Alejandro Queiruga's Presentation

## Smoothing Methods to Address the Numerical Stability of Peridynamics Near the Fracture Tip in Hydraulic Extension
### Alejandro F. Queiruga  
### George J. Moridis
#### June 7, 2018
#### Lawrence Berkeley National Lab

---

## Bottom Line Upfront

I found Peridynamics to be unsuitable:
1. It does not represent constant strain modes exactly.
2. The only acceptable formulation is unstable, especially around fracture tips.
3. We can fix the instability, but that doesn't change #1.
3. It is more computationally costly than the Finite Element Method.

Checkout my implementation with PeriFlakes:
[https://github.com/afqueiruga/PeriFlakes](https://github.com/afqueiruga/PeriFlakes)  
The Jupyter notebook replicates the graphs in this presentation.

Today, I'm going to talk about the path I took to make those
determinations, and discuss the stabilization around fracture tips.

---

## Outline

- The problem I was trying to solve
- The simulation I managed to program
- PeriFlakes
- *Briefly: constant strain problems*
- Thin crack benchmark problem
- Smoothing methods
- Comparisons including FEM

---

## Motivating Problem

Shale natural gas, Enhanced geothermal, Escape of gasses from nuclear
waste,  Bubbling of methane through oceanic sediments

---

## Peridynamic Representation

---

## Computational Framework

---

## Simulation Results

---

## Natural Fracture Interaction

---

## Videos are not results

It looks good, but not under actual scrutiny:

>Queiruga, A. F. and G. J. Moridis, "Numerical experiments on 
 > the convergence properties of state-based peridynamic laws and 
>  influence functions in two-dimensional problems." Computer 
 > Methods in Applied Mechanics and Engineering 322 (2017): 
 > 97-122.

**It is critical to verify and validate numerical methods on known
solutions and expected discretization refinement behavior.**

- Can't jump to dynamics and fracturing simulations immediately.

---

## Variety of Peridynamic Formulations

It is well known that bond-based does not work due to $\nu$ problem.

The are two classes of state-based models:

1. The dilation-state methods
  1. Silling's
  2. Oterkus's
2. The deformation gradient (F) correspondence methods
  1. Silling's original formulation
  1. Littlewood 2010's stabilization
  2. Silling 2017's stabilization (back-to-back in CMAME to Queiruga, 2017!)
  4. I tried post-processing smoothing

---

## Hyperparameters

Our goal is to obtain a computer program that solves the problems we
want accurately and efficiently.

In the Peridynamics scope, we have the following hyperparameters:

1. Formulation $\mathbf{t}$- discrete choice
2. Influence function $w$ - discrete choice
3. Relative support radius, $m=\delta/h$ - continuous
4. Stabilizer parameter $A$ - continuous
5. Smoothing function $w^s$ - discrete choice (including None)

These are **arbitrary**. We need to determine which combinations are
acceptable, and which of those are the best.

---

## PeriFlakes

[https://github.com/afqueiruga/PeriFlakes](https://github.com/afqueiruga/PeriFlakes)

- A cleaned up implementation of the code for *Queiruga, 2017*
- Released under GPL v3
- Written using my Cornflakes/Popcorn domain specific language
- Autogenerates code for new methods, including tangent matrices
- Jupyter notebook for all the plots in this presentation

- Databases of solutions archived at:
> Alejandro Francisco Queiruga. (2018). 2D Peridynamic Displacement Fields [Data set]. Zenodo. http://doi.org/10.5281/zenodo.1284634

(*I'm still working; it'll be 100% done by WCCM*)

## Code Generation

The file `peri_kernels_pop.py` contains object-oriented implementations like this:
```python
class Oterkus(StateBased):
    def m_expr(self):
        return self.w0I * alpha_I * rxabs
    def a_expr(self):
        return self.w0I * alpha_I * rxabs**2
    def b_expr(self):
        return self.w0I * alpha_I * (xn[0] * xn[1] )**2 * rxabs**2
    def theta_expr(self):
        return 2/m[0] * self.w0I * alpha_I * (ryabs-rxabs) *  (yn.T*xn)[0,0]
    def force(self):
        a = (c_K - 2*c_G*a_integ[0]/b_denom[0])/2
        b = c_G /( 2* b_denom[0])
        A = 2 * self.w0I* alpha_I  * ( a*2/m[0]*(yn.T*xn)[0,0]*theta[0] + b*(ryabs-rxabs) )
        return A * yn
```
		
- Vector symbolic notation using Sympy.
- Symbolic differentiation for tangent
- Generates a C module for matrix/vector contribution.

## Code Generation

We want to exhaustively search for methods:
```python
formulations = {
    'Silling' : StateBased,
    'Oterkus2': Oterkus,
    'Fbased':Fbased,
    'Fstab_Littlewood':Fstab_Littlewood,
    'Fstab_Silling':Fstab_Silling,
}
delta = i_delta[0]
weight_funcs = {
    'const' : lambda r : Piecewise( (1.0    , Le(r,delta)), (0.0,True)),
    'inv'   : lambda r : Piecewise( (delta/r, Le(r,delta)), (0.0,True)),
    'linear': lambda r : Piecewise( ((1.0-r / delta)**1.0, Le(r,delta)), (0.0,True)),
    'quadr' : lambda r : Piecewise( ((1.0-r / delta)**2.0, Le(r,delta)), (0.0,True)),
    'cubic' : lambda r : Piecewise( ((1.0-r / delta)**3.0, Le(r,delta)), (0.0,True)),
    'quarticA':lambda r: Piecewise( (1.0 - 6.0*(r/delta)**2 + 8.0*(r/delta)**3 - 3.0*(r/delta)**4, Le(r,delta)), (0.0,True))
}
for methodname,c in formulations.iteritems():
    for weightname,f in weight_funcs.iteritems():
       c(methodname+"_"+weightname, f).kernel()
```
This generates 30 different methods with only 240 lines of code!

---

## Constant Strain Modes

<img src="assets/figures/blocks.png" style="max-height: 30vh" >
<table><tr><td></td><td>isotropic</td><td>shear</td><td>uniaxial</td></tr><tr><td>Fbased</td><td>0.998518471486</td><td>0.88580762397</td><td>0.961508949756</td></tr><tr><td>Oterkus2</td><td>-0.0934425713562</td><td>0.992619961121</td><td>-0.133199906094</td></tr><tr><td>Silling</td><td>0.940215702548</td><td>-0.155718080685</td><td>-0.120152361872</td></tr></table>

- Best results out of an exhaustive search
- F-based has $O(h^1)$ convergence, which is slow!
- The dilation-based models do not converge to the correct result for
  all three cases.
- FEM is correct to within machine precision

## Dilation Method Search

<center>Silling's on the left, Oterkus's on the right.
<div align="center">
<img src="assets/figures/conv_ouchi_all.png" style="float:left; height:70vh"><img src="assets/figures/conv_oterkus_all.png" style="float:left; height:70vh">
</div>
</div>
</center>

## F-based Method Search

---

## Thin Pressurized Crack

We need to check for convergence to analytical solutions.

<img src="assets/figures/crack_diagram.png" style="float:left; height: 60vh" >
\begin{equation}
\hat{Z}\left(z = x+iy\right)=P\left(\sqrt{z^{2}-L^{2}}-z\right)
\end{equation}

\begin{equation}
u\_{x} = 
\frac{1}{2G}\left(\frac{2-4\nu}{2}\mathrm{Re}\left[\hat{Z}\right]-y\mathrm{Im}\left[Z\right]\right)
\end{equation}
\begin{equation}
u\_{y} = \frac{1}{2G}\left(\frac{4-4\nu}{2}\mathrm{Im}\left[\hat{Z}\right]-y\mathrm{Re}\left[Z\right]\right)
\end{equation}

Calculate an error to judge numerical solution $u^N$:
\begin{equation}
e = \frac{ \int \left\|u^N-u^a \right\|^2 \mathrm{d}x }{ \int
\left\|u^a \right\|^2 \mathrm{d}x }
\end{equation}

---

## Peridynamic Displacement Field

The F-based method's instability effects bond stretches, making the fracture branch
and jump nonphysically.

---

## Smoothing

My ad hoc approach was to post process the solution:
\begin{equation}
\mathbf{u}^s = \frac{\int\_\mathcal{H} w^s \alpha \mathbf{u} \mathrm{d}x'}{\int\_\mathcal{H} w^s \alpha \mathrm{d}x'}
\end{equation}

The smoothing kernel $w^s$ doesn't have to be the same $w$.

This worked to stabilize the fracture path in my original hydraulic fracturing simulations,
but the method never passed any other validations!

---

## The Stabilizers

We can add a stabilizer to the force:
\begin{equation}
\mathbf{T}\left<\xi\right> = w\alpha \sigma \mathbf{K}^{-1} \mathbf{\xi} + \mathbf{T}^s
\end{equation}

Silling's 2017 paper:
\begin{equation}
\mathbf{T}^s = w\alpha A \frac{18 K}{ w_0 \pi \delta^5} \left( \mathbf{Y} - \mathbf{F}\xi \right),
\quad w\_0 = \int_\mathcal{H} w\alpha \mathrm{d}x'
\end{equation}

Littlewood's 2010 paper:
\begin{equation}
\mathbf{T}^s = w\alpha A \frac{18 K}{\pi \delta^4}
\left( \frac{ \left|\mathbf{Y}\right| - \left|\mathbf{F}\xi\right| }{ \left|\xi\right| } \right)
\frac{\mathbf{Y}}{\left|\mathbf{Y}\right|}
\end{equation}

$A$ is a new penalty hyperparameter. It would be nice if it were
insensitive to the grid size.

---

## Smoothing solution:

Smoothing with a cubic function.

## Smoothing sweep

All of the functions are about the same.

---

## Silling's stabilizer:

Set $A = 2.0$

## Silling's stabilzer sweep

Looks grid-size invariant, best around 1.

---

## Littlewood's stabilizer:

Set $A=10.0$

## Littlewood's stabilizer sweep

Needed to sweep higher.

---

## Grid Convergence on Fracture

- The smoothing step is the best correction
- Silling's dilation is the best, but recall that it doesn't
  pass constant strain!

---

## What is our ultimate goal?

We're searching for a numerical method that:

- Passes all standard tests (continuum mechanics solutions)
- Solves the fracture problem efficiently
- Efficiently represents small changes in fractures
- Good tradeoff between accuracy and computational cost

Compare to a discretely meshed fracture solved by FEM:

- Meshed with GMSH, implemented with FEniCS
- Included in `PeriFlakes/fem/fem_crack.py`

Wall-clock time for matrix/vector assemble + solve.

---

## Error vs. Compute time

- FEM is an order of magnitude more accurate and faster.

---

## Minimal Representable Length

- Bond breaks are discrete events at fixed locations
- Peridynamic grid spacing is the minimal length a fracture can jump
- Same as a the mesh size in node-splitting or element-deleting FEM

---

## Minimal Representable Length

FEM solves a smaller length scale faster!

---

## Conclusions

- Need to do rigorous verifications of our numerical programs.
- Evaluate programs on a static crack problem with a known solution.
- Perform hyperparameter sweeps---every choice is arbitrary.
- Can fix the stability and improve the accuracy of the F-based
 Peridynamics model.
- Silling's dilation model is the best on the crack, but it doesn't
 pass all constant-strain tests.
- FEM is faster and more accurate.
- Peridynamics discretization isn't more efficient for a given
  fracture resolution.

We can churn out and test more Peridynamics models in PeriFlakes!

---

# Thank you

Email: afqueiruga@lbl.gov

Checkout my codes at  
[https://github.com/afqueiruga/PeriFlakes](https://github.com/afqueiruga/PeriFlakes)  
and look at the Jupyter notebook and Zenodo database.

I have another upcoming talk at WCCM in July:
> #2018172, Numerical Experiments on the Convergence Properties of State-based Peridynamic Laws and Influence Functions in Two-Dimensional Problems