Spring Design with Engineering Stress Verification

# Standard music wire diameters (mm) — partial table
MUSIC_WIRE_DIAMETERS = [
    0.10, 0.12, 0.14, 0.16, 0.18, 0.20, 0.22, 0.25,
    0.28, 0.30, 0.35, 0.40, 0.45, 0.50, 0.55, 0.60,
    0.65, 0.70, 0.80, 0.90, 1.00, 1.10, 1.20, 1.40,
    1.60, 1.80, 2.00, 2.20, 2.50, 2.80, 3.00, 3.50,
    4.00, 4.50, 5.00, 5.50, 6.00, 6.50, 7.00, 8.00,
    9.00, 10.00, 11.00, 12.00, 13.00, 14.00, 15.00, 16.00,
]

# Approximate tensile strength of music wire (Sut in MPa)
# Using the empirical formula: Sut = A / d^m
# For music wire (ASTM A228): A = 2211 MPa·mm^m, m = 0.145
def music_wire_tensile_strength(d_mm: float) -> float:
    """Approximate ultimate tensile strength of music wire.

    Based on curve fit to ASTM A228 data.
    Args:
        d_mm: Wire diameter in mm
    Returns:
        Ultimate tensile strength in MPa
    """
    A = 2211.0  # MPa·mm^m
    m = 0.145
    return A / (d_mm ** m)

def select_wire_diameter(d_required: float) -> float:
    """Select the nearest standard wire diameter >= required."""
    for d_std in MUSIC_WIRE_DIAMETERS:
        if d_std >= d_required:
            return d_std
    raise ValueError(
        f"Required diameter {d_required} mm exceeds "
        f"maximum standard size {MUSIC_WIRE_DIAMETERS[-1]} mm"
    )

import numpy as np

def wahl_factor(C: float) -> float:
    """Compute the Wahl correction factor.

    Args:
        C: Spring index (D/d), must be >= 3

    Returns:
        Wahl factor K_w
    """
    if C < 3:
        raise ValueError(f"Spring index C={C} is too low (min 3)")
    return (4*C - 1) / (4*C - 4) + 0.615 / C

def max_shear_stress(F: float, D: float, d: float) -> float:
    """Compute maximum shear stress with Wahl correction.

    Args:
        F: Applied force (N)
        D: Mean coil diameter (mm)
        d: Wire diameter (mm)

    Returns:
        Maximum shear stress (MPa)
    """
    C = D / d
    Kw = wahl_factor(C)
    tau = Kw * 8 * F * D / (np.pi * d**3)
    return tau

# Example: 50 N force, 20mm mean coil diameter, 2mm wire
tau = max_shear_stress(F=50, D=20, d=2.0)
print(f"Wahl factor: {wahl_factor(20/2):.3f}")
print(f"Max shear stress: {tau:.1f} MPa")

def goodman_fatigue_check(
    F_min: float, F_max: float,
    D: float, d: float,
    S_ut: float,
    shot_peened: bool = False,
    n_f_required: float = 1.5
) -> dict:
    """Perform Goodman fatigue analysis on a helical spring.

    Args:
        F_min: Minimum cyclic force (N)
        F_max: Maximum cyclic force (N)
        D: Mean coil diameter (mm)
        d: Wire diameter (mm)
        S_ut: Ultimate tensile strength of wire (MPa)
        shot_peened: Whether the spring is shot-peened
        n_f_required: Required fatigue safety factor

    Returns:
        Dictionary with stress values and safety factor
    """
    C = D / d
    Kw = wahl_factor(C)

    # Alternating and mean forces
    F_a = (F_max - F_min) / 2
    F_m = (F_max + F_min) / 2

    # Corrected alternating and mean shear stresses
    tau_a = Kw * 8 * F_a * D / (np.pi * d**3)
    tau_m = Kw * 8 * F_m * D / (np.pi * d**3)

    # Material properties
    S_su = 0.67 * S_ut    # Ultimate shear strength
    if shot_peened:
        S_se = 0.40 * S_ut  # Endurance limit (peened)
    else:
        S_se = 0.30 * S_ut  # Endurance limit (unpeened)

    # Goodman safety factor
    # tau_a/S_se + tau_m/S_su = 1/n_f
    # Solving for n_f:
    goodman_sum = tau_a / S_se + tau_m / S_su
    if goodman_sum > 0:
        n_f = 1.0 / goodman_sum
    else:
        n_f = float('inf')

    return {
        'F_a': F_a,
        'F_m': F_m,
        'tau_a_MPa': tau_a,
        'tau_m_MPa': tau_m,
        'S_se_MPa': S_se,
        'S_su_MPa': S_su,
        'Kw': Kw,
        'C': C,
        'n_f': n_f,
        'infinite_life': n_f >= n_f_required,
    }

# Example: spring cycling between 20N and 80N
result = goodman_fatigue_check(
    F_min=20, F_max=80,
    D=18, d=2.0,
    S_ut=music_wire_tensile_strength(2.0),
    shot_peened=False,
    n_f_required=1.5
)

