请问谁有求响应的程序呀？

erwin111

请问谁有求响应的程序呀？

xawbw

问题没有讲清楚！

erwin111

嘻嘻，简谐荷载的响应，作用在一般的梁或者框架上的

MVH

这个应该不难，不知道是否已经由数学模型了？

erwin111

数学模型？还没有。我只是想参考一下，想看看求响应具体用程序是怎么实现的。

MVH

原帖由 erwin111 于 2006-9-5 20:29 发表
数学模型？还没有。我只是想参考一下，想看看求响应具体用程序是怎么实现的。

还实现好好研究一下数学模型是怎么样的吧，模型搞清楚了，计算方面的问题就没那么困难了

欧阳中华

.

      响应计算也是包含时域和频域，时域里比较普通的是Wilsion-0法，频域里就是求传递函数喽。

     列段Wilsion-0法的程序：

        IMPLICIT REAL*8 (A-H,O-Z)
        COMMON /FILE/  IIN,IOUT,IELE
        DIMENSION GM(3),GK(3),GC(3),MAXA(3)
        DIMENSION U(3,2),UT(3,2),R(2),RT(2),D(3),W(3)
        IIN = 1
        IOUT= 2
        OPEN(IIN,FILE='WILSON.INP')
        OPEN(IOUT,FILE='WILSON.OUT')
        MAXA(1)= 1
        MAXA(2)= 2
        MAXA(3)= 4
        GM(1) = 2.
        GM(2) = 1.
        GM(3) = 0.
        GK(1) = 6.
        GK(2) = 4.
        GK(3) =-2.
        R(1)  = 0.
        R(2)  = 0.
        U(1,1)= 0.
        U(1,2)= 0.
        U(2,1)= 0.
        U(2,2)= 0.
        U(3,1)= 0.
        U(3,2)= 10.
        NWK = 3
        NWM = 3
        NN  = 2
        READ(IIN,'(I5)') KNW
        CALL WILSON(GK,GM,MAXA,GC,U,UT,R,RT,D,W,NWK,NWM,NN,KNW)
        STOP
        END

C-------------------------------------------------------------1998.06.10
C       This subroutine is used to calculate responce ,
C            using the Wilson or Newmark direct integral methods .
C
C  INPUT VARIABLES
C         A[NWK] => Stiffness matrix stored in compacted form
C         B[NWM] => Mass matrix stored in compacted form
C         C[NWK] => Damping matrix stored in compacted form
C     MAXA[NN+1] => Vector dontaining addresses of diagonal elements
C                   of stiffness matrix in  A
C        U[3,NN] => Inital displancment , velocity and acceleration
C             DT => Integral step of time
C          NTIME => Number of time integration
C             NN => Number of equation
C            NNM => NN + 1
C            NWK => Number of elements below skyline of stiffness matrix
C            NWM => Number of elements below skyline of mass matrix
C            IIN => Number of input device
C           IOUT => Number of output device
C
C   OUTPUT VARIABLES
C          D[NN] => Working storage
C          W[NN] => Working storage
C        U[3,NN] => Displancment , velocity and acceleration in t
C       UT[3,NN] => Displancment , velocity and acceleration in t+dt
C          R[NN] => Nodel load in t
C         RT[NN] => Nodel load in t+dt
C-----------------------------------------------------------------------
        SUBROUTINE WILSON(A,B,MAXA,C,U,UT,R,RT,D,W,NWK,NWM,NN,KNW)
        IMPLICIT REAL*8 (A-H,O-Z)
        COMMON /FILE/  IIN,IOUT,IELE
        DIMENSION A(1),B(1),C(1),MAXA(1)
        DIMENSION U(3,1),UT(3,1),R(1),RT(1),D(1),W(1)
C
        READ(IIN,'(I5,3F10.4)') NTIME,DT,ALPHA,BETA
        IF(NWM.EQ.NN) THEN
        DO 50 I=1,NN
        KL = MAXA(I)
   50   C(KL) = ALPHA*A(KL) + BETA*B(I)
        ELSE
        DO 60 I=1,NWK
   60   C(I) = ALPHA*A(I) + BETA*B(I)
        ENDIF
