U型渡槽动水压力荷载施加命令及说明

• 程序中需要用到ANSYS重启动，因为需提取前一步加速度结果以施加部分动水压力；
• 默认Y方向为重力方向，X方向为横槽向，Z方向为纵槽向；
• 需准备地震波文件；
• 需先将槽底面和槽壁面编为组元：bottom_areas,left_areas,right_areas，U型渡槽下部曲线在水平向既是槽壁，在竖向又是槽底；

pi=3.1415926

H_water=3.6   !槽内水深

R_water=2.3    !半槽宽

rou_water=1000   !水体密度

M_water=rou_water*(2*R_water*(H_water-R_water)+0.5*pi*R_water**2)

!控制点------截面槽底中心处节点
node_control=8450
coord_y_control=ny(node_control)
coord_x_control=nx(node_control)
coord_z_control=nz(node_control)
coord_y_water=coord_y_control+H_water

/prep7
!定义mass21单元
et,99,mass21
!定义弹簧单元
et,98,combin14
!初始实常数
nr=100

!读入地震波
*SET,NT,200      !地震加速度步数
*SET,DT,0.01      !地震加速度步长
*DIM,accel_x,,NT
*DIM,accel_y,,NT
*DIM,accel_z,,NT
*VREAD,accel_x,'accel_x','txt'
(F12.9)
*VREAD,accel_y,'accel_y','txt'
(F12.9)
*VREAD,accel_z,'accel_z','txt'
(F12.9)

!***********************************附加质量（开始）**********************************

!将附加质量的渡槽底面的所有节点编为一个数组

allsel

CMSEL,S,bottom_areas,AREA

NSLA,S,1

*get,n_bottom,node,,count

*dim,node_bottom,array,n_bottom

*get,nmin,node,,num,min

node_bottom(1)=nmin

*do,i,2,n_bottom

*get,nnum,node,nmin,nxth

nmin=nnum

node_bottom(i)=nmin

*enddo

allsel

!将附加质量的渡槽左壁面的所有节点编为一个数组

CMSEL,S,left_areas,AREA

NSLA,S,1

nsel,U,loc,x,coord_x_control-0.001,coord_x_control+0.001

*get,n_left,node,,count

*dim,node_left,array,n_left

*get,nmin,node,,num,min

node_left(1)=nmin

*do,i,2,n_left

*get,nnum,node,nmin,nxth

nmin=nnum

node_left(i)=nmin

*enddo

allsel

!将附加质量的渡槽右壁面的所有节点编为一个数组

CMSEL,S,right_areas,AREA

NSLA,S,1

nsel,U,loc,x,coord_x_control-0.001,coord_x_control+0.001

*get,n_right,node,,count

*dim,node_right,array,n_right

*get,nmin,node,,num,min

node_right(1)=nmin

*do,i,2,n_right

*get,nnum,node,nmin,nxth

nmin=nnum

node_right(i)=nmin

*enddo

allsel

!附加渡槽底面节点质量（竖向地震作用冲击动水压力）

*do,i,1,n_bottom

nnum=node_bottom(i)

m_wv=0.4*M_water/R_water*ARNODE(nnum)

type,99

nr=nr+1

r,nr,0,0,m_wv

real,nr

e,nnum

*enddo

!附加渡槽左壁面节点质量（横向地震作用冲击动水压力）

*if,(H_water/R_water),le,1.5,then

*do,i,1,n_left

nnum=node_left(i)

coord_y=ny(nnum)-coord_y_control

m_wh=0.5*M_water/R_water*(coord_y/H_water+0.5*(coord_y/H_water)**2)

m_wh=m_wh*sqrt(3)*tanh(sqrt(3)*R_water/H_water)*ARNODE(nnum)

!确保附加的节点位于水面以下

*if,coord_y,le,H_water,then

type,99

nr=nr+1

r,nr,m_wh,0,0

real,nr

e,nnum

*endif

*enddo

*else

*do,i,1,n_left

nnum=node_left(i)

m_wh=0.5*M_water/H_water*ARNODE(nnum)

coord_y=ny(nnum)-coord_y_control

