主题052：再入飞行器热防护辐射

摘要

再入飞行器从太空返回地球大气层时，由于高超声速飞行产生的气动加热，表面温度可达数千摄氏度，必须采用有效的热防护系统来保护飞行器结构和内部设备。本主题从辐射换热的基本原理出发，系统讲解再入飞行器热防护系统的辐射换热机制，包括高超声速气动加热原理、热防护材料的热辐射特性、烧蚀热防护机制、辐射冷却与主动冷却技术等。通过建立再入飞行器热环境模型，分析不同再入轨迹、不同飞行阶段的热流密度分布，设计热防护系统方案。结合Python仿真程序，模拟再入飞行过程中的表面温度变化、热流密度分布和热防护层厚度设计，为再入飞行器热防护系统设计提供科学依据。

关键词： 再入飞行器、气动加热、热防护、烧蚀、辐射冷却、高超声速

---

1. 再入飞行热环境概述

1.1 再入飞行的基本概念

再入飞行是指航天器从太空返回地球大气层并安全着陆的过程。这一过程涉及复杂的物理现象，其中热防护是最关键的技术挑战之一。

再入飞行的主要阶段：

（1）离轨阶段：

• 航天器通过反推发动机减速，脱离原轨道
• 进入再入走廊（Entry Corridor）
• 再入角通常控制在1°-2°范围内

（2）高超声速飞行阶段：

• 速度：约7-8 km/s（轨道速度）降至2-3 km/s
• 高度：120 km降至40 km
• 产生强烈的气动加热

（3）跨声速和亚声速阶段：

• 速度：降至声速以下
• 高度：40 km降至地面
• 气动加热减弱，可采用降落伞减速

（4）着陆阶段：

• 速度：降至安全着陆速度
• 着陆方式：伞降、滑翔、垂直着陆等

1.2 高超声速气动加热原理

再入飞行器面临的主要热载荷来自高超声速飞行时的气动加热。这种加热由以下机制产生：

（1）激波加热：

当飞行器以高超声速飞行时，前方形成强烈的弓形激波。激波将动能转化为热能，使气体温度急剧升高。

激波后温度（理想气体）：

$$T_2=T_1\cdot \frac{[2\gamma M_1^2-(\gamma -1)][(\gamma -1)M_1^2+2]}{(\gamma +1{)}^2{M}_1^2}$$

其中：

• $T_1$：激波前温度
• $M_1$：飞行马赫数
• $\gamma$：比热比（空气约1.4）

示例：

当$M_1=25$，$T_1=220$ K（高空大气温度）：

$$T_2\approx 220\times \frac{[2\times 1.4\times 625-0.4][0.4\times 625+2]}{5.76\times 625}$$

$$T_2\approx 220\times \frac{1749.6\times 252}{3600}\approx 27000\text{ K}$$

实际温度由于离解和电离效应会低于此值，但仍可达8000-10000 K。

（2）粘性耗散加热：

在边界层内，气体粘性导致动能转化为热能。边界层内的温度分布可用能量方程描述：

$$\rho c_{p}u\frac{\partial T}{\partial x}=\frac{\partial }{\partial y}\left(k\frac{\partial T}{\partial y}\right)+\mu {\left(\frac{\partial u}{\partial y}\right)}^2$$

其中最后一项代表粘性耗散。

（3）热流密度计算：

再入飞行器表面的热流密度是热防护设计的关键参数。常用的热流密度计算公式包括：

Fay-Riddell公式（驻点热流）：

$$q_w=0.76P{r}^{-0.6}({\rho }_e{\mu }_e{)}^{0.5}(\frac{d{u}_e}{dx}{)}^{0.5}(h_{aw}-h_w)$$

其中：

• $Pr$：普朗特数
• ${\rho }_e$：边界层外缘密度
• ${\mu }_e$：边界层外缘粘度
• $d{u}_e/dx$：速度梯度
• $h_{aw}$：绝热壁焓
• $h_w$：壁面焓

简化公式（工程估算）：

$$q_w=k\sqrt{\frac{\rho }{R_n}}{V}^3$$

其中：

• $k$：经验常数（约$1.83\times 10^{-4}$）
• $\rho$：大气密度
• $R_n$：鼻锥半径
• $V$：飞行速度

1.3 再入热环境特征

热流密度分布：

再入过程中，热流密度随时间和空间变化：

时间分布：

• 再入初期：热流密度逐渐增大
• 最大热流点：通常发生在30-50 km高度
• 再入后期：热流密度逐渐减小

空间分布：

• 驻点（Stagnation Point）：热流密度最大
• 迎风面：热流密度较高
• 背风面：热流密度较低
• 侧面：中等热流密度

典型再入热环境参数：

参数	载人飞船	航天飞机	再入弹头
再入速度	7.8 km/s	7.5 km/s	6.5 km/s
最大热流密度	300-500 W/cm²	100-150 W/cm²	500-1000 W/cm²
总加热量	10-20 MJ/kg	30-50 MJ/kg	5-10 MJ/kg
峰值温度	3000-4000 K	1500-2000 K	4000-5000 K

---

2. 热防护系统类型

2.1 被动热防护系统

被动热防护系统依靠材料本身的热容和隔热性能来吸收和阻隔热量。

（1）吸热式热防护（Heat Sink）：

利用高热容材料吸收热量，通过温升储存热能。

原理：
$$Q=m{c}_p\Delta T$$

其中：

• $m$：材料质量
• $c_p$：比热容
• $\Delta T$：温升

常用材料：

• 铜：$c_p=385$ J/kg·K，$\rho =8960$ kg/m³
• 铍：$c_p=1825$ J/kg·K，$\rho =1850$ kg/m³
• 石墨：$c_p=710$ J/kg·K，$\rho =2200$ kg/m³

特点：

• 结构简单、可靠性高
• 质量大，仅适用于短时加热
• 适用于再入弹头、小尺寸部件

（2）辐射式热防护（Radiative Cooling）：

利用高温表面的辐射散热来平衡气动加热。

原理：
$$q_{rad}=\epsilon \sigma T_w^4$$

当$q_{rad}=q_{aero}$时，达到平衡温度：

$$T_w={\left(\frac{q_{aero}}{\epsilon \sigma }\right)}^{1/4}$$

示例：

若$q_{aero}=100$ W/cm² = $10^6$ W/m²，$\epsilon =0.8$：

$$T_w={\left(\frac{10^6}{0.8\times 5.67\times 10^{-8}}\right)}^{1/4}\approx 2100\text{ K}$$

常用材料：

• 难熔金属（钨、钼、钽）：熔点>3000 K
• 陶瓷（碳化硅、氮化硅）：耐高温、低导热
• 碳-碳复合材料：高温强度好

特点：

• 可重复使用
• 需要耐高温材料
• 适用于航天飞机、可重复使用运载器

（3）烧蚀式热防护（Ablative Thermal Protection）：

通过材料的烧蚀（熔化、升华、分解）吸收热量，并形成隔热层保护内部结构。

烧蚀吸热机制：

1. 相变吸热： 材料熔化、升华吸收潜热
2. 热解吸热： 有机材料分解吸热
3. 质量引射： 分解气体向外喷射，阻挡热流
4. 炭化层隔热： 形成多孔炭化层，降低导热

常用材料：

• AVCOAT（阿波罗飞船）：环氧-酚醛树脂+二氧化硅微球
• PICA（星尘号）：酚醛浸渍碳烧蚀体
• SLA（航天飞机）：硅基低密度烧蚀材料
• 碳-酚醛：高密度烧蚀材料

烧蚀材料性能参数：

