高超声速飞行器磁控热防护霍尔电场数值方法研究∗

李开 柳军 刘伟强

(国防科学技术大学航天科学与工程学院,长沙410073)

高超声速飞行器磁控热防护霍尔电场数值方法研究∗

李开†柳军 刘伟强

(国防科学技术大学航天科学与工程学院,长沙410073)

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

作为一种新概念高超声速热防护手段,磁控热防护系统在实际应用中往往需要考虑霍尔效应的影响.为了在真实气体环境下求解霍尔电场,采用交替隐式近似因子分解法建立并验证了热化学非平衡流体域电场数值求解方法.分析了电场虚拟步进因子和收敛性的关系以及影响步进因子取值的因素,提出了当地变步进因子加速电场收敛方法.研究表明,存在一个最优的步进因子ap使得霍尔电场收敛速度最快,并且随网格尺度的减小和霍尔系数的增加,最优步进因子ap变大,电势场收敛速率变慢.对于局部加密网格而言,当地变步进因子法的电势收敛性明显优于常规的定步进因子法.

磁流体控制,热防护,霍尔效应,泊松方程

1 引言

热防护系统是保护高超声速飞行安全的关键子系统,一直以来都是高超声速飞行器系统设计的关键一环[1].由于电磁流动控制技术近年来得到了飞速发展[2−4],作为其在高超声速热防护领域的一个应用,磁控热防护技术显示出了巨大的潜在应用价值,得到了各国科研机构的广泛关注[5−7].

对于高超声速飞行器前缘磁控热防护而言,其核心在于在高超声速飞行器鼻锥内部添加磁铁,通过外加磁场和波后带电流体的相互作用,感应产生洛伦兹力推出激波、减速流体、偏转流线,从而降低到达壁面的热流,实现热防护.由于存在霍尔效应,磁场对电子的作用会使流动的法向和流向存在电势差,即感应霍尔电场.根据通用欧姆定律,感应电场的存在会使得电流出现三维效应,进而改变洛伦兹力,影响磁控热防护系统的效率.为分析霍尔效应对磁控热防护系统的影响,必须耦合求解热化学非平衡流场和霍尔电场.高超声速条件下,由于波后电离度较低,低磁雷诺数假设往往可以得到满足[8].此时,电磁场与非平衡流场的耦合计算可以通过求解电势泊松方程和包含电磁源项的非平衡流动方程实现[9].流体域霍尔电场求解的关键在于求解电势泊松方程.由于作为系数的电导率在激波位置存在间断,并且在大霍尔系数条件下系数矩阵刚性问题严重[10],高超声速流体域电势泊松方程的求解存在鲁棒性差、收敛性差等问题.Gaitonde等[11]和Wan等[12]分别基于有限差分法和有限体积法,采用虚拟时间步长以及近似因子分解-交替隐式方法实现了三维电势泊松方程的数值离散与求解.Bisek[13]基于有限体积法采用连续超松弛算法(successive over-relaxation,SOR)完成了三维电势泊松方程的迭代求解.吕浩宇等[14]、彭武等[15]同样基于有限差分法采用超松弛算法求解了二维电势泊松方程.Fujino等[16]采用基于Galerkin有限元素法求解二阶半离散电势泊松方程,与实验结果符合良好.胡海洋等[10]采用无虚拟时间步的LUSGS预处理BI-CGSTAB算法解决了大霍尔系数下泊松方程病态矩的求解效率低的问题.综上,目前见于文献的霍尔电场求解方法主要有三种:1)超松弛迭代显式算法;2)近似因子分解-交替隐式方法;3)有限元方法,如BI-CGSTAB,Galerkin方法等.相比而言,SOR算法为显式算法,应用简单且利于并行计算,但收敛效率低于隐式格式.有限元方法收敛性好,但在处理复杂壁面条件时存在困难.

