土石坝有限元分析

1、土石坝有限元分析1问题描述采用邓肯-张模型对土石坝施工过程和蓄水状态受力情况进行分析，选择通用有限元分析软件ANSYS作为研究平台，计算土石坝竣工期竖向方向沉降、水平方向沉降、最大主应力和最小主应力情况。坝体结构示意图如图1所示。本文主要完成以下工作：采用ANSYS内部参数化设计语言APDL编写邓肯-张模型计算材料弹性参数;使用中点增量法计算每步施工单元材料弹性参数；根据位移修正算法，编写专用程序对计算结果进行处理，获得坝体沉降云图。（本文针对每个步骤提供相应的APDL程序，方便后续研究人员进一步研究，也希望阅读本文的读者能够将自己的研究成果与大家分享。相关程序可能存在错误，笔者也未能完全认识

2、到，仅做参考）坝身填土/排水棱体淤泥质粘土uu岩全风化花岗图1坝体结构示意图计算所采用参数详见表1所示。表1各分区材料的物理力学参数材科湛容童了浮歩童fC/砒、<Pf址号(KN/rti1IIKN?m)KKKsn»&n19.3lOJtt72t>0.3i)5253<if).8223612茁t2W64nMOn.2n.4S4S弓n.717.279.H1WK)2Si*5柑<).2G.S6152finn.8均6010JR1H4T>n.3r>50435211.7注：马屮/7仿真分析结果如下图NOHALaOLUTIOWA33YS17.0Build17.0S

3、7E5-9999UYRSYS=0cmENWSKK-2L471B,2H42S=G263OT-.Qsasa.丄HT6"?-.1341S3-_CiH0fc3-.02J12a263&7图2竖向沉降位移云图图3水平方向位移云图AVGELEMENTSOLUTION3TEP=1SSUB=1TIME二1551DMX-,362925MU54015ATBtS艮合1它乱，右17.0Build17.0MX-5S4DL5-462012-3300DS-L930D5-£60口L,T-526013-396010-26007-1320030图4最大主应力云图和苗ELEKDHTSOLUTIONSTEP

4、-15SOB=1TIhE-15S3(fiVG)DMX=-3352SMN=-4巧吗2甘了SMX-32121呂-23&3d_7-7819丄起32121.6-46Q297-353&14-243SCI-13334a-40&L30-2&BS11-laesoa图5最小主应力云图仿真分析流程图如下图。图6仿真分析流程图2关键仿真分析过程2.1网格划分与单元组件创建当几何模型比较规则时，尽可能采用映射方式划分网格，网格分布规则，位移结果过渡光滑一些。一般情况下几何模型比较复杂，此时建议将截面网格尺寸设置小一些，可以设置为每次浇筑层厚度的四分之一。采用扫略的方式划分网格，扫略方向

5、可以设置少一些网格控制整体网格数量。有限元网格模型如下图。图7有限元网格模型坝体浇筑分为13步完成，每次浇筑层厚度为1m,根据竖向坐标选取浇筑层单元，创建单元组件，如图8和图9所示。AComporptrtMb冊的XCnnroncnts®NameCouncELIlament72E2Elernait72E3IlerTisit72E4El匕merit関ESElament6i3EbIlerTiGfitE7El已merit帥EEIlonsni60ESElernsitCOEWIlemmi6<lELIIlarwnt&E12Fieneit”1E12Ilansnt60|i|图8创建单元组

6、件图9单元示意图相关命令流程序如下ANSYSA17.0MA'f2&2Gl<i21=52;4«!单元分组vsel,s,loc,y,0,13.2alls,below,volucm,ebar,elemcm,nbar,node!选择坝体几何体!选择坝体单元和节点!创建单元组件ebar!创建节点组件nbarystep=13.2/13ytorl=0.2*do,i,1,13cmsel,s,nbarcmsel,s,ebar!浇筑层厚度!选择重叠区域范围!循环建立每步浇筑层组件nsel,r,loc,y,ystep*(i-1)-ytorl,ystep*i+ytorlesln,r,1

7、cm,e%i%,elem*enddo!组件名格式为exx2.2初始应力场计算初始应力场计算时，采用生死单元法抑制所有填筑层土体，仅保留地基土体处于激活状态。由于地基部分包含了两种土层：基岩和地表覆盖层。所以分为两步计算土体初始应力场。AVGELEMEWT30IAH打N3TEP-2SUB=1TIKE-28L冉呵DMX-.056951Z971DEMX=-lCl255.aM?L-251106-35341'-LS617'-iD5a72'-42126-2S5233-2011199-13-7745-74000-1.0255,9图10地基最小主应力分布云图相关命令流如下：!求解器!=

