• <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
    伴 娘
    免费看美女性在线毛片视频| 亚洲精品乱码久久久久久按摩| 国产三级在线视频| 国内精品一区二区在线观看| 在线播放无遮挡| 久久精品国产亚洲av天美| 日本免费一区二区三区高清不卡| 亚洲人成网站高清观看| 超碰97精品在线观看| 桃色一区二区三区在线观看| 成人性生交大片免费视频hd| 狂野欧美白嫩少妇大欣赏| 老司机影院成人| 国产精品美女特级片免费视频播放器| 麻豆成人午夜福利视频| 少妇熟女欧美另类| 午夜激情欧美在线| 久久99热6这里只有精品| 亚洲乱码一区二区免费版| 91aial.com中文字幕在线观看| 欧美一区二区国产精品久久精品| 国产成人aa在线观看| 亚洲人与动物交配视频| 亚洲av福利一区| 精品国内亚洲2022精品成人| 免费人成在线观看视频色| 国产伦在线观看视频一区| 黄色配什么色好看| 午夜福利在线观看免费完整高清在| ponron亚洲| 亚洲怡红院男人天堂| 中文精品一卡2卡3卡4更新| 国产亚洲午夜精品一区二区久久 | 亚洲在线观看片| 欧美成人午夜免费资源| 亚洲av中文字字幕乱码综合| 国产精品人妻久久久久久| 亚洲av中文字字幕乱码综合| 国产又黄又爽又无遮挡在线| 免费观看精品视频网站| 赤兔流量卡办理| 成人欧美大片| 大又大粗又爽又黄少妇毛片口| 国产在线男女| 久久精品国产亚洲av涩爱| 免费av不卡在线播放| 九九在线视频观看精品| 久久热精品热| 国产三级在线视频| 男女下面进入的视频免费午夜| 久久久久久国产a免费观看| 禁无遮挡网站| 久久久久久久久久久丰满| 看十八女毛片水多多多| 久久99热这里只有精品18| 亚洲精品影视一区二区三区av| 亚洲久久久久久中文字幕| 国产精品.久久久| 色噜噜av男人的天堂激情| 日日干狠狠操夜夜爽| 晚上一个人看的免费电影| 亚洲国产精品合色在线| 亚洲精品成人久久久久久| 一个人免费在线观看电影| 麻豆av噜噜一区二区三区| 国产伦理片在线播放av一区| 看非洲黑人一级黄片| 岛国在线免费视频观看| 久久精品国产鲁丝片午夜精品| 黄色日韩在线| 色吧在线观看| 亚洲成人久久爱视频| 99久久精品一区二区三区| 国内少妇人妻偷人精品xxx网站| 久久久久久国产a免费观看| 热99re8久久精品国产| 久久久久久久久中文| 午夜福利在线在线| 日本av手机在线免费观看| 国产精品伦人一区二区| 黄色日韩在线| 日韩欧美三级三区| 久久久国产成人免费| 成人一区二区视频在线观看| av在线观看视频网站免费| 深爱激情五月婷婷| 国产在线一区二区三区精 | 97在线视频观看| 欧美xxxx性猛交bbbb| 日韩在线高清观看一区二区三区| 村上凉子中文字幕在线| 亚洲国产精品成人综合色| 国产精品久久电影中文字幕| av在线蜜桃| 亚洲欧洲日产国产| 欧美成人精品欧美一级黄| 岛国毛片在线播放| 免费不卡的大黄色大毛片视频在线观看 | 成年女人永久免费观看视频| 青春草国产在线视频| 熟妇人妻久久中文字幕3abv| 美女国产视频在线观看| 午夜精品在线福利| 中文资源天堂在线| 26uuu在线亚洲综合色| 国产成人a∨麻豆精品| 乱码一卡2卡4卡精品| 寂寞人妻少妇视频99o| 九色成人免费人妻av| 少妇人妻精品综合一区二区| 久久久久久久久久久免费av| 亚洲av.av天堂| 2021少妇久久久久久久久久久| 夜夜爽夜夜爽视频| 久久精品91蜜桃| 免费观看a级毛片全部| 国产一区二区亚洲精品在线观看| 国产男人的电影天堂91| 国产男人的电影天堂91| 18+在线观看网站| 搞女人的毛片| 插逼视频在线观看| 国产 一区精品| 91av网一区二区| 国产成人精品婷婷| 成人漫画全彩无遮挡| 一夜夜www| 亚洲人与动物交配视频| 国产大屁股一区二区在线视频| 国产午夜精品一二区理论片| .国产精品久久| 国产av码专区亚洲av| 少妇被粗大猛烈的视频| 国国产精品蜜臀av免费| 一区二区三区乱码不卡18| 色视频www国产| 欧美97在线视频| 亚洲精品乱码久久久久久按摩| 在线观看美女被高潮喷水网站| 中文字幕制服av| 熟女电影av网| 欧美性猛交╳xxx乱大交人| 嘟嘟电影网在线观看| 免费av不卡在线播放| АⅤ资源中文在线天堂| 午夜福利成人在线免费观看| 日本免费一区二区三区高清不卡| 一个人观看的视频www高清免费观看| 丝袜喷水一区| 国产精品无大码| 午夜福利在线观看吧| 国产成人a区在线观看| 91午夜精品亚洲一区二区三区| 欧美性感艳星| 国产成人精品一,二区| 国产三级中文精品| 精品久久久久久久末码| 欧美97在线视频| 亚洲乱码一区二区免费版| 六月丁香七月| 国产成人午夜福利电影在线观看| 亚洲av.av天堂| 国产精品美女特级片免费视频播放器| 波野结衣二区三区在线| 九九热线精品视视频播放| 极品教师在线视频| 男的添女的下面高潮视频| 综合色丁香网| 一级毛片aaaaaa免费看小| 国产一区二区三区av在线| 欧美xxxx黑人xx丫x性爽| 青春草视频在线免费观看| 黄片无遮挡物在线观看| 日本欧美国产在线视频| 精品酒店卫生间| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 午夜日本视频在线| 性插视频无遮挡在线免费观看| 精品久久久噜噜| 亚洲成人精品中文字幕电影| 国模一区二区三区四区视频| 久久久欧美国产精品| 少妇丰满av| 亚洲四区av| 日韩成人伦理影院| 久久久久久久久久黄片| 99热这里只有是精品在线观看| 黄色日韩在线| 免费大片18禁| 99热精品在线国产| 精品一区二区免费观看| 久久久国产成人精品二区| 国产精品久久久久久久电影| av在线播放精品| 久久人妻av系列| 日本五十路高清| 好男人视频免费观看在线| 99久久精品一区二区三区| 欧美bdsm另类| 青春草国产在线视频| 人体艺术视频欧美日本| 中文字幕熟女人妻在线| 黄色欧美视频在线观看| 亚洲欧美成人综合另类久久久 | 少妇人妻一区二区三区视频| 精品人妻熟女av久视频| 精品国产一区二区三区久久久樱花 | 2021少妇久久久久久久久久久| 国产成人一区二区在线| 亚洲国产欧美人成| 中文字幕制服av| 一区二区三区高清视频在线| 国产高清有码在线观看视频| 成人午夜精彩视频在线观看| 色噜噜av男人的天堂激情| 成人av在线播放网站| 亚洲最大成人中文| АⅤ资源中文在线天堂| 亚洲精品一区蜜桃| 少妇人妻一区二区三区视频| 狂野欧美白嫩少妇大欣赏| 国产精品1区2区在线观看.| 97在线视频观看| 直男gayav资源| 国产三级中文精品| 99在线视频只有这里精品首页| 久久鲁丝午夜福利片| 久久精品影院6| 亚洲av电影不卡..在线观看| 人人妻人人澡欧美一区二区| 精品一区二区三区视频在线| 在线a可以看的网站| 狂野欧美激情性xxxx在线观看| 亚洲精品一区蜜桃| 深夜a级毛片| 好男人在线观看高清免费视频| 久久鲁丝午夜福利片| 亚洲欧洲日产国产| 欧美一级a爱片免费观看看| 精华霜和精华液先用哪个| 亚洲国产色片| 久久精品91蜜桃| 两个人视频免费观看高清| 日韩一区二区视频免费看| 99久久精品热视频| 午夜福利成人在线免费观看| 亚洲乱码一区二区免费版| 老司机影院成人| 人人妻人人看人人澡| 欧美日韩一区二区视频在线观看视频在线 | 日韩精品青青久久久久久| 3wmmmm亚洲av在线观看| 青青草视频在线视频观看| 亚洲精品一区蜜桃| 日本黄大片高清| 男女边吃奶边做爰视频| 男人和女人高潮做爰伦理| 日韩欧美在线乱码| eeuss影院久久| 97热精品久久久久久| 欧美97在线视频| 看免费成人av毛片| 亚洲av免费在线观看| 国产美女午夜福利| 午夜精品一区二区三区免费看| 国产高清不卡午夜福利| 黄片wwwwww| 日韩成人av中文字幕在线观看| 麻豆精品久久久久久蜜桃| 国产 一区 欧美 日韩| 岛国毛片在线播放| 久久精品久久久久久噜噜老黄 | 欧美高清成人免费视频www| 蜜臀久久99精品久久宅男| 久久精品综合一区二区三区| 男插女下体视频免费在线播放| 精品久久久久久久久久久久久| 中文字幕亚洲精品专区| 欧美高清成人免费视频www| 人妻夜夜爽99麻豆av| 青春草亚洲视频在线观看| 女的被弄到高潮叫床怎么办| 日韩,欧美,国产一区二区三区 | 日本三级黄在线观看| 免费观看在线日韩| 精品免费久久久久久久清纯| 免费看美女性在线毛片视频| 亚洲在线自拍视频| 免费看a级黄色片| 18禁裸乳无遮挡免费网站照片| 久久久久久久久中文| 亚洲欧美成人综合另类久久久 | 亚洲精品乱码久久久久久按摩| 亚洲五月天丁香| 大又大粗又爽又黄少妇毛片口| 亚洲精品,欧美精品| 亚洲欧美日韩卡通动漫| 内地一区二区视频在线| 国产免费男女视频| 五月伊人婷婷丁香| 免费观看性生交大片5| 99热全是精品| 白带黄色成豆腐渣| av免费在线看不卡| 久久婷婷人人爽人人干人人爱| 亚洲最大成人手机在线| 精品久久久噜噜| 国产av码专区亚洲av| 亚洲精品色激情综合| 亚洲欧美一区二区三区国产| 久久久久国产网址| 国产精品一区二区三区四区久久| 天堂影院成人在线观看| 七月丁香在线播放| 尤物成人国产欧美一区二区三区| 国产伦精品一区二区三区视频9| 如何舔出高潮| av国产免费在线观看| 成年免费大片在线观看| 国产精品国产三级专区第一集| 欧美色视频一区免费| 日本免费在线观看一区| 一本久久精品| av.在线天堂| 99国产精品一区二区蜜桃av| 欧美不卡视频在线免费观看| 欧美zozozo另类| 搞女人的毛片| 91在线精品国自产拍蜜月| 国产成人一区二区在线| 观看美女的网站| 精品欧美国产一区二区三| 色综合色国产| 少妇猛男粗大的猛烈进出视频 | 日韩欧美在线乱码| 亚洲精品国产成人久久av| 免费人成在线观看视频色| 97超视频在线观看视频| 欧美一级a爱片免费观看看| 美女黄网站色视频| 国产成人精品一,二区| av免费在线看不卡| 国产成人a区在线观看| 丰满人妻一区二区三区视频av| 亚洲不卡免费看| 日韩精品有码人妻一区| 亚洲国产精品久久男人天堂| 国产成人91sexporn| 亚洲精品乱码久久久v下载方式| 欧美bdsm另类| 99久久人妻综合| 26uuu在线亚洲综合色| 欧美成人免费av一区二区三区| 国产一区二区三区av在线| 国产亚洲一区二区精品| 91av网一区二区| 久久草成人影院| 亚洲成人精品中文字幕电影| 国产精品蜜桃在线观看| 国产一区二区在线观看日韩| 中文字幕免费在线视频6| 黄片wwwwww| 午夜福利成人在线免费观看| 久久这里只有精品中国| 直男gayav资源| 国产伦精品一区二区三区视频9| 国产精品久久电影中文字幕| 美女国产视频在线观看| 亚洲在线观看片| 亚洲国产欧美在线一区| 免费一级毛片在线播放高清视频| 国产精品爽爽va在线观看网站| 国产在线男女| 免费看美女性在线毛片视频| 看免费成人av毛片| 18禁在线无遮挡免费观看视频| 女人久久www免费人成看片 | 久久久久久九九精品二区国产| 成人性生交大片免费视频hd| 国产乱人视频| 色视频www国产| 国产精品无大码| 久久草成人影院| 五月伊人婷婷丁香| 久久久久久久久久久丰满| 亚洲国产精品专区欧美| 中文在线观看免费www的网站| 亚洲av.av天堂| 久久人妻av系列| 国产高清视频在线观看网站| 亚洲av中文字字幕乱码综合| 偷拍熟女少妇极品色| 欧美人与善性xxx| 国产黄片视频在线免费观看| 伦理电影大哥的女人| 亚洲精品影视一区二区三区av| 国产高清不卡午夜福利| 男人和女人高潮做爰伦理| 亚洲精品aⅴ在线观看| 国产精品一区www在线观看| 国产精品综合久久久久久久免费| 中文字幕精品亚洲无线码一区| 91午夜精品亚洲一区二区三区| 国产乱来视频区| 国产真实伦视频高清在线观看| 99久国产av精品| 成人国产麻豆网| 精品人妻一区二区三区麻豆| 欧美精品国产亚洲| 亚洲成av人片在线播放无| 少妇裸体淫交视频免费看高清| 亚洲在久久综合| 日日摸夜夜添夜夜爱| 国产精品1区2区在线观看.| 性色avwww在线观看| 永久网站在线| 非洲黑人性xxxx精品又粗又长| 欧美又色又爽又黄视频| 女的被弄到高潮叫床怎么办| 国产私拍福利视频在线观看| av专区在线播放| 日韩强制内射视频| 赤兔流量卡办理| 亚洲精品乱久久久久久| 国产色爽女视频免费观看| 欧美精品国产亚洲| 欧美日本视频| 免费大片18禁| 国产精品一区二区在线观看99 | 国产精品日韩av在线免费观看| av专区在线播放| 丰满人妻一区二区三区视频av| 亚洲av电影不卡..在线观看| 我的女老师完整版在线观看| a级一级毛片免费在线观看| 综合色av麻豆| 黄色一级大片看看| 国产精品久久视频播放| 日本黄大片高清| 日韩 亚洲 欧美在线| 欧美日本亚洲视频在线播放| 精品99又大又爽又粗少妇毛片| 久久久久久伊人网av| 国产色婷婷99| 欧美性感艳星| 国产精品一区二区在线观看99 | 建设人人有责人人尽责人人享有的 | 色哟哟·www| 亚洲精华国产精华液的使用体验| 一个人免费在线观看电影| 久久久久久久久久久免费av| 精品欧美国产一区二区三| 男人舔奶头视频| 啦啦啦啦在线视频资源| 国产一区二区在线观看日韩| 真实男女啪啪啪动态图| 精品久久久久久久久亚洲| 又爽又黄a免费视频| 日韩三级伦理在线观看| 亚洲成av人片在线播放无| 两个人视频免费观看高清| 欧美日韩综合久久久久久| 午夜免费男女啪啪视频观看| 97在线视频观看| 日韩国内少妇激情av| 亚洲高清免费不卡视频| 国产精品综合久久久久久久免费| 男人狂女人下面高潮的视频| 国产私拍福利视频在线观看| 亚洲人成网站在线观看播放| av视频在线观看入口| 日本免费在线观看一区| 久久人人爽人人片av| 26uuu在线亚洲综合色| 男女那种视频在线观看| 青春草亚洲视频在线观看| 成年版毛片免费区| 国产片特级美女逼逼视频| 自拍偷自拍亚洲精品老妇| 久久久久久久国产电影| 日韩三级伦理在线观看| 久久亚洲国产成人精品v| 69av精品久久久久久| 欧美日韩一区二区视频在线观看视频在线 | 日韩一本色道免费dvd| 国产一级毛片在线| 高清日韩中文字幕在线| 成人午夜精彩视频在线观看| 国产黄片视频在线免费观看| 插逼视频在线观看| 亚洲欧美一区二区三区国产| 国产 一区精品| 老女人水多毛片| 国产亚洲午夜精品一区二区久久 | 国产淫语在线视频| 观看美女的网站| 老司机影院毛片| 大香蕉97超碰在线| 成人毛片60女人毛片免费| 超碰97精品在线观看| 久久99热这里只有精品18| 亚洲丝袜综合中文字幕| av视频在线观看入口| 永久网站在线| 超碰97精品在线观看| 我的女老师完整版在线观看| kizo精华| 欧美成人一区二区免费高清观看| 久久久久久九九精品二区国产| 国产精品久久久久久久久免| 国产高清视频在线观看网站| 精品少妇黑人巨大在线播放 | 日韩人妻高清精品专区| 中文字幕免费在线视频6| 又粗又硬又长又爽又黄的视频| 色5月婷婷丁香| 日韩欧美精品免费久久| 国产精品国产三级专区第一集| av线在线观看网站| 一级毛片电影观看 | 美女xxoo啪啪120秒动态图| 老司机影院毛片| 午夜激情欧美在线| 亚洲伊人久久精品综合 | 久久久久久伊人网av| 国产又色又爽无遮挡免| 国模一区二区三区四区视频| 日日啪夜夜撸| 久久精品国产鲁丝片午夜精品| 成人特级av手机在线观看| 国语自产精品视频在线第100页| 亚洲国产精品成人久久小说| 看非洲黑人一级黄片| 啦啦啦观看免费观看视频高清| 日本av手机在线免费观看| 亚洲欧美精品自产自拍| 国产69精品久久久久777片| 国产精品国产三级国产av玫瑰| 免费av观看视频| 91在线精品国自产拍蜜月| 九草在线视频观看| 99九九线精品视频在线观看视频| 成人无遮挡网站| 国产熟女欧美一区二区| 麻豆成人午夜福利视频| 午夜福利网站1000一区二区三区| 日韩 亚洲 欧美在线| 亚洲精品,欧美精品| 欧美成人精品欧美一级黄| 婷婷六月久久综合丁香| 亚洲人成网站高清观看| 欧美性猛交黑人性爽| 久久精品国产亚洲av涩爱| videossex国产| 国内少妇人妻偷人精品xxx网站| 久久久久性生活片| 久久国产乱子免费精品| 亚洲天堂国产精品一区在线| 精品久久国产蜜桃| 麻豆成人av视频| 日本爱情动作片www.在线观看| 国产白丝娇喘喷水9色精品| 久久午夜福利片| 国产极品精品免费视频能看的| av视频在线观看入口| 插逼视频在线观看| 亚洲欧美日韩高清专用| 插逼视频在线观看| 亚洲最大成人手机在线| 91久久精品电影网| 国产乱来视频区| 美女大奶头视频| 老司机福利观看| 看非洲黑人一级黄片| 国产在线一区二区三区精 | or卡值多少钱| 欧美日本视频| 日日干狠狠操夜夜爽| 色网站视频免费| 国产精品久久久久久精品电影小说 | 久久鲁丝午夜福利片| 97热精品久久久久久| 国产精品一区二区性色av| 亚洲人与动物交配视频| 国产真实伦视频高清在线观看| 亚洲不卡免费看| 五月玫瑰六月丁香| 国产又色又爽无遮挡免| 天美传媒精品一区二区| 伦理电影大哥的女人| 午夜福利在线观看吧| 国产精品久久久久久久电影| 久久久国产成人免费| 日韩精品青青久久久久久| 特大巨黑吊av在线直播| 久久久久久伊人网av| 成人欧美大片| 校园人妻丝袜中文字幕| 国产黄片美女视频| 国产精品嫩草影院av在线观看| 日韩欧美 国产精品| 热99re8久久精品国产| 精品人妻视频免费看| 精品久久久久久成人av| 九九久久精品国产亚洲av麻豆| 亚洲av成人精品一区久久| 国产精品日韩av在线免费观看| 搞女人的毛片| 一级毛片我不卡| 男女啪啪激烈高潮av片| 久久久国产成人精品二区| 99久久精品热视频| 国产免费又黄又爽又色| 嫩草影院精品99| 亚洲在久久综合|