材料非线性有限元法

1、第四章材料非线性有限元法以上三章分别研究了线性弹性有限元法，材料非线性本构方程和非线性方程组解法，本章就可以研究材料非线性有限元法了。在材料非线性基本方程中，除第二章所述的本构方程外，与线性弹性一样，而非线性有限元法又归结为一系列线性弹性问题。因此，只要在第一章中改用第二章的本构方程，就可建立材料非线性有限元法的基本内容。4-1非线性弹性有限元法第二章提到，非线性弹性本构方程与形变理论弹塑性本构方程在形式上相同，所以与第章一样，这里也按塑性力学形变理论，研究非线性弹性有限元法，以便把二者统一起来。1非线性弹性基本方程为了便于以后直接引用，这里列出全量形式的非线性弹性(或形变理论弹塑性)基本方程

2、，并用矩矩阵表示。几何方程：本构方程：=D(2.13)pSH1+心平衡方程：在0内边界条件：=u(在上)utp(在上)U虚功方程：-J&ty+1ud0+1EutpldA=000A位能变分方程：n其中0=fW(eld/Q-ftp0-JtpdAaw(J)2非线性方程组的建立由于虚功方程本身不涉及材料性质，所以第一章由虚功方程得到的单元平衡方程(1.4)式8和总体平衡方程(1.1)0式9完全适用于非线性弹性(或形变理论弹塑性）问题。可见，只要把非线性弹性（或形变理论弹塑性）本构方程代入单元或总体的平衡方程，就可以建立非线性方程组。）1式3后，再把1）割线刚度方程仿照线性弹性有限元法，把（1.）1式3

3、后，再把1）31）3代入（1.）4式8便得单元割线刚度方程，即其中单元割线刚度矩阵k（.）1ttJbl/vseps而割线本构矩阵Dep）1式而割线本构矩阵Dep）1式3的推导，同样可得总体割线刚度方程IkK（u）1U=p仿照（即其中总体割线刚度矩阵，如（s1）式4所示。刚度矩阵ku)正td儿.)1LJe=1而总体节点载荷仍如（1）1式0所示。而总体节点载荷仍如（由（）式可知，总体割线刚度矩阵K）取决于各单元的等效应变；又由（s式可知，等效应变是由应变计算出来的；再由（3和（）式可知，应变与总体节点位移有关。可见，总体割线刚度矩阵K）是总体节点位移的函数，所s以总体割线刚度方程（）式是一个非线性

4、方程组。必须指出，建立非线性方程组（）式，只是为了说明非线性弹性（或形变理论弹塑性）有限元方程的非线性性质。实际求解时并不用（4.）4式。因为求解（4.）4式要用直接迭代法，而正如3指-出2，直接迭代法不但计算量太大，而且常常不收敛。（2）切线刚度矩阵由3-2-可-知3，-在6求解非线性方程组时，除上述直接迭代法外，都要用到切线刚度矩阵（至少要用到初始切线刚度矩阵e和）。为此，这0T0T里讨论一下建立非线性弹性（或形变理论弹塑性）有限元方程中的切线刚度矩阵问题。由（1.1）0和9（3.）1式1可知沁）=迟1力ee=1于是由（3.）1和0（4.）6式可得e=ldiddUdUdV由于由（e=ldi

5、ddUdUdV由于由（2.1），（61.）3和6（分别是d和8,有dj甬=）式，并考虑到符号o和8Id(&)Tepd,=b命=血所以把（）7式便得总体割线刚度矩阵，即所以把（）7式便得总体割线刚度矩阵，即K(u)T0式代入式e=l其中单元切线刚度矩阵1k（8）1=J时Id（JbI/v其中单元切线刚度矩阵1k（8）1=J时Id（JbI/vT8pe3）具有初应变理论或初应力的刚度方程仿照线性弹性有限元法，把形式上相同的（式便得单兀刚度方程，即0T）0式1代入（2.1）3式，0再,把（2.）1式3代入（qelk03=fe其中单元刚度矩阵k1和初应变，初应力节点载荷fe,0T8o,e仍分别如（o5）0

