Universal low-numerical dispersion finite difference time domain method

Abstract

The invention relates to a universal low numerical dispersion time domain finite difference method, which comprises the following steps of: based on a physical model structure of electromagnetic simulation engineering, selecting a time domain finite difference FDTD implicit or explicit subclass method which needs to be applied, and solving an optimal weighting factor and a correction factor; a stability CFL condition is verified; reconstructing an iterative formula that the electric field and the magnetic field advance along with time steps; the calculation main body carries out calculation solution, and the calculation main body involving implicit solution carries out calculation by using a tridiagonal matrix; on the basis of the steps, electric field parameters and magnetic field parameters in the whole area under a specific simulation instance are obtained, and therefore expected physical parameters are obtained through calculation. The isotropic sampling scheme is carried out in the stability condition limitation direction, the traditional sampling scheme is still kept in the remaining direction, and integration is carried out through a weighting factor and a correction factor. Wherein the values of the weighting factor and the correction factor are in one-to-one correspondence with the dispersion characteristics, so that the algorithm is ensured to have accurate parameter adaptability.

Description

Technical Field

[0001] The present invention relates to a universal low numerical dispersion finite-difference time-domain method, belonging to the technical field of electromagnetic simulation technology. Background Art

[0002] Since K.Y. Yee first proposed the finite-difference time-domain (FDTD) method in 1966, it has stood out among computational electromagnetics algorithms due to its significant advantages, including efficient time-domain computation and flexible boundary handling. In recent years, advances in computing power have broadened the application of the FDTD method to specific electromagnetic engineering problems, such as antenna radiation, resonant cavities, and RCS.

[0003] However, the analysis of fine structures usually requires a finer grid, which means that under the constraints of stability conditions, an extremely small time step must be set, which will inevitably lead to a decrease in computational efficiency. Therefore, in order to improve computational efficiency, some scholars have proposed ADI-FDTD, CN-FDTD, HIE-FDTD, WCS-FDTD, and SS-FDTD methods. These methods can weaken or even completely eliminate the restriction of the unit space step to the unit time step in the target direction through implicit operations during the calculation process, but they also lead to greater dispersion errors (the numerical phase velocity deviates from the exact phase velocity), which is very unfavorable for dealing with some phase-sensitive electromagnetic problems.

[0004] ADI-FDTD

[0005] This algorithm uses implicit differencing to relax the time step limit to several or even dozens of times that of explicit methods, significantly reducing the number of iterations. It is particularly suitable for large-scale problems (such as microwave devices and complex target scattering). However, due to its relaxed stability conditions, this method also has the worst dispersion error.

[0006] HIE-FDTD

[0007] This method is implemented by combining explicit and implicit difference schemes in different directions. Typically, an implicit scheme is used in directions where a fine mesh is required, while an explicit scheme is used in other directions. This has the advantage that the time step is no longer solely determined by the minimum mesh size in the computational domain, thereby relaxing the Courant stability condition and improving computational efficiency. For example, for electromagnetic models with fine structures in a certain direction, such as antennas with narrow slots or thin wires, the HIE-FDTD method can significantly reduce computational time while maintaining accuracy.

[0008] AA-HIE-FDTD

[0009] The AA method (artificial anisotropic parameter method) modifies the dielectric constant and permeability to make the speed of light in numerical simulations equal to the physical speed. Although this method can significantly reduce numerical dispersion, it still has anisotropic characteristics.

[0010] ID-HIE-FDTD

[0011] This method focuses on sampling and achieves isotropic dispersion through weighted approaches. Some researchers have incorporated the ID scheme into HIE-FDTD and ADI-FDTD algorithms, applying it to three different approaches. While this approach significantly reduces dispersion errors, it also increases computational complexity and reduces efficiency, hindering the rapid processing of complex electromagnetic problems.

[0012] In this context, the present invention proposes a new universal low numerical dispersion finite-difference time-domain method and system, which can be combined with various subclasses of FDTD algorithms to significantly reduce numerical dispersion and improve calculation accuracy. Summary of the Invention

[0013] Technical Problem: This paper addresses the dispersion error characteristics of weakly conditionally stable FDTD methods and proposes a universal low-numerical-dispersion finite-difference time-domain method. This method employs an isotropic sampling scheme in the directions constrained by the stability condition, while maintaining a traditional sampling scheme in the remaining directions. This scheme is then integrated using weighting and correction factors. The values ​​of these weighting and correction factors correspond one-to-one with the dispersion characteristics, ensuring precise parameter adaptability.

[0014] Technical solution: A universal low numerical dispersion finite-difference time-domain method of the present invention comprises the following steps:

