一个包括建模\分网\外程序调用优化\动力分析\模态提取的大型斜拉桥梁分析实例

!************************* ***桥有限元模型 ******************************

!**桥的整体有限元模型建立。该模型为鱼脊骨模型。先按照梁的自然节点（拉索位置）
!和主塔的自然节点（截面变化处、纵横梁交叉点）建立三维梁单元，然后分别读入前面
!建立的梁、塔的截面信息，生成空间整体三维有限元梁柱子体系。为了后面建立索的方便，
!索的两段的节点必须是按照等差数列的方式排列的。并且用参数获得索的第一个节点号、
!最后一个节点号，以及节点增加步长。然后使用循环语句分别生成八排索。

!**************************** 本程序的单位制 *******************************

!质量；吨； 力：牛顿； 长度：毫米； 时间：秒； 应力及弹性模量：兆帕
!密度：吨/立方毫米 能量：牛顿*毫米=0。001焦耳； 加速度：MM/S^-2=0.001g.

!****************************************************************************

finish
/clear,nostart
/filename,yhbridge,1
/title,yhbridge
/prep7

!NSKIP里存储的是本次调用在迭代序列中的序号，ANSYS朝FREQUENCE里保存数据时才用追加的形式保存
!而MATLAB朝OPTPARA里保存数据也是采用追加的方式保存。
!这样，当优化运行结束时，OPTPARA 里每行储存的数据就对应设计空间内一个相应的点，
!而FREQUENCE 则是该点对应的各阶频率 。

! NSLOOP 里按照逐行顺序，存储的是该次调用时各个循环指标的值（每个循环指标对应一个修正系数OPT）。
! NSLOOP 是按照刷新的方式存储的，所以里面的数据仅仅对应一个设计点（对应最后一次调用）。

!而NSKIP则是该次循环的序号，显然是各个指标的乘积。
*dim,NSLOOP,table,3,4
*tread,NSLOOP,NSLOOP,dat,'e:\ansyswork\yhdata\input'
NSKIP=1
*do,i,1,3,1
*do,j,1,4,1
NSKIP=NSKIP*NSLOOP(i,j)
*enddo
*enddo

! 由于是采用非追加的方式保存数据，所以，直接读出数据即可。
*dim,optpara,table,1,12,1
*tread,optpara,optpara,dat,'e:\ansyswork\yhdata\input',0
! opt1修正主梁弹性模量，opt2 修正主梁密度；opt3修正桥塔弹性模量，opt4修正桥塔密度；opt5修正索的弹性模量.
! pt6修正主塔垂直桥纵轴向水平约束；opt7 对应桥头水平约束的修正系数；opt8 对应辅助桥墩的水平约束修正系数。
!pt9 修正桥纵轴向水平约束；pt10 修正抗扭约束；pt11 修正4个桥墩竖向约束，pt12修正索幅在桥面上的宽度（桥面位置上下游索间距）

opt1=optpara(1,1)
opt2=optpara(1,2)
opt3=optpara(1,3)
opt4=optpara(1,4)
opt5=optpara(1,5)
opt6=optpara(1,6)
opt7=optpara(1,7)
opt8=optpara(1,8)
opt9=optpara(1,9)
opt10=optpara(1,10)
opt11=optpara(1,11)
opt12=optpara(1,12)

!定义主梁为第一类单元、实常数、材料组
et,1,189
mp,prxy,1,0.17
mp,ex,1,opt1*3.5e4
mp,dens,1,opt2*1.01*2.4e-9

!定义鱼刺（横向梁）为第二类单元 、实常数、材料组
et,2,4
mp,prxy,2,0.01
mp,ex,2,3.5e4
dens,2,0.002e-9

!定义桥塔为第三类单元、实常数、材料组
et,3,188
mp,prxy,3,0.17
mp,ex,3,opt3*0.95*3.5e4
dens,3,opt4*2.4e-9

!定义索为第四类单元。由于实常数、及弹性模量各个索有区别，所以在后面分别定义。
et,4,8

!定义主梁的弹性水平约束，采用 LINK 8 单元等效弹性约束
et,5,8
mp,prxy,5,0.3
mp,ex,5,2.0e5
dens,5,0.002e-9

!定义所有beam 4 单元的实常数。 注意实常数与单位有关。刚臂的实常数为R2。主梁当采用BEAM 4
!单元时，实在常数为R1
r,1,0.684e7,0.147e15,0.221e13,100,12000,,,,,,
r,2,1.0e6,0.33e14,0.33e14,100,100

!r,3,0.684e7,0,
!3 号实常数为主桥塔的垂直桥纵轴的水平弹性约束等效成拉压杆件的面积 。
!4 号实常数为两侧桥头的垂直桥纵轴的水平弹性约束等效成拉压杆件的面积 。
!5 号实常数为两侧辅助桥墩的垂直桥纵轴的水平弹性约束等效成拉压杆件的面积 。
r,3,opt6*52.5,0
r,4,opt7*1100,0
r,5,opt8*1100,0

!6 号实常数为两侧桥头的沿桥纵轴向弹性约束等效成拉压杆件的面积 。

r,6,opt9*1,0
r7=opt10*40000
r8=opt11*90000-2*r7

!7 号实常数为两侧桥头的外侧竖向抗扭转弹性约束等效成拉压杆件的面积 。
!7 号实常数为两侧辅助桥墩的外侧抗扭转弹性约束等效成拉压杆件的面积 。
r,7,r7,0
!8 号实常数为两侧辅助桥墩及桥头墩的中心竖向弹性约束等效成拉压杆件的面积 。
r,8,r8,0

!*************** 下面定义索的实常数 ***********************

!从CABLEDATA 里读出索的面积及初始应变。
!其中第一列为河侧索面积，单位为平方毫米。第二列为河侧索初始索力，单位为吨。第三列是岸侧索面积
!第四列是岸侧索初始索力。第五列是河侧索的水平投影，第六列是岸侧索的水平投影.当考虑索的非线性，
!采用 ERNST 修正法来修正索的弹性模量，要注意单位转化。设置索的容重gama，来作为后面计算索的弹性
!模量的数据。容重单位为KN/M^3,弹性模量单位为帕。
!定义数组 Eeqs 来存储后面生成的等效弹性模量。 第一列存储河侧索的等效弹性模量,第二列存储岸侧索的
!等效弹性模量。定义 COFF 来存储弹性模量等效时的系数。这样更容易看出是否出错。这个系数应该是略
!微小于1的正数。
!现在索的初始张拉力已经全部在后面用 CABLEF 来存储了，CABLEDATA 里的初始索力数据为**研究
!所的设计数据.

*dim,rcable,table,11,6,1
*tread,rcable,cabledata,txt,'e:\ansyswork\yhdata\input',

*dim,cablef,table,11,8,1
*tread,cablef,cablef,txt,'e:\ansyswork\yhdata\input',1

*set,gama,78000
*SET,e0,2.07e11

*dim,Eeqs,array,11,8,1
*dim,coff,array,11,8,1