6、和）式所示。但要强调，这里0的含义是单元初始切线刚度矩阵;e8中的初应变80或fe中的初应力0随迭代过程而变。o仿照线性弹性有限元法/样可得总体刚度方程，即TOC o 1-5 h z0Tkb=p-!f0Ta其中总体刚度矩阵K和总体初应变、初应力节点载荷FF在形式上均与线性0Ta弹性有限元法相同。3等效应力、等效应变关系由（4.）1-1-（-4.）1式6可知，要建立并求解非线性弹性（或形变理论弹塑性）有限元方程，关键是要具体知道材料的本构矩阵。而由（2.）1和（）式可知，只要（）和（）式中的函数关系b=（）是已知的，那么本构矩阵就是显式的。根据单一曲线假设，彳和列勺关系与单向拉伸时相同a三1丿=

7、E11一&再考虑体积不可压缩条件（V讦/2），则/1a三）=3G1-匕b其中eG）取决于所采用的简化模型。理想塑性（见图）：厂0eG）1-s1+See图4-1线性强化塑性（见图4-2）：图4-23（）=0,（）E（E+E）+（E-硃s,（_）s（4.20）E幂次强化塑性（见图4-3）：ee图4-3G)=1&n-1E(0n14.2)eeee4迭代公式的具体化由于非线性弹性（或形变理论弹塑性）有限元方程一般都写成全量形式，所以这里只相应的列出几种迭代类型解法的具体迭代公式。（）-法由（3、（）和（）式以及（3）和（3）式，有4.2)234e=BA4.2)234nn=i()2+6(2+2+2)2n3

8、iijj122331ni，j1k()e=BTD()BdVnTepnTeK(u)丄(Ae)Tk()eAenTnTe=14.2)R(u)=瓦(Ae)tJBt()dV4.2)nne=1TOC o 1-5 h zu=K(u)-1(PR(u)(4)nnTn2)初应变迭代法由(2.)1、(02.)1和3(3.9)、9(3.1)0式1可知=Dt(4)epee=D-1s(4.2)所以仿照E0=，D-!p法，并考虑到()和(F二(Ae)tBtDe0dVen所以仿照E0=，D-!p法，并考虑到()和(F二(Ae)tBtDe0dVen，1 HYPERLINK l bookmark16e=1e HYPERLINK l

9、 bookmark39U=K-1(P+F)n0Ten-1）1式0，有n，1e=BAeUnnE=2(E-E)2+6(E2+E2+E2)n3iijj1223i,j=31E0=-D(E)-!npnn3)初应力迭代法由(2.)1、(02和3（）1、（83.1）2式0可知所以仿照=DEeps=De0=DEp法，并考虑到（和6（）2式7，有F二(Ae)Ten，1e=1dVn，1Un=K0,1(P-Fn-i)E=BAeUnnC-E)+6(iijji,j=12+E2+E2/122331=-D(E)Enpnn非线性弹性手算例题为了熟悉非线性弹性（或形变理论弹塑性）有限元法及其非线性方程组的求解过程,这里以图所示

10、的弹塑性拉压超静定问题为例，用法、初应变迭代法、初应力迭代法进行手算。其中用法的求解作较详细的叙述，以便了解非线性弹性所或形变理论弹塑性的有限元分析的全过程，而用其他方法的求解只给出主要计算过程和计算结果。在该拉压超静定杆中划分的节点和单元如图（）中所示。单元和分别由线性强化材料和线性弹性材料制成，如图（）和（）所示。两个单元的截面积均为1212a，长度均为l，弹性模量均为E，单元的强化模量为EE/2。节点所受集中载荷P3A,，其中,为单元的屈服极限。2ss用-法求解用-法求解（）非线性有限元方程的形成首先去掉两固定端约束，用其约束反力P和P代替,13所以p=pppy（）123由于u=uuuT

11、123u=uuy（b）12u=uuy2所以100Ay1010010Ay21000设单元e内任意一点位移为u=a+xdu二a+xk所以e=b=冬二=Bubx-xji其中几何矩阵B=1-11l单元节点位移向量ue=UU1ij这里的单元都是单向应力状态，即J1=Q1J2=Q2)s)sDD1=E1epTT1尹所以由（2tGdV,(所以由（2tGdV,(A2)tfBtq2dVi）把（）（）式代入（）和（）式，并考虑到JdV=JdV=Al,有E1E10,ATTE1E1+E2E2TTTT0E2E2TT710R(u)二A(-1+210V丿1V丿l并考虑到)K（叫二式代入（3.）1式7，_110AE200_11

12、0AE2000110+T1011000011AE1T把（）、（）和（）Du=Du1E1E10,uATT1E1E1+E2E2”u1TTTT20E2E2uTT丿3P10、Du2DuT，有31P一A(-1”1+-2”一1)2(1)12121处理后的非线性有限元方程E2T3(1处理后的非线性有限元方程E2T3(1一2)E1+E2sT2)切线刚度迭代公式的建立Eu1-1二E1=r(111E2u12T由于1)P3最后由边界条件DU=Du3=0和P2=3A，用第一章所述的消行降阶法便得约束-二Ess所以无量纲迭代公式为C1nCs1(u)C1nCs(u)(1)佗2n18八(u)1)2n1)1811J1+7(u

13、)(2182nsC2nCs18(u)2nRC1G2n=nnACC11J1+7(u)(2182nsC2nCs18(u)2nRC1G2n=nnACCssE1TE2T(u)1(72n1)181(us)1)(2-n1)218s(u)2n18s(ET1(3nE2)+1AC(u)_(u)(u)/、181818sss()节点位移和单元应力的求解按迭代公式()()式，求节点位移(u2)n+118仿照用)0式仿照用)0式3有和单元应力Cn+1CC；+1C的迭代过程和计算结果如表所示：ss表%C)s(u)/2/(18)s(役)n%sCn/sC2/n拎s用初应变迭代法求解若令()式中E=E=E，则TT110AEK0

