Wind Load Effect Simulation on an Outdoor Gantry Crane based on TCWS

Topic: Equilibrium of Force Systems and Wind Loads

Subject: Fluid and Engineering Mechanics

Tool: QCAD & Scilab

By: Gani Comia, May 2025

• Gantry Crane Structure and Equilibrium of Force Systems

Gantry cranes are travelling cranes designed for lifting and moving heavy loads. This material handling equipment is generally used outdoors where it is not convenient to erect an overhead runway. Its structure consists of a girder or bridge, freestanding legs, electric motor that drives the steel wheel, and a rail system on the ground.  The bridge or girder is carried at the ends by the legs supported by trucks with wheels so that the crane can travel. The bridge carries a hoisting assembly unit. The crane is driven by the motor through a gear reduction shaft. [1]. The motor is referred to as a brake motor because it prevents movement caused by inertia when it is deenergized.

The objective of this article is to present an engineering simulation method of analyzing the effect of a wind loads due to typhoon to an outdoor gantry crane. The assumption is that there is a particular value of wind velocity that will make the crane to tip over at the wind direction.

Figure 1 shows a simple illustration of an outdoor gantry crane. On its side view is the free-body diagram (FBD) showing the corresponding concurrent forces acting on the structure. The following conventions will be used in the analysis: - \(W\) is the weight of the crane; \(F_{g h}\) is the horizontal induced wind load on the girder; \(F_{s h}\) and \(F_{s v}\) are the horizontal and vertical wind load acting on the support legs respectively; and \(R_{a c}\) and \(R_{b d}\) are the reaction forces at the four wheels, \(a\), \(b\), \(c\), and \(d\).

Figure 1. Outdoor Gantry Crane and Free-Body Diagram for the Forces System.

The solution is to determine as to whether the crane will be tipping over at the center of rotation, \(R_{a c}\), due to moment on the direction of the wind induced forces. Figure 2 illustrates two cases of wind velocity, \(V = 0\) and \(V > 0\). For the given domain of \(V\), there is a reaction force, \(R_{b d}\), where in there is a resultant moment from the coplanar force system that determines the rotation of crane from the center.

Figure 2. FDB of Gantry Crane due to Effect of Wind Velocity.

Case #1 is a condition where the wind velocity is \(V = 0\). With this there will be no induced forces and the force system present on FBD are crane weight, \(W\), and the two reaction forces, \(R_{a c}\) and \(R_{b d}\), on the wheel supporting the structure. In this condition, the equilibrium for concurrent force systems is obtained by determining the equations that produce a zero resultant. The resultant will be zero and equilibrium will exist when the following equations are satisfied:

$$\sum F_x = 0 \tag{1}$$

$$\sum F_y = 0 \tag{2}$$

Where:

\(F_x\) – forces parallel to x-axis

\(F_y\) – forces parallel to y-axis

With the two conditions of equilibrium, reaction forces can be determined using Equation (3).

$$R_{a c} = R_{b d} = \frac{W}{2} \tag{3}$$

Case #2 is the condition wherein the wind velocity acting on the crane is inducing resultant forces that make the structure to move by rotation. In this condition, the equilibrium in terms of moment summations about the center on its line of action can be used to calculate the resultant force. Hence the two other equations of equilibrium are

$$\sum M_A = 0 \tag{4}$$

$$\sum M_B = 0 \tag{5}$$

Where:

\(M_A\) – moment at center A

\(M_B\) – moment at center B

For the case #2 on Figure 2, reaction \(R_{a c}\) is used as the center of moment. For the increasing wind velocity, \(V_{wind}\), there exist a corresponding increase in wind load represented by \(F_{g h}\), \(F_{s h}\), and \(F_{s v}\) as they are functions of wind velocity in Equation 6.

$$F_{g h} = f_{g h} (V_{wind}), \quad F_{s h} = f_{s h} (V_{wind}), \quad F_{s v} = f_{s v} (V_{wind}) \tag{6}$$

$$F_{gh} \; \uparrow , F_{s h} \; \uparrow , F_{s v} \; \uparrow \quad \text{with respect to} \quad V_{wind} \tag{7}$$

• Wind Load Analysis [3]

Wind is a mass of air that moves in a mostly horizontal direction from an area of high pressure to an area with low pressure. High winds can be very destructive because they generate pressure against the surface of a structure. The intensity of this pressure is the wind load. The wind effect is dependent upon the size and shape of the structure. Calculating wind load is necessary for the design and construction of safer, more wind-resistant building and placement of objects such as antennas on top of the buildings.

