大变形问题的有限单元法讲义

1、第五章 大变形问题的有限单元法1. 弹性大变形问题的有限元法2. 弹性分支点稳定问题有限元分析3. 物质描述大变形增量问题的T.L 、U.L法2000.41哈尔滨建筑大学 王焕定教授制作1. 弹性大变形问题的有限元法 弹性大变形问题，需要考虑变形的非线性项和变形对平衡的影响。 若以初始自然平衡状态作初始位形，则物质描述的格林应变为式中线性部分非线性部分2000.42哈尔滨建筑大学 王焕定教授制作 为便于计算机编程，将张量转换为矩阵： 格林应变矩阵和张量的分量间有如下关系对应的克希荷夫应力矩阵和张量分量间关系为 引入两个算子矩阵2000.43哈尔滨建筑大学 王焕定教授制作式中再引入位移梯度向量的

2、记号3阶单位矩阵2000.44哈尔滨建筑大学 王焕定教授制作 在上述符号基础上，格林应变由位移表为则单元格林应变为其中线性应变矩阵B和有限元(I)一样设单元位移场和有限元(I)一样为非线性部分“应变矩阵”为2000.45哈尔滨建筑大学 王焕定教授制作式中G为如下93m的矩阵2000.46哈尔滨建筑大学 王焕定教授制作式中AL为如下69的矩阵2000.47哈尔滨建筑大学 王焕定教授制作由(AL)可见，格林应变-位移关系是非线性的。非线性部分因为所以 为用虚位移原理建立单元特性方程，还得建立应变增量和位移增量间的关系。对线性部分 综上所述，格林应变增量为验证如果记，则 。2000.48哈尔滨建筑大

3、学 王焕定教授制作 对弹性问题，在物质描述下本构关系为 在上述基础上，由虚位移原理可得式中 为单元结点力矩阵， 为单元等效结点荷载矩阵由于应变、应力以矩阵表示，因此弹性矩阵 应按下式并考虑应变矩阵定义来建立2000.49哈尔滨建筑大学 王焕定教授制作Pd和PE分别为直接和等效结点荷载矩阵，R为综合等效结点荷载矩阵。上式也可写为将克希荷夫应力表达式代入，可得 按集成规则集装后可得2000.410哈尔滨建筑大学 王焕定教授制作根据本构方程，则有式中K(U)是非对称的，为对非线性弹性问题 、 和 都是位移的函数。 根据非线性方程切线刚度矩阵的定义，可得又因B、G为已知矩阵2000.411哈尔滨建筑大

4、学 王焕定教授制作式中所以又因 ,所以2000.412哈尔滨建筑大学 王焕定教授制作由此可得式中 基于上述说明，可得2000.413哈尔滨建筑大学 王焕定教授制作如果引入如下记号则单元切线刚度矩阵为初应力或几何刚度矩阵线弹性刚度矩阵大位移刚度矩阵“结构”切线刚度矩阵为建立了切线刚度矩阵，用牛顿法等即可求解。2000.414哈尔滨建筑大学 王焕定教授制作 本节的讨论没涉及具体单元，因此具有普遍性。具体单元分析时，因形函数一般是自然坐标的函数，故需作坐标变换后代入相关公式，从而建立具体单元的切线刚度矩阵。需要指出的是 建议自行对各种单元自行推导切线矩阵。 本节只讨论了全量形式的弹性大变形分析，具体

5、求解步骤如讲义所示。为了保证收敛，拟用增量迭代法。对第二类稳定问题（极值点失稳问题）、弹塑性问题等，必须用 3. 所介绍的增量形式来解决。To26在我们、王勋成、谢贻权、徐次达等的有限元教材中，都有一些具体单元的切线刚度矩阵，需要时可供参考。2000.415哈尔滨建筑大学 王焕定教授制作2. 弹性分支点稳定问题有限元分析 对分支点稳定问题，关键是建立几何刚度矩阵。此时，以失稳前的平衡位置作初始位形，以失稳形态作现时位形。 和前述弹性大变形不同的是：大位移矩阵可忽略。应变仅包含失稳位移的非线性项。 分支点处相应失稳位移的综合荷载为零。 注意到上述差异，即可用上节结果解决分支点稳定的有限元分析。失

