稳态热传导问题

关键词： 稳态 线性 热传导

稳态热传导问题因其形式简单，常常用作有限元学习的入门范例。本节以二维稳态热传导问题为例，讲解如何使用 FEtch 系统求解线性椭圆型偏微分方程。

操作视频

控制方程

微分形式

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

这里考虑两类边界条件:

第一类边界条件，给定温度 \(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)\)

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

弱形式

运用变分原理，由微分方程 $\eqref{eq1}$ 两端乘以温度的变分，并在全域上积分，得 \(\int_{\Omega}\left(\frac{\partial}{\partial x}\left(k \frac{\partial u}{\partial x}\right)+\frac{\partial}{\partial y}\left(k \frac{\partial u}{\partial y}\right)+\ Q\right) \delta u d \Omega=0\) 对于二阶导数的部分，利用分部积分将其中的一阶导数转化到温度的变分上，得 \(\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} Q \delta u d \Omega+\int_{\Gamma} k \frac{\partial u}{\partial \mathbf{n}}\cdot \mathbf{n} \delta u d \Gamma\) 对于分部积分得到的边界积分项，代入边界条件，得 \(\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} Q \delta u d \Omega+\int_{\Gamma_2} q \delta u d \Gamma \tag{2}\label{eq2}\) 这就是温度场的最终的弱解形式。对式 $\eqref{eq2}$ 进行空间离散化后，可进一步写成矩阵形式 \(SU=F\tag{3}\label{eq3}\) 其中，$S$ 为刚度矩阵，$U$ 为温度向量，$F$ 为载荷向量。该式即为最终需要求解的线性方程组。

脚本编写

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

PDE 和 FBC 文件

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

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

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

MDI 文件

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

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

GCN 文件

GCN 文件给出计算流程信息，系统规定其文件名字要与 MDI 文件一致。作为线性椭圆方程，稳态热传导问题的计算流程非常简单，只需填写少量代码，具体内容如下。

NFE 文件

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

算例1

问题描述

矩形金属薄板，尺寸为 $1.0\ \mathrm{m}\times0.1\ \mathrm{m}$，热传导系数为 $80\ \mathrm{~W}/\mathrm{m}/{ }^{\circ}\mathrm{C}$，内热源密度为 $0$。左侧温度固定为 $0\ { }^{\circ}\mathrm{C}$，右侧固定为 $10\ { }^{\circ}\mathrm{C}$，上下两边绝热。求温度的分布。

前处理

填写材料参数文件

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

赋单元材料号

施加边界条件

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

网格剖分

根据所开发程序的计算要求，采用三角形线性单元进行网格剖分。共78个单元，62个节点。

计算结果

可以看到，温度场左低右高，呈线性分布。这与理论解完全一致，充分证明了算法和程序的有效性。

算例2

问题描述

矩形金属薄板，尺寸为 $0.6\ \mathrm{m}\times1.0\ \mathrm{m}$，热传导系数和内热源密度同算例 1。右侧和上方温度固定为 $0\ { }^{\circ}\mathrm{C}$，左侧绝热。底部施加恒定的热流，热流密度为 $1000\ \mathrm{~W}/\mathrm{m}^{-2}$。求稳定后温度场的分布。

前处理

填写材料参数文件

本算例与算例 1 的材料相同，差别仅在于包含了第二类边界条件。因此，只需要修改与边界单元相关的参数值。修改完成后的 mat 文件内容如下。

赋单元材料号

施加第一类边界条件

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

施加第二类边界条件

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

网格剖分

采用三角形线性单元进行网格剖分，共 524 个单元，295 个节点。

计算结果

可以看出，计算结果与理论解是一致的，进一步证明了算法和程序的有效性。