基于電流積分計(jì)算磁矢量勢修正的低磁雷諾數(shù)方法

丁明松 江濤 劉慶宗 董維中 高鐵鎖 傅楊?yuàn)W驍

(中國空氣動(dòng)力研究與發(fā)展中心計(jì)算空氣動(dòng)力研究所,綿陽 621000)

針對低磁雷諾數(shù)方法的適用性問題,分析了當(dāng)前低磁雷諾數(shù)條件應(yīng)用上存在分歧以及全磁流體力學(xué)方法在高超聲速領(lǐng)域局限性產(chǎn)生的原理.在低磁雷諾磁流體力學(xué)控制數(shù)值模擬方法的基礎(chǔ)上,基于感應(yīng)電流積分計(jì)算磁矢量勢,考慮截?cái)嘁蜃訉τ?jì)算域的縮減,提出了一種考慮感應(yīng)磁場修正的低磁雷諾數(shù)磁流體力學(xué)計(jì)算方法,并加以驗(yàn)證.結(jié)合RAM-C鈍錐體試驗(yàn)飛行狀態(tài)數(shù)值模擬,分析了“忽略感應(yīng)磁場”造成的計(jì)算偏差,探討了“低磁雷諾數(shù)假設(shè)”在高超聲速領(lǐng)域的使用原則.研究表明: 1)本文建立的修正計(jì)算方法,突破低磁雷諾數(shù)條件的限制,拓展了低磁雷諾數(shù)方法在高超聲速領(lǐng)域的適用性和應(yīng)用范圍,數(shù)值模擬結(jié)果可信度高,同時(shí)通過積分區(qū)域限制等方法使計(jì)算效率得到了較大的提升;2)高超聲速流動(dòng)過程中感應(yīng)磁場的影響,在宏觀上表現(xiàn)為對外加磁場的削弱和扭曲,一定程度上降低了磁控效果;本文計(jì)算條件下,“Rem <0.1”的低磁雷諾數(shù)條件可能過于保守,建議取為Rem <1.0,同時(shí)其特征電導(dǎo)率和特征尺度應(yīng)綜合考慮實(shí)際的等離子體分布.

1 引 言

早在二十世紀(jì)五六十年代,人們就開始了磁流體力學(xué) (magnetohydrodynamic,MHD)的研究.上世紀(jì)九十年代以來,隨著高超聲速飛行技術(shù)、計(jì)算機(jī)技術(shù)和便攜磁體制造業(yè)的發(fā)展,高超聲速磁流體力學(xué)控制掀起了新的研究熱潮.在飛行器不同部位加載合適的磁場發(fā)生裝置,通過動(dòng)量和能量的輸運(yùn)對高超聲速飛行器周圍的流動(dòng)狀態(tài)進(jìn)行控制,顯著地改進(jìn)和提升飛行器的氣動(dòng)力/熱特性、電磁通信性能、隱身與突防特性等,這就是高超聲速磁流體力學(xué)控制[1].高超聲速磁流體控制技術(shù)已經(jīng)成為一個(gè)多學(xué)科交叉融合的重要研究方向,其應(yīng)用前景十分廣泛,可用于飛行器氣動(dòng)熱環(huán)境管理、氣動(dòng)力操控、等離子體電子數(shù)密度控制、流速控制和邊界層控制等多個(gè)方面[2].由于它具備不改變飛行器外部結(jié)構(gòu)、無需機(jī)械傳動(dòng)、響應(yīng)快速準(zhǔn)確、兼容性強(qiáng)等優(yōu)點(diǎn),因而受到世界各國高度重視.

數(shù)值模擬是研究高超聲速磁流體力學(xué)控制的主要手段之一.從高超聲速磁流體力學(xué)的國內(nèi)外發(fā)展?fàn)顩r可以看出,當(dāng)前高超聲速M(fèi)HD數(shù)值模擬的方法主要有兩種[1,3]: 全MHD方法和低磁雷諾數(shù)MHD方法(或簡稱低磁雷諾數(shù)方法).

全MHD方法[4,5]是通過廣義歐姆定律和安培環(huán)路定理消去流動(dòng)控制方程和麥克斯韋方程組中電場E和感應(yīng)電流J,然后將流動(dòng)控制方程和磁擴(kuò)散方程(由法拉第定律得到)聯(lián)立構(gòu)成全MHD方程組,求解該方程組得到磁流體流場和感應(yīng)磁場分布.在此基礎(chǔ)上,可通過求解電流連續(xù)性方程和廣義歐姆定律,得到電場和電流分布.全MHD方法理論上適用于各類磁流體的數(shù)值計(jì)算分析,但在高超聲速飛行器流場(尤其是外流場)MHD數(shù)值模擬方面存在很大困難[1,2].

與其他領(lǐng)域不同,高超聲速飛行領(lǐng)域磁流體力學(xué)控制具有其自身的特點(diǎn).由于高超聲速飛行器外流場等離子體磁導(dǎo)率小、導(dǎo)電性一般較弱,電磁場變化的特征時(shí)間遠(yuǎn)小于流動(dòng)的特征時(shí)間,因此數(shù)值模擬飛行器高超聲速流場磁流體流動(dòng)時(shí)常采用“低磁雷諾數(shù)近似”: 在磁導(dǎo)率、電導(dǎo)率較小時(shí),感應(yīng)磁場相對于外加磁場很弱,基本上可以忽略,此時(shí)可以假設(shè)磁場未受流動(dòng)干擾,舍去電磁交叉項(xiàng),使數(shù)值模擬過程得到簡化,這就是低磁雷諾數(shù)MHD方法.對于穩(wěn)態(tài)(或準(zhǔn)穩(wěn)態(tài))高超聲速磁流體流動(dòng),采用低磁雷諾數(shù)方法,可以有效地避免偽磁場散度和低電導(dǎo)率帶來的剛性問題,減小了非必要的繁復(fù)計(jì)算,這是業(yè)內(nèi)通行的做法之一.國內(nèi)外很多高超飛行器磁流體力學(xué)控制方面的研究,如日本的Nagata等[6]、Otsu[7]、Fujino 等[8]和 Takahashi等[9],美國的Bisek等[10,11]和Lee等[12],意大利的Cristofolini等[13,14],國內(nèi)的李開等[15,16]、姚霄等[17]、潘勇[2]、張向洪[3]、何淼生等[18]、陳剛等[19]以及孫曉輝[20]等,在針對“低磁導(dǎo)率、低電導(dǎo)率特征”氣體的數(shù)值模擬時(shí),也幾乎都采用低磁雷諾數(shù)MHD方法.

