第五章：大变形问题的有限单元法.ppt

1、第五章 大变形问题的有限单元法,1. 弹性大变形问题的有限元法 2. 弹性分支点稳定问题有限元分析 3. 物质描述大变形增量问题的T.L 、U.L法,大变形时平衡方程和虚位移原理,变形体初始和现时位形如图所示，以欧拉应力表述平衡时,这是现时位形空间描述的平衡条件。,在外荷为保守力系,时,经推证得,上式乘以J，由于,由复合函数求导数，平衡方程改写为,平衡方程成为,对任意j恒有,证明,以j=1为例,同理可验证，对任意j恒有,上述拉格朗日和克希荷夫应力表示的平衡条件都是以初始位形作参考的物质描述。,利用克希荷夫应力和拉格朗日应力的关系，可将平衡条件改为,如果考虑到变形梯度和位移梯度间的关系，,对比小

2、变形情况，可见大变形时变形对平衡的影响，是通过变形或位移梯度表现出来的。,则克希荷夫应力表达的平衡方程,可改为,设现时位形微小虚位移在V内单值连续、在位移边界上为零。则外力总虚功为,考虑到位移边界虚位移为零和应力边界条件,空间描述的虚位移原理,将平衡方程引入，考虑到虚位移微小，则,利用格林公式，可得,也即虚位移原理的虚功方程为,若虚位移用虚速度、虚应变用虚应变率代替，则,柯西应变,虚功率方程,为建立物质描述虚功方程，先讨论能量共轭关系。,空间描述中,又因为,变形率张量是相对现时位形定义的柯西应变的速率,单位体积变形功率,因此克希荷夫应力和格林应变在能量上共轭。,由于变形率和欧拉应力张量是对称的

3、，因此,再利用欧拉和拉格朗日应力间关系，可得,因此拉格朗日应力和初始位形位移梯度在能量上共轭。,由此即可得到物质描述的虚功方程为,由虚位移原理可得,式中 为单元结点力矩阵， 为单元等效结点荷载矩阵,Pd和PE分别为直接和等效结点荷载矩阵，R为综合等效结点荷载矩阵。上式也可写为,按集成规则集装后可得,1. 弹性大变形问题的有限元法,弹性大变形问题，需要考虑变形的非线性项和变形对平衡的影响。,若以初始自然平衡状态作初始位形，则物质描述的格林应变为,式中,线性部分,非线性部分,为便于计算机编程，将张量转换为矩阵：,格林应变矩阵和张量的分量间有如下关系,对应的克希荷夫应力矩阵和张量分量间关系为,引入两

4、个算子矩阵,式中,再引入位移梯度向量的记号,3阶单位矩阵,在上述符号基础上，格林应变有位移表为,则单元格林应变为,其中线性应变矩阵B和线性分析一样,设单元位移场为,非线性部分“应变矩阵”为,式中G为如下93m的矩阵,式中AL为如下69的矩阵,由(AL)可见，格林应变-位移关系是非线性的。,非线性部分因为,为用虚位移原理建立单元刚度方程，还需有应变增量和位移增量间的关系。对线性部分,所以,综上所述，格林应变增量为,如果记,，则 。,对弹性问题，在物质描述下本构关系(克希荷夫应力和格林应变)为,在上述基础上，由虚位移原理可得,式中 为单元结点力矩阵， 为单元等效结点荷载矩阵,Pd和PE分别为直接和

5、等效结点荷载矩阵，R为综合等效结点荷载矩阵。上式也可写为,按集成规则集装后可得,将克希荷夫应力表达式代入，可得,式中K(U)是非对称的，为,对非线性弹性问题 、 和 都是位移的函数。,根据非线性方程切线刚度矩阵的定义，可得,根据本构方程，则有,又因,式中,所以,又因 ,所以,由此可得,式中,基于上述说明，可得,如果引入如下记号,初应力或几何刚度矩阵,线弹性刚度矩阵,大位移刚度矩阵,则单元切线刚度矩阵为,“结构”切线刚度矩阵为,建立了切线 刚度矩阵， 用牛顿法等 即可求解。,以上的讨论没涉及具体单元，因此具有普遍性。具体单元分析时，因形函数一般是自然坐标的函数，故需作坐标变换后代入相关公式，从而

