李 開，柳 軍，劉偉強

(國防科學技術大學航天科學與工程學院，長沙410073)

磁控熱防護系統高溫流場與電磁場耦合計算方法

李 開，柳 軍，劉偉強

(國防科學技術大學航天科學與工程學院，長沙410073)

為了研究真實氣體條件下霍爾效應對磁流體(MHD)控制熱防護效果的影響，建立了熱化學非平衡流場、外加磁場、霍爾電場的耦合計算方法。分析了耦合計算效率與電場更新間隔、電勢殘差收斂極限的關系。給出了采用非平衡霍爾系數模型時電場更新間隔和電勢收斂極限的較優(yōu)取值。研究表明，當霍爾系數較小(β=1.0)時，電場更新間隔大于100流場計算時間步時耦合計算效率較高，且導電壁面和絕緣壁面結論一致。當霍爾系數較大時，耦合計算時間過長，可適當增加電場迭代間隔和電勢收斂極限以提高耦合計算效率。

磁流體控制；熱防護；霍爾電場；耦合計算

0 引 言

臨近飛行器在高超聲速飛行時氣動熱問題日趨嚴重，其“熱障”問題已成為制約飛行器設計的瓶頸[1-3]。最近十年間高超聲速領域的新概念熱防護技術層出不窮。得益于近年來電磁流動控制技術、超導磁鐵技術的進步[4-5]，磁控熱防護技術的現實可行性逐漸增強，應用價值逐年提高，得到了各國研究者的普遍關注[6-7]。但值得注意的是，在典型軌道再入及超軌道再入飛行條件下，霍爾效應明顯，感應產生的霍爾電場會改變前緣流場電流分布，進而改變洛倫茲力，影響磁控熱防護系統的效率。磁控熱防護系統的霍爾效應研究涉及外加磁場、感應電場和非平衡流場三場的耦合計算，是該領域研究的難點[8]。

由于高超聲速飛行器前緣激波后電離度較低，低磁雷諾數假設往往可以得到滿足[9]。此時，外加磁場、感應電場與非平衡流場的耦合可以通過耦合求解電勢泊松方程和非平衡磁流體方程實現。由于大霍爾系數條件下電場求解剛性嚴重、收斂性差[10]，電磁場和熱化學非平衡流場的耦合計算效率較低。在現有文獻的研究中，耦合計算的研究可以分為兩大類：一類側重于構建精度較高的霍爾電場求解方法[11-13];另一類側重于分析霍爾效應對磁控系統的影響效果[10,14-15]。

盡管上述研究大都實現了流場和霍爾電場的耦合計算，但現有方法對于拓撲結構復雜、網格規(guī)律大的流體域耦合計算的效率仍較低。因此，非平衡流場和電場耦合計算效率還需仔細研究。在該松耦合過程中有兩個重要的控制參數，一是電場更新間隔S，即非平衡流場計算S時間步電場開始并完成一次迭代；二是電勢殘差收斂極限ε，即當電勢殘差小于該值時判定電場收斂。不同霍爾系數下電場的收斂性差別較大，給兩參數取值帶來困難。文獻[16]采用的電場更新間隔在10～100之間，具體數值取決于算例的磁場強度。文獻[17]在進行參數霍爾效應研究時，對于不同的霍爾系數采用了同樣的電場更新間隔(S=200)和同樣的電勢收斂極限(ε=1×10-4)?？梢姡F有文獻尚未考慮不同霍爾系數下電場收斂性對耦合效率的影響，也并未確定不同霍爾系數下電場更新間隔和收斂極限的合理取值方法。本文將首先建立三維熱化學非平衡流場和霍爾電場的數值方法以及流場-電磁場之間的耦合計算方法，而后分析在不同霍爾系數下電場更新間隔和收斂極限對耦合計算效率的影響，以便為高超聲速飛行器前緣磁控熱防護系統的霍爾效應研究提供參考。

1 數學模型

1.1 控制方程

針對高超聲速飛行器前緣流動特點，電磁場和流場耦合分析采用低磁雷諾數假設下的三維熱化學非平衡磁流體(Magnetohydrodynamic,MHD)控制方程，如式(1)所示。其中，熱化學非平衡模型采用Park雙溫度模型，化學反應動力學模型采用11組元20反應的Gupta模型。