C
        IF(KNW.EQ.1) THEN
        STA= 1.4
        A0 = 6./(STA*STA*DT*DT)
        A1 = 3./(STA*DT)
        A2 = 2.*A1
        F1 = 2.
        F2 = 2.
        F3 = STA*DT/2.
        F4 = A0/STA
        F5 =-A2/STA
        F6 = 1.-3./STA
        F7 = DT/2.
        F8 = F7
        A8 = DT*DT/6.
        ELSE
        DLT= 0.50
        ALF= 0.25
        STA= 1.00
        A0 = 1./(ALF*DT*DT)
        A1 = DLT/(ALF*DT)
        A2 = 1./(ALF*DT)
        F1 = 1./(2.*ALF) - 1.
        F2 = DLT/ALF - 1.
        F3 = DT*(DLT/ALF-2.)/2.
        F4 = A0
        F5 =-A2
        F6 =-F1
        F7 = DT*(1.-DLT)
        F8 = DLT*DT
        ENDIF
C
        DO 80 I=1,NWK
   80   A(I) = A(I) + A1*C(I)
        IF(NWM.EQ.NN) THEN
        DO 90 I=1,NN
        KL = MAXA(I)
   90   A(KL) = A(KL) + A0*B(I)
        ELSE
        DO 100 I=1,NWK
  100   A(I) = A(I) + A0*B(I)
        ENDIF
        CALL COLSOL(A,R,MAXA,NN,1)
C
        DO 200 I=1,NTIME
        READ(IIN,'(2E10.4)') (RT(J),J=1,NN)
        DO 130 J=1,NN
  130   D(J) = A0*U(1,J) + A2*U(2,J) + F1*U(3,J)
        CALL MULT(W,B,D,MAXA,NN,NWM)
        DO 140 J=1,NN
  140   R(J) = R(J) + STA*(RT(J)-R(J)) + W(J)
        DO 150 J=1,NN
  150   D(J) = A1*U(1,J) + F2*U(2,J) + F3*U(3,J)
        CALL MULT(W,C,D,MAXA,NN,NWK)
        DO 160 J=1,NN
  160   R(J) = R(J) + W(J)
C
        CALL COLSOL(A,R,MAXA,NN,2)
        DO 170 J=1,NN
        UT(3,J) = F4*(R(J)-U(1,J)) + F5*U(2,J) + F6*U(3,J)
        UT(2,J) = U(2,J) + F7*U(3,J) + F8*UT(3,J)
        UT(1,J) = U(1,J) + DT*U(2,J) + A8*(UT(3,J) + 2.*U(3,J))
        IF(KNW.EQ.2) UT(1,J)=R(J)
  170   CONTINUE
        DO 190 J=1,NN
        DO 180 K=1,3
  180   U(K,J)= UT(K,J)
  190   R(J)  = RT(J)
        WRITE(IOUT,'(I5,2F10.4)') I,UT(1,1),UT(1,2)
  200   CONTINUE
        RETURN
        END
C-----------------------------------------------------------------------
C       This subroutine is used to solve finite element static equations
C  in core , using compacted storage and column reduction scheme .
C
C  INPUT VARIABLES
C         A[NWK] => Stiffness matrix stored in compacted form
C          V[NN] => Right-hand-side load vecter
C      MAXA[NNM] => Vector dontaining addresses of diagonal elements
C                   of stiffness matrix in  A
C            NN  => Number of equations in the system
C           NWK  => Number of elements below skyline of matrix
C           NNM  => NN + 1
C            KK  => EQ.1 triangularization of matrix
C                   EQ.2 reduction and back-substitution of load vector
C           IOUT => Number of output device
C
C   OUTPUT VARIABLES
C         A[NWK] => D and L - factors of stiffness matrix
C          V[NN] => Displacement vector
C
C-----------------------------------------------------------------------
        SUBROUTINE COLSOL(A,V,MAXA,NN,KK)
        IMPLICIT REAL*8 (A-H,O-Z)
        COMMON /FILE/  IIN,IOUT,IELE
        DIMENSION A(1),V(1),MAXA(1)
