李玉川, 鲍麟
中国科学院大学工程科学学院, 北京 100049
2018年5月9日 收稿; 2018年5月28日 收修改稿
基金项目: 国家自然科学基金（11572319）资助
通信作者: 鲍麟, E-mail:lbao@ucas.ac.cn

摘要: 微粗糙壁面上高速气流的流动和传热特性是当前关注的问题。以半无穷长正弦波状板作为微粗糙壁面模型，使用OpenFOAM程序，数值研究不同来流参数和不同波幅的正弦波状板流动。结果表明，在波状壁流动中，壁面热流系数C_Q与总阻力系数C_D呈现出直接关联，C_Q/C_D是表征粗糙壁流动与传热特性的关键参数。对于附着型流动，C_Q/C_D仅取决于粗糙度与边界层厚度之比，即等效粗糙度$\tilde h$；而流动产生分离后，C_Q/C_D不再随粗糙度变化，仅与波状板相似变化后的位置参量有关。
关键词: 波状粗糙壁高速流热流-阻力系数比
Numerical study of the relationship between heat flux and drag in the flow on wavy rough-wall
LI Yuchuan, BAO Lin
School of Engineering Science, University of Chinese Academy of Sciences, Beijing 100049, China


Abstract: This study focuses on the flow and heat transfer characteristics of high-speed airflow over the rough-wall which has periodic roughness and small amplitude. By using the OpenFOAM program, flows on the semi-infinite sinusoidal wavy plates under different inflow conditions are numerically studied. The results show that there is a direct correlation between the wall heat flux coefficient C_Q and the total drag coefficient C_D. For the attached flow, C_Q/C_D only depends on the equivalent roughness $\tilde h$, which is the ratio of the roughness to the thickness of the boundary layer. When the flow separation occurs, C_Q/C_D no longer changes with roughness and only depends on the position variable $\tilde x$, which is defined by the similar transformation.
Keywords: wavy rough wallhigh-speed flowratio of heat flux coefficient to drag coefficient
非光滑表面的高速流动传热规律既是一个基础问题，也是航天工程领域广泛关注和研究的课题。以再入式航天飞船为例，其再入行星大气层时速度可高达20~30倍声速，相应产生高达数千K的高温，为使飞行器免受高温损坏，需采用烧蚀热防护技术，即在其表面涂上反应潜热高的热防护材料，依靠高温烧蚀过程，带走大量热流，以保护飞行器。由于典型材料烧蚀过程是非均匀的，会形成粗糙壁面，研究其气动加热规律成为航天工程中的重要问题之一。从20世纪70年代以来，以美国PANT计划^[1]为代表，针对粗糙壁面热流开展了大量研究，一是利用实验数据拟合得到若干经验公式计算粗糙壁面的热流，如Dirling^[2]，Phinney^[3]，Grabow^[4]，二是用各种数值方法计算流场和粗糙壁面的热流，如Yang等^[5]，李俊红等^[6]。这些工作共同揭示影响粗糙壁面热流主要参数有雷诺数、边界层厚度以及等效粗糙元高度等，为工程设计提供了重要理论基础。
由于前人的粗糙壁研究源于传统的大钝头式再入飞行器，其波阻远大于摩擦阻力，因此这些研究大多只关心粗糙壁的热流，而极少有对摩阻的讨论。当前正在研制的下一代航天飞行器将具备大气层内高速巡航能力，多采用扁平化的乘波体等气动外形，其波阻大大减小，而表面摩阻所占比重极大地上升，最多可达30%以上^[7]，因此，飞行器壁面流动的摩擦阻力研究越来越受到重视。另一方面，下一代飞行器将实现可重复使用功能，工程设计提出采用非烧蚀热防护系统，使飞行器气动外形保持不变。实际上，飞行器表面在极高温作用下，不可能完全不烧蚀，尤其是多次长时间飞行后，其壁面应当看作存在微烧蚀，即是微粗糙的。综合上述两方面，微粗糙表面在高速流动中的壁面摩擦和传热是当前值得研究的问题，对新型飞行器的气动布局以及热防护系统设计有应用价值。
回顾光滑平板层流边界层理论可知，壁面热流与摩阻由临近壁面的边界层特性决定，且它们之间存在雷诺比拟关系^[8]，即壁面热流与摩阻是成正比的，沿着平板流动方向，其比值是一个常数，仅与来流马赫数Ma_∞、壁面温度T_w、气体的Pr数和比热比γ有关。对微粗糙壁而言，高速气体流经非光滑表面时，不仅壁面摩阻与热流发生变化，此时雷诺比拟是否仍成立存在疑问，而且与一般低速流动不同，压强波动产生的波阻是超声速时特有的现象，热流又与波阻有什么关系，是需要讨论的问题。
本文基于OpenFOAM程序数值模拟波状粗糙壁面在高速来流中壁面阻力与热流，及其流场特征，结果发现简单雷诺比拟关系不适用于粗糙壁情况，应代之以壁面热流与阻力之间的比拟关系，并分析该比拟关系的主要影响参数。为了深入分析该问题中的基本规律，同时针对新型飞行器高空巡航的特点，本文计算仅考虑层流情况，流动介质是空气，且看作量热完全气体(即特定气体常数R=287 J/(kg·K)、比热比γ=1.4、Pr=0.72，以下不再说明)，但结果揭示出的流动规律和机理也可为进一步研究各类工程实际情况提供参考。
1 物理模型与数值方法在实际问题中，粗糙壁面形状是极其复杂的，即使是微粗糙的壁面，也不易用较为简单的形状来描述。早期实验中为了对其进行刻画，通常采用在光滑表面上均布不同形状和尺寸的粗糙元来模拟粗糙壁面，而为了使这些不同形状大小和排布的粗糙壁面具有相同的度量特征，Nikuradse^[9]提出等效砂砾粗糙度：排布着普通粗糙元的壁面热流，和紧密排布着高度为k_s的半球(三维)或半圆(二维)粗糙元的表面的热流相同，那么k_s即为此粗糙壁面的等砂砾粗糙度。Signal和Danberg^[10]对普通粗糙元到等效砂砾粗糙度的映射关系做了改进，在工程中得到广泛应用。
等效砂砾粗糙度概念的提出启示我们，可以利用一种标准粗糙元模型等效地模拟实际粗糙壁面形状，以获得壁面热流变化的规律，那么不失一般性，本文选取正弦波状外形作为标准粗糙元模型。
图 1为半无穷长的正弦波状板流动模型示意图，其中V_∞、ρ_∞、p_∞、T_∞分别为水平来流速度、密度、压强和温度，沿水平和竖直方向分别建立x轴和y轴，以波状板的前缘为原点O，板的壁面波形方程为y=εsin(2πx/λ)，其中，λ和ε分别为正弦波形的波长和波幅。由于波状板是半无穷长，因此，选取波长λ为特征长度，则无量纲坐标x=x/λ，y=y/λ，壁面波形方程可重写为y=$\frac{h}{2}\sin (2{\rm{ \mathit{ π} }}\bar x)$，其中，$h = \frac{{2\varepsilon }}{\lambda }$为无量纲的波幅，即相对粗糙度，本文仅考虑h＜0.1的微粗糙波状板的情况。为分析壁面传热，设波状板的壁面温度为常量T_w，定义壁温比Θ=T_w/T_∞，文中计算一律取T_∞=T_w=300 K，即Θ=1。
Fig. 1

	Download: JPG larger image
