一种新的玻尔兹曼方程可计算模型构造与分析1)
2021-11-10彭傲平李志辉吴俊林皮兴才蒋新宇
彭傲平 李志辉 吴俊林 皮兴才 蒋新宇
* (中国空气动力研究与发展中心超高速空气动力研究所,四川绵阳 621000)
† (北京航空航天大学国家计算流体力学实验室,北京 100191)
引言
近年来,随着对临近空间战略意义认识的加深,各航空航天大国发展了众多临近空间飞行器.为满足长航时、大载荷、高超声速巡航等需求,这类临近空间飞行器通常整体尺寸较大,各部件间差异显著,飞行过程会出现多物理多场耦合、多流区共存的多尺度复杂流动现象,传统的数值模拟方法如基于宏观连续介质假设的纳维−斯托克斯(N-S) 方程、模拟微观分子运动与碰撞的直接模拟蒙特卡洛(DSMC)方法难以满足多尺度一体化模拟需求,需要研究发展新的模拟手段[1-4].本文从介观气体分子速度分布函数满足的方程即玻尔兹曼方程[5-7]出发开展多尺度流动一体化模拟研究.
玻尔兹曼方程作为气体动理学理论的基本方程,能描述各流域气体流动输运现象,但它是一个高度非线性的积分−微分方程,其理论解和直接数值求解均难以实现[7].为利用玻尔兹曼方程能统一描述各流域气体流动的特性,众多学者提出了系列替代求解技术,其中应用最广泛的是模型方程方法,即用一个相对简单的表达式替换玻尔兹曼方程碰撞积分项[8],再对简化后的偏微分方程进行求解.
模型方程作为玻尔兹曼方程的近似,必须满足玻尔兹曼方程的基本属性,如总和不变量性质、H 定理、正确的输运系数等.第一个也是最著名的模型方程是巴特纳−格罗斯−克鲁克模型[8](BGK模型),但通过对其进行查普曼−恩斯科(Chapman−Enskog)展开分析得到的普朗特数(Pr)为1,而单原子气体的Pr约为2/3[9].为克服这一缺陷,国内外学者提出了一系列模型方程,包括椭球统计−BGK (ESBGK)模型[10-11]、沙克霍夫(Shakhov)模型[12-13]、刘(Liu)模型[14]、贝利(Belyi)模型[15]等,这些模型方程均假定其碰撞频率与分子速度无关.此外,也有学者发展了考虑碰撞频率与分子速度相关性的模型方程[16]以及适于多原子气体、考虑热力学非平衡的模型方程[17].
大量数学分析表明[16,18-20],上述模型方程均能满足玻尔兹曼方程的基本属性,且与原始玻尔兹曼方程的碰撞项具有类似的松弛特性,能保持与纳维−斯托克斯方程、直接模拟蒙特卡洛方法的收敛一致,因结构简单、贴近跨流域流体物理变化特点、可计算技术针对性强而得到广泛应用[21],如泰塔雷夫(Titarev)建立了物理空间和速度空间均采用非结构网格的内斯韦塔依(Nesvetay)求解器[22-23],针对返回式空间飞行器等高马赫数流动与直接模拟蒙特卡洛方法进行了系列对比分析.李志辉等[24-31]从研究描述各流域均适用的气体分子速度分布函数方程出发,吸收计算数学指数型积分求解原理,发展了气体分子运动论离散速度坐标法,研制经改进的高斯−埃尔米特积分方法等系列离散速度数值积分技术,消除气体分子速度分布函数对速度空间的连续依赖性,将玻尔兹曼模型方程化为在各个离散速度坐标点处具有非线性源项的双曲型守恒方程,基于计算流体力学数值求解技术,构造直接求解速度分布函数演化更新的气体动理论数值格式和大规模并行计算策略,建立了气体动理论统一算法(GKUA),开展了自由分子流到连续流复杂高超声速气动力/热问题大规模并行计算及其在航天工程中的应用研究,同时对跨流域多体干扰流动问题、跨流域非定常流动问题进行了研究,获得较好的成功.徐昆等[32-37]利用上述离散速度坐标法与BGK 类型格式相耦合,发展了统一气体动理学格式(UGKS),使之适于模拟稀薄流到连续流气体流动问题.
另一方面,从流动反映的物理现象来看,当流场中出现非平衡现象时,某点的速度分布函数已经不再是简单的麦克斯韦平衡态分布,而是呈现各种形态,如双峰分布、偏心分布等[38].这些分布与当地的气体分子速度、流场宏观速度、温度以及热流矢量、应力张量等宏观输运性质相关,这从玻尔兹曼方程的查普曼−恩斯科分析解[5,39]可以看出.因此在对玻尔兹曼方程碰撞积分进行物理分析、构造碰撞项的模型方程时应同时考虑这些因素,而常用的玻尔兹曼碰撞积分模型如沙克霍夫模型没有考虑应力张量的影响,ES-BGK 模型属于各向异性的类麦克斯韦分布,没有考虑热流矢量的影响.当模型方程不能全面反映这些因素时会导致模拟结果出现系统误差,例如在稀薄过渡流区流动的克努森层内黏性效应占主导,应考虑应力张量的影响,沙克霍夫模型与直接模拟蒙特卡洛方法模拟结果存在较大差别[40].
为此,本文将从玻尔兹曼方程出发,基于查普曼−恩斯科渐近分析解[9]的数学推导,构造一种同时考虑热流矢量和应力张量影响、满足玻尔兹曼方程高阶碰撞矩的综合型可计算碰撞模型,并在数学上分析其基本属性,建立与玻尔兹曼方程的联系,给出新模型与现有模型方程的关系,同时基于碰撞动力学确定该可计算模型碰撞松弛参数表达式.最后在气体动理论统一算法框架下,使用该新建玻尔兹曼方程可计算模型,通过模拟一维激波结构、二维近空间飞行环境平板和多体圆柱干扰流动来验证新模型的有效性和可靠性.
1 模型方程的构造
其中,i=x,y,z,Vx,Vy,Vz为气体分子速度分量;x,y,z为空间位置坐标,t为时间;S(f),(∂f/∂t)col为碰撞积分项.
在确定速度分布函数f后,可以通过对其取矩得到宏观流动参数,如密度 ρ,宏观速度分量Ux,Uy,Uz,温度T,应力张量分量Pij;热流矢量分量qx,qy,qz分别为
其中,R3为三维实数空间,dW=dVxdVydVz(下同),m为分子质量,n为分子数密度,kB为玻尔兹曼常数;i,j=x,y,z;C为分子热运动速度,Ci=Vi−Ui,etr为单位质量动量.此外流场压力p,黏性切应力张量分量 τij分别为
这里,符号 δij表示当i=j时为1,当i≠j时为0.
气体动理论统一算法所用玻尔兹曼方程可计算模型[24-31]可以写为
这里,ν 为气体分子碰撞松弛参数,fN为当地平衡态速度分布函数.常用的沙克霍夫模型[12]和贝利模型[15]的表达式分别为
这里fM为麦克斯韦速度分布函数,即
任何一个模型方程都必须满足玻尔兹曼方程碰撞积分的如下3 个基本属性[39].
(1)守恒性(质量、动量、能量守恒),即
(3)在平衡态时,有Sm(fM)=0.
此外,对于单原子气体,模型方程的查普曼−恩斯科分析得到的普朗特数Pr应约为2/3.
为研究构造新的玻尔兹曼方程可计算模型,首先分析其一阶查普曼−恩斯科分析解[9],其可以写为
显然,Φ (C)的第一项与常用的沙克霍夫模型[12]中第二项除系数外表达式一致,而 Φ (C)的第二项与贝利模型[15]的第二项除系数外表达式也一致.
从玻尔兹曼方程的查普曼−恩斯科分析解可以看出,流场中气体的速度分布函数与局部宏观流动参数如速度、热流、应力等相关,而沙克霍夫模型没有考虑应力的影响,贝利模型没有考虑热流的影响.当流场中局部区域的温度、速度梯度较大即热流和应力影响显著时,将出现较为严重的非平衡效应,仅考虑单一因素难以准确模拟复杂非平衡现象.因此,可以依据玻尔兹曼方程的查普曼−恩斯科分析解以及常用模型方程的表达式,构造如下同时考虑热流和应力影响的新型玻尔兹曼方程可计算模型
这里,A和B是与分子速度无关的常量,νN为碰撞松弛时间的倒数.显然,当A=1−Pr,B=0时对应沙克霍夫模型[12],A=0,B=(Pr−1)/Pr时对应贝利模 型[15].
2 新型玻尔兹曼模型方程的数学分析
在构造玻尔兹曼方程的一个新模型后,需要对其性质进行分析,验证其是否满足玻尔兹曼方程碰撞积分项的基本属性.
2.1 守恒性
由于
其中h=etr+p/ρ.上述3 个方程分别对应质量、动量和能量守恒方程.由此可知新模型方程满足玻尔兹曼方程碰撞积分项的守恒性.
在平衡态时热流和应力张量均为零,显然此时当地麦克斯韦分布是模型方程的解.
2.2 查普曼−恩斯科分析
为开展查普曼−恩斯科分析,首先将新模型方程的无量纲形式写成如下形式
其中 ξ 为与来流克努森数(Kn)同量级的小量,顶标 ∧表示无量纲量.按查普曼−恩斯科分析的流程[9],定义
进而黏性系数 µ、热传导系数 κ、普朗特数Pr为
可知对于沙克霍夫模型[12],A=1−Pr,B=0,νS=p/µ ;对于贝利模型[15],A=0,B=(Pr−1)/Pr,νVVB=Pr·p/µ,这与前述一致.通过选择合适的A,B可以得到正确的普朗特数Pr.
2.3 高阶矩分析
为了进一步验证新模型方程在数学上与玻尔兹曼方程的相容性,这里从麦克斯韦分子满足的玻尔兹曼方程高阶矩进行分析.由于麦克斯韦分子的特殊性,即其碰撞积分与相对速度无关,可以得到对应的高阶矩.
定义碰撞积分项的矩为
该结论与上一节一致,这进一步证明了新模型方程与玻尔兹曼方程的相容性,与沙克霍夫模型[12]和贝利模型[15]的一致性.
2.4 H 定理的证明
定义[39]
且从查普曼−恩斯科分析的过程可知分布函数f相对fM的扰动为一小量,因此有
即当上述不等式成立,且Kn较小时,新模型方程满足H 定理.从上面数学分析过程来看,νcoe在其取值范围内可以任意选取,但在实际物理过程中即分布函 数的碰撞松弛中 νcoe需要具有一定的物理意义.
2.5 碰撞松弛参数的确定
在气体动理学理论中,根据碰撞间隔理论[39]可以推导出碰撞时间间隔 η、黏性系数µ、压力p之间的关系为
显然,若将 η 直接视为玻尔兹曼模型方程中气体分子碰撞松弛参数 ν 的倒数,则上式与沙克霍夫模型一致,对应于 νcoe=1.但上述碰撞间隔没有考虑到分子碰撞后的速度残留[39],因此需要对其进行修正,通常将其取为扰动驰豫时间,即
表征为分布函数松弛至其初始值的 1/e所历经的时间(e 为自然对数的底),这里 λ 为平均自由程.
此外在直接模拟蒙特卡洛方法中,对于变径硬球(VHS)模型有[41]
上述 νcoe是通过微观分子碰撞间隔理论得到的,在气体动理学介观理论中,为了表征分布函数f松弛至的1/e 的物理意义,同时使 νcoe满足式(60),本文针对新模型方程取
对于单原子气体氩气(Ar),ω=0.81,此时νcoe=0 .819;对于空气,ω=0.77,νcoe=0.788 4.
3 数值分析
3.1 一维激波结构模拟分析
为验证本文新建模型的有效性,首先模拟了一维激波结构.设置激波马赫数为Mas=8,上游条件为:数密度n1=1 × 1023m−3,温度T1=200 K,速度V1=2105.738 m/s;下游条件:n2=3.821 × 1020m−3,T2=4174.414 K,V2=551.111 m/s;气体为氩气;计算中激波上下游取当地平衡分布,空间网格间距取上游气体分子平均自由程的0.05 倍.
图1 给出了采用沙克霍夫模型、本文新建模型、贝利模型模拟的激波轮廓线与直接模拟蒙特卡洛方法结果的对比,分别用“Shakhov”“new model”“VVB”和“DSMC”表示,其中直接模拟蒙特卡洛方法结果由贝德(Bird)的一维可视化直接模拟(DS1V)代码[41]计算得到,λ1表示激波上游气体子平均自由程,并将密度轮廓 (ρ−ρ1)/(ρ2−ρ1)=0.5 置于x/λ1=0 处.图2 给出了3 种模型模拟的一维激波内部无量纲热流和切应力变化曲线.从图1 可以看出本文新建模型能较好地反映Mas=8 的激波轮廓,验证了其有效性;对于几种模型而言,沙克霍夫模型与直接模拟蒙特卡洛方法结果符合最好,本文模型次之;差别主要在激波上游靠近x/λ1=0 区域,从图2 可以看出该区域无量纲切应力 τxx变化很小,而无量纲热流qx变化较大、且各模型结果区别也较大,仅考虑了热流项的沙克霍夫模型与直接模拟蒙特卡洛方法结果最为接近,说明在一维激波内部宏观切应力对内部参数分布影响很小,这可能是由于一维问题的宏观切应力仅出现在x方向且无量纲值相对热流较小,黏性效应不显著.此外模型方程模拟的温度轮廓在激波上游均出现了提前抬起的现象,这 可能是碰撞频率与分子速度无关造成的.
图1 不同模型模拟的激波轮廓与直接模拟蒙特卡洛方法对比Fig.1 Comparison of shock profiles simulated by GKUA with different models and DSMC
图2 不同模型模拟的无量纲热流和切应力对比Fig.2 Comparison of non-dimensional heat fluxes and shear stresses simulated by GKUA with different models
3.2 平板绕流模拟分析
为研究新模型方程模拟飞行器再入时在近空间环境黏性效应急剧增加的非平衡流场中黏性应力和热流影响较大的区域流动问题,首先考察平板绕流.通常对于这类问题,在克努森层内越靠近平板前缘,温度和速度梯度越大,非平衡效应越显著,对应的热流和黏性效应影响也越发严重.
来流条件:马赫数Ma∞=5,来流气体为氩气,来流温度和物面温度均为200 K,数密度为n∞=1.175 × 1020m−3.平板长度为1 m、宽0.04 m.计算中,为了揭示高马赫数附面层内流动与传热机制,直接模拟蒙特卡洛方法的背景网格间距取3.62 mm,统一算法贴近物面第一排网格间距为2.5 mm,均远小于来流分子平均自由程10 mm.统一算法模拟中除了所用新、旧碰撞模型外其他计算条件完全一样,物面采用与直接模拟蒙特卡洛方法一致的无穿透麦克斯韦散射分布,来流边界为来流无量纲麦克斯韦平衡分布,出口边界根据当地分子运动方向按特征边界条件处理.图3 给出了本文采用新模型方程模拟的温度分布云图(用GKUA_NewModel 表示的下半流场)与二维可视化直接模拟(DS2V)软件[42]模拟结果(用DSMC 表示的上半流场)的对比情况.从图中可以看出,两个结果几乎完全一致,而直接模拟蒙特卡洛方法尤其在尾部低速流场存在明显的统计波动.这说明新模型模拟的流场结果是可信的.
图3 本文模拟的温度云图与直接模拟蒙特卡洛方法结果对比Fig.3 Comparison of temperature contours simulated by GKUA with new model and DSMC
图4 给出了在x=0.1 m 截面上压力和温度沿y方向不同碰撞模型的统一算法模拟结果与直接模拟蒙特卡洛方法的对比,其中本文建立的新模型用“new model”表示,沙克霍夫模型[12]用“Shakhov”表示,贝利模型[15]用“VVB”表示.从图中可以看出,整体上使用新模型的统一算法结果与直接模拟蒙特卡洛方法模拟值相容,局部附面层区域y<0.15 m 时,相比沙克霍夫模型,新模型统一算法结果与直接模拟蒙特卡洛方法模拟值较为接近;而且还可看出对此黏性应力影响极为严重的附面层流动,使用贝利模型的统一算法结果与直接模拟蒙特卡洛方法模拟值符合更好;在0.15 m
图4 x=0.1 m 截面上压力和温度沿y 方向不同碰撞模型的统一算法计算值与直接模拟蒙特卡洛方法模拟结果的对比Fig.4 Comparison of pressure and temperature simulated by GKUA with different collision models and DSMC along the direction of y on the section of x=0.1 m
图5 和图6 分别给出了在x=0.5 m 和0.8 m 截面上压力和温度沿y方向不同碰撞模型的统一算法计算结果与直接模拟蒙特卡洛方法模拟值对比.随着x增大剪切层范围也扩大,在这些区域内流体黏性的影响较大,因而考虑了应力张量的碰撞模型,如使用3 种碰撞模型得到的统一算法计算分布变化趋势彼此一致,尤其是贝利模型与本文新提出的模型方程统一算法计算分布近乎完全重合,能较好模拟这些区域的局部流动,整体上都与直接模拟蒙特卡洛方法模拟结果吻合.而当y增大到剪切层外缘时,热流的影响开始显现加剧的某些区域,使用沙克霍夫模型得到的统一算法结果更吻合直接模拟蒙特卡洛方法模拟值.同样在当y进一步增大到远离附面层临近计算区域外缘,流动梯度仍然较大的近平衡态局部流场,几种模型得到的统一算法计算值均与直接模拟蒙特卡洛方法模拟值存在一定的系统误差,究其实质极可能缘于直接模拟蒙特卡洛方法完全使用麦克斯韦平衡态分布随机统计模拟且存在统计涨落,导致了两种方法体制的完全不一样,对于流场梯度较大的局部流动模拟,彼此结果间存在一定的计算偏差.
图5 x=0.5 m 截面上压力和温度沿y 方向不同碰撞模型统一算法计算结果与直接模拟蒙特卡洛方法模拟值比较Fig.5 Comparison of pressure and temperature simulated by GKUA with different collision models and DSMC along the direction of y on the section of x=0.5 m
图6 x=0.8 m 截面上压力和温度沿y 方向不同碰撞模型统一算法模拟结果与直接模拟蒙特卡洛方法的对比Fig.6 Comparison of pressure and temperature simulated by GKUA with different collision models and DSMC along the direction of y on the section of x=0.8 m
从图4～ 图6 可以看出,由于本文建立的玻尔兹曼方程新模型同时兼顾了流场中应力和热流对气体分子速度分布函数的影响,整体而言能较好地模拟全域流动.从统一算法计算x=0.1 m,0.5 m 和0.8 m 压力和温度变化来看,物面压力从1.5 Pa 降至0.9 Pa和0.8 Pa,温度从380 K 降至320 K 和300 K,相对于壁温200 K,物面温度跳跃逐渐降低,说明高马赫数绕流平板前缘强压缩流动在近空间飞行环境稀薄效应影响下逐渐变得稀疏,沿平板物面压力、温度不断减小,过渡到平板尾部低压常规流.
图7 给出了驻点线上不同碰撞模型得到的统一算法计算温度分布与直接模拟蒙特卡洛方法模拟值对比.平板头部压缩激波层是非平衡效应最为严重的区域,温度梯度和速度梯度均较大.从图中可以看出,在近空间飞行环境激波内部,本文建立的玻尔兹曼方程新模型统一算法计算结果始终与直接模拟蒙特卡洛方法模拟值符合较好.这进一步说明了需要同时考虑流场中黏性应力和热流对气体分子速度分布函数演化更新的影响.同时也看出基于玻尔兹曼方程可计算建模统一算法计算结果在驻点激波前部流动趋于近平衡态时温度高于直接模拟蒙特卡洛方法模拟值,这与一维激波结构中温度轮廓的变化规律类似.
图7 驻点线上不同碰撞模型统一算法计算的温度分布与直接模拟蒙特卡洛方法模拟值比较Fig.7 Comparison of temperature simulated by GKUA with different collision models and DSMC along the stagnation line
3.3 双圆柱干扰流动模拟分析
为验证新建模型统一算法对激波干扰的模拟能力,这里考察双圆柱干扰流动.来流条件和计算条件及边界条件与上节一致,圆柱直径为1 m,两圆柱圆心间距为4 m.在该条件下上下两个圆柱的头部脱体激波会相互作用,形成复杂的激波干扰流动.
图8 给出了分别使用3 种模型方程的统一算法计算温度分布云图与直接模拟蒙特卡洛方法模拟值的对比.两圆柱脱体激波相互作用后在x=−0.19m 的位置形成一个正激波类马赫杆,跨过马赫杆后温度和压力增加、速度降低.而马赫杆两端形成的斜激波仅对圆柱尾流产生影响,将导致单边圆柱法向力很小.从云图可以看出,使用贝利模型统一算法得到核心区温度分布与直接模拟蒙特卡洛方法模拟值最为接近,但马赫杆位置稍向后,而沙克霍夫模型差别最大,对马赫杆位置而言本文新建模型统一算法与直接模拟蒙特卡洛方法符合最好.表1 给出了3 种碰撞模型统一算法及直接模拟蒙特卡洛方法模拟得到的单边圆柱轴向力系数及法向力系数,以及不同模型所得轴向力系数与直接模拟蒙特卡洛方法模拟结果的偏差.从表中可以看出,本文新建立的碰撞模型统一算法计算结果与直接模拟蒙特卡洛方法模拟值最大偏差在4.6%以内,能较好地模拟气动力特性,精度优于常用的沙克霍夫模型.
图8 3 种模型方程计算得到的温度分布云图与直接模拟蒙特卡洛方法结果的对比Fig.8 Comparison of temperature contours simulated by GKUA with three collision models and DSMC
表1 不同模型/方法所得轴向力系数及其偏差与法向力系数Table 1 The axial and normal force coefficients by different models and DSMC method and errors between axial force coefficients
图9 给出了对称面(即y=0)上马赫杆附近沿x方向不同碰撞模型统一算法计算得到的压力、温度和速度分布与直接模拟蒙特卡洛方法模拟值的对比情况.本文新建碰撞模型较好地求解了二维绕流的激波位置及激波前后宏观参数的变化趋势,明显优于沙克霍夫模型,在激波内部模拟结果优于贝利模型.这是因为激波前后温度和速度梯度较大,但y=0 截面上仅存在x方向热流qx,而x和y方向的切应力 τxx和 τyy同时存在,速度梯度对分布函数的影响更大,使得本文新建碰撞模型能更好地模拟激波位置和内部参数变化.这与模拟一维激波结构所得结论不一致,可能是由于随着维度的增加多方向宏观切应力导致黏性效应更加显著,其对分子速度分布函数的影响增大,与热流共同作用影响激波位置 及其内部参数的分布.
图9 y=0 m 截面上压力、温度和速度沿 x 方向不同碰撞模型模拟结果与直接模拟蒙特卡洛方法的对比Fig.9 Comparison of pressure,temperature and velocity simulated by GKUA with different collision models and DSMC along the direction of x on the section of y=0 m
4 结论
本文通过分析玻尔兹曼方程的一阶查普曼−恩斯科分析解,结合沙克霍夫模型和贝利模型构造过程及数学推导与物理分析,建立了一种同时考虑热流矢量和应力张量影响、满足玻尔兹曼方程高阶碰撞矩的综合可计算碰撞模型,并在数学上分析了其守恒性、H 定理等基本属性,证明了新模型方程与玻尔兹曼方程的相容性,给出了新模型与现有模型的关系,通过碰撞动力学确定了碰撞松弛参数统一表达式.经过模拟一维激波结构、二维近空间飞行环境平板和多体圆柱干扰流动,并与直接模拟蒙特卡洛方法及不同碰撞模型方程的统一算法结果对比分析,得出以下结论:
(1)在稀薄效应显著、非平衡现象突出的流动区域,热流矢量和应力张量对当地气体分子速度分布函数的影响较大,在构造碰撞松弛模型时需要同时考虑两者的影响;
(2)在激波捕捉方面,流场中黏性效应显著的区域沙克霍夫模型结果与直接模拟蒙特卡洛方法模拟值差别较大,本文新建模型与贝利模型结果接近,在激波内部本文模型与直接模拟蒙特卡洛方法符合更好,而在黏性效应较小时仅考虑热流时结果与直接模拟蒙特卡洛方法符合更好,进一步说需要考虑多参数对分子速度分布函数的综合影响;
(3)本文新建模型的碰撞松弛参数中引入了与自然对数的底 e 相关的常数作为自由参数,间接表征了分布函数松弛的物理意义,虽然能满足模型方程的H 定理,但在流动的某些区域仍然与直接模拟蒙特卡洛方法存在差距.若能构造一个与当地非平衡参数相关联的碰撞松弛参数,将气体分子速度分布函数的松弛过程与当地流动参数的梯度耦合,应该能更好地模拟再入跨流域非平衡流动过程,这也是本文下一步拟开展的工作.
