Strain Energy Methods

Energy methods provide powerful alternative approaches to solving structural problems, especially for complex loading and indeterminate structures.

Strain Energy

When a body is loaded, work is done and stored as strain energy:

$$U = \frac{1}{2} \int \sigma \cdot \varepsilon \, dV$$

For different loading types:

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate

class StrainEnergy:
    """Calculate strain energy for various loading conditions."""

    @staticmethod
    def axial(P_N, L_mm, A_mm2, E_GPa):
        """
        Strain energy in an axially loaded member.

        U = P²L / (2AE)
        """
        E = E_GPa * 1e3  # MPa
        U = P_N**2 * L_mm / (2 * A_mm2 * E)  # N·mm = mJ
        return U

    @staticmethod
    def bending_uniform_moment(M_Nmm, L_mm, E_GPa, I_mm4):
        """
        Strain energy for uniform bending moment.

        U = M²L / (2EI)
        """
        E = E_GPa * 1e3
        U = M_Nmm**2 * L_mm / (2 * E * I_mm4)
        return U

    @staticmethod
    def bending_cantilever_end_load(P_N, L_mm, E_GPa, I_mm4):
        """
        Strain energy for cantilever with end load.

        M(x) = P(L-x), U = ∫M²/(2EI)dx = P²L³/(6EI)
        """
        E = E_GPa * 1e3
        U = P_N**2 * L_mm**3 / (6 * E * I_mm4)
        return U

    @staticmethod
    def bending_simply_supported_center_load(P_N, L_mm, E_GPa, I_mm4):
        """
        Strain energy for simply supported beam with center load.

        U = P²L³/(96EI)
        """
        E = E_GPa * 1e3
        U = P_N**2 * L_mm**3 / (96 * E * I_mm4)
        return U

    @staticmethod
    def torsion(T_Nmm, L_mm, G_GPa, J_mm4):
        """
        Strain energy in torsion.

        U = T²L / (2GJ)
        """
        G = G_GPa * 1e3
        U = T_Nmm**2 * L_mm / (2 * G * J_mm4)
        return U


# Example calculations
print("Strain Energy Examples")
print("=" * 50)

# Axial member
P = 50e3  # 50 kN
L = 2000  # mm
A = 500   # mm²
E = 200   # GPa

U_axial = StrainEnergy.axial(P, L, A, E)
print(f"\nAxial (P={P/1000}kN, L={L}mm, A={A}mm²):")
print(f"  Strain energy U = {U_axial:.2f} mJ = {U_axial/1000:.4f} J")

# Cantilever
I = 1e6   # mm⁴
U_cant = StrainEnergy.bending_cantilever_end_load(10e3, 1000, E, I)
print(f"\nCantilever (P=10kN, L=1000mm):")
print(f"  Strain energy U = {U_cant:.2f} mJ")

Castigliano's First Theorem

The partial derivative of strain energy with respect to a load gives the corresponding displacement:

$$\delta_i = \frac{\partial U}{\partial P_i}$$

For deflection at the point of load application:

$$\delta = \frac{\partial U}{\partial P}$$

def castigliano_cantilever(P_N, L_mm, E_GPa, I_mm4):
    """
    Use Castigliano's theorem to find cantilever tip deflection.

    U = P²L³/(6EI)
    δ = ∂U/∂P = PL³/(3EI)
    """
    E = E_GPa * 1e3

    # Strain energy
    U = P_N**2 * L_mm**3 / (6 * E * I_mm4)

    # Deflection (derivative of U with respect to P)
    delta = P_N * L_mm**3 / (3 * E * I_mm4)

    return {
        'U_mJ': U,
        'delta_mm': delta,
        'method': "Castigliano's Theorem"
    }

# Example
result = castigliano_cantilever(P_N=5000, L_mm=2000, E_GPa=200, I_mm4=2e6)
print(f"\nCastigliano - Cantilever Deflection:")
print(f"  Strain energy: {result['U_mJ']:.2f} mJ")
print(f"  Tip deflection: {result['delta_mm']:.4f} mm")

# Verify with standard formula
delta_standard = 5000 * 2000**3 / (3 * 200e3 * 2e6)
print(f"  Standard formula: {delta_standard:.4f} mm")

Deflection at Any Point (Dummy Load Method)

