2023年程序结构力学大作业

程序构造力学大作业结81孙玉进该程序可计算任意平面构造任意阶频率（包括重频），以及任意阶振型（包括重频对应正交振型）。计算时，振型不包括不包括单元固端型振型。‍

!****************************

moduleNumKind

!***************

implicitnone

integer(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._rkind

endmoduleNumKind

!**************

moduleTypeDef

!**************

useNumKind

implicitnonetype::typ_Joint

real(rkind)::x,y

integer(ikind)::GDOF(3)

endtypetyp_Jointtype::typ_Element

integer(ikind)::JointNo(2),GlbDOF(6)

real(rkind)::Length,CosA,SinA,EI,EA,mass

endtypetyp_Elementtype::typ_JointLoad

integer(ikind)::JointNo,LodDOF

real(rkind)::LodVal

endtypetyp_JointLoadtype::typ_ElemLoad

integer(ikind)::ElemNo,Indx

real(rkind)::Pos,LodVal

endtypetyp_ElemLoadcontains!===================================

subroutineSetElemProp(Elem,Joint)

!===================================

type(typ_Element),intent(inout)::Elem(:)

type(typ_Joint),intent(in)::Joint(:)

integer(ikind)::ie,NElem

real(rkind)::x1,x2,y1,y2

NElem=size(Elem,dim=1)

doie=1,NElem

x1=Joint(Elem(ie)%JointNo(1))%x

y1=Joint(Elem(ie)%JointNo(1))%y

x2=Joint(Elem(ie)%JointNo(2))%x

y2=Joint(Elem(ie)%JointNo(2))%y

Elem(ie)%Length=sqrt((x1-x2)**2+(y1-y2)**2)

Elem(ie)%CosA=(x2-x1)/Elem(ie)%Length

Elem(ie)%SinA=(y2-y1)/Elem(ie)%Length

Elem(ie)%GlbDOF(1:3)=Joint(Elem(ie)%JointNo(1))%GDOF(:)

Elem(ie)%GlbDOF(4:6)=Joint(Elem(ie)%JointNo(2))%GDOF(:)

enddo

return

endsubroutineSetElemProp!======================================

subroutineTransMatrix(ET,CosA,SinA)

!======================================

real(rkind),intent(out)::ET(:,:)

real(rkind),intent(in)::CosA,SinA

!ETcouldbe2x2,3x3or6x6dependingsize(ET)

ET=Zero

ET(1,1:2)=(/CosA,SinA/)

ET(2,1:2)=(/-SinA,CosA/)

if(size(ET,1)>2)ET(3,3)=One

if(size(ET,1)>3)ET(4:6,4:6)=ET(1:3,1:3)

return

endsubroutineTransMatrixendmoduleTypeDef

!*************

moduleBandMat

!*************

useNumKind

useTypeDef,only:typ_Element!仅用该模块中旳typ_Element

implicitnoneprivate!默认所有旳数据和过程为私有，增强封装性

public::SetMatBand,DelMatBand,VarBandSolvtype,public::typ_Kcol

real(rkind),pointer::row(:)

endtypetyp_Kcolcontains!==================================

subroutineSetMatBand(Kcol,Elem)

!==================================

!...[6-4-2]type(typ_KCol),intent(inout)::Kcol(:)

type(typ_Element),intent(in)::Elem(:)

integer(ikind)::minDOF,ELocVec(6)

integer(ikind)::ie,j,NGlbDOF,NElem

integer(ikind)::row1(size(Kcol,dim=1))

!row1是自动数组，子程序结束后将自动释放内存空间

NGlbDOF=size(Kcol,1)

NElem=size(Elem,1)

row1=NGlbDOF!先设始行码为大数

!确定各列始行码，放在数组row1(:)中

doie=1,NElem

ELocVec=Elem(ie)%GlbDOF

minDOF=minval(ELocVec,mask=ELocVec>0)

where(ELocVec>0)

row1(ELocVec)=min(row1(ELocVec),minDOF)

endwhere

enddo

!为各列旳半带宽分派空间并初始化

doj=1,NGlbDOF

allocate(Kcol(j)%row(row1(j):j))

Kcol(j)%row=Zero!清零

enddo

return

endsubroutineSetMatBand!===================================

subroutineDelMatBand(Kcol)

!===================================

!...[6-5-5]

type(typ_KCol),intent(inout)::Kcol(:)

integer(ikind)::j,NGlbDOF

NGlbDOF=size(Kcol,1)

doj=1,NGlbDOF

deallocate(Kcol(j)%row)

enddo

return

endsubroutineDelMatBand!=================================================

subroutineVarBandSolv(Disp,Kcol,GLoad)

!=================================================

type(typ_KCol),intent(inout)::Kcol(:)

real(rkind),intent(out)::Disp(:)

real(rkind),intent(in)::GLoad(:)

integer(ikind)::i,j,k,row1j,row_1,NCol

real(rkind)::Diag(size(Kcol,dim=1))

real(rkind)::s

!...[6-5-2]

NCol=size(Kcol,1)

Diag(1:NCol)=(/(Kcol(j)%row(j),j=1,NCol)/)

doj=2,NCol

row1j=lbound(Kcol(j)%row,1)

doi=row1j,j-1

row_1=max(row1j,lbound(Kcol(i)%row,1))

k=i-1

s=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)

