Exponential Illusions: Modeling Population Growth with Caution

Topic: Exponential Growth Model and Limitation

Subject: Numerical Methods

Tool: Scilab

By: Gani Comia, July 2025

• Ordinary Differential Equation (ODE) for Exponential Growth

The ODE representing exponential growth is a fundamental model frequently used in real-world application of predictive modeling, particularly in demographics, ecology, and urban planning, to forecast future population sizes based on past trends and other factors.

The population growth is represented by the differential equation shown in Equation (1). The equation signifies that the rate of population change is directly proportional to the current population size.

$$\frac{dP}{dt}=r \, P \tag{1}$$

Where:

\(dP/dt\) – rate of change of the population with respect to time

\(P\) – population size at any given initial time

\(r\) – growth rate, an assumed positive constant that influences population growth

Equation (1) has a general solution as shown in Equation (2).

$$P(t) = P_0 \; e^{r \, t}, \quad P(0) = P_0 \tag{2}$$

Where:

\(P(t)\) – population at any given time, \(t\)

\(P_0\) – population at time, \(t = 0\)

The ODE representing the population growth model has the following limitations if applied in predictive modelling. First, it predicts unlimited and unrealistic growth over the long period of time, which is not sustainable in the long run due to factors resulting to limitation of resources. Second, the growth rate considered constant does not account for changes in birth and death rates or other environmental factors.

Below is the list of models or techniques for estimating or describing the population aside from exponential growth.

1.  Exponential Growth
2. Logistic Growth
3. System Dynamics
4. Naive Bayes and Decision Tree
5. Grey Prediction

• Numerical Solution to ODE using Euler's Method

The Scilab scripts used for calculations implements the Euler method to numerically approximate the solution of a first-order ordinary differential equation modeling exponential growth. Euler method requires an ODE and initial condition in the form of Equation (3).

$$\frac{dP}{dt} = f(t,P) , \quad (t_0 , P_0) \tag{3}$$

The approximate solution using Euler’s method is shown in Equation (4).

$$\frac{\Delta P}{\Delta t} \approx \frac{P_{n+1} - P_n}{h} \approx f(t, P) \tag{4}$$

Where:

\(h\) – step size, \(h = \Delta t\)

\(P_{n+1}\) is the solution to the ODE and can be solved by rearranging the Equation (4) leading to Equation (5).

$$P_{n+1} = P_n + h \; f(t, P) \tag{5}$$

The \(f(t,P)\) requires finding the growth rate, \(r\), for the two conditions as shown in Equation (6).

$$\large r = \frac{\ln \left( {\frac{P_t}{P_0}} \right)}{t} \tag{6}$$

• Philippine Statistics Authority (PSA) Population Estimates

This article will use the exponential growth model to predict and estimates the population growth of the Philippines in comparison with a country’s agency (PSA) estimates. The agency is mandated to conduct national censuses and surveys including those on population.

Table 1 presents the population projections from the PSA up to the year 2055. These projections will be analyzed alongside the model based on exponential growth.

Table 1. Philippine Population Projection until 2055

• PSA Estimates and Exponential Growth Comparison

Figure 1 illustrates the difference between two methods of estimation. The exponential growth model deviates from the PSA estimates after 25 years. Exponential growth models can accurately predict population trends during the initial growth phase, but their reliability decreases as other factors cause deviations from exponential behavior.

Figure 1. Population Growth Model and PSA Estimates

Below is the Scilab script to produce the results shown in Figure 1.

• Scilab Script

// Copyright (C) 2025 - I.S.Comia
// Date of creation: Jul 6, 2025
// Predictive Modeling using Exponential Growth
clear;clc;

// Subroutine
// Population or Exponential Growth Model (1st-Order ODE)
function dPdt=f(t, P)
    global r;
    dPdt = r * P;
endfunction

function r=growthRate(P_0, P_t, t)
    // calculate the estimated growth rate
    // input:
    //      P_0 - quantity at t = 0
    //      P_t - quantity at a given time t
    //      t - time
    // output:
    //      r - estimated growth rate
    r = log(P_t./P_0)./t;
endfunction

// Main function
// Problem case scenario
// Phil. Statistic Authority (PSA) estimated population from year 2020 to 2050
// Primary parameters
actYear = [2020 2025 2030 2035 2040 2045 2050 2055];
// Population in millions
Pop = [109.2027 113.8631 118.8738 123.9636 128.8260 133.0245 136.2989 138.6727];
Yrs = actYear - 2020;                       // years covered
t = Yrs(2)-Yrs(1);
global r;
r = growthRate(Pop(1),Pop(2),t);            // growth rate estimate
T_initial = Yrs(1);                         // initial time
P_initial = Pop(1);                         // initial population (millions)

