热固耦合问题

关键词： 瞬态 线性 热传导 弹性力学 向后差分 热固耦合 解耦合

结构内温度场发生变化时，若受到外部约束或温度场不均匀，就会产生一定的应力，称为温度应力。温度应力的出现，已引起工程上普遍关注。在土木工程中，钢筋混凝土构筑物由于受到环境温度变化的影响，表面和内部会产生变形，若遇约束，会引起温度应力， 当应力达到一定值时，结构内部产生微观裂缝，甚至发展为裂缝；大体积混凝土地基和大坝结构由于水泥浇注期内水化热的作用，冷却收缩时温度应力若超过材料抗拉强度，也会呈现裂缝；道路和桥梁由于气候多变，在路面荷载和不均匀温度场共同作用下也有可能开裂受损。由此可见，不均匀温度场和温度应力计算具有重要意义，其结果可以直接为工程设计提供依据。

本节以二维弹性体的热固耦合问题为例，讲解如何使用 FEtch 系统对热固耦合问题进行解耦合求解，即如何先求温度场，再求位移场和应力场。

控制方程

微分形式

瞬态热传导

二维直角坐标系下，瞬态热传导服从如下偏微分方程 \(\rho c \frac{\partial T}{\partial t}-\frac{\partial}{\partial x}\left(k \frac{\partial T}{\partial x}\right)-\frac{\partial}{\partial y}\left(k \frac{\partial T}{\partial y}\right)=0 \ \left(\text{in } \Omega\right) \tag{1}\label{eq1}\) 其中，$\rho$ 表示材料密度，$c$ 为材料的比热，$T$ 表示温度，$k$ 是热传导系数，$x$ 和 $y$ 为二维直角坐标系下的坐标变量。

这里考虑两类边界条件:

第一类边界条件 \(T=\bar{T} \ \left(\text{on } \Gamma_1\right)\)

第二类边界条件 \(-k\frac{\partial T}{\partial \mathbf{n}} \cdot \mathbf{n}=-q\ \left(\text{on } \Gamma_2\right)\)

其中，$\bar{T}$ 和 $q$ 为边界上的温度和热流值。$\mathbf{n}$ 为边界上的单位外法向量。热流 $q$ 以进入区域为正。

初始条件为 \(T=T_0 \ \left(\text{in } \Omega\right)\)

其中，$T_0$ 为求解域内初始的温度分布。

热弹性力学

二维直角坐标系下，忽略惯性力和体力的作用，热弹性力学方程可写成如下形式 \(\frac{\partial \sigma_{x x}}{\partial x}+\frac{\partial \sigma_{x y}}{\partial y} = 0 \ \left(\text{in } \Omega\right) \tag{2}\label{eq2}\) \(\frac{\partial \sigma_{x y}}{\partial x}+\frac{\partial \sigma_{y y}}{\partial y} = 0 \ \left(\text{in } \Omega\right) \tag{3}\label{eq3}\) 几何方程为 \(\varepsilon_{xx}=\frac{\partial u}{\partial x}, \quad \varepsilon_{yy}=\frac{\partial v}{\partial y}, \quad \varepsilon_{x y}=\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\) 本构方程（平面应力）为 \(\left(\begin{array}{c}\sigma_{x x} \\ \sigma_{y y} \\ \sigma_{x y} \end{array}\right)=\frac{E}{(1+\nu)(1-\nu)}\left(\begin{array}{ccc} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & (1-\nu) / 2 \end{array}\right) \left(\begin{array}{l}\varepsilon_{x x} \\ \varepsilon_{y y} \\ \varepsilon_{x y}\end{array}\right) -\frac{E}{(1-\nu)}\alpha T\left(\begin{array}{l} 1\\ 1 \\ 0 \end{array}\right)\) 其中， $u$、$v$ 为位移 $\boldsymbol{u}$ 的分量。 $\varepsilon_{x x}$、$\varepsilon_{y y}$ 、 $\varepsilon_{x y}$ 为应变 $\boldsymbol{\epsilon}$ 的分量，$\sigma_{x x}$ 、$\sigma_{y y}$ 、$\sigma_{x y}$ 为应力 $\boldsymbol{\sigma}$ 的分量。参数 $E$ 为弹性模量， $\nu$ 为泊松比，$\alpha$ 为线膨胀系数。

考虑两类边界条件:

