Variable-Order Fracture Mechanics and its Application to Dynamic Fracture

This study presents the formulation, the numerical solution, and the validation of a theoretical framework based on the concept of variable-order mechanics and capable of modeling dynamic fracture in brittle and quasi-brittle solids. More specifically, the reformulation of the elastodynamic problem via variable and fractional order operators enables a unique and extremely powerful approach to model nucleation and propagation of cracks in solids under dynamic loading. The resulting dynamic fracture formulation is fully evolutionary hence enabling the analysis of complex crack patterns without requiring any a prior assumptions on the damage location and the growth path, as well as the use of any algorithm to track the evolving crack surface. The evolutionary nature of the variable-order formalism also prevents the need for additional partial differential equations to predict the damage field, hence suggesting a conspicuous reduction in the computational cost. Remarkably, the variable order formulation is naturally capable of capturing extremely detailed features characteristic of dynamic crack propagation such as crack surface roughening, single and multiple branching. The accuracy and robustness of the proposed variable-order formulation is validated by comparing the results of direct numerical simulations with experimental data of typical benchmark problems available in the literature.


Introduction
Fracture is one of the most commonly encountered mode of failure in structural systems across a broad spectrum of applications spanning the civil, mechanical, and aerospace engineering fields. The prevention of fracture-induced failure is a major concern of structural design and has historically motivated the development of theoretical and experimental methodologies to predict nucleation and propagation of structural damage. While the general topic of fracture mechanics is very complex in itself due to the coexistence of multiple physical processes occurring over multiple spatial scales, the specific topic of dynamic fracture is possibly even more challenging due to the occurrence of crack surface roughening, instabilities, and branching. Detailed discussions on the implications and modeling approaches for dynamic fracture can be found in many sources such as [1]. During the last few decades, the analysis of dynamic fracture has certainly largely benefited and made significant progress thanks to the rapid development of numerical methods. From a high level perspective, the approaches available for the analysis of damage can be divided into two categories, namely, discrete and continuum. This classification refers specifically to the modeling of the damage interface so that, while in both cases the solid is treated as a continuum, in the former class of approaches the displacement is modeled as a discontinuous field across the fracture surface. In the latter category, instead, the displacement is treated as continuous field everywhere (even across the crack surface), but the local value of the elastic energy is reduced by accounting for the softening of the material properties associated with the fracture-induced degradation. In the following, we briefly review some of the most accredited dynamic fracture models in order to clearly define the context in which our variable-order approach is defined.
Discrete approaches to the modeling of dynamic fracture include extended finite element methods (XFEM) [2][3][4][5], discontinuous cell methods (DCM) [6], cohesive interface element techniques [7][8][9], discontinous Galerkin methods [10,11], and lattice based models [12,13]. From a general perspective, these approaches are based either on linear elastic fracture mechanics (LEFM) [14] or on the cohesive zone model (CZM) [15]. Owing to its computational and multiscale analysis capabilities, the XFEM has quickly risen in popularity and it is currently one of the most widely used approaches. In XFEM, cracks are represented as discrete discontinuities that are embedded in the damaged elements by enriching the displacement field according to the method of partition of unity [16]. This approach implies that the front of the discontinuity (i.e. the crack) must be tracked explicitly. While several tracking algorithms have been proposed over time [4,5,17], the front tracking process is quite computationally intensive, particularly in three-dimensional problems involving complex crack topology. Another important limitation consists in the need for a branching criterion which is often ad-hoc and limited to two crack branches. Exceptions to this latter comment are the formulations based on either the DCM [6] or the interface elements [7,18], which however must be inserted in the model a priori hence posing the problem of knowing the location and path of propagation of damage. The front tracking limitation was addressed by the use of lattice models where the continuum is replaced by a system of rigid particles that interact via a network of linear and nonlinear springs. More recently, Silling [19,20] proposed an approach denominated peridynamics that models the solid medium as a nonlocal lattice of particles described via an integral formulation. During the last two decades this approach has received much attention and it has been used in many diverse applications. In the context of dynamic fracture, peridynamics has shown to be able to address several of the above shortcomings and could accurately capture crack intersections and branching in complex structures and materials. An important caveat of the lattice models derives from the fact that the springs stiffness is often defined heuristically and various elastic phenomena (e.g. the Poissons effect) cannot be reproduced exactly.
In the second category, the continuum approaches, we find the crack band method [21], nonlocal integral damage models [22], nonlocal stress based damage models [23], and the more recent class of phase-field models [24][25][26][27][28][29][30]. Phase-field models are undoubtedly the methods that have seen the highest popularity given their overall accuracy and ease of implementation. In phase-field models, sharp cracks are regularized by a diffused damage field while a variational approach is adopted to obtain evolutionary equations for both the displacement and the damage fields [24]. The formulation also includes a small and positive length scale parameter so that, in the limit for the parameter approaching zero, the phase-field representation of the crack converges to the original problem of a sharp crack [31]. The use of a phase-field regularization prevents the need for an explicit tracking of the crack surface discontinuity. It follows that the numerical implementation of the phase-field model is relatively straight-forward when compared to the previously mentioned discrete approaches. An important disadvantage of these models lies in their high computational cost, which follows from the need to solve a coupled system of partial differential equations for both the damage (phase) and the displacement fields [30]. This limitation becomes even more significant when the phase field approach is applied for fracture analysis in three-dimensional media. Additionally, phase-field models are subject to an artificial widening effect in the damaged area at the point of occurrence of instability [30,32], which is in contrast to the microbranching and crack surface roughening seen in experiments [33]. A detailed review on phase-field models can be found in [34].
In very recent years, variable-order fractional calculus (VO-FC) has emerged as a powerful mathematical tool to model a variety of discontinuities and nonlinear phenomena. Variable-order (VO) fractional operators are a natural extension of constant-order fractional operators that allow the differentiation and integration of functions to any real or complex valued order [35]. In VO operators, the order can be function of time, space, internal variables (e.g. system energy or stress) or even of external variables (e.g. temperature or external mechanical loads) [36]. As the VO-FC formalism allows updating the system's order depending on either its instantaneous or historical response, the corresponding model can evolve seamlessly to describe widely dissimilar dynamics without the need to modify the structure of the underlying governing equation. Thus, a very significant feature of VO-based physical models consists in their evolutionary nature; a property that can play a critical role in the simulation of nonlinear dynamical systems [36][37][38]. In recent years, many applications of VO-FC to practical real-world problems have been explored including, but not limited to, modeling of anomalous diffusion in complex structures with spatially and temporally varying properties [39], the response of nonlinear oscillators with spatially-varying constitutive law for damping [37], nonlinear dynamics with contacts and hysteresis [38], and even complex control systems [40]. The interested reader can find a comprehensive review of the applications of VO-FC in [41].
In this study, we present a theoretical and computational framework based on VO fractional operators and capable of effectively capturing the many features of dynamic fracture in brittle and quasi-brittle solids. We will show how the many unique capabilities of this framework build directly on the several remarkable properties afforded by these fascinating mathematical operators. Dynamic fracture is a quintessential evolutionary nonlinear dynamic problem that involves the propagation of nonlinearities and discontinuities through a system. VO fractional operators are uniquely equipped to model these complex class of dynamical problems. The VO framework presented in this paper builds upon the mathematical structure presented in [42] which focused on the modeling of the propagation of dislocations through lattices of particles using physics-informed order variations. More specifically, the VO model introduced in [42] leveraged an order variation law based on the relative displacements of particles within the lattice in order to capture the formation and annihilation of pair-wise bonds. The general strategy followed that outlined in [38,43] for physics-driven VO laws for discrete systems. The approach resulted in the formulation of evolutionary VO fractional differential equations capable of capturing the transition towards a nonlinear dynamic regime (associated with the motion of dislocations) without having to explicitly track the location of the dislocation. In this study, we extend this general approach to continuous systems by formulating a VO elastodynamic framework uniquely suited for the analysis of dynamic fracture and capable of detecting the formation and propagation of damage by means of a strain-driven order variation laws. The introduction of VO operators in the continuum elastodynamic formulation allows the governing equations to evolve (from linear to nonlinear) and adapt (by capturing discontinuities) based on both the local response and the underlying damage mechanism while eliminating the need for explicitly tracking the damage front. We will show that the resulting formulation is capable of capturing key features associated with the dynamic fracture mechanism such as roughening of the crack surface, crack instability, and crack branching without the need of any a-priori assumptions or ad-hoc criteria. Further, contrary to phase-field models, no additional partial differential equations are needed to predict the evolution of the damage field. Indeed, in the VO framework, the damage field evolves naturally guided by the variation of the order of the fractional operators that solely depends on the instantaneous response of the system. In the second part of this study, the VO dynamic fracture model is validated by applying it to the direct numerical simulation of three benchmark experiments available in the literature: 1) the Kalthoff-Winkler experiment that involves the impact shear loading of a doubly notched specimen [44], 2) the dynamic crack branching experiment [33], and 3) the John-Shah experiment that involves the impact loading of a pre-notched concrete slab [45].