!为了后面模拟施工过程的换索方便，将每根索都单独赋实常数。
!为了模拟有无初始索力的方便，在索的初始应变前都有个系数BEITA，
!当不考虑初始应变的时候，就给BEITA 赋 0 值。考虑，则为 1 。
!上游部分：
!定义北塔岸侧索的材料及实常数 。 编号从 11 开始到 21 结束
!定义北塔河侧索的材料及实常数 。 编号从 31 开始到 41 结束
!定义北塔河侧索的材料及实常数 。 编号从 51 开始到 61 结束
!定义南塔岸侧索的材料及实常数 。 编号从 71 开始到 81 结束
!下游部分：
!定义北塔岸侧索的材料及实常数 。 编号从 91 开始到 101 结束
!定义北塔河侧索的材料及实常数 。 编号从 111 开始到 121 结束
!定义北塔河侧索的材料及实常数 。 编号从 131 开始到 141 结束
!定义南塔岸侧索的材料及实常数 。 编号从 151 开始到 161 结束

*do,i,1,11,1
*set,coff1,(rcable(i,6)*0.001*gama)**2
*set,coff2,12*(cablef(i,1)*1.0e6/rcable(i,3))**3
coff(i,1)=1/((coff1*e0/coff2)+1)
Eeqs(i,1)=opt5*0.79*coff(i,1)*e0*1.0e-6
mp,prxy,i+10,0.3
mp,ex,i+10,Eeqs(i,1)
mp,dens,i+10,0.1e-11
mp,alpx,i+10,0.00001
*enddo

*do,i,1,11,1
*set,coff1,(rcable(i,5)*0.001*gama)**2
*set,coff2,12*(cablef(i,2)*1.0e6/rcable(i,1))**3
coff(i,2)=1/((coff1*e0/coff2)+1)
Eeqs(i,2)=opt5*0.79*coff(i,2)*e0*1.0e-6
mp,prxy,i+30,0.3
mp,ex,i+30,Eeqs(i,2)
mp,dens,i+30,0.1e-11
mp,alpx,i+30,0.00001
*enddo

*do,i,1,11,1
*set,coff1,(rcable(i,5)*0.001*gama)**2
*set,coff2,12*(cablef(i,3)*1.0e6/rcable(i,1))**3
coff(i,3)=1/((coff1*e0/coff2)+1)
Eeqs(i,3)=opt5*0.79*coff(i,3)*e0*1.0e-6
mp,prxy,i+50,0.3
mp,ex,i+50,Eeqs(i,3)
mp,dens,i+50,0.1e-11
mp,alpx,i+50,0.00001
*enddo

*do,i,1,11,1
*set,coff1,(rcable(i,6)*0.001*gama)**2
*set,coff2,12*(cablef(i,4)*1.0e6/rcable(i,3))**3
coff(i,4)=1/((coff1*e0/coff2)+1)
Eeqs(i,4)=opt5*0.79*coff(i,4)*e0*1.0e-6
mp,prxy,i+70,0.3
mp,ex,i+70,Eeqs(i,4)
mp,dens,i+70,0.1e-11
mp,alpx,i+70,0.00001
*enddo

*do,i,1,11,1
*set,coff1,(rcable(i,6)*0.001*gama)**2
*set,coff2,12*(cablef(i,5)*1.0e6/rcable(i,3))**3
coff(i,5)=1/((coff1*e0/coff2)+1)
Eeqs(i,5)=opt5*0.79*coff(i,5)*e0*1.0e-6
mp,prxy,i+90,0.3
mp,ex,i+90,Eeqs(i,5)
mp,dens,i+90,0.1e-11
mp,alpx,i+90,0.00001
*enddo

*do,i,1,11,1
*set,coff1,(rcable(i,5)*0.001*gama)**2
*set,coff2,12*(cablef(i,6)*1.0e6/rcable(i,1))**3
coff(i,6)=1/((coff1*e0/coff2)+1)
Eeqs(i,6)=opt5*0.79*coff(i,6)*e0*1.0e-6
mp,prxy,i+110,0.3
mp,ex,i+110,Eeqs(i,6)
mp,dens,i+110,0.1e-11
mp,alpx,i+110,0.00001
*enddo

*do,i,1,11,1
*set,coff1,(rcable(i,5)*0.001*gama)**2
*set,coff2,12*(cablef(i,7)*1.0e6/rcable(i,1))**3
coff(i,7)=1/((coff1*e0/coff2)+1)
Eeqs(i,7)=opt5*0.79*coff(i,7)*e0*1.0e-6
mp,prxy,i+130,0.3
mp,ex,i+130,Eeqs(i,7)
mp,dens,i+130,0.1e-11
mp,alpx,i+130,0.00001
*enddo

*do,i,1,11,1
*set,coff1,(rcable(i,6)*0.001*gama)**2
*set,coff2,12*(cablef(i,8)*1.0e6/rcable(i,3))**3
coff(i,8)=1/((coff1*e0/coff2)+1)
Eeqs(i,8)=opt5*0.79*coff(i,8)*e0*1.0e-6
mp,prxy,i+150,0.3
mp,ex,i+150,Eeqs(i,8)
mp,dens,i+150,0.1e-11
mp,alpx,i+150,0.00001
*enddo

!为了后面模拟施工过程的换索方便，将每根索都单独赋实常数。
!为了模拟有无初始索力的方便，在索的初始应变前都有个系数BEITA，
!当不考虑初始应变的时候，就给BEITA 赋 0 值。考虑，则为 1 。
!再赋实常数。
beita=1
*do,i,1,11,1
R,i+10,rcable(i,3),beita*cablef(i,1)/rcable(i,3)/Eeqs(i,1)
R,i+30,rcable(i,1),beita*cablef(i,2)/rcable(i,1)/Eeqs(i,2)
R,i+50,rcable(i,1),beita*cablef(i,3)/rcable(i,1)/Eeqs(i,3)
R,i+70,rcable(i,3),beita*cablef(i,4)/rcable(i,3)/Eeqs(i,4)
R,i+90,rcable(i,3),beita*cablef(i,5)/rcable(i,3)/Eeqs(i,5)
R,i+110,rcable(i,1),beita*cablef(i,6)/rcable(i,1)/Eeqs(i,6)
R,i+130,rcable(i,1),beita*cablef(i,7)/rcable(i,1)/Eeqs(i,7)
R,i+150,rcable(i,3),beita*cablef(i,8)/rcable(i,3)/Eeqs(i,8)
*enddo

!下面开始定义塔的几何模型。坐标起点从对称纵剖面、横剖面与基础承台的交点为坐标原点。
!下面定义塔脚坐标
k,1,8620,5000
k,2,-8620,5000

!下面定义塔第一横梁中心坐标。
k,3,6250,44758
k,4,-6250,44758
!下面定义两侧双排索从第一个索到第11个索的坐标。
k,5,6250,46100
k,6,-6250,46100

k,7,6250,48850
k,8,-6250,48850
k,9,6250,50450
k,10,-6250,50450
k,11,6250,51750
k,12,-6250,51750
k,13,6250,53000
k,14,-6250,53000
k,15,6250,54200
k,16,-6250,54200
k,17,6250,55400
k,18,-6250,55400
k,19,6250,56600
k,20,-6250,56600
k,21,6250,57800
k,22,-6250,57800
k,23,6250,59000
k,24,-6250,59000
!下面定义塔第二个横梁的坐标。
k,25,6250,59625
k,26,-6250,59625
!下面定义塔顶坐标。
k,27,6250,65000
k,28,-6250,65000

