辐射换热仿真-主题051_航天器热控系统辐射分析-主题051_航天器热控系统辐射分析

kkchenjj

52人浏览 · 2026-04-27 21:33:51

kkchenjj · 2026-04-27 21:33:51 发布

主题051：航天器热控系统辐射分析

摘要

航天器在太空极端环境中运行，热控系统是保证航天器正常工作的关键。太空环境具有高真空、强辐射、极端温度变化等特点，航天器与外界的换热几乎完全依靠辐射。本主题从辐射换热的基本原理出发，系统讲解航天器热控系统的辐射换热机制，包括太空热环境特性、航天器表面辐射特性、热控涂层和多层隔热材料的设计原理、辐射散热器的优化设计等。通过建立航天器热平衡模型，分析不同轨道、不同姿态下航天器的热环境，设计热控系统方案。结合Python仿真程序，模拟航天器在轨运行期间的热状态变化，计算温度分布和热流密度，为航天器热控系统设计提供科学依据。

关键词： 航天器、热控系统、辐射散热、太空环境、热平衡、热控涂层、多层隔热

在这里插入图片描述

1. 太空热环境概述

1.1 太空环境的基本特征

太空环境与地球表面环境截然不同，具有以下显著特征：

高真空环境：

近地轨道（LEO，200-2000 km）大气压力约 $10^{-6}$ - $10^{-9}$ Pa
高真空导致对流换热几乎为零
航天器与外界的换热主要依靠辐射

极端温度范围：

太阳直射区域：可达120°C以上
地球阴影区域：可低至-180°C以下
温差可达300°C以上

强辐射环境：

太阳辐射：1361 W/m²（太阳常数）
地球反照：约340 W/m²（平均）
地球红外辐射：约240 W/m²（平均）

微重力环境：

重力加速度约 $10^{-6}$ g
自然对流消失
相变换热（如热管）依靠毛细力

1.2 空间热流密度来源

航天器在轨运行时，接收的热流密度主要来自以下几个方面：

（1）太阳直接辐射（ $q_s$ ）：

太阳辐射是航天器最主要的外部热源。太阳辐射可视为温度为5778 K的黑体辐射，其光谱分布在紫外、可见光和近红外波段。

太阳直接辐射热流密度：

$qs=S0⋅f1⋅f2⋅cos⁡θq_s = S_0 \cdot f_1 \cdot f_2 \cdot \cos\theta$

其中：

$S_0 = 1361$ W/m²：太阳常数
$f_1$ ：日地距离修正因子（0.97-1.03）
$f_2$ ：太阳活动修正因子（0.9-1.1）
$θ\theta$ ：太阳入射角（表面法线与太阳光线的夹角）

（2）地球反照（ $q_a$ ）：

地球反射的太阳辐射称为地球反照。地球平均反照率约0.30，但随地区、季节和云层变化。

地球反照热流密度：

$qa=αE⋅S0⋅F12⋅cos⁡θsq_a = \alpha_E \cdot S_0 \cdot F_{12} \cdot \cos\theta_s$

其中：

$αE\alpha_E$ ：地球反照率（约0.30）
$F_{12}$ ：地球对航天器的角系数
$θs\theta_s$ ：太阳天顶角

对于近地轨道航天器，地球反照角系数可近似为：

$F12=sin⁡2ρF_{12} = \sin^2\rho$

其中 $ρ\rho$ 是地球对航天器的半张角：

$sin⁡ρ=RERE+h\sin\rho = \frac{R_E}{R_E + h}$

$R_E$ 是地球半径（6371 km）， $h$ 是轨道高度。

（3）地球红外辐射（ $q_{IR}$ ）：

地球吸收太阳辐射后向外发射红外辐射，这是航天器在阴影区的主要热源。

地球红外辐射热流密度：

$qIR=σTE4⋅F12q_{IR} = \sigma T_E^4 \cdot F_{12}$

其中：

$σ=5.67×10−8\sigma = 5.67 \times 10^{-8}$ W/m²·K⁴：斯蒂芬-玻尔兹曼常数
$T_E$ ：地球有效辐射温度（约255 K）
$F_{12}$ ：地球对航天器的角系数

（4）航天器内部热源（ $q_{int}$ ）：

航天器内部设备产生的热量，包括：

电子设备功耗
电池放电热量
人员代谢热（载人航天器）
推进系统热量

1.3 不同轨道的空间热环境

低地球轨道（LEO，200-2000 km）：

LEO轨道是最常用的航天器轨道，其热环境特点：

轨道周期约90-120分钟
光照期约60分钟，阴影期约30分钟
地球反照和红外辐射显著
大气阻力导致轨道衰减

热流密度变化规律：

光照期：太阳辐射 + 地球反照 + 地球红外
阴影期：仅地球红外辐射

地球同步轨道（GEO，35786 km）：

GEO轨道的热环境特点：

轨道周期24小时，与地球自转同步
地球半张角仅8.7°
地球反照和红外辐射较弱
几乎持续受太阳照射（除春分、秋分前后短暂地影期）

深空轨道：

深空探测任务的热环境特点：

远离地球，地球热源可忽略
仅太阳辐射为主要热源
温度极端（向阳面高温，背阳面极低温）
需要考虑行星/月球表面反照和红外辐射

月球轨道：

月球轨道热环境特点：

月球无大气，无反照
月球表面红外辐射显著（白天约400 K，夜晚约100 K）
月食期间经历极端温度变化

2. 航天器热平衡分析

2.1 航天器热平衡方程

航天器在轨运行时的热平衡方程为：

$Q_{in} + Q_{gen} = Q_{out} + Q_{stored}$

其中：

$Q_{in}$ ：外部输入热量（太阳辐射、地球反照、地球红外）
$Q_{gen}$ ：内部产生热量（设备功耗）
$Q_{out}$ ：向外散失热量（辐射散热）
$Q_{stored}$ ：储存热量（温度变化）

在稳态条件下（ $Q_{stored} = 0$ ）：

$Q_{in} + Q_{gen} = Q_{out}$

对于航天器表面微元面积 $d A$ ，热平衡方程可写为：

$αsqs+αaqa+εqIR+qint=εσT4\alpha_s q_s + \alpha_a q_a + \varepsilon q_{IR} + q_{int} = \varepsilon \sigma T^4$

其中：

$αs\alpha_s$ ：太阳吸收率
$αa\alpha_a$ ：反照吸收率（通常 $αa≈αs\alpha_a \approx \alpha_s$ ）
$ε\varepsilon$ ：红外发射率
$T$ ：表面温度

2.2 表面辐射特性参数

航天器表面的辐射特性主要由两个参数描述：

太阳吸收率（ $αs\alpha_s$ ）：

太阳吸收率表示表面对太阳辐射（0.3-2.5 μm）的吸收能力：

$αs=∫0.32.5αλEsλdλ∫0.32.5Esλdλ\alpha_s = \frac{\int_{0.3}^{2.5} \alpha_\lambda E_{s\lambda} d\lambda}{\int_{0.3}^{2.5} E_{s\lambda} d\lambda}$

其中 $EsλE_{s\lambda}$ 是太阳辐射光谱。

典型材料的太阳吸收率：

抛光铝：0.10-0.15
阳极氧化铝：0.25-0.40
白色涂料：0.20-0.30
黑色涂料：0.90-0.95
镀金表面：0.20-0.25

红外发射率（ $ε\varepsilon$ ）：

红外发射率表示表面在红外波段（通常2.5-50 μm）的发射能力：

$ε=∫2.550ελEbλ(T)dλ∫2.550Ebλ(T)dλ\varepsilon = \frac{\int_{2.5}^{50} \varepsilon_\lambda E_{b\lambda}(T) d\lambda}{\int_{2.5}^{50} E_{b\lambda}(T) d\lambda}$

典型材料的红外发射率：

