高超声速电离绕流的数值模拟

2024-04-22

高超声速电离绕流的数值模拟（共7篇）

篇1：高超声速电离绕流的数值模拟

钝锥高超声速绕流空间分离形态的数值模拟及拓扑分析

本文通过数值模拟和理论分析,系统研究了钝锥高超声速有迎角绕流时垂直于物面的截面流态沿轴向的演变过程.得出了高超声速情况下,背风分离区内的流态在大多数情况下貌似对称,其实是非对称的.这一重要结论.考察了周向网格疏密程度对流动结构分辩率的影响,发现对分离流动的数值模拟,必须在分离区内布置足够密的网格才能给出正确的结果.

作者：国义军宗文刚 Guo Yijun Zong Wengang 作者单位：中国空气动力研究与发展中心,四川绵阳,621000 刊名：空气动力学学报 ISTIC EI PKU英文刊名：ACTA AERODYNAMICA SINICA 年，卷(期)： 18(1) 分类号：V211.3 关键词：分离流动拓扑分析结构稳定性数值模拟

篇2：高超声速电离绕流的数值模拟

超声速高超声速球锥绕流的边界层稳定性特点初探

本文对一个球锥组合体,在超声速来流及零迎角条件下,对其边界层的稳定性特点进行了分析并讨论了转捩预测问题.结果表明转捩总是最先出现在来流马赫数不太大,或来流雷诺数较高并且具有较大的.周向速度或周向波数的情况下.

作者：袁湘江周恒赵耕夫 Yuan Xiangjiang Zhou Heng Hao Gengfu 作者单位：天津大学力学系,天津,300072 刊名：空气动力学学报 ISTIC EI PKU英文刊名：ACTA AERODYNAMICA SINICA 年，卷(期)：1999 “”(1) 分类号：V211.1 关键词：超声速边界层钝头锥稳定性

篇3：柱体绕流的自适应数值模拟方法

绕流问题在工程实际中常可遇到,如风对各种建筑物的绕流,河水流过桥墩,各种飞行器的设计,海洋石油工程中的开采平台、钻杆、水下输油管道等。在工业设备中绕流现象更是经常发生,如各类管壳式换热器。因此掌握流体绕物体流动的特性对工程实际和工业设备的设计非常重要;长期以来一直是学者们研究的热点问题,其中尤其以绕圆柱体和方柱体的流动最为常见和重要。这不仅因为它在工程技术中应用最广,而且研究它也是了解其他各种柱状钝体绕流的基础。绕方柱混合对流的研究是一个涉及到钝物体绕流和传热两方面相结合的课题。在过去的几十年里,该问题引起了众多的应用数学家、流体动力学家和传热学研究人员的极大关注。从学术研究的角度来看,绕方柱混合对流的基础是冷态条件下绕方柱流动流场的确定,它涉及到钝体绕流流动与分离、尾涡形成与变化的规律等涡旋动力学基础理论问题[1]。从工程应用的角度来看,深入研究由尖缘矩形横截面柱体所受到的绕流流体气动力载荷及其变化激起的流致振动问题是海洋工程和风工程等应用工程的需要。研究圆柱绕流问题在工程实际中也具有很重要的意义。如水流对桥梁、海洋钻井平台支柱、海底输运管线、桩基码头等的作用中,风对塔建筑、化工塔设备、高空电缆等的作用中,都有重要的工程应用背景。因此,对圆柱绕流进行深入研究,了解其流动机理和水动力学规律,不仅具有理论意义,还具有明显的社会经济效益。

