FEM之理论通俗有限元(2-1)---偏微分方程应用(Matlab/COMSOL)

本文介绍分别用Matlab，COMSOL， 以及自己手工 求解拉普拉斯方程。

1.Matlab

命令行输入 

1. pdetool

2. 启动界面，工具栏选择画长方形工具，画一个正方形，如图1-1

3. 点击菜单 Options--》Application--》Heat transfer

4. 点击菜单 Boundary--》 Boundary Mode， 选中图中左边一条边，温度设置为100， 选中右边一条边温度设为0.

5. 点击菜单 Solve--Solve PDE 如图1-2

图1-1：

图1-2：

Matlab 解偏微分方程最大的局限性是 只能求解平面不能解三维。

2. COMSOL

1. Space 选择2D

2. 求解类型选择Heat transfer-》Heat transfer in solid，求解类型选Stationary

3. 选中树形菜单Geometry，右键，选择rectangle，然后画出一个正方形

4. 选中树形菜单Material，选中一个Build-in的材料

5. 选中树形菜单Heat Transfer in Solids,右键选中 Temperature，然后选中长方形左边设置100,

6. 重复5，选中正方形右边设置0

7. 点击树形菜单Study1，右键计算 Compute selected step

结果看起来和Matlab 不太一样，是因为材料属性热导率不同，Matlab 默认边界温度为0。

3. 手工求解

先将正方形离散为三角单元，为三角单元建立形函数，然后利用伽辽金方法推导出积分公式，加入边界，建立可编程实现的有限元方程，求解线性方程，计算出结果。

三角单元离散化参考：FEM之单元(1)---三角单元介绍 。

假设三角单元温度场函数为：

T = [S1 S2 S3] [Ti Tj Tk]；

S1 S2 S3为 线性形函数；

S1=(a1+b1x+c1y)/2A；

S2=(a2+b2x+c2y)/2A；

S3=(a3+b3x+c3y)/2A；

Ti Tj Tk 为 三角形三个顶点上的温度；

A 为三角形面积；

T 表示了三角单元内部任意一点的温度。

下面最主要的一步是利用伽辽金方法为三角单元建立方程：

将形函数与偏微分方程乘积的积分余差为0.
其中

即为[S1 S2 S3]形函数，其中包含了变量x，y，上式可以拆成两个

对第一个式子进行解析，按照分步积分和链式求导法则，有如下表达式：

将该式子右边第二项移动到左边，就得到了要计算的表达式

其中 

求导无任何问题，简单的高等数学知识；

表达式

计算要利用一下 格林函数，将面积分转化成曲线积分，其实也是大学高等数学知识。最终的计算结果可以将积分表达式转化成用a，b，c表现的矩阵形式。

上面一步是很多教科书都没有写明的，一带而过，但这才是用伽辽金方法生成有限元方程的关键所在。

形成矩阵形式以后，就可以很方便的用编程实现了。详细推导过程可以参考《有限元分析--ANSYS理论与应用》P310--P312， 作者 Saeed Moavenl，这个系列的书都不错，深入浅出，通俗易懂。

下一步就是要将单个矩阵组装成整体矩阵表达式。知道了节点编号，就是对号入座就行了。假设节点数目为N，生成一个N*N的总体矩阵K，按节点编号，将单个单元标号的矩阵放入整体矩阵中，不同单元相同节点号直接相加即可，这表明了单元的连续性。之后就是边界条件，将已知边界的温度生成一个N的向量B，解方程 Kx=B，即可得到温度场分布。在解决实际问题时，由于矩阵特别大，在总刚组装以及求解线性方程组时会非常讲究，以后有时间讨论。

注：文中所使用的软件仅供学习使用。