'*************************************************************************
'* 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]. *
'* *
'* Quick Basic Release By J-P Moreau, Paris. *
'* (www.jpmoreau.fr) *
'*************************************************************************
'PROGRAM NDOF03
DEFDBL A-H, O-Z
DEFINT I-N
PI = 4# * ATN(1#)
CLS
PRINT
INPUT " How many degrees of freedom (d.o.f.): ", N
N2 = 2 * N
DIM M(N, N) AS DOUBLE
DIM K1(N, N) AS DOUBLE
DIM K2(N, N) AS DOUBLE
DIM A(N2, N2), X(N2), B(N2), B1(N2, 1)
DIM C9 AS INTEGER
DIM I1 AS DOUBLE
DIM M1 AS DOUBLE
PRINT
PRINT " Input Mass Matrix [M]:"
PRINT
FOR I = 1 TO N
FOR J = I TO N
PRINT " M("; I; ","; J; ") = "; : INPUT "", M(I, J)
NEXT J
NEXT I
PRINT
PRINT " Viscous Damping...: VIS"
PRINT " Structural Damping: STR"
PRINT
INPUT " Your choice: ", ANS$
IF ANS$ = "VIS" THEN GOTO 600
'case structural damping
C9 = 2
PRINT
PRINT " Input stiffness Matrix [K1]:"
FOR I = 1 TO N
FOR J = I TO N
PRINT " K1("; I; ","; J; ") = "; : INPUT "", K1(I, J)
NEXT J
NEXT I
PRINT
PRINT " Input Damping Matrix [K2]:"
FOR I = 1 TO N
FOR J = I TO N
PRINT " K2("; I; ","; J; ") = "; : INPUT "", K2(I, J)
NEXT J
NEXT I
PRINT
GOTO 900
600 'case viscous damping
C9 = 1
PRINT
PRINT " Input stiffness Matrix [K]:"
FOR I = 1 TO N
FOR J = I TO N
PRINT " K("; I; ","; J; ") = "; : INPUT "", K1(I, J)
NEXT J
NEXT I
PRINT
PRINT " Input Damping Matrix [C]:"
FOR I = 1 TO N
FOR J = I TO N
PRINT " C("; I; ","; J; ") = "; : INPUT "", K2(I, J)
NEXT J
NEXT I
PRINT
900 'Complete Matrices by symmetry
FOR I = 1 TO N
FOR J = I TO N
K1(J, I) = K1(I, J)
K2(J, I) = K2(I, J)
M(J, I) = M(I, J)
NEXT J
NEXT I
PRINT
PRINT " Input Excitation Vector:"
FOR I = 1 TO N
PRINT " F1("; I; ") = "; : INPUT " ", B(I)
NEXT I
PRINT
FOR I = N + 1 TO N2
J = I - N
PRINT " F2("; J; ") = "; : INPUT " ", B(I)
NEXT I
PRINT
INPUT " Number of d.o.f. to calculate: ", N3
PRINT
'calculate numerical response of d.o.f. #N3
PRINT
PRINT " Frequency Sweep"
PRINT " ---------------"
PRINT
1400 INPUT " Starting frequency = ", F1
INPUT " Ending frequency...= ", F2
PRINT
PRINT " Linear frequency step.....: LIN"
PRINT " Logarithmic frequency step: LOG"
PRINT
INPUT " Your choice: ", ANS$
PRINT
IF ANS$ = "LOG" THEN GOTO 1660
'case linear step
INPUT " Frequency step: ", F3
GOTO 1750
1660 'case log. step
IF F1 = 0 THEN
PRINT " CAUTION - Starting frequency must be > 0 for log. step."
GOTO 1400
ELSE
INPUT " Number of constant log. increments from F1 to F2: ", F3
F3 = 10 ^ ((LOG(F2) - LOG(F1)) / LOG(10#) / F3)
END IF
1750 'print header
PRINT
PRINT " Frequency (Hz) Displacement Phase (deg.)"
PRINT " --------------------------------------------"
F = F1 'Initialize frequency
1760 W = 2# * PI * F
FOR I = 1 TO N
I1 = I + N
FOR J = 1 TO 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)
NEXT J
NEXT I
'A=INV(A)
NOLD = N: N = N2: M = 0: GOSUB 1000 'A=INV(A)
N = NOLD
'X=A MPY B
FOR I = 1 TO N2
SUM = 0#
FOR J = 1 TO N2
SUM = SUM + A(I, J) * B(J)
NEXT J
X(I) = SUM
NEXT I
M1 = SQR(X(N3) ^ 2 + X(N3 + N) ^ 2)
IF M1 = 0 THEN
T1 = 0#: GOTO 1800
ELSE
O1 = X(N3) / M1
I1 = -X(N3 + N) / M1
IF O1 = -1# THEN
T1 = PI: GOTO 1800
ELSEIF O1 = 1# THEN
T1 = 0#: GOTO 1800
ELSE
IF I1 >= 0# THEN
'T1=ACOS(O1)
T1 = ATN(SQR(1 - O1 ^ 2) / O1)
ELSE
'T1=2*PI-ACOS(O1)
T1 = 2 * PI - ATN(SQR(1 - O1 ^ 2) / O1)
END IF
END IF
END IF
1800 T1 = T1 / PI / 180 'Convert phase in degrees
'print line of results
PRINT USING " ###.####"; F;
PRINT USING " ######.##"; M1;
PRINT USING " ###.# "; T1
IF ANS$ = "LIN" THEN
F = F + F3
ELSE
F = F * F3
END IF
IF F <= F2 THEN GOTO 1760 'next frequency
PRINT
END 'of main program
'*******************************************
'* SOLVING A LINEAR MATRIX SYSTEM AX = B1 *
'* 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 *
'* B1 MATRIX N*M *
'* --------------------------------------- *
'* OUTPUS: A INVERSE OF A N*N *
'* DET DETERMINANT OF A *
'* B1 SOLUTION MATRIX N*M *
'* --------------------------------------- *
'* NOTA - If M=0 inversion of A matrix *
'* only. *
'*******************************************
1000 'Subroutine MATINV
'Labels: 100, 200, 300
'Initializations
DET = 1#
FOR I = 1 TO N
PC(I) = 0#
PL(I) = 0#
CS(I) = 0#
NEXT I
'Main loop
FOR K = 1 TO N
'Searching greatest pivot
PV = A(K, K)
IK = K: JK = K
PAV = ABS(PV)
FOR I = K TO N
FOR J = K TO N
temp = ABS(A(I, J))
IF temp > PAV THEN
PV = A(I, J)
PAV = ABS(PV)
IK = I: JK = J
END IF
NEXT J
NEXT I
'Search terminated, the pivot is in location I=IK, J=JK.
'Memorizing pivot location:
PC(K) = JK
PL(K) = IK
'DETERMINANT DET is updated
'If DET=0, ERROR MESSAGE and STOP
'Machine dependant EPSMACH equals here 2E-16
IF IK <> K THEN
DET = -DET
END IF
IF JK <> K THEN
DET = -DET
END IF
DET = DET * PV
temp = ABS(DET)
IF temp < EPSMACH THEN
PRINT #2,
PRINT #2, " The determinant equals ZERO !!!"
END 'stop program
END IF
'POSITIONNING PIVOT IN K,K:}
IF IK <> K THEN
FOR I = 1 TO N
'EXCHANGE LINES IK and K of matrix A:}
TT = A(IK, I)
A(IK, I) = A(K, I)
A(K, I) = TT
NEXT I
END IF
IF M <> 0 THEN
FOR I = 1 TO M
TT = B1(IK, I)
B1(IK, I) = B1(K, I)
B1(K, I) = TT
NEXT I
END IF
'PIVOT is at correct line
IF JK <> K THEN
FOR I = 1 TO N
'EXCHANGE COLUMNS JK and K of matrix A:
TT = A(I, JK)
A(I, JK) = A(I, K)
A(I, K) = TT
NEXT I
END IF
'The PIVOT is at correct column.
'and is located in K,K.
'Column K of matrix A is stored in CS vector
'then column K is set to zero.
FOR I = 1 TO N
CS(I) = A(I, K)
A(I, K) = 0#
NEXT I
CS(K) = 0#
A(K, K) = 1#
'Line K of matrix A is modified:
temp = ABS(PV)
IF temp < EPSMACH THEN
PRINT #2,
PRINT #2, " PIVOT TOO SMALL - STOP"
END 'stop program
END IF
FOR I = 1 TO N
A(K, I) = A(K, I) / PV
NEXT I
IF M <> 0 THEN
FOR I = 1 TO M
B1(K, I) = B1(K, I) / PV
NEXT I
END IF
'Other lines of matrix A are modified:
FOR J = 1 TO N
IF J = K THEN
GOTO 100
END IF
FOR I = 1 TO N
'Line J of matrix A is modified:
A(J, I) = A(J, I) - CS(J) * A(K, I)
NEXT I
IF M <> 0 THEN
FOR I = 1 TO M
B1(J, I) = B1(J, I) - CS(J) * B1(K, I)
NEXT I
END IF
100 NEXT J
'Line K is ready
NEXT K 'end K loop
'MATRIX A INVERSION IS DONE - REARRANGEMENT OF MATRIX A
'EXCHANGE LINES
FOR I = N TO 1 STEP -1
IK = INT(PC(I))
IF IK = I THEN
GOTO 200
END IF
'EXCHANGE LINES I and PC(I) of matrix A:
FOR J = 1 TO N
TT = A(I, J)
A(I, J) = A(IK, J)
A(IK, J) = TT
NEXT J
IF M <> 0 THEN
FOR J = 1 TO M
TT = B1(I, J)
B1(I, J) = B1(IK, J)
B1(IK, J) = TT
NEXT J
END IF
'NO MORE EXCHANGE is NECESSARY
'Go to next line.
200 NEXT I 'End of loop i=N to 1 step -1
'EXCHANGE COLUMNS
FOR J = N TO 1 STEP -1
JK = INT(PL(J))
IF JK = J THEN
GOTO 300
END IF
'EXCHANGE COLUMNS J ET PL(J) of matrix A:
FOR I = 1 TO N
TT = A(I, J)
A(I, J) = A(I, JK)
A(I, JK) = TT
NEXT I
'NO MORE EXCHANGE is NECESSARY
'Go to next column.}
300 NEXT J
'REARRANGEMENT TERMINATED
RETURN 'End of subroutine 1000
'end of file ndof03.bas