第10章(非线性有限元)

公式号、图号等第十章非线性动力有限元法当机械结构受到较大的外载荷，或受到持续时间较短的冲击载荷作用时，结构会产生过大的变形,以至于必须考虑结构几何大变形对结构整体刚度及固有频率的影响，即所谓的几何非线性影响。另外,对于多数非线性动力学问题，还需要考虑材料非线性、接触非线性等方面的影响。非线性动力学分析求解的基本方程有如下形式0PIuM(4.141)式中，KuuCI为粘性效应项，考虑阻尼、粘塑、粘弹等效应。P为外部激励。对于考虑各种非线性效应的动力学问题求解，需要对动力学方程进行直接时间积分。即非线性动力有限元分析具有如下特点：（1）问题分析过程需要考虑时间积分效应，不必做模态分析，不必提取固有频率；（2）采用直接积分方法求解非线性动力学方程，需要对时间作积分计算，因此计算量远远大于线性模态动力学方法；（3）非线性动力学分析中可以施加不同类型的载荷，包括结点力、非零位移、单元载荷；（4）在每个时间步上，进行质量、阻尼、及刚度的集成，采用完整矩阵，不涉及质量矩阵的近似；（5）可以同时考虑几何、材料和接触等多种非线性效应。非线性动力有限元分析程序常采用隐式Hilber-Hughes-Taylor法进行时间积分运算。这种方法适于模拟非线性结构的动态问题，对于冲击、地震等激发的结构动态响应以及一些由于塑性或粘性阻尼造成的能量耗散，隐式算法特别有效。隐式积分方法需要对刚度矩阵求逆计算，并通过多次迭代求解增量步平衡方程。隐式Hilber-Hughes-Taylor时间积分算法为无条件稳定，对时间步长没有特别的限制。采用子空间法也可以对动力学平衡方程作时间积分运算。子空间法是提取模态分析得到的各阶特征模态，并采用与线性模态动力学分析方法相近的分析方式进行求解。对于带有微小非线性效应的问题，如材料小范围进行入屈服、结点转角不大的情况，子空间法效率比进接积分法要高。此外，非线性动力有限元分析还可以采用显式动态算法，如中心差分法。显式时间积分算法为有条件稳定，其临界稳定时间步长限制了时间步长的大小，与有限元模型最小单元尺寸、材料应力波速等有关。显式时间积分法适于模拟高速冲击、接触等问题。上述方法的选择需要综合考虑计算量、分析问题的规模、单元限制等多方面因素，需要丰富的有限元模拟的理论、经验和实践知识。以下以几何非线性问题和材料非线性问题为例介绍非线性有限元法，其中粘弹粘塑性非线性材料问题的分析是典型的非线性动力有限元的求解思想。9.1几何非线性问题的有限元法几何非线性问题一般是指物体经历大的刚体位移和转动，但固连于物体坐标系中的应变分量仍假设为小量,即大位移小应变情况。4.6.1.1几何非线性问题的牛顿迭代法由数值分析技术可知，求解非线性方程组的数值方法的常规方法是Newton-Raphson法，即牛顿迭代法，这是一种近似线性化迭代求解方法。对于非线性方程0)(x，具有一阶导数，在nx点作一阶泰勒级数展开，它在nx点的线性近似为d()()()()dnnnxxxxx(4.142)因此，非线性方程0)(x在nx附近似为线性方程：d()()()0dnnnxxxx(4.143)当d()0dnx时，由上式求得n步的修正项1d()/()dnnnXxx(4.144)Newton-Raphson方法的迭代公式为11nnnXxX(4.145)在几何非线性有限元法中，结构的刚度矩阵与其几何位置有关，平衡方程由变形后的位形描述，因此，结构的刚度矩阵是几何变形的函数。设变形为δ,结构的平衡方程式()0KδδR(4.146)为一个非线性方程组。记非线性方程()0KψδδR(4.147)用Newton-Raphson方法求()0ψδ的根时，迭代公式分别为11nnnδδδ(4.148)其中,1nδ满足下式1()TnnnKδRKδδ(4.149)式中,TnK称为切线刚度矩阵，表达式为d()()dTnnψδKδ(4.150)在每一个迭代步中，通过求解切线刚度矩阵TnK，进而用1nδ进行迭代求解，称为Newton-Raphson方法，又称切线刚度法。牛顿法的收敛性是好的。但是某些非线性问题中，使用牛顿法迭代时，若TnK出现奇异或病态，则对TK的求逆出现困难。关于这一点也可以采用其它修正办法，如引入阻尼因子。对于已经建立的有限元方程，设ψ表示内为和外力矢量的总和，有***d0TTTVVδψεσδR(4.151)式中,R为载荷列阵；*δ为虚位移；*ε为虚应变用应变的增量形式ddεBδ代入上式，消去*δ项，可以得到非线性问题的一般平衡方程式为()d0TVVψδBσR(4.152)该式不论位移或应变的大小与否均成立。在有限变形中，应变和位移之间的关系是非线性的，即B矩阵是δ的非线性函数。但是，近似地可将进行如下分解：0LBBB(4.153)式中,0B为线性应变分析的部分；LB为由非线性变形引起的，与δ有关。假定应力应变关系为线弹性，于是有00()σDεεσ(4.154)式中][D为材料的弹性矩阵；}{0为初应变列阵；}{0为初应力列阵对于式(4.152)的非线性平衡方程式，可用Newton-Raphson方法进行迭代求解。对该式微分，有dddddTTVVVVψBσBσ(4.155)不考虑初应变和初应力的影响，得dddσDεDBδ并且ddLBB这样可得ddddTLVVψBσKδ(4.156)这里0dTLVVKBDBKK(4.157)式中0K为通常的小位移的线性刚度矩阵。LK矩阵则是由于大位移引起，它可以写成00dTTTLLLLLVVKBDBBDBBDB(4.158)式(4.156)又可记成：0dddLTψKKKδKδ(4.159)式中ddTLVVKBσ(4.160)式中，K是关于应力水平的对称矩阵，称之为初应力矩阵或几何刚度矩阵。因此，用Newton-Raphson方法迭代求解几何非线性问题的步骤为：(1)用线弹性解作为1δ，即一次近似；(2)通过定义1()B求出1σ，求出1ψ；(3)确定切线刚度矩阵1TK；(4)211/TδψK，212δδδ；(5)重复上述迭代步骤，直至nψ足够小。在这里，没有考虑载荷R可能由于变形而发生的变化，即在这里假设了载荷不因变形而改变其大小和方向，否则是非保守力作用下的大变形问题，在此不做讨论。4.6.1.2典型单元的切线刚度矩阵求解具体的几何非线性问题时，必须计算单元的切线刚度矩阵。对于一般空间问题，无论位移和应变大小，都可以利用应变的基本定义写出位移和应变的关系式。用变形前的坐标),,(zyx做为自变量，可以用位移wvu,,定义如下大变形问题的应变分量表达式222)()()(21xwxvxuxux2221()()()2yvuvwyyyy222)()()(21zwzvzuzwzzwywzvyvzuyuywzvxy(4.161)zwxwzvxvzuxuzuxwyzywxwyvxvyuxuxvyuzx对于微小位移情况，可以略去二次以上的偏导数项，得到小变形时的应变公式。在有限变形中，假设应变仍为小量。应变和位移之间的关系为：0Lεεε(4.162)式中0ε为线性应变部分。对于非线性部分，可以写成：0000001122000TxTyxTzLyTTzyTTzzxTTyxθθθθεθCθθθθθθθθ(4.163)式中[][][]TxTyTzuvwxxxuvwyyyuvwzzzθθθ(4.164)式中C为96矩阵。根据θ的定义，可以将θ表示成任意一点位移的函数，引入形函数N后，可以得到eθGδ(4.165)对于(4.163)式进行微分，得11dddd22LεCθCθCθ(4.166)因此，deeLLεCGδBδ(4.167)B矩阵为0LBBB(4.168)这样得到0dTLVVKKBDB(4.169)另外，有dddddeTTTLVVVVKδBσGCσ(4.170)利用矩阵C和列阵θ的性质，得到dddxxyzxTeyxyyzzxzyzIIICσIIIθMGδIII(4.171)式中I为三阶单位矩阵，M是99的六个应力分量组成的矩阵。因此几何刚度矩阵为dTVVkGMG(4.172)故此，非线性三维单元的切线刚度矩阵为0TLkkkk(4.173)作为特例，可以直接写出三角形单元的上述有关表达式。由三角形单元的位移模式eeijmuNNNvfNδIIIδ(4.174)其中，2/)(ycxbaNiiii，2/)(ycxbaNjjjj，2/)(ycxbaNmmmm，Tmmjjiievuvuvu}{}{，式中的iiicba,,等由结点坐标确定，为三角形单元的面积。根据式(4.164)Txyuvuvxxyyθθθ把式(4.174)代入上式得00000010002000ijmijmeijmijmbbbbbbccccccθδ(4.175)由(4.165)式可以知道00000010002000ijmijmijmijmbbbbbbccccccG(4.176)根据定义，由式(4.163)确定的平面问题的C矩阵为000000TxTyTTyxuvxxuuyyuvuvyyxxθCθθθ(4.177)这样可以得到C的显式为001002iijjmmiijjmmiijjmmiijjmmiijjmmiijjmmiijjmmiijjmmbubububvbvbvcucucucvcvcvcucucucvcvcvbubububvbvbvC(4.178)故LBCG(4.179)而0B由线性问题给出，即000010002ijmijmiijjmmbbbccccbcbcbB(4.180)至此，0B、LB和G都是常数矩阵，只与单元结点坐标和结点位移有关。因此，线性刚度矩阵为：000TtkBDB(4.181)式中t为单元厚度。初始刚度矩阵为：00()TTTLLLLLtkBDBBDBBDB(4.182)几何刚度矩阵为：TtkGMG(4.183)式中00000000xxyxxyxyyxyyM(4.184)因此，对于几何非线性问题，平面三角形单元的切线刚度矩阵可以由式(4.173)求出。9.2材料非线性问题的有限元法材料非线性是指材料的本构方程是非线性的。一般主要分为两类:一类是非线性弹性问题，如橡胶、塑料、岩石等，在加载时应力应变关系等性质呈现非线性的物理现象，卸载时可逆。另一类是指材料的弹塑性问题。材料超过屈服极限后呈现出非线性。在机械结构分析中，常见的本构方程主要有线弹性和非线性弹性模型，其特点是应力仅应变的函数，加卸载规律相同。公式如下：klijklijC(4.185)其中，对于线弹性材料为常数，对于非线性弹性材料，ijklC是kl的函数。此外，还有超弹性模型、次弹性模型、弹塑性模型等。对于弹塑性模型，可以认为弹塑材料发生塑性变形时，其总应变可以分解为两部分：pijeijij(4.186)即总应变为弹性应变和塑性应变之和。加载时遵循一定规律，如P