(1)

式中：U為守恒變量，F、G、H分別為x、y、z方向的對流通量矢量，Fv、Gv、Hv為三方向上的黏性通量項，W為化學反應和振動能量源項矢量，具體表達式見文獻[18]。相對于熱化學非平衡黏流，上述方程增加了一個電磁源項通量SMFD，見下式

SMFD=

(2)

式中：J=(Jx,Jy,Jz)為電流密度矢量，E為電場強度矢量，B=(Bx,By,Bz)為磁感應強度矢量。γ表征不同非平衡模式間的電磁能量分配比，介于0和1之間,這里取γ=0.5[13]。磁感應強度B在給定磁場形態(tài)后即可確定，這里采用磁偶極子磁場[14]。考慮霍爾效應后，結合低磁雷諾數假設可以將對感應場強矢量E的求解轉化為對標量電勢φ的求解。由廣義歐姆定律(3)和電流守恒方程(4)

(3)

▽·J=0

(4)

可以得到關于φ的電勢泊松方程

(5)

(6)

式中：B為磁感應強度矢量B的模。σ為標量電導率，采用式(7)非平衡流電導率模型進行計算。β為霍爾系數。為全面分析不同霍爾系數下耦合效率，采用兩種方法確定β：1)式(7)非平衡模型[14]；2)均勻分布。

(7)

1.2 數值方法

磁流體流動方程(1)離散時，對流項差分采用基于MUSCL插值的AUSMPW格式，隱式求解采用了LU-SGS格式，并且對化學反應、振動以及電磁源項采用了隱式處理以削弱源項過大帶來的剛性從而提高收斂性。

電勢泊松方程(5)的離散基于單元中心有限體積法。采用交替隱式近似因子分解法(ADI-AF)，轉化為下式迭代求解

(8)

其中，系數a1i-1,b1i,c1i+1如下所示

(9)

此外，a2j-1,b2j,c2j+1,a3k-1,b3k,c3k+1形式類似，不再贅述。其中，

(10)

2 耦合方法與邊界條件

流場邊界條件：流動入口為自由來流；出口采用超聲速外推條件；壁面采用全催化等溫壁，壁溫取決于試驗工況?；魻栯妶鲞吔鐥l件：絕緣壁面,J·n=0; 導電壁面,φ=0 V; 其余邊界為Neumann邊界，▽φ·n=0。

3 仿真校驗

3.1 霍爾電場校驗

算例2為間隔電極霍爾效應算例，如圖3所示[13]?！纘的通道邊界內通入+x向的導電流體，電導率σ=1Ω-1m-1。有限寬度的平行電極沿周期性重復布置在通道兩側壁面上。由于通道無限長，圖中的進口和出口邊界設置為周期性邊界。為使通道內流體速度場對電勢場結果無影響，則必須滿足▽×(u×B)=0，因此可假定u=f(y)，B=f(z)，且通道內流動為充分發(fā)展的Poiseuille流動。其中，通道中心線上的流體速度umax=1m/s，通道半寬h=0.5m。

當B=0 T時，兩電極間只存在靜電場。圖4為無霍爾效應時的本文和文獻[13]的電勢等值線對比圖。圖5為有霍爾效應時(β=1.0、B=1T)時的本文與文獻的電流流線對比圖。從圖4～5可以看出，兩種情況下，本文計算結果與文獻結果吻合良好。

3.2 磁流體氣動熱校驗

氣動熱模擬的準確性是檢驗耦合非平衡磁流體計算方法正確性的重要標準。選用日本1996年發(fā)射的軌道再入試驗飛行器(Orbital reentry experiment, OREX)返回艙前體作為計算模型，如圖6所示。選取再入飛行時間在7471.5s(Ma=17.61、H=59.6km)的飛行工況，并采用了有限催化壁面模型(γ=5.0×10-3)進行了氣動熱數值模擬。圖7分別給出了三種外加磁場情況(B0=0.0、0.3、0.5 T)下的平動溫度和壁面熱流分布。通過與Fujino等[20]計算結果的比較可以看出，本文的激波脫體距離、壁面熱流計算結果和文獻[20]吻合良好。驗證了本文外加磁場作用下的非平衡流場及氣動熱數值模擬的準確性。值得一提的是，圖7(a)中平動溫度峰值的差異以及圖7(b)中當B0=0.3 T時在y/L=1.3(L為參考長度，L=1.0139m)附近的壁面熱流差異，可能與本文和文獻采用了不同的平動-振動松弛模型有關。

