高超声速湍流流动磁流体动力学控制机理*

2022-11-14罗仕超吴里银常雨

物理学报 2022年21期

罗仕超吴里银常雨

(中国空气动力研究与发展中心,超高速空气动力研究所,绵阳 621000)

基于低磁雷诺数假设建立完全气体湍流流场、磁场耦合模型.数值计算方法上,通过AUSMPW+格式和LUSGS 隐式处理方法求解磁流体动力学湍流流动方程,其中湍流模型采用Spalart-Allmaras 模型.分析了不同外加磁场条件下平板及压缩拐角湍流边界层流动控制效果.研究表明: 湍流边界层磁流体动力学流动控制效果与洛伦兹力大小正相关;外加磁场作用下,洛伦兹力的方向和流动方向相反,此时洛伦兹力起到减速的作用,减少了近壁面流体的动量,降低了边界层抵抗分离的能力;逆流向洛伦兹力减小了壁面的剪切应力,从而降低湍流流场壁面摩擦阻力系数,洛伦兹力对流体做负功,边界层内温度增加;磁相互作用位置对磁流体动力学分离区控制效果存在较大影响,工程应用中需配置合理的磁场布局方案.

1 引言

在超声速/高超声速飞行器的研制和使用过程中,均不可避免地会遇到由激波/边界层引起的各类问题,包括激波/边界层干扰、激波/激波干扰、进气道不起动等,这些可能导致推进系统工作效率下降、工作包线变窄,并使得飞行器气动阻力增加、局部出现高热流乃至烧蚀现象等[1,2].为此,有必要采取措施对流动进行控制.

变几何调节是目前研究较多的一类流动控制方法,其通过改变物面的几何位置、倾角等来实现对流动的控制.然而,由于变几何调节方法需要采用可转动或可平移的机械装置,这样不仅会使结构复杂、附加质量增加,还存在封严、热防护等问题.为此,各国学者一直在探索各种几何固定的流动控制方法,以在不改变飞行器几何外形的前提下对流动进行控制,避免变几何控制方法的不足.磁流体动力学(magnetohydrodynamic,MHD)流动控制是一项新型固定几何主动流动控制技术,可以用来解决高超声速飞行器面临的一系列关键问题,提升高超声速飞行器总体性能[3].电磁流动控制技术是采用“热电离”或人工电离方法提供等离子体流场,并在电磁激励的作用下对流场作用,从而实现对高超声速流场结构的控制.该固定几何流动控制方法具有非接触、响应快以及可重复利用等优点,并且可以根据来流状态进行参数调整,在控制方面有着很强的主动性,在高超飞行器设计领域具有潜在的广泛应用前景[4,5].

近20 年来,随着人工电离技术、超导磁体技术与计算机性能的突飞猛进,高超声速飞行器磁流体流动控制技术迎来了新的研究热潮[6-8],相关学者提出将MHD 技术应用于控制高超声速边界层的设想,并开展了大量的研究.2003 年,Bobashev 等[9]基于Ioffe 物理研究所激波风洞开展了高超声速边界层分离流动磁控实验研究.风洞实验段马赫数Ma=7.0,添加种子粒子后的电导率σe=1 700 S/m,磁场强度为1.5 T.结果表明,流场洛伦兹力与自由流方向相反时,流场由未分离状态转为分离状态.2004 年俄亥俄州立大学Meyer 等[10]在Ma=3.0超声速气流条件下,进行了洛伦兹力对边界层的影响实验研究,外加磁场强度为1.5 T.分析对比实验结果可知,加速洛伦兹力能够有效减少边界层内密度波动,使边界层流动趋于稳定,增强边界层抵抗分离的能力.普林斯顿大学Zaidi 等[11]在2006 年提出“雪橇(snowplow)”式边界层磁流体激励方式,用于边界层分离流动控制.东京工业大学Saito 等[12]在2008 年利用激波风洞,在Ma=1.5 超声速气流条件下,开展了磁流体边界层流动控制实验研究,实验结果表明,在减速洛伦兹力作用下,控制区域会出现显著的压升;加速洛伦兹力作用下,控制区域压力与无外加磁场工况基本一致;Saito 认为,加速洛伦兹力虽然可以减小局部静压,但是伴随产生的焦耳热效应会弱化加速洛伦兹力边界层分离控制效果.2013 年,Nagata 等[13]研究了外加磁场对双锥模型层流流动中激波-激波相互作用、激波-边界层相互作用的影响.磁场对分离泡及局部热流的特性有显著的影响,分离泡的大小随着磁场强度的增大而增大.2017 年,李益文等[14]建立了磁流体控制技术试验系统,采用电容耦合射频-直流组合放电对Ma=3.5 气流进行电离,在磁场作用下产生顺/逆气流方向的洛伦兹力控制流场.结果表明,在焦耳热和洛伦兹力的作用下,局部磁控流动减速时静压升高.2022 年,王江峰等[15]将电磁源项引入欧拉方程,建立了磁控进气道准一维简化模型,分析了不同外加磁场对典型进气道流场结构及性能影响.洛伦兹力作用下,分离区尺度减小,进气道内流结构显著改善.

