Dynamics of fracturing saturated porous media and self-organization of rupture
C. Peruzzo D. T. Cao E. Milanese P. Favia F. Pesavento F. Hussain B. A. Schrefler
 

    Nuh Gedik
    January 22, 2019
====================



The paper was presented as a general lecture at the EUROMECH Solid Mechanics Conference in Bologna.

E-mail addresses: mailcarloperuzzo@gmail.com (C.Peruzzo), toancaoduc@gmail.com (D.T.Cao), enrico.milanese@epfl.ch (E.Milanese), pietro.favia86@gmail.com (P.Favia), francesco.pesavento@dicea.unipd.it (F.Pesavento), fazlehussain@gmail.com (F.Hussain), bernhard.schrefler@dicea.unipd.it (B.A.Schrefler).



  Analytical solutions and a vast majority of numerical ones for fracture propagation in saturated porous media yield smooth behavior while experiments, field observations and a few numerical solutions reveal stepwise crack advancement and pressure oscillations. To explain this fact, we invoke self-organization of rupture observed in fracturing solids, both dry and fully saturated, when two requirements are satisfied: i) the external drive has a much slower timescale than fracture propagation; and ii) the increment of the external load (drive) is applied only when the internal rearrangement of fracture is over. These requirements are needed to obtain clean Self Organised Criticality (SOC) in quasi-static situations. They imply that there should be no restriction on the fracture velocity i.e. algorithmically the fracture advancement rule should always be independent of the crack velocity. Generally, this is not the case when smooth answers are obtained which are often unphysical. Under the above conditions hints of Self Organized Criticality are evident in heterogeneous porous media in quasi-static conditions using a lattice model, showing stepwise advancement of the fracture and pressure oscillations. We extend this model to incorporate inertia forces and show that this behavior still holds. By incorporating the above requirements in numerical fracture advancement algorithms for cohesive fracture in saturated porous continua we also reproduce stepwise advancements and pressure oscillations both in quasi-static and dynamic situations. Since dynamic tests of dry specimens show that the fracture advancement velocity is not constant we replicate such an effect with a model of a debonding beam on elastic foundation. This is the first step before introducing the interaction with a fluid.




§ INTRODUCTION

The term “self-organization of rupture” was first coined by Tzschicholz and Herrmann [1] when investigating pressure fluctuations and acoustic emission in hydraulic fracturing. According to Herrmann [2] it expresses an experimentally observed fact in disordered media where systems of cracks organize themselves to create a complex structure (often self-similar). Examples include localization of shear bands out of micro-cracks, formation of a drying cracks network, fractal cracks in hydraulic fracture or stress corrosion cracking and tectonic fault systems. By simulating hydraulic fracturing on a beam lattice model Tzschicholz and Herrmann found that the sequence of breaking events is organized in bursts, also called avalanches, that have power law distributed lifetime and quiescence intervals [3-10]. This indicates that self-organized criticality (SOC) [1] takes place or at least implied for samples of finite dimensions. Ignoring the requirements for self-organized criticality has serious consequences when designing the fracture advancement rule for analytical or numerical simulation of fracturing saturated homogeneous porous media - both in quasi-static and dynamic situations.

Fracturing in fully saturated porous media is of great importance in geophysics [11-15] and in non-conventional oil and gas extraction [16-18]. In geophysics both mechanical [12] action and pressure induced fracture [11, 19] are the case while in oil and gas extraction the latter is prevailing. A signature of hydraulic fracturing through lava is e.g. the magmatic dikes radiating from West Spanish Peak, Colorado, U.S. [20]. Despite the fact that in pressure induced fracture in geological situations fracture advancement speed can reach 1/5th of the speed of sound [11] and slugs of magma can rise in rocks up to a speed of 17m/s [21], inertia forces are usually neglected in simulations. Some solutions including inertia forces exist, however [22-25]. [22,23] show continuous fracture growth, [24] some pressure oscillations due to fracture toughness contrast between two adjacent layers, while [25] clearly exhibits stepwise behavior and pressure oscillations. This aspect will be addressed below and involves the use of algorithms which respect “self-organization of rupture”.

The many existing closed form and numerical solutions, see [25, 26] for extensive lists, are obtained mainly in a quasi-static context and a homogeneous setting. Stochastically heterogeneous fracturing media have been taken into account via statistical analysis in [1] as already mentioned, but without flow in the domain. Flow in the fracture and in the domain, has been investigated in [27] on a central force model for mechanical load and pressure induced fracture. Also, in these cases the fracture does not propagate in a continuous manner, but advances stepwise and the avalanches size (number of failing elements per loading step) follows a power law. However, it has also been found in [27] that in pressure induced fracture a power law exists only for low injection rates, while for higher injection rates the power law behavior breaks down, see Figure 1.


This hints at a transition in the behavior regime depending on the injection rate, which requires further research. 