第一类边界条件，给定位移 \(u=u_0,\ v=v_0\ \left(\text{on } \Gamma_1\right)\)

第二类边界条件，给定外力 \(T_{x}=\left(\sigma_{x x},\sigma_{x y}\right)\cdot \mathbf{n}=t_x,\ T_{y}=\left(\sigma_{x y},\sigma_{y y}\right)\cdot \mathbf{n}=t_y\ \left(\text{on } \Gamma_2\right)\)

其中，$u_0$、$v_0$为边界上的位移；$T_{x}$、$T_{y}$ 分别为边界力在 $x$ 和 $y$ 方向的分量。$\mathbf{n}$ 为边界上的单位外法向量。

弱形式

PDE 文件的编写依赖于待解偏微分方程的弱解方程。对热传导方程 $\eqref{eq1}$ 两边分别乘以温度的变分并在全域上积分，得 \(\int_{\Omega}\left[\rho c \frac{\partial T}{\partial t}-\frac{\partial}{\partial x}\left(k \frac{\partial T}{\partial x}\right)-\frac{\partial}{\partial y}\left(k \frac{\partial T}{\partial y}\right)\right] \delta T d \Omega=0\) 对于二阶导数的部分，利用分部积分将其中的一阶导数转化到温度的变分上，并继续代入边界条件，得 \(\int_{\Omega}\rho c \frac{\partial T}{\partial t} \delta T d \Omega+\int_{\Omega}\left(k \frac{\partial T}{\partial x} \frac{\partial \delta T}{\partial x}+k \frac{\partial T}{\partial y} \frac{\partial \delta T}{\partial y}\right) d \Omega=\int_{\Gamma_2} q \delta u d \Gamma \tag{4}\label{eq4}\) 该式即为温度场的弱形式。

根据虚位移原理，对力学平衡方程 $\eqref{eq2}$ 和 $\eqref{eq3}$ 两边分别乘以位移的变分并在求解区域内积分，得到其等效积分形式 \(\int_{\Omega}\left[\left(\frac{\partial \sigma_{x x}}{\partial x}+\frac{\partial \sigma_{x y}}{\partial y}\right) \delta u+\left(\frac{\partial \sigma_{x y}}{\partial x}+\frac{\partial \sigma_{y y}}{\partial y}\right) \delta v\right] d\Omega=0\) 对上式进行分部积分，并代入边界条件，得到弱形式 \(\int_{\Omega}\left(\sigma_{x x} \delta \varepsilon_{x x}+\sigma_{y y} \delta \varepsilon_{y y}+\sigma_{x y} \delta \varepsilon_{x y}\right)d\Omega\\=\int_{\Gamma_2}\left(t_{x} \delta u+t_{y} \delta v\right) d\Gamma\)

具体展开为 \(\int_{\Omega}\left(f_1\varepsilon_{x x} \delta \varepsilon_{x x}+f_1\varepsilon_{y y} \delta \varepsilon_{y y}+f_1\nu\varepsilon_{x x} \delta \varepsilon_{y y}+f_1\nu\varepsilon_{y y} \delta \varepsilon_{x x}+f_1f_2\varepsilon_{x y} \delta \varepsilon_{x y}\right)d\Omega\\=\int_{\Omega}\left(f_3\alpha T \delta \varepsilon_{x x}+f_3\alpha T \delta \varepsilon_{y y}\right)d\Omega+\int_{\Gamma_2}\left(t_{x} \delta u+t_{y} \delta v\right) d\Gamma \tag{5}\label{eq5}\) 其中， $f_1=\frac{E}{(1+\nu)(1-\nu)}$，$f_2=\frac{1-\nu}2$，$f_3=\frac{E}{(1-\nu)}$，均为与 $E$ 和 $\nu$ 相关的参数。式 $\eqref{eq5}$ 就是以位移为基本未知量的弱形式。这里温度 $T$ 是以已知量的形式出现的，作为荷载项出现在了方程的右端。

