Structural Dynamics¶

CalcpadCE worksheets in this section formulate structural dynamics problems as multi-degree-of-freedom systems with a flexibility-based mass-matrix formulation. The natural frequencies and mode shapes follow from the generalised eigenvalue problem \(\mathbf D \mathbf M \boldsymbol\phi = \lambda \boldsymbol\phi\), where \(\mathbf D\) is built by integrating bending and shear unit-load diagrams along the member and \(\mathbf M\) is the lumped mass matrix at the discretisation nodes. Time-history responses are then obtained from the Duhamel integral over each modal coordinate, with Rayleigh damping calibrated on the first two modes.

The steel pole free-vibration sheet treats a tapered cantilever of variable cross-section — rectangular, circular, elliptical, elliptical tube, trapezoidal, T or double-T — and integrates the area, second moment and shear area along the height before assembling the flexibility matrix and extracting the lowest natural frequencies, which are then compared with the empirical formula in ASCE/SEI 7-22. An animated SVG variant and a Plotly variant replay the mode shapes in the browser without re-running the numerical solution.

The beam impact sheet couples a steel ball drop with a simply supported reinforced concrete beam: the perfectly inelastic collision sets the post-impact velocity from momentum conservation, the contact duration follows the Hertzian impact contact-duration formula, and the time-history of the mid-span deflection is computed both as a single-degree-of-freedom approximation and as a full 11-node modal-superposition response with mass at every joint.

Beam Impact Analysis 🎬¶

Drop-weight impact of a steel ball on a simply supported reinforced concrete beam: SDOF and 11-node MDOF modal superposition with Rayleigh damping, Hertzian impact contact duration and Duhamel integration of the impulse force.

Code:

'Ball mass -'M_s = 2.1045t
'Ball material - steel
'Modulus of elasticity -'E_s = 206GPa
'Poisson′s ratio - 'ν_s = 0.3
'Mass density -'ρ_s = 7.85t/m^3
'Ball volume -'V_s = M_s/ρ_s' = 4/3πR³
'Ball radius -'R_s = cbrt(3/4*(V_s/π))|mm
'Height of bottom above the beam surface -'H = 2m
'Structure type - simply supported beam
'Beam length -'L = 12m
'Material - reinforced concrete C20/25
'Modulus of elasticity -'E = 20GPa
'Poisson′s ratio -'ν = 0.2
'Shear modulus -'G = E/(2*(1 + ν))
'Unit weight -'γ_b = 25kN/m^3
'Cross section - rectangular with dimensions:
'Width -'b = 350mm
'Height -'h = 650mm
'Area -'A = b*h|cm^2
'Second moment of area -'I = b*h^3/12
'Shear area -'A_Q = 5/6*A
'Self-weight -'g_b = A*γ_b
'Uniform load -'q = 10kN/m
'Gravity acceleration -'g = 9.80665m/s^2
'Uniform mass -'m = (g_b + q)/g
#hide
PlotSVG = 1
w = L
h_ = H + R_s + h + 0.5m
W = 250
H_ = h_*W/w
k = W/w
x1 = 0.5
y1 = (H + h)*k + 0.5
#def svg$ = '<svg viewbox="-10 '-R_s*k-20' 'w*k + 20' 'h_*k + 40'" xmlns="http://www.w3.org/2000/svg" version="1.1" style="font-family: Georgia Pro; font-size:9px; width:'W + 150'pt; height:'H_ + 250*H_/W'pt;">
#def thin_style$ = style="stroke:green; stroke-width:1; fill:none"
#def thick_style$ = style="stroke:green; stroke-width:2; fill:none"
#def ball_style$ = style="stroke:steelblue; stroke-width:0.5;  fill:url(#ball)"
#def solid_style$ = style="stroke:#888; stroke-width:1; fill:url(#concrete); fill-opacity: 0.5"
#def load_style$ = style="stroke:steelblue; stroke-width:0.5; fill:steelblue;"
#include svg_drawing.cpd
#show
#val
svg$
'<defs><radialGradient id="ball" cx="35%" cy="35%"><stop offset="0%" stop-color="lightcyan"/><stop offset="100%" stop-color="lightblue"/></radialGradient><pattern id="concrete" x="4" y="4" width="8" height="8" patternUnits="userSpaceOnUse"><rect x="0" y="0" width="8" height="8" fill="#eae9e8" /><circle cx="1" cy="1" r="1.2" fill="#ccb" /><circle cx="5" cy="2" r="1.6" fill="#dadad0" /><circle cx="4" cy="6" r="0.8" fill="#aa9" /><circle cx="3" cy="4" r="0.4" fill="#884" /><circle cx="7" cy="5" r="1.2" fill="#cacaba" /><circle cx="5" cy="3" r="0.9" fill="#fffded" /></pattern></defs>
circle$(L*k/2; -R_s*k; R_s*k; ball_style$)
rect$(0; H*k; L*k; h*k; solid_style$)
dimh$(0; L*k; (H+h+1m)*k;L='L'm)
dimv$(0.5*(L-1m)*k; 0; H*k;H='H'm)
texth$(0.5*(L+2*R_s+2.5m)*k; -R_s*k/2;M='M_s't)
#for i = 1 : 2
    #hide
    δ = w/30*k*sign(x1 - w/2*k)
    x2 = x1 - δ
    y2 = y1 - abs(δ)
    x3 = x1 + δ
    y3 = y1 + abs(δ)
    #show
    line$(x2; y3; x3; y3; thick_style$)
    line$(x2; y3; x1; y1; thin_style$)
    line$(x3; y3; x1; y1; thin_style$)
    #hide
    x1 = L*k - 0.5
    #show
#loop
#hide
δ = 10
δ_x = 0.2*δ
δ_y = -1.5*δ
x1 = L/2*k
y1 = -δ_y
#show
'<polygon points="'x1','y1 + δ_y' 'x1','y1' 'x1 - δ_x','y1 + δ_y/3' 'x1 + δ_x','y1 + δ_y/3' 'x1','y1'" load_style$ />
texth$(x1+4*δ_x;y1;g)
'</svg>
#equ
'<h4>Simple analytical solution</h4>
'The structure is reduced to a SDOF system for simplicity
'Dynamically equivalent mass -'M_d = 2*L/π*m
'Potential energy of the ball before dropping
E_p = M_s*g*H|kJ
#noc
'Kinetic energy immediately before the impact -'E_k = M_s*v_0^2/2
'The velocity at the moment before the impact is determined by the energy conservation law'E_k = E_p':
#equ
v_0 = sqrt(2*E_p/M_s)
'Perfectly inelastic collision model is assumed.
'Total mass after contact -'M_tot = M_s + M_d
'The velocity immediately after the contact is determined by the law of conservation of momentum:
v_1 = v_0*M_s/M_tot
'Structural stiffness for a vertical force applied at the middle point of the span
K = 48*E*I/L^3
'Deflection due to uniform load
z_0 = 5*(g_b + q)*L^4/(384*E*I)|mm
'Static displacement -'z_st = M_tot*g/K|mm
'Natural circular frequency -'ω_1 = sqrt(K/M_tot)
'Vibration period -'T_1 = 2*π/ω_1
'Dynamic factor
μ = 1 + sqrt(1 + (v_1*ω_1/g)^2)
'Dynamic displacement -'z_d = μ*z_st
'Dynamic force -'F_d = μ*M_s*g
'(without self-weight and uniform load)
'Simplified equation for the dynamic factor
μ_1 = 1 + v_1*ω_1/g
'The difference will be smaller for greater heights.
'Time before impact -'t_0 = sqrt(2*H/g)
'<h4>Elastic time history response of the structure as an SDOF system</h4>
'Damped vibration is assumed with factor -'ξ = 0.05
'Vibration amplitude -'A = z_d - z_st' or
A = v_1/ω_1|mm
#rad
'Theoretical equation of motion
y(t) = A*e^(-ξ*ω_1*t)*sin(ω_1*t)
'Solution by direct integration of the impulse load
'Duration of impulse transmission for a beam with infinite mass [1]
τ_L = 2.94*root((15/16*M_s*((1 - ν^2)/E + (1 - ν_s^2)/E_s))^2/(R_s*v_0); 5)|ms
'Duration of impulse transmission for a beam with finite mass  [2]
τ_L = 2.94*root((15/16*(M_s*M_d/(M_s + M_d))*((1 - ν^2)/E + (1 - ν_s^2)/E_s))^2/(R_s*v_0); 5)|ms
'The above values correspond well to the experimental data in [3], where the recorded durations are of a similar magnitude.
'The impulse force function will be determined by using the recommended expressions (9.20) - (9.22) in [1]
'The coefficient of restitution for perfectly inelastic collision is -'e. = 0
F(t) = M_s*v_0*(1 + e.)*(π/(2*τ_L))*sin(π/τ_L*t)*(abs(t) ≤ τ_L)|kN
'Impulse load diagram
'<!--'PlotWidth = 500','PlotHeight = 100'-->
$Plot{F(t) @ t = 0ms : 50ms}
'Maximal impulse load value -'F_max = F(τ_L/2)
'The equation of motion is expressed by the Duhamel′s integral
y_D(t) = 1/(M_tot*ω_1)*$Area{F(τ)*e^(-ξ*ω_1*(t - τ))*sin(ω_1*(t - τ)) @ τ = 0ms : min(t; τ_L)}|mm
'Static displacement for the midpoint of the beam
y_0(t) = z_0 + (z_st - z_0)*if(t < T_1/4; sin(2*π*t/T_1); 1)
'Time history of the midpoint displacement, [mm]
'<!--'PlotHeight = 200'-->
$Plot{t|-y_0(t) - y(t) & t|-y_0(t) - y_D(t) @ t = 0ms : 5s}
'<h4>Elastic time history response of the structure as an MDOF system</h4>
'Number of intermediate joints -'n_J = 11'(odd)
'Length of one segment -'Δx = L/(n_J + 1)
'Coordinate of joint <var>j</var> -'x_J(j) = Δx*j
'Shear forces due to unit vertical load <var>F</var><sub>j</sub> = 1 at joint <var>j</var>
V_1(x; j) = if(x < x_J(j); 1 - x_J(j)/L; -x_J(j)/L)
'<!--'PlotHeight = 100'-->
$Plot{V_1(x; 1) & V_1(x; 2) & V_1(x; 3) @ x = 0m : L}
'Bending moments due to unit vertical load <var>F</var><sub>j</sub> = 1 at joint <var>j</var>
M_1,max(j) = (x_J(j)/L - 1)*x_J(j)
M_1(x; j) = M_1,max(j)*if(x < x_J(j); x/x_J(j); (L - x)/(L - x_J(j)))
$Plot{-M_1(x; 1) & -M_1(x; 2) & -M_1(x; 3) @ x = 0m : L}
'Flexibility matrix
D(i; j) = $Area{M_1(x; i)*M_1(x; j) @ x = 0m : L}*(1/(E*I)) + $Area{V_1(x; i)*V_1(x; j) @ x = 0m : L}*(1/(G*A_Q))
#hide
D = symmetric(n_J)
#for i = 1 : n_J
    #for j = 1 : n_J
        x_min = min(x_J(i); x_J(j))
        D.(i; j) = D(i; j)*kN/mm
    #loop