!下面定义塔截面的方向关键点。
k,31,6250,44758,-10000
k,32,-6250,44758,-10000
k,33,6250,44758,-10000
k,34,6250,59625,-10000

!由于1、2两根短索接近垂直，与塔的夹角很小。直接将其朝轴线上简化误差比较大。
!为次，在该截面建立两个刚臂来解决。下面定义刚臂端点节点的关键点。

!北塔上游的刚臂
k,50,-6250,46000,126200-1200
k,51,-6250,46000,126200+1200
!北塔下游的刚臂
k,52,6250,46000,126200-1200
k,53,6250,46000,126200+1200
!南塔上游的刚臂
k,54,-6250,46000,386200-1200
k,55,-6250,46000,386200+1200
!南塔下游的刚臂
k,56,6250,46000,386200-1200
k,57,6250,46000,386200+1200

cm,ktower,kp
cmsel,s,ktower
!复制生成北桥塔的关键点，编号增加量为100，以便于控制。
kgen,2,1,34,1,0,0,126200,100,0,0
!移动生成南桥塔，编号增加量再增加100，以便于控制。
kgen,2,1,34,1,0,0,386200,0,0,1

!下面将塔的关键点用线连接起来，生成桥梁主塔的中心线。

*do,i,1,25,2
l,i,i+2
l,i+1,i+3
*enddo
l,3,4
l,25,26

numstr,line,101
*do,i,101,125,2
l,i,i+2
l,i+1,i+3
*enddo
l,103,104
l,125,126

!整个南塔北塔已经建模完毕。现在开始建立主梁模型。主粮模型按照单梁鱼脊骨模型
!建立。整个梁自然节点的Z坐标已经保存在数据文件 beamdata.txt 里面了。该数据的格式为：
!第一行和第一列分别表示行标号、列标号；按照列顺序分别为X、Y、Z坐标。读入该数据到
!数组参数 axis ，再用beamb、bdeep 两个标量参数来表示等效的索横向间距（桥面计算宽度）
!和斜拉索偏离主梁形心的间距，以考虑斜拉索建立的结构预应力对主梁产生的反拱。

*set,beamb,12000*opt12
*set,bdeep,368.25

!定义 axis 数组的维数。由于主梁在该计算模型中被分为51个节点，所以该数组
!为51*3的数组。

*dim,axis,table,51,3,1
*tread,axis,beamdata,txt,'e:\ansyswork\yhdata\input',

!以上 axis 里的是桥面的坐标及高程数据，但是实际上主梁的中性轴线距离桥面距离为613.754毫米，所以
!还需要把 axis 的第二列，也就是Y坐标向下平移这个距离，使得 axis 里是纵轴所在的坐标。

*do,i,1,11
axis(i,2)=axis(i,2)-613.754
*enddo

!定义cabup 数组的维数。该数组保存的是上游索的下节点位置。只保存Z坐标及Y坐标。X坐标由半幅梁宽得到。
!为51*3的数组。
*dim,cabup,table,51,3,1
cabup(0,0)=0
cabup(0,1)=1
cabup(0,2)=2
cabup(0,3)=3
*do,i,1,51,1
cabup(i,0)=i
cabup(i,1)=-1*beamb/2
cabup(i,2)=axis(i,2)-bdeep
cabup(i,3)=axis(i,3)
*enddo

!定义cabdown 数组的维数。该数组保存的是下游索的下节点位置。 只保存Z坐标及Y坐标。X坐标由半幅梁宽得到。
*dim,cabdown,table,51,3,1
cabdown(0,0)=0
cabdown(0,1)=1
cabdown(0,2)=2
cabdown(0,3)=3
*do,i,1,51,1
cabdown(i,0)=i
cabdown(i,1)=1*beamb/2
cabdown(i,2)=axis(i,2)-bdeep
cabdown(i,3)=axis(i,3)
*enddo
!定义bnup 数组的维数。该数组保存的是鱼脊椎上游方向节点位置。 只保存Z坐标及Y坐标。X坐标由半幅梁宽得到。
*dim,bnup,table,51,3,1
bnup(0,0)=0
bnup(0,1)=1
bnup(0,2)=2
bnup(0,3)=3
*do,i,1,51,1
bnup(i,0)=i
bnup(i,1)=-1*beamb/2
bnup(i,2)=axis(i,2)
bnup(i,3)=axis(i,3)
*enddo
!定义bndown 数组的维数。该数组保存的是鱼脊椎上游方向节点位置。 只保存Z坐标及Y坐标。X坐标由半幅梁宽得到。
*dim,bndown,table,51,3,1
bndown(0,0)=0
bndown(0,1)=1
bndown(0,2)=2
bndown(0,3)=3
*do,i,1,51,1
bndown(i,0)=i
bndown(i,1)=1*beamb/2
bndown(i,2)=axis(i,2)
bndown(i,3)=axis(i,3)
*enddo
!上面已经建立了梁的全部关键节点坐标。下面将利用上面各个参数数组里保存的
!相关坐标数据来生成鱼脊椎梁模型的线。生成的节点将从301开始编号。 300号节点
!将用来作为主梁的截面方向关键点。

!生成梁的轴线上关键节点。方向关键点K300。

k,300,0,30000,20000
*do,i,1,51,1
k,i+300,axis(i,1),axis(i,2),axis(i,3)
*enddo
!***********************************************************************
!此处采用 LINK 8 单元 生成所有构成边界条件的弹性约束。

!生成梁在桥墩处的垂直轴线方向的水平弹性约束，采用 LINK8 单元。
!spring 是边界的弹性杆长度，为了美观可以调节该参数。
spring=1000
! 355，356 ，357，358四个点分别在四个桥墩位置的正下方，用来生成竖向约束（模拟盆式支座）
k,355,axis(1,1),axis(1,2)-spring,axis(1,3)
k,356,axis(5,1),axis(5,2)-spring,axis(5,3)
k,357,axis(47,1),axis(47,2)-spring,axis(47,3)
k,358,axis(51,1),axis(51,2)-spring,axis(51,3)
! 361，362 ，363，364四个点分别在四个桥墩位置的侧向，用来生成水平约束（模拟盆式支座的摩擦力约束）
k,361,axis(1,1)+spring,axis(1,2),axis(1,3)
k,362,axis(5,1)+spring,axis(5,2),axis(5,3)
k,363,axis(47,1)+spring,axis(47,2),axis(47,3)
k,364,axis(51,1)+spring,axis(51,2),axis(51,3)

!生成梁在桥墩处的沿桥轴线方向的水平弹性约束，采用 LINK8 单元，用来生成水平约束（模拟盆式支座的摩擦力约束）。
k,371,axis(1,1),axis(1,2),axis(1,3)-spring
!生成梁在主塔处的沿桥垂直轴线方向的水平弹性约束，采用 LINK8 单元，用来生成水平约束（模拟桥塔的碰撞）。
k,372,axis(14,1)+spring,axis(14,2),axis(14,3)
k,373,axis(38,1)+spring,axis(38,2),axis(38,3)