盡管低磁雷諾數(shù)MHD方法,在高超聲速飛行領(lǐng)域得到了較為廣泛的應(yīng)用和發(fā)展,但由于它采用了“忽略感應(yīng)磁場”這一假設(shè),其應(yīng)用范圍受到限制.為了拓展低磁雷諾數(shù)MHD方法的應(yīng)用范圍,本文針對“忽略感應(yīng)磁場”這一核心假設(shè),探討低磁雷諾數(shù)MHD方法在高超聲速領(lǐng)域適用性問題,分析全MHD方法在高超聲速流動(dòng)領(lǐng)域局限性的產(chǎn)生原理,創(chuàng)建基于感應(yīng)電流積分計(jì)算感應(yīng)磁場對低磁雷諾數(shù)MHD方法進(jìn)行修正的數(shù)值模擬方法.在此基礎(chǔ)上結(jié)合典型計(jì)算狀態(tài),分析“忽略感應(yīng)磁場”對計(jì)算結(jié)果的影響,探討低磁雷諾數(shù)方法的適用性問題.

2 低磁雷諾數(shù)MHD方法適用性分析

采用低磁雷諾數(shù)MHD方法,忽略了感應(yīng)磁場的影響,必須滿足低磁雷諾數(shù)近似條件,即Rem≡μeσULp?1.這就對磁導(dǎo)率μe、電導(dǎo)率σ、速率U以及電磁流動(dòng)相互作用區(qū)域Lp等多個(gè)方面的流體特征提出了要求,客觀上限制了低磁雷諾數(shù)MHD方法的應(yīng)用范圍.

針對低磁雷諾數(shù)MHD方法的研究很多,但對于高超聲速流動(dòng)適用的低磁雷諾數(shù)條件,目前仍存在一些爭議,尚未見到統(tǒng)一的結(jié)論.2002 年,Poggie和Gaitonde[21]開展了鈍體飛行器磁流體力學(xué)控制研究,認(rèn)為Rem=3 時(shí)感應(yīng)磁場可以忽略,可采用低磁雷諾數(shù) MHD 方法;2007 年,Fujino 等[22]開展了再入飛行器磁控技術(shù)研究,認(rèn)為Rem<1 時(shí)可采用低磁雷諾數(shù)MHD方法;Khan等[23]采用定電導(dǎo)率方法開展了磁場對邊界層速度分布影響的研究,探討了低磁雷諾數(shù)MHD方法的適用性,發(fā)現(xiàn)當(dāng)Rem=2.5 時(shí)低磁雷諾數(shù)MHD方法基本失效,并給出了低磁雷諾數(shù)條件Rem<0.125 ,在此基礎(chǔ)上開展了二維鈍錐體Rem=0.179 時(shí)的磁流體數(shù)值模擬,發(fā)現(xiàn)采用低磁雷諾數(shù)MHD方法,對計(jì)算結(jié)果的影響小于 2.02%.2008 年,MacCormack[24]開展了強(qiáng)磁場下的磁流體力學(xué)控制數(shù)值模擬,認(rèn)為Rem<1時(shí)可采用低磁雷諾數(shù)MHD方法,并開展了Rem=0.17 時(shí)磁流體加速裝置研究,發(fā)現(xiàn)采用低磁雷諾數(shù)MHD方法對軸線速度造成的計(jì)算偏差很小;同期,田正雨[1]基于國外的低磁雷諾數(shù)條件的研究,給出了較為保守的低磁雷諾數(shù)條件Rem<0.1;2009 年,Kim[25]對高超聲速流動(dòng)的低磁雷諾數(shù)條件進(jìn)行了探討,認(rèn)為當(dāng)氣體電離度小于1%時(shí)可應(yīng)用低磁雷諾數(shù)MHD方法;2010年,Bocharov[26]綜合了Khan等[23]和MacCormack[24]等在低磁雷諾數(shù)MHD方法適用性方面的研究,給出低磁雷諾數(shù)條件Rem<0.5 ;2013 年,Fujino 和Ishikawa[8]開展了再入返回器大鈍體外形的磁流體數(shù)值模擬,計(jì)算狀態(tài) 64 km 時(shí),磁雷諾數(shù)Rem=29.12 (Lp取鈍體尺寸 0.4 m),感應(yīng)磁場最大僅相當(dāng)于外加磁場的5.4%,采用低磁雷諾數(shù)MHD方法(忽略感應(yīng)磁場)對流場結(jié)構(gòu)和表面熱流的影響很小,其偏差基本可以忽略.

由于低磁雷諾數(shù)條件尚未達(dá)成一致性的結(jié)論,為了保證低磁雷諾數(shù)MHD方法的絕對適用,可選取最為保守的限制條件,即Rem<0.1.但對于高超聲速流動(dòng),由于其流速較高,很容易就接近或超出這一條件.例如,對于特征電導(dǎo)率 500 S/m、特征長度 1 m和特征速度 6 km/s的磁流體狀態(tài),其特征磁雷諾數(shù)高達(dá)3.768.