[0015] Step 1. Based on the physical model structure of the electromagnetic simulation project, analyze whether it has a fine structure direction, and then select the implicit or explicit subclass method of finite-difference time-domain (FDTD) to be applied;

[0016] Step 2. Based on the selected finite-difference time-domain (FDTD) method, write the dispersion equation derived from the plane wave solution, apply the isotropic sampling scheme to the selected direction, and obtain the optimal weighting factor and correction factor; and verify the stability CFL condition;

[0017] Step 3. Based on the obtained optimal weighting factor, reconstruct the iterative formula of the electric field and magnetic field as the time step advances;

[0018] Step 4. Perform calculations and solve the problem based on the newly obtained iterative formula. The calculations involving implicit solutions are calculated using tridiagonal matrices.

[0019] Step 5. Based on the above steps, the electric field parameters and magnetic field parameters in the entire area under the specific simulation example are obtained, so as to obtain the desired physical parameters through calculation.

[0020] in,

[0021] In step 1, the selection of the implicit, explicit, or hybrid explicit-implicit subclass method of the finite-difference time-domain FDTD to be applied is specifically as follows: for ordinary objects, the ordinary FDTD method is selected; for methods with fine structures in one direction, the hybrid explicit-implicit finite-difference time-domain method HIE-FDTD method is adopted; for methods with fine structures in two directions, the weak conditionally stable finite-difference time-domain method WCS-FDTD is selected.

[0022] In step 1, the fine structure direction is determined by the meshing method based on the fine structure direction of the specific physical model. Assuming that the object has a fine structure in the y direction, the y direction is called the implicit direction, and the y-direction mesh size Δy is irrelevant to the stability condition. The x-direction mesh size Δx and the z-direction mesh size Δz do not participate in the stability condition.

[0023] The dispersion equation described in step 2 is the dispersion equation for one-step HIE:

[0024]

[0025] in:

[0026] k x 、k y 、k z are the numerical wave number components in the three directions in the rectangular coordinate system, where k z = kcosθ. k is the numerical wave number during simulation; where θ, are the elevation angle and azimuth angle in the spherical coordinate system,

[0027] c is the speed of light in vacuum, Δt is the time step of calculation, ω is the angular frequency of the wave; J x 、J y 、J z It is to simplify the algebraic expressions in the x, y, and z directions;

[0028] If the isotropic sampling ID sampling scheme is applied only to the x and z directions, then J in the dispersion equation x and J y Make the following corrections:

[0029]

[0030] The superscript ID indicates that an isotropic sampling scheme is applied in this direction, and α is the weighting factor.

[0031] Step 3: Based on the obtained optimal weighting factor, the iterative formula for reconstructing the electric and magnetic fields as time steps is specifically determined by applying the ID sampling scheme in the x and z directions. The maximum normalized phase velocity of the dispersion equation of the one-step HIE-FDTD method is R max , the minimum value is R min , let w=R max -R min， By traversing the parameters in MATLAB, we find the value of α when w is the minimum, which is called the optimal weighting factor. The dispersion equation of the universal low numerical dispersion finite-difference time-domain method is applied to the HIE case:

[0032]

[0033] The superscript ID indicates that an isotropic sampling scheme is applied in this direction.

[0034] In order to obtain the optimal weighting factor and correction factor, select As a reference, the dispersion equation is used, where θ, They are the elevation angle and azimuth angle in the spherical coordinate system, and the precise wave number k is obtained. exact The correction factor sc is determined:

[0035]

[0036] The correction factor sc acts on the dielectric constant and permeability ε′=sc×ε, μ′=sc×μ to obtain a dispersion curve with isotropy and almost close to the exact phase velocity; ε′ is the corrected dielectric constant, ε is the dielectric constant before correction, μ′ is the corrected permeability, and μ is the permeability before correction.

[0037] Step 4: Calculate and solve the new iterative formula based on the main body, and analyze the iterative formula of the main body, that is, the electric field and magnetic field changing with time. The calculation of the three-dimensional example includes the electric field value E in the x, y, and z directions. x 、E y 、E z and magnetic field value H x 、H y 、H z Applying the isotropic sampling ID scheme to the differences in the x and z directions yields:

[0038]

[0039]

[0040] The superscript n represents a specific time point in the simulation process, ε is the dielectric constant before correction, and μ is the magnetic permeability before correction.

[0041] Step 5: Based on the above steps, obtain the electric field parameters and magnetic field parameters in the entire area under the specific simulation example. After selecting the observation point or observation plane, obtain the desired physical parameters including but not limited to the transmission coefficient by calculation. Reflection coefficient Normalized phase velocity v p =ω / k.

