CN104679920B - Waveguide device microwave gas discharge method for numerical simulation - Google Patents

Abstract

The invention discloses a microwave gas discharge numerical simulation method of a waveguide device. The method establishes a fluid model of a microwave gas discharge in a waveguide structure, and performs accurate numerical simulations of electron density, electron fluid velocity and average electron energy in the process of gas breakdown. Through the analysis, the law of electric field intensity and electron density at different observation points with time can be obtained in the case of gas discharge in complex structure waveguide devices with different geometric shapes. Based on the time-domain spectral element method, the present invention can realize the modeling of waveguide devices with arbitrary complex structures; the generated mass matrix is a block diagonal matrix, and the inverse matrix can be directly obtained, which greatly reduces the calculation time; at the same time, the method of numerical simulation can Avoid repeated experiments to obtain the internal electric field distribution and electron density of waveguide devices, obtain the breakdown threshold, shorten the design cycle, and save design costs. It is very suitable for effective simulation design and numerical simulation of complex waveguide devices with discontinuous structures .

Description

Numerical Simulation Method of Microwave Gas Discharge in Waveguide Devices

technical field

The invention belongs to the simulation technology of electromagnetic characteristics involving multiple physical fields, in particular to a novel numerical simulation method for analyzing gas discharge breakdown of waveguide structure microwave devices.

Background technique

With the rapid development of space technology, the circuit integration of electronic equipment is getting higher and higher, and the operating frequency of microwave devices is constantly increasing; at the same time, with the continuous maturity of high-power pulse generator technology, its output power can reach the order of tens of Gw. The numerical analysis and computer simulation of the microwave breakdown and gas discharge phenomena of microwave devices with different gas-filled space waveguide structures are increasingly valued by scientific research institutions, and it is one of the hot problems in the field of space technology and high-power pulses. The research on this problem is helpful to the design and protection of microwave devices, and has high application value and practical significance. At present, the simulation analysis of such problems by numerical methods has the advantages of cost minimization, design optimization, cycle shortening, etc., and has become an important means of pre-production design and predictive simulation. When using numerical methods to simulate and analyze this type of multi-physics and complex nonlinear problems, how to efficiently perform fast full-wave analysis and accurately obtain key design thresholds is very important. Compared with the long period and high cost of experimental methods, the prediction by numerical methods faces difficulties such as complex structure of the solution object, high nonlinearity, long calculation time and low calculation accuracy.

At present, the analysis method of microwave breakdown and gas discharge mainly adopts the time-domain finite difference method in the differential equation method, such as literature 1. Yee, Jick H., et al, "Propagation of intense microwave pulses in air and in a waveguide, "IEEE Transactions on Antennas and Propagation, 39.9 (1991), pp:1421-1427. and finite volume integration method, such as literature 2. Adnane Hamiaz, et al," Finite VolumeTime Domain modeling of microwave breakdown and plasma formation in metallic aperture, "Computer Physics Communications 183(2012), pp:1634–1640. The above-mentioned published documents are all used to analyze the two-dimensional gas discharge problem with simple structure. Although the numerical distribution and prediction of the electric field and electron density of the simple problem can be given, but There are disadvantages such as simple model, low accuracy and long calculation time. The time-domain spectral element method uses a hexahedral unit discrete grid, which can well approximate the shape and structure of the target to be solved. The generated mass matrix is a diagonal matrix. The inversion process is simple and the calculation time is short. When analyzing microwave breakdown gas It has great advantages in multiphysics and nonlinear problems such as electric discharge. At present, there are no literature reports on the simulation of gas discharge phenomena in waveguide structures using the time-domain spectral element method.

Contents of the invention

The object of the present invention is to provide a numerical method for simulating the gas breakdown phenomenon in a complex microwave device with a waveguide structure based on the time-domain spectral element method. This method considers the interaction between the gaseous medium and the strong electromagnetic pulse in the waveguide device, It mainly includes avalanche ionization, electron-molecular collision, electron-ion recombination effect, and electron-molecular attachment effect, which can simulate and simplify the discharge physical process in waveguide devices. At the same time, it greatly shortens the development cycle and reduces research and development costs. The invention can be applied to various waveguides Threshold analysis and simulation of breakdown phenomena in devices.

