瞬态热传导问题

关键词： 瞬态 线性 热传导 向后差分

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

控制方程

微分形式

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

这里考虑三类边界条件:

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

第二类边界条件 \(-k\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)\) 其中，$\bar{u}$ 和 $q$ 为边界上的温度和热流值，$a$ 和 $b$ 为边界换热常数。$\mathbf{n}$ 为边界上的单位外法向量。热流 $q$ 以进入区域为正。

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

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

弱形式

运用虚位移原理求温度，由微分方程 $\eqref{eq1}$ 两端乘以温度的变分，并在全域上积分，得 \(\int_{\Omega}\left[\rho c \frac{\partial u}{\partial t}-\frac{\partial}{\partial x}\left(k \frac{\partial u}{\partial x}\right)-\frac{\partial}{\partial y}\left(k \frac{\partial u}{\partial y}\right)-\rho Q\right] \delta u d \Omega=0\) 对于二阶导数的部分，利用分部积分将其中的一阶导数转化到温度的变分上，得 \(\int_{\Omega} \rho c \frac{\partial u}{\partial t} \delta u d \Omega+\int_{\Omega}\left(k \frac{\partial u}{\partial x} \frac{\partial \delta u}{\partial x}+k \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 \frac{\partial u}{\partial \mathbf{n}}\cdot \mathbf{n} \delta u d \Gamma\) 对于分部积分得到的边界积分项，代入边界条件，得 \(\int_{\Omega} \rho c \frac{\partial u}{\partial t} \delta u d \Omega+\int_{\Omega}\left(k \frac{\partial u}{\partial x} \frac{\partial \delta u}{\partial x}+k \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 \dot{U}+SU=F\) 其中， $U$ 为温度向量，$\dot{U}$ 为温度对时间的导数向量，$M$ 为质量矩阵， $S$ 为刚度矩阵，$F$ 为载荷向量。采用向后差分法进行时域离散，可将上式改写为 \(M \frac{U^{t+\Delta t}-U^{t}}{\Delta t}+S U^{t+\Delta t}=F^{t+\Delta t}\) 整理得 \((M+S \Delta t) U^{t+\Delta t}=F^{t+\Delta t} \Delta t+M U^{t}\tag{3}\label{eq3}\) 上标 $t+\Delta t$ 代表当前时间步，$t$ 代表上一时间步。该式即为最终需要求解的线性方程组。

脚本编写

一个具体问题的有限元语言表述，通常需要包含 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 文件，声明程序具体使用的坐标系统、单元类型和积分方法。这里文件命名为 lp.mdi，具体内容如下。

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

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

GCN 文件

GCN 文件给出计算流程信息，系统规定其文件名字要与 MDI 文件一致。线性抛物方程计算流程比较简单，包含一个时间步的循环过程。只需填写少量代码，对应的 GCN 文件如下。

NFE 文件

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

算例1

问题描述

矩形金属薄板，尺寸为 $1\times0.1\ \mathrm{m}^2$，密度为 $8000\ \mathrm{kg}/\mathrm{m}^3$，比热容为 $450\ \mathrm{J}/\mathrm{kg}/^{\circ}\mathrm{C}$，热传导系数为 $80\ \mathrm{W}/\mathrm{m}/{ }^{\circ}\mathrm{C}$，内热源密度为 0。初始温度为 $0 { }^{\circ}\mathrm{C}$。左侧温度固定为 $0 { }^{\circ}\mathrm{C}$，右侧固定为 $10{ }^{\circ}\mathrm{C}$，上下两边绝热。求温度场随时间的变化。

前处理

填写材料参数文件

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

填写时间步长文件

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

100 10000

可知，时间步长为 100 s，计算终止时间为 10000 s。

赋单元材料

赋初始条件

施加第一类边界条件

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

网格剖分

根据所开发程序的计算要求，采用四边形线性单元进行网格剖分。

计算结果

温度场的演化过程

不同时刻温度场的分布

图中曲线从低到高依次为：100 s、2000 s、4000 s、6000 s、8000 s 和 1000 s。

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

算例2

问题描述

矩形金属薄板，同算例 1 的尺寸和材料相同。初始温度为 $0 { }^{\circ}\mathrm{C}$。左侧与环境换热，环境温度为 $0 { }^{\circ}\mathrm{C}$，换热系数为 $80 \mathrm{~W}/\mathrm{m}^2 /{ }^{\circ} \mathrm{C}$。右侧固定为 $10{ }^{\circ}\mathrm{C}$。上下两边绝热。求温度场随时间的变化。

前处理

填写材料参数文件

由于本算例同算例 1 的材料相同，只是变化了边界条件，因此，lp.mat 文件只需要修改边界单元的材料参数。修改后的文件内容如下。

填写时间步长文件

time.dat 文件内容如下：

600 60000

可知，时间步长为 600 s，计算终止时间为 60000 s。

赋单元材料

同算例1。

赋初始条件

同算例1。

施加第一类边界条件

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

施加第三类边界条件

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

网格剖分

同算例1。

计算结果

温度场的演化过程

左端点和中间点温度的变化

左端点为蓝线，中间点为红线。随着时间的增长，均逐渐逼近理论稳定值 $5 { }^{\circ}\mathrm{C}$ 和 $7.5 { }^{\circ}\mathrm{C}$ ，充分证明了算法和程序的有效性。