目前关于霍尔电场的研究往往仅限于Poisson方程求解方法的构建,对霍尔电场收敛性的研究较少.针对霍尔电场收敛性和步进因子ap的关系,Wan等[12]给出了其建议取值范围0.001—0.01;Gaitonde[11]则采用了Holst[17]提出的步进因子根据其上下限随时间步过渡变化的指数公式.但是,Wan给出的范围只针对其所采用的网格,不具有普适性.Holst的方法也必须首先确定步进因子的上下限,实际操作困难.因此,本文首先构造三维霍尔电场的隐式数值方法,而后进行霍尔电场收敛性的影响因素研究,以便为高超声速飞行器前缘磁控热防护系统的电磁场流场耦合计算分析提供参考.

2 数学模型

2.1 控制方程

高超声速飞行器定常绕流条件下,由于波后等离子体电离度往往低于1×10−4[18,19],此时,低磁雷诺数假设通常可以得到满足,并且外加稳定磁场时,磁场和感应电场的变化通常都可以忽略不计.因此,可以对安培-麦克斯韦公式、法拉第定律进行简化,将对感应场强矢量E以及感应电流密度J的求解转化为对标量电势ϕ的求解,结合广义欧姆定律

可以得到关于ϕ的电势泊松方程

其中u为流体速度矢量,B为磁感应强度矢量,电导率张量˜σ可以由下式求出:

其中,β为霍尔系数,并且为分析霍尔系数对步进因子取值、电场计算收敛性的影响,采用均布霍尔系数模型.σ为标量电导率,采用Fujino电导率模型计算[16]:

其中,e,ne,Me分别为电子电量、电子数密度、电子质量.为电子和其他组元的有效动量传输碰撞频率[20].

2.2 数值方法

由于电导率σ在激波处存在间断,作为系数的电导率张量在激波处也同样存在间断,这使得上述电势泊松方程(2)在求解时稳定性差.特别是在大霍尔系数条件下,系数矩阵很有可能会出现对角不占优的情况,使得方程刚性严重,影响迭代计算的收敛性.为此,我们采用能有效抑制刚性的近似因子分解-交替隐式方法.并且为了保证三维多分区网格在数值离散中的几何守恒性,将基于单元中心有限体积法对上述泊松方程进行离散.过程如下.

对于特定体积V,(1)式可以化为

将(6)式中右侧展开,可得到

其中|S|和(n1,n2,n3)分别是单元面的面积及法向矢量.令

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

对∇ϕ在三个计算坐标方向做分解并应用链式法则,则有

其中,ξ,η,ζ分别是i,j,k三方向的坐标变换系数.令XIC,ETC,ZTC如(10)式所示:

采用虚拟时间推进法进行迭代求解,则(6)式可以化为

应用中心差分对(6)式进行离散,得到(i,j,k)及其周围十八点的差分表达式.对于ϕi,j,k及临近的六个点ϕi+1,j,k,ϕi−1,j,k,ϕi,j+1,k,ϕi,j−1,k,ϕi,j,k+1,ϕi,j,k−1进行隐式处理ϕ(n+1)=ϕ(n)+∆ϕ(n),并令

a2j−1,b2j,c2j+1,a3k−1,b3k,c3k+1形式类似,不再赘述.令步进因子ap=1/∆t,则(11)式可以化为

令(13)式中右端项为RE,左端项采用交替隐式-近似因子分解法(ADI-AF),并定义算子Ai,Bj,Ck如(14)式所示:

则最终的电势方程迭代计算式为

霍尔电场边界条件为:绝缘壁面,J·n=0;导电壁面,ϕ=0;其余边界为NeuMann边界,∇ϕ·n=0.

3 数值验证

验证算例1取自参考文献[13],其控制方程为∆ϕ=xez.该方程拥有解析解ϕ=xez.若令则其控制方程可以化为和(2)式同样的简化形式.给定Dirichlet边界条件:ϕ=xez.图1给出了两套计算网格:单分区正方体网格和三分区球头网格.其中,正方体网格计算域为x=[0,1],y=[0,1],z=[0,1],球头半径为0.5m.从图2中两种网格数值解和解析解的对比可以看出,二者符合良好,验证了数值方法对不同网格的适应性.

