《有限元法的基本原理.docx》由会员分享，可在线阅读，更多相关《有限元法的基本原理.docx（20页珍藏版）》请在三一办公上搜索。
1、第二章有限单元法的基本原理作为一种比较成熟的数值计算方法，有限元的数学基础是变分原理。经过半个过 世纪的发展，它的数学基础已经比较完善。从数学角度分析，有限元法是以变分原理和剖分插值为基础的数值计算方法。它 广泛的应用于解算各种类型的偏微分方程，特别对椭圆型方程，因为椭圆型方程的边 值问题等价于适当的变分问题，即能量积分的级值问题。通过变分，导出相应的泛涵， 再把作用域从几何上剖分为足够小的单元，这样就能够用简单的图形去拟合复杂的边 界，用简单的初等函数去模拟单元的性质。在解算中先对每个单元进行分析，后在通 过连接单元的节点对作用域的整体进行分析，就是对泛涵求极值，从而把一个复杂的 偏微分方程
2、求解问题，变成解线形代数方程组的问题。尽管这样会出现大量的未知数， 由于采用了矩阵分析的方法，总体上很有规律，适合编制程序用计算机完成。通常的数学考虑包括这些：1）从古典变分方法原理去定义微分方程边值问题的 广义解以及在古典变分方法的框架对有限元进行理论分析。2）保证偏微分方程边值 问题的提法正确，即要求解存在、唯一和稳定，即保证数值解法是可靠的。3）有限 元中重要的一点是采用了分块多项式插值函数，因此，有限元的误差估计转化为插值 逼近的误差估计问题。4）有限元的收敛性和误差估计。由于本文是应用有限元的理论解决大地测量中的问题，因此，这里将不讨论上叙 问题，而是从固体力学的基本方程出发，通过虚
3、功原理建立起离散化的有限元方程。 另外，还以八节点六面体单元为例，简要叙述了实际中最常用的等参单元的概念及其 数值变化的一些公式。2.1弹性力学基本方程有限元法中经常要用到弹性力学的基本方程，这里写出这些方程的矩阵表达式。2-1-1、平衡方程对任意一点的受力情况分析，沿坐标轴方向x, y ,z分解得到平衡方程记为:d dx00d8y0d dz0d8y0d dxd dz000d dz0d8yd dzAo + F = 0o*xooyzTxyT+XZFXFyFZ其中A是微分算子，F是体积力向量。2-1-2、几何方程位移应变关系在微小位移和微小变形的情况下，略去位移倒数的高次幂，得到应变向量和位移向量
4、的几何关系:记为sXsysZYXyYyzLyd00dx0d08y00ddzdd08ydx0dddz8yd0ddz8yuVws = Bu其中B = AT2-1-3、物理方程应力-应变关系弹性力学中，应力应变关系也叫弹性关系，对于各向同性的弹性材料，有o = Ds弹性矩阵D为对称矩阵，它完全取决于弹性体材料的弹性模量E和泊松比vbXbyb z = E (IFT (1 + v)(1 - 2v)XyTyzT1v1 一 v000v1 一 v1v1 一 v0v1 一 vv1 一 v10000001 一 2v2(1-v)00000000001 一 2v2(1-v)000001 - 2vr _ n es X8
5、 ye* zY xy丫 yzY2(1 - v)2-1-4、力的边界条件在边界上S1已知弹性体单位面积上作用的面积力为P ,单位面积的内力为T :T-P 一xxTPyyTPzz进一步得:bTTcos( n, x)-P XxyxzxTbTcos(n, y)PyxyyzyTTbcos( n, z)PL zxzyz z其中cos(n,x), cos(n,y), cos(n,z)分别是边界面的法线与坐标轴的夹角。则有T = P2-1-5、几何边界条件在边界S2上弹性体的位移已知为Ur u r 一v=-Vww则有u = u2-1-6、弹性体的应变能应变能公式为：U(8 ) = 18 TD&综上，弹性力学方
6、程为：徵=/8 = Bub = D8T = Pu = u 2.2虚功原理和有限元方程的建立2-2-1虚功方程从平衡方程BTb + F = 0出发。这是一微分方程，同时它满足一定的边界条件。 为了求导出广义解，将它化为积分等式，取任意向量函数中3,y,z) = k w x，它在S1满足位移边界条件，用T乘方程的两端，在空间Q上求积分得到，0 = jjj T (Btb + F)dxdydz通过一系列数学变换，最终可得到:JJJ( B )t D( BU )dxdydz =J!J t Fdxdydz + J! TPdSQQS 2上述方程就是我们通常所说的虚功方程。中(X, y,z)称为虚位移，B中称为
7、虚应变8 *。因此上述方程的左端可表示成为JU (8*)tDbdxdydz，而上述等式就是虚功Q方程，它表示内力做功与外力相等。满足虚功方程同时又满足边界条件的向量函数U3, y,z)，即是微分方程的广义解。2-2-2离散化建立有限元方程1)单元剖分通常的单元有四面体、菱柱体和六面体单元。由于六面体等参元便于剖分拟合实际地形并且有较高的计算精度，在此我们以八 节点六面体单元为例进行说明。不妨假定区域Q是一多面体，做六面体剖分，单元e n(n=1,.,NE)，节点(六面体顶点)P 3, y ,Z )(I=1,.,NP)。 i i i i2)线性插值设在每个节点p.(七,y, z)上，位移u的值为
8、 vwi则对于任意单元e，由U在八个顶点上的值，恰好在单元内确定一线形插值。即 若令U - a x + a y + a z + a1234其中ai是 3*1列阵。它在顶点上满足U = a x + a y + a z + a (i=1,.,8) i 1 i 2 i 3 i 4 解出a.，再带回上式，并且经过适当整理可得到：U = N5 eN001N =0N000N1-iN0080N0800N8这里只是给出简要过程，在下一章我们将详细说明N的具体表达式及编程实现的方法。由几何方程和物理方程，应变e和应力。可表示为:e = B6e其中B = B .1而 B.(i = 1,.,8)是 6*3 矩阵I，
9、b = DB6eB8 8N0i-08x0dN0i-dy00dNi-B =dNdNdzi-0dydx0dNdNi-idzdydN0dNi-i-|_ dzdy _3）单元刚度矩阵与单元荷载向量剖分之后，虚功方程可表示为正 jjje *Tbdxdydz = 正 jjjTFdxdydz + 为n=1 en=1 en=1nnjjje *Tbdxdydz =jjjB6 * tDB6 dxdydz = 6 *t k 6eee e eAAAee故单元刚度矩阵为：k = i BrDBdV它是一个24*24的对称、非负矩阵，当它为零时，其力学意义是应变等于0,即 单元不发生形变，做刚性位移，它反映了单元的刚性，它
10、表明为了维持单元e的形变，我们需要在单元的节点上施加外力，使它达到平衡，这些外力是通过节点对单元作用， 称为等效节点力。为计算单元荷载向量需要计算积分fffo TFdxdydz 和 JJJ TPdSS2 nae容易知道jjjTFdxdydz = fff(N6* )tFdxdydz = 5* tfee其中F = Fe.FeJJJF N (x, y, z)dxdydz x iF = Jjf F N (x, y,z)dxdydzfff F N (x, y, z)dxdydzz iefffo TPdS =5*t ffNTPds =5*tFS 2 naeS2 nae其中Li L - 1F = F . F
11、eTFfff P N (x, y, z)dS s 2 naejJJ P N (x, y,z)dS snaeJJJ P N (x, y, z)dSLS 2 nae(I=1,.,8)这些复杂的积分运算，需要用到数值积分的方法，通常用高斯积分方法，具体实 现在下一章进行。4）总刚矩阵与总载荷向量的安装把单元刚度局阵，单元载荷向量的表达式代入等式AAA。每个单元刚度系数（24*24 矩阵），每个单元载荷向量（24*1列阵）按其脚标编号对号如座形成总刚矩阵和总荷 载向量。此时的总刚度矩阵是3NP阶方阵，总载荷向量是3NP维向量。叠加的结果为：芸5*t k 5=芸5*t F + Fnnnnnn = 1n
12、= 15*t（芸k ）5 = 5*t芸F + Fenenenn=1n=1F =芸F + F enn =1从而由上述方程得线形方程组5 *t （ K 5 - F ）= 0这里K =正k e nn = 1由于5 *是个任意NP维向量K5 * = F5）约束处理、解方程组解方程得到位移（1=1,NP），然后计算应 i带入边界条件，进行相关处理后，力。由于K的对称、带状等特征，在解方程时可以用一些特殊的方法，这里不加讨 论。这里从虚功原理出发推导了有限元的计算公式，另外，用加权余量法，变分原理？ 等同样可以推出这些公式。因为有限元的数学理论已经成熟，而且本文的目的是应用 这一数学工具得到华北地区的具体
13、情况，因此在这里没有展开讨论。具体可参见王瑁 成、姜礼尚的著作。 2.2等参单元的概念及空间等参元法选择恰当的插值函数来表达单元的位移分布，关系到有限元的计算精度和计算效 率，因此往往是研究的重点之一。等参元是目前大型有限元程序中应用最广泛的单元， 它不仅适用于各种曲线边界，而且能够构造出高精度的位移函数。所以在这里具体以 三维等参元为例给出其概念，并写出其变换规则。2-3-1等参单元的概念X = NXiN100 .N 800 _N =0N0 .0N01800N .00N118*8 J8 Z坐标变换和插值函数是以节点值为参数，并且参数的数目相同，采用的基函数也相同， 称为等参元。对于等参元，有
14、两套坐标系，整体坐标和自然坐标(局部坐标)。取如下插值函数将 下图的立方体映射成任意直边六面体的变换函数：其中：尤x= y_ z _X = y1 z1用试凑法求得N .:IN (r,s,t) = 8(1 + r r)(1 + s.s)(1 +11)IS2根据定义，坐标变换和插值函数是采用相同的基函数，得到位移函数为:U = N8 eU为位移函数，5 e为单元节点位移（沿整体坐标轴方向）2-3-2三维等参单元的一些导数变换在求应变矩阵时，要求位移对坐标的导数，下面推导位移对整体坐标和局部坐标 导数之间的关系。在空间问题中，位移u,v,w是整体坐标x,y,z的函数，而x,y,z又是局部坐标的 函数
15、，所以du dx du dy du dz=+drdx drdy drdz drdudu dxdu dydu dz=+dsdx d sdy dtdz dtdudu dxdu dydu dz一=+dtdx dtdy dtdz dt写成矩阵式:dudxdydzdududrdrdrdrdxdxdu=dxdydzduJ =Jdudsdsdsdsdydydudxdydzdududt _dtdtdt _dz _dz _其中J是雅可比矩阵:- &一araz一 csaz-初 空aray- csay-初 办一arax-csax-初 - _=-3 3 312 3J J JJJ1112J =JJ2122JJL 313
16、2令:J = J -1 =J11J必1J/3尸则有:31dududxddudu=Jdydsdudu云-0 _同理dvdv8w8w瓦济济8v8w=J， ，=J8s8y8s8v8v8w8w石_a _瓦将上面三个式子组成一个矩阵式，得出导数之间的变换式:8u8u云8T8u8u8y8s8u8u8z瓦8v8v云J 0 08r8vetc8v0 J 08一 一 瓦8v0 0 J8v云瓦8w8x8w8w8y宙8w8wbz J_8t -J，J，J 2的具体表达式为:1112138x8N8NJ=1 x +.+8 X118r8r 18r 88y8N8NJ1 y +.+ 8 y128r8r 18r88z8N8NJ1
17、z +.+耕z138r8r 18r 8同理可得J，J，J，J，J， J，将所有式子写成矩阵形式有:212223313233dNidrdNidsdNidtdN 2 drdN2ds dN 3 dtdN8 drx1y1z1dNxyz8222ds dN8xyzdt8882-3-3体积、面积微元的变换由矢量乘积计算体积微元：dV = dr (ds x dt) = Jdrdsdt面积微元，以r = 1为例：dS = ds x dt|在计算体积力(如重力)，表面力时，会用到上面的公式。具体推导过程这里从略。2-3-4八节点等参元的应变矩阵：- - aw一axcs 一aycs 一azav一axav一ayav-
18、azaw一办空ayaw一az- _OOO O O O O O O O1O O O O 1 O O O O8OO=-88y8z U 兴-O O O O O 1 O 1 O =V一*O 8 81 O O O O O O O OO 1 O 1 O O O O OO O 1 O O O 1 O OL88O导数之间变换有雅可比矩阵的逆矩阵得到:- - aw一araw一 角aw一azav一arav一 角av一azav一 禽如一角如一az- -J J J.11人12人13r nL-1J J Js100000000212223,J J Js000010000313233人人人yJ J Js0000000011
19、11213zJ J JV0 1 0 1 0 0 0 0 0人21人22人23xyJ J JV0 0 0 0 0 1 0 1 0313233yzJ J JV0 0 1 0 0 0 1 0 0.11人 12人 13zx 1J J J人21人22人23J J J313233将位移表达式U = N8 e对局部微分，得到：1 V 1 w 1 I 8 V 8 W 8- -R & o o o o o o N 禽W禽W金8 8 8只， R 只o o OW禽洲禽洲金o o O8/ 8 R欲苗欲右欲*。 o o o o o O欲*欲*欲*O 。欲*式欲* 。 。 O欲-a;1欲*也金o O=o o o O如一禽如
20、禽如金加禽所禽加金所禽揶一禽揶金综合上两式，可以得到 = B8 er八JJJ000000111213八八八000JJJ000212223八八八000000JJJB =八八八八八八313233JJJJJJ000212223111213八八000JJJJJJ八八八313233212223JJJ000JJJL 313233111213中的B矩阵为：dN i dr00dN2-udr00dN i ds00dN2-uds00dN1dt00dN.8- dt000dN 1 dr0.0dN-dr00dN 1 ds0.0dN 8- ds00dN 1 dt0.0dN2-dt000dN 1 dr.00dN81dr00dN 1 ds.00dN8ds00dN 1 dt.00dN8dt
