高超声速分离-再附流动峰值热流的模型理论

方 芳, 鲍 麟

(中国科学院大学, 北京 100049)

引 言

近年来, 近空间高超声速飞行器的研制受到航空航天领域的广泛关注[1-2]. 与传统的宇宙飞船、 航天飞机等具有钝头外形的再入式飞行器不同, 近空间飞行器需要在大气层内做长时间的巡航和机动飞行, 这就需要其具有高升阻比的气动外形和控制结构. 而这种非流线型的构型在高速飞行时将不可避免地产生强烈的激波-边界层干扰(shock-wave/boundary-layer interactions, SWBLI), 并常常伴随着流动的分离和再附. 由 SWBLI 引起的再附点附近的热流峰值可以达到无干扰时的 10～100 倍, 甚至是驻点热流的数倍之多[3-4]. 因此, 在高速飞行时, 飞行器将处于严峻而复杂的气动热环境中, 准确地预测峰值热流对飞行器的热防护设计尤为重要.

现有预测 SWBLI 引起的峰值热流的方法仍有不足. 一方面, 实验很难同时满足实际飞行时的焓、 Reynolds数等条件[5]. 另一方面, 尽管计算流体力学(CFD)已广泛用于飞行器设计, 但其对峰值热流的预测仍然存在非常大的不确定度[6-8]. 传统上工程主要依靠热流与压强的比拟关联式[9-11]间接预测峰值热流. 然而不同模型的关联式需要不同的拟合参数, 带有强烈的工程经验性质, 其物理本质也不明确. 因此, 对峰值热流的理论分析和有效预测显得很有必要.

在早期的研究中, 有两类理论方法可用于求解分离-再附流动的热流问题. 第一类是以Lees等[12]和 Nielsen等[13]为代表的动量积分方法， 它是 Karman-Pohlhausen 方法[14]的拓展和修正. Klineberg等[15]引入能量方程拓展了 Lees 的方法使其能够求解热流. Murphy[16]评估了上述几种方法, 发现这些方法只能预测弱激波干扰下的压强, 无法准确预测热流. 而且, 动量积分方法需要假设一个合理的速度剖面和温度剖面, 这是非常困难的. 另一类可求解分离-再附流动热流的方法是一种渐近展开方法, 即三层结构理论. Lighthill[17]最早提出两层结构的概念, 认为黏性只在边界层内靠近壁面的薄亚层起作用. 根据这一概念, Neiland[18]和Stewartson等[19]针对超声速分离流动， Messiter[20]针对不可压缩平板后缘附近流动， 分别独立地给出了三层结构理论. 随后, Brown等[21]将三层结构理论拓展到高超声速, Rizzetta等[22-23]引入了能量方程, 讨论了超声速压缩拐角流动的热流, Daniels[24-25]将三层结构理论类比到再附点附近. 尽管三层结构理论在多方面都得到了成功应用, 但理论本身还存在一些局限性, 例如: 引入极高Reynolds数假设(Re>108), 使得理论在较低Reynolds数下存在不足, 而且理论的控制方程是隐式椭圆偏微分方程, 其求解难度不亚于一般的 CFD 求解程序, 当出现显著的回流区时, 解可能会出现不稳定性, 难以解决大尺度分离的问题[26].

然而, 工程上更关心的近空间高超声速飞行器面临的流动Reynolds数较低(104～105), 且更关注大尺度分离引发的气动加热问题. 最近, 李邦明[27]针对再附点局部区域, 首次引入斜驻点流动模化再附点附近干扰区无黏流动, 虽然其最终仍须借助热流比拟关联式[9]来计算峰值热流, 但验证了构造这类非渐近的模型理论的可行性和有效性.

为了解决中低Reynolds数下高超声速大分离流动中的再附区气动加热问题, 本文结合三层结构理论的思想， 引入斜驻点流动模型, 用以模化再附点附近整个剪切层内部的流动, 构建了一个新的模型理论框架. 同时, 针对压缩拐角分离-再附流动模型, 本文给出了再附点区外缘以及斜驻点流动区参数的求解方案, 从而形成一套完整的求解压缩拐角流动峰值热流的方法. 本文仅考虑完全分离的定常二维层流流动, 并假设流体为量热完全气体.

