Electromagnetic thermal multi-physical field simulation method and system based on FDTD / FVTD algorithm

Abstract

The invention provides an electromagnetic thermal multi-physical field simulation method based on an FDTD / FVTD algorithm, and belongs to the technical field of electromagnetic simulation. Performing time domain finite difference discretization on the established spatial discrete model by adopting a Yee cellular structure, and calculating electromagnetic field distribution according to a Maxwell rotation equation; constructing an adaptive steady-state judgment criterion based on the instantaneous total energy change rate, and updating electromagnetic field distribution of the spatial discrete model by adopting an FDTD explicit method; calculating the electromagnetic loss power density of the spatial discrete model, and transmitting the electromagnetic loss power density into a thermal field for calculation; based on an FVTD method, thermal field distribution of the discrete model is calculated through an explicit or implicit format, and thermal field distribution information is obtained. The FDTD method is based on explicit time domain iteration, transient behaviors such as electromagnetic wave propagation and scattering can be visually captured, and the FDTD method has broadband response capability and complex medium adaptability and is especially suitable for the high-frequency electromagnetic problem. The FVTD method strictly ensures energy conservation, directly treats boundary heat flux and is suitable for non-uniform materials and complex heat transfer boundary conditions.

Description

Technical Field

[0001] The present invention belongs to the technical field of electromagnetic simulation, and in particular relates to an electromagnetic thermal multi-physics field simulation method and system based on FDTD / FVTD algorithm. Background Art

[0002] With the rapid development of high-frequency electronic devices, microwave communication systems, and high-power electromagnetic equipment, electromagnetic-thermal coupling has become a key factor restricting device performance and reliability. Traditional single-physics simulation methods (such as analyzing only the electromagnetic field or thermal field) have difficulty accurately predicting complex behaviors under multi-physics coupling, such as heat accumulation, material loss, and thermal stress distribution under transient electromagnetic excitation. Furthermore, multi-physics simulation of large-scale complex structures also faces the dual challenges of computational efficiency and accuracy. To address these issues, there is an urgent need to research and develop efficient and high-precision electromagnetic-thermal multi-physics simulation technologies. Summary of the Invention

[0003] The purpose of the present invention is to overcome the above-mentioned shortcomings and provide an electromagnetic-thermal multi-physics field simulation method and system based on the FDTD / FVTD algorithm. This solution uses the Finite Difference Time Domain (FDTD) method to solve the electromagnetic field, and based on the same set of structured grids, uses the Finite Volume Time Domain (FVTD) method to solve the thermal field, realizing electromagnetic-thermal one-way coupled simulation. At the same time, this solution introduces an adaptive steady-state judgment method in the electromagnetic field calculation and uses the stable biconjugate gradient iteration method in the thermal field calculation to achieve efficient solution of large-scale sparse linear equations, thereby improving the efficiency and accuracy of the calculation.

[0004] In order to achieve the above object, the present invention adopts the following technical solutions:

[0005] In a first aspect, the present invention provides an electromagnetic thermal multi-physics field simulation method based on the FDTD / FVTD algorithm, comprising the following steps:

[0006] A spatial discrete model of isotropic media is established in the passive region. The established spatial discrete model is discretized in the time domain using the Yee cell structure. The time-marching formulas of the electric and magnetic fields are constructed based on the Maxwell curl equations to calculate the electromagnetic field distribution.

[0007] An adaptive steady-state judgment criterion is constructed based on the instantaneous total energy change rate, and the FDTD explicit method is used to update the electromagnetic field distribution of the spatial discrete model. The calculation is terminated when the electromagnetic field reaches a steady state.

[0008] When the electromagnetic field reaches a steady state, the electromagnetic loss power density of the spatial discrete model is calculated and passed into the thermal field calculation;

[0009] Based on the FVTD method, the thermal field distribution of the discrete model is calculated by explicit or implicit format to obtain the thermal field distribution information.

[0010] The steps of establishing a spatial discrete model of an isotropic medium in a passive region, performing time-domain finite-difference discretization on the established spatial discrete model using a Yee cell structure, constructing a time-marching formula for the electric field and the magnetic field according to the Maxwell curl equation, and calculating the electromagnetic field distribution are as follows:

[0011] Construct a three-dimensional Yee cell structure and perform time-domain finite-difference discretization on the three-dimensional Yee cell structure, with the electric field component arranged on the cell edges and the magnetic field component arranged at the center of the surface;

[0012] The Maxwell curl equations are discretized by central difference to obtain the propulsion formulas for the electric and magnetic fields. The specific formulas are as follows:

[0013]

[0014] in, Indicates that at time n, The electric field E at the grid point x Component value; Indicates time, The magnetic field H at the grid point z Component value; CA(m) is The coefficients of the equation corresponding to the electromagnetic parameters at the grid points, The coefficients of the equations corresponding to the electromagnetic parameters at the grid points; σ(m) is the grid point The conductivity at the grid point, ε(m) is The dielectric constant at .

