二维稳态热传导基本方程的有限元方程

本文使用

创作并发布

本篇文章以笛卡尔坐标系下二维稳态热传导基本方程的有限元求解为例，所发展的方法也适用于理想流体流动、电场和对流-扩散等问题的求解。

符号说明

\vec{\cdot}——向量
\vec{i}, \vec{j}——二维笛卡尔坐标系基矢量
\{\cdot\}——列向量(向量分量的矩阵表示)[^1]
<\cdot>——行向量(向量分量的矩阵表示)[^1]
\boldsymbol{\nabla}——梯度算子 \{\nabla\}——梯度算子的矩阵表示
s——热源, 单位为 \frac{J}{m^2 \cdot s} (2D)
\vec q——热通量, 单位为 \frac{J}{m \cdot s} (2D)
\vec n——边界上的单位外法向
\boldsymbol{D} ——热传导张量, 单位为 \frac{W}{K} (2D)
[D] ——热传导张量分量的矩阵表达, 单位为 \frac{W}{K} (2D)
T ——温度, 单位为 K
\bar q ——给定边界通量, 单位为\frac{J}{m \cdot s} (2D)
\bar T ——给定边界温度, 单位为 K

热传导基本方程

张量记法

热传导基本方程的张量记法为[^2][^3]：

• \Omega域内的每一点满足能量守恒方程

\boldsymbol{\nabla} \cdot \vec{q} - s = 0 \tag{1}

• \Omega域内满足傅里叶定律

\vec{q} = -\boldsymbol{D}\boldsymbol{\nabla}T \tag{2}

• \Gamma_{q}边界上满足自然边界条件

q_n = \vec{q} \cdot \vec{n} = \bar{q} \tag{3}

• \Gamma_{T}边界上满足本质边界条件

  T = \bar{T} \tag{4}

矩阵记法

其实，如果采用矩阵记法，就隐式地采用了某种坐标系，这里我们选择最常用笛卡尔坐标系。那么则有

\boldsymbol{\nabla} = \frac{\partial }{\partial x}\vec i + \frac{\partial }{\partial y}\vec j \tag{5}

\vec{q} = q_x \vec i + q_y \vec j \tag{6}

\vec{n} = n_x \vec i + n_y \vec j \tag{7}

相应的分量的矩阵表达为：

\{ \nabla \} = \begin{bmatrix} \frac{\partial}{\partial x} \\ \frac{\partial}{\partial y} \\ \end{bmatrix} \tag{8}

\{q\} = \begin{bmatrix} q_x \\ q_y \\ \end{bmatrix} \tag{9}

\{n\} =\begin{bmatrix} n_x \\ n_y \\ \end{bmatrix} \tag{10}

热传导张量\boldsymbol{D}的分量的矩阵表达为：

[D] = \begin{bmatrix} d_{xx} & d_{xy} \\ d_{yx} & d_{yy} \\ \end{bmatrix} \tag{11}

那么式(1)-(4)用矩阵的形式可以表述如下：

• \Omega域内满足能量守恒方程[^4]

\{ \nabla \}^T\{q\} - s = 0 \tag{12}

• \Omega域内满足傅里叶定律

\{q\} = -[D]\{ \nabla T\} \tag{13}

• \Gamma_{n}边界上满足自然边界条件

q_n = \{q\}^{T} \{n\} = \bar{q} \tag{14}

• \Gamma_{e}边界上满足本质边界条件

  T = \bar{T} \tag{15}

弱形式

这里我们直接给出等价于式(1)-(4)的弱形式[^5]，感兴趣的读者可自行查阅文献[^2]。像上面一样，我们也会给出相应的矩阵表述，为了方便以后对有限元进行编程。

张量记法

对任意的\omega \epsilon U_0，寻找T \epsilon U，使其满足

\int_{\Omega} \boldsymbol{\nabla} \omega \cdot \vec{q} d\Omega = \int_{\Gamma_n} \omega \bar{q} d\Gamma-\int_{\Omega} \omega s d\Omega \tag{16}

其中：

U = \{ T(x,y)| T(x,y) \epsilon H^1, T_{\Gamma_{e}} = \bar{T} \}

U_0 = \{ \omega(x,y)| \omega(x,y) \epsilon H^1, \omega_{\Gamma_{e}} = 0 \}

矩阵记法

