999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

Simulations of energetic alpha particle loss in the presence of toroidal field ripple in the CFETR tokamak

2022-09-06 13:04:22YingfengXU徐穎峰DebingZHANG張德兵JialeCHEN陳佳樂andFangchuanZHONG鐘方川
Plasma Science and Technology 2022年10期

Yingfeng XU (徐穎峰),Debing ZHANG (張德兵),Jiale CHEN (陳佳樂) and Fangchuan ZHONG (鐘方川)

1 College of Science,Donghua University,Shanghai 201620,People's Republic of China

2 Member of Magnetic Confinement Fusion Research Centre,Ministry of Education,Shanghai 201620,People's Republic of China

3 Department of Physics,East China University of Science and Technology,Shanghai 200237,People's Republic of China

4 Institute of Plasma Physics,Chinese Academy of Sciences,Hefei 230031,People's Republic of China

Abstract Energetic alpha particle losses with the toroidal field ripple and the Coulomb collision in the CFETR tokamak have been simulated by using the orbit-following code GYCAVA for the steady-state and hybrid scenarios.The effects of the outer boundary and the ripple amplitude on alpha particle losses have been investigated.The loss fractions and heat loads of alpha particles in the hybrid scenario are much smaller than those in the steady-state scenario for a significant ripple amplitude.Some alpha particles in the plasma core are lost due to the ripple stochastic transport for a large ripple amplitude parameter.The heat loads with the last closed flux surface boundary are different from those with the wall boundary for the CFETR tokamak,which can be explained by typical alpha particle orbits.Discrete heat load spots have been observed in alpha particle loss simulations,which is due to the ripple well loss.The transition of the lost alpha particle behavior from the ripple stochastic diffusion to the ripple well trapping has been identified in our CFETR simulations.The Coulomb collision effect is responsible for this transition.

Keywords: CFETR,toroidal field ripple,orbit-following code,alpha particle loss

1.Introduction

The Chinese Fusion Engineering Test Reactor(CFETR)[1-6]is the next tokamak fusion device in the Chinese magnetic fusion energy development roadmap.The missions of CFETR include generating burning D-T plasmas in the steady-state and testing the self-sustainable burning state with a large fusion gain,Q=20-30.Energetic alpha particles produced by D-T fusion reaction play a significant role in magnetic confinement fusion plasmas.Fusion alpha particles can heat plasmas through the Coulomb collision with background ions and electrons.For the self-sustainable burning plasmas with Q=20-30 on CFETR,heating by alpha particles is the dominant plasma heating [5].Good confinement of alpha particles is beneficial to the alpha particle heating.Losses of alpha particles,which can be induced by the electromagnetic perturbations,will damage the plasma facing components of a tokamak and reduce the alpha particle heating efficiency of alpha particle.The toroidal field ripple,which comes from the discreteness of toroidal field coils,is one of the major toroidal asymmetric magnetic perturbations which can induce significant losses of alpha particles.

The pioneering numerical work of ripple-induced loss of alpha particles has been performed by the OMFC code [7].Numerical simulations and experimental studies of alpha particle losses induced by the ripple field have been performed in many tokamaks,such as TFTR [8-13],ITER [14-16],CFETR[17-19]and DEMO [20].In the previous simulations of alpha particle losses on CFETR [17-19],the orbit-following code ORBIT [21]was used and the last closed flux surface (LCFS)was chosen as the outer boundary.However,the realistic outer boundary of a tokamak is the first wall.The different choices of the outer boundary in orbit-following simulations may impact the loss fraction of alpha particles and the corresponding distribution of heat load.

A Monte Carlo orbit-following code GYCAVA [22,23]has been developed in the basis of the modern gyrokinetic theory [24]for the test particle simulation.Recently,another particle tracking code PTC [25]has been developed for studying the alpha particle behavior in a tokamak,especially for CFETR.The orbit code GYCAVA has been used for investigating the neutral beam injection (NBI) ion loss in the presence of the resonant magnetic perturbations [26]or the ripple field[27].In our previous numerical results of NBI ion loss on EAST computed by the orbit code GYCAVA and the NBI code TGCO [28],the outer boundary effect can impact the loss fraction and the heat load of NBI ions[29].Therefore,it is necessary to study the boundary effect on the alpha particle loss in the CFETR tokamak.The coordinates adopted in the orbit code GYCAVA can be chosen as the canonical variables based on the magnetic flux coordinates or the cylindrical coordinates.In contrast to the coordinates based on the magnetic flux coordinates,using the cylindrical coordinates in orbit-following codes can avoid the singularity problem of the safety factor at X points of the LCFS,then test charged particles in orbit-following simulations can move across the LCFS and reach the first wall.In the GYCAVA code,only the LCFS can be chosen as the outer boundary if the canonical variables are adopted,which is due to the limitation of the magnetic flux coordinates.However,both of the LCFS and the first wall can be chosen as the outer boundary if the cylindrical coordinates are adopted.

