材料有限元分析课件-4..ppt

1、 第四章第四章 非线性问题的有限元法非线性问题的有限元法4.1. 概述概述 线弹性体系的基本特点：线弹性体系的基本特点： 应变和位移关系是线性的（几何方程应变和位移关系是线性的（几何方程)； 应力和应变关系是线性的（本构方程应力和应变关系是线性的（本构方程)； 变形前应力和体力关系是线性的（平衡方程变形前应力和体力关系是线性的（平衡方程 ) 线弹性体系成立的前提：线弹性体系成立的前提： 结点位移无限小；结点位移无限小； 材料的应力和应变关系满足虎克定律；材料的应力和应变关系满足虎克定律； 加载时边界条件保持不变。加载时边界条件保持不变。0pLDLuT材料非线性材料非线性： 材料自身的应力和应变

2、关系是非线性的。如，纯金属及其绝大材料自身的应力和应变关系是非线性的。如，纯金属及其绝大多数合金材料、高分子材料等。多数合金材料、高分子材料等。应变应变屈服点屈服点.材料极限材料极限塑性应变塑性应变几何非线性几何非线性： 结构的位移造成体系受力状态发生显著变化。如，金属构结构的位移造成体系受力状态发生显著变化。如，金属构件的塑性变形和蠕变等。件的塑性变形和蠕变等。 大位移小应变（如，塑料的热成形）；大位移小应变（如，塑料的热成形）； 大位移大应变（如，金属的压力加工）。大位移大应变（如，金属的压力加工）。AB边界非线性边界非线性： 边界条件非线性变化。例如模锻毛坯的接触问题边界条件非线性变化。

3、例如模锻毛坯的接触问题4.2. 材料非线性材料非线性 共同特点：材料特性随温度和时间变化。共同特点：材料特性随温度和时间变化。（1）弹塑性问题）弹塑性问题 加载后，如果载荷恒定，则材料的变形不随时间变化加载后，如果载荷恒定，则材料的变形不随时间变化; 但另但另一方面，当载荷增加到某个值时，材料的变形会出现屈服现象。一方面，当载荷增加到某个值时，材料的变形会出现屈服现象。（2）粘弹性和粘塑性问题）粘弹性和粘塑性问题 加载后，材料的变形随时间变化。加载后，材料的变形随时间变化。蠕变蠕变应力恒定，应变应力恒定，应变随时间增加；随时间增加；松弛松弛应变恒定，应力随时间减小。应变恒定，应力随时间减小。t

4、 = hourlt = 0 Creep l1 l2 l1 = l2 t = 0t = hourStressRelaxation 4.2.1. 非线性方程组的求解非线性方程组的求解 离散化的非线性方程组一般可表示为离散化的非线性方程组一般可表示为阵，未知场的函数非线性方程组的系数矩解未知函数（场）的近似外载荷式中或)()24(0)()()() 14()(aKaQffaaKfaPaQaaK返回到P13(1) 直接迭代法直接迭代法后，迭代终止。，即小于某个规定值一直到误差的某种范数次近似解重复上述迭代，可得。的某个解作为通常可借用线弹性问题，式中得一次近似解，代入近似（试探）解初值为设)54( )4