对任意的\omega \epsilon U_0，寻找T \epsilon U，使其满足

\int_{\Omega} \{\nabla \omega\}^{T} \{q\} d\Omega = \int_{\Gamma_n} \omega \bar{q} d\Gamma-\int_{\Omega} \omega s d\Omega \tag{17}

上述弱形式(16)或(17)对线性材料和非线性材料均适用。如果材料本构满足傅里叶定律(13)，则有：

\int_{\Omega} \{\nabla \omega\}^{T} [D] \{\nabla T\} d\Omega =-\int_{\Gamma_n} \omega \bar{q} d\Gamma + \int_{\Omega} \omega s d\Omega \tag{18}

伽辽金有限元方程

采用有限元方法求解式(18)时，首先需要选择某种单元对\Omega在空间上进行离散。在二维问题中常用单元有三角形单元和四边形单元。这里我们不给出某种具体单元，仅仅是为了推导出通用的有限元方程。假设我们选择的单元含有n_{en}个节点数，离散后取其中一个单元e进行分析。

在单元e的内部，也满足式(18)，即[^6]：

\int_{\Omega^e} \{\nabla \omega^e\}^{T} [D] ^e \{\nabla T^e\} d\Omega =-\int_{\Gamma_n^e} \omega \bar{q} d\Gamma + \int_{\Omega^e} \omega^e s d\Omega \tag{19}

单元上温度T的近似函数为：

T^e(x, y) = \begin{bmatrix} N_1^e(x, y) & N_2^e(x, y) & \cdots & N_{n_{en}}^e(x, y) \end{bmatrix} \begin{bmatrix} T_1^e \\ T_2^e \\ \vdots\\ T_{n_{en}}^e \end{bmatrix} \tag{20}

由式(20)得到：

\{\nabla T^e\} =\begin{bmatrix} \frac{\partial N_1^e}{\partial x} & \frac{\partial N_2^e}{\partial x} & \cdots & \frac{\partial N_{n_{en}}^e}{\partial x}\\ \frac{\partial N_1^e}{\partial y} & \frac{\partial N_2^e}{\partial y} & \cdots & \frac{\partial N_{n_{en}}^e}{\partial y} \\ \end{bmatrix} \begin{bmatrix} T_1^e \\ T_2^e \\ \vdots\\ T_{n_{en}}^e \end{bmatrix} \tag{21}

对于伽辽金法，权函数取为形函数，且与未知数的个数相同。比如取第i个形函数作为权函数，将式(21)和权函数带入式(20)可得：

LHS = \int_{\Omega^e} [\frac{\partial N_i^e}{\partial x} \frac{\partial N_i^e}{\partial y} ] [D ^e] \begin{bmatrix} \frac{\partial N_1^e}{\partial x} & \frac{\partial N_2^e}{\partial x} & \cdots & \frac{\partial N_{n_{en}}^e}{\partial x}\\ \frac{\partial N_1^e}{\partial y} & \frac{\partial N_2^e}{\partial y} & \cdots & \frac{\partial N_{n_{en}}^e}{\partial y} \\ \end{bmatrix} \begin{bmatrix} T_1^e \\ T_2^e \\ \vdots\\ T_{n_{en}}^e \end{bmatrix} d\Omega \tag{22}

RHS = -\int_{\Gamma_n^e} N_i^e \bar{q} d\Gamma + \int_{\Omega^e} N_i^e s d\Omega \tag{23}

把n_{en}个方程都写出来，则有：

LHS = \int_{\Omega^e} \begin{bmatrix} \frac{\partial N_1^e}{\partial x} & \frac{\partial N_1^e}{\partial y} \\ \frac{\partial N_2^e}{\partial x} & \frac{\partial N_2^e}{\partial y} \\ \vdots \\ \frac{\partial N_{n_{en}}^e}{\partial x} & \frac{\partial N_{n_{en}}^e}{\partial y} \end{bmatrix} [D ^e] \begin{bmatrix} \frac{\partial N_1^e}{\partial x} & \frac{\partial N_2^e}{\partial x} & \cdots & \frac{\partial N_{n_{en}}^e}{\partial x}\\ \frac{\partial N_1^e}{\partial y} & \frac{\partial N_2^e}{\partial y} & \cdots & \frac{\partial N_{n_{en}}^e}{\partial y} \\ \end{bmatrix} \begin{bmatrix} T_1^e \\ T_2^e \\ \vdots\\ T_{n_{en}}^e \end{bmatrix} d\Omega \tag{24}

