问题描述

如图所示，一根长度为 LLL 的金属棒，左端有恒定热量 qqq 流入，右端保持恒定温度 TLT_LTL，有电流流过金属，恒定产热量为 QQQ，热导率 kkk 为常数求金属棒上的温度分布。

傅里叶定律

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

一维形式：
q=−kdTdx\boldsymbol{q}=-k\frac{dT}{dx} q=−kdxdT

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

Physical field	Flux	Force
Heat	q(Wm−2)\boldsymbol{q}(W\ m^{-2})q(W m−2)	−∇T-\nabla T−∇T
Electrics	I(Am−2)\boldsymbol{I}(A\ m^{-2})I(A m−2)	−∇ϕ-\nabla \phi−∇ϕ
Diffusion	J(molm−2s−1)\boldsymbol{J}(mol\ m^{-2}\ s^{-1})J(mol m−2 s−1)	−∇c-\nabla c−∇c

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

将问题描述转化为数学语言，因存在内部热源 Q(Wm−3)Q(W m^{-3})Q(Wm−3)，故控制方程为
∇⋅(−k∇T)=Q\nabla ·(-k\nabla T)=Q ∇⋅(−k∇T)=Q
一维形式为：
−kdT2d2x=Q,0<x<L-k\frac{dT^2}{d^2x}=Q, \ 0<x<L −kd2xdT2=Q, 0<x<L
边界条件为：
−kdTdx=q,x=0-k\frac{dT}{dx}=\boldsymbol{q}, \ x=0 −kdxdT=q, x=0T=TL,x=LT=T_L, \ x=L T=TL, 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)=TL+qk(L−x)+Q2k(L2−x2)T(x)=T_L+\frac{\boldsymbol{q}}{k}(L-x)+\frac{Q}{2k}(L^2-x^2) T(x)=TL+kq(L−x)+2kQ(L2−x2)

有限元求解思路

1. 单元划分

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

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

假设每一段内的温度分布都是线性变化的，满足线性函数 ϕk(x)\phi_k(x)ϕk(x)
则每一段内的温度分布可以近似为 T(xk)=ak⋅ϕk(x)T(x_k)=a_k·\phi_k(x)T(xk)=ak⋅ϕk(x)
aka_kak为系数，例如下图所示，每一段都是线性函数，但是其斜率是在变化的

总的温度分布可以找到一个近似表达式：
T(x)=∑i=1n+1ai⋅ϕi(x)T(x)=\sum_{i=1}^{n+1}a_i·\phi_i(x) T(x)=i=1∑n+1ai⋅ϕi(x)

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

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

2. 形函数的选择

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

3. 弱形式表达式
等式两边同时乘上试函数(test function) wi(x)w_i(x)wi(x)，再进行一次积分，则：
∫0Lwi(x)(−kdT2d2x−Q)dx=0\int_0^L w_i(x)(-k\frac{dT^2}{d^2x}-Q)dx=0 ∫0Lwi(x)(−kd2xdT2−Q)dx=0
在Galerkin处理中，试函数选择wi(x)=ϕi(x)w_i(x)=\phi_i(x)wi(x)=ϕi(x)

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

那么上式变成
∫0Lϕi(x)(−kdT2d2x)dx−∫0Lϕi(x)Qdx=\int_0^L \phi_i(x)(-k\frac{dT^2}{d^2x})dx-\int_0^L \phi_i(x)Qdx= ∫0Lϕi(x)(−kd2xdT2)dx−∫0Lϕi(x)Qdx=∫0LkdϕidxdTdxdx−ϕi(x)⋅kdTdx∣0L−∫0Lϕi(x)Qdx=0\int_0^L k\frac{d\phi_i}{dx}\frac{dT}{dx}dx-\phi_i(x) ·k \frac{dT}{dx}\bigg|_{0}^{L}-\int_0^L \phi_i(x)Qdx=0 ∫0LkdxdϕidxdTdx−ϕi(x)⋅kdxdT∣∣∣∣0L−∫0Lϕi(x)Qdx=0

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

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

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

即把T(x)=∑i=1n+1ai⋅ϕi(x)T(x)=\sum_{i=1}^{n+1}a_i·\phi_i(x)T(x)=∑i=1n+1ai⋅ϕi(x)代入
∫0LkdϕdxdTdxdx−ϕ(x)⋅kdTdx∣0L−∫0Lϕ(x)Qdx=0\int_0^L k\frac{d\phi}{dx}\frac{dT}{dx}dx-\phi(x) ·k \frac{dT}{dx}\bigg|_{0}^{L}-\int_0^L \phi(x)Qdx=0 ∫0LkdxdϕdxdTdx−ϕ(x)⋅kdxdT∣∣∣∣0L−∫0Lϕ(x)Qdx=0