随着CFD的发展,有关复杂几何外形的流场分析计算已经成为人们极为关心的问题。而合理设计并生成高质量的网格是CFD计算的前提条件。目前,处理复杂几何外形的CFD网格类型主要有如下三种:贴体的结构网格、非网格和Cartesian网格。事实上,对于结构的Cartesian网格由于其在网格生成方面的简易、快速等优点,在CFD发展的初期得到了广泛的应用,然而,由于其在处理固壁表面边界问题上的复杂性与低效性,很快又被贴体曲线网格所替代。近来,非结构的Cartesian网格由于采用的是四叉树的简单数据结构,易于结合网格自适应技术,其又重新引起了人们对Cartesian网格的普遍兴趣。然而,如何合理的处理复杂的物面边界条件仍然是Cartesian网格技术的一个关键问题。因此本文基于四叉树数据结构[2,3],发展了一种普遍适用于二维外形的Cartesian网格生成方法,运用对任意网格的切割细分算法,实现了针对几何外形的网格自动生成,使得网格生成更具灵活性和适用性;同时采用无网格方法[4,5,6,7,8]解决了Cartesian网格难以处理的物面边界问题,避免了通常采用的复杂物面切割方法[9]。此方法直接利用最小二乘法获得通量变量的值并且能够较好的处理物面附近网格点的复杂分布,具有简单明了、适用范围广等优点。因此本文将Cartesian网格方法和无网格方法相结合,利用Cartesian网格法处理计算区域内部网格点,而将无网格法用于处理物面边界条件。空间格式采用有限体积Roe格式,时间格式采用修正的四步Runge-Kutta时间推进方法,并且结合网格自适应,求解了二维复杂流场的Navier-Stokes 方程,实现了流场的自适应算法,与通常所采用的结构网格和非结构网格法相比,此方法大大改善了空间网格的分布,使得数值计算的分辨率、模拟效果及计算效率有了明显的提高。并将所发展的方法用于计算绕方柱流动和绕圆柱流动,所得的结果与相关文献[6,10—14等]的结果进行对比,进而验证了此方法在处理复杂物面边界条件下流动的合理性和有效性,也证明了Cartesian网格所具有的计算效率高、自适应能力强等优点。从而为Cartesain网格技术的未来发展提供了一种新的思路。

1 数值计算方法

1.1 二维Navier-Stokes方程组

当不考虑外加热和彻体力的影响时,直角坐标系下的二维可压非定常Navier-Stokes方程组的守恒积分形式为:

$\frac{\partial}{\partial t} \oint_{S} Q d S + \int^{\partial S} F \cdot n d l = 0 (1)$

式(1)中,Q为守恒向量,F为矢通量,S表示面积区域,∂S为S的边界,n为边界的外法向量,并且

$\begin{array}{l} Q = (ρ ρ u ρ v ρ E)^{Τ} ‚ \\ F = F_{c} - F_{v} = (F_{c x} - F_{v x}, F_{c y} - F_{v y}) (2) \end{array}$

矢通量F分解成对流矢通量Fc和黏性矢通量Fv两部分,各矢通量的具体表达式如下。

$F_{c x} = (\begin{array}{l} ρ u \\ (ρ u^{2} + p) \\ ρ u v \\ ρ u E + u p \end{array}) ‚ F_{c y} = (\begin{array}{l} ρ v \\ ρ u v \\ (ρ v^{2} + p) \\ ρ v E + v p \end{array}) 。$

其中 p, ρ, u, v, E 分别为压力,密度,直角坐标系下的速度分量和单位质量气体的总能量。

$F_{v x} = (\begin{array}{l} 0 \\ τ_{x x} \\ τ_{y x} \\ Π_{x} \end{array}), F_{v x} = (\begin{array}{l} 0 \\ τ_{x y} \\ τ_{y y} \\ Π_{y} \end{array}) (3)$

式(3)中

$\begin{array}{l} Π_{x} = u τ_{x x} + v τ_{x y} - q_{x}; \\ Π_{y} = u τ_{y x} + v τ_{y y} - q_{y} (4) \end{array}$

黏性应力项分别为:

${\begin{cases} τ_{x x} = 2 μ u_{x} - \frac{2}{3} μ (u_{x} + v_{y}) \\ τ_{y y} = 2 μ v_{y} - \frac{2}{3} μ (u_{x} + v_{y}) \\ τ_{x y} = τ_{y x} = μ (u_{y} + v_{x}) \end{cases} (5)$

热流量与温度梯度的关系符合Fourier定律,即

${\begin{cases} q_{x} = - k \frac{\partial Τ}{\partial x} \\ q_{y} = - k \frac{\partial Τ}{\partial y} \end{cases} (6)$

对于理想气体,有状态方程

${\begin{cases} p = ρ R Τ \\ h = c_{p} Τ \end{cases} (7)$

单位质量气体的总能量为

$E = \frac{p}{(γ - 1) ρ} + \frac{u^{2} + v^{2}}{2} (8)$

动力黏性系数μ是温度和压力的函数,在层流状态下通过Sutherland公式计算,即

$\frac{μ}{μ_{0}} \approx (\frac{Τ}{Τ_{0}})^{1.5} (\frac{Τ_{0} + Τ_{s}}{Τ + Τ_{s}}) (9)$

式(9)中,T0=273.16,对于空气,有

$μ_{0} = 1.716 \times 10^{- 5} Ρ a \cdot s, Τ_{s} = 124 Κ 。$

对于各向同性流体,导热系数k无方向性,仅随温度和压力变化,通过引入Pr数来确定,即

$k = \frac{μ c_{p}}{Ρ r} = \frac{μ γ R}{(γ - 1) Ρ r} (10)$

对于空气,在层流状态下可取Pr=0.72, cp和R分别为质量定压热容和气体常数。对于空间的离散采取的是Roe格式

$\tilde{F}_{i + 1 / 2} = \frac{1}{2} [F (U_{l}) + F (U_{r})] - \frac{1}{2} L^{- 1} | \tilde{A} | L (U_{r} - U_{l}) (11)$

对于时间的离散方式采用的是四步Runge-Kutta时间推进格式。

本文主要针对低雷诺数流动情形,采用层流模型计算。

1.2 Cartesian网格

1.2.1 数据结构

数据结构,在网格的生成和流场的计算中扮演着重要的角色,依靠其对网格数据的组织和管理,方便了网格的生成和数据的转换。本文采用的是树状数据结构,在Cartesian网格中,采用的基本数据结构是四叉树,即由一个父网格结点生成四个子网格结点。

图1 中表示了在四叉树数据结构的框架下的网格生成过程。

图1 中左边一列,是网格的深度值(Depth),表示当前网格的大小与细分层数,对于“根”网格,其深度值为0;右边是网格生成图(Grids),表示从“根”网格出发的网格细分和生成的过程;处于中间的是树状数据结构图(Tree),表示在树结构下,网格的生成过程。显然,没有再生成子结点的网格单元即是计算所需要的初始网格单元。

1.2.2 网格的分类

为了更好的处理所有的网格,可将所有网格划分为三个部分并定义如下:与物面相交的为物面网格(Wall Grids),剩余的网格中,在物面里面的是固体网格(Solid Grids),在流场计算中,固体网格是需要删除的,而物面外面的是流场网格(Flow Grids),面上的离散点表示固体的表面(如图2所示)。

1.2.3 网格自适应

合理的网格分布对于提高计算的效率和精度至关重要。本文采取压力梯度或者密度梯度的绝对值作为自适应参量进行网格自适应划分。阈值可以由下式给定:

$V_{s} = \frac{C R^{(n_{1} - n_{0})}}{Ν} (12)$

式(12)中Vs是阈值,n1为网格加密的深度值,n0为初始网格加密的深度值,N为参数,通常取100,CR为控制参数,CR的取值将直接影响网格的加密过程。对于任意一个网格,当其自适应参量大于给定的阀值时,其被标定为细分网格,进行自适应加密处理。相反地,当自适应参量小于给定的阀值,同时,网格深度小于或等于相邻网格的加密深度时,其被标定为粗化处理。

1.3 边界条件的处理

边界条件分为远场边界条件和物面边界条件。对于远场边界条件采用的类似于结构网格的处理方法,即通常采用的基于一维Riemann不变量的特征分析[19]来确定无黏流动变量的值,无需再进行特殊处理。由于采用Cartesian网格所得到的物面网格为非贴体网格,所以如何合理处理物面边界条件是Cartesian网格的一个关键。本文中采用的是基于最小二乘法的无网格方法来处理物面边界条件。

1.4 无网格法处理物面边界条件

为避免采用切割网格法生成的物面网格的复杂性,本文采用无网格法[5,8]来处理物面网格。无网格方法是新一代的计算方法,这种方法计算空间导数不需要借助于事先划定的网格,从而避免了高维拉普拉斯网格法中网格缠结和扭曲等问题[7]。

无网格法处理物面网格的具体应用过程为:

(1)对于一个物面网格找到一个离其中心点距离最近的点P,并且这个点在组成物面的有向线段上,如图3所示。

(2)按照由近及远的顺序在物面网格附近找到8个其它网格,并且这8个网格均是流场网格,如图4所示。

(3)可以假设P点附近区域的物理参数分布为:

$f (x, y) = a_{1} + a_{2} x + a_{3} y + a_{4} x y (13)$

因而根据8个网格中心点的参数并结合物面边界条件可以构成关于a1,a2,a3,a4的8个方程。这是一个超定方程组可以通过最小二乘法求解,便可以得到a1,a2,a3,a4的值,再代入式(13),即可得到f(x,y)的具体表达式。

(4)由以求得的式(13)的具体表达式可得物面网格的各条边界的参数值,这些参数值间接地满足物面流动条件。

对于物面处的无滑移边界条件[5,8]表达如下:

$\frac{\partial ρ}{\partial η} = 0 ‚ V_{ξ} = 0 ‚ V_{η} = 0, \frac{\partial ρ}{\partial η} = 0 (14)$

式(14)中 $\frac{\partial}{\partial η}$ 代表物面外法向方向导数,Vξ为物面切向速度分量,Vη为物面法向速度分量,如图5 所示

对于任一个物面网格,首先找到离其中心最近的物面点P。P点的法向和切向单位矢量可以根据其所在的有向线段求出,若P点恰为两条线段的交点,则其法向和切向单位矢量的值取为两条线段法向和切向的平均值。设已求得的物面法向单位矢量为

$\vec{η} = n_{x} \vec{i} + n_{y} \vec{j} (15)$

切向单位矢量为

$\vec{ξ} = m_{x} \vec{i} + m_{y} \vec{j} (16)$

然后进行坐标变换,将在x,y坐标系下的参数变换到以P点切向ξ和法向η构成的局部坐标系下,即可得到式(13)在ξ-η坐标系下的表达式为

$f (ξ ‚ η) = b_{1} + b_{2} ξ + b_{3} η + b_{4} ξ η (17)$

对于密度ρ将其在P处写为式(17)的形式

$ρ (ξ ‚ η) = b_{1} + b_{2} ξ + b_{3} η + b_{4} ξ η (18)$

由边界条件式(14)可知

$\frac{\partial ρ}{\partial η} |_{(0, 0)} = b_{3} + b_{4} ξ |_{(0 ‚ 0)} = b_{3} = 0 (19)$

将式(19)代入式(18)得

$ρ (ξ ‚ η) = b_{1} + b_{2} ξ + b_{4} ξ η (20)$

由上面的讨论可以看出,加入了边界条件后,未知数个数从4个降为3个。将已找到的邻近8个网格的中心点的密度值和坐标值分别代入式(20)即可构成8个方程3个未知数的方程组如式(21)。

$ρ_{i} (ξ, η) = b_{1} + b_{2} ξ + b_{4} ξ_{i} η_{i} ； i = 1, \dots, 8 (21)$

这是一个超定方程组,通过最小二乘法即可求解。求得密度分布的表达式之后,再将物面网格的边界的坐标代入其中就可以得到物面网格边界的密度值,这些值将在后面的流场计算中用到。

对于压强P,其求解过程与上面的对于密度ρ的求解过程完全一种样。

对于速度Vξ,利用边界条件有

$\begin{array}{l} V_{ξ} = b_{1} + b_{2} ξ + b_{3} η + b_{4} ξ η (22) \\ V_{ξ} = | (0, 0) = b_{1} = 0 (23) \end{array}$

将式(23)代入式(22)得到减少一个未知数的新函数

$V_{ξ} (ξ ‚ η) = b_{2} ξ + b_{3} η + b_{4} ξ η (24)$

将8个网格中心点的切向速度值和坐标值代入式(24),得到含有3个未知数的8个方程组

$V_{ξ_{i}} = b_{2} ξ_{i} + b_{3} η_{i} + b_{4} ξ_{i} η_{i} ； i = 1, \dots, 8 (25)$

这个方程组利用上面介绍的同样的方法求解。对于速度Vη,采用和Vξ同样的方法求解。最后,将得到的速度Vξ和Vη转换成x-y坐标系下的速度u,v的值。

2 结果与讨论

2.1 绕方柱流动

图6为边长为1的方柱中心位于坐标为(10.5,12.5),计算区域为41×25的初始网格图,来流条件为Ma=0.2,雷诺数分别为Re=20和Re=40。自适应网格图7、流线图8、9、10的结果表明当雷诺数 Re<50 时流动是稳定的,对于 Re=20,无量纲时间 t=100 左右达到定常解,对于 Re=50,t=200 左右达到定常解。

图8、图9、图10 分别为Re=20,Re=40,Re=50的方柱绕流流线显示。由此三图可以看出方柱尾流为两个对称的附体涡,为定常解,而且随雷诺数的增加,尾流对称涡的长度也在变大。当雷诺数 50<Re<55 时,方柱的尾迹涡型逐渐由定常附体涡向非定常涡旋脱落过渡。图11和图12为Re=55和Re=60的流线图,图中方柱尾流为非定常周期性涡旋脱落。因此,方柱绕流定常流动过渡到非定常流动的临界雷诺数 Rec应该在 50 和 55 之间,这与 Kelkar 等[12]采用线性稳定性分析得到的临界雷诺数Rec=53相符合。

图13为来流Ma=0.2,Re=100的方柱绕流的自适应网格图,图14为Re=100的方柱非定常方柱绕流的涡量云图,图15、图16分别Re=100和Re=250的非定常方柱绕流流线图,它们的特点是方柱尾流由定常的对称涡向非定常的涡旋脱落发展,涡旋交替由方柱的两个角点形成并脱落,并向周期性卡门涡街过渡。从表1、表2可以看出,Re=100和250方柱绕流的各项计算结果与文献 [12]符合的

2.2 绕圆柱流动

图17 为来流条件为Ma=0.2,Re=300的流体绕单位圆柱流动的初始网格图,从自适应网格图18,流线图19及涡量云图20上可以看出,流动过程为非定常过程,沿着来流方向,圆柱尾流涡旋的形成与脱落呈现出极强的规律性,周期性的涡旋从圆柱尾部交替脱落,形成了严格的卡门涡街。并且从图19的压力云图上可以看出对于压力梯度变化剧烈的地方,在图18的自适应网格上进行了相应的自适应加密,较好地捕捉了压力的布局情况。表3为与文献及实验结果的具体数据对比,可以看出符合的较好。

图21、图22为来流条件为Ma=0.2,Re=45的流体绕单位圆柱流动的涡量云图和流线图,从所得结果可以看出此流动过程为定常的。并且图22为与文献[16]相应结果的对比,可以看出符合的较好。

通过以上算例和分析可以看出Cartesian网格能够准确地根据压力梯度或者密度梯度的变化情况进行网格自适应加密或者粗化,使用无网格方法计算得到的结果与文献中使用贴体结构网格计算的结果或实验结果符合良好。因此本文所发展的方法能够较好的模拟绕方柱和绕圆柱等的定常和非定常的低雷诺数流动过程。这也说明本文所发展的网格生成技术、流场求解方法和物面边界处理的方法是成功、准确的。

3 结论

本文成功发展了基于四叉树数据结构的自适应网格生成及流场求解器。所发展方法的特点总结如下:

1) 基于四叉树数据结构发展了Cartesian网格的网格生成并与网格自适应相结合,改善了网格的空间分布,既保证了数值计算的分辨率和模拟效果又提高了计算效率。

2) 采用无网格法用来处理物面边界条件,使Cartesian网格亦能较好地解决具有复杂几何外形的物面边界问题。为Cartesian网格未来技术的发展提供了一种新的思路。

3) 通过柱体绕流的算例验证了所发展算法的有效性和合理性,计算结果证明了Cartesian网格所具有的计算效率高、自适应能力强等独特优点,为以后处理绕具有复杂几何外形物体的流动提供了新手段。下一步工作将此法发展到用于三维流场的计算。

