桩基节点单元3-5final.txt

是用apdl语言编的

!定义相关参量
*DIM,MK,ARRAY,1,4,, , ,     !定义水平抗力系数比例常数m的矩阵
*SET,MK(1,1,1) , 6000       !第四层土的m值
*SET,MK(1,2,1) , 5000       !第三层土的m值
*SET,MK(1,3,1) , 4000       !第二层土的m值
*SET,MK(1,4,1) , 2000       !第一层土的m值
*DIM,FQ,ARRAY,1,4,, , ,     !定义桩周单位面积摩擦力
*SET,FQ(1,1,1) , 40.53      !中性点以下的单位面积正摩擦力
*SET,FQ(1,2,1) , -33.78     !中性点以上第三层的单位面积负摩擦力
*SET,FQ(1,3,1) , -16.14  
*SET,FQ(1,4,1) , -4.8
*SET,QX, 100                !水平荷载
*SET,QZ, 100
*SET,b0, 1
*SET,h0, 0.5


!进入前处理器
/prep7
et,1,solid45                !桩身单元
et,2,combin14               !弹簧单元
et,3,surf154                !表面效应单元
keyopt,1,1,0
keyopt,2,3,0                !作为纵向三维弹簧
keyopt,3,4,1                !不包含中间节点
r,1
mp,ex,1,3.25e10             !桩的材料属性
mp,nuxy,1,0.3
mp,dens,1,2700

!建立整个桩模型的节点
n,1                         
n,6,1
fill,1,6
ngen,6,6,1,6,1,,,0.2
ngen,21,36,1,36,1,,0.5
/view,1,1,1,1               !以正等侧显示
/ang,1
/REP,FAST
GPLOT

!建立桩的几何实体模型
k,1
k,2,1
k,3,1,,1
k,4,,,1
k,5,,10,
k,6,1,10,
k,7,1,10,1
k,8,,10,1
/PNUM,kp,1
kplot
v,1,2,3,4,5,6,7,8
lsel,s,loc,y,0
lsel,a,loc,y,10
lesize,all,,,5
allsel
lsel,s,loc,y,5
lesize,all,,,20
vmesh,1
nummrg,all                  !压缩重复元素
numcmp,all                  !合并重复元素

!建立第一组弹簧模型,共7*16=102个弹簧
nsel,s,loc,x,1
nsel,r,loc,y,-0.1,8.1
nplot
ngen,2,1000,all,,,2
type,2
*do,i,1,17,1
 *if,i,le,4,then            !判断属于哪一层土
   kk=1
   *elseif,i,le,8,then
   kk=2
   *elseif,i,le,12,then
   kk=3
   *elseif,i,le,17,then
   kk=4
 *endif
 *do,j,6,36,6
 *get,vy,node,j+(i-1)*36,loc,y       !获取节点的Y坐标值
    r,kk+1,mk(1,kk,1)*abs(vy-8)*b0*h0/6      !根据Y值计算水平抗力系数，即弹簧的刚度
    real,kk+1
    e,j+(i-1)*36,1000+j+(i-1)*36
 *enddo
*enddo
allsel
!最大单元号602

!第二组弹簧单元,共7*16=102个弹簧
allsel
nsel,s,loc,x,0
nsel,r,loc,y,-0.1,8.1
nsel,r,loc,z,-0.1,1.1
nplot
ngen,2,2000,all,,,-2
type,2
*do,i,1,17,1
 *if,i,le,4,then
   kk=1
   *elseif,i,le,8,then
   kk=2
   *elseif,i,le,12,then
   kk=3
   *elseif,i,le,17,then
   kk=4
 *endif
 *do,j,1,31,6
 *get,vy,node,j+(i-1)*36,loc,y 
    r,kk+10,mk(1,kk,1)*abs(vy-8)*b0*h0/6
    real,kk+10
    e,j+(i-1)*36,2000+j+(i-1)*36
 *enddo
*enddo
allsel
!最大单元号704

!第三组弹簧单元,共7*16=102个弹簧
local,11,0,,,,,90
csys,0
allsel
nsel,s,loc,z,1
nsel,r,loc,y,-0.1,8.1
nsel,r,loc,x,-0.1,1.1
nplot
ngen,2,3000,all,,,,,2
type,2
*do,i,1,17,1
 *if,i,le,4,then
   kk=1
   *elseif,i,le,8,then
   kk=2
   *elseif,i,le,12,then
   kk=3
   *elseif,i,le,17,then
   kk=4
 *endif
 *do,j,31,36,1
 *get,vy,node,j+(i-1)*36,loc,y 
    r,kk+20,mk(1,kk,1)*abs(vy-8)*b0*h0/6
    real,kk+20
    esys,0
    e,j+(i-1)*36,3000+j+(i-1)*36
 *enddo
*enddo
allsel
!最大单元号码806

