[Python] How to implement Paris’ law crack growth simulation in Python?

I am trying to model fatigue crack growth in Python using Paris’ law:

da/dN = C (Delta K)^m

where:

a= crack length

N= load cycles

(stress intensity factor range)(deltaK)=Y.(delta sigma).(pi.a)^1/2

c,m= material constants

For a simple case, I set:c=1e-12

m=3

(delta sigma=100,mpa)

(geometry factor)Y=1.12

Initial crack length(0.001)m

Here is my first attempt in Python:

import numpy as np
import matplotlib.pyplot as plt

# Parameters
C, m = 1e-12, 3
sigma = 100e6 # Pa
Y = 1.12
a0 = 0.001 # m
Nmax = 1_000_000
dN = 1000

# Crack growth function
def crack_growth(a, dN):
dK = Y * sigma * np.sqrt(np.pi * a)
da = C * (dK**m) * dN
return a + da

# Iterative loop
a = a0
a_history = [a]
N_history = [0]

for N in range(0, Nmax, dN):
a = crack_growth(a, dN)
a_history.append(a)
N_history.append(N + dN)

plt.plot(N_history, a_history)
plt.xlabel("Load cycles (N)")
plt.ylabel("Crack length (m)")
plt.show()

1. 
   Is this iterative approach (Euler forward) numerically stable for crack growth problems, or should I use an ODE solver like scipy.integrate.solve_ivp?
   
   
2. 
   How can I better stop the simulation when the crack reaches a critical length (e.g., )?
   
   
3. 
   For variable amplitude loading (changing ), what would be the recommended way to implement it efficiently?
   

Continue reading...