The stepwise fracture growth and ensuing pressure oscillations are well known in the oil and gas industry since the late 80ties [17,18], Figure 2. There it has important implications for the hydrocarbon production [28]. An experiment of fluid injection in Colton Sandstone [29] evidenced that for high viscosity fluid and low constant injection rate a continuous fracture propagation was obtained (see for example Fig. 3) while high constant injection rate and low viscosity (field conditions for fracking) resulted in a stepwise advancement with ensuing pressure oscillations, Figure 2. Again, there is a transition in behavior depending on the injection rate as well as the viscosity. In [25] we conjectured that the transition from continuous to staccato advancement hides a stability problem, which has still to be proved. 

Stepwise advancement of a mechanically loaded hydrogel, both in liquid and in air has also been evidenced experimentally [30]. The stepwise behavior was mainly ignored in the mechanics and computational mechanics community. In this community, a first solution evidencing staccato fracture advancement was shown in a quasi-static setting with a standard Galerkin Finite element method with remeshing [31]. Only very recently other solutions capturing the stepwise advancement appeared; in the most recent papers [32,33] the Extended Finite Element Method is used.



We concentrate here our attention on the solid part and the interaction with the fluid. The behavior of the fluid, especially in the fracture will require further scrutiny. After addressing the issue of self-organization of rupture and the ensuing consequences for the design of the fracture advancement rule we shall extend the central force model of [27] to also include inertia forces and shall show solutions for quasi-static and dynamic analyses in stochastically heterogeneous saturated porous media. We shall consider both mechanical and hydraulic loading. Since it has been shown [5] that inclusion of inertia forces does not drive the system away from criticality, in the dynamic case we will omit the investigation of avalanche behavior, i. e. we shall limit us to single runs. Then we will show solutions obtained with the standard Galerkin Finite Element Method for fracturing homogeneous media in dynamics, based on an algorithm respecting the requirements for self-organization of rupture. Finally, we show a closed form solution for dry materials dynamics, based on a cohesive model. This is a first step towards investigating the above-mentioned conjecture. Concluding remarks will close the paper.



§ SELF-ORGANIZATION OF RUPTURE
 

When studying fracturing of heterogeneous media with and without fluid with methods of statistical physics in quasi-static situations, the following two requirements have to be satisfied, apart from heterogeneity [1, 3, 27]



  * External drive has a much slower timescale than fracture propagation;

  * increments of the external load (drive) are applied only when the internal rearrangement of fracture is over.

By satisfying these conditions, power law behavior has been found for the avalanche size probability distribution (breaking or damaging events in a time step), distribution of time duration of each avalanche, energy bursts, the probability distribution of finding two breaking events to be τ time steps apart. This last one corresponds to Omori's well-known law of aftershocks for earthquakes first formulated in 1894 [34] which has been verified from earthquake catalogues for aftershocks series ranging from a few hours to a couple of years after the main event [35]. Tzschichholz and Herrmann [1] conclude that one might expect the basic mechanism for burst sequences, or, respectively, aftershocks, to be universal due to self-organization of rupture. 
Also, faster loading than required above is possible, but then we may not be able to discern individual avalanches anymore. Clearly the avalanche behavior of the breaking events entails pressure oscillations with irregular intervals and amplitudes, which persist also in homogeneous media [1,36]. Satisfying condition i) and ii) implies that there should be no restriction on the fracture velocity i.e. algorithmically the fracture advancement rule should always be independent of the crack velocity. This is completely ignored when imposing velocity restriction to obtain solutions. Such velocity restrictions are imposed for instance when searching for analytical solutions, or with inadequate crack tip advancement/time stepping algorithms in numerical solutions. For the first case the statement “the solution near the tip of a propagating finite hydraulic fracture is captured by a stationary solution for a semi-finite crack moving at a constant velocity V –a solution akin to a traveling wave solution” is typical and a whole body of literature exists, based on this restriction [26]. For the second aspects, the most important papers are reviewed in [25]. Clearly a velocity restriction does not allow the system to behave in a physically correct way and we risk getting flawed results.

The consequences are evidenced by the following example, taken from [37]. Fluid is pumped in a borehole at constant rate Q equal to 0.0001 m3/s causing the fracture to advance.


The following analytical solutions [38] for crack length L, crack mouth opening displacement (CMOD), and for the pressure at the crack mouth (Pcm) respectively comes from an asymptotic solution and an approximation to it, based on some simplifying assumptions; in particular the fluid is incompressible, fracture impermeable, and linear elastic fracture propagation takes place. Further, the fluid completely occupies the fracture volume, hence affecting the tip velocity and excluding fluid lag in the tip region.


    L=0.65(GQ^3/μ(1-ν))^1/6t^2/3


    CMOD=2.14(μ(1-ν)Q^3/G)^1/6t^1/3


    p_cm=1.97(G^3 Q μ/(1-ν)^3 L^2)^1/4+S


