abaqus2用户单元子程序

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

另外，到目前为止，ABAQUS还只有隐式求解器ABAQUS/Standard支持用户自定义单元，而显式求解器ABAQUS/Explicit中还不支持这一功能。

20.1 非线性索单元

20.1.1 背景

钢索斜拉桥和斜拉索结构广泛应用于土木工程建筑上。索力的计算分析是设计和施工的关键环节。清华大学工程力学系在采用ABAQUS进行荆沙长江斜拉桥的计算机仿真分析（这个项目我们已在第15章“ABAQUS在土木工程中的应用（一）——荆州长江大桥南汊斜拉桥结构三维仿真分析”中讨论过）时，也曾进行了自行建立索单元的尝试。本节介绍的就是这方面的工作。

香港理工大学土木与结构工程系采用ABAQUS有限元软件进行计算，完成了香港Ting Kau斜拉桥和Tsing Ma悬索桥的结构计算和分析。对于钢索计算，他们采用梁单元进行模拟。由于梁单元含有弯曲刚度，计算的高阶频率值偏高，周期较低。

一般假设索是单向受拉力的构件。随着应变的非线性增加，索力呈非线性增加。尽管ABAQUS单元库中有500个以上的单元类型，但是，还没有索单元。本文发展了三维非线性索单元模型，形成ABAQUS的用户单元子程序，可以利用ABAQUS输入文件调入到具体的分析中。通过静态和动态例题的计算比较，索单元工作良好。

20-1

20.1.2 基本公式

在三维索单元计算中，如图20-1所示，坐标x和位移u的变量表达式为：

xjixjxiu (u,v,w) jiuju (x,y,z) i应变的公式为：

1Lxjiujiyjiv1jizjiwji2u222jivjiwji 公式（20-2）中，L为索的长度，索的张力为：

NAEN0

在公式（20-3）中，A为截面面积，E为弹性模量，N0为初始张力。

图20-1 索单元

在总体坐标系下，单元刚度矩阵为：

KeKK66KK 单元刚度矩阵中的子阵K分别由线性和非线性矩阵项组成：

20-2

(20-1)

(20-2)

(20-3)

(20-4)

K33KL,33KNL,33

(20-5)

在公式（20-5）中的KL和KNL均是3×3的对称矩阵，分别为：

x2jiEAKL3Lxjiyjiy2jixjizji100Nyjizji KNL10

L21zji1AL 2索单元的节点质量为：

m(20-6)

在公式（20-6）中，为密度。索单元的质量矩阵为：

结构的运动方程为：

Mem0

0m(20-7)

FextKu Mu(20-8)

公式（21-8）中Fext为作用在结构上的外力。在不断变化的索的变形中，求解运动方程，得到节点的位移值。

20.1.3 应用举例

图19-2 由五个单元组成的两端铰接的索杆结构

20-3

由5个单元组成的两端铰接的索杆结构，高5m长10m，6个节点号码依次为101~106，如图19-2所示。计算自由振动的频率和周期。 输入文件中的用户单元界面

ABAQUS输入文件（.inp）中的用户单元界面如下：

*HEADING

Two dimensional overhead hoist frame using 2 nodes self-developed truss element, Initial force N is defined in property(5) and referenced by user element SI Units

1-axis horizontal, 2-axis vertical ……

*USER ELEMENT, NODES=2, TYPE=U1, PROPERTIES=5, COORDINATES=3, VARIABLES=12 1,2,3

*UEL PROPERTY, ELSET=UTRUSS

1.963E-5, 2.0E11, 0.3, 7800, 10.0E5 *ELEMENT, TYPE=U1, ELSET=UTRUSS

……

计算结果和比较

表20-1列出了由用户索单元计算的图20-2所示结构的固有周期，并与应用ABAQUS梁单元B31的计算结果进行了比较。

索单元与梁单元前4阶模态的周期基本一致；索单元的第6~9阶模态与梁单元第7~10阶模态的周期基本一致。从第11阶模态开始，随着梁单元弯曲变形的增加，梁的弯曲刚度逐渐发挥作用并和轴向刚度耦合，与同阶模态的索单元相比，梁单元的振动周期显著降低，而频率高于索单元。

表20-1 ABAQUS用户索单元和梁单元B31计算的频率比较

20-4

振动模态 1 2 3 4 5 6 7 8 9 10 11 12 索单元固有周期 (Cycles/time) 112.42 112.42 213.83 213.83 249.51 294.32 294.32 345.99 345.99 474.60 653.22 767.91 梁单元固有周期 (Cycles/time) 112.48 112.48 213.95 213.95 222.55 222.75 294.48 294.48 346.19 346.19 422.37 423.69