The generic formula for wind load is

$$F = A \; P \; C_d \tag{8}$$

Where:

\(F\) – wind load or force, \(N\)

\(A\) – projected area of the object, \(m^2\)

\(P\) – wind pressure, \(N/m^2\)

\(C_d\) – drag coefficient

The generic formula is used for estimating the wind load on a specific object and it does not meet the building code requirements for planning a new construction.

The formula for wind pressure, \(P\), in SI units is

$$P = \text{0.613} \; V^2 \tag{9}$$

Where:

\(V\) – wind speed, \(m/s\)

The wind speed categorizes by the Philippine Atmospheric, Geophysical and Astronomical Services Administration (PAGASA) will be used as basis in calculation. PAGASA is the Philippine government agency responsible for providing weather forecasts, typhoon warnings, and other meteorological and astronomical services. Currently the Tropical Cyclone Wind Signal (TCWS) has been in use based on the adoption of best practices from other tropical cyclone warning centers.

The following are the corresponding wind threat based on the wind speed value.

TCWS #1: Wind Threat = 10.8 ~ 17.1 \(m/s\)

TCWS #2: Wind Threat = 17.2 ~ 24.4 \(m/s\)

TCWS #3: Wind Threat = 24.5 ~ 32.6 \(m/s\)

TCWS #4: Wind Threat = 32.7 ~ 51.2 \(m/s\)

TCWS #5: Wind Threat = 51.3 \(m/s\) or higher

\(C_d\) is the drag coefficient for the object subjected under wind pressure. Drag is the force that air exerts on the structure. Drag coefficient for the structure’s shape can be estimated roughly using the following standard:

\(C_d\) = 0.8 for short cylinder

\(C_d\) = 1.2 for long cylinder

\(C_d\) = 1.4 for shorter flat plate

\(C_d\) = 2.0 for long flat plate

• Wind Load Simulation for Gantry Crane

Figure 3 shows the effect of wind load on the gantry crane structure in terms of the equilibrium of force systems.

Figure 3. Wind Load Simulation using the Equilibrium of Force System on Gantry Crane.

The graph’s domain is based on the TCWS converted in terms of \(km/hr\). Corresponding wind speeds or velocities for the five categories of TCWS are plotted. Wheel reaction load for \(R_{b d}\) are shown for each TCWS. The red dash line which depicts \(R = 0\) is the critical wheel reaction load for \(R_{b d}\) on the FBD in Figure 2. With this \(R = 0\), the critical wind velocity is calculated to have a value of \(\text{127}\) \(km/hr\) under the TCWS #4.  A horizontal wind speed which of \(V_{wind} >\) \(\text{127}\) \(km/hr\) acting perpendicularly to structure will result to tripping over the crane to the wind direction.

With this, a conclusion can be made from the simulation stated as follows: - a wind speed of over \(\text{127}\) \(km/hr\) which is within the TCWS #4 will induce an equivalent critical wind load to affect the crane stability. The use of a device in gantry crane known as “Storm Lock” is highly recommended to increase its capacity for critical wind load.

Below is the Scilab script to generate the simulation plot on Figure 3. Parameters can be edited for the other cases.

• Scilab Script

// Copyright (C) 2025 - Gani Comia
// Date of creation: 24 Apr 2025
// Wind Load Effect Analysis on Outdoor Gantry Crane
clear;clc;

// gantry crane primary parameters for girders
lG = 1.35+22+0.35;                      // m, crane girder length
hG = 1;                                 // m, crane girder height
W = 14.5;                               // Tonne, crane weight
shapeG = 2;                             // crane girder as long flat plate
// gantry crane primary parameters for support legs
lS = 5.7;                               // m, support leg length
dS = 0.25;                              // m, support leg diameter
shapeS = 1.2;                           // leg support as long cylinder tube

// wind velocity domain
Vmps = 0:52;                            // m/s, wind speed, tropical cyclone
Vkph = Vmps*3.6;                        // km/hr, wind speed, tropical cyclone

// modules or functions
// wind load analysis using generic formula
// structure surface area
function A=area(l, h)
    // calculate girder area perpendicular to wind velocity
    // input argument:
        // l - surface length (m)
        // h - surface height (m)
    // ouput argument:
        // A - surface area (m^2)
    A = l.*h;
endfunction

