abaqus2用户单元子程序.docx
《abaqus2用户单元子程序.docx》由会员分享，可在线阅读，更多相关《abaqus2用户单元子程序.docx（24页珍藏版）》请在冰豆网上搜索。
abaqus2用户单元子程序
20ABAQUS用户单元子程序（UEL）
在这一章中将列举两个在这些年里发展过的ABAQUS/Standard用户单元子程序（UEL）。
第一个例子是一个非线性的索单元，我们的目的是通过这个比较简单的例子让读者了解用户单元子程序的基本开发过程；第二个例子是一个用于计算应变梯度理论的单元，应变梯度是当今比较热点的一个科研前沿问题，有各种理论，我们为了验证新的理论，需要数值结果与实验对照来进行评价，整个例子的目的是通过它说明用户子单元可以求解的问题范围很广，但是由于内容比较艰深，程序也很长，所以这个例子我们并没有给出最后的全部程序。
另外，到目前为止，ABAQUS还只有隐式求解器ABAQUS/Standard支持用户自定义单元，而显式求解器ABAQUS/Explicit中还不支持这一功能。
非线性索单元
20.1.1背景
钢索斜拉桥和斜拉索结构广泛应用于土木工程建筑上。
索力的计算分析是设计和施工的关键环节。
清华大学工程力学系在采用ABAQUS进行荆沙长江斜拉桥的计算机仿真分析（这个项目我们已在第15章“ABAQUS在土木工程中的应用
（一）——荆州长江大桥南汊斜拉桥结构三维仿真分析”中讨论过）时，也曾进行了自行建立索单元的尝试。
本节介绍的就是这方面的工作。
香港理工大学土木与结构工程系采用ABAQUS有限元软件进行计算，完成了香港TingKau斜拉桥和TsingMa悬索桥的结构计算和分析。
对于钢索计算，他们采用梁单元进行模拟。
由于梁单元含有弯曲刚度，计算的高阶频率值偏高，周期较低。
一般假设索是单向受拉力的构件。
随着应变的非线性增加，索力呈非线性增加。
尽管ABAQUS单元库中有500个以上的单元类型，但是，还没有索单元。
本文发展了三维非线性索单元模型，形成ABAQUS的用户单元子程序，可以利用ABAQUS输入文件调入到具体的分析中。
通过静态和动态例题的计算比较，索单元工作良好。
20.1.2基本公式
在三维索单元计算中，如图20-1所示，坐标x和位移u的变量表达式为：
（x,y,z）（u,v,w）（20-1）
应变的公式为：
（20-2）
公式（20-2）中，L为索的长度，索的张力为：
（20-3）
在公式（20-3）中，A为截面面积，E为弹性模量，N0为初始张力。
图20-1索单元
在总体坐标系下，单元刚度矩阵为：
（20-4）
单元刚度矩阵中的子阵K分别由线性和非线性矩阵项组成：
（20-5）
在公式（20-5）中的KL和KNL均是3×3的对称矩阵，分别为：
索单元的节点质量为：
（20-6）
在公式（20-6）中，
为密度。
索单元的质量矩阵为：
（20-7）
结构的运动方程为：
（20-8）
公式（21-8）中Fext为作用在结构上的外力。
在不断变化的索的变形中，求解运动方程，得到节点的位移值。
20.1.3应用举例
图19-2由五个单元组成的两端铰接的索杆结构
由5个单元组成的两端铰接的索杆结构，高5m长10m，6个节点号码依次为101~106，如图19-2所示。
计算自由振动的频率和周期。
输入文件中的用户单元界面
ABAQUS输入文件（.inp）中的用户单元界面如下：
*HEADING
Twodimensionaloverheadhoistframe
using2nodesself-developedtrusselement,
InitialforceNisdefinedinproperty（5）and
referencedbyuserelement
SIUnits
1-axishorizontal,2-axisvertical
……
*USERELEMENT,NODES=2,TYPE=U1,PROPERTIES=5,COORDINATES=3,VARIABLES=12
1,2,3
*UELPROPERTY,ELSET=UTRUSS
,,7800,
*ELEMENT,TYPE=U1,ELSET=UTRUSS
……
计算结果和比较
表20-1列出了由用户索单元计算的图20-2所示结构的固有周期，并与应用ABAQUS梁单元B31的计算结果进行了比较。
索单元与梁单元前4阶模态的周期基本一致；索单元的第6~9阶模态与梁单元第7~10阶模态的周期基本一致。
从第11阶模态开始，随着梁单元弯曲变形的增加，梁的弯曲刚度逐渐发挥作用并和轴向刚度耦合，与同阶模态的索单元相比，梁单元的振动周期显著降低，而频率高于索单元。
表20-1ABAQUS用户索单元和梁单元B31计算的频率比较
振动模态
索单元固有周期
（Cycles/time）
梁单元固有周期
（Cycles/time）
1
2
3
4
5
6
7
8
9
10
11
12
用户开发单元的缺点是不能采用ABAQUS的后处理进行显示，只能从数据文件（.dat）中读取结果。
另外，ABAQUS的接触算法等某些功能也无法应用。
20.1.4非线性索单元用户子程序
subroutineuel（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,mdload,pnewdt,jprops,njprop,period）
C
Include''
CAllcoordinatesinglobal
C
dimensionrhs（mlvarx,*）,amatrx（ndofel,ndofel）,
*svars（12）,energy（8）,props（5）,coords（mcrd,nnode）,
*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
dimensionsresid（6）,uji（3）,xji（3）,smatrx（3,3）
C
CMaterialproperties
area=props
（1）
e=props
（2）
anu=props（3）
rho=props（4）
CInitialtensionforceinuserelement
fn0=props（5）
C
CGeometry,stiffnessandmassparameters
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=*area*rho*alen
C
doi=1,3
uji（i）=u（i+3）-u（i）
xji（i）=coords（i,2）-coords（i,1）
enddo
strain=（xji
（1）*uji
（1）+xji
（2）*uji
（2）+xji（3）*uji（3）
*+*（uji
（1）**2+uji
（2）**2+uji（3）**2））/alen
tforce=e*area*strain+fn0
eal=e*area/alen**3
C
CStiffnessmatrixparameters
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）
smatrx（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/alen
C
do6k1=1,ndofel
sresid（k1）=
do2krhs=1,nrhs
rhs（k1,krhs）=
2continue
do4k2=1,ndofel
amatrx（k2,k1）=
4continue
6continue
C
if（lflags（3）.then
CNormalincrementation
if（lflags
（1）.then
C*Static
CElementstiffnessmatrix
doi=1,3
doj=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）
enddo
enddo
C
CReactionforce
if（lflags（4）.then
doi=1,3
force=ak*（u（i+3）-u（i））
dforce=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）
enddo
else
dok=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）
enddo
endif
C
elseif（lflags
（1）.then
C*Dynamic
alpha=params
（1）
beta=params
（2）
gamma=params（3）
C
dadu=（beta*dtime**2）
dvdu=gamma/（beta*dtime）
C
do14k1=1,ndofel
amatrx（k1,k1）=am*dadu
rhs（k1,1）=rhs（k1,1）-am*a（k1）
14continue
doi=1,3
doj=1,3
amatrx（i,j）=amatrx（i,i）++alpha）*smatrx（i,j）
amatrx（i+3,j+3）=amatrx（i+3,j+3）++alpha）*smatrx（i,j）
amatrx（i,j+3）=amatrx（i,j+3）-+alpha）*smatrx（i,j）
amatrx（i+3,j）=amatrx（i+3,j）-+alpha）*smatrx（i,j）
enddo
enddo
C
doi=1,3
force=ak*（u（i+3）-u（i））
sresid（i）=-force
sresid（i+3）=force
rhs（i,1）=rhs（i,1）-（+alpha）*sresid（i）
*-alpha*svars（i））
rhs（i+3,1）=rhs（i+3,1）-（+alpha）*sresid（i+3）
*-alpha*svars（i+3））
enddo
C
do16k1=1,ndofel
svars（k1+6）=svars（k1）
svars（k1）=sresid（k1）
16continue
endif
C
elseif（lflags（3）.then
CStiffnessmatrix
doi=1,3
doj=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）
enddo
enddo
C
elseif（lflags（3）.then
CMassmatrix
do40k1=1,ndofel
amatrx（k1,k1）=am
40continue
elseif（lflags（3）.then
C
CHalf-stepresidualcalculation
alpha=params
（1）
doi=1,3
force=ak*（u（i+3）-u（i））
sresid（i）=-force
sresid（i+3）=force
rhs（i,1）=rhs（i,1）-am*a（i）-+alpha）*sresid（i）
*+*alpha*（svars（i）+svars（i+6））
rhs（i+3,1）=rhs（i+3,1）-am*a（i+3）-+alpha）*sresid（i+3）
*+*alpha*（svars（i+3）+svars（i+9））
enddo
C
elseif（lflags（3）.then
CInitialaccelerationcalculation
do60k1=1,ndofel
amatrx（k1,k1）=am
60continue
doi=1,3
force=ak*（u（i+3）-u（i））
sresid（i）=-force
sresid（i+3）=force
rhs（i,1）=rhs（i,1）-sresid（i）
rhs（i+3,1）=rhs（i+3,1）-sresid（i+3）
enddo
C
do62k1=1,ndofel
svars（k1）=sresid（k1）
62continue
C
elseif（lflags（3）.then
COutputforperturbations
if（lflags
（1）.then
C*static
doi=1,3
force=ak*（u（i+3）-u（i））
dforce=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）
enddo
C
dokvar=1,nsvars
svars（kvar）=
enddo
doj=1,6
svars（j）=rhs（j,1）
enddo
elseif（lflags
（1）.then
C*Frequency
do90krhs=1,nrhs
doi=1,3
dforce=ak*（du（i+3,krhs）-du（i,krhs））
sresid（i）=-dforce
sresid（i+3）=dforce
rhs（i,krhs）=rhs（i,krhs）-sresid（i）
rhs（i+3,krhs）=rhs（i+3,krhs）-sresid（i+3）
enddo
90continue
dokvar=1,nsvars
svars（kvar）=
enddo
doj=1,6
svars（j）=rhs（j,1）
enddo
C
endif
endif
C
return
end
利用ABAQUS用户单元计算应变梯度塑性问题
我们在本节内容中主要介绍两种应变梯度理论，并在最后给出用这两种应变梯度理论编写的用户单元子程序的数值计算结果与实验的对比。
本节的目的在于让读者了解ABAQUS即使在面对如此复杂的理论问题时，也可以胜任。
20.2.1引言
研究应变梯度理论的意义
很多试验表明，当非均匀塑性变形特征长度在微米量级时，材料具有很强的尺度效应。
例如Fleck等在细铜丝的扭转试验中观察到，当铜丝的直径为12m时，无量纲的扭转硬化将增加至170m直径时的3倍；Stolken和Evans在薄梁弯曲试验中也观察到，当梁的厚度从100m减至12.5m时，无量纲的弯曲硬化也显著增加；而在单轴拉伸情况这种尺度效应并不存在。
在微米量级的尺度下，微观硬度试验与颗粒增强金属基复合材料中也观察到尺度效应，当压痕深度从10m减至1m时，金属的硬度增加了一倍；对于碳化硅颗粒加强的铝－硅基复合材料，Lloyd观察到当保持颗粒体积比为15%的条件下，将颗粒直径从16m减为7.5m后复合材料的强度显著增加。
由于在传统的塑性理论中，本构模型不包含任何尺度，所以它不能预测尺度效应。
然而，在工程实践中迫切需要处理微米量级的设计和制造问题，例如，厚度在1m或者更小尺寸下的薄膜；整个系统尺寸不超过10m的传感器、执行器和微电力系统（MEMS）；零部件尺寸小于10m的微电子封装；颗粒或者纤维的尺寸在微米量级的先进复合材料及微加工等等。
现在的设计方法，如有限元方法（FEM）和计算机辅助设计（CAD），都是基于经典的塑性理论，而它们在这一微小尺度不再适用。
另一方面，现在按照量子力学和原子模拟的方法在现实的时间和长度的尺度下处理微米尺度的结构依然很困难。
所以，建立连续介质框架下、考虑尺寸效应的本构模型就成为联系经典塑性力学和原子模拟之间必要的桥梁。
促使建立细观尺寸下连续介质理论的另一个目的是在韧性材料的宏观断裂行为和原子断裂过程之间建立联系。
在一系列值得注意的试验中，Elssner等测量了单晶铌/蓝宝石界面的宏观断裂韧度和原子分离功，使原子点阵或强界面分离所需要的力约为或者10Y（E为弹性模量，Y为拉伸屈服应力）。
而按照经典的塑性理论，Hutchinson指出裂纹前方最大应力水平只能达到4至5倍Y。
很明显这远远小于Elssner等在试验中观察到的结果，不足使原子分离。
考虑应变梯度的影响有望解释这一现象。
应变梯度理论简介
目前发展的应变梯度理论有很多种，包括CS理论（偶应力理论）、SG理论（拉伸和旋转梯度理论）、MSG理论（基于细观机制的应变梯度塑性理论）以及TNT理论（基于Taylor关系的非局部应变梯度理论）等。
我们利用ABAQUS用户单元主要进行了MSG和TNT两种理论应变梯度塑性的有限元分析，对于MSG塑性和TNT塑性都包括形变理论和流动理论，TNT塑性还包括了有限变形问题的形变和流动理论的分析。
现在对MSG理论和TNT理论加以简单的介绍。
20.2.2两种应变梯度理论
基于细观机制的MSG应变梯度塑性理论
基于位错机制的MSG应变梯度塑性理论是由位错理论出发的，它通过一个多尺度、分层次的框架由微观位错机制推导出了宏细观的应变梯度塑性理论。
这个理论相比于其它理论，物理机制更明确，构造方法系统，而且第一次提出了材料长度的表达式。
图20-3是MSG应变梯度塑性理论的原理图，在微观层次上塑性是由位错运动产生的，在细观层次上引入应变的梯度与微观的几何必须位错密度相关联，通过细观和微观的功等效由微观塑性推导出细观的本构理论。
为了在细观尺度下的应变梯度塑性和微尺度下的Taylor硬化关系之间建立联系，在MSG理论框架中采用如下的基本假设：
图20-3MSG理论中采用的多尺度框架
1）假设微尺度的流动应力由位错运动控制，并且遵守应变梯度律给出的Taylor硬化关系
（20-9）
2）微观尺度和细观尺度的联系是塑性功相等
（20-10）
3）在微尺度胞元中假设经典塑性的基本结构成立，其J2形变理论可以表示为
（20-11）
其中
，
。
微尺度的屈服条件为
（20-12）
基于以上的理论假设，应变梯度塑性MSG形变理论本构关系可以建立如下：
（20-13）
（20-14）
其中
（20-15）
式中：
（20-16）
为材料特征长度。
（20-17）
（20-18）
其中应变梯度表示为
（20-19）
（20-20）
基于Taylor关系的非局部应变梯度塑性理论（TNT理论）
MSG应变梯度理论物理机制明确、构造方法系统，那么为什么还要发展TNT理论呢因为前面提到的应变梯度塑性理论，无论CS、SG还是MSG，都是高阶理论，引入了高阶应力和附加的边界条件，而这些高阶应力和附加的边界条件都难以测量，难以想象，因此无法得到工业界的认可，难以走向实用。
经典的塑性问题控制方程是2阶，而在MSG理论中由于高阶应力的影响，其控制方程是4阶，这就增加了解决问题的复杂性。
非局部连续介质力学理论给我们以启发，采用应变的非局部加权积分来确定应变的梯度，这样就有了TNT——基于TAYLOR关系的非局部塑性理论。
这种理论既保持MSG理论的所有优点，同时与经典的塑性理论相比又不增加方程的阶数。
但也正是由于TNT理论的非局部性质，使其在求取解析解方面比较困难，也正因为如此，有限元解对TNT理论尤其重要。
TNT作为非局部塑性理论，有三个基本特点：
（1）流动应力遵从TAYLOR硬化关系，这是TNT塑性理论的出发点。
（2）应变梯度和几何必须位错密度是非局部量，表示为应变的加权平均。
这是TNT塑性理论的核心概念。
（3）TNT塑性理论保持经典塑性理论的基本结构。
这个特点使得TNT塑性理论具有很强的实用潜力。
和经典的塑性理论相比，TNT塑性理论的特别之处在于屈服条件的不同：
流动应力不仅依赖于应变，还同时依赖于应变的梯度。
这里应变的梯度是非局部量。
确定应变梯度的非局部积分如下：
将应变
在一点附近Taylor展开：
（20-21）
式中
为以x为坐标原点的局部坐标。
在包含x的表示体元内积分上式：
（20-22）
假定特征尺寸
足够小，可忽略
的高阶小量，于是梯度
可以表示为应变的积分形式：
（20-23）
从而由关系式（20-19）、（21-20）可以得出应变梯度
的值。
TNT形变理论的本构关系：
（20-24）
（20-25）
式中的
为由（20-16）式给出的材料长度，
为非局部变量，由（20-20）
升级会员 