平面8节点等参元子程序.doc

平面8节点等参元子程序 !平面8节点等参元程序 !清华大学土木工程系：陆新征 module Elem_Rect8 !八节点等参元 implicit none integer (kind(1)),parameter::ikind=(kind⑴) integer (kind(1))5parameter::rkind=(kind(O.dO)) type :: typ_Kcol real(rkind),pointer:: Row(:) end type typ_Kcol type :: typ_GValue !总体控制变量 integer(ikind) :: NNode, NEIem, NLoad, NMat, NSupport integer(ikind) :: NGIbDOF !整体 口 由度总数 integer(ikind) :: NGENS, NodeDOF.EIemNodeNo integer(ikind) :: Nlnt end type typ_GValue type Typ_Node !定义节点类型 real(rkind) :: coord(2) integer(ikind) :: GDOF(2) real(rkind) :: DISP(2) real(rkind) :: dDISP(2) !节点坐标 !整体H由度编码 !节点位移 !节点位移增量 real(rkind) :: dForce(2) !节点不平衡力 end type typ_Node Type typ_lntPoint !定义枳分点参数 real(rkind) :: EPS⑶ real(rkind) :: SIG(3) real(rkind) :: D(3,3) real(rkind) :: B(3,16) real(rkind) :: DET J end type Typ」ntPoint !应变 !应力 !本构矩阵 !几何矩阵 !雅克比行列式 type Typ_Rect8 !定义实体单元 integer(ikind) :: NodeNo(8) !节点编号  real(rkind) :: E real(rkind) :: u real(rkind) :: t real(rkind) :: EK(16,16) type(typ_intpoint) :: lntP(9) end type Typ_Rect8 type Typ_Load !定义荷载类型 integer(ikind) :: NodeNo real(rkind) :: Value(2) end type Typ_Load type Typ_Support!定义支座类型 integer(ikind) :: NodeNo integer(ikind) :: DOF real(rkind) :: Value end type Typ_Support !弹性模量!泊松比!单元厚度 !弹性模量 !泊松比 !单元厚度 !单元刚度矩阵 !积分点 !荷载作用节点编号 !荷载大小 !支座节点编号 !支座约束自由度 !支座位移大小 subroutine Get_Elem_K(GValue,Elem,Node) !核心程序，得到单元的刚度矩阵 implicit none type(typ_GValue) :: GValue type(typ_Node) :: Node(:) type(typ_Rec⑹::Elem(:) real*8 :: lnvJ(2,2) real*8:: Coord(8,2),lntPCoord(9,2),lntPwt(9),dN(2,8) integer:: l,ElemNo call Get_lntP_Prop(lntPCoord,lntPWt)!计算积分点信息 do ElemNo=1 ,size(Elem) Elem(ElemNo)%EK=0.0 do 1=1, GValue%ElemNodeNo; !得到节点处标 coord(l,:)=Node(Elem(ElemNo)%NodeNo(i))%coord; end do Elem(ElemNo)%EK=OdO DO 1=1 ,size(Elem(ElemNo)%lntP) !计算本构矩阵 call Get_D(Elem(ElemNo)%lntP(l)%D,Elem(ElemNo)%E,Elem(ElemNo)%u) !计算形函数对局部坐标导数 call Get_dN_dxi(dN,lntPCoord(i,1),lntPCoord(i,2)) lnvJ=matmul(dN, Coord) !得到雅克比行列式 Elem(ElemNo)%lntP(l)%DETJ=lnvJ(1,1)*lnvJ(2,2)-lnvJ(1,2)*lnvJ(2,1) ln