瞬态非线性热传导问题

关键词： 瞬态 非线性 热传导 向后差分 迭代 多种材料 第三类边界

本节以二维非线性瞬态热传导问题为例，讲解如何使用 FEtch 系统对非线性抛物型方程进行求解，以及如何自定义函数。

控制方程

微分形式

对于二维直角坐标系，瞬态非线性热传导服从如下偏微分方程 \(\rho c(u) \frac{\partial u}{\partial t}-\frac{\partial}{\partial x}\left(k(u) \frac{\partial u}{\partial x}\right)-\frac{\partial}{\partial y}\left(k(u) \frac{\partial u}{\partial y}\right)=\rho Q\ \left(\text{in } \Omega\right) \tag{1}\label{eq1}\) 边界条件为:

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

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

第三类边界条件 \(-k(u)\frac{\partial u}{\partial \mathbf{n}} \cdot \mathbf{n}=a(u-b)\ \left(\text{on } \Gamma_3\right)\) 其中，$u$ 表示温度，$\rho$ 表示材料密度，$c(u)$ 为材料的比热，$k(u)$ 是热传导系数，$Q$ 为热源密度，$\bar{u}$ 和 $q$ 为边界上的温度和热流值，$\mathbf{n}$ 为边界上的单位外法向量，$a$ 和 $b$ 为边界换热常数。$x$ 和 $y$ 为二维直角坐标系下的坐标变量，$t$ 为时间。

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

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

弱形式

运用变分原理，由方程 $\eqref{eq1}$ 得 \(\int_{\Omega}\left(\rho c(u) \frac{\partial u}{\partial t}-\frac{\partial}{\partial x}\left(k(u) \frac{\partial u}{\partial x}\right)-\frac{\partial}{\partial y}\left(k(u) \frac{\partial u}{\partial y}\right)-\rho Q\right) \delta u d \Omega=0\) 将其化为弱形式，可得 \(\int_{\Omega} \rho c(u) \frac{\partial u}{\partial t} \delta u d \Omega+\int_{\Omega}\left(k(u) \frac{\partial u}{\partial x} \frac{\partial \delta u}{\partial x}+k(u) \frac{\partial u}{\partial y} \frac{\partial \delta u}{\partial y}\right) d \Omega \\=\int_{\Omega} \rho Q \delta u d \Omega+\int_{\Gamma} k(u) \frac{\partial u}{\partial \mathbf{n}}\cdot \mathbf{n} \delta u d \Gamma\) 代入边界条件，可得 \(\int_{\Omega} \rho c(u) \frac{\partial u}{\partial t} \delta u d \Omega+\int_{\Omega}\left(k(u) \frac{\partial u}{\partial x} \frac{\partial \delta u}{\partial x}+k(u) \frac{\partial u}{\partial y} \frac{\partial \delta u}{\partial y}\right) d \Omega+\int_{\Gamma_3} au \delta u d \Gamma\\=\int_{\Omega} \rho Q \delta u d \Omega+\int_{\Gamma_2} q \delta u d \Gamma+\int_{\Gamma_3} ab \delta u d \Gamma \tag{2}\label{eq2}\)

离散化

式 $\eqref{eq2}$ 进行空间离散化后，可写成矩阵形式 \(M(U) \dot{U}+S(U)U=F\) 其中 $U$ 为温度向量，$\dot{U}$ 为温度对时间的导数向量，$M(U)$ 为质量矩阵， $S(U)$ 为刚度矩阵，$F$ 为载荷向量。采用向后差分法进行时域离散，可将上式改写为 \(M(U) \frac{U^{t+\Delta t}-U^{t}}{\Delta t}+S(U) U^{t+\Delta t}=F^{t+\Delta t}\) 整理得 \(\left(M(U)+S(U) \Delta t\right) U^{t+\Delta t}=F^{t+\Delta t} \Delta t+M(U) U^{t} \tag{3}\label{eq3}\) 其中，上标 $t+\Delta t$ 代表当前时间步，$t$ 代表上一时间步。$M(U)$ 和 $S(U)$ 均为 $U$ 的函数，因此，式 $\eqref{eq3}$ 是一个非线性方程组，无法直接求解，需要引入迭代过程。

最常用的处理方法是采用直接迭代法（也称 Picard 法），将式 $\eqref{eq3}$ 线性化，得 \(\left(M(U^i)+S(U^i) \Delta t\right) U^{i+1,t+\Delta t}=F^{t+\Delta t} \Delta t+M(U^i) U^{t} \tag{4}\label{eq4}\)

其中，上标 $i$ 表示 $t+\Delta t$ 时间步下的迭代步数。上式在每个时间步下都需要反复地迭代计算，直至误差 $\Vert U^{i+1}-U^i\Vert$ 小于某一设定的阈值。直接迭代法可以通过使用松弛因子，有效地提高收敛速度。松弛迭代法已被写入 NFE 算法模板库，是 FEtch 系统默认的非线性迭代算法。

