Numerical Experiments of State-based Peridynamics and Influence Functions in Classical Problems
Alejandro F. Queiruga
George J. Moridis
June 7, 2018
Lawrence Berkeley National Lab
Bottom Line Upfront
Peridynamics does not pass tests for a mechanics discretization:
- It does not represent constant strain modes exactly.
- The only converging formulation is unstable, especially around fracture tips.
- We can fix the instability, but that doesn't make #1 exact.
git clone my implementation PeriFlakes:
https://github.com/afqueiruga/PeriFlakes
The Jupyter notebooks replicate the graphs in this presentation.
Today, I'm going to talk about the constant strain modes.
Outline
- The problem I was trying to solve
- The simulation I managed to program
- How do we analyze this method?
- PeriFlakes
- Constant strain problems
- Boundary conditions
- Briefly: Linear fracture mechanics
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.
What is Peridynamics?
- We can describe linear mechanics as an operator $L$ in a continuum.
- PD defines an operator that depends on a damage $\alpha$: \begin{equation} \rho \ddot{\mathbf{u}} = L(\alpha)\mathbf{u} + \mathbf{b} \end{equation} with evolution laws for the damage, \begin{equation} \dot{\alpha} = \min{\left(f(\alpha,\mathbf{u}), 0\right)} \end{equation}
- Only considering static 2D plane strain linear elasticity: \begin{equation} 0 = L_{\alpha=1} \mathbf{u} \end{equation}
- Considering it only as a discretized operator.
What the properties of the Peridynamics specification of $L(\alpha)$?
How do we know our codes are right?
- Unit tests against certain solutions
- Convergence to known analytical solutions
- Convergence to manufactured solutions
- Convergence to super-refined self solutions
For a given method, some problems are in an exactly-representable solution space. E.g., Quadratic finite elements solve quadratic solutions exactly.
Which solutions are exactly representable by Peridynamics?
Constant Strain Problems
- Three basic modes of deformation with a constant gradient
- Fundamental solutions from basic mechanics
- FEM is correct to within machine precision
- This is the patch test.
These are Unit Tests in all of my codes! (PD, FEM, FDM, RKPM, etc.)
self.AssertTrue( error < 1.0e-12 )
Thin Pressurized Crack
Also check for convergence to more complicated analytical solutions.
\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 against the analytical solution $u^a$ to judge the 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}
Fixed-$\delta$ Convergence
Add more grid points keeping $w(r;\delta)$ fixed.
The support size grows with $m$
Periflakes/computation.ipynb
Larger $\delta$ vs $h$ means more points interact.This means the matrix fills up...
\begin{equation}
m=1.5
\end{equation}
Periflakes/computation.ipynb
This means the matrix fills up...
\begin{equation}
m=2.5
\end{equation}
Periflakes/computation.ipynb
This means the matrix fills up...
\begin{equation}
m=3.5
\end{equation}
Periflakes/computation.ipynb
This means the matrix fills up...
\begin{equation}
m=4.5
\end{equation}
Even if we don't assemble a matrix, this still shows data access patterns
of the program.
The number of nonzeros grows:
Even if we don't assemble a matrix, this still shows data access patterns
of the program.
Periflakes/computation.ipynb
Fixed-$m$ Convergence
- Changing $h$ but not $\delta$ changes the stencil.
- Refining $h$ but not $\delta$ leads to intractable computation!
- By keeping $m = \frac{\delta}{h}$ constant, the stencil stays the same.
- $m$ is the discretization-independent (hyper)parameter of the numerical method.
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:
- Formulation $\mathbf{T}\left[x\right]\left<\xi\right>$- discrete choice
- Influence function $w\left(r\right)$ - discrete choice
- Relative support radius, $m=\delta/h$ - continuous
- Boundary condition - discrete choice
- Stabilizer parameter $A$ - continuous
- 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.
Variety of Formulations
It is well known that bond-based does not work due to $\nu$ problem.
The are three classes of state-based models:
- The ordinary dilation-state methods
- Silling's
- Oterkus's
- The deformation gradient (F) correspondence methods
- Silling's original formulation
- Littlewood 2010/11's stabilization
- Silling 2017's stabilization (back-to-back in CMAME to Queiruga, 2017!)
- Post-processing smoothing
- Non-ordinary methods
- I haven't tried any yet.
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.
Influence Functions
$w(\mathbf{r})$ is another arbitrary choice!
\begin{equation} w(x'-x)=0 \quad\text{for}\quad x' \notin H(x) \end{equation}
I tried all of these (and more):
| Name | $w(r)$ in ${H}(x)$ |
|---|---|
| Const | 1 |
| Inv | $\delta/r$ |
| Linear | $1-\frac{r}{\delta}$ |
| Quadr | $\left(1-\frac{r}{\delta}\right)^2$ |
| Cubic | $\left(1-\frac{r}{\delta}\right)^3$ |
Different meshless methods have requirements on their radial kernel. Does PD have any?
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
Code Generation
The file peri_kernels_pop.py contains object-oriented implementations like this:
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:
formulations = {
'Silling' : StateBased,
'Oterkus2': Oterkus,
'Fbased':Fbased,
'Fstab_Littlewood':Fstab_Littlewood2010,
'Fstab_Littlewood2011':Fstab_Littlewood2011,
'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 36 different methods with only 240 lines of code!
Constant Strain Modes
Queiruga, 2017
Silling's Method Search
Queiruga, 2017
Oterkus's Method Search
Queiruga, 2017
F-based Method Search
Queiruga, 2017
Convergence Orders
How quickly do those methods converge to the true solution:
| isotropic | shear | uniaxial | |
| Fbased | 0.998518471486 | 0.88580762397 | 0.961508949756 |
| Oterkus2 | -0.0934425713562 | 0.992619961121 | -0.133199906094 |
| Silling | 0.940215702548 | -0.155718080685 | -0.120152361872 |
- Best results out of the exhaustive search
- F-based has $O(h^1)$ convergence, which is okay but slow!
- The dilation-based models do not converge to the correct result for
all three cases.
- Silling's method defines equivalance by a trial dilation
- Oterkus's method defines equivalence by a trial shear
- Cubic cone is the best influence function
Surface Effects
The $x$ component of Silling's with different $m$s:
PeriFlakes/surface_corrections.ipynb
Surface Effects
The $x$ component of Silling's with different $m$s:
PeriFlakes/surface_corrections.ipynb
Surface Corrections
Tried different ways of handling boundary conditions:
- No buffer zone \begin{equation} \mathbf{b} = \bar{\mathbf{t}}\frac{L}{A}, \quad \mathbf{u} = \bar{\mathbf{u}}, \quad\text{on} \,\Gamma_{bound} \end{equation}
- A buffer of size $\delta$ that integrates the PD operator, with the same BCs applied everywhere \begin{equation} \mathbf{b} = \bar{\mathbf{t}}\frac{L}{A}, \quad \mathbf{u} = \bar{\mathbf{u}} \quad \text{in}\, \Omega_{fict} \end{equation}
- A buffer with a finite difference for Neumann BCs: \begin{equation} \frac{u^f_y-u^o_y}{\Delta y} + \nu \frac{u^+_y-u^-_y}{2 \Delta x} = \frac{t_y(1-\nu^2)}{E} \end{equation}
- A buffer that integrates the PD operator, and then apples the BC constrain.
Surface Effects
- Boundary conditions are hard for all numerical methods, especially strong-form
- Every constant is computed numerically at each point in PeriFlakes
- Takes into account domain edge and $\alpha$ line-of-sight
- Described in Mandenci's 2014 book
- Le and Bobaru, 2018 tried many and found that ficticious nodes worked
Use a finite difference stencil: \begin{equation} \frac{u^f_y-u^o_y}{\Delta y} + \nu \frac{u^+_y-u^-_y}{2 \Delta x} = \frac{t_y(1-\nu^2)}{E} \end{equation} \begin{equation} \frac{u^f_x-u^o_x}{\Delta y} + \nu \frac{u^+_x-u^-_x}{2 \Delta x} = \frac{2 t_x (1+\nu)}{E} \end{equation}
The finite difference in vector notation:
Build the deformation gradient from differences, \begin{equation} \mathbf{F}=\frac{\left(\mathbf{u}^{f}-\mathbf{u}^{o}\right)}{\left|\mathbf{x}^{f}-\mathbf{x}^{o}\right|}\otimes\mathbf{n}+\frac{\left(\mathbf{u}^{+}-\mathbf{u}^{-}\right)}{\left|\mathbf{x}^{+}-\mathbf{x}^{-}\right|}\otimes\mathbf{e}_{t} \end{equation} and then apply the BC constraint by \begin{equation} 0=\lambda\mathrm{tr}\left[\varepsilon\right]\mathbf{n}+2G\varepsilon\mathbf{n}-\bar{\mathbf{t}} \end{equation} with $\varepsilon=\frac{1}{2}\left(\mathbf{F}^{T}+\mathbf{F}\right)$.
Equivalent to the stencil in Le and Bobaru, but rotates. Replace $-$ with $o$ for a corner.
Surface Corrections
PeriFlakes/surface_corrections.ipynb
F-based
- Adding the fictitious domain with same BCs still converges
- Stencil doesn't "catch" due to instability in operator.
PeriFlakes/surface_corrections.ipynb
Silling's dilation
- Adding the fictitious domain with same BCs still converges
- Stencil doesn't converge where it should. (Not 100% confident.)
PeriFlakes/surface_corrections.ipynb
Oterkus's dilation
- Adding the fictitious domain with same BCs still converges
- Stencil doesn't converge where it should. (Not 100% confident.)
PeriFlakes/surface_corrections.ipynb
Grid Convergence on Fracture
- The smoothing step is the best stabilization
- Silling's dilation is the best, but recall that it doesn't pass constant strain!
Periflakes/analysis.ipynb
Conclusions
- Need to do rigorous verifications of our numerical programs.
- Evaluate programs on fundamental problems looking for expected behavior.
- Perform hyperparameter sweeps---every choice is arbitrary.
| method | Uni. | Iso. | Shear | LFM | Stable |
|---|---|---|---|---|---|
| Silling | ❌ | ❎ | ❌ | ✅, #1 | ✅ |
| Oterkus | ❌ | ❌ | ❎ | ❌ | ✅ |
| Fbased | ❎ | ❎ | ❎ | ✅ | ❌ |
| Smoothed F | ❎ | ❎ | ❎ | ✅,#2 | ✅ |
- Can fix the stability and improve the accuracy of the F-based correspondence Peridynamics model.
- Silling's dilation model is the best on the crack, but it doesn't pass all constant-strain tests.
- BCs need work, but the error is bounded by $O(h^1)$.
- The operator itself has restrictions and needs more work.
Thank you
Email: afqueiruga@lbl.gov
Checkout my codes at
https://github.com/afqueiruga/PeriFlakes
and look at the Jupyter notebook and Zenodo database.
We can churn out and test more Peridynamics models in PeriFlakes!
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.
Periflakes/analysis.ipynb
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!
Periflakes/analysis.ipynb