2. 热弹性力学与稳态热传导基础

在热结构协同拓扑优化中,结构不仅要承受机械载荷,还要承受温度场引起的热载荷。理解热应力、热应变、热弹性本构方程以及稳态热传导方程,是建立热-结构耦合优化模型的前提。本章从物理量定义出发,逐步建立热弹性问题的数学描述,并说明这些方程在拓扑优化中的具体形式。

2.1 热应变与热膨胀系数

当温度变化时,材料内部质点间的平衡距离发生改变,产生热应变。对于各向同性材料,自由热膨胀不产生应力,热应变张量定义为:

ε_th = α · ΔT · I

其中 α 为线膨胀系数(单位 1/K),ΔT = T - T_ref 为相对参考温度的温升,I 为二阶单位张量。α 的物理含义是:单位温升引起的单位长度相对伸长量。例如,铝合金 α ≈ 23×10⁻⁶ /K,钢 α ≈ 12×10⁻⁶ /K,碳纤维复合材料沿纤维方向可低至 1×10⁻⁶ /K 量级。α 的数值差异直接决定了多材料结构中的热变形不协调程度。

需要特别强调符号约定:本文统一取温升为正(ΔT > 0),此时热应变为正(膨胀);若结构受到外部约束无法自由膨胀,则产生压应力。这一约定在后续所有公式和图示中保持一致。

2.2 热弹性本构方程

热弹性问题的本构方程将总应变分解为弹性应变和热应变两部分:

σ = D : (ε - ε_th) = D : (ε - α·ΔT·I)

其中 D 为弹性刚度张量,σ 为应力张量,ε 为总应变张量。该式表明:热应力并非由温度直接产生,而是由“总应变与热应变之差”即弹性应变驱动。若结构完全自由且温度均匀,ε = ε_th,则 σ = 0;若存在约束或温度梯度,弹性应变不为零,热应力随之产生。

在有限元离散形式下,热弹性本构关系写为:

σ = D (B u - α ΔT I)

其中 B 为应变-位移矩阵,u 为节点位移向量。对应的热载荷等效节点力为:

F_th = ∫ Bᵀ D α ΔT I dΩ

这一项将温度场转化为结构分析中的等效外载,是热-结构耦合的核心桥梁。

2.3 稳态热传导方程与边界条件

稳态热传导问题中,温度场不随时间变化,能量守恒退化为热流散度为零:

∇·(k ∇T) + Q = 0

其中 k 为导热系数(W/(m·K)),Q 为体积热源密度(W/m³)。该方程说明:进入微元体的净热流加上内部热源,必须为零。

热边界条件分为三类,它们在拓扑优化中的处理方式不同:

  • 第一类(Dirichlet)——给定温度: T = T̄,在边界 Γ_T 上直接约束节点温度。拓扑优化中通常固定在设计域边界,不随材料分布变化。
  • 第二类(Neumann)——给定热流: -k ∂T/∂n = q̄,在边界 Γ_q 上施加热流密度。热流方向沿外法线,正号表示流入。
  • 第三类(Robin)——对流换热: -k ∂T/∂n = h(T - T_∞),其中 h 为对流换热系数(W/(m²·K)),T_∞ 为环境温度。该条件将边界温度与环境温度耦合,是散热结构优化中最常见的边界。

辐射边界在高温场景不可忽略,其非线性形式为:

-k ∂T/∂n = ε_r σ_s (T⁴ - T_∞⁴)

其中 ε_r 为发射率,σ_s 为 Stefan-Boltzmann 常数。由于 T⁴ 非线性,辐射边界通常需要迭代求解。

2.4 热-结构耦合在拓扑优化中的数学描述

在热结构协同拓扑优化中,设计变量 ρ_e(单元相对密度)同时影响导热系数和弹性模量。常用的插值模型为:

k(ρ_e) = k_min + ρ_e^p (k_0 - k_min)
E(ρ_e) = E_min + ρ_e^p (E_0 - E_min)

其中 p 为惩罚因子(通常取 3),k_min 和 E_min 为空洞材料的小量。温度场和位移场通过以下耦合系统求解:

K_T(ρ) T = F_T
K_u(ρ) u = F_u + F_th(T, ρ)

其中 K_T 为热传导矩阵,K_u 为刚度矩阵,F_th 为热载荷。注意 F_th 依赖于温度场 T,而 T 又通过 k(ρ) 依赖于设计变量,因此这是一个双向耦合问题。在优化迭代中,通常先求解温度场,再计算热载荷并求解位移场,最后计算目标函数和灵敏度。

热-结构耦合拓扑优化求解链路 设计变量 ρₑ 单元相对密度 材料插值 k(ρ)=k_min+ρᵖ(k₀-k_min) E(ρ)=E_min+ρᵖ(E₀-E_min) 稳态温度场 T ∇·(k∇T)+Q=0 热载荷 F_th ∫BᵀDαΔT I dΩ 位移场 u K_u(ρ)u = F_u + F_th 热应力 σ σ = D(Bu - αΔT I) 目标函数 柔度 / 热应力指标 灵敏度分析 ∂c/∂ρₑ, ∂σ/∂ρₑ 设计变量更新 ρₑ ← ρₑ + Δρₑ 迭代循环直至收敛 热边界条件 T̄ / q̄ / h(T-T∞) / 辐射 作用于温度场求解

2.5 关键参数与量级参考

为便于建立数量级直觉,下表给出常见材料在室温附近的热弹性参数参考值。这些数值用于教学示例和初步估算,实际优化中应以具体材料数据为准。

材料 导热系数 k (W/(m·K)) 线膨胀系数 α (10⁻⁶/K) 弹性模量 E (GPa) 典型热应力来源
铝合金 6061 167 23.6 69 约束膨胀、温度梯度
结构钢 Q345 50 12.0 206 焊接热循环、约束
碳纤维/环氧(单向) 5~10(横向) 1~30(方向相关) 70~140 铺层方向热膨胀不匹配
氧化铝陶瓷 30 8.0 300 热冲击、脆性断裂
教学示例:取一根两端固定的铝合金杆,长度 L = 200 mm,温升 ΔT = 50 K,α = 23.6×10⁻⁶ /K,E = 69 GPa。自由热伸长量 ΔL = α·ΔT·L = 0.236 mm。由于两端固定,该伸长被完全约束,产生压应力 σ = -E·α·ΔT ≈ -81.4 MPa。若杆件同时承受 10 kN 拉伸机械载荷,截面 100 mm²,则机械应力为 100 MPa,叠加后实际应力约为 18.6 MPa(拉)。这说明热应力与机械应力可能相互抵消或叠加,拓扑优化中必须同时考虑。
注意:当温度梯度较大或材料导热系数较低时,稳态假设可能不成立,需要采用瞬态热传导分析。此外,辐射边界的 T⁴ 非线性会使问题变为非线性耦合,灵敏度推导需额外处理。在拓扑优化中,若忽略辐射,可能导致高温区域散热能力被高估,优化结果偏于保守或失真。

综上,热弹性力学与稳态热传导为热结构协同拓扑优化提供了完整的数学基础:热应变通过 α 和 ΔT 定义,热应力通过本构方程与弹性应变关联,温度场通过稳态热传导方程和边界条件确定,最终以等效热载荷的形式进入结构平衡方程。设计变量通过插值模型同时影响导热和刚度,形成双向耦合。理解这一链路,是后续建立优化模型、推导灵敏度和解释优化结果的前提。