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

资源描述

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

《第八章几何非线性问题的有限元法.docx》由会员分享，可在线阅读，更多相关《第八章几何非线性问题的有限元法.docx（41页珍藏版）》请在冰豆网上搜索。

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

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

引言

前而各章所讨论的问题都是在小变形假设的前提下进行的，即假泄物体所发生的位移远小于物体自身的几何尺寸，应变远小于2。

在此前提下，建立物体或微元体的平衡条件时可以不考虑物体的位這和形状（简称位形）的变化，因此在分析中不必区别变形前后位形的差别，且应变可用一阶无穷小的线性应变表达。

实际上，上述假设有时是不成立的。

即使实际应变可能是小的，且不超过材料的弹性极限，但如果需要精确地确怎位移，就必须考虑几何非线性，即平衡方程应该相对于变形后的位置得出，而几何关系应该计及二次项。

例如平板大挠度理论中，由于考虑了中面内的薄膜应力，求得的挠度比小挠度理论的结果有显著的减低。

再如在结构稳左性问题中，当载荷达到一左数值后，挠度比线性解答予示的结果更剧烈地增加，并且确实存在承载能力随继续变形而减低的现象。

在冷却塔、薄壁结构及其它比较细长的结构中，几何非线性分析都显得十分重要。

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

（1）大位移小应变问题。

一般工程结构所遇到的几何非线性问题大多属于这一类。

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

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

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

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

为了描述结构的变形需要设置一泄的参考系统。

一种做法是让单元的局部坐标系始终固定在结构发生变形之前的位宜，以结构变形前的原始位形作为基本的参考位形，这种分析方法称作总体的拉格朗日（Lagrange）列式法；另一种做法是让单元的局部坐标系跟随结构一起发生变位，分析过程中参考位形是不断被更新的，这种分析方法称作更新的拉格朗日列式法。

本章首先对几何非线性问题作一般性讨论，从中导岀经典的线性屈曲问题的公式：

然后建立平板大挠度问题和壳体的大位移（及大转动）分析的有限方法公式：

接着还给岀了大应变及大位移的一般公式，最后还详细讨论了杆系结构几何非线性问题的有关公式。

在讨论中我们采用总体的拉格朗日列式法，但对杆系结构，为应用方便我们给出了两种列式法的公式。

一般性讨论

理论基础

无论是对于何种几何非线性问题，虚功原理总是成立的。

由虚功原理，单元的虚功方程可以写成如下的形式

JJJ&F｛b｝S，-⑹｝"｛f｝「=0（）

其中｛f｝为单元节点力向量，｛$¥为单元的虚应变，0扌为节点虚位移向量。

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

上式中dpy表示单元节点位移的微分。

根据变分与微分运算在形式上的相似性，有以上两式中国］称为大位移情况下的增量应变矩阵，代表了单元应变增量与节点位移增量之间的关系。

在大位移情况下恨］应是肖点位移的函数。

若将上述应变增量矩阵分解为与节点位移无关的部分［禺］和与节点位移有关的部分

国（0}）］两部分组成，即

広卜阴+R］（）

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

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

