有限单元法基本思想-原理-数值计算过程

有限单元法学习报告在对力学问题分析求解过程中，方法可以概括为两种方法，一种为解析法，对具体问题具体分析，通过一定的推导用具体的表达式获得解答，由于实际工程中结构物的复杂性，此方法在处理工程问题是十分困难的；另一种是数值法，有限元法是其中一种方法，其数学逻辑严谨，物理概念清晰，又采用矩阵形式表达基本公式，便于计算机编程，因此在工程问题中获得广泛的应用。

有限元法基本原理是，将复杂的连续体划分为简单的单元体；将无限自由度问题化为有限自由度问题，因为单元体个数是有限的；将偏微分方程求解问题化为有限个代数方程组的求解问题。

通常以位移为基本未知量，通过虚功原理和最小势能原理来求解。

基本思想是先化整为零，即离散化整体结构，把整体结构看作是由若干个通过结点相连的单元体组成的整体；再积零为整，通过结点的平衡来建立代数方程组，最后计算出结果。

我将采用最简单的三结点三角形为基本单元体，解决弹性力学中的平面问题为例，解释有限单元法的基本原理、演示数值计算过程和一般性应用结论。

一、离散化解决平面问题时，主要单元类型包括三角形单元（三结点、六结点）和四边形单元（四结点矩形、四结点四边形、八结点四边形）等。

选用不同的单元会有不同的精度，划分的单元数越多，精度越高，但计算量也会越大。

因此在边界曲折，应力集中处单元的尺寸要小些，但最大与最小单元的尺寸倍数不宜过大。

在集中力作用点及分布力突变的点宜选为结点，不同厚度，不同材料不能划分在同一单元中。

三角形单元以内角接近60°为最好。

充分利用对称性与反对称性。

二、单元分析将一个单元上的所有未知量用结点位移表示，并将分布在单元上的外力等效到结点上。

1、位移函数选取：根据有限元法的基本思路，将连续体离散为有限的单元集合后，此时单元体满足连续性、均匀性、各向同性、完全线弹性假设。

单元与单元之间通过结点连接并传递力，位移法（应用最广）以结点位移δi=(u i v i)T为基本未知量，以离散位移场代替连续位移场。

单元体内的位移变化可以用位移函数（位移模式）来表示，因为有限元分析所得结果是近似结果，为了保证计算精度和收敛性，x位移函数应尽可能反应物体中的真实位移，即满足完备性和连续性的要求：①位移模式必须能反映单元的刚体位移。

②位移模式必须能反映单元的常量应变。

③位移模式应尽可能反应位移的连续性。

设三角形单元三个结点编号为i、j、m。

平面三角形单元位移函数选取为u=α1+α2x+α3y v=α4+α5x+α6y可以写成00u u yv v y ωω=-⎧⎨=+⎩的形式，00u v 、反映了单元的刚体平动，ω反映了单元的刚体转动，满足完备性和连续性的要求①。

采用插值法由单元结点位移列阵δe=()ii j j mm u v u v u v T计算α1、α2、α3、α4、α5、α6.，求出位移d=[u （x ，y), v （x ，y ）]。

6个未知量,6个代数方程,得d e=N δede=u v ⎛⎫ ⎪⎝⎭=00000i j m ijm N N N N N N ⎛⎫ ⎪⎝⎭()ii j j mm u v u v u v T式中N i =(a i +b i x+c i y)/2A ，a i =j j mmx y x y b i = -11imy y c i =11jmx x （i 、j 、m 轮换）A 为三角形面积，为避免A<0，i 、j 、m 按逆时针排列。

N 为形函数矩阵，形函数Ni 的性质有：①N i （x i ，y i ）=1 N i （x j ，y j ）=0 Ni （xm ，ym ）=0 ②N i （x ，y ）+N j （x ，y ）+N m （x ，y ）=1可推出三个形函数中，两个是独立的，反映了刚体平移。

令z=Ni ，在直接坐标系中画出Ni 、Nj 、Nm 的函数图形是以Ni （xi ，yi ）=1为高的四面体，所以结点位移影响单元的位移场，单元的位移场是线性分布的，相邻单元在公共边上的位移是连续的，单元相邻边的位移只取决于单元相邻公共边上的结点而与其他结点无关，无论以哪个单元计算相邻边的位移，结果一定相同。

形函数N i e决定了单元内的位移模式，反映了i 结点位移对单元内任意点位移的贡献率。

2、根据几何方程用单元结点位移表示单元应变：()00010002ij m ij m i i j j mm i ijjmm u x B b b v c c c u v u v u v y A cb c b c b u v y x ε⎛⎫∂ ⎪∂ ⎪⎛⎫ ⎪∂ ⎪== ⎪ ⎪∂ ⎪ ⎪⎝⎭⎪∂∂+ ⎪∂∂⎝⎭Te e B εδ=B 为几何矩阵B 可写为分块矩阵B =（B i B j B m ）T,B i =00ii ii b c c b ⎛⎫⎪⎪ ⎪⎝⎭，B 内所有元素与x ，y 无关，所以该单元内应变是常量，反映单元的常量应变，满足完备性和连续性的要求 ，这是一种常应变单元。

