Chapter 4: Fatigue and Fracture Mechanics

# Requirements:
# - Python 3.9+
# - matplotlib>=3.7.0
# - numpy>=1.24.0, <2.0.0

import numpy as np
import matplotlib.pyplot as plt

def generate_sn_curve(material='steel'):
    """Generate S-N (Wöhler) curve"""
    materials = {
        'steel': {'S_f': 200, 'b': -0.12, 'N_e': 1e6},
        'aluminum': {'S_f': 100, 'b': -0.10, 'N_e': 5e8}
    }
    props = materials[material]
    
    # Basquin equation: S = S_f * N^b
    N = np.logspace(3, 9, 100)
    S = props['S_f'] * (N / 2e3)**props['b']
    
    # Endurance limit for steel
    if material == 'steel':
        S[N > props['N_e']] = S[N > props['N_e']][0]
    
    return N, S

fig, ax = plt.subplots(figsize=(10, 6))
for material in ['steel', 'aluminum']:
    N, S = generate_sn_curve(material)
    ax.loglog(N, S, linewidth=2, label=material.capitalize())

ax.set_xlabel('Cycles to Failure (N)', fontsize=12)
ax.set_ylabel('Stress Amplitude (MPa)', fontsize=12)
ax.set_title('S-N Curves for Different Materials', fontsize=14, fontweight='bold')
ax.legend()
ax.grid(True, which='both', alpha=0.3)
plt.show()

def goodman_diagram(S_u=500, S_f=200, S_y=300):
    """Generate Goodman diagram"""
    sigma_m = np.linspace(0, S_u, 100)
    
    # Goodman line
    sigma_a_goodman = S_f * (1 - sigma_m / S_u)
    
    # Gerber parabola
    sigma_a_gerber = S_f * (1 - (sigma_m / S_u)**2)
    
    # Soderberg line
    sigma_a_soderberg = S_f * (1 - sigma_m / S_y)
    
    plt.figure(figsize=(10, 7))
    plt.plot(sigma_m, sigma_a_goodman, 'b-', linewidth=2, label='Goodman')
    plt.plot(sigma_m, sigma_a_gerber, 'r-', linewidth=2, label='Gerber')
    plt.plot(sigma_m, sigma_a_soderberg, 'g-', linewidth=2, label='Soderberg')
    plt.fill_between(sigma_m, 0, sigma_a_goodman, alpha=0.2)
    plt.xlabel('Mean Stress σ_m (MPa)', fontsize=12)
    plt.ylabel('Alternating Stress σ_a (MPa)', fontsize=12)
    plt.title('Goodman Diagram', fontsize=14, fontweight='bold')
    plt.legend()
    plt.grid(True, alpha=0.3)
    plt.show()

goodman_diagram()

def stress_intensity_factor(sigma, crack_length, geometry='center_crack'):
    """Calculate stress intensity factor K_I"""
    a = crack_length
    
    # Geometry factors
    if geometry == 'center_crack':
        Y = 1.0  # For infinite plate
    elif geometry == 'edge_crack':
        Y = 1.12
    elif geometry == 'semi_elliptical':
        Y = 0.73
    
    K_I = Y * sigma * np.sqrt(np.pi * a)
    return K_I

# Example: crack growth analysis
crack_lengths = np.linspace(1, 20, 100)  # mm
sigma = 100  # MPa
K_I = [stress_intensity_factor(sigma, a/1000, 'edge_crack') for a in crack_lengths]

plt.figure(figsize=(10, 6))
plt.plot(crack_lengths, K_I, 'b-', linewidth=2)
plt.axhline(30, color='r', linestyle='--', label='K_Ic = 30 MPa√m')
plt.xlabel('Crack Length (mm)', fontsize=12)
plt.ylabel('Stress Intensity Factor K_I (MPa√m)', fontsize=12)
plt.title('Stress Intensity vs Crack Length', fontsize=14, fontweight='bold')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

