有限元入门 Part 2.

微信扫一扫,分享到朋友圈

有限元入门 Part 2.

不知不觉变成年番了hhh 去年的Part 1.总结了变分法和有限元之间的关系,主要偏重于理论推导,本文Part 2.将偏重于如何具体实现。有兴趣的同学可以去回顾一下Part 1.~


Dr.Stein:有限元入门 Part 1. zhuanlan.zhihu.com

从Part 1.的最后两个公式出发,这次我们通过一个简单的MATLAB的程序体验一下如何具体实现有限元法(好像上传不了附件,但是如果有需要可以私信我邮箱,我可以发给你matlab script)。

  • 旋转轴对称的热弹性问题

问题描述 :一个空心圆柱体,内表面和外边面的温度均为均匀且温度差为 ,假设热膨胀的参考温度为内表面的温度。边界条件为上下表面在z方向上的位移为0。求该空心圆柱体均匀受热膨胀后的位移场,应变场和应力场。


我将把问题的求解划分为以下几个重要的步骤并一一进行详细讲解。

(1)划分网格

之前的文章主要在说如何用变分法将PDE转化为弱形式,而很少提起划分网格这一步。但实际上划分网格是一门艺术。大部分同学对网格的直观感受是:网格越密结果越精确。大多数情况下这个结论是正确的。但更重要的是,网格的单元类型和疏密经常会影响求解的收敛性。

大部分的有限元商业软件都会有自己的网格划分功能。除此之外,还有一些专门的网格生成器,比如我比较常用的Gmsh: https:// gmsh.info/ 。对于这个例子中的简单的矩形求解域的二维问题,我们可以使用最常见的三角形单元,并且用MATLAB就可以直接生成网格。不论用什么样的方式生成网格,我们通常用两个矩阵来描述一个网格,即节点矩阵N和连接矩阵T。拿三角形单元网格举例,节点矩阵是一个N乘以2的矩阵,矩阵的每一行是每个节点的x,y坐标。很明显N行代表一共有N个节点。而T矩阵是一个M乘以3的矩阵,每一行对应每一个单元,并且每一行的三个数是组成该单元的三个节点的标号。

对于这个旋转轴对称问题,我们可以利用旋转轴对称的特性把该问题从一个三维问题转化为在截面上的二维问题。(我这边有一个含有节点矩阵和连接矩阵的网格文件,同样,有需要可以给我留邮箱)

(2)构建单元刚度矩阵和载荷向量

单元刚度矩阵可以通过Part 1.中的倒数第二个公式的等式左边的积分求得。该积分的被积函数由四阶的stiffness tensor 和形函数的偏导数组成。这里我们将 用Voigt notation表示,即降阶为二阶矩阵[ D ]:基于线弹性力学,应力和应变的本构关系可通过下面的[ D ]矩阵描述: (对于旋转轴对称问题,四个元素分别对应rr, zz, , rz四个方向)

单元内的位移场和节点位移可以用下面的形函数矩阵 来描述:

其中 为该三角形单元的三个节点上的形函数。并且这里我们使用一阶线性的形函数,即 。根据形函数的特性,即形函数在自身节点处的值为1,在其余节点处的值为0,我们不难构建出三个形函数的表达式为:

其中 A 为该三角形单元的面积:

节点位移向量为 的向量(三个节点,每个节点两个自由度):

通过节点位移向量和形函数矩阵可以计算出单元内的位移场[ u ]:

接下来,我们构建应变和位移的关系: ,转换成矩阵的形式: ,其中矩阵[ C ]为偏导矩阵:

合并之前的公式可得: ,其中[ B ]为:

根据上一篇文章Part 1.的倒数第二个公式的等式左边,我们可以构建单元刚度矩阵(6 x 6):

这里我们为了计算方便做了一定的简化。理论上[ B ]矩阵是坐标的函数,但是这里我们用该三角形单元的中心( )来计算[ B ],这样被积分项就成为了常数。同时单元内的体积分可以简化为 乘以被积分项(A为三角形单元的面积)。这个假设在单元很小的时候是可以近似接受的。但是真正的有限元程序一般是用等参单元和雅克比矩阵去精确计算这个积分。

接下来我们需要计算载荷向量,即Part 1.的倒数第二个公式的等式右边。对于这个热膨胀问题,我们需要首先构建热应变向量(只有正应变而没有切应变):

类似于构建单元刚度矩阵,我们通过下式构建出载荷向量,同样用了 进行简化:

(3)Assemble单元刚度矩阵

得到单元刚度矩阵和单元载荷向量之后,我们需要将所有单元组合成全局刚度矩阵和载荷向量从而最终求解所有节点上的位移。假设网格一共有N个节点,那么这个问题的全局刚度矩阵将为2N x 2N大小。在程序里,我们可以添加一个遍历所有单元的循环,对于第m个单元,将对应的6 x 6的单元刚度矩阵放至[node, node],其中node为 为第m个单元的三个节点的index,可从网格文件的T矩阵中提取。类似地,将第m个单元的载荷向量放至全局载荷向量的[node]处。这样,我们就得到了全局刚度矩阵 和全局载荷向量

(4)施加边界条件

接下来,我们需要施加边界条件。我们知道如果一个问题没有边界条件,那么是得不到唯一解的。对于这个问题,需要施加的边界条件是上下表面的 ,是一个Dirichlet边界条件。假设节点i处于上下表面,我们需要把全局刚度矩阵[K]里的第2i行第2i列的元素设为1(每个节点有两个自由度,第一个对应于 ,第二个对应于 ),并把第2i行和第2i列的其余所有元素均设为0。对于载荷向量[Q],将其第2i个元素替换为需要施加的节点位移(对于这个问题,这个值为0)。最后,对于[Q]中除了第2i个元素的其余所有元素进行变换:例如对于第j个元素( ), 。其中 为施加边界条件前的原刚度矩阵和载荷向量的元素, 为所施加的第2i个自由度的位移值。对于本问题, ,所以不需要对其余的Q元素进行变换。

(5)求解

对于这个简单线性问题,我们可以在MATLAB中使用 ’\‘ 直接进行求解得到节点位移向量。

但现实中的FEM Code一般使用基于LU分解的Direct Solver或者类似于Conjugate Gradient Method的Iterative Solver进行求解。

Solutions to Linear Systems of Equations: Direct and Iterative Solvers www.comsol.com

(6)基于位移的结果计算应变和应力

解得了节点位移向量之后,我们可以使用形函数进行插值从而得到求解域上任意点的位移,并且使用下式计算应变场和应力场:

(7)结果

这个问题是可以计算出解析解的,所以我们可以用解析解去验证FEM的计算结果。

FEM位移场结果:


位移场的FEM结果和解析解的比较:


FEM应变场结果:


应变场的FEM结果和解析解的比较:


FEM应力场结果:


应力场的FEM结果和解析解的比较:


最后重申下,如果想参照我的MATLAB script和网格文件的同学可以私信我邮箱,或者大家可以自己用习惯的软件编写试试,ENJOY!

微信扫一扫,分享到朋友圈

有限元入门 Part 2.

软件定义网络(SDN)第二次实验报告

上一篇

世界首座高铁跨海大桥主塔封顶:设计时速350公里 1小时内福州到厦门

下一篇

你也可能喜欢

有限元入门 Part 2.

长按储存图像,分享给朋友