6、建立具体单元的切线刚度矩阵。,牛顿迭代法解大变形问题的具体步骤为：,牛顿迭代法解大变形问题的具体步骤为：,以上是全量形式的弹性大变形分析，为了保证收敛，拟用增量迭代法。,极值点失稳问题,弹塑性问题,2. 弹性分支点稳定问题有限元分析,研究弹性稳定问题时材料应力应变关系是弹性的，也即,。式中,对线弹性问题由弹性常数表示，,对非线性弹性问题由割线或切线系数表示。,稳定问题可分为两大类：分支点（第一类）稳定问题；极值点（第二类）稳定问题。分支点稳定问题存在稳定平衡状态和不稳定平衡状态的分界点，在此临界状态下，系统可在两种不同的变形形式下保持平衡。极值点（第二类）稳定问题变形形式整个过程中是不变的，但

7、当荷载达到临界值时，变形将迅速增加使系统丧失承载能力或破坏。只讨论第一类线性弹性稳定问题。,2. 弹性分支点稳定问题有限元分析,对于分支点稳定问题时，可设外荷是比例加载的。由于第一类稳定问题主要是解决临界荷载计算，因此失稳状态相对稳定平衡状态的位移仍然是微小的，故可略去上节所推得的大位移矩阵项。但因为现时失稳位形和初始稳定平衡是两种不同的变形形式，因此应变要考虑失稳位移的非线性项。临界状态时稳定平衡和不稳定平衡荷载不变，前者是没有与失稳位移对应的荷载的，由结构力学已知，稳定分析的关键是确定几何刚度矩阵。因此集装后的整体切线刚度方程为,式中,为整体线性刚度矩阵，,为整体几何刚度矩阵，,为不稳定平

8、衡位置相对稳定平衡位置的整体位移。,2. 弹性分支点稳定问题有限元分析,对分支点稳定问题，关键是建立几何刚度矩阵。此时，以失稳前的平衡位置作初始位形，以失稳形态作现时位形。,和前述弹性大变形不同的是：,大位移矩阵可忽略。,应变仅包含失稳位移的非线性项。,分支点处相应失稳位移的综合荷载为零。,失稳前变形是微小的。,2. 弹性分支点稳定问题有限元分析,由于是比例加荷，在基准荷载下初始位形（稳定平衡位置）时的应力可由线性分析求得，因此上节几何刚度矩阵式子里的S,与荷载比例系数成正比。,在比例加载时,为荷载系数，,为单元基准荷载下的几何刚度矩阵。,由此可见，第一类稳定问题最终归结为求解一个特征值问题。

9、,解得特征值后，即可得到临界荷载，一般只关心最小临界荷载。,设变形前单元长度为l，截面积为A，弹性模量为E。单元杆端位移矩阵为e,式中,2.1 桁架单元,。单元位移为,式中 为失稳位移（出平面的位移）。基于此，单元格林应变为（针对失稳前定义）,克希荷夫应力为,因分支点稳定关心的临界荷载，临界荷载时平衡有两重性，故临界状态应力等于失稳前状态的应力，也即,基于上述结果，单元几何刚度矩阵为,，失稳前应力为,设变形前单元长度为l，截面积和惯性矩为A、I，弹性模量为E。单元杆端位移矩阵为e,单元位移为,式中,2.2 梁单元,式中Ni和线性有限元梁单元一样。梁单元应变为,其中第一项为线性有限元里的线性项，

10、非线性的第二项为,基于上述结果，单元几何刚度矩阵为,象桁架单元说明一样，单元应力为,3.1 物质描述大变形增量问题的T.L法,对弹塑性、粘性-蠕变和施工力学等问题，介质的反应和变形的历史有关。对随时间变化的荷载，需要将时间变量离散成序列，以求解各时刻的响应。为此，都需要用增量法来解决。,从t到t+t的增量期间进行物质描述求解时，一般可选两种参考位形：初始和t时刻的位形。前者称为全拉格朗日(T.L)表述，后者称为修正拉格朗日(U.L)表述。,设从0到t时刻的全部反应、位形均已求得，现在的问题是，如何求t+t时刻的响应。,设t0、t和t+t的物理量分别用如下符号标记,对T.L法，介质位移是初始位形