[0042] Beneficial Effects: Under the same grid parameters, this method significantly outperforms traditional methods in computational accuracy, particularly when dealing with complex phase-sensitive problems such as electromagnetic scattering. The study systematically derived a dispersion formula for applying the universal low-dispersion finite-difference time-domain method to the HIE scenario, and fully verified its theoretical feasibility and practical effectiveness through simulation comparisons with traditional HIE schemes.

[0043] Further research revealed that the numerical dispersion error of HIE-FDTD in the fine structure direction remains constant, regardless of whether ID processing is used. Compared to conventional HIE methods that perform ID processing in all three directions, the proposed method significantly reduces the additional computational overhead associated with ID processing while maintaining dispersion performance, effectively simplifying the computational process and significantly improving computational efficiency.

[0044] To address the inherent dispersion error flaws of the FDTD algorithm, this paper proposes a novel, universal, low-numerical-dispersion finite-difference time-domain (FDTD) method and system. This system achieves isotropic dispersion by applying an ID sampling scheme in local (explicit) directions, significantly reducing dispersion error. Compared to existing technologies, this method achieves superior computational accuracy and dispersion performance while satisfying the CFL condition.

[0045] Compared to existing low-numerical-dispersion FDTD techniques, this method maintains the tridiagonal matrix solution form during implicit iterative solutions, resulting in the lowest computational complexity compared to other implicit ID methods. Comparison with FDTD methods reveals that the proposed method exhibits excellent isotropic dispersion characteristics and a more relaxed CFL condition. Three numerical examples demonstrate the advantages of high accuracy and computational efficiency, demonstrating significant advantages and promising applications in solving phase-sensitive problems. BRIEF DESCRIPTION OF THE DRAWINGS

[0046] Figure 1 This is a schematic diagram of dispersion error comparison.

[0047] Figure 2This is a schematic diagram comparing the three-dimensional dispersion performance of the three methods.

[0048] Figure 3 It is a calculation model of transmission and reflection coefficient.

[0049] Figure 4 are the transmission coefficient and the reflection coefficient. Figure 4 where a is the transmission coefficient of the periodic grating, Figure 4 Where b is the reflection coefficient of the periodic grating.

[0050] Figure 5 is a spherical scattering model.

[0051] Figure 6 It is the real part of the field value at the observation point and the imaginary part of the output value at the observation point. Figure 6 where a is the real part of the electric field value at the observation point, Figure 6 The b in it is the imaginary part of the electric field value at the observation point.

[0052] Figure 7 It is a flow chart of the invention implementation. DETAILED DESCRIPTION

[0053] The technical solutions in the embodiments of the present invention will be described clearly and completely below with reference to the accompanying drawings in the embodiments of the present invention.

[0054] like Figure 7 A universal low numerical dispersion finite-difference time-domain method of the present invention comprises the following steps:

[0055] Step 1. Based on the physical model structure of the electromagnetic simulation project, analyze whether it has a fine structure direction, and then select the implicit or explicit subclass method of finite-difference time-domain (FDTD) to be applied;

[0056] Step 2. Based on the selected finite-difference time-domain (FDTD) method, write the dispersion equation derived from the plane wave solution, apply the isotropic sampling scheme to the selected direction, and obtain the optimal weighting factor and correction factor; and verify the stability CFL condition;

[0057] Step 3. Based on the obtained optimal weighting factor, reconstruct the iterative formula of the electric field and magnetic field as the time step advances;

[0058] Step 4. Perform calculations and solve the problem based on the newly obtained iterative formula. The calculations involving implicit solutions are calculated using tridiagonal matrices.

[0059] Step 5. Based on the above steps, the electric field parameters and magnetic field parameters in the entire area under the specific simulation example are obtained, so as to obtain the desired physical parameters through calculation.

[0060] in,

[0061] The selection of the implicit, explicit, or hybrid implicit subclass method of the finite-difference time-domain FDTD method to be applied is specifically as follows: for common objects, the ordinary FDTD method is selected; for methods with fine structures in one direction, the hybrid implicit and explicit finite-difference time-domain method HIE-FDTD method is adopted; for methods with fine structures in two directions, the weak conditionally stable finite-difference time-domain method WCS-FDTD is selected.

[0062] The dispersion equation is a one-step HIE dispersion equation: preferably, the FDTD method to be applied is selected (particularly, the present invention takes the traditional one-step HIE-FDTD as an example).

