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

1、第八章几何非线性问题的有限元法8.1 引言前面各章所讨论的问题都是在小变形假设的前提下进行的，即假定物体所发生的位移远小于物体自身的几何尺寸，应变远小于1。在此前提下，建立物体或微元体的平衡条件时可以不考虑物体的位置和形状(简称位形)的变化，因此在分析中不必区别变形前后位形的差 别，且应变可用一阶无穷小的线性应变表达。实际上，上述假设有时是不成立的。即使实际应变可能是小的， 且不超过材料的弹性极限，但如果需要精确地确定位移，就必须考虑几何非线性，即平衡方程应该相对于变形后的位置得出，而几何关系应该计及二次项。例如平板大挠度理论中， 由于考虑了中面内的薄膜应力，求得的挠度比小挠度理论的结果有显著

2、的减低。再如在结构稳定性问题中，当载荷达到一定数值后，挠度比线性解答予示的结果更剧烈地增加，并且确实存在承载能力随继续变形而减低的现象。在冷却塔、薄壁结构及其它比较细长的结构中，几何非线性分析都显得十分重要。几何非线性问题可以分为以下几种类型：(1)大位移小应变问题。一般工程结构所遇到的几何非线性问题大多属于这一类。例 如高层建筑或高耸构筑物以及大跨度网壳等结构的分析常需要考虑到结构大位移的影响。(2)大位移大应变问题，如金属压力加工中所遇到的问题就属于这一类型。(3)结构的变形引起外载荷大小、方向或边界支承条件的变化等。结构的平衡实际上是在结构发生变形之后达到的，对于几何非线性问题来说，平衡

3、方程必须建立在结构变形之后的状态上。为了描述结构的变形需要设置一定的参考系统。一种做法是让单元的局部坐标系始终固定在结构发生变形之前的位置，以结构变形前的原始位形作为基本的参考位形，这种分析方法称作总体的拉格朗日(Lagrange)列式法；另一种做法是让单元的局部坐标系跟随结构一起发生变位，分析过程中参考位形是不断被更新的，这种分析方法称作更新的拉格朗日列式法。本章首先对几何非线性问题作一般性讨论，从中导出经典的线性屈曲问题的公式；然后建立平板大挠度问题和壳体的大位移(及大转动)分析的有限方法公式；接着还给出了大应变及大位移的一般公式，最后还详细讨论了杆系结构几何非线性问题的有关公式。在讨论中

4、我们采用总体的拉格朗日列式法， 但对杆系结构，为应用方便我们给出了两种列式法的公式。8.2 一般性讨论8.2.1 理论基础虚功原理总是成立的。 由虚功原理，单元的虚功方程无论是对于何种几何非线性问题, 可以写成如下的形式eT eeT edvF 0(8.1)ve其中F为单元节点力向量，为单元的虚应变,e为节点虚位移向量。增量形式的应变一位移关系可表示为上式中d e表示单元节点位移e的微分。根据变分与微分运算在形式上的相似性，有(8.3)以上两式中B称为大位移情况下的增量应变矩阵，代表了单元应变增量与节点位移增量之间的关系。在大位移情况下B应是节点位移的函数。若将上述应变增量矩阵分解为与节点位移无

5、关的部分Bo和与节点位移有关的部分Bl（）两部分组成，即(8.4)B BoBl此时Bo也就是一般线性分析时的应变矩阵。将式（8.3）代入（8.1）,并考虑到节点虚位移的任意性，可将单元的平衡方程写T(8.5)B dv F 0v按照式（8.5）可以对整个结构建立有限元列式，这种列式方法可称为全量列式方式，在几何非线性分析中，按照这种列式方法得到的单元和结构刚度矩阵一般是非对称的，于求解不利。因此，在分析非线性问题时大多采用增量列式法。以下就着重介绍这一方法。式（8.5）所示的平衡方程可以写成微分的形式Tee八d（ B ）dv d F 0（8.6）v由于在几何非线性问题中，应变矩阵 B和应力e都是

