newmark.for

工程中微分方程数值积分经常需要使用的Newmark法

－－－－－－－－－－－－－－－－－－－－－－－－－  
  
Module New_mark  
  
Contains  
  
Subroutine Newmark (N, dt, Nt, K, M, C, R, U, V, A)  
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
! 输入:  
! N 整体矩阵维数  
! dt 时程积分的步长  
! Nt 时程积分的总步数  
! K 整体刚度矩阵(N,N)  
! M 整体质量矩阵(N,N)  
! C 整体阻尼矩阵(N,N)  
! R 荷载矩阵(N,Nt)  
! 输出:  
! U 输出位移矩阵(N,Nt)  
! V 输出速度矩阵(N,Nt)  
! A 输出加速度矩阵(N,Nt)  
!  
!  Zhaoxin，2000,10.9  
!      [e-mail]  
! 参考文献：李杰，李国强，地震工程学导论，地震出版社  
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
USE NUMERICAL_LIBRARIES  
Use dflib  
Integer  N, Nt, LDA,IPVT(N),p,q  
Real(8)  dt,Gama,Alpha,Bete,A0,A1,A2,A3,A4,A5,A6,A7,aa  
Real(8), Dimension(N, N) :: K,M,C,TK,Ke,MINV  
Real(8), Dimension(N, Nt) :: R,U,V,A,RR  
TYPE (qwinfo)  winfo  
LOGICAL(4)     result  
  
OPEN (UNIT=1, FILE='USER',   
     *      TITLE='TRANSIENT ANALYSIS OF NEWMARK METHOD')  
winfo%TYPE = QWIN$SET  
winfo%X = 80  
winfo%Y = 0  
winfo%H = 34  
winfo%W = 60  
result=SETWSIZEQQ(1, winfo)  
!  CALL DLINRG (N, M, N, MINV, N)  
  
Gama=0.005  
Alpha=0.5+Gama  
Bete=0.25*(0.5+Alpha)**2  
A0=1/(Bete*dt**2)  
A1=Alpha/(Bete*dt)  
A2=1/(Bete*dt)  
A3=1/(2*Bete)-1  
A4=Alpha/Bete-1  
A5=dt/2*(Alpha/Bete-2)  
A6=dt*(1-Alpha)  
A7=Alpha*dt  
  
Ke=K+A0*M+A1*C  
  
do i=1,N  
if (Ke(i,i).eq.0) then  
  
aa=i  
aa=i  
end if  
  
end do  
  
LDA=N  
CALL DLFTSF (N, Ke, LDA, TK, LDA, IPVT)  
  
Do i=2, Nt  
  
WRITE (1,'("***********   LOAD STEP: ",I5,"   ***********")') i  
  
RR(1:N,i)= R(1:N,i)+Matmul(M,A0*U(1:N,i-1)+A2*V(1:N,i-1)  
     * +A3*A(1:N,i-1))+ Matmul(C,A1*U(1:N,i-1)  
     * +A4*V(1:N,i-1)+A5*A(1:N,i-1))  
  
CALL DLFSSF (N, TK, LDA, IPVT, RR(1:N,i), U(1:N,i))  
  
A(1:N,i)=A0*(U(1:N,i)-U(1:N,i-1))-A2*V(1:N,i-1)-A3*A(1:N,i-1)  
V(1:N,i)= V(1:N,i-1)+A6*A(1:N,i-1)+A7*A(1:N,i)  
  
End do  
  
WRITE (1,'("***********   SOLUTION IS DONE!  ***********")')  
  
End Subroutine Newmark  
  
End Module