4 計算工況

為了提高耦合計算效率，需要確定電磁場參數相對于流動參數的最優(yōu)的更新方式，即研究電場更新間隔S、電場殘差收斂極限ε對耦合計算效率的影響。其中耦合計算效率以平均單步迭代時間ta，總收斂時長tc表示。

計算模型網格和第2.2節(jié)相同。采用OREX飛行器7461.5s飛行工況(Ma=20.04，H=63.60km)。駐點磁感應強度B0=0.2 T。流場計算采用當地時間步長，庫朗特數取0.2。鑒于霍爾系數不同對電場收斂性影響巨大，從而會對耦合計算效率產生重要影響。這里根據霍爾系數不同設計了兩組算例：1)均布霍爾系數β=1.0，步進因子ap=0.002；2)Fujino霍爾系數模型，步進因子ap=0.006。前者霍爾系數相對較小，電場收斂快；后者霍爾系數相對較大(β≈5.0～10.0)，電場收斂較慢。耦合計算以忽略霍爾效應的流場收斂解作為初場。

5 結果與分析

表1給出了均布霍爾系數β=1.0時不同電場更新間隔下的單步平均迭代時間和收斂時的總迭代時間。從表1可以看出，和S=10相比，電場更新間隔S≥100時，單步迭代時間和收斂總時間大幅減小，而后隨著電場更新間隔的增加，單步平均迭代時間和收斂總時間雖略微減少但整體變化不大。導電壁面和絕緣壁面規(guī)律一致。

表1 不同電場更新間隔下的迭代時間對比Table 1 Iteration time under different updating intervals of electric field (β=1.0)

圖8(a)和(b)分別為導電壁面不同電場更新間隔下流場和電場的收斂曲線。從圖8可以看出，電場更新間隔對流場收斂曲線形狀影響較小，隨電場更新間隔的增加收斂步數略微減少。電場更新間隔對電勢場收斂影響較大，并且電場更新間隔越大，相同電場迭代步數流場收斂程序越高，相應的電勢收斂所需的電場虛擬時間步越少。

表2給出了采用非平衡霍爾系數模型時不同電場更新間隔下的單步平均迭代時間和收斂總時間。表2還對比了兩種電勢收斂殘差極限下的結果?？梢钥闯?，電勢殘差收斂極限為10-4時，電場收斂耗時較多，單步平均迭代時間較長(S=10時約為58s)，耦合計算時間呈現先降后升的趨勢，最優(yōu)的電場更新間隔為1000。而當電勢殘差收斂極限為10-3時，在S=10～5000范圍內，電場更新間隔越大，單步平均迭代時間和總收斂時間越小，同β=1.0時規(guī)律一致。另外，當S=1000時，兩個收斂極限(ε=10-3、10-4)下的總收斂時間相近，但電場間隔越小二者收斂時間相差越大，S=10時總收斂時間10-4可達10-3的110倍，此時采用相對較大的收斂殘差極限(10-3)可以極大地提高耦合計算效率。

表2 不同電場更新間隔下的迭代時間對比Table 2 Iteration time under different updating intervals of electric field (nonequilibrium Hall parameter model)

5 結 論

以高速飛行器磁控熱防護系統為研究對象，針對其在真實氣體條件下高溫流場和電磁場的耦合求解問題，建立了非平衡流場、外加磁場和霍爾電場的松耦合計算方法，在此基礎上分析了耦合計算效率與電場更新間隔、電勢殘差收斂極限的關系。研究表明，霍爾系數較小(β=1.0)時，電場更新間隔S=100相比S=10單步迭代時間和收斂總時間大幅減少。S>100后隨電場更新間隔的增加略微減少，計算效率變化不大，并且導電壁面和絕緣壁面規(guī)律一致?；魻栂禂递^大時，電場收斂緩慢，耦合計算時間過長，可以通過適當增加電場迭代間隔以及提高電勢收斂極限的方法有效提高耦合計算效率。采用霍爾系數非平衡模型時，電場更新間隔的建議取值為1000，電勢收斂極限可取10-3。

