一维热传导的有限元求解基础与COMSOL弱形式实现

问题描述

如图所示，一根长度为 $L$ 的金属棒，左端有恒定热量 $q$ 流入，右端保持恒定温度 $T_L$，有电流流过金属，恒定产热量为 $Q$，热导率 $k$ 为常数求金属棒上的温度分布。

傅里叶定律

傅里叶定律，也称热传导定律，描述热传导过程。其微分形式如下：
$$\mathbf{q}=-k\nabla T$$
其中，
$\mathbf{q}$ 为局部热流密度，$W{m}^{-2}$
$k$ 为材料热导率，$W{m}^{-1}{K}^{-1}$
$\nabla T$ 为温度梯度，$K{m}^{-1}$

一维形式：
$$\mathbf{q}=-k\frac{dT}{dx}$$

As an aside
物理场中flux与force很有意思

Physical field	Flux	Force
Heat	$\mathbf{q}(W{m}^{-2})$	$-\nabla T$
Electrics	$\mathbf{I}(A{m}^{-2})$	$-\nabla \varphi$
Diffusion	$\mathbf{J}(mol{m}^{-2}{s}^{-1})$	$-\nabla c$

后两者分别是Ohm定律与Fick定律

将问题描述转化为数学语言，因存在内部热源 $Q(W{m}^{-3})$，故控制方程为
$$\nabla \cdot (-k\nabla T)=Q$$
一维形式为：
$$-k\frac{d{T}^2}{d^{2}x}=Q,0<x<L$$
边界条件为：
$$-k\frac{dT}{dx}=\mathbf{q},x=0$$$$T=T_L,x=L$$
求解方程即可得到温度分布

解析解

二阶常微分方程，可直接采用公式求解或者手动求解。
也可采用MATLAB的符号运算求解

% 创建符号变量，位置坐标x，内部单位产热量Q，材料热导率k，右端恒温TL，总长L，C1,C2
% 为积分常数, q为左端恒定流入热量
syms x C1 C2 Q k TL L q;

% 傅里叶热传导定律：dT^2=-Q/k*d^2x。进行积分
dT = int(-(Q/k),x)+C1;
T = int(dT,x)+C2;

% 代入边界条件 dT(0)=-q/k, T(L)=TL
% subs(T,x,0)表示将变量x替换为值0；solve(eqn,var),要求eqn右侧为0
s = solve(subs(dT,x,0)+q/k,subs(T,x,L)-TL,C1,C2);
s.C1;
s.C2;
T = subs(T,{C1,C2},{s.C1,s.C2})

可得
$$T(x)=T_L+\frac{\mathbf{q}}{k}(L-x)+\frac{Q}{2k}(L^2-x^2)$$

有限元求解思路

1. 单元划分

对与长为L的线段，我们可以将其划分为许多小单元，如下图：

每一单元可以表示为 e k = { x : x k ≤ x ≤ x k + 1 } e_k=\left\{x:x_k\leq x\leq x_{k+1}\right\} ek​={x:xk​≤x≤xk+1​}

假设每一段内的温度分布都是线性变化的，满足线性函数 ${\varphi }_k(x)$
则每一段内的温度分布可以近似为 $T(x_k)=a_k\cdot {\varphi }_k(x)$
$a_k$为系数，例如下图所示，每一段都是线性函数，但是其斜率是在变化的

总的温度分布可以找到一个近似表达式：
$$T(x)=\sum _{i=1}^{n+1}{a}_i\cdot {\varphi }_i(x)$$

分的越细，则该近似表达式越接近于解析解。

那么下一步的问题就是如何确定 ${\varphi }_i(x)$, 该函数称为形函数（shape function），即描述T的形状的函数。

2. 形函数的选择

注意温度的近似表达式需要满足原控制方程
$$-k\frac{d{T}^2}{d^{2}x}=Q$$
假如取${\varphi }_i(x)=bx+c$, 显然其二阶导并不存在！存在二阶导形式的函数将使得问题变得更加复杂。
由此，引进弱形式，即将原控制方程进行一次积分，变为一阶微分方程。

