!匀质圆环,在原程序上加了静水压力，故本命令流适合水土分算 fini /cle

*dim,wxn,array,2000 !!定义名为wxn的数组，2000行1列，下同 *dim,wyn,array,2000 *dim,xn,array,2000 *dim,yn,array,2000 *dim,fxn,array,2000 *dim,fyn,array,2000 *dim,aa1,array,2000 *dim,bb1,array,2000 *dim,aa2,array,2000 *dim,bb2,array,2000 *dim,jx1,array,2000 *dim,jy1,array,2000 *dim,jx2,array,2000 *dim,jy2,array,2000

定义参数并赋值 nodesum=120 单元数 pi=3.1415927 height=0.35 环厚0.35米 width=1.0 环宽 area=height*width 面积 inertia=width*height*height*height/12 惯性矩 emod=34500*1000000 弹性模量，按C50混凝土计 radius=2.925 计算半径 density=25*1000/10 材料密度 wradius=radius+0.5 弹簧节点所处半径 tankang=12*1000000*width 弹簧刚度,MPa/m 赋荷载值 ptop=200*1000*width 顶板压力，单位按KN计 pltop=80*1000*width 侧向压力 plbot=50*1000*width 侧向沿竖向增加的侧压，即墙底水平压力为pltop+plbot pbot=140*1000*width 底板压力 water=8.0*1000*width 顶板处水位埋深

/prep7 et,1,beam3 定义单元类型,弹模,密度，实常数 mp,ex,1,emod mp,dens,1,density mp,prxy,1,0.2 r,1,area,inertia,height et,2,combin39 定义弹簧单元 keyopt,2,1,0 0为沿与加载相同的路径卸载 keyopt,2,2,0 按照所定义的受压曲线加压力荷载 keyopt,2,3,0 沿节点X轴方向平动 keyopt,2,4,3 二维轴向单元 keyopt,2,6,0 基本单元输出

*do,i,1,nodesum,1 xn(i)=radius*cos(pi/2+(i-1)*(2*pi/nodesum)) 定义管片节点I的x坐标 yn(i)=radius*sin(pi/2+(i-1)*(2*pi/nodesum)) 定义管片节点I的y坐标 n,i,xn(i),yn(i) 生成管片节点I *enddo

*do,i,1,nodesum,1 wxn(i)=wradius*cos(pi/2+(i-1)*(2*pi/nodesum)) 定义弹簧节点I的x坐标 wyn(i)=wradius*sin(pi/2+(i-1)*(2*pi/nodesum)) 定义弹簧节点I的y坐标 n,i+nodesum,wxn(i),wyn(i) 生成弹簧节点I *enddo

生成单元，类型为1，实常数为1 type,1 real,1 *do,i,1,nodesum-1,1 e,i,i+1 *enddo e,nodesum,1 生成弹簧单元 len=radius*(2*pi/nodesum) wytankang=len*tankang r,2,-1,-wytankang,0,0,1,wytankang/1000000 type,2 real,2

*do,i,1,nodesum,1 e,i,i+nodesum *enddo

*do,i,2,nodesum,1 aa1(i)=abs(xn(i)-xn(i-1)) 取两节点X坐标差的绝对值，并赋给数组aa1 bb1(i)=abs(yn(i)-yn(i-1)) jx1(i)=xn(i)*xn(i-1) 取两节点X坐标的乘积，并赋给数组jx1 jy1(i)=yn(i)*yn(i-1) *enddo aa1(1)=abs(xn(1)-xn(nodesum)) bb1(1)=abs(yn(1)-yn(nodesum)) jx1(1)=xn(1)*xn(nodesum) jy1(1)=yn(1)*yn(nodesum)

*do,i,1,nodesum-1,1 aa2(i)=abs(xn(i)-xn(i+1)) 相邻两节点坐标差的绝对值 bb2(i)=abs(yn(i)-yn(i+1)) jx2(i)=xn(i)*xn(i+1) 相邻两节点x坐标乘积 jy2(i)=yn(i)*yn(i+1) *enddo aa2(nodesum)=abs(xn(nodesum)-xn(1)) bb2(nodesum)=abs(yn(nodesum)-yn(1)) jx2(nodesum)=xn(nodesum)*xn(1) jy2(nodesum)=yn(nodesum)*yn(1)

*do,i,1,nodesum,1 *if,yn(i),lt,0,then 如yn(i)小于0（x轴以下的节点），则

a1=aa1(i) a2=aa2(i) fyn(i)=pbot*(a1+a2)/2

*elseif,yn(i),gt,0 如yn(i)大于0（x轴以上的节点），则 a1=aa1(i) a2=aa2(i) fyn(i)=-ptop*(a1+a2)/2 *else fyn(i)=0 *endif *enddo

*do,i,1,nodesum,1 *if,xn(i),lt,0,then 如xn(i)小于0(y轴左侧节点)，则 b1=bb1(i) b2=bb2(i) px=pltop+plbot*(radius-yn(i))/(2*radius) fxn(i)=px*(b1+b2)/2

*elseif,xn(i),gt,0 如xn(i)大于0(y轴右侧节点)，则 b1=bb1(i) b2=bb2(i) px=-(pltop+plbot*(radius-yn(i))/(2*radius)) fxn(i)=px*(b1+b2)/2 *else fxn(i)=0 *endif *enddo

进入求解层 /solu d,nodesum+1,ux,0,,nodesum+nodesum,1 弹簧节点x,y轴进行约束 d,nodesum+1,uy,0,,nodesum+nodesum,1 d,1+nodesum/2,ux,0 将圆环最下侧节点进行水平约束 ACEL,0,10,0 *do,i,1,nodesum,1 对管片节点施加节点力 f,i,fx,fxn(i) f,i,fy,fyn(i) *enddo

施加静水压力 *do,i,1,nodesum,1 ni=nelem(i,1) 得到单元I的起始节点号 nj=nelem(i,2) 得到单元I的末节点号 water1=(water-ny(ni)+radius)*10 单元I处i节点静水压力 water2=(water-ny(nj)+radius)*10 sfbeam,i,1,pres,-water1,-water2 *enddo NROPT,FULL,, !采用全牛顿－拉普森法进行求解 Allsel !选择所有内容 Outres,all,all !输出所有内容 solve finish 进入后处理 /post1 allsel,all esel,s,type,,1 /AUTO,1 相当于fit /post1 /pnum,node,0 /RGB,INDEX,100,100,100, 0 /RGB,INDEX, 80, 80, 80,13 /RGB,INDEX, 60, 60, 60,14 /RGB,INDEX, 0, 0, 0,15 /REPLOT /COLOR,PBAK,OFF ! 关闭背景颜色

! 画变形结果图 /TITLE,deform diagram pldisp,1 ! 保存为图片 /image,save,deform_1,jpg

建立单元表 etable,imoment,smisc,6 etable,jmoment,smisc,12 etable,ishear,smisc,2 etable,jshear,smisc,8 etable,iforce,smisc,1 etable,jforce,smisc,7 !更新单元表 ETABLE,REFL !画弯矩分布图 /TITLE,Bending moment diagram plls,imoment,jmoment,-1 /image,save,Bending_moment_1,jpg ! 画剪力分布图