Broadband uniaxial anisotropic complete matching layer absorption boundary method

Abstract

The invention discloses a broadband uniaxial anisotropic complete matching layer absorption boundary method, and belongs to the field of computational electromagnetism. The method comprises the following steps: initializing electromagnetic field component coefficients (Ex, Ey, Hz), an excitation source and W-UPML parameters; constructing a matrix equation through an AH domain differential operator to solve a magnetic field component coefficient, and converting an AH domain field value into a time domain field value in combination with a domain inverse transformation formula; and high-efficiency calculation and boundary absorption of a broadband electromagnetic field are realized by utilizing a Hermite orthogonal basis function and an attenuation factor b. Compared with a traditional method, the time step size stability limitation is avoided through AH domain transformation, and the calculation efficiency is remarkably improved; the low-frequency wave absorption performance is enhanced by the introduced attenuation factor, and the reflection error is reduced by 20dB compared with that of the traditional UPML. The method is suitable for the fields of electromagnetic simulation, seismic wave simulation and the like, and effective technical support is provided for high-precision and broadband electromagnetic field analysis.

Description

Technical Field

[0001] The invention relates to a broadband uniaxial anisotropic perfect matching layer absorption boundary method, which belongs to the field of computational electromagnetics. Background Art

[0002] The finite-difference time-domain (FDTD) method, due to its computational simplicity and ease of implementation, is widely used to address electromagnetic problems in various complex media, providing a powerful analytical tool for electromagnetic theory. However, when simulating electromagnetic problems, the traditional FDTD method is constrained by the Cauchy (Courant Friedrichs Lewy) stability condition for its time and space steps, resulting in reduced computational efficiency. To overcome the time-step limitations imposed by these stability conditions, various unconditionally stable FDTD methods have been proposed. These include the alternating direction implicit finite-difference time-domain (ADI-FDTD) method and the associated Hermite finite-difference time-domain (AH-FDTD) method. Although the ADI-FDTD method achieves unconditional stability, increasing the time step increases its numerical dispersion error, thereby reducing computational accuracy. The AH-FDTD algorithm can not only break through the limitations of stability conditions, but also avoid large dispersion errors when using a large time step. Therefore, the AH-FDTD algorithm can efficiently solve electromagnetic problems in complex media.

[0003] In addition, when using the FDTD method for numerical simulation of electromagnetic wave propagation, due to limited computing resources, FDTD calculations can only be performed in a limited area. In order to simulate free space or an infinite area, it is necessary to set an absorption boundary at the cutoff of the calculation area. The uniaxial anisotropy perfectly matched layer (UPML) absorption boundary is widely used in anisotropic dispersive media and has the advantages of high precision, fast computational efficiency, and simple implementation. However, studies have found that the absorption effect of UPML on low-frequency waves is not ideal. To address this problem, the present invention adopts an improved wideband UPML (Wideband UPML, abbreviated as W-UPML) absorption boundary to achieve a good absorption effect on low-frequency waves.

[0004] Purpose of the Invention

[0005] The purpose of the present invention is to provide a broadband uniaxial anisotropic perfectly matched layer absorbing boundary method. By introducing an attenuation factor, the problem of poor attenuation of low-frequency waves is solved, which has certain guiding significance for the fields of electromagnetic field simulation, seismic wave simulation, etc.

[0006] The purpose of the present invention is achieved through the following technical solutions:

[0007] A broadband uniaxial anisotropic perfectly matched layer absorbing boundary method is proposed, and the specific steps are as follows:

[0008] Step 1: Input parameters into the electromagnetic simulation calculation area model;

[0009] Step 2: Initialize the electromagnetic field component coefficients (E x ,E y ,H z ), current density source term J y 、Hermite polynomial H q (t); and set the parameter σ of the W-UPML absorbing boundary η , κ η , coefficient C η ; Get the interface position η0 of the W-UPML layer close to the simulation area; the thickness d of the W-UPML absorbing boundary; the unit matrix β of the Q order; the differential operator α of the AH domain;

[0010] Step 3: Add an excitation source in the y direction Among them, t c , t d and f c is the excitation source parameter, t c is the pulse center time, t d is the pulse width parameter, f c is the frequency;

