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

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

图8.1a 图8.1b

可见，工程实际中分枝屈曲现象实为罕见，它仅出现在完全无结构缺陷，完全沿轴向加压的绝对直杆及完整空球壳在均匀外压的情况下。分枝屈曲现象虽然罕见，但实际中有不少结构屈曲状态接近分枝屈曲，而分枝屈曲的计算工作量又远小于计算极限屈曲的工作量，况且，不少作者得出结论，一些中等非线性的屈曲状态，可以用线性屈曲问题特征向量的线性组合近似得到。因此线性屈曲理论还是有其实际价值。
屈曲的含义可简述为：结构处于一种平衡状态，载荷增量为一个微量，其位移增量很大。用方程来表达这种物理现象，则由总体拉格朗日列式法建立的结构刚度方程（8.21）变成为
（8.22）
根据式（8.18）和（8.20），有
（8.23）
将式（8.23）代入（8.22）得
（8.24）
在线性屈曲情况下，屈曲前结构处于原始位形的线性平衡状态，因此上式中的大位移矩阵应为零，此时式（8.24）简化为
（8.25）
由式（8.16）可以看出，并不明显地包含位移增量，在小变形情况下，该矩阵与应力水平成正比。由于屈曲前线性假设，多数情况下，应力与外载也为线性关系，因此，如果令某一参考载荷对应的初应力刚度阵为，令屈曲极值载荷为，与参考裁荷有关系，称载荷比例因子，为极值裁荷的比例因子。极值载荷时的初应力刚度阵为
（8.26）
代入式（8.25），则有
（8.27）
可以看出式（8.27）是一个广义特征值方程，也是经典弹性稳定理论的最后控制方程。
实际求解时可按以下步骤进行：
（1）按线弹性问题的有限元法形成各单元的刚度矩阵，并用常规方法组装成结构刚度矩阵，即。
（2）对结构施加参考载荷，并求解有限元方程，进而可求得各单元节点应力。
（3）对有膜应力存在的单元按式（8.16）形式初应力矩阵，并组装成结构初应力矩阵，即。
（4）求解广义特征值问题，即式（8.27），得到最低几阶特征值及对应的特征模态。理论上它们都是平衡模态，当载荷达到分枝载荷时，平衡由一种模态“跳列”另一种模态，这也就是分枝屈曲的含义。
（5）从分枝载荷中挑选一个最小的载荷，即，作为极值载荷或临界载荷。
值得指出的是，由式（8.27）所表征的线性屈曲问题，是建立在下述假设的基础上的，即假设线性应变刚度矩阵在屈曲前不产生明显的变化，且初应力矩阵简单地与应力水平成正比。如前所述，这在实际问题中是很罕见的。由式（8.27）所确定的极值载荷只能是近似的。由于线性屈曲理论存在精度差，以及适用范围窄的限制，所以，在一般情况下，应当用切线刚度矩阵来研究这个问题。当时，发生随遇平衡。显然，这里应该用逐次逼近的方法进行求解。关于非线性屈曲问题的有限元解法，读者可参考其它文献。
8.4 板的大挠度及线性屈曲
8.4.1 基本问题
在板的大挠度问题中，板中的内力除弯曲内力外，还有薄膜内力，如图8.2所示，在这种情况下，横向位移会引起薄膜应力，因此在板的大挠度问题中，面内变形和弯曲变形不再象小挠度板那样能分别处理，二者是相互耦合的。

图8.2 图8.3

和以前一样，平板的应变可用中面的位移来描述。如果令平面与中面一致，则广义应变和广义应力向量可以表示为
（8.28）
（8.29）
式中上标p和b分别表示面内和弯曲的分量。
由图8.3可以看出，挠度使中面在x和y方向产生附加伸长和附加角变形，因此广义应变可表达成
（8.30）
式中第一项是我们熟知的线性项，而第二项则是非线性项。
如果所考虑的材料是线弹性的，则平板的弹性矩阵由平面应力弹性矩阵和弯曲弹性矩阵组成，即
（8.31）
平板单元的位移可以采用适当的形函数，通过节点位移来表达，写成
（8.32）
为方便起见，将节点位移向量也分为面内的和弯曲的两部分，即