图 1 正弦波状粗糙壁面示意图 Fig. 1 Schematic diagram of the sinusoidal rough wall 图 1 正弦波状粗糙壁面示意图 Fig. 1 Schematic diagram of the sinusoidal rough wall -->

选取来流密度ρ_∞、速度V_∞和温度T_∞分别为特征密度、特征速度和特征温度，可定义波状特征雷诺数为$R{e_\lambda } = \frac{{{\rho _\infty }{V_\infty }\lambda }}{{{\mu _\infty }}}$，来流马赫数为Ma_∞=$\frac{{{V_\infty }}}{{\sqrt {\gamma R{T_\infty }} }}$，其中，γ、R和μ_∞分别为空气的比热比、特定气体常数以及动力黏性系数。所以，来流状态由Re_λ和Ma_∞描述。若飞行器飞行高度大于40 km，马赫数取4~12，λ小于0.1 m，可估计Re_λ范围约为5×10³~1×10⁵。
本文将利用数值计算方法求解上述粗糙壁模型的流动问题，主要探讨壁面热流、摩阻与压阻之间的关系，为此，引入定义无量纲壁面压强系数C_p、摩阻系数C_f与热流率系数C_h ^[8]：

$\begin{array}{l}{C_{\rm{p}}} = \frac{{p - {p_\infty }}}{{\frac{1}{2}{\rho _\infty }V_\infty ^2}},\\{C_{\rm{f}}} = \frac{{{\tau _{\rm{w}}}}}{{\frac{1}{2}{\rho _\infty }V_\infty ^2}},\\{C_{\rm{h}}} = \frac{{ - {q_{\rm{w}}}}}{{{\rho _\infty }{V_\infty }\left( {{h_{{\rm{aw}}}} - {h_{\rm{w}}}} \right)}}.\end{array}$

(1)

已知高速气流的波状板流动的控制方程为二维定常、可压缩、黏性Navier-Stokes方程(以下简称N-S方程)：

