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

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

8.1引言

前面各章所讨论的问题都是在小变形假设的前提下进行的，即假定物体所发生的位移远小r物体自

身的几何尺寸，应变远小于1。在此前提下，建立物体或微元体的平衡条件时可以不考虑物体的位置和

形状（简称位形）的变化，因此在分析中不必区别变形前后位形的差异，且应变可用一阶无穷小的线性

应变表达。

实际上，上述假设有时是不成立的。即使实际应变可能是小的，且不超过材料的弹性极限，但如果

需要精确地确定位移，就必须考虑几何非线性，即平衡方程应该相对于变形后的位置得出，而几何关系

应该计及二次项。例如平板大挠度理论中，由于考虑了中面内的薄膜应力，求得的挠度比小挠度理论的

结果有显著的减低。再如在结构稳定性问题中，当载荷到达一定数值后，挠度比线性解答予示的结果更

剧烈地增加，并且确实存在承载能力随继续变形而减低的现象。在冷却塔、薄壁结构及其它比拟细长的

结构中，几何非线性分析都显得十分重要。

几何非线性问题可以分为以下几种类型：

（1）大位移小应变问题。一般工程结构所遇到的几何非线性问题大多属于这一类。例如高层建筑

或高耸构筑物以及大跨度网壳等结构的分析常需要考虑到结构大位移的影响。

（2）大位移大应变问题，如金属压力加工中所遇到的问题就属于这一类型。

（3）结构的变形引起外载荷大小、方向或边界支承条件的变化等。

结构的平衡实际上是在结构发生变形之后到达的，对于几何非线性问题来说，平衡方程必须建立在结构

变形之后的状态上。为了描述结构的变形需要设置一定的参考系统。一种做法是让单元的局部坐标系始

终固定在结构发生变形之前的位置，以结构变形前的原始位形作为根本的参考位形，这种分析方法称作

总体的拉格朗日（Lagrange）列式法；另一种做法是让单元的局部坐标系跟随结构一起发生变位，分析

过程中参考位形是不断被更新的，这种分析方法称作更新的拉格朗日列式法。

本章首先对几何非线性问题作一般性讨论，从中导出经典的线性屈曲问题的公式；然后建立平板大

挠度问题和壳体的大位移（及大转动）分析的有限方法公式；接着还给出了大应变及大位移的一般公式,

最后还详细讨论了杆系结构几何非线性问题的有关公式。在讨论中我们采用总体的拉格朗日列式法，但

对杆系结构，为应用方便我们绐出了两种列式法的公式。

一般性讨论

8.2.1理论根底

无论是对于何种几何非线性问题，虚功原理总是成立的。由虚功原理，单元的虚功方程可以写成如

下的形式

n卜｝"｛"小-w｝"但丫=。（）

V

其中｛尸｝为单元节点力向量，为单元的虚应变，加｝为节点虚位移向量。

增量形式的应变一位移关系可表示为

=［即｛邸（）

上式中”伤『表示单元节点位移份｝'的微分。根据变分与微分运算在形式上的相似性，有

m服｝'（8.3）

以上两式中限］称为大位移情况下的增量应变矩阵，代表了单元应变增量与节点位移增量之间的关系。

在大位移情况下®］应是节点位移的函数。

假设将上述应变增量矩阵分解为与节点位移无关的局部［B。］和与节点位移有关的局部仍乂伤}）］两

局部组成，即

同=同+回［0

此时［纥］也就是一般线性分析时的应变矩阵。

将式0代入（），并考虑到节点虚位移标］的任意性，可将单元的平衡方程写成

