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

    Theoretical analyses on the one-dimensional charged particle transport in a decaying plasma under an electrostatic field

    2023-10-11 07:55:44YaoTingWang汪耀庭XinLiSun孫鑫禮LanYueLuo羅嵐月ZiMingZhang張子明HePingLi李和平DongJunJiang姜東君andMingShengZhou周明勝
    Chinese Physics B 2023年9期
    關(guān)鍵詞:東君和平

    Yao-Ting Wang(汪耀庭), Xin-Li Sun(孫鑫禮), Lan-Yue Luo(羅嵐月), Zi-Ming Zhang(張子明),He-Ping Li(李和平), Dong-Jun Jiang(姜東君), and Ming-Sheng Zhou(周明勝)

    Department of Engineering Physics,Tsinghua University,Beijing 100084,China

    Keywords: transport of charged particles, decaying plasma, low-pressure plasma, theoretical analysis,

    1.Introduction

    Wall-confined decaying plasmas widely exist in the applications of plasma immersion ion implantation,[1–4]afterglow plasmas[5–7]and ion extraction process of laser-induced plasmas.[8–14]In these fields of engineering and academic research,ion extraction flux is one of the key parameters,which is determined by the evolution of sheath and the propagation of waves in plasmas.In this paper,we focus on the ion extraction process of the laser-induced plasmas, and the spatiotemporal evolutions of the decaying plasmas under an externally applied electrostatic field are analyzed theoretically.

    For a typical ion extraction process, as shown in Fig.1,the metal vapor is formed starting from the pool of the molten metal; then, a weakly-ionized metal plasma is induced by a laser pulse during a very short period; and then, an electrostatic field is imposed between the electrodes, which results in the extraction of the charged particles until all of them are collected on the electrode surfaces.Two parameters, i.e., the ion extraction flux (Γi), and ion extraction time(text) defined as the time required from the start of the ion extraction process to the time point when 90% of ions are extracted, are used to describe the features of the ion extraction process in previous studies.[12]Up to now, many efforts have been made to improveΓi, such as improvements of extracting electrode configurations (e.g., Π type[10,13,15,16]and M type[15,17–19]) and weakening of the shielding effect of high voltage sheath (e.g., RF resonance method[11,20–24]),etc.Analytical description,[7,9,25–27]fluid simulation,[28,29]and particle-in-cell (PIC) simulation[15,16,30,31]are the three major theoretical methods for studying the ion extraction process.The fluid method consumes fewer computing resources,but it provides less accurate descriptions of the whole process and is even not applicable to some cases.[16]Fewer assumptions are employed in the PIC method which gives a more accurate description compared with the fluid model, however,it consumes much more computing resources.For analytical description,it can provide a clear physical image of the ion extraction process provided that a reasonable simplification can be made.Thus, it is of great significance to establish an analytical model to analyze and predict the physical processes during the decaying of a laser-induced plasma.

    Physically, for the singly- or doubly-bounded decaying plasmas, the sheath dynamics can operate in a supersonic regime or a subsonic regime depending on the relative speed between the sheath front movement and the ion rarefaction wave(IRW)propagation.[30]For the sheath dynamics in a supersonic regime,the theoretical analysis on the ion extraction characteristics in a decaying plasma was firstly conducted by Chen.[9]Chen’s work shows that, during the whole ion extraction process, firstly the sheath will expand rapidly until the balance between the Langmuir flux and the Bohm flux is obtained; then, the density of the bulk plasma will drop slowly with almost unchanged sheath thickness; and finally,the plasma will collapse quickly when the density and the thickness of the bulk plasma reach sufficient small values.Based on Chen’s work, Okano[25]derived a simplified formula forΓithat is applicable over a wide initial plasma density range to evaluatetext.While for the sheath dynamics in a subsonic region, the IRW propagation is also a dominant factor contributing to the ion fluxes.Combining sheath expansion and IRW propagation, Murakami and Nishihara[26]developed an analytical model for describing the sheath expansion and IRW propagation in a decaying plasma bounded by a single negatively biased electrode.In addition,the whole process is generally divided into two phases,i.e.,the ion matrix sheath relaxation and the Langmuir sheath expansion,in the plasma immersion ion implantation (PIII) with a decaying single-wall confined plasma;[1,2]and correspondingly,the formulas for the ion flux or the governing equations for sheath dynamics in the whole implantation process have also been presented.Focusing on the ion extraction process in the decaying plasmas confined between two electrodes,recent studies showed that: the ion flux to the anode could not be ignored and was greatly influenced by the initial electron temperature of plasmas;[30]and the whole ion extraction process could be divided into three stages, i.e., the electron oscillation stage, sheath expansion and IRW propagation stage,and the subsequent rapid decaying stage.[16,30]However, up to now, the ion flux to the anode during the whole process and the fluxes of the charged particles during the electron oscillation stage are not included in the analytical model.Therefore,the goal of this paper is to establish a complete analytical model for describing the transport of the charged particles during the whole decaying process of the laser-induced plasma bounded between two electrodes under the externally applied electrostatic field.This paper is organized as follows:descriptions on the physical model and the PIC code are presented in Section 2, then, theoretical analyses and simulation validations are provided in Section 3, and then, discussions are presented in Section 4, and finally, the major conclusions are summarized in Section 5 with a brief discussion on future research.

    Fig.1.Schematic diagram of the ion extraction device.

    2.Descriptions on the physical model and the PIC code

    In this study, a one-dimensional(1D)ion extraction system is considered as shown in Fig.2.A laser-induced plasma consisting of argon ions(Ar+)and electrons(e)is suddenly generated att=0 and uniformly distributed between a grounded anode (?x= ?L) and a negatively biased cathode fixed at-?U0(?x=0).Once an electrostatic field is imposed,the ion extraction process starts at ?t >0.It is assumed that:(i) the plasma considered here is collisionless because the species mean free path is much larger than the characteristic length of the system due to the typical low chamber pressure(~10-5Torr,1 Torr=1.33322×102Pa); (ii)the particles are absorbed by the metal plates when they reach the electrodes,and the subsequent sputtering and secondary electron emission effects are neglected; (iii) the initial velocity distributions of the electrons and ions are Maxwellian corresponding to their initial temperatures, ?Te,0and ?Ti,0,respectively.

    Fig.2.Schematic diagram for studying the 1D transport of charged particles in a decaying plasma.

    In this study, a 1D code, EDIPIC (the electrostatic direct implicit particle-in-cell code),[32]is employed to simulate the evolutional properties of the bounded decaying plasmas.EDIPIC is a well-developed PIC code for studying the ion extraction process and its accuracy and applicability have been well-verified in previous studies.[30,31]In this paper,considering the needs for both accuracy and consumption of the computing resources,the spatial and time intervals are selected as follows: the grid size(Δ?x)is 0.1 times the Debye length(λD),that is, Δ?x=0.1λD; the time interval is determined to make sure that an electron will not travel across a grid during one time interval(Δ?t),that is ?veΔ?t <Δ?x,where ?veis the maximum electron speed considered in the PIC simulation; in particular, ?veequals 20 times of the electron thermal velocity in this study.The initial macroparticles in a cell are set to be 256,i.e.,Np=256.The operating parameters for a typical case studied in the following sections are listed in Table 1.Here,it should be emphasize that: since in previous PIC simulations of the ion extraction process,the variation ranges for ?Te,0, ?Ti,0,?L,and ?n0are 0.1 eV~1.0 eV,0.01 eV~0.1 eV,1.0 cm~3.0 cm,and(0.1~10)×1016m-3,respectively,[15,30,31,33,34]the values of ?Te,0=0.5 eV, ?Ti,0=0.02 eV, ?L=2.0 cm, and ?n0=1.0×1016m-3are used in this study;while for the externally applied voltage,the value of ?U0=300 V is used according to Ref.[35]for the extraction of argon ions.

    Table 1.Typical simulation parameters used in this paper.

    3.Theoretical analyses and numerical validations

    3.1.Normalization for the physical parameters and equations

    In this section, all parameters and the related governing equations will be normalized before we start the analyses.The specific normalization process is as follows(for simplicity,the variables with and without “^” represent the parameters before and after normalization):

    wheresrefers to the sheath thickness,xis the spatial variable,ωpiis the ion plasma frequency,urefers to the macroscopic ion velocity,csis the ion acoustic velocity,nrepresents the ion or electron number density,φrefers to the electric potential,eis the elementary charge,kis the Boltzmann constant, andLstands for the distance between electrodes.Thus,the normalized formula for the Langmuir flux and the Poisson equation are expressed as

    According to Refs.[26,36],the sheath motion function,the ion density,and velocity profiles in an IRW are as follows:

    wherex0andt0are the starting point of IRW in space and time.

    3.2.Theoretical analyses of the whole ion extraction process

    Based on the previous studies,[16,30]the whole ion extraction process can be divided into three stages, i.e., the electron oscillation and ion matrix sheath extraction stage,sheath expansion and ion rarefaction wave(IRW)propagation stage,and the plasma collapse stage.However,the specific criterion for indicating the starting and/or ending of each stage is still not clear.In this sub-section, the criteria for differentiating various stages during the decaying of the laser-induced plasmas under an externally applied electrostatic field,as shown in Fig.2,are derived theoretically;and then,the expressions for the ion extraction flux at each stage are obtained for predicting the characteristics of the ion extraction process.

    3.2.1.Electron oscillation and ion matrix sheath extraction stage

    When an external electric field is suddenly imposed on a plasma,the electrons will immediately respond to the electric field to form an electron-free region,which is the so-called ion matrix sheath.In the previous studies on PIII,[1,2]the thickness of the initial matrix near the cathode is thought to be equal to (2U0)1/2.However, when the plasma is confined by two electrodes, electrons will oscillate between electrodes as shown in Fig.3.Here,we define the ion matrix sheath as the largest electron-free region near the cathode and anode, and thus, the ion matrix sheath extraction process will occur near both electrodes.In Ref.[37], it was shown that the thickness of the initial ion matrix sheath near the cathode is

    where the subscripts“m”and“ca”represent the matrix sheath and cathode, respectively.Because the normalized electron motion equation can be expressed as

    the initial ion matrix sheath thickness near the anode is

    where the subscript“an”stands for the anode.

    In previous studies on the evolution of the ion matrix sheath,[1]the ions are assumed to be uniformly distributed in space, and then, the temporal evolution of the ion density in the matrix sheath is expressed as

    While for a stable Langmuir sheath, the spatial distribution of the ion density can be approximated as follows(a detailed derivation is presented in Appendix A):

    wherexrefers to the distance to the sheath edge.Since the ion matrix sheath extraction process is also the formation of the Langmuir sheath, it is reasonable to assume that the ion density in the vicinity of the cathode can be expressed by Eq.(9),and a time-dependent variableγ(t)can be introduced into Eq.(10) to describe the spatiotemporal evolution of the ion matrix sheath as

    Fig.3.Schematics of the species number density distributions during the electron oscillation process when the electrons move towards the anode(a)and the cathode(b),where ue represents the motion of electrons.

    After a coordinate transformation,the evolution of the ion matrix sheath near the cathode can be simplified as

    For the ion matrix sheath near the anode,during the first stage, the plasma potential will be higher than that of the anode.[16,30]Thus, the evolution of the matrix sheath will be similar to the situation near the cathode.According to Eq.(7),during the first stage of the ion extraction process, the range of the electron oscillation in the matrix sheath near the cathode and the anode are[sm,ca-sm,an,sm,ca]and[L-sm,an,L],respectively.Therefore, it is reasonable to believe that the ion density evolutions in the ion matrix sheath near the anode([L-sm,an,L])is similar to that in the region of[sm,ca-sm,an,sm,ca]near the cathode.Finally,the evolution of the ion matrix sheath near the anode can be expressed as

    In this subsection, comparisons are conducted between the theoretically predicted and PIC simulated results for the ion density distributions in the ion matrix sheath near the electrodes for the typical case listed in Table 1.It is seen from Fig.4 that a good agreement is obtained at different time points during the electron oscillation and ion matrix sheath extraction stage.Based on the preceding discussions,the criterion indicating the end of the first stage withγ=1 can be obtained from Eq.(11)as

    Correspondingly, the ion extraction flux can be obtained provided that the ion density evolution is known.Thus, the ion extraction fluxes near the cathode and anode during the first stage, i.e., the electron oscillation and ion matrix sheath extraction stage,can be expressed as

    in the above equations, the Hyp2F1 is the hypergeometric function.

    Fig.4.Comparisons of the theoretically predicted (blue line) and PIC simulation (red line) results of the ion density distributions in the ion matrix sheath both near the cathode and anode at t=0.47(a),0.91(b),2.00(c),and 3.09(d).

    3.2.2.Sheath expansion and IRW propagation stage

    Following the electron oscillation and ion matrix sheath extraction stage, the sheath expansion and the IRW propagation become the dominant phenomena during the ion extraction process.In the second stage, on the one hand, a moving Langmuir sheath is formed in the vicinity of the cathode,and moves towards the bulk plasma; on the other hand, the IRW starts to propagate from the initial matrix sheath edge near the anode, while the remained ions in the matrix sheath will be extracted until the IRW reaches the anode.For simplicity,the IRWs excited near the cathode and anode are noted as IRW-C and IRW-A,respectively,hereafter in this study.

    Firstly,the ion extraction flux on the cathode will be considered, as shown in Fig.5(a).According to Ref.[26], the propagating speed of the sheath edge towards the bulk plasma(vs= ds/dt) will gradually decrease with the increase of the sheath thickness.Therefore,compared tocs,the sheath expansion process may include two states:supersonic sheath expansion state whenvsexceedscs,and subsonic sheath expansion state whenvsis less thancs.For the supersonic sheath expansion state,since the sheath expansion speed is higher thancs,the IRW cannot be excited, and thus, makes no contributions to the sheath evolution.Therefore,the normalized sheath motion function during the supersonic state can be expressed as

    Since the supersonic sheath expansion state ends whenvs=1,the corresponding sheath thickness and the time point can be expressed as

    Here, it should be emphasized that, if the thickness of the initial ion matrix sheath exceedsssup,end,ca,there exists no supersonic sheath expansion state.For universality and the sake of discussion,we assume the existence of a supersonic sheath expansion process in the following analysis.

    During the subsonic state,the IRW will be excited and its front will overtake the sheath edge.Thus,there is a two-stage distribution for the ion density, as shown in Fig.5(b).In this case, the IRW will influence the ion density and speed at the sheath edge.Therefore, the normalized sheath motion equation is written as

    Fig.5.Schematics of the ion density distributions near the cathode during the supersonic(a)and subsonic(b)sheath expansion states.

    Secondly,in the vicinity of the anode,after the first stage ends, the sheath will gradually become a neutral region with the potential of the bulk plasma getting closer to that of the anode, as shown in Fig.6.Therefore, the IRW-A will propagate toward both the bulk plasma and the anode,and the corresponding time period when the IRW-A reaches the anode can be expressed as

    After the IRW reaches the anode,due to the properties of the IRW,the ion flux to the anode can be expressed as

    While before the IRW-A reaches the anode, the ions in the vanishing sheath will be slowly extracted.Considering the continuity of the ion fluxes,the ion extraction flux at this substage can be expressed as

    In addition, although the starting point of the IRW-A wheren=exp(-1) will move towards the anode after the front of IRW-A reaches the anode, the starting point of IRW-A will be set tox=Lfor simplicity in the following analysis since the distance between this starting point and the anode surface is relative small compared with the IRW travelling distance.Based on Eq.(19), it is shown that the IRW affects the ion speed and density at the sheath boundary;and thus,affects the sheath expansion process near the cathode.Equations (17)–(22) offer a comprehensive illustration of the ion extraction fluxes before the convergence of the two IRW fronts.Figure 7 gives a schematic description of the evolution of the ion density distributions before and after the IRW-A and IRW-C meet.It is seen from Fig.7 that, after the meet of the fronts of IRW-A and IRW-C,the two IRWs continue to propagate in their original directions as before.In the interaction region,the ion density and speed show parabolic and linear distributions,respectively.[7]During the whole second stage,i.e.,the sheath expansion and IRW propagation stage, the sheath advancing process near the cathode is only affected by the high voltage sheath and the IRW-C,while the extraction process on the anode is divided into two sub-stages.After the front of the IRW-A meets the sheath near the cathode, sheath properties near the cathode start to be affected by both IRW-A and IRWC, thus, the criterion for the beginning of the third stage and the corresponding sheath thickness at this time will be

    Fig.6.Schematic diagram of the ion density distributions near the anode after the electron oscillation and matrix sheath extraction stage ends.

    Fig.7.Schematic diagram of the ion density distributions before and after the IRW-A and IRW-C meet.

    3.2.3.Plasma collapse stage

    After the disappearance of the IRW-A,the ion density and speed at the sheath edge start to be affected by IRW-A.It is seen from Fig.8 that, in the third stage of the ion extraction process, i.e., the plasma collapse stage, the Langmuir sheath continues to expand towards the bulk plasma, while the ion density at the sheath edge drops quickly.The ion density in the interaction region typically exhibits a parabolic distribution; however, for the sake of simplicity, it can be approximated as being uniformly distributed,as discussed in detail in Appendix B.With the front of IRW-C as the boundary point,the ion velocity presents a piecewise linear distribution.

    Fig.8.Schematics of the ion density(a)and speed(b)distributions in the plasma collapse stage.

    For the sheath motion near the cathode, the sheath function is written as

    Thus,the ion density and velocity at the sheath edge can be expressed as

    From the analysis in Appendix C,the ion flux to the anode can be written as

    After the disappearance of IRW-A,it can be assumed that du/dxis constant,i.e.,

    In summary,the ion fluxes both to the cathode and anode in the third stage can be described by Eqs.(24)–(27).

    3.2.4.PIC modeling validations

    Based on the preceding discussions,there exist three successive stages in a typical ion extraction process, which are the electron oscillation and ion matrix sheath extraction stage(stage I),sheath expansion and IRW propagation stage(stage II),and plasma collapse stage(stage III).During stage II of the ion extraction process,considering the relationship betweenvsandcs,the extraction process near the cathode can be furtherly divided into stage II-1C and stage II-2C;while for the ion extraction process near the anode, it can be divided into stage II-1A and stage II-2A.In stage II-1A, the remaining ions in the matrix sheath will be extracted, while in stage II-2C, the ion extraction flux to the anode will remain unchanged and be equal to exp(-1).The formulas and the governing equations for the whole ion extraction process are summarized in Table 2.

    For validating the analytical model developed in this study(Table 2),the PIC simulations for the typical case listed in Table 1 are conducted.Comparisons of the ion extraction fluxes to the cathode and anode between the theoretically predicted results using the equations in Table 2 and the PIC modeling results are shown in Fig.9.It is seen from Fig.9 that: (i)The analytical model gives the same time-dependent trend of the ion extraction fluxes both to the cathode and the anode as that of the PIC modeling results fromt=0 totext;and in particular,two high flux peaks to the cathode and anode in stage I and a constant ion extraction flux to the anode in stage II-2A are both predicted by the analytical model and the PIC simulation.(ii)The percentages of the ions extracted to the cathode and anode in each stage are calculated theoretically and numerically, and the corresponding maximum relative discrepancies in each stage between the analytical model and the PIC simulation are all less than 3%.(iii) The ion extraction time obtained from the PIC simulation istext=335.13 when 90%of the ions are extracted;while att=text,the analytical model predicts that 87.8%of the ions are extracted with a relative discrepancy of 2.2%compared to that of the PIC modeling result.In summary,the preceding comparisons of the trend of the ion extraction flux,percentages of the ions extracted to the anode and cathode in each stage, and the percentage of the ions extracted att=textbetween the theoretical analysis and the PIC simulation have shown a good agreement.This indicates that the analytical model established in this paper can provide a good description on the whole ion extraction process.

    Table 2.List of the formulas and governing equations for describing the whole ion extraction process.

    4.Discussions on the characteristics of the ion extraction process

    As shown in Fig.9,the ion extraction fluxes to the cathode and anode both exhibit considerable peaks during stage I,which means that a large number of ions will be quickly extracted in the first several ion plasma periods; while in stage III,the magnitudes of the ion extraction fluxes are quite small.Therefore, in order to shorten the ion extraction time,it would be very effective to maximize the number of ions extracted in the first stage and reduce the number of ions extracted in the third stage.For the convenience of description,we define the ratio of the number of ions extracted in the first stage, as well as the total number of ions extracted beforet=tdis,anto the initial total number of ions asr1andr2, respectively.In the following two subsections, we will discuss the influence of plasma parameters on these two ratios.

    4.1.Proportion of ions extracted during stage I

    At the end of the first stage of the ion extraction process,there exist two Langmuir sheaths near the cathode and anode.The total number of ions in a Langmuir sheath is quite small compared with that in the initial matrix sheath.Thus,the ions left in the Langmuir sheath can be neglected when giving a rough estimation of the ions extracted during stage I.Therefore,r1can be approximately expressed as

    It is seen clearly thatr1increases with the increase of the magnitude of the externally applied voltage,and decrease with increasing in the initial plasma density and plasma width.Meanwhile,as for the initial electron temperature,it does not appear in Eq.(28),that is to say,it has almost no effects on the ion extraction characteristics in stage I for the ion extraction process.

    4.2.Proportion of remaining ions before the starting of stage III

    In this subsection,the number of the left ions at the beginning of the third stage will be roughly estimated.During stage II-2C, the moving speed of the sheath edge is lower thancs,and drops quickly with the increase of the sheath thickness.Neglecting the ions in the Langmuir sheath and assuming a uniform ion density distribution in the bulk plasma,r2can be roughly expressed as

    Letting(sdis,an-ssup,e,ca)/(tdis,an-tsup,e,ca)=α,the value ofαreflects the position of the sheath edge relative to the front of the IRW.The smallerαis, the more the front of the IRW exceeds the sheath edge.It is seen clearly thatαis affected by the externally applied voltage and the initial electron temperature.As the applied voltage increases, the sheath advances faster, and the value ofαalso increases.Moreover, the increase in the electron temperature will lead to a faster propagation of IRW,and a decrease inα.Thus,r2can be re-written as

    Based on the foregoing analysis,it is seen that the values ofsm,an/L,tm,e/L,tsup,e,ca/L,andssup,e,ca/Lall merely depend on the externally applied voltage.Therefore, ifαbeing constant,r2will be only affected by the applied voltage.Here,the values ofr2are calculated using Eq.(31) whenαranges from 0 to 1.0 with an interval of 0.1 andU0varies from 100 to 3000,while keeping other parameters being unchanged for the typical case listed in Table 1.The averaged value ofr2and its variation region enclosed byα=0.0 as the lower boundary andα=1.0 as the upper boundary are shown in Fig.10.For the typical case, the estimated averaged value ofr2obtained from Eq.(31)is 66.59%,which is in good agreement with the predicted value of 69.78%by the analytical model and that of 71.31%by the PIC modeling.

    Based on the preceding discussions,the analytical model listed in Table 2 can be employed to study the influences of the key physical parameters,including the externally applied voltage, initial electron temperature, initial plasma density, and plasma width, on the values ofr1andr2.It is shown from Fig.11 that: (i)The value ofr1is significantly positively correlated with ?U0,and negatively correlated with ?n0and ?L,while it remains almost unchanged with ?Te.(ii)The applied voltage has a significant impact onr2, while the variations of other parameters such as ?n, ?L, and ?Tehave minimal effects on it.It should be noted thatr2reflects the proportion of ions extracted before the plasma collapse, and the larger value ofr2is,the less ions need to be extracted in stage III,which is not directly related to the ion extraction time.Although the increasing in the initial electron temperature results in a higher ion acoustic velocity and shorter time,tdis,an,when the IRW-A and sheath meets,and consequently,leads to the reduction of the ion extraction time,text,as predicted in previous study,e.g.,Ref.[30], the ratios of the extracted ions are nearly independent of the initial electron temperature.The modeling results presented in Fig.11 suggest that the increase of ?U0will be most helpful for increasingr2,but in practice,the increase of ?U0is often limited by some factors such as the sputtering effect or secondary ionization caused by excessively high voltage.Therefore,novel methods need to be developed to improver2for the ion extraction process.

    5.Concluding remarks

    In this paper,a complete analytical model for describing the spatiotemporal evolutions of the ion extraction process in a decaying plasma bounded by a negatively biased cathode and a grounded anode is proposed,and also validated by comparing with the PIC simulation results.The evolutions of the ion density distribution, the sheath expansion and IRW propagation during the whole ion extraction process are summarized in Fig.12.The major conclusions are as follows:

    (i) The entire ion extraction process can be divided into three stages with the formation of the Langmuir sheath and the disappearance of the IRW-A after reaching the cathode sheath as the dividing points,i.e.,the electron oscillation and ion matrix sheath extraction stage,sheath expansion and IRW propagation stage,and plasma collapse stage.In stage II,depending on the relative values of the sheath expansion velocity and the ion acoustic velocity, the extraction process near the cathode can be further divided into two sub-stages,in which the sheath expansion is supersonic (stage II-1C) and/or subsonic (stage II-2C).While the ion extraction process near the anode can also be divided into two sub-stages according to whether the front of IRW reaches the anode or not.In the first sub-stage(stage II-1A),the remaining ions in the matrix sheath continue to be extracted,and the ion flux to the anode will be constant in the second sub-stage(stage II-2A).

    Fig.12.Schematic diagram of the spatiotemporal evolutions of the ion number density distribution, sheath expansion, and IRW propagation during the whole ion extraction process.

    (ii) Under the typical case studied in this paper, the ion extraction fluxes both to the cathode and anode obtained from the analytical model show a great agreement with the PIC simulation results.Considering the percentage of ions extracted att=text,the relative discrepancy between the analytical and numerical results is less than 2.2%.

    (iii)The ratio of ions extracted during the electron oscillation and matrix sheath extraction process increases with the increase of the externally applied voltage,and decreases with the increase of the initial plasma number density and plasma width;while the initial electron temperature shows little effect on the ion extractions features in stage I.The remaining ions before the plasma collapse are only significantly related to the externally applied voltage.The higher the applied voltage,the more ions are extracted before the plasma collapse,while the other parameters have little influences on this aspect.

    In practice, there may exist certain limitations when increasing the applied voltage to improve the ion extraction efficiency, and thus, the remaining ions before the plasma collapse cannot be efficiently reduced.Therefore, new methods based on the analytical model still need to be investigated in future research.In addition,since different electrode configurations,including Π type,M type,etc.,are often employed in actual applications,it is also indispensable to develop the twodimensional analytical model to describe the characteristics of the ion extraction process with complex electrode configurations facing various practical applications.

    Appendix A:Derivation of the ion density distribution in a Langmuir sheath

    As can be seen in Fig.A1, ignoring the electrons in a Langmuir sheath,the sheath function will be reduced to

    whereu0is the normalized ion velocity entering the Langmuir sheath.For a high voltage sheath,series expansion can be conducted forφ ?1,i.e.,

    By employing the same assumptions as the derivation for the Langmuir flux in Ref.[38]that the potential and the electric field are both equal to 0 at the sheath edge, the electric potential distribution in the Langmuir sheath is

    Assumingu0=1, the ion density distribution in a stable Langmuir sheath is then obtained as

    Fig.A1.Schematic diagram of the ion density, electron density, and electric potential distributions in a Langmuir sheath.

    Appendix B:Derivation of the ion density distribution in the interaction region of two IRWs

    In this section,we will consider such a system presented as Fig.B1, where region II refers to the interaction region,while regions I and III are the non-interacting regions only influenced by IRW-C and IRW-A,respectively.Att=0,IRW-C propagating in the positive direction will be excited atx=x1,while IRW-A propagating in the opposite direction will be excited atx=x2.

    Fig.B1.Schematic diagram of the ion density distribution after the meet of the two IRW fronts.

    Here we define the point with zero ion velocity asxc.Before the disappearance of IRW-A,xccan be expressed as

    According to Ref.[7], the ion number density in the interaction region exhibits a parabolic distribution.Considering the symmetry,the ion density in region II can be expressed as

    whereaandbare the parameters describing the parabolic relationship.

    After a traveling time oft, the front of IRW-A arrives atx=x4=x2-t, while the front of IRW-C arrives atx=x3=x1+t.The ion density and its derivate atx=x3can be obtained with the property of the IRW as

    Substituting Eq.(B3)into Eq.(B2),aandbcan be expressed as

    After the disappearance of the IRW-A (t=tdis,an),xccan be written as

    The ion number density and its derivate atx=x3are

    Substituting Eq.(B6)into Eq.(B5),aandbcan be re-written as

    The ratio of the difference between the ion density atx=xc[(n(xc)]andx=x3[n(x3)]ton(xc), denoted asη, can be employed to evaluate the uniformity of the ion density distribution in the interaction region,i.e.,

    Before the disappearance of the IRW-A,

    The rough range ofηis

    After the appearance of the IRW-C,

    The corresponding rough range ofηis

    Therefore,it is reasonable to approximate that the spatial distribution of the ion number density in the interaction region is uniform.

    Appendix C:Derivation of the ion extraction flux with a uniformly distributed ion number density and a linearly distributed ion velocity

    After the disappearance of the IRW-A, the ion number density and velocity exhibit a parabolic and piecewise linear distribution, respectively.From the analysis in Appendix B,it is known that the ion number density is almost uniformly distributed in the interaction region.Thus, in this section,the plasmas with uniformly distributed ion number density and linear distributed ion velocity will be considered.For the plasma near the anode,the electric field can be neglected and the ion thermal pressure can also be ignored due to the low temperature of ions.Therefore, the governing equations describing the ion motions can be written as

    Combining Eqs.(C1)and(C2),we can get that

    wherefis an arbitrary function ofx.Thus, the ion number density evolution will be

    wheren0refers to the ion number density att=0.So,the ion flux to the anode atx=Lwill be

    Considering the initial value of the ion number density and the ion velocity,equation(C5)can be re-written as

    Since the ion velocities atx=Landt=0 are both equal to 1,equation(C6)can be simplified as

    It can be seen from Eq.(C7) that the ion extraction flux on the anode is approximately proportional tot-2in this stage.

    Acknowledgement

    Project supported by the National Natural Science Foundation of China(Grant No.11775128).

    猜你喜歡
    東君和平
    和平之路
    小說的互文與改寫——讀東君短篇小說《與楊志共飲》
    都市(2022年12期)2022-03-04 09:11:56
    山茶花
    李東君
    中國書畫(2017年2期)2017-06-05 14:15:23
    溫暖的情境
    ——讀李東君的新作
    李東君作品
    博弈·和平
    特別文摘(2016年18期)2016-09-26 16:42:36
    期盼和平
    李東君作品選
    美術(shù)界(2014年7期)2014-04-29 00:44:03
    和平
    小說月刊(2014年2期)2014-04-18 14:06:40
    日韩有码中文字幕| 91久久精品国产一区二区成人 | 1024香蕉在线观看| 国产 一区 欧美 日韩| 琪琪午夜伦伦电影理论片6080| 亚洲欧洲精品一区二区精品久久久| 久久人人精品亚洲av| 两性午夜刺激爽爽歪歪视频在线观看| av黄色大香蕉| av女优亚洲男人天堂 | 国产激情久久老熟女| 精品免费久久久久久久清纯| 美女午夜性视频免费| 91老司机精品| 十八禁人妻一区二区| 精品久久久久久久末码| 午夜福利成人在线免费观看| netflix在线观看网站| 999精品在线视频| 国产伦精品一区二区三区视频9 | 欧美乱妇无乱码| 亚洲男人的天堂狠狠| 久久精品国产综合久久久| 三级国产精品欧美在线观看 | 日日摸夜夜添夜夜添小说| 黄色成人免费大全| 亚洲片人在线观看| 国产成人精品无人区| 国内精品久久久久久久电影| 免费一级毛片在线播放高清视频| 三级国产精品欧美在线观看 | 99国产精品99久久久久| 我的老师免费观看完整版| 色综合亚洲欧美另类图片| 亚洲九九香蕉| 亚洲欧美激情综合另类| 一级毛片高清免费大全| 好男人电影高清在线观看| 欧美黑人欧美精品刺激| 禁无遮挡网站| 国产高清三级在线| 日韩中文字幕欧美一区二区| 好男人电影高清在线观看| 日日夜夜操网爽| 亚洲精品美女久久久久99蜜臀| 中文字幕人成人乱码亚洲影| 一区二区三区高清视频在线| 免费看光身美女| 久久天躁狠狠躁夜夜2o2o| 一区二区三区国产精品乱码| 九色成人免费人妻av| 亚洲午夜精品一区,二区,三区| av女优亚洲男人天堂 | 国内久久婷婷六月综合欲色啪| 十八禁人妻一区二区| a级毛片a级免费在线| 国产高清激情床上av| 色哟哟哟哟哟哟| 国内精品一区二区在线观看| 欧美最黄视频在线播放免费| 人妻夜夜爽99麻豆av| 宅男免费午夜| 免费观看的影片在线观看| 色综合欧美亚洲国产小说| 欧美日本视频| 亚洲成人精品中文字幕电影| 日本a在线网址| 久久久久性生活片| 91在线观看av| 亚洲国产精品sss在线观看| 午夜福利免费观看在线| 桃色一区二区三区在线观看| 亚洲专区字幕在线| av国产免费在线观看| 国产精品女同一区二区软件 | 久久这里只有精品中国| 亚洲五月天丁香| 美女扒开内裤让男人捅视频| 国语自产精品视频在线第100页| 啦啦啦免费观看视频1| 午夜亚洲福利在线播放| 成人三级做爰电影| 亚洲午夜精品一区,二区,三区| 国产真人三级小视频在线观看| 91av网站免费观看| 久久久久国产一级毛片高清牌| 成人永久免费在线观看视频| 亚洲成人中文字幕在线播放| 免费搜索国产男女视频| 亚洲第一电影网av| 成人特级黄色片久久久久久久| 一进一出好大好爽视频| 国产精华一区二区三区| 中文字幕久久专区| 久久精品人妻少妇| 久久中文字幕一级| 精品日产1卡2卡| 亚洲va日本ⅴa欧美va伊人久久| 亚洲欧美激情综合另类| 国产精品野战在线观看| 欧美日本亚洲视频在线播放| 亚洲欧美日韩卡通动漫| 国产高清videossex| 欧美日韩综合久久久久久 | 国产v大片淫在线免费观看| 国产伦在线观看视频一区| 手机成人av网站| avwww免费| 欧美黑人巨大hd| 校园春色视频在线观看| 变态另类成人亚洲欧美熟女| 亚洲18禁久久av| 99在线视频只有这里精品首页| 亚洲人成网站高清观看| 午夜福利高清视频| 国产精品精品国产色婷婷| 热99re8久久精品国产| 在线国产一区二区在线| 欧美在线黄色| 日韩高清综合在线| 精品国内亚洲2022精品成人| 久久欧美精品欧美久久欧美| 男女之事视频高清在线观看| 丰满人妻熟妇乱又伦精品不卡| 黄色成人免费大全| 变态另类丝袜制服| 免费在线观看成人毛片| 少妇熟女aⅴ在线视频| 国产一区二区在线观看日韩 | 午夜精品久久久久久毛片777| 99精品在免费线老司机午夜| 人人妻人人澡欧美一区二区| 亚洲黑人精品在线| 黑人欧美特级aaaaaa片| 大型黄色视频在线免费观看| 少妇裸体淫交视频免费看高清| 夜夜夜夜夜久久久久| 国产免费男女视频| 一进一出好大好爽视频| 国产精品美女特级片免费视频播放器 | 午夜福利高清视频| 18禁美女被吸乳视频| 成人亚洲精品av一区二区| 老司机在亚洲福利影院| 国产成年人精品一区二区| 国产真实乱freesex| 欧美高清成人免费视频www| 脱女人内裤的视频| 99国产综合亚洲精品| 亚洲一区二区三区不卡视频| 国产综合懂色| 亚洲va日本ⅴa欧美va伊人久久| 特大巨黑吊av在线直播| 日韩大尺度精品在线看网址| 怎么达到女性高潮| 网址你懂的国产日韩在线| 免费av毛片视频| 后天国语完整版免费观看| 日韩欧美在线乱码| 国产精品久久久久久精品电影| 午夜精品在线福利| 国产精品99久久久久久久久| 日韩欧美三级三区| 亚洲人与动物交配视频| 国产黄片美女视频| 色综合站精品国产| 国内久久婷婷六月综合欲色啪| 我的老师免费观看完整版| 欧美性猛交╳xxx乱大交人| 特大巨黑吊av在线直播| 午夜视频精品福利| 亚洲国产欧美一区二区综合| 免费一级毛片在线播放高清视频| 国产精品av视频在线免费观看| 青草久久国产| 老司机深夜福利视频在线观看| 1024香蕉在线观看| av片东京热男人的天堂| 国产亚洲av高清不卡| 国产亚洲欧美在线一区二区| 欧洲精品卡2卡3卡4卡5卡区| 此物有八面人人有两片| 亚洲五月天丁香| 日日摸夜夜添夜夜添小说| 巨乳人妻的诱惑在线观看| 桃红色精品国产亚洲av| 又黄又粗又硬又大视频| 天堂√8在线中文| 成年女人毛片免费观看观看9| 女同久久另类99精品国产91| 99国产综合亚洲精品| 伦理电影免费视频| 亚洲国产精品999在线| 欧美xxxx黑人xx丫x性爽| 国产久久久一区二区三区| 亚洲 欧美一区二区三区| xxxwww97欧美| 真人做人爱边吃奶动态| 动漫黄色视频在线观看| 久久精品aⅴ一区二区三区四区| 亚洲狠狠婷婷综合久久图片| 色精品久久人妻99蜜桃| 久久草成人影院| 国产亚洲精品av在线| 身体一侧抽搐| 国产精品免费一区二区三区在线| а√天堂www在线а√下载| 桃色一区二区三区在线观看| 中文字幕av在线有码专区| 国产极品精品免费视频能看的| 热99在线观看视频| 色播亚洲综合网| 精品久久蜜臀av无| 国产爱豆传媒在线观看| 在线观看午夜福利视频| 超碰成人久久| 国产成人福利小说| 天堂影院成人在线观看| 午夜福利欧美成人| 久99久视频精品免费| 国产精品99久久99久久久不卡| 欧美日韩中文字幕国产精品一区二区三区| 亚洲男人的天堂狠狠| www.自偷自拍.com| 一本一本综合久久| 国产精品一区二区精品视频观看| 香蕉丝袜av| 亚洲人成网站高清观看| 久久精品国产亚洲av香蕉五月| 中文字幕精品亚洲无线码一区| 国产精品98久久久久久宅男小说| 在线观看美女被高潮喷水网站 | 免费大片18禁| 国产欧美日韩精品一区二区| 免费电影在线观看免费观看| 欧美另类亚洲清纯唯美| 人人妻人人澡欧美一区二区| 国产精品美女特级片免费视频播放器 | 日韩三级视频一区二区三区| 国产 一区 欧美 日韩| 精品久久久久久,| 亚洲成av人片免费观看| 婷婷精品国产亚洲av| 精品久久久久久久毛片微露脸| 精品久久久久久久久久久久久| 亚洲精品在线观看二区| 丝袜人妻中文字幕| 久久久精品大字幕| 婷婷亚洲欧美| 久久精品亚洲精品国产色婷小说| 男人舔女人下体高潮全视频| 午夜精品在线福利| 好看av亚洲va欧美ⅴa在| 一区福利在线观看| 给我免费播放毛片高清在线观看| 日日夜夜操网爽| av天堂在线播放| 一个人观看的视频www高清免费观看 | 国产成人啪精品午夜网站| 小说图片视频综合网站| 国产精品日韩av在线免费观看| 在线看三级毛片| 亚洲一区二区三区不卡视频| 成人国产一区最新在线观看| 听说在线观看完整版免费高清| 亚洲成人久久爱视频| 亚洲av片天天在线观看| 国产三级在线视频| 亚洲aⅴ乱码一区二区在线播放| 久久草成人影院| 亚洲色图av天堂| 中文字幕人成人乱码亚洲影| 麻豆av在线久日| 淫妇啪啪啪对白视频| 一进一出抽搐gif免费好疼| 久久精品综合一区二区三区| 99热6这里只有精品| 日韩精品中文字幕看吧| 性色av乱码一区二区三区2| 国产激情偷乱视频一区二区| www.999成人在线观看| 欧美一级毛片孕妇| 国产男靠女视频免费网站| av女优亚洲男人天堂 | 亚洲成人久久爱视频| 国产精品久久久久久久电影 | 午夜成年电影在线免费观看| 美女免费视频网站| 国产成人欧美在线观看| 亚洲精品一区av在线观看| 久久久国产成人免费| 99视频精品全部免费 在线 | 一二三四社区在线视频社区8| 成人18禁在线播放| 热99在线观看视频| 欧美中文综合在线视频| 国产高潮美女av| 久久中文字幕一级| 久久精品91无色码中文字幕| 亚洲五月婷婷丁香| 久久久精品大字幕| 国产av在哪里看| 亚洲aⅴ乱码一区二区在线播放| 日日摸夜夜添夜夜添小说| 午夜a级毛片| 国产精品美女特级片免费视频播放器 | 国产精品乱码一区二三区的特点| 少妇丰满av| 这个男人来自地球电影免费观看| 黑人欧美特级aaaaaa片| 人妻夜夜爽99麻豆av| 亚洲九九香蕉| 国内久久婷婷六月综合欲色啪| 成人特级黄色片久久久久久久| 久久精品影院6| 午夜激情福利司机影院| 网址你懂的国产日韩在线| 日本精品一区二区三区蜜桃| 久久精品夜色国产| 黄色配什么色好看| 九九久久精品国产亚洲av麻豆| 亚洲高清免费不卡视频| 天堂网av新在线| 国产精品久久久久久精品电影| 欧美另类亚洲清纯唯美| 欧美bdsm另类| 欧美另类亚洲清纯唯美| 精华霜和精华液先用哪个| 欧美3d第一页| 国产av在哪里看| 精品无人区乱码1区二区| 成人毛片a级毛片在线播放| 午夜福利视频1000在线观看| 国产探花极品一区二区| 欧美精品国产亚洲| 久久人人爽人人爽人人片va| 观看美女的网站| 久久人人爽人人片av| 少妇高潮的动态图| 国产av不卡久久| 日韩一区二区三区影片| 大香蕉久久网| 97在线视频观看| 别揉我奶头 嗯啊视频| 91午夜精品亚洲一区二区三区| 最新中文字幕久久久久| 热99re8久久精品国产| 精品熟女少妇av免费看| 一个人观看的视频www高清免费观看| 18+在线观看网站| 国产午夜精品久久久久久一区二区三区| 美女cb高潮喷水在线观看| 日韩视频在线欧美| 成人三级黄色视频| 天天躁日日操中文字幕| 直男gayav资源| 人人妻人人澡人人爽人人夜夜 | 亚洲丝袜综合中文字幕| 深夜a级毛片| 六月丁香七月| 国产又黄又爽又无遮挡在线| 男人舔奶头视频| 国产免费视频播放在线视频 | 国产成人一区二区在线| 99久久精品热视频| 亚洲国产精品成人久久小说| 熟妇人妻久久中文字幕3abv| 色综合亚洲欧美另类图片| 最近的中文字幕免费完整| 一个人看的www免费观看视频| av又黄又爽大尺度在线免费看 | 黄色日韩在线| 久久精品综合一区二区三区| 高清毛片免费看| 最新中文字幕久久久久| 亚洲图色成人| 亚洲一级一片aⅴ在线观看| 久久久a久久爽久久v久久| 久久久久久久久大av| 国产大屁股一区二区在线视频| 亚洲综合色惰| 国产精品麻豆人妻色哟哟久久 | 熟女人妻精品中文字幕| 国产三级中文精品| 啦啦啦啦在线视频资源| 久久精品国产自在天天线| 欧美成人免费av一区二区三区| 欧美bdsm另类| 中文在线观看免费www的网站| 高清日韩中文字幕在线| 亚洲三级黄色毛片| 男女国产视频网站| 国产伦在线观看视频一区| 搡老妇女老女人老熟妇| 国产精品1区2区在线观看.| 一区二区三区免费毛片| 国产极品精品免费视频能看的| 国产精品熟女久久久久浪| 日本黄大片高清| 99在线视频只有这里精品首页| ponron亚洲| 少妇熟女欧美另类| 日本午夜av视频| 国产精品久久电影中文字幕| 亚洲av男天堂| 黄色欧美视频在线观看| 亚洲最大成人av| 熟妇人妻久久中文字幕3abv| 嫩草影院入口| 亚洲第一区二区三区不卡| 欧美色视频一区免费| 日韩av不卡免费在线播放| a级毛色黄片| 哪个播放器可以免费观看大片| 亚洲人成网站在线播| 亚洲av成人精品一二三区| 中文字幕人妻熟人妻熟丝袜美| 免费不卡的大黄色大毛片视频在线观看 | 秋霞伦理黄片| 观看美女的网站| 亚洲色图av天堂| 一本久久精品| 22中文网久久字幕| 国产69精品久久久久777片| 亚洲成人中文字幕在线播放| 亚洲国产精品国产精品| 亚洲av福利一区| 国产在线男女| 国产黄片视频在线免费观看| 人体艺术视频欧美日本| 欧美色视频一区免费| 午夜日本视频在线| 成人亚洲精品av一区二区| 日本wwww免费看| 最近最新中文字幕大全电影3| 一边亲一边摸免费视频| 国产成人一区二区在线| 精品久久久久久久末码| 欧美精品一区二区大全| 婷婷色av中文字幕| 日本wwww免费看| 亚洲在线自拍视频| 久久久精品94久久精品| 亚洲国产精品合色在线| 亚洲va在线va天堂va国产| 日韩欧美 国产精品| 国产视频内射| 免费播放大片免费观看视频在线观看 | 三级男女做爰猛烈吃奶摸视频| 午夜福利视频1000在线观看| 色吧在线观看| 高清在线视频一区二区三区 | av线在线观看网站| 又粗又硬又长又爽又黄的视频| 午夜视频国产福利| 欧美激情在线99| 最近中文字幕2019免费版| 好男人在线观看高清免费视频| 欧美又色又爽又黄视频| 欧美高清性xxxxhd video| 日韩三级伦理在线观看| 欧美日本视频| 久久午夜福利片| 国产高潮美女av| 免费黄色在线免费观看| www日本黄色视频网| 乱码一卡2卡4卡精品| 国产黄a三级三级三级人| 中国美白少妇内射xxxbb| 久久韩国三级中文字幕| 午夜福利网站1000一区二区三区| 麻豆乱淫一区二区| 日韩制服骚丝袜av| 看片在线看免费视频| 最近中文字幕高清免费大全6| av福利片在线观看| 熟女人妻精品中文字幕| 少妇被粗大猛烈的视频| 夜夜爽夜夜爽视频| 天美传媒精品一区二区| 国产精品久久久久久久电影| 免费黄色在线免费观看| 热99在线观看视频| 国产精品野战在线观看| 欧美zozozo另类| 永久网站在线| 成年免费大片在线观看| 国产精品久久久久久精品电影| 美女被艹到高潮喷水动态| 视频中文字幕在线观看| 久久精品91蜜桃| 99热6这里只有精品| 精品国产三级普通话版| 日韩一区二区视频免费看| 午夜亚洲福利在线播放| 免费观看在线日韩| 成年免费大片在线观看| 最近最新中文字幕免费大全7| 人人妻人人看人人澡| 亚洲欧美日韩卡通动漫| 国产免费一级a男人的天堂| 麻豆av噜噜一区二区三区| 亚洲伊人久久精品综合 | 99久久中文字幕三级久久日本| 赤兔流量卡办理| 91精品伊人久久大香线蕉| 身体一侧抽搐| 如何舔出高潮| 最近最新中文字幕免费大全7| kizo精华| 丰满少妇做爰视频| 汤姆久久久久久久影院中文字幕 | 一级毛片aaaaaa免费看小| 一夜夜www| 永久网站在线| 极品教师在线视频| av又黄又爽大尺度在线免费看 | 欧美日韩综合久久久久久| 三级国产精品欧美在线观看| 亚洲五月天丁香| 在线a可以看的网站| 一个人观看的视频www高清免费观看| av天堂中文字幕网| 欧美成人一区二区免费高清观看| 欧美不卡视频在线免费观看| 亚洲av熟女| 看免费成人av毛片| 一区二区三区四区激情视频| 伦理电影大哥的女人| 欧美成人a在线观看| 在线观看66精品国产| 亚洲va在线va天堂va国产| 色网站视频免费| 午夜老司机福利剧场| 你懂的网址亚洲精品在线观看 | 一级av片app| 亚洲精品成人久久久久久| 久99久视频精品免费| 免费av不卡在线播放| 一本久久精品| 1000部很黄的大片| 51国产日韩欧美| 日日啪夜夜撸| 啦啦啦啦在线视频资源| 精华霜和精华液先用哪个| 中文字幕制服av| 亚洲av福利一区| 久久鲁丝午夜福利片| 国内精品宾馆在线| 少妇被粗大猛烈的视频| 麻豆一二三区av精品| 国产成人a∨麻豆精品| 亚洲精品日韩在线中文字幕| 国产av一区在线观看免费| 激情 狠狠 欧美| 你懂的网址亚洲精品在线观看 | 亚洲精品乱久久久久久| 日本熟妇午夜| 国产精品美女特级片免费视频播放器| 日韩制服骚丝袜av| 亚洲性久久影院| 老司机影院成人| 久久精品夜夜夜夜夜久久蜜豆| 亚洲av中文字字幕乱码综合| 亚洲国产最新在线播放| 村上凉子中文字幕在线| 中文字幕久久专区| 欧美成人a在线观看| 国产精品国产三级专区第一集| av在线老鸭窝| 日本爱情动作片www.在线观看| 日韩欧美 国产精品| av线在线观看网站| 超碰av人人做人人爽久久| 我的老师免费观看完整版| 欧美高清性xxxxhd video| 舔av片在线| 国产一级毛片七仙女欲春2| 久久久久久久国产电影| 亚洲最大成人av| 天堂中文最新版在线下载 | 国模一区二区三区四区视频| 村上凉子中文字幕在线| 国产淫语在线视频| 欧美成人一区二区免费高清观看| 91aial.com中文字幕在线观看| 国产成人a区在线观看| 日韩av在线大香蕉| 亚洲真实伦在线观看| 一级爰片在线观看| 亚洲欧美日韩卡通动漫| 日日啪夜夜撸| 天堂网av新在线| 六月丁香七月| 91aial.com中文字幕在线观看| 亚洲丝袜综合中文字幕| 日韩高清综合在线| 日韩中字成人| 18+在线观看网站| 国产精品嫩草影院av在线观看| 国产不卡一卡二| 国产精品不卡视频一区二区| 亚洲欧美日韩高清专用| 亚洲欧美成人综合另类久久久 | 午夜免费激情av| 五月伊人婷婷丁香| 亚洲精华国产精华液的使用体验| 午夜精品在线福利| 亚洲怡红院男人天堂| 内射极品少妇av片p| 最近的中文字幕免费完整| 国产一区亚洲一区在线观看| 国语自产精品视频在线第100页| 亚洲一区高清亚洲精品| 成年女人永久免费观看视频|