篇4：附属多圆柱绕流的数值模拟

附属多圆柱绕流的数值模拟

利用FLUENT软件对圆柱周围等分布置6个小直径圆柱的流动问题进行了数值模拟.通过数值模拟,着重研究了附属圆柱直径及主从圆柱之间相对距离的`变化对升力系数、拖曳力系数和涡脱落频率的影响作用.

作者：杨安梅吕林董国海滕斌 YANG An-mei LU Lin DONG Guo-hai TENG Bin 作者单位：大连理工大学,大连,116023刊名：中国海洋平台 ISTIC英文刊名：CHINA OFFSHORE PLATFORM年，卷(期)：24(6)分类号：O351.2关键词：FLUENT 圆柱绕流附属小圆柱升力系数阻力系数涡脱落频率

篇5：高超声速电离绕流的数值模拟

对高超声速圆球模型飞行流场进行数值模拟,分别采用空气完全气体模型、平衡气体模型以及热化学非平衡11组元气体模型求解非定常轴对称N-S方程组.使用有限差分时间相关法捕捉激波,得到了定常流场的解.差分方程隐式部分采用了LU-SGS方法以避免矩阵运算,对化学反应和振动能量源项采用预处理矩阵以解决刚性问题.由计算结果处理得到的.阴影图和干涉条纹图与再入物理弹道靶实验照片进行了对比分析,验证了实验中圆球飞行流场大部分区域接近于平衡状态.

作者：柳军乐嘉陵杨辉作者单位：柳军(国防科技大学航天与材料工程学院,湖南,长沙,410073)

乐嘉陵,杨辉(中国空气动力研究与发展中心,四川,绵阳,621000)

篇6：高超声速电离绕流的数值模拟

二维高超声速吸气式飞行器一体化流场数值模拟

高超声速吸气式飞行器多采用机体/发动机一体化设计,使用CFD方法,对二维高超声速吸气式飞行器机体/发动机一体化流场进行了数值模拟,在来流马赫数6条件下,分别研究了发动机关闭,发动机通流与发动机点火三种状态下飞行器受到的`气动力系数的变化,以评估该飞行器的纵向气动性能,并同时对三种工作状态下飞行器机体/发动机一体化流场的流动特征进行了分析.

