土石壩有限元綜述綜述

1、土石壩有限元分析1 問題描述采用鄧肯-張模型對土石壩施工過程和蓄水狀態(tài)受力情況進行分析，選擇通用有限元分析軟件ANSYS作為研究平臺，計算土石壩竣工期豎向方向沉降、水平方向沉降、最大主應力和最小主應力情況。壩體結構示意圖如圖1所示。本文主要完成以下工作：采用ANSYS內部參數(shù)化設計語言APDL編寫鄧肯-張模型計算材料彈性參數(shù)；使用中點增量法計算每步施工單元材料彈性參數(shù)；根據(jù)位移修正算法，編寫專用程序對計算結果進行處理，獲得壩體沉降云圖。（本文針對每個步驟提供相應的APDL程序，方便后續(xù)研究人員進一步研究，也希望閱讀本文的讀者能夠將自己的研究成果與大家分享。相關程序可能存在錯誤，筆者也未能完全認

2、識到，僅做參考）全風化花崗巖淤泥質粘土壩身填土排水棱體圖1 壩體結構示意圖計算所采用參數(shù)詳見表1所示。仿真分析結果如下圖。圖2 豎向沉降位移云圖圖3 水平方向位移云圖圖4 最大主應力云圖圖5 最小主應力云圖仿真分析流程圖如下圖。建立幾何模型劃分網格針對每次澆筑層，建立單元組件計算地基初始彈性參數(shù)地基初始應力場計算使用EB模型計算彈性參數(shù)中點增量法計算彈性參數(shù)激活下一層土初始應力場計算使用EB模型計算彈性參數(shù)否填至壩頂是結束圖6 仿真分析流程圖2 關鍵仿真分析過程2.1 網格劃分與單元組件創(chuàng)建當幾何模型比較規(guī)則時，盡可能采用映射方式劃分網格，網格分布規(guī)則，位移結果過渡光滑一些。一般情況下幾何模型

3、比較復雜，此時建議將截面網格尺寸設置小一些，可以設置為每次澆筑層厚度的四分之一。采用掃略的方式劃分網格，掃略方向可以設置少一些網格，控制整體網格數(shù)量。有限元網格模型如下圖。圖7 有限元網格模型壩體澆筑分為13步完成，每次澆筑層厚度為1m，根據(jù)豎向坐標選取澆筑層單元，創(chuàng)建單元組件，如圖8和圖9所示。圖8 創(chuàng)建單元組件圖9 單元示意圖相關命令流程序如下：!單元分組!= vsel,s,loc,y,0,13.2! 選擇壩體幾何體alls,below,volu! 選擇壩體單元和節(jié)點cm,ebar,elem! 創(chuàng)建單元組件ebarcm,nbar,node! 創(chuàng)建節(jié)點組件nbarystep=13.2/13!

4、 澆筑層厚度ytorl=0.2! 選擇重疊區(qū)域范圍*do,i,1,13! 循環(huán)建立每步澆筑層組件cmsel,s,nbarcmsel,s,ebarnsel,r,loc,y,ystep*(i-1)-ytorl,ystep*i+ytorlesln,r,1cm,e%i%,elem! 組件名格式為exx*enddo2.2 初始應力場計算初始應力場計算時，采用生死單元法抑制所有填筑層土體，僅保留地基土體處于激活狀態(tài)。由于地基部分包含了兩種土層：基巖和地表覆蓋層。所以分為兩步計算土體初始應力場。圖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! 保存參數(shù)信

6、息!第2次計算fini/soluallsmatnew,emntr,1,1! 根據(jù)eb模型計算土體彈性參數(shù)cmsel,s,volu4alls,below,voluemcalc! 根據(jù)中點增量法計算土體實際參數(shù)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模型計算土體彈性參數(shù)matnew,emntr,1,1cmsel,s,ebarekill,allcmsel,s,volu3alls,below,voluealive,allesel,s,live! 按照中點增量法計算土體彈性參數(shù)emcalctime,2allssave,case1_2,db,model

8、solveparsav,all,parms2.3 鄧肯-張模型（Duncan-Chang EB Model）鄧肯-張模型根據(jù)單元的應力狀態(tài)來評估彈性參數(shù)。當且時，單元處于卸荷狀態(tài)，彈性模量用；否則，單元處于加荷狀態(tài)，彈性模量用表示。為歷史最大偏應力；為歷史最大應力水平。卸荷或重復加載時的回彈模量采用下式計算：（1）加荷時，彈性模量采用下式計算：（2）（3）（4）式中：和分別為最大和最小主應力；為標準大氣壓；為破壞比；為卸荷時彈性系數(shù)；為加荷時彈性系數(shù)；為破壞剪應力；為剪應力極限值。切線體積模量采用下式計算：（5）式中，為體積模量系數(shù)；為體積模量指數(shù)。摩擦角隨圍壓變化公式如下：（6）相關命令流如