1 模型理论的构建

前人的研究指出, 由 SWBLI 引起的不同类型的分离-再附流动并没有本质的区别[26]. 无论是激波入射边界层, 还是压缩拐角流动, 或者是正激波诱导的分离和再附, 均是由边界层在分离点脱出形成自由剪切流, 并撞向壁面形成再附流动, 其压力和热流分布也相似. 李邦明[27]认为所有的再附流动是相似的, 可以用统一的理论来描述. 根据这个观点, 本文将压缩拐角分离-再附流动分成两个部分: 包含热流峰值点的再附点局部区域和外部无黏流动. 首先建立再附点局部区域的模型理论.

1.1 再附点局部区域的建模

再附区流动的特征很大程度上取决于分离产生的自由剪切层的特性. 边界层在分离点脱出后形成的自由剪切层一般只受到惯性力和黏性力作用, 经过不断演化, 层内流动更接近于类 Couette 流. 这种剪切层携带均匀涡量, 撞击壁面形成再附流动, 非常符合斜驻点流动的特征, 即由驻点流(西门子流动)[28]加入均匀剪切分量得到.

图1(a) 是斜驻点流动的示意图, 平板壁面处于Descartes坐标系y=0 处, 流体处在y>0 区域,θ为斜驻点的倾斜角, 即零流线与平板的夹角. 流动在远场的流函数可写为[29-31]

(1)

即西门子流动和均匀剪切流的叠加. 易知

(2)

其中,a表示x=0处的流向速度梯度,b表示涡量,a，b为常量. 由于均匀剪切流并非边界层型的流动, 斜驻点流动形成的边界层厚度可以用相应的西门子流动的边界层厚度表征, 即为

(3)

由经典斜驻点理论 (1) 式看出, 斜驻点流动速度在远场趋于无穷. 在实际流动中, 远场应为均匀流. 因此, 该模型需要修正. 根据三层结构理论[19], 流动在分离点或再附点等奇性点附近具有三层结构, 即有黏有旋的下层、 无黏有旋的主层和无黏无旋的上层, 见图 1(b). 这意味着, 附着后的剪切层的绝大部分可以当作无黏流来处理, 黏性只在剪切层靠近壁面的薄亚层起作用, 而剪切层的外部是无黏势流. 对于紧靠壁面的黏性层, 实际上是剪切层再附之后发展出来的新的边界层.

(a) Sketch of oblique stagnation-point flow

因此, 可以建立描述再附点下游局部区域的“斜驻点-三层结构”模型, 如图 2 所示. 靠近壁面的下层是有黏有旋的斜驻点边界层, 简称为“黏性层”, 其厚度可设为cδδ0, 其中cδ为常数; 中间层是无黏有旋的斜驻点边界层外流, 可称之为“无黏剪切层”, 其厚度为δe; 上层是无黏无旋的势流, 决定了无黏剪切层外缘的流动参数, 如压强pe, 温度Te, 密度ρe, 速度ue等.

图2 斜驻点-三层结构模型Fig. 2 Oblique-stagnation-point-triple-deck model

斜驻点-三层结构模型反映了再附点下游局部区域的流场特征, 同样也反映了其传热特征. 高速流动的能量方程[32]为

图3 剪切层传热示意图Fig. 3 Sketch of heat transfer in shear layer

1.2 斜驻点-三层结构模型的数学描述

基于再附点下游局部区域流动物理特性, 建立了“斜驻点-三层结构”模型, 下面将给出模型中黏性层和无黏剪切层即斜驻点部分的数学描述.

黏性层的厚度相对于流向尺度而言是一个小量, 因此可仿照经典边界层的量级分析方法, 给出黏性层的控制方程

(4)