[0063] Preferably, the grid division method is determined according to the fine structure direction of the specific physical model. Assuming that the object has a fine structure in the y direction, the y direction is called the implicit direction, and the y-direction grid size Δy is irrelevant to the stability condition; the x-direction grid size Δx and the z-direction grid size Δz do not participate in the stability condition.

[0064] Preferably, the dispersion equation of the original solution is written out. The dispersion equation of the one-step HIE is:

[0065]

[0066] in:

[0067] where k x 、k y 、k z is the numerical wave number in the three directions in the rectangular coordinate system, c is the speed of light in vacuum, Δt is the time step of the calculation, and ω is the angular frequency of the wave.

[0068] Step 4. Apply the ID sampling scheme only in the x and z directions, then J in the dispersion equation x and J y Make the following corrections:

[0069]

[0070] The superscript ID indicates that an isotropic sampling scheme is applied in this direction, and α is the weighting factor.

[0071] Preferably, the maximum normalized phase velocity of the dispersion equation obtained by applying the universal low numerical dispersion finite-difference time-domain method to the HIE case is defined as R max , the minimum value is R min , let w=R max -R min， The α obtained by traversing the parameters through MATLAB to find the minimum value of w is called the optimal weighting factor.

[0072] Preferably, select As a reference, the dispersion equation at the time of θ can be used to obtain the exact wave number k exact The correction factor sc is determined:

[0073]

[0074] Preferably, applying the correction factor to the dielectric constant and the permeability ε′=sc×ε, μ′=sc×μ can obtain a dispersion curve with isotropy and almost close to the exact phase velocity.

[0075] Preferably, the calculation subject is analyzed, and the calculation of the three-dimensional example includes the electric field value E in the x, y, and z directions. x 、E y 、E z and magnetic field value H x 、H y 、H z Applying the ID scheme to the differences in the x and z directions yields:

[0076]

[0077] The superscripts indicate specific time points, for example, E x n+1 Represents the electric field value in the x direction at time n+1.

[0078] Preferably, the parts involving implicit solution are solved according to the tridiagonal matrix idea, and the variables for explicit solution can be directly iterated.

[0079] Numerical stability (CFL condition) (Step 2. Based on the selected finite-difference time-domain (FDTD) method, write the dispersion equation derived from the plane wave solution, apply the isotropic sampling scheme to the selected direction, and obtain the optimal weighting factors and correction factors; and verify the stability CFL condition.)

[0080] Through matrix analysis, we can determine whether the stability condition of the new algorithm will change after processing the x and z directions.

[0081] The calculation process of one-step HIE-FDTD is expressed as follows:

[0082] Au n+1 =Bu n+1 ,

[0083] where u=[E x ,E y ,E z ,H x ,H y ,H z ].

[0084]

[0085] in:

[0086]

[0087] Let the growth factor be p. Then, in order to satisfy the stability condition, det|Ap-B|=0. Solve the resulting equation as a quartic equation with variable p:

[0088] p1=p2=1;

[0089]

[0090] in:

[0091] M=4ab+Dy 2 +2Gx 2 +2Gz 2 ;

[0092] N=8(Dy 2 +Gx 2 +Gz 2 )(Gx 2 +Gz 2 +4ab)) / 16;

[0093]

[0094] To ensure the stability of the algorithm, |p|<1 is required, so the polynomial in the square root must be negative. It is observed that the numerators of the two unknown solutions only differ in their signs, M+N<0; M+N<0, the sum of the two gives M<0.

[0095] The solution is:

[0096]

[0097] Dispersion analysis (Step 2. Based on the selected finite-difference time-domain (FDTD) method, write the dispersion equation derived from the plane wave solution, apply the isotropic sampling scheme to the selected direction, and obtain the optimal weighting factors and correction factors; and verify the stability CFL condition.)

[0098] In fact, regardless of the FDTD method, the sine function term containing the spatial step in the dispersion equation is obtained by the interaction of the plane wave function and the derivative operation. In the HIE case, since the y direction is implicit, only the x and z terms in the dispersion equation need to be changed; therefore, the Euler formula is still used to process the displacement operator, which will be obtained:

[0099]

[0100] Here D xis a formal notation that represents the partial derivative operation in the x direction. Similarly, performing the same operation in the z direction gives:

[0101]

[0102] The above derivation reveals that the shift operator is equivalent to multiplying the cosine function when applied to the difference. It is convenient to deduce that the dispersion equation applied to HIE-FDTD by this method has a similar form to that of the traditional HIE-FDTD, and only the following transformation is required:

[0103]

[0104] Where J represents the trigonometric function term generated by the difference in the dispersion equation; when α is 0, the dispersion equation degenerates into the traditional form;

