張澤宇，黃永華*，梁益濤，耑 銳，張 亮，卜劭華

（1.上海交通大學(xué)制冷與低溫工程研究所，上海 200240；

2.上海宇航系統(tǒng)工程研究所，上海 201108）

0 引言

在空間微重力環(huán)境中，表面張力相對于自然對流占主導(dǎo)作用，因此流體的流動(dòng)和傳熱特性與地面常重力條件下存在顯著區(qū)別，這對低溫推進(jìn)劑的在軌貯存和管理技術(shù)等帶來了挑戰(zhàn)。進(jìn)行微重力實(shí)驗(yàn)的常規(guī)研究手段有落塔[1]、拋物飛機(jī)[2]和探空火箭等，其中最完美的手段是在空間站等真實(shí)微重力場所進(jìn)行實(shí)驗(yàn)。落塔實(shí)驗(yàn)通常整個(gè)過程都不到5 s，而其他手段費(fèi)時(shí)費(fèi)力，成本高昂。此外，相比于常溫流體，進(jìn)行低溫流體微重力實(shí)驗(yàn)挑戰(zhàn)性更大。近20年來，磁重力補(bǔ)償成為這一領(lǐng)域的新方向，特別是物理科學(xué)和生命科學(xué)實(shí)驗(yàn)，如對液氦[3]、液氫[4]、液氧[5]、水[6]甚至玻璃等介質(zhì)的相關(guān)實(shí)驗(yàn)。在力平衡理論方面，從麥克斯韋方程出發(fā)，引入場向量G=grad（B2）（其中B為磁感應(yīng)強(qiáng)度）對磁重力補(bǔ)償?shù)臈l件進(jìn)行描述，以便分析磁力場分布特征[7]，從而指導(dǎo)所需的磁場設(shè)計(jì)。特別需要指出的是，由于液氧本身具有順磁性特征，且其磁化系數(shù)高出其他低溫流體（液氫、液氮等）幾個(gè)數(shù)量級，同時(shí)液氧又正好是低溫推進(jìn)劑的氧化劑，因而最適合用來進(jìn)行磁補(bǔ)償微重力實(shí)驗(yàn)。

國際上相關(guān)研究機(jī)構(gòu)已對低溫流體的磁補(bǔ)償微重力模擬開展了一定的實(shí)驗(yàn)研究。如美國國家強(qiáng)磁場研究中心（NHFML）在1999年首先利用超導(dǎo)磁體對直徑Φ10 mm的球形區(qū)域內(nèi)進(jìn)行了0.1 g水平的液氦液滴微重力補(bǔ)償實(shí)驗(yàn)觀察[8]。日本筑波大學(xué)磁場研究中心于2002年通過數(shù)值模擬的方法對超導(dǎo)線圈的結(jié)構(gòu)進(jìn)行了分析設(shè)計(jì)，搭建了一套磁補(bǔ)償微重力實(shí)驗(yàn)臺，在Φ10 mm×20 mm區(qū)域的柱形區(qū)域內(nèi)實(shí)現(xiàn)了0.01 g水平的微重力補(bǔ)償，并進(jìn)行了水等流體的微重力實(shí)驗(yàn)[9]。法國原子能總署（CEA）于2004年搭建了針對液氫的微重力補(bǔ)償實(shí)驗(yàn)臺，并對液氫進(jìn)行了磁補(bǔ)償微重力環(huán)境下臨界流動(dòng)實(shí)驗(yàn)[10]，而后又對實(shí)驗(yàn)臺進(jìn)行了改進(jìn)，在Φ30 mm×50 mm的柱形區(qū)域內(nèi)達(dá)到了0.04 g的微重力補(bǔ)償水平，并對液氧開展了微重力池沸騰[11]及變重力加速度[12]等相關(guān)實(shí)驗(yàn)。