[1] 孟松鶴，楊強，霍施宇，等. 一體化熱防護技術現狀和發(fā)展趨勢[J]. 宇航學報, 2013, 34(10): 1295-1302. [Meng Song-he, Yang Qiang, Huo Shi-yu, et al. State of arts and trend of integrated thermal protection systems [J]. Journal of Astronautics, 2013, 34(10): 1295-1302.]

[2] 李鋒，艾邦成，姜貴慶. 一種熱平衡等溫機制的新型熱防護及相關技術[J]. 宇航學報, 2013, 34(12): 1644-1650. [Li Feng, Ai Bang-cheng, Jiang Gui-qing. A new thermal protection technology based on heat balance isothermal mechanism [J]. Journal of Astronautics, 2013, 34(12): 1644-1650.]

[3] 潘勇,王江峰，伍貽兆. 非結構網格高超聲速MHD流場逆風格式數值模擬[J]. 宇航學報, 2008, 29(1): 104-109. [Pan Yong, Wang Jiang-feng, Wu Yi-zhao. Numerical simulation of hypersonic MHD flows using upwind scheme on unstructured grids [J]. Journal of Astronautics, 2008, 29(1): 104-109.]

[4] 張康平,丁國昊,田正雨,等.磁流體動力學控制二維擴壓器流場數值模擬研究[J].國防科技大學學報, 2009, 31(6):39-41. [Zhang Kang-ping, Ding Guo-hao, Tian Zheng-yu, et al. Numerical simulation of magnetohydrodynamic (MHD) control on 2D diffuser′s flow field [J]. Journal of National University of Defense Technology, 2009, 31(6):39-41.]

[5] 田正雨,張康平,潘沙,等.磁流體動力學斜激波控制數值模擬分析[J]. 力學季刊, 2008, 29(1):72-77. [Tian Zheng-yu, Zhang Kang-ping, Pan Sha, et al. Numerical investigation and analysis for MHD oblique shock control[J]. Chinese Quarterly of Mechanics， 2008, 29(1):72-77. ]

[6] Bityurin V A, Bocharov A N. Study of catalytic effects at reentry vehicle [C].The 52nd Aerospace Sciences Meeting, National Harbor, Maryland, USA, January 13-17, 2014.

[7] Cristofolini A, Borghi C A, Neretti G, et al. MHD interaction around a blunt body in a hypersonic unseeded air flow[C].The 18th AIAA/3AF International Space Planes and Hypersonic Systems and Technologies Conference, Tours, France, September 24-28, 2012.

[8] 劉深深,桂業(yè)偉,唐偉,等.一種多場耦合數據傳遞新方法[J]. 宇航學報, 2016, 37(1):61-67. [Liu Shen-shen, Gui Ye-wei, Tang Wei, et al. A new data transfer method in fluid-thermal-structrure coupling problems [J]. Journal of Astronautics, 2016, 37(1): 61-67.]

[9] 呂浩宇, 李椿萱.Hall效應對磁流體壓縮管道電磁流動特性的影響[J].科學通報, 2010, 55(12):1182-1188. [Lv Hao-yu, Li Chun-xuan. Influence of Hall effects on characteristics of magnetohydrodynamic converging channel [J]. Chin. Sci. Bull, 2010, 55(12):1182-1188.]

[10] 胡海洋, 楊云軍, 周偉江. 大霍爾系數下電離氣體與磁場相互作用規(guī)律數值研究[J]. 力學學報，2011, 43(3): 453-460. [Hu Hai-yang, Yang Yun-jun, Zhou Wei-jiang. Numerical research on the interaction between ionized gas and magnetic field under high Hall parameter [J]. Chinese Journal of Theoretical and Applied Mechanics, 2011, 43(3):453-460.]

