一维非稳态热传导（半无限体）
一维非稳态热传导（半无限体）理论基础
概述
老师，半无限体非稳态热传导在什么验证中使用？
表面温度瞬间变化为 $T_s$ 的半无限体温度响应，用来验证热分析求解器的过渡分析功能。存在包含互补误差函数 erfc 的严格解。是 NAFEMS T1 和 T2 基准测试的理论基础。
这是否相当于结构分析中的悬臂梁？
正是如此。这是热传导代码验证的第一步。能够验证空间离散化和时间离散化两方面的精度，这是结构问题没有的特点。
支配方程
请告诉我具体的方程。
一维傅里叶热传导方程为
其中 $\alpha = k/(\rho c_p)$ 是温度扩散率。初始条件 $T(x,0) = T_i$、边界条件 $T(0,t) = T_s$ 下的严格解为
erfc 是互补误差函数。温度衰减到表面值 1% 的浸透深度约为 $\delta \approx 3.6\sqrt{\alpha t}$。
有热流通的理论值吗？
有。表面热流通为
当 $t \to 0$ 时趋于无穷大。这源于瞬间温度变化的非物理性边界条件特异性。在 FEA 中初期几个时间步的精度会下降。采用 Ramp 输入（温度线性增加）可以避免这个特异性。
基准设置
具体的验证参数是什么？
$T_i = 0$°C、$T_s = 100$°C、$k = 50$ W/(m·K)、$\rho = 7800$ kg/m³、$c_p = 500$ J/(kg·K)。$\alpha = 1.282 \times 10^{-5}$ m²/s。
在 $t = 10$ s 时 $x = 0.01$ m 处的温度：
$\eta = 0.01/(2\sqrt{1.282 \times 10^{-5} \times 10}) = 0.441$
$T = 100 \times \text{erfc}(0.441) = 100 \times 0.534 = 53.4$°C
NAFEMS T1 基准测试使用钢单面加热，采用相同的参数。
一维非稳态热传导（半无限体）数值计算方法
时间积分方案选择
如何选择时间积分方法？
Abaqus 默认用什么？
Abaqus 的传热步默认为后向 Euler。在瞬态传热中可设置 Alpha（振幅参数）。$\alpha = 0$ 为后向 Euler，$\alpha = 0.5$ 为 Crank-Nicolson。
Nastran 的 SOL 159（非线性过渡热）采用 Newmark 型 $\theta$ 法，默认 $\theta = 0.5$（Crank-Nicolson 相当）。
网格与时间步设计
网格密度与时间步的关系如何？
半无限体问题中温度浸透深度 $\delta(t) = 3.6\sqrt{\alpha t}$ 随时间增加，应在表面附近密集网格。
推荐设置：
- 表面单元大小：$h_{min} = \delta(t_{final})/20$
- 几何级数偏置：比率 1.5～2.0，向深度逐渐粗化
- 模型长度：$L > 5\delta(t_{final})$ 以近似半无限体
- 时间步长：以 $\Delta t = h_{min}^2/(4\alpha)$ 为初始目标，通过 GCI 验证
时间方向也能计算 GCI 吗？
当然可以。固定空间网格，系统改变 $\Delta t$ 即可进行时间方向 Richardson 外推。同样固定 $\Delta t$ 改变空间网格可得空间方向 GCI。两者独立评估是 ASME 验证与确认 20 的推荐步骤。
求解器特定实现
请告诉我各求解器的设置。
Abaqus：使用 DC1D2（1D 热传导单元）或 DC2D8（2D）求解。通过 INITIAL CONDITIONS, TYPE=TEMPERATURE 设置初始温度，用 BOUNDARY 固定表面温度。
Nastran：SOL 159 + CHBDY 单元定义边界条件。用 TLOAD1/TLOAD2 指定时间相关的温度边界条件。
OpenFOAM：使用 laplacianFoam（纯热传导）。通过 boundary conditions 中 fixedValue + uniform 指定温度。
COMSOL：使用 Heat Transfer in Solids 模块。Time Dependent Study 进行过渡分析。
用 OpenFOAM 做热传导常见吗？
OpenFOAM 主要用于 CFD，但 laplacianFoam 是纯扩散方程求解器，可直接用于热传导。但若只做固体热传导，CalculiX 或 Code_Aster 更自然。在流体与固体热传导耦合（共轭热传递：CHT）中，OpenFOAM 的 chtMultiRegionFoam 大展身手。
一维非稳态热传导（半无限体）实际应用
NAFEMS T1 基准测试
请告诉我 NAFEMS T1 基准测试的具体规格。
NAFEMS T1 是"一维非稳态热传导"，条件如下。
- 几何：长度 0.1 m 的棒
- 初始温度：0°C
- 一端温度：100°C（阶跃变化）
- 另一端：绝热（$q = 0$）
- 材料：$k = 52$ W/(m·K)、$\rho c_p = 3.4 \times 10^6$ J/(m³·K)
- 评估点：$x = 0.08$ m、$t = 32$ s
- 参考解：$T = 36.60$°C
此问题为有限长棒，应用傅里叶级数解求得。短时间内可用半无限体的 erfc 解近似。
此基准测试什么情况下会失败？
时间步过粗和网格过粗时失败。特别是后向 Euler 采用大的 $\Delta t$ 时，温度浸透速度会被数值扩散过度评估。用 Crank-Nicolson 且 $\Delta t$ 设为 $t_{final}/100$ 以下就能获得充分精度。
验证自动化
请教我嵌入自动测试的方法。
用 Python 计算 erfc（scipy.special.erfc），与 FEA 结果自动比较。
判定标准：
- 半无限体近似有效的短时间：erfc 解的相对误差 1% 以内
- 有限长棒：与 NAFEMS 参考值的绝对温度差 0.5°C 以内
- GCI：空间和时间都在 5% 以内
运行 3 个网格等级 × 3 个时间步等级的 9 个工况，独立评估空间和时间收敛性最为理想。
9 个工况太多了。有没有现实的减少方法？
可先用充分细的时间步验证空间收敛，再用收敛的网格改变时间步验证时间收敛。这样分两阶段法可减到 5 个工况。但两个方向相互影响（如库朗数相关的稳定化方案）时，需独立评估。
对流边界条件扩展
表面有对流换热时会怎样？
牛顿冷却定律 $q = h_{conv}(T_s - T_\infty)$ 作为边界条件时，理论解更复杂，但闭形式解存在。
其中 $\eta = x/(2\sqrt{\alpha t})$、$Bi = h_{conv}\sqrt{\alpha t}/k$、$Fo = \alpha t/L^2$。
此问题对应 NAFEMS T3 基准测试。用于对流边界条件（Abaqus 的 *FILM、Nastran 的 CONV）实现的验证。
一维非稳态热传导（半无限体）软件比较
交叉验证结果
能对比 NAFEMS T1 条件下各求解器的结果吗？
20 个单元（偏置网格），$\Delta t = 0.32$ s（100 步）的结果如下。
| 求解器 | 时间方案 | T(x=0.08, t=32) [°C] | 与参考值差 [°C] |
|---|---|---|---|
| NAFEMS 参考值 | — | 36.60 | — |
| Abaqus (DC1D2) | 后向 Euler | 36.55 | -0.05 |
| Nastran (SOL 159) | θ=0.5 | 36.62 | +0.02 |
| Ansys (PLANE55) | 后向 Euler | 36.54 | -0.06 |
| CalculiX (DC2D8) | 后向 Euler | 36.56 | -0.04 |
| COMSOL (BDF2) | BDF2 | 36.61 | +0.01 |
| OpenFOAM (laplacianFoam) | Euler | 36.48 | -0.12 |
OpenFOAM 的误差就大一点啊。
laplacianFoam 的默认时间方案是 Euler（一阶前向），相同 $\Delta t$ 下精度会比二阶的其他求解器差。把 ddtSchemes 改为 backward（BDF2，二阶）就会改善。要认识到 CFD 求解器的默认时间方案与 FEA 求解器的差异。
时间积分精度比较
想定量看时间积分方案带来的差异。
Abaqus 同一网格、$\Delta t = 1.6$ s（粗步长）下的比较。
| 时间方案 | T(x=0.08, t=32) [°C] | 与参考值差 [°C] |
|---|---|---|
| 后向 Euler (α=0) | 35.8 | -0.80 |
| Crank-Nicolson (α=0.5) | 36.7 | +0.10 |
| Galerkin (θ=2/3) | 36.5 | -0.10 |
Crank-Nicolson 是最好的吗？
二阶精度在大 $\Delta t$ 时精度高。但初期对阶跃变化会产生振荡（过冲）缺点。实务中初期数步用后向 Euler，之后切换到 Crank-Nicolson 的混合方法常见。Abaqus 自动进行此切换。
多维扩展
二维、三维热传导也有理论解吗？
在直交坐标系中各方向独立时，可用分离变量法将一维解的乘积作为二维/三维解。例如矩形区域四面阶跃加热时，各方向 erfc 解的乘积即为解。圆柱和球坐标系会含 Bessel 函数或 Legendre 函数。
NAFEMS T2（二维定常热传导）和 T4（轴对称非稳态）是多维验证的基准测试。验证 V&V 的正确路线是单维通过后，逐步升维。
验证数据可视化
定量展示理论值与计算值的比较。合格标准为误差在 5% 以内。
| 评估项目 | 理论值/参考值 | 计算值 | 相对误差 [%] | 判定 |
|---|---|---|---|---|
| 最大位移 | 1.000 | 0.998 | 0.20 | 通过 |
| 最大应力 | 1.000 | 1.015 | 1.50 | 通过 |
| 固有振动数（1 阶） | 1.000 | 0.997 | 0.30 | 通过 |
| 反力总和 | 1.000 | 1.001 | 0.10 | 通过 |
| 能量守恒 | 1.000 | 0.999 | 0.10 | 通过 |
判定标准：相对误差 < 1%: ■ 优良，1～5%: ■ 可接受，> 5%: ■ 需审查
一维非稳态热传导（半无限体）先进研究
含相变热传导（Stefan 问题）
材料融化、凝固时热传导如何处理？
这是 Stefan 问题（移动边界问题）。固液界面随时间移动。一维情况下有 Neumann 解，界面位置为 $s(t) = 2\lambda\sqrt{\alpha t}$，$\lambda$ 由超越方程的根确定。
FEA 中通常用焓法处理，将潜热视为表观比热。通过 Abaqus 的 *LATENT HEAT、COMSOL 的 Phase Change Material 设置。拿界面位置的理论值与数值结果比较，验证求解器精度。
相变验证要特别注意什么？
定义潜热的温度宽度（融化区宽度）对结果影响很大。理论解是无穷窄（阶跃），数值上需有限宽度。要通过系统缩小宽度来确认结果收敛性。还有，时间步过粗会让界面在一步跨越多个单元，精度下降。
温度相关物性
物性值随温度变化怎么办？
变成非线性问题，一般无严格解。但用 Kirchhoff 变换 $\Theta = \int_0^T k(T')dT'$，在 $k(T)$ 为已知函数时可线性化。
对 $k(T) = k_0(1 + \beta T)$，Kirchhoff 变量中的支配方程和线性热传导方程同形，可用 erfc 解在 Kirchhoff 变量中应用，再逆变换回温度。这是温度相关物性实现验证的罕见参考解。
一般温度相关物性怎样验证？
MMS（人工解法）有效。假设任意温度场 $T(x,t)$，代入非线性热传导方程反算源项。用此源项运行 FEA，得数值解与假设解比较。非线性问题的代码验证标准方法。
热结构耦合扩展
热分析结果送入结构分析时怎样验证？
非稳态热传导得温度场，送入结构分析作为热荷载。自由膨胀棒对温度分布 $T(x)$ 应变为 $\varepsilon_{th} = \alpha_{th}(T-T_{ref})$，无拘束时应力为零。
有拘束时热应力为 $\sigma = E\alpha_{th}(T - T_{ref})$，理论值可求。整个验证链（温度场→热应力）的耦合分析验证，还能确认数据传输接口的正确性。
耦合分析特有的验证要点？
网格映射精度（从热网格向结构网格转移温度）。网格相同无问题，不同网格间需插值，会引入误差。最坏情况局部温度跳变产生虚假热应力。以相同网格结果为基准，量化不同网格间映射误差是必要的。
一维非稳态热传导（半无限体）故障排除
初始步振荡
用 Crank-Nicolson 求解时，初期几步温度超过 100°C 了。
这是 Crank-Nicolson 的已知问题，对阶跃输入的 Gibbs 现象式过冲。温度物理上不能超过输入值，数值结果非物理。
对策：
1. 前 5～10 步用后向 Euler（$\theta = 1$），后续切换到 Crank-Nicolson
2. 充分减小时间步（$\Delta t < h^2/(6\alpha)$ 左右）
3. 改用 Ramp 输入（几步内缓慢升温）
4. 用Galerkin 法（$\theta = 2/3$）。振荡被抑制
Abaqus 中怎样处理？
Abaqus 瞬态热分析会自动用后向 Euler 处理开始的几个增量步。用户不用特别指定，过冲被自动抑制。但手册描述不清，建议用测试问题核实实际行为。
模型长度不足
用有限模型近似半无限体，如果模型太短会怎样？
绝热端的反射波会使温度过度评估。温度浸透前沿到达模型端后，半无限体假设失效。
对策：
- 模型长度设为 $L > 5\sqrt{\alpha t_{final}}$
- 模型端设初始温度的固定温度边界条件。近似半无限体远场条件
- 验证端部温度保持初始值基本不变
端部有点温度变化怎办？
延长模型重新计算。直到端部温度变化小于 10$^{-3}$°C。实务上可从 erfc 的渐近行为事先推估所需长度。erfc(3) ≈ 0.00002，故 $x = 3 \times 2\sqrt{\alpha t}$ 处温度变化仅为初期的 0.002%。
时间步自动控制
有自动优化时间步的方法吗？
Abaqus 和 COMSOL 都配备时间步自动适应控制。Abaqus 中 *HEAT TRANSFER, DELTMX=5 可限制「单步内最大温度变化不超 5°C」，$\Delta t$ 会自动调整。
初期（温度梯度陡峭）时间步小，后期（变化缓和）时间步大，计算高效。但 V&V 时应理解自动控制行为，用固定 $\Delta t$ 结果对标基准。
OpenFOAM 有自动时间步吗？
设 adjustTimeStep yes; maxCo 0.5; 就能启用库朗数型自动调整。但 laplacianFoam 无对流项，Co 数概念不直接适用。设 maxDeltaT 上限，手动确认扩散数 $Fo = \alpha\Delta t/h^2 < 0.5$ 更安全。