An Efficient Time Domain Domain Decomposition Algorithm Based on the Spectral Element Method in Time Domain

Abstract

The present invention discloses an efficient time-domain domain decomposition algorithm based on the time-domain spectral element method. This algorithm performs sub-region modeling on the electromagnetic space, constructs the corresponding numerical electric field model in each region respectively, and imposes the electric field transmission conditions on the interfaces between adjacent regions. At the same time, according to the time-domain spectral element method, the numerical electric field models in each region are discretized in space and time, and finally an efficient time-domain domain decomposition algorithm with sub-region and explicit solution characteristics is formed. In the spatial discretization of the present invention, the technology of classifying and separating the electric field inside the region and the electric field on the region interface is adopted to realize the sub-region solution process with completely independent unknowns; in the time discretization, the central difference scheme is adopted to ensure that the algorithm in each region has an efficient explicit solution format. The present invention not only has a flexible modeling method, but also can achieve efficient time-domain electromagnetic simulation of finite periodic structures through multiple numerical extensions of typical modeling regions, with a fast simulation speed and small computer memory occupancy.

Description

Technical Field

[0001] The present invention belongs to the technical field of electromagnetic simulation, and particularly relates to an efficient time-domain domain decomposition algorithm based on the spectral element method in the time domain. Background Art

[0002] Most modern artificial electromagnetic materials have periodic structures, and with the increase in the number of unit structure periods, actual artificial electromagnetic materials usually have electrically large sizes, which greatly increases the difficulty of achieving efficient numerical electromagnetic simulation. To solve this problem, periodic boundary conditions are usually set around the unit structure to minimize the unknowns, thereby reducing memory consumption and improving calculation efficiency. However, the periodic boundary conditions can only simulate the electromagnetic characteristics of materials in the case of infinite periodic structures, while actual materials must be finite periodic structures, and the electromagnetic coupling effect at the material boundary cannot be ignored.

[0003] For the efficient electromagnetic simulation of electrically large finite periodic structures, the domain decomposition method (DDM) is an effective solution. Using the domain decomposition method, sub-region modeling can be performed on several unit structures with typical boundary conditions, and according to the periodic characteristics of the material, the typical regions are extended a finite number of times to construct the electromagnetic simulation of the finite periodic structure.

[0004] The domain decomposition method can be implemented through frequency-domain algorithms or time-domain algorithms. When the domain decomposition method is combined with frequency-domain algorithms, usually only the electric field Helmholtz equation needs to be solved, without involving the solution of the magnetic field, so the unknowns are few and the calculation speed is fast. However, due to the frequency-sweeping characteristics of frequency-domain full-wave simulation, it usually results in low efficiency when obtaining broadband frequency responses using this algorithm, while time-domain algorithms can obtain frequency responses in a wide frequency band range with a relatively small computational cost. However, when the existing domain decomposition method is combined with time-domain algorithms, the electric field Maxwell equation and the magnetic field Maxwell equation need to be solved separately, resulting in a doubling of the unknowns. Therefore, combining the domain decomposition method with a time-domain algorithm that only calculates the electric field Helmholtz equation is an effective way to achieve efficient electromagnetic simulation of broadband frequency responses for electrically large finite periodic structures. Summary of the Invention

[0005] The object of the present invention is to propose a new time-domain domain decomposition algorithm based on the spectral element method in the time domain, which can achieve efficient time-domain electromagnetic simulation of finite periodic structures. This algorithm performs sub-region processing on the electromagnetic structure according to the periodic characteristics, obtains the numerical models of each region based on the time-domain electric field Helmholtz equation and the time-domain electric field transmission conditions on the domain interfaces, then uses the electric field classification and separation technology and the central difference scheme to discretize the electric fields in each region in space and time, and finally realizes the efficient time-domain electromagnetic simulation of finite periodic structures through multiple periodic numerical extensions of the typical regions.

[0006] The technical solution for achieving the object of the present invention is as follows: An efficient time-domain domain decomposition algorithm based on the time-domain spectral element method, comprising the following steps:

[0007] Step 1, perform sub-region modeling on the electromagnetic space. On the basis of fully considering regional associations, construct the numerical model of the electric field for each region to form the weak solution expression of the electric field in each region.

[0008] Step 2, in each region, use the central difference format to discretize the weak solution expression of the electric field in time to obtain its corresponding time-domain iterative expression of the electric field.

[0009] Step 3, in each region, adopt the classification and separation technique to achieve an orderly arrangement of the discrete spatial electric field.

[0010] Step 4, utilize the orderly arrangement characteristics of the electric field to establish an independent solution mode for the electric field in sub-regions.

[0011] Step 5, load the excitation source and set the total number of calculation time steps. During the time stepping process, first calculate the electric field on the regional interface, and then calculate the internal electric field of each sub-region until the set number of time steps is completed.

[0012] Compared with the prior art, the significant advantages of the present invention are as follows:

[0013] 1. Compared with the frequency-domain domain decomposition method:

[0014] The frequency-domain domain decomposition method can only obtain the electromagnetic response at one frequency point each time. To obtain the broadband frequency response, it is necessary to calculate point by point for each frequency point, so the efficiency is low. While the present invention is based on the time-domain algorithm. After calculating the time-domain waveform of the signal, through a single fast Fourier transform, the frequency response in a very wide frequency band range can be obtained. Therefore, for the calculation of broadband frequency response, the time-domain domain decomposition method is more efficient.

[0015] 2. Compared with the existing time-domain domain decomposition method:

[0016] (1) The novel time-domain domain decomposition method of the present invention is based on the electric field Helmholtz equation, while the existing time-domain domain decomposition method is based on the electric field Maxwell equation and the magnetic field Maxwell equation, and the number of unknowns is twice that of the present invention.

[0017] (2) The novel time-domain domain decomposition method of the present invention only needs to calculate the electric field Helmholtz equation once at each time step, while the existing time-domain domain decomposition method needs to calculate the electric field Maxwell equation and the magnetic field Maxwell equation once at each time step, and the calculation efficiency is low.

[0018] 3. The present invention combines all the advantages of the domain decomposition method and the time-domain spectral element method:

[0019] (1) Sub-region modeling is adopted, and a coarse-grid curvilinear hexahedron is used for modeling. The modeling method is flexible, the fitting degree of the model structure is high, and the number of calculated variables formed is small.

[0020] (2) On the basis of using the high-order orthogonal basis functions of the time-domain spectral element method, an electric field classification and separation technique is adopted, making the electric field inside each region independent of the electric field at the region interface. Parallel computing is carried out, perfectly inheriting the explicit solution characteristics of the time-domain spectral element method, with fast calculation speed and high calculation efficiency.

[0021] (3) It is suitable for the use of the periodic numerical continuation technique. When modeling, only a small number of physical models of typical periodic structures need to be constructed. Through numerical continuation and matrix reuse, an efficient electromagnetic simulation of electrically large periodic targets is formed, greatly reducing the memory occupation of the computer.

[0022] The present invention will be further described in detail below with reference to the accompanying drawings. Description of the Drawings

[0023] Figure 1 It is a flow chart of a new time-domain domain decomposition algorithm based on the time-domain spectral element method.

[0024] Figure 2 It is a schematic diagram of a waveguide structure.

[0025] Figure 3 It is a schematic diagram of the S-parameters of a waveguide with 4 repeating structures.

[0026] Figure 4 It is a schematic diagram of the S-parameters of a waveguide with 8 repeating structures.

[0027] Figure 5 It is a schematic diagram of the S-parameters of a waveguide with 12 repeating structures.

[0028] Figure 6 It is a comparison chart of the calculation time between the present invention and the traditional time-domain spectral element method.

[0029] Figure 7 It is a comparison chart of the memory occupation between the present invention and the traditional time-domain spectral element method. Specific Embodiments

[0030] Combined with Figure 1 , an efficient time-domain electromagnetic simulation of a finite periodic structure is carried out using a new time-domain domain decomposition algorithm based on the time-domain spectral element method, including the following steps:

