有限元法
一种将连续体离散化,以求解各种力学问题的数值方法,又称有限单元法或有限元素法。
20世纪50年代M.J.特纳和L.C.托普等人把求解杆件结构的方法推广到求解连续体力学问题并在数学上采用矩阵表示法,对有限元法的早期发展起了重要作用。1960年R.W.克拉夫首先使用“有限元”这一名称。60年代,由于电子计算机的广泛使用,有限元法得到很快的发展。经过二十多年的发展,有限元法已经成为处理力学、物理、工程等计算问题的有效方法之一。
从应用数学的角度看,有限元法基本概念的提出,可以追溯到R.库朗1943年的工作,他采用三角形单元组成分区近似函数,并用最小势能原理,讨论了柱体的扭转问题。由于当时没有求解大型联立方程的计算工具,这种方法在长期内没有实际应用。中国刘徽早在公元3世纪时就用割圆术求圆周率。这种把圆周分割成有限个单元,用有限来逼近无限的思想,可说是现代有限元法的早期萌芽。
基本原理 有限元法的离散化,可从物理和数学两个角度来理解。从物理角度看,一个连续体可以近似地用有限个在节点处互相连接的单元所组成的组合体来代表,因而可把连续体的分析问题变成单个单元的分析和所有单元的组合问题。从数学角度看,一个连续域可以分割为有限个子域,每个子域的场函数是只包含有限个参数的简单场函数,用这些子域的场函数的集合,就能近似地代表整个连续域的场函数。于是,求解连续场函数的微分方程就转化为求有限个待定参数的代数方程组。总之,离散化就是化无限为有限,从而达到化难为易的目的。
推导方法 有限元法所用的基本推导方法有三种:
①直接法 它根据物理概念,直接建立待求问题的方程式,就象在解杆系结构力学问题中用位移法或力法建立方程一样。这种方法虽有概念明确、易于理解的优点,但一般只能处理比较简单的问题。
②变分方法 根据能量泛函的驻值条件,建立求解问题的方程式。这种方法可处理各种复杂问题。在固体力学领域里,能量变分原理是有限元法获得发展的一个重要基础。
③加权残数法 如已知问题的微分方程,但无可用的变分泛函表达形式,可以用加权残数法处理。伽辽金法就是加权残数法的一种。
固体力学问题有多种解法。基本未知量如选取为节点位移,称为位移法或刚度法;如选取为节点力,称为力法或柔度法;如选取一部分为节点位移,另一部分为节点力,则称为混合法。这三种方法中,位移法有易于实现自动计算的优点。因此,在有限元法中,以位移法的应用最广。混合法应用较晚,但在板壳问题和接触问题中已经显示出它的优点。
采用不同的变分原理,可以得到不同类型的有限元模型,其中有基于最小势能原理的协调型模型,有基于最小余能原理的平衡型模型,有基于广义变分原理的杂交模型和混合模型等。
计算过程 现以弹性力学平面问题为例,说明有限元位移法的计算过程。它包含三个主要步骤:离散化、单元特性分析和整体组合分析。
①离散化 图1是一个重力坝及其基础的横截面图。可按弹性力学平面问题进行计算。用有限元法须先把原来的连续体分割成三角形网格,使它变成由有限个在节点处用铰相连的三角形单元所组成的离散体系(图2)。外界的载荷也等效地转移到节点上,在位移为零的边界节点处设置支座。
图2
②单元特性分析 图 3是上述重力坝离散体系中的一个典型的三角形单元。在弹性力学平面问题中,单元的每个节点有沿x、y轴的两个位移分量u、υ,一个单元有三个节点,共有六个节点位移分量。单元的节点位移矢量矩阵{δ}可记为:
(1)
这就是问题的基本未知量。与基本未知量对应的物理量是图4所示的节点力分量。节点力矢量矩阵可记为:
(2)
节点力矢量与节点位移矢量之间有如下关系:
(3)
式中[k]称为单元的刚度矩阵,其物理意义是:矩阵中的每个元素表示单位节点位移引起的节点力。单元特性分析的主要任务就是求单元的刚度矩阵[k]。
图3
图4
单元特性分析的步骤是:首先确定单元内的位移插值函数,即把单元内任一点的位移用节点位移{δ}表示;然后写出单元的总势能;最后应用势能驻值条件,即可导出刚度方程(3),从而得出单元的刚度矩阵[k]。
③整体组合分析 得到单元刚度矩阵[k]后,再将各个单元组合成整个结构,求出整个结构的刚度矩阵[K],这就是整体组合分析的主要任务。
整体分析的步骤是:首先对所有单元采用同一个公共坐标,经过坐标转换,建立在公共坐标系中的单元刚度矩阵;然后,把各个单元刚度矩阵进行组装,得出结构的整体刚度矩阵[K],从而得出整个结构的结点力{F}和节点位移{δ}之间的刚度方程:
(4)
由式(4)可解出结点位移,进而可求出各单位的应力。
单元类型举例 随着有限元法的发展,各种类型的新单元不断出现,它们在单元的几何形状,基本未知量的性质和个数,以及插值函数的形式等方面都各有特点。现将位移法的几种常用单元形式列于表中。
表中序号1、2是一维单元,3、4、5、6是二维单元,7是三维单元。此外,二维或三维的等参单元也是常用的。为了适应不同问题的需要,近年来还发展出奇异元、无限元、边界元、伸缩元等新型单元。如果要求解的问题规模较大,为了提高计算效率,常把整个结构划分为若干个子结构,先对子结构进行特性分析,然后再做整体组合分析。这样,子结构实质上就是大型的单元。如果问题规模很大,而计算机容量有限,还可以把结构化为多层的子结构。
在其他领域的应用 经过二十多年的发展,用有限元法解线性范围内的固体力学和工程结构方面的问题已经相当成熟。除了在线性静力分析上的应用外,在其他领域里的应用也在不断发展,这些领域是:
①动力分析 在动力分析中,常用有限元法来求结构的自然频率和振型。例如,对于无阻尼的自由振动问题,结构的运动方程为:
[K]{δ}=-[Μ]{},(5)
式中[K]和{δ}的含义与静力问题相同。只是位移{δ}是时间t的函数,[Μ]是结构的质量矩阵。在自由振动中,位移可表为:
{δ}={δ}sinωt,(6)
而加速度为:
{}=-ω{δ}sinωt,(7)
ω为结构自由振动时的圆频率。将式(6)和(7)代入(5)得:
([K]-ω[Μ]) {δ}={0},(8)
这就是一个典型的求特征值问题。
②非线性问题 根据引起非线性反应的原因,结构分析中的非线性问题可分两类,即物理非线性和几何非线性问题(见计算结构力学)。这两类非线性问题的解法基本一致,一般均采用迭代法或增量法。不管用哪一种方法,每一步计算都相当于一次线性分析。
③稳定问题 研究结构的稳定性,必须考虑几何非线性的影响。若失稳前结构的反应是线性的,可以只考虑几何刚度影响。这时结构的刚度矩阵[K]可由下式给出:
[K]=[K]+[K],(9)
式中[K]为弹性刚度矩阵,[K]为几何刚度矩阵,在载荷增量{F}的作用下,结构的平衡方程为:
{F}=[K]{δ}。(10)
对于稳定问题,实际上只要考虑由一个给定的平衡状态开始的第一步增量,而且这一步的载荷增量为零,于是由式(10)可得:
([K]+[K]){δ}={0},(11)
所以确定失稳时的临界载荷和屈曲形态问题也是一个典型的特征值问题。
④流体力学问题 有限元法在固体力学中的成就,启发人们把它用于流体力学,例如飞行器的设计,要设计好飞行器,须知道空气对飞行器的作用力。由于飞行器的形状复杂,用原有流体力学的数值方法来处理这类问题就有困难。有限元法则可以灵活地采用各种形状和尺寸的单元,因而能够描述复杂的形状,为解决这类问题开拓新的途径。但流体力学毕竟与固体力学有不少差别。差别之一是流体的位移很大。因此,流体力学方程(如欧拉方程组或纳维-斯托克斯方程)就有很多非线性项,例如,流场中的加速度就由几个非线性项表示。此外,流体力学中的边界条件也很不同。因此,把有限元法用于流体力学就具有自己的特点。用数学术语来说就是,如果流体力学问题能化成线性椭圆型方程,用有限元法就比较合适。近几年来,已有人开始用有限元法研究跨声速飞行器问题,也有人用它来解用双曲型方程描述的波动问题。这些发展趋势引起研究者较大的兴趣。
发展方向 目前有待研究的主要课题可归纳为如下几个方面:①对有限元法的数学基础和变分原理的进一步研究;②对新型单元,如奇异元、无限元等的研究;③对非线性问题的研究;④对有限元法与其他离散化方法或解析方法综合应用的研究;⑤对在力学领域以外应用的研究;⑥对效能和自动化程度更高的有限元软件系统的研究等。
参考书目 O,C.Zienkiewiez, The Finite Element Method,McGraw-Hill,London,1977. R.H.Gallagher,Finite Element Analysis,Fundamentals,Prentice-Hall,Englewood Cliffs,New Jersey,1975. 龙驭球编:《有限元法概论》,人民教育出版社,北京,1978。 《1980年全国计算力学会议文集》,北京大学出版社,北京,1981。