事實(shí)上,國內(nèi)外有很多研究在具體應(yīng)用低雷諾數(shù)方法時(shí),沒有嚴(yán)格討論其適用性問題.例如,2005年,Otsu[7]在開展鈍體磁流體數(shù)值模擬研究時(shí)采用低磁雷諾數(shù)MHD方法,其磁雷諾數(shù)Rem=3.26 (特征電導(dǎo)率取 400 S/m);2007 年,Yoshino等[27]在開展沿再入軌道的磁流體數(shù)值模擬時(shí)采用低磁雷諾數(shù)MHD方法,飛行高度75 km時(shí)磁雷諾數(shù)為 12.66.2012 年,Nagata 等[6]在開展鈍柱體磁流體力學(xué)控制阻力特性研究時(shí)采用低磁雷諾數(shù)MHD方法,其磁雷諾數(shù)Rem=2.04 (特征電導(dǎo)率取 250 S/m);2015 年,Takahashi等[9]在開展了火星探測器磁流體研究,經(jīng)測算其磁雷諾數(shù)Rem高達(dá)196.9 (特征電導(dǎo)率取波后峰值電導(dǎo)率1600 S/m),數(shù)值模擬時(shí)仍采用低磁雷諾數(shù)MHD方法;2016年,國內(nèi)學(xué)者在開展返回艙OREX磁控?zé)岱雷o(hù)系統(tǒng)時(shí)采用低磁雷諾數(shù)MHD方法[28],經(jīng)本文測算其數(shù)值模擬狀態(tài)的磁雷諾數(shù)Rem高達(dá)11.6.

造成這種現(xiàn)狀的主要原因可能包括兩個(gè)方面:

1)國外在開展低磁雷諾數(shù)條件探討時(shí),一般通過對比低磁雷諾數(shù)MHD方法和全MHD方法.而全MHD方法在高超聲速領(lǐng)域的應(yīng)用受到很大的限制,這將在某種程度上影響數(shù)值模擬的精準(zhǔn)度.

全MHD方程組中的磁擴(kuò)散方程包含磁擴(kuò)散項(xiàng)νe?2B,其磁擴(kuò)散率νe=1/(σμe) ,這里B為磁場磁感應(yīng)強(qiáng)度.由于高超聲速流動(dòng)的氣體磁導(dǎo)率μe、電導(dǎo)率σ一般較低,磁擴(kuò)散率高,其磁擴(kuò)散項(xiàng)大,甚至高出磁擴(kuò)散方程中其他項(xiàng)幾個(gè)數(shù)量級,此時(shí)方程存在巨大的剛性,極大地影響了計(jì)算的穩(wěn)定性和收斂性.尤其是高超聲速流場中很多區(qū)域氣體電導(dǎo)率接近于零,其磁擴(kuò)散項(xiàng)趨近于無窮大,全MHD方程組幾乎無法直接求解,需要進(jìn)行近似處理.

從物理的角度來看,電導(dǎo)率越小,磁場與流動(dòng)的相互作用應(yīng)該越弱,當(dāng)電導(dǎo)率趨近于零時(shí),高超聲速磁流體流動(dòng)就退化為一般的高超聲速流動(dòng),使物理現(xiàn)象和問題得到簡化.此時(shí)全MHD方法反而出現(xiàn)數(shù)值模擬方面的困難,這在一定程度上反映了這種基于磁擴(kuò)散方程的計(jì)算方法的局限性.

這種局限性產(chǎn)生機(jī)制,可結(jié)合全MHD方法的基本原理進(jìn)行分析.由完整的磁流體力學(xué)控制基本方程[1,2]可知,在等離子體電磁流體相互作用過程中,電磁場對流場動(dòng)量輸運(yùn)F和能量輸運(yùn)EM分別為:

這里U為氣流速度.可見,電流J和電場E在電磁流體相互作用過程中,占據(jù)極其重要的地位,直接表征了電磁場對流動(dòng)的干擾.而在全MHD方法中電流J和電場E由磁場的變化得到[4]:

對于數(shù)值模擬來說,物理量偏導(dǎo)數(shù)的數(shù)值模擬精準(zhǔn)度一般要低于物理量本身,例如熱流和摩擦阻力的計(jì)算精準(zhǔn)度,一般就低于溫度和速度等參量的模擬精度,因此?×B的計(jì)算很容易出現(xiàn)數(shù)值誤差.而對于地球大氣,磁導(dǎo)率非常小μe≈4π×10?7h/m.當(dāng)磁場B的一階偏導(dǎo)數(shù)由于某種原因(如數(shù)值誤差、迭代中間差量等)出現(xiàn)較小的非物理擾動(dòng)時(shí),該擾動(dòng)反映在電流J和電場E上就會(huì)被放大106倍甚至更高.尤其是低電導(dǎo)率或零電導(dǎo)率時(shí),電場E就面臨著“極小值除以極小值”甚至“零除以零”的極端情況.結(jié)合(1)式可以看出,這將直接影響電磁流動(dòng)核心的動(dòng)量與能量傳輸過程.可見,全MHD方法“以磁場的變化表征電場、電流變化”基本思路,在低磁導(dǎo)率、低電導(dǎo)率的高超聲速飛行領(lǐng)域存在缺陷,這將極大地影響計(jì)算的穩(wěn)定性和有效性.

此外,由于磁擴(kuò)散方程的數(shù)值解強(qiáng)烈依賴于磁場邊界條件,因此不恰當(dāng)邊界條件可能帶來計(jì)算的準(zhǔn)確性問題.求解經(jīng)典麥克斯韋方程組,分界面上電磁場的連續(xù)性條件可寫為[29]:

這里H2,H1分別為界面兩邊的磁場強(qiáng)度矢量,B2,B1分別為界面兩邊的磁感應(yīng)強(qiáng)度矢量,js為界面電流密度,n為界面方向矢量.由于界面電流很難有效確定,真實(shí)的磁場交界面的切向分量很難直接給出.因此全MHD方法磁場氣固界面一般都采用近似條件[1?3]: 對于完全導(dǎo)體界面 dB/dn=0 ;對于絕緣壁面B=Bw.對比(3)式可知,這兩種邊界均不能完全涵蓋分界面電磁場的連續(xù)性條件.磁場邊界的不準(zhǔn)確,將可能引入非物理的數(shù)值模擬誤差.

