土石坝有限元分析.docx

《土石坝有限元分析.docx》由会员分享，可在线阅读，更多相关《土石坝有限元分析.docx（52页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。

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

2、识到，仅做参考）全风化花岗岩淤泥质粘土坝身填土排水棱体图1 坝体结构示意图计算所采用参数详见表1所示。仿真分析结果如下图。图2 竖向沉降位移云图图3 水平方向位移云图图4 最大主应力云图图5 最小主应力云图仿真分析流程图如下图。建立几何模型划分网格针对每次浇筑层，建立单元组件计算地基初始弹性参数地基初始应力场计算使用EB模型计算弹性参数中点增量法计算弹性参数激活下一层土初始应力场计算使用EB模型计算弹性参数否填至坝顶是结束图6 仿真分析流程图2 关键仿真分析过程2.1 网格划分与单元组件创建当几何模型比较规则时，尽可能采用映射方式划分网格，网格分布规则，位移结果过渡光滑一些。一般情况下几何模型

3、比较复杂，此时建议将截面网格尺寸设置小一些，可以设置为每次浇筑层厚度的四分之一。采用扫略的方式划分网格，扫略方向可以设置少一些网格，控制整体网格数量。有限元网格模型如下图。图7 有限元网格模型坝体浇筑分为13步完成，每次浇筑层厚度为1m，根据竖向坐标选取浇筑层单元，创建单元组件，如图8和图9所示。图8 创建单元组件图9 单元示意图相关命令流程序如下：!单元分组!= vsel,s,loc,y,0,13.2! 选择坝体几何体alls,below,volu! 选择坝体单元和节点cm,ebar,elem! 创建单元组件ebarcm,nbar,node! 创建节点组件nbarystep=13.2/13!

4、 浇筑层厚度ytorl=0.2! 选择重叠区域范围*do,i,1,13! 循环建立每步浇筑层组件cmsel,s,nbarcmsel,s,ebarnsel,r,loc,y,ystep*(i-1)-ytorl,ystep*i+ytorlesln,r,1cm,e%i%,elem! 组件名格式为exx*enddo2.2 初始应力场计算初始应力场计算时，采用生死单元法抑制所有填筑层土体，仅保留地基土体处于激活状态。由于地基部分包含了两种土层：基岩和地表覆盖层。所以分为两步计算土体初始应力场。图10 地基最小主应力分布云图相关命令流如下：!求解器!= fini/solu!第1步：激活基岩部分!= anty

5、pe,0nropt,fullrescontrol,define,all,last! 仅通过最后子步重启动分析outres,all,last! 保存在最后子步保存所有结果acel,9.806! 施加重力加速度载荷cmsel,s,ebar! 杀死所有坝体填筑层单元ekill,allcmsel,s,volu3! 杀死地表覆盖层土体单元alls,below,voluekill,allcmsel,s,volu4! 计算基岩初始应力场alls,below,volumyinismyinisalls! 计算time,1savesolverstnew! 提取应力结果parsav,all,parms! 保存参数信

6、息!第2次计算fini/soluallsmatnew,emntr,1,1! 根据eb模型计算土体弹性参数cmsel,s,volu4alls,below,voluemcalc! 根据中点增量法计算土体实际参数time,1! 计算allssave,case1_1,db,modelsolveparsav,all,parms!第2步：激活地表土层!= fini/soluantype,restart,1,continue! 重启动分析parres,parmsalls! 杀死坝体填筑层cmsel,s,ebarekill,allcmsel,s,volu3! 计算地表覆盖层土体初始应力alls,below,v

7、oluealive,allmyinistime,2! 求解allssave,case1_2,db,modelsolverstnewparsav,all,parms!第2次计算fini/soluantype,restart,1,continue! 重启动分析parres,parmsalls! 按照EB模型计算土体弹性参数matnew,emntr,1,1cmsel,s,ebarekill,allcmsel,s,volu3alls,below,voluealive,allesel,s,live! 按照中点增量法计算土体弹性参数emcalctime,2allssave,case1_2,db,model

8、solveparsav,all,parms2.3 邓肯-张模型（Duncan-Chang EB Model）邓肯-张模型根据单元的应力状态来评估弹性参数。当且时，单元处于卸荷状态，弹性模量用；否则，单元处于加荷状态，弹性模量用表示。为历史最大偏应力；为历史最大应力水平。卸荷或重复加载时的回弹模量采用下式计算：（1）加荷时，弹性模量采用下式计算：（2）（3）（4）式中：和分别为最大和最小主应力；为标准大气压；为破坏比；为卸荷时弹性系数；为加荷时弹性系数；为破坏剪应力；为剪应力极限值。切线体积模量采用下式计算：（5）式中，为体积模量系数；为体积模量指数。摩擦角随围压变化公式如下：（6）相关命令流如

9、下：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%*1e4rf=rf%index%dens0=dens2%index%*afun,degpa=1e5! 标准大气压p1=-arrs3(num)! 最大主应力p3=-arrs1(num)! 最小主应力*if,p3,lt,0.1*pa,then! 对土体应力进行检查，避免出现拉力状态p3=0.1*