!梁在桥墩处边上的的垂直方向的抗扭转弹性约束，采用 LINK8 单元。但是不需要再次建立线
!因为后面建立索端刚臂时在编号含有1，5，47，51 处时，恰好不是索端而是桥墩位置，在后
!面划分网格的时候注意把这些打算作为边界条件弹性约束的线挑选出来划分成边界条件的单元
!而非索端刚臂！
!***********************************************************************

! 生成梁上游的全部脊椎骨端节点 。节点编号从401开始
*do,i,1,51,1
k,i+400,bnup(i,1),bnup(i,2),bnup(i,3)
*enddo
!生成梁下游的全部脊椎骨端节点 。节点编号从501开始
*do,i,1,51,1
k,i+500,bndown(i,1),bndown(i,2),bndown(i,3)
*enddo

! 生成梁上游的全部索端节点 。节点编号从601开始 ，

k,1+600,cabup(1,1),cabup(1,2)-spring,cabup(1,3)

*do,i,2,4,1
k,i+600,cabup(i,1),cabup(i,2),cabup(i,3)
*enddo

k,5+600,cabup(5,1),cabup(5,2)-spring,cabup(5,3)

*do,i,6,46,1
k,i+600,cabup(i,1),cabup(i,2),cabup(i,3)
*enddo

k,47+600,cabup(47,1),cabup(47,2)-spring,cabup(47,3)

*do,i,48,50,1
k,i+600,cabup(i,1),cabup(i,2),cabup(i,3)
*enddo

k,51+600,cabup(51,1),cabup(51,2)-spring,cabup(51,3)

!生成梁下游的全部索端节点 。节点编号从701开始 ，

k,1+700,cabdown(1,1),cabdown(1,2)-spring,cabdown(1,3)

*do,i,2,4,1
k,i+700,cabdown(i,1),cabdown(i,2),cabdown(i,3)
*enddo

k,5+700,cabdown(5,1),cabdown(5,2)-spring,cabdown(5,3)

*do,i,6,46,1
k,i+700,cabdown(i,1),cabdown(i,2),cabdown(i,3)
*enddo

k,47+700,cabdown(47,1),cabdown(47,2)-spring,cabdown(47,3)

*do,i,48,50,1
k,i+700,cabdown(i,1),cabdown(i,2),cabdown(i,3)
*enddo

k,51+700,cabdown(51,1),cabdown(51,2)-spring,cabdown(51,3)

!生成梁的轴线 。编号从300开始编号。 编号为从300到349号。

numstr,line,300
*do,i,1,50,1
l,i+300,i+301
*enddo

!由前面描述的设计数据，301，305，347，351 分别是桥主梁轴线在桥墩位置处的坐标。
numstr,line,361
l,361,301
l,362,305
l,363,347
l,364,351

!生成桥纵轴向水平弹性约束，编号为370号。
numstr,line,370
l,301,371
l,314,372
l,338,373

!生成桥竖向水平弹性约束，编号为380号。
numstr,line,380
l,301,355
l,305,356
l,347,357
l,351,358

!生成全部脊椎骨上的鱼刺梁。 编号从400开始编号。
numstr,line,400
*do,i,1,51,1
l,i+300,i+400
l,i+300,i+500
*enddo

!生成全部上游索端刚臂梁。 编号从600开始编号。
numstr,line,600
*do,i,1,51,1
l,i+400,i+600
*enddo

!生成全部下游索端刚臂梁。 编号从700开始编号。
numstr,line,700
*do,i,1,51,1
l,i+500,i+700
*enddo

!以上已经全部建立完毕梁单元和塔单元。下面将建立索单元。

!主梁北塔北侧上游索端点编号我们从最短索到最长索存到 BCABLE1 。
!主梁北塔南侧上游索端点编号我们从最短索到最长索存到 BCABLE2 。
!主梁南塔北侧上游索端点编号我们从最短索到最长索存到 BCABLE3 。
!主梁南塔南侧上游索端点编号我们从最短索到最长索存到 BCABLE4 。

!主梁北塔北侧下游索端点编号我们从最短索到最长索存到 BCABLE5 。
!主梁北塔南侧下游索端点编号我们从最短索到最长索存到 BCABLE6 。
!主梁南塔北侧下游索端点编号我们从最短索到最长索存到 BCABLE7 。
!主梁南塔南侧下游索端点编号我们从最短索到最长索存到 BCABLE8 。

!下面开始定义以上12个参数的维数。
*dim,BCABLE1,array,11,1,1
*dim,BCABLE2,array,11,1,1
*dim,BCABLE3,array,11,1,1
*dim,BCABLE4,array,11,1,1
*dim,BCABLE5,array,11,1,1
*dim,BCABLE6,array,11,1,1
*dim,BCABLE7,array,11,1,1
*dim,BCABLE8,array,11,1,1

!边界条件的说明：
!北塔北侧 1 号桥墩的节点号码为1（301、401、501、601、701）
!北塔北侧 2 号桥墩的节点号码为5（305、405、505、605、705）
!北塔 3 号桥墩（主塔自身）的节点号码为14（314、414、514、614、714）

!中心对称点的节点号码为 26 。
!北塔 4 号桥墩（主塔自身）的节点号码为38 。
!南塔南侧 5 号桥墩的节点号码为47 。
!南塔南侧 5 号桥墩的节点号码为51 。

BCABLE1(1)=613,612,611,610,609,608,607,606
BCABLE1(9)=604,603,602

BCABLE2(1)=615,616,617,618,619,620,621,622
BCABLE2(9)=623,624,625

BCABLE3(1)=637,636,635,634,633,632,631,630
BCABLE3(9)=629,628,627

BCABLE4(1)=639,640,641,642,643,644,645,646
BCABLE4(9)=648,649,650
*do,i,1,11,1
BCABLE5(i,1)=BCABLE1(i,1)+100
BCABLE6(i,1)=BCABLE2(i,1)+100
BCABLE7(i,1)=BCABLE3(i,1)+100
BCABLE8(i,1)=BCABLE4(i,1)+100
*enddo

!下面再来使用查询函数KP来得到桥塔上索节点的编号，并把他们保存到TCABLE 数组总去。
!北塔上游塔上索节点从最短索到最长索的编号我们将收集来存在数组 TCABLE1 。
!北塔下游塔上索节点从最短索到最长索的编号我们将收集来存在数组 TCABLE2 。
!南塔上游塔上索节点从最短索到最长索的编号我们将收集来存在数组 TCABLE3 。
!南塔下游塔上索节点从最短索到最长索的编号我们将收集来存在数组 TCABLE4 。

*dim,tcable1,array,11,1,1
*dim,tcable2,array,11,1,1
*dim,tcable3,array,11,1,1
*dim,tcable4,array,11,1,1
*dim,cabcoord,array,11,1,1

! 北桥塔下游索中心点的坐标为（6250，，126200）
! 北桥塔上游索中心点的坐标为（-6250，，126200）
! 南桥塔索下游中心点的坐标为（6250，，386200）
! 南桥塔索上游中心点的坐标为（-6250，，386200）
! 南桥塔索中心点的坐标见下列数组 cabcoord 中。

cabcoord(1)=59000,57800,56600,55400,54200,53000,51750,50450
cabcoord(9)=48850,46000,46000