enddo

s=sum(Diag(row1j:j-1)*Kcol(j)%row(row1j:j-1)**2)

Diag(j)=Diag(j)-s

enddo

Disp(:)=GLoad(:)

!...[6-5-3节旳代码：其中GP换为Disp]

doj=2,NCol

row1j=lbound(Kcol(j)%row,1)

Disp(j)=Disp(j)-sum(Kcol(j)%row(row1j:j-1)*Disp(row1j:j-1))

enddo

!...[6-5-4节旳代码：其中GP换为Disp]

Disp(:)=Disp(:)/Diag(:)

doj=NCol,1,-1

row1j=lbound(Kcol(j)%row,1)

Disp(row1j:j-1)=Disp(row1j:j-1)-Disp(j)*Kcol(j)%row(row1j:j-1)

enddo

return

endsubroutineVarBandSolv

endmoduleBandMat!****************************

moduleDispMethod

!****************************

useNumKind

useTypeDef

useBandMat

implicitnonecontains

!========================

subroutineSolveDisp(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)::NGlbDOF

NGlbDOF=size(Disp,dim=1)

allocate(GLoad(NGlbDOF))

allocate(Kcol(NGlbDOF))callSetMatBand(Kcol,Elem)

callGLoadVec(GLoad,Elem,JLoad,ELoad,Joint)

callGStifMat(Kcol,Elem)

callVarBandSolv(Disp,Kcol,GLoad)

callDelMatBand(Kcol)

deallocate(GLoad)

return

endsubroutineSolveDisp

!========================

subroutineGStifMat(Kcol,Elem)

!========================

type(typ_Kcol),intent(inout)::Kcol(:)

type(typ_Element),intent(in)::Elem(:)

!...[6-4-3]

integer(ikind)::ie,j,JGDOF,NElem

real(rkind)::EK(6,6),ET(6,6)

integer(ikind)::ELocVec(6)

NElem=size(Elem,1)

doie=1,NElem

callEStifMat(EK,Elem(ie)%Length,Elem(ie)%EI,Elem(ie)%EA)

callTransMatrix(ET,Elem(ie)%CosA,Elem(ie)%SinA)

EK=matmul(transpose(ET),matmul(EK,ET))

ELocVec=Elem(ie)%GlbDOF

doj=1,6

JGDOF=ELocVec(j)

if(JGDOF==0)cycle

where(ELocVec<=JGDOF.and.ELocVec>0)

Kcol(JGDOF)%row(ELocVec)=Kcol(JGDOF)%row(ELocVec)+EK(:,j)