[0015] The Courant stability condition is used to determine the time step when the Maxwell curl equation is discretized by central difference. The specific formula is as follows:

[0016]

[0017] Where Δx, Δy, and Δz represent the grid space steps in the x, y, and z coordinate directions, respectively; Δt represents the time step; and c represents the propagation speed of electromagnetic waves in a specific medium.

[0018] In the step of constructing an adaptive steady-state judgment criterion based on the instantaneous total energy change rate, updating the electromagnetic field distribution of the spatial discrete model using the FDTD explicit method, and terminating the calculation when the electromagnetic field reaches a steady state, the adaptive steady-state judgment criterion is as follows:

[0019] For scattering problems, the electromagnetic field is calculated by the transient energy within the closed surface S;

[0020] For radiation problems, the electromagnetic field is calculated by adding the power term of the excitation source;

[0021] Set the stable convergence parameter and terminate the calculation when the relative wave pair energy of adjacent periods is less than the stable convergence parameter.

[0022] For scattering problems, the formula for calculating the electromagnetic field through the transient energy in the closed surface S is expressed as follows:

[0023]

[0024] Where n is the index of the time step, i is the index variable, M represents the total number of discrete surface elements, subscript m represents the index of the discrete surface element, E m (iΔt) is the electric field vector of the mth discrete element at the moment (iΔt), H m (iΔt) is the magnetic field vector of the mth discrete element at the moment (iΔt), Δs m represents the area vector associated with the mth discrete surface element.

[0025] For the radiation problem, the formula for calculating the electromagnetic field by adding the excitation source power term is expressed as follows:

[0026]

[0027] Where N and N′ represent the number of grid edges where the excitation electric field and excitation magnetic field are located, respectively, and n and n′ are the corresponding indices.

[0028] When the electromagnetic field reaches a steady state, the electromagnetic loss power density of the spatial discrete model is calculated and passed into the thermal field calculation step. The total electromagnetic loss power density of the electromagnetic field in the medium can be expressed as:

[0029]

[0030] in, is the conductive loss, also known as Joule heat, σ is the conductivity of the medium, and |E| is the amplitude of the electric field strength; is the polarization loss, ω represents the angular frequency of the electromagnetic wave, ε0 is the vacuum dielectric constant, and ε″ is the dielectric loss factor of the medium; is the magnetization loss, μ0 is the vacuum permeability, and μ″ is the magnetic loss factor of the medium.

[0031] In the step of obtaining thermal field distribution information by calculating the thermal field distribution of the discrete model in an explicit or implicit format based on the FVTD method, the law of conservation of energy is applied to the cells. The specific formula is as follows:

[0032]

[0033] in are the internal energy of the cell at time t1 and t0 respectively, and are the heat flowing from the six surfaces of the cell from t0 to t1, is the heat generated by the heat source in the cell from t0 to t1.

[0034] The stable biconjugate gradient method is used to calculate the thermal field distribution of the implicit format discrete model.

[0035] In a second aspect, the present invention provides an electromagnetic thermal multi-physics field simulation system based on the FDTD / FVTD algorithm, comprising:

[0036] The first module is used to establish a spatial discrete model of isotropic media in the passive region. The Yee cell structure is used to perform time-domain finite-difference discretization on the established spatial discrete model. The time-marching formula of the electric and magnetic fields is constructed based on the Maxwell curl equation to calculate the electromagnetic field distribution.

[0037] The second module is used to construct an adaptive steady-state judgment criterion based on the instantaneous total energy change rate, and uses the FDTD explicit method to update the electromagnetic field distribution of the spatial discrete model. The calculation is terminated when the electromagnetic field reaches a steady state.

[0038] The third module is used to calculate the electromagnetic loss power density of the spatial discrete model after the electromagnetic field reaches a steady state, and pass it into the thermal field calculation;

[0039] The fourth module is used to calculate the thermal field distribution of the discrete model based on the FVTD method through explicit or implicit format to obtain thermal field distribution information.

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

[0041] The present invention provides an electromagnetic thermal multi-physics field simulation method based on the FDTD / FVTD algorithm, comprising the following steps: establishing a spatial discrete model of an isotropic medium in a passive region, discretizing the established spatial discrete model using a Yee cell structure in the time domain (FDTD), constructing time-marching formulas for the electric and magnetic fields based on the Maxwell curl equations, and calculating the electromagnetic field distribution; constructing an adaptive steady-state judgment criterion based on the instantaneous total energy change rate, updating the electromagnetic field distribution of the spatial discrete model using the FDTD explicit method, and terminating the calculation when the electromagnetic field reaches a steady state; after the electromagnetic field reaches a steady state, calculating the electromagnetic loss power density of the spatial discrete model and passing it into the thermal field calculation; and calculating the thermal field distribution of the discrete model using the FVTD method using an explicit or implicit format to obtain thermal field distribution information. The FDTD method, based on explicit time-domain iteration, can intuitively capture transient behaviors such as electromagnetic wave propagation and scattering, has wideband response capabilities and adaptability to complex media, making it particularly suitable for high-frequency electromagnetic problems; the FVTD method strictly ensures energy conservation and directly processes boundary heat flux, making it suitable for inhomogeneous materials and complex heat transfer boundary conditions.