!确保附加的节点位于水面以下

*if,coord_y,le,H_water,then

type,99

nr=nr+1

r,nr,m_wh,0,0

real,nr

e,nnum

*endif

*enddo

*endif

!附加渡槽右壁面节点质量（横向地震作用冲击动水压力）

*if,(H_water/R_water),lt,1.5,then

*do,i,1,n_right

nnum=node_right(i)

coord_y=ny(nnum)-coord_y_control

m_wh=0.5*M_water/R_water*(coord_y/H_water+0.5*(coord_y/H_water)**2)

m_wh=m_wh*sqrt(3)*tanh(sqrt(3)*R_water/H_water)*ARNODE(nnum)

!确保附加的节点位于水面以下

*if,coord_y,le,H_water,then

type,99

nr=nr+1

r,nr,m_wh,0,0

real,nr

e,nnum

*endif

*enddo

*else

*do,i,1,n_right

nnum=node_right(i)

m_wh=0.5*M_water/H_water*ARNODE(nnum)

coord_y=ny(nnum)-coord_y_control

!确保附加的节点位于水面以下

*if,coord_y,le,H_water,then

type,99

nr=nr+1

r,nr,m_wh,0,0

real,nr

e,nnum

*endif

*enddo

*endif

!***********************************附加质量（结束）**********************************

!***********************************施加弹簧（开始）**********************************

g_=9.81

temp=(H_water-R_water)/R_water

M1=M_water*(0.571-1.276/((1+temp)**0.627)*(tanh(0.331*temp))**0.932)

w1_quadratic=1.323+0.228*(tanh(1.505*temp))**0.768-0.105*(tanh(1.105*temp))**4.659

w1_quadratic=w1_quadratic*g_/R_water

K1=M1*w1_quadratic

h1=H_water*(1-(temp**0.664)*(0.394+0.097*sinh(1.534*temp))/cosh(1.534*temp))

!计算槽壁施加弹簧点的x坐标

*if,h1,ge,R_water,then

coord_x_spring_0=R_water

*else

coord_x_spring_0=sqrt(R_water**2+(R_water-h1)**2)

*endif

!寻找按高度施加弹簧的节点

coord_y_spring_0=coord_y_control+h1

allsel

CMSEL,S,left_areas,AREA

CMSEL,A,right_areas,AREA

NSLA,S,1

node_find=node(coord_x_spring_0,coord_y_spring_0,coord_z_control)

coord_y_spring_1=ny(node_find)

NSEL,R,loc,Y,coord_y_spring_1-0.001,coord_y_spring_1+0.001

*get,n_spring,node,,count

*dim,node_spring,array,n_spring

*do,i,1,n_spring

nmin=node(0,0,100000)!选择距离(0,0,100000)最近的点

node_spring(i)=nmin

NSEL,U,node,,nmin!去掉刚选出的节点，重新选择

*enddo

allsel

*get,node_max,node,,num,max

dist1=0

*do,i,1,n_spring,2

nnum1=node_spring(i)

node1_x=nx(nnum1)

node1_y=ny(nnum1)

node1_z=nz(nnum1)

nnum2=node_spring(i+1)

node2_x=nx(nnum2)

*if,i,eq,n_spring-1,then

dist2=0

*else

nnum3=node_spring(i+2)

node3_z=nz(nnum3)

dist2=node1_z-node3_z

*endif

!创建中间节点

node_max=node_max+1

n,node_max,(node1_x+node2_x)*0.5,node1_y,node1_z

d,node_max,uy

d,node_max,uz

type,99

nr=nr+1

r,nr,M1*0.5*(dist2+dist1),0,0

real,nr

e,node_max

!创建弹簧（模拟对流动水压力）

type,98

nr=nr+1

r,nr,K1*0.5*(dist2+dist1)

real,nr

e,nnum1,node_max

e,node_max,nnum2

dist1=dist2

*enddo

!***********************************施加弹簧（结束）**********************************

!**********************施加地震加速度和面荷载形式的动水压力（开始）******************************

!查找槽底槽壁节点对应的渡槽截面中间点

allsel