材料	密度 (kg/m³)	烧蚀热 (MJ/kg)	最高使用温度 (K)
AVCOAT	500-600	10-15	3500
PICA	250-300	15-20	3800
碳-酚醛	1400-1600	8-12	4000
硅橡胶	400-500	5-8	3000

2.2 主动热防护系统

主动热防护系统通过工质流动带走热量。

（1） transpiration cooling（发汗冷却）：

通过多孔壁面向边界层注入冷却剂（气体或液体），形成隔热层。

原理：

• 冷却剂吸收壁面热量
• 在边界层内形成低温膜层
• 改变边界层温度分布

冷却效率：
$$\eta =\frac{q_{nocooling}-q_{withcooling}}{q_{nocooling}}$$

典型效率可达80-95%。

（2）薄膜冷却（Film Cooling）：

在壁面切向注入冷却剂，形成保护性气膜。

与发汗冷却的区别：

• 发汗冷却：通过多孔壁面垂直注入
• 薄膜冷却：通过缝隙切向注入

（3）再生冷却（Regenerative Cooling）：

冷却剂在壁面内部的冷却通道中流动，吸收热量后用于推进系统或其他用途。

常用于：

• 液体火箭发动机燃烧室
• 高超声速飞行器前缘

2.3 热防护系统的选择

热防护系统的选择取决于以下因素：

（1）热环境 severity：

• 低热流密度（<100 W/cm²）：辐射式、吸热式
• 中热流密度（100-500 W/cm²）：烧蚀式、主动冷却
• 高热流密度（>500 W/cm²）：烧蚀式、发汗冷却

（2）任务类型：

• 一次性任务：烧蚀式
• 可重复使用：辐射式、主动冷却

（3）质量约束：

• 烧蚀式：质量效率低（单位面积质量大）
• 辐射式：质量效率高
• 主动冷却：需要额外冷却剂质量

（4）成本考虑：

• 烧蚀式：成本低，不可重复使用
• 辐射式：成本高，可重复使用

---

3. 烧蚀热防护机理

3.1 烧蚀过程的热物理模型

烧蚀是一个复杂的热-化学-力学耦合过程，涉及多种物理现象：

（1）热解区（Pyrolysis Zone）：

温度升高导致有机粘结剂分解，产生气体和炭残留物。

热解反应：
$$\text{聚合物}\to \text{炭}+\text{热解气体}$$

热解速率可用Arrhenius方程描述：

$$\frac{d\alpha }{dt}=A(1-\alpha {)}^n\exp \left(-\frac{E_a}{RT}\right)$$

其中：

• $\alpha$：热解转化率
• $A$：指前因子
• $E_a$：活化能
• $n$：反应级数

（2）炭化区（Charring Zone）：

热解后形成的多孔炭化层具有以下特性：

• 高孔隙率（50-80%）
• 低导热系数（0.1-0.5 W/m·K）
• 高辐射发射率（0.8-0.95）

炭化层的隔热作用：
$$q={\lambda }_{char}\frac{T_{surface}-T_{pyrolysis}}{{\delta }_{char}}$$

（3）表面烧蚀（Surface Ablation）：

表面材料通过以下机制被移除：

• 熔化： 无机填料（二氧化硅）熔化并被气流剪切带走
• 升华： 碳材料直接升华
• 氧化： 表面碳与边界层中的氧气反应
• 机械剥蚀： 气流剪切力导致材料剥落

3.2 烧蚀质量损失计算

质量烧蚀率：

$$\dot{m}=\frac{q_{net}}{Q^{\ast }}$$

其中：

• $\dot{m}$：质量烧蚀率（kg/m²·s）
• $q_{net}$：净热流密度（W/m²）
• $Q^{\ast }$：有效烧蚀热（J/kg）

有效烧蚀热：

有效烧蚀热是衡量烧蚀材料性能的关键参数，包括：

$$Q^{\ast }=Q_{phase}+Q_{pyrolysis}+Q_{gas}+Q_{blockage}$$

其中：

• $Q_{phase}$：相变潜热（熔化、升华）
• $Q_{pyrolysis}$：热解吸热
• $Q_{gas}$：热解气体升温吸热
• $Q_{blockage}$：质量引射阻塞效应

典型烧蚀材料的$Q^{\ast }$值：

• AVCOAT：10-15 MJ/kg
• PICA：15-20 MJ/kg
• 碳-碳：20-30 MJ/kg

3.3 炭化层厚度增长模型

炭化层厚度随时间增长：

$$\frac{d{\delta }_{char}}{dt}=\frac{\dot{m}}{{\rho }_{char}}$$

其中：

• ${\delta }_{char}$：炭化层厚度
• ${\rho }_{char}$：炭化层密度

热防护层厚度设计：

所需热防护层厚度：

$${\delta }_{TPS}={\delta }_{char}+{\delta }_{virgin}+{\delta }_{safety}$$

其中：

• ${\delta }_{char}$：炭化层厚度（随时间增长）
• ${\delta }_{virgin}$：原始材料厚度（隔热用）
• ${\delta }_{safety}$：安全裕度（通常20-30%）

---

4. 辐射冷却设计

4.1 辐射平衡温度

对于辐射式热防护，表面温度由辐射平衡决定：

$$q_{aero}=\epsilon \sigma T_w^4+q_{conduction}$$

在稳态条件下（$q_{conduction}\approx 0$）：

$$T_w={\left(\frac{q_{aero}}{\epsilon \sigma }\right)}^{1/4}$$

温度限制：

辐射冷却的表面温度受材料耐温限制：

• 镍基超合金：约1300 K
• 铌合金：约1700 K
• 碳-碳复合材料：约2300 K
• 陶瓷（ZrB₂、HfC）：>2500 K

4.2 辐射冷却材料

（1）难熔金属：

材料	熔点 (K)	密度 (kg/m³)	发射率
钨 (W)	3695	19300	0.4-0.6
钼 (Mo)	2896	10200	0.3-0.5
钽 (Ta)	3290	16600	0.3-0.5
铌 (Nb)	2750	8570	0.3-0.4

（2）超高温陶瓷（UHTCs）：

材料	熔点 (K)	密度 (kg/m³)	发射率
ZrB₂	3310	6090	0.7-0.85
HfB₂	3380	11200	0.7-0.85
HfC	4160	12400	0.5-0.7
TaC	4270	14600	0.4-0.6

（3）碳-碳复合材料：

特点：

• 高温强度保持率高
• 热导率低（2-10 W/m·K）
• 发射率：0.8-0.9
• 最高使用温度：2300 K（氧化环境）/ 3000 K（惰性环境）

4.3 辐射散热器设计

对于需要排散内部热量的辐射散热器：

散热功率：

$$Q_{rad}=\epsilon \sigma A{T}^4$$

散热器面积：

$$A=\frac{Q_{internal}}{\epsilon \sigma T^4}$$

设计考虑：

• 散热器应避开高气动加热区
• 可采用可展开散热器增大面积
• 表面发射率优化（高发射率涂层）

---

5. 热防护系统设计方法

5.1 热分析流程

热防护系统设计的一般流程：

（1）确定热环境：

• 再入轨迹分析
• 气动加热计算
• 热流密度分布

（2）选择热防护类型：

• 根据热环境 severity
• 考虑任务要求（一次性/可重复使用）
• 质量、成本约束

（3）材料选择：

• 耐温性能
• 烧蚀性能（如适用）
• 力学性能
• 工艺性

（4）厚度设计：

• 热传导分析
• 烧蚀深度计算（如适用）
• 安全裕度

（5）验证试验：

• 等离子体风洞试验
• 电弧加热器试验
• 飞行试验

5.2 热传导分析

一维热传导方程：

