March 2, 2021

# Variable-order fracture mechanics and its application to dynamic fracture

### Evolutionary governing equations

We briefly discuss the general strategy leading to the formulation of evolutionary governing equations based on VO Riemann–Liouville (VO-RL) derivatives of constants38,44. Then, we apply these operators to formulate an evolutionary elastodynamic framework suitable for the modeling of dynamic fracture. Some background and discussion of the fundamental properties of the VO-RL operators used in this study are provided in Supplementary Note 1.

A particularly interesting property of fractional-order Riemann-Liouville operators stems out of their behavior when applied to the fixed-order derivative of a constant. It is found that this fractional-order derivative is not equal to zero unless the order converges to an integer. While this is an unexpected and maybe even unsettling property of such operators, at least in view of classical integer order calculus, we will show that this property has extremely valuable implications for modeling physical systems exhibiting highly nonlinear and discontinued behavior. Mathematically, the RL derivative of a constant A0(in mathbbR) to a constant fractional-order α0(in mathbbR^+) defined on the interval ((a,t)in mathbbR) is given as35

$$D_t^alpha _0A_0=fracA_0,(t-a)^-alpha _0Gamma (1-alpha _0)$$

(2)

where Γ() is the Gamma function. Note that, although apparently nonintuitive, this is merely an intrinsic property of the RL operator. The use of this property was originally outlined and extended to variable-order in38,44 where it was applied to the modeling of highly nonlinear mechanisms in dynamical systems. More specifically, the properties offered by the VO-RL derivative of a constant creates a unique opportunity to formulate governing equations in an evolutionary form. In the following, we briefly review these characteristics in order to lay the necessary foundation for the development of the VO elastodynamic formulation.

Consider a function α(t) constructed using a continuous real-valued function κ(t) in the following fashion38

$$alpha (t)=exp (-kappa _0kappa (t))$$

(3)

where the function κ(t) is some function designed to capture the desired physical mechanism of interest and the one producing the order variation. Specific details on the selection of this function in the context of fracture mechanics will be provided when addressing the VO dynamic fracture formulation. We emphasize that, while the characteristic function κ(t) introduced above was defined to be a function of time t, the functional dependence can be extended to include any other dependent or independent variables depending on the specific physical problem. Further, (kappa _0in mathbbR^+) is a scaling factor that allows calibrating the order variation on the scale of the characteristic response of the physical system. A detailed discussion of the procedure to determine the value of κ0 is outlined in Supplementary Note 2 along with an illustrative example. For a given κ0, the limiting behavior of α(t) is

