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

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

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

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

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

=0（）

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

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

>

将式（）代入（），则有

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

VV

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

〃何=［刃忖（）

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

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

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

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

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

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

〃何=［切［科何（）

将式（）代入式（）后得

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

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

ffj阳何dv=（fJJ［Bor［Z）K>

+（J『［Bj［D血>

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

若记

V

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

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

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

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

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

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

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

I

上式中［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

根据式（）和（），有

将式（）代入（）得

在线性屈曲情况下，屈曲前结构处于原始位形的线性平衡状态，因此上式中的大位移矩阵［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平面与中而一致，则广义应变和广义应力向量可以表示为

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

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

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

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

什组成，即

（）

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

II

<

4=[“]｛砂

W

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

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

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

应变矩阵屈的计算

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

首先应

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

变的非线性项引起的，它可以通过g"

对参数｛》"

｝取微分来导出。

下面就来推导［b什的表

达式。

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

O亦一勿亦一&

釦一&

O亦¥

=

■■

亦一办鈿¥

/•・w・、

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

OW

{&

}"

監匸⑹旳

其中

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

将式（）取微分，有

〃《£

}=£

〃［人］{&

}+g［勿{&

}

注意到矩阵［A］和{&

}具有两个有用的性质:

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

于是有

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

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

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

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

[b^\=[a]g]

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

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

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

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

卧[卯kJ

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

最后，利用式（）可求出初应力矩阵[©

L为此，将式（）中的[bJ进行微分，得

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

利用式（），有

N,

6/［annv•=

最后可将初应力矩阵写成

式中

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

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

板的屈曲问题

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

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

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

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

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

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

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

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

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

根据格林

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

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

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

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

这时，•和①是原来平行于坐标轴的线段的伸长率，而Yxy.Yyz和；

^表示原来平行于坐标轴的两垂直线段夹角的变化。

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

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

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

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

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

应变矩阵［bJ的推导

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

｛£

｝=｛筍｝+｛勺｝

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

⑷是一个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｝「，由上式得

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

旧卜陶+［紅］+阳

杆系结构的大位移分析

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

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

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

原理和实施步骤。

桁架单元的切线刚度矩阵

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

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

du

Hx

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

这种附加的轴向应变与单元在变位过程中所转过的角度&

有关，此时单元的轴向应变可表示为

du1,d\\y八

£

.=—+—（一广（）

dx2dx

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

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

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

u=a.+a.x

12（）

v=a3+a4x

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

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

于是

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

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

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

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

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

-1

将上式与式（）比较可得

「00

0-1

00

01

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

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

另外应注意到,

由于［灯］和［©

］均为对称矩阵，所以桁杆的切线刚度矩阵旧］也是对称矩阵，由［約］所

组装的结构刚度矩阵［K」也是对称的，这对求解是很有利的。

刚架单元的切线刚度矩阵

在大位移的情况下，仍可假肚刚架单元的位移函数与小位移情况下的相同，即单元的轴向位移U是局部坐标X的线函数，而横而位移”则是X的三次函数。

用式子表示如下

u=q+a2x

v=+aAx+

由杆两端的节点位移条件，可从上式解得于是单元上任意点的位移可用杆端位移

表示成

{/}=<

%MW

［吩

1-3（歹+2（汀

X

00-00

L

兰0

03（-）2-2（-）3x（--l）（-）LLLL}

1-丄

01-3（-）2+2（-）3x（l--）20

LLL

哙2-吩吩-咛

于是［N］可表示成

M=

而式（）可写成

假若单元的曲率仍能用横向位移V的二阶导数近似地表示，则代表弯曲所引起的轴向应变的项仍与线性分析时相同。

而由于横向位移引起的单元轴向应仍应按式（）il•算，此时有

式中代表由线位移所引起的轴向应变，/代表由弯曲所引起的轴向应变。

而

为非线性项。

由式（）和（）得

将上式取微分，并注意到［nJ和［nJ不含有

由式（）和（）可得

{£

l}=

将上式取微分，有

dx

L2

6x6x2

~73+ir

3x22x

由式（）可得

〃｛£

｝=〃低｝+〃伍｝

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

d{£

}=

d{3}e+

将上式与式（）比较，并注意到式（）的关系有

［讣

将式（）

将国）］代入式（）就得到线性分析时的单元刚度矩阵［&

］,如式（）所示。

和（）代入式（），就可以得到单元的大位移矩阵

A

A-

-~B

0-/

-刁0

、ic

EA

fL

-B

B0

■<

Jo

0-

■A

B

0A

-c

-A

A0

C

-C0

牙2

—、

-AB0一"

AC

-AB

—TB「

0AB

-BC

（C

—•>

AB

—2

-AC

^BC0-AC

dx+L

A=Aj2（Vj-Vi）+AlBl0i-A[C0

B=AeOj-*）+-BiC®

>

C=A】C]（i,—匕）+dC0—0■

4x

3L-L

1}

2x

3x2

13

B严1

G=

q◎分别为单元/端和丿•端的截而转角。

在计算式（）中的积分时，可以应用下列积分公式

A^dx=—B;

dx=—L

J7AC//.v=丄

Joloj

Jo5厶Jo15

“1rL1

[A^Bxdx=—一\B,adx=—L

Jo10J（）30

以及

^B,dx=冬

Jo1135/?

「B^C.dx=—

Jo「140

A:

B~dx=—

Jo1135厶

超C4=o

CB^dx=—L

Joi35

Jo1135L2