圆柱坐标系的热传导
圆柱坐标系热传导的理论基础
圆柱坐标系的热传导方程
老师，配管和轴的热解析在笛卡尔坐标下不能很好处理吗？
圆柱形状在笛卡尔坐标下表达浪费很大，网格复杂化。使用圆柱坐标 $(r, \phi, z)$ 可利用对称性降低维数。一维径向定常热传导方程为
与笛卡尔坐标的傅里叶定律相比，多了 $1/r$ 项。
体积元素与 $r$ 成正比，这项反映了面积的变化。这也是临界绝热半径出现的物理根源。
解析解（内部发热无）
当 $\dot{q}_v = 0$ 时的解为
温度分布为对数函数，不同于平板的线性分布。热流量为
$\ln$ 出现是因为面积随 $r$ 变化的缘故。
正确。截面积 $A = 2\pi r L$ 与 $r$ 成正比，所以热流密度 $q'' = q/A$ 与 $r$ 成反比。内表面的热流密度最大。
内部发热情况
有均匀内部发热 $\dot{q}_v$ 时的解为
在中心 $r=0$ 处达到最高温度 $T_{\max} = T_s + \dot{q}_v R^2/(4k)$。
电线的焦耳发热就是这种情况。
正是如此。AWG18铜线（直径1.02mm）通10A电流产生约10 W/m的热。$\dot{q}_v \approx 1.2 \times 10^7$ W/m$^3$，可用该式估算中心温度升高。
圆柱坐标拉普拉斯算子的来源
圆柱坐标的热传导方程为∂²T/∂r² + (1/r)∂T/∂r + (1/r²)∂²T/∂θ² + ∂²T/∂z² + q̇/λ = 0。该拉普拉斯算子中的(1/r)∂T/∂r项源自坐标变换的雅可比行列式。Lamé在1833年用弹性论整备了圆柱坐标拉普拉斯算子，其后20年才被用于Fourier热方程。
圆柱坐标系热传导的数值计算方法
轴对称FEM模型
用FEM求解圆柱热传导时，需要建立3D模型吗？
若沿周向对称，2D轴对称模型就足够了。计算成本可降低到1/100以下。
| 模型类型 | 自由度数 | 计算时间 | 精度 |
|---|---|---|---|
| 3D全体 | ～100万 | 分钟级 | 基准 |
| 3D 1/4对称 | ～25万 | 秒级 | 等效 |
| 2D轴对称 | ～1000 | 瞬间 | 等效 |
2D轴对称足够的情况下还用3D就太浪费了。
利用对称性是计算效率的基础。在Ansys Mechanical中使用PLANE55（KEYOPT(3)=1），在Abaqus中使用DCAX4（轴对称4节点单元）。
有限差分法离散化
圆柱坐标的中心差分需要小心处理。
在 $r = 0$（中心轴）处 $1/r$ 会出现奇点，需应用洛必达法则：
中心轴的奇点需要特殊处理。
FEM也是如此。轴对称单元会自动正确处理中心轴，但用3D模型在轴上放置楔形单元时需要格外注意。
多层圆柱求解
多层圆柱（管+保温材料+覆盖层等）各层可用解析解直接连接。
包含层间接触热阻时，分母中加入 $1/(2\pi r_i h_{c,i} L)$。
和电路中串联电阻的求和一样。
正是热阻网络的思想。即使5层以上的多层结构也能手算，这是圆柱坐标的一大优势。
对数温度分布的求解步骤
无限长空心圆柱（内半径ri、外半径ro）的定常温度分布为T(r)=Ti + (Ti−To)·ln(r/ro)/ln(ri/ro)，呈对数型。该求解过程为变量分离→常微分方程积分→2个边界条件确定常数，共3个步骤完成。日本大学入试中有类似题目，国内热工学教科书必在章首例题中讲解该问题。
圆柱坐标系热传导的工程应用
配管放热计算
请讲一下圆柱热传导在工程中的典型应用。
最常见的是工业装置配管的热损失计算。举蒸汽管（外径114.3mm、SUS304、保温层50mm）的例子。
| 层 | $r$ [mm] | $k$ [W/(m K)] | 热阻 [m K/W] |
|---|---|---|---|
| 管壁 | 48.6→57.15 | 16.3 | 0.00160 |
| 保温材料 | 57.15→107.15 | 0.05 | 2.016 |
| 内表面对流 | — | — | 0.00033 |
| 外表面对流 | — | — | 0.0149 |
全热阻 $R_{\text{total}} = 2.033$ m K/W。保温材料占全热阻的99%。
管壁的热阻几乎可以忽略。
金属管的热阻比保温材料小1000倍，实际上没有温降。因此管壁用铜还是不锈钢对放热量影响微乎其微。
对数平均半径
圆柱单位长度热阻 $R = \ln(r_2/r_1)/(2\pi k)$ 可用等效平板 $R = t/(k \cdot A_{lm})$ 表示。$A_{lm}$ 是对数平均面积
用对数平均面积就能套用平板公式。
当 $r_2/r_1 < 2$ 时，用算术平均 $(r_1+r_2)/2$ 替代，误差在4%以内。简易计算常采用算术平均。
网格生成的注意事项
圆柱轴对称网格的要点列如下。
- 径向采用偏置网格（内表面部分细化）
- 薄壁管径向至少3个单元
- 多层结构中，层间节点共享或使用Tied Contact
偏置网格是因为内表面热流密度更大吧。
正是。热流密度 $q'' \propto 1/r$，内表面温度梯度陡峭，需要细密网格。
核反应堆燃料棒的温度分布
核反应堆中UO₂燃料棒（直径8～10 mm、λ≈3 W/m·K）的内部发热用圆柱热传导分析，中心温度可比表面高出最多1800℃。1979年的三里岛事故中燃料棒中心超过熔点（2865℃）的部分，事故后温度复现解析被整理在NRC报告书（NUREG-0772）中。
圆柱坐标系热传导的软件对比
商用工具中的实现
各工具如何设置圆柱热传导分析？
比较轴对称模型的设置方法。
| 工具 | 单元类型 | 设置方法 |
|---|---|---|
| Ansys Mechanical | PLANE55 (KEYOPT(3)=1) | 将分析类型设置为2D轴对称 |
| Abaqus | DCAX4, DCAX8 | 在Part定义时选择Axisymmetric |
| COMSOL | 2D Axisymmetric | 模型向导中选择坐标系 |
| STAR-CCM+ | 2D Axisymmetric mesh | 在网格模型中指定轴对称 |
APDL实现示例
给出蒸汽配管放热计算的例子。
```
/PREP7
ET,1,PLANE55,,,1 ! 轴对称
! 管壁 (SUS304)
MP,KXX,1,16.3
! 保温材料
MP,KXX,2,0.05
CYL4,0,0,48.6,0,57.15,90 ! 管壁
CYL4,0,0,57.15,0,107.15,90 ! 保温材料
AGLUE,ALL
ESIZE,1
AMESH,ALL
/SOL
SFL,内面LINE,CONV,500,180 ! 蒸汽侧 h=500, T=180℃
SFL,外面LINE,CONV,10,25 ! 外侧 h=10, T=25℃
SOLVE
```
内表面 $h$ 是500，外表面是10，数值差异很大。
蒸汽的强制对流 $h = 100$～$1000$ W/(m$^2$ K)，外表面自然对流 $h = 5$～$15$ 左右。所以温降主要发生在保温层。
Abaqus设置
Abaqus用*FILM条件定义对流。
```
*STEP
*HEAT TRANSFER, STEADY STATE
*FILM
inner_surface, F2, 180., 500.
outer_surface, F2, 25., 10.
*END STEP
```
记法虽然不同，但物理设置是一样的。
正是。求解器的记法各异，但离散化依据傅里叶定律是共通的。结果也应该相同。
Abaqus轴对称单元CAX4的威力
Abaqus（Dassault Systèmes）的轴对称实心单元CAX4可用2D网格高精度求解圆柱形状热传导。相比完全3D模型，该方法可将节点数削减99%以上。石油化工装置厚壁配管（壁厚超100mm）的热应力分析采用此方法，相关案例在2019年Abaqus用户会议上有所介绍。
圆柱坐标系热传导的先端研究
各向异性圆柱热传导
CFRP圆柱的 $k$ 因方向而异对吧。
圆柱坐标各向异性热传导方程为
CFRP的缠绕方式（螺旋、环向、0度）会造成 $k_r$、$k_\phi$、$k_z$ 差异很大。
| 方向 | $k$ [W/(m K)] |
|---|---|
| 纤维方向 | 5～10 |
| 垂直方向 | 0.5～1 |
| 厚度方向 | 0.3～0.7 |
各向异性达10倍以上。
CFRP压力容器（储氢罐、火箭马达壳体）的热分析必须考虑各向异性。在Ansys中把材料坐标系设定为圆柱坐标，分别定义KXX/KYY/KZZ。
包含接触问题的多层圆柱
烧结和压入的圆柱中，接触压力随温度变化，接触热阻也会变动。需要结构-热的耦合分析。
温度升高导致热膨胀改变接触压力，接触热阻随之改变，又影响温度分布，形成恶性循环。
正是这样。Ansys Mechanical和Abaqus都支持Sequential Coupled Thermal-Stress分析来处理此耦合。还有更高级的Fully Coupled分析可一次求解两个物理过程。
旋转体热传导
旋转机械（汽轮机转子、电动机转子）中，旋转会产生离心应力，同时还有摩擦热和风损。需要将热传导和流体对流在旋转坐标系中耦合（CHT分析）。
汽轮机转子是典型的圆柱坐标问题。
在Ansys CFX或STAR-CCM+中设置Rotating Reference Frame可进行CHT分析。转子温度分布直接影响蠕变寿命，需要高精度分析。
旋转圆柱中的热传导与离心力耦合
高速旋转的汽轮机转子中，离心力产生的弹性变形与热膨胀耦合，求解温度分布时需同时求解热传导和弹性方程。GE Aviation在1980年代为T700发动机汽轮机转子建立了轴对称热-结构耦合分析模型，该手法成为今日Ansys Coupled Field分析的原型。
圆柱坐标系热传导的故障排查
常见故障与对策
请讲圆柱热传导分析的注意点。
整理一下常见的故障。
1. 2D与3D结果不一致
原因：2D模型忘记设置轴对称。默认采用平面应力/平面应变，其中没有 $r$ 的效果。
对策：在Ansys PLANE55中设KEYOPT(3)=1，Abaqus中选择轴对称单元，COMSOL在建模时选"2D Axisymmetric"。
2. 多层结构温度分布不连续
各层接合处有问题？
原因：层间节点未共享或Tied Contact/Bonded Contact设置遗漏。
对策：在公共面处合并节点，或正确设置接触条件。如要引入层间接触热阻，则设定Gap Conductance。
3. 内表面对流系数估算错误
管内 $h$ 值因流速、流体物性、管径而变化很大。
| 流体条件 | $h$ [W/(m$^2$ K)] |
|---|---|
| 空气·自然对流 | 5～25 |
| 空气·强制对流 | 25～250 |
| 水·强制对流 | 500～10,000 |
| 蒸汽凝结 | 5,000～50,000 |
| 沸腾 | 2,500～75,000 |
水和空气相差100倍以上。
用Dittus-Boelter公式或Gnielinski公式计算是基本做法。$\text{Nu} = 0.023 \text{Re}^{0.8} \text{Pr}^{0.4}$（Dittus-Boelter）。Re数算错会导致 $h$ 数值偏差一个数量级。
4. 单位制混淆
SI制（m）和mm制时 $k$ 的单位不同。$k = 16.3$ W/(m K) 换成mm制是 $k = 0.0163$ W/(mm K)。Ansys Mechanical的Workbench环保境默认mm制，需要特别注意。
单位制真的很容易出错。
养成用理论解检算的习惯。$q = 2\pi k L \Delta T / \ln(r_2/r_1)$ 手工计算一下，和结果对比就能发现单位制错误。
内外径单位错误导致的重大事故
配管圆柱热传导分析中，ri和ro应输入mm但误输为m，导致壁厚计算放大1000倍，这类人为错误在现场屡见不鲜。鉴于1999年NASA火星气候轨道器坠毁事故（磅力与牛顿单位混用），多家求解器厂商在输入单位验证处添加了明确的单位确认对话框。
关联主题
细的内容
错误