温度场有限元计算的研究（1）

研究温度场计算的相关理论知识，主要是通过研究电器产品的数学建模过程，得到PDE和它的边界条件，即构成一个定解问题，再将该定解问题通过变分原理或者Galerkin方法得到一个离散的积分，最后离散化得到代数方程组。

传热学入门

主要是掌握基本方程和几种边界条件。

基本概念

热能传递有三种基本方式：热传导，热对流，热辐射。

热传导指的是物体各部分之间不发生相对位移时，依靠分子、原子及自由电子等微观粒子的热运动而产生的热能传递，遵循的是傅里叶定律。热对流是指由于流体的宏观运动而引起的流体各部分之间发生相对的位移，冷、热流体相互掺混所导致的热量传递过程，仅能发生在流体当中。对流传热的基本计算式是牛顿冷却公式。物体通过电磁波来传递能量的方式称为热辐射。

针对电器产品，我们首要考虑的是电器产品的热传导。

热传导方程及其边界条件

热传导的基本定律是傅里叶定律，热传导的基本方程和边界条件如下:

$\left\{ \begin{aligned} &\Omega:\nabla\cdot\lambda\nabla T = -G\\ &\Gamma_1:T=T_g\\ &\Gamma_2:\lambda\dfrac{\partial T}{\partial n} = q\\ &\Gamma_3:\lambda\dfrac{\partial T}{\partial n} = h(T-T_0) \end{aligned} \right.$

其中， $\lambda$ 是导热系数( $W/(m\cdot k)$ )， $G$ 是内热源强度( $W/m^3$ )， $q$ 是热流密度( $W/m^2$ )， $h$ 是对流换热系数( $W/(m^2\cdot k)$ )， $T_0$ 是周围介质温度( $k$ )。

第一类边界条件：边界上的温度或温度分布函数已知。

第二类边界条件：边界上的热流密度（热通量） $q$ 已知（热传导边界条件）。

第三类边界条件：边界与周围环境介质间的对流传热系数 $h$ 和介质温度 $T_0$ 已知（对流边界条件）。

有限元离散

为了简化计算，工程师提出了将温度场等效为热路的手段，这种方法虽然简单，但是不精确。想要准确的计算空间中的温度分布还是要求解场的PDE。PDE的解析解基本上不可能求出，所以需要介绍将PDE通过有限元方法离散成代数方程组，也就是形如 $[S][A] = [F]$ 的矩阵形式的方法。

有限元的基本原理是，通过分网将空间划分为成百上千个单元(比如说二维的三角形单元，三维的四面体单元)，每一个单元都满足同一个PDE。对每一个单元单独分析，即在该单元的范围内将PDE离散为代数方程组，再通过装配得到整个求解域内的代数方程组，最终求解装配后的大型代数方程组得到每一个节点的待求物理量的值。

尽管变分过程或加权余量法对于有限元离散是必要的，但是这个过程并不是通过编程来实现，而且对数学差的工科生（正是在下）而言，这种级别的微积分过程确实有难度。因此，对新手而言，直接根据PDE的形式，找到对应的有限元离散格式，编程计算也就够了。

变分过程

太麻烦了，要学泛函，建议不学这玩意，直接参考加权余量法中的Galerkin方法。反正笔者没学明白（逃。

事实上，通过Galerkin方法离散得到的代数方程组和通过变分得到的代数方程组是一样的。

加权余量法

参见'The finite element method in engineering' P187 \textbf{5.9 WEIGHTED RESIDUAL APPROACH}

加权余量法和变分法一样，都是离散化PDE的手段，但是使用加权余量法是不需要了解变分过程的，应该说这是一种变分的替代手段，这里主要介绍伽辽金法。如果学过数值分析的话，可以感觉到，加权余量法和最小二乘法原理上有相通之处。

一般意义下的加权余量法

$u$ 是空间中待求的物理量，该物理量在空间中分布的算子方程为 $D(u)=T(u)-f=0$ 。假设在空间中有 $n$ 个固定点，这些固定点的坐标已知，它们的 $u$ 是未知的，那么，可以认为只要求解出这 $n$ 个固定点的 $u$ ，那么空间中任意位置的 $u$ ，都可以通过这 $n$ 个固定点的 $u$ 插值得到。

插值完成后，空间中任意一点的待求量 $\hat{u}=\sum_{j=1}^{n}u_jN_j$ ，其中 $u_j$ 就是已知的每个节点的待求量； $N_j$ 是插值后的形函数，只和位置有关； $\hat{u}$ 被称为试探解。

尽管完成了插值，但是等号右侧的每个固定点的 $u_j$ 是未知的。我们无法保证每一点的试探解 $\hat{u}$ 都能够满足算子方程，除非 $\hat{u}$ 是真解，但是我们希望 $\hat{u}$ 和 $u_j$ 能够尽量接近真解。为了求解出 $u_j$ ，我们提出一种策略，那就是将空间中所有通过插值得到的试探解 $\hat{u}$ 代入到算子方程 $T(u)-f$ 当，如果它们的和能够满足算子方程，那么我们认为固定点处的 $u_j$ 是接近于真实的 $u$ 的。很明显，这是对整个求解域的积分过程。

$\int_{\Omega}(T(\hat{u})-f)d\Omega=0$

$\int_{\Omega}\sum_{j=1}^{n}(u_jT(N_j)-f)d\Omega=0$

但是这个策略本身是有问题的，因为这一组积分只能得到一组关于 $u_j$ 的关系式，也就是说，对于n个未知节点而言需要n组关系式来求解。如果能够定义n个权函数，保证上述积分在加权的情况下都能够满足要求，我们认为求出的 $u_j$ 接近真实解。

$\int_{\Omega}\sum_{j=1}^{n}W_i(u_jT(N_j)-f)d\Omega=0,\:i=1,2,...,n$

伽辽金方法

权函数有多种选取方法，Galerkin方法是应用最广泛也是最容易记住的一种方法。Galerkin方法直接取插值后的形函数 $N_i$ 作为权函数。

$\int_{\Omega}\sum_{j=1}^{n}N_i(u_jT(N_j)-f)d\Omega=0,\:i=1,2,...,n$

加权余量法在有限元当中的应用

之前提到的加权余量法，都是建立在已知空间中 $N$ 个节点的情况下，但是有限元方法是基于单元分析-装配-求解这个过程的，加权余量法仅仅是应用于单元分析当中，将PDE离散成代数方程组，然后装配。

我们可以思考一个问题，假设存在一个二维模型，三角形分网后存在 $n$ 个已知坐标的节点。一般理解上，加权余量法应该直接对这 $n$ 个节点插值，加权积分然后离散，但是事实上在有限元方法中，加权余量法是对每个单元插值，加权积分然后离散，最后装配。这个过程为什么是成立的？也就是说，有限元分析中的加权余量法，其装配的原理是什么？

为此笔者给出自己的理解。对于加权余量法而言，权函数的选取并没有固定的要求，只要满足线性独立的条件既是可以的。对于 $N$ 个节点而言，重要的是获取 $N$ 组线性独立的权函数，至于整个大型方程组中的权函数的形式并不是一个值得关心的问题。将每个单元通过Galerkin方法离散，使得每个单元的误差满足Galerkin方法的要求，并且将单元中的权函数按照节点对应的位置进行装配后，得到新的 $N$ 组权函数，这就够了。此时的大型线性方程组应该是不满足通常意义上对n个节点插值的Galerkin方法的。

温度场有限元中的伽辽金方法

温度场求解域内的PDE为 $\Omega:\nabla\cdot\lambda\nabla T+G=0$ ，设 $\hat{T}=\sum_{j=1}^{n}N_jT_j$ ，通过Galerkin方法积分得到：

$\tag{1} \int_{\Omega}N_i(\nabla\cdot\lambda\nabla\hat{T}+G)d\Omega=0,\:i=1,2,...,n$

由Green变换（不会，但是这个变换在Laplace方程和Poission方程的有限元求解中都能用上，建议记住）得到：
$\tag{2} \int_{\Omega}N_i(\nabla\cdot\lambda\nabla\hat{T})d\Omega= -\int_{\Omega}\lambda\nabla N_i\cdot \nabla\hat{T}d\Omega+ \int_{\Gamma}N_i\lambda\dfrac{\partial\hat{T}}{\partial n}dS$

将(2)代入(1)中，得到：
$\tag{3} \int_{\Omega}\lambda\nabla N_i\cdot \nabla\hat{T}d\Omega- \int_{\Omega}N_iGd\Omega- \int_{\Gamma}N_i\lambda\dfrac{\partial\hat{T}}{\partial n}dS=0$

分析(3)中的第三项，即对边界积分的项，根据三种不同的边界条件分类讨论：

第一类边界条件：在强迫边界条件下这一项是未知的，不能在单元内求解，所以直接忽略该项，装配完成以后在求解过程中直接给该节点赋值，计算出其它点的结果即可。

第二类边界条件： $\lambda\dfrac{\partial\hat{T}}{\partial n}=q_w$ ，
$\int_{\Gamma_2}N_i\lambda\dfrac{\partial\hat{T}}{\partial n}dS= \int_{\Gamma_2}N_iq_wdS$

第三类边界条件： $\lambda\dfrac{\partial\hat{T}}{\partial n}=k(T-T_f)$ ，
$\int_{\Gamma_3}N_i\lambda\dfrac{\partial\hat{T}}{\partial n}dS= \int_{\Gamma_3}N_ik(\hat{T}-T_f)dS$

将 $\hat{T}=\sum_{j=1}^{n}N_jT_j$ 代入(3)：

$\int_{\Omega}\lambda\nabla N_i\cdot \nabla\hat{T}d\Omega=\sum_{j=1}^{n}T_j\int_{\Omega}\lambda\nabla Ni\cdot\nabla N_jd\Omega,\:i = 1,2,...,n$

$\int{\Gamma_3}N_ik\hat{T}dS= \sum_{j=1}^{n}T_j\int_{\Gamma_3}N_ikN_jdS,\:i=1,2,...,n$

综上可得：
$\sum_{j=1}^{n}T_j\int_{\Omega}\lambda\nabla Ni\cdot\nabla N_jd\Omega+ \sum_{j=1}^{n}T_j\int_{\Gamma_3}kN_iN_jdS= \int_{\Omega}N_iGd\Omega+\int_{\Gamma_2}N_iq_wdS,\:i=1,2,...,n$

离散为代数方程组

考虑三种模型的代数方程组离散：二维平面温度场，二维轴对称温度场，三维温度场。这里给出的是上面的方程求解后的结论，需要编程求解的话直接参考这里即可。

一些必须的几何参数

参考颜威利《电气工程电磁场数值分析》。

二维场中的三角形单元：假设三个节点分别编号为1,2,3，坐标分别为 $(x_1,y_1)$ , $(x_2,y_2)$ , $(x_3,y_3)$ 。

$p_1=x_2y_3-y_2x_3$ ; $p_2=x_3y_1-y_3x_1$ ; $p_3=x_1y_2-y_1x_2$ ;

$q_1=y_2-y_3$ ; $q_2=y_3-y_1$ ; $q_3=y_1-y_2$ ;

$r_1=x_3-x_2$ ; $r_2=x_1-x_3$ ; $r_3=x_2-x_1$ ;

$\bigtriangleup=\dfrac{1}{2}(q_1r_2-q_2r_1)$ —— 三角形面积。

二维场中的线单元：假设两个节点分别编号为1,2，坐标分别为 $(x_1,y_1)$ , $(x_2,y_2)$ 。

$d_l=\sqrt{(x_2-x_1)^2+(y_2-y_1)^2}$ —— 线单元长度

三维场中的四面体单元：假设四个节点分别编号为1,2,3,4，坐标分别为 $(x_1,y_1)$ , $(x_2,y_2)$ , $(x_3,y_3)$ , $(x_4,y_4)$ 。

$p_1 = \left|\begin{array}{ccc} x_2 & y_2 & z_2\\ x_3 & y_3 & z_3\\ x_4 & y_4 & z_4 \end{array}\right|$ , $p_2 = -\left|\begin{array}{ccc} x_1 & y_1 & z_1\\ x_3 & y_3 & z_3\\ x_4 & y_4 & z_4 \end{array}\right|$ , $p_3 = \left|\begin{array}{ccc} x_1 & y_1 & z_1\\ x_2 & y_2 & z_2\\ x_4 & y_4 & z_4 \end{array}\right|$ , $p_4 = -\left|\begin{array}{ccc} x_1 & y_1 & z_1\\ x_2 & y_2 & z_2\\ x_3 & y_3 & z_3 \end{array}\right|$

$q_1 = -\left|\begin{array}{ccc} 1 & y_2 & z_2\\ 1 & y_3 & z_3\\ 1 & y_4 & z_4 \end{array}\right|$ , $q_2 = \left|\begin{array}{ccc} 1 & y_1 & z_1\\ 1 & y_3 & z_3\\ 1 & y_4 & z_4 \end{array}\right|$ , $q_3 = -\left|\begin{array}{ccc} 1 & y_1 & z_1\\ 1 & y_2 & z_2\\ 1 & y_4 & z_4 \end{array}\right|$ , $q_4 = \left|\begin{array}{ccc} 1 & y_1 & z_1\\ 1 & y_2 & z_2\\ 1 & y_3 & z_3 \end{array}\right|$

$r_1 = -\left|\begin{array}{ccc} x_2 & 1 & z_2\\ x_3 & 1 & z_3\\ x_4 & 1 & z_4 \end{array}\right|$ , $r_2 = \left|\begin{array}{ccc} x_1 & 1 & z_1\\ x_3 & 1 & z_3\\ x_4 & 1 & z_4 \end{array}\right|$ , $r_3 = -\left|\begin{array}{ccc} x_1 & 1 & z_1\\ x_2 & 1 & z_2\\ x_4 & 1 & z_4 \end{array}\right|$ , $r_4 = \left|\begin{array}{ccc} x_1 & 1 & z_1\\ x_2 & 1 & z_2\\ x_3 & 1 & z_3 \end{array}\right|$

$s_1 = -\left|\begin{array}{ccc} x_2 & y_2 & 1\\ x_3 & y_3 & 1\\ x_4 & y_4 & 1 \end{array}\right|$ , $s_2 = \left|\begin{array}{ccc} x_1 & y_1 & 1\\ x_3 & y_3 & 1\\ x_4 & y_4 & 1 \end{array}\right|$ , $s_3 = -\left|\begin{array}{ccc} x_1 & y_1 & 1\\ x_2 & y_2 & 1\\ x_4 & y_4 & 1 \end{array}\right|$ , $s_4 = \left|\begin{array}{ccc} x_1 & y_1 & 1\\ x_2 & y_2 & 1\\ x_3 & y_3 & 1 \end{array}\right|$

$Ve=\dfrac{1}{6}\left|\begin{array}{cccc} 1 & x_1 & y_1 & z_1\\ 1 & x_2 & y_2 & z_2\\ 1 & x_3 & y_3 & z_3\\ 1 & x_4 & y_4 & z_4 \end{array}\right|$ ——四面体体积

三维场中的三角形单元：和二维场的三角形单元差不多，需要用到的参数只有一个 $\bigtriangleup$ 。 $\bigtriangleup=\dfrac{1}{2}\left|\vec{AB}\times\vec{AC}\right|$ .

二维平面温度场

温度场求解和磁场求解有个很大的区别就是边界条件的处理。对于磁场来说，基本上只要考虑第一类边界条件，第二类边界条件为0，但是温度场的第二类边界条件和第三类边界条件是普遍存在的。为了处理边界条件，区域内使用三角形单元，边界用线单元，两种类型的单元装配成整个系统进行计算。

三角形单元： $S_{ij}^e = \dfrac{\lambda}{4\bigtriangleup}(r_ir_j+q_iq_j)$ ， $F_i^e = \dfrac{1}{3}G_e\bigtriangleup$

线单元：第二类边界条件： $S_{ij}^l = 0$ ， $F_i^l = \dfrac{1}{2}q_ld_l$

第三类边界条件： $S_{ij}^l = \dfrac{hd_l}{6}(1+\delta_{ij}^l)$ ， $\delta_{ij}^l = \left\{ \begin{aligned} 1 \quad j=i\\ 0 \quad j\neq i\\ \end{aligned} \right.$ ， $F_i^l = \dfrac{1}{2}hT_0d_l$

二维轴对称温度场

关于轴对称温度场的离散方法很多，不同的离散方法会导致些许误差。笔者最初参考的教材是'The Finite Element Method in Engineering'，求出来的结果根本不对。后来检查发现这本书关于第二类边界条件和第三类边界条件离散方法存在一些问题。按照这本书的离散方法，三角形单元和线单元分别求解之后，两种单元的系数矩阵量纲是不同的，右侧向量的量纲也不同，这样的话装配得到的结果自然是错误的。后来找到一篇很古老的文章，《轴对称热传导有限元格式》，该文批判了当时的各种方法（90年代），然后提出了一种自己的计算方法，按照这篇文章给出的全新方法测试，发现和COMSOL的计算结果误差还是挺大的，而且很麻烦，所以就参考了这篇文章引用的文献[1]和文献[2]的方法。

三角形单元： $S_{ij}^e = \dfrac{r\lambda}{4\bigtriangleup}(r_ir_j+q_iq_j)$ ， $F_i^e =\dfrac{\pi G_e\bigtriangleup}{6}(3r+r_i)$

其中， $r=\dfrac{r_i+r_j+r_k}{3}$ .

线单元：第二类边界条件：& $S_{ij}^l=0$ ， $F_i^l=\dfrac{\pi qd_l}{2r'+2r_i}$

第三类边界条件： $S_{ij}^l=\dfrac{\pi hd_l}{6}(2r'+2r_i)$ ， $F_i^l=\dfrac{\pi hT_0d_l}{3}(2r'+r_i)$

其中 $r'=\dfrac{r_i+r_j}{2}$

三维温度场

三维温度场的求解域内使用四面体单元，边界使用三角形单元。

四面体单元： $S_{ij}^e=\dfrac{\lambda}{36V^e}(q_iq_j+r_ir_j+s_is_j)$ ， $F_i^e=\dfrac{G_eV^e}{4}$

三角形单元：第二类边界条件：& $S_{ij}^e=0$ ， $F_i^e=\dfrac{q\bigtriangleup}{3}$

第三类边界条件：& $S_{ij}^e=\dfrac{h\bigtriangleup}{12}(1+\delta_{ij})$ ， $\delta_{ij}^l = \left\{\begin{aligned} 1 \quad j=i\\ 0 \quad j\neq i\\ \end{aligned} \right.$ ， $F_i^e=\dfrac{hT_0\bigtriangleup}{3}$

COMSOL分网数据的解析

求解之前需要先建模和分网，COMSOL的分网数据含义是很清晰的，它由如下几部分构成：

（1）存在的节点数目，以及各个节点的坐标；

（2）存在几种类型的单元；

（3）各种单元包含了几个节点（四面体单元4个，三角形单元3个，线单元2个，顶点单元1个），所有该单元包括的节点坐标编号，以及每个单元所在的求解域。

值得注意的是这里的节点坐标编号和求解域编号是从0开始的，求解域编号为COMSOL建模后生成的求解域编号-1，所以想要提取某个边界单元的话直接利用分网信息中标出的域即可。

最后编辑于：2019.03.14 22:22:32

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 215,384评论 6赞 497
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 91,845评论 3赞 391
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 161,148评论 0赞 351
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 57,640评论 1赞 290
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 66,731评论 6赞 388
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 50,712评论 1赞 294
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 39,703评论 3赞 415
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 38,473评论 0赞 270
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 44,915评论 1赞 307
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 37,227评论 2赞 331
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 39,384评论 1赞 345
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 35,063评论 5赞 340
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 40,706评论 3赞 324
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 31,302评论 0赞 21
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 32,531评论 1赞 268
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 47,321评论 2赞 368
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 44,248评论 2赞 352