磁場力非均勻度對液氧磁補償微重力自由界面的影響

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

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

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

0 引言

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

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

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

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

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

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

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

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

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

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

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

圖1 磁場分布及Z軸方向上磁感應強度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磁補償區(qū)域均勻性表征圖Fig.2 Characterization of uniformity of Φ40 mm×60 mm magnetic compensation region

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

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

2.1 網(wǎng)格劃分及無關性驗證

利用COMSOL Multiphysics軟件對磁場、流場及多相流三種物理場在二維軸對稱坐標下利用有限元方法進行了耦合計算。其中流場為層流，多相流采用水平集（Level Set）方法進行計算，物理場離散方法為有限元方法，模型中心軸線為軸對稱邊界條件，其余三個邊界均設為壁面邊界條件，且液體與壁面之間的接觸角為15°。計算容差10-5，計算時間步長為計算殘差控制的變時間步長。表面張力用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ū)域進行了網(wǎng)格無關性驗證。為了能定量比較不同計算條件下得到的氣液相界面，引入RMS均方根計算不同計算條件下氣液相界面曲線之間的平均距離：

對Φ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)格單元進行劃分；對Φ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)格單元進行網(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)格單元進行網(wǎng)格劃分，并將計算結(jié)果分別與最密的網(wǎng)格劃分計算的結(jié)果進行對比，如圖4所示。（a）為對Φ17 mm×25.5 mm的不同網(wǎng)格單元大小的計算結(jié)果與0.000 1 m網(wǎng)格單元的計算結(jié)果之間的均方根；（b）為Φ24 mm×36 mm的不同網(wǎng)格單元大小的計算結(jié)果與0.000 3 m網(wǎng)格單元的計算結(jié)果之間的均方根；（c）為Φ40 mm×60 mm的不同網(wǎng)格單元大小的計算結(jié)果與0.000 3 m網(wǎng)格單元的計算結(jié)果之間的均方根；（d）為Φ70 mm×105 mm不同網(wǎng)格單元大小的計算結(jié)果與0.000 5 m網(wǎng)格單元計算結(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的計算區(qū)域進行劃分比較合理。

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

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

2.2 計算結(jié)果

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

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

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

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

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

2.3 結(jié)果分析

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

圖9 四種區(qū)域范圍內(nèi)液體沿壁面的爬升過程及磁補償條件與真實微重力條件偏離情況圖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越小意味著磁補償條件下的液面曲線與空間微重力條件下的液面曲線越相似，即更準確的補償效果。由圖9可見，在0.05%的非均勻度條件下，液氧在磁補償環(huán)境下的自由界面流動能夠很好吻合空間微重力條件（10-6g）下真實自由界面流動，特別是液面爬升階段和形成穩(wěn)定氧氣枕階段的時間標簽都取得了較好的同步對應。

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

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

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

2.4 規(guī)律拓展分析

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

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

3 結(jié)論

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