For the numerical solution, we have here adapted the cohesive model [31] to approach as much as possible linear elastic fracture and adopted a finite domain. Further the axial symmetry of the problem is broken by introducing a notch from which the crack enucleates. All the other simplifying assumptions of [38, 39] have been taken into account, Table 1. The numerical and analytical solutions are presented in Fig. 4. The crack mouth opening displacement matches pretty well, but there are some discrepancies with the fracture advancement although not too serious. The salient point is the non-regular stepwise advancement in the numerical solution. This is due to the fact that we use a crack-tip advancement/ time stepping algorithm shown in section 4 which respects the self-organization of rupture i.e. it imposes no tip velocity restriction. The pressure makes the difference.




After an initial period needed for starting, influencing the solution due to the initial condition chosen in the numerical model, the average pressure follows the exact solution, equation (3). While all analytical solutions [26] show smooth pressures versus time, in Figure 5 clearly downward pressure jumps are observed at each step of the fracture advancement which have been explained by Biot’s theory (second partial scenario of [27]): if flow is specified in the fluid continuity equation (second Biot equation, next section) the flow effect is transmitted to the solid through the pressure coupling term in the effective stress (first Biot equation). The solid is loaded and upon rupture produces a sudden increase of the volumetric strain and an ensuing drop in pressure. In quasi-static situations fluid lag does not matter much in this reasoning because it has been found [40] that in most practical cases fluid lags, if it exists at early stages, diminishes with the propagating fracture and may be ignored. In dynamics however, differences between the velocities of crack tip and fluid front in the fracture may be expected and may hence influence the phenomenon described. These pressure drops are observed both in experiments and in the field and make the solution relevant for steering fracking operation [28], and for investigating volcanic [15] and subduction tremor in geophysics [13]. In such situations, smooth solutions are perfectly useless and unphysical. In the next two sections we will show that these pressure jumps persist also in dynamic solutions.



§ DYNAMICS IN HETEROGENEOUS POROUS MEDIA


The central force model previously developed for quasi-static situations [27] will be briefly recalled in the next sub-section and extended to consider inertia forces. In the same section, we will also explain how the concept of the self-organization of rupture is built inside this model. In sub-section 3.2 we will present two applications. The first one with a central flux imposed (hydraulic loading), the second one with a pressure imposed at the centre (mechanical loading). These examples which respect the requirement ii) of section 2 show that:

	
  * there are still the pressure rises and drops at fracture, as in the quasi static case;
	
  * accounting for inertia will lead to a more physical description of the pressure oscillations and the fracture advancement.





 §.§ The model

The model is the continuous-damage fully saturated model (CDFSM) of [27] based on the generalized Biot’s theory [41,42] under the following assumptions:

	
  * linear elastic stress–strain relations;
	
  * small strains;
	
  * incompressible solid and fluid phases. 

The extension of the CDFSM to dynamics is here summarized. The velocity of the fluid v^w is defined taking the velocity of the solid skeleton v^s as a reference and introducing the relative velocity v^ws

    v^ws=v^w-v^s

The water acceleration a^w will be [43]:

    a^w=a^s+D^sv^ws/Dt+(v^s+v^ws)·v^ws=a^s+a^ws+(v^s+v^ws)·v^ws

The last term in (2) is called convective term and a^ws is the acceleration of the fluid with respect to the water. Following Zienkiewicz et al. [44], those terms are not taken into account. With these assumptions the linear momentum balance equation is as follows:

    -ρa^s +ρg +∇·σ =0

where ρ^s,ρ^w,ρ are respectively the intrinsic density (ID) of the solid phase, the ID of water and ρ=(1-n)ρ^s+nρ^w is the averaged density of the multiphase systems. n is the porosity. a^s is the acceleration of the solid matrix. ∇ is the divergence operator and σ the total stress tensor; g is the acceleration related to gravitational effects.
Terzaghi’s effective stress principle holds and it is written in the form:

    σ=σ'-δ p^w

with σ' the effective stress tensor, p^w the pressure and δ the Kronecker delta. The continuity of the solid and fluid phases, together with Darcy’s law for fluid flow in a porous medium, yield:

    ∂ε/∂t+[k/μ^w(-p^w+ρ^w(g-a^w)) ]=0

where ε is the volumetric strain, k and μ^w are the fluid intrinsic permeability and dynamic viscosity respectively (both assumed constant in space and time). Considering the water acceleration in the water mass balance equation will lead to a so called “dynamic seepage contribution”. Its effect is negligible, [44], and will not be considered.

Discretizing Eqs. (3) and (5) by means of the Galerkin Finite Element method, the coupled system of equations assumes the same form as in [27] except for the equilibrium equation of the mixture solid plus fluid (first equation) that now accounts for inertia term (6). The second equation is the continuity equation of the fluid (7).

    ∫_Ω^B^T σ' dΩ-Qp̅_w +Mü̅̈=f_u


    Hp̅_w +Q^T u̇̅̇=f_p