用户开发单元的缺点是不能采用ABAQUS的后处理进行显示，只能从数据文件（.dat）中读取结果。另外，ABAQUS的接触算法等某些功能也无法应用。

20.1.4 非线性索单元用户子程序

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,mdload,pnewdt,jprops,njprop,period) C

Include 'aba_param.inc' C All coordinates in global C

dimension rhs(mlvarx,*), amatrx(ndofel, ndofel), * svars(12), energy(8), props(5), coords(mcrd,nnode),

20-5

* 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) C

C Material properties area=props(1) e =props(2) anu =props(3) rho =props(4)

C Initial tension force in user element fn0 =props(5) C

C 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*alen C

do i=1,3

uji(i)=u(i+3)-u(i)

xji(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

20-6

tforce=e*area*strain+fn0 eal=e*area/alen**3 C

C 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) 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

do 6 k1=1,ndofel sresid(k1)= 0.0d0 do 2 krhs=1,nrhs rhs(k1,krhs)=0.0d0 2 continue

do 4 k2=1,ndofel amatrx(k2,k1)=0.0d0 4 continue 6 continue C

if (lflags(3).eq.1) then C Normal incrementation

if (lflags(1).eq.1.or.lflags(1).eq.2) then C *Static

C Element stiffness matrix do i=1,3

20-7

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 do C

C Reaction force

if (lflags(4).ne.0) then do i=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) 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 if C

else if(lflags(1).eq.11.or.lflags(1).eq.12) then C *Dynamic

20-8

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)

amatrx(i+3,j+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) amatrx(i+3,j)=amatrx(i+3,j)-(1.0d0 + alpha)*smatrx(i,j) end do end do C

do i=1,3

force =ak*(u(i+3)-u(i)) sresid(i)= -force sresid(i+3)= force

rhs(i,1)=rhs(i,1)-((1.0d0+alpha)*sresid(i) * -alpha*svars(i))

rhs(i+3,1)=rhs(i+3,1)-((1.0d0+alpha)*sresid(i+3) * -alpha*svars(i+3)) end do C

20-9

do 16 k1=1,ndofel svars(k1+6)=svars(k1) svars(k1)=sresid(k1) 16 continue end if C

else if(lflags(3).eq.2) then C 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 do C

else if(lflags(3).eq.4) then C Mass matrix do 40 k1=1,ndofel amatrx(k1,k1)=am 40 continue

else if(lflags(3).eq.5) then C

C Half-step residual calculation alpha = params(1) do i=1,3

force =ak*(u(i+3)-u(i)) sresid(i)= -force sresid(i+3)= force

20-10

rhs(i,1)=rhs(i,1)-am*a(i)-(1.0d0+alpha)*sresid(i) * +0.5d0*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 C

else if(lflags(3).eq.6) then C Initial acceleration calculation do 60 k1=1,ndofel amatrx(k1,k1)=am 60 continue do i=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) end do C

do 62 k1=1,ndofel svars(k1)=sresid(k1) 62 continue C

else if(lflags(3).eq.100) then C Output for perturbations

if (lflags(1).eq.1.or.lflags(1).eq.2) then C *static do i=1,3

force =ak*(u(i+3)-u(i))

dforce =ak*(du(i+3,1)-du(i,1))

20-11

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 C

do kvar=1,nsvars svars(kvar)=0.0d0 end do do j=1,6

svars(j)=rhs(j,1) end do

else if(lflags(1).eq.41) then C *Frequency

do 90 krhs=1,nrhs do i=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) end do

90 continue do kvar=1,nsvars svars(kvar)=0.0d0 end do do j=1,6

svars(j)=rhs(j,1) end do C

20-12

end if end if C

return end

20.2 利用ABAQUS用户单元计算应变梯度塑性问题

我们在本节内容中主要介绍两种应变梯度理论，并在最后给出用这两种应变梯度理论编写的用户单元子程序的数值计算结果与实验的对比。本节的目的在于让读者了解ABAQUS即使在面对如此复杂的理论问题时，也可以胜任。

20.2.1 引言

研究应变梯度理论的意义

