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

    HL-2A 高約束先進運行模式等離子體電流剖面集成模擬*

    2021-12-16 07:59:06張洪銘吳靜李佳鮮姚列明徐江城吳巖占劉旗艷郭鵬程
    物理學報 2021年23期
    關(guān)鍵詞:位形托卡馬克等離子體

    張洪銘 吳靜 李佳鮮 姚列明 徐江城吳巖占 劉旗艷 郭鵬程

    1) (電子科技大學物理學院,成都 611731)

    2) (核工業(yè)西南物理研究院聚變研究所,成都 610225)

    磁約束托卡馬克裝置等離子體物理參數(shù),如等離子體電流、約束磁場、安全因子的實時監(jiān)測對高約束模式(H 模)的穩(wěn)態(tài)運行至關(guān)重要.本文介紹了基于HL-2A 實驗數(shù)據(jù)和集成建模手段,重建先進高約束運行模式(H 模)下等離子體位形,電流密度和安全因子參數(shù)剖面.通過搭建動力學平衡位形集成模擬平臺,結(jié)合工作流快速處理手段,采用真實實驗數(shù)據(jù)和集成模擬模型,得到了HL-2A 高約束模式放電條件的離子和電子溫度、密度,電流密度和安全因子q 剖面,重建了高約束模式放電物理圖像和托卡馬克內(nèi)部磁面位形和等離子體邊界參數(shù)連續(xù)分布,計算了等離子體歐姆電流、自舉電流和射頻電流等剖面及各自份額,給出了先進運行模式下安全因子q 徑向分布.通過集成模擬,發(fā)現(xiàn)HL-2A 該炮號參數(shù)條件下臺基寬度約7.52 cm,壓強梯度的方向在ρ(r/a)=0.1 改變,在ρ=0.7 附近數(shù)值達到最大,可能是負剪切產(chǎn)生的內(nèi)部輸運壘位形(ITB)引起的.通過等離子體磁場和電流等剖面重建和實時監(jiān)測,可以評估HL-2A 高約束放電實驗的品質(zhì),為HL-2M 高比壓(βn)先進模式的穩(wěn)態(tài)運行提供參考.

    1 引言

    在磁約束聚變裝置中,高約束模式[1](H 模)相比于低約束模式(L 模),能大幅度提高等離子體性能[2].在高約束模式下,等離子體達到同等水平的運行參數(shù)所需的裝置規(guī)模要小,可大大提高未來聚變能的經(jīng)濟性.磁約束托卡馬克出現(xiàn)H 模首先在德國的磁約束托卡馬克裝置ASDEX 發(fā)現(xiàn),一直是各個聚變裝置的重要研究內(nèi)容,后來通過不斷開展放電物理實驗及理論模擬研究,逐漸掌握低約束模式到高約束模式(L-H 模)的轉(zhuǎn)換機理及對等離子體性能的影響.H 模放電的典型特征之一是壓強剖面(溫度、密度相關(guān))在等離子體邊界附近具有一定臺基,即邊界附近壓強剖面突然變陡峭,壓強梯度迅速增大.臺基的高度和寬度與等離子體性能緊密相關(guān),臺基高度越高,也更容易得到芯部參數(shù)更高的溫度密度剖面.臺基高度和寬度受到動力氣球模和剝離氣球模不穩(wěn)定性的限制[3],可由EPED程序進行預測[4,5].邊界附近壓強梯度的變化較大產(chǎn)生較大自舉電流,影響等離子體電流的剖面分布,進而影響安全因子q的分布.電流剖面和安全因子q徑向分布是研究H 模以及比H 模先進運行模式(參數(shù)更高)的重要分析指標.但是電流分布,尤其是自舉電流分布及安全因子q分布、芯部數(shù)據(jù),難以從實驗診斷實時獲取,一般需要通過輸運模擬程序模擬分析得到.研究人員需要選取物理模型(粒子和流體模型),進行集成建模,將診斷獲取的參數(shù)進行整合,重建電流分布剖面和安全因子剖面,對托卡馬克穩(wěn)態(tài)運行進行實時監(jiān)控.

    磁約束聚變研究涉及很多前沿科學和工程技術(shù)難題,在磁約束托卡馬克裝置中,放電實驗主要依靠診斷設備獲取相關(guān)信號(等離子體溫度和密度,電子溫度和密度,旋轉(zhuǎn)速度等),并轉(zhuǎn)換為便于理解和分析的實驗數(shù)據(jù).由于診斷技術(shù)、診斷空間、經(jīng)費投入等因素,使得實驗獲取的診斷數(shù)據(jù)難以滿足放電實驗實時控制需求.比如,溫度密度剖面數(shù)據(jù),長期難以得到完整的剖面數(shù)據(jù),由于分辨率低,獲取的剖面無法反描述H 模放電的臺基參數(shù)等.因此,需要利用有限的實驗數(shù)據(jù)和理論模型分析,建立綜合的集成模擬平臺,獲取更豐富的放電信息,更加全面評估放電品質(zhì).

    美國通用原子公司(General Atom,GA)開發(fā)了基于DIII-D 實驗數(shù)據(jù)的磁面重建程序EFIT[6,7],可以獲得磁場位形,利用外部的電磁診斷設備(磁通環(huán)和磁探針)獲取磁信號,對實驗數(shù)據(jù)進行全面深入分析,準確重建等離子體平衡位形和相關(guān)的物理參數(shù).由于缺乏芯部診斷數(shù)據(jù),EFIT 平衡反演程序難以給出準確的電流及壓強剖面分布.可以通過增加內(nèi)部診斷手段,利用可見光譜的動態(tài)斯塔克效應(motional Stark effect,MSE)得到等離子體電流剖面分布[8]、等離子體壓強(可見光譜診斷等手段得到的等離子體離子密度和溫度)的診斷數(shù)據(jù),利用EFIT 完成了更為準確的磁面重建和剖面參數(shù)分析.EFIT 應用于磁約束托卡馬克裝置,成為分析實驗數(shù)據(jù)的主要工具,可以重建等離子體邊界和磁場位形參數(shù),也能給出等離子體壓強、電流密度及安全因子分布剖面.隨著聚變研究的深入,基于EFIT 程序無法提供更準確的剖面分布,尤其是不同電流份額及分布(比如自舉電流)條件下的磁面重建結(jié)果,難以滿足高約束模式H 模,尤其是先進運行模式的深入研究需求.

    托卡馬克裝置DIII-D 發(fā)展了集成模擬框架程序OMFIT (one modeling framework for integrated tasks)[9-11].通過開發(fā)工作流,根據(jù)目的將不同程序進行集成運行,可以開展動力學平衡位形重建,磁流體穩(wěn)定性研究,高次諧波快波電流驅(qū)動優(yōu)化,磁流體和微觀湍流的相互作用以及自洽的平衡輸運模擬等物理內(nèi)容研究.DIII-D 托卡馬克裝置研究人員提出設想并使用該程序進行了實踐工作,充分顯示了OMFIT 框架在托卡馬克聚變研究的眾多領(lǐng)域的超強能力[12].

    HL-2A 裝置作為國內(nèi)實現(xiàn)偏濾器位形高約束模式(H 模)放電實驗的托卡馬克裝置,使中國成為歐美日之后成功實現(xiàn)高約束模式運行的國家[13,14].HL-2A 具有相對較高的加熱/驅(qū)動功率,開展了與ITER(國際熱核聚變實驗堆)物理相關(guān)研究.該研究基于HL-2A 高約束模式相關(guān)物理實驗,利用OMFIT 集成模擬平臺,搭建了包含EFIT,ONETWO[15,16],TGYRO[17,18]等核心程序的具備開展動力學平衡位形重建的工作流,通過實驗診斷數(shù)據(jù)給出的離子/電子的溫度/密度剖面,開展了高約束模式(H 模)放電實驗的集成模擬分析,重建了HL-2A高約束模式下等離子體的磁面位形,給出了準確的等離子體邊界,提供了完整的溫度演化、密度剖面以及安全因子剖面.基于以上物理參數(shù),本文分析了等離子體電流成分,分別給出歐姆電流、自舉電流、驅(qū)動電流占總電流的份額和剖面分布.基于OMFIT 模擬程序和HL-2A H 模實驗數(shù)據(jù),考慮物理模型自洽性,基于物理量之間的數(shù)據(jù)智能傳遞耦合,得到了托卡馬克等離子體電流剖面和磁面的重建.本研究結(jié)合HL-2A 的實驗結(jié)果、整合多個實驗反演程序和模擬程序進行集成建模分析(OMFIT),整合托卡馬克裝置大量的物理數(shù)據(jù),得到實時托卡馬克裝置物理整體剖面參數(shù)并進行分析.

    本文首先介紹理論分析模型和建立集成模擬程序OMFIT 框架及文中研究所需主要模塊,描述HL-2A 高約束放電集成模擬工作流設計,然后選取HL-2A 實驗高約束模式診斷物理結(jié)構(gòu)參數(shù),介紹了基于OMFIT 框架的整合,對實驗數(shù)據(jù)和動力學平衡反演程序計算,介紹了重建后的HL-2A磁面位形,電流(包含自舉電流份額)和安全因子q剖面,最后通過等離子體磁場和電流等剖面重建和實時監(jiān)測,評估了HL-2A 高約束放電實驗的品質(zhì).

    2 集成模擬模型

    2.1 理論模型

    作為綜合模擬平臺,集成模擬程序OMFIT 能夠用于各種物理程序之間的耦合組成復雜的流程模塊以分析物理問題,其特點在于其對托卡馬克實驗數(shù)據(jù)的實時整合.目前該平臺已經(jīng)整合了物理模擬程序,如EFIT,EPED,ONETWO和TGYRO 等.在自洽的平衡與輸運預測模型基礎(chǔ)上,利用OMFIT集成模擬平臺,搭建了包含EFIT,ONETWO,TGYRO 等核心程序的具備開展動力學平衡位形重建的建模過程,采用智能處理手段工作流針對HL-2A 高約束模式H 模放電進行了研究.能量源和能量損失用平衡反演程序ONETWO 可以計算,包括輔助加熱和輻射損耗以及離子和電子之間的能量交換.湍流輸運和新經(jīng)典輸運特性用模擬代碼TGYRO,TGLF[19]和NEO[20]分別計算,得到能量、粒子和環(huán)向角動量的輸運通量.TGYRO 通過匹配能量通量和目標通量來計算穩(wěn)態(tài)溫度分布,目標通量由相關(guān)的能量源和損失項的體積積分給出,磁場分布由二維(2-D)EFIT 平衡反演程序求解.

    本文采用固定邊界求解理想磁流體力學平衡Grad-Shafranov(G-S)方程來獲得平衡剖面演化特性.磁流體方程組為(1)式—(5)式,其中(1)式為粒子守恒方程,(2)式為動量守恒方程,(3)式為能量守恒方程,(4)式為等離子體在電磁場中的安培定律,(5)式為法拉第定律.(6)式為受力平衡條件.以上為理想條件下的磁流體MHD 方程組,用來描述磁約束等離子體中范圍較大的粒子運動狀態(tài).等離子體的每個組分的粒子處于平衡態(tài),粒子的速度分布函數(shù)基本符合麥克斯韋分布,可以描述等離子體在電磁場中的平衡分布:

    模擬條件假設為:磁流體重力相對于電磁力的約束很小,等離子體在磁場中達到平衡條件:?t v=0,(2) 式可以寫成(7)式,假設等離子體的宏觀運動速度為0 m/s,即v=0 m/s,寫成(8)式,分別用J和B來乘等式兩邊,可以得到(10)式:

    磁場的洛倫茲力和熱壓力之間的平衡可以從(8) 式得到,在輸運的尺度上作為力的平衡方程,包含電流密度、壓力和電磁力之間的作用.在時間尺度上作為與時間無關(guān)的平均速度為零的平衡理想磁流體方程的解.通過求解該方程可以得到等離子體的位形(包含壓強、磁場、粒子密度、溫度等)以及約束等離子體的外場(磁場,電場)隨空間坐標變化.由 (9) 式可以得到壓強相等的磁面輪廓,由 (10) 式可得到磁面上的電流密度及徑向電流導致的電磁場力平衡方程.從動量守恒方程出發(fā),假設等離子體的宏觀運動速度趨近于平衡,推導出理想條件下的MHD 平衡方程.考慮壓力張量是非負的,得到的平衡物理量能描述線性理想磁流體不穩(wěn)定性.作為壓力、電磁作用和電流密度之間的平衡的解,從中能判斷出約束等離子體的外場條件和絕對等離子體的形狀.在柱坐標系f(R,φ,Z)中,用函數(shù)F(F=RBf)和磁通函數(shù)ψ來表達磁場B.φ作為大環(huán)方向的角度,ψ作為等離子體的極向磁通被準確的用來描述磁面,其作為封閉環(huán)形磁面內(nèi)的極向磁通量的總和.在磁約束托卡馬克等離子體,圓截面的等離子體方程通常采用柱坐標系,方向之間有如(11)式 的關(guān)系.R為徑向參數(shù),φ為角向參數(shù),Z為軸向參數(shù).由磁場B無旋性,如(12)式所示.考慮磁流體的對稱性,那么在軸對稱平衡中,將 (11) 式兩邊同時求散度.由于軸對稱平衡中 0=?·(?ψ×?φ),(13) 式變?yōu)?14)式:

    考慮到對稱性,標量函數(shù)F和ψ都與環(huán)向角度φ無關(guān),作為Z和R的函數(shù).將求解磁場B的(11) 式兩邊同時求旋度,考慮?×(?T)=0和?×(?×t)=?(?·t)-?2t,(15) 式可以推導出(16)式,代入安培定律,得到等離子體電流密度J.從(8)式—(10)式可知,由于磁面是平衡的等壓強磁面,等離子體壓強P只和極向磁通函數(shù)ψ有關(guān),即P=P(ψ).由于P的函數(shù)只考慮ψ,由(10)式得到(18)式,將(17)式中的等離子體電流密度代入(18)式,于是得到(19)式.通過前面的推導,已經(jīng)得到J和B的式,代入MHD 式中.通過求解磁流體平衡方程,可以得到等離子體壓強與等離子體電流分布,得到平衡位形參大小數(shù),可在托卡馬克裝置中描述等離子體平衡電流位形和磁場位形分布:

    將(19)式和(20)式代入(8)式中可得:

    其中(?φ·?φ=R-2).比較等式(19),得到MHD磁流體平衡方程:

    2.2 HL-2A 高約束放電集成模擬

    基于磁流體平衡方程(23),將HL-2A 初始條件和邊界條件下,電子溫度、離子溫度和電子密度等進行擬合處理,代入模擬程序,具體的工作流如圖1所示.采用EFIT 程序計算得到HL-2A 平衡位形重建的平衡物理參數(shù),代入集成模擬程序,將實驗診斷數(shù)據(jù)(電子溫度Te,離子溫度Ti,電子密度ne和旋轉(zhuǎn)信息Vr)導入動力學剖面計算程序,計算外部驅(qū)動電流、歐姆電流和自舉電流,給出演化電流剖面信息,計算源和損失項.源項包括粒子源、能量源、外部電流驅(qū)動電源和動量源,損失項為輻射損失.粒子源主要來自中性束注入氫及同位素、充氣和壁材料濺射.能量源主要包括輔助加熱,例如中性束和波加熱.動量源主要來自非對稱中性束注入引發(fā)的等離子體旋轉(zhuǎn).最終得到源和損失項的參數(shù)以及新的平衡參數(shù)、壓強信息與所有電流份額組成的平衡電流剖面信息.

    圖1 HL-2A 高約束模式(H 模)放電模擬循環(huán)示意圖Fig.1.Design and analysis of HL-2A high confinement model (H model) discharge integrated simulation workflow.

    當輸出端存在電流演化時,利用EFIT 反演平衡參數(shù)剖面來更新的平衡參數(shù)剖面,即壓力梯度和極向電流梯度,得到新的平衡參數(shù),包含壓強信息與所有電流份額組成的平衡電流剖面信息,代入到EFIT 程序,求解G-S 方程,得到自洽的二維平衡方程.將自洽的平衡和源項的參數(shù)輸入到計算程序,利用源和損失項以及更新的平衡參數(shù),進行動理學剖面的演化,采取臺基區(qū)固定,演化芯部的動力學剖面,得到新的芯部的動理學剖面,繼而將新平衡動力學剖面(電子溫度Te,離子溫度Ti,電子密度ne和旋轉(zhuǎn)信息Vr)輸入到程序中,進行迭代,直到相鄰迭代之間的誤差最小,表明得到了收斂自洽的解.通過比較實驗數(shù)據(jù)與模擬的結(jié)果,完成HL-2A 磁面和安全因子與電流剖面重建.

    2.3 HL-2A 裝置參數(shù)

    本模擬采用HL-2A 幾何位形和物理參數(shù):大半徑R=1.65 m、小半徑a=0.4 m、環(huán)向磁場BT=2.8 T、環(huán)向等離子體電流IP=480 kA,具有偏濾器位形,圓截面,H 模放電實驗.HL-2A 具有相對較高的加熱功率,其電流加熱系統(tǒng)和輔助加熱包括歐姆加熱、離子回旋共振加熱、中性束加熱和電子回旋共振加熱以及低雜波系統(tǒng).HL-2A 裝置的電子回旋系統(tǒng)總規(guī)模為5 MW,包括68 GHz/3 MW/1 s 的一次諧波系統(tǒng)和140 GHz/2 MW/3 s 的二次諧波系統(tǒng).HL-2A 中性束注入系統(tǒng)注入功率為1.1 MW,束能量為44—55 keV、脈沖寬度為2 s,包括4 套大功率離子源、中性束線部件、中性束電源和診斷測量系統(tǒng)[22].

    基于HL-2A 實驗參數(shù)進行先進運行模式(高約束模式,H 模)動力學平衡位形重建,將實驗測量的離子和電子溫度、離子密度、旋轉(zhuǎn)角速度剖面參數(shù)作為平衡物理量導入到EFIT 平衡反演程序中.電子密度和溫度,離子溫度和旋轉(zhuǎn)角速度的參數(shù)剖面輸入將是集成模擬分析重點.在對放電參數(shù)選取的時候,選擇了高約束模式(H 模)的電子密度和溫度、離子溫度,通過數(shù)據(jù)庫選擇臺基明顯的H 模放電參數(shù),以及輸運壘比較明顯的物理參數(shù),即#37012 炮放電參數(shù):等離子體電流為153799 A,磁場為1.3 T,等離子體位置坐標小半徑a=0.376 m,大半徑R=1.53 m,采用兩束中性束加熱,加熱功率分別為600 keV和725 keV.該次放電炮號的放電情況如圖2所示,在882 s 后出現(xiàn)了高約束H模.加熱參數(shù)(加熱功率、電子溫度、離子溫度等)通過相關(guān)的計算得到:比壓βn=2.3,離子溫度Ti=1.58 keV,電子溫度Te=0.51 keV.設置好集成模擬程序初始等離子體位形重建物理量和邊界等離子體參數(shù)后,就可以進入到反演電流剖面的計算過程.

    圖2 HL-2A 炮號#37012 放電參數(shù)Fig.2.Shot #37012 parameters of HL-2A discharge.

    3 分析和討論

    3.1 HL-2A 實驗數(shù)據(jù)

    本文選取托卡馬克裝置HL-2A 參數(shù)#37012放電參數(shù)如下:模擬采用初始放電脈沖的平頂段處的電子溫度、離子溫度以及電子密度.從HL-2A 導入的放電數(shù)據(jù)繪制的曲線如圖3所示,通過將高能中性原子注入來將電子和離子進行加熱,其中離子加熱為主,離子溫度高于電子溫度,選取芯部電子和離子溫度(大半徑R=1.65 m),對噪聲進行平滑處理.電子密度的分布是等離子體芯部區(qū)域(ρ=0)到邊緣區(qū)域(ρ=±1,+1 代表右側(cè)弱場側(cè),—1 代表左側(cè)強場側(cè))的數(shù)據(jù),在處理時采用右半部分弱場側(cè)(0 <ρ< 1)的參數(shù).通過建立絕對坐標的位置和磁軸的位置相對應變換關(guān)系,得到磁面的位置來確定邊緣位置,分別對應歸一化后的0 <ρ< 1位置.選取離散的數(shù)據(jù)集來構(gòu)造一個解析函數(shù),讓原來的離散函數(shù)盡可能接近擬合曲線來對實驗數(shù)據(jù)進行擬合,通過迭代返回在x(磁零點)處的結(jié)果,同時在0<ρ<1 平均生成離散點,將數(shù)據(jù)進行3 次樣條數(shù)據(jù)插值,流程圖如圖4所示.能量源和能量損失采用ONETWO 程序計算,將電子溫度、離子溫度與電子密度剖面導入ONETWO模塊.TGYRO 由ONETWO和EFIT 計算所得到的已歸一化的參數(shù)剖面進行迭代,計算離子溫度和電子溫度、電子密度和等離子體環(huán)向旋轉(zhuǎn)速度及安全因子q的穩(wěn)態(tài)物理參數(shù)的剖面,從而計算芯部的輸運通量.計算選取的HL-2A 的設計參數(shù)為:βn=2.4,2.3;ne=2.95×1019m—3;Ti=1.53,1.58 keV;Te=0.53,0.51 keV;PNBI1=600 kW;PNBI2=725 kW;fbs/fohm=34.6/53.3;fbeam=12.1%;a=0.376 m;R=1.667 m;Ip=153799 A;B=1.3 T.

    圖3 HL-2A 托卡馬克平衡物理參數(shù)實驗數(shù)據(jù) (a)等離子體中離子溫度Ti徑向分布;(b)電子溫度Te徑向分布;(c)電子密度ne徑向分布Fig.3.(a) Radial distribution profile of ion temperature Tiin HL-2A Tokamak plasma;(b) radial distribution profile of electron temperature Tein HL-2A Tokamak plasma;(c) radial distribution profile of electron density nein HL-2A tokamak plasma.

    圖4 HL-2A 實驗數(shù)據(jù)處理,物理參數(shù)擬合徑向分布和OMFIT 耦合流程Fig.4.Flow chart for getting fitting curve processing by using discrete and interpolation of HL-2A experimental data.

    3.2 結(jié)果和討論

    基于HL-2A 實驗數(shù)據(jù)放電數(shù)據(jù),進行實驗擬合并將等離子體參數(shù)耦合到OMFIT 集成模擬平臺,求解MHD 磁流體平衡方程,計算電流密度(帶電粒子)和能量的輸運信息,建立邊緣和芯部溫度徑向分布函數(shù),給出邊緣溫度和計算出的臺基溫度剖面和電子密度徑向分布.芯部的新經(jīng)典輸運和湍流輸運由動理學剖面的演化計算,通過得到芯部的動理學剖面,固定臺基區(qū)的參數(shù)進行迭代獲得.

    圖5給出了重建后的電子溫度、離子溫度和電子密度,通過集成模擬重建電流剖面曲線與安全因子參數(shù)的剖面.通過整合等離子體離子溫度,電子溫度和密度,等離子體旋轉(zhuǎn)角速度,移動斯塔克效應測量的法拉第旋轉(zhuǎn)角等參數(shù),重建先進的高約束模式下等離子體整體輪廓參數(shù)(壓強、電流密度、安全因子)演化物理圖像.

    圖5 采用集成模擬程序OMFIT 重建的 (a) HL-2A 的電子溫度Te和離子溫度Ti剖面;(b)電子密度ne剖面Fig.5.Reconstruction of HL-2A parameters by OMFIT integrated simulation code:(a) Electron temperature Teand ion temperature Tiprofiles;(b) electron density ne profile.

    圖6所示為穩(wěn)態(tài)等離子體參數(shù)剖面的實驗擬合與集成模擬結(jié)果比較.當達到穩(wěn)態(tài)時,模擬得到的電子溫度、離子溫度、電子密度和旋轉(zhuǎn)因子,與#37012 次放電采用的集成模擬程序OMFIT 重建的該參數(shù)剖面進行對比.在芯部(ρ <0.8)區(qū)域,湍流輸運和新經(jīng)典輸運采用TGYRO 計算得到,邊界附近處(0.8<ρ<1.0)的溫度采取實驗數(shù)據(jù)作為初始值進行迭代計算,芯部的電子溫度初始值為0.6 keV,芯部的離子溫度初始值為1.2 keV.從圖6可以看到,通過模擬計算得到的電子和離子溫度空間分布、電子密度和旋轉(zhuǎn)速度等模擬結(jié)果與實驗曲線在演化過程中有少許差異,實驗數(shù)據(jù)擬合誤差及模擬模型的近似有關(guān),但整體趨勢基本一致.

    圖6 HL-2A 托卡馬克等離子體電子密度(a)、電子溫度(b)、離子溫度(c)和旋轉(zhuǎn)角速度(d)隨徑向剖面的集成模擬和實驗擬合結(jié)果對比Fig.6.HL-2A Tokamak electron density (a),electron temperature(b),ion temperature (c),rotational angular velocity,and(d) versus radial profile by OMFIT and experimental fitting in shot #37012.

    圖7所示為HL-2A 托卡馬克磁面位形、壓強分布、安全因子q、壓強梯度和極向電流梯度等物理參數(shù)剖面.將剖面參數(shù)導入集成建模工作流中,通過輸運程序ONETWO和TGYRO 與平衡程反演序EFIT,迭代得到重建平衡參數(shù)剖面。通過反演平衡結(jié)果,可以重建電流密度與安全因子剖面.采用EFIT 平衡反演程序,結(jié)合偏振干涉儀測量的法拉第旋轉(zhuǎn)角(pitch angele),反演獲得安全因子q分布,最終計算電流密度分布(Jdistributtion).從圖7可以看到,壓強梯度在芯部(ρ=0)下降,后又上升,在ρ=0.7 附近有一個最大值,對應著較強的內(nèi)部輸運壘(internal transport barrier,ITB),內(nèi)部輸運壘的形成有利于高約束模式的維持,起到凈空雜質(zhì)的作用,相當于在芯部和邊緣之間建立了一個密度和溫度梯度勢壘.

    圖7 集成模擬程序OMFIT 得到的HL-2A 托卡馬克等離子體磁場位形結(jié)構(gòu)圖 (a);等離子體壓強(單位Pa)徑向剖面分布(b);壓強梯度P’(無量綱)徑向分布(c);安全因子q(無量綱)隨著徑向變化分布剖面(d);極向電流梯度FF′(A/m3)徑向分布剖面(e)Fig.7.HL-2A Tokamak plasma magnetic field configuration obtained by integrated simulation code OMFIT (a);pressure profile(unit Pa) (b);profile distribution of the pressure gradient(Dimensionless quantity) (c);profile distribution of safety factor q(Dimensionless quantity) obtained (d);profile distribution of the polar current gradient FF′(unit A/m3)(e).

    HL-2A 較強的內(nèi)部輸運壘所形成的臺基區(qū)(pedestal region)[21-22],可以從離子溫度0.7 <ρ< 0.9之間形成的密度快速下降區(qū)域得到(如圖5所示).磁場的x點在z=—0.4 m 附近,該結(jié)果有利于設計偏濾器的位形,引導最后一個閉合磁面的帶電粒子流動來轟擊偏濾器靶板.磁面的位形接近圓形,壓力圖在邊緣附近有梯度的變化.由實驗結(jié)果圖5(a)(b)可以看到,與低約束模式(L 模)對比,邊緣密度、溫度剖面變陡相當于壓強剖面中增加了一個臺階,因而在同樣的加熱功率下,總儲能增大,能量約束時間隨之延長.約束改善源于局域性的反常輸運的減弱,即輸運壘的建立,或邊緣附近湍流輸運的減弱,使邊緣區(qū)的密度、電子和離子溫度剖面都變陡.

    高約束模式(H 模)首次在德國托卡馬克裝置ASDEX 上發(fā)現(xiàn),具有重要的物理意義[23].高約束模式(H 模式)相對于低約束模式(L 模),粒子的約束密度極限和能量約束時間都有非常大的提升.當能量連續(xù)注入到等離子體中時,邊緣剪切流會在注入功率達到某一閾值時增強,然后在邊界區(qū)域形成輸運壘,等離子體就從L 模轉(zhuǎn)換成H 模.國際上其他的有較大加熱功率的聚變裝置都發(fā)現(xiàn)了H 模約束,在球馬克和仿星器上也都有發(fā)現(xiàn),說明H 模的出現(xiàn)與偏濾器或限制器的位形和加熱方式無關(guān)[24],并且與約束裝置類型無關(guān).

    當?shù)入x子體達到H 模以后,邊界的溫度和密度不斷上升,從而形成了稱為臺基的結(jié)構(gòu),即邊界輸運壘(ETB),此刻的等離子體的能量約束水平已大大提高,約束性能也大幅提高.文獻給出了ITER臺基區(qū)的預測圖[25],可得到臺基區(qū)的壓強和密度溫度預測,臺基結(jié)構(gòu)處的計算可以用修正tanh 函數(shù)[26,27]來擬合:

    其中,XSYM為臺基中心的徑向坐標,X為臺基區(qū)密度的徑向坐標,w是臺基的寬度,a是小半徑,為0.376 m,電子密度ne=2.95×1019m-3,A+B為臺基的高度[28].如圖3(b)所示,HL-2A 的該炮號產(chǎn)生的臺基區(qū)的寬度在0.8 <ρ< 1.0,臺基寬度約w=0.376 m× 0.2=7.52 cm.

    托卡馬克聚變裝置中臺基的寬度一般僅僅為幾毫米到幾厘米,但對等離子體的約束有及其大的作用.在等離子體的芯部隨著壓力梯度的持續(xù)增長達到某一臨界值的時候,湍流運輸和各種不穩(wěn)定模式(例如捕獲電子模TEM)都會出現(xiàn)閾值效應的爆發(fā)導致梯度變平,然后芯部等離子體的相關(guān)參數(shù)不能無限的增長.等離子體約束性能隨著臺基高度增大而加強,使得等離子體的聚變功率提高[29].

    H 模具有較長的能量約束時間,較高的自舉電流份額和較大的粒子約束密度,被視為未來ITER實驗聚變堆的基本運行模式.研究人員對等離子體的能量約束時間和相關(guān)的參數(shù)進行了許多實驗[30],在L 模下的約束時間可以用(26)式定標率來估計:

    在H 模具有高臺基時,其H 模約束下約束時間可根據(jù)(27)式來估計:

    式中,Ip為等離子體電流,BT為環(huán)向磁場強度,M為原子質(zhì)量,R為大半徑,HL-2A 的大半徑為1.65 m,a為小半徑,HL-2A 的小半徑為0.376 m,κ為等離子體拉長比,Pt是等離子體得到的凈的加熱功率.有了ITER98 的定標率之后[31],可以推導出能量約束增強因子,為托卡馬克裝置獲得能量的約束時間和定標率所預測的能量時間之間的比例,體現(xiàn)約束狀態(tài)的好壞.

    圖8所示為采用集成模擬得到的電流密度剖面(電流份額)與安全因子q剖面.從圖8可以看出,安全因子q最小值為1.8,該值較低可能引發(fā)危險性的低nMHD 模式,從而導致因為MHD 不穩(wěn)定性的放電破裂.圖8(a)給出了中性束(beam)電流,自舉電流(BS),歐姆電流(Ohmic),射頻電流(Rf)和全電流(Tot)剖面分布,中性束加熱效果良好,可以看出在中心區(qū)(0 <ρ< 0.4)所占的份額最多,該區(qū)域中性束電流密度增大,中性束電流密度在中心區(qū)強所占份額最多,隨著徑向分布演化,在邊緣區(qū)域(0.7 <ρ< 1.0)逐漸減弱到趨近0,同時自舉電流密度(BS)份額大幅增大,而托卡馬克穩(wěn)態(tài)運行需要大幅提高自舉電流的份額.通過計算得到的電流份額和空間分布可以對HL-2A 實驗控制和等離子體位形設計提供參考.

    圖8 (a)磁約束托卡馬克裝置shot #37012 由集成模擬程序OMFIT 計算的整體電流密度、包含中性束電流密度、自舉電流密度、歐姆電流密度和射頻電流密度及份額分布曲線;(b)OMFIT 程序和平衡反演程序EFIT 計算的安全因子q 剖面對比圖Fig.8.(a) Magnetic confinement Tokamak device shot #37012 total current densities calculated by the integrated simulation program OMFIT,including neutral beam current density,bootstrap current density,Ohmic current density,and RF current density with percentages of respective current densities in total current densities;(b) comparison of the safety factor q profiles calculated by the OMFIT program and the balance inversion code EFIT.

    通過自舉電流分布可知,在0.2 <ρ< 1.0 區(qū)域,自舉電流在整體電流的份額在上升,而歐姆電流和中性束的電流份額在下降,對比電流驅(qū)動剖面與安全因子q剖面發(fā)現(xiàn),自舉電流的改變與q的剖面的梯度的變化有聯(lián)系.總電流由歐姆電流Iohm、自舉電流IBS、中性束驅(qū)動電流INBI和射頻感應電流IRF組成.歐姆電流、自舉電流、中性束驅(qū)動電流和射頻感應電流份額通過計算分別為53.2%,34.6%,12.0%和0.2%.通過與圖3對比,發(fā)現(xiàn)中性束NBI注入以后,在1.6 m<R<1.7 m 處離子溫度上升明顯,證明加熱芯部效果良好,中性束粒子在注入過程中慢化的過程更加接近芯部.電子密度受到離子加熱的影響,在該區(qū)域溫度也快速上升,形成臺基區(qū)的內(nèi)部輸運壘位形(ITB).

    從圖7(c)壓強梯度隨徑向變化規(guī)律可以看到,中性束注入初期,電子溫度很快升高,形成一個中空的電流分布,后期的高功率NBI 使等離子體壓強梯度迅速加大,中心區(qū)的粒子密度很快提高,這兩種效應都使自舉電流比例加大,負磁剪切位形得以長時間維持,這種運行模式能同時具有中心負磁剪切和非常峰化的密度分布,具有較長的約束時間和較高的約束因子H98.在0.1 <ρ< 0.2 區(qū)間壓強梯度方向發(fā)生變化(由正變負),邊緣溫度的增加使等離子體壓強梯度在邊緣區(qū)不斷增加,如圖7(c)所示,溫度分布在邊緣區(qū)變陡,極向旋轉(zhuǎn)速度在此區(qū)域明顯加大,于是在此區(qū)域形成邊緣輸運壘ETB(edge transport barrier).在穩(wěn)態(tài)托卡馬克JT-60U,TFTR和DIII-D 都觀測到約束改善模式H 模式和各種負剪切位形放電.內(nèi)部輸運壘位形(ITB)均是以負磁剪切位形為基礎(chǔ)的改善約束模式,負磁剪切位形是在等離子體中心形成負的磁剪切區(qū),安全因子q的分布函數(shù)不再是單調(diào)下降分布,具有高的自舉電流份額,如圖8(a)所示,自舉電流大幅提升.在電流份額中,自舉電流所占份額對未來自持穩(wěn)態(tài)放電至關(guān)重要.一般可通過中性束加熱和低雜波驅(qū)動等手段,大幅增加自舉電流份額.

    穩(wěn)定狀態(tài)的托卡馬克等離子體中的可能產(chǎn)生扭曲模不穩(wěn)定性,需有環(huán)向磁場BT和具有上限的等離子體電流.在環(huán)形等離子體中,這個極限電流可用等離子體邊界的安全因子q> 1 滿足扭曲不穩(wěn)定性的穩(wěn)定條件,q稱作安全因子,Bp是等離子體電流所產(chǎn)生的極向磁場.安全因子可由兩種磁場來計算:

    其中,r為位置坐標,可由r/a=ρ來計算,BT是環(huán)向磁場,Bp是極向磁場.

    通過圖8(b)中q剖面的變化,模擬得到的安全因子q剖面并非單調(diào)分布曲線.在0.2 <ρ< 0.9區(qū)域,通過實驗數(shù)據(jù)集成模擬得到的q剖面出現(xiàn)了波動,集成模擬得到的q最小值略小于2.0,在ρ=0處q≈ 2.0,在芯部略有波動.實驗擬合得到安全因子q分布單調(diào)遞增,中間無明顯漲落趨勢變化.集成模擬得到的安全因子q的剖面比單獨采用EFIT由磁診斷數(shù)據(jù)擬合反演得到的q的剖面更加趨于合理.EFIT 是利用外部的電磁診斷設備(磁通環(huán)和磁探針)獲取的磁信號,對實驗數(shù)據(jù)進行了全面、深入的分析,準確重建等離子體平衡位形和相關(guān)的形狀參數(shù).由于缺乏芯部診斷數(shù)據(jù),難以給出較為豐富和更加準確的電流剖面及壓強剖面分布.

    集成模擬OMFIT 是對芯部的數(shù)據(jù)進行計算,能較為豐富地給出電流剖面及壓強剖面分布,即集成模擬得到的q剖面比單獨采用EFIT 由磁診斷數(shù)據(jù)反演得到的安全因子q分布更加合理.由圖8(b)可知,EFIT 得到的q變化趨勢平緩,而結(jié)合等離子體參數(shù),輸運參數(shù)集成模擬得到的安全因子剖面q分布更加豐富.在0.3 <ρ< 0.7 區(qū)域,安全因子增加明顯.

    4 結(jié)論

    以EFIT,ONETWO,TGYRO 等核心模擬程序為基礎(chǔ),搭建了具備開展動力學平衡位形重建的集成模擬平臺OMFIT.針對HL-2A 先進高約束模式(H 模)放電實驗的數(shù)據(jù),采用集成建模手段,重建了HL-2A 放電實驗的磁面位形,給出準確的等離子體邊界物理量剖面,反演了完整的離子和電子溫度徑向分布、電子密度徑向分布、等離子體旋轉(zhuǎn)角速度徑向分布以及安全因子q徑向分布.

    基于實驗數(shù)據(jù)和模擬重建的結(jié)果,給出等離子體電流成分構(gòu)成,包括歐姆電流、自舉電流、驅(qū)動電流占總電流的份額,并給出了這些電流詳細剖面分布.基于集成模擬計算的HL-2A 高約束模式放電集成模型動力學剖面(電子溫度、離子溫度、電子密度、旋轉(zhuǎn))與實驗診斷擬合圖像基本一致,驗證了該集成模型.采用OMFIT 集成模擬平臺,整合了物理模擬和數(shù)據(jù)重建程序EFIT,ONETWO,TGYRO 等,得到豐富的q剖面以及電流剖面參數(shù).在托卡馬克等離子體運行過程中,對于安全因子q測量,研究人員通常采用運動斯塔克效應光譜診斷來獲取法拉第旋轉(zhuǎn)角和徑向電場分布.但是對于可見光譜診斷,面向芯部的光譜測量遇到了芯部等離子體溫度高,原子被完全電離,光譜信號弱,邊緣等離子體受到來自偏濾器和外壁內(nèi)壁的反射雜散光干擾,因此獲取的安全因子圖像誤差較大,這對穩(wěn)態(tài)托卡馬克反饋控制安全有潛在風險.從集成模擬思路出發(fā),基于HL-2A 實驗數(shù)據(jù)和物理模擬模型整合,重建了非單調(diào)變化安全因子q的徑向分布,比單獨通過EFIT 平衡反演程序得到的單調(diào)變化安全因子徑向分布更為豐富.下一步計劃對集成模擬程序部分模塊進行實時耦合和自洽性驗證,得到HL-2M 相應的等離子體參數(shù)的實時剖面信息演化,實現(xiàn)對穩(wěn)態(tài)運行實時反饋控制.

    猜你喜歡
    位形托卡馬克等離子體
    中間支撐剛度對雙跨梁屈曲穩(wěn)定性的影響
    英原型聚變堆ST40實現(xiàn)1億℃等離子體高溫
    國外核新聞(2022年4期)2022-02-08 14:20:36
    連續(xù)磁活動對等離子體層演化的影響
    基于低溫等離子體修飾的PET/PVC浮選分離
    等離子體種子處理技術(shù)介紹
    基于旋量理論的四自由度抓取機械手奇異位形分析
    EAST托卡馬克上截面效應對電荷交換復合光譜測量結(jié)果的影響
    核技術(shù)(2016年4期)2016-08-22 09:05:30
    基于可操作度的機器人最優(yōu)初始位形研究
    大眾科技(2015年11期)2015-11-24 01:57:16
    基于最優(yōu)初始位形的冗余度機器人可操作度優(yōu)化*
    反射內(nèi)存網(wǎng)絡在托卡馬克裝置快控制器中的應用
    久久精品人人爽人人爽视色| 真人一进一出gif抽搐免费| 91成年电影在线观看| 在线天堂中文资源库| 长腿黑丝高跟| 国产av一区二区精品久久| 久热这里只有精品99| 国产av一区二区精品久久| 免费高清在线观看日韩| 99精品久久久久人妻精品| 欧美成人免费av一区二区三区| 中文字幕高清在线视频| 欧美一区二区精品小视频在线| 亚洲欧美日韩无卡精品| av电影中文网址| 国产精品日韩av在线免费观看 | 亚洲狠狠婷婷综合久久图片| 性色av乱码一区二区三区2| 久久久久久久精品吃奶| 国产精品久久久久久人妻精品电影| 日韩av在线大香蕉| 精品熟女少妇八av免费久了| 一级作爱视频免费观看| 久久香蕉激情| 一边摸一边抽搐一进一小说| a级毛片在线看网站| 国产精品久久久久久人妻精品电影| 少妇裸体淫交视频免费看高清 | 精品欧美一区二区三区在线| 欧美日韩亚洲国产一区二区在线观看| 欧美另类亚洲清纯唯美| 日本在线视频免费播放| av天堂久久9| 99精品在免费线老司机午夜| 色播在线永久视频| 又黄又爽又免费观看的视频| 欧美精品亚洲一区二区| 国产蜜桃级精品一区二区三区| 俄罗斯特黄特色一大片| 在线av久久热| 欧美国产精品va在线观看不卡| av天堂在线播放| 欧美激情极品国产一区二区三区| 欧美色视频一区免费| 九色国产91popny在线| 老司机在亚洲福利影院| 国产免费av片在线观看野外av| 亚洲精品中文字幕一二三四区| 桃色一区二区三区在线观看| 超碰成人久久| 日韩 欧美 亚洲 中文字幕| 日本免费a在线| 黄色片一级片一级黄色片| 亚洲欧洲精品一区二区精品久久久| 又黄又爽又免费观看的视频| 日本 av在线| 久久精品国产亚洲av高清一级| 久久久久久亚洲精品国产蜜桃av| 免费在线观看视频国产中文字幕亚洲| 一二三四在线观看免费中文在| 手机成人av网站| videosex国产| 亚洲黑人精品在线| 国产精品 国内视频| 99精品久久久久人妻精品| 99久久精品国产亚洲精品| 91老司机精品| 老司机在亚洲福利影院| 国产亚洲精品第一综合不卡| 欧美日本中文国产一区发布| 满18在线观看网站| 最近最新免费中文字幕在线| 69精品国产乱码久久久| 国产极品粉嫩免费观看在线| 精品国产一区二区久久| 亚洲成a人片在线一区二区| 精品午夜福利视频在线观看一区| 国产高清videossex| 91大片在线观看| 国产精品一区二区免费欧美| 久久 成人 亚洲| 欧美乱码精品一区二区三区| 岛国在线观看网站| 国产精品综合久久久久久久免费 | 别揉我奶头~嗯~啊~动态视频| 亚洲第一电影网av| 不卡一级毛片| 天堂√8在线中文| 精品欧美国产一区二区三| 欧美久久黑人一区二区| 999久久久国产精品视频| 啦啦啦 在线观看视频| 99精品久久久久人妻精品| 日日夜夜操网爽| 999久久久国产精品视频| 成人免费观看视频高清| 少妇的丰满在线观看| 老司机深夜福利视频在线观看| 91成人精品电影| 国产精品精品国产色婷婷| 首页视频小说图片口味搜索| 窝窝影院91人妻| 亚洲午夜精品一区,二区,三区| 亚洲精品在线观看二区| 日韩欧美三级三区| av福利片在线| 精品高清国产在线一区| 婷婷精品国产亚洲av在线| 大陆偷拍与自拍| 亚洲一区高清亚洲精品| www日本在线高清视频| 国产xxxxx性猛交| 国产成人精品无人区| 日本三级黄在线观看| 国产精品乱码一区二三区的特点 | 亚洲国产中文字幕在线视频| 日韩中文字幕欧美一区二区| 久久精品影院6| av超薄肉色丝袜交足视频| bbb黄色大片| 日韩欧美国产在线观看| 最近最新中文字幕大全电影3 | 国产精品九九99| xxx96com| 曰老女人黄片| 婷婷精品国产亚洲av在线| 宅男免费午夜| 色综合站精品国产| 亚洲无线在线观看| 亚洲性夜色夜夜综合| 国产精品久久视频播放| 非洲黑人性xxxx精品又粗又长| 很黄的视频免费| 91国产中文字幕| 亚洲人成伊人成综合网2020| 欧美国产日韩亚洲一区| 最新美女视频免费是黄的| 黄色毛片三级朝国网站| 亚洲va日本ⅴa欧美va伊人久久| 亚洲人成网站在线播放欧美日韩| 伊人久久大香线蕉亚洲五| 1024香蕉在线观看| 日韩视频一区二区在线观看| 午夜福利在线观看吧| 黄色片一级片一级黄色片| 制服人妻中文乱码| 欧美乱色亚洲激情| 亚洲精品在线美女| 男女下面进入的视频免费午夜 | 丁香欧美五月| 啦啦啦 在线观看视频| 99久久国产精品久久久| 亚洲色图av天堂| 99精品欧美一区二区三区四区| 一区在线观看完整版| 丝袜美足系列| 亚洲第一av免费看| 久久久久九九精品影院| 亚洲国产欧美网| 久久国产乱子伦精品免费另类| 久久久久久久久久久久大奶| 制服人妻中文乱码| 91成人精品电影| 亚洲aⅴ乱码一区二区在线播放 | 欧美中文综合在线视频| 国产高清视频在线播放一区| 性色av乱码一区二区三区2| 人妻久久中文字幕网| 侵犯人妻中文字幕一二三四区| 男女床上黄色一级片免费看| 久久久久久亚洲精品国产蜜桃av| 9热在线视频观看99| 丰满的人妻完整版| 久99久视频精品免费| 又黄又粗又硬又大视频| 亚洲视频免费观看视频| 非洲黑人性xxxx精品又粗又长| 无限看片的www在线观看| 日本 av在线| 国产亚洲av嫩草精品影院| 免费一级毛片在线播放高清视频 | 成人18禁在线播放| 国产精品影院久久| 免费看十八禁软件| 麻豆一二三区av精品| 波多野结衣av一区二区av| av超薄肉色丝袜交足视频| 精品久久久久久久久久免费视频| 免费少妇av软件| 久久精品成人免费网站| 亚洲午夜理论影院| 一进一出好大好爽视频| 国内精品久久久久久久电影| 久久国产精品影院| 90打野战视频偷拍视频| 国产精品一区二区在线不卡| 日本五十路高清| a级毛片在线看网站| 亚洲欧美日韩高清在线视频| 人人妻人人澡欧美一区二区 | 午夜视频精品福利| 免费观看人在逋| 欧美国产日韩亚洲一区| 午夜免费成人在线视频| 亚洲五月天丁香| 亚洲色图综合在线观看| 在线观看一区二区三区| 国产极品粉嫩免费观看在线| 麻豆成人av在线观看| АⅤ资源中文在线天堂| 欧美日韩亚洲国产一区二区在线观看| 老熟妇仑乱视频hdxx| 欧美日韩乱码在线| 黄色a级毛片大全视频| 国产野战对白在线观看| 青草久久国产| 搞女人的毛片| av有码第一页| 欧美国产日韩亚洲一区| 人人妻人人澡欧美一区二区 | 国产欧美日韩一区二区三区在线| 两个人视频免费观看高清| 50天的宝宝边吃奶边哭怎么回事| 美女午夜性视频免费| 国产亚洲av嫩草精品影院| 国产亚洲精品av在线| 久久亚洲精品不卡| 日本一区二区免费在线视频| 欧美精品亚洲一区二区| 国产精品爽爽va在线观看网站 | 亚洲情色 制服丝袜| 日韩欧美国产一区二区入口| a在线观看视频网站| 韩国av一区二区三区四区| 欧美黑人精品巨大| 国产区一区二久久| 欧美日韩黄片免| 国产伦人伦偷精品视频| 亚洲欧美日韩高清在线视频| 一区二区日韩欧美中文字幕| 亚洲人成伊人成综合网2020| 黄色 视频免费看| 男女之事视频高清在线观看| 国产1区2区3区精品| 久久国产精品男人的天堂亚洲| 禁无遮挡网站| 午夜激情av网站| 亚洲精品粉嫩美女一区| 亚洲欧洲精品一区二区精品久久久| 此物有八面人人有两片| 真人一进一出gif抽搐免费| 一本大道久久a久久精品| 亚洲男人天堂网一区| 久99久视频精品免费| 免费一级毛片在线播放高清视频 | 在线十欧美十亚洲十日本专区| 国产精品乱码一区二三区的特点 | 好男人在线观看高清免费视频 | 美女高潮喷水抽搐中文字幕| 国内毛片毛片毛片毛片毛片| 日本撒尿小便嘘嘘汇集6| 91精品国产国语对白视频| 欧美亚洲日本最大视频资源| 午夜免费观看网址| 免费人成视频x8x8入口观看| 日韩欧美在线二视频| 欧美精品啪啪一区二区三区| 国产成人系列免费观看| 免费看美女性在线毛片视频| 国产成人av教育| 日日夜夜操网爽| 纯流量卡能插随身wifi吗| 国产一区二区激情短视频| 国产极品粉嫩免费观看在线| 免费一级毛片在线播放高清视频 | 精品熟女少妇八av免费久了| 日日摸夜夜添夜夜添小说| 涩涩av久久男人的天堂| 中亚洲国语对白在线视频| 精品高清国产在线一区| 99精品欧美一区二区三区四区| 一卡2卡三卡四卡精品乱码亚洲| 每晚都被弄得嗷嗷叫到高潮| 亚洲五月婷婷丁香| 国产乱人伦免费视频| 丝袜美腿诱惑在线| 国产1区2区3区精品| 欧美老熟妇乱子伦牲交| 午夜成年电影在线免费观看| 亚洲第一av免费看| 最新在线观看一区二区三区| 国产高清视频在线播放一区| 亚洲色图 男人天堂 中文字幕| 成人永久免费在线观看视频| 免费高清在线观看日韩| 国产激情欧美一区二区| 午夜免费观看网址| 久久久久国产一级毛片高清牌| 伊人久久大香线蕉亚洲五| 日韩大码丰满熟妇| 国产精品1区2区在线观看.| 国产av精品麻豆| 怎么达到女性高潮| 亚洲,欧美精品.| 精品欧美一区二区三区在线| 精品久久久久久,| 欧美日韩一级在线毛片| 亚洲三区欧美一区| 中亚洲国语对白在线视频| 两性午夜刺激爽爽歪歪视频在线观看 | 亚洲男人的天堂狠狠| 日韩精品青青久久久久久| 国产色视频综合| 日韩精品中文字幕看吧| 一夜夜www| 国内精品久久久久精免费| 国产高清激情床上av| 亚洲成国产人片在线观看| 亚洲专区中文字幕在线| 如日韩欧美国产精品一区二区三区| 老司机在亚洲福利影院| 国产精品久久久久久人妻精品电影| 黑人操中国人逼视频| 亚洲熟妇中文字幕五十中出| av天堂久久9| 国产精品香港三级国产av潘金莲| 久久精品国产99精品国产亚洲性色 | 91麻豆精品激情在线观看国产| 狂野欧美激情性xxxx| 亚洲五月色婷婷综合| 变态另类丝袜制服| 久久人妻av系列| 国产av一区二区精品久久| 久久草成人影院| 涩涩av久久男人的天堂| 非洲黑人性xxxx精品又粗又长| 每晚都被弄得嗷嗷叫到高潮| 波多野结衣一区麻豆| 亚洲va日本ⅴa欧美va伊人久久| 久久狼人影院| 国产男靠女视频免费网站| 最近最新免费中文字幕在线| 久久 成人 亚洲| 久久国产精品男人的天堂亚洲| 一本久久中文字幕| 中文字幕av电影在线播放| 色在线成人网| 欧美成人一区二区免费高清观看 | 丰满人妻熟妇乱又伦精品不卡| 一区二区三区激情视频| 久久精品人人爽人人爽视色| 亚洲免费av在线视频| 欧美激情久久久久久爽电影 | 国产精品永久免费网站| 亚洲精品中文字幕在线视频| 国产精品爽爽va在线观看网站 | 国产精华一区二区三区| 麻豆一二三区av精品| 国产又色又爽无遮挡免费看| 在线观看免费日韩欧美大片| 色播在线永久视频| 亚洲激情在线av| 午夜老司机福利片| 久久香蕉精品热| 91av网站免费观看| 国产人伦9x9x在线观看| 久久精品国产亚洲av香蕉五月| 精品日产1卡2卡| 国产成人系列免费观看| 中出人妻视频一区二区| 午夜a级毛片| 91麻豆av在线| 亚洲片人在线观看| 国产精品永久免费网站| 国产免费av片在线观看野外av| 国产激情久久老熟女| 男人舔女人下体高潮全视频| 熟女少妇亚洲综合色aaa.| 免费少妇av软件| 成年女人毛片免费观看观看9| 精品人妻1区二区| 亚洲人成伊人成综合网2020| 男女午夜视频在线观看| 中出人妻视频一区二区| 午夜福利欧美成人| 国产精品 欧美亚洲| 男人舔女人下体高潮全视频| 国内久久婷婷六月综合欲色啪| 日韩精品免费视频一区二区三区| √禁漫天堂资源中文www| 午夜福利欧美成人| 亚洲国产中文字幕在线视频| 国产精品久久久人人做人人爽| 国产片内射在线| 婷婷六月久久综合丁香| 天堂影院成人在线观看| 99热只有精品国产| 窝窝影院91人妻| 啦啦啦 在线观看视频| 中文字幕av电影在线播放| 日韩三级视频一区二区三区| 色综合婷婷激情| 亚洲av电影在线进入| 午夜日韩欧美国产| 少妇裸体淫交视频免费看高清 | 91麻豆av在线| 久久狼人影院| av天堂久久9| 精品人妻在线不人妻| 午夜老司机福利片| 丝袜在线中文字幕| bbb黄色大片| 午夜福利高清视频| 免费在线观看视频国产中文字幕亚洲| 两个人免费观看高清视频| 精品久久蜜臀av无| 国产高清激情床上av| 国产高清videossex| 欧美午夜高清在线| 看免费av毛片| 亚洲av日韩精品久久久久久密| 丝袜美足系列| 免费观看人在逋| 日韩精品中文字幕看吧| 久久天躁狠狠躁夜夜2o2o| 精品少妇一区二区三区视频日本电影| 亚洲五月婷婷丁香| 999久久久国产精品视频| 午夜老司机福利片| 88av欧美| 久久久久久久久免费视频了| 美女高潮到喷水免费观看| 成人精品一区二区免费| 国产精品av久久久久免费| 亚洲一区二区三区不卡视频| 日韩av在线大香蕉| 国产亚洲av高清不卡| 成年人黄色毛片网站| 丝袜在线中文字幕| 非洲黑人性xxxx精品又粗又长| 91成人精品电影| 中出人妻视频一区二区| 精品国产美女av久久久久小说| 久久午夜综合久久蜜桃| 欧美激情极品国产一区二区三区| netflix在线观看网站| 日韩欧美免费精品| 精品福利观看| 久久九九热精品免费| 成人av一区二区三区在线看| 三级毛片av免费| 91成人精品电影| 成熟少妇高潮喷水视频| 亚洲视频免费观看视频| 国产精品野战在线观看| 丝袜美腿诱惑在线| 国产精品免费视频内射| 久久香蕉激情| 人人妻人人澡欧美一区二区 | 欧美黑人欧美精品刺激| 国内精品久久久久精免费| 国产一区二区激情短视频| av电影中文网址| 亚洲av熟女| 亚洲五月色婷婷综合| 欧美成狂野欧美在线观看| 亚洲全国av大片| 亚洲成a人片在线一区二区| 亚洲国产看品久久| 国产精品一区二区在线不卡| av有码第一页| 色婷婷久久久亚洲欧美| 亚洲人成电影观看| 日韩精品青青久久久久久| 国产精品亚洲av一区麻豆| 日本 av在线| 九色国产91popny在线| 国产亚洲欧美在线一区二区| 亚洲自偷自拍图片 自拍| a级毛片在线看网站| 激情视频va一区二区三区| 在线观看免费午夜福利视频| 久久国产乱子伦精品免费另类| 法律面前人人平等表现在哪些方面| 宅男免费午夜| 99久久99久久久精品蜜桃| 日本a在线网址| 黑人巨大精品欧美一区二区mp4| 亚洲色图 男人天堂 中文字幕| 人人妻,人人澡人人爽秒播| 人妻久久中文字幕网| 91麻豆精品激情在线观看国产| 成人欧美大片| 欧美日韩亚洲综合一区二区三区_| 九色亚洲精品在线播放| 麻豆国产av国片精品| 国产野战对白在线观看| svipshipincom国产片| 日本精品一区二区三区蜜桃| 亚洲av成人一区二区三| a级毛片在线看网站| 国产精品99久久99久久久不卡| 男女下面插进去视频免费观看| 国产免费男女视频| 天堂动漫精品| 欧美在线一区亚洲| 看片在线看免费视频| 女生性感内裤真人,穿戴方法视频| 热re99久久国产66热| 宅男免费午夜| 夜夜躁狠狠躁天天躁| АⅤ资源中文在线天堂| 精品国内亚洲2022精品成人| 深夜精品福利| 国产成人一区二区三区免费视频网站| 亚洲精品一区av在线观看| 中文字幕人成人乱码亚洲影| 欧美黄色淫秽网站| 精品久久久久久,| 69精品国产乱码久久久| 伊人久久大香线蕉亚洲五| 亚洲成a人片在线一区二区| 在线视频色国产色| 一卡2卡三卡四卡精品乱码亚洲| 亚洲国产精品999在线| 99在线视频只有这里精品首页| 国产精品综合久久久久久久免费 | 欧美成人免费av一区二区三区| 久久人人精品亚洲av| 欧美一级毛片孕妇| 久久婷婷人人爽人人干人人爱 | 午夜福利一区二区在线看| 色尼玛亚洲综合影院| 久久婷婷人人爽人人干人人爱 | 亚洲视频免费观看视频| 正在播放国产对白刺激| 在线永久观看黄色视频| 国产av一区在线观看免费| 午夜福利一区二区在线看| 日本在线视频免费播放| 久久天堂一区二区三区四区| √禁漫天堂资源中文www| 国产一卡二卡三卡精品| 久久影院123| 国产精品,欧美在线| 成人18禁在线播放| 色综合亚洲欧美另类图片| 成人18禁高潮啪啪吃奶动态图| 国产蜜桃级精品一区二区三区| 正在播放国产对白刺激| 欧美午夜高清在线| 人人澡人人妻人| 久久久水蜜桃国产精品网| 丁香六月欧美| 90打野战视频偷拍视频| 99国产精品一区二区三区| 免费高清视频大片| 91精品国产国语对白视频| 亚洲九九香蕉| 久久精品亚洲熟妇少妇任你| 满18在线观看网站| 亚洲欧美日韩另类电影网站| 母亲3免费完整高清在线观看| 一边摸一边做爽爽视频免费| 国产色视频综合| 亚洲伊人色综图| 婷婷丁香在线五月| 丰满的人妻完整版| 人妻丰满熟妇av一区二区三区| 天堂√8在线中文| 99久久99久久久精品蜜桃| 亚洲精品中文字幕一二三四区| 欧美精品亚洲一区二区| 欧美乱妇无乱码| 免费在线观看影片大全网站| 极品人妻少妇av视频| 男女午夜视频在线观看| 成人三级做爰电影| 国产精品,欧美在线| 黄色a级毛片大全视频| 国产精品九九99| 丝袜美足系列| 日本免费a在线| 亚洲国产精品合色在线| av有码第一页| 国产伦一二天堂av在线观看| 亚洲七黄色美女视频| 夜夜看夜夜爽夜夜摸| 女警被强在线播放| 日本欧美视频一区| 高清在线国产一区| 男女床上黄色一级片免费看| 亚洲 欧美一区二区三区| 亚洲av日韩精品久久久久久密| 亚洲中文日韩欧美视频| 成人18禁高潮啪啪吃奶动态图| 88av欧美| 中文字幕高清在线视频| 91av网站免费观看| 九色亚洲精品在线播放| 后天国语完整版免费观看| 国产野战对白在线观看| 欧美成人免费av一区二区三区| 午夜福利18| 亚洲第一欧美日韩一区二区三区| 久久亚洲精品不卡| 亚洲欧洲精品一区二区精品久久久| 国产欧美日韩精品亚洲av| 女人爽到高潮嗷嗷叫在线视频| 亚洲欧美一区二区三区黑人| 在线观看66精品国产| 国产精品久久视频播放| 悠悠久久av| 国产区一区二久久| 日韩欧美一区视频在线观看|