• <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
    √禁漫天堂资源中文www| 日本猛色少妇xxxxx猛交久久| 丰满迷人的少妇在线观看| 国产亚洲一区二区精品| 日本黄色日本黄色录像| 人人妻人人澡人人看| 精品国产国语对白av| 精品久久久噜噜| 亚洲激情五月婷婷啪啪| 最近中文字幕高清免费大全6| 丰满饥渴人妻一区二区三| 国产成人精品福利久久| 国产精品国产三级国产av玫瑰| 丝袜脚勾引网站| 免费看不卡的av| 国产国语露脸激情在线看| 亚洲四区av| 久久毛片免费看一区二区三区| 国产精品一国产av| 精品久久久精品久久久| 99视频精品全部免费 在线| videos熟女内射| freevideosex欧美| 18禁在线无遮挡免费观看视频| 久久毛片免费看一区二区三区| 午夜福利网站1000一区二区三区| 免费观看无遮挡的男女| 午夜视频国产福利| 久久久久视频综合| 午夜福利视频在线观看免费| 五月玫瑰六月丁香| h视频一区二区三区| 亚洲精品国产色婷婷电影| 亚洲精华国产精华液的使用体验| 丰满迷人的少妇在线观看| 国产综合精华液| 亚洲精品乱久久久久久| 欧美精品人与动牲交sv欧美| 日韩成人伦理影院| 18禁观看日本| av电影中文网址| 欧美日韩在线观看h| 久久精品国产鲁丝片午夜精品| 亚洲精品乱久久久久久| 成人午夜精彩视频在线观看| 黑人猛操日本美女一级片| 人人妻人人爽人人添夜夜欢视频| 国产av国产精品国产| 久久久久久久久久人人人人人人| 美女cb高潮喷水在线观看| 国产女主播在线喷水免费视频网站| 51国产日韩欧美| 成人手机av| 精品99又大又爽又粗少妇毛片| 草草在线视频免费看| 精品亚洲成a人片在线观看| 亚洲精品色激情综合| 日韩精品免费视频一区二区三区 | 蜜桃久久精品国产亚洲av| 色婷婷av一区二区三区视频| 国产高清三级在线| 午夜久久久在线观看| 久久av网站| 国产成人免费无遮挡视频| 黄色一级大片看看| 亚洲精品国产av成人精品| 中文天堂在线官网| 日韩一区二区三区影片| 国产精品无大码| 极品少妇高潮喷水抽搐| 又粗又硬又长又爽又黄的视频| 亚洲精华国产精华液的使用体验| 中文字幕av电影在线播放| 狠狠婷婷综合久久久久久88av| 97精品久久久久久久久久精品| 精品久久久久久久久亚洲| av黄色大香蕉| 91久久精品国产一区二区三区| 看免费成人av毛片| 最近手机中文字幕大全| 亚洲av国产av综合av卡| 狠狠精品人妻久久久久久综合| 久久久久视频综合| 国产伦理片在线播放av一区| 亚洲成人手机| 人成视频在线观看免费观看| 国产69精品久久久久777片| 少妇人妻 视频| 久热这里只有精品99| 一区二区三区四区激情视频| 男女免费视频国产| 妹子高潮喷水视频| 日韩大片免费观看网站| 中国美白少妇内射xxxbb| 精品国产国语对白av| 欧美97在线视频| 最近中文字幕高清免费大全6| 欧美xxⅹ黑人| 亚洲成人一二三区av| 免费不卡的大黄色大毛片视频在线观看| 亚洲人成网站在线观看播放| 91精品伊人久久大香线蕉| 日韩人妻高清精品专区| 中文字幕免费在线视频6| 高清av免费在线| 亚洲中文av在线| 欧美精品亚洲一区二区| 又大又黄又爽视频免费| 考比视频在线观看| av国产久精品久网站免费入址| 日韩av在线免费看完整版不卡| 老熟女久久久| 黄片播放在线免费| 欧美xxⅹ黑人| 亚洲美女搞黄在线观看| 亚洲欧洲精品一区二区精品久久久 | 最黄视频免费看| 亚洲精品视频女| 日本黄大片高清| 看十八女毛片水多多多| 婷婷成人精品国产| 亚洲精品国产色婷婷电影| 久久人妻熟女aⅴ| 午夜精品国产一区二区电影| 狂野欧美激情性xxxx在线观看| 天堂中文最新版在线下载| 中文字幕免费在线视频6| 青春草亚洲视频在线观看| 大话2 男鬼变身卡| 国产在线一区二区三区精| 亚洲国产成人一精品久久久| 亚洲精品中文字幕在线视频| 欧美xxxx性猛交bbbb| 欧美 亚洲 国产 日韩一| 久久 成人 亚洲| 观看美女的网站| 亚洲综合色惰| 成人免费观看视频高清| 如何舔出高潮| 在线观看美女被高潮喷水网站| 欧美精品人与动牲交sv欧美| 新久久久久国产一级毛片| 综合色丁香网| 婷婷色麻豆天堂久久| 精品一品国产午夜福利视频| 一本大道久久a久久精品| 国产极品粉嫩免费观看在线 | 亚洲中文av在线| 亚洲精品456在线播放app| 男女边吃奶边做爰视频| 免费观看在线日韩| 亚洲精品国产av蜜桃| 中国国产av一级| 日韩人妻高清精品专区| 亚洲第一av免费看| 插阴视频在线观看视频| 国产精品秋霞免费鲁丝片| 国产在线免费精品| 飞空精品影院首页| 最近中文字幕高清免费大全6| 街头女战士在线观看网站| kizo精华| 麻豆精品久久久久久蜜桃| 欧美国产精品一级二级三级| 国产免费一级a男人的天堂| 欧美精品高潮呻吟av久久| 精品久久久精品久久久| 亚洲av.av天堂| 久久久久久久久久人人人人人人| 亚洲美女搞黄在线观看| 国产成人a∨麻豆精品| 久久午夜福利片| 成人国产麻豆网| videosex国产| 久久久久久久久久久丰满| 国模一区二区三区四区视频| 亚洲综合色惰| 大话2 男鬼变身卡| 一区二区日韩欧美中文字幕 | 亚洲精品日韩av片在线观看| 99国产精品免费福利视频| 欧美日韩亚洲高清精品| 毛片一级片免费看久久久久| 精品久久久久久电影网| 肉色欧美久久久久久久蜜桃| 免费人成在线观看视频色| 欧美xxxx性猛交bbbb| 黄片播放在线免费| 又大又黄又爽视频免费| 人体艺术视频欧美日本| 日本欧美视频一区| 日韩 亚洲 欧美在线| 成人18禁高潮啪啪吃奶动态图 | 日韩精品有码人妻一区| 日韩欧美一区视频在线观看| 最近手机中文字幕大全| 亚洲精品日韩在线中文字幕| 午夜精品国产一区二区电影| 91久久精品国产一区二区三区| 色94色欧美一区二区| 亚洲美女视频黄频| 久热这里只有精品99| 18禁在线无遮挡免费观看视频| 亚洲av中文av极速乱| 国产乱人偷精品视频| 国产亚洲午夜精品一区二区久久| 伦理电影大哥的女人| 最近的中文字幕免费完整| 草草在线视频免费看| a级片在线免费高清观看视频| 久久综合国产亚洲精品| 精品国产露脸久久av麻豆| 一边摸一边做爽爽视频免费| 在线观看三级黄色| 蜜桃在线观看..| 国产精品蜜桃在线观看| 亚洲国产精品国产精品| 精品一区二区三卡| 成年女人在线观看亚洲视频| 日日摸夜夜添夜夜爱| 国产精品国产三级专区第一集| 免费观看在线日韩| 国产成人一区二区在线| 丝袜喷水一区| 亚洲精品456在线播放app| 欧美 亚洲 国产 日韩一| 国产精品无大码| 成年人免费黄色播放视频| 在线亚洲精品国产二区图片欧美 | 黄色视频在线播放观看不卡| 十分钟在线观看高清视频www| 国产av精品麻豆| 又黄又爽又刺激的免费视频.| 精品少妇内射三级| 亚洲av男天堂| 成人亚洲欧美一区二区av| 伊人亚洲综合成人网| 爱豆传媒免费全集在线观看| 午夜免费鲁丝| 麻豆精品久久久久久蜜桃| 十八禁高潮呻吟视频| 日本黄色日本黄色录像| 2022亚洲国产成人精品| 久久人人爽人人爽人人片va| av不卡在线播放| 丝袜在线中文字幕| 免费看不卡的av| 日韩av不卡免费在线播放| 亚洲天堂av无毛| 男女高潮啪啪啪动态图| 99九九在线精品视频| 日韩一区二区三区影片| 91成人精品电影| 亚洲精品久久成人aⅴ小说 | 亚洲欧美精品自产自拍| 成人黄色视频免费在线看| 欧美少妇被猛烈插入视频| 亚洲四区av| 亚洲人成网站在线播| 久久久精品免费免费高清| 免费观看在线日韩| 精品一区二区免费观看| 校园人妻丝袜中文字幕| 最近的中文字幕免费完整| 免费少妇av软件| 精品酒店卫生间| 一个人看视频在线观看www免费| 美女主播在线视频| 青青草视频在线视频观看| 狂野欧美激情性bbbbbb| 国产日韩欧美在线精品| 全区人妻精品视频| 国产国语露脸激情在线看| 欧美性感艳星| 日韩 亚洲 欧美在线| 亚洲精品一区蜜桃| 国产色婷婷99| 91精品三级在线观看| 两个人的视频大全免费| 永久网站在线| 桃花免费在线播放| 丝袜美足系列| 国产精品一二三区在线看| 亚洲丝袜综合中文字幕| 欧美亚洲 丝袜 人妻 在线| 久久人人爽人人爽人人片va| 婷婷色综合www| 特大巨黑吊av在线直播| 午夜免费观看性视频| 午夜福利网站1000一区二区三区| 亚洲国产精品专区欧美| 午夜免费鲁丝| 亚洲国产色片| 国产欧美日韩一区二区三区在线 | 新久久久久国产一级毛片| 国产日韩欧美在线精品| 国产白丝娇喘喷水9色精品| 亚洲精品视频女| 搡女人真爽免费视频火全软件| 18禁在线播放成人免费| 免费黄色在线免费观看| 免费少妇av软件| 国产成人精品在线电影| 中文欧美无线码| 在线观看三级黄色| 少妇猛男粗大的猛烈进出视频| 天堂俺去俺来也www色官网| 亚洲精品久久成人aⅴ小说 | 久久热精品热| 亚洲精品日本国产第一区| 91精品三级在线观看| 免费观看a级毛片全部| 亚洲精品成人av观看孕妇| 日韩成人伦理影院| 国产精品一区二区在线不卡| 亚洲精品一区蜜桃| 成人国产麻豆网| 午夜激情av网站| av.在线天堂| 日本av手机在线免费观看| 自线自在国产av| 在线观看免费日韩欧美大片 | 狠狠婷婷综合久久久久久88av| 亚洲成人手机| 欧美最新免费一区二区三区| 黄色毛片三级朝国网站| 在线播放无遮挡| 一级a做视频免费观看| 18禁动态无遮挡网站| 亚洲成色77777| 国产高清有码在线观看视频| 人妻 亚洲 视频| av又黄又爽大尺度在线免费看| 特大巨黑吊av在线直播| 国产黄色免费在线视频| 少妇 在线观看| 成人免费观看视频高清| 国产成人免费无遮挡视频| 80岁老熟妇乱子伦牲交| 18禁裸乳无遮挡动漫免费视频| 日日啪夜夜爽| 国产精品一国产av| 国产无遮挡羞羞视频在线观看| 插逼视频在线观看| 搡老乐熟女国产| 久久精品国产a三级三级三级| 99热国产这里只有精品6| 熟女人妻精品中文字幕| 国产精品一区二区在线观看99| 狂野欧美激情性xxxx在线观看| 婷婷成人精品国产| 国产精品无大码| 一本一本综合久久| 国产精品.久久久| 有码 亚洲区| 欧美成人午夜免费资源| 国产一区二区在线观看av| 51国产日韩欧美| 久久久久久久久久人人人人人人| 日本av手机在线免费观看| 99re6热这里在线精品视频| 大片电影免费在线观看免费| 亚洲精品乱码久久久v下载方式| 中文天堂在线官网| 最近2019中文字幕mv第一页| 成年美女黄网站色视频大全免费 | 国产深夜福利视频在线观看| 日韩精品有码人妻一区| 韩国av在线不卡| 欧美日韩综合久久久久久| 在线观看国产h片| 久久久精品免费免费高清| 亚洲综合色惰| 亚洲av综合色区一区| 欧美日韩一区二区视频在线观看视频在线| √禁漫天堂资源中文www| h视频一区二区三区| 亚洲国产av新网站| 少妇熟女欧美另类| 免费观看av网站的网址| 精品一区在线观看国产| 日韩在线高清观看一区二区三区| 久久久久视频综合| 亚洲av不卡在线观看| 99国产综合亚洲精品| 在线观看美女被高潮喷水网站| 国产色婷婷99| 日日撸夜夜添| 国产深夜福利视频在线观看| 汤姆久久久久久久影院中文字幕| 狂野欧美白嫩少妇大欣赏| 亚洲av二区三区四区| 欧美日韩国产mv在线观看视频| 国产精品国产av在线观看| 少妇人妻久久综合中文| 王馨瑶露胸无遮挡在线观看| 五月开心婷婷网| 国产精品三级大全| 久久久精品免费免费高清| 极品人妻少妇av视频| 久久精品久久久久久久性| 日本黄大片高清| 丝袜脚勾引网站| 人妻一区二区av| 日韩三级伦理在线观看| 极品少妇高潮喷水抽搐| 亚洲精品乱码久久久久久按摩| 国产 精品1| 精品国产一区二区三区久久久樱花| 天美传媒精品一区二区| 日韩欧美一区视频在线观看| 久久午夜综合久久蜜桃| 日韩精品免费视频一区二区三区 | 国产成人精品无人区| 久久97久久精品| 亚洲欧美成人精品一区二区| 久久久精品94久久精品| 亚洲国产毛片av蜜桃av| 色吧在线观看| 日韩电影二区| 久久久久久久久大av| 国产69精品久久久久777片| a级毛片在线看网站| 水蜜桃什么品种好| 午夜福利在线观看免费完整高清在| 一二三四中文在线观看免费高清| 精品人妻熟女毛片av久久网站| 18禁裸乳无遮挡动漫免费视频| 日韩电影二区| 精品午夜福利在线看| 欧美成人午夜免费资源| 中文欧美无线码| 18禁动态无遮挡网站| 一区二区三区精品91| 99热全是精品| 91午夜精品亚洲一区二区三区| 久久精品久久久久久久性| 国产免费一级a男人的天堂| 日韩欧美精品免费久久| 女人久久www免费人成看片| videossex国产| 国产在线免费精品| 精品少妇久久久久久888优播| 久久精品人人爽人人爽视色| 日日摸夜夜添夜夜爱| 一个人看视频在线观看www免费| 日本vs欧美在线观看视频| 日本欧美视频一区| 各种免费的搞黄视频| 内地一区二区视频在线| 日本色播在线视频| av在线app专区| 成人国产麻豆网| 日韩一区二区三区影片| 亚洲欧洲国产日韩| 成年av动漫网址| 我的女老师完整版在线观看| 香蕉精品网在线| 国产成人91sexporn| 一个人免费看片子| av免费在线看不卡| 天天影视国产精品| 日韩不卡一区二区三区视频在线| 精品一区二区三卡| 亚洲精品色激情综合| 精品人妻偷拍中文字幕| 国模一区二区三区四区视频| 日本黄色日本黄色录像| 色视频在线一区二区三区| 久热这里只有精品99| 啦啦啦在线观看免费高清www| 狠狠精品人妻久久久久久综合| 国产国语露脸激情在线看| 内地一区二区视频在线| 精品国产一区二区久久| 国产精品麻豆人妻色哟哟久久| 人人妻人人澡人人看| 久久久久久久大尺度免费视频| 成人国语在线视频| 国产亚洲av片在线观看秒播厂| 一级毛片电影观看| 欧美变态另类bdsm刘玥| 卡戴珊不雅视频在线播放| 亚洲国产精品999| 久久国产亚洲av麻豆专区| 一个人看视频在线观看www免费| 日韩欧美精品免费久久| 夜夜看夜夜爽夜夜摸| 亚洲人与动物交配视频| 天天躁夜夜躁狠狠久久av| a 毛片基地| 久久免费观看电影| 日本欧美国产在线视频| 如日韩欧美国产精品一区二区三区 | 热99国产精品久久久久久7| 亚洲一区二区三区欧美精品| 亚洲人与动物交配视频| 欧美bdsm另类| 亚洲不卡免费看| 桃花免费在线播放| 国产免费又黄又爽又色| 少妇人妻 视频| 精品久久久噜噜| 国产精品无大码| 国产亚洲av片在线观看秒播厂| 久久久久精品久久久久真实原创| 青春草视频在线免费观看| 国产精品免费大片| 大香蕉久久网| 男男h啪啪无遮挡| 欧美+日韩+精品| 最后的刺客免费高清国语| 成人亚洲精品一区在线观看| 亚洲欧美中文字幕日韩二区| 视频在线观看一区二区三区| 国内精品宾馆在线| 国产精品偷伦视频观看了| 自线自在国产av| 久久精品久久久久久噜噜老黄| 久久狼人影院| 黄色怎么调成土黄色| 人人妻人人爽人人添夜夜欢视频| 一级毛片电影观看| 欧美激情国产日韩精品一区| 有码 亚洲区| av有码第一页| 国产高清有码在线观看视频| 国产成人精品无人区| 熟女人妻精品中文字幕| 少妇人妻 视频| 亚洲av在线观看美女高潮| 国产国语露脸激情在线看| 秋霞伦理黄片| 青春草亚洲视频在线观看| 美女脱内裤让男人舔精品视频| 国产亚洲精品久久久com| 日本免费在线观看一区| 热99国产精品久久久久久7| 中文天堂在线官网| 内地一区二区视频在线| 亚洲怡红院男人天堂| 久热这里只有精品99| 国产精品三级大全| 亚洲精品美女久久av网站| 精品一区二区三区视频在线| 男的添女的下面高潮视频| 热99国产精品久久久久久7| 久久久亚洲精品成人影院| 麻豆精品久久久久久蜜桃| 国产一区二区在线观看日韩| 中文字幕久久专区| 一级毛片电影观看| 人人妻人人澡人人爽人人夜夜| 大香蕉97超碰在线| 黄片无遮挡物在线观看| 超碰97精品在线观看| 欧美日韩视频精品一区| 久久精品国产鲁丝片午夜精品| 国产精品一区二区在线不卡| 伊人久久精品亚洲午夜| 亚洲欧美中文字幕日韩二区| 欧美亚洲 丝袜 人妻 在线| 午夜福利视频在线观看免费| 在线观看三级黄色| 国产成人精品福利久久| 午夜老司机福利剧场| 大香蕉久久成人网| 亚洲人成网站在线播| 国产片内射在线| 一个人看视频在线观看www免费| 精品一区二区三区视频在线| 国产精品国产三级专区第一集| 老女人水多毛片| 视频区图区小说| 久久毛片免费看一区二区三区| 国产高清国产精品国产三级| 美女脱内裤让男人舔精品视频| 久久国内精品自在自线图片| 日本-黄色视频高清免费观看| 丝袜在线中文字幕| 国产欧美日韩综合在线一区二区| 亚洲人成网站在线播| 成人毛片a级毛片在线播放| 亚洲人与动物交配视频| 午夜av观看不卡| 久久久久久人妻| 男女国产视频网站| 欧美激情国产日韩精品一区| av在线播放精品| 日韩av在线免费看完整版不卡| 亚洲精品中文字幕在线视频| 色哟哟·www| 亚洲久久久国产精品| 久久99热这里只频精品6学生| 久久久久久久久久久免费av| 嫩草影院入口| 街头女战士在线观看网站| 中文字幕制服av| 女人精品久久久久毛片| 国产成人freesex在线| 精品国产国语对白av| 纵有疾风起免费观看全集完整版| 国产在线一区二区三区精| 亚洲av福利一区| 国产精品一区二区在线观看99| 一本一本综合久久| 美女大奶头黄色视频| av线在线观看网站| 少妇的逼好多水| 大陆偷拍与自拍| 简卡轻食公司| 在线看a的网站| 丝瓜视频免费看黄片| 日韩三级伦理在线观看| 久久久国产精品麻豆|