很多试验表明，当非均匀塑性变形特征长度在微米量级时，材料具有很强的尺度效应。例如Fleck等在细铜丝的扭转试验中观察到，当铜丝的直径为12m时，无量纲的扭转硬化将增加至170m直径时的3倍；Stolken和Evans在薄梁弯曲试验中也观察到，当梁的厚度从100m减至12.5m时，无量纲的弯曲硬化也显著增加；而在单轴拉伸情况这种尺度效应并不存在。在微米量级的尺度下，微观硬度试验与颗粒增强金属基复合材料中也观察到尺度效应，当压痕深度从10m减至1m时，金属的硬度增加了一倍；对于碳化硅颗粒加强的铝－硅基复合材料，Lloyd观察到当保持颗粒体积比为15%的条件下，将颗粒直径从16m减为7.5m后复合材料的强度显著增加。

由于在传统的塑性理论中，本构模型不包含任何尺度，所以它不能预测尺度效应。然而，在工程实践中迫切需要处理微米量级的设计和制造问题，例如，厚度在1m或者更小尺寸下的薄膜；整个系统尺寸不超过10m的传感器、执行器和微电力系统（MEMS）；零部件尺寸小于10m的微电子封装；颗粒或者纤维的尺寸在微米量级的先进复合材料及微加工等等。现在的设计方法，如有限元方法（FEM）和计算机辅助设计（CAD），都是基于经典的塑性理论，而它们在这一微小尺度不再适用。另一方面，现在按照量子力学和原子模拟的方法在现实的时间和长度的尺度下处理微米尺度的结

20-13

构依然很困难。所以，建立连续介质框架下、考虑尺寸效应的本构模型就成为联系经典塑性力学和原子模拟之间必要的桥梁。

促使建立细观尺寸下连续介质理论的另一个目的是在韧性材料的宏观断裂行为和原子断裂过程之间建立联系。在一系列值得注意的试验中，Elssner等测量了单晶铌/蓝宝石界面的宏观断裂韧度和原子分离功，使原子点阵或强界面分离所需要的力约为0.03E或者10（Y为拉伸屈服应力）。而按照经典的塑性理论，HutchinsonYE为弹性模量，指出裂纹前方最大应力水平只能达到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-14

图20-3 MSG理论中采用的多尺度框架

1）假设微尺度的流动应力由位错运动控制，并且遵守应变梯度律给出的Taylor硬化关系

~Yf2~l

ij(20-9)

2）微观尺度和细观尺度的联系是塑性功相等

VcelldV~ij~ijijijk Vcell ijk(20-10)

3）在微尺度胞元中假设经典塑性的基本结构成立，其J2形变理论可以表示为

~2~~ ij3~ ije32(20-11)

~其中e

~~~ijij，23~ij~ij。微尺度的屈服条件为

~~ e(20-12)

基于以上的理论假设，应变梯度塑性MSG形变理论本构关系可以建立如下：

20-15

其中 

ij2ij3

22ijklijkijkYffijk

(20-13)

(20-14)

Yf

2l

(20-15)

式中： 

为材料特征长度。 

2l18refb2(20-16)

ikppjjkppi72ijk2ijkkjikij142ijkmnikjmnjkimn154ikjpjkippmn4(20-17) (20-18)

其中应变梯度表示为 

ijkuk,ijik,jjk,iij,k

1ijkijk4(20-19)

(20-20)

基于Taylor关系的非局部应变梯度塑性理论（TNT理论）

MSG应变梯度理论物理机制明确、构造方法系统，那么为什么还要发展TNT理论呢？因为前面提到的应变梯度塑性理论，无论CS、SG还是MSG，都是高阶理论，引入了高阶应力和附加的边界条件，而这些高阶应力和附加的边界条件都难以测量，难以想象，因此无法得到工业界的认可，难以走向实用。

经典的塑性问题控制方程是2阶，而在MSG理论中由于高阶应力的影响，其控制方程是4阶，这就增加了解决问题的复杂性。非局部连续介质力学理论给我们以启发，

20-16

采用应变的非局部加权积分来确定应变的梯度，这样就有了TNT —— 基于TAYLOR关系的非局部塑性理论。这种理论既保持MSG理论的所有优点，同时与经典的塑性理论相比又不增加方程的阶数。但也正是由于TNT理论的非局部性质，使其在求取解析解方面比较困难，也正因为如此，有限元解对TNT理论尤其重要。

TNT作为非局部塑性理论，有三个基本特点：

（1）流动应力遵从TAYLOR硬化关系，这是TNT塑性理论的出发点。

（2）应变梯度和几何必须位错密度是非局部量，表示为应变的加权平均。这是TNT塑性理论的核心概念。

