程序结构力学大作业

程序结构力学大作业结81 孙玉进该程序可计算任意平面结构任意阶频率（包括重频），以及任意阶振型（包括重频对应正交振型）。计算时，振型不包括不包括单元固端型振型。!*module NumKind!*implicit noneinteger (kind(1),parameter : ikind = kind(1), rkind = kind(0.D0)real (rkind),parameter : Zero = 0._rkind, One = 1._rkind, Two = 2._rkind, &Three= 3._rkind, Four = 4._rkind, Five = 5._rkind, &Six = 6._rkind, Seven= 7._rkind, Eight = 8._rkind, &Nine = 9._rkind, Ten =10._rkind, Twelve=12._rkindend module NumKind!*module TypeDef!*use NumKindimplicit nonetype : typ_Jointreal (rkind) : x,yinteger (ikind) : GDOF(3)end type typ_Jointtype : typ_Elementinteger (ikind) : JointNo(2),GlbDOF(6)real (rkind) : Length,CosA,SinA,EI,EA,massend type typ_Elementtype : typ_JointLoadinteger (ikind) : JointNo,LodDOFreal (rkind) : LodValend type typ_JointLoadtype : typ_ElemLoadinteger (ikind) :ElemNo,Indxreal (rkind) : Pos,LodValend type typ_ElemLoadcontains!=subroutine SetElemProp (Elem, Joint)!=type (typ_Element),intent(in out) : Elem(:)type (typ_Joint),intent(in) : Joint(:)integer(ikind):ie,NElemreal(rkind):x1,x2,y1,y2NElem=size(Elem,dim=1)do ie=1,NElemx1=Joint(Elem(ie)%JointNo(1)%xy1=Joint(Elem(ie)%JointNo(1)%yx2=Joint(Elem(ie)%JointNo(2)%xy2=Joint(Elem(ie)%JointNo(2)%yElem(ie)%Length=sqrt(x1-x2)*2+(y1-y2)*2)Elem(ie)%CosA=(x2-x1)/Elem(ie)%LengthElem(ie)%SinA=(y2-y1)/Elem(ie)%LengthElem(ie)%GlbDOF(1:3)=Joint(Elem(ie)%JointNo(1)%GDOF(:)Elem(ie)%GlbDOF(4:6)=Joint(Elem(ie)%JointNo(2)%GDOF(:)end doreturnend subroutine SetElemProp!=subroutine TransMatrix (ET, CosA,SinA)!=real(rkind),intent(out) : ET(:,:)real(rkind),intent(in) : CosA,SinA! ET could be 2x2, 3x3 or 6x6 depending size(ET)ET = ZeroET(1,1:2) = (/ CosA, SinA /)ET(2,1:2) = (/-SinA, CosA /)if (size(ET,1) 2) ET(3,3) = Oneif (size(ET,1) 3) ET(4:6,4:6) = ET(1:3,1:3)returnend subroutine TransMatrixend module TypeDef!*module BandMat!*use NumKinduse TypeDef,only : typ_Element ! 仅用该模块中的typ_Elementimplicit noneprivate ! 默认所有的数据和过程为私有，增强封装性public : SetMatBand, DelMatBand, VarBandSolvtype,public : typ_Kcolreal (rkind),pointer : row(:)end type typ_Kcolcontains!=subroutine SetMatBand (Kcol, Elem)!=! .6-4-2type (typ_KCol),intent(in out) : Kcol(:)type (typ_Element),intent(in) : Elem(:)integer (ikind) : minDOF,ELocVec(6)integer (ikind) : ie,j,NGlbDOF,NEleminteger (ikind) : row1(size(Kcol,dim=1)! row1是自动数组，子程序结束后将自动释放内存空间NGlbDOF = size(Kcol,1)NElem = size(Elem,1)row1 = NGlbDOF ! 先设始行码为大数! 确定各列始行码，放在数组row1(:)中do ie=1,NElemELocVec = Elem(ie)%GlbDOFminDOF = minval (ELocVec,mask = ELocVec 0)where (ELocVec 0)row1(ELocVec) = min(row1(ELocVec), minDOF)end whereend do! 为各列的半带宽分配空间并初始化do j=1,NGlbDOFallocate ( Kcol(j)%row(row1(j):j) )Kcol(j)%row = Zero ! 清零end doreturnend subroutine SetMatBand!=subroutine DelMatBand (Kcol)!=!.6-5-5type (typ_KCol), intent(in out) : Kcol(:)integer (ikind) : j,NGlbDOFNGlbDOF = size(Kcol,1)do j=1,NGlbDOFdeallocate ( Kcol(j)%row )end doreturnend subroutine DelMatBand!=subroutine VarBandSolv (Disp, Kcol,GLoad)!=type (typ_KCol), intent(in out) : Kcol(:)real (rkind), intent(out) : Disp(:)real (rkind), intent(in) : GLoad(:)integer (ikind) : i,j,k,row1j,row_1,NColreal(rkind) : Diag(size(Kcol,dim=1)real(rkind) :s!.6-5-2NCol=size(Kcol,1)Diag(1:NCol)=(/(Kcol(j)%row(j),j=1,NCol)/)do j=2,NColrow1j=lbound(Kcol(j)%row,1)do i=row1j,j-1row_1=max(row1j,lbound(Kcol(i)%row,1)k=i-1s= sum(Diag(row_1:k)*Kcol(i)%row(row_1:k)*Kcol(j)%row(row_1:k)Kcol(j)%row(i)=(Kcol(j)%row(i)-s)/Diag(i)end dos=sum(Diag(row1j:j-1)*Kcol(j)%row(row1j:j-1)*2)Diag(j)=Diag(j)-send doDisp(:) = GLoad(:)!. 6-5-3节的代码：其中GP换为Disp do j=2,NColrow1j=lbound(Kcol(j)%row,1)Disp(j)=Disp(j)-sum(Kcol(j)%row(row1j:j-1)*Disp(row1j:j-1)end do!. 6-5-4节的代码：其中GP换为Disp Disp(:)=Disp(:)/Diag(:)do j=NCol,1,-1row1j=lbound(Kcol(j)%row,1)Disp(row1j:j-1)=Disp(row1j:j-1)-Disp(j)*Kcol(j)%row(row1j:j-1)end doreturnend subroutine VarBandSolvend module BandMat!*module DispMethod!*use NumKinduse TypeDefuse BandMatimplicit nonecontains!=subroutine SolveDisp(Disp,Elem,Joint,JLoad,ELoad)!=real(rkind),intent(out) : Disp(:)type(typ_Element),intent(in) : Elem(:)type(typ_Joint) ,intent(in) : Joint(:)type(typ_JointLoad),intent(in) : JLoad(:)type(typ_Elemload),intent(in) : ELoad(:)type(typ_Kcol),allocatable : Kcol(:) real(rkind),allocatable : GLoad(:) integer(ikind) : NGlbDOFNGlbDOF=size(Disp,dim=1)allocate(GLoad(NGlbDOF)allocate(Kcol(NGlbDOF)call SetMatBand(Kcol,Elem)call GLoadVec(GLoad,Elem,JLoad,ELoad,Joint)call GStifMat(Kcol,Elem)call VarBandSolv(Disp,Kcol,GLoad)call DelMatBand(Kcol)deallocate(GLoad) returnend subroutine SolveDisp!=subroutine GStifMat(Kcol,Elem)!=type(typ_Kcol),intent(in out) : Kcol(:)type(typ_Element),intent(in) : Elem(:)! .6-4-3integer (ikind) : ie,j,JGDOF,NElemreal (rkind) : EK(6,6),ET(6,6)integer (ikind) : ELocVec(6)NElem=size(Elem,1)do ie=1,NElemcall EStifMat (EK,Elem(ie)%Length,Elem(ie)%EI,Elem(ie)%EA)call TransMatrix(ET,Elem(ie)%CosA,Elem(ie)%SinA)EK=matmul(transpose(ET),matmul(EK,ET)ELocVec=Elem(ie)%GlbDOFdo j=1,6JGDOF=ELocVec(j)if(JGDOF=0) cyclewhere (ELocVec0)Kcol(JGDOF)%row(ELocVec)=Kcol(JGDOF)%row(ELocVec)+EK(:,j)end whereend doend doreturnend subroutine GStifMat!=subroutine GLoadVec(GLoad,Elem,JLoad,ELoad,Joint)!=real (rkind),intent(out) : GLoad(:)type (typ_Element), intent(in) : Elem(:)type (typ_Joint), intent(in) : Joint(:)type (typ_JointLoad), intent(in) :JLoad(:)type (typ_ElemLoad), intent(in) :ELoad(:)integer (ikind) : ie,NJLoad,NELoadreal (rkind) : F0(6),ET(6,6)NJLoad=size(JLoad,dim=1)NELoad=size(ELoad,dim=1)GLoad(:)=Zerodo ie=1,NJLoadGLoad(Joint(JLoad(ie)%JointNo)%GDOF(JLoad(ie)%LodDOF)=GLoad(Joint(JLoad(ie)%JointNo)%GDOF(JLoad(ie)%LodDOF)+JLoad(ie)%LodValend dodo ie=1,NELoadcall EFixendF (F0, ELoad(ie)%Indx,ELoad(ie)%Pos,ELoad(ie)%LodVal,Elem(ELoad(ie)%ElemNo)call TransMatrix(ET,Elem(ELoad(ie)%ElemNo)%CosA,Elem(ELoad(ie)%ElemNo)%SinA)where (Elem(ELoad(ie)%ElemNo)%GlbDOF0)GLoad(Elem(ELoad(ie)%ElemNo)%GlbDOF)=GLoad(Elem(ELoad(ie)%ElemNo)%GlbDOF)-matmul(transpose(ET),F0(:)end where end doreturnend subroutine GLoadVec!=subroutine EStifMat(EK,Elen,EI,EA)!=real(rkind),intent(out) : EK(:,:)real(rkind),intent(in) : Elen,EI,EAinteger(ikind) : i,jEK=ZeroEK(1,1) = EA/ElenEK(1,4) = -EA/ElenEK(2,2) = Twelve*EI/Elen*ThreeEK(2,3) = Six*EI/Elen*TwoEK(2,5) = -Twelve*EI/Elen*ThreeEK(2,6) = Six*EI/Elen*TwoEK(3,3) = Four*EI/ElenEK(3,5) = -Six*EI/Elen*TwoEK(3,6) = Two*EI/ElenEK(4,4) = EA/ElenEk(5,5) = Twelve*EI/Elen*ThreeEK(5,6) = -Six*EI/Elen*TwoEK(6,6) = Four*EI/Elendo j=1,6do i=j+1,6EK(i,j)=EK(j,i)end doend doreturnend subroutine EStifMat!=subroutine EFixendF(F0,Indx,a,q,Elem)!=real(rkind), intent(out): F0(:)real (rkind), intent(in) :a,qinteger (ikind), intent(in): Indxtype (typ_Element), intent(in) : Elemreal :ll=Elem%lengthselect case (Indx)case(1)F0(1:3)=(/Zero,-q*(a*l)/Two*(Two-Two*a*Two+a*Three),-q*(a*l)*Two/Twelve*(Six-Eight*a+Three*a*Two)/)F0(4:6)=(/Zero,-q*(a*l)*a*Two/Two*(Two-a),q*(a*l)*Two*a/Twelve*(Four-Three*a)/)case(2)F0(1:3)=(/Zero,-q*(One-Two*a+a*Two)*(One+Two*a),-q*(a*l)*(One-Two*a+a*Two)/)F0(4:6)=(/Zero,-q*a*Two*(Three-Two*a),q*a*Two*(One-a)*l/)case(3)if(aOne/Two)thenF0(1)=Elem%EA*q/lF0(4)=-F0(1)elseF0(1)=-Elem%EA*q/lF0(4)=-F0(1)end ifF0(2) = ZeroF0(3) = ZeroF0(5) = ZeroF0(6) = Zero case(4)if(a1)thenMFreq=MFreq-1Freq(k-FreqStart+1)=Freq(k-FreqStart)cycleend ifif(kFreqStart) freq_l=freq(k-FreqStart)docall CJ0(J01,freq_l,Elem)call CJk(Jk,freq_l,Kcol,Elem)J_l=J01+Jkif(J_l=k) exitfreq_l=freq_ufreq_u=freq_u*Twoend dodofreqm=(freq_l+freq_u)/Twocall CJ0(J0,freqm,Elem)call CJk(Jk,freqm,Kcol,Elem)Jm=J0+Jkif(Jm=k)thenfreq_u=freqmelsefreq_l=freqmend ifif(freq_u-freq_l)Tol*(One+freq_l) thencall CJ0(J0,freq_u,Elem)call CJk(Jk,freq_u,Kcol,Elem)Jm=J0+JkMFreq=Jm-k+1exitend ifend doFreq(k+1-FreqStart)=(freq_u+freq_l)/2end do returnend subroutine GetFreq!-subroutine Mode(Elem,Kcol,FreqStart,NFreq,Tol)!获取振型!-type(typ_Element),intent(in) :Elem(:)type(typ_Kcol) ,intent(in out) :Kcol(:)real(rkind) ,intent(in) :Tol integer(ikind) ,intent(in) :FreqStart,NFreqreal(rkind) ,allocatable :det0(:),det(:),GLoad(:),diag(:),GLoad2(:),Freq(:),Mdet(:,:) !记录重频振型integer(ikind) ,allocatable :MFreq(:)!记录重频已出现次数real(rkind) :de,s,freq0,arf,EDisp(6)integer(ikind) :NGlbDOF,NElem,i,j,row1j,row_1,NCol,k,t,rNGlbDOF=size(Kcol,1)NCol=size(Kcol,1)NElem=size(Elem,1)r=min(NFreq,NGlbDOF)allocate(det0(NGlbDOF)allocate(det(NGlbDOF)allocate(GLoad(NGlbDOF)allocate(diag(NGlbDOF)allocate(MFreq(NFreq+1)allocate(GLoad2(NGlbDOF)allocate(Freq(NFreq)MFreq=Zerocall GetFreq(Freq,Elem,Kcol,FreqStart,NFreq,Tol)do i=1,NFreq-1 !计算重频对应的已出现次数if(Freq(i+1)-Freq(i)0) then !重频振型正交化do t=1,MFreq(k)call EGLoadVec(GLoad,Elem,det,Freq(k)call EGLoadVec(GLoad2,Elem,Mdet(1:NGlbDOF,t),Freq(k)arf=dot_product(Mdet(1:NGlbDOF,t),GLoad)/dot_product(Mdet(1:NGlbDOF,t),GLoad2)det=det-Mdet(1:NGlbDOF,t)*arfend doend if!逆幂迭代call SetMatBand(Kcol,Elem)call GStifMat(Kcol,Elem,Freq(k)diag(:)=(/(Kcol(j)%row(j),j=1,NCol)/)do j=2,NColrow1j=lbound(Kcol(j)%row,1)do i=row1j,j-1row_1=max(row1j,lbound(Kcol(i)%row,1)s=sum(diag(row_1:i-1)*Kcol(i)%row(row_1:i-1)*Kcol(j)%row(row_1:i-1)Kcol(j)%row(i)=(Kcol(j)%row(i)-s)/diag(i)end dos=sum(diag(row1j:j-1)*Kcol(j)%row(row1j:j-1)*2)diag(j)=diag(j)-send do do j=2,NColrow1j=lbound(Kcol(j)%row,1)GLoad(j)=GLoad(j)-sum(Kcol(j)%row(row1j:j-1)*GLoad(row1j:j-1)end doGLoad(:)=GLoad(:)/Diag(:)do j=NCol,1,-1row1j=lbound(Kcol(j)%row,1)GLoad(row1j:j-1)=GLoad(row1j:j-1)-GLoad(j)*Kcol(j)%row(row1j:j-1)end dode=GLoad(1)do i=2,NGlbDOF !标准化if(abs(GLoad(i)abs(de) de=GLoad(i) end dodo i=1,NGlbDOFGLoad(i)=GLoad(i)/deend doif(dot_product(GLoad-det0,GLoad-det0)0) Mdet(1:NGlbDOF,MFreq(k+1)=GLoad !记录重频振型，用于后续正交化exitend ifdet0=GLoadfreq0=Freq(k)Freq(k)=Freq(k)-1/de !牛顿外插频率if(Freq(k)freq0+Tol*(One+freq0) Freq(k)=freq0 !牛顿导护if(Freq(k)freq0-Tol*(One+freq0) Freq(k)=freq0end dowrite(8,*) Freq(k),0 do t=1,NElem call ElemDisp2(EDisp,GLoad,Elem(t)write(8,*) EDisp end doend doEDisp=0do k=r+1,NFreqwrite(8,*) Freq(k),0do t=1,NElemwrite(8,*) EDispend doend dodeallocate(diag)deallocate(det0)deallocate(GLoad)call DelMatBand (Kcol)returnend subroutine Mode!计算单元固端频率数subroutine CJ0(J0,freq,Elem)integer(ikind) ,intent(out) :J0real(rkind) ,intent(in) :freqtype(typ_Element) ,intent(in) :Elem(:)integer(ikind) :NElem,Ja,Jb,i,jreal(rkind) :nu,lambda,sgNElem=size(Elem,dim=1)J0=0do i=1,NElemnu = freq*Elem(i)%Length*sqrt(Elem(i)%mass/Elem(i)%EA)Ja = int(nu/Pi)lambda= Elem(i)%Length*sqrt(sqrt(freq*2*Elem(i)%mass/Elem(i)%EI)sg = sign(One,One-cosh(lambda)*cos(lambda)j = int(lambda/Pi)Jb = j-(1-(-1)*j*sg)/2J0 = J0+Ja+Jbend doreturnend subroutine CJ0!计算Jksubroutine CJk(Jk,freq,Kcol,Elem)integer(ikind) ,intent(out) :Jkreal(rkind) ,intent(in) :freqtype(typ_Kcol) ,intent(in out) :Kcol(:)type(typ_Element),intent(in) :Elem(:)integer(ikind) :NGlbDOFreal(rkind),allocatable :diag(:)integer(ikind) :i,j,row1j,row_1,NColreal(rkind) :sNGlbDOF=size(Kcol,1)allocate(diag(NGlbDOF)call SetMatBand(Kcol,Elem)call GStifMat(Kcol,Elem,freq)NCol=size(Kcol,1)diag(:)=(/(Kcol(j)%row(j),j=1,NCol)/)do j=2,NColrow1j=lbound(Kcol(j)%row,1)do i=row1j,j-1row_1=max(row1j,lbound(Kcol(i)%row,1)s=sum(diag(row_1:i-1)*Kcol(i)%row(row_1:i-1)*Kcol(j)%row(row_1:i-1)Kcol(j)%row(i)=(Kcol(j)%row(i)-s)/diag(i)end dos=sum(diag(row1j:j-1)*Kcol(j)%row(row1j:j-1)*2)diag(j)=diag(j)-send do Jk=count(mask=diagZero,dim=1)deallocate(diag)returnend subroutine CJk!整体动力刚度矩阵subroutine GStifMat(Kcol,Elem,freq)type(typ_Kcol),intent(in out) : Kcol(:)type(typ_Element),intent(in) : Elem(:)real(rkind),intent(in) : freqinteger(ikind) : ie,j,JGDOF,NElemreal(rkind) : EK(6,6),ET(6,6)integer(ikind) : ELocVec(6)NElem=size(Elem,1)do ie=1,NElemcall EStifMat(EK,Elem,ie,freq)call TransMatrix(ET,Elem(ie)%CosA,Elem(ie)%SinA)EK=matmul(transpose(ET),matmul(EK,ET)ELocVec=Elem(ie)%GlbDOFdo j=1,6JGDOF=ELocVec(j)if(JGDOF=0) cyclewhere (ELocVec0)Kcol(JGDOF)%row(ELocVec)=Kcol(JGDOF)%row(ELocVec)+EK(:,j)end whereend doend doreturnend subroutine GStifMat!等效整体荷载向量，用于逆幂迭代!=subroutine EGLoadVec(GLoad,Elem,det,freq)!=real (rkind),intent(out) : GLoad(:)type (typ_Element), intent(in) : Elem(:)real (rkind),intent(in) : freqreal (rkind),intent(in) : det(:)integer (ikind) : ie,NElem,ireal (rkind) : F0(6),EK(6,6)GLoad(:)=ZeroNElem=size(Elem,1)do ie=1,NElemF0=Zerodo i=1,6if(Elem(ie)%GlbDOF(i)0) F0(i)=det(Elem(ie)%GlbDOF(i)end docall DEStifMat(EK,Elem,ie,freq)where (Elem(ie)%GlbDOF0)GLoad(Elem(ie)%GlbDOF)=GLoad(Elem(ie)%GlbDOF)+matmul(EK,F0)end where end doreturnend subroutine EGLoadVec!单元动力刚度矩阵!-subroutine EStifMat(EK,Elem,ie,freq)!-real(rkind),intent(out) :EK(6,6)type(typ_Element),intent(in) :Elem(:)integer(ikind),intent(in) :iereal(rkind),intent(in) :freqreal(rkind) :EI,EA,Length,mreal(rkind) :critreal(rkind) :B1,B2,T,R,Q,H,S,Creal(rkind) :nu,lambdainteger(ikind) :i,jEK = Zero EI = Elem(ie)%EIEA = Elem(ie)%EALength= Elem(ie)%Lengthm = Elem(ie)%massnu = freq*Length*sqrt(m/EA)lambda= Length*(freq*2*m/EI)*0.25crit = One-cosh(lambda)*cos(lambda)B1=nu/tan(nu)B2=nu/sin(nu)T = lambda*3*(sin(lambda)*cosh(lambda)+cos(lambda)*sinh(lambda)/critR = lambda*3*(sin(lambda)+sinh(lambda)/critQ = lambda*2*(sin(lambda)*sinh(lambda)/critH = lambda*2*(cosh(lambda)-cos(lambda)/critS = lambda*(sin(lambda)*cosh(lambda)-cos(lambda)*sinh(lambda)/critC = lambda*(sinh(lambda)-sin(lambda)/critEK(1,1) = B1*EA/LengthEK(1,4) = -B2*EA/LengthEK(2,2) = T*EI/Length*3EK(2,3) = Q*EI/Length*2EK(2,5) = -R*EI/Length*3EK(2,6) = H*EI/Length*2EK(3,3) = S*EI/LengthEK(3,5) = -H*EI/Length*2EK(3,6) = C*EI/LengthEK(4,4) = B1*EA/LengthEk(5,5) = T*EI/Length*3EK(5,6) = -Q*EI/Length*2EK(6,6) = S*EI/Lengthdo j=1,6do i=j+1,6EK(i,j)=EK(j,i)en