# Simulating Crack Propagation

In this notebook, we will examine some of the equations and numerical techniques for simulating **dynamic crack propagation**. The context will be **earthquake source mechanics** but much of the approach can be generalised to other problems.

## Basic Equations

From the notes on Earthquakes, we derived a frequency space expression for the stress change on a slipping crack under the quasidynamic assumption.

\begin{equation}
\hat{\tau} = \hat{\tau}_0 - \frac{i\mu\omega}{2c}\hat{\delta}-\frac{\mu|k|}{2}\hat{\delta}.
\end{equation}

The corresponding equation in real space is

\begin{equation}
\tau=\tau_0-\frac{\mu}{2c}V - \frac{\mu}{2\pi}\int\limits_{-\infty}^\infty \frac{d\delta/dx'}{x-x'}dx'.
\end{equation}

where $V=d\delta/dt$ is called the slip velocity.

## Part 1: Computing Slip-induced Stress Changes

To solve dynamic crack problems, we need a method to compute the integral term above. Recall, this is the stress change on the crack resulting from slip, $\delta$, after all seismic waves have dissipated. It is only realised at very long time, and is sometimes referred to as the 'quasistatic term'.

One way to solve this term, for arbitrary slip function, is by directly evaluating the integral. 

Let's see how that works. ***First, execute the cell below to create a slip profile.***


```python
%matplotlib inline
import numpy as np
from matplotlib import pyplot as plt
x = np.linspace(0,1.e3,128)      # spatial coordinate (1 km)
c,a,umax = 500,250,1             # center and aperture of crack, maximum slip
delta = np.array([np.sqrt(a**2 - (xi-c)**2)/a*umax if abs(xi-c)<a else 0 for xi in x ])   # elliptical slip profile
f,ax = plt.subplots(1,1)
ax.plot(x,delta,'k-')
ax.set_xlabel('x [m]')
ax.set_ylabel('slip [m]')
plt.show()
```

Now, compute the integral **at each $x$** using the Rectangular rule.


```python
xm = 0.5*(x[:-1]+x[1:])              # cell mid-points
ddelta = np.diff(delta)/np.diff(x)   # slip gradient at midpoints
dts = 0.*x                           # empty vector to compute quasistatic stress change at midpoints
dx = x[1]-x[0]
mu = 30.e9                           # typical shear modulus of rock, 30 GPa
for i,xi in enumerate(x):           # loop over each value xm and compute dts
    dts[i] = -mu/2./np.pi*np.sum(ddelta/(xi-xm))*dx     # Rectangular rule

f,ax = plt.subplots(1,1)
ax.plot(x,delta,'k-')
ax.set_xlabel('x [m]')
ax.set_ylabel('slip [m]')
ax1=ax.twinx()
ax1.plot(x,dts/1.e6,'b-')
ax1.set_ylabel('stress change [MPa]')
plt.show()
```

An alternative way to compute the quasistatic stress term is to recognise that the integral is formulated as a convolution. Therefore, we can ***multiply*** the Fourier transforms of the terms $1/x$ and $\delta$ with appropriate frequency space terms corresponding to the spatial derivative $d/dx$. 

\begin{equation}
\int\limits_{-\infty}^\infty \frac{d\delta/dx'}{x-x'}dx' \Longleftrightarrow \pi\,|k|\,\hat{\delta}
\end{equation}

The algorithm then becomes

1. Compute (absolute values of) the spatial wavenumbers, $|k|$. This needs only be done once.
2. Take the Discrete Fourier Transform (DFT) of the slip vector, $\delta$. This needs to be done each timestep, because $\delta$ evolves with time.
3. Multiply $|k|$ and $\hat{\delta}$.
4. Take the Inverse Discrete Fourier Transform (IDFT) of the product.

We will use the Fast Fourier Transform algorithm to compute the DFTs and IDFTs.


