激波内部结构的数值求解方法1)

朱清波 周文元 杨庆春 徐 旭

* (北京航空航天大学宇航学院,北京 102206)

† (上海空间推进研究所,上海 201112)

引言

在气体动力学中,激波常被简化为零厚度的间断面,经过激波,流动参数发生阶跃,熵不可逆地增加,这对大部分情况是很好的近似.然而真实的宏观流动中不会出现数学意义上的不连续,激波也不例外,从波前向波后状态的过渡必然是连续的,而激波厚度则是一个有限值.图1 展示了一个典型驻定正激波的内部流动参数变化曲线,其中Ma1和Ma2分别为波前和波后马赫数,δ 为基于最大速度梯度定义的激波厚度.

图1 激波内部结构示意图Fig.1 Schematic diagram of the internal structure of shock waves

研究激波内部流动参数变化过程的问题被称为激波结构问题,它不仅是一个重要的科学问题,还具有工程实用价值.在稀薄气体动力学[1]、分子动理论[2]、非平衡态热力学[3]和天体物理[4]等基础科学领域中,激波结构问题促进了对激波、输运现象以及流体力学理论基础的理解.在航空航天[5-6]、光学[7]以及核[8]等技术领域中,激波结构的正确解析对解决一些工程问题至关重要.例如航天器再入[5-6]时,高层大气极其稀薄导致激波厚度较大,准确模拟激波内部结构是外流场仿真的关键之一;又比如,激波会使光发生偏折,正确描述激波结构有助于修正它对光学系统 (风洞试验的光学测量装置、超声速飞行器的光学观瞄设备与机载激光武器等) 产生的不利影响[7].此外,驻定正激波是一个一维定常且边界条件极其简单的强非平衡流动,这使其成为验证气体动力学模型及相关数值方法准确性与可靠性的理想场景,再加上丰富的实验数据[9-10],激波结构长期被用作严苛而敏感的标准测试算例[11-14].

激波结构的求解有解析和数值两种方式.一般来说,包括激波结构问题在内的大部分流动问题非常复杂,难以求得解析解.但Becker[15]发现当气体的Prandtl 数等于3/4 这个特殊值时,激波内部流动所满足的能量方程可以被简化,总焓沿流向保持不变,进而获得了第一个对任意强度激波有效的解析解.Becker 解随后被不断改进,从最初的只适用于输运系数为常数的气体,到适用于硬球和Maxwell 分子模型气体[16]、软球分子模型气体[17]和van der Waals 气体[18]等.尽管适用范围不断扩大,这些解析解对Prandtl 数的要求却从未放宽.另一个被广泛使用的解析解是Mott-Smith[19]提出的著名的双峰分布函数解,然而它只是一定假设下的近似解,且从原理上看不能很好地描述弱激波.实际上,即使对强激波,Mott-Smith 解也不被实验和直接模拟Monte Carlo(DSMC) 计算的结果支持[10].总而言之,目前的解析方法仍存在较多局限性,只能解决一些特殊条件下的简单问题,而大部分激波结构问题较为复杂,需借助数值方法求解.

随着计算科学的发展,数值技术在流体力学中发挥着越来越重要的作用,激波结构问题也受益于此.按所采用的数学物理模型的类型来划分,激波结构的数值模拟方法主要有3 种:

(1) 微观方法,即基于还原论观点、从第一性原理出发直接模拟分子运动和碰撞的粒子方法,主要包括以分子动力学[20](MD) 为代表的确定论方法和以DSMC[21-22]为代表的概率方法;

(2) 介观方法,求解基于速度分布函数描述的Boltzmann 方程[23-24],或其简化方程/模型[25],包括著名的BGK 模型方程[26]和离散速度模型[27-28](DVM) 等;

(3) 宏观方法,求解基于宏观量描述的连续介质流体力学方程/模型,包括经典和改进的Navier-Stokes 方程[11,29-30]、以Burnett 方程为代表的高阶流体力学方程[31-32]、Grad 矩方程及其变体[14,33]、广义流体力学[34](generalized hydrodynamics) 方程、非线性耦合本构关系[35](NCCR) 模型以及扩展/广延热力学[3](extended thermodynamics) 模型等.

