4. 热力耦合控制方程与数值方法

金属切削仿真之所以比一般结构分析复杂，根本原因在于温度场与变形场互相“喂养”：塑性变形和摩擦生热改变温度，温度又通过热软化、热膨胀改变流动应力和应变分布。要得到可信的切削力、切屑形态和刀具温度，必须把这两类物理场放进同一套控制方程中求解。本节围绕热弹塑性本构、能量守恒与动量守恒的耦合形式、有限元离散以及显式/隐式求解策略展开，重点说明“方程怎么建、离散怎么做、求解器怎么选”。

一、热弹塑性本构：把温度写进应力—应变关系

切削区材料经历大应变、高应变率和显著温升，本构模型必须同时反映硬化、率相关和热软化。工程上常用 Johnson-Cook 形式：

σ = [A + B·εⁿ] · [1 + C·ln(ε̇/ε̇₀)] · [1 − ((T − T₀)/(T_m − T₀))^m]

其中 A 为初始屈服应力，B、n 描述应变硬化，C 描述应变率敏感度，m 描述热软化。以常见切削钢的示例参数（仅作教学假设）为例：A=553 MPa、B=600 MPa、n=0.234、C=0.0134、m=1.0，熔点 T_m=1733 K，参考温度 T₀=293 K。当切削区温度从 293 K 升到 900 K 时，热软化因子从 1 降到约 0.58，意味着同样应变下流动应力下降四成以上——这正是切削力随速度升高而不再线性增长的原因。

热弹塑性本构还需给出热膨胀与弹性模量随温度的变化。总应变分解为弹性、塑性、热应变三部分：ε = εᵉ + εᵖ + εᵗ，其中 εᵗ = α(T)(T − T_ref)。塑性部分采用关联流动法则和等向硬化，屈服函数 f = σ_eq − σ_y(εᵖ, ε̇ᵖ, T) ≤ 0。温度进入屈服面，使本构成为“热—力”双向耦合的入口。

二、动量守恒与能量守恒的耦合形式

动量守恒（准静态忽略惯性时）写作：∇·σ + f = 0；能量守恒写作：ρc·∂T/∂t = ∇·(k∇T) + Q_plastic + Q_friction。两式的耦合项非常明确：

• 力→热：塑性功转化为热，Q_plastic = β·σ_eq·ε̇ᵖ，β 为 Taylor-Quinney 系数，通常取 0.85～0.95；刀屑摩擦功 Q_friction = τ_f·v_slip。
• 热→力：温度通过本构的 σ_y(T) 和热应变 εᵗ 反馈到动量方程。

这种双向耦合意味着不能先算完力场再算温度场，而要在每个时间步内交替或同时更新。若忽略热→力反馈，只做“热源单向加载”，得到的切削力会偏高，切屑形态也会偏“硬”。

三、有限元离散：空间与时间两条线

空间离散把连续体划分为单元，位移场 u = N·u_e，温度场 T = N·T_e。代入动量方程得到离散平衡方程：M·ü + F_int(u,T) = F_ext；代入能量方程得到热平衡方程：C·Ṫ + K_t·T = Q(u,ε̇ᵖ)。其中 F_int 包含本构积分得到的应力，Q 包含塑性功和摩擦热。切削仿真中，切屑与工件之间的大滑移接触通常用罚函数或拉格朗日乘子处理，接触面同时传递法向力和摩擦热，是耦合最活跃的区域。

时间离散决定求解器类型。显式方法用中心差分，不需要组装整体刚度矩阵，但稳定步长受最小单元尺寸和波速限制。以典型切削钢为例，E≈210 GPa、ρ≈7800 kg/m³，波速 c≈5190 m/s；若最小单元尺寸为 5 μm，则 Δt ≤ 5×10⁻⁶ / 5190 ≈ 9.6×10⁻¹⁰ s，即约 1 ns。要模拟 1 ms 的切削过程，需要约百万个时间步，计算量巨大但每步极快。隐式方法用 Newton-Raphson 迭代，步长可放大到 10⁻⁵～10⁻⁴ s，但每步需要更新切线刚度并处理收敛问题，在大变形、接触状态突变时容易不收敛。

四、显式与隐式求解策略的选择

实际工程中很少“二选一”，更常见的是分阶段混合：

• 切屑形成阶段：用显式求解，因为涉及材料分离、大变形、自接触和绝热剪切，隐式几乎无法稳定收敛。质量缩放可适度放大稳定步长，但需控制动能与内能之比小于 5%～10%，否则惯性效应会污染切削力。
• 稳态温度场与刀具应力：把显式得到的切削力、热流密度作为边界条件，切换到隐式求解稳态热传导和刀具热应力，效率更高。
• 耦合方式：强耦合（每步同时更新温度与位移）精度高但代价大；弱耦合（交替迭代，每若干步交换一次场变量）适合热惯性较大的工况。切削区温度梯度极陡，建议采用强耦合或至少每步交换。

五、一个可复现的示例设置

以二维正交切削为例（教学假设参数）：工件材料为 45 钢，刀具为硬质合金，切削速度 120 m/min，进给 0.15 mm/r，切削深度 2 mm，刀具前角 6°、后角 6°。显式分析步总时长 2 ms，初始时间增量 1×10⁻⁹ s，质量缩放因子 100，单元类型 CPE4RT（四节点热力耦合平面应变单元）。接触摩擦系数 0.3，摩擦生热分配系数 0.5。输出变量包括：切削力 F_c、进给力 F_f、切削区最高温度、切屑厚度比。若 F_c 在 0.8～1.2 ms 内稳定在 620～680 N 区间，且最高温度稳定在 850～950 K，可认为该阶段结果进入可接受范围；若温度持续上升不收敛，应检查热边界条件或摩擦热分配是否合理。

热力耦合仿真的核心不是“把两个场都打开”，而是理解耦合项如何进入方程、离散后如何传递、求解器如何取舍。本构决定温度对强度的反馈强度，能量方程决定热源分布，离散方式决定计算代价，显式/隐式选择决定你能算多长、算多准。把这四层关系理清，后续的切屑形态、刀具磨损和残余应力分析才有可靠基础。