[0105] Considering that the purpose of α is to achieve isotropic dispersion, the maximum normalized phase velocity of the dispersion equation is defined as R by applying the universal low numerical dispersion finite-difference time-domain method to the HIE case. max , the minimum value is R min , let w=R max -R min ; Obviously, the value of w changes with the change of α, and the α that makes w take the minimum value is the optimal weighting factor.

[0106] Iterative formula for the low numerical dispersion finite-difference time-domain method. (Step 4. Perform computational solution based on the newly obtained iterative formula. The computational body involving implicit solution uses tridiagonal matrix calculations.)

[0107] In the traditional FDTD method, the magnetic field and electric field surround each other and satisfy the law of electromagnetic induction. We then perform approximate differentiation through central difference. To facilitate the subsequent representation of the ID operation, we first define the shift operator as follows:

[0108] Hu + F(...,u,...)=F(...,u+Δu / 2,...)

[0109] Hu - F(...,u,...)=F(...,u-Δu / 2,...)

[0110] H is the shift operator, F represents the scalar electric field or scalar magnetic field, the superscripts + and - represent positive and negative shifts, and u represents the three coordinate components x, y, and z in the direct coordinate system. Next, the shift operator is used to represent the process of replacing derivatives with differences:

[0111]

[0112] Here, the superscript (1) indicates that the sampling scheme is the original FDTD scheme;

[0113] The spatial derivative of the new scheme is a weighted combination of the original FDTD scheme and the isotropic sampling scheme, and is of the following form:

[0114]

[0115] in, represents the derivative in the u direction, the superscript (2) represents the vertex scheme of one of the ID schemes, and α is the weighting coefficient when using the vertex scheme;

[0116] In the x-direction, the relationship between the vertex solution and the original solution has rotation symmetry and is expressed as follows:

[0117]

[0118] To simplify the expression, a new shift operator is defined as follows:

[0119]

[0120]

[0121] Finally, the differential operation of the new algorithm will consist of three parts:

[0122]

[0123] The spatial coordinates for the difference are (i, j, k). The three rows from top to bottom in the above formula represent the original coordinate points, the one-shifted points, and the two-shifted points respectively; the one-shifted and two-shifted points refer to the number of shift operators. Overall, isotropic sampling is achieved by expanding the sampling points involved in the difference.

[0124] Applying the above transformation to HIE-FDTD, that is, performing ID processing on the spatial differences in the x and z directions, but not on the spatial differences in the y direction, the updated equation of the electromagnetic field is obtained as follows:

[0125]

[0126] This method still maintains the tridiagonal matrix form when solving implicitly, and can be solved iteratively in a simple manner.

[0127] The purpose of the present invention is to provide a universal low numerical dispersion finite-difference time-domain method to achieve efficient and accurate simulation analysis with low numerical dispersion performance.

[0128] Example:

[0129] In order to make the above-mentioned objects, features and advantages of the present invention more clear and understandable, the present invention is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0130] First, the traditional one-step HIE-FDTD method is selected, and then its dispersion equation is given:

[0131]

[0132] in:

[0133] where k x 、k y 、k z is the numerical wave number in the three directions of the rectangular coordinate system, c is the speed of light in vacuum, Δt is the time step of calculation, and ω is the angular frequency of the wave; J x 、J y 、J z In order to simplify the algebraic expressions in the x, y, and z directions, the specific parameters are set as follows: Δx = Δz = 5Δy = λ / 10 = 1 cm, λ is the wavelength, the time step and frequency are Δt = Δy / 2c = 3.33 ps and the source frequency f = 3 GHz; ω = 2πf = 1.885×10 10 rad / s. Numerical wave number k=62.83rad / m,k x 、k y 、k z Parameter traversal will be performed.

[0134] If the isotropic sampling ID sampling scheme is applied only to the x and z directions, then J in the dispersion equation x and J y Make the following corrections:

[0135]

[0136] The superscript ID indicates that an isotropic sampling scheme is applied in this direction, and α is the weighting factor;

[0137] The maximum normalized phase velocity of the dispersion equation in the HIE case is defined as R max , the minimum value is R min , let w=R max -R min， The α obtained by traversing the parameters in MATLAB when w is the minimum is called the optimal weighting factor. In this example, the α found is 0.169, and the minimum error w is 6.7×10 -4 The dispersion equation at this time is:

[0138]

[0139] Select Among them, They are the elevation angle and azimuth angle in the spherical coordinate system respectively, and the dispersion equation at t is used as a reference to obtain the precise wave number k exact The correction factor sc is determined:

[0140]

[0141] Applying the correction factor to the dielectric constant and permeability ε′=sc×ε, μ′=sc×μ, we can obtain a dispersion curve with isotropy and almost accurate phase velocity; ε′ is the corrected dielectric constant, ε is the dielectric constant before correction, μ′ is the corrected permeability, and μ is the permeability before correction.

[0142] Through matrix analysis, we can determine whether the stability condition of the new algorithm will change after processing the x and z directions.

[0143] The calculation process of one-step HIE-FDTD is expressed as follows:

[0144] Au n+1 =Bu n+1 ,

[0145] where u=[E x ,E y ,E z ,H x ,H y ,H z ].

[0146]

[0147] in:

[0148]

[0149] Let the growth factor be p. Then, in order to satisfy the stability condition, det|Ap-B|=0. Solve the resulting equation as a quartic equation with variable p:

[0150] p1=p2=1;

[0151]

[0152] in:

[0153] M=4ab+Dy 2 +2Gx 2 +2Gz 2 ;

[0154] N=8(Dy 2 +Gx 2 +Gz 2)(Gx 2 +Gz 2 +4ab)) / 16;

[0155]

[0156] To ensure the stability of the algorithm, |p|<1 is required, so the polynomial in the square root must be negative. It is observed that the numerators of the two unknown solutions only differ in their signs, M+N<0; M+N<0, the sum of the two gives M<0.

[0157] The solution is:

[0158]

[0159] In the traditional FDTD method, the magnetic field and electric field surround each other and satisfy the law of electromagnetic induction. We then perform approximate differentiation through central difference. To facilitate the subsequent representation of the ID operation, we first define the shift operator as follows:

[0160] Hu + F(...,u,...)=F(...,u+Δu / 2,...)

[0161] Hu - F(...,u,...)=F(...,u-Δu / 2,...)

[0162] H is the shift operator, F represents the scalar electric field or scalar magnetic field, the superscripts + and - represent positive and negative shifts, and u represents the three coordinate components x, y, and z in the direct coordinate system. Next, the shift operator is used to represent the process of replacing derivatives with differences:

[0163]

[0164] Here, the superscript (1) indicates that the sampling scheme is the original FDTD scheme;

[0165] The spatial derivative of the new scheme is a weighted combination of the original FDTD scheme and the isotropic sampling scheme, and is of the following form:

[0166]

[0167] in, represents the derivative in the u direction, the superscript (2) represents the vertex scheme of one of the ID schemes, and α is the weighting coefficient when using the vertex scheme;

[0168] In the x-direction, the relationship between the vertex solution and the original solution has rotation symmetry and is expressed as follows:

[0169]

[0170] To simplify the expression, a new shift operator is defined as follows:

[0171]

[0172] Finally, the differential operation of the new algorithm will consist of three parts:

[0173]

[0174] The spatial coordinates for the difference are (i, j, k). The three parts of the above formula from top to bottom represent the original coordinate point, the one-shift point, and the two-shift point respectively; the one-shift and two-shift refer to the number of shift operators. Overall, isotropic sampling is achieved by expanding the sampling points involved in the difference.

[0175] Applying the above transformation to HIE-FDTD, that is, performing ID processing on the spatial differences in the x and z directions, but not on the spatial differences in the y direction, the updated equation of the electromagnetic field is obtained as follows:

[0176]

[0177] This method still maintains the tridiagonal matrix form when solving implicitly, and can be solved iteratively in a simple manner.

[0178] Example 1 shows the normalized dispersion error curve and the three-dimensional dispersion plane. Figure 1 and Figure 2 .v p is the numerical phase velocity; in order to quantitatively observe the dispersion error of the method of applying the universal low numerical dispersion finite-difference time-domain method to the HIE case and the traditional method. The present invention uses the normalized phase velocity and the normalized phase velocity error for comparison, and the normalized phase velocity and the normalized phase velocity error are defined as v respectively. p / c and (1-v p / c); the grid size of the discrete cells is Δx = Δz = Δy = λ / 10 = 1 cm, λ is the wavelength, the time step and frequency are Δt = Δy / 2c and the source frequency f = 3 GHz, respectively; to demonstrate the dispersion characteristics, a comparison between the universal low numerical dispersion finite-difference time-domain method applied to the HIE case and the traditional HIE-FDTD method is given through MATLAB simulation under the same parameters.