从现有的研究来看,耦合磁场、等离子体流场、湍流的复杂多场高超声速流动控制研究难度极大,国内外尚未建立完善的理论体系和研究手段.耦合电磁场的高超声速湍流流场数值模拟研究方面,国内已经在该领域开展了部分研究工作,技术水平与国外相当.但与国外已开展研究相比,多场耦合模型研究的系统性、深入性,以及相关的计算流体动力学计算方法等研究手段还比较薄弱;考虑电磁效应的湍流模型及其数值计算方法有待进一步研究.电磁场耦合高超声速流动机理研究方面,针对多场耦合流动控制机理、规律研究开展较少,亦缺乏相应的理论基础研究成果;多场耦合模型相互作用机理有待进一步深入研究.

基于低磁雷诺数假设下的MHD 模型,建立了耦合电磁场的完全气体湍流流动MHD 数值模拟方法,针对典型构型开展等离子体流场与磁场耦合效应数值模拟研究,较为系统地分析不同外加磁场条件下湍流边界层流动控制效果,为发展电磁流动控制技术提供理论支撑.

2 完全气体湍流流场耦合外加磁场数值模拟方法

高超声速飞行器的实际飞行环境中,流体介质通常以低电导率(σmax≤102S/m)为特征,较低的电导率决定了低磁雷诺数.此时,感应磁场相对于外加磁场可以忽略不计,考虑电磁效应的Navier-Stokes(N-S)方程,即在低磁雷诺数假设下,通过在N-S 方程中添加电磁源项的方式来表征电磁场对导电流体运动的影响.在笛卡尔坐标系中,包含电磁源项SMFD的三维守恒形式的N-S 方程[16]为

其中,Q为守恒变量;F,G,H分别为x,y,z三方向对流通量项;Fv,Gv,Hv分别为x,y,z三方向黏性通量项.

qx,qy,qz分别为x,y,z三方向热流密度.外加磁场后附加了一个电磁源项,即

其中,J为电流密度矢量,B为磁感应强度矢量.完全气体无量纲的状态方程和总能、内能关系式为

式中,u,v,w分别为x,y,z三方向速度分量;ρ表示气体密度,T表示气体温度,Cv为定容比热;et为气体总能,e为气体内能,为平均分子量;Rg为气体常数,ℜ为阿伏伽德罗常数.

应力张量各分量为

热流密度为

在(7)式中,µ=µL+µT,其中µL为层流黏性系数(单位为N·s/m2),可以用Sutherland 公式求出:

µT为湍流黏性系数(由湍流模型计算),为自由流黏性系数,作为流场黏性系数无量纲的参考量.k为热传导系数,定义为

式中,Cp为等压比热,PrL与PrT分别为层流与湍流普朗特数.

从(3)式可以看出,考虑电磁场后,动量方程增加了洛伦兹力相关项J×B.能量方程增加了一个源项J·E,该源项由两部分组成,第一项为流动功,第二项为耗散焦耳热,具体表示为

忽略霍尔效应的假设[17,18]时,该电磁能量源项为0,流动功全部转化为焦耳热耗散,即J·E=0.此时,相对于常规的N-S 方程,额外引入的量只有电流密度矢量J和磁感应强度矢量B.不考虑感应磁场后,B由外加磁场给定.为使方程封闭,还必须引入广义欧姆定律和电导率模型.电导率可以退化为一个标量,且对于轴对称的外加磁场,可认为是无感应电场,并且激波层内的感应电流密度只有周向分量[19],如下式所示:

其中,σ表示流场电导率.