[11] Gaitonde D V, Poggie J. Elements of a numerical procedure for 3-D MGD flow control analysis [C]. 40th AIAA Aerospace Sciences Meeting and Exhibit, Reno, USA, January 14-17, 2002.

[12] Wan T, Candler G V, Macheret S O, et al. Three-dimensional simulation of the electric field and magnetohydrodynamic power generation during reentry [J]. AIAA Journal, 2009, 47(6):1327-1336.

[13] Bisek N J. Numerical study of plasma-assisted aerodynamic control for hypersonic vehicles [D]. Michigan: The University of Michigan, 2010.

[14] Fujino T, Matsumoto Y, Kasahara J, et al. Numerical studies of magnetohydrodynamic flow control considering real wall electrical conductivity [J]. Journal of Spacecraft and Rockets, 2007, 44(3):625-632.

[15] 呂浩宇,李椿萱,董海濤. 三維超聲速磁流體發(fā)生器的流動特性[J]. 中國科學G輯,2009,39(3):435-445. [Lv Hao-yu, Li Chun-xuan, Dong Hai-tao. Characterization of the three-dimensional supersonic flow for the MHD generator [J]. Science in China Series G-Physics Mechan. Astron. , 2009, 39(3):435-445.]

[16] Bisek N J, Gosse R, Poggie J. Computational study of impregnated ablator for improved magnetohydrodynamic heat shield [J]. Journal of Spacecraft and Rockets, 2013, 50(5): 927-935.

[17] Otsu H, Konigorski D, Abe T. Influence of Hall effect on electrodynamic heat shield system for reentry vehicles [J]. AIAA Journal, 2010, 48(10):2177-2186.

[18] 柳軍.熱化學非平衡流及其輻射現象的實驗和數值計算研究[D].長沙：國防科技大學，2004. [Liu Jun. Experimental and numerical research on thermo-chemical nonequilibrium flow with radiation phenomenon [D]. Changsha: National University of Defense Technology, 2004.]

[19] Gnoffo P A, Gupta R N, Shinn J L. Conservation equations and physical models for hypersonic air flows in thermal and chemical nonequilibrium [R]. NASA, TP-2867, 1989.

[20] Fujino T, Ishikawa M. Numerical simulation of control of plasma flow with magnetic field for thermal protection in earth reentry flight [J]. IEEE Transactions on Plasma Science, 2006,34(2):409-420.

通信地址：湖南省長沙市國防科技大學航天科學與工程學院(410073)

電話：15616246843

E-mail:likai898989@126.com

(編輯：牛苗苗)

Numerical Methods of Coupling High Temperature Flow Field with Electro Magnetic Field for Magnetohydrodynamic Heat Shield System

LI Kai, LIU Jun, LIU Wei-qiang

(College of Aerospace Science and Engineering, National University of Defense Technology, Changsha 410073, China)

Under the real flight conditions, the Hall effect usually cannot be neglected which may affect the performance of the magnetohydrodynamics(MHD) thermal protection. Therefore, the coupled numerical simulations are conducted for the thermochemical nonequilibrium flow field, the externally applied magnetic field and the induced electric field around a typical hypersonic vehicle. After that, the relations between the coupled calculation efficiency and the updating interval of the electric field as well as the converging residual limit of the electric potential are analyzed. Results show that the coupled calculation efficiency is higher under the condition that the updating interval is greater than 100 flow steps during a relatively smaller Hall parameter. However, during a larger Hall parameter, the coupling efficiency is lower but can be improved by increasing the updating interval of the electric field and potential residual limit. The suggested updating interval and potential residual limit are also given while the nonequilibrium Hall parameter model is utilized.

MHD flow control; Thermal protection; Hall electric field; Coupled computation

2016-09-26;

2017-02-27

國家自然科學基金(90916018)；湖南省自然科學基金重點資助項目(13JJ2002)

V211.1

A

1000-1328(2017)05-0474-07

10.3873/j.issn.1000-1328.2017.05.005

李 開(1989-)，男，博士生，主要從事磁流體流動控制，高溫氣體動力學以及飛行器熱防護系統設計方面的研究。