作者：金亮吴先宇罗世彬王振国 JIN Liang WU Xianyu LUO Shibin WANG Zhenguo 作者单位：国防科学技术大学航天与材料工程学院,长沙,410073刊名：弹箭与制导学报 PKU英文刊名：JOURNAL OF PROJECTILES, ROCKETS, MISSILES AND GUIDANCE年，卷(期)：28(6)分类号：V235.213关键词：高超声速吸气式一体化数值模拟

篇7：高超声速电离绕流的数值模拟

关键词：高超声速,数值模拟,热流,温度梯度,网格准则

未来高超声速飞行器外形较为复杂,尖化前缘[1]及流动干扰特征尤为明显.为了实施低冗余度无烧蚀或微烧蚀热防护设计,必须准确把握其气动热特性.目前高超声速气动热预测多采用以可压缩边界层理论为基础的工程预测方法,该类方法针对简单外形大面积区域已发展较为成熟,结合无黏流动和轴对称比拟思想[2]的流线跟踪方法[3]更进一步拓展了其应用范围.然而对于复杂外形,尤其是存在流动干扰区域,工程方法的预测精度仍旧难以令人满意.

数值模拟方法对于复杂流动的热环境预测有着先天的优势,它通过求解统一的控制方程而得到流场内的全部信息,从而在理论上非常适合于复杂外形的计算.然而由于高超声速问题自身的复杂性以及数值方法尚不完善等原因,目前可应用于工程型号计算的高超声速气动热环境数值模拟技术还存在较多问题有待改进[4,5,6],其中热环境计算的网格尺度选择问题尤为突出[7,8].