*do,i,1,11,1
num=12-i
tcable1(num)=kp(-62500,cabcoord(i),126200)
tcable2(num)=kp(62500,cabcoord(i),126200)
tcable3(num)=kp(-62500,cabcoord(i),386200)
tcable4(num)=kp(62500,cabcoord(i),386200)
*enddo

!生成全部塔上连接1、2两号索的小刚臂梁。 编号从801开始。
numstr,line,801
l,50,tcable1(1)
l,51, tcable1(1)
l,52,tcable2(1)
l,53,tcable2(1)
l,54,tcable3(1)
l,55,tcable3(1)
l,56,tcable4(1)
l,57,tcable4(1)

!下面用上面 TCABLE 、BCABLE 里面的数据来建立索单元。索单元的编号规则如下：
!从1001号开始，每排索面为一组，
!第一组：北塔北上游 。号码为1001-1011
!第二组：北塔南上游 。号码为1101-1111
!第三组：南塔北上游 。号码为1201-1211
!第四组：南塔南上游 。号码为1301-1311
!第五组：北塔北下游 。号码为1401-1411
!第六组：北塔南下游 。号码为1501-1511
!第七组：南塔北下游 。号码为1601-1611
!第八组：南塔南下游 。号码为1701-1711

!建立第一组索单元。北塔北上游
numstr,line,1001
l,50,bcable1(1)
l,50,bcable1(2)
*do,i,3,11,1
l,tcable1(i),bcable1(i)
*enddo
!建立第二组索单元 。北塔南上游 。
numstr,line,1101
l,51,bcable2(1)
l,51,bcable2(2)
*do,i,3,11,1
l,tcable1(i),bcable2(i)
*enddo
!建立第三组索单元 。南塔北上游。
numstr,line,1201
l,54,bcable3(1)
l,54,bcable3(2)
*do,i,3,11,1
l,tcable3(i),bcable3(i)
*enddo
!建立第四组索单元。 南塔南上游。
numstr,line,1301
l,55,bcable4(1)
l,55,bcable4(2)
*do,i,3,11,1
l,tcable3(i),bcable4(i)
*enddo
!建立第五组索单元 。北塔北下游 。
numstr,line,1401
l,52,bcable5(1)
l,52,bcable5(2)
*do,i,3,11,1
l,tcable2(i),bcable5(i)
*enddo
!建立第六组索单元 。北塔南下游 。
numstr,line,1501
l,53,bcable6(1)
l,53,bcable6(2)
*do,i,3,11,1
l,tcable2(i),bcable6(i)
*enddo
!建立第七组索单元 。南塔北下游 。
numstr,line,1601
l,56,bcable7(1)
l,56,bcable7(2)
*do,i,3,11,1
l,tcable4(i),bcable7(i)
*enddo
!建立第八组索单元 。南塔南下游 。
numstr,line,1701
l,57,bcable8(1)
l,57,bcable8(2)
*do,i,3,11,1
l,tcable4(i),bcable8(i)
*enddo
SAVE
!整个前处理单元建立完毕

!************************************************************************************************

!永和桥网格划分

!************************************************************************************************
!下面将塔的截面全部读进内存。
sectype,1,beam,mesh !指定读入的截面类型在后面使用中编号
secoffset,shrc !指定截面在梁纵轴上的偏移量
secread,'yhtsec1','sect','e:\ansyswork\yhbridge\seclib',mesh !读入截面。如果截面保存在其他路径，可以采用绝对路径的方法确定

sectype,2,beam,mesh !指定读入的截面类型在后面使用中编号
secoffset,shrc!指定截面在梁纵轴上的偏移量
secread,'yhtsec2','sect','e:\ansyswork\yhbridge\seclib',mesh !读入截面。如果截面保存在其他路径，可以采用绝对路径的方法确定

sectype,3,beam,mesh !指定读入的截面类型在后面使用中编号
secoffset,shrc !指定截面在梁纵轴上的偏移量
secread,'yhtsec3','sect','e:\ansyswork\yhbridge\seclib',mesh !读入截面。如果截面保存在其他路径，可以采用绝对路径的方法确定

sectype,4,beam,mesh !指定读入的截面类型在后面使用中编号
secoffset,shrc !指定截面在梁纵轴上的偏移量
secread,'yhtsec4','sect','e:\ansyswork\yhbridge\seclib',mesh !读入截面。如果截面保存在其他路径，可以采用绝对路径的方法确定

sectype,5,beam,mesh !指定读入的截面类型在后面使用中编号
secoffset,cent !指定截面在梁纵轴上的偏移量
secread,'yhbsec1','sect','e:\ansyswork\yhbridge\seclib',mesh !读入截面。如果截面保存在其他路径，可以采用绝对路径的方法确定

!对于188\189单元，不需要实常数，但是如果不给其指定一个空的实常数号，会自动把第一个实常数号给它分配。产生警告信息。

!***********索的分布编号说明*************
!北塔岸侧索：1001-1011；1401-1411
!北塔河侧索：1101-1111；1501-1511
!南塔河侧索：1201-1211；1601-1411
!南塔岸侧索：1301-1311；1701-1711
!按照从上游到下游，从北到南的顺序！

numstr,elem,1
*do,i,1,11
lsel,s,line,,1000+i
latt,i+10,i+10,4
esize,,1
lmesh,all
*enddo

numstr,elem,101
*do,i,1,11
lsel,s,line,,1100+i
latt,i+30,i+30,4
esize,,1
lmesh,all
*enddo

numstr,elem,201
*do,i,1,11
lsel,s,line,,1200+i
latt,i+50,i+50,4
esize,,1
lmesh,all
*enddo

numstr,elem,301
*do,i,1,11
lsel,s,line,,1300+i
latt,i+70,i+70,4
esize,,1
lmesh,all
*enddo

numstr,elem,401
*do,i,1,11
lsel,s,line,,1400+i
latt,i+90,i+90,4
esize,,1
lmesh,all
*enddo

numstr,elem,501
*do,i,1,11
lsel,s,line,,1500+i
latt,i+110,i+110,4
esize,,1
lmesh,all
*enddo

numstr,elem,601
*do,i,1,11
lsel,s,line,,1600+i
latt,i+130,i+130,4
esize,,1
lmesh,all
*enddo

numstr,elem,701
*do,i,1,11
lsel,s,line,,1700+i
latt,i+150,i+150,4
esize,,1
lmesh,all
*enddo

! 再给四个桥墩位置的横桥向水平弹性约束划分单元，同样采用LINK 8 单元。
! 采用编号控制，其LINE 的编号为361-364

numstr,elem,801

lsel,s,line,,361,361,1,0
latt,5,4,5
esize,,1
lmesh,all

lsel,s,line,,362,362,1,0
latt,5,5,5
esize,,1
lmesh,all
lsel,s,line,,363,363,1,0
latt,5,5,5
esize,,1
lmesh,all
lsel,s,line,,364,364,1,0
latt,5,4,5
esize,,1
lmesh,all

! 再给桥墩位置的桥纵轴向水平弹性约束划分单元，同样采用LINK 8 单元。
! 采用编号控制，其LINE 的编号为370

lsel,s,line,,370,370,1,0
latt,5,6,5
esize,,1
lmesh,all

! 再给住塔位置的桥垂直轴向水平弹性约束划分单元，同样采用LINK 8 单元。
! 采用编号控制，其LINE 的编号为371,372

