大变形问题的有限单元法

大变形问题的有限单元法第一页，共七十三页，2022年，8月28日

大变形时平衡方程和虚位移原理变形体初始和现时位形如图所示，以欧拉应力表述平衡时这是现时位形空间描述的平衡条件。在外荷为保守力系时第二页，共七十三页，2022年，8月28日经推证得上式乘以J，由于由复合函数求导数，平衡方程改写为平衡方程成为对任意j恒有第三页，共七十三页，2022年，8月28日证明以j=1为例同理可验证，对任意j恒有第四页，共七十三页，2022年，8月28日上述拉格朗日和克希荷夫应力表示的平衡条件都是以初始位形作参考的物质描述。利用克希荷夫应力和拉格朗日应力的关系，可将平衡条件改为第五页，共七十三页，2022年，8月28日如果考虑到变形梯度和位移梯度间的关系，对比小变形情况，可见大变形时变形对平衡的影响，是通过变形或位移梯度表现出来的。则克希荷夫应力表达的平衡方程可改为第六页，共七十三页，2022年，8月28日设现时位形微小虚位移在V内单值连续、在位移边界上为零。则外力总虚功为考虑到位移边界虚位移为零和应力边界条件空间描述的虚位移原理第七页，共七十三页，2022年，8月28日将平衡方程引入，考虑到虚位移微小，则利用格林公式，可得也即虚位移原理的虚功方程为若虚位移用虚速度、虚应变用虚应变率代替，则柯西应变虚功率方程第八页，共七十三页，2022年，8月28日为建立物质描述虚功方程，先讨论能量共轭关系。空间描述中又因为变形率张量是相对现时位形定义的柯西应变的速率单位体积变形功率因此克希荷夫应力和格林应变在能量上共轭。第九页，共七十三页，2022年，8月28日

由于变形率和欧拉应力张量是对称的，因此再利用欧拉和拉格朗日应力间关系，可得因此拉格朗日应力和初始位形位移梯度在能量上共轭。由此即可得到物质描述的虚功方程为第十页，共七十三页，2022年，8月28日由虚位移原理可得式中为单元结点力矩阵，为单元等效结点荷载矩阵Pd和PE分别为直接和等效结点荷载矩阵，R为综合等效结点荷载矩阵。上式也可写为按集成规则集装后可得第十一页，共七十三页，2022年，8月28日1.弹性大变形问题的有限元法

弹性大变形问题，需要考虑变形的非线性项和变形对平衡的影响。若以初始自然平衡状态作初始位形，则物质描述的格林应变为式中线性部分非线性部分第十二页，共七十三页，2022年，8月28日为便于计算机编程，将张量转换为矩阵：格林应变矩阵和张量的分量间有如下关系对应的克希荷夫应力矩阵和张量分量间关系为引入两个算子矩阵第十三页，共七十三页，2022年，8月28日式中再引入位移梯度向量的记号3阶单位矩阵第十四页，共七十三页，2022年，8月28日第十五页，共七十三页，2022年，8月28日在上述符号基础上，格林应变有位移表为则单元格林应变为其中线性应变矩阵B和线性分析一样设单元位移场为非线性部分“应变矩阵”为第十六页，共七十三页，2022年，8月28日式中G为如下9×3m的矩阵第十七页，共七十三页，2022年，8月28日式中AL为如下6×9的矩阵第十八页，共七十三页，2022年，8月28日由(AL)可见，格林应变-位移关系是非线性的。非线性部分因为为用虚位移原理建立单元刚度方程，还需有应变增量和位移增量间的关系。对线性部分第十九页，共七十三页，2022年，8月28日所以综上所述，格林应变增量为如果记，则。对弹性问题，在物质描述下本构关系(克希荷夫应力和格林应变)为第二十页，共七十三页，2022年，8月28日在上述基础上，由虚位移原理可得式中为单元结点力矩阵，为单元等效结点荷载矩阵Pd和PE分别为直接和等效结点荷载矩阵，R为综合等效结点荷载矩阵。上式也可写为按集成规则集装后可得第二十一页，共七十三页，2022年，8月28日将克希荷夫应力表达式代入，可得式中K(U)是非对称的，为对非线性弹性问题、和都是位移的函数。根据非线性方程切线刚度矩阵的定义，可得第二十二页，共七十三页，2022年，8月28日根据本构方程，则有又因式中所以第二十三页，共七十三页，2022年，8月28日又因,所以由此可得式中第二十四页，共七十三页，2022年，8月28日基于上述说明，可得如果引入如下记号初应力或几何刚度矩阵线弹性刚度矩阵大位移刚度矩阵第二十五页，共七十三页，2022年，8月28日则单元切线刚度矩阵为“结构”切线刚度矩阵为建立了切线刚度矩阵，用牛顿法等即可求解。以上的讨论没涉及具体单元，因此具有普遍性。具体单元分析时，因形函数一般是自然坐标的函数，故需作坐标变换后代入相关公式，从而建立具体单元的切线刚度矩阵。第二十六页，共七十三页，2022年，8月28日牛顿迭代法解大变形问题的具体步骤为：