$$\rho c_p\frac{\partial T}{\partial t}=\frac{\partial }{\partial x}\left(k\frac{\partial T}{\partial x}\right)$$

边界条件：

• 表面：$-k\frac{\partial T}{\partial x}{|}_{x=0}=q_{aero}-\epsilon \sigma T_s^4$
• 背面：绝热或对流换热

数值求解：

采用有限差分法或有限元法求解。对于烧蚀问题，需要考虑移动边界。

5.3 典型热防护系统设计

（1）载人飞船返回舱：

热环境：

• 最大热流密度：300-500 W/cm²
• 总加热量：10-20 MJ/kg

设计方案：

• 底部：AVCOAT烧蚀材料，厚度50-70 mm
• 侧面：AVCOAT，厚度30-40 mm
• 背部：低密度烧蚀材料，厚度20-30 mm

（2）航天飞机轨道器：

热环境：

• 最大热流密度：100-150 W/cm²
• 总加热量：30-50 MJ/kg

设计方案：

• 鼻锥和翼前缘：碳-碳复合材料
• 底部：高温可重复使用表面隔热瓦（HRSI，Li₂SiO₃）
• 侧面：低温可重复使用表面隔热瓦（LRSI）
• 背部：柔性可重复使用表面隔热毡（FRSI）

（3）高超声速飞行器：

热环境：

• 前缘热流密度：500-1000 W/cm²
• 持续时间：数百秒至数千秒

设计方案：

• 前缘：主动冷却（再生冷却或发汗冷却）
• 机身：辐射式热防护（超高温陶瓷或碳-碳）

---

6. Python仿真案例分析

6.1 案例1：再入气动加热计算

本案例计算不同再入轨迹下的气动加热环境。

模拟内容：

1. 计算再入过程中的速度、高度变化
2. 计算热流密度随时间变化
3. 分析不同再入角的影响
4. 计算总加热量

关键公式：

再入动力学：

$$\frac{dV}{dt}=-\frac{C_{D}A}{2m}\rho V^2-g\sin \gamma$$

$$\frac{d\gamma }{dt}=\left(\frac{C_{L}A}{2m}\rho V-\frac{g}{V}\right)\cos \gamma$$

$$\frac{dh}{dt}=-V\sin \gamma$$

其中：

• $V$：速度
• $\gamma$：飞行路径角
• $h$：高度
• $C_D$：阻力系数
• $C_L$：升力系数
• $\rho$：大气密度

热流密度：

$$q=k\sqrt{\frac{\rho }{R_n}}{V}^3$$

模拟结果分析：

弹道式再入（无升力）：

• 再入角：1.5°
• 最大热流密度：约400 W/cm²
• 最大热流高度：约50 km
• 总加热量：约15 MJ/kg

滑翔式再入（有升力）：

• 升阻比：1.0
• 最大热流密度：约150 W/cm²
• 再入时间延长，但总加热量增加

6.2 案例2：烧蚀热防护计算

本案例模拟烧蚀热防护层的温度和烧蚀深度。

模拟内容：

1. 计算热防护层内温度分布
2. 计算炭化层厚度增长
3. 计算表面烧蚀率
4. 评估热防护层厚度是否足够

模拟方法：

一维热传导（含烧蚀）：

$$\rho c_p\frac{\partial T}{\partial t}=\frac{\partial }{\partial x}\left(k\frac{\partial T}{\partial x}\right)+\dot{m}{c}_{pg}\frac{\partial T}{\partial x}$$

其中最后一项代表热解气体的对流换热。

炭化层厚度：

$${\delta }_{char}(t)={\int }_0^t\frac{\dot{m}(\tau )}{{\rho }_{char}}d\tau$$

模拟结果分析：

AVCOAT材料（厚度50 mm）：

• 表面温度：2500-3000 K
• 背面温度：<400 K（安全）
• 炭化层厚度：15-20 mm
• 剩余厚度：30-35 mm（安全裕度充足）

6.3 案例3：辐射冷却分析

本案例分析辐射式热防护的性能。

模拟内容：

1. 计算不同材料的最大允许热流密度
2. 分析辐射平衡温度
3. 计算向内部结构的热传导
4. 优化材料厚度和发射率

模拟结果分析：

碳-碳复合材料：

• 最高使用温度：2300 K
• 对应最大热流密度（ε=0.85）：
  $$q_{max}=0.85\times 5.67\times 10^{-8}\times 2300^4\approx 150{\text{ W/cm}}^2$$

超高温陶瓷（ZrB₂）：

• 最高使用温度：2500 K
• 对应最大热流密度（ε=0.8）：
  $$q_{max}=0.8\times 5.67\times 10^{-8}\times 2500^4\approx 220{\text{ W/cm}}^2$$

6.4 案例4：热防护系统综合设计

本案例综合设计一个再入飞行器的热防护系统。

模拟内容：

1. 根据再入轨迹确定热环境
2. 选择热防护类型和材料
3. 计算热防护层厚度
4. 验证设计是否满足要求

设计案例：

飞行器参数：

• 质量：5000 kg
• 底部直径：3 m
• 再入速度：7.5 km/s
• 再入角：1.5°

热环境：

• 最大热流密度：350 W/cm²
• 总加热量：18 MJ/kg
• 再入时间：约600 s

设计方案：

• 底部：AVCOAT，厚度60 mm
• 侧面：AVCOAT，厚度40 mm
• 背部：低密度烧蚀材料，厚度25 mm

验证结果：

• 背面温度：<350 K（满足要求）
• 剩余厚度：>50%（安全裕度充足）
• 总质量：约800 kg（可接受）

---

7. 结果讨论与工程应用

7.1 关键发现总结

通过上述仿真分析，可以得出以下关键结论：

（1）再入热环境极其严酷

再入飞行器面临的热流密度可达数百至上千W/cm²，表面温度可达数千摄氏度。热防护系统是再入飞行器设计的关键。

（2）烧蚀式热防护适用于高焓值任务

烧蚀材料通过相变、热解、质量引射等机制有效吸收热量，适用于载人飞船、再入弹头等高焓值任务。

（3）辐射式热防护适用于可重复使用运载器

辐射冷却系统可重复使用，适用于航天飞机、高超声速飞行器等任务，但受材料耐温限制。

（4）热防护设计需要综合考虑

热防护系统的选择需要在热性能、质量、成本、可靠性之间权衡。没有一种热防护系统适用于所有任务。

7.2 热防护系统设计要点

基于仿真分析，总结热防护系统设计的关键要点：

（1）准确预测热环境是基础

热环境预测的准确性直接影响热防护设计的可靠性。需要考虑气动加热、辐射加热、催化效应等多种因素。

（2）材料选择至关重要

热防护材料的耐温性能、烧蚀性能、力学性能直接影响热防护系统的有效性。需要根据具体任务选择合适材料。

（3）安全裕度必须充足

热防护设计必须留有足够的安全裕度（通常20-30%），以应对不确定性因素（如再入轨迹偏差、材料性能分散性等）。

（4）地面试验验证必不可少

热防护系统必须通过等离子体风洞、电弧加热器等地面设备验证，才能用于飞行任务。

7.3 新技术发展趋势

热防护技术正在向以下方向发展：

新型烧蚀材料：

• 3D编织碳-碳复合材料
• 纳米改性烧蚀材料
• 可编程烧蚀材料（智能响应）

可重复使用热防护：

• 超高温陶瓷基复合材料（CMCs）
• 难熔金属合金
• 陶瓷-金属梯度材料

主动热防护：

• 微通道冷却技术
• 相变材料储能
• 磁流体动力学（MHD）热防护

多功能热防护：

• 结构/热防护一体化
• 隐身/热防护一体化
• 能量收集/热防护一体化

