变形固体中的热传递

在之前的博客文章中，我们介绍了关于静止固体共轭传热的应用。静止固体简化了将要求解的热方程，并且通常是温度场的精确近似值。今天，我们将向大家介绍传热和固体力学耦合时材料热弹性效应的相关物理场。

材料和空间坐标系

在介绍物理场之前，我们先简要回顾一下 COMSOL Multiphysics 中使用的坐标系。当几何非线性时，固体力学接口会区分材料坐标系和空间坐标系。材料坐标系以初始状态 \mathbf{X}= (X, Y, Z) 的坐标表示物理量，而空间坐标系以当前状态 \mathbf{x}=(x, y, z) 的坐标表示物理量。

下图为两个承受压应变的正方形示例。该正方形长 10cm，其左下角最初位于 (X, Y) = (1~\textrm{cm}, 1~\textrm{cm}) ，之后其左右两侧受到边界载荷的压缩而发生变形。这种变形会使正方形上几乎所有点的位置发生改变。例如，左下角移至新的位置 (x, y) = (1.54~\textrm{cm},0.82~\textrm{cm}) 。

以材料坐标表示的变形正方形，左侧为初始状态，右侧为最终状态。

以空间坐标表示的变形正方形，左侧为初始状态，右侧为最终状态。

材料坐标始终及时指向同一粒子，该粒子最初位于指定点 (X, Y, Z)。在该坐标系中规定了固体力学的动量方程。另一方面，空间坐标系中的一点 (x, y, z) 是指当前状态下可能位于此处的任何粒子，在该坐标系中规定了热方程。

在这两个坐标系中，与体积相关的物理量具有不同的值。例如，在没有任何质量源的情况下，材料坐标的密度在变换之前和之后保持恒定，而空间坐标的密度随体积变化而变化。因此，为了耦合材料坐标（结构力学）上的方程和在空间坐标（传热）上的另一个方程，需要在每个坐标上合理评估这些值。下表列出了一些从材料坐标到空间坐标的热物理量的转换。这些转换涉及变形梯度 \mathbf{F} = {\partial \mathbf{x}} /{\partial \mathbf{X}} 及其行列式 J 。两者均使用由固体力学接口计算的位移场进行评估。

从材料坐标到空间坐标的热物理量转换。

这些转换还反映了应力和应变会通过修改几何结构（如空间坐标所示）影响传热。例如，延伸的边界更可能通过辐射接收更多的热量（Q_\mathrm{r} > Q_{\mathrm{r}, 0}）， 如下图所示。

在固体上表面，初始状态（左）和拉伸上表面（右）后受到的辐射热通量。

另一个示例是空间坐标中的导热系数表达式，通常使用初始状态值 \bold{k}_0 ，以及和固体应变有关的数量  \mathbf{F} 和 J 。

固体变形后空间坐标上导热系数的变化。

与温度相关的应力和应变

在材料坐标中定义固体力学方程。通过线性动量平衡方程和应力-应变关系将位移 \mathbf{u}、第二类皮奥拉-基尔霍夫应力张量 \mathbf{S} ，以及弹性应变张量 \mathbf{E}_\mathrm{el} 相关联。

(1)

(2)

其中， \mathbf{C} 是弹性张量，通常由杨氏模量和泊松系数定义。如果是碳钢 1020，则可能由温度决定。

碳钢 1020 的杨氏模量取决于温度。

下列公式说明，在没有任何塑性效应的情况下，弹性应变张量 \mathbf{E}_\mathrm{el} 通过热应变张量 \mathbf{E}_\mathrm{th} 承载温度依赖性：

(3)

(4)

(5)

热膨胀系数 \boldsymbol{\alpha} ，表征材料由于温度变化而收缩和膨胀的能力。通常它是标量，但可能更常采用张量形式。下表列出了各向同性 \boldsymbol{\alpha} 的一些典型值。

某些材料的热膨胀系数。

此外，\boldsymbol{\alpha} 本身可能取决于温度，如下例所示。

碳钢 1020 的热膨胀系数，取决于温度。

从这些示例中可以看出， \boldsymbol{\alpha} 的值通常是 10^-5 K^-1。因此，为了使 \mathbf{E}_\mathrm{th} 变得显著，与参考温度状态具有高温差是必要的。例如，铝需要达到高于参考温度约 500K 才能显示出仅 1.2％ 的热伸长率。

采用 1:1 变形比例的加热至 500K 受约束铝梁的热膨胀示例。

注意，在等式（3）-（5）的公式中，从总应变中减去热应变。由于较小的 值，热应变值通常是小应变的适当近似值。适用于较大热应变的更精确的乘法公式如下所示，后文不再赘述。该公式用于 COMSOL Multiphysics 中的超弹性材料。