（8.33）
而

这样，形函数也应分成
（8.34）
通过上述规定，除了非线性应变项以外，其余各项都和标准的线性分析相同，不再重复。
8.4.2 应变矩阵的计算
为了确定单元的切线刚度矩阵，我们先来讨论单元的应变矩阵的计算。首先应注意到
（8.35）
式中

这里，和分别是平面单元和弯曲单元按线性分析所得到的标准矩阵，而是应变的非线性项引起的，它可以通过对参数取微分来导出。下面就来推导的表达式。
在式（8.30）中，非线性应变分量可以方便地写成
（8.36）
式中
（8.37）
(8.38)
w的一阶导数可以用弯曲节点位移表示成
（8.39）
其中
（8.40）
可见只与单元的节点坐标有关。
将式（8.36）取微分，有
（8.41）
注意到矩阵和具有两个有用的性质：
1．设是一个任意向量，则有
（8.42）
于是有
（8.43）
2．设是一个任意向量，则有
（8.44）
利用上述性质1，并注意到式（8.39），可将式（8.41）写成
（8.45）
观察上式，由应变矩阵的定义，可见
（8.46）
8.4.3 单元切线刚度矩阵的计算
根据式（8.18），单元的切线刚度矩阵由三部分组成，即
（8.47）
由第二章和第四章给出的刚度矩阵，可将线性小变形的刚度矩阵写成下式
（8.48）
将式（8.35）中的和代入式（8.15），可得大位移矩阵
（8.49）
最后，利用式（8.16）可求出初应力矩阵。为此，将式（8.35）中的进行微分，得

把上式代入式（8.16），并利用式（8.29）和（8.46），给出

利用式（8.44），有

最后可将初应力矩阵写成
（8.50）
式中
（8.51）
是平板对称形式的初应力矩阵。
将式（8.48）至（8.50）代入（8.47），便得到板的切线刚度矩阵，其余计算步骤已如8.2.2节所述，不再重复。
8.4.4 板的屈曲问题
板的线性屈曲问题，可从式（8.25）出发作为特征值问题考虑。对于只在平面内承受压力载荷的情况，先求问题的弹性解，然后用式（8.51）求出，而由线弹性方法求得，此时，在完成结构刚度矩阵的组装后得到了广义特征值问题
（8.52）
解之，可得平板线性屈曲问题的解。
上述平板的几何非线性有限元理论，也适用于由平板单元组成的壳体结构的几何非线性分析，只是多了一步将局部坐标系中的切线刚度矩阵及其它有关量转换到整体坐标系中的工作。
值得指出的是，通常对壳体结构进行线性屈曲的分析结果远远大于实际的失稳载荷。因此，确定变形对屈曲的影响是十分重要的，也就是说，为了得到正确的结果，必须进行完全的非线性分析。
8.5 三维单元的大应变和大位移公式
上节平板大挠度问题中，非线性应变一位移关系是在特定的情况下建立的。根据格林（Green）应变的定义，现在可以从一般方法出发推导出几何非线性理论的有限元计算公式，而不管应变或位移是大的还是小的。用直角坐标表示的格林应变为
（8.53）
其它四个应变只要通过轮换，即可得到。
在小位移的情况下，忽略二次项，就得到线性应变公式。这时，和是原来平行于坐标轴的线段的伸长率，而和表示原来平行于坐标轴的两垂直线段夹角的变化。但在应变较大时，上述几何意义就不再成立。
现在我们来推导三维应力状态下非线性的的一般计算公式。稍加改变即可得到一维和二维情况下的计算公式。对于板壳问题，可在特定条件下忽略某些项后就可导出上节中的结果。对杆系结构我们将在下节详细讨论。
8.5.1 应变矩阵的推导
变形体某点的应变应是线性应变和非线性应变之和，即
（8.54）
式中
（8.55）
按照式（8.53），非线性应变可用下式表示
（8.56）
式中

