---
tags: Mathematical modeling and optimization
title: 2.1 Analyzing mechanical systems
---
2.1 Analyzing mechanical systems
===
Analyzing mechanical systems is a central task in engineering design and problem-solving. Engineers employ various methods and concepts to understand the behavior of these systems, ensuring they meet desired performance criteria and safety standards. In this discussion, we'll explore a comprehensive approach to analyzing structural systems, starting with the Free Body Diagrams technique and progressing to systems with multiple objects.
### 2.1.1 Free Body Diagrams (FBD)
A Free Body Diagram is a graphical representation used to isolate a body from its environment, allowing us to analyze the forces acting on that body. It's a simplified sketch that includes all the external forces acting on the body, such as applied forces, reactions, and gravitational forces.
Free Body Diagrams are essential because they provide a clear visual representation of the forces at play in a given system. This simplification allows engineers to apply the principles of equilibrium and Newton's laws of motion to analyze and solve complex engineering problems.
__Components of a Free Body Diagram:__
- __Body Representation__:
Begin by drawing a simple, labeled representation of the body in question. This can be a point particle, a rigid body, or a specific component of a larger system.
- __Forces__:
Identify and draw all the external forces acting on the body. These include applied forces, reactions, and gravitational forces. Clearly label each force vector with appropriate symbols.
- __Coordinate System__:
Establish a coordinate system to define the direction of forces accurately. This makes it easier to analyze and resolve forces into their respective components.
A simple example demonstrates a box resting on an inclined plane.
<img src="https://live.staticflickr.com/65535/53558688349_76441b0151_z.jpg">
where $W$ is the weight of box, $F_N$ is the normal force to the incline, and $F_f$ is the friction of the surface, acting in the opposite direction of the movement.
By analyzing the Free Body Diagram, we can apply Newton's laws of motion and equilibrium equations to determine:
- Whether the box is at rest or in motion.
- The minimum angle of inclination before the box starts sliding (if friction is present).
- The forces acting on the box and their magnitudes and directions.
__Forces acting on a single body__
In force analysis based on FBDs, Newton's second law is applied to relate the forces depicted in the diagram to the resulting motion of the body
Mathematically, this can be expressed as:
$$ F_i=m\ a_i\qquad,$$
where
- $F_i$ is the net force acting on the object in each direction $i$,
- $m$ is the mass of the object, and
- $a_i$ is its acceleration in the corresponding direction.
In the case of force equilibrium, the net force $\displaystyle F_i = 0$, resulting in no acceleration of the object. The object stands still or moves at a constant velocity. (Newton's first law)
The acceleration leads to the state change of an object. The velocity is updated as $\displaystyle u_i(t+ \Delta t) = u_i(t) + a_i\ \Delta t$, and the position is updated as $\displaystyle x_i(t+\Delta t) = x_i(t) + \frac{1}{2} a_i\ (\Delta t)^2$ will describe the trajectory of the object repectively.
### 2.1.2 Multi-body Analysis
Using the free body diagram as a foundation, we can break a complex system into several elements and study the dynamic interactions between interconnected bodies. The concept of [multi-body simulation](https://en.wikipedia.org/wiki/Multibody_simulation) involves applying numerical methods to model each element and evaluate the performance of the complex system.
The general steps can be followed by:
1. Free Body Diagrams provide a detailed understanding of the forces acting on __individual bodies within a system__.
2. While FBDs focus on static equilibrium, multi-body analysis extends this concept to __dynamic scenarios__, considering factors such as motion, inertia, and momentum.
3. Multi-body analysis allows engineers to __model complex systems with interconnected components__, simulating interactions such as collisions, joints, and constraints.
4. Multi-body simulations involve numerical integration methods to __solve differential equations governing the motion and state of bodies__, integrating concepts from calculus and numerical analysis.
__A Scenario of Multibody Analysis:__
A two-cart system is a simple mechanical system consisting of two carts connected by a spring and a damper. This system is often used to study concepts related to oscillations, vibrations, and control theory.
<img src="https://live.staticflickr.com/65535/53597465562_9ddbbb8769.jpg">
<br/>
<br/>
<br/>
The related components are:
- Carts: The system comprises two carts, often represented as mass points, denoted as Cart 1 and Cart 2.
- Spring: A spring connects the two carts, providing a restoring force proportional to the displacement between them.
- Damper: A damper, or damping element, is often included to dissipate energy from the system and dampen oscillations.
The dynamics of the system are described as follows:
- The motion of the two-cart system is governed by Newton's laws of motion and Hooke's law for the spring.
- The spring force exerted on each cart is proportional to the displacement between them, according to Hooke's law: $F_s=−\kappa\cdot \Delta x$, where $\kappa$ is the spring constant and $\Delta x$ is the displacement.
- The damping force exerted by the damper is typically proportional to the velocity of each cart, according to the equation: $F_d=−c\cdot \dot{x}$, where $c$ is the damping coefficient and $\dot{x}$ is the velocity.
The dynamics can be implemented in Python as:
```
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
# Define system parameters
m1 = 1.0 # Mass of first cart
m2 = 0.5 # Mass of second cart
k = 10.0 # Spring constant
c = 0.5 # Damping coefficient
def dynamics(t, state):
# Unpack state variables
x1, x1_dot, x2, x2_dot = state
# Compute accelerations
x1_ddot = (k*(x2 - x1) - c*x1_dot) / m1
x2_ddot = (-k*(x2 - x1) - c*x2_dot) / m2
return [x1_dot, x1_ddot, x2_dot, x2_ddot]
# Initial conditions
initial_state = [0.0, 0.0, 0.5, 0.0] # [x1, x1_dot, x2, x2_dot]
# Time span and integration
t_span = np.linspace(0, 10, 1000)
sol = solve_ivp(dynamics, [t_span[0], t_span[-1]], initial_state, t_eval=t_span)
# Plot results
plt.figure(figsize=(10, 6))
plt.plot(sol.t, sol.y[0], label='Cart 1')
plt.plot(sol.t, sol.y[2], label='Cart 2')
plt.xlabel('Time')
plt.ylabel('Displacement')
plt.title('Two-Cart System Simulation')
plt.legend()
plt.grid(True)
plt.show()
```
<img src="https://live.staticflickr.com/65535/53586195081_fa6128f7a7_z.jpg">
#### <ins>Scenario -- double pendulum system
A double pendulum system is one of the fundamental elemtens in robotics. Two swinging objects connected by a rod are completely driven by gravity. The trajectory of the second mass is chaotic and determined by initial conditions.
<img src="https://blogs.mathworks.com/images/pick/doublependulum_animation.gif" width="45%">
The given length for rods are 1m and 0.5m, with the second mass of 90° offset to the first mass.
__Compare the trajectories of two different parameter sets and visualize them in Python.__
- parameter I : mass 1 = 1kg, mass 2 = 3 kg
- parameter II: mass 1 = 1kg, mass 2 = 2.9 kg
[__Sample solution__](https://cocalc.com/share/public_paths/384f8ab69241a3e91fcaf08fee267697896a3356)
### 2.1.3 Structural analysis and 1-D Beam model
The concept of FBD can also be applied to infinitesimal elements. Discretizing an object into several small elements, the interaction of these local elements represents the macroscale state or dynamics. This is the main idea behind Finite Element/Volume/Difference analysis.
<img src="https://live.staticflickr.com/65535/53613770525_4b82ac60f4_z.jpg" width=80%>
<br/>
Consider an individual finite volume within a structure. The forces and moments acting on this element can be represented mathematically using the equilibrium of physical properties, such that
$$
\begin{align}
\displaystyle
m\cdot \frac{d\ \vec{v}}{d t} = &\sum_{}F\\
I \cdot \frac{d \vec{\omega}}{dt} = & \sum M
\end{align}
$$
#### __1-D Beam Model in FE Framework__
For instance, for a 1-D beam element subjected to axial, shear, and bending loads, the equilibrium states can be written as:
<img src="https://live.staticflickr.com/65535/53612463917_16d28afbc7_z.jpg" width=80%>
<br/>
<br/>
- Axial Equilibrium: $\vec{F}(x+\Delta x)\ +\ \vec{F}(x)=0$
- Shear Equilibrium: $\vec{V}(x+\Delta x)\ +\ \vec{V}(x)=0$
- Bending Moment Equilibrium: $\vec{M}(x+\Delta x) \ + \ \vec{M}(x)=0$
Where:
- $\vec{F}$, $\vec{V}$, and $\vec{M}$ are the axial force, shear force, and bending moment on the selected element
- $x$ and $(x+\Delta x)$ are the boundary positions of this selected element.
These equations serve as the basis for analyzing the forces and moments within each finite element, allowing engineers to derive the governing equations for structural behavior.
The mathematical description of such beams can be derived using the [Euler-Bernoulli beam theory](https://en.wikipedia.org/wiki/Euler%E2%80%93Bernoulli_beam_theory), which assumes small deformations and neglects shear deformation. The __Euler-Bernoulli beam equation__ describes the relationship between bending moment, curvature, and material properties of the beam. It can be derived from the equilibrium of forces and moments acting on an infinitesimal element of the beam, resulting in the differential equation:
$$ \frac{d^2}{dx^2}(EI\frac{dw^2}{dx^2})=q(x)$$
Where:
- $w(x)$ is the deflection of the beam at position $x$.
- $E$ is the Young's modulus of the material.
- $I$ is the moment of inertia of the beam's cross-section.
- $q(x)$ is the applied load per unit length.
The boundary conditions of the Euler-Bernoulli beam equation depend on the physical constraints. A summary of all boundary conditions is listed below:
| Boundary Condition | Displacement ($w$) | Rotation ($\frac{dw}{dx}$) | Bending Moment ($M = EI·\frac{d^2w}{dx^2}$) | Shear Force ($V = -EI·\frac{d^3w}{dx^3}$) |
|--------------------|------------------|------------------|----------------------------------|-------------------------------|
| **Fixed** | `w = 0` | `dw/dx = 0` | Non-zero | Non-zero |
| **Pinned** | `w = 0` | Unconstrained | `M = 0` | Non-zero |
| **Free** | Unconstrained | Unconstrained | `M = 0` | `V = 0` |
| **Guided/Roller** | Unconstrained | `dw/dx = 0` | Non-zero | `V = 0` |
Numerically, the fourth-order derivative of the variable ww can be approximated as
$$\frac{dw^4}{dx^4} \approx \frac{1}{\Delta x^4}\left( w_{i-2}-4\ w_{i-1} + 6\ w_i - 4 \ w_{i+1} + w_{i+2}\right)$$
The same idea can be applied to the boundary conditions. Thus
$$
\begin{align}
& \frac{dw}{dx} \mid_{x_0}= \frac{1}{\Delta x}\left( w_{1} - w_0\right) \qquad && \frac{dw}{dx} \mid_{x_n}= \frac{1}{\Delta x}\left( w_{n-1} - w_n\right)\\
\\
& \frac{d^2w}{dx^2} \mid_{x_1}= \frac{1}{\Delta x^2}\left( w_2 - 2w_{1} + w_0\right) \qquad && \frac{d^2w}{dx^2} \mid_{x_{n-1}}= \frac{1}{\Delta x^2}\left( w_{n-2}-2w_{n-1} - w_n\right)\\
\\
& \frac{d^3w}{dx^3} \mid_{x_0}= \frac{1}{\Delta x^3}\left( −w_0+3w_1−3w_2+w_3\right) && \frac{d^3w}{dx^3} \mid_{x_n}= \frac{1}{\Delta x^3}\left( -w_{n-4}+2w_{n-3}-2w_{n-1} + w_n\right)
\end{align}
$$
##### __A Scenario of 1D Beam Analysis:__
Consider a simply supported beam of length $L$ subjected to a uniform load $q(x)$. We aim to simulate the deflection of the beam using the finite element method and solve the Euler-Bernoulli beam equation numerically.
<img src="https://live.staticflickr.com/65535/54458293400_14638b7ebe_h.jpg" width= 60%>
<br/>
<br/>
The corresponding boundary conditions for a simply supported (pinned) beam are:
$\qquad \displaystyle w(0) = 0$ (no displacement), $\qquad w^{\prime\prime}(0) = 0$ (no moment)
as well as for the other boundary side.
$\qquad \displaystyle w(l) = w^{\prime\prime}(l) = 0$
```
import numpy as np
import matplotlib.pyplot as plt
# Constants
L = 10.0 # Length of the beam (m)
E = 2.1e11 # Young's modulus of the material (Pa)
I = 8.33e-6 # Moment of inertia of the beam's cross-section (m^4)
q = -1000.0 # Uniform load (N/m)
# Discretization
n_points = 100 # Number of points for discretization
dx = L / (n_points - 1) # Discretization step
x = np.linspace(0, L, n_points) # Discretized positions along the beam
# Construct coefficient matrix
A = np.zeros((n_points, n_points))
A[0, 0] = 1 # Fixed end boundary condition: w(0) = 0
A[1, 0:3] = np.array([1, -2, 1]) / dx**4
A[-2, -4:-1] = np.array([1, -2, 1]) / dx**4
A[-1, -1] = 1 # Free end boundary condition: dw/dx(L) = 0
for i in range(2, n_points - 2):
A[i, i-2:i+3] = np.array([1, -4, 6, -4, 1]) / dx**4
# Construct right-hand side vector
b = np.zeros(n_points)
b[2:-3] = q / (E * I)
# Solve linear system
w = np.linalg.solve(A, b)
# Plot deflection
plt.plot(x, w)
plt.xlabel('Position along the beam (m)')
plt.ylabel('Deflection (m)')
plt.title('Deflection of Cantilevered Beam under Uniform Load')
plt.grid(True)
plt.show()
```
<img src="https://live.staticflickr.com/65535/54458995003_2ec02e2d65_o.png">
#### <ins>Scenario -- non-uniform load on cantilever
Switching from a beam supported on both sides to a cantilever, where one side is free without any support. The fixed side is attached to a wall, where no displacement or rotation is allowed.
<img src="https://live.staticflickr.com/65535/54458949794_afebdc4b19_o.png" width="70%">
What is the deflection if only the outer part is loaded?
[__Sample solution__](https://cocalc.com/share/public_paths/bab139350fd29803677b376bd51979baded50ba8)
#### __Dynamic 1-D Beam Model in the FE Framework__
Extending the static condition, the equation of motion for the beam in a dynamic state can be written as:
$\displaystyle \qquad \qquad \displaystyle ρA\frac{\partial^2 w}{\partial t^2}−EI\frac{\partial^4 w}{\partial x^4}=q(x,t)$
Where:
- $w(x,t)$ is the displacement of the beam as a function of position $x$ and time $t$.
- $\rho$ is the mass density of the beam material.
- $A$ is the cross-sectional area of the beam.
- $E$ is the Young's modulus of the beam material.
- $I$ is the moment of inertia of the beam's cross-sectional area.
- $q(x,t)$ is the external force per unit length applied to the beam as a function of position and time.
Applying an external force $F(x,t)$ only on the left half section of a cantilever beam ($x \in R, \quad L/2 \leq x \leq L$ ) with the sinusiod load $F$:
$\displaystyle \qquad \qquad \qquad F(x,t)=F_0\ \sin(\omega \ t)$
Where:
- $F_0$ is the amplitude of the external force.
- $\omega$ is the angular frequency of the external force.
The finite difference method is used to discretize the beam equation in both space and time. By discretizing the spatial and temporal derivatives, we obtain a system of equations that can be solved numerically. In the code, the second-order central differences are used for the spatial derivatives, and the first-order forward difference is used for the time derivative.
The resulting discrete equation for the acceleration $a(x,t)$ at each spatial grid point and time step can be derived from the Euler-Bernoulli beam equation and the applied external force. For example, the term representing the acceleration due to the external force can be discretized as follows:
$\qquad \qquad \displaystyle \frac{F(x_i,t_n)}{\rho\ A}=\frac{F_i^n}{\rho \ A}=\frac{F_0\sin(\omega\ t_n)}{\rho\ A}$
Where $F_i^n$ represents the external force applied at spatial grid point $x_i$ and time step $t_n$.
```
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import time
rho = 100.0 # Mass density
A = 0.05 # Cross-sectional area
L = 1.0 # Length of the beam (m)
E = 2.1e8 # Young's modulus of the material (Pa)
I = 8.33e-6 # Moment of inertia of the beam's cross-section (m^4)
#q = -100.0 # Uniform load (N/m)
# Discretization
n_points = 21 # Number of points for discretization
dx = L / (n_points - 1) # Discretization step
x = np.linspace(0, L, n_points) # Discretized positions along the beam
# Numerical parameters
endTime = 4e-2 # Duration of simulation
timeStep = 4e-5
number_of_iteration = int(endTime/timeStep)
# Load definition
F0 = -100 # Amplitude of external force
omega = 1e2 * np.pi # Frequency of external force
timeSequence = np.arange(0, endTime, timeStep)
ForceSequence = F0* np.sin(omega * timeSequence)
# Stability condition (Courant condition)
C = E*I/(rho*A) * timeStep**2/dx**4
print(f"Courant number: {C}")
if C > 30:
print(f"C= {C}")
raise ValueError("Stability condition not satisfied: c * dt / dx > 1")
# Initialize displacement
w = np.zeros(n_points)
w_prev = np.zeros(n_points)
w_next = np.zeros(n_points)
# Initialize force array
q = np.zeros(n_points)
# Set up the figure and axis for animation
fig, (ax0, ax1) = plt.subplots(2, 1, figsize=(10, 9))
ax0.set_title('Dynamic 1D Beam simulation on fixed end cantilever')
ax0.set_xlabel('time [s]')
ax0.set_ylabel('loading force [Nt/m]')
ax0.plot(timeSequence, ForceSequence, "--", color="black")
point, = ax0.plot([], 'o', color='red', markersize=10)
ax0.set_ylim(-100, 100)
ax0.grid()
#ax1.plot()
time_template = 'time = %.4fs'
time_text = ax0.text(0.05, 0.9, '', transform=ax0.transAxes)
line, = ax1.plot(x, w)
ax1.set_xlabel('Position x')
ax1.set_ylabel('Displacement w [mm]')
ax1.set_ylim(-5, 5)
ax1.grid()
time_template = 'time = %.4fs'
time_text = ax1.text(0.05, 0.9, '', transform=ax1.transAxes)
# Function to update the plot for each frame
def update(frame):
global w, w_prev, w_next
currentTime = frame * timeStep
# External force applied at the middle section
q[int(0.5 * n_points):] = F0 * np.sin(omega * currentTime)
# Compute the next time step
for i in range(2, n_points-2):
w_next[i] = 2*w[i] - w_prev[i] \
- (timeStep**2.0 / (rho * A)) * ((E * I) / dx**4) * (w[i+2] - 4*w[i+1] + 6*w[i] -4* w[i-1] + w[i-2]) \
+ timeStep**2.0 * q[i] / (rho * A)
# boundary condition for fixed point
w_next[0] = 0 # w = 0
w_next[1] = w_next[0] # dw/dx=0 at x = 0
# boundary condition for free point
w_next[-2] = 2* w_next[-3] - w_next[-4] # d2w/dx2 = 0 Use forward difference:
w_next[-1] = 2* w_next[-2] - w_next[-3] # d3w/dx3 = 0 Use forward difference:
#w_next[-1] = 0
#w_next[-2] = w_next[-1]
# Update previous and current time steps
w_prev = w.copy()
w = w_next.copy()
if frame % 10 == 1:
time_text.set_text(time_template % currentTime)
# Update the plot
line.set_ydata(w*1e3)
point.set_data([[currentTime], [F0 * np.sin(omega * currentTime)]])
return line, time_text, point
```
<img src="https://live.staticflickr.com/65535/54464844863_7f0946bb63_o.gif">
<br/>
<br/>
So far, the modeling is carried out in the one-dimensional regime, both in stationary and transient types. The implementation into scripts can be easily carried out from scratch. For analysis in higher dimension, the complexity increases exponentially. Simple scripts are not easy to implement and it is also computational inefficient. For such cases, simulation software for finite element calculations is rather recommended.
---
[__BACK TO CONTENT__](https://hackmd.io/@SamuelChang/H1LvI_eYn)
Sign in with Wallet
Connect another wallet