// wind pressure
function P=windPressure(V)
    // calculate the induced wind pressure
    // input argument:
        // V - wind velocity (m/s)
    // ouput argument:
        // P - wind pressure (N/m^2 or Pa)
    P = 0.613*V.^2;
endfunction

// drag coefficient
function Cd=dragCoefficient(shape)
    // calculate the theortical drag coefficient
    // input argument: structure shape
        // shape = 2 for long flat plate
        // shape = 1.4 for shorter flat plate
        // shape = 1.2 for long cylinder tube
        // shape = 0.8 for short cylinder tube
    // output argument:
        // Cd - drag coefficient (no unit)
    Cd = shape;
endfunction

// induced wind load
function Fw=windLoad(A, P, Cd)
    // calculate the induced wind load
    // input argument:
        // A - surface area (m^2)
        // P - wind pressure (N/m^2)
        // Cd - drag coefficient (no unit)
    // output argument:
        // Fw - wind load (N)
    Fw = A.*P.*Cd;               
endfunction

// main function
// wind load calculation at the girder
A_g = area(lG,hG);
P_g = windPressure(Vmps);
Cd_g = dragCoefficient(shapeG);
Fgh = windLoad(A_g,P_g,Cd_g);
Fgh = Fgh/9806.65;                      // Tonne, wind load    

// wind load calculation at the support legs
A_s = area(lS,dS);                      // m^2, surface area
// angle of support leg
phi = atand(((3.5/2)-0.168)/5.495);     // degrees, angle
Vmps_s = Vmps./cosd(phi);               // m/s, wind speed perpendicular to legs
P_s = windPressure(Vmps_s);
Cd_s = dragCoefficient(shapeS); 
Fsr = windLoad(A_s,P_s,Cd_s);

// force component at legs
Fsh = Fsr*cosd(phi); Fsh = 3*Fsh/9806.65;   // N to tonnes
Fsv = Fsr*sind(phi); Fsv = 3*Fsv/9806.65;   // N to tonnes

// critial load
Rbd = (W*(3.5/2)+Fsv*((3.5/2)+0.937)-Fgh*(6.25+0.375)-Fsh*(2.9+0.375))/3.5;

// Plotting calculation results
clf;
fig = gcf()
fig.figure_size = [700,700]
plot(Vkph,Rbd,"b-","linewidth",4) 
note = "$\LARGE R_{wheel} = f(V_{wind})$"
legend(note,with_box=%F)
plot([0,200],[0,0],"r--","linewidth",1.8)
title("Wind Load Effect on Outdoor Gantry Crane","fontsize",4)
xlabel("Wind Velocity, Tropical Cyclone, V, (kph)","fontsize",3.5)
ylabel("Wheel Reaction Load, R, (Tonne)","fontsize",3.5)
xgrid(color("green"),1,7)
xstring(10,0,"Critical Wheel Reaction Load = 0")
xstring(120,5.0,"https://gani-mech-toolbox.blogspot.com")

// Wind velocity limit of tropical cyclone
Vwind = [10.8,17.2,24.5,32.7,51.3];         // m/s, wind velocity
Vwind = Vwind*3.6;                          // kph, wind velocity 

// tropical cyclone signal number
nVwind = length(Vwind)

for i = 1:nVwind
    R(i) = interp1(Vkph,Rbd,Vwind(i))
    vLx = [Vwind(i) Vwind(i)]; vLy = [-10 R(i)]
    hLx = [0 Vwind(i)]; hLy = [R(i) R(i)]
    plot(vLx,vLy,"k--")
    plot(hLx,hLy,"k--")
    plot(Vwind(i),R(i),"marker","d","markerFaceColor","red")
    xstring(Vwind(i)+8,-9.8,["Signal No. ",string(i)],-90)
end

// critical wind velocity
format(5)
Rcrit = 0
Vcrit = interp1(Rbd,Vkph,Rcrit)
plot(Vcrit,Rcrit,"marker","S","markerFaceColor","red")
xstring(Vcrit,Rcrit,["Critical Velocity (kph) =",string(Vcrit)])

Feel free to comment for inquiry, clarification, correction or suggestion for improvement. Drop your email to make a request to the author.

References