其中， 下标 e 表示黏性层外缘即无黏剪切层中的物理量. 与经典的边界层方程不同的是, 斜驻点边界层的外流是带有涡量的无黏剪切流, 满足

对应动量和能量方程各增加一项.

鉴于斜驻点黏性流动方程 (4) 与经典边界层方程类似, 下面将拓展自相似理论, 给出方程 (4) 的数学求解方案. 第一, 类似 Howarth-Dorodnitsyn 变换[34-35], 引入坐标变换

(5)

并假设∂η/∂x=0, 其中Cρ=ρe/ρ. 由于再附点附近ρe,Te等均达到峰值, 其流向的变化可以忽略, 均看作常数, 因而δe也是常数. 第二, 参考 Tooke等[36]的不可压缩斜驻点流动解法, 假设黏性层内流函数仍然具有与无黏剪切层类似的驻点分量加剪切分量的形式, 即

则黏性层内的流动速度为

(6)

第三, 如上所述, 斜驻点边界层外缘并非均匀流, 再附点局部区域存在沿流向的温度梯度, 由速度沿流向线性分布可知, 能量方程的压缩项和耗散项沿流向呈二次函数分布, 因此对于壁温为Tw的等温壁面可假设温度分布也为二次函数分布, 即

(7)

从而可得无黏剪切层的流动速度为

(8)

边界层位移厚度的参数cδ可类比于平板边界层中的参考温度法, 取为

其中,c1=-1.27,c2=0.15,c3=1.57. 将 (6)～(8) 式代入到方程 (4) 中, 并按 (5) 式进行坐标变换, 最终得到

(9)

其中,

利用黏性层与无黏剪切层之间的匹配和壁面条件, 并假定ψ(ξ,0)=0, 可得方程的边界条件为

(10)

假设在再附点局部区域C仅是坐标η的函数, 即C=C(η), 则方程 (9) 和 (10) 与ξ无关, 即流动在再附点局部是自相似的. 然而此处再附点的自相似与经典的边界层自相似又有所不同. 可以看到, 再附点附近流动的动量方程化成了两个常微分方程 9(a) 和 9(b), 说明再附点附近流动是两个自相似流动即西门子流动和均匀剪切流叠加而成. 而能量传输是3个自相似过程的叠加.

要定解方程(9) 和 (10), 须给出黏性层外缘的流动条件, 从而确定方程中的无量纲参数Hr，Dr， ΔTr. 然而, 黏性层外缘为无黏剪切层, 其流动状态仍不易获得. 为了解决这个问题, 可以借鉴黏性激波层理论[37], 将控制方程描述的区域从黏性层扩展到整个剪切层. 再附点附近的剪切层一般由上游边界层分离脱出形成, 其厚度与上游边界层相当, 相对于流向尺度而言仍是一个小量. 因此, 方程 (9) 和 (10) 同样适用于整个剪切层. 而剪切层外缘的流动条件可以通过无黏理论给出, 较易获得, 具体过程将在下一节讨论. 此时有

最终, 对方程(9) 采用数值迭代法进行求解, 可给出无量纲壁面温度梯度Θ0′(0), 那么再附点后的峰值热流为

(11)

1.3 压缩拐角流动再附点局部的外缘条件

上节建立了斜驻点-三层结构模型的理论求解方案, 在求解高超声速分离-再附流动时, 仍须给出斜驻点流动外缘的无黏流动条件, 本节选取压缩拐角流动为例, 见图 4. 压缩拐角流动是一类典型的分离-再附流动, 图 4(a) 给出了二维中等角度压缩拐角模型的激波结构示意图. 流动再附发生在 3 区上游, 热流峰值位于再附点压缩波系的根部. 需要注意的是, 3 区可能很小, 以至于只有一个点. 为了应用前文给出的模型理论, 需要通过来流Mach数Ma∞, Reynolds数ReL, 温度T∞及壁温Tw, 拐角角度θw得到再附点附近的流动参数, 即由 0 区的参数求解 3 区的参数. 3 区的流动参数可分为两类: 第一类是无黏剪切层外缘的压强pe, 温度Te, 密度ρe, 速度ue, 不妨称之为“外缘参数”; 第二类是无黏剪切层的厚度δe以及a,b, 不妨称之为“剪切层参数”. 下面将分别说明如何得到这两类参数.

