Thermodynamically-consistent derivation and computation of twinning and fracture in brittle materials by means of phase-field approaches in the finite element method

A theoretical-computational framework is proposed for predicting the failure behavior of two anisotropic brittle materials, namely, single crystal magnesium and boron carbide. Constitutive equations are derived, in both small and large deformations, by using thermodynamics in order to establish a fully coupled and transient twin and crack system. To study the common deformation mechanisms (e.g., twinning and fracture), which can be caused by extreme mechanical loading, a monolithically-solved Ginzburg--Landau-based phase-field theory coupled with the mechanical equilibrium equation is implemented in a finite element simulation framework for the following problems: (i) twin evolution in two-dimensional single crystal magnesium and boron carbide under simple shear deformation; (ii) crack-induced twinning for magnesium under pure mode I and mode II loading; and (iii) study of fracture in homogeneous single crystal boron carbide under biaxial compressive loading. The results are verified by a steady-state phase-field approach and validated by available experimental data in the literature. The success of this computational method relies on using two distinct phase-field (order) parameters related to fracture and twinning. A finite element method-based code is developed within the Python-based open-source platform FEniCS. We make the code publicly available and the developed algorithm may be extended for the study of phase transformations under dynamic loading or thermally-activated mechanisms, where the competition between various deformation mechanisms is accounted for within the current comprehensive model approach.


Introduction
Understanding and predicting anisotropic fracture and damage evolution in brittle materials have been long-standing problems in engineering designs. Owing to the advent of novel modeling techniques and the advances in computational capabilities, the usage of accurate and robust numerical methods plays a key role in situations where purely experimental approaches are of high cost and not always readily accessible (e.g., high-energy in situ X-ray computed microtomography [1], in situ electron backscattered diffraction (EBSD) [2], and micro/nano-mechanical testing [3]). In the literature, the simulation of fracture in solids at the atomic scale is commonly treated by molecular dynamics [4], density functional theory [5], or lattice static models that are based on spring networks [6]. Despite addressing nonlinearities at the crack tip, avoiding singularity-related issues, and considering bond breaking between atoms [7,8,9], there exist challenges in using atomistic models to cover the timeand length-scales necessary to analyze the structural response at the macroscale needed for engineering applications.
Conventionally, there are two main categories of numerical approaches that are employed to provide explicit simulations of material failure: (i) discrete crack models (e.g., the discrete element method [10,11], the extended finite element method (XFEM) [12], the cohesive zone method [13], and the cohesive segment methods [14]) in which the displacement field is allowed to be discontinuous across the fracture surfaces, and (ii) smeared (continuum) crack models (e.g., damage models [15], and diffuse interface models [16,17]) that consider a continuous displacement everywhere, assuming gradually decreasing stiffness to model the degradation process. Regardless of showing much success in modeling crack propagation [18], the discrete crack models need additional criteria based on stress, strain energy density, energy release rate, or virtual crack closure techniques [19] to predict the crack initiation (nucleation), growth, and branching in dynamic fracture problems [20]. Further, the sharp representation of cracks requires remeshing algorithms or using the partition of unity method [21], both having their own difficulties in tracking the multiple crack fronts in complex three-dimensional morphologies [22,23].
In the smeared crack approach, regularizing strong discontinuities caused by strain localizations within a finite and thin band leads to a precise approximation of the crack topology [24]. The gradient damage model [25,26], physical/mechanical community-based phase-field fracture models [27,28,29] that traces back to the reformulation of Griffith's principle [30], and peridynamics [31,32,33], which may be regarded as generalized non-local continuum mechanics, fall within this category. Replacing partial differential equations in the phase-field model by integrals in peridynamics allows for topologicallycomplex fractures such as intersecting and branching to be handled in both two and three dimensions [34]. Coupling smeared and discrete crack approaches, for example, the element deletion method [35], the combined non-local damage and cohesive zone method [36,37], and the thick level-set method [38] have also shown promising results in modeling fracture. In the thick level-set method, a discontinuous crack description is surrounded by continuous strain-softening regions, defined by a level-set function to separate the undamaged from the damaged zone [39]. However, the dependence of the results on the finite element mesh and the convergence path of the solutions, for a mesh size tending to zero, results in numerical errors [40].

Phase-field approach
As an alternative approach, the phase-field model has been widely used recently in the context of phase transition processes, ranging from solidification [41] and phase transformation in solids [42] to the modeling of ferroelectric materials [43]. Having the capability to model the microstructural evolution, phase-field modeling has been successfully adopted in the simulation of martensitic phase transformations [44,45], reconstructive phase transformations [46], phase transformations in liquids [47], dislocations [48], twinning [49], damage [50], and their interactions [51,52,53]. Initiated with the celebrated work by Francfort and Marigo on the variational approach to brittle fracture [54], where the total energy is minimized simultaneously with respect to the crack geometry and the displacement field, the concept of applying the phase-field method in fracture mechanics has gained significant interest in the literature [55,56,57,58,59,60,61,62,63,64,65,66,67]. Due to the thermodynamic driving forces, the evolution of interfaces (e.g., merging and branching of multiple cracks) is predicted with no additional effort [68]. Also, being quantitative, material-specific, and simple to couple to other calculations (e.g., stress or temperature [69]) makes phase-field modeling a powerful and flexible method for studying the fracture of single crystalline [70] and polycrystalline materials [71] as well as in granular solids [72]. The high computational cost in the phase-field model due to resolving the gradient term by using sufficiently refined mesh in the damaged zone is straightforwardly tackled by parallel implementations [73] and adaptive remeshing [74].
Using the phase-field model to study the failure mechanisms in brittle materials has recently received increasing attention [75,76,77,78,79]. Brittle solids may fail along grain boundaries or cracks propagate along the constituent phases in the case of geomaterials [80,81]. A fourth-order model for the phase-field approximation in brittle materials leads to a unified model to simulate the mechanics of damage and failure in concrete [82], where an explicit Hilber-Hughes-Taylor-α [83] method is used with a phase-field approach given by a single-well energy potential to describe the fracture behavior. This method is well-established in the mechanics community. In the physics community, on the other hand, the phase-field models are commonly derived by adapting the phase transition formalism of Landau and Ginzburg [84]. For example, Aranson et al. [85] combined elastic equilibrium with the Ginzburg-Landau equation, which accounted for the dynamics of defects, to study the crack propagation in brittle amorphous solids. Another Ginzburg-Landau based phase-field approach restricted to mode III fracture (antiplane shear) was proposed by Karma et al. [27] and Hakim and Karma [16] in the twoand three-dimensional settings, respectively. Considering fracture as a solid-gas transformation, the double-well energy potential appeared in phase-field modeling of damage [58]. Some of the disadvantages of the double-well potential, such as crack widening and lateral growth during crack propagation, can be eliminated by using a single-well term [86]; however, the realistic shape of the stress-strain curves obtained from the experiments or atomistic simulations are more difficult to capture by the single-well free energy density [57].
In numerical implementations, nonlinear problems with a strong coupling between the equilibrium equation and the phase-field parameter is solved through two approaches: 1. A staggered solution scheme is based on decoupling balance equations and the phase-field problem into the system of two equations that are solved in a subsequent manner [87,88]. The implementation is more modular and two smaller systems to solve is faster than two systems together; computational time increases exponentially. The method is robust due to giving rise to two convex minimization problems, but depending on the coupling (and application), a significant amount of staggered iterations may be required at a fixed loading step, thus resulting in a higher computational cost [89].