Material and Methods
We briefly discuss the general strategy leading to the formulation of evolutionary governing equations based on VO Riemann-Liouville (VO-RL) derivatives of constants [38,43]. 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 Information (SI).
Evolutionary governing equations via VO-RL derivatives of constants: 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 discontinues behavior. Mathematically, the RL derivative of a constant A 0 ∈ R to a constant fractional-order α 0 ∈ R + defined on the interval (a, t) ∈ R is given as [35]: where Γ(·) is the Gamma function. Note that, although apparently non intuitive, this is merely an intrinsic property of the RL operator. The use of this property was originally outlined and extended to variable-order in [38,43] where it was applied to the modeling of highly nonlinear mechanisms in dynamical systems. More specifically, the properties offered by the variable-order Riemann-Liouville (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 fashion [38]: 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, κ 0 ∈ R + 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 the SI along with an illustrative example. For a given κ 0 , the limiting behavior of α(t) is: Now, we can indicate the VO-RL derivative of the constant A 0 to the order α(t) on the interval (a, t) as a D α(t) t f (t) or, in the interest of a more compact notation, as D It appears that, when the VO-RL operator is applied to a constant under the conditions in Eq. (2), 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 types of nonlinearities or discontinuities results in an implicit reformulation of the underlying system dynamics. This change in the underlying dynamics can be captured in the order α(t) via the function κ(t). It immediately follows that a change in the order α(t) results in an implicit reformulation of the equations of motion 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 behavior [38], and motion of edge-dislocations in lattice structures [42]. 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 Fig. (S3) in SI) is given in the well-known form: 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: where ∂Ω u and ∂Ω t denote the portion of the boundary where essential and natural boundary conditions are applied, respectively. u and t denote the externally applied displacement and traction at ∂Ω u and ∂Ω t , respectively. While u 0 and v 0 denote the displacement and velocity fields at t = 0. The stress developed in the medium upon damage is defined in the following fashion: where C 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 this study, the damage variable d is defined such that d = 0 indicates the undamaged state, while d = 1 indicates a fully damaged state. We note that the same stress-damage-strain constitutive relation has also been adopted in several classical dynamic fracture formulations [6,29,30].
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: where κ 0 is the previously introduced scaling factor. ε u ∈ R + is the material parameter defining the ultimate tensile strain limit governing the onset of damage. ε(x, t) is the maximum principal strain that occurs at a given point x until the instant t. More specifically: whereε(x, τ ) denotes the maximum principal strain component (i.e., maximum of different eigen strain values ε x (x, τ ),ε y (x, τ ),ε z (x, τ )) 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. (8), the damage variable can be written as: where d 0 = 1 indicates the maximum possible damage. Before discussing the specific role of the two terms in Eq. (10), we explain the different parameters introduced in the equation. ε R is defined as: where l t = 2EG f /σ 2 u is the characteristic material length for an isotropic solid having Young's modulus E, fracture energy G f , and elastic strength σ u [46]. l f 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, l f 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, Fig. (S3) of SI). In order to guarantee the insensitivity of the results to the specific choice of the numerical mesh adopted, it is necessary that l f < l t [6,46]. 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.
It follows from Eqs. (8,10) that d(x, t) = 0 for ε(x, t) ≤ ε u and d(x, t) → 1 when ε(x, t) ε u . Thus, it appears that, when the maximum principal strain ε(x, t) exceeds the critical strain limit ε u , 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. (10). While term I in Eq. (10) 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 (ε) 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. (5,10) 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.
The VO dynamic fracture formalism deserves some additional remarks. First, note that the constitutive relations defined in Eq. (7) result in an identical tensile and compressive fracture behaviour 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 energy [26,32]. In this study, we incorporate this asymmetric behavior via the maximum principal strain based damage criterion, which follows from the well known Rankine criterion [30]. 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. (8) allows to model the asymmetric damage behaviour. Further, the definition of the parameter ε in Eq. (9) based on its past history, along with the conditionḋ(x, t) ≥ 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 (t 1 ) ⊆ Γ c (t 2 ) ∀ t 1 < t 2 . Additionally, the use of the strain-history based parameter ε 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 fracture [26,28]. Following this staggered 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, note that the damage value d = 0 obtained for ε = ε u indicates a zero crack-tip opening displacement. As ε increases, the damage increases leading to an increase in the crack-tip displacement. Since the crack-tip opening displacement is directly linked to the strain value, it follows that the ratio of the instantaneous crack-tip opening displacement to the maximum crack-tip opening displacement is equal to the ratio of the instantaneous strain to the maximum strain (obtained for d = 1). Further, given the strainbased definition of the damage variable in Eq. (10), it follows that the ratio of the crack-tip displacement (δ c ) to the critical crack-tip displacement (δ 0 ) can be given as the ratio of the damage parameters d/d 0 . This formulation allows expressing different traction separation laws (also called softening laws) directly in terms of the damage variable d without the need for additional characteristic functions similar to [6]. In this study, we use two different traction separation laws expressed in terms of the damage variable as: Equation (12a) is a linear law for brittle materials [29], while Eq. (12b) is the Cornelissen law [47] for concrete (a quasi-brittle material). This latter law was obtained experimentally in [47] and the coefficients were found to be η 1 = 3 and η 2 = 6.93. Note that, in both cases, ψ(d) is bounded so that ψ(d) ∈ [0, 1] for d ∈ [0, 1]. 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 in [43], albeit only for discrete systems.