%同（{""一{尸}'=00

V

按照式（）可以对整个结构建立有限元列式，这种列式方法可称为全量列式方式，在几何非线性分析中，

按照这种列式方法得到的单元和结构刚度矩阵•般是非对称的，于求解不利。因此，在分析非线性问题

时大多采用增量列式法。以下就着重介绍这一方法。

式（）所示的平衡方程可以写成微分的形式

JJjd（同—）办々{*'=（）0

V

由于在儿何非线性问题中，应变矩阵限］和应力匕『都是节点位移的函数，因此有

”（［研{a}e）=”［4才+回d{a}e（）

将式0代入0,那么有

JJjd\B］{o-\dv十.同d{cr}ldv=d［F\（）

VV

单元内部的应力增量与应变增量存在确定的关系，这种关系可以用增量形式表示为

d{cr}e=［Dp{£}r0

式中称为应力-应变关系矩阵，或称为材料的本构关系矩阵。如果材料属于线性弹性的，［。］将是

一个常数矩阵。并且，对于线性弹性材料来说有

■=同（卜『-{q}'）+Me（）

上式中{%}'和{5}分别为单元材料中可能存在的初应变和初应力。

将式0代入式0就可以得到应力增量与单元节点位移增量之间的关系

”{叶=［“而{"（）

招式0代入式0后得

47『=网间+［4］）丽，0

于是，式（）左端中的第二项便可表示为

JJJ回人=皿仍—。加+皿M［DIBL］dv

VVV（）

。r

+ffi/fW］dv+JJJ［BL］［D\BL］dv））d［^

VV

假设记

k0］=册与H。］限如（8.14）

V

k0］是与单元节点位移无美，它就是一般线性分析时的单元刚度矩阵。式（）右端第二层括号内的项可

记为

h］=JJJ（回『［。加.］+,j向综］+®］）小（）

V

》』称为单元的初位移矩阵或大位移矩阵，表示单元位置的变动对单元刚度矩阵的影响。

现在再来看式0左端的第一项。考虑到式0的关系并注意到［B。］与节点位移无关，因此对•节点

位移的微分等于零，对于一个确定的有限元分析模型，式0左端的第一项可一般地写成

洞'{垃小=e瓦卜仍（）

JJJJJJd［BLY{a}dv=

VV

上式中卜」称为单元的初应力矩阵或几何刚度矩阵，它表示单元中存在的应力对单元刚度矩阵的影响。

由上式0和式（），并考虑到式（）和。的关系，有

（庶］+应］+瓦］）4/,=4丹°o

假设记

k/=k（）］+k/+kJ（）

就称为单元的切线刚度矩阵。此时，有增量形式的单元刚度方程

ee

［kTyl{6}=d{F}0

由此可以看出，单元切线刚度矩阵心」代表了单元于某种变形位置时的瞬时刚度，或者说代表了单元节

点力与节点位移之间的瞬时关系。

有了单元切线刚度矩阵就可以按照常规的方法，即单元集成法组装结构的切线刚度矩阵，即有

［勺］=£k］［）

并进而得到结构的增量刚度方程

及卜⑻二小}0

前面在推导式0时，假定载荷{F『与变形无关。但有些情况并非如此。例如，作用于特大变形结

构I.的压力载荷，与变形有关的气动载荷便是这样。在这种情况下，式（）应计及载荷相对于〃砂}的

微分项，本书后面的推导中均不考虑这一影响。

8.2.2求解方法

对于实际应用，载荷增量不可能取成微分的形式，总是一个有限值。于是，按式()求得的位移增

量使结构偏移了其真实的平衡位置。为了解决这一问题，可以根据当时的结构位移情况按式()求各单

元上作用的节点力，并继而求得各节点合力。然后将外载荷与上述节点合力之差，即节点的不平衡力，

作为一种载荷施加于结构，由此求得节点位移的修正值。上述过程也可以反复屡次。

综上所述，总体的拉格朗日列式方法的一次完整的迭代步骤一般可归纳如下：

(1)按线性分析得到节点位移的初值份｝。

(2)形成局部坐标系中的单元切线刚度矩阵［加］,并按式0计算单元的节点力｛尸｝"

(3)将吟］和｛F｝"转换到整体坐标系。

(4)对所有单元重复(2)至(3)的步骤。生成结构的切线刚度矩阵［KJ和节点力合力仍｝。

(5)计算节点不平衡力｛〃｝产｛八—Z｛处"

(6)求解结构刚度方程［勺］△⑻।=｛%,得节点位移增量△伤｝1

(7)将一份,叠加到节点位移向量伤｝中,即伤｝2=同+丽］。

(8)收敛条件判断，如果不满足那么反回到步骤⑵。

上述在总载荷下进行迭代的方法有时会遇到困难。在非线性程度较高的问题中可能收敛较慢，此外,

当解答非唯一时，有可能得到实际上不需要的那个解。在这种情况下，可采用节中所介绍的增量法求解,

并得到每•增量步的非线性解.。如迭代中再带有自平衡校正，并采用小的载荷增量，通常一步运算就能

足够精确地得到该增量步的解。

以上两节所介绍的增量形式的总体拉格朗日列式法，在结构的非线性分析中应用十分广泛，有关计

算公式及求解方法对板、壳或杆件体系的非线性分析都同样适用。由上面的分析也可以看出，米用总体

的拉格朗日方法求解非线性问题的关键是形成单元的切线刚度矩阵。

8.3屈曲问题

非线性分析，尤其是几何非线性分析在很多情况下是估算一个结构在失去稳定性前所能承受的最大

教荷。这是结构屈曲问题的研究目标，是固体力学的一个重要分支，也是工程实践中经常出现的问题。

小位移线性理论假设在结构受载变形过程中忽略了结构的位移变化，因此在加载的各个阶段总是认

为结构在未加载的原始位形上产生平衡，当屈曲发生时，结构位形突然跳到另一个平衡位置，图8.1(a)

为线性屈曲的示意图。2为裁有比例因子，其含义梢后会讲到，它与位移S在屈曲前为线性关系，当载

荷到达极限值［图中分枝点)时结构失积，曲线改变，结构平衡转向另一模态。这就是线性屈曲

也称分枝屈曲。

严格说来，结构的平衡实际上是在结构发生变形之后到达的，因此，从加载一开始就出现了几何非

线性的特性，图8.1(b)为非线性屈曲的示意图，当载荷比例因子增加时，曲线是非线性的，一直

到达极限，这种在结构发生变形•直到失稳，在变形后的位形上考虑平衡•直到达极限的方法称非线性

屈曲或极限屈曲。

图图

可见，工程实际中分枝屈曲现象实为罕见，它仅出现在完全无结构缺陷，完全沿轴向加压的绝对直

杆及完整空球壳在均匀外压的情况下。分枝屈曲现象虽然罕见，但实际中有不少结构屈曲状态接近分枝

屈曲，而分枝屈曲的计算工作量又远小于计算极限屈曲的工作量，况且，不少作者得出结论，i些中等

非线性的屈曲状态，可以用线性屈曲问题特征向量的线性组合近似得到。因此线性屈曲理论还是有其实

际价值。

屈曲的含义可简述为：结构处于一种平衡状态，载荷增量为一个微量，其位移增量很大.用方程来

表达这种物理现象，那么由总体拉格朗日列式法建立的结构刚度方程（）变成为

应卜{5}=00

根据式（）和（），有

及］=区］-及］+区10

将式0代入0得

（及］+及［+［K』）d{b}=00

在线性屈曲情况下，屈曲前结构处于原始位形的线性平衡状态，因此上式中的大位移矩阵［KJ应为

零，此时式（）简化为

（及］+及皿&}=0（）

由式（）可以看出，［K/并不明显地包含位移增量d{b},在小变形情况下，该矩阵与应力水平成

正比。由于屈曲前线性假设，多数情况下，应力与外载也为线性关系，因此，如果令某一参考载荷F'.对

应的初应力刚度阵为［K；］,令屈曲极值载荷为产「，与参考裁荷有/，=元/'关系，4称载荷比例因子,

才为极值裁荷的比例因子。极值载荷时的初应力刚度阵为

依］=尤应］⑹

代入式〔），那么有

（及］+牙［《］以伤}=01）

可以看出式0是一个广义特征值方程，也是经典弹性稳定理论的最后控制方程。

实际求解时可按以下步骤进行：

（1）按线弹性问题的有限元法形成各单元的刚度矩阵［尢』，并用常规方法组装成结构刖度矩阵,

即K]=Zko]0

