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

    Effect of edge transport barrier on required toroidal field for ignition of elongated tokamak?

    2019-04-13 01:14:32CuiKunYang楊翠坤MingShengChu朱名盛andWenFengGuo郭文峰
    Chinese Physics B 2019年4期
    關鍵詞:文峰

    Cui-Kun Yang(楊翠坤),Ming-Sheng Chu(朱名盛),and Wen-Feng Guo(郭文峰)

    1Institute of Plasma Physics,Chinese Academy of Sciences,Hefei 230031,China

    2University of Science and Technology of China,Hefei 230031,China

    1.Introduction

    Ohmic heating is the easiest and the most effective method of heating the plasma.[1]Therefore using Ohmic heating to achieve plasma ignition is the most desirable ignition path.To realize this goal,we are faced with two challenging technical issues.First,Ohmic heating loses its effectiveness when the plasma temperature gets too high.Therefore,we have to employ large current to reach ignition temperature before Ohmic heating loses its effectiveness.Second,consequently,a strong toroidal field is needed to maintain the plasma equilibrium and stability and also the plasma transport must be reduced to assure that the ignition condition is achievable.Ohmic ignition of the tokamak has been proposed by many authors.[2–6]The critical issue is that the marginal toroidal field(BT)for a plasma with elongation κ~1 is required to be as high as about 20 T.[3]This is at or beyond the limit of the existing technologies.It was proposed in Ref.[4]that the marginal BTcan be reduced by elongating the tokamak.It was shown that with elongation~5,the marginal BTcan be reduced to~6 T.Then it is well within the limit of existing technology.

    The zero-dimensional(0D)system models,when properly bench-marked against full transport codes,have been found to be extremely useful in design studies.If we generate 0D models based on the calibrated one-dimensional(1D)profile models,we can obtain some extra information regarding the effects of profile on the 0D models.Especially,if the profiles are modified from the profile used in the bench mark,then we can use this to obtain perturbative information regarding the effect of these modified profiles.We apply this method to the study of Ohmic ignition of an elongated tokamak.In previous study,a 0D ignition model was generated by using a set of 1D profiles(This set of profiles is defined as ‘I’model hereafter.)The ignition phenomenon is analyzed by using the 0D model.Particular cases were then bench-marked against transport code runs,establishing its validity.However,the study was made relatively early,before the identification of the possibility of an edge transport barrier.[7,8]It is natural to inquire whether the presence of the edge barrier can help us reduce the required marginal BTfor ignition.This is the main purpose of the present paper.

    In the present paper,we study this question by generalizing the 1D profiles(‘I’model)that were used in Ref.[4]to include an edge transport barrier.Most of the analyses made in Ref.[4]can be readily repeated because the structure of the problem is not changed.A seemingly counter intuitive conclusion is reached:the effect of an edge transport barrier would increase the marginal BTfor ignition.Of course,in existing experiments,the achieved edge transport barrier is relatively small.The resulting small increase can be absorbed into the design margins.But more extreme profile modifications have been proposed.Also,we try to understand the physics reason underlying the result of this new discovery.

    To obtain our conclusion regarding the effect of edge barrier,we implement two possible edge barrier models.The first edge barrier model is that the edge barrier effect on the pro-file is assumed to be mainly limited to the edge;it does not penetrate all the way into the center.This is close to the case where the extreme edge is cleaned by using a lithium wall.We call this profile model the ‘DC’profile model.The second edge barrier model is the model that is induced by the H-mode.In this type of edge transport barrier,the whole profile can be considered to be lifted up by an amount that is equal to the amount caused by the edge barrier.This is called the‘H’profile model.Both types of profile modifications imply‘improved edge confinement’.Although the profile modifications improve the edge confinement to different extents,they lead to similar results of‘higher’marginal BTrequired for ignition.One hypothesis is that the edge transport barrier actually results in broadening the profiles of the various physical quantities such as density,temperature,and current.Although broadness and peakedness are intuitively clear concepts,to be more precise,we define the broadness br of a profile ξ (ρ)as

    The peakedness is pk=1?br.(It is quite clear that br=1and pk=0 for a flat profile with ξ =ξ(0).For a profile that is peaked at the center ξ= ξ (0)δ(0),with δ(0)being the Kroneker delta function,then br=0,pk=1).We note that the broadness defined here is closely related to the average of the profile across the plasma cross-section,and for a large aspect ratio torus it is also the same as the average of the profile across the plasma volume.It is specifically defined here to make our discussion more precise.

    To further strengthen our hypothesis,we go on testing more general profiles and we obtain the following conclusion:a broader profiles will lead to a higher required marginal BT,whereas a more peaked profile can bring about a lower required marginal BT.This conclusion is further confirmed by using a set of profiles that are more general than those used in Ref.[4].Therefore,introducing an edge transport barrier is just a means in broadening the plasma profile.This phenomenon is due to the concentration of the fusion energy production needed for ignitionin near the center of the plasma.Broadening the profiles actually moves the fusion energy input away from the center,which is required for ignition.A fictitious model in which the energy input profile could be modified is examined.It is found that the trend of the required margin BTinfluenced by the transport barrier can be substantially modified.Application of this result can lead to the benefits expected from the transport barrier.

    The rest of this paper is organized as follows.In Section 2,we briefly review the model and method of analysis of Ref.[4].In Section 3,we analyze the effect of profile variation on ignition.First we introduce two different models to study the effect of incorporating the edge transport barriers.In Subsection 3.1 we study ‘DC’type of profiles,in which the effect of the edge barrier does not propagate to the plasma center.In Subsection 3.2,we study the effect of edge barriers in which the elevated edge barrier propagate all the way into the plasma center,resulting in an elevated plasma ‘H’profile typical in H-mode type of discharge.In Subsection 3.3,we modify the profiles in Ref.[4]to include a larger set of profiles,either broader ‘GB’or more peaked ‘GP’,but without a transport barrier.We show that,in general,a broader profile leads to a more stringent requirement for BT;whereas a more peaked profile will reduce the required BT.In Section 4,we introduce a fictitious model,in which the fusion power input could be adjusted to be concentrated in or moved away from the center.We show in this model that when the fusion power is much less concentrated in the center,we can possibly utilize the edge pedestal more to our advantage,i.e.,reducing the required BT.A brief summary and discussion are given in Section 5.

    2.Summary of model and analysis in Ref.[4]

    Here we follow the calibrated model given by Chu et al.in Ref.[4]to formulate the Ohmic ignition problem.The local power generation and power loss of the plasma are related by

    In Eq.(2),P?,Pf,Pb,Pe,and Piare the Ohmic power,fusion input,bremsstrahlung loss,electron transport loss,ion transport loss of the plasma,respectively.Here,the MKS unit is used.The P,n,and T are in units of W,m?3,and keV,respectively.In the density and temperature range of interest,these power input and loss components can be expressed as P?=CjT?3/2,Pf=Cfn2T3,Pb=Cbn2T1/2,and Pe=CeT,Pi=Cin2T1/2.Using Ohmic heating to achieve ignition relies on high plasma density and relatively low(in comparison to other heating processes)plasma temperature.Therefore we assume that the electron and ion temperature are equilibrated on the flux surface.The C coefficients depend on other characteristics of the system,such as j,a,A,R0,κ,q,...,which denote the current density,minor radius,aspect ratio,major radius,elongation,safety factor,...).Because,in general,each of n,T,j,q,κ,...is a function of flux surface.Equation(2)depends on flux surface.To create a 0D model out of Eq.(2),we take the average of Eq.(2)over the plasma volume,i.e.,we define ignition as

    Here,the over bar represents the average over plasma vol-To facilitate analytic parametric study,it is assumed that the large aspect ratio,and uniform plasma elongation limit of the flux surface are valid;and the major variation across the flux surfaces comes from the density,temperature and current profiles. A library of profiles for the physical variables ξ={n,T,j}=(density,temperature,current density)can be represented as

    Here, ρ2=x2/a2+y2/a2κ2, ξ0={n0,T0,j0}and p={pn,pT,pj}.(Here the lower case p is used to denote the exponent of the profile)The subscript I specifies that this is the model used in Ref.[4].Then the volume averaged power generation and loss areA profile average function is defined as

    Utilizing the F function,the C coefficients are written in terms of one part that depend strongly on the profile,the K function;and another part that depends very mildly on the profile.We obtain

    with

    with

    with

    with

    with

    Here,Zeff,,K,q,and μ denote the effective ion charge,κBT/qA,2κ2/1+κ2,the safety factor,and the ion mass number,respectively,and BTis the toroidal field,Ftris the fraction of trapped particles,and α is the fraction of α particle energy that is absorbed locally by the plasma.

    If all the’s values are known,then the value of P be comes a function ofand.At marginal ignition,the functionsatify the condition of a saddle point(see the panel on the left side in Fig.1).The system parameters used here are α =1,μ =2.5,Zeff=1,a=0.25,q=2,κ =5,and A=3.(In the following,we will regard this as the standard case)Therefore,the marginal point for ignition in the)space is defined as

    These conditions are solved analytically,thus resulting in a constraint on thecoefficients as follows:

    Once equation(7)is satisfied,the marginal values of density and temperature can be obtained in terms of theas given below.

    Concurrently,the condition P?=Pe=Pf=(Pb+Pi)was also satisfied.Equation(7)can be expressed alternatively as a constraint on the physical values which are the input values to the problem.If we regard a,A,Zeff,κ,μ ...as known variables,and BTas the only unknown variable,then its marginal value is given by

    Here,

    With BTas defined in Eq.(9),we can readily evaluate the numerical values of the C’s and consequentlyAll other physical quantities such as total current,density,and temperature profiles are also obtained for an elongated tokamak(see the red curves in Fig.2).Therefore BTor equation(9)can be regarded as the key to the ignition condition.Usually the Ohmic ignition is considered as a technologically challenging route to achieving the ignition because for a circular tokamak with κ=1,the required BT(see Fig.1)is around 20 T,at or beyond the practical engineering limit.

    One of the main contributions of Ref.[4]is that it shows that at large κ the required toroidal field is scaled as 1/κ.The result of the model is also calibrated against a transport code with pn=pT=1,and pj=1.5,for a plasma with α=1,μ=2.5,Zeff=1,and a=0.25.Although at κ=1,the required toroidal field is≈ 20 T,whereas for κ =5,the required toroidal field for ignition is lowered to≈6 T.

    We note here that the profiles adopted in this model,Eq.(4),do not include the effect of a transport barrier at the plasma edge.Therefore,our goal of this paper is to study the effect of an edge transport barrier.[7,8]

    3.Effect of edge transport barriers and more general types of profile

    There has been two different methods of producing an edge transport barrier.The first kind is usually created by extreme edge cleaning,for instance by using lithium coated wall,or other pumping methods.In this type of edge transport barrier,the profile change is more limited to the plasma edge.The effect of this type of transport barrier is to be discussed in Subsection 3.1.The second type of edge transport barrier is created by an H mode through large power input.The edge confinement improvement is most likely due to the excitation of large poloidal flow in the edge region which suppresses the anomalous transport.This type of edge barrier actually gives rise to a substantial increase in the total energy content. The effect of this type of transport barrier is discussed in Subsection 3.2.In principle,these two types of edge-barriers are probably not mutually exclusive.But here,we regard them as two pure cases.We note that both types of transport barrier lead to the broadening of the profiles of the physical variables,such as density,temperature and current density. They lead to a similar higher marginal BTrequired for igntion.We also subsequently invoke an even wider class of profiles than those employed in Ref.[4]and summarized in Section 2 for the case without a transport barrier,in which the profile can be either broader or more peaked.We show that the effects of profile changes in all cases invoked in this section have one common similarity.That is,in general,a broader profile will cause the required marginal BT,to increase,which is the cases discussed in Subsections 3.1 and 3.2.On the other hand,if the profile is more peaked than that used in Section 2,then the required marginal BTcan be reduced.We suggest that this is due to the fact that the ignition condition is facilitated by concentrating the fusion energy input in the plasma center.

    3.1.Effect of discharge cleaning type of edge transport barrier on BT

    In discharge cleaning induced edge transport barrier,we assume that it gives rise to a pedestal in profile,but this change does not propagate into the plasma center.For this purpose,we propose to use the function fDCto model the effect of this discharge-cleaning ‘DC’type of transport barrier.Here

    We see that at h=0,fDCreduces to fI(see Eq.(4)).At 0

    These averages modify the definitions of the K functions,but do not change the relationship betweenand K.With the modified values of K,we can still use the same method outlined in Section 2 to obtain the expressions for the various quantities at ignition,i.e.,the formulae are the same,although their numerical values have changed,because the values of K have changed,thus thehave changed.

    The results of our study indicate that with fDCtype of profile,the resultant BTis always larger than that of the original fIprofile.We show in Fig.1 the comparison of the saddle point structure of the averaged power in the(n T)plane.Figure 1(a)shows the contour of averaged power of original profiles using fI; figure 1(b)displays the averaged power for a plasma with an fDCtype of edge barrier with ‘h=0.2’.We observe that the edge pedestal moves the ignition point towards higher averaged density and higher averaged temperature.

    Fig.1.(a)Original profile ‘I’model and(b)discharge cleaning profile‘DC’model(see Eq.(10))with ‘h=0.2’,showing that bothare shifted to higher values in ‘DC’profile than in ‘I’profile.

    3.2.H mode-induced edge transport barrier

    In the H mode,the plasma develops an edge transport barrier which propagates into the center of the plasma.This results[7]in an improved overall plasma confinement.A simplified but useful assumption is made for this model as follows:the whole profile increases by the same amount as the pedestal.This actually also broadens the plasma profile.We can model this type of edge transport barrier as

    Quite obviously,if h=0,fHdecreases to fI(see Eq.(4)).At h<1,the edge value increases by a value of hp,and the value in the center is(1+h)p.If h=1,the profile does not become a constant as in the case of Subsection 3.1.The profile average of fHis

    To facilitate our understanding of these two different types of profiles,we show in Fig.2 their comparison with the original profile.It is interesting to note that for fIn=(1?ρ2),br=pk=0.5;for the DC profile with fDCn=(1?(1?h)ρ2),br=0.5(1+h),pk=0.5(1?h);and for the H profile with fHn=(1?ρ2+h),br=(0.5+h)/(1+h),pk=0.5/(1+h).Therefore,the increase in the edge pedestal will always result in a broader profile.

    The effect of fHtype of profile in affecting various quantities at ignition is shown in Fig.3.The system parameters are the same as those used in Fig.1.Figure 3(a)shows the values of the marginal BT, figure 3(b)the values of the total current IP, figure 3(c)the values of the density,and figure 3(d)the values of averaged β.The red line is for the original case,and the blue line for the case with an H-type edge barrier with‘h=0.2’.The black curves refer to the Troyon β limits.[9]

    Fig.2.Comparison of the profiles in the ‘I’profile with that in the ‘DC’model and in the ‘H’model.Note that both DC and H models have non-zero values at the edge.Also for H profile,central value exceeds 1.Solid red curves represent density profiles,the dotted blue curves refer to temperature profiles,and green curves indicate current density profiles.

    Fig.3.Physical parameters at ignition for‘I’model(red)and ‘H’model(blue)for aspect ratio=3(solid)and 7(dashed)plasma versus elongation κ.For‘H’model,h=0.2(see Eq.(12))is used.Panel(a)shows marginal toroidal field BTin unit T,panel(b)exhibits toroidal current IPin MA,panel(c)refers to average densitn 1020/m3and panel(d)indicates volume averaged β values.The black curves represent Tryon limits.

    We see that both types of edge transport barriers at this modest barrier height result in higher required BT,higher total current,higher density,and higher β value.We believe that this is due to the fact that both types of transport barriers result in broadening the profile.

    3.3.Effects of more general types of profiles(broader or more peaked)on marginal BBBTTT

    In this subsection,we use an even more extended set of profiles than the ’I’models which are summarized in Section 2 to show that the peakedness of profile will lead to the improvement of ignition condition or the reduction of the toroidal field required for ignition.A broader profile will lead to an increase in the required BT.We consider the following profiles:

    We note that this set of functions is the same as the original set if l=1.If l>1,then it is broader;if l<1,then it is more peaked.The average of this profile over the plasma cross-section gives

    Fig.4.Comparisons among variations of marginal BTrequired for all models at ignition relative to the ‘I’model for an aspect ration A=3 plasma,with edge parameter h(see Eqs.(10)and(12)).In order to accommodate all models in the same figure,we also used h to denote l±1.(see Eq.(14))for the ‘GB’(broader)and ‘GP’(more peaked)models.Panel(a)is the change of BTin unit of T,and panel(b)is the change relative to that of the ‘I’model(see Eq.(4)).

    Fig.5.(a)and(c)Plots of different components of power versus squared minor radius.(b)and(d)total sum of power input,power loss,and net power.Values for standard ‘I’profile are plotted in panels(a)and(b),and values for the ‘H’profile with A=3,κ =5 are plotted in panel(c)and panel(d).

    or

    Here,B is the Beta function,i.e.,

    and it is readily evaluated.Therefore,we have the generalized FGfunction as given below.

    We note that this generalization of the profile function is different from that of the pedestal function introduced in Subsections 3.1 and 3.2,shown in Fig.4 is the comparison of the effect among the different models defined up to now in affecting the marginal BTas a function of the barrier height h.It is seen that if the barrier height is increased, then the marginal BTis raised;similarly,if the profile is broader,then the marginal ignition BTis raised.Because raising the barrier height is to broaden the profile,we conclude that,in general,broadening the profile leads to raising the BTand vice versa.We attribute this to the fact that most of the fusion power input required for the ignition is located in the center of the plasma.

    Figure 5 shows the distributions of various components of power generation(loss)versus variable ρ2.It is clear that the fusion power is the most concentrated component in the plasma center.Therefore,as expected,the main cause of ignition is due to the fast rise of fusion power production near the plasma center which keeps the central temperature rising,and eventually leads to ignition.Raising the edge density and temperature everywhere by the same amount would not be conducible to the ignition.The interesting question is“Can we change the situation if the profile of fusion power generation is modified?”

    4.Effects of modified energy deposition profile

    To further investigate the basic cause of‘failure to reduce the required BTin plasma in which the edge density and temperature are raised’,in this section we introduce a fictitious energy deposition model to clarify this issue. Because we hope to attribute the result to energy production being too concentrated in the center of the plasma,a model in which the energy deposition profile is more flexible should be useful.From Fig.5 we see that only the fusion power is very concentrated in the plasma center,therefore we propose to consider the following modification model for the fusion energy production

    Here,w represents the replacement of center energy generation by averaged energy generation.It is seen that if w=0,we recover the original distribution of energy generation.Ifw>0,the energy distribution is switched from the central region near the magnetic axis towards the outside;whereas if w<0,the energy generation becomes more concentrated in the plasma center.This modified energy generation can modify only the Kfcoefficients.In this case,

    With this modified model,when w increases,the minimum BTincreases with w.However,the effect of edge transport barrier on BTis reduced.As a matter of fact,consider specifically the case of w=0.7,then the fusion and the Ohmic power deposition profile will become the same,then increasing the edge pedestal actually becomes favorable.At this time the power balance as a function of volume becomes more favorable at the mid portion of the plasma,and unfavorable in the center.The various components of the energy deposition profiles for w=0.7 are shown in Fig.6.With this value of w,we should expect the edge barriers to have a reduced BTfor the‘DC’and‘H’model relative to the ‘I’model.

    Fig.6.(a)Different components of power versus squared minor radius for modified fictitious plasma ignition model with w=0.7(seeEq.(17))profile with A=3,and κ=5;(b)total sum of power input,power loss,and net power.By comparing with Fig.5,we see that the fusion power input is now the same as the Ohmic power input,and less peaked than the Bremsstrahlung or the ion loss.Also the net power is not peaked in the center but in the mid-range of plasma radius.At this time,the edge pedestal leads to a reduced required BT.

    Figure 7(a)shows the plots of BTrequired for ignition versus fusion redistribution parameter w for the various models.We can see that as expected,when the fusion power becomes less peaked in the center,the required toroidal fields all increase.Figure 7(b)shows the plots of the required BTrelative to the original model versus fusion redistribution parameter w.We can see that as expected,when w increases,the required BTvalue for each of the ‘DC’and ‘H’model gradually becomes smaller.At the same time,the BTfor model‘GB’also becomes smaller,while that of model‘GP’becomes larger than that of the original model.As expected,the values of BTand delta BTchanging over actually occur at w=0.45 before reaching w=0.7.

    Therefore we conclude that the counter intuitive effect of the edge pedestal on the required minimum BTis due to the fact that the fusion power deposition is extremely concentrated in the plasma center.

    Fig.7.(a)Plots of BTrequired versus fusion power redistribution parameter w for various models,showing that as the fusion deposition power moves away from the center,the required BTincreases for all models.(b)Plots of difference in required BTrelative to that of the original model used in Ref.[4],versus fusion power redistribution parameter w,indicating that as w increases,values of required BTfor‘DC’,‘H’,and ‘GB’models are all improved with respect to those of original model and ‘GP’model becomes less favorable.The values of BTand delta BTchanging over actually occur around w=0.45 before reaching 0.7.

    5.Summary and conclusions

    Ohmic ignition remains an attractive option to achieve ignition in fusion devices.[6]The major challenge is the large toroidal field required.It was proposed in a previous research[4]that this major difficulty could be circumvented by using an elongated tokamak.The plasma profiles employed in the previous study did not include the presence of a transport barrier.In this work,we examined a wide class of different profiles,including the possibility of a transport barrier as shown in Table 1.

    The comparisons among different types of profiles are shown in Fig.8.

    We found that broader profiles would lead to higher BT:more peaked profiles would lead to smaller values of BT.This phenomenon is not restricted to whether or not an edge transport barrier exists.We also generalize the ignition problem to the case where there are different fusion energy deposition profiles,and show that the relationship between the required BTand the profile peakedness is a direct consequence of the fusion power generation concentrating near the plasma center.When a sufficient amount of fusion power generation is moved away from the plasma center towards the outer parts of the plasma radius,an edge transport barrier then leads to a reduced BTrelative to profile without a barrier.However,moving the fusion power away from the center,in general,leads to a higher BTregardless of profile models.Only when the power input becomes less concentrated at the plasma center,then the transport barriers show advantages with respect to profiles without a transport barrier.

    Table 1.Profiles and broadness for various models.

    Fig.8.Comparisons among different plasma profiles studied with h=0.2 and l=1+h=1.2 for GB model and l=1+h=0.8 for GP model.

    Acknowledgment

    This work was done on the ShenMa High Performance Computing Cluster at the Institute of Plasma Physics,Chinese Academy of Sciences.

    [1]Li J G 2016 Physics 45 88(in Chinese)

    [2]Coppi B 1977 Comments Plasma Phys.Control.Fusion 3 47

    [3]Wagner C E 1981 Phys.Rev.Lett.46 654

    [4]Chu M S,Harvey R W,Politzer P A and Yamagishi T 1985 Nucl.Fusion 25 835

    [5]Cohn D R and Bromberg L 1986 J.Fusion Energy 5 161

    [6]https://wikivisually.com/wiki/Ignitor

    [7]Wagner F,Fussmann G,Grave T,et al.,1984 Phys.Rev.Lett.53 1453

    [8]Behn R,Alfier A,Yu S,Zhuang G,Pasqualotto R,Nielsen P,Martin Y and team T C V 2007 Plasma Phys.Control.Fusion 49 1289

    [9]Troyon F,Gruber T,Saurenmann H,Semenzato S and Succi S 1984 Plasma Phys.Control.Fusion 26 1A 209

    猜你喜歡
    文峰
    Design of broadband achromatic metasurface device based on phase-change material Ge2Sb2Te5
    近30年國內外籃球研究現(xiàn)狀與發(fā)展趨勢
    體育師友(2022年1期)2022-04-17 10:42:34
    留暉(荊文峰)
    寶藏(2020年12期)2021-01-21 02:15:54
    胡文峰博士簡介
    中東省蘭陵縣尚巖鎮(zhèn)初級中學紅色文峰文學社
    中學時代(2020年1期)2020-02-17 10:38:38
    “期貨交易”型走私犯罪的刑法分析
    犯罪研究(2019年3期)2019-11-27 19:28:54
    小畫廊
    小主人報(2018年17期)2018-09-12 07:21:34
    文峰街
    重慶與世界(2016年6期)2016-10-09 06:27:10
    油庫是我家
    洗腦
    国产 精品1| 国产一区有黄有色的免费视频| 高清在线视频一区二区三区| 免费人成在线观看视频色| 午夜福利,免费看| 国产在线视频一区二区| 中文欧美无线码| 精品午夜福利在线看| 国产黄片视频在线免费观看| 国产精品一二三区在线看| 天天影视国产精品| 亚洲av日韩在线播放| 国产成人91sexporn| 欧美日韩视频高清一区二区三区二| 国产精品一区www在线观看| 久久久a久久爽久久v久久| 国产免费现黄频在线看| 久久精品久久久久久噜噜老黄| 日本av手机在线免费观看| 国产老妇伦熟女老妇高清| 精品亚洲乱码少妇综合久久| 九九在线视频观看精品| 各种免费的搞黄视频| 成人无遮挡网站| 亚洲av综合色区一区| 蜜桃在线观看..| 男的添女的下面高潮视频| xxx大片免费视频| 最近手机中文字幕大全| 蜜桃久久精品国产亚洲av| 欧美日本中文国产一区发布| 国产熟女午夜一区二区三区 | 建设人人有责人人尽责人人享有的| 一级毛片aaaaaa免费看小| 亚洲第一区二区三区不卡| 极品人妻少妇av视频| 啦啦啦视频在线资源免费观看| 蜜桃国产av成人99| 99国产精品免费福利视频| 久久99精品国语久久久| 日本爱情动作片www.在线观看| 三级国产精品片| 80岁老熟妇乱子伦牲交| 久久99蜜桃精品久久| av天堂久久9| 亚洲人与动物交配视频| 亚洲丝袜综合中文字幕| 国产极品天堂在线| 国产黄片视频在线免费观看| 欧美 日韩 精品 国产| 欧美老熟妇乱子伦牲交| 视频在线观看一区二区三区| 91久久精品国产一区二区三区| 精品人妻熟女av久视频| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | 少妇熟女欧美另类| a 毛片基地| 日韩视频在线欧美| 亚洲精品久久久久久婷婷小说| 亚洲av福利一区| 最近手机中文字幕大全| 日韩一区二区视频免费看| 日本黄大片高清| 亚洲精品色激情综合| 性色avwww在线观看| 久久精品国产a三级三级三级| 亚洲,欧美,日韩| 91aial.com中文字幕在线观看| 中国国产av一级| 人妻少妇偷人精品九色| 国产不卡av网站在线观看| 久久免费观看电影| 中文字幕最新亚洲高清| 人妻系列 视频| 国产精品国产av在线观看| 在线看a的网站| 国产欧美另类精品又又久久亚洲欧美| 亚洲国产色片| av视频免费观看在线观看| 男的添女的下面高潮视频| 国产精品人妻久久久影院| av不卡在线播放| av在线老鸭窝| 亚洲精品亚洲一区二区| 午夜免费鲁丝| 高清毛片免费看| 国产精品一区二区在线观看99| 亚洲美女视频黄频| 人妻人人澡人人爽人人| 成人综合一区亚洲| 日韩一本色道免费dvd| 一级a做视频免费观看| 欧美xxxx性猛交bbbb| 国产成人午夜福利电影在线观看| 日产精品乱码卡一卡2卡三| 青春草亚洲视频在线观看| 国产视频内射| 女性生殖器流出的白浆| 久久久久视频综合| 国产片特级美女逼逼视频| 在线观看免费视频网站a站| 欧美精品一区二区大全| 97在线视频观看| 少妇熟女欧美另类| av有码第一页| 一级毛片 在线播放| 夜夜骑夜夜射夜夜干| 欧美一级a爱片免费观看看| 永久网站在线| 亚洲国产精品一区三区| 99热这里只有精品一区| 免费高清在线观看日韩| 2018国产大陆天天弄谢| 久久精品国产亚洲网站| 成人免费观看视频高清| 桃花免费在线播放| 五月天丁香电影| 久久精品熟女亚洲av麻豆精品| 成人毛片a级毛片在线播放| 高清av免费在线| 熟女人妻精品中文字幕| 亚洲欧美成人精品一区二区| 精品国产乱码久久久久久小说| 亚洲精品av麻豆狂野| 人妻 亚洲 视频| 久久久午夜欧美精品| 欧美性感艳星| 久久久久久久久久久丰满| 插逼视频在线观看| 激情五月婷婷亚洲| 狠狠精品人妻久久久久久综合| 国产探花极品一区二区| 久久久久久久久大av| 99九九在线精品视频| 高清在线视频一区二区三区| 边亲边吃奶的免费视频| 国产高清国产精品国产三级| 亚洲精品456在线播放app| 日韩,欧美,国产一区二区三区| 制服诱惑二区| 人妻人人澡人人爽人人| av专区在线播放| freevideosex欧美| 在线观看免费高清a一片| 十八禁高潮呻吟视频| 久久人妻熟女aⅴ| 久久久国产一区二区| 一级黄片播放器| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | 国产一区有黄有色的免费视频| 人成视频在线观看免费观看| 91在线精品国自产拍蜜月| av在线app专区| 亚洲欧美一区二区三区黑人 | 视频区图区小说| 最后的刺客免费高清国语| 麻豆乱淫一区二区| 人妻制服诱惑在线中文字幕| 99热国产这里只有精品6| 黑人猛操日本美女一级片| 狂野欧美白嫩少妇大欣赏| 国产毛片在线视频| 国产熟女午夜一区二区三区 | 久久久亚洲精品成人影院| 青青草视频在线视频观看| 欧美日韩亚洲高清精品| 制服丝袜香蕉在线| 最近中文字幕高清免费大全6| 亚洲国产精品国产精品| 国产精品成人在线| 久久国内精品自在自线图片| 色5月婷婷丁香| 少妇人妻 视频| 老司机影院成人| 在线播放无遮挡| 成人国产麻豆网| 精品一区二区三区视频在线| 少妇高潮的动态图| 亚洲内射少妇av| 国产成人aa在线观看| 国产毛片在线视频| 日韩中文字幕视频在线看片| 亚洲欧美日韩另类电影网站| 最近中文字幕2019免费版| 777米奇影视久久| 韩国av在线不卡| 性高湖久久久久久久久免费观看| 五月玫瑰六月丁香| 女的被弄到高潮叫床怎么办| 久久久久久久久久久久大奶| 亚洲,一卡二卡三卡| 欧美精品高潮呻吟av久久| 久久久午夜欧美精品| 午夜福利视频精品| 亚洲三级黄色毛片| 在线观看免费日韩欧美大片 | 人成视频在线观看免费观看| 国产精品久久久久久av不卡| 考比视频在线观看| 一本色道久久久久久精品综合| 狂野欧美白嫩少妇大欣赏| 性高湖久久久久久久久免费观看| 成人影院久久| 精品亚洲成国产av| 欧美bdsm另类| 一边摸一边做爽爽视频免费| 老司机影院毛片| 国产精品免费大片| 嘟嘟电影网在线观看| 国产有黄有色有爽视频| 国产精品女同一区二区软件| 男女无遮挡免费网站观看| 亚洲熟女精品中文字幕| 国产片特级美女逼逼视频| 在线精品无人区一区二区三| 99久久中文字幕三级久久日本| 熟女人妻精品中文字幕| 黄色一级大片看看| 制服人妻中文乱码| 日本欧美视频一区| 成人国产麻豆网| 国产精品一区www在线观看| 看免费成人av毛片| 18禁在线播放成人免费| 国产淫语在线视频| 日韩不卡一区二区三区视频在线| 韩国高清视频一区二区三区| 久久久久久久久久久久大奶| 中文字幕精品免费在线观看视频 | 男人操女人黄网站| 亚洲精品日本国产第一区| 人人妻人人爽人人添夜夜欢视频| 亚洲一级一片aⅴ在线观看| 国产片特级美女逼逼视频| 亚洲国产日韩一区二区| 黑人猛操日本美女一级片| 亚洲图色成人| 91国产中文字幕| 夜夜骑夜夜射夜夜干| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | 99久久精品一区二区三区| 亚洲天堂av无毛| 国产亚洲精品久久久com| 午夜福利,免费看| 中国三级夫妇交换| 18禁动态无遮挡网站| 亚洲欧美一区二区三区黑人 | 国产男人的电影天堂91| 熟妇人妻不卡中文字幕| 蜜臀久久99精品久久宅男| 99九九在线精品视频| 中国美白少妇内射xxxbb| 我的女老师完整版在线观看| 美女视频免费永久观看网站| 纵有疾风起免费观看全集完整版| 亚洲婷婷狠狠爱综合网| 91精品三级在线观看| 一级片'在线观看视频| 纯流量卡能插随身wifi吗| 国产成人精品久久久久久| 国产黄片视频在线免费观看| 大又大粗又爽又黄少妇毛片口| 日本午夜av视频| 女性被躁到高潮视频| 国产在线视频一区二区| 在线观看国产h片| 国产日韩欧美在线精品| 69精品国产乱码久久久| 日韩亚洲欧美综合| 国产欧美日韩一区二区三区在线 | 亚洲精品456在线播放app| 交换朋友夫妻互换小说| 春色校园在线视频观看| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | 久久鲁丝午夜福利片| 欧美精品国产亚洲| av黄色大香蕉| 看非洲黑人一级黄片| 水蜜桃什么品种好| 日本黄大片高清| 卡戴珊不雅视频在线播放| av专区在线播放| 成人国产av品久久久| 夫妻性生交免费视频一级片| 一区二区三区四区激情视频| 国产不卡av网站在线观看| 国产精品人妻久久久影院| 国产欧美日韩一区二区三区在线 | 丝袜在线中文字幕| 人成视频在线观看免费观看| 在线 av 中文字幕| 欧美日韩综合久久久久久| 人人妻人人添人人爽欧美一区卜| kizo精华| 国产国拍精品亚洲av在线观看| 18禁观看日本| 午夜av观看不卡| 亚洲国产精品成人久久小说| 在线看a的网站| 性色avwww在线观看| www.av在线官网国产| 久久午夜综合久久蜜桃| 秋霞伦理黄片| 欧美精品国产亚洲| 欧美一级a爱片免费观看看| 亚洲伊人久久精品综合| 久久精品国产a三级三级三级| 国产精品一国产av| 免费看av在线观看网站| 大香蕉久久成人网| 国产日韩欧美视频二区| 91午夜精品亚洲一区二区三区| 热99国产精品久久久久久7| .国产精品久久| 久久青草综合色| 亚洲无线观看免费| 伦理电影免费视频| 亚洲精品一二三| 九色成人免费人妻av| 国产精品麻豆人妻色哟哟久久| 少妇丰满av| 亚洲av综合色区一区| freevideosex欧美| 麻豆乱淫一区二区| 只有这里有精品99| 国产视频内射| 国产午夜精品久久久久久一区二区三区| 亚洲人与动物交配视频| 男女边摸边吃奶| 国产片特级美女逼逼视频| 夜夜看夜夜爽夜夜摸| 精品少妇内射三级| 黄色一级大片看看| 91精品国产国语对白视频| 亚洲性久久影院| 男女国产视频网站| a级片在线免费高清观看视频| 亚洲内射少妇av| 免费少妇av软件| 97超视频在线观看视频| 天天影视国产精品| 少妇的逼好多水| 欧美一级a爱片免费观看看| 99久久中文字幕三级久久日本| 乱人伦中国视频| 人人妻人人澡人人看| 女的被弄到高潮叫床怎么办| 欧美bdsm另类| 亚洲图色成人| 一级毛片 在线播放| 伊人久久国产一区二区| 美女国产视频在线观看| 伊人久久国产一区二区| 最近最新中文字幕免费大全7| 自线自在国产av| 精品国产乱码久久久久久小说| 欧美成人精品欧美一级黄| 欧美精品一区二区免费开放| 亚洲精品中文字幕在线视频| 日本午夜av视频| 欧美亚洲日本最大视频资源| av播播在线观看一区| 国产视频内射| 亚洲熟女精品中文字幕| 久久久久久久国产电影| 亚洲美女视频黄频| 精品亚洲成a人片在线观看| 日本wwww免费看| 中文字幕免费在线视频6| 一区在线观看完整版| 婷婷色av中文字幕| 26uuu在线亚洲综合色| 香蕉精品网在线| 制服诱惑二区| 特大巨黑吊av在线直播| 久久精品久久久久久噜噜老黄| 久热久热在线精品观看| 亚洲精品日韩av片在线观看| 狠狠精品人妻久久久久久综合| 美女主播在线视频| 亚洲人成网站在线播| 婷婷色综合大香蕉| 三上悠亚av全集在线观看| 日韩伦理黄色片| 纯流量卡能插随身wifi吗| 99热国产这里只有精品6| 自线自在国产av| 亚洲美女视频黄频| 午夜福利视频在线观看免费| 亚洲综合色网址| 久久精品熟女亚洲av麻豆精品| 亚洲内射少妇av| 日韩欧美一区视频在线观看| 精品一区在线观看国产| 色5月婷婷丁香| 天美传媒精品一区二区| 人妻制服诱惑在线中文字幕| 精品卡一卡二卡四卡免费| 99国产综合亚洲精品| 插阴视频在线观看视频| 成人毛片60女人毛片免费| 黄色视频在线播放观看不卡| 国产极品天堂在线| 欧美成人精品欧美一级黄| 精品国产乱码久久久久久小说| 99热这里只有精品一区| 最近最新中文字幕免费大全7| 欧美精品高潮呻吟av久久| 欧美精品国产亚洲| 国产成人午夜福利电影在线观看| 美女xxoo啪啪120秒动态图| 人人妻人人爽人人添夜夜欢视频| 免费观看性生交大片5| 黄色一级大片看看| 国语对白做爰xxxⅹ性视频网站| 日本-黄色视频高清免费观看| 亚洲精品日本国产第一区| 欧美精品人与动牲交sv欧美| 国产乱来视频区| 午夜av观看不卡| 久久久久久人妻| 桃花免费在线播放| 精品一品国产午夜福利视频| 在线观看国产h片| 91精品三级在线观看| 91精品国产国语对白视频| 99精国产麻豆久久婷婷| 亚洲精品日本国产第一区| 男女啪啪激烈高潮av片| 超碰97精品在线观看| av福利片在线| 极品少妇高潮喷水抽搐| 日韩av免费高清视频| 18+在线观看网站| 观看美女的网站| 岛国毛片在线播放| 一本大道久久a久久精品| 边亲边吃奶的免费视频| 欧美丝袜亚洲另类| 亚洲国产av影院在线观看| 国产高清三级在线| 激情五月婷婷亚洲| 精品人妻熟女毛片av久久网站| 午夜福利视频精品| 欧美xxⅹ黑人| 国产 精品1| 国产黄频视频在线观看| 美女视频免费永久观看网站| 美女国产视频在线观看| 精品一区二区三卡| 狂野欧美激情性xxxx在线观看| 十八禁网站网址无遮挡| 国产成人a∨麻豆精品| 全区人妻精品视频| 女性被躁到高潮视频| a级毛片在线看网站| 亚洲五月色婷婷综合| 亚洲精品,欧美精品| 飞空精品影院首页| 美女国产视频在线观看| 久久狼人影院| 丰满迷人的少妇在线观看| 精品少妇黑人巨大在线播放| 亚洲高清免费不卡视频| 久久女婷五月综合色啪小说| 久久久久久久精品精品| 国产有黄有色有爽视频| 97精品久久久久久久久久精品| 亚洲美女黄色视频免费看| 中文字幕免费在线视频6| 亚洲一级一片aⅴ在线观看| 秋霞伦理黄片| 久久精品久久精品一区二区三区| 一级毛片黄色毛片免费观看视频| 亚洲人与动物交配视频| 欧美bdsm另类| 久久国产精品男人的天堂亚洲 | 精品人妻熟女av久视频| 日本色播在线视频| 精品人妻一区二区三区麻豆| 黑人猛操日本美女一级片| 一个人看视频在线观看www免费| 久久午夜福利片| 我的老师免费观看完整版| 国产亚洲一区二区精品| 日韩制服骚丝袜av| 青春草视频在线免费观看| 久久午夜综合久久蜜桃| 老司机影院毛片| 国产 精品1| 哪个播放器可以免费观看大片| av黄色大香蕉| 成人影院久久| 久久久久久久久久人人人人人人| a 毛片基地| 观看av在线不卡| 午夜影院在线不卡| 观看美女的网站| 高清午夜精品一区二区三区| 午夜激情福利司机影院| 日韩精品有码人妻一区| 久久精品熟女亚洲av麻豆精品| 日韩成人av中文字幕在线观看| 久久久久久久久久人人人人人人| 免费黄色在线免费观看| 亚洲怡红院男人天堂| 午夜免费观看性视频| 亚洲,一卡二卡三卡| 日本黄大片高清| 国产无遮挡羞羞视频在线观看| 日本vs欧美在线观看视频| 午夜视频国产福利| 久久韩国三级中文字幕| 天天躁夜夜躁狠狠久久av| 在线亚洲精品国产二区图片欧美 | 91精品三级在线观看| 美女cb高潮喷水在线观看| 亚洲精品,欧美精品| av在线app专区| 亚洲欧洲精品一区二区精品久久久 | 日本黄色日本黄色录像| 日韩亚洲欧美综合| av免费观看日本| 日本欧美视频一区| 免费观看性生交大片5| 寂寞人妻少妇视频99o| 亚洲av综合色区一区| 各种免费的搞黄视频| 伦精品一区二区三区| 夜夜骑夜夜射夜夜干| 国产精品成人在线| 亚洲av成人精品一区久久| 亚洲国产精品一区三区| 夜夜看夜夜爽夜夜摸| 热re99久久精品国产66热6| 久久久久久久久大av| 熟女av电影| 午夜91福利影院| 视频在线观看一区二区三区| 日日啪夜夜爽| 搡老乐熟女国产| 熟女电影av网| 热re99久久精品国产66热6| 国产精品国产三级国产av玫瑰| 国产成人一区二区在线| 我要看黄色一级片免费的| 一级,二级,三级黄色视频| 欧美日韩亚洲高清精品| 国产成人91sexporn| 精品一区二区三区视频在线| 中文欧美无线码| 婷婷成人精品国产| 成人国语在线视频| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | 免费av不卡在线播放| 国产黄色视频一区二区在线观看| 2022亚洲国产成人精品| 黑人巨大精品欧美一区二区蜜桃 | 91精品国产国语对白视频| 久久久久久久久久久久大奶| 黑人巨大精品欧美一区二区蜜桃 | 亚洲国产最新在线播放| 午夜视频国产福利| 一级毛片我不卡| 国产精品偷伦视频观看了| 亚洲国产欧美日韩在线播放| 国产亚洲欧美精品永久| 女性生殖器流出的白浆| 国内精品宾馆在线| 欧美bdsm另类| 国产午夜精品一二区理论片| 啦啦啦中文免费视频观看日本| 日本-黄色视频高清免费观看| 亚洲国产欧美日韩在线播放| 精品久久国产蜜桃| 精品人妻熟女毛片av久久网站| 美女国产高潮福利片在线看| 国产一区二区在线观看日韩| 久久久国产欧美日韩av| 啦啦啦中文免费视频观看日本| 欧美精品人与动牲交sv欧美| 嘟嘟电影网在线观看| 精品久久国产蜜桃| 国产不卡av网站在线观看| 高清欧美精品videossex| 男女啪啪激烈高潮av片| 精品99又大又爽又粗少妇毛片| 国产亚洲最大av| 看免费成人av毛片| 少妇人妻精品综合一区二区| 69精品国产乱码久久久| 少妇的逼水好多| 激情五月婷婷亚洲| av.在线天堂| 三上悠亚av全集在线观看| 亚洲精品色激情综合| 亚洲中文av在线| 日韩av免费高清视频| 久久精品国产亚洲av涩爱| 熟女人妻精品中文字幕| 中文欧美无线码| h视频一区二区三区| 欧美性感艳星| 欧美精品一区二区大全| 涩涩av久久男人的天堂| 精品一区在线观看国产| 午夜福利在线观看免费完整高清在| 国产成人aa在线观看| 亚洲国产精品专区欧美| 午夜影院在线不卡| 欧美另类一区| 成年美女黄网站色视频大全免费 | 亚洲精华国产精华液的使用体验| 免费观看无遮挡的男女| 久久99精品国语久久久|