• <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
    一级a爱片免费观看的视频| 欧美性感艳星| 俺也久久电影网| 夜夜夜夜夜久久久久| 黄色欧美视频在线观看| 欧美+日韩+精品| 国产亚洲精品久久久久久毛片| 亚洲精品亚洲一区二区| 免费不卡的大黄色大毛片视频在线观看 | 99在线人妻在线中文字幕| 欧美激情在线99| 国产精品美女特级片免费视频播放器| 淫妇啪啪啪对白视频| 色综合色国产| 日韩中字成人| 美女xxoo啪啪120秒动态图| 亚洲av免费高清在线观看| 一卡2卡三卡四卡精品乱码亚洲| 成人精品一区二区免费| 一区二区三区高清视频在线| 丰满的人妻完整版| 校园人妻丝袜中文字幕| 午夜福利在线观看吧| 一个人观看的视频www高清免费观看| 久久韩国三级中文字幕| 欧美在线一区亚洲| 成人漫画全彩无遮挡| 色综合站精品国产| 听说在线观看完整版免费高清| 免费观看人在逋| 久久久久国内视频| 久久九九热精品免费| 噜噜噜噜噜久久久久久91| 午夜福利在线在线| 久久久久久伊人网av| 国产精品久久久久久精品电影| 亚洲国产欧洲综合997久久,| 亚洲七黄色美女视频| 高清毛片免费观看视频网站| 一区二区三区免费毛片| 欧美成人a在线观看| 国产欧美日韩一区二区精品| 女人被狂操c到高潮| 国产单亲对白刺激| 91精品国产九色| 国产美女午夜福利| 中国美女看黄片| 男女视频在线观看网站免费| 成人午夜高清在线视频| 久久天躁狠狠躁夜夜2o2o| 国产亚洲精品久久久com| 国产麻豆成人av免费视频| 欧美3d第一页| 日本三级黄在线观看| 国产精品一区二区免费欧美| 国产伦精品一区二区三区四那| 特级一级黄色大片| 熟女电影av网| 黄片wwwwww| 免费观看人在逋| 桃色一区二区三区在线观看| 久久精品国产亚洲av涩爱 | 亚洲综合色惰| 免费人成视频x8x8入口观看| 九九爱精品视频在线观看| 免费看光身美女| 日本与韩国留学比较| 露出奶头的视频| 乱系列少妇在线播放| 人妻丰满熟妇av一区二区三区| 国产精品久久久久久久电影| 精品乱码久久久久久99久播| 成人无遮挡网站| 午夜精品国产一区二区电影 | 寂寞人妻少妇视频99o| 久久久久久久久大av| 免费搜索国产男女视频| 国产高清三级在线| 欧美+亚洲+日韩+国产| 日韩强制内射视频| 如何舔出高潮| 亚洲成人av在线免费| 九九久久精品国产亚洲av麻豆| 亚洲高清免费不卡视频| 成人综合一区亚洲| 国产精品福利在线免费观看| 欧美成人a在线观看| 天堂影院成人在线观看| 亚洲中文字幕日韩| 亚洲精品国产成人久久av| 99久久久亚洲精品蜜臀av| 亚洲18禁久久av| 在现免费观看毛片| 欧美高清成人免费视频www| 欧美又色又爽又黄视频| 国产精品爽爽va在线观看网站| 麻豆乱淫一区二区| 欧美绝顶高潮抽搐喷水| 亚洲综合色惰| 99热6这里只有精品| 97超碰精品成人国产| 精品一区二区三区人妻视频| 国产精品爽爽va在线观看网站| 美女高潮的动态| 欧美绝顶高潮抽搐喷水| 亚州av有码| 男人和女人高潮做爰伦理| 中文字幕久久专区| 我的老师免费观看完整版| 男女边吃奶边做爰视频| 亚洲欧美精品综合久久99| 波野结衣二区三区在线| 精品午夜福利在线看| 国产人妻一区二区三区在| 成人精品一区二区免费| 自拍偷自拍亚洲精品老妇| 搡老妇女老女人老熟妇| 亚洲不卡免费看| 国产视频一区二区在线看| 国产男人的电影天堂91| 性欧美人与动物交配| 看免费成人av毛片| 久久久久久久久久黄片| 欧美色欧美亚洲另类二区| 国产免费男女视频| 久久中文看片网| 精品一区二区三区视频在线| 成人永久免费在线观看视频| 我的女老师完整版在线观看| 国产麻豆成人av免费视频| 欧美人与善性xxx| 美女免费视频网站| av视频在线观看入口| 久久精品国产99精品国产亚洲性色| 国产欧美日韩一区二区精品| 干丝袜人妻中文字幕| 日韩在线高清观看一区二区三区| 男女做爰动态图高潮gif福利片| 国产精品精品国产色婷婷| 91久久精品国产一区二区三区| 久久久久国内视频| 别揉我奶头 嗯啊视频| 国产伦一二天堂av在线观看| 国产男靠女视频免费网站| 日韩三级伦理在线观看| 秋霞在线观看毛片| 国产在线男女| 麻豆久久精品国产亚洲av| 最新中文字幕久久久久| 综合色丁香网| 午夜激情福利司机影院| 精品日产1卡2卡| 国产精品一区二区三区四区免费观看 | 午夜影院日韩av| 国产亚洲91精品色在线| 国产真实乱freesex| 久久久久精品国产欧美久久久| 又黄又爽又免费观看的视频| 日日啪夜夜撸| 成人特级av手机在线观看| 日本一二三区视频观看| 午夜久久久久精精品| 国产成人freesex在线 | 97超视频在线观看视频| 最后的刺客免费高清国语| 久久草成人影院| 国产午夜精品论理片| 我的老师免费观看完整版| 国产精品国产三级国产av玫瑰| 久久国内精品自在自线图片| 啦啦啦观看免费观看视频高清| 在线观看免费视频日本深夜| 欧美精品国产亚洲| 春色校园在线视频观看| 欧美人与善性xxx| 波多野结衣高清无吗| av视频在线观看入口| 99久国产av精品| 美女被艹到高潮喷水动态| 村上凉子中文字幕在线| 欧美一区二区亚洲| 美女 人体艺术 gogo| 日日摸夜夜添夜夜添小说| 变态另类成人亚洲欧美熟女| 久久人人爽人人爽人人片va| 精品人妻熟女av久视频| 国产精品野战在线观看| 97人妻精品一区二区三区麻豆| 国产精品精品国产色婷婷| 一级黄片播放器| 精品国内亚洲2022精品成人| 波多野结衣高清无吗| 91麻豆精品激情在线观看国产| 美女 人体艺术 gogo| 国产精品电影一区二区三区| 欧美一区二区精品小视频在线| 免费高清视频大片| 久久九九热精品免费| 国产在视频线在精品| 国产乱人视频| 俄罗斯特黄特色一大片| 极品教师在线视频| 网址你懂的国产日韩在线| 亚洲精品色激情综合| 麻豆乱淫一区二区| 不卡视频在线观看欧美| 国产精品日韩av在线免费观看| 亚洲av美国av| 精品久久久久久久久av| 一级黄片播放器| 在线国产一区二区在线| 欧美又色又爽又黄视频| 男人和女人高潮做爰伦理| 国产91av在线免费观看| 麻豆一二三区av精品| 欧美+亚洲+日韩+国产| 免费高清视频大片| 少妇人妻精品综合一区二区 | 久久久精品欧美日韩精品| 女人被狂操c到高潮| 午夜福利在线观看免费完整高清在 | 免费av不卡在线播放| 久久精品夜夜夜夜夜久久蜜豆| 精品一区二区三区视频在线| av专区在线播放| 国内揄拍国产精品人妻在线| 中国美女看黄片| 免费av观看视频| 最好的美女福利视频网| 国产成人影院久久av| 国产色爽女视频免费观看| 白带黄色成豆腐渣| 精品午夜福利在线看| 欧美成人免费av一区二区三区| 色综合色国产| 久久久午夜欧美精品| 久久6这里有精品| 久久99热6这里只有精品| 看片在线看免费视频| 日韩欧美一区二区三区在线观看| 给我免费播放毛片高清在线观看| 免费一级毛片在线播放高清视频| 久久久成人免费电影| 欧美xxxx黑人xx丫x性爽| 国产精品久久电影中文字幕| 日日啪夜夜撸| 一夜夜www| 国产综合懂色| 色综合色国产| 日韩三级伦理在线观看| 少妇猛男粗大的猛烈进出视频 | 中文字幕精品亚洲无线码一区| 看片在线看免费视频| 日韩人妻高清精品专区| 99久久无色码亚洲精品果冻| 国产精品一区二区免费欧美| 婷婷精品国产亚洲av在线| 一区二区三区免费毛片| 一本精品99久久精品77| 国产精品亚洲一级av第二区| 国产精品人妻久久久影院| 一级毛片久久久久久久久女| 免费黄网站久久成人精品| 男女啪啪激烈高潮av片| 99九九线精品视频在线观看视频| 国内精品久久久久精免费| 一区二区三区高清视频在线| 国产精品一区二区三区四区久久| 别揉我奶头 嗯啊视频| av天堂中文字幕网| 一级毛片电影观看 | 黄色日韩在线| 亚洲性夜色夜夜综合| 午夜激情福利司机影院| 国产免费一级a男人的天堂| 久久久成人免费电影| 久久国内精品自在自线图片| h日本视频在线播放| 日本精品一区二区三区蜜桃| 最近在线观看免费完整版| 免费av不卡在线播放| 男女做爰动态图高潮gif福利片| 最新在线观看一区二区三区| 日韩,欧美,国产一区二区三区 | 欧美人与善性xxx| 国产欧美日韩精品亚洲av| 成人性生交大片免费视频hd| 12—13女人毛片做爰片一| 久久99热这里只有精品18| 麻豆av噜噜一区二区三区| 国产伦精品一区二区三区四那| 色5月婷婷丁香| 成年女人永久免费观看视频| 麻豆乱淫一区二区| 12—13女人毛片做爰片一| 欧美性猛交╳xxx乱大交人| 欧美3d第一页| 亚洲av中文av极速乱| 日本免费a在线| 男人狂女人下面高潮的视频| 此物有八面人人有两片| 女人十人毛片免费观看3o分钟| 成年女人毛片免费观看观看9| 少妇熟女aⅴ在线视频| 国产不卡一卡二| 美女高潮的动态| 精品福利观看| 亚洲18禁久久av| 麻豆精品久久久久久蜜桃| 成人永久免费在线观看视频| 国内精品美女久久久久久| 一级毛片久久久久久久久女| 国产色爽女视频免费观看| 嫩草影视91久久| 成人av在线播放网站| 国产亚洲精品av在线| 天堂√8在线中文| 久久久久性生活片| 波多野结衣高清无吗| 别揉我奶头~嗯~啊~动态视频| 午夜视频国产福利| 免费高清视频大片| 桃色一区二区三区在线观看| 精品久久久久久久人妻蜜臀av| 十八禁网站免费在线| 99热6这里只有精品| 欧美一区二区精品小视频在线| 深爱激情五月婷婷| 人妻久久中文字幕网| 少妇人妻一区二区三区视频| 亚洲欧美清纯卡通| 白带黄色成豆腐渣| 精品久久久久久久久av| 亚洲在线自拍视频| 3wmmmm亚洲av在线观看| 亚洲欧美精品综合久久99| 女同久久另类99精品国产91| 俄罗斯特黄特色一大片| 少妇的逼水好多| 久久精品国产自在天天线| 麻豆av噜噜一区二区三区| 成人亚洲欧美一区二区av| 成年女人永久免费观看视频| 精品久久久久久成人av| 国产一区二区在线观看日韩| 可以在线观看的亚洲视频| 免费看光身美女| 高清日韩中文字幕在线| 亚洲图色成人| 欧美又色又爽又黄视频| 午夜免费激情av| 人人妻人人澡人人爽人人夜夜 | 欧美日本视频| 国产蜜桃级精品一区二区三区| 国产精品久久久久久av不卡| 亚洲成人av在线免费| 亚洲内射少妇av| 欧美日本视频| 国产成人91sexporn| 日本撒尿小便嘘嘘汇集6| 日日摸夜夜添夜夜添小说| 亚洲熟妇熟女久久| 人妻久久中文字幕网| 大香蕉久久网| 欧美xxxx黑人xx丫x性爽| 国产激情偷乱视频一区二区| 午夜a级毛片| 免费观看的影片在线观看| 熟妇人妻久久中文字幕3abv| 久久99热这里只有精品18| 国产黄色视频一区二区在线观看 | 插阴视频在线观看视频| 在线国产一区二区在线| 激情 狠狠 欧美| 国产单亲对白刺激| 悠悠久久av| 男女之事视频高清在线观看| 99久久九九国产精品国产免费| 哪里可以看免费的av片| 在线播放国产精品三级| 嫩草影视91久久| 97热精品久久久久久| 18+在线观看网站| 久久99热这里只有精品18| 免费在线观看成人毛片| 寂寞人妻少妇视频99o| 校园人妻丝袜中文字幕| 联通29元200g的流量卡| 欧美一区二区国产精品久久精品| 免费黄网站久久成人精品| 亚洲精品亚洲一区二区| 午夜免费激情av| 又爽又黄a免费视频| 99久久中文字幕三级久久日本| 麻豆精品久久久久久蜜桃| 成人亚洲欧美一区二区av| 九九热线精品视视频播放| 成人精品一区二区免费| 日韩三级伦理在线观看| 国产91av在线免费观看| 男人舔奶头视频| 精品无人区乱码1区二区| 亚洲精品一卡2卡三卡4卡5卡| 九九热线精品视视频播放| 国产久久久一区二区三区| 亚洲一区二区三区色噜噜| 日本免费a在线| 久久精品国产鲁丝片午夜精品| 久久久久久九九精品二区国产| 成人午夜高清在线视频| 久久久久国产网址| 成人特级av手机在线观看| 亚洲第一电影网av| 国内精品美女久久久久久| 淫秽高清视频在线观看| 亚洲在线自拍视频| www日本黄色视频网| 亚洲av电影不卡..在线观看| 男人舔奶头视频| 一本久久中文字幕| 精品久久久久久久久久免费视频| 免费高清视频大片| 亚洲图色成人| 人人妻人人澡欧美一区二区| 国产精品嫩草影院av在线观看| av在线蜜桃| 久久精品91蜜桃| 亚洲美女视频黄频| 免费高清视频大片| 国产一区二区在线观看日韩| 久久精品夜夜夜夜夜久久蜜豆| 亚洲七黄色美女视频| 亚洲一区高清亚洲精品| 男女那种视频在线观看| 亚洲人成网站在线观看播放| 精品不卡国产一区二区三区| 99国产极品粉嫩在线观看| 亚洲精品日韩在线中文字幕 | 淫秽高清视频在线观看| 精品99又大又爽又粗少妇毛片| 婷婷亚洲欧美| 免费电影在线观看免费观看| 美女内射精品一级片tv| 级片在线观看| 身体一侧抽搐| 搡老岳熟女国产| 亚洲国产欧美人成| 午夜爱爱视频在线播放| 国产精品一区二区免费欧美| 亚洲av二区三区四区| 亚洲成人中文字幕在线播放| 国产爱豆传媒在线观看| 一区福利在线观看| 我要看日韩黄色一级片| 欧美3d第一页| 国语自产精品视频在线第100页| avwww免费| aaaaa片日本免费| 中文资源天堂在线| www日本黄色视频网| 精品久久久久久久末码| 免费看光身美女| 高清日韩中文字幕在线| 欧美xxxx性猛交bbbb| 国产精品不卡视频一区二区| 床上黄色一级片| 桃色一区二区三区在线观看| 伦理电影大哥的女人| 亚洲欧美日韩无卡精品| 成人午夜高清在线视频| 级片在线观看| 91久久精品电影网| 国产 一区精品| 欧美日韩乱码在线| 亚洲七黄色美女视频| 如何舔出高潮| 不卡视频在线观看欧美| 夜夜看夜夜爽夜夜摸| 十八禁网站免费在线| 能在线免费观看的黄片| 久久亚洲国产成人精品v| 自拍偷自拍亚洲精品老妇| 久久6这里有精品| 久久久久久久亚洲中文字幕| 国产精品久久久久久亚洲av鲁大| 国产麻豆成人av免费视频| 精品熟女少妇av免费看| 国产成人a∨麻豆精品| 人妻夜夜爽99麻豆av| 性插视频无遮挡在线免费观看| 久久亚洲国产成人精品v| a级毛片免费高清观看在线播放| 国产精品99久久久久久久久| 偷拍熟女少妇极品色| 波多野结衣巨乳人妻| 99久久中文字幕三级久久日本| 极品教师在线视频| 男人舔奶头视频| 亚洲国产精品成人久久小说 | 永久网站在线| 性色avwww在线观看| 久久鲁丝午夜福利片| 日韩高清综合在线| 国产精品1区2区在线观看.| 丝袜美腿在线中文| 一进一出抽搐gif免费好疼| 国产精华一区二区三区| 亚洲电影在线观看av| 国产精品一二三区在线看| 大型黄色视频在线免费观看| 麻豆成人午夜福利视频| 欧美日韩精品成人综合77777| 麻豆乱淫一区二区| 此物有八面人人有两片| 好男人在线观看高清免费视频| 精品一区二区三区视频在线观看免费| 亚洲婷婷狠狠爱综合网| 国产探花在线观看一区二区| а√天堂www在线а√下载| 国产成人一区二区在线| 中文字幕熟女人妻在线| 国产又黄又爽又无遮挡在线| 亚洲性夜色夜夜综合| 亚洲中文日韩欧美视频| 国产av不卡久久| 国产大屁股一区二区在线视频| eeuss影院久久| 久久草成人影院| 亚洲精品在线观看二区| 亚洲七黄色美女视频| 中文字幕人妻熟人妻熟丝袜美| 免费看光身美女| 97人妻精品一区二区三区麻豆| 给我免费播放毛片高清在线观看| 99久久精品一区二区三区| av中文乱码字幕在线| 性插视频无遮挡在线免费观看| 美女被艹到高潮喷水动态| 久久精品国产自在天天线| 熟女人妻精品中文字幕| 久久久久久久久久成人| 久久久久久久久大av| 熟女人妻精品中文字幕| 波野结衣二区三区在线| 九九热线精品视视频播放| 亚洲国产精品久久男人天堂| 日韩欧美 国产精品| 天美传媒精品一区二区| 国产成人a区在线观看| 午夜影院日韩av| 国产亚洲精品久久久com| 亚洲欧美成人精品一区二区| 少妇被粗大猛烈的视频| 久久精品国产亚洲网站| 免费av不卡在线播放| 国产成人影院久久av| 日韩成人伦理影院| 男女之事视频高清在线观看| 亚洲图色成人| 99久久中文字幕三级久久日本| 女人被狂操c到高潮| 国产亚洲精品久久久久久毛片| 深爱激情五月婷婷| 少妇熟女aⅴ在线视频| 亚洲国产欧美人成| 男女那种视频在线观看| 99在线人妻在线中文字幕| 久久草成人影院| 日日啪夜夜撸| 男人狂女人下面高潮的视频| 91午夜精品亚洲一区二区三区| 六月丁香七月| 亚洲无线观看免费| 久久精品影院6| 色尼玛亚洲综合影院| 欧美激情国产日韩精品一区| 99在线人妻在线中文字幕| 亚洲精品成人久久久久久| 日韩精品中文字幕看吧| 久久久久久久久大av| 国产麻豆成人av免费视频| 亚洲成人av在线免费| 在线观看美女被高潮喷水网站| 免费看光身美女| 亚洲无线在线观看| 亚洲精品乱码久久久v下载方式| 又黄又爽又免费观看的视频| 成年女人永久免费观看视频| 欧美极品一区二区三区四区| 精品人妻偷拍中文字幕| 91在线精品国自产拍蜜月| 国产欧美日韩精品一区二区| or卡值多少钱| 国产精品一区二区三区四区免费观看 | 成人欧美大片| 97人妻精品一区二区三区麻豆| 乱系列少妇在线播放| 少妇熟女aⅴ在线视频| 国产精品一区二区三区四区免费观看 | 悠悠久久av| 久久午夜福利片| 日本撒尿小便嘘嘘汇集6| 国产真实乱freesex| 精品午夜福利视频在线观看一区| 久久久久国产精品人妻aⅴ院| 日本 av在线| 99riav亚洲国产免费| 久久久久久久亚洲中文字幕| 亚洲中文字幕日韩| 有码 亚洲区| 高清毛片免费看| 精品人妻视频免费看| 成年免费大片在线观看| 美女大奶头视频| 欧美高清性xxxxhd video|