Beam Bending Theory

Beams are structural elements that carry loads perpendicular to their axis. Understanding beam bending is crucial for designing everything from bridges to aircraft wings.

The Bending Problem

When a beam bends:

• The top fibers compress (shorter)
• The bottom fibers stretch (longer)
• The neutral axis remains unstressed

The Flexure Formula

The bending stress at distance $y$ from the neutral axis:

$$\sigma = \frac{My}{I}$$

Where:

• $\sigma$ = Bending stress (Pa)
• $M$ = Bending moment (N·m)
• $y$ = Distance from neutral axis (m)
• $I$ = Moment of inertia (m⁴)

$$\sigma_{max} = \frac{Mc}{I} = \frac{M}{Z}$$

Where $c$ is the distance to the extreme fiber and $Z = I/c$ is the section modulus.

Moment of Inertia

The moment of inertia quantifies a cross-section's resistance to bending:

$$I = \int y^2 \, dA$$

import numpy as np

def rectangular_section(b, h):
    """
    Properties of rectangular cross-section.

    Parameters:
    -----------
    b : float - Width (mm)
    h : float - Height (mm)

    Returns:
    --------
    dict with area, I, Z, and centroid
    """
    A = b * h
    I = (b * h**3) / 12
    c = h / 2
    Z = I / c

    return {
        'area_mm2': A,
        'I_mm4': I,
        'c_mm': c,
        'Z_mm3': Z,
        'description': f'Rectangle {b}×{h} mm'
    }

def circular_section(d):
    """Properties of circular cross-section."""
    r = d / 2
    A = np.pi * r**2
    I = (np.pi * d**4) / 64
    c = d / 2
    Z = I / c

    return {
        'area_mm2': A,
        'I_mm4': I,
        'c_mm': c,
        'Z_mm3': Z,
        'description': f'Circle Ø{d} mm'
    }

def hollow_circular_section(d_outer, d_inner):
    """Properties of hollow circular (tube) cross-section."""
    A = np.pi * (d_outer**2 - d_inner**2) / 4
    I = np.pi * (d_outer**4 - d_inner**4) / 64
    c = d_outer / 2
    Z = I / c

    return {
        'area_mm2': A,
        'I_mm4': I,
        'c_mm': c,
        'Z_mm3': Z,
        'description': f'Tube Ø{d_outer}/Ø{d_inner} mm'
    }

def i_section(b_f, t_f, h_w, t_w):
    """
    Properties of I-section (wide flange).

    Parameters:
    -----------
    b_f : float - Flange width (mm)
    t_f : float - Flange thickness (mm)
    h_w : float - Web height (mm)
    t_w : float - Web thickness (mm)
    """
    h_total = h_w + 2 * t_f

    # Area
    A = 2 * b_f * t_f + h_w * t_w

    # Moment of inertia (about horizontal centroidal axis)
    # Two flanges + web
    I_flanges = 2 * (b_f * t_f**3 / 12 + b_f * t_f * ((h_w + t_f) / 2)**2)
    I_web = t_w * h_w**3 / 12
    I = I_flanges + I_web

    c = h_total / 2
    Z = I / c

    return {
        'area_mm2': A,
        'I_mm4': I,
        'c_mm': c,
        'Z_mm3': Z,
        'h_total_mm': h_total,
        'description': f'I-section: flanges {b_f}×{t_f}, web {h_w}×{t_w} mm'
    }

# Compare sections
print("Cross-Section Properties Comparison")
print("=" * 60)

sections = [
    rectangular_section(50, 100),
    circular_section(60),
    hollow_circular_section(60, 50),
    i_section(b_f=100, t_f=10, h_w=180, t_w=8)
]

for sec in sections:
    print(f"\n{sec['description']}:")
    print(f"  Area: {sec['area_mm2']:.1f} mm²")
    print(f"  I: {sec['I_mm4']:.2e} mm⁴")
    print(f"  Z: {sec['Z_mm3']:.2e} mm³")

Bending Stress Analysis

def bending_stress_analysis(M_kNm, section_props, material_yield_MPa):
    """
    Analyze bending stress in a beam section.

    Parameters:
    -----------
    M_kNm : float - Bending moment in kN·m
    section_props : dict - From section property functions
    material_yield_MPa : float - Yield strength

    Returns:
    --------
    dict with stress analysis results
    """
    M = M_kNm * 1e6  # Convert to N·mm

    I = section_props['I_mm4']
    c = section_props['c_mm']

    # Maximum bending stress
    sigma_max = M * c / I  # MPa (since M in N·mm, c in mm, I in mm⁴)

    # Factor of safety
    FoS = material_yield_MPa / abs(sigma_max)

    return {
        'sigma_max_MPa': sigma_max,
        'sigma_tension_MPa': sigma_max if M > 0 else -sigma_max,
        'sigma_compression_MPa': -sigma_max if M > 0 else sigma_max,
        'factor_of_safety': FoS,
        'is_safe': FoS >= 1.5
    }