[0042] Furthermore, an adaptive steady-state judgment method was introduced in the electromagnetic field calculation. It periodically monitors the relative fluctuations of the instantaneous total energy in the calculation domain, dynamically judges whether the electromagnetic field has reached a steady state, avoids invalid iterations or premature termination, and significantly saves calculation time while ensuring accuracy.

[0043] Furthermore, a stable biconjugate gradient iteration method is introduced in the thermal field calculation. It improves the efficiency and stability of solving asymmetric matrices for large-scale sparse linear equations generated by implicit FVTD discretization, and is suitable for efficient solution of complex three-dimensional thermal fields.

[0044] Furthermore, in the present invention, electromagnetic field and thermal field calculations share the same set of structured grids, which improves coupling efficiency and reduces memory requirements. BRIEF DESCRIPTION OF THE DRAWINGS

[0045] Figure 1 Flow chart of the method of the present invention;

[0046] Figure 2 Schematic diagram of the Yee cell structure in the present invention;

[0047] Figure 3 Schematic diagram of the closed boundary surface S in the scattering problem of the present invention;

[0048] Figure 4 Schematic diagram of the closed boundary surface S in the radiation problem of the present invention;

[0049] Figure 5 Schematic diagram of cell energy change in the present invention;

[0050] Figure 6Schematic diagram of heat exchange on a solid boundary in the present invention;

[0051] Figure 7 Schematic diagram of heat exchange on the third type of boundary in the present invention;

[0052] Figure 8 Schematic diagram of the simulation model structure of Example 2 of the present invention;

[0053] Figure 9 Schematic diagram of the change of the stability coefficient over time in Example 2 of the present invention;

[0054] Figure 10 is the electromagnetic loss power P of the dielectric substrate in Example 2 of the present invention d Schematic diagram of calculation results;

[0055] Figure 11 Schematic diagram of thermal field simulation results of the dielectric substrate surface in Example 2 of the present invention;

[0056] Figure 12 Schematic diagram of the simulation results of the maximum temperature point of the quality substrate in Example 2 of the present invention. DETAILED DESCRIPTION

[0057] In order to further understand the content of the present invention, the present invention is described in detail below in conjunction with the accompanying drawings and specific embodiments. It should be understood that the embodiments are only for explaining the present invention and are not intended to limit it.

[0058] Example 1

[0059] An electromagnetic thermal multi-physics field simulation method based on FDTD / FVTD algorithm includes the following steps:

[0060] S1: Establish a spatial discrete model of isotropic media in the passive region, perform time-domain finite-difference discretization on the established spatial discrete model using the Yee cell structure, construct the time-marching formula of the electric and magnetic fields based on the Maxwell curl equations, and calculate the electromagnetic field distribution;

[0061] S2: Based on the instantaneous total energy change rate, an adaptive steady-state judgment criterion is constructed, and the electromagnetic field distribution of the spatial discrete model is updated using the FDTD explicit method. The calculation is terminated when the electromagnetic field reaches a steady state.

[0062] S3: When the electromagnetic field reaches a steady state, the electromagnetic loss power density of the spatial discrete model is calculated and passed into the thermal field calculation;

[0063] S4: Based on the FVTD method, the thermal field distribution of the discrete model is calculated using an explicit or implicit format to obtain the thermal field distribution information.

[0064] Specifically, in S1, a spatial discrete model of an isotropic medium is established in the passive region. The established spatial discrete model is discretized in the time domain using the Yee cell structure. The time-marching formulas of the electric and magnetic fields are constructed based on the Maxwell curl equations to calculate the electromagnetic field distribution. The specific method is as follows:

[0065] The finite-difference time-domain method (FDTD) is a numerical method widely used in computational electromagnetics. It is primarily used to solve the time-varying Maxwell equations and simulate the propagation, scattering, and radiation behavior of electromagnetic waves in complex media or structures. Due to its intuitive time-domain characteristics, wide bandwidth, and adaptability to complex media, the FDTD method has become an important tool in electromagnetic field simulation, particularly suitable for transient problems and high-frequency applications.

[0066] Considering a passive region in space, assuming that the medium parameters do not change with time and are isotropic, the Maxwell curl equation is:

[0067]

[0068] Where, E is the electric field intensity, D is the electric flux density, H is the magnetic field intensity, B is the magnetic flux density, J is the current density, J m is the magnetic flux density. The medium constitutive relation is:

[0069] D=εE (7)

[0070] B=μH (8)

[0071] J=σE (9)