(2)对结构施加参考载荷厂｝,并求解有限元方程火监｝二｛尸｝,进而可求得各单元节点应力。

(3)对有膜应力存在的单元按式()形式初应力矩阵并组装成结构初应力矩阵，即

(4)求解广义特征值问题，即式()，得到最低几阶特征值片及对应的特征模态。理论上它们都是

平衡模态，当载荷到达分枝载荷看尸「时，平衡由一•种模态“跳列”另•种模态，这也就是分枝屈曲的

含义。

(5)从分枝载荷中挑选一个最小的载荷，即耳京=4nin9■，作为极值载荷或临界载荷,

值得指出的是，由式()所表征的线性屈曲问题，是建立在下述假设的根底上的，即假设线性应变

刚度矩阵在屈曲前不产生明显的变化，且初应力矩阵简单地与应力水平成正比。如前所述，这在实际问

题中是很罕见的。由式。所确定的极值载荷只能是近似的。由于线性屈曲理论存在精度差，以及适用

范围窄的限制，所以，在一般情况下，应当用切线刚度矩阵[K/来研究这个问题。当[给卜⑻三。时,

发生随遇平衡。显然，这里应该用逐次逼近的方法进行求解。关于非线性屈曲问题的有限元解法，读者

可参考其它文献。

8.4板的大挠度及线性屈曲

8.4.1根本问题

在板的大挠度问题中，板中的内力除弯曲内力外，还有薄膜内力，如下图，在这种情况下，横向位

移会引起薄膜应力，因此在板的大挠度问题中，面内变形和弯曲变形不再象小挠度板那样能分别处理，

二者是相互耦合的。

和以前一样，平板的应变可用中面的位移来描述。如果令工-)，平面与中面一致，那么广义应变和

广义应力向量可以表示为

cPd2wd2wd2w]..