2.
A monolithic solution where all variables are solved at once (simultaneously) [90]. In some cases, for example highly coupled systems, the monolithic solution is more efficient as a result of (in total) less Newton-Raphson iterations [91]. To the best of our knowledge, no studies have focused on solving the Ginzburg-Landau based phase-field problem concerned with predicting the twinning and fracture behavior of brittle materials by using a monolithic scheme; this is addressed in the current paper.

Goals and outlook
In the present study, we seek to extend the Ginzburg-Landau phase-field approach to predict fracture and twinning in single crystal anisotropic brittle materials (e.g., magnesium, Mg, and boron carbide, B 4 C). The advantage of the Ginzburg-Landau approach over the incremental energy minimization method used in Clayton and Knap [70] is that material parameters associated with time scales for interfacial motion enter the model by calibrating it with the most recent Molecular Dynamics simulations [92], making it a more general model for studying the deformation mechanisms of intrinsically brittle materials (e.g., single crystalline Mg and B 4 C). This is an important consideration for deformation twinning since the propagation speed of twin boundaries can be difficult to measure, and could even be supersonic if the driving stress is sufficiently large [93]. In addition, this work focuses on derivation of governing equations and solving them monolithically in order to increase the accuracy for applications with strong coupling between mechanics and damage. We solve the nonlinear and coupled differential equations by using the open-source parallel computing platform FEniCS [94]. By following the works on local stress concentrations in nanoscale defect-free volumes or by high pressures [95], as well as shear arising from twinning [96], we develop a nonlinear phase-field theory for elasticity along with anisotropic surface energy [97]. To address this, the governing equations for both small and large deformations are derived because finite rotations may occur under some loading conditions even at small strains, which necessitate considering both regimes in crystallographic theory [98]. A new decomposition for the strain energy density based on [99] is proposed to reproduce the experimentally-observed crack propagation under compressive loading, and the simulated results are compared with analytical solutions. In these comparisons, a double-well energy potential is considered for studying the fracture behavior (e.g., crack initiation, growth, and propagation) and twinning in anisotropic brittle solids.
The remainder of this paper is outlined as follows. In Section 2, we briefly describe the chosen materials. Continuum mechanics and thermodynamically sound derivation of equations are shown in Section 3. Variational formulation and the finite element method is explained in Section 3.5. Results and representative material properties along with the discussion of phase-field simulations are reported in Section 4. The conclusions of the study are drawn in Section 5.

Materials
The focus of the present paper is to model the deformation behavior of magnesium and boron carbide single crystals. The low ductility of these materials suggests to comprehend them as brittle elastic materials with large driving forces for dislocation glide, leading to other mechanisms such as phase transformations [100], deformation twinning [101,102], and fracture [103].

Magnesium
Having low mass density (∼23% of steel and 66% of aluminum), high strength, and durability for a wide range of temperatures in high performance automotive and aerospace applications, magnesium and its alloys have attracted considerable attention in recent years [104,105]. Mg alloys tend to be brittle due to their limited number of dislocation systems [106]. As a result of possessing low-symmetry crystallographic structure and larger critical resolved shear stress (CRSS), twinning is the dominant deformation mode in magnesium as it is subjected to twinning-favored loads stretched along [0001] [107,108,109], resulting in transitions in the material behavior at high strain rates [110]. A previous study indicated that the formation of intersecting twins may improve the ductility of Mg alloys [111]. Therefore, understanding and predicting the twinning behavior during plastic deformation of magnesium is critical towards the realization of next-generation lightweight metallic materials for application in automotive and defense applications. In order to investigate twinning in magnesium, various techniques such as high-resolution transmission electron microscopy [112], visco-plastic self-consistent polycrystal models [113], elasto-plastic self-consistent polycrystal models [114], molecular dynamics simulations [115], crystal plasticity models [116], and quasi-static phase-field models [49] have been employed. In this paper, the fracture and twinning behaviors of single crystal magnesium are studied using an advanced time-dependent phase-field theory by numerically solving engineering problems.

Boron carbide
As a result of possessing hardness above 30 GPa, low mass density (2.52 g /cm 3 ), and high Hugoniot elastic limit (17)(18)(19)(20), boron carbide (B 4 C) has received considerable attention in ballistic applications [117]. Due to its high melting point and thermal stability [118], favorable abrasion resistance [119], and high temperature semiconductivity [120], boron carbide excels 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 [121], stress-induced phase transformations [122,123], and various fracture behaviors [124] when subjected to mechanical stresses exceeding their elastic limit. The key failure mechanisms in boron carbide (e.g., cleavage fracture and twinning) are commonly studied experimentally using numerous characterization techniques (e.g., transmission electron microscopy [125] and Raman spectroscopy [126]). 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 [127,128]. Finite deformation continuum models, such as cohesive zone models for fracture [129] and crystal plasticity [130] 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 armor ceramics by understanding the important plastic deformation and brittle fracture mechanisms that govern its high rate performance. The current framework does not incorporate slip for B 4 C due to having very large resistance to dislocation glide in certain directions at low temperatures [131].

Formulation
In this section, we develop a model for single twinning and fracture systems in solids based on thermodynamical derivation. The present approach extends the model of Clayton and Knap [49,70,132,133] by accounting for the time-evolution of order parameters towards an equilibrium state for predicting the twinning and crack paths in anisotropic single crystal materials adequately. This allows the study of spatio-temporal fluctuations of order parameters (e.g., twinning and fracture variables), as well as the nanoscale dynamics that govern various pattern forming phenomena [134]. Moreover, the interfaces, their propagation, and interactions, which are the most important features governing the formation of microstructures in materials, can be studied via this newly implemented approach. It is worth mentioning that the present work does not address plastic slip distinct from the motion of twinning partials inherent in deformation twinning. We refer to [135,136] and references therein for further details on ways to simultaneously address plastic slip and twinning. In what follows, we are interested in elastic twinning (i.e., the reversible nature of the corresponding deformation of the crystal) in which a twin appears and grows in the crystal lattice due to the presence and to the increase of an external load, and escapes from the crystal upon removal of the load [137].