边界条件对于计算有重要的作用,边界条件的处理对计算的稳定性和数值误差有直接影响.本文考虑低磁雷诺数假设下的磁流体动力学控制方程,并且忽略霍尔效应.因此,流场的初始条件只需要给定外加磁场并且对流场进行初始化.计算过程中,磁场固定不变,边界条件只包括流动边界.本文以自由来流条件为初场,而速度则在壁面上给定无滑移条件,计算域其余部分则赋值为来流.计算区域一般包括远场/入口边界、出口边界、对称边界以及壁面边界.并行计算中,各分区还包括分区交界面边界.

远场/入口边界为自由来流条件:Φin=Φ∞.出口边界一般位于超声速区域,采用线性外插法:ΦN=2ΦN1-ΦN2.对称面边界,采用零梯度对称条件:=0,Φ为原始变量,Φin为远场/入口边界原始变量,Φ∞为自由流原始变量,ΦN1和ΦN2分别为第一层和第二层虚拟网格原始变量.

壁面边界具体包括:

1)速度无滑移条件:uWall=vWall=wWall=0;

2)等温壁条件:TW=TConst;

基于以上数学物理模型,我们采用有限体积法对方程(1)进行离散求解,对流项采用AUSMPW+(advection upstream splitting method by pressurebased weight functions)格式离散,黏性项采用二阶中心差分法进行离散,隐式求解基于LUSGS(lowerupper symmetric Gauss Seidel)方法,对电磁源项进行隐式处理以削弱源项过大而产生的刚性问题,从而可以提高程序的收敛性,湍流模型则选用SA(Spalart-Allmarars)方程.图1 是完全气体湍流流场、外加电磁场和霍尔电场耦合计算程序流程图,程序主要通过控制各参数来实现各类计算子程序模块的调用.

图1 完全气体磁流体湍流流场计算示意图(µt 为湍流黏性系数,kt 为湍流热传导系数,Fstep 为耦合计算过程流场迭代步数)Fig.1.Schematic of coupling computations of the magnetohydrodynamic turbulent flows based on the perfect gas model(µt is the turbulent viscosity coefficient,kt is the turbulent heat transfer coefficient,Fstep is the number of iteration steps of the flow field in the coupling calculation process).

3 数值计算方法验证及网格无关性

3.1 超声速平板湍流流动验证算例

为考核本文数值计算方法对湍流边界层预测的准确性,选用超声速平板湍流流动计算模型进行验证[20],计算工况见表1.此外,平板长度L=1.0 m,计算区域的高度为0.8 m,网格规模为101×69×3(流向×法向×周向).壁面附近法向网格加密,近壁面网格最小尺度y+＜1.0.

表1 超声速平板湍流流动计算自由流条件Table 1.Freestream conditions for the supersonic flat plate flow.

图2 为绝热壁面条件下SA 湍流模型计算得到的壁面律(速度剖面,x=1 m 处)和理论上的壁面律、参考文献[21]k-ω模型计算结果的对比.结果表明,SA 湍流模型能很好地模拟无分离平板湍流流动壁面律.图3 给出了绝热壁面条件下SA 湍流模型计算得到的壁面摩阻系数(Cf)与范德莱斯特(Van Driest)的理论值、参考文献[21]k-ω湍流模型计算值的比较.可以看出,在选定的网格间距下,本文SA 湍流模型计算得到的超声速平板湍流壁面摩阻系数与理论解、参考文献k-ω湍流模型计算值符合较好,本文数值方法能够较准确地模拟无分离平板湍流边界层流动.

图2 绝热壁面下SA 湍流模型计算得到的壁面律(速度剖面,x=1 m 处)和其他结果[21]的对比Fig.2.Comparison of the wall law calculated by the SA turbulence model(velocity profile,at x=1 m)under the adiabatic wall conditions with other results[21].

图3 绝热壁条件下,壁面热流分布对比Fig.3.Comparison of wall heat flow distribution under the adiabatic wall.

3.2 低磁雷诺数MHD 层流流动验证算例

选取圆球算例来验证低磁雷诺数MHD 层流流动计算模块的可靠性.表2 整理了圆球绕流算例计算的工况.

表2 圆球绕流算例计算自由流条件Table 2. Freestream conditions for the hemisphere flow.

外加磁场与流场相互作用产生洛伦兹力,激波层内流动减速,激波脱体距离增加.图4 给出了当外加磁场By=6.472 T 时,低磁雷诺数MHD 流场计算程序计算得到的流场压力轮廓图及文献[22]的结果,通过对比可知,本文计算结果与文献[22]结果符合良好.图5 是当外加磁场By=6.472 T 时无量纲压力沿驻点线分布的对比图.可以看出,本文MHD 程序计算结果与文献[22]计算结果一致,验证了本文耦合磁场的低磁雷诺数MHD 层流流动计算模块的可靠性.

