讨论有限元法求解瞬态响应的Wilson-θ程序
本人正在求解转子系统的瞬态响应问题，需要用FORTRAN编写的用于有限元法求解瞬态响应的Wilson-θ程序，C也可以。[ 本帖最后由 gh688 于 2008-5-25 21:47 编辑 ] .
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 .
贡献段程序代码，还得花费9个振动币？...
回复 3楼 的帖子
你的振动币确实太少了！积分和威望很高:lol
致谢
非常感谢！以后一定把这种共享精神发扬下去 给你增加振动币，呵呵
页:
[1]