!******************************************************************************* !* Resonance Frequencies of a bending beam * !* - Modal Mass and Stiffness * !* - Deformation modes and Maximum Strain * !* --------------------------------------------------------------------------- * !* SAMPLE RUN: * !* * !* Fixed-Free Beam, M=1 * !* Fixed-Supported, M=2 * !* Fixed-Fixed, M=3 * !* Free-Free, M=4 * !* Free-Supported, M=5 * !* Supported-Supported, M=6 * !* * !* Limit Conditions, M= 3 * !* * !* Young Modulus: 2E11 * !* * !* Volumic Mass.: 7800 * !* * !* Beam Width....: 0.4 * !* Beam Thickness: 0.5 * !* Beam Length...: 10 * !* * !* * !* Frequency (Hz) = 26.025373 * !* Frequency (Hz) = 71.739944 * !* Frequency (Hz) = 140.638973 * !* Frequency (Hz) = 232.483361 * !* Frequency (Hz) = 347.290041 * !* * !* * !* Do you want the modes, modal masses & Stiffnesses (y/n): y * !* * !* How many modes (Maximum 5): 5 * !* * !* How many points for deviation, slope & max. strain: 11 * !* * !* Do you want automatic divisions (y/n): y * !* * !* * !* MODE #1 * !* ------- * !* * !* OMEGA= 163.5222 FREQUENCY= 26.0254 * !* * !* Modal Mass= 5.40002743637149E+0006 Modal Stiffness= 1.44394161046433E+0011 * !* * !* X Deviation Slope Strain (x10^6) * !* ------------------------------------------------ * !* 0.0000 0.0000 0.0000 22771.74 * !* 1.0000 0.1925 0.3498 12232.17 * !* 2.0000 0.6304 0.4927 2225.41 * !* 3.0000 1.1155 0.4494 -6194.01 * !* 4.0000 1.4814 0.2635 -11846.94 * !* 5.0000 1.6164 0.0000 -13841.17 * !* 6.0000 1.4814 -0.2635 -11846.94 * !* 7.0000 1.1155 -0.4494 -6194.01 * !* 8.0000 0.6304 -0.4927 2225.41 * !* 9.0000 0.1925 -0.3498 12232.17 * !* 10.0000 -0.0000 -0.0000 22771.74 * !* * !* MODE #2 * !* ------- * !* * !* OMEGA= 450.7554 FREQUENCY= 71.7399 * !* * !* Modal Mass= 1.64384471757779E+009 Modal Stiffness= 3.33997017195899E+014 * !* * !* X Deviation Slope Strain (x10^6) * !* ------------------------------------------------ * !* 0.0000 0.0000 0.0000 61624.92 * !* 1.0000 0.4554 0.7516 14015.32 * !* 2.0000 1.2058 0.6202 -24480.45 * !* 3.0000 1.5043 -0.0778 -40796.43 * !* 4.0000 1.0338 -0.8261 -29766.62 * !* 5.0000 0.0000 -1.1411 -0.00 * !* 6.0000 -1.0338 -0.8261 29766.62 * !* 7.0000 -1.5043 -0.0778 40796.43 * !* 8.0000 -1.2058 0.6202 24480.46 * !* 9.0000 -0.4554 0.7516 -14015.31 * !* 10.0000 -0.0000 0.0000 -61624.92 * !* * !* MODE #3 * !* ------- * !* * !* OMEGA= 883.6607 FREQUENCY= 140.6390 * !* * !* Modal Mass= 6.30221918268558E+011 Modal Stiffness= 4.92112747787204E+017 * !* * !* X Deviation Slope Strain (x10^6) * !* ------------------------------------------------ * !* 0.0000 0.0000 0.0000 120907.45 * !* 1.0000 0.7701 1.1128 -6282.87 * !* 2.0000 1.5079 0.1215 -77726.96 * !* 3.0000 0.8687 -1.2982 -47991.93 * !* 4.0000 -0.6284 -1.3976 39638.73 * !* 5.0000 -1.4060 -0.0000 85988.24 * !* 6.0000 -0.6284 1.3976 39638.73 * !* 7.0000 0.8687 1.2982 -47991.92 * !* 8.0000 1.5079 -0.1215 -77726.96 * !* 9.0000 0.7701 -1.1128 -6282.87 * !* 10.0000 -0.0000 -0.0000 120907.44 * !* * !* MODE #4 * !* ------- * !* * !* OMEGA=1460.7360 FREQUENCY= 232.4834 * !* * !* Modal Mass= 2.62456494665287E+014 Modal Stiffness= 5.60016486006314E+020 * !* * !* X Deviation Slope Strain (x10^6) * !* ------------------------------------------------ * !* 0.0000 0.0000 0.0000 199859.15 * !* 1.0000 1.0745 1.2736 -58760.62 * !* 2.0000 1.3192 -0.9913 -120007.61 * !* 3.0000 -0.4227 -1.9219 45103.94 * !* 4.0000 -1.3935 0.3075 139911.06 * !* 5.0000 -0.0000 1.9969 0.02 * !* 6.0000 1.3935 0.3075 -139911.05 * !* 7.0000 0.4227 -1.9219 -45103.98 * !* 8.0000 -1.3192 -0.9913 120007.59 * !* 9.0000 -1.0745 1.2736 58760.66 * !* 10.0000 -0.0000 0.0000 -199859.09 * !* * !* MODE #5 * !* ------- * !* * !* OMEGA=2182.0877 FREQUENCY= 347.2900 * !* * !* Modal Mass= 1.14990481414791E+017 Modal Stiffness= 5.47527941977509E+023 * !* * !* X Deviation Slope Strain (x10^6) * !* ------------------------------------------------ * !* 0.0000 0.0000 0.0000 298555.55 * !* 1.0000 1.3218 1.1293 -144271.18 * !* 2.0000 0.6736 -2.2318 -91130.35 * !* 3.0000 -1.3394 -0.7648 201616.08 * !* 4.0000 -0.2202 2.4118 33178.41 * !* 5.0000 1.4146 0.0000 -211057.80 * !* 6.0000 -0.2202 -2.4118 33178.38 * !* 7.0000 -1.3394 0.7648 201616.09 * !* 8.0000 0.6736 2.2318 -91130.32 * !* 9.0000 1.3218 -1.1293 -144271.20 * !* 10.0000 -0.0000 -0.0000 298555.50 * !* * !* --------------------------------------------------------------------------- * !* REFERENCE: "Mécanique des vibrations linéaires By M. Lalanne, * !* P. Berthier, J. Der Hagopian, Masson, Paris 1980" [BIBLI 16]. * !* * !* F90 Release By J-P Moreau, Paris. * !* (www.jpmoreau.fr) * !******************************************************************************* !NOTE: For a bending beam of constant section S, the deviation v(x) is given !by: ! EI d4v/dx4 + rho*S d2v/dt2 - Tex = 0 (1) !where: ! E = Young's modulus of beam material (steel: 2E11 Pa) ! I = Moment of inertia of S (rectangle: B*H^3/12) ! rho = Volumic Mass (steel: 7800 kg/m3) ! Tex = External force per length unit ! x = beam abscissa from x=0 to x=L ! ! d4v/dx4: 4th partial derivative with respect to x, ! d2v/dt2: 2nd partial derivative with respect to time, t. ! !We seek solutions of v(x,t) having the form: v(x) * f(t). This leads to solving ! ! d2f(t)/dt2 + w2 f(t) = 0 (2) !and ! d4v(x)/dx4 - rho.S/(EI) w2 v(x) = 0 (3) !with ! w = pulsation of force excitation ! !The solution of (2) has the form: f(t) = A sin(wt) + B cos(wt) (4) !The solution of (3) is obtained by considering v(x) = V e^rx (5) ! !the characreristic equation of which is: r^4 - rho.S w2 /(EI) = 0 (6) ! !The complex roots of (6) are: beta, -beta, j.beta and -j.beta, with ! ! beta = 4th Root of (rho.s w2 / (EI)) (7) ! !From (5) and the the 4 complex roots of (6), we can write: ! !v(x) = C sin(beta.x) + D cos(beta.x) + E sh(beta.x) + F ch(beta.x) (8) ! !The resonance pulsations are given by: wn = Xn^2 /L^2 * sqrt(EI/(rho.S)) (9) ! !The Xn values are presented below for various limit conditions: ! ! Limit Conditions X1^2 X2^2 X3^2 X4^2 X5^2 !----------------------------------------------------------------------- !Fixed-free: 1+ch(X)*cos(X)=0 3.516 22.03 61.69 120.9 199.8 !Sup.-Sup.: sin(X)=0 9.869 39.47 88.82 157.9 246.7 !Fixed-Fixed or Free-Free: !1 - ch(X) cos(X) = 0 22.37 61.67 120.9 199.8 298.5 !Fixed-Sup. or Free-Sup.: !tan(X) = tanh(X) 15.41 49.96 104.2 178.2 272.0 !----------------------------------------------------------------------- ! !The max. strain is given by: sigma = E * H/2 * d2vi(x) / dx2 ! !------------------------------------------------------------------------------- PROGRAM BEAM1 IMPLICIT REAL*8 A-H, O-Z DIMENSION A1(5),X(25) REAL*8 I1,I2,L,L0,J,K1,M1 CHARACTER*3 Ans PRINT *,' ' PRINT *,'Fixed-Free Beam, M=1' PRINT *,'Fixed-Supported, M=2' PRINT *,'Fixed-Fixed, M=3' PRINT *,'Free-Free, M=4' PRINT *,'Free-Supported, M=5' PRINT *,'Supported-Supported, M=6' PRINT *,' ' WRITE(*,10,advance='no'); READ *, M PI = 4.D0 * DATAN(1.D0) IF (M == 1) THEN A1(1) = 3.5160152D0 A1(2) = 22.034491D0 A1(3) = 61.697214D0 A1(4) = 120.90191D0 A1(5) = 199.85953D0 ELSE IF (M == 2.OR.M == 5) THEN A1(1) = 15.418205D0 A1(2) = 49.964862D0 A1(3) = 104.24769D0 A1(4) = 178.26972D0 A1(5) = 272.03097D0 ELSE IF (M == 3.OR.M == 4) THEN A1(1) = 22.373285D0 A1(2) = 61.672822D0 A1(3) = 120.90339D0 A1(4) = 199.85944D0 A1(5) = 298.55553D0 ELSE IF (M == 6) THEN A1(1) = 9.8696044D0 A1(2) = 39.478417D0 A1(3) = 88.826439D0 A1(4) = 157.91367D0 A1(5) = 246.74011D0 END IF PRINT *,' ' WRITE(*,20,advance='no'); READ *, E PRINT *,' ' WRITE(*,30,advance='no'); READ *, R PRINT *,' ' WRITE(*,40,advance='no'); READ *, B WRITE(*,50,advance='no'); READ *, H WRITE(*,60,advance='no'); READ *, L PRINT *,' ' J = B * H**3 / 12.D0 !Moment of Inertia S = B * H !Section DO I = 1, 5 B1 = DSQRT(A1(I)) F = A1(I) / (2.D0 * PI * L**2) * DSQRT(E * J / R / S) WRITE(*,190) F END DO PRINT *,' ' WRITE(*,70,advance='no'); READ *, Ans IF (Ans == 'n') STOP PRINT *,' ' WRITE(*,80,advance='no'); READ *, N1 PRINT *,' ' WRITE(*,90,advance='no'); READ *, N2 PRINT *,' ' WRITE(*,100,advance='no'); READ *, Ans IF (Ans == 'n') THEN DO I = 1, N2 WRITE(*,110,advance='no') I; READ *, X(I) END DO ELSE L0 = L / (N2 - 1) X(1) = 0.D0 DO I = 2, N2 X(I) = X(I - 1) + L0 END DO END IF DO I = 1, N1 B1 = DSQRT(A1(I)) B2 = B1 / L O1 = DCOS(B1); I1 = DSIN(B1) O2 = DCOS(2.D0 * B1); I2 = DSIN(2.D0 * B1) C1 = (DEXP(B1) + DEXP(-B1)) / 2.D0 S1 = (DEXP(B1) - DEXP(-B1)) / 2.D0 C2 = (DEXP(2.D0 * B1) + DEXP(-2.D0 * B1)) / 2.D0 S2 = (DEXP(2.D0 * B1) - DEXP(-2.D0 * B1)) / 2.D0 A = 1.D0 IF (M == 1) THEN B = -(I1 + S1) / (O1 + C1) C = -1.D0 D = -B ELSE IF (M == 2.OR.M == 3) THEN B = (-I1 + S1) / (O1 - C1) C = -1.D0 D = -B ELSE IF (M == 4) THEN B = (I1 - S1) / (-O1 + C1) C = 1.D0 D = B ELSE IF (M == 5) THEN B = -(I1 + S1) / (O1 + C1) C = 1.D0 D = B ELSE IF (M == 6) THEN B = 0.D0 C = 0.D0 D = 0.D0 END IF T1 = A ** 2 / 2 * (L - 1 / 2.D0 / B2 * I2) T2 = B ** 2 / 2 * (L + I2 / 2.D0 / B2) T3 = C ** 2 / 2 * (S2 / 2.D0 / B2 - L) T4 = D ** 2 / 2 * (S2 / 2.D0 / B2 + L) T5 = A * B / 2.D0 / 2.D0 / B2 * (1.D0 - O2) T6 = C * D / 2.D0 / 2.D0 / B2 * (C2 - 1.D0) T7 = DEXP(B1) / 2.D0 / 2.D0 / B2 * (A * C + A * D) * (I1 - O1) T7 = T7 + DEXP(-B1) / 2.D0 / 2.D0 / B2 * (A * D - A * C) * (I1 + O1) T7 = T7 + A * C / 2.D0 / B2 T8 = DEXP(B1) / 2.D0 / 2.D0 / B2 * (B * C + B * D) * (I1 + O1) T8 = T8 + DEXP(-B1) / 2.D0 / 2.D0 / B2 * (B * D - B * C) * (I1 - O1) T8 = T8 - B * C / 2.D0 / B2 T = T1 + T2 + T3 + T4 + T5 + T6 + T7 + T8 M1 = R * S * T K1 = M1 * B2 ** 4 * E * H ** 2 / 12.D0 / R W = DSQRT(K1 / M1) F = W / 2.D0 / PI PRINT *,' ' WRITE(*,200) I PRINT *,' ----------' PRINT *,' ' WRITE(*,210) W, F PRINT *,' ' PRINT *,' Modal Mass=', M1, ' Modal Stiffness=', K1 PRINT *,' ' PRINT *,' X Deviation Slope Strain (x 10^6) ' PRINT *,' -------------------------------------------------' DO K = 1, N2 Z = X(K) I1 = DSIN(B2 * Z); O1 = DCOS(B2 * Z) S1 = (DEXP(B2 * Z) - DEXP(-B2 * Z)) / 2.D0 C1 = (DEXP(B2 * Z) + DEXP(-B2 * Z)) / 2.D0 Y = A * I1 + B * O1 + C * S1 + D * C1 P = B2 * (A * O1 - B * I1 + C * C1 + D * S1) M1 = B2**2 * (-A * I1 - B * O1 + C * S1 + D * C1) M1 = E * H * M1 / 2.D0 WRITE(*,220) Z, Y, P, M1 END DO print *,' ' Pause END DO 10 FORMAT(' Limit Conditions, M= ') 20 FORMAT(' Young''s Modulus: ') 30 FORMAT(' Volumic Mass...: ') 40 FORMAT(' Beam Width....: ') 50 FORMAT(' Beam Thickness: ') 60 FORMAT(' Beam Length...: ') 70 FORMAT(' Do you want the modes, modal masses & Stiffnesses (y/n): ') 80 FORMAT(' How many modes (Maximum 5): ') 90 FORMAT(' How many points for deviation, slope & max. strain: ') 100 FORMAT(' Do you want automatic divisions (y/n): ') 110 FORMAT(' X(',I2,') = ') 190 FORMAT(' Frequency (Hz) = ',F11.6) 200 FORMAT(' MODE #',I1) 210 FORMAT(' OMEGA=', F9.4, ' FREQUENCY=', F9.4) 220 FORMAT(F10.4,F10.4,F10.4,F12.2) END !end of file beam1.f90