是一个6×9的矩阵。是一个9×1的列向量，而

余类推。
容易验证上述定义的正确性，并重新建立上节中矩阵和的两个性质。我们再次有
（8.57）
如果用形函数和节点位移向量来表示，则可写出
（8.58）
式中的表达式和上节中的式（8.40）类似，为形函数对坐标的一阶导数所组成，因而只与节点坐标有关。将上式代入式（8.57）得

因此由应变矩阵的定义即有
（8.59）
8.5.2 单元切线刚度矩阵的推导
注意到式（8.4），即，由式（8.14）积分得到，而由式（8.15）可求得，又由式（8.16），则有
（8.60）
由式（8.59）

注意到不含，由上式得
（8.61）
将式（8.61）代入（8.60）的右端，则
（8.62）
同上一节一样，利用和的第二个性质，得
（8.63）
式中

由式（8.58）可得
（8.64）
将式（8.64）代入（8.63）则
（8.65）
再将式（8.65）代入（8.62），使得到具有对称形式的初应力矩阵
（8.66）
至此，三维单元的切线刚度矩阵就可由下式求得

8.6 杆系结构的大位移分析
对杆系结构进行大位移分析时，可采用总体的拉格朗日列式法或更新的拉格朗日列式法。对总体的拉格朗日列式法，我们将推导出桁杆单元和梁单元的切线刚度矩阵及其显式，其余计算均如8.2.2节所述。而对于更新的拉格朗日列式法，则以刚架单元为例，重点介绍其计算原理和实施步骤。
8.6.1 桁架单元的切线刚度矩阵
考虑一个横截而积为A，弹性模量为E，长度为L的桁杆单元，它在发生变位前后的位置如图8.4所示。

图8.4

在小位移的情况下，桁杆单元上某一点的轴向应变为
（8.67）
如果节点的位移比较大，则由于横向位移v会使单元发生附加的轴向伸缩。这种附加的轴向应变与单元在变位过程中所转过的角度有关，此时单元的轴向应变可表示为
（8.68）
式中左端的第二项是考虑大位移时附加的非线性项。显然上式也可由格林应变公式直接写出。
在大位移分析中，桁杆单元的轴向和横向位移函数可精确地取x的一次函数，即取
（8.69）
同以前线性分析时的单元分析过程一样，由杆两端节点位移条件可解出，于是单元上任意点的位移可用杆端节点位移表示为
（8.70）
式中
（8.71）
为桁杆单元的形函数矩阵。将式（8.70）代入式（8.68）再考虑到式（8.71），有
（8.72）
于是

（8.73）
将上式与式（8.2）比较，并考虑到式（8.4），有
（8.74）
(8.75)
在式（8.75）中
（8.76）
为单元在变位过程中发生的转角，如图8.4所示，于是式（8.75）可以写成
（8.77）
将式（8.74）代入式（8.14），可得线性分析时桁杆单元的刚度矩阵，如式（5.6）所示。
将式（8.74）和（8.77）代入式（8.15）可得桁杆单元的大位移矩阵如下

（8.78）
由式（8.78）可以看出，当单元处于变位以前的原始位置时，，因而大位移矩阵。
下面来推导桁杆单元的初应力矩阵。由式（8.77）并考虑式（8.76）的关系，有
（8.79）
式（8.79）与点有坐标无关，而对于桁杆单元来说，是一个常数，因此在积分时这些项均可提到积分号之外。若记为桁杆单元的轴向力，则有