$\begin{array}{l}\frac{{\partial \rho }}{{\partial t}} + \frac{{\partial \left( {\rho u} \right)}}{{\partial x}} + \frac{{\partial \left( {\rho v} \right)}}{{\partial y}} = 0,\\\frac{{\partial \left( {\rho u} \right)}}{{\partial t}} + \frac{{\partial \left( {\rho {u^2}} \right)}}{{\partial x}} + \frac{{\partial \left( {\rho uv} \right)}}{{\partial y}} = - \frac{{\partial p}}{{\partial x}} + \frac{{\partial {\tau _{xx}}}}{{\partial x}} + \frac{{\partial {\tau _{yx}}}}{{\partial y}},\\\frac{{\partial \left( {\rho v} \right)}}{{\partial t}} + \frac{{\partial \left( {\rho uv} \right)}}{{\partial x}} + \frac{{\partial \left( {\rho {v^2}} \right)}}{{\partial y}} = - \frac{{\partial p}}{{\partial y}} + \frac{{\partial {\tau _{xy}}}}{{\partial x}} + \frac{{\partial {\tau _{yy}}}}{{\partial y}},\\\frac{{\partial \left( {\rho E} \right)}}{{\partial t}} + \frac{{\partial \left( {\rho Eu} \right)}}{{\partial x}} + \frac{{\partial \left( {\rho Ev} \right)}}{{\partial y}} = \frac{{\partial \left( {pu} \right)}}{{\partial x}} + \frac{{\partial \left( {pv} \right)}}{{\partial y}} - \\\frac{{\partial {q_x}}}{{\partial x}} - \frac{{\partial {q_y}}}{{\partial y}} + \Phi .\end{array}$

(2)

式中：ρ、p、T分别表示流体密度、压强和温度；u、v分别为水平方向的速度分量和竖直方向的速度分量；E和e分别表示单位质量流体具有的总能和内能，$E = e + \frac{{\left( {{u^2} + {v^2}} \right)}}{2}$；流体内应力和热流率分别满足关系式：${\tau _{xx}} = \frac{4}{3}\mu \frac{{\partial u}}{{\partial x}} - \frac{2}{3}\mu \frac{{\partial v}}{{\partial y}}$，${\tau _{yy}} = \frac{4}{3}\mu \frac{{\partial v}}{{\partial y}} - \frac{2}{3}\mu \frac{{\partial u}}{{\partial x}}$, ${\tau _{xy}} = {\tau _{yx}} = \mu \frac{{\partial u}}{{\partial y}} + \mu \frac{{\partial v}}{{\partial x}}$, ${q_x} = - k\frac{{\partial T}}{{\partial x}}, {q_y} = - k\frac{{\partial T}}{{\partial y}}$, μ和k分别为动力黏度系数和导热系数；Φ为耗散函数，形式为$\Phi = \frac{{\partial \left( {u{\tau _{xx}} + v{\tau _{xy}}} \right)}}{{\partial x}} + \frac{{\partial \left( {u{\tau _{yx}} + v{\tau _{yy}}} \right)}}{{\partial y}}$。
有限体积法(finite volume methods，简称FVM)是计算流体力学中常用的数值方法，能够较好地保证流动的物理守恒性，因此，通常认为FVM类方法计算的壁面热流或摩阻的精度较高，本文利用OpenFOAM开源FVM程序，选取二阶Kurganov-Tadmor中心格式^[11]求解上述N-S方程。具体格式详见文献[11]，在此不再赘述。重要的是该格式满足TVD特性和最大值原理，同时模拟流场具有高解析度，且与时间步长大小无关，因此，是一种理想的高速流动的计算格式。此程序的具体算例验证将在下一节中给出。
2 数值程序的验证高速粗糙壁流动中压强梯度较大，这要求计算程序具有更高的流场解析度，为了验证本文采用的OpenFOAM程序的可靠性，首先，计算高马赫数光滑平板边界层前缘的黏性干扰区，并将壁面压强、摩阻与热流计算结果与已有的黏性干扰理论值作对比；其次，模拟超声速正弦波状板流动，在距离前缘充分远的区域，检验流场的数值解是否与理论上的二维定常超声速波形壁理论解一致；最后，对正弦波状板的高速流动计算结果的时间步长和网格密度无关性做检验。
首先，黏性干扰区算例的参数选取为Ma_∞=8，单位雷诺数Re_L=2×10⁵ m^-1。依据Hayes和Probstein^[12]得出的强干扰区摩阻和热流的分布公式为：

$\begin{array}{l}{C_{\rm{f}}} = {c_0}{{\bar \chi }^{1/2}}\sqrt {\frac{C}{{R{e_x}}}} ,\\{C_{\rm{h}}} = \frac{1}{2}{c_0}{{\bar \chi }^{1/2}}\sqrt {\frac{C}{{R{e_x}}}} ,\\p = {g_0}\bar \chi {p_\infty }.\end{array}$

(3)

式中：c₀≈0.35(T_w/T₀+4.51)，$\bar \chi = Ma_\infty ^3\sqrt {\frac{C}{{R{e_x}}}} $，g₀≈0.1(T_w/T₀+0.176), $C = \frac{{{\mu _{\rm{w}}}{T_\infty }}}{{{\mu _\infty }{T_{\rm{w}}}}}$。图 2选取算例中0.5＜x＜1的区间，对应6＜χ＜9，此时为强干扰区域，将热流率系数和摩阻系数与强干扰区理论值作对比。结果表明，C_f和C_h理论值和计算值基本一致，偏差5%以内，而C_p最大偏差为10%，即p的偏差4%，因此，程序计算有压力梯度流场的热流和切应力是较为可靠的。
Fig. 2

	Download: JPG larger image