// Secondary parameters
T_final = Yrs($);                           // years, final time
dt = 1/12;                                  // month, step size,
N = T_final/dt;                             // number of steps

// Euler method solution to IVP ODE
// Initial conditions for time and P-value
T(1) = T_initial;
P(1) = P_initial;
// Solution
for k = 1:N
    dPdt = f(T(k), P(k));
    T(k+1) = T(k) + dt;
    P(k+1) = P(k) + dPdt * dt;
end

// Plotting of ODE solution
clf;
fig = gcf()
fig.figure_size = [700, 700]
plot(T,P,"b-","LineWidth",4)
title("Population Growth Model, Year 2020~2055","FontSize",4)
xlabel("Time (years)","FontSize",3.5)
ylabel("Population (millions)","FontSize",3.5)
xgrid(7,1)

// Population estimate for T_value~5 years or Year 2025
idT = max(find(T <= 5.05));
T_value = T(idT);
P_value = P(idT);
mprintf("For T_value: %4.2f, P_value: %6.3f \n", T_value, P_value)

// Plot of information for Year 2025
plot(T_value,P_value,"rs","MarkerFaceColor","r","MarkerSize",13)
ann_1=xstring(T_value+0.75,P_value-1,["Year:2025,Pop:",string(round(P_value)),"M"])
ann_1.font_size = 2
ann_2 = xstring(20,102.5,"https://gani-mech-toolbox.blogspot.com")
ann_2.font_size = 2
ann_3 = xstring(T_initial-2.5,P_initial-2.5,"Year: 2020")
ann_3.font_size = 2

T_line = [T_value T_value T_initial-5];
P_line = [100 P_value P_value];
plot(T_line,P_line,"r--","LineWidth",1)

// PSA population estimate
plot(Yrs,Pop,"ko","MarkerFaceColor","k","MarkerSize",8)
ann_4 = ["Exponential Growth Model";"Est Pop @ Year 2025";"";"PSA Pop Estimate"]
ann_5 = legend(ann_4,2,with_box=%F)
ann_5.font_size = 2

ax = gca()
ax.data_bounds = [-5 100; 40 150]

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

References

1. James, G., et al., Modern Engineering Mathematics. Addison-Wesley Publishing Company, Inc. 1993.
2. “Philippine Population is Projected to be around 138.67 Million by 2055 under Scenario 2”. Philippine Statistics Authority. January 31, 2024. https://psa.gov.ph/content/philippine-population-projected-be-around-13867-million-2055-under-scenario-2
3. ChatGPT. AI Tools. Accessed on July 6, 2025.
4. Google AI Overview. AI Tools. Accessed on July 9, 2025. 

C++ and Python for Numerical Solution to Spring-Mass-Damper Model

Topic: SDOF Spring-Mass-Damper System
Subject: Numerical Methods & Mechanical Vibration
Tools: C++, CodeBlocks, GitHub Copilot, LabPlot, Python, VS Code, and QCAD

By: Gani Comia, May 2025

• Spring-Mass-Damper Model and Tools for Numerical Solution

Numerical solutions, such as Euler's method or Runge-Kutta method, can be used to determine the free response of a damped single-degree-of-freedom system. These frequently employed numerical techniques are utilized to solve vibration problems for which closed-form analytical solutions are not readily obtainable. Runge-Kutta simulations provide a higher degree of accuracy and quality compared to those generated by the Euler method.

C++ and Python programming language will be used in this article to solve numerically the response from a Spring-Mass-Damper Model using the Runge-Kutta method of 4^th-Order. One may refer to an article in this blog entitled “Simple Harmonic Motion Solution using Runge-Kutta RK4” for the discussion of RK4 numerical method as applied to 2^nd-Order ODE.

As observed, free oscillations in most systems eventually decay and the motion ceases entirely. This phenomenon can be modelled by adding a damper to the spring-mass model of the mechanical systems. The equation of motion of the spring-mass-damper model is represented by

$$m \, y''(t) + c \, y'(t) + k \, y(t) = 0  \tag{1}$$

Where:

\(m\) – mass, \(kg\)

\(c\) – damping coefficient, \(kg/s\)

\(k\) – spring constant, \(N/m\)

\(y’’(t)\) – acceleration as a function of time, \(m/s^2\)

\(y’(t)\) – velocity as a function of time, \(m/s\)