{2}=/>8..E.--------------------------Z--------()

odx2dy2dxd)\

匕｝=巴=牝N,N"%MMJ-O

<yv

式中上标p和b分别表示面内和弯曲的分量。

由图可以看出，挠度棒使中面在X和y方向产生附加伸长和附加角变形，因此广义应变可表达成

式中第一项为哪一项我们熟知的线性项，而第二项那么是非线性项。

如果所考虑的材料是线弹性的，那么平板的弹性矩阵曰平面应力弹性矩阵和弯曲弹性矩阵

［则组成，即

网=，》］(8.31)

平板单元的位移可以采用适当的形函数，通过节点位移来表达，写成

u

,…［N网’0

为方便起见，将节点位移向量也分为面内的和弯曲的两局部，即

区｝」倒n

阴｝J

而

k)=k次极｝小啖管)［

%oy

这样，形函数也应分成

一♦.°

通过上述规定，除了非线性应变项b”以外，其余各项都和标准的线性分析相同，不再重复。

8.4.2应变矩阵向的计算

为了确定单元的切线刚度矩阵［的］，我们先来讨论单元的应变矩阵［可的计算。首先应注意到

［司=闻+闯0

式中

出］』［第］。］屈］［0［明

。|_0［B种J［o0

这里，和忸彳］分别是平面单元和弯曲单元按线性分析所得到的标准矩阵，而是应变的非线性

项引起的，它可以通过,八对参数｛*'｝取微分来导出。下面就来推导口:］的表达式。

在式。中，非线性应变分量可以方便地写成

_

_。

，

一

_&如

_

生

11更1⑷

P>o—K

u

<=，>=

L-'、<-*

丝

J2《2

如⑪

av-v<I

ay&

.

中