（3）TNT塑性理论保持经典塑性理论的基本结构。这个特点使得TNT塑性理论具有很强的实用潜力。和经典的塑性理论相比，TNT塑性理论的特别之处在于屈服条件的不同：流动应力不仅依赖于应变，还同时依赖于应变的梯度。这里应变的梯度是非局部量。

确定应变梯度的非局部积分如下： 将应变ij在一点附近Taylor展开： 

ijxξijxij,mm0ξ

2(20-21)

式中ξ为以x为坐标原点的局部坐标。在包含x的表示体元内积分上式： 

VcellxξijkdVijxkdVij,mVcellVcellkmdV

(20-22)

假定特征尺寸l足够小，可忽略l的高阶小量，于是梯度ij,k可以表示为应变的积分形式：

ij,kijxξijxmdVkmdVVVcellcell从而由关系式（20-19）、（21-20）可以得出应变梯度的值。

1(20-23)

TNT形变理论的本构关系： 

kk3Kkk

(20-24)

20-17



ijTNT流动理论的本构关系：

2reff2l3ij

(20-25)

式中的l为由（20-16）式给出的材料长度，为非局部变量，由（20-20）式给出。



kkkk3K23ij6klklreflij2ij22e3fferefpp(20-26)

(20-27)

加载时1，卸载时0。

20.2.3 ABAQUS用户单元的使用

如上文所述，由于MSG理论和TNT理论本身的复杂性，用它们做解析解比较困难。只是对于很少几个简单问题才有解析解。为比较理论与实验以及解析解的符合程度，必须用有限元计算加以验证。因为考虑应变梯度的本构关系与经典理论完全不同，所以在有限元实现中不能利用现有的程序。但是由于ABAQUS具有的用户子程序功能，为我们实现应变梯度的有限元计算提供了很方便的条件。

为了保证用户子程序能够完成计算应变梯度问题的功能，编写时必须遵守ABAQUS规定的法则。具体包括INCLUDE声明、命名约定、重新定义变量、编译和链接、测试和调试、用户可用及不可用的通道号、中断分析等。我们编写的关于应变梯度的用户子程序中，用的是9节点矩形单元，由于节点变量包括应变梯度和高阶应力的分量，故在编程过程中不能利用原有的几何关系，也不能利用ABAQUS提供的前后处理程序，但可以方便地利用它的求解器。我们在用户子程序中计算出有限元的刚度矩阵和右端项，再利用ABAQUS的求解器解线形代数方程，在用户子程序UEL中需要保存的变量，如计算出的广义应变、广义应力等，都必须存在ABAQUS指定的变量SVARS中，下一次调用UEL时再从变量SVARS中读取出需要的变量值。关于用户子程序的调试可通过读写输出文件来获得信息，子程序中可以将调试信息写在ABAQUS的消息文件（.msg）（通道号为7）或数据文件（.dat）（通道号为6）中，读这些文件可得到调试信息，从而验证程序是否已满足了计算的要求。

20-18

20.2.4 有限元计算的结果

MSG理论的有限元计算结果

我们编写了应变梯度本构的平面问题和轴对称问题的用户子程序。利用轴对称程序计算了微小尺度下的浅压头压痕问题，并对微尺度压痕的试验结果进行了参数拟合；对于平面问题的程序，进行了静态裂纹的分析，分析过程主要针对I型裂纹。

5.04.54.03.53.02 MSG Deformation Theory Experimental Data(H/H0)2.52.01.51.00.50.0012345-16781/h (m)

图20-4 MSG理论拟合微压痕试验结果

图20-4是利用MSG理论计算微压痕的曲线：坐标横轴是压痕深度的倒数，坐标纵轴是无量纲化的硬度的平方，其中H是微压痕的硬度，H0是h取值很大时的压痕硬度，它与压痕的深度h无关。图中的试验点是McElhaney等人在1998年对多晶铜所作的结果。可见，利用MSG理论计算的结果曲线在从1/10个微米到几个微米很大的范围内与试验结果符合的非常好，计算拟合的材料特征长度l =1.52也符合多晶铜由（20-16）式的计算结果以及试验对铜的估算值，这表明MSG理论能够比较准确地反映微米到亚微米量级材料的塑性行为。

下面讨论利用MSG本构的平面问题程序计算I型静止裂纹的结果。整个计算区域为1000l，远场施加的是弹性位移边界条件，外加K场强度为20。图20-5给出了双对数坐标下裂尖的奇异性曲线，即裂纹延长线上（实际上取的是最靠近裂纹面的一列高斯