!第四组弹簧单元,共7*16=102个弹簧
allsel
nsel,s,loc,z,0
nsel,r,loc,y,-0.1,8.1
nsel,r,loc,x,-0.1,1.1
nplot
ngen,2,4000,all,,,,,-2
type,2
*do,i,1,17,1
 *if,i,le,4,then
   kk=1
   *elseif,i,le,8,then
   kk=2
   *elseif,i,le,12,then
   kk=3
   *elseif,i,le,17,then
   kk=4
 *endif
 *do,j,1,6,1
 *get,vy,node,j+(i-1)*36,loc,y 
    r,kk+20,mk(1,kk,1)*abs(vy-8)*b0*h0/6
    real,kk+20
    esys,0
    e,j+(i-1)*36,4000+j+(i-1)*36
 *enddo
*enddo
allsel
!最大单元号码为908

!底层弹簧
allsel
csys,0
ngen,2,5000,1,36,1,,-1
type,2
r,1000,mk(1,1,1)*8*b0*h0/6*20
real,1000
*do,i,1,36,1
  e,i,5000+i
*enddo
allsel
gplot

!第一组表面效应单元
allsel
r,1001
real,1001
local,13,0,,,,90
type,3
*do,i,1,16,1
 *do,j,6,30,6
  nsel,s,node,,j+(i-1)*36
  nsel,a,node,,j+6+(i-1)*36
  nsel,a,node,,j+42+(i-1)*36
  nsel,a,node,,j+36+(i-1)*36
  esys,13
  esurf
 *enddo
*enddo

!第二组表面效应单元
allsel
local,13,0,,,,90
type,3
*do,i,1,16,1
 *do,j,1,29,6
  nsel,s,node,,j+(i-1)*36
  nsel,a,node,,j+6+(i-1)*36
  nsel,a,node,,j+42+(i-1)*36
  nsel,a,node,,j+36+(i-1)*36
  esys,13
  esurf
 *enddo
*enddo

!第三组表面效应单元
allsel
local,14,0,,,,90
type,3
*do,i,1,16,1
 *do,j,1,5,1
  nsel,s,node,,j+(i-1)*36
  nsel,a,node,,j+1+(i-1)*36
  nsel,a,node,,j+37+(i-1)*36
  nsel,a,node,,j+36+(i-1)*36
  esys,14
  esurf
 *enddo
*enddo
csys,0

!第四组表面效应单元
allsel
local,15,0,,,,90
type,3
*do,i,1,16,1
 *do,j,31,35,1
  nsel,s,node,,j+(i-1)*36
  nsel,a,node,,j+1+(i-1)*36
  nsel,a,node,,j+37+(i-1)*36
  nsel,a,node,,j+36+(i-1)*36
  esys,15
  esurf
 *enddo
*enddo
save
finish

!进入求解器
/solu
csys,0
allsel
nsel,s,loc,x,3                
nsel,a,loc,x,-2
nsel,a,loc,z,3
nsel,a,loc,z,-2
d,all,all                   !约束边界条件
nsel,a,loc,y,-1
d,all,all
nsel,a,loc,y,0
d,all,ux
d,all,uz
ACEL,0,9.81,0,              !施加重力加速度

!桩顶表面平均分布载荷
allsel
nsel,s,loc,y,10
nsel,r,loc,x,-0.1,1.1
nsel,r,loc,z,-0.1,1.1
nplot
esln,s
sfe,all,6,pres,,5000        !桩顶表面分布荷载
nsel,r,loc,x,1
f,all,fx,QX                  !水平载荷
nsel,s,loc,y,10
nsel,r,loc,z,1
f,all,fz,QZ                  !水平载荷
allsel
gplot

!第四层土的表面正摩擦力
csys,0
allsel
nsel,s,loc,y,0,1.5
esln,s
esel,r,type,,3
sfe,all,2,pres,,FQ(1,1,1)
/PSF,PRES,TANx,2,0,1
gplot

!第三层土的表面负摩擦力
csys,0
allsel
nsel,s,loc,y,2.5,3.5
esln,s
esel,r,type,,3
sfe,all,2,pres,,FQ(1,2,1)
/PSF,PRES,TANx,2,0,1
gplot

!第二层土的表面负摩擦力
csys,0
allsel
nsel,s,loc,y,4.5,5.5
esln,s
esel,r,type,,3
sfe,all,2,pres,,FQ(1,3,1)
/PSF,PRES,TANx,2,0,1
gplot

!第一层土的表面负摩擦力
csys,0
allsel
nsel,s,loc,y,6.5,7.5
esln,s
esel,r,type,,3
sfe,all,2,pres,,FQ(1,4,1)
/PSF,PRES,TANx,2,0,1
gplot

allsel
gplot
time,1
autot,on
nlgeom,on
nropt,full
csys,0
solve
finish

/POST1
!PLNSOL,U,Y,1,1
!PLNSOL,S,Z,1,1

PATH,path1,2,30,20, 
PPATH,P51X,1
PATH,STAT   
 
/PBC,PATH,1 
!/REPLOT 
/PBC,PATH,0 
  
AVPRIN,0,0, 
 
PDEF,SY,S,Y,AVG 
/PBC,PATH, ,0   
 
PRANGE,1,1,10,YG,   
PLPATH,SY