验证算例2为间隔电极霍尔效应算例,如图3所示[13].±Y的通道边界内通入+X方向的导电流体,电导率σ=1Ω−1·m−1.有限宽度的平行电极沿周期性重复布置在通道两侧壁面上.由于通道无限长,图3中的Inlet和Outlet边界设置为周期性边界.为使通道内流体速度场对电势场结果无影响,则必须满足∇×(u×B)=0,因此可假定u=f(y),B=f(z),且通道内流动为充分发展的Poiseuille流动,如(16)式所示:

图1 (网刊彩色)验证算例1的两套网格(a)单分区立方体;(b)三分区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 (网刊彩色)验证算例1的数值解和解析解对比(a)立方体网格;(b)球头网格Fig.2.(color on line)CoMparison of nuMerical and analytical resu lts:(a)Cube Mesh;(b)quarter sphere head Mesh.

图3 (网刊彩色)间隔电极霍尔效应算例示意图Fig.3.(color online)ScheMatic view of segMented electrodes.

其中,umax为通道中心线上的流体速度,umax=1m/s.h为通道半宽,h=0.5m.

当B=0 T时,两电极间只存在静电场.图4和图5分别给出了无霍尔效应时本文和文献[17]的电势等值线、电流流线对比.图6给出了有霍尔效应时(β=1.0,B=1 T)时的本文与文献[13]的电流流线对比.可以看出,两种情况下,本文计算结果与文献[13]的结果符合良好.

图4 (网刊彩色)无霍尔效应时静电场电势等值线对比图(a)本文结果;(b)文献[13]结果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 无霍尔效应时静电场电流线对比图(a)本文结果;(b)文献[13]结果Fig.5.E lectric cu rrent streaMlines w ithout Magnetic field:(a)Ou r resu lts;(b)Ref.[13].

图6 有霍尔效应(β=1.0)时静电场电流线对比图(a)本文结果;(b)文献[13]结果Fig.6.E lectric current streaMlines w ith Hall eff ect(β=1.0):(a)Our resu lts;(b)Ref.[13].

4 霍尔电场收敛性分析

4.1 步进因子a p

由(15)式可以得到i向的系数矩阵如(7)式所示.可以看出,步进因子ap在迭代计算中的主要作用是保证迭代矩阵对角占有.因此,ap的合适取值应和a,b,c三个系数的大小有关.ap取值过大会造成收敛速度缓慢;过小则无法保证系数矩阵对角占优,使得计算发散.

根据(10)式和(12)式可以看出,系数b1∝ξnσ1i,n=x,y,z;i=1,2,3.因此ap的取值与网格尺度以及电导率、霍尔系数都有关系.网格越密,ξn越小,相应的系数a,b,c越小,ap可以取值更小.这里选用验证算例1立方体网格进行计算,分析不同网格尺度下ap的取值对电场收敛性的影响.

图7为三套网格下(10×10×10,30×30×30,50×50×50)不同ap取值对应的电势收敛曲线.可以看出,对应每套网格都存在一个最优的步进因子ap使得电势收敛速度最快;随着网格尺度的减小,最优步进因子减小,三套网格对应的最优步进因子分别约为0.125,0.015,0.006,与之前理论分析一致.

