Python Project: Beam Analyzer
In this capstone lesson, we'll build a comprehensive beam analysis tool that combines everything we've learned. This tool will calculate reactions, generate SFD/BMD diagrams, compute deflections, and check stress limits.
Project Structure
beam_analyzer/
├── beam.py # Main Beam class
├── sections.py # Cross-section definitions
├── materials.py # Material properties
├── loads.py # Load types
├── solver.py # Analysis engine
└── visualization.py # Plotting functionsComplete Implementation
import numpy as np
import matplotlib.pyplot as plt
from dataclasses import dataclass, field
from typing import List, Optional, Callable
from enum import Enum
# ============================================================
# MATERIALS MODULE
# ============================================================
@dataclass
class Material:
"""Material properties."""
name: str
E: float # Young's modulus (GPa)
G: float # Shear modulus (GPa)
Sy: float # Yield strength (MPa)
density: float # kg/m³
@classmethod
def steel(cls):
return cls("Steel", E=200, G=77, Sy=250, density=7850)
@classmethod
def aluminum(cls):
return cls("Aluminum 6061-T6", E=70, G=26, Sy=275, density=2700)
@classmethod
def stainless_steel(cls):
return cls("Stainless Steel 304", E=193, G=77, Sy=215, density=8000)
# ============================================================
# SECTIONS MODULE
# ============================================================
@dataclass
class Section:
"""Cross-section properties."""
name: str
A: float # Area (mm²)
Ix: float # Moment of inertia about x-axis (mm⁴)
Iy: float # Moment of inertia about y-axis (mm⁴)
Zx: float # Section modulus x (mm³)
Zy: float # Section modulus y (mm³)
c_top: float # Distance to top fiber (mm)
c_bot: float # Distance to bottom fiber (mm)
@classmethod
def rectangle(cls, b: float, h: float):
"""Create rectangular section."""
A = b * h
Ix = b * h**3 / 12
Iy = h * b**3 / 12
c = h / 2
Zx = Ix / c
Zy = Iy / (b/2)
return cls(f"Rectangle {b}x{h}", A, Ix, Iy, Zx, Zy, c, c)
@classmethod
def circular(cls, d: float):
"""Create circular section."""
A = np.pi * d**2 / 4
I = np.pi * d**4 / 64
c = d / 2
Z = I / c
return cls(f"Circle Ø{d}", A, I, I, Z, Z, c, c)
@classmethod
def i_section(cls, bf: float, tf: float, hw: float, tw: float):
"""Create I-section (wide flange)."""
h = hw + 2 * tf
A = 2 * bf * tf + hw * tw
Ix = (bf * h**3 - (bf - tw) * hw**3) / 12
Iy = (2 * tf * bf**3 + hw * tw**3) / 12
c = h / 2
Zx = Ix / c
Zy = Iy / (bf / 2)
return cls(f"I-Section {bf}x{h}", A, Ix, Iy, Zx, Zy, c, c)
# ============================================================
# LOADS MODULE
# ============================================================
class LoadType(Enum):
POINT = "point"
DISTRIBUTED = "distributed"
MOMENT = "moment"
@dataclass
class Load:
"""Base load class."""
load_type: LoadType
position: float # mm from left end
@dataclass
class PointLoad(Load):
"""Point load (force)."""
magnitude: float # N (positive = downward)
def __init__(self, position: float, magnitude: float):
super().__init__(LoadType.POINT, position)
self.magnitude = magnitude
@dataclass
class DistributedLoad(Load):
"""Uniformly or linearly distributed load."""
start: float # mm
end: float # mm
w_start: float # N/mm at start
w_end: float # N/mm at end
def __init__(self, start: float, end: float, w_start: float, w_end: float = None):
super().__init__(LoadType.DISTRIBUTED, start)
self.start = start
self.end = end
self.w_start = w_start
self.w_end = w_end if w_end is not None else w_start
@dataclass
class PointMoment(Load):
"""Applied moment."""
magnitude: float # N·mm (positive = CCW)
def __init__(self, position: float, magnitude: float):
super().__init__(LoadType.MOMENT, position)
self.magnitude = magnitude
# ============================================================
# MAIN BEAM CLASS
# ============================================================
class Beam:
"""Complete beam analysis class."""
def __init__(self, length: float, support_type: str = "simply_supported"):
"""
Initialize beam.
Parameters:
-----------
length : float - Beam length in mm
support_type : str - 'simply_supported', 'cantilever', 'fixed_fixed',
'propped_cantilever'
"""
self.L = length
self.support_type = support_type
self.material: Optional[Material] = None
self.section: Optional[Section] = None
self.loads: List[Load] = []
self.reactions = {}
self._solved = False
def set_material(self, material: Material):
"""Set beam material."""
self.material = material
self._solved = False
return self
def set_section(self, section: Section):
"""Set beam cross-section."""
self.section = section
self._solved = False
return self
def add_point_load(self, position: float, magnitude: float):
"""Add point load (kN, converted internally to N)."""
self.loads.append(PointLoad(position, magnitude * 1000))
self._solved = False
return self
def add_distributed_load(self, start: float, end: float,
w_start: float, w_end: float = None):
"""Add distributed load (N/mm)."""
self.loads.append(DistributedLoad(start, end, w_start, w_end))
self._solved = False
return self
def add_moment(self, position: float, magnitude: float):
"""Add moment (kN·m, converted to N·mm)."""
self.loads.append(PointMoment(position, magnitude * 1e6))
self._solved = False
return self
def solve(self):
"""Calculate reactions."""
if self._solved:
return self.reactions
# Sum of forces and moments from loads
total_force = 0
total_moment_A = 0 # About left support
for load in self.loads:
if isinstance(load, PointLoad):
total_force += load.magnitude
total_moment_A += load.magnitude * load.position
elif isinstance(load, DistributedLoad):
length = load.end - load.start
avg_w = (load.w_start + load.w_end) / 2
W = avg_w * length
total_force += W
# Centroid
if load.w_start == load.w_end:
centroid = load.start + length / 2
else:
centroid = load.start + length/3 * (load.w_start + 2*load.w_end) / (load.w_start + load.w_end)
total_moment_A += W * centroid
elif isinstance(load, PointMoment):
total_moment_A += load.magnitude
# Calculate reactions based on support type
if self.support_type == "simply_supported":
R_B = total_moment_A / self.L
R_A = total_force - R_B
self.reactions = {'R_A': R_A, 'R_B': R_B, 'M_A': 0, 'M_B': 0}
elif self.support_type == "cantilever":
R_A = total_force
M_A = -total_moment_A
self.reactions = {'R_A': R_A, 'M_A': M_A}
elif self.support_type == "fixed_fixed":
# Simplified: assume UDL for fixed-fixed
# More complex analysis needed for general case
R_A = total_force / 2
R_B = total_force / 2
M_A = -total_moment_A / 2 + total_force * self.L / 8
M_B = total_moment_A / 2 - total_force * self.L / 8
self.reactions = {'R_A': R_A, 'R_B': R_B, 'M_A': M_A, 'M_B': M_B}
self._solved = True
return self.reactions
def shear_at(self, x: float) -> float:
"""Calculate shear force at position x."""
if not self._solved:
self.solve()
V = 0
# Reaction at left
if x > 0:
V += self.reactions.get('R_A', 0)
# Point loads
for load in self.loads:
if isinstance(load, PointLoad) and load.position < x:
V -= load.magnitude
elif isinstance(load, DistributedLoad) and x > load.start:
x_eff = min(x, load.end)
length = x_eff - load.start
if load.w_start == load.w_end:
V -= load.w_start * length
else:
ratio = length / (load.end - load.start)
w_at_x = load.w_start + ratio * (load.w_end - load.w_start)
avg_w = (load.w_start + w_at_x) / 2
V -= avg_w * length
return V
def moment_at(self, x: float) -> float:
"""Calculate bending moment at position x."""
if not self._solved:
self.solve()
M = 0
# Fixed-end moment
M += self.reactions.get('M_A', 0)
# Reaction moment
if x > 0:
M += self.reactions.get('R_A', 0) * x
# Loads
for load in self.loads:
if isinstance(load, PointLoad) and load.position < x:
M -= load.magnitude * (x - load.position)
elif isinstance(load, PointMoment) and load.position < x:
M += load.magnitude
elif isinstance(load, DistributedLoad) and x > load.start:
x_eff = min(x, load.end)
length = x_eff - load.start
if load.w_start == load.w_end:
W = load.w_start * length
centroid = load.start + length / 2
else:
ratio = length / (load.end - load.start)
w_at_x = load.w_start + ratio * (load.w_end - load.w_start)
avg_w = (load.w_start + w_at_x) / 2
W = avg_w * length
centroid = load.start + length / 2
M -= W * (x - centroid)
return M
def deflection_at(self, x: float) -> float:
"""
Calculate deflection using numerical integration.
"""
if self.material is None or self.section is None:
raise ValueError("Material and section must be set")
E = self.material.E * 1e3 # MPa
I = self.section.Ix
# Double integration of M/(EI) for deflection
# Use numerical integration
n = 200
xi = np.linspace(0, x, n)
M_vals = np.array([self.moment_at(xii) for xii in xi])
# First integration: slope
curvature = M_vals / (E * I)
slope = np.trapz(curvature, xi)
# This is a simplified approach; full analysis requires
# boundary condition handling
dx = x / (n - 1) if n > 1 else 0
# For simply supported, integrate twice
if self.support_type == "simply_supported":
# Use conjugate beam or direct formula
# Simplified: δ ≈ M_max * L² / (10 * E * I) at center
idx_max = n // 2
M_center = self.moment_at(self.L / 2)
return M_center * self.L**2 / (10 * E * I) * (4 * x / self.L - 4 * (x/self.L)**2)
elif self.support_type == "cantilever":
# δ(x) from double integration
xi_full = np.linspace(0, self.L, n)
M_full = np.array([self.moment_at(xii) for xii in xi_full])
curv_full = M_full / (E * I)
# Integrate twice from fixed end
slope_arr = np.zeros(n)
defl_arr = np.zeros(n)
dx = self.L / (n - 1)
for i in range(1, n):
slope_arr[i] = slope_arr[i-1] + curv_full[i-1] * dx
defl_arr[i] = defl_arr[i-1] + slope_arr[i-1] * dx
# Find deflection at x
idx = int(x / self.L * (n - 1))
return defl_arr[min(idx, n-1)]
return 0
def max_stress(self) -> dict:
"""Calculate maximum bending stress."""
if self.section is None:
raise ValueError("Section must be set")
n = 100
x = np.linspace(0, self.L, n)
M = np.array([self.moment_at(xi) for xi in x])
M_max_pos = np.max(M)
M_max_neg = np.min(M)
M_max = M_max_pos if abs(M_max_pos) > abs(M_max_neg) else M_max_neg
# σ = M * c / I
sigma_top = M_max * self.section.c_top / self.section.Ix
sigma_bot = M_max * self.section.c_bot / self.section.Ix
return {
'M_max_Nmm': M_max,
'M_max_kNm': M_max / 1e6,
'sigma_top_MPa': sigma_top,
'sigma_bot_MPa': sigma_bot,
'sigma_max_MPa': max(abs(sigma_top), abs(sigma_bot))
}
def check_design(self, FoS: float = 2.0) -> dict:
"""Check if beam design is adequate."""
if self.material is None:
raise ValueError("Material must be set")
stress = self.max_stress()
sigma_allow = self.material.Sy / FoS
return {
'sigma_max_MPa': stress['sigma_max_MPa'],
'sigma_allow_MPa': sigma_allow,
'FoS_actual': self.material.Sy / stress['sigma_max_MPa'],
'is_safe': stress['sigma_max_MPa'] <= sigma_allow
}
def plot(self, show_deflection: bool = True, num_points: int = 200):
"""Generate complete beam analysis plots."""
if not self._solved:
self.solve()
x = np.linspace(0, self.L, num_points)
V = np.array([self.shear_at(xi) for xi in x])
M = np.array([self.moment_at(xi) for xi in x])
n_plots = 4 if show_deflection and self.material and self.section else 3
fig, axes = plt.subplots(n_plots, 1, figsize=(12, 3*n_plots))
# 1. Loading diagram
ax1 = axes[0]
ax1.axhline(y=0, color='black', lw=3)
# Supports
if self.support_type == "simply_supported":
ax1.scatter([0, self.L], [0, 0], s=300, marker='^', color='green', zorder=5)
elif self.support_type == "cantilever":
ax1.fill_between([-self.L*0.02, 0], [-20, -20], [20, 20], color='gray', alpha=0.5)
# Loads
for load in self.loads:
if isinstance(load, PointLoad):
sign = 1 if load.magnitude > 0 else -1
ax1.annotate('', xy=(load.position, 0),
xytext=(load.position, sign * 40),
arrowprops=dict(arrowstyle='->', lw=2, color='red'))
ax1.text(load.position, sign * 45,
f'{abs(load.magnitude)/1000:.1f}kN',
ha='center', fontsize=9, color='red')
elif isinstance(load, DistributedLoad):
xs = np.linspace(load.start, load.end, 15)
for xi in xs:
ax1.annotate('', xy=(xi, 0), xytext=(xi, 25),
arrowprops=dict(arrowstyle='->', lw=1, color='blue', alpha=0.7))
# Reactions
ax1.annotate('', xy=(0, 0), xytext=(0, -35),
arrowprops=dict(arrowstyle='->', lw=2, color='green'))
ax1.text(0, -45, f"R_A={self.reactions.get('R_A', 0)/1000:.2f}kN",
ha='center', fontsize=9, color='green')
if self.support_type == "simply_supported":
ax1.annotate('', xy=(self.L, 0), xytext=(self.L, -35),
arrowprops=dict(arrowstyle='->', lw=2, color='green'))
ax1.text(self.L, -45, f"R_B={self.reactions.get('R_B', 0)/1000:.2f}kN",
ha='center', fontsize=9, color='green')
ax1.set_xlim(-self.L*0.05, self.L*1.05)
ax1.set_ylim(-60, 60)
ax1.set_title('Loading Diagram', fontsize=11)
ax1.axis('off')
# 2. Shear Force Diagram
ax2 = axes[1]
ax2.fill_between(x, 0, V/1000, where=(V >= 0), alpha=0.3, color='blue')
ax2.fill_between(x, 0, V/1000, where=(V < 0), alpha=0.3, color='red')
ax2.plot(x, V/1000, 'k-', lw=1.5)
ax2.axhline(y=0, color='k', lw=0.5)
ax2.set_ylabel('V (kN)')
ax2.set_title(f'Shear Force Diagram | Vmax = {np.max(np.abs(V))/1000:.2f} kN')
ax2.grid(True, alpha=0.3)
# 3. Bending Moment Diagram
ax3 = axes[2]
ax3.fill_between(x, 0, M/1e6, where=(M >= 0), alpha=0.3, color='green')
ax3.fill_between(x, 0, M/1e6, where=(M < 0), alpha=0.3, color='orange')
ax3.plot(x, M/1e6, 'k-', lw=1.5)
ax3.axhline(y=0, color='k', lw=0.5)
ax3.set_ylabel('M (kN·m)')
ax3.set_title(f'Bending Moment Diagram | Mmax = {np.max(np.abs(M))/1e6:.2f} kN·m')
ax3.grid(True, alpha=0.3)
# 4. Deflection (if material and section are set)
if show_deflection and self.material and self.section:
ax4 = axes[3]
y = np.array([self.deflection_at(xi) for xi in x])
ax4.fill_between(x, 0, -y, alpha=0.3, color='purple')
ax4.plot(x, -y, 'purple', lw=1.5)
ax4.axhline(y=0, color='k', lw=0.5)
ax4.set_xlabel('Position (mm)')
ax4.set_ylabel('δ (mm)')
ax4.set_title(f'Deflection | δmax = {np.max(np.abs(y)):.3f} mm')
ax4.grid(True, alpha=0.3)
else:
axes[-1].set_xlabel('Position (mm)')
plt.tight_layout()
plt.show()
return {'V': V, 'M': M, 'x': x}
def summary(self):
"""Print analysis summary."""
if not self._solved:
self.solve()
print("=" * 60)
print("BEAM ANALYSIS SUMMARY")
print("=" * 60)
print(f"\nGeometry:")
print(f" Length: {self.L} mm ({self.L/1000:.2f} m)")
print(f" Support: {self.support_type}")
if self.material:
print(f"\nMaterial: {self.material.name}")
print(f" E = {self.material.E} GPa")
print(f" Sy = {self.material.Sy} MPa")
if self.section:
print(f"\nSection: {self.section.name}")
print(f" A = {self.section.A:.1f} mm²")
print(f" Ix = {self.section.Ix:.2e} mm⁴")
print(f" Zx = {self.section.Zx:.2e} mm³")
print(f"\nReactions:")
for key, value in self.reactions.items():
if 'M' in key:
print(f" {key} = {value/1e6:.3f} kN·m")
else:
print(f" {key} = {value/1000:.3f} kN")
if self.section:
stress = self.max_stress()
print(f"\nStress Analysis:")
print(f" Mmax = {stress['M_max_kNm']:.3f} kN·m")
print(f" σmax = {stress['sigma_max_MPa']:.2f} MPa")
if self.material:
check = self.check_design()
print(f" FoS = {check['FoS_actual']:.2f}")
print(f" Status: {'SAFE ✓' if check['is_safe'] else 'UNSAFE ✗'}")
print("=" * 60)
# ============================================================
# EXAMPLE USAGE
# ============================================================
if __name__ == "__main__":
# Create a beam
beam = Beam(length=6000, support_type="simply_supported")
# Set properties
beam.set_material(Material.steel())
beam.set_section(Section.i_section(bf=150, tf=10, hw=280, tw=8))
# Add loads
beam.add_point_load(2000, 25) # 25 kN at 2m
beam.add_point_load(4000, 15) # 15 kN at 4m
beam.add_distributed_load(0, 6000, 5, 5) # 5 N/mm UDL
# Analyze
beam.summary()
beam.plot()