将上式与式（8.16）比较可得
（8.80）
这就是桁杆单元的初应力矩阵，它计及了桁杆单元的轴向力对杆端横向位移的影响。
最后将式（5.6），（8.78）和（8.80）代入式（8.18），便得到桁杆单元的切线刚度矩阵。另外应注意到，由于均为对称矩阵，所以桁杆的切线刚度矩阵也是对称矩阵，由所组装的结构刚度矩阵也是对称的，这对求解是很有利的。
8.6.2 刚架单元的切线刚度矩阵
在大位移的情况下，仍可假定刚架单元的位移函数与小位移情况下的相同，即单元的轴向位移u是局部坐标x的线函数，而横而位移v则是x的三次函数。用式子表示如下

（8.81）
由杆两端的节点位移条件，可从上式解得，于是单元上任意点的位移可用杆端位移表示成
（8.82）
式中
（8.83）
若记

（8.84）
于是可表示成

而式（8.82）可写成
（8.85）
假若单元的曲率仍能用横向位移v的二阶导数近似地表示，则代表弯曲所引起的轴向应变的项仍与线性分析时相同。而由于横向位移引起的单元轴向应仍应按式（8.68）计算，此时有
(8.86)
式中代表由线位移所引起的轴向应变，代表由弯曲所引起的轴向应变。而
（8.87）
为线性分析时的应变项。
（8.88）
为非线性项。
由式（8.85）和（8.87）得

将上式取微分，并注意到和不含，有
（8.89）
由式（8.85）和（8.88）可得

将上式取微分，有
（8.90）
式中

（8.91）
由式（8.86）可得

将式（8.89）和（8.90）代入上式，得

将上式与式（8.2）比较，并注意到式（8.4）的关系有
（8.92）
(8.93)
将式（8.84）代入（8.92），得

（8.94）
将代入式（8.14）就得到线性分析时的单元刚度矩阵，如式（5.37）所示。将式（8.93）
和（8.94）代入式（8.15），就可以得到单元的大位移矩阵
（8.95）
式中
（8.96）
而
（8.97）
分别为单元i端和j端的截面转角。
在计算式（8.95）中的积分时，可以应用下列积分公式
（8.98）
以及
（8.99）
以下我们来推导初应力矩阵。对刚架单元来说，轴向应力由两部分组成，可表示为
（8.100）
其中代表单元由于拉伸或压缩而引起的轴向应力，在单元上任意点都相同，可以将这一常数应力用表示。代表由弯曲变形而引起的轴向应力。由式（8.93）可知中第二行元素全部为零，这样，在式（8.16）中将不发生作用，这说明单元的初应力矩阵与它的弯曲应力无关。若记中的第一行为

注意到矩阵中的项只是x的函数，由式（8.16）得

将式（8.91）代入上式，引入式（8.97）的关系，同时注意到为单元的轴向力，则有