11、坐标的函数,坐标 ，密度,面积和体积,设有限元分析时单元形状描述为,又设有限元分析时单元位移场为,式中雅可比矩阵J为,象线性有限元等参元分析一样，由于形函数一般是对自然坐标 定义的，因此有限元分析中的对坐标求导等，应象线性有限元一样进行转换。,在上述记号下，格林应变为,时刻t t +t的应变增量为,式中,t +t时刻,t时刻,由于，u (t时刻)在增量步内已知，因此,同大变形有限元，将张量转换为矩阵，则,引入大变形所用算子记号，则有,由于u在增量步内已知，因此,B和BL都是已知的。又若记,则有,因为 ，因此,综上可得,为进行有限元列式，还需讨论克希荷夫应力。设t和t+t时刻的应力分量分别为,基

12、于上述分析，利用t+t时刻的虚位移原理虚功方程,则象应变分析一样，可将 分成 。同样换为矩阵表示，则有,再次强调，t时刻及其前的量都是已知的，因此变分为零。基于此,将此结果代入虚功方程，可得单元刚度方程,式中 和 是对初始位形定义的，t+t时刻的体积力和表面力，它们是已知的。,式中,或,将其按集成规则集装后可得,再引入如下记号,将 和 的表达式代入，可得,或,利用这些关系，非线性平衡方程可写为,为求解上述方程，尚需解决如下两方面问题,首先假设,然后将S和E的关系线性化。根据本构关系,则有,为使其线性化，设( t时刻的材料性质矩阵),在做了上述两方面处理后，可得,由此出发，用非线性方程的相关解法

13、，即可解决大变形非线性（材料非线性）问题。,将其代回非线性平衡方程，可得,3.2 物质描述大变形增量问题的U.L法,因t+t的位移是用t时刻位形为基准度量的，因此,在t,t+t间隔内，以t时刻位形为参考位形，其增量位移为,象T.L法一样，设单元和位移的描述为,但需指出的是，式中形函数N是t时刻单元自然坐标的函数。在计算 等导数时，要先作坐标变换(Xi应换为xi)。,类似地，用矩阵来表示则有,在时刻t和t+t的格林应变是以t时刻位形定义的，因而它们分别为,式中算子矩阵象T.L法一样，但应将Xi换为xi。,式中算子符号象T.L法一样，但应将Xi换为xi。,象T.L法一样可导得,关于应力的处理也和T

14、.L法一样，对t 时位形定义的t和t+t时刻的克希荷夫应力分别为,象T.L法一样由虚位移原理虚功方程可导得,象T.L法一样推导，引入如下符号定义,式中,其中 分别为t时刻位形定义的单元体积、表面应力、体力和表面力。,则可得,象T.L法一样，为求解上述方程也需解决线性化问题。首先讨论S的计算。因为,因此,可改写为,由第四章已知,式中,可得,由如下两式消去 并利用 ,且注意到Vij对称、ij反对称,又由于 、 和 ，因此,上式最后一项将使本构张量不对称，对金属类不可压缩介质，这一项可略去，也即,在有限的克希荷夫应力和格林应变增量之间仍认为,这就是U.L法的本构关系线性化。,与T.L法一样，除本构关

15、系线性化外，还需解决几何方面的线性化。因为,为对上述三项积分作几何方面的线性化,这一关系可有两种矩阵表达方式：其一是把 写成六维向量形式,以B代替 进行几何非线性的线性化,其另一方式是把 写成九维向量形式,引入算子矩阵 ，则有,第一个积分用,作第二、三个积分时,基于上述讨论，积分 时，作如下处理,式中 为由 组成的66矩阵， 矩阵为,结 束,则t+ t时刻的非线性平衡方程,可改写成,要提高精度可减小t或采用适当的修正技术。,两种方法的比较 在大变形增量问题有限元计算中，用那种方法，主要根据计算的效率。,在线性化后的总刚度矩阵中， 法的 里包含 ，而它包含位移效应，因此它是满阵。而在 法中与其对应的 在线性化后，用的是B ，它不包含位移效应，所以计算比 法节省机时。,在 法中，形函数的导数仅与初始坐标有关，因此只要第一步计算后存起来即可。而 法，因为是对