6、节点位移的函数，因此有-T e二 e - T ed( B ) d B B d(8.7)将式(8.7)代入(8.6),则有eT . eed B dv B d dv d F(8.8)vv单元内部的应力增量与应变增量存在确定的关系，这种关系可以用增量形式表示为(8.9)式中D称为应力一应变关系矩阵，或称为材料的本构关系矩阵。如果材料属于线性弹性的,D将是一个常数矩阵。并且，对于线性弹性材料来说有D(o e)(8.10)上式中 0e和 0e分别为单元材料中可能存在的初应变和初应力。将式（8.2）代入式（8.9）就可以得到应力增量与单元节点位移增量之间的关系eed D Bd(8.(11)将式（8.4）代

7、入式（8.11）后得e _d D ( BoBl )d(8.(12)日ZE,式（8.8）左端中的第二项便可表示为e，一dv (BovBL TvT D Bo dvB0 dv(BovBL TvT D BL dvD BL dv)d若记koTBoB0 dvk0是与单元节点位移无关，它就是一般线性分析时的单元刚度矩阵。式（ 层括号内的项可记为8.13)(8.(13)(8.14)右端第二kL ( Bo T D BLvBl T D B0 BL T D BL )dv(8.(15)kL称为单元的初位移矩阵或大位移矩阵,表示单元位置的变动对单元刚度矩阵的影响。现在再来看式（8.8）左端的第一项。考虑到式（8.4）的

8、关系并注意到Bo与节点位移无关，因此对节点位移的微分等于零，对于一个确定的有限元分析模型，式（ 第一项可一般地写成8.8）左端的T ed B dv d BL vveedv k d(8.(16)上式中k称为单元的初应力矩阵或几何刚度矩阵,它表示单元中存在的应力对单元刚度矩阵的影响。由上式（8.16）和式（8.13）,并考虑到式（16.14） （8.15）的关系，有若记(kok kL )d(8.(17)kTkokkL(8.(18)kT就称为单元的切线刚度矩阵。此时，有增量形式的单元刚度方程kT d(8.19)由此可以看出，单元切线刚度矩阵kT 代表了单元于某种变形位置时的瞬时刚度，或者说代表了单元

9、节点力与节点位移之间的瞬时关系。有了单元切线刚度矩阵就可以按照常规的方法，即单元集成法组装结构的切线刚度矩阵，即有KTkT（ 8.20）并进而得到结构的增量刚度方程KT dd F（ 8.21 ）前面在推导式（8.8）时，假定载荷F e与变形无关。但有些情况并非如此。例如，作用于特大变形结构上的压力载荷，与变形有关的气动载荷便是这样。在这种情况下，式 （ 8.8）应计及载荷相对于d 的微分项，本书后面的推导中均不考虑这一影响。8.2.2 求解方法对于实际应用，载荷增量不可能取成微分的形式，总是一个有限值。于是， 按式 （ 8.21 ）求得的位移增量使结构偏移了其真实的平衡位置。为了解决这一问题，

10、可以根据当时的结构位移情况按式（8.5）求各单元上作用的节点力，并继而求得各节点合力。然后将外载荷与上述节点合力之差，即节点的不平衡力，作为一种载荷施加于结构，由此求得节点位移的修正值。上述过程也可以反复多次。综上所述，总体的拉格朗日列式方法的一次完整的迭代步骤一般可归纳如下：（ 1）按线性分析得到节点位移的初值1 。（ 2） 形成局部坐标系中的单元切线刚度矩阵kT ， 并按式 （ 8.5） 计算单元的节点力F e 。e（ 3）将kT 和 F 转换到整体坐标系。（ 4）对所有单元重复（2）至（3）的步骤。生成结构的切线刚度矩阵KT 1和节点力合力 F 1 。5）计算节点不平衡力1 F 11。6

11、）求解结构刚度方程K T 111 ，得节点位移增量（7）将1 叠加到节点位移向量1 中，即211 。（8）收敛条件判断，如果不满足则反回到步骤（2） 。上述在总载荷下进行迭代的方法有时会遇到困难。在非线性程度较高的问题中可能收敛较慢，此外，当解答非唯一时，有可能得到实际上不需要的那个解。在这种情况下，可采用7.2.4 节中所介绍的增量法求解，并得到每一增量步的非线性解。如迭代中再带有自平衡校正，并采用小的载荷增量，通常一步运算就能足够精确地得到该增量步的解。以上两节所介绍的增量形式的总体拉格朗日列式法，在结构的非线性分析中应用十分广泛，有关计算公式及求解方法对板、壳或杆件体系的非线性分析都同样

12、适用。由上面的分析也可以看出，采用总体的拉格朗日方法求解非线性问题的关键是形成单元的切线刚度矩阵。8.3屈曲问题非线性分析，尤其是几何非线性分析在很多情况下是估算一个结构在失去稳定性前所能 承受的最大载荷。这是结构屈曲问题的研究目标，是固体力学的一个重要分支，也是工程实践中经常出现的问题。小位移线性理论假设在结构受载变形过程中忽略了结构的位移变化，因此在加载的各个阶段总是认为结构在未加载的原始位形上产生平衡，当屈曲发生时，结构位形突然跳到另一个平衡位置。图8.1（a）为线性屈曲的示意图。为裁荷比例因子，其含义稍后会讲到，它与位移 在屈曲前为线性关系，当载荷达到极限值（图中分枝点）时结构失稳，曲

13、线改变，结构平衡转向另一模态。这就是线性屈曲也称分枝屈曲。严格说来，结构的平衡实际上是在结构发生变形之后达到的，因此，从加载一开始就出现了几何非线性的特性， 图8.1（b）为非线性屈曲的示意图， 当载荷比例因子增加时，曲线是非线性的，一直达到极限，这种在结构发生变形一直到失稳，在变形后的位形上考虑平衡一直达到极限的方法称非线性屈曲或极限屈曲。图 8.1a图 8.1b可见，工程实际中分枝屈曲现象实为罕见， 压的绝对直杆及完整空球壳在均匀外压的情况下。它仅出现在完全无结构缺陷，完全沿轴向加分枝屈曲现象虽然罕见，但实际中有不少结构屈曲状态接近分枝屈曲，而分枝屈曲的计算工作量又远小于计算极限屈曲的工作

14、量，况且，不少作者得出结论， 一些中等非线性的屈曲状态，可以用线性屈曲问题特征向量的线性组合近似得到。因此线性屈曲理论还是有其实际价值。屈曲的含义可简述为：结构处于一种平衡状态， 载荷增量为一个微量，其位移增量很大。用方程来表达这种物理现象，则由总体拉格朗日列式法建立的结构刚度方程（8.21）变成为KT d0（8.22）根据式（8.18）和（8.20）,有KTK0 K KL（8.23）将式（8.23）代入（8.22）得（K0 K KL ）d 0（8.24）在线性屈曲情况下，屈曲前结构处于原始位形的线性平衡状态，因此上式中的大位移矩K。阵KL应为零，此时式（8.24）简化为(8.25)8.16）

15、可以看出，K 并不明显地包含位移增量d ，在小变形情况下，该矩阵与应力水平成正比。由于屈曲前线性假设，多数情况下，应力与外载也为线性关系，因此，如果令某一参考载荷F r 对应的初应力刚度阵为K r ，令屈曲极值载荷为F c ，与参考裁荷有Fc cFr关系，称载荷比例因子，c为极值裁荷的比例因子。极值载荷时的初应力刚度阵为Kc c Kr（ 8.26）代入式（8.25） ，则有（ K0 c Kr ）d 0（ 8.27）可以看出式（8.27）是一个广义特征值方程，也是经典弹性稳定理论的最后控制方程。实际求解时可按以下步骤进行：（ 1）按线弹性问题的有限元法形成各单元的刚度矩阵k0 ，并用常规方法组装

16、成结构刚度矩阵，即K0k0 。（ 2）对结构施加参考载荷Fr ，并求解有限元方程K0Fr ，进而可求得各单元节点应力。（ 3）对有膜应力存在的单元按式（8.16）形式初应力矩阵kr ，并组装成结构初应力矩阵，即K rkr 。（4）求解广义特征值问题，即式（8.27）,得到最低几阶特征值c及对应的特征模态。理论上它们都是平衡模态，当载荷达到分枝载荷cFr时，平衡由一种模态“跳列”另一种模态，这也就是分枝屈曲的含义。（ 5） 从分枝载荷中挑选一个最小的载荷，即 Fmrincmin F r ， 作为极值载荷或临界载荷。值得指出的是，由式（8.27）所表征的线性屈曲问题，是建立在下述假设的基础上的，即

17、假设线性应变刚度矩阵在屈曲前不产生明显的变化，且初应力矩阵简单地与应力水平成正比。 如前所述，这在实际问题中是很罕见的。由式 （ 8.27） 所确定的极值载荷只能是近似的。由于线性屈曲理论存在精度差，以及适用范围窄的限制，所以，在一般情况下，应当用切线刚度矩阵K T 来研究这个问题。当K T d 0 时，发生随遇平衡。显然，这里应该用逐次逼近的方法进行求解。关于非线性屈曲问题的有限元解法，读者可参考其它文献。8.4 板的大挠度及线性屈曲8.4.1 基本问题在板的大挠度问题中，板中的内力除弯曲内力外，还有薄膜内力，如图8.2 所示，在这种情况下，横向位移会引起薄膜应力，因此在板的大挠度问题中，面

18、内变形和弯曲变形不再象小挠度板那样能分别处理，二者是相互耦合的。图8.2图8.3和以前一样，平板的应变可用中面的位移来描述。 义应变和广义应力向量可以表示为如果令y平面与中面2w2x2 TW zx y(8.28)phNxbxNyN2yMxMyMxyT(8.29)式中上标p和b分别表木面内和弯曲的分量。x和y方向产生附加伸长和附加角变形,因此广义由图8.3可以看出，挠度w使中面在: 应变可表达成y x2w2x2w2-y22wx y-(2wW)2 X w)2 y wp 0b0(8.30)式中第一项是我们熟知的线性项，而第二项则是非线性项。如果所考虑的材料是线弹性的，则平板的弹性矩阵由平面应力弹性矩

19、阵Dp和弯曲弹性矩阵Db组成，即D Dp0Db(8.31 )平板单元的位移可以采用适当的形函数，通过节点位移来表达，写成(8.32)为方便起见，将节点位移向量也分为面内的和弯曲的两部分，即(8.33)UiTViWi(-)ixT(-)iy这样，形函数也应分成通过上述规定，除了非线性应变项NiNip00bNi(8.34)以外,其余各项都和标准的线性分析相同，不再重复。8.4.2 应变矢I阵B的计算为了确定单元的切线刚度矩阵我们先来讨论单元的应变矩阵B的计算。首先应注意到式中BoB0p0BBo0BbBlBl0Bb00这里，B0p和Bo分别是平面单元和弯曲单元按线性分析所得到的标准矩阵，而变的非线性项

20、引起的，它可以通过L对参数 b取微分来导出。下面就来推导达式。在式（8.30）中，非线性应变分量可以方便地写成0 x1 cW02 ywwyx(8.35)Bb是应Bl的表1A 2(8.36)式中(8.37)wx wyew的一阶导数可以用弯曲节点位移b表示成(8.38)其中wxwy(8.39)bNibN2GNixbN2yy(8.40)可见G只与单元的节点坐标有关。将式(8.36)取微分，有c11dL-dAAd(8.41)22注意到矩阵 A和具有两个有用的性质：1.设xXi,X2 T是一个任意向量，则有wd(一)0x w0d(一)ywwd(-) d(-)yxXiX2Xi0d(-w)Xi00X2d(-

21、xw、0x2 dX2Xi)X2Xi(8.42)于是有d AAd(8.43)2.设yH，V3个任意向量，则有d AT yw d(一) xd(利用上述性质wd( 一 )yd(w) y w) xyiy2y3i,并注意到式（8.39）,可将式（yiy38.41)y3y2写成LpAd观察上式，由应变矩阵的定义，可见Bb A G8.4.3单元切线刚度矩阵 片的计算根据式（8.18）,单元的切线刚度矩阵由三部分组成，即kTk0 k kL由第二章和第四章给出的刚度矩阵,k0将式（8.35）中的B0和Bl代入式0B°p T最后，利用式（8.16）w d(一) xw d() y可将线性小变形的刚度矩阵写

22、成下式k°p00kb8.15),Dp Bb可求出初应力矩阵d Bl可得大位移矩阵对称b Tb dvBb Dp Bbk 。为此，将式（8.35）中的0T d Bb把上式代入式（8.16）,并利用式（8.29)和(8.46),给出利用式（8.44）,有0G Td ATNxNyNxydvMxM yMxyyiy3Bly3dy2(8.44)(8.45)(8.46)(8.47)(8.48)(8.49)进行微分,Nx d AT NyNxyNxNxyNxyNyNxNxyNxyNybeG d b最后可将初应力矩阵写成式中kb0 kb(8.50)NxNxyNxyNyG dv(8.51 )是平板对称形式的

23、初应力矩阵。（8.47）,便得到板的切线刚度矩阵,其余计算步骤已如8.2.2将式（8.48）至（8.50）代入 节所述，不再重复。8.4.4板的屈曲问题板的线性屈曲问题，可从式（8.25）出发作为特征值问题考虑。对于只在平面内承受压力载荷的情况，先求问题的弹性解，然后用式（8.51）求出kb ,而由线弹性方法求得 kb ,此时kL0 ,在完成结构刚度矩阵的组装后得到了广义特征值问题（K； c Kb ）d 0（8.52）解之，可得平板线性屈曲问题的解。上述平板的几何非线性有限元理论，也适用于由平板单元组成的壳体结构的几何非线性 分析，只是多了一步将局部坐标系中的切线刚度矩阵及其它有关量转换到整体

24、坐标系中的工 作。值得指出的是，通常对壳体结构进行线性屈曲的分析结果远远大于实际的失稳载荷。因此，确定变形对屈曲的影响是十分重要的，也就是说，为了得到正确的结果，必须进行完全 的非线性分析。8.5三维单元的大应变和大位移公式上节平板大挠度问题中，非线性应变一位移关系是在特定的情况下建立的。根据格林 （Green）应变的定义，现在可以从一般方法出发推导出几何非线性理论的有限元计算公式， 而不管应变或位移是大的还是小的。用直角坐标表示的格林应变为xux12(u)2 xv 2 ()x(-w 2 )xuvu uvvw wxyyxx yxyx y其它四个应变只要通过轮换xyzx,uvwu ，即可得到。在

25、小位移的情况下，忽略二次项，就得到线性应变公式。这时,(8.53)x, y和Z是原来平行于坐标轴的线段的伸长率，而 xy, yz和7*表示原来平行于坐标轴的两垂直线段夹角的变化。但在应变较大时，上述几何意义就不再成立。现在我们来推导三维应力状态下非线性的B和的一般计算公式。稍加改变即可得到一维和二维情况下的计算公式。对于板壳问题，可在特定条件下忽略某些项后就可导出上节中的结果。对杆系结构我们将在下节详细讨论。8.5.1应变矢I阵BL的推导变形体某点的应变应是线性应变0和非线性应变L之和,式中yzzxxyvwzywuxzuvyx(8.(54)(8.(55)按照式（8.53）,非线性应变可用下式表

26、示式中0T zT yA是一个6X9的矩阵。是一个0T zT y0T y 0T z 0T x9X 10T y 0T z 0T x0T zT yT x0(8.56)0T zT yT x0的列向量，而余类推。容易验证上述定义的正确性，并重新建立上节中矩阵A和的两个性质。我们再次有d l -d 2AdAd(8.57)如果 用形函数N和节点位移向量则可写出(8.(58)式中G的表达式和上节中的式（8.40）类似，为形函数对坐标的一阶导数所组成，因而只Bl8.5.2单元切线刚度矩阵%的推导注意到式（8.4）,即BBoBl8.14）积分得到k0 ,而由式（8.15）可与节点坐标有关。将上式代入式（8.57）

27、得因此由应变矩阵的定义即有求得kL ,又由式（8.16）,则有d Bl Tvedv(8.60)由式（8.59）BlT ATe 注意到G不含 ，由上式得d BlG Td AT(8.61)将式（8.61）代入（8.60）的右端，则同上一节一样，利用 A和的第二个性质,G Td ATedv(8.(62)d ATxIxy I xzIyIyzIzI(8.(63)式中xI对M xyIyI称xzIyzIzI由式（8.58）可得ed G d（8.64）将式（8.64）代入（8.63）则d ATe M G d e（8.65）再将式（8.65）代入（8.62）,使得到具有对称形式的初应力矩阵k G T M Gdv

28、（8.66）v至此，三维单元的切线刚度矩阵kT就可由下式求得kTko 。 k8.6杆系结构的大位移分析对杆系结构进行大位移分析时，可采用总体的拉格朗日列式法或更新的拉格朗日列式 法。对总体的拉格朗日列式法，我们将推导出桁杆单元和梁单元的切线刚度矩阵及其显式， 其余计算均如8.2.2节所述。而对于更新的拉格朗日列式法，则以刚架单元为例，重点介绍 其计算原理和实施步骤。8.6.1 桁架单元的切线刚度矩阵考虑一个横截而积为 A,弹性模量为巳长度为L的桁杆单元，它在发生变位前后的位 置如图8.4所示。图8.4在小位移的情况下，桁杆单元上某一点的轴向应变为(8.67)dux dx如果节点的位移比较大，则

29、由于横向位移v会使单元发生附加的轴向伸缩。这种附加的轴向应变与单元在变位过程中所转过的角度有关，此时单元的轴向应变可表示为(8.68)du 1/dv（一）dx 2 dx式中左端的第二项是考虑大位移时附加的非线性项。显然上式也可由格林应变公式直接写 出。12x在大位移分析中，桁杆单元的轴向和横向位移函数可精确地取x的一次函数，即取(8.(69)34x同以前线性分析时的单元分析过程一样,由杆两端节点位移条件可解出单元上任意点的位移可用杆端节点位移表示为(8.(70)式中为桁杆单元的形函数矩阵。将式（8.70）代入式L(8.68)xL再考虑到式8.71),有e)2(8.(71)(8.(72)10-0

30、L(8.73)将上式与式（8.2)比较,并考虑到式8.4),有在式（8.75）Bo(8.74)Bl(8.75)为单元在变位过程中发生的转角,如图8.4所示,1L (VjVi)(8.76)于是式（8.75）可以写成Bl(8.77)将式（8.74）代入式（8.14）,可得线性分析时桁杆单元的刚度矩阵，如式（5.6）所示。将式（8.74）和（8.77）代入式（8.15）可得桁杆单元的大位移矩阵如下22kLEAL 0(8.78)由式（8.78）可以看出，当单元处于变位以前的原始位置时,0,因而大位移矩阵kL0。卜面来推导桁杆单元的初应力矩阵。由式（8.77）8.76）的关系，有Bl T01L 1 /、

31、10 L(vj vi)L21Luiviujvj(8.79)式（8.79）与点有坐标无关，而对于桁杆单元来说,个常数，因此在积分时这些项均可提到积分号之外。若记N A为桁杆单元的轴向力，则有TN 0101d BldvdL 0 00001 01将上式与式（8.16）比较可得0000N0101k (8.80)L000001 01这就是桁杆单元的初应力矩阵，它计及了桁杆单元的轴向力对杆端横向位移的影响。最后将式（5.6）, （8.78）和（8.80）代入式（8.18）,便得到桁杆单元的切线刚度矩阵。另外应注意到，由于k0 , kL和k均为对称矩阵，所以桁杆的切线刚度矩阵kT也是对称矩阵，由kT所组装的

32、结构刚度矩阵KT也是对称的，这对求解是很有利的。8.6.2刚架单元的切线刚度矩阵在大位移的情况下，仍可假定刚架单元的位移函数与小位移情况下的相同，即单元的轴向位移u是局部坐标x的线函数，而横而位移 v则是x的三次函数。用式子表示如下(8.81 )u12x23v34乂5乂 6乂由杆两端的节点位移条件，可从上式解得1 6 ,于是单元上任意点的位移可用杆端位移表不成式中23x(1 Lx)23(：)2若记N1N23(X)2X 0 L2(：)3x(1：)23(-x)22(x)3而式可表不成(8.(82) 成(8.(83)(8.(84)3/X Xx(l 1)( l)(8.(85)N1N2(8.85)du1

33、 ,dv、2dua-()2vx2 dxdxbd2vd 2vy 2y 2dxdx式中a代表由线位移所引起的轴向应变,1(dX)20 l (8.86)0b代表由弯曲所引起的轴向应变。而N1N2假若单元的曲率仍能用横向位移V的二阶导数近似地表示，则代表弯曲所引起的轴向应变的项仍与线性分析时相同。而由于横向位移引起的单元轴向应仍应按式(8.68)计算，此时有du(8.87)dX d2v y 2 dx为线性分析时的应变项。(8.88)为非线性项。由式(8.85)和(8.87)得d Ni dx d2 N2y2-dx将上式取微分，并注意到N1和N2不含e,有由式（8.85）和（8.88）可得将上式取微分，有

34、式中d N2 s dx由式（8.86）可得将式（8.89）将上式与式（d N1 dx d d2N2 dy2-dx(8.89)N2dx1 ,d(N22( dx0ed N2 dx(8.90)6x6x214x3x206x6x23x22xL27rLL2L2L3L2L和（8.90）代入上式,(8.91 )得d N1dxd y一8.2）比较，并注意到式2N2 dx2 (8.4)B。的关系有d N1 dxBld2 Nydx2se2(8.92)(8.93)0将式（8.84）代入（8.92）,得Bo6 小2X2- 3X八 6 2X2 03x、2(1 )y (2 )y 02(1 一)y 一(1 一)yL2LL L

35、L2 LLL(8.94)k0 ,如式(5.37)所示。将式(8.93)将B0代入式(8.14)就得到线性分析时的单元刚度矩阵8.94) 代入式8.15), 就可以得到单元的大位移矩阵0 A BA00EA L B00kL 一八一 一L 00 A BA00C000 AA0B00AA0C0C00dx LC00(8.95)式中0020Al0AB00020A0aC000-2AB 0 A2B 0 AB000-2AB 0 ABC 0 AC0ACBCdx 0AC2C2 .AA1 (vjvi)A1B1 iA1c1jBABiVvi)Bi2 iB1C1jCA1C1 (vjvi)B1C1iC12j(8.96)AiBi

36、Ci6xL21丝L2xL6x2"Lr3x27r3x21r(8.97)i, j分别为单元i端和j端的截面转角。以及L 2 .6L 2 .2 1L 2 .2 1A dxn BdxLC1 dxL015L01150115L1L1L1n ABdxB1C1dxLn AGdx010030010在计算式(8.95)中的积分时，可以应用下列积分公式(8.98)L 472L 42 1L 42 1Adx3B14 dx一 LC14dx一 L035L30'350135L 39L 39L31n A3B1dx2n A3Gdx-2n A B；dx035L035L0140L 3LL3LL31n B3GdxB&

37、#163;3dxn AC；dx014001400140L 2_2o A2B2dx3L _ 22B2C12dxLL 22o A2C12dx3035L0210035LL 2A2BC1dx 0L2AB12c1dx1L2. ABC12dx1001400140k 。对刚架单元来说，轴向应力由两部分组成，可表示为(8.99)以下我们来推导初应力矩阵(8.100)其中a代表单元由于拉伸或压缩而引起的轴向应力，在单元上任意点都相同，可以将这一常数应力用表示。b代表由弯曲变形而引起的轴向应力。由式（8.93）可知BL中第二行元素全部为零，这样，在式（8.16）中b将不发生作用，这说明单元的初应力矩阵与它的弯曲应

38、力无关。若记BL中的第一行为一e -BL1S S注意到矩阵 S中的项只是x的函数，由式（8.16）得T ,LT ,aL-T, T-T,d BLdv A d BL1dx A S d S dxI o L1ov将式（8.91）代入上式，引入式（8.97）的关系，同时注意到 A N为单元的轴向力，则 有Td BlvdvAiB0eTAi B 0dx将式(8.98)N0a2ABi 0 A2AGaCiABiB20A1B1BiCi0a2AB0A2AGAGBiCi0AGCi2k d中的积分公式代入上式便得到单元的初应力矩阵,36dx d ek也0 30L 03L3L4L2363L3LL2(8.101)363L3

39、LL2363L3L4L2最后将式(5.37), (8.95)和(8.101)代入式(8.18)，便得到了刚架单元的切线刚度矩阵。同桁杆单元一样，由于 k0 , kL和k 均为对称矩阵，所以刚架单元的切线刚度矩阵也 是对称的。8.6.3更新的拉格朗日列式法前面我们在推导桁杆单元和刚架单元的切线刚度矩阵时，采用了总体的拉格朗日列式法。其特点是，所有静力学和运动学方面的变量，例如单元刚度矩阵， 单元节点位移和单元 节点应力等，都是以t=0时刻的位形，即变形时的位形为参照系统。这里t 0, t,2 t,只是用来描述物体位形变化发展过程的一个标志量，而并不一定具有真正的时间意义。上述各量，即单元的刚度矩

40、阵，节点位移和节点应力等量的描述和计算，也可以用物体变形过程中某一时刻 t的位形为参照系统，进而推算 t t时刻时物体的位形。由于 t时刻 的位形和坐标值随计算而变化，所以称为更新的拉格朗日列式法。由此产生的非线性分析方法可称作带流动坐标(或变更坐标)的迭代法。这一方法对于杆系结构的大位移分析特别显 示出其优越性。尤其是在杆件发生比较大的转动时，采用这一方法更为适宜。 这是因为通过变更局部坐标系可以方便地描述单元的刚体转动，从而较容易地确定变位后的单元在变形后的结构中所发挥的作用。我们知道结构坐标系中的单元杆端力是由结构坐标系中的单元刚度矩阵与该单元杆端 位移向量的乘积得到的。 而结构坐标系中

41、的单元刚度矩阵不仅取决于单元本身的属性，还与该单元所处的方位有关。 在线性的位移分析中，由于节点位移引起的单元方位的变化十分微小而可以忽略，此时在计算结构的位移和节点合力时仍然能利用变形前方位的单元刚度矩 阵。如果结构的位移比较大，或者说在非线性问题中，单元方位的这种变化就不能忽略，此 时在计算节点合力时就应该使用变形后方位的单元刚度矩阵。在进行结构的大位移分析时， 可以将按线性分析所得到的节点位移作为结构位移的第一 次近似值。根据上述节点位移可以对单元刚度矩阵进行修改，从而反映单元在变位后的位置上所起的作用。在带有流动坐标的迭代法中所采用的做法是：根据已求的节点位移修改单元的坐标转换矩阵 T

42、 ,从而达到修正结构坐标系中单元刚度矩阵的目的。根据修正后的各单元刚度矩阵乃至刚度方程，可以计算出节点合力。按照上述结构位移的第一项近似值算出的 节点合力与节点所受到的外载荷并不相等，也就是说此时节点的平衡条件未被满足。这是因为按线性分析所得到的节点位移并不代表结构真正的平衡位置。于是，在这样的位置上结构当然也就无法保持平衡。现将原结构的等效节点载荷与上述节点合力之差称为节点的不平衡 力。为了求得结构真正的平衡位置，可以将不平衡力作为一组新的外载荷施加于上述已发生 变形的结构上，求得节点位移的修正值又可以重新修改各单元的坐标转换矩阵，并进而得到新的节点合力和节点不平衡力，继而再将新的节点不平衡

43、力施加于变形以后的结构。重复上述过程一般总可以使节点不平衡力减小到可以被忽略的水平，此时的节点位移所对应的便是结构在发生大位移之后真正的平衡位置。按照上述平衡位置可以计算结构在大位移情况下的杆件内力。以上就是带有流动坐标的迭代法分析大位移问题的基本思路。以下讨论流动坐标系中的杆端位移和杆端力的计算方法。图8.5 (a)给出了在结构坐标系xy中的一个未变形的刚架单元，它在弯形后的位形如图8.5b所示。xy是该单元的流动(a)坐标，它的坐标原点位于变形后的杆端，x轴沿变形后杆端节点的连线方向。(b)图8.5对比图8.5a、b可得XiXoUj Ui yi yo Vj Vi(8.(102)(8.(10

44、3), yiarctan 一Xi于是，刚架单元在流动坐标系 x y中的节点位移可表示为1''''11/ 22.21UiViVjUj l lo (xiyi) lo''i i (0 )j j (0)因此，刚架单元在 x y中的节点位移向量为T(8.104)8.105)0 0 i uj 0 j此时在上述局部坐标系中的杆端力可表示为（Fk其中k'为局部坐标系x'y'中的单元刚度矩阵，仍可按式（5.37）计算。由于采用了流动坐标系，方向角就成为杆端位移的函数，可由式（8.102）计算。这样，单元的坐标转换矩阵T 和通过 T 作用于

45、k' 得到的结构坐标系中的单元刚度矩阵k也都成为杆端位移的函数。综上所述，在采用带有流动坐标的迭代法时，一个典型的迭代过程包括以下步骤。（ 1）为外载荷作用下的结构假定一组节点位移。（ 2）根据结构坐标下的节点位移向量，确定单元两端的位置，建立单元的局部坐标系。（ 3）计算上述局部坐标系中的杆端位移向量' 。（ 4）形成局部坐标系中的单元刚度矩阵k ' 和杆端力向量F ' 。，一 ' 一 一. （ 5）将 k 和 F 转换到结构坐标系，得到k 和 F 。（ 6）对所有单元重复（2）至（5）的步骤。生成结构刚度矩阵K k 和节点合力F。（ 7）计算节点不平

46、衡力R R F 。（ 8）求解结构方程KR 得节点位移增量。（ 9）将叠加到节点位移向量中。（ 10）收敛条件判断，如果不满足则返回到步骤（2） 。以上的过程可概括为如下的迭代公式却是一个值得研究的问8.106)这样通Ki i1 R ki iII ii1般来说作第一次迭代运算时节点位移常有利于加快收敛速度。但当结构的非线性程度很高时如何假定 题。例如对于图8.6所示结构，按线性分析方法无法确定结构的位移或者说结构的位移将成为无穷大，此时就需要为载荷作用点C先假定一个适当大小的竖向位移，然后才能开始迭代过程。而在某些情况下需要把外载荷分成数级逐步施加到结构上。图8.6迭代过程的收敛可根据上述不平

47、衡力进行判断，当不平衡力和外载荷的比率减小到一个给定的限度时则可以认为迭代过程已达到了收敛，这样的收敛条件称为力收敛条件。 在结构的大位移分析中一般采用位移收敛率条件效果更好一些，因此这种收敛条件的应用更为普遍。位移收敛条件可由不同的形式提出，若记Nn J T n(8.107)式中 n是经过n次迭代运算得到的节点位移向量，则位移收敛条件的一种形式是Nn Nn 1Ner(8.108)这里er是精度要求，可以根据工程的要求和问题的性质而定，一般可取10-6到10-2。作为一种更加严格的收敛判断，位移收敛条件也可以要求每一个节点自由度上的位移满足条件|L| e(8.109)有时，位移收敛条件和力收敛条件可以同时运用。 结构的非线性分析需要经过反复迭代 运算，需要的时间一般要比线性分析时多得多。 这样，收敛条件的正确与灵活运用就十分重 要，以便在获得满意的计算精度的同时节约计算时间，这里需要一定的实践经验作为基础。上述带有流动坐标的迭代法对于板和壳的大位移分析也可以适用，只是单元刚度矩阵的形式不同，此时