3. 弱形式表达式
等式两边同时乘上试函数(test function) $w_i(x)$，再进行一次积分，则：
$${\int }_0^L{w}_i(x)(-k\frac{d{T}^2}{d^{2}x}-Q)dx=0$$
在Galerkin处理中，试函数选择$w_i(x)={\varphi }_i(x)$

接下来，利用分步积分法消去二阶导项：
补充：$\int udv=uv-\int vdu$

那么上式变成
$${\int }_0^L{\varphi }_i(x)(-k\frac{d{T}^2}{d^{2}x})dx-{\int }_0^L{\varphi }_i(x)Qdx=$$$${\int }_0^{L}k\frac{d{\varphi }_i}{dx}\frac{dT}{dx}dx-{\varphi }_i(x)\cdot k\frac{dT}{dx}{|}_0^L-{\int }_0^L{\varphi }_i(x)Qdx=0$$

此即为弱形式表达式。可见在进行一次积分之后，消除了二阶项，但同时放宽了解要求，因此称之为weak form。

4. 试函数
为什么要引入试函数？直接两边积分不可以吗？
$${\int }_0^L{w}_i(x)(-k\frac{d{T}^2}{d^{2}x}-Q)dx=0$$
对于精确解，应当满足$-k\frac{d{T}^2}{d^{2}x}-Q=0$，同时也满足以上积分函数。
但是需要注意，我们在这里引入了$T(x)$的近似表达式，所以不可能完全满足，例如$-k\frac{d{T}^2}{d^{2}x}-Q=0.0001$，尽管存在误差（残差），但该误差是我们可以接受的。那么这个时候试函数将起到作用，即使得 $w_i(x)=0$，整个积分函数就是满足条件的。

5. 求解
将构造的近似表达式$T(x)$带入弱形式表达式中，最终求得每一单元内的温度分布。

即把$T(x)={\sum }_{i=1}^{n+1}{a}_i\cdot {\varphi }_i(x)$代入
$${\int }_0^{L}k\frac{d\varphi }{dx}\frac{dT}{dx}dx-\varphi (x)\cdot k\frac{dT}{dx}{|}_0^L-{\int }_0^L\varphi (x)Qdx=0$$

那么对于某一个单元$x_i\le x\le x_{i+1}$内
$$T(x)=a_i{\varphi }_i(x)+a_{i+1}{\varphi }_{i+1}(x)$$
选取试函数使得其满足（每一个节点的$a_i$值就是该节点的温度）
$$T(x_i)=a_i;T(x_{i+1})=a_{i+1}$$
令
$${\varphi }_i(x_i)=1;{\varphi }_i(x_{i+1})=0$$$${\varphi }_{i+1}(x_i)=0;{\varphi }_{i+1}(x_{i+1})=1$$
由此解得
$${\varphi }_i(x)=\frac{x_{i+1}-x}{x_{i+1}-x_i};{\varphi }_{i+1}(x)=\frac{x-x_i}{x_{i+1}-x_i}$$$$\frac{d{\varphi }_i}{dx}=-\frac{1}{x_{i+1}-x_i};\frac{d{\varphi }_{i+1}}{dx}=\frac{1}{x_{i+1}-x_i}$$
为了方便起见，只划分两个单元，$e_1(0,L/2),e_2(L/2,L)$
将$T(x)=a_1{\varphi }_1(x)+a_2{\varphi }_2(x)$到第一个单元内的弱表达式：
试函数为${\varphi }_1$时，
$$\sum _{j=1}^{2}k\left[{\int }_0^{L/2}\frac{d{\varphi }_1}{dx}\frac{d{\varphi }_j}{dx}\right]a_j-{\varphi }_1(x)\cdot k\frac{dT}{dx}{|}_0^{L/2}-{\int }_0^{L/2}{\varphi }_1(x)Qdx$$
试函数为${\varphi }_2$时，
$$\sum _{j=1}^{2}k\left[{\int }_0^{L/2}\frac{d{\varphi }_2}{dx}\frac{d{\varphi }_j}{dx}\right]a_j-{\varphi }_2(x)\cdot k\frac{dT}{dx}{|}_0^{L/2}-{\int }_0^{L/2}{\varphi }_2(x)Qdx$$
代入
$${\varphi }_1(x)=\frac{L/2-x}{L/2};{\varphi }_2(x)=\frac{x}{L/2}$$$$\frac{d{\varphi }_1}{dx}=-\frac{1}{L/2};\frac{d{\varphi }_2}{dx}=\frac{1}{L/2}$$
可以得到矩阵形式的线性方程：
$$\frac{2k}{L}\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]\left[\begin{array}{c}{a}_1\\ a_2\end{array}\right]-\frac{QL}{4}\left[\begin{array}{c}1\\ 1\end{array}\right]-\left[\begin{array}{c}q\\ 0\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right]$$
对第二个单元进行同样的操作，可得
$$\frac{2k}{L}\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]\left[\begin{array}{c}{a}_2\\ a_3\end{array}\right]-\frac{QL}{4}\left[\begin{array}{c}1\\ 1\end{array}\right]-\left[\begin{array}{c}0\\ 0\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right]$$
进行矩阵组装，可得
$$\frac{2k}{L}\left[\begin{array}{ccc}1& -1& 0\\ -1& 2& -1\\ 0& -1& 1\end{array}\right]\left[\begin{array}{c}{a}_1\\ a_2\\ a_3\end{array}\right]-\frac{QL}{4}\left[\begin{array}{c}1\\ 2\\ 1\end{array}\right]=\left[\begin{array}{c}q\\ 0\\ 0\end{array}\right]$$
最终解得$a_1,a_2,a_3$，也就是每一个节点的值