式

迦

。

私

电

。

力

如

生

小

。

]=

f、

加

j一

小

[刎

丁

W的一阶导数可以用弯曲节点位移物叩表示成

r

加

，-

改

=<—

如

a

kI

其中

0

可见［G］只与单元的节点坐标有关。

将式0取微分，有

扉a｝=3"⑷网｛处

注意到矩阵［A］和｛0｝具有两个有用的性质：

1.设卜｝二卜,X』'是一个任意向量，那么有

条

dzo

(

\

a再

。dz

(

x¥

-M

沙

d/dz

f

i(

x¥)x

于是有

"同例=同的｝0

2.设｛），｝=%,刈,以「是一个任意向量，那么有

'、

包

而3

<、

dzxdz

o(

一

I-(一

x¥)x

静

必

y第"

_•

aor3=I\I

〃A/)，=4，yl

I1

——Fx1)

X为d

z\dzz

(ayl((%」

一

一

x-zxJ%x¥)一

通<

XI

利用上述性质1,并注意到式。，可将式（）写成

南―砥［4同/的0

观察上式，由应变矩阵的定义，可见

［昨同G］0

843单元切线刚度矩阵卜」的计算

根据式0,单元的切线刚度矩阵由三局部组成，即

k］=［%］+k］+瓦］0

由第二章和第四章给出的刚度矩阵，可将线性小变形的刚度矩阵写成下式

0

将式0中的同］和仍/代入式0,可得大位移矩阵

瓦］=1［第y-］疝回I用产0

最后，利用式（）可求出初应力矩阵［勺］,为此，将式（）中的也』进行微分，得

、

00

()

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

N、

*

0

忆》固=聊时水了°UMydv

（）J以

%

心

利用式（），有

N,

小川『N卜.MN“d{0}=«N”

N.N,

MyN”N.

最后可将初应力矩阵写成

00

()

0k

式中

MM,

网=川同N；$[G]c/v0

是平板对称形式的初应力矩阵。

将式（）至。代入1）,便得到板的切线刚度矩阵，其余计算步骤已如节所述，不再重复。

8.4.4板的屈曲问题

板的线性屈曲问题，可从式0出发作为特征值问题考虑。对于只在平面内承受压力载荷的情况，

先求问题的弹性解，然后用式0求出卜外，而由线弹性方法求得[%],此时kJ=0,在完成结构刚

度矩阵的组装后得到了广义特征值问题

（福+元阳4⑻=（）0

解之，可得平板线性屈曲问题的解。

上述平板的几何非线性有限元理论，也适用于由平板单元组成的壳体结构的几何非线性分析，只是

多了一步将局部坐标系中的切线刚度矩阵及其它有关量转换到整体坐标系中的工作。

值得指出的是，通常对壳体结构进行线性屈曲的分析结果远远大于实际的失稳载荷。因此，确定变

形对屈曲的影响是十分重要的，也就是说，为了得到正确的结果，必须进行完全的非线性分析。

8.5三维单元的大应变和大位移公式

上节平板大挠度问题中，非线性应变一位移关系是在特定的情况下建立的。根据格林（Green）应

变的定义，现在可以从一般方法出发推导出几何非线性理论的有限元计算公式，而不管应变或位移是大

的还是小的。用直角坐标表示的格林应变为

du1初、)QcR叭2

一+—(—)■+(―)■+(—广

dx2dxdxdx

0

dudvdududvdvdwdw

九尸区+说+------------1--------------1-------------

dxdydxdydxdy

其它四个应变只要通过轮换xryfz—-u―他一＞〃，即可得到。

在小位移的情况下，忽略二次项，就得到线性应变公式。这时，j,£,和田是原来平行于坐标轴的

线段的伸长率，而八y,八z和7h表示原来平行于坐标轴的两垂直线段夹角的变化。但在应变较大时，上

述几何意义就不再成立。

现在我们来推导三维应力状态下非线性的®Wk"的一般计算公式。稍加改变即可得到一维和二维

情况下的计算公式。对于板壳问题，可在特定条件下忽略某些项后就可导出上节中的结果。对杆系结构

我们将在下节详细讨论。

8.5.1应变矩阵区］的推导

变形体某点的应变应是线性应变也＞}和非线性应变{〃}之和，即

{£}={4}+{4}（）

式中

I包

,&

,史

,

,力

,包

_,&

1aa

v+vv()

-¥

a-z

aw

1-

-4-az

cv?

-

av

、

按照式0,非线性应变可用下式表示

MF00

0｛以0限｝:

00没｝7・Q

＞=扣场}()

0｛"｛疗帆｝)

阳,0初}『

0

式中

也｝，o0

0M0

但}]