In this work,simulation results of alpha particle losses in the CFETR steady-state and hybrid scenarios with the toroidal field ripple and the Coulomb collision computed by the orbitfollowing code GYCAVA are presented.The LCFS or the first wall has been used as the loss/outer boundary in our alpha particle loss simulations.In order to compute to alpha particle behaviors on the wall,the cylindrical coordinates were used in our all simulations.The ripple effect and the outer boundary effect on loss fractions,loss regions and heat loads of alpha particles have been investigated.

The remaining part of the paper is organized as follows.In section 2,Hamilton’s equations of motion with ripple and the Coulomb collision for an alpha particle adopted in the orbit code GYCAVA are given.The magnetic equilibriums of CFETR in the steady-state and hybrid scenarios,their plasma profiles,the birth profiles of alpha particles and the toroidal field ripple are presented in section 3.The numerical results of alpha particle losses with ripple and collision in the CFETR tokamak are given in section 4.The main conclusions are given in section 5.

2.Hamilton’s equations of motion with ripple and Coulomb collision used in GYCAVA

The GYCAVA code was used for studying the toroidal field ripple effect on alpha particle losses.In contrast to the version of GYCAVA used in [27],in which some approximations related to ripple have been adopted,more accurate Hamilton’s equations with ripple are adopted in the present version of GYCAVA.The guiding-center Hamilton’s equations with ripple are introduced below.

The total magnetic field B includes the axisymmetric equilibrium magnetic field B0and the toroidal ripple field δBrip,that is

The equilibrium magnetic field is written as

Here,ψ and φ are the poloidal magnetic flux and the toroidal angle,respectively.

For simplicity,only the toroidal component of the ripple field is kept in the GYCAVA code,because the ripple effect on energetic ion losses is mainly due to the toroidal component of the ripple field[30,31].The toroidal ripple field can be expressed as

Here,δ is the toroidal field ripple parameter and N is the number of the toroidal field coils.δ is a function of R,Z and is independent of φ.

Then the total magnetic field including the ripple field can also be expressed as

The guiding-center Hamilton’s equations with ripple in the coordinates (X,v‖) can be expressed as Here,H is the Hamiltonian and the subscript s denotes the species of a particle.msis the mass and esis the electric charge.B*is defined asB*=B+×bandB‖*=B* ·b.b=B/B and B=|B0+δBrip|.

The Hamilton’s equations(8)and(9)can be rewritten in terms of the cylindrical coordinates as

Here,Jijare the components of the Poisson matrix.Note that the ripple effect is included in the Poisson matrix and the Hamiltonian.

The Coulomb collision between energetic ions and background particles was included in the GYCAVA code.The Coulomb collision contains the slowing down part and the pitch angle scattering part.The slowing down of energetic ions due to drag of background particles can be written as

with νsis the slowing down collision rate,expressed as

Zmis defined asFor energetic alpha particles in D-T plasmas,Zm?1.67,when the deuterium and tritium densities are equal.vcis the critical velocity,which can be expressed as

Here,Zf,mfare mass and charge number of an energetic ion,respectively.The subscript e denotes the electron.

The pitch angle scattering of an energetic ion is governed by

Here,Δt is the time step and λ=v‖/v is the particle pitch.νdis the pitch angle scattering rate,written as

Here,Zeffis the effective charge number.In our orbit-following simulations on CFETR,the default effective charge number is chosen as Zeff=2.45.Note that the slowing down collision rate is independent of Zeffbut the pitch angle scattering rate linearly depends on Zeff.

3.Simulation setup

The magnetic configurations,plasma profiles,initial distributions of energetic alpha particles on CFETR are given in this section.These data were used as input data in our alpha loss simulations performed by the GYCAVA code.

