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

    EXL-50U 球形環(huán)中快離子磁場波紋損失的優(yōu)化模擬研究*

    2023-11-24 05:06:12郝保龍李穎穎陳偉郝廣周顧翔孫恬恬王嵎民董家齊袁保山彭元凱石躍江謝華生劉敏勝ENNTEAM
    物理學(xué)報 2023年21期
    關(guān)鍵詞:位形局域波紋

    郝保龍 李穎穎 陳偉 郝廣周 顧翔 孫恬恬王嵎民 董家齊 袁保山 彭元凱 石躍江謝華生 劉敏勝 ENN TEAM

    1) (核工業(yè)西南物理研究院,成都 610041)

    2) (河北省緊湊型聚變重點實驗室,廊坊 065001)

    3) (新奧科技發(fā)展有限公司,廊坊 065001)

    EXL-50U 裝置高參數(shù)等離子體的實現(xiàn)對中性束注入(NBI)加熱的依賴非常敏感,期望NBI 快離子約束良好并通過碰撞慢化把能量傳給背景等離子體.本文基于集成模擬給出的平衡位形、快離子分布和裝置波紋度數(shù)據(jù)對快離子波紋損失開展了模擬研究.發(fā)現(xiàn)快離子波紋損份額約為37%,局域熱斑約 0.6 MW/m2,對裝置實驗運行來說不可接受.其優(yōu)化方案包括移動等離子體位置和加FI(鐵素體鋼插件)降低波紋度,增大 Ip(等離子體電流)以及優(yōu)化 NBI 角度.結(jié)果顯示必須控制波紋度分布且增大 Ip 到 600 kA 以上,才能使快離子損失降低到 3%—4%,局域熱斑降低一個量級.本文總結(jié)了裝置設(shè)計時快離子波紋損失評估的方法,包括相空間快離子分布和波紋損失區(qū)重合度,全要素慢化時間尺度粒子跟蹤.還總結(jié)了降低波紋損失的工程和物理途徑,為集成模擬迭代優(yōu)化和裝置運行提供模擬支持.

    1 引言

    中國磁約束核聚變經(jīng)過半個多世紀(jì)的實驗和物理研究取得了巨大成就,以及不斷刷新參數(shù)記錄的進展.除專業(yè)科研院所外,高校和企業(yè)在理論模擬和裝置建設(shè)運行上也出現(xiàn)了多點開花的有利局面,為我國聚變能開發(fā)利用實現(xiàn)從無到有、從有到優(yōu)打下堅實基礎(chǔ).傳統(tǒng)托卡馬克的大環(huán)徑比為3—4或稍高,EAST 約4.2.小環(huán)徑比(~2 或更低)裝置更緊湊,沒有中心螺線管,有望實現(xiàn)更高的等離子體比壓,即更高的約束效率,位于廊坊新奧科技有限公司的EXL-50U 即屬此類磁約束核聚變裝置,表1 是其主機參數(shù).EXL-50U 設(shè)計有~7 MW輔助加熱和電流驅(qū)動系統(tǒng),包括50 keV/1.5 MW的NBI 和2 MW 的離子回旋波,經(jīng)過前期集成模擬設(shè)計,背景離子溫度可達Ti(0)~40 eV,如圖1所示.EXL-50U 的參數(shù)實現(xiàn)主要依靠NBI 芯部加熱,因此,NBI 快離子的良好約束是實現(xiàn)裝置目標(biāo)的前提.與傳統(tǒng)托卡馬克一樣,分立的有限數(shù)目縱場線圈(TF)引入磁場波紋擾動可引起快離子顯著損失.理論上,任何破壞環(huán)對稱約束平衡位形的電磁擾動都會引起快離子約束變差,損失是否可接受的判據(jù)是加熱功率的損失不影響裝置目標(biāo)實現(xiàn),且損失粒子局域沉積的熱斑不影響裝置安全[1,2].

    圖1 EXL-50U 集成模擬設(shè)計中背景電子密度ne、溫度Te和離子溫度分布TiFig.1.Distribution of bulk electron density ne,electron temperature Te and bulk ion temperature Ti in EXL-50U integrated modeling.

    表1 EXL-50U 與其他托卡馬克裝置主機參數(shù)對比Table 1.Main parameters comparison of EXL-50U and other tokomak facilities.

    磁場波紋擾動引起的快離子損失主要有兩種機制: 波紋磁阱俘獲和隨機波紋擴散.波紋磁阱是裝置局域的波紋擾動較大時,如外中平面附近,形成二級磁阱,平行速度較小的捕獲粒子經(jīng)過該區(qū)域時就環(huán)向俘獲在該局域磁阱內(nèi),通過徑向漂移快速損失掉或經(jīng)碰撞和漂移逃出該損失區(qū),該損失機制和位形、波紋度分布和快離子分布有關(guān),損失的時間尺度在幾個極向回彈周期內(nèi).當(dāng)波紋度較小不足以形成局域磁阱時,有限的波紋擾動也會使捕獲粒子極向投影不閉合,香蕉軌道轉(zhuǎn)折點在平衡位置附近小位移振蕩,當(dāng)香蕉軌道轉(zhuǎn)折點在TF 周期擾動下形成的相空間共振島相互重疊超過隨機判據(jù)時,即發(fā)生連續(xù)周期運動之間解耦,形成隨機波紋擴散,碰撞效應(yīng)可降低該隨機擴散閾值.通行粒子因沿磁力線回旋變換會平均該擾動,因此僅有捕獲粒子有波紋損失.快離子波紋損失是裝置設(shè)計和實驗運行需要開展的必要課題,因為縱場波紋擾動是始終存在的[3,4].

    本文基于EXL-50U 集成模擬設(shè)計參數(shù),利用中性束沉積計算模塊NUBEAM 得到NBI 快離子分布函數(shù)和碰撞率,在哈密頓導(dǎo)心軌道跟蹤程序Orbit 中構(gòu)建平衡疊加波紋擾動后的總磁場位形,一個慢化時間內(nèi)跟蹤蒙卡樣本粒子,統(tǒng)計損失信息.NUBEAM 是20 世紀(jì)70 年代開發(fā)維護至今的經(jīng)典成熟NBI 模塊,可單機或集成在隨時演化輸運程序中運行,可計算NBI 和聚變產(chǎn)物的功率沉積、電流驅(qū)動、動量輸運和加料等.在軸對稱位形中跟蹤蒙卡粒子,可計算多種成分的初始和經(jīng)典穩(wěn)態(tài)慢化分布,也可考慮反常輸運效應(yīng),如Sawtooth,Fishbone 和Ripple loss,但相應(yīng)的模塊粗糙,需要設(shè)定人為經(jīng)驗系數(shù),如Ripple loss 中可以直接指定捕獲粒子的損失時間和隨機擴散閾值,定量分析受限.近年在NUBEAM 中實現(xiàn)的RF-kick,MHDkick 和CGM 等簡化模型在實驗分析中得到不斷發(fā)展和完善[5-8].

    程序Orbit 是20 世紀(jì)80 年代開發(fā)升級至今的經(jīng)典成熟快粒子跟蹤程序,基于直線磁面坐標(biāo)系讀入軸對稱數(shù)值平衡或解析構(gòu)建圓截面平衡,在二維(2D)樣條插值構(gòu)建的背景場中疊加磁場波紋擾動得到總磁場位形,讀入仿星器三維平衡版本為Orbit-3D.在總磁場位形基礎(chǔ)上,Orbit 可選擇讀入 δB=?×αB0或 δB=?×(ξ×B0) 形式的磁流體力學(xué)不穩(wěn)定性(MHD)擾動,其中α是擾動函數(shù)的標(biāo)量形式,ξ是矢量形式的擾動位移,粒子的空間坐標(biāo)為(ψp,θ,ζ),分別是極向磁通、極向角和廣義環(huán)向角.Orbit中粒子導(dǎo)心哈密頓量表達為H=B2/2+μB+Φ,ρ//=v///B是歸一化的平行磁場速度,B是總磁場強度,μ=/2B是粒子磁矩,v⊥是粒子垂直磁場速度,Φ是靜電勢.粒子的速度由能量E、磁矩μ和螺距角正負號決定,Orbit 程序求解跟蹤粒子哈密頓導(dǎo)心軌道時考慮碰撞修正,導(dǎo)心方程中表達磁場擾動的是總磁場強度的空間導(dǎo)數(shù) (?B/?ζ,?B/?θ,?B/ψp),而無需矢量分量,Orbit 中根據(jù)磁場散度為零的條件和不同函數(shù)形式構(gòu)建擾動場[9-11].

    本文第2 節(jié)根據(jù)EXL-50U 裝置的波紋場工程數(shù)據(jù),在Orbit 程序中數(shù)值重建,判定波紋損失區(qū),分析相空間中NBI 快離子初始分布函數(shù)和波紋損失區(qū)的重合度.第3 節(jié)全要素計算,即在總磁場中跟蹤大量樣本粒子分布一個慢化時間,統(tǒng)計粒子信息,針對波紋損失計算結(jié)果,提出并論證多個優(yōu)化方向,對EXL-50U 設(shè)計建設(shè)提供數(shù)值支持.第4節(jié)總結(jié)本文工作并討論一般托卡馬克裝置波紋損失評估流程.

    2 EXL-50U 裝置波紋場數(shù)值重建和NBI 快離子分布函數(shù)

    2.1 縱場波紋擾動重建

    EXL-50U 環(huán)向磁場由12 柄瘦D 形線圈實現(xiàn),對比其他托卡馬克裝置,柄數(shù)較少,有較大磁場波紋擾動.定義在極向截面的波紋度為δ(r,θ)=(Bmax-Bmin)/(Bmax+Bmin),Bmax和Bmin分別是同一極向位置不同環(huán)向角處的磁場最大值和最小值.工程設(shè)計中沒考慮包層鐵素體鋼材料,僅有縱場線圈時的波紋度分布見圖2(b),波紋度在等離子體外中平面Rmax~1.4 m 附近最大,約3%,波紋擾動幅度從芯部到邊界呈自然指數(shù)增大.托卡馬克縱場磁體設(shè)計的一般要求是TF 線圈在16—24 柄,等離子體約束區(qū)波紋度最大值即外中平面處小于0.5%,EXL-50U 由于特殊的主機和等離子體尺寸設(shè)計不滿足此要求.

    圖2 EXL-50U 縱場波紋度分布 (a) 解析函數(shù)實現(xiàn)值;(b) 工程設(shè)計值Fig.2.Distribution of toroidal field ripple perturbation amplitude in EXL-50U: (a) Ripple data by analytical equation;(b) engineering data in design.

    在Orbit 程序中引入波紋擾動時,在磁面坐標(biāo)系中疊加縱場波紋后的總磁場可表達為B(ψp,θ,φ)=B0(ψp,θ)[1+δcos(Nφ)],其中B0(ψp,θ) 是平衡場.本文采用解析公式δ(R,Z)=δ0exp{[(R-Rrip)2+bripZ2]0.5/wrip}來擬合工程設(shè)計值,該解析函數(shù)的初始形式是針對TFTR 類似圓截面TF 線圈產(chǎn)生的一系列波紋度同心標(biāo)準(zhǔn)圓等高線提出的,后在ITER 等裝置的D 形線圈波紋度幾何分布發(fā)展為包含橢圓形變和等高線圓心位移的現(xiàn)在形式,在一系列裝置上應(yīng)用成熟,包括EAST 和CFETR[12-14].基于EXL-50U 裝置波紋度工程數(shù)據(jù),數(shù)學(xué)含義和擬合系數(shù)分別為波紋度最小值δ0=4.19×10-8,波紋度不同值的等高線幾何圓心位置Rrip=1.77×10-3+0.106Z2m,橢圓率brip=0.297 和波紋擾動特征空間尺度wrip=0.103 m.波紋擾動的解析函數(shù)實現(xiàn)和工程設(shè)計值對比見圖2—4.由圖2 可知解析函數(shù)重建的波紋度幾何輪廓和幅值較大精度還原了工程設(shè)計;由圖3 和圖4 的不同縱橫截面對比可知解析函數(shù)重建的波紋場與工程數(shù)據(jù)有一定差距,圖3 中不同Z截面的波紋度基本不變,這是EXL-50U 特殊的TF 線圈引起的,即瘦長D 形而非傳統(tǒng)托卡馬克的胖D 形或圓形.本文解析函數(shù)重建波紋場的整體符合程度稍差于傳統(tǒng)托卡馬克[3,15],但仍精確重建了等離子體約束區(qū)的波紋場,即LCFS (最后一個閉合磁面)內(nèi)的擬合精確度足夠高,10-5— 10-4,鑒于該函數(shù)形式有清晰的幾何含義,如表2 所列,優(yōu)于高階多項式擬合,本文工作沿用.

    圖3 EXL-50U 裝置波紋度不同Z 截面工程數(shù)據(jù)和數(shù)值擬合結(jié)果對比 (a) Z=0 m;(b) Z=0.3 m;(c) Z=0.6 m;(d) Z=0.9 m Fig.3.Ripple comparison between engineering design and fitting curve in different Z plane of EXL-50U: (a) Z=0 m;(b) Z=0.3 m;(c) Z=0.6 m;(d) Z=0.9 m.

    圖4 EXL-50U 裝置波紋度不同R 截面工程數(shù)據(jù)和數(shù)值擬合結(jié)果對比 (a) R=0.5 m;(b) R=0.7 m;(c) R=0.9 m;(d) R=1.1 mFig.4.Ripple comparison between engineering design and fitting curve in different R plane of EXL-50U: (a) R=0.5 m;(b) R=0.7 m;(c) R=0.9 m;(d) R=1.1 m.

    表2 EXL-50U 與其他托卡馬克裝置縱場波紋數(shù)據(jù)擬合結(jié)果對比Table 2.Ripple field fitting parameters comparison of EXL-50U and other tokomak facilities.

    2.2 NBI 快離子分布函數(shù)

    利用NUBEAM 程序計算得到的快離子初始分布函數(shù)一般用于長時間尺度的粒子跟蹤損失,包含相空間梯度信息的穩(wěn)態(tài)慢化分布用于分析波粒非線性相互作用.本文基于EXL-50U 集成模擬中的平衡和背景參數(shù)剖面計算得到的NBI 初始分布函數(shù)如圖5 和圖6 所示,可以看出,NBI 有良好的芯部沉積和微弱的穿透損失.不同于傳統(tǒng)超導(dǎo)托卡馬克,EXL-50U 的NBI 束線路徑不經(jīng)過高場側(cè)靶板,為標(biāo)準(zhǔn)的切向注入,快離子pitch(螺距)角較大,由于整體安全因子剖面較高,平衡位形中快離子分布有捕獲粒子份額約30%.NBI 設(shè)計運行束能量在20—50 keV,進窗口功率在0.25—1.5 MW,NBI 系統(tǒng)建成后的束能和束流強度對應(yīng)關(guān)系類似伏安特性曲線.本文評估快離子損失的初始能量設(shè)定為45 keV/1.0 MW,平衡位形的縱場方向都是俯視逆時針,等離子體電流俯視順時針.

    圖5 EXL-50U 上NBI 快離子初始分布的粒子空間位置俯視圖(a)和極向投影(b)Fig.5.Initial distribution of beam ions in EXL-50U: (a) Bird's view;(b) poloidal cross section.

    圖6 Orbit 程序中讀入的NBI 快離子初始分布函數(shù) (a) RZ 空間分布;(b)粒子密度的極向磁通分布;(c)能量分布;(d) pitch 角分布Fig.6.Initial distribution of NBI fast ions read by Orbit code: (a) Particle location in RZ coordinate;(b) particle density distribution in poloidal flux;(c) energy distribution;(d) pitch angle distribution.

    2.3 波紋損失區(qū)

    在裝置設(shè)計初期,工程和物理設(shè)計迭代速度快,無法通過耗時的慢化時間尺度內(nèi)全要素粒子跟蹤進行快離子損失計算.為直觀高效討論EXL-50U 上快離子波紋損失,本節(jié)采用相空間損失區(qū)域和粒子分布重合度來評估.圖7 是在平衡位形疊加磁場波紋后的總磁場位形中刻畫得到的波紋損失區(qū),其中波紋磁阱俘獲區(qū)由磁阱形成條件ε|sinθ|<qNδ確定,無碰撞波紋隨機擴散區(qū)的確定依據(jù)是簡化判據(jù)GWB(Goldston-White-Boozer)δs=(ε/Nπq)1.5/ρLq',其中ε是當(dāng)?shù)氐目v橫比倒數(shù),ρL是拉莫爾回旋半徑,q'=dq/dr是安全因子導(dǎo)數(shù),快離子分布相空間對應(yīng)的波紋幅度大于GWB 判據(jù)即認為捕獲粒子轉(zhuǎn)折點在此被隨機擴散損失掉,因此波紋隨機擴散損失和快離子能量、pitch 角、安全因子剖面等有關(guān)系.由于EXL-50U 裝置的整體波紋度較大,且等離子體約束區(qū)位置靠外,波紋損失區(qū)基本占滿空間,意味著捕獲粒子基本都會損失[16].

    初始分布的NBI 快離子能量可認為是單值分布,即決定粒子軌道的3 個物理量中固定E,可以得到一個二維平面(Pζ,μB0/E),其中B0是磁軸處磁場強度,為Orbit 程序中的磁場強度單位;Pζ是正則環(huán)向角動量.在軸對稱位形和能量守恒下可以在(Pζ,μB0/E)平面內(nèi)得到由一系列曲線分割開的不同軌道類型區(qū)域,如圖8(a),其中 T-C 指約束捕獲粒子軌道,T-L 指損失捕獲粒子軌道,P±-L指損失通行粒子軌道,P±-C 指約束通行粒子軌道,下標(biāo)的正負號指螺距角符號[10,17];橫軸題Pζ/ψw是被LCFS 處環(huán)向磁通歸一化的正則環(huán)向角動量.該平面中還有一些特殊的軌道沒有標(biāo)記出,如土豆軌道和極向靜止軌道,不同平衡位形和快離子分布下的軌道類型區(qū)域占比不同,但拓撲結(jié)構(gòu)不變.圖8(b)為波紋損失區(qū)的分布,同圖7 類似,捕獲粒子空間幾乎被波紋損失區(qū)占滿.捕獲粒子的香蕉軌道轉(zhuǎn)折點落入波紋損失區(qū)即表示快離子會損失,快離子損失可能會反復(fù)經(jīng)歷兩個損失通道.隨機碰撞會改變快離子運動的螺距角,可能讓快離子提前進入或者逃出波紋損失區(qū),因此下文通過判斷快離子是否位于波紋損失區(qū)來判斷粒子損失與否時,忽略了隨機碰撞,一般來說現(xiàn)有裝置參數(shù)水平下,碰撞波紋隨機擴散是主要的快離子波紋損失通道;聚變堆參數(shù)下,無碰撞波紋隨機擴散是主要的損失通道,因此圖8 中只用GWB 判據(jù),即無碰撞波紋隨機擴散標(biāo)記的波紋損失區(qū)來評估和快離子密度分布的重合情況.

    圖8(c)和圖8(d)是原始設(shè)計下兩個離子源初始快離子的相空間分布.由圖8 可知,NBI 沉積芯部峰化,因兩個離子源注入角度不同,快離子分布中捕獲粒子份額不同,co-perp (較垂直)源的捕獲粒子份額較多,初始軌道損失份額較多.根據(jù)相空間的波紋損失區(qū)和快離子分布重疊情況,可以判斷目前的設(shè)計參數(shù)下,NBI 快離子的波紋損失較顯著,主要原因是波紋損失區(qū)太大,且快離子分布中捕獲粒子份額太大,為后續(xù)定量的全要素計算結(jié)果和優(yōu)化方向提供清晰物理圖像的支持.

    3 NBI 快離子初始分布慢化時間尺度的導(dǎo)心軌道跟蹤

    3.1 基于集成模擬設(shè)計中平衡和快離子分布的計算結(jié)果

    本節(jié)計算輸入即基于前述集成模擬平衡位形和快離子分布函數(shù),在一個慢化時間內(nèi)跟蹤波紋擾動疊加后的總磁場位形中NBI 快離子初始分布樣本,碰撞頻率分布由NBUEAM 程序基于背景等離子體參數(shù)計算得到,計算取徑向中間位置為典型值,初始能量的能量慢化率νε0=0.015 s-1,螺距角散射率為νpa0=0.05 s-1.一個計算步長 Δt后,Orbit 程序?qū)αW蛹磿r能量Ei和螺距角λi=v///v的修正分別為

    其中T(eV) 是等離子體溫度[18].程序中對NBI 快離子初始分布跟蹤一個慢化時間后的分布信息如圖9 所示,對比初始分布圖6 可知,粒子能量在一個慢化時間內(nèi)從初始能量慢化到接近背景溫度,pitch 角分布從特殊的典型集中值到基本各向同性,符合物理實際.

    考慮縱場波紋擾動和碰撞效應(yīng)的粒子損失份額如表3 所列,其中捕獲粒子份額為在平衡位形疊加波紋場后的總磁場位形中得到,損失粒子份額基本接近捕獲粒子份額,印證了前述分析中波紋損失區(qū)基本占滿捕獲粒子相空間的結(jié)論.損失粒子信息見圖10,縱場俯視順時針,損失粒子因梯度漂移和曲率漂移方向集中在外中平面附近偏下,初始軌道損失后的損失份額隨著碰撞和波紋擾動的時間累計逐漸增加.Orbit 程序定義LCFS 為損失邊界,由于裝置空間有限且射頻天線等第一壁部件需要盡可能貼合LCFS 以獲得良好的功率耦合,因此LCFS 與第一壁形貌和距離接近,因此可以根據(jù)該處的粒子沉積位置、數(shù)量和能量信息評估局域熱負荷.圖11 是EXL-50U 集成模擬參數(shù)下,假設(shè)coperp 離子源對應(yīng)0.5 MW 的快離子損失在LCFS上形成的局域熱斑,峰值約0.2 MW/m2,一共有12 個熱斑,和裝置TF 線圈數(shù)目對應(yīng),位置位于相鄰兩柄線圈中間.考慮到co-tang 離子源損失的粒子沉積在相近位置,熱負荷疊加,總的熱斑峰值約0.4 MW/m2,因此對裝置運行安全十分不利.NBI快離子損失份額過大意味著輔助加熱效率低,且考慮局域熱斑危害需要進一步優(yōu)化快離子約束.按一般研究經(jīng)驗,假設(shè)EXL-50U 第一壁部件的熱負荷安全閾值在0.5 MW/m2,考慮到NBI 多個離子源的波紋損失局域熱負荷疊加,粗估EXL-50U 的NBI 快離子波紋損失份額上限為15%作為優(yōu)化計算的參照.

    圖10 Co-perp 快離子分布下的(a)損失粒子極向位置分布,(b)損失時間和能量記錄,(c)損失份額的隨時演化Fig.10.(a) Poloidal distribution of lost particle,(b) lost time and energy record,(c) time evolution of loss fraction for co-perp beam ion distribution.

    圖11 一個慢化時間后損失的NBI 快離子局域沉積在LCFS 處得到的熱負荷Fig.11.Heat load at the last closed flux surface due to NBI fast ions loss after a slowing down time.

    表3 基于集成模擬平衡位形和快離子分布的波紋損失全要素計算結(jié)果Table 3.Ripple loss results of full calculation based on integrated modeling equilibrium and beam ion distribution.

    3.2 降低波紋度: 移動等離子體位置和加FI(鐵素體鋼插件)

    前述計算的快離子份額過大的顯著原因是裝置波紋度較大,且NBI 快離子的捕獲粒子份額較大.托卡馬克裝置的波紋度是從芯部到邊界e 指數(shù)增大,由圖2 可知,等離子體約束區(qū)位于波紋度輪廓中心的偏低場側(cè),往高場側(cè)移動等離子體即可降低波紋度分布,即減小LCFS 的最大徑向位置Rmax.圖12 是LCFS的Rmax移動到1.32 m 時的波紋損失區(qū),對比圖7 和圖8中Rmax在1.35 m 的分布可知,LCFS 移動3 cm,等離子體約束區(qū)的波紋度就有顯著下降.為定量對比移動LCFS的Rmax對NBI快離子波紋損失影響,保持等離子體電流Ip=500 kA 固定,僅移動Rmax的計算結(jié)果如表4 所列.為方便對比,表中計算輸入的快離子分布來自2.2 節(jié)平衡位形Rmax~1.35 m 兩個離子源同時注入,快離子在不同LCFS 位置下的RZ平面分布不同,表4 僅用于定量對比說明LCFS 移動可顯著降低波紋度和波紋損失.

    表4 EXL-50U 中固定 Ip=500 kA,不同Rmax時波紋損失計算結(jié)果Table 4.Ripple loss results of different Rmax of LCFS equilibrium with Ip=500 kA in EXL-50U.

    工程上解決托卡馬克TF 線圈尺寸和柄數(shù)有限時波紋度較大的主要辦法是引入FI,在模擬和實驗上證實效果顯著[19].EXL-50U 引入FI 后可顯著降低裝置波紋度,以波紋度最大的外中平面處Rmax~1.226 m 為例,波紋度可從0.63%降低到0.34%.工程上FI 降低波紋度的效果取決于所加材料的多少,聚變堆由于空間有限,FI 效果有限.EXL-50U 的FI 材料有充分設(shè)計空間,效果顯著,圖13 是固定平衡位形下加FI 前后的波紋損失區(qū)對比.

    3.3 優(yōu)化NBI 注入角度和束能

    EXL-50U 初始設(shè)計中NBI 的兩個離子源注入的切向半徑Rtang分別為0.428 m 和0.607 m,在赤道面水平注入,如圖5 所示,快離子分布中捕獲粒子份額較大.本節(jié)優(yōu)化NBI 注入角度包括在裝置窗口空間允許范圍內(nèi)掃描切向半徑Rtang,保證NBI 束線穿過窗口時法蘭中心不動且沿垂直方向移動離子源位置Zelev,注入布局改變即改變NBI角度,優(yōu)化的前提條件是保證NBI 功率主要沉積在芯部.通常,NBI 角度小更容易加熱芯部,尤其是高密度聚變堆在有限束能工程能力限制下,NBI不能過于切向,否則功率無法芯部沉積.同時,NBI角度小引入的捕獲粒子份額多,不利于快離子約束和慢化加熱,需要數(shù)值模擬找到最優(yōu)方案.表5 是利用NUBEAM 程序掃描不同NBI 角度時的捕獲粒子份額,可知垂直方向移動離子源基本不改變捕獲粒子份額,其與NBI 整體束線長度和等離子體尺寸有關(guān).改變NBI 幾何布局會影響沉積剖面,一般情況下NBI 芯部沉積優(yōu)先級最高,本文優(yōu)化模擬中對離軸沉積算例舍棄.如圖14 所示,切向半徑在Rtang~1 m 附近時就有NBI 沉積大幅度偏離芯部,因此物理上最優(yōu)的NBI 注入角度在Rtang~0.8 m 附近.

    圖14 不同NBI 注入角度下的NBI 沉積路徑的極向截面 (a)切向半徑Rtang=0.8 m,離子源垂直中平面抬升Zelev=0 m;(b)切向半徑Rtang=1.0 m,離子源垂直中平面抬升Zelev=0.6 mFig.14.Cross section of NBI deposition trajectory with different NBI geometry: (a) Beam tangency radius Rtang=0.8 m,elevation of beam ion source above midplane Zelev=0 m;(b) beam tangency radius Rtang=1.0 m,elevation of beam ion source above midplane Zelev=0.6 m.

    表5 不同NBI 角度下的快離子分布中的捕獲粒子份額Table 5.Trapped particle faction of beam ions with different NBI geometry.

    為了探索優(yōu)化波紋度后的NBI 合適運行區(qū)間,表6 列出了Ip=500 kA 時不同NBI 切向半徑和束能下的快離子損失結(jié)果,其中NBI 切向半徑Rtang=0.428 m,0.607 m 和0.807 m 時的捕獲粒子份額分別為0.514,0.433 和0.352.快離子損失份額隨著束能增大而增大,物理圖像是快離子能量增大,回旋半徑增大,進而粒子運動的軌道寬度增大,初始軌道損失和波紋損失顯著依賴軌道寬度.

    表6 平衡位形LCFS Rmax~1.16 m,Ip=500 kA 時不同NBI 角度和束能Enb 下的快離子損失份額Table 6.Loss faction of NBI fast ions with different NBI geometry and beam energy in LCFS Rmax~1.16 m,Ip=500 kA equilibrium.

    基于表6 結(jié)果可知,原EXL-50U 設(shè)計中的NBI 兩個離子源的切向半徑僅有Rtang=0.607 m的co-tang 源的損失份額在允許范圍內(nèi).Rtang=0.428 m 的離子源損失份額過大,降低束能意味著降低束功率,雖然可使熱負荷降低,但集成模擬研究顯示裝置的設(shè)計目標(biāo)敏感依賴NBI 加熱功率,降低束能得不償失.增大NBI 切向半徑到Rtang=0.807 m 后可以允許滿功率,但對裝置窗口條件要求高,且快離子損失份額仍有優(yōu)化空間,需要探索其他方向,在工程和物理上折中選擇,避免某一個方向極端.

    3.4 提高運行電流Ip 降低快離子損失

    快離子初始軌道損失來源于快離子初始位置太靠近損失邊界,且軌道寬度過大,以至于無法完成完整的極向運動周期即損失掉,一般參數(shù)下的捕獲粒子軌道寬度比通行粒子大幾倍甚至一個量級,降低初始軌道損失需要使快離子分布剖面盡可能聚芯,遠離損失邊界且捕獲粒子成分少.捕獲粒子的軌道寬度近似為Δwidth(8/ε),其中qρL為極向磁場強度對應(yīng)的粒子回旋半徑,和電流Ip反比.降低快離子波紋損失的途徑類似,即使等離子體盡可能處在低波紋度區(qū)域,快離子分布剖面盡可能聚芯,減少捕獲粒子份額.現(xiàn)在托卡馬克裝置設(shè)計之初即盡可能降低波紋度,波紋磁阱區(qū)域一般很小.波紋隨機擴散損失中,捕獲粒子軌道轉(zhuǎn)折點的擾動位移近似為 Δr(Nπ/|sinθb|)1/2(q/ε)2/3ρLδcos(NΔφb),其中θb是轉(zhuǎn)折點的極向角,Δφb是相鄰轉(zhuǎn)折點的環(huán)向角度差,該表達形式同樣敏感依賴軌道寬度.

    初始軌道損失和波紋損失的物理圖像顯示,在工程優(yōu)化前提下,提高等離子體電流是效果最顯著的方法.表7 計算結(jié)果是基于LCFS 都在Rmax<1.18 m,稍有不同的是,不同Ip下對應(yīng)的平衡位形不同的,在TF+FI 波紋場中,采用同一NUBEAM快離子分布計算結(jié)果,包含原始設(shè)計NBI 切向半徑的45 keV 兩個離子源.表7 中同一快離子分布的捕獲粒子份額沒有明顯規(guī)律,是因為NBI 束線和等離子體相對空間位置和安全因子q剖面沒有單調(diào)規(guī)律,而非純解析對比捕獲粒子份額對損失的影響,只是用來說明Ip提高可顯著降低損失,對比波紋損失區(qū)也會發(fā)現(xiàn)隨著Ip增大其明顯減少.Ip在600 kA 以上的快離子損失份額在允許范圍,以Ip=700 kA 的平衡位形為例,單個NBI 離子源快離子分布下的損失份額在Rtang=0.607 m 時為5.3%,Rtang=0.707 m 時為4.3%,局域熱負荷比原始設(shè)計低一個量級,損失份額滿足加熱效率和第一壁熱負荷安全.

    表7 不同Ip 下的NBI 快離子損失粒子份額Table 7.Trapped particle faction of NBI fast ions with different Ip.

    4 總結(jié)和討論

    本文基于EXL-50U 集成模擬設(shè)計中的平衡位形、NBI 快離子分布和裝置波紋度,利用導(dǎo)心程序Orbit 在一個慢化時間內(nèi)跟蹤粒子損失份額,并評估局域熱負荷,結(jié)果顯示原始設(shè)計損失過大.優(yōu)化措施包括往高場側(cè)移動LCFS 和加FI 降低波紋度,增大Ip,優(yōu)化NBI 角度,計算結(jié)果證實優(yōu)化措施效果顯著,同時在多種優(yōu)化方法下的快離子份額可降到4.3%,局域熱斑降低一個量級.本文給出的局域熱斑評估方法有幾點限制,即定義損失邊界為LCFS 會高估損失,忽略有限回旋半徑會低估損失,現(xiàn)有裝置的LCFS 和第一壁形貌的貼合度不如聚變堆,如果有第一壁部件突出、形變,會有更局域的熱斑.此外,第一壁部件的熱負荷引起的溫升和等離子體運行時間有關(guān),EXL-50U 的運行時間約1 s,因此同樣的材料能承受的熱負荷比長脈沖或穩(wěn)態(tài)運行時高,本文假定的熱負荷安全閾值完全能確保裝置運行時的安全.

    優(yōu)化快離子損失采取的措施需要在工程和物理上統(tǒng)籌折中、反復(fù)迭代,本文計算僅限于固定平衡和背景參數(shù)的粒子損失,更自洽的模擬需要多次集成模擬迭代,本文結(jié)論為更進一步設(shè)計和實驗運行提供參考.數(shù)值模擬和放電實驗一樣,都不是一勞永逸,能一次性從物理邏輯和數(shù)值程序都自洽給出結(jié)果,需要多個小方向和集成模擬之間的互相迭代改進.

    快離子損失的另一個主要損失通道是MHD不穩(wěn)定性引起的損失和再分布,類似本文通過相空間中快離子密度和損失區(qū)的重合度來評估快離子輸運,MHD 擾動模式信息也可以用來在相空間二維平面(Pζ,E)中刻畫隨機區(qū)域,通過對比快離子在該平面的密度分布評估是否損失或者局域再分布.快離子KAM 面破缺位置位于(Pζ,E)邊界,且有ωnPζ-nE=constant 斜線連接快離子主分布區(qū)和損失邊界,即表示MHD 會引起快離子直接損失.如果KAM 破缺標(biāo)記的隨機區(qū)域位于芯部,即表示有局域的再分布.

    聚變裝置中快離子損失是否在合理范圍內(nèi)主要有兩大判據(jù),第一是損失的輔助加熱或聚變自加熱功率是否影響背景參數(shù)性能的實現(xiàn);第二是損失的快離子沉積在第一壁部件上,是否引起的局域熱斑損失工程部件,如射頻波天線,或者雜質(zhì)問題.聚變界對快離子損失份額有個經(jīng)典的說法,即ITER概念設(shè)計在21 世紀(jì)10 年代指出α 粒子損失份額超過5%即對第一壁造成破壞性損傷,超過25%無法完成聚變點火,第一壁熱斑的承受閾值在0.5 MW/m2,背景輻射約占0.2 MW/m2,允許的快離子損失沉積熱斑約0.3 MW/m2,該說法一直被廣泛采用,出處見參考文獻 [20,21].ITER 第一壁可承受熱斑閾值0.5 MW/m2在文獻 [22]中也再次明確,說明ITER 更新設(shè)計中保留了該閾值.考慮到EXL-50U 非穩(wěn)態(tài)運行和工程實際,第一壁熱斑安全閾值假定0.5 MW/m2合理.國內(nèi)聚變堆設(shè)計中的第一壁熱負荷閾值約1.0 MW/m2,具體數(shù)值還要根據(jù)設(shè)計更新組織專家討論.值得指出的是,聚變堆聚焦的是α 粒子損失份額,由于α 粒子分布天然比NBI 快離子分布剖面聚芯,裝置波紋度是從芯部到邊界e 指數(shù)增大,NBI 從低場側(cè)邊界開始沉積,因此NBI 快離子波紋損失份額一般要遠比α 粒子大.因此,除了第一壁局域熱斑閾值可參考,NBI 快離子損失份額不能簡單限定在5%,需要根據(jù)具體裝置設(shè)計計算.簡而言之,判定NBI 快離子損失份額是否可接受,一是看損失的NBI 功率是否影響運行目標(biāo)實現(xiàn),二是看損失粒子是否打壞部件.通常,降低快離子波紋損失有以下辦法.

    1) 優(yōu)化降低等離子體區(qū)域的裝置波紋度,如增加TF 線圈柄數(shù)、加鐵素體鋼插件.往高場側(cè)移動等離子體,第一壁部件需要跟隨設(shè)計.一般來說FI 效果取決于所加材料量的多少,在聚變堆設(shè)計階段,因第一壁空間限制,FI 效果有限,EXL-50U第一壁空間充足,FI 降低波紋度的效果較好.

    2) 保證NBI 功率芯部沉積的前提下優(yōu)化NBI角度,一般較切向注入有更少的捕獲粒子份額,但需要窗口空間,超導(dǎo)裝置的窗口狹長,優(yōu)化NBI 角度的空間有限,需要主機參數(shù)設(shè)計初期即完成物理和工程的快速迭代,盡可能在集成模擬階段發(fā)揮NBI 加熱能力.

    3) 優(yōu)化背景等離子體參數(shù),使電子密度ne聚芯、峰化分布,以使NBI 沉積聚芯,ne剖面基本決定了NBI 的沉積剖面.快離子越靠近芯部,一般約束越好,此外NBI 離軸沉積也有電流驅(qū)動、輸運壘和抑制MHD 不穩(wěn)定性等方面的可能,需要滿足基本加熱目前的情況下綜合考慮.

    4) 增大ne以縮短慢化時間,減少因碰撞和時間累積產(chǎn)生的粒子損失,以使快離子更快的慢化.增大裝置縱場,減小快離子回旋半徑.降低束能以降低慢化時間和波紋擾動步長,同時意味著降低功率,限制參數(shù)提高.

    5) 增大等離子體運行電流.在工程設(shè)計基本確定時,增大Ip是效果最顯著的優(yōu)化方案,可有效降低快離子軌道寬度和波紋擾動步長,但需要更高的電流驅(qū)動能力.在主機參數(shù)設(shè)計時,增加中心螺線管伏秒數(shù)比增加昂貴、復(fù)雜的輔助系統(tǒng)能力收益大,是托卡馬克類裝置作為電物理設(shè)備產(chǎn)生感應(yīng)電流的核心.

    猜你喜歡
    位形局域波紋
    中間支撐剛度對雙跨梁屈曲穩(wěn)定性的影響
    基于NACA0030的波紋狀翼型氣動特性探索
    小波紋的童話
    局域積分散列最近鄰查找算法
    電子測試(2018年18期)2018-11-14 02:30:34
    為什么水面波紋蕩漾
    學(xué)與玩(2017年5期)2017-02-16 07:06:26
    基于旋量理論的四自由度抓取機械手奇異位形分析
    PET成像的高分辨率快速局域重建算法的建立
    基于可操作度的機器人最優(yōu)初始位形研究
    大眾科技(2015年11期)2015-11-24 01:57:16
    基于最優(yōu)初始位形的冗余度機器人可操作度優(yōu)化*
    基于局域波法和LSSVM的短期負荷預(yù)測
    電測與儀表(2015年7期)2015-04-09 11:39:50
    国产视频一区二区在线看| 最近视频中文字幕2019在线8| 最新中文字幕久久久久 | 免费大片18禁| 夜夜爽天天搞| 精品熟女少妇八av免费久了| 日韩欧美国产在线观看| 成人国产综合亚洲| 精品久久久久久,| 噜噜噜噜噜久久久久久91| 18禁裸乳无遮挡免费网站照片| 别揉我奶头~嗯~啊~动态视频| 午夜免费激情av| 国产精品久久久av美女十八| 黑人巨大精品欧美一区二区mp4| 色在线成人网| 久久婷婷人人爽人人干人人爱| 丰满人妻一区二区三区视频av | 91av网一区二区| 长腿黑丝高跟| 亚洲国产欧美一区二区综合| 成人av在线播放网站| 中文字幕精品亚洲无线码一区| 色av中文字幕| 成人亚洲精品av一区二区| 中文在线观看免费www的网站| 国产人伦9x9x在线观看| 精品久久久久久久久久久久久| 啦啦啦观看免费观看视频高清| 中文字幕高清在线视频| 色老头精品视频在线观看| 成年版毛片免费区| 国产亚洲av高清不卡| 12—13女人毛片做爰片一| 日韩成人在线观看一区二区三区| 一进一出好大好爽视频| 欧美性猛交黑人性爽| 午夜亚洲福利在线播放| 在线观看午夜福利视频| 午夜激情福利司机影院| 熟女少妇亚洲综合色aaa.| 精品久久久久久久久久免费视频| 欧美日韩国产亚洲二区| 色哟哟哟哟哟哟| 午夜福利欧美成人| 国产三级黄色录像| 色播亚洲综合网| 老司机午夜福利在线观看视频| 好看av亚洲va欧美ⅴa在| 99国产极品粉嫩在线观看| 精品久久久久久久毛片微露脸| 在线观看免费午夜福利视频| 真人做人爱边吃奶动态| 女人高潮潮喷娇喘18禁视频| 精品免费久久久久久久清纯| 国产v大片淫在线免费观看| 国产亚洲精品一区二区www| 女人高潮潮喷娇喘18禁视频| 亚洲欧美日韩卡通动漫| 日本在线视频免费播放| 精品国产超薄肉色丝袜足j| 国产精品,欧美在线| 日韩高清综合在线| 久久午夜综合久久蜜桃| 日本黄大片高清| 国产精品一区二区三区四区久久| 99在线视频只有这里精品首页| 国产成人av激情在线播放| 动漫黄色视频在线观看| 女生性感内裤真人,穿戴方法视频| 成人性生交大片免费视频hd| 在线国产一区二区在线| 亚洲美女黄片视频| 亚洲精华国产精华精| 国产久久久一区二区三区| 白带黄色成豆腐渣| 99热只有精品国产| bbb黄色大片| 中文字幕av在线有码专区| 久久精品人妻少妇| 免费看美女性在线毛片视频| 校园春色视频在线观看| 一级黄色大片毛片| 天堂√8在线中文| 久99久视频精品免费| 国产精品野战在线观看| 久久人人精品亚洲av| 偷拍熟女少妇极品色| 亚洲一区高清亚洲精品| 亚洲精品一卡2卡三卡4卡5卡| 久久天堂一区二区三区四区| 在线观看日韩欧美| 国产精品99久久久久久久久| 男女做爰动态图高潮gif福利片| 亚洲国产欧美一区二区综合| 亚洲aⅴ乱码一区二区在线播放| а√天堂www在线а√下载| 亚洲欧美日韩高清专用| 校园春色视频在线观看| 欧美乱码精品一区二区三区| avwww免费| 在线观看午夜福利视频| 国产伦精品一区二区三区四那| 日韩人妻高清精品专区| 男女之事视频高清在线观看| 丰满人妻熟妇乱又伦精品不卡| 国产亚洲精品久久久久久毛片| 精品乱码久久久久久99久播| 中文字幕精品亚洲无线码一区| 国产午夜精品论理片| 国产成人系列免费观看| 日本一二三区视频观看| 国产精品电影一区二区三区| 搡老岳熟女国产| 日韩av在线大香蕉| 亚洲在线观看片| 色哟哟哟哟哟哟| 九九热线精品视视频播放| 两人在一起打扑克的视频| 亚洲专区字幕在线| 国产精品一及| 一a级毛片在线观看| 国产高清视频在线观看网站| 日本撒尿小便嘘嘘汇集6| 国产免费男女视频| 久久中文看片网| 一级作爱视频免费观看| 在线观看舔阴道视频| 啦啦啦免费观看视频1| 国产精品久久久久久人妻精品电影| 午夜成年电影在线免费观看| 精品一区二区三区视频在线观看免费| 亚洲精品粉嫩美女一区| 99久久综合精品五月天人人| 变态另类成人亚洲欧美熟女| 琪琪午夜伦伦电影理论片6080| 国产视频内射| 日韩精品中文字幕看吧| 亚洲熟妇中文字幕五十中出| 三级毛片av免费| 亚洲va日本ⅴa欧美va伊人久久| 久久精品aⅴ一区二区三区四区| 亚洲五月婷婷丁香| 国产伦精品一区二区三区四那| 亚洲自拍偷在线| 亚洲黑人精品在线| 国产av麻豆久久久久久久| av天堂在线播放| 最近最新中文字幕大全免费视频| 男女那种视频在线观看| 村上凉子中文字幕在线| 麻豆av在线久日| 怎么达到女性高潮| 国产一区在线观看成人免费| av女优亚洲男人天堂 | av视频在线观看入口| 国产av一区在线观看免费| 12—13女人毛片做爰片一| 免费看美女性在线毛片视频| 亚洲第一电影网av| 人人妻,人人澡人人爽秒播| 欧美在线一区亚洲| 可以在线观看毛片的网站| 2021天堂中文幕一二区在线观| 亚洲av熟女| 91老司机精品| 18禁黄网站禁片免费观看直播| 国模一区二区三区四区视频 | 欧美日韩乱码在线| 亚洲欧洲精品一区二区精品久久久| 身体一侧抽搐| 亚洲九九香蕉| 天堂影院成人在线观看| 免费看a级黄色片| 久久天躁狠狠躁夜夜2o2o| 99热这里只有精品一区 | 九色成人免费人妻av| 免费看十八禁软件| 黄色丝袜av网址大全| 国产高清视频在线播放一区| 一个人免费在线观看电影 | 露出奶头的视频| 亚洲av电影不卡..在线观看| 波多野结衣高清无吗| cao死你这个sao货| 亚洲中文日韩欧美视频| 在线免费观看的www视频| 热99re8久久精品国产| 欧美精品啪啪一区二区三区| 床上黄色一级片| 久久精品国产99精品国产亚洲性色| 女同久久另类99精品国产91| 视频区欧美日本亚洲| 精品久久久久久成人av| 此物有八面人人有两片| 99久久无色码亚洲精品果冻| 亚洲午夜理论影院| 精品国产三级普通话版| 免费在线观看成人毛片| 亚洲精品在线美女| 国产伦人伦偷精品视频| 成人欧美大片| 18禁国产床啪视频网站| 很黄的视频免费| 亚洲欧美精品综合一区二区三区| 免费看美女性在线毛片视频| 成人欧美大片| 国产精品,欧美在线| 免费在线观看成人毛片| 亚洲 欧美 日韩 在线 免费| 亚洲aⅴ乱码一区二区在线播放| 国产精品电影一区二区三区| 美女黄网站色视频| 精品久久久久久成人av| 亚洲自拍偷在线| 久久午夜亚洲精品久久| 黄色成人免费大全| 91在线精品国自产拍蜜月 | 大型黄色视频在线免费观看| 日韩高清综合在线| 国产精品国产高清国产av| 变态另类丝袜制服| 国产人伦9x9x在线观看| avwww免费| 无人区码免费观看不卡| 久久精品91蜜桃| 国产精品女同一区二区软件 | 亚洲国产精品999在线| 国产精品久久视频播放| 99久久99久久久精品蜜桃| 精品国产三级普通话版| 18禁黄网站禁片午夜丰满| 欧美精品啪啪一区二区三区| 噜噜噜噜噜久久久久久91| www.精华液| 日韩人妻高清精品专区| 男人的好看免费观看在线视频| 国产亚洲精品综合一区在线观看| 激情在线观看视频在线高清| 99精品欧美一区二区三区四区| 老司机午夜福利在线观看视频| 午夜福利在线观看免费完整高清在 | 午夜影院日韩av| www.www免费av| 午夜日韩欧美国产| 精品国产美女av久久久久小说| 亚洲人成电影免费在线| av福利片在线观看| 香蕉丝袜av| svipshipincom国产片| 成人特级黄色片久久久久久久| 色吧在线观看| 偷拍熟女少妇极品色| 91麻豆av在线| 亚洲精品在线观看二区| 午夜精品在线福利| 一级黄色大片毛片| 一夜夜www| 久久热在线av| 男人舔女人下体高潮全视频| 欧美日韩一级在线毛片| 色老头精品视频在线观看| 国产一区二区在线av高清观看| 一卡2卡三卡四卡精品乱码亚洲| 身体一侧抽搐| 99riav亚洲国产免费| 国产伦在线观看视频一区| 国产精品一区二区三区四区久久| 日韩欧美国产在线观看| 国产野战对白在线观看| 1024手机看黄色片| 又黄又爽又免费观看的视频| 九九热线精品视视频播放| www.熟女人妻精品国产| av国产免费在线观看| 久久久国产成人免费| 婷婷精品国产亚洲av| 黄色视频,在线免费观看| 久久久久久久久免费视频了| 久久久成人免费电影| 夜夜夜夜夜久久久久| 久久精品91蜜桃| 亚洲欧美日韩卡通动漫| 小蜜桃在线观看免费完整版高清| 免费高清视频大片| 九色国产91popny在线| 成人三级黄色视频| 亚洲成a人片在线一区二区| 日韩中文字幕欧美一区二区| 免费人成视频x8x8入口观看| 国产欧美日韩精品一区二区| 日韩欧美三级三区| 亚洲国产精品久久男人天堂| 欧美av亚洲av综合av国产av| 狂野欧美白嫩少妇大欣赏| 欧美日韩综合久久久久久 | 三级毛片av免费| 夜夜看夜夜爽夜夜摸| 欧美成人性av电影在线观看| 熟妇人妻久久中文字幕3abv| 久久久国产精品麻豆| 国产成人av激情在线播放| 国产蜜桃级精品一区二区三区| 精品日产1卡2卡| 好男人电影高清在线观看| 男人的好看免费观看在线视频| 免费无遮挡裸体视频| 少妇裸体淫交视频免费看高清| 丁香六月欧美| 欧美黑人欧美精品刺激| 久久久精品大字幕| 亚洲av成人av| 9191精品国产免费久久| 波多野结衣高清作品| 色尼玛亚洲综合影院| 精华霜和精华液先用哪个| 亚洲人成伊人成综合网2020| 波多野结衣高清无吗| 欧美日韩福利视频一区二区| 99久久久亚洲精品蜜臀av| 又大又爽又粗| av福利片在线观看| 这个男人来自地球电影免费观看| 麻豆成人av在线观看| 久久久久久久精品吃奶| 白带黄色成豆腐渣| 日本 欧美在线| 国产毛片a区久久久久| 极品教师在线免费播放| 久久精品亚洲精品国产色婷小说| 日本免费a在线| 中文字幕高清在线视频| 噜噜噜噜噜久久久久久91| 看片在线看免费视频| 九色成人免费人妻av| 亚洲欧美日韩高清在线视频| 久久久水蜜桃国产精品网| 天堂网av新在线| 在线免费观看的www视频| av黄色大香蕉| 丰满人妻一区二区三区视频av | 国产1区2区3区精品| 国产精品一区二区三区四区久久| 国产精品一区二区三区四区久久| 国产精品一区二区免费欧美| 国产成年人精品一区二区| 亚洲欧洲精品一区二区精品久久久| 精品国产三级普通话版| 人妻久久中文字幕网| av天堂在线播放| 精品一区二区三区四区五区乱码| 99国产极品粉嫩在线观看| 一进一出抽搐动态| 国产精品影院久久| 两人在一起打扑克的视频| 人人妻,人人澡人人爽秒播| 非洲黑人性xxxx精品又粗又长| 亚洲专区中文字幕在线| 美女 人体艺术 gogo| 老熟妇仑乱视频hdxx| 精品一区二区三区视频在线 | 19禁男女啪啪无遮挡网站| 亚洲在线自拍视频| 一二三四社区在线视频社区8| 久久久精品大字幕| 国产精品九九99| 亚洲 欧美一区二区三区| 亚洲av成人av| 小蜜桃在线观看免费完整版高清| 成在线人永久免费视频| 99国产精品一区二区三区| 一夜夜www| 欧美精品啪啪一区二区三区| 嫩草影院入口| 国产亚洲av高清不卡| 久久久久国内视频| 午夜福利在线观看免费完整高清在 | 久久久久久久久久黄片| 人妻丰满熟妇av一区二区三区| 亚洲精品美女久久久久99蜜臀| 久久久久久久久免费视频了| 19禁男女啪啪无遮挡网站| 日日摸夜夜添夜夜添小说| 亚洲中文av在线| 亚洲中文日韩欧美视频| 91av网站免费观看| 成年版毛片免费区| 午夜免费激情av| 精品一区二区三区视频在线观看免费| 成人av在线播放网站| 999久久久国产精品视频| 三级毛片av免费| 亚洲真实伦在线观看| 成人鲁丝片一二三区免费| 搡老岳熟女国产| 欧美极品一区二区三区四区| 黄色女人牲交| 两个人看的免费小视频| 亚洲熟妇熟女久久| 国产成人精品久久二区二区免费| 国产伦人伦偷精品视频| 国产精品98久久久久久宅男小说| 禁无遮挡网站| a在线观看视频网站| 成年女人看的毛片在线观看| 麻豆成人午夜福利视频| 此物有八面人人有两片| 亚洲精品456在线播放app | 精品国产乱子伦一区二区三区| 亚洲国产欧洲综合997久久,| 中文字幕av在线有码专区| 黄色视频,在线免费观看| 欧美绝顶高潮抽搐喷水| 免费搜索国产男女视频| 亚洲成人免费电影在线观看| 日韩有码中文字幕| 欧美日本亚洲视频在线播放| 人妻久久中文字幕网| 国产蜜桃级精品一区二区三区| 欧美大码av| 老汉色av国产亚洲站长工具| 欧美高清成人免费视频www| 免费看美女性在线毛片视频| 精品福利观看| 午夜精品在线福利| 神马国产精品三级电影在线观看| 国产v大片淫在线免费观看| 久久亚洲精品不卡| 九色国产91popny在线| 哪里可以看免费的av片| 搡老熟女国产l中国老女人| 日日干狠狠操夜夜爽| 精品福利观看| 1024手机看黄色片| 2021天堂中文幕一二区在线观| 色综合欧美亚洲国产小说| 久久久成人免费电影| 亚洲人成网站高清观看| 一个人免费在线观看的高清视频| 一二三四在线观看免费中文在| 69av精品久久久久久| 免费av不卡在线播放| 俺也久久电影网| 我的老师免费观看完整版| 99久久国产精品久久久| www日本黄色视频网| 久久久久精品国产欧美久久久| 亚洲av美国av| 亚洲国产欧美网| 午夜亚洲福利在线播放| 在线a可以看的网站| 1024手机看黄色片| 国产三级在线视频| 99热这里只有精品一区 | 久久草成人影院| 国内揄拍国产精品人妻在线| 亚洲欧美日韩高清在线视频| 国产精品精品国产色婷婷| 岛国在线观看网站| 黄片大片在线免费观看| 国产乱人视频| 欧美极品一区二区三区四区| 亚洲18禁久久av| 亚洲av片天天在线观看| 亚洲欧美日韩卡通动漫| 午夜福利在线观看免费完整高清在 | 久久久国产精品麻豆| 免费在线观看视频国产中文字幕亚洲| 丝袜人妻中文字幕| 日韩人妻高清精品专区| 精品日产1卡2卡| 丰满人妻一区二区三区视频av | 色老头精品视频在线观看| 久久久久九九精品影院| 国产三级中文精品| 精品日产1卡2卡| 在线免费观看的www视频| 波多野结衣巨乳人妻| 亚洲av第一区精品v没综合| 欧美日韩综合久久久久久 | 999精品在线视频| 午夜激情欧美在线| 男女下面进入的视频免费午夜| 九色国产91popny在线| a级毛片a级免费在线| 国产毛片a区久久久久| 两性夫妻黄色片| 国产男靠女视频免费网站| 一卡2卡三卡四卡精品乱码亚洲| 久久久久久久久久黄片| 亚洲无线观看免费| 成人无遮挡网站| 成人av一区二区三区在线看| 亚洲第一欧美日韩一区二区三区| 国产精品,欧美在线| 人人妻人人看人人澡| av视频在线观看入口| 精品国产超薄肉色丝袜足j| aaaaa片日本免费| 国产真人三级小视频在线观看| 成人鲁丝片一二三区免费| 午夜福利高清视频| 一级a爱片免费观看的视频| 久久精品影院6| 一边摸一边抽搐一进一小说| 国产成人av激情在线播放| 欧美av亚洲av综合av国产av| 一本久久中文字幕| 波多野结衣巨乳人妻| 12—13女人毛片做爰片一| 性色av乱码一区二区三区2| 国产激情欧美一区二区| 此物有八面人人有两片| a在线观看视频网站| 天天一区二区日本电影三级| 亚洲中文av在线| 成熟少妇高潮喷水视频| 欧美日本亚洲视频在线播放| 成人av在线播放网站| 精品福利观看| 制服丝袜大香蕉在线| 国产亚洲欧美98| 精品不卡国产一区二区三区| 岛国在线免费视频观看| netflix在线观看网站| 久久精品影院6| 亚洲在线观看片| 黄色女人牲交| 精品午夜福利视频在线观看一区| 18禁观看日本| 成人无遮挡网站| 久久久久久久久免费视频了| 久久国产乱子伦精品免费另类| 别揉我奶头~嗯~啊~动态视频| 一个人看的www免费观看视频| 精品国产乱码久久久久久男人| 一本精品99久久精品77| 美女午夜性视频免费| 国产69精品久久久久777片 | 日本成人三级电影网站| 中文字幕熟女人妻在线| 国产成人av激情在线播放| 97超视频在线观看视频| 法律面前人人平等表现在哪些方面| 亚洲中文字幕一区二区三区有码在线看 | 国产乱人视频| 嫁个100分男人电影在线观看| 激情在线观看视频在线高清| 国产精品女同一区二区软件 | 人人妻人人看人人澡| 老汉色av国产亚洲站长工具| 级片在线观看| 午夜a级毛片| 亚洲精品456在线播放app | 亚洲av成人不卡在线观看播放网| 亚洲成av人片免费观看| 白带黄色成豆腐渣| 舔av片在线| 99精品在免费线老司机午夜| 国产私拍福利视频在线观看| 久久草成人影院| 婷婷精品国产亚洲av在线| 久久久色成人| 日本黄色视频三级网站网址| 国产精品久久久久久人妻精品电影| 欧美日韩国产亚洲二区| 十八禁网站免费在线| 久久国产精品人妻蜜桃| 久久中文看片网| 亚洲天堂国产精品一区在线| 国产欧美日韩一区二区三| 国产成人啪精品午夜网站| 国产不卡一卡二| 成人亚洲精品av一区二区| 国产野战对白在线观看| 三级毛片av免费| 老熟妇仑乱视频hdxx| 1024手机看黄色片| 亚洲精品在线美女| 男人的好看免费观看在线视频| 麻豆成人av在线观看| 国产精品久久视频播放| 国产精品99久久99久久久不卡| 69av精品久久久久久| 免费看光身美女| 成人午夜高清在线视频| 国产熟女xx| 久久精品综合一区二区三区| 90打野战视频偷拍视频| 亚洲,欧美精品.| 法律面前人人平等表现在哪些方面| 国产不卡一卡二| 在线播放国产精品三级| 国产精品免费一区二区三区在线| 欧美大码av| 国产精品一区二区三区四区久久| 国产aⅴ精品一区二区三区波| 久久99热这里只有精品18| 成人18禁在线播放| 黄色女人牲交| 观看免费一级毛片| 久久久水蜜桃国产精品网| 99国产精品99久久久久| 俄罗斯特黄特色一大片| 国产伦一二天堂av在线观看| 亚洲av熟女| 少妇熟女aⅴ在线视频| 岛国在线免费视频观看| 欧美日韩精品网址| 久久天躁狠狠躁夜夜2o2o| 日本成人三级电影网站| 成年人黄色毛片网站| 精品电影一区二区在线| 亚洲人与动物交配视频| 日本与韩国留学比较|