在求得了温度和位移之后，利用本构方程和几何方程可以得到应力与位移、温度的关系 $\boldsymbol{\sigma}(\boldsymbol{u}，T)$。应用最小二乘法，以待求应力 $\boldsymbol{\sigma}$ 与由位移和温度求得的应力 $\boldsymbol{\sigma}_0(\boldsymbol{u},T)$ 的偏差为目标泛函，即 \(F(\boldsymbol{\sigma})=\int_{\Omega}[\boldsymbol{\sigma-\sigma_0}(\boldsymbol{u},T)]^{2} d \Omega\) 对 $F(\boldsymbol{\sigma})$ 取极值，由变分原理知需满足 $\delta F(\sigma) = 0$。于是有 \(\delta F(\boldsymbol{\sigma})=\int_{\Omega} 2 \delta \sigma^{\mathrm{T}}[\boldsymbol{\sigma}-\boldsymbol{\sigma}_0(\boldsymbol{u},T)] \mathrm{d} \Omega=0\) 整理得 \(\int_{\Omega} \boldsymbol{\sigma} \delta \boldsymbol{\sigma} d\Omega=\int_{\Omega} \boldsymbol{\sigma_{0}}(\boldsymbol{u},T) \delta \boldsymbol{\sigma} d\Omega \tag{6}\label{eq6}\)

该式即为应力场的弱形式。

脚本编写

求解策略

理论上，热固耦合问题通常耦合强度较弱，可以采用解耦合的策略，对温度场、位移场和应力场依次求解。具体步骤如下：

1. 温度场是完全独立的。对式 $\eqref{eq4}$ 进行空间离散化后，采用向后差分法进行时域离散。继续求解线性方程组，即可求得温度场。
2. 将温度场作为已知量代入式 $\eqref{eq5}$ 右侧，进行空间离散化后求解线性方程组，即可求得位移场。
3. 同样地，继续将温度场和位移场作为已知量代入式 $\eqref{eq6}$ 右侧，即可继续求解应力场。

位移场和应力场完全依赖于温度场的结果，用户可以根据需要自行决定何时计算它们，如每个时间步都计算，或者在温度场稳定后再计算。

一个具体问题的有限元语言表述，通常需要包含 5 种类型的文件：PDE、FBC、MDI、GCN 和 NFE。下面依次给出详细介绍。

PDE 和 FBC 文件

温度场

根据有限元语言的语法规则，将方程 $\eqref{eq4}$ 的体积分项写入 PDE 文件如下。

这里将 shap 行和 gaus 行的参数设为待定，参数值 %1、%2 和 %3 将从 MDI 文件中取得。

同样地，将方程 $\eqref{eq4}$ 的边界积分项写入 FBC 文件如下。

位移场

基于有限元语言语法，将式 $\eqref{eq5}$ 的体积分部分，编写成了disp.pde 文件，具体内容如下。

这里将 shap 行和 gaus 行的参数设为待定，参数值 %1、%2 和 %3 将从 MDI 文件中取得。

对于单元载荷向量中的边界积分部分，用 FBC 文件来描述（FBC 的命名习惯上常取对应 PDE 的文件名，disp.fbc）。将方程 $\eqref{eq1}$ 的边界积分项写入 FBC 文件如下。

应力场

应力场只有体积分项而没有边界积分项，所以只需填写 PDE 文件。对应式 $\eqref{eq2}$，利用有限元语言语法编写的 PDE 文件如下（命名为 sdisp.pde）

注意，后面会看到，由于应力场计算直接使用了位移场的数据作为输入数据，所以，这里应力场 PDE 文件中 mate 行材料参数的定义必须和位移场 PDE 文件 mate 行的定义完全相同。

MDI 文件

MDI 文件负责组织不同场的 PDE 和 FBC 文件，声明程序具体使用的坐标系统、单元类型和积分方法。这里文件命名为 tds.mdi，具体内容如下。

注意，同一个场的 PDE 和 FBC 的单元类型和积分方法需要保证匹配，否则会导致计算精度降低，甚至计算失败的情况发生。

对于弹性力学问题，理论上四边形单元的计算精度要明显好于三角形单元。因此，这里配合2点高斯积分，采用了四边形单元。

MDI 文件声明的单元类型和积分方式将以参数的形式传递给 PDE 和 FBC 文件中的的 %1、%2 和 %3。如果后续想要使用其他的单元类型或者积分方法，只需修改 MDI 文件即可。

GCN 文件

在这个例子中，有位移场和应力场两个场。除了对应的 PDE 和 FBC 文件，还要给出每个场使用的算法文件和计算的命令流程，这些信息由 GCN 文件给出。系统要求 GCN 与 MDI 文件同名，具体内容如下。