$$alpha (t)to left{beginarrayllinfty &kappa (t),<,0 0&kappa (t),>,0endarrayright.$$

(4)

Now, we can indicate the VO-RL derivative of the constant A0 to the order α(t) on the interval (a, t) as (,_aD_t^alpha (t)f(t)) or, in the interest of a more compact notation, as (D_t^alpha (t)f(t)). Equations (2)–(4) lead to

$$,_aD_t^alpha (t)f(t)equiv D_t^alpha (t)A_0=left{beginarrayllmathopmathrmlimlimits_alpha (t)to infty fracA_0,(t-a)^-alpha (t)Gamma [1-alpha (t)]=0&kappa (t)le 0 mathopmathrmlimlimits_alpha (t)to 0fracA_0,(t-a)^-alpha (t)Gamma [1-alpha (t)]=A_0&kappa (t), >,0endarrayright.$$

(5)

It appears that, when the VO-RL operator is applied to a constant under the conditions in Eq. (3), a discontinuous (switch-like) behavior can be captured simply following a change in sign of the function κ(t). It is exactly this switching behavior that can be exploited to simulate the occurrence of certain nonlinear and discontinuous dynamical properties of mechanical systems. More specifically, consider defining VO operators as part of a governing equation such that its variation can capture changes in the properties of the systems such as, for example, a change in stiffness (e.g., bilinear stiffness) or the occurrence of geometric discontinuities (e.g., dislocations in a lattice or crack in a continuum). In all these cases, the response of the system changes from initially linear to, potentially, highly nonlinear. The onset of either type of nonlinearities or discontinuities results in an implicit reformulation of the underlying system dynamics which can be captured in the order α(t) via the function κ(t) embedded in it. In other terms, the governing equations describing a system can be implicitly reformulated via a change in the order α(t), following a change in the underlying physical mechanisms dominating the response of the system. This characteristic was illustrated to formulate evolutionary equations to model contact dynamics, hysteretic behavior38, and motion of edge-dislocations in lattice structures43. In the present study, we extend this unique behavior of the VO-RL operator to simulate the initiation and propagation of cracks in solids. Such behavior is achieved by proper integration of the VO-RL operators in the elastodynamic formulation.

### VO elastodynamic formulation

The strong form of the governing equation for a solid having a volume Ω (see Supplementary Fig. 4) is given in the well-known form

$$nabla cdot boldsymbolsigma +boldsymbolf=rho ddotboldsymbolu$$

(6)

where σ denotes the stress field, u denotes the displacement field, f denotes the externally applied force, and ρ denotes the density of the solid. The bold-face is used to indicate either vectors or tensors. The above equations of motion are subject to the following boundary (BC) and initial (IC) conditions:

$$,textBC,:left{beginarrayllboldsymbolu(boldsymbolx,t)=overlineboldsymbolu(boldsymbolx,t)&boldsymbolxin partial Omega _u hatboldsymbolncdot boldsymbolsigma =overlineboldsymbolt(boldsymbolx,t)&boldsymbolxin partial Omega _tendarrayright.$$

(7)

$$,textIC,:left{beginarrayllboldsymbolu(boldsymbolx,0)=boldsymbolu_0(boldsymbolx)&boldsymbolxin Omega dotboldsymbolu(boldsymbolx,0)=boldsymbolv_0(boldsymbolx)&boldsymbolxin Omega endarrayright.$$

(8)

where ∂Ωu and ∂Ωt denote the portion of the boundary where essential and natural boundary conditions are applied, respectively. (overline{boldsymbolu}) and (overlineboldsymbolt) denote the externally applied displacement and traction at ∂Ωu and ∂Ωt, respectively. While u0 and v0 denote the displacement and velocity fields at t = 0. The stress developed in the medium upon damage is defined in the following fashion

$$boldsymbolsigma =psi (d)mathbbC:boldsymbolvarepsilon$$

(9)

where (mathbbC) is the classical fourth-order elasticity tensor and ε is the symmetric displacement-gradient strain tensor. ψ(d) is a degradation function of the damage variable d [0, 1] such that ∂ψ/∂d ≤ 0; this latter condition originates from the thermodynamic consideration that the degradation function must lead to a decrease in the elastic energy with an increase in damage size. In addition, the degradation function must be bounded in the interval ψ(d)  [0, 1]. The condition on the lower-bound of the degradation function (that is ψ(d) = 0) guarantees the convexity of the potential energy of the system ((=frac12psi (d)boldsymbolvarepsilon :mathbbC:boldsymbolvarepsilon )); this latter aspect being essential for a well-posed formulation and finite element implementation. In this study, the damage variable d is defined such that d = 0 indicates the undamaged state, while d = 1 indicates a fully damaged state.

In the VO dynamic fracture formulation, we adopt a strain-based criterion to detect the onset of damage. More specifically, damage at a given point occurs when the maximum principal strain at the given point exceeds a critical strain derived from the elastic strength of the material. The VO-RL formalism presented previously allows us to define the characteristic function κ(x, t) which allows detecting the onset of damage following the strain-driven physical law. More specifically, we define a VO α(x, t) in the following manner:

$$alpha (boldsymbolx,t)=exp left[-kappa _0underbraceleft(fracoverlinevarepsilon (boldsymbolx,t)-varepsilon _uvarepsilon _uright)_kappa (boldsymbolx,t)right]$$

(10)

where κ0 is the previously introduced scaling factor. (varepsilon _uin mathbbR^+) is the material parameter defining the ultimate tensile strain limit governing the onset of damage. (overlinevarepsilon (boldsymbolx,t)) is the maximum principal strain that occurs at a given point x until the instant t. More specifically

$$overlinevarepsilon (boldsymbolx,t)=mathopmax limits_tau in [0,t]tildevarepsilon (boldsymbolx,tau )=mathopmax limits_tau in [0,t]left[max lefttildevarepsilon _x(boldsymbolx,tau ),tildevarepsilon _y(boldsymbolx,tau ),tildevarepsilon _z(boldsymbolx,tau )rightright]$$

(11)

where (tildevarepsilon (boldsymbolx,tau )) denotes the maximum principal strain component (i.e., maximum of different eigen strain values (tildevarepsilon _x(boldsymbolx,tau ),tildevarepsilon _y(boldsymbolx,tau ),tildevarepsilon _z(boldsymbolx,tau ))) at the point x and at a given time instant τ. Recall that a change in the sign of the argument κ(x, t) within the exponential of the VO results in a reformulation of the underlying governing equations. Exploiting the previously described property of VO-RL operators and defining a physics-driven variation of the order according to Eq. (10), the damage variable can be written as

$$d(boldsymbolx,t)=underbraceD_t^alpha (boldsymbolx,t)d_0_mathrmTerm mathrmI-underbraceD_t^alpha (boldsymbolx,t)varepsilon _uleft[frac1overlinevarepsilon (boldsymbolx,t)exp left(-frac(overlinevarepsilon (boldsymbolx,t)-varepsilon _u)varepsilon _Rright)right]_mathrmTerm mathrmII$$

(12)

where d0 = 1 indicates the maximum possible damage. Before discussing the specific role of the two terms in Eq. (12), we explain the different parameters introduced in the equation. εR is defined as

$$varepsilon _R=2varepsilon _uleft(1-fracl_fl_tright)$$

(13)

where (l_t=2EG_f/sigma _u^2) is the characteristic material length for an isotropic solid having Young’s modulus E, fracture energy Gf, and elastic strength σu47. lf determines a characteristic physical dimension of the area within which the crack is localized and, in numerical implementations, it is directly related to the size of the elements used for the spatial discretization of the domain. In other terms, lf dictates the width of the crack path at a given point, that is the distance perpendicular to the crack path at the same point, within which the damage varies between its extreme values. Further, the parameter εR governs the damage evolution rate that determines the level of damage via term II (see also, Supplementary Fig. 4). In order to guarantee the insensitivity of the results to the specific choice of the numerical mesh adopted, it is necessary that lf < lt6,47. The latter condition also follows from the fact that the size of the elements used to simulate the crack must be smaller than the characteristic material length for accurate resolution of the crack path. Additionally, the condition lf < lt, ensures that the damage variable is bounded, i.e., d [0, 1], which is essential to guarantee the boundedness of the degradation function. Recall that the boundedness of the degradation function guarantees the convexity of the potential energy of the system.

It follows from Eqs. (10) and (12) that d(x, t) = 0 for (overlinevarepsilon (boldsymbolx,t)le varepsilon _u) and d(x, t) → 1 when (overlinevarepsilon (boldsymbolx,t)gg varepsilon _u). Thus, it appears that, when the maximum principal strain (overlinevarepsilon (boldsymbolx,t)) exceeds the critical strain limit εu, the damage is initiated at that particular point. The specific value of the damage variable is determined by the combined effects of the two terms in Eq. (12). While term I in Eq. (12) sets the value of the maximum damage, term II allows for an exponential interpolation of the damage between its extreme values (0 and 1) depending on the amount by which the maximum principal strain ((overlinevarepsilon )) exceeds the critical strain (εu) by using the parameter ϵR. Note that the evolution of both these terms is guided by the VO-RL derivatives. More specifically, the VO-RL operator allows detecting the onset of damage in the solid driven by the VO α(x, t). This leads to an automatic reformulation of the underlying governing equations via Eqs. (6) and (12) in order to account for the occurrence and evolution of damage. Remarkably, the resulting evolutionary VO model does not require any front tracking algorithm nor additional criteria to capture the characteristic features of the dynamic crack mechanism such as roughening and branching. This latter comment will be more evident when contrasting the above formulation with the numerical results presented below.

Remarkably, the VO description of the damage field presented in Eq. (12) leads to a direct interpretation of the approach from the perspective of a CZM with general softening laws. This latter interpretation also allows establishing the length-scale (lf) insensitivity of the VO formulation. In the following, we provide a brief discussion on these aspects and refer the reader to Supplementary Note 3 for further details. The VO approach allows an accurate approximation of the different softening laws in terms of the damage variable. By leveraging Eq. (12), an analytical expression for the ratio of the displacement at a certain level of damage with the displacement at the critical case (d = 1) can be written in terms of the damage variable and other material parameters. The resulting ratio can be substituted within analytically or experimentally identified softening laws. Note that the above procedure to derive the softening laws in the VO approach differentiates itself markedly from phase-field approaches where it is typically not possible to obtain a closed-form expression of the displacement ratio in terms of the damage variable. To overcome this latter aspect, a specific form of the softening law must be assumed a priori in phase-field approaches29,30,32. On the contrary, in the VO formulation, rational softening laws arise naturally without any a priori assumption. Rational softening laws have been connected to the ability to suppress mesh-bias in numerical implementations48, hence the fact that they naturally emerge in the VO formulation highlights another key feature of the methodology. The following two softening laws are used in this study

$$psi (d)=(1-d)left[frac1+C_0+0.5C_0^21+C_0+0.5C_0^2-dright]$$

(14)

$$psi (d)=left[left[1+left(fraceta _1d(C_0+0.5C_0^2)1+C_0+0.5C_0^2-dright)^3right]exp left[fraceta _2(C_0+0.5C_0^2)1+C_0+0.5C_0^2-dright]-frac(1+eta _1^3)(C_0+0.5C_0^2)1+C_0+0.5C_0^2-dexp (-eta _2)right]$$

(15)

where C0 = εu/εR, η1 = 3 and η2 = 6.93. Equation (14) is the VO approximation of the linear law for brittle materials, while Eq. (15) is the VO approximation of Cornelissen’s law49 for concrete (a quasi-brittle material). Note that, in both cases, ψ(d) is a rational function of d and is bounded so that ψ(d)  [0, 1] for d [0, 1]. Also, C0 = εu/εR = 2(1 − lf/lt) and thus, in the limit lf → 0 or lflt (that is when the crack is a sharp discontinuity), the softening laws are independent of lf. We have shown in Supplementary Note 3 that both the analytical softening laws are reproduced very well by the approximated VO laws in Eqs. (14) and (15) for lflt. This result confirms that the proposed fractional-order approach is essentially equivalent to the CZM in a 1D setting. In Supplementary Note 3 we show that, for a well-posed model satisfying the thermodynamic constraint lflt, the displacement jump (and consequently, both strain and stress) do not depend on lf. These characteristics indicate that the proposed fractional-order approach is insensitive to the length-scale parameter lf subject to the thermodynamic constraint lflt, at least in the 1D scenario. For more complex settings, the lf insensitivity was established by comparing the temporal variations of quantities such as the elastic strain energy, the crack length, and the crack-tip velocity, for different values of lf (see Supplementary Note 3).

The VO dynamic fracture formalism deserves some additional remarks. First, note that the constitutive relations defined in Eq. (9) result in an identical tensile and compressive fracture behavior which is not generally true when modeling the failure of brittle and quasi-brittle solids. Several researchers have captured the asymmetric tensile/compressive damage by performing a spectral decomposition of the strain energy density and by degrading only the positive strain energy26,32. In this study, we incorporate this asymmetric behavior via the maximum principal strain-based damage criterion, which follows from the well known Rankine criterion30. In other terms, as described previously, the crack is allowed to nucleate only when the maximum principal strain exceeds the critical tensile strain (εu). This specific feature ensured via the VO defined in Eq. (10) allows modeling the asymmetric damage behavior. Further, the definition of the parameter (overlinevarepsilon ) in Eq. (11) based on its past history, along with the condition (dotd(boldsymbolx,t)ge 0), ensure irreversibility of the system. More specifically, these conditions ensure that the length of the crack, denoted by Γc, is monotonically increasing, that is Γc(t1)  Γc(t2) t1 < t2. In addition, the use of the strain-history-based parameter (overlinevarepsilon ) leads to simpler numerical implementations as it allows for an operator split algorithm within a given time-step, wherein the displacement field and the damage field are updated in a staggered manner. The same concept, albeit using a strain energy-based history variable, is often employed in phase-field models of dynamic fracture26,28. Following this staggering numerical implementation, the computation of the damage field is a purely algebraic operation and it does not require the minimization of an additional potential function which, in the case of phase-field models, corresponds to the crack surface density function. The most immediate consequence is that the VO approach reduces the computational cost of dynamic fracture when compared to phase-field models. Finally, we highlight that the VO elastodynamic formulation described above does not account for contact conditions, such as those that occur when the free surfaces of a crack come in contact when subjected to compressive loads. Note that this is not a limitation of the methodology but merely a decision of the authors to focus this work on aspects concerning crack initiation and propagation. Indeed, the contact problem is typically not addressed in classical treatments of dynamic fracture. However, the VO formulation can easily account for contact dynamics by simply adding dedicated terms in the VO derivative. The case of contact via VO operators was previously treated by the authors in44, albeit only for discrete systems.