/************************************************************************ * 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.80 0.0 * * 0.3140 1041.65 0.0 * * 0.3145 1824.42 0.0 * * 0.3150 4347.56 0.0 <-- 3rd resonance * * 0.3155 2248.74 -0.0 * * 0.3160 1151.46 -0.0 * * 0.3165 757.64 -0.0 * * 0.3170 560.61 -0.0 * * 0.3175 443.06 -0.0 * * 0.3180 365.13 -0.0 * * 0.3185 309.74 -0.0 * * 0.3190 268.36 -0.0 * * 0.3195 236.29 -0.0 * * 0.3200 210.72 -0.0 * * * * --------------------------------------------------------------------- * * 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 // Labels: e600,e900,e1400,e1660,e1750,e1760,e1800. #define NMAX 13 #define MACH_EPS 1.2e-16 #define PI 4.0*atan(1.0) #define ABS fabs typedef double MAT[NMAX][NMAX]; typedef double MAT1[2*NMAX][2]; typedef double MAT2[2*NMAX][2*NMAX]; typedef double VEC2[2*NMAX]; MAT M, K1, K2; MAT2 A; VEC2 X, B; MAT1 B1; int C9,I,I1,J,J1,N,N3; double DET,F,F1,F2,F3, M1,O1,S1,SUM,T1,tmp,W; char ANS[4]; /****************************************** * SOLVING A LINEAR MATRIX SYSTEM AX = B * * with Gauss-Jordan method using full * * pivoting at each step. During the pro- * * cess, original AA and BB matrices are * * destroyed to spare storage location. * * --------------------------------------- * * INPUTS: AA MATRIX N*N * * BB MATRIX N*M * * --------------------------------------- * * OUTPUS: AA INVERSE OF AA N*N * * DET DETERMINANT OF AA * * BB SOLUTION MATRIX N*M * * --------------------------------------- * * NOTA - If M=0 inversion of AA matrix * * only. * ******************************************/ void MATINV(int N,int M,MAT2 AA,MAT1 BB,double *DET) { int PC[2*NMAX],PL[2*NMAX]; double CS[2*NMAX]; double PV,PAV,temp,TT; int I,IK,J,JK,K; // Initializations : *DET= 1.0; for (I=1; I<=N; I++) { PC[I]= 0; PL[I]= 0; CS[I]= 0.0; } // Main loop for (K=1; K<=N; K++) { // Searching greatest pivot : PV=AA[K][K]; IK=K; JK=K; PAV= ABS(PV); for (I=K; I<=N; I++) for (J=K; J<=N; J++) { temp = ABS(AA[I][J]); if (temp > PAV) { PV=AA[I][J]; PAV= ABS(PV); IK=I; JK=J; } } // 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 1e-20 if (IK!=K) *DET=-*DET; if (JK!=K) *DET=-*DET; *DET=*DET*PV; temp= ABS(*DET); if (temp < MACH_EPS) { printf("\n The determinant equals ZERO !!!\n"); return; } // POSITIONNING PIVOT IN K,K: if(IK!=K) for (I=1; I<=N; I++) { // EXCHANGE LINES IK and K of matrix AA: TT=AA[IK][I]; AA[IK][I]=AA[K][I]; AA[K][I]=TT; } if(M!=0) for (I=1; I<=M; I++) { TT=BB[IK][I]; BB[IK][I]=BB[K][I]; BB[K][I]=TT; } // PIVOT is at correct line if(JK!=K) for (I=1; I<=N; I++) { // EXCHANGE COLUMNS JK and K of matrix AA: TT=AA[I][JK]; AA[I][JK]=AA[I][K]; AA[I][K]=TT; } // The PIVOT is at correct column. // and is located in K,K. // Column K of matrix AA is stored in CS vector // then column K is set to zero. for (I=1; I<=N; I++) { CS[I]=AA[I][K]; AA[I][K]= 0.0; } CS[K]= 0.0; AA[K][K]= 1.0; // Line K of matrix AA is modified: temp= ABS(PV); if(temp < MACH_EPS) { printf("\n PIVOT TOO SMALL - STOP\n"); return; } for (I=1; I<=N; I++) AA[K][I]=AA[K][I]/PV; if (M!=0) for (I=1; I<=M; I++) BB[K][I]=BB[K][I]/PV; // Other lines of matrix AA are modified: for (J=1; J<=N; J++) { if (J==K) J++; for (I=1; I<=N; I++) // Line J of matrix AA is modified: AA[J][I]=AA[J][I]-CS[J]*AA[K][I]; if (M!=0) for (I=1; I<=M; I++) BB[J][I]=BB[J][I]-CS[J]*BB[K][I]; } // Line K is ready } //End of K loop // MATRIX AA INVERSION IS DONE - REARRANGEMENT OF MATRIX AA // EXCHANGE LINES for (I=N; I>0; I--) { IK=PC[I]; if (IK==I) goto fin1; // EXCHANGE LINES I and PC(I) of matrix AA: for (J=1; J<=N; J++) { TT=AA[I][J]; AA[I][J]=AA[IK][J]; AA[IK][J]=TT; } if (M!=0) for (J=1; J<=M; J++) { TT=BB[I][J]; BB[I][J]=BB[IK][J]; BB[IK][J]=TT; } // NO MORE EXCHANGE is NECESSARY // Go to next line fin1: ; } // for i // EXCHANGE COLUMNS: for (J=N; J>0; J--) { JK=PL[J]; if (JK==J) goto fin2; // EXCHANGE COLUMNS J ET PL(J) of matrix AA: for (I=1; I<=N; I++) { TT=AA[I][J]; AA[I][J]=AA[I][JK]; AA[I][JK]=TT; } fin2: ; // NO MORE EXCHANGE is NECESSARY // Go to next column. } // REARRANGEMENT TERMINATED. return; } //Matinv() void main() { printf("\n How many degrees of freedom (d.o.f.): "); scanf("%d", &N); printf("\n Input Mass Matrix [M]:\n\n"); for (I=1; I<=N; I++) for (J=I; J<=N; J++) { printf(" M(%d,%d) = ",I,J); scanf("%lf", &tmp); M[I][J] = tmp; } printf("\n"); printf(" Viscous Damping...: VIS\n"); printf(" Structural Damping: STR\n\n"); printf(" Your choice: "); scanf("%s", ANS); if (ANS[0] == 'V') goto e600; // case structural damping C9 = 2; printf("\n Input stiffness Matrix [K1]:\n"); for (I=1; I<=N; I++) for (J=I; J<=N; J++) { printf(" K1(%d,%d) = ",I,J); scanf("%lf", &tmp); K1[I][J] = tmp; } printf("\n Input Damping Matrix [K2]:\n"); for (I=1; I<=N; I++) for (J=I; J<=N; J++) { printf(" K2(%d,%d) = ",I,J); scanf("%lf", &tmp); K2[I][J] = tmp; } printf("\n"); goto e900; e600:; // case viscous damping C9 = 1; printf("\n Input stiffness Matrix [K]:\n"); for (I=1; I<=N; I++) for (J=I; J<=N; J++) { printf(" K(%d,%d) = ",I,J); scanf("%lf", &tmp); K1[I][J] = tmp; } printf("\n Input Damping Matrix [C]:\n"); for (I=1; I<=N; I++) for (J=I; J<=N; J++) { printf(" C(%d,%d) = ",I,J); scanf("%lf", &tmp); K2[I][J] = tmp; } e900:; // Complete Matrices by symmetry for (I=1; I<=N; I++) for (J=I; J<=N; J++) { K1[J][I] = K1[I][J]; K2[J][I] = K2[I][J]; M[J][I] = M[I][ J]; } printf("\n Input Excitation Vector:\n"); for (I=1; I<=N; I++) { printf(" F1(%d) = ", I); scanf("%lf", &tmp); B[I]=tmp; } printf("\n"); for (I=N+1; I<=2*N; I++) { J = I - N; printf(" F2(%d) = ", J); scanf("%lf", &tmp); B[I]=tmp; } printf("\n Number of d.o.f. to calculate: "); scanf("%d", &N3); printf("\n\n"); // calculate numerical response of d.o.f. #N3 printf(" Frequency Sweep\n"); printf(" ---------------\n\n"); e1400: printf(" Starting frequency: "); scanf("%lf", &F1); printf(" Ending frequency..: "); scanf("%lf", &F2); printf("\n"); printf(" Linear frequency step.....: LIN\n"); printf(" Logarithmic frequency step: LOG\n"); printf("\n Your choice: "); scanf("%s", ANS); printf("\n"); if (ANS[2] == 'G') goto e1660; // case linear step printf(" Frequency step: "); scanf("%lf", &F3); goto e1750; e1660:; // case log. step if (F1 == 0.0) { printf(" CAUTION - Starting frequency must be > 0 for log. step.\n"); goto e1400; } else { printf(" Number of constant log. increments from F1 to F2: "); scanf("%lf", &F3); F3 = pow(10.0, ((log(F2) - log(F1)) / log(10.0) / F3)); } e1750:; // write header printf("\n"); printf(" Frequency (Hz) Displacement Phase (deg.)\n"); printf(" --------------------------------------------\n"); F = F1; // Initialize frequency e1760: W = 2.0 * PI * F; for (I=1; I<=N; I++) { I1 = I + N; for (J=1; J<=N; J++) { J1 = J + N; A[I][J] = K1[I][J] - W * W * M[I][J]; A[I1][J1] = A[I][J]; if (C9 == 1) A[I1][J] = K2[I][J] * W; else A[I1][J] = K2[I][J]; A[I][J1] = -A[I1][J]; } } // A=INV(A) MATINV(2*N,0,A,B1,&DET); // X=A MPY B for (I=1; I<=2*N; I++) { SUM = 0.0; for (J=1; J<=2*N; J++) SUM += A[I][J] * B[J]; X[I] = SUM; } M1 = sqrt(X[N3]*X[N3] + X[N3+N]*X[N3+N]); if (M1 == 0.0) T1 = 0.0; else { O1 = X[N3] / M1; S1 = -X[N3+N] / M1; if (O1 == -1.0) T1 = PI; else if (O1 == 1.0) T1 = 0.0; else { if (S1 >= 0.0) // T1=ACOS(O1) T1 = atan(sqrt(1.0 - O1*O1) / O1); else // T1=2*PI-ACOS(O1) T1 = 2.0 * PI - atan(sqrt(1.0 - O1*O1) / O1); } } T1 = T1 / PI / 180.0; //Convert phase in degrees // write line of results printf(" %8.4f %9.2f %5.1f\n", F, M1, T1); if (ANS[2] == 'G') F *= F3; else F += F3; if (F <= F2) goto e1760; // next frequency printf("\n"); } // end of file ndof03.cpp