自20世纪80年代末期,国内外大量研究者对于热流计算的网格问题展开了研究,研究的重点主要集中于法向网格尺度的选择.Klopfer等[7]指出,热流计算结果与物面附近网格尺寸有密切关系,并且对于钝头体给出了通过减小网格法向间距改善热流计算的方法.Hoffmann等[8]详细描述了高超音速黏性绕流的表面热流与计算网格的关系,并且给出了计算球锥外形时网格尺寸与来流马赫数和雷诺数的关系.Cheatwood等[9]在LAURA软件中加入了准一维网格自适应进程,以保证物面法向网格雷诺数在1的量级.李君哲[10]研究了多个物理量对物面第一层网格高度确定的影响.谢锦睿[11]通过对边界层相似解的分析,指出需要综合温度边界层与速度边界层的分布特征来确定法向网格尺度.程晓丽等[12]在理论分析和数值试验的基础上,提出了基于分子自由程的法向网格划分标准.

本文结合钝双锥试验模型研究网格尺度对热流计算的影响,提出一种基于当地温度梯度的法向网格划分准则,为工程型号具体应用提供指导.

1 数学模型及算法

在笛卡儿坐标系下,采用非定常可压缩雷诺平均Navier-Stokes方程

其中,Q为守恒变量:F,G,H为3个方向的无黏通量向量;Fv,Gv,Hv为3个方向的黏性通量向量.