for key, val in result.items():
    if isinstance(val, float):
        print(f"  {key}: {val:.3f}")
    else:
        print(f"  {key}: {val}")

import matplotlib.pyplot as plt
import numpy as np

def plot_goodman_diagram(result: dict, title: str = "Goodman Diagram"):
    """Plot the Goodman diagram showing the operating point
    relative to the failure boundary.

    Args:
        result: Output from goodman_fatigue_check()
    """
    S_se = result['S_se_MPa']
    S_su = result['S_su_MPa']
    tau_a = result['tau_a_MPa']
    tau_m = result['tau_m_MPa']

    fig, ax = plt.subplots(1, 1, figsize=(8, 6))

    # Goodman line (failure boundary)
    tau_m_line = np.linspace(0, S_su, 100)
    tau_a_line = S_se * (1 - tau_m_line / S_su)
    ax.plot(tau_m_line, tau_a_line, 'r-', linewidth=2,
            label='Goodman line (failure)')

    # Yield line: tau_a + tau_m = S_sy
    S_sy = 0.45 * (S_se / 0.30)  # Back-calculate S_ut
    tau_a_yield = np.linspace(0, S_sy, 100)
    tau_m_yield = S_sy - tau_a_yield
    ax.plot(tau_m_yield, tau_a_yield, 'b--', linewidth=1.5,
            label='Yield line')

    # Safe region (shaded)
    ax.fill_between(tau_m_line, 0, tau_a_line,
                     alpha=0.1, color='green', label='Safe region')

    # Operating point
    ax.plot(tau_m, tau_a, 'ko', markersize=10, zorder=5)
    ax.annotate(
        f'Operating point\n'
        f'n_f = {result["n_f"]:.2f}',
        xy=(tau_m, tau_a),
        xytext=(tau_m + S_su*0.05, tau_a + S_se*0.1),
        fontsize=10,
        arrowprops=dict(arrowstyle='->', color='black'),
        bbox=dict(boxstyle='round,pad=0.3', facecolor='yellow',
                  alpha=0.8),
    )

    # Load line from origin through operating point
    if tau_m > 0:
        slope = tau_a / tau_m
        tau_m_load = np.linspace(0, S_su, 100)
        tau_a_load = slope * tau_m_load
        mask = tau_a_load <= S_se * 1.1
        ax.plot(tau_m_load[mask], tau_a_load[mask], 'g:',
                linewidth=1, label='Load line')

    ax.set_xlabel('Mean Shear Stress τ_m (MPa)', fontsize=12)
    ax.set_ylabel('Alternating Shear Stress τ_a (MPa)', fontsize=12)
    ax.set_title(title, fontsize=14)
    ax.legend(loc='upper right')
    ax.set_xlim(0, S_su * 1.1)
    ax.set_ylim(0, S_se * 1.1)
    ax.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.savefig("goodman_diagram.png", dpi=150)
    plt.show()

# Plot for our example
# plot_goodman_diagram(result, "Spring Fatigue — Goodman Diagram")

import numpy as np

# === DESIGN REQUIREMENTS ===
F_max = 80.0     # N, maximum working load
F_min = 20.0     # N, minimum load (installed)
delta_w = 15.0   # mm, working deflection (F_max - F_min)
OD_max = 22.0    # mm, maximum outer diameter
material = "music_wire"
G = 81_700.0     # MPa, shear modulus (music wire)
n_s = 1.3        # Static safety factor
n_f = 1.5        # Fatigue safety factor

# === STEP 1: Spring Rate ===
k = (F_max - F_min) / delta_w  # N/mm
print(f"Required spring rate: k = {k:.2f} N/mm")

# === STEP 2: Wire Diameter and Coil Geometry ===
# Start with a trial spring index C = 8 (good starting point)
C_trial = 8.0

# Iterate: select wire diameter, compute coil geometry
# The design equation for wire diameter from stress limit:
# d >= (8 * K_w * F_max * C / (pi * tau_allow))^(1/3)

# Trial wire diameters
for d_trial in [1.6, 1.8, 2.0, 2.2, 2.5]:
    D = C_trial * d_trial
    OD = D + d_trial

    if OD > OD_max:
        continue

    # Number of active coils
    N_a = G * d_trial**4 / (8 * D**3 * k)

    if N_a < 3 or N_a > 20:
        continue

    # Free length (squared-and-ground ends)
    # Solid height
    L_s = (N_a + 2) * d_trial
    # Deflection at F_max
    delta_max = F_max / k
    # Clash allowance (15% of working deflection)
    delta_clash = 0.15 * delta_w
    L_f = L_s + delta_max + delta_clash

    # Check stress
    Kw = wahl_factor(C_trial)
    tau_max = Kw * 8 * F_max * D / (np.pi * d_trial**3)
    S_ut = music_wire_tensile_strength(d_trial)
    S_sy = 0.45 * S_ut
    tau_allow = S_sy / n_s

    # Check fatigue
    fatigue = goodman_fatigue_check(
        F_min, F_max, D, d_trial, S_ut,
        shot_peened=False, n_f_required=n_f
    )

    print(f"\n--- Trial: d = {d_trial} mm ---")
    print(f"  D = {D:.1f} mm, OD = {OD:.1f} mm")
    print(f"  N_a = {N_a:.1f} coils")
    print(f"  L_f = {L_f:.1f} mm, L_s = {L_s:.1f} mm")
    print(f"  tau_max = {tau_max:.1f} MPa")
    print(f"  tau_allow = {tau_allow:.1f} MPa")
    print(f"  Static OK: {tau_max < tau_allow}")
    print(f"  S_ut = {S_ut:.0f} MPa")
    print(f"  Fatigue n_f = {fatigue['n_f']:.2f}")
    print(f"  Infinite life: {fatigue['infinite_life']}")