抛光金属：0.02-0.10
阳极氧化铝：0.80-0.90
白色涂料：0.85-0.90
黑色涂料：0.90-0.95

吸收-发射比（ $αs/ε\alpha_s/\varepsilon$ ）：

吸收-发射比是热控设计的关键参数：

比值高（>1）：表面温度升高
比值低（<1）：表面温度降低
比值=1：温度由外部热流决定

典型热控涂层的 $αs/ε\alpha_s/\varepsilon$ ：

光学太阳反射器（OSR）：0.05-0.15
白色涂料：0.20-0.35
黑色涂料：0.90-1.00
镀金表面：2.0-4.0

2.3 航天器温度估算

在稳态条件下，航天器表面温度可由热平衡方程估算：

$\left[\frac{\alpha_s q_s + \alpha_a q_a + \varepsilon q_{IR} + q_{int}}{\varepsilon \sigma}\right]^{1/4}$

简化估算示例：

假设某航天器表面：

太阳吸收率 $αs=0.3\alpha_s = 0.3$
红外发射率 $ε=0.8\varepsilon = 0.8$
接收太阳辐射 $q_s = 1361$ W/m²
地球反照 $q_a = 300$ W/m²
地球红外 $q_{IR} = 200$ W/m²
内部热流 $q_{int} = 50$ W/m²

则表面温度：

$\left[\frac{0.3 \times 1361 + 0.3 \times 300 + 0.8 \times 200 + 50}{0.8 \times 5.67 \times 10^{-8}}\right]^{1/4}$

$\left[\frac{408.3 + 90 + 160 + 50}{4.536 \times 10^{-8}}\right]^{1/4}$

$\left[1.56 \times 10^{10}\right]^{1/4} \approx 353 \text{ K} \approx 80°C$

这个温度对于大多数航天器设备来说过高，需要通过热控措施降低。

3. 热控涂层与表面处理

3.1 热控涂层的分类与原理

热控涂层是航天器热控系统的基础，通过调节表面的太阳吸收率和红外发射率来控制温度。

（1）光学太阳反射器（OSR）：

OSR是目前最高效的热控涂层之一，通常由石英玻璃或熔融二氧化硅制成，背面镀银或铝。

特性：

太阳吸收率：0.05-0.15
红外发射率：0.80-0.85
吸收-发射比：0.06-0.18

工作原理：

石英玻璃对太阳辐射（短波）几乎透明
金属镀层反射大部分太阳辐射
石英玻璃在红外波段强吸收并发射

（2）白色涂料：

白色涂料通过散射反射太阳辐射，同时在红外波段高发射。

特性：

太阳吸收率：0.20-0.30
红外发射率：0.85-0.90
吸收-发射比：0.22-0.35

常用类型：

氧化锌涂料
二氧化钛涂料
硅酸钾涂料

（3）黑色涂料：

黑色涂料用于需要吸收太阳辐射或高发射率的场合。

特性：

太阳吸收率：0.90-0.95
红外发射率：0.90-0.95
吸收-发射比：0.95-1.05

（4）金属表面：

抛光金属表面具有低太阳吸收率和低红外发射率。

特性：

太阳吸收率：0.10-0.30
红外发射率：0.02-0.10
吸收-发射比：1.0-10.0

应用：

镀金表面：用于需要保温的部件
镀铝表面：用于辐射屏蔽
不锈钢表面：用于结构件

3.2 热控涂层的性能退化

热控涂层在太空环境中会发生性能退化，主要因素包括：

紫外辐射损伤：

紫外光子破坏有机粘结剂
导致涂层变黄、吸收率增加
在LEO轨道尤为严重

带电粒子辐射：

高能质子、电子损伤涂层结构
导致光学性能变化
在范艾伦辐射带影响显著

原子氧剥蚀：

LEO轨道原子氧与涂层反应
导致涂层质量损失、厚度减小
对有机涂料影响严重

微流星体撞击：

高速微粒撞击表面
造成局部损伤、表面粗糙化
增加太阳吸收率

性能退化模型：

太阳吸收率随时间变化的经验公式：

$αs(t)=αs0+Δαs⋅(1−e−t/τ)\alpha_s(t) = \alpha_{s0} + \Delta\alpha_s \cdot (1 - e^{-t/\tau})$

其中：

$αs0\alpha_{s0}$ ：初始太阳吸收率
$Δαs\Delta\alpha_s$ ：最大增量（通常0.05-0.15）
$τ\tau$ ：时间常数（通常1-5年）
$t$ ：在轨时间

3.3 多层隔热材料（MLI）

多层隔热材料是航天器热控的重要手段，用于减少辐射换热。

结构：

MLI由多层反射屏和间隔层组成：

反射屏：镀铝聚酰亚胺薄膜或铝箔（厚度约6-25 μm）
间隔层：涤纶网或玻璃纤维网（防止反射屏接触）
层数：通常10-30层

隔热原理：

多层隔热通过以下机制减少热流：

每层反射屏反射大部分辐射
多层反射形成辐射屏蔽
间隔层减少传导换热
高真空环境消除对流

等效热导率：

MLI的等效热导率极低：

$W/m⋅K\lambda_{eff} = 10^{-4} - 10^{-3} \text{ W/m·K}$

相比之下：

静止空气：0.026 W/m·K
聚氨酯泡沫：0.02-0.03 W/m·K
不锈钢：15 W/m·K

热流计算：

通过MLI的热流密度：

$\frac{\sigma(T_h^4 - T_c^4)}{\frac{1}{\varepsilon_h} + \frac{1}{\varepsilon_c} - 1 + N(\frac{2}{\varepsilon} - 1)}$

其中：

$T_h, T_c$ ：热端和冷端温度
$εh,εc\varepsilon_h, \varepsilon_c$ ：内外表面发射率
$N$ ：反射屏层数
$ε\varepsilon$ ：反射屏发射率

简化公式（当 $N$ 较大时）：

$\approx \frac{\sigma(T_h^4 - T_c^4)}{(N+1)(\frac{2}{\varepsilon} - 1)}$

4. 辐射散热器设计

4.1 辐射散热器的工作原理

辐射散热器是航天器排散废热的主要手段，通过增大辐射面积、提高发射率来增强散热能力。

散热功率计算：

辐射散热器的散热功率：

$Qrad=εσAT4Q_{rad} = \varepsilon \sigma A T^4$

其中：

$ε\varepsilon$ ：散热器表面发射率
$A$ ：散热器面积
$T$ ：散热器温度

散热器面积估算：

假设需要排散热量 $Q$ ，散热器温度 $T$ ，则所需面积：

$\frac{Q}{\varepsilon \sigma T^4}$

示例：

需要排散1000 W热量，散热器温度300 K，发射率0.85：

$\frac{1000}{0.85 \times 5.67 \times 10^{-8} \times 300^4}$

$\frac{1000}{0.85 \times 5.67 \times 10^{-8} \times 8.1 \times 10^9}$

$\frac{1000}{390.6} \approx 2.56 \text{ m}^2$

4.2 散热器类型

（1）体装式散热器：

直接利用航天器结构表面作为散热器。

特点：

结构简单、质量小
散热面积受限
受外部热流影响大

应用：低功耗设备、小型航天器

（2）展开式散热器：

通过展开机构增大散热面积。

特点：

发射时收拢，入轨后展开
散热面积大
增加系统复杂度和质量

类型：

刚性散热器面板
柔性散热器（可卷绕）
可展开散热翼

（3）热管散热器：

利用热管将热量从热源传输到散热器。

特点：

高效传热
等温性好
可远距离传输热量

工作原理：

蒸发段吸收热量，工质蒸发
蒸汽流向冷凝段
冷凝段释放热量，工质冷凝
毛细芯将液体回流到蒸发段

（4）环路热管（LHP）和毛细泵环（CPL）：

先进的相变传热系统，用于大功率散热。

