The study of diffuse interface propagation of dynamic failure in advanced ceramics using the phase-field approach

Dynamic failure in advanced ceramic materials is a complex mechanical phenomenon, which is challenging to study in experimental and computational works. This work models and investigates the crack initiation, growth, propagation, and branching in two-dimensional single crystalline boron carbide at the nanometer length-scale and picosecond time-scale. A phase-field approach is used by carefully selecting the interfacial energy allowing a monolithic numerical solution method capturing strong coupling between mechanics and damage. Specifically, the transient Ginzburg–Landau equation is coupled with the balance of momentum including geometric nonlinearities and solved by using the open-source computing platform FEniCS. The monolithic computation ensures necessary accuracy for dynamic fracture under pure mode I loading, and also demonstrates the capability of crack branching depending on the loading rate. The effect of crack regularization length and kinetic coefficient on the crack interface profile and crack tip velocity is also studied. Altogether, the results are important for modeling anisotropic fracture in advanced ceramics and designing materials with desired dynamic failure characteristics.


Introduction
Being the most important failure mode in solids and structures, the evolution of cracks as well as their impact on material designs are of great importance to both scientists and engineers equally.According to Cox et al. [1], all fracture is dynamic at some length-scale and time-scale.Even under quasi-static loading, the crack may propagate at a speed comparable to that of the mechanical waves [2].Therefore, the unstable propagation of cracks under external dynamic loading necessitates the consideration of inertial effects to analyze the transient behavior of structures or the interactions between stress waves and cracks [3].To date, numerous experimental, analytical, and numerical studies have been considered to investigate the dynamic failure of brittle materials.Experimentally, the split-Hopkinson pressure bar [4][5][6], ballistic impact [7], plate impact [8], laser spall [9], laser shock [10], and coupling of real-time visualization techniques (e.g., ultra high-speed imaging) with laboratory-based experiments [11] have been employed to probe the fracture behavior of advanced ceramics during dynamic loading; however, the energy flux to the crack tip and its connection to the crack velocity, as one of the main difficulties in the theory of high speed dynamic crack propagation, cannot be directly measured from experimental tests [12].
Atomistic simulations, largely represented by molecular dynamics (MD) simulations and density functional theory (DFT), are the bridge between continuum and atomic description of dynamic fracture.Examples are the method of lattice dynamics modeling [13] and tight binding methods [14].Despite having successful connection between atomistic approaches and continuum-scale theory and experimentation [15], these atomistic viewpoints often suffers from limitation of modeling time, which leads to a loading rate much higher than the rate of practical interest [16], and are forbiddingly expensive if extended to realistic length scales [17].
Among the theoretical studies for understanding the underlying physics and mechanics of dynamic fracture [18,19], the Griffith criterion was one of the earliest description which was based on an energy balance equation [20]; however, practical engineering examples are usually too complex to be solved using analytical techniques.Therefore, numerical methods play a vital role in understanding dynamic failure.For example, discrete crack models, such as the discrete element method [21,22], the extended finite element method (XFEM) [23][24][25], the cohesive zone method [26], and the cohesive segment method [27], allow the displacement field to be discontinuous across the fracture surfaces.Regardless of showing much success in modeling crack propagation [28], the discrete crack models need additional energy-based criteria to simulate dynamic fracture problems [29,30].In addition, the discrete crack model requires remeshing algorithms or using the partition of unity method [31], both having their own difficulties in tracking the multiple crack fronts in complex three dimensional morphologies [32,33].In the smeared (continuum) crack model, including the gradient damage model [34,35] and diffuse interface models [36][37][38][39][40], the displacement is continuous everywhere and stresses gradually decrease to model the degradation process.The coupling of discrete and smeared crack approaches (e.g., the element deletion method [41] and the thick level-set method [42]) have also shown promising results in modeling fracture.However, there are drawbacks as the dependence of the results on the finite element meshes and the convergence path of the solutions result in numerical errors [43].
The phase-field model is an alternative approach that has been successfully adopted in the simulation of martensitic phase transformations [44,45], melting [46,47], dislocations [48,49], twinning [50,51], fracture [52][53][54][55][56][57][58][59][60][61], and their interactions [62][63][64][65].One of the main advantages of the phase-field approach to study the failure mechanisms in brittle materials over the previous methods is predicting the evolution of interfaces (e.g., merging and branching of multiple cracks) with no additional effort [66][67][68], making phase-field modeling a powerful and flexible method for studying the fracture of single-crystalline [69,70], polycrystalline [71,72], and granular materials [73,74].In the physics community, the phase-field models are commonly derived by adapting the phase transition formalism of Landau and Ginzburg [75] to study dynamic crack propagation in brittle solids [76,77].The advantage of the Ginzburg-Landau approach over the incremental energy minimization method [69] is that material parameters associated with time scales for interfacial motion enter the model, making it a more general model for studying the dynamic failure of brittle materials.
Numerically, the coupled nonlinear equilibrium and the Ginzburg-Landau equations are solved via a staggered or monolithic scheme.The former is based on decoupling balance equations and the phase-field problem into the system of two equations that are solved in a subsequent manner [78,79].Giving rise to two convex minimization problems, the method is robust; however, a significant amount of staggered iterations may be required at a fixed loading step, resulting in a higher computational cost [80].In the monolithic scheme, all solution variables are solved simultaneously [81], and this is more efficient for strongly coupled systems of equations as a result of less Newton-Raphson iterations [82].
The current paper extends the Ginzburg-Landau phase-field approach to predict dynamic fracture in single crystal boron carbide (B 4 C) under mode I loading.Here, the focus of this work is on derivation of governing equations and solving them monolithically in order to increase the accuracy for applications with strong coupling between mechanics and damage growth.Governing equations for the phase-field method is motivated differently than balance equations and there is a computational difficulty to implement monolithic schemes.Hence, the literature lacks studies focused on solving the Ginzburg-Landau based phase-field problem for predicting the dynamic crack branching of brittle materials by using a monolithic scheme.By following the works on local stress concentrations in nanoscale defect-free volumes or by high pressures [83], we develop a highly nonlinear phase-field theory for elasticity with considerations for anisotropic surface energy [84].In this way, numerical problems are circumvented owing to governing equations motivated by thermodynamics, and we manage to solve these nonlinear and coupled differential equations by exploiting the open-source parallel computing platform FEniCS [85,86].
The remainder of this paper is outlined as follows: In Section 2, we present continuum mechanics and thermodynamically sound derivation of equations and their variational formulation for the finite element method.Then in Section 3, results and representative material properties are reported along with the discussion of phase-field simulations.The conclusions of the study are drawn in Section 4.