将式（8.98）中的积分公式代入上式便得到单元的初应力矩阵，即
（8.101）
最后将式（5.37），（8.95）和（8.101）代入式（8.18）,便得到了刚架单元的切线刚度矩阵。同桁杆单元一样，由于均为对称矩阵，所以刚架单元的切线刚度矩阵也是对称的。
8.6.3 更新的拉格朗日列式法
前面我们在推导桁杆单元和刚架单元的切线刚度矩阵时，采用了总体的拉格朗日列式法。其特点是，所有静力学和运动学方面的变量，例如单元刚度矩阵，单元节点位移和单元节点应力等，都是以t=0时刻的位形，即变形时的位形为参照系统。这里只是用来描述物体位形变化发展过程的一个标志量，而并不一定具有真正的时间意义。
上述各量，即单元的刚度矩阵，节点位移和节点应力等量的描述和计算，也可以用物体变形过程中某一时刻t的位形为参照系统，进而推算时刻时物体的位形。由于t时刻的位形和坐标值随计算而变化，所以称为更新的拉格朗日列式法。由此产生的非线性分析方法可称作带流动坐标（或变更坐标）的迭代法。这一方法对于杆系结构的大位移分析特别显示出其优越性。尤其是在杆件发生比较大的转动时，采用这一方法更为适宜。这是因为通过变更局部坐标系可以方便地描述单元的刚体转动，从而较容易地确定变位后的单元在变形后的结构中所发挥的作用。
我们知道结构坐标系中的单元杆端力是由结构坐标系中的单元刚度矩阵与该单元杆端位移向量的乘积得到的。而结构坐标系中的单元刚度矩阵不仅取决于单元本身的属性，还与该单元所处的方位有关。在线性的位移分析中，由于节点位移引起的单元方位的变化十分微小而可以忽略，此时在计算结构的位移和节点合力时仍然能利用变形前方位的单元刚度矩阵。如果结构的位移比较大，或者说在非线性问题中，单元方位的这种变化就不能忽略，此时在计算节点合力时就应该使用变形后方位的单元刚度矩阵。
在进行结构的大位移分析时，可以将按线性分析所得到的节点位移作为结构位移的第一次近似值。根据上述节点位移可以对单元刚度矩阵进行修改，从而反映单元在变位后的位置上所起的作用。在带有流动坐标的迭代法中所采用的做法是：根据已求的节点位移修改单元的坐标转换矩阵，从而达到修正结构坐标系中单元刚度矩阵的目的。根据修正后的各单元刚度矩阵乃至刚度方程，可以计算出节点合力。按照上述结构位移的第一项近似值算出的节点合力与节点所受到的外载荷并不相等，也就是说此时节点的平衡条件未被满足。这是因为按线性分析所得到的节点位移并不代表结构真正的平衡位置。于是，在这样的位置上结构当然也就无法保持平衡。现将原结构的等效节点载荷与上述节点合力之差称为节点的不平衡力。为了求得结构真正的平衡位置，可以将不平衡力作为一组新的外载荷施加于上述已发生变形的结构上，求得节点位移的修正值又可以重新修改各单元的坐标转换矩阵，并进而得到新的节点合力和节点不平衡力，继而再将新的节点不平衡力施加于变形以后的结构。重复上述过程一般总可以使节点不平衡力减小到可以被忽略的水平，此时的节点位移所对应的便是结构在发生大位移之后真正的平衡位置。按照上述平衡位置可以计算结构在大位移情况下的杆件内力。以上就是带有流动坐标的迭代法分析大位移问题的基本思路。
以下讨论流动坐标系中的杆端位移和杆端力的计算方法。图8.5（a）给出了在结构坐标系中的一个未变形的刚架单元，它在弯形后的位形如图8.5b所示。是该单元的流动坐标，它的坐标原点位于变形后的杆端，轴沿变形后杆端节点的连线方向。

(a) （b）
图8.5

对比图8.5a、b可得
（8.102）
于是，刚架单元在流动坐标系中的节点位移可表示为
（8.103）
因此，刚架单元在中的节点位移向量为
　　　　　　　　　　　　　　　　　　 (8.104)
此时在上述局部坐标系中的杆端力可表示为
（8.105）
其中为局部坐标系中的单元刚度矩阵，仍可按式（5.37）计算。
由于采用了流动坐标系，方向角就成为杆端位移的函数，可由式（8.102）计算。这样，单元的坐标转换矩阵和通过作用于得到的结构坐标系中的单元刚度矩阵也都成为杆端位移的函数。
综上所述，在采用带有流动坐标的迭代法时，一个典型的迭代过程包括以下步骤。
（1）为外载荷作用下的结构假定一组节点位移。
（2）根据结构坐标下的节点位移向量，确定单元两端的位置，建立单元的局部坐标系。
（3）计算上述局部坐标系中的杆端位移向量。
（4）形成局部坐标系中的单元刚度矩阵和杆端力向量。
（5）将和转换到结构坐标系，得到和。
（6）对所有单元重复（2）至（5）的步骤。生成结构刚度矩阵和节点合力。
（7）计算节点不平衡力。
（8）求解结构方程得节点位移增量。
（9）将叠加到节点位移向量中。
（10）收敛条件判断，如果不满足则返回到步骤（2）。
以上的过程可概括为如下的迭代公式
（8.106）
一般来说作第一次迭代运算时节点位移可取为线性分析所得的节点位移值，这样通常有利于加快收敛速度。但当结构的非线性程度很高时如何假定却是一个值得研究的问题。
例如对于图8.6所示结构，按线性分析方法无法确定结构的位移或者说结构的位移将成为无穷大，此时就需要为载荷作用点C先假定一个适当大小的竖向位移，然后才能开始迭代过程。而在某些情况下需要把外载荷分成数级逐步施加到结构上。

