Multi-stage time-stepping numerical simulation method of wave equation-based visco-acoustic model

Abstract

The application discloses a kind of multistage time steping numerical simulation method of wave equation-based lossy model, which comprises the following steps: modeling a model with ultra-thin lossy structure; the model is meshed, and the solution domain is spatially dispersed; derive the second-order wave equation containing lossy medium, introduce the model boundary condition, and carry out Galerkin test; the second-order wave equation is combined with spatial dispersion to convert into a lossy global matrix partial differential equation; all unknowns of the partial differential equation are written as Taylor expansion, and the partial differential equation is substituted into the Taylor expansion; the order of the high-order partial derivative term is reduced; the unknowns in the Taylor expansion after order reduction are solved using a multistage stepping scheme; the stability condition of multistage time stepping is given; the acceleration ratio is obtained through theoretical analysis, and the test verification is carried out. The application can not only use flexible spatial dispersion method, but also realize high-precision and stable numerical solution, especially for ultra-thin lossy structure with high computational efficiency, significantly reducing the calculation cost.

Description

Technical Field

[0001] This invention relates to a transient and efficient simulation method for lossy models, and in particular to a numerical simulation method for a multi-stage time-stepping lossy model based on the wave equation. Background Technology

[0002] The wave equation can be derived from Maxwell's equations. However, the wave equation only needs to store the unknowns of the electric field, while Maxwell's equations need to store the unknowns of both the electric and magnetic fields, which will result in greater memory overhead. Meanwhile, both [C.-Y.Tian, ​​Y.Shi, KMShum, and CHChan, "Wave equation-based discontinuous Galerkin time-domain method for co-simulation of electromagnetic-circuit systems," IEEE Trans. Antennas Propag., vol. 68, no. 4, pp. 3026–3036, Apr. 2020.] and [C.-Y.Tian, ​​Y.Shi, and CHChan, "Interior penalty discontinuous Galerkin time-domain method based on wave equation for 3-D electromagnetic modeling," IEEE Trans. Antennas Propag., vol. 65, no. 12, pp. 7174–7184, Dec. 2017.] show that, under the same numerical analysis algorithm, the wave equation has a higher simulation efficiency than the Maxwell equation.

[0003] High-conductivity lossy model materials are widely used due to their low internal resistance, high current carrying capacity, significant skin effect, and rapid attenuation of electromagnetic waves. For example, coating the outer surface of certain equipment with conductive composite materials can absorb radar waves and reduce target detection signals. However, these coatings are usually very thin, forming a cross-scale feature with the surrounding structure. This poses a challenge to transient full-wave simulation, as the global time step is limited by these lossy structures, necessitating the exploration of more efficient simulation strategies. Summary of the Invention

[0004] The purpose of this invention is to provide a numerical simulation method for a multi-level time-stepping lossy model based on the wave equation. This method can not only describe the model with infinite accuracy by increasing the order of the Taylor expansion, but also significantly improve computational efficiency and reduce computational cost for lossy models, especially high conductivity structures with ultrathin structures.

[0005] The technical solution to achieve the purpose of this invention is as follows:

[0006] A numerical simulation method for a multi-stage time-stepping lossy model based on the wave equation includes the following steps:

[0007] Step 1: Model the lossy model with ultrathin structure, perform mesh generation, spatially discretize the solution domain, and derive mesh and node information;

[0008] Step 2: Derive the second-order wave equation containing lossy media, introduce the model boundary conditions, and perform Galerkin tests to transform the frequency domain second-order wave equation into a lossy global matrix partial differential equation by combining spatial discretization.

[0009] Step 3: Write all the unknowns of the depleted global matrix partial differential equation as Taylor expansions, and substitute the depleted global matrix partial differential equation into the Taylor expansions to reduce the order of higher-order partial derivatives.

[0010] Step 4: Solve the unknowns in the Taylor expansion after order reduction using a multi-step approach. Solve the unknowns in the ultrathin lossy structure using a small time step, while solve the unknowns in other regions using a large time step, to obtain the numerical simulation values ​​of the lossy model.

[0011] Step 5: List the stability conditions for multi-level time stepping and perform critical multi-level time stepping tests;

[0012] Step 6: Based on the numerical discretization results of the lossy structure obtained in Step 4, the speedup ratio is obtained through theoretical analysis, and the simulation resource overhead test is performed based on the time consumption in Step 5. The speedup ratio of the overhead test is compared with the speedup ratio of the theoretical analysis. If the difference between the speedup ratio of the overhead test and the speedup ratio of the theoretical analysis does not meet the requirements, the numerical simulation needs to be repeated.