1. E.A. Avallone and T. Baumeister III. Marks Standard Handbook for Mechanical Engineers. 9^th Ed. McGraw-Hill International Edition. 1987.
2. F.L. Singer. Engineering Mechanics. 2^nd Ed. Harper International Edition, New York, Evanston & London. 1970. 
3. Joseph Quinones. “How to Calculate Wind Load”, last updated July 24, 2023. https://www.wikihow.com/Calculate-Wind-Load#Calculating-Wind-Load-Using-the-Generic-Formula 
4. “Tropical Cyclone Wind Signal.” Philippine Atmospheric, Geophysical and Astronomical Services Administration, March 23, 2022. https://www.pagasa.dost.gov.ph/learning-tools/tropical-cyclone-wind-signal 

Design of a Circular Flat Plate Cover for Vacuum Chamber using Acrylic Material

Topic: Flat Plate & Maximum Principal Stress Theory
Subject: Machine Design
Tool: QCAD & Scilab

By: Gani Comia, Apr. 2025

• Circular Flat Plate Application

This topic will present the solution for the design of an acrylic material used as a cover in the vacuum chamber. Figure 1 illustrates an application wherein the thickness of circular plate cover is calculated for the considerable factor-of-safety in order to overcome the induced stress from a vacuum pressure.

Figure 1. Vacuum Chamber with Circular Plate Cover Application.

For a circular plate under uniform load and simply supported at the edge, the maximum stress at the center is

$$\sigma = \frac{3\,(3 + \nu)\,p \,r^2 }{8\,t^2} \tag{1}$$

Where:

\(\sigma\) – maximum bending stress at the center, \(MPa\)

\(\nu\) – Poisson’s ratio of the material

\(p\) - uniform load or pressure, \(MPa\)

\(r\) - radius of circular plate, \(mm\)

\(t\) – thickness of plate, \(mm\)

• Maximum Principal Stress Theory

This failure theory is accredited to W.J.M. Rankine. It gives good predictions and often used for brittle materials. It presumes that when the maximum principal stress exceeds a certain limiting value, for this case the ultimate strength, \(\sigma_u\), for brittle materials failure occurs. This design stress for uniaxial load with shear is given by the Equation (2)

$$\sigma_d = \frac{\sigma}{2} + \sqrt{\left( \frac{\sigma}{2} \right)^2 + \tau^2} \tag{2}$$

Where:

\(\sigma_d\) – design stress, \(MPa\)

\(\sigma\) – normal stress, \(MPa\)

\(\tau\) – shear stress, \(MPa\)

The normal stress from Equation (2) is the maximum bending stress, \(\sigma\), as in this case of circular plate cover. Shear stress, \(\tau\), is calculated from the force due to the difference between atmospheric pressure and vacuum acting on the area around the chamber opening circumference with the given thickness of the cover plate. Shear stress is defined as

$$\tau = \frac{F}{A}, \quad F = \pi\,r^2\,p, \quad A = 2\,\pi\,r\,t \tag{3}$$

Where:

\(F\) – induced force due to pressure difference, \(N\)

\(A\) – area under shear at the cover around the circumference of vacuum chamber, \({mm}^2\)

For a static loading of a brittle material the limiting stress is taken as the ultimate strength and the design stress is expressed as

$$\sigma_d = \frac{\sigma_u}{N} \tag{4}$$

Where:

\(\sigma_u\) – ultimate strength, \(MPa\)

\(N\) – factor-of-safety

Our objective for the design is to determine the thickness of the circular plate cover which will give a considerable factor-of-safety while addressing other factors such as stress concentration due to the effect of attaching vacuum fittings or other devices and the availability of material based on standard or off-the-shelf thickness.  

• Scilab Script

// Copyright (C) 2025 - Gani Comia
// Date of creation: 2 Apr 2025
// Circular Plate, Uniform Load, Simply Supported Edge
clear;clc;
// module or functions
// (1) maximum stress at the center due to pressure
function sigma_fp=flatPlate(p, r, t, nu)
    // calculate the stress on the flat plate.
    // input argument:
        // p - MPa or N/mm^2, uniform load
        // r - mm, radius of circular plate
        // t - mm, plate thickness
        // nu - Poisson's ratio
    // output argument:
        // sigma_fp - MPa, maximum principal stress at the center
    sigma_fp = (3*(3+nu).*p.*r.^2)./(8*t.^2);
endfunction