2005 年,MacCormack[30]采用全 MHD方法開展鈍體磁流體數(shù)值模擬,得到的結(jié)論是磁場使頭部區(qū)域 熱 流上升;2006 年,Khan 等[31]采用全MHD方法開展鈍體磁流體數(shù)值模擬,其中部分狀態(tài)磁雷諾數(shù)僅為0.01,此時(shí)感應(yīng)(誘導(dǎo))磁場遠(yuǎn)小于外加磁場,但其計(jì)算結(jié)果與低磁雷諾數(shù)MHD方法結(jié)果存在明顯差異,這一計(jì)算結(jié)果的有效性也值得商榷;2009年,國內(nèi)學(xué)者采用全MHD方法開展了RAM-C磁流體數(shù)值模擬[32],得到的結(jié)論同樣是磁場使頭部熱流上升.事實(shí)上,在飛行器頭部加載一定的磁場可以有效降低駐點(diǎn)熱流,這一點(diǎn)在試驗(yàn)中已經(jīng)得到了驗(yàn)證.例如意大利航天中心開展了鈍體磁控?zé)岱雷o(hù)試驗(yàn),磁場可使表面熱流下降46%[13,28].近兩年國防科技大學(xué)柳軍等開展了鈍體的磁控?zé)岱雷o(hù)試驗(yàn)研究,試驗(yàn)結(jié)果也表明外加磁場能有效地降低端頭熱流.

由此可見,在探討低磁雷諾數(shù)條件時(shí),全MHD方法在高超聲速領(lǐng)域的局限性,有可能會(huì)影響其定量結(jié)論的取得.

2)開展低磁雷諾數(shù)條件定量分析[12,23,24,30,31,33]時(shí),氣體電導(dǎo)率計(jì)算多采用定電導(dǎo)率方法.

定電導(dǎo)率方法,即假定全流場的氣體電導(dǎo)率均為某一固定值,此時(shí)等離子體特征電導(dǎo)率及其作用區(qū)域的特征尺度較為確定,因此磁雷諾數(shù)提取相對簡單.而對于高超聲速流動(dòng),其混合氣體電導(dǎo)率與流動(dòng)結(jié)構(gòu)緊密相關(guān).一般來說,僅激波后很小區(qū)域內(nèi)氣體電導(dǎo)率相對較高,達(dá)到102—103S/m量級,其他流場絕大部分區(qū)域電導(dǎo)率均小于10 S/m,甚至接近于0.其等效等離子體特征電導(dǎo)率遠(yuǎn)小于流場中峰值電導(dǎo)率,等效的等離子體特征厚度遠(yuǎn)小于飛行器特征尺度(如球頭半徑等),且這兩者的直接提取相對困難.此時(shí),如果采用流場中峰值電導(dǎo)率和飛行器特征尺度計(jì)算磁雷諾數(shù),就會(huì)使計(jì)算得到的磁雷諾數(shù)遠(yuǎn)遠(yuǎn)大于真實(shí)的等效結(jié)果.

由此可見,低磁雷諾數(shù)MHD方法在低磁雷諾數(shù)條件的定量表征及其應(yīng)用上,仍存在不確定性.對于具有低電導(dǎo)率特征的高超聲速磁流體模擬,選擇偏向保守的低磁雷諾數(shù)條件(Rem<0.1 )來判斷低磁雷諾數(shù)MHD方法是否適用,其必要性仍有待商榷.反之,如果完全忽視低磁雷諾數(shù)條件的限制,則有可能由于感應(yīng)磁場不可忽略,引入較大的計(jì)算偏差.而常用的考慮“感應(yīng)磁場”的全MHD方法在高超聲速領(lǐng)域存在較大的局限性,影響了計(jì)算的有效性.因此,有必要建立適用于高超聲速領(lǐng)域的考慮感應(yīng)磁場的磁流體數(shù)值模擬方法.

3 修正方法的建立

由于低磁雷諾數(shù)MHD方法在流動(dòng)控制方程方面沒作任何簡化,其形式與完整的磁流體力學(xué)方程組[1,2]中的流動(dòng)方程一致.因此如果在低磁雷諾數(shù)MHD方法中,進(jìn)一步考慮感應(yīng)磁場的影響,那么這種方法就能從理論上突破低磁雷諾數(shù)條件的限制,在更大程度上反映完整磁流體力學(xué)方程的物理本質(zhì),因此適用范圍更加廣泛.

按照磁場矢量疊加原理,總的物理磁場可寫為

這里Bext為外加磁場,由機(jī)載磁場發(fā)生裝置產(chǎn)生;Bin為感應(yīng)磁場(或稱誘導(dǎo)磁場),由等離子體中感應(yīng)電流誘導(dǎo)產(chǎn)生.對于高超聲速磁流體力學(xué)控制來說,磁場對流場動(dòng)量和能量的輸運(yùn)需要時(shí)間的累積,因此其控制過程一般可近似為定常(或準(zhǔn)定常)狀態(tài).此時(shí)其收斂狀態(tài)的流場中感應(yīng)電流J和感應(yīng)磁場Bin的分布不隨(或近似不隨)時(shí)間變化,滿足(或近似滿足)電流穩(wěn)恒和靜磁條件.因此r處的感應(yīng)磁場Bin可由畢奧-薩伐爾定律積分給出:

這里為r′和V′為體積積分的位置矢量和區(qū)間.

由畢奧-薩伐爾定律可直接推導(dǎo)出安培環(huán)路定理和磁場高斯定律[29],因此基于感應(yīng)電流磁誘導(dǎo)積分得到的Bin,自動(dòng)滿足了安培環(huán)路定理和磁場高斯定律.此時(shí),按照磁場唯一性原理[29],這一Bin應(yīng)與直接求解麥克斯韋方程中安培環(huán)路定理和磁場高斯定律等同,這從理論上說明了該方法的可行性.

基于畢奧-薩伐爾定律積分計(jì)算Bin,在省去了繁瑣的磁擴(kuò)散方程或磁矢量勢泊松方程組離散與迭代求解過程的同時(shí),還可以考慮邊界處可能存在的面電流影響(可將其直接列入積分),避免磁場邊界或磁矢量勢邊界難以確定的問題.因此其物理意義更加明確,計(jì)算更加簡便.