[0072] J m =σ m H (10)

[0073] Among them, ε is the dielectric constant of the medium, μ is the magnetic permeability coefficient, σ is the conductivity, σ m is the magnetic permeability.

[0074] In the rectangular coordinate system, equations (1) and (2) are expressed as the following scalar equations:

[0075]

[0076] The above equation is discretized by FDTD difference, and the three-dimensional space to be solved is meshed. The mesh space steps in the x, y and z coordinate directions are Δx, Δy and Δz respectively. The time is discretized with Δt as the time step. The discrete values ​​in the time and space domains can be expressed as:

[0077] f(x,y,z,t)=f(iΔx,jΔy,kΔz,nΔt)=f n (i,j,k) (13)

[0078] Here, i, j, k represent the electric field components taken at integer grid points, i.e., i, j, k = 0, 1, 2…, N.

[0079] For electromagnetic field calculation problems, the field components at each point in space are distributed according to Yee cells. Figure 2 As shown, all electric field components are located at the cell edges, while all magnetic field components are distributed at the center points of the six faces. Each magnetic field component is surrounded by four adjacent electric field components, and each electric field component is also surrounded by four adjacent magnetic field components. The electric and magnetic fields are sampled alternately in time, with the sampling intervals differing by half a time step. This allows the Maxwell curl equations to be discretized into explicit difference equations, which can be solved iteratively over time.

[0080] The time-marching formulas for the electric and magnetic fields in FDTD are obtained by discretizing the first terms of (7) and (8) using the above method:

[0081]

[0082] in, Indicates that at time n, The electric field E at the grid point x Component value; Indicates time, The magnetic field H at the grid point z Component value; CA(m) is The coefficients of the equation corresponding to the electromagnetic parameters at the grid point, CB(m) are The coefficients of the equations corresponding to the electromagnetic parameters at the grid points; σ(m) is the grid point The conductivity at the grid point, ε(m) is The dielectric constant at .

[0083]

[0084]

[0085] in, express time, Magnetic field H at the grid point x Component value; Indicates that at time n, The electric field E at the grid point z Component value; CP(m) is The coefficients of the equation corresponding to the electromagnetic parameters at the grid points, CQ(m) are The coefficient of the equation corresponding to the electromagnetic parameters at the grid point); σ m (m) is The magnetic permeability at the grid point, μ(m) is Magnetic permeability at a grid point;

[0086] Similarly, the time advancement formulas for the components of the electric field and magnetic field in other directions can be obtained.

[0087] Furthermore, numerical convergence and stability are prerequisites for the application of the FDTD algorithm. In order to prevent the numerical overflow of the field value as the number of iterations increases, the time step and the unit grid size in the FDTD algorithm must satisfy a certain relationship. This relationship is also called the Courant stability condition and is given by the following formula:

[0088]

[0089] Where c represents the propagation speed of electromagnetic waves in a specific medium.

[0090] Specifically, in S2, an adaptive steady-state judgment criterion is constructed based on the instantaneous total energy change rate, the FDTD explicit method is used to update the electromagnetic field distribution of the spatial discrete model, and the iterative calculation of the electromagnetic field distribution is dynamically terminated; the specific method is as follows:

[0091] Since the FDTD method is a time-domain numerical algorithm, the time-domain numerical algorithm realizes the process from transient state to steady state of the system through time-domain iteration. In order to avoid the algorithm falling into an infinite iterative loop, and also to avoid the algorithm terminating the calculation prematurely before the electromagnetic field reaches a steady state, thereby affecting the accuracy of the results, it is necessary to set a suitable iteration termination condition for the algorithm.

[0092] To address the above problem, the time domain iteration process is terminated by periodically monitoring the changes in the instantaneous total energy in the calculation area when the instantaneous total energy meets the preset convergence threshold.

[0093] like Figure 3 As shown in the figure, for the scattering problem, according to Poynting's theorem, a closed surface S is taken for the scattering target within the FDTD total field boundary, then the energy in the volume V surrounded by the surface S satisfies the following formula:

[0094]

[0095] in is the power entering the closed boundary S, represents the power generated by the external and internal sources in the closed area V, represents the external conduction current density, Represents the external magnetic current density. is the power dissipated due to material losses, represents the ohmic loss current density, represents the magnetic loss current density; It represents the power stored in the closed area V; D is the electric displacement vector, and B is the magnetic induction intensity.

[0096] Then the transient energy Egy(t) stored in the finite volume V at time t can be obtained by integrating the left side of the equation over time, that is,

[0097]

[0098] Since there is no external source in the closed volume V, Substituting this condition into the above formula and discretizing the time variable, we get:

[0099]

[0100] Where n is the index of the time step, i is the index variable, M represents the total number of discrete surface elements, subscript m represents the index of the discrete surface element, E m (iΔt) is the electric field vector of the mth discrete element at the moment (iΔt), H m (iΔt) is the magnetic field vector of the mth discrete element at the moment (iΔt), Δs m represents the area vector associated with the mth discrete surface element.

