（有问题的）弹簧单元做的桩土接触分析

　　下面是我自己用弹簧单元做的桩和土的接触分析，但是不成功，不知道是什么原因，希望能用高人给指点一二！！！不胜感激～！！！！！
　　　程序有点冗长了！！！！不好意思！！！因为学的时间不是很长　！！！！
　　　工程背景是这样的：桩是扩顶的，上面４０cm扩大头，直径是０.７m，下面桩体的直径是０.４m，桩长１５m　，
　　　土体一共分了４层！！！具体参数见程序中！！！！

*DIM,MK,ARRAY,1,4,,,,
*SET,MK(1,1,1),3000
*SET,MK(1,2,1),2000
*SET,MK(1,3,1),3000
*SET,MK(1,4,1),5000
*DIM,FQ,ARRAY,1,5,,,,
*SET,FQ(1,1,1),-4.46
*SET,FQ(1,2,1),-14.79
*SET,FQ(1,3,1),36.95
*SET,FQ(1,4,1),78.54
*SET,FQ(1,5,1),82.94
*SET,b1,0.55
*SET,h1,0.05
*SET,b0,0.314
*SET,h0,0.1
/PREP7
et,1,plane42
et,2,SOLID45
et,3,combin14
et,4,surf154
keyopt,2,1,0
keyopt,3,3,0
keyopt,4,4,0
r,1
mp,ex,2,2.55e10
mp,nuxy,2,0.2
mp,dens,2,2500
n,1,
n,2,0.05,
n,20,-0.05,0,0
csys,5
fill,2,20,17
ngen,2,18,3,19,1,,180
ngen,4,36,2,37,1,0.05
ngen,147,145,1,145,1,,,-0.1
nsel,s,loc,z,0
nplot
ngen,4,100036,110,145,1,0.05
nsel,s,loc,z,0
ngen,9,100253,all,,,,,0.05
allsel
gplot
csys,0
k,1,
k,2,0.35,0,0
k,3,0.35,0.4,0
k,4,0,0.4,0
k,5,0.2,0,0
k,6,0,-14.6,0
k,7,0.2,-14.6,0
a,1,5,2,3,4
a,1,5,7,6
aglue,1,2
numcmp,all
lsel,s,loc,y,0.4
lesize,all,,,7
lsel,s,loc,y,0.2
lesize,all,,,8
lsel,s,,,1,7,6
lesize,all,,,4
lsel,s,,,2
lesize,all,,,3
lsel,s,,,6,8,2
lesize,all,,,146
allsel
lplot
amesh,1,2,1
TYPE, 2
EXTOPT,ESIZE,9,0,
EXTOPT,ACLEAR,1
EXTOPT,ATTR,0,0,0
MAT,_Z2
REAL,_Z4
ESYS,0
vrotat,all,,,,,,1,6,360
nummrg,all
numcmp,all
csys,5
!下层侧面弹簧
nsel,s,loc,x,0.2
nsel,r,loc,z,0,-14.6
nplot
ngen,2,100000,all,,,0.4
type,3
wpro,,90.000000,
wpro,,90.000000,
*do,i,1,147,1
*if,i,le,13,then
kk=1
*elseif,i,le,53,then
kk=2
*elseif,i,le,103,then
kk=3
*elseif,i,le,147,then
kk=4
*endif
*do,j,110,145,1
*get,vy,node,j+(i-1)*145,loc,z
r,kk+1,mk(1,kk,1)*abs(vy+0.4)*b0*h0/36
real,kk+1
e,j+(i-1)*145,100000+j+(i-1)*145
*enddo
*enddo
allsel
!底层弹簧单元
nsel,s,loc,z,-14.6
nsel,r,loc,x,0,0.21
nplot
ngen,2,200000,all,,,,,-0.2
type,3
r,10,mk(1,4,1)*15*b0*h0/36*150
real,10
*do,i,21172,21315,1
e,i,200000+i
*enddo
allsel
gplot
!上底层弹簧
nsel,s,,,110,145,1
nsel,a,,,21316,21351,1
nsel,a,,,21497,21532,1
nsel,a,,,21714,21749,1
nplot
ngen,2,300000,all,,,,,-0.4
type,3
r,10,mk(1,1,1)*0.4*b1*h1/36*8
real,10
*do,i,110,145,1
e,i,300000+i
*enddo
*do,i,21316,21351,1
e,i,300000+i
*enddo
*do,i,21497,21532,1
e,i,300000+i
*enddo
*do,i,21714,21749,1
e,i,300000+i
*enddo
allsel
!上侧面弹簧
nsel,s,,,21714,21749,1
nsel,a,,,21967,22002,1
nsel,a,,,22220,22255,1
nsel,a,,,22473,22508,1
nsel,a,,,22726,22761,1
nsel,a,,,22979,23014,1
nsel,a,,,23232,23267,1
nsel,a,,,23340,23375,1
nsel,a,,,23412,23447,1
ngen,2,400000,all,,,0.7
type,3
*do,j,21714,21749,1
*get,vy,node,j,loc,z
r,2,mk(1,1,1)*abs(vy-0.4)*b1*h1/36
real,2
e,j,400000+j
*enddo
*do,j,21967,22002,1
*get,vy,node,j,loc,z
r,2,mk(1,1,1)*abs(vy-0.4)*b1*h1/36
real,2
e,j,400000+j
*enddo
*do,j,22220,22255,1
*get,vy,node,j,loc,z
r,2,mk(1,1,1)*abs(vy-0.4)*b1*h1/36
real,2
e,j,400000+j
*enddo
*do,j,22473,22508,1
*get,vy,node,j,loc,z
r,2,mk(1,1,1)*abs(vy-0.4)*b1*h1/36
real,2
e,j,400000+j
*enddo
*do,j,22726,22761,1
*get,vy,node,j,loc,z
r,2,mk(1,1,1)*abs(vy-0.4)*b1*h1/36
real,2
e,j,400000+j
*enddo
*do,j,22979,23014,1
*get,vy,node,j,loc,z
r,2,mk(1,1,1)*abs(vy-0.4)*b1*h1/36
real,2
e,j,400000+j
*enddo
*do,j,23232,23267,1
*get,vy,node,j,loc,z
r,2,mk(1,1,1)*abs(vy-0.4)*b1*h1/36
real,2
e,j,400000+j
*enddo
*do,j,23340,23375,1
*get,vy,node,j,loc,z
r,2,mk(1,1,1)*abs(vy-0.4)*b1*h1/36
real,2
e,j,400000+j
*enddo
*do,j,23412,23447,1
*get,vy,node,j,loc,z
r,2,mk(1,1,1)*abs(vy-0.4)*b1*h1/36
real,2
e,j,400000+j
*enddo
allsel
wpro,,90.000000,
wpro,,90.000000,
csys,0
allsel
r,8
real,8
local,13,0,,,,270
type,4