// (2) maximu shear stress at the chamber opening
function tau_fp=shearPlate(p, r, t)
    // calculate the shear at the plate around the chamber opening.
    // input argument:
        // p - MPa or N/mm^2, uniform load
        // r - mm, radius of circular plate
        // t - mm, plate thickness
    // output argument:
        // tau_fp - MPa, maximum shear stress    
    c = 2*%pi*r;            // mm, circumference
    a_shear = c.*t;         // mm^2, area under shear
    a_normal = %pi*r.^2;    // mm^2, area perpendicular to pressure
    F = p.*a_normal;        // N, force due to pressue
    tau_fp = F./a_shear;
endfunction

// (3) maximum principal stress theory
function sigma_d=max_Pst(sigma, tau)
    // calculate the design stress based on max principal stress theory
    // input argument:
        // sigma - MPa, principal stress (uniaxial)
        // tau - MPa, shear stress
    // output argument:
        // sigma_d - MPa, design stress based on MPST
    sigma_d = (sigma/2)+sqrt((sigma/2).^2 + tau.^2);
endfunction

// (4) factor of safety
function N=FoS(sigma_d, sigma_u)
    // calculate the factor of safety for brittle material
    // input argument:
        // sigma_d - MPa, design stress
        // sigma_u - MPa, ultimate strength
    // output argument:
        // N - no unit, factor-of-safety 
    N = sigma_u./sigma_d;   
endfunction

// main function
// (5) primary parameters
p_atm = 760;                // Torr, atmospheric pressure
p_vac = 0.001;              // Torr, vacuum pressure
r = 500/2;                  // mm, radius of flat circular cover
nu = 0.35;                  // Poisson's ratio of acrylic or
                            // polymethyl methacrylate of PMMA
sigma_u = 75;               // MPa, ultimate tensile strength of PMMA

// (6) secondaty parameters
p_eqv = p_atm-p_vac;            // Torr, equivalent pressue
p_eqv = p_eqv*0.000133322;      // MPa, equivalent pressure

// (7) thickness domain under examination
t = linspace(1,75,200)

// (8) calculation
sigma_fp = flatPlate(p_eqv,r,t,nu)
tau_fp = shearPlate(p_eqv,r,t)
sigma_d = max_Pst(sigma_fp,tau_fp)
N = FoS(sigma_d,sigma_u)

// (9) visualization of results
clf;
fig = gcf()
fig.figure_size = [700,700];

plot(t,N,"b-","linewidth",3.5)
title("Thickness of Vacuum Chamber Cover","fontsize",4.5)
xlabel("Design Thickness (t), mm","fontsize",4)
ylabel("Factor-of-Safety (N)","fontsize",4)
xgrid(3,0)
note="$\LARGE N = f(t),\quad\text{with Max Principal Stress Theory}$"
legend(note,with_box = %F)
xstring(50,53,"https://gani-mech-toolbox.blogspot.com")

// (10) plot of standard thickness and factor-of-safety
std_thk = [25 50 75]
Nstdthk = length(std_thk)
format(5)

for i = 1:Nstdthk
    N(i) = interp1(t,N,std_thk(i))
    vLx = [std_thk(i) std_thk(i)]; vLy = [0 N(i)];
    hLx = [0 std_thk(i)]; hLy = [N(i) N(i)];
    plot(vLx,vLy,"r--")
    plot(hLx,hLy,"r--")
    plot(std_thk(i),N(i),"marker","o","markerFaceColor","red")
    xstring(5,N(i),["$\large N = $",string(N(i))])
    xstring(std_thk(i)+3.5,0,["$\large t = $",string(std_thk(i))],-90)
end

ax = gca()
ax.data_bounds = [0 0; 80 60];

• Scilab Output (Figure 2)

Figure 2. Factor-of-Safety as a function of thickness for the design of vacuum chamber cover.

A vacuum chamber cover with the thickness of \(\text{50 mm}\) would be reasonable enough for the consideration such that its weight will provide an initial compression set for gasket in between the chamber and the cover, but a toggle clamp or any clamping mechanism is still necessary for a leak free sealing, and the provision of fittings for inlet/outlet vacuum port and instruments at the chamber cover leads to an increase in stress concentration that in effect reduce the factor-of-safety.

Feel free to comment for inquiry, clarification, correction or suggestion for improvement. Drop your email to make a request to the author.

References

1. M.F. Spotts. Design of Machine Elements. 6th Ed. Englewood Cliffs, N.J. 07632. Prentice-Hall, Inc. 1985.
2. Virgil Moring Faires. Design of Machine Elements. 4^th Ed. The Macmillan Company, New York. 1968. 

