有关有限元的其他一些经典abaqus2用户单元子程序

1、20ABAQUS 用户单元子程序(UEL)在这一章中将列举两个在这些年里发展过的 ABAQUS/Standard 用户单元子程序（UEL）。第一个例子是一个非线性的索单元，的目的是通过这个比较简单的例子让读者了解用户单元子程序的基本开发过程；第二个例子是一个用于计算应变梯度理论的单元，应变梯度是比较热点的一个科研前沿问题，有各种理论，为了验证新的理论，需要数值结果与实验对照来进行评价，整个例子的目的是通过它说明用户子单元可以求解范围很广，但是由于内容比较艰深，程序也很长，所以这个例子并没有给出最后的全部程序。另外，到目前为止，ABAQUS 还只有隐式求解器 ABAQUS/Standard 支持

2、用户自定义单元，而显式求解器 ABAQUS/Explic it 中还不支持这能。20.1非线性索单元20.1.1 背景钢索斜拉桥和斜结构广泛应用于土木工程建筑上。索力的计算分析是设计和施工程力学系在采用 ABAQUS 进行工的关键环节。真分析（这个项目长江斜拉桥的计算机仿已在第 15 章“ABAQUS 在土木工程中的应用（一）荆州长江大桥斜拉桥结构三维仿真分析”中过）时，也曾进行了自行建立索单元的尝试。本节介绍的就是这方面的工作。理工大学土木与结构工程系采用 ABAQUS 有限元进行计算，完成了Ting Kau 斜拉桥和 Tsing Ma 悬索桥的结构计算和分析。对于钢索计算，他们采用元进行模

3、拟。由于元含有弯曲刚度，计算的高阶频率值偏高，周期较低。一般假设索是单向受拉力的构件。随着应变的非线性增加，索力呈非线性增加。尽管 ABAQUS 单元库中有 500 个以上的单元类型，但是，还没有索单元。本文发展了三维非线性索单元模型，形成 ABAQUS 的用户单元子程序，可以利用 ABAQUS 输入文件调入到具体的分析中。通过静态和动态例题的计算比较，索单元工作良好。20-120.1.2 基本公式在三维索单元计算中，如图 20-1 所示，坐标 x 和位移 u 的变量表达式为：i(x,y,z)(u,v,w)(20-1)u u ujiji应变的公式为：v z w v2 w 1 x1 u yu2j

4、i2(20-2)L ji jiji jijijijiji2公式（20-2）中，L 为索的长度，索的张力为：N AE N0在公式（20-3）中，A 为截面面积，E 为弹性模量，N0 为初始张力。(20-3)在总体坐标系下，单元刚度矩阵为：K K Ke(20-4) KK66单元刚度矩阵中的子阵 K 分别由线性和非线性矩阵项组成：K33 KL,33 KNL, 33(20-5)20-2图 20-1 索单元在公式（20-5）中的 KL 和 KNL 均是 33 的对称矩阵，分别为：x21010 x yx zjiji jiji ji N L EA K y2KNL0y zji ji LjiL31z2ji索单元

5、的节点质量为：m 1 AL(20-6)2在公式（20-6）中， 为密度。索单元的质量矩阵为： m0 M e(20-7) 0m结构的运动方程为：Mu F ext Ku(20-8)公式（21-8）中 Fext 为作用在结构上的外力。在不断变化的索的变形中，求解运动方程，得到节点的位移值。20.1.3 应用举例由 5 个单元组成的两端铰接的索杆结构，高 5m 长 10m，6 个节点号码依次为101106，如图 19-2 所示。计算振动的频率和周期。20-3图 19-2 由五个单元组成的两端铰接的索杆结构输入文件中的用户单元界面ABAQUS 输入文件（.inp）中的用户单元界面如下：*HEADINGT

6、wo dimenal overhead hoist frameusing 2 nodef-developed truement,Initial force N is defined in property(5) andreferenced by user element SI Units1-axis horizontal, 2-axis vertical*USERELEMENT,NODES=2,TYPE=U1,PROPERTIES=5,COORDINATES=3, VARIABLES=121,2,3*UEL PROPERTY, ELSET=UTRUSS 1.963E-5, 2.0E11, 0.