图4 By=6.472 T 时的流场压力轮廓图(a)本文计算;(b)文献[22]的结果Fig.4.Pressure profile of flow field with By=6.472 T:(a)Paper results;(b)results from Ref.[22].

图5 By=6.472 T,压力沿驻点线线压力分布对比Fig.5.Comparison of the pressure distribution along axial line with By=6.472 T.

3.3 网格无关性分析

选用轨道再入试验飞行器(orbital reentry experiment,OREX)验证程序的网格无关性.令试验速度U∞=3873.4 m/s,高度h=51.99 km,OREX外形见图6.建立了六套网格进行网格收敛性分析,网格雷诺数分别为0.8,1.6,5.5,11.0,28.0,56.0;图7 对比了六套网格下的壁面热流计算结果,基于自由流参数网格雷诺数ReΔn,∞≤1.6 时,相应网格下的驻点热流峰值计算结果相对误差不超过1%,程序网格无关性得到了有效验证.

图6 OREX 外形Fig.6.Geometry for the OREX.

图7 不同网格下壁面热流计算结果对比Fig.7.Heat flux distributions under different normal grid spacing at the wall.

4 平板磁流体湍流流动

选取典型的低磁雷诺数MHD 平板流动算例,以分析不同外加磁场条件下平板湍流边界层流动控制效果.超声速平板磁流体湍流流动计算工况如表3 所列,入口边界为给定的自由来流条件.此外,外加均布磁场作用于整个平板流动计算域,磁场方向垂直于平板,磁感应强度Bmax=0—1.3 T,此时洛伦兹力方向与流动方向相反.外加磁场平板湍流计算网格规模为100×80×3(流向×法向×周向),近壁面附近法向网格加密,第一层网格的高度为2.0×10-7m.

表3 超声速平板磁流体湍流流动计算自由流条件Table 3. Freestream conditions for the flat plate magnetohydrodynamic turbulent flows.

首先进行耦合磁场的平板层流流动算例分析,将计算结果与文献[23]中平板层流边界层速度剖面进行比较.图8 为不同磁感应强度下的平板层流边界层速度剖面图,纵轴表示平板无量纲法向距离,无量纲参考长度为平板长度L,y为壁面法向距离,计算结果和文献[23]符合良好(误差在1%以内).结果表明,本文程序能准确可靠地模拟低磁雷诺数下的平板磁流体层流流动.

图8 x=0.06 m 时,不同磁感应强度下的层流边界层速度分布的对比Fig.8.Comparison of laminar velocity profiles under the different magnetic field strength with x=0.06 m.

图9 和图10是x=0.06 m 位置不同外加磁场作用下,平板湍流边界层速度、温度剖面的对比图.外加磁场作用下,洛伦兹力的方向和流动方向相反,此时洛伦兹力起到减速作用,边界速度型面更陡峭;此外,洛伦兹力对流体做负功,在外加磁场作用下边界层内温度增加;湍流边界层流动控制效果与洛伦兹力大小正相关.图11 给出了外加磁场作用下层流/湍流摩擦阻力系数的对比,给定x=0.04 m为流场强制转捩点,转捩点前计算流场全为层流流动,转捩点后计算流场全为湍流.磁场与流场相互作用产生的逆流向洛伦兹力减小了壁面的剪切应力,从而降低层流/湍流流场壁面摩擦阻力系数.施加1.0 T 磁场时,x=0.03 m 站位层流流场局部壁面摩擦系数减小25.5%;x=0.06 m 站位湍流流场局部壁面摩擦系数最大减小7.2%.

图9 x=0.06 m 时,不同磁感应强度下的湍流边界层速度型面对比Fig.9.Comparison of turbulent velocity profiles under the different magnetic field strength with x=0.06 m.

图10 x=0.06 m 时,不同磁感应强度下的湍流边界层温度剖面对比Fig.10.Comparison of the turbulent temperature profiles under the different magnetic field strength with x=0.06 m.

图11 x=0.06 m 时,不同磁感应强度下湍流摩阻系数对比Fig.11.Comparison of turbulent skin friction coefficients under the different magnetic field strength with x=0.06 m.

5 高超声速压缩拐角磁流体湍流流动