9、下：num=arg1! 單元編號index=arg2! 土體區(qū)域索引! 參數(shù)賦值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! 對土體應力進行檢查，避免出現(xiàn)拉力狀態(tài)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步有限元分析彈性參數(shù)按照下式計算。（7）（8）式中：、分別為初始彈性模量和泊松比；、分別為根據(jù)第1步結果計算的彈性模量和泊松比。相關命令流如下：*create,emcalc,mac*get,ecount,elem,0,count*if,ecount,ne,0,theneitem=elnext(0)*do,i,1,ecount! 根據(jù)中點增量法計算實際彈性參數(shù)(相關

13、參數(shù)為初始彈性模量和本次迭代求得的彈性模量)etii=e0(eitem,1)*0.25+em(eitem,1)*0.75muii=e0(eitem,2)*0.25+em(eitem,2)*0.75! 保存實際彈性模量到數(shù)組中，供下次迭代時使用e0(eitem,1)=etiie0(eitem,2)=muii!定義材料，并賦予給單元mp,ex,eitem,etiimp,nuxy,eitem,muiimpchg,eitem,eitemeitem=elnext(eitem)*enddo*endif*end2.5 初始應力場計算通常采用下面方法確定新填筑層的初始應力狀態(tài)：（9）（10）（11）式中，為新

14、填土層的重度；為單元形心在土體表面以下的深度；為土體靜止側壓力系數(shù)；為此種材料的內摩擦角。相關命令流如下：yitem=centry(eitem)str1=dens2%matid%*(ndymax-yitem)*10str3=dens2%matid%*(ndymax-yitem)*10*(0.95-sin(fail0%matid%)2.6 位移修正算法新填土層位移修正算法如下：（12）式中，為新填土層厚度；為土層內部節(jié)點距離新填土層表面距離；為有限元計算位移。相關命令流如下：!位移修正法!= 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!定義參數(shù)!= !(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!變量監(jiān)視!= *dim,el_moni,30,8! s1,s3,p,q,ex,nuxy,str_max,s_maxemntr=1027! 監(jiān)視單元iiloops=0! 全局索引變量!鄧肯模型(05/29/2016)!= *create,dunc2,macnum=arg1! 單元編號index=arg2! 土

26、體區(qū)域索引! 參數(shù)賦值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! 對土體應力進行檢查，避免出現(xiàn)拉力狀態(tài)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! 當前區(qū)域材料編號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! 根據(jù)中點增量法計算實際彈性參數(shù)(相關參數(shù)為初始彈性模量和本次迭代求得的彈性模量)etii=e0(eitem,1)*0.25+em(eitem,1)*0.75muii=e0(eitem,2)*0.25+em(eitem,2)*0.75! 保存實際彈性模量到數(shù)組中，供下次迭代時使用e0(eitem,1)=etiie0(eitem,2)=muii!定義材料，并賦予給單元mp,ex,eitem,etiimp,nuxy,eitem,muiimpchg,eitem,eitemeitem=elnext(eitem

34、)*enddo*endif*end!提取初始狀態(tài)應力!= *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!計算土體彈性參數(shù)(05/29/2016)!= *create,matnew,mac

35、! 更新地基土體材料參數(shù)*do,rziloop,1,4cmsel,s,volu%rziloop%alls,below,voluelread,rziloop,0*enddo*end!前處理!= /prep7!定義單元類型!= et,1,200,6et,2,185!定義材料參數(shù)!= 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,

37、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,5000!定義數(shù)組參數(shù)!= *get,ecount,elem,0,count*dim,arrs1,ecount*dim

38、,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,-

39、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!針對不同區(qū)域定義初始材料屬性!= mp,ex,5001,1e9mp,nuxy,5001,0.3mp,dens,5001,dens21mp,ex,5002,1e9mp,nuxy,5002,0.3mp,dens,5002,dens22mp,ex,5003,1e9mp,nuxy,5003,0.3mp,dens,5003,dens23mp,ex,5004,1e9mp,nuxy,5004,0.3mp,de

40、ns,5004,dens24!定義體組件!= vsel,s,1,2vsel,a,4,5vsel,a,7cm,volu1,volualls,below,voluemod,all,mat,5001vsel,s,3cm,volu2,volualls,below,voluemod,all,mat,5002vsel,s,6,12,3cm,volu3,volualls,below,voluemod,all,mat,5003vsel,s,8vsel,a,10,11cm,volu4,volualls,below,voluemod,all,mat,5004alls!求解器!= fini/solu!第1步：激活基

41、巖部分!= antype,0nropt,fullrescontrol,define,all,lastoutres,all,lastacel,9.806cmsel,s,ebarekill,allcmsel,s,volu3alls,below,voluekill,allcmsel,s,volu4alls,below,volumyinismyinisallstime,1savesolverstnewparsav,all,parms!第2次計算fini/soluallsmatnew,emntr,1,1cmsel,s,volu4alls,below,voluemcalctime,1allssave,ca

42、se1_2,db,modelsolveparsav,all,parms!第2步：激活地表土層!= fini/soluantype,restart,1,continueparres,parmscmsel,s,ebarekill,allcmsel,s,volu3alls,below,voluealive,allmyinistime,2allssolverstnewparsav,all,parms!第2次計算fini/soluantype,restart,1,continueparres,parmsallsmatnew,emntr,1,1cmsel,s,ebarekill,allcmsel,s,vo

43、lu3alls,below,voluealive,allesel,s,liveemcalctime,2allssave,case1_2,db,modelsolveparsav,all,parms!分層施工!= *do,iloop,1,13parsav,all,parms!第4步：激活第i層填筑層!= fini/soluantype,restart,iloop+1,continueparres,parmscmsel,s,ebarekill,allcmsel,s,volu3alls,below,voluealive,all*if,iloop,ne,1,thenesel,none*do,ii,1,iloop-1cmsel,a,e%ii%*end