def paris_law_integration(a_0, sigma, Y, K_Ic, C=1e-11, m=3):
    """Integrate Paris law to predict crack growth"""
    a = a_0
    N = 0
    a_history = [a]
    N_history = [N]
    
    while True:
        # Calculate K
        K_max = Y * sigma * np.sqrt(np.pi * a)
        K_min = 0  # R = 0
        Delta_K = K_max - K_min
        
        # Check failure
        if K_max >= K_Ic:
            break
        
        # Paris law: da/dN
        da_dN = C * (Delta_K * 1e6)**m  # Convert to Pa√m
        
        # Update crack length
        dN = 100  # Cycle increment
        a = a + da_dN * dN
        N = N + dN
        
        a_history.append(a)
        N_history.append(N)
        
        if N > 1e7:  # Maximum cycles
            break
    
    return np.array(N_history), np.array(a_history)

# Simulation
N, a = paris_law_integration(a_0=0.001, sigma=100, Y=1.12, K_Ic=30e6)

plt.figure(figsize=(10, 6))
plt.plot(N, a*1000, 'b-', linewidth=2)
plt.xlabel('Cycles (N)', fontsize=12)
plt.ylabel('Crack Length (mm)', fontsize=12)
plt.title('Fatigue Crack Growth Prediction', fontsize=14, fontweight='bold')
plt.grid(True, alpha=0.3)
plt.show()

print(f"Cycles to failure: {N[-1]:.0f}")

def coffin_manson(N_f, epsilon_f=0.5, c=-0.6):
    """Coffin-Manson equation for LCF"""
    # Δε_p/2 = ε_f' * (2*N_f)^c
    delta_epsilon_p = epsilon_f * (2 * N_f)**c
    return delta_epsilon_p

N_values = np.logspace(1, 5, 100)
epsilon_p = coffin_manson(N_values)

plt.figure(figsize=(10, 6))
plt.loglog(N_values, epsilon_p, 'b-', linewidth=2)
plt.xlabel('Cycles to Failure (N_f)', fontsize=12)
plt.ylabel('Plastic Strain Amplitude', fontsize=12)
plt.title('Coffin-Manson Relationship (LCF)', fontsize=14, fontweight='bold')
plt.grid(True, which='both', alpha=0.3)
plt.show()

def miners_rule(stress_levels, cycles_applied):
    """Calculate cumulative fatigue damage"""
    # S-N curve parameters
    S_f, b = 200, -0.12
    
    damage = 0
    for S, n in zip(stress_levels, cycles_applied):
        # Calculate N_f from S-N curve
        N_f = (S / S_f)**(1/b) * 2e3
        
        # Add damage fraction
        damage += n / N_f
    
    return damage

# Example: Variable amplitude loading
stress_levels = [150, 200, 250]  # MPa
cycles_applied = [50000, 20000, 5000]

damage = miners_rule(stress_levels, cycles_applied)
remaining_life_fraction = 1 - damage

print(f"Cumulative Damage: {damage:.3f}")
print(f"Remaining Life: {remaining_life_fraction*100:.1f}%")
if damage >= 1:
    print("Failure predicted!")

def calculate_KIc(P_Q, B, W, a):
    """Calculate K_Ic from compact tension test"""
    # P_Q: Load at 5% secant
    # B: Thickness, W: Width, a: Crack length
    
    alpha = a / W
    
    # Geometry function for CT specimen
    f_alpha = (2 + alpha) / (1 - alpha)**1.5 * \
              (0.886 + 4.64*alpha - 13.32*alpha**2 + 14.72*alpha**3 - 5.6*alpha**4)
    
    K_Q = (P_Q / (B * np.sqrt(W))) * f_alpha
    
    # Validity checks (ASTM E399)
    valid = True
    if a < 0.45*W or a > 0.55*W:
        valid = False
    
    return K_Q, valid

# Example
P_Q = 15000  # N
B, W, a = 25e-3, 50e-3, 25e-3  # m

K_Ic, valid = calculate_KIc(P_Q, B, W, a)
print(f"K_Ic = {K_Ic/1e6:.1f} MPa√m")
print(f"Test valid: {valid}")