压缩拐角是一类经典的激波/边界层干扰问题,虽然外形比较简单,但局部的激波、边界层分离等流动现象复杂,给数值模拟带来了很大的挑战.选取Settles 和Dodson[24]实验数据库中的34°高超声速压缩拐角为研究对象,进行激波/边界层干扰磁流体湍流数值模拟研究.自由流条件如表4 所列,高超声速压缩拐角计算网格数为149× 180×4(流向×法向×周向),近壁面区域、拐角区域网格加密,近壁面区域网格最小尺度y+＜1.0.

表4 34°压缩拐角磁流体湍流流动自由流条件Table 4. Freestream conditions for the 34° ramp magnetohydrodynamic turbulent flows.

高超声速拐角流动分离控制系统包含外加磁场系统以及人工电离系统(电子枪)等设备.由于来流工况与电离设备的制约,绕流介质所能获得的电导率较低,属于低磁雷诺数磁流体流动.文献[25]指出,当电子束电流密度为 5 mA/cm2,电子束能为30 keV 时产生了约 5.0 S/m 的电导率,为简化计算,电导率的分布简化处理为均匀分布.图12 是压缩拐角流动控制系统局部受力机理分析原理图,在外加均布磁场条件下对压缩拐角湍流流动控制局部受力进行机理分析.U为流体的速度矢量,B为磁感应强度矢量.外加磁场和导电流体相互作用产生一个垂直X-Y平面向外的感应电流J,电流和磁场相互作用产生一个逆流向的洛伦兹力F,逆流向洛伦兹力F主要作用是减速流体.

图12 34°压缩拐角MHD 流动控制系统局部受力机理分析Fig.12.Analysis of local force mechanism of 34° the ramp compression corner MHD flow control system.

MHD 控制压缩拐角湍流边界层分离,宏观上表现为分离泡尺度变化,分离激波和再附激波位置改变,本质上是外加磁场与导电流场相互作用产生的洛伦兹力改变边界层速度,从而拐角点附近压升的幅度发生变化.为研究电磁相互作用位置对压缩拐角流动MHD 控制效果影响,选定控制区域位于平板区域,共设计了4 组计算算例,具体电磁控制参数如表5 所列.其中,无外加磁场条件下,34°高超声速压缩拐角湍流流场分离点位于xS=0.550 m,再附点位于xR=0.569 m;x1为磁相互作用位置的起点,x2为磁相互作用位置的终点;控制区域局部电导率高度取10 mm.

表5 电磁流动控制参数Table 5. Electromagnetic flow control parameters.

图13(a)是平板区域不同MHD 作用位置壁面压力分布对比结果.可清楚地看到: 无控制时流动发生小分离,分离泡(压力平台区)尺度较小;外加磁场作用下,压力平台区域明显增大,Case 1 控制效果最好;不同MHD 作用位置对分离流动控制效果影响很大,在不同位置施加MHD 控制时,分离泡长度最大相差45.2%.特定控制区域内,MHD区域越靠近原分离点上游,控制效果越显著,压力平台区显著增大.图13(b)为平板区域不同MHD作用位置壁面摩阻系数对比结果,随着MHD 区域前移,分离点向上游移动,再附点向下游移动,分离区尺寸变大.外加磁场作用下,控制区域逆流向洛伦兹力减小了壁面的剪切应力,从而降低了湍流流场壁面摩擦阻力系数.

图13 MHD 作用位置对壁面压力系数及摩阻系数分布影响(a)壁面压力系数分布;(b)摩阻系数分布Fig.13.Effect of MHD position on wall pressure coefficient and skin friction coefficient distribution:(a)Wall pressure distribution;(b)skin friction coefficient distribution.

图14(a)是不同MHD 作用位置时,x=0.4 m站位截面速度分布,可看出施加逆流向洛伦兹力能明显减速边界层,使近壁面区速度型面更加陡峭,减少了局部控制区域近壁面流体的动量,从而降低了边界层抵抗分离的能力.图14(b)为不同MHD作用位置时,x=0.4 m 站位截面温度分布,等温壁面边界条件下,逆流向洛伦兹力,对流体做负功,边界层内温度增加.

图14 MHD 作用位置对x=0.4 m 站位截面速度及温度分布影响(a)速度剖面;(b)温度剖面Fig.14.Effect of MHD position on velocity and temperature distribution at x=0.4 m station section:(a)Profiles of velocities;(b)profiles of temperatures.