第10章(非线性有限元)

1、公式号、图号等第十章 非线性动力有限元法当机械结构受到较大的外载荷，或受到持续时间较短的冲击载荷作用时，结构会产生过大的变形, 以至于必须考虑结构几何大变形对结构整体刚度及固有频率的影响，即所谓的几何非线性影响。另外, 对于多数非线性动力学问题，还需要考虑材料非线性、接触非线性等方面的影响。 非线性动力学分析求解的基本方程有如下形式 (4.141)式中，为粘性效应项，考虑阻尼、粘塑、粘弹等效应。P为外部激励。对于考虑各种非线性效应的动力学问题求解，需要对动力学方程进行直接时间积分。即非线性动力有限元分析具有如下特点：（1）问题分析过程需要考虑时间积分效应，不必做模态分析，不必提取固有频率；（2

2、）采用直接积分方法求解非线性动力学方程，需要对时间作积分计算，因此计算量远远大于线性模态动力学方法；（3）非线性动力学分析中可以施加不同类型的载荷，包括结点力、非零位移、单元载荷；（4）在每个时间步上，进行质量、阻尼、及刚度的集成，采用完整矩阵，不涉及质量矩阵的近似；（5）可以同时考虑几何、材料和接触等多种非线性效应。非线性动力有限元分析程序常采用隐式Hilber-Hughes-Taylor法进行时间积分运算。这种方法适于模拟非线性结构的动态问题，对于冲击、地震等激发的结构动态响应以及一些由于塑性或粘性阻尼造成的能量耗散，隐式算法特别有效。隐式积分方法需要对刚度矩阵求逆计算，并通过多次迭代求解

3、增量步平衡方程。隐式Hilber-Hughes-Taylor时间积分算法为无条件稳定，对时间步长没有特别的限制。采用子空间法也可以对动力学平衡方程作时间积分运算。子空间法是提取模态分析得到的各阶特征模态，并采用与线性模态动力学分析方法相近的分析方式进行求解。对于带有微小非线性效应的问题，如材料小范围进行入屈服、结点转角不大的情况，子空间法效率比进接积分法要高。此外，非线性动力有限元分析还可以采用显式动态算法，如中心差分法。显式时间积分算法为有条件稳定，其临界稳定时间步长限制了时间步长的大小，与有限元模型最小单元尺寸、材料应力波速等有关。显式时间积分法适于模拟高速冲击、接触等问题。上述方法的选择

4、需要综合考虑计算量、分析问题的规模、单元限制等多方面因素，需要丰富的有限元模拟的理论、经验和实践知识。以下以几何非线性问题和材料非线性问题为例介绍非线性有限元法，其中粘弹粘塑性非线性材料问题的分析是典型的非线性动力有限元的求解思想。9.1 几何非线性问题的有限元法几何非线性问题一般是指物体经历大的刚体位移和转动，但固连于物体坐标系中的应变分量仍假设为小量, 即大位移小应变情况。4.6.1.1 几何非线性问题的牛顿迭代法 由数值分析技术可知，求解非线性方程组的数值方法的常规方法是Newton-Raphson法，即牛顿迭代法，这是一种近似线性化迭代求解方法。对于非线性方程，具有一阶导数，在点作一阶

5、泰勒级数展开，它在点的线性近似为 (4.142)因此，非线性方程在附近似为线性方程： (4.143)当时，由上式求得步的修正项 (4.144)Newton-Raphson方法的迭代公式为 (4.145)在几何非线性有限元法中，结构的刚度矩阵与其几何位置有关，平衡方程由变形后的位形描述，因此，结构的刚度矩阵是几何变形的函数。设变形为, 结构的平衡方程式 (4.146)为一个非线性方程组。记非线性方程 (4.147)用Newton-Raphson方法求的根时，迭代公式分别为 (4.148)其中, 满足下式 (4.149)式中, 称为切线刚度矩阵，表达式为 (4.150) 在每一个迭代步中，通过求解

6、切线刚度矩阵，进而用进行迭代求解，称为Newton-Raphson方法，又称切线刚度法。牛顿法的收敛性是好的。但是某些非线性问题中，使用牛顿法迭代时，若出现奇异或病态，则对的求逆出现困难。关于这一点也可以采用其它修正办法，如引入阻尼因子。对于已经建立的有限元方程，设表示内为和外力矢量的总和，有 (4.151)式中, R为载荷列阵；为虚位移；为虚应变用应变的增量形式代入上式，消去项，可以得到非线性问题的一般平衡方程式为 (4.152)该式不论位移或应变的大小与否均成立。在有限变形中，应变和位移之间的关系是非线性的，即B矩阵是的非线性函数。但是，近似地可将进行如下分解： (4.153)式中, 为线