但這種方法也存在著明顯的缺陷.在感應(yīng)磁場的每一步耦合迭代更新過程,對于流場中每一個(gè)流體微元都需要進(jìn)行全場的感應(yīng)電流的磁誘導(dǎo)積分.假設(shè)流場中N個(gè)微元,那么每一步迭代耦合的計(jì)算量將會(huì)達(dá)到N2量級.隨網(wǎng)格量增大,計(jì)算量將迅速增加.對于高超聲速復(fù)雜流動(dòng)問題,計(jì)算網(wǎng)格很容易達(dá)到千萬(107)量級,此時(shí)計(jì)算量將高達(dá)1014量級,一般的計(jì)算機(jī)系統(tǒng)很難承受.

為了減小計(jì)算量和計(jì)算時(shí)間,必須對這種基于畢奧-薩伐爾定律直接積分的方法進(jìn)行改進(jìn).這里由磁場高斯定律?·Bin=0 ,引入磁矢量勢Ain的概念:

然后由畢奧-薩伐爾定律給出磁矢量勢,(5)式變形可得:

這里算符?是對r的微分算符,與r′無關(guān),有恒等式:

因此(7)式可寫為

這里磁矢量勢Ain可寫為

可見,只要基于感應(yīng)電流積分計(jì)算式(10)式得到磁矢量勢,然后就可由(6)式得到感應(yīng)磁場Bin分布.相比于(5)式,(10)式的積分計(jì)算要簡單得多,計(jì)算量也要小得多.

由于高超聲速磁流體力學(xué)控制磁場與流動(dòng)相互作用的主要區(qū)域,相比于全流場空間的占比較小.因此,在計(jì)算感應(yīng)磁場Bin時(shí),可以根據(jù)磁流體流動(dòng)特征對計(jì)算范圍進(jìn)行縮減,進(jìn)一步減小計(jì)算量.

1)縮減需要計(jì)算感應(yīng)磁場Bin的區(qū)域

流場中并不是所有區(qū)域都需要計(jì)算感應(yīng)磁場Bin,圖1(a)給出了文獻(xiàn)[34]計(jì)算的流場電導(dǎo)率,可以看出很大區(qū)域的電導(dǎo)率趨近于零.當(dāng)流場某區(qū)域氣流的電導(dǎo)率很低甚至接近于零時(shí),流體幾乎不受常規(guī)磁場的宏觀作用,此時(shí)該區(qū)域感應(yīng)磁場存在與否,對流動(dòng)干擾很小.換句話說,如果某區(qū)域的流體受常規(guī)磁場的影響很小,那么該區(qū)域的感應(yīng)磁場可以不用計(jì)算.由于磁相互作用數(shù)可以在一定程度上表征磁場與流動(dòng)耦合作用的整體效果,這里引入單位磁相互作用數(shù)Q′m的概念,令磁相互作用數(shù)中的特征磁感應(yīng)強(qiáng)度和特征長度分別為1 T和1 m,即有=σ/ρu,這里ρ,u為氣體密度和速度大小.因此,需要計(jì)算感應(yīng)磁場Bin的區(qū)域,可以通過氣體電導(dǎo)率或者單位磁相互作用數(shù)高于某個(gè)值進(jìn)行判定:

為了盡可能地覆蓋明顯影響磁流體力學(xué)控制的絕大部分區(qū)域,同時(shí)盡可能減少需要積分的計(jì)算區(qū)域,截?cái)嚯妼?dǎo)率σc或者截?cái)鄦挝淮畔嗷プ饔脭?shù)Qc應(yīng)取值恰當(dāng).例如對于高超聲速流動(dòng),其激波后電導(dǎo)率為 102—103S/m,σc可取為 0.1—1.0 S/m或者更小.為了減小不同計(jì)算狀態(tài)的人為干預(yù),這里引入截?cái)嘁蜃应?

其中σ0為磁流體作用區(qū)域典型電導(dǎo)率,為磁流體作用區(qū)域典型單位磁相互作用數(shù).在迭代初期,截?cái)嘁蜃应强蛇x用較大的值,如η=10?1—10?3,以減小計(jì)算量,加快收斂;在迭代后期,η可選用較小的值,如η=10?4—10?6,以提高計(jì)算精度.

圖1 流場中電導(dǎo)率和環(huán)形感應(yīng)電流分布[34] (a)電導(dǎo)率;(b) 電流密度Fig.1.Distribution of electronic conductivity and annular electric current: (a) Conductivity;(b) current.

2)縮減計(jì)算某一點(diǎn)感應(yīng)磁場Bin所需要的感應(yīng)電流積分區(qū)域

流場中某一點(diǎn)的感應(yīng)磁場,由全場每個(gè)微元電流在該點(diǎn)產(chǎn)生的感應(yīng)磁場疊加得到.由圖1(b)可以看出,高超聲速流場中感應(yīng)電流主要集中于流場與電磁場相互作用的小部分區(qū)域,很大區(qū)域電流趨近于零.因此,也可以引入類似于(11)式和(12)式的條件,來判別計(jì)算感應(yīng)磁場Bin時(shí)需要積分的空間區(qū)域:

這里Jc為截?cái)嚯娏髅芏却笮?J0為磁相互作用區(qū)典型電流密度.

3)縮減相互作用區(qū)間

此外由(5)式還可以看出,由電流微元產(chǎn)生的感應(yīng)磁場與距離的平方呈反比關(guān)系.可見由穩(wěn)恒(或近似穩(wěn)恒)的感應(yīng)電流產(chǎn)生的磁場,其主要作用區(qū)域應(yīng)該為有限的空間,將其定為

這里L(fēng)c為磁流體作用區(qū)域的特征尺度;L為感應(yīng)磁場作用區(qū)間尺度,超出L流場空間,不參與該外加磁場作用區(qū)域的感應(yīng)磁場計(jì)算;λ為放大因子,可結(jié)合具體的物理實(shí)際進(jìn)行設(shè)置,本文數(shù)值模擬時(shí)λ=5.

4)耦合策略

由于流動(dòng)收斂的速度遠(yuǎn)遠(yuǎn)小于感應(yīng)電流生成誘導(dǎo)磁場的響應(yīng)速度,因此在實(shí)際的數(shù)值模擬過程中,先采用不考慮感應(yīng)磁場的低磁雷諾數(shù)MHD方法迭代計(jì)算直至結(jié)果基本收斂,然后每隔500—5000步進(jìn)行一次感應(yīng)磁場修正,迭代計(jì)算直至結(jié)果重新收斂.

