有限元法的基本原理.doc

第二章 有限单元法的基本原理作为一种比较成熟的数值计算方法，有限元的数学基础是变分原理。经过半个过世纪的发展，它的数学基础已经比较完善。从数学角度分析，有限元法是以变分原理和剖分插值为基础的数值计算方法。它广泛的应用于解算各种类型的偏微分方程，特别对椭圆型方程，因为椭圆型方程的边值问题等价于适当的变分问题，即能量积分的级值问题。通过变分，导出相应的泛涵，再把作用域从几何上剖分为足够小的单元，这样就能够用简单的图形去拟合复杂的边界，用简单的初等函数去模拟单元的性质。在解算中先对每个单元进行分析，后在通过连接单元的节点对作用域的整体进行分析，就是对泛涵求极值，从而把一个复杂的偏微分方程求解问题，变成解线形代数方程组的问题。尽管这样会出现大量的未知数，由于采用了矩阵分析的方法，总体上很有规律，适合编制程序用计算机完成。通常的数学考虑包括这些：1）从古典变分方法原理去定义微分方程边值问题的广义解以及在古典变分方法的框架对有限元进行理论分析。2）保证偏微分方程边值问题的提法正确，即要求解存在、唯一和稳定，即保证数值解法是可靠的。3）有限元中重要的一点是采用了分块多项式插值函数，因此，有限元的误差估计转化为插值逼近的误差估计问题。4）有限元的收敛性和误差估计。由于本文是应用有限元的理论解决大地测量中的问题，因此，这里将不讨论上叙问题，而是从固体力学的基本方程出发，通过虚功原理建立起离散化的有限元方程。另外，还以八节点六面体单元为例，简要叙述了实际中最常用的等参单元的概念及其数值变化的一些公式。2.1 弹性力学基本方程有限元法中经常要用到弹性力学的基本方程，这里写出这些方程的矩阵表达式。2-1-1、 平衡方程对任意一点的受力情况分析，沿坐标轴方向x, y ,z分解得到平衡方程 记为： 其中A 是微分算子，F 是体积力向量。2-1-2、几何方程-位移 应变关系在微小位移和微小变形的情况下，略去位移倒数的高次幂，得到应变向量和位移向量的几何关系：记为 其中 2-1-3、物理方程-应力-应变关系弹性力学中，应力应变关系也叫弹性关系，对于各向同性的弹性材料，有弹性矩阵 D 为对称矩阵，它完全取决于弹性体材料的弹性模量E和泊松比 v 。2-1-4、力的边界条件在边界上S1已知弹性体单位面积上作用的面积力为， 单位面积的内力为：进一步得： 其中cos(n,x), cos(n,y), cos(n,z) 分别是边界面的法线与坐标轴的夹角。则有 2-1-5、几何边界条件在边界S2上弹性体的位移已知为则有 2-1-6、弹性体的应变能应变能公式为： 综上，弹性力学方程为： 2.2 虚功原理和有限元方程的建立2-2-1虚功方程从平衡方程出发。这是一微分方程，同时它满足一定的边界条件。为了求导出广义解，将它化为积分等式，取任意向量函数，它在S1满足位移边界条件，用乘方程的两端，在空间上求积分得到，通过一系列数学变换，最终可得到：上述方程就是我们通常所说的虚功方程。称为虚位移，称为虚应变。因此上述方程的左端可表示成为，而上述等式就是虚功方程，它表示内力做功与外力相等。满足虚功方程同时又满足边界条件的向量函数，即是微分方程的广义解。2-2-2离散化建立有限元方程1) 单元剖分通常的单元有四面体、菱柱体和六面体单元。由于六面体等参元便于剖分拟合实际地形并且有较高的计算精度，在此我们以八节点六面体单元为例进行说明。不妨假定区域是一多面体，做六面体剖分，单元 (n=1,NE)，节点（六面体顶点）(I=1,NP)。2）线性插值设在每个节点上，位移U的值为则对于任意单元，由U在八个顶点上的值，恰好在单元内确定一线形插值。即若令其中是3*1 列阵。它在顶点上满足 (i=1,8)解出，再带回上式，并且经过适当整理可得到： 这里只是给出简要过程，在下一章我们将详细说明N的具体表达式及编程实现的方法。 由几何方程和物理方程，应变和应力可表示为： ， 其中 而是6*3矩阵3）单元刚度矩阵与单元荷载向量剖分之后，虚功方程可表示为 AAA 故单元刚度矩阵为：它是一个24*24的对称、非负矩阵，当它为零时，其力学意义是应变等于0，即单元不发生形变，做刚性位移，它反映了单元的刚性，它表明为了维持单元e的形变，我们需要在单元的节点上施加外力，使它达到平衡，这些外力是通过节点对单元作用，称为等效节点力。为计算单元荷载向量需要计算积分 和容易知道其中 （I=1,8）其中 （I=1,8）这些复杂的积分运算，需要用到数值积分的方法，通常用高斯积分方法，具体实现在下一章进行。4）总刚矩阵与总载荷向量的安装把单元刚度局阵，单元载荷向量的表达式代入等式AAA。每个单元刚度系数(24*24矩阵)，每个单元载荷向量（24*1列阵）按其脚标编号对号如座形成总刚矩阵和总荷载向量。此时的总刚度矩阵是3NP阶方阵，总载荷向量是3NP维向量。叠加的结果为：即 故这里 ，由于是一个任意NP维向量，从而由上述方程得线形方程组5）约束处理、解方程组带入边界条件，进行相关处理后，解方程得到位移(I=1, , NP)，然后计算应力。 由于K的对称、带状等特征，在解方程时可以用一些特殊的方法，这里不加讨论。这里从虚功原理出发推导了有限元的计算公式，另外，用加权余量法，变分原理？等同样可以推出这些公式。因为有限元的数学理论已经成熟，而且本文的目的是应用这一数学工具得到华北地区的具体情况，因此在这里没有展开讨论。具体可参见王瑁成、姜礼尚的著作。2.2 等参单元的概念及空间等参元法选择恰当的插值函数来表达单元的位移分布，关系到有限元的计算精度和计算效率，因此往往是研究的重点之一。等参元是目前大型有限元程序中应用最广泛的单元，它不仅适用于各种曲线边界，而且能够构造出高精度的位移函数。所以在这里具体以三维等参元为例给出其概念，并写出其变换规则。2-3-1等参单元的概念坐标变换和插值函数是以节点值为参数，并且参数的数目相同，采用的基函数也相同，称为等参元。对于等参元，有两套坐标系，整体坐标和自然坐标（局部坐标）。取如下插值函数将下图的立方体映射成任意直边六面体的变换函数：其中： 用试凑法求得: 根据定义，坐标变换和插值函数是采用相同的基函数，得到位移函数为：U为位移函数，为单元节点位移（沿整体坐标轴方向）2-3-2三维等参单元的一些导数变换在求应变矩阵时，要求位移对坐标的导数，下面推导位移对整体坐标和局部坐标导数之间的关系。在空间问题中，位移u,v,w是整体坐标x,y,z的函数，而x,y,z又是局部坐标的函数，所以写成矩阵式：其中J是雅可比矩阵：令：则有：同理将上面三个式子组成一个矩阵式，得出导数之间的变换式： ，的具体表达式为：同理可得， ，将所有式子写成矩阵形