[0179] Figure 1The normalized phase velocity error curves for conventional HIE-FDTD and the modified new algorithm are shown. Clearly, the normalized phase velocities of the new algorithm are more isotropic, and after the correction, the new algorithm can achieve nearly exact isotropic dispersion. This confirms the correctness of the proposed method and demonstrates that for the implicit HIE-FDTD method in the y-direction, excellent correction can be achieved simply by introducing ID operations in the x and z directions.

[0180] Figure 2 The superiority of the new algorithm's dispersion can be intuitively appreciated through the three-dimensional plot. The yellow surface shows the dispersion of the method using the universal low-dispersion finite-difference time-domain method applied to the HIE case, the green surface shows the dispersion plot of the conventional FDTD method, and the blue surface shows the dispersion plot of the traditional HIE-FDTD method. The yellow surface is a plane at all angles, and the value of each point on the plane is 1. The blue and green surfaces, however, fluctuate with the propagation angle, and their values ​​are all less than 1, confirming the superiority of the proposed method.

[0181] Example 2 simulates a 70 nm thick (relative dielectric constant ε r =2) grating structure, and the results calculated by rigorous coupled wave analysis (RCWA) are used as the benchmark. Figure 3 and Figure 4 .

[0182] A Gaussian pulse source was used to generate plane wave excitation. The grid size was set to Δx = Δz = Δy = 1 nm, and the time step was Δt = Δx / 2c = 1.67 × 10 -18 Mur absorbing boundary conditions were applied in the propagation direction of the excitation source, while periodic boundary conditions were set along the x and y directions. Furthermore, the recording plane was appropriately positioned to capture both the transmitted and reflected electric fields. Specifically, the wave source was located 308 nm to the left of the dielectric grating, the transmitted electric field recording plane was located 470 nm to the right of the dielectric grating, and the reflected electric field was located 80 nm to the left of the wave source. The reflected electric field recording plane Fourier transforms the recorded reflected and transmitted electric fields, as well as the incident electric field. The ratios of these values ​​allow the electric field reflection and transmission coefficients to be calculated under the broadband excitation source.

[0183] The proposed method is highly consistent with the results obtained from rigorous coupled wave analysis. Further observation revealed that the proposed method and the ID-FDTD method exhibited almost identical errors. However, under fine grid conditions, the proposed method required less CPU time and memory than the ID-FDTD method. For a simulation involving 4000 time steps, the proposed method took 118.19 seconds, a significant improvement compared to the 151.17 seconds required by the ID-FDTD method. This also confirms the high efficiency of the proposed method.

[0184] Example 3 simulates the scattering characteristics of a lossy dielectric sphere to verify the effectiveness of the proposed algorithm in dealing with complex media with energy dissipation and phase sensitivity problems. Figure 5 and Figure 6 .

[0185] Parameters of simulation settings: cell size Δx = 5 Δy = Δz = 1.67 × 10 -7 The number of grid cells in the x, y, and z directions is 96. The CFL condition is Δt = Δx / 3c. The relative permittivity of the sphere is 6-0.1j, and its radius is 12 grid cells, equivalent to 16 microns. The observation line is set 36 grid cells from the center of the sphere along the propagation direction of the incident wave.

[0186] Since the ID-FDTD method agrees very well with the exact Mie series, we consider the results obtained by the ID-FDTD method as the reference solution. The scattered field values ​​along the observation line are calculated as follows. Figure 5 、 Figure 6 As shown in the figure, along the observation line, the proposed method shows a high degree of agreement with the reference solution in the real part, comparable to the traditional hybrid implicit explicit finite-difference time-domain (HIE-FDTD) method. In the imaginary part, the proposed method demonstrates a clear advantage, producing more accurate results. In electromagnetic field problems, the phase of the electric field is determined by the ratio of the real to imaginary parts. Therefore, the proposed method is particularly advantageous in dealing with phase-sensitive problems.

[0187] This document uses specific examples to illustrate the principles and implementation methods of the present invention. The above examples are only intended to help understand the method and core concept of the present invention. At the same time, those skilled in the art will find that the specific implementation methods and application scopes may vary based on the concept of the present invention. In summary, the contents of this specification should not be construed as limiting the present invention.

Claims

1. A universal low numerical dispersion finite-difference time-domain method, characterized by The following steps are involved: Step 1. Based on the physical model structure of the electromagnetic simulation project, analyze whether it has a fine structure direction, and then select the implicit or explicit subclass method of finite-difference time-domain (FDTD) to be applied; Step 2. Based on the selected finite-difference time-domain (FDTD) method, write the dispersion equation derived from the plane wave solution, apply the isotropic sampling scheme to the selected direction, and obtain the optimal weighting factor and correction factor; And verify the stability CFL conditions; Step 3. Based on the obtained optimal weighting factor, reconstruct the iterative formula of the electric field and magnetic field as the time step advances; Step 4. Perform calculations and solve the problem based on the newly obtained iterative formula. The calculations involving implicit solutions are calculated using tridiagonal matrices. Step 5. Based on the above steps, the electric field parameters and magnetic field parameters in the entire area under the specific simulation example are obtained, so as to obtain the desired physical parameters through calculation.