控制方程采用有限体积法进行离散.为使格式达到2阶精度,选用minmod限制器[13]通过MUSCL方法进行插值,无黏通量选择Roe格式[14],黏性通量采用2阶中心格式离散,时间方向则采用隐式LU-SGS方法[15]推进.

2 法向网格尺度影响

这里采用钝双锥层流试验模型[16]进行法向网格尺度研究.对称面及壁面网格如图1所示,模型长L=122.24mm,头部曲率半径Rn=3.835mm,前锥半锥角12.84°,后锥半锥角7°.

来流条件为:马赫数Ma∞=9.86,攻角α=10°,单位雷诺数Re∞=1.730×106 m-1,来流温度T∞=51.66K,壁面温度Tw=300K.计算共采用7套计算网格,具体法向网格尺度见表1.

不同计算网格均进行了200 000步隐式迭代以保证流场的充分收敛.由于不同计算网格得到的流场基本一致,这里仅给出网格VII计算的前锥部分压力和热流分布云图,见图2.

图3给出了不同母线的热流计算结果.其中φ=0°为背风线,φ=90°为侧线,φ=180°为迎风线.可以看到,随着法向网格的加密,计算结果收敛于恒定值,且该恒定值与试验数据吻合较好,这说明了本文采用的数值方法具有较好的计算精度以及网格收敛性.

从上述计算结果不难看出:对于给定外形和来流参数,存在一最大法向网格尺度δmax.当计算网格法向尺度大于δmax时,热流计算结果将偏离最终的网格无关解.故计算网格应首先满足

按照上述网格准则,近壁法向网格应尽可能减小以获得更准确的计算结果.然而在实际工程应用中这一点难以实现,原因在于过小的网格尺度将严重降低数值模拟的收敛速度,以至于在可接受的时间范围内难以得到收敛的数值解.

图4给出了7套计算网格数值模拟的密度残差收敛曲线.可以明显看到,较粗计算网格很快实现了残差收敛,而随着计算网格的加密,残差收敛速度呈非线性降低.

本算例共采用200000步隐式迭代才实现了最密网格Ⅶ热流的完全收敛.而在实际工程应用中,由于计算模型的复杂性以及计算资源的限制,如此充足的计算时间往往难以保证,所以有必要对较密网格随迭代时间的热流收敛特性进行研究.