8、fini/solu!第1步：激活基岩部分!=antype,0nropt,fullrescontrol,define,all,lastoutres,all,lastacel,9.806!仅通过最后子步重启动分析!保存在最后子步保存所有结果!施加重力加速度载荷cmsel,s,ebar杀死所有坝体填筑层单元ekill,allcmsel,s,volu3alls,below,voluekill,all杀死地表覆盖层土体单元cmsel,s,volu4alls,below,volumyinismyinis计算基岩初始应力场allstime,1savesolverstnewparsav,all,parms!第

9、2次计算fini/solu计算!提取应力结果!保存参数信息allsmatnew,emntr,1,1cmsel,s,volu4alls,below,voluemcalctime,1allssave,case1_1,db,model!根据eb模型计算土体弹性参数!根据中点增量法计算土体实际参数!计算solveparsav,all,parms!第2步：激活地表土层!=fini/soluantype,restart,1,continueparres,parmsallscmsel,s,ebarekill,allcmsel,s,volu3alls,below,voluealive,allmyinistim

10、e,2allssave,case1_2,db,modelsolverstnewparsav,all,parms!第2次计算fini/soluantype,restart,1,continueparres,parmsallsmatnew,emntr,1,1cmsel,s,ebarekill,all!重启动分析!杀死坝体填筑层!计算地表覆盖层土体初始应力!求解!重启动分析!按照EB模型计算土体弹性参数cmsel,s,volu3alls,below,voluealive,allesel,s,liveemcalc!按照中点增量法计算土体弹性参数time,2allssave,case1_2,db,mod

11、elsolveparsav,all,parms2.3 邓肯-张模型(Duncan-ChangEBModel)邓肯-张模型根据单元的应力状态来评估弹性参数。当G-)<G-)13且s<s13maxmax5)弹性模量用E表示。t时，单元处于卸荷状态，弹性模量用e；否则，单元处于加荷状态,urG-b)为历史最大偏应力；s为历史最大应力水平。13maxmax卸荷或重复加载时的回弹模量采用下式计算：E=KPururannur1)加荷时，弹性模量E采用下式计算:tE二KPta1_R£fb1_G3f2)Gb)2ccos申+2bsin申13f1一sin申3)G_b)Rf13ult4)式中：

12、b和b分别为最大和最小主应力；p为标准大气压;13a荷时弹性系数；K为加荷时弹性系数；G-b)为破坏剪应力；13f为破坏比；K为卸urG-b)为剪应力极限13ult值。切线体积模量采用下式计算：式中，k为体积模量系数；m为体积模量指数。摩擦角9随围压b3变化公式如下:A9lg06)相关命令流如下num=arg1index=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%if,