(a) Shock wave structure

“外缘参数”可以通过无黏理论和黏性干扰理论得到. 1 区的流动参数可由强/弱干扰理论[38]给出, 流动从 1 区到 2 区实际经过了分离激波-边界层干扰区, 假设分离点位置已知, 则可由自由干扰理论[39]给出 2 区流动参数; 2 区到 3 区须经过再附激波-边界层干扰区, 目前没有相关理论可以描述. 若并不关心干扰区内的流动变化, 可以将 2 区到 3 区之间的压缩波系简化成两道斜激波, 如图 4(b) 所示, 2 区和 3 区之间形成 6 区, 设 6 区的气流角度为

θ6=(θw-θ2)α+θ2， 0<α<1

(12)

则α决定了压缩波系的强度, 而α由来流参数和几何参数决定. 本文所使用的计算和实验数据表明, 在较广的流动参数范围内,α对流动参数的变化不敏感, 且有α=0.8. 下文将显示两道斜激波简化计算是合理的. 综上所述, 1～ 6 区的流动参数均可以通过无黏激波理论得到.

“剪切层参数”可以通过基本流动原理和近似得到. 具体地, 假设从流动分离到 3 区, 剪切层内的流量守恒, 近似有ρeueδe=ρuuuδu, 其中， 下标u表示上游干扰起点U处的参数. 则 3 区剪切层的厚度δe为

(13)

由式(8),ue=ax+bδe, 对于图 4 所示的压缩拐角流动, 剪切层撞向壁面的角度θ很小, 有a<

(14)

由式(2)可知a的值, 其中倾斜角θ=θw-θ2. 至此, “剪切层参数”均可以得到.

2 结果与讨论

2.1 理论结果与数值结果和实验结果的对比验证

为了验证前文的理论, 本文使用 DSMC 方法模拟了一系列的高超声速压缩拐角流动. DSMC 方法在中低Reynolds数高超声速流动中是公认可靠的, 已在相当多的研究当中得到应用[40-43]. 计算采用 Bird[44]开发的 DS2V 程序. 计算气体为空气(79%N2,21%O2)且采用量热完全气体模型. 分子碰撞采用变径硬球模型, 壁面采用完全漫反射模型. 根据不同算例的来流参数, 所用的模拟分子数约为 107～108. 计算终止条件为流场达到定常状态, 即计算结果不随时间变化. 达到定常时不同算例所对应的物理流动时间为 20～30 ms, 约为 30～50 倍的流动特征时间. 流动参数为Ma∞=7～11,ReL=(1,2)×104,θw=25°,T∞=60 K,Tw=300 K. 另外, 理论得到的峰值热流也与前人的 CFD 结果[45]、 实验结果[46-48]进行了对比. 各算例的流动参数见表 1.

表1 各算例的流动参数、 峰值热流及其对应的理论结果Table 1 Flow parameters, peak heat flux coefficients and theoretical results

针对不同的实验条件和数值模拟条件, 利用本文的理论框架预测其峰值热流系数St=qw/(ρ∞u∞cp(T0-Tw)), 其中，T0为来流总温. 理论预测的结果与数值模拟、 实验结果的对比情况如图 5 所示, 具体数值见表 1. 图中横坐标为数值模拟或实验结果, 纵坐标为理论预测结果; 实心点表示实验条件, 空心点表示数值模拟条件. 从图中可以看出, 理论预测与数值结果符合得较好, 说明本文的理论在中低Reynolds数的条件下是定量合理的. 相对于实验结果, 理论预测偏高 30% 左右. 这可能是由于实验中存在着真实气体效应和三维效应, 而理论中并未考虑这些因素的影响. 如 Nompelis等[49]指出, 振动非平衡效应可以使热流降低约 20%. 而三维效应显著减小分离区的大小[50-51], 这显然会导致热流的变化. 由下文的讨论可知, 减小分离区会降低峰值热流的大小.