图 2 平板黏性干扰区的Cf、Cp、Ch计算值与理论值的对比 Fig. 2 Comparison between the calculated and theoretical values of Cf, Cp, and Ch along a flat plate with viscous interference 图 2 平板黏性干扰区的Cf、Cp、Ch计算值与理论值的对比 Fig. 2 Comparison between the calculated and theoretical values of Cf, Cp, and Ch along a flat plate with viscous interference -->

其次，无穷长超声速定常波形壁的二维流动问题中的压力系数的理论值^[13]为C_p=$\frac{{4{\rm{ \mathit{ π} }}\varepsilon }}{{\lambda \sqrt {Ma_\infty ^2 - 1} }}\cos (2{\rm{ \mathit{ π} }}x/\lambda )$，当来流雷诺数较大，而x也较大时，本问题和它是极接近的。图 3给出Ma_∞=4，Re_λ=10⁵，h=0.02，15＜x＜19的压力系数的计算值和理论值对比，其最大偏差为7%，即压强p的偏差3%，因此，程序在计算波形壁面高速流动中的压强是较可靠的。
Fig. 3

	Download: JPG larger image
图 3 无穷长波状板二维超声速流动的压强分布的理论值与计算值对比 Fig. 3 Comparison between the calculated and theoretical values of the pressure distribution along a wavy plate in 2-D supersonic flow 图 3 无穷长波状板二维超声速流动的压强分布的理论值与计算值对比 Fig. 3 Comparison between the calculated and theoretical values of the pressure distribution along a wavy plate in 2-D supersonic flow -->

最后，选取不同的网格和时间步长计算一个波状板模型的流动，其参数为Ma_∞=4，Re_λ=2×10⁴，波状板的波幅(粗糙度) h=0.02，板长为20λ，计算显示结果具有良好的收敛性。图 4(a)和4(b)分别展示不同网格密度以及不同CFL数下的结果对比，其中CFL数$C = \frac{{{V_\infty }\Delta t}}{{\Delta x}}$，Δt为时间步长，Δx为最小网格尺度。为比较不同算例之间的细节差异，图中仅画出一个周期波形的壁面上无量纲热流率系数C_h、摩阻系数C_f和压强系数C_p。结果表明，当网格数大于500×200时，计算结果不随网格加密而变化；当CFL数小于0.01时，计算结果已趋于一致。因此在后续的计算中选取网格数为500×200，CFL数为0.01。
Fig. 4

	Download: JPG larger image
图 4 不同计算参数下的波状板计算Ch、CQ、Cp曲线 Fig. 4 The calculated curves of Ch, CQ, and Cp along the wavy plate with different calculation parameters 图 4 不同计算参数下的波状板计算Ch、CQ、Cp曲线 Fig. 4 The calculated curves of Ch, CQ, and Cp along the wavy plate with different calculation parameters -->

3 算例及分析3.1 沿波状壁面的热流-阻力关联关系在光滑壁面的平板层流边界层中，沿壁面的压强是不变的，但壁面热流和摩阻沿平板的变化规律相似^[14]：C_f∝Re_x^1/2，C_h∝Re_x^1/2，因此，存在雷诺比拟关系式${C_{\rm{h}}}/{C_{\rm{f}}} = \frac{1}{2}P{r^{ - 2/3}}$。对于粗糙壁面而言，由于壁面形状凹凸不平，流动随之压缩或膨胀，产生壁面压强波动，更重要的是，超声速流动中，这种压强波动在一个周期内的平均作用效果并不是零，而是对壁面产生了阻力作用，即波阻作用，这意味着高速气流受粗糙壁的作用包括摩擦阻力和压差阻力(或称波阻)两部分，这两种阻力都对其做功，进而影响气流的总能量以及传入壁面的热量。再考虑式(2)中能量方程右端显示，导热项同时受到压力做功以及耗散作用的联合影响，在粗糙壁面附近，不仅气流的摩擦剪切作用强，耗散作用显著，而且压差带来的阻力是不可忽略的，压力功也非常显著。由此推测，粗糙壁面的气动加热量应当与壁面总阻力存在关联，热流系数与摩阻系数的比拟关系有可能不再满足。
本小节通过计算结果分析波状板y=$\frac{h}{2}\sin (2{\rm{ \mathit{ π} }}\bar x)$的壁面热流-摩阻比及热流-(总)阻力比沿x方向分布是否存在特定规律，以此判断粗糙壁的气动加热究竟是与摩阻直接关联，还是与总阻力存在相关性。需要指出，由于流体的黏性(或热传导)效应响应时间比压力响应时间慢得多，所以，周期性波状板面上的热流、摩阻和压强分布虽然都具有周期性，但它们的相位不是完全同步的，为了分析热流与阻力之间的关联，我们选取一个波形周期内的物理量积分平均值作为研究参量，比如平均压阻为压强在一个周期壁面上积分的水平分量的平均值。因此，结合式(1)，分别定义单个正弦波周期内壁面上的平均摩阻系数C_Df、平均压阻系数C_Dp、平均阻力系数C_D和平均热流系数C_Q如下：