6、稳前变形是微小的。2000.416哈尔滨建筑大学 王焕定教授制作 设变形前单元长度为l，截面积为A，弹性模量为E。单元杆端位移矩阵为e式中2.1 桁架单元。单元位移为式中 为失稳位移。基于此，单元格林应变为2000.417哈尔滨建筑大学 王焕定教授制作 克希荷夫应力为 因分支点稳定关心的临界荷载，临界荷载时平衡有两重性，故临界状态应力等于失稳前状态的应力，也即 基于上述结果，单元几何刚度矩阵为 ，失稳前应力为2000.418哈尔滨建筑大学 王焕定教授制作 设变形前单元长度为l，截面积和惯性矩为A、I，弹性模量为E。单元杆端位移矩阵为e单元位移为式中2.2 梁单元2000.419哈尔滨建筑大学

7、王焕定教授制作 式中Ni和有限元(I)一样。梁单元应变为其中第一项为有限元(I)里的线性项，非线性的第二项为 基于上述结果，单元几何刚度矩阵为 象桁架单元说明一样，单元应力为同结构力学2000.420哈尔滨建筑大学 王焕定教授制作 有限元(I)里的二维问题单元位移为2.3 二维单元(薄板稳定)但失稳时的位形，将有出平面的位移w，和杆单元一样，应变需考虑w及其非线性项非线性项2000.421哈尔滨建筑大学 王焕定教授制作 设出平面位移为面内弹性应力为 对出平面位移，其G矩阵为 对应的M矩阵为2000.422哈尔滨建筑大学 王焕定教授制作对各种具体单元，将Nw的具体形函数代入G矩阵的表达式，即可积

8、分得到具体单元的几何刚度矩阵。 基于上述结果，单元几何刚度矩阵为 对于分支点稳定问题，结构力学已经指出，在比例加载下，最终归结为一个特征值问题 解得特征值后，即可得到临界荷载，一般只关心最小临界荷载。2000.423哈尔滨建筑大学 王焕定教授制作3.1 物质描述大变形增量问题的T.L法 对弹塑性、粘性-蠕变和施工力学等问题，介质的反应和变形的历史有关。对随时间变化的荷载，需要将时间变量离散成序列，以求解各时刻的响应。为此，都需要用增量法来解决。 从t到t+t的增量期间进行物质描述求解时，一般可选两种参考位形：初始和t时刻的位形。前者称为全拉格朗日(T.L)表述，后者称为修正拉格朗日(U.L)表

9、述。 设从0到t时刻的全部反应、位形均已求得，现在的问题是，如何求t+t时刻的响应。2000.424哈尔滨建筑大学 王焕定教授制作未知 设t0、t和t+t的物理量分别用如下符号标记 对T.L法，介质位移是初始位形坐标的函数 坐标 ，密度 面积和体积 设有限元分析时单元形状描述为又设有限元分析时单元位移场为2000.425哈尔滨建筑大学 王焕定教授制作式中雅可比矩阵J为 象有限元(I)等参元分析一样，由于形函数一般是对自然坐标 定义的，因此有限元分析中的对坐标求导等，应象有限元(I)一样进行转换。To15To362000.426哈尔滨建筑大学 王焕定教授制作 在上述记号下，格林应变为 时刻t +

10、t和t的应变增量为式中2000.427哈尔滨建筑大学 王焕定教授制作 要强调指出的是，u在增量步内已知，因此 同大变形有限元，将张量转换为矩阵，则引入大变形所用算子记号，则有2000.428哈尔滨建筑大学 王焕定教授制作 由于u在增量步内已知，因此B和BL都是已知的。又若记则有 因为 ，因此综上可得2000.429哈尔滨建筑大学 王焕定教授制作 为进行有限元列式，还需讨论克希荷夫应力。设t和t+t时刻的应力分量分别为 基于上述分析，利用t+t时刻的虚位移原理虚功方程则象应变分析一样，可将 分成 。同样换为矩阵表示，则有2000.430哈尔滨建筑大学 王焕定教授制作 再次强调，t时刻及其前的量都

