李開(kāi) 柳軍 劉偉強(qiáng)

(國(guó)防科學(xué)技術(shù)大學(xué)航天科學(xué)與工程學(xué)院,長(zhǎng)沙410073)

高超聲速飛行器磁控?zé)岱雷o(hù)霍爾電場(chǎng)數(shù)值方法研究?

李開(kāi)?柳軍 劉偉強(qiáng)

(國(guó)防科學(xué)技術(shù)大學(xué)航天科學(xué)與工程學(xué)院,長(zhǎng)沙410073)

(2016年9月16日收到;2017年1月22日收到修改稿)

作為一種新概念高超聲速熱防護(hù)手段,磁控?zé)岱雷o(hù)系統(tǒng)在實(shí)際應(yīng)用中往往需要考慮霍爾效應(yīng)的影響.為了在真實(shí)氣體環(huán)境下求解霍爾電場(chǎng),采用交替隱式近似因子分解法建立并驗(yàn)證了熱化學(xué)非平衡流體域電場(chǎng)數(shù)值求解方法.分析了電場(chǎng)虛擬步進(jìn)因子和收斂性的關(guān)系以及影響步進(jìn)因子取值的因素,提出了當(dāng)?shù)刈儾竭M(jìn)因子加速電場(chǎng)收斂方法.研究表明,存在一個(gè)最優(yōu)的步進(jìn)因子ap使得霍爾電場(chǎng)收斂速度最快,并且隨網(wǎng)格尺度的減小和霍爾系數(shù)的增加,最優(yōu)步進(jìn)因子ap變大,電勢(shì)場(chǎng)收斂速率變慢.對(duì)于局部加密網(wǎng)格而言,當(dāng)?shù)刈儾竭M(jìn)因子法的電勢(shì)收斂性明顯優(yōu)于常規(guī)的定步進(jìn)因子法.

磁流體控制,熱防護(hù),霍爾效應(yīng),泊松方程

1 引言

熱防護(hù)系統(tǒng)是保護(hù)高超聲速飛行安全的關(guān)鍵子系統(tǒng),一直以來(lái)都是高超聲速飛行器系統(tǒng)設(shè)計(jì)的關(guān)鍵一環(huán)[1].由于電磁流動(dòng)控制技術(shù)近年來(lái)得到了飛速發(fā)展[2?4],作為其在高超聲速熱防護(hù)領(lǐng)域的一個(gè)應(yīng)用,磁控?zé)岱雷o(hù)技術(shù)顯示出了巨大的潛在應(yīng)用價(jià)值,得到了各國(guó)科研機(jī)構(gòu)的廣泛關(guān)注[5?7].

對(duì)于高超聲速飛行器前緣磁控?zé)岱雷o(hù)而言,其核心在于在高超聲速飛行器鼻錐內(nèi)部添加磁鐵,通過(guò)外加磁場(chǎng)和波后帶電流體的相互作用,感應(yīng)產(chǎn)生洛倫茲力推出激波、減速流體、偏轉(zhuǎn)流線,從而降低到達(dá)壁面的熱流,實(shí)現(xiàn)熱防護(hù).由于存在霍爾效應(yīng),磁場(chǎng)對(duì)電子的作用會(huì)使流動(dòng)的法向和流向存在電勢(shì)差,即感應(yīng)霍爾電場(chǎng).根據(jù)通用歐姆定律,感應(yīng)電場(chǎng)的存在會(huì)使得電流出現(xiàn)三維效應(yīng),進(jìn)而改變洛倫茲力,影響磁控?zé)岱雷o(hù)系統(tǒng)的效率.為分析霍爾效應(yīng)對(duì)磁控?zé)岱雷o(hù)系統(tǒng)的影響,必須耦合求解熱化學(xué)非平衡流場(chǎng)和霍爾電場(chǎng).高超聲速條件下,由于波后電離度較低,低磁雷諾數(shù)假設(shè)往往可以得到滿足[8].此時(shí),電磁場(chǎng)與非平衡流場(chǎng)的耦合計(jì)算可以通過(guò)求解電勢(shì)泊松方程和包含電磁源項(xiàng)的非平衡流動(dòng)方程實(shí)現(xiàn)[9].流體域霍爾電場(chǎng)求解的關(guān)鍵在于求解電勢(shì)泊松方程.由于作為系數(shù)的電導(dǎo)率在激波位置存在間斷,并且在大霍爾系數(shù)條件下系數(shù)矩陣剛性問(wèn)題嚴(yán)重[10],高超聲速流體域電勢(shì)泊松方程的求解存在魯棒性差、收斂性差等問(wèn)題.Gaitonde等[11]和Wan等[12]分別基于有限差分法和有限體積法,采用虛擬時(shí)間步長(zhǎng)以及近似因子分解-交替隱式方法實(shí)現(xiàn)了三維電勢(shì)泊松方程的數(shù)值離散與求解.Bisek[13]基于有限體積法采用連續(xù)超松弛算法(successive over-relaxation,SOR)完成了三維電勢(shì)泊松方程的迭代求解.呂浩宇等[14]、彭武等[15]同樣基于有限差分法采用超松弛算法求解了二維電勢(shì)泊松方程.Fujino等[16]采用基于Galerkin有限元素法求解二階半離散電勢(shì)泊松方程,與實(shí)驗(yàn)結(jié)果符合良好.胡海洋等[10]采用無(wú)虛擬時(shí)間步的LUSGS預(yù)處理BI-CGSTAB算法解決了大霍爾系數(shù)下泊松方程病態(tài)矩的求解效率低的問(wèn)題.綜上,目前見(jiàn)于文獻(xiàn)的霍爾電場(chǎng)求解方法主要有三種:1)超松弛迭代顯式算法;2)近似因子分解-交替隱式方法;3)有限元方法,如BI-CGSTAB,Galerkin方法等.相比而言,SOR算法為顯式算法,應(yīng)用簡(jiǎn)單且利于并行計(jì)算,但收斂效率低于隱式格式.有限元方法收斂性好,但在處理復(fù)雜壁面條件時(shí)存在困難.