lsel,s,line,,371,372,1,0
latt,5,3,5
esize,,1
lmesh,all

numstr,elem,1000
!开始给鱼脊椎骨主梁划分单元
lsel,s,line,,300,349,1,0
latt,1,1,1,,300,,5
esize,2000
lmesh,all

numstr,elem,2000
!开始给鱼脊刚性横梁主梁划分单元 .
lsel,s,line,,400,501,1,0
latt,2,2,2,,,,
esize,,1
lmesh,all

!开始给鱼脊椎骨下的刚性刚臂梁划分单元 .
lsel,s,line,,601,649,1,0
lsel,a,line,,701,749,1,0
lsel,u,line,,704,704,1,0
lsel,u,line,,604,604,1,0
lsel,u,line,,646,646,1,0
lsel,u,line,,746,746,1,0
latt,2,2,2,,,,
esize,,1
lmesh,all

!开始给边界条件的扭转弹簧划分单元中间两个辅助墩.
lsel,s,line,,704,704,1,0
lsel,a,line,,604,604,1,0
lsel,a,line,,646,646,1,0
lsel,a,line,,746,746,1,0
latt,5,7,5,,,,
esize,,1
lmesh,all

!开始给边界条件的扭转弹簧划分单元。 两个桥头桥墩.
lsel,s,line,,700,700,1,0
lsel,a,line,,600,600,1,0
lsel,a,line,,650,650,1,0
lsel,a,line,,750,750,1,0
latt,5,7,5,,,,
esize,,1
lmesh,all

!开始给边界条件的竖向弹簧划分单元。
lsel,s,line,,380,380,1,0
lsel,a,line,,381,381,1,0
lsel,a,line,,382,382,1,0
lsel,a,line,,383,383,1,0
latt,5,8,5,,,,
esize,,1
lmesh,all

!开始给塔上的小刚臂梁划分单元 。
lsel,s,line,,801,808,1,0
latt,2,2,2,,,,
esize,,1
lmesh,all

numstr,elem,4000

!先给北塔分配属性和单元划分，区别仅仅是关键点及线的编号增量为100。
lsel,s,line,,101,101,1,0
latt,3,25,3,,131,,1
esize,2000
lmesh,all

lsel,s,line,,102,102,1,0
latt,3,25,3,,132,,1
esize,2000
lmesh,all

lsel,s,line,,103,125,2,0
latt,3,25,3,,131,,3
esize,2000
lmesh,all

lsel,s,line,,104,126,2,0
latt,3,25,3,,132,,3
esize,2000
lmesh,all

lsel,s,line,,127,127,1,0
latt,3,25,3,,133,,4
esize,2000
lmesh,all

lsel,s,line,,128,128,1,0
latt,3,25,3,,134,,4
esize,2000
lmesh,all

numstr,elem,5000
!再给南塔分配属性及单元划分
lsel,s,line,,1,1,1,0
latt,3,25,3,,31,,1
esize,2000
lmesh,all

lsel,s,line,,2,2,1,0
latt,3,25,3,,32,,1
esize,2000
lmesh,all

lsel,s,line,,3,25,2,0
latt,3,25,3,,31,,3
esize,2000
lmesh,all

lsel,s,line,,4,26,2,0
latt,3,25,3,,32,,3
esize,2000
lmesh,all

lsel,s,line,,27,27,1,0
latt,3,25,3,,33,,4
esize,2000
lmesh,all

lsel,s,line,,28,28,1,0
latt,3,25,3,,34,,4
esize,2000
lmesh,all
numstr,elem,4000
SAVE

!ANSYS 在 分网的时候，最后一次操作生成单元之后，相应单元的原来的几何模型（线、面）并没有被删除。还可以显示。
allsel
eplo

!*****************************************************************

!永和桥的边界条件

!*****************************************************************
!永和桥的模型边界条件、成桥荷载、及成桥状态静力分析。寻找边界分别是有 AXIS里面存储的梁轴线坐标去寻找的。
! 由于边界条件对振型和频率影响十分显著，所以四个桥墩全部采用弹性边界条件。然后根据动力参数进行修正。
/solu
csys,0

!给北塔北第一个桥墩(桥头)设置边界条件（横向水平约束）
nsel,s,loc,z,axis(1,3)-10,axis(1,3)+10
nsel,r,loc,y,axis(1,2)-10,axis(1,2)+10
nsel,r,loc,x,axis(1,1)+spring-10,axis(1,1)+spring+10
d,all,uy
d,all,ux
d,all,uz

!给北塔北第一个桥墩(桥头)设置边界条件（沿桥纵轴向水平约束）
nsel,s,loc,z,axis(1,3)-spring-10,axis(1,3)-spring+10
nsel,r,loc,y,axis(1,2)-10,axis(1,2)+10
nsel,r,loc,x,axis(1,1)-10,axis(1,1)+10
d,all,uy
d,all,ux
d,all,uz
!给北塔北第一个桥墩(桥头)设置边界条件（沿桥两侧竖直方向抗扭转弹簧边界条件）
nsel,s,loc,z,cabup(1,3)-10,cabup(1,3)+10
nsel,r,loc,y,cabup(1,2)-spring-10,cabup(1,2)-spring+10
nsel,r,loc,x,cabup(1,1)-10,cabup(1,1)+10
d,all,uy
d,all,ux
d,all,uz

nsel,s,loc,z,cabdown(1,3)-10,cabdown(1,3)+10
nsel,r,loc,y,cabdown(1,2)-spring-10,cabdown(1,2)-spring+10
nsel,r,loc,x,cabdown(1,1)-10,cabdown(1,1)+10
d,all,uy
d,all,ux
d,all,uz

!给北塔北第二个桥墩(辅助墩)设置边界条件（横向水平约束）
nsel,s,loc,z,axis(5,3)-10,axis(5,3)+10
nsel,r,loc,y,axis(5,2)-10,axis(5,2)+10
nsel,r,loc,x,axis(5,1)+spring-10,axis(5,1)+spring+10
d,all,uy
d,all,ux
d,all,uz

!给北塔北第一个桥墩(辅助墩)设置边界条件（沿桥两侧竖直方向抗扭转弹簧边界条件）
nsel,s,loc,z,cabup(5,3)-10,cabup(5,3)+10
nsel,r,loc,y,cabup(5,2)-spring-10,cabup(5,2)-spring+10
nsel,r,loc,x,cabup(5,1)-10,cabup(5,1)+10
d,all,uy
d,all,ux
d,all,uz

nsel,s,loc,z,cabdown(5,3)-10,cabdown(5,3)+10
nsel,r,loc,y,cabdown(5,2)-spring-10,cabdown(5,2)-spring+10
nsel,r,loc,x,cabdown(5,1)-10,cabdown(5,1)+10
d,all,uy
d,all,ux
d,all,uz

!给北塔设置边界条件
nsel,s,loc,z,axis(14,3)-10,axis(14,3)+10
nsel,r,loc,y,5000-10,5000+10
d,all,all

!给南塔设置边界条件
nsel,s,loc,z,axis(38,3)-10,axis(38,3)+10
nsel,r,loc,y,5000-10,5000+10
d,all,all