nsel,s,loc,x,coord_x_control-0.001,coord_x_control+0.001

nsel,r,loc,y,coord_y_control-0.001,coord_y_control+0.001

*dim,node_bottom_near,array,n_bottom,2

*do,i,1,n_bottom

nnum=node_bottom(i)

nnum_near=nnear(nnum)

node_bottom_near(i,1)=nnum_near

*enddo

*dim,node_left_near,array,n_left,2

*do,i,1,n_left

nnum=node_left(i)

nnum_near=nnear(nnum)

node_left_near(i,1)=nnum_near

*enddo

*dim,node_right_near,array,n_right,2

*do,i,1,n_right

nnum=node_right(i)

nnum_near=nnear(nnum)

node_right_near(i,1)=nnum_near

*enddo

/solu

ANTYPE,TRANS

TRNOPT,FULL

OUTRES,ALL,ALL

allsel

!第一步计算

t=1

time,DT*t

acel,accel_x(t),accel_y(t),accel_z(t)

kbc,0

nsubst,1

rescontrol,DEFINE,last,last !最后时间步重启

parsav,all,aqueduct_parameter,sav !参数保存

solve

!rescontrol,file_summary

finish

*do,t,2,5

!从结果提取控制点加速度

/POST1

SET,t-1,1

*do,i,1,n_bottom

*get,node_bottom_near(i,2),node,node_bottom_near(i,1),a,x

*enddo

*do,i,1,n_left

*get,node_left_near(i,2),node,node_left_near(i,1),a,y

*enddo

*do,i,1,n_right

*get,node_right_near(i,2),node,node_right_near(i,1),a,y

*enddo

parsav,all,aqueduct_parameter,sav !参数保存

!再次进入求解模块

/solu

ANTYPE,,REST,t-1,last,0

parres,new,aqueduct_parameter,sav !恢复参数

!依据渡槽宽高比施加槽底动水压力荷载（横向地震作用下槽底的冲击动水压力）

*if,(H_water/R_water),lt,1.5,then

*do,i,1,n_bottom

a_wh=node_bottom_near(i,2)

nnum=node_bottom(i)

coord_x_pbh=coord_x_control-nx(nnum)

p_bh=0.5*M_water/R_water*a_wh

p_bh=p_bh*0.5*sqrt(3)*sinh(sqrt(3)*coord_x_pbh/H_water)/cosh(sqrt(3)*R_water/H_water)

F,nnum,FY,p_bh*ARNODE(nnum)

*enddo

*else

*do,i,1,n_bottom

nnum=node_bottom(i)

coord_x_pbh=coord_x_control-nx(nnum)

p_bh=0.5*M_water/H_water/R_water*coord_x_pbh

F,nnum,FY,p_bh*ARNODE(nnum)

*enddo

*endif

!左槽壁（竖向地震槽壁动水压力）

sign_left=(nx(node_left(1))-coord_x_control)/abs(nx(node_left(1))-coord_x_control)

*do,i,1,n_left

a_wv=node_left_near(i,2)

nnum=node_left(i)

coord_y=ny(nnum)-coord_y_control

*if,coord_y,le,H_water,then

p_wv=0.4*M_water/R_water*a_wv*cos(pi*0.5*coord_y/H_water)

F,nnum,FX,p_wv*ARNODE(nnum)*sign_left

*endif

*enddo

!右槽壁（竖向地震槽壁动水压力）

sign_right=(nx(node_right(1))-coord_x_control)/abs(nx(node_right(1))-coord_x_control)

*do,i,1,n_right

a_wv=node_right_near(i,2)

nnum=node_right(i)

coord_y=ny(nnum)-coord_y_control

*if,coord_y,le,H_water,then

p_wv=0.4*M_water/R_water*a_wv*cos(pi*0.5*coord_y/H_water)

F,nnum,FX,p_wv*ARNODE(nnum)*sign_right

*endif

*enddo

time,DT*t

acel,accel_x(t),accel_y(t),accel_z(t)

nsubst,1

allsel

solve

finish

*enddo

!**********************施加地震加速度和面荷载形式的动水压力（结束）******************************