[0013] Compared with the prior art, the significant advantages of this invention are:

[0014] (1) The method of the present invention adopts multi-level time stepping technology, which can solve the problem of time multi-scale caused by the time stability condition of lossy and lossless models, and significantly improve the computational efficiency.

[0015] (2) For the lower-order terms in the higher-order expansion, the lower-order expansion can be used directly by taking advantage of the stacking relationship; while the higher-order partial derivatives can be reduced to the lower-order terms, avoiding the complex differentiation of the higher-order terms and the consumption of unnecessary storage resources.

[0016] (3) Compared with the global unified time stepping technique, no new unknowns are generated and memory consumption is almost not increased.

[0017] (4) This invention combines the relationship between the local coordinate system and the reference coordinate system of PML at any rotation angle to derive a method for handling the boundary of the lossy model of PML+ABC at any angle, which can be more flexibly applied to the simulation of irregular lossy models. Attached Figure Description

[0018] Figure 1 This is a schematic diagram of a high-power radio frequency absorbing device with an ultra-thin lossy structure.

[0019] Figure 2 This is a schematic diagram showing the electric field strength at the field probe of a high-power RF absorbing device (comparison of global unified time stepping and multi-level time stepping methods).

[0020] Figure 3 This is a schematic diagram showing the reflection coefficient of a high-power radio frequency absorbing device (comparison of global unified time stepping and multi-level time stepping methods). Detailed Implementation

[0021] This embodiment provides a numerical simulation method for a multi-level time-stepping lossy model based on the wave equation, including the following steps:

[0022] First, a lossy structure is modeled. Generally, due to the high conductivity and small skin depth of lossy structures, only a very thin lossy structure is needed to achieve the function. The ultra-thin structure forms a multi-scale structure with the dielectric substrate and air. To facilitate flexible mesh generation, a more flexible tetrahedron can be used to model the local region, and then uniformly converted into a hexahedron for spatial discretization of the solution domain. The mesh number, the corresponding material number in the mesh, the corresponding node number in the mesh, and the node spatial coordinates are exported to complete the model construction.

[0023] Secondly, the wave equation of the lossy model can be written as:

[0024]

[0025] Where μ is the magnetic permeability, ε is the permittivity, σ is the electrical conductivity, t is time, and E is the unknown electric field. Based on the expression for heat power density P in the electrothermal relationship... d =σ|E| 2It can be seen that this model is lossy, with electric field energy being converted into heat power and dissipated. To ensure good absorption, a PML (perfectly matched layer) is used for absorption at the boundary, and an ABC (absorption boundary condition) is used to truncate the computational threshold. The frequency domain wave equation containing PML can be written as:

[0026]

[0027] in The matching matrix for the corner region. ξ takes values ​​x, y, z, σ ξ Let ω represent the conductivity distribution in different directions, ω be the angular frequency, and ε₀ be the vacuum permittivity. If a reference coordinate system is used to set up the PML, then according to the mapping relationship between the reference coordinate system (u, v, w) and the global coordinate system (x, y, z), we know that:

[0028]

[0029]

[0030] Therefore, if the PML parameters are set using a reference coordinate system, then the vertex matrix [Λ] xyz The relationship to be expressed as the reference vertex matrix [Λ] is as follows:

[0031] [Λ] xyz =Q T [Λ]Q (6)

[0032] The frequency domain wave equation in a reference coordinate system with PML can then be written as:

[0033]

[0034] Expanding the unknown electric field in the above equations and performing Galerkin tests, we obtain:

[0035]

[0036] Where the subscripts ei and ej represent the numbers of the unknowns, N ej To expand the basis functions, N ei To test the basis functions. Additionally, the solution domain outside the PML is truncated using ABC, which can be represented as:

[0037]

[0038] Substituting the ABC boundary into the above expansion equation, we get:

[0039]

[0040] Will Substitute into the above equation, and multiply both sides by S.x S y S z ,have to:

[0041]

[0042] in, Therefore, we can conclude that:

[0043]

[0044] And because σ x σ y σ z The three terms and σ cannot all be non-zero at the same time, as they belong to different regions. Therefore, the last term can be simplified directly:

[0045]

[0046]

[0047] in

[0048]

[0049] Multiplying both sides by the permeability, we transition to the time domain:

[0050]

[0051] We obtain a compact form of the second-order wave matrix equation:

[0052]

[0053] in

[0054]

[0055] [C] ij =[T] ij =με∫N ei ·N ej dv (25)

[0056]