4 數(shù)值計(jì)算分析

4.1 無限長直流導(dǎo)線產(chǎn)生磁場

為了保證電流數(shù)值離散積分過程的正確性,首先對電流積分計(jì)算磁場的程序模塊進(jìn)行校驗(yàn),這是本文低磁雷諾數(shù)MHD方法修正的基礎(chǔ).

對于通電電流為I的無限長直導(dǎo)線,導(dǎo)線外距離為r的任意一點(diǎn)的磁場B由安培環(huán)路定理計(jì)算:

這里L(fēng)為半徑為r的圓周,S為圓面,圓面方向?yàn)殡娏鞣较?由于對稱性因此有磁感應(yīng)強(qiáng)度大小B=μeI/2πr,其方向符合右手系.

采用半徑為1 mm的直導(dǎo)線,電流I=1.0A.假設(shè)電流在導(dǎo)線截面內(nèi)均勻分布,則電流密度約為J=3.18×105A/m2,這與圖1(b)流場感應(yīng)電流密度在同一數(shù)量級.為了保證“無限長”的假設(shè)近似成立,本文取直導(dǎo)線長為100 m.采用圓柱形網(wǎng)格,網(wǎng)格點(diǎn)設(shè)為 1 0001×51×181 ,沿導(dǎo)線長方向的網(wǎng)格點(diǎn)均勻分布,其他兩個(gè)方向非均勻分布.

圖2給出了不同位置的磁場感應(yīng)強(qiáng)度與理論值的比較.由圖2(a)可以看出,在導(dǎo)線的端頭(x=0.0 m)附近,數(shù)值積分的磁感應(yīng)強(qiáng)度明顯小于理論值,這是因?yàn)閤<0.0m 的部分被人為截?cái)?沒有參與數(shù)值積分,因此不滿足“近似無限長條件”.隨著x增大,這種“人為截?cái)唷碑a(chǎn)生影響逐漸減小,數(shù)值結(jié)果逐漸趨近于理論值,長直導(dǎo)線中點(diǎn) (x=50.0 m),數(shù)值結(jié)果幾乎與理論值完全重合.圖2(b)中也有類似的現(xiàn)象: 當(dāng)r較小時(shí),數(shù)值結(jié)果與理論值幾乎完全重合,隨r進(jìn)一步增大,r逐漸接近甚至大于導(dǎo)線的長度,此時(shí)“人為截?cái)唷碑a(chǎn)生的影響逐漸增大,不再滿足“近似無限長條件”,數(shù)值結(jié)果與理論值偏離.

由此可見,數(shù)值積分計(jì)算結(jié)果符合理論預(yù)期.這說明本文基于電流積分的磁場計(jì)算程序模塊能較為準(zhǔn)確地給出空間磁場分布,滿足感應(yīng)磁場計(jì)算精度要求.

4.2 三維鈍柱體外形考慮感應(yīng)磁場修正的磁流體數(shù)值模擬

為了進(jìn)一步校驗(yàn)本文修正計(jì)算方法的有效性,開展鈍柱體外形高超聲速磁流體數(shù)值模擬,該算例有低磁雷諾數(shù)方法和全MHD方法計(jì)算結(jié)果[33].計(jì)算外形頭部半徑 0.01 m,柱體長 0.075 m.計(jì)算來流為混合電離氣體,來流溫度 9000.0 K,來流壓力7726.0 Pa,來 流 密 度 0.01429 kg/m3,來 流 速 度6710.0 m/s,壁面溫度 12000.0 K.采用偶極子磁場,磁場中心位于頭部球心,磁特征感應(yīng)強(qiáng)度B0=4.0T,磁場特征長度為 0.01 m.采用固定電導(dǎo)率方法,全場混合氣體電導(dǎo)率均為 4800 S/m,特征磁相互作用數(shù)為69.

圖3和圖4分別給出了本文計(jì)算的外加磁場和感應(yīng)磁場分布,并與文獻(xiàn)[33]給出的結(jié)果進(jìn)行對比,圖中數(shù)值乘以轉(zhuǎn)化為以特斯拉為單位的量.由圖可以看出,本文修正方法計(jì)算得到的感應(yīng)磁場,相對于外加磁場量值較小;除局部細(xì)節(jié)外,整體分布變化規(guī)律和量值均與全MHD方法的結(jié)果[33]相符合,這說明本文修正方法能較好計(jì)算磁流體流動(dòng)過程中的感應(yīng)磁場分布.

為了進(jìn)一步分析感應(yīng)磁場細(xì)節(jié)差異產(chǎn)生的原因,圖5給出了不同計(jì)算方法得到的流場壓力分布.可以看出,在磁場作用下本文計(jì)算得到的流場壓力分布規(guī)律與文獻(xiàn)[33]基本一致,但本文給出的激波脫體距離明顯大于文獻(xiàn)結(jié)果.這是由于文獻(xiàn)[33]沒有給出混合電離氣體的成分,也沒有明確氣體狀態(tài)計(jì)算方法,而本文采用的等效比熱比的氣體計(jì)算方法,可能與文獻(xiàn)[33]中的氣體模型不一致.事實(shí)上,文獻(xiàn)[33]給出的該狀態(tài)磁場使激波脫體距離增大的倍數(shù)并不確定,為 3—7.5倍[33],圖2—圖5中文獻(xiàn)結(jié)果為增大3倍的結(jié)果,而本文計(jì)算結(jié)果為磁場使激波脫體距離增大約4—5倍.本文計(jì)算激波脫體距離較遠(yuǎn),導(dǎo)致了感應(yīng)磁場分布細(xì)節(jié)上的差異.

圖4 感應(yīng)磁場 (a)本文 BX;(b)文獻(xiàn) BX;(c)本文 BY;(d)文獻(xiàn) BYFig.4.Induced magnetic field: (a) BX of this paper;(b) BX[33];(c) BY of this paper;(d) BY[33].