endwhere

enddo

enddo

return

endsubroutineGStifMat

!========================

subroutineGLoadVec(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,NELoad

real(rkind)::F0(6),ET(6,6)

NJLoad=size(JLoad,dim=1)

NELoad=size(ELoad,dim=1)

GLoad(:)=Zero

doie=1,NJLoad

GLoad(Joint(JLoad(ie)%JointNo)%GDOF(JLoad(ie)%LodDOF))=GLoad(Joint(JLoad(ie)%JointNo)%GDOF(JLoad(ie)%LodDOF))+JLoad(ie)%LodVal

enddo

doie=1,NELoad

callEFixendF(F0,ELoad(ie)%Indx,ELoad(ie)%Pos,ELoad(ie)%LodVal,Elem(ELoad(ie)%ElemNo))

callTransMatrix(ET,Elem(ELoad(ie)%ElemNo)%CosA,Elem(ELoad(ie)%ElemNo)%SinA)

where(Elem(ELoad(ie)%ElemNo)%GlbDOF>0)

GLoad(Elem(ELoad(ie)%ElemNo)%GlbDOF)=GLoad(Elem(ELoad(ie)%ElemNo)%GlbDOF)-matmul(transpose(ET),F0(:))

endwhere

enddo

return

endsubroutineGLoadVec!========================

subroutineEStifMat(EK,Elen,EI,EA)

!========================

real(rkind),intent(out)::EK(:,:)

real(rkind),intent(in)::Elen,EI,EA

integer(ikind)::i,j

EK=Zero

EK(1,1)=EA/Elen

EK(1,4)=-EA/Elen

EK(2,2)=Twelve*EI/Elen**Three

EK(2,3)=Six*EI/Elen**Two

EK(2,5)=-Twelve*EI/Elen**Three

EK(2,6)=Six*EI/Elen**Two

EK(3,3)=Four*EI/Elen

EK(3,5)=-Six*EI/Elen**Two

EK(3,6)=Two*EI/Elen

EK(4,4)=EA/Elen

Ek(5,5)=Twelve*EI/Elen**Three

EK(5,6)=-Six*EI/Elen**Two

EK(6,6)=Four*EI/Elen

doj=1,6

doi=j+1,6

EK(i,j)=EK(j,i)

enddo

enddo

return

endsubroutineEStifMat!========================

subroutineEFixendF(F0,Indx,a,q,Elem)

!========================

real(rkind),intent(out)::F0(:)

real(rkind),intent(in)::a,q

integer(ikind),intent(in)::Indx

type(typ_Element),intent(in)::Elem

real::l

l=Elem%length

selectcase(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(a<One/Two)then

F0(1)=Elem%EA*q/l

F0(4)=-F0(1)

else

F0(1)=-Elem%EA*q/l

F0(4)=-F0(1)

endif

F0(2)=Zero

F0(3)=Zero

F0(5)=Zero

F0(6)=Zero

case(4)

if(a<One/Two)then

F0(2)=Twelve*Elem%EI*q/l**Three

F0(3)=Six*Elem%EI*q/l**Two

F0(5)=-F0(2)

F0(6)=F0(3)

else

F0(2)=-Twelve*Elem%EI*q/l**Three

F0(3)=-Six*Elem%EI*q/l**Two

F0(5)=-F0(2)

F0(6)=F0(3)

endif

F0(1)=Zero

F0(4)=Zero

endselect

return

endsubroutineEFixendF

!========================

subroutineElemDisp(EDisp,Disp,Elem)

!========================

real(rkind),intent(out)::EDisp(:)

real(rkind),intent(in)::Disp(:)

type(typ_Element),intent(in)::Elem

integer(ikind)::i

doi=1,6

if(Elem%GlbDOF(i)==0)then

EDisp(i)=Zero

else

EDisp(i)=Disp(Elem%GlbDOF(i))

endif

enddo

return

endsubroutineElemDisp

!========================

subroutineElemForce(EForce,ie,Disp,Elem,ELoad)

!========================

real(rkind),intent(out)::EForce(:)

real(rkind),intent(in)::Disp(:)

type(typ_Element),intent(in)::Elem(:)

type(typ_ElemLoad),intent(in)::ELoad(:)

integer(ikind),intent(in)::ie

real(rkind)::EK(6,6),ET(6,6),EP(6),F0(6),EDisp(6)

integer(ikind)::NELoad,i

NELoad=size(ELoad,1)

EP(:)=0

doi=1,NELoad

if(ELoad(i)%ElemNo==ie)then

callEFixendF(F0,ELoad(i)%Indx,ELoad(i)%Pos,ELoad(i)%LodVal,Elem(ELoad(i)%ElemNo))

EP(:)=-F0(:)

endif

enddo

callEStifMat(EK,Elem(ie)%Length,Elem(ie)%EI,Elem(ie)%EA)

callTransMatrix(ET,Elem(ie)%CosA,Elem(ie)%SinA)

callElemDisp(EDisp,Disp,Elem(ie))

EForce(:)=matmul(EK,matmul(ET,EDisp))-EP(:)

return

endsubroutineElemForce

endmoduleDispMethod!************************

moduleVibration

!**************************

useNumKind

useTypeDef

useBandMat

implicitnone

real(rkind)::Pi=3.contains!-----------------------------------------------------

subroutineGetFreq(Freq,Elem,Kcol,FreqStart,NFreq,Tol)!二分法获取频率

!-----------------------------------------------------

real(rkind),intent(inout)::Freq(:)

type(typ_Element),intent(in)::Elem(:)

type(typ_Kcol),intent(inout)::Kcol(:)

integer(ikind),intent(in)::FreqStart,NFreq

real(rkind),intent(in)::Tol

real(rkind)::freq_l,freq_u,freqm

integer(ikind)::k

integer(ikind)::J01,J02,J_l,J_u,Jk,J0,Jm,MFreq

MFreq=1

freq_l=One

freq_u=Ten

dok=FreqStart,FreqStart+NFreq-1

if(MFreq>1)then

MFreq=MFreq-1

Freq(k-FreqStart+1)=Freq(k-FreqStart)

cycle

endif

if(k>FreqStart)freq_l=freq(k-FreqStart)

do

callCJ0(J01,freq_l,Elem)

callCJk(Jk,freq_l,Kcol,Elem)

J_l=J01+Jk

if(J_l<k)exit

freq_l=freq_l/Two

enddo

do

callCJ0(J02,freq_u,Elem)

callCJk(Jk,freq_u,Kcol,Elem)

J_u=J02+Jk

if(J_u>=k)exit

freq_l=freq_u

freq_u=freq_u*Two

enddo

do

freqm=(freq_l+freq_u)/Two

callCJ0(J0,freqm,Elem)

callCJk(Jk,freqm,Kcol,Elem)

Jm=J0+Jk

if(Jm>=k)then

freq_u=freqm

else

freq_l=freqm

endif

if((freq_u-freq_l)<Tol*(One+freq_l))then

callCJ0(J0,freq_u,Elem)

callCJk(Jk,freq_u,Kcol,Elem)

Jm=J0+Jk

MFreq=Jm-k+1

exit

endif

enddo

Freq(k+1-FreqStart)=(freq_u+freq_l)/2

enddo

return

endsubroutineGetFreq

!--------------------------------

subroutineMode(Elem,Kcol,FreqStart,NFreq,Tol)!获取振型

!--------------------------------

type(typ_Element),intent(in)::Elem(:)

type(typ_Kcol),intent(inout)::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,r

NGlbDOF=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=ZerocallGetFreq(Freq,Elem,Kcol,FreqStart,NFreq,Tol)

doi=1,NFreq-1!计算重频对应旳已出现次数

if((Freq(i+1)-Freq(i))<Two*Tol*(1+Freq(i+1)))MFreq(i+1)=MFreq(i)+1

enddoallocate(Mdet(NGlbDOF,sum(MFreq)))

write(8,*)NFreqdok=1,rcallrandom_number(det(:))!振型初始化

det0=-One

doif(MFreq(k)>0)then!重频振型正交化

dot=1,MFreq(k)

callEGLoadVec(GLoad,Elem,det,Freq(k))

callEGLoadVec(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)*arf

enddo

endif

!逆幂迭代

callSetMatBand(Kcol,Elem)

callGStifMat(Kcol,Elem,Freq(k))

diag(:)=(/(Kcol(j)%row(j),j=1,NCol)/)

doj=2,NCol

row1j=lbound(Kcol(j)%row,1)

doi=row1j,j-1

row_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)

enddo

s=sum(diag(row1j:j-1)*Kcol(j)%row(row1j:j-1)**2)

diag(j)=diag(j)-s

enddo

doj=2,NCol

row1j=lbound(Kcol(j)%row,1)

GLoad(j)=GLoad(j)-sum(Kcol(j)%row(row1j:j-1)*GLoad(row1j:j-1))

enddo

GLoad(:)=GLoad(:)/Diag(:)

doj=NCol,1,-1

row1j=lbound(Kcol(j)%row,1)

GLoad(row1j:j-1)=GLoad(row1j:j-1)-GLoad(j)*Kcol(j)%row(row1j:j-1)

enddo

de=GLoad(1)

doi=2,NGlbDOF!原则化

if(abs(GLoad(i))>abs(de))de=GLoad(i)

enddo

doi=1,NGlbDOF

GLoad(i)=GLoad(i)/de

enddo

if(dot_product(GLoad-det0,GLoad-det0)<Tol*Tol)then

if(MFreq(k+1)>0)Mdet(1:NGlbDOF,MFreq(k+1))=GLoad!记录重频振型，用于后续正交化

exit

endif

det0=GLoad

freq0=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)=freq0