\(y(t)\) – displacement as a function of time, \(m\)

Subject to initial conditions

$$y(0) = y_0 \quad \text{and} \quad y'(0) = v_0  \tag{2}$$

Where:

\(y_0\) – initial position, \(m\), at \(t = 0 \; s\)

\(v_0\) – initial velocity, \(m/s\), at \(t = 0 \, s\)

Equation (1) is also referred to as a damped single-degree-of-freedom (SDOF) system. In analyzing the damped SDOF system, the non-dimensional number, \(\zeta\), called the damping ratio dictates the type of its response. Damping ratio is defined in Equation (3)

$$\zeta = \frac{c}{2 \, \sqrt{k \, m}}  \tag{3}$$

Damping ratio defines the types of motion of the response as follows

1. Underdamped: \(0 < \zeta < 1\)
2. Overdamped: \(\zeta > 1\)
3. Critically Damped:  \(\zeta = 1\)

• Application

Figure 1 shows a sample application where the response in terms of displacement and velocity as a function of time are to be solved.

Figure 1. Schematic of a SDOF Spring-Mass-Damper Vibration Model.

The given sample of damped system can be represented by the Equation (4).

$$3 \, y''(t) + y'(t) + 2 \, y(t) = 0  \tag{4}$$

With initial conditions as follows

$$y(0) = 0 \quad \text{and} \quad y'(0) = 0.25  \tag{5}$$

• Numerical Solution Using C++, Codeblocks, Github Copilot, and LabPlot.

Figure 2 is the plot of the numerical solution or the response on the system in terms of displacement and velocity. The code is a modified solution through the assistance of Github Copilot (AI). Codeblocks is used as the Integrated Development Environment (IDE) for the C++. The calculation results from C++ are exported as .csv file. LabPlot is the application used to import the .csv file to make the plot as shown in Figure 2.

Figure 2. Plot of Numerical Solution to Damped System using C++ and LabPlot.

• C++ Code

#include <iostream>
#include <fstream>
#include <vector>
#include <cmath>
#include <iomanip>

// Define the system of ODEs
// y1' = y2
// y2' = f(t, y1, y2) = -(2/3)*y1-(1/3)*y2

double f1(double t, double y1, double y2) {
    return y2; // y1' = y2
}
double f2(double t, double y1, double y2) {
    // Example: y'' = -(2/3)*y1-(1/3)*y2 (Simple Harmonic Oscillator)
    return -(2.0/3.0)*y1-(1.0/3.0)*y2; // y2' = -(2/3)*y1-(1/3)*y2
}
// Runge-Kutta 4th Order Method
void rungeKutta4(double t0, double y10, double y20, double tEnd, double h, const std::string& filename) {
    double t = t0;
    double y1 = y10;
    double y2 = y20;
    // Open the CSV file for writing
    std::ofstream file(filename);
    if (!file.is_open()) {
        std::cerr << "Error: Could not open file " << filename << " for writing.\n";
        return;
    }
    // Write the header to the CSV file
    file << "t,y1,y2\n";
    // Write the initial conditions
    file << t << "," << y1 << "," << y2 << "\n";
    while (t < tEnd) {
        // Compute RK4 coefficients for y1
        double k1_y1 = h * f1(t, y1, y2);
        double k1_y2 = h * f2(t, y1, y2);
        double k2_y1 = h * f1(t + h / 2.0, y1 + k1_y1 / 2.0, y2 + k1_y2 / 2.0);
        double k2_y2 = h * f2(t + h / 2.0, y1 + k1_y1 / 2.0, y2 + k1_y2 / 2.0);
        double k3_y1 = h * f1(t + h / 2.0, y1 + k2_y1 / 2.0, y2 + k2_y2 / 2.0);
        double k3_y2 = h * f2(t + h / 2.0, y1 + k2_y1 / 2.0, y2 + k2_y2 / 2.0);
        double k4_y1 = h * f1(t + h, y1 + k3_y1, y2 + k3_y2);
        double k4_y2 = h * f2(t + h, y1 + k3_y1, y2 + k3_y2);
        // Update y1 and y2
        y1 += (k1_y1 + 2 * k2_y1 + 2 * k3_y1 + k4_y1) / 6.0;
        y2 += (k1_y2 + 2 * k2_y2 + 2 * k3_y2 + k4_y2) / 6.0;
        // Update time
        t += h;
        // Write the results to the CSV file
        file << t << "," << y1 << "," << y2 << "\n";
    }
    // Close the file
    file.close();
    std::cout << "Results saved to " << filename << "\n";
}