13、str_max(num),gt,p1-p3,and,s_max(num),gt,s,thenet=kur*pa*(p3/pa)*nur1e4rf=rf%index%dens0=dens2%index%*afun,degpa=1e5!标准大气压p1=-arrs3(num)!最大主应力!对土体应力进行检查，避免出现拉力状态p3=-arrs1(num)!最小主应力*if,p3,lt,0.1*pa,thenp3=0.1*pa*endiffail=fail0-ffail*log10(p3/pa)str=2*(c*cos(fail)+p3*sin(fail)/(1-sin(fail)s=(p1-p3)/s

14、tr*if,s,gt,0.95,thens=0.95*endifs_max(num)=s*elseif,str_max(num),le,p1-p3,and,s_max(num),gt,s,thenei=k*pa*(p3/pa)*net=ei*(1-rf*s)*2str_max(num)=p1-p3*elseif,str_max(num),le,p1-p3,and,s_max(num),le,s,thenei=k*pa*(p3/pa)*net=ei*(1-rf*s)*2str_max(num)=p1-p3s_max(num)=s*endifbt=kb*pa*(p3/pa)*mmu=(3*bt-e

15、t)/(6*bt)*if,mu,ge,0.49,thenmu=0.49*elseif,mu,lt,0.01,thenmu=0.01*endifmp,ex,num,etmp,nuxy,num,mump,dens,num,dens0mpchg,num,num2.4中点增量法将载荷分为若干级载荷增量，对每级载荷增量作两次有限元计算，第2步有限元分析弹性参数按照下式计算。E二0.25E+0.75E(7)0mp二0.25p+0.75p(8)0m式中：E、p分别为初始弹性模量和泊松比；E、p分别为根据第1步结果计算00mm的弹性模量和泊松比。相关命令流如下：*create,emcalc,mac*get,e

16、count,elem,0,count*if,ecount,ne,0,theneitem=elnext(0)*do,i,1,ecount!根据中点增量法计算实际弹性参数(相关参数为初始弹性模量和本次迭代求得的弹性模量)etii=e0(eitem,1)*0.25+em(eitem,1)*0.75muii=e0(eitem,2)*0.25+em(eitem,2)*0.75!保存实际弹性模量到数组中，供下次迭代时使用e0(eitem,1)=etiie0(eitem,2)=muii!定义材料，并赋予给单元mp,ex,eitem,etiimp,nuxy,eitem,muiimpchg,eitem,eite

17、meitem=elnext(eitem)*enddo*endif*end2.5初始应力场计算通常采用下面方法确定新填筑层的初始应力状态：c1=yh(9)二KYh(10)30K=0.95一sin9(11)0式中，丫为新填土层的重度；h为单元形心在土体表面以下的深度；K为土体静止侧0压力系数；9为此种材料的内摩擦角。相关命令流如下：yitem=centry(eitem)str1=dens2%matid%*(ndymax-yitem)*10str3=dens2%matid%*(ndymax-yitem)*10*(0.95-sin(fail0%matid%)2.6位移修正算法新填土层位移修正算法如下：

18、U=2yu/(h+y)(12)式中，h为新填土层厚度；y为土层内部节点距离新填土层表面距离；u为有限元计算位移。相关命令流如下!位移修正法!=alls*get,ndnum,node,0,count*dim,ndu,ndnum,3!将地基以下土体位移置零!=nditem=0*do,i,1,ndnumnditem=ndnext(nditem)ndu(nditem,1)=0ndu(nditem,2)=0ndu(nditem,3)=0*enddo!计算累积位移，同时进行位移修正!=*do,iloop,1,13lcdef,1,iloop+1lcdef,2,iloop+2lcase,2lcoper,sub

19、,1!累加下层土体位移allsesel,none*do,i,iloop,13cmsel,a,e%i%*enddoesel,invensle,s,allcm,ndcm1,node*get,ndset1,node,0,countnditem=0*do,i,1,ndset1nditem=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%il

20、oop%nsel,s,extnsel,r,loc,z,4.9,5.1cmsel,r,ndcm1cm,nd_line1,nodecmsel,s,e%iloop%nsel,s,extnsel,r,loc,z,4.9,5.1cmsel,u,ndcm1cm,nd_line2,nodecmsel,s,nd_line1*get,ndnum_1,node,0,countar_line1=*dim,ar_line1,ndnum_1,4nditem=0*do,i,1,ndnum_1nditem=ndnext(nditem)ar_line1(i,1)=nditemar_line1(i,2)=nx(nditem)a

21、r_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),thenitem=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_l

22、ine1(j,4)ar_line1(j,1)=itemar_line1(j,2)=itemxar_line1(j,3)=itemyar_line1(j,4)=itemz*endif*enddo*enddotb_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,countar_line2=*dim,ar_line2,ndnum_2,4nditem=0*

23、do,i,1,ndnum_2nditem=ndnext(nditem)ar_line2(i,1)=nditemar_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),thenitem=ar_line2(i,1)itemx=ar_line2(i,2)itemy=ar_line2(i,3)itemz=ar_line2(i,4)ar_line2(i,1)=ar_li

24、ne2(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)=itemar_line2(j,2)=itemxar_line2(j,3)=itemyar_line2(j,4)=itemz*endif*enddo*enddotb_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)!=!修