!给南塔北第二个桥墩(辅助墩)设置边界条件（横向水平约束）
nsel,s,loc,z,axis(47,3)-10,axis(47,3)+10
nsel,r,loc,y,axis(47,2)-10,axis(47,2)+10
nsel,r,loc,x,axis(47,1)+spring-10,axis(47,1)+spring+10
d,all,uy
d,all,ux
d,all,uz

!给南塔南第二个桥墩(辅助墩)设置边界条件（沿桥两侧竖直方向抗扭转弹簧边界条件）
nsel,s,loc,z,cabup(47,3)-10,cabup(47,3)+10
nsel,r,loc,y,cabup(47,2)-spring-10,cabup(47,2)-spring+10
nsel,r,loc,x,cabup(47,1)-10,cabup(47,1)+10
d,all,uy
d,all,ux
d,all,uz

nsel,s,loc,z,cabdown(47,3)-10,cabdown(47,3)+10
nsel,r,loc,y,cabdown(47,2)-spring-10,cabdown(47,2)-spring+10
nsel,r,loc,x,cabdown(47,1)-10,cabdown(47,1)+10
d,all,uy
d,all,ux
d,all,uz

!给南塔第一个桥墩(桥头)设置边界条件（横向水平约束）
nsel,s,loc,z,axis(51,3)-10,axis(51,3)+10
nsel,r,loc,y,axis(51,2)-10,axis(51,2)+10
nsel,r,loc,x,axis(51,1)+spring-10,axis(51,1)+spring+10
d,all,uy
d,all,ux
d,all,uz

!给南塔第一个桥墩(桥头)设置边界条件（沿桥两侧竖直方向抗扭转弹簧边界条件）
nsel,s,loc,z,cabup(51,3)-10,cabup(51,3)+10
nsel,r,loc,y,cabup(51,2)-spring-10,cabup(51,2)-spring+10
nsel,r,loc,x,cabup(51,1)-10,cabup(51,1)+10
d,all,uy
d,all,ux
d,all,uz

nsel,s,loc,z,cabdown(51,3)-10,cabdown(51,3)+10
nsel,r,loc,y,cabdown(51,2)-spring-10,cabdown(51,2)-spring+10
nsel,r,loc,x,cabdown(51,1)-10,cabdown(51,1)+10
d,all,uy
d,all,ux
d,all,uz

!给四个桥墩设置竖向边界条件（模拟盆式支座）
nsel,s,loc,z,axis(1,3)-10,axis(1,3)+10
nsel,r,loc,y,axis(1,2)-spring-10,axis(1,2)-spring+10
nsel,r,loc,x,axis(1,1)-10,axis(1,1)+10
d,all,uy
d,all,ux
d,all,uz

nsel,s,loc,z,axis(5,3)-10,axis(5,3)+10
nsel,r,loc,y,axis(5,2)-spring-10,axis(5,2)-spring+10
nsel,r,loc,x,axis(5,1)-10,axis(5,1)+10
d,all,uy
d,all,ux
d,all,uz

nsel,s,loc,z,axis(47,3)-10,axis(47,3)+10
nsel,r,loc,y,axis(47,2)-spring-10,axis(47,2)-spring+10
nsel,r,loc,x,axis(47,1)-10,axis(47,1)+10
d,all,uy
d,all,ux
d,all,uz

nsel,s,loc,z,axis(51,3)-10,axis(51,3)+10
nsel,r,loc,y,axis(51,2)-spring-10,axis(51,2)-spring+10
nsel,r,loc,x,axis(51,1)-10,axis(51,1)+10
d,all,uy
d,all,ux
d,all,uz

!给两个主桥塔处设置垂直桥纵轴线方向水平边界条件（模拟挡块约束）
nsel,s,loc,z,axis(14,3)-10,axis(38,3)+10
nsel,r,loc,y,axis(14,2)-10,axis(38,2)+10
nsel,r,loc,x,axis(14,1)+spring-10,axis(38,1)+spring+10
d,all,uy
d,all,ux
d,all,uz

SAVE
allsel
finish

!************************************************************************************************

!动力分析之模态分析

!************************************************************************************************
allsel

/SOLU
ANTYPE,MODAL
MODOPT,LANB,20
MXPAND,20,,,0
LUMPM,0
PSTRES,on
SOLVE
save
FINISH

/POST1
*GET,OMG1,MODE,1,FREQ
*GET,OMG2,MODE,2,FREQ
*GET,OMG3,MODE,3,FREQ
*GET,OMG4,MODE,4,FREQ
*GET,OMG5,MODE,5,FREQ
*GET,OMG6,MODE,6,FREQ
*GET,OMG7,MODE,7,FREQ
*GET,OMG8,MODE,8,FREQ
*GET,OMG9,MODE,9,FREQ
*GET,OMG10,MODE,10,FREQ

*GET,OMG11,MODE,11,FREQ
*GET,OMG12,MODE,12,FREQ
*GET,OMG13,MODE,13,FREQ
*GET,OMG14,MODE,14,FREQ
*GET,OMG15,MODE,15,FREQ
*GET,OMG16,MODE,16,FREQ
*GET,OMG17,MODE,17,FREQ
*GET,OMG18,MODE,18,FREQ
*GET,OMG19,MODE,19,FREQ
*GET,OMG20,MODE,20,FREQ

*dim,frequence,array,1,20,1

frequence(1,1,1)=omg1
frequence(1,2,1)=omg2
frequence(1,3,1)=omg3
frequence(1,4,1)=omg4
frequence(1,5,1)=omg5
frequence(1,6,1)=omg6
frequence(1,7,1)=omg7
frequence(1,8,1)=omg8
frequence(1,9,1)=omg9
frequence(1,10,1)=omg10
frequence(1,11,1)=omg11
frequence(1,12,1)=omg12
frequence(1,13,1)=omg13
frequence(1,14,1)=omg14
frequence(1,15,1)=omg15
frequence(1,16,1)=omg16
frequence(1,17,1)=omg17
frequence(1,18,1)=omg18
frequence(1,19,1)=omg19
frequence(1,20,1)=omg20

!*cfopen,iteration,dat,'e:\ansyswork\yhdata\output',APPEND,
!*vwrite,NSKIP,NSLOOP(1:1),NSLOOP(1,2),NSLOOP(1,3),NSLOOP(1,4),NSLOOP(2,1),NSLOOP(2,2),NSLOOP(2,3),NSLOOP(2,4),NSLOOP(3,1),NSLOOP(3,2),NSLOOP(3,3),NSLOOP(3,4)
!(F3.0,' ',f8.0,' ',f8.0,' ',f8.0,' ',f8.0,' ',f8.0,' ',f8.0,' ',f8.0,' ',f8.0,' ',f8.0,' ',f8.0,' ',f8.0,' ',f8.0)
!*CFCLOS

*cfopen,optparament,dat,'e:\ansyswork\yhdata\output',APPEND,
*vwrite,NSKIP,optpara(1,1),optpara(1,2),optpara(1,3),optpara(1,4),optpara(1,5),optpara(1,6),optpara(1,7),optpara(1,8),optpara(1,9),optpara(1,10),optpara(1,11),optpara(1,12)
(F3.0,' ',f10.4,' ',f10.4,' ',f10.4,' ',f10.4,' ',f10.4,' ',f10.4,' ',f10.4,' ',f10.4,' ',f10.4,' ',f10.4,' ',f10.4,' ',f10.4)
*CFCLOS

