'*********************************************************************************** '* 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]. * '* * '* Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '*********************************************************************************** 'Program Ndof01 DEFDBL A-H, O-Z DEFINT I-N DIM C9 AS INTEGER DIM V9 AS INTEGER DIM F9 AS INTEGER DIM C6 AS INTEGER DIM C3 AS INTEGER DIM M2 AS DOUBLE PI = 4# * ATN(1#) TINY = .000000000001# M9 = 0: K9 = 0: C9 = 0: V9 = 0: I9 = 0: F9 = 0 200 CLS : 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 INPUT " Number of elements: ", N DIM C(N) AS INTEGER DIM T(5, 5), T1(5, 5), T2(5, 5), A(5, 5), V1(5), V2(5) DIM B(2, 2), E1(2), E2(2) PRINT FOR I = 1 TO N PRINT " Kind of element "; I; ": "; : INPUT "", C(I) IF C(I) = 1 THEN K9 = K9 + 1 ELSEIF C(I) = 2 THEN M9 = M9 + 1 ELSEIF C(I) = 3 THEN C9 = C9 + 1 ELSEIF C(I) = 4 THEN V9 = V9 + 1 ELSEIF C(I) = 5 THEN I9 = I9 + 1 ELSEIF C(I) = 6 THEN F9 = F9 + 1 ELSE PRINT " Unknown Element." END END IF NEXT I IF F9 <> 0 THEN GOTO 680 PRINT PRINT " There is no excitation Force, Redo!" GOTO 200 680 PRINT IF M9 = 0 THEN GOTO 770 'input M9 masses DIM M1(M9) AS DOUBLE FOR I = 1 TO M9 PRINT " Mass #"; I; " = "; : INPUT "", M1(I) NEXT I PRINT 770 IF K9 = 0 THEN GOTO 840 'input K9 springs DIM K1(K9) AS DOUBLE FOR I = 1 TO K9 PRINT " Spring #"; I; " = "; : INPUT "", K1(I) NEXT I PRINT 840 IF C9 = 0 THEN GOTO 910 'input C9 viscous Dampers DIM C1(C9) FOR I = 1 TO C9 PRINT " Viscous Damper #"; I; " = "; : INPUT "", C1(I) NEXT I PRINT 910 IF V9 = 0 THEN GOTO 1000 'input V9 springs + viscous Dampers in parallel DIM K2(V9) AS DOUBLE DIM C2(V9) FOR I = 1 TO V9 PRINT " Spring + Viscous Damper #"; I; " K = "; : INPUT "", K2(I) PRINT " Spring + Viscous Damper #"; I; " C = "; : INPUT "", C2(I) NEXT I PRINT 1000 IF I9 = 0 THEN GOTO 1090 'input I9 springs with structural Damping DIM K3(I9) AS DOUBLE DIM E3(I9) FOR I = 1 TO I9 PRINT " Spring with structural Damping #"; I; " K = "; : INPUT "", K3(I) PRINT " Spring with structural Damping #"; I; " E = "; : INPUT "", E3(I) NEXT I PRINT 1090 'input F9 sinusoidal forces DIM F1(F9), F2(F9) FOR I = 1 TO F9 PRINT " Excitation Force #"; I; " F.COS(PHI) = "; : INPUT "", F1(I) PRINT " Excitation Force #"; I; " F.SIN(PHI) = "; : INPUT "", F2(I) NEXT I 1180 PRINT INPUT " For which node number do you want the response: ", 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 INPUT " Code = ", C8 1300 PRINT PRINT " Frequency Sweep" PRINT " ---------------" INPUT " Starting Frequency: ", F3 INPUT " Ending Frequency..: ", F4 INPUT " Frequency Step....: ", F5 PRINT 'end of data section 'Main frequency loop FOR F = F3 TO F4 STEP F5 M9 = 0: K9 = 0: C9 = 0: V9 = 0: I9 = 0: F9 = 0: C6 = 1 IF C6 = N1 THEN 'T2=identity matrix FOR I = 1 TO 5 FOR J = 1 TO 5 IF J = I THEN T2(I, J) = 1# ELSE T2(I, J) = 0# END IF NEXT J NEXT I END IF W = 2 * PI * F W2 = W * W 'T=0# FOR I = 1 TO 5 FOR J = 1 TO 5 T(I, J) = 0# NEXT J NEXT I FOR I = 1 TO N 'A=identity matrix FOR K = 1 TO 5 FOR J = 1 TO 5 IF J = K THEN A(K, J) = 1# ELSE A(K, J) = 0# END IF NEXT J NEXT K C3 = C(I) 'Select according to kind of element IF C3 = 1 THEN GOSUB 2500 ELSEIF C3 = 2 THEN GOSUB 2450 ELSEIF C3 = 3 THEN GOSUB 2550 ELSEIF C3 = 4 THEN GOSUB 2600 ELSEIF C3 = 5 THEN GOSUB 2650 ELSEIF C3 = 6 THEN GOSUB 2700 END IF IF I = 1 THEN 'T=A FOR K = 1 TO 5 FOR J = 1 TO 5 T(K, J) = A(K, J) NEXT J NEXT K ELSE 'T1=A MPY T FOR I1 = 1 TO 5 FOR J = 1 TO 5 Sum = 0# FOR K = 1 TO 5 Sum = Sum + A(I1, K) * T(K, J) T1(I1, J) = Sum NEXT K NEXT J NEXT I1 'T=T1 FOR K = 1 TO 5 FOR J = 1 TO 5 T(K, J) = T1(K, J) NEXT J NEXT K END IF C6 = C6 + 1 IF C6 = N1 THEN 'T2=T FOR K = 1 TO 5 FOR J = 1 TO 5 T2(K, J) = T(K, J) NEXT J NEXT K END IF NEXT I IF C8 = 1 THEN 'Fixed-Fixed I1 = 2: J1 = 1 ELSEIF C8 = 2 THEN 'Fixed-Free I1 = 1: J1 = 1 ELSEIF C8 = 3 THEN 'Free-Free I1 = 1: J1 = 2 ELSEIF 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 ABS(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 FOR I = 1 TO 5 V1(I) = 0# NEXT I V1(5) = 1# IF J1 <> 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 FOR I = 1 TO 5 Sum = 0# FOR J = 1 TO 5 Sum = Sum + T2(I, J) * V1(J) NEXT J V2(I) = Sum NEXT I M2 = SQR(V2(2) ^ 2 + V2(4) ^ 2) IF ABS(M2) < TINY THEN T3 = 0# ELSE C4 = V2(2) / M2 S4 = -V2(4) / M2 IF C4 = -1# THEN T3 = PI ELSEIF C4 = 1# THEN T3 = 0# ELSE IF S4 >= 0# THEN 'T3=ACOS(C4) T3 = ATN(SQR(1 - C4 ^ 2) / C4) ELSE 'T3=2*PI-ACOS(C4) T3 = 2 * PI - ATN(SQR(1 - C4 ^ 2) / C4) END IF END IF END IF 'Convert phase in degrees T3 = T3 / PI * 180 'print current results PRINT USING " Freq=###.###"; F; PRINT USING " Displacement=#######.#"; M2; PRINT USING " Phase=###.# Deg."; T3 NEXT F 'end main loop PRINT PRINT " Change limit conditions: LIM" PRINT " Change frequency sweep : SWE" PRINT " Change response node...: NOD" PRINT " Exit program...........: EXI" PRINT INPUT " Your answer: ", A\$ IF A\$ = "LIM" THEN GOTO 1200 IF A\$ = "SWE" THEN GOTO 1300 IF A\$ = "NOD" THEN GOTO 1180 END 'of main program 2450 'Mass Matrix M9 = M9 + 1 A(1, 2) = -M1(M9) * W2 A(3, 4) = A(1, 2) RETURN 2500 'Stiffness Matrix K9 = K9 + 1 A(2, 1) = 1# / K1(K9) A(4, 3) = A(2, 1) RETURN 2550 'Viscous damping Matrix C9 = C9 + 1 A(4, 1) = -1# / C1(C9) / W A(2, 3) = -A(4, 1) RETURN 2600 'Spring + viscous damper in parallel Matrix 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 2650 'Spring with structural damping Matrix I9 = I9 + 1 A(2, 1) = 1# / (K3(I9) * (1# + E3(I9) ^ 2)) A(4, 3) = A(2, 1) A(4, 1) = -E3(I9) / (K3(I9) * (1# + E3(I9) ^ 2)) A(2, 3) = -A(4, 1) RETURN 2700 'Force Matrix F9 = F9 + 1 A(1, 5) = -F1(F9) A(3, 5) = -F2(F9) RETURN 'end of file ndof01.bas