目前關(guān)于霍爾電場(chǎng)的研究往往僅限于Poisson方程求解方法的構(gòu)建,對(duì)霍爾電場(chǎng)收斂性的研究較少.針對(duì)霍爾電場(chǎng)收斂性和步進(jìn)因子ap的關(guān)系,Wan等[12]給出了其建議取值范圍0.001—0.01;Gaitonde[11]則采用了Holst[17]提出的步進(jìn)因子根據(jù)其上下限隨時(shí)間步過(guò)渡變化的指數(shù)公式.但是,Wan給出的范圍只針對(duì)其所采用的網(wǎng)格,不具有普適性.Holst的方法也必須首先確定步進(jìn)因子的上下限,實(shí)際操作困難.因此,本文首先構(gòu)造三維霍爾電場(chǎng)的隱式數(shù)值方法,而后進(jìn)行霍爾電場(chǎng)收斂性的影響因素研究,以便為高超聲速飛行器前緣磁控?zé)岱雷o(hù)系統(tǒng)的電磁場(chǎng)流場(chǎng)耦合計(jì)算分析提供參考.

2 數(shù)學(xué)模型

2.1 控制方程

高超聲速飛行器定常繞流條件下,由于波后等離子體電離度往往低于1×10?4[18,19],此時(shí),低磁雷諾數(shù)假設(shè)通常可以得到滿足,并且外加穩(wěn)定磁場(chǎng)時(shí),磁場(chǎng)和感應(yīng)電場(chǎng)的變化通常都可以忽略不計(jì).因此,可以對(duì)安培-麥克斯韋公式、法拉第定律進(jìn)行簡(jiǎn)化,將對(duì)感應(yīng)場(chǎng)強(qiáng)矢量E以及感應(yīng)電流密度J的求解轉(zhuǎn)化為對(duì)標(biāo)量電勢(shì)?的求解,結(jié)合廣義歐姆定律

可以得到關(guān)于?的電勢(shì)泊松方程

其中u為流體速度矢量,B為磁感應(yīng)強(qiáng)度矢量,電導(dǎo)率張量?σ可以由下式求出:

其中,β為霍爾系數(shù),并且為分析霍爾系數(shù)對(duì)步進(jìn)因子取值、電場(chǎng)計(jì)算收斂性的影響,采用均布霍爾系數(shù)模型.σ為標(biāo)量電導(dǎo)率,采用Fujino電導(dǎo)率模型計(jì)算[16]:

其中,e,ne,Me分別為電子電量、電子數(shù)密度、電子質(zhì)量.為電子和其他組元的有效動(dòng)量傳輸碰撞頻率[20].

2.2 數(shù)值方法

