摘要
激波/边界层干扰会导致飞行器表面出现显著的气动加热，而壁温升高后会反过来影响原有流场结构。为研究共轭传热情况下激波/湍流边界层干扰流场特性，采用数值模拟方法分析了不同壁面材料、气动加热时间和来流总温下的流场结构、壁面参数和固体域温度分布情况。结果表明，存在共轭传热时的分离区尺寸介于等温壁和绝热壁结果之间；随着壁面材料热扩散系数的减小和气动加热时间的增加，流场分离区尺寸增加，壁面温度峰值升高，热流峰值减小；相同材料和时间下，增大总温使分离区减小，说明改变来流总温引起的雷诺数变化对分离特性占主导地位。激波/边界层干扰使再附区下方的固体域温度明显上升，电木壁面内部的高温区扩散很慢，呈片状分布。
Abstract
Shock wave/boundary layer interaction results in significant aerodynamic heating on the surface of the aircraft. The original flow field structure is also affected by the increase in wall temperature. To investigate the flow field characteristics of shock wave/boundary layer interaction under CHT (conjugate heat transfer) conditions, numerical simulations were performed. The flow field structure, wall parameters and solid domain temperature distribution were analyzed under different wall materials, aerodynamic heating time and total temperature. The numerical results indicate that: the length of the separation bubble obtained using the CHT method is between the results for isothermal and adiabatic walls. As the thermal diffusivity of the wall material decreases and the aerodynamic heating time increases, the length of the separation bubble increases. At the same time, the wall temperature peak increases, and the heat flux peak decreases. For the same material and heating time, an increase in total temperature reduces the length of the separation bubble. This suggests that the change in Reynolds number caused by the total temperature dominates the separation characteristics. Additionally, for the solid domain at the corresponding position of the reattachment zone downstream of the interaction, a significant temperature rise is observed. The high-temperature zone inside the bakelite wall diffuses very slowly and exhibits a flaky distribution.
激波/边界层干扰（shock wave/ boundary layer interaction，SWBLI）广泛存在于高超声速飞行器内外流场中[1-3]。激波强度较大时会导致边界层发生分离，形成大尺度分离区，严重降低进气道的流量和抗反压能力。此外，SWBLI也会显著增加壁面热流，从而严重威胁飞行器结构可靠性。SWBLI导致的气动热问题被认为是实现长航时高速飞行所面临的最严峻挑战之一[4]。
有关的SWBLI的数值模拟研究大多采用理想壁面条件[5]，如绝热壁和等温壁。Yang等[6]基于绝热壁假设，通过求解雷诺平均纳维-斯托克斯（Reynolds-averaged Navier-Stokes， RANS）方程研究了SWBLI分离区特征尺度，提出了一种广义标度相关性来表征干扰强度对分离区流向特征长度的影响。Lopez和Harris[7]针对冷却壁条件下的SWBLI问题开展了大涡模拟（large eddy simulation，LES）研究，发现强冷却壁了分离区长度的显著缩短和低频压力波动的增加。同样基于等温壁假设，Zuo等[8]和Pasha等[9]的RANS研究表明干扰区尺度随壁面冷却和雷诺数增加而减小，马赫数和壁面温度的增加会降低摩阻和斯坦顿数。朱轲等[10]通过直接数值模拟（direct numerical simulation，DNS）发现随着壁温升高，摩阻因壁面附近速度梯度降低而减小，分离区尺寸增大。Bernardini等[11]和Volpiani等[12-13]分别针对不同壁温条件下的SWBLI进行了较为系统的DNS研究。
然而，SWBLI实际上存在对壁面的气动加热效应，并不满足等温和绝热假设，干扰特性会与气动加热时间和壁面材料等因素相关。由于固体域的传热时间尺度显著大于流体域，如果考虑气流与壁面之间的传热会极大地增加实验和高精度数值模拟的成本，因此有关气动加热问题的研究通常采用基于共轭传热（conjugate heat transfer，CHT）的RANS方法[14]。
CHT方法应用于处理同时包含流动换热和固体热传导的问题[15]，考虑了流场结构和物理量分布对固体域温度分布的影响，以及固体域温度改变后对流场的影响。通过与实验数据对比，Dechaumphai等[16]和Zhao等[17]验证了CHT方法的准确性和可靠性。John等[18]使用CHT方法对放置于高超声速气流中的圆柱形前缘进行了研究，对比了不同耦合方式下的计算效率。黄杰[19]采用CHT方法对高速飞行器头部扰流问题进行了流热固耦合计算，发现结构壁面温升对流场热流有显著影响。尹亮等[20]对头锥逆向喷流问题进行了CHT研究，发现逆向喷流结构显著降低了固体结构温度。周印佳等[21]利用CHT方法研究了超高温陶瓷材料的热物性非线性对高超声速流场传热过程的影响，发现材料热导率对气动热环境和表面温度分布影响显著。
目前，考虑共轭传热的SWBLI研究仍然较少。Muller等[22]将CHT与LES结合，研究了SWBLI的流场特性，结果表明，气流与壁面之间的传热会显著影响热流、温度、摩阻分布和湍流特性，然而由于计算量限制，计算总时长仅有约0.1 s。Peetala等[23]采用CHT方法对双楔压缩拐角问题进行了数值模拟研究，发现基于CHT的计算结果要比等温和绝热壁结果更为准确。
本文采用CHT方法对SWBLI问题进行了RANS方法研究，分析了不同壁面材料、气动加热时间和来流总温下的SWBLI流场结构以及壁面压力、温度和热流分布规律，最后研究了固体域温度分布情况。
1 数值方法及验证
1.1 数值方法
由于采用共轭传热方法，需要同时对流体域和固体域进行求解。针对流体域，采用二维可压缩Navier-Stokes方程作为控制方程：
(1)
其中，U为守恒型状态变量，F、G、Fv、Gv分别为流向、法向的通量及相应黏性项，t为时间变量，x、y分别为流向和法向的空间坐标。各项具体表达式如下：
(2)
(3)
(4)
(5)
(6)
其中，u、v分别为流向、法向速度分量，p、ρ、T分别为压力、密度和温度，κ为气体热传导系数，μ为动力黏性系数，采用Sutherland公式计算：
(7)
式中，μ0=1.716×10-5 N·s/m2，T0=273.11 K，Ts=110.56 K。
对流项采用Roe格式[24]用于捕捉激波，空间离散采用二阶迎风格式，黏性项使用二阶中心差分格式，通过引入理想气体模型封闭方程组。湍流模型采用Menter[25]发展的剪切应力传输（shear stress transfer，SST）k-ω模型，湍动能K、耗散率ω的输运方程为：
(8)
(9)
其中，GK、Gω分别为湍动能和耗散率的生成项，ΓK、Γω分别为湍动能和耗散率的有效扩散项，Dω为交叉扩散项，式（3）和式（4）中的最后一项表示湍流引起的湍动能和耗散率的耗散。
对于固体域，采用二维非定常热传导方程进行计算：
(10)
式中，ρs、cs、λs分别为固体材料的密度、比热、热导率，选择铜、铝、不锈钢、电木进行计算，不同材料的热物性参数如表1所示。其中，热扩散系数为热导率与密度和比热的乘积之比，用于表征物体温度变化的快慢。热传导方程采用中心差分离散，采用经典一阶Runge-Kutta方法对时间进行离散。
表1壁面材料参数
Tab.1Material parameters of the flat
共轭传热问题涉及固体外部对流换热和固体内部热传导的耦合。由于固体的热容远大于流体，近似认为当流场达到稳态时，固体导热尚未开始，因此首先按环境温度设置初始温度，计算稳态流场，再以此为初始条件进行非定常耦合传热计算。由于流场发展相比固体传热过程的特征时间非常小，可以近似认为在固体域传热进行的一个时间步长内，流场的特性没有发生太大改变，流场的计算是瞬时稳定的，即松耦合方法。
具体求解思路如图1所示：
1）初始化流体域和固体域，其中流体域以初始壁面温度Tw作为Dirichlet边界条件，固体域按照初始温度进行初始化，通过稳态计算获得初始流场和耦合界面的热流值qw，结果作为瞬态耦合计算的初始条件。
2）在共轭传热计算起始点，进行瞬态求解。将耦合界面的热流值qw作为Neumann边界条件，在固体域求解传热方程，经过固体域一个时间步长Δt后，获得重新分布的固体域温度场，得到此时的壁面温度Tw。
3）将Tw作为Dirichlet边界条件，在流体域进行瞬态计算，经过流体域一个时间步长Δt后，得到此时的壁面热流值qw。
4）循环此过程，直到完成给定的计算时间。
图1共轭传热方法示意图
Fig.1Schematic of conjugate heat transfer method
1.2 仿真模型及工况
设置的计算域如图2所示，气流方向由左至右，来流马赫数为6.0，总压1.0 MPa。计算域上边界设置斜劈用于产生斜激波，气流转折角为θ=10°。为了避免斜劈表面边界层引起的激波强度变化，将斜劈所在上边界设置为滑移壁面。入口边界采用压力远场边界条件，出口边界设置为超声速压力出口。在此基础上，本文研究了不同壁面材料、时间和来流总温下的SWBLI流场特性，算例工况的具体设置如表2所示，其中等温壁温度设置为300 K。
图2SWBLI仿真模型示意图
Fig.2Geometry and boundary conditions of the SWBLI model
1.3 数值方法验证
通过NASA Langley圆柱壳前缘高超声速流动实验[26]验证CHT方法的准确性。实验中，来流马赫数为6.47，来流静温为241.5 K，静压为648.1 Pa，圆管材料为321不锈钢，内径为25.4 mm，外径为38.1 mm。图3（a）和图3（b）分别展示了仿真结果的流体域马赫数云图和固体域温度云图，计算得到的流场形态和激波位置符合实验结果，圆管最前端温度最高，这是流场驻点气动加热的结果。图3（c）为5 s时的圆管表面温度分布，实验数据为约5 s时的数据[20]。可以发现，仿真结果与实验数据较为一致，证明本文采用的共轭传热方法能够较好地应用于高超声速流动和传热问题的计算分析。
表2算例工况设置
Tab.2Cases condition setup
为了验证计算SWBLI问题的准确性，对Kussoy和Horstman的实验[27]进行了数值模拟。实验中，来流马赫数为8.18，静温为81 K，静压为430 Pa，实验模型如图4（a）所示，斜劈前缘与平板前缘的距离x0设置在120~140 cm的范围内，距平板高度为10.16 cm，斜劈上表面与水平面之间的角度θc和斜劈角度θm均设置为10°。图4（b）和图4（c）分别展示了实验和仿真中沿流向的壁面压力和热流分布，可以看出数值模拟结果与实验数据吻合较好，证明本文采用的数值方法能够获得具有较高可信度的SWBLI仿真结果。
图3CHT验证算例计算结果
Fig.3Validation of CHT method
为了检验网格无关性，选取三套不同疏密程度的网格进行先期计算，网格数量如表3所示，验证情况如图5所示。所有计算网格下壁面第一层网格的无量纲壁面高度y+值均小于1。三套网格计算所获得的壁面压力和热流结果如图5（a）、图5（b）所示，对比发现，所选三套网格均能较准确地预测壁面的压力分布，当网格量在中等及以上时，计算结果几乎不受网格尺寸影响，综合考虑计算精度及效率，采用中等网格进行计算。
图4SWBLI验证算例计算结果
Fig.4Calculation results of the SWBLI verification example
表3网格数量
Tab.3Flow field grids
图5网格与迭代次数无关性验证
Fig.5Comparison of wall parameter predicted on different grid levels and iterations
选择上述中等网格进行共轭传热模拟，时间步长为0.005 s，分别选取迭代次数10、30、50、70进行迭代次数验证。壁面温度结果如图5（c）所示，在迭代50次以上时，计算结果几乎不变，为了保证精度和节约计算资源，后续的计算采用迭代次数50。
2 结果与讨论
2.1 SWBLI流场结构分析
图6显示了不同壁面边界条件（等温壁、绝热壁和耦合壁）下SWBLI流场的压力云图。其中，耦合壁指考虑共轭传热的算例，流场结果取来流总温为1 000 K、t=10 s时的仿真数据。为方便分析，将流场中的声速线标记在云图中，使用白色虚线表示。如图6所示，所有工况下的流场结构基本一致，气流经过入射激波后压力骤升，随后在膨胀波作用下压力下降，膨胀波影响SWBLI流场的区域位于干扰区下游。入射激波导致的压升通过壁面附近的亚声速区向上游传播，形成分离区并抬升上游的边界层，形成一道分离激波，引起下游压力的增加，但从图中可知，压力升高并不明显。随后，入射激波与分离激波相交产生明显的高压区，入射激波在透射后继续向壁面方向传播。然而由于当地马赫数逐渐减小，透射激波逐渐弯曲，并最终在声速线位置消失，过程中不断向上方反射出膨胀波，削弱了激波相交诱导的高压区。抬升后的边界层在绕过分离区后撞击壁面，气流的突然偏折形成一道再附激波，引起下游压力大范围的显著增加。对比不同壁面边界条件下的流场可以发现，绝热壁情况下流场分离程度最大，激波相交产生的高压区范围远大于其他工况，而铜、铝、钢和电木的压力分布情况与等温壁下的结果近似一致，流场分离程度均较小。
图6不同壁面边界条件下SWBLI流场的压力云图
Fig.6Pressure contours of SWBLI with different wall boundary conditions
图7给出了不同壁面边界条件下SWBLI流场的马赫数云图，其中耦合壁算例同样取总温为1 000 K、t=10 s时的流场数据。图7中标记的黑色虚线为沿流向速度分量为0的等值线，用以描述分离区特征。如上所述，绝热壁条件下分离区尺度最大，而其他工况下尺度相当，其中壁面材料为电木时的分离区要稍大于其他材料。分离区形态可以近似看作一个扁平的三角形，从分离区顶点到分离起始点的距离要显著大于其到再附点的距离。为定量分析激波诱导的分离区特性，图8（a）给出了不同壁面边界条件下沿流向的壁面摩阻曲线，用三种不同符号表示分离区的三个特征位置。实心圆表示干扰作用的起始点，实心正方形表示分离点，实心三角形表示再附点，分离点和再附点分别为壁面摩阻曲线的两个零点。可以发现，不同壁面材料下干扰起始点与分离点位置略有不同，而再附点位置几乎相同，且位于无粘激波入射位置附近。将图8（a）中的分离点和再附点位置提取出来，统计得到不同壁面边界条件下流场分离、再附位置以及分离区长度，结果如图8（b）所示。结合表1可知，分离区长度与壁面材料的热扩散系数呈负相关，热扩散系数越大分离程度越接近等温壁结果，反过来越接近绝热壁结果。即使对于电木这种热扩散系数很小的材料，在气动加热10 s的情况下流场特性仍然更趋近于等温壁结果。
图7不同壁面边界条件下SWBLIs流场的马赫数云图
Fig.7Mach contours of SWBLIs with different wall boundary conditions
图8不同壁面边界条件下分离区变化情况
Fig.8Variations in the separation zone under different wall boundary conditions
对于电木材料，将来流总温为1 000 K时不同时间下的分离点和再附点位置提取出来，得到图9（a）。对比发现，随着时间的增加，分离区长度逐渐增大，这是由于气动加热效应使壁面温度升高，从而加重了流动分离[8]。值得注意的是，分离点随时间增加近似线性地向上游移动，而再附点位置几乎不变，因此分离区长度的增加量主要取决于分离点的前移程度。
图9（b）展示了电木壁面在不同总温条件下t=10 s时的分离特征。可以发现，随着来流总温的增加，分离区长度并没有增大，而是减小，主要体现在分离点向下游移动，而再附点近似不变。这是因为总温增加后虽然促进了气动加热，但也改变了来流雷诺数。将图9（b）中四种不同总温工况的雷诺数ReL和边界层厚度δ0与其对应分离区长度Lsep表示在图10中，发现Lsep与δ0·（ReL）0.5近似成正比。
图9不同时刻和总温下分离区变化情况
Fig.9Variations in the separation zone of different time and total temperature
2.2 SWBLI壁面参数分析
2.2.1 压力分布
图11（a）显示了不同壁面材料在总温为1 000 K、t=10 s时的壁面压力曲线。其中，干扰起始点、分离点和再附点分别在图中用对应曲线颜色的实心圆、正方形和三角形进行标记。可以观察到，压力存在两次明显的上升过程，自干扰起始点，压力第一次上升，上升速率逐渐变缓，所有工况会在x=340 mm附近达到近似同一个平台值。根据自由干扰理论，从干扰起始点到分离点为自由干扰区，该区域的压力分布仅与来流条件有关。在下游，伴随着流动再附，压力经历第二次增加，并在下游x=390 mm附近达到峰值，随后开始在上游膨胀波的作用下逐渐下降（参考图6）。对比不同材料下的压力分布，发现电木的压力分布曲线与其他工况存在一定差异，这是由于电木热扩散系数远小于其他材料，在相同时间下壁温更高，从而影响流场，使干扰起始点和分离点向上游移动，分离区增大。
图10分离区长度与δ0·（ReL）0.5的关系
Fig.10Influence of the δ0· (ReL) 0.5 on the separation length
图11（b）和图11（c）分别给出了壁面材料为电木时，同一总温（T0=1 000 K）不同时间下的压力分布曲线以及同一时间（t=10 s）不同总温下的压力分布曲线。壁面压力随时间的变化主要体现在压力第一次升高的起始位置，如上所述，随着时间增加，干扰起始点前移。在压力二次上升后，不同时间下的结果基本一致。对于不同来流总温下的压力分布，干扰起始点和下游压力峰值均有一定的差异，随着总温升高，干扰起始点会向下游移动，后移程度与总温变化呈现出较为明显的非线性特征，壁面压力达到峰值的流向位置基本不变，但数值略有下降。
图11壁面压力分布
Fig.11Distributions of wall pressure
2.2.2 温度和热流分布
图12（a）展示了不同材料在总温为1 000 K、t=10 s时的壁面温度分布曲线，同样标记了干扰起始点、分离点和再附点的位置。如图12（a）所示，在干扰开始后，壁面温度略有下降，随后在分离点下游一定位置处开始升高。壁面温度在达到干扰区温度后，会出现一个快速上升的过程。对比图11，发现这一转折点对应的流向位置刚好在压力平台位置附近。此外，与压力分布曲线的规律相似，壁面温度在再附点下游达到温度峰值，但是峰值对应的流向位置更靠近上游一些。整体来看，在同一时间，壁面材料热扩散系数越小，同一流向位置处壁温越高，沿流向方向的壁面温差也越大。图12（b）给出了不同材料在总温为1 000 K、t=10 s 时的壁面热流分布曲线，其趋势与壁温的变化趋势相同。耦合壁条件下的热流均低于等温壁的结果，且壁面材料热扩散系数越大，同一时间的壁面热流越接近于等温壁情况。考虑到不同壁面材料下干扰起始点壁面温度和热流不同，对壁面温度和热流使用干扰起始点的对应参数进行无量纲处理，得到分离点和再附点的无量纲温度和热流，结果如表4所示。电木壁面分离点处无量纲温度最低，但再附点处无量纲温度最高，而无量纲热流的规律相反。同一时间，随壁面材料热扩散系数减小，再附点处无量纲温度增加，无量纲热流减小。
图12不同材料壁面温度和热流分布
Fig.12Distributions of wall temperature and heat flux for different wall materials
表4无量纲温度和热流
Tab.4Nondimensional temperature and heat flux
图13（a）和图13（b）分别给出了电木壁面条件下总温为1 000 K时壁面温度和热流分布随时间的变化曲线。可以发现，壁面温度随时间不断升高，但升高的趋势放缓，壁面热流则随时间下降，同样下降速率会变小。在再附点下游，壁面温度和热流的变化更为明显，这是由于再附区上方的气流经过分离激波和再附激波压缩后温度急剧升高，使得这部分壁面被显著加热，而随着壁温升高，气流与壁面之间的温差降低，从而导致热流下降。
如上所述，温度和热流的变化速率会随时间的增加而减小。根据图13，壁面温度和热流曲线的峰值点位置均在x=390 mm附近，为此提取不同壁面材料下x=390 mm处的温度和热流作为对应的峰值，得到其随时间的变化情况，结果分别如图14（a）和图14（b）所示。可知，随着时间增加，金属材料（铜、铝、钢）的峰值温度和热流近似线性变化，且热扩散系数越大，拟合曲线斜率的值越小。这是由于较大热扩散系数的材料热量扩散得更快，导致壁面温度和热流上升速率较慢。而对于热扩散系数较小的非金属材料（电木），峰值温度和热流与时间近似呈指数关系，温度随时间的变化趋于稳定，热流随时间减小，其变化率衰减。
图15（a）和图15（b）分别给出了不同来流总温下t=10 s时电木壁面的温度和热流沿流向的变化曲线。可以发现，壁面温度和热流的变化规律一致，随着总温升高，壁面温度呈现整体性升高的趋势，干扰区上下游的温差和热流差也随总温升高而增大，温度和热流的峰值位置略向下游移动。壁面热流和温度变化规律的相似性表明，增大来流总温虽然会使壁温更快地升高，但壁面与外界气流之间的温差仍然会增大。
图13电木不同时间壁面温度和热流分布
Fig.13Distributions of wall temperature and heat flux for different time of bakelite
图14壁面温度和热流最大值分布
Fig.14Distribution of maximum wall temperature and heat flux
图15电木不同总温壁面温度和热流分布
Fig.15Distributions of wall temperature and heat flux for total temperature
2.3 SWBLI固体域导热特性
图16分别给出了不同来流总温、不同材料下固体域的温度云图。为便于分析流场结构与固体域温度分布的联系，在固体域温度上方给出了马赫数云图，其结果取自同一材料、总温为1 000 K时的流场数据。如图16所示，对于同一材料和相同时间，随着总温升高，再附区下方固体域的温度明显增大。对比不同材料下的结果，发现热扩散系数最大的铝壁面高温范围远大于其他工况，铜和铁壁面次之。在来流总温较小时，高温区近似以x=0.385 m为中心向固体域扩散，而随着总温增大，扩散中心逐渐向下游移动，这与图13中温度峰值位置的变化规律一致。对于热扩散系数最小的电木壁面，热量更多地集中在壁面附近，高温区呈现片状，局部温度上升更为明显。
图16不同材料在不同总温下固体域温度云图及对应流体域马赫数云图
Fig.16Temperature contours of solid filed for different materials under different total temperature with corresponding Mach contours
图17给出了电木壁面在来流总温为1 000 K时，不同时间（t=2 s、4 s、6 s、8 s和10 s）下固体域温度分布云图，以及在t=0时流场的压力云图。随着时间推移，固体域内温度升高的区域同样位于壁面附近，高温区向下方扩展的速度很慢，温度分布变化情况与改变总温的规律相似。因此，对于不同总温和气动加热时间的干扰流场，固体域温度分布可能是相似的。而对比不同壁面材料的结果（见图16），温度的空间分布存在明显差异。因此可以推断，如果两种壁面的热扩散系数不同，则难以找出一个来流总温和时间使它们的固体域温度分布具有相似性。
图17不同时刻电木固体域温度云图及对应流体域压力云图
Fig.17Temperature contours of bakelite solid filed for different time with corresponding pressure contours
3 结论
现有的SWBLI数值模拟研究大多采用等温或绝热壁面条件，并没有考虑结构传热。本文采用共轭传热方法对不同壁面材料下高超声速SWTBLI问题进行数值模拟研究，得到的结论如下：
1）壁面材料为热扩散系数较大的铜、铝和钢时，气动加热时间t=10 s下的分离区长度仍与等温壁结果相当，而对于热扩散系数较小的电木材料，分离区长度明显增大，但小于绝热壁结果。
2）对于相同来流总温和气动加热时间，随着壁面材料热扩散系数的减小，分离点位置明显前移，再附点位置稍向下游移动，分离区尺度增大，壁面压力峰值不变，温度峰值升高，热流峰值相应减小。
3）对于相同来流总温和壁面材料，随着气动加热时间增加，分离区特性及壁面参数分布的变化情况与减小壁面热扩散系数的规律一致。此外，电木壁面的温度和热流峰值随时间的变化率逐渐减小，而其他壁面材料下温度和热流峰值近似线性变化。
4）对于相同壁面材料和气动加热时间，随着总温升高，分离点和再附点分别向下游和上游移动，分离区减小，表明改变来流总温引起的雷诺数变化对分离特性占主导地位。增加总温导致壁面压力峰值略有下降，温度和热流峰值均上升。
5）SWTBLI使再附区下方的固体域温度显著升高，高温区近似以温度峰值对应的位置为中心向固体域内部扩散。对于热扩散系数较小的材料，高温区紧贴壁面，呈现片状。不同总温和气动加热时间下固体域的温度分布可能具有相似性。