[0101] like Figure 4 As shown, for radiation problems, the excitation source is usually placed in the computational domain, so the external source power term in equation (18) is no longer zero. The power generated by the discrete source can be obtained by multiplying the port voltage and the port current:

[0102]

[0103] Substituting the above formula into formula (18) and numerically discretizing the time variable, we can obtain:

[0104]

[0105] Where N and N′ represent the number of grid edges where the excitation electric field and excitation magnetic field are located, respectively, and n and n′ are the corresponding indices.

[0106] To quantitatively analyze the transition from transient to steady state in a time-domain iterative algorithm, this can be achieved by setting steady-state convergence parameters. For time-harmonic field signal excitation, since it is an energy-infinite signal, the total energy within the computational domain typically does not decay to infinitesimal values. In such cases, the relative energy fluctuations between two adjacent cycles can be used as a basis for determining whether the energy-infinite signal has reached steady state. Therefore, the following stability coefficient ξ(t) is used:

[0107]

[0108] The iterative calculation is terminated when the relative wave pair energy of adjacent periods is less than the stable convergence parameter.

[0109] Specifically, in S3, when the electromagnetic field reaches a steady state, the electromagnetic loss power density of the spatial discrete model is calculated and passed to the thermal field calculation. The specific method is as follows:

[0110] The key to electromagnetic-thermal coupling calculations is determining how much of the lost electromagnetic energy in a medium is converted into heat. When a lossy medium is in an electromagnetic field, the electromagnetic energy lost within the medium can be quantitatively calculated using the electromagnetic loss power. Depending on how the electromagnetic energy is lost, the corresponding electromagnetic loss power calculation formula also varies. When the applied electric field is time-harmonic and the conductivity of the medium is not zero, the total electromagnetic loss power density of the electromagnetic field in the medium can be expressed as:

[0111]

[0112] in, is the conductive loss, also known as Joule heat, σ is the conductivity of the medium, and |E| is the amplitude of the electric field strength; is the polarization loss, ω represents the angular frequency of the electromagnetic wave, ε0 is the vacuum dielectric constant, and ε″ is the dielectric loss factor of the medium; is the magnetization loss, μ0 is the vacuum permeability, and μ″ is the magnetic loss factor of the medium.

[0113] Without considering the magnetization loss, the above formula becomes:

[0114]

[0115] Where E is the root mean square of E. For the time-harmonic field, the components of E in the x, y, and z directions can be calculated using the following formula:

[0116]

[0117] Where, Δt e represents the sampling time interval, is a complex exponential term used to project the time domain signal into the frequency domain.

[0118] Since the electrothermal calculations share the same grid, and the electromagnetic loss P d (i, j, k) is located at the center of the Yee cell, so it is necessary to introduce a linear interpolation method to obtain the center position of the Yee cell. value. The calculation formulas for each component are:

[0119]

[0120] Specifically, in S4, based on the FVTD method, the thermal field distribution of the discrete model is calculated by explicit or implicit format to obtain the thermal field distribution information. The specific method is as follows:

[0121] The finite volume time domain (FVTD) method is an efficient numerical method for solving heat transfer problems. Its core idea is to discretize the calculation area into non-overlapping control volumes (cells). Discretization is achieved by integrating the energy conservation equation within each control volume, thereby directly dealing with the flux exchange of heat flow at the cell boundary and strictly ensuring local and global energy conservation.

[0122] like Figure 5 The cell shown is fixed in rectangular coordinate space. The law of conservation of energy is applied to the cell. That is, the increase in energy in the cell from time t0 to time t1 is equal to the sum of the heat flowing from the cell surface and the heat generated inside the cell during this period. The formula is expressed as:

[0123]

[0124] in are the internal energy of the cell at time t1 and t0 respectively, and are the heat flowing from the six surfaces of the cell from t0 to t1, is the heat generated by the heat source in the cell from t0 to t1.

[0125]

[0126] Where C is the specific heat capacity of the cellular material, ρ is the density of the cellular material, are the temperature values ​​of the center of the cell at time t1 and t0 respectively.

[0127] According to Fourier's law of heat conduction, the amount of heat conducted through a given unit cross section per unit time is proportional to the temperature change rate and cross-sectional area in the direction perpendicular to the cross section, and the direction of heat transfer is opposite to the direction of temperature increase. The formula is expressed as:

[0128]

[0129] Where k is the thermal conductivity of the cellular material.

[0130]

[0131] Among them, s p is the power density of the heat source in the cell.

[0132] Substituting equations (32)-(35) into equation (31), we obtain:

[0133]

[0134] Where ρ is the material density; C is the material specific heat capacity; are the temperature values ​​of the cell center at time t1 and t0 respectively; Δt is the time step t1-t0, q e ,q w ,q s ,q n ,q b and q f is the heat flux density through the six surfaces of the cell (the amount of heat conducted through a given unit cross-section per unit time); s p is the power density of the heat source in the cell.