。。没｝『{0}={My}

0泣｝,机｝,

泣｝70依｝，

跟『向F0

［A］是一个6X9的矩阵。｛夕｝是一个9X1的列向量，而

0y=包包如

［&dxdx_

余类推。

容易验证上述定义的正确性，并重新建立上节中矩阵［A］和｛/的两个性质。我们再次有

；皿网+：［A］cl｛0｝=［A］d｛0｝（）

如果⑹用形函数［N］和节点位移向量&｝来表示，那么可写出

别=冏伤y0

式中［G］的表达式和上节中的式（）类似，为形函数对坐标的一阶导数所组成，因而只与节点坐标有关。

将上式代入式0得

，/｛4｝=同⑷⑻

因此由应变矩阵的定义即将

间=［八同0

8.5.2单元切线刚度矩阵卜」的推导

注意到式（），即郎闻+闻,由式（）积分得到［%］,而由式（）可求得心』，又由式（）,

那么有

kk⑻"=JJJ胭『｛（j｝edv（）

V

由式0

也『=皿川

注意到［G］不含同"由上式得

d[Bj=[G]rd[A]r0

将式（）代入（）的右端，那么

瓦皿r="[G冏硝才办0

同上一节一样，利用同和网的第二个性质，得

对

44.｛疗=%/ayI称49｝=[“卜｛。｝0

式中

由式（）可得

d\0\=\G\l\6\0

将式（）代入（）那么

小4nb『=[M][G卜附0

再将式（）代入（），使得到具有对称形式的初应力矩阵

比卜啊MG出（）

V

至此，三维单元的切线刚度矩阵卜"就可由下式求得

[即]=—一]+一]

8.6杆系结构的大位移分析

对杆系结构进行大位移分析时，可采用总体的拉格朗日列式法或更新的拉格朗日列式法,对总体的

拉格朗口列式法，我们将推导出桁杆单元和梁单元的切线刚度矩阵及其显式，其余计算均如节所述。而

对于更新的拉格朗日列式法，那么以刚架单元为例，重点介绍其计算原理和实施步骤。

8.6.1桁架单元的切线刚度矩阵

考虑一个横截而积为A,弹性模量为E,长度为L的桁杆单元，它在发生变位前后的位置如下图。

在小位移的情况下，桁杆单元上某一点的轴向应变为

du

=——(8.67)

xdx

如果节点的位移比拟大，那么由于横向位移v会使单元发生附加的轴向伸缩。这种附加的轴向应变

与单元在变位过程中所转过的角度6有关.此时单元的轴向应变可表示为

式中左端的第二项是考虑大位移时附加的非线性项。显然上式也可由格林应变公式直接写出,

在大位移分析中，桁杆单元的轴向和横向位移函数可精确地取x的一次函数，即取

〃=«+ax

2►

U=+a4X

同以前线性分析时的单元分析过程一样，由杆两端节点位移条件可解出6于是单元上任意

点的位移可用杆端节点位移表示为

{/}={:=网伤}，（）

式中

[1--0-0

NLLO

L[]」=\XX

01--0-

LL\

为桁杆单元的形函数矩阵。将式0代入式0再考虑到式（），有

0

于是

0"{邸。1。B的