[0057] Based on the basis functions (Φ) of the mapping space (ξ, η, ζ) i Φ i With divergence and curl compatibility, we obtain the degeneracy matrix:

[0058]

[0059] [C] ij =[T] ij =με∫(J -1 Φ i)·(J -1 Φ j )|J|dξdηdζ (30)

[0060]

[0061] Then, the compact form of the second-order wave matrix equation obtained in (21) is transformed into a system of equations, and the first-order partial derivatives of each unknown in the system of equations can be written as:

[0062]

[0063]

[0064] To facilitate the demonstration of the derivation process, all unknowns are reduced to the transient scheme. After the first-order Taylor expansion, the first-order partial derivatives of each unknown are substituted in:

[0065]

[0066] After performing a second-order Taylor expansion, substituting it into the first-order Taylor expansion, we obtain:

[0067]

[0068] Similarly, depending on the required precision, Taylor expansions of any order can be performed on the unknowns. For lower-order terms within higher-order expansions, the lower-order expansions can be directly used by leveraging the stacking relationship; and higher-order partial derivatives can be reduced to lower-order terms, avoiding complex differentiation of higher-order terms and unnecessary storage resource consumption.

[0069] Next, the unknowns in the reduced-order Taylor expansion are solved using a multi-stage stepping scheme, while the unknowns in the ultrathin lossy structure are solved using a small time step Δt. s The solution is obtained by using a large time step Δt for the unknowns in the lossless region. lThe solution involves finding the time step ratio N between the lossless and lossy regions. When solving for the unknowns in both the lossy and lossless regions, the vectors are divided according to different time steps. When solving for the lossless region, the partial derivatives of the lossy region are still used; similarly, when solving for the lossy region, the partial derivatives of the lossless region are also applied. Since the unknowns are solved using different time steps, the vectors solved at smaller time steps are interpolated to ensure consistency in the global time solution. The interpolation process is as follows: the unknowns after the first hour step in the lossy region are obtained by using the partial derivatives of the unknowns from the previous large time step and all regions; the unknowns after the m-th hour step (2≤m≤N, and are integers) are obtained by using the unknowns after the (m-1)-th hour step and all regions. Due to the layering relationships in Taylor expansions of different orders and the coupling relationship between partial derivatives and unknowns, when solving higher-order expansions, it is only necessary to store the results of each polynomial in the previous order expansion.

[0070] The stability condition of multi-level time stepping is related not only to the size of the mesh in different regions, but also to the material parameters of each region (such as dielectric constant, magnetic permeability and electrical conductivity).

[0071] For the long-term stability of the lossless region, the time step is mainly affected by the fluctuation term and needs to satisfy the following empirical formula:

[0072]

[0073] For the long-term stability of the lossy region, the time step is mainly affected by the lossy term and needs to satisfy the following empirical formula:

[0074]

[0075] Where C is an empirical constant, related to the partitioning type and spatial dimension, l min Let c0 be the minimum grid size of the region, c0 be the speed of light in vacuum, and μ be the minimum grid size of the region. r To solve for the permeability of the region, ε r To solve for the dielectric constant of the region, it can be seen that in the lossy region, where the conductivity is very high (≥10^4 S / m), the choice of time step will be completely dominated by the conductivity.

[0076] The initial time step ratio N0 can also be selected using the above formula:

[0077]

[0078] Furthermore, due to the influence of the partitioning type, numerical discretization method, spatial dimension, etc., Δt l , Δt sThe actual time step ratio N should be selected as a critical value through multiple tests.

[0079] Finally, since the relationship between the unknown quantity and the machine solution time is approximately linear, the theoretical speedup of the proposed multi-stage time-stepping method can be written as:

[0080]

[0081] Where T 全局 T 局部 DOFs 全局 DOFs 有耗 DOFs 无耗 These figures represent the computation time consumed by a globally unified time step, the computation time consumed by a multi-level time step, the number of all unknowns globally, the number of unknowns in the costly region, and the number of unknowns in the costless region. Since the computational matrix has a relatively uniform dimension, it does not incur excessive computational memory consumption. The proposed method and its acceleration effect can be verified through costly model simulations.

[0082] Example