第二十七页，共七十三页，2022年，8月28日牛顿迭代法解大变形问题的具体步骤为：

第二十八页，共七十三页，2022年，8月28日以上是全量形式的弹性大变形分析，为了保证收敛，拟用增量迭代法。极值点失稳问题弹塑性问题第二十九页，共七十三页，2022年，8月28日2.弹性分支点稳定问题有限元分析研究弹性稳定问题时材料应力应变关系是弹性的，也即

。式中

对线弹性问题由弹性常数表示，对非线性弹性问题由割线或切线系数表示。稳定问题可分为两大类：分支点（第一类）稳定问题；极值点（第二类）稳定问题。分支点稳定问题存在稳定平衡状态和不稳定平衡状态的分界点，在此临界状态下，系统可在两种不同的变形形式下保持平衡。极值点（第二类）稳定问题变形形式整个过程中是不变的，但当荷载达到临界值时，变形将迅速增加使系统丧失承载能力或破坏。只讨论第一类线性弹性稳定问题。第三十页，共七十三页，2022年，8月28日2.弹性分支点稳定问题有限元分析

对于分支点稳定问题时，可设外荷是比例加载的。由于第一类稳定问题主要是解决临界荷载计算，因此失稳状态相对稳定平衡状态的位移仍然是微小的，故可略去上节所推得的大位移矩阵项。但因为现时失稳位形和初始稳定平衡是两种不同的变形形式，因此应变要考虑失稳位移的非线性项。临界状态时稳定平衡和不稳定平衡荷载不变，前者是没有与失稳位移对应的荷载的，由结构力学已知，稳定分析的关键是确定几何刚度矩阵。因此集装后的整体切线刚度方程为式中

为整体线性刚度矩阵，

为整体几何刚度矩阵，

为不稳定平衡位置相对稳定平衡位置的整体位移。

第三十一页，共七十三页，2022年，8月28日2.弹性分支点稳定问题有限元分析对分支点稳定问题，关键是建立几何刚度矩阵。此时，以失稳前的平衡位置作初始位形，以失稳形态作现时位形。和前述弹性大变形不同的是：大位移矩阵可忽略。应变仅包含失稳位移的非线性项。分支点处相应失稳位移的综合荷载为零。失稳前变形是微小的。第三十二页，共七十三页，2022年，8月28日2.弹性分支点稳定问题有限元分析由于是比例加荷，在基准荷载下初始位形（稳定平衡位置）时的应力可由线性分析求得，因此上节几何刚度矩阵式子里的S

与荷载比例系数成正比。

在比例加载时

为荷载系数，

为单元基准荷载下的几何刚度矩阵。

由此可见，第一类稳定问题最终归结为求解一个特征值问题。解得特征值后，即可得到临界荷载，一般只关心最小临界荷载。第三十三页，共七十三页，2022年，8月28日设变形前单元长度为l，截面积为A，弹性模量为E。单元杆端位移矩阵为δe式中2.1桁架单元。单元位移为式中为失稳位移（出平面的位移）。基于此，单元格林应变为（针对失稳前定义）第三十四页，共七十三页，2022年，8月28日克希荷夫应力为因分支点稳定关心的临界荷载，临界荷载时平衡有两重性，故临界状态应力等于失稳前状态的应力，也即基于上述结果，单元几何刚度矩阵为，失稳前应力为第三十五页，共七十三页，2022年，8月28日设变形前单元长度为l，截面积和惯性矩为A、I，弹性模量为E。单元杆端位移矩阵为δe单元位移为式中2.2梁单元第三十六页，共七十三页，2022年，8月28日式中Ni和线性有限元梁单元一样。梁单元应变为其中第一项为线性有限元里的线性项，非线性的第二项为基于上述结果，单元几何刚度矩阵为象桁架单元说明一样，单元应力为第三十七页，共七十三页，2022年，8月28日3.1物质描述大变形增量问题的T.L法对弹塑性、粘性-蠕变和施工力学等问题，介质的反应和变形的历史有关。对随时间变化的荷载，需要将时间变量离散成序列，以求解各时刻的响应。为此，都需要用增量法来解决。从t到t+Δt的增量期间进行物质描述求解时，一般可选两种参考位形：初始和t时刻的位形。前者称为全拉格朗日(T.L)表述，后者称为修正拉格朗日(U.L)表述。设从0到t时刻的全部反应、位形均已求得，现在的问题是，如何求t+Δt时刻的响应。第三十八页，共七十三页，2022年，8月28日未知设t0、t和t+Δt的物理量分别用如下符号标记对T.L法，介质位移是初始位形坐标的函数坐标，密度面积和体积设有限元分析时单元形状描述为又设有限元分析时单元位移场为第三十九页，共七十三页，2022年，8月28日式中雅可比矩阵J为象线性有限元等参元分析一样，由于形函数一般是对自然坐标定义的，因此有限元分析中的对坐标求导等，应象线性有限元一样进行转换。第四十页，共七十三页，2022年，8月28日在上述记号下，格林应变为时刻t－t+Δt的应变增量为式中t+Δt时刻