!面1
*do,i,1,146,1
*do,j,110,118,1
nsel,s,node,,j+(i-1)*145
nsel,a,node,,j+1+(i-1)*145
nsel,a,node,,j+146+(i-1)*145
nsel,a,node,,j+145+(i-1)*145
esys,13
esurf
*enddo
*enddo

!面2

*do,i,1,146,1
*do,j,119,127,1
nsel,s,node,,j+(i-1)*145
nsel,a,node,,j+1+(i-1)*145
nsel,a,node,,j+146+(i-1)*145
nsel,a,node,,j+145+(i-1)*145
esys,13
esurf
*enddo
*enddo

!面3

*do,i,1,146,1
*do,j,128,136,1
nsel,s,node,,j+(i-1)*145
nsel,a,node,,j+1+(i-1)*145
nsel,a,node,,j+146+(i-1)*145
nsel,a,node,,j+145+(i-1)*145
esys,13
esurf
*enddo
*enddo

!面4

*do,i,1,146,1
*do,j,137,144,1
nsel,s,node,,j+(i-1)*145
nsel,a,node,,j+1+(i-1)*145
nsel,a,node,,j+146+(i-1)*145
nsel,a,node,,j+145+(i-1)*145
esys,13
esurf
*enddo
*enddo

!面5
allsel
*do,i,1,146,1
nsel,s,node,,110+(i-1)*145
nsel,a,node,,110+35+(i-1)*145
nsel,a,node,,110+180+(i-1)*145
nsel,a,node,,110+145+(i-1)*145
esys,13
esurf
*enddo

!面6

allsel

local,14,0,,,,270
*do,i,1,6,1
*do,j,21714,21722,1
nsel,s,node,,j+(i-1)*253
nsel,a,node,,j+1+(i-1)*253
nsel,a,node,,j+254+(i-1)*253
nsel,a,node,,j+253+(i-1)*253
esys,14
esurf
*enddo
*enddo

!面7