C
C       PERFORM L*D*L(T) FACTORIZATION OF STIFFNESS MATRIX
C
        IF(KK-2)40,150,150
   40   DO 140 N=1,NN
        KN = MAXA(N)
        KL = KN + 1
        KU = MAXA(N+1) - 1
        KH = KU - KL
        IF(KH) 110,90,50
   50   K = N - KH
        IO = 0
        KLT = KU
        DO 80 J=1,KH
        IO = IO + 1
        KLT = KLT - 1
        KI = MAXA(K)
        ND =MAXA(K+1) - KI - 1
        IF(ND) 80,80,60
   60   KK = MIN(IO,ND)
        C = 0.
        DO 70 L=1,KK
   70   C = C + A(KI+L)*A(KLT+L)
        A(KLT)=A(KLT) - C
   80   K = K + 1
   90   K = N
        B = 0.
        DO 100 KK=KL,KU
        K = K - 1
        KI = MAXA(K)
        C = A(KK)/A(KI)
        B = B + C*A(KK)
  100   A(KK) = C
        A(KN) = A(KN) - B
  110   IF(A(KN)) 120,120,140
  120   WRITE(IOUT,2000)N,A(KN)
        STOP
  140   CONTINUE
        RETURN
C
C       REDUCE RIGHT-HAND-SIDE LOAD VECTOR
C
  150   DO 180 N=1,NN
        KL = MAXA(N) + 1
        KU = MAXA(N+1) - 1
        IF(KU-KL) 180,160,160
  160   K = N
        C = 0.
        DO 170 KK=KL,KU
        K = K - 1
  170   C = C + A(KK)*V(K)
        V(N) = V(N) - C
  180   CONTINUE
C
C       BACK-SUBSTITUTE
C
        DO 200 N=1,NN
        K = MAXA(N)
  200   V(N) = V(N)/A(K)
        IF(NN.EQ.1) RETURN
        N = NN
        DO 230 L=2,NN
        KL = MAXA(N) + 1
        KU = MAXA(N+1) - 1
        IF(KU-KL) 230,210,210
  210   K = N
        DO 220 KK=KL,KU
        K = K - 1
  220   V(K) = V(K) - A(KK)*V(N)
  230   N = N - 1
        RETURN
C
2000   FORMAT(//5x,'STOP !  Stiffness matrix not positive definite'//
     &  15x,'Nonpositive pivot for equation = ',i4//
     &  15x,'                         Pivot = ',E20.12//)
        END
C-----------------------------------------------------------------------
C   PROGRAM TO EVALUATE PRODUCT OF B TIMES RR AND STORE RESULT IN TT
C-----------------------------------------------------------------------
        SUBROUTINE MULT(TT,B,RR,MAXA,NN,NWM)
        IMPLICIT REAL*8 (A-H,O-Z)
        DIMENSION TT(1),B(1),RR(1),MAXA(1)
        IF(NWM.GT.NN) GOTO 20
        DO 10 I=1,NN
   10   TT(I) = B(I)*RR(I)
        RETURN
C
   20   DO 40 I=1,NN
   40   TT(I) = 0.
        DO 100 I=1,NN
        KL = MAXA(I)
        KU = MAXA(I+1) - 1
        II = I + 1
        CC = RR(I)
        DO 100 KK=KL,KU
        II = II - 1
  100   TT(II) = TT(II) + B(KK)*CC
        IF(NN.EQ.1) RETURN
        DO 200 I=2,NN
        KL = MAXA(I) + 1
        KU = MAXA(I+1) - 1
        IF(KU-KL) 200,210,210
  210   II = I
        AA = 0.
        DO 220 KK=KL,KU
        II = II - 1
  220   AA = AA + B(KK)*RR(II)
        TT(I) = TT(I) + AA
  200   CONTINUE
        RETURN
        END
C-----------------------------------------------------------------------
C     Ntot , Ncom =  285   88

输入文件：

1
12,0.28,0.,0.
0.,10.
0.,10.
0.,10.
0.,10.
0.,10.
0.,10.
0.,10.
0.,10.
0.,10.
0.,10.
0.,10.
0.,10.


输出文件：
    1     .0069     .4002
    2     .0598    1.5422
    3     .2242    3.0271
    4     .5588    4.4110
    5    1.0783    5.3466
    6    1.7292    5.6472
    7    2.3902    5.3203
    8    2.8993    4.5457
    9    3.1002    3.6024
   10    2.8932    2.7722
   11    2.2729    2.2505
   12    1.3411    2.0942


     好久以前的东西了，有问题自己研究吧... ...

NASA

龙格库塔法、newmark法、打靶法也都是很常见的方法吧

zfhnlg

NASA 发表于 2006-9-7 21:07

                               
登录/注册后可看大图

龙格库塔法、newmark法、打靶法也都是很常见的方法吧

对拟合点打耙法，求周期解，您了解吗？

sszy0336