爱华网虽然采用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.7hjfmyfile,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为最后求的刚度矩阵的逆矩阵！逆矩阵的求解算法采用的是高斯约旦消元法！
采用高斯约旦消元法，易于实现，但是求解的速度比较慢，对于大型的整体刚度矩阵的求逆，比较耗时，对其计算机的配置要求比较高！希望能后续读者能采用更好更简洁的算法，提高的计算的效率，节约计算的时间！