The technical scheme that realizes the object of the present invention is:

The first step is to use Ansys software to establish the geometric dissection model of the device to be analyzed. According to the geometric dimensions of the complex waveguide structure, use computer-aided design tools to model, and use the hexahedron element based on the GLL basis function to dissect the target model to obtain The geometric information of the target and set the position information of the excitation source and apply the excitation source;

In the second step, the initial state of the numerical simulation is given through the given atmospheric pressure, and the initial value of the electron density in the waveguide device is set;

In the third step, according to the geometric information of the device to be analyzed obtained in the first step and the initial state set in the second step, a set of differential equations describing the nonlinear gas discharge phenomenon is established by using the time-domain spectral element method, and the selected The GLL basis function approximates the five unknown variables of the unknown electric field, magnetic field, electron density, electron velocity, and average electron energy, and then substitutes them into the partial differential equations composed of Maxwell's equations and electron hydrodynamic equations, and finally selects the GLL basis The function is used as a weighting function, through the Galerkin test, so that the margin of each differential equation is zero in the sense of weighted average, thereby converting the continuous differential equations into a matrix equation, and obtaining the unknown information of each GLL point in the space ;

In the fourth step, according to the matrix equations obtained in the third step, the central difference method is used to calculate the magnetic field strength, electron velocity, electric field strength, electron density, and average electron energy of each node in the space at each moment, and obtain the The regularity of the above variables over time;

The fifth step is to update the ionization parameters according to the average electron energy obtained in the fourth step, and repeat steps 3 to 4 until the preset calculation time is reached;

The sixth step is to output the electromagnetic field value and electron density calculated at each moment, and obtain the time-varying law of the electromagnetic field and electron density in the waveguide device. From this, various electromagnetic characteristic parameters can be calculated, and the microwave gas discharge phenomenon of the waveguide device can be completed. analysis process.

Compared with the prior art, the present invention has the remarkable advantages of:

(1) The fluid model of microwave breakdown gas discharge in the waveguide device was established, and the numerical simulation calculation was carried out for the electromagnetic field, electron density, electron velocity, and average electron energy of each node inside the waveguide device in the gas discharge, which can be given quantitatively For the unknown variable at any point to be sought, using the waveguide device microwave gas discharge numerical simulation method of the present invention can avoid obtaining the gas breakdown threshold of the waveguide device through repeated experiments, shorten the design cycle, save design costs, and realize Effective protection simulation and pre-production design of complex waveguide devices with structures;

(2) A numerical simulation method for microwave gas discharge of a waveguide device proposed by the present invention is based on the time-domain spectral element method, and the hexahedron subdivision unit can be used to model waveguide devices with any complex structure, so as to realize the geometric information of three-dimensional complex structures Modeling and discretization; at the same time, the mass matrix generated by the spectral element method can be directly inverted as a block diagonal matrix, which greatly reduces the calculation time, reduces the calculation cost, and shortens the design cycle;

(3) In the microwave gas discharge numerical model of the waveguide device proposed by the present invention, the interaction between the gaseous medium in the waveguide device and the strong electromagnetic pulse is considered, mainly including avalanche ionization, electron-molecular collision, electron-ion recombination effect, and electron-molecular adhesion effect. It can simulate the discharge physical process in the simplified waveguide device, and adopts the improved electron energy distribution function solution method. At each time step, according to the obtained average electron energy, the ionization parameters of each point are updated, and the Maxwell equation and the electron are solved. The strong coupling of fluid mechanics equations ensures the accuracy of ionization parameter updates and improves the accuracy of fluid models.

Description of drawings

Figure 1 is a schematic cross-sectional view of a waveguide structure with a slot.

Fig. 2 is a diagram showing the change of electric field intensity with time at the observation point at the midpoint of the first square cone.

Fig. 3 is a graph showing the change of electron density with time at the observed point at the midpoint of the first square cone.

Fig. 4 is a flowchart of the present invention.

detailed description

The present invention will be described in further detail below in conjunction with the accompanying drawings.