# Example: Cantilever beam with point load
# Simply supported beam with center load
span = 4000  # mm
load = 20    # kN
M_max = load * span / 4 / 1000  # kN·m for simply supported with center load

section = i_section(b_f=150, t_f=12, h_w=276, t_w=8)

print(f"Beam Analysis: Simply Supported, Center Point Load")
print(f"Span: {span} mm, Load: {load} kN")
print(f"Maximum moment: {M_max:.2f} kN·m")
print(f"\nSection: {section['description']}")

result = bending_stress_analysis(M_max, section, material_yield_MPa=250)

print(f"\nResults:")
print(f"  Max bending stress: {result['sigma_max_MPa']:.2f} MPa")
print(f"  Tension (bottom): {result['sigma_tension_MPa']:.2f} MPa")
print(f"  Compression (top): {result['sigma_compression_MPa']:.2f} MPa")
print(f"  Factor of Safety: {result['factor_of_safety']:.2f}")
print(f"  Safe design: {'Yes' if result['is_safe'] else 'No - Increase section!'}")

Stress Distribution Visualization

import matplotlib.pyplot as plt

def visualize_bending_stress(b, h, M_kNm, y_points=50):
    """
    Visualize bending stress distribution in a rectangular beam.
    """
    # Section properties
    I = b * h**3 / 12
    c = h / 2

    # Stress at various y positions
    y = np.linspace(-c, c, y_points)
    M = M_kNm * 1e6  # N·mm
    sigma = M * y / I

    # Create visualization
    fig, axes = plt.subplots(1, 3, figsize=(14, 6))

    # 1. Cross-section
    ax1 = axes[0]
    rect = plt.Rectangle((-b/2, -h/2), b, h,
                         facecolor='lightblue', edgecolor='black', lw=2)
    ax1.add_patch(rect)
    ax1.axhline(y=0, color='red', linestyle='--', lw=2, label='Neutral Axis')
    ax1.set_xlim(-b, b)
    ax1.set_ylim(-h*0.7, h*0.7)
    ax1.set_aspect('equal')
    ax1.set_xlabel('Width (mm)')
    ax1.set_ylabel('Height (mm)')
    ax1.set_title('Cross-Section')
    ax1.legend()
    ax1.grid(True, alpha=0.3)

    # 2. Stress distribution
    ax2 = axes[1]
    ax2.fill_betweenx(y, 0, sigma, alpha=0.3,
                      color='blue' if M_kNm > 0 else 'red')
    ax2.plot(sigma, y, 'b-' if M_kNm > 0 else 'r-', lw=2)
    ax2.axvline(x=0, color='k', lw=0.5)
    ax2.axhline(y=0, color='red', linestyle='--', lw=2)

    # Annotate max stresses
    sigma_max = max(abs(sigma))
    ax2.annotate(f'σmax = {max(sigma):.1f} MPa\n(Tension)',
                xy=(max(sigma), c), xytext=(max(sigma)*1.1, c*0.8),
                fontsize=10, color='blue')
    ax2.annotate(f'σmin = {min(sigma):.1f} MPa\n(Compression)',
                xy=(min(sigma), -c), xytext=(min(sigma)*1.1, -c*0.8),
                fontsize=10, color='red')

    ax2.set_xlabel('Bending Stress σ (MPa)')
    ax2.set_ylabel('Distance from NA (mm)')
    ax2.set_title('Stress Distribution')
    ax2.grid(True, alpha=0.3)

    # 3. Beam deflection shape (qualitative)
    ax3 = axes[2]
    x_beam = np.linspace(0, 1000, 100)
    # Parabolic deflection for simply supported beam
    y_deflection = -10 * (x_beam/500 - (x_beam/500)**2) * 4
    ax3.plot(x_beam, y_deflection, 'b-', lw=3)
    ax3.plot(x_beam, np.zeros_like(x_beam), 'k--', lw=1, alpha=0.5)

    # Supports
    ax3.scatter([0, 1000], [0, 0], s=200, marker='^', color='green', zorder=5)

    ax3.set_xlabel('Beam Length')
    ax3.set_ylabel('Deflection')
    ax3.set_title('Deflected Shape (Exaggerated)')
    ax3.set_aspect('equal')
    ax3.set_ylim(-20, 10)

    plt.tight_layout()
    plt.show()