[0011] Differential operator α based on AH domain and coefficient C of W-UPML absorbing boundary η , through the matrix equation AH z =B solves the magnetic field component coefficient H of the Associated Hermite (AH) domain z , where A is a banded sparse coefficient matrix and B represents the excitation term;

[0012] Step 4: Change step H z Substitute into the AH domain electric field component equation and calculate the AH domain electric field component coefficient E x 、E y ;

[0013] Step 5: Obtain the time domain field value at the observation point through the domain inverse transformation formula:

[0014]

[0015] Among them, U represents the electromagnetic field component in the time domain, q is the order, and U q (r) represents the electromagnetic field component coefficient in the AH domain, is the q-order AH orthogonal basis function, is a q-order Hermite polynomial, t represents time, Express Perform q-order derivative operation;

[0016] Step 6: The time domain field values ​​obtained by domain inverse transformation are used to further analyze the electromagnetic field waveform and frequency domain characteristics of the observation point, and to evaluate the reflection error and broadband absorption performance of the W-UPML absorption boundary.

[0017] Furthermore, the parameters input in step 1 include:

[0018] Calculate the area size N x ×N y , where N x is the number of grids in the x-direction, N y is the number of grid cells in the y-direction; the spatial steps Δx and Δy, where x is the abscissa and y is the ordinate; the time step Δt; the magnetic permeability μ0 and the dielectric constant ε0 in vacuum; the number of layers NUPML of the absorbing boundary; and the related parameter κ of the W-UPML absorbing boundary. ηmax , σ ηmax , m, where κ ηmax Take an integer, κ ηmax The value range of σ is [1,60], ηmax The value of σ depends on opt , σ opt =(m+1) / (150πΔz),σ ηmax / σ opt The value range of is (0,12], the value range of m is [1,20]; the simulation calculation time T f ; The total order of the Hermite polynomial Q, Q ≥ 0 and is an integer; time scale factor l; observation point; field source parameters.

[0019] Furthermore, the parameter σ of the W-UPML absorbing boundary in step 2 is η , κ η , coefficient C η , specifically:

[0020] σ η =σ ηmax |η-η0| m / d m

[0021] κ η =1+(κηmax -1)|η-η0| m / d m

[0022] Where η = x, y, η0 is the interface position of the W-UPML layer close to the simulation area, σ ηmax is the maximum conductivity parameter of the absorption boundary, κ η is the dielectric constant adjustment parameter of the absorbing edge, κ ηmax is κ η The maximum value of m is the polynomial order that controls the spatial distribution of the absorption boundary parameters, d is the thickness of the W-UPML absorption boundary, and the value range of Δη is (0,λ12], where λ is the wavelength of the electromagnetic wave.

[0023] C η =k η β+σ η (bβ+ε0α) -1

[0024] Where η = x, y, β is the Q-order identity matrix, α is the differential operator in the AH domain, and b is the attenuation factor.

[0025] Furthermore, the differential operator α in the AH domain is:

[0026]

[0027] Wherein, Q ≥ 0 and is an integer.

[0028] Furthermore, the attenuation factor b gradually decays from a maximum value to zero within the boundary, so that the absorption efficiency of the absorption boundary for low-frequency electromagnetic waves is higher than that of a conventional UPML boundary.

[0029] Furthermore, the matrix equation AH in step 3 is z =B is solved by LU decomposition and pursuit method.

[0030] Furthermore, the AH domain electric field component coefficient E in step 4 x 、E y satisfy:

[0031]

[0032] Where Δx and Δy represent the distances between adjacent computational grids in the x and y directions, respectively; i represents the i-th computational grid in the x direction; and j represents the j-th computational grid in the y direction. x 、C y is the coefficient of the W-UPML absorbing boundary; E x | i,j 、E y | i,j 、Hz | i,j are the expansion coefficients of the x-direction electric field component, the y-direction electric field component, and the z-direction magnetic field component in the AH domain at the grid point (i, j), respectively. y | i,j represents the expansion coefficient of the excitation source in the y direction at the grid point (i, j) in the AH domain.

[0033] The time domain field values ​​obtained by the domain inverse transformation in the present invention can be used for: (1) analyzing the electromagnetic field waveform and frequency domain characteristics of the observation point to verify the simulation accuracy; (2) evaluating the reflection error and broadband absorption performance of the W-UPML absorption boundary. Compared with the traditional UPML, the reflection error of the W-UPML absorption boundary is reduced by 20dB; (3) providing data support for electromagnetic field simulation, seismic wave simulation, etc., breaking through the stability limitations of the traditional FDTD method and improving computational efficiency.