The invention provides a numerical simulation method for microwave gas discharge of a waveguide device. The specific steps of the present invention will be further described in detail below in conjunction with the accompanying drawings, taking the device with a waveguide structure with a slot shown in FIG. 1 as an example.

See Figure 1 for a schematic cross-sectional view of the waveguide structure with slots. The geometric dimensions of the model are as follows: BJ-100 waveguide model is used, the interior is filled with argon gas, the pressure is 3.9Torr, the long side of the waveguide a=22.86mm, the short side of the waveguide b=10.16mm, The total length of the waveguide l=27.45mm, the center distance of the resonant element is 9.15mm, the discharge electrodes are a pair of cuboids, the side length a1=1.0mm; the diaphragm width m=8.0mm, the height b=10.16mm, the thickness is 1.0mm, the middle The gap height h1=1.0mm. Analyzing the microwave gas discharge numerical simulation method of band slot waveguide structure shown in Fig. 1 according to the present invention, its specific operation steps are as follows:

In the first step, according to the geometric dimensions of the waveguide structure with slots shown in Figure 1, the Ansys software is used to model it, and the hexahedral element based on the GLL basis function is used to divide the waveguide structure with slots, that is, the waveguide structure with slots is obtained The number of unknown quantities and the coordinate information of the nodes, as for the boundary information of the waveguide structure with slot, can be dealt with by setting the four-sided rectangular waveguide wall as ideal metal and the two sides of the transmission direction as absorbing boundary conditions according to the structure and size of the waveguide with slot. And set the location information of the excitation source;

In the second step, according to the atmospheric pressure of the given working environment, the initial state of the numerical simulation is given, and the initial value of the electron density in the waveguide device is set to 10 ⁶ m ^-3 ;

In the third step, according to the geometric information of the waveguide structure with slots and the initial state set in the second step, the time-domain spectral element method is used to establish the differential equations describing the nonlinear gas discharge phenomenon, that is, the selected GLL basis function is unknown The five unknown variables of electric field, magnetic field, electron density, electron velocity and average electron energy are approximately expanded, and then respectively substituted into the partial differential equations composed of Maxwell's equations and electron hydrodynamic equations. Finally, the GLL basis function is selected as the weighting function. Through the Galerkin test, the margin of the equation in the sense of weighted average is zero, thus converting the continuous differential equations into matrix equations, as shown in formula (1):

Among them, T represents the mass matrix, S represents the stiffness matrix, and the subscripts of the mass matrix and the stiffness matrix correspond to the expansion basis and the test basis respectively. R and F represent boundary integral matrix and source vector matrix, and ν _i , ν _a , ν _c and Q _l are ionization parameters representing ionization rate, attachment rate, collision rate and energy loss rate, respectively;

The fourth step is to calculate the magnetic field strength H, electron velocity u, electric field strength E, electron density n, and average electron energy of each node in the space at each moment using the central difference method based on the matrix equations obtained in the third step The law of the above variables changing with time at different positions is obtained. The display format of each time step iteration of each variable is shown in the following formula (2):

Among them, Δt represents the time step, and the superscript k represents the number of time steps;

The fifth step is to update the ionization parameters according to the average electron energy obtained in the fourth step:

1. Obtain the ionization rate according to the following formula (3):

where ν _i is the ionization rate, N _gas is the gas number density, m _e is the electron mass, ε _e is the electron energy, σ _i is the ionization cross-sectional area, and f(ε _e ) is the electron energy distribution function;

2. The electron energy distribution function required to solve the integral term in the above formula is expressed as the average electron energy The function:

in ξ ₁ = 3/(2χ), ξ ₂ = 5/(2χ), χ = 6.5 corresponds to the case of argon filling,

3. According to the above-mentioned functional relationship between the ionization parameters and the average electron energy and the average electron energy value calculated in the fourth step, use the interpolation method to obtain the ionization parameters of each position point in the updated space, and substitute them for the iterative calculation of the next cycle , repeat steps 3 to 4 until the preset calculation time is reached;