7、3, 7800, 10.0E5*ELEMENT, TYPE=U1, ELSET=UTRUSS计算结果和比较表 20-1ABAQUS列出了由用户索单元计算的图 20-2 所示结构的固有周期，并与应用元 B31 的计算结果进行了比较。元前 4 阶模态的周期基本一致；索单元的第 69 阶模态与索单元与元第710 阶模态的周期基本一致。从第 11 阶模态开始，随着元弯曲变形的增加，梁的弯曲刚度逐渐发挥作用并和轴向刚度耦合，与同阶模态的索单元相比，期显著降低，而频率高于索单元。元的振动周表 20-1 ABAQUS 用户索单元和元 B31 计算的频率比较20-4用户开发单元的缺点是不能采用 ABAQUS

8、的后处理进行显示，只能从数据文件（.dat）中结果。另外，ABAQUS 的接触算法等某些功能也无法应用。非线性索单元用户子程序subroutine uel(rhs, amatrx, svars, energy, ndofel, nrhs, nsvars,props, nprops, coords, mcrd, nnode, u,du,v,a,jtype,time,dtime,kstep,kinc,jelem,params,ndload,jdltyp,adlmag,predef,npredf,* lflags,mlvarx,ddlmag,mdlonewdt,jprops,njprop,perio

9、d)CInclude aba_param.incAll coordinates in globalCCdimenrhs(mlvarx,*), amatrx(ndofel, ndofel),svars(12), energy(8), props(5), coords(mcrd,nnode),u(ndofel), du(mlvarx,*), v(ndofel), a(ndofel), time(2),20-5振动模态索单元固有周期(Cycles/time)元固有周期(Cycles/time)1112.42112.482112.42112.483213.83213.954213.83213.9552

10、49.51222.556294.32222.757294.32294.488345.99294.489345.99346.1910474.60346.1911653.22422.3712767.91423.69*params(3),jdltyp(mdload,*), adlmag(mdload,*),ddlmag(mdload,*), predef(2,npredf,nnode),lflags(*), jprops(*)Cdimensresid(6), uji(3), xji(3), smatrx(3,3)CCMaterial propertiesarea=props(1)e=props(2)

11、anu =props(3)rho =props(4) Initial ten fn0 =props(5)Cforce in user elementCCGeometry, stiffness and mass parametersdxdy=coords(1,2)-coords(1,1)=coords(2,2)-coords(2,1)dz=coords(3,2)-coords(3,1)alen=sqrt(dx*dx+dy*dy+dz*dz)ang =atan(dy/dx)ak =area*e/alenam =0.5d0*area*rho*alenC=1,3uji(i)=u(i+3)-u(i) x