enddo

write(8,*)Freq(k),0

dot=1,NElem

callElemDisp2(EDisp,GLoad,Elem(t))

write(8,*)EDisp

enddo

enddo

EDisp=0

dok=r+1,NFreq

write(8,*)Freq(k),0

dot=1,NElem

write(8,*)EDisp

enddo

enddo

deallocate(diag)

deallocate(det0)

deallocate(GLoad)

callDelMatBand(Kcol)

returnendsubroutineMode

!计算单元固端频率数

subroutineCJ0(J0,freq,Elem)

integer(ikind),intent(out)::J0

real(rkind),intent(in)::freq

type(typ_Element),intent(in)::Elem(:)

integer(ikind)::NElem,Ja,Jb,i,j

real(rkind)::nu,lambda,sg

NElem=size(Elem,dim=1)

J0=0

doi=1,NElem

nu=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)/2

J0=J0+Ja+Jb

enddo

return

endsubroutineCJ0!计算Jk

subroutineCJk(Jk,freq,Kcol,Elem)

integer(ikind),intent(out)::Jk

real(rkind),intent(in)::freq

type(typ_Kcol),intent(inout)::Kcol(:)

type(typ_Element),intent(in)::Elem(:)

integer(ikind)::NGlbDOF

real(rkind),allocatable::diag(:)

