Saturday, 5 April 2025

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. 4th Ed. The Macmillan Company, New York. 1968. 

Sunday, 23 March 2025

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. 2nd 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.

Thursday, 20 March 2025

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. 2nd Ed. Addison-Wesley Publishing Company. 1990.

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, LabP...