Figure 1 shows the magnetic configurations in the CFETR steady-state and hybrid scenarios.The directions of the plasma current and the equilibrium magnetic field are both counter-clockwise from the top view.It means that the magnetic drift of an alpha particle is in the upward direction.The minor radius is a=2.2 m.The magnetic field and the major radius at the axis are=5.7 T and=8.0 m for the steady-state scenario,and=6.1 T and=7.6 m for the hybrid scenario.The plasma currents for the steady-state and hybrid scenarios are=11 MAand=13 MA,respectively.The safety factor profiles in the CFETR steady-state and hybrid scenarios are shown in figure 2.ρtis the square root of the normalized toroidal magnetic flux.We can see from figure 2 that the safety factor in the hybrid scenario is smaller than that in the steady-state scenario.The safety factors at the axis for the steady-state and hybrid scenarios are=3.9and=1.4,respectively.The safety factors at 95% of the poloidal magnetic flux are=7.3and=5.8,respectively.

Figure 1.The magnetic configurations (magenta) in the CFETR steady-state (a) and hybrid (b) scenarios.The black lines denote the shape of the first wall and the blue lines denote the LCFS.

Figure 2.The safety factor profiles in the CFETR steady-state(black) and hybrid (blue) scenarios.

Figure 3 shows the density and temperature profiles,including the densities of electron and deuterium,and the temperatures of electron and bulk ion in the CFETR steadystate [32]and hybrid [33]scenarios.The tritium density profile is the same as the deuterium one.These profiles were obtained by integrated modeling described in [32,33].

Figure 3.The electron density(a)and temperature(b),and deuterium density(c)and ion temperature(d)profiles in the CFETR steady-state(black) and hybrid (blue) scenarios.The tritium density profile is the same as the deuterium one.

The source of fusion alpha particles is given bySα=nDnT()DT,which determines the birth distribution of alpha particles.Here,nD,Tare the densities of deuterium and tritium,respectively.()DTis the D-T fusion rate,which is given by ()DT=3.68 × 10?18exp (?19.94) m3s?1for T ≤25 keV [34,35].Here,the unit of ion temperature T is keV.Using the density and temperature profiles and the equations above,we can compute the source or birth distribution profiles of alpha particles.Figure 4 shows the radial birth distributions of fusion alpha particles in the CFETR steady-state and hybrid scenarios.The energy of the initial alpha particles is Eα=3.5 MeV.

Figure 4.The radial birth distributions of fusion alpha particles in the CFETR steady-state (black) and hybrid (blue) scenarios.

Figure 5.The contour lines(black)of the toroidal field ripple and the ripple well boundaries (red) in the CFETR steady-state (a) and hybrid (b) scenarios.The black bold lines denote the shape of the first wall and the blue lines denote the LCFS.

Figure 6.The particle and power fractions of lost alpha particles as a function of the ripple amplitude parameter δ0 in the CFETR steady-state and hybrid scenarios with the LCFS and wall boundaries.Here,=1.57 × 10?5.

The ripple parameter δ can be fitted with the following analytic function [11,36]

For the CFETR tokamak,the ripple-related coefficients are chosen as N=16,δ0=1.57×10?5,Rrip=6.01 ?0.062Z2,brip=0.021,wrip=0.63 m [18,19].The contour plots of the toroidal field ripple δ and the ripple well boundaries in the CFETR steady-state and hybrid scenarios are shown in figure 5.The ripple well regions are mainly located on the right side of the ripple well boundaries.The maximal ripples within the LCFS in the steady-state and hybrid scenarios are both about 0.4%.Both of the maximal ripple values are located at R=9.4 m and Z=1 m,that is,the top right of the LCFS.The ripple well boundary is determined bywhere ?is the inverse aspect ratio and θ is the geometric poloidal angle.

In all our alpha particle loss simulations,the Coulomb collision effect was included.The total simulation time is about 1.3 s,which is sufficient for an alpha particle slowing down and becoming a thermal ion.About 20 000 test particles were used for simulating loss fractions,loss regions and heat loads of alpha particles.

4.Numerical results

The numerical results of the alpha particle loss in the CFETR tokamak with ripple and collision computed by the GYCAVA code are presented in this section.