相对于Boltzmann 方程,控制气体宏观运动的流体力学方程是一种降阶的描述,自由度更少,因而更简单.例如在Chapman-Enskog 理论[36]中,Navier-Stokes 方程和Burnett 方程分别是对Boltzmann 方程的一阶和二阶近似.因此,第3 类方法原则上应归入第2 类.但从介观Boltzmann 方程到宏观流体力学方程的粗化并不容易,两者之间的关系在数学上也远非显而易见 (详见Hilbert 第6 问题[37]).更重要的是,流体力学方法基于与分布函数描述迥然不同的局部统计平均量 (分布函数的矩) 描述,体系相对独立而内容又极其丰富,故单独列为一类.

前述3 类方法中,基于流体力学方程的数值模拟最受青睐,因为它在具有良好精度的同时还有较高的计算效率,其宏观描述也为人们乐于接受,因此相关算法得到高度发展.具体到激波结构的仿真计算中,有限差分法常被使用,这类通用的计算流体力学技术十分成熟,此处不再赘述.

由于以下3 个原因,绝大多数激波结构问题可归结为一维定常正激波的结构问题.

(1) 激波的典型厚度非常小,流体穿过激波的时间极短,一般仅需数次分子碰撞即可完成从波前到波后状态的过渡,弛豫速度远快于其他宏观非平衡流动现象.因此,对于大部分激波结构问题,运动随时间的发展不重要,非定常效应可以忽略.

(2) 平面激波流动可以被视为法向的正激波流动叠加切向的平行流动.因此,二维或三维平面激波的结构问题可被简化为法向上的一维正激波结构问题.

(3) 曲面激波波面的最小曲率半径一般远大于激波厚度,所以曲率对激波内部结构的影响通常可以忽略,曲面激波的局部可被近似为平面激波,进而分解为法向上的正激波和切向的平行流.

毫无疑问,一维定常正激波是检验气体动力学模型 (尤其是稀薄气体动力学模型) 的理想算例: 极小的厚度意味着Knudsen 数较大,间断分子效应显著;强耗散意味着远离热力学平衡,非平衡效应显著;一维定常的特性和简单的边界条件还降低了问题的复杂性.定常正激波的结构问题,其本质是求解一组流体力学常微分方程.对这类问题,直接数值积分比有限差分法更简单高效.然而对激波结构应用空间推进的数值积分仍存在一些问题,例如沿流向推进计算的发散问题,本文的目的是解决这些问题,并建立在流体力学框架内数值求解激波结构的一般程序.

由于关注的问题具有普遍性,本文将只讨论最基础的情况,即单组分简单气体的情况,多组分气体和等离子体的情况暂不考虑,以避免不必要的细节把问题复杂化.控制方程将采用简单的Navier-Stokes 方程,这不会导致本文的讨论失去一般性,因为宏观连续的流体系统都具有Navier-Stokes 形式的守恒方程,区别仅在于用怎样的本构方程封闭它们,而这不影响本文的主要结论.

1 控制方程与边界条件

一维气体流动的质量、动量和能量守恒方程分别可写为

其中,t,x,ρ,V,p,u,ht,τ 和q分别为时间、流向坐标、密度、速度、压强、比内能、比总焓、黏性应力以及热流密度.需要说明的是,为保持与气体运动论领域的习惯一致,本文约定当 τ 表现为压应力时其符号为正,这与流体力学中法向应力的定义相反.式(1)～式(3) 直接来源于三大守恒定律,也可从Boltzmann 方程的矩方程——Maxwell 输运方程导出,未对 τ 和q作任何假设,因此被所有连续介质气体动力学模型的一维形式所遵守.

若采用Navier-Stokes 和Fourier 线性本构关系,有

式中,T是温度,µ 和 κ 分别为动力黏度和导热系数.将式 (4) 代入式(1)～式(3) 即得到一维可压Navier-Stokes 方程组.

考虑如图1 所示的坐标系固定在波面上的驻定正激波情况,由于流动是定常的,时间偏导项为0,所有参数仅与位置x有关,式(1)～式(3) 可简化为常微分方程,将所得常微分方程对x积分,得到一维定常正激波内部流动的控制方程

式中,Jm,JP和JE是积分常数,不随x改变.实际上它们分别是质量、动量和能量的通量,并且因为在波前波后无穷远处 τ 和q趋于0,根据式(5)～式(7)显然有

下标1 和2 分别表示上游和下游无穷远处的状态.

进一步引入理想气体状态方程

其中R为气体常数.式(4)～式(7) 与式(11) 组成一个封闭的微分方程组,从中消去 τ,q,p以及 ρ 可得

该系统满足以下渐近边界条件