由此可見，盡管國外已對低溫流體開展了一定的磁補(bǔ)償微重力實(shí)驗(yàn)，但都只能達(dá)到0.01 g量級的微重力水平（可能與其所利用的超導(dǎo)磁體結(jié)構(gòu)單一有關(guān)），且基本上未涉及對磁場均勻性，特別是關(guān)于磁場均勻性是如何影響磁重力補(bǔ)充水平方面的闡述。另一方面，國內(nèi)在磁補(bǔ)償實(shí)現(xiàn)低溫流體微重力方面尚未開展相關(guān)研究。本文將利用磁場、流體速度場等耦合的數(shù)值模擬方法，建立大范圍、高水平的液氧微重力磁補(bǔ)償環(huán)境，對不同磁場力均勻度條件下的氣液兩相界面在容器內(nèi)的爬升過程進(jìn)行模擬，并且與空間真實(shí)微重力環(huán)境（10-6g）進(jìn)行對比，以確定液氧磁補(bǔ)償微重力實(shí)驗(yàn)所依賴的超導(dǎo)磁體產(chǎn)生磁場力水平的技術(shù)要求，為后續(xù)的超導(dǎo)磁體微重力實(shí)驗(yàn)臺的設(shè)計(jì)與搭建提供數(shù)據(jù)支持。

1 磁補(bǔ)償微重力環(huán)境的實(shí)現(xiàn)

在磁場環(huán)境中，順磁性物質(zhì)液氧受到磁體積力的作用，因而如果液氧所受到的磁場力的方向與重力方向相反，理論上便能夠得到等效的微重力環(huán)境。根據(jù)磁化理論，磁性流體在磁場中所受到的磁體積力fm為：

根據(jù)力平衡，則實(shí)現(xiàn)利用磁場力對重力g進(jìn)行完全磁補(bǔ)償?shù)膶?shí)現(xiàn)條件為：

式中：χ為磁化系數(shù)；μ0為真空磁導(dǎo)率；ρ為流體的密度。對于液氧，其磁化系數(shù)χ=3.45×10-3，密度ρ=1 141.8 kg/m3，則對應(yīng)1 g所需的G0=8.15 T2/m。因而，若要得到均勻的磁場力環(huán)境，需要使得磁場強(qiáng)度B平方的梯度在一定范圍內(nèi)盡可能均勻。微重力的實(shí)現(xiàn)水平與該梯度均勻度密切相關(guān)，梯度越均勻，微重力水平越高。

通過超導(dǎo)磁體所實(shí)現(xiàn)的空間磁場力分布是不均勻的，均勻度越差意味著離理想的微重力環(huán)境相差越大。本文提出以下非均勻度指標(biāo)對磁場力的非均勻度分布進(jìn)行定量分析：

式中：εz為磁場力在Z軸方向上補(bǔ)償?shù)姆蔷鶆蚨?，εz越低意味著磁場力在Z軸方向上與重力越接近，即越接近空間微重力環(huán)境。εr為磁場力在R軸方向上補(bǔ)償?shù)姆蔷鶆蚨取o定的工質(zhì)，對磁感應(yīng)強(qiáng)度B進(jìn)行二階Tylor展開并利用麥克斯韋方程div（B→）=0和rot（B→B→）=0得到：

式中：B為磁感應(yīng)強(qiáng)度；R為磁補(bǔ)償區(qū)域的特征長度。對于給定的磁體結(jié)構(gòu)和實(shí)驗(yàn)工質(zhì)，不同的補(bǔ)償范圍對應(yīng)著一組εz和εr的值。因此選取不同的補(bǔ)償區(qū)域作為對象時(shí)，都可以計(jì)算出各自的εz，用于評估研究磁場力補(bǔ)償?shù)姆蔷鶆蚨葘ξ⒅亓?shí)驗(yàn)的影響。

圖1為比選后構(gòu)建的某一個(gè)超導(dǎo)磁體線圈加載電流時(shí)在其內(nèi)部及周邊區(qū)域所形成的磁場分布及Z軸方向上磁場強(qiáng)度B和物理量G的分布曲線。超導(dǎo)線圈由兩組線圈組成，外線圈內(nèi)徑250 mm、外徑263.03 mm、高250.44 mm，內(nèi)線圈內(nèi)徑266.72 mm、外徑272.03 mm、高924.48 mm，內(nèi)外線圈的工作電流均為135 A。所得到的磁場能夠在Φ40 mm×60 mm的柱形區(qū)域內(nèi)實(shí)現(xiàn)0.3%的磁場力不均勻度。由圖1（b）可以看到，在Z軸坐標(biāo)-0.263 m附近的區(qū)域內(nèi)存在著一片與重力方向相反的均勻磁場力區(qū)域，該區(qū)域即為液氧的磁補(bǔ)償理想?yún)^(qū)域。圖2為超導(dǎo)磁體在Φ40 mm×60 mm區(qū)域內(nèi)磁補(bǔ)償?shù)摩舲和εr的分布。由圖2（a）能夠看到εz沿補(bǔ)償點(diǎn)呈中心對稱分布。