t时刻第四十一页，共七十三页，2022年，8月28日由于，u(t时刻)在增量步内已知，因此同大变形有限元，将张量转换为矩阵，则引入大变形所用算子记号，则有第四十二页，共七十三页，2022年，8月28日由于u在增量步内已知，因此B和BL都是已知的。又若记则有因为，因此综上可得第四十三页，共七十三页，2022年，8月28日为进行有限元列式，还需讨论克希荷夫应力。设t和t+Δt时刻的应力分量分别为基于上述分析，利用t+Δt时刻的虚位移原理虚功方程则象应变分析一样，可将分成。同样换为矩阵表示，则有第四十四页，共七十三页，2022年，8月28日再次强调，t时刻及其前的量都是已知的，因此变分为零。基于此将此结果代入虚功方程，可得单元刚度方程式中和是对初始位形定义的，t+Δt时刻的体积力和表面力，它们是已知的。式中第四十五页，共七十三页，2022年，8月28日t时刻应力引起的等效结点力矩阵t+Δt时刻荷载引起的等效结点力矩阵或将其按集成规则集装后可得再引入如下记号将和的表达式代入，可得第四十六页，共七十三页，2022年，8月28日几何或非线性应变增量刚度矩阵或利用这些关系，非线性平衡方程可写为为求解上述方程，尚需解决如下两方面问题第四十七页，共七十三页，2022年，8月28日首先假设然后将ΔS和ΔE的关系线性化。根据本构关系则有为使其线性化，设(t时刻的材料性质矩阵)第四十八页，共七十三页，2022年，8月28日在做了上述两方面处理后，可得由此出发，用非线性方程的相关解法，即可解决大变形非线性（材料非线性）问题。将其代回非线性平衡方程，可得第四十九页，共七十三页，2022年，8月28日3.2物质描述大变形增量问题的U.L法因t+Δt的位移是用t时刻位形为基准度量的，因此在[t,t+Δt]间隔内，以t时刻位形为参考位形，其增量位移为象T.L法一样，设单元和位移的描述为但需指出的是，式中形函数N是t时刻单元自然坐标的函数。在计算

等导数时，要先作坐标变换(Xi应换为xi)。第五十页，共七十三页，2022年，8月28日类似地，用矩阵来表示则有在时刻t和t+Δt的格林应变是以t时刻位形定义的，因而它们分别为式中算子矩阵象T.L法一样，但应将Xi换为xi。第五十一页，共七十三页，2022年，8月28日式中算子符号象T.L法一样，但应将Xi换为xi。象T.L法一样可导得关于应力的处理也和T.L法一样，对t时位形定义的t和t+Δt时刻的克希荷夫应力分别为第五十二页，共七十三页，2022年，8月28日象T.L法一样由虚位移原理虚功方程可导得象T.L法一样推导，引入如下符号定义式中其中分别为t时刻位形定义的单元体积、表面应力、体力和表面力。几何或非线性应变增量刚度矩阵荷载引起的等效结点力矩阵第五十三页，共七十三页，2022年，8月28日则可得t时刻应力引起的等效结点力矩阵t+Δt时刻的非线性平衡方程象T.L法一样，为求解上述方程也需解决线性化问题。首先讨论ΔS的计算。因为第五十四页，共七十三页，2022年，8月28日因此可改写为由第四章已知式中第五十五页，共七十三页，2022年，8月28日可得由如下两式消去并利用,且注意到Vij对称、Ωij反对称又由于、和，因此第五十六页，共七十三页，2022年，8月28日上式最后一项将使本构张量不对称，对金属类不可压缩介质，这一项可略去，也即在有限的克希荷夫应力和格林应变增量之间仍认为这就是U.L法的本构关系线性化。与T.L法一样，除本构关系线性化外，还需解决几何方面的线性化。因为第五十七页，共七十三页，2022年，8月28日为对上述三项积分作几何方面的线性化这一关系可有两种矩阵表达方式：其一是把写成六维向量形式以B代替进行几何非线性的线性化第五十八页，共七十三页，2022年，8月28日

第五十九页，共七十三页，2022年，8月28日其另一方式是把写成九维向量形式

引入算子矩阵，则有第六十页，共七十三页，2022年，8月28日

第一个积分用

作第二、三个积分时

基于上述讨论，积分时，作如下处理式中为由组成的6×6矩阵，矩阵为第六十一页，共七十三页，2022年，8月28日结束则t+Δt时刻的非线