Results
To demonstrate the accuracy and the robustness of the VO fracture mechanics framework, we apply it to perform numerical simulations of three classical benchmark experiments available in the literature. The first two examples refer to the Kalthoff-Winkler experiment [44] and the classical crack-branching experiment [33]; both pertain to the dynamic fracture of brittle solids. The last example consists of mixed-mode fracture of a concrete slab under an impact load that was analyzed experimentally in [45]. In all three examples, the numerical results were obtained using a plane strain elastic model. The computational domain was discretized using uniform quadrilateral elements and the dynamic solution was computed using an explicit Newmark solver. Further, a lumped mass matrix was used in the dynamic solver to suppress high-frequency elastic oscillations (noises) and to ensure conservation of energy. The time-step (∆t) used in the dynamic solver, was determined using the Courant-Friedrichs-Lewy (CFL) condition. For a more conservative scheme, we used ∆t = 0.9∆t 0 = 0.9h/c, where h denotes the size of an element within the mesh and c denotes the speed of compressional waves in the medium under consideration. Further details on the numerical implementation are provided in SI. Further, videos of the growing crack front in the three benchmark cases are provided as multimedia supplementary information.
Kalthoff-Winkler experiment: The classical Kalthoff-Winkler experiment [44] consists of an unrestrained doubly notched specimen subject to an impact load, as illustrated in Fig. (1a). Following the original experimental setup [44], the specimen was made of maraging steel with the following material properties: E = 190 GPa, σ u = 844 MPa, G f = 22.2 N/mm, ν = 0.3, and ρ = 8000 kg/m 3 [6]. The characteristic material length corresponding to the aforementioned material properties is obtained as l t = 0.012 m. It was observed experimentally that, for lower strain rates (v 0 = 16.5 m/s), brittle failure occurs and the cracks nucleate from the edges of both notches at an angle of about 70 • with respect to the horizontal axis (which coincides with the line of symmetry). In the following, we numerically analyze this benchmark problem using the VO dynamic fracture model.
Given the symmetry of both the specimen and the test conditions in the original experiment, we modeled only the upper-half of the specimen in order to reduce the computational cost. The vertical component of the displacement field was set as zero (u y =0) at the line of symmetry indicated in Fig. (1a) to impose the symmetric boundary conditions. To model the impact load, a velocity boundary condition was applied at the nodes corresponding to the impact zone and the impulse was kept constant throughout the dynamic simulation. Further, the linear softening law (see Eq. (12a)) was used to model the degradation in the elastic energy upon damage development. The damage pattern generated using the VO model is presented in Fig. (1b,1c) for two different mesh configurations. The element size is taken as h = 0.5 mm in Fig. (1b), and h = 0.25 mm in Fig. (1c). As evident from Figs. (1b,1c), although a sharper crack is obtained in the case of the finer mesh, the overall crack propagation features are insensitive to the specific mesh configuration. For both mesh configurations, the crack develops in a direction that forms approximately 70 • with the horizontal axis. The average angle from the initial crack tip to the point where the crack intersects the top boundary is obtained as 72 • which matches well with the experimental result in [44]. Further, the crack intersects the top boundary of the specimen at a time instant of 75 µs for both the mesh configurations; this time corresponds to an average crack propagation velocity of c = 1064 m/s. Note that c < 0.6c R , where c R (= 2745 m/s) denotes the Rayleigh wave speed in the medium, as observed commonly in experiments [33]. For completeness, we highlight that unlike the results presented in [5,6] and obtained using different FEM and DCM models, we did not observe any spurious cracks developing from the bottom right corner and travelling towards the tip of the notch.
Dynamic crack branching: In this benchmark problem, we model a pre-cracked specimen loaded dynamically in tension as illustrated in Fig. (2a). This problem has been widely adopted in the literature to study dynamic crack branching, both experimentally [33] and numerically [4][5][6]30,32]. The specific material parameters used in this simulation were E = 32 GPa, σ u = 3.1 MPa, G f = 3 J/m 2 , ν = 0.2 and ρ = 2450 kg/m 3 [6], which correspond to a glass type material. The associated characteristic length is found to be l t = 0.02 m. Further, given the brittle nature of the material in this experiment, the linear softening law was used to account for the degradation in the elastic strain energy upon damage. A uniform and constant traction of magnitude σ 0 = 1 MPa is applied instantaneously to the rectangular specimen on its top and bottom surfaces at the initial time step and is held constant throughout the simulation. All other surfaces of the specimen are left free. This load condition is such that crack branching occurs in the specimen. The simulations obtained via the VO formulation are presented in Fig. (2b-2d). Figure (  The results obtained via the VO model lead to the following remarks. First, as discussed in [33], upon the onset of instability several microcracks develop from the principal propagating crack branch and interact with one another simultaneously. This process ultimately leads to a roughening of the crack surface. As evident from the inset in Fig. (2b), the VO formulation is able to capture in great detail the roughening of the crack surface due to the emerging microbranches. Similar to the experiments conducted in [33], these microbranches vary in size and the larger ones develop into full fledged branches. The remaining microbranches are arrested as a result of dynamic interaction with the growing ones. This set of results also highlights some of the advantages of the VO model over phase-field models which invariably capture an artificial widening effect in the damaged area at the point of occurrence of instability [30,32]. Also extremely remarkable, unlike classical dynamic fracture models that capture two branches [5,6,30,32], the VO model predicts four branches nucleating from the point of instability. This result closely matches the experimental results in [33], where it was demonstrated that the number of branches developed can vary between two and four.Further, we emphasize that unlike classical discrete approaches to dynamic fracture [2,3,5,8] we did not impose any additional criteria within the VO model to facilitate the crack branching behavior; the branching and roughening occurs naturally as a result of the local response field. Similar to phase field (variational) approaches, the VO dynamic fracture model leads to full crack identification without the support of additional branching conditions. Finally, as shown in [33], the number of branches exceeds four when also considering unsuccessful (i.e. not fully developed) branches. This feature is also captured by the VO model wherein we see that a number of unsuccessful branches nucleate from the principal branches.
John-Shah experiment: Another benchmark problem to test the performance of the VO fracture mechanics framework involves the three-point bending of concrete beams subject to impact loading [45]. The geometry and boundary conditions for the specimen involved in this test are illustrated in Fig. (3a). In this classical benchmark problem, a pre-built notch (offset from the mid-span axis of symmetry) is used to study mixed-mode fracture in concrete beams. It was observed by John and Shah that the parameter γ = l 1 /l 2 (see Fig. (3a)), that controls the placement of the notch, plays a critical role in determining what failure mode and damage pattern would occur in the specimen after the impact. Indeed, there exits a critical value γ c such that, for γ < γ c , the crack nucleates from the notch tip while, for γ > γ c , the crack nucleates from the mid-span. In addition, there exists an intermediate value of γ close to and less than γ c wherein both cracks develop. The experimentally determined value of γ c was γ c = 0.77 [45].
We simulated this benchmark problem using the VO framework. The specific material parameters used in the simulation were E = 34 GPa, σ u = 1 MPa, G f = 31.1 J/m 2 , ν = 0.2 and ρ = 2400 kg/m 3 . The degradation in the strain energy upon damage was modeled using the Cornelissen softening law for concrete (see Eq. (12b)). The impact velocity is given by the linear ramp [17]: where v 0 = 0.06 m/s and t 0 = 196 µs. Using the above material properties, geometry, and loading conditions we simulated the dynamic three-point bending for three different values of γ ∈ {0.72, 0.76, 0.79} corresponding to three different notch locations. In all the three cases, the computational domain was uniformly discretized using elements of size h = 0.635 mm. Note that the characteristic material length corresponding to the material properties of the concrete specimen is obtained as l t = 2.1 m. We merely observe that for quasi-brittle materials like concrete, the characteristic length-scale is generally too large when compared to the dimensions of laboratory specimens and hence, the condition l f < l t is virtually meaningless for quasi-brittle materials. The crack obtained for the three different cases are presented in Figs. (3b-3d).
Overall, the results of the three numerical experiments compare very well with the experimental results. In particular, as in the experimental results, the crack nucleates from the tip of the notch for γ = 0.72 < γ c and from the mid-span of the beam for γ = 0.79 > γ c , and propagates towards the top surface of the beam. Further, similar to the experiment conducted in [45], a transition state is observed for the case γ = 0.76 wherein cracks propagate from both the notch tip and mid-span towards the top surface. It follows that the estimate of the critical notch location γ c obtained via the VO dynamic fracture model lies in (0.76, 0.79), which is in good agreement with the experimental value of γ c = 0.77.

Conclusions
This study presented a novel elastodynamic formulation based on variable order fractional operators and capable to provide accurate estimates of dynamic fracture in brittle and quasi-brittle solids. From a mathematical perspective, the peculiar properties of the VO-RL operator enable capturing the behavior of highly nonlinear systems with evolving discontinuities, such as those involved in the nucleation and propagation of cracks in solids. We showed that an apparently unsettling property of the RL operator, that is the non-vanishing derivative of a constant, can have very useful implications to model dynamic fracture. Furthermore, the ability of VO operators to update their order as a function of either dependent or independent variables, results in governing equations that can evolve in real time to capture growing cracks without requiring modifications to the fundamental governing equations. Even more remarkable, and certainly in stark contrast with more traditional approaches to dynamic fracture, is the fact that VO governing equations do not require a priori assumptions or additional conditions to detect characteristic aspects of dynamic fracture such as crack nucleation, crack surface roughening, crack instability and branching. In other terms, the nonlinear and discontinuous dynamic behavior associated with fracture naturally emerges based on the instantaneous response of the system. Further, given the many recent advances in the formulation of fractional order mechanics as a comprehensive approach to nonlocal elasticity, it can be envisioned that the present VO elastodynamic framework could be easily integrated in a fully fractional formulation hence leading to a powerful tool for dynamic fracture analysis of nonlocal media.