#loop
#show
D'mm/kN
'Mass matrix
d_M.j = m*Δx/t
#hide
d_M = vector(n_J)
#for j = 1 : n_J
    d_M.j = m*Δx/t
#loop
'Middle joint number -'j_m = (n_J + 1)/2
d_M.(j_m) = d_M.(j_m) + M_s/t
M = vec2diag(d_M)
#def n$ = 7
#show
M
'Total mass of the structure -'sum(d_M)'t
'Eigenvalues
M_sq = sqrt(M)
C = copy(M_sq*D*M_sq; symmetric(n_J); 1; 1)
λ = eigenvals(C*10^-3; -n$)
'Natural circular frequences -'ω = sqrt(1/λ)'<i>s</i>⁻¹
'Natural vibration frequences -'f = ω/(2*π)*Hz
'Natural vibration periods -'T = 1/f
'Eigenvectors
V = transp(eigenvecs(C*10^-3; -n$))
Φ = inverse(M_sq)*V
X = stack(matrix(1; 3); Φ; matrix(1; 3))
#def X$(m$) = i*Δx|spline(i + 1; m$; X)
$Plot{X$(1) & X$(2) & X$(3) @ i = 0 : n_J + 1}
'Modal masses -'m_Φ = diag2vec(transp(Φ)*M*Φ)*t
'Rayleigh damping model is assumed
β = 2*ξ/(ω.1 + ω.2)', 'α = β*ω.1*ω.2
ξ(ω) = α/(2*ω) + β*ω/2
'Modal damping factors -'ξ_Φ = ξ(ω)
'<!--'PlotWidth = 300'-->
$Plot{ξ(ω_) & ω.1|ξ_Φ.1 & ω.2|ξ_Φ.2 & 0|0 & ω.3|ξ_Φ.3 & ω.4|ξ_Φ.4 & ω.5|ξ_Φ.5 @ ω_ = ω.1/5 : ω.n$}
'Damped natural frequences
ω_D = ω*sqrt(1 - ξ_Φ^2)*s^-1
'Dynamic load vector
F_Φ(i; t) = Φ.(j_m; i)*F(t)
'The equations of modal dynamic displacements are expressed by the Duhamel′s integral
y_Φ(i; t) = 1/(m_Φ.i*ω_D.i)*$Area{F_Φ(i; τ)*e^(-ξ_Φ.i*ω.i*s^-1*(t - τ))*sin(ω_D.i*(t - τ)) @ τ = 0ms : min(t; τ_L)}|mm
'Joint displacements
y_J(j; t) = $Sum{Φ.(j; i)*y_Φ(i; t) @ i = 1 : n$}
'<!--'PlotWidth = 500','PlotHeight = 200','PlotAdaptive = 0'-->
'Comparison of time history records of the midpoint displacements for SDOF and MDOF systems, [mm]
$Plot{t|-y_0(t) - y_J(j_m; t) & t|-y_0(t) - y_D(t) @ t = 0ms : 5s}
#def y$(j$) = t|-y_J(j$; t)
'Time history records for the amplitudes of separate joints, [mm] <!--Turn adaptive plot off for best results-->
$Plot{y$(1) & y$(2) & y$(3) & y$(4) & y$(5) @ t = 0ms : 5s}
#hide
Y = matrix(5; n_J + 2)*mm
Δt = 2*τ_L/5
t(i) = i*Δt
d(ξ) = 3.2*(ξ^4 - 2*ξ^3 + ξ)
$Repeat{$Repeat{Y.(i; j + 1) = d(j*Δx/L)*y_0(t(i)) + y_J(j; t(i)) @ j = 1 : n_J} @ i = 1 : 5}
#def Y$(i$) = j*Δx|-spline(j + 1; row(Y;i$))
'<!--'PlotHeight = 100','PlotAdaptive = 1'-->
n_T = 400
#show
#val
'<p>Beam deflections for the first five time steps at Δt = 'Δt' ms</p>
#equ
$Plot{Y$(1) & Y$(2) & Y$(3) & Y$(4) & Y$(5) @ j = 0 : n_J + 1}
#hide
Δt = 1.1*T.1/(n_T - 1)
t(i) = (i - 1)*Δt
Z = matrix(n_T; n_J + 2)*mm
$Repeat{$Repeat{Z.(i; j + 1) = d(j*Δx/L)*y_0(t(i)) + y_J(j; t(i)) @ j = 1 : n_J} @ i = 1 : n_T}
n_T = n_T/3
k_t = 20
t_0 = t_0/k_t
n_0 = ceiling(t_0/Δt)
t_0 = (k_t - 1)*t_0
N = n_0 + n_T
z_0(t) = H - g/2*(t_0 + t)^2
z_0 = z_0(0s)
PlotHeight = 100*(1 + z_0/100mm)
k_R = PlotWidth*mm/L
#show
'Animation of beam elastic response (slowed down)
#val
'<style>.fr{display:none;} circle{stroke-width: 0;} .fr circle.Series3{stroke-width:'6*R_s*k_R'px!Important;stroke-opacity:1; fill:none;}</style>
#for k = 1 : N
    #hide
    k_0 = max(k - n_0; 1)
    z = row(Z; k_0)
    j_m = ceiling(n_J/2)
    z_b(k) = if(k ≤ n_0; z_0((k - 1)*Δt); -spline(j_m; z)) + R_s*k_R
    #show
    '<div class="fr" id="beamimpact-fr-'k'">
    $Plot{j*Δx|-spline(j + 1; z) & 0m|-100mm & j_m*Δx|z_b(k) & L|z_0 + 50mm @ j = 0 : n_J + 1}
    '</div>
#loop
'<script>(function(){$("#beamimpact-fr-1").show();var i=1;var fr=$("#beamimpact-fr-1");setInterval(function(){fr.hide();if(++i>'N')i=1;fr=$("#beamimpact-fr-"+i);fr.show();}, 20);})();</script>
#equ
'[1] Harris C. M., Piersol A.G., HARRIS’ SHOCK AND VIBRATION HANDBOOK, Fifth Edition, McGraw-Hill 2002, ISBN 0-07-137081-1
'[2] Qing Peng, Xiaoming Liu, Yueguang Wei, Elastic impact of sphere on large plate, Journal of the Mechanics and Physics of Solids, Volume 156, 2021, 104604, ISSN 0022 - 5096, <a href="https://doi.org/10.1016/j.jmps.2021.104604">https://doi.org/10.1016/j.jmps.2021.104604</a>
'[3] Hong Hao and Thong M. Pham, Performance of RC Beams with or without FRP Strengthening Subjected to Impact Loading, Proceedings of the 2nd World Congress on Civil, Structural, and Environmental Engineering (CSEE’17) Barcelona, Spain – April 3– 4, 2017 ISSN:2371 - 5294 <a href="https://www.researchgate.net/publication/316618143_Performance_of_RC_Beams_with_or_without_FRP_Strengthening_Subjected_to_Impact_Loading">DOI:10.11159/icsenm17.1</a>
'[4] Gugan, D. “Inelastic collision and the Hertz theory of impact.” American Journal of Physics 68 (2000): 920-924., <a href="http://www.oxfordcroquet.com/tech/gugan/index.asp">http://www.oxfordcroquet.com/tech/gugan/index.asp</a>

