Zhihang ZHAO(赵志航),Xinlao WEI(魏新劳),Shuang SONG(宋爽),Lin CUI(崔林),Kailun YANG(杨凯伦) and Zhonghua ZHANG(张中华)
1 Key Laboratory of Engineering Dielectrics and Application of Ministry of Education,Harbin University of Science and Technology,Harbin 150080,People’s Republic of China
2 No.703 Research Institute of CSIC(China Shipbuilding Industry Corporation),Harbin 150078,People’s Republic of China
3 Yunnan Electric Test&Research Institute Group Co.Ltd,Kunming 650217,People’s Republic of China
Abstract In this paper,an improved air discharge fluid model under non-uniform electric field is constructed based on the plasma module COMSOL Multiphysics with artificial stability term,and the boundary conditions developed in the previous paper are applied to the calculation of photoionization rate.Based on the modified model,the characteristics of low temperature subatmospheric air discharge under 13 kV direct current voltage are discussed,including needle-plate and needle-needle electrode structures.Firstly,in order to verify the reliability of the model,a numerical example and an experimental verification were carried out for the modified model respectively.Both verification results show that the model can ensure the accuracy and repeatability of the calculation.Secondly,according to the calculation results of the modified model,under the same voltage and spacing,the reduced electric field under low temperature subatmosphere pressure is larger than that under normal temperature and atmospheric pressure.The high electric field leads to the air discharge at low temperature and sub atmospheric pressure entering the streamer initiation stage earlier,and has a faster propagation speed in the streamer development stage,which shortens the overall discharge time.Finally,the discharge characteristics of the two electrode structures are compared,and it is found that the biggest difference between them is that there is a pre-ionization region near the cathode in the needle-needle electrode structure.When the pre-ionization level reaches 1013 cm−3,the propagation speed of the positive streamer remains unchanged throughout the discharge process,and is no longer affected by the negative streamer.The peak value of electric field decreases with the increase of pre-ionization level,and tends to be constant during streamer propagation.Based on the previous paper,this paper constructs the air discharge model under non-uniform electric field,complements with the previous paper,and forms a relatively complete set of air discharge simulation system under low temperature and sub atmospheric pressure,which provides a certain reference for future research.
Keywords:low temperature and sub atmospheric pressure,artificial stability term,reduced electric field,pre-ionization,simulation system
Streamer discharge,as one of the main forms of gas discharge,has always been a fundamental problem in electrical engineering.With the development of science and technology,the relationship between disciplines is closer,which makes the gas discharge involved in more and more fields.As the author mentioned in[1],with the further development of aviation industry,on the premise of safe operation,it has become an urgent problem for aviation and high-voltage engineering to study the air discharge characteristics of aircraft in low-temperature sub atmospheric pressure environment,which has certain significance for promoting the development of aviation technology.
Streamer discharge is a transient process with drastic changes.So far,experiment is still an important method to study the streamer discharge process[2–5].However,the existing experimental methods and existing experimental data still cannot fully clarify the micro-mechanism of the discharge process,and it is difficult to obtain the net charge density,electric field distribution,streamer radius and other microscopic parameters in the discharge channel.Therefore,numerical simulation has become an important method to promote the development of gas discharge theory.By comparing the similarities and differences between the simulation results and the experimental results,optimizing the selection of parameters,and improving the calculation model,it is helpful to establish a streamer discharge model that is closer to the real situation,to deepen the understanding of the discharge process and better promote the development of gas discharge research.
With the improvement of computer technology,it is possible to realize the numerical simulation of gas discharge physical process by using high-performance computer and optimized numerical algorithm.The preliminary study shows that the correct numerical simulation can obtain the results basically consistent with the experimental results[6,7].Streamer discharge is a multiscale system,including many processes,such as chemical reaction,component transport and fluid coupling.It is a great challenge to simulate the above processes with high accuracy and high efficiency.Firstly,streamer discharge process is highly nonlinear.Due to the interaction between the streamer head and the space charge layer in front of the streamer,high spatial accuracy is required at the junction.In addition,the one-dimensional model cannot accurately describe the arc structure in front of streamer.Therefore,Morrow R[8]and others combined twodimensional Poisson equation with one-dimensional transport equation to build a 1.5-dimensional model to simulate the radial streamer,but this method cannot simulate the discharge of negative polarity streamer and surface streamer.The Schaffer-Gummel(SG)method introduced in[9,10]can simulate the spark and glow discharge in two-dimensional axisymmetric geometry.This method can also analyze the dynamic process of positive and negative streamers[11,12].In[13,14],the improved method was applied to the needle-plate discharge process and the accurate simulation was realized.With the further improvement of Pechereau F[15],this method has been applied to the calculation of streamer medium contact problem.Secondly,there are many spatial scales in the process of streamer propagation.Under atmospheric pressure,the propagation distance of streamer between electrodes is generally centimeter level,while the charge layer thickness of streamer head is only tens of microns.Therefore,the streamer model needs to span four orders of magnitude in spatial scale.The most common solution to this problem is to use unstructured mesh and encrypt the mesh appropriately on the streamer propagation path.In order to further improve the computational efficiency,[16,17]use the latest structured adaptive mesh.The streamer process also has multiple time scales.The propagation time scale of streamer between electrodes is generally nanosecond,and the minimum time scale for simulating streamer transport process at atmospheric pressure is 10−14s.Therefore,the streamer model still needs to span multiple orders of magnitude in time scale.The above difficulties can be overcome by using semi-implicit propulsion[18,19]and asynchronous time propulsion[20].Finally,streamer discharge has strong nonlocality.It is mainly reflected in the distribution of electric field and photoionization,which makes parallel acceleration very difficult.Especially in the case of medium,the change of boundary conditions makes the local electric field have a huge difference.There are three ways to solve the nonlocality problem:one is to use the local energy approximation,the second is to modify the ionization reaction rate coefficient,the third is the hybrid simulation of particle and fluid.Li C[21]and others applied particle model and fluid model to simulate streamer head and streamer body respectively.This method integrates Monte Carlo model and fluid model,and obtains a lot of research in the follow-up[22,23].The particle Monte Carlo model is often used alone[24,25]to study runaway electrons and streamer bifurcation.
In recent years,the commercial finite element software COMSOL Multiphysics has stood out among many streamer simulation software due to its unique advantages.Its advantages are mainly reflected as follows:(1)the simple and easyto-understand operation interface allows beginners to get started quickly and get some initial results;(2)based on the powerful finite element method,it can deal with extremely complex geometric figures;(3)modular settings can help beginners to reduce programming and other complex operations,and the partial differential equation module is more versatile strong,can realize the free organization between models;(4)adaptive mesh can increase the number of mesh to a certain extent;(5)the built-in post-processing system also makes the calculation results clearer and more visible.
Although COMSOL Multiphysics has many advantages,its defects are still very significant.These defects make it impossible to achieve the most ideal tool for plasma simulation under atmospheric pressure.Firstly,the plasma module of COMSOL Multiphysics uses standard Galerkin method,which is suitable for the convection dominated region,but not for the transport dominated region of the streamer head,cannot guarantee the calculation results are true and effective.Secondly,COMSOL Multiphysics with ‘black box’ mode is different from other programmable software,its built-in module cannot be modified by human,so that users can not start when they encounter problems.Finally,COMSOL Multiphysics was unable to make full use of the parallel structure of the computer due to its own reasons,making its calculations slow and memory consumption huge.
Aiming at the above defects,the weak form partial differential equation is introduced into COMSOL Multiphysics software as artificial stability term,as a solution to overcome numerical oscillation and pseudo-diffusion problems in plasma simulation,and the influence of photoionization reaction is considered.The second order scattering boundary condition proposed in[1]is used in photoionization calculation condition.The choice of chemical reaction rate and transport coefficient is an important problem in the theory of air discharge,which is often obtained by determining the electron energy distribution.Figure 1 shows the relationship between the electron energy calculated by the two-term approximation module of COMSOL Multiphysics Boltzmann equation and the electron energy distribution function.Air discharge belongs to non-equilibrium plasma discharge at low temperature,and its electron energy distribution is different from Maxwell distribution and Druyvesteyn distribution.Therefore,it is necessary to accurately calculate the electron energy distribution function to determine whether the model is reliable.In this paper,the electronic energy distribution function is calculated by the open-source electronic energy distribution function calculation software BOLSIG+,which has been widely recognized in the world[26].
Figure 1.Comparison of Maxwellian,Druyvesteyn and two-term Boltzmann distribution functions for different mean electron energies.
This paper is divided into three parts.In the first part,the particle transport equation,Poisson equation and photoionization term used in the model are introduced.Then,the corresponding boundary conditions are introduced,part of the mesh generation is given,and the mesh quality is evaluated.Finally,the input form of artificial stability term in COMSOL Multiphysics is derived in weak form.In the second part,the modified model is verified by examples and experiments to ensure its credibility.Then,the modified model is used to calculate the low temperature sub atmospheric streamer discharge under non-uniform electric field.In the third part,the corresponding conclusions are given.
The model studied in this paper is a two-dimensional axisymmetric configuration based on the COMSOL Multiphysics 5.6®package.First,a relatively complete set of plasma equations was established,including transport equations,Poisson equations and Helmholtz equations.
The electron transport equation is:
whereneis the electron density(cm−3),tis time(s),μeis the electron mobility(cm2V−1s−1),Deis the electron diffusion coefficient(cm2s−1),E is the electric field(kV cm−1),Seis the electron yield dominated by plasma chemical reaction(cm−3s−1),Sphis the electron yield dominated by photoionization(cm−3s−1),photoionization reaction provides seed electrons for streamer discharge,forming secondary electron avalanche and promoting discharge formation.For the convenience of comparison with the benchmark example,the reaction rate and transport coefficient in section 3.1 are calculated by the formula given in[13],but the electron diffusion,electron mobility and reaction rate coefficient used in sections 3.2 and 3.3 are obtained by the calculation software BOLSIG+.
Under the nanosecond discharge time scale,it can be assumed that ions and neutral particles do not move,and only the reaction equation excluding transport term can be solved:
whereniis the ion density(cm−3),divided intoni+andni−,Diis the ion diffusivity(cm2s−1),Siis the ion yield dominated by chemical reactions(cm−3s−1),Sphis the ion yield dominated by photoionization(cm−3s−1).
The Poisson equation is:
whereeis the elementary charge,e=1.602176565×10−19C,ε is the dielectric constant of air.
The Helmholtz equation for calculating the photoionization term is:
wherepO2is the partial pressure of oxygen,which is 150 Torr at standard atmospheric pressure.The parameters λjandAjare consistent with the data in[1].Equation(5)is the new boundary condition ‘second order scattering boundary condition’ developed in the previous paper[1].For the detailed introduction of boundary conditions,please refer to[1],which will not be repeated here.
All the calculations involved in this model are performed on the workstation with dual 2.9 GHz Intel(R)Xeon processors in 64 bit,256 GB ram and windows 7 operating system.The geometry of needle-plate electrode and needle-needle electrode involved in the model is shown in figure 2,where the radius of curvature of the needle is 175 μm.
Figure 2.Simulation geometry of the streamer discharge models.(a)Needle-plate electrode,(b)needle-needle electrode.
Accurate boundary conditions and initial values are the key to solve the equation.The boundary conditions corresponding to figures 2(a)and(b)are given in tables 1 and 2.The discharge gap is 1 cm and the voltage source is 13 kV DC.The ambient temperature is 223 K and the air pressure is 170 Torr,which all meet the external environmental conditions of the aircraft(such as airplane)11 km from the ground.In order to skip the initial stage of streamer discharge,it is assumed that there is an initial plasma(electrons and positive ions)in the calculation region.The electron density and positive ion density show Gaussian distribution in both radial and axial directions[13]:
Table 1.Boundary conditions for needle-plate electrode.
Table 2.Boundary conditions for needle-needle electrode.
wheren0=1014cm−3,σr=0.01 cm,σz=0.025 cm.The value ofz0depends on the electrode position in figure 2.The initial condition of negative ionni−is assumed to be zero.
Figure 3 shows the fineness of the mesh generation of the model through some meshes.In numerical calculation,accurate mesh generation is the guarantee of model convergence.
Figure 3.Mesh generation.
The quality of mesh elements is an important factor in meshing[27].Its value is a scalar value between 0 and 1.The larger the value is,the higher the quality of mesh elements is.According to table 3,the minimum unit mass of the needleplate electrode structure is 0.7645,and the average unit mass is 0.9512.In the needle-needle electrode structure,the minimum unit mass is 0.7821 and the average unit mass is 0.9733.In order to improve the convergence of the model,the adaptive mesh is used based on the basic mesh.For the transient process such as streamer discharge,the adaptive mesh is undoubtedly a very economical and effective method,which can get more accurate results.However,the benefit of this method largely depends on the basic mesh.Since the problem must be solved on the basic mesh for a period,the basic mesh does not need to be too refined under the premise of ensuring the quality of the mesh.
Table 3.Mesh statistics.
In order to better illustrate the weak form of the artificial stability term,this section takes the transport equation as an example to illustrate the derivation process of the weak form in detail.First,we rewrite the electron transport equation(1)into the following form:
By multiplying the ‘test function’uon the left and right sides of the equation and integrating the finite element volume,we can get the following results:
The first two terms of the above formula are integrated by parts,and simplify to obtain a weak form equation that can be directly input into COMSOL Multiphysics:
In the above formula,the first line is the weak form partial differential equation in the plasma computational domain,and the second line represents the boundary conditions naturally satisfied by it.If the drift diffusion process described in equation(9)is dominated by the diffusion process,the standard Galerkin method adopted by COMSOL Multiphysics can obtain reliable and stable solutions.If the solution region is dominated by the drift process,then equation(9)shows the characteristics of hyperbolic equation,and the numerical method is no longer stable.The results of low-quality mesh show violent oscillation,and the pseudo diffusion solution will be formed in high-quality mesh.
Figure 4 uses the built-in plasma module of COMSOL Multiphysics to reproduce the calculation results of the needle-plate discharge in the[13]at two different grid sizes of 50 and 15 μm.In figures 4(a)and(b),50 μm low-quality mesh is used,while in figures 4(c)and(d),15 μm high-quality mesh is used.
Figures 4(a)and(c)compare the axial electron density distribution(blue line)of streamer channel in two mesh sizes using COMSOL Multiphysics built-in plasma module with that in the benchmark example(red line).In the low-quality mesh,as shown in figure 4(b),the results of COMSOL Multiphysics are consistent with the benchmark in streamer propagation speed,but there is a severe nonphysical numerical oscillation in the whole discharge channel as shown in figure 4(a),which makes it impossible to analyze the calculation results in depth.If we only rely on the refined mesh,although there is no numerical oscillation in the calculation results of axial electron density in figure 4(c),we can find that the streamer propagation distance in figure 4(d)lags the benchmark example,which indicates that COMSOL Multiphysics standard finite element method cannot capture the streamer discharge process in the benchmark example,and the calculation results are diffuse streamer.The above two results show that it is impossible to correctly calculate streamer discharge process by using the default numerical format and calculation method of COMSOL Multiphysics.
Figure 4.Failed cases when repeating the benchmark proposed by[13]using COMSOL with conventional methods.Two types of numerical instability appear:(a),(b)oscillations in a coarse mesh and(c),(d)diffusion in a refined mesh.
This numerical instability is mainly due to the mismatch between the computational grid and the algorithm.The standard Galerkin method in the finite element method is the same as the standard SG method in the finite difference method.In the atmospheric pressure range,the accuracy of streamer discharge calculation will be significantly reduced.If the electron mobility and the electron diffusion rate conform to the Einstein relationship,the following relationship should be ensured in order to ensure the accuracy of the calculation results[14]:
wherekrepresents thekth mesh node,ΔEis the electric field gradient(kV cm−1),his the mesh size(cm),Teis the electron temperature(V).The above formula shows that the potential difference between two adjacent nodes should be much smaller than the electron temperature.For the numerical calculation in the atmospheric pressure range,if this condition is satisfied,high mesh density is required,which makes highperformance calculation almost impossible.Unless the calculation area is very small,the mesh conditions specified in equation(10)cannot be realized in COMSOL Multiphysics.Therefore,artificial stability can be introduced to improve the calculation accuracy.Streamline upwind/Petrov Galerkin is one of the most widely used artificial stability methods[28,29].Through the weak formal equation,the artificial stability term is added to each finite element in COMSOL Multiphysics:
In the above formula,the meanings of each item are as follows:
wherenelis the total number of finite elements,nehis the electron density in a finite element of characteristic sizeh,uhis the corresponding test function;Гexand Гeyare the components of the electron flux Г=(−De∇ne+neμeE)in thexandydirections;αxand αyare the parameters related to the Péclet number in thexandydirections,respectively:
The stable solution can be obtained by adding equation(11)to the weak form of electron transport equation in the form of weak form equation.The verification of the modified model will be introduced in detail in section 3.1 through two aspects of example verification and experimental verification.
3.1.1.Calculation benchmark verification.Needle-plate discharge is a classical plasma discharge structure.A lot of research has been carried out by early scholars,mainly using finite volume method to solve the transport equation,Poisson equation and photoionization term.Among them,Kulikovsky and others from Russian Academy of Sciences carried out numerical calculation research on needle-plate discharge,and constructed a series of numerical simulation benchmark examples which are easy to verify.Reference[13]calculated the breakdown process of a needle-plate discharge structure with a gap of 1 cm under a DC voltage of 13 kV,considering the ionization reaction,three-body and two-body attachment reactions,positive ion electron recombination and positive and negative ion recombination reaction,and successfully reproducing the avalanche-streamer transition phenomenon and the streamer propagation process.The parameters and equations in[13]are very complete and detailed,and have been reproduced by many research groups.It is widely recognized as one of the important benchmark examples for numerical calculation of atmospheric pressure low temperature plasma.
In this section,the modified model is used,and the same governing equations,transport parameters and reaction system as in[13]are used for calculation.The results are shown in figures 5 and 6.
Figure 5.Comparison of the electron density between the modified model and Kulikovsky’s benchmark case[13].(a)The axial distribution of electron density in the benchmark case,(b)the axial distribution of electron density in the modified model,(c)distribution cloud map of electron density in the benchmark case,(d)distribution cloud map of electron density in the modified model.
Figure 6.Comparison of the electric field between the modified model and Kulikovsky’s benchmark case[13].(a)The axial distribution of electric field in the benchmark case,(b)the axial distribution of electric field in the modified model,(c)distribution cloud map of electric field in the benchmark case,(d)distribution cloud map of electric field in the modified model.
Figures 5(a)and(c)are the cloud images of the axial electron density distribution and the spatial distribution of the electron density at different times obtained in[13].Figures 5(b)and(d)are the corresponding results obtained by using the modified model.Figures 6(a)and(c)are the cloud images of axial electric field distribution and spatial distribution of electric field at different times obtained in[13].Figures 6(b)and(d)are the corresponding results obtained by using the modified model.The calculation results of the modified model are qualitatively consistent with those of the benchmark example,but there is a certain difference in quantity.In order to better quantitatively analyze the difference between the two,this paper provides the corresponding data and error percentage in table 4.
Table 4.Quantitative analysis between calculation results and benchmark example results.
It is worth noting that the streamer diameter calculated by the modified model is slightly smaller than that calculated by the reference example,and the reason for this difference may be the different calculation methods of photoionization.Kulikovsky[13]adopted the simplified ‘integral ring’method,which only calculates the photoionization source term near the streamer head with strong electric field,and simplifies the global integration to the annular volume near the streamer head.In this paper,we use the Helmholtz equation photoionization model,the accuracy of which has been widely verified[30],and use the new boundary condition ‘second order scattering boundary condition’developed in[1]to make the calculation more efficient and accurate.The red lines of the benchmark example in figures 4(a)and(c)corresponds to the third lines(19 ns time)on the left in figures 5(a)and(b).Comparing with figures 4 and 5,there is no numerical oscillation and pseudo diffusion in the calculation results of the modified model.
3.1.2.Verification of experimental results.Although the results of[13]can be used as a good benchmark for numerical calculation,they have not been fully verified by experiments.Based on overcoming the numerical oscillation and pseudo diffusion problems caused by the standard finite element method,the modified model is compared with the existing experimental results.
Intensified Charge Coupled Device is a common tool in streamer discharge experiments.Reference[31]used it to study the optical properties of needle plate structure with 1.6 cm gap under the action of nanosecond pulse overvoltage(43 kV)at atmospheric pressure,as shown in figure 7(a).The follow-up study[32]also showed that when the voltage rise time is significantly less than the characteristic time required to produce the electric field shielding effect,the local electric field is significantly higher than the ionization threshold,so the large diameter inverted cone discharge structure is produced.
Under the same conditions as[31],this paper attempts to reproduce the inverted cone discharge structure by using the uncorrected plasma module.However,due to the‘black box’,the model does not converge and the calculation cannot be completed.Therefore,this paper only shows the calculation results using the modified COMSOL Multiphysics plasma module in figure 7(b).
Figure 7.Comparison between overvoltage discharge experiment and the modified model.(a)Experimental measurement of discharge morphology,(b)the radiation intensity of the second positive band system calculated by modified model(Abel transformation).
In the experiment,the main source of optical radiation is N2(C3Пu),the second positive band system of nitrogen molecule.In order to directly compare with the experimental results,the density of N2(C3Пu)is calculated by using the modified model,and the calculated results are transformed by Abel transformation to consider the thickness of the dischargeregion.The calculated results are normalized between 0 and 1.As shown in figure 7,both experiments and calculations capture the local strong radiation in the anode and cathode regions,and the maximum discharge diameter near 1 cm ofzaxis is 1.2 cm.The calculated results of the modified model are still in good agreement with the experimental results.
Both numerical and experimental results show that the modified model can effectively overcome the numerical oscillation and streamer pseudo diffusion problems in discharge calculation,ensure the accuracy and repeatability of plasma calculation,and improve the reliability of air streamer discharge calculation results at low temperature and sub atmospheric pressure.
In this paper,the needle-plate streamer discharge at low temperature and sub atmospheric pressure is studied by using the modified model,mainly focusing on the electron density,electric field,and propagation process.However,due to the detailed introduction of needle-plate streamer discharge at room temperature and atmospheric pressure,the calculation results of needle-plate streamer discharge at low temperature and sub atmospheric pressure will not be repeated.
In order to compare with the needle-plate electrode,this section still analyzes the discharge characteristics of the streamer under the needle-needle electrode structure from the three aspects of electron density,electric field,and streamer propagation speed.
The first is the electron density.It can be seen from figure 8 that the electron density of the positive streamer head is higher than that of the negative streamer head.The main reason is that the negative streamer develops in the same direction as the electron migration,while the positive streamer develops in the opposite direction.Therefore,the electron density gradient of the positive streamer head is larger,and the electric field generated by space charge is higher,which leads to the ionization reaction and photoionization reaction of the positive streamer head more intense than that of the negative streamer.In order to further analyze the streamer propagation speed,this paper selects three times with the same time interval,which aret=11 ns,t=12 ns,t=13 ns,and gives the annotation in figure 8.11–12 ns is called the first-time interval,and 12–13 ns is called the second-time interval.The analysis shows that the streamer propagation distance of the first-time interval is significantly less than that of the second-time interval,that is,the streamer propagation speed of the second-time interval is higher than that of the first-time interval.The reason of different streamer propagation speeds in the same time interval will be analyzed in section 3.3.
Figure 8.Axial profiles of the electron density.
The second is the electric field.According to figure 9,in the initial stage of streamer,the electric field at the streamer head is mainly affected by the high background field,and the space charge almost does not exist,showing Laplacian field distribution.When the streamer enters the propagation stage,the streamer head gradually moves away from the high electric field region near the tip of the needle.Due to the existence of space charge,the streamer head presents a Possion field distribution.According to the above,due to the large electron density gradient in the positive streamer head,the electric field in the space charge layer is high,so the peak electric field always stays in the positive streamer head region.The anode potential is positive relative to the plasma potential,so it attracts electrons near the anode tip,repels positive ions,makes the electron density exceed the ion density,and finally forms an electron space charge layer near the anode tip,which is called ‘electron sheath’.Similarly,because the cathode potential is negative relative to the plasma potential,the cathode tip attracts positive ions and repels electrons,so that the ion density exceeds the electron density,and finally forms an ion space charge layer near the cathode tip,which is called ‘ion sheath’.According to the electric field distribution,the thickness of the ion sheath is obviously larger than that of the electron sheath because the electron mobility is much higher than that of the ion mobility.Therefore,the electric field in the ion sheath is higher than that in the electron sheath,which is consistent with the distribution of the electric field in the sheath at room temperature and atmospheric pressure.Three times with the same time interval are also marked in figure 9,which are 12 ns,13 ns and 14 ns respectively.In addition to the difference of propagation speed,the electric field of streamer head decreases in both time intervals.The reason for the decrease of electric field in streamer head will also be analyzed in section 3.3.
Figure 9.Axial profiles of the electric field.
The final one is the speed of streamer propagation.In order to better explain the streamer propagation time,the neutral particle density N at low temperature and sub atmospheric pressure is calculated according to formula(17).
wherekBis the Boltzmann constant,kB=1.380649×10−23(J/K).According to the calculation,N=7.23×1018cm−3,which is 1/3 of the neutral particle density at normal temperature and pressure.The decrease of neutral particle density means that there is a larger reduced electric field at low temperature and sub atmospheric pressure under the same gap and applied voltage.According to[33],the breakdown threshold of reduced electric field in air at atmospheric pressure is 120 Td.Therefore,it can reach the threshold faster and enter the initial stage of streamer earlier at low temperature and sub atmospheric pressure.In the streamer propagation stage,the increase of reduced field makes the ionization reaction at the streamer head more intense and improves the streamer propagation speed.The propagation time at low temperature and sub atmospheric pressure is smaller than that at normal temperature and pressure,that is to say,the propagation speed at low temperature and sub atmospheric pressure is greater than that at normal temperature and pressure,which is consistent with that in[1].Figures 10 and 11 show the spatial electron density distribution and electric field distribution att=1 ns,10 ns and 15 ns,respectively.The whole streamer discharge process can be better understood by the spatial distribution at different times.
Although the electric field distribution of the two electrode structures belongs to non-uniform electric field,there are still differences.It can be seen from figure 10 that in the needleneedle electrode structure,the negative streamer discharge near the cathode tip due to high electric field produces a certain range of pre-ionization region,and the streamer propagation speed and peak electric field change greatly before and after the positive streamer enters the region.The preionization region near the cathode is the main reason for the difference between the two electrode structures.
Figure 10.Space electron density distribution at different times.
In order to study the discharge process more conveniently,firstly,according to the three different times shown in figures 10 and 11,the discharge process under the needleneedle electrode is divided into three different stages.(a)The time range of the stage is 0–2 ns,which belongs to the initial stage of streamer,and the positive streamer is gradually generated near the needle electrode.(b)The time range of the stage is 2–12.5 ns,which belongs to the streamer propagation stage,and the positive and negative streamers begin to propagate,and gradually approach.(c)The time range of the stage is 12.5–15 ns,which belongs to the positive and negative streamer connection stage,and the positive streamer gradually enters the negative streamer pre-ionization region,while the plasma parameters are greatly different before and after the connection.
Secondly,according to figure 11,in the whole process of propagation,the head of positive streamer is always the area where the peak electric field is located,so this paper uses the position change of the peak electric field to represent the process of positive streamer propagation.It can be seen from figure 12 that in stage(a),the streamer propagation speed is relatively fast,reaching 700 km s−1due to the high electric field near the tip;when the positive streamer enters stage(b),the streamer propagation also enters a stable state,with the propagation speed of about 520 km s−1,because the electric field decreases and tends to be stable after it is far away from the high electric field area of the tip;when the positive streamer enters stage(c),it is affected by the negative streamer,and the positive streamer accelerates to the cathode in the pre-ionization area generated by the negative streamer,and the speed reaches nearly 1280 km s−1,which is 2.5 times of that in stage(b).
Figure 11.Space electric field distribution at different times.
Figure 12.Position of the maximum electric field Emax along the axis of symmetry as a function of time.
As shown in figure 13,with the growth of time,the peak value of electric field also presents different changing rules in three stages in addition to changes in its position.In stage(a),due to the application of external voltage,the electric field increases sharply,reaching the maximum value of 160 kV cm−1in the whole discharge process att=1 ns;when the positive streamer starts to propagate into stage(b),the peak value of electric field decreases gradually due to the gradual departure from the high electric field region near the tip,and almost tends to the fixed value of 90 kV cm−1at 6–10 ns.At the end of stage(b),the peak value of electric field increases slightly due to the influence of negative streamer,and it is worth noting that the peak value of electric field fluctuates slightly because of space charge;when the positive streamer enters stage(c),the discharge is easier due to the pre-ionization region generated by negative streamer,and the required electric field decreases,and the peak value of electric field decreases greatly.
Figure 13.Position of the value of Emax along the axis of symmetry as a function of time.
It can be seen from the above discussion that stage(c)is a unique discharge process in the needle-needle electrode structure,which is different from the needle-plate electrode structure due to the existence of pre-ionization region caused by negative streamer.The level of pre-ionization determines the streamer discharge characteristics in stage(c).In order to better discuss the effect of pre-ionization on stage(c),background ionization is temporarily used instead of photoionization in the model.Although the background ionization is for the whole computational domain,this paper only focuses on stage(c).
Firstly,pre-ionization is divided into three levels:ne=np=1011cm−3,ne=np=1012cm−3,ne=np=1013cm−3,nn=0 cm−3.The variations of streamer propagation speed and peak electric field at three different levels are discussed.As shown in figure 14,with the increase of preionization level,the speed change of positive streamer before and after entering pre-ionization region tends to be smooth.When the pre-ionization level reachesne=np=1013cm−3,the positive streamer propagation speed almost keeps the same throughout the discharge process.
Figure 14.Time evolutions of the position of the maximum electric field Emax along the axis of symmetry for three different preionization levels: ne=np=1011 cm−3,ne=np=1012 cm−3,ne=np=1013 cm−3,nn=0 cm−3.
It can be seen from figure 15 that regardless of the preionization level,the variation trend of the peak electric field in the whole discharge process is almost the same.It is worth noting that the peak value of electric field tends to a fixed value in the period of 6–10 ns in stage(b),which gradually decreases with the increase of pre-ionization level,and is 90 kV cm−1,80 kV cm−1,and 70 kV cm−1respectively.There is no need for high electric field to maintain streamer propagation when the pre-ionization level is increased.
Figure 15.Time evolutions of the value of Emax along the axis of symmetry for three different pre-ionization levels: ne=np=1011 cm−3,ne=np=1012 cm−3,ne=np=1013 cm−3,nn=0 cm−3.
In this paper,based on the COMSOL Multiphysics plasma module with artificial stability term,and applying the boundary conditions developed in the previous paper to the calculation of photoionization rate,an improved air discharge model under non-uniform electric field is constructed.The model is complementary to the previous paper and forms a complete set of numerical simulation system of air discharge at low temperature and sub atmospheric pressure.The main conclusions are as follows:
(1)In this paper,based on the streamline upwind/Petrov Galerkin method,an artificial stability term is added to the COMSOL Multiphysics plasma module,and an improved modified model is constructed.The calculation results of the modified model can keep good consistency in both the example verification and the experimental verification,which proves that the model can effectively overcome the numerical oscillation and streamer pseudo diffusion problems caused by the algorithm defects of COMSOL Multiphysics built-in module,and ensure the accuracy and repeatability of the calculation.
(2)The modified model was used to simulate the streamer discharge process under the needle-needle electrode,and the one-dimensional axial distribution and two-dimensional spatial distribution of the electron density and electric field were obtained.Because the density of neutral particles under low temperature and sub-atmospheric pressure is only 1/3 of that under normal temperature and pressure,it has a greater reduced electric field under the same voltage and spacing.The increase of reduced electric field shortens the time of streamer through the gap.Therefore,the streamer propagation speed at low temperature and sub atmospheric pressure is faster than that at normal temperature and pressure.
(3)The process of air discharge is divided into three stages:the initial stage,the propagation stage,and the connection stage.It can be seen from the comparison that the biggest difference between the two electrode structures lies in the unique positive and negative streamer connection stage of the needle-needle electrode structure.In the connection stage,due to the existence of pre-ionization region caused by negative streamer,the propagation speed of positive streamer in the pre-ionization region is doubled after the positive and negative streamers are connected.However,with the increase of pre-ionization level,the propagation speed changes before and after the connection tend to be smooth.When the pre-ionization level reachesne=np=1013cm−3,the speed of positive streamer almost remains unchanged during the whole discharge process.In the connection stage,the peak value of the electric field of the needle-needle electrode is lower than that of the needle plate electrode.The main reason is that the pre-ionization leads to the decrease of the required electric field.The higher the pre-ionization level is,the lower the required electric field is.However,regardless of the pre-ionization level,the peak value of electric field tends to a fixed value at 6–10 ns in the propagation stage.The constant value decreases with the increase of preionization level.
Due to the limitation of paper length,this paper only discusses the electron density,electric field,and streamer propagation speed in plasma parameters.The study of electron temperature,net charge density,streamer radius and other parameters will be reflected in future papers.
First,we would like to thank the National Key RESEARCH and Development Program of the Ministry of Science and Technology ‘Life Prediction and Operation Risk Assessment of UHV Equipment under long-term Service conditions(No.2017YFB0902705)’for supporting this work.Second,thanks to the No.703 Research Institute of CSIC(China Shipbuilding Industry Corporation)and Yunnan Electric Test&Research Institute Group CO.,Ltd for assistance in this paper.The authors would also like to thank Dr Yifei ZHU for helpful discussions.
Plasma Science and Technology2021年7期