注意：stress 求解器是一种特殊化了的求解器，专门用于求解已知场的梯度场。它与 str.nfe 算法模板文件互相匹配。stress 求解器不必初始化（即不使用 #start*），而是直接使用与之耦合的 a 场的网格和材料数据作为输入数据。因此，后面我们会看到，GiD 前处理模板中没有为 b 场添加菜单选项。

NFE 文件

根据 GCN 文件 defi 段的定义，初始化后，系统基于模板文件自动生成了各个场的 NFE 文件。因为 b 场和 c 场耦合了其他场的计算结果，用户需要根据具体情况进行一定的修改。

温度场

a 场使用了 parb 算法，系统自动生成了以下的 tdsa.nfe 文件。为了使后处理时的图例显示更加规范，这里稍微改动了一下 gidpost 语句。由于 a 场是独立的，并未使用其他场的数据，因此该 NFE 文件其他语句完全可以不作任何改动。

位移场

通过已知温度求位移，b 场使用了 ell 算法，系统基于模板文件自动生成了 tdsb.nfe 文件。由于 b 场使用了 a 场的计算结果，因此，需添加 coef 段和相关的数据读取语句。修改后文件的具体内容如下。

应力场

通过已知位移和温度求应力，c 场使用了 str 算法，系统基于模板文件自动生成了 tdsc.nfe 文件。由于 c 场同时使用了 a 场和 b 场的计算结果，因此，需修改 coef 段和相关的数据读取语句。修改后文件的具体内容如下。

注意，NFE 文件的 coef 行，要保证与 PDE 文件中 coef 行耦合场变量的数目和顺序保持一致，否则会出现信息传递错误。另外，要特别关注临时中间文件 unod 和 unodb 等的读写操作，读和写内容必须相互匹配。

算例

边长为 $0.4\ \mathrm{m}\times0.3\ \mathrm{m}$ 的长方形薄板，上下边的位移竖向约束，横向自由，两侧边自由。初始温度为 $0 { }^{\circ}\mathrm{C}$，底部突然施加 $10 { }^{\circ}\mathrm{C}$ 的固定温度，上边和两侧绝热，求 $100\ \mathrm{s}$ 内薄板的温度、变形和应力的变化情况。具体参数值为 $\rho=5000\ \mathrm{kg}/\mathrm{m}^3$，$c=10\ \mathrm{J}/\mathrm{kg}/^{\circ}\mathrm{C}$，$k=100\ \mathrm{W}/\mathrm{m}/{ }^{\circ}\mathrm{C}$，$E＝1.0×10^{10}\ \mathrm{N/m^2}$，$\nu = 0.3$，$\alpha=1.0\times10^{-5}\ /^{\circ}\mathrm{C}$。

前处理

根据对称性，取薄板的右半边进行建模。

填写材料参数文件

系统自动生成了 tds.mat 文件，用于设置材料参数。对于本算例，修改完成后的文件内容如下。由于应力场计算直接使用了位移场的数据作为输入数据，所以没有相关的单元和参数选项。

填写时间步文件

由 tds.gcn 文件可知，用户自定义了时间步长文件 time.dat，给出时间步长和终止时间。文件内容如下：

1 100

可知， 时间步长为 1 s，结束时间为 100 s。

赋单元材料

温度场

位移场

赋初始条件

只有温度场存在初始值。

施加第一类边界条件

温度场

位移场

其中：
u-I 为边界类型，-1 表示约束，1 表示自由
u-D 为相应的约束值或施加节点力

由于有不同的位移边界条件交接，角点位置上的位移边界条件可能出现歧义。为此，继续明确定义角点处的边界条件。

网格剖分

根据所开发程序的计算要求，网格剖分采用 4 节点等参单元。共使用了 150 个单元，176 个节点。

计算结果

对计算结果稍加整理，效果如下。

温度的演化过程和薄板变形情况

随着传热的进行，温度自下而上逐渐增大。薄板从下到上逐渐膨胀。为了更好地显示变形效果，动画中的位移放大了 800 倍。

位移的演化过程

水平位移

竖向位移

薄板从下到上逐渐膨胀。最终，水平向位移成线性分布，而竖向由于受到约束，不发生位移变化。

应力的演化过程

x 方向正应力

y 方向正应力

切应力

进一步分析可以发现数值解与解析解符合得很好，进而证明算法和程序的有效性。这里不再给出具体过程，用户可自行验证。