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.

Dynamic Amplitude Response of an Undamped System under Harmonic Force

Topic: Harmonically Excited Vibration

Subject: Mechanical Vibration

Tool: Scilab & QCAD

By: Gani Comia, Mar. 2025

The undamped system under harmonic force is a dynamic system subjected to external force or excitation. This excitation is called the forcing or excitation function. The excitation function is usually time-dependent and may be harmonic or nonharmonic. 

This article presents a dynamic response of a single degree of freedom system under the harmonic motions of the form

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

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

If the force \(F(t)\) acts on an undamped spring-mass system with \(\phi = 0\), 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

The dynamic amplitude of the undamped spring-mass system is a vibration response which refers to the maximum displacement under dynamic conditions subjected to an oscillatory force or motion. This can be expressed as

$$Y = \large \frac{\delta_{st}}{1 - \left(\frac{\omega}{\omega_n}\right)^2} \tag{4}$$

Where:
\(Y\) – dynamic amplitude
\(\delta_{st}\) – static deflection
\(\omega\) – excitation angular frequency
\(\omega_n\) – natural angular frequency

The static deflection and natural angular frequency are defined as follows

$$\delta_{st} = \frac{F_0}{k}\;, \quad \omega_n = \sqrt{\frac{k}{m}} \tag{5}$$

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

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