其中V2和T2可根据来流参数V1,T1以及正激波前后参数的关系求出.

2 打靶法

式 (12)～式(13) 构成一个封闭的一阶常微分方程组,沿x对其积分可获得V和T的沿程分布.然而边界条件 (14) 给在无穷远处,无法在数值积分中直接使用,人们很自然地想到将无限域截断成有限的计算域,再用打靶法[6,30,38]将边值问题转化为初值问题求解,其具体步骤如下.

(1) 选一略小于V1的速度值,将其对应的位置作为积分的起点,坐标指定为x=0,该速度值记为V(0).

(2) 选一略大于T1的温度值作为T(0) 的猜测值.

(3) 以V(0) 和T(0) 为初值,对式(12)～式(13) 构成的系统进行数值积分,积分沿着流动的方向进行,向下游推进足够远后截断.

(4) 若第3 步积分出的V和T的曲线像图1 给出的分布曲线那样逐渐趋平,则结束计算并输出结果;若积分过程中发散或出现参数明显远离波后值的趋势,则中断计算、修正T(0) 的值并回到第3 步.

第1 步的目的是对计算域的上游边界进行人工截断;第3 步可采用各种数值积分算法,包括但不限于Runge-Kutta 法;第4 步的目的是对与V(0) 相对应的初值T(0) 进行迭代修正.不同文献使用的算法在细节上可能有所差异,例如文献 [6,30,38] 直接求解的变量是p和 ρ,积分的起始点则选在p=(p1+p2)/2处,并从该点分别向上下游推进,这些文献的方法与本节的方法没有本质区别.

2.1 打靶法的失效

用前述的打靶法对单原子气体中Ma1=2 的正激波结构进行了求解,数值积分采用自适应步长的4 级Runge-Kutta 法,计算结果如图2 所示.结果表明,如果计算程序输入的参数组合经过精细调整,打靶法能给出看似合理的解,然而这依赖于对T(0) 值和下游截断位置极其谨慎的选择和微调,T(0) 稍有改变或计算域继续向下游延伸就可能导致发散.具体到图2 的算例,在大约x/λ1=12 (λ 表示分子平均自由程) 的位置之前,随着计算的推进,各参数快速趋近于波后值,并在波后值附近保持了一段距离的“稳定”,此时参数变化的趋势以及参数分布曲线的形态是正常的,在该位置终止计算可以得到看似收敛的解.但如果继续向下游积分,原本将要趋平的参数会快速偏离波后值,出现发散.

图2 打靶法计算的单原子气体激波结构 (Ma1=2)Fig.2 Profile of a Ma1=2 shock in a monoatomic gas calculated with the shooting method

合格的算法不应过于依赖人工介入,然而在激波结构问题中,打靶法很难像求解一般的两点边值问题那样使用二分法、割线法之类的迭代方法自动修正初值,因为第4 步极其依赖人工判断.此外,自动修正T(0) 需要第3 步中有一个明确的截断标准,但实际上计算域的长度很难确定,因为在计算中发现,无论T(0) 如何接近其准确值,只要向下游推进足够远,发散总是不可避免;过短的计算域无法保证覆盖整个激波结构,过长的计算域则导致初值必须非常准确才能保证收敛,因此对下游边界的截断和T(0)的迭代都依赖手动操作.高度的敏感性暗示着不稳定性的存在,需要在深入分析系统的动力学性质的基础上提出改进算法.

2.2 积分方向的重要性

激波内任意位置气体的状态可以表示为相空间中的一个点,那么激波压缩过程可以表示为相空间中一条光滑轨线.分别记波前和波后状态点为1 点和2 点,求解激波结构实际就是要找到一条连接1,2 两点的积分曲线,将V-T相平面中这样一条曲线叫做激波过程线 (即图3 中的曲线S).已知上下游无穷远处速度和温度梯度趋于0,即1,2 两点使式 (12)～式(13) 的右侧均为0,因此它们都是奇点,激波过程线的斜率 dT/dV在这两点处是 0/0 型未定式.

图3 V-T 平面中的相轨迹Fig.3 Phase portrait in the V-T plane

约定x增加的方向/流动方向/压缩进行的方向为激波过程的正向.为了更直观地了解系统的动力学特性,根据式 (12)～式(13) 所确定的方向场/斜率场,画出V-T平面中的相轨迹图,如图3 所示,其中箭头表示正向.相轨迹即图中带箭头的实线,是随着x的增加、相点的运动轨迹,每一条相轨迹都代表系统一个可能的解;相轨迹上每一点处的斜率与斜率场给定的一致,即相点的运动方向由方向场决定.相轨迹图非常有用,它为微分方程建立了一种类似于流线图的图案,揭示了系统的定性性质,通过它可以对解的行为形成图形化的认识.