那么对于某一个单元xi≤x≤xi+1x_i\leq x\leq x_{i+1}xi≤x≤xi+1内
T(x)=aiϕi(x)+ai+1ϕi+1(x)T(x)=a_i\phi_i(x)+a_{i+1}\phi_{i+1}(x) T(x)=aiϕi(x)+ai+1ϕi+1(x)
选取试函数使得其满足（每一个节点的aia_iai值就是该节点的温度）
T(xi)=ai;T(xi+1)=ai+1T(x_i)=a_i; T(x_{i+1})=a_{i+1} T(xi)=ai;T(xi+1)=ai+1
令
ϕi(xi)=1;ϕi(xi+1)=0\phi_i(x_i)=1; \phi_i(x_{i+1})=0 ϕi(xi)=1;ϕi(xi+1)=0ϕi+1(xi)=0;ϕi+1(xi+1)=1\phi_{i+1}(x_i)=0; \phi_{i+1}(x_{i+1})=1 ϕi+1(xi)=0;ϕi+1(xi+1)=1
由此解得
ϕi(x)=xi+1−xxi+1−xi;ϕi+1(x)=x−xixi+1−xi\phi_i(x)=\frac{x_{i+1}-x}{x_{i+1}-x_i}; \phi_{i+1}(x)=\frac{x-x_i}{x_{i+1}-x_i} ϕi(x)=xi+1−xixi+1−x;ϕi+1(x)=xi+1−xix−xidϕidx=−1xi+1−xi;dϕi+1dx=1xi+1−xi\frac{d\phi_i}{dx}=-\frac{1}{x_{i+1}-x_i};\frac{d\phi_{i+1}}{dx}=\frac{1}{x_{i+1}-x_i} dxdϕi=−xi+1−xi1;dxdϕi+1=xi+1−xi1
为了方便起见，只划分两个单元，e1(0,L/2),e2(L/2,L)e_1(0,L/2), e_2(L/2, L)e1(0,L/2),e2(L/2,L)
将T(x)=a1ϕ1(x)+a2ϕ2(x)T(x)=a_1\phi_1(x)+a_{2}\phi_{2}(x)T(x)=a1ϕ1(x)+a2ϕ2(x)到第一个单元内的弱表达式：
试函数为ϕ1\phi_1ϕ1时，
∑j=12k[∫0L/2dϕ1dxdϕjdx]aj−ϕ1(x)⋅kdTdx∣0L/2−∫0L/2ϕ1(x)Qdx\sum_{j=1}^2k\left[\int_0^{L/2} \frac{d\phi_1}{dx}\frac{d\phi_j}{dx}\right]a_j-\phi_1(x) ·k \frac{dT}{dx}\bigg|_{0}^{L/2}-\int_0^{L/2} \phi_1(x)Qdx j=1∑2k[∫0L/2dxdϕ1dxdϕj]aj−ϕ1(x)⋅kdxdT∣∣∣∣0L/2−∫0L/2ϕ1(x)Qdx
试函数为ϕ2\phi_2ϕ2时，
∑j=12k[∫0L/2dϕ2dxdϕjdx]aj−ϕ2(x)⋅kdTdx∣0L/2−∫0L/2ϕ2(x)Qdx\sum_{j=1}^2k\left[\int_0^{L/2} \frac{d\phi_2}{dx}\frac{d\phi_j}{dx}\right]a_j-\phi_2(x) ·k \frac{dT}{dx}\bigg|_{0}^{L/2}-\int_0^{L/2} \phi_2(x)Qdx j=1∑2k[∫0L/2dxdϕ2dxdϕj]aj−ϕ2(x)⋅kdxdT∣∣∣∣0L/2−∫0L/2ϕ2(x)Qdx
代入
ϕ1(x)=L/2−xL/2;ϕ2(x)=xL/2\phi_1(x)=\frac{L/2-x}{L/2}; \phi_{2}(x)=\frac{x}{L/2} ϕ1(x)=L/2L/2−x;ϕ2(x)=L/2xdϕ1dx=−1L/2;dϕ2dx=1L/2\frac{d\phi_1}{dx}=-\frac{1}{L/2};\frac{d\phi_2}{dx}=\frac{1}{L/2} dxdϕ1=−L/21;dxdϕ2=L/21
可以得到矩阵形式的线性方程：
2kL[1−1−11][a1a2]−QL4[11]−[q0]=[00]\frac{2k}{L}\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}\begin{bmatrix} a_1 \\ a_2 \end{bmatrix}-\frac{QL}{4}\begin{bmatrix} 1 \\ 1 \end{bmatrix}-\begin{bmatrix} q \\ 0 \end{bmatrix}=\begin{bmatrix} 0 \\ 0 \end{bmatrix} L2k[1−1−11][a1a2]−4QL[11]−[q0]=[00]
对第二个单元进行同样的操作，可得
2kL[1−1−11][a2a3]−QL4[11]−[00]=[00]\frac{2k}{L}\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}\begin{bmatrix} a_2 \\ a_3 \end{bmatrix}-\frac{QL}{4}\begin{bmatrix} 1 \\ 1 \end{bmatrix}-\begin{bmatrix} 0 \\ 0 \end{bmatrix}=\begin{bmatrix} 0 \\ 0 \end{bmatrix} L2k[1−1−11][a2a3]−4QL[11]−[00]=[00]
进行矩阵组装，可得
2kL[1−10−12−10−11][a1a2a3]−QL4[121]=[q00]\frac{2k}{L}\begin{bmatrix} 1 & -1 & 0 \\ -1 & 2 & -1\\ 0 & -1 &1 \end{bmatrix}\begin{bmatrix} a_1 \\ a_2 \\ a_3 \end{bmatrix}-\frac{QL}{4}\begin{bmatrix} 1 \\2\\ 1 \end{bmatrix}=\begin{bmatrix} q \\ 0\\0 \end{bmatrix} L2k⎣⎡1−10−12−10−11⎦⎤⎣⎡a1a2a3⎦⎤−4QL⎣⎡121⎦⎤=⎣⎡q00⎦⎤
最终解得a1,a2,a3a_1,a_2,a_3a1,a2,a3，也就是每一个节点的值

