• <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?
    av卡一久久| 国产精品一区二区在线观看99 | 尤物成人国产欧美一区二区三区| 日本撒尿小便嘘嘘汇集6| 亚洲av第一区精品v没综合| 噜噜噜噜噜久久久久久91| 婷婷色综合大香蕉| 欧美精品一区二区大全| 少妇人妻一区二区三区视频| 美女国产视频在线观看| 亚洲18禁久久av| av在线天堂中文字幕| 欧美性感艳星| 卡戴珊不雅视频在线播放| 免费搜索国产男女视频| 亚洲婷婷狠狠爱综合网| 国产黄色小视频在线观看| 国产高清有码在线观看视频| 亚洲va在线va天堂va国产| 国产精品久久久久久精品电影小说 | 免费在线观看成人毛片| 国产69精品久久久久777片| 成人二区视频| 精品日产1卡2卡| 精品久久久噜噜| 国产伦一二天堂av在线观看| 欧美日韩在线观看h| 成人特级黄色片久久久久久久| 日本黄大片高清| 22中文网久久字幕| 国产乱人视频| 五月伊人婷婷丁香| 日韩精品有码人妻一区| 久久久久久久久久黄片| 日本免费a在线| 嫩草影院新地址| 高清毛片免费观看视频网站| 中文字幕精品亚洲无线码一区| 99热全是精品| 中出人妻视频一区二区| 国产亚洲av片在线观看秒播厂 | 一卡2卡三卡四卡精品乱码亚洲| 精品人妻一区二区三区麻豆| 久久久久久久久中文| 久久久久国产网址| 伦理电影大哥的女人| 欧美不卡视频在线免费观看| 国内少妇人妻偷人精品xxx网站| 26uuu在线亚洲综合色| 国产精品国产三级国产av玫瑰| 国产黄a三级三级三级人| 精品久久久久久久久av| 久久久久久大精品| 卡戴珊不雅视频在线播放| 1024手机看黄色片| 久久久国产成人免费| 欧美色欧美亚洲另类二区| 最近最新中文字幕大全电影3| 中文在线观看免费www的网站| 99久久人妻综合| 特大巨黑吊av在线直播| 国产视频内射| 欧美另类亚洲清纯唯美| АⅤ资源中文在线天堂| 国内少妇人妻偷人精品xxx网站| 日本免费一区二区三区高清不卡| 亚洲精品乱码久久久久久按摩| 国内揄拍国产精品人妻在线| 久久久久久国产a免费观看| 精品久久久久久久末码| 少妇熟女欧美另类| 最近最新中文字幕大全电影3| 久久久精品大字幕| 国产伦一二天堂av在线观看| 天天躁日日操中文字幕| a级毛片a级免费在线| 18禁在线无遮挡免费观看视频| 国产麻豆成人av免费视频| 国产av一区在线观看免费| 国产午夜精品一二区理论片| 女同久久另类99精品国产91| 精品人妻视频免费看| 国产不卡一卡二| 九九热线精品视视频播放| 一级毛片我不卡| 中文字幕精品亚洲无线码一区| 久99久视频精品免费| 国产精品蜜桃在线观看 | 午夜爱爱视频在线播放| av在线老鸭窝| 国语自产精品视频在线第100页| 人妻夜夜爽99麻豆av| 久久99蜜桃精品久久| 日韩av不卡免费在线播放| 久久精品国产亚洲av香蕉五月| 精品久久国产蜜桃| av视频在线观看入口| av在线播放精品| 亚洲av第一区精品v没综合| 一级毛片电影观看 | av又黄又爽大尺度在线免费看 | 亚洲欧美日韩无卡精品| 午夜a级毛片| 精品一区二区三区人妻视频| 色噜噜av男人的天堂激情| 亚洲国产欧美在线一区| 99热网站在线观看| 久久精品国产自在天天线| 久久热精品热| 12—13女人毛片做爰片一| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 国产探花极品一区二区| 日本成人三级电影网站| 精品不卡国产一区二区三区| 国产极品精品免费视频能看的| 亚洲经典国产精华液单| 免费看光身美女| 国产一级毛片七仙女欲春2| 亚洲欧美精品专区久久| 最新中文字幕久久久久| 亚洲综合色惰| 久久国产乱子免费精品| 国产午夜精品久久久久久一区二区三区| 日韩欧美精品v在线| eeuss影院久久| 精品午夜福利在线看| 国产一区二区激情短视频| 少妇裸体淫交视频免费看高清| 国产成人影院久久av| 桃色一区二区三区在线观看| 秋霞在线观看毛片| 精品日产1卡2卡| 精品久久久噜噜| 国产成人午夜福利电影在线观看| 成人一区二区视频在线观看| 九草在线视频观看| 国产白丝娇喘喷水9色精品| 日韩成人伦理影院| 久久久久国产网址| 2022亚洲国产成人精品| 国产精品国产三级国产av玫瑰| 亚洲国产欧美人成| 精品久久国产蜜桃| 亚洲av一区综合| 两个人视频免费观看高清| 国内精品宾馆在线| 国产69精品久久久久777片| 少妇高潮的动态图| 国产一区二区亚洲精品在线观看| 亚洲欧美精品综合久久99| 在线a可以看的网站| 亚洲av成人av| 欧美丝袜亚洲另类| 免费无遮挡裸体视频| 在线免费观看的www视频| 国产淫片久久久久久久久| 国产精品人妻久久久久久| 两个人的视频大全免费| 日韩精品有码人妻一区| 欧美日韩精品成人综合77777| 少妇被粗大猛烈的视频| 亚洲欧美日韩无卡精品| 26uuu在线亚洲综合色| 日日撸夜夜添| 激情 狠狠 欧美| 欧美成人a在线观看| 国产人妻一区二区三区在| 亚洲三级黄色毛片| 最新中文字幕久久久久| 日本一二三区视频观看| 欧美人与善性xxx| 国产成人91sexporn| 亚洲av免费高清在线观看| 日韩在线高清观看一区二区三区| 久久精品国产亚洲av香蕉五月| 麻豆成人午夜福利视频| 国语自产精品视频在线第100页| 欧美精品国产亚洲| 男人狂女人下面高潮的视频| 欧美日韩一区二区视频在线观看视频在线 | 中文在线观看免费www的网站| 久久热精品热| 免费搜索国产男女视频| 亚洲欧美成人精品一区二区| 在线免费观看的www视频| 非洲黑人性xxxx精品又粗又长| 欧美不卡视频在线免费观看| 亚洲最大成人av| 午夜福利成人在线免费观看| 亚洲精品亚洲一区二区| 18禁裸乳无遮挡免费网站照片| 国产伦精品一区二区三区四那| 亚洲成人中文字幕在线播放| 简卡轻食公司| 美女脱内裤让男人舔精品视频 | 国产色爽女视频免费观看| 午夜福利在线在线| 一本久久中文字幕| 欧美激情久久久久久爽电影| 长腿黑丝高跟| 日本在线视频免费播放| 男人舔奶头视频| av卡一久久| 免费观看精品视频网站| 日本五十路高清| 深爱激情五月婷婷| 色播亚洲综合网| 一进一出抽搐动态| 日本在线视频免费播放| 日韩成人伦理影院| 你懂的网址亚洲精品在线观看 | 色综合亚洲欧美另类图片| 亚洲最大成人av| 欧美色欧美亚洲另类二区| 内射极品少妇av片p| 久久久久久久久久黄片| 亚洲av电影不卡..在线观看| 日韩精品青青久久久久久| 天天躁日日操中文字幕| 内射极品少妇av片p| 精品午夜福利在线看| 免费看日本二区| 免费观看在线日韩| 久久韩国三级中文字幕| 看黄色毛片网站| 国产高清激情床上av| 少妇人妻精品综合一区二区 | 国产视频首页在线观看| 中文字幕人妻熟人妻熟丝袜美| 最好的美女福利视频网| 国产成人影院久久av| 亚洲性久久影院| 日韩欧美三级三区| 99热这里只有精品一区| 欧美日韩在线观看h| 国产精品嫩草影院av在线观看| 亚洲在久久综合| 长腿黑丝高跟| 最近最新中文字幕大全电影3| 国产真实伦视频高清在线观看| 禁无遮挡网站| 黄色日韩在线| 免费看美女性在线毛片视频| 免费一级毛片在线播放高清视频| 美女黄网站色视频| 直男gayav资源| 只有这里有精品99| 简卡轻食公司| 插逼视频在线观看| 美女内射精品一级片tv| 亚洲在久久综合| 国产一区二区在线av高清观看| 成年av动漫网址| 蜜臀久久99精品久久宅男| 日韩欧美精品免费久久| 三级经典国产精品| 天美传媒精品一区二区| 一级毛片aaaaaa免费看小| 久久99热这里只有精品18| 99热这里只有是精品在线观看| 亚洲欧美精品自产自拍| 国产成人影院久久av| 国产成人freesex在线| 成年女人永久免费观看视频| 中文在线观看免费www的网站| 久久国产乱子免费精品| 麻豆av噜噜一区二区三区| 99久久久亚洲精品蜜臀av| 性插视频无遮挡在线免费观看| 国产精品人妻久久久久久| 黑人高潮一二区| 晚上一个人看的免费电影| 日韩人妻高清精品专区| 精品欧美国产一区二区三| 色综合色国产| 最近中文字幕高清免费大全6| 国产精品电影一区二区三区| 久久久久久久久久黄片| 免费黄网站久久成人精品| 亚洲真实伦在线观看| 国产精品嫩草影院av在线观看| 久久久久久久久久久免费av| 最近中文字幕高清免费大全6| 麻豆一二三区av精品| 免费电影在线观看免费观看| 尤物成人国产欧美一区二区三区| 亚洲欧美成人综合另类久久久 | 亚洲在久久综合| 97在线视频观看| 十八禁国产超污无遮挡网站| 色综合亚洲欧美另类图片| 蜜桃亚洲精品一区二区三区| 自拍偷自拍亚洲精品老妇| 国产真实乱freesex| 九九在线视频观看精品| 亚洲精品影视一区二区三区av| 国产成人一区二区在线| 18禁在线播放成人免费| 日日摸夜夜添夜夜添av毛片| 亚洲最大成人av| 啦啦啦韩国在线观看视频| 久久久成人免费电影| 国产精品伦人一区二区| 蜜臀久久99精品久久宅男| 亚洲国产色片| 插逼视频在线观看| 午夜免费激情av| 国产黄片美女视频| 亚洲欧美日韩卡通动漫| 国产成人freesex在线| 中国国产av一级| 99久久精品一区二区三区| 免费电影在线观看免费观看| 99热这里只有精品一区| av免费观看日本| 亚洲成a人片在线一区二区| 边亲边吃奶的免费视频| 午夜福利视频1000在线观看| 在线天堂最新版资源| 亚洲七黄色美女视频| 国产日韩欧美在线精品| 91久久精品电影网| 一边亲一边摸免费视频| 麻豆国产av国片精品| 国产色婷婷99| 99视频精品全部免费 在线| 看非洲黑人一级黄片| 欧美zozozo另类| 亚洲av电影不卡..在线观看| 精品不卡国产一区二区三区| 精品久久久久久久久av| 国产69精品久久久久777片| 一级黄色大片毛片| 亚洲精品久久国产高清桃花| 国产熟女欧美一区二区| 插阴视频在线观看视频| 欧美区成人在线视频| 嫩草影院入口| 国语自产精品视频在线第100页| 一区二区三区四区激情视频 | 亚洲精品日韩av片在线观看| 久久精品夜色国产| 亚洲精品久久久久久婷婷小说 | 久久草成人影院| 2021天堂中文幕一二区在线观| 精品无人区乱码1区二区| 日韩一本色道免费dvd| 成人av在线播放网站| 不卡一级毛片| 简卡轻食公司| 国产三级中文精品| 99热全是精品| 亚洲乱码一区二区免费版| 麻豆精品久久久久久蜜桃| 日韩一本色道免费dvd| 久99久视频精品免费| 色哟哟哟哟哟哟| 99久国产av精品国产电影| 国产精品久久久久久久电影| 国产亚洲精品久久久com| 欧美日韩一区二区视频在线观看视频在线 | 精品日产1卡2卡| 亚洲成人中文字幕在线播放| 亚洲久久久久久中文字幕| 少妇丰满av| 极品教师在线视频| 久久精品久久久久久噜噜老黄 | 色综合色国产| 中国美女看黄片| av福利片在线观看| 美女cb高潮喷水在线观看| 国产熟女欧美一区二区| 亚洲欧美日韩高清专用| 男女啪啪激烈高潮av片| 国产亚洲欧美98| 日产精品乱码卡一卡2卡三| 午夜激情福利司机影院| 村上凉子中文字幕在线| 国产精品伦人一区二区| 久久精品夜夜夜夜夜久久蜜豆| 99热这里只有精品一区| 久久久久久伊人网av| 在线播放国产精品三级| 91麻豆精品激情在线观看国产| 亚洲第一电影网av| 亚洲精品色激情综合| 黄色欧美视频在线观看| a级毛片免费高清观看在线播放| 看黄色毛片网站| 激情 狠狠 欧美| 亚洲无线观看免费| 嫩草影院新地址| 黄色一级大片看看| 能在线免费看毛片的网站| 熟女电影av网| 亚洲欧美日韩高清专用| 欧美+亚洲+日韩+国产| 最近的中文字幕免费完整| 男人舔奶头视频| 99热这里只有是精品在线观看| 国产午夜精品久久久久久一区二区三区| 免费看av在线观看网站| 国内揄拍国产精品人妻在线| 老女人水多毛片| 国产高清视频在线观看网站| 小蜜桃在线观看免费完整版高清| 青青草视频在线视频观看| 欧美人与善性xxx| 日日干狠狠操夜夜爽| 国产白丝娇喘喷水9色精品| 一边摸一边抽搐一进一小说| 国内精品美女久久久久久| 久99久视频精品免费| 日日啪夜夜撸| 嘟嘟电影网在线观看| 午夜a级毛片| 又爽又黄a免费视频| 免费不卡的大黄色大毛片视频在线观看 | 性色avwww在线观看| 免费观看在线日韩| 一进一出抽搐动态| 成人三级黄色视频| 高清毛片免费看| 我的女老师完整版在线观看| 国产精品电影一区二区三区| 18禁在线播放成人免费| 看十八女毛片水多多多| 久久九九热精品免费| 国产爱豆传媒在线观看| 免费看av在线观看网站| 男女做爰动态图高潮gif福利片| 国产 一区 欧美 日韩| 欧美精品国产亚洲| 国产亚洲精品av在线| 亚洲七黄色美女视频| av天堂中文字幕网| 欧美色视频一区免费| 12—13女人毛片做爰片一| 非洲黑人性xxxx精品又粗又长| 亚洲乱码一区二区免费版| 国产不卡一卡二| 欧美成人免费av一区二区三区| 精华霜和精华液先用哪个| 欧美另类亚洲清纯唯美| 黑人高潮一二区| 99精品在免费线老司机午夜| 久久久精品欧美日韩精品| 欧美日本视频| 天美传媒精品一区二区| 波野结衣二区三区在线| 日本欧美国产在线视频| 亚洲真实伦在线观看| 人妻久久中文字幕网| 精品一区二区免费观看| 中文在线观看免费www的网站| av又黄又爽大尺度在线免费看 | 日本一本二区三区精品| 国产黄a三级三级三级人| 亚洲无线观看免费| 美女内射精品一级片tv| 91久久精品电影网| 免费看日本二区| www.色视频.com| 久久久精品大字幕| 偷拍熟女少妇极品色| 亚洲色图av天堂| 国产亚洲av嫩草精品影院| 又粗又爽又猛毛片免费看| 国产一级毛片在线| 美女国产视频在线观看| 赤兔流量卡办理| 91在线精品国自产拍蜜月| av福利片在线观看| 熟女人妻精品中文字幕| 成人三级黄色视频| 久久久a久久爽久久v久久| 国产一区二区激情短视频| 欧美激情国产日韩精品一区| 插阴视频在线观看视频| 欧美又色又爽又黄视频| 亚洲精品乱码久久久v下载方式| 久久久久久久久久黄片| 中国美白少妇内射xxxbb| 国产精品久久久久久av不卡| 欧美精品一区二区大全| 国产精品蜜桃在线观看 | 亚洲欧美精品综合久久99| 99久久中文字幕三级久久日本| 日韩一区二区视频免费看| 一夜夜www| 插逼视频在线观看| 少妇熟女欧美另类| 菩萨蛮人人尽说江南好唐韦庄 | 日本熟妇午夜| 美女cb高潮喷水在线观看| 波多野结衣高清无吗| 亚洲久久久久久中文字幕| 国产亚洲91精品色在线| 一个人看视频在线观看www免费| 国产成人影院久久av| 此物有八面人人有两片| 日本在线视频免费播放| 久久久精品94久久精品| 看十八女毛片水多多多| 熟女人妻精品中文字幕| www日本黄色视频网| 中文字幕人妻熟人妻熟丝袜美| 欧美一区二区国产精品久久精品| 国产免费男女视频| 18禁裸乳无遮挡免费网站照片| 亚洲激情五月婷婷啪啪| 国产91av在线免费观看| 久久精品影院6| 伊人久久精品亚洲午夜| 精品欧美国产一区二区三| 免费看a级黄色片| 亚洲国产精品久久男人天堂| 精品久久国产蜜桃| 免费无遮挡裸体视频| 国产淫片久久久久久久久| 97人妻精品一区二区三区麻豆| 中文字幕人妻熟人妻熟丝袜美| 亚洲欧洲日产国产| 黄色一级大片看看| 久99久视频精品免费| 国内精品久久久久精免费| 丰满乱子伦码专区| 久久久久久久久久久免费av| 狂野欧美白嫩少妇大欣赏| 国产不卡一卡二| 一级毛片aaaaaa免费看小| 亚洲av熟女| 寂寞人妻少妇视频99o| 日本黄色片子视频| 国产精品国产三级国产av玫瑰| 日日摸夜夜添夜夜添av毛片| 99riav亚洲国产免费| 久久欧美精品欧美久久欧美| 又爽又黄a免费视频| 久久久精品欧美日韩精品| 国产探花极品一区二区| 99在线视频只有这里精品首页| 亚洲精品成人久久久久久| 精品久久久噜噜| 国产黄色小视频在线观看| 亚洲第一区二区三区不卡| 欧美3d第一页| av免费在线看不卡| 国产午夜精品论理片| 美女cb高潮喷水在线观看| 亚洲欧美日韩无卡精品| 中文字幕久久专区| 色播亚洲综合网| 人妻夜夜爽99麻豆av| 久久草成人影院| 99在线人妻在线中文字幕| 伊人久久精品亚洲午夜| 国产真实乱freesex| 亚洲欧洲国产日韩| 精品久久久久久久人妻蜜臀av| 久久午夜福利片| 波多野结衣高清无吗| 国产综合懂色| 国产精品女同一区二区软件| 日本黄大片高清| 成人亚洲精品av一区二区| 精品久久久久久久久av| 国产精品1区2区在线观看.| 日韩一区二区三区影片| 日韩三级伦理在线观看| av在线亚洲专区| 日韩高清综合在线| 亚洲乱码一区二区免费版| 啦啦啦啦在线视频资源| 亚洲成人久久爱视频| 免费大片18禁| 国产色婷婷99| 偷拍熟女少妇极品色| 夫妻性生交免费视频一级片| 搡老妇女老女人老熟妇| 麻豆成人av视频| 亚洲精品国产av成人精品| 日本-黄色视频高清免费观看| 狂野欧美白嫩少妇大欣赏| 午夜a级毛片| 亚洲欧美日韩高清在线视频| 搡老妇女老女人老熟妇| 国产免费一级a男人的天堂| 久久精品国产99精品国产亚洲性色| 99久久中文字幕三级久久日本| 亚洲精华国产精华液的使用体验 | 九九久久精品国产亚洲av麻豆| 天堂网av新在线| 观看美女的网站| 午夜福利在线观看免费完整高清在 | 中文字幕人妻熟人妻熟丝袜美| 听说在线观看完整版免费高清| 精品不卡国产一区二区三区| 麻豆精品久久久久久蜜桃| 97超视频在线观看视频| 中文字幕av在线有码专区| 国模一区二区三区四区视频| 久久99蜜桃精品久久| 黄色欧美视频在线观看| 卡戴珊不雅视频在线播放| 亚洲精品亚洲一区二区| 老师上课跳d突然被开到最大视频| 欧美日韩乱码在线| 国产女主播在线喷水免费视频网站 | 在线观看美女被高潮喷水网站| 中文亚洲av片在线观看爽| 久久久国产成人精品二区| 天堂√8在线中文| 国产精品1区2区在线观看.| 99热这里只有是精品在线观看|