integer(ikind)::i,j,row1j,row_1,NCol

real(rkind)::s

NGlbDOF=size(Kcol,1)

allocate(diag(NGlbDOF))

callSetMatBand(Kcol,Elem)

callGStifMat(Kcol,Elem,freq)

NCol=size(Kcol,1)

diag(:)=(/(Kcol(j)%row(j),j=1,NCol)/)

doj=2,NCol

row1j=lbound(Kcol(j)%row,1)

doi=row1j,j-1

row_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)

enddo

s=sum(diag(row1j:j-1)*Kcol(j)%row(row1j:j-1)**2)

diag(j)=diag(j)-s

enddo

Jk=count(mask=diag<Zero,dim=1)

deallocate(diag)

return

endsubroutineCJk!整体动力刚度矩阵

subroutineGStifMat(Kcol,Elem,freq)

type(typ_Kcol),intent(inout)::Kcol(:)

type(typ_Element),intent(in)::Elem(:)

real(rkind),intent(in)::freq

integer(ikind)::ie,j,JGDOF,NElem

real(rkind)::EK(6,6),ET(6,6)

integer(ikind)::ELocVec(6)

NElem=size(Elem,1)

doie=1,NElem

callEStifMat(EK,Elem,ie,freq)

callTransMatrix(ET,Elem(ie)%CosA,Elem(ie)%SinA)