int main() {
    // Initial conditions
    double t0 = 0.0;    // Initial time
    double y10 = 0.0;   // Initial value of y (y1)
    double y20 = 0.25;   // Initial value of y' (y2)
    double tEnd = 20.0; // End time
    double h = 0.1;     // Step size
    // Output file name
    std::string filename = "rk4_2nd_ode_output.csv";
    // Solve the ODE and save results to the CSV file
    rungeKutta4(t0, y10, y20, tEnd, h, filename);
    return 0;
}

• Numerical Solution Using Python, Github Copilot, and VS Code.

Figure 3 is the plot for the response of the same damped system also in terms of displacement and velocity. Github Copilot (Artificial Intelligence) is used to generate the Python script with minor modification. VSCode is used as IDE to run the script and plot the results.

Figure 3. Plot of Numerical Solution to Damped System using Python Script.

• Python Script

import numpy as np
import matplotlib.pyplot as plt

def runge_kutta_4_system(funcs, y0, t0, t_end, h):
    """

    Solve a system of first-order ODEs using the 4th-order Runge-Kutta method (RK4).
    Parameters:
        funcs (list): List of functions representing the system of ODEs.
        y0 (list): Initial values for the system [y1(0), y2(0), ...].
        t0 (float): Initial time.
        t_end (float): End time.
        h (float): Step size.
    Returns:
        tuple: Arrays of t and y values for each variable.
    """
    t_values = np.arange(t0, t_end + h, h)
    y_values = np.zeros((len(t_values), len(y0)))
    y_values[0] = y0
    for i in range(1, len(t_values)):
        t = t_values[i - 1]
        y = y_values[i - 1]
        k1 = [h * f(t, *y) for f in funcs]
        k2 = [h * f(t + h / 2, *(y + np.array(k1) / 2)) for f in funcs]
        k3 = [h * f(t + h / 2, *(y + np.array(k2) / 2)) for f in funcs]
        k4 = [h * f(t + h, *(y + np.array(k3))) for f in funcs]
        y_values[i] = y + (np.array(k1) + 2 * np.array(k2) + 2 * np.array(k3) + np.array(k4)) / 6
    return t_values, y_values

# Example application
if __name__ == "__main__":
    # Define the system of ODEs for 3y'' + y' + 2y = 0
    # Let y1 = y and y2 = -(2.0/3.0)*y1-(1.0/3.0)*y2, then:

    # y1' = y2
    # y2' = -(2.0/3.0)*y1-(1.0/3.0)*y2
    def f1(t, y1, y2):
        return y2
    def f2(t, y1, y2):
        return -(2.0/3.0)*y1-(1.0/3.0)*y2

    # Initial conditions
    y0 = [0, 0.25]  # y(0) = 0, y'(0) = 0.25
    t0 = 0       # Initial time
    t_end = 20   # End time
    h = 0.1      # Step size
    # Solve the system
    t_values, y_values = runge_kutta_4_system([f1, f2], y0, t0, t_end, h)
    # Extract y1 (the solution to the original ODE)
    y1_values = y_values[:, 0]
    y2_values = y_values[:, 1]

    # Plot the solution
    plt.figure(figsize=(8, 6))
    plt.plot(t_values, y1_values, linewidth=3, label="y(t)", color="blue")
    plt.plot(t_values, y2_values, linestyle='--', linewidth=2, label="y'(t)", color="red")
    plt.title("Using Python to Solve 2nd-Order ODE with RK4")
    plt.xlabel("Time, t")
    plt.ylabel("Displacement, y(t) & Velocity, y'(t)")
    plt.grid(True)
    plt.ylim(-0.3, 0.3)
    plt.legend()
    plt.show()

The damping ratio for the given application can be solved as shown

$$\zeta = \frac{c}{2 \, \sqrt{k \, m}} = \frac{1}{2 \sqrt{(2)(3)}} = 0.204$$

With the damping ratio, \(\zeta = 0.204\), the motion of the response of the system is considered underdamped.

C++ and Python programming languages using Codeblocks and VS Code as IDEs, and LabPlot for plotting application are some powerful engineering tools to analyze a system of mechanical vibration. AI through the use of Github Copilot provided the applicable codes for the solution.

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

References

1. D.I.Inman. Engineering Vibration. 2^nd Ed. Prentice Hall International, Inc., New Jersey 07458. 2001.
2. GitHub Copilot. AI Tools. Accessed on May 8, 2025 within the VS Code environment. 

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