图8.6

迭代过程的收敛可根据上述不平衡力进行判断，当不平衡力和外载荷的比率减小到一个给定的限度时则可以认为迭代过程已达到了收敛，这样的收敛条件称为力收敛条件。在结构的大位移分析中一般采用位移收敛率条件效果更好一些，因此这种收敛条件的应用更为普遍。位移收敛条件可由不同的形式提出，若记

（8.107）
式中是经过n次迭代运算得到的节点位移向量，则位移收敛条件的一种形式是
（8.108）
这里er是精度要求，可以根据工程的要求和问题的性质而定，一般可取10-6到10-2。作为一种更加严格的收敛判断，位移收敛条件也可以要求每一个节点自由度上的位移满足条件
（8.109）
有时，位移收敛条件和力收敛条件可以同时运用。结构的非线性分析需要经过反复迭代运算，需要的时间一般要比线性分析时多得多。这样，收敛条件的正确与灵活运用就十分重要，以便在获得满意的计算精度的同时节约计算时间，这里需要一定的实践经验作为基础。
上述带有流动坐标的迭代法对于板和壳的大位移分析也可以适用，只是单元刚度矩阵的形式不同，此时与杆端力和杆端位移相对应的则是一般单元的节点力和节点位移。采用流动坐标对结构进行大位移分析时，所有本质性的非线性特性均是通过坐标转换计及的。只要单元划分得足够小，带有流动坐标的迭代法可以为结构大位移问题提供满意的解答。最后，需要指出的是，在采用带有流动坐标的迭代法时，也有将单元的被应力矩阵也考虑在内的做法，特别是在杆件的轴向力或者是板、壳单元的中面应力相当大的时候。
8.6.4 算例
采用带有流动坐标的迭代法对图8.7所示结构进行大位移分析。设在图示载荷作用下，杆件材料仍将处在弹性工作阶段，弹性模量，杆件的模截面面积A=4cm2。

图8.7
该结构因为C点处的铰很靠近AB连线，结构在图示裁荷作用下将表现出明显的非线性特性。根据对称性，该结构只有一个线位移自由度，即C点的竖向位移，可记为。将CB杆看作一个单元，其长度为，将该单元局部坐标系的原点设在C点，单元的方向角。
结构的刚度矩阵显然应等于CB杆单元刚度矩阵的两倍，即

其中为流动的局部坐标系的方向角。采用带有流动坐标的迭代法分析大位移问题时，将跟随结构的变位而发生变化，即有

由式（8.103）可知单元的伸长量

此时局部坐标系中的杆端力和结构的节点合力F分别为

设采用位移收敛条件

带有流动坐标迭代法的求解过程如表8.1所示。

迭代次数

F（kN）
R—F
（kN）

0
0
0
20
-31.269
-31.269
1
5.369
-4279.284
-2711.470
2691.470
16.882
-14.387
2
1.314
-1050.809
-339.859
319.859
7.646
-6.741
3
0.361
-288.981
-50.328
30.328
2.500
-4.241
4
0.175
-139.621
-17.182
-2.818
-0.420
-4.661
5
0.202
-161.247
-21.434
1.434
0.203
-4.458
6
0.188
-150.619
-19.414
-0.586
-0