Solving an Undamped System Under Harmonic Force Using Laplace Transform

Topic: Harmonically Excited Vibration
Subject: Mechanical Vibration
Tool: Euler Math Toolbox, QCAD, & Scilab

By: Gani Comia, Mar. 2025

• Undamped System Under Harmonic Force - Equation of Motion

The undamped system under harmonic force is a dynamic system subjected to external force or excitation. Consider the harmonic excitations of the form

$$F(t) = F_0 \cos(\omega t) \tag{1}$$

Where:
\(F(t)\) – excitation of forcing function
\(F_0\) – constant static force
\(\omega\) – excitation angular frequency
\(t\) - time

If the force \(F(t)\) acts on an undamped spring-mass system, the equation of motion is

$$m\,y'' + k\,y = F(t) \tag{2}$$

or

$$m\,y'' + k\,y = F_0 \cos(\omega t) \tag{3}$$

Where:
\(m\) – mass
\(k\) – spring constant
\(y’’\) – acceleration of mass
\(y\) – displacement from rest position

Under a harmonic excitation, the response of the system will also be harmonic. If the frequency of excitation \(\omega\) coincides with the natural frequency \(\omega_n\) of the system, the response of the system will be very large. The condition, \(\omega = \omega_n\), is what we called resonance. Resonance is to be avoided to prevent system’s failure as in this case mechanical damage.

$$\omega = \omega_n \quad \left(\omega_n = \sqrt{k / m}\right) \tag{4} $$

Where:
\(\omega_n\) – natural angular frequency

• Laplace Transform of Second-Order Linear Differential Equation

Laplace transform is one of the methods to solve for the response for the undamped system under harmonic force. Consider a general second-order linear differential equation

$$a\,\frac{d^2 y}{d t^2} + b\, \frac{dy}{dt} + c\,y = u(t) \quad \left(t \geq 0 \right) \tag{5}$$

Subject to initial conditions

$$t = 0, \quad y(0) = y_0, \quad \frac{dy}{dt}(0) = \nu_0 \tag{6}$$

Such differential equations may model the dynamics of some system for which the variable \(y(t)\) determines the response of the system to the forcing or excitation term \(u(t)\). The terms system input and system output are also frequently used for \(u(t)\) and \(y(t)\) respectively.

• Laplace and Inverse Laplace Transform of an Undamped System Under Harmonic Force

Taking the Laplace transform of each term of Equation (3)

$$m\, \mathcal{L} \left\{ \frac{d^2 y}{d t^2} \right\} + k\, \mathcal{L} \{y\} = F_0 \, \mathcal{L} \{ \cos(\omega t)\} \tag{7}$$

With initial conditions

$$t = 0, \quad y(0) = \frac{F_0}{k - m\,\omega^2}, \quad \frac{dy}{dt}(0) = 0 \tag{8}$$

Leads to

$$m \, \left( s^2 \, Y(s) - \frac{s\,F_0}{k - m \, \omega^2} \right) + k\,Y(s) = \frac{s \, F_0}{\omega^2 + s^2} \tag{9}$$

By rearranging and solving for \(Y(s)\) gives

$$Y(s) = - \frac{s\,F_0}{m\,\omega^4 + \omega^2 \, (m\,s^2 - k) - k\,s^2} \tag{10}$$

The inverse Laplace transform solves the system output \(y(t)\) for \(\omega > 0\)

$$\mathcal{L^{-1}} \{Y(s)\} = \mathcal{L^{-1}} \left\{ - \frac{s\,F_0}{m\,\omega^4 + \omega^2 \, (m\,s^2 - k) - k\,s^2} \right\} \quad (\omega > 0) \tag{11}$$

The steady-state response \(y(t)\) is

$$y(t) = \frac{F_0 \, \cos(\omega \, t)}{m \, \omega^2 - k} \tag{12}$$

The response of the system can be identified to be of three types.

Case 1. When \(0 < {\omega}/{\omega_n} < 1\), the harmonic response of the system \(y(t)\) is said to be lagging the force.

Case 2. When \({\omega}/{\omega_n} > 1\), the response of the system to a harmonic force of very high frequency is close to zero. The response of the system is said to be leading the force.

Case 3. When \({\omega}/{\omega_n} = 1\), the response is called resonance and the amplitude of \(y(t)\) becomes infinite. This phenomenon can result on large displacements and stresses.