7、性应变分析的部分； 为由非线性变形引起的，与有关。假定应力应变关系为线弹性，于是有 (4.154)式中 为材料的弹性矩阵； 为初应变列阵；为初应力列阵对于式(4.152)的非线性平衡方程式，可用Newton-Raphson方法进行迭代求解。对该式微分，有 (4.155)不考虑初应变和初应力的影响，得并且这样可得 (4.156)这里 (4.157)式中为通常的小位移的线性刚度矩阵。矩阵则是由于大位移引起，它可以写成 (4.158)式(4.156)又可记成： (4.159)式中 (4.160)式中，是关于应力水平的对称矩阵，称之为初应力矩阵或几何刚度矩阵。 因此，用Newton-Raphson方法

8、迭代求解几何非线性问题的步骤为： (1) 用线弹性解作为，即一次近似； (2) 通过定义求出，求出； (3) 确定切线刚度矩阵；(4) ， ； (5) 重复上述迭代步骤，直至足够小。在这里，没有考虑载荷R可能由于变形而发生的变化，即在这里假设了载荷不因变形而改变其大小和方向，否则是非保守力作用下的大变形问题，在此不做讨论。4.6.1.2 典型单元的切线刚度矩阵求解具体的几何非线性问题时，必须计算单元的切线刚度矩阵。对于一般空间问题，无论位移和应变大小，都可以利用应变的基本定义写出位移和应变的关系式。用变形前的坐标做为自变量，可以用位移定义如下大变形问题的应变分量表达式 (4.161) 对于微小

9、位移情况，可以略去二次以上的偏导数项，得到小变形时的应变公式。在有限变形中，假设应变仍为小量。应变和位移之间的关系为： (4.162)式中为线性应变部分。对于非线性部分，可以写成： (4.163)式中 (4.164)式中C为矩阵。根据的定义，可以将表示成任意一点位移的函数，引入形函数N后，可以得到 (4.165)对于(4.163)式进行微分，得 (4.166)因此， (4.167)B矩阵为 (4.168)这样得到 (4.169) 另外，有 (4.170)利用矩阵C和列阵的性质，得到 (4.171)式中I为三阶单位矩阵，M是的六个应力分量组成的矩阵。因此几何刚度矩阵为 (4.172)故此，非线性

10、三维单元的切线刚度矩阵为 (4.173)作为特例，可以直接写出三角形单元的上述有关表达式。由三角形单元的位移模式 (4.174)其中，式中的等由结点坐标确定，为三角形单元的面积。 根据式(4.164)把式(4.174)代入上式得 (4.175) 由(4.165)式可以知道 (4.176)根据定义，由式(4.163)确定的平面问题的C矩阵为 (4.177)这样可以得到C的显式为(4.178)故 (4.179)而由线性问题给出，即 (4.180) 至此，、和G都是常数矩阵，只与单元结点坐标和结点位移有关。因此，线性刚度矩阵为： (4.181)式中为单元厚度。初始刚度矩阵为： (4.182)几何刚度

11、矩阵为： (4.183)式中 (4.184)因此，对于几何非线性问题，平面三角形单元的切线刚度矩阵可以由式(4.173)求出。9.2 材料非线性问题的有限元法材料非线性是指材料的本构方程是非线性的。一般主要分为两类: 一类是非线性弹性问题，如橡胶、塑料、岩石等，在加载时应力应变关系等性质呈现非线性的物理现象，卸载时可逆。另一类是指材料的弹塑性问题。材料超过屈服极限后呈现出非线性。在机械结构分析中，常见的本构方程主要有线弹性和非线性弹性模型，其特点是应力仅应变的函数，加卸载规律相同。公式如下： (4.185)其中，对于线弹性材料为常数，对于非线性弹性材料，是的函数。此外，还有超弹性模型、次弹性模

12、型、弹塑性模型等。对于弹塑性模型，可以认为弹塑材料发生塑性变形时，其总应变可以分解为两部分： (4.186)即总应变为弹性应变和塑性应变之和。加载时遵循一定规律，如Prandtl-Reuss方程，而卸载时为弹性。对于应力足够大时的金属、土壤、岩石等材料，都有此类特征。在非线性动力有限元分析时，具有黏性特性的模型十分重要。对于黏弹性模型，一般包括松驰（指突加应变作用下应力逐渐减少）和蠕变（指突加应力作用下应变会逐渐增加）。典型的松驰模型是如下Maxwell模型： (4.187)蠕变模型如Voigt-Kelvin模型为： (4.188)而弹/粘塑性模型是指材料的塑性变形与时间有关，本构方程中出现非