```python
k = np.fft.fftfreq(len(x), x[1]-x[0])      # get the wave numbers - do this once
deltahat = np.fft.fft(delta)               # FFT of slip vector - do this each time step
dts_FFT = -mu/2.*np.fft.ifft(abs(k)*deltahat)*2*np.pi   # inverse transform of product

f,ax = plt.subplots(1,1)
ax.plot(x,delta,'k-')
ax.set_xlabel('x [m]')
ax.set_ylabel('slip [m]')
ax1=ax.twinx()
ax1.plot(x,dts/1.e6,'b--',label='direct')
ax1.plot(x,dts_FFT/1.e6,'g-', label = 'FFT')
ax1.set_ylabel('stress change [MPa]')
ax1.legend()
plt.show()
```

Not the same, but similar. The differences arise due to a periodic boundary condition when implementing the FFT method. This assumes that slip is not localised to the crack, but rather is an infinitely tesselating pattern (an infinite number of slipping cracks, all affecting each other). One way to counteract this effect is to ***zero pad*** the slip vector.

***Change the crack width from 250m to 25m and note how agreement between the FFT and direct evaluation improve***.

***Define some other slip functions and recomputing the quasistatic stress change.***

- ** A sinusoid.**
- ** A square function.**
- ** A Gaussian. **

***What causes a stress singularity?***

What motivates an FFT method - why should we bother? The main reason is some quite significant gains in **computational efficiency**.

The code below simulates many repetitions of the direct integration and FFT methods for different vector lengths, $N$.

***Execute the cell below to analyse scaling performance for the two methods***


```python
from time import time

def dts_direct(x,u,reps):
    ''' return time to evaluate direct integral REPS times
    '''
    # start the timer
    tstart = time()
    
    # outside loop, calculations done once at beginning of simulation
    xm = 0.5*(x[:-1]+x[1:])     
    mu = 30.e9               
    dx = (x[1]-x[0])*-mu/2./np.pi                     
    dts = 0.*x                           
    for i in range(reps):
        # inside loop, calculations that change each time
        ddelta = np.diff(delta)/np.diff(x)   
        for i,xi in enumerate(x):           
            dts[i] = np.sum(ddelta/(xi-xm))*dx     
            
    # end the timer
    tend = time()
    
    return (tend-tstart)/reps
    
def dts_FFT(x,u,reps):
    ''' return time to evaluate FFT integral REPS times
    '''
    # start the timer
    tstart = time()
    
    # outside loop, calculations done once at beginning of simulation
    mu = 30.e9                           # typical shear modulus of rock, 30 GPa
    absk = -mu/2.*abs(np.fft.fftfreq(len(x), x[1]-x[0]))*2*np.pi      # get the wave numbers
    for i in range(reps):
        # inside loop, calculations that change each time        
        deltahat = np.fft.fft(delta)               # FFT of slip vector
        dts_FFT = np.fft.ifft(absk*deltahat)
        
    # end the timer
    tend = time()
    
    return (tend-tstart)/reps
    

# loop over different size vectors and record execution times
n = 2**np.array(range(6,12))
time_direct = []
time_FFT = []
for ni in n:
    x = np.linspace(0,1.,ni)      # spatial coordinate (1 km)
    delta = np.random.rand(ni)    # slip function (not important what)
    time_FFT.append(dts_FFT(x,delta,1000))     # have to do heaps of REPS to register on system timing
    time_direct.append(dts_direct(x,delta,100))  # not for the direct integration though, so slow...

# plot scalings
f,ax = plt.subplots(1,1)
ax.plot(n, time_direct, 'ks',label='direct integration')
ax.plot(n, time_FFT, 'ko', mfc = 'w',label='FFT')
ax.legend()
ax.set_xlabel('vector length')
ax.set_ylabel('execution time [s]')
ax.set_xscale('log')
ax.set_yscale('log')
plt.show()    
```

*Estimate the gradients for curves fitted to the points above. Hence propose scaling laws for the two methods.*

*Do some Googling and see how your answer compares.*

## Part 2: Solving for Crack Propagation

To solve the dynamic crack problem, we first need to formulate some evolution equations for the variables of interest. These will be **friction**, $f(x,t)$, and **slip**, $\delta(x,t)$.

Recall, we introduced a **linear slip-weakening** friction formulation

\begin{equation}
f=\begin{cases} 
      f_s - (f_s-f_d)\delta/\delta_c & \delta\leq \delta_c \\
      f_d & \delta > \delta_c 