5、4()( )() 34()(a )24(1110001010rnnrnneaaeefKanaaKKfKa 直接迭代法求解非线性方程组的局限：只适用于求解与变直接迭代法求解非线性方程组的局限：只适用于求解与变形历史无关的非线性问题（例如，弹塑性问题）。形历史无关的非线性问题（例如，弹塑性问题）。的形成与历史无关。的选取和因为 )( 10naKa 直接迭代法收敛的几何含义（图中的直接迭代法收敛的几何含义（图中的 为标量，非线性为标量，非线性系统是单自由度的）系统是单自由度的）(2) NewtonRaphson迭代法迭代法 目的：进一步提高近似解的精度和解的收敛速度。目的：进一步提高近似解的精度和解

6、的收敛速度。 定的收敛条件为止。进行迭代，直到满足给对结合其中内，得为切线矩阵，将其代入式中公式展开附近按在将项存在，并且有式的假设)64()84()( )84()()()( )74( )()74(0)()( )()64( )( )24(111111 nnnTnTnnTnnTnTnnnnnnnnnnaPPaKKfPKKaaKdadPdadadadaaTayloraaaaaa回到P14 N-R法中的初始近似解法中的初始近似解 ，可简单的设为，可简单的设为 ；这样，；这样，的的初值初值在非线性问题中就是弹性刚度矩阵。在非线性问题中就是弹性刚度矩阵。 N-R法求解过程的几何表示法求解过程的几何表示0

7、a00anTK0TK(3) 修正的修正的NewtonRaphson法法 目的：克服每迭代一次都需重新生成并求逆切线矩阵的麻烦。目的：克服每迭代一次都需重新生成并求逆切线矩阵的麻烦。 思路之一：思路之一：)104()()()84()94(100fPKaKKnTnTnT得代入令 式式(4-10)是以牺牲收敛速度为代价来换取计算量的减少。是以牺牲收敛速度为代价来换取计算量的减少。思路之二：思路之二： 迭代若干次（例如迭代若干次（例如 m 次）后，更新次）后，更新 为为 ，再进行，再进行以后的迭代。以后的迭代。TKmTK(4) 增量法增量法 如果初始状态下，解向量如果初始状态下，解向量 和载荷向量和载

8、荷向量 f 均为定值，则均为定值，则用增量法求解非线性方程组可以得到较好的收敛解。用增量法求解非线性方程组可以得到较好的收敛解。 令式令式(4-2)中的载荷项中的载荷项 ，得增量方程，得增量方程的。和是针对或式中的足标利用欧拉公式改写或微分，有上式对 1 )144()()( )134()134()( )124(0 )114(0)(101101000fmmfKfaKaafaKddafddaKfddadadPfaPmmTmmTmmTT0ff解得解得(4-14)中的中的 ，再用，再用“修正的修正的”欧拉公式改进之，即欧拉公式改进之，即01ma) 10()1 ()()()154()(0111mmmmT

9、mTmmTmmaaaaKKfKaa式中 按按(4-14)或或(4-15)计算出计算出来的来的 （如图）一般会（如图）一般会导致解的漂移。为克服解的导致解的漂移。为克服解的漂移现象，可以将漂移现象，可以将N-R法或法或mN-R法用于每一增量步。法用于每一增量步。1ma例如，用例如，用(4-8)修正修正(4-14)，得，得)164(0)()()(110110111nmmnTmnmmnmaKfaPfaP由上式解出由上式解出)184( 1 )()()174()()(111110111011111nmnmnmmmmmTmnTmnmmnTnmaaanamaanKKfaPKa次迭代的改进值的第则最后可得次解

10、），（即增量法的第如果开始迭代时取次迭代的改进值的第是式中回到14，17，40 为缩减为缩减(4-17)的计算量，采用修正的的计算量，采用修正的NewtonRaphson方法，方法，此时此时。就是实际上，则在精确解，即存程增量步结束时，增量方若进一步假设，当前一式中只求解一次，则有如果)144()214( )214()( 0 )114()204()()()174()194()()()(01101011011101mmTmmmmmmmmmTmmmTmTmnTfKafPfPKaaaKKK (4-17)和和(4-20)被称为考虑平衡校正后的迭代算法，其几何意被称为考虑平衡校正后的迭代算法，其几何意义

11、如图所示义如图所示变斜率切线变斜率切线对应对应(4-17)恒斜率直线恒斜率直线对应对应(4-20)(5) 加速收敛的方法（以加速收敛的方法（以Aitken法为例）法为例） 考查单自由度非线性系统，无考查单自由度非线性系统，无Aitken加速的加速的mN-R迭代和有迭代和有Aitken加速的加速的mN-R迭代求解示意图如迭代求解示意图如(a)、(b)所示：所示：特点：切线和割线交替出现当经过迭代当经过迭代12次后，从次后，从(b)图得到的两次迭代之差值为图得到的两次迭代之差值为局部割线刚度矩阵起始切线刚度矩阵式中smTmmmTmsKKaaKaK )( )224()()(1101011111111

12、111111111111 1 )()()234( )234()( )( mmmmsmTmmmTmsmaaaabaaaaKKaaKaKaAitken所以为加速因子。，称可知，比值、或图由于是可得，即令来加速非线性解的收敛法利用于是，增量法公式于是，增量法公式(4-18)可改写成可改写成加速迭代其中 ), 5 , 3 , 1( ), 4 , 2 , 0(1 )244(111111111naaanaaaaanmnmnmnnmnnmnm 公式公式(4-24)表明，表明，Aitken加速收敛法的特点是：迭代和加加速收敛法的特点是：迭代和加速交叉进行。速交叉进行。将将(4-24)推广到推广到N个自由度的系

13、统个自由度的系统), 5 , 3 , 1(), 4 , 2 , 0(1 )254(1,11,111111naaaniaaanminminmnininnmnnmnm为个元素第是对角矩阵，该矩阵的其中， 为避免因上式分母项为避免因上式分母项 的值很小时，计算量剧的值很小时，计算量剧增的情况出现，特对增的情况出现，特对(4-25)进行修正，用标量代替进行修正，用标量代替(4-25)中的对角中的对角矩阵矩阵nminmiaa1,11,n加速收敛法。即为修正的此时的AitkennaaaaaaananminmiTnminminmTnminmin)254(), 5 , 3 , 1()()()(), 4 , 2

