ANSYS邓肯-张材料模型

楼主给的在ANSYS上实现邓肯－张模型的方法很有用，但其中还有几点需要修正的，这也是楼上的兄弟们有疑问的原因。我把楼主的代码运行了一下，然后对照作了修改，现在上传一下，有问题的兄弟可以仔细对照一下，在这里我对其中几个比较明显的问题说明一下：
1.MP命令不能直接给单元加材料，这是对的。在这里，楼主遗漏了一下命令：MPCHG，具体见下面的修改过的代码。
2.关于密度的问题。这些要在宏中定义，每修改一种材料(即调用一次邓肯－张子程序)就要修改一次材料的密度，其他有关材料的问题可以类推。
3.关于施加重力的问题。要在调用宏后，在同一个循环中重新定义一下重力。
以下是我修改过的楼主的代码，希望对兄弟们有所帮助。
!用APDL得到初步成果，贴于此供感兴趣的朋友参考，不当之处敬请指正，
!欢迎加以完善。
!基本思路：
!邓肯-张模型的关键点是材料的弹性模量随大小主应力差
!及小主应力(围压)的变化而变化，用APDL实现之的基本思路是：
!给每个单元定义一个材料号，分级施加荷载，在每个荷载步结束时提取出各
!单元的大小主应力，据此计算出下个荷载步的弹性模量Et,修改各单元之MP，
!用于下一步计算。
!以下是一个简单算例，copy出去可直接运行。
!!!常规三轴试验模拟
!**********************************************************
FINISH
/CLEAR
/TITLE,Numerical Simulation of three axes testing of soils
/PREP7
*dim,SUy,array,50 !Settlement records
*dim,MaxPs,array,120 !Max history p1-p3
*dim,MaxDs,array,120 !Max history Ds
!*dim,EEt,array,50 !Et of elememt
!!!Duncan-Chang Model
!!!Symbols：c-粘滞力，Fai-内摩擦角，Sf-破坏强度(p1-p3)f，
!Ds-应力水平，Pa-大气压,P3-围压
!**********************************************************************
*CREATE,Duncan-Chang !Creat Macro file
*afun,deg !Unit of angle
*set,Pa,1e5
*set,P1,-ArrS3(i) !注意：岩土工程中应力为拉负压正
*set,P3,-ArrS1(i)
*if,P3,LT,0.1*Pa,then
P3=0.1*Pa !围压最小取值
*endif
Sf0=2*(c0*cos(Fai)+P3*sin(Fai))/(1-sin(Fai)) !Mohr-Coulomb破坏强度(p1-p3)f
Ds=(P1-P3)/Sf0 !应力水平，
*if,Ds,GT,0.95,then
Ds=0.95 !应力水平最大取值
*endif
!判断加卸荷，如果(P1-P3)小于历史最大值视为卸荷－再加荷过程
*if,MaxPs(i),LT,P1-P3,then
Ei=k0*Pa*(P3/Pa)**n0
Et0=Ei*(1-Rf*Ds)**2 !加荷情况的切线模量
MaxPs(i)=P1-P3 !保存历史最大应力
*elseif,MaxPs(i),GE,P1-P3
Et0=Kur*Pa*(P3/Pa)**n0 !卸荷模量
*endif
mp,ex,i,Et0 !修改材料i的Et
mp,nuxy,i,Mu0
mp,dens,I,1800 ！重新定义第I种材料的密度
*END
!***********************单元类型***********************************
et,1,42

!平面四节点单元
KEYOPT,1,3,2 !平面应变
!*********************以下定义材料初始模量*************************
mp,ex,1,3.728e7 !砂土的弹性模量
mp,nuxy,1,0.33
mp,dens,1,1800

!******************************************************************
blc4,0,0,0.08,0.15 !建立几何模型
/PNUM,AREA,1
/REPLOT
aesize,all,0.01 !网格划分
mat,1
amesh,all
nsel,s,loc,y,0 !边界条件
d,all,Uy !底边界竖向约束
nsel,s,loc,x,0
d,all,Ux !左侧边界水平向约束
nsel,all
/replot
fini
/SOLU
time,0.01 !施加围压
ACEL,0,9.8,0
sfl,all,pres,2e5 !200kPa
solve
!********分级施加荷载，实现非线性计算,荷载增量10kPa,共50级********
*DO,ti,1,2 !取出计算结果，修改弹性模量
/POST1
*get,SUy(ti),node,29,u,y !Settlement record of time ti
ETABLE,EtabS1,S,1 !取各单元第一主应力
ETABLE,EtabS3,S,3 !取各单元第三主应力
*dim,ArrS1,array,120
*dim,ArrS3,array,120
*do,Num,1,120 !Num为单元编号
*get,ArrS1(Num),elem,Num,etab,EtabS1 !将单元结果存入数组
*get,ArrS3(Num),elem,Num,etab,EtabS3
*enddo
/PREP7
c0=0
Fai=35
Rf=0.7
k0=400
n0=0.6
Mu0=0.33
Kur=326.7 !修改砂土单元的Et,单元号1-120
*do,i,1,120 !各单元循环计算
*use,Duncan-Chang,c0,Fai,Rf,k0,n0,Mu0,Kur !调用Duncan-Chang宏文件
mpchg,i,i
!！！！！！注意，这个命令是把第I种材料施加给第I个单元
*enddo
/SOLU !EEt(ti)=ET !保存第120单元之ET
time,ti　
sfl,3,pres,2e5+1e4*ti !施加荷载，增量1e4
ACEL,0,9.8,0 ！重新施加重力
solve !对ti级荷载情况求解
*ENDDO