13、齐次的时间微分。下面以弹/粘塑性材料的变形分析为例，说明非线性动力有限元的基本步骤。即弹/粘塑性材料的变形过程中要考虑时间效应，材料开始屈服后，塑性流动、应力和应变均与时间有关。4.6.2.1 弹/粘塑性问题的基本表达式 假设总应变分离成弹性应变和粘塑性应变，即 (4.189)其中表示对时间的求导。总应力率取决于弹性应变率，有 (4.190)式中，D是弹性矩阵。粘塑性性质为 (4.191)式中，是单向屈服应力，是硬化参数的函数。设粘塑性应变率仅取决于当前的应力，如下式所示： (4.192)式中，是塑性势，是流动参数，项对于是正单调增量函数，即 (4.193)在这里只讨论的情况。函数的两种常用形

14、式为 (4.194)或 (4.195)式中，为常数。4.6.2.2 粘塑性应变增量和应力增量 对于式(4.192)所表示的应变定律，定义时间间隔出现的应变增量为，有 (4.196)其中为常数， 。上式中的可用如下近似公式表达 (4.197)式中, . 取决于应力水平，可导出显式公式，具体内容参见有关文献欧文，塑性有限元。对于应力增量，有 (4.198)用位移增量表示时，有 (4.199)4.6.2.3 弹/粘塑性平衡方程 在任何瞬时的平衡方程应满足 (4.200)在时间增量中，平衡方程由增量形式给出，即 (4.201)在时步中出现的位移增量能为 (4.202) (4.203)这里，是切线刚度矩

15、阵，使用如下形式，即 (4.204)把位移增量代回式(4.199)，可以得到应力增量. 再有 (4.205) (4.206)且有 (4.207) (4.208)应力增量的计算是基于增量平衡方程(4.201)的线性化形式，累积所有这样的应力增量得到总应力是不正确的，并不会真实地满足平衡方程(4.200)。为此，可以采用计算残余平衡力的方法进行迭代求解，即 (4.209) 对于几何非线性问题，由位移计算得出。然后，残余平衡力迭加到下一个时间步的载荷增量上。在弹/粘塑性分析中，时间步长的选择十分重要，限于篇幅，这里不再加以讨论。9.3 材料非线性问题的动力有限元法在材料非线性问题的动力有限元分析中，

16、首先考虑材料的非线性本构关系，此外还要考虑是否包含大位移、大应变等因素，在实际工程分析中，需要对具体不同问题具体处理，目前没有统一的方法。在这里只限于考虑小应变条件下的材料非线性问题，且对有限元动平衡方程进行时间上的数值积分，以实现动力分析。根据虚功原理，可以导出结构材料在时刻的动平衡方程。平衡方程与材料本身的性质无关。对于每一个结点应满足相应的平衡方程： (4.211)其中为内阻力： (4.212)体力一致力为(为作用的体力矢量) (4.213)式中，为时刻的形函数。惯性力为(为密度) (4.214)阻尼力为(为阻尼系数阵)： (4.215)边界面力的一致力为(为面力矢量)： (4.216)

17、式中, 为受面力作用的边界与单元边界重合的部分。利用Gauss-Legendre乘积方法，依据单元形函数，可以对上述各式进行数值计算。设材料的粘弹性本构方程中各部分如下： 1) 弹性应变 (4.217)其中为弹性矩阵，设为常数矩阵。 2) 粘弹性应变增量粘弹性应变率可以表示为 (4.218)即 (4.219)在时刻为 (4.220)且有 (4.221)采用如下插分格式 (4.222)其中，为插分系数，如0，0.5，1等。则粘弹性应变增量与应力增量之间的关系为 (4.223) 3)粘塑性应变增量当满足一定的屈服准则时，材料产生粘塑性变形。采用类似的插分方程 (4.224)其中用Taylor级数来

18、近似 (4.225)其中由屈服准则和粘塑性流动准则确定 (4.226) 采用Ducker-Prager准则： (4.227)式中由材料粘聚力和内摩擦角定；为应力第一不变量；为应力偏量的第二不变量，当取关联的流动法则时有 (4.228)这样，由式(4.224)和式(4.225)得到 (4.229) 4)总应变增量由应变分解可知 (4.230)其中总应变增量可以由结点位移得到 (4.231)式中, 为位移应变矩阵；为单元结点位移增量。由式(4.223)、式(4.229)和式 (4.230)可得应力增量为： (4.232)是粘弹粘塑性模量 (4.233)5) 动态有限元刚度矩阵在任一时刻，可以证明，该时刻的切线刚度矩阵为 (4.234) 3. 显式时间积分法任一时刻的动平衡方程可以写成如下矩阵形式： (4.235)式中为总质量矩阵，为总阻尼矩阵，为内力的总矢量，为作用的体力、面力的一致结点力矢量。为简便起见，用中心差分公式对上述动平衡方程进行离散化。加速度、速度分别为： (4.236) (4.237)通过差分处理，在时刻