abaya.txt

源码公开－线弹性动力分析直接积分 声明： 贴出源程序只是为了给有这方面需要的弟兄提供参考

发信人: ABAYA (ABAYA), 信区: FEA
标  题: 源码公开－线弹性动力分析直接积分
发信站: 同舟共济站 (Sat Sep 28 18:26:19 2002)

声明：
贴出源程序只是为了给有这方面需要的弟兄提供参考，并
非用于和做这方面工作的弟兄讨论用。
本人不对程序的正确性提供任何保证，不对程序的可读性
承担任何义务，也不承担对使用者在使用过程中的任何问
题提供解答的义务。
同意本声明的弟兄可以以任何的方式修改、使用
和传播本程序。否则，请勿以任何方式使用本程序。

NEWMARK BETA源程序（FORTRAN）。程序是一个独立的文件（
F77格式），使用者可将其拷入一个文件并改名为.FOR即可。
如：NEWMARK.FOR。程序中的矩阵和向量均为满存储。

    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)
!
!           赵昕，[e-mail]
!               同济大学建筑工程系，2000,11.9
!   参考文献：李杰，李国强，地震工程学导论，地震出版社
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    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