前言
依据水工抗震规范中关于渡槽动水压力的部分编一个用于ANSYS渡槽模型动水压力施加的命令流，是我研究生时一直想要做的一件事，原因嘛主要是想对比一下规范提供的方法和ANSYS声学流体单元模拟水体这两种方法的结果。但是由于时间关系，毕业以前并没有完成，所幸毕业工作之后还有点空闲时间，就把这个命令编了一下。不过我估计几乎不会有人用这个命令了以及规范的这种方法了，毕竟从理论上来说，用声学单元模拟得到的结果肯定是比规范的方法精确的，而且操作也很简单，不涉及编命令。但是不管怎么样，还是放上来充个数吧。
命令依据NB 35047-2015《水电工程水工建筑物抗震设计规范》中附录B 渡槽槽体内动水压力计算部分编写，分为两个，一个用于矩形渡槽，另一个用于U型渡槽。
矩形渡槽的动水荷载命令说明
命令运行前的准备工作
保证模型的坐标系为：Y方向为重力方向，X方向为横槽向，Z方向为纵槽向；
- 计算中需要用到ANSYS重启动技术，因为施加的部分动水荷载的数值需要由渡槽截面槽底中心处的水平或竖直加速度响应值计算得到，所以命令中加入了导入地震波的命令，需要提前准备好地震波文件accel_x.txt、accel_y.txt、accel_z.txt并记得修改步长DT与总步数NT；
- 需先将槽底面和槽壁面编为组元：bottom_areas, left_areas,right_areas
修改相关参数（红色）。
pi=3.1415926 !圆周率
H_water=2 !槽内水深
l_water=2 !半槽宽
rou_water=1000 !水体密度
M_water=H_water*l_water*rou_water*2 ! 1m长度的水体质量
node_control=158 !控制点------截面槽底中心处的节点号
coord_y_control=ny(node_control) !得到控制点的y坐标
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
*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
*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/l_water*ARNODE(nnum)
type,99
nr=nr+1
r,nr,0,0,m_wv
real,nr
e,nnum
*enddo
!附加渡槽左壁面节点质量（横向地震作用冲击动水压力）
*if,(H_water/l_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/l_water*(coord_y/H_water+0.5*(coord_y/H_water)**2)
m_wh=m_wh*sqrt(3)*tanh(sqrt(3)*l_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/l_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/l_water*(coord_y/H_water+0.5*(coord_y/H_water)**2)
m_wh=m_wh*sqrt(3)*tanh(sqrt(3)*l_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=sqrt(5/2)*H_water/l_water
M1=2*rou_water*H_water*l_water*(sqrt(5/2)*l_water/H_water*tanh(temp)/3)
K1=M1*g_/l_water*sqrt(5/2)*tanh(temp)
h1=H_water*(1-(cosh(temp)-2)/temp/sinh(temp))
!寻找按高度施加弹簧的节点
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_control,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
type,99
nr=nr+1
r,nr,M1*0.5*(dist2+dist1),0,0
real,nr
e,node_max
d, node_max,uy
d, node_max,uz
!创建弹簧（模拟对流动水压力）
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) !在前面选出的节点中寻找与nnum节点最接近的节点号，并赋值给nnum_near
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,NT
!从结果提取控制点加速度
/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,1,0
parres,new,aqueduct_parameter,sav !恢复参数
!依据渡槽宽高比施加槽底动水压力荷载（横向地震作用下槽底的冲击动水压力）
*if,(H_water/l_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/l_water*a_wh
p_bh=p_bh*0.5*sqrt(3)*sinh(sqrt(3)*coord_x_pbh/H_water)/cosh(sqrt(3)*l_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/l_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/l_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/l_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
!**********************施加地震加速度和面荷载形式的动水压力并计算（结束）******************************