球壳的定常热传导
球壳定常热传导的理论基础
球壳与圆筒的区别
老师，球壳的定常热传导与圆筒的情况有什么不同？两者都是"圆形"，好像应该很相似啊。
最大的差别是传热面积的增长方式。在圆筒中，面积与半径成正比（$A = 2\pi r L$），但在球体中，面积与半径的平方成正比（$A = 4\pi r^2$）。
面积增长方式不同，温度分布也会改变吗？
完全正确。在圆筒中，温度变化为$\ln r$（对数函数），而在球壳中为$1/r$。直观地说，球壳向外侧扩展时面积增长更快，所以温度降低衰减也更快。
| 形状 | 面积 $A(r)$ | 温度分布 | 热阻的形式 |
|---|---|---|---|
| 平板 | 常数 | $T \propto r$（线性） | $L/(kA)$ |
| 圆筒 | $\propto r$ | $T \propto \ln r$ | $\ln(r_2/r_1)/(2\pi k L)$ |
| 球壳 | $\propto r^2$ | $T \propto 1/r$ | $(1/r_1 - 1/r_2)/(4\pi k)$ |
也就是说，面积的"幂次"越高，温度降低就越缓和？那在实际工程中，球壳主要用在哪些场合呢？
主要有三种应用场景。
- 液化天然气球形罐（直径40米级）的隔热设计 ― 预测蒸发气体损失
- 核反应堆压力容器的肉厚方向温度分布 ― 热应力评估的输入
- 球形燃料芯块（高温气冷堆HTGR的TRISO燃料）的中心温度估算
这些场合都能使用球对称假设，通过一维解析解快速给出设计的初步估计，这正是球壳模型的优势所在。
支配方程与解析解
请教我球壳的支配方程。圆筒的时候有$1/r$项对吧。
球坐标中，定常、径向一维热传导方程是这样的。
与圆筒的$\frac{1}{r}\frac{d}{dr}\left(kr\frac{dT}{dr}\right)=0$相比，$r$的指数从1变成了2。这正反映了面积与$r^2$成正比的事实。
如果$k$是常数，怎么积分呢？
当$k$为常数时，$\frac{d}{dr}\left(r^2 \frac{dT}{dr}\right) = 0$积分两次得
代入边界条件 $T(r_1) = T_1$、$T(r_2) = T_2$，即可得到温度分布的解析解。
可以看作圆筒中的$\ln r$分布被替换成了$1/r$的形式，这是圆筒和球壳的直观对应。
$1/r$分布具体的图形是什么样呢？
举例说，内径100mm、外径200mm的钢球壳，$T_1 = 500$°C、$T_2 = 100$°C时，
- $r = 100$ mm: $T = 500$°C（内表面）
- $r = 133$ mm（肉厚1/3处）：$T = 300$°C
- $r = 200$ mm: $T = 100$°C（外表面）
对于平板，肉厚1/3处的温度为$(500-100)\times(1-1/3)+100 = 367$°C，可见球壳在内表面附近的温度降低更快。这也意味着内表面附近的热流密度更大。
球壳的热阻
通过球壳的热流量怎么计算？能像圆筒那样用热阻的概念吗？
完全可以。用傅里叶定律 $q = -kA(r)\frac{dT}{dr}$，代入$A(r) = 4\pi r^2$并积分
用$Q = \Delta T / R$的形式整理，球壳的导热热阻为
与圆筒不同的是，这里不含长度$L$。这是因为球是全向闭合的形状，不需要管长这样的参数。
如果内外表面有对流，直接加上对流热阻就行吗？
对的。球面的对流热阻为 $R_{\text{conv}} = 1/(4\pi r^2 h)$，从内表面流体温度$T_{\infty,1}$到外表面流体温度$T_{\infty,2}$
这就是对流+导热+对流的串联电路，与欧姆定律的结构完全相同。
多层球壳模型
液化天然气罐那样夹着隔热材料的多层情况，怎么计算呢？
对于$n$层球壳，各层热导率为$k_i$，半径为$r_0, r_1, \ldots, r_n$，总导热热阻为
包含对流的总传热系数$U$，以内表面$A_1 = 4\pi r_0^2$为基准面积
这样就能参数化地计算隔热材料厚度变化时的热损失了。
正是这样。以直径40米的液化天然气球形罐为例，钢壳30mm（$k = 50$ W/(m·K)）加珍珠岩隔热材500mm（$k = 0.04$ W/(m·K)），隔热材的热阻约为钢壳的2500倍。用热阻法能定量识别哪一层是控制因素，这正是热阻法的强大之处。
内部发热的情况
像核反应堆燃料球那样内部发热的情况，该怎么处理？
有均匀体积热源$\dot{q}$ [W/m³]的实心球（半径$R$、表面温度$T_s$），支配方程为
在中心对称条件 $dT/dr|_{r=0} = 0$ 和$T(R) = T_s$下求解
中心温度为 $T_{\max} = T_s + \dot{q} R^2 / (6k)$。与圆筒的情况（分母为$4k$）相比，球的中心温度上升更小，因为表面积更大散热更快。
直径约1mm的高温气冷堆TRISO燃料粒子，用这个式子就能估算中心温度了。
正确。以UO₂内核为例（$k \approx 3$ W/(m·K)、$\dot{q} \approx 2 \times 10^8$ W/m³、$R = 0.25$ mm），$\Delta T_{\max} = 2\times10^8 \times (0.25\times10^{-3})^2 / (6 \times 3) \approx 0.7$°C。粒子很小，所以温度升幅很小。但实际上需要加上多层被覆（如SiC层）的接触热阻。
球壳热传导与地球科学
球壳定常热传导理论的渊源可追溯到19世纪的泊松（1820年代）和开尔文勋爵。开尔文曾把地球视为一个冷却的球体，试图用定常热传导的$1/r$分布从地球年龄推断出来（约1亿年），但他没有考虑放射性衰变导致的内部热源，因此严重低估了。现代地球物理仍在用球壳热传导作为地幔对流模型的初始近似。
球壳定常热传导的数值计算方法
轴对称FEM模型
球壳热传导的有限元分析，一定要建立3D球体模型吗？
如果有球对称性，用二维轴对称模型就足够了。这相当于画出断面的半圆，然后围绕对称轴旋转。计算成本比3D模型低1000倍以上。
| 模型类型 | 单元数量 | 自由度数 | 计算时间 |
|---|---|---|---|
| 3D全体（整个球） | 50万以上 | 50万以上 | 分钟级 |
| 3D 1/8对称 | 6万以上 | 6万以上 | 秒至分钟 |
| 2D轴对称 | 200以上 | 200以上 | 瞬时 |
| 1D解析解 | 无 | 无 | 手算 |
用轴对称模型要注意什么？
有三点需要注意。
- 单元类型的选择：Ansys Mechanical选用PLANE55（线性）/PLANE77（二阶），KEYOPT(3)=1（轴对称）。Abaqus选用DCAX4/DCAX8。
- 对称轴上的边界条件：$r=0$轴上自动应用绝热条件（$\partial T/\partial r = 0$），但某些软件需要显式设置。
- 网格偏置：内表面温度梯度陡，要细分；外表面梯度缓，可粗化。通常内表面单元尺寸设为外表面的1/3到1/5。
有限差分法的离散化
如果不用有限元，想用有限差分自己编程（用Python），该怎么做？
球坐标定常热传导方程展开后
用中心差分离散化，在节点$i$（位置$r_i$）
整理成三对角矩阵的线性方程组
$r_i$越小（接近内表面），$\Delta r / r_i$的影响越大，所以等间距网格在内表面附近精度下降。用非均匀网格或变量变换$s = 1/r$可以改善精度。
$s = 1/r$的变换很有趣，变换后会怎样？
做变换$s = 1/r$，球壳方程变成$\frac{d^2T}{ds^2} = 0$。也就是说，在$s$坐标空间中，温度分布变成线性的。用等间距的$s$网格离散化，中心差分会给出精确解。这是球壳热传导特有的技巧。
网格设计指南
有限元网格的单元尺寸怎么决定？
球壳定常热传导的网格设计参考如下。
| 项目 | 推荐值 | 理由 |
|---|---|---|
| 肉厚方向分层数 | 10层以上 | 准确再现$1/r$分布 |
| 内表面偏置比 | 1:3～1:5 | 内表面温度梯度陡峭 |
| 单元类型 | 二阶单元推荐 | $1/r$分布用线性单元近似较粗糙 |
| 多层界面 | 节点重合 | 界面处温度梯度不连续 |
| 宽高比 | 5以下 | 过度扁平的单元精度下降 |
怎么验证网格收敛性？
用三个或以上不同密度的网格进行计算，检查关键点的温度是否收敛。因为球壳有解析解，直接与解析解对比是最有效的验证方法。通常肉厚方向5、10、20层进行计算，比较内表面热流，变化率小于1%即可判定为收敛。
验证用基准问题
有现成的标准基准问题用来验证我的有限元模型是否正确吗？
有个典型的基准问题。
| 参数 | 数值 |
|---|---|
| 内径 $r_1$ | 0.1 m |
| 外径 $r_2$ | 0.2 m |
| 热导率 $k$ | 50 W/(m·K) |
| 内表面温度 $T_1$ | 500°C |
| 外表面温度 $T_2$ | 100°C |
解析解：$r = 0.15$ m处的温度
热流量 $Q = 4\pi \times 50 \times 400 / (1/0.1 - 1/0.2) = 50{,}265$ W。有限元结果与解析解相差在0.1%以内即通过验证。
球壳定常热传导的实务应用
液化天然气罐的隔热设计
液化天然气球形罐的隔热设计中，球壳热传导计算实际是怎么用的？
液化天然气球形罐（Moss型，直径40米级）的设计要求日蒸发气体（BOG）损失限制在贮存容量的0.05%以下。用多层球壳热传导计算来确定满足这一要求的隔热材厚度。
- 钢壳：9% Ni钢，厚30mm，$k = 15$ W/(m·K)
- 隔热材：珍珠岩粉末，厚500mm，$k = 0.04$ W/(m·K)
- 内表面条件：液化天然气温度-162°C，$h_1 = 500$ W/(m²·K)（沸腾）
- 外表面条件：大气温度+35°C，$h_2 = 10$ W/(m²·K)（自然对流）
用多层球壳公式计算，侵入热量约80 kW，BOG率约0.04%，满足要求。如果隔热材减为400mm，侵入热量超100 kW，会超过要求值。
手算出80kW后，还有必要再用有限元验证吗？
实际的罐不是完全球对称的。支撑结构（裙座）、焊缝处隔热材的缝隙、上下温度差导致的自然对流等三维效应，一维模型无法捕捉。有限元的作用是评估这些局部效应的影响。手算是"设计框架"，有限元是"详细检讨"，两者分工明确。
核反应堆压力容器的温度评估
核反应堆的设计中，球壳热传导怎么用？
压水堆（PWR）的压力容器下半球是半球壳形状。正常运行时的肉厚方向温度分布可用球壳定常解评估。然后据此计算热应力，与压力膜应力叠加评估应力强度。
更重要的是加压热冲击（PTS）评估。紧急停堆时冷却水注入，容器内表面急速冷却，内表面产生拉应力，已有缺陷可能扩展。这类非定常问题的初始条件就是定常温度分布。
边界条件的设置
球壳热传导分析的边界条件设置有什么需要注意的？
| 边界条件类型 | 数学表达 | 物理含义 | 注意事项 |
|---|---|---|---|
| 第1类（Dirichlet） | $T(r_i) = T_0$ | 表面温度指定 | 最简单，但实际多含对流 |
| 第2类（Neumann） | $-k\frac{dT}{dr}\bigg|_{r_i} = q_0$ | 热流密度指定 | 绝热条件是$q_0 = 0$的特例 |
| 第3类（Robin） | $-k\frac{dT}{dr}\bigg|_{r_i} = h(T - T_\infty)$ | 对流换热 | $h$值影响大。需确认文献值 |
| 第4类（辐射） | $-k\frac{dT}{dr}\bigg|_{r_i} = \epsilon\sigma(T^4 - T_\infty^4)$ | 热辐射 | 非线性。高温（数百°C以上）主导 |
对流换热系数$h$的值如何确定，这在工程上似乎是大问题。
确实，这需要经验。球体外表面自然对流可用Churchill & Chu相关式（$\text{Nu} = 2 + 0.589 \text{Ra}^{1/4}/[1+(0.469/\text{Pr})^{9/16}]^{4/9}$），但结果对$h$的假定线性敏感。$h$偏差20%，侵入热量也偏差20%。必须进行灵敏度分析，评估$h$的不确定性对结果的影响。
质量保证检查表
球壳热传导分析的审核，要看哪些关键点？
Moss型液化天然气球形罐的隔热设计
Moss型液化天然气罐（由挪威Kvaerner公司1970年代开发）采用直径40米的球形铝罐，外覆珍珠岩隔热材。隔热材厚500mm能将侵入热量控制在约80kW，每天约30吨液化天然气蒸发（BOG率0.04%）。这些蒸发气体作为锅炉燃料消耗，并非浪费，但若BOG率过高会导致罐内压力升高，安全阀作动。
球壳定常热传导的软件比较
商业工具的实现
能求解球壳定常热传导的CAE软件有哪些选择？
主要的商业软件及球壳建模特点对比如下。
| 软件 | 轴对称单元 | 多层支持 | 辐射耦合 | 脚本 |
|---|---|---|---|---|
| Ansys Mechanical | PLANE55/77 | 按材料号分层 | SURF151/152 | APDL / PyAnsys |
| Abaqus | DCAX4/DCAX8 | Section分配 | *RADIATION | Python脚本 |
| COMSOL | 2D轴对称 | GUI域分割 | 表面到环境 | MATLAB LiveLink |
| Code_Aster | AXIS_FOURIER | 材料分配 | ECHANGE_PAROI | Python / comm文件 |
APDL实现示例
能展示一个用Ansys APDL求解球壳热传导的最小脚本吗？
球壳（$r_1 = 0.1$ m、$r_2 = 0.2$ m、$k = 50$ W/(m·K)）定常温度场的APDL示例。
/PREP7
ET,1,PLANE55,,,1 ! 轴对称热传导单元
MP,KXX,1,50 ! 热导率 50 W/(m·K)
! 建立球壳断面（半圆）
K,1, 0.1, 0 ! 内径点（下端）
K,2, 0.2, 0 ! 外径点（下端）
K,3, 0, 0.2 ! 外径点（上端）
K,4, 0, 0.1 ! 内径点（上端）
LARC,1,4,0,0.1 ! 内面圆弧
L,4,3 ! 上端直线
LARC,3,2,0,0.2 ! 外面圆弧
L,2,1 ! 下端直线
AL,1,2,3,4 ! 定义面
ESIZE,,20 ! 肉厚方向20分层
AMESH,ALL
/SOLU
ANTYPE,STATIC
! 边界条件
NSEL,S,LOC,X,0.1 ! 内表面节点
D,ALL,TEMP,500 ! 内表面温度 500°C
NSEL,S,LOC,X,0.2 ! 外表面节点
D,ALL,TEMP,100 ! 外表面温度 100°C
ALLSEL
SOLVE$r = 0.15$ m处的温度输出为233.3°C即通过验证。
Abaqus的设置
Abaqus怎么设置？
Abaqus的INP文件设置示例。
*HEADING
Spherical shell steady-state conduction
*NODE
! 轴对称模型节点坐标 (r, z)
1, 0.1, 0.0
...
*ELEMENT, TYPE=DCAX4, ELSET=SHELL
! 轴对称4节点扩散单元
1, 1, 2, 12, 11
...
*MATERIAL, NAME=STEEL
*CONDUCTIVITY
50.0
*SOLID SECTION, ELSET=SHELL, MATERIAL=STEEL
*STEP, NAME=STEADY
*HEAT TRANSFER, STEADY STATE
*BOUNDARY
INNER_NODES, 11, 11, 500.0
OUTER_NODES, 11, 11, 100.0
*OUTPUT, FIELD
*NODE OUTPUT
NT
*ELEMENT OUTPUT
HFL
*END STEP关键是单元类型为DCAX4（轴对称4节点扩散单元），自由度编号11代表温度，与结构单元的1-6不同。
开源软件的选择
没有商业软件许可的情况下，有免费的工具吗？
| 软件 | 许可证 | 球壳支持 | 备注 |
|---|---|---|---|
| Code_Aster | GPL | 轴对称单元 | 法国EDF开发。文档多为法文，但英文版完整 |
| FEniCS / FEniCSx | LGPL | Python API灵活定义 | 变分形式直接编程。适合研究 |
| Elmer | GPL | 2D/3D热传导 | 芬兰CSC开发。有GUI |
| 自编Python实现 | ― | 有限差分/SciPy | 教学最简洁。用解析解即时验证 |
球壳定常热传导的故障排除
温度分布与理论解不符
老师，我的球壳有限元模型与解析解差5%，哪里有问题？
球壳偏差5%很大。逐项检查。
- 单元类型设置错误：忘了轴对称（KEYOPT(3)=1），变成平面应力/应变模型。这是最常见的原因。$r^2$的面积效应没有反映，结果完全错误。
- 对称轴位置：Ansys APDL中$r$轴对应$X$轴，球断面应在YZ平面。若画在XY平面就模型没有旋转。Abaqus也是$r=X$、$z=Y$。
- 线性单元精度不足：4节点线性单元用$1/r$分布近似，特别是薄壳（$r_2/r_1 > 3$）误差大。换8节点二阶单元或增加分层数。
除了这些，还有其他常见的错误吗？
热流密度的不连续
球壳沿半径方向输出热流密度，单元边界处波动很大…
这是有限元的特性：温度连续（C⁰），但梯度（即热流）在单元边界不连续。对策三种。
- 节点平均：相邻单元的值取平均。大多数软件默认执行。
- 用二阶单元：梯度近似精度高，不连续幅度减小。
- 加密网格：足够细则不连续宽度忽略不计。
后处理中输出"Unaveraged"（未平均）热流密度，相邻单元差超5%表示网格不足。
多层模型的接触热阻
多层球壳模型中，层间接触热阻怎么处理？手算假设完全接触，实际有缝隙。
考虑层间接触热阻$R_c$ [m²·K/W]有三种方法。
- 等效薄层：在接触面插入虚拟薄层，厚度$\delta$（如0.1mm），导热率$k_c = \delta / R_c$。与手算兼容。
- 接触单元：Ansys的CONTA/TARGE单元，Abaqus的TIE + GAP CONDUCTANCE定义接触热导。
- 直接导热系数：COMSOL"薄层"功能直接输入接触导热系数$h_c = 1/R_c$ [W/(m²·K)]。
接触面粗糙度（打磨、涂脂等）使$R_c$变化范围$10^{-5}$～$10^{-2}$ m²·K/W。接触面压力越大，微凸顶端压扁，$R_c$减小。
常见的单位错误
前辈说"单位错误浪费了三天"。球壳计算特别容易出什么单位错？
| 现象 | 原因 | 对策 |
|---|---|---|
| 热流量大1000倍 | 半径用mm，$k$用W/(m·K) | 全用SI（m）。或mm系统用mW/(mm·K) |
| 温度分布呈线性 | 轴对称KEYOPT忘设 | 与解析解对比马上发现 |
| 多层界面温度断跳 | 节点不共享（网格分离） | Merge nodes或Tie constraint连接 |
| 外表面热流为零 | 对流$h$单位错误 | 确认$h$为W/(m²·K)。mm系统乘10⁻³ |
最后回到起点——解析解对比是最有效的检验。
完全同意。解析解存在的问题是验证自己模型的宝贵机会。球壳的$T(r) = 1/r$分布验证合格后，再进行复杂的3D模型，是最稳妥的做法。