11、是已知的，因此变分为零。基于此将此结果代入虚功方程，可得单元刚度方程式中 和 是对初始位形定义的，t+t时刻的体积力和表面力，它们是已知的。式中2000.431哈尔滨建筑大学 王焕定教授制作t时刻应力引起的等效结点力矩阵t+t时刻荷载引起的等效结点力矩阵或 将其按集成规则集装后可得再引入如下记号将 和 的表达式代入，可得2000.432哈尔滨建筑大学 王焕定教授制作几何或非线性应变增量刚度矩阵或利用这些关系，非线性平衡方程可写为 为求解上述方程，尚需解决如下两方面问题2000.433哈尔滨建筑大学 王焕定教授制作 首先假设 然后将S和E的关系线性化。根据本构关系则有为使其线性化，设( t时刻的

12、材料性质矩阵)2000.434哈尔滨建筑大学 王焕定教授制作 在做了上述两方面处理后，可得由此出发，用非线性方程的相关解法，即可解决大变形非线性（材料非线性）问题。将其代回非线性平衡方程，可得 讲义上给出了T.L法的求解步骤，可供大家编程序参考。2000.435哈尔滨建筑大学 王焕定教授制作3.2 物质描述大变形增量问题的U.L法 因t+t的位移是用t时刻位形为基准度量的，因此 在t,t+t间隔内，以t时刻位形为参考位形，其增量位移为象T.L法一样，设单元和位移的描述为但需指出的是，式中形函数N是t时刻单元自然坐标的函数。在计算 等导数时，要先作坐标变换(Xi应换为xi)。To262000.4

13、36哈尔滨建筑大学 王焕定教授制作 类似地，用矩阵来表示则有 在时刻t和t+t的格林应变是以t时刻位形定义的，因而它们分别为式中算子矩阵象T.L法一样，但应将Xi换为xi。2000.437哈尔滨建筑大学 王焕定教授制作再次强调，式中算子符号象T.L法一样，但应将Xi换为xi。 基于上述说明，象T.L法一样可导得 关于应力的处理也和T.L法一样，对t 时位形定义的t和t+t时刻的克希荷夫应力分别为2000.438哈尔滨建筑大学 王焕定教授制作 象T.L法一样由虚位移原理虚功方程可导得 象T.L法一样推导，引入如下符号定义式中其中 分别为t时刻位形定义的单元体积、应力表面、体力和表面力。几何或非线

14、性应变增量刚度矩阵荷载引起的等效结点力矩阵2000.439哈尔滨建筑大学 王焕定教授制作则可得t时刻应力引起的等效结点力矩阵t+ t时刻的非线性平衡方程 象T.L法一样，为求解上述方程也需解决线性化问题。首先讨论S的计算。因为2000.440哈尔滨建筑大学 王焕定教授制作因此可改写为 由第四章已知Chap4 64式中Chap4 44Chap4 262000.441哈尔滨建筑大学 王焕定教授制作可得 由如下两式消去 并利用 ,且注意到Vij对称、ij反对称GoTo 49又由于 、 和 ，因此GoTo 502000.442哈尔滨建筑大学 王焕定教授制作 上式最后一项将使本构张量不对称，对金属类不可压缩介质，这一项可略去，也即在有限的克希荷夫应力和格林应变增量之间仍认为这就是U.L法的本构关系线性化。 与T.L法一样，除本构关系线性化外，还需解决几何方面的线性化。因为2000.443哈尔滨建筑大学 王焕定教授制作为对上述三项积分作几何方面的线性化，可设这一关系可有两种矩阵表达方式：其一是其另一方式是2000.444哈尔滨建筑大学 王焕定教授制作引入讲义上式(5. 102)算子矩阵 ，则有 基于上述讨论，积分 时，作如下处理式中 为由 组成的66矩阵， 矩阵