/********************************************************************************** * 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]. * * * * C++ Release By J-P Moreau, Paris. * * (www.jpmoreau.fr) * **********************************************************************************/ #include #include #define NMAX 25 #define TINY 1e-12 //Labels: e200,e680,e770,e840,e910,e1000,e1090,e1180,e1200,e1300 typedef double MAT5[6][6]; typedef double VEC[NMAX]; typedef double VEC5[6]; int N,N1, I,I1, ITMP, J,J1, K, I9, M9, K9, C3, C6, C8, C9, V9, F9; double C4,D,F,F3,F4,F5, M2, PI, S4, Sum, T3,TMP, W,W2; int C[NMAX]; MAT5 T,T1,T2,A; VEC5 V1,V2; double B[3][3]; double E1[3], E2[3]; VEC C1,C2,E3,F1,F2,K1,K2,K3,M1; char Ans[4]; double Sqr(double x) { return x*x; } void S2450() { //Mass Matrix M9++; A[1][2] = -M1[M9] * W2; A[3][4] = A[1][2]; } void S2500() { //Stiffness Matrix K9++; A[2][1] = 1.0 / K1[K9]; A[4][3] = A[2][1]; } void S2550() { //Viscous damping Matrix C9++; A[4][1] = -1.0 / C1[C9] / W; A[2][3] = -A[4][1]; } void S2600() { //Spring + viscous damper in parallel Matrix V9++; A[2][1] = K2[V9] / (Sqr(K2[V9]) + Sqr(C2[V9]) * W2); A[4][3] = A[2][1]; A[4][1] = -C2[V9] * W / (Sqr(K2[V9]) + Sqr(C2[V9]) * W2); A[2][3] = -A[4][1]; } void S2650() { //Spring with structural damping Matrix I9++; A[2][1] = 1.0 / (K3[I9] * (1.0 + Sqr(E3[I9]))); A[4][3] = A[2][1]; A[4][1] = -E3[I9] / (K3[I9] * (1.0 + Sqr(E3[I9]))); A[2][3] = -A[4][1]; } void S2700() { //Force Matrix F9++; A[1][5] = -F1[F9]; A[3][5] = -F2[F9]; } void main() { PI = 4.0*atan(1.0); M9 = 0; K9 = 0; C9 = 0; V9 = 0; I9 = 0; F9 = 0; e200: printf("\n"); printf(" Kind of elements\n"); printf(" ----------------\n"); printf(" 1; Spring\n"); printf(" 2; Mass\n"); printf(" 3; Viscous Damper^n"); printf(" 4; Spring + Viscous Damper in parallel\n"); printf(" 5; Spring with structural Damping\n"); printf(" 6; Sinusoidal Force\n\n"); printf(" Number of elements; "); scanf("%d", &N); printf("\n"); for (I=1; I<=N; I++) { printf(" Kind of element %d: ", I); scanf("%d", &ITMP); C[I]=ITMP; if (C[I] == 1) K9++; else if (C[I] == 2) M9++; else if (C[I] == 3) C9++; else if (C[I] == 4) V9++; else if (C[I] == 5) I9++; else if (C[I] == 6) F9++; else { printf(" Unknown Element.\n\n"); return; } } if (F9 != 0) goto e680; printf("\n There is no excitation Force, Redo!\n\n"); goto e200; e680: printf("\n"); if (M9 == 0) goto e770; //input M9 masses for (I=1; I<=M9; I++) { printf(" Mass #%d = ", I); scanf("%lf", &TMP); M1[I]=TMP; } printf("\n"); e770: if (K9 == 0) goto e840; //input K9 springs for (I=1; I<=K9; I++) { printf(" Spring #%d = ", I); scanf("%lf", &TMP); K1[I]=TMP; } printf("\n"); e840: if (C9 == 0) goto e910; //input C9 viscous Dampers for (I=1; I<=C9; I++) { printf(" Viscous Damper #%d = ", I); scanf("%lf", &TMP); C1[I]=TMP; } printf("\n"); e910: if (V9 == 0) goto e1000; //input V9 springs + viscous Dampers in parallel for (I=1; I<=V9; I++) { printf(" Spring + Viscous Damper #%d K = ", I); scanf("%lf", &TMP); K2[I]=TMP;; printf(" Spring + Viscous Damper #%d C = ", I); scanf("%lf", &TMP); C2[I]=TMP; } printf("\n"); e1000: if (I9 == 0) goto e1090; //input I9 springs with structural Damping for (I=1; I<=I9; I++) { printf(" Spring with structural Damping #%d K = ", I); scanf("%lf", &TMP); K3[I]=TMP; printf(" Spring with structural Damping #%d E = ", I); scanf("%lf", &TMP); E3[I]=TMP; } printf("\n"); e1090:; //input F9 sinusoidal forces for (I=1; I<=F9; I++) { printf(" Excitation Force #d F.COS(PHI) = ", I); scanf("%lf", &TMP); F1[I]=TMP; printf(" Excitation Force #d F.SIN(PHI) = ", I); scanf("%lf", &TMP); F2[I]=TMP; } e1180: printf("\n"); printf(" For which node number do you want the response: "); scanf("%d", &N1); e1200: printf("\n"); printf(" Fixed-Fixed System, code=1.\n"); printf(" Fixed-Free System, code=2.\n"); printf(" Free-Fixed System, code=3.\n"); printf(" Free-Free System, code=4.\n\n"); printf(" Code = "); scanf("%d", &C8); e1300: printf("\n"); printf(" Frequency Sweep\n"); printf(" ---------------\n"); printf(" Starting Frequency; "); scanf("%lf", &F3); printf(" Ending Frequency..; "); scanf("%lf", &F4); printf(" Frequency Step....; "); scanf("%lf", &F5); printf("\n"); //end of data section //Main frequency loop F=F3; while (F <= F4) { M9 = 0; K9 = 0; C9 = 0; V9 = 0; I9 = 0; F9 = 0; C6 = 1; if (C6 == N1) //T2=identity matrix for (I=1; I<6; I++) for (J=1; J<6; J++) if (J==I) T2[I][J] = 1.0; else T2[I][J] = 0.0; W = 2 * PI * F; W2 = W * W; //T=0 for (I=1; I<6; I++) for (J=1; J<6; J++) T[I][J] = 0.0; for (I=1; I<=N; I++) { //A=identity matrix for (K=1; K<6; K++) for (J=1; J<6; J++) if (J==K) A[K][J] = 1.0; else A[K][J] = 0.0; C3 = C[I]; //Select according to kind of element if (C3 == 1) S2500(); else if (C3 == 2) S2450(); else if (C3 == 3) S2550(); else if (C3 == 4) S2600(); else if (C3 == 5) S2650(); else if (C3 == 6) S2700(); if (I == 1) // T=A for (K=1; K<6; K++) for (J=1; J<6; J++) T[K][J] = A[K][J]; else { // T1=A MPY T for (I1=1; I1<6; I1++) for (J=1; J<6; J++) { Sum = 0.0; for (K=1; K<6; K++) { Sum += A[I1][K] * T[K][J]; T1[I1][J] = Sum; } } // T=T1 for (K=1; K<6; K++) for (J=1; J<6; J++) T[K][J] = T1[K][J]; } C6++; if (C6 == N1) // T2=T for (K=1; K<6; K++) for (J=1; J<6; J++) T2[K][J] = T[K][J]; } //I loop if (C8 == 1) { //Fixed-Fixed I1 = 2; J1 = 1; } else if (C8 == 2) { //Fixed-Free I1 = 1; J1 = 1; } else if (C8 == 3) { //Free-Free I1 = 1; J1 = 2; } else if (C8 == 4) { //Free-Fixed I1 = 2; J1 = 2; } 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 (fabs(D) < TINY) { printf(" CAUTION, SYSTEM IS SINGULAR!\n\n"); goto e1200; // Change limit conditions } 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; I<6; I++) V1[I] = 0.0; V1[5] = 1.0; if (J1 != 1) { V1[2] = E2[1]; V1[4] = E2[2]; } else { V1[1] = E2[1]; V1[3] = E2[2]; } // V2=T2 MPY V1 for (I=1; I<6; I++) { Sum = 0.0; for (J=1; J<6; J++) Sum += T2[I][J] * V1[J]; V2[I] = Sum; } M2 = sqrt(Sqr(V2[2]) + Sqr(V2[4])); if (fabs(M2) < TINY) T3 = 0.0; else { C4 = V2[2] / M2; S4 = -V2[4] / M2; if (C4 == -1.0) T3 = PI; else if (C4 == 1.0) T3 = 0.0; else if (S4 >= 0.0) // T3=ACOS(C4) T3 = atan(sqrt(1.0 - C4*C4) / C4); else // T3=2*PI-ACOS(C4) T3 = 2.0 * PI - atan(sqrt(1.0 - C4*C4) / C4); } // Convert phase in degrees T3 = T3 / PI * 180.0; // printf current results printf(" Freq=%7.3f Displacement= %7.1f Phase=%5.1f Deg.\n", F, M2, T3); F += F5; } //while (F<=F4) //end frequency loop printf("\n"); printf(" Change limit conditions; LIM\n"); printf(" Change frequency sweep ; SWE\n"); printf(" Change response node...; NOD\n"); printf(" Exit program...........; EXI\n\n"); printf(" Your answer: "); scanf("%s", Ans); if (Ans[0] == 'L') goto e1200; if (Ans[0] == 'S') goto e1300; if (Ans[0] == 'N') goto e1180; } //end of file ndof01.cpp