探讨关于粘弹性人工边界在ANSYS中的实现

最近几天一直在看粘弹性人工边界的理论，对比一篇论文做了一个实例，与大家一同分享，欢饮拍砖。

原文：基于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)返回的节点影响面积，这样才能使刚度系数和位移与弹性恢复力正比，使阻尼系数和速度与阻尼力正比。