Case 2. When \(\mathbf{ \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 \(\mathbf{ \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.

• Application

Consider a sample application illustrated in Figure 1, let us analyze its dynamic amplitude response in terms of excitation frequency based on the three types of responses.

Figure 1. Undamped System under Harmonic Force – Sample Problem

Below is the calculation of dynamic amplitude for the illustrated application with the use of Scilab script.

• Scilab Script

// Copyright (C) 2025 - Gani Comia
// Date of creation: 20 Mar 2025
// Dynamic Amplitude from Undamped System under Harmonic Force
clear;clc;
// (1) dynamic amplitude function
function Y = dynAmp(w)
    Y = delta./(1 - (w./wn).^2)
endfunction

// (2) primary parameters
m = 100/9.8                     // kg, mass
k = 2000                        // N/m, spring stiffness or constant
F0 = 25                         // N, static force (constant)
delta = F0./k                   // m, static amplitude or deflection

// (3) secondary parameters
wn = sqrt(k./m)                 // Hz, natural angular frequency
mprintf("Natural frequency, w_n = %3.1f",wn)

// (4) frequency domain
w_l = linspace(0,wn-0.05,100)    // left side of natural frequency
w_r = linspace(wn+0.05,60,100)   // right side of natural frequency

// (5) calculation of dynamic amplitude
Y_l = dynAmp(w_l)
Y_r = dynAmp(w_r)

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

// (7) plot of results
plot(w_l,Y_l,"b-","linewidth",3.5)
plot(w_r,Y_r,"m-","linewidth",3.5)
plot([wn wn],[-1,1],"r--","linewidth",2)
note1 = "Vibration Response of an Undamped System w/ Harmonic Force"
title(note1,"fontsize",3.75)
ylabel("$\Large\text{Dynamic Amplitude (m),}\;\;\mathbf{Y(\omega)}$")
xlabel("$\Large\text{Mechanical Frequency (Hz),}\;\;\mathbf{\omega}$")
xgrid(3,0)

// (8) additional plot information
note2 = "$\LARGE m\ddot{y}(t) + k y(t) = F_0\,cos(\omega\,t)$"
xstring(34,0.425,note2)
note3 = "$\LARGE Y(\omega\;<\;\omega_n)$"
note4 = "$\LARGE Y(\omega\;>\;\omega_n)$"
note5 = "$\LARGE Y(\omega\;=\;\omega_n)$"
legend(note3,note4,note5,with_box=%F)
note6 = "$\Large \text{Resonance}$"
xstring(16,0.4,note6,-90)
xstring(13.5,-0.8,note6,-90)

// (9) axes properties
ax = gca()
ax.x_location = 'origin' 
ax.data_bounds = [0 -1; 60 1];

• Scilab Output (Figure 2)

Figure 2. Vibration Response to Undamped System with Harmonic Excitation

Figure 2 shows the dynamic amplitude in three types of responses to harmonic excitation. For the illustrated problem, a suitable forcing frequency \(\omega\) can be chosen based on the requirements. Resonance frequency shows an infinite displacement.

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

Reference

1. Singiresu S. Rao. Mechanical Vibrations. 2^nd Ed. Addison-Wesley Publishing Company. 1990.

Bearing Selection using Static, Dynamic Load Rating, Rating Life, and Survival Rate

Topic: Rolling Contact Bearings

Subject: Machine Design

Tool: Scilab & QCAD

By: Gani Comia, Mar. 2025

Let us use an illustrated problem shown in Figure 1 to demonstrate the basic process of selecting an appropriate deep groove ball bearing that satisfy the minimum requirements for an application. The design on Figure 1 requires a deep groove ball bearing with a fixed bore diameter, \(d\), subjected to a pure constant and steady radial load, \(F_r\). The assembly requires the outer ring of the bearing to rotate at 1200 \(rpm\).

Figure 1. Deep Groove Ball Bearing Sample Design Requirements.

For the application with a simple pure radial load, a deep groove ball bearing will be used. For this example, a Timken brand will be chosen to select an appropriate bearing number that will satisfy requirements for static, dynamic load rating and bearing rating life.

• Static Load Rating

The basic static load rating, \(C_{0r}\), is based on a maximum contact stress within a non-rotating bearing at the center of the most heavily loaded rolling element and raceway contact. Equation (1) is the requirement for the static load rating.

$$C_{0r} > P_r \tag{1}$$

Where:

\(C_{0r}\) – static load rating, \(kN\)

\(P_r\) – equivalent radial load, \(kN\)

The equivalent radial load is defined in Equation (2)

$$P_r = C_1 V_1 F_r \tag{2}$$

Where:

\(C_1\) – shock and impact factors, \(C_1 = 1.0\) for a constant and steady load

\(V_1\) – race rotation factor, \(V_1 = 1.0\) for inner ring rotation and \(V_1 = 1.2\) for outer ring rotation

Calculating the equivalent radial load of the given design problem will arrived at

$$P_r = (1) (1.2) (1\;\text{kN}) = 1.2 \; \text{kN}$$

From the given \(P_r\), a bearing with static load rating, \(C_{0r}\), is chosen together with the information of bore diameter, \(d = 10 \; mm\). In reference to Timken’s catalog for deep groove ball bearing, there are three possible bearing numbers that satisfy the two requirements and they are Bearing no. 6000, 6200, and 6300 as shown in Table 1.

Table 1. Basic Information of Timken Bearing no. 6000, 6200, and 6300.

• Dynamic Load Rating and Rating Life

Dynamic load rating, designated as \(C\), is defined as the radial load under which a population of bearings will achieve an \(L_{10}\) life of one million revolutions. Bearing life is defined as the length of time, or number of revolutions, until a fatigue spall of \(6 \; {mm}^2\) develops. The rating life, \(L_{10}\), is the life that 90 percent of a group of identical bearings will complete or exceed before a fatigue spall develops. The \(L_{10}\) life is associated with 90 percent reliability for a single bearing under a certain load. It has been calculated using the dynamic equivalent radial load, \(P_r\), and the dynamic load rating, \(C\), based on one million cycles.

$$L_{10} = \left(\frac{C}{P_r}\right)^e \; \left(\frac{1\, \times \,{10}^6}{60 \, n} \right) \tag{3}$$

Where:

\(L_{10}\) – rating life, \(hr\)

\(C\) – dynamic load rating, \(kN\)

\(P_r\) – equivalent radial load, \(kN\)

\(n\) – bearing’s inner or outer race rotation, \(rpm\)

\(e\) -constant, \(e = 3\) for ball bearings and \(e = 10/3\) for roller bearings

Using the rating life defined in Equation (3), the appropriate bearing number can be selected based on the number of hours it can be run until it completes or exceeds the fatigue limit. Below is the calculation for the illustrated problem with the use of Scilab script.

• Scilab Script

// Copyright (C) 2025 - Gani Comia
// Date of creation: 5 Mar 2025
// Deep Groove Ball Bearing Selection
clear;clc;
// (1) static load rating calculation
// (1-1) equivalent radial load function
function Pr=eqRadLoad(Fr)
    C1 = 1.0;                   // constant and steady load factor
    // V1 = 1.0;                // inner ring rotation factor
    V1 = 1.2;                   // outer ring rotation factor
    Pr = C1.*V1.*Fr;            // kN, equivalent radial load
endfunction

Fr = 1;                         // kN, pure radial load
Pr = eqRadLoad(Fr);             // kN, equivalent radial load
mprintf("Equivalent Radial Load, Pr = %3.2f kN\n", Pr);
mprintf("Required Static Load Rating, Cor > %3.2f kN\n", Pr);

// (1-2) design decision
// for the given d = 10 and equivalent radial load for
// deep groove ball bearing, 6000, 6200 and 6300 can be used since
// their static load rating, Cor > Pr.

// (2) dynamic load rating and bearing life calculation
// bearing life, L10 (hrs) as a function of dynamic load ratings
function L10=bLife(Cr, Pr, n)
    e = 3.0;                    // for ball bearings
    L10 = ((Cr./Pr).^e).*(1e6/(60*n));
endfunction

// (2-1) dynamic load ratings for Bearing No. 6000, 6200, and 6300
Cr = [4.60 5.10 8.10];          // kN, dynamic load ratings                 
nCr = length(Cr);               // number of bearing samples
Pr = 0.1:0.005:2.0;             // domain of equivalent radial load
n = 1200;                       // rpm, design rotation

// (2-2) calculation and plot
clf;
f = gcf();
f.figure_size = [700,700];

// (2-3) calculation of bearing life, L10
for i = 1:nCr
    L10(i,:) = bLife(Cr(i),Pr,n);    
end

// (2-4) plot of L10 and Pr
plot("nl",Pr,L10(1,:),"b-","linewidth",3);
plot("nl",Pr,L10(2,:),"r-","linewidth",3);
plot("nl",Pr,L10(3,:),"m-","linewidth",3);
title("Timken Deep Groove Ball Bearing","fontsize",4.0);
xlabel(["$\mathbf{P_r}$",",Equiv Radial Load, kN"],"fontsize",3.5);
ylabel(["$\mathbf{L_{10}}$",",Bearing Life, hr"],"fontsize",3.5);
note1 = "$\Large\text{Bearing no. 6000}$";
note2 = "$\Large\text{Bearing no. 6200}$";
note3 = "$\Large\text{Bearing no. 6300}$";
legend([note1,note2,note3],with_box=%F);
xgrid(3,0);

// (2-5) plot of L10 for Prad = 1.2 kN for 3 bearing nos.
Prad = 1.2;                     // kN, equivalent radial load
for j = 1:nCr
    L10(j) = bLife(Cr(j),Prad,n);
end

// (2-5-1) vertical line or Prad
vX = [Prad Prad]; vY = [1e2 L10(3)];
plot(vX,vY,"k--");
// (2-5-2) horizontal line for L10 for 3 bearing nos.
for k = 1:nCr
    hX = [0 Prad]; hY = [L10(k) L10(k)];
    plot(hX,hY,"k--");
    plot(Prad,L10(k),"marker","o","markerFaceColor","black");
end
// (2-5-3) L10 labels
format(7);
xstring(0.1,L10(1)-375,["$\large \mathbf {L_{10}=}$",string(L10(1))]);
xstring(0.1,L10(2),["$\large \mathbf {L_{10}=}$",string(L10(2))]);
xstring(0.1,L10(3),["$\large \mathbf {L_{10}=}$",string(L10(3))]);
xstring(1.4,1e6,["@ n = ", string(n), "rpm"]);
xstring(0.4,5e6,"https://gani-mech-toolbox.blogspot.com");
mprintf("Bearing no. 6300 has L10 = %6.1f hrs\n", L10(3));

• Scilab Output (Figure 2)

Figure 2. Bearing’s Rating Life for an Equivalent Radial Load.

Though the three bearing numbers satisfy the bore dimensions and static load rating, the suitable bearing can be further trimmed down on the requirements of applicable rating life as shown in Figure 2.

• Bearing Survival

If the machine is assembled with a total number of \(N\) bearings and each has the same reliability, \(R\), then the reliability of the group of bearings or assembly is

$$R_N = R^N \tag{4}$$

Where:

\(R_N\) – reliability of the group of \(N\) bearings

\(R\) – reliability of each bearing, for \(L_{10}\), \(R = 0.90 \;\; \text{or} \;\; 90\,\text{%}\)

\(N\) – total number of bearings in the assembly

Let us assume that there are two similar bearings in the machine assembly, then its reliability is

$$R_N = (0.9)^2 = 0.81 \;\; \text{or} \;\; 81\,\text{%}$$

The presented factors as basis for bearing selection are some of the important considerations. Actual design might need to take additional considerations into account and it is encouraged the reader to refer to the rolling contact bearing manufacturers information for more details.

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

References

1. Timken Engineering Manual, Catalog No. 10424.
2. Timken Deep Groove Ball Bearing, Catalog No. 10857.
3. M.F. Spotts. Design of Machine Elements. 6th Ed. Englewood Cliffs, N.J. 07632. Prentice-Hall, Inc. 1985.
4. J.E. Shigley and C.R. Mischke. Mechanical Engineering Design. 5th Ed. New York, New York 10020. McGraw-Hill Publishing Company. 1989.