Nonlinear phase-field theory
We describe a model for a single fracture system in solids based on thermodynamical derivations.The present approach accounts for the time-evolution of the fracture order parameter towards an equilibrium state for predicting the crack paths in anisotropic single crystal materials.This allows the study of spatio-temporal fluctuations of damage variables, as well as the nanoscale dynamics that govern various pattern forming phenomena [87,88].Moreover, the interfaces, their propagation, and interactions, which are the most important features governing the formation of microstructures in materials [89][90][91][92], is studied via this newly implemented approach.

Material
As a result of possessing high hardness (>30 GPa), low mass density (2520 kg/m 3 ), and high Hugoniot elastic limit-the elastic wave generated under elastic-plastic-shock phenomena (17-20 GPa) [8], boron carbide (B 4 C) has received considerable attention in ballistic and shock applications [93].Due to its high melting point and thermal stability [94], extreme abrasion resistance [95], and high temperature semiconductivity [96], boron carbide is widely used in refractory, nuclear, and novel electronic applications, respectively; however, its performance is hindered by one or more of a number of inelastic deformation mechanisms, including deformation twinning [97], stress-induced phase transformations [98,99], and various anisotropic fracture behaviors [100] when subjected to mechanical stresses exceeding their elastic limit.In the literature, the key failure mechanisms in boron carbide (e.g., cleavage fracture [101] and twinning [102]) are commonly studied experimentally using numerous characterization techniques (e.g., transmission electron microscopy [103] and Raman spectroscopy [104]).Fracture in the form of shear failure, cavitation, and cleavage has been confirmed from atomic simulation results, either via first principles or molecular dynamics simulations [105,106].Finite deformation continuum models, such as cohesive zone models for fracture [107] and crystal plasticity [108] have also been used to investigate inelastic deformation in single and polycrystalline boron carbide.The present time-evolved phase-field model seeks to engineer the next generation of anisotropic boron carbide-based ceramics by understanding the important plastic deformation and brittle fracture mechanisms that govern its high rate performance.