[0034] Technical Effects

[0035] Compared with the prior art, the present invention has the following beneficial effects:

[0036] 1) The present invention provides a method for realizing a broadband uniaxial anisotropic perfectly matched layer absorption boundary. This method transforms the Maxwell equations into the AH domain for calculation. This method does not involve a time step when calculating the electromagnetic field component coefficients in the entire calculation region, thus breaking through the limitation of stability conditions and improving calculation efficiency.

[0037] 2) The present invention provides a method for realizing a broadband uniaxial anisotropic perfectly matched layer absorption boundary. Due to the introduction of the attenuation factor b, the method is more effective in absorbing low-frequency waves compared with the conventional UPML absorption boundary.

[0038] In summary, the time domain field values ​​obtained by the domain inverse transformation of the present invention can be used for: (1) analyzing the electromagnetic field waveform and frequency domain characteristics of the observation point to verify the simulation accuracy; (2) evaluating the reflection error and broadband absorption performance of the W-UPML absorption boundary. Compared with the traditional UPML, the reflection error of the W-UPML absorption boundary is reduced by 20dB; (3) providing data support for electromagnetic field simulation, seismic wave simulation, etc., breaking through the stability limitations of the traditional FDTD method and improving computational efficiency. BRIEF DESCRIPTION OF THE DRAWINGS

[0039] Figure 1 It is a schematic flow diagram of the present invention;

[0040] Figure 2 is a schematic diagram of a calculation model in an embodiment of the present invention;

[0041] Figure 3 This is a comparison diagram of the time domain magnetic field waveform at the observation point between the method of the present invention and the traditional FDTD method;

[0042] Figure 4 is a graph of relative reflection errors of different absorption boundaries at observation points in an embodiment of the present invention;

[0043] Figure 5 2 is a comparison of the L2 errors of the absorption boundary at different frequencies in an embodiment of the present invention. DETAILED DESCRIPTION

[0044] The present invention will be described in further detail below with reference to the accompanying drawings and specific embodiments. Figure 1 As shown, the steps are as follows:

[0045] Step 1: Input model file;

[0046] Step 2: Initialize parameters and set parameters;

[0047] Step 3: Add an excitation source in the y direction and calculate the magnetic field component coefficient H of the Associated Hermite (AH) domain z ;

[0048] Step 4: Calculate the electric field component coefficient E in the AH domain x 、E y ;

[0049] Step 5: Obtain the time domain field value at the observation point through domain inverse transformation.

[0050] Specifically, the Maxwell equations for the W-UPML absorbing boundary are expressed as:

[0051]

[0052] in, represents the electric field vector, represents the magnetic field vector, j is the imaginary unit, ω is the angular frequency, ε0 ​​is the vacuum dielectric constant, μ0 is the vacuum permeability, is the differential operator, is the anisotropic feature matching matrix, which can be written as:

[0053]

[0054] Among them, s x 、s y are the relative permittivity tensor and relative permeability tensor of the anisotropic medium in the x-direction and y-direction, respectively, expressed as:

[0055]

[0056] Where η = x, y, κ η , σ η is the parameter of W-UPML, and b is the attenuation factor.

[0057] The present invention only considers the two-dimensional TE in simple lossless media Z In the wave case, equation (1) can be written as:

[0058]

[0059] Among them, E x 、E y represent the electric field in the x and y directions respectively, H z represents the magnetic field in the z direction.

[0060] According to the operation rules of AH linear operator Equations (4)-(6) can be directly converted to the AH domain

[0061]

[0062] Where i represents the i-th computational grid on the horizontal axis, j represents the j-th computational grid on the vertical axis, Δη (η=x,y) represents the distance between discrete points in space, C η is the W-UPML coefficient related to the coordinate grid and is calculated as:

[0063] C η =κ η β+σ η (bβ+αε0) -1 (10)

[0064] Among them, β is the Q-order identity matrix, Q is the total order of the Hermite polynomial, and α is the differential operator in the AH domain, which is specifically expressed as:

[0065]