圖1 磁場分布及Z軸方向上磁感應(yīng)強(qiáng)度B和物理量G的分布圖Fig.1 Magnetic field distribution and distribution of magnetic induction intensity and physical quantity G in Z-direction

圖2 Φ40 mm×60 mm磁補(bǔ)償區(qū)域均勻性表征圖Fig.2 Characterization of uniformity of Φ40 mm×60 mm magnetic compensation region

2 不同磁場力非均勻度環(huán)境下的自由液面爬升模擬

在上述模型和評價(jià)指標(biāo)的基礎(chǔ)上，選取不同大小的磁場力補(bǔ)償區(qū)域，對區(qū)域內(nèi)的液氧在0.1 MPa下的箱內(nèi)自由液面爬升過程進(jìn)行數(shù)值模擬，由此評估不同的磁場力非均勻度（0.05%、0.10%、0.30%、1.0%）工況與對應(yīng)的空間真實(shí)微重力條件工況的特性差異。

2.1 網(wǎng)格劃分及無關(guān)性驗(yàn)證

利用COMSOL Multiphysics軟件對磁場、流場及多相流三種物理場在二維軸對稱坐標(biāo)下利用有限元方法進(jìn)行了耦合計(jì)算。其中流場為層流，多相流采用水平集（Level Set）方法進(jìn)行計(jì)算，物理場離散方法為有限元方法，模型中心軸線為軸對稱邊界條件，其余三個(gè)邊界均設(shè)為壁面邊界條件，且液體與壁面之間的接觸角為15°。計(jì)算容差10-5，計(jì)算時(shí)間步長為計(jì)算殘差控制的變時(shí)間步長。表面張力用CSF模型[13]：

式中：σ為表面張力系數(shù)，取13.2 mN/m；δ為界面delta函數(shù)；κ為界面曲率；n為界面單位向量。

模擬模型網(wǎng)格劃分如圖3所示，其中（a）為流場區(qū)域網(wǎng)格的放大，流場區(qū)域和磁體結(jié)構(gòu)由四邊形網(wǎng)格劃分，其他區(qū)域由三角形網(wǎng)格劃分。由于不同大小的流場區(qū)域需要不同的網(wǎng)格數(shù)量，因此分別對Φ17 mm×25.5 mm、Φ24 mm×36 mm、Φ40 mm×60 mm、Φ70 mm×105 mm的流場區(qū)域進(jìn)行了網(wǎng)格無關(guān)性驗(yàn)證。為了能定量比較不同計(jì)算條件下得到的氣液相界面，引入RMS均方根計(jì)算不同計(jì)算條件下氣液相界面曲線之間的平均距離：

對Φ17 mm×25.5 mm的區(qū)域分別按0.000 75 m、0.000 5 m、0.000 375 m、0.000 25 m和0.000 1 m的網(wǎng)格單元進(jìn)行劃分；對Φ24 mm×36 mm和Φ40 mm×60 mm的區(qū)域分別按0.001 m、0.000 75 m、0.000 5 m、0.000 4 m和0.000 3 m的網(wǎng)格單元進(jìn)行網(wǎng)格劃分；對Φ70 mm ×105 mm的區(qū)域按0.003 m、0.002 m、0.001 m、0.000 75 m和0.000 5 m的網(wǎng)格單元進(jìn)行網(wǎng)格劃分，并將計(jì)算結(jié)果分別與最密的網(wǎng)格劃分計(jì)算的結(jié)果進(jìn)行對比，如圖4所示。（a）為對Φ17 mm×25.5 mm的不同網(wǎng)格單元大小的計(jì)算結(jié)果與0.000 1 m網(wǎng)格單元的計(jì)算結(jié)果之間的均方根；（b）為Φ24 mm×36 mm的不同網(wǎng)格單元大小的計(jì)算結(jié)果與0.000 3 m網(wǎng)格單元的計(jì)算結(jié)果之間的均方根；（c）為Φ40 mm×60 mm的不同網(wǎng)格單元大小的計(jì)算結(jié)果與0.000 3 m網(wǎng)格單元的計(jì)算結(jié)果之間的均方根；（d）為Φ70 mm×105 mm不同網(wǎng)格單元大小的計(jì)算結(jié)果與0.000 5 m網(wǎng)格單元計(jì)算結(jié)果之間的均方根。結(jié)果表明，分別選取0.000 25 m、0.000 5 m、0.000 5 m、0.001 m的網(wǎng)格單元對Φ17 mm×25.5 mm、Φ24 mm×36 mm、Φ40 mm×60 mm、Φ70mm×105mm的計(jì)算區(qū)域進(jìn)行劃分比較合理。