Order parameter
The main desired feature of the proposed model is to introduce an order parameter η assigned to each material point X to represent fracture, where η = 0 indicates undamaged material, η = 1 fully damaged material, and η ∈ (0, 1) denotes partially degraded material.This variable is commonly assumed to be at least C 2 -continuous with respect to X according to the diffuse interface theory [109,110].The fracture order parameter varies in time and is subjected to, in general, time-dependent boundary conditions.In addition, the global irreversibility constraint of crack evolution is satisfied locally by ensuring a positive variational derivative of the crack surface function and a positive evolution of the crack phase field [76].With this fundamental constraint, cyclic loading can be performed [73]; however, considering fatigue is not in the scope of the current study.

Kinematics
We use standard continuum mechanics notation and understand a summation over repeated indices.The reference, stress-relaxed intermediate, and current configurations are denoted by B 0 , B t , and B, respectively.For the balance of momentum, the computational domain will be B 0 with its closure ∂B 0 .On Neumann boundaries, ∂B 0 N i , the gradient of the solution is known by the given traction vector and on Dirichlet boundaries, ∂B 0 Di , the solution itself (displacement) is given.The motion from the reference position X to the current (deformed) position x = X + u is given by the displacement tensor of rank one, u = u(X, t), as a function in space, X, and time, t.The deformation gradient, F = ∇ x, is a compatible (∇ × F = 0) non-singular two-point tensor; ∇ is the gradient operator in the undeformed state.The right Cauchy-Green deformation tensor and the Jacobian determinant read respectively.Then, the strain measure is nonlinear, called Green-Lagrange strain tensor, as Therefore, we model the geometric nonlinearities accurately.Even if the material model is elastic, the governing equation is nonlinear.A fourth-degree potential ϕ(a, η) allows an analytical solution for a propagating interface, leading to correct introduction of the interfacial stress [111] ϕ(a, η where a is a constant material parameter.This potential should be monotone for 0 ≤ η ≤ 1 (a is chosen between 0 and 6), and satisfy the antisymmetry and the following conditions [112] ϕ(a, 0 To satisfy the antisymmetry condition (i.e., ϕ(a, 1 − η) = 1 − ϕ(a, η) [113]) one must impose a = 3 and obtain where φ(η) is the continuous interpolation function used for dynamic fracture interface simulations.

Free energy and dissipation inequality
The dissipation inequality is derived by following the thermodynamics of irreversible phenomena [114].Considering the balance of energy and subtracting the balance of kinetic energy constructed from the balance of momentum, we obtain the balance of internal energy for each material point as Here, ρ 0 is the mass density in the initial configuration, U is the specific (per unit mass) internal energy, Q is the (heat) flux term across the orthogonal direction, r is the (known) specific volumetric heat supply rate per unit mass, the right-hand side, P : ḞT , is a production term, and P is the first Piola-Kirchhoff stress tensor.We begin with a rather general model and define entropy flux as the heat flux times inverse of temperature, T , as well as an entropy supply term as the internal energy supply term multiplied by 1/T .Hence, we obtain the balance of entropy where specific (per mass) entropy, s, and the entropy production, Σ , need to be defined.The second law of thermodynamics asserts that the entropy production is zero for reversible and positive for irreversible processes.By replacing the supply term in Eq. ( 6) and using divergence theorem, we obtain If we model the heat flux by the Fourier's law where κ T is the positive-definite heat conductivity tensor, then the entropy production's last term is always positivewe emphasize that T is in kelvin.By substituting U with the specific Helmholtz free energy, ψ = U − T s, and assuming that the heat conduction and other thermomechanical processes are mutually independent [115], the mechanical dissipation inequality reads In thermoelasticity, the free energy depends on strain and temperature.Next, we extend the model for the free energy in such a way that the phase-field approach is acquired.We prescribe the constitutive relation for ψ by assuming dependency also on the fracture order parameter and its first derivative, ∇η.Substituting ψ into Eq.( 10), using an integration by parts, and regrouping by means of variables in the energy F, T , η, and ∇η, we obtain Since this inequality holds for every process, we finalize the derivation of the constitutive equations for reversible processes generating zero entropy production We may write the inequality, with a so-called thermodynamic flux, ρ 0 X , conjugate to a so-called thermodynamic force, η.Hence, we model the thermodynamic flux depending on the thermodynamic force in such a way that their multiplication remains positive.Therefore, this model is for irreversible processes.The simplest model is a linear one by choosing a constant mobility parameter L η , which is also called the Ginzburg-Landau (evolution) equation.L η is a positive (kinetic) coefficient steering an energy dissipation in the form of heat during crack propagation [36].
Although taking into account an anisotropic kinetic coefficient, which depends on various free energy functional parameters (e.g., temperature or interface orientation), is required to accurately describe the interface kinetics [116], the kinetic energy here does not consider the velocity of the damage order parameter.

Governing equations
According to the strength and fracture toughness of B 4 C being constant from room temperature to 1500 • K [117], an isothermal process is assumed.Therefore, we aim for solving the balance of momentum and evolution equations for the phase-field parameter.Also, we emphasize that the model of the system accounts for inertial effects, since it is essential for modeling unstable crack propagation in brittle materials [118].We neglect the deformation due to weight and set the gravitational body force to zero.In other words, for an isothermal and dynamic case, the deformation is caused by the mechanical loading on boundaries such that the governing equations become where displacement, u, and phase-field variable are unknowns to be solved.We stress that the whole formulation is acquired by choosing the Helmholtz free energy adequately.We propose to use the following free energy where g(η) is the degradation function, ψ e and ψ ∇ are the mechanical energy and interfacial energy per unit mass, respectively.The degradation function indicates that the mechanical energy of the structure at each position degrades by the order parameter The constant ζ ensures a minimal residual stiffness for fully fractured materials.This constant is also required to ensure the numerical stability and to avoid the singularity of the first part of the energy defined in Eq. ( 16) [119].
The quadratic degradation of elastic energy has likewise been used in a number of other phase-field and gradient damage models [34,120].For the elastic strain energy density in the initial configuration, we use a quadratic energy description leading to a linear material model The fourth-order stiffness tensor, C, depends on the elasticity tensor of the perfect (virgin, without damage) material C(η = 0) and the fracture order parameter as For the interfacial energy density ψ ∇ , the following decomposition is used For cleavage fracture, which is the primary failure mode in boron carbide, we choose the terms in Eq. ( 20) as where B = Υ /l is the ratio of fracture surface energy Υ and crack thickness or regularization length l, ω 0 = Υl is a material constant, β is the cleavage anisotropy factor, and M is the orientation of the cleavage plane, which is known a priori [78,79].The cleavage plane can be a plane of low surface energy or low intrinsic strength in the crystal [121].The parameter β penalizes fracture on planes not normal to M so that β = 0 results in isotropic damage.The regularization length is taken as the cohesive process zone for shear failure [122] where µ 0/2π is the theoretical shear failure strength, and ν 0 = (3k 0 −2µ 0 ) /(6k 0 +2µ 0 ) [123].By using the aforementioned material modeling and strain definition in Eq. ( 18), the governing Eqs. ( 15) of displacement read as with For the phase-field method, in the case of a homogeneous material ρ 0 = const.|X , we have

Variational formulation
Governing equations are replaced by their weak forms.We follow the standard techniques, so-called variational formulation, for generating weak forms and then solve them numerically by means of the finite element method [124].The space discretization is incorporated by approximating fields, u and η, by spanning over nodal values after a triangulation of the computational domain, Ω , with its closure, ∂Ω , into finite elements.For the sake of a simpler notation, we skip a notational change for approximated fields, since their analytical and discrete representations never occur in the same formulation.We emphasize that all unknowns, { u, η } , are solved in a monolithic manner, therefore, the Hilbertian-Sobolev space, H 1 , on a triangulation τ , as We use a discretization using Lagrange elements and generate piecewise continuous polynomials that are adequate for approximation in H 1 .This triangulation is denoted T and consists of non-overlapping triangles, τ .These standard FEM elements are of order a = 2 and b = 1 such that we use linear elements for the phase-field and quadratic elements for displacement.As is common in the Galerkin approach, the same space from above is used in defining test functions For time discretization, we use a Euler backwards scheme for the order parameter and velocity.For example, for the current unknown value, η = η(t), we utilize the solution from one time step before, η 0 = η(t − ∆t), and hence Using a constant time step, ∆t, the second time derivative of the displacement can be written as where u 0 i and u 00 i denote the displacement solutions at previous time steps t n−1 and t n−2 , respectively.Multiplying governing equations by test functions, generating integral forms, and then integrating by parts where necessary, we obtain where a traction vector, t, is given on Neumann boundaries, ∂Ω N .Analogously, we construct a weak form for the phase-field governing equation, where we use integration by parts only for the second derivative in η term The implementation solves the nonlinear weak form For the order parameter, we consider the homogeneous Neumann boundary conditions as a result of assuming a phase-independent energy of the external surface.This assumption leads to the orthogonality of the fracture regions to the boundary because their conjugate is proportional to the normal gradient.The derivation of evolution equations for damage and plastic kinematic variables and the different interpretations about the boundary conditions in the gradient variational approach to damage are studied in [125,126] and the references therein.The solution method is based on a Newton-Raphson solver.In each iteration, the weak form depending on unknowns, P = {u, η}, and their corresponding test functions, δ P = {δu, δη}, is linearized by an expansion around the unknowns from the last iteration, P, in order to calculate P + ∆ P. Here, the solution is for the increment, ∆ P, where the problem is linear by cutting the expansion at second order terms F( P + ∆ P, δ P) = F( P, δ P) + ∇ P F( P, δ P) where the derivative in unknowns is established by the Gateaux (or directional) derivative This standard formulation is beneficial but cumbersome for implementation if the latter derivative is determined manually.Therefore, on a higher level, we use open-source packages from SyFi in FEniCS [127,128] allowing to obtain this part via the symbolic derivative such that the weak form's nonlinearity may be complicated, yet the implementation remains the same.

Numerical simulations demonstrating the instability of a tensile loaded fast moving crack in a single crystal boron carbide
The weak forms in Eq. ( 32) are nonlinear and coupled.We present a monolithic implementation of this transient, fully coupled system of equations by using open-source packages from the FEniCS Project [129].We emphasize that the literature is often suggesting a staggered scheme, see for example [130][131][132] for implementations in FEniCS.The staggered solution solves many smaller problems than one larger, which is faster since the computational cost increases exponentially [133,134].However, in a staggered algorithm, several iterations are necessary for solving one time step in order to ensure that coupling between unknowns is fulfilled with the chosen accuracy.Generally speaking, for highly-coupled systems, a monolithic approach is more feasible.Herein, the implementation for this study involves using a monolithic approach, where displacement and phase fields are solved at once.As aforementioned, the linearization is done automatically by means of a symbolic derivation which is proven to be more reliable [135,136].The conjugate gradient method with a Jacobi preconditioner from PETSc packages has been employed for solving the nonlinear equations.The simulation has been performed by a computing node using Intel Xeon E7-4850, in total 64 cores each with the 40 MB cache, equipped with 256 GB memory in total, running Linux Kernel 5 Ubuntu 20.04.The code is written in Python, although the FEniCS software wraps the formulation to a C++ code and solves as a compiled program.Therefore, yet efficient in developing the code, all computation is running in parallel such that large scale problems are possible to solve under the current implementation.
A popular 2-D example under different loading conditions is demonstrated next in order to simulate unstable crack propagation observed in ceramic boron carbide.The results are adequate, qualitatively and quantitatively.The material properties used in the simulations are shown in Table 1 for B 4 C.The mesh size, h, is such that 10 elements are considered at the interface to resolve the sharp variation along the interface width.Experience has shown that this size provides sufficient accuracy without over resolving the crack.Analysis of the computational solution of this example improves understanding a rapid failure at length-and time-scales not feasible to demonstrate experimentally.These extreme conditions are so challenging to detect experimentally, we must rely on accurate multiphysics simulations.Hence, we stress the importance of the choice of a monolithic solution strategy.
In this example, we model a pre-notched rectangular plate loaded dynamically in tension.This phenomena has been extensively studied by experimental approaches [141,142].Regarding phase-field models, the mode I branching problem has been treated in Borden et al. [40] and Bleyer et al. [143].Both investigations used the small strain framework for precisely locating when and where branching occurs; however, the material is exposed to large strain before and during fracture at the nanoscale [144].The geometry and loading conditions are depicted in Fig. 1.
The order parameter distribution along the x-axis is illustrated in Fig. 2. The mechanical loading of amplitude σ t is established by a traction vector vertical to top and bottom boundaries, σ t n, where n is the surface normal outward from the continuum body.This traction is applied at the initial time step and held constant throughout the  simulation.We investigate three cases with the normal traction amplitude chosen as 1.5 GPa, 1.7 GPa, and 1.9 GPa.The regularization length is l = 0.5 nm, and plane strain is assumed.The highlighted region shows that the value of the fracture order parameter is increasing until reaching to the final value η = 1, where the actual crack tip forms at t = 1.3 ps, and after that the crack starts to grow and propagate along the crack length.Next, distribution of the phase-field variable is analyzed.At different time steps under the same loading scenario, the damage field along different vertical slices is plotted in Fig. 3.The kinetic coefficient parameter used in the phase-field model, which is previously unknown for boron carbide, is assumed to be L η = 1000 Pa −1 s −1 and the crack regularization length is l = 0.5 nm.Under the loading σ t = 1.5 GPa, at x = 1 nm distance from the notch tip, the vertical distribution of damage in Fig. 3(a) provides an insight that the crack reaches this distance between 5 ps and 10 ps.Also, there is a small plateau for the interface profile at t = 22 ps which is related to a small deviation of the crack path.At x = 2 nm, the interface profiles demonstrate the motion of the crack more clearly.By moving through the length of the specimen, x = 4 nm, it is clear that the crack thickness also increases by comparing the distance between every two points on order parameter profiles at each time step in Fig. 3(d).In addition, there is a small drop in the peak point of the damage variable at t = 22 ps which shows the existence of branching.
Increasing the tensile stress by using another traction of σ t = 1.9 GPa leads to a completely different response in Fig. 4. In comparison to the previous case, the crack tip is formed faster.At x = 2 nm in Fig. 4, for the time instant t = 5 ps, the peak point of the interface profile is around three times higher than in Fig. 3(b).In addition, the branching happens quicker as well since the distance of bifurcation is smaller in the case of σ t = 1.5 GPa.A crack branching angle of 45 • followed by a straight extension is observed, which is in good agreement with the literature [141,145].In addition, the crack thickness is larger than the initial regularization length (l = 0.5 nm) for the simulation results, showing that the crack tends to grow in its thickness as it propagates towards the right side of the domain.Furthermore, the damage band widening is more visible for this case in comparison with the loading case of σ t = 1.5 GPa (Fig. 3), which is consistent with experiments and other simulations [146].
In order to better understand the effect of stress magnitude under dynamic loading, the change of crack tip location with time is presented in Fig. 5.In order to determine the crack tip position, we consider the position of top right node that has reached the phase field threshold of η max = 0.995 with the origin taken at the pre-notch tip (inset of Fig. 5) and by following the upper-branch tip.Having the coordinates of the nodes with order parameter values greater than η max , the crack tip velocity can be calculated from a linear fitting to three points of the (t n , α n ) curve, where α n is the total crack length [147].For σ t = 1.5 GPa, the evolution is almost linear and no branching is observed.It is also seen that crack tip approaches to a limiting point when increasing the loading.This procedure occurs for earlier times up to t = 7.5 ps and for higher loading scenarios.The initial plateau for different values of the imposed uniaxial stress is related to the time when the crack tip forms, which occurs sooner for higher stress values.
Finally, the effect of interface kinetic coefficients L η on the damage variable at different times is illustrated in Fig. 6.This result is important to predict the evolution of arbitrary morphologies and complex microstructures in dynamic failure of brittle materials without explicitly tracking the positions of interfaces.In this example, the time step remains constant, the crack regularization length is assumed to be l = 0.5 nm, and the plate is subjected to a uniaxial tension stress of σ t = 1.3 GPa.It is shown that increasing the interface kinetic coefficient leads to decreasing the formation and propagation of crack's speed.In addition, there cannot be seen any crack kinking or branching for higher kinetic coefficients, even at later time steps.This result is in good agreement with the expression for the chemical potential on cracked surface defined in [148] due to the existence of the inertial effects and tip splitting for high applied tension.Also, the crack speed is limited by the sound speed for fast growth and all dissipation takes place in the interface region [149], resulting in crack widening perpendicular to the propagation direction which is more visible for lower kinetic coefficients especially before a branching event (i.e., increasing of fracture surface energy).

Conclusion
A robust finite element procedure for solving a coupled system of equilibrium and time-dependent Ginzburg-Landau equations has been motivated by using thermodynamically-sound derivation of governing equations.Use of the variational procedure and thermodynamic consistency of the model ensures that it has a strict relaxational behavior for the free energy; hence, the models are more than a phenomenological description of an interfacial problem when compared to the available literature [150].The dissipation and time scales associated with growth kinetics are also derived and addressed in this manuscript as complementing previous works [69].The model has been used for studying the evolution of fracture in anisotropic single crystal boron carbide at finite strains.For the first time, a monolithic strategy has been employed for solving the coupled mechanical equilibrium and order parameters evolution equations under extreme conditions.As a challenge in computational mechanics, nanometer length scale and picosecond time scale have been demonstrated in simulations.The computational procedures and numerical algorithms are implemented using the open-source platform FEniCS.The present nonlinear finite element code has been developed and used to study the prediction of paths derived from dynamic crack growth under uniaxial tensile stress loading in single crystal boron carbide.The numerical results for all the problems are in agreement with the available experimental data in the literature.The current contribution opens up new possibilities for multi-scale fracture models.In the future, our finite element based phase-field model can be applied for studies of phase transformations (e.g., amorphization) and interaction between plasticity and fracture under high strain-rate loading.As a next step, the current model could be combined with discrete localized plastic flow (e.g., shear band and dislocation pileups) and thermally-activated mechanisms (e.g., melting) to capture the behavior of brittle materials in laser spall experiments.

Code availability
The Python code, generated during the current study, is part of the FEniCS project available at http://www.fenic sproject.org/download,and an example for the computational implementation is available in [151] to be used under the GNU Public license [152].

CRediT authorship contribution statement
Benhour Amirian: Developed the model, Wrote the code, Designed and performed all simulations, Formal analysis, Writing -original draft, Discussed the results.Bilen Emek Abali: Allocated the computational resources, Writing -review & editing, Discussed the results.James David Hogan: Supervision, Funding acquisition, Writing -review & editing, Discussed the results.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.The geometry and boundary conditions for the dynamic crack branching example.The specimen contains a crack and is subjected to a symmetric traction load.

Fig. 2 .
Fig. 2. The dynamic crack profiles versus the crack notch's length x at different time steps for regularization length of l = 0.5 nm under a uniaxial tension stress of σ t = 1.5 GPa.The inset shows the highlighted order parameter's distribution in yellow for t = 1 ps, 1.1 ps, and t = 1.3 ps.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 3 .
Fig. 3. Symmetric damage interfaces under a uniaxial tension stress of 1.5 GPa at different time steps along different vertical slices {(x, y) ∈ R | x = i}.(a) i = 1 nm; (b) i = 2 nm; (c) i = 3 nm; (d) i = 4 nm.The magnified simulations results for each case clearly indicates the crack widening before splitting into two distinct cracks.

Fig. 4 .
Fig. 4. Symmetric damage interfaces under a uniaxial tension stress of 1.9 GPa at different time steps along different vertical slices {(x, y) ∈ R | x = i}.(a) i = 1 nm; (b) i = 2 nm; (c) i = 3 nm; (d) i = 4 nm.The magnified simulations results for each case clearly indicates the crack widening before splitting into two distinct cracks.

Fig. 5 .
Fig. 5. Evolution of crack tip location for different loading.The inset shows the way of determining the crack tip location at different time steps through the domain's central line in the x direction.

Fig. 6 .
Fig.6.Distribution of the damage variable in a square after loading with a size of 50 nm, extracted from the original domain before loading (refer to the Fig.1).Results are shown for five different times of t = 5 ps, 10 ps, 15 ps, 20 ps, and 25 ps (designated on the left), and for four different kinetic coefficients L η under a uniaxial tension stress of σ t = 1.3 GPa and crack regularization length of l = 0.5 nm.

Table 1
Material properties and model constants for B 4 C.