${C_{{\rm{Df}}}} = \frac{{\int_{{x_0} - \lambda /2}^{{x_0} + \lambda /2} \tau \cdot {\rm{d}}x}}{{\frac{1}{2}{\rho _\infty }V_\infty ^2 \cdot \lambda }} = \int_{{{\bar x}_0} - 1/2}^{{{\bar x}_0} + 1/2} {{C_{\rm{f}}}} {\rm{d}}\bar x,$

(4)

${C_{{\rm{Dp}}}} = \frac{{\int_{{x_0} - \lambda /2}^{{x_0} + \lambda /2} {p \cdot \tan \theta {\rm{d}}x} }}{{\frac{1}{2}{\rho _\infty }V_\infty ^2 \cdot \lambda }} = \int_{{{\bar x}_0} - 1/2}^{{{\bar u}_0} + 1/2} {{C_{\rm{p}}}\tan \theta {\rm{d}}\bar x} ,$

(5)

${C_{\rm{D}}} = {C_{{\rm{Df}}}} + {C_{{\rm{Dp}}}},$

(6)

${C_{\rm{Q}}} = \frac{{\int_{{x_0} - \lambda /2}^{{x_0} + \lambda /2} {q \cdot {\rm{d}}x/\cos \theta } }}{{{\rho _\infty }{V_\infty }\left( {{h_{{\rm{aw}}}} - {h_{\rm{w}}}} \right) \cdot \lambda }} = \int_{{{\bar x}_0} - 1/2}^{{{\bar x}_0} + 1/2} {\frac{{{C_{\rm{h}}}}}{{\cos \theta }}{\rm{d}}\bar x} .$

(7)

式中：θ为波状板表面倾角，$\theta = \frac{{{\rm{d}}\bar y}}{{{\rm{d}}\bar x}} = {\rm{ \mathit{ π} }}h\cos (2{\rm{ \mathit{ π} }}\bar x)$，x₀为正弦波形中心，x₀=x₀/λ。
计算来流Ma_∞=4和Re_λ=2×10⁴的波状板流动，图 5显示该算例中平均热流-摩阻系数比C_Q/C_Df和平均热流-阻力系数比C_Q/C_D随x的变化规律，x代表单个正弦波的中心位置。结果显示：对于较小的h(0.01≤h≤0.03)，C_Q/C_Df在板前端沿的变化非常显著，随后变化趋于平缓，而C_Q/C_D沿x基本不变，这说明热流与总阻力呈现出直接相关性，对于粗糙壁流动，分析热流-阻力系数比更有意义。从另一个角度看，光滑平板上的雷诺比拟可以看作此关联关系的一个特例，即压阻为零的特例。
Fig. 5

	Download: JPG larger image
图 5 h较小时CQ/CDf和CQ/CD随x的分布(Ma∞=4，Reλ=2×104) Fig. 5 Distributions of CQ/CDf and CQ/CD along x at the small h values (Ma∞=4, Reλ=2×104) 图 5 h较小时CQ/CDf和CQ/CD随x的分布(Ma∞=4，Reλ=2×104) Fig. 5 Distributions of CQ/CDf and CQ/CD along x at the small h values (Ma∞=4, Reλ=2×104) -->

图 6给出更大范围的粗糙度(0.01≤h≤0.1)下，C_Q/C_D随x的变化规律，可以看到，虽然h≥0.04后，C_Q/C_D随x有明显变化，但是，几条C_Q/C_D曲线绝大部分重合在一起，即平均热流-阻力系数比几乎不随h改变而改变，体现出另一种特殊的热流-阻力系数相关关系。
Fig. 6

	Download: JPG larger image
图 6 不同h下的CQ/CD随x的变化曲线(Ma∞=4，Reλ=2×104) Fig. 6 Distribution of CQ/CD along x at various h values(Ma∞=4, Reλ=2×104) 图 6 不同h下的CQ/CD随x的变化曲线(Ma∞=4，Reλ=2×104) Fig. 6 Distribution of CQ/CD along x at various h values(Ma∞=4, Reλ=2×104) -->

通过对比不同h时的流场结果，我们发现两种C_Q/C_D关联规律实际上对应着不同的流动结构。图 7(a)、7(b)分别给出h=0.03和h=0.06时波状壁面附近的流场结构。可以看到h=0.03时，流场没有出现分离现象，与光滑平板边界层的雷诺比拟关系类似，热流-阻力系数比C_Q/C_D几乎不随x变化；h=0.06时，由于壁面波幅较大，在波谷区域较强的逆压梯度使得壁面流动分离，形成明显的分离涡结构，此时，分离区内壁面热流和阻力分布与层流边界层结果存在很大差异，但是分离流场的结构是非常相似的，与h关系不大，对应的，流动分离后平均阻力系数等壁面参数应基本不变。
Fig. 7

	Download: JPG larger image
