高超声速飞行器以7km/s以上的速度再入大气层时,遭受强烈的气动加热,飞行器物面辐射平衡温度高达2 000K以上,无法使内部装置正常工作。为此再入飞行器表面多采用烧蚀材料并通过高温气体与物面材料发生的气固热化学反应和质量交换,吸收大部分气动加热能量,使飞行器内壁保持允许的温度[1]。烧蚀产物引射到气体边界层中,影响到流场的化学反应和热力学特性,并进一步影响飞行器的光电特性。关于战略弹头等钝头体烧蚀流场特性,国内外已开展了大量深入的工作[2-8],取得了富有成效的结果。本文在前人研究基础上,将烧蚀流场的研究扩充到其他高超声速飞行器外形或部件模型。
对再入飞行器和吸气式高超飞行器而言,控制面和发动机进气道区域的流场是需要重点研究的,因为在这些区域很容易发生激波/边界层干扰,影响控制面的效率和发动机的性能[9]。这些区域的流动可以简化为平板后接压缩面的压缩拐角流动。平板边界层从平板区逐渐发展,在拐角前会受到压缩面产生的斜激波的干扰。如果该斜激波足够强,边界层会发生分离,并在拐角区域形成回流区。分离使边界层变厚且压缩来流,产生分离激波。在再附点附近,气流受到压缩面压缩产生一系列压缩波,这些压缩波汇聚成一道激波,和分离激波相互作用,影响到压缩面壁面特性。在高焓来流条件下,高温气体效应将会对流场中的波系和分离区产生影响[9]。本文重点研究存在烧蚀的压缩拐角流动特性。在自由流Ma=10~30,总焓值6~55MJ/kg范围,对15°、18°、24°压缩拐角模型,分别进行无烧蚀的壁面催化与非催化条件和石墨烧蚀条件下的流场计算,分析各类条件下压缩拐角的流场特性和壁面条件的影响。
1 热化学模型石墨材料表面上发生的反应主要有碳的氧化反应,碳催化的O原子复合反应和碳的升华反应(生成C,C2,C3),具体反应式和速率常数确定详见文献[3]。据此可以得到表面反应生成的各组元净质量流率。
对存在石墨烧蚀的高温空气流场计算时,考虑16个组元,29个化学反应,详见表 1。反应式中,催化物M1、M2、M3、M4、M5、M6为各反应的催化物,就是流场中的某些组元。对反应1~反应21,速率系数采用文献[10]给出的方法,其余反应采用文献[11]的方法计算。
| 序号 | 反应方程式 |
| 1 | N2+M1⇄N+N+M1 |
| 2 | O2+M2⇄O+O+M2 |
| 3 | C2+M3⇄C+C+M3 |
| 4 | CN+M4⇄C+N+M4 |
| 5 | N2+e-⇄N+N+e- |
| 6 | O+e-⇄O++e-+e- |
| 7 | N+e-⇄N++e-+e- |
| 8 | N2+O⇄NO+N |
| 9 | NO+O⇄O2+N |
| 10 | CO+C⇄C2+O |
| 11 | CO+O⇄O2+C |
| 12 | CO+N⇄CN+O |
| 13 | N2+C⇄CN+N |
| 14 | CN+O⇄NO+C |
| 15 | CN+C⇄C2+N |
| 16 | CO+C2⇄C3+O |
| 17 | C3+N⇄CN+C2 |
| 18 | C3+C⇄C2+C2 |
| 19 | O+N⇄NO++e- |
| 20 | N+N⇄N2++e- |
| 21 | CO2+M6⇄CO+O+M6 |
| 22 | CO2+N⇄CN+O2 |
| 23 | CO+NO⇄CO2+N |
| 24 | CO2+O⇄CO+O2 |
| 25 | 2CO⇄CO2+C |
| 26 | N2+CO2⇄N+N+CO2 |
| 27 | O2+CO2⇄O+O+CO2 |
| 28 | C2+CO2⇄C+C+CO2 |
| 29 | CN+CO2⇄C+N+CO2 |
热力非平衡条件下,系统有两个温度,平动/转动温度T和振动温度Tv,振动松弛时间采用Park修正的Millikan和White振动松弛模型计算。双温度近似下,不同的反应采用不同的控制温度。气体热力学特性和输运特性的计算方法详见文献[10]。
2 流动控制方程与壁面边界条件 2.1 控制方程与计算方法控制方程为时间相关的二维非平衡流N-S方程,求流场的定常解。无粘通量项采用AUSMPW+格式[12]离散,粘性项采用中心差分离散。流体力学方程和热化学动力学方程耦合求解时,采用全耦合隐式处理无粘通量项和源项,隐式处理方法为LU-SGS方法[13]。控制方程与求解的详细内容参见文献[14]。
2.2 烧蚀壁面边界条件与求解壁面处s组元的质量平衡关系为
|
(1) |
采用等离子体准中性假设后,可不含电子的质量平衡关系,式(1)共有15个方程。式中n为壁面单位法矢量,