[0135] like Figure 6 As shown, taking the grid in the x direction as an example, for solid material grids, the heat flux density at the interface between adjacent grids can be expressed as:

[0136]

[0137]

[0138] Similarly, the results in the y and z directions can be obtained.

[0139] For the third type of boundary conditions, the heat flux density is determined by the surface heat transfer coefficient h between the object and the surrounding fluid (usually air) on the boundary surface and the temperature T of the surrounding fluid environment. air Given. Figure 7 For example, there is a heat flux on the boundary:

[0140]

[0141] The heat flux density in other directions can be obtained in the same way.

[0142] When the heat flux q is taken as the temperature variable and When , it is the explicit FVTD heat transfer method; the explicit FVTD heat transfer method needs to meet the time stability condition, the specific expression is

[0143]

[0144] When the heat flux q is taken as the temperature variable and When , it is the implicit FVTD heat transfer method; the implicit FVTD heat transfer method is unconditionally stable, but it needs to solve a large-scale linear algebraic equation system in terms of calculation, so it is necessary to use a suitable method to solve it.

[0145] Methods for solving algebraic equations fall into two main categories: direct and iterative methods. The BiCGSTAB (Stabilized Biconjugate Gradient) method is an iterative algorithm used in numerical heat transfer to solve large, sparse linear systems. It is particularly well-suited for nonsymmetric matrices. Proposed by H.A. van der Vorst in 1992, the BiCGSTAB method aims to improve upon the traditional BiCG method, enhancing its convergence and stability. By introducing polynomial acceleration techniques, BiCGSTAB reduces oscillations during iterations, thereby accelerating convergence.

[0146] The implicit FVTD heat transfer discrete equation is expressed in matrix form:

[0147] AT=b (45)

[0148] Where A is the coefficient matrix, which is generated by the discretized heat transfer equation and is affected by thermal parameters and geometric parameters; T is the temperature vector, which represents the temperature field to be determined; b is a known vector containing known conditions, such as heat sources and boundary conditions.

[0149] The solution steps of the stable biconjugate gradient method are as follows:

[0150] 1. Initialization:

[0151] 1) Select the initial guess T0;

[0152] 2) Calculate the residual r0 = b-AT0;

[0153] 3) Settings

[0154] 2. Iterative process (for k = 0, 1, 2, ...):

[0155] calculate If k=0, set p k =r k Otherwise, calculate β k =(ρ k / ρ k-1 )(α k-1 / ω k-1 ); Update p k =r k +β k (p k-1 -ω k-1 v k-1 );

[0156] Calculate v k =Ap k ;

[0157] calculate

[0158] Calculate s=r k -α k v k

[0159] Calculate t = As;

[0160] Calculate ω k =(t T s) / (t T t);

[0161] Update solution T k+1 =T k +α k p k +ω k s

[0162] Update residual r k+1 =s-ω k t;

[0163] Check convergence conditions and stop if they are met.

[0164] Where T0 is the initial temperature; T k is the temperature solution of step k; r0 is the initial residual vector; is the initial residual projection vector; r k is the residual vector after the k-th iteration, which represents the error between the current temperature solution and the true temperature solution; ρ k is the residual vector r k and the initial residual projection vector The inner product of p reflects the projection change of the residual; k is the search direction vector of the kth step; β k is the BiCG direction update coefficient of the kth step; v k is the coefficient matrix A and the search direction p k The product of is used to calculate α k ; α k is the step size factor, controlling the search direction p k The movement amount is adjusted to adjust the update amplitude of the solution to reduce the residual; s is the intermediate residual vector, indicating the movement along the direction p k The residual error after the move; t is the product of the coefficient matrix A and the intermediate residual s, which is used to calculate ω k ;ω k To stabilize the step size, the algorithm stability is enhanced by minimizing the norm of the intermediate residual s.

[0165] Example 2

[0166] Based on the method of the present invention, the frequency selective surface unit is composed of two parts: a cross-shaped gap metal layer and an FR-4 dielectric substrate. The specific method is as follows:

[0167] like Figure 8 As shown, the electromagnetic parameters of the metal layer are: ε r =1,μ r =1,σ=1e10,σ m =0; the electromagnetic parameters of the dielectric substrate are: ε r =4.5, μ r =1,σ=0.004,σ m = 0. The thermal parameters of the metal layer are: thermal conductivity: 0.3W / K / m, specific heat capacity: 1369J / K / kg, density: 1900kg / m 3 ; Thermal parameters of dielectric substrate are: thermal conductivity: 400W / K / m, specific heat capacity: 385J / K / kg, density: 8960kg / m 3 The metal layer is 1 mm thick, and the dielectric substrate is 3 mm thick. A plane wave is incident perpendicularly on the FSS element in the +z direction, with the electric field polarized in the x direction. Assume the ambient temperature is 20°C and the convection heat transfer coefficient h = 5. The plane wave excitation source is a sinusoidal wave. The excitation source frequency is 2.185 GHz, and the amplitude is 4000 V / m. The model is meshed with a uniform structure of 0.5 mm × 0.5 mm × 0.5 mm.