*cfopen,frequence,dat,'e:\ansyswork\yhdata\output',APPEND,
*vwrite,NSKIP,frequence(1,1,1),frequence(1,2,1),frequence(1,3,1),frequence(1,4,1),frequence(1,5,1),frequence(1,6,1),frequence(1,7,1),frequence(1,8,1),frequence(1,9,1),frequence(1,10,1),frequence(1,11,1),frequence(1,12,1),frequence(1,13,1),frequence(1,14,1),frequence(1,15,1),frequence(1,16,1),frequence(1,17,1),frequence(1,18,1)
(F3.0,' ',f10.8,' ',f10.8,' ',f10.8,' ',f10.8,' ',f10.8,' ',f10.8,' ',f10.8,' ',f10.8,' ',f10.8,' ',f10.8,' ',f10.8,' ',f10.8,' ',f10.8,' ',f10.8,' ',f10.8,' ',f10.8,' ',f10.8,' ',f10.8)
*CFCLOS
*SET,object,(OMG1/0.176-1)**2+(OMG2/0.427-1)**2+(OMG3/0.618-1)**2+(OMG6/1.105-1)**2+(OMG11/1.347-1)**2
object=object/5
SAVE

!**************************************************************************************************************************************8

! 提取出模态位移，保存下来，以供MATLAB 调用

!****************************************************************************************************************************************

*dim,Beamnode,array,51,1
*dim,towernode,array,4,1
allsel
*do,i,1,51
Beamnode(i,1)=node(axis(i,1),axis(i,2),axis(i,3))
*enddo
towernode(1,1)=node(6250,65000,axis(14,3))
towernode(2,1)=node(-6250,65000,axis(14,3))
towernode(3,1)=node(6250,65000,axis(38,3))
towernode(4,1)=node(-6250,65000,axis(38,3))

! modenum 是选择存储的模态数
! xmode 是X方向的模态位移
! ymode 是Y方向的模态位移
modenum=20
*dim,xmode,array,51,modenum
*dim,ymode,array,51,modenum
*dim,zmode,array,51,modenum
*dim,summode,array,51,modenum
*dim,rotzmode,array,51,modenum
*dim,towermodez,array,4,modenum
*dim,towermodex,array,4,modenum

*do,j,1,modenum
set,1,j
*do,i,1,51
*get,xmode(i,j),node,Beamnode(i,1),u,x
*get,ymode(i,j),node,Beamnode(i,1),u,y
*get,zmode(i,j),node,Beamnode(i,1),u,z
*get,summode(i,j),node,Beamnode(i,1),u,sum
*get,rotzmode(i,j),node,Beamnode(i,1),rot,z
*enddo
*enddo

*do,i,1,51
*do,j,1,modenum
rotzmode(i,j)=rotzmode(i,j)*12000
*enddo
*enddo


*do,j,1,modenum
set,1,j
*do,i,1,4
*get,towermodez(i,j),node,towernode(i,1),u,z
*get,towermodex(i,j),node,towernode(i,1),u,x
*enddo
*enddo
FILENAME='xmode'
*cfopen,%FILENAME%%NSKIP%,dat,'e:\ansyswork\yhdata\output'
*vwrite,sequ,xmode(1,1),xmode(1,2),xmode(1,3),xmode(1,4),xmode(1,5),xmode(1,6),xmode(1,7),xmode(1,8),xmode(1,9),xmode(1,10),xmode(1,11),xmode(1,12),xmode(1,13),xmode(1,14),xmode(1,15),xmode(1,16),xmode(1,17),xmode(1,18)
(F3.0,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12)
*CFCLOS
FILENAME='ymode'
*cfopen,%FILENAME%%NSKIP%,dat,'e:\ansyswork\yhdata\output'
*vwrite,sequ,ymode(1,1),ymode(1,2),ymode(1,3),ymode(1,4),ymode(1,5),ymode(1,6),ymode(1,7),ymode(1,8),ymode(1,9),ymode(1,10),ymode(1,11),ymode(1,12),ymode(1,13),ymode(1,14),ymode(1,15),ymode(1,16),ymode(1,17),ymode(1,18)
(F3.0,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12)
*CFCLOS

FILENAME='zmode'
*cfopen,%FILENAME%%NSKIP%,dat,'e:\ansyswork\yhdata\output'
*vwrite,sequ,zmode(1,1),zmode(1,2),zmode(1,3),zmode(1,4),zmode(1,5),zmode(1,6),zmode(1,7),zmode(1,8),zmode(1,9),zmode(1,10),zmode(1,11),zmode(1,12),zmode(1,13),zmode(1,14),zmode(1,15),zmode(1,16),zmode(1,17),zmode(1,18)
(F3.0,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12)
*CFCLOS
FILENAME='summode'
*cfopen,%FILENAME%%NSKIP%,dat,'e:\ansyswork\yhdata\output'
*vwrite,sequ,summode(1,1),summode(1,2),summode(1,3),summode(1,4),summode(1,5),summode(1,6),summode(1,7),summode(1,8),summode(1,9),summode(1,10),summode(1,11),summode(1,12),summode(1,13),summode(1,14),summode(1,15),summode(1,16),summode(1,17),summode(1,18)
(F3.0,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12)
*CFCLOS
FILENAME='rotzmode'
*cfopen,%FILENAME%%NSKIP%,dat,'e:\ansyswork\yhdata\output'
*vwrite,sequ,rotzmode(1,1),rotzmode(1,2),rotzmode(1,3),rotzmode(1,4),rotzmode(1,5),rotzmode(1,6),rotzmode(1,7),rotzmode(1,8),rotzmode(1,9),rotzmode(1,10),rotzmode(1,11),rotzmode(1,12),rotzmode(1,13),rotzmode(1,14),rotzmode(1,15),rotzmode(1,16),rotzmode(1,17),rotzmode(1,18)
(F3.0,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12)
*CFCLOS
FILENAME='towermodez'
*cfopen,%FILENAME%%NSKIP%,dat,'e:\ansyswork\yhdata\output'
*vwrite,sequ,towermodez(1,1),towermodez(1,2),towermodez(1,3),towermodez(1,4),towermodez(1,5),towermodez(1,6),towermodez(1,7),towermodez(1,8),towermodez(1,9),towermodez(1,10),towermodez(1,11),towermodez(1,12),towermodez(1,13),towermodez(1,14),towermodez(1,15),towermodez(1,16),towermodez(1,17),towermodez(1,18)
(F3.0,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12)
*CFCLOS
FILENAME='towermodex'
*cfopen,%FILENAME%%NSKIP%,dat,'e:\ansyswork\yhdata\output'
*vwrite,sequ,towermodex(1,1),towermodex(1,2),towermodex(1,3),towermodex(1,4),towermodex(1,5),towermodex(1,6),towermodex(1,7),towermodex(1,8),towermodex(1,9),towermodex(1,10),towermodex(1,11),towermodex(1,12),towermodex(1,13),towermodex(1,14),towermodex(1,15),towermodex(1,16),towermodex(1,17),towermodex(1,18)
(F3.0,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12,' ',f16.12)
*CFCLOS