Figure 6 shows the relation between the particle and power losses of alpha particles and the ripple amplitude parameter δ0in the CFETR steady-state and hybrid scenarios with the LCFS and wall boundaries.δ0=0 denotes the case that the ripple effect is not included.In the case without the ripple effect,the alpha particle loss is very small.It is because most of alpha particles are deposited in the core region and in this case the alpha particle loss is only due to the magnetic drift and collision.In the cases with the ripple effect,that is,δ0≠0,we can see from figure 6 that the ripple effect enhance the alpha particle loss.The loss fraction of alpha particles in the steady-state scenario is much larger than that in the hybrid scenario forδ0≥,which is due to the fact that the safety factor in the steady-state scenario is larger than that in the hybrid scenario (see figure 2).According to Goldston et al's theory[31],the ripple stochastic diffusion is positively related to the safety factor.The loss fractions of alpha particles increase rapidly with δ0increasing,especially for the steady-state scenario.The main ripple loss mechanisms are the ripple stochastic diffusion and the ripple well loss.In contrast to the cases with the LCFS boundary,the alpha particle loss fractions are smaller in the cases with the wall boundary,because the wall is outside the LCFS.The alpha particle loss fractions in the steady-state scenario are much larger than those in the hybrid scenario,especially for the large δ0.In the case ofδ0=≡ 1.57 × 10?5,the particle loss fractions of alpha particles loss are small (about 2%) in the two scenarios with the LCFS boundary.However,in the case ofδ0==3.93 × 10?5,the particle loss fraction is 17.9% in the steady-state scenario with the LCFS boundary,which is very closed to the previous numerical results of alpha ripple loss[19]computed by the ORBIT code.Their numerical results indicated that the particle loss fraction is 17% in the case of δ0=4×10?5(see table 2 of reference[19]).Figure 7 shows birth distributions and final positions of lost alpha particles for different ripple amplitude parameters δ0in the CFETR steady-state scenario with the LCFS boundary.We can see that the loss regions(initial positions of lost ions)of alpha particles are mainly at the edge of plasma.With δ0increasing,the loss region of alpha particles becomes wider.Because the loss of alpha particles for δ0=0 is very small,that is,the first orbit loss is small,the loss of alpha particles is closely related to the ripple and Coulomb collision effects for δ0≠0.Forδ0=,some alpha particles in the results on birth distributions and final positions of lost alpha particles are similar to those in the steady-state scenario.

Figure 7.Contours of birth distribution and final positions (blue) of lost alpha particles for different ripple amplitude parameters δ0 in the CFETR steady-state scenario with the LCFS boundary.

Figure 8.Contours of birth distribution and final positions (blue) of lost alpha particles for different ripple amplitude parameters δ0 in the CFETR steady-state scenario with the wall boundary.

Figures 9 and 10 show the 2D distributions of heat loads induced by lost alpha particles for different ripple amplitude parameters δ0in the CFETR steady-state scenario with the LCFS and wall boundaries,respectively.In the case with the LCFS boundary and without the ripple effect shown in figure 9(a),the heat loads of lost alpha particles are located on the LCFS above the mid-plane.However,in the case with the wall boundary and without the ripple effect shown in figure 10(a),heat loads of lost alpha particles are mainly located on the lower divertor.It is related to the fact that plasmas or the LCFS are more close to the lower divertor in contrast to the upper divertor.In the case with the LCFS boundary andδ0=δc0fetrshown in figure 9(b),most of heat loads of lost alpha particles are near the mid-plane,which is mainly induced by the ripple stochastic loss.In the case with plasma core are lost due to the ripple stochastic loss.The final positions of lost alpha particles are mainly located at the LCFS above the mid-plane.Figure 8 shows birth distributions and final positions of lost alpha particles in the steady-state scenario with the wall boundary.We can see that the final positions of lost alpha particles are on the wall.The loss regions of alpha particles are similar to those with the LCFS boundary.

Figure 9.Heat loads induced by lost alpha particles for different ripple amplitude parameters δ0 in the CFETR steady-state scenario with the LCFS boundary.

Alpha particle loss simulations in the hybrid scenario with different boundaries have also been performed.The larger δ0shown in figures 9(c) and (d),we can see some localized heat loads at θ=70°-80°.The number of heat load spots is 16,which is the same at the toroidal field coil number.In addition,the poloidal position of heat load agrees with the upward direction of the magnetic drift of an alpha particle.Therefore,these heat loads are induced by the ripple well loss.

In the case with the wall boundary and0δ=shown in figure 10(b),we also can see 16 heat load spots at the position that the poloidal angle is θ ?70°,which are different from the heat loads in the case with the LCFS boundary shown in figure 9(b).It means that the ripple well loss occurs outside the LCFS and makes the heat loads become localized.Comparing the heat loads in figure 10(b) with those in figure 9(b),we can obtain the result that losses of alpha particles are due to the combination effect of the ripple stochastic diffusion and the ripple well loss in the case with the wall boundary.The other evidence of this result is given by a typical orbit of a lost alpha particle shown in figure 12(c)below.In the cases with the wall boundary,the maximum value of the heat loads is about 0.035 MW m?2forδ0=,but it is about 1 MW m?2forδ0=.