# === STEP 3: Select Best Design ===
# Choose d = 2.0 mm based on results above
d = 2.0                          # mm (standard music wire)
D = C_trial * d                  # 16.0 mm mean coil diameter
OD = D + d                       # 18.0 mm outer diameter
ID = D - d                       # 14.0 mm inner diameter
N_a = G * d**4 / (8 * D**3 * k)  # Active coils
N_a = round(N_a)                 # Round to nearest integer

# Recalculate actual spring rate with rounded N_a
k_actual = G * d**4 / (8 * D**3 * N_a)
print(f"\n=== FINAL DESIGN ===")
print(f"Wire diameter:     d = {d} mm")
print(f"Mean coil diameter: D = {D} mm")
print(f"Outer diameter:    OD = {OD} mm")
print(f"Inner diameter:    ID = {ID} mm")
print(f"Active coils:      N_a = {N_a}")
print(f"Spring rate:       k = {k_actual:.3f} N/mm")

# Physical dimensions
L_s = (N_a + 2) * d
delta_max = F_max / k_actual
delta_clash = 0.15 * (F_max - F_min) / k_actual
L_f = L_s + delta_max + delta_clash
pitch = (L_f - 2 * d) / N_a

print(f"Solid height:      L_s = {L_s:.1f} mm")
print(f"Free length:       L_f = {L_f:.1f} mm")
print(f"Pitch:             p = {pitch:.2f} mm")
print(f"Total coils:       N_t = {N_a + 2}")

# === STEP 4: Final Verification ===
Kw = wahl_factor(D / d)
tau_max = Kw * 8 * F_max * D / (np.pi * d**3)
S_ut = music_wire_tensile_strength(d)
S_sy = 0.45 * S_ut
n_static = S_sy / tau_max

fatigue_final = goodman_fatigue_check(
    F_min, F_max, D, d, S_ut,
    shot_peened=False, n_f_required=n_f
)

print(f"\n=== STRESS VERIFICATION ===")
print(f"Wahl factor:       K_w = {Kw:.3f}")
print(f"Max shear stress:  tau = {tau_max:.1f} MPa")
print(f"Shear yield:       S_sy = {S_sy:.0f} MPa")
print(f"Static safety:     n_s = {n_static:.2f}")
print(f"Fatigue safety:    n_f = {fatigue_final['n_f']:.2f}")
print(f"Infinite life:     {fatigue_final['infinite_life']}")

# === STEP 5: Generate 3D Geometry ===
spring_3d = compression_spring(
    d=d,
    D=D,
    N_a=N_a,
    L_f=L_f,
    end_type="squared_ground"
)

# Export production files
cq.exporters.export(spring_3d, "spring_verified.step")
cq.exporters.export(spring_3d, "spring_verified.stl")

print("\nExported: spring_verified.step, spring_verified.stl")

def extension_spring_design(
    F_max: float,       # Maximum load (N)
    F_i: float,         # Initial tension (N)
    delta_max: float,   # Maximum deflection (mm)
    D: float,           # Mean coil diameter (mm)
    d: float,           # Wire diameter (mm)
    G: float = 81_700,  # Shear modulus (MPa)
) -> dict:
    """Design an extension spring with initial tension.

    Returns dictionary of design parameters.
    """
    # Spring rate (not counting initial tension)
    k = (F_max - F_i) / delta_max

    # Active coils
    N_a = G * d**4 / (8 * D**3 * k)

    # Body length (coils touching)
    L_body = (N_a + 1) * d  # +1 for end coil

    # Hook stress concentration factor (typical)
    # Hooks have higher stress than the body
    K_hook = (4 * (D/d)**2 - (D/d) - 1) / (4 * (D/d) * ((D/d) - 1))

    # Maximum body stress
    C = D / d
    Kw = wahl_factor(C)
    tau_body = Kw * 8 * F_max * D / (np.pi * d**3)

    # Hook bending stress (the critical point)
    sigma_hook = K_hook * 16 * F_max * D / (np.pi * d**3)

    return {
        'k_N_per_mm': k,
        'N_a': N_a,
        'L_body_mm': L_body,
        'tau_body_MPa': tau_body,
        'sigma_hook_MPa': sigma_hook,
        'K_hook': K_hook,
    }