20-19

点，1.0143）的等效应力对应距裂尖距离的曲线。横坐标是无量纲化的点到裂尖的距离，纵坐标是无量纲化的等效应力。图中给出了利用MSG形变理论和MSG流动理论两种结果，同时也给出了同一问题经典塑性理论的计算结果。容易看出MSG形变理论和流动理论的计算结果相差的非常小，这表明，虽然裂纹尖端场不是严格的比例加载，但是利用MSG形变理论的计算结果还是足够准确的。图中也可以看出，三种理论计算的裂尖塑性区大小完全相同，说明塑性区大小与应变梯度效应无关；MSG理论给出的曲线在约0.3l处与经典塑性理论分开，MSG理论的结果比经典理论高很多，在0.1l处（对应于铜约是0.4微米）等效应力是经典理论的两倍还多。这在某种程度上可以解释裂尖的断裂问题。

10 MSG Flow Theory MSG Deformation Theory Classical Plasticity Theorye/Y1KI= 20Yl0.10.011/20.11101001000r /l图20-5 MSG理论计算I型裂纹尖端场的奇异性曲线

理论的有限元计算结果

对于TNT理论同样我们编写了平面问题和轴对称问题的用户子程序。利用轴对称程序计算了微尺度下的浅压头压痕问题，利用平面问题程序计算了静态裂纹I型问题。

20-20

5.04.54.03.53.02=0.3, l=1.52(H/H0)2.52.01.51.00.50.0012345-1 TNT Deformation Theory McElhaney et al.'s Experiment6781/h (m)图20-6 TNT理论拟合微压痕试验结果

图20-6是应用TNT小变形形变理论预测的多晶铜的压痕硬度：坐标横轴是压痕深度的倒数，坐标纵轴是无量纲化的硬度的平方。从图20-6中可以看出：理论预测的压痕硬度的平方与压痕深度的倒数成线形关系；压痕深度从几微米到零点几微米的范围内，理论预测值与实验符合的都非常好。这说明TNT非局部塑性理论能很好的描述材料在微米和亚微米量级时的机械行为。经典塑性理论由于不包含尺度，所以所预测的压痕硬度与压痕深度无关，在图上是一条水平线。

下面讨论利用TNT本构的平面问题程序计算I型静止裂纹的结果。边界条件和网格划分与MSG理论相同。图20-7给出了双对数坐标下裂尖的奇异性曲线，图中的虚线为经典塑性理论的预测值，实线为TNT理论的预测值。可见TNT曲线在弹性部分和经典塑性理论完全重合，在塑性部分的外面部分也和经典塑性理论一致；从r = 0.3~0.4l开始，TNT塑性理论的预测值开始超过经典塑性理论的预测值。经典塑性理论的双对数曲线的奇异性为N/(N+1)，与HRR场一致，而TNT非局部塑性理论的奇异性曲线的奇异性甚至超过了1/2。在r = 0.1l处TNT塑性理论的预测值比经典塑性理论的预测值高出一倍多，这为解释韧性材料中的解理断裂提供了一种新的解释机理。在TNT塑性理论的曲线与经典理论的曲线分开的r = 0.3~0.4l处，对于铜l = 4~5微米，r约为1.5微米，是位

20-21

错间距的10多倍，连续介质塑性理论仍然是适用的。整体的来看，这条曲线可以把裂纹尖端场区分为4个区：最靠近裂尖的无位错区或稀疏位错区，这个区里连续介质塑性理论不适用；TNT塑性理论控制区；HRR场控制区；和弹性K场。

10 TNT Deformation Theory Classical Plasticity TheoryKI=20Yle/Y11/2III0.10.01III0.11101001000r/l

图20-7 TNT形变理论计算I型裂纹尖端场奇异性曲线

10 TNT Flow Theory MSG Flow Theorye/Y1KI=20Yl0.10.011/20.11101001000r/l

图20-8 TNT流动理论与MSG流动理论奇异性曲线比较

20-22

图20-8是TNT流动塑性理论与MSG流动塑性理论计算裂纹尖端场的结果比较。很明显，两种理论预测的曲线完全重合。

本章主要内容引自：

[1] 黄东平，庄茁，发展ABAQUS用户单元-非线性索单元，2000年ABAQUS用户年会，清华大学，2000

[2] 邱信明，黄克智，利用ABAQUS用户单元计算应变梯度塑性问题，2000年ABAQUS用户年会，清华大学，2000。

20-23