图 7 h=0.03和h=0.06时壁面附近的流场图 Fig. 7 The flow field diagram near the wall when h=0.03 and 0.06 图 7 h=0.03和h=0.06时壁面附近的流场图 Fig. 7 The flow field diagram near the wall when h=0.03 and 0.06 -->

上述计算结果表明，在波状板流动中，壁面热流与阻力呈现出明显的关联关系，平均热流-阻力系数比C_Q/C_D是衡量流动与传热特性的关键参数。在相对粗糙度h很小时，该比值近似与x无关，随h增大而减少；而对略大的粗糙度h，波形板面上的波谷内存在分离涡结构，对应的C_Q/C_D随x变化而近似与h无关。
计算中同样发现，当h很大或者x很大时，流场变得异常紊乱，壁面上出现多个涡聚集，壁面热流和阻力剧烈抖动，此时需要考虑非定常、转捩及湍流等复杂因素。相对于光滑壁面，粗糙壁面会使边界层提前转捩，PANT转捩准则^[15]是常用的粗糙壁面转捩准则：

$R{e_{\theta tr}} = 215{\left( {\frac{k}{\theta }\frac{{{T_{\rm{e}}}}}{{{T_{\rm{w}}}}}} \right)^{ - 0.7}},\left( {Ma \ne 1} \right).$

(8)

式中：Re_θtr是以动量边界层厚度θ作为特征尺度的雷诺数；k为粗糙元尺度，即为本文中2ε。定性上讲，波状板的粗糙元尺度越大，位置离前缘越远，流动越趋向转捩乃至湍流。按本文设定的来流条件，设最大板长为x=1 m，即Re_x约为10⁵~10⁶，则动量边界层厚度θ量级约为10^-3 m^[16]，代入式(8)得到临界转捩粗糙元高度约为k≈0.01 m，即h=0.1。因此，本文仅讨论h＜0.1时的波状板流动，将其看作层流流动，并且不包含出现局部紊乱流场时的情况，下面将集中探讨定常层流流动条件下的C_Q/C_D的具体变化规律。
3.2 热流-阻力关联关系的归一化波状模型的粗糙壁流动中，平均热流-阻力系数比C_Q/C_D不仅与x和h有关，还受Ma_∞、Re_λ影响。在不同的Ma_∞和Re_λ下，热流-阻力关联关系是否一致存在，是本小节讨论的问题。计算选取Ma_∞为4~14、Re_λ为5×10³~1×10⁵等共10组算例，每组算例中，相对粗糙度范围为0.01≤h≤0.1。
3.2.1 未分离的波状板流动选取微小粗糙度的算例进行讨论，即波形板流动未出现明显分离时，例如图 5中对应的h≤0.03的曲线，C_Q/C_D几乎不随x变化。再考虑到不同来流参数的影响，此时，热流-阻力关联式可表示为C_Q/C_D=f(h, Ma_∞, Re_λ)。
回顾大钝头体层流粗糙壁研究，Dirling等^[2]提出粗糙壁与光滑壁面的热流之比仅取决于Re_∞^1/2·k/δ^*，即表面粗糙度与边界层位移厚度之比和来流雷诺数的乘积，该雷诺数由钝头体的头部曲率半径R定义，Re_∞=ρ_∞V_∞R/μ_∞。同理，对于波状粗糙模型的流动问题，不考虑x影响的条件下，与其说影响壁面热流和阻力的参数是相对粗糙度h，不如说是粗糙度与边界层厚度之比。对于微粗糙板的高速边界层厚度并无一般解析解，但是不妨采用可压缩平板边界层理论来估算，由于本文中气体比热比、Pr数及壁温比均固定不变，则波状板上一个周期波长λ的壁面上，特征边界层厚度δ_λ应满足^[17]：$\delta _ { \lambda } / \lambda \sim F ( M a _ { \infty } ) / \sqrt { R e _ { \lambda } }$，F(Ma_∞)是关于Ma_∞的函数。进一步，在本文考虑的马赫数4~12的高超声速流动范围内，对于靠近板前缘的有限个波形内，黏性干扰效应较显著，平板的黏性干扰理论^[17]已经给出：无论是压强系数C_p还是摩阻系数C_f均由关联参数V=$\frac { M a _ { \infty } } { \sqrt { R e _ { x } } } \sqrt { C }$控制，其中$C = \frac{{{\rho _{\rm{w}}}{\mu _{\rm{w}}}}}{{{\rho _{\rm{e}}}{\mu _{\rm{e}}}}}$，由此推测，高速气流经过微粗糙壁面时，边界层厚度具有类似的形式：$\delta _ { \lambda } / \lambda \propto M a _ { \infty } / \sqrt { R e _ { \lambda } }$。因此，估算波状板上的粗糙度和波形特征边界层厚度的比值为

