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

    Effects of fast ions produced by ICRF heating on the pressure at EAST

    2020-03-09 13:22:02ZhenZHENG鄭振NongXIANG項農ChengYANG楊程andYingfengXU徐穎峰
    Plasma Science and Technology 2020年2期

    Zhen ZHENG(鄭振),Nong XIANG(項農),Cheng YANG(楊程) and Yingfeng XU (徐穎峰)

    1 Institute of Plasma Physics, Chinese Academy of Sciences, Hefei 230031, People’s Republic of China

    2 University of Science and Technology of China, Hefei 230026, People’s Republic of China

    3 Center for Magnetic Fusion Theory, Chinese Academy of Sciences, Hefei 230031, People’s Republic of China

    Abstract Ion cyclotron resonance heating (ICRH), which can produce fast ions, is an important auxiliary heating method at EAST. To analyze the effect of ICRH-induced fast ions on the plasma pressure at EAST,simulations are performed using TRANSP and TORIC codes.It is found that the ICRF-induced fast ion pressure cannot be negligible when the ICRF power is sufficiently high.The magnitude of the total ion pressure can be raised up to 60%of the total pressure as the input power rises above 3 MW. The pressure profile is also significantly modified when the resonant layer is changed.It is shown that by changing the wave frequency and antenna position,the total ion pressure profile can be broadened,which might provide an option for profile control at EAST.

    Keywords: ICRH, equilibrium, fast ion, pressure, minority heating

    1. Introduction

    Plasma equilibrium [1-5], which is comprised of the plasma geometry and important profiles of the plasma pressure,current,and safety factor, plays an important role in analyzing the experimental data, magneto-hydrodynamic (MHD) activities,and transport process. The fast ions produced by neutral beam injection(NBI)and ion cyclotron resonance heating(ICRH)can trigger many interesting physics, e.g. edge-localized modes [6]and Alfvén eigenmodes [7]. In addition, they can contribute to the plasma pressure greatly and affect the plasma equilibrium.There are several methods and numerical codes developed for equilibrium reconstruction. In this paper, the equilibrium is obtained using EFIT code [1, 2].

    Recently, a variety of high-performance experiments have been successfully performed on the EAST tokamak[6-10]. These experiments aimed to explore how the plasma would behave under the high-performance conditions. ICRH is one of important auxiliary heating methods at EAST,which can produce considerable fast ions.Fast ions can contribute to the plasma pressure at EAST, which can influence the MHD activities and plasma confinement. The effects of fast ions produced by NBI on the EAST equilibrium have been discussed in the previous work [11]. In this work, we discuss how ICRH-induced fast ions affect the plasma pressure when high ICRH power is used at EAST.

    The injection of waves in the ion cyclotron range of frequencies is a well-established method of heating and current drive in the tokamak. At EAST, the minority heating scheme is chosen to heat ions. The wave power is mainly absorbed through collisionless wave-particle interactions by hydrogen ions,which are then slowed down via electron drag or collisions with the background ions. More recently, much measured information related to ICRH has been obtained on EAST, indicating the existence of a considerable number of fast ions [12, 13]. To analyze the diagnostic data and understand the physical processes involved, a reliable equilibrium is needed which should include contribution of fast ions. In this paper, the TRANSP /TORIC [14, 15] codes are adopted to calculate the power deposition during ICRF and the resulting fast ion pressure. To demonstrate the fast ion effects on plasma pressure, the dependence of the fast ion pressure on the input power is discussed.

    It is well known that the bootstrap current strongly depends on the plasma pressure gradient, and many MHD instabilities such as the neoclassical tearing mode(NTM)and ballooning modes are determined by the bootstrap current gradient. The importance of the bootstrap current for tokamaks has already been pointed out[16,17].It has been shown that a high bootstrap current would improve prospects for energy production[18].Moreover,a high bootstrap fraction is beneficial for the high-performance phase of ITER and CFETR[19,20].Therefore,whether the plasma pressure can be modified by fast ions is of interest. It is shown that the pressure profile can be modified by changing the position of the power deposition and the fast ion pressure, ultimately promoting the bootstrap and resulting in hollow current and reverse q profile.To achieve this aim,we simulate the effects of different antenna positions, resonance layer positions, and ICRF powers on the plasma equilibrium.

    This paper is organized as follows. The numerical tools,including TRANSP and TORIC, are introduced in section 2.Section 3 shows the simulation results of the effects of fast ions induced by ICRF waves on the dependence of the RF power, RF frequency, and antenna positions. The summary and discussion are presented in section 4.

    2. Physical model and simulation setup

    The equilibrium used in this paper is based on the previous work on kinetic equilibrium reconstruction for the H-mode plasma on the EAST tokamak [11]. The electron and ion temperature profiles, as a function of the square root of the normalized toroidal flux, ρ, are shown in figure 1(a). The electron density profile is shown in figure 1(b). For this discharge, the toroidal magnetic field is about 2.25 T.The plasma consists of hydrogen(H),deuterium(D),and impurity ions.The H ion concentration is 3.4%. The equilibrium and fundamental cyclotron resonance for H and the second harmonic for D ions with the wave frequency of 34 MHz are plotted in figure 1(c). There are two ICRF antennas located at windows B and I, whose frequencies andφn (the toroidal mode number) spectrum are nearly the same. They are both in the equatorial plane. Since the fast ion pressure is our primary concern and the fast ion pressure produced by the ICRF heating mainly depends on the input power,wave frequency and antenna position as shown later.We assume that there is one ICRH antenna launching a wave with the toroidal mode numberφn = 23 for the sake of simplicity.

    The TRANSP code has been applied to compute fast ion distribution functions. It requires the plasma equilibrium and kinetic profiles. TRANSP then applies an equilibrium solver taking the plasma pressure and bootstrap current into account.When TRANSP calculates the total fast ion pressure, a fullwave code is needed to analyze physical quantities related to fast ions,such as the wave distribution and power deposition.

    Figure 1.Measured and fitted profiles of (a) electron density, (b)electron temperature, and ion temperature. (c) Position of the main resonance layer for fRF = 34 MHz on EAST; n is the harmonic number of ICRH.

    Figure 2.2D profiles of power density to (a) minority ions and (b)bulk ions,in W/cm3.(c)Power density to electrons absorbed by the fast wave and IBW.

    TORIC has been widely used to study the ion cyclotron wave propagation and power deposition in the tokamak plasma[21-24].The physics for ICRF are well described in the TORIC solver including the wave propagation, mode conversion, and power deposition in a tokamak.It has been successfully applied to the design and modeling work on the ICRF heating scenarios on EAST [25, 26]. The Maxwell’s wave equation for ICRF in arbitrary axisymmetric toroidal geometry,including the tokamak equilibrium, can be solved in TORIC:

    Figure 3.(a) Hydrogen ion energy distribution function (normalized) at ρ = 0.025. (b) Effective temperature.

    where E is the wave electric feild,represents the antenna current,andis the response current of the plasma derived by solving the linear Vlasov equation based on the fniite Larmor radius assumption [15, 27].

    In TRANSP,the FPPRF[28]module computes the minority ion phase space distribution. This module uses the bouncedaveraged Fokker Planck equationin whichμ is the magnetic moment,C is the bounced-averaged and linearized collision operator, Q is the complete bounce-averaged quasilinear operator, S represents sources and sinks, and R is a radial transport operator.

    In this work, we focus on the effect of the fast ion produced by ICRF heating on the plasma pressure. We choose the equilibrium of #62585@3.8 s because of the good measurements of temperature and density profiles, which can be used for the equilibrium.We use TORIC in TRANSP to give out the wave distribution and power deposition; then the fast ion pressure can be obtained by FPPRF.

    3. Simulation results

    3.1. Fast ion pressure for 1.6 MW RF power

    Firstly, we simulate the case with fRF= 34 MHz and 1.6 MW of RF power. The cyclotron resonance of the hydrogen ions is located at R = 1.9 m. Figures 2(a) and (b) illustrate the power deposition to the minority ions and deuteron ions.The width of the resonance region observed in figures 2(a) and (b) is conditioned by the Doppler shift from equationThe simulation shows that most of the power is deposited in the core region. About 77.7% power is transferred to the minority ions via H fundamental cyclotron resonance. The fraction of power directly absorbed by the electrons via transit time magnetic pumping and Landau damping is 7.2%.The power to deuteron ions is about 11.2% via D harmonic cyclotron resonance.The remaining 3.9%of the power goes to the impurities.It can be seen that the electron absorption is very weak.Figure 2(c) indicates that most of the power damping by electrons is from the fast wave.Nevertheless,the ion heating is very effective, as about 90% of total absorbed power goes to ions,and the minority ion heating dominates,as the minority ions are effectively heated in the core region.This can be seen clearly in figure 3.

    The heating of ICRF dramatically increases the perpendicular energy of minority ions through the collisionless,wave-particle resonant interaction in the resonance region,which will cause a part of the minority ions to become anisotropic. As a result, a large minority H tail is formed (see figure 3(a)) and fast ions are generated. These fast ions deliver their energy to the bulk ions via collision, heating the plasma.The ions with energies over the critical energy [29] =Ecrabout 33 keV, will primarily transfer their energy to electrons.Fast ions will contribute to the total pressure, which is what we are concerned with in this paper.This tail temperature can be estimated by Stix’s isotropic model [30]. In this model, when E approaches infinity, Teffapproaches the asymptotic value127 keV.This shows a good agreement between the theoretical and simulation models when E < 14 keV (see figure 3(b)). The effective tail temperatures for the two models are different but have the same order when E > 14 keV;a possible reason is that the theoretical curve is from a simple isotropic model,which is less precise than the anisotropic model used in the simulations as mentioned in section 2.

    Once the H distribution function is obtained, the fast ion pressure is then determined. As can be seen from figure 4,the fast ion pressure peak appears in the core region where the power is mainly deposited.The contribution of fast ions to the total pressure in the core region is about 40%,while it is less than 5%in the edge region.Clearly,the pressure is unaffected for ρ > 0.4. Such a weak effect of fast ions on the total pressure might be due to the low input power.Since the ICRF system at EAST has the capability of a 12 MW power source,it is necessary to find out the fast ion contribution to the total pressure with higher input power.

    Figure 4.Fast ion pressure (red line), bulk ion (D) pressure (black line) and total pressure (magenta line).

    Figure 5.(a)Power density profiles with different RF powers.(b)Fast ion pressures,background pressures,and total pressures with different RF powers. Black lines represent the background pressures, and magenta lines represent the total pressures. Red, green, and blue lines represent the fast ion pressures at 3.2, 6.4, and 10 MW, respectively.

    3.2. Dependence of fast ion pressure on RF power

    The simulations are carried out with a scan of RF powers from 3.2 to 10 MW. The profiles of fast ion pressure and power density are plotted in figures 5(a) and (b). It can be found that the profiles of both pressure and power density increase with the increase in RF power, as expected. The powers are all deposited in the core region, and the background pressures are all the same. Then we examinenamely the ratio of the fast ion pressure to the background pressure. At the peaks,is 0.59, 0.85, and 1.0,respectively, for the cases in which the wave power is 3.2,6.4, and 10 MW. Clearly, theincreases with the RF input power.When the input power is larger than 3.2 MW,the maximum value ofis larger than 0.5, which means that the fast ion pressure is significant. The value ofincreases with the input power even in the outer region, for example at ρ = 0.5.decreases rapidly when moving outward to the plasma edge. The fast ion effects almost disappear for ρ > 0.6, even for 10 MW input power.

    The location of the peak of the fast ion pressure should be noted.The position of the power density in the core region is within ρ = 0.3, just as depicted in figure 2(a). The power density and the fast ion pressure are close from the axis.However, the fast ions have their own banana orbits, and a flux surface average should be done on the fast ion pressure first. Also, the pressure is independent of the volume. As a result, the locations of the peaks of the fast ion pressure and power density are different but still close to the axis, within 0.1 < ρ < 0.3.

    Figure 6.The hydrogen ion distribution function in velocity coordinate at ρ = 0.475 with 10 MW. Two energy levels, 0.25 and 0.5 MeV, are shown.

    The contours of hydrogen ion distribution function in the velocity space coordinate at ρ = 0.475 is shown in figure 6.We can see that the velocity distribution is anisotropic and that the contour has a‘rabbit-ear’shape.The‘ears’are inside the passing-trapped boundary. The ‘ear’ shape of the contours shown in figure 6 is similar to the previous results[31-33]. The anisotropy and ear-like features suggest that a large number of trapped particles are produced. The tips of‘ears’are related to the trapped particles whose banana tips lie on the resonance layer. The ‘ears’ are a result of the fast ion resonance localization. The ICRF wave gives primarily particles perpendicular energy through wave-particle interaction,and thus resonant passing particles are turned to trapped particles. The ICRF heating trends to produce energetic trapped particles with banana tips near the resonant layer.

    The simulation results of the on-axis heating indicate that the equilibrium in the core region is significantly modified as the input power exceeds 3.2 MW. This is because the cyclotron resonances are located in the core region where fast ions are produced.It can be expected that as the resonance is moved outward by changing the wave frequency and antenna position, off-axis heating might occur [34, 35]. As a result,the total pressures might be modified in the outer region,which is discussed in the following sections.

    3.3. Dependence of fast ion pressure on RF frequency

    The hydrogen resonance location along the major radius from the central region can be shifted to the low-field side by changing the RF wave frequency from 34 to 30 MHz with 1.6 MW RF power. Figure 7 shows the resonance layers of different frequencies.

    Figure 7.Resonance layers at 34 MHz (magenta), 32 MHz (green)and 30 MHz (red).

    Figure 8 shows the power deposition to minority ions for the wave frequencies of 32 and 30 MHz. Generally, the fast wave is launched from different positions of the antenna,propagating toward the core region and being enhanced. In the off-axis heating case, the wave meets the resonance layer first. The wave resonates with the minority ions in the mid-ρ region and deposits the power here.With the resonance layer moving to the low-field side,more power will be deposited in the mid-ρ region. Consequently, the fast ion pressure profile is significant only off-axis.

    The power density and fast ion pressure profiles for different incident frequencies are illustrated in figure 9. It is observed in figure 9(a)that as the frequency decreases from 34 to 32 MHz, the peaks of the power density profiles are still located in the core,but the peak of the latter is lower and more wave energy is deposited in the mid-ρ region.Accordingly,the peaks of the fast pressure profiles for the two cases appear in the core region while the fast pressure is broader for the lowerfrequency case. As the frequency becomes 30 MHz, a significant portion of the wave energy is deposited in the outer region, as shown in figure 8(b). As can be seen in figure 9(a),the peak of the power density profiles moves towards the mid-ρ region.Meanwhile the fast ion pressure peaks around ρ = 0.3 because the ion density increases rapidly with ρ.

    In summary,as the frequency decreases,more power will be deposited in the outer region due to the fact that the cyclotron resonance moves from the core region to the outer region. Accordingly, the peak of the fast ion pressure also moves outwards and the fast ion pressure is greatly broadened. Such modification of the fast ion pressure should be more significant with a higher input power.Therefore,we will illustrate how the ion pressure can be modified by varying the cyclotron resonances.

    We choose the case in which the frequency is 30 MHz,and increase the power to 6.4 MW to find out how the fast ion pressure will change with a higher power. The power density and fast ion pressure are shown in figure 10. It can be found that both the power density and the pressure increase dramatically. The pressure broadens in the region of ρ = 0.2-0.3,and the pressure in ρ = 0.7 also rises.Generally,the profile is broadened in the region of ρ = 0.1-0.8. It is as expected to broaden the profile off-axis. By choosing appropriate frequencies with more antennas, the fast ion pressure can be more important in the outer region.

    3.4. Dependence of fast ion pressure on antenna positions

    Figure 8.2D profiles of power density to minority ions for (a) 32 and (b) 30 MHz, in W/cm3.

    Figure 9.(a) Power density profiles at different frequencies. (b) Fast ion pressure profiles at different frequencies.

    Figure 10.(a) Power density and (b) fast ion pressure profile at the frequency of 30 MHz and power of 6.4 MW.

    Figure 11.Distribution of the electric field as the wave is launched from the antenna at the upper right side.

    As just shown, the fast ion pressure strongly relies on how the incident wave arrives at the resonant layer where the power is deposited.When the wave reaches the resonant layer at the outer region,the fast ion pressure could be much broadened as shown in figure 10(b), in which the resonance layer does not pass the core region.Thus we can imagine that for a given resonant layer,as shown in figure 1,the wave power could be absorbed in the outer region if the antenna is moved to the upper right side.Considering the geometries of the plasma and the antenna, the wave might not focus in the core region but propagate towards the off-axis region and then turn downward. Sometimes the wave reflects in the lower right side. During this process, the wave may penetrate the resonance layer once or several times,or never reach it, depending on where the wave is launched,φn ,and the location of the resonance layer.The footprint where the wave reaches the resonance layer is the power density region;as it is weak when the wave reaches the resonance layer at the second time, the power is mainly deposited in the first pass.Consequently, the power density and fast ion pressure profiles are expected off-axis. Figure 11 gives an example of how the wave propagates from the upper right side.

    Figure 12.Profiles with different antenna positions. Red curves represent antenna position I, blue curves represent position II, and green curves represent position III. (a) Three antenna positions chosen to launch the fast wave. (b) Power density profiles. (c) Fast ion pressure profiles.

    To study the effect of different antenna positions on the plasma, three positions (I, II, III) are set as depicted in figure 2(a). The RF wave frequency is 34 MHz and the RF power is 1.6 MW. Figures 12(b)-(c) illustrate the power density and the fast ion pressure.It can be found that the peak of the power density profile moves towards the outside and decreases as the antenna position changes from the middle plane to the upper side and the profile broadens in position III.The peak of the pressure profile is found to be related to that of the power density.In position I,the peaks of profiles are in the central region,within ρ < 0.4.As the antenna goes up to II, the fast ion pressure in the core region drops dramatically and the peak moves from ρ = 0 to ρ = 0.15. At ρ = 0.18,= 0.2.It can be noticed that the pressure rises and gets broader in the region 0.25 < ρ < 0.5. In position III, the profile becomes much lower but broadened and flattened. At ρ = 0.18,= 0.1. In the off-axis region, compared with position II, the profile rises little. In addition, all the profiles of the three positions are within ρ < 0.5. It is expected that the profile in the edge region will rise under a high-power injection, as mentioned before.

    We choose the case at antenna position III, with RF powers scanning from 1.6 to 10 MW,to find out how the fast ion pressure will change with a higher power.The results are shown in figure 13. In antenna position III, the results are similar to those in antenna position I. Clearly, the pressure increases with power for all the profiles. The power density profiles, compared with the results in figure 12(c), rise apparently in the edge region. At ρ = 0.1,is 0.09,0.15, 0.26, and 0.37, respectively, which is much lower but flatter than that in antenna position I case.At ρ = 0.5,is 0.05,0.13,0.32,and 0.5,respectively,which is close to the result in antenna position I.In addition, the profiles rise even in the edge region with ρ > 0.6. Generally speaking, the pressure profiles become flatter in the outer region.

    4. Summary and discussion

    In this work, with different RF frequencies, RF powers, and antenna positions,the contribution of fast ions induced by ICRF to the plasma pressure is studied using TRANSP ad TORIC codes. Firstly, we simulate the case with fRF= 34 MHz and 1.6 MW of RF power. The ion distribution function is also given,which shows that numerous fast ions are generated by the heating of ICRF. Next, we simulate on-axis minority heating cases with different RF powers,and it is found that the fast ion pressure in the edge region cannot be affected by the on-axis heating. When the power is over 3.2 MW for the on-axis heating, the contribution of fast ions induced by ICRF cannot be neglected. Then we simulate cases with different RF frequencies,and it is found that the fast ion pressure profile can be broadened through off-axis heating and, moreover, that the fast ion pressure in the edge region can be increased.Finally,by changing the antenna position,we find that the fast ion pressure can also be broadened and even made flatter as the antenna position moves up. We also discover how the fast ion pressure changes at different RF powers in the antenna position III.

    In our simulations, we use the RF power 10 MW. This is actually too high for EAST,and it may cause the fast ions to not be constrained by the magnetic field, eventually leading to fast ion losses to the first wall.In addition,the absorption efficiency of ICRF at EAST is about 40%-50%, and the technique of changing the antenna position is not feasible right now. However, the simulations still provide some reasonable predictions and give us some idea as to how the waves propagate from different launch locations and how this affects the profile of the fast ion pressure in the resonance region.This will be helpful,for example, in optimizing the antenna design in the future. Changing the frequency and the antenna position will broaden the peak of the pressure profile. By choosing appropriate RF frequencies and antenna positions, the purpose of controlling the fast ion pressure profile can be achieved,ultimately generating a broadened bootstrap current. A broadened bootstrap profile in the core region is beneficial to high-bootstrap and high-βNoperation in future devices such as ITER and CFETR.The width of the improved core confinement region should be large enough to avoid the development of a narrow domain with a too localized and very steep pressure gradient [36].

    Acknowledgments

    The authors would like to thank Doctors Guosheng Xu,Damao Yao, and Yongliang Li for helpful discussions.This work was supported by the National Key R&D Program of China (Nos. 2017YFE0300400 and 2017YFE0300406),and National Natural Science Foundation of China (Nos.11575239 and 11775265).

    Figure 13.(a)Power densities with different RF powers in antenna position III.(b)Fast ion pressure profiles(solid lines)and total pressure profiles (dash lines) at different RF powers in antenna position III; black line represents the background pressure.

    ORCID iDs

    Zhen ZHENG (鄭振)https://orcid.org/0000-0001-6866-1903

    国产片内射在线| 亚洲成色77777| 国产不卡av网站在线观看| 久久狼人影院| 成年av动漫网址| 纯流量卡能插随身wifi吗| 国产亚洲精品久久久com| 人人妻人人添人人爽欧美一区卜| 亚洲一级一片aⅴ在线观看| 天天躁夜夜躁狠狠久久av| 午夜福利,免费看| 国产免费一区二区三区四区乱码| 国产欧美另类精品又又久久亚洲欧美| 最新的欧美精品一区二区| 亚洲人成77777在线视频| 午夜免费观看性视频| 久久久久久久久久久丰满| 另类精品久久| 国产黄色视频一区二区在线观看| 国精品久久久久久国模美| 亚洲国产欧美在线一区| 一区在线观看完整版| 一级毛片黄色毛片免费观看视频| 免费看不卡的av| 99九九在线精品视频| 男女啪啪激烈高潮av片| 极品少妇高潮喷水抽搐| 狂野欧美激情性xxxx在线观看| 免费看光身美女| 曰老女人黄片| 搡女人真爽免费视频火全软件| av专区在线播放| 热re99久久国产66热| 成人毛片60女人毛片免费| 热re99久久精品国产66热6| 亚洲国产精品专区欧美| 中国美白少妇内射xxxbb| 91精品国产国语对白视频| 22中文网久久字幕| 欧美成人精品欧美一级黄| 国产欧美另类精品又又久久亚洲欧美| 亚洲人成77777在线视频| 最近中文字幕2019免费版| 久久久欧美国产精品| 在线亚洲精品国产二区图片欧美 | 亚洲精品456在线播放app| 日本欧美国产在线视频| 啦啦啦中文免费视频观看日本| 亚洲内射少妇av| www.色视频.com| 午夜日本视频在线| 亚洲经典国产精华液单| 日韩成人伦理影院| 久久国产精品男人的天堂亚洲 | av女优亚洲男人天堂| 久久韩国三级中文字幕| 丰满迷人的少妇在线观看| 大码成人一级视频| 国产极品粉嫩免费观看在线 | 大码成人一级视频| 久久精品久久久久久噜噜老黄| 日产精品乱码卡一卡2卡三| 男的添女的下面高潮视频| 肉色欧美久久久久久久蜜桃| www.av在线官网国产| 一级毛片aaaaaa免费看小| 国产日韩欧美亚洲二区| 国产深夜福利视频在线观看| 九九爱精品视频在线观看| 三上悠亚av全集在线观看| 国产欧美日韩一区二区三区在线 | 制服丝袜香蕉在线| 国产成人精品在线电影| 日韩熟女老妇一区二区性免费视频| 国产成人精品福利久久| 免费大片黄手机在线观看| 亚洲无线观看免费| 色5月婷婷丁香| 日韩一区二区视频免费看| 久久狼人影院| 日韩电影二区| 天天影视国产精品| 国产爽快片一区二区三区| 亚洲欧美日韩另类电影网站| 男人爽女人下面视频在线观看| 精品人妻熟女av久视频| 日韩 亚洲 欧美在线| 亚洲欧美中文字幕日韩二区| 一区二区三区免费毛片| 成人无遮挡网站| 在线观看免费高清a一片| 精品卡一卡二卡四卡免费| 午夜老司机福利剧场| 王馨瑶露胸无遮挡在线观看| 久久久久久久精品精品| 亚洲欧洲精品一区二区精品久久久 | 18禁动态无遮挡网站| 亚洲欧洲精品一区二区精品久久久 | 高清视频免费观看一区二区| 秋霞在线观看毛片| 国产片内射在线| 视频在线观看一区二区三区| 少妇高潮的动态图| 久久鲁丝午夜福利片| 国产成人一区二区在线| a级片在线免费高清观看视频| 18+在线观看网站| h视频一区二区三区| 欧美97在线视频| 日韩视频在线欧美| 国产高清三级在线| av黄色大香蕉| 精品亚洲乱码少妇综合久久| 国产在线一区二区三区精| 亚洲美女搞黄在线观看| 亚洲美女搞黄在线观看| 精品亚洲成a人片在线观看| 亚洲少妇的诱惑av| 国产乱人偷精品视频| av黄色大香蕉| 国产亚洲精品第一综合不卡 | 超碰97精品在线观看| 国产无遮挡羞羞视频在线观看| a级毛色黄片| 少妇人妻精品综合一区二区| 精品99又大又爽又粗少妇毛片| 欧美日韩亚洲高清精品| 又黄又爽又刺激的免费视频.| 王馨瑶露胸无遮挡在线观看| 嫩草影院入口| 一区二区三区免费毛片| av不卡在线播放| a级毛片黄视频| 色婷婷久久久亚洲欧美| 国产亚洲精品久久久com| 国内精品宾馆在线| 亚洲四区av| 一本一本综合久久| 一本一本综合久久| 简卡轻食公司| 99热这里只有是精品在线观看| 丰满饥渴人妻一区二区三| 欧美3d第一页| 18禁观看日本| 欧美精品高潮呻吟av久久| 午夜福利视频在线观看免费| 精品一区二区三卡| 午夜视频国产福利| 成人免费观看视频高清| 亚洲欧美一区二区三区黑人 | 美女内射精品一级片tv| 热99国产精品久久久久久7| 欧美精品亚洲一区二区| 丰满饥渴人妻一区二区三| 国产精品蜜桃在线观看| 天天影视国产精品| 赤兔流量卡办理| 在线精品无人区一区二区三| 国产白丝娇喘喷水9色精品| 一级毛片黄色毛片免费观看视频| 男的添女的下面高潮视频| 久久久午夜欧美精品| 69精品国产乱码久久久| 亚洲精品,欧美精品| 观看av在线不卡| 成人国产麻豆网| 一级毛片 在线播放| 亚洲综合色网址| 久久久精品94久久精品| 国产综合精华液| 国产在视频线精品| 精品人妻熟女av久视频| 欧美3d第一页| 国产深夜福利视频在线观看| 又黄又爽又刺激的免费视频.| 精品人妻在线不人妻| 精品少妇黑人巨大在线播放| 黑丝袜美女国产一区| 22中文网久久字幕| 亚洲精品久久久久久婷婷小说| 一区在线观看完整版| 色婷婷av一区二区三区视频| 久久久久久伊人网av| 少妇的逼好多水| 黄片播放在线免费| 久久青草综合色| 大片免费播放器 马上看| 免费日韩欧美在线观看| 中文乱码字字幕精品一区二区三区| 午夜老司机福利剧场| 国产成人aa在线观看| 久久精品久久久久久久性| 狠狠精品人妻久久久久久综合| av黄色大香蕉| 欧美变态另类bdsm刘玥| 久久精品国产亚洲av涩爱| 国产男女超爽视频在线观看| a级毛片在线看网站| 日本色播在线视频| 国产成人免费观看mmmm| 国产亚洲欧美精品永久| 亚洲国产精品一区三区| 精品卡一卡二卡四卡免费| a 毛片基地| 久热久热在线精品观看| 蜜桃久久精品国产亚洲av| 午夜福利视频在线观看免费| 蜜桃在线观看..| 色94色欧美一区二区| 交换朋友夫妻互换小说| 亚洲熟女精品中文字幕| 日韩中文字幕视频在线看片| 美女xxoo啪啪120秒动态图| videossex国产| 国产精品一区www在线观看| 特大巨黑吊av在线直播| 久久国产精品男人的天堂亚洲 | 午夜免费观看性视频| 久久久久国产网址| 制服丝袜香蕉在线| 如日韩欧美国产精品一区二区三区 | 亚洲精品久久久久久婷婷小说| 人人妻人人澡人人看| 一级片'在线观看视频| 国产亚洲最大av| 国语对白做爰xxxⅹ性视频网站| 国产日韩欧美视频二区| 久久精品国产a三级三级三级| 亚洲精品久久成人aⅴ小说 | 久久久久久久久久久免费av| 满18在线观看网站| 亚洲精品日韩av片在线观看| 99久久中文字幕三级久久日本| 亚洲精品久久午夜乱码| 久久人人爽人人片av| 午夜影院在线不卡| 亚洲精品美女久久av网站| 欧美精品一区二区大全| 99热国产这里只有精品6| 日韩av免费高清视频| 亚洲国产最新在线播放| 啦啦啦啦在线视频资源| 亚洲精品久久成人aⅴ小说 | 婷婷色综合大香蕉| 一级a做视频免费观看| 最后的刺客免费高清国语| 国产亚洲一区二区精品| 亚洲第一av免费看| 欧美激情极品国产一区二区三区 | 中文天堂在线官网| 精品国产乱码久久久久久小说| 欧美性感艳星| 久久久久国产精品人妻一区二区| 免费黄网站久久成人精品| 久久毛片免费看一区二区三区| 亚洲人成77777在线视频| 国产精品国产三级专区第一集| 大码成人一级视频| 国产成人免费无遮挡视频| 久久人人爽人人片av| 丁香六月天网| www.av在线官网国产| 婷婷色麻豆天堂久久| 日韩大片免费观看网站| 天美传媒精品一区二区| 老司机影院成人| 老熟女久久久| 日韩在线高清观看一区二区三区| 国产免费福利视频在线观看| 国产国语露脸激情在线看| av天堂久久9| 欧美 亚洲 国产 日韩一| 晚上一个人看的免费电影| 在线观看免费视频网站a站| 色网站视频免费| 欧美日韩精品成人综合77777| 97在线视频观看| 精品久久久精品久久久| 如何舔出高潮| 日本黄大片高清| 成人影院久久| 热re99久久精品国产66热6| 国产在线视频一区二区| 少妇人妻 视频| 插阴视频在线观看视频| 丝袜在线中文字幕| 日本黄色日本黄色录像| 草草在线视频免费看| 国模一区二区三区四区视频| 精品人妻一区二区三区麻豆| 亚洲av二区三区四区| 人妻制服诱惑在线中文字幕| 久久久久网色| 久久久精品区二区三区| 一级爰片在线观看| 成年美女黄网站色视频大全免费 | 日韩欧美精品免费久久| 高清在线视频一区二区三区| 黑丝袜美女国产一区| 久久久久视频综合| 久久国内精品自在自线图片| 秋霞伦理黄片| 亚洲精品aⅴ在线观看| 国产男女内射视频| 亚洲国产色片| 一区二区日韩欧美中文字幕 | 97超碰精品成人国产| freevideosex欧美| 搡老乐熟女国产| 如何舔出高潮| 免费高清在线观看日韩| 精品一区二区免费观看| 日韩欧美精品免费久久| 又黄又爽又刺激的免费视频.| 母亲3免费完整高清在线观看 | 91精品国产九色| 人人妻人人爽人人添夜夜欢视频| 26uuu在线亚洲综合色| 久久国产亚洲av麻豆专区| 99热全是精品| 99国产精品免费福利视频| 男人操女人黄网站| 久久97久久精品| 97在线人人人人妻| 人体艺术视频欧美日本| 日本vs欧美在线观看视频| 大香蕉97超碰在线| 欧美变态另类bdsm刘玥| 街头女战士在线观看网站| 视频中文字幕在线观看| 亚洲精品乱码久久久v下载方式| 精品国产露脸久久av麻豆| 日韩av免费高清视频| 午夜精品国产一区二区电影| 亚洲欧美一区二区三区黑人 | 赤兔流量卡办理| 亚洲天堂av无毛| 王馨瑶露胸无遮挡在线观看| 国产色爽女视频免费观看| 三级国产精品欧美在线观看| av免费在线看不卡| 2021少妇久久久久久久久久久| 亚洲精品久久成人aⅴ小说 | 伦理电影免费视频| 国产无遮挡羞羞视频在线观看| 五月开心婷婷网| 国产成人免费观看mmmm| 性色avwww在线观看| 中文字幕久久专区| 美女大奶头黄色视频| 18禁在线播放成人免费| 亚洲精品乱码久久久v下载方式| 国产成人精品久久久久久| 欧美少妇被猛烈插入视频| 晚上一个人看的免费电影| 久久久国产精品麻豆| 中国三级夫妇交换| 老司机亚洲免费影院| 国产精品人妻久久久久久| 婷婷色综合大香蕉| 亚洲综合精品二区| 中文字幕久久专区| av播播在线观看一区| 91精品国产国语对白视频| 国产一区有黄有色的免费视频| 视频在线观看一区二区三区| 国产综合精华液| 老女人水多毛片| 母亲3免费完整高清在线观看 | 久久精品国产亚洲网站| 亚洲国产精品999| 精品久久久久久电影网| 亚洲av免费高清在线观看| 国产不卡av网站在线观看| 免费观看a级毛片全部| 高清不卡的av网站| 精品一区二区三卡| 日本与韩国留学比较| 中文字幕av电影在线播放| 国产色爽女视频免费观看| 国产精品 国内视频| 久久精品久久精品一区二区三区| 日本91视频免费播放| 日韩av在线免费看完整版不卡| 国产黄片视频在线免费观看| 亚洲国产毛片av蜜桃av| 国产黄频视频在线观看| 插阴视频在线观看视频| 亚洲国产毛片av蜜桃av| 黄片无遮挡物在线观看| 亚洲国产最新在线播放| 少妇的逼水好多| 日韩精品有码人妻一区| 国产在视频线精品| 免费人成在线观看视频色| av播播在线观看一区| 2021少妇久久久久久久久久久| 91精品三级在线观看| 久久久久久伊人网av| 丝袜脚勾引网站| 亚洲久久久国产精品| 亚洲经典国产精华液单| 精品一区二区免费观看| 熟女人妻精品中文字幕| 秋霞伦理黄片| 国产国语露脸激情在线看| av不卡在线播放| 尾随美女入室| 国产欧美日韩综合在线一区二区| 99久久精品国产国产毛片| 亚洲无线观看免费| 国产欧美另类精品又又久久亚洲欧美| 国产精品久久久久久精品古装| 人人妻人人澡人人爽人人夜夜| 国产精品一区二区在线不卡| 老熟女久久久| 2022亚洲国产成人精品| 亚洲婷婷狠狠爱综合网| 国产成人午夜福利电影在线观看| 亚洲怡红院男人天堂| 狂野欧美白嫩少妇大欣赏| 国产成人精品在线电影| 日本猛色少妇xxxxx猛交久久| 如何舔出高潮| 久久热精品热| av电影中文网址| 国产av一区二区精品久久| 视频区图区小说| 久久久久精品性色| 成人亚洲精品一区在线观看| 多毛熟女@视频| 亚洲av在线观看美女高潮| 国产黄频视频在线观看| 能在线免费看毛片的网站| 免费看不卡的av| 桃花免费在线播放| 另类精品久久| 狠狠婷婷综合久久久久久88av| 欧美性感艳星| 国产淫语在线视频| 亚洲欧美成人综合另类久久久| 亚洲国产精品成人久久小说| 久久鲁丝午夜福利片| 国产男女内射视频| 狠狠精品人妻久久久久久综合| 日韩一本色道免费dvd| 大片免费播放器 马上看| 99久久精品一区二区三区| 97超视频在线观看视频| 精品99又大又爽又粗少妇毛片| 两个人的视频大全免费| av卡一久久| 亚洲,欧美,日韩| 久久人人爽人人爽人人片va| 丝袜在线中文字幕| 边亲边吃奶的免费视频| 午夜91福利影院| 高清视频免费观看一区二区| 国国产精品蜜臀av免费| 午夜影院在线不卡| 亚洲欧洲国产日韩| 国产精品国产三级国产专区5o| 看十八女毛片水多多多| 欧美性感艳星| 蜜桃在线观看..| 国产成人精品婷婷| 99精国产麻豆久久婷婷| 亚洲精品日韩在线中文字幕| 2021少妇久久久久久久久久久| 国产白丝娇喘喷水9色精品| 欧美少妇被猛烈插入视频| 新久久久久国产一级毛片| 欧美日韩成人在线一区二区| 日本欧美国产在线视频| www.av在线官网国产| 国产精品欧美亚洲77777| 亚洲国产毛片av蜜桃av| 最黄视频免费看| 中文字幕制服av| 日韩精品有码人妻一区| 久久精品国产a三级三级三级| 亚洲国产精品专区欧美| 丁香六月天网| 亚洲欧美成人综合另类久久久| 51国产日韩欧美| 全区人妻精品视频| 国产成人av激情在线播放 | 国产一区二区三区综合在线观看 | 国产精品免费大片| 日日啪夜夜爽| 婷婷色综合大香蕉| 国产乱来视频区| 久久人妻熟女aⅴ| 黄色毛片三级朝国网站| 国产又色又爽无遮挡免| av免费在线看不卡| 中文字幕亚洲精品专区| 国产精品成人在线| 狂野欧美白嫩少妇大欣赏| 国产精品一区二区三区四区免费观看| 黑人猛操日本美女一级片| 久久精品国产鲁丝片午夜精品| 观看av在线不卡| 久久99热6这里只有精品| 制服丝袜香蕉在线| 久久韩国三级中文字幕| 麻豆成人av视频| 男女国产视频网站| 国产日韩欧美在线精品| 日韩不卡一区二区三区视频在线| 国产亚洲一区二区精品| 男的添女的下面高潮视频| 插阴视频在线观看视频| 色94色欧美一区二区| 亚洲,一卡二卡三卡| 国产精品国产av在线观看| 啦啦啦在线观看免费高清www| 丝瓜视频免费看黄片| 国产一区二区三区av在线| a级毛色黄片| 男女啪啪激烈高潮av片| 在线 av 中文字幕| 免费观看在线日韩| 一本久久精品| 亚洲精品美女久久av网站| 久久ye,这里只有精品| 精品一品国产午夜福利视频| 我的女老师完整版在线观看| 乱人伦中国视频| 日韩不卡一区二区三区视频在线| 国产亚洲一区二区精品| 男女边摸边吃奶| 午夜激情av网站| 寂寞人妻少妇视频99o| 美女cb高潮喷水在线观看| 免费av不卡在线播放| 国产成人a∨麻豆精品| 国产成人精品福利久久| 久久国产亚洲av麻豆专区| 看免费成人av毛片| 91在线精品国自产拍蜜月| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | 啦啦啦中文免费视频观看日本| 熟女av电影| 欧美少妇被猛烈插入视频| 午夜免费观看性视频| 有码 亚洲区| 欧美人与善性xxx| 九色成人免费人妻av| 又黄又爽又刺激的免费视频.| 一本—道久久a久久精品蜜桃钙片| av在线app专区| 久久久久精品久久久久真实原创| 免费看av在线观看网站| 国产精品偷伦视频观看了| 一边摸一边做爽爽视频免费| 22中文网久久字幕| 国产成人aa在线观看| 22中文网久久字幕| 丝袜脚勾引网站| 久久免费观看电影| 日韩一区二区三区影片| 最新的欧美精品一区二区| 51国产日韩欧美| 最新的欧美精品一区二区| 国产精品一二三区在线看| 亚洲精品日韩在线中文字幕| 人妻系列 视频| 国产极品天堂在线| 狂野欧美激情性bbbbbb| 人成视频在线观看免费观看| 在线 av 中文字幕| 国产永久视频网站| 午夜激情福利司机影院| freevideosex欧美| 99热6这里只有精品| 免费高清在线观看视频在线观看| 一边摸一边做爽爽视频免费| 蜜臀久久99精品久久宅男| 久久99热6这里只有精品| 亚洲熟女精品中文字幕| 晚上一个人看的免费电影| 69精品国产乱码久久久| 亚洲人成网站在线播| 亚洲美女视频黄频| 伊人久久国产一区二区| 五月玫瑰六月丁香| 国产淫语在线视频| 欧美少妇被猛烈插入视频| 另类亚洲欧美激情| 久久午夜综合久久蜜桃| 99re6热这里在线精品视频| 黄色怎么调成土黄色| 精品酒店卫生间| av国产久精品久网站免费入址| 99精国产麻豆久久婷婷| 夜夜骑夜夜射夜夜干| 欧美三级亚洲精品| 日韩亚洲欧美综合| 日韩一本色道免费dvd| 丰满饥渴人妻一区二区三| 日本色播在线视频| 久久久久久久国产电影| 久久久久久伊人网av| 欧美精品高潮呻吟av久久| 亚洲内射少妇av| 免费看不卡的av| 中国国产av一级| 少妇高潮的动态图| 性色avwww在线观看| 日本黄大片高清| a 毛片基地| 中文欧美无线码| 纯流量卡能插随身wifi吗| 菩萨蛮人人尽说江南好唐韦庄|