可以看出，迭代过程的引入是非线性问题较线性问题计算量明显增大的主要原因。因此，如何在保证计算稳定性的同时提高收敛速度，一直是非线性求解领域永恒的话题。

材料本构模型

非线性问题远比线性问题要复杂。由于非线性参数的函数形式千变万化，所以通常需要针对具体的问题专门地设计程序。对于 FEtch 系统而言，由于脚本文件是完全自主可控的，只要基本框架搭好，后面需要改动的地方极少。

为后续算例的计算需要，这里考虑铸钢这种非线性材料。它的比热和密度的乘积为

导热系数为

可见，铸钢的材料参数包含了两个依赖于自变量温度 $u$ 的函数，$\rho c(u)$ 和 $k(u)$ 。用户需认真体会有限元语言是如何定义和使用这些函数的。

脚本编写

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

PDE 和 FBC 文件

在这个例子中，只有一个物理场，需要给出温度场对应的 PDE 文件和 FBC 文件。根据有限元语言的语法规则，将方程 $\eqref{eq2}$ 的体积分项写入 PDE 文件如下。

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

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

MDI 文件

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

注意，同一个场的 PDE 和 FBC 的单元类型和积分方法需要保证匹配，否则会导致计算精度降低，甚至计算失败的情况发生。MDI 文件声明的单元类型和积分方式将以参数的形式传递给 PDE 和 FBC 文件中的的 %1、%2 和 %3。如果后续想要使用其他的单元类型或者积分方法，只需修改 MDI 文件即可。

GCN 文件

GCN 文件给出计算流程信息，系统规定其文件名字要与 MDI 文件一致。非线性抛物方程的计算流程较为复杂，包含两大循环过程，时间步循环和迭代步循环。对应的 GCN 文件和说明如下。

NFE 文件

GCN 文件使用了 nparb 算法，系统会基于模板文件自动生成以下的 npa.nfe 文件。如果用户有其他需要，可以在此基础上进一步地修改。

算例

这里给出一个经典算例，在 FEPG 的文献资料中经常被引用。《 FEPG 有限元应用深入剖析》一书的内容就是围绕本算例展开的。

一钢铸件及其砂模的横截面尺寸如图所示。砂模的导热系数是 $2.5\times10^{-2}$，比热和密度的乘积是 $7.112\times10^{-2}$；铸钢的材料参数已写入程序内部，材料号标记为 2，可直接调用。铸钢的初始温度为 $2875 { }^{\circ} \mathrm{C}$，砂模的初始温度为 $80 { }^{\circ} \mathrm{C}$；砂模的外边界为对流边界条件，其中对流系数为 $0.014$，环境温度为 $80 { }^{\circ} \mathrm{C}$。求 3 个时刻内铸钢及砂模的温度场的演化过程。

前处理

填写材料参数文件

系统自动生成了 np.mat 文件，用于设置材料参数。参数初始值全部为 PDE 和 FBC 文件给出的默认值。可根据具体问题进行修改。对于本算例，修改完成后的文件内容如下。

注意，对于铸钢材料，由于它的参数是由自定义函数计算得到的，仅需要保证标记参数 itag 的值为 2 即可，其他的参数并不发挥作用，可随意取值。

填写时间步长文件

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

0.001 3

可知， 初始时间步长为 0.001，结束时间为 3。

赋单元材料

赋初始条件

为避免歧义，继续明确定义交界位置的初始条件。

施加第三类边界条件

第二类边界条件以线单元的形式施加

网格剖分

根据所开发程序的计算要求，网格剖分采用 4 节点等参单元。网格大小为 $0.5\times 0.5$，共使用了 864 个单元，925 个节点。

计算结果

通过自动改变时间步长，变化范围从 0.001 到 0.059，总共使用了 211 个时间步。对计算结果稍加整理，效果如下。

温度场的演化过程

可以看到，伴随着传热过程，铸钢温度逐渐下降，向着全域等温的稳态发展。整个变化过程是平滑和稳定的，在一定程度上证明了算法和程序的有效性。

小结

对本章的计算思路做一下总结。对于瞬态非线性热传导问题，当材料的物性参数随温度变化时，离散方程式的系数也随温度变化，故每次迭代计算时需要更新物性参数。在每个时间步里温度场收敛了，才能进入下一时间步的迭代。在某个时间步内的迭代计算思路为：

1. 根据上一次迭代计算得到的温度（第一次迭代之前应该是初始温度） 计算（更新）物性参数；
2. 根据新物性参数计算温度场；
3. 计算当前时间步长内结果是否收敛；
4. 若收敛，进入下一时间步长迭代计算；若不收敛，继续更新物性参数，迭代计算温度场，直到收敛。

以上过程的实现需要 PDE、NFE 和 GCN 这 3 种文件的相互协调配合。建议用户仔细研读本章用到的脚本文件的源代码，充分理解后再去尝试求解其他的瞬态非线性问题，稳扎稳打，避免踩坑。