'***************************************************************************** '* Step by step solution of system [M] X" + [C] X' + [K] X = F(t) * '* by the "Wilson-Theta" Method * '* ------------------------------------------------------------------------- * '* SAMPLE RUN: * '* * '* 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 ndof03.bas). * '* * '* * '* How many degrees of freedom (d.o.f.): 3 * '* * '* Mass Matrix: * '* M(1,1) = 2 * '* M(1,2) = 0 * '* M(1,3) = 0 * '* M(2,2) = 1 * '* M(2,3) = 0 * '* M(3,3) = 3 * '* * '* Damping Matrix: * '* 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 * '* * '* Stiffness Matrix: * '* K(1,1) = 3 * '* K(1,2) = -2 * '* K(1,3) = 0 * '* K(2,2) = 3 * '* K(2,3) = -1 * '* K(3,3) = 1 * '* * '* Starting time (sec.) = 0 * '* Ending time......... = 60 * '* Time increment...... = 0.1 * '* * '* Starting motion X(1) = 0 * '* Starting motion X(2) = 0 * '* Starting motion X(3) = 0 * '* * '* Starting speed V(1) = 0 * '* Starting speed V(2) = 0 * '* Starting speed V(3) = 0 * '* * '* Starting acceleration G(1) = 0 * '* Starting acceleration G(2) = 0 * '* Starting acceleration G(3) = 0 * '* * '* Theta = 1.4 * '* * '* Number of force componant: 1 * '* * '* Force (maximum) = 1000 * '* Force frequency = 0.315 * '* * '* Output file ndof04.txt contains: * '* * '* Time= 0 Force= 0 * '* Node Displacement Speed Acceleration * '* --------------------------------------- * '* 1 0 0 0 * '* 2 0 0 0 * '* 3 0 0 0 * '* Time= .1 Force= 196.6306946154201 * '* 1 .1630633481454764 4.891900444364293 97.83800888728585 * '* 2 1.055011797717105D-03 3.165035393151314D-02 .6330070786302628 * '* 3 1.134321720457676D-06 3.402965161373028D-05 6.805930322746056D-04* '* Time= .2 Force= 385.5839922773965 * '* 1 1.342109522392924 20.69568389433056 218.2376601120394 * '* 2 1.215582989602615D-02 .238073481154732 3.495455465834115 * '* 3 1.681963252616181D-05 3.684703693299332D-04 6.008221322049452D-03* '* Time= .3 Force= 559.4822581021671 * '* 1 4.697087083250554 48.3460760314658 334.7701826306654 * '* 2 6.533331077055958D-02 .9444046906348329 10.6311687237679 * '* 3 1.203901950183077D-04 2.069765070002037D-03 2.801767269139262D-02* '* Time= .4 Force= 711.5356772092854 * '* 1 11.37323365004051 86.85373580923384 435.3830129246953 * '* 2 .2351846769611723 2.67517316826032 23.98420082874183 * '* 3 5.752348065811248D-04 8.104924572310808D-03 9.268551735478277D-02* '* Time= .5 Force= 835.8073613682702 * '* 1 22.36456143833293 134.2632113840701 512.8064985720298 * '* 2 .6577455128974099 6.127268700129398 45.05770980863972 * '* 3 2.103809379023869D-03 2.501311216092159D-02 .2454782344174329 * '* Time= .6 Force= 927.4451533346613 * '* 1 38.43691496520708 188.0038581094827 562.0064359362229 * '* 2 1.545162720612546 12.11509334076331 74.69878300403855 * '* 3 6.349156254191256D-03 6.506027021230676D-02 .5554649266102704 * '* --- --------------------- -------------------- ----------------- * '* Time= 59.50000000000058 Force=-998.8898749620234 * '* 1 1311.501557220317 -4027.91127756982 -5149.756315555141 * '* 2 -2554.79215720466 9594.998126819037 10332.61120309285 * '* 3 184.3997346660952 -884.9594932352101 -754.7119091791997 * '* Time= 59.60000000000058 Force=-988.651744737743 * '* 1 885.6672742053009 -4461.718119533076 -3526.380523709975 * '* 2 -1550.158169105559 10432.39282918031 6415.282844132672 * '* 3 92.74860312996083 -941.8793641546515 -383.6855092096296 * '* Time= 59.70000000000058 Force=-939.8119510859298 * '* 1 424.805036107236 -4726.111877690296 -1761.494639434426 * '* 2 -481.8064257360606 10865.00250051768 2236.91058261469 * '* 3 -2.712865272797816 -960.9010483129747 3.251826043166988 * '* Time= 59.80000000000058 Force=-854.2774316987006 * '* 1 -53.54968238474825 -4810.343067407215 76.87084509603785 * '* 2 608.7489562817452 10874.81093036808 -2040.741985606796 * '* 3 -98.14045939900252 -941.1883184623499 391.0027709693283 * '* Time= 59.90000000000058 Force=-735.387860780239 * '* 1 -531.1317036387142 -4710.618045059347 1917.629601861331 * '* 2 1679.008274813626 10460.1947945006 -6251.580731742692 * '* 3 -189.6819550233697 -883.4183703547816 764.3961911820371 * '* Time= 60.00000000000058 Force=-587.7852522915433 * '* 1 -989.652416589454 -4430.266778496569 3689.395729394218 * '* 2 2687.136442441895 9636.034476434012 -10231.62562958909 * '* 3 -273.6278425865369 -789.7596957445547 1108.7773010225 * '* * '* Maximum Displacement (node #1) = 2432.1840040250 * '* * '* ------------------------------------------------------------------------- * '* REFERENCE: "Mécanique des vibrations linéaires By M. Lalanne, * '* P. Berthier, J. Der Hagopian, Masson, Paris 1980" [BIBLI 16].* '* * '* Qwick Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '***************************************************************************** 'Program ndof04 DEFDBL A-H, O-Z DEFINT I-N PI = 4# * ATN(1#) OPEN "ndof4.txt" FOR OUTPUT AS #1 CLS PRINT INPUT " How many degrees of freedom (d.o.f.): ", N DIM M(N, N) AS DOUBLE DIM K(N, N) AS DOUBLE DIM K1(N, N) AS DOUBLE DIM B1(N, 2) DIM A(N, N), C(N, N), X0(N), V0(N), G0(N), F0(N), F1(N), F2(N) DIM X1(N), G1(N), X2(N) 'Data section PRINT PRINT " Mass Matrix:" FOR I = 1 TO N FOR J = I TO N PRINT " M("; I; ","; J; ") = "; : INPUT "", M(I, J) NEXT J NEXT I PRINT PRINT " Damping Matrix" FOR I = 1 TO N FOR J = I TO N PRINT " C("; I; ","; J; ") = "; : INPUT "", C(I, J) NEXT J NEXT I PRINT PRINT " Stiffness Matrix:" FOR I = 1 TO N FOR J = I TO N PRINT " K("; I; ","; J; ") = "; : INPUT "", K(I, J) NEXT J NEXT I 'Complete Matrices by symmetry FOR I = 1 TO N FOR J = I TO N K(J, I) = K(I, J) C(J, I) = C(I, J) M(J, I) = M(I, J) NEXT J NEXT I PRINT INPUT " Starting time (sec.) = ", T0 INPUT " Ending time......... = ", T3 INPUT " Time increment...... = ", D PRINT FOR I = 1 TO N PRINT " Starting motion X("; I; ") = "; : INPUT "", X0(I) NEXT I PRINT FOR I = 1 TO N PRINT " Starting speed V("; I; ") = "; : INPUT "", V0(I) NEXT I PRINT FOR I = 1 TO N PRINT " Starting acceleration G("; I; ") = "; : INPUT "", G0(I) NEXT I PRINT INPUT " Theta = ", T2 PRINT INPUT " Number of force componant: ", N1 PRINT INPUT " Force (maximum) = ", Fmax INPUT " Force frequency = ", freq F = Fmax * SIN(2# * PI * freq * T0) 'end of data section PRINT #1, PRINT #1, " Time="; T0; " Force="; F PRINT #1, " Node Displacement Speed Acceleration" PRINT #1, " ---------------------------------------" FOR I = 1 TO N PRINT #1, " "; I; " "; X0(I); " "; V0(I); " "; G0(I) NEXT I PRINT FOR I = 1 TO N F0(I) = 0# NEXT I F0(N1) = F A9 = T2 * D A0 = 6# / A9 ^ 2 A1 = 3# / A9 A2 = 2# * A1 A3 = A9 / 2# A4 = A0 / T2 A5 = -A2 / T2 A6 = 1# - 3# / T2 A7 = D / 2# A8 = D ^ 2 / 6# 'Construct K+A0*M+A1*C '--------------------- 500 FOR I = 1 TO N F1(I) = 0# NEXT I T1 = T0 + D F = Fmax * SIN(2# * PI * freq * T1) PRINT #1, " Time="; T1; " Force="; F F1(N1) = F FOR I = 1 TO N F2(I) = 0# NEXT I FOR I = 1 TO N F2(I) = F0(I) + T2 * (F1(I) - F0(I)) FOR J = 1 TO N K1(I, J) = K(I, J) + A0 * M(I, J) + A1 * C(I, J) F2(I) = F2(I) + M(I, J) * (A0 * X0(J) + A2 * V0(J) + 2# * G0(J)) F2(I) = F2(I) + C(I, J) * (A1 * X0(J) + 2# * V0(J) + A3 * G0(J)) NEXT J NEXT I 'K1=INV(K1) FOR I = 1 TO N FOR J = 1 TO N A(I, J) = K1(I, J) NEXT J NEXT I M = 0: GOSUB 1000 FOR I = 1 TO N FOR J = 1 TO N K1(I, J) = A(I, J) NEXT J NEXT I 'X2=K1 MPY F2 FOR I = 1 TO N SUM = 0# FOR J = 1 TO N SUM = SUM + K1(I, J) * F2(J) NEXT J X2(I) = SUM NEXT I FOR I = 1 TO N G1(I) = A4 * (X2(I) - X0(I)) + A5 * V0(I) + A6 * G0(I) V1(I) = V0(I) + A7 * (G1(I) + G0(I)) X1(I) = X0(I) + D * V0(I) + A8 * (G1(I) + 2# * G0(I)) IF ABS(X1(1)) > XMAX THEN XMAX = X1(1) NEXT I FOR I = 1 TO N PRINT #1, " "; I; " "; X1(I); " "; V1(I); " "; G1(I) NEXT I FOR I = 1 TO N X0(I) = X1(I): V0(I) = V1(I): G0(I) = G1(I) NEXT I T0 = T1 IF T1 < T3 THEN GOTO 500 'next time step PRINT PRINT " Results in file ndof4.txt..." PRINT #1, " Maximum Displacement (node #1) = "; XMAX CLOSE #1 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 ndof4.bas