[0083] The implementation effect is demonstrated using a high-power radio frequency absorbing device that integrates a novel 0.1µm thick ceramic material onto a high thermal conductivity SiC substrate 11 as an example. Figure 1 As shown, the RF absorbing device also includes a conductive microstrip line 13 and a metal via 10. The ultrathin lossy structure 12 formed by the novel ceramic material has a melting point greater than 3400℃, which allows it to maintain good thermal reliability when absorbing high-power microwaves. Simultaneously, the novel ceramic material possesses an ultra-high conductivity of 350,000 S / m and a relative permittivity of 4.2. The relative permittivity of the SiC substrate is 9.7. The dimensional parameters of the high-power RF absorbing device are shown in Table 1. In the numerical simulation, a 0.3 mm thick PML absorbing layer was set in the x, y, and z directions of the model, and the outermost layer adopted an ABC cutoff computational domain. A modulated Gaussian source with a center frequency of 15 GHz is injected from port A. The total simulation time is 200 ps. Compared with the globally unified time-stepping method, which uses a maximum time step of 0.167 fs, the multi-level time-stepping technique used in this invention can achieve a time step of 100 times that of the globally unified time step in the lossless region, using a time step of 16.7 fs. However, a time step of 0.167 fs is still used in the lossy region. Due to the high conductivity of the new ceramic material, the skin depth of electromagnetic waves is relatively small, and the electromagnetic waves attenuate quickly, requiring an extremely small time step to ensure the convergence of this method. Figure 2 The transient time-varying electric field was obtained using two methods for the field probe. Figure 3The reflection coefficients of PortA, simulated by two methods, show significant return loss below 22 GHz. Figure 2 and Figure 3 Transient simulations show a high degree of consistency between the two methods, verifying that the method proposed in this patent does not reduce simulation accuracy compared to the global time-stepping method. Furthermore, as shown in Table 2, the multi-level time-stepping scheme reduces computation time by more than 90% while maintaining the same level of memory consumption. This further highlights the superior computational efficiency of the multi-level time-stepping lossy model numerical simulation method based on the wave equation.

[0084] Table 1 Structural Dimensions of High-Power RF Absorbing Devices

[0085]

[0086] Table 2 Computational resource consumption during the simulation of high-power RF absorbing devices

[0087]

[0088] This invention proposes a multi-level time-stepping numerical simulation method for lossy models based on the wave equation. In ultrathin structures with extremely high conductivity, a sufficiently small time step ensures the stability of the lossy structure, while in the lossless region, the time step can be amplified according to the stability relationship, achieving higher computational efficiency. Simultaneously, by utilizing Taylor expansions of different orders with stacked relationships, high-order accuracy numerical simulations can be conveniently achieved. Compared with globally unified time-stepping techniques, no new unknowns are generated, and memory consumption is almost negligible. This invention not only employs flexible spatial discretization methods but also achieves high-precision and stable numerical solutions, exhibiting particularly high computational efficiency for ultrathin lossy structures and significantly reducing computational costs.

[0089] Obviously, those skilled in the art can make various modifications and variations to the embodiments of the present invention without departing from the spirit and scope of the embodiments of the present invention. Thus, if these modifications and variations to the embodiments of the present invention fall within the scope of the claims of the present invention and their equivalents, the present invention also intends to include these modifications and variations.

Claims

1. A numerical simulation method for a multi-stage time-stepping lossy model based on the wave equation, characterized in that, Including the following steps: Step 1: Model the lossy model with ultrathin structure, perform mesh generation, spatially discretize the solution domain, and derive mesh and node information; Step 2: Derive the second-order wave equation containing lossy media, introduce the model boundary conditions, and perform Galerkin tests to transform the frequency domain second-order wave equation into a lossy global matrix partial differential equation by combining spatial discretization. Step 3: Write all the unknowns of the depleted global matrix partial differential equation as Taylor expansions, and substitute the depleted global matrix partial differential equation into the Taylor expansions to reduce the order of higher-order partial derivatives. Step 4: Solve the unknowns in the Taylor expansion after order reduction using a multi-step approach. Solve the unknowns in the ultrathin lossy structure using a small time step, while solve the unknowns in other regions using a large time step, to obtain the numerical simulation values ​​of the lossy model. Step 5: List the stability conditions for multi-level time stepping and perform critical multi-level time stepping tests; Step 6: Based on the number of unknowns in the lossy model obtained in Step 4, calculate the theoretical speedup ratio and obtain the test speedup ratio based on the test in Step 5. If the difference between the test speedup ratio and the theoretical speedup ratio does not meet the requirements, the numerical simulation needs to be repeated.

2. The numerical simulation method for a multi-level time-stepping lossy model based on the wave equation according to claim 1, characterized in that, The lossy model is modeled using tetrahedrons, and then uniformly converted into hexahedrons for spatial discretization of the solution domain. The grid number, the corresponding material number in the grid, the corresponding node number in the grid, and the node spatial coordinates are exported to complete the model construction.

