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

    Simulations of NBI fast ion loss in the presence of toroidal field ripple on EAST

    2021-09-10 09:26:10YingfengXU徐穎峰YoujunHU胡友俊XiaodongZHANG張曉東XingyuanXU徐行遠(yuǎn)LeiYE葉磊XiaotaoXIAO肖小濤andZhenZHENG鄭振
    Plasma Science and Technology 2021年9期
    關(guān)鍵詞:張曉東徐行

    Yingfeng XU (徐穎峰),Youjun HU (胡友俊),Xiaodong ZHANG (張曉東),Xingyuan XU (徐行遠(yuǎn)),4,Lei YE (葉磊),Xiaotao XIAO (肖小濤) and Zhen ZHENG (鄭振),4

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

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

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

    4 University of Science and Technology of China,Hefei 230026,People's Republic of China

    Abstract NBI fast ion losses in the presence of the toroidal field ripple on EAST have been investigated by using the orbit code GYCAVA and the NBI code TGCO.The ripple effect was included in the upgraded version of the GYCAVA code.It is found that loss regions of NBI fast ions are mainly on the low field side near the edge in the presence of ripple.For co-current NBIs,the synergy effect of ripple and Coulomb collision on fast ion losses is dominant,and fast trapped ions located on the low field side are easily lost.The ripple well loss and the ripple stochastic loss of fast ions have been identified from the heat loads of co-current NBI fast ions.The ripple stochastic loss and the collisioninduced loss are much larger than the ripple well loss.Heat loads of lost fast ions are mainly localized on the right side of the radio frequency wave antennas from the inside view toward the first wall.For counter-current NBIs,the first orbit loss due to the magnetic drift is the dominant loss channel.In addition,fast ion loss fraction with ripple and collision for each NBI linearly increases with the effective charge number,which is related to the pitch angle scattering effect.

    Keywords: fast ion loss,ripple,NBI,heat load,collision

    1.Introduction

    Fast particles play an important role in tokamak fusion plasmas.Fusion reaction can produce energetic alpha ions.Fast particles can also be produced by radio frequency wave(RFW) heating and neutral beam injection (NBI).Fast ion losses not only reduce the heating efficiency,but also destroy the first wall of tokamaks and produce impurities.

    The discreteness of toroidal field coils destroys the axisymmetry of the tokamak.This toroidal asymmetric perturbation is called the toroidal field ripple.The toroidal field ripple can induce energetic or fast ion losses and may cause the localized heat load on the first wall,especially RFW antennas.The resulting heat load on the first wall limits the toroidal field ripple amplitude and the number of toroidal field coils in a tokamak fusion reactor.Simulations of energetic ion losses with the toroidal field ripple have been performed by using various orbit codes,such as OFMC [1,2],ORBIT [3],and VENUS-LEVIS[4].Ripple-induced losses of NBI fast ions have been studied experimentally and numerically in many tokamaks,such as JT-60U [5–7],JET [8,9],TFTR [10,11],and EAST [12,13].Previous simulations of the ripple loss on EAST were performed by using the orbit code ORBIT.The canonical variables used by ORBIT for computing the guiding-center orbits are based on the magnetic flux coordinates,in which the safety factor is singular at the X points.Therefore,the last closed flux surface (LCFS),instead of the first wall,was used as the outer/loss boundary in the previous ripple loss simulations on EAST.In addition,the 2D distributions of heat loads and loss regions have not been given in the previous ripple-induced fast ion loss simulations on EAST.

    We have recently developed a Monte Carlo orbit-following code GYCAVA [14,15],which can be used for investigating the behavior of test particles in a tokamak.It is of interest to numerically investigate NBI ion losses and their 2D heat loads in the presence of ripple on EAST,with the first wall being the loss boundary by the GYCAVA code.

    In this work,numerical simulations of NBI ion losses on EAST in the presence of the toroidal field ripple and collision by the orbit code GYCAVA and the NBI code TGCO are presented.The TGCO code has been benchmarked with the NBI code NUBEAM.The GYCAVA code has been integrated with TGCO[15].The toroidal field ripple effect has been added in the new version of the code GYCAVA.The cylindrical coordinates were used in our GYCAVA simulations for avoiding the problem of the singularity of the safety factor at the X points.The first wall without and with RFW antennas has been used as the outer boundary in our orbit simulations.The ripple and collision effects on fast ion behaviors of NBI ions on EAST have been studied.The TGCO code is used for computing the initial distribution of NBI ions,which is an input of the orbit-following code GYCAVA.Previous work on the ripple loss of NBI ions on JT-60U [5] has pointed out that the radial electric field has a small effect on the heat load position induced by ripple-trapped NBI ions.Therefore,the effect of the radial electric fieldErhas not been included in our simulations for simplicity.The radial electric field effect on fast ion losses will be investigated in future.

    The remaining part of this paper is organized as follows.The equations of motion in the presence of toroidal field ripple and collision used in the GYCAVA code are presented in section 2.In section 3,the simulation setup,including the equilibrium of EAST,the plasma profiles,the birth distributions of NBI ions and the toroidal field ripple,is given.The numerical results of ripple- and collision-induced NBI ion losses on EAST are presented and discussed in section 4.Section 5 summarizes the main results.

    2.Equations of motion in the presence of ripple and collision used in GYCAVA

    For investigating the ripple effect on NBI ion losses,we have added the ripple effect in the new version of the GYCAVA code.The guiding-center equations of motion including the toroidal field ripple and the Coulomb collision used in GYCAVA are introduced as follows.

    The axisymmetric equilibrium magnetic field for a tokamak can be expressed as

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

    The toroidal field ripple can be defined as

    Here,BTminandBTmaxare the minimum and maximum values of the toroidal magnetic field which are computed at two points with the same radial and vertical positions,but different toroidal positions.Therefore,the magnitude of ripple δ is a function ofRandZ.

    The ripple effect on fast ion losses,including the ripple well loss and the ripple stochastic loss,is mainly due to the toroidal component of the ripple field [16,17].Therefore,only the toroidal ripple field is kept in our GYCAVA code for simplicity.The magnetic field including toroidal field ripple isThe toroidal field perturbation related to the ripplecan be written as

    We can see thatδBripis minimal whenfor the fixedRandZ.Here,Nis the number of the toroidal field coils.Then the total magnetic field can be written as

    The guiding-center equations of motion including the toroidal field ripple can be written as

    Here,Xis the guiding-center position of a charged particle andv‖is the parallel velocity.msandesare the mass and charge of the charged particles,respectively.B0.The guiding-center Hamiltonian is expressed as

    Note thatB?B0+δBrip.

    The guiding-center equations of motion (6),(7) can be explicitly expressed in terms of the Possion matrix and the cylindrical coordinates as

    The guiding-center orbits were computed in GYCAVA simulations.However,the finite Larmor radius effect can impact the loss fractions [18] and the hit points [19] of fast ions.Therefore,the particle coordinates,which are obtained by the guiding-center transform from the guiding-center coordinates,are used for judging whether a fast ion is lost or not in GYCAVA.

    The Coulomb collision of fast ions and background particles including the slowing down and the pitch angle scattering is contained in the GYCAVA code [15,20].The slowing down of fast ions due to electron and ion drag can be expressed as

    Here,νsis the slowing down collision rate,written as

    neis the electron density.eandefare the electron and fast ion electric charges,respectively.Zmis defined as=ZmThe subscriptjis the background ion species.Λln is the Coulomb logarithm.vcis the critical velocity,expressed as

    Here,me,mfare the electron and fast ion masses,respectively.Teis the electron temperature.

    The pitch angle scattering of fast ions is computed by a Monte Carlo method [21].The new pitch of particle,λnew,is computed from the old one,λold,by

    Here,λ=v‖/vand Δtis the time step used in orbit simulations.νdis the pitch angle scattering rate,expressed as

    Here,Zeffis the effective charge number,defined asIn our simulations,the default effective charge number is chosen asZeff=2.The NBI fast ions and bulk ions on EAST are both deuterium.Therefore,Zm?1.

    3.Simulation setup

    In this section,the magnetic configuration of EAST,plasma profiles,the birth distributions of fast ions produced by different NBIs,and the ripple field on EAST are presented.

    Figure 1.The magnetic configuration (magenta line) and the shape of the first wall(black line)and the LCFS(blue line)on EAST.The red line denotes the shape of the RFW antennas.

    Figure 2.The safety factor profile for EAST discharge #62585.

    The parameters of EAST discharge#62585 were used in our NBI ion loss simulations.The magnetic equilibrium configuration of this discharge and the first wall with the antennas of RFW are shown in figure 1.The direction of the plasma toroidal current is counter-clockwise from the top view,and the direction of the magnetic field is clockwise,which means that the magnetic drift direction of ions is downward.The major radius at the axis isRaxis=1.91 m and the minor radius isa=0.46 m.The magnetic field at the axis isBaxis=2.2 T and the plasma current isIp=0.56 MA.Figure 2 shows the safety factor profile.ρtdenotes the square root of the normalized toroidal magnetic flux.The range of the safety factor is 1.6–7.6 with the value at 95% of the poloidal magnetic flux isq95=4.8.The density and temperature profiles were used as inputs of the TGCO code for simulating the deposition of NBIs on EAST.Figure 3 shows the plasma profiles,including the electron density,the temperatures of electron and ion for EAST discharge#62585.The maximum electron density isne0=6.8 ×1019m?3.The maximum electron and ion temperatures areTe0=2.2 keV andTi0=2.3 keV,respectively.

    Figure 3.The electron density (a),electron temperature (b) and ion temperature (c) profiles for EAST discharge #62585.

    Figure 4.The initial spatial distributions of different NBI ions from the top view: (a) co-perp NBI,(b) co-tang NBI,(c) ctr-tang NBI,(d) ctr-perp NBI,and (e) the distribution for the co-perp NBI in the cross-section.

    The parameters of NBIs used in our simulations were chosen as those of NBIs [13] on EAST.Four NBIs with different tangential radii were used in most simulations.Two NBIs are injected in the co-current direction,while the other two NBIs are in the counter-current direction.The tangent radii of the co-current tangential and perpendicular NBIs arerespectively.The tangent radii of the counter-current tangential and perpendicular NBIs are0.606 m,respectively.The full energy of beams are chosen asEb=50 keV in all cases of our simulations.The species of NBI ions is deuterium on EAST.NBI deposition on EAST is simulated by the NBI code TGCO.This code is used for computing the birth distributions of fast ions produced by four NBIs on EAST.Figure 4 shows the birth distributions of fast ions produced by NBIs.Figures 4(a)–(d) show the birth spatial distributions for four NBIs from the top view,respectively.Figure 4(e) shows the birth spatial distributions for the co-current perpendicular NBI in the poloidal crosssection.In fact,the spatial distributions for different NBIs in the poloidal cross-section are similar.

    The ripple parameter δ can be fitted with an analytic function [13,22,23],expressed as

    For the EAST tokamak,the coefficients related to the toroidal field ripple are chosen asN=16,δ0=1.267 ×10?4,Rrip=1.714 ?0.181Z2,brip=0.267,wrip=0.149 m [12,13].Figure 5 shows the contour lines of the toroidal field ripple δ.The maximum value of the ripple δ within the LCFS is 0.76%.The position of maximal δ is located atR=2.32 m andZ=0.04 m,that is,the rightmost of the LCFS.

    The poloidal cross-section of plasmas can be divided by the ripple well boundary[16]into two regions,the ripple well region and the reminder.The ripple well boundary is given by

    Here,?is the inverse aspect ratio,that is,?=r/R.ris the minor radius.θ is the geometric poloidal angle.The ripple well region is the plasma region with α<1,in which toroidal wells exist.The ripple well boundary for EAST is shown in figure 5.For the regions enclosed by the ripple well boundary,α>1 and toroidal wells do not exist.The plasma region outside the ripple well boundary is the ripple well region for EAST.In the ripple well region,fast ions with very small pitch can be trapped in the toroidal wells and then they are lost through the magnetic drift in the vertical direction.The loss positions of fast ions due to the ripple well loss are vertically above or below the ripple well region,which depends on the direction of the magnetic drift or the toroidal magnetic field.

    4.Numerical results

    Orbits of two trapped 50 keV fast ions in the EAST tokamak in the presence of the toroidal field ripple are shown in figure 6.From figure 6(a),we can see that the fast ion with the initial pitchv‖/v=0.2 finally hits the first wall near the mid-plane.This fast ion loss is due to the ripple stochastic loss[17].Ripple makes turning points of the fast trapped ion become stochastic.From figure 6(b),we can see that the fast ion with a very small pitch (v‖/v=0.05) finally hits the first wall vertically below the initial position of the fast ion.This fast ion loss is due to the ripple well loss,which is induced by the magnetic drift and the ripple well.The magnetic drift is downward in this case,and therefore the lost position is below the mid-plane.

    In figures 7–10,fast ion loss simulations were performed in four cases: (i)without ripple and collision,(ii)with ripple,(iii) with collision and (iv) with collision and ripple.In these simulations,the 3D wall (the wall with RFW antennas) has been used as the outer boundary.The 3D wall used in GYCAVA includes two LHW antennas [24] and two ICRF antennas[25]on EAST.After the upgrade of EAST in 2020,LHW antennas are located at ports B and E of EAST,respectively.Two ICRF antennas are located at ports N and I,respectively.The parameter setup of RFW antennas in GYCAVA roughly agrees with that on EAST.

    Figure 5.The contour plot(black)of the toroidal field ripple and the boundaries of the ripple well(red)on the EAST tokamak,as well as the shapes of the LCFS (blue) and the wall (black).

    Figure 6.Orbits (red line) of two trapped 50 keV fast ions in the EAST tokamak in the presence of the toroidal field ripple with the initial pitch v‖/v=0.2 (a) and v‖/v=0.05 (b).The initial radial position is ?ψ = 0.8.Here,?ψ is the normalized poloidal magnetic flux.The cross and plus symbols denote the initial position of a fast ion and the magnetic axis position,respectively.

    The loss fraction at the final time computed from our simulations,which almost reaches a steady state,can represent the loss fraction of fast ions for continuous NBI,if the NBI ion source,full magnetic fields and plasma profiles remain unchanged.The heat loads computed from our simulations also represent the heat loads for continuous NBI.

    Figure 7.Time histories of particle (a) and power (b) loss of co-current perpendicular NBI fast ions in the four cases: without ripple and collision (cyan),with ripple (red),with collision (blue),with collision and ripple (black).

    Figure 8.Time histories of particle(a)and power(b)loss of co-current tangential NBI fast ions in the four cases:without ripple and collision(cyan),with ripple (red),with collision (blue),with collision and ripple (black).

    Time histories of particle and power loss fraction of cocurrent perpendicular(co-perp)and co-current tangential(cotang)NBI fast ions are shown in figures 7 and 8,respectively.We can see that the first orbit losses are small for co-perp and co-tang NBIs.In the cases with collision,fast ion loss increases with time at the early stage and reaches a steady state finally.The pitch angle scattering is responsible in collision-induced loss,especially for the co-perp NBI.In contrast to the co-tang NBI,the ripple-induced loss for the co-perp NBI is more significant.This is because trapped fast ions tend to lose due to ripple in contrast to passing fast ions,and the co-perp NBI produces more trapped ions than the co-tang NBI.The power loss with collision is smaller than the particle loss,which is due to the slowing down effect.

    Time histories of particle and power loss fraction of counter-current tangential (ctr-tang) and counter-current perpendicular (ctr-perp) NBI fast ions are shown in figures 9 and 10,respectively.We find that the first orbit loss is very large for counter-current NBIs.The first orbit loss is 41%for ctr-tang NBI and 31% for ctr-perp NBI,respectively.This is because most NBI fast ions are deposited in the low field side in contrast to the high field side.The ripple-induced loss for the ctr-perp NBI is more significant than the ctr-tang NBI,which is similar to the result of co-current NBIs.In the cases with collision and ripple,fast ion loss for the ctr-tang NBI is similar to that for the ctr-perp NBI.The particle loss fraction is about 55% and the power loss fraction is about 50%.

    Figure 9.Time histories of particle(a)and power(b)loss of counter-current tangential NBI ions in the four cases:without ripple and collision(cyan),with ripple (red),with collision (blue),with collision and ripple (black).

    Figure 10.Time histories of particle (a) and power (b) loss of counter-current perpendicular NBI ions in the four cases: without ripple and collision (cyan),with ripple (red),with collision (blue),with collision and ripple (black).

    We have also performed fast ion loss simulations with the 2D wall boundary.It is found that the fast ion loss with the 3D wall boundary is slightly larger than that with the 2D wall boundary.This is because the RFW antennas are included in the 3D wall boundary.

    Birth distributions and striking positions of lost beam ions for different NBIs are shown in figures 11 and 12.In figure 11,ripple and collision have not been included and the 2D wall has been used as the outer boundary.Fast ion losses are only due to the first orbit loss.We can see that the loss region,that is,the birth distribution of lost beam ions,is determined by the NBI directions.Loss regions are on the high field side near the edge for co-current NBIs,while loss regions of beam ions are on the low field side for countercurrent NBIs.In figure 12,ripple and collision have been included and the 3D wall has been used as the outer boundary.We can find that loss regions for four different NBIs are mainly on the low field side.The differences of the loss regions shown in figures 11 and 12 are due to the fact that fast ion losses are enhanced by the synergy effect of ripple and collision.Ripple makes trapped fast ions easily lose and trapped fast ions are mainly on the low field side.Collisioninduced pitch angle scattering can enhance the ripple loss.The loss region with the 2D wall boundary is similar to that with the 3D wall boundary.However,final positions of lost fast ions are significantly affected by the RFW antennas.

    The heat loads induced by lost NBI ions with ripple and collision for different NBIs are shown in figures 13–14.In figure 13,the 2D wall boundary has been used and the 3D wall boundary has been used in figure 14.

    From figure 13,we can find that heat loads induced by lost NBI ions are on the first wall below the mid-plane and the lower divertor.The positions of these heat loads are affected by the magnetic field direction.The toroidal magnetic field used in our simulations is in the clockwise direction from the top view.This means that the magnetic drift direction of a fast ion is downward and thus lost fast ions hit the wall below the mid-plane.Heat loads for counter-current NBIs are much more significant than those for co-current NBIs.This is because the first orbit loss for counter-current NBIs is much larger,which is due to the fact that fast ions born on the low field side are much larger than those on the high field side.

    From figures 13(a) and (b),we can see that most of heat loads are on the first wall near the mid-plane for co-current NBIs.These heat loads are mainly due to the synergy effect of ripple and collision.From figure 6(a),we know that the locations of these heat loads are the same as the striking points of lost trapped fast ions due to the ripple stochastic loss.In addition,the pitch angle scattering effect induced by collision can enhance the ripple stochastic loss significantly.Note that the effect of the first orbit loss on the heat load for co-current NBIs is small in contrast to the synergy effect of ripple and collision,which is shown in figure 7 and 8.Near θ=?60°(θ is the poloidal angle);we also can see that some of localized heat loads in the toroidal direction.From figure 6(b),we know that these heat loads are mainly due to the ripple well loss.The number of the localized heat loads near θ=?60°is the same as the number of the toroidal field coils,that is,N=16.Ripple well-induced losses are on the first wall vertically below the ripple well region.Fast ions with very small pitch,which are produced by collisioninduced pitch angle scattering,can be trapped in the toroidal wells.Then they are lost through the magnetic drift in the vertical direction.The fact that the loss position is below the Figure 11.Contours of birth distribution and striking positions(blue)of lost NBI ions in the absence of ripple and collision with the 2D wall boundary.Figure 12.Contours of birth distribution and striking positions(blue)of lost NBI ions in the presence of ripple and collision with the 3D wall boundary.mid-plane is also related to the fact that the magnetic drift is downward in our simulations.Ripple well loss is a rapid process.The timescale of the ripple well loss isa/vdwithvdthe magnetic drift.In contrast to the ripple well loss,the ripple stochastic loss and the collision-induced loss are dominant.

    Figure 13.Heat loads induced by lost fast ions divided by the NBI injected power in the presence of ripple and collision with the 2D wall boundary.

    Figure 14.Heat loads induced by lost fast ions divided by the NBI injected power in the presence of ripple and collision with the 3D wall boundary.

    Figure 15.The particle and power loss fractions of fast ions as a function of the effective charge number for different NBIs in the presence of ripple and collision.

    For counter-current NBIs,we find that heat loads induced by lost fast ions produced by tangential NBI are on the first wall below the mid-plane.For the tangential NBI,significant heat loads are on the lower divertor,which is shown in figure 13(c).Meanwhile heat loads of fast ions produced by the perpendicular NBI are mainly on the first wall near the mid-plane,which is shown in figure 13(d).These heat loads are due to the first orbit loss and the loss induced by ripple and collision.For counter-current NBIs,the effect of the first orbit loss on heat loads is dominant in contrast to the ripple and collision effects.

    From figure 14,we can see that heat loads induced by lost fast ions near the mid-plane become more localized in contrast to figure 13,which is caused by the presence of the antennas of RFW.Heat loads are mainly located on the right side of all antennas of RFW from the inside view toward the first wall,especially on their edge.For co-current NBIs,heat loads are on the right side of all RFW antennas.Heat loads are mainly on the lower right corner of the ICRH2 antenna for ctr-tang NBI,while they are mainly on the lower right corner of the LHCD2 antenna for ctr-perp NBI.

    The original counter-current NBIs on EAST will become the co-current NBIs after the upgrade of the EAST tokamak.The tangent radii of the new co-current tangential and perpendicular NBIs are0.606 m,respectively.Figure 15 shows the relation between the particle and power loss fractions of fast ions and the effective charge number for different NBIs in the presence of ripple and collision.The 3D wall boundary was used as the outer boundary.The simulation time for each case is about 50 ms so that the loss fraction can reach a steady state.It is found from figure 15 that the loss fractions of particle or power for each NBI linearly increase with the effective charge number.Note that the slowing down effect and the pitch angle scattering effect are included in the Coulomb collision.The slowing down effect is independent of the effective charge number,but the pitch angle scattering effect is proportional to the effective charge number.Therefore,the fast ion loss fractions induced by the effect of the pitch angle scattering increase with the effective charge number increasing.In contrast to the counter-current NBIs,the fast ion loss fractions of co-current NBIs are much smaller,especially for the cocurrent tangential NBIs.

    5.Summary

    NBI fast ion losses in the presence of the ripple field on the EAST tokamak have been simulated by the orbit-following code GYCAVA and the NBI code TGCO.The ripple effect was included in the new version of the orbit code GYCAVA.The walls without and with RFW antennas have been used as the outer or loss boundary in our orbit-following simulations.

    The effects of the magnetic drift,ripple and Coulomb collision on fast ion losses of NBI have been investigated.Loss regions of fast ions are mainly on the low field side near the edge in the presence of ripple for all NBIs.For co-current NBIs,fast trapped ions located on the low field side near the edge are easily lost due to ripple and collision,and the NBI fast ion loss induced by the synergy effect of ripple and collision is the dominant loss channel.For counter-current NBIs,the first orbit loss is the dominant loss channel,and fast ions on the low field side near the edge tend to lose.

    The ripple well loss and the ripple stochastic loss of fast ions have been identified from heat loads of fast ions for cocurrent NBIs.The lost fast ions due to the ripple well loss have very small pitch and they can be produced by the pitch angle scattering.The ripple stochastic loss and the collisioninduced loss are much larger than the ripple well loss.Heat loads induced by lost fast ions are mainly localized on the right side the RFW antennas with the wall including RFW antennas from the inside view toward the first wall.Heat loads for co-current NBIs are near the mid-plane,while heat loads for counter-current NBIs are on the lower right corner of RFW antennas and the lower divertor.

    The loss fraction of particle or power with ripple and collision for each NBI linearly increases with the effective charge number increasing.This is due to the fact that the collision-induced pitch angle scattering effect is proportional to the effective charge number.

    In this work,we focus on simulations of NBI ion losses with ripple.The simulation results will be compared with experimental observations on EAST in future.

    Acknowledgments

    The authors appreciate the support from the EAST team.Numerical simulations were performed on the ShenMa High Performance Computing Cluster in the Institute of Plasma Physics,Chinese Academy of Sciences.The work was supported by National Natural Science Foundation of China(No.11 775 265).

    ORCID iDs

    猜你喜歡
    張曉東徐行
    一剪梅·中秋賞月
    狂奔向你:秒速5厘米
    偷吃一枚小月亮
    花火彩版A(2021年1期)2021-09-10 07:22:44
    張曉東《篆書國有歲以團(tuán)扇》
    張曉東治療原發(fā)性痛經(jīng)經(jīng)驗
    窗花
    小河的漣漪
    我們的愿望
    竹外疏花
    家居廊(2019年10期)2019-09-10 07:22:44
    《花亂開》
    欧美三级亚洲精品| 夜夜夜夜夜久久久久| 1024香蕉在线观看| 免费看a级黄色片| 白带黄色成豆腐渣| 国产在线精品亚洲第一网站| 成年免费大片在线观看| 久99久视频精品免费| 一本久久中文字幕| 色综合站精品国产| 亚洲精品美女久久av网站| 国产91精品成人一区二区三区| 午夜影院日韩av| 日韩中文字幕欧美一区二区| 无限看片的www在线观看| 免费一级毛片在线播放高清视频| 黄色丝袜av网址大全| 亚洲一区高清亚洲精品| 成人av在线播放网站| 午夜福利在线观看吧| 国产高清有码在线观看视频| 国产一区二区在线av高清观看| 国产免费av片在线观看野外av| 一进一出抽搐gif免费好疼| 国产视频一区二区在线看| 欧美精品啪啪一区二区三区| 亚洲成av人片在线播放无| 久久久精品大字幕| 色尼玛亚洲综合影院| 亚洲欧美精品综合久久99| 成人av在线播放网站| 欧美成狂野欧美在线观看| 日韩高清综合在线| 国产伦人伦偷精品视频| 99精品欧美一区二区三区四区| 久久婷婷人人爽人人干人人爱| 18禁黄网站禁片免费观看直播| 国产精品久久久av美女十八| 国产三级中文精品| 日韩精品中文字幕看吧| 欧美日本视频| 两个人看的免费小视频| 精品久久久久久久久久久久久| 午夜a级毛片| 动漫黄色视频在线观看| 亚洲av成人一区二区三| 亚洲欧美日韩高清在线视频| 亚洲av成人精品一区久久| 中出人妻视频一区二区| 亚洲人成电影免费在线| 成人无遮挡网站| 国产精品久久久久久人妻精品电影| 欧美性猛交╳xxx乱大交人| АⅤ资源中文在线天堂| 日本精品一区二区三区蜜桃| 亚洲国产精品999在线| 精品电影一区二区在线| 麻豆久久精品国产亚洲av| 白带黄色成豆腐渣| 久久精品国产综合久久久| 欧美国产日韩亚洲一区| 精品99又大又爽又粗少妇毛片 | 国产成人一区二区三区免费视频网站| 18美女黄网站色大片免费观看| 日日干狠狠操夜夜爽| 噜噜噜噜噜久久久久久91| 88av欧美| 亚洲av美国av| 1024香蕉在线观看| 久久精品国产清高在天天线| 亚洲美女视频黄频| 男女做爰动态图高潮gif福利片| 国产又色又爽无遮挡免费看| 国产不卡一卡二| 午夜福利欧美成人| 91av网一区二区| 精品久久久久久,| 婷婷精品国产亚洲av在线| 亚洲电影在线观看av| 全区人妻精品视频| 日本黄色视频三级网站网址| 欧美日韩综合久久久久久 | 久久久久久国产a免费观看| 香蕉av资源在线| 91麻豆精品激情在线观看国产| 美女午夜性视频免费| 少妇人妻一区二区三区视频| 伊人久久大香线蕉亚洲五| 色哟哟哟哟哟哟| 久久香蕉精品热| 亚洲专区国产一区二区| 老熟妇乱子伦视频在线观看| 麻豆一二三区av精品| 波多野结衣巨乳人妻| 日日夜夜操网爽| 精品久久久久久久久久久久久| 国产久久久一区二区三区| www.精华液| 亚洲第一电影网av| 高清在线国产一区| 亚洲国产精品久久男人天堂| 亚洲精品456在线播放app | 一本一本综合久久| 久久久久国内视频| 亚洲欧美日韩卡通动漫| 香蕉丝袜av| 成年女人永久免费观看视频| 亚洲国产欧洲综合997久久,| 免费电影在线观看免费观看| 亚洲美女视频黄频| 两个人的视频大全免费| 淫妇啪啪啪对白视频| 国产真人三级小视频在线观看| 亚洲精品中文字幕一二三四区| 亚洲国产精品久久男人天堂| 在线免费观看的www视频| 国产高清视频在线播放一区| av在线蜜桃| 1000部很黄的大片| 亚洲成人精品中文字幕电影| 啦啦啦免费观看视频1| 久久天堂一区二区三区四区| 亚洲国产欧洲综合997久久,| 国产三级黄色录像| 日本精品一区二区三区蜜桃| 午夜免费成人在线视频| 亚洲国产精品sss在线观看| 精品国产三级普通话版| 中文字幕人妻丝袜一区二区| 一区二区三区高清视频在线| 欧美日本亚洲视频在线播放| 国产高清激情床上av| 国产精品爽爽va在线观看网站| 最近最新免费中文字幕在线| 夜夜夜夜夜久久久久| 国产欧美日韩精品一区二区| 中文亚洲av片在线观看爽| 日韩精品青青久久久久久| 国产av不卡久久| 午夜激情福利司机影院| 国产精品一区二区三区四区免费观看 | 国产免费男女视频| 宅男免费午夜| 成人av在线播放网站| 国产午夜精品论理片| 最新中文字幕久久久久 | 国产伦一二天堂av在线观看| 国产精品,欧美在线| 成人av在线播放网站| 亚洲专区字幕在线| 91麻豆精品激情在线观看国产| 别揉我奶头~嗯~啊~动态视频| 亚洲国产欧美网| 欧美又色又爽又黄视频| 在线观看美女被高潮喷水网站 | 色精品久久人妻99蜜桃| 高潮久久久久久久久久久不卡| 亚洲精品在线观看二区| 这个男人来自地球电影免费观看| 在线观看舔阴道视频| 国产精品香港三级国产av潘金莲| 欧美乱妇无乱码| 91在线观看av| 日本精品一区二区三区蜜桃| 成人三级做爰电影| 每晚都被弄得嗷嗷叫到高潮| 黄色视频,在线免费观看| 精品99又大又爽又粗少妇毛片 | 国产97色在线日韩免费| 亚洲熟女毛片儿| 成人特级av手机在线观看| 欧美成人免费av一区二区三区| 国产亚洲精品久久久久久毛片| 国产精品九九99| 性欧美人与动物交配| 在线观看66精品国产| 99在线人妻在线中文字幕| xxxwww97欧美| av在线蜜桃| 日本免费a在线| 国产高清videossex| 久久久久久人人人人人| 成人午夜高清在线视频| 可以在线观看毛片的网站| 亚洲国产精品999在线| 午夜免费激情av| 欧美高清成人免费视频www| 久久精品aⅴ一区二区三区四区| 精品一区二区三区视频在线 | www日本黄色视频网| 给我免费播放毛片高清在线观看| 老司机深夜福利视频在线观看| 国产免费av片在线观看野外av| 亚洲成人久久爱视频| 少妇熟女aⅴ在线视频| 欧美黑人欧美精品刺激| 99久久综合精品五月天人人| www.自偷自拍.com| 青草久久国产| 色老头精品视频在线观看| 亚洲精品色激情综合| 欧美黄色淫秽网站| 在线观看免费视频日本深夜| 亚洲精品中文字幕一二三四区| 国产高清三级在线| 国产男靠女视频免费网站| 日本 欧美在线| 国产伦精品一区二区三区四那| 亚洲国产看品久久| 国产视频内射| 悠悠久久av| 一夜夜www| 天堂网av新在线| 亚洲美女黄片视频| av女优亚洲男人天堂 | 18美女黄网站色大片免费观看| 中文字幕高清在线视频| 国产高清激情床上av| 国产麻豆成人av免费视频| 激情在线观看视频在线高清| 极品教师在线免费播放| 国产在线精品亚洲第一网站| 亚洲国产精品999在线| 两性夫妻黄色片| 国产伦一二天堂av在线观看| 亚洲中文av在线| 国产成人一区二区三区免费视频网站| 欧美日韩精品网址| or卡值多少钱| 亚洲 欧美一区二区三区| 午夜a级毛片| 在线免费观看不下载黄p国产 | 99riav亚洲国产免费| 91在线观看av| 变态另类丝袜制服| 级片在线观看| 国内揄拍国产精品人妻在线| 又爽又黄无遮挡网站| 国内精品美女久久久久久| 欧美激情在线99| 国产精品野战在线观看| 久久婷婷人人爽人人干人人爱| 免费在线观看成人毛片| 我的老师免费观看完整版| 亚洲国产精品成人综合色| 国产精品野战在线观看| 国产高清视频在线观看网站| 成人亚洲精品av一区二区| 大型黄色视频在线免费观看| 十八禁人妻一区二区| 久久久久免费精品人妻一区二区| 亚洲精品国产精品久久久不卡| 日本在线视频免费播放| 啦啦啦免费观看视频1| 国产精品美女特级片免费视频播放器 | 久久久久久国产a免费观看| 成人18禁在线播放| 国产不卡一卡二| 午夜激情福利司机影院| 每晚都被弄得嗷嗷叫到高潮| 久久久精品欧美日韩精品| 国产综合懂色| 久久久成人免费电影| 国产日本99.免费观看| 成人特级黄色片久久久久久久| 亚洲欧美一区二区三区黑人| 国产精品野战在线观看| 美女黄网站色视频| 国产高清videossex| 亚洲成人精品中文字幕电影| 国产真实乱freesex| 变态另类成人亚洲欧美熟女| 国产精品自产拍在线观看55亚洲| 99久国产av精品| 999久久久国产精品视频| 精品久久久久久久末码| 成年版毛片免费区| 亚洲av熟女| 中文资源天堂在线| 色综合站精品国产| 国产av麻豆久久久久久久| 国产免费av片在线观看野外av| 国产97色在线日韩免费| 久久久久久九九精品二区国产| 亚洲精品久久国产高清桃花| 国语自产精品视频在线第100页| 成年女人永久免费观看视频| 久久国产精品人妻蜜桃| 免费av毛片视频| 老司机午夜福利在线观看视频| 亚洲专区中文字幕在线| 午夜精品久久久久久毛片777| 欧美日韩综合久久久久久 | 韩国av一区二区三区四区| АⅤ资源中文在线天堂| 天堂网av新在线| 国产精品久久久人人做人人爽| 首页视频小说图片口味搜索| 日本免费一区二区三区高清不卡| 午夜福利视频1000在线观看| 国产黄片美女视频| 日韩av在线大香蕉| 婷婷亚洲欧美| 欧美在线一区亚洲| 老司机深夜福利视频在线观看| 久久久国产精品麻豆| 99热这里只有是精品50| 美女高潮喷水抽搐中文字幕| 国产精品亚洲av一区麻豆| 国产爱豆传媒在线观看| av女优亚洲男人天堂 | 久久久久性生活片| 此物有八面人人有两片| 久久这里只有精品19| 12—13女人毛片做爰片一| 成人无遮挡网站| 婷婷亚洲欧美| 免费人成视频x8x8入口观看| 91九色精品人成在线观看| 久久久久久国产a免费观看| 淫妇啪啪啪对白视频| 香蕉av资源在线| 国产激情久久老熟女| 国产成人系列免费观看| 久久久国产成人精品二区| 波多野结衣高清无吗| 亚洲国产色片| ponron亚洲| 亚洲欧美日韩高清专用| 999精品在线视频| 俺也久久电影网| 国产精品98久久久久久宅男小说| 亚洲七黄色美女视频| 欧美日韩国产亚洲二区| 国产精品av久久久久免费| 中文字幕高清在线视频| 久久久精品欧美日韩精品| 真人做人爱边吃奶动态| 综合色av麻豆| 在线观看舔阴道视频| 久久久久久大精品| 国产亚洲精品综合一区在线观看| 老汉色av国产亚洲站长工具| 真实男女啪啪啪动态图| 国产成年人精品一区二区| 久久久久精品国产欧美久久久| 99re在线观看精品视频| 熟女少妇亚洲综合色aaa.| 亚洲国产精品sss在线观看| 国产亚洲欧美98| 亚洲精品在线观看二区| 国产97色在线日韩免费| 亚洲国产精品999在线| 成人鲁丝片一二三区免费| 午夜福利在线在线| 最新在线观看一区二区三区| av天堂在线播放| 国产麻豆成人av免费视频| 国产精品综合久久久久久久免费| 欧美色视频一区免费| 身体一侧抽搐| 成人av在线播放网站| 悠悠久久av| 国产淫片久久久久久久久 | 欧美最黄视频在线播放免费| 国产成年人精品一区二区| 日本一本二区三区精品| АⅤ资源中文在线天堂| 一区二区三区国产精品乱码| 亚洲av日韩精品久久久久久密| 国产1区2区3区精品| 亚洲精品美女久久av网站| 啦啦啦观看免费观看视频高清| 亚洲国产欧洲综合997久久,| 欧美一级a爱片免费观看看| 亚洲国产欧洲综合997久久,| 国产成人av激情在线播放| 97超视频在线观看视频| 最近最新中文字幕大全电影3| 日韩欧美一区二区三区在线观看| 欧美中文综合在线视频| avwww免费| 亚洲国产精品久久男人天堂| 久久99热这里只有精品18| 亚洲男人的天堂狠狠| 最新中文字幕久久久久 | 亚洲国产精品999在线| 少妇的丰满在线观看| 精华霜和精华液先用哪个| 色综合欧美亚洲国产小说| 午夜影院日韩av| 女人被狂操c到高潮| 无遮挡黄片免费观看| 国产成人欧美在线观看| 午夜a级毛片| 成人18禁在线播放| 99精品欧美一区二区三区四区| 国内精品久久久久精免费| www.999成人在线观看| 精品一区二区三区av网在线观看| 一二三四社区在线视频社区8| 夜夜躁狠狠躁天天躁| 国产高清视频在线观看网站| 99热精品在线国产| 一个人免费在线观看的高清视频| 91九色精品人成在线观看| 制服丝袜大香蕉在线| 欧美日韩瑟瑟在线播放| 午夜福利在线观看吧| 亚洲国产高清在线一区二区三| 久久午夜亚洲精品久久| 国内精品美女久久久久久| 一级作爱视频免费观看| 亚洲真实伦在线观看| 少妇的丰满在线观看| 一a级毛片在线观看| svipshipincom国产片| 久久精品国产清高在天天线| 国产成人av激情在线播放| 亚洲精品国产精品久久久不卡| 国产成人精品久久二区二区免费| 国产精品自产拍在线观看55亚洲| 老司机深夜福利视频在线观看| cao死你这个sao货| 成人特级黄色片久久久久久久| tocl精华| 男人舔奶头视频| 欧美3d第一页| 99久久精品热视频| 夜夜躁狠狠躁天天躁| 国产亚洲欧美在线一区二区| 搡老熟女国产l中国老女人| 亚洲七黄色美女视频| 黄色日韩在线| 欧美色视频一区免费| 国产三级中文精品| 精品午夜福利视频在线观看一区| 黄色女人牲交| 99国产极品粉嫩在线观看| 久久久久性生活片| 男人舔女人下体高潮全视频| 真实男女啪啪啪动态图| 久久久国产精品麻豆| 又黄又粗又硬又大视频| 欧美黑人巨大hd| 成年女人看的毛片在线观看| 18禁黄网站禁片午夜丰满| 亚洲真实伦在线观看| 亚洲成人久久爱视频| 村上凉子中文字幕在线| 亚洲avbb在线观看| 丰满人妻一区二区三区视频av | 国内精品一区二区在线观看| a在线观看视频网站| 日韩欧美 国产精品| 美女 人体艺术 gogo| 日本与韩国留学比较| 亚洲欧美精品综合久久99| 国产欧美日韩一区二区三| 人妻丰满熟妇av一区二区三区| 午夜福利免费观看在线| 国产精品一区二区三区四区免费观看 | 欧美日韩福利视频一区二区| 在线视频色国产色| 欧美大码av| 麻豆国产97在线/欧美| 亚洲一区二区三区不卡视频| 久久亚洲真实| 日本黄大片高清| 老司机午夜十八禁免费视频| 男人舔女人下体高潮全视频| 亚洲欧美日韩高清在线视频| 久久久精品欧美日韩精品| 国产精品一区二区三区四区免费观看 | 精品一区二区三区四区五区乱码| 日韩 欧美 亚洲 中文字幕| 制服人妻中文乱码| 久久久久久国产a免费观看| 搞女人的毛片| 男插女下体视频免费在线播放| 男人舔女人的私密视频| 国产v大片淫在线免费观看| 午夜免费观看网址| 免费在线观看日本一区| 黄频高清免费视频| 夜夜看夜夜爽夜夜摸| 热99re8久久精品国产| 国产一区二区在线观看日韩 | 亚洲国产日韩欧美精品在线观看 | 亚洲第一电影网av| 亚洲中文av在线| 欧美精品啪啪一区二区三区| 成年人黄色毛片网站| 最好的美女福利视频网| 成人18禁在线播放| 欧美一区二区国产精品久久精品| 亚洲午夜精品一区,二区,三区| 欧美乱妇无乱码| 欧美日本亚洲视频在线播放| 午夜久久久久精精品| 亚洲国产欧美人成| 成人国产综合亚洲| 久久久国产成人免费| 亚洲国产精品999在线| 最新美女视频免费是黄的| 国产三级黄色录像| 白带黄色成豆腐渣| 欧美日韩福利视频一区二区| 少妇丰满av| 国内精品久久久久精免费| 免费观看人在逋| 亚洲中文字幕日韩| 国产97色在线日韩免费| 国产精品免费一区二区三区在线| 中文资源天堂在线| 午夜免费激情av| 99国产极品粉嫩在线观看| 日本在线视频免费播放| 亚洲av日韩精品久久久久久密| 曰老女人黄片| 国产淫片久久久久久久久 | 国产黄片美女视频| 18禁黄网站禁片午夜丰满| 在线观看一区二区三区| 国产伦精品一区二区三区四那| 久久久久久久精品吃奶| 1000部很黄的大片| 成人亚洲精品av一区二区| 天堂√8在线中文| 757午夜福利合集在线观看| 午夜精品在线福利| 天天躁狠狠躁夜夜躁狠狠躁| 男女那种视频在线观看| 亚洲九九香蕉| 国产精品久久久久久久电影 | 老司机深夜福利视频在线观看| 禁无遮挡网站| 国产精品久久久av美女十八| 中文字幕久久专区| 日韩av在线大香蕉| 国产精品亚洲美女久久久| 99精品在免费线老司机午夜| 国产免费av片在线观看野外av| 欧美国产日韩亚洲一区| 国产私拍福利视频在线观看| 日本一本二区三区精品| 99在线人妻在线中文字幕| 久久热在线av| 99热只有精品国产| 欧美最黄视频在线播放免费| 亚洲国产欧洲综合997久久,| 国产 一区 欧美 日韩| 国产精品,欧美在线| 欧美国产日韩亚洲一区| 亚洲乱码一区二区免费版| 免费看日本二区| 国产伦精品一区二区三区视频9 | 久久久久亚洲av毛片大全| 在线免费观看不下载黄p国产 | 国产精品免费一区二区三区在线| 精品无人区乱码1区二区| 日本五十路高清| 最近在线观看免费完整版| 999精品在线视频| 三级毛片av免费| 999久久久国产精品视频| 精品一区二区三区视频在线 | 国产三级中文精品| 女警被强在线播放| 亚洲成人精品中文字幕电影| 法律面前人人平等表现在哪些方面| 免费av毛片视频| 欧美大码av| 国产午夜精品久久久久久| 搡老熟女国产l中国老女人| 日本a在线网址| 麻豆成人av在线观看| www.www免费av| 国产精品综合久久久久久久免费| 国产成人av教育| 亚洲精品在线观看二区| 精品午夜福利视频在线观看一区| 日韩欧美一区二区三区在线观看| 国内精品久久久久精免费| 精品久久久久久久毛片微露脸| 国产伦在线观看视频一区| 母亲3免费完整高清在线观看| 亚洲成人免费电影在线观看| 别揉我奶头~嗯~啊~动态视频| 久久欧美精品欧美久久欧美| 少妇的丰满在线观看| 欧美午夜高清在线| 99精品在免费线老司机午夜| 亚洲九九香蕉| 白带黄色成豆腐渣| a在线观看视频网站| 97碰自拍视频| 国产精品一区二区免费欧美| 91在线精品国自产拍蜜月 | 亚洲精品美女久久av网站| 国产精品98久久久久久宅男小说| 美女 人体艺术 gogo| 在线免费观看的www视频| 欧美绝顶高潮抽搐喷水| 久久午夜亚洲精品久久| 色噜噜av男人的天堂激情| 日韩欧美在线乱码| 不卡av一区二区三区| 亚洲欧美日韩高清专用| 日本 av在线| 亚洲中文日韩欧美视频| 青草久久国产| www日本黄色视频网| 亚洲av成人一区二区三| 这个男人来自地球电影免费观看|