第09章 几何非线性

1第9章几何非线性在以前各章所讨论的问题中，都是基于小变形的假设。它包含了两方面内容：一是假定物体所发生的位移远小于物体自身的几何尺度，在此前提下，建立结构或微元体的平衡条件时可以不考虑物体的位置和形状（简称位形）的变化，因此分析中不必区分变形前和变形后的位形，即如我们通常习惯上所用的以变形前描述位形去表达变形后的平衡位形；二是假定在加载和变形过程中的应变可用一阶微量的线性应变进行度量，即应变与位移成一阶线性关系。实际上我们会碰到很多不符合小变形假设的问题，例如板、壳等薄壁结构在一定载荷作用下，尽管应变很小，甚至远未超过弹性极限，但是位移很大，材料的线元素有较大的转动；又如压杆失稳后的变形分析，这时必须考虑变形的位形对平衡的影响，即平衡条件必须建立在变形后的位形上；再如薄板的大挠度问题理论中，由于考虑了中面内薄膜力的影响，可能使得按小挠度理论分析得到的挠度有很大程度的缩减。本章讨论的就是放弃小变形假设的非线性问题即几何非线性问题，它分为小变形几何非线性和大变形几何非线性问题。9.1小变形几何非线性有限元方程的建立与求解回顾前面介绍过的建立有限元单元平衡方程时采用应变能泛函求极值。单元体的应变能等于外力所做的功与单元变形所产生的应变能之和。如果用F表示内力和外力矢量的总和，那么由位移增量表示的应变能泛函可以写成如下表达式ddVd=-∫TTεσuFP（9-1）式中F代表所有载荷列阵；du为位移列阵；dε为应变列阵,σ为应力矩阵。如果我们用应变的向量形式写成位移增量和应变增量的关系dd=εBu将上式代入式（9-1），并对应变能泛函求极值，得到非线性问题的一般平衡方程式()dV=-=∫ysTuBF0（9-2）上式中的积分运算是应用通常的方法，由各个元素的积分对于节点平衡所作的贡献总合而成。式（9.2）的推导，仅应用了最小势能原理，没有涉及材料、几何是线性还是非线性的问题，因此不论位移（应变）是大还是小，都应完全适用。在线性有限元中，由于应变与位移之间及应力与应变之间的关系均是线性的，即==εBuσDε将上述关系代入式（9-2），因此可以得到dV-=∫TBDBuF02现在考察小变形几何非线性问题，应变和位移的关系是非线性的，因此矩阵B是u的函数，为简便起见，可以写成()=+0LBBBu（9-3）式中0B作为线性应变分析的矩阵项；LB是大位移应变矩阵项，并取决于u，它是由非线性变形引起的。一般地说，LB是位移列阵u的线性函数。而应力应变关系还仍然认为是线性弹性关系，则有()=+00σDε-εσ式中D为材料弹性矩阵；0ε为初应变列阵；0σ为初应力列阵。很显然，式（9-2）的解可以通过迭代方法求得。如果采用牛顿-拉斐逊（Newton-Raphson）方法，必须寻求dψ和du之间的关系。通过（9-2）式，取ψ的微分，于是有()dddVddV=+∫∫yssTTuBB（9-4）如不考虑初应变和初应力的影响，得到dddse=D=DBu由式（9-3）可知dd=LBB所以式（9-4）可以改写为()dddVdd=+=+∫sysTLΒKuKKu（9-5）式中dV==+∫T0LKBDBKK，且0K表示通常的小位移动线性刚度矩阵，即dV=∫T000KBDB，()dV=++∫TTTL0LLLL0KBDBBDBBDB，表示大位移所引起的影响矩阵，一般称大位移矩阵或初位移矩阵；dddV=∫ssTLKuB，表示应力水平的对称矩阵，称它为初应力矩阵或几何刚度矩阵。这样，完成了小应变几何非线性问题的非线性方程组的建立。下面我们应用大挠度板单元的例题对小变形几何非线性问题给予具体的说明。在大挠度范围内，当平板承受横向载荷时，板中的内力除弯曲内力外，还有薄膜内力。众所周知，横向位移可以引起薄膜应变。所以在大挠度的情况下，平面变形和弯曲变形不再认为是互不相关，而是相互耦合的。如果OXY平面与平板中面重合，平板应变可以用中面位移描述。令ε表示中面应变和曲率列阵，σ表示薄膜内力和弯曲内力列阵。于是有3222222⎡⎤∂∂∂⎡⎤==⎢⎥⎣⎦∂∂∂∂⎣⎦⎡⎤⎡⎤==⎣⎦⎣⎦TTplbxyzTTplbNxNyNxyxyxyεσ（9-6）式中下标“pl”表示中面的薄膜应变，“b”为中面的弯曲应变。由板壳大挠度理论可知，平板挠度w对于中面线所在x和y方向产生附加伸长和附加角变形为221212uwxxvwyyuvwwyxyxeeg∂∂⎛⎞=+⎜⎟∂∂⎝⎠⎛⎞∂∂=+⎜⎟∂∂⎝⎠∂∂∂∂=++∂∂∂∂xyxy将上式代入式（9-6）可得2222222121200002uxwvxywuvyyx∂⎧⎫⎪⎪∂⎧⎫⎪⎪∂⎛⎞∂⎪⎪⎜⎟⎪⎪∂⎝⎠⎪⎪⎪⎪∂⎪⎪⎪⎪⎛⎞∂∂∂⎪⎪⎪⎪⎜⎟+⎪⎪∂⎪⎪⎝⎠∂∂⎧⎫⎧⎫⎪⎪⎪⎪⎪⎪⎪⎪=+=+⎨⎬⎨⎬⎨⎬⎨⎬∂∂∂⋅⎪⎪⎪⎪⎪⎪⎪⎪⎩⎭⎩⎭∂∂⎪⎪⎪⎪∂⎪⎪⎪⎪∂⎪⎪⎪⎪∂⎪⎪⎪⎪⎪⎪⎪⎪⎩⎭∂⎪⎪⎪⎪∂∂⎩⎭0Lplpl0bεεεε（9-7）式中右端第一项表示线性应变，第二项就是给出的非线性项。如果所考虑到材料是弹性的，于是平板弹性矩阵D是平面应力弹性矩阵pD和弯曲弹性矩阵bD所组成，即⎡⎤=⎢⎥⎣⎦pbD0D0D式中21111hmmmm⎡⎤⎢⎥=⎢⎥-⎢⎥⎢⎥⎣⎦pD00-0024()321121Ehmmmm⎡⎤⎢⎥=⎢⎥-⎢⎥⎢⎥⎣⎦bD0101-002平板弯曲的位移模式可以利用适当的形函数通过节点位移列阵来表示，即[]uvw=TeNu为方便起见，节点位移列阵也可以区分为平面的和弯曲的两类：wwuvwxy⎡⎤∂∂=⎢⎥∂∂⎣⎦Tiiiiiiu同样形函数也可以表示为⎡⎤=⎢⎥⎢⎥⎣⎦piibiN0N0N至此，我们可以计算B矩阵：=+0LBBB这里⎡⎤⎡⎤=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦pb0L0Lb0B00BBB=000B（9-8）式中p0B，b0B就是平面元素和弯曲元素按线性分析所得到的标准矩阵，而bLB可以由eLp通过节点弯曲位移列阵bu来确定。也就是说，bLB是由于应变的非线性项而引起的。在式（9-7）中，应变分量的非线性部分可方便地写成011022wwxwxwyywwyx⎡⎤∂⎢⎥∂⎧⎫∂⎢⎥⎪⎪⎢⎥∂∂⎪⎪==⎨⎬⎢⎥∂∂⎪⎪⎢⎥∂⎪⎪⎢⎥∂∂⎩⎭⎢⎥∂∂⎢⎥⎣⎦eLpCθ（9-9）式中00Twwxywwyx∂∂⎡⎤⎢⎥∂∂⎢⎥=∂∂⎢⎥⎢⎥∂∂⎣⎦C，wwxy⎡⎤∂∂==⎢⎥∂∂⎣⎦eTbθGu，xxyy⎡⎤∂∂⎢⎥∂∂⎢⎥=⎢⎥∂∂⎢⎥∂∂⎢⎥⎣⎦LLbbjibbjiNNGNN，矩阵G完全由单元坐标所定义。可以证明：a）ddq⋅=⋅CCq5b）若设[]123yyy=Ty的列阵，则1332yyddyy⎡⎤⋅=⎢⎥⎣⎦TCyθ由上述结果，我们可以推导bLB的表达式。取式（9-8）的微分，可得1122ddddd=⋅+⋅==eLbpεCCCCGuqqq于是由LB定义，可以直接写出bLB=CG由于对非线性方程求解经常采用牛顿-拉斐逊方法，需要求出切线刚度矩阵，如式（9-5）所示的sTK=K+K。现在计算板单元的切线刚度矩阵TK，对于线性的小变形的刚度矩阵⎡⎤⎢⎥⎢⎥⎣⎦pl00b0K0K=0K在以前已有定义。对于初始位移矩阵，可以把式（9-8）代入式（9-5）得到()()dxdy⎡⎤⎢⎥=⎢⎥⎣⎦∫TTLplplbbplb0LLL0KBDBBDB对称最后，利用式（9-5）可以算出几何刚度矩阵：⎡⎤⎢⎥⎣⎦ssb00K=0K式中FFdxdyFF⎡⎤=⎢⎥⎢⎥⎣⎦∫sNxNxybTNxyNyKGG牛顿-拉斐逊迭代法如同第7章所讨论的，其基本步骤是：首先按小变形线性理论算出位移u的一阶近似值，它是平面应力和弯曲应力问题的耦合解。然后，按所求得的位移近似值，利用式（9-7）求得应变值，包括对线性部分和非线性部分的贡献，通过材料的弹性矩阵决定相应的应力值。再由式（9-2）就可以算出()yu的一阶近似值y1，为了下一次迭代，算出平板的切线刚度矩阵，接下去就可以利用牛顿-拉斐逊方法进行多次迭代，直到yn足够小为止。9.2大变形几何非线性的几何描述6涉及大变形几何非线性问题的有限单元法中，结构的位移随加载过程不断变化，因此必须建立参考位形，以参考位形表征其变化过程，目前对此基本上采用两种不同的表达格式，第一种格式中所有静力学和运动学变量总是参考于初始位形，即在整个分析过程中参考位形保持初始位形不变，这种格式称为完全的拉格朗日（Lagrangian）格式；另一种格式中所有静力学和运动学的变量参考于每一载荷或时间步长开始时的位形，即在分析过程中参考位形是不断地被更新的，这种格式称为修正的或更新的拉格朗日格式。按上述两种格式考察，可以发现大变形几何非线性问题的几何描述必须建立两个坐标系，并对两种描述格式予以定义。设运动物体中的任一点在t=0时刻（初始位形）的空间位置用空间笛卡儿坐标系的一组坐标Xi来表示，而该点在t时刻的空间位置用一组坐标ix表示（如图9-1），则x3,X3x1,X1x2,X2•Xi•xiV0tt=00V图9-1坐标描述(),iijxxXt=（9-10）或□□□□(),iijXXxt=□□□□□□□□□□□□□□□□□□□□□（9-11）如果对于物体内所有点都能给出这样的方程，那么整个物体的运动也就被确定了。由于我们所考察的物体运动和变形等参数都是随时间连续变化的，因此在度量物体运动和变形时需要选定一个特定时刻的状态作为基准，即称为参考状态。为了正确定义描述，我们规定初始时刻任意点的坐标iX称为物质坐标。同一质点在t时刻的坐标xi称为空间坐标，同时定义以物质坐标为自变量，描述物体运动和变形的方法（即如式（9-10）所示）称为拉格朗日描述，或物质描述。反之，如果以空间坐标作为自变量（或参考系）加以描述（即如式（9-11）所示）称为欧拉描述，或称空间描述。假设物体的运动和变形是时间的单值连续函数，从几何观点看，式（9-10）表示一个物体从初始状态所占据的区域V0到现时（t时刻）状态所占据区域V的映射。由于是单值连续假定，这个映射是一一对应的。即在V0内处处有71111232221233331230xxxXXXxxxJXXXxxxXXX∂∂∂∂∂∂∂∂∂=≠∂∂∂∂∂∂∂∂∂定义J为变形梯度矩阵的行列式，xX∂∂ij为变形梯度。如果初始状态中由OABC组成的平行六面体体元，其中的三条棱是OA、OB、OC，对应线元分量是dXp、Xdl、XΔm（如图9.2所示），在映射下（如式（9-10）），它们的现时线元是oa、ob和oc，这些线元的分量()(),,dxxXdXtxXt=+-iijjij123123xxxxdXdXdXdXXXXXxxXXxxXXdd⎫∂∂∂∂=++=⎪∂∂∂∂⎪⎪∂⎪=⎬∂⎪⎪∂⎪Δ=Δ∂⎪⎭iiiijjiijjiijj同理（9-12）X1x1x2x3X2X3AOCBaobc0t=t而由它们构成的六面体体元的体积为123123123dxdxdxdVxxxedxxxxxxdddd==ΔΔΔΔijkijk图9-2变形前后的体元8xxxedXXXXXXd∂∂∂=Δ∂∂∂jikijkplmplmJedXXXJdVd=Δ=plmplm0（9-13）式中eijk为排列张量。按排列张量性质，当i、j、k按偶排列时，取值为1；奇排列时取值为－1；其余为0。而dV0为123123123dXdXdXdVXXXedXXXXXXdddd==ΔΔΔΔ0plmplm由此可知，对于不可压缩物体1dVJdV==0；如物质运动，遵守质量守恒定律，则dVJdVrr==00。下面讨论变形过程中面元的大小和方程的变化。取任意两点间的线元（图9.3所示）dXi、Xdi，变形后的现时状态为dxi、xdi，并用NdAi0及ndAi表示变形前后的面元法线方向和面积大小,则变形前后的有向面元分别是dXi和Xdi、dxi和xdi的矢量积，可表示为：X1x1x2x3X2X3ABOabdA0odAnN0t=tNdAedXXndAedxxdd==i0ijkjkiijkjk将上述第二式两端左乘以xX∂∂ip，则由式（9-12）可得xxxxndA