由圖5還可以看出,本文修正方法計(jì)算得到的壓力與一般的低磁雷諾數(shù)MHD方法模擬結(jié)果,除局部細(xì)節(jié)外的分布規(guī)律基本一致,與文獻(xiàn)[33]“低磁雷諾數(shù)MHD方法和全MHD方法結(jié)果之間的差異”類似.這里以磁場尺度和來流速度為特征量,計(jì)算磁雷諾數(shù)Rem≈0.4.可見,此時(shí)采用低磁雷諾數(shù)假設(shè)忽略感應(yīng)磁場,對流場結(jié)構(gòu)的影響不大,這與圖3和圖4顯示的“感應(yīng)磁場相對于外加磁場的量值較小”相互印證.

4.3 Ram-C鈍錐“忽略感應(yīng)磁場”對磁流體數(shù)值模擬影響的分析

文獻(xiàn)[35]為了分析電導(dǎo)率模擬準(zhǔn)確性問題,基于低磁雷諾數(shù)方法探討了國內(nèi)外常見的多種電導(dǎo)率模型差異及其影響,其中采用電導(dǎo)率模型M6 時(shí),流場中峰值電導(dǎo)率高達(dá) 6000 S/m,其磁雷諾數(shù)遠(yuǎn)高于0.1,低磁雷諾數(shù)方法的適用性有待商榷.本文結(jié)合這一典型飛行狀態(tài)與計(jì)算條件,開展一般低磁雷諾數(shù)方法和本文修正方法的數(shù)值對比計(jì)算,分析“忽略感應(yīng)磁場”對磁流體控制數(shù)值模擬的影響,探討低磁雷諾數(shù)條件適用性.

采用RAM-C鈍錐體外形[35,36],計(jì)算飛行高度 71 km,飛行速度 7650 m/s,等溫壁面條件,壁面溫度 1500 K,完全非催化壁面條件.磁場配置采用磁偶極子磁場,磁感應(yīng)特征強(qiáng)度B0=0.5T ,特征長度r0=0.1524 m,磁偶極子方向?yàn)橹苯亲鴺?biāo)橫軸負(fù)方向.采用電導(dǎo)率模型M6[35],其形式為σ=1.56×10?2T1.5/ln(1.23×107T1.5/n0e.5),T為溫度,ne為電子數(shù)密度.采用11組分Park反應(yīng)模型計(jì)算等離子體分布.

圖6給出了采用M6模型計(jì)算得到的流場中電導(dǎo)率分布,圖6(a)為采用一般低磁雷諾數(shù)MHD方法模擬得到的流場電導(dǎo)率分布云圖,圖6(b)為采用修正方法和一般低磁雷諾數(shù)MHD方法計(jì)算得到的駐點(diǎn)線電導(dǎo)率分布比較.可以看出,在磁場作用下流場峰值電導(dǎo)率可達(dá)5000 S/m,但取得峰值電導(dǎo)率的等離子區(qū)域厚度較薄,波后較大區(qū)域內(nèi)氣體電導(dǎo)率僅為1000 S/m左右.采用不同方法計(jì)算,磁場對激波脫體距離的外推效果存在差別,采用修正方法得到的激波脫體距離,比采用一般低磁雷諾數(shù)MHD方法結(jié)果稍小.

圖5 不同方法計(jì)算的流場壓力分布 (a)文獻(xiàn)低磁雷諾數(shù)MHD方法;(b)文獻(xiàn)全MHD方法;(c)本文低磁雷諾數(shù)MHD方法;(b)本文修正方法Fig.5.Distribution of pressure in the flow computed by different method: (a) Low Rem method[33];(b) full MHD method[33];(c) low Rem method of this paper;(d)improved method of this paper.

圖6 采用電導(dǎo)率模型M6計(jì)算的流場電導(dǎo)率分布 (a)全場云圖;(b)駐點(diǎn)線參數(shù)分布Fig.6.Distribution of electronic conductivity using M6: (a) Full contour map;(b) parameters along stagnation line.

圖7給出了修正方法計(jì)算得到的感應(yīng)磁場分布.對比圖4可以看出,感應(yīng)磁場的整體分布與4.2節(jié)鈍柱體的感應(yīng)磁場分布大體相似,但二者又有明顯的差異.造成差異的原因可能有兩個(gè)方面:一是計(jì)算外形差別;二是電導(dǎo)率分布存在差別,尤其是激波波前位置,4.2節(jié)采用定電導(dǎo)率方法,其流場波前氣體電導(dǎo)率為4800 S/m,此時(shí)波前流動(dòng)也會(huì)產(chǎn)生較強(qiáng)的環(huán)形感應(yīng)電流;而本節(jié)狀態(tài)的流場波前電導(dǎo)率接近于零,波前環(huán)形電流也接近于零.這種感應(yīng)電流分布的差異必然影響感應(yīng)磁場的分布.

為了衡量感應(yīng)磁場的相對大小,圖8進(jìn)一步給出了外加磁場分布和全磁場分布,這里全磁場是外加磁場與感應(yīng)磁場疊加的結(jié)果.結(jié)合圖7可以看出,除局部細(xì)節(jié)外,感應(yīng)磁場方向大體上與外加磁場相反,感應(yīng)磁場的影響在某種程度上相當(dāng)于對外加磁場作用效果的削弱.這是符合磁電阻擴(kuò)散基本原理的,由楞次定律可知,具有導(dǎo)電性的等離子體在磁場中運(yùn)動(dòng)時(shí),其感應(yīng)磁場方向總是對抗磁通量的變化方向,因此整體上表現(xiàn)為對原磁場的某種削弱,整體影響幅度不大,其x分量最大值相當(dāng)于外加磁場x方向最大值的7.9%左右,其y分量最大值相當(dāng)于外加磁場y方向最大值的4.6%左右.

圖7 修正方法計(jì)算的鈍錐 RAM-C 感應(yīng)磁場 (a) Bx 分量;(b) By 分量Fig.7.Induced magnetic field of RAM-C using improved method: (a) Component Bx ;(b) component By.