COMSOL求解

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

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

解析解为
T(x)=300+10(1−x)+10(1−x2)T(x)=300+10(1-x)+10(1-x^2)T(x)=300+10(1−x)+10(1−x2)
以下采用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. 弱形式偏微分方程
已知弱形式为
∫0LkdϕidxdTdxdx−ϕi(x)⋅kdTdx∣0L−∫0Lϕi(x)Qdx=0\int_0^L k\frac{d\phi_i}{dx}\frac{dT}{dx}dx-\phi_i(x) ·k \frac{dT}{dx}\bigg|_{0}^{L}-\int_0^L \phi_i(x)Qdx=0 ∫0LkdxdϕidxdTdx−ϕi(x)⋅kdxdT∣∣∣∣0L−∫0Lϕi(x)Qdx=0
我们把非积分项移到右边
∫0L(kdϕidxdTdx−ϕi(x)Q)dx=ϕi(x)⋅kdTdx∣0L\int_0^L (k\frac{d\phi_i}{dx}\frac{dT}{dx}-\phi_i(x)Q)dx=\phi_i(x) ·k \frac{dT}{dx}\bigg|_{0}^{L} ∫0L(kdxdϕidxdT−ϕi(x)Q)dx=ϕi(x)⋅kdxdT∣∣∣∣0L
由此可以写出在COMSOL中的弱表达式，注意试函数用test表示

右侧非积分项写进弱贡献中，左端x=0时,
−ϕi(x)⋅kq-\phi_i(x) ·k q −ϕi(x)⋅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

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

1. comsol 中的弱形式问题——一维稳态导热
弱形式推导的一般方法写出问题的偏微分方程乘以试函数并对方程进行积分采用分部积分进行微分降级得到方程的弱形式一维稳态导热微分方程的弱形式 1.问题描述如图所示,一根长度为 L=1m的金属棒,左 ...
COMSOL弱形式解微分方程
COMSOL中建模在COMSOL中,根据模型向导,选择一维稳态弱形式偏微分方程(w)(数学->偏微分方程接口),因变量可以设置为T. 对于几何模型,通过线段间隔命令建立一段线段,左端点为1, ...
利用有限元法求解一维热传导问题
我们都知道热传导方程是大名鼎鼎的傅里叶提出的,而对于此方程的解也可以用傅里叶变换的方法得到,但笔者已经记不起该怎么求解了(尽管上学期才刚学过数理方法啊).但是无伤大雅,因为我已经掌握了(其实还没有)对 ...
python有限元传热求解_二维稳态热传导基本方程的有限元求解(2)
四节点矩形单元在二维稳态热传导基本方程的有限元求解(1)这篇文章中,我们仅仅给出了有限元单元方程的一种比较标准的推导步骤,并未涉及某种具体的单元.且在式(20)中,单元上温度的近似函数表示成节点 ...
有限元_一维拉杆问题的求解方法（解析法、差分法、试函数法）
一维拉杆问题的求解方法(解析法.差分法.试函数法) 为改善求解精度,就可以取多个基底函数.多个待定系数.
Python——求解一维热传导
1.问题求一维热传导方程混合问题: 求解其的数值解,取N=10,h=0.1,计算到K=36为止. 2.程序 #python import numpy as np import matplotlib. ...
【CAE】优秀的开源有限元求解器
工业有限元求解器是工业仿真中最重要,最核心的,目前市场上有较多的开源求解器,下面列举几个比较有名的开源求解器,仅供各位对国产CAE感兴趣的朋友参考.借鉴,学习CAE软件开发的框架,思路等. 1. Op ...
一维有限元法matlab,有限元matlab研究.ppt
* 有限元方法--介绍有限元方法是数值求解偏微分方程边值问题的一种方法,此方法首先于20世纪50年代初由工程师提出,并用于求解简单的结构问题.有限元方法是这一种系统的数值方法,并奠定其数学基础,是在 ...
MATLAB模拟导热过程,一维热传导MATLAB模拟.doc
PAGE 昆明学院 2015 届毕业设计(论文) 设计(论文)题目一维热传导问题的数值解法及其MATLAB模拟子课题题目无姓名伍有超学号 201117030225 所属系物 ...
热传递 matlab,一维热传导MATLAB模拟.pdf
昆明学院 2015 届毕业设计(论文) 设计(论文)题目一维热传导问题的数值解法及其 MATLAB 模拟子课题题目无姓名伍有超学号 201117030225 所属系物理科学 ...

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

目录