[0168] The FDTD method is used to simulate the electromagnetic field of the model, and the adaptive steady-state judgment method is used to monitor the changes in electromagnetic field energy during the simulation process. Figure 9 As shown, at 2.7 ns, the stability coefficient ξ(t) is lower than the set convergence threshold of -20 dB, and it is judged that the electromagnetic field has reached the steady-state condition.

[0169] After the electromagnetic field reaches a steady state, the electromagnetic loss power density P in the model is calculated from the electric field root mean square E d , the calculation results are as follows Figure 10 shown.

[0170] The FVTD method is used to simulate the thermal field of the model and the heat transfer time is calculated to be 1200s. The thermal field simulation results on the dielectric plate are as follows: Figure 11 shown.

[0171] Further, such as Figure 12The figure shows the temperature variation over time at the maximum point of the thermal field on the dielectric substrate. Compared with the commercial multi-physics simulation software COMSOL, it can be seen that after 1200 seconds of calculation, the calculation result of this method is 206.28℃, and the calculation result of COMSOL software is 205.13℃. The relative error between the two results is only 0.56%, which shows that the accuracy of this method is consistent with that of existing commercial multi-physics calculation software. For large-scale electromagnetic thermal multi-physics coupling problems, the method of the present invention can better utilize parallel computing resources and improve computing efficiency. At the same time, its flexible algorithm and model definition method enable it to more accurately describe the complex coupling relationship between electromagnetics and heat. COMSOL may be limited by software performance and algorithms when dealing with such large-scale coupling problems.

[0172] Example 3:

[0173] An electromagnetic thermal multi-physics field simulation system based on FDTD / FVTD algorithm, including

[0174] The first module is used to establish a spatial discrete model of isotropic media in the passive region. The Yee cell structure is used to perform time-domain finite-difference discretization on the established spatial discrete model. The time-marching formula of the electric and magnetic fields is constructed based on the Maxwell curl equation to calculate the electromagnetic field distribution.

[0175] The second module is used to construct an adaptive steady-state judgment criterion based on the instantaneous total energy change rate, and uses the FDTD explicit method to update the electromagnetic field distribution of the spatial discrete model. The calculation is terminated when the electromagnetic field reaches a steady state.

[0176] The third module is used to calculate the electromagnetic loss power density of the spatial discrete model after the electromagnetic field reaches a steady state, and pass it into the thermal field calculation;

[0177] The fourth module is used to calculate the thermal field distribution of the discrete model based on the FVTD method through explicit or implicit format to obtain thermal field distribution information.

[0178] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit it. Although the present invention has been described in detail with reference to the above embodiments, ordinary technicians in the field should understand that the specific implementation methods of the present invention can still be modified or replaced by equivalents. Any modification or equivalent replacement that does not depart from the spirit and scope of the present invention should be covered by the scope of protection of the claims of the present invention.

Claims

1. An electromagnetic thermal multi-physics field simulation method based on FDTD / FVTD algorithm, characterized in that: The steps include: A spatial discrete model of isotropic media is established in the passive region. The established spatial discrete model is discretized in the time domain using the Yee cell structure. The time-marching formulas of the electric and magnetic fields are constructed based on the Maxwell curl equations to calculate the electromagnetic field distribution. An adaptive steady-state judgment criterion is constructed based on the instantaneous total energy change rate, and the FDTD explicit method is used to update the electromagnetic field distribution of the spatial discrete model. The calculation is terminated when the electromagnetic field reaches a steady state. When the electromagnetic field reaches a steady state, the electromagnetic loss power density of the spatial discrete model is calculated and passed into the thermal field calculation; Based on the FVTD method, the thermal field distribution of the discrete model is calculated by explicit or implicit format to obtain the thermal field distribution information.

2. The electromagnetic thermal multi-physics field simulation method based on the FDTD / FVTD algorithm according to claim 1, characterized in that: The steps of establishing a spatial discrete model of an isotropic medium in a passive region, performing time-domain finite-difference discretization on the established spatial discrete model using a Yee cell structure, constructing a time-marching formula for the electric field and the magnetic field according to the Maxwell curl equation, and calculating the electromagnetic field distribution are as follows: Construct a three-dimensional Yee cell structure and perform time-domain finite-difference discretization on the three-dimensional Yee cell structure, with the electric field component arranged on the cell edges and the magnetic field component arranged at the center of the surface; The Maxwell curl equations are discretized by central difference to obtain the propulsion formulas for the electric and magnetic fields. The specific formulas are as follows: in, Indicates that at time n, The electric field E at the grid point x Component value; Indicates time, The magnetic field H at the grid point z Component value; CA(m) is The coefficients of the equation corresponding to the electromagnetic parameters at the grid point, CB(m) are The coefficients of the equations corresponding to the electromagnetic parameters at the grid points; σ(m) is the grid point The conductivity at the grid point, ε(m) is The dielectric constant at .