[0066] Substituting (7) and (8) into (9) and eliminating the electric field component, we can obtain the five-point equation related only to the magnetic field:

[0067]

[0068] (12) is written in the form of a matrix equation, specifically:

[0069] AH z =B (13)

[0070] Where A is a banded sparse coefficient matrix, H z is the unknown vector to be solved, consisting of the magnetic field coefficients of all orders at all points, and B represents the excitation term. The matrix equation is decomposed into LU, and the magnetic field coefficients in the AH domain are solved using the pursuit method. Substituting the magnetic field coefficients into (7) and (8) yields the electric field expansion coefficients in the AH domain.

[0071] Finally, the time domain results can be reconstructed through the AH domain electromagnetic field expansion coefficients and basis functions. The specific reconstruction formula is as follows:

[0072]

[0073] Among them, U represents the electromagnetic field component in the time domain, U q (r) represents the electromagnetic field component coefficient in the AH domain, is the q-order AH orthogonal basis function, is a Hermite polynomial of order q.

[0074] The effect of the present invention is described below by experiment:

[0075] Experiment: Simulation calculation of parallel plate waveguide.

[0076] According to the method steps of the present invention, Figure 2 As shown in the figure, the entire calculation area in the experiment is a 700×700 grid with a grid size of 0.015m×0.015m, that is, Δx=Δy=0.015m. A cosine modulated Gaussian pulse is selected as the excitation source in the y direction, and the specific expression is:

[0077]

[0078] Among them, t d =1 / (2f c ), t c =4t d , f c =1GHz. The observation point A is set at the 2 grid positions after the W-UPML layer. The time step Δt = 250ps, the order of the Hermite polynomial Q = 60, and the time scale factor l = 5.12×10 -10 The entire simulation time is T f =10.81ns. W-UPML absorbing edge parameter κ ηmax =5,σ ηmax =σ opt , m=4,NUPML=8,the attenuation factor gradually decays from the maximum value to zero within the boundary, and the maximum value is b max =0.4. The magnetic field waveform at the observation point calculated by the method of the present invention is shown in Figure 3 .from Figure 3 It can be seen that the calculation results of the method of the present invention are consistent with those of the FDTD method, which verifies the correctness of the method of the present invention. Figure 4 is the relative reflection error of different absorption boundaries at the observation point, and the calculation formula is specifically expressed as:

[0079]

[0080] Among them, H z (t) is the magnetic field value at the observation point under different boundary conditions, is the reference magnetic field value, is the maximum absolute value of the reference magnetic field. Figure 4 It can be seen that the maximum reflection error of the W-UPML absorption boundary is lower than -65dB, which is about 20dB lower than the reflection error of the conventional UPML absorption boundary, indicating that the W-UPML absorption boundary has better absorption performance. The attenuation factor b can greatly improve the performance of the absorption boundary. Figure 5 To compare the L2 error of the absorbing boundary at different frequencies, the calculation formula of the L2 error is as follows:

[0081]

[0082] Depend on Figure 5 It can be seen that as the frequency decreases, the accuracy of the W-UPML absorption boundary is higher than that of the conventional UPML absorption boundary, which shows that the W-UPML absorption boundary can effectively improve the problem of poor attenuation of low-frequency waves.

[0083] The foregoing description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Those skilled in the art will readily appreciate that various modifications and variations of the present invention are possible. Any modifications, equivalent substitutions, or improvements made within the spirit and principles of the present invention are intended to be within the scope of protection of the present invention.

Claims