特点：

传输距离远（可达数米）
传热功率大（可达数千瓦）
可实现温度控制

4.3 散热器优化设计

散热器朝向优化：

散热器应尽量避免接收外部热流：

背向太阳
背向地球
与太阳成较大角度

表面发射率优化：

提高发射率可增强散热：

使用高发射率涂层（黑色涂料、阳极氧化铝）
发射率可达0.90-0.95

散热器温度优化：

提高散热器温度可显著增加散热功率：

$\propto T^4$

温度从300 K提高到350 K，散热功率增加：

$350/300)^4 = 1.88$

散热功率增加88%。

散热器面积优化：

在质量和体积约束下优化散热器面积：

采用蜂窝结构增加比表面积
使用展开机构增大面积
优化散热器形状（ fins、肋片）

5. 航天器热控系统设计

5.1 热控系统的组成

航天器热控系统通常包括以下部分：

被动热控措施：

热控涂层（OSR、涂料、金属表面）
多层隔热材料（MLI）
热管
导热填料
辐射散热器

主动热控措施：

电加热器
百叶窗（可变发射率装置）
热开关
机械泵驱动回路
蒸气压缩制冷

热控系统的选择原则：

优先采用被动措施（简单、可靠）
主动措施用于需要精确控温的场合
考虑质量、功耗、可靠性权衡

5.2 热设计流程

航天器热设计的一般流程：

（1）任务分析：

确定轨道参数
分析热环境
识别极端工况

（2）热分析建模：

建立几何模型
划分热节点
计算角系数和辐射耦合

（3）热平衡计算：

计算各节点温度
分析稳态和瞬态工况
确定极端温度

（4）热控方案设计：

选择热控措施
确定涂层和隔热方案
设计散热器

（5）热试验验证：

热真空试验
热平衡试验
验证分析模型

5.3 典型航天器热控方案

通信卫星：

通信卫星热控特点：

大功率电子设备
需要精确温度控制
南北向散热板

热控方案：

南北散热板：OSR + 热管
东西板：MLI包覆
对地板：MLI包覆
对天板：OSR
内部：热管网络 + 辐射器

载人飞船：

载人飞船热控特点：

人员热舒适要求
生命保障系统
返回大气层再入热

热控方案：

返回舱：防热结构 + 烧蚀材料
轨道舱：MLI + 流体回路
主动控温：空调系统

深空探测器：

深空探测器热控特点：

极端温度环境
放射性同位素热源（RHU/RTG）
长期运行

热控方案：

多层隔热（30-50层）
放射性同位素加热
百叶窗控温
高发射率散热器

6. Python仿真案例分析

6.1 案例1：航天器表面温度计算

本案例计算不同轨道、不同表面特性下航天器的温度。

模拟内容：

计算LEO轨道航天器在光照期和阴影期的温度
分析不同热控涂层（OSR、白漆、黑漆）的温度差异
计算地球反照和红外辐射对温度的影响
对比不同轨道高度（300 km、500 km、800 km）的热环境

关键公式：

热平衡方程：

$\left[\frac{\alpha_s q_s + \alpha_a q_a + \varepsilon q_{IR} + q_{int}}{\varepsilon \sigma}\right]^{1/4}$

地球反照角系数：

$F12=sin⁡2ρ=(RERE+h)2F_{12} = \sin^2\rho = \left(\frac{R_E}{R_E + h}\right)^2$

模拟结果分析：

LEO 400 km轨道：

光照期：太阳辐射1361 W/m²，地球反照约300 W/m²，地球红外约200 W/m²
阴影期：仅地球红外约200 W/m²

不同涂层温度对比：

涂层	$αs\alpha_s$	$ε\varepsilon$	光照期温度	阴影期温度
OSR	0.10	0.85	-20°C	-30°C
白漆	0.25	0.88	15°C	-25°C
黑漆	0.90	0.90	85°C	-20°C
镀金	0.25	0.05	180°C	50°C

6.2 案例2：多层隔热材料性能分析

本案例模拟多层隔热材料的热性能。

模拟内容：

计算不同层数MLI的热流密度
分析反射屏发射率对隔热性能的影响
模拟MLI在极端温差下的性能
计算等效热导率

模拟方法：

多层辐射换热网络法：

$\frac{\sigma(T_h^4 - T_c^4)}{R_{total}}$

其中总热阻：

$Rtotal=1εh+1εc−1+∑i=1N−1(2εi−1)R_{total} = \frac{1}{\varepsilon_h} + \frac{1}{\varepsilon_c} - 1 + \sum_{i=1}^{N-1}\left(\frac{2}{\varepsilon_i} - 1\right)$

模拟结果分析：

层数对隔热性能的影响：

层数	热流密度 (W/m²)	等效热导率 (W/m·K)
5	25.3	$5.0 \times 10^{-3}$
10	13.8	$2.7 \times 10^{-3}$
20	7.2	$1.4 \times 10^{-3}$
30	4.9	$9.6 \times 10^{-4}$

（假设： $T_h = 300$ K， $T_c = 100$ K， $ε=0.05\varepsilon = 0.05$ ）

结论：

层数增加可显著降低热流
但层数超过20后，边际效益递减
实际工程通常采用15-25层

6.3 案例3：辐射散热器设计优化

本案例设计并优化辐射散热器。

模拟内容：

计算不同面积、温度、发射率下的散热功率
优化散热器朝向以最小化外部热流吸收
设计热管散热器系统
分析散热器在轨温度变化

优化目标：

在给定质量约束下最大化散热功率：

$max⁡Qrad=εσAT4\max Q_{rad} = \varepsilon \sigma A T^4$

约束条件：

质量约束： $\leq m_{max}$
体积约束： $\leq V_{max}$
温度约束： $\leq T_{max}$

模拟结果分析：

散热器面积与散热功率：

面积 (m²)	温度 (K)	发射率	散热功率 (W)
1.0	300	0.85	390
2.0	300	0.85	780
2.0	320	0.85	1010
2.0	320	0.95	1130

优化结论：

提高温度对散热功率影响最大（ $T^4$ 关系）
提高发射率也有显著效果
面积增加线性提高散热功率

6.4 案例4：航天器在轨热状态仿真

本案例综合仿真航天器在轨运行期间的热状态。

模拟内容：

建立简化的航天器几何模型
计算一个轨道周期内的温度变化
分析不同姿态对温度的影响
评估热控系统的有效性

航天器模型：

简化立方体航天器：

尺寸：1m × 1m × 1m
轨道：LEO 400 km，圆轨道
周期：93分钟
光照期：60分钟
阴影期：33分钟

表面特性：

向阳面：OSR（ $αs=0.10,ε=0.85\alpha_s=0.10, \varepsilon=0.85$ ）
背阳面：MLI包覆
侧面：白漆（ $αs=0.25,ε=0.88\alpha_s=0.25, \varepsilon=0.88$ ）

内部热源：500 W

瞬态热平衡方程：

$c_p \frac{dT}{dt} = Q_{in}(t) + Q_{gen} - Q_{out}(t)$

模拟结果分析：

温度变化曲线：

光照期开始：温度约-20°C
光照期结束：温度约+40°C
阴影期结束：温度约-10°C
一个周期内温度波动约60°C

热控措施效果：

无热控：温度范围-100°C至+150°C
有热控：温度范围-20°C至+40°C
热控系统有效限制了温度波动

7. 结果讨论与工程应用

7.1 关键发现总结

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

（1）热控涂层的选择至关重要

不同涂层的温度差异可达100°C以上。OSR是最有效的散热涂层，可使表面温度降低50-80°C。镀金表面虽然太阳吸收率低，但红外发射率也低，不适合散热。

（2）多层隔热效果显著

20层MLI可将热流密度降低至无隔热时的5%以下。但层数超过20后，边际效益递减，需要权衡质量和隔热性能。