圖3 數(shù)值模擬模型網(wǎng)格劃分圖Fig.3 Mesh generation of unmerical simulation model

圖4 不同計(jì)算模型的網(wǎng)格無關(guān)性檢驗(yàn)圖Fig.4 Grid independence test of different computing models

2.2 計(jì)算結(jié)果

圖5 、圖6、圖7和圖8分別為補(bǔ)償區(qū)域Φ17 mm×25.5 mm、Φ24 mm×36 mm、Φ40 mm×60 mm、Φ70 mm×105 mm（對應(yīng)的磁場力非均勻度分別為0.05%、0.10%、0.30%、1.0%）條件下，磁場力補(bǔ)償（a）和空間零重力環(huán)境（b）液氧氣液相界面的瞬態(tài)變化。

圖5 0.05%非均勻度條件時(shí)磁補(bǔ)償和真實(shí)微重力環(huán)境下氧氣液相界面演化圖Fig.5 Evolution of oxygen liquid interface under magnetic compensation and real microgravity at 0.05%inhomogeneity

圖6 0.10%非均勻度條件時(shí)磁補(bǔ)償和真實(shí)微重力環(huán)境下氧氣液相界面演化圖Fig.6 Evolution of oxygen liquid interface under magnetic compensation and real microgravity at 0.10%inhomogeneity

圖7 0.30%非均勻度條件時(shí)磁補(bǔ)償和真實(shí)微重力環(huán)境下氧氣液相界面演化圖Fig.7 Evolution of oxygen liquid interface under magnetic compensation and real microgravity at 0.30%inhomogeneity

圖8 1.0%非均勻度條件時(shí)磁補(bǔ)償和真實(shí)微重力環(huán)境下氧氣液相界面演化圖Fig.8 Evolution of oxygen liquid interface under magnetic compensation and real microgravity at 1.0%inhomogeneity

2.3 結(jié)果分析

基于圖5～8中的計(jì)算結(jié)果，將四種計(jì)算條件下的液面沿固體壁面的爬升過程以及磁補(bǔ)償條件與對應(yīng)的空間微重力條件下的氣液相界面曲線平均距離進(jìn)行了對比，如圖9所示。

圖9 四種區(qū)域范圍內(nèi)液體沿壁面的爬升過程及磁補(bǔ)償條件與真實(shí)微重力條件偏離情況圖Fig.9 The climbing process of liquid along the wall in four regions and the deviation of magnetic compensation condition from the real microgravity condition

RMS越小意味著磁補(bǔ)償條件下的液面曲線與空間微重力條件下的液面曲線越相似，即更準(zhǔn)確的補(bǔ)償效果。由圖9可見，在0.05%的非均勻度條件下，液氧在磁補(bǔ)償環(huán)境下的自由界面流動(dòng)能夠很好吻合空間微重力條件（10-6g）下真實(shí)自由界面流動(dòng)，特別是液面爬升階段和形成穩(wěn)定氧氣枕階段的時(shí)間標(biāo)簽都取得了較好的同步對應(yīng)。

0.10%的非均勻度條件下，液氧在磁補(bǔ)償環(huán)境下的自由界面流動(dòng)無論是液面爬升階段還是形成氧氣枕階段都得到了一定程度的還原，因此0.1%的非均勻度可視為能基本滿足要求的磁補(bǔ)償環(huán)境條件，優(yōu)于此技術(shù)指標(biāo)的工況將取得更佳的微重力模擬水平。