[0031] Step 1: Model by region according to the characteristics of the periodic structure, discretize the physical model using curvilinear hexahedron elements, construct a numerical model for each region using the time-domain electric field Helmholtz equation, and impose the time-domain electric field transmission condition on the adjacent region interfaces.

[0032] In each region, the time-domain electric-field Helmholtz equation can be expressed as:

[0033]

[0034] where is the Hamiltonian operator, E is the electric-field strength, ε is the permittivity, and μ is the permeability. The imposed time-domain electric-field transmission conditions at the interface between adjacent regions are:

[0035]

[0036] where the subscript a is the identifier of this region, b is the identifier of the adjacent region, E a is the tangential electric-field strength at the interface in this region, E b is the tangential electric-field strength at the interface in the adjacent region, Z a , Z b are the wave impedances of the two regions, Z 0 is the wave impedance in vacuum, is the unit normal vector of the interface, denotes the partial derivative with respect to time t.

[0037] Using high-order orthogonal basis functions as the expansion basis functions for the electric field in Equation (1), and the same basis functions as the test basis functions to perform the Galerkin transformation on Equation (1). At the same time, adding the electric-field transmission conditions at the interface between adjacent regions in Equation (2) to this transformation process, the weak-solution expression form of the electric field in each region is obtained:

[0038]

[0039] where

[0040]

[0041] T = ε∫∫∫N am ·N an dV

[0042]

[0043]

[0044]

[0045]

[0046] In the above equation, e a represents all the unknown electric-field quantities to be solved in this region, e ai and e birespectively represent the unknown electric field quantities to be solved on the regional interface of this region and the adjacent region. Matrix T and matrix S are volume integral matrices, and matrix T 1 、T 2 、S 1 、S 2 are surface integral matrices. In the integral terms for forming the elements of each matrix, N am 、N an respectively represent the m-th test basis function and the n-th expansion basis function of this region, and N bn represents the n-th expansion basis function of the adjacent region on the interface.

[0047] Step 2: Use the central difference scheme to discretize Equation (3) in time. Obtain the time-domain iterative format of the electric field with Δt as the time step:

[0048]

[0049] where Δt represents the time step, and the superscripts n - 1, n, and n + 1 are identifiers for three adjacent time steps. The n + 1 moment is the current moment for solving the unknown quantity. Obviously, the unknown quantity at the current moment solved in Equation (4) includes those of both this region and the adjacent region. Therefore, Equation (4) does not satisfy the computational characteristics of domain decomposition.

[0050] Step 3: Classify and separate the left-hand side term of Equation (4), that is, first classify and divide the electric field in each region into the internal electric field of the region and the electric field on the regional interface, and then reorder the unknown electric field quantities in each region in the order of the internal electric field first and then the interface electric field. After classification and separation, the left-hand side term of Equation (4) can be expressed as:

[0051]

[0052] where e ar represents the unknown electric field quantity to be solved inside this region, and e ai and e bi respectively represent the unknown electric field quantities to be solved on the regional interface of this region and the adjacent region.

[0053] Step 4: Aggregate the unknown electric field quantities on the interface to form a solution region specifically for the unknown electric field quantities inside each region and a solution region specifically for the unknown electric field quantities on the interface respectively. For example, assume a model is decomposed into region a and region b. Then, classify and separate the unknown electric field quantities in their respective regions, and also perform corresponding classification and separation on the elements of each related matrix. The finally formed global equations can be simplified and expressed as:

[0054]

[0055] where,

[0056] T′ ai = T ai + 0.5ΔtT 1ai

[0057] T″ ai = 0.5ΔtT 2ai

[0058] T′ bi = T bi + 0.5ΔtT 1bi

[0059] T″ bi = 0.5ΔtT 2bi

[0060] b ar = (2T ar - Δt 2 S ar )e a n - T ar e ar n-1

[0061] b br = (2T br - Δt 2 S br )e b n - T br e br n-1

[0062] b ai = (2T ai - Δt 2 S ai - Δt 2 S 1ai )e a n +(0.5ΔtT 1ai - T ai )e ai n-1 - Δt 2 S 2ai e bi n + 0.5ΔtT 2ai e bi n-1