（3）散热器设计需要综合考虑

散热器面积、温度和发射率都影响散热功率。提高温度效果最显著，但受设备耐温限制。实际设计需要在质量、体积和散热能力之间权衡。

（4）轨道热环境变化剧烈

LEO轨道航天器在一个周期内温度波动可达60-100°C。热控系统需要适应这种周期性变化，保持设备温度在允许范围内。

7.2 热控系统设计要点

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

（1）热环境分析是基础

准确分析不同轨道、不同姿态下的热环境是热控设计的前提。需要考虑极端工况（最大/最小热流）和典型工况。

（2）被动措施优先

在满足要求的前提下，优先采用被动热控措施（涂层、MLI、热管），因其简单可靠、不消耗功率。

（3）主动措施精确控温

对于需要精确控温的设备（如光学仪器、电子设备），采用主动热控措施（电加热器、百叶窗）。

（4）热试验验证必不可少

热分析模型存在不确定性，必须通过热真空试验验证。试验数据也用于修正分析模型。

7.3 新技术发展趋势

航天器热控技术正在向以下方向发展：

智能热控：

可变发射率涂层（电致变色、热致变色）
自适应热控系统
基于微流控的主动热控

高效传热：

纳米流体
脉动热管
环路热管技术改进

轻量化设计：

新型隔热材料（气凝胶、真空绝热板）
3D打印散热器
可展开柔性散热器

深空探测热控：

放射性同位素电源（RPS）
斯特林制冷机
相变储能材料

附录：Python仿真程序说明

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

A.1 程序结构

run_simulation.py
├── 案例1：航天器表面温度计算
│   ├── LEO轨道热环境计算
│   ├── 不同涂层温度对比
│   └── 轨道高度影响分析
├── 案例2：多层隔热材料性能分析
│   ├── 不同层数热流计算
│   ├── 等效热导率分析
│   └── 极端温差性能
├── 案例3：辐射散热器设计优化
│   ├── 散热功率计算
│   ├── 朝向优化
│   └── 参数敏感性分析
└── 案例4：航天器在轨热状态仿真
    ├── 轨道热环境建模
    ├── 瞬态温度计算
    └── 热控效果评估

A.2 关键算法

地球角系数计算：

def earth_view_factor(h, R_E=6371e3):
    """计算地球对航天器的角系数"""
    sin_rho = R_E / (R_E + h)
    return sin_rho**2

热平衡温度计算：

def surface_temperature(alpha_s, epsilon, q_s, q_a, q_IR, q_int=0):
    """计算表面热平衡温度"""
    sigma = 5.67e-8
    q_in = alpha_s * q_s + alpha_s * q_a + epsilon * q_IR + q_int
    T = (q_in / (epsilon * sigma))**0.25
    return T

多层隔热热流计算：

def mli_heat_flux(T_hot, T_cold, N, epsilon=0.05):
    """计算通过MLI的热流密度"""
    sigma = 5.67e-8
    R_total = (2/epsilon - 1) * (N + 1)
    q = sigma * (T_hot**4 - T_cold**4) / R_total
    return q

A.3 可视化输出

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

温度对比图：不同涂层在不同工况下的温度
热流分析图：MLI热流密度随层数变化
散热器性能图：散热功率随参数变化
在轨温度曲线：一个轨道周期内的温度变化

A.4 使用说明

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

pip install numpy matplotlib scipy

运行程序：

python run_simulation.py