25、正i层填土位移!=cmsel,s,e%iloop%nsle,s,allcmsel,u,ndcm1*get,ndset1,node,0,countnditem=0*do,i,1,ndset1nditem=ndnext(nditem)h1=tb_line2(nx(nditem)-tb_line1(nx(nditem)z1=tb_line2(nx(nditem)-ny(nditem)*if,h1,eq,0,thenndu(nditem,1)=0ndu(nditem,2)=0ndu(nditem,3)=0*elsefactor1=2*z1/(h1+z1)ndu(nditem,1)=ux(nditem)

26、*factor1ndu(nditem,2)=uy(nditem)*factor1ndu(nditem,3)=uz(nditem)*factor1*endif*enddo!显示当前步累积变形!=allsdof,ux,uy,uz*do,i,1,ndnumdnsol,i,u,x,ndu(i,1),ndu(i,2),ndu(i,3)*enddo*if,iloop,ne,13,thenesel,none*do,i,iloop+1,13cmsel,a,e%i%*enddoesel,invensle,s,all*endif/titleplns,u,y/image,save,MidRst-%iloop%,pn

27、g!=*enddoallsdof,ux,uy,uz*do,i,1,ndnumdnsol,i,u,x,ndu(i,1),ndu(i,2),ndu(i,3)*enddoalls/title/gline,1,-1/dscal,1,1plns,u,y3仿真分析命令流文件!=!模型1土石坝分析!=fini/clear/filn,case1!定义参数!=!(1)坝身填土kur1=720kb1=400k1=600m1=0.3n1=0.5c1=2.5fail01=30ffail1=0rf1=0.8dens11=1010dens21=1893!(2)排水棱体kur2=1200kb2=640k2=840m2=0.

28、2n2=0.45c2=0fail02=45ffail2=5rf2=0.7dens12=1275dens22=2256!(3)淤泥质粘土kur3=600kb3=250k3=500m3=0.2n3=0.56c3=1.5fail03=20ffail3=0rf3=0.8dens13=981dens23=1727!(4)全风化花岗岩kur4=1000kb4=640k4=840m4=0.3n4=0.5c4=0.4fail04=35ffail4=2rf4=0.7dens14=1039dens24=1960变量监视*dim,el_moni,30,8emntr=1027iiloops=0!s1,s3,p,q,e

29、x,nuxy,str_max,s_max!监视单元!全局索引变量邓肯模型(05/29/2016)*create,dunc2,macnum=arg1index=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%if,s,gt,0.95,thens=0.951e4rf=rf%index%dens0=dens2%index%*afun,degpa=1e5p1=-arrs3(num)p

30、3=-arrs1(num)!标准大气压!最大主应力!最小主应力*if,p3,lt,0.1*pa,thenp3=0.1*pa!对土体应力进行检查，避免出现拉力状态*endiffail=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,thenet=kur*pa*(p3/pa)*nur*elseif,str_max(num),gt,p1-p3,and,s_max(num),le,s,t

31、henei=k*pa*(p3/pa)*net=ei*(1-rf*s)*2s_max(num)=s*elseif,str_max(num),le,p1-p3,and,s_max(num),gt,s,thenei=k*pa*(p3/pa)*net=ei*(1-rf*s)*2str_max(num)=p1-p3*elseif,str_max(num),le,p1-p3,and,s_max(num),le,s,thenei=k*pa*(p3/pa)*net=ei*(1-rf*s)*2str_max(num)=p1-p3s_max(num)=s*endifbt=kb*pa*(p3/pa)*mmu=(3*

32、bt-et)/(6*bt)*if,mu,ge,0.49,thenmu=0.49*elseif,mu,lt,0.01,thenmu=0.01*endifmp,ex,num,etmp,nuxy,num,mump,dens,num,dens0mpchg,num,num*end!主应力结果更新(02/04/1600:33:18)!=*create,rstnew,macfini/post1set,lastetable,s1,s,1etable,s3,s,3*vget,arrs1(1),elem,1,etab,s1*vget,arrs3(1),elem,1,etab,s3*end!初始应力迭代*creat

33、e,elread,macmatid=arg1key1=arg2当前区域材料编号是否计算初始应力场*get,ecount,elem,0,count*if,ecount,ne,0,theneitem=elnext(0)*do,i,1,ecount*if,key1,eq,1,thenyitem=centry(eitem)str1=dens2%matid%*(ndymax-yitem)*10str3=dens2%matid%*(ndymax-yitem)*10*(0.95-sin(fail0%matid%)arrs3(eitem)=-str1arrs1(eitem)=-str3dunc2,eitem,

34、matid*if,emntr,eq,eitem,theniiloops=iiloops+1el_moni(iiloops,1)=p1el_moni(iiloops,2)=p3el_moni(iiloops,3)=str_max(eitem)el_moni(iiloops,4)=p1-p3el_moni(iiloops,5)=s_max(eitem)el_moni(iiloops,6)=sel_moni(iiloops,7)=etel_moni(iiloops,8)=mu*endife0(eitem,1)=ete0(eitem,2)=mueitem=elnext(eitem)*elsedunc2

35、,eitem,matid*if,emntr,eq,eitem,theniiloops=iiloops+1el_moni(iiloops,1)=p1el_moni(iiloops,2)=p3el_moni(iiloops,3)=str_max(eitem)el_moni(iiloops,4)=p1-p3el_moni(iiloops,5)=s_max(eitem)el_moni(iiloops,6)=sel_moni(iiloops,7)=etel_moni(iiloops,8)=mu*endifem(eitem,1)=etem(eitem,2)=mueitem=elnext(eitem)*en

36、dif*enddo*endif*end!中点增量法!=*create,emcalc,mac*get,ecount,elem,0,count*if,ecount,ne,0,theneitem=elnext(0)*do,i,1,ecount!根据中点增量法计算实际弹性参数(相关参数为初始弹性模量和本次迭代求得的弹性模量)etii=e0(eitem,1)*0.25+em(eitem,1)*0.75muii=e0(eitem,2)*0.25+em(eitem,2)*0.75!保存实际弹性模量到数组中，供下次迭代时使用e0(eitem,1)=etiie0(eitem,2)=muii!定义材料，并赋予给单

37、元mp,ex,eitem,etiimp,nuxy,eitem,muiimpchg,eitem,eitemeitem=elnext(eitem)*enddo*endif*end!提取初始状态应力!=*create,myinis,maccm,ecm,elemnsle,s,all*get,ndymax,node,0,mxloc,ycmsel,s,ecmcmsel,s,volu1eslv,relread,1,1cmsel,s,ecmcmsel,s,volu2eslv,relread,2,1cmsel,s,ecmcmsel,s,volu3eslv,relread,3,1cmsel,s,ecmcmsel,

38、s,volu4eslv,relread,4,1*end!计算土体弹性参数(05/29/2016)!=*create,matnew,mac!更新地基土体材料参数*do,rziloop,1,4cmsel,s,volu%rziloop%alls,below,voluelread,rziloop,0*enddo*end前处理!=/prep7!定义单元类型!=et,1,200,6et,2,185!定义材料参数!=mp,ex,5000,1e9mp,nuxy,5000,0.3mp,dens,5000,1010!建立几何模型k,2,2.75*13.2,13.2k,3,kx(2)+4,ky(2)k,4,kx(3

39、)+2.5*9.2,ky(3)-9.2k,10,kx(4)+2*4,0k,5,kx(4)+2,ky(4)k,6,kx(5)+2*4,0k,7,kx(5)+2*1,ky(5)-1k,8,kx(7)+8,ky(7)k,9,kx(6)+8,ky(6)l,1,2l,2,3l,3,4l,4,10l,1,10al,alla,10,6,5,4a,7,6,9,8aplorect,-50,kx(9)+50,-10,0rect,-50,kx(9)+50,-26.4,-10aglue,allwpro,90kwpa,7asbw,allkwpa,4asbw,allnumc,allwpcsys,1,0wpro,90asb

40、w,allkwpa,9asbw,all!划分网格!=esize,3lsel,s,1,3,2lesize,all,9lsel,s,20,26,3lesize,all,1lsel,s,7,11,2lsel,a,15lesize,all,3lsel,s,2,4,2lsel,a,5,18,13lesize,all,20allsmat,5000msha,0,2dmshk,1amesh,1,5,1amesh,7amap,14,1,20,22,9amesh,6,8,2amesh,9,10amesh,13esize,3vext,all,0,0,15aclear,allnumc,elemmpcopy,1,500

41、0!定义数组参数!=*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.2alls,below,volucm,ebar,elemcm,nbar,nodeystep=13.2/13ytorl=0.2*do,i,1,13cmsel,s,nbarcmsel,s,ebarnsel,r,loc,y,ystep*(i-1)-ytorl,ystep*i+ytorlesln,r,1cm,e%i%,elem*enddo!施加载荷和边界条件!=nsel,s,loc,y,-26.4d,all,all,0nsel,s,loc,x,-50nsel,a,loc,x,131d,all,ux,0nsel,s,loc,z,0nsel,a,loc,z,15d,all,uz,0alls!针对不同区域定义初始材料属性!=mp,ex,5001,1e9mp,nuxy,5001,0.3mp,dens,5001,dens21mp,ex,5002,1e9mp,nuxy,5002,0.3mp,dens,5002,dens22mp,ex,5003,