COMSOL求解

以下代入具体数值进行求解

Parameters	Value
$L(m)$	1
$k(W{m}^{-1}{K}^{-1})$	1
$Q(W{m}^{-3})$	20
$T_L(K)$	300
$q(W{m}^{-2})$	10

解析解为
$$T(x)=300+10(1-x)+10(1-x^2)$$
以下采用COMSOL求解数值解

传热模块

首先绘制线段，其次设定固体1

添加热源节点，并选中域

设定右侧边界条件，添加温度1

设定左侧边界条件，添加热通量1

随后进行求解可得

PDE模块

1. 系数形式偏微分方程
建议阅读COMSOL的博客加深理解
https://cn.comsol.com/blogs/brief-introduction-weak-form/
https://cn.comsol.com/blogs/implementing-the-weak-form-in-comsol-multiphysics/
https://cn.comsol.com/blogs/discretizing-the-weak-form-equations/
https://www.comsol.com/blogs/strength-weak-form/

选择系数形式偏微分方程，因变量设为T，单位K，源项单位为W/m^3

构建几何，将各参数值填入

右端选择狄利克雷边界条件

左端选择通量/源条件

稳态求解器，求解

2. 弱形式偏微分方程
已知弱形式为
$${\int }_0^{L}k\frac{d{\varphi }_i}{dx}\frac{dT}{dx}dx-{\varphi }_i(x)\cdot k\frac{dT}{dx}{|}_0^L-{\int }_0^L{\varphi }_i(x)Qdx=0$$
我们把非积分项移到右边
$${\int }_0^L(k\frac{d{\varphi }_i}{dx}\frac{dT}{dx}-{\varphi }_i(x)Q)dx={\varphi }_i(x)\cdot k\frac{dT}{dx}{|}_0^L$$
由此可以写出在COMSOL中的弱表达式，注意试函数用test表示

右侧非积分项写进弱贡献中，左端x=0时,
$$-{\varphi }_i(x)\cdot kq$$
因此

右侧x=L时，温度恒定，计算弱贡献时需引入拉格朗日乘子u

也可以直接采用狄利克雷边界条件
最后计算求解

Ref:
https://en.wikipedia.org/wiki/Thermal_conduction#Fourier’s_law
Book, The finite element method _ basic concepts and applications with MATLAB, MAPLE