M is the mass matrix,B  is the strain operator,Q  the coupling matrix,H  the permeability matrix, σ'  the effective stress tensor responsible for deformation and strength,p̅_w  and u̅ respectively the nodal values of the pressure and of the displacements and f_u and f_p the loading terms. Dot means differentiation with respect to time.
∫_Ω^B^T σ' dΩ=Ku̅, Mü̅̈ and the other matrices have been derived for the three-layer element used in the following and are shown in Figure 6.


The system of equations (6) and (7) can be conveniently rewritten as:

    [ M 0; 0 0 ][    ü̅̈; p̈̅̈_w ]
    +
    [   0   0; Q^T   0 ][  u̇̅̇; ṗ̅̇_w ]
    +
    [  K -Q;  0  H ][   u̅; p̅_w ]
    =
    [ f_u; f_p ]
    ↓
    Mä+Cȧ +Ka= f

The time integration of the system (8) is carried out by the “Truncated Taylor series collocation algorithm GN22” [45]. The solution for the unknowns a=[   u̅ p̅_w ] at a generic time step t_n+1=t_n+Δ t is given by:

    a_n+1=-A^-1(f_n+1+Câ̇_n+1+Mâ̈_n+1)   
    
    A=2/β_2Δt^2M+2β_1/β_2ΔtC+K
    â̇_n+1=-2β_1/β_2Δta_n+(1-2β_1/β_2)ȧ_n+(1-β_1/β_2)Δtä_n
    â̈_n+1=-2/β_2Δt^2a_n-2/β_2Δtȧ_n-1-β_2/β_2ä_n

where β_2 ≥β_1 ≥ 0.5 are parameters for time integration.

In the CDFSM continuous damage is applied to the truss elements as follows. At the beginning of the simulation a uniform Young’s modulus of 100 MPa is assigned to each element. The spatial distribution of defect in the solid matrix (heterogeneous behavior) is taken in to account by randomly extracting a stress threshold from the range (0,1) MPa. Each threshold has the same probability of being extracted (uniform distribution). During the simulation, when the stress in a truss exceeds the local threshold, the elastic modulus of the truss is reduced by the factor (1-D) , thus considering damage. D is the fixed damage amount. The first extreme value D = 0 means no damage and the other extreme value D = 1 means that the truss is completely broken and does no more contribute to the global stiffness. To represent gradual damage, D is set to 0.1: each time that a truss is damaged, its new Young’s modulus value is 90% of the previous one. Each truss can be damaged up to thirty times before the final failure, thus conveying an asymptotically decreasing damage rule. The threshold is updated after each damage event (annealed behavior).

No limits are imposed on the number of elements that could be damaged within a time step. The next time step starts only when, taken an element i , the stress σ_i is less than the threshold σ^max_i for each i , see Figure 7.


The concept of the self-organization of rupture at the beginning of the section 2 is expressed in dynamics only by the requirement ii); while the requirement i) is implicitly satisfied. 


 §.§ Numerical applications

The two examples of this section have the same domain, boundary and initial condition shown in Figure 8. The examples differ from each other because of the central boundary condition: a constant pressure is imposed p(t)=1MPa in the example of 3.2.1, while a constant flux f(t)=3*10^-2mm^3/s is imposed in the example 3.2.2.

Physically one can think about those examples as a plane strain representation of a soil specimen where a fluid is either pumped with a constant and a variable fluid flow in the centre. No mechanical load is applied. In both cases the permeability is kept constant, either in the porous matrix and in the fracture path. 


  §.§.§ Water pressure imposed at the centre of the domain

This loading case is of mechanical type and is applied on the first Biot equation (8). In Figure 9, continuous line, both pressure rises and drops can be observed. This corresponds to the full scenario of [27], i.e. the second partial scenario of section 2 has to be completed by a first partial scenario which reads “if a load, pressure, or displacement boundary condition is applied suddenly (all these conditions acting on the equilibrium of the solid-liquid mixture), then the fluid bears initially almost all the induced load because its immediate response is undrained (rigid and non-flowing). Then through the coupling with the fluid, the overpressures decrease and the solid gets loaded”. This situation produces a pressure rise upon rupture. However, the sudden changes of displacements would induce also some pressure drop in this first partial scenario due to the volumetric strain and pressure rise in the second partial scenario; therefore, the effects of both scenarios will appear in general (whole scenario). In quasi-static case the two scenarios are well separated for mode I fracturing. In dynamics, this is not the case. In fact, when the sample is locally subjected to a sudden damaging process, or to a fracture advancement, part of the strain energy is released and converted to kinetic energy. Strain and displacement waves of finite velocity appear. During the simulation, we observe damage at multiple points, so we have multiple generation of strain waves the interaction of which can create conditions for new damage or new fracture advancement. This is influencing also the interaction between fluid and solid. The presence of these waves changes the behavior of the sample. Neglecting continuous damage would lead to a fast consolidation problem shown in Figure 9 with dashed lines while the continuous one represents the solution with continuous damage. Figure 10 shows the time evolution of the system at different time stations in terms of Young modulus (above) and stresses (below).




  §.§.§ Imposed Flux