M［可{刃加-耐=0（）

按照式（）可以对整个结构建立有限元列式，这种列式方法可称为全量列式方式，在几何非线性分析中，按照这种列式方法得到的单元和结构刚度矩阵一般是非对称的，于求解不利。

因此，在分析非线性问题时大多采用增量列式法。

以下就着重介绍这一方法。

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

出〃（［可{巧川-d{F}"=0（）

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

d（［可何）=丽］{巧+［可d{b}「（）

将式（）代入（），则有

JJJ〃［秋dv+出［可d{a}edv=d{F}e（）

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

〃何=［刃忖（）

式中［D］称为应力一应变关系矩阵，或称为材料的本构关系矩阵。

如果材料属于线性弹性的,

［D］将是一个常数矩阵。

并且，对于线性弹性材料来说有

W=bk®一M

（）

上式中｛勺｝‘和｛b（）y分別为单元材料中可能存在的初应变和初应力。

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

〃何=［切［科何（）

将式（）代入式（）后得

〃何=［拆国）］+如）〃紂（）

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

ffj阳何dv=（fJJ［Bor［Z）K>+（J『［Bj［D血>

+nmqi垃mm［b』mem））〃耐（）

若记

kJ是与单元节点位移无关，它就是一般线性分析时的单元刚度矩阵。

式（）右端第二层括号内的项可记为

孔卜MM［dIbl］+［8』［DlB.］+［Bj［DlB^dv（）

［kJ称为单元的初位移矩阵或大位移矩阵，表示单元位苣的变动对单元刚度矩阵的影响。

现在再来看式（）左端的第一项。

考虑到式（）的关系并注意到［禺］与节点位移无关,因此对节点位移的微分等于零，对于一个确泄的有限元分析模型，式（）左端的第一项可一般地写成

JJJ〃［列｛/dv訂JJ创dv=\kc紂（）

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

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

（［灯］+阳+［俎］）〃｛疔=〃耐（）

若记

［心］就称为单元的切线刚度矩阵。

此时，有增量形式的单元刚度方程

旧0｛疔=〃｛町（）

由此可以看出，单元切线刚度矩阵［心］代表了单元于某种变形位宜时的瞬时刚度，或者说代表了单元节点力与节点位移之间的瞬时关系。

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

K小工旧］（）

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

［KT\l\d｝=d｛F｝（）

前而在推导式（）时，假立载荷｛F『与变形无关。

但有些情况并非如此。

例如，作用于特大变形结构上的压力载荷，与变形有关的气动载荷便是这样。

在这种情况下，式（）应计及载荷相对于〃0｝的微分项，本书后面的推导中均不考虑这一影响。

求解方法

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

于是，按式（）求得的位移增量使结构偏移了其真实的平衡位置。

为了解决这一问题，可以根据当时的结构位移情况按式（）求各单元上作用的巧点力，并继而求得各节点合力。

然后将外载荷与上述节点合力之差，即节点的不平衡力，作为一种载荷施加于结构，由此求得卩点位移的修正值。

上述过程也可以反复多次。

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

（1）按线性分析得到节点位移的初值

（2）形成局部坐标系中的单元切线刚度矩阵［知］,并按式（）计算单元的节点力｛F｝「。

（3）将［心］和｛F｝'转换到整体坐标系。

（4）对所有单元重复

（2）至（3）的步骤。

生成结构的切线刚度矩阵［K?

］和节点力合力｛必。

（5）计•算节点不平衡力初｝严｛必。

（6）求解结构刚度方程得节点位移增量△｛砒。

（7）将△｛》£叠加到节点位移向疑｛础中，即｛J｝2=｛J｝,+A｛J｝,.

（8）收敛条件判断，如果不满足则反回到步骤

（2）。

上述在总载荷下进行迭代的方法有时会遇到困难。

在非线性程度较髙的问题中可能收敛较慢，此外，当解答非唯一时，有可能得到实际上不需要的那个解。

在这种情况下，可采用节中所介绍的增量法求解，并得到每一增虽:

步的非线性解。

如迭代中再带有自平衡校正，并采用小的载荷增量，通常一步运算就能足够精确地得到该增量步的解。

以上两节所介绍的增疑形式的总体拉格朗日列式法,在结构的非线性分析中应用十分广泛，有关计算公式及求解方法对板、壳或杆件体系的非线性分析都同样适用。

由上而的分析也可以看岀，采用总体的拉格朗日方法求解非线性问题的关键是形成单元的切线刚度矩阵。

屈曲问题

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

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

小位移线性理论假设在结构受载变形过程中忽略了结构的位移变化，因此在加载的各个阶段总是认为结构在未加载的原始位形上产生平衡，当屈曲发生时，结构位形突然跳到另一个平衡位巻。

图⑻为线性屈曲的示意图。

兄为裁荷比例因子，其含义稍后会讲到，它与位移5在屈曲前为线性关系，当载荷达到极限值（图中分枝点）时结构失稳，曲线改变,结构平衡转向另一模态。

这就是线性屈曲也称分枝屈曲。

严格说来，结构的平衡实际上是在结构发生变形之后达到的，因此，从加载一开始就岀现了几何非线性的特性，图（b）为非线性屈曲的示意图，当载荷比例因子增加时，曲线是非线性的，一直达到极限，这种在结构发生变形一直到失稳，在变形后的位形上考虑平衡一直达到极限的方法称非线性屈曲或极限屈曲。

可见，工程实际中分枝屈曲现象实为罕见，它仅出现在完全无结构缺陷，完全沿轴向加压的绝对直杆及完整空球壳在均匀外压的情况下。

分枝屈曲现象虽然罕见，但实际中有不少结构屈曲状态接近分枝屈曲，而分枝屈曲的计算工作量又远小于计算极限屈曲的工作量，况且，不少作者得出结论，一些中等非线性的屈曲状态，可以用线性屈曲问题特征向疑的线性组合近似得到。

因此线性屈曲理论还是有其实际价值。

屈曲的含义可简述为：

结构处于一种平衡状态，载荷增量为一个微量，其位移增量很大。

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

[K#/0}=O

根据式（）和（），有

[K」=[Kj+[K订+[K」

（）

将式（）代入（）得

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

（［爲］+［心］”0}=0（〉

由式（）可以看岀，［K"］并不明显地包含位移增量〃0},在小变形情况下，该矩阵与应力水平成正比。

由于屈曲前线性假设，多数情况下，应力与外载也为线性关系，因此，如果令某一参考载荷对应的初应力刚度阵为［K；］,令屈曲极值载荷为与参考裁荷有

Fc=XFr关系，几称载荷比例因子，才为极值裁荷的比例因子。

极值载荷时的初应力刚度阵为

K；卜刃K；］（）

代入式（），则有

（［心］+/阳）〃0}=0（）

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

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

（1）按线弹性问题的有限元法形成各单元的刚度矩阵［心］・并用常规方法组装成结构刚度矩阵，即［K（）］=》k）］。

（2）对结构施加参考载荷并求解有限元方程［K.^}={Fr}^进而可求得各单元节点应力。

（3）对有膜应力存在的单元按式（）形式初应力矩阵\k；\,并组装成结构初应力矩阵,即

（4）求解广义特征值问题，即式（），得到最低几阶特征值石及对应的特征模态。

理论上它们都是平衡模态，当载荷达到分枝载荷人；7"时，平衡由一种模态“跳列”另一种模态，这也就是分枝屈曲的含义。

（5）从分枝载荷中挑选一个最小的载荷，即F；in=^minFr,作为极值载荷或临界载荷。

值得指出的是，由式（）所表征的线性屈曲问题，是建立在下述假设的基础上的，即假

设线性应变刚度矩阵在屈曲前不产生明显的变化，且初应力矩阵简单地与应力水平成正比。

如前所述，这在实际问题中是很罕见的。

由式（）所确左的极值载荷只能是近似的。

由于线性屈曲理论存在精度差，以及适用范用窄的限制，所以，在一般情况下，应当用切线刚度矩阵［Kj来研究这个问题。

当［Kr}l{3}=0时，发生随遇平衡。

显然，这里应该用逐次逼近的方法进行求解。

关于非线性屈曲问题的有限元解法，读者可参考其它文献。

板的大挠度及线性屈曲

基本问题

在板的大挠度问题中，板中的内力除弯曲内力外，还有薄膜内力，如图所示，在这种情况下，横向位移会引起薄膜应力，因此在板的大挠度问题中，面内变形和弯曲变形不再象小挠度板那样能分别处理，二者是相互耦合的。

和以前一样，平板的应变可用中面的位移来描述。

如果令x-y平面与中而一致，则广义应变和广义应力向量可以表示为

{b}=：

={NxNy

d2wd2\v

dxdy

（）

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

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

dudxw

-2

dyd

砒O

rul

1>IJ"Obo££广-J

式中第一项是我们熟知的线性项，而第二项则是非线性项。

如果所考虑的材料是线弹性的，则平板的弹性矩阵［d］由平而应力弹性矩阵［zr］和弯曲禅性矩阵［»什组成，即

（）

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

（）

<4=[“]｛砂

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

而

这样，形函数也应分成

（）

通过上述规定，除了非线性应变项｛$:

｝以外，其余各项都和标准的线性分析相同，不再重复。

应变矩阵屈的计算

为了确定单元的切线刚度矩阵［kT］.我们先来讨论单元的应变矩阵国］的计算。

首先应

注意到

式中

阪卜［崩+血］

（）

这里，［酬］和肉］分別是平而单元和弯曲单元按线性分析所得到的标准矩阵・而［璃是应

变的非线性项引起的，它可以通过g"对参数｛》"｝取微分来导出。

下面就来推导［b什的表

达式。

在式（）中.非线性应变分量可以方便地写成

人

1-2

空&酬¥

r1-」

——X

O鈿一勿空去亦一&O亦¥

-I丿\」

O亦一勿亦一&

釦一&O亦¥

■■

亦一办鈿¥

/•・w・、

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

{&}"監匸⑹旳

其中

6N：

6N；~dx

莎

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

将式（）取微分，有

〃《£}=£〃［人］{&}+g［勿{&}

注意到矩阵［A］和{&}具有两个有用的性质:

1.设{x}=［xpx2Jz是一个任意向量，则有

oXI

空去avv¥

zlxznx

r——丿

avaro

z（\

d[A^x}=<

av¥

z（\

于是有

（）

d[A^}=[A}l{0}

2.设{y}=[）w）J是一个任意向量，则有

〃{&}（）

儿J2

>1）?

亦一&亦¥

z<\/<\

32V-y

•123yyy亦T囂T&/（Xz（xJJ

/l\

d[A]{y}=<

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

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

（）

[b^\=[a]g]

单元切线刚度矩阵[灯]的汁算

根据式（），单元的切线刚度矩阵由三部分组成，即

[kT]=[k0]+[ka]+[kL]

（）

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

卧[卯kJ

将式（）中的和代入式（），可得大位移矩阵

2卩订01剧亂血1

（）

最后，利用式（）可求出初应力矩阵[©L为此，将式（）中的[bJ进行微分，得

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

［曲训备d［町

N、・

M、

利用式（），有

N,Nxv

6/［annv•=

最后可将初应力矩阵写成

式中

（）

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

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

板的屈曲问题

板的线性屈曲问题，可从式（）出发作为特征值问题考虑。

对于只在平面内承受压力载荷的情况，先求问题的弹性解，然后用式（）求出网］,而由线弹性方法求得［剧，此时［紅］=0,在完成结构刚度矩阵的组装后得到了广义特征值问题

（此］+才比弘/0}=0（）

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

上述平板的几何非线性有限元理论，也适用于由平板单元组成的壳体结构的几何非线性分析，只是多了一步将局部坐标系中的切线刚度矩阵及其它有关量转换到整体坐标系中的工作。

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

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

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

上廿平板大挠度问题中，非线性应变一位移关系是在特左的情况下建立的。

根据格林

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

用直角坐标表示的格林应变为

du1

=——+—

dx2

dudv一+一+dydx

dudu內內d\vdxdydxdyoxdy

其它四个应变只要通过轮换xTyTzTX/UTWT”，即可得到。

在小位移的情况下，忽略二次项，就得到线性应变公式。

这时，•和①是原来平行于坐标轴的线段的伸长率，而Yxy.Yyz和；^表示原来平行于坐标轴的两垂直线段夹角的变化。

但在应变较大时，上述几何意义就不再成立。

现在我们来推导三维应力状态下非线性的恨加1［妬］的一般计算公式。

稍加改变即可得到一维和二维情况下的计算公式。

对于板壳问题，可在特龙条件下忽略某些项后就可导出上节中的结果。

对杆系结构我们将在下节详细讨论。

应变矩阵［bJ的推导

变形体某点的应变应是线性应变｛吕）｝和非线性应变｛巧｝之和，即

（）

｛£｝=｛筍｝+｛勺｝

<§¥<§-azav一ax

加-axcg-內鈿一&+、■+•7+

av_az亦一办加一勿

、

按照式（），非线性应变可用下式表示

式中

仅卩

｛叮

Ikf

h}J

（）

研

耐

{&F

{叮

{須

阿］{&}=<（&}

〔{加

⑷是一个6X9的矩阵。

｛&｝是一个9X1的列向崑而

dvdwdxdxdx

余类推。

容易验证上述左义的正确性，并重新建立上节中矩阵［q］和｛&｝的两个性质。

我们再次有如｝=|幽］｛&｝+g［A｝1｛0｝=［A\l｛0｝（）

如果｛0｝用形函数［N］和节点位移向疑0｝来表示，则可写岀

®｝=［G］｛疔（）

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

将上式代入式（）得

冷｝=［A］［G00｝

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

［b」=［a］［g］（）

单元切线刚度矩阵［匕］的推导

注意到式（），即屈=個］+阴，由式（）积分得到［心］,而由式（）可求得［紅］,

又由式（），则有

kJ/仞｛b｝d

由式（）

注意到［G］不含0｝「，由上式得

将式（）代入（）的右端，则

d[Bj=[G]d[A]

化00『=JJJW〃［町何dv

同上一节一样，利用［A］和｛0｝的第二个性质，得

bj对

〃⑷{刃='TxyIbJ利r.Js

式中

由式（）可得

将式（）代入（）则

bj对rxJbj称

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

g阿［MlGk

（）

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

旧卜陶+［紅］+阳

杆系结构的大位移分析

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

对总体的拉格朗日列式法，我们将推导出桁杆单元和梁单元的切线刚度矩阵及其显式,其余il算均如肖所述。

而对于更新的拉格朗日列式法，则以刚架单元为例，重点介绍其计算

原理和实施步骤。

桁架单元的切线刚度矩阵

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

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

如果节点的位移比较大，则由于横向位移V会使单元发生附加的轴向伸缩。

这种附加的轴向应变与单元在变位过程中所转过的角度&有关，此时单元的轴向应变可表示为

du1,d\\y八

£.=—+—（一广（）

dx2dx

式中左端的第二项是考虑大位移时附加的非线性项。

显然上式也可由格林应变公式直接写出。

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

u=a.+a.x

12（）

v=a3+a4x

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

单元上任意点的位移可用杆端节点位移表示为

（）

式中

1--0

[N]"

01——

{/}=f=[W

再考虑到式（）,

（）

为桁杆单元的形函数矩阵。

将式（）代入式（）

有

于是

0d0Y+0-丄0丄{J}r0-丄0丄”0}。

LLLL

将上式与式（）比较，并考虑到式（），有

-丄0

一丄

-Pi

在式（）中

0为单元在变位过程中发生的转角，如图所示，于是式（）可以写成

（）

将式（）代入式（），可得线性分析时桁杆单元的刚度矩阵，如式（）所示。

将式（）和（）代入式（）可得桁杆单元的大位移矩阵如下

一&

-6

一少

0一&

-0-01

oe2

由式（）可以看岀，当单元处于变位以前的原始位置时,

0=0,因而大位移矩阵化］=0。

下而来推导桁杆单元的初应力矩阵。

由式（）并考虑式o

广

■

1（i

1/、

-1

>7倍_气）

[}

U：

—

J丿丿

的关系，有

［心

式（）与点有坐标无关，而对于桁杆单元来说，｛b『=b是一个常数，因此在积分时

这些项均可提到积分号之外。

若记N=gA为桁杆单元的轴向力，则有

-1

将上式与式（）比较可得

「00

0-1

（）

这就是桁杆单元的初应力矩阵，它计及了桁杆单元的轴向力对杆端横向位移的影响。

最后将式（），（）和（）代入式（），便得到桁杆单元的切线刚度矩阵。

另外应注意到,

展开阅读全文