!*********************************************************************************** !* Response of a N d.o.f. Mass-Spring System with Damping to a sinusoidal Force * !* By Transfer Matrices Method * !* ------------------------------------------------------------------------------- * !* EXPLANATION: * !* Find response of Mass 2m in 3 d.o.f. undamped system below: * !* Nodes: 1 2 3 4 5 6 7 * !* * !* --> x1 --> x2 --> x3 * !* \| k ------ 2k ------ k ------ * !* \|--/\/\/\--| 2m |--/\/\/\--| m |--/\/\/\--| 3m | * !* \| ------ ------ ------ * !* --> F1(t) * !* * !* The differential equations of the system motions is: * !* | 2m 0 0 | | x1' | | 3k -2k 0 | | x1 | | 0 | * !* | 0 m 0 | | x2' | + |-2k 3k -k | | x2 | = | 0 | * !* | 0 0 3m | | x3' | | 0 -k k | | x3 | | 0 | * !* * !* | X1 | jwt * !* Noting X = | X2 | e the matrix form of equations becomes: * !* | X3 | * !* * !* | 3k - 2mw2 -2k 0 | | X1 | * !* | -2k 3k - mw2 -k | | X2 | = 0 * !* | 0 -k k - 3mw2 | | X3 | * !* * !* The matrix determinant is null for w1 = 0.10517 sqr(k/m), w2 = 0.80861 sqr(k/m) * !* and w3 = 3.9195 sqr(k/m), hence the resonance fresuencies 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 | * !* * !* ------------------------------------------------------------------------------- * !* SAMPLE RUN: * !* * !* Kind of elements * !* ---------------- * !* 1: Spring * !* 2: Mass * !* 3: Viscous Damper * !* 4: Spring + Viscous Damper in parallel * !* 5: Spring with structural Damping * !* 6: Sinusoidal Force * !* * !* Number of elements: 7 * !* * !* Kind of element 1: 1 * !* Kind of element 2: 2 * !* Kind of element 3: 1 * !* Kind of element 4: 2 * !* Kind of element 5: 1 * !* Kind of element 6: 2 * !* Kind of element 7: 6 * !* * !* Mass #1 = 2 * !* Mass #2 = 1 * !* Mass #3 = 3 * !* * !* Spring #1 = 1 * !* Spring #2 = 2 * !* Spring #3 = 1 * !* * !* Excitation Force #1 F.COS(PHI) = 1000 * !* Excitation Force #1 F.SIN(PHI) = 0 * !* * !* For which node number do you want the response: 3 * !* * !* Fixed-Fixed System, code=1. * !* Fixed-Free System, code=2. * !* Free-Fixed System, code=3. * !* Free-Free System, code=4. * !* * !* Code = 2 * !* * !* Frequency Sweep * !* --------------- * !* Starting Frequency: 0.30 * !* Ending Frequency..: 0.32 * !* Frequency Step....: 0.001 * !* * !* Freq= 0.300 Displacement= 96.1 Phase= 0.0 Deg. * !* Freq= 0.301 Displacement= 101.2 Phase= 0.0 Deg. * !* Freq= 0.302 Displacement= 107.1 Phase= 0.0 Deg. * !* Freq= 0.303 Displacement= 114.0 Phase= 0.0 Deg. * !* Freq= 0.304 Displacement= 122.2 Phase= 0.0 Deg. * !* Freq= 0.305 Displacement= 132.1 Phase= 0.0 Deg. * !* Freq= 0.306 Displacement= 144.1 Phase= 0.0 Deg. * !* Freq= 0.307 Displacement= 159.3 Phase= 0.0 Deg. * !* Freq= 0.308 Displacement= 178.8 Phase= 0.0 Deg. * !* Freq= 0.309 Displacement= 204.7 Phase= 0.0 Deg. * !* Freq= 0.310 Displacement= 240.9 Phase= 0.0 Deg. * !* Freq= 0.311 Displacement= 294.8 Phase= 0.0 Deg. * !* Freq= 0.312 Displacement= 383.9 Phase= 0.0 Deg. * !* Freq= 0.313 Displacement= 558.2 Phase= 0.0 Deg. * !* Freq= 0.314 Displacement= 1051.8 Phase= 0.0 Deg. * !* Freq= 0.315 Displacement= 12214.5 Phase= 0.0 Deg. <-- 3rd resonance * !* Freq= 0.316 Displacement= 1226.3 Phase=180.0 Deg. * !* Freq= 0.317 Displacement= 574.1 Phase=180.0 Deg. * !* Freq= 0.318 Displacement= 370.7 Phase=180.0 Deg. * !* Freq= 0.319 Displacement= 271.5 Phase=180.0 Deg. * !* Freq= 0.320 Displacement= 212.8 Phase=180.0 Deg. * !* * !* Example #2: * !* Nodes: 1 2 3 4 5 6 7 * !* --> x1 --> x2 k --> x3 * !* \| k *----* 2k *----*--/\/\/\--*----* * !* \|--/\/\/\--| 2m |--/\/\/\--| m | ---* | 3m | * !* \| *----* *----*---| |----*----* * !* --> F1(t) ---* * !* E * !* * !* Number of elements: 7 * !* * !* Kind of element 1: 1 * !* Kind of element 2: 2 * !* Kind of element 3: 1 * !* Kind of element 4: 2 * !* Kind of element 5: 5 * !* Kind of element 6: 2 * !* Kind of element 7: 6 * !* * !* Mass #1 = 2 * !* Mass #2 = 1 * !* Mass #3 = 3 * !* * !* Spring #1 = 1 * !* Spring #2 = 2 * !* * !* Spring with structural Damping #1 K = 1 * !* Spring with structural Damping #1 E = 0.05 * !* * !* Excitation Force #1 F.COS(PHI) = 1000 * !* Excitation Force #1 F.SIN(PHI) = 0 * !* * !* For which node number do you want the response: 3 * !* * !* Fixed-Fixed System, code=1. * !* Fixed-Free System, code=2. * !* Free-Fixed System, code=3. * !* Free-Free System, code=4. * !* * !* Code = 2 * !* * !* Frequency Sweep * !* --------------- * !* Starting Frequency: 0.30 * !* Ending Frequency..: 0.32 * !* Frequency Step....: 0.001 * !* * !* Freq= 0.300 Displacement= 95.7 Phase= 3.5 Deg. * !* Freq= 0.301 Displacement= 100.6 Phase= 3.9 Deg. * !* Freq= 0.302 Displacement= 106.3 Phase= 4.5 Deg. * !* Freq= 0.303 Displacement= 113.0 Phase= 5.1 Deg. * !* Freq= 0.304 Displacement= 120.9 Phase= 5.8 Deg. * !* Freq= 0.305 Displacement= 130.4 Phase= 6.7 Deg. * !* Freq= 0.306 Displacement= 141.8 Phase= 7.8 Deg. * !* Freq= 0.307 Displacement= 156.0 Phase= 9.1 Deg. * !* Freq= 0.308 Displacement= 173.9 Phase= 10.8 Deg. * !* Freq= 0.309 Displacement= 197.2 Phase= 13.0 Deg. * !* Freq= 0.310 Displacement= 228.3 Phase= 15.9 Deg. * !* Freq= 0.311 Displacement= 271.8 Phase= 20.1 Deg. * !* Freq= 0.312 Displacement= 334.9 Phase= 26.5 Deg. * !* Freq= 0.313 Displacement= 429.1 Phase= 37.0 Deg. * !* Freq= 0.314 Displacement= 557.6 Phase= 55.2 Deg. * !* Freq= 0.315 Displacement= 644.2 Phase= 84.1 Deg. <-- 3rd resonance * !* Freq= 0.316 Displacement= 562.8 Phase=-65.6 Deg. * !* Freq= 0.317 Displacement= 422.0 Phase=-45.6 Deg. * !* Freq= 0.318 Displacement= 317.1 Phase=-34.2 Deg. * !* Freq= 0.319 Displacement= 247.5 Phase=-27.3 Deg. * !* Freq= 0.320 Displacement= 200.3 Phase=-22.8 Deg. * !* * !* Note: with a moderate structural damping added to spring #3 (eta=5%), the 3rd * !* resonance frequency is practically unchanged, but the response is much * !* lower and the phases are completely different. * !* ------------------------------------------------------------------------------- * !'* 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. * !* (www.jpmoreau.fr) * !*********************************************************************************** Program Ndof01 INTEGER, PARAMETER :: NMAX = 25 INTEGER C9, V9, F9, C6, C3, C8 REAL*8 M2, PI, TINY REAL*8 C1(NMAX), C2(NMAX), K1(NMAX), K2(NMAX), M1(NMAX) REAL*8 K3(NMAX), E3(NMAX) INTEGER C(NMAX) REAL*8 T(5, 5), T1(5, 5), T2(5, 5), A(5, 5), V1(5), V2(5) REAL*8 B(2, 2), E1(2), E2(2), F1(5), F2(5) REAL*8 C4,D,F,F3,F4,F5,S4,SUM,T3, W,W2 CHARACTER*3 ANS PI = 4.D0 * DATAN(1.D0) TINY = 1.D-12 M9 = 0; K9 = 0; C9 = 0; V9 = 0; I9 = 0; F9 = 0 PRINT *,' PI=', PI 200 PRINT *, ' ' PRINT *, ' Kind of elements' PRINT *, ' ----------------' PRINT *, ' 1: Spring' PRINT *, ' 2: Mass' PRINT *, ' 3: Viscous Damper' PRINT *, ' 4: Spring + Viscous Damper in parallel' PRINT *, ' 5: Spring with structural Damping' PRINT *, ' 6: Sinusoidal Force' PRINT *, ' ' WRITE(*,10,advance='no'); READ *, N PRINT *, ' ' DO I = 1, N WRITE(*,20,advance='no') I; READ *, C(I) IF (C(I) == 1) THEN K9 = K9 + 1 ELSE IF (C(I) == 2) THEN M9 = M9 + 1 ELSE IF (C(I) == 3) THEN C9 = C9 + 1 ELSE IF (C(I) == 4) THEN V9 = V9 + 1 ELSE IF (C(I) == 5) THEN I9 = I9 + 1 ELSE IF (C(I) == 6) THEN F9 = F9 + 1 ELSE PRINT *, ' Unknown Element.' STOP END IF END DO IF (F9.NE.0) GOTO 680 PRINT *, ' ' PRINT *, ' There is no excitation Force, Redo!' GOTO 200 680 PRINT *, ' ' IF (M9 == 0) GOTO 770 !input M9 masses DO I = 1, M9 WRITE(*,30,advance='no') I; READ *, M1(I) END DO PRINT *, ' ' 770 IF (K9 == 0) GOTO 840 !input K9 springs DO I = 1, K9 WRITE(*,40,advance='no') I; READ *, K1(I) END DO PRINT *, ' ' 840 IF (C9 == 0) GOTO 910 !input C9 viscous Dampers DO I = 1, C9 WRITE(*,50,advance='no') I; READ *, C1(I) END DO PRINT *, ' ' 910 IF (V9 == 0) GOTO 1000 !input V9 springs + viscous Dampers in parallel DO I = 1, V9 WRITE(*,60,advance='no') I; READ *, K2(I) WRITE(*,70,advance='no') I; READ *, C2(I) END DO PRINT *, ' ' 1000 IF (I9 == 0) GOTO 1090 !input I9 springs with structural Damping DO I = 1, I9 WRITE(*,80,advance='no') I; READ *, K3(I) WRITE(*,90,advance='no') I; READ *, E3(I) END DO PRINT *, ' ' 1090 CONTINUE !input F9 sinusoidal forces DO I = 1, F9 WRITE(*,100,advance='no') I; READ *, F1(I) WRITE(*,110,advance='no') I; READ *, F2(I) END DO 1180 PRINT *, ' ' WRITE(*,120,advance='no'); READ *, N1 1200 PRINT *, ' ' PRINT *, ' Fixed-Fixed System, code=1.' PRINT *, ' Fixed-Free System, code=2.' PRINT *, ' Free-Fixed System, code=3.' PRINT *, ' Free-Free System, code=4.' PRINT *, ' ' WRITE(*,130,advance='no'); READ *, C8 1300 PRINT *, ' ' PRINT *, ' Frequency Sweep' PRINT *, ' ---------------' WRITE(*,140,advance='no'); READ *, F3 WRITE(*,150,advance='no'); READ *, F4 WRITE(*,160,advance='no'); READ *, F5 PRINT *, ' ' !end of data section !Main frequency loop DO F = F3, F4, F5 M9 = 0; K9 = 0; C9 = 0; V9 = 0; I9 = 0; F9 = 0; C6 = 1 IF (C6==N1) THEN !T2=identity matrix DO I = 1, 5 DO J = 1, 5 IF (J==I) THEN T2(I, J) = 1.D0 ELSE T2(I, J) = 0.D0 END IF END DO END DO END IF W = 2.D0 * PI * F W2 = W * W ! T=0 DO I = 1, 5 DO J = 1, 5 T(I, J) = 0.D0 END DO END DO DO I = 1, N !A=identity matrix DO K = 1, 5 DO J = 1, 5 IF (J == K) THEN A(K, J) = 1.D0 ELSE A(K, J) = 0.D0 END IF END DO END DO C3 = C(I) !Select according to kind of element IF (C3 == 1) THEN CALL S2500(K9,A,K1) ELSE IF (C3 == 2) THEN CALL S2450(M9,A,M1,W2) ELSE IF (C3 == 3) THEN CALL S2550(C9,A,C1,W) ELSE IF (C3 == 4) THEN CALL S2600(V9,A,K2,C2,W,W2) ELSE IF (C3 == 5) THEN CALL S2650(I9,A,K3,E3) ELSE IF (C3 == 6) THEN CALL S2700(F9,A,F1,F2) END IF IF (I == 1) THEN !T=A DO K = 1, 5 DO J = 1, 5 T(K, J) = A(K, J) END DO END DO ELSE !T1=A MPY T DO I1 = 1, 5 DO J = 1, 5 Sum = 0.D0 DO K = 1, 5 Sum = Sum + A(I1, K) * T(K, J) T1(I1, J) = Sum END DO END DO END DO !T=T1 DO K = 1, 5 DO J = 1, 5 T(K, J) = T1(K, J) END DO END DO END IF C6 = C6 + 1 IF (C6 == N1) THEN !T2=T DO K = 1, 5 DO J = 1, 5 T2(K, J) = T(K, J) END DO END DO END IF END DO !I loop IF (C8 == 1) THEN !Fixed-Fixed I1 = 2; J1 = 1 ELSE IF (C8 == 2) THEN !Fixed-Free I1 = 1; J1 = 1 ELSE IF (C8 == 3) THEN !Free-Free I1 = 1; J1 = 2 ELSE IF (C8 == 4) THEN !Free-Fixed I1 = 2; J1 = 2 END IF B(1, 1) = T(I1, J1) B(1, 2) = T(I1, J1 + 2) B(2, 1) = T(I1 + 2, J1) B(2, 2) = T(I1 + 2, J1 + 2) E1(1) = -T(I1, 5) E1(2) = -T(I1 + 2, 5) D = B(1, 1) * B(2, 2) - B(1, 2) * B(2, 1) !Zero Divide Protection IF (DABS(D) < TINY) THEN PRINT *, ' CAUTION, SYSTEM IS SINGULAR!' GOTO 1200 !Change limit conditions END IF E2(1) = (B(2, 2) * E1(1) - B(1, 2) * E1(2)) / D E2(2) = (B(1, 1) * E1(2) - B(2, 1) * E1(1)) / D !V1=0 DO I = 1, 5 V1(I) = 0.D0 END DO V1(5) = 1.D0 IF (J1.NE.1) THEN V1(2) = E2(1); V1(4) = E2(2) ELSE V1(1) = E2(1); V1(3) = E2(2) END IF !V2=T2 MPY V1 DO I = 1, 5 Sum = 0.D0 DO J = 1, 5 Sum = Sum + T2(I, J) * V1(J) END DO V2(I) = Sum END DO M2 = DSQRT(V2(2)**2 + V2(4)**2) IF (DABS(M2) < TINY) THEN T3 = 0.D0 ELSE C4 = V2(2) / M2 S4 = -V2(4) / M2 IF (C4 == -1.D0) THEN T3 = PI ELSE IF (C4 == 1.D0) THEN T3 = 0.D0 ELSE IF (S4 >= 0.D0) THEN T3=DACOS(C4) !T3 = DATAN(DSQRT(1.D0 - C4**2) / C4) ELSE T3=2.D0*PI-DACOS(C4) !T3 = 2.D0 * PI - DATAN(DSQRT(1.D0 - C4**2) / C4) END IF END IF END IF !Convert phase in degrees T3 = T3 / PI * 180 !PRINT *, current results WRITE(*,180) F, M2, T3 END DO !end of frequency loop PRINT *, ' ' PRINT *, ' Change limit conditions: LIM' PRINT *, ' Change frequency sweep : SWE' PRINT *, ' Change response node...: NOD' PRINT *, ' Exit program...........: EXI' PRINT *, ' ' WRITE(*,181,advance='no'); READ *, ANS IF (ANS == 'LIM') GOTO 1200 IF (ANS == 'SWE') GOTO 1300 IF (ANS == 'NOD') GOTO 1180 STOP 10 FORMAT(' Number of elements: ') 20 FORMAT(' Kind of element ',I2,': ') 30 FORMAT(' Mass #',I2,' = ') 40 FORMAT(' Spring #',I2,' = ') 50 FORMAT(' Viscous Damper #',I2,' = ') 60 FORMAT(' Spring + Viscous Damper #',I2,' K = ') 70 FORMAT(' Spring + Viscous Damper #',I2,' C = ') 80 FORMAT(' Spring with structural Damping #',I2,' K = ') 90 FORMAT(' Spring with structural Damping #',I2,' E = ') 100 FORMAT(' Excitation Force #',I2,' F.COS(PHI) = ') 110 FORMAT(' Excitation Force #',I2,' F.SIN(PHI) = ') 120 FORMAT(' For which node number do you want the response: ') 130 FORMAT(' Code = ') 140 FORMAT(' Starting Frequency: ') 150 FORMAT(' Ending Frequency..: ') 160 FORMAT(' Frequency Step....: ') 180 FORMAT('Freq=',F7.3,' Displacement=',F9.1,' Phase=',F5.1,' Deg.') 181 FORMAT(' Your answer: ') END !of main program Subroutine S2450(M9,A,M1,W2) !Mass Matrix REAL*8 A(5,5), M1(*), W2 M9 = M9 + 1 A(1, 2) = -M1(M9) * W2 A(3, 4) = A(1, 2) RETURN End Subroutine S2500(K9,A,K1) !Stiffness Matrix REAL*8 A(5,5), K1(*) K9 = K9 + 1 A(2, 1) = 1.D0 / K1(K9) A(4, 3) = A(2, 1) RETURN End Subroutine S2550(C9,A,C1,W) !Viscous damping Matrix INTEGER C9 REAL*8 A(5,5), C1(*), W C9 = C9 + 1 A(4, 1) = -1.D0 / C1(C9) / W A(2, 3) = -A(4, 1) RETURN End Subroutine S2600(V9,A,K2,C2,W,W2) !Spring + viscous damper in parallel Matrix INTEGER V9 REAL*8 A(5,5), K2(*), C2(*), W, W2 V9 = V9 + 1 A(2, 1) = K2(V9) / (K2(V9)**2 + C2(V9)**2 * W2) A(4, 3) = A(2, 1) A(4, 1) = -C2(V9) * W / (K2(V9)**2 + C2(V9)**2 * W2) A(2, 3) = -A(4, 1) RETURN End Subroutine S2650(I9,A,K3,E3) !Spring with structural damping Matrix REAL*8 A(5,5), K3(*), E3(*) I9 = I9 + 1 A(2, 1) = 1.D0 / (K3(I9) * (1.D0 + E3(I9)**2)) A(4, 3) = A(2, 1) A(4, 1) = -E3(I9) / (K3(I9) * (1.D0 + E3(I9)**2)) A(2, 3) = -A(4, 1) RETURN End Subroutine S2700(F9,A,F1,F2) !Force Matrix INTEGER F9 REAL*8 A(5,5), F1(*), F2(*) F9 = F9 + 1 A(1, 5) = -F1(F9) A(3, 5) = -F2(F9) RETURN End !end of file ndof01.f90