Figure 11 shows the heat loads of lost alpha particles in the hybrid scenario with the wall boundary.The poloidal and toroidal distributions of heat loads shown in figure 11 are similar to those in the steady-state scenario shown in figure 10.The difference between the results of these two scenarios is that the heat loads in the hybrid scenario are much smaller than those in the steady-state scenario.

Figure 10.Heat loads induced by lost alpha particles for different ripple amplitude parameters δ0 in the CFETR steady-state scenario with the wall boundary.

Figure 12.Guiding-center orbits(red)of trapped alpha particles in the CFETR steady-state scenario and in the four cases:(a)the case without collision and ripple,(b) the case without collision and in the presence of ripple δ0=,(c) the case with collision and ripple δ0=,and (d) the case with collision and ripple δ0=.The plus symbol denotes the position of the magnetic axis and the cross symbol denotes the initial position of an alpha particle.

Figure 13.Evolutions of ΔE/E0 in the four cases which are the same as figure 12.Here,ΔE=E ?E0.

In order to well understand the 2D heat load distributions of lost alpha particles,typical orbits of lost alpha particles are studied and discussed below.Guiding-center orbits of trapped alpha particles in the CFETR steady-state scenario in the four cases are shown in figure 12.Four cases include the case without collision and ripple,the case without collision and in the presence of ripple(δ0=),the case with collision and ripple (δ0=),and the case with collision and ripple(δ0=).The corresponding ΔE/E0(ΔE=E ?E0)evolutions of these alpha particles are shown in figure 13.These energetic alpha particles have the same initial phase space position.The initial energy of alpha particles is E0=3.5 MeV and the initial pitch is=? 0.2.In the case without collision and ripple shown in figure 12(a),we can see that the guiding-center orbit of alpha particle is an unperturbed banana orbit.In the case without collision and in the presence of ripple shown in figure 12(b),we can find that the ripple stochastic effect makes the alpha particle move radially and hit the first wall finally.The numerical error of the particle energy of alpha particle in the cases without collision,shown in figures 13(a) and (b),is on the order of 10?4.The simulation time shown in figure 13(a)is 1.3 s,that is,5.4×104poloidal cycles for the trapped alpha particle.It means that the numerical energy conservation is accurate in our long-time simulation of the alpha particle orbit without the collision effect.

In the case with the collision and ripple shown in figure 12(c),it is found that the alpha particles move radially due to the ripple stochastic diffusion effect in the early stage and then move vertically upward due to the ripple well loss mechanism outside the LCFS.In the other words,the alpha particle loss is due to the ripple stochastic loss when the LCFS boundary is used,whereas the alpha particle loss is due to the ripple stochastic transport and the ripple well loss when the wall boundary is used.The radial position of the ripple well loss agrees with the analytic ripple well region shown in figure 5(a).The orbit shown in figure 12(c) is a typical orbit of a lost alpha particle in the CFETR magnetic configuration.The typical orbit can explain the heat loads of lost alpha particles with different outer boundaries shown in figures 9(b)and 10(b).By comparing the cases without and with collision shown in figures 12(b)and(c),we can find that the Coulomb collision effect plays an important role in the transition of the lost alpha particle behavior from the ripple stochastic diffusion to the ripple well trapping.In the cases with collision,the alpha particle energy decreases rapidly with time due to the collision-induced slowing down effect,which can be seen from figures 13(c) and (d).

As shown in figure 12(d),the alpha particle orbit in the case withδ0=has the similar behavior of that in the case with 0 0cfetrδ=δshown in figure 12(c).The transition from the ripple stochastic transport to the ripple well trapping occurs within the LCFS in contrast to the orbit shown in figure 12(c).It is because the larger ripple amplitude parameter δ0has been used in the former case.These typical orbits can also explain the heat loads of lost alpha particles shown in figures 9(c) and 10(c).In addition,we can see from figures 13(c) and (d) that the alpha particle loss occurs more rapidly in the case withδ0=.

5.Summary

Energetic alpha particle losses with the toroidal field ripple and the Coulomb collision have been numerically studied by the orbit-following code GYCAVA in the CFETR steadystate and hybrid scenarios.The effects of the outer boundary and the ripple amplitude on alpha particle losses have been investigated.The LCFS or the first wall has been used as the loss boundary for simulating loss fractions,loss regions and heat loads of alpha particles.

