faintmm 发表于 2006-9-25 12:45

求助!WILSON法!急!

请教各位高手
我做的是复合板问题,
我编FORTRAN程序的步骤应该是怎样的呢?主要是输入数据的做法!请高手赐教!
还有在总装后进行微分方程组求解的时候,我找不到WISON法的FORTRAN程序,请有这个程序的高手给我一个!
llcbjjedut@163.com

boyi0706 发表于 2006-9-25 19:35

我这没fortran版的，只有matlab版的

欧阳中华 发表于 2006-9-25 20:56

.
找段十年前做的程序你研究一下吧，掌握问题还是要下功夫的...

      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
CINPUT VARIABLES
C         A => Stiffness matrix stored in compacted form
C         B => Mass matrix stored in compacted form
C         C => Damping matrix stored in compacted form
C   MAXA => Vector dontaining addresses of diagonal elements
C                   of stiffness matrix inA
C      U => 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 => Working storage
C         W => Working storage
C         U => Displancment , velocity and acceleration in t
C       UT => Displancment , velocity and acceleration in t+dt
C            R => Nodel load in t
C          RT => 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
Cin core , using compacted storage and column reduction scheme .
C
CINPUT VARIABLES
C         A => Stiffness matrix stored in compacted form
C          V => Right-hand-side load vecter
C      MAXA => Vector dontaining addresses of diagonal elements
C                   of stiffness matrix inA
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 => D and L - factors of stiffness matrix
C          V => 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

21172485 发表于 2006-9-25 22:10

这个程序不难，关于有限元程序设计的书里
我见过那里有FORTRAN的源程序，你找找看

faintmm 发表于 2006-9-26 13:48

非常感谢

在输入数据时要注意什么吗?

欧阳中华 发表于 2006-9-27 17:48

.
    将输入数据打入相应的数据文件就可以了... ...

faintmm 发表于 2006-10-8 18:50

谢谢!

欧阳老师,请问您有关于这个方法原理的资料吗?我找到一些都很粗略!
还有您程序中的参数的意义可以告诉我吗?

欧阳中华 发表于 2006-10-23 13:40

.
    《结构计算程序设计》张汝清、董明编，重庆出版社，1988年版P371，7.3直接积分法... ...

navigator 发表于 2006-10-24 00:27

对着书看看应该还好理解

faintmm 发表于 2006-10-24 21:02

谢谢!
我这个人性子很急呵呵!

hanxiao 发表于 2006-10-30 10:52

可以到源代码专区去找找看

faintmm 发表于 2006-11-2 21:41

IF(KK-2)40,150,150
这个语句是什么意思呢?
很多老书上都有!一直不明白!

欧阳中华 发表于 2006-11-3 13:11

.
    若KK-2小于0，转向40语句，等于0或大于0，转向150语句。

    一般FORTRAN书上都应该说明的.

faintmm 发表于 2006-11-3 15:11

噢!
可以这样用啊!
谢谢

faintmm 发表于 2006-11-11 15:07

请问欧阳老师
这个一维压缩存储是个什么概念呢?
在什么书上可以找到!?
页: [1] 2
查看完整版本: 求助!WILSON法!急!