14、 , 0(11,11,1,11,111,11,4.2.2. 材料非线性本构关系材料非线性本构关系4.2.2.1. 材料的弹塑性行为材料的弹塑性行为 单调加载单调加载理想弹塑性硬化塑性非线性弹性和塑性当材料发生应变硬化（加工硬化）时，有当材料发生应变硬化（加工硬化）时，有)262()(pss即：加载过程中，材料的下一步屈服与前一步应变有关。 反向加载反向加载 针对硬化材料，如果在一针对硬化材料，如果在一个方向加载进入塑性后，卸载个方向加载进入塑性后，卸载并在反方向加载，直至进入新并在反方向加载，直至进入新的塑性。的塑性。混合硬化且运动硬化各向同性硬化则卸载应力反向加载的屈服应力正向加载的屈服应力

15、设 2 2 0111101111110ssrrsssrrsrss循环加载循环加载一次循环循环松弛循环硬化循环蠕变（棘轮效应）的过程）。的过程）。下一个下一个循环到循环到前前环加载的过程（即从当环加载的过程（即从当加载分支：完成一个循加载分支：完成一个循次加载的屈服应力次加载的屈服应力第第次反转加载的应力次反转加载的应力第第图中图中 sisisiriii 通常，循环加载条件与材料特性的关系通常，循环加载条件与材料特性的关系 循环加载条件循环加载条件材料特性材料特性等幅应变控制等幅应变控制 循环硬（软）化循环硬（软）化不等幅应变控制不等幅应变控制循环松弛循环松弛不等幅应力控制不等幅应力控制 循环蠕

16、变（棘轮效应）循环蠕变（棘轮效应）4.2.2.2. 塑性力学的基本法则塑性力学的基本法则(1)屈服条件屈服条件对于初始各向同性材料，其开始进入塑性流动的条件为对于初始各向同性材料，其开始进入塑性流动的条件为式中式中 k “硬化硬化”参数，参数， 应力向量阵列，应力向量阵列，Y(k) 单向屈服应力单向屈服应力)274( 0)()(),(kYfkF(2) 流动法则流动法则假设假设 塑性应变增量与塑性势（能）有关塑性应变增量与塑性势（能）有关则则式中式中塑性应变增量塑性应变增量 ； 与材料硬化法则有关的参数；与材料硬化法则有关的参数； Q 塑性势函数塑性势函数 )284( dQppd返回到26，33

17、， 35，38对于稳定的应变硬化材料，如果存在关联塑性，即对于稳定的应变硬化材料，如果存在关联塑性，即Q = F其中其中 F 后继屈服函数（后继加载函数或加载曲面，即与加载历后继屈服函数（后继加载函数或加载曲面，即与加载历史有关的屈服函数）史有关的屈服函数）此时，公式此时，公式(4-28)变成变成)294( dFp(3) 硬化法则硬化法则公式公式(4-29)中，后继屈服函数的一般表达式中，后继屈服函数的一般表达式对于理想弹塑性材料，因无硬化效应，所以后继屈服函数对于理想弹塑性材料，因无硬化效应，所以后继屈服函数应与初始屈服函数应与初始屈服函数 F( ) 一致，即一致，即)304( 0) , ,