由于電導(dǎo)率σ在激波處存在間斷,作為系數(shù)的電導(dǎo)率張量在激波處也同樣存在間斷,這使得上述電勢(shì)泊松方程(2)在求解時(shí)穩(wěn)定性差.特別是在大霍爾系數(shù)條件下,系數(shù)矩陣很有可能會(huì)出現(xiàn)對(duì)角不占優(yōu)的情況,使得方程剛性嚴(yán)重,影響迭代計(jì)算的收斂性.為此,我們采用能有效抑制剛性的近似因子分解-交替隱式方法.并且為了保證三維多分區(qū)網(wǎng)格在數(shù)值離散中的幾何守恒性,將基于單元中心有限體積法對(duì)上述泊松方程進(jìn)行離散.過(guò)程如下.

對(duì)于特定體積V,(1)式可以化為

將(6)式中右側(cè)展開(kāi),可得到

其中|S|和(n1,n2,n3)分別是單元面的面積及法向矢量.令

(i=1,2,3).則有

對(duì)??在三個(gè)計(jì)算坐標(biāo)方向做分解并應(yīng)用鏈?zhǔn)椒▌t,則有

其中,ξ,η,ζ分別是i,j,k三方向的坐標(biāo)變換系數(shù).令XIC,ETC,ZTC如(10)式所示:

采用虛擬時(shí)間推進(jìn)法進(jìn)行迭代求解,則(6)式可以化為

應(yīng)用中心差分對(duì)(6)式進(jìn)行離散,得到(i,j,k)及其周?chē)它c(diǎn)的差分表達(dá)式.對(duì)于?i,j,k及臨近的六個(gè)點(diǎn)?i+1,j,k,?i?1,j,k,?i,j+1,k,?i,j?1,k,?i,j,k+1,?i,j,k?1進(jìn)行隱式處理?(n+1)=?(n)+??(n),并令

a2j?1,b2j,c2j+1,a3k?1,b3k,c3k+1形式類(lèi)似,不再贅述.令步進(jìn)因子ap=1/?t,則(11)式可以化為

令(13)式中右端項(xiàng)為RE,左端項(xiàng)采用交替隱式-近似因子分解法(ADI-AF),并定義算子Ai,Bj,Ck如(14)式所示:

則最終的電勢(shì)方程迭代計(jì)算式為

霍爾電場(chǎng)邊界條件為:絕緣壁面,J·n=0;導(dǎo)電壁面,?=0;其余邊界為NeuMann邊界,??·n=0.

3 數(shù)值驗(yàn)證

驗(yàn)證算例1取自參考文獻(xiàn)[13],其控制方程為??=xez.該方程擁有解析解?=xez.若令則其控制方程可以化為和(2)式同樣的簡(jiǎn)化形式.給定Dirichlet邊界條件:?=xez.圖1給出了兩套計(jì)算網(wǎng)格:單分區(qū)正方體網(wǎng)格和三分區(qū)球頭網(wǎng)格.其中,正方體網(wǎng)格計(jì)算域?yàn)閤=[0,1],y=[0,1],z=[0,1],球頭半徑為0.5m.從圖2中兩種網(wǎng)格數(shù)值解和解析解的對(duì)比可以看出,二者符合良好,驗(yàn)證了數(shù)值方法對(duì)不同網(wǎng)格的適應(yīng)性.

驗(yàn)證算例2為間隔電極霍爾效應(yīng)算例,如圖3所示[13].±Y的通道邊界內(nèi)通入+X方向的導(dǎo)電流體,電導(dǎo)率σ=1??1·m?1.有限寬度的平行電極沿周期性重復(fù)布置在通道兩側(cè)壁面上.由于通道無(wú)限長(zhǎng),圖3中的Inlet和Outlet邊界設(shè)置為周期性邊界.為使通道內(nèi)流體速度場(chǎng)對(duì)電勢(shì)場(chǎng)結(jié)果無(wú)影響,則必須滿足?×(u×B)=0,因此可假定u=f(y),B=f(z),且通道內(nèi)流動(dòng)為充分發(fā)展的Poiseuille流動(dòng),如(16)式所示:

圖1 (網(wǎng)刊彩色)驗(yàn)證算例1的兩套網(wǎng)格(a)單分區(qū)立方體;(b)三分區(qū)1/4球頭Fig.1.(color on line)Two Meshes for validation case 1:(a)Cube of single segMent;(b)quarter sphere head of th ree segMents.

