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

    Hybrid simulation of q=1 high-order harmonics driven by passing energetic particles in tokamak plasmas

    2023-12-18 03:54:34ShengLIU劉勝ZhenzhenREN任珍珍WeihuaWANG汪衛(wèi)華WeiSHEN申偉JinhongYANG楊錦宏andHongweiNING寧洪偉
    Plasma Science and Technology 2023年12期
    關鍵詞:劉勝

    Sheng LIU (劉勝) ,Zhenzhen REN (任珍珍) ,Weihua WANG (汪衛(wèi)華),* ,Wei SHEN(申偉) ,Jinhong YANG(楊錦宏) and Hongwei NING(寧洪偉)

    1 Institutes of Physical Science and Information Technology,Anhui University,Hefei 230601,People’s Republic of China

    2 School of Physics and Optoelectronic Engineering,Anhui University,Hefei 230601,People’s Republic of China

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

    Abstract High-order harmonics q(ψs)=1 energetic particle modes(EPMs)have been observed in toroidal plasmas experiments with neutral beam injection.To investigate these phenomena,linear properties and nonlinear dynamics of these EPMs driven by passing energetic particles(EPs)are studied via the global hybrid kinetic-magnetohydrodynamic code M3D-K.Simulation results demonstrate that passing EPs’effects on high mode-number harmonics(q(ψs)=m/n=2/2,3/3,4/4) instability are more obvious than the q(ψs)=1/1 mode,especially when q-profile is sufficiently flat in the core region.Furthermore,the effects of the pitch angle Λ0 and beam ion pressure Phot/Ptotal on the features of high n components are also analyzed specifically.It is found that there exists only one resonant condition for these EPMs.In the nonlinear phase,these high mode-number harmonics can induce significant energetic ions redistribution and chirping up phenomena,which differs from the classical fishbone excited by passing EPs.These discoveries are conducive to better apprehend the underlying physical mechanisms of the highorder harmonics driven by passing EPs.

    Keywords: high-order harmonics,passing energetic particles,wave-particle resonance,tokamak

    1.Introduction

    In magnetically-confined fusion devices,energetic particles(EPs) physics is a critical issue in achieving high confinement performance and steady-state operation.To increase the temperature of plasma,electron cyclotron resonance heating,ion cyclotron resonance heating and neutral beam injection(NBI) can be applied in tokamaks.In the process,a large number of EPs can be produced,which could interact with magnetohydrodynamics (MHD) activities or drive various instabilities including fishbone modes and Alfvén eigenmodes[1-6].Those instabilities can contribute to EPs loss/redistribution.In the Poloidal Divertor Experiment (PDX),fishbone instability excited by trapped particles was first observed under the perpendicular NBI[7].Subsequently,in the Princeton beta experiment (PBX) with tangential NBI,the fishbone-like internal kink modes excited by passing particles were also reported [8,9].Since then,in the last few decades,the experimental phenomenon of energetic ions acting on the q(ψs)=1/1 mode has been widely found in the tokamak experiments,for instance,Tokamak Fusion Test Reactor(TFTR) [10],Joint European Torus (JET) [11],HL-2A[12-14],and DIII-D [15].Furthermore,in comparison to q(ψs)=1/1 EPM,the instability of the high mode-number harmonics is also unstable which has been observed in MAST[16],EAST [17],ASDEX [18] and HT-7 [19] during NBI heating.Therefore,the study of high-order energetic particle modes (EPMs) excited by EPs in tokamaks is of great significance.

    Theoretically,we can better apprehend the physical mechanisms of these MHD modes[20-24],including internal kinks,sawteeth,long livedmodes,and fishbones.In particular,the past related research reports have mainly focused on m=1,n=1 mode,which is resonantly excited by energetic ions generated from perpendicular or tangential NBI[25-29].When trapped particles excite the q(ψs)=1/1 mode,and the wave-particle resonance can be specifically expressed as ω=ωφ.When passing particles drive the q(ψs)=1/1 mode,and the wave-particle resonance can be specifically expressed as ω=ωφor ω=ωφ-ωθ[23,29-33],where ωφand ωθare the toroidal transit frequency and the poloidal transit frequency,respectively.m=n >1 EPMs are studied in recent years based on the experimental observations of these modes [34,35].Zhang et al numerically investigated q(ψs)=1/1,2/2 EPMs driven by trapped particles in tokamak[36].Up to now,the simulation results of the trapped particles affecting the high-order harmonics EPMs have been well reported[37].Nevertheless,it should be stressed that the impact of the passing EPs on the q(ψs)=1(n ≥1)harmonics were similarly observed in some experiments [38,39],and few simulation studies have been taken into account so far.The physical process of these n >1 EPMs can reduce the core plasma confinement,causing significant energetic-ion redistribution.To effectively achieve steady-state operation in tokamak plasmas,investigating the physical mechanism of high order harmonics driven by passing EPs is essential.

    In this work,we concentrate our attention on the characteristics of the high-order harmonics EPMs,and numerically analyze the passing EPs impact on the mode excitation via M3D-K code [40].The rest of this work is organized as follows.In section 2,the simulation code and associated parameters are briefly introduced.In sections 3 and 4,the linear and nonlinear numerical results with passing EPs are presented,respectively.Ultimately,the summary and discussion are in section 5.

    2.Numerical model and basic parameters

    2.1.Physical analytical model

    M3D-K is a global hybrid kinetic-magnetohydrodynamic(MHD) code,we mainly apply this code to numerically calculate the extended MHD equations and the drift-kinetic equations,respectively,which can better simulate the physical phenomenon of q(ψs)=1 EPMs in toroidal plasmas[25].The plasmas in this code are composed of the EPs and the background plasma.The EPs are demonstrated by driftkinetic equations,whose solution is δ? particle-in-cell method.The background plasma contains electrons and ions,which is considered as a solitary fluid in the MHD equations.The finite element method is used to calculate these equations.Up to now,M3D-K code was widely used to simulate MHD instabilities excited by EPs,such as fishbone mode,RSAE,EPM and TAE [41-47].

    Figure1.(a) The pressure p equilibrium profiles and (b) q profile.

    2.2.Equilibrium profiles and parameters setup

    In this simulation,the main parameters are listed as follows based on HL-2A like conditions:circular cross-section,?-1=R0/a=3.6667,elongation κ=1,triangularity δ=0,B0=1.37 T,βhot/βtotal=0.5,the central total beta is fixed at βtotal=1.2%including the beta of both bulk plasma and EPs,E0=37.78 keV,where vA=B0/(μ0ρ0)1/2is Alfvén speed,τA=R0/vAis Alfvén time and ωA=vA/R0is Alfvén frequency.

    Figure 1(a) shows the radial profile of total pressure,which is as follows:

    where P0is pressure at the magnetic axis,the magnetic poloidal flux Ψ is a radial variable,the edge of the plasma Ψ=1,the center of the plasma Ψ=0.Correspondingly,for the q(Ψs)=1 EPMs,figure 1(b) is the spatial profiles of the safety factor with q0=0.9814.

    The beam ion is a slowing down distribution in velocity space and Gaussian distribution in pitch angle space(Λ=μB0/E).The EPs distribution is as follows:

    where ΔΨ=0.4,ΔΛ=0.3,Λ0=0.7,c is a normalization factor,H is the step function,vc=0.962vAis the critical velocity,which is as follows:

    Figure2.(a)Growth rate from toroidal mode number n=1 to n=4 and (b) corresponding mode frequency.

    3.Linear simulation results

    3.1.Features of the high-order harmonics

    In the linear simulation part of this work,Λ0=0.7 is chosen for analyzing the high-order harmonics,and the width of the pitch angle is set to be ΔΛ=0.3.Due to the impact of the pitch angle width,the distribution of pitch angle has a certain variation range,and the energy particles include passing particles and trapped particles in the beam ion distribution.In addition,according to wave-particle resonant interaction,we found that the high-order harmonics are driven by passing particles with Λ0=0.7,as shown in figure 4.These are the reasons why passing particles are analyzed with near perpendicular neutral beam injection in our work.

    In the following the numerical results of the EPMs driven by passing EPs are presented.Initially,the numerical simulations without NBI are carried out,that is,MHD simulations are firstly performed,and these results show that MHD mode presents a steady state.Secondly,when NBI heating is included,according to the analysis of perturbed distribution δE,we find that the high order EPMs will become more unstable under the effects of passing EPs.As is shown in figure 2,passing EPs’effects on the q(ψs)=2/2,3/3,4/4 harmonics instability are more significant than those on the q(ψs)=1/1 mode,indicating that these high n component instabilities become more dominant.Furthermore,the mode frequency becomes higher when the toroidal mode number n increases.

    It is shown in figure 3 that linear mode numbers of q(ψs)=1 EPMs are n=1,2,3,4,separately.U is the velocity stream function,which is associated with the plasma velocityν=R2??⊥U× ??+?χ+ν???,where the variables such as toroidal angle ? and compressible component χ are included.Due to the energetic ions’effect,the mode structures for these high order harmonics show slightly twisted feature,which is different from the typical internal kink mode.The mode structures are inside the q(ψs)=1 rational surfaces,which are denoted by the black circles.

    Due to wave-particle resonant interaction,generally,the free energy of the radial variation of the EPs distribution excites the EPMs instability.For passing particles,the resonant condition expression is as follows [48]:

    where p is an integer.Previous analyzes have initially focused on the q(ψs)=1/1 mode,satisfying the primary resonances with p=-1 for the EPMs branch of low-frequency or p=0 for the EPMs branch of high-frequency[31].In this work,we find that high frequency branch resonances can also be applied to high n harmonics.Figure 4 shows the perturbed fast ion energy δE in the phase space P?-E.The pitch angle parameter is fixed when we plot resonant conditions of different harmonics,where Λ is defined as Λ=μB0/E.When Λ is fixed,the resonance curves of the different harmonics are fitted and plotted in figure 4.To track the resonant locations,the significant changes of particle energy in the linear phase are investigated in figure 4,which shows the phase space locations of the particles with∣δE∣>0.125Emax,whereδEis the particle energy change,Emaxis the maximum particle energy change.As a result,there are many particles with significant energy changes located along the resonant curves in figure 4.The magnetic moment values of these resonant particles vary in figure 4 with different energies.In addition,figure 4 shows that trapped particles are mainly centered in the regions of larger P?,while passing EPs are mainly centered in the regions of smaller P?.For the passing and trapped particles,we use the horizontal dashed lines as the dividing line of the different types of particles in figure 4.Significantly,the different harmonics with q(ψs)=1/1,2/2,3/3,4/4 modes have only one branch of resonance condition,with the resonance corresponding to p=0,-1,-2,-3 respectively.Furthermore,we observe that large δE for different harmonic modes are all located below the dashed line.Therefore,these simulations show that the resonance of waves and particles is caused by the passing EPs when Λ0=0.7.Because toroidal angular momentum can be expressed asP?=eψ+mDν‖RB?B,which is related to ψ.Compared to the q(ψs)=1/1 mode,n=2,3,4 components have a slightly smaller toroidal angular momentum P?,which indicates that locations of passing particle resonances are nearer to the core region of the tokamak plasmas.Generally,for energetic particle driven instabilities,the correlative growth rate expression isγ∝+[36],where P?is canonical toroidal angular momentum.This formula shows that the drive related to the special gradient of distribution is proportional to toroidal mode number n,which partially explains why n=2,3,4 components instabilities are more unstable compared to the n=1 component.Furthermore,it is shown in the following that the excitation of high order harmonics depends on the shape of q profile and various related parameters.

    Figure3.The velocity stream function U.(a)The(1,1)harmonic,(b)the(2,2)harmonic,(c)the(3,3)harmonic and(d)the(4,4)harmonic.

    Figure4.The resonant conditions of different harmonics: including (a) the (1,1) harmonic,(b) the(2,2)harmonic,(c) the (3,3) harmonic and(d)the(4,4)harmonic.The horizontal dashed lines in figures 4(a)-(d)show the dividing line for the different types of particles.Above the dotted line is the trapped particles region,and below the line is the passing particles region.The color bar shows the amplitude of the energy change δE of energetic particles.

    3.2.Effects of important parameters on high-order harmonics

    Figure5.(a) Growth rate of the q(ψs)=1 EPMs from pitch angles Λ0=0.3 to Λ0=1.0 and (b) corresponding mode frequency.

    Figure6.Resonant condition: (a) Λ0=0.8,the mode frequency is 0.090128ωA,(b) Λ0=0.9,the mode frequency is 0.02387ωA and(c) Λ0=1.0,the mode frequency is 0.02285ωA for the q(ψs)=1/1 mode.

    Figure7.(a) Growth rate of the q(ψs)=1 harmonics from the function of Phot/Ptotal=0.2 to 0.8 and (b) corresponding mode frequency.

    Figure8.Resonant condition: (a) Phot/Ptotal=0.3,the mode frequency of the n=1 mode is 0.1092926045ωA and (b) Phot/Ptotal=0.6,the mode frequency of the n=1 mode is 0.0840079ωA for Λ0=0.7.

    To testify to the impact of the central pitch angle (Λ0),it is shown in figure 5 that the pitch angle varies from 0.3 to 1.0,which affects the mode frequency and instability of different harmonics.The q(ψs)=1(n ≥1)harmonics are stable in Λ0≤0.3.However,with Λ0increasing,the energy particles effect on the stability of these EPMs becomes more obvious.Furthermore,the linear frequency of these EPMs decreases as Λ0increases.The reason is that the transit frequency of passing EPs is determined by particle parallel velocity,associated with pitch angle Λ[49].According to the changes of the perturbed distribution δE,we mark the most prominent locations of the wave-particle resonance under different Λ0.For the n=1 component,Λ0=0.8,0.9,1.0 are chosen to analyze as the mode frequency changes significantly in this region.When Λ0=0.8,the modes are mainly excited by passing particles,while Λ0=0.9,1.0,the modes are mainly excited by trapped particles.The resonant results are shown in figures 6(a)-(c),when Λ0=0.8,the m=1,n=1 mode has only one resonant condition with p=0,and passing particles energy changes significantly below the red dashed line.However,when Λ0=0.9,the m=1,n=1 mode has two resonant conditions with p=0,1.When Λ0=1.0,the m=1,n=1 mode has also one branch resonant conditions with p=0.For both the EPMs with Λ0=0.9 and Λ0=1.0,particle energy changes are all located above the red dashed line,indicating that the impact of the trapped particles on the n=1 EPMs is dominant when Λ0≥0.9.

    Figure9.(a)Spatial profiles of safety factor varying from q0=0.82 to q0=0.9814,(b)the corresponding linear growth rate for different toroidal mode numbers and (c) corresponding mode frequency.

    As is illustrated in figure 7,with fixed total pressure,the effects of the EPs pressure fraction Phot/Ptotalon the excitation of the q(ψs)=1 (n ≥1) harmonics are examined,which become unstable when the EPs pressure over a certain threshold.As the beam ion pressure fraction Phot/Ptotalincreases,the growth rate of these EPMs gradually increases.However,the mode frequency of these EPMs decreases with the beam ion pressure increasing,which is mainly associated with the resonance locations.Taking q(ψs)=1/1 mode for instance,figures 8(a) and (b) respectively show the resonant conditions of different beam ion pressure fractions Phot/Ptotal,when Phot/Ptotalis smaller,the resonance locations of the passing particles are closer to the core.

    In the previous study,it is found that the instability of the q(ψs)=1 modes is associated with δq=|q-1| [50].To testify the dependence of the q(ψs)=1(n ≥1)harmonics on the q-profile,we change q0from 0.82 to 0.9814 in figure 9(a).Figure 9(b)shows that the stability of these EPMs is sensitive to the safety factor and the flattened region of q-profile.For a smaller q0(q0=0.82),the n=1 component grows linearly at the rate of γτA≈0.005,but the high-order EPMs are stable.With q0gradually increasing and q-profile becoming more flattened,the instability of the high mode-number harmonics becomes very unstable,which is even more dominant when q0>0.9.Furthermore,it is found that the mode frequency changes of the q(ψs)=2/2,3/3,4/4 components are different from the q(ψs)=1/1 mode in figure 9(c),whose mode frequencies with different q0are almost the same.

    Figure10.(a)Growth rate of the q(ψs)=1 harmonics with different ΔΛ=0.2 and 0.3 and (b) corresponding mode frequency.

    Figure11.Time evolution of the kinetic energy for different toroidal mode numbers.

    In some experiments in HL-2A and EAST with flat q profile in the core region,the high order harmonics do not emerge.There are two possible explanations for this puzzle.Firstly,according to the theoretical work by Hastie and Hender [50],the stability criterion for higher m modes without energetic particle effects is given bywhereδq=∣q-1 ∣,r1is the radial location of q=1 surface.As a result,the stability of higher m modes depends on various parameters including the beta value,the location of q=1 surface,and the mode number.Secondly,in the following,we show that a few parameters related to EPs can also affect the linear growth rate of high order harmonics,including Λ,Phot/Ptotal,etc.As a result,although the condition of flat q profile is satisfied,the excitation of high order harmonics still depends on various related parameters,which explains why in some experiments the high order harmonics do not dominate with flat q profile.

    Figure12.Fourier spectrogram including (a) the (1,1) mode,(b) the (2,2) mode,(c) the (3,3) mode,(d) the (4,4) mode.

    Figure13.Time evolution of 1D distribution function with (a) n=1,(b) n=2,(c) n=3,(d) n=4.

    To better explore the impact of the pitch angle width on the high-order harmonics,Λ0=0.7 is chosen for analysis.The linear growth rates and mode frequencies are shown in figure 10,for two different ΔΛ: ΔΛ=0.2 (blue dashed line)and ΔΛ=0.3 (red dashed line),and the radial width is Δψ=0.4(ψmax-ψmin).The instability of different modes decreases with increasing the width of the pitch angle (ΔΛ),and the q(ψs)=1 high mode-number harmonics are still dominant compared to the n=1 component.However,the frequencies of different modes are nearly independent of ΔΛ.

    Figure14.Time evolution of 2D distribution function with (a) n=1,(b) n=2,(c) n=3,and (d) n=4.

    4.Nonlinear simulation results

    4.1.Saturation level analysis of different modes

    In the nonlinear simulation part of this work,the kinetic energy evolution of different toroidal modes without MHD non-linearity is shown in figure 11.Here we choose Λ0=0.7,which is consistent with the linear simulations.Significantly,these high n components driven by passing EPs have larger saturation level in the nonlinear phase.Generally,the EPMs saturation is associated with the flattening of the particle distribution [25].In [37],it was proposed that the saturation levels of the q(ψs)=2/2,3/3,4/4 harmonics excited by trapped EPs are smaller than that of the q(ψs)=1/1 harmonic in the nonlinear phase.Nevertheless,the high mode-number harmonics driven by passing EPs are always dominant both in the linear and nonlinear saturation phases.These results are different from the previous work [37],in which the high mode-number harmonics driven by trapped EPs are just dominant in the linear phase.

    4.2.Characteristics of the mode frequency

    The phenomenon of frequency chirping in the nonlinear phase is caused by the wave-particle resonance.Figure 12 shows the frequency evolution of the high-order EPMs.For the q(ψs)=1/1 mode,the frequency shifts downward about δω/ω=13%.Nevertheless,the frequency of the high-order EPMs shifts upward.For the q(ψs)=2/2 mode,the frequency starts to stay around a constant value with no significant change,and then shifts upward about δω/ω=20% at around t ≈2500τA,and nearly no variation afterwards.Furthermore,the frequency evolution of the n=3 and n=4 parts is similar,where the chirping range is even larger,and shifts upward about δω/ω=41%at around t ≈1500τAand δω/ω=40%at around t ≈1200τArespectively.These physical phenomena of frequency chirping are due to the radial flattening of the EPs distribution,which may induce the resonance island shift in the phase space and cause the loss of energetic ions.

    4.3.Evolution of energetic particles’distribution function

    Due to the high-order harmonics instability,EPs generate redistribution in phase space.Figures 13(a)-(d) show the 1D distribution function ?(P?) of passing EPs.Firstly,in comparison with the nonlinear later phase,the redistribution levels are relatively small in the early initial saturation phase.Then,with the nonlinear time evolution of these EPMs,the redistribution of the beam ions expands outwards/inwards radially and has a large flattening region,located nearer the core region in the radial direction.Moreover,figures 13(b)-(d)show that the redistribution regions of the q(ψs)=2/2,3/3,4/4 harmonics are relatively larger than that of the q(ψs)=1/1 mode.There are two possible explanations for these results:(1)the ratios of the linear growth rate to the mode frequency are relative larger for the high n components.For instance,the n=1 component is γ/ω=0.093 while the n=3 component is γ/ω=0.175;(2)the saturation levels of the high mode-number harmonics are much higher,especially for the saturation peak of kinetic energy when n=3,which reaches the highest level.Consequently,the high mode-number harmonics can induce a more obvious redistribution of energetic ions.

    Figures 14(a)-(d) show the 2D distribution function ?(P?) of passing EPs at E ≈0.3,0.4,0.46,0.5,which is marked by the horizontal black dashed.In addition,we plot the vertical black dashed at P?=0.4 in order to better observe redistribution.The movements of these particles in the phase space P?-E are determined by the formula[28].Here,the redistribution of EPs includes the first saturated phase and the later nonlinear phase.It is found that EPs have migrated from the core to the edge in the tokamak experiment.Especially with the nonlinear evolution of these EPMs,the levels of redistribution for the n=2,3,4 modes are further enhanced.Thus,the significant redistribution of energetic ions for the q(ψs)=1/1 mode is relatively smaller.As a result,energetic ions are transported in the phase space P?-E,indicating that the high-order EPMs can induce the redistribution of energy ions.

    5.Summary and discussion

    In summary,the linear and nonlinear evolutions of the q(ψs)=1 high mode-number harmonics driven by passing EPs have been systematically studied via the global kinetic-MHD code M3D-K.In the linear phase,the impacts of the important parameters are first analyzed in this work.With the central pitch angle or the energetic particle pressure increasing,the instability of these EPMs becomes stronger.Specifically,the transition of the wave-particle resonance condition from passing particles dominant to trapped particles dominant is observed when the central pitch angle exceeds a certain range.Furthermore,numerical results show that passing EPs’ effects on high mode-number harmonics instability are more significant than that on the q(ψs)=1/1 mode.The high-order EPMs driven by passing EPs satisfy different resonant conditions under different pitch angles.

    Additionally,the nonlinear features of these EPMs have also been investigated.It is found that frequencies of the n >1 components chirp up,which is different from the classical fishbone.Ultimately,because high n components have a relatively larger saturated level,we find that strong redistribution of the passing EPs can be induced by the high-order harmonics.The instability phenomena of these EPMs driven by passing EPs can be commonly observed in HL-2A experiments.Therefore,the related discoveries in this work are conducive to guiding future tokamak experiments,especially in controlling the instability of high-order harmonics modes.

    Acknowledgments

    This work is supported by National Key R&D Program of China (Nos.2019YFE03050002,2018YFE0310400,and 2022YFE03040002),and National Natural Science Foundation of China (Nos.12005003 and 11975270),and Science Foundation of Institute of Plasma Physics,Chinese Academy of Sciences (No.DSJJ-2022-04).

    ORCID iDs

    猜你喜歡
    劉勝
    河北滿城漢墓出土的飲酒器具
    東方收藏(2024年3期)2024-05-19 10:17:33
    中山靖王,玩轉“障眼法”
    評論:用意象和象征建造當代漢詩的北方迷樓
    詩選刊(2023年1期)2023-05-30 06:19:22
    Anderson localization of a spin–orbit coupled Bose–Einstein condensate in disorder potential
    一袋棉花
    劉勝的詩
    詩選刊(2019年8期)2019-08-12 02:29:36
    職場沒有敵人
    投資
    雜文選刊(2017年9期)2017-09-06 00:44:14
    記劉勝金縷玉衣的第三次修復
    文物春秋(2015年6期)2015-05-30 10:48:04
    伴 娘
    最近最新免费中文字幕在线| 亚洲精品乱码久久久v下载方式| 美女被艹到高潮喷水动态| 一个人看的www免费观看视频| 午夜福利在线在线| 久久久久久大精品| 一级黄片播放器| 成人美女网站在线观看视频| 国产毛片a区久久久久| 欧美日本亚洲视频在线播放| 极品教师在线免费播放| 国产亚洲欧美98| 国产激情偷乱视频一区二区| 999久久久精品免费观看国产| 国产不卡一卡二| 亚洲av一区综合| 免费在线观看亚洲国产| 成人精品一区二区免费| 老熟妇仑乱视频hdxx| 亚州av有码| 中国美女看黄片| 成人亚洲精品av一区二区| 久久久久免费精品人妻一区二区| 日韩精品青青久久久久久| АⅤ资源中文在线天堂| 可以在线观看的亚洲视频| 真实男女啪啪啪动态图| 国产一区二区三区视频了| 国产精品一及| 成年人黄色毛片网站| 又爽又黄a免费视频| 免费无遮挡裸体视频| 国产男靠女视频免费网站| 国产三级黄色录像| 老司机午夜十八禁免费视频| netflix在线观看网站| ponron亚洲| 亚洲18禁久久av| 午夜免费激情av| 国产精品永久免费网站| 久久国产精品影院| 亚洲精品色激情综合| 性色avwww在线观看| 日本在线视频免费播放| 亚洲经典国产精华液单 | 久久久色成人| 久久人妻av系列| 久久国产精品影院| 国产av不卡久久| 日韩亚洲欧美综合| 18+在线观看网站| 大型黄色视频在线免费观看| 国产精品爽爽va在线观看网站| 久99久视频精品免费| 啦啦啦韩国在线观看视频| 国内精品美女久久久久久| 波多野结衣高清作品| 亚洲最大成人手机在线| 真人做人爱边吃奶动态| 亚洲精品一区av在线观看| 久久久久久久久久黄片| 别揉我奶头~嗯~啊~动态视频| 欧美一级a爱片免费观看看| 日本黄大片高清| 亚洲18禁久久av| 亚洲国产色片| 午夜福利在线观看免费完整高清在 | 久久久久久久精品吃奶| 日本在线视频免费播放| 每晚都被弄得嗷嗷叫到高潮| 中文字幕熟女人妻在线| 蜜桃亚洲精品一区二区三区| 黄色配什么色好看| eeuss影院久久| av视频在线观看入口| 亚洲欧美精品综合久久99| 亚洲最大成人中文| 午夜a级毛片| 99久久无色码亚洲精品果冻| 一本精品99久久精品77| 欧美+亚洲+日韩+国产| 欧美性猛交黑人性爽| 欧美xxxx性猛交bbbb| 午夜视频国产福利| 91av网一区二区| 日日干狠狠操夜夜爽| 精品久久久久久成人av| 999久久久精品免费观看国产| 国产国拍精品亚洲av在线观看| 国产真实乱freesex| 可以在线观看毛片的网站| 一区二区三区高清视频在线| 美女xxoo啪啪120秒动态图 | 色哟哟·www| 国产美女午夜福利| 国产精品不卡视频一区二区 | 级片在线观看| 亚洲精品亚洲一区二区| 亚洲七黄色美女视频| 亚洲不卡免费看| 一级黄色大片毛片| 亚洲国产欧美人成| 精品99又大又爽又粗少妇毛片 | 精品久久久久久久久久久久久| 男人和女人高潮做爰伦理| 午夜精品久久久久久毛片777| 久久九九热精品免费| 一个人免费在线观看电影| 男人的好看免费观看在线视频| 国产欧美日韩一区二区精品| 午夜免费男女啪啪视频观看 | 全区人妻精品视频| a级毛片免费高清观看在线播放| 免费在线观看亚洲国产| 国产精品乱码一区二三区的特点| 欧美黑人欧美精品刺激| 亚洲成a人片在线一区二区| 特级一级黄色大片| 少妇被粗大猛烈的视频| 午夜福利视频1000在线观看| 免费人成在线观看视频色| 亚洲 国产 在线| 国产精品自产拍在线观看55亚洲| 男人狂女人下面高潮的视频| 夜夜看夜夜爽夜夜摸| 国产精品三级大全| 校园春色视频在线观看| 午夜福利在线观看免费完整高清在 | 别揉我奶头 嗯啊视频| 婷婷亚洲欧美| 色综合站精品国产| 精品久久久久久成人av| 久久热精品热| 一个人观看的视频www高清免费观看| 黄色丝袜av网址大全| 麻豆久久精品国产亚洲av| 欧美一区二区国产精品久久精品| 久久久国产成人精品二区| 老熟妇乱子伦视频在线观看| 天堂√8在线中文| 国产中年淑女户外野战色| 亚洲,欧美,日韩| 久久精品国产自在天天线| 97热精品久久久久久| 免费在线观看日本一区| 3wmmmm亚洲av在线观看| 国产真实乱freesex| 好看av亚洲va欧美ⅴa在| 久久精品国产自在天天线| 美女cb高潮喷水在线观看| 日日夜夜操网爽| 亚洲成人精品中文字幕电影| 色播亚洲综合网| 久久人妻av系列| 97碰自拍视频| 久久久久国内视频| 国产精品电影一区二区三区| 一区福利在线观看| 亚洲成人久久爱视频| 亚洲欧美日韩东京热| 亚洲激情在线av| 别揉我奶头~嗯~啊~动态视频| 国产成+人综合+亚洲专区| 丰满的人妻完整版| 欧美色欧美亚洲另类二区| 少妇高潮的动态图| 免费高清视频大片| 欧美日韩亚洲国产一区二区在线观看| 在线观看av片永久免费下载| 天天躁日日操中文字幕| 亚洲成av人片在线播放无| 欧美bdsm另类| 精品欧美国产一区二区三| 好看av亚洲va欧美ⅴa在| 久久精品国产亚洲av香蕉五月| 亚洲国产色片| 婷婷亚洲欧美| 国产一区二区三区在线臀色熟女| 中文字幕人妻熟人妻熟丝袜美| 一进一出抽搐gif免费好疼| 国产精品女同一区二区软件 | 国产亚洲精品久久久久久毛片| 久久久久久久午夜电影| 精品久久久久久久久久久久久| 91字幕亚洲| 欧美成人性av电影在线观看| 美女被艹到高潮喷水动态| 五月伊人婷婷丁香| 亚洲人成伊人成综合网2020| 午夜两性在线视频| 69人妻影院| 又黄又爽又刺激的免费视频.| 久久久久亚洲av毛片大全| 一个人观看的视频www高清免费观看| 国产探花在线观看一区二区| 国产精品久久久久久精品电影| 亚洲成人久久爱视频| 激情在线观看视频在线高清| 淫秽高清视频在线观看| 成人美女网站在线观看视频| 日本 av在线| 国产真实乱freesex| 天堂√8在线中文| 国产真实乱freesex| 欧美日韩瑟瑟在线播放| 久久久久久大精品| 欧美一级a爱片免费观看看| 亚洲国产精品sss在线观看| 亚洲最大成人手机在线| .国产精品久久| 精品久久国产蜜桃| 乱码一卡2卡4卡精品| 嫩草影视91久久| 一个人观看的视频www高清免费观看| 99久久精品一区二区三区| 精品久久久久久久久亚洲 | 成人一区二区视频在线观看| 中文字幕人成人乱码亚洲影| 很黄的视频免费| 少妇高潮的动态图| 在线观看美女被高潮喷水网站 | 国产主播在线观看一区二区| 成人av在线播放网站| 中文资源天堂在线| 黄片小视频在线播放| 亚洲成a人片在线一区二区| 日韩中文字幕欧美一区二区| 欧美在线一区亚洲| 在线播放无遮挡| 欧美黄色片欧美黄色片| 精品99又大又爽又粗少妇毛片 | 乱人视频在线观看| 日本一本二区三区精品| 久久99热这里只有精品18| 在线看三级毛片| 国产午夜福利久久久久久| 国产黄色小视频在线观看| 亚洲在线自拍视频| 亚洲国产精品999在线| 欧美日韩乱码在线| av专区在线播放| 久久热精品热| 日本免费a在线| 夜夜躁狠狠躁天天躁| 欧美bdsm另类| 久久草成人影院| 日本免费一区二区三区高清不卡| 免费观看的影片在线观看| 欧美+日韩+精品| 欧美黑人欧美精品刺激| 国产在视频线在精品| 亚洲欧美日韩无卡精品| 国产日本99.免费观看| 动漫黄色视频在线观看| 嫩草影院入口| 亚洲av五月六月丁香网| 国产高清有码在线观看视频| 国产伦在线观看视频一区| 亚洲av成人av| 精品久久久久久久久久免费视频| 国产精品av视频在线免费观看| x7x7x7水蜜桃| 日韩人妻高清精品专区| 男人和女人高潮做爰伦理| 亚洲色图av天堂| 极品教师在线免费播放| 午夜日韩欧美国产| 亚洲av美国av| 亚洲第一电影网av| 在线观看一区二区三区| 欧美黄色淫秽网站| 好男人电影高清在线观看| 欧美性猛交╳xxx乱大交人| 成人性生交大片免费视频hd| 亚洲不卡免费看| 国内精品一区二区在线观看| 亚洲国产精品999在线| 久久久久久久久久黄片| 精品一区二区三区视频在线| 亚洲在线自拍视频| 久久久久久九九精品二区国产| 露出奶头的视频| 一级av片app| 国产黄片美女视频| 宅男免费午夜| 一本一本综合久久| 日韩欧美在线乱码| 美女大奶头视频| 最近最新免费中文字幕在线| 高潮久久久久久久久久久不卡| 日韩欧美国产一区二区入口| 69人妻影院| 国产又黄又爽又无遮挡在线| 91久久精品电影网| 日本a在线网址| 国产激情偷乱视频一区二区| 波多野结衣高清作品| 级片在线观看| 国内揄拍国产精品人妻在线| 精品久久久久久成人av| 久久精品国产自在天天线| 桃色一区二区三区在线观看| 波多野结衣高清无吗| 久久久国产成人免费| 全区人妻精品视频| 美女cb高潮喷水在线观看| 一个人看的www免费观看视频| 性色av乱码一区二区三区2| 天堂av国产一区二区熟女人妻| 亚洲自偷自拍三级| 真实男女啪啪啪动态图| 久久这里只有精品中国| 欧美日韩亚洲国产一区二区在线观看| 日韩欧美精品v在线| 18禁黄网站禁片午夜丰满| 亚洲成人免费电影在线观看| 男人的好看免费观看在线视频| 一个人看的www免费观看视频| 亚洲av第一区精品v没综合| 精品一区二区三区av网在线观看| 91午夜精品亚洲一区二区三区 | 最近在线观看免费完整版| 亚洲不卡免费看| 精品不卡国产一区二区三区| 国产精品自产拍在线观看55亚洲| 蜜桃久久精品国产亚洲av| 久久精品国产亚洲av涩爱 | 搡老妇女老女人老熟妇| 欧美黄色淫秽网站| 日韩欧美在线乱码| 最新中文字幕久久久久| 精品一区二区三区视频在线| 老司机深夜福利视频在线观看| 精品一区二区免费观看| 国内精品久久久久久久电影| 成人国产综合亚洲| 国产欧美日韩精品亚洲av| 国产视频内射| 波野结衣二区三区在线| 亚洲精品粉嫩美女一区| 一级av片app| 亚洲中文字幕一区二区三区有码在线看| or卡值多少钱| 免费黄网站久久成人精品 | 在线免费观看不下载黄p国产 | 亚洲av成人av| 可以在线观看毛片的网站| 国产欧美日韩精品亚洲av| 亚洲在线自拍视频| 男女床上黄色一级片免费看| 九色国产91popny在线| 九九久久精品国产亚洲av麻豆| 国产精品久久电影中文字幕| 国产免费av片在线观看野外av| 好看av亚洲va欧美ⅴa在| 国产精品爽爽va在线观看网站| 波多野结衣巨乳人妻| 久久精品影院6| 亚洲av成人精品一区久久| 久久久久久久久久黄片| 看免费av毛片| 麻豆国产97在线/欧美| 亚洲国产精品久久男人天堂| 黄色日韩在线| 成人一区二区视频在线观看| 自拍偷自拍亚洲精品老妇| 成人无遮挡网站| 成人一区二区视频在线观看| av国产免费在线观看| 国产亚洲精品久久久com| 91九色精品人成在线观看| 国产精品免费一区二区三区在线| 69av精品久久久久久| 免费在线观看日本一区| 亚洲美女黄片视频| 男人的好看免费观看在线视频| 欧美高清成人免费视频www| 日韩欧美一区二区三区在线观看| 在线免费观看不下载黄p国产 | 亚洲va日本ⅴa欧美va伊人久久| 亚洲精品粉嫩美女一区| 国内久久婷婷六月综合欲色啪| 嫩草影视91久久| 亚洲欧美清纯卡通| 国产探花极品一区二区| 变态另类丝袜制服| 国内久久婷婷六月综合欲色啪| 久久欧美精品欧美久久欧美| 性色avwww在线观看| 特级一级黄色大片| 国产av不卡久久| 直男gayav资源| 精品久久久久久久久亚洲 | 757午夜福利合集在线观看| 日本a在线网址| 亚洲最大成人中文| 亚洲精品久久国产高清桃花| 欧美三级亚洲精品| bbb黄色大片| 日本五十路高清| 亚洲成a人片在线一区二区| 亚洲成人精品中文字幕电影| 亚洲美女搞黄在线观看 | 美女免费视频网站| 高清在线国产一区| 深夜精品福利| 国产伦精品一区二区三区四那| 国产亚洲欧美在线一区二区| 国产亚洲精品av在线| 国产黄色小视频在线观看| 黄色一级大片看看| 亚洲五月婷婷丁香| 每晚都被弄得嗷嗷叫到高潮| 国产视频内射| 久久精品夜夜夜夜夜久久蜜豆| 又爽又黄a免费视频| 国产成年人精品一区二区| 日日夜夜操网爽| 中文字幕熟女人妻在线| 两性午夜刺激爽爽歪歪视频在线观看| 色精品久久人妻99蜜桃| 人妻夜夜爽99麻豆av| 三级毛片av免费| 国产av在哪里看| 欧美3d第一页| 久久99热这里只有精品18| 日本一本二区三区精品| a级毛片a级免费在线| 特级一级黄色大片| 国模一区二区三区四区视频| 日韩欧美 国产精品| 免费在线观看成人毛片| 精华霜和精华液先用哪个| 欧美国产日韩亚洲一区| 天堂√8在线中文| 久久久久久久久中文| 观看美女的网站| 亚洲经典国产精华液单 | 免费看a级黄色片| 色综合欧美亚洲国产小说| 国产高清激情床上av| 一本精品99久久精品77| 久久久久久国产a免费观看| 一本综合久久免费| 毛片一级片免费看久久久久 | 国产野战对白在线观看| 国产成人啪精品午夜网站| 久久久久久久久久黄片| 欧美高清成人免费视频www| 69人妻影院| 女同久久另类99精品国产91| 国产美女午夜福利| 日韩中字成人| 久久人妻av系列| 国产一区二区亚洲精品在线观看| 国产成+人综合+亚洲专区| 天天躁日日操中文字幕| 两个人的视频大全免费| 两人在一起打扑克的视频| 五月伊人婷婷丁香| 中文在线观看免费www的网站| 长腿黑丝高跟| 99久久99久久久精品蜜桃| 性色av乱码一区二区三区2| 又紧又爽又黄一区二区| 国产大屁股一区二区在线视频| 我的女老师完整版在线观看| 麻豆一二三区av精品| 偷拍熟女少妇极品色| 99久久九九国产精品国产免费| 免费观看的影片在线观看| 欧美一区二区精品小视频在线| 国内精品久久久久久久电影| 久久亚洲精品不卡| 午夜精品一区二区三区免费看| 国产精品自产拍在线观看55亚洲| 亚洲经典国产精华液单 | 久久国产乱子免费精品| 亚洲第一电影网av| 琪琪午夜伦伦电影理论片6080| 免费看光身美女| 一区二区三区免费毛片| 在线观看av片永久免费下载| 成人亚洲精品av一区二区| 男女床上黄色一级片免费看| 一本综合久久免费| 午夜免费激情av| 国产爱豆传媒在线观看| 国产精品伦人一区二区| 91麻豆精品激情在线观看国产| 真人做人爱边吃奶动态| 最近中文字幕高清免费大全6 | 嫩草影视91久久| 成人亚洲精品av一区二区| 精品午夜福利在线看| 免费av观看视频| 天天躁日日操中文字幕| 精品人妻偷拍中文字幕| 欧美成人免费av一区二区三区| 亚州av有码| 91狼人影院| 99riav亚洲国产免费| 国产精品久久久久久亚洲av鲁大| 日本在线视频免费播放| 久久久久久大精品| 在线观看66精品国产| 别揉我奶头~嗯~啊~动态视频| 美女大奶头视频| 国产精品一区二区免费欧美| 日韩欧美免费精品| 成年女人看的毛片在线观看| 首页视频小说图片口味搜索| 亚洲av免费在线观看| 91狼人影院| 俺也久久电影网| av在线观看视频网站免费| 欧洲精品卡2卡3卡4卡5卡区| 欧美日韩中文字幕国产精品一区二区三区| 欧美激情久久久久久爽电影| 精品欧美国产一区二区三| 国产真实伦视频高清在线观看 | 色哟哟哟哟哟哟| 亚洲精品456在线播放app | 简卡轻食公司| avwww免费| 免费在线观看日本一区| 欧美zozozo另类| 日韩欧美 国产精品| 成人亚洲精品av一区二区| 亚洲熟妇熟女久久| 欧美激情久久久久久爽电影| 简卡轻食公司| 国产精品av视频在线免费观看| 真人做人爱边吃奶动态| 中亚洲国语对白在线视频| 亚洲精品亚洲一区二区| 国产成人啪精品午夜网站| 久久久久久国产a免费观看| 亚洲色图av天堂| 狂野欧美白嫩少妇大欣赏| 久久精品久久久久久噜噜老黄 | 亚洲自偷自拍三级| 国产成+人综合+亚洲专区| 禁无遮挡网站| 无遮挡黄片免费观看| 观看免费一级毛片| 99国产综合亚洲精品| 夜夜看夜夜爽夜夜摸| 欧美黄色淫秽网站| 亚洲熟妇中文字幕五十中出| 久久久久国产精品人妻aⅴ院| 简卡轻食公司| 久久婷婷人人爽人人干人人爱| 日韩欧美在线乱码| 男插女下体视频免费在线播放| 国产av在哪里看| 又黄又爽又刺激的免费视频.| 一夜夜www| 成人毛片a级毛片在线播放| 欧美高清性xxxxhd video| 综合色av麻豆| 在线观看66精品国产| 午夜激情欧美在线| 黄色女人牲交| 亚洲在线观看片| 欧美性猛交黑人性爽| 国产又黄又爽又无遮挡在线| 中文字幕av在线有码专区| 亚洲熟妇熟女久久| 亚洲成av人片在线播放无| 国产精品自产拍在线观看55亚洲| 国产精品野战在线观看| 国产精品自产拍在线观看55亚洲| 国产亚洲精品av在线| 日韩欧美免费精品| 天天一区二区日本电影三级| 免费搜索国产男女视频| 中文字幕高清在线视频| 亚洲欧美日韩东京热| av在线天堂中文字幕| 成人永久免费在线观看视频| 女人十人毛片免费观看3o分钟| av福利片在线观看| 欧美乱妇无乱码| 免费看a级黄色片| 亚洲avbb在线观看| 亚洲精品日韩av片在线观看| 成熟少妇高潮喷水视频| 村上凉子中文字幕在线| 男人狂女人下面高潮的视频| 亚洲成人精品中文字幕电影| 久久草成人影院| 他把我摸到了高潮在线观看| 久久久久久久亚洲中文字幕 | 精品午夜福利视频在线观看一区| 久久久成人免费电影| 天堂网av新在线| 亚洲av电影在线进入| 免费看光身美女| 俺也久久电影网| 国产三级中文精品| 欧美极品一区二区三区四区| 欧美午夜高清在线| a级一级毛片免费在线观看| 一夜夜www| 国产精品久久久久久久电影| 亚洲成人久久爱视频| av国产免费在线观看| 老女人水多毛片| 成年女人永久免费观看视频| 亚洲国产精品合色在线| 毛片女人毛片| 久久欧美精品欧美久久欧美| 亚洲avbb在线观看| 长腿黑丝高跟|