第八章几何非线性问题的有限元法

1第八章几何非线性问题的有限元法8.1引言前面各章所讨论的问题都是在小变形假设的前提下进行的，即假定物体所发生的位移远小于物体自身的几何尺寸，应变远小于1。在此前提下，建立物体或微元体的平衡条件时可以不考虑物体的位置和形状（简称位形）的变化，因此在分析中不必区别变形前后位形的差别，且应变可用一阶无穷小的线性应变表达。实际上，上述假设有时是不成立的。即使实际应变可能是小的，且不超过材料的弹性极限，但如果需要精确地确定位移，就必须考虑几何非线性，即平衡方程应该相对于变形后的位置得出，而几何关系应该计及二次项。例如平板大挠度理论中，由于考虑了中面内的薄膜应力，求得的挠度比小挠度理论的结果有显著的减低。再如在结构稳定性问题中，当载荷达到一定数值后，挠度比线性解答予示的结果更剧烈地增加，并且确实存在承载能力随继续变形而减低的现象。在冷却塔、薄壁结构及其它比较细长的结构中，几何非线性分析都显得十分重要。几何非线性问题可以分为以下几种类型：（1）大位移小应变问题。一般工程结构所遇到的几何非线性问题大多属于这一类。例如高层建筑或高耸构筑物以及大跨度网壳等结构的分析常需要考虑到结构大位移的影响。（2）大位移大应变问题，如金属压力加工中所遇到的问题就属于这一类型。（3）结构的变形引起外载荷大小、方向或边界支承条件的变化等。结构的平衡实际上是在结构发生变形之后达到的，对于几何非线性问题来说，平衡方程必须建立在结构变形之后的状态上。为了描述结构的变形需要设置一定的参考系统。一种做法是让单元的局部坐标系始终固定在结构发生变形之前的位置，以结构变形前的原始位形作为基本的参考位形，这种分析方法称作总体的拉格朗日（Ｌagrange）列式法；另一种做法是让单元的局部坐标系跟随结构一起发生变位，分析过程中参考位形是不断被更新的，这种分析方法称作更新的拉格朗日列式法。本章首先对几何非线性问题作一般性讨论，从中导出经典的线性屈曲问题的公式；然后建立平板大挠度问题和壳体的大位移（及大转动）分析的有限方法公式；接着还给出了大应变及大位移的一般公式，最后还详细讨论了杆系结构几何非线性问题的有关公式。在讨论中我们采用总体的拉格朗日列式法，但对杆系结构，为应用方便我们给出了两种列式法的公式。8.2一般性讨论8.2.1理论基础无论是对于何种几何非线性问题，虚功原理总是成立的。由虚功原理，单元的虚功方程可以写成如下的形式0veeTeeTFdv（8.1）其中F为单元节点力向量，e为单元的虚应变，e为节点虚位移向量。增量形式的应变一位移关系可表示为eedBd（8.2）上式中ed表示单元节点位移e的微分。根据变分与微分运算在形式上的相似性，有2eeB(8.3)以上两式中B称为大位移情况下的增量应变矩阵，代表了单元应变增量与节点位移增量之间的关系。在大位移情况下B应是节点位移的函数。若将上述应变增量矩阵分解为与节点位移无关的部分0B和与节点位移有关的部分)(LB两部分组成，即LBBB0（8.4）此时0B也就是一般线性分析时的应变矩阵。将式（8.3）代入（8.1），并考虑到节点虚位移的任意性，可将单元的平衡方程写成0veeTFdvB（8.5）按照式（8.5）可以对整个结构建立有限元列式，这种列式方法可称为全量列式方式，在几何非线性分析中，按照这种列式方法得到的单元和结构刚度矩阵一般是非对称的，于求解不利。因此，在分析非线性问题时大多采用增量列式法。以下就着重介绍这一方法。式（8.5）所示的平衡方程可以写成微分的形式0)(eveTFddvBd（8.6）由于在几何非线性问题中，应变矩阵B和应力e都是节点位移的函数，因此有eTeeTdBBdBd)(（8.7）将式（8.7）代入（8.6），则有eeTvevFddvdBdvBd（8.8）单元内部的应力增量与应变增量存在确定的关系，这种关系可以用增量形式表示为eedDd（8.9）式中D称为应力一应变关系矩阵，或称为材料的本构关系矩阵。如果材料属于线性弹性的，D将是一个常数矩阵。并且，对于线性弹性材料来说有eeeeD00)(（8.10）3上式中e0和e0分别为单元材料中可能存在的初应变和初应力。将式（8.2）代入式（8.9）就可以得到应力增量与单元节点位移增量之间的关系eedBDd（8.11）将式（8.4）代入式（8.11）后得eLedBBDd)(0（8.12）于是，式（8.8）左端中的第二项便可表示为veLvTLTLvvLTTevTddvBDBdvBDBdvBDBdvBDBdvdB))((0000（8.13）若记dvBDBkTv000(8.14)0k是与单元节点位移无关，它就是一般线性分析时的单元刚度矩阵。式（8.13）右端第二层括号内的项可记为vLTLTLLTLdvBDBBDBBDBk)(00（8.15）Lk称为单元的初位移矩阵或大位移矩阵，表示单元位置的变动对单元刚度矩阵的影响。现在再来看式（8.8）左端的第一项。考虑到式（8.4）的关系并注意到0B与节点位移无关，因此对节点位移的微分等于零，对于一个确定的有限元分析模型，式（8.8）左端的第一项可一般地写成eevTLevTdkdvBddvBd（8.16）上式中k称为单元的初应力矩阵或几何刚度矩阵，它表示单元中存在的应力对单元刚度矩阵的影响。由上式（8.16）和式（8.13），并考虑到式（8.14）和（8.15）的关系，有eeLFddkkk)(0（8.17）若记LTkkkk0（8.18）Tk就称为单元的切线刚度矩阵。此时，有增量形式的单元刚度方程eeTFddk（8.19）4由此可以看出，单元切线刚度矩阵Tk代表了单元于某种变形位置时的瞬时刚度，或者说代表了单元节点力与节点位移之间的瞬时关系。有了单元切线刚度矩阵就可以按照常规的方法，即单元集成法组装结构的切线刚度矩阵，即有TTkK（8.20）并进而得到结构的增量刚度方程FddKT（8.21）前面在推导式（8.8）时，假定载荷eF与变形无关。但有些情况并非如此。例如，作用于特大变形结构上的压力载荷，与变形有关的气动载荷便是这样。在这种情况下，式（8.8）应计及载荷相对于d的微分项，本书后面的推导中均不考虑这一影响。8.2.2求解方法对于实际应用，载荷增量不可能取成微分的形式，总是一个有限值。于是，按式（8.21）求得的位移增量使结构偏移了其真实的平衡位置。为了解决这一问题，可以根据当时的结构位移情况按式（8.5）求各单元上作用的节点力，并继而求得各节点合力。然后将外载荷与上述节点合力之差，即节点的不平衡力，作为一种载荷施加于结构，由此求得节点位移的修正值。上述过程也可以反复多次。综上所述，总体的拉格朗日列式方法的一次完整的迭代步骤一般可归纳如下：（1）按线性分析得到节点位移的初值1。（2）形成局部坐标系中的单元切线刚度矩阵Tk，并按式（8.5）计算单元的节点力eF。（3）将Tk和eF转换到整体坐标系。（4）对所有单元重复（2）至（3）的步骤。生成结构的切线刚度矩阵1TK和节点力合力1F。（5）计算节点不平衡力eFF11。（6）求解结构刚度方程111TK，得节点位移增量1。（7）将1叠加到节点位移向量1中，即112。（8）收敛条件判断，如果不满足则反回到步骤（2）。上述在总载荷下进行迭代的方法有时会遇到困难。在非线性程度较高的问题中可能收敛较慢，此外，当解答非唯一时，有可能得到实际上不需要的那个解。在这种情况下，可采用7.2.4节中所介绍的增量法求解，并得到每一增量步的非线性解。如迭代中再带有自平衡校正，并采用小的载荷增量，通常一步运算就能足够精确地得到该增量步的解。以上两节所介绍的增量形式的总体拉格朗日列式法，在结构的非线性分析中应用十分广5泛，有关计算公式及求解方法对板、壳或杆件体系的非线性分析都同样适用。由上面的分析也可以看出，采用总体的拉格朗日方法求解非线性问题的关键是形成单元的切线刚度矩阵。8.3屈曲问题非线性分析，尤其是几何非线性分析在很多情况下是估算一个结构在失去稳定性前所能承受的最大载荷。这是结构屈曲问题的研究目标，是固体力学的一个重要分支，也是工程实践中经常出现的问题。小位移线性理论假设在结构受载变形过程中忽略了结构的位移变化，因此在加载的各个阶段总是认为结构在未加载的原始位形上产生平衡，当屈曲发生时，结构位形突然跳到另一个平衡位置。图8.1(a)为线性屈曲的示意图。为裁荷比例因子，其含义稍后会讲到，它与位移在屈曲前为线性关系，当载荷达到极限值（图中分枝点）时结构失稳，曲线改变，结构平衡转向另一模态。这就是线性屈曲也称分枝屈曲。严格说来，结构的平衡实际上是在结构发生变形之后达到的，因此，从加载一开始就出现了几何非线性的特性，图8.1(b)为非线性屈曲的示意图，当载荷比例因子增加时，曲线是非线性的，一直达到极限，这种在结构发生变形一直到失稳，在变形后的位形上考虑平衡一直达到极限的方法称非线性屈曲或极限屈曲。图8.1a图8.1b可见，工程实际中分枝屈曲现象实为罕见，它仅出现在完全无结构缺陷，完全沿轴向加压的绝对直杆及完整空球壳在均匀外压的情况下。分枝屈曲现象虽然罕见，但实际中有不少结构屈曲状态接近分枝屈曲，而分枝屈曲的计算工作量又远小于计算极限屈曲的工作量，况且，不少作者得出结论，一些中等非线性的屈曲状态，可以用线性屈曲问题特征向量的线性组合近似得到。因此线性屈曲理论还是有其实际价值。屈曲的含义可简述为：结构处于一种平衡状态，载荷增量为一个微量，其位移增量很大。用方程来表达这种物理现象，则由总体拉格朗日列式法建立的结构刚度方程（8.21）变成为0dKT（8.22）根据式（8.18）和（8.20），有LTKKKK0（8.23）将式（8.23）代入（8.22）得0)(0dKKKL（8.24）在线性屈曲情况下，屈曲前结构处于原始位形的线性平衡状态，因此上式中的大位移矩阵LK应为零，此时式（8.24）简化为00dKK（8.25）6由式（8.16）可以看出，K并不明显地包含位移增量d，在小变形情况下，该矩阵与应力水平成正比。由于屈曲前线性假设，多数情况下，应力与外载也为线性关系，因此，如果令某一参考载荷rF对应的初应力刚度阵为rK，令屈曲极值载荷为cF，与参考裁荷有rccFF关系，称载荷比例因子，c为极值裁荷的比例因子。极值载荷时的初应力刚度阵为rccKK（8.26）代入式（8.25），则有0)(0dKKrc（8.27）可以看出式（8.27）是一个广义特征值方程，也是经典弹性稳定理论的最后控制方程。实际求解时可按以下步骤进行：（1）按线弹性问题的有限元法形成各单元的刚度矩阵0k，并用常规方法组装成结构刚度矩阵，即00kK。（2）对结构施加参考载荷rF，并求解有限元方程rFK0，进而可求得各单元节点应力。（3）对有膜应力存在的单元按式（8.16）形式初应力矩阵