1. A broadband uniaxial anisotropic perfectly matched layer absorbing boundary method, characterized by: The specific steps are as follows: Step 1: Input parameters into the electromagnetic simulation calculation area model; Step 2: Initialize the electromagnetic field component coefficients (E x ,E y ,H z ), current density source term J y 、Hermite polynomial H q (t); and set the parameter σ of the W-UPML absorbing boundary η , κ η , coefficient C η ; Get the interface position η0 of the W-UPML layer close to the simulation area; the thickness d of the W-UPML absorbing boundary; the unit matrix β of the Q order; the differential operator α of the AH domain; Step 3: Add an excitation source in the y direction Among them, t c , t d and f c is the excitation source parameter, t c is the pulse center time, t d is the pulse width parameter, f c is the frequency; Differential operator α based on AH domain and coefficient C of W-UPML absorbing boundary η , through the matrix equation AH z =B to solve for the magnetic field component coefficient H z , where A is a banded sparse coefficient matrix and B represents the excitation term; Step 4: Change step H z Substitute into the AH domain electric field component equation and calculate the AH domain electric field component coefficient E x 、E y ; Step 5: Obtain the time domain field value at the observation point through the domain inverse transformation formula: Among them, U represents the electromagnetic field component in the time domain, q is the order, and U q (r) represents the electromagnetic field component coefficient in the AH domain, is the q-order AH orthogonal basis function, is a q-order Hermite polynomial, t represents time, Express Perform q-order derivative operation; Step 6: The time domain field values ​​obtained by domain inverse transformation are used to further analyze the electromagnetic field waveform and frequency domain characteristics of the observation point, and to evaluate the reflection error and broadband absorption performance of the W-UPML absorption boundary.

2. The broadband uniaxial anisotropic perfectly matched layer absorbing boundary method according to claim 1, characterized in that: The parameters input in step 1 include: Calculate the area size N x ×N y , where N x is the number of grids in the x-direction, N y is the number of grid cells in the y-direction; the spatial steps Δx and Δy, where x is the abscissa and y is the ordinate; the time step Δt; the magnetic permeability μ0 and the dielectric constant ε0 in vacuum; the number of layers NUPML of the absorbing boundary; and the related parameter κ of the W-UPML absorbing boundary. ηmax , σ ηmax , m, where κ ηmax Take an integer, κ ηmax The value range of σ is [1,60], ηmax The value of σ depends on opt ,σ opt =(m+1) / (150πΔz),σ ηmax / σ opt The value range of is (0,12], the value range of m is [1,20]; the simulation calculation time T f ; The total order of the Hermite polynomial Q, Q ≥ 0 and is an integer; time scale factor l; observation point; field source parameters.

3. The broadband uniaxial anisotropic perfectly matched layer absorbing boundary method according to claim 1, characterized in that: The parameter σ of the W-UPML absorbing boundary in step 2 is η , κ η , coefficient C η , specifically: s η =s ηmax |n-n0| m / d m k η =1+(k ηmax -1)|η-η0| m / d m Where η = x, y, η0 is the interface position of the W-UPML layer close to the simulation area, σ ηmax is the maximum conductivity parameter of the absorption boundary, κ η is the dielectric constant adjustment parameter of the absorbing edge, κ ηmax is κ η The maximum value of m is the polynomial order that controls the spatial distribution of the absorption boundary parameters, d is the thickness of the W-UPML absorption boundary, and the value range of Δη is (0,λ12], where λ is the wavelength of the electromagnetic wave. C η =k η b+s η (bβ+ε0α) -1 Where η = x, y, β is the Q-order identity matrix, α is the differential operator in the AH domain, and b is the attenuation factor.

4. The broadband uniaxial anisotropic perfectly matched layer absorbing boundary method according to claim 3, characterized in that: The differential operator α in the AH domain is: Wherein, Q ≥ 0 and is an integer.

5. The broadband uniaxial anisotropic perfectly matched layer absorbing boundary method according to claim 3, characterized in that: The attenuation factor b gradually decays from a maximum value to zero within the boundary, so that the absorption efficiency of the absorption boundary for low-frequency electromagnetic waves is higher than that of a conventional UPML boundary.

6. The broadband uniaxial anisotropic perfectly matched layer absorbing boundary method according to claim 1, characterized in that: The matrix equation AH in step 3 z =B is solved by LU decomposition and pursuit method.

7. The broadband uniaxial anisotropic perfectly matched layer absorbing boundary method according to claim 1, characterized in that: The electric field component coefficient E in the AH domain in step 4 x 、E y satisfy: Where Δx and Δy represent the distances between adjacent computational grids in the x and y directions, respectively; i represents the i-th computational grid in the x direction; and j represents the j-th computational grid in the y direction. x 、C y is the coefficient of the W-UPML absorbing boundary; E x | i,j 、E y | i,j 、H z | i,j are the expansion coefficients of the x-direction electric field component, the y-direction electric field component, and the z-direction magnetic field component in the AH domain at the grid point (i, j), respectively. y | i,j represents the expansion coefficient of the excitation source in the y direction at the grid point (i, j) in the AH domain.

Images