图5 理论预测与数值模拟、 实验结果的对比 Fig. 5 Comparison of theoretical prediction with the results from simulations and experiments

2.2 与驻点边界层理论的关系

由前文的理论分析知, 再附点附近的流动可以看作是驻点流和均匀剪切流的叠加, 且已验证在流动倾斜度较大(即θ较小)的压缩拐角再附流动中是合理的. 本节将讨论另一种极限情况, 即流动接近垂直于壁面, 即正驻点流动.

当再附流动接近驻点流时, 剪切分量为零, 流动不带有涡量, 有b=0, 则Hr=0. 此时, 方程(9)(c) 和 (9)(d) 为齐次方程且具有齐次边界条件, 为零解. 方程 (9) 可简化为

这与驻点热流的结果[52]一致. 由此可见, 本文的模型理论可以退化到经典的驻点边界层理论, 因此, 在大倾斜度的再附流动中也应当具有一定的合理性.

2.3 分离区大小和边界层厚度对峰值热流的影响

控制分离区的大小和边界层的厚度是流动控制中的重要手段[53], 因而探讨分离区的大小和边界层的厚度对峰值热流的影响规律很有必要. 分离区的大小或边界层厚度的改变将影响再附点区的外缘参数和剪切层参数, 从而改变热流. 在本文的理论中, 可以直接设置分离点的位置或干扰起点U处的边界层厚度, 来研究其对峰值热流的影响.

图6 为不同的Mach数和Reynolds数条件下压缩拐角流动的峰值热流随着分离区大小的变化曲线. 其中， 横坐标Lu为无量纲的上游影响距离, 即图 4 中O点到U点的距离与平板长度之比, 表征分离区的大小; 纵坐标是通过理论得到的峰值热流系数. 当分离区极大时(Lu>0.8), 由于分离点太靠近前缘, 分离激波边界层干扰与前缘激波边界层干扰相互影响, 前文的压缩拐角求解方案不适用, 因此不予讨论. 一般情况下(Lu<0.8), 减小分离区的大小可以降低峰值热流, 且Mach数越低, 效果越显著.

图6 分离区大小对峰值热流的影响 Fig. 6 Influence of the separation length on the peaking heat flux

图 7 为不同的Mach数和Reynolds数条件下压缩拐角流动的峰值热流随着边界层分离时的厚度变化曲线. 其中, 横坐标为干扰起点的边界层厚度, 纵坐标是通过理论得到的峰值热流系数. 在本文所使用的算例中, 分离时的实际边界层厚度为 0.07～0.14. 从图中可以看出, 增大干扰起点的边界层厚度有利于降低峰值热流. 相对而言, 在较低Mach数和较低Reynolds数条件下, 效果更为显著.

图7 干扰起点的边界层厚度对峰值热流的影响 Fig. 7 Influence of the boundary layer thickness on the peaking heat flux

3 结论

本文提出了再附点附近流动的斜驻点-三层结构模型理论, 并据此给出了预测高超声速压缩拐角分离-再附流动的峰值热流的工程理论. 理论预测的压缩拐角再附热流峰值与 DSMC 计算结果大致相符, 验证了斜驻点流动模型理论在中低Reynolds数参数范围内的合理性和有效性. 并且, 理论可以退化至经典的驻点边界层理论, 表明模型理论具有可拓展性. 最后利用该理论探讨了压缩拐角流动中分离区的大小和干扰起点的边界层厚度对峰值热流的影响, 发现减小分离区的大小或增大干扰起点的边界层厚度均可以有效降低峰值热流, 尤其对于中低Reynolds数情况, 效果更为显著.

致谢作者得到中国科学院大学余永亮教授和王智慧教授的热心帮助和指导, 在此表示衷心的感谢.