abaqus2用户单元子程序.docx

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

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

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

4、度，索的张力为： (20-3)在公式20-3中，A为截面面积，E为弹性模量，N0为初始张力。图20-1 索单元在总体坐标系下，单元刚度矩阵为： (20-4)单元刚度矩阵中的子阵K分别由线性和非线性矩阵项组成： (20-5)在公式20-5中的KL和KNL均是33的对称矩阵，分别为： 索单元的节点质量为： (20-6)在公式20-6中，为密度。索单元的质量矩阵为： (20-7)结构的运动方程为： (20-8)公式21-8中Fext为作用在结构上的外力。在不断变化的索的变形中，求解运动方程，得到节点的位移值。20.1.3 应用举例图19-2 由五个单元组成的两端铰接的索杆结构由5个单元组成的两端铰接

5、的索杆结构，高5m长10m，6个节点号码依次为101106，如图19-2所示。计算自由振动的频率和周期。输入文件中的用户单元界面ABAQUS输入文件.inp中的用户单元界面如下：*HEADINGTwo dimensional overhead hoist frame using 2 nodes self-developed truss element,Initial force N is defined in property(5) and referenced by user elementSI Units1-axis horizontal, 2-axis vertical*USER ELE

6、MENT, NODES=2, TYPE=U1, PROPERTIES=5, COORDINATES=3, VARIABLES=121,2,3*UEL PROPERTY, ELSET=UTRUSS*ELEMENT, TYPE=U1, ELSET=UTRUSS计算结果和比拟表20-1列出了由用户索单元计算的图20-2所示结构的固有周期，并与应用ABAQUS梁单元B31的计算结果进行了比拟。索单元与梁单元前4阶模态的周期根本一致；索单元的第69阶模态与梁单元第710阶模态的周期根本一致。从第11阶模态开始，随着梁单元弯曲变形的增加，梁的弯曲刚度逐渐发挥作用并和轴向刚度耦合，与同阶模态的索单元相比，梁

7、单元的振动周期显著降低，而频率高于索单元。表20-1 ABAQUS用户索单元和梁单元B31计算的频率比拟振动模态索单元固有周期(Cycles/time)梁单元固有周期(Cycles/time)123456789101112用户开发单元的缺点是不能采用ABAQUS的后处理进行显示，只能从数据文件.dat中读取结果。另外，ABAQUS的接触算法等某些功能也无法应用。20. 非线性索单元用户子程序 subroutine uel(rhs, amatrx, svars, energy, ndofel, nrhs, nsvars, * props, nprops, coords, mcrd, nnode,

8、 u,du,v,a,jtype,time,dtime, * kstep,kinc,jelem,params,ndload,jdltyp,adlmag,predef,npredf, * lflags,mlvarx,ddlmag,mdload,pnewdt,jprops,njprop,period)C Include aba_param.incC All coordinates in global C dimension rhs(mlvarx,*), amatrx(ndofel, ndofel), * svars(12), energy(8), props(5), coords(mcrd,nnod

9、e), * u(ndofel), du(mlvarx,*), v(ndofel), a(ndofel), time(2), * params(3),jdltyp(mdload,*), adlmag(mdload,*), * ddlmag(mdload,*), predef(2,npredf,nnode),lflags(*), * jprops(*) C dimension sresid(6), uji(3), xji(3), smatrx(3,3)CC Material properties area=props(1) e =props(2) anu =props(3) rho =props(

10、4)C Initial tension force in user element fn0 =props(5)CC Geometry, stiffness and mass parameters dx =coords(1,2)-coords(1,1) dy =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/alen am =0.5d0*area*rho*alenC do i=1,3 uji(i)=u(i+3)-u(i) xji

11、(i)=coords(i,2)-coords(i,1) end do strain=(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 eal=e*area/alen*3CC Stiffness matrix parameters smatrx(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

12、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/alenC do 6 k1=1,ndofel do 2 krhs=1,nrhs2 continue do 4 k2=1,ndofel4 continue6 continueC if (lflags(3).eq.1) thenC Normal incrementati

13、on if (lflags(1).eq.1.or.lflags(1).eq.2) thenC *StaticC Element stiffness matrix do i=1,3 do 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,j) end do end doCC Reaction force if (lflags(4).ne.0) then do i=1,3 force =ak*(u(i+3)-u(i) dforce=

14、ak*(du(i+3,1)-du(i,1) sresid(i)=-dforce sresid(i+3)=dforce rhs(i,1)=rhs(i,1)-sresid(i) rhs(i+3,1)=rhs(i+3,1)-sresid(i+3) end do else do k=1,3 force =ak*(u(k+3)-u(k) sresid(k)=-force sresid(k+3)=force rhs(k,1)=rhs(k,1)-sresid(k) rhs(k+3,1)=rhs(k+3,1)-sresid(k+3) end do end ifC else if(lflags(1).eq.11.or.lflags(1).eq.12) then C *Dynamic alpha =params(1) beta =params(2) gamma =params(3)C dadu =1.0d0/(beta*dtime*2) dvdu =gamma/(beta*dtime)C do 14 k1=1,ndofel amatrx(k1,k1)=am*dadu rhs(k1,1)=rhs(k1,1)-am*a(k1)14 continue do i=1,3 do j=1,3 amatrx(i,j)=amatrx(i,i)+(1.0d0 + alpha)*smatrx(i,j) ama