非线性有限元分析.docx
《非线性有限元分析.docx》由会员分享，可在线阅读，更多相关《非线性有限元分析.docx（20页珍藏版）》请在冰点文库上搜索。
非线性有限元分析
非线性有限元分析
1概述
在科学技术领域内，对于许多力学问题和物理问题，人们已经得到了它们所应遵循的基本方程（常微分方程或偏微分方程）和相应的定解条件（边界条件）。
但能够用解析方法求出精确解的只是少数方程性质比较简单，并且几何形状相当规则的问题。
对于大多数工程实际问题，由于方程的某些特征的非线性性质，或由于求解区域的几何形状比较复杂，则不能得到解析的答案。
这类问题的解决通常有两种途径。
一是引入简化假设，将方程和几何边界简化为能够处理的情况，从而得到问题在简化状态下的解答。
但是这种方法只是在有限的情况下是可行的，因为过多的简化可能导致误差很大甚至是错误的解答。
因此人们多年来一直在致力于寻找和发展另一种求解途径和方法——数值解法。
特别是五十多年来，随着电子计算机的飞速发展和广泛应用，数值分析方法已成为求解科学技术问题的主要工具。
已经发展的数值分析方法可以分为两大类。
一类以有限差分法为代表，主要特点是直接求解基本方程和相应定解条件的近似解。
其具体解法是将求解区域划分为网格，然后在网格的结点上用差分方程来近似微分方程，当采用较多结点时，近似解的精度可以得到改善。
但是当用于求解几何形状复杂的问题时，有限差分法的精度将降低，甚至发生困难。
另一类数值分析方法是首先建立和原问题基本方程及相应定解条件相等效的积分提法，然后再建立近似解法并求解。
如果原问题的方程具有某些特定的性质，则它的等效积分提法可以归结为某个泛函的变分，相应的近似解法实际上就是求解泛函的驻值问题。
诸如里兹法，配点法，最小二乘法，伽辽金法，力矩法等都属于这一类方法。
但此类方法也只能局限于几何形状规则的问题，原因在于它们都是在整个求解区域上假设近似函数，因此，对于几何形状复杂的问题，不可能建立合乎要求的近似函数。
1960年，R.W.CLOUGH发表了有限单元法的第一篇文献“TheFiniteElementMethodinPlaneStressAnalysis”，这同时也标志着有限单元法（FEM）的问世。
有限单元法的基本思想是将连续的求解区域离散为一组有限个，且按一定方式相互联接在一起的单元的组合体。
由于单元能按不同的联结方式进行组合，且单元本身又可以有不同形状，因此可以模型化几何形状复杂的求解域。
并且可以利用在每一个单元内假设的近似函数来分片地表示全求解域上待求的未知场函数，从而使一个连续的无限自由度问题变成离散的有限自由度问题。
现已证明，有限单元法是基于变分原理的里兹法的另一种形式，从而使里兹法分析的所有理论基础都适用于有限单元法，确认了有限单元法是处理连续介质问题的一种普遍方法。
利用变分原理建立有限元方程和经典里兹法的主要区别是有限单元法假设的近似函数不是在全求解域而是在单元上规定的，而且事先不要求满足任何边界条件，因此可以用来处理很复杂的连续介质问题。
在短短四十余年的时间里，有限单元的分析方法已经迅速地发展为适合于使用各种类型计算机解决复杂工程问题的一种相当普及的方法。
如今，有限元广泛地应用于各个学科门类，已经成为工程师和科研人员用于解决实际工程问题，进行科学研究不可或缺的有力工具。
有限单元法的应用范围已由弹性力学平面问题扩展到空间问题，板壳问题，由静力平衡问题扩展到稳定问题，动力问题和波动问题。
分析的对象从弹性材料扩展到塑性，粘弹性，粘塑性和复合材料等，从固体力学扩展到流体力学，传热学等连续介质力学领域。
在工程分析中的作用已从分析和校核扩展到优化设计并和计算机辅助设计技术相结合。
各种各样商业化的大型通用有限元软件层出不穷，不断推陈出新。
可以预见，随着现代力学，计算数学，计算机技术等学科的发展，有限单元法作为一个具有巩固理论基础和广泛应用范围的数值分析工具，必将得到进一步的完善和发展。
2非线性问题的类型和求解特点
2．1非线性问题的类型
2.1.1线性分析的含义
在有限元分析中的线性假设包含下列含义：
即结点位移为无限小量，材料为线弹性，加载时边界条件的性质保持不变。
于是，静力平衡方程可以表示为：
（2.1）
其中，
为刚度矩阵，
为荷载矢量。
由于
和
的元素为常数，故位移响应
是荷载矢量
的线性函数。
也就是说，如果
变为
，则
变为
，其中，
为常数。
这就是所谓的线性有限元分析。
如果上述假设中的任何一条不能得到满足，那么就属于非线性有限元分析。
2.1.2非线性分析的必要性
结构力学问题，从本质上讲都是非线性的，线性假设只是实际工程问题的一种简化。
当然，任何实际工程问题的求解都避免不了适当地简化，简化是否合理主要应根据求解效果和实际经验来判断。
对于目前工程实际中的很多问题，如地震作用下结构的弹塑性动力响应，高层建筑抗风，大跨度网壳结构动力稳定性，索膜结构找形荷载与裁减分析，大型桥梁风致振动等问题的研究，仅仅假设为线性问题是很不够的，常常需要进一步考虑为非线性问题。
因此，对各种工程结构的非线性分析就是必不可少且日趋重要了。
对于结构力学的非线性问题来说，有限单元法是最为有效的数值分析方法。
2.1.3非线性问题的类型
通常，把非线性问题分为两大类，即分为几何非线性和材料非线性。
但从建立基本方程和程序设计的方便出发，又可分为三种类型：
1．材料非线性：
非线性效应仅由应力应变关系的非线性引起，位移分量仍假设为无限小量，故仍可采用工程应力和工程应变来描述，即仅材料为非线性。
非线性的应力应变关系是结构非线性的常见原因，许多因素都可以影响材料的应力应变性质，包括加载历史（如在弹塑性响应状况下），环境状况（如温度），加载的时间总量（如在蠕变响应状况下）等。
2．几何非线性：
如果结构经受大变形，则变化了的几何形状可能会引起结构的非线性响应，这又可以分为两种情形：
第一种情形，大位移小应变。
只是物体经历了大的刚体平动和转动，固连于物体坐标系中的应变分量仍假设为无限小。
此时的应力应变关系则根据实际材料和实际问题可以是线性的也可以是非线性的。
第二种情形，大位移大应变。
也即最一般的的情况，此时结构的平动位移，转动位移和应变都不再是无限小量，本构关系也是非线性的。
3．状态非线性：
除以上两种非线性问题之外，还有一种非线性问题，即由于系统刚度和边界条件的性质随物体的运动发生变化所引起的非线性响应。
例如，一根只能拉伸的钢索可能是松散的,也可能是绷紧的；轴承套可能是接触的,也可能是不接触的；冻土可能是冻结的,也可能是融化的。
这些系统的刚度和边界条件由于系统状态的改变在不同的值之间突然变化。
状态改变也许和载荷直接有关，也可能由某种外部原因引起。
最为典型的就是接触问题，接触是状态非线性类型中一个特殊而重要的子集。
通常情况下，状态非线性问题可以在上述材料非线性和几何非线性类型中的每一种同时出现，从而使得问题的分析变得更为复杂。
2．2非线性问题的求解特点
2.2.1非线性分析的基本问题
非线性分析的基本问题是求出在当前荷载作用下的平衡状态。
如果作用的荷载被描述成时间的函数，则物体有限元离散系统的平衡方程可以表示为：
（2.2）
其中，矢量
由t时刻外荷载的结点力分量所构成，而矢量
则表示t时刻的单元应力所引起的结点力分量。
平衡方程（2.2）应针对t时刻的几何位形建立，并应计入所有的非线性效应。
如果是动力分析，矢量
中还应当包括惯性力和阻尼力。
在求解非线性问题时，（2.2）式应在全部加载历史中成立。
变量t的引入并不意味着一定是动力问题。
在静力分析中，t不具有真实“时间”的含义，它的不同取值只是表示相应于不同位形的不同的荷载水平。
但是，在动力分析或具有时间效应的静力分析中，变量t就有了它本来的“时间”的含义。
2.2.2非线性方程组的增量逐步解法
对于许多工程结构，我们所关心的常常是在特定的荷载水平下，或相应的时间物体中的应力和变形。
实际问题根据其解法可以分为两大类型。
第一类问题无需计算中间变形过程，可直接求解在给定荷载下的平衡位形。
但是，如果问题的几何性质或材料性质与路径相关或与时间相关，即该问题依赖于变形历史，则中间变形过程的计算是不可缺少的，这就是第二类问题。
从本质上来说，非线性问题是第二类问题。
此时，往往采用增量分析的方法。
增量逐步解法的基本思想是：
假定t时刻的解为已知，要求t＋Δt时刻的解，其中，Δt是适当选择的时间增量。
在t＋Δt时刻，式（2.2）写成为：
（2.3）
这里，左上标表示为t＋Δt时刻的量。
由于t时刻的解为已知，因此，可以写为：
（2.4）
式中，
表示t到t＋Δt时间间隔内，由于单元内应力增量所引起的结点力增量矢量。
这一矢量可以近似表示为：
（2.5）
式中，
为相应于t时刻材料和几何条件的切线刚度矩阵。
为Δt时间间隔中的结点位移增量，现在它还是未知的。
将式（2.4）和（2.5）代入式（2.3）中，得到：
（2.6）
上式中只有位移增量
为未知，一旦解出，即可算得t＋Δt时刻的位移：
（2.7）
根据
，就容易算出t＋Δt时刻的应力及
，
，于是马上可以着手下一步的计算。
但要指出的是，式（2.5）是一个近似表达式，因此t＋Δt时刻的解也是近似的，如果急于求成的作下去，最终结果可能出现不可忽视的重大误差以致于达到荒谬的地步。
解决这一困难的办法是以花费计算时间为代价，即在t到t＋Δt时步中进行足够次数的迭代，以保证最终的解获得足够的精度。
2.2.3引入修正Newton－Raphson迭代格式的增量逐步解法
现在更多采用的方法是在每一个荷载增量步中，使用Newton－Raphson迭代法或修正的Newton－Raphson迭代法。
由于后者不需要每次迭代时都计算切线刚度矩阵，因此在实际中具有更广泛的应用。
现对该方法做简单的介绍。
在t时刻到t＋Δt时刻的时步中，修正Newton－Raphson法的迭代公式可以表示为：
（2.8）
（2.9）
其中，i表示迭代步数，依次取1，2，3，…，其迭代所用的初始值正是t时刻的解，即：
（2.10）
式（2.8）的右端项：
称为第i步迭代前的不平衡荷载。
在迭代过程中，
随i的增加而逐步接近
。
因此，我们可事先对不平衡荷载的模给定一个精度指标，每次迭代后检查不平衡荷载是否小于该指标。
若满足精度，则在求出
之后转入下一时步的计算，否则继续迭代，直到满足精度要求为止。
3材料非线性问题的有限单元法
3．1材料非线性问题概述
在所有的非线性分析问题中，材料非线性问题的处理相对简单，不需要重新列出整个问题的表达格式，只要将材料本构关系线性化，就可将线性问题的表达格式推广用于非线性分析。
一般来说，通过试探和迭代的过程求解一系列线性问题，如果在最后阶段，材料的状态参数被调整得满足材料的非线性本构关系，则最终可得到问题的解答。
材料非线性问题可以分为两种类型。
一类是不依赖于时间的弹塑性问题，其特点是当荷载作用以后，材料变形立即发生，并且不再随时间而变化。
另一类是依赖于时间的粘（弹，塑）性问题，其特点是荷载作用以后，材料不仅立即发生变形，而且变形随时间继续变化。
在荷载保持不变的条件下，由于材料粘性而继续增长的变形称之为蠕变；另一方面，在变形保持不变的条件下，由于材料粘性而使应力衰减称之为松弛。
显然，后一类材料非线性问题在求解时更为困难一些。
3．2材料非线性本构关系
限于篇幅，本文仅讨论最为常见的弹塑性非线性本构关系。
弹塑性材料进入塑性的特征是当荷载卸去以后存在不可恢复的永久变形，因而在涉及卸载的情况下，应力应变之间不再存在惟一的对应关系，这是区别于非线弹性材料的基本属性。
以材料的单向受力情况为例，只是在加载时应力应变呈现非线性关系，还不足以判定材料是非线性弹性还是弹塑性。
但是一经卸载立即发生两者的区别，非线性弹性材料将沿原路径返回，而弹塑性材料将依据不同的加载历史卸载后产生不同的永久变形。
任何一种弹塑性材料都应当满足塑性力学的四条基本准则，这里对此作简单的介绍：
1．初始屈服条件：
规定了材料开始塑性变形的应力状态。
在有限元分析中，通常采用V.Mises准则。
2．流动准则：
规定塑性应变增量的分量和应力分量以及应力增量分量之间的关系。
3．硬化准则：
规定材料进入塑性变形后的后继屈服函数。
对于理想弹塑性材料，因无硬化效应，后继屈服函数和初始屈服函数一致；对于硬化材料，通常又有各向同性硬化准则，随动硬化准则和混合硬化准则三种不同的准则。
4．加载，卸载准则：
用以判别从一塑性状态出发是继续塑性加载还是弹性卸载，这是计算中判定是否继续塑性变形以及决定是采用弹塑性本构关系还是弹性本构关系所必须的。
各种类型的弹塑性材料可以从对各自的后继屈服函数进行微分出发，进而推导出各自相应的应力应变的增量关系，这里不一一列举。
需要进一步说明的是，对于处于高温条件下工作的结构，必须考虑温度对本构关系的影响。
比如随着温度的升高，屈曲极限有所降低，材料硬化特性也有所减少，并逐渐接近理想塑性材料，同时材料常数E，μ，α等也随温度变化而有所变化。
至于长期工作在高温条件下的结构还必须考虑蠕变的效应。
3．3材料非线性问题的有限元表达格式
对于弹塑性材料，由于材料和结构的弹塑性行为与加载以及变形的历史有关。
因此，在进行结构的弹塑性分析时，通常将荷载分成若干个增量，然后对于每一荷载增量，将弹塑性方程线性化，从而使弹塑性分析这一非线性问题分解为一系列的线性问题。
按照这种思想，首先建立增量形式的荷载条件和位移条件，进而建立增量形式的虚位移原理，即增量形式的最小势能原理，最终即可得到基于增量形式的有限元表达格式。
系统平衡方程形式同前（3.5）式，其中切线刚度矩阵
在这里是系统的弹塑性刚度矩阵。
弹塑性增量有限元分析在将加载过程划分为若干增量步以后，对于每一个增量步应包含下列三个算法步骤：
1．线性化弹塑性本构关系，并形成增量有限元方程。
2．求解有限元方程。
注意在求解过程中每个增量步或每次迭代时弹塑性刚度矩阵都可能发生局部的变化。
3．积分本构方程，决定新的应力状态，检查平衡条件，并决定是否进行新的迭代。
上述每一步骤的算法方案和数值方法，以及荷载增量步长的选择都关系到整个求解过程的稳定性，精度和效率。
这里尤其需要注意的是非线性方程组求解方案的选择。
通常可以采用以下几种求解方案：
无迭代的增量解法，具有变刚度迭代（N-R迭代）的增量解法和具有常刚度迭代（mN-R迭代）的增量解法。
变刚度迭代具有良好的收敛性，允许采用较大的时间步长，但每次迭代都要重新形成和分解新的刚度矩阵。
而采用常刚度迭代可以节省上述计算费用，缺点是收敛速度较慢，特别在接近荷载的极限状况时，因此经常需要同时采用加速迭代的措施。
具体采用何种求解方案，应根据具体问题的特点，综合考虑精度和效率两方面因素。
对于除弹塑性以外的材料非线性问题，例如热弹塑性—蠕变问题，粘弹塑性问题等，由于同时涉及独立于时间和依赖于时间的两类非弹性变形以及本构方程的高度非线性，无论是其本构方程的建立和它的积分方法，还是非线性方程组的求解方法都远比通常的弹塑性分析困难得多。
但还是有很多共性的方面，这里不再展开详述。
4几何非线性问题的有限单元法
4．1几何非线性问题概述
在某一固体力学问题中，如果假定物体所发生的位移远小于物体自身的几何尺度，应变远小于1，那么此问题就称作满足“小变形假定”。
在此前提下，建立物体或微元体的平衡条件时可以不考虑物体的位置和形状（简称位形）的变化。
因此分析中不必区分变形前和变形后的位形。
同时在加载和变形过程中的应变可用一阶无穷小的线性应变进行度量。
但是在实际中，我们往往会遇到很多不符合小变形假定的问题，例如板壳等薄壁结构的屈曲问题。
此时必须考虑变形对平衡的影响，即平衡条件应建立在变形后的位形上，同时应变表达式也应包括位移的二次项。
这样一来，平衡方程和几何关系都将是非线性的。
这种由于大平动和大转动引起的非线性问题称为几何非线性问题。
几何非线性问题还有另外一种类型，例如金属的成型，橡皮型材料受荷载作用，都可能会出现很大的应变，这时除了采用非线性的平衡方程和几何关系以外，还需要引入相应的应力应变关系，尽管对于后一问题材料通常还处于弹性状态。
当然大多数大应变问题是和材料的非弹性性质联系在一起的。
这类几何非线性问题即通常所说的大平动，大转动，大应变问题。
4．2几何非线性问题的有限元表达格式
早期几何非线性有限元分析基本上是线性分析的扩展，针对各个具体问题分别进行分析。
而近年来，基于非线性连续介质力学原理的有限元分析取得很大发展，得到了统一的一般非线性分析的表达格式。
基于非线性连续介质力学，首先应当对大变形情况下的应变和应力进行度量。
这是因为在非线性问题中，由于存在的大位移，大应变而导致有限变形，使得原来传统的小变形下的Cauchy方程不再适用。
此时，根据连续体在不同的位形下坐标的变换，对变形前后物体上某一线段变形的度量可以采用两种不同的应变度量方式。
即用变形前坐标表示的Green应变张量和用变形后坐标表示的Almansi应变张量。
在大变形问题中，是用从变形后的物体内截取出的微元体来建立平衡方程和与之等效的虚功原理的。
因此，在从变形后物体内截取出的微元体上面定义的应力张量称为Euler应力张量。
如果用于变形前的位形，可以具体定义另外两种应力张量：
Lagrange应力张量和Kirchhoff应力张量。
此外，在连续介质力学中还定义了一种其分量不随材料刚体转动而变化的速率型应力张量，Jaumann应力速率张量。
在涉及几何非线性问题的有限单元法中，通常都采用增量分析的方法。
为了得到方程的解答，所有的变量都应参考某一已经求得的平衡位形。
在实际分析中，通常有以下两种选择：
1．全Lagrange格式（TotalLagrangeFormulation，简称T.L.格式），这种格式中所有变量以时间0的位形作为参考位形。
2．更新Lagrange格式（UpdatedLagrangeFormulation，简称U.L.格式），这种格式中所有变量以时间t的位形作为参考位形。
因为在求解过程中参考位形是不断改变的，所以称之为更新的Lagrange格式。
由以上两种格式导出的求解方程在理论上是等效的，如若采用数学上相一致的本构关系，它们将产生相同的结果。
但在求解的有限元矩阵方程本身和求解步骤上仍有一定的差别。
在通用的有限元程序中，通常同时包括这两种格式，使用时可以根据所分析问题及材料本构关系的具体特点和形式选择最有效的格式。
为进一步说明非线性分析的特点，下表列出按非线性问题的不同分类所采用的不同描述方法和应力应变。
表1非线性问题分类
分析类型
特点
描述方法
应力和应变
仅材料
非线性
平动位移，转动位移和应变无限小，
应力应变关系是非线性的。
仅材料非线性
（M.N.O.）
工程应力
工程应变
大平动
大转动
小应变
线元的平动位移和转动位移充分大，
但线元的伸长和线元之间的角度改变无限小，
应力应变关系是线性的或非线性的。
全Lagrange描述
（T.L.）
Kirchhoff应力
Green应变
更新Lagrange描述
（U.L.）
Cauchy应力
Almansi应变
大平动
大转动
大应变
线元的伸长和线元之间的角度改变充分大，
线元的平动位移和转动位移也可以充分大，
应力应变关系是线性的或非线性的。
全Lagrange描述
（T.L.）
Kirchhoff应力
Green应变
更新Lagrange-Jaumann
描述
（U.L.J.）
Jaumann应力率
Almansi应变率
4．3几何非线性问题有限元方程的求解
对于几何非线性有限元的求解，一般采用等参元对求解域进行离散。
两种表达格式T.L.和U.L.都可应用，关键在于对求解方程的线性化处理。
因为无论是T.L.格式还是U.L.格式，都是基于线性化处理后的虚位移原理建立的有限元矩阵方程，该矩阵方程仅是对于每一时间步长所应求解的非线性方程的近似。
由于系统的非线性性质，线性化处理带来的误差将可能导致解的漂移或不稳定。
因此，仍需采用基于Newton-Raphson迭代格式或修正Newton-Raphson迭代格式的增量逐步解法求解方程组。
在实际分析中，两种格式用于求解的时间一般情况下相差不多，究竟选择哪种格式通常取决于所采用的本构关系的具体形式。
也就应当在求解之初便首先区分是大应变还是小应变，选择格式已在表1中列出，此处不再详述。
和材料非线性问题相比，几何非线性问题有着更为复杂多样的荷载—位移路径，如在荷载控制下的疾速通过和位移控制下的疾速通过。
因此，荷载增量步长的自动选择就显得格外重要。
近些年来，广泛应用的一类荷载增量步长的自动选择方法是“广义弧长法”。
在广义弧长法中，用于调节荷载增量和位移增量在弧长ΔL中作用的比例因子α对弧长法的总体性能有很大的影响。
一般采用的比例因子有：
α＝1的球面弧长法，α＝0的柱面弧长法，α＝Sp的椭圆弧长法。
对于不同的结构和荷载情况，很难说以上α不同取值的三种情况中哪一种具有绝对的优势，但是α＝0的柱面弧长法具有较好的普遍适应性。
5杆索非线性有限元理论
在结构非线性有限元分析中，最为重要也最为基本的是建立精度适合的各种有限单元列式，并在基于某些假定的基础上推导出其单元刚度矩阵和有限元求解方程。
下由于目前的研究和应用中已经出现了相当多的上述单元的理论和模型，限于篇幅并基于应用角度，每种单元仅选择一种最为常用，精度也较高的单元加以介绍。
5．1非线性有限杆单元理论
5.1.1基本假定
1.杆单元只能承受轴向力；
2.杆单元的应力应变关系符合虎克定律；
3.杆单元位移变形为大位移小应变。
5.1.2刚度矩阵及有限元方程
假定单元位移函数线性插值：
（5.1）
（5.2）
在局部坐标系ξ中，应力应变关系为：
在局部和整体坐标系关系中，转轴时应力增量和应变增量的变换矩阵T为：
式中：
l,m,n是方向余弦。
由此可得：
（5.3）
式中：
，
表示局部坐标系下单元的应变和应力；
，
表示整体坐标系下单元的应力和应变。
应变增量和位移增量的关系可用[B]矩阵表示，则有限元矩阵可表示为：
（5.4）
式中：
SE，σE分别为Kirchhoff应力矩阵和向量。
将积分式展开，得到线性刚度矩阵，非线性刚度矩阵和内力项的矩阵表达式：
（5.5）
（5.6）
（5.7）
按虚位移原理的矩阵列式为：
（5.8）
上式即为有限元基本方程。
5．2非线性有限索单元理论
索结构在大跨结构中已得到广泛的应用。
随着连续长索的不断应用，对于索力学模型的精度要求也越来越高。
初期的研究以解析法为基础，对较为简单理想的外荷载和边界条件作了分析。
随着计算技术的提高，提出并采用了考虑大变形的各种离散模型，主要有：
两节点直线杆单元模型，以等效弹性模量来考虑垂度影响；两节点抛物线索单元模型，以及为了提高分析精度采用内插节点的多节点索单元（三节点，四节点，五节点索单元）模型和采用B样条基构建的索单元模型。
下面简要介绍悬链线索单元模型。
5.2.1基本假定
1.索为理想柔索，不受压且无弯曲刚度；
2.满足大变形，小应变要求；
3.索中外荷载沿索长均匀分布。
5.2.2刚度矩阵及有限元方程
作几何非线性分析时，索单元的切线刚度可按下述方法计算。
如图1中为一个索单元，其中i点的位移是⊿1，⊿2，⊿3；j点的位移是⊿4，⊿5，⊿6；节点力由原来的F01，F02，F03，F04，F05，F06增加到F1，F2，F3，F4，F5，F6。
此时，节点力及节点位移间的平衡方程式如下：
（假设
）
（5
升级会员 