圖2 (網(wǎng)刊彩色)驗(yàn)證算例1的數(shù)值解和解析解對(duì)比(a)立方體網(wǎng)格;(b)球頭網(wǎng)格Fig.2.(color on line)CoMparison of nuMerical and analytical resu lts:(a)Cube Mesh;(b)quarter sphere head Mesh.

圖3 (網(wǎng)刊彩色)間隔電極霍爾效應(yīng)算例示意圖Fig.3.(color online)ScheMatic view of segMented electrodes.

其中,umax為通道中心線上的流體速度,umax=1m/s.h為通道半寬,h=0.5m.

當(dāng)B=0 T時(shí),兩電極間只存在靜電場(chǎng).圖4和圖5分別給出了無(wú)霍爾效應(yīng)時(shí)本文和文獻(xiàn)[17]的電勢(shì)等值線、電流流線對(duì)比.圖6給出了有霍爾效應(yīng)時(shí)(β=1.0,B=1 T)時(shí)的本文與文獻(xiàn)[13]的電流流線對(duì)比.可以看出,兩種情況下,本文計(jì)算結(jié)果與文獻(xiàn)[13]的結(jié)果符合良好.

圖4 (網(wǎng)刊彩色)無(wú)霍爾效應(yīng)時(shí)靜電場(chǎng)電勢(shì)等值線對(duì)比圖(a)本文結(jié)果;(b)文獻(xiàn)[13]結(jié)果Fig.4.(color on line)E lectric potential contou rs for the segMented electrode channelw ithou tMagnetic field:(a)Our results;(b)Ref.[13].

圖5 無(wú)霍爾效應(yīng)時(shí)靜電場(chǎng)電流線對(duì)比圖(a)本文結(jié)果;(b)文獻(xiàn)[13]結(jié)果Fig.5.E lectric cu rrent streaMlines w ithout Magnetic field:(a)Ou r resu lts;(b)Ref.[13].

圖6 有霍爾效應(yīng)(β=1.0)時(shí)靜電場(chǎng)電流線對(duì)比圖(a)本文結(jié)果;(b)文獻(xiàn)[13]結(jié)果Fig.6.E lectric current streaMlines w ith Hall eff ect(β=1.0):(a)Our resu lts;(b)Ref.[13].

4 霍爾電場(chǎng)收斂性分析

4.1 步進(jìn)因子a p

由(15)式可以得到i向的系數(shù)矩陣如(7)式所示.可以看出,步進(jìn)因子ap在迭代計(jì)算中的主要作用是保證迭代矩陣對(duì)角占有.因此,ap的合適取值應(yīng)和a,b,c三個(gè)系數(shù)的大小有關(guān).ap取值過(guò)大會(huì)造成收斂速度緩慢;過(guò)小則無(wú)法保證系數(shù)矩陣對(duì)角占優(yōu),使得計(jì)算發(fā)散.

根據(jù)(10)式和(12)式可以看出,系數(shù)b1∝ξnσ1i,n=x,y,z;i=1,2,3.因此ap的取值與網(wǎng)格尺度以及電導(dǎo)率、霍爾系數(shù)都有關(guān)系.網(wǎng)格越密,ξn越小,相應(yīng)的系數(shù)a,b,c越小,ap可以取值更小.這里選用驗(yàn)證算例1立方體網(wǎng)格進(jìn)行計(jì)算,分析不同網(wǎng)格尺度下ap的取值對(duì)電場(chǎng)收斂性的影響.

圖7為三套網(wǎng)格下(10×10×10,30×30×30,50×50×50)不同ap取值對(duì)應(yīng)的電勢(shì)收斂曲線.可以看出,對(duì)應(yīng)每套網(wǎng)格都存在一個(gè)最優(yōu)的步進(jìn)因子ap使得電勢(shì)收斂速度最快;隨著網(wǎng)格尺度的減小,最優(yōu)步進(jìn)因子減小,三套網(wǎng)格對(duì)應(yīng)的最優(yōu)步進(jìn)因子分別約為0.125,0.015,0.006,與之前理論分析一致.

圖7 (網(wǎng)刊彩色)不同網(wǎng)格尺度下不同a p取值的電場(chǎng)收斂曲線(a)10×10×10;(b)30×30×30;(c)50×50×50Fig.7.(color on line)Residual cu rves of electric potential under various stepping factors for diff erent grid densities:(a)10×10×10;(b)30×30×30;(c)50×50×50.