图7 (网刊彩色)不同网格尺度下不同a p取值的电场收敛曲线(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 霍尔系数的影响

选用日本1996年发射的轨道再入试验飞行器(orbital reentry experimental,OREX)返回舱前体作为计算模型.计算域和网格分别如图8所示.选取再入飞行时间在7461.5 s飞行工况(Ma=20.04,H=63.60 km).驻点磁感应强度B0=0.2 T,霍尔系数β=1.0,5.0,10.0.采用全催化、绝缘壁面条件.

图9给出了不同霍尔系数下不同ap的电势收敛曲线.ap取值小于图中的最小值时即发散.可以看出,霍尔系数β=1.0,5.0,10.0对应的最优ap分别约为1.0×10−3,2.0×10−3,4.0×10−3,并且当收敛至电势残差为1.0×10−4时,需要的虚拟迭代步数分别为900,8900,31000步.可见,网格一定,对应某个霍尔系数存在一个最优的步进因子ap,并且随着霍尔系数的增加,最优步进因子ap增加,电势场收敛速率变慢.

图8 (网刊彩色)轨道再入试验飞行器的几何模型、网格及边界示意图Fig.8.(color online)GeoMetry and boundary conditions of OREX.

图9 (网刊彩色)不同霍尔系数下不同a p的电势收敛曲线(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 变步进因子加速法

由上节分析可知,步进因子的最优取值和网格尺度、霍尔系数相关.不同网格尺度、霍尔系数下,最优步进因子差别较大.气动热计算网格不均匀程度高并且霍尔系数也随着外加磁感应强度、飞行状态的改变变化较大,这些都给固定步进因子的取值带来困难.因此,实际流场各点的ap可以取不同值,以在当地网格尺度和当地霍尔系数下都达到较优的收敛效果.由于ap的主要作用是在计算中保证元素b1,b2,b3对角占优,故参考计算流体力学中当地时间步长的做法,令

其中,ap0为常系数.图10给出了验证算例1局部加密立方体网格不同ap0下的电势收敛曲线.其中,图10(a)为X向加密网格,图10(b)为三向均加密网格.图中实心和空虚曲线分别代表全局定ap和当地变ap收敛曲线.可以看出,对于局部加密网格而言,变ap条件下的电势收敛性明显优于定ap.尽管变ap条件下同样需要确定最优的步进因子系数ap0,但其范围相对比较固定.较优的取值范围位于0.01—0.50之间,具体取值要根据实际网格和霍尔系数确定.

图10 (网刊彩色)不同步进因子下的电势收敛曲线(a)X向加密;(b)三个方向加密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 结论

针对高速飞行器磁控热防护系统霍尔电场的求解问题,采用交替隐式近似因子分解法,建立了热化学非平衡条件下的电势Poisson方程数值方法,并采用标准算例完成了程序验证.而后从迭代矩阵性质出发,分析了步进因子ap和电场收敛性的关系以及影响步进因子取值的因素(网格尺度、霍尔系数),并参考计算流体力学中当地时间步长加速收敛的方法,提出了当地变步进因子加速电场收敛方法.研究表明:

1)在网格和霍尔系数一定的情况下,存在一个最优的步进因子ap使得霍尔电场收敛速度最快;

2)随网格尺度的减小和霍尔系数的增加,最优步进因子ap变大,电势场收敛速率变慢;

3)对于局部加密网格而言,当地变步进因子方法的电势收敛性明显优于常规的定步进因子法.

[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)[朱艳菊,江月松,华厚强,张崇辉,辛灿伟2014物理学报63 244101]

[2]Y in J F,You Y X,LiW,Hu T Q 2014 Acta Phys.Sin.63 044701(in Chinese)[尹纪富,尤云祥,李巍,胡天群2014物理学报63 044701]

[3]Zhao G Y,Li Y H,Liang H,Hua W Z,Han MH 2015 Acta Phys.Sin.64 015101(in Chinese)[赵光银,李应红,梁华,化为卓,韩孟虎2015物理学报64 015101]

[4]Yu H Y 2014 Acta Phys.Sin.63 047502(in Chinese)[于红云2014物理学报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科学通报55 1182]

[9]Li K,Liu W Q 2016 Acta Phys.Sin.65 064701(in Chinese)[李开,刘伟强2016物理学报65 064701]

[10]Hu H Y,Yang Y J,Zhou W J 2011 Chin.J.Theo.App.Mechan.43 453(in Chinese)[胡海洋,杨云军,周伟江2011力学学报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中国科学G辑39 435]

[15]Peng W,He Y G,Fang G F,Fan X T 2013 Acta Phys.Sin.62 020301(in Chinese)[彭武,何怡刚,方葛丰,樊晓腾2013物理学报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)[张康平,丁国昊,田正雨,潘沙,李桦2009国防科技大学学报31 39]

[19]T ian Z Y,Zhang K P,Pan S,Li H 2008 Chin.Quar.Mechan.29 72(in Chinese)[田正雨,张康平,潘沙,李桦2008力学季刊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

∗湖南省自然科学基金(批准号:13JJ2002)和国家自然科学基金(批准号:90916018)资助的课题.

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

©2017中国物理学会C h inese P hysica l Society

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