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

    Conductivity effects during the transition from collisionless to collisional regimes in cylindrical inductively coupled plasmas

    2022-06-01 07:55:50WeiYANG楊唯FeiGAO高飛andYounianWANG王友年
    Plasma Science and Technology 2022年5期
    關(guān)鍵詞:王友

    Wei YANG (楊唯), Fei GAO (高飛) and Younian WANG (王友年)

    1 College of Science, Donghua University, Shanghai 201620, People’s Republic of China

    2 School of Physics, Dalian University of Technology, Dalian 116024, People’s Republic of China

    Abstract A numerical model is developed to study the conductivity effects during the transition from collisionless to collisional regimes in cylindrical inductively coupled argon plasmas at pressures of 0.1-20 Pa.The model consists of electron kinetics module, electromagnetics module, and global model module.It allows for self-consistent description of non-local electron kinetics and collisionless electron heating in terms of the conductivity of homogeneous hot plasma.Simulation results for non-local conductivity case are compared with predictions for the assumption of local conductivity case.Electron densities and effective electron temperatures under non-local and local conductivities show obvious differences at relatively low pressures.As increasing pressure, the results under the two cases of conductivities tend to converge, which indicates the transition from collisionless to collisional regimes.At relatively low pressures the local negative power absorption is predicted by non-local conductivity case but not captured by local conductivity case.The two-dimensional (2D) profiles of electron current density and electric field are coincident for local conductivity case in the pressure range of interest, but it roughly holds true for non-local conductivity case at very high pressure.In addition,an effective conductivity with consideration of non-collisional stochastic heating effect is introduced.The effective conductivity almost reproduces the electron density and effective electron temperature for the non-local conductivity case,but does not capture the non-local relation between electron current and electric field as well as the local negative power absorption that is observed for nonlocal conductivity case at low pressures.

    Keywords: inductively coupled plasmas, conductivity effects, electron kinetics, plasma parameters, electromagnetic wave characteristics, electrodynamics

    1.Introduction

    Low pressure inductively coupled plasmas (ICPs) are widely utilized for ion sources and material processing.Their popularity for both scientific and industrial applications depends on the capability to achieve high density even for relatively low power and high plasma uniformity, as well as on the possibility of independent control of the ion energy at the plasma boundary.

    In the past few decades, substantial experimental and theoretical efforts on the electron kinetics and the electron heating mechanisms have been made in order to understand low pressure ICPs and optimize the discharge characteristics in terms of the individual requirements [1, 2].The common feature of low pressure ICPs is their non-equilibrium nature.In particular,electrons are not only far from equilibrium with ions and neutrals, but even not in equilibrium within their own ensemble [3].Non-local electrodynamics and non-local electron kinetics as well as anomalous skin effect and local negative power absorption are distinctive features of the lowpressure regime [4-8].The differences between non-local electrodynamics and non-local electron kinetics have been introduced in detail by Kolobov et al[8].The electron energy probability function (EEPF) can deviate from a Maxwellian distribution in these discharges (non-local electron kinetics)[9, 10].Due to the smaller collision frequency of electrons with neutrals than excitation frequency, the local relation between electron current and electric field (Ohm’s law) does not hold true (non-local electrodynamics).In this case, the electron current relies on the distribution of the electric field in the whole plasma and size of the discharge device[11-13].It is due to the spatial dispersion of conductivity caused by electron thermal motion.Therefore, a non-local conductivity is required to determine the electric field penetration into the plasma.

    The accurate description of the discharge properties under low pressures relies on kinetic simulations, where the EEPF,non-local conductivity and plasma density need to be self-consistently calculated.The non-local conductivities for homogeneous plasmas were derived in planar-type [14, 15] and solenoidal-type [16] ICP reactors by Yoon et al.Due to assuming a constant electron density and a Maxwellian EEPF,only the averaged electron Boltzmann equation under linear approximation and the 2D Maxwell equations were solved in their models.Therefore, the interaction between electron properties and electric field cannot be self-consistently considered.Regarding the complexity of theoretical studies in cylindrical geometry,a self-consistent kinetic model with the calculations of nonlocal conductivity, EEPF, and non-uniform profiles of electron density in a one-dimensional (1D) slab geometry was developed by Kaganovich et al [12].Their model contains the averaged electron Boltzmann equation for EEPF, the Maxwell equation for the electric field, the quasi-neutrality condition for electrostatic potential, and the fluid conservation equations for ion density profile.Subsequently,the kinetic model was used to study the influences of non-local conductivity and non-Maxwellian EEPF on plasma power absorption and plasma density profiles in a 1D planar inductively coupled argon discharge under the non-local regime (pressure ≤10 mTorr) [17, 18].Recently, Kushner et al reported the influence of non-local conductivity on plasma absorption power in a 2D planar-type argon ICP at 10 mTorr,where electron Monte Carlo simulation instead of electron Boltzmann equation was used to calculate the electron current density [19].

    The aforementioned theoretical modellings in the collisionless regime have been well established, while the consideration of inhomogeneous plasmas in a 2D cylindrical geometry is very cumbersome.The fluid equations need to be solved and meanwhile the EEPF as a function of not only the kinetic energy but the potential energy needs to be taken into account.It might limit the studies to some fixed discharge conditions regarding the computational expense.To the best of our knowledge,the effects of plasma conductivity have not been fully reported over a wide pressure range that involves the transition from collisionless to collisional regimes.As mentioned earlier, plasma conductivity couples plasma kinetics to the electromagnetic fields.In this study,a numerical model based on previous work[20]that has been validated is developed to study the conductivity effects during the transition from collisionless to collisional regimes in cylindrical inductively coupled Ar plasmas.In the model, the electron kinetics are addressed by the electron Boltzmann equation, electromagnetic wave properties are obtained by the Maxwell’s equations,and plasma density is calculated from global model.The non-local conductivity,local conductivity, and effective conductivity are respectively considered in the model and their effects on plasma parameters and electromagnetic wave characteristics are compared under a wide pressure range of 0.1-20 Pa.The non-local conductivity is computed by the electron Boltzmann equation, while the local conductivity is obtained based on the Ohm’s law.The effective conductivity is introduced based on the Ohm’s law but with consideration of non-collisional stochastic heating effect by defining a stochastic collision frequency.In principle, the nonlocal conductivity model should be valid no matter which heating mechanisms (i.e.collisional, collisionless or both of them) dominate the discharge in a given situation.

    This manuscript is structured as follows.Section 2 shows a brief description of the model.Verification and validation of the model are shown in section 3.Simulation results and discussions are presented in section 4, where the non-local conductivity effects are compared with local conductivity case and with effective conductivity case.Section 5 presents the conclusions.

    2.Model description

    In this work, we consider a cylindrical ICP source with a solenoidal-type coil, and the structure of the chamber is the same as that of the driver chamber of a two-chamber radiofrequency (RF) ICP in the experiments by Li et al [21].The radius R is 6 cm and the length H is 14 cm.The model is based on the assumption of homogeneous plasma, since the discharge is assumed to be under inductive H mode to produce a thin sheath near the walls.

    The model mainly includes three modules of the electron kinetics module (EKM), electromagnetics module (EMM),and global model module(GMM).Harmonic electromagnetic field is generated by the alternating current in the coil by solving the frequency domain wave equation (i.e.EMM),where all physical quantities are assumed to be axisymmetric.For local model,electron current density is obtained based on Ohm’s law where the conductivity is local.For non-local model, the determination of electron current density needs a non-local conductivity that is deduced from the EKM (linearized electron Boltzmann equation for perturbed distribution function) with EMM.Electric field from the EMM is used in the EKM (electron Boltzmann equation under quasilinear approximation for isotropic distribution function) to produce normalized EEPF.The normalized EEPF is used to calculate electron impact rate coefficients and other related parameters for the use in the GMM.Particle balance equations for densities of particle species, and power balance equation for density of electrons are solved in the GMM.It is worth noting that only the parts involved in this work as well as variations and improvements in terms of our previous model are shown here and more details can be found in our previous work[20].

    2.1.Electromagnetics module

    The solution of the inductive electric field can be obtained by the mode analysis method using the Fourier series in a rectangular coordinate system(x,z) [20]

    whereknandqmcan refer to the previous work [20], and the expansion coefficientEnmof the electric field is given as[20]

    whereμ0,ω,I0,N, andzkare respectively the vacuum permeability,excitation frequency,current amplitude in the coil,number of turns of the coil,and location of thekth turn in the axial direction.k0iswith speed of light c.The conductivityσwill be discussed in the EKM.

    Based on the solution of electric field, the electron current densityJpcan be expressed in terms ofEyvia a conductivityσ, and it is given as [20]

    2.2.Electron kinetics module

    In this part, the non-local conductivity, local conductivity,and effective conductivity, as well as normalized EEPF will be presented.

    2.2.1.Non-local conductivity.In the non-local limit, where the electron-neutral collision frequencyνenis of the same order as (or even smaller than) the excitation frequencyω,collisionless heating mechanism might compete with or even dominate over collisional heating in the low-pressure regime.The Ohm’s law for electron current densityJpis therefore considered to be not valid anymore.To account for collisionless electron heating,Jpis calculated based on the perturbed distribution functionf1,and it is given by

    whereeis the elementary charge,neis the electron density,andvyis the velocity component in the y direction.

    The non-local conductivity can be self-consistently computed by equations (3) and (4).σnis thenth coefficient of non-local conductivity in the Fourier series,and it is given as [20]

    whereεis the electron kinetic energy,meis the electron mass,is the electron velocity, andf0is isotropic distribution function(i.e.the normalized EEPF).Ψ andΦ can refer to the previous work [20].

    2.2.2.Local conductivity.When the electron-neutral collision frequencyνenis much larger than the excitation frequencyωat higher pressures, the conductivity is defined by the Ohm’s law.It indicates that the collisional heating is the only manner that the electromagnetic field transfers random thermal energy to the plasma.In cases where the local relation between electron current and electric field is valid,the local conductivity for an arbitrary EEPF is given by[17, 22, 23]

    The expansion coefficientEnm(see equation (2)) of electric fieldEyand electron current densityJp(see equation(3))mentioned in the EMM are therefore determined based on the three cases of conductivities.

    2.2.4.Normalized EEPF.The normalized EEPFf0can be solved by

    where the energy diffusion coefficientDeis deduced in appendix in detail and is expressed as

    Table 1.Reactions considered in the model.

    where the first term on the RHS has the same form as the term on the LHS of equation (8).DeeandVeein equation (10) are respectively given as [25]

    whereveeis the frequency of electron-electron collision [25].

    The normalization condition foris required.Equation (8) forf0can be written in the form of an ordinary nonlinear differential equation[25],and solved using the tridiagonal speed-up method.

    2.3.Global model module

    Global model is also called zero-dimensional model, since it ignores the spatial variation of plasma properties.As known,the global model usually includes particle balance equations for each species except for the background gas, the quasineutrality condition, the equation of state for total particle number density, and the power balance equation.Therefore,the global model can be used for the calculation of electron temperature and densities of all particle species.In this study,the normalized EEPFf0(effective electron temperatureTeff)is calculated using the electron Boltzmann equation under quasilinear approximation, so the global model is only used for the calculation of the densities of particle species.We discard one of the particle balance equations due to a reduced unknown,i.e.electron temperature in the global model.In our previous work [20], only the power balance equation that is used to determine the electron density is considered.In this work, the particle balance is also included to calculate densities of other particle species, because the Ar metastable is also considered.Even if including Ar metastable has very small effects on the plasma parameters and electromagnetic field properties in this work, the consideration of particle balance equation will help to extend the model to more complex gas discharges for future work.The particle balance equation other than background gas is given as [1, 26]

    whereniis the density of the ith reactant.For the electronimpact volume reactions,kis calculated based on the normalized EEPF [1]

    For the surface reaction for neutrals and positive ions,the rate coefficients can be found in [26].In this work, the reaction processes of Ar gas are shown in table 1.The cross sections of electron-impact elastic scattering, excitation, de-excitation and ionization are in [27-29].

    The density of background gas is determined by the equation of state.Meanwhile, the discharge is assumed to satisfy the quasi-neutrality condition.These two criteria are respectively given by

    whereTgasis the gas temperature andkBis the Boltzmann constant.In this work, equation (17) is simplified asne=

    The power balance equation in our previous work[20]is based on the assumption of particle balance by equating the total surface particle loss to the total volume ionization[1].In fact, it is not quite appropriate, since the EEPF is not strictly Maxwellian and the electron temperature is not determined by the particle balance equation but by the electron Boltzmann equation under quasilinear approximation.It leads to the underestimation of electron density predicted by simulations[20]as increasing pressure when comparing with experiments[21], due to the overestimation of electron loss rate at relatively high pressures.In this work, the power balance equation is utilized in a more general form.The absorption power equals to the bulk power losses caused by electronneutral collisions, and surface power losses due to charged species flowing to the walls [26]:

    where the absorption powerPabsis obtained from the EMM.In this work, we only have one kind of ions, i.e.Ar+, so the ion densityniequals to electron densityne.PVis the bulk power loss of an electron per unit volume andPSis the surface power loss of an ion on the walls per unit volume, respectively.They are given as [26]

    3.Verification and validation of the model

    3.1.Verification of the model

    In order to verify the model,it is compared with the non-local approach to the solution of the Boltzmann equation developed by Kolobov et al [34].In their work, electron kinetics in a cylindrical ICP sustained by a long solenoidal coil was studied for the near-collisionless regime when the electron mean free path was large compared to the radius of the chamber.For simplicity, they prescribed 1D profiles of the fields and skin depth rather than self-consistently solving the Maxwell equations with electron kinetics and the plasma density.The amplitude of magnetic field (B0= 2 G), skin depth(δ= 1 cm) related to plasma density, radius of the chamber(R = 5 cm), and RF frequency (f= 13.56 MHz) are the input parameters in their model.The former two quantities are self-consistently calculated in our model.In order to present more reasonable comparison, we try to find the discharge condition capable of making the results of field and plasma density closer to the input ones used in the model of Kolobov et al and the electron mean free path larger than the radius of the chamber [34].In our simulation, the same radius and RF frequency are used.The length of the chamber is 15 cm, the pressure is 0.5 Pa and the absorption power is 80 W.Under the discharge condition,the calculated skin depthδis around 0.95 cm and the amplitude of electric fieldE0is 4 V cm?1that corresponds to the amplitude of magnetic fieldB0of around 1.9 G in our model.

    Figure 1 shows that the normalized EEPF predicted by our model is compared with that by the model of Kolobov et al[34].In their model,the electron-electron collisions were excluded.The normalized EEPFs with and without consideration of electron-electron collisions are predicted by our model.Note that two models achieve reasonable agreement with each other.The model of Kolobov et al predicts a more remarkable depletion at the electron energy higher than 18 eV and more cold electrons with the energy lower than 2.5 eV.The differences between the two models can be due to the different collision processes considered in the two models.The collision processes included in the model of Kolobov et al were not mentioned [34].If more inelastic collision processes were included in the model of Kolobov et al, it would lead to more depletion in the energetic electrons and enhancement in low-energy ones.As expected, our model without consideration of electron-electron collisions agrees better with the model of Kolobov et al.The inclusion of the electron-electron collisions decreases cold electrons with the energy lower than 2.5 eV, and slightly increases energetic electrons with the energy higher than 18 eV.In general,electron-electron collisions tend to Maxwellize the EEPF.

    Figure 1.Normalized EEPF versus electron energy.Predictions by the model in this work with and without electron-electron collisions and by the model of [34] without electron-electron collisions.Reprinted with permission from [34].Copyright (1997) by the American Physical Society.

    In addition,the effects of chamber lengths(i.e.5,10,15,and 20 cm)have also been studied,even if they are not shown here.The simulated EEPF shows that the energetic electrons in inelastic energy range slightly increase with increasing chamber length, because the increased electron energy relaxation length decreases electron energy loss via the inelastic collisions.It leads to slightly higher effective electron temperature at larger chamber length.Shorter chamber length, e.g.5 cm, makes its scale smaller than the energy relaxation length, but the increased electron-electron collisions due to the increased plasma density still tend to Maxwellize the EEPF.This is because in the model the smaller chamber length produces larger power deposition density and thus higher plasma density for fixed total absorption power.

    3.2.Validation of the model

    In the previous study [20], the EEPF (electron kinetics)as well as the electron temperature and electron density (basic plasma properties) predicted by the model, was validated with the experiments[21]under the same discharge conditions.Here,we validate the model against experimental measurements in a cylindrical ICP discharge, where the data on electromagnetic wave characteristics under different discharge conditions are available[7].The experiment is composed of a Pyrex glass tube(Φ 26.3 cm (ID)×45 cm×0.5 cm) and a lower stainless-steel chamber (Φ 45 cm×25 cm×0.5 cm).The ICP source was sustained by a one-turn RF coil made of a copper strip(Φ 32 cm×6 cm×0.05 cm),and the vertical distance between the lower side of the copper strip and the upper side of the stainless-steel chamber was 14.9 cm.The plasma source was shielded with an aluminium cage to suppress RF electromagnetic interference.It is worth noting that the plasma can be very non-uniform especially along the axial direction for this kind of coil position and chamber size.The radial distribution of the RF magnetic field(axial component)was measured using the RF magnetic probe at the center of strip coil along the axial direction.The azimuthal electric field was obtained based on the measured magnetic field.

    Figures 2(a) and (b) show the 1D radial distributions of electric fields versus power obtained from the experiments and simulations, respectively.Both experiments and simulations show that the electric field rapidly decays from the boundary near the coil to the center of chamber.The electric field decreases faster in the experiment than that in the simulation.This is because the measurements have been performed at the center of strip coil along the axial direction where the plasma density is almost the maximum.Nevertheless, the uniform plasma density assumed in the simulations can possibly make the volume-averaged density smaller than the one in the experiments.Therefore, the lower plasma density in the simulations leads to the deeper penetration and slower decay of the electric field into the plasma.As expected,the lower power facilitates the penetration of electric field into the plasma, which is captured by both simulations and experiments.Note that the higher power produces the nonmonotonic field decay,which is typically the anomalous skin effect caused by the collisionless electron thermal motion.The phenomenon has also been observed and discussed in the experiments of an ICP source driven by a planar coil [35].The non-monotonic evolution of the electric field appears at higher power in the simulation, due to higher power required in the simulation to produce the comparable plasma density to that in the experiment.

    Figure 2.Radial distribution of normalized electric field under different absorption powers at 0.1 Pa.Experimental data (a)reprinted with the permission from [7], copyright (2015), AIP Publishing LLC and predictions by the model in this work (b).

    4.Results and discussion

    4.1.Non-local conductivity effects

    In this section, the effects of non-local and local conductivities on plasma parameters and electromagnetic wave characteristics are compared.The non-local and local conductivities are respectively obtained by solving the electron Boltzmann equation and by using the local approximation.In the simulations, the absorption power is 50 W, the RF excitation frequency is 13.56 MHz, and the pressure varies from 0.1 to 20 Pa.The conductivity is expressed asσin the following figures.

    The characteristic lengths play a key role in the definition of non-local and local electrodynamics [8].If the characteristic size of plasma is much larger than electron mean free pathλbut is comparable to electron energy relaxation lengthλε, in the collisional (but nonlocal) regime the isotropic part of the electron energy distributionf0at a given point relies not only on the electric field at this point but also on plasma properties in the vicinity of the point of the size, while the anisotropic partf1is a local function of the field (i.e.local electrodynamics).Ifλbecomes comparable or even larger than the characteristic size of the plasma.In this nearly collisionless regime, thef1at a point is determined not only by the local value of the electric field at this point,but also by the profile of the electric field in the vicinity of the point of the size along the electron trajectory (i.e.nonlocal electrodynamics) [8].

    Figure 3.Electron energy relaxation length and mean free path versus pressure.

    The electron density and effective electron temperature versus pressure under non-local and local conductivities are shown in figure 4.In figure 4(a), electron density monotonously increases with increasing pressure, which achieves qualitatively agreement with the previous experiments [21].The electron density of non-local conductivity is higher than that of local conductivity at relatively low pressures.It also shows one of the important reasons why the fluid simulation is not valid anymore at very low pressures.The fluid simulations are based on the Ohm’s law where only the collisional heating considered is not enough to sustain discharge in lowpressure regime.For example, the electron density of nonlocal conductivity is far higher than that of local conductivity at 0.1 Pa in our simulation.As increasing pressure, the electron densities under two cases of conductivities tend to converge, which indicates the transition from non-local to local electrodynamics regime.As expected, the effective electron temperature of non-local conductivity is lower than that of local conductivity at relatively low pressures,as shown in figure 4(b).Similar to the variation of electron density,the effective electron temperatures under two cases of conductivities also tend to converge as increasing pressure.The simulation results indicate that the consideration of the nonlocal conductivity in a low-pressure regime is very important for the development of a numerical model.

    Figure 4.Electron density(a)and effective electron temperature(b)versus pressure under non-local and local conductivities.The excitation frequency is 13.56 MHz and the absorption power is 50 W.

    The effective electron temperature is obtained based on the integral of normalized EEPF.Figure 5 shows the normalized EEPF versus electron energy at different pressures under non-local and local conductivities.Two cases of conductivities result in a big difference in the normalized EEPF,especially at 0.1 Pa, where the local conductivity produces more electrons in the inelastic energy range and fewer electrons in the elastic energy range.The shape of EEPF is determined by the expansion coefficients in the Fourier series of the electric field, according to equation (9), where the energy diffusion coefficient is proportional to the square of the expansion coefficients of electric field.The effects of collisionless heating under non-local conductivity case allow the electromagnetic wave to couple energy to plasma more efficiently.As a result, the electric field required to deposit a given power into the plasma is smaller than that needed when only considering the collisional heating under local conductivity case.The simulation results of electric fields will be shown later.The stronger electric field under local conductivity case leads to a larger energy diffusion coefficientDein the simulations and therefore more energetic electrons penetrate into the RF field region.Consequently, local conductivity produces more electrons in the inelastic energy range.

    Figure 5.Normalized EEPF versus electron energy at different pressures under non-local and local conductivities.The discharge conditions refer to figure 4.

    In the simulations, the energy diffusion coefficientDeis smaller than the energy diffusion coefficient related to electron-electron collisionsDee(see equation (11)) at low electron energy.Therefore, the electron-electron interactions are the only mechanism of energy exchange for the low energy electrons that cannot reach the region of strong electric field where the heating occurs.Deeunder non-local conductivity case is larger than that under local conductivity case due to higher electron density (see figure 4(a)).It leads to more low energy electrons and meanwhile tends to Maxwellize the EEPF in elastic energy range at low pressure of 0.1 Pa under non-local conductivity case.As increasing pressure, the collisionless heating is less efficient and the thermalization of EEPF occurs.The normalized EEPFs evolve into a Maxwellian distribution and are almost identical at 20 Pa for two cases of conductivities.This is because the model-predictedDeandDeefor two cases of conductivities are almost the same at 20 Pa.

    Figure 6 shows the 2D profile of the absorption power density versus pressure under non-local conductivity case.The absorption power density is taken the absolute values and plotted on a log scale to emphasize the negative power absorption.The local negative power absorption in the plasma is predicted by the model with consideration of non-local conductivity at low pressures, for example, 0.1 and 0.5 Pa,and it still exists at 1 Pa even if it is not shown in figure 6.It is well known that the negative power absorption is caused by the thermal motion of electrons.For the case of non-local conductivity, the phase difference between electric field and electron current is not constant.The energetic electrons originating from the discharge edge can be out of phase with an electric field outside the skin layer (close to the center), and lose energy to electric field, thus resulting in negative power absorption.For example, negative power absorption appears close to the discharge center outside the skin layer at 0.1 and 0.5 Pa.The simulation results confirm that there is a spatial region in the discharge where the electrons, instead of taking energy from the field,give back energy to the electromagnetic wave.The size and the number of negative and positive power absorption region depend on the power, which has been experimentally confirmed by Cunge et al[36]and Ding et al [7].At low power, there is one area of positive and one area of negative power absorption in the radial direction.As increasing power,the number of negative and positive power absorption region increases.Due to a low absorption power of 50 W used in our simulations, there is only one area of positive and one area of negative power absorption.

    At higher pressures,for example,3 and 20 Pa,the region of negative power absorption vanishes.It implies that the dominant power transfer mechanism is collisional heating at higher pressures.Even if homogeneous plasma is assumed in the model, the important physical property, i.e.the local negative power absorption, is still captured by the model.It has been demonstrated in the work of Godyak and Kolobov,when comparing the 1D profile of the absorption power density obtained in the experiments with that in the calculation assuming homogeneous plasma [6].The position and amplitude of positive and negative power densities predicted by the model can be very possibly different from the reality in the experiments where the plasma is inhomogeneous.The uniformity of plasma is better at lower pressure and power.It is well known that the skin depth inversely relies on the electron density [1].The magnitude and profile of plasma density can affect the magnitude and profile of the electromagnetic fields, electron current, and absorption power density.The absorption power density can be underestimated at the center of the discharge chamber and overestimated at the boundary of the discharge chamber.In addition, the power density profile shrinks and its value becomes larger towards the boundary near the coil, due to the increase in electron density (see figure 4(a)) as increasing pressure.

    The 2D profile of the absorption power density versus pressure under local conductivity case is shown in figure 7.Very different from the non-local conductivity case,the local negative power absorption is not captured by the local conductivity case at low pressures.For the case of local conductivity, the phase difference between electron current and electric field is constant.Thus,the profile of absorption power relies only on the magnitudes of electric field and local conductivity.At higher pressures, for example, 20 Pa, the 2D profile of the absorption power density under local conductivity case is almost identical with that under non-local conductivity case.It indicates that the collisional heating dominates the discharge at higher pressures.Similar to the non-local conductivity case, the amplitude of absorption power density increases with increasing pressure under local conductivity case due to the increased electron density.As expected,the difference of amplitude in the absorption power density under two cases of conductivities significantly decreases with increasing pressure.For example, for nonlocal and local conductivities the amplitudes of absorption power density are 0.62 W cm?3and 0.28 W cm?3at 0.1 Pa and 1.6 W cm?3and 1.8 W cm?3at 20 Pa.

    Figure 6.2D profiles of absorption power density versus pressure under non-local conductivity.The discharge conditions refer to figure 4.Power density is taken the absolute values and plotted on a log scale to emphasize the negative power absorption.

    Figure 7.2D profiles of absorption power density versus pressure under local conductivity.The discharge conditions refer to figure 4.

    The 2D profiles of electric field and electron current density versus pressure under non-local conductivity are shown in figure 8.As increasing pressure, the 2D profiles of electric field and electron current density shrink towards the boundary near the coil,and the amplitude of the electric field decreases that is consistent with variation of the effective electron temperature with pressure, but the amplitude of the electron current density increases due to the increase in electron density.The non-local relation between electron current density and electric field is more pronounced at lower pressures represented by very different 2D profiles of them.The skin depthsδpredicted by the model are respectively 1.5,1.0, 0.7 and 0.4 cm at 0.1, 0.5, 3 and 20 Pa.At lower pressures,more electrons contribute to the electron current outside the skin layer and near the plasma bulk, which is a non-local behavior.As increasing pressure,the electron current tends to be determined by the local electric field.Due to homogeneous plasma assumed in the model, the 2D profiles of electron current density and electric field can be slightly different from the reality where the electron density usually has maximum in the discharge center and decreases at the boundary,especially for higher pressures.

    The 2D profiles of electric field and electron current density versus pressure for local conductivity case are shown in figure 9.It is clearly shown that the electron current is determined by the local electric field in the investigated pressure range due to the assumption of Ohm’s law between electron current density and electric field.The 2D profiles of electron current density and electric field are identical to each other.The skin depthsδpredicted by the model are respectively 2.3,1.2,0.7 and 0.4 cm at 0.1, 0.5, 3 and 20 Pa.Note thatδunder local conductivity case with lower electron density (see figure 4(a)) is larger than that under non-local conductivity case (mentioned in figure 8) at lower pressures.As increasing pressure to 20 Pa,the skin depths are the same(i.e.0.4 cm)under two cases of conductivities due to the same electron density.

    Figure 8.2D profiles of induced electric field (first row) and plasma current density (second row) versus pressure under non-local conductivity.The discharge conditions refer to figure 4.

    The amplitudes of electric field under local and non-local(see figure 8) conductivities are very closer with each other at higher pressures and very different at lower pressures.The amplitude of the electric field under local conductivity case is far higher than that under non-local conductivity case at lower pressures.As mentioned in figure 5, for the case of non-local conductivity the effects of warm-electron and long mean free path cause the electromagnetic wave to couple energy to the plasma more efficiently.The lower electric fields are therefore required to deposit the same total power.It leads to the higher electron density and lower effective electron temperature under non-local conductivity case (see figure 4).For example, at 0.1 Pa,the amplitudes of electric field under local and non-local conductivities are respectively 12.4 V cm?1(here) and 5.02 V cm?1(see figure 8), therefore the effective electron temperatures are respectively 8.6 eV and 5.4 eV (see figure 4(b)).The consideration of non-local conductivity by solving electron Boltzmann equation in the model is very important for the prediction of non-local relation between electron current and electric field as well as the local negative power absorption in low-pressure regime.

    In order to estimate the variation of heating mechanisms with pressure,the characteristic frequencies of plasmas versus pressure are presented in figure 10.These frequencies include electron-neutral collision frequencyνen,stochastic collision frequencyνstoc, and effective collision frequencyνeff=νen+νstoc.In addition, the value of the excitation frequencyωis also shown.Bothνenandνstocare obtained based on the results of the model with consideration of non-local conductivity.The collision frequencyνenalmost linearly increases with increasing pressure because of its proportionality to the density of neutral particles.νenis larger than excitation frequencyωat higher pressures (≥6 Pa), thus the Ohmic heating dominates the discharge.The stochastic collision frequencyνstocshows a slight increase with pressure.Comparingνenwithνstoc, the Ohmic heating is expected to be dominant at higher pressures asνenexceedsνstocby almost one order of magnitude.Due to almost negligible effect of stochastic heating at higher pressures,νefftends to be equivalent toνen.

    Figure 9.2D profiles of induced electric field (first row) and plasma current density (second row) versus pressure under local conductivity.The discharge conditions refer to figure 4.

    Figure 10.Characteristic frequencies of electrons and excitation frequency (ω )versus pressure.The characteristic frequencies include the electron-neutral collision frequency (ν en ),the stochastic heating frequency (νs toc ),and the effective collision frequency (ν eff ).The discharge conditions refer to figure 4.

    Towards the lower pressures (≤3 Pa), the influence of stochastic heating is important.As revealed by the simulation results of absorption power density (see figure 6), the local negative power absorption vanishes at 3 Pa.At 0.8 Pa,νenandνstocare the same magnitude.At the lowest pressure limit of 0.1 Pa, because ofνstocfar higher thanνen,νeffsignificantly deviates fromνenand tends to be equivalent toνstoc.Therefore,the evaluation of the characteristic frequencies shows a transition of the heating mechanism from non-collisional stochastic heating to collision heating as increasing pressure from 0.1 to 20 Pa.It is necessary to consider stochastic heating in low-pressure regime.

    4.2.Effective conductivity effects

    Due to the importance of collisionless heating caused by the thermal motion of electrons at low pressures, the effective conductivity (equation (7)) is introduced based on the local conductivity (equation (6)) but accounts for the stochastic heating similar to collisional heating by defniing a stochastic collision frequencyνstoc.In order to evaluate the feasibility of effective conductivity at low pressures, its effects on plasma parameters and electromagnetic wave characteristics are compared with non-local conductivity case.

    Figure 11.Electron density (a) and effective electron temperature (b) versus pressure under non-local and effective conductivities.The discharge conditions refer to figure 4.

    Figure 12.Normalized EEPF versus electron energy at different pressures under non-local and effective conductivities.The discharge conditions refer to figure 4.

    The electron density and effective electron temperature versus pressure under effective and non-local conductivities are shown in figure 11.Different from the comparison between local and non-local conductivities,i.e.figure 4,where very pronounced discrepancies of effective electron temperature and electron density at lower pressures are presented, the effective conductivity almost reproduces the effective electron temperature and electron density under non-local conductivity case in the investigated pressure range.The addition of collisionless heating to collisional heating by introducing the stochastic collision frequency effectively increases electron density compared with the case of only collisional heating taken into account in low-pressure regime.Figure 12 shows the normalized EEPF versus electron energy at different pressures under effective and non-local conductivities.As expected, the normalized EEPF under non-local conductivity case is also almost reproduced by the effective conductivity.It confirms that the effective collision frequency instead of electronneutral collision frequency taken into account in the solver of local conductivity achieves a more accurate description of electron density and effective electron temperature in low-pressure regime.

    The 2D profile of absorption power density versus pressure under effective conductivity case is shown in figure 13.The effective conductivity almost reproduces the skin depth and the 2D profile of electric field under non-local conductivity case at low pressures of 0.1 and 0.5 Pa, and therefore the effective electron temperature and electron density (see figure 11).However,it does not capture the local negative power absorption that is revealed by the non-local conductivity at low pressures.Meanwhile, the non-local relation between electron current and electric field at low pressures is also not captured by the effective conductivity,as shown in figure 14.This is because despite the effective collision frequencyνeffinstead of the collisional frequencyνenused in the calculation of conductivity, the determination of electron current density still depends on the Ohm’s law.The nature of the local relation between electron current and electric field is not changed for an effective conductivity case.

    Figure 13.2D profiles of absorption power density versus pressure under effective conductivity.The discharge conditions refer to figure 4.

    Figure 14.2D profiles of induced electric field (first row) and plasma current density (second row) versus pressure under effective conductivity.The discharge conditions refer to figure 4.

    5.Conclusions

    Plasma conductivity is a fundamental parameter defining the electrical characteristics of a plasma device.It couples plasma kinetics to electromagnetic fields.In this study, a numerical model for homogeneous plasma is developed to study the conductivity effects in cylindrical inductively coupled Ar discharges at pressures of 0.1-20 Pa.The discharges experience a transition from collisionless to collisional regimes as the pressure increases from 0.1 to 20 Pa.The model consists of an electron kinetics module, electromagnetics module, and global model module.The non-local conductivity, local conductivity, and effective conductivity are respectively considered in the model and their effects on plasma parameters and electromagnetic wave characteristics are compared.The non-local conductivity is addressed by the electron Boltzmann equation.The local conductivity is defined by relating electron current to the electric field using Ohm’s law.The effective conductivity is introduced based on the Ohm’s law but with consideration of non-collisional stochastic heating effect by defining a stochastic collision frequency.

    Simulation results under non-local conductivity cases are first compared with predictions under the assumption of local conductivity case (Ohm’s law).The lower effective electron temperature, as well as higher electron density, is obtained for non-local conductivity case than that for local conductivity case at relatively low pressures.This is because the collisionless heating effects under non-local conductivity case allow the electromagnetic wave to couple energy to the plasma more efficiently.Consequently,the electric field required to deposit a given power into the plasma is smaller than that needed when only considering the collisional heating under local conductivity case.The weaker electric field under non-local conductivity case leads to lower effective electron temperature and therefore higher electron density at lower pressures.As increasing pressure, the effective electron temperature and electron density,respectively,tend to converge under two cases of conductivities.It displays a transition of electron kinetics from non-local to local regime.The local conductivity predicts that there is no negative power absorption and the absorption power density monotonically decays from the boundary near the coil towards the discharge center in the investigated pressure range.However,the non-local conductivity predicts a local negative power absorption region and two peaks respectively for positive and negative power absorption regions at relatively low pressures.Almost the same 2D profiles of the absorption power density for non-local conductivity as that for local conductivity are obtained at relatively high pressures.This difference between two cases of conductivities at relatively low pressures can be understood by the prediction of the 2D profiles of electron current density and electric field.More spread-out 2D profile of electron current density than that of electric field under non-local conductivity case indicates a large fraction of electrons coming from the skin layer contributing to the electron current outside the skin layer.The 2D profiles of electron current density and electric field are coincident for local conductivity case in the pressure range of interest,but it only roughly holds true for non-local conductivity case at very high pressure.It verifies that the local relation between electron current and electric field is not valid anymore in the low-pressure regime.The non-local conductivity model is valid no matter which heating mechanisms(i.e.collisional,collisionless or both of them) dominate the discharge in a given situation.

    In addition, an effective conductivity with consideration of non-collisional stochastic heating effect is introduced.The effective conductivity almost reproduces the electron density,effective electron temperature,and normalized EEPF under nonlocal conductivity case in the investigated pressure range.However, the 2D profiles of electric field and electron current density are coincident and there is no negative power absorption in the investigated pressure range.Therefore, the non-local relation between electron current and electric field is not captured by the effective conductivity.This is because the introduction of effective conductivity is still based on the Ohm’s law and therefore the nature of the local relation between electron current and electric field is not changed.Even so, for the calculations of effective electron temperature and electron density,the effective conductivity can be alternative to non-local conductivity under the non-local regime in complex modellings such as multi-dimensional fluid simulations for simplicity.

    In this work,the effect of the ponderomotive force is not investigated, because it has been proven to be negligibly small at the industrial frequency 13.56 MHz [37].In future work, the model will be developed to describe the molecular gas discharges, such as hydrogen gas discharge that are extensively used in negative hydrogen ion sources operated at lower RF frequency (1 MHz).In that case, the magnetic(ponderomotive)force can be rather large,because it is almost inversely proportional to the RF frequency [34, 38].

    Acknowledgments

    This work was sponsored by National Natural Science Foundation of China (Nos.12105041, 11935005 and 12035003), Fundamental Research Funds for the Central Universities (No.2232020D-40) and Shanghai Sailing Program (No.20YF1401300).

    Appendix.Derivation of the energy diffusion coefficient in the kinetic equation

    In order to clarify the derivation of energy diffusion coefficientDεin equation (8) of the main text, the deduction is shown in detail here.A quasilinear equation for isotropic distribution functionf0is given as [12, 20]

    whereEandBare respectively the electric field and magnetic field.For arbitrary physical quantities pertaining to simple harmonic oscillation with timeA(t)=Ae-iωtandB(t)=Be-iωt,by averaging their product over the oscillation periodT=we have

    Therefore, the kinetic equation forf0averaged over the RF period is give as

    Then,using the orthogonality conditions[20],the kinetic equation forf0averaged over the space is given as

    Substituting equation (A6) into (A5), we have

    where the flux in the velocity spaceΓ is given as

    Using the spherical coordinates in the velocity space

    we have

    Averaging the?v·Γover azimuthal angle in the velocity space, we obtain

    where

    Substituting equation (A15) into (A14), we have

    where the functionΨ can be found in [20].

    Substituting equation (A16) into (A12) and combining equation (A8), we can obtain

    where

    For convenience,the RHS of equation(A17)is expressed as a function of energyε.It is given as

    where the energy diffusion coefficientDein equation (8) of the main text is obtained as follows

    Finally, the effect of magnetic field is cancelled out, because the Lorentz force term (v×B) is averaged over azimuthal angle in the velocity space (see equation (A12)).

    ORCID iDs

    猜你喜歡
    王友
    3D fluid model analysis on the generation of negative hydrogen ions for negative ion source of NBI
    Fundamental study towards a better understanding of low pressure radio-frequency plasmas for industrial applications
    High energy electron beam generation during interaction of a laser accelerated proton beam with a gas-discharge plasma
    Multi-layer structure formation of relativistic electron beams in plasmas
    Influence of magnetic filter field on the radiofrequency negative hydrogen ion source of neutral beam injector for China Fusion Engineering Test Reactor
    Numerical investigation of radio-frequency negative hydrogen ion sources by a three-dimensional fluid model?
    Modulation of the plasma uniformity by coil and dielectric window structures in an inductively coupled plasma
    Time-resolved radial uniformity of pulse-modulated inductively coupled O2/Ar plasmas?
    Experimental investigation of the electromagnetic effect and improvement of the plasma radial uniformity in a large-area,very-high frequency capacitive argondischarge
    Spatio-temporal measurements of overshoot phenomenon in pulsed inductively coupled discharge?
    最近2019中文字幕mv第一页| 国产精品久久久久成人av| 亚洲av免费高清在线观看| 久久精品夜色国产| 午夜影院在线不卡| 美女脱内裤让男人舔精品视频| 亚洲熟女精品中文字幕| 亚洲精品久久午夜乱码| 欧美日韩精品网址| 欧美97在线视频| 欧美最新免费一区二区三区| av片东京热男人的天堂| videossex国产| 国产成人欧美| 国产精品.久久久| 成人手机av| 夫妻午夜视频| 少妇 在线观看| 亚洲国产看品久久| 寂寞人妻少妇视频99o| av免费在线看不卡| 最近手机中文字幕大全| 一级爰片在线观看| 国产成人精品久久二区二区91 | 亚洲一级一片aⅴ在线观看| 91精品国产国语对白视频| 亚洲av综合色区一区| 亚洲人成网站在线观看播放| 91午夜精品亚洲一区二区三区| 免费不卡的大黄色大毛片视频在线观看| 99久久综合免费| 亚洲国产成人一精品久久久| 日韩熟女老妇一区二区性免费视频| 国产精品成人在线| 久久久久人妻精品一区果冻| 青春草国产在线视频| 午夜福利乱码中文字幕| 美女大奶头黄色视频| 免费不卡的大黄色大毛片视频在线观看| 国语对白做爰xxxⅹ性视频网站| 久久精品久久久久久久性| 考比视频在线观看| 国产精品成人在线| 高清不卡的av网站| videos熟女内射| 亚洲av电影在线进入| 日韩伦理黄色片| av国产精品久久久久影院| 天堂俺去俺来也www色官网| 久久毛片免费看一区二区三区| 亚洲精品aⅴ在线观看| xxx大片免费视频| 亚洲精品国产av成人精品| 精品亚洲乱码少妇综合久久| 国产成人精品在线电影| 日本wwww免费看| 日韩制服骚丝袜av| 午夜老司机福利剧场| 欧美成人午夜免费资源| av.在线天堂| 女人精品久久久久毛片| www日本在线高清视频| 国产精品麻豆人妻色哟哟久久| 国产精品av久久久久免费| 亚洲国产av影院在线观看| 91精品伊人久久大香线蕉| 天天躁夜夜躁狠狠躁躁| 亚洲在久久综合| 男女下面插进去视频免费观看| 亚洲av中文av极速乱| 成人午夜精彩视频在线观看| 熟女av电影| 99精国产麻豆久久婷婷| 亚洲av国产av综合av卡| 大话2 男鬼变身卡| 久久精品亚洲av国产电影网| 日韩一卡2卡3卡4卡2021年| 国产精品欧美亚洲77777| 久久国产亚洲av麻豆专区| 夜夜骑夜夜射夜夜干| 人妻系列 视频| 韩国精品一区二区三区| 美女中出高潮动态图| 热99久久久久精品小说推荐| 一级,二级,三级黄色视频| 在线观看免费日韩欧美大片| 人人妻人人澡人人看| 亚洲一级一片aⅴ在线观看| 久久影院123| av在线播放精品| 丝袜美足系列| 狂野欧美激情性bbbbbb| 国产精品一二三区在线看| 国产日韩欧美视频二区| 女的被弄到高潮叫床怎么办| 高清在线视频一区二区三区| 91成人精品电影| 久久精品aⅴ一区二区三区四区 | 欧美老熟妇乱子伦牲交| 精品人妻熟女毛片av久久网站| 看十八女毛片水多多多| 欧美激情高清一区二区三区 | 欧美日韩一区二区视频在线观看视频在线| 亚洲欧美一区二区三区黑人 | 99久久精品国产国产毛片| 水蜜桃什么品种好| 超碰97精品在线观看| 91成人精品电影| 久久亚洲国产成人精品v| 丝袜美腿诱惑在线| 欧美xxⅹ黑人| 国产黄色免费在线视频| 婷婷色综合大香蕉| 亚洲一级一片aⅴ在线观看| 亚洲欧洲精品一区二区精品久久久 | 天天影视国产精品| 亚洲视频免费观看视频| 欧美激情极品国产一区二区三区| 天堂中文最新版在线下载| av在线观看视频网站免费| 色吧在线观看| 免费av中文字幕在线| 看免费av毛片| 亚洲av免费高清在线观看| 热re99久久精品国产66热6| av在线老鸭窝| 久久久精品免费免费高清| 人人妻人人澡人人爽人人夜夜| 纵有疾风起免费观看全集完整版| 国产精品国产三级专区第一集| 免费播放大片免费观看视频在线观看| 有码 亚洲区| 观看美女的网站| 色婷婷久久久亚洲欧美| 韩国高清视频一区二区三区| 伊人久久国产一区二区| 精品人妻在线不人妻| 99精国产麻豆久久婷婷| 丝袜人妻中文字幕| 久久这里只有精品19| 免费av中文字幕在线| 性色avwww在线观看| 在线免费观看不下载黄p国产| 国产xxxxx性猛交| 午夜福利,免费看| 99国产综合亚洲精品| 国产黄色免费在线视频| 七月丁香在线播放| √禁漫天堂资源中文www| 2021少妇久久久久久久久久久| 成人二区视频| 精品视频人人做人人爽| av.在线天堂| 日韩一本色道免费dvd| 九草在线视频观看| 天堂中文最新版在线下载| 亚洲国产av影院在线观看| 久久人人爽av亚洲精品天堂| 成年人午夜在线观看视频| 亚洲,欧美精品.| 纵有疾风起免费观看全集完整版| 亚洲欧洲国产日韩| 多毛熟女@视频| 春色校园在线视频观看| 中文字幕亚洲精品专区| 日韩欧美一区视频在线观看| av.在线天堂| 亚洲欧美成人综合另类久久久| www.精华液| 精品国产国语对白av| 午夜激情av网站| 国产色婷婷99| 熟女电影av网| 成人国产麻豆网| 天堂8中文在线网| 在现免费观看毛片| 老熟女久久久| 国产精品二区激情视频| 天天操日日干夜夜撸| 一级毛片 在线播放| 日本av手机在线免费观看| 三级国产精品片| 久久这里有精品视频免费| 亚洲成色77777| 亚洲美女黄色视频免费看| 啦啦啦在线观看免费高清www| 丝袜美足系列| 女人被躁到高潮嗷嗷叫费观| 18+在线观看网站| 亚洲av欧美aⅴ国产| 成人国产麻豆网| 亚洲国产欧美网| 国产乱人偷精品视频| 捣出白浆h1v1| 欧美人与性动交α欧美精品济南到 | 亚洲国产精品999| 大陆偷拍与自拍| 日韩中字成人| 久久综合国产亚洲精品| 亚洲五月色婷婷综合| 男女无遮挡免费网站观看| 免费观看在线日韩| av免费在线看不卡| 嫩草影院入口| 久久久久久人人人人人| 91精品三级在线观看| 欧美亚洲 丝袜 人妻 在线| 欧美另类一区| 国产片特级美女逼逼视频| 在线观看一区二区三区激情| 女性生殖器流出的白浆| 久久久久精品久久久久真实原创| 国产免费视频播放在线视频| 欧美 日韩 精品 国产| 人人妻人人爽人人添夜夜欢视频| 亚洲av中文av极速乱| www.精华液| 日韩中文字幕欧美一区二区 | av卡一久久| 韩国av在线不卡| 一区福利在线观看| 欧美精品人与动牲交sv欧美| 久久久久视频综合| 99久国产av精品国产电影| 亚洲精品一区蜜桃| 精品国产一区二区三区久久久樱花| 久久这里只有精品19| 丝袜美足系列| 亚洲国产精品国产精品| 国产精品av久久久久免费| 国产精品秋霞免费鲁丝片| 高清不卡的av网站| 在线亚洲精品国产二区图片欧美| 亚洲国产精品999| 欧美精品一区二区免费开放| 在线观看美女被高潮喷水网站| 三级国产精品片| 亚洲久久久国产精品| 最黄视频免费看| 国产av国产精品国产| 天天躁夜夜躁狠狠躁躁| 亚洲国产日韩一区二区| 26uuu在线亚洲综合色| av卡一久久| 亚洲精品一区蜜桃| 丝袜美腿诱惑在线| av片东京热男人的天堂| 91精品伊人久久大香线蕉| 亚洲美女搞黄在线观看| 欧美最新免费一区二区三区| 国产极品粉嫩免费观看在线| 91在线精品国自产拍蜜月| 亚洲精品国产av蜜桃| 十八禁网站网址无遮挡| 最黄视频免费看| 一区二区三区四区激情视频| 久久久久久人妻| 久久女婷五月综合色啪小说| 成人国产麻豆网| 大片免费播放器 马上看| 乱人伦中国视频| 国产精品.久久久| 少妇人妻 视频| 色播在线永久视频| 男女午夜视频在线观看| 欧美日韩综合久久久久久| 国产成人精品婷婷| av女优亚洲男人天堂| 人妻 亚洲 视频| 国产精品 国内视频| 欧美精品av麻豆av| 国产精品一二三区在线看| 亚洲精品美女久久久久99蜜臀 | 纯流量卡能插随身wifi吗| 黄色配什么色好看| 久久久久精品久久久久真实原创| 好男人视频免费观看在线| 国产野战对白在线观看| 成人亚洲精品一区在线观看| 亚洲av男天堂| 亚洲欧美一区二区三区久久| 精品国产露脸久久av麻豆| 青草久久国产| 色视频在线一区二区三区| 如何舔出高潮| 亚洲中文av在线| 精品一区二区免费观看| 久久午夜福利片| 少妇精品久久久久久久| 久久99一区二区三区| 男女边吃奶边做爰视频| 下体分泌物呈黄色| 久久毛片免费看一区二区三区| 欧美日韩成人在线一区二区| 国产一区二区 视频在线| 三上悠亚av全集在线观看| 90打野战视频偷拍视频| 中国三级夫妇交换| 丝袜美足系列| 久久久久久久亚洲中文字幕| 最近最新中文字幕大全免费视频 | 成人午夜精彩视频在线观看| 成人毛片a级毛片在线播放| 色哟哟·www| 国产精品国产三级国产专区5o| 午夜免费男女啪啪视频观看| 在线观看免费日韩欧美大片| 啦啦啦中文免费视频观看日本| 女人高潮潮喷娇喘18禁视频| 午夜福利网站1000一区二区三区| 日韩av免费高清视频| 2021少妇久久久久久久久久久| 午夜福利视频在线观看免费| 男女下面插进去视频免费观看| 欧美另类一区| 超碰97精品在线观看| 欧美激情极品国产一区二区三区| 亚洲美女视频黄频| 王馨瑶露胸无遮挡在线观看| av线在线观看网站| 日本免费在线观看一区| 国产片内射在线| 国产片特级美女逼逼视频| 80岁老熟妇乱子伦牲交| 街头女战士在线观看网站| 免费观看性生交大片5| 丝袜在线中文字幕| 99re6热这里在线精品视频| 26uuu在线亚洲综合色| 精品国产国语对白av| 国产精品久久久av美女十八| 1024视频免费在线观看| 青青草视频在线视频观看| 女的被弄到高潮叫床怎么办| 五月开心婷婷网| 色吧在线观看| 黄色视频在线播放观看不卡| 国产伦理片在线播放av一区| 香蕉国产在线看| 亚洲精品av麻豆狂野| 高清在线视频一区二区三区| 国产免费福利视频在线观看| 黑人猛操日本美女一级片| 超碰97精品在线观看| 免费高清在线观看日韩| 女人精品久久久久毛片| 国产黄色视频一区二区在线观看| 久久久欧美国产精品| 日本色播在线视频| 午夜日韩欧美国产| 美女午夜性视频免费| 男女边摸边吃奶| 精品国产国语对白av| 成年女人在线观看亚洲视频| 国产 一区精品| 久久久久网色| 欧美黄色片欧美黄色片| 少妇人妻久久综合中文| 亚洲欧美一区二区三区国产| 国产一区二区三区综合在线观看| 一级毛片 在线播放| 最近的中文字幕免费完整| 国产精品二区激情视频| 香蕉国产在线看| 精品久久蜜臀av无| 国产黄色视频一区二区在线观看| 欧美精品高潮呻吟av久久| 日韩av不卡免费在线播放| 99九九在线精品视频| av视频免费观看在线观看| 少妇人妻 视频| 国产成人aa在线观看| 亚洲中文av在线| 97在线视频观看| 自线自在国产av| 亚洲精品视频女| 久久韩国三级中文字幕| 亚洲欧美色中文字幕在线| 婷婷色麻豆天堂久久| 啦啦啦在线免费观看视频4| 十八禁高潮呻吟视频| 日韩人妻精品一区2区三区| 国产成人a∨麻豆精品| 91成人精品电影| 国产极品天堂在线| 亚洲综合精品二区| 宅男免费午夜| 中国三级夫妇交换| 激情视频va一区二区三区| 国产福利在线免费观看视频| 国产av国产精品国产| 秋霞伦理黄片| 亚洲av成人精品一二三区| 欧美精品亚洲一区二区| 十八禁网站网址无遮挡| 日本vs欧美在线观看视频| 国产黄色视频一区二区在线观看| www.av在线官网国产| 成人影院久久| 亚洲国产精品国产精品| 婷婷色麻豆天堂久久| 美国免费a级毛片| 亚洲精品一二三| 亚洲精品乱久久久久久| 日韩一区二区三区影片| 在线天堂最新版资源| 日本av手机在线免费观看| 久久狼人影院| 最新中文字幕久久久久| 少妇 在线观看| 国产熟女欧美一区二区| 看免费av毛片| 日韩视频在线欧美| 看非洲黑人一级黄片| 欧美日韩精品网址| 丝袜美腿诱惑在线| 一区二区三区四区激情视频| 亚洲第一av免费看| 2022亚洲国产成人精品| 国语对白做爰xxxⅹ性视频网站| 国产成人一区二区在线| 大码成人一级视频| 国产精品久久久久久久久免| av网站免费在线观看视频| 日韩欧美精品免费久久| 欧美变态另类bdsm刘玥| 亚洲av成人精品一二三区| 日韩 亚洲 欧美在线| 国产又色又爽无遮挡免| 日韩一卡2卡3卡4卡2021年| 少妇熟女欧美另类| 欧美日韩一区二区视频在线观看视频在线| 黄色怎么调成土黄色| 午夜福利乱码中文字幕| 视频在线观看一区二区三区| 天天躁日日躁夜夜躁夜夜| 日韩一区二区视频免费看| 国产精品秋霞免费鲁丝片| 少妇 在线观看| 久久人妻熟女aⅴ| 日日摸夜夜添夜夜爱| 一区福利在线观看| 亚洲av.av天堂| 三上悠亚av全集在线观看| 最新的欧美精品一区二区| 婷婷色综合www| 亚洲中文av在线| 丰满乱子伦码专区| 国产日韩欧美亚洲二区| 两性夫妻黄色片| 久久久久国产精品人妻一区二区| 黑丝袜美女国产一区| 青青草视频在线视频观看| 午夜福利乱码中文字幕| 国产97色在线日韩免费| 啦啦啦啦在线视频资源| 电影成人av| 日本欧美国产在线视频| 九草在线视频观看| 人人澡人人妻人| 中文字幕亚洲精品专区| 麻豆av在线久日| 又黄又粗又硬又大视频| 波野结衣二区三区在线| 久久久久久久精品精品| kizo精华| 又粗又硬又长又爽又黄的视频| 1024视频免费在线观看| 卡戴珊不雅视频在线播放| 在线观看免费视频网站a站| 最新中文字幕久久久久| 国产 精品1| 国产精品一区二区在线观看99| 日韩成人av中文字幕在线观看| 日韩欧美一区视频在线观看| 十分钟在线观看高清视频www| 久久久久久久国产电影| 只有这里有精品99| 下体分泌物呈黄色| 肉色欧美久久久久久久蜜桃| 亚洲 欧美一区二区三区| 国产精品国产av在线观看| www.精华液| 国产一区有黄有色的免费视频| 少妇人妻 视频| 欧美成人精品欧美一级黄| 丝袜人妻中文字幕| 最近的中文字幕免费完整| 国产福利在线免费观看视频| 国产av国产精品国产| 97在线人人人人妻| av免费观看日本| 国产一级毛片在线| 黄片无遮挡物在线观看| 欧美日韩精品网址| 高清黄色对白视频在线免费看| 老司机影院成人| 国产成人精品久久二区二区91 | a级毛片在线看网站| 99久国产av精品国产电影| 亚洲av电影在线进入| 看免费av毛片| 日韩在线高清观看一区二区三区| 一级片'在线观看视频| 久久久a久久爽久久v久久| 国产国语露脸激情在线看| 18禁动态无遮挡网站| 制服诱惑二区| 久久久精品区二区三区| 国产黄色免费在线视频| 国产精品免费视频内射| 亚洲第一区二区三区不卡| 亚洲情色 制服丝袜| 久久久国产欧美日韩av| 人体艺术视频欧美日本| 久久 成人 亚洲| 日韩,欧美,国产一区二区三区| 亚洲情色 制服丝袜| 国产又色又爽无遮挡免| 伦理电影免费视频| 精品人妻一区二区三区麻豆| 国产黄频视频在线观看| 亚洲欧美一区二区三区久久| 人妻 亚洲 视频| 中文字幕人妻熟女乱码| 丰满少妇做爰视频| 国产亚洲av片在线观看秒播厂| 国产有黄有色有爽视频| 日韩大片免费观看网站| av卡一久久| 国产精品久久久av美女十八| 在线看a的网站| 精品第一国产精品| 秋霞伦理黄片| xxx大片免费视频| 亚洲欧美色中文字幕在线| 成人国产av品久久久| 大香蕉久久成人网| 秋霞在线观看毛片| 丝袜人妻中文字幕| 亚洲三区欧美一区| 亚洲精品自拍成人| 亚洲欧美日韩另类电影网站| 日韩一本色道免费dvd| 最近2019中文字幕mv第一页| 久久久亚洲精品成人影院| 咕卡用的链子| 国产成人精品福利久久| xxxhd国产人妻xxx| 亚洲国产精品999| 午夜免费鲁丝| 在线天堂中文资源库| 久久国产亚洲av麻豆专区| 五月开心婷婷网| 美国免费a级毛片| 午夜福利在线观看免费完整高清在| 国产精品熟女久久久久浪| 国产福利在线免费观看视频| 成年人午夜在线观看视频| 亚洲av欧美aⅴ国产| 亚洲精品成人av观看孕妇| 久久精品亚洲av国产电影网| 一区二区三区激情视频| 国产 一区精品| 亚洲国产成人一精品久久久| 麻豆av在线久日| 国产精品一二三区在线看| 国产一区二区三区av在线| 女的被弄到高潮叫床怎么办| 久久精品久久久久久噜噜老黄| 丝袜在线中文字幕| 国产国语露脸激情在线看| 极品少妇高潮喷水抽搐| 日韩三级伦理在线观看| 亚洲第一av免费看| 最近最新中文字幕免费大全7| 天天影视国产精品| 2018国产大陆天天弄谢| 国产成人午夜福利电影在线观看| 青春草亚洲视频在线观看| 亚洲欧美成人精品一区二区| 天堂8中文在线网| 性高湖久久久久久久久免费观看| 大片免费播放器 马上看| 欧美日韩一区二区视频在线观看视频在线| 曰老女人黄片| 只有这里有精品99| 999精品在线视频| 亚洲国产av影院在线观看| 日本91视频免费播放| 2018国产大陆天天弄谢| 1024香蕉在线观看| 在线观看美女被高潮喷水网站| 伦理电影大哥的女人| 美国免费a级毛片| 亚洲一码二码三码区别大吗| 一二三四在线观看免费中文在| 亚洲精品,欧美精品| 亚洲精品中文字幕在线视频| 男人操女人黄网站| 自线自在国产av| 日韩三级伦理在线观看| 成人午夜精彩视频在线观看| 欧美日韩视频高清一区二区三区二| 精品酒店卫生间| 最近手机中文字幕大全| 老鸭窝网址在线观看| 青春草国产在线视频| 男人操女人黄网站| 亚洲国产精品国产精品| 大码成人一级视频| 欧美变态另类bdsm刘玥| 综合色丁香网| 99精国产麻豆久久婷婷| 日本欧美国产在线视频|