图5给出了20000迭代步时不同母线的热流计算结果比较.可以看到,除网格Ⅵ和网格Ⅶ外,其他网格均实现了完全收敛.同时,网格Ⅵ和网格Ⅶ在所有高热流驻点区域内也得到了较好的结果,仅在大面积低热流区域出现了一定的偏差,且这种偏差随着热流的降低而增大.

从以上的计算规律不难看出,在计算时间(迭代步数)限定的前提下,对于给定外形和来流参数,存在一最小法向网格尺度δmin.当计算网格法向尺度小于δmin时,计算热流难以实现充分收敛.故计算网格还应满足

总之,在气动热的具体计算中,法向网格过粗或过密均不利于准确热流结果的获得.

需要说明的是,虽然过粗和过密网格均难以实际应用,但其产生偏差的机理却完全不同.

粗网格由于网格尺度过大,其空间分辨率难以准确模拟壁面温度梯度,因此会产生严重的空间离散误差进而影响热流计算.这一点可以由不同流动区域中粗网格的热流计算表现来证明.在高热流区域,壁面温度梯度较大,因此数值离散中模拟壁面温度型的截断误差较大,最终导致空间离散误差过大,从而产生较为明显的计算偏差.而在低热流区域,壁面温度梯度较小,对于相同壁面网格,其模拟温度型的截断误差较小,从而降低了空间离散误差,保证了计算精度.

密网格由于网格尺度过小,严重限制了迭代时间步长,从而在限定计算时间的前提下难以得到热流收敛解.同时,黏性项的收敛速度依赖于当地速度梯度和温度梯度.因此在背风区、大面积等速度梯度和温度梯度较小的黏性显著区域,收敛速度会进一步恶化.这就解释了密网格在低热流区域收敛较慢的原因.

3 法向网格准则

在工程应用中,总是希望能够找到一种合适的网格尺度满足式(1)和式(2),以同时保证计算精度和计算效率.

首先需要确定的是,所选择的网格尺度应偏向于δmin还是δmax.由前文计算结果可以看到,较密网格(网格Ⅴ～Ⅶ)在热流收敛区域内对于空间精度的提高并不明显,而其对低热流区域的收敛速度却影响较大.与之相对的是,较粗网格(网格Ⅰ～Ⅳ)在全热流范围内都有着较好的收敛速度,而在大部分热流范围内都可以满足热流计算精度.鉴于上述两点考虑,本文主要给出网格准则δmax,δmin的确定可参考文献[12].

为寻找这种网格尺度,首先对壁面处的热流通量数值计算过程进行分析.壁面数值格式一般存在多种离散形式,这里仅采用一类较为常用的基于虚拟网格技术的壁面热流通量格式,具体形式见图6,其中T1为壁面第一层网格格心处温度,T-1为壁面虚网格格心处温度,Kw为壁面处气体热传导系数,Qw为壁面热流,△y为壁面第一层网格高度.

此种情况下,热流通量可表达为

其中

在冷壁情况下,由近壁处温度边界层特性可知,随△y的增大,T1值将迅速增大,采用二阶精度计算壁面温度梯度的误差也随之迅速增大.

为将温度梯度误差控制在合理范围,就需要对壁面距离作出限制.这里不直接针对网格距离给出准则,而是将网格距离限制转化为对内点温度的限制,即不允许T1值超过一定限度.而由于直接确定T1值思路仍然不够清晰,这里利用式(4)通过对T-1值定义计算准则来间接实现对T1值的限制.

对于该类壁面数值格式,根据工程应用经验,这里给出限制准则

基于这种准则可以很容易得到内点温度及法向距离限制,即

由式(7)确定的δmax在相同壁温下随当地热流的变化而不同,热流较高时对应临界网格尺度较小,热流较低时对应临界网格尺度较大,这一变化趋势与之前数值计算结果是吻合的.

下面以真实计算状态来确定具体网格尺度,并与计算结果进行比较.

首先考虑驻点处计算结果,壁面参数为

由式(7)确定的网格准则为

可以看到网格Ⅰ不满足求得的网格临界准则,网格Ⅱ基本符合,而其他网格均满足,这与计算结果基本一致.

再考虑弹体低热流区计算结果,壁面参数为