从图3 可以看出,1 点是不稳定结点,2 点是鞍点,激波过程线S是从1 点发出并最终抵达2 点的异宿轨.在相平面中创建曲线M和N分别代表µ=0和 κ=0 情况下的激波过程线.在V-T相平面中,曲线M的方程可通过令式 (12) 右侧中括号内的项等于0 获得;同理曲线N的方程也可通过令式 (13) 右侧中括号内的项等于0 获得.由各自的曲线方程易知,M是抛物线;如果气体是定比热的,N也是抛物线.将M和N围成的区域命名为L.文献 [39] 证明了 ρ-p图中的激波过程线位于曲线M和N之间;类似地,很容易证明V-T图中的激波过程线S也落在M和N之间,即区域L中.

如图3 所示,从区域L内任意一点出发进行正向积分,除非该点严格位于所要求解的目标轨线S上,否则相点会先靠近2 点,但随后又迅速远离,始终无法抵达2 点,这也与图2 算例发散的具体表现(流动参数先趋近波后值,随后又迅速偏离) 吻合.实际计算中,即使初值点刚好落在目标轨线上,不可避免的舍入误差也会导致积分曲线偏离目标轨线,并被导出L区域.波后点是鞍点,这是所有正向求解随着数值积分向x=+∞ 推进最终都会发散的根本原因.有理由认为,所有包含正向推进的激波结构计算方法所获得的“收敛”结果,都是在发散之前的“合适”位置强行终止计算所产生的假象,这种人为截断既不合理,也具有很大的随意性.

3 逆向推进计算方法

前已述及,正向推进必然导致数值积分曲线偏离激波过程线、相点偏离2 点.这个问题可用逆向推进解决,即从区域L内、2 点附近选择一初值点,从该点出发向上游积分至流动参数几乎不再变化为止.图3 显示逆向积分曲线必然汇聚于1 点,所以该求解思路是可行的,但仍存在一个问题: 对初值不敏感固然是优点,却也导致逆向求解无法与打靶法搭配使用,因为打靶法靠迭代确定合适的初值,从而逼近精确解,而逆向求解会使初值误差衰减、相点无限趋近于1 点,无法为初值的修正提供反馈信息.因此需要一种新的初值确定方法.

3.1 基于L’Hôpital 法则的初值确定方法

虽然2 点的参数无法直接作为初值使用,但如果已知激波过程线的终点斜率 (dT/dV)2,就能以较高的精度在2 点附近确定一初值点 (V0,T0),其中V0的值由人为指定且略大于V2,T0可按Euler 格式给出

将式 (13) 和式 (12) 两端分别相除,整理后得激波过程线在V-T图中斜率的表达式

式 中JPm和JEm是常数:JPm≡JP/Jm,JEm≡JE/Jm.注意理想气体的比内能u是且仅是温度T的函数,且du=cvdT,其中cv是定容比热;µ 和 κ 一般也是温度T的函数.

2.2 节提到,1,2 两点是奇点,直接将波前或波后参数代入式 (16) 会使f/g成为 0/0 型未定式,给两端点处激波过程线斜率的确定带来困难,这个问题可用L’Hôpital 法则解决.首先,激波过程线S位于曲线M和N之间,所以其端点斜率也介于M和N的端点斜率之间;又因为曲线M和N为或者近似为抛物线,它们的端点斜率为有限值,所以S的端点斜率存在且为有限值,进而未定式f/g的极限也存在且为有限值,因此L’Hôpital 法则适用.

将式 (16) 中的 ζ 乘到等号左边,然后对等式两端同时取极限——沿激波过程线向端点i(i=1,2)逼近,在这个过程中T可以视为V的函数 (反之亦然).将f/g的分子分母同时对V求全导,其极限的值保持不变

其中下标V和T表示偏导,fV=JPm-V,fT=cv,gV=2V-JPm,gT=R.值得一提的是,这里允许气体是变比热的,即不要求cv是常数,它完全可以是T的函数.

式 (20) 可整理成关于端点斜率 (dT/dV)i的一元二次方程

其具有一对相异实根