The sixth step is to output the calculated electric field value and electron density at each moment, and obtain the time-varying law of electric field strength and electron density in the device with a waveguide structure with a gap.

According to the method of the present invention, the device with a waveguide structure with a slot is simulated, and the electric field intensity and electron density at the observation point at the midpoint of the first square cone change with time as shown in Fig. 2 and Fig. 3 . It can be seen from the figure that at the midpoint of a pair of square cones, the field strength increases sharply, and when the field strength threshold of gas breakdown is reached, the electron concentration increases exponentially. With the increase of electron concentration, plasma is generated at the center of the square cone, making The incident electric field attenuates, and when the electron concentration increases from 10 ⁶ m ^-3 to 2.12×10 ¹⁹ m ^-3 it tends to be stable, and finally the electron concentration of the generated plasma is 2.12×10 ¹⁹ m ^-3 . The simulation results can well explain a series of physical phenomena such as the surge in electron concentration and the attenuation of the incident electric field after the microwave gas breakdown of the complex waveguide structure with gaps. The changing trend of exponential increase can quantitatively give the changing law of unknown variables at the observation point, which has high guiding significance for the design and protection of waveguide structure microwave devices.

Claims (3)

1. A waveguide device microwave gas discharge numerical simulation method is characterized in that the steps are as follows: In the first step, according to the geometric size of the complex waveguide structure, the geometric subdivision model of the target is established, and the target model is subdivided by the hexahedral element based on the Guass-Lobatto-Legendre basis function to obtain the geometric information of the target, and then the excitation source is set. position information and apply stimulus; The second step is to set the initial value of the electron density in the waveguide device through the given atmospheric pressure, and give the initial state of the microwave gas discharge numerical simulation of the waveguide device; In the third step, according to the geometric information of the target obtained in the first step and the initial value set in the second step, the time-domain spectral element method is used to establish the differential equations describing the nonlinear gas discharge phenomenon, that is, the selected GLL base Function for unknown magnetic field strength H, electron velocity u, electric field strength E, electron density n, average electron energy The five unknown variables are approximately expanded, and then respectively substituted into the partial differential equations composed of Maxwell's equations and electrohydrodynamic equations. Finally, the GLL basis function is selected as the weighting function, and through the Galerkin test, each equation in the sense of weighted average The margin of is zero, thus converting the continuous differential equations into matrix equations, as shown in the following formula (1): <mrow><mfenced open = "{" close = ""><mtable><mtr><mtd><mrow><msub><mi>T</mi><mrow><mi>h</mi><mi>i</mi><mi>h</mi><mi>j</mi></mrow></msub><mfrac><mrow><mi>d</mi><mi>H</mi></mrow><mrow><mi>d</mi><mi>t</mi></mrow></mfrac><mo>+</mo><msub><mi>S</mi><mrow><mi>h</mi><mi>i</mi><mi>e</mi><mi>j</mi></mrow></msub><mo>&amp;CenterDot;</mo><mi>E</mi><mo>=</mo><mn>0</mn></mrow></mtd></mtr><mtr><mtd><mrow><msub><mi>T</mi><mrow><mi>e</mi><mi>i</mi><mi>e</mi><mi>j</mi></mrow></msub><mfrac><mrow><mi>d</mi><mi>E</mi></mrow><mrow><mi>d</mi><mi>t</mi></mrow></mfrac><mo>+</mo><msub><mi>S</mi><mrow><mi>e</mi><mi>i</mi><mi>n</mi><mi>j</mi><mi>u</mi><mi>k</mi></mrow></msub><mo>&amp;CenterDot;</mo><mi>n</mi><mi>u</mi><mo>-</mo><msub><mi>S</mi><mrow><mi>e</mi><mi>i</mi><mi>h</mi><mi>j</mi></mrow></msub><mi>H</mi><mo>+</mo><msub><mi>R</mi><mrow><mi>e</mi><mi>i</mi><mi>e</mi><mi>j</mi></mrow></msub><mi>E</mi><mo>=</mi>mo><mn>0</mn></mrow></mtd></mtr><mtr><mtd><mrow><mfrac><mrow><mi>d</mi><mi>n</mi></mrow><mrow><mi>d</mi><mi>t</mi></mrow></mfrac><mo>-</mo><mrow><mo>(</mo><msub><mi>v</mi><mi>i</mi></msub><mo>-</mo><msub><mi>v</mi><mi>a</mi></msub><mo>)</mo></mrow><mi>n</mi><mo>=</mo><mn>0</mn></mrow></mtd></mtr><mtr><mtd><mrow><msub><mi>T</mi><mrow><mi>u</mi><mi>i</mi><mi>u</mi><mi>j</mi></mrow></msub><mfrac><mrow><mi>d</mi><mi>u</mi></mrow><mrow><mi>d</mi><mi>t</mi></mrow></mfrac><mo>+</mo><mrow><mo>(</mo><msub><mi>v</mi><mi>i</mi></msub><mo>+</mo><msub><mi>v</mi><mi>c</mi></msub><mo>)</mo></mrow><msub><mi>T</mi><mrow><mi>u</mi><mi>i</mi><mi>u</mi><mi>j</mi></mrow></msub><mi>u</mi><mo>-</mo><msub><mi>S</mi><mrow><mi>u</mi><mi>i</mi><mi>e</mi><mi>j</mi></mrow></msub><mi>E</mi><mo>=</mo><mn>0</mn></mrow></mtd></mtr><mtr><mtd><mrow><msub><mi>T</mi><mrow><mi>U</mi><mi>i</mi><mi>U</mi><mi>j</mi></mrow></msub><mfrac><mrow><mi>d</mi><msub><mover><mi>&amp;epsiv;</mi><mo>&amp;OverBar;</mo></mover><mi>e</mi></msub></mrow><mrow><mi>d</mi><mi>t</mi></mrow></mfrac><mo>+</mo><msub><mi>v</mi><mi>i</mi></msub><msub><mi>T</mi><mrow><mi>U</mi><mi>i</mi><mi>U</mi><mi>j</mi></mrow></msub><msub><mover><mi>&amp;epsiv;</mi><mo>&amp;OverBar;</mo></mover><mi>e</mi></msub><mo>-</mo><msub><mi>S</mi><mrow><mi>U</mi><mi>i</mi><mi>e</mi><mi>j</mi><mi>u</mi><mi>k</mi></mrow></mi>msub><mi>E</mi><mi>u</mi><mo>+</mo><mi>F</mi><mo>&amp;CenterDot;</mo><msub><mi>Q</mi><mi>l</mi></msub><mo>=</mo><mn>0</mn></mrow></mtd></mtr></mtable></mfenced><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>1</mn><mo>)</mo></mrow></mrow> Among them, T represents the mass matrix, S represents the stiffness matrix, and the subscripts of the mass matrix and the stiffness matrix correspond to the expansion basis and the test basis respectively; R and F represent the boundary integral matrix and source vector matrix, ν _i , ν _a , ν _c and _Ql are ionization parameters representing ionization rate, attachment rate, collision rate and energy loss rate, respectively; In the fourth step, according to the matrix equations (1) obtained in the third step, the central difference method is used to calculate the magnetic field strength, electron velocity, electric field strength, electron density, and average electron energy of each node in the space at each moment, and obtain The law of the above variables changing with time at different locations; ① Expand the ordinary differential equations shown in formula (1) using the central difference scheme, transform them into the form of algebraic equations, and finally give explicit iterative formulas for five unknown quantities, as shown in formula (2): <mrow><mfenced open = "{" close = ""><mtable><mtr><mtd><mrow><msup><mi>H</mi><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msup><mo>=</mo><msup><mi>H</mi><mi>k</mi></msup><mo>-</mo><msup><msub><mi>T</mi><mrow><mi>h</mi><mi>i</mi><mi>h</mi><mi>j</mi></mrow></msub><mrow><mo>-</mo><mn>1</mn></mrow></msup><msub><mi>S</mi><mrow><mi>h</mi><mi>i</mi><mi>e</mi><mi>j</mi></mrow></msub><mo>&amp;CenterDot;</mo><msup><mi>E</mi><mi>k</mi></msup><mo>&amp;CenterDot;</mo><mi>&amp;Delta;</mi><mi>t</mi></mrow></mtd></mtr><mtr><mtd><mrow><msup><mi>E</mi><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msup><mo>=</mo><msup><mi>E</mi><mi>k</mi></msup><mo>+</mo><msup><msub><mi>T</mi><mrow><mi>e</mi><mi>i</mi><mi>e</mi><mi>j</mi></mrow></msub><mrow><mo>-</mo><mn>1</mn></mrow></msup><mrow><mo>(</mo><mrow><msub><mi>S</mi><mrow><mi>e</mi><mi>i</mi><mi>n</mi><mi>j</mi><mi>u</mi><mi>k</mi></mrow></msub><mo>&amp;CenterDot;</mo><msup><mi>n</mi><mi>k</mi></msup><msup><mi>u</mi><mi>k</mi></msup><mo>-</mo><msub><mi>S</mi><mrow><mi>e</mi><mi>i</mi><mi>h</mi><mi>j</mi></mrow></msub><mo>&amp;CenterDot;</mo><msup><mi>H</mi><mi>k</mi></msup><mo>+</mo><msub><mi>R</mi><mrow><mi>e</mi><mi>i</mi><mi>e</mi><mi>j</mi></mrow></msub><msup><mi>E</mi><mi>k</mi></msup></mrow><mo>)</mo></mrow><mi>&amp;Delta;</mi><mi>t</mi></mrow></mtd></mtr><mtr><mtd><mrow><msup><mi>n</mi><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msup><mo>=</mo><mrow><mo>&amp;lsqb;</mo><mrow><mn>1</mn><mo>+</mo><mi>&amp;Delta;</mi><mi>t</mi><mrow><mo>(</mo><mrow><msub><mi>v</mi><mi>i</mi></msub><mo>-</mo><msub><mi>v</mi><mi>a</mi></msub></mrow><mo>)</mo></mrow></mrow><mo>&amp;rsqb;</mo></mrow><msup><mi>n</mi><mi>k</mi></msup></mrow></mtd></mtr><mtr><mtd><mrow><msup><mi>u</mi><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msup><mo>=</mo><mrow><mo>&amp;lsqb;</mo><mrow><mn>1</mn><mo>-</mo><mrow><mo>(</mo><mrow><msub><mi>v</mi><mi>i</mi></msub><mo>+</mo><msub><mi>v</mi><mi>c</mi></msub></mrow><mo>)</mo></mrow><mi>&amp;Delta;</mi><mi>t</mi></mrow><mo>&amp;rsqb;</mo></mrow><msup><mi>u</mi><mi>k</mi></msup><mo>+</mo><msubsup><mi>T</mi><mrow><mi>u</mi><mi>i</mi><mi>u</mi><mi>j</mi></mrow><mrow><mo>-</mo><mn>1</mn></mrow></msubsup><msub><mi>S</mi><mrow><mi>u</mi><mi>i</mi><mi>e</mi><mi>j</mi></mrow></msub><msup><mi>E</mi><mi>k</mi></msup><mi>&amp;Delta;</mi><mi>t</mi></mrow></mtd></mtr><mtr><mtd><mrow><msup><msub><mover><mi>&amp;epsiv;</mi><mo>&amp;OverBar;</mo></mover><mi>e</mi></msub><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msup><mo>=</mo><mrow><mo>&amp;lsqb;</mo><mrow><mn>1</mn><mo>-</mo><msub><mi>v</mi><mi>i</mi></msub><mi>&amp;Delta;</mi><mi>t</mi></mrow><mo>&amp;rsqb;</mo></mrow><msup><msub><mover><mi>&amp;epsiv;</mi><mo>&amp;OverBar;</mo></mover><mi>e</mi></msub><mi>k</mi></msup><mo>+</mo><msup><msub><mi>T</mi><mrow><mi>U</mi><mi>i</mi><mi>U</mi><mi>j</mi></mrow></msub><mrow><mo>-</mo><mn>1</mn></mrow></msup><mrow><mo>(</mo><mrow><msub><mi>S</mi><mrow><mi>U</mi><mi>i</mi><mi>e</mi><mi>j</mi><mi>u</mi><mi>k</mi></mrow></msub><msup><mi>E</mi><mi>k</mi></msup><mo>&amp;CenterDot;</mo><msup><mi>u</mi><mi>k</mi></msup><mo>-</mo><mi>F</mi><mo>&amp;CenterDot;</mo><msub><mi>Q</mi><mi>l</mi></msub></mrow><mo>)</mo></mrow><mi>&amp;Delta;</mi><mi>t</mi></mrow></mtd></mtr></mtable></mfenced><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>2</mn><mo>)</mo></mrow></mrow> Among them, Δt represents the time step, and the superscript k represents the number of time steps; ② At each time step, the value obtained in the previous time step is substituted into the above explicit iteration formula (2) in order to obtain the variable value of the next time step and save it; The fifth step is to update the ionization parameters according to the average electron energy obtained in the fourth step, and repeat steps 3 to 4 until the preset calculation time is reached; ① Obtain the ionization rate according to the following formula (3): <mrow><msub><mi>v</mi><mi>i</mi></msub><mo>=</mo><msub><mi>N</mi><mrow><mi>g</mi><mi>a</mi><mi>s</mi></mrow></msub><mrow><mo>&amp;lsqb;</mo><mrow><munderover><mo>&amp;Integral;</mo><mn>0</mn><mi>&amp;infin;</mi></munderover><msqrt><mfrac><mrow><mn>2</mn><msub><mi>&amp;epsiv;</mi><mi>e</mi></msub></mrow><msub><mi>m</mi><mi>e</mi></msub></mfrac></msqrt><msub><mi>&amp;sigma;</mi><mi>i</mi></msub><mi>f</mi><mrow><mo>(</mo><msub><mi>&amp;epsiv;</mi><mi>e</mi></msub><mo>)</mo></mrow><msub><mi>d&amp;epsiv;</mi><mi>e</mi></msub></mrow><mo>&amp;rsqb;</mo></mrow><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>3</mn><mo>)</mo></mrow></mrow> where ν _i is the ionization rate, N _gas is the gas number density, m _e is the electron mass, ε _e is the electron energy, σ _i is the ionization cross-sectional area, and f(ε _e ) is the electron energy distribution function; ② The electron energy distribution function required to solve the integral term in the above formula is expressed as the average electron energy The function: <mrow><mi>f</mi><mrow><mo>(</mo><msub><mi>&amp;epsiv;</mi><mi>e</mi></msub><mo>)</mo></mrow><mo>=</mo><msub><mi>c</mi><mn>1</mn></msub><msup><msub><mi>&amp;epsiv;</mi><mi>e</mi></msub><mfrac><mn>1</mn><mn>2</mn></mfrac></msup><mi>exp</mi><mrow><mo>(</mo><mrow><mo>-</mo><msub><mi>c</mi><mn>2</mn></msub><msup><msub><mi>&amp;epsiv;</mi><mi>e</mi></msub><mi>&amp;chi;</mi></msup></mrow><mo>)</mo></mrow><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>4</mn><mo>)</mo></mrow></mrow> in ξ ₁ =3/(2χ),ξ ₂ =5/(2χ), χ is determined by the gas composition in the waveguide device; ③ According to the above-mentioned functional relationship between the ionization parameter and the average electron energy and the average electron energy value calculated in the fourth step, use the interpolation method to obtain the ionization parameters of each position point in the updated space, and substitute them for the iterative calculation of the next cycle; The sixth step is to output the calculated magnetic field strength, electric field strength and electron density at each moment, and obtain the time-varying law of the electromagnetic field and electron density in the waveguide device, from which various electromagnetic characteristic parameters can be calculated. 2. The waveguide device microwave gas discharge numerical simulation method according to claim 1, characterized in that: in the first step, the geometric information of the target refers to the number of unknown quantities and the coordinate information and boundary information of nodes. 3. waveguide device microwave gas discharge numerical simulation method according to claim 1 is characterized in that: in the second step, by given atmospheric pressure p (Torr), the initial value of electron density in the waveguide device is set