ANSYS邓肯张材料模型

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