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

    Ion-focused propagation of a relativistic electron beam in the self-generated plasma in atmosphere

    2022-06-29 08:54:16JianHongHao郝建紅BiXiXue薛碧曦QiangZhao趙強(qiáng)FangZhang張芳JieQingFan范杰清andZhiWeiDong董志偉
    Chinese Physics B 2022年6期
    關(guān)鍵詞:張芳趙強(qiáng)

    Jian-Hong Hao(郝建紅) Bi-Xi Xue(薛碧曦) Qiang Zhao(趙強(qiáng))Fang Zhang(張芳) Jie-Qing Fan(范杰清) and Zhi-Wei Dong(董志偉)

    1School of Electrical and Electronic Engineering Department,North China Electric Power University,Beijing 102206,China

    2Institute of Applied Physics and Computational Mathematics,Beijing 100094 China

    Keywords: ion-focused regime (IFR), charge neutralization, electron avalanche, relativistic electron beam(REB),particle-in-cell Monte Carlo collision(PIC-MCC)

    1. Introduction

    During long-range propagation of an artificial relativistic electron beam(REB)in space environment,multiple practical applications in active space experiments can be achieved by probing secondary particles of interactions between primary beam electrons with ambient neutrals,including tracing magnetic field lines and analyzing beam–particle interactions.[1–5]Due to strong space charge effects of intense REBs, it is difficult to achieve long-range propagation in the order of several hundreds of kilometers in vacuum. In addition,owing to degradation of beam quality,REBs are unable to meet need of practical applications after long-range propagation.Therefore,how to increase beam transport efficiency in low-pressure atmosphere during long-range propagation has always been the focus of associated research.[6]As the ambient neutral density of atmosphere varies in a complicated manner at different altitudes,propagation properties of REBs can be affected by multiple factors.Previous studies have discussed long-range propagations of REBs in atmosphere from different perspectives thoroughly.[7–12]We also gave the Monte–Carlo(MC)simulation,particle-in-cell(PIC)simulation and theoretical study on the influence of ambient neutral density, intensity and direction of geomagnetic field, and the ion channel oscillation on the long-range propagation of REBs.[13–15]

    After emission of an REB into low-pressure atmosphere environment, the impact ionization between primary beam electrons and ambient particles produces positive ions and secondary electrons in the beam.Under the influence of Coulomb repulsion, secondary electrons in the beam will be expelled transversely and the heavier positive ions will be left on the beam path. Accordingly, an ion-focused regime (IFR) is established. In the low-pressure atmosphere, ambient particle density varies at a relatively slow rate, meaning that the IFR could maintain for a longer time in such environment. The IFR also has multiple applications in devices such as plasmaassisted microwave sources.[16–24]However,the REB and ambient neutral gas require to fulfill some conditions for generation of an IFR.Previous studies indicated that in neutral gas with a relatively low pressure, it is impossible to accumulate enough positive ions for the self-focused propagation of an REB in a short period of time. In the neutral gas with a relatively high pressure,the amounts of both secondary electrons and ions increase rapidly,leading to accumulation of multiple secondary electrons in the beam. Such a phenomenon may lead to various instabilities,including the two-stream instability and resistive hose instability,which will result in a significant decrease in the beam propagation efficiency.[25,26]Therefore, it is necessary to provide a more detailed description of the connection between physical characteristics of the IFR and parameters of REBs.Studies mentioned above mainly concentrated on device physical problems, such as the propagation of REBs in an enclosed environment. Under such conditions,transverse velocity of secondary electrons will be lost at the drift tube wall. However, for the REB which propagates in an open environment such as the Earth atmosphere,secondary electrons are able to move transversely to a further position and a reverse return current will be generated on the surface of the beam. In PIC simulations, this difference should be taken into account by setting appropriate boundary conditions and dimensions of the simulation region. For the longitudinal motion of secondary electrons in the beam and other complex physical processes,associated analytical studies typically neglect these effects for ease of calculation. A PIC simulation code can take into account the influence of more physical fields, thus a particle-in-cell Monte Carlo collision(PICMCC)method is used in this study to calculate the propagation of REBs in low-pressure neutral gas and to analyze characteristics of the IFR generated by the beam.

    In the present study, our aim is to explore the possibility of long-range propagation of REBs by utilizing the IFR in the atmosphere. First, the basic physical processes during long-range propagation of REBs and envelope equations of the sheet beam are briefly introduced in Section 2. Long-range propagation of REBs in low-pressure neutral gas is numerically investigated by using a PIC-MCC code. The establishment of IFR and the corresponding influence on the long-range propagation of REBs are analyzed by combination with associated theories. Finally, possible solutions for the realization of long-range propagation of REBs in the atmosphere are discussed.

    2. The calculation model

    2.1. Collision processes and physical model

    2.1.1. Electron impact ionization and electron avalanche

    In accordance with previous studies, parameters of the REB and neutral gas used in this study are listed in Table 1.Neutral gas is composed of 70%of N2and 30%of O2to simulate the Earth atmosphere. During the propagation of the REB in the neutral gas, gas atoms become ionized because of collisions with beam electrons and ions,and secondary electrons are produced. In the impact ionization,the ion density can be modeled as[26]

    whereniandnbare the densities of ions and beam electrons,respectively;tbis the propagation time of the beam after it is injected into the gas, andτeis the characteristic electron impact ionization time expressed as follows:

    Note thatngis the density of neutral gas;σbis ionization crosssection;vis the velocity of electrons. Ionization time of the REB with parameters in Table 1 is approximately 1ns, which is much less than the longitudinal and transverse ion oscillation time.Thus,ions can be considered immobile in this study.

    Table 1. Physical parameters of the beam and plasma.

    2.1.2. Electron avalanche

    After the ionization between beam electrons and neutral gas atoms, secondary electrons will be accelerated and expelled transversely out of the beam by the REB.Motion of secondary electrons will lead to further ionization of background gas,which can be expressed by

    wherensecis the density of secondary electrons;teandtsare the effective electron avalanche time and the effective secondary escape time, respectively. Before the beam is totally charge neutralized,tsbecomes relatively small, which means that secondary electrons will quickly escape from the beam transversely. In this case, electron avalanching effects can be usually neglected.

    2.2. The PIC-MCC model

    Compared to the previous PIC research on the beam transportation in an enclosed environment, the present study is mainly concentrated on the influence of IFR on the REB dynamics in the atmosphere environment. Without any shielding boundaries in low-pressure neutral gas,it has been proved that the significant presence of numerical noise may result in decrease of the calculation accuracy in a fully electromagnetic model. Considering the cost in computational capacity when calculating the long-range propagation, a PIC-MCC quasielectromagnetic model in the so-called 2D3V (2D in spatial and 3D in velocity)formulation is used to simulate the propagation of the REB in neutral gas. A snapshot of PIC-MCC simulations is shown in Fig.1,whereLxandLyare the transverse and longitudinal dimensions. In addition, the moving window technology is employed in PIC simulations to meet the need in propagation distance of the beam,whereas a cylindrical coordinate system cannot be collocated with the moving window. Therefore, a two-dimensional Cartesian coordinate system is set to depict the particle distribution, in which thexandyaxes are along the propagation direction and transverse direction of the REB, respectively. Thezaxis is auxiliary direction. A sheet beam model is used accordingly.The beam electrons are uniformly distributed in space, and the beam is infinitely long in the longitude direction. Fromtb=0,the beam is continuously emitted from the left boundary(x=0)along the positive direction ofx-axis, which indicates thattbcorresponds to the beam propagation time. The beam parameters are listed in Table 1. The total size and grids of the simulation satisfy the Courant–Friedrichs–Lewy(CFL)condition.[27]At the left boundary, we use a Dirichlet (i.e.,conducting) boundary condition on the electrostatic potentialψx=0=0. At the right boundary,we use a Neumann boundary conditionψ′x=Lx=0. The top and bottom boundaries are set to be periodic.

    Fig. 1. Snapshot of simulation particle positions: (a) distribution of beam electrons and secondary ions, (b) distribution of beam electrons and secondary electrons.

    To investigate the collisions between particles, a Monte Carlo collision technique is used, including the following reactions: (1)beam electron,secondary electron impact ionizations

    (2)beam electron,secondary electron elastic collisions

    (3)beam electron,secondary electron attachments

    Ionization cross sections of low-energy electrons (E <10 keV) are obtained from experiments, and the cross sections of high-energy electrons(E >10 keV)are calculated by Coulomb collision cross section.[28,29]In the PIC-MCC simulations, “macroparticles” represent several particles whose physical characteristics are close to each other. During the long-range propagation, electron avalanche can lead to a dramatic increase of secondary electrons, which further results in a massive reduction of computational speed. Null reactions“combination”can decrease the amount of the macroparticles in the area with a high particle density by increasing the weights of the macroparticles. In this way, the computation efficiency and ensuring precision can be improved to satisfy the need of our calculation.

    2.3. The analytical model

    The envelope equations in the paraxial axial approximation can describe the envelope variations of REBs during longrange propagation. If the energy loss caused by the beam–atmosphere collisions is ignored,envelope equations of a sheet beam can be written as[30]

    whereIbandIpare the beam current and return current,respectively.Previous studies usually assume that the beam is charge neutral and the current neutralization is negligible. In order to analyze the influence of beam neutralization on the long-range propagation,neutralization factors,which vary with propagation time and longitudinal positions, are obtained in accordance with the data in PIC simulations. By taking neutralization factors into envelope equations,the results can be verified by comparing the PIC results with the results obtained by envelope equations.

    3. Charge and current neutralizations in IFR

    3.1. Charge neutralization

    Establishment processes of IFR can be described by the transverse velocity variations of secondary electrons at different propagation timetb. Figure 2 shows the transverse velocity variations of secondary electrons whentb=9 ns andtb=57 ns. In Fig. 2(a), as the beam propagates for a short time,the beam has not been totally charge neutralization.Most secondary electrons are expelled transversely by the beam. At this time, most secondary electrons located at the positiveyaxis transport along the positive direction or conversely.Whentb=57 ns,beam body of the REB has been totally charge neutralization,leading to the accumulation of excess positive ions in the beam. These ions can effectively restrict the radical motion of newly generated secondary electrons, or even pull the other secondary electrons back to the beam to neutralize the beam. Consequently, transverse velocities of secondary electrons are more homogenously distributed along theyaxis with a much smaller peak amplitude,as depicted in Fig.2(b).

    Fig. 2. Transverse velocities of secondary electrons at different transverse positions: (a)tb=9 ns,(b)tb=57 ns.

    Distribution of the REB is shown in Fig.3 and densities of various particles in the beam are shown in Fig.4. At the beam head, as the duration of impact ionization is relatively short,density of secondary particles is lower than beam electrons,which indicates that the IFR has not been well established.The beam expands transversely under the influence of space charge effects and initial emittance, resulting in a density of beam electrons lower than the initial density. As the coordinatexdecreases, the duration of impact ionization improves gradually, leading to the increase of secondary particle densities. Furthermore, in accordance with Eq. (6),feincreases with the ion densityni.At the beam neck,whenfesatisfies the so-called Budker condition,i.e.,fe>1/γ2,the corresponding phenomenon called Bennett pinch occurs because of the transverse pinch force supposed by the IFR.The transverse dimension can be reduced to initial size gradually. At the junction part between the beam neck and beam body,the beam can be over pinched(fe>1). The significantly over pinched part of the beam may even break up, which eventually interrupts the efficient beam propagation. In the region of beam body,densities of various secondary particles exceed the beam electron density and the transverse dimension of the beam is close to the initial size,which indicates that the IFR has been fully established.

    Fig.3. Distribution of the REB when tb=57 ns.

    Fig.4. Relationship between the densities of various particles and longitudinal positions.

    3.2. Current neutralization

    During the long-range propagation of the REB in the plasma environment, a longitudinal electric fieldEzis generated by net currentInetnear the beam. The return currentIp,which flows in opposite directions against the beam,is driven byEz;le=c/ωpis the so-called electromagnetic skin depth,which is also the thickness of the sheath near the edge of the beam to which the net currentIpis confined;[31,32]ωpis the plasma frequency. The relationship between the skin depthleand the transverse dimension of the REB can be modeled as[33]

    Variations of the above currents in PIC-MCC simulations are shown in Fig.5,whereInet=Ib+Ip. At the front of the beam head,beam electrons are accelerated because of the longitudinal space charge effects of other beam electrons,and the first peak beam current is formed. In the beam neck area,with the Bennett pinch of the REB, beam density increases gradually with the decrease of the coordinatexof the beam electrons,and finally the second peak of beam current is formed. In the connection area between the beam neck and beam body, the beam current gradually approaches the initial beam current.During the above processes of the change in beam current,the return current is generated by the motion of secondary electrons to suppress the effect of beam current variation on generated plasma environment. The return currentIpcan be written as[32]

    Fig.5. Variations of Ib,-Ip,and Inet in the beam.

    Previous studies use the electric diffusion time of the plasma to describe the growth period of return current,which can be written as[31]

    Here,reis the classical electron radius,andveis the effective collision frequency which can be estimated by the formula[34]

    wherengis the density of neutral gas. Figure 6 shows the variation ofτewith the propagation time. With the increase oftb, transverse and longitudinal dimension of the beam increase because of the space charge effects. Thus, the density of beam electrons decreases. Reduction ofneleads to the increase ofτe. This regularity can be shown both in PIC-MCC simulations and the analytical model.

    Fig.6. Variation of τe with the propagation time tb.

    The formulas of the two neutralization factors are given as Eqs. (6) and (7). Combined the densities of various particles and distributions of various currents in Figs. 4 and 5,variations of neutralization factors with longitudinal positions are shown in Fig. 7. During the long-range propagation of a 1 MeV REB in neutral gas of 10 Pa,impact ionizations between the beam and ambient particles make the charge neutralization factor beyond the Budker condition in a short duration(<1 ns)and reach its peak in the connection area between the beam neck and beam body. Ions and secondary electrons in the beam change the value offevarying around 1 in the beam body, which means that the subsequent part of the beam is totally charge neutralized. Under the influence of the longitudinal motion of both beam electrons and secondary electrons,fmreaches its maximum at the tail of the beam.

    Fig.7. Variations of neutralization factors fe and fm with the longitudinal positions.

    4. Propagation properties of the REB in IFR

    4.1. Influence of IFR on the longitudinal motion of the REB

    After the REB is emitted into the vacuum or the neutral gas, the distribution of longitudinal velocities of beam electrons shows certain features. This phenomenon is more evident in long-range propagation of REBs,but only a few studies have been focused on such influence of the IFR. Figure 8 shows distributions of longitudinal velocities of beam electrons whentb=15 ns andtb=57 ns.Whentb=15 ns,the duration of impact ionizations between beam electrons and ambient particles is not long enough for the establishment of IFR.Under the influence of longitudinal space-charge effects,beam electrons which are just emitted into the gas will be decelerated by other beam electrons. On the contrary,beam electrons at the front of the beam head are accelerated by other beam electrons. Without the effect of any other external forces,the longitudinal space-charge effects can cause the increase of longitudinal dimension, resulting in a low transport efficiency.Whentb= 57 ns, at the beam head and neck, because the density of ions is not high enough to affect the propagation of REB, distribution of longitudinal velocities in this area is similar totb=15 ns. In the connection area between the beam neck and body,a peak of velocity appears in the position wherefe>1 in Fig.7,which indicates that the excessive ions in the beam attract the beam electrons,leading to the transverse over pinch and longitudinal acceleration. For beam electrons in the beam body, the above discussion in Subsection 3.1 indicates that the beam is neutralized under the influence of IFR,which means that the longitudinal space-charge effects are offset accordingly. Consequently, the longitudinal velocities are close to the initial velocity of the beam.It can be concluded that IFR can neutralize the longitudinal space-charge effects and limit the longitudinal expansion of the beam efficiently. Moreover,combining the previous calculation results,[15]when the REB propagates in the natural plasma environment in the atmosphere,ion channel always continues to oscillate periodically over the entire propagation,leading to the synchronous oscillation of longitudinal velocities of beam electrons. When the beam propagates in the self-generated plasma in neutral gas,as the density of secondary particles varies significantly at different longitudinal positions,the IFR will not oscillate at a certain frequency. Hence,the longitudinal velocity keeps close to the initial velocity,as depicted in Fig.8(b). In these two conditions, although both the ion channel and the IFR can focus the propagation of the REB similarly, amplitude of the longitudinal velocity gradually increases during the long-range propagation in the natural plasma, which eventually leads to significant instabilities as well as a low transport efficiency.This phenomenon does not appear when the REB propagates in the self-generated plasma.

    Fig. 8. Distributions of longitudinal velocities of beam electrons: (a)tb=15 ns,(b)tb=57 ns.

    4.2. Influence of IFR on the transverse motion of the REB

    Comparison of the beam envelopes between PIC-MCC simulations and envelope equations is shown in Fig.9. At the beam head,complicated variations of various particle densities result in the certain fluctuations of charge neutralization factorfe, some errors are expected in the calculated result of envelope equations,while the changing rules of these two methods are consistent. There is an obvious transverse expansion at the beam head. Whenfesatisfies the Budker condition,the transverse expansion of the beam is limited by IFR and the beam envelope gradually approaches the initial value. In the beam body area,beam envelope tends to be stable under the common effects of multiple forces. The equilibrium transverse dimension is close to the initial transverse dimension.

    Fig.9. Comparison of beam envelopes calculated by PIC-MCC simulations and envelope equations.

    5. Problems faced by the long-range propagation of the REB in the earth atmosphere

    5.1. Instabilities of the REB in IFR

    According to Section 4, IFR can contribute to the longrange propagation of REB by limiting the transverse and longitudinal expansion. However, after the REB has propagated in neutral gas for a period,electron avalanche appears because of the continuous occurrence of impact ionizations and may disrupt the stable propagation of the REB. Figure 10 shows the distribution of beam electrons whentb= 75 ns. Compared with Fig.3,the obvious transverse and longitudinal expansion can be observed, and a significant instability occurs at the tail of the beam. With the increase oftb, this instability will spread through the entire beam and totally disrupt the propagation. Figure 11 shows the variation of secondary electron densities in different longitudinal positions withtb. The above discussion in Sections 3 and 4 indicates that electron avalanche begins when the beam has finished charge neutralization in the beam body area,resulting in the exponential increase and accumulation of secondary electrons in the beam.Multiple secondary electrons may lead to a two-stream instability with beam electrons and further other instabilities. The increase of secondary electrons also improves the return current as well as the current neutralization factorfm. Current neutralization may weaken the focused effect of IFR and lead to the occurrence of resistive hose instability. Multiple factors mentioned above eventually result in the macro instability in Fig.10.

    Fig.10. Distribution of the REB when tb=75 ns.

    Combined the results of both theoretical analyses and PIC-MCC simulations, in order to avoid the occurrence of such instabilities, the decrease of pulse width Δtof the REB can be taken into account to reduce the duration of impact ionization at the beam tail.Distributions of REBs from Δt=72 ns to Δt=65 ns whentb=75 ns are shown in Fig. 12. With the decrease of Δt,the influence of electron avalanche recedes gradually. When Δt <65 ns, the effect of electron avalanche at the beam tail can be basically ignored. It can be concluded that reduction of the pulse width can be used to avoid the disruption of beam stability by electron avalanche.

    Fig. 11. Variation of secondary electron densities with tb at different longitudinal positions.

    Fig. 12. Distributions of REBs when tb =75 ns: (a) Δt =72 ns, (b)Δt=70 ns,(c)Δt=68 ns,and(d)Δt=65 ns.

    5.2. Reduction of beam quality during long-rang propagation

    As the REB propagates to a longer distance in the neutral gas, beam-atmosphere interactions will significantly reduce the beam quality. Figure 13 presents the beam electron density at the beam head,beam fluence and charge neutralization factorfein the beam body area. Under the influence of impact ionizations,elastic collisions,and space charge effects of the beam, the beam head will expand consistently and the beam electron density decreases correspondingly. As a result,number of ions generated by the beam head is not enough for the establishment of the stable IFR in the beam body. The beam quality reduces gradually, as shown by the decrease of fluence andfe.

    A possible solution is to use a train of pulses for longrange propagation in the atmosphere to handle the problem of beam quality and the reduction during long-range propagation.[32]The early emitted pulses are used to generate an IFR in the atmosphere with a similar mechanism to the beam head in the above experience. By sacrificing several early emitted pulses, a stable IFR can be established in space environment and the long-range propagation of other REB pulses may be achieved.

    Fig.13. Reduction of beam quality during long-rang propagation.

    6. Conclusion

    In this study, we mainly focus on the establishment and influence of IFR during propagation of the REB in neutral gas and further explore the possibility of long-range propagation.A PIC-MCC method is used to calculate charge neutralization factorfeand current neutralization factorfm. The result shows that afterfesatisfies the Budker condition, transverse dimension of the beam head which has already expanded will be reduced to the initial size and form the beam neck. The beam may even be over-pinched in the connection area between beam neck and beam body. In the beam body area,the beam can be totally charge neutralized, which indicates that the IFR has been established. The PIC-MCC simulation results show that the beam is well focused by IFR in the beam body. Compared with the propagation of the REB in the natural plasma environment,longitudinal velocities of beam electrons will not oscillate periodically when the beam propagates in the self-generated plasma in neutral gas. However,effect of electron avalanche exacerbates gradually with time after the beam has been totally charge neutralization, leading to a significant current neutralization at the beam tail. Macro instabilities may also occur and interrupt efficient beam propagation. Decrease of the pulse width of the REB can reduce the influence of electron avalanche,whereas the beam quality during long-range propagation may be improved by using a train of REB pulses. Further improvement of beam quality during long-range propagation will be the next step in our study.

    In accordance with the results in the present study,it can be concluded that the IFR can effectively limit the transverse and longitude expansion of the REB and improve the transport efficiency during long-range propagation.

    Acknowledgement

    Project supported by the Joint Funds of the National Natural Science Foundation of China (Grant Nos. 61372050 and U1730247).

    猜你喜歡
    張芳趙強(qiáng)
    打電話
    雜文選刊(2023年10期)2023-10-12 21:08:09
    打電話
    Physical properties of relativistic electron beam during long-range propagation in space plasma environment?
    贏棋有獎(jiǎng)
    花鳥剪紙
    Seasonal phenology of the heterotrophic dino flagellate Noctiluca scintillans (Macartney) in Jiaozhou Bay and adjacent coastal Yellow Sea, China*
    趙強(qiáng)
    另類碰瓷
    東方劍(2018年11期)2018-02-20 07:56:28
    打電話
    小說月刊(2017年8期)2017-08-16 20:33:29
    我愿做“頂缸爸爸”,父母卻要她打胎才結(jié)婚
    一区二区三区精品91| 亚洲,欧美精品.| 日日撸夜夜添| 免费观看av网站的网址| 巨乳人妻的诱惑在线观看| 国产无遮挡羞羞视频在线观看| 建设人人有责人人尽责人人享有的| 侵犯人妻中文字幕一二三四区| 国产精品嫩草影院av在线观看| 91精品国产国语对白视频| 啦啦啦中文免费视频观看日本| 国产片内射在线| 一级毛片电影观看| 中文精品一卡2卡3卡4更新| 国产黄频视频在线观看| 捣出白浆h1v1| 侵犯人妻中文字幕一二三四区| 看免费av毛片| 日本猛色少妇xxxxx猛交久久| 亚洲色图 男人天堂 中文字幕| 久久久精品国产亚洲av高清涩受| 国产男人的电影天堂91| www.精华液| 最近最新中文字幕免费大全7| 两性夫妻黄色片| 日韩大码丰满熟妇| 亚洲av日韩在线播放| 一区二区日韩欧美中文字幕| 亚洲精品中文字幕在线视频| 日韩欧美精品免费久久| 又粗又硬又长又爽又黄的视频| 69精品国产乱码久久久| 亚洲欧美一区二区三区黑人| 国产福利在线免费观看视频| 欧美日韩亚洲综合一区二区三区_| 国产成人精品无人区| 午夜日本视频在线| 黑丝袜美女国产一区| 国产精品国产av在线观看| 午夜免费鲁丝| 天堂8中文在线网| 一级片'在线观看视频| 亚洲av国产av综合av卡| 国产精品香港三级国产av潘金莲 | 亚洲国产看品久久| 一二三四在线观看免费中文在| 啦啦啦在线免费观看视频4| 亚洲国产av影院在线观看| 精品酒店卫生间| 九九爱精品视频在线观看| 卡戴珊不雅视频在线播放| 十分钟在线观看高清视频www| 亚洲精品在线美女| 日日爽夜夜爽网站| 色婷婷av一区二区三区视频| 国产av码专区亚洲av| 免费在线观看视频国产中文字幕亚洲 | 日本av免费视频播放| 成人手机av| 日本wwww免费看| 日韩一卡2卡3卡4卡2021年| 成人国产麻豆网| 国产亚洲午夜精品一区二区久久| 两个人看的免费小视频| av国产久精品久网站免费入址| 无限看片的www在线观看| 丝袜脚勾引网站| 妹子高潮喷水视频| 女的被弄到高潮叫床怎么办| 秋霞伦理黄片| 亚洲精品中文字幕在线视频| 日日撸夜夜添| 亚洲天堂av无毛| 麻豆av在线久日| 国产精品av久久久久免费| 高清视频免费观看一区二区| 国产成人午夜福利电影在线观看| 久久亚洲国产成人精品v| 久久精品aⅴ一区二区三区四区| 亚洲专区中文字幕在线 | 99久久人妻综合| 另类精品久久| 韩国高清视频一区二区三区| 无限看片的www在线观看| 黄色视频在线播放观看不卡| 日韩成人av中文字幕在线观看| 亚洲国产av新网站| 日韩制服骚丝袜av| 操美女的视频在线观看| 欧美另类一区| 亚洲国产av新网站| 操出白浆在线播放| 日本一区二区免费在线视频| av在线观看视频网站免费| 男女之事视频高清在线观看 | 在线观看一区二区三区激情| 日日啪夜夜爽| 狠狠精品人妻久久久久久综合| 超碰成人久久| 亚洲欧美一区二区三区久久| 国产男女内射视频| 久久久精品国产亚洲av高清涩受| 日韩不卡一区二区三区视频在线| 精品国产乱码久久久久久小说| 日韩 欧美 亚洲 中文字幕| 中文字幕制服av| 欧美日韩av久久| 老司机亚洲免费影院| 精品少妇内射三级| 国产99久久九九免费精品| 国产日韩一区二区三区精品不卡| 人人妻人人澡人人看| 女性生殖器流出的白浆| 大香蕉久久成人网| 国产精品二区激情视频| 成年动漫av网址| 伊人久久大香线蕉亚洲五| 婷婷成人精品国产| 国产日韩欧美视频二区| 老司机深夜福利视频在线观看 | 天堂8中文在线网| 久久久久精品久久久久真实原创| 欧美xxⅹ黑人| 日韩一区二区视频免费看| e午夜精品久久久久久久| 黄片无遮挡物在线观看| 男女免费视频国产| 人体艺术视频欧美日本| 深夜精品福利| 国产在线一区二区三区精| 丝瓜视频免费看黄片| 如何舔出高潮| 欧美亚洲日本最大视频资源| 久久精品久久久久久久性| 日韩 亚洲 欧美在线| 99久久综合免费| 美女扒开内裤让男人捅视频| 欧美精品av麻豆av| 久久热在线av| 国产片内射在线| 2018国产大陆天天弄谢| 最新的欧美精品一区二区| 久久精品国产亚洲av涩爱| av视频免费观看在线观看| 丝瓜视频免费看黄片| 国产熟女午夜一区二区三区| 你懂的网址亚洲精品在线观看| 老司机影院毛片| 新久久久久国产一级毛片| 亚洲av成人不卡在线观看播放网 | 成年人午夜在线观看视频| 久久午夜综合久久蜜桃| 99久久精品国产亚洲精品| 精品国产一区二区久久| 成人国产麻豆网| 免费高清在线观看视频在线观看| 精品国产乱码久久久久久男人| 精品人妻熟女毛片av久久网站| 亚洲人成网站在线观看播放| 亚洲国产欧美日韩在线播放| 黄色 视频免费看| 亚洲自偷自拍图片 自拍| 天美传媒精品一区二区| 亚洲av福利一区| 久久久久精品久久久久真实原创| 久久精品亚洲av国产电影网| 国产精品av久久久久免费| 五月开心婷婷网| 久久毛片免费看一区二区三区| 欧美97在线视频| 欧美精品高潮呻吟av久久| 99久久人妻综合| 日韩中文字幕欧美一区二区 | 满18在线观看网站| 晚上一个人看的免费电影| 亚洲熟女毛片儿| 午夜福利影视在线免费观看| 国产人伦9x9x在线观看| 99热国产这里只有精品6| 国产极品天堂在线| 久久人妻熟女aⅴ| 国产探花极品一区二区| 国产精品国产av在线观看| 天天躁狠狠躁夜夜躁狠狠躁| 人妻人人澡人人爽人人| 久久国产亚洲av麻豆专区| 久久久国产一区二区| 日韩制服骚丝袜av| 男人操女人黄网站| 日韩av在线免费看完整版不卡| 欧美激情高清一区二区三区 | 久久久久视频综合| 看非洲黑人一级黄片| 汤姆久久久久久久影院中文字幕| 午夜福利网站1000一区二区三区| 成人亚洲欧美一区二区av| 你懂的网址亚洲精品在线观看| 在线精品无人区一区二区三| 男女无遮挡免费网站观看| 天天躁夜夜躁狠狠躁躁| 如何舔出高潮| 中文字幕亚洲精品专区| 国产免费视频播放在线视频| 夫妻午夜视频| 欧美日韩亚洲国产一区二区在线观看 | 色网站视频免费| 午夜福利一区二区在线看| 亚洲精品一区蜜桃| videos熟女内射| 久热爱精品视频在线9| 午夜福利影视在线免费观看| 9色porny在线观看| 国产人伦9x9x在线观看| 免费av中文字幕在线| 久久久久精品人妻al黑| 午夜免费男女啪啪视频观看| 尾随美女入室| h视频一区二区三区| 涩涩av久久男人的天堂| 丝袜在线中文字幕| 街头女战士在线观看网站| 母亲3免费完整高清在线观看| 日韩av不卡免费在线播放| 大话2 男鬼变身卡| 亚洲欧美一区二区三区黑人| 制服丝袜香蕉在线| 天天操日日干夜夜撸| 国产精品 国内视频| 汤姆久久久久久久影院中文字幕| 亚洲,欧美精品.| 波多野结衣av一区二区av| 亚洲国产毛片av蜜桃av| 国产午夜精品一二区理论片| 国产精品香港三级国产av潘金莲 | 久久综合国产亚洲精品| 天堂8中文在线网| www.熟女人妻精品国产| 各种免费的搞黄视频| 亚洲免费av在线视频| 日本午夜av视频| 欧美变态另类bdsm刘玥| 80岁老熟妇乱子伦牲交| 搡老岳熟女国产| 精品一区二区三卡| 国产日韩欧美在线精品| 在线观看人妻少妇| 亚洲少妇的诱惑av| 免费少妇av软件| 国产一区二区三区综合在线观看| 精品国产乱码久久久久久男人| 18在线观看网站| 久久久久精品久久久久真实原创| 可以免费在线观看a视频的电影网站 | 永久免费av网站大全| 永久免费av网站大全| 亚洲成色77777| 夫妻性生交免费视频一级片| 国产在线视频一区二区| 国产日韩欧美在线精品| 久久av网站| 国产av精品麻豆| 亚洲av国产av综合av卡| av.在线天堂| 日韩,欧美,国产一区二区三区| 一区二区三区乱码不卡18| 一级毛片电影观看| 2018国产大陆天天弄谢| 看免费av毛片| 国产成人av激情在线播放| 亚洲国产最新在线播放| 观看美女的网站| 久久青草综合色| 亚洲婷婷狠狠爱综合网| 97人妻天天添夜夜摸| 18禁动态无遮挡网站| 新久久久久国产一级毛片| 亚洲精品久久成人aⅴ小说| 男人添女人高潮全过程视频| 欧美黑人精品巨大| 国产精品.久久久| 美女高潮到喷水免费观看| 精品免费久久久久久久清纯 | 纯流量卡能插随身wifi吗| 黄色 视频免费看| 久久国产精品男人的天堂亚洲| 久久精品国产亚洲av涩爱| 黄色一级大片看看| 久久99一区二区三区| 少妇人妻久久综合中文| 久久青草综合色| 久久久久久免费高清国产稀缺| 国产深夜福利视频在线观看| 午夜福利乱码中文字幕| 亚洲成人av在线免费| 亚洲精品久久久久久婷婷小说| 亚洲av欧美aⅴ国产| 99热国产这里只有精品6| 男女国产视频网站| av电影中文网址| 欧美日韩亚洲高清精品| 亚洲人成网站在线观看播放| 最近最新中文字幕大全免费视频 | 成人亚洲精品一区在线观看| 久久99精品国语久久久| 性少妇av在线| 国产精品 欧美亚洲| 在线观看国产h片| 9191精品国产免费久久| 国产精品国产三级国产专区5o| 亚洲av成人不卡在线观看播放网 | 国产一区有黄有色的免费视频| 国产av码专区亚洲av| 国产精品成人在线| 日本vs欧美在线观看视频| 激情五月婷婷亚洲| 亚洲欧美清纯卡通| 欧美日韩亚洲高清精品| 考比视频在线观看| 最近中文字幕2019免费版| av免费观看日本| 一边亲一边摸免费视频| 热re99久久国产66热| 国产无遮挡羞羞视频在线观看| 在线观看人妻少妇| 久久鲁丝午夜福利片| 国产麻豆69| 亚洲精品av麻豆狂野| a 毛片基地| 男女午夜视频在线观看| 欧美激情极品国产一区二区三区| 日韩大码丰满熟妇| 少妇人妻久久综合中文| 久久久久久久久免费视频了| 一本大道久久a久久精品| 久久久久久久国产电影| 久热这里只有精品99| 中文字幕人妻丝袜制服| 一区二区av电影网| 国产片内射在线| 在现免费观看毛片| 伊人久久国产一区二区| 热99久久久久精品小说推荐| 在线天堂最新版资源| 亚洲av成人精品一二三区| 综合色丁香网| 少妇人妻精品综合一区二区| 老熟女久久久| 国产一区二区激情短视频 | 七月丁香在线播放| 看免费av毛片| 男人爽女人下面视频在线观看| 999精品在线视频| 午夜福利乱码中文字幕| 精品国产乱码久久久久久男人| 久久精品久久久久久噜噜老黄| 看非洲黑人一级黄片| 丰满迷人的少妇在线观看| 亚洲精品日本国产第一区| 91成人精品电影| 性高湖久久久久久久久免费观看| 这个男人来自地球电影免费观看 | a级毛片黄视频| 九九爱精品视频在线观看| 在线精品无人区一区二区三| 亚洲国产欧美一区二区综合| 婷婷色综合大香蕉| 日本爱情动作片www.在线观看| 美国免费a级毛片| 亚洲国产成人一精品久久久| 亚洲av日韩精品久久久久久密 | 亚洲欧美精品综合一区二区三区| 国产xxxxx性猛交| 欧美日韩亚洲综合一区二区三区_| 亚洲熟女精品中文字幕| 黄色怎么调成土黄色| 久久久久国产精品人妻一区二区| 天天影视国产精品| 成人亚洲欧美一区二区av| 久久久精品国产亚洲av高清涩受| 午夜免费鲁丝| 老司机在亚洲福利影院| 男人添女人高潮全过程视频| 亚洲激情五月婷婷啪啪| 在线观看国产h片| 亚洲精品中文字幕在线视频| 看免费成人av毛片| 一本色道久久久久久精品综合| 亚洲欧美中文字幕日韩二区| 久久久欧美国产精品| 日韩人妻精品一区2区三区| 黄色怎么调成土黄色| 男人舔女人的私密视频| 国产有黄有色有爽视频| 国产亚洲最大av| 欧美 日韩 精品 国产| 美女脱内裤让男人舔精品视频| 女人精品久久久久毛片| 老司机亚洲免费影院| av视频免费观看在线观看| 91aial.com中文字幕在线观看| 久久97久久精品| 女人高潮潮喷娇喘18禁视频| 久久人妻熟女aⅴ| 国产精品三级大全| 如何舔出高潮| 精品免费久久久久久久清纯 | 亚洲av欧美aⅴ国产| 丰满乱子伦码专区| 爱豆传媒免费全集在线观看| 超碰成人久久| 午夜免费男女啪啪视频观看| 男女边摸边吃奶| 国产成人免费观看mmmm| 国产精品一二三区在线看| 国产黄色视频一区二区在线观看| 精品人妻一区二区三区麻豆| 成年美女黄网站色视频大全免费| 一本色道久久久久久精品综合| 老司机深夜福利视频在线观看 | 51午夜福利影视在线观看| 男女之事视频高清在线观看 | av视频免费观看在线观看| 婷婷色麻豆天堂久久| 亚洲欧美一区二区三区久久| 一级爰片在线观看| 久久免费观看电影| 亚洲精品久久午夜乱码| 男人舔女人的私密视频| 国产免费现黄频在线看| 少妇人妻 视频| 丝袜脚勾引网站| 亚洲成人国产一区在线观看 | 国产熟女欧美一区二区| 国产精品秋霞免费鲁丝片| 丝袜美足系列| 伊人久久国产一区二区| 国产黄色免费在线视频| 黑人猛操日本美女一级片| 国产无遮挡羞羞视频在线观看| 黄色一级大片看看| 日本av免费视频播放| 国产人伦9x9x在线观看| 交换朋友夫妻互换小说| 国产亚洲欧美精品永久| av不卡在线播放| 亚洲国产看品久久| 久久天躁狠狠躁夜夜2o2o | 天天操日日干夜夜撸| 啦啦啦在线观看免费高清www| 九色亚洲精品在线播放| 考比视频在线观看| 国产一卡二卡三卡精品 | 久久久国产一区二区| 国产精品一区二区精品视频观看| 国产不卡av网站在线观看| 青草久久国产| 亚洲综合精品二区| 国产成人系列免费观看| 色精品久久人妻99蜜桃| 中文字幕高清在线视频| 国产又色又爽无遮挡免| 街头女战士在线观看网站| 国产精品女同一区二区软件| 国产在线一区二区三区精| 在线天堂最新版资源| 免费久久久久久久精品成人欧美视频| 亚洲精品在线美女| 亚洲国产精品一区二区三区在线| 人人澡人人妻人| 男女午夜视频在线观看| 黑丝袜美女国产一区| 日本一区二区免费在线视频| 校园人妻丝袜中文字幕| 中文字幕亚洲精品专区| 久久久国产精品麻豆| 亚洲伊人久久精品综合| 老司机深夜福利视频在线观看 | 夫妻性生交免费视频一级片| 又大又黄又爽视频免费| 国产精品秋霞免费鲁丝片| 精品少妇一区二区三区视频日本电影 | 日本猛色少妇xxxxx猛交久久| 黄网站色视频无遮挡免费观看| 夫妻性生交免费视频一级片| bbb黄色大片| 欧美日韩综合久久久久久| 日韩大码丰满熟妇| 亚洲av欧美aⅴ国产| 2018国产大陆天天弄谢| 高清av免费在线| 久久人人97超碰香蕉20202| 19禁男女啪啪无遮挡网站| 女性生殖器流出的白浆| 99精国产麻豆久久婷婷| 亚洲天堂av无毛| 成人18禁高潮啪啪吃奶动态图| 天天躁夜夜躁狠狠躁躁| 欧美在线一区亚洲| 午夜日韩欧美国产| 不卡视频在线观看欧美| 日韩中文字幕视频在线看片| 亚洲av综合色区一区| 美女高潮到喷水免费观看| 亚洲美女视频黄频| 国产黄色视频一区二区在线观看| 久久久久精品人妻al黑| 在线观看免费午夜福利视频| 水蜜桃什么品种好| www.av在线官网国产| 久久久久视频综合| 水蜜桃什么品种好| 亚洲精品美女久久av网站| 可以免费在线观看a视频的电影网站 | 国产精品久久久久成人av| 超碰97精品在线观看| 两性夫妻黄色片| 国产成人欧美| 久久精品久久久久久噜噜老黄| 乱人伦中国视频| 最新的欧美精品一区二区| av电影中文网址| 精品视频人人做人人爽| 亚洲成国产人片在线观看| 在线亚洲精品国产二区图片欧美| 亚洲精品,欧美精品| 丝袜脚勾引网站| 欧美日韩一区二区视频在线观看视频在线| 丝袜美足系列| 十八禁网站网址无遮挡| 中文字幕最新亚洲高清| 涩涩av久久男人的天堂| 亚洲美女黄色视频免费看| 婷婷色av中文字幕| 国产精品成人在线| 水蜜桃什么品种好| 在线天堂中文资源库| 久久影院123| 男女午夜视频在线观看| 国产一区二区三区av在线| 中文乱码字字幕精品一区二区三区| 国产日韩欧美在线精品| 这个男人来自地球电影免费观看 | 午夜福利影视在线免费观看| 成人影院久久| 免费黄色在线免费观看| 一级毛片黄色毛片免费观看视频| 国产激情久久老熟女| 久久久久久人人人人人| 人人妻人人添人人爽欧美一区卜| 亚洲欧洲国产日韩| 中文字幕亚洲精品专区| 男人操女人黄网站| 久久久精品区二区三区| 在线观看免费日韩欧美大片| 亚洲欧美日韩另类电影网站| 日韩人妻精品一区2区三区| 午夜激情av网站| 又大又爽又粗| 久久精品久久久久久久性| 欧美精品一区二区免费开放| 建设人人有责人人尽责人人享有的| 1024视频免费在线观看| 九色亚洲精品在线播放| 老鸭窝网址在线观看| 一二三四在线观看免费中文在| 久久精品国产亚洲av高清一级| 欧美 亚洲 国产 日韩一| 欧美日韩亚洲综合一区二区三区_| 国产免费视频播放在线视频| 在线精品无人区一区二区三| 色视频在线一区二区三区| 亚洲第一青青草原| 国产亚洲最大av| 蜜桃在线观看..| 中文字幕最新亚洲高清| 1024香蕉在线观看| 黑人猛操日本美女一级片| 一区二区日韩欧美中文字幕| 亚洲国产中文字幕在线视频| 成年动漫av网址| 国产激情久久老熟女| 一区二区日韩欧美中文字幕| √禁漫天堂资源中文www| 在线观看免费视频网站a站| 高清av免费在线| 欧美黑人精品巨大| 国产亚洲午夜精品一区二区久久| 国产色婷婷99| 一边亲一边摸免费视频| 亚洲精品久久成人aⅴ小说| 国产精品欧美亚洲77777| 国产精品久久久人人做人人爽| √禁漫天堂资源中文www| 午夜福利一区二区在线看| 一级毛片电影观看| 国产不卡av网站在线观看| 国产精品一区二区在线不卡| 蜜桃在线观看..| 两个人看的免费小视频| 国产伦理片在线播放av一区| 国产一区二区 视频在线| 十八禁网站网址无遮挡| 亚洲色图综合在线观看| av有码第一页| 丰满少妇做爰视频| 精品福利永久在线观看| 老司机影院成人| 国产成人精品无人区| 在线免费观看不下载黄p国产| 亚洲美女黄色视频免费看| 老司机靠b影院| 一区二区av电影网| 在线观看人妻少妇| 中文精品一卡2卡3卡4更新| 亚洲成av片中文字幕在线观看|