ANSYS求结构整体刚度矩阵逆矩阵APDL命令流

Fluidmach

　　虽然采用matlab等相关软件对于大型矩阵的求解非常方便，但是这里涉及到ANSYS与matlab接口的处理问题，比较繁琐，本文编制一个直接在ANSYS中提取结构整体刚度矩阵和对整体刚度求逆的APDL程序，能在ANSYS中直接操作刚度矩阵，比较实用。

　　利用ANSYS求结构刚度矩阵逆矩阵，模型采用单层球面网壳!

　　相关的APDL命令流如下所示。

　　finish

　　/clear

　　/filname,modal

　　/prep7

　　et,1,beam188

　　mp,ex,1,2.1e11

　　mp,prxy,1,0.28

　　mp,dens,1,7850

　　*afun,deg

　　!用户界面设计，输入几何参数

　　multipro,'start',4

　　*cset,1,3,f,'Enter the f',6.8

　　*cset,4,6,span,'Enter the span',30

　　*cset,7,9,Kn,'Enter the radial number',24

　　*cset,10,12,Nx,'Enter the node circle number',3

　　*cset,61,62,'please enter the geometry parameters'

　　multipro,'end'

　　!计算节点坐标位置，并定义节点

　　csys,2

　　r=(f*f+span*span/4)/(2*f)

　　Dalfa=atn(span/2/sqrt(r*r-(span/2)*(span/2)))/Nx

　　n,1,r,0,90

　　*do,i,1,Nx

　　*do,j,1,Kn

　　x=r

　　y=(j-1)*360/Kn

　　z=90-i*Dalfa

　　n,1+Kn*(i-1)+j,x,y,z,

　　*enddo

　　*enddo

　　!定义单元连接

　　ro=0.219/2

　　ri=(0.219-0.010*2)/2

　　sectype,1,beam,ctube,,,

　　secdata,ri,ro,8

　　type,1

　　mat,1

　　secnum,1

　　!定义环向杆连接

　　*do,i,1,Nx

　　*do,j,1,Kn-1,1

　　e,1+Kn*(i-1)+j,1+Kn*(i-1)+j+1

　　*enddo

　　e,1+Kn*i,1+i*Kn-Kn+1

　　*enddo

　　!定义径向杆连接

　　*do,j,1,Kn,1

　　e,1,1+j,

　　*enddo

　　*do,i,1,Nx-1

　　*do,j,1,Kn

　　e,1+Kn*(i-1)+j,1+Kn*i+j

　　*enddo

　　*enddo

　　!定义斜杆连接

　　*do,i,1,Nx-1

　　*do,j,1,Kn-1,1

　　e,1+Kn*(i-1)+j,1+Kn*i+j+1,1

　　*enddo

　　e,1+i*Kn,1+i*Kn+1

　　*enddo

　　!定义边界约束节点和节点荷载

　　*do,i,1,1+Kn*Nx

　　*if,i,le,1+Kn*(Nx-1),then

　　f,i,fz,-1000

　　*else

　　d,i,all

　　*endif

　　*enddo

　　finish

　　/solu

　　solve

　　/aux2

　　file,modal,full

　　hbmat,modal,txt,,ascII,stiff,yes !转换hb.file文件hb.txt

　　finish

　　*DIM,CONTLINE,,5 !定义一维数组

　　*VREAD,CONTLINE(1),MODAL,TXT,,,5,,,1 !跳过第1行后读入5个数据

　　(5F14.0)

　　PTRCRD=CONTLINE(2) !保存列指针总行数

　　INDCRD=CONTLINE(3) !保存行索引总行数

　　VALCRD=CONTLINE(4) !保存矩阵元素总行数

　　RHSCRD=CONTLINE(5) !保存右边项总行数

　　*VREAD,CONTLINE(1),MODAL,TXT,,,4,,,2 !跳过第2行后读入4个数据

　　(A3,11X,4F14.0)

　　NROW=CONTLINE(2)$NCOL=CONTLINE(3) !保存刚度矩阵的行列数

　　STRLINE=$CONTLINE= !删除数组

　　*IF,RHSCRD,EQ,0,THEN !如果无右边项取LS0=4行，否则取LS0=5

　　LS0=4$*ELSE$LS0=5$*ENDIF

　　*DIM,POINTR,,PTRCRD !定义列指针数组

　　*DIM,ROWIND,,INDCRD !定义行索引数组

　　*DIM,VALUES,,VALCRD !定义矩阵元素值数组

　　*DIM,RHSVAL,,RHSCRD !定义右边项元素值数组

　　*VREAD,POINTR(1),MODAL,TXT,,,PTRCRD,,,LS0

　　(F14.0) !读入列指针数据

　　*VREAD,ROWIND(1),MODAL,TXT,,,INDCRD,,,LS0+PTRCRD

　　(F14.0) !读入行索引数据

　　*VREAD,VALUES(1),MODAL,TXT,,,VALCRD,,,LS0+PTRCRD+INDCRD

　　(D25.15) !读入矩阵元素数据

