平面8節(jié)點(diǎn)等參元完整程序.doc

平面8節(jié)點(diǎn)等參元完整程序module Elem_Rect8 ! 八節(jié)點(diǎn)等參元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 !整體自由度總數(shù)integer(ikind) : NGENS, NodeDOF,ElemNodeNointeger(ikind) : NIntend type typ_GValuetype Typ_Node !定義節(jié)點(diǎn)類型real(rkind) : coord(2) !節(jié)點(diǎn)坐標(biāo)integer(ikind) : GDOF(2) !整體自由度編碼real(rkind) : DISP(2) !節(jié)點(diǎn)位移real(rkind) : dDISP(2) !節(jié)點(diǎn)位移增量real(rkind) : dForce(2) !節(jié)點(diǎn)不平衡力end type typ_Node!=Type typ_IntPoint !定義積分點(diǎn)參數(shù)real(rkind) : EPS(3) !應(yīng)變r(jià)eal(rkind) : SIG(3) !應(yīng)力real(rkind) : D(3,3) !本構(gòu)矩陣real(rkind) : B(3,16) !幾何矩陣real(rkind) : DETJ !雅克比行列式end type Typ_IntPointtype Typ_Rect8 !定義實(shí)體單元integer(ikind) : NodeNo(8) !節(jié)點(diǎn)編號(hào)real(rkind) : E !彈性模量real(rkind) : u !泊松比real(rkind) : t !單元厚度real(rkind) : EK(16,16) !單元?jiǎng)偠染仃噒ype(typ_intpoint) : IntP(9) !積分點(diǎn)!.end type Typ_Rect8type Typ_Load ! 定義荷載類型integer(ikind) : NodeNo !荷載作用節(jié)點(diǎn)編號(hào)real(rkind) : Value(2) !荷載大小end type Typ_Loadtype Typ_Support ! 定義支座類型integer(ikind) : NodeNo !支座節(jié)點(diǎn)編號(hào)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) !節(jié)點(diǎn)坐標(biāo)，積分點(diǎn)坐標(biāo)，權(quán)函數(shù)，形函數(shù)導(dǎo)數(shù)integer : I,ElemNocall Get_IntP_Prop(IntPCoord,IntPWt) ! 計(jì)算積分點(diǎn)信息do ElemNo=1,size(Elem)Elem(ElemNo)%EK=0.0do I=1, GValue%ElemNodeNo; !得到節(jié)點(diǎn)坐標(biāo)coord(I,:)=Node(Elem(ElemNo)%NodeNo(i)%coord; end doElem(ElemNo)%EK=0d0 DO I=1,size(Elem(ElemNo)%IntP)! 計(jì)算本構(gòu)矩陣call Get_D(Elem(ElemNo)%IntP(I)%D,Elem(ElemNo)%E,Elem(ElemNo)%u) ! 計(jì)算形函數(shù)對(duì)局部坐標(biāo)導(dǎo)數(shù)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); ! 對(duì)雅克比行列式求逆dN = matmul(InvJ,dN) ! 得到形函數(shù)對(duì)整體坐標(biāo)的導(dǎo)數(shù)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) ! 得到本構(gòu)矩陣，平面應(yīng)力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) ! 得到高斯積分點(diǎn)的參數(shù)implicit none real*8 :IntPCoord(:,:),IntPWt(:) ! 高斯積分點(diǎn)坐標(biāo)，高斯積分點(diǎn)權(quán)重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) ! 得到形函數(shù)對(duì)局部坐標(biāo)系的導(dǎo)數(shù)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) ! 計(jì)算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 ! 數(shù)據(jù)輸入模塊 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) !讀入節(jié)點(diǎn)坐標(biāo)信息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) ! 讀入結(jié)點(diǎn)坐標(biāo)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(*,*) 結(jié)點(diǎn)信息讀入完畢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) ! 結(jié)點(diǎn)坐標(biāo)變換TransNode=(/1,7,5,3,8, 6,4, 2/) ! 注意單元節(jié)點(diǎn)編號(hào)匹配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 ! 罰函數(shù)allocate(GLoad(NGlbDOF)allocate(GDisp(NGlbDOF)!查找?guī)抎o 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) ! 從整體位移向量得到結(jié)點(diǎn)位移向量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)*Kcol(