Rendered Output:

Ball mass - M_s = 2.1 t

Ball material - steel

Modulus of elasticity - E_s = 206 GPa

Poisson′s ratio - ν_s = 0.3

Mass density - ρ_s = 7.85 tm³

Ball volume - V_s = M_sρ_s = 2.1 t7.85 t ∕ m³ = 0.268 m³ = 4/3πR³

Ball radius - R_s = 3 34 · V_sπ = 3 34 · 0.268 m³3.14 = 400 mm

Height of bottom above the beam surface - H = 2 m

Structure type - simply supported beam

Beam length - L = 12 m

Material - reinforced concrete C20/25

Modulus of elasticity - E = 20 GPa

Poisson′s ratio - ν = 0.2

Shear modulus - G = E2 · ( 1 + ν ) = 20 GPa2 · ( 1 + 0.2 ) = 8.33 GPa

Unit weight - γ_b = 25 kNm³

Cross section - rectangular with dimensions:

Width - b = 350 mm

Height - h = 650 mm

Area - A = b · h = 350 mm · 650 mm = 2275 cm²

Second moment of area - I = b · h³12 = 350 mm · ( 650 mm ) ³12 = 8009895833 mm⁴

Shear area - A_Q = 56 · A = 56 · 2275 cm² = 1895.83 cm²

Self-weight - g_b = A · γ_b = 2275 cm² · 25 kN ∕ m³ = 5.69 kN ∕ m

Uniform load - q = 10 kNm

Gravity acceleration - g = 9.81 ms²

Uniform mass - m = g_b + qg = 5.69 kN ∕ m + 10 kN ∕ m9.81 m ∕ s² = 1.6 t ∕ m

Simple analytical solution

The structure is reduced to a SDOF system for simplicity

Dynamically equivalent mass - M_d = 2 · Lπ · m = 2 · 12 m3.14 · 1.6 t ∕ m = 12.22 t

Potential energy of the ball before dropping

E_p = M_s · g · H = 2.1 t · 9.81 m ∕ s² · 2 m = 41.28 kJ

Kinetic energy immediately before the impact - E_k = M_s · v₀²2

The velocity at the moment before the impact is determined by the energy conservation law E_k = E_p :

v₀ = 2 · E_pM_s = 2 · 41.28 kJ2.1 t = 6.26 m ∕ s

Perfectly inelastic collision model is assumed.

Total mass after contact - M_tot = M_s + M_d = 2.1 t + 12.22 t = 14.33 t

The velocity immediately after the contact is determined by the law of conservation of momentum:

v₁ = v₀ · M_sM_tot = 6.26 m ∕ s · 2.1 t14.33 t = 0.92 m ∕ s

Structural stiffness for a vertical force applied at the middle point of the span

K = 48 · E · IL³ = 48 · 20 GPa · 8009895833 mm⁴ ( 12 m ) ³ = 4449.94 kN ∕ m

Deflection due to uniform load

z₀ = 5 · ( g_b + q ) · L⁴384 · E · I = 5 · ( 5.69 kN ∕ m + 10 kN ∕ m ) · ( 12 m ) ⁴384 · 20 GPa · 8009895833 mm⁴ = 26.44 mm

Static displacement - z_st = M_tot · gK = 14.33 t · 9.81 m ∕ s²4449.94 kN ∕ m = 31.57 mm

Natural circular frequency - ω₁ = KM_tot = 4449.94 kN ∕ m14.33 t = 17.62 s^-1

Vibration period - T₁ = 2 · πω₁ = 2 · 3.1417.62 s^-1 = 0.356 s

Dynamic factor

μ = 1 + 1 + (v₁ · ω₁g)² = 1 + 1 + (0.92 m ∕ s · 17.62 s^-19.81 m ∕ s²)² = 2.93

Dynamic displacement - z_d = μ · z_st = 2.93 · 31.57 mm = 92.58 mm

Dynamic force - F_d = μ · M_s · g = 2.93 · 2.1 t · 9.81 m ∕ s² = 60.52 kN

(without self-weight and uniform load)

Simplified equation for the dynamic factor

μ₁ = 1 + v₁ · ω₁g = 1 + 0.92 m ∕ s · 17.62 s^-19.81 m ∕ s² = 2.65

The difference will be smaller for greater heights.

Time before impact - t₀ = 2 · Hg = 2 · 2 m9.81 m ∕ s² = 0.639 s

Elastic time history response of the structure as an SDOF system

Damped vibration is assumed with factor - ξ = 0.05

Vibration amplitude - A = z_d − z_st = 92.58 mm − 31.57 mm = 61.01 mm or

A = v₁ω₁ = 0.92 m ∕ s17.62 s^-1 = 52.21 mm

Theoretical equation of motion

y ( t ) = A · e^{-ξ · ω₁ · t} · sin ( ω₁ · t )

Solution by direct integration of the impulse load

Duration of impulse transmission for a beam with infinite mass [1]

τ_L = 2.94 · 5 (1516 · M_s · (1 − ν²E + 1 − ν_s²E_s))²R_s · v₀ = 2.94 · 5 (1516 · 2.1 t · (1 − 0.2²20 GPa + 1 − 0.3²206 GPa))²400 mm · 6.26 m ∕ s = 3.93 ms

Duration of impulse transmission for a beam with finite mass [2]

τ_L = 2.94 · 5 (1516 · M_s · M_dM_s + M_d · (1 − ν²E + 1 − ν_s²E_s))²R_s · v₀ = 2.94 · 5 (1516 · 2.1 t · 12.22 t2.1 t + 12.22 t · (1 − 0.2²20 GPa + 1 − 0.3²206 GPa))²400 mm · 6.26 m ∕ s = 3.69 ms

The above values correspond well to the experimental data in [3], where the recorded durations are of a similar magnitude.

The impulse force function will be determined by using the recommended expressions (9.20) - (9.22) in [1]

The coefficient of restitution for perfectly inelastic collision is - e. = 0

F ( t ) = M_s · v₀ · ( 1 + e. ) · π2 · τ_L · sin(πτ_L · t) · ( | t | ≤ τ_L )

Impulse load diagram

Maximal impulse load value - F_max = F (τ_L2) = F (3.69 ms2) = 5613.66 kN

The equation of motion is expressed by the Duhamel′s integral

y_D ( t ) = 1M_tot · ω₁ · min ( t; τ_L ) ∫0 ms F ( τ ) · e^{-ξ · ω₁ · ( t − τ )} · sin ( ω₁ · ( t − τ ) ) dτ

Static displacement for the midpoint of the beam

