土石坝有限元分析

土石坝有限元分析
1 •问题描述
采用邓肯-张模型对土石坝施工过程和蓄水状态受力情况进行分析，选择通用有限元分 析软件ANSYS 作为研究平台，计算土石坝竣工期竖向方向沉降、水平方向沉降、最大主应 力和最小主应力情况。

坝体结构示意图如图
1所示。

本文主要完成以下工作：采用 ANSYS
内部参数化设计语言 APDL 编写邓肯-张模型计算材料弹性参数；使用中点增量法计算每步施 工单元材料弹性参数；
根据位移修正算法，编写专用程序对计算结果进行处理，
获得坝体沉
降云图。

(本文针对每个步骤提供相应的 APDL 程序，方便后续研究人员进一步研究，也希
望阅读本文的读者能够将自己的研究成果与大家分享。

相关程序可能存在错误， 笔者也未能
表
1各分区材料的物理力学参数
披号 (KN/m' >
i KN/m' )
K
岛
审f
() &
Cz (t/m*)
()
①
1H P 93
400
0.3 <1 5 35
fl.H
2256
12.75
1.W MH tun n.2 0.4S
n
务
(1 7 ③
17,27
w
c.2 n.Sfr i 5
n n.fi
ts.w
1<1,W
1
“
MO K40 0.3 “5
0.4
35
2
H.7
仿真分析结果如下图。

图2竖向沉降位移云图
urvkz. soLnrrcM
QTEZMT
UX
EJLVY3)
M¥5-0
PHM ・・471a sra ■-・iDt43fc
AMK -a M«L2»
-,lj6^3S -・ 7 60 ?7T
-a 01：i51S
-#03524e
.907211 图3水平方向位移云图
NQML aoLurrwa RSr3^0 EKS --21471B SM« ■-.Z144^B X ■・C2£397 AHSYS ReiMae lT.g
Build 1 " . Li
--214426
-.inn
=T™is
-.1J41S3
■■口血
M
-.M712 汕
7
MfSYS R B 1«4?B IJ.O
Euild 1丁・0
rDZ^gil
图4最大主应力云图
f g”
TJTTEti
-113m
YOPIM
-2M0L7
-14A504
图5最小主应力云图
hv- FLF^rrr
STET^-15 3UB 二 1
TIHE-15
tMK
3 帼—
吕
MX -32121,3
Ml^va He L ea se 17 ・Q
E-LLld 17-0
仿真分析流程图如下图。

图6仿真分析流程图
2 •关键仿真分析过程
2.1网格划分与单元组件创建
当几何模型比较规则时，尽可能采用映射方式划分网格，网格分布规则，位移结果过渡光滑一些。

一般情况下几何模型比较复杂，此时建议将截面网格尺寸设置小一些，可以设置为每次浇筑层厚度的四分之一。

采用扫略的方式划分网格，扫略方向可以设置少一些网格，
控制整体网格数量。

有限元网格模型如下图。

图7有限元网格模型
坝体浇筑分为13步完成，每次浇筑层厚度为1m，根据竖向坐标选取浇筑层单元，仓U 建单元组件，如图8和图9所示。

图8创建单元组件
图9单元示意图
相关命令流程序如下: ! 单元分组
================ vsel,s,loc,y,0,13.2 alls,below,volu cm,ebar,elem cm,n bar, node ystep=13.2/13 ytorl=0.2 *do,i,1,13
cmsel,s ,n
bar
cmsel,s,ebar
n sel,r,loc,y,ystep*(i-1)-ytorl,ystep*i+ytorl esl n,r,1 cm,e%i%,elem !
组件名格式为exx *enddo
2.2初始应力场计算
初始应力场计算时，采用生死单元法抑制所有填筑层土体， 仅保留地基土体处于激活状 态。

由于地基部分包含了两种土层：基岩和地表覆盖层。

所以分为两步计算土体初始应力场。

图10地基最小主应力分布云图
相关命令流如下:
! 求解器
!==================================================
选择坝体几何体 选择坝体单元和节点 创建单元组件ebar 创建节点组件nbar
浇筑层厚度 选择重叠区域范围 循环建立每步浇筑层组件
AVG EI.EKRMT SOWTTOH
Y L7
S-M N
Jr STgP-2
SUB =1
TiMK-2 81
(AVS)
LKX -..096951
am --2S7m
EKX =-102 w. a
I'H ?|：I3 h
21:523^6
-26SJJJ
-2J2--L1
-2014 0^ -1577415
/solu
! 第1 步：激活基岩部
分!=================== == antype,0 nropt,full rescontrol,define,all,last
outres,all,last acel,,9.806
cmsel,s,ebar ekill,all
cmsel,s,volu3 alls,below,volu ekill,all
cmsel,s,volu4 alls,below,volu myinis myinis
alls time,1 save solve rstnew parsav,all,parms
! 第2 次计算fini
/solu
alls matnew,emntr,1,1
cmsel,s,volu4 alls,below,volu emcalc
time,1 alls
save,case1_1,db,,model 仅通过最后子步重启动分析保
存在最后子步保存所有结果施
加重力加速度载荷
杀死所有坝体填筑层单元
杀死地表覆盖层土体单元
计算基岩初始应力场
计算
提取应力结果
保存参数信息
根据eb 模型计算土体弹性参数
根据中点增量法计算土体实际参数计算
parsav,all,parms
! 第2 步：激活地表土
层!=================== == fini
/solu
antype,,restart,1,,continue parres,,parms
alls
cmsel,s,ebar
ekill,all
cmsel,s,volu3
alls,below,volu
ealive,all
myinis
time,2
alls
save,case1_2,db,,model solve
rstnew
parsav,all,parms
! 第2 次计算
fini
/solu
antype,,restart,1,,continue parres,,parms
alls matnew,emntr,1,1
cmsel,s,ebar
ekill,all
重启动分析
杀死坝体填筑层
计算地表覆盖层土体初始应力
求解
重启动分析
按照EB 模型计算土体弹性参数
cmsel,s,volu3 alls,below,volu ealive,all
按照中点增量法计算土体弹性参数
time,2
alls
save,case1_2,db,,model solve
parsav,all,parms
2.3 邓肯-张模型(Duncan-Chang EB Model )
邓肯-张模型根据单元的应力状态来评估弹性参数。

当6 -乙”：；「1-匚3 max 且S^s^x
时，单元处于卸荷状态，弹性模量用
E ur ；否则，单元处于加荷状态，弹性模量用 E t 表示。

G -匚3 max 为历史最大偏应力；
s max 为历史最大应力水平。

卸荷或重复加载时的回弹模量采用下式计算:
式中：G 和二3分别为最大和最小主应力； P a 为标准大气压；R f 为破坏比；K ur 为卸
荷时弹性系数；K 为加荷时弹性系数； 二1-；「3 f 为破坏剪应力；匚1-；「3ult 为剪应力极限
值。

切线体积模量采用下式计算:
B t 式中，K b 为体积模量系数； m 为体积模量指数。

esel,s,live
emcalc E ur
加荷时，弹性模量 E t 采用下式计算:
E t = KP a
I Pa 丿
2ccos 「亠 2；「3S in :
R 「上亠
W ult
(1)
(2)
(3)
(4)
=K ur
巳
巳
二1 -二3 f
二
K P
摩擦角「随围压匚3变化公式如下:
(6) 相关命令流如下:
num=arg1 in dex=arg2 单兀编号土体区域索引
!参数赋值
failO=failO%i ndex% ffail=ffail%i ndex%
kur=kur%i ndex%
kb=kb%i ndex% k=k%i ndex% m=m%i ndex% n=n%in dex%
nur=n%in dex% c=c%i ndex%*1e4 rf=rf%i ndex% den s0=de
ns2%i ndex%
*afu n, deg
pa=1e5 p1=-arrs3( num) p3=-
arrs1( num) 标准大气压最大主应力最小主应力
*if,p3,lt,0.1*pa,then
p3=0.1*pa
*en dif
对土体应力进行检查，避免出现拉力状态
fail=fail0-ffail*log10(p3/pa)
str=2*(c*cos(fail)+p3*si n(fail))/(1-si
n(fail)) s=(p1-p3)/str
*if,s,gt,0.95,then
s=0.95
*en dif
*if,str_max (nu m),gt,p1-p3,a nd,s_max (nu m),gt,s,the n et=kur*pa*(p3/pa)** nur
*elseif,str_max (nu m),gt,p1-p3,a nd,s_max (nu
m),le,s,the n ei=k*pa*(p3/pa)** n et=ei*(1-rf*s)**2
s_max( nu m)=s
*elseif,str_max (nu m),le,p1-p3,a nd,s_max (nu m),gt,s,the n ei=k*pa*(p3/pa)** n
et=ei*(1-rf*s)**2
str_max (num )=p1-p3
*elseif,str_max( nu m),le,p1-p3,a nd,s_max( nu m),le,s,the n ei=k*pa*(p3/pa)** n
et=ei*(1-rf*s)**2
str_max (num )=p1-p3
s_max( nu m)=s
*en dif bt=kb*pa*(p3/pa)**m mu=(3*bt-et)/(6*bt)
*if,mu,ge,0.49,the n mu=0.49
*elseif,mu,lt,O.O1,then mu=0.01
*en dif mp,ex, nu m,et
mp, nu xy, num,mu mp,de ns,nu m,de nsO mpchg, num,num
2.4中点增量法
性参数按照下式计算。

E= 0.25E 。

0.75E m
=0.25% 0.75l m
的弹性模量和泊松比。

相关命令流如下: *create,emcalc,mac
*get,eco un t,elem,0,co unt
*if,eco un t, ne,0,the n
eitem=e In ext(0)
*do,i,1,eco unt
!根据中点增量法计算实际弹性参数
(相关参数为初始弹性模量和本次迭代求
得的弹性模量) etii=e0(eitem,1)*0.25+em(eitem,1)*0.75
muii=e0(eitem,2)*0.25+em(eitem,2)*0.75
!保存实际弹性模量到数组中，供下次迭代时使用
e0(eitem,1)=etii
e0(eitem,2)=muii
! 定义材料，并赋予给单元
mp,ex,eitem,etii
mp, nu xy,eitem,muii
mpchg,eitem,eitem
将载荷分为若干级载荷增量，
对每级载荷增量作两次有限元计算, 第2步有限元分析弹
(7) (8) 式中：E o 、％分别为初始弹性模量和泊松比; E m 、J m 分别为根据第1步结果计算
eitem=e In ext(eitem)
*enddo
*en dif
*end
2.5初始应力场计算
通常采用下面方法确定新填筑层的初始应力状态：
■■-1 = h ( 9)
3 = K0 h( 10)
K0 =0.95-sin「( 11)式中，为新填土层的重度；h为单元形心在土体表面以下的深度；K0为土体静止侧
压力系数；:为此种材料的内摩擦角。

相关命令流如下：
yitem=ce ntry(eitem)
str1=de ns2%matid%*( ndymax-yitem)*10
str3=de ns2%matid%*( ndymax-yitem)*10*(0.95-si n(fail0%matid%))
2.6位移修正算法
新填土层位移修正算法如下：
u = 2yu/(h y) ( 12)式中，h为新填土层厚度；y为土层内部节点距离新填土层表面距离；u为有限元计算位移。

相关命令流如下:
! 位移修正法
!================================================== alls
*get,ndnum,node,0,count *dim,ndu,,ndnum,3
! 将地基以下土体位移置零
!================================================== nditem=0
*do,i,1,ndnum nditem=ndnext(nditem) ndu(nditem,1)=0 ndu(nditem,2)=0 ndu(nditem,3)=0 *enddo
! 计算累积位移，同时进行位移修正
!================================================== *do,iloop,1,13
lcdef,1,iloop+1
lcdef,2,iloop+2
lcase,2
lcoper,sub,1
! 累加下层土体位移
alls
esel,none
*do,i,iloop,13 cmsel,a,e%i%
*enddo
esel,inve nsle,s,all cm,ndcm1,node
*get,ndset1,node,0,count
nditem=0
*do,i,1,ndset1 nditem=ndnext(nditem) ndu(nditem,1)=ndu(nditem,1)+ux(nditem) ndu(nditem,2)=ndu(nditem,2)+uy(nditem) ndu(nditem,3)=ndu(nditem,3)+uz(nditem) *enddo
! 记录各个土层上下边界!================================================== cmsel,s,e%iloop%
nsel,s,ext nsel,r,loc,z,4.9,5.1 cmsel,r,ndcm1 cm,nd_line1,node
cmsel,s,e%iloop% nsel,s,ext nsel,r,loc,z,4.9,5.1 cmsel,u,ndcm1 cm,nd_line2,node
cmsel,s,nd_line1 *get,ndnum_1,node,0,count ar_line1=
*dim,ar_line1,,ndnum_1,4
nditem=0
*do,i,1,ndnum_1 nditem=ndnext(nditem) ar_line1(i,1)=nditem ar_line1(i,2)=nx(nditem) ar_line1(i,3)=ny(nditem) ar_line1(i,4)=nz(nditem)
*enddo
*do,i,1,ndnum_1-1
*do,j,i+1,ndnum_1 *if,ar_line1(i,2),gt,ar_line1(j,2),then
item=ar_line1(i,1)
itemx=ar_line1(i,2)
itemy=ar_line1(i,3)
itemz=ar_line1(i,4)
ar_line1(i,1)=ar_line1(j,1)
ar_line1(i,2)=ar_line1(j,2)
ar_line1(i,3)=ar_line1(j,3)
ar_line1(i,4)=ar_line1(j,4)
ar_line1(j,1)=item
ar_line1(j,2)=itemx
ar_line1(j,3)=itemy
ar_line1(j,4)=itemz *endif
*enddo
*enddo
tb_line1= *dim,tb_line1,table,ndnum_1,1 *vfun,tb_line1(1,0),copy,ar_line1(1,2)
*vfun,tb_line1(1,1),copy,ar_line1(1,3) !========================================== ========
cmsel,s,nd_line2 *get,ndnum_2,node,0,count ar_line2= *dim,ar_line2,,ndnum_2,4 nditem=0 *do,i,1,ndnum_2 nditem=ndnext(nditem) ar_line2(i,1)=nditem ar_line2(i,2)=nx(nditem) ar_line2(i,3)=ny(nditem) ar_line2(i,4)=nz(nditem)
*enddo
*do,i,1,ndnum_2-1
*do,j,i+1,ndnum_2 *if,ar_line2(i,2),gt,ar_line2(j,2),then
item=ar_line2(i,1) itemx=ar_line2(i,2) itemy=ar_line2(i,3) itemz=ar_line2(i,4)
ar_line2(i,1)=ar_line2(j,1)
ar_line2(i,2)=ar_line2(j,2)
ar_line2(i,3)=ar_line2(j,3)
ar_line2(i,4)=ar_line2(j,4)
ar_line2(j,1)=item ar_line2(j,2)=itemx ar_line2(j,3)=itemy ar_line2(j,4)=itemz
*endif
*enddo
*enddo
tb_line2=
*dim,tb_line2,table,ndnum_2,1
*vfun,tb_line2(1,0),copy,ar_line2(1,2)
*vfun,tb_line2(1,1),copy,ar_line2(1,3)
!==================================================
! 修正i 层填土位移
!================================================== cmsel,s,e%iloop% nsle,s,all
cmsel,u,ndcm1
*get,ndset1,node,0,count
nditem=0
*do,i,1,ndset1 nditem=ndnext(nditem) h1=tb_line2(nx(nditem))-tb_line1(nx(nditem))
z1=tb_line2(nx(nditem))-ny(nditem) *if,h1,eq,0,then
ndu(nditem,1)=0
ndu(nditem,2)=0
ndu(nditem,3)=0
*else
factor1=2*z1/(h1+z1)
ndu(nditem,1)=ux(nditem)*factor1
ndu(nditem,2)=uy(nditem)*factor1
ndu(nditem,3)=uz(nditem)*factor1
*endif
*enddo
! 显示当前步累积变形
!==================================================
alls
dof,ux,uy,uz
*do,i,1,ndnum dnsol,i,u,x,ndu(i,1),ndu(i,2),ndu(i,3)
*enddo
*if,iloop,ne,13,then
esel,none *do,i,iloop+1,13
cmsel,a,e%i%
*enddo
esel,inve
nsle,s,all
*endif
/title plns,u,y
/image,save,MidRst-%iloop%,png !========================================= =========
*enddo
alls
dof,ux,uy,uz *do,i,1,ndnum
dnsol,i,u,x,ndu(i,1),ndu(i,2),ndu(i,3)
*enddo
alls
/title
/gline,1,-1
/dscal,1,1
plns,u,y
3．仿真分析命令流文
件!================================================== ! 模型1——土石坝分析!================================================== fini
/clear /filn,case1
! 定义参数!================================================== ! (1) 坝身填土
kur1=720 kb1=400 k1=600 m1=0.3 n1=0.5 c1=2.5
fail01=30 ffail1=0 rf1=0.8 dens11=1010 dens21=1893
! (2) 排水棱体kur2=1200 kb2=640 k2=840 m2=0.2 n2=0.45 c2=0 fail02=45 ffail2=5 rf2=0.7 dens12=1275 dens22=2256
! (3) 淤泥质粘土kur3=600 kb3=250 k3=500 m3=0.2 n3=0.56 c3=1.5 fail03=20 ffail3=0 rf3=0.8 dens13=981 dens23=1727
! (4) 全风化花岗岩kur4=1000 kb4=640 k4=840 m4=0.3 n4=0.5 c4=0.4 fail04=35 ffail4=2 rf4=0.7 dens14=1039 dens24=1960
变量监视
*dim,el_moni,,30,8 emntr=1027 iiloops=0 ! s1,s3,p,q,ex,nuxy,str_max,s_max ! 监视单元
! 全局索引变量
! 邓肯模型(05/29/2016) !===================== *create,dunc2,mac
num=arg1 index=arg2 ! 单元编号
! 土体区域索引
! 参数赋值fail0=fail0%index% ffail=ffail%index% kur=kur%index% kb=kb%index% k=k%index% m=m%index% n=n%index% nur=n%index% c=c%index%*1e4 rf=rf%index% dens0=dens2%index%
*afun,deg
pa=1e5
p1=-arrs3(num) p3=-arrs1(num) ! 标准大气压! 最大主应力! 最小主应力
*if,p3,lt,0.1*pa,then
p3=0.1*pa
*endif
! 对土体应力进行检查，避免出现拉力状态
fail=fail0-ffail*log10(p3/pa) str=2*(c*cos(fail)+p3*sin(fail))/(1-sin(fail)) s=(p1-p3)/str *endif
*if,str_max(num),gt,p1-p3,and,s_max(num),gt,s,then et=kur*pa*(p3/pa)**nur
*if,s,gt,0.95,then s=0.95
*elseif,str_max(num),gt,p1-p3,and,s_max(num),le,s,then ei=k*pa*(p3/pa)**n et=ei*(1-rf*s)**2 s_max(num)=s
*elseif,str_max(num),le,p1-p3,and,s_max(num),gt,s,then ei=k*pa*(p3/pa)**n et=ei*(1-rf*s)**2 str_max(num)=p1-p3
*elseif,str_max(num),le,p1-p3,and,s_max(num),le,s,then ei=k*pa*(p3/pa)**n et=ei*(1-rf*s)**2 str_max(num)=p1-p3 s_max(num)=s
*endif
bt=kb*pa*(p3/pa)**m
mu=(3*bt-et)/(6*bt)
*if,mu,ge,0.49,then
mu=0.49
*elseif,mu,lt,0.01,then
mu=0.01
*endif
mp,ex,num,et
mp,nuxy,num,mu
mp,dens,num,dens0
mpchg,num,num
*end
! 主应力结果更新( 02/04/16
00:33:18 ) !==================================================
*create,rstnew,mac
fini
/post1
set,last
etable,s1,s,1
etable,s3,s,3
*vget,arrs1(1),elem,1,etab,s1
*vget,arrs3(1),elem,1,etab,s3
当前区域材料编号 是否计算初始应力场
*get,ecount,elem,0,count
*if,ecount,ne,0,then eitem=elnext(0) *do,i,1,ecount
*if,key1,eq,1,then yitem=centry(eitem) str1=dens2%matid%*(ndymax-yitem)*10
str3=dens2%matid%*(ndymax-yitem)*10*(0.95-sin(fail0%matid%))
arrs3(eitem)=-str1
arrs1(eitem)=-str3 dunc2,eitem,matid
*if,emntr,eq,eitem,then
iiloops=iiloops+1
el_moni(iiloops,1)=p1
el_moni(iiloops,2)=p3 el_moni(iiloops,3)=str_max(eitem)
el_moni(iiloops,4)=p1-p3 el_moni(iiloops,5)=s_max(eitem)
el_moni(iiloops,6)=s el_moni(iiloops,7)=et el_moni(iiloops,8)=mu
*endif
e0(eitem,1)=et e0(eitem,2)=mu eitem=elnext(eitem)
*else
dunc2,eitem,matid
*if,emntr,eq,eitem,then
iiloops=iiloops+1
el_moni(iiloops,1)=p1 el_moni(iiloops,2)=p3
el_moni(iiloops,3)=str_max(eitem) el_moni(iiloops,4)=p1-p3
el_moni(iiloops,5)=s_max(eitem) el_moni(iiloops,6)=s
el_moni(iiloops,7)=et el_moni(iiloops,8)=mu
*endif
em(eitem,1)=et em(eitem,2)=mu eitem=elnext(eitem)
*endif
*enddo
*endif
*end
! 中点增量法 !==================================================
*create,emcalc,mac
*end
初始应力迭代
*create,elread,mac
matid=arg1
key1=arg2
*get,ecount,elem,0,count
*if,ecount,ne,0,then
eitem=elnext(0)
*do,i,1,ecount
! 根据中点增量法计算实际弹性参数( 相关参数为初始弹性模量和本次迭代求得的
弹性模量)
etii=e0(eitem,1)*0.25+em(eitem,1)*0.75 muii=e0(eitem,2)*0.25+em(eitem,2)*0.75
! 保存实际弹性模量到数组中，供下次迭代时使用e0(eitem,1)=etii e0(eitem,2)=muii
! 定义材料，并赋予给单元mp,ex,eitem,etii mp,nuxy,eitem,muii mpchg,eitem,eitem
eitem=elnext(eitem)
*enddo
*endif *end ! 提取初始状态应
力!================================================== *create,myinis,mac cm,ecm,elem nsle,s,all
*get,ndymax,node,0,mxloc,y
cmsel,s,ecm
cmsel,s,volu1
eslv,r
elread,1,1
cmsel,s,ecm
cmsel,s,volu2
eslv,r
elread,2,1
cmsel,s,ecm
cmsel,s,volu3
eslv,r
elread,3,1
cmsel,s,ecm
cmsel,s,volu4
eslv,r
elread,4,1
*end
! 计算土体弹性参数
(05/29/2016) !==================================================
*create,matnew,mac
! 更新地基土体材料参数
*do,rziloop,1,4 cmsel,s,volu%rziloop% alls,below,volu elread,rziloop,0
*enddo
*end
前处理
!================================================== /prep7
! 定义单元类型
!================================================== et,1,200,6
et,2,185
! 定义材料参数
!================================================== mp,ex,5000,1e9
mp,nuxy,5000,0.3
mp,dens,5000,1010
! 建立几何模型
k,2,2.75*13.2,13.2
k,3,kx(2)+4,ky(2) k,4,kx(3)+2.5*9.2,ky(3)-9.2 k,10,kx(4)+2*4,0 k,5,kx(4)+2,ky(4) k,6,kx(5)+2*4,0 k,7,kx(5)+2*1,ky(5)-1 k,8,kx(7)+8,ky(7) k,9,kx(6)+8,ky(6)
l,1,2
l,2,3
l,3,4
l,4,10
l,1,10 al,all a,10,6,5,4 a,7,6,9,8 aplo rect,-50,kx(9)+50,-10,0 rect,-50,kx(9)+50,-26.4,-10 aglue,all wpro,,90 kwpa,7 asbw,all kwpa,4
asbw,all numc,all wpcsys,1,0
wpro,,,90
asbw,all
kwpa,9
asbw,all
! 划分网格
!================================================== esize,3
lsel,s,,,1,3,2
lesize,all,,,9 lsel,s,,,20,26,3 lesize,all,,,1 lsel,s,,,7,11,2 lsel,a,,,15 lesize,all,,,3 lsel,s,,,2,4,2
lsel,a,,,5,18,13 lesize,all,,,20
alls
mat,5000 msha,0,2d mshk,1 amesh,1,5,1 amesh,7
amap,14,1,20,22,9
amesh,6,8,2
amesh,9,10
amesh,13
esize,,3
vext,all,,,0,0,15
aclear,all numc,elem mpcopy,,1,5000
! 定义数组参数
!==================================================
*get,ecount,elem,0,count
*dim,arrs1,,ecount
*dim,arrs3,,ecount *dim,s_max,,ecount
*dim,str_max,,ecount *dim,e0,,ecount,2
*dim,em,,ecount,2
! 单元分组
!================================================== vsel,s,loc,y,0,13.2 alls,below,volu cm,ebar,elem cm,nbar,node
ystep=13.2/13
ytorl=0.2
*do,i,1,13
cmsel,s,nbar
cmsel,s,ebar nsel,r,loc,y,ystep*(i-1)-ytorl,ystep*i+ytorl
esln,r,1 cm,e%i%,elem
*enddo
! 施加载荷和边界条件
!================================================== nsel,s,loc,y,-26.4 d,all,all,0
nsel,s,loc,x,-50 nsel,a,loc,x,131
d,all,ux,0
nsel,s,loc,z,0 nsel,a,loc,z,15
d,all,uz,0
alls
! 针对不同区域定义初始材料属性
!================================================== mp,ex,5001,1e9 mp,nuxy,5001,0.3
mp,dens,5001,dens21
mp,ex,5002,1e9
mp,nuxy,5002,0.3
mp,dens,5002,dens22
mp,ex,5003,1e9
mp,nuxy,5003,0.3
mp,dens,5003,dens23
mp,ex,5004,1e9
mp,nuxy,5004,0.3
mp,dens,5004,dens24
! 定义体组件
!================================================== vsel,s,,,1,2 vsel,a,,,4,5
vsel,a,,,7
cm,volu1,volu
alls,below,volu
emod,all,mat,5001
vsel,s,,,3
cm,volu2,volu
alls,below,volu
emod,all,mat,5002
vsel,s,,,6,12,3
cm,volu3,volu
alls,below,volu
emod,all,mat,5003
vsel,s,,,8
vsel,a,,,10,11
cm,volu4,volu
alls,below,volu
emod,all,mat,5004
alls
! 求解器
!================================================== fini
/solu
! 第1 步：激活基岩部分!================================================== antype,0 nropt,full
rescontrol,define,all,last outres,all,last
acel,,9.806
cmsel,s,ebar ekill,all
cmsel,s,volu3 alls,below,volu ekill,all
cmsel,s,volu4 alls,below,volu myinis myinis
alls
time,1 save solve rstnew parsav,all,parms
! 第2 次计算fini
/solu
alls matnew,emntr,1,1
cmsel,s,volu4 alls,below,volu emcalc
time,1 alls save,case1_2,db,,model solve
parsav,all,parms
第2 步：激活地表土层
!================================================== fini
/solu
antype,,restart,1,,continue
parres,,parms
cmsel,s,ebar
ekill,all
cmsel,s,volu3
alls,below,volu ealive,all myinis
time,2
alls
solve
rstnew
parsav,all,parms
! 第2 次计算
fini
/solu
antype,,restart,1,,continue
parres,,parms
alls
matnew,emntr,1,1
cmsel,s,ebar
ekill,all
cmsel,s,volu3
alls,below,volu
ealive,all
esel,s,live
emcalc
time,2
alls
save,case1_2,db,,model
solve
parsav,all,parms
! 分层施工!================================================== *do,iloop,1,13 parsav,all,parms
! 第4 步：激活第i 层填筑
层!================================================== fini
/solu
antype,,restart,iloop+1,,continue parres,,parms
cmsel,s,ebar
ekill,all
cmsel,s,volu3 alls,below,volu ealive,all
*if,iloop,ne,1,then
esel,none
*do,ii,1,iloop-1
cmsel,a,e%ii%
*enddo
ealive,all
*endif
cmsel,s,e%iloop% ealive,all myinis
time,iloop+2
alls
solve
rstnew
parsav,all,parms
! 第2 次计算
fini
/solu
antype,,restart,iloop+1,,continue parres,,parms
alls
matnew,emntr,1,1
cmsel,s,ebar
ekill,all
cmsel,s,volu3
alls,below,volu
ealive,all
esel,none
*do,ii,1,iloop
cmsel,a,e%ii%
*enddo
ealive,all
esel,s,live
emcalc
time,iloop+2
alls
save,case1_%iloop+2%,db,,model
solve
parsav,all,parms
*enddo
!================================================== ! 结果后处理
!================================================== fini
/post1
显示位置设置
!================================================== /view,1,,,1
/ang,1
/auto,1
/gline,1,-1
/dscal,1,1
! 位移修正法
!================================================== alls
*get,ndnum,node,0,count *dim,ndu,,ndnum,3
! 将地基以下土体位移置零
!================================================== nditem=0
*do,i,1,ndnum nditem=ndnext(nditem) ndu(nditem,1)=0 ndu(nditem,2)=0 ndu(nditem,3)=0 *enddo
! 计算累积位移，同时进行位移修正
!================================================== *do,iloop,1,13
lcdef,1,iloop+1
lcdef,2,iloop+2
lcase,2
lcoper,sub,1
! 累加下层土体位移
alls
esel,none
*do,i,iloop,13 cmsel,a,e%i%
*enddo
esel,inve nsle,s,all cm,ndcm1,node
*get,ndset1,node,0,count
nditem=0
*do,i,1,ndset1 nditem=ndnext(nditem) ndu(nditem,1)=ndu(nditem,1)+ux(nditem) ndu(nditem,2)=ndu(nditem,2)+uy(nditem) ndu(nditem,3)=ndu(nditem,3)+uz(nditem) *enddo
! 记录各个土层上下边界
!================================================== cmsel,s,e%iloop%
nsel,s,ext nsel,r,loc,z,4.9,5.1 cmsel,r,ndcm1 cm,nd_line1,node
cmsel,s,e%iloop%
nsel,s,ext nsel,r,loc,z,4.9,5.1 cmsel,u,ndcm1 cm,nd_line2,node
cmsel,s,nd_line1 *get,ndnum_1,node,0,count
ar_line1= *dim,ar_line1,,ndnum_1,4
nditem=0
*do,i,1,ndnum_1 nditem=ndnext(nditem) ar_line1(i,1)=nditem ar_line1(i,2)=nx(nditem) ar_line1(i,3)=ny(nditem) ar_line1(i,4)=nz(nditem)
*enddo
*do,i,1,ndnum_1-1
*do,j,i+1,ndnum_1 *if,ar_line1(i,2),gt,ar_line1(j,2),then
item=ar_line1(i,1) itemx=ar_line1(i,2) itemy=ar_line1(i,3) itemz=ar_line1(i,4)
ar_line1(i,1)=ar_line1(j,1) ar_line1(i,2)=ar_line1(j,2) ar_line1(i,3)=ar_line1(j,3)
ar_line1(i,4)=ar_line1(j,4)
ar_line1(j,1)=item ar_line1(j,2)=itemx ar_line1(j,3)=itemy ar_line1(j,4)=itemz
*endif
*enddo
*enddo
tb_line1= *dim,tb_line1,table,ndnum_1,1 *vfun,tb_line1(1,0),copy,ar_line1(1,2)
*vfun,tb_line1(1,1),copy,ar_line1(1,3) !=========================================== =======
cmsel,s,nd_line2 *get,ndnum_2,node,0,count ar_line2= *dim,ar_line2,,ndnum_2,4 nditem=0
*do,i,1,ndnum_2 nditem=ndnext(nditem) ar_line2(i,1)=nditem ar_line2(i,2)=nx(nditem) ar_line2(i,3)=ny(nditem) ar_line2(i,4)=nz(nditem)
*enddo
*do,i,1,ndnum_2-1
*do,j,i+1,ndnum_2 *if,ar_line2(i,2),gt,ar_line2(j,2),then
item=ar_line2(i,1) itemx=ar_line2(i,2) itemy=ar_line2(i,3) itemz=ar_line2(i,4)
ar_line2(i,1)=ar_line2(j,1)
ar_line2(i,2)=ar_line2(j,2)
ar_line2(i,3)=ar_line2(j,3)
ar_line2(i,4)=ar_line2(j,4) ar_line2(j,1)=item ar_line2(j,2)=itemx ar_line2(j,3)=itemy
ar_line2(j,4)=itemz
*endif
*enddo
*enddo
tb_line2=
*dim,tb_line2,table,ndnum_2,1 *vfun,tb_line2(1,0),copy,ar_line2(1,2)
*vfun,tb_line2(1,1),copy,ar_line2(1,3) !=========================================== =======
! 修正i 层填土位移
!================================================== cmsel,s,e%iloop%
nsle,s,all
cmsel,u,ndcm1
*get,ndset1,node,0,count
nditem=0
*do,i,1,ndset1 nditem=ndnext(nditem) h1=tb_line2(nx(nditem))-tb_line1(nx(nditem)) z1=tb_line2(nx(nditem))-ny(nditem)
*if,h1,eq,0,then ndu(nditem,1)=0 ndu(nditem,2)=0 ndu(nditem,3)=0
*else factor1=2*z1/(h1+z1) ndu(nditem,1)=ux(nditem)*factor1
ndu(nditem,2)=uy(nditem)*factor1 ndu(nditem,3)=uz(nditem)*factor1 *endif
*enddo
显示当前步累积变形
!================================================== alls dof,ux,uy,uz
*do,i,1,ndnum dnsol,i,u,x,ndu(i,1),ndu(i,2),ndu(i,3)
*enddo
*if,iloop,ne,13,then
esel,none
*do,i,iloop+1,13
cmsel,a,e%i%
*enddo
esel,inve nsle,s,all *endif
/title
plns,u,y
/image,save,MidRst-%iloop%,png
!==================================================
*enddo
alls
dof,ux,uy,uz
*do,i,1,ndnum dnsol,i,u,x,ndu(i,1),ndu(i,2),ndu(i,3) *enddo
alls
/title
/gline,1,-1
/dscal,1,1
plns,u,y。