14、T-1121011法的求解，由(ssssssssr1-r1-10uPFae11E1121u_P+F,l22E2011uPF1133JE3丿约束处理后可得AE2u二3A+Fl2se2因为只有单元具有塑性所以初应变节点载荷2只有由单元提供。由于F_(Ai)tBtE(e0)1dVe10-1、101,E(e0)1Al_1，001J0_1E(e0)1A_i-1所以由图4）可知Fe2二由图4）可知Fe2二E(eo)iA2(i一)E2iEssssssssssssssss因此(e0)(e0)1=E11ee1一Es于是由（）、（）于是由（）、（）（）和（）0式9可得无量纲迭代公式，即(E(E0)1n1Es1ns

15、n12An12As(E0)1n12esssssssssssssssssTOC o 1-5 h z(U)_3丄(F)lE22Ass_1(1+也)2lE2(u)n_一2nlE00按迭代公式（）（）式，求节点位移（U2）n（伦）和单元应力，n,，：,的迭代过sss程和计算结果如表所示：表30)1/nKs(F)/e2y(2A,)s叫）ss,2/n，s用初应力迭代法求解仿照用初应变迭代法求解，有1-1-10uPFAE11，1121VuP-F卜l22，20-11uPF1133，300约束处理后可得AE2u(,1-,i)Ae所以F，(1-1)A()2e于是由()、()、()和()式可得无量纲迭代公式，即11

16、门丄(u)n-1，(12n1)TOC o 1-5 h z24lEss(1)(u)cn12n1()22lEss(F)1(1)2n1n1c_n1(2A2ss(u)_3(F),、lE22Ass按迭代公式()()式，求节点位移(u)/(lc)和单元应力s/s、s/s的迭2nsnsns代过程和计算结果如表4-所3示：表(F)/2咲2A)s(u)/2/(lE)ss2/(U)/(2n(ie)ss4-3弹塑性有限元法在4-中1讨论了基于形变理论的弹塑性有限元法。这里将讨论基于流动理论的弹塑性有限元法。为了跟踪加载历史求出位移、应变和应力的全量，基于流动理论弹塑性有限元方程只能取增量形式。00同非线性弹性(或形

17、变理论弹塑性)问题一样，第一章由虚功方程得到的平衡方程也完全适用于流动理论的弹塑性问题。所以由(完全适用于流动理论的弹塑性问题。所以由(1)和1(3.)59式便得总体增量0000刚度方程，即K(u)dU=dPT其中总体切线刚度矩阵(参考(4.)1式1)K(u)区(Ae)Tk(G)eAe(TTe1单元切线刚度矩阵(参考(4.1)2式)k(G)e=JBTD(g)BdV(TepTe(4式中的切线弹塑性矩阵D(G)，按()式(正则屈服面)或(epT和(2.)8式8(非正则屈服面)计算。几种常见材料模型的切线弹塑性矩阵已在2-中给出显示形式。由于流动理论弹塑性有限元方程只能取增量形式，所以在第三章所述的

18、解法中，只有增量类型的解法才能用来对它进行求解，其中常用的是法、一次迭代增量类型的解法才能用来对它进行求解，其中常用的是法、一次迭代0000法、初应变增量法和初应力增量法。i逐步求解公式的具体化这里列出一求得G的第增量步的具体公式。i)4式6，并仿照非法对比()4式6，并仿照非法对比(线性弹塑性有限元法，l法逐步求解公式为:k(g)e=BTD(g)BdViTepiTeK(U)迟(Ae)Tk(G)eAeiTiTe1UK(u)-iAPiiTi門+i門+卩一次迭代法:同样，仿照非线性弹塑性有限元法，一次迭代法:同样，仿照非线性弹塑性有限元法，一次迭代法逐步求解公0000式为:式为:伙Qi)t=JBt

19、D(g)BdVepiTK(U)区(Ae)Tk(G)eAeiTiTe1R(u)二,(Ae)tJBtgdViie1TOC o 1-5 h z，U=K(u)_1(PR(u)()iiTi+1i=+,()i+1ii()初应变增量法仿照中的初应变迭代法，并考虑到()和()3.11式6，有( HYPERLINK l bookmark37,F=瓦(Ae)tJBtD,0dV(5ei1i1e=1e，U=K-i(AP+,F)()i0Tiei1UUi+i=Ui+，Ui0000As0=-D()-i,ipii对于限增量形式，则由（模型，若把塑性应变增量对于限增量形式，则由（模型，若把塑性应变增量ds直接看作是初应变增量d

20、s0，并写成有p）3和5（2.1）3式6可得：）、2（60000AsoAso=As=iLr1）Gr1）T，piHQiQiii初应变增量法：仿照式，有：（中1的初应力迭代法，并考虑到（3.1）2初应变增量法：仿照式，有：（,F=区(Ae)tJBt,0dVi-1i-1e=1，Ui=K0-1（，鬥厂，尸丄1）Ui+i=Ui+，Ui,0=,-D,siii2模型和1）单元弹、塑性状态的确定（模型由于同一增量步的不同单元可能处于不同的状态一弹性或塑性，而不同状态下单元的应力增量应按不同公式计算，所以要计算（）（6式中的单元应力增量，必须参照中的加卸载准则，确定第增量步单元的i弹、塑性状态。为了确定第增量步

21、单元的状态，可先假定该单元在此增量步内不产生新的塑性变形弹性、卸载或中性变载），于是有：=+，=+D，sei+1ieiii据此可以计算第步终了时的加载函数值(f*)或屈服函数值(f)。为了计算方便，可TOC o 1-5 h zi+1i+1把加载函数写成另一种形式。由()、()和()式可知，等向强化加载函数值为：11(f*)，(S2S2S22S22S22S2)-(Q)2()i12112233122331i13si1随动强化(线性)加载函数值为：11(f*)，-(S2S2S22S22S22S2)-Q2()i12112233122331i13S理想塑性加载函数(即屈服函数)值为：(f*)，(s2s2

22、s2)2(s22s22S2)i12112233i1122331i1(4.6)-92)26k(q)-k21i11i1参照2-中2的加卸载准则，若(f*)0i1则应力s在加载面内(弹性或卸载)或在加载面上(中性变载)，说明上述假定是正确ei+1的。我们把这种单元简称为弹性单元。若(f*)0i+1则应力Q在加载面外，这当然是不可能的，说明上述假定是错误的。但这种情况下的ei1单元只有两种可能状态：全塑性或弹塑性，而单元的确切状态又取决于第步终了时(或第步开始时)的单元状态。若第步终了时单元是塑性的，即(f*)，0i则第步单元是全塑性的(有一种塑性状态进入另一种塑性状态，即完全加载状态)，我们把这种单

23、元简称为塑性单元。若第步终了时单元是弹性的，即(f*)0i则第步单元是弹塑性的(有弹性状态进入塑性状态，即部分加载状态)，我们把这种单元简称为弹塑性单元(或过度单元)。()单元应力增量Ac的计算对于弹性单元，显然有：iAc，DAs(6ii对于塑性单元，当单元应变增量Ae较小时，有：iAc，D(Q)As()iepiTi当单元应变增量A较大时，单元应力增量应当用切线模量按增量法求下列积分的近似i值：Aj=,AsiD()d()i0ep对于弹塑性单元，单元应力增量由弹性和塑性两部分组成。若用表示弹性应力增量与总应力增量A之比参看图()和相应的一维问题图()则可通过满足加i载条件确定这个比例因子。确定出

24、以后，当单元应变增量De较小时，有iTOC o 1-5 h z=+rDAs()iiiA二rD+(1-r)D()As(7iepiTi当单元应变增量As较大时，单元应力增量应当用切线模量按增量法求下列积分的近似.值：显然，当A显然，当A二JrAsi|.0二rDAs+.时0，此式即为(Dds+JAeiD()dsrAs.印丁,AiD()dsrAs.ep丁7式。()比例因子的确定假定用A()比例因子的确定假定用ASe或ASe表示弹性应力偏张量，于是把(.jij6式中的S或S换(a)(a)为sij+rse或s+r为sij+rse或s+re，并令其值为零便得ijijijrA其中A、B和C,Mises模型等向

25、强化材料为TOC o 1-5 h zA(3訂)2+(s；?)2+(s#3)2+2sf2)2+(s；3)2+()21Bse(s)+se(s)+se(s)+2se(s)+se(s)+se(s)11111i2222i3333i1212i2323i3131i2C(s2+s2+s2)+2(s2+s2+s2)一()2112233i122313i3siMises模型随动强化(线化)材料为A(e)2+(e)2+(e)2+2e)2+(e)2+(e)21112233122331Be()+e()+e()+2e()+e()+e()11111i2222i3333i1212i2323i3131i2C(2+2+2)+2(2

26、+2+2)Q2112233i122313i3（4.75）（4.76（4.75）（4.76）（4.77）（4.78）（4.79）（4.80）（4.81）（4.82）4.83）（4.84）A(s)2+(s?)2+(s#3)2+2(As(2)2+(sg)2+(s#)2118。2(Qe)2BAse(s)+Ase(s)+Ase(s)+se(s)+Ase(s)+Ase(s)11111i2222i3333i1212i2323i3131i+aAqeQk-18a(Q)1i1iC(s2+s2+s2)+2(s2+s2+s2)+a(q)12k18a(Q)bk2112233i122313i1i1i（4）积分近似值的计算

27、（4.71）和（4.74）式中的积分很难积得显式。为此，我们可以把积妥对应的塑性变形部分再分为若干子增量，用分段线性的计算结果去逼近积分结果。由于（4.71）式是（4.74）的r=0的特殊情况，所以这里只讨论（4.74）式中积分的近似计算问题。若把（1-厂）么，,或（1-r）aq）再分为N个相等的子增量，（例如图4i-6中N=5）,则有4.85)图4-6由于初始子增量步开始时的单元应力为y=+r4.85)图4-6由于初始子增量步开始时的单元应力为y=+rdKa,.t4.86)第j子增量步开始时的单元应力为yj=+IdQj-ij亞epT4.87)所以ac=rbh,+“tQj)5，.epj=1T=

28、rb+r“tQj)，Nepj=1（5）第i增量步的单元计算过程在第I增量步,已知的是ii4.88)要计算的00是&和耳+1。若用0表示判定单元应变增量a,.大小的标准则计算过程如图4-7所示。3.Tresca模型由于Tresca模型的屈服面比较复杂，所以这里只限于平面应力问题的Euler-Cauchy解法。epTepT由25中的和tj可知，Tresca模型弹塑性矩阵各元素均与主应力方向epTepT00关。当单元开始进入塑性状态以后，主方向随载荷的增加而改变，因而每一增量步都要重新对它进行计算。为了使应力满足加载条件，可以按下述插值方法对它进行计算第i增量步的主方向,+1值。(1)按计算主应力和

29、塑性功用t()(塑性单元)和b(弹性单元)以iepiT及它们的组合(弹塑性单元)形成Ik(u)，求得位移增量AU以后，有iTias=IbIaAu/tAg二()As.tepiTig=j+Agtitg，=1G，+g)+2(g，)(sc)1i21122i21122ii12ii+g，)一1(a，-g22i21122i1“，=arctgi2)c2-s2)+2S)(sc)i12ii2L)12L-L丿n22iW二W二AWpipipi2）按计算主应力和塑性功仿照（1）中的计算，有iA讣Ag”.=b(，)as”.IIiepiTi”g.=5+Ag.IG二1(L+g)+丄(L_g)c2-S2)+2(a)sc)1i2