To find deflection at a point where no load exists, apply a dummy load Q at that point:

$$\delta = \frac{\partial U}{\partial Q}\bigg|_{Q=0}$$

def castigliano_midspan_deflection():
    """
    Find midspan deflection of cantilever using dummy load method.

    Cantilever with end load P, find deflection at midspan.
    """
    import sympy as sp

    # Symbolic variables
    P, Q, L, E, I, x = sp.symbols('P Q L E I x', positive=True, real=True)

    # Moment expression (with dummy load Q at L/2)
    # For 0 ≤ x ≤ L/2: M = P(L-x) + Q(L/2-x)
    # For L/2 < x ≤ L: M = P(L-x)

    # Strain energy (integrating M²/2EI)
    M1 = P*(L-x) + Q*(L/2 - x)  # 0 to L/2
    M2 = P*(L-x)                 # L/2 to L

    U1 = sp.integrate(M1**2 / (2*E*I), (x, 0, L/2))
    U2 = sp.integrate(M2**2 / (2*E*I), (x, L/2, L))
    U_total = U1 + U2

    # Derivative with respect to Q, then set Q=0
    delta_mid = sp.diff(U_total, Q)
    delta_mid = delta_mid.subs(Q, 0)
    delta_mid = sp.simplify(delta_mid)

    print("Castigliano's Theorem - Dummy Load Method")
    print("-" * 50)
    print(f"Midspan deflection of cantilever with end load:")
    print(f"  δ_mid = {delta_mid}")

    # Numerical example
    result = delta_mid.subs([(P, 10000), (L, 2000), (E, 200e3), (I, 1e6)])
    print(f"\nFor P=10kN, L=2000mm, E=200GPa, I=10⁶mm⁴:")
    print(f"  δ_mid = {float(result):.4f} mm")

    # Compare with tip deflection
    delta_tip = P * L**3 / (3 * E * I)
    delta_tip_val = float(delta_tip.subs([(P, 10000), (L, 2000), (E, 200e3), (I, 1e6)]))
    print(f"  δ_tip = {delta_tip_val:.4f} mm")
    print(f"  Ratio: δ_mid/δ_tip = {float(result)/delta_tip_val:.4f}")

try:
    castigliano_midspan_deflection()
except ImportError:
    print("Note: SymPy required for symbolic calculations")
    print("Install with: pip install sympy")

Numerical Integration for Complex Cases

def strain_energy_numerical(M_func, L_mm, E_GPa, I_mm4, n_points=100):
    """
    Calculate strain energy by numerical integration.

    U = ∫(M²/2EI)dx

    Parameters:
    -----------
    M_func : callable - Function returning moment at position x
    """
    E = E_GPa * 1e3
    x = np.linspace(0, L_mm, n_points)
    M = np.array([M_func(xi) for xi in x])

    # Integrand: M²/(2EI)
    integrand = M**2 / (2 * E * I_mm4)

    # Numerical integration using trapezoidal rule
    U = np.trapz(integrand, x)

    return {
        'U_mJ': U,
        'x': x,
        'M': M,
        'integrand': integrand
    }

# Example: Simply supported beam with UDL
def moment_ss_udl(x, w, L):
    """Moment for simply supported beam with UDL."""
    return w * x * (L - x) / 2

L = 4000  # mm
w = 10    # N/mm
E = 200   # GPa
I = 5e6   # mm⁴

M_func = lambda x: moment_ss_udl(x, w, L)
result = strain_energy_numerical(M_func, L, E, I)

print(f"\nSimply Supported Beam with UDL (w={w} N/mm, L={L}mm):")
print(f"  Strain energy (numerical): {result['U_mJ']:.2f} mJ")

# Analytical formula: U = w²L⁵/(240EI)
U_analytical = w**2 * L**5 / (240 * E * 1e3 * I)
print(f"  Strain energy (analytical): {U_analytical:.2f} mJ")

Statically Indeterminate Structures

Energy methods are powerful for solving indeterminate structures.

