A Novel Approach for Direct Numerical Simulation of Hydraulic Fracture Problems

Abstract

Hydraulic fracturing is a non-linear and multi-physics problem involving the break up of a solid medium due to the action of hydrodynamic forces. Fluid and solid mechanics are involved at the same time together with fracture mechanics. Despite its relevance in many scientific and engineering fields, the theoretical and numerical description of hydraulic fracturing remains a challenging matter and the capabilities of existing models for applications are still limited. In this context, we propose a novel numerical approach to the Direct Numerical Simulation of hydraulic fracturing based on the Navier–Stokes equations coupled with peridynamic theory of solid mechanics through a multi-direct Immersed Boundary Method. The main advantage of this approach consists in the reliable crack-detection and tracking capabilities of peridynamics together with the capability of the Immersed Boundary Method of managing no-slip and no-penetration boundary conditions on complex and time-evolving interfaces. A massive-parallel solver based on this model has been implemented. We present a detailed theoretical description of the proposed methodology as well as the results of an extensive validation campaign for the new solver. Different benchmarking tests are provided together with the qualitative results of a simulation reproducing the fracture of a solid structure in a laminar, unsteady flow.

Appendix: Interface Detection Algorithm

In a Lagrangian framework, a discretized continuum is represented as a set of finite-size material points. In this frame, the detection and the tracking of surfaces and interfacial regions is not a trivial problem. A possible solution consists in tracking the interfaces by means of a set of Lagrangian markers, distributed on the interfaces themselves and moving solidly with them. In the peridynamic frame, the natural choice for interface markers consists in a specific subset of peridynamic particles selected within the overall set of material particles used for the solid discretization, namely interfacial particles. In other words, the Lagrangian points used for interface detection and tracking coincide with the Lagrangian positions of a specific set of solid particles, located in the proximity of the interface. This choice is particularly convenient for different reasons. The former relies in the fact that it does not require the usage of an additional set of dummy Lagrangian tracking points. Moreover, interfaces are automatically tracked by evolving the Lagrangian governing equations of peridynamics (9). Finally, as new interfaces are generated due to the formation of cracks, the set of Lagrangian markers can be easily extended. Indeed, when the size of the gap between solid fragments created by a crack overcomes a threshold value, the solid particles located at the opposite sides of the gap are automatically detected as interfacial ones and their Lagrangian position added to the list of interface markers. The list of interfacial particles is updated at every time-step of the solid solver cycle. All the external loads acting on solid object are assumed to be applied to the set of the interfacial particles, only.

Fig. 10

Two-dimensional schematic representation of the geometrical criterion used for interface detection with $N_s=8$, $N_v=2$, $r_{avg,t} = 2{\varDelta }_s$. The space around each material particle is subdivided into eight angular sectors. a A material particle is located on the interface if two or more adjacent angular sectors do not contain any other material particles interacting with the considered one. b A material point is located in the solid interior in all the other cases. Panels (c) and (d) show a sketch of the outcomes provided by the interface detection algorithm applied to the plate described in Sect. 3.1. The black particles are interfacial ones while the grey particles are the interior ones. The simulation parameters are the same used in Sect. 3.1 in except the applied load that is $\sigma =24\,\hbox {MPa}$. c Crack path at $t = 25\,\upmu \hbox {s}$. d Crack path at $t = 50\,\upmu \hbox {s}$

The numerical detection of interfacial particles is achieved via a geometrical criterion. The basic idea consists in assuming that a solid particle is located on an interface if exists a region of the space surrounding it that is completely depleted of other material points interacting with the considered one. The description of the algorithm is limited to the two-dimensional case for a clearer explanation. The numerical procedure is based on a loop over the whole set of peridynamic particles used for the discretization of the solid media. For each particle, ${\varvec{X}}_h$, a circular region of space of radius $R_s=\delta$ centered on the particle itself is divided into $N_s$ circular sectors of amplitude ${\varDelta }\theta _s = 2\pi /N_s$, as sketched in Fig. 10a, b. The radius of the sectors used by the algorithm is finite and is set to be equal to the peridynamic horizon. Indeed, cracks form due to the rupture of bonds that, by definition, always occurs within the horizon of a peridynamic particle. In this frame, the characteristic length of the smallest detectable interface is of the order of the peridynamic horizon. A second loop is used to scan the particles, ${\varvec{X}}_l$, pertaining to the family of the central one, ${\varvec{X}}_h$. Each of these material particles is associated with one of the circular sectors partitioning the horizon of ${\varvec{X}}_h$ and the status of the pairwise interaction ${\varvec{X}}_h \leftrightarrow {\varvec{X}}_l$ is checked by means of the parameter $\lambda _{h,l}$ defined in Sect. 2.1. For each circular sector, k, a mean relative distance, $r_{avg,k}$, between the central particle, ${\varvec{X}}_h$, and the particles ${\varvec{X}}_l$ for which $\lambda _{h,l}=0$ is computed:

$$\begin{aligned} r_{avg,k}=\frac{\sum _{l=1}^{N_{k}} ||({\varvec{X}}_l - {\varvec{X}}_h)|| (1-\lambda _{h,l})}{\sum _{l=1}^{N_{k}}(1-\lambda _{h,l})}. \end{aligned}$$

(27)

The summation is taken separately over each of the sectors partitioning the horizon of ${\varvec{X}}_h$ with $N_k$ the number of peridynamic particles pertaining simultaneously to the horizon of ${\varvec{X}}_h$ and to the sector k. The interface detection criterion is based on the assumption that the central particle, ${\varvec{X}}_h$, is an interfacial particle if one or both of the following conditions are satisfied:

• at least $N_{v}$ arbitrary consecutive sectors are depleted of particles, as sketched in Fig. 10a;

• all the pairwise interactions, ${\varvec{X}}_h \leftrightarrow {\varvec{X}}_l$, in at least $N_{v}$ arbitrary consecutive sectors have been suppressed and the distances $r_{avg,k}$ in each of these sectors are greater than a prescribed threshold value, $r_{avg,t}$.

In all the other cases the particle is assumed to be part of the solid interior, as sketched in Fig. 10b. The second condition of the two above avoids the spurious detection of interfacial particles in regions where the material is only partially damaged. Moreover, a proper setting of $r_{avg,t}$ can mitigate the problem of inaccurate evaluation of fluid forces in the narrow fluid gaps between solid fragments, with a size smaller than that of a given threshold. The values of $N_s$, $N_v$ and $r_{avg,t}$ should be carefully chosen since they can affect the accuracy of the surface detection procedure. In a peridynamic frame, the setting of these values have to be related to the resolution of the peridynamic discretization, $m=\delta /{\varDelta }_s$. In our experience $N_s=8$, $N_v=2$ and $r_{avg,t} = 3{\varDelta }_s$ is the proper setting for an accurate detection of the interfaces in the range $m=2$ to 4. Figure 10c, d shows a sketch of the outcomes provided by the interface detection algorithm applied to the plate described in Sect. 3.1. The simulation parameters are the same used in Sect. 3.1 in except the applied load that is $\sigma =24\,\hbox {MPa}$. The criterion can be extended to the three-dimensional case by using thirty-two spherical sectors and an analogous scanning procedure. The computational cost of the two-dimensional algorithm is of the order of $O({ NMN}_s)$ with N the total number of peridynamic particles used for the discretization of solids and M the average number of particles in each peridynamic family.