邓肯-张模型的关键点是材料的弹性模量随大小主应力差及小主应力(围压)的变化而变化， 用 APDL 实现之的基本思路是：给每个单元定义一个材料号，分级施加荷载，在每个荷载步 结束时提取出各单元的大小主应力，据此计算出下个荷载步的弹性模量 Et,修改各单元之 MP，用于下一步计算。 以下是一个简单算例，copy 出去可直接运行。 !常规三轴试验模拟 !by taomingxing,NWPU !2003.7.16 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 Sf=2*(c*cos(Fai)+P3*sin(Fai)/(1-sin(Fai) !Mohr-Coulomb 破坏强度(p1-p3)f Ds=(P1-P3)/Sf !应力水平， *if,Ds,GT,0.95,then Ds=0.95 !应力水平最大取值 *endif !判断加卸荷，如果(P1-P3)小于历史最大值视为卸荷再加荷过程 *if,MaxPs(i),LT,P1-P3,then Ei=k*Pa*(P3/Pa)*n Et=Ei*(1-Rf*Ds)*2 !加荷情况的切线模量 MaxPs(i)=P1-P3 !保存历史最大应力 *elseif,MaxPs(i),GE,P1-P3 Et=Kur*Pa*(P3/Pa)*n !卸荷模量 *endif mp,ex,i,Et !修改单元 i 的 Et mp,nuxy,i,Mu *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 !8cm X 15cm /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 !施加围压 sfl,all,pres,2e5 !200kPa solve !分级施加荷载，实现非线性计算 !荷载增量 10kPa,共 50 级 *DO,ti,1,50 !取出计算结果，修改弹性模量 /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 !修改砂土单元的 Et,单元号 1-120 c=0 $Fai=35 $Rf=0.7 $k=400 $n=0.6 $Mu=0.33 $Kur=326.7 *do,i,1,120 !各单元循环计算 *use,Duncan-Chang,c,Fai,Rf,k,n,Mu,Kur !调用 Duncan-Chang 宏文件 *enddo !EEt(ti)=ET !保存第 120 单元之 ET /SOLU time,ti sfl,3,pres,2e5+1e4*ti !施加荷载，增量 1e4 solve !对 ti 级荷载情况求解 *ENDDO