def propped_cantilever_reaction():
    """
    Solve propped cantilever using Castigliano's theorem.

    Cantilever with prop at free end under UDL.
    Find reaction at prop (R).

                w (UDL)
    ↓↓↓↓↓↓↓↓↓↓↓↓
    ================
    |              △
    Fixed          Prop (R)
    |-----L-------|
    """
    print("Propped Cantilever - Castigliano's Method")
    print("-" * 50)

    # At prop, deflection = 0
    # ∂U/∂R = 0 (minimum potential energy)

    # Moment: M(x) = R(L-x) - w(L-x)²/2
    # where R is the prop reaction (unknown)

    # For compatibility: deflection at prop = 0
    # This leads to: R = 3wL/8

    import sympy as sp
    w, L, x, R, E, I = sp.symbols('w L x R E I', positive=True)

    M = R*(L-x) - w*(L-x)**2/2

    # Strain energy
    U = sp.integrate(M**2 / (2*E*I), (x, 0, L))

    # For zero deflection at prop: ∂U/∂R = 0
    dU_dR = sp.diff(U, R)
    R_solution = sp.solve(dU_dR, R)[0]

    print(f"Reaction at prop: R = {R_solution}")
    print(f"                    = {sp.simplify(R_solution)}")

    # Numerical example
    w_val = 10  # N/mm
    L_val = 3000  # mm

    R_val = float(R_solution.subs([(w, w_val), (L, L_val)]))
    print(f"\nFor w={w_val} N/mm, L={L_val} mm:")
    print(f"  Prop reaction R = {R_val:.2f} N = {R_val/1000:.3f} kN")
    print(f"  Total load wL = {w_val * L_val / 1000:.2f} kN")
    print(f"  Fixed end reaction = {w_val * L_val - R_val:.2f} N")

try:
    propped_cantilever_reaction()
except ImportError:
    print("SymPy required for this example")

Maxwell-Betti Reciprocal Theorem

The deflection at point A due to load at B equals deflection at B due to same load at A:

$$\delta_{AB} = \delta_{BA}$$

def demonstrate_reciprocity(L_mm, E_GPa, I_mm4):
    """
    Demonstrate Maxwell-Betti reciprocal theorem.
    """
    E = E_GPa * 1e3
    P = 10000  # N

    print("Maxwell-Betti Reciprocal Theorem")
    print("-" * 50)
    print(f"Simply supported beam, L={L_mm}mm")
    print(f"Load P={P/1000}kN")

    # Point A at L/3, Point B at 2L/3
    a = L_mm / 3
    b = 2 * L_mm / 3

    # Deflection at A due to load at B
    # Using superposition formula for point load
    def deflection_at_x_from_load_at_a(x, a_load, L, P, E, I):
        """Deflection at x from point load P at a_load."""
        b_load = L - a_load
        if x <= a_load:
            return P * b_load * x * (L**2 - b_load**2 - x**2) / (6 * L * E * I)
        else:
            return P * a_load * (L - x) * (2*L*x - a_load**2 - x**2) / (6 * L * E * I)

    # δ_AB: deflection at A (L/3) due to load at B (2L/3)
    delta_AB = deflection_at_x_from_load_at_a(a, b, L_mm, P, E, I_mm4)

    # δ_BA: deflection at B (2L/3) due to load at A (L/3)
    delta_BA = deflection_at_x_from_load_at_a(b, a, L_mm, P, E, I_mm4)

    print(f"\nLoad at B (x={b:.0f}mm), measure at A (x={a:.0f}mm):")
    print(f"  δ_AB = {delta_AB:.6f} mm")

    print(f"\nLoad at A (x={a:.0f}mm), measure at B (x={b:.0f}mm):")
    print(f"  δ_BA = {delta_BA:.6f} mm")

    print(f"\nReciprocity verified: δ_AB ≈ δ_BA")
    print(f"  Difference: {abs(delta_AB - delta_BA):.10f} mm")

demonstrate_reciprocity(L_mm=3000, E_GPa=200, I_mm4=5e6)

Summary of Energy Methods

Key Takeaways

• Strain energy is stored elastic energy: $U = \int \sigma \varepsilon \, dV / 2$
• Castigliano's theorem: $\delta = \partial U / \partial P$
• Dummy load method: Add fictitious load where deflection needed
• Indeterminate structures: Use $\partial U / \partial R = 0$ for redundant reactions
• Maxwell-Betti: $\delta_{AB} = \delta_{BA}$ (reciprocity)

Next lesson: We'll build a complete beam analysis tool in Python.