$\frac{{2\varepsilon }}{{{\delta _\lambda }}} = \frac{{2\varepsilon /\lambda }}{{{\delta _\lambda }/\lambda }} \propto \frac{{h\sqrt {R{e_\lambda }} }}{{M{a_\infty }}},$

定义等效粗糙度${\tilde h}$为

$\tilde h = h\sqrt {R{e_\lambda }} /M{a_\infty }.$

(9)

则波状粗糙板的流动和传热特性应取决于参数${\tilde h}$的控制，即平均热流-阻力系数比可写为

${C_{\rm{Q}}}/{C_{\rm{D}}} = f\left( {\tilde h} \right).$

为验证此结论，选择计算Ma_∞为4~14、Re_λ为5×10³~5×10⁴的算例，波状板长均为10λ。当粗糙度很小时，C_Q/C_D几乎不随x变化，对应取波状板中间部位的第5、6正弦波壁面对应C_Q/C_D的平均值为研究对象。图 8给出9组波状板算例中，小粗糙度下的C_Q/C_D与${\tilde h}$的对应关系，结果显示不同算例近似落在同一条曲线上，这说明波状粗糙壁上平均热流-阻力系数比C_Q/C_D仅取决于一个组合参数，即等效粗糙度${\tilde h}$。采用3次多项式拟合，可近似得到

$\begin{array}{*{20}{c}}{{C_{\rm{Q}}}/{C_{\rm{D}}} = f\left( {\tilde h} \right) = 0.20{{\tilde h}^3} - 0.36{{\tilde h}^2} - }\\{0.0018\tilde h + 0.62.}\end{array}$

(10)

Fig. 8

	Download: JPG larger image
图 8 不同Ma∞和Reλ下波状板附着流动CQ/CD随${\tilde h}$的变化规律 Fig. 8 The law of the CQ/CD change with ${\tilde h}$ at different Ma∞ and Reλ values along the wavy plate in attached flow 图 8 不同Ma∞和Reλ下波状板附着流动CQ/CD随${\tilde h}$的变化规律 Fig. 8 The law of the CQ/CD change with ${\tilde h}$ at different Ma∞ and Reλ values along the wavy plate in attached flow -->

式(10)由数值计算结果拟合得出，并无可对比的实验结果。在此仅从数值验证角度估算其误差范围：由第2节，壁面热流系数和摩阻系数最大偏差约为5%，而压强系数的最大偏差为10%，因此，可大致估计C_Q/C_D最大偏差约为15%~20%，可以供工程预研或工程估算时参考。
3.2.2 存在分离的波状板流动当流动在壁面附近出现分离时，壁面平均热流-阻力系数比呈现新的规律，C_Q/C_D等随x逐渐增大，而近似与h无关。例如图 5中，对固定的x₀，h>0.04时几条曲线几乎重合，取它们的C_Q/C_D平均值为研究对象，可以给出C_Q/C_D相对于x的变化曲线。对于不同的来流参数，该曲线显然是不同的，即C_Q/C_D=f(x, Ma_∞, Re_λ)。考虑到本文计算的是半无穷长的波状板模型流动，即便是对于分离流动情况，同样具有一定的自相似性特征。借鉴经典的平板边界层Blasius解^[18]结论——自相似变换中尺度的缩放因子是Re_x^1/2，引入波状板尺度的相似变换

$\tilde x = \bar x/\sqrt {R{e_\lambda }} .$

(11)

图 9给出波状板流动出现分离后，不同Ma_∞和Re_λ下，8组算例中C_Q/C_D与等效位置${\tilde x}$的对应关系，同样可以看到，不同算例点近似组成一条曲线。因此，存在分离的波状板流动中，壁面平均热流-阻力系数比可看作仅受${\tilde x}$的影响，同样，通过多项式拟合得到

${C_{\rm{Q}}}/{C_{\rm{D}}} = g\left( {\tilde x} \right) = - 31.7{{\tilde x}^2} + 4.34\tilde x + 0.30.$

(12)

Fig. 9

Download: JPG larger image

图 9 不同Ma∞和Reλ下波状板分离流动的CQ/CD随等效位置${\tilde x}$的变化规律 Fig. 9 The law of the CQ/CD change with the equivalent position ${\tilde x}$ at different Ma∞ and Reλ values along the wavy plate in separation flow 图 9 不同Ma∞和Reλ下波状板分离流动的CQ/CD随等效位置${\tilde x}$的变化规律 Fig. 9 The law of the CQ/CD change with the equivalent position ${\tilde x}$ at different Ma∞ and Reλ values along the wavy plate in separation flow -->

