平面8节点等参元完整源程序.doc

可执行程序 Planar_8_Nodes.f90!=! 平面8节点等参元完整程序! 清华大学土木工程系：陆新征！程序说明 !=module Elem_Rect8 ! 八节点等参元implicit noneinteger (kind(1),parameter :ikind=(kind(1)integer (kind(1),parameter :rkind=(kind(0.d0)type : typ_Kcolreal(rkind),pointer : Row(:)end type typ_Kcoltype : typ_GValue !总体控制变量integer(ikind) : NNode, NElem, NLoad, NMat, NSupportinteger(ikind) : NGlbDOF !整体自由度总数integer(ikind) : NGENS, NodeDOF,ElemNodeNointeger(ikind) : NIntend type typ_GValuetype Typ_Node !定义节点类型real(rkind) : coord(2) !节点坐标integer(ikind) : GDOF(2) !整体自由度编码real(rkind) : DISP(2) !节点位移real(rkind) : dDISP(2) !节点位移增量real(rkind) : dForce(2) !节点不平衡力end type typ_Node!=Type typ_IntPoint !定义积分点参数real(rkind) : EPS(3) !应变real(rkind) : SIG(3) !应力real(rkind) : D(3,3) !本构矩阵real(rkind) : B(3,16) !几何矩阵real(rkind) : DETJ !雅克比行列式end type Typ_IntPointtype Typ_Rect8 !定义实体单元integer(ikind) : NodeNo(8) !节点编号real(rkind) : E !弹性模量real(rkind) : u !泊松比real(rkind) : t !单元厚度real(rkind) : EK(16,16) !单元刚度矩阵type(typ_intpoint) : IntP(9) !积分点!.end type Typ_Rect8type Typ_Load ! 定义荷载类型integer(ikind) : NodeNo !荷载作用节点编号real(rkind) : Value(2) !荷载大小end type Typ_Loadtype Typ_Support ! 定义支座类型integer(ikind) : NodeNo !支座节点编号integer(ikind) : DOF !支座约束自由度real(rkind) : Value !支座位移大小end type Typ_Supportcontainssubroutine Get_Elem_K(GValue,Elem,Node) ! 核心程序，得到单元的刚度矩阵implicit nonetype(typ_GValue) : GValuetype(typ_Node) : Node(:)type(typ_Rect8) : Elem(:)real*8 : InvJ(2,2) ! 雅克比矩阵及其逆矩阵real*8 : Coord(8,2),IntPCoord(9,2),IntPwt(9)real*8 : dN(2,8) !节点坐标，积分点坐标，权函数，形函数导数integer : I,ElemNocall Get_IntP_Prop(IntPCoord,IntPWt) ! 计算积分点信息do ElemNo=1,size(Elem)Elem(ElemNo)%EK=0.0do I=1, GValue%ElemNodeNo; !得到节点坐标coord(I,:)=Node(Elem(ElemNo)%NodeNo(i)%coord; end doElem(ElemNo)%EK=0d0 DO I=1,size(Elem(ElemNo)%IntP)! 计算本构矩阵call Get_D(Elem(ElemNo)%IntP(I)%D,Elem(ElemNo)%E,Elem(ElemNo)%u) ! 计算形函数对局部坐标导数call Get_dN_dxi(dN,IntPCoord(i,1),IntPCoord(i,2)InvJ=matmul(dN, Coord)! 得到雅克比行列式Elem(ElemNo)%IntP(I)%DETJ=InvJ(1,1)*InvJ(2,2)-InvJ(1,2)*InvJ(2,1)InvJ=matinv(InvJ); ! 对雅克比行列式求逆dN = matmul(InvJ,dN) ! 得到形函数对整体坐标的导数call Get_B (Elem(ElemNo)%IntP(I)%B,dN) ! 得到几何矩阵Elem(ElemNo)%EK=Elem(ElemNo)%EK+&matmul(matmul(transpose(Elem(ElemNo)%IntP(I)%B),&Elem(ElemNo)%IntP(I)%D),Elem(ElemNo)%IntP(I)%B)*&Elem(ElemNo)%IntP(I)%DetJ*IntPWt(i)*Elem(ElemNo)%tend doend doreturnend subroutinesubroutine Get_D(D,e,v) ! 得到本构矩阵，平面应力implicit nonereal*8 : e,vreal*8 : D (:,:)D=0.d0D(1,1)=1D0; D(2,2)=1d0D(1,2)=v; D(2,1)=vD(3,1:2)=0d0; D(1:2,3)=0d0;D(3,3)=(1.d0-v)/2.d0;D=E/(1.d0-v*v)*D;returnend subroutine Get_D subroutine Get_IntP_Prop(IntPCoord,IntPWt) ! 得到高斯积分点的参数implicit none real*8 :IntPCoord(:,:),IntPWt(:) ! 高斯积分点坐标，高斯积分点权重real*8 : z,x(3),y(3),w(3)integer : i,jz=.2d0*sqrt(15.d0) x=(/-1.d0,0.d0,1.d0/); y=(/1.d0,0.d0,-1.d0/)w=(/5.d0/9.d0,8.d0/9.d0,5.d0/9.d0/)do i=1,3do j=1,3IntPCoord(i-1)*3+j,1)=x(j)*zIntPCoord(i-1)*3+j,2)=y(i)*zIntPWt(i-1)*3+j)=w(i)*w(j)end doend doreturnend subroutine Get_IntP_Propsubroutine Get_dN_dxi(dN_dxi,xi,eta) ! 得到形函数对局部坐标系的导数implicit nonereal*8 : dN_dxi(:,:), xi,etareal*8: x1, x2, x3, x4, x5 ,x6 ,x7 ,x8 x1=.25*(1.-eta); x2=.25*(1.+eta); x3=.25*(1.-xi); x4=.25*(1.+xi) dN_dxi(1,1)=x1*(2.*xi+eta)dN_dxi(1,2)=-8.*x1*x2dN_dxi(1,3)=x2*(2.*xi-eta)dN_dxi(1,4)=-4.*x2*xidN_dxi(1,5)=x2*(2.*xi+eta)dN_dxi(1,6)=8.*x2*x1dN_dxi(1,7)=x1*(2.*xi-eta)dN_dxi(1,8)=-4.*x1*xidN_dxi(2,1)=x3*(xi+2.*eta)dN_dxi(2,2)=-4.*x3*etadN_dxi(2,3)=x3*(2.*eta-xi)dN_dxi(2,4)=8.*x3*x4dN_dxi(2,5)=x4*(xi+2.*eta)dN_dxi(2,6)=-4.*x4*etadN_dxi(2,7)=x4*(2.*eta-xi)dN_dxi(2,8)=-8.*x3*x4 returnend subroutine Get_dN_dxisubroutine Get_B(B,dN_dx) ! 得到几何矩阵implicit nonereal*8 : B(:,:), dN_dx(:,:) integer:k,l,m,n , ih,nod; real*8 : x,y,zB=0. do m=1,8k=2*m; l=k-1; x=dN_dx(1,m); y=dN_dx(2,m)B(1,l)=x; B(3,k)=x; B(2,k)=y; B(3,l)=yend do returnend subroutine Get_Bfunction matinv(A) result (B) ! 计算2x2逆矩阵real(rkind) ,intent (in): A(:,:)real(rkind) , pointer:B(:,:)real(rkind) : xinteger : NN=size(A,dim=2)allocate(B(N,N)B(1,1)=A(2,2)B(2,2)=A(1,1)B(1,2)=-A(1,2)B(2,1)=-A(2,1)x=A(1,1)*A(2,2)-A(1,2)*A(2,1)B=B/xreturnend function matinvend modulemodule Data_Input ! 数据输入模块 use Elem_Rect8 implicit nonecontainssubroutine DataInput(GValue,Node,Elem,Load,Support)type(typ_GValue) : GValuetype(typ_Node),pointer : Node(:)type(typ_Rect8),pointer : Elem(:)type(typ_Load),pointer : Load(:)type(typ_Support),pointer : Support(:)call ReadGValue(GValue) !读入整体控制变量call ReadNode(GValue,Node) !读入节点坐标信息call ReadElem(GValue,Elem) !读入单元原始信息call ReadMaterial(GValue,Elem)call ReadLoad(GValue,Load) !读入荷载信息 call ReadSupport(GValue,Support) !读入支座信息returnend subroutine DataInputsubroutine ReadGValue(GValue) type(typ_GValue) : GValue GValue%NGENS=3GValue%NodeDOF=2GValue%ElemNodeNo=8GValue%NInt=3returnend subroutine ReadGValuesubroutine ReadNode(GValue,Node) ! 读入结点坐标type(typ_GValue) : GValuetype(typ_Node),pointer : Node(:)integer(ikind) : I,Jopen(55,file=node.txt)read(55,*) GValue%NNodeGValue%NGlbDOF=2*GValue%NNodeallocate(Node(GValue%NNode)do I=1, GValue%NNoderead(55,*) J,Node(I)%Coord(1:2)end doclose(55)write(*,*) 结点信息读入完毕returnend subroutine ReadNodesubroutine ReadElem(GValue,Elem) ! 读入单元信息type(typ_GValue) : GValuetype(typ_Rect8),pointer : Elem(:)integer(ikind) : I,J,Kinteger(ikind) : TransNode(8),N(8),N1(8) ! 结点坐标变换TransNode=(/1,7,5,3,8, 6,4, 2/) ! 注意单元节点编号匹配open(55,file=Element.txt)read(55,*) GValue%NElemallocate(Elem(GValue%NElem)do I=1,GValue%NElemread(55,*) J, Ndo J=1,size(N)K=TransNode(J)N1(K)=N(J) end doElem(I)%NodeNo=N1end doclose(55)write(*,*) 单元信息读入完毕returnend subroutine ReadElemsubroutine ReadMaterial(GValue,Elem) ! 读入材料信息type(typ_GValue) : GValuetype(typ_Rect8),pointer : Elem(:)integer (ikind) : NElemreal(rkind) : E, muinteger(ikind) : I,Jinteger(ikind),allocatable : K(:)open(55,file=Material.txt)read(55,*) E, muElem(:)%E=EElem(:)%u=mu Elem(:)%t=0.5close(55)write(*,*) 材料信息读入完毕returnend subroutine ReadMaterialsubroutine ReadLoad(GValue,Load) ! 读入荷载信息type(typ_GValue) : GValuetype(typ_Load),pointer : Load(:)integer(ikind) : I,Jopen(55,file=Load.txt)read(55,*) GValue%NLoadallocate(Load(GValue%NLoad)do I=1, GValue%NLoadread(55,*) J, Load(I)%NodeNo,Load(I)%Value(1:2)end doclose(55)write(*,*) 荷载信息读入完毕returnend subroutine ReadLoadsubroutine ReadSupport(GValue,Support) ! 读入支座信息type(typ_GValue) : GValuetype(typ_Support),pointer : Support(:)integer(ikind) : I,Jopen(55,file=Support.txt)read(55,*) GValue%NSupportallocate(Support(GValue%NSupport)do I=1, GValue%NSupportread(55,*) J, Support(I)%NodeNo, Support(I)%DOF, Support(I)%Valueend doclose(55)write(*,*) 支座信息读入完?quot;returnend subroutine ReadSupportend modulemodule Mat_solve ! 矩阵求解模块use Elem_Rect8implicit noneinteger : NGlbDOFcontainssubroutine Matsolve(GValue,Elem,Node,Load, Support)type(typ_GValue) : GValuetype(typ_Node),pointer : Node(:)type(typ_Rect8),pointer : Elem(:)type(typ_Load),pointer : Load(:)type(typ_Support),pointer : Support(:)type(typ_Kcol),allocatable : Kcol(:)real(rkind),allocatable : GLoad(:), GDisp(:)real(rkind) : Penatlyinteger(ikind) : BandWidth(2*GValue%NNode)integer(ikind) : ELocVec(16)integer(ikind) : I,J,Kinteger(ikind) : MinDOFNumNGlbDOF=2*GValue%NNodePenatly=1.0 ! 罚函数allocate(GLoad(NGlbDOF)allocate(GDisp(NGlbDOF)!查找带宽do I=1,NGlbDOFBandWidth(I)=Iend dodo I=1,GValue%NElemdo J=1,size(Elem(I)%NodeNo)ELocVec(J*2-1)=2*Elem(I)%NodeNo(J)-1ELocVec(J*2 )=2*Elem(I)%NodeNo(J) end doMinDOFNum=minval(ELocVec)do J=1,size(ElocVec)BandWidth(ELocVec(J)=min(MinDOFNum,BandWidth(ELocVec(J) end doend dowrite(*,*) 完成带宽查找allocate(Kcol(NGlbDOF)do I=1, NGlbDOFallocate(Kcol(I)%Row(BandWidth(I):I)Kcol(I)%Row=0.d0end dodo I=1, GValue%NElemdo J=1,size(Elem(I)%NodeNo)ELocVec(J*2-1)=2*Elem(I)%NodeNo(J)-1ELocVec(J*2 )=2*Elem(I)%NodeNo(J) end dodo J=1,size(ElocVec)do K=1,Jif(ELocVec(J)ELocVec(K) thenKcol(ELocVec(J)%row(ELocVec(K)=&Kcol(ELocVec(J)%row(ELocVec(K)+Elem(I)%EK(J,K)elseKcol(ELocVec(K)%row(ELocVec(J)=&Kcol(ELocVec(K)%row(ELocVec(J)+Elem(I)%EK(J,K)end ifend doend doPenatly=max(Penatly,maxval(abs(Elem(I)%EK)end dowrite(*,*) 完成总刚集成GLoad=0.d0do I=1, size(Load)J=Load(I)%NodeNoGLoad(J*2-1:J*2)=GLoad(J*2-1:J*2)+Load(I)%Value(:) end dodo I=1,GValue%NSupportJ=2*(Support(I)%NodeNo-1)+Support(I)%DOFGLoad(J)=GLoad(J)+Support(I)%Value*Penatly*1.0D8Kcol(J)%Row(J)=KCol(J)%Row(J)+Penatly*1.0d8end dowrite(*,*) 完成支座集成call BandSolv(Kcol,Gload,GDisp)call Get_Node_dDisp(GValue,Node,GDisp)do I=1,NGlbDOFdeallocate(Kcol(I)%Row)end dodeallocate(Kcol)deallocate(GDisp)deallocate(GLoad)returnend subroutinesubroutine Get_Node_dDisp(GValue,Node,GDisp) ! 从整体位移向量得到结点位移向量type(typ_GValue) : GValuetype(typ_Node) : Node(:)real(rkind) : GDisp(:)integer(ikind) : I,Jdo I=1, GValue%NNodeNode(I)%dDisp(1:2)=GDisp(I*2-1:I*2)Node(I)%Disp=Node(I)%Disp+Node(I)%dDispend doreturnend subroutine Get_Node_dDisp!-subroutine BandSolv(Kcol,GLoad,diag)!-type (typ_Kcol) : Kcol(:);real(rkind) : GLoad(:),diag(:);integer : row1,ncol,row,j,iereal(rkind) : sncol=NGlbDOFdiag(1:ncol)=(/(Kcol(j)%row(j),j=1,ncol)/)do j=2,ncolrow1=lbound(Kcol(j)%row,1)do ie=row1,j-1row=max(row1,lbound(Kcol(ie)%row,1)s=sum(diag(row:ie-1)*Kcol(ie)%row(row:ie-1)*Kcol(j)%row(row:ie-1)Kcol(j)%row(ie)=(Kcol(j)%row(ie)-s)/diag(ie)end dos=sum(diag(row1:j-1)*Kcol(j)%row(row1:j-1)*2)diag(j)=diag(j)-send dodo ie=2,ncolrow1=lbound(Kcol(ie)%row,dim=1)GLoad(ie)=GLoad(ie)-sum(Kcol(ie)%row(row1:ie-1)*GLoad(row1:ie-1)end doGLoad(:)=GLoad(:)/diag(:)do j=ncol,2,-1row1=lbound(Kcol(j)%row,dim=1)GLoad(row1:j-1)=GLoad(row1:j-1)-GLoad(j)*