4.2 霍爾系數(shù)的影響

選用日本1996年發(fā)射的軌道再入試驗(yàn)飛行器(orbital reentry experimental,OREX)返回艙前體作為計(jì)算模型.計(jì)算域和網(wǎng)格分別如圖8所示.選取再入飛行時(shí)間在7461.5 s飛行工況(Ma=20.04,H=63.60 km).駐點(diǎn)磁感應(yīng)強(qiáng)度B0=0.2 T,霍爾系數(shù)β=1.0,5.0,10.0.采用全催化、絕緣壁面條件.

圖9給出了不同霍爾系數(shù)下不同ap的電勢(shì)收斂曲線.ap取值小于圖中的最小值時(shí)即發(fā)散.可以看出,霍爾系數(shù)β=1.0,5.0,10.0對(duì)應(yīng)的最優(yōu)ap分別約為1.0×10?3,2.0×10?3,4.0×10?3,并且當(dāng)收斂至電勢(shì)殘差為1.0×10?4時(shí),需要的虛擬迭代步數(shù)分別為900,8900,31000步.可見(jiàn),網(wǎng)格一定,對(duì)應(yīng)某個(gè)霍爾系數(shù)存在一個(gè)最優(yōu)的步進(jìn)因子ap,并且隨著霍爾系數(shù)的增加,最優(yōu)步進(jìn)因子ap增加,電勢(shì)場(chǎng)收斂速率變慢.

圖8 (網(wǎng)刊彩色)軌道再入試驗(yàn)飛行器的幾何模型、網(wǎng)格及邊界示意圖Fig.8.(color online)GeoMetry and boundary conditions of OREX.

圖9 (網(wǎng)刊彩色)不同霍爾系數(shù)下不同a p的電勢(shì)收斂曲線(a)β=1.0;(b)β=5.0;(c)β=10.0Fig.9.(color on line)Residual curvesunder variousstepping factors for diff erent HallparaMeters:(a)β=1.0;(b)β=5.0;(c)β=10.0.

5 變步進(jìn)因子加速法

由上節(jié)分析可知,步進(jìn)因子的最優(yōu)取值和網(wǎng)格尺度、霍爾系數(shù)相關(guān).不同網(wǎng)格尺度、霍爾系數(shù)下,最優(yōu)步進(jìn)因子差別較大.氣動(dòng)熱計(jì)算網(wǎng)格不均勻程度高并且霍爾系數(shù)也隨著外加磁感應(yīng)強(qiáng)度、飛行狀態(tài)的改變變化較大,這些都給固定步進(jìn)因子的取值帶來(lái)困難.因此,實(shí)際流場(chǎng)各點(diǎn)的ap可以取不同值,以在當(dāng)?shù)鼐W(wǎng)格尺度和當(dāng)?shù)鼗魻栂禂?shù)下都達(dá)到較優(yōu)的收斂效果.由于ap的主要作用是在計(jì)算中保證元素b1,b2,b3對(duì)角占優(yōu),故參考計(jì)算流體力學(xué)中當(dāng)?shù)貢r(shí)間步長(zhǎng)的做法,令

其中,ap0為常系數(shù).圖10給出了驗(yàn)證算例1局部加密立方體網(wǎng)格不同ap0下的電勢(shì)收斂曲線.其中,圖10(a)為X向加密網(wǎng)格,圖10(b)為三向均加密網(wǎng)格.圖中實(shí)心和空虛曲線分別代表全局定ap和當(dāng)?shù)刈僡p收斂曲線.可以看出,對(duì)于局部加密網(wǎng)格而言,變ap條件下的電勢(shì)收斂性明顯優(yōu)于定ap.盡管變ap條件下同樣需要確定最優(yōu)的步進(jìn)因子系數(shù)ap0,但其范圍相對(duì)比較固定.較優(yōu)的取值范圍位于0.01—0.50之間,具體取值要根據(jù)實(shí)際網(wǎng)格和霍爾系數(shù)確定.

圖10 (網(wǎng)刊彩色)不同步進(jìn)因子下的電勢(shì)收斂曲線(a)X向加密;(b)三個(gè)方向加密Fig.10.(color on line)Residual curves under various stepp ing factor for d iff erent refined Meshes:(a)Refined in the X direction;(b)refined in three directions.

6 結(jié)論