附录：Python仿真程序说明

本主题配套的Python仿真程序包含以下功能模块：

A.1 程序结构

run_simulation.py
├── 案例1：再入气动加热计算
│   ├── 再入轨迹计算
│   ├── 热流密度计算
│   └── 不同再入角对比
├── 案例2：烧蚀热防护计算
│   ├── 温度分布计算
│   ├── 炭化层厚度增长
│   └── 烧蚀深度预测
├── 案例3：辐射冷却分析
│   ├── 辐射平衡温度
│   ├── 材料性能对比
│   └── 最大热流密度
└── 案例4：热防护系统综合设计
    ├── 热环境分析
    ├── 材料选择
    └── 厚度设计验证

A.2 关键算法

再入轨迹计算：

def reentry_dynamics(V0, gamma0, h0, m, CD, A, dt=0.1):
    """计算再入轨迹"""
    # 积分运动方程
    # 返回速度、高度、路径角随时间变化

热流密度计算：

def heat_flux(V, rho, Rn, k=1.83e-4):
    """计算驻点热流密度"""
    return k * np.sqrt(rho / Rn) * V**3

烧蚀深度计算：

def ablation_depth(q, Q_star, rho, t):
    """计算烧蚀深度"""
    m_dot = q / Q_star  # 质量烧蚀率
    return (m_dot / rho) * t  # 烧蚀深度

A.3 可视化输出

程序生成以下可视化结果：

1. 再入轨迹图：高度-速度、高度-时间曲线
2. 热流密度图：热流密度随时间/高度变化
3. 温度分布图：热防护层内温度分布
4. 烧蚀分析图：炭化层厚度、烧蚀深度随时间变化
5. 材料对比图：不同材料的性能对比

A.4 使用说明

运行程序前请确保安装以下依赖：

pip install numpy matplotlib scipy

运行程序：

python run_simulation.py