(6)

(7)

(8)

变形固体的热方程

热方程是从热力学第一定律推导出的能量平衡方程。对于固体，在空间坐标上采用以下形式：

(9)

耦合项  W_\sigma 是由于固体压缩或膨胀而产生的热源，其定义为：

(10)

其中，由于 \boldsymbol{\alpha}  与温度无关，可以简化为：

(11)

其中， \boldsymbol{\alpha}  是和 \mathbf{E}_\mathrm{th}
一样的热膨胀系数 。如上表所示， \boldsymbol{\alpha} 的低值必须通过足够高的 T {\mathrm{d} \mathbf{S}} / {\mathrm{d}t} 值作为补偿，以使  W_\sigma  成为重要的热源，即：

• 在高温下
• 通过快速和高度变化的应力

下述为传热与固体力学之间进行多物理场耦合的四个关键贡献：

1. 应变和应力对材料或空间坐标中的热量和边界热通量的影响
2. 弹性矩阵的温度依赖性
3. 通过热应变张量的弹性应变张量的温度依赖性
4. 热源 W_\sigma ，对应于固体中的热弹性阻尼

接下来，我们将说明最后两个耦合贡献，并通过几个建模示例展示如何在 COMSOL Multiphysics 中处理它们。

示例 1：涡轮静叶片中的热应力

我们之前的博客文章中曾详细地描述了如何模拟涡轮静叶片中的热应力。在这里，为了显示 J_\mathrm{th} 的效果，我们只展示了结果。因为这是稳态模型，所以可以忽略热弹性阻尼 W_\sigma。

在材料坐标中表示叶片表面的温度场。

由于处于高温环境，与最初静叶片形状的参考温度 300K 相比，温度场的值介于 870K 和 1100K 之间。如此高的温度使材料更容易发生热变形。平均热膨胀系数和温度约为 1.2·10 ^-5 K ^-1 和 1070K， \mathbf{E}_\mathrm{th} 大约是 0.9％。

由于热效应而产生的体积膨胀，大变形为  \Delta V/V_0 = J_\mathrm{th}-1 （J_\mathrm{th}
  在等式（8）中被引入）。对于小应变，它仍然是一个很好的近似值，约有 2.80％ 的膨胀。在后处理中，实际的体积膨胀为 2.76％。

静叶片的温度场和变形，为了提高可见度，将绘图放大了 3 倍。

示例2：传热支架的瞬态分析

在结构力学模块模型库和案例库中都可以找到支架瞬态分析模型。在此模型中，支架臂根据快速的时间相关负载移动。因此，应该会发生很小的温度变化。

现有模型忽略了这些热效应，因此，我们需要在固体传热接口中增加新的传热。

然后，添加以下两个多物理场特征，以耦合 固体传热 和 固体力学 接口：

• 热膨胀
  • 会修正应用于整个支架域的热应变张量 \mathbf{E}_\mathrm{th} 和热弹性热源 W_\sigma 。
• 温度耦合
  • 将通过固体传热接口和固体力学接口计算出的温度变量耦合在一起

该研究还可以扩展到 30ms，以观察更多的负载周期。

从各处的 20 °C等温曲线开始，小的温度变化导致可忽略的热应变张量。现在，由于快速的应力变化，对热效应的主要贡献是热弹性热源。

随时间变化的支架温度曲线，为了提高可见度，将绘图放大了10倍。

在支架的极端温度之间可以观察到约 0.8K 的差异。如预期的那样，加热和冷却过程位于应力更重要且应力变化更大的拐角处。

通过求解热方程和动量平衡方程，可以数值计算变形固体中的热传递。根据实际情况，我们对两个坐标系加以区分：

1. 制定运动方程的材料坐标
2. 热方程的空间坐标

两个坐标中与体积有关的数量具有不同的值，并且需要相互转换，特别是对于特定的能量和密度。

这两个控制方程均包含耦合项，这些耦合项使固体运动取决于温度，而传热则取决于固体变形。如前两个示例所示，COMSOL Multiphysics 中提供了合适的功能可以描述这些过程。

当温度保持在参考状态附近并且应力没有太快变化时，这些耦合效应可以忽略不计。否则，应将它们添加到模型的公式中。

如果要深入研究该主题，您可以下载此处提到的模型相关文件，并通过下面的链接阅读一些相关的博客文章。

更多资源

• 模型下载：
  • 涡轮静叶片的热应力分析
  • 支架教程模型
• 之前的博客文章：
  • 共轭传热
  • 模拟涡轮静叶片中的热应力

编者注：为了与 COMSOL Multiphysics 5.1 版本保持一致，此博客文章已于 2015 年 7 月 23 日更新。