上面等式右面的雅各比逆矩阵,我们已经求得,形函数对st坐标的导数也已经求得,因此,上面引入的矩阵就可以求出来,然后再按照元素对应的规律将上述矩阵以此存放到B矩阵中。
在计算质量矩阵的时候,首先计算单元的额协调质量矩阵,然后将行质量集中到对角元,让后将非对角元素置为零,具体可参考帖子:
单元的阻尼矩阵可按照比例阻尼求解,即刚度矩阵和质量矩阵的线性组合。
求解刚度矩阵的主程序为
subroutine KKmartix(KK,coords,DMatx,mcrd,nnode,jelem)
INCLUDE 'ABA_PARAM.INC'
double precision coords(mcrd,nnode),dmatx(3,3)
double precision INTEGERPOINT(3),wgt(3)
double precision KK(2*nnode,2*nnode)
double precision B(3,2*NNODE),J(2,2)
double precision DETJ,G,H,COFF
double precision NPXY(2,NNODE)
INTEGERPOINT(1)= 0.D0
INTEGERPOINT(2)= 0.7745966692D0
INTEGERPOINT(3)=-0.7745966692D0
wgt(1)=0.8888888889D0
wgt(2)=0.5555555556D0
wgt(3)=0.5555555556D0
B=0.d0
J=0.d0
DETJ=0.d0
G=0.d0
H=0.d0
COFF=0.d0
NPXY=0.d0
do I=1,3
do II=1,3
G=INTEGERPOINT(I)
H =INTEGERPOINT(II)
call CalBJ(B,J,DETJ,G,H,MCRD,NNODE,JELEM,COORDS)
COFF=DETJ*WGT(I)*WGT(II)
KK=KK+MATMUL(MATMUL(TRANSPOSE(B),DMatx),B)*COFF
enddo
enddo
return
end
求解质量矩阵的主程序为
subroutine MMmartix(MM,COORDS,PROPS,MCRD,NNODE,JELEM)
INCLUDE 'ABA_PARAM.INC'
INTEGER MCRD,NNODE,JELEM
DOUBLE PRECISION MM(2*NNODE,2*NNODE),COORDS(MCRD,NNODE)
DOUBLE PRECISION NN(2,16)
DOUBLE PRECISION INTEGERPOINT(3),WGT(3)
DOUBLE PRECISION G,H,DETJ,DENSITY,COFF
DOUBLE PRECISION J(2,2),PROPS(*)
MM=0.D0
wgt=0.D0
JJ=0.D0
NN=0.D0
G=0.D0
H=0.D0
j=0.D0
DETJ=0.D0
Coff=0.D0
DENSITY=0.D0
INTEGERPOINT(1)= 0.D0
INTEGERPOINT(2)= 0.7745966692D0
INTEGERPOINT(3)=-0.7745966692D0
wgt(1)=0.8888888889D0
wgt(2)=0.5555555556D0
wgt(3)=0.5555555556D0
DENSITY=PROPS(3)
DO I=1,3
DO II=1,3
G =INTEGERPOINT(I)
H =INTEGERPOINT(II)
CALL CalDETJ(J,DETJ,G,H,MCRD,NNODE,JELEM,COORDS)
CALL CalNN(NN,G,H)
COFF=DETJ*WGT(I)*WGT(II)*DENSITY
MM =MM+(MATMUL(TRANSPOSE(NN),NN))*COFF
ENDDO
ENDDO
return
end
至此,单元的刚度矩阵、质量矩阵和阻尼矩阵已经求得,uel程序的核心部分已经完成,下一个帖子会更新自编CPE8/CPS8 uel程序的计算结果。