　　*VREAD,RHSVAL(1),MODAL,TXT,,,RHSCRD,,,LS0+PTRCRD+INDCRD+VALCRD

　　(D25.15) !读入右边项元素数据

　　*DIM,SMATR,,NROW,NCOL !定义矩阵行列数，满矩阵存储的矩阵

　　*DO,ICOL,1,NCOL !以列数循环

　　STACOL=POINTR(ICOL) !得到当前列指针(元素的列号)

　　ENDCOL=POINTR(ICOL+1) !得到下一列指针

　　*DO,IROW,STACOL,ENDCOL-1 !以当前列中的非零元素个数循环

　　TRUEROW=ROWIND(IROW) !得到当前元素的行号

　　SMATR(TRUEROW,ICOL)=VALUES(IROW) !按行列号将元素值保存到矩阵中

　　*ENDDO

　　*ENDDO !结束两个循环

　　*DO,IROW,1,NROW !形成上三角元素，进而得到满矩阵

　　*DO,ICOL,1,NCOL

　　SMATR(IROW,ICOL)=SMATR(ICOL,IROW)

　　*ENDDO

　　*ENDDO

　　POINTR=$ROWIND=$VALUES=$RHSVAL=$ICOL=$IROW=$LS0=$STACOL= !以下为删除临时变量和数组变量

　　ENDCOL=$TRUEROW=$TOTCRD=$PTRCRD=$INDCRD=$VALCRD=$RHSCRD=

　　!*cfopen,F:\Syntax-editor_V1.7\hjf\myfile,txt

　　!*CFCLOS

　　!对结构刚度矩阵求逆矩阵

　　FINISH

　　*DIM,T_SMATR,ARRAY,NROW,NCOL !定义中间存贮矩阵

　　*DIM,R_SMATR,ARRAY,NROW,NCOL !最终所求的刚度矩阵的逆矩阵

　　*DIM,IS,ARRAY,NROW

　　*DIM,JS,ARRAY,NCOL

　　*do,i,1,NROW

　　*do,j,1,NCOL

　　T_SMATR(i,j)=SMATR(i,j)

　　*enddo

　　*enddo

　　*do,k,1,NROW

　　FLAG=0.00

　　*do,i,k,NROW

　　*do,j,k,NCOL

　　*if,abs(SMATR(i,j)),gt,FLAG,then

　　FLAG=abs(SMATR(i,j))

　　IS(k)=i

　　JS(k)=j

　　*else

　　*endif

　　*enddo

　　*enddo

　　!*if,dd+1.0,eq,1.0,then

　　! l=0

　　!*else

　　!*endif

　　*do,j,1,NCOL

　　t=SMATR(k,j)

　　SMATR(k,j)=SMATR(IS(k),j)

　　SMATR(IS(k),j)=t

　　*enddo

　　*do,i,1,NROW

　　t=SMATR(i,k)

　　SMATR(i,k)=SMATR(i,JS(k))

　　SMATR(i,JS(k))=t

　　*enddo

　　SMATR(k,k)=1/SMATR(k,k)

　　*do,j,1,NCOL

　　*if,j,ne,k,then

　　SMATR(k,j)=SMATR(k,j)*SMATR(k,k)

　　*else

　　*endif

　　*enddo

　　*do,i,1,NROW

　　*if,i,ne,k,then

　　*do,j,1,NCOL

　　*if,j,ne,k,then

　　SMATR(i,j)=SMATR(i,j)-SMATR(i,k)*SMATR(k,j)

　　*else

　　*endif

　　*enddo

　　*else

　　*endif

　　*enddo

　　*do,i,1,NROW

　　*if,i,ne,k,then

　　SMATR(i,k)=-SMATR(i,k)*SMATR(k,k)

　　*else

　　*endif

　　*enddo

　　*enddo

　　*do,k,NROW,1,-1

　　*do,j,1,NCOL

　　t=SMATR(k,j)

　　SMATR(k,j)=SMATR(JS(k),j)

　　SMATR(JS(k),j)=t

　　*enddo

　　*do,i,1,NROW

　　t=SMATR(i,k)

　　SMATR(i,k)=SMATR(i,IS(k))

　　SMATR(i,IS(k))=t

　　*enddo

　　*enddo

　　*do,i,1,NROW

　　*do,j,1,NCOL

　　R_SMATR(i,j)=0

　　*do,l,1,NROW

　　R_SMATR(i,j)=R_SMATR(i,j)+SMATR(i,l)*T_SMATR(l,j)

　　*enddo

　　*enddo

　　*enddo

　　R_SMATR为最后求的刚度矩阵的逆矩阵!逆矩阵的求解算法采用的是高斯约旦消元法!

　　采用高斯约旦消元法，易于实现，但是求解的速度比较慢，对于大型的整体刚度矩阵的求逆，比较耗时，对其计算机的配置要求比较高!希望能后续读者能采用更好更简洁的算法，提高的计算的效率，节约计算的时间!



转自：http://blog.sina.com.cn/s/blog_5e94ccea01010ck5.html