y₀ ( t ) = z₀ + ( z_st − z₀ ) · {if t < T₁4: sin(2 · π · tT₁)
else: 1

Time history of the midpoint displacement, [mm]

Elastic time history response of the structure as an MDOF system

Number of intermediate joints - n_J = 11 (odd)

Length of one segment - Δx = Ln_J + 1 = 12 m11 + 1 = 1 m

Coordinate of joint j - x_J ( j ) = Δx · j

Shear forces due to unit vertical load F_j = 1 at joint j

V₁ ( x; j ) = {if x < x_J ( j ) : 1 − x_J ( j ) L
else: -x_J ( j ) L

Bending moments due to unit vertical load F_j = 1 at joint j

M_1,max ( j ) = (x_J ( j ) L − 1) · x_J ( j )

M₁ ( x; j ) = M_1,max ( j ) · {if x < x_J ( j ) : xx_J ( j )
else: L − xL − x_J ( j )

Flexibility matrix

D ( i; j ) = ( L∫0 m M₁ ( x; i ) · M₁ ( x; j ) dx) · 1E · I + ( L∫0 m V₁ ( x; i ) · V₁ ( x; j ) dx) · 1G · A_Q

D = 0.02160.03780.04890.05520.05740.0560.05140.04430.0350.02420.0124 0.03780.07040.0930.1060.1110.1090.10.08640.06850.04740.0242 0.04890.0930.1280.1490.1580.1550.1440.1240.09880.06850.035 0.05520.1060.1490.1790.1930.1930.180.1560.1240.08640.0443 0.05740.1110.1580.1930.2140.2170.2050.180.1440.10.0514 0.0560.1090.1550.1930.2170.2270.2170.1930.1550.1090.056 0.05140.10.1440.180.2050.2170.2140.1930.1580.1110.0574 0.04430.08640.1240.1560.180.1930.1930.1790.1490.1060.0552 0.0350.06850.09880.1240.1440.1550.1580.1490.1280.0930.0489 0.02420.04740.06850.08640.10.1090.1110.1060.0930.07040.0378 0.01240.02420.0350.04430.05140.0560.05740.05520.04890.03780.0216 mm/kN

Mass matrix

d_M.j = m · Δxt = 1.6 t ∕ m · 1 mt = 1.6

M = 1.60000000000 01.6000000000 001.600000000 0001.60000000 00001.6000000 000003.700000 0000001.60000 00000001.6000 000000001.600 0000000001.60 00000000001.6

Total mass of the structure - sum ( ⃗d_M ) = 19.7 t

Eigenvalues

M_sq = √ M = 1.260000000000 01.26000000000 001.2600000000 0001.260000000 00001.26000000 000001.9200000 0000001.260000 00000001.26000 000000001.2600 0000000001.260 00000000001.26

C = copy ( M_sq · D · M_sq; symmetric ( n_J ) ; 1; 1 ) = copy ( M_sq · D · M_sq; symmetric ( 11 ) ; 1; 1 ) = 0.03450.06050.07810.08830.09180.1360.08220.07080.0560.03870.0198 0.06050.1130.1490.170.1780.2650.160.1380.110.07580.0387 0.07810.1490.2040.2380.2520.3780.230.1990.1580.110.056 0.08830.170.2380.2870.3090.4690.2870.250.1990.1380.0708 0.09180.1780.2520.3090.3430.5290.3280.2870.230.160.0822 0.1360.2650.3780.4690.5290.8390.5290.4690.3780.2650.136 0.08220.160.230.2870.3280.5290.3430.3090.2520.1780.0918 0.07080.1380.1990.250.2870.4690.3090.2870.2380.170.0883 0.0560.110.1580.1990.230.3780.2520.2380.2040.1490.0781 0.03870.07580.110.1380.160.2650.1780.170.1490.1130.0605 0.01980.03870.0560.07080.08220.1360.09180.08830.07810.06050.0345

⃗λ = eigenvals ( C · 10^-3; -7 ) = [0.00261 0.000137 3.32×10^-5 9.33×10^-6 4.78×10^-6 2.17×10^-6 1.51×10^-6]

Natural circular frequences - ⃗ω = 1⃗λ = [19.57 85.55 173.53 327.32 457.58 678.76 814.79] s⁻¹

Natural vibration frequences - ⃗f = ⃗ω2 · π · Hz = ⃗ω2 · 3.14 · Hz = [3.11 Hz 13.61 Hz 27.62 Hz 52.09 Hz 72.83 Hz 108.03 Hz 129.68 Hz]

Natural vibration periods - ⃗T = 1⃗f = [0.321 s 0.0734 s 0.0362 s 0.0192 s 0.0137 s 0.00926 s 0.00771 s]

Eigenvectors

V = transp ( eigenvecs ( C · 10^-3; -7 ) ) = 0.09510.2040.2750.3540.380.4080.391 0.1840.3540.4030.3540.235-1.23×10^-12-0.143 0.2610.4080.3193.59×10^-13-0.237-0.408-0.337 0.320.3540.0714-0.354-0.3911.46×10^-130.273 0.3580.204-0.195-0.354-0.03420.4080.264 0.5640-0.484-1.44×10^-140.4249.64×10^-14-0.372 0.358-0.204-0.1950.354-0.0342-0.4080.264 0.32-0.3540.07140.354-0.391-2.73×10^-130.273 0.261-0.4080.319-3.5×10^-13-0.2370.408-0.337 0.184-0.3540.403-0.3540.2351.37×10^-12-0.143 0.0951-0.2040.275-0.3540.38-0.4080.391

Φ = inverse ( M_sq ) · V = 0.07520.1610.2170.280.3010.3230.309 0.1450.280.3190.280.186-9.74×10^-13-0.113 0.2060.3230.2522.84×10^-13-0.187-0.323-0.266 0.2530.280.0565-0.28-0.3091.15×10^-130.216 0.2830.161-0.155-0.28-0.0270.3230.209 0.2930-0.251-7.46×10^-150.225.01×10^-14-0.193 0.283-0.161-0.1550.28-0.027-0.3230.209 0.253-0.280.05650.28-0.309-2.16×10^-130.216 0.206-0.3230.252-2.77×10^-13-0.1870.323-0.266 0.145-0.280.319-0.280.1861.09×10^-12-0.113 0.0752-0.1610.217-0.280.301-0.3230.309

X = stack ( matrix ( 1; 3 ) ; Φ; matrix ( 1; 3 ) ) = 0000000 0.07520.1610.2170.280.3010.3230.309 0.1450.280.3190.280.186-9.74×10^-13-0.113 0.2060.3230.2522.84×10^-13-0.187-0.323-0.266 0.2530.280.0565-0.28-0.3091.15×10^-130.216 0.2830.161-0.155-0.28-0.0270.3230.209 0.2930-0.251-7.46×10^-150.225.01×10^-14-0.193 0.283-0.161-0.1550.28-0.027-0.3230.209 0.253-0.280.05650.28-0.309-2.16×10^-130.216 0.206-0.3230.252-2.77×10^-13-0.1870.323-0.266 0.145-0.280.319-0.280.1861.09×10^-12-0.113 0.0752-0.1610.217-0.280.301-0.3230.309 0000000

Modal masses - ⃗m_Φ = diag2vec ( transp ( Φ ) · M · Φ ) · t = [1 t 1 t 1 t 1 t 1 t 1 t 1 t]

Rayleigh damping model is assumed

β = 2 · ξ⃗ω₁ + ⃗ω₂ = 2 · 0.0519.57 + 85.55 = 0.000951 , α = β · ⃗ω₁ · ⃗ω₂ = 0.000951 · 19.57 · 85.55 = 1.59

ξ ( ω ) = α2 · ω + β · ω2

Modal damping factors - ⃗ξ_Φ = ξ ( ⃗ω ) = [0.05 0.05 0.0871 0.158 0.219 0.324 0.389]

Damped natural frequences

⃗ω_D = ⃗ω · √ 1 − ⃗ξ_Φ^⊙ 2 · s^-1 = [19.54 s^-1 85.44 s^-1 172.87 s^-1 323.2 s^-1 446.43 s^-1 642.14 s^-1 750.77 s^-1]

Dynamic load vector

F_Φ ( i; t ) = Φ_{j_m,i} · F ( t )

The equations of modal dynamic displacements are expressed by the Duhamel′s integral

y_Φ ( i; t ) = 1⃗m_Φ.i · ω_D.i · min ( t; τ_L ) ∫0 ms F_Φ ( i; τ ) · e^{-⃗ξ_Φ.i · ⃗ω_i · s^-1 · ( t − τ )} · sin ( ω_D.i · ( t − τ ) ) dτ

Joint displacements

y_J ( j; t ) = 7∑i= 1Φ_j,i · y_Φ ( i; t )

Comparison of time history records of the midpoint displacements for SDOF and MDOF systems, [mm]

Time history records for the amplitudes of separate joints, [mm]

Beam deflections for the first five time steps at Δt = 1.48 ms

Animation of beam elastic response (slowed down)

[1] Harris C. M., Piersol A.G., HARRIS’ SHOCK AND VIBRATION HANDBOOK, Fifth Edition, McGraw-Hill 2002, ISBN 0-07-137081-1

[2] Qing Peng, Xiaoming Liu, Yueguang Wei, Elastic impact of sphere on large plate, Journal of the Mechanics and Physics of Solids, Volume 156, 2021, 104604, ISSN 0022 - 5096, https://doi.org/10.1016/j.jmps.2021.104604

[3] Hong Hao and Thong M. Pham, Performance of RC Beams with or without FRP Strengthening Subjected to Impact Loading, Proceedings of the 2nd World Congress on Civil, Structural, and Environmental Engineering (CSEE’17) Barcelona, Spain – April 3– 4, 2017 ISSN:2371 - 5294 DOI:10.11159/icsenm17.1

[4] Gugan, D. “Inelastic collision and the Hertz theory of impact.” American Journal of Physics 68 (2000): 920-924., http://www.oxfordcroquet.com/tech/gugan/index.asp

Free Vibrations of Steel Pole 🎬¶

Same tapered steel pole as the static variant, with the lower mode shapes animated in pure SVG by frame switching driven by a small inline jQuery script.

Code:

"Free vibrations of steel pole with variable cross-section of arbitrary shape
'<hr/>
#rad
'<h3>Static scheme</h3>
'Cantilever with length -'L = 110m
'<h3>Material</h3>
'<h4>Steel</h4>
'Elastic modulus -'E = 206GPa
'Poisson′s ratio - 'ν = 0.3
'Shear modulus -'G = E/(2*(1 + ν))
'Mass density -'γ_s = 7.85*t/m^3
'<!--Section shape 'shape = 3', 'Precision = 10^-4'-->
'<h3>Cross-section</h3>
'Section height:
'- at bottom -'h_b = 3000mm
'- at top -'h_t = 750mm
'- difference -'Δh = h_t - h_b
'- as a function of distance to base -'h(x) = h_b + Δh*x/L
'Section width:
'- at bottom -'w_b = 3000mm
'- at top -'w_t = 750mm
'- difference -'Δw = w_t - w_b
'- as a function of distance to base -'w(x) = w_b + Δw*x/L
'<p>Cross-section shape -
#if shape ≡ 0
    '<b>Rectangular</b></p>
#else if shape ≡ 1
    '<b>Circular</b>
#else if shape ≡ 2
    '<b>Elliptical</b></p>
#else if shape ≡ 3
    '<b>Elliptical tube</b></p>
    'Thickness -'t_b = 20mm', 't_t = 9mm
#else if shape ≡ 4
    '<b>Trapezoidal</b></p>
    'Width at top edge -'b_2 = 90mm
#else if shape ≡ 5
    '<b>T-profile</b></p>
    'Top flange -'b_2 = 70mm', 'h_2 = 20mm
#else if shape ≡ 6
    '<b>Double-T profile</b></p>
    'Bottom flange -'b_1 = 70mm', 'h_1 = 20mm
    'Top flange -'b_2 = 90mm', 'h_2 = 15mm
#end if
'Function of cross-section outline
#if shape ≡ 0
    b(x; z) = w(x)
#else if shape ≡ 1
    b(x; z) = 2*sqrt((h(x)/2)^2 - (z - h(x)/2)^2)
#else if shape ≡ 2
    b(x; z) = 2*w(x)/h(x)*sqrt((h(x)/2)^2 - (z - h(x)/2)^2)
#else if shape ≡ 3
    t(x) = t_b + (t_t - t_b)*x/L
    b_i(x; z) = w(x) - 2*t(x)
    h_i(x; z) = h(x) - 2*t(x)
    b_1(x; z) = 2*w(x)/h(x)*sqrt((h(x)/2)^2 - (z - h(x)/2)^2)
    b_2(x; z) = if((z > t(x))*(z < (h(x) - t(x))); 2*b_i(x; z)/h_i(x; z)*sqrt((h_i(x; z)/2)^2 - (z - h(x)/2)^2); 0m)
    b(x; z) = b_1(x; z) - b_2(x; z)
#else if shape ≡ 4
    b(x; z) = w(z) + z*(b_2 - w(x))/h(x)
#else if shape ≡ 5
    b(x; z) = if(z > h(x) - h_2; b_2; w(x))
#else if shape ≡ 6
    b(x; z) = if(z > h(x) - h_2; b_2; if(z < h_1; b_1; w(x)))
#end if
'<h4>Cross section properties</h4>
δ = 0.01mm
#if shape ≡ 3
    'Area -'A(x) = $Integral{b_1(x; z) @ z = 0mm : h(x)} - $Integral{b_2(x; z) @ z = t(x) : h(x) - t(x)}|cm^2
    'First moment of area -'S(x) = $Integral{b_1(x; z)*z @ z = 0mm : h(x)} - $Integral{b_2(x; z)*z @ z = t(x) : h(x) - t(x)}|cm^3
    'Geometrical center-'z_c(x) = S(x)/A(x)|mm
    'Second moment of area
    I_1(x) = $Integral{b_1(x; z)*(z - z_c(x))^2 @ z = 0mm : h(x)}
    I_2(x) = $Integral{b_2(x; z)*(z - z_c(x))^2 @ z = t(x) : h(x) - t(x)}|cm^4
    I(x) = I_1(x) - I_2(x)
    'First moment of area under z
    S_1(x; z) = $Integral{b_1(x; ζ)*(z_c(x) - ζ) @ ζ = δ : z}
    S_2(x; z) = $Integral{b_2(x; ζ)*(z_c(x) - ζ) @ ζ = t(x) + δ : z}
    'Shear area
    A_Q1(x) = I_1(x)^2/$Integral{S_1(x; z)^2/b_1(x; z) @ z = δ : h(x) - δ}
    A_Q2(x) = I_2(x)^2/$Integral{S_2(x; z)^2/b_2(x; z) @ z = t(x) + δ : h(x) - t(x) - δ}
    A_Q(x) = A_Q1(x) - A_Q2(x)
#else
    'Area -'A(x) = $Area{b(x; z) @ z = 0mm : h(x)}|cm^2
    'First moment of area -'S(x) = $Area{b(x; z)*z @ z = 0mm : h(x)}|cm^3
    'Center of gravity -'z_c(x) = S(x)/A(x)|mm
    'Second moment of area -'I(x) = $Area{b(x; z)*(z - z_c(x))^2 @ z = 0mm : h(x)}|cm^4
    'First moment of area under z -'S_1(x; z) = $Area{b(x; ζ)*(z_c(x) - ζ) @ ζ = δ : z}
    'Shear area -'A_Q(x) = I(x)^2/$Area{S_1(x; z)^2/b(x; z) @ z = δ : h(x) - δ}
#end if
'<!--'PlotWidth = 200','PlotHeight = 200','PlotSVG = 1'-->
#def PlotSection$(x$)
    #if shape ≡ 3
        $Plot{0mm|z_c(x$) & b_1(x$; z)/2|z & -b_1(x$; z)/2|z & b_2(x$; z)/2|z & -b_2(x$; z)/2|z & h(x$)/2 - z|z_c(x$) @ z = 0mm : h(x$)}
    #else
        $Plot{0mm|z_c(x$) & b(x$; z)/2|z & -b(x$; z)/2|z & h(x$)/2 - z|z_c(x$) @ z = 0mm : x$}
    #end if
#end def
PlotSection$(0mm)
PlotSection$(L)
'<h3>Mass</h3>
'Distributed mass along height -'m(x) = A(x)*γ_s
'<h3>Solution</h3>
'Number of nodes -'n_J = 11
'Length of one segment -'Δx = L/n_J
'Elevation of node <var>j</var> -'x(j) = Δx*j
'Bending due to horizontal force <var>F</var><sub>j</sub> = 1 at node <var>j</var>
M(j; x) = max(x(j) - x; 0m)
'Flexibility matrix
D(i; j) = $Integral{M(i; x)*M(j; x)/I(x) @ x = 0m : L}*(1/E) + $Integral{1/A_Q(x) @ x = 0m : L}*(1/G)
#hide
D = symmetric(n_J)
#for i = 1 : n_J
    #for j = 1 : n_J
        x_min = min(x(i); x(j))
        D.(i; j) = ($Integral{M(i; x)*M(j; x)/I(x) @ x = 0m : x_min}/E + $Integral{1/A_Q(x) @ x = 0m : x_min}/G)*kN/mm
    #loop
#loop
#show
D'mm/kN
'Mass matrix
d_M.j = $Integral{m(x) @ x = x(j) - Δx/2 : x(j) + Δx/2}
#hide
d_M = vector(n_J)
#for j = 1 : n_J - 1
    d_M.j = $Integral{m(x) @ x = x(j) - Δx/2 : x(j) + Δx/2}/t
#loop
d_M.n_J = $Integral{m(x) @ x = L - Δx/2 : L}/t
M = vec2diag(d_M)
#show
M
'Total mass of structure -'M_tot = sum(d_M)'t
'Eigenvalues
M_sq = sqrt(M)
C = copy(M_sq*D*M_sq; symmetric(n_J); 1; 1)
λ = eigenvals(C*10^-3; -3)
'Natural frequencies
ω = sqrt(1/λ)
'Vibration frequencies
f = ω/(2*π)*Hz
'Vibration periods
T = 1/f
'Mode shapes
V = transp(eigenvecs(C*10^-3; -3))
Φ = inverse(M_sq)*V
X = stack(matrix(1; 3); Φ)
#hide
PlotHeight = 400
PlotWidth = 120
n_T = [80; 60; 30]
#show
$Plot{spline(i + 1; 1; X)|i*Δx & spline(i + 1; 2; X)|i*Δx & spline(i + 1; 3; X)|i*Δx @ i = 0 : n_J}
'Comparison with ASCE SEI 7/22 <span class="ref">(C26.11-13)</span>
λ_1 = 1.9*exp(-4*(h_t/h_b)) + 6.65/(0.9 + (t_t/t_b)^0.666)
'<p class="ref">(C26.11-12)</p>
'Fundamental natural frequency
n_1 = λ_1/(2*π*L^2)*sqrt(E*I(0m)/m(0m))|Hz
'Fundamental period of vibrations
T_1 = 1/n_1
'Animation for mode'i = 2
#val
#rad
'<style>[id^="steelpole-fr-"]{display:none;}</style>
#for k = 0 : n_T.i
    '<div id="steelpole-fr-'k'">
    $Plot{spline(j + 1; i; X)*cos((k/n_T.i*2*π))|j*Δx & -0.6|0m & 0.6|0m @ j = 0 : n_J}
    '</div>
#loop
'<script>var i='n_T.i';f=$("#steelpole-fr-0");setInterval(function(){f.hide();if(++i>'n_T.i')i=1;f=$("#steelpole-fr-"+i);f.show();},'1000/n_T.i*T.i');</script>

Rendered Output:

Free vibrations of steel pole with variable cross-section of arbitrary shape

Static scheme

Cantilever with length - L = 110 m

Material

Steel

Elastic modulus - E = 206 GPa

Poisson′s ratio - ν = 0.3

Shear modulus - G = E2 · ( 1 + ν ) = 206 GPa2 · ( 1 + 0.3 ) = 79.23 GPa

Mass density - γ_s = 7.85 tm³

Cross-section

Section height:

- at bottom - h_b = 3000 mm

- at top - h_t = 750 mm

- difference - Δh = h_t − h_b = 750 mm − 3000 mm = -2250 mm

- as a function of distance to base - h ( x ) = h_b + Δh · xL

Section width:

- at bottom - w_b = 3000 mm

- at top - w_t = 750 mm

- difference - Δw = w_t − w_b = 750 mm − 3000 mm = -2250 mm

- as a function of distance to base - w ( x ) = w_b + Δw · xL

Cross-section shape - Elliptical tube

Thickness - t_b = 20 mm , t_t = 9 mm

Function of cross-section outline

t ( x ) = t_b + ( t_t − t_b ) · xL

b_i ( x; z ) = w ( x ) − 2 · t ( x )

h_i ( x; z ) = h ( x ) − 2 · t ( x )

b₁ ( x; z ) = 2 · w ( x ) h ( x ) · (h ( x ) 2)² − (z − h ( x ) 2)²

b₂ ( x; z ) = {if ( z > t ( x ) ) · ( z < h ( x ) − t ( x ) ) : 2 · b_i ( x; z ) h_i ( x; z ) · (h_i ( x; z ) 2)² − (z − h ( x ) 2)²
else: 0 m

b ( x; z ) = b₁ ( x; z ) − b₂ ( x; z )

Cross section properties

δ = 0.01 mm

Area - A ( x ) = h ( x ) ∫0 mm b₁ ( x; z ) dz − h ( x ) − t ( x ) ∫t ( x ) b₂ ( x; z ) dz

First moment of area - S ( x ) = h ( x ) ∫0 mm b₁ ( x; z ) · z dz − h ( x ) − t ( x ) ∫t ( x ) b₂ ( x; z ) · z dz

Geometrical center- z_c ( x ) = S ( x ) A ( x )

Second moment of area

I₁ ( x ) = h ( x ) ∫0 mm b₁ ( x; z ) · ( z − z_c ( x ) ) ² dz

I₂ ( x ) = h ( x ) − t ( x ) ∫t ( x ) b₂ ( x; z ) · ( z − z_c ( x ) ) ² dz

I ( x ) = I₁ ( x ) − I₂ ( x )

First moment of area under z

S₁ ( x; z ) = z∫δ b₁ ( x; ζ ) · ( z_c ( x ) − ζ ) dζ

S₂ ( x; z ) = z∫t ( x ) + δ b₂ ( x; ζ ) · ( z_c ( x ) − ζ ) dζ

Shear area

A_Q1 ( x ) = I₁ ( x ) ² h ( x ) − δ∫δ S₁ ( x; z ) ²b₁ ( x; z ) dz

A_Q2 ( x ) = I₂ ( x ) ² h ( x ) − t ( x ) − δ∫t ( x ) + δ S₂ ( x; z ) ²b₂ ( x; z ) dz

A_Q ( x ) = A_Q1 ( x ) − A_Q2 ( x )

Mass

Distributed mass along height - m ( x ) = A ( x ) · γ_s

Solution

Number of nodes - n_J = 11

Length of one segment - Δx = Ln_J = 110 m11 = 10 m

Elevation of node j - x ( j ) = Δx · j

Bending due to horizontal force F_j = 1 at node j

M ( j; x ) = max ( x ( j ) − x; 0 m )

Flexibility matrix

D ( i; j ) = ( L∫0 m M ( i; x ) · M ( j; x ) I ( x ) dx) · 1E + ( L∫0 m 1A_Q ( x ) dx) · 1G

D = 0.009110.02190.03460.04740.06010.07290.08560.09840.1110.1240.137 0.02190.07310.1290.1850.2420.2980.3540.410.4660.5220.578 0.03460.1290.2630.4030.5430.6830.8240.9641.11.241.38 0.04740.1850.4030.6730.9521.231.511.792.072.352.63 0.06010.2420.5430.9521.431.932.432.923.423.914.41 0.07290.2980.6831.231.932.733.564.385.26.026.85 0.08560.3540.8241.512.433.564.846.157.478.7810.1 0.09840.410.9641.792.924.386.158.1710.2412.314.37 0.1110.4661.12.073.425.27.4710.2413.3916.6419.89 0.1240.5221.242.353.916.028.7812.316.6421.6726.88 0.1370.5781.382.634.416.8510.114.3719.8926.8835.23 mm/kN

Mass matrix

d_M.j = x ( j ) + Δx2∫x ( j ) − Δx2 m ( x ) dx = 1.06 t

M = 13.010000000000 011.43000000000 009.9400000000 0008.550000000 00007.26000000 000006.0800000 0000004.990000 00000004000 000000003.1200 0000000002.330 00000000000.904

Total mass of structure - M_tot = sum ( ⃗d_M ) = 71.62 t

Eigenvalues

M_sq = √ M = 3.610000000000 03.38000000000 003.1500000000 0002.920000000 00002.7000000 000002.4700000 0000002.230000 00000002000 000000001.7700 0000000001.530 00000000000.951

C = copy ( M_sq · D · M_sq; symmetric ( n_J ) ; 1; 1 ) = copy ( M_sq · D · M_sq; symmetric ( 11 ) ; 1; 1 ) = 0.1190.2670.3940.50.5850.6480.690.710.7080.6830.469 0.2670.8351.381.832.22.482.672.772.782.71.86 0.3941.382.613.714.615.315.86.086.155.994.15 0.51.833.715.757.518.889.8710.4810.6910.57.31 0.5852.24.617.5110.4212.8214.6115.7616.2716.1111.3 0.6482.485.318.8812.8216.619.5821.622.6422.6816.05 0.692.675.89.8714.6119.5824.1527.5129.4729.9821.45 0.712.776.0810.4815.7621.627.5132.7136.1737.627.33 0.7082.786.1510.6916.2722.6429.4736.1741.7744.8933.39 0.6832.75.9910.516.1122.6829.9837.644.8950.5839.03 0.4691.864.157.3111.316.0521.4527.3333.3939.0331.84

⃗λ = eigenvals ( C · 10^-3; -3 ) = [0.195 0.0173 0.00341]

Natural frequencies

⃗ω = 1⃗λ = [2.27 7.6 17.13]

Vibration frequencies

⃗f = ⃗ω2 · π · Hz = ⃗ω2 · 3.14 · Hz = [0.361 Hz 1.21 Hz 2.73 Hz]

Vibration periods

⃗T = 1⃗f = [2.77 s 0.826 s 0.367 s]

Mode shapes

V = transp ( eigenvecs ( C · 10^-3; -3 ) ) = 0.008960.03160.0687 0.03480.1120.217 0.07580.2190.357 0.130.3240.402 0.1940.3960.298 0.2640.4040.0531 0.3360.326-0.239 0.4030.151-0.414 0.459-0.104-0.312 0.495-0.3910.105 0.372-0.470.48

Φ = inverse ( M_sq ) · V = 0.002480.008750.019 0.01030.03320.0643 0.0240.06950.113 0.04430.1110.138 0.07180.1470.111 0.1070.1640.0216 0.150.146-0.107 0.2010.0754-0.207 0.26-0.0589-0.176 0.324-0.2560.069 0.391-0.4940.505

X = stack ( matrix ( 1; 3 ) ; Φ ) = 000 0.002480.008750.019 0.01030.03320.0643 0.0240.06950.113 0.04430.1110.138 0.07180.1470.111 0.1070.1640.0216 0.150.146-0.107 0.2010.0754-0.207 0.26-0.0589-0.176 0.324-0.2560.069 0.391-0.4940.505

Comparison with ASCE SEI 7/22 (C26.11-13)

λ₁ = 1.9 · exp(-4 · h_th_b) + 6.650.9 + (t_tt_b)^0.666 = 1.9 · exp(-4 · 750 mm3000 mm) + 6.650.9 + (9 mm20 mm)^0.666 = 5.17

(C26.11-12)

Fundamental natural frequency

n₁ = λ₁2 · π · L² · E · I ( 0 m ) m ( 0 m ) = 5.172 · 3.14 · ( 110 m ) ² · 206 GPa · I ( 0 m ) m ( 0 m ) = 0.367 Hz

Fundamental period of vibrations

T₁ = 1n₁ = 10.367 Hz = 2.72 s

Animation for mode i = 2

Free Vibrations of Steel Pole Plotly 🎬¶

Tapered steel pole free-vibration analysis with the mode shapes animated in the browser by Plotly, smoothly interpolating between the positive and negative amplitude of each mode.

Code:

#rad
'<h3>Static scheme</h3>
'Cantilever with length -'L = 110m
'<h3>Material</h3>
'<h4>Steel</h4>
'Elastic modulus -'E = 206GPa
'Poisson′s ratio - 'ν = 0.3
'Shear modulus -'G = E/(2*(1 + ν))
'Mass density -'γ_s = 7.85*t/m^3
'<!--Section shape 'shape = 3', 'Precision = 10^-4'-->
'<h3>Cross-section</h3>
'Section height:
'- at bottom -'h_b = 3000mm
'- at top -'h_t = 750mm
'- difference -'Δh = h_t - h_b
'- as a function of distance to base -'h(x) = h_b + Δh*x/L
'Section width:
'- at bottom -'w_b = 3000mm
'- at top -'w_t = 750mm
'- difference -'Δw = w_t - w_b
'- as a function of distance to base -'w(x) = w_b + Δw*x/L
'<p>Cross-section shape -
#if shape ≡ 0
    '<b>Rectangular</b></p>
#else if shape ≡ 1
    '<b>Circular</b>
#else if shape ≡ 2
    '<b>Elliptical</b></p>
#else if shape ≡ 3
    '<b>Elliptical tube</b></p>
    'Thickness -'t_b = 20mm', 't_t = 9mm
#else if shape ≡ 4
    '<b>Trapezoidal</b></p>
    'Width at top edge -'b_2 = 90mm
#else if shape ≡ 5
    '<b>T-profile</b></p>
    'Top flange -'b_2 = 70mm', 'h_2 = 20mm
#else if shape ≡ 6
    '<b>Double-T profile</b></p>
    'Bottom flange -'b_1 = 70mm', 'h_1 = 20mm
    'Top flange -'b_2 = 90mm', 'h_2 = 15mm
#end if
'Function of cross-section outline
#if shape ≡ 0
    b(x; z) = w(x)
#else if shape ≡ 1
    b(x; z) = 2*sqrt((h(x)/2)^2 - (z - h(x)/2)^2)
#else if shape ≡ 2
    b(x; z) = 2*w(x)/h(x)*sqrt((h(x)/2)^2 - (z - h(x)/2)^2)
#else if shape ≡ 3
    t(x) = t_b + (t_t - t_b)*x/L
    b_i(x; z) = w(x) - 2*t(x)
    h_i(x; z) = h(x) - 2*t(x)
    b_1(x; z) = 2*w(x)/h(x)*sqrt((h(x)/2)^2 - (z - h(x)/2)^2)
    b_2(x; z) = if((z > t(x))*(z < (h(x) - t(x))); 2*b_i(x; z)/h_i(x; z)*sqrt((h_i(x; z)/2)^2 - (z - h(x)/2)^2); 0m)
    b(x; z) = b_1(x; z) - b_2(x; z)
#else if shape ≡ 4
    b(x; z) = w(z) + z*(b_2 - w(x))/h(x)
#else if shape ≡ 5
    b(x; z) = if(z > h(x) - h_2; b_2; w(x))
#else if shape ≡ 6
    b(x; z) = if(z > h(x) - h_2; b_2; if(z < h_1; b_1; w(x)))
#end if
'<h4>Cross section properties</h4>
δ = 0.01mm
#if shape ≡ 3
    'Area -'A(x) = $Integral{b_1(x; z) @ z = 0mm : h(x)} - $Integral{b_2(x; z) @ z = t(x) : h(x) - t(x)}|cm^2
    'First moment of area -'S(x) = $Integral{b_1(x; z)*z @ z = 0mm : h(x)} - $Integral{b_2(x; z)*z @ z = t(x) : h(x) - t(x)}|cm^3
    'Geometrical center-'z_c(x) = S(x)/A(x)|mm
    'Second moment of area
    I_1(x) = $Integral{b_1(x; z)*(z - z_c(x))^2 @ z = 0mm : h(x)}
    I_2(x) = $Integral{b_2(x; z)*(z - z_c(x))^2 @ z = t(x) : h(x) - t(x)}|cm^4
    I(x) = I_1(x) - I_2(x)
    'First moment of area under z
    S_1(x; z) = $Integral{b_1(x; ζ)*(z_c(x) - ζ) @ ζ = δ : z}
    S_2(x; z) = $Integral{b_2(x; ζ)*(z_c(x) - ζ) @ ζ = t(x) + δ : z}
    'Shear area
    A_Q1(x) = I_1(x)^2/$Integral{S_1(x; z)^2/b_1(x; z) @ z = δ : h(x) - δ}
    A_Q2(x) = I_2(x)^2/$Integral{S_2(x; z)^2/b_2(x; z) @ z = t(x) + δ : h(x) - t(x) - δ}
    A_Q(x) = A_Q1(x) - A_Q2(x)
#else
    'Area -'A(x) = $Area{b(x; z) @ z = 0mm : h(x)}|cm^2
    'First moment of area -'S(x) = $Area{b(x; z)*z @ z = 0mm : h(x)}|cm^3
    'Center of gravity -'z_c(x) = S(x)/A(x)|mm
    'Second moment of area -'I(x) = $Area{b(x; z)*(z - z_c(x))^2 @ z = 0mm : h(x)}|cm^4
    'First moment of area under z -'S_1(x; z) = $Area{b(x; ζ)*(z_c(x) - ζ) @ ζ = δ : z}
    'Shear area -'A_Q(x) = I(x)^2/$Area{S_1(x; z)^2/b(x; z) @ z = δ : h(x) - δ}
#end if
'<!--'PlotWidth = 200','PlotHeight = 200'-->
#def PlotSection$(x$)
    #if shape ≡ 3
        $Plot{0mm|z_c(x$) & b_1(x$; z)/2|z & -b_1(x$; z)/2|z & b_2(x$; z)/2|z & -b_2(x$; z)/2|z & h(x$)/2 - z|z_c(x$) @ z = 0mm : h(x$)}
    #else
        $Plot{0mm|z_c(x$) & b(x$; z)/2|z & -b(x$; z)/2|z & h(x$)/2 - z|z_c(x$) @ z = 0mm : x$}
    #end if
#end def
PlotSection$(0mm)
PlotSection$(L)
'<h3>Mass</h3>
'Distributed mass along height -'m(x) = A(x)*γ_s
'<h3>Solution</h3>
'Number of nodes -'n_J = 11
'Length of one segment -'Δx = L/n_J
'Elevation of node <var>j</var> -'x(j) = Δx*j
'Bending due to horizontal force <var>F</var><sub>j</sub> = 1 at node <var>j</var>
M(j; x) = max(x(j) - x; 0m)
'Flexibility matrix
D(i; j) = $Integral{M(i; x)*M(j; x)/I(x) @ x = 0m : L}*(1/E) + $Integral{1/A_Q(x) @ x = 0m : L}*(1/G)
#hide
D = symmetric(n_J)
#for i = 1 : n_J
    #for j = 1 : n_J
        x_min = min(x(i); x(j))
        D.(i; j) = ($Integral{M(i; x)*M(j; x)/I(x) @ x = 0m : x_min}/E + $Integral{1/A_Q(x) @ x = 0m : x_min}/G)*kN/mm
    #loop
#loop
#show
D'mm/kN
'Mass matrix
d_M.j = $Integral{m(x) @ x = x(j) - Δx/2 : x(j) + Δx/2}
#hide
d_M = vector(n_J)
#for j = 1 : n_J - 1
    d_M.j = $Integral{m(x) @ x = x(j) - Δx/2 : x(j) + Δx/2}/t
#loop
d_M.n_J = $Integral{m(x) @ x = L - Δx/2 : L}/t
M = vec2diag(d_M)
#show
M
'Total mass of structure -'M_tot = sum(d_M)'t
'Eigenvalues
M_sq = sqrt(M)
C = copy(M_sq*D*M_sq; symmetric(n_J); 1; 1)
λ = eigenvals(C*10^-3; -3)
'Natural frequencies
ω = sqrt(1/λ)
'Vibration frequencies
f = ω/(2*π)*Hz
'Vibration periods
T = 1/f
'Mode shapes
V = transp(eigenvecs(C*10^-3; -3))
Φ = inverse(M_sq)*V
X = stack(matrix(1; 3); Φ)
#hide
PlotHeight = 400
PlotWidth = 120
#show
$Plot{spline(i + 1; 1; X)|i*Δx & spline(i + 1; 2; X)|i*Δx & spline(i + 1; 3; X)|i*Δx @ i = 0 : n_J}
'Comparison with ASCE SEI 7/22 <span class="ref">(C26.11-13)</span>
λ_1 = 1.9*exp(-4*(h_t/h_b)) + 6.65/(0.9 + (t_t/t_b)^0.666)
'<p class="ref">(C26.11-12)</p>
'Fundamental natural frequency
n_1 = λ_1/(2*π*L^2)*sqrt(E*I(0m)/m(0m))|Hz
'Fundamental period of vibrations
T_1 = 1/n_1
'<script src="https://cdnjs.cloudflare.com/ajax/libs/d3/3.5.17/d3.min.js"></script>
'<script src="https://cdn.plot.ly/plotly-2.26.0.min.js" charset="utf-8"></script>
#hide
y = range(1; n_J; 1)*L/n_J
#show
#def PlotlyAnimateMode$(n$)
    #val
    '<h4>Animation of mode n$ - T = 'T.n$'s</h4>
    '<div id="Moden$"></div>
    #hide
    x = col(X; n$)
    x_lim = round(norm_i(x)*10)/10
    #post
    '<script>
    'const yn$ = 'y';
    'var xn$ = 'x';
    'var tn$ = 'T.n$*500';
    'var layout = {autosize:false,width:300,height:600,xaxis:{range:[-'x_lim','x_lim'],tick0:-'x_lim',dtick: 0.1,ticklen:'2*x_lim'},yaxis:{dtick:10},margin:{l:25,r:10,b:50,t:10,pad:2},};
    'var data = [{x:xn$,y:yn$,line:{shape:"spline"},}];
    'Plotly.newPlot("Moden$",data,layout);
    'function animateModen$(){
    'for (let i = 0; i < xn$.length; i++) {xn$[i] = -xn$[i];}
    'Plotly.animate("Moden$",{data:[{x:xn$,y:yn$,line:{shape:"spline"},}],traces:[0],layout:layout},
    '{
    'transition:{duration:tn$,easing:"cos"},
    'frame:{duration:tn$, redraw: false},
    'mode: "next"
    '}).then(()=>requestAnimationFrame(animateModen$)).catch(()=>{});}
    'requestAnimationFrame(animateModen$);
    '</script>
    #equ
#end def
PlotlyAnimateMode$(1)

Free Vibrations of Steel Pole¶

Tapered cantilever pole with selectable cross-section shape; flexibility matrix built by integration of bending and shear along the height, three lowest natural frequencies, with the fundamental period checked against ASCE SEI 7/22.

Code:

#rad
'<h3>Static scheme</h3>
'Cantilever with length -'L = 110m
'<h3>Material</h3>
'<h4>Steel</h4>
'Elastic modulus -'E = 206GPa
'Poisson′s ratio - 'ν = 0.3
'Shear modulus -'G = E/(2*(1 + ν))
'Mass density -'γ_s = 7.85*t/m^3
'<!--Section shape 'shape = 3', 'Precision = 10^-4'-->
'<h3>Cross-section</h3>
'Section height:
'- at bottom -'h_b = 3000mm
'- at top -'h_t = 750mm
'- difference -'Δh = h_t - h_b
'- as a function of distance to base -'h(x) = h_b + Δh*x/L
'Section width:
'- at bottom -'w_b = 3000mm
'- at top -'w_t = 750mm
'- difference -'Δw = w_t - w_b
'- as a function of distance to base -'w(x) = w_b + Δw*x/L
'<p>Cross-section shape -
#if shape ≡ 0
    '<b>Rectangular</b></p>
#else if shape ≡ 1
    '<b>Circular</b>
#else if shape ≡ 2
    '<b>Elliptical</b></p>
#else if shape ≡ 3
    '<b>Elliptical tube</b></p>
    'Thickness -'t_b = 20mm', 't_t = 9mm
#else if shape ≡ 4
    '<b>Trapezoidal</b></p>
    'Width at top edge -'b_2 = 90mm
#else if shape ≡ 5
    '<b>T-profile</b></p>
    'Top flange -'b_2 = 70mm', 'h_2 = 20mm
#else if shape ≡ 6
    '<b>Double-T profile</b></p>
    'Bottom flange -'b_1 = 70mm', 'h_1 = 20mm
    'Top flange -'b_2 = 90mm', 'h_2 = 15mm
#end if
'Function of cross-section outline
#if shape ≡ 0
    b(x; z) = w(x)
#else if shape ≡ 1
    b(x; z) = 2*sqrt((h(x)/2)^2 - (z - h(x)/2)^2)
#else if shape ≡ 2
    b(x; z) = 2*w(x)/h(x)*sqrt((h(x)/2)^2 - (z - h(x)/2)^2)
#else if shape ≡ 3
    t(x) = t_b + (t_t - t_b)*x/L
    b_i(x; z) = w(x) - 2*t(x)
    h_i(x; z) = h(x) - 2*t(x)
    b_1(x; z) = 2*w(x)/h(x)*sqrt((h(x)/2)^2 - (z - h(x)/2)^2)
    b_2(x; z) = if((z > t(x))*(z < (h(x) - t(x))); 2*b_i(x; z)/h_i(x; z)*sqrt((h_i(x; z)/2)^2 - (z - h(x)/2)^2); 0m)
    b(x; z) = b_1(x; z) - b_2(x; z)
#else if shape ≡ 4
    b(x; z) = w(z) + z*(b_2 - w(x))/h(x)
#else if shape ≡ 5
    b(x; z) = if(z > h(x) - h_2; b_2; w(x))
#else if shape ≡ 6
    b(x; z) = if(z > h(x) - h_2; b_2; if(z < h_1; b_1; w(x)))
#end if
'<h4>Cross section properties</h4>
δ = 0.01mm
#if shape ≡ 3
    'Area -'A(x) = $Integral{b_1(x; z) @ z = 0mm : h(x)} - $Integral{b_2(x; z) @ z = t(x) : h(x) - t(x)}|cm^2
    'First moment of area -'S(x) = $Integral{b_1(x; z)*z @ z = 0mm : h(x)} - $Integral{b_2(x; z)*z @ z = t(x) : h(x) - t(x)}|cm^3
    'Geometrical center-'z_c(x) = S(x)/A(x)|mm
    'Second moment of area
    I_1(x) = $Integral{b_1(x; z)*(z - z_c(x))^2 @ z = 0mm : h(x)}
    I_2(x) = $Integral{b_2(x; z)*(z - z_c(x))^2 @ z = t(x) : h(x) - t(x)}|cm^4
    I(x) = I_1(x) - I_2(x)
    'First moment of area under z
    S_1(x; z) = $Integral{b_1(x; ζ)*(z_c(x) - ζ) @ ζ = δ : z}
    S_2(x; z) = $Integral{b_2(x; ζ)*(z_c(x) - ζ) @ ζ = t(x) + δ : z}
    'Shear area
    A_Q1(x) = I_1(x)^2/$Integral{S_1(x; z)^2/b_1(x; z) @ z = δ : h(x) - δ}
    A_Q2(x) = I_2(x)^2/$Integral{S_2(x; z)^2/b_2(x; z) @ z = t(x) + δ : h(x) - t(x) - δ}
    A_Q(x) = A_Q1(x) - A_Q2(x)
#else
    'Area -'A(x) = $Area{b(x; z) @ z = 0mm : h(x)}|cm^2
    'First moment of area -'S(x) = $Area{b(x; z)*z @ z = 0mm : h(x)}|cm^3
    'Center of gravity -'z_c(x) = S(x)/A(x)|mm
    'Second moment of area -'I(x) = $Area{b(x; z)*(z - z_c(x))^2 @ z = 0mm : h(x)}|cm^4
    'First moment of area under z -'S_1(x; z) = $Area{b(x; ζ)*(z_c(x) - ζ) @ ζ = δ : z}
    'Shear area -'A_Q(x) = I(x)^2/$Area{S_1(x; z)^2/b(x; z) @ z = δ : h(x) - δ}
#end if
'<!--'PlotWidth = 200','PlotHeight = 200'-->
#def PlotSection$(x$)
    #if shape ≡ 3
        $Plot{0mm|z_c(x$) & b_1(x$; z)/2|z & -b_1(x$; z)/2|z & b_2(x$; z)/2|z & -b_2(x$; z)/2|z & h(x$)/2 - z|z_c(x$) @ z = 0mm : h(x$)}
    #else
        $Plot{0mm|z_c(x$) & b(x$; z)/2|z & -b(x$; z)/2|z & h(x$)/2 - z|z_c(x$) @ z = 0mm : x$}
    #end if
#end def
PlotSection$(0mm)
PlotSection$(L)
'<h3>Mass</h3>
'Distributed mass along height -'m(x) = A(x)*γ_s
'<h3>Solution</h3>
'Number of nodes -'n_J = 11
'Length of one segment -'Δx = L/n_J
'Elevation of node <var>j</var> -'x(j) = Δx*j
'Bending due to horizontal force <var>F</var><sub>j</sub> = 1 at node <var>j</var>
M(j; x) = max(x(j) - x; 0m)
'Flexibility matrix
D(i; j) = $Integral{M(i; x)*M(j; x)/I(x) @ x = 0m : L}*(1/E) + $Integral{1/A_Q(x) @ x = 0m : L}*(1/G)
#hide
D = symmetric(n_J)
#for i = 1 : n_J
    #for j = 1 : n_J
        x_min = min(x(i); x(j))
        D.(i; j) = ($Integral{M(i; x)*M(j; x)/I(x) @ x = 0m : x_min}/E + $Integral{1/A_Q(x) @ x = 0m : x_min}/G)*kN/mm
    #loop
#loop
#show
D'mm/kN
'Mass matrix
d_M.j = $Integral{m(x) @ x = x(j) - Δx/2 : x(j) + Δx/2}
#hide
d_M = vector(n_J)
#for j = 1 : n_J - 1
    d_M.j = $Integral{m(x) @ x = x(j) - Δx/2 : x(j) + Δx/2}/t
#loop
d_M.n_J = $Integral{m(x) @ x = L - Δx/2 : L}/t
M = vec2diag(d_M)
#show
M
'Total mass of structure -'M_tot = sum(d_M)'t
'Eigenvalues
M_sq = sqrt(M)
C = copy(M_sq*D*M_sq; symmetric(n_J); 1; 1)
λ = eigenvals(C*10^-3; -3)
'Natural frequencies
ω = sqrt(1/λ)
'Vibration frequencies
f = ω/(2*π)*Hz
'Vibration periods
T = 1/f
'Mode shapes
V = transp(eigenvecs(C*10^-3; -3))
Φ = inverse(M_sq)*V
X = stack(matrix(1; 3); Φ)
#hide
PlotHeight = 400
PlotWidth = 120
#show
$Plot{spline(i + 1; 1; X)|i*Δx & spline(i + 1; 2; X)|i*Δx & spline(i + 1; 3; X)|i*Δx @ i = 0 : n_J}
'Comparison with ASCE SEI 7/22 <span class="ref">(C26.11-13)</span>
λ_1 = 1.9*exp(-4*(h_t/h_b)) + 6.65/(0.9 + (t_t/t_b)^0.666)
'<p class="ref">(C26.11-12)</p>
'Fundamental natural frequency
n_1 = λ_1/(2*π*L^2)*sqrt(E*I(0m)/m(0m))|Hz
'Fundamental period of vibrations
T_1 = 1/n_1

Spotted an error? Edit these examples.