RHS = -\int_{\Gamma_n^e} \begin{bmatrix} N_1^e \\ N_1^e \\ \vdots\\ N_{n_{en}}^e \end{bmatrix} \bar{q} d\Gamma + \int_{\Omega^e} \begin{bmatrix} N_1^e \\ N_1^e \\ \vdots\\ N_{n_{en}}^e \end{bmatrix} s d\Omega \tag{25}

由于单元节点上的温度独立于空间坐标，所以可以提取到积分号外，则有：

LHS = \int_{\Omega^e} \begin{bmatrix} \frac{\partial N_1^e}{\partial x} & \frac{\partial N_1^e}{\partial y} \\ \frac{\partial N_2^e}{\partial x} & \frac{\partial N_2^e}{\partial y} \\ \vdots \\ \frac{\partial N_{n_{en}}^e}{\partial x} & \frac{\partial N_{n_{en}}^e}{\partial y} \end{bmatrix} [D ^e] \begin{bmatrix} \frac{\partial N_1^e}{\partial x} & \frac{\partial N_2^e}{\partial x} & \cdots & \frac{\partial N_{n_{en}}^e}{\partial x}\\ \frac{\partial N_1^e}{\partial y} & \frac{\partial N_2^e}{\partial y} & \cdots & \frac{\partial N_{n_{en}}^e}{\partial y} \\ \end{bmatrix} d\Omega \begin{bmatrix} T_1^e \\ T_2^e \\ \vdots\\ T_{n_{en}}^e \end{bmatrix} \tag{26}

进一步，借助矩阵符号，式(26)和(25)可写成紧凑格式，如下：

LHS = \int_{\Omega^e} [B^e]^T [D ^e] [B^e] d\Omega \{d\}^e \tag{27}

RHS = -\int_{\Gamma_n^e} <{N^e}>^T \bar{q} d\Gamma + \int_{\Omega^e} <{N^e}>^T s d\Omega \tag{28}

其中：

[B^e] = \begin{bmatrix} \frac{\partial N_1^e}{\partial x} & \frac{\partial N_2^e}{\partial x} & \cdots & \frac{\partial N_{n_{en}}^e}{\partial x}\\ \frac{\partial N_1^e}{\partial y} & \frac{\partial N_2^e}{\partial y} & \cdots & \frac{\partial N_{n_{en}}^e}{\partial y} \\ \end{bmatrix}

\{d^e\} = \begin{bmatrix} T_1^e \\ T_2^e \\ \vdots\\ T_{n_{en}}^e \end{bmatrix}

<\boldsymbol{N}^e> = \begin{bmatrix} N_1^e(x, y) & N_2^e(x, y) & \cdots & N_{n_{en}}^e(x, y) \end{bmatrix}

记

[K^e] = \int_{\Omega^e} [B^e]^T [D ^e] [B^e] d\Omega

\{f^e\} =-\int_{\Gamma_n^e} <{N^e}>^T \bar{q} d\Gamma + \int_{\Omega^e} <{N^e}>^T s d\Omega

则有限元方程为：

[K^e] \{d^e\} = \{f^e\} \tag{29}

至此，使用有限元方法求解二维稳态热传导基本方程的细节还没有完全介绍完，比如，单元形函数的推导、单元刚度直接组装法、边界条件的处理等。因此本文章会不断被更新。

[^1]: R.M. Barannon. Functional and Structured Tensor Analysis for Engineers(draft)[M]. 2003, P21.
[^2]: Jacob Fish, Ted Belytschko. A First Course in Finite Elements[M]. 2007, P142.
[^3]: Eduardo W. V. Chaves. Notes on Continuum Mechanics[M]. 2013, P328-P329.
[^4]: 上标T表示转置，除此之外，T代表温度
[^5]: 推导过程会用到格林公式，而格林公式的推导基于散度定理，散度定理的推导又基于格林定理
[^6]: 其实这里把\bar q写成\{q\}^{T} \{n\}会更好一些，因为选取的单元的任意性，并不一定包含给定通量\bar q的边界，这样也能方便理解总体有限元方程中的通量向量只有给定通量的边界有贡献