搜索
有限元法有限元法是现代工程数值分析中应用广泛的一种方法，是解决工程实际问题的一种有力的数值计算工具。昀初，这种方法是被用来研究复杂飞机结构中的应力，是将弹性理论、计算数学和计算机软件有机地结合于一起的一种数值分析技术。后来由于这一方法的灵活、快速和有效性，使其迅速发展成为求解各领域中的数理方程的一种通用的近似计算方法。目前，它在许多科学领域和实际工程问题中都得到广泛的应用，因此，在工科院校和工业界受到普遍的重视。在求解工程技术领域的实际问题时，建立基本方程和边界条件相对来说还是比较容易的，但是由于其几何形状、材料特性和外部荷载的复杂性和不规则性，使得求其解析解答确实十分的困难。因此，寻求近似解法就成了必由之路。经过多年的探索，近似算法发展了许多种。但常用的数值分析方法往往或多或少地存在着局限性。例如：我们在前面第二章介绍的有限差分法计算模型可给出其基本方程的逐点近似值（差分网格中的点），但对于不规则的几何形状和不规则的特殊边界条件差分法就难于应用了；而加权残值在应用上很大程度依赖于权函数和试函数的正确选择，需要相当多的经验积累才能做到。有限元法是将求解区域看作由许多小的在节点处互相连接的子域（单元）所构成，其分析模型是给出基本方程的（子域）分片近似解。由于单元（子域）可以被分割成各种形状和大小不同的尺寸，所以它能适应复杂的几何形状、复杂的材料持性和复杂的边界条件。“有限元法”这个名称第一次出现在1960年，是由美国的克拉夫（Clough）在一篇关于平面弹性问题的论文中首先使用的。但是有限元法分析的概念却可以追溯到更远的20世纪40年代。1943年，Courant第一次在他的论文中，取定义在三角形区域上的分片连续函数，利用昀小势能原理研究了圣维南扭转问题。然而，这一方法当时并没有引起人们的注意，因而发展缓慢，几乎过了十年才有人再次采用这些离散化的概念。1956年Turner，Clough，Martin和Topp等人，第一次给出了采用三角形单元所求得的平面应力问题的真正解答，利用弹性理论的方程求出了三角单元的特性，并第一次介绍了我们所熟知的确定单元特性的直接刚度法。在1960年Clough进—步处理了平面弹性问题之后，工程师们开始认识了有限元法的功效，此后有限元法在工程界获得了广泛的应用。到20世纪70年代以后，随着计算机和软件技术的发展，有限元法也随之迅速地发展起来，进入了有限元法的鼎盛时期，并开始对该方法进行全面而深入的研究。涉及的内容有：有限元法在数学和力学领域所依据的理论；单元的划分原则、形状函数的选取及协调性；有限元法所涉及的各种数值计算方法及其误差、收敛性和稳定性；计算机程序设计技术等；以及其向其他各领域的推广。关于有限元的渊源，不同领域的学者持有不同看法。数学家们认为有限元法是偏微分方程定解问题数值分析方法中的一种新进展；而力学家们认为其是基于变分原理的近似解法的革新和发展；工程结构专家更倾向认为有限元法是结构力学的矩阵分析法在连续介质力学中的推广与应用。事实上，各个领域的学者们在各自不同的时期，从不同的角度出发，独立或者一起发展了有限元法，并在各自领域中作出了积极的贡献。再加上有限元法有成熟的大型软件系统支持，使其已成为非常受欢迎的、应用极广的一种数值计算方法。特别一提的是我国在有限元法领域的发展。我国是在20世纪50年代末期成功研制了第一台大型电子计算机（130型），当时著名数学家冯康教授以大型水坝结构的应力分析为开端，进行了大量的研究探讨工作，独立于西方创立了一套现代化和系统化的求解微分方程的近似方法，并取名为“基于变分原理的差分格式”。这一方法的内容在实质上就是当时国际上差不多同一时期提出的所谓的“有限元法”，并先于西方建立了严密的理论体系，是国际公认的当代计算数学的一项重大成就，不同的只是我国是从数学方面提出有限元法的。到目前为止，有限元法已被应用于固体力学、流体力学、热传导、电磁学、声学、生物力学等各个领域，能求解由杆、梁、板、壳、块体等各类单元构成的弹性（线性和非线性）、弹塑性或塑性问题（包括静力和动力问题）；能求解各类场分布问题（流体场、温度场、电磁场等的稳态和瞬态问题）；还能求解水流管路、电路、润滑、噪声以及固体、流体、温度等相互作用的问题。此外，有限元法有比较固定的一套分析顺序，对于不同的工程结构，往往可以使用同一个计算程序来解决，便于求解过程规范化，有高度的通用性。相关的有限元程序发展也很快，目前国外有名的主要有限元软件有：ASKA（结构分析自动系统），NASTRAN（NASA结构分析程序），SAFE（有限元结构分析程序），SAP系列（结构分析程序），ADINA，MARC，ANSYS等。有些程序还具备了前后处理功能，不仅解题的速度提高，还极大地方便了使用者，这对有限元法的普及与应用必然起到很大的促进作用。§1有限元的直观方法和基本概念本节首先从直观方法来阐述有限元法的基本概念和求解过程，这种直观法也是有限元早期发展的开端，这将有利于初学者对这一方法的思路和物理意义有一个直观的了解。杆、梁是工程技术人员所熟悉的标准结构构件，早在40年代，Hrenikoff,Mchenry,Newmark等就想用杆或梁的组合体来模拟弹性体的性质，直观地使用节点连结离散构件来逼近连续结构，并利用结构力学中的位移法或力法建立有限元计算格式，这种方法叫作直观有限元法，简称为直接法。直接法的优点是易于理解，但只能用于简单的问题。§1.1有限元法的直接法——杆的分析（一）单元的刚度矩阵在平面桁架分析中，我们经常会碰到杆元件，如图1.1所示。图1.1杆单元以及坐标系示意图记杆的两个端点分别为i和j，其节点位移和节点力分别可表示为⎥⎥⎦⎤⎢⎢⎣⎡=⎥⎦⎤⎢⎣⎡=iyixiiiiffvu][,][fu；⎥⎥⎦⎤⎢⎢⎣⎡=⎥⎥⎦⎤⎢⎢⎣⎡=jyjxjjjjffvu][,][fu（1.1）其中iivu,（或jjvu,）分别表示节点i（或j）沿坐标轴x方向和y方向的节点位移分量；yiixff,（或yjjxff,）为相应的x方向和y方向的节点力分量。并且规定：节点位移和节点力的符号与坐标轴yx,取向一致为正，反之为负。首先分析杆的应变-位移关系。杆的长度l可表示为222)()(ijijyyxxl−+−=（1.2）对其两边进行微分运算，得到)dd(sin)dd(cosdijijyyxxl−+−=θθ（1.3）其中，lyylxxijij)(sin,)(cos−=−=θθ。由于杆件受力变形，节点i的坐标将从),(iiyx改变为),(iiiivyux++，即iiiivyux==d,d（1.4）同理jjjjvyux==d,d（1.5）杆的应变为)(sin)(cosdijijvvluulll−+−==θθε（1.6）单元内的轴向力（规定拉力为正）为)](sin)([cosijijvvuulAEAEN−+−==θθε（1.7）这里A和E分别表示杆的横截面积和弹性模量。节点力与轴力的关系可表示为xyjuiujvivxjfyjfyifxifθijθθθθsin,cos;sin,cosNfNfNfNfyjxjyixi==−=−=（1.8）结合（1.7）和（1.8），可以得到力与位移间关系表示式⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡−−−−−−−−=⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡jjiiyjxjyixivuvulAEffffθθθθθθθθθθθθθθθθθθθθθθθθ22222222sinsincossinsincossincoscossincoscossinsincossinsincossincoscossincoscos（1.9）或记为如下的矩阵形式eee][][][UKF=（1.10）其中单元内的节点力矩阵、位移矩阵和单元刚度矩阵e][K如下⎥⎦⎤⎢⎣⎡=][][][jieffF，⎥⎦⎤⎢⎣⎡=][][][jieuuU，ejjjiijiie⎥⎦⎤⎢⎣⎡=][][][][][kkkkK（1.11）并且⎥⎦⎤⎢⎣⎡=22][βαβαβαlAEiik，⎥⎦⎤⎢⎣⎡−−−−==22][][βαβαβαlAEjiijkk，⎥⎦⎤⎢⎣⎡=22][βαβαβαlAEjjk，θβθαsin,cos==（1.12）我们注意到，上式中的单元刚度矩阵e][K具有对称性。在图1.1中，如果我们建立局部坐标系}{yox，x轴平行于杆的轴线，y轴垂直于杆的轴线。相应的节点位移、节点力均以上面加“–”的符号来表示（如图1.2）。图1.2杆单元以及局部坐标系示意图类似于前面的推导，可得到在局部坐标系中的节点力与位移间的关系为xyjuiujvivxjfyjfiyfjxfij⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡−−=⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡jjiijyjxiyixvuvulAEffff0000010100000101（1.13）或记为如下矩阵形式eee][][][UKF=（1.14）由图1.1和1.2，在局部坐标系和整体坐标系之间，不难写出力的转换关系θθθθcossin,sincosyixiiyyixiixffffff+−=+=（1.15a）θθθθcossin,sincosyjxjjyyjxjjxffffff+−=+=（1.15b）或写成eee][][][FTF=（1.16）其中⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡−−=θθθθθθθθcossin00sincos0000cossin00sincos][eT（1.17）为转换矩阵。对于节点位移，在局部坐标与整体坐标之间，也存在着类似的关系：eee][][][UTU=（1.18）将（1.16）和（1.18）代入（1.14）,得到eeeee][][][][][UTKFT=（1.19）由于转换矩阵e][T为正交矩阵，有Tee][][1TT=−。对上式两边同时左乘1][−eT，得到eeeTee][][][][][UTKTF=（1.20）将上式与（1.10）相比较不难看出，在整体坐标系中的单元刚度矩阵e][K可以由局部坐标系中的刚度矩阵e][K得到，即eeTee][][][][TKTK=（1.21）由此可见，在建立杆单元刚度分析过程中，只需在局部坐标系中进行，而且所有杆单元的局部坐标系中的刚度矩阵都是类似的表达式（如（1.13）式），即而后通过杆单元的转换矩阵（对应于不同的倾角θ），即可得到整体坐标中的刚度矩阵。这为处理大量的杆单元带来了便利以及便于程序化处理。（二）节点平衡方程与整体刚度矩阵在桁架结构中，往往是由多个杆件在节点处连结而成。现我们从中考察一个节点i，如图1.3所示，设环绕节点i共有三个单元，即ipimij,,，节点承受的沿整体坐标}{xoy的载荷分量记为yixiFF,。图1.3节点i的平衡根据力的平衡，作用于杆单元的节点力与作用于该节点的节点力，其大小相等而方向相反。以杆ij为例，作用于杆单元的节点力是Tyijxijijff],[][=f，而作用于节点i的节点力为Tyijxijff],[−−。为了分析节点的平衡，将节点i作为脱离体取出来（如图1.3）。节点的平衡方程为∑==pmjexiexifF,,，∑==pmjeyieyifF,,（1.22）显然与节点i无关的单元不计入上述求和。根据力与位移的关系（1.10），可知杆单元ij施加与节点i的节点力为]][[]][[][jijiiiyijxijijffukukf+=⎥⎥⎦⎤⎢⎢⎣⎡=（1.23）相应地，可以得到单元im和ip施加与节点i的节点力分别为]][[]][[][mimiiiyimximimffukukf+=⎥⎦⎤⎢⎣⎡=，]][[]][[][pipiiiyipxipipffukukf+=⎥⎥⎦⎤⎢⎢⎣⎡=（1.24）将（1.23）和（1.24）代入平衡方程（1.22）得到][][][][][,,)(epmjeieieiiiukukF∑∑=+=（1.25）xiFyiFipjmxiFyiFixijfximfxipfyijfyimfyipf其中TyixiiFF}{][=F。对于每个节点都有类似于如上的平衡方程，对于全部节点，Ni,,2,