圖8 鈍錐 RAM-C 外加磁場和修正方法計(jì)算的全磁場分布 (a)全磁場 Bx 分量;(b)全磁場 By 分量;(c)外加磁場 Bx 分量;(b) 外加磁場 By 分量Fig.8.Total magnetic field computed using improved method and externally applied magnetic field of RAM-C: (a) Total Bx ;(b) total By ;(c) externally Bx ;(c) externally By.

表1給出了不同條件或方法計(jì)算得到的阻力系數(shù),可以看出,采用修正方法計(jì)算得到的阻力系數(shù)與一般低磁雷諾數(shù)MHD方法計(jì)算結(jié)果的差異在1%左右,其磁控效果沒有本質(zhì)差別.

表1 鈍錐 RAM-C 阻力系數(shù)Table 1.Drag coefficient of RAM-C.

圖9進(jìn)一步給出了不同計(jì)算方法得到的表面熱流分布.由圖可以看出,采用修正方法計(jì)算得到的熱流與采用一般低磁雷諾數(shù)MHD方法的結(jié)果局部存在差異,但整體分布趨于一致,這說明磁控?zé)岱雷o(hù)效果沒有本質(zhì)差異.在X=0 m 附近,兩者差異稍大,采用修正方法的熱流的磁控?zé)岱雷o(hù)效率比一般低磁雷諾數(shù)MHD方法的下降不到10%.可見,考慮感應(yīng)磁場影響在一定程度上降低了局部區(qū)域的磁控?zé)岱雷o(hù)效果.

圖9 修正方法和一般低磁雷諾數(shù)MHD方法計(jì)算得到的熱流分布比較Fig.9.Heat flux computed using Low Rem method or improvbed method.

這里結(jié)合磁雷諾數(shù)進(jìn)行分析,按傳統(tǒng)的低磁雷諾數(shù)的計(jì)算方法,以流場中峰值電導(dǎo)率、磁場特征尺度和氣體來流流速為特征量,則磁雷諾數(shù)Rem≈7.3,此時(shí)不符合低磁雷諾數(shù)條件要求.如果以波后較大區(qū)域電導(dǎo)率1000 S/m和頭部等離子體層厚度0.1 m為特征量,則磁雷諾數(shù)Rem≈0.96.兩者差異較大,后者綜合考慮了等離子體的實(shí)際分布特征、磁擴(kuò)散方向以及電磁相互作用區(qū)間,因而更加準(zhǔn)確,但其結(jié)果仍然比保守的低磁雷諾數(shù)條件(Rem?0.1 )高得多.結(jié)合圖8和表1 “忽略感應(yīng)磁場對計(jì)算結(jié)果影響較小”的結(jié)論可以看出,對于本文這種不含人工電離的高超聲速空氣流場計(jì)算狀態(tài),“Rem<0.1 ”的低磁雷諾數(shù)條件過于保守,可擴(kuò)展為Rem<1.0 ,同時(shí)其特征電導(dǎo)率和特征尺度應(yīng)綜合考慮實(shí)際的等離子體分布.

圖10進(jìn)一步給出兩種修正方法計(jì)算得到的熱流和殘差收斂曲線.Modified Method 1 為直接積分計(jì)算畢奧-薩伐爾定律積分的方法,Modified Method 2為基于電流積分計(jì)算磁矢量勢同時(shí)考慮截?cái)嘁蜃訉τ?jì)算區(qū)域進(jìn)行縮減的方法.由圖可以看出,這兩種方法計(jì)算得到的熱流幾乎完全重合,計(jì)算精度相當(dāng).本文 Modified Method 2 的計(jì)算時(shí)間較短,效率較高.

圖10 采用不同修正方法計(jì)算得到的熱流和殘差收斂曲線 (a)熱流;(b)殘差Fig.10.Heat flux and residual error computed using different modified method: (a) Heat flux;(b) residual error.

5 結(jié) 論

(1)本文在詳細(xì)探討低磁雷諾數(shù)方法適用性以及全MHD方法局限性及其原理的基礎(chǔ)上,基于感應(yīng)電流積分計(jì)算磁矢量勢得到感應(yīng)磁場,提出了一種低磁雷諾數(shù)MHD修正計(jì)算方法,突破了“低磁雷諾數(shù)假設(shè)”的限制,增強(qiáng)了低磁雷諾數(shù)MHD方法的適用性,拓寬了應(yīng)用范圍.考核驗(yàn)算表明,該方法可較為準(zhǔn)確地計(jì)算模擬磁流體流動(dòng)中的感應(yīng)磁場,數(shù)值模擬結(jié)果與理論分析及文獻(xiàn)結(jié)果相符,具有較高的可信度.

(2)針對典型計(jì)算狀態(tài),開展了常規(guī)低磁雷諾數(shù)MHD方法和修正方法的數(shù)值對比計(jì)算,探討了感應(yīng)磁場對磁流體力學(xué)控制數(shù)值模擬的影響,分析了采用常規(guī)低磁雷諾數(shù)MHD方法在流場特性、氣動(dòng)力/熱特性等方面造成的計(jì)算偏差.研究結(jié)果表明: 1)高超聲速流動(dòng)過程中感應(yīng)磁場的影響在宏觀上表現(xiàn)為對外加磁場的削弱和扭曲,一定程度上降低了磁控效果;2)文獻(xiàn)[35]采用常規(guī)的低磁雷諾數(shù)MHD方法“忽略感應(yīng)磁場”造成的計(jì)算偏差,對其計(jì)算結(jié)論的取得不構(gòu)成實(shí)質(zhì)性影響;3)相比于全域直接積分計(jì)算畢奧-薩劃爾定律,本文基于電流積分計(jì)算磁矢量勢同時(shí)考慮截?cái)嘁蜃訉τ?jì)算區(qū)域進(jìn)行縮減的耦合計(jì)算方法,在計(jì)算時(shí)間和效率方面具有明顯優(yōu)勢;4)本文計(jì)算條件下,“Rem<0.1”的低磁雷諾數(shù)條件過于保守,建議取為Rem<1.0,同時(shí)其特征電導(dǎo)率和特征尺度應(yīng)綜合考慮實(shí)際的等離子體分布.