*do,i,1,6,1
*do,j,21723,21731,1
nsel,s,node,,j+(i-1)*253
nsel,a,node,,j+1+(i-1)*253
nsel,a,node,,j+254+(i-1)*253
nsel,a,node,,j+253+(i-1)*253
esys,14
esurf
*enddo
*enddo

!面8

*do,i,1,6,1
*do,j,21732,21740,1

nsel,s,node,,j+(i-1)*253
nsel,a,node,,j+1+(i-1)*253
nsel,a,node,,j+254+(i-1)*253
nsel,a,node,,j+253+(i-1)*253
esys,14
esurf
*enddo
*enddo

!面9

*do,i,1,6,1
*do,j,21741,21748,1
nsel,s,node,,j+(i-1)*253
nsel,a,node,,j+1+(i-1)*253
nsel,a,node,,j+254+(i-1)*253
nsel,a,node,,j+253+(i-1)*253
esys,14
esurf
*enddo
*enddo

!面10
allsel
*do,i,1,6,1
nsel,s,node,,21714+(i-1)*253
nsel,a,node,,21714+35+(i-1)*253
nsel,a,node,,21714+288+(i-1)*253
nsel,a,node,,21714+253+(i-1)*253
esys,14
esurf
*enddo
!面11
allsel
*do,j,23232,23266,1
nsel,s,node,,j
nsel,a,node,,j+1
nsel,a,node,,j+109
nsel,a,node,,j+108
esys,14
esurf
*enddo
!面12
allsel
*do,j,23340,23374,1
nsel,s,node,,j
nsel,a,node,,j+1
nsel,a,node,,j+73
nsel,a,node,,j+72
esys,14
esurf
*enddo
allsel
nsel,s,node,,23232
nsel,a,node,,23340
nsel,a,node,,23375
nsel,a,node,,23267
esys,14
esurf
allsel
nsel,s,node,,23340
nsel,a,node,,23412
nsel,a,node,,23447
nsel,a,node,,23375
esys,14
esurf

!边界条件
/solu
csys,5
allsel
!侧面约束
nsel,s,loc,x,1.05
nsel,r,loc,z,0,0.4
d,all,all
allsel
nsel,s,loc,x,0.6
nsel,r,loc,z,0,-14.6
d,all,all

!下底面约束
allsel
nsel,s,loc,z,-14.8
d,all,all

!上底面约束
allsel
nsel,s,loc,z,-0.4
nsel,r,loc,x,0.23,0.35
nplot
d,all,all
allsel
nsel,s,,,300110,300145,1
nplot
d,all,all

!桩下底节点的约束
nsel,s,loc,z,-14.6
nsel,r,loc,x,0,0.21
nplot
d,all,ux
d,all,uz

!桩上下底节点的约束

nsel,s,loc,z,0
nsel,r,loc,x,0.19,0.36
nplot
d,all,ux
d,all,uz
csys,0
ACEL,0,-9.81,0
allsel
eplot

！施加外部载荷

!桩顶表面平均载荷

allsel
csys,5
nsel,s,loc,z,0.4
nsel,r,loc,x,0,0.34
nsel,r,loc,y,0,360
nplot
esln,s
sfe,all,6,pres,,2000
sfe,all,14,pres,,2000
sfe,all,22,pres,,2000
sfe,all,30,pres,,2000
allsel
eplot

!土对桩的表面摩擦力
!5
csys,5
allsel
nsel,s,loc,z,-10.3,-14.6
esln,s
esel,r,type,,4
sfe,all,2,pres,,-FQ(1,5,1)
/psf,pres,TANx,2,0,1
eplot
!4
allsel
nsel,s,loc,z,-5.7,-10.1
esln,s
esel,r,type,,4
sfe,all,2,pres,,-FQ(1,4,1)
/psf,pres,TANx,2,0,1
eplot
!3
allsel
nsel,s,loc,z,-4.0,-5.5
esln,s
esel,r,type,,4
sfe,all,2,pres,,-FQ(1,3,1)
/psf,pres,TANx,2,0,1
eplot
!2
allsel
nsel,s,loc,z,-1.7,-3.8
esln,s
esel,r,type,,4
sfe,all,2,pres,,-FQ(1,2,1)
/psf,pres,TANx,2,0,1
eplot
!1
allsel
nsel,s,loc,z,0.4,-1.5
esln,s
esel,r,type,,4
sfe,all,2,pres,,-FQ(1,1,1)
/psf,pres,TANx,2,0,1
eplot
allsel
eplot

time,1
autot,on
nlgeom,on
nropt,full
csys,5
solve
finish