球壳稳态热传导
理论与物理
球壳与圆柱的区别
老师，球壳的稳态热传导和圆柱坐标的情况有什么不同？两者都是“圆形”，看起来话题会很相似。
最大的区别在于传热面积的增长方式。圆柱的面积与半径成正比（$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)$ |
原来如此，面积的“幂次”越高，温度下降就越平缓啊。球壳在实际工作中用在什么场合呢？
主要有三个典型应用。
- LNG球形储罐（直径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$分布类比，形式变成了$\ln r \to 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，所以球壳在内表面附近的温度下降更急剧。这也意味着内表面附近的热流密度更高。
球壳的热阻
通过球壳的热流量怎么计算？能像圆柱一样使用热阻的概念吗？
可以。将$A(r) = 4\pi r^2$代入傅里叶定律 $q = -kA(r)\frac{dT}{dr}$ 并积分，得到
整理成 $Q = \Delta T / R$ 的形式，则球壳的传导热阻为
与圆柱不同，球壳的特点是没有出现长度$L$。因为球是全方向封闭的形状，所以不需要相当于管长的参数。
如果内表面和外表面有对流，是不是只要加上对流热阻就可以了？
是的。球面的对流热阻是 $R_{\text{conv}} = 1/(4\pi r^2 h)$，所以从内侧流体温度$T_{\infty,1}$到外侧流体温度$T_{\infty,2}$有
这是对流+传导+对流的串联回路。和电路欧姆定律的结构完全一样。
多层球壳模型
像LNG储罐那样夹有隔热材料的情况，是要按多层计算吧？
对于$n$层球壳，各层热导率为$k_i$，半径为$r_0, r_1, \ldots, r_n$时，总传导热阻为
包含对流在内的总传热系数$U$，若取内表面积$A_1 = 4\pi r_0^2$为基准面积，则
使用这个公式，就可以参数化地计算改变隔热材料厚度时的热损失了，对吧。
是的。以LNG球形储罐（直径40m）为例，钢壳30mm（$k = 50$ W/(m·K)）+珍珠岩隔热材料500mm（$k = 0.04$ W/(m·K)）时，隔热材料的热阻约占钢壳的约2,500倍。热阻法的优势在于能够定量地把握哪个层是主导因素。
存在内部发热的情况
像反应堆燃料球那样内部发热的情况会怎样呢？
对于具有均匀体积发热率$\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$，所以球的中心温升较小。因为表面积大，更容易散热。
那么高温气冷堆的TRISO燃料颗粒（直径约1mm）就可以用这个公式估算中心温度了。
是的。对于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亿年，但由于不知道放射性衰变产生的内部发热，导致严重低估）。即使在现代，球壳热传导仍被用作地球地幔对流模型的初始近似。
球坐标拉普拉斯算子的推导
- 拉普拉斯算子$\nabla^2 T$（球坐标）：将直角坐标的$\frac{\partial^2 T}{\partial x^2}+\frac{\partial^2 T}{\partial y^2}+\frac{\partial^2 T}{\partial z^2}$转换为球坐标$(r,\theta,\phi)$，得到 $$\nabla^2 T = \frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial T}{\partial r}\right) + \frac{1}{r^2\sin\theta}\frac{\partial}{\partial\theta}\left(\sin\theta\frac{\partial T}{\partial\theta}\right) + \frac{1}{r^2\sin^2\theta}\frac{\partial^2 T}{\partial\phi^2}$$ 若球对称（$T$仅为$r$的函数），则$\theta$项和$\phi$项消失，剩下$\frac{1}{r^2}\frac{d}{dr}\left(r^2\frac{dT}{dr}\right) = 0$。
- $r^2$的物理意义：源于半径$r$处球面的面积$4\pi r^2$。当恒定热流量$Q$通过时，热流密度$q'' = Q/(4\pi r^2)$与$r^2$成反比减小。
假设条件与适用范围
- 稳态：无时间变化（$\partial T/\partial t = 0$）。经过足够长时间、蓄热项可忽略后的温度场。
- 球对称：温度仅为$r$的函数（与$\theta$, $\phi$无关）。存在局部热源或非均匀对流时不适用。
- 各向同性·均质：各层内热导率$k$恒定。温度依赖性较大时（钢的$k$在0〜800°C下从50→25 W/(m·K)变化）需要进行非线性分析。
- 傅里叶定律：$q = -k \nabla T$。极低温（超流氦）或飞秒激光加热时需要非傅里叶模型。
量纲分析与主要参数
| 变量 | SI单位 | 典型值·注意事项 |
|---|---|---|
| 温度 $T$ | K 或 °C | 包含辐射时必须使用绝对温度[K] |
| 热导率 $k$ | W/(m·K) | 钢≈50、隔热材料≈0.03〜0.05、UO₂≈3 |
| 内径 $r_1$ | m | 球壳要求$r_1 > 0$（实心球另有解） |
| 球壳热阻 $R$ | K/W | $(1/r_1 - 1/r_2)/(4\pi k)$ ― 不含长度$L$ |
| 体积发热率 $\dot{q}$ | W/m³ | 核燃料≈$10^8$、焦耳发热≈$10^5$〜$10^7$ |
数值解法与实现
轴对称FEM模型
用FEM求解球壳热传导时，必须建立3D球体模型吗？
如果球对称，2D轴对称模型就足够了。相当于绘制一个半圆截面，并绕对称轴旋转而成的模型。计算成本是3D模型的1/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。
有限差分法的离散化
如果不想用FEM，想自己用有限差分法实现呢？我想用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$网格上离散化，中心差分就能给出精确解。这个变换是球壳热传导特有的技巧。
网格设计指南
用FEM划分网格时，根据什么标准决定单元尺寸呢？
对于球壳稳态热传导，可以参考以下指南。
| 项目 | 推荐值 | 理由 |
|---|---|---|
| 壁厚方向分割数 | 10层以上 | 为了准确再现$1/r$分布 |
| 内表面侧偏置比 | 1:3〜1:5 | 内表面温度梯度陡峭 |
| 单元类型 | 推荐二次单元 | $1/r$分布用线性单元近似较粗糙 |
| 多层界面 | 节点对齐 | 界面处温度梯度不连续 |
| 长宽比 | 5以下 | 过度扁平的单元会导致精度下降 |
相关主题
なった
詳しく
報告