需要说明，虽然流动出现局部分离涡后，波状板上C_Q/C_D具有相似性规律，但流动在何处产生分离仍是一个难题。定性上讲，波状板来流Ma_∞越低，Re_λ越高，流动越容易分离，并且波状粗糙度越大，位置离前缘越远，流动越趋向于分离。然而暂时无法定量地给出流动开始分离的位置，这仍有待深入研究。
4 结论本文选取正弦波状粗糙元作为典型粗糙壁面外形，通过数值模拟分析在高速来流中波状板壁面热流与阻力系数的关联关系。结果表明：层流流动条件下，波状粗糙壁流动中壁面热流与摩阻没有明显的相关性，而与总阻力有直接关联，平均热流-阻力系数比C_Q/C_D是衡量流动与传热特性的关键参数。对于粗糙度非常小的波状板，流动具有附着型边界层流动特征，C_Q/C_D近似与x无关，仅取决于粗糙度与边界层厚度之比，即等效粗糙度${\tilde h}$；而波状板流动产生分离后，壁面附近流动存在具有相似的分离涡结构，对应的C_Q/C_D仅取决于相似变换尺度，即等效位置${\tilde x}$。综上：

$\begin{array}{l}{C_{\rm{Q}}}/{C_{\rm{D}}} = \\\left\{ {\begin{array}{*{20}{l}}{0.20{{\tilde h}^3} - 0.36{{\tilde h}^2} - 0.0018\tilde h + 0.62,未分离,}\\{ - 31.7{{\tilde x}^2} + 4.34\tilde x + 0.30,已分离.}\end{array}} \right.\end{array}$

无论从前人历史研究的结论看，还是从波状粗糙板流动本身的分析角度看，定量确定流动开始分离的位置仍然是一个非常困难的问题，本文仅讨论壁面热流-阻力的关系，分离条件仍有待下一步深入研究。
本文得到王智慧副教授的诸多指导和启发，在具体工作中也得到博士生李贤冬、罗健的热心帮助，在此谨表感谢。
参考文献

[1]	Abbett M J, Anderson A D, Cooper L, et al. Passive nosetip technology (PANT) program. Volume 20. Investigation of flow phenomena over reentry vehicle nosetips[R]. Acurex Corp/Aerotherm Mountain View CA, 1975.
[2]	Dirling J, Swain C, Stokes T. The effect of transition and boundary layer development on hypersonic reentry shape change[C]//10th Thermophysics Conference. 1975: 673.
[3]	Phinney R E. Mechanism for heat-transfer to a rough blunt body[J]. Letters in Heat and Mass Transfer, 1974, 1(2): 181-186. Doi:10.1016/0094-4548(74)90156-8
[4]	Grabow R M, White C. Surface roughness effects nosetip ablation characteristics[J]. AIAA Journal, 2014, 13(5): 605-609.
[5]	Yang S H, Lin G P, Xin S. Comparative study on the numerical computation of convective heat transfer over rough surfaces[J]. Journal of Aerospace Power, 2011, 26(3): 570-575.
[6]	李俊红, 张亮, 俞继军, 等. 高超声速可压缩流中粗糙壁热流研究[J]. 计算物理, 2017, 34(2): 165-174. Doi:10.3969/j.issn.1001-246X.2017.02.006
[7]	Gnoffo P A. Planetary-entry gas dynamics[J]. Annual Review of Fluid Mechanics, 1999, 31(1): 459-494. Doi:10.1146/annurev.fluid.31.1.459
[8]	Anderson J D. Fundamentals of aerodynamics[M]. 6th ed. NewYork: McGraw-Hill Education, 2017.
[9]	Nikuradse J. Laws of flow in rough pipes[M]. Washington: National Advisory Committee for Aeronautics, 1950.
[10]	Sigal A, Danberg J E. New correlation of roughness density effect on the turbulent boundary layer[J]. AIAA Journal, 1990, 28(3): 554-556. Doi:10.2514/3.10427
[11]	Kurganov A, Tadmor E. New high-resolution central schemes for nonlinear conservation laws and convection-diffusion equations[J]. Journal of Computational Physics, 2000, 160(1): 241-282. Doi:10.1006/jcph.2000.6459
[12]	Hayes W, Probstein R F. Hypersonic flow theory[M]. New York: Academic Press, 1966.
[13]	童秉纲, 孔祥言, 邓国华. 气体动力学[M]. 2版. 北京: 高等教育出版社, 2012.
[14]	Van Driest E R. Investigation of laminar boundary layer in compressible fluids using the Crocco method[J]. Technical Report Archive and Image Library, 1952, 10(1): 15-31.
[15]	姜贵庆. 高速气流传热与烧蚀热防护[M]. 北京: 国防工业出版社, 2003.
[16]	赵学端, 廖其奠. 黏性流体力学[M]. 北京: 机械工业出版社, 1983.
[17]	Anderson J D. Hypersonic and high-temperature gas dynamics[M]. 2nd ed. Reston: AIAA Publications, 2006.
[18]	庄礼贤, 尹协远, 马晖扬. 流体力学[M]. 2版. 合肥: 中国科学技术大学出版社, 2009.