• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    Thermophysical properties of pure gases and mixtures at temperatures of 300–30 000K and atmospheric pressure:thermodynamic properties and solution of equilibrium compositions

    2022-01-10 14:51:44ZhongyuanCHI池中源WeijunZHANG張衛(wèi)軍andQiangdaYANG楊強大
    Plasma Science and Technology 2021年12期

    Zhongyuan CHI (池中源) , Weijun ZHANG (張衛(wèi)軍) and Qiangda YANG (楊強大)

    School of Metallurgy, Northeastern University, Shenyang 110819, People’s Republic of China

    Abstract The equilibrium compositions and thermodynamic properties (density, enthalpy, etc at constant pressure)of plasma of pure gases and mixtures under local thermodynamic nonequilibrium have been calculated in this paper.The homotopy Levenberg–Marquardt algorithm was proposed to accurately solve nonlinear equations with singular Jacobian matrices, and is constructed by the Saha equation and Guldberg–Waage equation combined with mass conservation, the electric neutrality principle and Dalton’s partial pressure law,to solve the problem of dependence on the initial value in the process of iteration calculation.In this research, the equations at a higher temperature were solved and used as the auxiliary equations, and the homotopy control parameters’ sequence of the homotopy equations was selected by equal ratios.For auxiliary equations, the iterative initial value was obtained by assuming that there were only the highestvalence atomic cations and electrons at this temperature,and the plasma equilibrium composition distribution with the required accuracy was ultimately solved under the current conditions employing the Levenberg–Marquardt algorithm.The control parameter sequence was arranged according to the geometric sequence and the homotopy step was gradually shortened to ensure continuity of the homotopy process.Finally, the equilibrium composition and thermodynamic properties of pure N2,Mg(30%)–CO2(70%)and Mg(40%)–CO(50%)–N2(10%)mixture plasma at atmospheric pressure were calculated and the calculation process of some specified temperatures was shown and analyzed.The calculation accuracy of equilibrium composition is higher than other findings in the literature.The results for the thermodynamic properties are in good agreement with data reported by the literature.

    Keywords:plasma equilibrium compositions,thermodynamic properties,homotopy Levenberg–Marquardt algorithm (HLMA), auxiliary equations

    1.Introduction

    Plasma, as a new clean and important energy source, has the advantages of high temperature, high enthalpy, high energy concentration, and controllable compositions.Hence, it has been widely used in metallurgy [1], electric-arc furnaces [2],semiconductor processing [3], material surface modification[4],environmental protection[5],and other fields.Plasma is a mixture of neutral molecules and/or atoms produced by a series of chemical reactions.It is composed of electrons,neutral molecules, atoms, and their cations with high conductivity, heat conduction, and fluidity.The thermodynamic properties are the necessary characteristic parameters for analyzing the temperature distribution in a plasma region.For further studies on plasma application, the exact composition distribution in thermal plasma is the basic condition for calculating or obtaining thermodynamic properties.However, it is very difficult to directly solve a nonlinear problem such as the plasma equilibrium composition model iteratively,and the difficulty lies in the determination of the initial iteration values.In this paper, a new homotopy Levenberg–Marquardt algorithm (HLMA) is proposed to solve the nonlinear problem based on the particularity of higher temperature and the composition distribution of some plasma is calculated;the corresponding thermodynamic parameters are obtained at the same time.

    For the calculation of plasma thermomechanical parameters, the handing method is simple and easy if the particle number density and internal partition function are given[6, 7], which will not be described in too much detail here.The models expressing the plasma composition have two forms.The first one is composed of the mass conservation equation,Dalton’s law, and the charge conservation equation combined with the Saha equation and the Guldberg–Waage equation.The other one is obtained by solving the minimum Gibbs free energy of the plasma system.These two forms are equivalent in gas[8,9].The first method directly presents the physical and chemical process of plasma formation, the Saha and Guldberg–Waage equations are in the form of the equations, and abstract mathematical model is intuitive and simple.

    For modeling the plasma composition, some researchers have done a lot of studies.Cressault [10] calculated the plasma composition and the net emission coefficient (NEC)for Ar–H2–He thermal plasma in local thermodynamic equilibrium at atmospheric pressure.Murphy and Arundelli [11]determined the equilibrium composition of thermal plasma using the principle of minimization of the Gibbs free energy.Wang [12] determined the equilibrium composition and thermodynamic properties of carbon and water in the temperature range of 400–30 000 K and the pressure range of 0.1–10 atm by the method of Gibbs free energy minimization.Gleizes and Chervy[13]calculated the composition of a twotemperature SF6plasma in the pressure range of 0.1–1.6 MPa using the generalized Saha and Guldberg–Waage equations.Colombo and Ghedini[14]obtained the chemical equilibrium composition and determined the thermodynamic and transport properties in nonequilibrium argon, oxygen, and nitrogen thermal plasmas.Gordana [15] calculated the equilibrium plasma composition of a U-shaped argon-stabilized DC arc based on minimization of the Gibbs free energy.

    In both forms, they eventually converge to a system of nonlinear equations; Gordon and Mcbride [16] minimization is the most popular method to solve the plasma composition model, and Lagrange multipliers along with steepest descent Newton–Raphson iteration are utilized to solve the nonlinear equation system numerically [12–14].Godinet al[17] proposed a methodology that uses the chemical basis concept to reduce to a minimum the order of the system of equations to solve.Frolov and Ivanov [18] solved these equations by the modified Newton method for increments of logarithms of unknown values.Maria[19]and Zhonget al[20]constructed a new distribution function according to the maximum entropy principle and the statistical theory methods for gases with physical and chemical processes, to avoid solving the Saha equations.However, these methods have an obvious disadvantage in that they are sensitive to the iteration initialization and the globally optimal solution cannot be obtained if no reasonable initial values are given.

    In order to find the global optimal solution,the HLMA is proposed, and this method combines the advantages of the Levenberg–Marquardt algorithm (LMA) and homotopy analysis method(HAM).The LMA is an iterative method located between the Gauss–Newton method and the steepest descent method, and can be used to deal with singular nonlinear equations [21]; meanwhile, the adaptive step size also improves the chance of escaping local optimality during iteration.The HAM is a global convergence algorithm and is not constrained by the initial values[22].These methods play important roles in solving nonlinear equations.Fan [23]introduced the linear search for the approximate Levenberg–Marquardt (L-M) step, improved the convergence speed, and reduced the calculation load.Chen [24] established a highorder accelerating modified L-M method by adjusting the approximate L-M step one more time at every iteration.The convergence order of this new method was shown to be fourth.However, the LMA has the problem that the initial values are difficult to determine like other iterative methods.The HAM has been developed for the solution of nonlinear ordinary differential equations with a fractional order [25].Abbasbandy [26] applied the HAM to solve the nonlinear heat transfer problem by adjusting and controlling the convergence region.However,the HAM has insufficient ability to deal with singular nonlinear equations.When the value oftkapproaches 1, the homotopy equations are singular, and how to adjust solutions to ensure the continuity of the iterative process is the key issue.

    In this paper, the HLMA is proposed to solve singular nonlinear equations of modeling the plasma compositions.In this new algorithm, the advantage of the HAM with no constraint for the initial value is combined with the advantage of rapid convergence of the LMA realizing the solution of singular nonlinear equations.In other words, the HAM provides an exact iteration initial value for the LMA, and the LMA quickly solves the homotopy equations.For the HLMA,two important tasks must be taken into account, i.e., the auxiliary equations and the homotopy control parameters’sequences [tk].The equations expressing the plasma composition at a higher temperature which could be solved accurately were chosen as the auxiliary equations,and the control parameter sequences [tk] were constructed using the equilateral ratios.To demonstrate the accurate calculation of the plasma compositions by the HLMA, three example cases are described and discussed in detail in the third part.

    Four parts are covered in this paper.Firstly, the mechanism of the plasma composition model under local thermodynamic nonequilibrium and thermal nonequilibrium parameter θ=Te/This elaborated based on the conservation equation of mass, the conservation equation of charge, Dalton’s law, and the Shah and Guldberg–Waage equations.Based on the theory of statistical thermodynamics, the calculation model of thermodynamic parameters is established.Secondly, the solution procedure of the plasma composition model is deduced by the HLMA to solve the auxiliary equations and set up the control parameter sequences [tk].Thirdly, in order to verify the general applicability of the algorithm, examples of different atomic species and different auxiliary temperatures are given to solve the problem,and the example calculations are shown to demonstrate the feasibility and accuracy of the HLMA method.Finally, the equilibrium compositions of pure N2, Mg–CO–N2and Mg-CO2mixture plasmas in specific proportions are calculated, and then the corresponding thermodynamic parameters are obtained.

    2.Model of equilibrium compositions and thermodynamic properties

    2.1.Equilibrium compositions

    Plasma is a mixture of neutral molecules formed after a series of chemical reactions.The total pressure of the mixture is the sum of the gas pressure of each component,

    wherepis the total pressure of the plasma;neandniare the number densities of electrons and theith particle, respectively;kBis the Boltzmann constant; andTeandThare the electron and heavy temperatures of the plasma, respectively.θ=Te/This the thermal nonequilibrium parameter.

    The ratio of the number densities of two different elements is constant:

    wherenziis the number density of theith atom element andnz1is the number density of the first atom element.

    The plasma is neutral on the macro level:

    whereziis the electric charge of theith particle including electrons.

    The ionization reactionar→a(r+1)++e and dissociation reactionab→a+bare the main chemical reactions in plasma.The chemical equation and the chemical equilibrium relationship among the concentrations of different molecules can be expressed by the Saha equation equation (4) and Guldberg–Waage equation (5),

    where e is an electron,meis the mass of an electron,andhis the Planck constant.AB represents a polyatomic molecule, and A and B are the products of the decomposition AB.nAB,nA,nB,mAB,mA, andmBare the number density and mass of AB, A,and B, respectively.QrandQr+1are the partition functions of the atoms and their ions,respectively,QAB,QA,andQBare the partition functions of the AB,A,and B molecules,respectively,andEI,r+1andEdare the ionization and dissociation energies of the associated chemical reaction,respectively.The calculation of the partition coefficient can be referred to in the relevant literature [14, 27, 28], but will not be introduced here.

    The model can be expressed by equations (1)–(5).They include one charge conservation equation, one pressure balance equation,Ne?1 mass conservation equations(Neis the number of atomic types and these equations do not exist ifNe=1), andm?Ne?1 chemical equations.

    The values of the right-hand side of equations (4) or (5)are significantly different between the different chemical reactions.For example, in the air plasma, whenT=3000 K,their maximum value divided by their minimum value may equal nearly 108which is the main reason for the singularity of the equation system.Moreover, the coefficients of other equations are less than or around 10 while the coefficients of the chemical reactions are in the range of 10?5–1010.The Saha and Guldberg–Waage equations can be expressed by equation (7) in terms of vectors so that equation (6) can be expressed as equation (8):

    2.2.Partition function

    The atom partition function is the sum of electron level states,diatomic molecules are the result of interactions between electronic, vibrational, and rotational energy levels, and a polyatomic molecule is the product of simple harmonic partition functions ignoring the contribution of rotational energy levels.

    For atoms and their ions,

    wherewiand εiare the degeneracy (statistical weight)and energy level of theith energy level of the atom, respectively.Atomic energy levels are based on the latest spectral data from the National Institute of Standards and Technology(NIST) [29].

    For diatomic molecules and their ions,

    where εe,εv,εr,ge,gv,andgrare the electron vibrational,rotational energy level, and the corresponding degeneracy of the molecule, respectively.

    The vibration energy of thevth vibrational state andeth electronic state can be computed as

    where ωe, ωexe, ωeye,Be,De, αe, and βeare the spectroscopic constants of theeth electronic state of a diatomic molecule.BecauseDvis much smaller thanBv,theDvterm is ignored when calculating the rotational energy levels in order to simplify the calculation.

    For polyatomic molecules,

    For linear molecules, vibrational degrees of freedomm=3n?5, nonlinear moleculem=3n?6, andnis the number of atoms in the molecule.

    2.3.Thermodynamic properties

    Once compositions and partition functions are known, it is possible to compute the thermodynamic properties.According to statistical thermodynamics, and given the number densities, the average density equation (14), internal energy equation (15), specific enthalpy equation (16), and constantpressure specific heat equation (17) have been listed.

    where ρ,e,h,andcpare the density,internal energy,enthalpy,and constant-pressure specific heat of plasma, respective.n,m, andEare the number density, mass, and formation enthalpy of a molecule, respectively, andQintrepresents the partition functions.It can be seen from equations (14)–(17)that the thermodynamic parameters are the summation of the data of different particles and can be calculated directly.

    3.Homotopy Levenberg–Marquardt algorithm to solve plasma equilibrium compositions

    3.1.HLMA

    The HAM is based on both the homotopy in the topology and the Maclaurin series, and is a simple analytical technique for linear and nonlinear problems[24,25].This method does not require small or large parameters in the perturbation techniques, but it cannot be applied to singular nonlinear equation systems including the current equation system.In this paper,the LMA was combined with the HAM to solve singular nonlinear equations expressing the plasma composition model.The LMA was responsible for solving the homotopy equations, while the HAM was used to construct the homotopy equations and provide the iterative initial values for the LMA.The homotopy equations are expressed as

    where H(x,t) represents the homotopy equations,tis the homotopy sequence factor, F(x) represents the solving equations, and G(x) represents the auxiliary equations that need to be constructed or assigned.Ift=1, H(x,t)=G(x),and ift=0, H(x,t)=F(x).Essentially, the process of the HAM method in dealing with equations is to conduct the transformation from H(x,1)=G(x) to H(x,0)=F(x) by solving a series of homotopy equations H(x,tk),wheretkis thekth homotopy sequence factor and 0

    3.2.Solution of auxiliary equations and construction of homotopy equations

    For the plasma composition model, there is a common sense conclusion that only the highest-valence atomic cations and electron are present in plasma at a higher temperature such as 30 000 K.The common sense conclusion can also be used as the initial value of the iterative methods and the solution satisfying requirements can be obtained through the iterative process.Hence, equations at a higher temperature were chosen as the auxiliary equations.Such equations had the advantages of a highly precise solution and the same structure as solving the system at any temperature.If the solution at one temperature can be obtained, we can have solutions at any temperature below that temperature via the HLMA.The homotopy equations can be expressed as

    where FTmax(x) represents the plasma composition equations at a higher temperature.The data array corresponding to the number densities of electron and the highest-valence atomic cations is existing or nonzero and all those corresponding to the number densities of other particles are equal to zero or do not exist.This data array serves as a reasonable initial value for iteration and calculation of the exact solution of the auxiliary equations by the LMA.Otherwise, we can only reconstruct the auxiliary equations by raising the temperature ofTmaxsuch asTmax=Tmax+2000.Every homotopy equation is still in a singular nonlinear system that cannot be solved analytically.However,the LMA can easily solve these equations if a reasonable initial value is given.When the ratio of HT(x,tk?1)/HT(x,tk) is close to 1, the solution of HT(x,tk?1)can be used as the initial value of HT(x,tk).In this way,we achieved the combination of the homotopy method and LMA where the homotopy method determined the homotopy equations by adjustingtkand the LMA solved these equations.

    3.3.Determination of the homotopy sequence

    Another key factor of the HLMA is to determine the control parameter sequence [tk] such as 0=t0

    The last solution can be obtained byNtime iterations of the LMA calculation,whereNis the number of the homotopy step.The equations of HT(x,tk) are also singular nonlinear equations.For a given plasma, the maximum value occurs in a fixed chemical reaction that does not vary with temperature and is marked asSa,max.Sa,maxis the maximum singular value and this reaction is the maximum singular traction.Similarly,the maximum singular valueSt,maxof the solving equation can be obtained from the same reaction.In this paper, the ratioRmaxof the maximum singular coefficientSa,maxof FTmax(x)=0 toSt,maxof FT(x)=0 was selected as the reference value to determine the control parameter sequence[tk].In regard to the homotopy step number,the larger number is the better in principle.Hence,we selectedN=(Tmax?T)/50 in the initial period in this paper.When using the LMA algorithm to solve HT(x,tk) =0, it is necessary to calculate the Jacobian matrix J.J is calculated by

    The element is calculated by

    In the actual calculation,ifxjis not equal to 0,dx=1e–5 ×xj; else dx=1e – 5.

    The control parameter sequence[tk]is constructed via the equal-scale model, and the proportionality coefficientqand the homotopy coefficienttkare expressed as

    wheretk=1 whenk=0 andtk=0 whenk=N.Whenkis small, the difference oftk?tk?1is large, and with the increase ink, this difference gradually decreases.The step numberNis the key parameter to adjust the homotopy sequence so that each step of LMA calculation can be completed successfully.If LMA fails,Nmust be dynamically adjusted to reduce the difference oftk?tk?1to make sure that the solution of HT(x,tk?1) will be a reasonable initial value for HT(x,tk).This process is iterated until the whole homotopy process is accomplished.In this paper,Nais defined as the number of successful solutions by the LMA.When theNa+1 LMA calculation fails, ifNa<10,Nis reset toN+Na; else,Nis reset toN—Na/2.The proportionality coefficientqis calculated as follows:

    The values ofqandtkare recalculated by equations(23)and (24).The purpose of adjusting step N is to reduce the number of LMA calculations and avoid wasting computing resources..

    3.4.Procedure of the HLMA

    Figure 1.Calculation flow chart of the HLMA for solving the plasma composition.

    The procedure of the HLMA can be divided into two parts that are executed successively.The first part is to construct and solve the auxiliary equations,and the second part is to adjust the step numberNand achieve the solution of a series of homotopy equations by the LMA; the conditions of convergence of the LMA are||F||

    The calculation process of the HLMA is shown in figure 1 and the detailed procedures for solving the plasma composition model are as follows:

    (1) Determine the maximum temperatureTmaxand homotopy step numberNmaxand construct the auxiliary equations G(x);set the calculation accuracyepsand setNa=0,tc=0 at the same time.

    (2) Assuming that only electrons and the highest-valence atomic cations exist in the plasma, construct and solve the linear equations containingNe+ 1 variables according to the law of charge conservation, Dalton’s law of partial pressure,and the assumption that the ratio of the total number densities of atoms is constant; then the exact solution or a reasonable iterative initial value under the assumption is obtained.

    (3) Solve the auxiliary equations and obtain their variancep.Ifp

    (4) Find the maximum singular valuesSa,maxof the auxiliary equations andSt,maxof the solving equations,and calculateRmax=Sa,max/St,max.

    (5) Calculate the proportional coefficientqand setk=0.

    (6) Compute the homotopy coefficienttkand construct the homotopy equations H(x,tk).

    (7) Solve the homotopy equations H(x,tk) by the LMA,with two conditions of convergence: (a) ||F||

    The physician s failure to recognize that he is subordinate to his godfather, and that despite his privileged position as godson he too is mortal, leads to his downfall. Note that he is not satisfied with wealth and power, but now strives to marry a princess and win a crown. What is considered a legitimate30 goal in tales of magic becomes a mark of hubris31 in this tale. (198)

    (8) Calculate their variancep.If notp

    (9) IfNa<10, setNmax=Nmax+Na.Otherwise, setNmax=Nmax—Na/2,Rmax=(1 ?tc)Rmaxand go to step 5.

    (10) If notk=Nmax, setk=k+1,tc=tkand go back to step 6.Otherwise, go to step 11.

    (11) Finished.

    3.5.Cases

    According to the above statement, the HLMA can be used to solve the plasma equilibrium composition at arbitrary temperatureT,number of atomic elementsNe,and total pressureP.Three example cases are selected whereT=6000 K,P=0.1 MPa,Ne=2, and auxiliary temperatureTmax=7000 K in Case 1;T=10 000 K,P=0.1 MPa,Tmax=30 000 K, andNe=3 in Case 2; andT=12 000 K,P=0.1 MPa,Tmax=30 000 K, andNe=4 in Case 3.The computational accuracy ||F||<1e?10.They are used to demonstrate the accuracy and feasibility of HLMA for the plasma equilibrium composition model.

    Case 1: Mg(30%)–CO2(70%) mixture plasma, temperatureTe=10 000 K, and auxiliary temperatureTmax=30 000 K.There are 19 particles and 15 chemical reactions considered for this plasma, including eight neutral particles,ten ions and five dissociation reactions, and ten ionization reactions.These particles contain one triatomic molecule and four diatomic molecules,and the other particles are electrons,monatomic molecules and their ions.Table 1 shows the list of molecules in the plasma.The possible chemical reactions are listed in table 2,and table 3 lists the solutions for the solving equations (t=1) and the auxiliary equations (t=0).It can be seen that the values of some solutions vary greatly,such as CO; the value of the solution changes from 1.43e?35to 8.79e6, with the difference reaching 41 orders of magnitude.A total of 710 homotopy calculations were carried out.

    Table 1.List of molecules in the Mg–CO2 mixture plasma.

    Table 2.List of possible chemical reactions in the Mg(30%)–CO2(70%) mixture plasma.

    Table 3.Solutions of auxiliary and solving equations for the Mg(30%)–CO2(70%) mixture plasma.

    Table 4.List of molecules in the Mg(40%)–CO(50%)–N2(10%) mixture plasma.

    The variations of the control parameter, the precision,and the number density of each particle with the homotopy step number are shown in figure 2.It is found that the transition from the auxiliary solution to the exact solution is a continuous and smooth gradual process.The fluctuation of the control parametertkchanges from violent to gentle and finally becomes close to zero.Hence,we can judge that the precision of each HLMA step meets the calculation requirements.

    Figure 2.Variations of the control parameter,the precision,and the number density of each particle with the homotopy step number for the Mg(30%)–CO2(70%) mixture plasma.

    Table 6 gives the solutions of the auxiliary function(t=0)and the solution function (t=1)in the calculation of Case 2.The solution accuracy of the auxiliary function ||F||=1.13e?16,and the solving function solution precision||F|| = 1.83e?18, which is better.A total of 362 HLMA calculations were performed during the calculation, and the variations of the control parameter, the precision, and the number density of each particle with each LMA are shown in figure 3.At this temperature, the number densities of the polyatomic molecules are less than 1 cm?3.It can also be seen that only electrons, atoms, and their ions exist in the plasmaand there are no polyatomic molecules when the temperature is higher, which is consistent with our assumption for determining the auxiliary equations.

    Table 5.List of possible chemical reactions in the Mg(40%)–CO(50%)–N2(10%) mixture plasma.

    Table 6.Solutions of auxiliary and solving equations for Mg(40%)–CO(50%)–N2(10%) mixture plasma.

    Table 7.List of molecules in the N2 plasma.

    Table 8.List of possible chemical reactions in the N2 plasma.

    Table 9.Solutions of the auxiliary and solving equations for the N2 plasma.

    Figure 3.Variations of the control parameter,the precision,and the number density of each particle with the homotopy step number for the Mg(40%)–CO(50%)–N2(10%) mixture plasma.

    Case 3: pure N2plasma, temperatureTe= 6000 K and auxiliary temperatureTmax= 7000 K.The auxiliary temperature of this example is 7000 K, which is not consistent with the so-called higher temperature,but it is possible to calculate the equilibrium compositions if we have number densities that are higher than the calculated temperature.The lists of molecules and reactions for pure N2are shown in tables 7 and 8, respectively.The values of equilibrium compositions atT=7000 K andT=6000 K have been listed in table 9.In the cases ofT=6000 K orT=7000 K,the reaction N?N++e is dominant in the plasma, and the numbers of N++, N+++and N2+are very few.In figure 4, the number densities of electrons and N+particles appear to coincide.

    It can be seen from the above figures that the control parameters vary from 0 to 1.With the progress of thehomotopy process, the difference between adjacent control parameters decreases gradually, and the minimum value reaches 1e?10.The change of control parameters ensures the continuity of the homotopy process.In the process of calculation, each variable is continuously transitioning from the initial value to the final solution, which ensures the correctness of the iterative direction.The most important factor is that there is no need to worry about the calculation process falling into the local optimum and being difficult to escape.In the HLMA calculation, if the local optimum is trapped, the step size of control parameters can be adjusted appropriately,and the LMA can jump out of the local optimal solution by automatically adjusting the range of the trust region.This is also the advantage of the HLMA algorithm to solve nonlinear problems, which can not only ensure the correctness of the optimization direction but also escape from the local optimal region smoothly.

    Figure 4.Variations of the control parameter,the precision,and the number density of each particle with the homotopy step number for the N2 plasma.

    In the process of batch calculation, it is not necessary to select an auxiliary function at a higher temperature, such as 30 000 K,as long as the temperature is higher than the current solution under the same pressure.In the practical application process, the author summed up an immature experience, and suggested that ifT> 15 000 K, thenTmax–T=1000 K; ifT>8000 K, thenTmax–T=500 K; ifT>5000 K, thenTmax–T=100 K;ifT>3000 K,thenTmax–T=50 K;elseTmax–T=10 K.These suggestions can effectively shorten the gap between the auxiliary equations and the solving equations, ensure that the solving variables are in a controllable range, and also greatly reduce the number of homotopy calculations, providing work efficiency.

    4.Results and discussion

    4.1.Equilibrium composition

    In this paper,the equilibrium composition number densities are calculated in the temperature range of 300–30 000 K for pure N2, Mg(30%)–CO2(70%) and Mg(40%)–CO(50%)–N2(10%)mixture plasmas at atmospheric pressure;the calculation results are shown in figures 5–7.It can be seen from the figures that the main particles of plasma at higher temperature are electrons and the highest atomic cations(e,C3+,O3+,Mg2+,N3+),and it is shown that polyatomic molecules are completely decomposed and monatomic molecules are deeply ionized.This conclusion also proves the correctness of choosing equilibrium compositions at higher temperatures as the auxiliary function.With the decrease of temperature, low-order atomic cations gradually appear until the monatomic molecules are precipitated, and the value of number density increases rapidly.Molecular ionization and dissociation reactions are dominant at lower temperatures.In the actual calculation process, in order to reduce the amount of calculation,the plasma temperature as an auxiliary function can be solved by the HLMA algorithm as long as it is higher than the solution, rather thanTmax=30 000 K as shown in Case 3.Here, the authors make an empirical conclusion, the calculation precision of the auxiliary function||F||<1e?10,which can effectively reduce the number of LMAs.The equilibrium compositions of pure N2at temperatures of 300–3000 K and pressure of 0.1 MPa are displayed in figure 8.We found that the composition proportion of N2+is obviously higher than in the case of θ=1, and the larger θ, the greater the proportion.

    Figure 5.Equilibrium compositions of Mg(30%)–CO2(70%) mixture plasma at temperatures of 300–3000 K and pressure of 0.1 MPa.(a)300–8000 K, (b) 8000–30 000 K.

    Figure 6.Equilibrium compositions of Mg(40%)–CO(50%)-N2(10%) mixture plasma at temperatures of 300–3000 K and pressure of 0.1 MPa.(a) 300–8000 K, (b) 8000–30 000 K.

    In these calculations, other producing polyatomic molecules are considered, such as O2, C2, CO2, MgO, NO, CN,and CO.At the same time, the ionization reaction of an O2molecule with a strong oxidation property is also considered.The calculation of the partition function takes into account the various electronic states of the molecule, which effectively improves the accuracy of the composition number density,and is especially important for the calculation of spectral data.The electronic state energy level data of monatomic molecules are obtained from the NIST Atomic Spectra Database[29].For CO,the electronic states X1Σ+,A1Π and B1Σ+are mentioned, and their spectral data come from Farrenq [30],Ostrowska [31] and Eidelsberg [32], respectively.For CO+,X2Σ+[33], A2Πi[33], B2Σ+[34] and C2Δ [35] electronic states are used to calculate.For O2[36]and O2+[37],we refer to X3Σ?g,B3Σ?utwo states and X2Πg, a4Πg, A2Πg, b4Σ?gfour states.C2have eight electron states are obtained from Kirby[38], and the four electronic states MgO have been found from the NIST-JANAF thermochemical tables [39].For CN,X2Σ+, A2Πi, and B2Σ+three electronic states have been taken into account; the spectral data are available from Peteret al[40].For NO, the electronic states of X2Πr, A2Σ+and D2Σ+are given by Reddy [41], and the states of B2Σrand C2Πrare referred to[41].The spectral data of the X1Σ+state for NO+have been listed by Peteret al[40].The calculations of the partition function for CO2have been discussed by Grayet al[42]; only the vibrational level of the ground state is considered in this paper, and the vibrational constants are 2349 cm?1, 1337 cm?1and 667 cm?1(2).

    Figure 7.Equilibrium compositions of pure N2 at temperatures of 300–3000 K and pressure of 0.1 MPa.(a)300–8000 K,(b)8000–30 000 K.

    4.2.Thermodynamic properties

    Thermodynamic properties (including density, enthalpy, etc)are important conditions for calculating the temperature of a system according to the law of conservation of energy, and they are also the basis for calculating the temperature field in computer aided engineering (CAE).When the equilibrium composition number density and particle partition function are known, the mass density and enthalpy can be calculated by equations (14) and (16).For pure N2, the results for mass density and enthalpy are in good agreement with data reported by Colomboet al[14].The comparisons are shown in figures 9 and 10.For the Mg–CO2and Mg–CO–N2mixture plasmas, there are no comparable calculations, and their calculations are shown in figures 11 and 12.

    With the increase of temperature, the thermal motion behavior of particles becomes more active, and the force on the boundary becomes stronger.Under the same pressure,the total number density of plasma decreases.At the same time,the energy level of the particle jumps from the ground state to the excited state, which enlarges the internal energy and partition function of the particle.As the internal energy of the molecule enlarges, the enthalpy of the plasma also increases,and the energy density of the unit volume of the plasma becomes higher.Due to the small volume of arc plasma in industrial applications, the time and energy required to raise the plasma temperature are far less than those of conventional fossil energy, which also lays the foundation for the application of high-temperature and ultra-high-temperature plasma.

    Figure 9.Comparison of mass density in N2 plasma at temperatures of 300–30 000 K and pressure of 0.1 MPa.The literature result was calculated by Colombo et al [14].

    Figure 10.Comparison of enthalpy in N2 plasma at temperatures of 300–30 000 K and pressure of 0.1 MPa.The literature result was calculated by Colombo et al [14].

    5.Conclusions

    In this paper,an equilibrium plasma composition model under local thermodynamic nonequilibrium conditions was established by the Saha and Guldberg–Waage equations combined with the mass conservation equation,the charge conservation equation, and Dalton’s law.The homotopy Levenberg–Marquardt algorithm(HLMA)was proposed for the first time and used to solve the plasma composition model.The main conclusions are as follows.

    The auxiliary equations constructed at a higher temperature could be accurately solved by the LMA and the initial value was chosen according to the assumption that the number density ratio between electrons and the highestvalence atomic cations is constant.

    The HLMA could solve the singular nonlinear equations when the auxiliary equations were obtained and the difference in values oftk?tk?1between two adjacent control parameters was small.

    HLMA could solve the plasma composition model at arbitrary temperatures and arbitrary thermal nonequilibrium parameter θ, number of atomic elements, and pressure with high accuracy and convergence of the homotopy process.The calculation accuracy of three examples is ||F||<1e?12.

    Figure 11.Mass density of Mg(30%)–CO2(70%)and Mg(40%)–CO(50%)–N2(10%)mixture plasmas at temperatures of 300–30 000 K and pressure of 0.1 MPa.

    Figure 12.Specific enthalpy of Mg(30%)–CO2(70%)and Mg(40%)–(CO50%)–N2(10%) mixture plasmas at temperatures of 300–30 000 K and pressure of 0.1 MPa.

    Acknowledgments

    This study was supported by the National Key Research and Development Program of China(No.2017YFA0700300),the Fundamental Research Funds for the Central Universities(No.N2025032),and the Liaoning Provincial Natural Science Foundation (No.2020-MS-362).

    ORCID iDs

    国产伦一二天堂av在线观看| 亚洲国产中文字幕在线视频| 桃红色精品国产亚洲av| 精品免费久久久久久久清纯| 久久精品综合一区二区三区| 久久久国产精品麻豆| 一个人免费在线观看电影 | 国内毛片毛片毛片毛片毛片| 91av网站免费观看| 欧美zozozo另类| 丝袜人妻中文字幕| 久久久精品国产亚洲av高清涩受| 久久这里只有精品19| www.精华液| 看免费av毛片| 老司机深夜福利视频在线观看| 在线观看免费视频日本深夜| 999精品在线视频| 国产三级在线视频| 黄色女人牲交| 国产成人精品久久二区二区免费| 看免费av毛片| 国产午夜精品论理片| 国产精品一区二区精品视频观看| 亚洲 欧美 日韩 在线 免费| 人人妻人人澡欧美一区二区| 两个人免费观看高清视频| 亚洲av成人av| 少妇裸体淫交视频免费看高清 | 亚洲欧洲精品一区二区精品久久久| 国产精品电影一区二区三区| 国产伦一二天堂av在线观看| 精品少妇一区二区三区视频日本电影| 又粗又爽又猛毛片免费看| 亚洲精品久久成人aⅴ小说| 亚洲熟妇熟女久久| 亚洲一码二码三码区别大吗| 国产亚洲av嫩草精品影院| 97超级碰碰碰精品色视频在线观看| 国内久久婷婷六月综合欲色啪| 成人精品一区二区免费| 国产黄片美女视频| 黄色 视频免费看| 国产真人三级小视频在线观看| 久久久水蜜桃国产精品网| 久久精品国产亚洲av高清一级| 好男人在线观看高清免费视频| 国产v大片淫在线免费观看| 国产精品一区二区免费欧美| 不卡av一区二区三区| 999精品在线视频| 97人妻精品一区二区三区麻豆| 亚洲熟妇中文字幕五十中出| 日本a在线网址| av福利片在线| www.精华液| 91九色精品人成在线观看| 日韩中文字幕欧美一区二区| 女人高潮潮喷娇喘18禁视频| 欧美黑人精品巨大| 亚洲精品美女久久av网站| 日韩欧美免费精品| 这个男人来自地球电影免费观看| 国产真实乱freesex| 精品一区二区三区四区五区乱码| 亚洲精品国产一区二区精华液| 一级毛片精品| 精品一区二区三区av网在线观看| 国产人伦9x9x在线观看| 99久久99久久久精品蜜桃| 国产视频一区二区在线看| 嫩草影视91久久| 亚洲国产看品久久| 午夜福利免费观看在线| 亚洲国产欧美一区二区综合| 99久久国产精品久久久| 国产激情久久老熟女| 在线国产一区二区在线| 女同久久另类99精品国产91| 2021天堂中文幕一二区在线观| 18禁国产床啪视频网站| 一个人免费在线观看电影 | 可以在线观看的亚洲视频| www.精华液| 国产伦在线观看视频一区| 夜夜看夜夜爽夜夜摸| 国产在线观看jvid| 天堂动漫精品| а√天堂www在线а√下载| 精品熟女少妇八av免费久了| 亚洲av电影不卡..在线观看| 亚洲中文日韩欧美视频| 午夜两性在线视频| 亚洲乱码一区二区免费版| 波多野结衣高清作品| 亚洲国产精品成人综合色| 精品无人区乱码1区二区| 两个人看的免费小视频| 国产精品国产高清国产av| 国产精品久久久久久亚洲av鲁大| 黄片小视频在线播放| 在线观看美女被高潮喷水网站 | 欧美日韩一级在线毛片| 午夜成年电影在线免费观看| 男男h啪啪无遮挡| 精品国产亚洲在线| 一a级毛片在线观看| 久久久精品欧美日韩精品| 给我免费播放毛片高清在线观看| 人妻久久中文字幕网| 亚洲国产欧美网| 精品国产超薄肉色丝袜足j| 国产高清激情床上av| a级毛片a级免费在线| 又大又爽又粗| 欧美中文日本在线观看视频| 脱女人内裤的视频| 午夜两性在线视频| 啦啦啦免费观看视频1| 九九热线精品视视频播放| 嫁个100分男人电影在线观看| 欧美成人免费av一区二区三区| 国产1区2区3区精品| 亚洲人成网站在线播放欧美日韩| 亚洲国产精品sss在线观看| 91成年电影在线观看| 黄片大片在线免费观看| 婷婷丁香在线五月| 欧美日韩国产亚洲二区| 大型黄色视频在线免费观看| 国产精品一及| 午夜精品一区二区三区免费看| 精品久久久久久久人妻蜜臀av| 身体一侧抽搐| 国产主播在线观看一区二区| 国产亚洲精品久久久久5区| 黄色片一级片一级黄色片| 久久人妻av系列| 欧美日韩精品网址| 精品久久久久久,| 色综合站精品国产| 欧美色视频一区免费| 国产aⅴ精品一区二区三区波| 色老头精品视频在线观看| 成人三级黄色视频| 国产精品精品国产色婷婷| av福利片在线| 精品一区二区三区av网在线观看| 国产精品免费视频内射| 成人国产一区最新在线观看| 深夜精品福利| 久久天躁狠狠躁夜夜2o2o| 亚洲欧美精品综合久久99| 中文资源天堂在线| 又粗又爽又猛毛片免费看| www.999成人在线观看| 嫁个100分男人电影在线观看| 欧美中文综合在线视频| 久久久久国产精品人妻aⅴ院| 亚洲美女黄片视频| 亚洲av电影在线进入| 别揉我奶头~嗯~啊~动态视频| 欧美乱色亚洲激情| 欧美日韩黄片免| 日韩欧美三级三区| 免费看十八禁软件| 国产精品久久电影中文字幕| 亚洲最大成人中文| 国产男靠女视频免费网站| 淫秽高清视频在线观看| 麻豆国产av国片精品| 亚洲人成网站高清观看| 国产精品免费视频内射| 国产高清激情床上av| 日韩av在线大香蕉| 老司机午夜福利在线观看视频| netflix在线观看网站| www.自偷自拍.com| 琪琪午夜伦伦电影理论片6080| 亚洲色图 男人天堂 中文字幕| 国产aⅴ精品一区二区三区波| av福利片在线观看| 国产一区在线观看成人免费| 国产高清有码在线观看视频 | 搞女人的毛片| av中文乱码字幕在线| 亚洲av美国av| or卡值多少钱| 午夜精品久久久久久毛片777| 免费观看精品视频网站| 亚洲 欧美 日韩 在线 免费| 国产精品一及| 桃红色精品国产亚洲av| 九色国产91popny在线| 亚洲熟女毛片儿| 丰满人妻一区二区三区视频av | 国产亚洲精品一区二区www| 午夜视频精品福利| 日韩三级视频一区二区三区| 午夜精品久久久久久毛片777| 中文字幕久久专区| 我的老师免费观看完整版| 免费在线观看完整版高清| 国产精华一区二区三区| 久久国产精品人妻蜜桃| 亚洲成人免费电影在线观看| 久久精品国产清高在天天线| 欧美成人一区二区免费高清观看 | 50天的宝宝边吃奶边哭怎么回事| 三级国产精品欧美在线观看 | 九九热线精品视视频播放| 久久久精品国产亚洲av高清涩受| 丁香六月欧美| 国产高清视频在线播放一区| 国产av不卡久久| 久9热在线精品视频| 丁香六月欧美| 亚洲全国av大片| 欧美性猛交黑人性爽| 非洲黑人性xxxx精品又粗又长| 麻豆一二三区av精品| 久久精品人妻少妇| 亚洲精品粉嫩美女一区| 一进一出抽搐动态| 天堂av国产一区二区熟女人妻 | 99国产综合亚洲精品| 国产97色在线日韩免费| 久久久久免费精品人妻一区二区| 亚洲成人精品中文字幕电影| 亚洲av第一区精品v没综合| 91麻豆av在线| 岛国在线观看网站| av欧美777| 99riav亚洲国产免费| 日本一本二区三区精品| 亚洲中文字幕一区二区三区有码在线看 | 伦理电影免费视频| 亚洲成人免费电影在线观看| www.999成人在线观看| 国产成人影院久久av| 蜜桃久久精品国产亚洲av| 国产精品免费视频内射| 精品乱码久久久久久99久播| 久久香蕉激情| 99久久99久久久精品蜜桃| 亚洲乱码一区二区免费版| av视频在线观看入口| 亚洲片人在线观看| 一边摸一边抽搐一进一小说| 最近最新中文字幕大全电影3| 我的老师免费观看完整版| 欧美日韩中文字幕国产精品一区二区三区| 亚洲人与动物交配视频| 国产成人精品久久二区二区免费| 亚洲美女黄片视频| av视频在线观看入口| 国产高清videossex| 一区二区三区激情视频| 亚洲精品美女久久久久99蜜臀| 国产av又大| 国产1区2区3区精品| 色在线成人网| 亚洲av熟女| 法律面前人人平等表现在哪些方面| 麻豆成人av在线观看| 亚洲av成人av| 久久这里只有精品19| 少妇熟女aⅴ在线视频| 他把我摸到了高潮在线观看| 精品欧美国产一区二区三| 国产99久久九九免费精品| 久久欧美精品欧美久久欧美| 性色av乱码一区二区三区2| 成人18禁在线播放| 精品一区二区三区视频在线观看免费| 国产蜜桃级精品一区二区三区| 日韩av在线大香蕉| 久久伊人香网站| 欧美三级亚洲精品| 九色成人免费人妻av| 精品一区二区三区四区五区乱码| 成人高潮视频无遮挡免费网站| 一a级毛片在线观看| 久久精品成人免费网站| 高潮久久久久久久久久久不卡| 伊人久久大香线蕉亚洲五| 黄频高清免费视频| 亚洲欧洲精品一区二区精品久久久| 久久精品91蜜桃| 麻豆国产97在线/欧美 | 欧美乱色亚洲激情| 最好的美女福利视频网| 日本 欧美在线| 亚洲av成人精品一区久久| 久久欧美精品欧美久久欧美| 成人午夜高清在线视频| 日本一本二区三区精品| 神马国产精品三级电影在线观看 | 亚洲国产精品合色在线| 国产精品一及| 国产真实乱freesex| 午夜日韩欧美国产| 五月伊人婷婷丁香| 日本一二三区视频观看| 看片在线看免费视频| 丁香欧美五月| 中国美女看黄片| 国产精品电影一区二区三区| 少妇人妻一区二区三区视频| 一本久久中文字幕| 亚洲成av人片在线播放无| 欧美黑人精品巨大| 欧美日韩一级在线毛片| 日本熟妇午夜| 男女床上黄色一级片免费看| 欧美丝袜亚洲另类 | 最近最新中文字幕大全电影3| 亚洲专区中文字幕在线| 亚洲av熟女| 人妻丰满熟妇av一区二区三区| 亚洲自拍偷在线| 午夜福利视频1000在线观看| 国产精品野战在线观看| 午夜福利欧美成人| 久久久国产成人免费| 不卡av一区二区三区| 亚洲18禁久久av| 中出人妻视频一区二区| 全区人妻精品视频| 欧美性猛交╳xxx乱大交人| 亚洲天堂国产精品一区在线| 精品人妻1区二区| 亚洲五月天丁香| 午夜日韩欧美国产| 久久精品国产清高在天天线| 老司机午夜十八禁免费视频| 老熟妇乱子伦视频在线观看| 日韩欧美三级三区| 久久久国产成人精品二区| 精品一区二区三区视频在线观看免费| 黄色视频,在线免费观看| 观看免费一级毛片| 国产在线观看jvid| 男人的好看免费观看在线视频 | 叶爱在线成人免费视频播放| 正在播放国产对白刺激| 久久久久亚洲av毛片大全| 欧美成人一区二区免费高清观看 | 欧美日韩瑟瑟在线播放| 长腿黑丝高跟| 精品久久久久久成人av| 国产爱豆传媒在线观看 | 一级黄色大片毛片| 亚洲狠狠婷婷综合久久图片| 高清在线国产一区| 国产男靠女视频免费网站| 日本一二三区视频观看| 亚洲最大成人中文| 黄色视频不卡| 国内久久婷婷六月综合欲色啪| 又粗又爽又猛毛片免费看| 一边摸一边抽搐一进一小说| 国产午夜精品论理片| 亚洲精品美女久久久久99蜜臀| 亚洲avbb在线观看| 老熟妇仑乱视频hdxx| 亚洲国产欧美一区二区综合| 看黄色毛片网站| 亚洲国产看品久久| 日本一区二区免费在线视频| 在线a可以看的网站| 嫩草影院精品99| 久久久久精品国产欧美久久久| 午夜福利在线在线| 午夜福利视频1000在线观看| 免费在线观看视频国产中文字幕亚洲| 久久精品影院6| 欧美国产日韩亚洲一区| 国产男靠女视频免费网站| 日本一二三区视频观看| www日本在线高清视频| 国产午夜精品久久久久久| 国产激情久久老熟女| 久久欧美精品欧美久久欧美| a在线观看视频网站| 不卡av一区二区三区| 国产亚洲欧美在线一区二区| 国产av一区在线观看免费| 村上凉子中文字幕在线| 99精品久久久久人妻精品| 日本成人三级电影网站| 91av网站免费观看| 亚洲 欧美一区二区三区| 在线永久观看黄色视频| www.自偷自拍.com| 桃红色精品国产亚洲av| 超碰成人久久| 国产一区二区在线av高清观看| 丰满人妻熟妇乱又伦精品不卡| 欧美日本亚洲视频在线播放| 久久久国产成人免费| 亚洲欧美精品综合一区二区三区| 亚洲人成77777在线视频| 91大片在线观看| 一级a爱片免费观看的视频| 国产成人一区二区三区免费视频网站| 国产精品亚洲一级av第二区| 搡老妇女老女人老熟妇| 看片在线看免费视频| 精品乱码久久久久久99久播| 校园春色视频在线观看| 99国产极品粉嫩在线观看| tocl精华| 精品少妇一区二区三区视频日本电影| 午夜亚洲福利在线播放| 一区二区三区国产精品乱码| 久久久久久久久久黄片| 欧美极品一区二区三区四区| 人人妻人人澡欧美一区二区| 操出白浆在线播放| 成人特级黄色片久久久久久久| 欧美zozozo另类| 中文资源天堂在线| 九色成人免费人妻av| 亚洲欧洲精品一区二区精品久久久| 美女午夜性视频免费| 老司机午夜十八禁免费视频| 亚洲成人久久性| 99riav亚洲国产免费| 又黄又爽又免费观看的视频| 天天添夜夜摸| 久久 成人 亚洲| 国产精品,欧美在线| 校园春色视频在线观看| 欧美日本亚洲视频在线播放| 无人区码免费观看不卡| 国产一级毛片七仙女欲春2| www.自偷自拍.com| 欧美成人性av电影在线观看| 国产亚洲精品综合一区在线观看 | 一卡2卡三卡四卡精品乱码亚洲| 日本五十路高清| 一本一本综合久久| 国产精品免费视频内射| 狠狠狠狠99中文字幕| 日韩欧美在线乱码| 麻豆国产97在线/欧美 | 亚洲全国av大片| 欧美成狂野欧美在线观看| 免费无遮挡裸体视频| 久久久久久国产a免费观看| 欧美一级a爱片免费观看看 | 久久久久精品国产欧美久久久| 日韩精品免费视频一区二区三区| 亚洲一区高清亚洲精品| 又粗又爽又猛毛片免费看| 久久香蕉激情| 精品一区二区三区视频在线观看免费| 欧美成人一区二区免费高清观看 | 亚洲中文字幕一区二区三区有码在线看 | 99re在线观看精品视频| 国产区一区二久久| 精品人妻1区二区| 久久久精品国产亚洲av高清涩受| 日韩欧美 国产精品| 两性夫妻黄色片| 在线观看66精品国产| 色精品久久人妻99蜜桃| 岛国在线观看网站| a级毛片a级免费在线| 免费在线观看日本一区| www.999成人在线观看| 两个人看的免费小视频| 99久久精品热视频| 欧美极品一区二区三区四区| 搡老岳熟女国产| 又爽又黄无遮挡网站| 国产av在哪里看| 国产亚洲欧美98| 国产三级在线视频| 精品一区二区三区四区五区乱码| 伊人久久大香线蕉亚洲五| 可以在线观看的亚洲视频| 久久久久久九九精品二区国产 | 亚洲中文字幕一区二区三区有码在线看 | 床上黄色一级片| 国产激情欧美一区二区| 久久精品影院6| 国产探花在线观看一区二区| 免费高清视频大片| 亚洲欧美日韩高清在线视频| 欧美黄色片欧美黄色片| 久久国产精品人妻蜜桃| 免费在线观看完整版高清| 成人18禁高潮啪啪吃奶动态图| 88av欧美| 亚洲在线自拍视频| 中文字幕人妻丝袜一区二区| 国产又黄又爽又无遮挡在线| 1024香蕉在线观看| 中文在线观看免费www的网站 | 法律面前人人平等表现在哪些方面| 男女下面进入的视频免费午夜| 老司机福利观看| 午夜福利欧美成人| 亚洲男人的天堂狠狠| 亚洲片人在线观看| www.精华液| 免费av毛片视频| 蜜桃久久精品国产亚洲av| 国产精品av视频在线免费观看| 午夜福利成人在线免费观看| av欧美777| 色尼玛亚洲综合影院| 他把我摸到了高潮在线观看| 国产av一区在线观看免费| 丰满的人妻完整版| 正在播放国产对白刺激| 国产精品,欧美在线| xxx96com| 亚洲熟妇中文字幕五十中出| 中文字幕久久专区| 国产亚洲欧美在线一区二区| 在线国产一区二区在线| 91av网站免费观看| 少妇熟女aⅴ在线视频| 母亲3免费完整高清在线观看| 久久国产精品影院| 男插女下体视频免费在线播放| 日本黄色视频三级网站网址| 久久午夜亚洲精品久久| 18禁观看日本| 久久久精品国产亚洲av高清涩受| 亚洲真实伦在线观看| 少妇的丰满在线观看| 日韩 欧美 亚洲 中文字幕| 亚洲国产中文字幕在线视频| 香蕉av资源在线| 欧美一级a爱片免费观看看 | 精品高清国产在线一区| av在线天堂中文字幕| 超碰成人久久| 老鸭窝网址在线观看| 熟女电影av网| 国产熟女午夜一区二区三区| 一区二区三区高清视频在线| 99久久国产精品久久久| 村上凉子中文字幕在线| 欧美另类亚洲清纯唯美| 国产精品久久久人人做人人爽| 日本一本二区三区精品| 国产一区二区在线av高清观看| 欧美日本视频| 最近最新中文字幕大全免费视频| 一二三四在线观看免费中文在| 又大又爽又粗| 老汉色∧v一级毛片| 日韩欧美精品v在线| 可以免费在线观看a视频的电影网站| 日韩精品免费视频一区二区三区| 精品久久久久久久久久久久久| 精品乱码久久久久久99久播| 一级毛片精品| 91成年电影在线观看| 国产免费男女视频| 日韩av在线大香蕉| 亚洲 欧美 日韩 在线 免费| 欧美日韩一级在线毛片| av中文乱码字幕在线| 成人三级做爰电影| 他把我摸到了高潮在线观看| 怎么达到女性高潮| 中文字幕熟女人妻在线| 亚洲欧美日韩东京热| 1024视频免费在线观看| 三级国产精品欧美在线观看 | 热99re8久久精品国产| 天天躁夜夜躁狠狠躁躁| 最新在线观看一区二区三区| 午夜福利欧美成人| 国产aⅴ精品一区二区三区波| 亚洲av成人av| 男男h啪啪无遮挡| 亚洲欧美一区二区三区黑人| 精品少妇一区二区三区视频日本电影| 99在线视频只有这里精品首页| 亚洲精品在线观看二区| 欧美另类亚洲清纯唯美| 88av欧美| netflix在线观看网站| 人成视频在线观看免费观看| 动漫黄色视频在线观看| 三级毛片av免费| 中文字幕人成人乱码亚洲影| 午夜久久久久精精品| 久久久水蜜桃国产精品网| 亚洲专区字幕在线| 一边摸一边抽搐一进一小说| 天天躁夜夜躁狠狠躁躁| 婷婷丁香在线五月| 久久精品国产清高在天天线| 亚洲av中文字字幕乱码综合| 亚洲七黄色美女视频| 国产精品久久久久久久电影 | 丁香欧美五月| 亚洲天堂国产精品一区在线| 国产精品日韩av在线免费观看| 亚洲va日本ⅴa欧美va伊人久久| www.www免费av| 成人国产一区最新在线观看| 特大巨黑吊av在线直播| 亚洲男人天堂网一区| 国产成人精品久久二区二区免费| 国产成人精品无人区| 色哟哟哟哟哟哟|