壁面的能量平衡关系为
|
(2) |
其中, qw为垂直于壁面向上的热流, qcond, w为来自固体材料的垂直于壁面向上的热传导热流。这里考虑准定常烧蚀情况,将qcond, w设为零。式(2)中第三项为从壁面辐射的能量,ε是材料的辐射系数,σ是斯蒂芬-玻尔兹曼常数。式(2)中第4项为气固交界面处由于表面反应从固体进入交界面的能流

存在壁面质量引射时,壁面法向速度由总质量流率确定
|
(3) |
上述方程中出现的壁面密度根据混合气体的状态方程确定,而壁面压力根据法向压力梯度为零确定。
联立求解各组元的壁面质量平衡关系式和壁面能量平衡关系式共16个非线性方程,即可确定壁面处的各组元质量分数和壁面温度。为求解方便,也可先假定壁温,求解质量平衡关系得到质量分数后再通过能量平衡关系式求解温度,之后进行迭代。非线性方程组的求解可采用牛顿迭代法。
另外,为对比分析壁面条件对流场特性的影响,本文还进行了三种无烧蚀壁面条件下的流场计算,分别是:等壁温(300K)的全催化和非催化条件,辐射平衡壁温的非催化条件。全催化壁面处质量分数为当地温度压力下的平衡值,非催化壁面处组元质量分数根据梯度为零确定。辐射平衡壁温根据能量平衡关系-qw=εσTw4确定。无烧蚀时壁面处速度为零。
3 计算结果与分析 3.1 计算程序的验证作为对本文计算程序模拟存在激波/边界层干扰和流动分离的高超声速流场有效性的验证,作者曾对Holden[15]在Calspan 48-inch的激波风洞实验条件下的二维压缩拐角流场进行了计算。风洞自由流马赫数14.1、温度89K、密度5.27×10-4kg/m3,计算得到的压缩拐角壁面压强和热流分布与文献[16]给出的实验测量结果和计算结果一致,结果对比图详见文献[17]。作为对本文含烧蚀非平衡流场计算程序的验证,曾对再入条件(飞行速度10km/s、高度65km)下的球头石墨烧蚀流场进行计算,所得壁面热流等结果与文献[3]的计算结果一致性也较好,结果对比图详见文献[8]。
3.2 算例条件与计算网格取上述的Holden在激波风洞中开展实验的二维压缩拐角实验模型(模型总长0.6096m,其中前段平板长0.3048m;模型宽0.6096m)。自由流密度设为风洞条件(5.27×10-4kg/m3),但自由流温度设为300K,以体现高超声速流场的高温效应。在斜面倾角为15°、18°、24°三种情况下,分别计算了马赫数为10、14、20、30时的石墨烧蚀流场。作为对比,还进行了无烧蚀流场的计算,壁面分别设定为壁温300K条件下的全催化壁、非催化壁、非催化辐射平衡壁温。
采用131×81(流向×法向)的网格。流向方向上,在头部、拐角处、压力峰值处均进行了加密。法向方向在壁面处进行了加密,第一层网格高度为6.096×10-5m,法向外边界与平板前缘产生的激波大致平行。具体网格布置如图 1,图中网格进行了稀疏显示处理。
|
| 图 1 计算网格 Fig. 1 Computational mesh |
3.3 流场结构特点分析
高超声速压缩拐角流场中的波系结构如图 2[17]。该图基于马赫数7.7条件下24°压缩拐角流场的压力等值线绘制。该条件下流动发生了分离,流场中除前缘激波和主激波外,还有较强的分离激波和再附激波。对未发生流动分离的情况,在拐角处前方,气流经过前缘激波后参数也已不均匀,斜面压缩在此产生的激波不再是斜激波,而是有一定的弯曲,且强度弱于主激波。这一段弯曲激波与流动分离时的再附激波作用类似,本文也将其命名为“再附激波”,方便后面的统一分析。前缘激波、分离激波、再附激波都可能穿透主激波或和主激波相交、融合,造成相交融合点附近的压力温度升高并引起峰值热流的急剧升高。具体的波系结构除了与自由流条件和压缩面角度有关,与拐角处平板的边界层状态也密切相关。
本文诸自由流条件下,15°压缩拐角流动均未出现分离,18°压缩拐角流动仅Ma=10的烧蚀条件下有轻微分离,24°压缩拐角流动在Ma=30时无分离,在其它较小Ma均出现明显分离。这里以24°压缩拐角、Ma=14算例为代表,在图 3中给出了全催化壁(用Cat.表示,壁温300K)和烧蚀条件下(Abaltion)流场的压强等值线和分离区流线。图 4、图 5分别给出了壁面压强系数和热流分布,图 6给出了分离点附近的壁面摩阻分布。该算例流场化学反应微弱,等壁温条件下,壁面的催化特性对壁面热流和分离区范围都没有明显影响,烧蚀条件下的分离区范围比无烧蚀的辐射平衡壁温条件稍有增大。相对于低壁温条件,无烧蚀的辐射平衡壁温和壁面烧蚀条件下壁面热流有所降低。但流动的激波层增厚,分离区增大,这是高壁温导致边界层增厚从而影响激波形状及其与边界层相互作用特性的结果。分离特性的不同还将进一步影响到下游的流动特性,分离区增大后分离激波、再附激波与主激波相交于更远位置,斜面上压力、摩阻和热流峰值点均有所后移。
|
| 图 3 Ma=14时24°压缩拐角流场压强分布与回流区流线 Fig. 3 Pressure distribution and streamline in recirculation zone at Ma=14 |
|
| 图 4 Ma=14时24°压缩拐角壁面压强系数分布 Fig. 4 Surface pressure distribution at Ma=14 |
|
| 图 5 Ma=14时24°压缩拐角壁面热流分布 Fig. 5 Surface heat transfer distribution at Ma=14 |
|
| 图 6 Ma=14时24°压缩拐角分离点附近摩阻系数分布 Fig. 6 Friction distribution near the separation point at Ma=14 |
斜面倾角和壁面条件都相同的情况下,随着马赫数增大分离可能性减小。这是因为一方面马赫数增大后,压缩面形成的主激波激波角减小,激波更贴近物面,逆压梯度的范围减小;另一方面马赫数增大后边界层内的亚声速区范围也有所减小,逆压梯度向前传播的的范围也减小。马赫数越大,分离区越小,压缩面上压力峰值点越靠前。表 2给出了24°压缩拐角流动存在分离的各算例的分离起始位置x/L,其中L为平板长度,x为离平板前缘距离。
| x/L | ||||
| 壁面条件马赫数 | 催化壁壁温300K | 非催化壁温300K | 无烧蚀辐射平衡壁温 | 壁面石墨烧蚀 |
| 10 | 0.6588 | 0.6612 | 0.5136 | 0.5018 |
| 14 | 0.7383 | 0.7408 | 0.5852 | 0.5554 |
| 20 | 0.8712 | 0.8732 | 0.7731 | 0.7717 |
3.4 Ma=30时24°压缩拐角流场热化学参数分布
24°压缩拐角在Ma=30时的流场是激波最强、振动激发和化学反应也最强的情况。下面以24°压缩拐角在Ma=30、烧蚀壁面条件下流场热化学参数分布为例进行分析。图 7给出了流场的压力等值线,图 8分别给出了平动、振动温度等值线。由图 8(a)可见在斜面前端出现高平动温度区,这是经过前缘激波压缩升温后的气流和拐角区较厚的边界层内升温了的气流,又经过“再附激波”的进一步压缩和升温所致。而由图 7可见压力峰值区出现在主激波后,自由来流经过强主激波的压缩和升温,压力达到全场峰值,但由于自由来流温度较低(300K),温度仍低于“再附激波”后的气流。振动温度的峰值区在主激波后,这与振动激发随压力升高而加强有关。主激波较强的部分下游存在一条高温带,高温带中振动激发和化学反应明显要比周围流场剧烈,原子质量分数的高值区和振动温度的高值区范围接近。图 9、图 10分别给出了流场的O、NO+质量分数分布。
|
| 图 7 Ma=30时24°压缩拐角流场压强分布 Fig. 7 Pressure distribution at Ma=30 |
|
| 图 8 Ma=30时24°压缩拐角流场平动与振动温度分布 Fig. 8 Translational temperature and vibrational temperature distribution at Ma=30 |
|
| 图 9 Ma=30时24°压缩拐角流场O质量分数分布 Fig. 9 CO distribution at Ma=30 |
|
| 图 10 Ma=30时24°压缩拐角流场NO+质量分数分布 Fig. 10 CNO+ distribution at Ma=30 |
烧蚀产物基本上局限在边界层内,质量分数最高的烧蚀组元依次是CO、C3、C2、CN。图 11、图 12分别给出了CO、C3组元的质量分数分布。另外,C2质量分数峰值达到0.0574,CN值为0.0230,C峰值为0.0211,CO2峰值在10-3量级。烧蚀产物在接近拐角斜面尾部时达到峰值,图 13给出了x/L=1.9时的烧蚀组元质量分数沿斜面法向分布。
|
| 图 11 Ma=30时24°压缩拐角流场CO质量分数分布 Fig. 11 CCO distribution at Ma=30 |
|
| 图 12 Ma=30时24°压缩拐角流场C3质量分数分布 Fig. 12 CC3 distribution at Ma=30 |
|
| 图 13 Ma=30时24°压缩拐角流场烧蚀组元法向分布 Fig. 13 Ablation species normal distribution at Ma=30 |
图 14给出了四种壁面条件下的壁面热流分布。相对于低壁温(300K)的催化和非催化条件,无烧蚀辐射平衡壁温和烧蚀条件下壁温有较大升高,所以拐角斜面上的热流有明显降低。其中低壁温条件下催化壁热流高于非催化壁,是因为算例条件下主激波后有较强的非平衡化学反应,催化壁条件下壁面处存在扩散热流。烧蚀条件热流又高于无烧蚀的辐射平衡壁温条件,也是因为烧蚀壁面处的放热化学反应。
|
| 图 14 Ma=30时24°压缩拐角壁面热流分布 Fig. 14 Distribution of surface heat flux at Ma=30 |
3.5 壁面条件对流场特性影响的分析
在3.3节中已知壁面条件对流场总体结构的影响特点。相对于低壁温条件,无烧蚀的辐射平衡壁温和壁面烧蚀条件下斜面上的热流率降低,但分离可能性和分离区范围增大,并且“再附激波”和主激波的相交位置后移,导致壁面压力和热流峰值区后移。无分离时的影响特点类似,但程度较小。
另外,无烧蚀的辐射平衡壁温和烧蚀条件下的激波层厚度和激波强度明显大于低壁温条件,并且烧蚀条件还比辐射平衡壁温条件下略大。这是壁温升高或烧蚀使边界层增厚,因而前缘激波和分离激波增强,激波也被推离物面所致。这里以斜面上压力峰值点位置为例,在图 15、图 16中给出了24°压缩拐角在Ma=14(存在流动分离)和Ma=30(无分离)条件下温度沿壁面法向的分布。斜面上压力峰值点位于再附激波和主激波相交区域。壁面条件不影响压力峰值,但高壁温条件下该处温度峰值更高;并且壁面条件的影响在存在流动分离的Ma=14算例上表现比无分离的Ma=30算例更明显。
|
| 图 15 Ma=14时24°压缩拐角斜面上压力峰值点处温度沿法向分布 Fig. 15 Normal distribution of temperature at pmax at Ma=14 |
|
| 图 16 Ma=30时24°压缩拐角斜面上压力峰值点处温度沿法向分布 Fig. 16 Normal distribution of temperature at pmax at Ma=30 |
由于激波层厚度的差异,低壁温和高壁温条件(无烧蚀的辐射平衡壁温和壁面烧蚀条件)下的组元分布也表现出相应的明显差异。不过无烧蚀的辐射平衡壁温和烧蚀两种壁面条件下组元分布差异主要体现在边界层内。烧蚀产物扩散到边界层内后进一步影响边界层内的化学反应和扩散,使高温空气的中性分子组元N2、O2增加,原子组元N、O和主要的离子组元NO+、N2+减少。图 17、图 18给出了斜面上压力峰值点处的N、NO、O组元和NO+、N2+组元质量分数沿法向分布,可以看出烧蚀的上述影响。
|
| 图 17 Ma=30时24°压缩拐角斜面上压力峰值点处N、NO、O质量分数沿法向分布 Fig. 17 Normal distribution of CN, CNO, and CO at pmax at Ma=30 |
|
| 图 18 Ma=30时24°压缩拐角斜面上压力峰值点处CNO+、CN2+质量分数沿法向分布 Fig. 18 Normal distribution of CNO+and CN2+ at pmax at Ma=30 |
4 结论
本文对15°、18°、24°压缩拐角模型,在自由流Ma=10~30,总焓值6~55MJ/kg范围,采用16组元双温度模型进行了无烧蚀和石墨烧蚀条件下的流场计算,得到如下结论:
1) 压缩拐角处主激波后升压在壁面附近的边界层内前传,一定条件下逆压梯度可能导致平板边界层分离。随着压缩拐角度数增大激波增强,分离可能性增大。虽然自由流马赫数升高导致边界层增厚,但压缩面的主激波角变小,激波更贴近物面,分离可能性减小。本文15°和18°压缩拐角流动未出现分离,而存在分离的24°压缩拐角流动,随着马赫数升高分离区减小。
2) 相对于低壁温条件,无烧蚀的辐射平衡壁温和壁面烧蚀条件下壁面热流有所降低。但流动的激波层增厚,分离区增大,这是高壁温导致边界层增厚从而影响激波形状及其与边界层相互作用特性的结果。分离特性的不同还将进一步影响到下游的流动特性,斜面上压力、摩阻和热流峰值点均有所后移。可见对这类存在较强激波-边界层相互作用的流动,壁面条件对流动特性的影响要远大于球头、球锥等外形。
3) 壁面烧蚀产物扩散到边界层后进一步影响边界层内的化学反应和扩散。相对于无烧蚀的辐射平衡壁温条件,烧蚀条件下边界层内高温空气的中性分子组元含量有所升高,原子和离子组元含量则有所下降。不过这两种壁面条件下的流场结构和边界层外流动特性没有明显差异,烧蚀的影响基本局限于边界层内。
本文针对准定常烧蚀流场,亦即壁面温度最高、烧蚀产物含量最大的情况下开展烧蚀对流场特性的影响研究,比较具有代表性。考虑到不同壁温条件下烧蚀产物组成不同,非定常过渡过程中烧蚀对流场特性的影响会有不同表现。进一步的工作需要耦合固体热传导方程,针对非定常烧蚀流场开展详细研究。
| [1] |
Le Jialing, Gao Tiesuo, Zeng Xuejun.
Reentry physics[M]. Beijing: National Defense Industry Press, 2005 .
(in Chinese) 乐嘉陵, 高铁锁, 曾学军. 再入物理[M]. 北京: 国防工业出版社, 2005 . |
| [2] | Keenan J A, Candler G V. Simulation of ablation in Earth atmosphere entry. AIAA-93-2789[R]. Reston:AIAA, 1993. doi:10.2514/6.1993-2789. |
| [3] | Keenan J A, Candler G V. Simulation of graphite sublimation and oxidation under re-entry conditions. AIAA-94-2083[R]. Reston:AIAA, 1994. doi:10.2514/6.1994-2083. |
| [4] | Suzuki K, Kubot H, Fujita K, et al. Chemical nonequilibrium ablation analysis of MUSES-C super-orbital reentry capsule. AIAA-97-2481[R]. Reston:AIAA, 1997. doi:10.2514/6.1997-2481. |
| [5] |
Gao T S, Dong W Z, Zhang Q Y. The computation and analysis for the hypersonic flow over reentry vehicles with ablation[J].
Acta Aerodynamica Sinica, 2006, 24(1):41–45.
(in Chinese) 高铁锁, 董维中, 张巧芸. 高超声速再入体烧蚀流场计算分析[J]. 空气动力学学报, 2006, 24(1) : 41–45. |
| [6] | Tissera S, Titarev V, Drikakis D. Chemically reacting flows around a double-cone including ablation effects. AIAA 2010-1285[R]. Reston:AIAA, 2010. doi:10.2514/6.2010-1285. |
| [7] | Bianchi D, Nasuti F, Onofri M. Aerothermodynamic analysis of reentry flows with coupled ablation. AIAA 2011-2273[R]. Reston:AIAA, 2011. doi:10.2514/6.2011-2273. |
| [8] |
Zhang W, Zeng M, Xiao L F, et al. Numerical study for the effects of ablation and pyrolysis on the hypersonic reentry flow[J].
Journal of National University of Defense Technology, 2014, 36(4):41–48.
(in Chinese) 张威, 曾明, 肖凌飞. 碳-酚醛材料烧蚀热解对再入流场特性影响的数值研究[J]. 国防科技大学学报, 2014, 36(4) : 41–48. DOI:10.11887/j.cn.20140400 |
| [9] | Mallinson S G, Gai S L, Mudford N R. High enthalpy, hypersonic compression corner flow[J]. AIAA Journal, 1996, 34(6):1130–1139. DOI:10.2514/3.13202 |
| [10] | Park C, Jaffe R L, Partridge H. Chemical-kinetic parameters of hyperbolic Earth entry[J]. Journal of Thermophysics and Heat Transfer, 2001, 15(1):76–90. DOI:10.2514/2.6582 |
| [11] | Blottner F G. Prediction of electron density in the boundary layer on entry vehicles with ablation[R]. N71-21113, 1971. |
| [12] | Kim K H, Kim C, Rho O H. Accurate computations of hypersonic flows using AUSMPW+ scheme and shock-aligned grid tehnique. AIAA-98-2442[R]. Reston:AIAA, 1998. doi:10.2514/6.1998-2442. |
| [13] | Stoll P, Gerlinger P, Bruggermann D. Domain decomposition for an implicit LU-SGS scheme using overlapping grids. AIAA-97-0770[R]. Reston:AIAA, 1997. doi:10.2514/6.1997-0770. |
| [14] |
Zeng M. Numerical rebuilding of free-stream measurement and analysis of nonequilibrium effects in high-enthalpy tunnel[D]. Beijing:Institute of Machanics, CAS, 2007. (in Chinese) 曾明.高焓风洞流场测量的数值重建与非平衡效应的数值分析[D].北京:中国科学院力学研究所, 2007. http://www.irgrid.ac.cn/handle/1471x/8637?mode=full&submit_simple=Show+full+item+record |
| [15] | Holden M S, Moselle J. Theoretical and experimental studies of the shock wave-boundary layer interaction on compression surfaces in hypersonic flow[R]. ARL 70-0002, 1970. |
| [16] | Rudy D H, Thomas J L, Kumar A. Computation of laminar viscous-inviscid interactions in high-speed internal flows[R]. N 91-21087, 1991. |
| [17] | Wang D F. Study of hypersonic viscous intereaction including high-temperature gas effects[D]. Changsha:National University of Defense Technology. 2011. (in Chinese) |