針對(duì)高速飛行器磁控?zé)岱雷o(hù)系統(tǒng)霍爾電場(chǎng)的求解問(wèn)題,采用交替隱式近似因子分解法,建立了熱化學(xué)非平衡條件下的電勢(shì)Poisson方程數(shù)值方法,并采用標(biāo)準(zhǔn)算例完成了程序驗(yàn)證.而后從迭代矩陣性質(zhì)出發(fā),分析了步進(jìn)因子ap和電場(chǎng)收斂性的關(guān)系以及影響步進(jìn)因子取值的因素(網(wǎng)格尺度、霍爾系數(shù)),并參考計(jì)算流體力學(xué)中當(dāng)?shù)貢r(shí)間步長(zhǎng)加速收斂的方法,提出了當(dāng)?shù)刈儾竭M(jìn)因子加速電場(chǎng)收斂方法.研究表明:

1)在網(wǎng)格和霍爾系數(shù)一定的情況下,存在一個(gè)最優(yōu)的步進(jìn)因子ap使得霍爾電場(chǎng)收斂速度最快;

2)隨網(wǎng)格尺度的減小和霍爾系數(shù)的增加,最優(yōu)步進(jìn)因子ap變大,電勢(shì)場(chǎng)收斂速率變慢;

3)對(duì)于局部加密網(wǎng)格而言,當(dāng)?shù)刈儾竭M(jìn)因子方法的電勢(shì)收斂性明顯優(yōu)于常規(guī)的定步進(jìn)因子法.

[1]Zhu Y J,Jiang Y S,Hua H Q,Zhang C H,X in C W 2014 Acta Phys.Sin.63 244101(in Chinese)[朱艷菊,江月松,華厚強(qiáng),張崇輝,辛燦偉2014物理學(xué)報(bào)63 244101]

[2]Y in J F,You Y X,LiW,Hu T Q 2014 Acta Phys.Sin.63 044701(in Chinese)[尹紀(jì)富,尤云祥,李巍,胡天群2014物理學(xué)報(bào)63 044701]

[3]Zhao G Y,Li Y H,Liang H,Hua W Z,Han MH 2015 Acta Phys.Sin.64 015101(in Chinese)[趙光銀,李應(yīng)紅,梁華,化為卓,韓孟虎2015物理學(xué)報(bào)64 015101]

[4]Yu H Y 2014 Acta Phys.Sin.63 047502(in Chinese)[于紅云2014物理學(xué)報(bào)63 047502]

[5]Bityu rin V A,Bocharov A N 2014 52nd Aerospace Sciences Meeting National Harbor,Mary land,January 13–17,2014

[6]Bisek N J,Gosse R,Poggie J 2013 J.Spacecraft Rockets 50 927

[7]C ristofolini A,BorghiC A,NerettiG,Battista F,Schettino A,Trifoni E,Filipp is F D,Passaro A,Baccarella D 2012 18th A IAA/3AF In ternational Space P lanes and Hypersonic SysteMs and Techno logies Conference Tours,France,Sep teMber 24–28,2012

[8]Lv H Y,Lee C H 2010 Chin.Sci.Bu ll.55 1182(in Chinese)[呂浩宇,李椿萱2010科學(xué)通報(bào)55 1182]

[9]Li K,Liu W Q 2016 Acta Phys.Sin.65 064701(in Chinese)[李開(kāi),劉偉強(qiáng)2016物理學(xué)報(bào)65 064701]

[10]Hu H Y,Yang Y J,Zhou W J 2011 Chin.J.Theo.App.Mechan.43 453(in Chinese)[胡海洋,楊云軍,周偉江2011力學(xué)學(xué)報(bào)43 453]

[11]Gaitonde D V,Poggie J 2002 40th A IAA Aerospace Sciences Meeting and Exhibit Reno,Nevada,January 14–17,2002

[12]W an T,Cand ler G V,Macheret SO,Shneider MN 2009 A IAA J.47 1327

[13]Bisek N J 2010 Ph.D.D issertation(Michigan:University of Michigan)

[14]Lv H Y,Lee C H,Dong H T 2009 Sci.Sin.Phys.Mechan.Astron.39 435(in Chinese)[呂浩宇,李椿萱,董海濤2009中國(guó)科學(xué)G輯39 435]

[15]Peng W,He Y G,Fang G F,Fan X T 2013 Acta Phys.Sin.62 020301(in Chinese)[彭武,何怡剛,方葛豐,樊曉騰2013物理學(xué)報(bào)62 020301]