10、pa*endiffail=fail0-ffail*log10(p3/pa)str=2*(c*cos(fail)+p3*sin(fail)/(1-sin(fail)s=(p1-p3)/str*if,s,gt,0.95,thens=0.95*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,thenei=k*pa*(p3/pa)*net=ei*(1-rf*s)*2s_max(num)=s*elsei

11、f,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-et)/(6*bt)*if,mu,ge,0.49,thenmu=0.49*elseif,mu,lt,0

12、.01,thenmu=0.01*endifmp,ex,num,etmp,nuxy,num,mump,dens,num,dens0mpchg,num,num2.4 中点增量法将载荷分为若干级载荷增量，对每级载荷增量作两次有限元计算，第2步有限元分析弹性参数按照下式计算。（7）（8）式中：、分别为初始弹性模量和泊松比；、分别为根据第1步结果计算的弹性模量和泊松比。相关命令流如下：*create,emcalc,mac*get,ecount,elem,0,count*if,ecount,ne,0,theneitem=elnext(0)*do,i,1,ecount! 根据中点增量法计算实际弹性参数(相关

13、参数为初始弹性模量和本次迭代求得的弹性模量)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,eitemeitem=elnext(eitem)*enddo*endif*end2.5 初始应力场计算通常采用下面方法确定新填筑层的初始应力状态：（9）（10）（11）式中，为新

14、填土层的重度；为单元形心在土体表面以下的深度；为土体静止侧压力系数；为此种材料的内摩擦角。相关命令流如下：yitem=centry(eitem)str1=dens2%matid%*(ndymax-yitem)*10str3=dens2%matid%*(ndymax-yitem)*10*(0.95-sin(fail0%matid%)2.6 位移修正算法新填土层位移修正算法如下：（12）式中，为新填土层厚度；为土层内部节点距离新填土层表面距离；为有限元计算位移。相关命令流如下：!位移修正法!= alls*get,ndnum,node,0,count*dim,ndu,ndnum,3!将地基以下土体位

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

16、*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%iloop%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,n

17、dcm1cm,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)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

18、,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_line1(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

19、_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*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)=

20、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_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)=itemar_line2(j

21、,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)!= ! 修正i层填土位移!= cmsel,s,e%iloop%nsle,s,allcmsel,u,ndcm1*get,ndset1,node,0,countnditem=0*do,i,1,ndset1nditem=ndnext(nditem

22、)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)*factor1ndu(nditem,2)=uy(nditem)*factor1ndu(nditem,3)=uz(nditem)*factor1*endif*enddo! 显示当前步累积变形!= allsdof,ux,uy,uz

23、*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%,png!= *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

24、,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.2n2=0.45c2=0fail02=45ffail2=5rf2=0.7dens12=1275dens22=2256!(3) 淤泥质粘土kur3=600kb3=250k3=500m3=0.2n3=0.56c3=

25、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,8! s1,s3,p,q,ex,nuxy,str_max,s_maxemntr=1027! 监视单元iiloops=0! 全局索引变量!邓肯模型(05/29/2016)!= *create,dunc2,macnum=arg1! 单元编号index=arg2! 土

26、体区域索引! 参数赋值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%*1e4rf=rf%index%dens0=dens2%index%*afun,degpa=1e5! 标准大气压p1=-arrs3(num)! 最大主应力p3=-arrs1(num)! 最小主应力*if,p3,lt,0.1*pa,then! 对土体应力进行检查，避免出现拉力状态p3=0.1*pa*endiffail=fail0-ffail*log1

27、0(p3/pa)str=2*(c*cos(fail)+p3*sin(fail)/(1-sin(fail)s=(p1-p3)/str*if,s,gt,0.95,thens=0.95*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,thenei=k*pa*(p3/pa)*net=ei*(1-rf*s)*2s_max(num)=s*elseif,str_max(num),le,p1-p3,and,s

28、_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-et)/(6*bt)*if,mu,ge,0.49,thenmu=0.49*elseif,mu,lt,0.01,thenmu=0.01*endifmp,ex,nu

29、m,etmp,nuxy,num,mump,dens,num,dens0mpchg,num,num*end!主应力结果更新（02/04/16 00: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! 初始应力迭代!= *create,elread,macmatid=arg1! 当前区域材料编号key1=arg2! 是否计算初始应力场*get,ecount,elem,0,count*

30、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,matid*if,emntr,eq,eitem,theniiloops=iiloops+1el_moni(iiloops,1)=p1el_moni

31、(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,eitem,matid*if,emntr,eq,eitem,theniiloops=iiloops+1el_moni(iiloops,1)=p1

32、el_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)*endif*enddo*endif*end!中点增量法!= *create,emcalc,mac*get,ecount,elem,0,count*if

33、,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,eitemeitem=elnext(eitem

34、)*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,s,volu4eslv,relread,4,1*end!计算土体弹性参数(05/29/2016)!= *create,matnew,mac

35、! 更新地基土体材料参数*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,1k,2,2.75*13.2,13.2k,3,kx(2)+4,ky(2)k,4,kx(3)+2.5*9.2,ky(3)-9.2k,10,kx(4)+2*4,0k,5,kx(4)+2,ky(4)k,6,k

36、x(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,90asbw,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,e