\end{cases}
\end{equation}

Therefore, friction evolution is governing directly by slip evolution

\begin{equation}
\frac{df}{dt}=\begin{cases} 
      - (f_s-f_d)V/\delta_c & \delta\leq \delta_c \\
      0 & \delta > \delta_c 
\end{cases}
\end{equation}

where $V=d\delta/dt$. Furthermore, we can rearrange the stress equation introduced at the beginning of this notebook in terms of $V=d\delta/dt$ explicitly

\begin{equation}
\frac{d\delta}{dt}=\frac{2c}{\mu}\left(\tau_0-f\sigma_n - \frac{\mu}{2\pi}\int\limits_{-\infty}^\infty \frac{d\delta/dx'}{x-x'}dx'\right).
\end{equation}

where we are **requiring** that fault stress, $\tau$, at all times is **equal** to the (evolving) frictional strength, $f\sigma_n$.

Thus, we have a pair of **coupled ODEs** that can be solved according to whichever timestepping method you prefer. 

For example, an **improved Euler step**:

1. Specify friction, $f_0$, and slip, $\delta_0$, at time $t$.
2. Compute velocity, $V_0$, at time $t$, using $f_0$ and $\delta_0$.
2. Compute predictor friction, $f_1^P$, at time, $t+dt$, using $V_0$.
3. Compute predictor slip, $\delta_1^P$, at time, $t+dt$, using $V_0$.
4. Compute predictor velocity, $V_1^P$, at time $t+dt$, using $f_1^P$ and $\delta_1^P$.
5. Compute corrector velocity, $V_1^C=(V_0+V_1^P)/2$.
6. Compute corrector friction, $f_1^C$, and slip, $\delta_1^C$, at time $t+dt$, using $V_1^C$.

These equations have been implemented for you in the Python module `qk2.py`.

***Execute the cell below to create and visualise an initial stress state for your fault (shear and normal stress, static and dynamic frictions).***


```python
# CELL 1: CREATE INITIAL STRESS STATE

from qk2 import*

# create 1D mesh (power of 2 spacing for the FFT methods)
    # 5 km long fault
L = 5.e3
    # power of 2 grid points
x = np.linspace(0,L,2**10+1)
x = 0.5*(x[1:]+x[:-1])

# assign initial stress state
    # normal stress
sn = 30.e6
    # static and dynamic friction coefficients
fs,fd = 0.8, 0.6
    # constant, just below the static strength but above dynamic
s = (fs-0.01)*sn+0.*x             # the '+0.*x' is just a hacky way of getting the vector the right length...

# zero pad the stress vector on both ends to avoid FFT reflection problems
s[np.where(abs(x-L/2.)>L/8.)]=0

# create a fault object with this initial stress state
flt = Fault(x, s, fs=fs, fd=fd, normal=sn, dc=1e-2, dtn=30000, tmax=10.)

# plot the initial stress state
flt.show_stress(xlim = [L/4,3*L/4])
```

***Execute the cell below to trigger an earthquake on your fault. Visualise how slip, stress, and friction evolve over time.***


```python
# CELL 2: SIMULATE AN EARTHQUAKE (CRACK PROPAGATION) AND VIEW THE OUTPUT

# trigger an earthquake halfway along the fault
flt.trigger(x=L/2.)

# have a look at the earthquake
flt.show_eq(xlim = [L/4,3*L/4])
```

***Referring to the top left plot, explain why the crack (earthquake) stops propagating where it does.***

*your answer here*

***Referring to the right middle plot, identify periods of nucleation, propagation and arrest.***

*your answer here*

***Is the final slip profile consistent with expressions derived in the notes?***

*your answer here*

***Is the length of the nucleation period consistent with numbers derived in the notes?***

*your answer here*

***Experiment with different initial stress configurations and triggering locations.***

- **What happens if we don't trigger at the center of the positive stress drop region?**
- **What happens if we have multiple regions of positive and negative stress drop?**
- **What happens if the stress drop is larger?**

Hint: getting the right time output in these simulations is tricky. Experiment with the parameters `dtn` and `tmax`, while inspecting the right middle plot. You want to set these parameters not too much larger/longer than the length of the earthquake.


```python
# CELL 3: EXPERIMENT
# create 1D mesh (power of 2 spacing for the FFT methods)
    # 5 km long fault
from qk2 import*
L = 5.e3
    # power of 2 grid points
x = np.linspace(0,L,2**10+1)
x = 0.5*(x[1:]+x[:-1])

# assign initial stress state
    # normal stress
sn = 30.e6
    # static and dynamic friction coefficients
fs,fd = 0.8, 0.6
    # constant, just below the static strength but above dynamic
s = (fs-0.01)*sn+0.*x             # the '+0.*x' is just a hacky way of getting the vector the right length...

# zero pad the stress vector on both ends to avoid FFT reflection problems
s[np.where(abs(x-L/2.)>L/8.)]=0

# create a fault object with this initial stress state
flt = Fault(x, s, fs=fs, fd=fd, normal=sn, dc=1e-2, dtn=30000, tmax=10.)

# plot the initial stress state
    # OPTIONAL
#flt.show_stress(xlim = [L/4,3*L/4])

# trigger an earthquake halfway along the fault
flt.trigger(x=2.5e3)

# have a look at the earthquake
flt.show_eq(xlim = [L/4,3*L/4])
```

## Timestep control

The physics of crack nucleation and propagation occur on substantially **different time scales**, which presents a challenge to time stepping methods. Should we set the timestep small, to **accurately capture** rapid evolution during crack propagation, but take an **unnecessarily large number** of timesteps during nucleation? Or is it better to **step rapidly** through nucleation, and **resolve poorly** the rupture phase?

Ideally what we'd like is a timestep that changes after each improved Euler update, in a way that responds to the system evolution.

### The Proportional-Integral-Derivative (PID) controller

This a control feedback mechanism used to **adjust a process toward some desirable setpoint** on the basis of a computed error value, $e(t)$. The key is that the imposed adjustment (in this, changing the timestep) depends both on the size of the error, $e$ (the **proportional** part), its recent history, $\int e dt$, (the **integral** part) and how rapidly it is changing, $de/dt$, (the **derivative** part).

Fortunately, an improved Euler step embeds an error estimate within itself.

\begin{equation}
e(t) = \max\left\lvert\frac{U_1^C-U_1^P}{U_1^C}\right\rvert
\end{equation}

The timestep is then updated by multiply by a ratio that depends on the PID terms

\begin{equation}
dt_{i+1} = dt_i \left(\frac{\epsilon_r}{e_i}\right)^\alpha \left(\frac{e_{i-1}}{e_i}\right)^\beta \left(\frac{e_{i-1}^2}{e_{i-2}e_{i}}\right)^\gamma
\end{equation}

where $\epsilon_r$ is the relative error tolerance.

Taking the log of both sides

\begin{equation}
\log(dt_{i+1}) = \log(dt_i)+\alpha \left[\log\epsilon_r-\log e_i\right] + \beta \left[\log e_{i-1}-\log e_i\right]+\gamma\left[-\log e_{i-2}+2\log e_{i-1}-\log e_i\right].
\end{equation}

We can recognize:

- $\log\epsilon_r-\log e_i$, the integral term, the total error at time $t$.
- $\log e_{i-1}-\log e_i$, the proportional term, the accumulated error at this time step.
- $-\log e_{i-2}+2\log e_{i-1}-\log e_i$, the derivative term, the rate at which the error is changing.

The powers $\alpha$=0.17, $\beta$=0.245 and $\gamma$=0.05 are selected based on analysis and experience.

483db2d0e82b6a5f5877c01aff893e4e45321e70
ipynb
Jupyter Notebook
eq-sim.ipynb
ddempsey/ENGSC741_2019
dd0060c6c1042cdc435456c6fe14059cc5020240
MIT
eq-sim.ipynb
ddempsey/ENGSC741_2019
dd0060c6c1042cdc435456c6fe14059cc5020240
MIT
eq-sim.ipynb
ddempsey/ENGSC741_2019
dd0060c6c1042cdc435456c6fe14059cc5020240
MIT
Qwen/Qwen-72B
1. YES
2. YES