Order parameters
The main desired feature of the proposed model is to introduce an order parameter η assigned to each material point X for the description of the twinning, η = 0 within the parent (original or untwinned phase at (X, t)) elastic crystal, whereas η = 1 denotes the twin. The twin boundary zone is determined by the diffuse interval, 0 < η < 1. Another order parameter, ξ, is used to represent fracture, where ξ = 0 indicates undamaged ("virgin") material, ξ = 1 fully damaged material, and ξ ∈ (0, 1) partially degraded material. Both of these state variables are commonly assumed to be at least C 2 -continuous with respect to X according to the diffuse interface theory [138,139]. They also vary in time and are subjected, in general, to time-dependent boundary conditions. In addition, the global irreversibility constraint of crack evolution is satisfied by ensuring locally a positive variational derivative of the crack surface function and a positive evolution of the crack phase field [85]. Without this fundamental constraint, no cyclic loading can be performed [61,72]; however, cyclic loading is not in the scope of the current study.

Kinematics
We use standard continuum mechanics notation and understand a summation over repeated indices. A continuum body as in Figure 1, composed of many grains, is considered in the reference frame, This framework is simply the initial configuration-before loading.
For modeling purposes, we introduce a stress-relaxed intermediate frame, B * , and a current frame, B.
Since the formulation is established in the reference frame, it is a material system where coordinates denote particles. In this frame, the mass density is given as a function in X, we circumvent solving the balance of mass. For the balance of momentum, the computational domain will be B 0 with its closure ∂B 0 . On Neumann boundaries, ∂B 0 N i , gradient of the solution (traction vector) is known and on Dirichlet boundaries, ∂B 0 Di , the solution (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 = ∇ 0 x, is multiplicatively decomposed as where F E is the recoverable elastic deformation work-conjugate to Piola stress, analogously, F η is the deformation associated with structural defects, twinning in the current study, evolving within the material. It is important to highlight that twinning is distinguished from plastic slip. First, the former occurs by collective motion of defects that preserves the particular orientation relationship between the twin and original phases [140]. In addition, twinning is unidirectional while usually slip is not [141]. In contrast to F , which always satisfies compatibility conditions ∇ 0 × F = 0, the deformation maps F E and F η are generally not integrable to a vector field, when considered individually, as a result of existing crystal defects [142,143]. In other words, these two deformations are generally anholonomic (∇ 0 × F η = 0 and ∇ × (F E ) −1 = 0) [144,145]. Additional structural changes may be included in this theory for representation of other defects, such as point defects [146] or dislocation slips [147,148]. The kinematics for twinning in simple shear is given as [149] F η (η) where s and m are the orthogonal unit vectors in the directions of twinning and normal to the twinning plane, respectively; and γ 0 is the magnitude of the maximum twinning shear. All functions are defined in the reference frame, since we have a material system. For small deformation, this distinction is negligible. The continuous interpolation function φ(η) is obtained from a general representative function ϕ(a, η) within a fourth-degree potential [44] defined as where a is a constant parameter-in order to ensure that ϕ(a, η) is a monotonous function, a is chosen between 0 and 6. The functions ϕ(a, η) and φ(η) should be monotone for 0 ≤ η ≤ 1 and satisfy the following conditions [44] ϕ(a, 0) = 0, ϕ(a, 1) = 1, , which obeys the antisymmetry condition, i.e., φ(1 − η) = 1 − φ(η) [150]. As usually assumed, the plastic deformation is deviatoric such that the volume remains the same, det F η = 1. Therefore, the Jacobian determinant reads, J = det F = det F E . With the right Cauchy-Green deformation tensors, C ij = F ki F kj and C E ij = F E ki F E kj , we obtain the Green-Lagrange total and elastic strains, respectively. They are now used, as follows: In the case of small deformations, by using the standard linearization approach, we obtain With the same approach, by discarding all nonlinear terms (including displacement gradient multiplied by s or m), we determine from Eq. (6), in the case of small deformations, By inserting Eq. (7), we obtain In this way, we aim for using F and F η in large deformation simulations; analogously, ε and ε η in small deformation simulations. A comparison will lead to the significance of nonlinearity in twinning simulations.

Constitutive equations
In order to derive the constitutive equations, we follow thermodynamics of irreversible phenomena and refer to [151] for historical remarks. By starting with the balance of energy and subtracting the balance of momentum, for an arbitrary volume V 0 in the undeformed configuration B 0 , we obtain the balance of internal energy: with the specific internal energy, u in J/kg, its (heat) flux term, Q in W/m 2 , across the direction N (surface normal outward the volume), a specific volumetric heat supply rate, r in W/kg, and a production term defined by Piola stress, P in N/m 2 , and deformation gradient rate, among others see [152,Sect. 2.4] for a straight-forward derivation. Temperature T in K, and specific entropy s in J/(K kg), are related by the global entropy balance equation, in the undeformed configuration, with the assumption that the entropy flux is 1/T times heat flux leading to where the entropy production, Σ, is zero for reversible and positive for irreversible processes. This assertion, Σ ≥ 0, is the Second Law of Thermodynamics. Since T ≥ 0, we immediately acquire T Σ ≥ 0. By replacing the supply term in Eq. (10) and using Gauß-Ostrogradskiy (divergence) theorem, we obtain Now, by using the Helmholtz free energy, ψ = U − T s, and using temperature as an independent thermodynamic parameter instead of entropy, the entropy production related term reads As the latter is positive, several restrictions are possible for constitutive relations. We choose the simplest possible constitutive relation for the heat flux. The linear dependency is called Fourier conduction law, where κ is the (symmetric and positive-definite) thermal conductivity tensor, it may depend on temperature [136] but is constant in temperature gradient. Hence, we reduce T Σ ≥ 0 to the following inequality We start modeling by assuming the free energy dependency on order parameters and their first derivatives, ψ = ψ(F , T, η, η ,i , ξ, ξ ,i ), leading to By using the latter in Eq. (14), regrouping by means of variables in the energy F , T , η, η ,i , ξ, and ξ ,i , we obtain This model is the simplest one and we obtain after subsequent Gauß-Ostrogradskiy theorems The boundary conditions for the evolution of the order parameters are obtained by assuming a phase-independent energy of the external surface ∂V 0 of V 0 With this assumption, fractured or twinned regions are always orthogonal to the boundary because their conjugate is proportional to the normal gradient. We refer to [153] for different interpretations about the boundary conditions in the damage gradient approach. However, in the phase-field theory, one can use a stricter approach by introducing a generalized surface force conjugated to the rate of change of the order parameter. This generalized surface force balances the terms that appears due to the dependence of the free energy on the gradient of the order parameter, and this will give us more general boundary conditions [154]. The inequality (17) has to hold for any process, the first two terms may only be zero, since T • and F • are in general not restricted and calculated by balance of entropy and momentum, respectively. Hence, we acquire the well-known relations where the second term is known as Castigliano's theorem. In the case of small deformations, the latter differentiation reads σ ij = ρ 0 ∂ψ/∂ε ji . The Second Law of Thermodynamics is satisfied for all processes by choosing (mobility) parameters L η and L ξ in order to achieve Ginzburg-Landau (evolution) equations which reduce Eq. (17) to The L η and L ξ are positive kinetic coefficients for twinning and fracture evolution, respectively.

Governing equations
For a given temperature, i.e. for an isothermal process, the balance of entropy is fulfilled and we aim for solving the balance of momentum and evolution equations for order parameters. We make further assumptions and model the system without inertia and body forces. In other words, for an isothermal, quasi-static case, the deformation is caused by the mechanical loading on boundaries such that the governing equations become As usual in mechanics, we search for displacement, u, from the balance of momentum in Eq. (22) 1 . Equation (22) Technically, ψ ∇ represents the gradient energy per mass. The mechanical energy degrades by the order parameter ξ denoting the microporosity (ξ = 1 means fracture) of the structure in each position. This field function, ξ, is used to obtain a degradation function phenomenologically The constant ζ ensures a minimal residual stiffness for fully fractured materials. The quadratic degradation of elastic energy has likewise been used in a number of other phase-field and gradient damage models [155,156,157,158]. The reflection or rotation of the reference frame of the crystal lattice commensurate with twinning should be taken into account for anisotropic elastic constants [159]. By using Green-Lagrange strains in Eq. (5), for the deformation energy density, W M , we use a quadratic energy description The stiffness tensor, C, has minor and major symmetries, C ijkl = C jikl = C klij . Also it depends on the twin stiffness C T and initial stiffness C P by the phase-field approach, The order parameter, η, is used to determine the amount of each phase in various position. Elastic coefficients of the fully twinned crystals, η = 1, are related to those of the untwinned state, η = 0, by where Q is the reorientation matrix transforming the original lattice to twin lattice within a centrosymmetric structure Type I and type II twins differ in reflections or rotations of the lattice vectors in the twin and parent phase. In the case of homogeneous materials, ∇ρ 0 = 0, we simplify the notation and use this for the gradient energy density, W ∇ = ρ 0 ψ ∇ , and use the following decomposition: The first term consists of a standard double-well potential [160,161,154] where A = 12Γ /l characterizes the energy barrier between two stable phases (minima), relating to the equilibrium energy per unit area, Γ, and thickness, l, of an unstressed interface [49]; ι(ξ) is a coupling degradation function which degrades with the fracture parameter ξ. It is assumed that ι(ξ) = g(ξ), meaning that the twin boundary energy and the elastic deformation energy degrade with damage according to the same quadratic function. The regularization length is taken as the cohesive process zone for shear failure [162] l = 16πΥ where Υ is the fracture surface energy, µ0 /2π is the theoretical shear failure strength, and ν 0 = (3k0 − 2µ0) /(6k0 + 2µ0) [163]. The second term on the right-hand side of Eq. (29) follows from the Cahn-Hilliard formalism [138] where κ ij = κ 0 ι(ξ)δ ij is a diagonal tensor of rank two, and κ 0 = 3Γl /4 is a gradient energy parameter. For cleavage fracture, which is the primary failure mode in boron carbide, we choose the terms in Eq. (29) as follows where B = Υ /h is the ratio of fracture surface energy and crack thickness, ω 0 = Υh is a material constant, β is the cleavage anisotropy factor, and M is a unit vector in material coordinates that is normal the cleavage plane [87,88]. The cleavage plane can be a plane of low surface energy or low intrinsic strength in the crystal [164]. Generally, there is no predefined relation between cleavage and twinning planes. Orientations of these planes are specified a priori and may or may not coincide [133]. The parameter β penalizes fracture on planes not normal to M so that β = 0 results in isotropic damage. This formulation has been used in recent continuum models of fracture as a result of its ability to converge to the correct surface energy of a singular surface when the twin boundary thickness tends to zero [165,166].
By using the aforementioned material modeling and strain definition in Eq. (2), the governing equations (22) read for displacement For phase-fields, in the case of a homogeneous material ρ 0 = const.
By using Eqs. (1), (2), we have For a better analogy, we misused the notation for the same degradation function in Eq. (24) as follows: Finally, we obtain the following equations to solve numerically, and

Variational formulation
We follow the standard techniques for generating weak forms to solve numerically by means of the finite element method [167]. The space discretization is incorporated by approximating fields, u, η, ξ, by spanning over nodal values after a triangulation of the computational domain, Ω, with its closure, ∂Ω, into finite elements. For simplicity, 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 n , with the polynomial order, n, as follows: We use mixed spaces for quantities, depending on the problem, displacement (linear or quadratic) and phase-fields (linear) standard Lagrange finite element. On each node, unknowns read 2 + 1 + 1 = 4 degrees of freedom (DOF) in 2-D and 3 + 1 + 1 = 5 (DOF) in 3-D space. As usual in the Galerkin approach, we use the same space for test functions, δu, δη, δξ , where they vanish on Dirichlet boundariesV = δu, δη, δξ ∈ [H 1 (Ω)] DOF : δu, δη, δξ To ensure the irreversibility of the order parameter rates in our simulations,η andξ, we restrict them to be zero if they become negative. For the time discretization, we use Euler backwards scheme for order parameters, for example where ∆t is the time step. For simplicity we use constant time steps. This method is implicit, hence for real valued problems stable, and converges to the correct solution. 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 .
where we have employed the fact that N i η ,i vanishes at the boundary of Ω in connection with Eq. (17). Analogously, we obtain The implementation solves the nonlinear weak form after a symbolic derivation and Newton-Raphson iterations.

Multiphysics simulation
The weak forms in Eq. (48) are nonlinear and coupled. We have implemented a transient, fully coupled solution strategy by using open-source packages from the FEniCS Project [168]. We refer to [169,170,171] for implementations in FEniCS by means of a staggered scheme. Herein, the implementation uses a monolithic approach, where displacement and phase fields are solved at once. Staggered solution solves many smaller problems than one larger, which is faster since the computational cost increases exponentially. However, in a staggered algorithm, several iterations are necessary for solving one time step in order to ensure that coupling between unknowns are fulfilled. Generally speaking, for highly-coupled systems, a monolithic approach is more feasible. For the linearization, we use a standard Newton-Raphson approach. The linearization is done automatically by means of a symbolic derivation that allows the user to write the weak form without going through the error-prone linearization process by hand [172,173]. 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 very efficiently. In short, the problem-specific parts of the computer code used to perform the simulations have been generated automatically from a high-level description that resembles closely the notation used in this work.
Examples under different loading conditions in two-dimensional samples are demonstrated next in order to simulate deformation mechanisms observed in metallic magnesium and ceramic boron carbide. The results are adequate, qualitatively and quantitatively. The material properties used in the simulations are shown in Table 1 for Mg and B 4 C. Five independent second-order elastic constants [174] are listed by using Voigt notation, with indices I, J from 1 to 6, as follows: The bulk modulus listed below for each undamaged material is obtained by [175] k = (C 11 + C 12 ) C 33 − 2C 2 13 According to the primary inelastic mechanisms for each material mentioned before, three different problems are simulated and discussed in the following to represent the degenerate cases: (i) Twin propagation in two-dimensional single crystals magnesium and boron carbide in Sect. 4.1; Fracture is suppressed in this condition by assuming ξ = 0.
(ii) Analysis of twinning induced by a crack in magnesium under pure mode I or mode II loading in Sect. 4.2; Similar to the previous case, fracture is not calculated, ξ = 0. These examples demonstrate that we generate knowledge about mechanical deformations in very small length-scales and extreme loading rates causing a twin or crack initiation and propagation at very small time scales. These extreme conditions are challenging to observe experimentally, where we rely on accurate multiphysics simulations as presented herein.

Twin growth and propagation
Nucleation and evolution of deformation twinning in a single crystal of magnesium (in Sect. 4.1.1) and boron carbide (in Sect. 4.1.2) are presented in a two-dimensional domain in plain strain conditions to be depicted in Fig. 2. For validation, the model is initially solved for elastically isotropic pure magnesium single crystals with the properties listed in Table 1. The isotropic elastic approximation appears reasonable because magnesium single crystals are not strongly anisotropic elastically [181]. A circular twin nucleus, η = 1, of initial radius a = 3 nm is embedded in a rectangular domain with a surrounding parent material, η = 0. The domain is of 40 × 40 nm in size for the magnesium simulations. The initial radius of the twin embryo is set to 3 nm as a result of the fact that a bifurcation from circular to elliptical shape occurs for a radius of 3.2 nm, corresponding to the analytical sharp interface solution [182]. The lattice orientation vectors are in the form where θ denotes the orientation of the habit plane. Also, according to the following matrix-form gradient coefficient both isotropic (κ 11 = κ 22 = κ 0 ) and anisotropic ( κ 11 2 = 2κ 22 = κ 0 ) twin boundary surface energies are employed in different simulations in order to explore their effects as well as for validation purposes. The following simple shear with Dirichlet boundary conditions on ∂B D for top and bottom boundaries are used where Λ = 0.08 is the magnitude of applied shear for all simulations in the following section. The twin growth to the boundary is inhibited by the displacement boundary conditions. The order parameter gradients also vanish at the boundaries due to the Neumann boundary conditions defined in Eq. (18). Figure 3 shows contour plots demonstrating the spatial distributions of numerical results for the growth of a circular twin embryo in a single crystal magnesium with an orientation of the habit plane θ = 0.

Twin embryo propagation and growth in single crystal magnesium
The embryo is undergoing a simple shear at 8% displacement prescribed on the top. Parameters of interest include the twin order parameter (i, ii), y displacement (iii, iv), and shear stress (v, vi). Each image pair considers both small (left side) and large strains (right side), as well as isotropic (a, b) and anisotropic surface energies (c, d). For this case, there is no significant difference in the simulation results between linear (left side) and nonlinear (right side) elasticity. The results are shown at time instants of t = 50 ps and t = 500 ps to show the evolution of the twin's morphology. The mesh of the rectangular domain includes 160,000 linear triangular elements. By using a standard h-convergence, we have chosen this particular mesh to deliver mesh insensitive results. The 1 0 1 1 plane and {1 0 1 2} direction are considered as the primary twinning system in magnesium [183].
First, the evolution of the twin order parameter is shown in Fig. 3(a, b)(i, ii) under simple shear with the boundary conditions defined at t = 50 ps and t = 500 ps for small and large deformation with isotropic twin boundary energy. As can be seen, the twin embryo grows until it is repelled by the rigid outer boundaries, where the order parameter is set to zero. Under these numerical conditions, a small orientation of the twin evolution is realized due to the difference in the driving force for twinning, which is a factor of (F η ) −1 . We emphasize that the twin morphology at the final stage is in qualitative agreement with the (static) phase-field results [49] and molecular dynamics simulations [115], thus serving to verify the set of results in Fig. 3(a, b)(i, ii) for the proposed time-dependent phase-field model.
Second, the distribution of the displacement in the y direction for the domain under simple shear loading for small and large strains at different times are depicted in Fig. 3(a, b)(iii, iv). The positive and negative displacement values indicate that the left and right sides of the twinned boundary regions are under compressive and tensile loading, respectively. This distribution in a simple shear is not possible for a homogeneous material with prescribed vanishing y displacement on top. Herein, we stress that the parent to twin phase change introduces a heterogeneity in stiffness parameters across the twin boundaries. By means of this relatively simple simulation, we gather an insight into the material response. Moreover, the range of displacement magnitudes at the very last time instant are lower than those at initial times as a result of inhibiting by the boundaries. The corresponding evolution of the shear stress for small and large strains with consideration of the isotropic surface energy at various times are illustrated in Fig. 3(a, b)(v, vi). Investigating the shear stress distribution improves our knowledge of the redistribution of high local stress, resulting from twinning [184], and this provides new insights into demonstrating the driving force for the propagation and growth of twin within a small region in the microstructure. Third, in Fig. 3(a, b)(v, vi), one component of the stress tensor (shear) is shown. Again, in a homogeneous material under simple shear conditions, a constant shear stress is created. The distribution is caused by the twinning, as visible by the shape compared to the twin distribution.
Fourth, the effect of anisotropic surface energy is studied in Fig. 3(c, d)(i -vi). For the twin order parameter in Fig. 3(c)(i, ii), the equilibrium shape of the twin embryo under small strains is wider in the horizontal direction (parallel to the habit plane) and flatter in the vertical direction at t = 50 ps as compared to the isotropic energy case shown in Fig. 3(a)(i, ii). This behavior has been observed previously in the time-independent phase-field approach [49], where the results are in agreement with those in this study, thus providing an additional confirmation of this phenomenon. After completing its growth in the horizontal direction, Fig. 3(c)(i, ii), the twin begins to grow in width for later times, t = 500 ps, Fig. 3(d)(i, ii). This behavior is correlated to the surface energy anisotropy ratio κ11 /κ22. Moreover, the twin interface thickness has a lower value in the direction normal to the habit plane for the anisotropic surface energy scenario depicted in Fig. 3(c, d)(i, ii) as compared with the isotropic case from Fig. 3(a, b)(i, ii). This phenomenon is related to the contribution of the core and elastic energies to the total surface energy of the interface [185]. The displacement for the anisotropic case in Fig. 3(c, d)(iii, iv) is lower than in the isotropic one in Fig. 3(a, b)(iii, iv). Finally, the variation of shear stresses for anisotropic surface energies at various time instants under small and large strains are depicted in Fig. 3(c, d)(v, vi). Considering the results at t = 500 ps, Fig. 3(c, d)(vi), the maximum and minimum shear stress values for the current simulations are within a 7% difference of the results obtained in [49] by means of a static simulation, demonstrating the significance of inertial terms in extreme loading conditions. For both isotropic and anisotropic surface energies, the magnitude of the shear stress within the twinning region decreases as a function in time and, eventually, becomes negative. This observation is consistent with experimental results for single crystal magnesium under simple shear loading [186].
For the next set of simulation examples in Fig. 4, the same boundary conditions and numerical setup from Fig. 3 are considered for θ = π /6. The layout of the figure is similar to that of Fig. 3 with θ = 0 where (a, b)(i -vi) and (c, d)(i -vi) are the simulation results for the order parameter, displacement, and shear stress under small and large strains at t = 50 ps and t = 500 ps for isotropic and anisotropic surface energies, respectively. For the isotropic surface energy case in Fig. 4(a)(i, ii) at t = 50 ps, the twin is smaller as a consequence of less driving force under the same shear loading of 8% as compared with Fig. 3(a)(i, ii). Further, the twin area fraction at t = 500 ps shown in Fig. 4(b)(i, ii) is much smaller than the case, when the orientation of the habit plane is aligned with the shear loading direction (previously in Fig. 3(a)(i, ii)). In the case of large strains, the twin tends to grow more prominently in the direction of the habit plane when θ = π /6 than when θ = 0 ( Fig. 4(a)(ii)). The displacement contours shown in Fig. 4(a -d)(iii, iv) indicate that the upper and lower sides of the twin's interface are under tensile and compressive loading, respectively, which is similar to Fig. 3(a -d)(iii, iv). The displacement in the vertical direction ( Fig. 4(a)(iv)) is ∼17% greater than that for the small deformation case depicted in Fig. 4(a)(iii), and the maximum shear stress under large strain conditions ( Fig. 4(a)(vi)) is 5% greater than that for the small deformation case (Fig. 4(a)(v)). At t = 500 ps, the twin embryo has a greater thickness for small deformations (Fig. 4(b)(i)) as compared with its growth in length in the direction of the habit plane for the case of large deformations (Fig. 4(b)(ii)), until it is prohibited by the boundaries. The displacement at the end of the simulation is around 17% larger for small strains (Fig. 4(b)(iii)) as compared with the large deformation result (Fig. 4(b)(iv)). Lastly, the spatial variations of shear stress at t = 50 ps and t = 500 ps are depicted in Fig. 4(a, b)(v, vi). As can be seen, the minimum and maximum shear stress values happen in the twinned region and matrix, respectively. The heterogeneous stress distribution around the twins is due to a sudden change in the stresses within the twin interface [187].
Next, the phase-field results for the anisotropic surface energy and θ = π /6 are shown in Fig. 4(c, d)(i-vi). Considering the distribution of the twin order parameter for small strains, Fig. 4(c, d)(i, ii), the twin boundaries tend to be expanded parallel to the habit plane when compared with the isotropic case because the elongation in the direction of s is favored due to a decreasing contribution of the gradient energy term [49]. Pointing to Fig. 4(c)(iii, iv), the maximum displacement values for large deformations are 20% higher than those in the small deformation case from Fig. 4(a)(iii, iv). At the tip of the twin, the shear stress is maximum and ∼10% larger for large strain conditions (Fig. 4(a)(vi)) as compared to the small deformation case (Fig. 4(a)(v)). For the same boundary conditions, the results are depicted for t = 500 ps in Fig. 4(b, d). Here, the twin embryo has a different equilibrium shape than what was shown in Fig. 3 for θ = 0. Namely, the twin is rotated in such a way that one axis in the reference coordinate is aligned to the direction s of twinning shear, as shown in Fig. 4(d)(i, ii). The twin interface also has a lower thickness in the direction normal to the habit plane due to the various contributions of the core and elastic energies to the interface energy [185]. For the displacement contour, the values are 30% larger for the anisotropic energy (Fig. 4(d)(iii, iv)) as compared to the isotropic case, while the difference in shear stress for small and large strains is negligible.

Twin embryo propagation and growth in single crystal boron carbide
For the first time in the literature, the numerical results obtained from phase-field approach are validated with the high-resolution transmission electron microscopy (HRTEM) [188] for the twinning propagation, growth, and interactions in B 4 C. The ubiquitous existence of twins and stacking faults in pressureless sintered and hot-pressed B 4 C, reported in the previous literature [120,189], has motivated studies of their impact [190,191]. This is important because it is widely accepted that existing nanotwins, ranging from t = 1 nm up to t = 30 nm in width for milled and unmilled samples [192], would enhance the strength and hardness of boron carbide [193] by arresting twin boundary slip within the nanotwins [193]. As a result, the presented results opens a number of interesting possibilities for simulating and controlling microstructure pattern development in materials experiencing extreme mechanical loading [194,195]. Given the lack of true images of the twin interfaces in boron carbide [121,196] and the difficulty in experimentally tracking the twin growth process, the present continuum mechanics model will provide insight into the deformation behavior of pre-existed twinned B 4 C, which have been largely neglected in previous works [49,197]. In addition, the morphology of mature twins will be affected by the early stages of the twin nucleus evolution, which necessitates a comprehensive model as herein. In this light, understanding how twins are formed and then developing effective strategies for incorporating twin boundaries into polycrystalline microstructures constitute an attractive approach for enhancing the mechanical response of ceramics. To address this, we conducted numerical simulations using the proposed phase-field model in a boron carbide single crystal.
The combination of growth of a single twin embryo is measured along two critical directions, including twin thickening through twin boundary (TB) migration and twin tip (TT) propagation. The simulation results are then compared with experiments in Fig. 5. Shear strains are applied by displacing all the boundary regions, while the bottom side is fixed. A time step of ∆t = 1 fs is chosen for solving the problem. The dimensions of the simulation domain are 40 nm × 40 nm in the X and Y directions, and contains 160,000 linear triangular elements. One circular twin embryo with a radius of 5 nm is inserted at the center of a square containing the perfect B 4 C crystal lattice, using the Eshelby method as in [198]. The magnitude of applied shear Λ is set to 0.3, which is maximum at the top and zero at the bottom. Additional simulations showed that choosing a shear magnitude lower than 0.3 leads to shrinking and disappearing of the twin.  [193]. (c) and (d) reproduced with permission from [193]. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.) Schematics of the simulation result for an initially circular twin embryo in boron carbide at t = 1 ps and t = 2 ps with the X-axis along the [1101] direction are shown in Figs. 5(a) and 5(b), respectively. Figure 5(a) depicts the "twin tip", which occurs in the primary direction of twin growth, and the "twin boundary", which occurs in a direction perpendicular to the twin growth. Under shear loading, the size and shape of the initial circular twin has changed until reaching a stable configuration. Similar to other ceramics such as calcite, the twin was contracted at the beginning of loading, and this has been shown to be related to the stress reversal [199]. Next, the twin embryo's shape and growth direction at t = 2 ps from Fig. 5(b) is compared with the high-resolution transmission electron microscope images (Fig. 5(c)) and density functional theory results (Fig. 5(d)) [193]. The shape and angle of the twin embryo obtained from the numerical simulations are in good agreement with the previously published results, showing the symmetric twin with an inclination angle of 73.1 • to 73.3 • .
Following this basic validation for boron carbide with results under restrictions of experimental limitations in the literature, the change of the twin size (e.g., length and thickness) and twin interactions in a single crystal boron carbide are explored in order to measure the velocity of twin tips and boundaries (Fig. 6). Being an important parameter for indicating the twin boundary propagation as a key plasticity mechanism, the present findings have important implications for studying the morphology of twins. In order to accomplish this endeavor, the velocities are calculated by tracking the mid points (η = 0.5) on the twin tip and twin boundary interfaces with respect to time. Currently, there is no such statistical data on twin boundary velocity for single crystalline boron carbide, and so we make an attempt to provide some new insights. Considering only one nucleus in the center of the domain, the twin boundary (red colored) and twin tip (blue colored) velocities are shown in Fig. 6(a). The distribution of the twin order parameter at different steps along with the direction for twin tip and twin boundary are also shown in the inset, where the applied shear loading of 0.3 is in the [1101] direction. By choosing ∆t = 1 fs as the time step, the initial circular nucleus shrinks in size until reaching to a stable shape. After that time, the twin starts to grow in the direction of 73.6 • with respect to the loading direction. In this case, the twin tip and twin boundary velocities are larger at the beginning of the loading in comparison with later time instants due to the detwinning process [199] and larger space for unconfined propagation. In addition, the average of twin tip velocities (2.71 ± 0.86 nm /ps) are larger than twin boundaries (2.91 ± 0.37 nm /ps) as a result of having a larger aspect ratio. For the two nuclei scenario shown in Fig. 6(b), the average of twin boundary velocities of the middle embryo (2.76 ± 0.48 nm /ps) are larger than the single twin case because of the tendency of the middle twin to interact with the twin at the top of the inset (termed as Twin #2). The variation of the twin tip velocity is also smaller than the single twin case on the basis of the fast growth of the twin's aspect ratio. Moreover, Twin #2 has a lower aspect ratio, indicating that the two twins will have a wedge shape in the case of interaction between each other. The spreading of a wedge shaped twin has been seen for other ceramics as a result of rapid load drop associated with the twinning process [200]. When placing Twin #3 at the bottom right of the specimen near the fixed boundary conditions (Fig. 6(c)), the average twin boundary velocities of Twins #1 and #2 are increased. This fact is likely a consequence of increasing the twins' aspect ratio, which can be related to the high tendency of twins to interact. Moreover, Twin #3 grows in the direction perpendicular to other embryos because of arresting at the boundary in the scenario depicted for the three twin systems in Fig. 6(c). By adding another embryo close to the fixed boundary condition in a four twin system (Fig. 6(d)), all the twins' aspect ratio has decreased, with Twin #2 by ∼30% in both length and width. Furthermore, the embryo in the middle tended to connect to the nucleus at the top of the domain as a result of the proximity of Twin #2 with the shear loading. Altogether, adding more twin nuclei leads to decreasing the twin boundary velocity of Twin #1, which may be caused by the local stress created from other nuclei to restrict the movement of the boundary.

Fracture-induced twinning in single crystal magnesium
The next example seeks to evaluate the current phase-field approach for studying twinning at a crack tip in magnesium. This simulation is motivated by the urge in understanding the sequence and competition between twinning and fracture, which is difficult to unravel experimentally (e.g., via nanoindentation tests [3]). This is important because far less attention has been given to nucleation and propagation of twins at crack tips and it could offer valuable information on the deformation twinning processes and help to elucidate the role of nanotwinning in crack propagation [201]. In this subsection, a stationary pre-existing crack is considered by a thin notch in a two-dimensional geometry for studying twinning under mode I and mode II cracking. The numerical setup is shown in Fig. 7. An initially square domain of size 100 nm by 100 nm with a pre-existing edge crack of length 50 nm and thickness 4 nm with a rounded tip of radius 2 nm is considered for simulations under a plain strain condition. The crack is assigned a finite radius to alleviate extreme deformations due to singular stress fields at the tip [202]. conditions. The origin of the (X, Y ) coordinate system is at the crack tip, with positive X downward and positive Y to the right.
For boundary conditions, the crack surface is a free surface with a zero Neumann boundary condition. Along each external boundary condition except for the crack surface, the displacements for pure mode I or mode II loading are imposed as in [203]. The orientation of the twin system, s and m, is chosen such that the resolved shear stress is maximum (i.e., θ = 1.2 rad for mode I and θ = 0 rad for mode II). In addition, a small twin nucleus with a radius of 0.8 nm at the crack tip is considered as the initial condition for the twin order parameter. The phase-field results for mode I loading are illustrated in Fig. 8, where a contour of the twin order parameter is plotted. It is clear that the twin growth to the external boundaries is prohibited by the imposed displacement boundary conditions. By progressing in time, the 1 0 1 1 {1 0 1 2} twin band is nucleated at the crack tip and develops at an externally applied strain of 5% due to the stress concentration. The shape and angle of the twin of 69.4 • at t = 110 ps are in agreement with the atomistic simulation results of tensile twinning in single crystal magnesium [204], where a value of 69 • has been reported. The mode II case is shown in Fig. 9 for the twin order parameter at various time instants. Similar to the mode I case, the twin nucleates at the crack tip and starts to grow until it is inhibited by the right boundary condition. As expected, the twin system is aligned in a direction that has the maximum resolved shear stress (θ = 0). These results are in qualitative agreement with the stationary phase-field model under similar boundary conditions [132]. This needle-shaped lenticular twin, which has also been observed in [205], suggests that twin growth occurs by extension of a fast twin tip followed by a coordinated slower migration of the boundaries [206]. The subject of crack growth in the literature has mainly focused on mode I fracture because opening mode crack growth is preferred before that under mixed mode or pure shear mode conditions [207]. It is recognized that, even under pure shear loading, local tensile stresses at the tip result in crack growth under mode I conditions [208]. However, cracks can grow in brittle materials under mode II loading when the ratio between the critical stress intensity factors, KIIc /KIc, is low [209]. It is also motivated that at a sufficiently high confining pressure, the crack is assumed to extend along a smooth curved path that maximizes K II [210]. In heterogeneous brittle solids, the different microstructural inhomogenities (e.g., voids and microcraks) result in a large process regions at the crack tip, and this may lead to macroscopic mode II failure under compressive loads [211]. The study of crack initiation and propagation of mode II fracture is, thus, important in order to better understand the behavior of cracks in brittle solids.
Classically in phase-field modeling in the literature [155], it is assumed that for compressive deformation states, crack growth does not take place. To deal with this, a common technique is to decompose the strain energy density into tensile and compressive parts using a spectral decomposition [87], or a hydrostatic-deviatoric approach [165]; however, both of these decompositions have disadvantages that have yet to be addressed. Specifically, regarding the spectral decomposition, the force-displacement curve shows unphysical stiffening in the fully-cracked specimen [212]. For the hydrostatic-deviatoric method, there are limitations for compression-dominated loading (e.g., the material is allowed to crack in volumetric expansion and shear, but not in volumetric compression) [62]. In addition, both of these popular decompositions can only be used for isotropic materials [99]. Nevertheless, boron carbide has strong anisotropic elasticity ( Emax E min = 8.11, where E max and E min are the general maximum and minimum Young's modulus, respectively) [213]. This analysis is important because the plastic deformation in nanograined boron carbide is assumed to be dominated by intergranular fracture [214] and these new results can be employed toward guiding material design for B 4 C under extreme dynamic loading.

Crack initiation and propagation under biaxial compressive stress in single crystal boron carbide
Consider the biaxial compression test of a single crystal B 4 C specimen with a single pre-existing notch under plain strain condition as shown in Fig. 10. The dimensions of the square domain are those of Fig. 7, and the material parameters are the same as those mentioned in Table 1. Additionally, the fracture surface energy (Υ) and cleavage anisotropy factor (β) are set to 3.27 and 100 J /m 2 (or 0 for isotropic damage), respectively [180]. In the simulations, a total of 323,460 triangular elements were used to discretize the domain with a finer mesh assigned to critical zones. A high confining stress is chosen such that the opening stress intensity factors at the tip of the crack in any direction is zero. The stress parallel to the crack plane is assumed to be larger than the stress values normal to the crack plane (σ XX > σ Y Y ). The initial condition for the time-dependent fracture order parameter and time step are set to ξ(t = 0) = 0.01 and ∆t = 0.5 fs, respectively. In the simulations, all the frictional effects on the crack surfaces are disregarded. The crack evolution process under these numerical conditions is depicted in Fig. 11. As shown, biaxial compression first leads to the initiation of the crack from the tip of the notch (Fig. 11(a)). The range of the fracture order parameter indicates that the crack is not fully formed at t = 0.5 ps. At t = 0.75 ps (Fig. 11(b)), the crack kinks as two single straight branched cracks at a small angle. By progressing in time to t = 0.9 ps, two anti-symmetric cracks begin to propagate toward the top and bottom boundaries due to the larger compressive normal stress parallel with the crack plane ( Fig. 11(c)). In addition, the crack grows in incrementally small steps that are consistent with experimental observations for other brittle materials [215,216]. At the last time frame of t = 1 ps, the propagation path of cracks in single crystal B 4 C is shown (Fig. 11(d)). As can be seen, the crack patterns follow a curvilinear path described by a function ax b . The crack paths reported analytically in [217] and measured experimentally in [215] support this computational result herein. In the studied experiments, b has been found in the interval of 1.43 to 1.58 for pre-fractured specimens of gypsum under uniaxial and biaxial compression [210]. From the analytical model, the exponent was required to be equal to 1.5 in order to be independent of the crack extension length [217]. For boron carbide in this study, the exponent b is obtained as 1.66 ± 0.06 and this is in reasonable agreement with the predicted theory for brittle materials. The curvature parameter a is equal to 0.49 ± 0.08 (nm) −0.66 and the angle of the branched kink is 73.1 • , which is in good agreement with the value (70 • ) reported in [218]. Finally, the homogeneous damage distribution (ξ ∼ 0.5) in Fig. 11 is also due to the hydrostatic nature of the loading. As it is shown in [58], the damage initiation criteria in the current phase-field potential for fracture is fulfilled at infinitesimal load; however, at the crack free surface where the load is not applied, the color is dark blue, which indicates no damage, as expected.

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 thermodynamically consistency of the model ensures that it has a strict relaxational behaviour of the free energy; hence, the models are more than a phenomenological description of an interfacial problem as was done previously in the literature [219]. The dissipation and time scales associated with growth kinetics are also derived and addressed in our paper. The model has been used for studying the evolution of twinning deformation and fracture in anisotropic single crystal magnesium and boron carbide at finite strains. The formulation considers distinct order parameters for fracture and twinning. 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 continuum mechanics, nanometer length scale and picosecond time scale have been used in simulations in this paper.
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: (i) the growth and propagation of deformation twinning in single crystal magnesium and boron carbide, (ii) fracture-induced twinning in single crystal magnesium under pure mode I and mode II loading, and (iii) the prediction of the crack path under biaxial compressive stress loading in single crystal boron carbide. The numerical results for all the problems are in agreement with the available experimental data and analytical solutions in the literature. It has been demonstrated through numerical simulations that the proposed model delivers adequate results matching qualitatively a variety of observed phenomena, including the growth of existing twin embryos, the effect of pre-existing cracks on the twin path under various loading, and the propagation of cracks under compression for highly anisotropic boron carbide. 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 [220]) 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 [221]) and thermally-activated mechanisms (e.g., melting [222]) to capture the behavior of brittle materials in laser spall experiments.

Declaration of Competing Interests
The authors declare no competing financial interests or personal relationships.