# Visualize
visualize_bending_stress(b=50, h=100, M_kNm=5)

Parallel Axis Theorem

When calculating I about an axis that doesn't pass through the centroid:

$$I = I_c + Ad^2$$

Where:

• $I_c$ = Moment of inertia about centroidal axis
• $A$ = Area
• $d$ = Distance from centroid to new axis

def composite_section_I(sections):
    """
    Calculate moment of inertia for a composite section.

    Parameters:
    -----------
    sections : list of dicts
        Each: {'A': area, 'I_c': centroidal I, 'y_c': centroid y-coordinate}
    """
    # Find composite centroid
    total_A = sum(s['A'] for s in sections)
    y_bar = sum(s['A'] * s['y_c'] for s in sections) / total_A

    # Calculate I using parallel axis theorem
    I_total = 0
    for s in sections:
        d = s['y_c'] - y_bar
        I_total += s['I_c'] + s['A'] * d**2

    return {
        'total_area': total_A,
        'centroid_y': y_bar,
        'I_total': I_total
    }

# Example: T-section
# Flange: 150×20 mm at top
# Web: 20×180 mm below flange

flange = {
    'A': 150 * 20,
    'I_c': 150 * 20**3 / 12,
    'y_c': 180 + 10  # Center of flange from bottom
}

web = {
    'A': 20 * 180,
    'I_c': 20 * 180**3 / 12,
    'y_c': 90  # Center of web from bottom
}

result = composite_section_I([flange, web])

print("T-Section Analysis:")
print(f"  Total area: {result['total_area']:.0f} mm²")
print(f"  Centroid from bottom: {result['centroid_y']:.2f} mm")
print(f"  Moment of inertia: {result['I_total']:.2e} mm⁴")

Section Selection for Design

def select_beam_section(M_max_kNm, sigma_allow_MPa, section_type='rectangular'):
    """
    Select minimum section size for given moment.

    Parameters:
    -----------
    M_max_kNm : float - Maximum bending moment
    sigma_allow_MPa : float - Allowable stress
    section_type : str - 'rectangular' or 'circular'
    """
    M = M_max_kNm * 1e6  # N·mm

    # Required section modulus
    Z_required = M / sigma_allow_MPa

    print(f"Design Moment: {M_max_kNm} kN·m")
    print(f"Allowable Stress: {sigma_allow_MPa} MPa")
    print(f"Required Z: {Z_required:.0f} mm³")
    print()

    if section_type == 'rectangular':
        # For b = h/2 (typical proportions)
        # Z = bh²/6 = (h/2)h²/6 = h³/12
        h_required = (12 * Z_required) ** (1/3)
        b_required = h_required / 2

        # Round up to standard sizes
        h_std = np.ceil(h_required / 10) * 10
        b_std = np.ceil(b_required / 10) * 10

        sec = rectangular_section(b_std, h_std)
        print(f"Rectangular section (b = h/2):")
        print(f"  Calculated: {b_required:.1f} × {h_required:.1f} mm")
        print(f"  Selected: {b_std:.0f} × {h_std:.0f} mm")
        print(f"  Actual Z: {sec['Z_mm3']:.0f} mm³")

    elif section_type == 'circular':
        # Z = πd³/32
        d_required = (32 * Z_required / np.pi) ** (1/3)
        d_std = np.ceil(d_required / 5) * 5

        sec = circular_section(d_std)
        print(f"Circular section:")
        print(f"  Calculated diameter: {d_required:.1f} mm")
        print(f"  Selected diameter: {d_std:.0f} mm")
        print(f"  Actual Z: {sec['Z_mm3']:.0f} mm³")

# Example
select_beam_section(M_max_kNm=15, sigma_allow_MPa=150, section_type='rectangular')

Key Takeaways

• Flexure formula: $\sigma = My/I$ — stress varies linearly with distance from neutral axis
• Maximum stress: Occurs at extreme fibers: $\sigma_{max} = Mc/I = M/Z$
• Section modulus Z = I/c — higher Z means lower stress for same moment
• I-sections are efficient because material is placed far from the neutral axis
• Parallel axis theorem: $I = I_c + Ad^2$ — for composite sections

Next lesson: We'll calculate beam deflections using integration methods.