"""
主题052：再入飞行器热防护辐射 - Python仿真程序
================================================================================
本程序包含以下案例：
  1. 再入气动加热计算
  2. 烧蚀热防护计算
  3. 辐射冷却分析
  4. 热防护系统综合设计

作者：仿真教学团队
日期：2026年
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, PillowWriter
import matplotlib
matplotlib.use('Agg')  # 使用非交互式后端

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

# 物理常数
SIGMA = 5.67e-8  # 斯蒂芬-玻尔兹曼常数 W/m²·K⁴
R_E = 6371e3  # 地球半径 m
g0 = 9.81  # 重力加速度 m/s²

print("=" * 80)
print("主题052：再入飞行器热防护辐射 - Python仿真")
print("=" * 80)
print("\n本程序包含以下案例：")
print("  1. 再入气动加热计算")
print("  2. 烧蚀热防护计算")
print("  3. 辐射冷却分析")
print("  4. 热防护系统综合设计")
print("\n" + "=" * 80)


def atmosphere_density(h):
    """
    计算大气密度 (指数模型)
    
    参数:
    - h: 高度 (m)
    
    返回:
    - rho: 大气密度 (kg/m³)
    """
    rho0 = 1.225  # 海平面密度 kg/m³
    H = 8500  # 标高 m
    rho = rho0 * np.exp(-h / H)
    return rho


def atmosphere_temperature(h):
    """
    计算大气温度 (简化模型)
    
    参数:
    - h: 高度 (m)
    
    返回:
    - T: 大气温度 (K)
    """
    # 简化模型：对流层温度递减
    T0 = 288.15  # 海平面温度 K
    lapse_rate = 0.0065  # 温度递减率 K/m
    T = T0 - lapse_rate * h
    return max(T, 220)  # 最低温度限制


def reentry_trajectory(V0, gamma0, h0, m, CD, A, CL=0, dt=1.0, t_max=2000):
    """
    计算再入轨迹
    
    参数:
    - V0: 初始速度 (m/s)
    - gamma0: 初始飞行路径角 (rad, 负值表示下降)
    - h0: 初始高度 (m)
    - m: 飞行器质量 (kg)
    - CD: 阻力系数
    - A: 参考面积 (m²)
    - CL: 升力系数 (默认0，弹道式再入)
    - dt: 时间步长 (s)
    - t_max: 最大计算时间 (s)
    
    返回:
    - t: 时间数组 (s)
    - V: 速度数组 (m/s)
    - h: 高度数组 (m)
    - gamma: 飞行路径角数组 (rad)
    - rho: 大气密度数组 (kg/m³)
    """
    # 初始化
    t = [0]
    V = [V0]
    h = [h0]
    gamma = [gamma0]
    rho = [atmosphere_density(h0)]
    
    # 数值积分
    i = 0
    while h[i] > 0 and t[i] < t_max:
        # 当前状态
        Vi = V[i]
        hi = h[i]
        gammai = gamma[i]
        rhoi = atmosphere_density(hi)
        
        # 重力加速度 (随高度变化)
        gi = g0 * (R_E / (R_E + hi))**2
        
        # 气动力
        D = 0.5 * rhoi * Vi**2 * CD * A  # 阻力
        L = 0.5 * rhoi * Vi**2 * CL * A  # 升力
        
        # 运动方程
        dV_dt = -D/m - gi * np.sin(gammai)
        dgamma_dt = (L/(m*Vi) - gi/Vi * np.cos(gammai)) * np.cos(gammai)
        dh_dt = -Vi * np.sin(gammai)
        
        # 更新状态
        V.append(Vi + dV_dt * dt)
        gamma.append(gammai + dgamma_dt * dt)
        h.append(hi + dh_dt * dt)
        t.append(t[i] + dt)
        rho.append(rhoi)
        
        i += 1
        
        # 防止无限循环
        if i > 100000:
            break
    
    return np.array(t), np.array(V), np.array(h), np.array(gamma), np.array(rho)


def stagnation_heat_flux(V, rho, Rn, k=1.83e-4):
    """
    计算驻点热流密度 (Sutton-Graves公式)
    
    参数:
    - V: 速度 (m/s)
    - rho: 大气密度 (kg/m³)
    - Rn: 鼻锥半径 (m)
    - k: 经验常数
    
    返回:
    - q: 热流密度 (W/m²)
    """
    q = k * np.sqrt(rho / Rn) * V**3
    return q


def radiation_equilibrium_temperature(q, epsilon):
    """
    计算辐射平衡温度
    
    参数:
    - q: 热流密度 (W/m²)
    - epsilon: 表面发射率
    
    返回:
    - T: 辐射平衡温度 (K)
    """
    T = (q / (epsilon * SIGMA))**0.25
    return T


# ==============================================================================
# 案例1：再入气动加热计算
# ==============================================================================
print("\n" + "=" * 80)
print("案例1：再入气动加热计算")
print("=" * 80)

# 再入参数
V0 = 7500  # 初始速度 m/s (约7.5 km/s)
h0 = 120e3  # 初始高度 m (120 km)
m = 5000  # 质量 kg
CD = 1.5  # 阻力系数
A = 10  # 参考面积 m² (底部直径约3.6m)
Rn = 1.0  # 鼻锥半径 m

# 不同再入角
gamma_angles = [-1.0, -1.5, -2.0]  # 度

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

results = {}
for gamma_deg in gamma_angles:
    gamma0 = np.radians(gamma_deg)
    
    print(f"\n【再入角 {gamma_deg}° 计算中...】")
    
    # 计算再入轨迹
    t, V, h, gamma, rho = reentry_trajectory(V0, gamma0, h0, m, CD, A, dt=2.0)
    
    # 计算热流密度
    q = stagnation_heat_flux(V, rho, Rn)
    
    # 存储结果
    results[gamma_deg] = {
        't': t, 'V': V, 'h': h, 'q': q, 'rho': rho
    }
    
    # 找到最大热流
    q_max_idx = np.argmax(q)
    q_max = q[q_max_idx]
    h_qmax = h[q_max_idx]
    t_qmax = t[q_max_idx]
    
    # 计算总加热量 (积分)
    Q_total = np.trapezoid(q, t)  # J/m²
    
    print(f"  再入时间: {t[-1]:.1f} s")
    print(f"  最大热流密度: {q_max/1e4:.1f} W/cm² = {q_max/1e6:.2f} MW/m²")
    print(f"  最大热流高度: {h_qmax/1e3:.1f} km")
    print(f"  总加热量: {Q_total/1e6:.1f} MJ/m²")

# 1.1 高度-速度曲线
ax1 = axes[0, 0]
for gamma_deg, data in results.items():
    ax1.plot(data['V']/1e3, data['h']/1e3, '-', label=f'γ={gamma_deg}°', linewidth=2)
ax1.set_xlabel('速度 (km/s)', fontsize=11)
ax1.set_ylabel('高度 (km)', fontsize=11)
ax1.set_title('再入轨迹：高度-速度曲线', fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(alpha=0.3)
ax1.invert_yaxis()

# 1.2 热流密度随时间变化
ax2 = axes[0, 1]
for gamma_deg, data in results.items():
    ax2.plot(data['t'], data['q']/1e4, '-', label=f'γ={gamma_deg}°', linewidth=2)
    # 标注最大热流点
    q_max_idx = np.argmax(data['q'])
    ax2.plot(data['t'][q_max_idx], data['q'][q_max_idx]/1e4, 'o', markersize=8)
ax2.set_xlabel('时间 (s)', fontsize=11)
ax2.set_ylabel('热流密度 (W/cm²)', fontsize=11)
ax2.set_title('驻点热流密度随时间变化', fontsize=12, fontweight='bold')
ax2.legend()
ax2.grid(alpha=0.3)

# 1.3 热流密度随高度变化
ax3 = axes[1, 0]
for gamma_deg, data in results.items():
    ax3.semilogy(data['h']/1e3, data['q']/1e4, '-', label=f'γ={gamma_deg}°', linewidth=2)
ax3.set_xlabel('高度 (km)', fontsize=11)
ax3.set_ylabel('热流密度 (W/cm²)', fontsize=11)
ax3.set_title('热流密度随高度变化', fontsize=12, fontweight='bold')
ax3.legend()
ax3.grid(alpha=0.3, which='both')
ax3.invert_xaxis()

# 1.4 不同鼻锥半径的影响
ax4 = axes[1, 1]
gamma_deg = -1.5  # 选择中间再入角
data = results[gamma_deg]
Rn_values = [0.5, 1.0, 2.0, 3.0]  # m

for Rn in Rn_values:
    q = stagnation_heat_flux(data['V'], data['rho'], Rn)
    ax4.plot(data['t'], q/1e4, '-', label=f'Rn={Rn}m', linewidth=2)
    q_max_idx = np.argmax(q)
    print(f"\n  鼻锥半径{Rn}m: 最大热流{q[q_max_idx]/1e4:.1f} W/cm²")

ax4.set_xlabel('时间 (s)', fontsize=11)
ax4.set_ylabel('热流密度 (W/cm²)', fontsize=11)
ax4.set_title(f'鼻锥半径对热流密度的影响 (γ={gamma_deg}°)', fontsize=12, fontweight='bold')
ax4.legend()
ax4.grid(alpha=0.3)

plt.tight_layout()
plt.savefig('case1_reentry_heating.png', dpi=150, bbox_inches='tight')
plt.close()

print("\n✓ 案例1结果已保存: case1_reentry_heating.png")


# ==============================================================================
# 案例2：烧蚀热防护计算
# ==============================================================================
print("\n" + "=" * 80)
print("案例2：烧蚀热防护计算")
print("=" * 80)

# 烧蚀材料参数 (AVCOAT类似材料)
ablative_material = {
    'name': 'AVCOAT类烧蚀材料',
    'rho': 550,  # 密度 kg/m³
    'rho_char': 300,  # 炭化层密度 kg/m³
    'cp': 1500,  # 比热容 J/kg·K
    'k': 0.15,  # 导热系数 W/m·K
    'Q_star': 12e6,  # 有效烧蚀热 J/kg
    'T_pyrolysis': 600,  # 热解温度 K
    'T_surface_max': 3500,  # 最高表面温度 K
    'epsilon': 0.9  # 表面发射率
}

# 使用案例1的再入数据 (γ=-1.5°)
data = results[-1.5]
t = data['t']
q = data['q']

# 热防护层参数
thickness_TPS = 0.06  # 热防护层厚度 m (60 mm)
n_nodes = 50  # 节点数
dx = thickness_TPS / (n_nodes - 1)

# 材料属性
rho = ablative_material['rho']
cp = ablative_material['cp']
k = ablative_material['k']
Q_star = ablative_material['Q_star']
rho_char = ablative_material['rho_char']
epsilon = ablative_material['epsilon']

# 初始化温度场
T = np.ones(n_nodes) * 300  # 初始温度 300 K
T_history = [T.copy()]
T_surface_history = [T[0]]
T_back_history = [T[-1]]
char_depth_history = [0]
ablation_depth_history = [0]

# 时间步进 (使用较小的时间步长)
dt = 0.5  # s
n_steps = len(t)

print(f"\n  正在进行烧蚀热防护计算...")
print(f"  热防护层厚度: {thickness_TPS*1000:.0f} mm")
print(f"  时间步数: {n_steps}")

# 瞬态热传导计算
for i in range(1, n_steps):
    # 当前热流
    q_current = q[i]
    
    # 辐射散热
    q_rad = epsilon * SIGMA * T[0]**4
    
    # 净热流
    q_net = max(0, q_current - q_rad)
    
    # 计算烧蚀率
    m_dot = q_net / Q_star  # kg/m²/s
    ablation_rate = m_dot / rho  # m/s
    
    # 更新炭化层厚度
    if i == 1:
        char_depth = 0
    else:
        # 简化的炭化层增长模型
        if T[0] > ablative_material['T_pyrolysis']:
            char_growth_rate = m_dot / rho_char
            char_depth = char_depth_history[-1] + char_growth_rate * dt
        else:
            char_depth = char_depth_history[-1]
    
    # 有限差分法求解热传导
    T_new = T.copy()
    
    # 表面边界条件 (考虑烧蚀)
    T_new[0] = T[0] + (2*dt/(rho*cp*dx)) * (q_net - k*(T[0]-T[1])/dx)
    
    # 内部节点
    for j in range(1, n_nodes-1):
        T_new[j] = T[j] + (k*dt/(rho*cp*dx**2)) * (T[j+1] - 2*T[j] + T[j-1])
    
    # 背面边界条件 (绝热)
    T_new[-1] = T[-1] + (2*k*dt/(rho*cp*dx**2)) * (T[-2] - T[-1])
    
    # 更新温度
    T = T_new.copy()
    
    # 限制温度 (物理约束)
    T = np.clip(T, 300, 4000)
    
    # 记录历史
    if i % 10 == 0:  # 每10步记录一次
        T_history.append(T.copy())
        T_surface_history.append(T[0])
        T_back_history.append(T[-1])
        char_depth_history.append(char_depth)
        ablation_depth_history.append(ablation_depth_history[-1] + ablation_rate * dt * 10)

T_history = np.array(T_history)
T_surface_history = np.array(T_surface_history)
T_back_history = np.array(T_back_history)
char_depth_history = np.array(char_depth_history)
ablation_depth_history = np.array(ablation_depth_history)

print(f"  计算完成")
print(f"  最高表面温度: {np.max(T_surface_history):.0f} K ({np.max(T_surface_history)-273:.0f}°C)")
print(f"  最高背面温度: {np.max(T_back_history):.0f} K ({np.max(T_back_history)-273:.0f}°C)")
print(f"  最终炭化层厚度: {char_depth_history[-1]*1000:.1f} mm")
print(f"  最终烧蚀深度: {ablation_depth_history[-1]*1000:.1f} mm")
print(f"  剩余厚度: {(thickness_TPS - ablation_depth_history[-1])*1000:.1f} mm")

# 创建可视化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 2.1 表面和背面温度随时间变化
ax1 = axes[0, 0]
t_plot = np.linspace(0, t[-1], len(T_surface_history))
ax1.plot(t_plot, T_surface_history, 'r-', label='表面温度', linewidth=2)
ax1.plot(t_plot, T_back_history, 'b-', label='背面温度', linewidth=2)
ax1.axhline(y=ablative_material['T_pyrolysis'], color='orange', linestyle='--', 
           label=f"热解温度 {ablative_material['T_pyrolysis']}K")
ax1.set_xlabel('时间 (s)', fontsize=11)
ax1.set_ylabel('温度 (K)', fontsize=11)
ax1.set_title('热防护层温度变化', fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(alpha=0.3)

# 2.2 温度分布 (不同时刻)
ax2 = axes[0, 1]
x = np.linspace(0, thickness_TPS*1000, n_nodes)  # mm

# 选择几个时刻显示
time_indices = [0, len(T_history)//4, len(T_history)//2, 3*len(T_history)//4, len(T_history)-1]
colors = plt.cm.viridis(np.linspace(0, 1, len(time_indices)))

for idx, color in zip(time_indices, colors):
    ax2.plot(x, T_history[idx], '-', color=color, 
            label=f't={t_plot[idx]:.0f}s', linewidth=2)

ax2.set_xlabel('深度 (mm)', fontsize=11)
ax2.set_ylabel('温度 (K)', fontsize=11)
ax2.set_title('热防护层内温度分布', fontsize=12, fontweight='bold')
ax2.legend()
ax2.grid(alpha=0.3)

# 2.3 炭化层和烧蚀深度随时间变化
ax3 = axes[1, 0]
ax3.plot(t_plot, char_depth_history*1000, 'g-', label='炭化层厚度', linewidth=2)
ax3.plot(t_plot, ablation_depth_history*1000, 'r-', label='烧蚀深度', linewidth=2)
ax3.axhline(y=thickness_TPS*1000, color='black', linestyle='--', label=f'初始厚度 {thickness_TPS*1000:.0f}mm')
ax3.set_xlabel('时间 (s)', fontsize=11)
ax3.set_ylabel('深度 (mm)', fontsize=11)
ax3.set_title('炭化层和烧蚀深度变化', fontsize=12, fontweight='bold')
ax3.legend()
ax3.grid(alpha=0.3)

# 2.4 不同热防护层厚度对比
ax4 = axes[1, 1]
thickness_values = [40, 50, 60, 70]  # mm

for thickness_mm in thickness_values:
    thickness = thickness_mm / 1000  # m
    
    # 简化的背面温度估算 (假设线性温度分布)
    # 实际应该重新运行热传导计算
    T_back_approx = 300 + (np.max(T_surface_history) - 300) * (1 - np.exp(-thickness_mm/30))
    
    ax4.bar(thickness_mm, T_back_approx, width=5, alpha=0.7, 
           label=f'{thickness_mm}mm')
    print(f"  厚度{thickness_mm}mm: 估算背面温度 {T_back_approx:.0f}K")

ax4.axhline(y=350+273, color='red', linestyle='--', label='允许温度上限 (350°C)')
ax4.set_xlabel('热防护层厚度 (mm)', fontsize=11)
ax4.set_ylabel('背面温度 (K)', fontsize=11)
ax4.set_title('不同厚度热防护层背面温度', fontsize=12, fontweight='bold')
ax4.legend()
ax4.grid(alpha=0.3, axis='y')

plt.tight_layout()
plt.savefig('case2_ablative_tps.png', dpi=150, bbox_inches='tight')
plt.close()

print("\n✓ 案例2结果已保存: case2_ablative_tps.png")


# ==============================================================================
# 案例3：辐射冷却分析
# ==============================================================================
print("\n" + "=" * 80)
print("案例3：辐射冷却分析")
print("=" * 80)

# 辐射冷却材料参数
radiative_materials = {
    '碳-碳复合材料': {
        'T_max': 2300,  # K
        'epsilon': 0.85,
        'k': 5,  # W/m·K
        'rho': 1800,  # kg/m³
        'color': '#2F4F4F'
    },
    'ZrB₂超高温陶瓷': {
        'T_max': 2500,
        'epsilon': 0.8,
        'k': 60,
        'rho': 6000,
        'color': '#8B4513'
    },
    '碳化硅 (SiC)': {
        'T_max': 1900,
        'epsilon': 0.9,
        'k': 120,
        'rho': 3200,
        'color': '#708090'
    },
    '钨合金': {
        'T_max': 3000,
        'epsilon': 0.5,
        'k': 170,
        'rho': 19300,
        'color': '#C0C0C0'
    }
}

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 3.1 不同材料的最大允许热流密度
ax1 = axes[0, 0]

material_names = []
q_max_values = []

for name, props in radiative_materials.items():
    T_max = props['T_max']
    epsilon = props['epsilon']
    q_max = epsilon * SIGMA * T_max**4
    material_names.append(name)
    q_max_values.append(q_max / 1e4)  # W/cm²
    
    print(f"\n  {name}:")
    print(f"    最高使用温度: {T_max} K")
    print(f"    发射率: {epsilon}")
    print(f"    最大允许热流密度: {q_max/1e4:.1f} W/cm²")

x_pos = np.arange(len(material_names))
bars = ax1.bar(x_pos, q_max_values, color=[radiative_materials[name]['color'] 
                                          for name in material_names], 
               edgecolor='black', alpha=0.8)
ax1.set_xticks(x_pos)
ax1.set_xticklabels(material_names, rotation=15, ha='right')
ax1.set_ylabel('最大允许热流密度 (W/cm²)', fontsize=11)
ax1.set_title('不同材料的最大允许热流密度', fontsize=12, fontweight='bold')
ax1.grid(alpha=0.3, axis='y')

# 添加数值标签
for bar, q_val in zip(bars, q_max_values):
    height = bar.get_height()
    ax1.annotate(f'{q_val:.0f}',
                xy=(bar.get_x() + bar.get_width() / 2, height),
                xytext=(0, 3), textcoords="offset points",
                ha='center', va='bottom', fontsize=9)

# 3.2 辐射平衡温度随热流密度变化
ax2 = axes[0, 1]
q_range = np.linspace(1e4, 5e5, 100)  # W/m²

for name, props in radiative_materials.items():
    epsilon = props['epsilon']
    T_eq = radiation_equilibrium_temperature(q_range, epsilon)
    ax2.plot(q_range/1e4, T_eq, '-', label=name, linewidth=2, 
            color=props['color'])
    # 标注最高温度限制
    ax2.axhline(y=props['T_max'], color=props['color'], 
               linestyle='--', alpha=0.5)

ax2.set_xlabel('热流密度 (W/cm²)', fontsize=11)
ax2.set_ylabel('辐射平衡温度 (K)', fontsize=11)
ax2.set_title('辐射平衡温度随热流密度变化', fontsize=12, fontweight='bold')
ax2.legend(fontsize=9)
ax2.grid(alpha=0.3)

# 3.3 发射率对辐射冷却的影响
ax3 = axes[1, 0]
epsilon_range = np.linspace(0.3, 0.95, 50)
q_values = [50, 100, 200, 400]  # W/cm²

for q in q_values:
    T_eq = radiation_equilibrium_temperature(q * 1e4, epsilon_range)
    ax3.plot(epsilon_range, T_eq, '-', label=f'q={q} W/cm²', linewidth=2)

ax3.set_xlabel('表面发射率', fontsize=11)
ax3.set_ylabel('辐射平衡温度 (K)', fontsize=11)
ax3.set_title('发射率对辐射平衡温度的影响', fontsize=12, fontweight='bold')
ax3.legend()
ax3.grid(alpha=0.3)

# 3.4 向内部的热传导
ax4 = axes[1, 1]

# 简化的热传导分析
# 假设材料厚度10mm，计算背面温度
thickness = 0.01  # m
q_surface = 200e4  # W/m²

for name, props in radiative_materials.items():
    k = props['k']
    epsilon = props['epsilon']
    
    # 表面辐射平衡温度
    T_surface = radiation_equilibrium_temperature(q_surface, epsilon)
    
    # 简化的背面温度估算 (稳态一维热传导)
    # q = k * (T_surface - T_back) / thickness
    T_back = T_surface - q_surface * thickness / k
    
    # 绘制温度分布示意图
    x = np.linspace(0, thickness*1000, 10)  # mm
    T_dist = T_surface - (T_surface - T_back) * x / (thickness*1000)
    
    ax4.plot(x, T_dist, '-', label=f'{name} (T_back={T_back:.0f}K)', 
            linewidth=2, color=props['color'])

ax4.set_xlabel('深度 (mm)', fontsize=11)
ax4.set_ylabel('温度 (K)', fontsize=11)
ax4.set_title('辐射冷却材料内温度分布 (q=200 W/cm²)', fontsize=12, fontweight='bold')
ax4.legend(fontsize=8)
ax4.grid(alpha=0.3)

plt.tight_layout()
plt.savefig('case3_radiative_cooling.png', dpi=150, bbox_inches='tight')
plt.close()

print("\n✓ 案例3结果已保存: case3_radiative_cooling.png")


# ==============================================================================
# 案例4：热防护系统综合设计
# ==============================================================================
print("\n" + "=" * 80)
print("案例4：热防护系统综合设计")
print("=" * 80)

# 设计案例：载人飞船返回舱
print("\n【设计案例：载人飞船返回舱热防护系统】")

# 飞行器参数
spacecraft = {
    'name': '载人飞船返回舱',
    'mass': 5000,  # kg
    'diameter': 3.0,  # m
    'bottom_area': np.pi * (3.0/2)**2,  # m²
    'side_area': np.pi * 3.0 * 2.0,  # m² (假设高度2m)
    'back_area': np.pi * (3.0/2)**2,  # m²
}

# 再入参数
reentry_params = {
    'V0': 7500,  # m/s
    'h0': 120e3,  # m
    'gamma': np.radians(-1.5),  # rad
}

# 计算热环境
t, V, h, gamma, rho = reentry_trajectory(
    reentry_params['V0'], reentry_params['gamma'], 
    reentry_params['h0'], spacecraft['mass'], CD, A, dt=2.0
)
q = stagnation_heat_flux(V, rho, Rn=1.0)

# 热环境统计
q_max = np.max(q)
q_avg = np.mean(q)
Q_total = np.trapezoid(q, t)

print(f"\n  热环境分析:")
print(f"    最大热流密度: {q_max/1e4:.1f} W/cm²")
print(f"    平均热流密度: {q_avg/1e4:.1f} W/cm²")
print(f"    总加热量: {Q_total/1e6:.1f} MJ/m²")
print(f"    再入时间: {t[-1]:.1f} s")

# 热防护方案设计
print(f"\n  热防护方案设计:")

# 底部 (最高热流)
q_bottom = q_max  # 驻点热流
design_bottom = {
    'location': '底部',
    'q_max': q_bottom,
    'material': 'AVCOAT',
    'thickness': 0.060,  # m
    'area': spacecraft['bottom_area'],
    'mass': 0  # 待计算
}
design_bottom['mass'] = design_bottom['thickness'] * design_bottom['area'] * ablative_material['rho']

# 侧面 (中等热流)
q_side = q_max * 0.6  # 约为驻点的60%
design_side = {
    'location': '侧面',
    'q_max': q_side,
    'material': 'AVCOAT',
    'thickness': 0.040,  # m
    'area': spacecraft['side_area'],
    'mass': 0
}
design_side['mass'] = design_side['thickness'] * design_side['area'] * ablative_material['rho']

# 背部 (低热流)
q_back = q_max * 0.3  # 约为驻点的30%
design_back = {
    'location': '背部',
    'q_max': q_back,
    'material': '低密度烧蚀材料',
    'thickness': 0.025,  # m
    'area': spacecraft['back_area'],
    'mass': 0
}
design_back['mass'] = design_back['thickness'] * design_back['area'] * 300  # 低密度材料密度约300 kg/m³

total_mass = design_bottom['mass'] + design_side['mass'] + design_back['mass']

print(f"\n    底部: {design_bottom['material']}, 厚度{design_bottom['thickness']*1000:.0f}mm, "
      f"面积{design_bottom['area']:.1f}m², 质量{design_bottom['mass']:.1f}kg")
print(f"    侧面: {design_side['material']}, 厚度{design_side['thickness']*1000:.0f}mm, "
      f"面积{design_side['area']:.1f}m², 质量{design_side['mass']:.1f}kg")
print(f"    背部: {design_back['material']}, 厚度{design_back['thickness']*1000:.0f}mm, "
      f"面积{design_back['area']:.1f}m², 质量{design_back['mass']:.1f}kg")
print(f"\n    热防护系统总质量: {total_mass:.1f} kg")
print(f"    占飞行器质量比例: {total_mass/spacecraft['mass']*100:.1f}%")

# 创建可视化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 4.1 热防护系统质量分布
ax1 = axes[0, 0]
locations = [design_bottom['location'], design_side['location'], design_back['location']]
masses = [design_bottom['mass'], design_side['mass'], design_back['mass']]
colors = ['#FF6B6B', '#4ECDC4', '#45B7D1']

wedges, texts, autotexts = ax1.pie(masses, labels=locations, autopct='%1.1f%%',
                                     colors=colors, startangle=90,
                                     textprops={'fontsize': 11})
ax1.set_title(f'热防护系统质量分布 (总质量{total_mass:.0f}kg)', fontsize=12, fontweight='bold')

# 4.2 不同区域热流密度和厚度
ax2 = axes[0, 1]
x_pos = np.arange(len(locations))
width = 0.35

# 热流密度 (左轴)
ax2_twin = ax2.twinx()
bars1 = ax2.bar(x_pos - width/2, [design_bottom['q_max']/1e4, design_side['q_max']/1e4, 
                                 design_back['q_max']/1e4], 
               width, label='最大热流密度', color='#FF6B6B', alpha=0.8)
ax2.set_ylabel('最大热流密度 (W/cm²)', fontsize=11, color='#FF6B6B')
ax2.tick_params(axis='y', labelcolor='#FF6B6B')

# 厚度 (右轴)
bars2 = ax2_twin.bar(x_pos + width/2, [design_bottom['thickness']*1000, 
                                       design_side['thickness']*1000, 
                                       design_back['thickness']*1000], 
                    width, label='热防护层厚度', color='#4ECDC4', alpha=0.8)
ax2_twin.set_ylabel('热防护层厚度 (mm)', fontsize=11, color='#4ECDC4')
ax2_twin.tick_params(axis='y', labelcolor='#4ECDC4')

ax2.set_xticks(x_pos)
ax2.set_xticklabels(locations)
ax2.set_title('不同区域热流密度和热防护层厚度', fontsize=12, fontweight='bold')
ax2.grid(alpha=0.3, axis='y')

# 添加图例
lines1, labels1 = ax2.get_legend_handles_labels()
lines2, labels2 = ax2_twin.get_legend_handles_labels()
ax2.legend(lines1 + lines2, labels1 + labels2, loc='upper right')

# 4.3 设计验证：背面温度估算
ax3 = axes[1, 0]

# 简化的背面温度估算
# 假设稳态热传导: q = k * (T_surface - T_back) / thickness
designs = [design_bottom, design_side, design_back]
k_material = ablative_material['k']

T_back_estimates = []
for design in designs:
    q_loc = design['q_max']
    thickness = design['thickness']
    # 简化的背面温度估算 (假设表面温度3000K)
    T_surface = 3000  # K
    T_back = T_surface - q_loc * thickness / k_material
    T_back_estimates.append(T_back)
    print(f"\n  {design['location']}: 估算背面温度 {T_back:.0f}K ({T_back-273:.0f}°C)")

bars = ax3.bar(locations, T_back_estimates, color=colors, edgecolor='black', alpha=0.8)
ax3.axhline(y=350+273, color='red', linestyle='--', linewidth=2, label='允许温度上限 (350°C)')
ax3.set_ylabel('估算背面温度 (K)', fontsize=11)
ax3.set_title('热防护设计验证：背面温度估算', fontsize=12, fontweight='bold')
ax3.legend()
ax3.grid(alpha=0.3, axis='y')

# 添加数值标签
for bar, T_val in zip(bars, T_back_estimates):
    height = bar.get_height()
    ax3.annotate(f'{T_val-273:.0f}°C',
                xy=(bar.get_x() + bar.get_width() / 2, height),
                xytext=(0, 3), textcoords="offset points",
                ha='center', va='bottom', fontsize=10)

# 4.4 不同任务类型的热防护方案对比
ax4 = axes[1, 1]

mission_types = ['载人飞船', '航天飞机', '再入弹头', '高超声速飞行器']
max_heat_flux = [400, 150, 800, 600]  # W/cm²
tps_types = ['烧蚀式', '辐射式', '烧蚀式', '主动冷却']
tps_masses = [800, 5000, 200, 1000]  # kg (估算值)

x_pos = np.arange(len(mission_types))
width = 0.6

bars = ax4.bar(x_pos, tps_masses, width, color=['#FF6B6B', '#4ECDC4', '#45B7D1', '#96CEB4'],
               edgecolor='black', alpha=0.8)

# 在柱状图上标注热流密度和热防护类型
for i, (bar, q, tps_type) in enumerate(zip(bars, max_heat_flux, tps_types)):
    height = bar.get_height()
    ax4.annotate(f'{q} W/cm²\n{tps_type}',
                xy=(bar.get_x() + bar.get_width() / 2, height),
                xytext=(0, 3), textcoords="offset points",
                ha='center', va='bottom', fontsize=9)

ax4.set_xticks(x_pos)
ax4.set_xticklabels(mission_types, rotation=15, ha='right')
ax4.set_ylabel('热防护系统质量 (kg)', fontsize=11)
ax4.set_title('不同任务类型的热防护方案对比', fontsize=12, fontweight='bold')
ax4.grid(alpha=0.3, axis='y')

plt.tight_layout()
plt.savefig('case4_tps_design.png', dpi=150, bbox_inches='tight')
plt.close()

print("\n✓ 案例4结果已保存: case4_tps_design.png")


# ==============================================================================
# 生成GIF动画：再入过程热环境变化
# ==============================================================================
print("\n" + "=" * 80)
print("生成GIF动画：再入过程热环境变化")
print("=" * 80)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# 使用γ=-1.5°的数据
data = results[-1.5]
t_anim = data['t']
h_anim = data['h']
V_anim = data['V']
q_anim = data['q']

# 降采样以加快动画生成
skip = 10
t_anim = t_anim[::skip]
h_anim = h_anim[::skip]
V_anim = V_anim[::skip]
q_anim = q_anim[::skip]

def init():
    ax1.clear()
    ax2.clear()
    return []

def update(frame):
    ax1.clear()
    ax2.clear()
    
    # 左图：高度-速度曲线
    ax1.plot(V_anim[:frame+1]/1e3, h_anim[:frame+1]/1e3, 'b-', linewidth=2)
    ax1.scatter([V_anim[frame]/1e3], [h_anim[frame]/1e3], c='red', s=100, zorder=5)
    ax1.set_xlim(0, 8)
    ax1.set_ylim(0, 120)
    ax1.set_xlabel('速度 (km/s)', fontsize=11)
    ax1.set_ylabel('高度 (km)', fontsize=11)
    ax1.set_title('再入轨迹', fontsize=12, fontweight='bold')
    ax1.grid(alpha=0.3)
    ax1.invert_yaxis()
    
    # 添加当前状态标注
    ax1.annotate(f't={t_anim[frame]:.0f}s\nV={V_anim[frame]/1e3:.1f}km/s\nh={h_anim[frame]/1e3:.0f}km',
                xy=(V_anim[frame]/1e3, h_anim[frame]/1e3),
                xytext=(V_anim[frame]/1e3+0.5, h_anim[frame]/1e3+10),
                fontsize=10, arrowprops=dict(arrowstyle='->', color='red'))
    
    # 右图：热流密度-时间曲线
    ax2.fill_between(t_anim[:frame+1], 0, q_anim[:frame+1]/1e4, alpha=0.3, color='red')
    ax2.plot(t_anim[:frame+1], q_anim[:frame+1]/1e4, 'r-', linewidth=2)
    ax2.scatter([t_anim[frame]], [q_anim[frame]/1e4], c='blue', s=100, zorder=5)
    ax2.set_xlim(0, t_anim[-1])
    ax2.set_ylim(0, np.max(q_anim)/1e4 * 1.1)
    ax2.set_xlabel('时间 (s)', fontsize=11)
    ax2.set_ylabel('热流密度 (W/cm²)', fontsize=11)
    ax2.set_title('驻点热流密度', fontsize=12, fontweight='bold')
    ax2.grid(alpha=0.3)
    
    # 添加当前热流标注
    ax2.annotate(f'q={q_anim[frame]/1e4:.1f} W/cm²',
                xy=(t_anim[frame], q_anim[frame]/1e4),
                xytext=(t_anim[frame]+20, q_anim[frame]/1e4+20),
                fontsize=10, arrowprops=dict(arrowstyle='->', color='blue'))
    
    return []

# 创建动画
print("  正在生成GIF动画...")
anim = FuncAnimation(fig, update, frames=len(t_anim), init_func=init, blit=False)
writer = PillowWriter(fps=15)
anim.save('reentry_heating_animation.gif', writer=writer)
plt.close()

print("\n✓ GIF动画已保存: reentry_heating_animation.gif")


# ==============================================================================
# 总结输出
# ==============================================================================
print("\n" + "=" * 80)
print("所有案例运行完成！")
print("=" * 80)

print("\n生成的文件：")
print("  - case1_reentry_heating.png")
print("  - case2_ablative_tps.png")
print("  - case3_radiative_cooling.png")
print("  - case4_tps_design.png")
print("  - reentry_heating_animation.gif")

print("\n" + "=" * 80)
print("主题052：再入飞行器热防护辐射 - 仿真完成")
print("=" * 80)