[0063] b bi = (2T bi - Δt 2 S bi - Δt 2 S1bi )e b n +(0.5ΔtT 1bi -T bi )e bi n-1 -Δt 2 S 2bi e ai n +0.5ΔtT 2bi e ai n-1

[0064] In the above formula, the subscripts ar and br respectively represent the identifiers of the electric field related quantities inside the a-region and the b-region, the subscripts ai and bi respectively represent the identifiers of the electric field related quantities on the interface between the a-region and the b-region, and the subscripts a and b respectively represent the identifiers of all the electric field related quantities in the a-region and all the electric field related quantities in the b-region.

[0065] The global coefficient matrix of Equation (6) has an obvious block diagonal property and can be split into three region-characteristic equations, which are respectively expressed as:

[0066] [T ar [e ar n+1 =[b ar (7-1)

[0067] [T br [e br n+1 =[b br (7-2)

[0068]

[0069] Among them, Equation (7-1) is used to solve the unknown electric field in region a, Equation (7-2) is used to solve the unknown electric field in region b, and Equation (7-3) is used to solve the unknown electric field on the interface between region a and region b. The solutions of these three equations are completely independent, and only the electric field information of the first two time steps in the adjacent regions at the interface element is used when forming the right vector b of the equations.

[0070] And so on. As the number of regions increases, independent equations for solving the internal electric field and the interface electric field will also be formed in the newly added regions; if the electromagnetic model is a periodic structure, numerical continuation can be performed on the repeated regions of the periodic structure, and the electric field transmission conditions are imposed on the newly generated region interfaces between the extended regions and the regions to be extended, so as to obtain the numerical model of the extended regions. The matrices of the extended regions can be reused, thus greatly reducing the memory occupation of the computer.

[0071] ​​Step 5: Load the excitation source and set the total number of computational time steps. In each time step, first calculate the electric field on the regional interface, and then calculate the electric field inside each region until the set number of time steps is completed, and post-process the obtained data.

[0072] Taking a rectangular waveguide structure filled with a full-height dielectric as an example, the present invention is compared with the traditional spectral element time domain method (SETD).

[0073] Embodiment

[0074] As Figure 2 shown in the waveguide structure. Among them, the size of the waveguide cross-section is 20mm×5mm, the filled dielectric block is located in the middle of the waveguide, its size is 20mm×5mm×10mm, and the relative dielectric constant of the dielectric block is ε r = 2. 16 layers of PML absorbing boundaries are set at both ends of the waveguide, and the thickness of each layer of PML is 5mm. The excitation source is a modulated Gaussian pulse with a center frequency of 9.5 GHz and a signal bandwidth of 2 GHz, and the source application surface is located in Region 1. Two observation points are set, located on both sides of the dielectric block in Region 2 respectively.

[0075] Use the new time-domain domain decomposition algorithm based on the spectral element time domain method to perform efficient time-domain electromagnetic analysis on this waveguide structure. The specific steps are as follows:

[0076] Step 1, according to the structural characteristics of the waveguide, divide it into 3 typical regions, as Figure 2 shown. Among them, Region 1 and Region 3 contain PML absorption layers, and Region 2 is a discontinuous structure layer that can be periodically extended. Use computer-aided design tools to model the typical regions, and discretize the model using hexahedral elements to obtain the geometric parameter information of all elements.

[0077] Step 2, according to the new time-domain domain decomposition method based on the spectral element time domain method proposed by the present invention, respectively construct the solution equations for the unknown electric field quantities inside the 3 typical regions and the solution equations for the unknown electric field quantities on the 2 regional interfaces. Set the calculation to 10,000 steps, and record the time-domain electric field values of the 2 observation points at each time step until the set number of time steps is completed.

[0078] Step 3, perform a fast Fourier transform on the time-domain waveforms recorded at the 2 observation points to obtain the corresponding frequency-domain information, and then calculate the S11 parameter and S21 parameter according to the definition of the scattering parameters.

[0079] Step 4, use the periodic numerical extension technique to extend Region 2. Using the new time-domain domain decomposition method based on the spectral element time domain method proposed by the present invention, calculate the S11 parameter and S21 parameter of the waveguide when Region 2 is extended 3 times, 7 times, 11 times, 15 times, and 19 times respectively.

[0080] Step 5: According to the actual extension times of Region 2, establish corresponding full-scale models, and calculate their corresponding S11 parameters and S21 parameters respectively using the time-domain spectral element method and CST commercial simulation software, and compare with the calculation results in Step 4. Figures 3 to 5 The waveguide S parameters obtained by the above three algorithms are compared when Region 2 is extended 3 times, 7 times, and 11 times, and the calculation results of the three algorithms are in complete agreement. Figure 6 and Figure 7 The differences in calculation time and memory occupancy between the traditional time-domain spectral element method (SETD) and the novel time-domain domain decomposition method (SETD-DDM) based on the time-domain spectral element method proposed in the present invention are compared. The comparison results confirm the advantages of the fast simulation speed and less computer memory occupancy of the present invention.

[0081] Thus far, the accuracy and efficiency of the time-domain electromagnetic simulation of the finite periodic structure of the present invention have been verified.

Claims

1. An efficient time-domain domain decomposition algorithm based on the spectral element method in the time domain, characterized in that, it includes the following steps: Step 1: Perform sub-region modeling on the electromagnetic space. On the basis of fully considering regional correlations, construct the numerical model of the electric field in each region to form the weak solution expression of the electric field in each region; Use the time-domain electric field Helmholtz equation to construct the basic numerical model of the electric field in each region: where ▽ is the Hamiltonian operator, E is the electric field strength, ε is the permittivity, and μ is the magnetic permeability of the electromagnetic wave propagation medium; Introduce the electric field transmission condition for the region interface. The tangential continuity of the electric field on the region interface is ensured by the time-domain electric field transmission condition: where the subscript a represents this region, and the subscript b represents the adjacent region, E a is the tangential electric field strength of this region on the interface, E b is the tangential electric field strength of the adjacent region on this interface, Z a and Z b are the wave impedances of the two adjacent regions respectively, Z 0 is the wave impedance in vacuum, is the unit normal vector of the interface, represents the partial derivative with respect to time t; Perform the Galerkin transformation on Equation (1) using high-order orthogonal basis functions; at the same time, add the electric field transmission condition on the interface between adjacent regions in Equation (2) to this transformation process to obtain the weak solution expression form of the electric field in each region: where In the above formula, e a represents all the electric fields to be solved in this region, e ai and e bi respectively represent the electric fields to be solved on the regional interface of this region and the adjacent region. The matrices T and S are volume integral matrices, and the matrix T 1 , T 2 , S 1 , S 2 are surface integral matrices. In the integral terms for forming each matrix element, N am , N an respectively represent the m-th test basis function and the n-th expansion basis function in this region, and N bn represents the n-th expansion basis function on the interface of the adjacent region; Step 2: In each region, use the central difference scheme to discretize the weak solution expression of the electric field in time to obtain its corresponding time-domain iteration expression of the electric field; Step 3: In each region, adopt the classification and separation technique to achieve an orderly arrangement of the discrete spatial electric field; Step 4: Utilize the orderly arrangement characteristics of the electric field to establish an independent solution mode for the electric field in sub-regions; Step 5: Load the excitation source and set the total number of calculation time steps. During the time stepping process, first calculate the electric field on the region interface, and then calculate the internal electric field of each sub-region until the set number of time steps is completed.

2. The efficient time-domain domain decomposition algorithm based on the spectral element method in the time domain according to claim 1, characterized in that, in Step 2, use the central difference scheme to discretize Equation (3) in time; obtain the time-domain iteration format of the electric field with Δt as the time step: where Δt represents the time step, and the superscripts n - 1, n, n + 1 are identifiers used to represent three adjacent time steps, and the n + 1 moment is the current moment for solving the unknown quantity.

3. The efficient time-domain domain decomposition algorithm based on the spectral element method in the time domain according to claim 2, characterized in that, in Step 3, use the classification and separation method for the spatially discrete electric field in each sub-region, that is, first classify the internal electric field of the region and the electric field on the region interface, then sort them in the order of the internal electric field first and the interface electric field later, and finally separate the electric fields on all region interfaces from each region to form an independent interface electric field region, and finally achieve independent solution of each region.

4. The efficient time-domain domain decomposition algorithm based on the spectral element method in the time domain according to claim 3, characterized in that, in Step 4, according to the classified and sorted electric field, recombine Equation (4) in each sub-region, and the left term of the recombined Equation (4) is expressed as: Among them, e ar represents the electric field to be solved inside this region, e ai and e bi respectively represent the electric fields to be solved on the region interface of this region and the adjacent region.

5. The efficient time-domain domain decomposition algorithm based on the spectral element method in the time domain according to claim 4, characterized in that, In step 4, the electric fields in each region are combined and analyzed, and the electric fields at the interfaces are summarized to form a global system of equations with the characteristic of sub-region calculation; that is, if a model is decomposed into region a and region b, in each region, according to the sorting method of the electric field in Equation (5), Equation (4) is recombined, which is manifested as the corresponding classification and separation of each matrix element related to Equation (4), and is reflected in the global matrix after summarizing each region; the finally formed global system of equations has the characteristic of sub-region calculation, and is simplified and expressed in the form of a block matrix vector multiplication as: Wherein, T a ′ i = T ai + 0.5ΔtT 1ai T a ″ i = 0.5ΔtT 2ai T b ′ i = T bi + 0.5ΔtT 1bi T b ″ i = 0.5ΔtT 2bi b ar = (2T ar - Δt 2 S ar )e a n - T ar e ar n-1 b br = (2T br - Δt 2 S br )e b n - T br e br n-1 b ai =(2T ai -Δt 2 S ai -Δt 2 S 1ai )e a n +(0.5ΔtT 1ai -T ai )e ai n-1 -Δt 2 S 2ai e bi n +0.5ΔtT 2ai e bi n-1 b bi = (2T bi - Δt 2 S bi - Δt 2 S 1bi )e b n + (0.5ΔtT 1bi - T bi )e bi n-1 -Δt 2 S 2bi e ai n +0.5ΔtT 2bi e ai n-1 In the above formula, the subscripts ar and br respectively represent the identifiers of the quantities related to the internal electric field in region a and region b, the subscripts ai and bi respectively represent the identifiers of the quantities related to the electric field at the interface between region a and region b, and the subscripts a and b respectively represent all the quantities related to the electric field in region a and all the quantities related to the electric field in region b.

6. The efficient time-domain domain decomposition algorithm based on the time-domain spectral element method according to claim 5, characterized in that, In step 4, the global system of equations has the characteristic of sub-region calculation, that is, Equation (6) is split into 3 systems of equations with the characteristic of independent sub-region solution, and are respectively expressed as: [T ar [e ar n+1 = [[b ar (7-1)[T br [[e br n+1 = [[b br (7-2)​​ Wherein, Equation (7-1) is used to solve the internal electric field in region a, Equation (7-2) is used to solve the internal electric field in region b, and Equation (7-3) is used to solve the electric field at the interface between region a and region b; the solutions of these 3 systems of equations are completely independent, except that the electric field information at the first 2 time steps in the interface elements of their respective adjacent regions is used when forming the right-hand vector b of the system of equations; If the electromagnetic model is a periodic structure, then the numerical continuation is performed on the repeated region of the periodic structure, and the electric field transmission condition is re-imposed on the newly generated region interface between the extended region and the extended region, so as to obtain the numerical model of the extended region.

7. The efficient time-domain domain decomposition algorithm based on the time-domain spectral element method according to claim 1, characterized in that, In step 5, the independent interface electric field region is solved first, and then the internal electric field of each sub-region is calculated; since the electric field transmission condition is introduced at the region interface, the independently solved interface electric field must satisfy the continuity condition of the tangential electric field; then, the obtained interface electric field is used as the excitation source for each region, and further the electric field values inside each region are obtained.

Images