0.30%的非均勻度條件下，液氧在磁補(bǔ)償環(huán)境下的自由界面流動(dòng)能夠部分還原空間微重力條件（10-6g）下的自由界面流動(dòng)情況。在液面爬升階段，磁補(bǔ)償環(huán)境較好地還原了液氧液面沿固體壁面的爬升，但未能完全還原液氧底部隨著爬升過程而出現(xiàn)的振蕩現(xiàn)象。在形成穩(wěn)定氧氣枕階段，磁補(bǔ)償環(huán)境沒能正確還原空間微重力環(huán)境下的液氧氣泡?？梢?，雖然不均勻度僅僅從0.10%惡化到0.30%，所帶來的微重力模擬水平卻出現(xiàn)了明顯下滑。

1.0%以內(nèi)的非均勻度條件下，液氧在磁補(bǔ)償環(huán)境下的箱內(nèi)自由界面流動(dòng)基本上不再能夠還原空間微重力條件（10-6g）下的自由界面流動(dòng)。無論是在液面爬升階段還是形成穩(wěn)定氧氣枕階段，液氧液面均與空間微重力環(huán)境下的真實(shí)爬升速度、形狀出現(xiàn)了明顯的差異。由于非均勻磁場力的干擾，氧氣泡會在液體內(nèi)不斷地變形運(yùn)動(dòng)而無法形成穩(wěn)定的氧氣枕。

2.4 規(guī)律拓展分析

為了研究更一般性的規(guī)律，通過改變磁補(bǔ)償區(qū)域的大小，對0.05%～1.0%共11個(gè)工況的不同非均勻度磁場力下的氣液相界面變化過程進(jìn)行了模擬。圖10為曲線平均距離RMS在5 s計(jì)算時(shí)間內(nèi)的平均值隨磁場力非均勻度的變化情況，反映了真實(shí)微重力環(huán)境與磁補(bǔ)償微重力環(huán)境下計(jì)算得到的液面曲線間的差異隨磁場力非均勻度的變化情況。空間微重力環(huán)境下與磁補(bǔ)償環(huán)境下氣液相界面之間的相似程度僅在0.10%以內(nèi)的補(bǔ)償情況下有了明顯的降低，說明在0.10%以內(nèi)的磁場力非均勻度條件下，液體沿壁面的爬升過程以及包裹氣體形成氣泡的過程開始得到較好的還原，而在更低的非均勻度條件下（0.05%）則能夠得到很高的真實(shí)微重力狀態(tài)還原度。

圖10 平均RMS隨磁場力非均勻度的變化Fig.10 Variation of mean RMS with inhomogeneity of magnetic field force

3 結(jié)論

本文對順磁性流體液氧在磁場力作用下實(shí)現(xiàn)微重力懸浮的自由界面流動(dòng)過程進(jìn)行了磁場、流場、多相流耦合的數(shù)值模擬，研究了磁場力的非均勻度對液氧磁補(bǔ)償微重力自由界面的影響。提出了評估磁場力非均勻性的定量指標(biāo)和評估方法，通過與空間真實(shí)均勻微重力環(huán)境的比對，發(fā)現(xiàn)磁補(bǔ)償微重力環(huán)境的制造水平對磁場力的非均勻度很敏感。在0.05%的非均勻度磁補(bǔ)償條件下，不管是液氧液面爬升過程還是形成氣枕的過程，均可實(shí)現(xiàn)很高的真實(shí)微重力狀態(tài)還原度；在0.1%的非均勻度磁補(bǔ)償條件下，磁補(bǔ)償微重力環(huán)境初步具備還原空間微重力環(huán)境的能力；非均勻度惡化到0.3%時(shí)，僅能對液氧的液面爬升過程進(jìn)行可接受的模擬，其液體底部的振蕩過程和氧氣枕形成過程都出現(xiàn)明顯的偏差；當(dāng)進(jìn)一步惡化到1.0%的非均勻度時(shí)，所有過程均出現(xiàn)了明顯差異，甚至氣液相界面會無法穩(wěn)定，基本判定已不適合用于空間微重力真實(shí)狀態(tài)模擬。上述研究結(jié)論有助于確定液氧磁補(bǔ)償微重力實(shí)驗(yàn)所依賴的超導(dǎo)磁體產(chǎn)生磁場力水平的技術(shù)要求，為后續(xù)的超導(dǎo)磁體微重力實(shí)驗(yàn)臺的設(shè)計(jì)與搭建提供數(shù)據(jù)支持。