"""
主题051：航天器热控系统辐射分析 - 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⁴
S0 = 1361  # 太阳常数 W/m²
R_E = 6371e3  # 地球半径 m

print("=" * 80)
print("主题051：航天器热控系统辐射分析 - Python仿真")
print("=" * 80)
print("\n本程序包含以下案例：")
print("  1. 航天器表面温度计算")
print("  2. 多层隔热材料性能分析")
print("  3. 辐射散热器设计优化")
print("  4. 航天器在轨热状态仿真")
print("\n" + "=" * 80)


def earth_view_factor(h, R_E=R_E):
    """
    计算地球对航天器的角系数
    
    参数:
    - h: 轨道高度 (m)
    - R_E: 地球半径 (m)
    
    返回:
    - F_12: 地球对航天器的角系数
    """
    sin_rho = R_E / (R_E + h)
    return sin_rho**2


def orbital_period(h, R_E=R_E):
    """
    计算圆轨道周期 (开普勒第三定律)
    
    参数:
    - h: 轨道高度 (m)
    - R_E: 地球半径 (m)
    
    返回:
    - period: 轨道周期 (s)
    """
    mu = 3.986e14  # 地球引力常数 m³/s²
    a = R_E + h  # 轨道半长轴
    period = 2 * np.pi * np.sqrt(a**3 / mu)
    return period


def leo_thermal_environment(h, albedo=0.3, T_E=255):
    """
    计算LEO轨道热环境
    
    参数:
    - h: 轨道高度 (m)
    - albedo: 地球反照率
    - T_E: 地球有效辐射温度 (K)
    
    返回:
    - q_s: 太阳辐射热流密度 (W/m²)
    - q_a: 地球反照热流密度 (W/m²)
    - q_IR: 地球红外辐射热流密度 (W/m²)
    - F_12: 地球角系数
    - period: 轨道周期 (s)
    - eclipse_ratio: 阴影期比例
    """
    # 太阳辐射
    q_s = S0
    
    # 地球角系数
    F_12 = earth_view_factor(h, R_E)
    
    # 地球反照 (简化计算)
    q_a = albedo * S0 * F_12
    
    # 地球红外辐射
    q_IR = SIGMA * T_E**4 * F_12
    
    # 轨道周期
    period = orbital_period(h, R_E)
    
    # 阴影期比例 (简化计算，假设圆轨道)
    # 阴影角 = 2*arcsin(R_E/(R_E+h))
    eclipse_angle = 2 * np.arcsin(R_E / (R_E + h))
    eclipse_ratio = eclipse_angle / (2 * np.pi)
    
    return q_s, q_a, q_IR, F_12, period, eclipse_ratio


def surface_temperature(alpha_s, epsilon, q_s, q_a, q_IR, q_int=0, in_sunlight=True):
    """
    计算表面热平衡温度
    
    参数:
    - alpha_s: 太阳吸收率
    - epsilon: 红外发射率
    - q_s: 太阳辐射热流密度 (W/m²)
    - q_a: 地球反照热流密度 (W/m²)
    - q_IR: 地球红外辐射热流密度 (W/m²)
    - q_int: 内部热流密度 (W/m²)
    - in_sunlight: 是否在光照期
    
    返回:
    - T: 表面温度 (K)
    """
    if in_sunlight:
        q_in = alpha_s * q_s + alpha_s * q_a + epsilon * q_IR + q_int
    else:
        # 阴影期无太阳辐射和地球反照
        q_in = epsilon * q_IR + q_int
    
    T = (q_in / (epsilon * SIGMA))**0.25
    return T


def mli_heat_flux(T_hot, T_cold, N, epsilon_screen=0.05):
    """
    计算通过多层隔热材料的热流密度
    
    参数:
    - T_hot: 热端温度 (K)
    - T_cold: 冷端温度 (K)
    - N: 反射屏层数
    - epsilon_screen: 反射屏发射率
    
    返回:
    - q: 热流密度 (W/m²)
    """
    # 多层隔热总热阻
    R_total = (2 / epsilon_screen - 1) * (N + 1)
    q = SIGMA * (T_hot**4 - T_cold**4) / R_total
    return q


def mli_effective_conductivity(T_hot, T_cold, N, epsilon_screen=0.05, thickness=0.01):
    """
    计算MLI等效热导率
    
    参数:
    - T_hot: 热端温度 (K)
    - T_cold: 冷端温度 (K)
    - N: 反射屏层数
    - epsilon_screen: 反射屏发射率
    - thickness: MLI总厚度 (m)
    
    返回:
    - lambda_eff: 等效热导率 (W/m·K)
    """
    q = mli_heat_flux(T_hot, T_cold, N, epsilon_screen)
    T_avg = (T_hot + T_cold) / 2
    lambda_eff = q * thickness / (T_hot - T_cold)
    return lambda_eff


def radiator_power(A, epsilon, T):
    """
    计算辐射散热器散热功率
    
    参数:
    - A: 散热器面积 (m²)
    - epsilon: 表面发射率
    - T: 散热器温度 (K)
    
    返回:
    - Q: 散热功率 (W)
    """
    return epsilon * SIGMA * A * T**4


def radiator_area(Q, epsilon, T):
    """
    计算所需散热器面积
    
    参数:
    - Q: 需要排散的热量 (W)
    - epsilon: 表面发射率
    - T: 散热器温度 (K)
    
    返回:
    - A: 所需面积 (m²)
    """
    return Q / (epsilon * SIGMA * T**4)


# ==============================================================================
# 案例1：航天器表面温度计算
# ==============================================================================
print("\n" + "=" * 80)
print("案例1：航天器表面温度计算")
print("=" * 80)

# 定义热控涂层
coatings = {
    'OSR': {'alpha_s': 0.10, 'epsilon': 0.85, 'color': '#E0F7FA'},
    '白漆': {'alpha_s': 0.25, 'epsilon': 0.88, 'color': '#FFFFFF'},
    '黑漆': {'alpha_s': 0.90, 'epsilon': 0.90, 'color': '#212121'},
    '镀金': {'alpha_s': 0.25, 'epsilon': 0.05, 'color': '#FFD700'},
    '阳极氧化铝': {'alpha_s': 0.35, 'epsilon': 0.85, 'color': '#C0C0C0'}
}

# 不同轨道高度
altitudes = [300e3, 400e3, 500e3, 600e3, 800e3]  # m

# 存储结果
results_temp = {}

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

# 1.1 LEO 400 km轨道温度计算
h = 400e3
q_s, q_a, q_IR, F_12, period, eclipse_ratio = leo_thermal_environment(h)

print(f"\n【LEO {h/1000:.0f} km轨道热环境】")
print(f"  太阳辐射: {q_s:.1f} W/m²")
print(f"  地球反照: {q_a:.1f} W/m²")
print(f"  地球红外: {q_IR:.1f} W/m²")
print(f"  地球角系数: {F_12:.4f}")
print(f"  轨道周期: {period/60:.1f} 分钟")
print(f"  阴影期比例: {eclipse_ratio*100:.1f}%")

# 计算不同涂层的温度
temps_sunlight = []
temps_eclipse = []
coating_names = []

for name, props in coatings.items():
    T_sun = surface_temperature(props['alpha_s'], props['epsilon'], 
                                q_s, q_a, q_IR, in_sunlight=True)
    T_ecl = surface_temperature(props['alpha_s'], props['epsilon'], 
                                q_s, q_a, q_IR, in_sunlight=False)
    temps_sunlight.append(T_sun - 273.15)  # 转换为°C
    temps_eclipse.append(T_ecl - 273.15)
    coating_names.append(name)
    
    print(f"\n  {name} (αs={props['alpha_s']:.2f}, ε={props['epsilon']:.2f}):")
    print(f"    光照期温度: {T_sun-273.15:.1f}°C")
    print(f"    阴影期温度: {T_ecl-273.15:.1f}°C")
    print(f"    温差: {T_sun-T_ecl:.1f} K")

# 绘制温度对比图
ax1 = axes[0, 0]
x = np.arange(len(coating_names))
width = 0.35
bars1 = ax1.bar(x - width/2, temps_sunlight, width, label='光照期', color='#FFB74D', edgecolor='black')
bars2 = ax1.bar(x + width/2, temps_eclipse, width, label='阴影期', color='#64B5F6', edgecolor='black')
ax1.set_xlabel('热控涂层', fontsize=11)
ax1.set_ylabel('温度 (°C)', fontsize=11)
ax1.set_title(f'LEO {h/1000:.0f} km轨道不同涂层温度对比', fontsize=12, fontweight='bold')
ax1.set_xticks(x)
ax1.set_xticklabels(coating_names, rotation=15, ha='right')
ax1.legend()
ax1.grid(axis='y', alpha=0.3)
ax1.axhline(y=0, color='gray', linestyle='--', alpha=0.5)

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

# 1.2 轨道高度对温度的影响
ax2 = axes[0, 1]
altitudes_km = np.array(altitudes) / 1000

for name, props in coatings.items():
    temps_sun = []
    temps_ecl = []
    for h in altitudes:
        q_s, q_a, q_IR, F_12, period, eclipse_ratio = leo_thermal_environment(h)
        T_sun = surface_temperature(props['alpha_s'], props['epsilon'], 
                                    q_s, q_a, q_IR, in_sunlight=True)
        T_ecl = surface_temperature(props['alpha_s'], props['epsilon'], 
                                    q_s, q_a, q_IR, in_sunlight=False)
        temps_sun.append(T_sun - 273.15)
        temps_ecl.append(T_ecl - 273.15)
    
    ax2.plot(altitudes_km, temps_sun, 'o-', label=f'{name} (光照)', 
             color=props['color'], markeredgecolor='black', markersize=6)
    ax2.plot(altitudes_km, temps_ecl, 's--', label=f'{name} (阴影)', 
             color=props['color'], markeredgecolor='black', markersize=5, alpha=0.6)

ax2.set_xlabel('轨道高度 (km)', fontsize=11)
ax2.set_ylabel('温度 (°C)', fontsize=11)
ax2.set_title('不同轨道高度下的表面温度', fontsize=12, fontweight='bold')
ax2.legend(fontsize=8, ncol=2)
ax2.grid(alpha=0.3)
ax2.axhline(y=0, color='gray', linestyle='--', alpha=0.5)

# 1.3 吸收-发射比与温度的关系
ax3 = axes[1, 0]
alpha_range = np.linspace(0.05, 0.95, 50)
epsilon_values = [0.05, 0.3, 0.6, 0.9]

for eps in epsilon_values:
    temps = []
    for alpha in alpha_range:
        T = surface_temperature(alpha, eps, q_s, q_a, q_IR, in_sunlight=True)
        temps.append(T - 273.15)
    ax3.plot(alpha_range/eps, temps, '-', label=f'ε = {eps:.2f}', linewidth=2)

ax3.set_xlabel('吸收-发射比 (αs/ε)', fontsize=11)
ax3.set_ylabel('温度 (°C)', fontsize=11)
ax3.set_title('吸收-发射比对表面温度的影响', fontsize=12, fontweight='bold')
ax3.legend()
ax3.grid(alpha=0.3)
ax3.axhline(y=0, color='gray', linestyle='--', alpha=0.5)
ax3.axvline(x=1.0, color='red', linestyle=':', alpha=0.5, label='αs/ε = 1')

# 1.4 内部热源对温度的影响
ax4 = axes[1, 1]
q_int_range = np.linspace(0, 500, 50)  # W/m²

for name, props in coatings.items():
    temps = []
    for q_int in q_int_range:
        T = surface_temperature(props['alpha_s'], props['epsilon'], 
                               q_s, q_a, q_IR, q_int=q_int, in_sunlight=True)
        temps.append(T - 273.15)
    ax4.plot(q_int_range, temps, '-', label=name, linewidth=2, color=props['color'])

ax4.set_xlabel('内部热流密度 (W/m²)', fontsize=11)
ax4.set_ylabel('温度 (°C)', fontsize=11)
ax4.set_title('内部热源对表面温度的影响', fontsize=12, fontweight='bold')
ax4.legend()
ax4.grid(alpha=0.3)
ax4.axhline(y=0, color='gray', linestyle='--', alpha=0.5)

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

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


# ==============================================================================
# 案例2：多层隔热材料性能分析
# ==============================================================================
print("\n" + "=" * 80)
print("案例2：多层隔热材料性能分析")
print("=" * 80)

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

# 2.1 层数对热流密度的影响
ax1 = axes[0, 0]
N_range = np.arange(1, 51)
T_hot = 300  # K
T_cold = 100  # K
epsilon_screen = 0.05

heat_fluxes = [mli_heat_flux(T_hot, T_cold, N, epsilon_screen) for N in N_range]

ax1.semilogy(N_range, heat_fluxes, 'b-', linewidth=2)
ax1.set_xlabel('反射屏层数', fontsize=11)
ax1.set_ylabel('热流密度 (W/m²)', fontsize=11)
ax1.set_title(f'MLI热流密度随层数变化 (ε={epsilon_screen})', fontsize=12, fontweight='bold')
ax1.grid(alpha=0.3, which='both')

# 标注关键点
for N in [5, 10, 20, 30]:
    q = mli_heat_flux(T_hot, T_cold, N, epsilon_screen)
    ax1.plot(N, q, 'ro', markersize=8)
    ax1.annotate(f'N={N}\nq={q:.1f} W/m²', xy=(N, q), xytext=(N+3, q*1.5),
                fontsize=9, arrowprops=dict(arrowstyle='->', color='red'))
    print(f"  {N}层MLI热流密度: {q:.2f} W/m²")

# 2.2 等效热导率
ax2 = axes[0, 1]
thickness = 0.02  # m (2 cm)
lambda_effs = [mli_effective_conductivity(T_hot, T_cold, N, epsilon_screen, thickness) 
               for N in N_range]

ax2.semilogy(N_range, lambda_effs, 'g-', linewidth=2)
ax2.set_xlabel('反射屏层数', fontsize=11)
ax2.set_ylabel('等效热导率 (W/m·K)', fontsize=11)
ax2.set_title('MLI等效热导率随层数变化', fontsize=12, fontweight='bold')
ax2.grid(alpha=0.3, which='both')

# 对比其他材料
materials = {
    '静止空气': 0.026,
    '聚氨酯泡沫': 0.025,
    '不锈钢': 15
}
for name, lambda_val in materials.items():
    ax2.axhline(y=lambda_val, color='red', linestyle='--', alpha=0.5)
    ax2.text(35, lambda_val*1.2, name, fontsize=9, color='red')

# 2.3 反射屏发射率的影响
ax3 = axes[1, 0]
epsilon_range = np.linspace(0.02, 0.15, 50)
N_values = [5, 10, 20, 30]

for N in N_values:
    heat_fluxes_eps = [mli_heat_flux(T_hot, T_cold, N, eps) for eps in epsilon_range]
    ax3.semilogy(epsilon_range, heat_fluxes_eps, '-', label=f'N={N}', linewidth=2)

ax3.set_xlabel('反射屏发射率', fontsize=11)
ax3.set_ylabel('热流密度 (W/m²)', fontsize=11)
ax3.set_title('反射屏发射率对隔热性能的影响', fontsize=12, fontweight='bold')
ax3.legend()
ax3.grid(alpha=0.3, which='both')

# 2.4 温差对热流密度的影响
ax4 = axes[1, 1]
T_cold_fixed = 100  # K
T_hot_range = np.linspace(150, 400, 50)
N_values = [5, 10, 20, 30]

for N in N_values:
    heat_fluxes_temp = [mli_heat_flux(T_hot, T_cold_fixed, N, epsilon_screen) 
                        for T_hot in T_hot_range]
    ax4.plot(T_hot_range - 273.15, heat_fluxes_temp, '-', label=f'N={N}', linewidth=2)

ax4.set_xlabel('热端温度 (°C)', fontsize=11)
ax4.set_ylabel('热流密度 (W/m²)', fontsize=11)
ax4.set_title('温差对MLI热流密度的影响', fontsize=12, fontweight='bold')
ax4.legend()
ax4.grid(alpha=0.3)

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

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

print(f"\n【关键计算结果】")
print(f"20层MLI热流密度: {mli_heat_flux(300, 100, 20, 0.05):.2f} W/m²")
print(f"等效热导率: {mli_effective_conductivity(300, 100, 20, 0.05, 0.02):.2e} W/m·K")
print(f"相比静止空气改善: {0.026/mli_effective_conductivity(300, 100, 20, 0.05, 0.02):.0f}倍")


# ==============================================================================
# 案例3：辐射散热器设计优化
# ==============================================================================
print("\n" + "=" * 80)
print("案例3：辐射散热器设计优化")
print("=" * 80)

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

# 3.1 散热器面积与散热功率
ax1 = axes[0, 0]
A_range = np.linspace(0.5, 5, 50)  # m²
T_rad = 320  # K
epsilon_rad = 0.88

powers = [radiator_power(A, epsilon_rad, T_rad) for A in A_range]

ax1.plot(A_range, powers, 'b-', linewidth=2)
ax1.set_xlabel('散热器面积 (m²)', fontsize=11)
ax1.set_ylabel('散热功率 (W)', fontsize=11)
ax1.set_title(f'散热器面积与散热功率 (T={T_rad}K, ε={epsilon_rad})', 
              fontsize=12, fontweight='bold')
ax1.grid(alpha=0.3)

# 标注关键点
for A in [1, 2, 3, 4]:
    Q = radiator_power(A, epsilon_rad, T_rad)
    ax1.plot(A, Q, 'ro', markersize=8)
    ax1.annotate(f'{Q:.0f}W', xy=(A, Q), xytext=(A+0.1, Q+50),
                fontsize=9)
    print(f"  面积{A:.0f}m²散热器功率: {Q:.0f} W")

# 3.2 散热器温度与散热功率
ax2 = axes[0, 1]
T_range = np.linspace(250, 400, 50)  # K
A_fixed = 2.0  # m²

powers_temp = [radiator_power(A_fixed, epsilon_rad, T) for T in T_range]

ax2.plot(T_range - 273.15, powers_temp, 'r-', linewidth=2)
ax2.set_xlabel('散热器温度 (°C)', fontsize=11)
ax2.set_ylabel('散热功率 (W)', fontsize=11)
ax2.set_title(f'散热器温度与散热功率 (A={A_fixed}m²)', 
              fontsize=12, fontweight='bold')
ax2.grid(alpha=0.3)

# 标注T^4关系
ax2.annotate('Q ∝ T⁴', xy=(100, 1500), fontsize=14, color='blue', fontweight='bold')

# 3.3 发射率对散热功率的影响
ax3 = axes[1, 0]
epsilon_range = np.linspace(0.5, 0.95, 50)
A_values = [1.0, 2.0, 3.0, 4.0]

for A in A_values:
    powers_eps = [radiator_power(A, eps, T_rad) for eps in epsilon_range]
    ax3.plot(epsilon_range, powers_eps, '-', label=f'A={A}m²', linewidth=2)

ax3.set_xlabel('表面发射率', fontsize=11)
ax3.set_ylabel('散热功率 (W)', fontsize=11)
ax3.set_title(f'发射率对散热功率的影响 (T={T_rad}K)', fontsize=12, fontweight='bold')
ax3.legend()
ax3.grid(alpha=0.3)

# 3.4 散热器朝向优化（考虑外部热流）
ax4 = axes[1, 1]

# 散热器朝向角度（0°=面向太阳，90°=侧向太阳，180°=背向太阳）
angles = np.linspace(0, 180, 50)  # 度

# 计算不同朝向下接收的外部热流
h = 400e3
q_s, q_a, q_IR, F_12, period, eclipse_ratio = leo_thermal_environment(h)

# 简化模型：散热器接收的热流随角度变化
# 假设散热器法线与太阳夹角为theta
q_external = []
for theta_deg in angles:
    theta = np.radians(theta_deg)
    # 太阳辐射分量
    q_surf = q_s * np.cos(theta) if theta < np.pi/2 else 0
    # 地球辐射（假设始终可见）
    q_earth = (q_a + q_IR) * F_12
    q_external.append(q_surf + q_earth)

# 计算净散热功率（假设散热器温度350K，发射率0.88）
T_rad_opt = 350  # K
epsilon_opt = 0.88
A_opt = 2.0  # m²

net_powers = []
for i, theta_deg in enumerate(angles):
    # 散热器发射功率
    Q_emit = epsilon_opt * SIGMA * A_opt * T_rad_opt**4
    # 吸收的外部热流
    Q_abs = 0.25 * q_external[i] * A_opt  # 假设太阳吸收率0.25
    # 净散热
    Q_net = Q_emit - Q_abs
    net_powers.append(Q_net)

ax4.plot(angles, net_powers, 'g-', linewidth=2)
ax4.set_xlabel('散热器法线与太阳光线夹角 (°)', fontsize=11)
ax4.set_ylabel('净散热功率 (W)', fontsize=11)
ax4.set_title('散热器朝向对散热性能的影响', fontsize=12, fontweight='bold')
ax4.grid(alpha=0.3)
ax4.axvline(x=90, color='red', linestyle='--', alpha=0.5, label='侧向太阳')
ax4.axvline(x=180, color='blue', linestyle='--', alpha=0.5, label='背向太阳')
ax4.legend()

# 标注最佳朝向
max_idx = np.argmax(net_powers)
ax4.plot(angles[max_idx], net_powers[max_idx], 'r*', markersize=15)
ax4.annotate(f'最佳朝向\n{angles[max_idx]:.0f}°\n{net_powers[max_idx]:.0f}W', 
            xy=(angles[max_idx], net_powers[max_idx]), 
            xytext=(angles[max_idx]-40, net_powers[max_idx]+100),
            fontsize=10, arrowprops=dict(arrowstyle='->', color='red'))

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

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

print(f"\n【关键计算结果】")
print(f"散热器面积2m²、温度320K、发射率0.88时的散热功率: {radiator_power(2, 0.88, 320):.0f} W")
print(f"排散1000W所需面积 (T=320K, ε=0.88): {radiator_area(1000, 0.88, 320):.2f} m²")


# ==============================================================================
# 案例4：航天器在轨热状态仿真
# ==============================================================================
print("\n" + "=" * 80)
print("案例4：航天器在轨热状态仿真")
print("=" * 80)

# 航天器参数
spacecraft_mass = 500  # kg
spacecraft_cp = 800  # J/kg·K (复合材料比热容)
spacecraft_area = 6  # m² (1m×1m立方体表面积)
internal_heat = 500  # W

# 轨道参数
h = 400e3  # m
q_s, q_a, q_IR, F_12, period, eclipse_ratio = leo_thermal_environment(h)

# 表面特性（简化模型）
alpha_s_eff = 0.20  # 等效太阳吸收率
epsilon_eff = 0.85  # 等效红外发射率

# 时间参数
dt = 10  # s
t_total = period * 2  # 仿真2个轨道周期
t = np.arange(0, t_total, dt)
n_steps = len(t)

# 初始化温度
T_initial = 273.15  # K (0°C)
T = np.zeros(n_steps)
T[0] = T_initial

# 热容
C = spacecraft_mass * spacecraft_cp

# 瞬态热平衡计算
print("\n  正在进行瞬态热平衡计算...")
for i in range(1, n_steps):
    # 当前时刻在轨道周期中的位置
    t_orbit = t[i] % period
    
    # 判断是否在阴影期
    eclipse_duration = period * eclipse_ratio
    in_eclipse = t_orbit < eclipse_duration
    
    # 计算外部热流
    if in_eclipse:
        Q_external = epsilon_eff * q_IR * spacecraft_area * F_12
    else:
        Q_solar = alpha_s_eff * q_s * spacecraft_area * 0.25  # 假设25%面积受照
        Q_albedo = alpha_s_eff * q_a * spacecraft_area * F_12
        Q_IR = epsilon_eff * q_IR * spacecraft_area * F_12
        Q_external = Q_solar + Q_albedo + Q_IR
    
    # 辐射散热
    Q_radiation = epsilon_eff * SIGMA * spacecraft_area * T[i-1]**4
    
    # 热平衡方程
    dT_dt = (Q_external + internal_heat - Q_radiation) / C
    T[i] = T[i-1] + dT_dt * dt

print(f"  计算完成，温度范围: {np.min(T)-273.15:.1f}°C 至 {np.max(T)-273.15:.1f}°C")

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

# 4.1 温度随时间变化
ax1 = axes[0, 0]
t_minutes = t / 60
ax1.plot(t_minutes, T - 273.15, 'b-', linewidth=1.5)
ax1.set_xlabel('时间 (分钟)', fontsize=11)
ax1.set_ylabel('温度 (°C)', fontsize=11)
ax1.set_title('航天器在轨温度变化 (2个轨道周期)', fontsize=12, fontweight='bold')
ax1.grid(alpha=0.3)
ax1.axhline(y=0, color='gray', linestyle='--', alpha=0.5)

# 标注阴影期
for i in range(2):
    start = i * period / 60
    end = start + period * eclipse_ratio / 60
    ax1.axvspan(start, end, alpha=0.2, color='blue', label='阴影期' if i == 0 else '')
ax1.legend()

# 标注温度极值
T_max_idx = np.argmax(T)
T_min_idx = np.argmin(T)
ax1.plot(t_minutes[T_max_idx], T[T_max_idx]-273.15, 'r^', markersize=10)
ax1.plot(t_minutes[T_min_idx], T[T_min_idx]-273.15, 'bv', markersize=10)
ax1.annotate(f'最高: {T[T_max_idx]-273.15:.1f}°C', 
            xy=(t_minutes[T_max_idx], T[T_max_idx]-273.15),
            xytext=(t_minutes[T_max_idx]+5, T[T_max_idx]-273.15+5),
            fontsize=9, arrowprops=dict(arrowstyle='->', color='red'))
ax1.annotate(f'最低: {T[T_min_idx]-273.15:.1f}°C', 
            xy=(t_minutes[T_min_idx], T[T_min_idx]-273.15),
            xytext=(t_minutes[T_min_idx]+5, T[T_min_idx]-273.15-10),
            fontsize=9, arrowprops=dict(arrowstyle='->', color='blue'))

# 4.2 一个轨道周期内的温度变化（平均）
ax2 = axes[0, 1]
n_orbits = int(t_total / period)
period_steps = int(period/dt)
T_avg = np.zeros(period_steps)
for i in range(n_orbits):
    start_idx = int(i * period / dt)
    end_idx = start_idx + period_steps
    if end_idx <= len(T):
        T_avg += T[start_idx:end_idx]
T_avg /= n_orbits

t_orbit = np.arange(0, period, dt) / 60
# 确保数组维度一致
min_len = min(len(t_orbit), len(T_avg))
ax2.plot(t_orbit[:min_len], T_avg[:min_len] - 273.15, 'g-', linewidth=2)
ax2.set_xlabel('轨道周期内时间 (分钟)', fontsize=11)
ax2.set_ylabel('平均温度 (°C)', fontsize=11)
ax2.set_title('一个轨道周期内的平均温度变化', fontsize=12, fontweight='bold')
ax2.grid(alpha=0.3)
ax2.axhline(y=0, color='gray', linestyle='--', alpha=0.5)

# 标注光照期和阴影期
eclipse_min = period * eclipse_ratio / 60
ax2.axvspan(0, eclipse_min, alpha=0.2, color='blue', label='阴影期')
ax2.axvspan(eclipse_min, period/60, alpha=0.2, color='yellow', label='光照期')
ax2.legend()

# 4.3 不同内部热功率下的温度
ax3 = axes[1, 0]
internal_powers = [200, 500, 800, 1200]  # W
colors = ['blue', 'green', 'orange', 'red']

for Q_int, color in zip(internal_powers, colors):
    T_sim = np.zeros(n_steps)
    T_sim[0] = T_initial
    
    for i in range(1, n_steps):
        t_orbit = t[i] % period
        eclipse_duration = period * eclipse_ratio
        in_eclipse = t_orbit < eclipse_duration
        
        if in_eclipse:
            Q_external = epsilon_eff * q_IR * spacecraft_area * F_12
        else:
            Q_solar = alpha_s_eff * q_s * spacecraft_area * 0.25
            Q_albedo = alpha_s_eff * q_a * spacecraft_area * F_12
            Q_IR = epsilon_eff * q_IR * spacecraft_area * F_12
            Q_external = Q_solar + Q_albedo + Q_IR
        
        Q_radiation = epsilon_eff * SIGMA * spacecraft_area * T_sim[i-1]**4
        dT_dt = (Q_external + Q_int - Q_radiation) / C
        T_sim[i] = T_sim[i-1] + dT_dt * dt
    
    ax3.plot(t_minutes, T_sim - 273.15, '-', color=color, 
            label=f'Q_int={Q_int}W', linewidth=1.5)

ax3.set_xlabel('时间 (分钟)', fontsize=11)
ax3.set_ylabel('温度 (°C)', fontsize=11)
ax3.set_title('不同内部热功率下的温度变化', fontsize=12, fontweight='bold')
ax3.legend()
ax3.grid(alpha=0.3)
ax3.axhline(y=0, color='gray', linestyle='--', alpha=0.5)

# 4.4 不同涂层（发射率）下的温度
ax4 = axes[1, 1]
epsilon_values = [0.6, 0.75, 0.85, 0.92]

for eps, color in zip(epsilon_values, colors):
    T_sim = np.zeros(n_steps)
    T_sim[0] = T_initial
    
    for i in range(1, n_steps):
        t_orbit = t[i] % period
        eclipse_duration = period * eclipse_ratio
        in_eclipse = t_orbit < eclipse_duration
        
        if in_eclipse:
            Q_external = eps * q_IR * spacecraft_area * F_12
        else:
            Q_solar = alpha_s_eff * q_s * spacecraft_area * 0.25
            Q_albedo = alpha_s_eff * q_a * spacecraft_area * F_12
            Q_IR = eps * q_IR * spacecraft_area * F_12
            Q_external = Q_solar + Q_albedo + Q_IR
        
        Q_radiation = eps * SIGMA * spacecraft_area * T_sim[i-1]**4
        dT_dt = (Q_external + internal_heat - Q_radiation) / C
        T_sim[i] = T_sim[i-1] + dT_dt * dt
    
    ax4.plot(t_minutes, T_sim - 273.15, '-', color=color, 
            label=f'ε={eps:.2f}', linewidth=1.5)

ax4.set_xlabel('时间 (分钟)', fontsize=11)
ax4.set_ylabel('温度 (°C)', fontsize=11)
ax4.set_title('不同发射率下的温度变化', fontsize=12, fontweight='bold')
ax4.legend()
ax4.grid(alpha=0.3)
ax4.axhline(y=0, color='gray', linestyle='--', alpha=0.5)

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

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

print(f"\n【关键计算结果】")
print(f"轨道周期: {period/60:.1f} 分钟")
print(f"阴影期比例: {eclipse_ratio*100:.1f}%")
print(f"温度波动范围: {np.min(T)-273.15:.1f}°C 至 {np.max(T)-273.15:.1f}°C")
print(f"平均温度: {np.mean(T)-273.15:.1f}°C")


# ==============================================================================
# 生成GIF动画：航天器在轨温度变化
# ==============================================================================
print("\n" + "=" * 80)
print("生成GIF动画：航天器在轨温度变化")
print("=" * 80)

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

# 准备动画数据
n_frames = 100
t_anim = np.linspace(0, period, n_frames)
T_anim = np.interp(t_anim, t[:int(period/dt)], T[:int(period/dt)])

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

def update(frame):
    ax1.clear()
    ax2.clear()
    
    current_time = t_anim[frame]
    current_temp = T_anim[frame]
    
    # 左图：温度-时间曲线
    ax1.plot(t_anim[:frame+1]/60, T_anim[:frame+1]-273.15, 'b-', linewidth=2)
    ax1.scatter([current_time/60], [current_temp-273.15], c='red', s=100, zorder=5)
    ax1.set_xlim(0, period/60)
    ax1.set_ylim(np.min(T)-273.15-5, np.max(T)-273.15+5)
    ax1.set_xlabel('时间 (分钟)', fontsize=11)
    ax1.set_ylabel('温度 (°C)', fontsize=11)
    ax1.set_title('航天器在轨温度变化', fontsize=12, fontweight='bold')
    ax1.grid(alpha=0.3)
    ax1.axhline(y=0, color='gray', linestyle='--', alpha=0.5)
    
    # 标注阴影期
    eclipse_min = period * eclipse_ratio / 60
    ax1.axvspan(0, eclipse_min, alpha=0.2, color='blue', label='阴影期')
    ax1.axvspan(eclipse_min, period/60, alpha=0.2, color='yellow', label='光照期')
    
    # 右图：航天器示意图
    # 绘制立方体
    from mpl_toolkits.mplot3d import Axes3D
    ax2.set_xlim(-1.5, 1.5)
    ax2.set_ylim(-1.5, 1.5)
    ax2.set_aspect('equal')
    
    # 根据温度设置颜色
    temp_norm = (current_temp - np.min(T)) / (np.max(T) - np.min(T))
    color = plt.cm.RdYlBu_r(temp_norm)
    
    # 绘制简化的航天器
    rect = plt.Rectangle((-0.5, -0.5), 1, 1, facecolor=color, edgecolor='black', linewidth=2)
    ax2.add_patch(rect)
    
    # 绘制太阳
    sun_x = 1.2
    sun_y = 1.0
    sun = plt.Circle((sun_x, sun_y), 0.25, color='yellow', ec='orange', linewidth=2)
    ax2.add_patch(sun)
    ax2.annotate('太阳', xy=(sun_x, sun_y), ha='center', va='center', fontsize=10, fontweight='bold')
    
    # 绘制地球
    earth = plt.Circle((0, -1.8), 0.6, color='blue', ec='darkblue', linewidth=2)
    ax2.add_patch(earth)
    ax2.annotate('地球', xy=(0, -1.8), ha='center', va='center', fontsize=10, fontweight='bold', color='white')
    
    # 判断是否在阴影期
    in_eclipse = current_time < period * eclipse_ratio
    phase = "阴影期" if in_eclipse else "光照期"
    
    ax2.set_title(f'航天器热状态: {phase}\n温度: {current_temp-273.15:.1f}°C', 
                 fontsize=12, fontweight='bold')
    ax2.axis('off')
    
    # 添加图例
    sm = plt.cm.ScalarMappable(cmap=plt.cm.RdYlBu_r, 
                               norm=plt.Normalize(vmin=np.min(T)-273.15, vmax=np.max(T)-273.15))
    sm.set_array([])
    cbar = plt.colorbar(sm, ax=ax2, shrink=0.6)
    cbar.set_label('温度 (°C)', fontsize=10)
    
    return []

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

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


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

print("\n生成的文件：")
print("  - case1_spacecraft_temperature.png")
print("  - case2_mli_performance.png")
print("  - case3_radiator_design.png")
print("  - case4_orbit_thermal_simulation.png")
print("  - spacecraft_thermal_animation.gif")

print("\n" + "=" * 80)
print("主题051：航天器热控系统辐射分析 - 仿真完成")
print("=" * 80)