In this second case a constant water flux is applied on the second Biot equation (8). Looking at Figure 11, we observe pressure both rises and drops: similarly to the previous section 3.2.1 we have the full scenario of [27]. The solution of the problem without the damage algorithm, dashed line in Figure 11, shows an increase of the pressure which slowly reaches a constant mean value (not shown in the Figure). Figure 12 shows the time evolution of the system at different time stations in terms of Young modulus (left) and stresses (right).


When the damage occurs at some position we mostly expect a pressure drop because of the new space created. This concept is proved in Figure 13, by comparing the state of the system before and after a pressure drop in the two pairs between steps 210 and 220 and steps 340 and 350. On the right-hand side, we see the evolution of pressure in some selected points (A-D), where the drops in the above intervals can be clearly seen.

Looking at the bottom pairs of pictures of Figure 13, the location of point 73 of Figure 11 and the related pressure evolution, we observe that at this node the overpressure is going to zero because of the nearby formation of a large fracture channel. On the contrary, at node 217 of Figure 11 (located in an un-fractured region, Fig. 13), the pressure is globally not decreasing even if it is still below the solution obtained without considering the fracture process. From Figure 12, we also note that on the boundary of an existing fracture, a new one has been nucleated and grows in a different direction from the previous one.

Note that the fracture width, unknown a priori, is bigger than the size of one cell of the lattice.





§ DYNAMICS IN HOMOGENEOUS POROUS MEDIA

In this section, we show a few results of two loading cases for fast dynamic fracturing in homogeneous saturated porous media taken from [25]. The following fracture advancement algorithm has been used [31] which respects the two requirements of section 2: a fracture arises when the Rankine criterion is not satisfied at a point. The fracture is then advanced by an increment Δ s in the direction normal to the maximum principal stress and the mesh is consequently re-built. 

In 2D the fracture follows directly the direction normal to the maximum principal stress while in 3D the fracture follows the face of the element around the fracture tip which is closest to the direction normal to the maximum principal stress; the fracture tip becomes a curve in space (front). If during the advancement a new node is created at the front the resulting elements for the fluid in the crack are triangular in 2D and tetrahedral in 3D. At each time t_n, j successive tip (front) advancements are possible within the same time step; see Fig. 14. Their number in general depends on the chosen time step increment Δ t, the adopted crack length increment Δ s, and the variation of the applied loads. This algorithm works also in dynamics and in [25] it has been shown that still several advancements within a time step happen despite the fast loading. For another algorithm, which avoids remeshing but respects the two conditions of section 2, see [32].

he two loading cases are shown in Figure 15. The saturated porous medium sample is a rectangle of 0.2 m height and 2 m length. In case of mechanical loading the sample is loaded in traction by two uniform vertical velocities with magnitude 2.35x102 m/s applied in opposite directions to the left end of the top and bottom edges. In case of pressure loading the same cross section is subjected to a fast pressure increase of 5x10^11 Pa/s at the crack mouth reaching a maximum value of 100 MPa. Vertical and horizontal displacements are constrained at the right edge, and the boundary is impervious and adiabatic. For material parameters see [25]. In both cases stepwise advancement and pressure oscillations are obtained but the frequency content of these oscillations differ completely between the two loading cases: while for mechanical loading a broad spectrum appears with peaks at 2 kHz, 12 kHz and 18 kHz , in case of hydraulic loading a single peak prevails, see also [13, 14].

We show in for mechanical loading Figure 16 a few significant snapshots for the pressure wave propagation at respectively (a) 0.02025 s; (b) at 0.02035 s; and (c) at 0.0205 s.

For pressure loading the pressure wave contours are depicted in Figure 17 at (a) 0.02005 s; (b) 0.02025 s; and (c) 0.0203 s. Also from these graphs the different responses can be appreciated. Note that the chosen time instants differ from those of [25] and Figure 17 is much better resolved. Outgoing waves parallel to the crack can be noticed in this case.




§ DYNAMIC SOLUTION OF A DEBONDING BEAM ON ELASTIC FOUNDATION

Fineberg et al. [46] experimentally studied instability of dynamic fracture on plastic polymethylmethacrylate samples (PMMA). These experiments reveal strong velocity oscillations. Tvergaard and Needleman [47] investigated numerically dynamic crack growth in a porous ductile material and found strong oscillations of J versus time at and after the fracture initiation. Avalanche behavior has been found in fracturing heterogeneous dry materials in quasi static conditions also [3, 27]. We investigate here dynamic fracturing on a dry material by adapting the model shown in Williams [48] which was used to assess the energy dissipated by a single debonding fiber in the cohesive zone using different cohesive forces laws. The mode is based on the deformation of a Euler-Bernoulli beam. We have extended the original model taking into account the inertia force and this first formulation is a precursor for a full model including also the interaction with the fluid phase. With this model, it is possible to show that the velocity of crack tip advancement is not constant but changes in time due to dynamics.