The loss fractions and heat loads of alpha particles in the hybrid scenario are much smaller than those in the steady-state scenario for a significant ripple amplitude δ0.The loss fraction of alpha particles increases rapidly and the loss region becomes wider with δ0increasing.The loss of alpha particles is closely related to the ripple and Coulomb collision effects.The loss regions of alpha particles are mainly at the edge of plasma.Some alpha particles in the plasma core are lost due to the ripple stochastic transport for a large ripple amplitude δ0.

For the CFETR tokamak,heat loads with the LCFS boundary are mainly near the mid-plane,which is due to the ripple stochastic loss.However,heat loads with the wall boundary are localized and mainly located at the position that the poloidal angle is θ=70°-80°.This position of heat load is related to the direction of the magnetic drift.Discrete heat load spots are due to the ripple well loss.

The heat loads of lost alpha particles with different outer boundaries can be explained by typical alpha particle orbits.The transition of the lost alpha particle behavior from the ripple stochastic diffusion to the ripple well trapping has been identified in the CFETR tokamak.The transition can occur inside or outside the LCFS,which is dependent of the ripple amplitude parameter δ0.The Coulomb collision effect is responsible for the transition from the ripple stochastic diffusion to the ripple well trapping.

Acknowledgments

The authors appreciate the support from the CFETR team.Numerical simulations were performed on the ShenMa High Performance Computing Cluster in Institute of Plasma Physics,Chinese Academy of Sciences.The work was jointly supported by National Natural Science Foundation of China(Nos.12175034,12005063),the National Key R&D Program of China (No.2019YFE03030001) and the Fundamental Research Funds for the Central Universities (No.2232022G-10).

ORCID iDs

主站蜘蛛池模板: 日本一区二区三区精品视频| 欧美日韩高清| 欧美三级自拍| 一级片一区| 青青青草国产| 中文字幕天无码久久精品视频免费 | 97视频在线观看免费视频| 精品人妻无码中字系列| 国产人成乱码视频免费观看| 久久性妇女精品免费| 欧美视频在线播放观看免费福利资源 | 综1合AV在线播放| 秋霞国产在线| 亚洲区一区| 中日韩一区二区三区中文免费视频| 久久国产精品影院| 麻豆精品久久久久久久99蜜桃| 91网红精品在线观看| 国产成人AV综合久久| 99这里精品| 在线观看亚洲人成网站| 亚洲成年人片| 最新国产网站| 亚洲丝袜中文字幕| 91麻豆国产视频| 伊人国产无码高清视频| 国产99在线观看| 一级毛片a女人刺激视频免费| 在线亚洲小视频| 97精品国产高清久久久久蜜芽| 国产三级国产精品国产普男人 | www.日韩三级| 成年人视频一区二区| 天堂成人在线| 国产精品黄色片| 伊人久久久久久久久久| 国产AV无码专区亚洲A∨毛片| 久久中文无码精品| 亚洲精品第五页| 狠狠色综合久久狠狠色综合| 美女视频黄又黄又免费高清| 一级在线毛片| 免费看的一级毛片| 亚洲永久精品ww47国产| 欧美日韩在线亚洲国产人| 无码福利视频| 国产aⅴ无码专区亚洲av综合网| 免费A级毛片无码无遮挡| 狠狠ⅴ日韩v欧美v天堂| 日韩大片免费观看视频播放| 高清大学生毛片一级| 国产成人亚洲无码淙合青草| 99999久久久久久亚洲| 亚洲精品天堂在线观看| 在线播放国产一区| 久久久久久久久18禁秘| 好紧太爽了视频免费无码| 又爽又黄又无遮挡网站| 久久久噜噜噜| 国产欧美精品一区aⅴ影院| 美女扒开下面流白浆在线试听| 日韩午夜片| 免费国产高清视频| 国产中文一区二区苍井空| 欧美一级视频免费| 色综合激情网| 成人毛片在线播放| av大片在线无码免费| 青青草国产在线视频| 98精品全国免费观看视频| 精品免费在线视频| 动漫精品中文字幕无码| 色香蕉影院| 国产小视频在线高清播放 | 国产激情第一页| 国产日本一区二区三区| 亚洲天堂免费观看| 日韩AV无码免费一二三区| 欧美中文字幕一区二区三区| 国模沟沟一区二区三区| 四虎国产在线观看| 亚洲精品亚洲人成在线|