利用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://s1/middle/4af3952e49659e1bc31c0&690