We consider a fiber of total length L as shown in Figure 18, bonded to a fixed support for a length L-L_o. The black line of Figure 18 represents the configuration of half of a crack of initial length L_o; the remaining bonded part, of length L-L_o is the process zone and the fully bonded zone.

From the length L_o to L-L_o, the behavior is described by the equation of a beam on an elastic foundation (10):


    EJ∂^4v(x,t)/∂x^4+f(v)=q(x,t)-ρ A∂^2v(x,t)/∂t^2

where v(x,t) is the vertical displacement of the beam axis, i.e. the crack mouth opening (Figure 18), and it is a function of the spatial coordinate x and the time t. ρ A is the beam weight per unit of length and represents the inertia of the upper part of the solid body with regard to the fracture. The flexural rigidity EJ of the beam represents the stiffness of the homogeneous medium above the fracture plane. q(x,t) is the far field stress due to the external load applied to the body and in the following it is assumed constant in space (i.e. a uniformly distributed load is considered). Finally, f(v) is the function expressing the cohesive forces. At the beginning of the process an initial fracture of length Lo is already formed, so there the cohesive forces are set to 0 (i.e.f(v)=0 ) in equation (10). This leads to the following form of this equation which holds for the region between the origin x= 0 and x=L_o:


    EJ∂^4v(x,t)/∂x^4=q(x,t)-ρ A∂^2v(x,t)/∂t^2

Among the various formulations for the cohesive traction law proposed by Williams [48], we used the one shown in Figure 19 with three different values of threshold. Moreover, at the interface between the two regions proper continuity conditions should be assured. These are described by the following four equations defined at x=x_interface for each time instant (x_interface = L_o at the beginning):

    ∂^iv^L(x_interface,t)/∂x^i=∂^iv^R(x_interface,t)/∂x^i  ,     i=0,...,3

in which the superscripts L and R indicate the displacement v, respectively on the Left and on the Right side of the point at the interface. The coordinate x_interface is changing in time considering that fracture is a moving boundary problem.

The final set of governing equations is formed by eqs. (10), (11) and (12) with the following boundary conditions defined on the left and right side respectively:


    ∂^2v(0,t)/∂x^2=0,  ∂^3v(0,t)/∂x^3=0   ∀ t
    
    ∂^2v(L,t)/∂x^2=0,  ∂^3v(L,t)/∂x^3=0   ∀ t

which correspond to the free-end conditions. In this first application, we consider free boundary at the right-hand side which can produce wave reflection. In the next developments, we will consider   also a transmitting boundary. Note that here we have a fourth order non-homogeneous differential equation compared to a homogeneous second order one as usual in seismic analysis in a continuum [44].

The applied load is considered constant in space and increasing in time from zero up to a constant value. The external load is distributed for the whole length L of the system and it is directed according to the positive vertical axis, (Figure 18).


tarting from an initial condition corresponding to zero velocity and zero displacement, we obtain the fracture length vs time depicted in Figure 19 for different thresholds in the constitutive law. It appears clearly that, due to the inertia forces, the fracture speed is not constant but alternates periods of fast growth with periods of slow down.

Figure 20 shows the evolution of crack opening during the simulation in three distinct time instants: at the beginning of the process for x_interf=L_o, (A), at an intermediate time station (B) and close to the complete debonding (C).




§ CONCLUSIONS

By combining insight from statistical physics with that from fracture mechanics, both analytical and computational, a new light has been shed on simulation of fracture advancement. We have identified two requirements coming from statistical physics dealing with self-organization of rupture which should be included in the fracture advancement rule [2, 3]. These requirements are:  i) the external drive has a much slower timescale than fracture propagation; and ii) the increment of the external load (drive) is applied only when the internal rearrangement of fracture is over. Actually ii) is a consequence of i), and they constitute a "quasi-static" simulation. We show that the second requirement works also well in dynamics situations. These requirements should hence be taken into account when designing experiments and algorithms if we want to observe clean SOC. Ignoring them gives rise to conceptually defective algorithms yielding unphysical results. But experimentally and numerically speaking they are a choice: one can also load faster and then not be able to discern individual avalanches anymore. Most numerical algorithms and exact solutions do this and yield smooth results both for tip advancement and pressure, when carefully carried out experiments and field observations show the opposite. Finally, we show with a debonding beam on elastic foundation model that for a dry material the velocity of crack tip advancement is not constant but changes in time due to dynamics. In this model, the interaction with a fluid will be introduced. 



§ ACKNOWLEDGEMENTS

The authors thank Professor H.J. Herrmann for most useful comments. The work of B.A. Schrefler, has been carried out with the support of the Technische Universität München - Institute for Advanced Study, funded by the German Excellence Initiative and the TÜV SÜD Foundation. C. Peruzzo and F. Pesavento were supported by the project 734370-BESTOFRAC (“Environmentally best practices and optimisation in hydraulic fracturing for shale gas/oil development”)-H2020-MSCA-RISE-2016.