30、1122i21122ii12ii-G，)C2-S2)+20，)(sc)22iii21122i2c”二-G+)-G2i21122i21112ii4.89)4.90)4.91)4.92)4.93)4.94)4.95)4.96)4.97)4.98)4.99)4.100)Wpi=W=AWpipi4.101)3）按插值方法计算主方向由i、；、Q1），、Q2）、Q1）、(G)、G(W)2spi和（W）按一种插值方法求，使得用它求得的应力满足加载条件。spi边6，有i+1i+1W，)_(g)-g一pi1t丄一一gy一g)-gW，)+gW，)对于图26中的，gW，)gy.spi1，4.102)pipi0000

31、44弹塑性手算例题我们仍以42中的弹塑性拉压问题为例，用Euler-Cauchy法、Euler一次迭代法、初应变增量法和初应力增量法进行手算。为了反映应力与加载历史的相关性，除Euler-Cauchy法外，均用叙述方式。1)用Euler-Cauchy法求解仿照42中用Newton-Raphson法的求解，约束处理后可得E(1)E(1)0fAuAP、AT-E(1)TE(1)+E(2)E(2)约束处理后可得E(1)E(1)0=1AP，lT0TT-ETTET2Au32AP3(1)+E(2)kulTTR)=AP222若以总载荷丁3A“s为参考载荷并按两个增量步取A”广3，A”则有i=1:ET(1)=E

