最近几天一直在看粘弹性人工边界的理论，对比一篇论文做了一个实例，与大家一同分享，欢饮拍砖。原文：基于ABAQUS的粘弹性边界单元及在重力坝抗震分析中的应用本帖：ANSYS算例程序如下：/PREP7!*******************************************************************************L=800 !水平长度
原文：基于ABAQUS的粘弹性边界单元及在重力坝抗震分析中的应用
本帖：ANSYS算例
程序如下：
/PREP7
!*******************************************************************************
L=800 !水平长度
H=400 !竖向深度
E=1.323E+10 !弹性模量
density=2.7E+03 !密度
mu=0.25 !泊松比
dxyz=20 !网格尺寸
G = E/(2.*(1.+mu)) !剪切模量
alfa = E*(1-mu)/((1.+mu)*(1.-2.*mu)) !若计算平面应力，此式需要修改
Cp=sqrt(alfa/density) !压缩波速
Cs=sqrt(G/density) !剪切波速
R=(sqrt(L*L/4.+H*H)+L/2+H)/3 !波源到边界点等效长度
KbT=0.5*G/R
KbN=1.0*G/R
CbT=density*Cs
CbN=density*Cp
!*******************************************************************************
ET, 1, plane42,,,2 !按平面应变计算
et, 2, combin14, ,, 2 !切向
et, 3, combin14, ,, 2 !法向
r, 2, KbT, CbT
r, 3, KbN, CbN
MP, EX, 1, E
MP, PRXY, 1, mu
MP, DENS, 1, density
rectng,-L/2.,L/2,0.,H
asel, all
aesize, all, dxyz
mshape,0,2D
mshkey,1
amesh, all
!*******************************************************************************
PI = Acos(-1)
TNUM = 400
btime = 0.005
etime = 2.00
DTIME = 0.005
!****************************
nsel,s,loc,y,0 !选中需要加等效荷载的底边界点
*get,NBOTM,node,,count !得到选中的结点数，存入NBOTM
*dim,BOTMU,array,TNUM,3, , , , !储存底边界位移和速度
*dim,BOTMF,array,TNUM, , , , , !储存底边界等效节点荷载
*do,T,1,100,1
BOTMU(T,1)=SIN(4.*PI*T*DTIME)-0.5*SIN(8.*PI*T*DTIME)
BOTMU(T,2)=4.*PI*COS(4.*PI*T*DTIME)-4.*PI*COS(8.*PI*T*DTIME)
BOTMU(T,3)=-16.*PI*PI*SIN(4.*PI*T*DTIME)+32.*PI*PI*SIN(8.*PI*T*DTIME)
*enddo
*do,T,101,TNUM,1
BOTMU(T,1)=0
BOTMU(T,2)=0
BOTMU(T,3)=0
*enddo
*do,T,1,TNUM,1
BOTMF(T) = 1.*(KbT*BOTMU(T,1)+CbT*BOTMU(T,2)+density*Cs*BOTMU(T,2))
*enddo
!****************************
nsel,s,loc,x,-L/2 !选择左边界上的全部节点
*get,NNODE,node,,count !求得已选中节点的数量
*dim,DTIME1,array,NNODE,,,,, !储存入射波从底边界到侧边界任意节点的传播时间
*dim,DTIME2,array,NNODE,,,,, !储存反射波从地表到侧边界任意节点的传播时间
*dim,VERTU1,array,TNUM,NNODE,,,, !储存入射波从底边界传播到侧边界任意节点而引起的位移
*dim,VERTV1,array,TNUM,NNODE,,,, !储存入射波从底边界传播到侧边界任意节点而引起的速度
*dim,VERTU2,array,TNUM,NNODE,,,, !储存反射波从地表传播到侧边界任意节点而引起的位移
*dim,VERTV2,array,TNUM,NNODE,,,, !储存反射波从地表传播到侧边界任意节点而引起的速度
*dim,VERTFN,array,TNUM,NNODE,,,, !储存侧边界的法向等效荷载
*dim,VERTFT,array,TNUM,NNODE,,,, !储存侧边界的切向等效荷载
*do,inum,2,NNODE,1
DTIME1(INUM)=((INUM-1.)*dxyz)/Cs
DTIME2(INUM)=(2.*H-(INUM-1.)*dxyz)/Cs
*enddo
!*
*DO,T,1,TNUM,1
*do,inum,2,NNODE,1
ITIME1 = T*DTIME-DTIME1(INUM)
ITIME2 = T*DTIME-DTIME2(INUM)
*if,ITIME1,LT,0.,THEN
VERTU1(T,inum) = 0.
VERTV1(T,inum) = 0.
*elseif,ITIME1,GE,0.,and,ITIME1,LE,0.5,THEN
VERTU1(T,inum) = SIN(4.*PI*ITIME1)-0.5*SIN(8.*PI*ITIME1)
VERTV1(T,inum) = 4.*PI*COS(4.*PI*ITIME1)-4.*PI*COS(8.*PI*ITIME1)
*else
VERTU1(T,inum) = 0
VERTV1(T,inum) = 0
*endif
*if,ITIME2,LT,0.,THEN
VERTU2(T,inum) = 0.
VERTV2(T,inum) = 0.
*elseif,ITIME2,GE,0.,and,ITIME2,LE,0.5,THEN
VERTU2(T,inum) = SIN(4.*PI*ITIME2)-0.5*SIN(8.*PI*ITIME2)
VERTV2(T,inum) = 4.*PI*COS(4.*PI*ITIME2)-4.*PI*COS(8.*PI*ITIME2)
*else
VERTU2(T,inum) = 0
VERTV2(T,inum) = 0
*endif
*enddo
*ENDDO
!*
*DO,T,1,TNUM,1
*do,inum,2,NNODE,1
VERTFN(T,inum) = 1.*(KbN*VERTU1(T,inum)+CbN*VERTV1(T,inum)+KbN*VERTU2(T,inum)+CbN*VERTV2(T,inum))
VERTFT(T,inum) = 1.*CbT*(VERTV1(T,inum)-VERTV2(T,inum))
*enddo
*ENDDO
!*******************************************************************************
!以下建立底边界法向和切向弹簧阻尼单元
nsel,s,loc,y,0.
*get,np,node,,count !得到选中的结点数，存入np
*get,npmax,node,,num,maxd !得到已经定义的最大结点数，存入npmax
*do,ip,1,np
npnum=node((ip-1)*dxyz-L/2.,0.,0.)
x=nx(npnum)
y=ny(npnum)
z=nz(npnum)
npmax=npmax+1
n,npmax,x,y-dxyz/20,z !定义底边界法向结点以便与边界点形成法向单元
type,3
real,3
e,npnum,npmax
d,npmax,all,0. !约束新生成的点
npmax=npmax+1
n,npmax,x-dxyz/20,y,z !定义底边界切向结点以便与边界点形成切向单元
type,2
real,2
e,npnum,npmax
d,npmax,all,0. !约束新生成的点
*enddo
!以下建立左边界法向和切向弹簧阻尼单元
nsel,s,loc,x,-L/2
*get,np,node,,count !得到选中的结点数，存入np
*get,npmax,node,,num,maxd !得到已经定义的最大结点数，存入npmax
*do,ip,2,np !侧边界最下面一个点按底边界上处理
npnum=node(-L/2,(ip-1)*dxyz,0.)
x=nx(npnum)
y=ny(npnum)
z=nz(npnum)
npmax=npmax+1
n,npmax,x-dxyz/20,y,z !定义左边界法向结点以便与边界点形成法向单元
type,3
real,3
e,npnum,npmax
d,npmax,all,0. !约束新生成的点
npmax=npmax+1
n,npmax,x,y-dxyz/20,z !定义左边界切向结点以便与边界点形成切向单元
type,2
real,2
e,npnum,npmax
d,npmax,all,0. !约束新生成的点
*enddo
!以下建立右边界法向和切向弹簧阻尼单元
nsel,s,loc,x,L/2
*get,np,node,,count !得到选中的结点数，存入np
*get,npmax,node,,num,maxd !得到已经定义的最大结点数，存入npmax
*do,ip,2,np !侧边界最下面一个点按底边界上处理
npnum=node(L/2,(ip-1)*dxyz,0.)
x=nx(npnum)
y=ny(npnum)
z=nz(npnum)
npmax=npmax+1
n,npmax,x+dxyz/20,y,z !定义右边界法向结点以便与边界点形成法向单元
type,3
real,3
e,npnum,npmax
d,npmax,all,0. !约束新生成的点
npmax=npmax+1
n,npmax,x,y-dxyz/20,z !定义右边界切向结点以便与边界点形成切向单元
type,2
real,2
e,npnum,npmax
d,npmax,all,0. !约束新生成的点
*enddo
allsel,all
/pnum,type,1
/number,1
eplot
finish
!*******************************************************************************
/solu
ANTYPE,trans
TRNOPT,FULL
LUMPM,0
*do,T,1,TNUM,1
Time,T*Dtime
!****************************
*do,inum,1,NBOTM,1 !选中需要加等效荷载的底边界点
npnum=node(-L/2+(inum-1)*dxyz,0.,0.)
F,NPNUM,FX,ARNODE(npnum)*BOTMF(T)
*enddo
!****************************
*do,inum,2,NNODE,1 !选中需要加等效荷载的左边界点
npnum=node(-L/2,(inum-1)*dxyz,0.)
F,NPNUM,FX,ARNODE(npnum)*VERTFN(T,inum)
F,NPNUM,FY,ARNODE(npnum)*VERTFT(T,inum)
*enddo
!****************************
*do,inum,2,NNODE,1 !选中需要加等效荷载的右边界点
npnum=node(L/2,(inum-1)*dxyz,0.)
F,NPNUM,FX,ARNODE(npnum)*VERTFN(T,inum)
F,NPNUM,FY,-1.*ARNODE(npnum)*VERTFT(T,inum)
*enddo
allsel,all
SOLVE
*ENDDO 后处理文件：
finish
/POST26 !进入时间历程处理器
!-----------------选取节点定义变量--------------------------
NA=node(-L/2,H,0.) !选择左边界顶部节点
NB=node(0,H,0.) !选择上边界中间节点
NC=node(L/2,H,0.) !选择右边界顶部节点
ND=node(L/2,0,0.) !选择右边界底部节点
NE=node(0,0,0.) !选择底边界中间节点
NF=node(L/2,0,0.) !选择左边界底部节点
!-----------------提取位移----------------------------------
/AXLAB,X,Time(s) !设置X轴标识
/AXLAB,Y,Displacement Z (m) !设置Y轴标识
nsol,2,NA,u,x,NAux !提取底部左边端节点的uy向位移
plvar,2 !显示2号变量随时间变化曲线
/ui,copy,save,bmp,,,,port
nsol,3,NB,u,x,NBux
plvar,3
/ui,copy,save,bmp,,,,port
nsol,4,NC,u,x,NCux
plvar,4
/ui,copy,save,bmp,,,,port
nsol,5,ND,u,x,NDux
plvar,5
/ui,copy,save,bmp,,,,port
nsol,6,NE,u,x,NEux
plvar,6
/ui,copy,save,bmp,,,,port
nsol,7,NF,u,x,NFux
plvar,5,6,7
/ui,copy,save,bmp,,,,port
!-----------------提取数据---------------------------------
prvar,2,3,4,5,6,7
结果1
结果2
结果3
结果4
结果5
结果6
底部三个观察点对比：
激励位移时程：
本算例跟许多文献中算例取值是基本相同的，但是得出来的结果显然是有问题的，不知道是什么问题，仍然在查找原因。
哈哈，终于找到原因了，是一个很低级的错误。
这个错误就在于KbT,KbN,CbT,CbN的设置，两者都必须要在定义的时候乘以dxyz，或者通过ARNODE(NODENUM)返回的节点影响面积，这样才能使刚度系数和位移与弹性恢复力正比，使阻尼系数和速度与阻尼力正比。