3. The numerical simulation method for a multi-level time-stepping lossy model based on the wave equation according to claim 1, characterized in that, PML is used for absorption at the model boundary. The second-order wave equation in the frequency domain containing PML is: Where μ is the magnetic permeability, ε is the permittivity, σ is the conductivity, and E is the unknown electric field. For the matching matrix, ξ takes values ​​x, y, z, σ ξ Let ω represent the conductivity distribution in different directions, and ω be the angular frequency. Here, x, y, z represent the density operator, and x, y, z represent the coordinates in the global coordinate system.

4. The numerical simulation method for a multi-level time-stepping lossy model based on the wave equation according to claim 3, characterized in that, The parameters of PML are set using a reference coordinate system, and the second-order wave equation in the frequency domain is: Wherein, the mapping matrix between the reference coordinate system and the global coordinate system QQ T =I, where u, v, and w are coordinates in the reference coordinate system.

5. The numerical simulation method for a multi-stage time-stepping lossy model based on the wave equation according to claim 4, characterized in that, The second-order wave equation in the frequency domain is transformed into a lossy global matrix partial differential equation by combining spatial discretization: in, e is an unknown quantity. Considering the divergence and curl compatibility between the mapping space (ξ,η,ζ) and the original space (u,v,w), the coefficient matrices in the equation are as follows: [C] ij =with∫(J -1 F i )·(J -1 F j )|J|dξdηdζ Where the subscripts i and j represent the numbers of the basis functions Φ, J is the Jacobian matrix, and n is the normal vector.

6. The numerical simulation method for a multi-level time-stepping lossy model based on the wave equation according to claim 5, wherein in step 3, all unknowns of the lossy global matrix partial differential equation are written as Taylor expansions, and the lossy global matrix partial differential equation is substituted into the Taylor expansion to perform order reduction processing of higher-order partial derivative terms, the specific process being as follows: The costly global matrix partial differential equation is transformed into a system of equations, and the first-order partial derivatives of each unknown in the system of equations are written as: Restore all unknowns to the transient scheme, perform a first-order Taylor expansion, and substitute the first-order partial derivatives of each unknown: After applying the second-order Taylor expansion and substituting it into the first-order Taylor expansion, we get: Similarly, depending on the required precision, Taylor expansions of any order are performed on the unknowns. For lower-order terms within higher-order expansions, the lower-order expansions are directly used based on the stacking relationship; higher-order partial derivatives are reduced to lower orders.

7. The numerical simulation method for a multi-level time-stepping lossy model based on the wave equation according to claim 1, characterized in that, In step 4, when the unknowns in the ultrathin lossy structure are solved using small time steps, while the unknowns in other regions are solved using large time steps, the specific stepping process is as follows: the unknowns in the ultrathin lossy structure are solved using small time steps Δt. s The solution is obtained by using a large time step Δt for the unknowns in the lossless region. l The solution involves finding the ratio of time steps between the lossless and lossy regions, where N is the number of steps. Solving for the lossless region requires using partial derivatives of the lossy region; similarly, solving for the lossy region requires using partial derivatives of the lossless region. The interpolation stepping process is as follows: the unknown quantity after stepping through the first hour in the lossy region is obtained by using the partial derivatives of the unknown quantities from the previous large time step and all regions; the unknown quantity after stepping through the m-th hour is obtained by using the unknown quantity after stepping through the (m-1)-th hour and all regions, where 2 ≤ m ≤ N and is an integer.

8. The numerical simulation method for a multi-stage time-stepping lossy model based on the wave equation according to claim 4, characterized in that, The stability conditions for multi-level time stepping in step 5 include: For the lossless region, the time step satisfies: For the lossy region, the time step satisfies: Where C is an empirical constant, related to the partitioning type and spatial dimension, l min Let c0 be the minimum grid size of the region, c0 be the speed of light in vacuum, and μ be the minimum grid size of the region. r To solve for the permeability of the region, ε r To determine the dielectric constant of the solution region; the initial time step ratio N0 is: 。 9. The numerical simulation method for a multi-stage time-stepping lossy model based on the wave equation according to claim 4, characterized in that, The theoretical speedup ratio in step 6 is: Where T 全局 T 局部 DOFs 全局 DOFs 有耗 DOFs 无耗 These represent the computation time consumed by a globally unified time step, the computation time consumed by a multi-level time step, the number of all unknowns to be solved globally, the number of unknowns to be solved in a costly region, and the number of unknowns to be solved in a costless region.

Images