其中 γ 为比热比,Pr为Prandtl 数.不难看出这对根的符号相反.1 点是不稳定结点,有无数对积分曲线从该点发出,其中一对在1 点处的斜率取正根;包括激波过程线在内的其他积分曲线在1 点处相切,斜率取负根.2 点是鞍点,从该点发出和抵达该点的积分曲线各有一对,两个根分别为这两对积分曲线在2 点处的斜率,其中负根是所需要的激波过程线的终点斜率.

3.2 一般求解程序

综上,激波内部结构逆向求解的一般程序可总结为以下4 步:

(1) 根据波前参数计算出波后参数,尤其是波后马赫数Ma2、波后速度V2以及波后温度T2;

(2) 将波后参数代入式 (22) 中负根的表达式,获得激波过程线的终点斜率 (dT/dV)2;

(3) 选一略大于V2的速度值V0,将其与其他所需参数代入式 (15),计算出相应的T0;

(4) 以 (V0,T0) 为初值点,使用负步长 (即逆流向推进) 对式(12)～式(13) 构成的系统进行数值积分;当V与V1、T与T1的差值满足精度要求时结束计算.

通过这4 步可获得V和T的沿程分布,其他参数可基于它们进一步计算获得而无需参与数值积分,例如将V和T的值代入式(5) 和式(11) 可求出ρ和p的值.

为了在精度和效率之间取得最佳平衡,建议在激波结构问题中使用变步长的数值积分方法.这不仅是因为激波厚度随激波强度的变化很大 (从无限弱激波的无穷大厚度到强激波仅数倍分子平均自由程的厚度),还因为速度梯度和温度梯度在激波内部不同位置差异巨大.自适应步长可避免在参数变化剧烈的区域网格过疏或参数变化缓慢的地方网格过密.

3.3 算例验证

为验证逆向推进法的有效性,对单原子气体中Ma1=1.01～100的正激波进行了求解,数值积分方法为步长自适应的4 级Runge-Kutta 法.不失一般性地,假设气体是定比热的,且其输运系数随温度的变化规律符合硬球分子模型的描述,即 µ 和κ 正比于T1/2.图4 给出了代表性工况的沿程参数分布.

图4 逆向推进法计算的单原子气体激波结构Fig.4 Shock structures in a monatomic gas calculated with the backward marching method

图4 中展示的流动参数已经过归一化处理,x=0 的位置被统一规定在V=(V1+V2)/2 处,坐标无量纲化所采用的分子平均自由程 λ 由硬球模型给定[1]

进一步地,将打靶法 (shooting method,SM) 与逆向推进法 (backward marching method,BMM) 计算的单原子气体中Ma1=2 的正激波内部的速度与温度分布进行了对比,如图5 所示.

图5 打靶法与逆向推进法对激波结构的计算结果的对比Fig.5 Comparison between the shock profiles calculated with the shooting method and the backward marching method

以上结果表明,逆向推进法可以正确可靠地求解激波内部流动结构,明显优于打靶法.

4 结论

针对激波结构问题,根据一维定常正激波控制方程的方向场画出了激波内部流动的相轨迹图,基于其拓扑结构分析了系统的动力学性质,阐述了积分方向对求解该问题的重要性,并指出在激波结构问题中,采用正向推进的传统打靶法无法避免发散.为解决该问题,提出采用逆向推进计算方法,并给出与之配合使用的初值确定方法.在非常宽的马赫数范围 (1.01～100) 内测试了逆向推进法,结果表明,该方法可以正确有效地求解激波结构问题,且与传统的打靶法相比,其具有以下优势:

(1) 计算无条件收敛;

(2) 对数值计算中不可避免的初值误差、舍入误差等不敏感,且误差随着计算的进行不断降低;

(3) 无需对初值进行迭代,积分一次即获得最终解,计算量小.

应当指出,黏性应力与热流密度的表达式虽然会影响讨论的定量细节,但不会使激波系统内在的定性性质发生根本变化,因此本文的求解方法及主要结论的适用范围不限于只具有线性本构的经典Navier-Stokes 方程.原则上,只要本构关系不违背热力学第二定律,即等效黏度和等效导热系数在任意位置处都非负,逆向推进法可以向各种具有Navier-Stokes 形式守恒方程 (式(1)～式(3)) 的复杂流体系统推广,如改进的Navier-Stokes 方程、以Burnett 方程为代表的高阶流体力学方程和Grad 矩方程等.对多原子气体中考虑体积黏度或不透明气体中考虑辐射传热的激波结构问题,新方法亦适用.