Fluidmach 发表于 2016-3-8 11:00

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

  虽然采用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
页: [1]
查看完整版本: ANSYS求结构整体刚度矩阵逆矩阵APDL命令流