3、根据物理方程用单元结点位移表示单元应力：e e D σε= 2101011002E D μμμμ⎛⎫ ⎪⎪= ⎪- ⎪- ⎪⎝⎭D 为弹性矩阵e e e e D DB S σεδδ=== S 为应力矩阵S=DB 中，每一个元素都是常数，所以e σ的每一个分量与单元内x ，y 位置无关，这是一种常应力单元。

因为在三结点三角形单元中，位移函数中含有坐标的一次项，其误差为()2o x ∆，而应力、应变是常量，其误差为()o x ∆，比位移精度低。

4、根据虚功原理用单元结点位移表示单元结点力单元在结点处受力，单元会发生变形，因此单元在结点处所受到的力与单元结点位移肯定有关系。

单元间通过结点的相互作用成为整体，因此每一单元的受力——位移关系找出来，整体的受力——位移关系也就出来了。

记单元节点力为()e ijm F F F F =T，单元结点虚位移为()*e***i j m =Tδδδδ单元内应力为()e xy xy σσστ=T， 单元内虚应变()****x yxy εεεγ=T根据虚功原理，()()**TTe e e AF dxdy t δεσ=∙⎰⎰，可得e T e AF B DBdxdy t δ=⋅⋅⎰⎰因为B 、D 中元素都是常数，e T e e F B DBtA K δδ==，K=B TDB tA 为单元刚度矩阵。

K 为6行6列矩阵可写为()T i T j i jm T m B K tA B D B B B B ⎛⎫ ⎪== ⎪ ⎪⎝⎭iiij im jij j jm mimjmm k k k k k k k k k ⎛⎫ ⎪ ⎪ ⎪⎝⎭，xxxy ij ij ij i j yxyy ijij k k k B DB At At k k ⎛⎫== ⎪ ⎪⎝⎭，xyij k 表示j 结点处发生y 方向的单位位移时所引起的i 结点处x 方向的结点力。

不同类型不同形式的单元，只有弹性矩阵D 和几何矩阵B 不同，计算子块矩阵的公式相同，平面问题中，影响刚度矩阵K 的只有几何矩阵B 。

K 的性质有：K 中每个元素表示个单元结点沿坐标方向发生单位位移时所引起的结点力。

②K 为对称矩阵。

③单元做刚体位移时，单元内不产生应变应力，结点力为0，所以K 中每行每列元素之和为0，所以0K=，所以只根据e e F K δ=无法求得唯一解。

5、根据虚功等效原则计算等效结点力根据有限元的基本方法，单元内任意点的位移、应变、应力等最终都要用结点位移来表示，所以作用在物体上的外力也要用结点位移表示。

为了计算等效结点力，在任意的虚位移上，使原载荷与等效载荷虚功相等。

设外力为p f ，结点虚位移为*e δ，则任意点虚位移为**e e d N δ=，等效节点载荷为e L F ,有*eT *eT e p L d f t F δ= e T L p F N f t =（集中力）同理得e T L S F N f ds t =⋅⎰（面力），e L A F Nfdxdy t =⋅⎰⎰（体力）。

三、整体分析将结构的所有单元通过结点连接起来，形成一个整体的离散结构以代替实际的连续体，以形成以结点位移为未知量的整体结构的有限元代数方程组，最后求得结点位移。

对结点受力分析：结点受到与之相关的单元给它的反作用力和外载荷的等效结点力，这两组力坐标轴方向相反，所以应该相等，即i Li eeF F =∑∑，设有n 个结点，每个结点建立两个方向的方程，不考虑外界约束时，共2n 个方程，2n 个未知量（,,...ix iy jx δδδ），为了建立这个代数方程组，建立整个弹性体的结点力和结点位移的关系式L K F δ=，K （2n ×2n ）为整理刚度矩阵，δ为整体结点位移列阵，F L 整体结点载荷列阵。

为了求整体刚度矩阵，要找到它与已求得的单元刚度矩阵的关系，在整体中对结点编码，设整体刚度矩阵中某元素为Kij ，意为j 个结点在x 或y 方向发生位移引起i 个结点x 或y 方向的结点力，找到同时用到i 与j 结点的单元，并用与之对应的单元刚度矩阵中的元素ksm 相加得到Kij ，整体刚度矩阵也是奇异矩阵，必须考虑边界约束条件，消除K 的奇异性，才能求解结点位移。

再由单位结点等效载荷得到整体结点载荷列阵F L 。

这样K 、F L 已知，求解代数方程，解出整体结点位移列阵δ，得到相应的单元结点位移δe。

δe得到了，相应的d e、σe、εe等就得到了。