EK=matmul(transpose(ET),matmul(EK,ET))

ELocVec=Elem(ie)%GlbDOF

doj=1,6

JGDOF=ELocVec(j)

if(JGDOF==0)cycle

where(ELocVec<=JGDOF.and.ELocVec>0)

Kcol(JGDOF)%row(ELocVec)=Kcol(JGDOF)%row(ELocVec)+EK(:,j)

endwhere

enddo

enddo

return

endsubroutineGStifMat!等效整体荷载向量，用于逆幂迭代

!========================

subroutineEGLoadVec(GLoad,Elem,det,freq)

!========================

real(rkind),intent(out)::GLoad(:)

type(typ_Element),intent(in)::Elem(:)

real(rkind),intent(in)::freq

real(rkind),intent(in)::det(:)

integer(ikind)::ie,NElem,i

real(rkind)::F0(6),EK(6,6)

GLoad(:)=Zero

NElem=size(Elem,1)

doie=1,NElem

F0=Zero

doi=1,6

if(Elem(ie)%GlbDOF(i)>0)F0(i)=det(Elem(ie)%GlbDOF(i))

enddo

callDEStifMat(EK,Elem,ie,freq)

where(Elem(ie)%GlbDOF>0)

GLoad(Elem(ie)%GlbDOF)=GLoad(Elem(ie)%GlbDOF)+matmul(EK,F0)

endwhere

enddo

return

endsubroutineEGLoadVec

!单元动力刚度矩阵

!----------------------------------

subroutineEStifMat(EK,Elem,ie,freq)

!----------------------------------

real(rkind),intent(out)::EK(6,6)

type(typ_Element),intent(in)::Elem(:)

integer(ikind),intent(in)::ie

real(rkind),intent(in)::freq

real(rkind)::EI,EA,Length,m

real(rkind)::crit

real(rkind)::B1,B2,T,R,Q,H,S,C

real(rkind)::nu,lambda

integer(ikind)::i,j

EK=Zero

EI=Elem(ie)%EI

EA=Elem(ie)%EA

Length=Elem(ie)%Length

m=Elem(ie)%mass

nu=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))/crit

R=lambda**3*(sin(lambda)+sinh(lambda))/crit

Q=lambda**2*(sin(lambda)*sinh(lambda))/crit

H=lambda**2*(cosh(lambda)-cos(lambda))/crit

S=lambda*(sin(lambda)*cosh(lambda)-cos(lambda)*sinh(lambda))/crit

C=lambda*(sinh(lambda)-sin(lambda))/critEK(1,1)=B1*EA/Length

EK(1,4)=-B2*EA/Length

EK(2,2)=T*EI/Length**3

EK(2,3)=Q*EI/Length**2

EK(2,5)=-R*EI/Length**3

EK(2,6)=H*EI/Length**2

EK(3,3)=S*EI/Length

EK(3,5)=-H*EI/Length**2

EK(3,6)=C*EI/Length

EK(4,4)=B1*EA/Length

Ek(5,5)=T*EI/Length**3

EK(5,6)=-Q*EI/Length**2

EK(6,6)=S*EI/Lengthdoj=1,6

doi=j+1,6

EK(i,j)=EK(j,i)

enddo

enddo

return

endsubroutineEStifMat!单元动力刚度矩阵旳导数矩阵