[16]Fu jino T,MatsuMoto Y,Kasahara J,Ishikawa M2007 J.Spacecraft Rocket 44 625

[17]Holst T 2000 Progress in Aerospace Sci.36 1

[18]Zhang K P,D ing G H,T ian Z Y,Pan S,Li H 2009 J.National Univ.Defense Tech.31 39(in Chinese)[張康平,丁國(guó)昊,田正雨,潘沙,李樺2009國(guó)防科技大學(xué)學(xué)報(bào)31 39]

[19]T ian Z Y,Zhang K P,Pan S,Li H 2008 Chin.Quar.Mechan.29 72(in Chinese)[田正雨,張康平,潘沙,李樺2008力學(xué)季刊29 72]

[20]Gnoff o P A,Gup ta R N,Shinn J L 1989 NASA TP–2867

(Received 16 Sep teMber 2016;revised Manuscrip t received 22 January 2017)

PACS:47.40.Ki,47.85.L–,52.30.Cv,41.20.GzDOI:10.7498/aps.66.084702

*Project supported by the Natural Science Foundation of Hunan Province,China(G rant No.13JJ2002)and the National Natu ral Science Foundation of China(G rant No.90916018).

?Corresponding author.E-Mail:LiKai898989@126.com

N uMerical so lu tion p rocedu re for H all electric fi eld of the hyperson ic Magnetohyd rodynaMic heat sh ield system?

Li Kai?Liu Jun Liu Wei-Qiang

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

MagnetohydrodynaMic(MHD)heat shield systeMis a novel-concept thermal protection technique for hypersonic vehicles,which hasbeen proved by lotsof researchersw ith both nuMericaland experiMentalMethods.Most of researchers neglect the Hall eff ect in their researches.However,in the hypersonic reentry process,the Hall eff ect is soMetiMes so significant that the electric current distribution in the shock layer can be changed by the induced electric field.Consequently,the Lorentz force aswellas the Joule heat is varied,and thus the effi ciency of the MHD heat shield systeMis aff ected.

In order to analyze the influence of Halleff ect,the induced electric field must be taken into consideration.According to the weakly-ionized characteristics of hypersonic flow post bow shock,the Magneto-Reynolds number is assuMed to be sMall.Therefore,the Maxwell equations are siMp lified w ith the generalized Ohm’s law,and the induced electric field is governed by the potential Possion equation.Numericalmethods are hence established to solve the Hall electric field equations in the therMocheMical nonequilibriuMflow field.The electric potential Poisson equation is of significant rigidity and diffi cult to solve for two reasons.One is that the coeffi cientmatrix may not be diagonally doMinant when the Hall parameter is large in the shock layer,and the other is that thismatrix including the electric conductivity is discontinuous across the shock.In this paper,a virtual stepping factor is included to strengthen the diagonal doMinance and iMp rove the coMputational stability.Moreover,approximate factor and alternating direction iMp licit method are eMp loyed for further iMp roving the stability.W ith these Methods,a FORTRAN code is w ritten and validated by coMparing the nuMerical resultsw ith the analytical ones aswell as results available froMprevious references.A fter that,relation between the convergence property and the virtualstepping factor is revealed by theoreticalanalysisand numerical simulations.Based on thesework,a local variable stepping factorMethod is p roposed to accelerate the iterating process.Resu lts show that the convergence property is closely related to theMesh density and Hall paraMeter,and there exists a best stepping factor for a particu larmesh aswell as a particular Hall parameter.Since the best stepping factor varies a lot for diff erentMeshes and diff erent Hall paraMeter,its appropriate value is hard to choose.The best value of stepping factor coeffi cient still exists in the local step factor Method,but its value range is relatively sMaller.More iMportantly,the local stepping factor method yields better convergence property than the regular constant one when eMp loying a locally refined Mesh.

magnetohydrodynaMic flow control,thermal protection,Hall eff ect,Poisson equation

10.7498/aps.66.084702

?湖南省自然科學(xué)基金(批準(zhǔn)號(hào):13JJ2002)和國(guó)家自然科學(xué)基金(批準(zhǔn)號(hào):90916018)資助的課題.

?通信作者.E-Mail:LiKai898989@126.com

?2017中國(guó)物理學(xué)會(huì)C h inese P hysica l Society

http://w u lixb.iphy.ac.cn