!************************************************************************* !* Response of a N d.o.f. Mass-Spring System with Damping * !* to a sinusoidal Force By a Direct Method * !* --------------------------------------------------------------------- * !* EXPLANATION: * !* Let us consider the 3 d.o.f. Mass-Spring system with viscous damping: * !* * !* --> x1 --> x2 k --> x3 * !* \| k *----* 2k *----*--/\/\/\--*----* * !* \|--/\/\/\--| 2m |--/\/\/\--| m | ---* | 3m | * !* \| *----* *----*---| |----*----* * !* --> F1(t) ---* * !* c * !* | 2 0 0 | * !* The Mass Matrix M is | 0 1 0 | assuming m=1. * !* | 0 0 3 | * !* * !* | 3 -2 0 | * !* The Stiffness Matrix K is | -2 3 -1 | assuming k=1. * !* | 0 -1 1 | * !* * !* | 0 0 0 | * !* The Damping Matrix is | 0 0 0 | assuming c=0.5. * !* | 0 0 0.5 | * !* * !* The resonance fresuencies (neglecting c) are: * !* F1 = w1/2pi = 0.05161 sqr(k/m) * !* F2 = w2/2pi = 0.1431 sqr(k/m) * !* F3 = w3/2pi = 0.3151 sqr(k/m) * !* with the corresponding modal vectors: * !* | 1 | | 1 | | 1 | * !* phi1 = | 1.395 | phi2 = | 0.6914 | phi3 = | -2.420 | * !* | 2.038 | |-0.4849 | | 0.2249 | * !* * !* (See program ndof01.bas). * !* --------------------------------------------------------------------- * !* SAMPLE RUN: * !* * !* How many degrees of freedom (d.o.f.): 3 * !* * !* Input Mass Matrix [M]: * !* M(1,1) = 2 * !* M(1,2) = 0 * !* M(1,3) = 0 * !* M(2,2) = 1 * !* M(2,3) = 0 * !* M(3,3) = 3 * !* * !* Viscous Damping...: VIS * !* Structural Damping: STR * !* * !* Your choice: VIS * !* * !* Input Stiffness Matrix [K]: * !* K(1,1) = 3 * !* K(1,2) = -2 * !* K(1,3) = 0 * !* K(2,2) = 3 * !* K(2,3) = -1 * !* K(3,3) = 1 * !* * !* Input Damping Matrix [C]: * !* C(1,1) = 0 * !* C(1,2) = 0 * !* C(1,3) = 0 * !* C(2,2) = 0 * !* C(2,3) = 0 * !* C(3,3) = 0.5 * !* * !* Input Excitation Vector: * !* F1(1) = 1000 * !* F1(2) = 0 * !* F1(3) = 0 * !* F2(1) = 0 * !* F2(2) = 0 * !* F2(3) = 0 * !* * !* Number of d.o.f. to calculate: 3 * !* * !* Frequency Sweep * !* --------------- * !* Starting Frequency: 0.31 * !* Ending Frequency..: 0.32 * !* * !* Linear frequency step.....: LIN * !* Logarithmic frequency step: LOG * !* * !* Your choice: LIN * !* * !* Frequency Step....: 0.0005 * !* * !* * !* Frequency (Hz) Displacement Phase (deg.) * !* -------------------------------------------- * !* 0.3100 240.66 0.0 * !* 0.3105 264.73 0.0 * !* 0.3110 294.72 0.0 * !* 0.3115 333.08 0.0 * !* 0.3120 383.85 0.0 * !* 0.3125 454.17 0.0 * !* 0.3130 557.90 0.0 * !* 0.3135 725.79 0.0 * !* 0.3140 1041.61 0.0 * !* 0.3145 1824.30 0.0 * !* 0.3150 4347.33 0.0 <-- 3rd resonance * !* 0.3155 2248.97 0.0 * !* 0.3160 1151.54 0.0 * !* 0.3165 757.68 0.0 * !* 0.3170 560.64 0.0 * !* 0.3175 443.08 0.0 * !* 0.3180 365.14 0.0 * !* 0.3185 309.74 0.0 * !* 0.3190 268.37 0.0 * !* 0.3195 236.30 0.0 * !* 0.3200 210.73 0.0 * !* * !* --------------------------------------------------------------------- * !* REFERENCE: "Mécanique des vibrations linéaires By M. Lalanne, * !* P. Berthier, J. Der Hagopian, Masson, Paris 1980". * !* [BIBLI 16]. * !* * !* F90 Release By J-P Moreau, Paris * !* (in simple precision). * !* (www.jpmoreau.fr) * !************************************************************************* PROGRAM NDOF03 INTEGER, PARAMETER :: NMAX = 12 REAL*4 M(NMAX, NMAX), K1(NMAX, NMAX), K2(NMAX, NMAX) REAL*4 A(2*NMAX, 2*NMAX), X(2*NMAX), B(2*NMAX), B1(2*NMAX, 1) INTEGER C9 REAL*4 I1, M1 CHARACTER*3 ANS PI = 4. * ATAN(1.) PRINT *,' ' WRITE(*,10,advance='no'); READ *, N N2 = 2 * N PRINT *,' ' PRINT *,' Input Mass Matrix [M]:' DO I = 1, N DO J = I, N WRITE(*,20,advance='no') I,J; READ *, M(I,J) END DO END DO PRINT *,' ' PRINT *,' Viscous Damping...: VIS' PRINT *,' Structural Damping: STR' PRINT *,' ' WRITE(*,30,advance='no'); READ *, ANS IF (ANS=='VIS') GOTO 600 !case structural damping C9 = 2 PRINT *,' ' PRINT *,' Input stiffness Matrix [K1]:' DO I = 1, N DO J = I, N WRITE(*,40,advance='no') I,J; READ *, K1(I,J) END DO END DO PRINT *,' ' PRINT *,' Input Damping Matrix [K2]:' DO I = 1, N DO J = I, N WRITE(*,41,advance='no') I,J; READ *, K2(I,J) END DO END DO PRINT *,' ' GOTO 900 600 CONTINUE !case viscous damping C9 = 1 PRINT *,' ' PRINT *,' Input stiffness Matrix [K]:' DO I = 1, N DO J = I, N WRITE(*,50,advance='no') I,J; READ *, K1(I,J) END DO END DO PRINT *,' ' PRINT *,' Input Damping Matrix [C]:' DO I = 1, N DO J = I, N WRITE(*,51,advance='no') I,J; READ *, K2(I,J) END DO END DO 900 CONTINUE !Complete Matrices by symmetry DO I = 1, N DO J = I, N K1(J, I) = K1(I, J) K2(J, I) = K2(I, J) M(J, I) = M(I, J) END DO END DO PRINT *,' ' PRINT *,' Input Excitation Vector:' DO I = 1, N WRITE(*,60,advance='no') I; READ *, B(I) END DO PRINT *,' ' DO I = N + 1, N2 J = I - N WRITE(*,61,advance='no') J; READ *, B(I) END DO PRINT *,' ' WRITE(*,70,advance='no'); READ *, N3 PRINT *,' ' !calculate numerical response of d.o.f. #N3 PRINT *,' ' PRINT *,' Frequency Sweep' PRINT *,' ---------------' PRINT *,' ' 1400 CONTINUE WRITE(*,80,advance='no'); READ *, F1 WRITE(*,81,advance='no'); READ *, F2 PRINT *,' ' PRINT *,' Linear frequency step.....: LIN' PRINT *,' Logarithmic frequency step: LOG' PRINT *,' ' WRITE(*,30,advance='no'); READ *, ANS PRINT *,' ' IF (ANS=='LOG') GOTO 1660 !case linear step WRITE(*,85,advance='no'); READ *, F3 GOTO 1750 1660 CONTINUE !case log. step IF (F1==0.) THEN PRINT *,' CAUTION - Starting frequency must be > 0 for log. step.' GOTO 1400 ELSE WRITE(*,86,advance='no'); READ *, F3 F3 = 10**((LOG(F2) - LOG(F1)) / LOG(10.) / F3) END IF 1750 CONTINUE !print header PRINT *,' ' PRINT *,' Frequency (Hz) Displacement Phase (deg.)' PRINT *,' --------------------------------------------' F = F1 !Initialize frequency 1760 W = 2. * PI * F DO I = 1, N I1 = I + N DO J = 1, N J1 = J + N A(I, J) = K1(I, J) - W * W * M(I, J) A(I1, J1) = A(I, J) IF (C9==1) THEN A(I1, J) = K2(I, J) * W ELSE A(I1, J) = K2(I, J) END IF A(I, J1) = -A(I1, J) END DO END DO ! A=INV(A) CALL MATINV(N2,0,A,B1,DET) ! X=A MPY B DO I = 1, N2 SUM = 0. DO J = 1, N2 SUM = SUM + A(I, J) * B(J) END DO X(I) = SUM END DO M1 = SQRT(X(N3)**2 + X(N3 + N)**2) IF (M1==0) THEN T1 = 0. ELSE O1 = X(N3) / M1 I1 = -X(N3 + N) / M1 IF (O1==-1.) THEN T1 = PI ELSE IF (O1==1.) THEN T1 = 0. ELSE IF (I1 >= 0.) THEN T1=ACOS(O1) ELSE T1=2.*PI-ACOS(O1) END IF END IF END IF T1 = T1 / PI / 180. !Convert phase in degrees !print line of results WRITE(*,90) F, M1, T1 IF (ANS=='LIN') THEN F = F + F3 ELSE F = F * F3 END IF IF (F <= F2) GOTO 1760 !next frequency PRINT *,' ' 10 FORMAT(' How many degrees of freedom (d.o.f.): ') 20 FORMAT(' M(',I2,',',I2,') = ') 30 FORMAT(' Your choice: ') 40 FORMAT(' K1(',I2,',',I2,') = ') 41 FORMAT(' K2(',I2,',',I2,') = ') 50 FORMAT(' K(',I2,',',I2,') = ') 51 FORMAT(' C(',I2,',',I2,') = ') 60 FORMAT(' F1(',I2,') = ') 61 FORMAT(' F2(',I2,') = ') 70 FORMAT(' Number of d.o.f. to calculate: ') 80 FORMAT(' Starting frequency = ') 81 FORMAT(' Ending frequency...= ') 85 FORMAT(' Frequency step: ') 86 FORMAT(' Number of constant log. increments from F1 to F2: ') 90 FORMAT(' ',F8.4,' ',F9.2,' ',F5.1) END !of main program !******************************************* !* SOLVING A LINEAR MATRIX SYSTEM AX = B * !* with Gauss-Jordan method using full * !* pivoting at each step. During the pro- * !* cess, original A and B matrices are * !* destroyed to spare storage location. * !* --------------------------------------- * !* INPUTS: A MATRIX N*N * !* B MATRIX N*M * !* --------------------------------------- * !* OUTPUS: A INVERSE OF A N*N * !* DET DETERMINANT OF A * !* B SOLUTION MATRIX N*M * !* --------------------------------------- * !* NOTA - If M=0 inversion of A matrix * !* only. * !******************************************* SUBROUTINE MATINV(N,M,AA,BB,DET) integer, parameter :: NMAX = 12 REAL*4,PARAMETER :: EPSMACH = 2.E-12 REAL*4 AA(2*NMAX,*),BB(2*NMAX,*) INTEGER PC(N),PL(N) REAL*4 CS(N) REAL*4 PV,PAV,DET,TT !Initializations : DET=1. DO I=1,N PC(I)=0 PL(I)=0 CS(I)=0. END DO !main loop DO K=1,N !Searching greatest pivot : PV=AA(K,K) IK=K JK=K PAV=ABS(PV) DO I=K,N DO J=K,N IF (ABS(AA(I,J)).GT.PAV) THEN PV=AA(I,J) PAV=ABS(PV) IK=I JK=J ENDIF ENDDO ENDDO !Search terminated, the pivot is in location I=IK, J=JK. !Memorizing pivot location: : PC(K)=JK PL(K)=IK !Determinant DET is actualised !If DET=0, ERROR MESSAGE and STOP !Machine dependant EPSMACH equals here 1.DE-20 IF (IK.NE.K) DET=-DET IF (JK.NE.K) DET=-DET DET=DET*PV IF (ABS(DET).LT.EPSMACH) THEN !Error message and Stop PRINT 10 STOP ENDIF !POSITIONNING PIVOT IN K,K: IF(IK.NE.K) THEN DO I=1,N !EXCHANGE LINES IK and K: TT=AA(IK,I) AA(IK,I)=AA(K,I) AA(K,I)=TT ENDDO ENDIF IF (M.NE.0) THEN DO I=1,M TT=BB(IK,I) BB(IK,I)=BB(K,I) BB(K,I)=TT ENDDO ENDIF !Pivot is at correct line IF(JK.NE.K) THEN DO I=1,N !Exchange columns JK and K of matrix AA TT=AA(I,JK) AA(I,JK)=AA(I,K) AA(I,K)=TT ENDDO ENDIF !Pivot is at correct column and located in K,K !Store column K in vector CS !then set column K to zero DO I=1,N CS(I)=AA(I,K) AA(I,K)=0.D0 ENDDO ! CS(K)=0. AA(K,K)=1. !Modify line K : IF(ABS(PV).LT.EPSMACH) THEN WRITE(*,*) ' PIVOT TOO SMALL - STOP' STOP ENDIF DO I=1,N AA(K,I)=AA(K,I)/PV ENDDO IF (M.NE.0) THEN DO I=1,M BB(K,I)=BB(K,I)/PV ENDDO ENDIF !Modify other lines of matrix AA: DO J=1,N IF (J.EQ.K) CONTINUE DO I=1,N !Modify line J of matrix AA : AA(J,I)=AA(J,I)-CS(J)*AA(K,I) ENDDO IF (M.NE.0) THEN DO I=1,M BB(J,I)=BB(J,I)-CS(J)*BB(K,I) ENDDO ENDIF ENDDO !Line K is ready. ENDDO !End of K loop !The matrix AA is inverted - Rearrange AA !Exchange lines DO I=N,1,-1 IK=PC(I) IF (IK.EQ.I) CONTINUE !EXCHANGE LINES I AND PC(I) OF AA: DO J=1,N TT=AA(I,J) AA(I,J)=AA(IK,J) AA(IK,J)=TT ENDDO IF (M.NE.0) THEN DO J=1,M TT=BB(I,J) BB(I,J)=BB(IK,J) BB(IK,J)=TT ENDDO ENDIF !NO MORE EXCHANGE NEEDED !GO TO NEXT LINE ENDDO !EXCHANGE COLUMNS DO J=N,1,-1 JK=PL(J) IF (JK.EQ.J) CONTINUE !EXCHANGE COLUMNS J AND PL(J) OF AA : DO I=1,N TT=AA(I,J) AA(I,J)=AA(I,JK) AA(I,JK)=TT ENDDO !NO MORE EXCHANGE NEEDED !GO TO NEXT COLUMN ENDDO !REARRANGEMENT TERMINATED. RETURN 10 FORMAT(//' DETERMINANT EQUALS ZERO, NO SOLUTION!'/) END ! end of file ndof03.f90