!----------------------------------

subroutineDEStifMat(EK,Elem,ie,freq)

!----------------------------------

real(rkind),intent(out)::EK(6,6)

type(typ_Element),intent(in)::Elem(:)

integer(ikind),intent(in)::ie

real(rkind),intent(in)::freq

real(rkind)::EI,EA,Length,m

real(rkind)::crit

real(rkind)::T,R,Q,H,S,C,B11,B21,T1,R1,Q1,H1,S1,C1

real(rkind)::nu,lambda

integer(ikind)::i,j

EK=Zero

EI=Elem(ie)%EI

EA=Elem(ie)%EA

Length=Elem(ie)%Length

m=Elem(ie)%mass

nu=freq*Length*sqrt(m/EA)

lambda=Length*(freq**2*m/EI)**0.25crit=One-cosh(lambda)*cos(lambda)

B11=nu/freq/tan(nu)-nu*nu/freq*(One/tan(nu)/tan(nu)+1)

B21=nu/freq/sin(nu)-nu*nu/freq*cos(nu)/sin(nu)/sin(nu)

T=lambda**3*(sin(lambda)*cosh(lambda)+cos(lambda)*sinh(lambda))/crit

R=lambda**3*(sin(lambda)+sinh(lambda))/crit

Q=lambda**2*(sin(lambda)*sinh(lambda))/crit

H=lambda**2*(cosh(lambda)-cos(lambda))/crit

S=lambda*(sin(lambda)*cosh(lambda)-cos(lambda)*sinh(lambda))/crit

C=lambda*(sinh(lambda)-sin(lambda))/critT1=One/Two/freq*lambda*(Three*T/lambda-T*S/lambda+Two*lambda**3*cos(lambda)*cosh(lambda)/crit)

R1=One/Two/freq*lambda*(Three*R/lambda-R*S/lambda+lambda**3*(cos(lambda)+cosh(lambda))/crit)

Q1=One/Two/freq*lambda*(Two*Q/lambda-Q*S/lambda+T/lambda)

H1=One/Two/freq*lambda*(Two*H/lambda-H*S/lambda+R/lambda)

S1=One/Two/freq*lambda*(S/lambda-S*S/lambda+Two*Q/lambda)

C1=One/Two/freq*lambda*(C/lambda-C*S/lambda+H/lambda)EK(1,1)=B11*EA/Length

EK(1,4)=-B21*EA/Length

EK(2,2)=T1*EI/Length**3

EK(2,3)=Q1*EI/Length**2

EK(2,5)=-R1*EI/Length**3

EK(2,6)=H1*EI/Length**2

EK(3,3)=S1*EI/Length

EK(3,5)=-H1*EI/Length**2

EK(3,6)=C1*EI/Length

EK(4,4)=B11*EA/Length

Ek(5,5)=T1*EI/Length**3

EK(5,6)=-Q1*EI/Length**2

EK(6,6)=S1*EI/Lengthdoj=1,6

doi=j+1,6

EK(i,j)=EK(j,i)

enddo

enddo

return

endsubroutineDEStifMat!========================

subroutineElemDisp2(EDisp,Disp,Elem)

!========================

real(rkind),intent(out)::EDisp(:)

real(rkind),intent(in)::Disp(:)

type(typ_Element),intent(in)::Elem

integer(ikind)::i

doi=1,6

if(Elem%GlbDOF(i)==0)then

EDisp(i)=Zero

else

EDisp(i)=Disp(Elem%GlbDOF(i))

endif

enddo

return

endsubroutineElemDisp2endmoduleVibration

!*****************************

programSM_90

!*****************************

useNumKind

useTypeDef

useDispMethod

useVibrationimplicitnonetype(typ_Element),allocatable::Elem(:)

type(typ_Joint),allocatable::Joint(:)

type(typ_JointLoad),allocatable::JLoad(:)

type(typ_Elemload),allocatable::ELoad(:)

real(rkind),allocatable