2. The universal low numerical dispersion finite-difference time-domain method according to claim 1, characterized in that: In step 1, the selection of the implicit, explicit, or hybrid explicit-implicit subclass method of the finite-difference time-domain FDTD to be applied is specifically as follows: for ordinary objects, the ordinary FDTD method is selected; for methods with fine structures in one direction, the hybrid explicit-implicit finite-difference time-domain method HIE-FDTD method is adopted; for methods with fine structures in two directions, the weak conditionally stable finite-difference time-domain method WCS-FDTD is selected.

3. The universal low numerical dispersion finite-difference time-domain method according to claim 2, characterized in that: In step 1, the fine structure direction is determined by the meshing method based on the fine structure direction of the specific physical model. Assuming that the object has a fine structure in the y direction, the y direction is called the implicit direction, and the y-direction mesh size Δy is irrelevant to the stability condition. The x-direction mesh size Δx and the z-direction mesh size Δz do not participate in the stability condition.

4. The universal low numerical dispersion finite-difference time-domain method according to claim 1, characterized in that: The dispersion equation described in step 2 is the dispersion equation for one-step HIE: in: k x 、k y 、k z are the numerical wave number components in the three directions in the rectangular coordinate system, where k z = kcosθ. k is the numerical wave number during simulation; where θ, are the elevation angle and azimuth angle in the spherical coordinate system, c is the speed of light in vacuum, Δt is the time step of calculation, ω is the angular frequency of the wave; J x 、J y 、J z It is to simplify the algebraic expressions in the x, y, and z directions; If the isotropic sampling ID sampling scheme is applied only to the x and z directions, then J in the dispersion equation x and J y Make the following corrections: The superscript ID indicates that an isotropic sampling scheme is applied in this direction, and α is the weighting factor.

5. The universal low numerical dispersion finite-difference time-domain method according to claim 1, characterized in that: Step 3: Based on the obtained optimal weighting factor, the iterative formula for reconstructing the electric and magnetic fields as time steps is specifically determined by applying the ID sampling scheme in the x and z directions. The maximum normalized phase velocity of the dispersion equation of the one-step HIE-FDTD method is R max , the minimum value is R min , let w=R max -R min， By traversing the parameters in MATLAB, we find the value of α when w is the minimum, which is called the optimal weighting factor. The dispersion equation of the universal low numerical dispersion finite-difference time-domain method is applied to the HIE case: The superscript ID indicates that an isotropic sampling scheme is applied in this direction.

6. The universal low numerical dispersion finite-difference time-domain method according to claim 1 or 4, characterized in that: The optimal weighting factor and correction factor are obtained, and the As a reference, the dispersion equation is used, where θ, They are the elevation angle and azimuth angle in the spherical coordinate system, and the precise wave number k is obtained. exact The correction factor sc is determined:

7. The universal low numerical dispersion finite-difference time-domain method according to claim 6, characterized in that: The correction factor sc acts on the dielectric constant and permeability ε′=sc×ε, μ′=sc×μ to obtain a dispersion curve with isotropy and almost close to the exact phase velocity; ε′ is the corrected dielectric constant, ε is the dielectric constant before correction, μ′ is the corrected permeability, and μ is the permeability before correction.

8. The universal low numerical dispersion finite-difference time-domain method according to claim 1, characterized in that: Step 4: Calculate and solve the new iterative formula based on the main body, and analyze the iterative formula of the main body, that is, the electric field and magnetic field changing with time. The calculation of the three-dimensional example includes the electric field value E in the x, y, and z directions. x 、E y 、E z and magnetic field value H x 、H y 、H z Applying the isotropic sampling ID scheme to the differences in the x and z directions yields: The superscript n represents a specific time point in the simulation process, ε is the dielectric constant before correction, and μ is the magnetic permeability before correction.

9. The universal low numerical dispersion finite-difference time-domain method according to claim 1, characterized in that: Step 5: Based on the above steps, obtain the electric field parameters and magnetic field parameters in the entire area under the specific simulation example. After selecting the observation point or observation plane, obtain the desired physical parameters including but not limited to the transmission coefficient by calculation. Reflection coefficient Normalized phase velocity v p =ω / k.

Images