18、(pF)314( 0)() , ,(FFp对于硬化材料对于硬化材料 各向同性硬化法则各向同性硬化法则各向同性硬化特点各向同性硬化特点材料进入塑性变形后，加载屈服面（后继屈服面）在各材料进入塑性变形后，加载屈服面（后继屈服面）在各方向上均匀向外扩张，其形状、中心和方位保持不变。方向上均匀向外扩张，其形状、中心和方位保持不变。 例如，当例如，当 时，初始屈服轨迹与后继屈服轨迹之间时，初始屈服轨迹与后继屈服轨迹之间的关系的关系03加载屈服面初始屈服面)(31 )334( , )(21 )324( 0)(3),(2/1222222zyxmmzzmyymxxxyzxyzzyxspssssssF，式中根据

19、根据Von.Mises流动法则流动法则(4-28)，各向同性硬化后继屈服函数，各向同性硬化后继屈服函数（加载屈服面）的通式为（加载屈服面）的通式为)(sp式式(4-32)中的中的 是加载时的后继屈服应力，它是等效塑是加载时的后继屈服应力，它是等效塑性应变性应变 的函数。其中的函数。其中)344( 322/1pTpppddd的值可由材料单轴拉伸试验的曲线获得，定义的值可由材料单轴拉伸试验的曲线获得，定义 为材料的塑性模量（硬化系数），它与弹性模量为材料的塑性模量（硬化系数），它与弹性模量 E 和切线模量和切线模量 之间存在关系之间存在关系)(ps)354( pspddEddE )364( EEE

20、EEp各向同性硬化法则适用的材料各向同性硬化法则适用的材料 单调加载，且单调加载，且 。11rs 运动硬化法则运动硬化法则 运动硬化法则特点运动硬化法则特点 材料进入塑性后，加载屈服面在应力空间作刚性移动，其形状、材料进入塑性后，加载屈服面在应力空间作刚性移动，其形状、大小和方位均保持不变。大小和方位均保持不变。 此时的后继屈服函数可表示此时的后继屈服函数可表示为为 F( , ) = 0 (4-37) 式中式中 是加载曲面中心在应是加载曲面中心在应力空间的移动张量力空间的移动张量(法向量），与法向量），与材料的硬化特性和变形历史有关。材料的硬化特性和变形历史有关。 运动硬化法则适用的材料运动硬

21、化法则适用的材料 单调加载时，单调加载时， ； 卸载时，卸载时， 。11rs0112ssr初始屈服面加载屈服面 混合硬化法则混合硬化法则 混合硬化的后继屈服函数混合硬化的后继屈服函数的参数等效应变各向同性硬化中类似于（与运动硬化有关）间的移动张量加载曲面中心在应力空式中 )384( 0),( pF混合硬化法则适用的材料混合硬化法则适用的材料 并且主要用于反向加载和循环加载的场合。并且主要用于反向加载和循环加载的场合。011112 ssrrs卸载时单调加载时(4)加载、卸载准则加载、卸载准则由于材料变形过程中，其应力状态是在变化的，因此，由于材料变形过程中，其应力状态是在变化的，因此，用加载、卸

22、载准则来判断材料从当前状态出发，下一步是继用加载、卸载准则来判断材料从当前状态出发，下一步是继续加载还是弹性卸载，并以此确定后继计算是采用弹塑性本续加载还是弹性卸载，并以此确定后继计算是采用弹塑性本构方程还是采用弹性本构方程。由构方程还是采用弹性本构方程。由(427)知知 0 , 0 0 , 0 0 , 0 dfFdfFdfF若载由塑性加载转向弹性卸则若继续塑性加载则若则 对于理想弹塑性材料，继续塑性加载；对于硬化材料：中性变载（即仍保持塑性状态，但不发生新的塑性流动）4.2.2.3. 应力应变关系（本构关系）应力应变关系（本构关系） 当应力产生无穷小增量时，假设应变由弹性应变和塑性应变两当应

23、力产生无穷小增量时，假设应变由弹性应变和塑性应变两部分组成，即部分组成，即因因 弹性应变弹性应变 塑性应变塑性应变)394 （pedddQddDdpe1)404 1（QdDd返回到32 当材料发生塑性屈服时，应力状态处在式当材料发生塑性屈服时，应力状态处在式(4-27)所表示的屈服面所表示的屈服面上，对上，对(4-27)微分，得微分，得)434( 1 )424( 0 )414( 0 dFAAdFdFdFdFdFdFdFdFdFTxyxyzxzxyzyzyyyyxx式中或)464( )454( )404( 0 )424()444()444( QDFADFQDDDdDdAQDFDdFQDFDdFd

24、FTTepepTTTTT其中，得，并代回从上式解出，有代入将为了消去参数为了消去参数 ，用，用 左乘式左乘式(4-40)两端，得两端，得DFT 式式(4-45)中的中的 称为称为弹塑性矩阵弹塑性矩阵。该矩阵只有在材料是关联。该矩阵只有在材料是关联 塑性的情况下才对称。此外，对于理想塑性材料（这时塑性的情况下才对称。此外，对于理想塑性材料（这时 A = 0），），矩阵矩阵 仍有定义。仍有定义。epDepD4.2.3. 粘塑性问题粘塑性问题4.2.3.1. 粘塑性材料的本构方程粘塑性材料的本构方程 对于具有粘塑性的材料，在应力空间中，其总应变速率等于对于具有粘塑性的材料，在应力空间中，其总应变速率

25、等于弹性应变速率与粘塑性应变速率之和。即弹性应变速率与粘塑性应变速率之和。即 粘塑性材料的屈服条件在形式上与塑性材料相同粘塑性材料的屈服条件在形式上与塑性材料相同（4-27），即即)474( vpe)484( 0)(),(),(kYfkFvpvp返回43根据粘塑性流动法则，材料的粘塑性应变速率可表示为根据粘塑性流动法则，材料的粘塑性应变速率可表示为屈服函数动率）粘度系数（控制塑性流塑性势函数式中 ),( )504( 0F 00F )()( ),( )494( )( vpvpvpFFFFQQQF)524( )( )514( 11QFDDDee为粘塑性材料的本构方程返回到42，434.2.3.2.

26、 蠕变蠕变 蠕变的特征：在常应力条件下，材料的变形与时间和温度有关。蠕变的特征：在常应力条件下，材料的变形与时间和温度有关。设设 蠕变应变为蠕变应变为 ，则蠕变应变率为，则蠕变应变率为c)534( ),(cct = hourlt = 0 Creep4.2.4. 温度对材料非线性本构方程的影响温度对材料非线性本构方程的影响 考虑温度的影响，则非线性材料的应变增量（用张量形式）考虑温度的影响，则非线性材料的应变增量（用张量形式）应为应为蠕变应变增量塑性应变增量；温度应变增量；弹性应变增量；式中 )544( cklpklkleklcklklpkleklklddddddddd 此时的应力应变关系可表示

27、为此时的应力应变关系可表示为列子矩阵行第弹性矩阵中的第式中 )554( )( jiDdddddDdeklijpklckleklklkleklijij上述应变增量的具体表示式：上述应变增量的具体表示式：(1)弹性应变增量弹性应变增量 考虑到温度对弹性模量考虑到温度对弹性模量 E 和泊松比和泊松比 的影响，有的影响，有量变化引起的弹性应变增由引起的弹性应变增量由温度增量列元素行第中的第表示矩阵足标弹性柔度张量，式中eklijklijeklekleklijkleklijeklijijeklijeklCddddjiCjiCdddCdCd , )564( (2) 塑性应变增量塑性应变增量)574( kl

28、klpklQdd(3) 温度应变增量温度应变增量 常数温度膨胀系数式中daltaKronecker 31 )584( 2klklkldd(4) 蠕变应变增量蠕变应变增量 平均正应力分量用张量表示的偏斜应力等效应变等效应力式中 )(31 , 23 ,23 )594( 23 3322111/21/2mijmijijijijeijijklecklsssdtsdtdad4.2.5. 弹塑性问题的有限元求解弹塑性问题的有限元求解(增量法增量法) 由于材料和结构的弹塑性行为与加载及变形历史有关，所以，由于材料和结构的弹塑性行为与加载及变形历史有关，所以，通常把载荷分解成若干个增量，针对每一个载荷增量，线性

29、化弹塑通常把载荷分解成若干个增量，针对每一个载荷增量，线性化弹塑性方程，从而将非线性问题转化成一系列线性问题（即按载荷步求性方程，从而将非线性问题转化成一系列线性问题（即按载荷步求解）。解）。 应用应用(4-16)，即，即段载荷步足标表示步迭代的切线矩阵第应变矩阵式中 1 )614( )( )604(0)()( 1111011mnKBdBaPaKfaPnTnmTnmnmmnTmnm然后利用迭代法可求得然后利用迭代法可求得弹塑性矩阵式中 )( 11111111111epnmnmepnmnmnmnmnmnmnmDDaaa4.2.6. 粘塑性问题的有限元求解粘塑性问题的有限元求解(增量法增量法) 按

30、时间间隔按时间间隔(增量增量) 求解。求解。 假设：在假设：在 时刻已求得结点位移时刻已求得结点位移 和应力和应力 ，且载荷向，且载荷向量量 已知，则已知，则(1) 应变增量应变增量 由式由式(4-49)表示的应变率法则，可求得表示的应变率法则，可求得 内产生的应变增量内产生的应变增量 式中式中mmmttt1mtmammf)624( mmmmvpmvpCt)634( mmmHtC参数 ,mvpmHmt(2)应力增量应力增量由由(4-51)的增量形式和的增量形式和(4-47)，得，得)644( )(mvpmmemDD用位移增量表示总应变增量用位移增量表示总应变增量应变矩阵式中 )654( BaB

31、mmm将将(4-62)和和(4-65)代入代入(4-64)，得，得)674( )()(D )664( )t( 111mmmvpmmmmmmCDDDCIaBD式中返回45（3）平衡方程）平衡方程 在任一瞬间在任一瞬间 ，结点应力和等效结点载荷之间满足平衡方程，结点应力和等效结点载荷之间满足平衡方程，用积分式表示用积分式表示)684( 0 )(mmTmfdBmt对于对于 ，有，有mt)694( 0 )(mmTmfdB将将(4-66)代入代入(4-69)，可得时间步长，可得时间步长 的位移增量的位移增量mt)704( )( )694( )( )704( )( 1dBDBKfdtDBQQKammTmmTmmmvpmTmmmmTm式中最后将最后将(4-70)代回代回(4-66)，得应力增量，得应力增量 ，从而有，从而有mmmmmmmaaa11 利用利用(4-64)、(4-65)得得)724( )714( 11mvpmvpmvpmmmmvpDaB于是有(4-72)即为即为“平衡状态平衡状态”下求解粘塑性问题的增量公式。下求解粘塑性问题的增量公式。