32、T(2)=E,(AP2)i=2A“S,(Au)2-=1Ssi=2:TE比“E(1)=当一-2(“(AP)=A“2i(Au)佗s于是无量纲逐步求解公式为(Au)2/(i二1)(i二2)T)er=0s0r=0f*S+仏)=0，r=-B+亍图一（）(u)(u)+(A0r=0f*S+仏)=0，r=-B+亍图一（）(u)(u)+(Au)2ill=2-i2illssU）O(1)illOs(u)2illl2”1ls(u)2_i!l丿(i二1)(i二2)V）ill=2illOlssW）(l)和单元应力O(1)O、O(2)Osi+1si+1s按逐步求解公式（W）(l)和单元应力O(1)O、O(2)Osi+1si

33、+1sssssssss00表44iAi(u)/(l,)2rs(u)/(l,)2i+1s61)/、i+1s(2)/、i+1s1231.0001.0001.000-1.0002130.6671.6671.333-1.667（2）用Euler一次迭代法求解参照（只）和（S）式，由（4.51）（4.55）式可得Euler一次迭代法的逐步求解公式，即(p2)ij(p)+(AP)(Au)2(Au)2il，s(E+E(2)TTi(P)2i+1As(一(2)issssssss00ssssssss00(u)(u)(Au)2_i+12_i+2_il，l，l，sss若以总载荷P3A2解过程和计算结果如下。i1:为参

34、考载荷，并按两个增量步取Ai若以总载荷P3A2解过程和计算结果如下。i1:为参考载荷，并按两个增量步取Ai3,A23,则逐步求(u)2i0,(P)0,(AP)=AP2Al，212112s(P)(P)+(AP)2A，222121s(1)1s(2)010s（ET1）（ET2）E（按线性弹性计算）(Au)-1l,(E+E(2)ATT1(P)(61)(2)22一11ssssssss00ssssssss001l,s(u)(u)(Au)1l,s2_22_1-+l,l,ss由于单元的屈服条件为1而此时单元的应力为(1)(1)(u)22l,1ssssssss00ssssssss00可见恰好满足屈服条件，所以在第1增量步单元处于弹性状态。-1,-1,ssssssssssss00i2:(AP)A,PAg2222s(P)i2:(AP)A,PAg2222s(P)(P)(AP)3Ag232222(G(1)2(uAi,(G(2)(u)2_2-2lsET(1)ET(2)E按线性弹性计算）(Au)22ls(P)(G一G(2)22(ET1)eT2)2