12、ji(i)=coords(i,2)-coords(i,1) end dostrain=(xji(1)*uji(1)+xji(2)*uji(2)+xji(3)*uji(3)+0.5*(uji(1)*2+uji(2)*2+uji(3)*2)/alen tforce=e*area*strain+fn0*20-6eal=e*area/alen*3CCStiffness matrix parameterssmatrx(1,1)=eal*xji(1)*2+tforce/alen smatrx(1,2)=eal*xji(1)*xji(2) smatrx(1,3)=eal*xji(1)*xji(3) smat

13、rx(2,1)=smatrx(1,2) smatrx(2,2)=eal*xji(2)*2+tforce/alen smatrx(2,3)=eal*xji(2)*xji(3) smatrx(3,1)=smatrx(1,3) smatrx(3,2)=smatrx(2,3)smatrx(3,3)=eal*xji(3)*2+tforce/alenCdo 6 k1=1,ndofelsresid(k1)= 0.0d0do 2 krhs=1,nrhsrhs(k1,krhs)=0.0d02continuedo 4 k2=1,ndofelamatrx(k2,k1)=0.0d0continue46Ccontinu

14、eif (lflags(3).eq.1) thenCNormal incremenionif (lflags(1).eq.1.or.lflags(1).eq.2) thenCC*SicElement stiffness matrix=1,3do j=1,320-7amatrx(i,j)=smatrx(i,j)amatrx(i,j+3)=-smatrx(i,j) amatrx(i+3,j)=-smatrx(i,j) amatrx(i+3,j+3)=smatrx(i,j) end doend doCCReaction forceif (lflags(4).ne.0) then=1,3force =

15、ak*(u(i+3)-u(i)dforce=ak*(du(i+3,1)-du(i,1) sresid(i)=-dforce sresid(i+3)=dforcerhs(i,1)=rhs(i,1)-sresid(i)rhs(i+3,1)=rhs(i+3,1)-sresid(i+3)end doelsedo k=1,3force =ak*(u(k+3)-u(k)sresid(k)=-forcesresid(k+3)=forcerhs(k,1)=rhs(k,1)-sresid(k) rhs(k+3,1)=rhs(k+3,1)-sresid(k+3) end doend ifCelse if(lfla

16、gs(1).eq.11.or.lflags(1).eq.12) then*Dynamicalpha =params(1)C20-8beta =params(2)gamma =params(3)Cdadu =1.0d0/(beta*dtime*2)dvdu =gamma/(beta*dtime)Cdo 14 k1=1,ndofelamatrx(k1,k1)=am*dadu rhs(k1,1)=rhs(k1,1)-am*a(k1) continue=1,3do j=1,3amatrx(i,j)=amatrx(i,i)+(1.0d0 + alpha)*smatrx(i,j) amatrx(i+3,j

17、+3)=amatrx(i+3,j+3)+(1.0d0 + alpha)*smatrx(i,j)amatrx(i,j+3)=amatrx(i,j+3)-(1.0d0 + alpha)*smatrx(i,j)14amatrx(i+3,j)=amatrx(i+3,j)-(1.0d0 + alpha)*smatrx(i,j)end doend doC=1,3force =ak*(u(i+3)-u(i)sresid(i)= -forcesresid(i+3)= forcerhs(i,1)=rhs(i,1)-(1.0d0+alpha)*sresid(i)-alpha*svars(i)rhs(i+3,1)=

18、rhs(i+3,1)-(1.0d0+alpha)*sresid(i+3)-alpha*svars(i+3)end do*Cdo 16 k1=1,ndofel20-9svars(k1+6)=svars(k1)svars(k1)=sresid(k1) continueend if16Celse if(lflags(3).eq.2) thenStiffness matrix=1,3do j=1,3 amatrx(i,j)=smatrx(i,j) amatrx(i,j+3)=-smatrx(i,j) amatrx(i+3,j)=-smatrx(i,j) amatrx(i+3,j+3)=smatrx(i

19、,j) end doend doCCelse if(lflags(3).eq.4) thenCMass matrixdo 40 k1=1,ndofelamatrx(k1,k1)=amcontinue40else if(lflags(3).eq.5) thenCCHalf-step residual calculationalpha = params(1)=1,3force =ak*(u(i+3)-u(i) sresid(i)= -force sresid(i+3)= forcerhs(i,1)=rhs(i,1)-am*a(i)-(1.0d0+alpha)*sresid(i)20-10*+0.5

20、d0*alpha*(svars(i)+svars(i+6)rhs(i+3,1)=rhs(i+3,1)-am*a(i+3)-(1.0d0+alpha)*sresid(i+3)+0.5d0*alpha*(svars(i+3)+svars(i+9)end do*Celse if(lflags(3).eq.6) thenInitial acceleration calculation do 60 k1=1,ndofel amatrx(k1,k1)=amcontinue=1,3force =ak*(u(i+3)-u(i) sresid(i)= -force sresid(i+3)= forcerhs(i

21、,1)=rhs(i,1)-sresid(i)C60rhs(i+3,1)=rhs(i+3,1)-sresid(i+3)end doCdo 62 k1=1,ndofelsvars(k1)=sresid(k1)continue62Celse if(lflags(3).eq.100) thenOutput for perturbationsif (lflags(1).eq.1.or.lflags(1).eq.2) thenCC*sic=1,3force =ak*(u(i+3)-u(i)dforce =ak*(du(i+3,1)-du(i,1) sresid(i)=-dforce20-11sresid(

22、i+3)=dforcerhs(i,1)=rhs(i,1)-sresid(i) rhs(i+3,1)=rhs(i+3,1)-sresid(i+3) end doCdo kvar=1,nsvarssvars(kvar)=0.0d0 end dodo j=1,6 svars(j)=rhs(j,1) end doelse if(lflags(1).eq.41) then*Frequencydo 90 krhs=1,nrhs=1,3Cdforce =ak*(du(i+3,krhs)-du(i,krhs)sresid(i)=-dforcesresid(i+3)=dforcerhs(i,krhs)=rhs(

23、i,krhs)-sresid(i)rhs(i+3,krhs)=rhs(i+3,krhs)-sresid(i+3)end do90continuedo kvar=1,nsvars svars(kvar)=0.0d0 end dodo j=1,6 svars(j)=rhs(j,1)end doCend if20-12end ifCreturnend20.2利用 ABAQUS 用户单元计算应变梯度塑性问题在本节内容中主要介绍两种应变梯度理论，并在最后给出用这两种应变梯度理论编写的用户单元子程序的数值计算结果与实验的对比。本节的目的在于让读者了解 ABAQUS 即使在面对如此复杂的理论问题时，也可以胜

24、任。20.2.1 引言研究应变梯度理论的意义很多试验表明，当非均匀塑性变形特征长度在微米量级时，材料具有很强的尺度效应。例如 Fleck 等在细铜丝的扭转试验中观察到，当铜丝的直径为 12m 时，无量纲的扭转硬化将增加至 170m 直径时的 3 倍；Stolken 和 Evans 在薄梁弯曲试验中也观察到，当梁的厚度从 100m 减至 12.5m 时，无量纲的弯曲硬化也显著增加；而在单轴拉伸情况这种尺度效应并不存在。在微米量级的尺度下，微观硬度试验与颗粒增强金属基复合材料中也观察到尺度效应，当压痕深度从 10m 减至 1m 时，金属的硬度增加了一倍；对于碳化硅颗粒加强的铝硅基复合材料，Lloy

25、d 观察到当保持颗粒体积比为 15%的条件下，将颗粒直径从 16m 减为 7.5m 后复合材料的强度显著增加。由于在传统的塑性理论中，本构模型不包含任何尺度，所以它不能尺度效应。然而，在工程实践中迫切需要处理微米量级的设计和制造问题，例如，厚度在 1m 或者更小尺寸下的薄膜；整个系统尺寸不超过 10m 的传感器、执行器和微电力系统（MEMS）；零尺寸小于 10m 的微电子封装；颗粒或者的尺寸在微米量级的先进复合材料及微加工等等。现在的设计方法，限元方法（FEM）和计算机辅助设计（CAD），都是基于经典的塑性理论，而它们在这一微小尺度不再适用。另一方面，现在按照量子力学和原子模拟的方法在现实的时

26、间和长度的尺度下处理微米尺度的结构依然很。所以，建立连续介质框架下、考虑尺寸效应的本构模型就成为联系经典塑性力学和原子模拟之间必要的桥梁。20-13促使建立细观尺寸下连续介质理论的另一个目的是在韧性材料的宏观断裂行为和原子断裂过程之间建立联系。在一系列值得注意的试验中，Eler 等测量了单晶铌/蓝宝石界面的宏观断裂韧度和原子分离功，使原子点阵或强界面分离所需要的力约为0.03E 或者10（Y E 为弹性模量，Y 为拉伸屈服应力）。而按照经典的塑性理论，Hhinson裂纹前方最大应力水平只能达到 4 至 5 倍Y。很明显这远远小于 Eler 等在试验中观察到的结果，使原子分离。考虑应变梯度的影响

27、有望解释这一现象。应变梯度理论简介目前发展的应变梯度理论有很多种，包括 CS 理论（偶应力理论）、SG 理论（拉伸和旋转梯度理论）、MSG 理论（基于细观机制的应变梯度塑性理论）以及理论（基于 Taylor 关系的非局部应变梯度理论）等。利用 ABAQUS 用户单元主要进行了MSG 和形变理论和两种理论应变梯度塑性的有限元分析，对于 MSG 塑性和塑性都包括理论，塑性还包括了有限变形问题的形变和理论的分析。现在对 MSG 理论和理论加以简单的介绍。20.2.2 两种应变梯度理论基于细观机制的 MSG 应变梯度塑性理论基于位错机制的 MSG 应变梯度塑性理论是由位错理论出发的，它通过一个多尺度、

28、分层次的框架由微观位错机制推导出了宏细观的应变梯度塑性理论。这个理论相比于其它理论，物理机制更明确，构造方法系统，而且第一次提出了材料长度的表达式。图 20-3是 MSG 应变梯度塑性理论的原理图，在微观层次上塑性是由位错运动产生的，在细观层次上引入应变的梯度与微观的几何必须位错密度相关过细观和微观的功等效由微观塑导出细观的本构理论。为了在细观尺度下的应变梯度塑性和微尺度下的 Taylor 硬化关系之间建立联系，在MSG 理论框架中采用如下的基本假设：20-141）假设微尺度的化关系应力由位错运动控制，并且遵守应变梯度律给出的 Taylor 硬 Yf 2 l(20-9)2）微观尺度和细观尺度的

29、联系是塑性功相等 dV V(20-10)ijijijijijkijkcellVcell3）在微尺度胞元中假设经典塑性的基本结立，其 J2 形变理论可以表示为 23 (20-11)ijije其中 e ， 。微尺度的屈服条件为322ij ijij ij3e (20-12)基于以上的理论假设，应变梯度塑性 MSG 形变理论本构关系可以建立如下：20-15图 20-3 MSG 理论中采用的多尺度框架ij 2ij 3(20-13) l f f 22(20-14)ijkijkijkYijk其中 Yf l2(20-15)式中：2l 182 b(20-16)ref为材料特征长度。 2 72 1(20-17)i

30、jkijkkjikijik ppjjk ppi4 1 542(20-18)ijkmnikjmnjk imnik jpjk ippmn4其中应变梯度表示为ijk uk,ij ik, j jk,i ij,k(20-19)1 (20-20)ijk ijk4基于 Taylor 关系的非局部应变梯度塑性理论（理论）MSG 应变梯度理论物理机制明确、构造方法系统，那么为什么还要发展理论呢？因为前面提到的应变梯度塑性理论，无论 CS、SG 还是 MSG，都是高阶理论，引入了高阶应力和附加的边界条件，而这些高阶应力和附加的边界条件都难以测量，难以想象，因此无法得到工业界的认可，难以实用。经典的塑性问题控制方程

31、是 2 阶，而在 MSG 理论中由于高阶应力的影响，其控制20-16方程是 4 阶，这就增加了解决问题的复杂性。非局部连续介质力学理论给以启发，采用应变的非局部积分来确定应变的梯度，这样就有了 基于 TAYLOR关系的非局部塑性理论。这种理论既保持 MSG 理论的所有优点，同时与经典的塑性理论相比又不增加方程的阶数。但也正是由于理论的非局部性质，使其在求取理论尤其重要。解方面比较，也正因为如此，有限元作为非局部塑性理论，有三个基本特点：（1）应力遵从 TAYLOR 硬化关系，这是塑性理论的出发点。（2）应变梯度和几何必须位错密度是非局部量，表示为应变的平均。这是塑性理论的概念。（3）塑性理论保

32、持经典塑性理论的基本结构。这个特点使得塑性理论具有很强的实用潜力。和经典的塑性理论相比，塑性理论的特别之处在于屈服条件的不同：应力不仅依赖于应变，还同时依赖于应变的梯度。这里应变的梯度是非局部量。确定应变梯度的非局部积分如下：将应变 ij 在一点附近 Taylor 展开：ij x ij x ij,mm 0 2式中 为以 x 为坐标原点的局部坐标。在包含 x 的表示体元内积分上式：ij x k dV ij x k dV ij,m k mdV(20-21)(20-22)VcellVcellVcell假定特征尺寸l 足够小，可忽略l 的高阶小量，于是梯度ij,k 可以表示为应变的积分形式：1 x x

33、 dV dV (20-23)ij,kijijmkmVcellVcell从而由关系式（20-19）、（21-20）可以得出应变梯度 的值。形变理论的本构关系：kk 3Kkk20-17(20-24)f 2 l2 ref(20-25)ijij3式中的l 为由（20-16）式给出的材料长度， 为非局部变量，由（20-20）式给出。理论的本构关系： kk 3K kk(20-26)36 2 lijkl klref 2 (20-27)ijij2 fp fp23eeref加载时 1，卸载时 0 。20.2.3 ABAQUS 用户单元的使用如上文所述，由于 MSG 理论和理论本身的复杂性，用它们做解比较。只是对

34、于很少几个简单问题才有解。为比较理论与实验以及解的符合程度，必须用有限元计算加以验证。因为考虑应变梯度的本构关系与经典理论完全不同，所以在有限元实现中不能利用现有的程序。但是由于 ABAQUS 具有的用户子程序功能，为我们实现应变梯度的有限元计算提供了很方便的条件。为了保证用户子程序能够完成计算应变梯度问题的功能，编写时必须遵守 ABAQUS规定的法则。具体包括 INCLUDE、命名约定、重新定义变量、编译和、测试和调试、用户可用及不可用的通道号、中断分析等。编写的关于应变梯度的用户子程序中，用的是 9 节点矩形单元，由于节点变量包括应变梯度和高阶应力的分量，故在编程过程中不能利用原有的几何关

35、系，也不能利用 ABAQUS 提供的前后处理程序，但可以方便地利用它的求解器。在用户子程序中计算出有限元的刚度矩阵和右端项，再利用 ABAQUS 的求解器解线形代数方程，在用户子程序 UEL 中需要保存的变量，如计算出的广义应变、广义应力等，都必须存在 ABAQUS 指定的变量 SVARS 中，下一次调用 UEL 时再从变量 SVARS 中出需要的变量值。关于用户子程序的调试可通过读写输出文件来获得信息，子程序中可以将调试信息写在 ABAQUS 的消息文件（.msg）（通道号为 7）或数据文件（.dat）（通道号为 6）中，读这些文件从而验证程序是否已满足了计算的要求。到调试信息，20-182

36、0.2.4 有限元计算的结果MSG 理论的有限元计算结果编写了应变梯度本构的平面问题和轴对称问题的用户子程序。利用轴对称程序计算了微小尺度下的浅压头压痕问题，并对微尺度压痕的试验结果进行了参数拟合；对于平面问题的程序，进行了静态裂纹的分析，分析过程主要针对 I 型裂纹。图 20-4 MSG 理论拟合微压痕试验结果图 20-4 是利用 MSG 理论计算微压痕的曲线：坐标横轴是压痕深度的倒数，坐标纵轴是无量的硬度的平方，其中 H 是微压痕的硬度，H0 是 h 取值很大时的压痕硬度，它与压痕的深度 h 无关。图中的试验点是 McElhaney 等人在 1998 年对多晶铜所作的结果。可见，利用 MS

37、G 理论计算的结果曲线在从 1/10 个微米到几个微米很大的范围内与试验结果符合的非常好，计算拟合的材料特征长度 l =1.52 也符合多晶铜由（20-16）式的计算结果以及试验对铜的估算值，这表明 MSG 理论能够比较准确地反映微米到米量级材料的塑性行为。下面利用 MSG 本构的平面问题程序计算 I 型裂纹的结果。整个计算区域为 1000l，远场施加的是弹性位移边界条件，外加 K 场强度为 20。图 20-5 给出了双对数坐标下裂尖的奇异性曲线，即裂纹延长线上（实际上取的是最靠近裂纹面的一列20-19点， 1.0143 ）的等效应力对应距裂尖距离的曲线。横坐标是无量的点到裂尖的距离，纵坐标是

38、无量的等效应力。图中给出了利用 MSG 形变理论和 MSG理论两种结果，同时也给出了同一问题经典塑性理论的计算结果。容易看出 MSG 形变理论和理论的计算结果相差的非常小，这表明，虽然裂纹尖端场不是严格的比例加载，但是利用 MSG 形变理论的计算结果还是足够准确的。图中也可以看出，三种理论计算的裂尖塑性区大小完全相同，说明塑性区大小与应变梯度效应无关；MSG 理论给出的曲线在约 0.3l 处与经典塑性理论分开，MSG 理论的结果比经典理论高很多，在 0.1l 处（对应于铜约是 0.4 微米）等效应力是经典理论的两倍还多。这在某种程度上可以解释裂尖的断裂问题。1011/2KI= 20Yl0.10.010.11101001000r /l图 20-5M SG 理论计算 I 型裂纹尖端场的奇异性曲线理论的有限元计算结果对于理论同样编写了平面问题和轴对称问题的用户子程序。利用轴对称程序计算了微尺度下的浅压头压痕问题，利用平面问题程序