Consider a sample application illustrated in Figure 1, let us analyze the response \(y(t)\) in terms of excitation \(F(t)\) based on the Case 1 and 2 responses. Case 3 will not be possible to plot. For this let us assume an excitation frequency, \(\omega\), of less than and more than the natural frequency, \(\omega_n\).

Figure 1. Application of an Undamped System under Harmonic Force

Below is the Scilab script for plotting both the system input and output for the illustrated problem.

• Scilab Script

// Copyright (C) 2025 - Gani Comia
// Date of creation: 22 Mar 2025
// Undamped System under Harmonic Force using Laplace Transform
clear;clc;
// (1) solution y(t) from m y''(t) + k y(t) = F0 cos(wt)
function y=position(t)
    y = -(F0.*cos(omega.*t))./((m.*omega.^2)-k)
endfunction

// (2) excitation of forcing function
function F=force(t)
    F = F0.*cos(omega.*t)
endfunction

// (3) primary parameters from m y''(t) + k y(t) = F0 cos(wt)
m = 100/9.8;                        // kg, mass (where W = 100 N)
k = 2000;                           // N/m, spring constant
F0 = 25;                            // N, constant static force

// (4) secondary parameters
omega_n = sqrt(k./m);               // Hz, natural angular frequency
mprintf("omega_n (Hz) = %3.1f\n", omega_n)

clf;
// (5) figure properties
fig = gcf()
fig.figure_size = [700,800];

// (6) for excitation or forcing function
omega = omega_n
omega_t = linspace(0,2*%pi,100);
t = omega_t./omega;

// (6-1) plot of excitation force
subplot(3,1,1)
F = force(t)
plot(t,F,"m-","linewidth",4)
title("$\Large\text{Excitation Function}\quad F(t)=F_0\cos(\omega t)$")
ylabel("$\Large \mathbf{F(t)}$")
xlabel("$\Large \text{time,}\; t$")
legend("$\LARGE\omega = \omega_n$",4,with_box=%F)
xgrid(3,0)
// (6-2) axes properties
ax = gca()
ax.x_location = 'origin' 
ax.data_bounds = [0 -30; 0.5 30]

// (7) for Case 1. 0 < omega < omega_n
omega = 0.95*omega_n
omega_t = linspace(0,2*%pi,100);
t = omega_t./omega;

// (7-1) plot of a sampe case 1
subplot(3,1,2)
y = position(t)
plot(t,y,"b-","linewidth",4)
title("$\Large\text{Displacement at}\quad \mathbf{\omega < \omega_n}$")
ylabel("$\Large \mathbf{y(t)}$")
xlabel("$\Large \text{time,}\; t$")
legend("$\LARGE \omega = \text{0.95} \; \omega_n$",4,with_box=%F)
xgrid(3,0)
ax = gca()
ax.x_location = 'origin'
ax.data_bounds = [0 -0.2; 0.5 0.2]

// (8) for Case 2. omega > omega_n
omega = 1.05*omega_n
omega_t = linspace(0,2*%pi,100);
t = omega_t./omega;

// (8-1) plot of a sampe case 2
subplot(3,1,3)
y = position(t)
plot(t,y,"r-","linewidth",4)
title("$\Large\text{Displacement at}\quad \mathbf{\omega > \omega_n}$")
ylabel("$\Large \mathbf{y(t)}$")
xlabel("$\Large \text{time,}\; t$")
legend("$\LARGE \omega = \text{1.05} \; \omega_n$",1,with_box=%F)
xstring(0.31,-0.2,"https://gani-mech-toolbox.blogspot.com")
xgrid(3,0)
ax = gca()
ax.x_location = 'origin'
ax.data_bounds = [0 -0.2; 0.5 0.2] 

• Scilab Output (Figure 2)

Figure 2. Forcing Function and Responses to Undamped System with Harmonic Excitation

The harmonic steady-state response \(y(t)\) for case 1 is said to be lagging with the forcing function \(F(t)\). As for the case 2 the response \(y(t)\) is said to be leading with external harmonic force \(F(t)\).

Feel free to comment for inquiry, clarification, or suggestion for improvement. Drop your email to request the soft copy of the file.

References

1. Singiresu S. Rao. Mechanical Vibrations. 2^nd Ed. Addison-Wesley Publishing Company. 1990.
2. G. James, D. Burley, P. Dyke, J. Searl, N. Steele, J. Wright. Advanced Modern Engineering Mathematics. Addison-Wesley Publishing Company, Inc. 1993.