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

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

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

0 引言

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

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

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

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

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

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

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

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

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

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

圖1為比選后構(gòu)建的某一個(gè)超導(dǎo)磁體線圈加載電流時(shí)在其內(nèi)部及周邊區(qū)域所形成的磁場(chǎng)分布及Z軸方向上磁場(chǎng)強(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。所得到的磁場(chǎng)能夠在Φ40 mm×60 mm的柱形區(qū)域內(nèi)實(shí)現(xiàn)0.3%的磁場(chǎng)力不均勻度。由圖1（b）可以看到，在Z軸坐標(biāo)-0.263 m附近的區(qū)域內(nèi)存在著一片與重力方向相反的均勻磁場(chǎng)力區(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)呈中心對(duì)稱分布。

圖1 磁場(chǎng)分布及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 不同磁場(chǎng)力非均勻度環(huán)境下的自由液面爬升模擬

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

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

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

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

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

對(duì)Φ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)行劃分；對(duì)Φ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)格劃分；對(duì)Φ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)行對(duì)比，如圖4所示。（a）為對(duì)Φ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)格單元對(duì)Φ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)格無(wú)關(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（對(duì)應(yīng)的磁場(chǎng)力非均勻度分別為0.05%、0.10%、0.30%、1.0%）條件下，磁場(chǎng)力補(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ì)算條件下的液面沿固體壁面的爬升過(guò)程以及磁補(bǔ)償條件與對(duì)應(yīng)的空間微重力條件下的氣液相界面曲線平均距離進(jìn)行了對(duì)比，如圖9所示。

圖9 四種區(qū)域范圍內(nèi)液體沿壁面的爬升過(guò)程及磁補(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可見(jiàn)，在0.05%的非均勻度條件下，液氧在磁補(bǔ)償環(huán)境下的自由界面流動(dòng)能夠很好吻合空間微重力條件（10-6g）下真實(shí)自由界面流動(dòng)，特別是液面爬升階段和形成穩(wěn)定氧氣枕階段的時(shí)間標(biāo)簽都取得了較好的同步對(duì)應(yīng)。

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

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

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

2.4 規(guī)律拓展分析

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

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

3 結(jié)論

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