9
1
Tzschichholz, F., Herrmann, H.J., 1995. 
Simulations of pressure fluctuations and acoustic emission in hydraulic fracturing,
Phys. Rev. E 51 (3), 1961–1970.
2
Herrmann, H.J. 1990.
Fracture patterns and scaling laws
, Physica A 163, 359-372
3
Zapperi, S., Vespignani, A. and Stanley, H.E., 1997. 
Plasticity and avalanche behaviour in microfracturing phenomena.
Nature, 388(6643), p.658. (1997).
4
Alava, M.J., Nukala, P.K. and Zapperi, S., 2008. 
Fracture size effects from disordered lattice models.
International journal of fracture, 154(1-2), p.51. 
5
Salerno, K.M., Craig E. Maloney, C.E., Robbins M.O., 2012. 
Avalanches in Strained Amorphous Solids: Does Inertia Destroy Critical Behavior.
Physical Review Letters 109, 105703, 1-5.
6
Dahmen, K.A., Ben-Zion, Y. and Uhl, J.T., 2011. 
A simple analytic theory for the statistics of avalanches in sheared granular materials.
Nature Physics, 7(7), p.554. 
7
Hidalgo, R.C., Kun, F. and Herrmann, H.J., 2001. 
Bursts in a fiber bundle model with continuous damage.
Physical Review E, 64(6), p.066122.
8
Hidalgo, R.C., Kun, F., Kovács, K. and Pagonabarraga, I., 2009. 
Avalanche dynamics of fiber bundle models.
Physical Review E, 80(5), p.051108. 
9
Kun, F., Zapperi, S. and Herrmann, H.J., 2000. Damage in fiber bundle models.
The European Physical Journal B-Condensed Matter and Complex Systems, 17(2), pp.269-279. 
10
Raischel, F., Kun, F. and Herrmann, H.J., 2008. Continuous damage fiber bundle model for strongly disordered materials.
Physical Review E, 77(4), p.046102. 
11
Phillips, W.J., 1972. Hydraulic fracturing and mineralization.
Journal of the Geological Society 128, 337-359.
12
Nolet, G., 2009. Slabs do not go gently.
Science 324 (5931), 1152-1153.
13
Schwartz, S.Y., Rokosky, J.M., 2007. 
Slow slip events and seismic tremor at circum-pacific subduction zones.
Reviews of Geophysics 45, RG3004.
14
Burlini, L., Di Toro G., Meredith, P., 2009. 
Seismic tremor in subduction zones: Rock physics evidence.
Geophysical Research Letters 36, L08305.
15
Burlini, L., and Di Toro G., 2008. 
Volcanic Symphony in the Lab.
Science 322 (5899), 207-208.
16
Morita, N., Black, A.D., and Guh, G.-F., 1990. 
Theory of Lost Circulation Pressure.
Presented at the SPE Annual Technllical Conference and Exhibition, New Orleans, USA, 23–26 September. SPE-20409-MS. http://dx.doi.org/10.2118/20409-MS.
17
Fuh, G.-F., Morita, N., Boyd, P. A. et al., 1992. 
A New Approach to Preventing Lost Circulation While Drilling. 
Presented at the SPE Annual Technical Conference and Exhibition, Washington, DC, 4–7 October. SPE-24599-MS. http://dx.doi.org/10.2118/24599-MS
18
Okland, D., Gabrielsen, G.K., Gjerde, J., Koen, S., Williams, E.L., 2002. 
The Importance of Extended Leak-Off Test Data for Combatting Lost Circulation.
Society of Petroleum Engineers. doi:10.2118/78219-MS. 
19
Candela, T., Brecht Wassing, B., ter Heege, J., Buijze, L., 2018. 
How earthquakes are induced. 
Science360 (6389), 598-599.
20
West Spanish Peak, Colorado, U.S., Wikipedia.
21
Rosen, J., 2016. 
Crystal clocks. 
Science, 6314 (354), 822-825.
22
Khoei, A.R., Barani, O.R., Mofid, M., 2011. 
Modeling of dynamic cohesive fracture propagation in porous saturated media.
International Journal for Numerical and Analytical Methods in Geomechanics 35, 1160–1184.
23
Kostov, N.M., Ning, J., Gosavi, S.V., Gupta P., Kulkarniand, K., Sanz P.F., 2015. 
Dynamic Hydraulic Fracture Modeling for Wellbore Integrity Prediction in a Porous Medium. SIMULIA Community Conference, 1-9. www.3ds.com/simulia.
24
Ouchi, H., John T. Foster, J.T., Sharma, M.M., 2017. 
Effect of reservoir heterogeneity on the vertical migration of hydraulic Fractures.
J. Petroleum Science and Engineering, 151, 384-408
25
Cao, T.D., Hussain F., Schrefler, B.A., 2018. 
Porous media fracturing dynamics: stepwise crack advancement and fluid pressure oscillations. J. Mech. Phys. Solids 111, 113-133.
26
Lecampion, B., Bunger A.,  Zhang, X., 
Numerical methods for hydraulic fracture propagation: A review of recent trends.
J.  Natural Gas Science and Engineering, 49, 66-83
27
Milanese E., Yilmaz O., Molinari J.F., Schrefler B.A., 2016. 
Avalanches in dry and saturated disordered media at fracture.
Physical Review E 93 (4), 043002. http://dx.doi.org/10.1103/PhysRevE.93.043002.
28
Soliman M.Y., Wigwe, M., Alzahabi, A., Pirayesh E., Stegent, N., 2014. 
Analysis of Fracturing Pressure Data in Heterogeneous Shale Formations.
Hydraulic Fracturing Journal 1(2), 8-12.
29
L'homme, T.P., de Pater, C.J., Helferich, P.H., 2002. 
Experimental study of hydraulic fracture initiation in Colton Sandstone
, SPE/ISRM 78187, SPE/ISRM Rock Mechanics Conference, Irving, Texas.
30
Pizzocolo F., Huyghe, J.M.R.J., Ito, K., 2013. 
Mode I crack propagation in hydrogels is stepwise.
Engineering Fracture Mechanics 97, 72-79.
31
Schrefler, B.A., Secchi, S., Simoni L., 2006. 
On adaptive refinement techniques in multifield problems including cohesive fracture.
Computer Methods in Applied Mechanics and Engineering 195, 444–461.
32
Vahab, M., Akhondzadehb, Sh., Khoei A.R., Khalili, N. 2018. 
An X-FEM investigation of hydro-fracture evolution in naturally-layered domains.
Engineering Fracture Mechanics 191,187-204. 
33
Ding, J., Remij, E.W., Remmers J.J.C., Huyghe, J.M., 2018. 
An investigation of the stepwise Mode I crack propagation in hydrogels.
Journal of Applied Mechanics (in press)
34
Omori, F., 1894. J. Coll. Sci. Imper. Univ. Tokyo, 7, 111.
35
Shaw, B.E., 1993. 
Generalized Omori law for aftershocks from a simple dynamics
, Geophys. Res. Lett. 20, 907-910.
36
Cao, D.T., Milanese, E., Remij, E.W., Rizzato, P., Simoni, L.,  Huyghe, J.M., Hussai, F.,Schrefler B.A., 2017. 
Interaction between crack tip advancement and fluid flow in fracturing saturated porous media
, Mechanics Research Communications, 80, 24-37.
37
Secchi, S., Simoni, .L, Schrefler, B.A., 2007. 
Numerical procedure for discrete fracture propagation in porous materials. 
Int J Numer Anal. Methods Geomech 31:109–370
38
Spence, D. A., Sharp, P. (1985). 
Self-similar solutions for elasto-hydrodynamic cavity flow. 
Proceedings of the Royal Society of London, 400, 289–313. 
39
Emerman SH, Turcotte DL, Spence DA.
Transport of magma and hydrothermal solutions by laminar and turbulent fluid fracture. 
Physics of the Earth and Planetary Interiors 1986; 41:249–259.
40
Chen, B., Cen, S., Barrow, A.R., Owen, D.R.J, Li C, 2018. 
Numerical investigation of the fluid lag during hydraulic fracturing.
Engineering Computations, online first, DOI 10.1108/EC-02-2018-0087
41
Biot M.A., 1941. 
General theory of three-dimensional consolidation.
J. Appl. Phys. 12, 155–164.
42
Biot M.A., 1956. 
General solution of the equation of elasticity and consolidation for a porous material. 
ASME J. Appl. Mech. 23, 91–96. 
43
Lewis R.W., Schrefler B.A., 1998. 
The Finite Element Method in the Static and Dynamic Deformation and Consolidation of Porous Media
(John Wiley and Sons, New York).
44
Zienkiewicz O.C., Chan A.H.C., Pastor M., Paul D.K., Shiomi T., 1990. 
Static and dynamic behaviour of soils: a rational approach for quantitative solutions. I. Fully-saturated problems.
Proc. R. Soc. London. A, 429, 285-309.
45
Zienkiewicz, O.C., 2000. 
The Finite Element Method – Volume 1: The Basis
– 5th edition.
46
Fineberg, J., Gross, S.P., Marder, M., Swinney, H.L., 1991. 
Instability in dynamic fracture.
Phys. Rev. Lett. 67 (4), 457–460.
47
Tvergaard, V., Needleman, A., 1993. 
An analysis of the brittle-ductile transition in dynamic crack growth.
International Journal of Fracture, 59, 53-67.
48
Williams J.G., Hadavinia, H., 2002. 
Analytical solutions for cohesive zone models
, JMPS 50, 809–825.




