土石坝有限元分析

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

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

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

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

相关程序可能存在错误，笔者也未能完全认识到，仅做参考）图1 坝体结构示意图计算所采用参数详见表1所示。

仿真分析结果如下图。

坝身填土排水棱体淤泥质粘土全风化花岗岩图2竖向沉降位移云图图3水平方向位移云图图4最大主应力云图图5最小主应力云图仿真分析流程图如下图。

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

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

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

有限元网格模型如下图。

图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! 浇筑层厚度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地基最小主应力分布云图相关命令流如下：!求解器!==================================================/solu!第1步：激活基岩部分!================================================== antype,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! 保存参数信息!第2次计算fini/soluallsmatnew,emntr,1,1! 根据eb模型计算土体弹性参数cmsel,s,volu4alls,below,voluemcalc! 根据中点增量法计算土体实际参数time,1! 计算allssave,case1_1,db,,modelparsav,all,parms!第2步：激活地表土层!==================================================fini/soluantype,,restart,1,,continue! 重启动分析parres,,parmsalls! 杀死坝体填筑层cmsel,s,ebarekill,allcmsel,s,volu3! 计算地表覆盖层土体初始应力alls,below,voluealive,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,2 allssave,case1_2,db,,model solveparsav,all,parms2.3 邓肯-张模型（Duncan -Chang EB Model ）邓肯-张模型根据单元的应力状态来评估弹性参数。

当()()max 3131--σσσσ<且maxs s <时，单元处于卸荷状态，弹性模量用ur E ；否则，单元处于加荷状态，弹性模量用t E 表示。

()max 31-σσ为历史最大偏应力；m ax s 为历史最大应力水平。

卸荷或重复加载时的回弹模量采用下式计算：urn a a ur ur P P K E ⎪⎪⎭⎫⎝⎛=3σ（1）加荷时，弹性模量t E 采用下式计算：()2313131⎥⎥⎦⎤⎢⎢⎣⎡---⎪⎪⎭⎫⎝⎛=f fna a t R P KP E σσσσσ （2）()ϕϕσϕσσsin 1sin 2cos 2-331-+=c f（3）()()ultffR 3131σσσσ--=（4）式中：1σ和3σ分别为最大和最小主应力；a P 为标准大气压；f R 为破坏比；ur K 为卸荷时弹性系数；K 为加荷时弹性系数；()f 31-σσ为破坏剪应力；()ult 31-σσ为剪应力极限值。

切线体积模量采用下式计算：ma ab t P P K B ⎪⎪⎭⎫⎝⎛=3σ（5）式中，b K 为体积模量系数；m 为体积模量指数。

摩擦角ϕ随围压3σ变化公式如下：⎪⎪⎭⎫⎝⎛∆=a P 30lg -σϕϕϕ （6）相关命令流如下： 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 *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,then s=0.95 *endif*if,str_max(num),gt,p1-p3,and,s_max(num),gt,s,then et=kur*pa*(p3/pa)**nur*elseif,str_max(num),gt,p1-p3,and,s_max(num),le,s,then ei=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,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 *endifbt=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 *endifmp,ex,num,et mp,nuxy,num,mu mp,dens,num,dens0mpchg,num,num2.4 中点增量法将载荷分为若干级载荷增量，对每级载荷增量作两次有限元计算，第2步有限元分析弹性参数按照下式计算。

m E E E 75.025.00+=（7） m μμμ75.025.00+=（8）式中：0E 、0μ分别为初始弹性模量和泊松比；m E 、m μ分别为根据第1步结果计算的弹性模量和泊松比。

相关命令流如下：*create,emcalc,mac *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 *end2.5 初始应力场计算通常采用下面方法确定新填筑层的初始应力状态：h γσ=1 （9） h K γσ03=（10） ϕsin 95.00-=K（11）式中，γ为新填土层的重度；h 为单元形心在土体表面以下的深度；0K 为土体静止侧压力系数；ϕ为此种材料的内摩擦角。