3. The electromagnetic thermal multi-physics field simulation method based on FDTD / FVTD algorithm according to claim 2, characterized in that: The Courant stability condition is used to determine the time step when performing central difference discretization of Maxwell's curl equations. The specific formula is as follows: Where Δx, Δy, and Δz represent the grid space steps in the x, y, and z coordinate directions, respectively; Δt represents the time step; and c represents the propagation speed of electromagnetic waves in a specific medium.

4. The electromagnetic thermal multi-physics field simulation method based on FDTD / FVTD algorithm according to claim 1, characterized in that: In the step of constructing an adaptive steady-state judgment criterion based on the instantaneous total energy change rate, updating the electromagnetic field distribution of the spatial discrete model using the FDTD explicit method, and terminating the calculation when the electromagnetic field reaches a steady state, the adaptive steady-state judgment criterion is as follows: For scattering problems, the electromagnetic field is calculated by the transient energy within the closed surface S; For radiation problems, the electromagnetic field is calculated by adding the power term of the excitation source; Set the stable convergence parameter and terminate the calculation when the relative wave pair energy of adjacent periods is less than the stable convergence parameter.

5. The electromagnetic thermal multi-physics field simulation method based on FDTD / FVTD algorithm according to claim 4, characterized in that: For scattering problems, the formula for calculating the electromagnetic field through the transient energy in the closed surface S is expressed as follows: Where n is the index of the time step, i is the index variable, M represents the total number of discrete surface elements, subscript m represents the index of the discrete surface element, E m (iΔt) is the electric field vector of the mth discrete element at the moment (iΔt), H m (iΔt) is the magnetic field vector of the mth discrete element at the moment (iΔt), Δs m represents the area vector associated with the mth discrete surface element.

6. The electromagnetic thermal multi-physics field simulation method based on FDTD / FVTD algorithm according to claim 5, characterized in that: For the radiation problem, the formula for calculating the electromagnetic field by adding the excitation source power term is expressed as follows: Where N and N′ represent the number of grid edges where the excitation electric field and excitation magnetic field are located, respectively, and n and n′ are the corresponding indices.

7. The electromagnetic thermal multi-physics field simulation method based on FDTD / FVTD algorithm according to claim 1, characterized in that: When the electromagnetic field reaches a steady state, the electromagnetic loss power density of the spatial discrete model is calculated and passed into the thermal field calculation step. The total electromagnetic loss power density of the electromagnetic field in the medium can be expressed as: in, is the conductive loss, also known as Joule heat, σ is the conductivity of the medium, and |E| is the amplitude of the electric field strength; is the polarization loss, ω represents the angular frequency of the electromagnetic wave, ε0 is the vacuum dielectric constant, and ε″ is the dielectric loss factor of the medium; is the magnetization loss, μ0 is the vacuum permeability, and μ″ is the magnetic loss factor of the medium.

8. The electromagnetic thermal multi-physics field simulation method based on FDTD / FVTD algorithm according to claim 1, characterized in that: In the step of obtaining thermal field distribution information by calculating the thermal field distribution of the discrete model in an explicit or implicit format based on the FVTD method, the law of conservation of energy is applied to the cells. The specific formula is as follows: in are the internal energy of the cell at time t1 and t0 respectively, and are the heat flowing from the six surfaces of the cell from t0 to t1, is the heat generated by the heat source in the cell from t0 to t1.

9. The electromagnetic thermal multi-physics field simulation method based on FDTD / FVTD algorithm according to claim 8, characterized in that: The stable biconjugate gradient method is used to calculate the thermal field distribution of the implicit format discrete model.

10. The electromagnetic thermal multi-physics field simulation system based on FDTD / FVTD algorithm according to any one of claims 1 to 9, characterized in that: include: The first module is used to establish a spatial discrete model of isotropic media in the passive region. The Yee cell structure is used to perform time-domain finite-difference discretization on the established spatial discrete model. The time-marching formula of the electric and magnetic fields is constructed based on the Maxwell curl equation to calculate the electromagnetic field distribution. The second module is used to construct an adaptive steady-state judgment criterion based on the instantaneous total energy change rate, and uses the FDTD explicit method to update the electromagnetic field distribution of the spatial discrete model. The calculation is terminated when the electromagnetic field reaches a steady state. The third module is used to calculate the electromagnetic loss power density of the spatial discrete model after the electromagnetic field reaches a steady state, and pass it into the thermal field calculation; The fourth module is used to calculate the thermal field distribution of the discrete model based on the FVTD method through explicit or implicit format to obtain thermal field distribution information.

Images