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

    電磁驅(qū)動高能量密度動力學實驗的一維磁流體力學多物理場數(shù)值模擬平臺:SSS-MHD*

    2023-11-07 11:24:08孫承緯趙繼波羅斌強谷卓偉王桂吉張旭平陳學秒周中玉張紅平王剛?cè)A孫奇志文尚剛譚福利趙劍衡莫建軍蔡進濤金云聲趙小明劉倉理
    爆炸與沖擊 2023年10期
    關(guān)鍵詞:構(gòu)形算例方程組

    孫承緯,陸 禹,趙繼波,羅斌強,谷卓偉,王桂吉,張旭平,陳學秒,周中玉,李 牧,袁 紅,張紅平,王剛?cè)A,孫奇志,文尚剛,譚福利,趙劍衡,莫建軍,蔡進濤,金云聲,賀 佳,種 濤,趙小明,劉倉理

    (1. 中國工程物理研究院流體物理研究所,四川 綿陽 621999;2. 上海激光等離子體研究所,上海 201800;3. 深圳技術(shù)大學先進材料測試技術(shù)研究中心,廣東 深圳 518118;4. 深圳技術(shù)大學大數(shù)據(jù)與互聯(lián)網(wǎng)學院,廣東 深圳 518118;5. 中國工程物理研究院化工材料研究所,四川 綿陽 621999;6. 中國工程物理研究院應(yīng)用電子學研究所,四川 綿陽 621999;7. 中國工程物理研究院,四川 綿陽 621999)

    極端物理學顧名思義是極端條件下的物理學,但按當今科學前沿的實際理解是高能量密度狀態(tài)下的物理學,即高能量密度物理[1]。高能量密度狀態(tài)的公認條件是受載過程中物質(zhì)內(nèi)能的增量大于0.1 MJ/cm3,即其熱力學(靜水)壓力超過100 GPa[2],或者波長約1 μm 的輻射能束對物質(zhì)的輻照度大于3 PW/cm2。數(shù)十年前高能量密度狀態(tài)顯然只能在天體、行星、地球數(shù)千千米深部等環(huán)境以及核武器物理過程中存在,但產(chǎn)生高功率短脈沖激光的啁啾技術(shù)和美國Sandia 實驗室100 TW 級功率低阻抗強電流驅(qū)動器Z 機器的問世,使這種狀態(tài)已可在不少實驗室內(nèi)實現(xiàn)。研究高能量密度狀態(tài)下物質(zhì)、輻射及其相互作用的多學科交叉的極端物理學前沿領(lǐng)域,對慣性聚變、天體物理、行星物理和核武器物理研究具有重要現(xiàn)實意義。

    電磁驅(qū)動下宏觀物體的高能量密度動力學實驗,包括極端材料動力學、流體動力學行為和高溫高密度等離子體產(chǎn)生等,均具有重要的技術(shù)背景和應(yīng)用價值[3-4]。從驅(qū)動原理上看,這種實驗可分為2 類:磁驅(qū)動和磁壓縮。磁驅(qū)動系利用導體構(gòu)形表面上強驅(qū)動電流脈沖與自生磁場作用的電磁力加載,實現(xiàn)材料準等熵(斜波)壓縮或發(fā)射超高速固體飛片(統(tǒng)稱isentropic compression experiment, ICE),可用于沖擊壓縮、驅(qū)動金屬套筒內(nèi)爆、形成高溫高壓等離子體等。磁壓縮系利用炸藥爆轟或電磁驅(qū)動金屬套筒的低速內(nèi)爆運動,把較大體積中的磁通量壓縮到小區(qū)域中,形成很高磁場和磁壓,把該小區(qū)域內(nèi)的材料(尤其是體積較大的低密度材料)準等熵壓縮到高密度。利用炸藥爆轟為動力的磁通量壓縮裝置,又稱為磁通量壓縮(MC-1)發(fā)生器。這2 類實驗均涉及極端條件下實驗構(gòu)形的力學與物理學耦合問題,難以從理論方程組上完善考慮,但能通過數(shù)值模擬途徑用計算解決。

    1980 年代,美國Los Alamos 國家實驗室研制了一維輻射磁流體力學編碼RAVEN[5],具有一維任意拉格朗日-歐拉(arbitrary Lagrangian-Eulerian, ALE)功能,可以處理帶真空磁腔的負載構(gòu)形。RAVEN 編碼利用算子分裂方式,解決了流體彈塑性材料、雙流雙溫等離子體和高溫輻射物質(zhì)等強電磁加載下多物理場耦合模擬問題。目前,公認的最先進磁流體力學(magnetohydrodynamics, MHD)構(gòu)形計算編碼是美國Sandia 國家實驗室的多介質(zhì)多物理場ALE 編碼ALEGRA 系列。1990 年代初,Sandia 實驗室將高精度沖擊動力學有限差分(歐拉)編碼CTH 與改進的拉格朗日有限元編碼結(jié)合,得到在無結(jié)構(gòu)單元有限元網(wǎng)格上表述的ALE 編碼,即ALEGRA 最初的版本。經(jīng)過30 多年的開發(fā),ALEGRA 已成為功能強大的計算輻射磁流體-沖擊動力學-多物理場耦合的ALE 有限元編碼,得到了廣泛應(yīng)用[6]。

    本文中敘述的一維沖擊爆轟和磁流體力學多物理場計算編碼SSS-MHD,可對電磁驅(qū)動高能量密度動力學實驗開展磁流體力學多物理場數(shù)值模擬,分析強電流或者炸藥驅(qū)動的瞬態(tài)電磁場加載、受載實驗構(gòu)形和材料樣品中發(fā)生的彈塑性磁流體力學運動和熱力學變化等過程。SSS-MHD 編碼的基礎(chǔ)是1980 年代流體物理研究所研發(fā)的沖擊爆轟動力學編碼SSS[7],后又擴展了激光輻照效應(yīng)[8-9]、等離子體電磁驅(qū)動的計算功能[10-11]。隨著磁驅(qū)動和磁壓縮實驗的開展[12-13],SSS 編碼也向MHD 擴展,并于2012 年提出SSS-MHD 編碼的初稿[14-16],10 年來經(jīng)過許多實例計算核對和修改完善后基本定稿[17-22]。圖1 所示的各類實驗裝置基本概括了具有重要技術(shù)應(yīng)用背景的電磁驅(qū)動實驗類型,均可用SSS-MHD 編碼進行不同程度的模擬。

    圖1 適合SSS-MHD 編碼模擬的各種類型電磁驅(qū)動高能量密度物理實驗Fig. 1 Kinds of the magnetically-driven high-energy-density physics experiments suitably simulated with the SSS-MHD code

    多物理場耦合的數(shù)值研究,主要體現(xiàn)在多個物理系統(tǒng)方程組在離散化計算中進行恰當?shù)模〝?shù)學上完備的)相互耦合計算。SSS-MHD 編碼主要理論框架是彈塑性磁流體力學方程組,自身就是力學與電磁學兩大系統(tǒng)的耦合,進一步還能擴展到與等離子體物理和輻射流體力學等學科的耦合。由于使用場景和目的不同,磁流體力學方程組形式多樣。圖1 的各類實驗涉及低頻脈沖強電流和強磁場與導電體的短時間非相對論性相互作用,在電磁力和能量表達式中磁場的份額遠高于電場。這種意義下的磁流體力學近似,只需要2 個電磁學變量,如電流密度和磁感應(yīng)強度,相關(guān)的方程式為磁擴散方程和歐姆定律。上述耦合體現(xiàn)于運動方程中的電磁力、能量方程中的焦耳熱,同時磁擴散方程與介質(zhì)運動速度有關(guān),這種架構(gòu)為SSS-MHD 編碼結(jié)構(gòu)的簡化創(chuàng)造了有利條件。

    本文組成如下。第1 節(jié),敘述彈塑性磁流體力學多物理場耦合基本方程組;第2 節(jié),建立作為SSSMHD 編碼基礎(chǔ)的拉格朗日一維方程組的具體形式;第3 節(jié),敘述SSS-MHD 編碼的創(chuàng)新結(jié)構(gòu)和功能;第4 節(jié),展示平面構(gòu)形ICE 實驗的SSS-MHD 模擬,包括鋁、鉭、PBX 炸藥的準等熵斜波壓縮和鋁飛片高速發(fā)射實驗;第5 節(jié),展示圓柱面構(gòu)形的固體套筒電磁內(nèi)爆實驗及炸藥內(nèi)爆磁通量壓縮(MC-1)發(fā)生器實驗的模擬計算;第4~5 節(jié)將SSS-MHD 編碼計算與Sandia 實驗室Z 裝置的3 個高水平實驗和計算實例進行比較,而且模擬流體物理研究所CQ、CJ 系列裝置的4 個典型實驗;第6 節(jié)的結(jié)束語中指出,上述算例表明SSS-MHD 模擬與實驗數(shù)據(jù)以及ALEGRA 計算數(shù)據(jù)的相對偏差基本不超過5%。本文中的數(shù)值模擬還以豐富的物理量剖面圖顯示實驗過程全面的物理圖像。編碼SSS-MHD 的創(chuàng)新特色以及原有的沖擊動力學編碼的牢固基礎(chǔ),使其模擬通用性強、功能靈活全面、計算速度快、精度高。編碼SSS-MHD 的開發(fā)不但能為極端物理學、極端材料動力學實驗提供有力的數(shù)值模擬平臺,而且有助于多維MHD 多物理場編碼的研發(fā)。

    1 彈塑性磁流體力學基本方程組

    連續(xù)介質(zhì)力學運動的基本規(guī)律為質(zhì)量、動量和能量守恒定律,從一般的連續(xù)介質(zhì)動力學方程組可以方便地導出磁流體力學的基本方程組[26]。首先,在動量方程和能量方程中分別增加與電磁力和電阻產(chǎn)生熱量(焦耳熱)相關(guān)的項,把向量形式的連續(xù)介質(zhì)動力學方程組改寫如下,依次為連續(xù)性方程(質(zhì)量守恒)、運動方程(動量守恒)和能量方程(能量守恒):

    將式(2)全部點乘以u,再與式(3)相加,得到以總比能E=e+u·u/2 表達的能量方程:

    式(4)是SSS-MHD 編碼實際使用的能量方程,其優(yōu)點是簡化了張量運算。后文中,能量方程中無外體力F項,外源項Q寫成WL即單位質(zhì)量介質(zhì)吸收的激光輻射功率。

    在沖擊動力學或爆炸力學中,對連續(xù)介質(zhì)通常采用流體彈塑性的力學模型加以描述。將應(yīng)力張量Σ分解為 Σ=-pI+T,寫成分量式:

    式中:比容v=1/ρ 。SSS-MHD 編碼計算中材料EOS 基本均采用完全物態(tài)方程數(shù)據(jù)庫SESAME[27]。SESAME 庫中包括熱力學和低溫等離子體性質(zhì),為SSS-MHD 編碼多物理計算提供了重要基礎(chǔ)。由式(5)~(6)和未寫出的本構(gòu)關(guān)系,溫度以及應(yīng)力張量各個分量都可以通過微元的比內(nèi)能、比容、粒子速度分量及其導數(shù)或積分表達,此時未知的力學函數(shù)尚有e、 ρ 和u,式(1)~(3)恰好是2 個標量和1 個矢量方程,只要其中的電磁量J和B能借助于電動力學方程組耦合求解,即可實現(xiàn)方程組數(shù)學上的封閉。

    考慮電動力學之始,先假定介質(zhì)中不存在自由電荷,無需考慮位移電流和電介質(zhì)的極化或磁化,并得到電場強度E和磁感應(yīng)強度B都是無源的: ?·E=0 , ?·B=0 。麥克斯韋方程組還有如下定律。

    法拉第定律:

    安培定律:

    再增加廣義歐姆定律:

    前面方程組已含有粒子速度u為未知函數(shù),式(7)~(9)提供了求解E、J和B等3 個電磁量的方程,其中 μ 為介質(zhì)的磁導率。由于安培定律,式(1)~(3)中僅出現(xiàn)磁場B,在該意義下系統(tǒng)被稱為磁流體力學方程組。進一步化簡,可從式(7)~(9)中消去E和J,得到?jīng)Q定B的磁擴散方程:

    至此,式(1)~(6)、(8)、(10)以及未寫明的介質(zhì)材料本構(gòu)關(guān)系,組成了一般情況下多維彈塑性磁流體力學的封閉方程組,而且是用矢量、張量和微分算子表述的一般形式,提供了寫出任意正交曲線坐標系中分量形式微分方程組的方便途徑。

    2 一維拉格朗日形式彈塑性磁流體力學方程組

    根據(jù)第1 節(jié)中的基本方程組,可寫出一維情形下平面、柱面、球面坐標系中歐拉形式的彈塑性磁流體力學方程組(分量式),如下。連續(xù)性方程:

    運動方程:

    能量方程:

    在一維磁流體力學情形中,幾何指數(shù)N只能取0 或1,SSS-MHD 編碼的電磁學方程組(式(9)~(10))的分量形式如下。安培定律:

    磁擴散方程組:

    式中:Cθ=rNBθ稱為磁感應(yīng)強度的廣義 θ 分量,在柱對稱情形中有重要意義。

    在歐拉-拉格朗日形式的轉(zhuǎn)換中,時間t保持不變,拉格朗日空間坐標采用質(zhì)量坐標M,即計算構(gòu)形自最左邊界點到格點處的累計質(zhì)量。歐拉空間坐標r成為格點拉格朗日坐標M的位置函數(shù),位于r處格點的質(zhì)量 ?M則是厚度為dr、底面積分別是單位面積(平面)、r× 單位寬度(柱面)和r2(球面)的(單位弧度和立體角)體積中包含的介質(zhì)質(zhì)量。因此,2 個坐標系之間的轉(zhuǎn)換關(guān)系是:

    為了加以區(qū)分,將拉格朗日形式方程組中的r改寫為R,粒子速度u改寫為U,徑(縱)向應(yīng)力 σrr簡寫為 σ 。每個拉格朗日格點的質(zhì)量保持不變體現(xiàn)了質(zhì)量守恒,也省去了連續(xù)性方程。為此,需要添加一個速度方程,即U等于R的時間偏導數(shù)。

    由于高速度、高壓力問題參數(shù)范圍的特點,沖擊爆轟編碼SSS 采用了“沖擊波單位制”,其基本單位是質(zhì)量g、長度cm 和時間μs,因而速度單位是cm/μs 即10 km/s,壓力單位是g/(cm · μs2)=100 GPa。前面推導的電磁學方程組采用國際單位制(SI),因此在耦合方程組中電磁量相關(guān)項需添加修正系數(shù),達到與力學量單位的自洽。磁導率 μ=μ′μ0,其中真空磁導率 μ0=4π×10-7H/m , μ′為介質(zhì)的相對磁導率。

    經(jīng)過上述推導和轉(zhuǎn)換,得到SSS-MHD 編碼使用的拉格朗日形式一維彈塑性磁流體力學方程組如下。

    連續(xù)性方程和速度方程:

    運動方程:

    能量方程:

    安培定律(Cθ≡rNBθ):

    磁擴散方程組:

    對式(18)~(23)做離散時,時間離散格式是簡單的前向差分,空間離散中格點序號自左向右排列,除了每一格的位置、速度、電流密度定義在格子中點以外,其他力學量和電磁量均定義在格子右點。因此,速度方程和安培定律是中心差分格式,其余方程都是前向差分格式。整個差分方程組的具體形式可參看文獻[5]。

    3 SSS-MHD 編碼的結(jié)構(gòu)和功能

    3.1 計算編碼的結(jié)構(gòu)

    SSS-MHD 編碼的結(jié)構(gòu)立足于SSS 程序原有的框架[7],保持彈塑性流體動力學計算的各個子程序模塊構(gòu)架,進行適應(yīng)MHD 計算的改造。首先,在導程序中增加MHD 建模的內(nèi)容,在總卡(GENERL)和區(qū)卡(COMPNT)中引入與MHD 計算相關(guān)的數(shù)組、變量和指數(shù)等,然后,在專門的MHD 建模子程序(MHDZON)里再對有關(guān)區(qū)卡賦以MHD 負載構(gòu)形的類型指數(shù)、構(gòu)形參數(shù)以及初始條件等,完成MHD 計算建模。接著,在相關(guān)子程內(nèi)部和計算循環(huán)節(jié)點處加入控制MHD 計算的語句,如在計算循環(huán)開始處設(shè)定邊界條件狀態(tài)及參數(shù)的同時,檢查和調(diào)整空腔和磁腔的狀態(tài)及參數(shù)。

    隨著力學計算流程的進行,在運動方程子程序(VELOC)、能量方程子程序(ENERGY)和熱輸運計算子程序(HETCON)中,使用上時刻電磁變量值分別做洛倫茲力和焦耳熱計算。力學計算完成之后,流程進入專設(shè)的電磁計算程序包(MHDBLK,后詳)。最后是數(shù)據(jù)更新(UPDATE)、后處理等,進入下一個時間循環(huán),直至計算時間t達到設(shè)定時刻ts。編碼的計算流程見圖2。

    圖2 SSS-MHD 編碼的計算流程Fig. 2 Flowchart of the SSS-MHD code

    原編碼SSS 是一個框架式的通用程序,具有多種物理場耦合計算的功能,可以包容或擴展各種物態(tài)方程和物性數(shù)據(jù)庫,加入有關(guān)局部性物理過程的計算(例如激光與物質(zhì)的熱和沖量耦合、平衡電離的Saha 方程、動態(tài)斷裂的NAG 模型和含能材料的燃燒爆炸反應(yīng)等),適合平、柱、球?qū)ΨQ幾何的動力學問題,設(shè)有多種邊界條件和初始條件可供選擇,能方便地實現(xiàn)沖擊或爆轟的加載模擬[7],這些特點已全部為SSS-MHD 編碼沿用。下面各小節(jié)將具體闡述SSS-MHD 計算編碼的關(guān)鍵內(nèi)容和特色。

    3.2 MHD 負載構(gòu)形的統(tǒng)一建模

    在磁流體力學計算模型的建模中,由于力學與電磁量具有不同方向的坐標分量,首先需要設(shè)定力學構(gòu)形和相應(yīng)的加載磁場類型,然后賦予各個MHD 計算區(qū)域及格點處的電磁學物性參數(shù)和相應(yīng)初始值,必要時在整體模型邊界處補充電磁學邊界條件。其中的實質(zhì)性工作是負載構(gòu)形及磁場的分類和表述。

    圖3 是SSS-MHD 編碼規(guī)定的實驗負載構(gòu)形種類,N是一維連續(xù)性方程式(11)中的幾何指數(shù)。這樣的一維幾何中介質(zhì)速度只有徑(縱)向分量,但電流和磁場可以有2 個方向,即軸向(z)和角向(θ)分量。磁擴散相應(yīng)地有2 個方向的方程式,還可以用2 個驅(qū)動電路分別提供z、θ 方向的加載電流。一維磁流體力學原則上可以計算z、θ 及其組合的螺旋磁力線相關(guān)問題。負載構(gòu)形是各種材料的多層平行平板或者同軸圓柱殼層組成的幾何結(jié)構(gòu),其間可以包括空腔(有磁場的空腔稱為磁腔)。建模依據(jù)主要的電磁驅(qū)動機制,對構(gòu)形加載或被其壓縮的磁場(磁力線)有軸向(z)或角向(θ)以及兩者組合(z、θ 螺旋)等3 種形式。螺旋磁場中軸向與角向分量各自成體系、并行不悖,雖然圖中未畫出這個類別。從計算體系來說,只有圖3 的4 種基本構(gòu)形,其中平面z和θ 等2 種構(gòu)形實際相同,構(gòu)形的名稱表示其加載電流方向。構(gòu)形的板、殼可以是多層組合,也沒有具體畫出。

    圖3 SSS-MHD 編碼中實驗負載構(gòu)形的種類(紅線表示回流導體)[29]Fig. 3 Types of experimental configurations in the SSS-MHD code(red curves standing for return conductors)[29]

    圖1 中各類實驗的構(gòu)形和磁場均可抽象為圖3 的類型。負載構(gòu)形中的磁場可以是加載電流感生的(持續(xù)勵磁),也可以是初始時由外界給定的(瞬態(tài)種子磁場或定常勵磁);推動構(gòu)形運動的可以是加載電流自身形成的洛倫茲力(磁驅(qū)動),也可以利用炸藥爆轟或沖擊的力量(磁壓縮)。

    3.3 空腔和磁腔

    SSS-MHD 編碼沿用了SSS 編碼的空腔功能,并擴展到帶磁場的磁腔情形。計算模型通常為順序排列的介質(zhì)分區(qū)所組成的單連通區(qū),如果隊列中存在真空間隔并不妨礙歐拉形式計算,但對于拉格朗日計算則發(fā)生“零質(zhì)量格點”的困難,通常需采取分區(qū)分階段計算合成或?qū)嵭幸痪SALE 策略,這樣可能使得編碼復雜、機時增多,計算精度也受影響。

    SSS 編碼采取物質(zhì)區(qū)與空腔區(qū)一體化拉格朗日計算的方法。圖4 表示把空腔(區(qū)I)處理為只有3 個零質(zhì)量格點(J至J+2)的特殊區(qū),不參加物質(zhì)區(qū)運算,但其左右兩岸點J和J+3 分別具有左右物質(zhì)區(qū)岸點位置坐標R的數(shù)據(jù),空腔間隙等于R(J+3)-R(J) 。若空腔間隙大于零,左右兩岸點J和J+3 處均取自由邊界條件,空腔中間兩點J+1 和J+2 分別是J和J+3 點的鬼點(影子點)。如果空腔間隙為零或負值,意味著空腔閉合甚至兩岸碰撞過沖,需要退到前一時刻,適當減小時間步長重算,使得重算的空腔間隙在規(guī)定誤差下為零,并使點J和點J+3 的R值保持相等,成為同一個右岸介質(zhì)內(nèi)點。如果空腔原來是閉合的,但表觀重合的左、右岸點相互作用的拉力已變?yōu)槌^一定閾值的斥力,則空腔將被自動拉開。對于空腔開、閉變動的暫態(tài)振蕩過程,程序給予致穩(wěn)處理。大量算例表明,SSS 編碼的空腔一體化計算功能比較成功,步驟簡便,精度較高。

    圖4 SSS 編碼中空腔區(qū)及鄰區(qū)中格點的設(shè)置Fig. 4 Meshes in the cavity and nearby in the code SSS

    磁流體力學計算中磁空腔較常見,其力學運動可沿用SSS 編碼的處理方法,但需增加電磁學計算,并保持兩岸點磁場及磁力對稱相等的處理。經(jīng)過改進,SSS-MHD 編碼實現(xiàn)了一維磁流體力學多連通區(qū)的拉格朗日一體化計算,具備了對于較復雜實驗及裝置進行數(shù)值模擬的能力。

    磁腔處理方法基于電磁學一維問題的特點并利用了SSS 編碼的基礎(chǔ)。如圖3 所示,運用積分形式的安培定律和法拉第定律可以證明,實驗構(gòu)形中磁腔界面處的磁場總是與相應(yīng)電流分量正交的,磁腔左右岸導體表面的電流只能是面電流??涨唬ɑ螂娊橘|(zhì))中磁場的生成,無外乎饋入加載電流或初始種子磁場以及構(gòu)形導體中產(chǎn)生感應(yīng)電流等幾類原因,都可由磁擴散方程組計算確定。由于不考慮靜電電荷和位移電流,電介質(zhì)的電磁學性質(zhì)等同于磁腔。如果最右端的回流導體質(zhì)量較大或可不加考慮,則可簡化為最右端的單純磁腔(力學上是自由邊界,不是空腔)。此時的計算模型稱為單邊驅(qū)動,加載電流線只要一條,還需設(shè)定右端磁腔條件來計算該腔的磁通量。圓柱面θ 構(gòu)形可能存在軸線處的中心磁腔,但只能是環(huán)形電流驅(qū)動的均勻軸向磁場,如MC-1 裝置的情形。

    3.4 磁腔磁場的計算

    圖3 表明在平面一維情形(N= 0)、以一對無限大平行電極板為界的平磁腔(即陽陰極間隙AK-gap)中,加載電流可以有2 個方向,分別產(chǎn)生均勻磁場Bθ或Bz,磁力線均為直線;圓柱面一維(N= 1)幾何具有θ 或z等2 種構(gòu)形,同軸的無限長雙圓柱筒的間隙或單柱筒內(nèi)部均為圓(環(huán))磁腔,可近似推廣到雙多邊形柱筒之間的環(huán)形磁腔(環(huán)腔)情形。因加載電流有θ 或z方向,環(huán)腔中磁力線則為軸向直線或同軸環(huán)線。按照圖3 畫出的磁力線圍線,由安培定律給出平直腔和圓(環(huán))腔的磁場與構(gòu)形加載電流的關(guān)系,得到平直磁腔和圓(環(huán))磁腔中的磁場分別為:

    式中:I為構(gòu)形的總加載電流,i為該電流經(jīng)過單位寬度表面的平均電流(線密度),w為腔寬度或環(huán)腔的周線長度,下標表示電流或磁場的方向。按照一維幾何的原意,平直腔是無限寬的,總電流I無限大,式(24)是用有限寬度平均電流近似無限寬的理想狀況,帶來的誤差修正下面說明。式(25)適用于圓腔,由于非圓形環(huán)腔也是封閉的,可以近似采用圓環(huán)腔來處理,從而環(huán)腔的直線段部分相當于理想化的平直腔,雖然其電流密度因需按周線長度平均計算而大為下降。

    柱筒構(gòu)形按式(25)的近似很接近實際的柱面二維情況,誤差不大。然而,理想平面一維的式(24)與實際構(gòu)形的三維電磁場以及實驗結(jié)果的差異就較顯著。式(24)第一個等號右邊的因子 (I/w) 即以磁腔實際寬度w計算的平均加載電流密度,相當于用電流密度 (I/w) 、寬度無限的理想磁腔做近似。而如圖3所示,實際的平面一維“單條片”構(gòu)形(磁驅(qū)動實驗的一對加載電極板)是“折疊”后供加載電流來回的狹長金屬片,其間的AK-gap 即是磁腔。從能量角度看,式(24)的近似顯然是對磁腔中電流密度、磁場及磁壓的一種高估。陸禹[30]利用商業(yè)電磁學軟件[31]進行了二維平面電磁場計算,以頻率0.35 MHz 的正弦電流加載于一對銅電極板(截面寬2 mm、厚1 mm),兩板間的磁腔截面寬度(w)為2 mm、間隙(g)為0.1 mm。此無限長“單條片構(gòu)形”電磁響應(yīng)的二維平面高頻計算結(jié)果表示,加載電流(幅值)密度分布明顯不均勻,集中于有限尺度磁腔截面的4 個角點,證實磁腔中部電流密度低于邊部。這種計算提供了對有限尺度平直磁腔中電流密度分布不均勻性的定量估計。定義K為等效電流因數(shù),即利用商業(yè)電磁學軟件模擬計算得到的平面磁腔中部幅值電流密度與全寬度平均幅值電流密度的比值。將式(24)(不計下標)修正為:

    式中:修正后的電流密度iK=KI/w。平面一維磁腔的等效電流因數(shù)K是一個小于1 的正數(shù),與磁腔幾何參數(shù)w、g以及加載電流參數(shù)有關(guān),如果w→∞,g→0 ,則K→1 。

    嚴格說來,上述等效電流因數(shù)只反映實驗開始時的情況,由于運動中構(gòu)形動態(tài)電感的變化需要進行難以做到的至少是二維的磁擴散計算,才能對電流分布作更適當?shù)目紤]。動態(tài)電感的增大降低了加載電流的幅度,也可看作為運動過程中因數(shù)K的繼續(xù)下降,然而實際上MHD 計算包括了這部分的相互作用(詳見3.6 節(jié)),使得K動態(tài)變化的影響大為減弱,絕大部分計算中只需沿用初始的K值,即可得到全程與實驗數(shù)據(jù)相符的良好結(jié)果。這就是SSS-MHD 編碼處理平面磁腔的方法。

    本節(jié)中論述的平面磁腔磁場近似計算方法對于SSS-MHD 編碼的磁擴散計算具有重要意義,因為磁腔中的磁場Bz或Cθ都是常值,根據(jù)式(24)~(26)可由構(gòu)形加載總電流的即時值直接確定。該總電流(時間函數(shù))可以作為輸入數(shù)據(jù)給定,也可從電路方程組計算獲得。因而,這些磁腔的磁場值可以用作為磁擴散方程計算中可靠的定標值,避免了隱式計算的麻煩和誤差。該問題在理論上有更嚴謹?shù)姆椒?,如Lemke 等[32]利用Poisson 方程計算了含有四邊形環(huán)形磁腔的多連通區(qū)平面二維磁流體力學問題,但并不能得到有普遍意義的結(jié)果。

    3.5 顯式計算格式

    SSS-MHD 編碼主要用來模擬MHD 實驗構(gòu)形的力學和電磁運動,進而研究樣品材料的物態(tài)方程和材料動力學問題。為了滿足數(shù)值計算要求的穩(wěn)定性條件(如SSS 編碼中的Courant 條件),提高計算精度,并可防止鄰近格點尺度差距過大引起的剛性問題,希望采取一體化的拉格朗日顯式運算。力學方程組可以從給定邊界條件的左端點開始做顯式計算。把磁擴散方程同樣離散為顯式有限差分格式,但其邊界條件不在左、右端點給定,需要通過迭代使得本時刻的整體磁場分布、磁通量、反電動勢(構(gòu)形的動態(tài)電感)與加載或勵磁電流值達到高度自洽,才能獲得較高的計算精度,具體做法見3.7 節(jié)的說明。

    避免計算不穩(wěn)定性和剛性問題的重要措施是編碼中介質(zhì)區(qū)可選擇多種方法進行分格。計算值的振蕩往往發(fā)生于幾何或物理因素造成鄰近格子尺度或物理量急劇變化的區(qū)域,如內(nèi)爆聚心階段的中心范圍、點爆炸和散心爆轟開始階段的核心區(qū)域、高速碰撞的接觸面附近、脈沖電流的集膚層、高功率激光的吸收層和燒蝕層等。如果把這些區(qū)域做小尺度分格,并平滑過渡到周圍正常尺度分格的鄰區(qū),則有可能避免振蕩,得到物理上合理的計算結(jié)果。SSS-MHD 編碼中除了通常的等步長、等質(zhì)量分格方式外,還有在指定分區(qū)中使格子尺度(或質(zhì)量)按一定比例增大或縮小的分格方式。激光吸收層處理和從冷問題開始的等離子體演化等SSS 編碼計算,就是成功實例[8-9]。

    3.6 SSS-MHD 計算與驅(qū)動電路方程的耦合

    負載構(gòu)形的動態(tài)電感變化與驅(qū)動器放電過程之間存在重要的相互作用,若要理解這些復雜過程,首先應(yīng)把MHD 計算與驅(qū)動器電路計算進行耦合,建立準確預測實際加載電流的能力。與驅(qū)動電路的耦合是MHD 計算能否解決實際問題的一個關(guān)鍵,也是SSS-MHD 編碼的重要特色。這個要求至今仍是對強電磁驅(qū)動實驗計算模擬能力的重大挑戰(zhàn)。

    簡單實驗裝置中驅(qū)動電流的傳輸及成形通常由RLC 電路描述,其計算是常微分方程組的初值問題。復雜的大型裝置使用脈沖形成介質(zhì)線技術(shù),可由一階偏微分方程組描述,并可離散化為線性代數(shù)方程組求解。在與磁流體力學計算耦合方面,兩者沒有實質(zhì)性差別。下面以當前使用的集中參數(shù)RLC 電路為例,說明與MHD 計算的耦合的方法。先把RLC 電路寫為:

    式中:第1 個等號左邊為電容器組的剩余電壓,Q和C0分別為其剩余電荷量和初始電容值;U0和I分別為充電電壓和負載電流;L0、LL、R和VS依次為外電路(負載構(gòu)形)的固定電感、其他可變電感、固定電阻和開關(guān)電壓降(伏安特性),均在數(shù)據(jù)文件中輸入;ε 為負載構(gòu)形對于驅(qū)動電路的反電動勢,同一構(gòu)形可以有z、θ 等2 個驅(qū)動電路和2 套電路參數(shù)。

    驅(qū)動電路與負載構(gòu)形耦合的關(guān)鍵是計算反電動勢ε,即構(gòu)形路端電壓(電感性)VAK和沿其電流回路的電阻性電壓降VLP之和。根據(jù)積分形式的法拉第定律,VAK為構(gòu)形總磁通量 Φ 的時間導數(shù),即VAK=dΦ/dt;根據(jù)歐姆定律,有VLP=ηlLPI(或相應(yīng)的積分式),lLP為構(gòu)形電流回路長度。任何時刻,電容器組剩余電壓減去各項電壓降之后應(yīng)等于負載構(gòu)形的反電動勢。從物理上說,VAK就是構(gòu)形磁通量變化在圍繞其磁場的導體圍線中感應(yīng)的路端電壓。因此,電路耦合不但要求解電路方程,更需要求解磁擴散方程組得到構(gòu)形中的磁場分布,進行分區(qū)和全構(gòu)形的磁通量計算。

    3.7 SSS-MHD 計算程序包MHDBLK

    力學計算之后就可計算在上一時刻應(yīng)力、電磁力和各種能量作用下、在已發(fā)生力學變形并且即時負載電流作用下的構(gòu)形中電磁場的變化,也就是求解本時刻的電路方程和磁擴散方程組,這就是SSSMHD 編碼中體現(xiàn)MHD 多場耦合的程序包MHDBLK 的工作(圖5)。它的功能有3 項:(1)電路方程子程序(DLRK),依據(jù)上一時刻(j)已迭代n次的反電動勢ε 和負載電流I值解算電路方程(調(diào)用子程CURRNT 可設(shè)定該方程形式),得到負載電流n+1 次迭代值In+1;(2)磁擴散方程組子程序(MAGDFU),依據(jù)算出的In+1,求得磁場定標點處值,計算全部構(gòu)形范圍的磁場分布;(3)依據(jù)得到的磁場分布,計算負載構(gòu)形n+1 次迭代的總磁通量及反電動勢ε 的迭代值。若ε 迭代收斂不符合要求,則把n+1 次I和ε 的值輸入DLRK 子程序,進行下一次迭代計算,直到第n次和第n+1 次的反電動勢差值小于千分之一,才輸出本時刻j+1 的負載電流和磁場分布的確定值。

    圖5 SSS-MHD 編碼中計算程序包MHDBLK 的組成Fig. 5 Flow chart for the routine package MHDBLK in the SSS-MHD code

    若無法獲得驅(qū)動器電路方程的資料,加載電流也可選擇指定電流饋入方式,如實驗測量的實際負載電流或者另行估算的負載電流。此時,指定電流的數(shù)據(jù)以電流函數(shù)(CURFOM)形式饋入,直接把本時刻j+1 的指定電流值輸入MAGDFU,算得磁場后直接輸出,與電路無關(guān),不進行迭代。SSS-MHD 編碼整體為顯式運算,但通過同一時刻的負載電流與磁腔磁場的自洽需通過迭代計算,可得到與隱式運算相同的效果,而且精度更高。按照這種方法,把上時刻j到本時刻j+1 的迭代擴大到整個力學計算范圍也是可行的,但增加如此多的計算量對提高精度不一定有意義。

    4 SSS-MHD 編碼應(yīng)用于平面構(gòu)形磁驅(qū)動實驗的典型算例

    SSS-MHD 編碼針對各類實驗研究以及國外高水平實驗結(jié)果,進行了大量驗算和改進工作。本節(jié)列舉的平面構(gòu)形磁驅(qū)動實驗(ICE)算例分為2 類:材料樣品準等熵壓縮和高速金屬飛片發(fā)射。

    圖6 是上述2 類磁驅(qū)動實驗構(gòu)形的原理示意圖,圖6(a)~(b)都是單條片(strip-line)構(gòu)形;圖6(c)是方柱筒構(gòu)形,其芯柱(陰極)和側(cè)面板(陽電極)相當于4 個單條片,它們并聯(lián)構(gòu)成了四邊筒形立體結(jié)構(gòu)。樣品自由面或窗口界面速度歷史的測量采用激光干涉技術(shù),如VISAR(適用于任意反射面的速度干涉儀)、PDV(光子多普勒測速儀)等。當材料樣品受到很強壓縮脈沖時,其自由表面會發(fā)生微噴射現(xiàn)象,干擾激光干涉測速過程,這種情況必須采用力學阻抗適當?shù)墓鈱W單晶體片(LiF、藍寶石等)作為窗口,嚴密覆蓋樣品表面以抑制噴射,測量的數(shù)據(jù)則是樣品/窗口界面的速度歷史。

    圖6 磁驅(qū)動準等熵壓縮實驗的典型負載構(gòu)形Fig. 6 Typical loading configurations for magnetically-driven isentropic compression experiments (ICE)

    4.1 磁驅(qū)動平面材料樣品準等熵壓縮實驗(算例1~4)

    算例1 是對美國Sandia 實驗室Z 裝置實驗Z-1220[33]的驗算,該實驗采用類似圖6(c)的6061-T6 鋁合金4 層方筒構(gòu)形,將一個側(cè)面(陽電極板)加工為4 段樣品,厚度h分別為0.616、0.868、1.118 和1.371 mm,相當于一個4 階的臺階靶。每段樣品外表面粘有鍍增透膜的LiF 光學窗口。樣品與窗口的界面速度用VISAR 測量。實驗加載電流峰值約18 MA,上升沿約300 ns。美國Z-1220 實驗數(shù)據(jù)以及ALEGRA-1D 編碼與SSS-MHD 編碼的計算模擬都示于圖7。

    圖7 對于美國Sandia 實驗室磁驅(qū)動等熵壓縮實驗Z-1220[33]的SSS-MHD 模擬(算例1)Fig. 7 SSS-MHD simulations for theisentropic compression experiment Z-1220[33]at the Sandia National Laboratories, USA (example 1)

    SSS-MHD 計算是依據(jù)圖7(a)的加載電流實驗波形進行的,此例的矩形環(huán)腔周長5.8 cm,換算得知平均電流密度峰值約3.1 MA/cm,等效電流因數(shù)K為1。4 種厚度樣品的2 種數(shù)值模擬速度曲線與實驗數(shù)據(jù)的符合程度都較好。SSS-MHD 計算給出的材料樣品最高壓力為51 GPa,文獻[33]給出的最高縱向應(yīng)力為54 GPa。2 種模擬的差別在于SSS-MHD 計算在彈塑性轉(zhuǎn)變階段速度起跳過早,其原因可能是該計算的樣品材料是純鋁而不是Z 實驗的6061-T6 鋁合金,此外還需要考察計算中輸入電流波形初始50 ns 段落的細節(jié)以及鋁樣品中磁場擴散的情況,進而或可做出相應(yīng)改進。對于模擬速度曲線末期與實驗數(shù)據(jù)的偏差,主要是因為鋁的彈塑性行為與模擬選用的材料動態(tài)本構(gòu)模型之間存在差異。

    算例2~4 是對流體物理研究所CQ 系列脈沖功率裝置上相關(guān)實驗[14,16,29]的驗算,材料樣品包括重金屬鉭、輕金屬鋁和以HMX 為基的塑料黏結(jié)炸藥(PBX)JO-9 159,實驗參數(shù)及結(jié)果列于表1。

    表1 磁驅(qū)動準等熵壓縮實驗算例2~4 的主要參數(shù)[14,16]Table 1 Parameters of examples 2-4 for magnetically-driven isentropic compression experiments[14,16]

    算例2~4 加載電流和樣品表面速度的實驗和計算結(jié)果示于圖8。由于CQ 系列裝置早期的實驗樣品寬度較窄,算例2(Ta)和算例4(JO-9 159)的等效電流因數(shù)均取為0.6,算例3(Al)的等效電流因數(shù)取為0.4;后續(xù)電流較高的算例中,該因數(shù)大多為0.7 以上。SSS-MHD 編碼具有2 種加載電流輸入方式:電路(CKT)方程組耦合計算(見圖5)和實驗負載電流(CUR)數(shù)據(jù)饋入,圖8 將這2 種計算方式做了比較。早期工作都是使用CKT 計算的[14-16],隨著強電流測量技術(shù)的提高,本文中補充了CUR 計算,得到了更全面的認識。

    圖8 對于流體物理研究所CQ 裝置磁驅(qū)動準等熵壓縮實驗的SSS-MHD 模擬(算例2~4)Fig. 8 SSS-MHD simulations for magnetically-driven quasi-isentropic compression experiments at IFP (examples 2-4)

    將圖8 中3 種材料樣品臺階靶速度波形處理后,可獲得它們的準等熵壓縮線,其高端壓力如表1 右端所示。算例2 中鉭金屬樣品的密度高,同時該例的加載電流更強,使得鉭的等熵線高端達80 GPa 以上。根據(jù)圖8(a)可知,如果加載電流密度提高3 倍多,等熵壓縮壓力可提高一個量級,能得到重金屬1 TPa以下范圍內(nèi)的等熵線,正如美國Z 裝置上開展的實驗[34]。當然,要把磁驅(qū)動加載電流強度提高若干倍,對驅(qū)動器硬件和實驗技術(shù)都極為困難。

    表1 中提及了電流計算(CUR)與電路計算(CKT)方式的比較。SSS-MHD 編碼CKT 計算采用的集中參數(shù)RLC 電路可以描述很多實際驅(qū)動電路,但電路中若使用伏安特性復雜的氣體開關(guān),則RLC 電路方程就不夠準確。圖8(d)中CKT 計算優(yōu)于CUR 方式,原因可能是該算例驅(qū)動器CQ-1.5 裝置采用了伏安性能簡單的爆炸開關(guān)[35],正適合于CKT 計算的簡單RLC 電路描述(早期實驗使用的電流測量技術(shù)也存在較大誤差)。其余2 例的驅(qū)動器CQ-4 電路采用氣體開關(guān),雖然目前簡單的CKT 計算也可達到偏差低于10%的精度,但明顯不如CUR 方式。本文電路耦合CKT 計算的原理是正確的,但對具體電路的描述還有待改進。

    4.2 磁驅(qū)動高速金屬飛片實驗(算例5)

    磁驅(qū)動高速金屬飛片實驗利用等熵壓縮的磁力斜波,以極高壓力(數(shù)百吉帕以至太帕)驅(qū)動作為電極板一部分的金屬飛片高速射出,鋁、銅飛片的速度很容易達到10 km/s 以上[21],2011 年美國Sandia 實驗室已能把亞毫米厚度的鋁飛片發(fā)射到45 km/s[36],開創(chuàng)了沖擊壓縮實驗的新境界,對極端物理學研究具有重要意義。飛片運動前期由強磁壓驅(qū)動,電流衰減的后期則由飛片后側(cè)電弧燒蝕等離子體推進,這類電磁熱力強耦合的復雜實驗必須采取多物理場耦合計算進行設(shè)計或核算。

    算例5 是Lemke 等[36]2011 年所做系列實驗之“11 mm-2 s”的模擬,該系列實驗研究了美國Sandia實驗室Z 裝置驅(qū)動高速金屬飛片的能力,給定了用于直至TPa 范圍的沖擊壓縮標準加載手段。實驗“11 mm-2 s”是文獻[36]中唯一報道了加載電流波形的具體實驗。該實驗采用單條片構(gòu)形,材料為鋁,電極板寬度為11 mm、長度為25 mm,構(gòu)形的盲孔為矩形凹槽,槽底即矩形金屬飛片,AK-gap 即磁腔間隙為1 mm。電極板凹槽頂面直接覆蓋7.5 mm 厚的銅靶板,凹槽深度即飛片的飛行距離約為7.8 mm。實驗測量飛片后自由面速度波形,并提供其剖面參數(shù)及高速碰撞過程的信息。計算采用的等效電流因數(shù)K為文獻[36]給定值0.76(本文自算值為0.73)。圖9 為SSS-MHD 編碼的模擬結(jié)果及其比較。

    圖9 對于Sandia 實驗室Z 裝置鋁飛片實驗“11 mm-2 s”的SSS-MHD 模擬(算例5)Fig. 9 SSS-MHD simulations for the Al flyer experiment “11 mm-2 s”on the Z machine at the Sandia National Laboratories, USA (example 5)

    圖9(a)顯示SSS-MHD 計算的“11 mm-2 s”飛片速度與原實驗及ALEGRA-2D 編碼模擬相當一致,但在3.2 μs 后有些偏高。圖9(b)顯示在磁力推動后期(3.05 μs)磁壓力高達285 GPa。圖9(c)表示飛行后期3.24 μs 時飛片自由面后固體區(qū)溫度過高(1 000 K 左右),導致密度略低于常值,這個偏差與圖9(a)后期計算速度偏高的情況相符。雖然這里沒有進行飛片中熔化邊界的計算,但從圖9(c)中3 條剖面曲線斜率變化的拐點可以判斷,飛片后自由面固態(tài)區(qū)厚度約為0.1 mm,符合Z 裝置飛片實驗通常的數(shù)據(jù)。圖9(b)是與圖9(c)時間不同的剖面,但與文獻[36]中另一實驗“11 mm-1 s”剖面的數(shù)據(jù)相當一致,提供了有依據(jù)的物理圖像。

    5 SSS-MHD 編碼應(yīng)用于圓筒構(gòu)形電磁內(nèi)爆和MC-1 發(fā)生器實驗的典型算例

    本節(jié)計算的金屬套筒電磁內(nèi)爆實驗和CJ 型MC-1 發(fā)生器實驗,可使材料樣品中斜波壓力達到太帕量級,或可使構(gòu)形軸線處的壓縮磁場達到700~1 000 T(相當于磁壓200~400 GPa),對于低密度物質(zhì)的壓縮具有特殊意義。

    5.1 固體金屬套筒高速電磁內(nèi)爆實驗(算例6)

    磁驅(qū)動圓筒構(gòu)形的內(nèi)爆實驗可以高效實現(xiàn)樣品材料的斜波壓縮或套筒的高速度高對稱性內(nèi)爆運動,其負載構(gòu)形即是作為陰、陽電極的內(nèi)外兩層同軸金屬圓筒,其間為圓環(huán)磁腔。這種構(gòu)形的主體—作內(nèi)爆運動的內(nèi)筒稱為套筒(liner);外筒稱為回流筒(也可由若干金屬柱組成),做速度不快的外向擴展運動,對內(nèi)爆影響不大,然而它回流的電路功能不可或缺。由于套筒內(nèi)爆運動是會聚的,其電磁驅(qū)動力和內(nèi)聚壓力越來越強,電流有效作用時間和效率都明顯高于平面構(gòu)形。金屬套筒內(nèi)爆速度有可能超過100 km/s,內(nèi)爆壓力可達到數(shù)太帕,成為高能量密度物理宏觀樣品實驗的主要手段。

    算例6 是對于美國Z 機器磁驅(qū)動Al/Cu 復合套筒超高壓內(nèi)爆壓縮實驗[34]的SSS-MHD 數(shù)值模擬。圖10(a)是該實驗裝置的截面示意圖,圖中的陽極—回流筒(厚度0.45 mm、內(nèi)半徑13 mm)和陰極—內(nèi)套筒的驅(qū)動層(厚度1 mm、外半徑3.43 mm)由6061-T6 鋁合金制造。外筒的膨脹速度歷史由裝置外圍排列的VISAR 探針測量。內(nèi)套筒的內(nèi)層—樣品層為銅材料,厚度0.53 mm、內(nèi)半徑1.9 mm(此系列還進行了鉭和鋁樣品材料的內(nèi)爆壓縮實驗)。套筒的有效高度為10 mm,內(nèi)窺PDV 探針由6 根PDV 光纖組成,封裝于壁厚0.15 mm、外半徑0.35 mm 的鉑管之中。圖10(b)顯示SSS-MHD 編碼的計算速度與實驗及ALERGRA-1D 模擬結(jié)果相符。當內(nèi)爆中止即銅套筒撞上鉑管時(約3.016 μs),波形優(yōu)化設(shè)計的加載電流達到約16.5 MA,此時鋁驅(qū)動層外半徑約2 mm,電流線密度高達13 MA/cm,磁壓力超過1 TPa,并且壓縮過程是準等熵的。圖10(c)顯示銅樣品中最高壓力為1.25 TPa,最高密度達22.4 g/cm3,可看出SSS-MHD編碼計算的剖面臺階與ALEGRA-2D 模擬結(jié)果略有差別,主要原因是介質(zhì)密度接觸間斷的處理存在缺陷,有可能是本文計算中空間格子尺度偏大等原因造成的,需要改進。

    圖10 對于Sandia 實驗室Z 裝置上金屬套筒電磁內(nèi)爆實驗的SSS-MHD 模擬(算例6)Fig. 10 SSS-MHD simulation of magnetically-driven liner implosions on the Z-machine at the Sandia National Laboratories, USA (example 6)

    電磁內(nèi)爆壓縮實驗的意義在于獲得太帕等級的物態(tài)方程數(shù)據(jù),關(guān)鍵在于如何較準確地處理實驗數(shù)據(jù)。文獻[34]中經(jīng)過精確的物態(tài)方程計算,得到圖10(c)中銅樣品的最高壓力、密度標定值的偏差均小于1%,是否意味著已具有狀態(tài)方程參考價值,這些重要問題需要繼續(xù)從實驗和模擬方面進行探討。

    5.2 圓柱形炸藥內(nèi)爆磁通量壓縮實驗(算例7)

    MC-1 發(fā)生器利用圓柱形炸藥驅(qū)動金屬套筒內(nèi)爆,把套筒空腔內(nèi)種子磁場的初始磁通量壓縮到軸線周圍小區(qū)域中,從而得到高磁場和高磁壓力[37]。MC-1 發(fā)生器與5.1 節(jié)套筒電磁內(nèi)爆的基本不同點在于MC-1 套筒是圖3 中的圓柱形(N=1)θ 構(gòu)形,而電磁內(nèi)爆套筒屬于z構(gòu)形。MC-1 套筒的環(huán)形“勵磁”電流物理上是注入套筒腔內(nèi)的種子磁場Bz0在內(nèi)爆套筒中感生的渦電流。從圓柱電磁構(gòu)形的分析可知,z構(gòu)形產(chǎn)生的環(huán)形磁場驅(qū)動套筒內(nèi)爆具有較高效率,θ 構(gòu)形壓縮腔內(nèi)軸向磁通量則有較高效率。

    CJ-100 型MC-1 裝置是中國工程物理研究院流體物理所研制的小型MC-1 發(fā)生器[38-39],圖11 是該裝置的截面示意圖和實物照片。CJ-100 型裝置的套筒材料為304 不銹鋼,外直徑100 mm、厚度1.5 mm、高度210 mm。炸藥為RHT-901(m(TNT)/m(RDX)=40/60),尺寸為內(nèi)直徑100 mm、外直徑210 mm、高度65 mm,總炸藥量不到3 kg。

    圖11 CJ-100 型MC-1 發(fā)生器的示意圖和實物照片F(xiàn)ig. 11 Schematics and picture of the CJ-100 type MC-1 generator

    算例7 是CJ-100 型MC-1 裝置磁場壓縮實驗Shot20150630[39]的SSS-MHD 驗算。該實驗套筒內(nèi)部的種子磁場Bz0約為5.5 T。該實驗中軸線附近磁場增長歷史以及套筒內(nèi)、外半徑變化的模擬計算示于圖12(a),在磁探針停止工作之前實測磁場與計算值十分符合。

    圖12 CJ-100 型MC-1 發(fā)生器磁通量壓縮實驗 及其性能的SSS-MHD 模擬計算(算例7)Fig. 12 SSS-MHD simulations for the magnetic flux compression experiment and the performances of the CJ-100 type MC-1 generator (example 7)

    計算的套筒內(nèi)層回彈半徑約2.6 mm,此時計算的峰值磁場約為1 450 T。但實測磁場峰值為690 T,從圖12(a)看出此時套筒腔內(nèi)磁場區(qū)的半徑值約3.8 mm,此值相當于實驗中磁場測點的半徑值。圖12(b)是不同種子磁場Bz0之下SSS-MHD 計算的CJ-100 型發(fā)生器達到的峰值磁場和套筒回彈半徑,可作為實驗設(shè)計的參考。圖12(b)中還畫出了該型號發(fā)生器以前2 發(fā)實驗EX1、EX2 的壓縮磁場,以及對相應(yīng)磁場測點半徑值的估計,數(shù)據(jù)與本例相近。從圖12(b)可知,設(shè)計適當?shù)姆N子磁場值,可得到折中的壓縮磁場值和可利用磁場的空間范圍。

    6 結(jié) 束 語

    SSS-MHD 編碼是拉格朗日形式一維平面和圓柱面幾何、多介質(zhì)、多組分、多連通區(qū)的彈塑性磁流體力學多物理場計算編碼,是一維沖擊爆轟編碼SSS 向磁流體力學擴展而形成,具有25 個子程序和函數(shù),共約6 000 條Fortran 語句。SSS-MHD 編碼具有框架式結(jié)構(gòu),便于靈活增加和擴充介質(zhì)類型、物性數(shù)據(jù)及物理功能。目前,除了SSS 編碼原有的功能外,SSS-MHD 編碼主要面向強電磁驅(qū)動的高能量密度動力學實驗,特別是為極端材料動力學實驗提供數(shù)值模擬平臺。

    本文中,第1~3 節(jié)闡述了SSS-MHD 編碼的理論基礎(chǔ)、程序結(jié)構(gòu)和特色,第4~5 節(jié)給出了該編碼模擬強電磁驅(qū)動實驗的7 個算例,包括平面構(gòu)形材料磁驅(qū)動準等熵壓縮實驗、磁驅(qū)動發(fā)射高速金屬飛片實驗、套筒電磁內(nèi)爆材料壓縮實驗和MC-1 發(fā)生器磁通量壓縮實驗??傮w看來,SSS-MHD 模擬與實驗及其他方法數(shù)值模擬結(jié)果符合較好,相對偏差均在5%以下。這些算例的實驗形式廣泛多樣,物理變量范圍寬廣,實驗測量和模擬計算精度較高,為考核SSS-MHD 編碼的實際能力提供了合適的試題,同時說明該編碼已達到強電磁驅(qū)動實驗一維數(shù)值模擬平臺的要求,可以推廣應(yīng)用。

    通過與該領(lǐng)域先進實驗及數(shù)值模擬結(jié)果的核對,為SSS-MHD 編碼的繼續(xù)改進和擴展指出了方向。例如,補充流體力學模型,適應(yīng)內(nèi)爆動力學實驗模擬的需要;在基本方程組層面引入更普遍的非線性連續(xù)介質(zhì)力學本構(gòu)理論,以適應(yīng)大變形高應(yīng)變率計算的需要;改進單溫等離子體為多溫多流模型,向簡單輻射磁流體力學計算發(fā)展;建立更精細準確的驅(qū)動器電路和關(guān)鍵器件模型,實現(xiàn)更準確的MHD 與驅(qū)動電路的耦合計算;改進高功率激光與物質(zhì)耦合計算模型,建立與SOP (streaked opticapyrometer)等先進測溫技術(shù)配合的功能,進一步擴大對高能量密度物理實驗的模擬能力。

    本文工作得到了流體物理研究所諸多磁驅(qū)動、磁壓縮實驗數(shù)據(jù)的支持,謹向參加有關(guān)實驗的仝延錦、唐小松、宋振飛、程誠、李建明、匡學武、吳剛、稅榮杰、胥超、鄧順益、馬驍?shù)韧局乱灾孕母兄x!

    猜你喜歡
    構(gòu)形算例方程組
    深入學習“二元一次方程組”
    雙星跟飛立體成像的構(gòu)形保持控制
    《二元一次方程組》鞏固練習
    通有構(gòu)形的特征多項式
    一類次臨界Bose-Einstein凝聚型方程組的漸近收斂行為和相位分離
    對一個幾何構(gòu)形的探究
    基于振蕩能量的低頻振蕩分析與振蕩源定位(二)振蕩源定位方法與算例
    互補問題算例分析
    基于CYMDIST的配電網(wǎng)運行優(yōu)化技術(shù)及算例分析
    非自治耗散Schr?dinger-Boussinesq方程組緊致核截面的存在性
    97热精品久久久久久| 99久久无色码亚洲精品果冻| 超碰av人人做人人爽久久| 亚洲五月天丁香| 91精品伊人久久大香线蕉| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 国产三级中文精品| 看十八女毛片水多多多| 成人欧美大片| 国产成人a区在线观看| 国产成人a区在线观看| 黄色欧美视频在线观看| 欧美日本视频| 亚洲经典国产精华液单| 国产亚洲av嫩草精品影院| 国产精品国产三级国产av玫瑰| 日日撸夜夜添| 中文精品一卡2卡3卡4更新| 听说在线观看完整版免费高清| 一区二区三区免费毛片| 亚洲精品国产av成人精品| 最近2019中文字幕mv第一页| 欧美高清性xxxxhd video| 国产又黄又爽又无遮挡在线| 99热这里只有是精品50| 欧美变态另类bdsm刘玥| 国产伦在线观看视频一区| 国语自产精品视频在线第100页| 人人妻人人看人人澡| 黄色欧美视频在线观看| 网址你懂的国产日韩在线| 色网站视频免费| 欧美区成人在线视频| 在线a可以看的网站| 亚洲国产精品合色在线| av天堂中文字幕网| 国产成人一区二区在线| 最近中文字幕2019免费版| 22中文网久久字幕| 亚洲欧美日韩无卡精品| 26uuu在线亚洲综合色| 久久婷婷人人爽人人干人人爱| 久久久久久久久久黄片| 男人舔女人下体高潮全视频| 精品久久久久久久久亚洲| 神马国产精品三级电影在线观看| 国产精品国产高清国产av| 天天一区二区日本电影三级| 九九在线视频观看精品| 久久亚洲国产成人精品v| 亚洲丝袜综合中文字幕| 亚洲在线观看片| 看片在线看免费视频| 嘟嘟电影网在线观看| 亚洲va在线va天堂va国产| 国产高清有码在线观看视频| 黑人高潮一二区| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 国产精品麻豆人妻色哟哟久久 | 亚洲aⅴ乱码一区二区在线播放| 亚洲高清免费不卡视频| 麻豆精品久久久久久蜜桃| 高清在线视频一区二区三区 | 91精品一卡2卡3卡4卡| 国产成人a区在线观看| 国产v大片淫在线免费观看| 中文乱码字字幕精品一区二区三区 | 91精品一卡2卡3卡4卡| 天天一区二区日本电影三级| 国产av码专区亚洲av| 免费av观看视频| 国产午夜精品一二区理论片| 亚洲,欧美,日韩| 亚洲无线观看免费| 国产伦一二天堂av在线观看| 欧美最新免费一区二区三区| 22中文网久久字幕| 男人舔女人下体高潮全视频| 国模一区二区三区四区视频| 亚洲美女搞黄在线观看| 欧美日韩综合久久久久久| 热99在线观看视频| 91精品伊人久久大香线蕉| 蜜桃久久精品国产亚洲av| 亚洲性久久影院| 汤姆久久久久久久影院中文字幕 | 日本av手机在线免费观看| 国产精品一二三区在线看| 国产成人午夜福利电影在线观看| 日日啪夜夜撸| 校园人妻丝袜中文字幕| 国产91av在线免费观看| 久久鲁丝午夜福利片| 婷婷色综合大香蕉| 亚洲高清免费不卡视频| 中文在线观看免费www的网站| 成人亚洲欧美一区二区av| 婷婷六月久久综合丁香| 亚洲高清免费不卡视频| 久久精品国产鲁丝片午夜精品| 成年版毛片免费区| 在线观看一区二区三区| av卡一久久| 国产黄片视频在线免费观看| 最近手机中文字幕大全| 天堂影院成人在线观看| 寂寞人妻少妇视频99o| 99久久精品热视频| 亚洲最大成人手机在线| 国产欧美日韩精品一区二区| 亚洲精品乱码久久久久久按摩| 亚洲熟妇中文字幕五十中出| 国产淫片久久久久久久久| 99久久成人亚洲精品观看| 免费av不卡在线播放| 国产不卡一卡二| 卡戴珊不雅视频在线播放| 精品久久久久久成人av| 欧美区成人在线视频| 菩萨蛮人人尽说江南好唐韦庄 | 亚洲成av人片在线播放无| 亚洲精品自拍成人| 最近中文字幕高清免费大全6| 久久草成人影院| 国产免费男女视频| 男女下面进入的视频免费午夜| 精品免费久久久久久久清纯| 亚洲人成网站在线播| 国产精品一区二区在线观看99 | 床上黄色一级片| 简卡轻食公司| 久久久久免费精品人妻一区二区| 久久精品久久久久久噜噜老黄 | 亚洲精品乱码久久久久久按摩| 久久久精品94久久精品| 国产视频内射| 国产精品1区2区在线观看.| 色综合亚洲欧美另类图片| 99久国产av精品国产电影| 久久久精品大字幕| 中文字幕av成人在线电影| 国产一区二区亚洲精品在线观看| 人妻系列 视频| 国产伦一二天堂av在线观看| 成人国产麻豆网| 亚洲国产最新在线播放| 久久久久久久亚洲中文字幕| 男的添女的下面高潮视频| 亚洲国产欧洲综合997久久,| 国产片特级美女逼逼视频| 亚洲av中文字字幕乱码综合| 久久国产乱子免费精品| 一区二区三区乱码不卡18| 国产又黄又爽又无遮挡在线| 日韩欧美 国产精品| 欧美性猛交黑人性爽| 日本黄色片子视频| av.在线天堂| 亚洲熟妇中文字幕五十中出| 久久精品久久久久久久性| 亚洲国产精品合色在线| 久久精品影院6| 在线免费观看不下载黄p国产| 高清视频免费观看一区二区 | 午夜亚洲福利在线播放| 如何舔出高潮| 99视频精品全部免费 在线| 亚洲精品色激情综合| 麻豆成人午夜福利视频| 亚洲国产日韩欧美精品在线观看| 国产 一区精品| 欧美一区二区亚洲| 九九爱精品视频在线观看| 亚洲成人av在线免费| 中文精品一卡2卡3卡4更新| 久99久视频精品免费| 免费av观看视频| av在线观看视频网站免费| 视频中文字幕在线观看| 中文字幕人妻熟人妻熟丝袜美| 99视频精品全部免费 在线| 亚洲精品色激情综合| 午夜爱爱视频在线播放| 日日摸夜夜添夜夜添av毛片| 亚洲欧美日韩无卡精品| 国内精品宾馆在线| 日韩欧美在线乱码| 精品一区二区免费观看| 最近的中文字幕免费完整| 国产又色又爽无遮挡免| 久久99蜜桃精品久久| 最近最新中文字幕免费大全7| 精品国产一区二区三区久久久樱花 | 亚洲欧美日韩卡通动漫| 亚洲国产最新在线播放| 男人狂女人下面高潮的视频| 国产伦精品一区二区三区四那| 国产老妇女一区| 国产亚洲5aaaaa淫片| 亚洲av成人精品一二三区| 黑人高潮一二区| 中国国产av一级| 国产人妻一区二区三区在| 又爽又黄a免费视频| 免费看光身美女| 亚洲av中文av极速乱| 村上凉子中文字幕在线| 免费播放大片免费观看视频在线观看 | 亚洲国产精品sss在线观看| 九九热线精品视视频播放| 国产黄色视频一区二区在线观看 | 五月伊人婷婷丁香| 狠狠狠狠99中文字幕| 老司机福利观看| 成人欧美大片| 久久这里有精品视频免费| 亚洲无线观看免费| 精品人妻偷拍中文字幕| 国产伦精品一区二区三区视频9| 亚洲一级一片aⅴ在线观看| 狠狠狠狠99中文字幕| a级毛色黄片| 成人欧美大片| 男的添女的下面高潮视频| 欧美精品国产亚洲| 久久国内精品自在自线图片| 午夜福利在线在线| 七月丁香在线播放| 亚洲av熟女| 久久久国产成人精品二区| 国产在线男女| 久久久久久九九精品二区国产| 婷婷色综合大香蕉| 国产亚洲午夜精品一区二区久久 | 一个人免费在线观看电影| 午夜久久久久精精品| 男女啪啪激烈高潮av片| 亚洲欧美日韩卡通动漫| 天堂av国产一区二区熟女人妻| 精品国产三级普通话版| 成人性生交大片免费视频hd| 夫妻性生交免费视频一级片| 久久99热这里只频精品6学生 | 久久久久久伊人网av| 国产乱来视频区| 在线观看66精品国产| 国产一区二区在线av高清观看| 免费看av在线观看网站| 亚洲精品aⅴ在线观看| av在线播放精品| 亚洲精品乱码久久久v下载方式| 日本免费一区二区三区高清不卡| 嫩草影院新地址| 久久久精品大字幕| 亚洲av成人精品一二三区| 国产乱来视频区| 亚洲电影在线观看av| 我要搜黄色片| 99久久精品一区二区三区| 日韩欧美精品v在线| 插阴视频在线观看视频| 精品少妇黑人巨大在线播放 | 91精品一卡2卡3卡4卡| 在线免费观看的www视频| 亚洲av免费高清在线观看| 日韩国内少妇激情av| 国产成人福利小说| 十八禁国产超污无遮挡网站| 亚洲美女视频黄频| 国产v大片淫在线免费观看| ponron亚洲| 天天躁日日操中文字幕| 日本熟妇午夜| 美女xxoo啪啪120秒动态图| 久久久精品大字幕| 国产成人免费观看mmmm| 亚洲av不卡在线观看| 久久精品人妻少妇| 99在线视频只有这里精品首页| 一个人看视频在线观看www免费| 变态另类丝袜制服| 国产亚洲av片在线观看秒播厂 | 免费播放大片免费观看视频在线观看 | 成人一区二区视频在线观看| 美女黄网站色视频| 女的被弄到高潮叫床怎么办| 夜夜看夜夜爽夜夜摸| 在线观看美女被高潮喷水网站| 亚洲av成人精品一区久久| 91av网一区二区| 久久久国产成人免费| .国产精品久久| 听说在线观看完整版免费高清| 国产 一区 欧美 日韩| 免费看a级黄色片| 大香蕉97超碰在线| 精品国产露脸久久av麻豆 | 国产在视频线精品| 老司机影院毛片| 免费看美女性在线毛片视频| 精品久久久久久久人妻蜜臀av| 菩萨蛮人人尽说江南好唐韦庄 | 亚洲真实伦在线观看| 1024手机看黄色片| 国产极品精品免费视频能看的| 少妇猛男粗大的猛烈进出视频 | 高清毛片免费看| 亚洲自拍偷在线| 麻豆成人午夜福利视频| 日本免费在线观看一区| 最后的刺客免费高清国语| 亚洲最大成人中文| 一区二区三区四区激情视频| 不卡视频在线观看欧美| 偷拍熟女少妇极品色| 国产探花在线观看一区二区| 亚洲熟妇中文字幕五十中出| 日韩视频在线欧美| 国产精品一区二区三区四区久久| 中文在线观看免费www的网站| 欧美bdsm另类| 国产精品伦人一区二区| 日本色播在线视频| 久久人人爽人人片av| 国产免费一级a男人的天堂| 午夜老司机福利剧场| 大香蕉久久网| 国产精品爽爽va在线观看网站| 国产色婷婷99| 亚洲熟妇中文字幕五十中出| 国产淫语在线视频| 美女国产视频在线观看| 精品久久久久久成人av| 少妇丰满av| 日本黄色片子视频| 久久99蜜桃精品久久| 欧美xxxx性猛交bbbb| 我的老师免费观看完整版| 一级毛片我不卡| 免费无遮挡裸体视频| av视频在线观看入口| 伦精品一区二区三区| 国产视频首页在线观看| 97超碰精品成人国产| 亚洲熟妇中文字幕五十中出| 国产色爽女视频免费观看| 99热网站在线观看| 99久久成人亚洲精品观看| 最新中文字幕久久久久| 九九爱精品视频在线观看| 亚洲自拍偷在线| 日本黄色片子视频| 深夜a级毛片| 桃色一区二区三区在线观看| 亚洲av中文av极速乱| 免费黄色在线免费观看| 寂寞人妻少妇视频99o| 一本久久精品| 日本黄色视频三级网站网址| av卡一久久| 色综合色国产| or卡值多少钱| 色5月婷婷丁香| 一级毛片电影观看 | 亚洲一级一片aⅴ在线观看| 国产精品久久久久久久久免| 亚洲欧美精品自产自拍| 久久久久网色| 精品国产露脸久久av麻豆 | 久久99热6这里只有精品| 中文天堂在线官网| 秋霞在线观看毛片| 国产精品福利在线免费观看| 国产精品国产三级国产专区5o | 中文字幕人妻熟人妻熟丝袜美| 亚洲怡红院男人天堂| 欧美成人精品欧美一级黄| 国产一区二区在线观看日韩| 内地一区二区视频在线| 久久国产乱子免费精品| 中文字幕亚洲精品专区| 亚洲国产精品专区欧美| 十八禁国产超污无遮挡网站| 长腿黑丝高跟| kizo精华| 天美传媒精品一区二区| 国产成人freesex在线| 亚洲欧美日韩卡通动漫| 亚洲av二区三区四区| 亚洲人成网站在线播| 欧美高清成人免费视频www| 亚洲怡红院男人天堂| 热99re8久久精品国产| 你懂的网址亚洲精品在线观看 | 亚洲国产欧美人成| 欧美丝袜亚洲另类| 国语自产精品视频在线第100页| 有码 亚洲区| 亚洲图色成人| 国产三级在线视频| 久久精品影院6| 久久久久久久国产电影| 日本免费一区二区三区高清不卡| 亚洲va在线va天堂va国产| 亚洲人成网站在线观看播放| 在线免费观看的www视频| 亚洲人与动物交配视频| 国产av码专区亚洲av| 中国国产av一级| 内地一区二区视频在线| 蜜臀久久99精品久久宅男| 性色avwww在线观看| 久久亚洲国产成人精品v| 国产伦精品一区二区三区视频9| 非洲黑人性xxxx精品又粗又长| 免费无遮挡裸体视频| 欧美一区二区国产精品久久精品| 直男gayav资源| 只有这里有精品99| 午夜福利视频1000在线观看| 日韩国内少妇激情av| 男女国产视频网站| 中文字幕人妻熟人妻熟丝袜美| 国产三级在线视频| 18禁裸乳无遮挡免费网站照片| 99久久精品国产国产毛片| 99久久九九国产精品国产免费| 精品国产三级普通话版| 国产在视频线精品| 看十八女毛片水多多多| 两个人视频免费观看高清| 成人欧美大片| 美女黄网站色视频| 级片在线观看| 少妇丰满av| 日本黄色片子视频| 欧美日韩综合久久久久久| 欧美性感艳星| av福利片在线观看| 草草在线视频免费看| 高清av免费在线| 欧美极品一区二区三区四区| 一区二区三区高清视频在线| 久久午夜福利片| 日韩在线高清观看一区二区三区| 国产精品人妻久久久影院| 午夜激情福利司机影院| 午夜福利在线观看免费完整高清在| 国产高清有码在线观看视频| 国产激情偷乱视频一区二区| 亚洲综合色惰| 欧美极品一区二区三区四区| 亚洲精品,欧美精品| 99久国产av精品国产电影| 国语自产精品视频在线第100页| 国产精品女同一区二区软件| 国产欧美日韩精品一区二区| 七月丁香在线播放| 亚洲av电影不卡..在线观看| 亚洲激情五月婷婷啪啪| 精品欧美国产一区二区三| 在线观看66精品国产| 美女cb高潮喷水在线观看| av.在线天堂| 中文字幕精品亚洲无线码一区| 日韩欧美精品v在线| 蜜桃久久精品国产亚洲av| 国产精品,欧美在线| 久久精品熟女亚洲av麻豆精品 | 久久精品影院6| 欧美高清成人免费视频www| 日本av手机在线免费观看| 网址你懂的国产日韩在线| 久久国内精品自在自线图片| 亚洲精品日韩在线中文字幕| 国产精品蜜桃在线观看| 99久久中文字幕三级久久日本| 可以在线观看毛片的网站| 偷拍熟女少妇极品色| 欧美不卡视频在线免费观看| 国语对白做爰xxxⅹ性视频网站| 亚洲精品aⅴ在线观看| 久久99热6这里只有精品| 免费播放大片免费观看视频在线观看 | 亚洲欧美日韩东京热| 男女视频在线观看网站免费| 麻豆久久精品国产亚洲av| 国产视频内射| 久久99热这里只有精品18| 最近最新中文字幕大全电影3| 精品久久久久久久久av| 麻豆成人av视频| 18禁在线无遮挡免费观看视频| 日本免费在线观看一区| 国产中年淑女户外野战色| 3wmmmm亚洲av在线观看| 99热这里只有是精品50| 亚洲av一区综合| 中文字幕亚洲精品专区| 日日啪夜夜撸| 能在线免费看毛片的网站| 熟妇人妻久久中文字幕3abv| 97在线视频观看| 秋霞在线观看毛片| 我的老师免费观看完整版| av视频在线观看入口| 久久这里有精品视频免费| 亚洲av成人av| 午夜精品国产一区二区电影 | 国产成人福利小说| 精品99又大又爽又粗少妇毛片| 2022亚洲国产成人精品| 亚洲欧美中文字幕日韩二区| 高清午夜精品一区二区三区| 成人午夜高清在线视频| 久久久久久久久久成人| 色哟哟·www| 国产免费又黄又爽又色| 免费电影在线观看免费观看| 老师上课跳d突然被开到最大视频| 亚洲伊人久久精品综合 | 人体艺术视频欧美日本| 嫩草影院入口| 亚洲国产精品合色在线| 校园人妻丝袜中文字幕| 久久久久久久久中文| 国产av不卡久久| 18+在线观看网站| 久久午夜福利片| 日韩一本色道免费dvd| 久久草成人影院| 亚州av有码| 亚洲无线观看免费| 99久国产av精品国产电影| 日本与韩国留学比较| 99久久人妻综合| 国产成人福利小说| 欧美激情久久久久久爽电影| 麻豆成人午夜福利视频| 国内揄拍国产精品人妻在线| 国产欧美日韩精品一区二区| 丝袜美腿在线中文| 国产精品综合久久久久久久免费| 老司机福利观看| 天堂√8在线中文| 中文字幕av在线有码专区| 岛国在线免费视频观看| 国产免费一级a男人的天堂| 成人亚洲精品av一区二区| 精品久久久久久久久av| .国产精品久久| 国产乱来视频区| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 国产麻豆成人av免费视频| 午夜a级毛片| 国产激情偷乱视频一区二区| 亚洲婷婷狠狠爱综合网| 亚洲精品,欧美精品| 国产一区有黄有色的免费视频 | 禁无遮挡网站| 亚洲精品成人久久久久久| 国产免费又黄又爽又色| 卡戴珊不雅视频在线播放| 国产三级中文精品| 国产 一区 欧美 日韩| 校园人妻丝袜中文字幕| 欧美一区二区精品小视频在线| 午夜精品国产一区二区电影 | 中文天堂在线官网| 蜜桃久久精品国产亚洲av| 国产视频首页在线观看| 亚洲电影在线观看av| 18禁动态无遮挡网站| 欧美一区二区国产精品久久精品| 欧美日韩综合久久久久久| 干丝袜人妻中文字幕| 日韩三级伦理在线观看| 少妇高潮的动态图| 少妇熟女欧美另类| av又黄又爽大尺度在线免费看 | 欧美一区二区亚洲| 国产免费一级a男人的天堂| 九九爱精品视频在线观看| 少妇的逼好多水| 18+在线观看网站| 亚洲av一区综合| 亚洲天堂国产精品一区在线| 人体艺术视频欧美日本| 老司机影院成人| 国产一级毛片七仙女欲春2| 亚洲欧美一区二区三区国产| 男人舔奶头视频| 国产69精品久久久久777片| 国产精品福利在线免费观看| 少妇人妻精品综合一区二区| 日本与韩国留学比较| 热99re8久久精品国产| 精品人妻视频免费看| 老司机福利观看| 亚洲精品乱久久久久久| 欧美变态另类bdsm刘玥| 亚洲av熟女| 九九在线视频观看精品| 成人性生交大片免费视频hd| 久久99蜜桃精品久久| 国产精品永久免费网站| 亚洲内射少妇av| 国产伦理片在线播放av一区| 69av精品久久久久久| 国产成人精品久久久久久| 午夜激情欧美在线| 亚洲欧美精品专区久久| 别揉我奶头 嗯啊视频| 色哟哟·www| 欧美激情久久久久久爽电影|