A novel phase-field approach to brittle damage mechanics of gradient metamaterials combining action formalism and history variable

Metamaterials response is generally modeled by generalized continuum based theories. Their inherent substructure leads to a necessity for higher-order theories

thermodynamical modeling of several "phases" within the bulk of material.Hence, often the name "phase-field" is used for the damage parameter and it connects two length scales without any mathematical difficulties [7].Especially in fracture mechanics, the phase-field is applied as a diffusive fracture modeling [8][9][10][11].
The phase-field is also introduced as an order parameter.From a mathematical point of view, the mechanism is based on a potential [12] developed in [13] as a minimization problem or derived by a gamma-convergence as in [14].From a physics point of view, it models microcracks as a bump function such that a differential equation is formally proposed in form of an evolution equation.Loss of bound between molecules, which we call the crack in the microscale, is measured as a jump in the displacement.Balance on singular surfaces is used to bring the phase-field modeling in consistence with configurational forces [15,16].Both viewpoints result in the same parameter, either called a damage or phase-field variable, although historically developed by using different visions [17].
A governing equation-in form of an evolution or balance equation-is necessary for the phase-field.By using a thermodynamical formulation, a balance equation is deduced in [18].Obviously, there are many applications for damage mechanics and the parameters in the formulation are yet to be understood by comparing simulations to experiments or other fine scale simulations by using benchmarks [19,20].Among others, rate dependency [21] and finite speed [22] in crack propagation for brittle materials have been studied.Especially the tip zone and localization of stress [23,24] are investigated with respect to the inherent structure.Hence, they are justified as a size effect [25][26][27].One possible way of modeling this size effect [28,29] relies on strain gradient elasticity or so-called generalized mechanics [30,31], since additional forces on the fracture tip play a significant role in the damage evolution [32].
For a consistent formulation, there are two challenges: formulation and computation.In the formulation, a mathematical chain of reasoning necessitates a postulate of an energy formulation.For first-order (in space derivatives) theories, a model based on a postulated energy is validated by well-established measurements like tensile, compression, or shear tests.Generalization of the energy formulation to gradient metamaterials is indeed possible, even in the case of damage variables, but an experimental validation possibility is not yet available.Especially damage related higher order terms are challenging to comprehend.Therefore, it is ambiguous to discuss whether a postulated energy is correct.
A possible remedy is a straight forward variational approach leading to the weak form.Technically, only the weak form is necessary for computations in generalized damage mechanics.For such a formulation including higher order terms in damage and fracture mechanics, a well-posed methodology is introduced as a non-standard (hemi-)variational approach [50][51][52][53].As a result, governing equations are obtained for the displacement field coupled with Karush-Kuhn-Tucker (KKT) type conditions.Both provide the evolution equations for displacement and damage fields that are assumed as the fundamental kinematical descriptors.The advantage of using this hemi-variational approach is that the KKT conditions and damage evolution equations are derived directly from a dissipation functional, in a systematic manner.Therefore, this formulation yields a robust criterion for loading-unloading-reloading conditions based upon purely mechanical considerations.The damage evolution in such a formulation is an algebraic equation.In fact, a phase-field approach is expected to model more sophisticated systems with crack branching as a contrast to the (hemi-)variational approach.Thus, we aim for a fully variational approach in this work.For the computational implementation, especially in generalized mechanics, different approaches are present [54][55][56][57][58].One necessary condition is a conforming numerical scheme [59][60][61] providing the adequate smoothness in fields.
In this work, we aim at a straightforward variational approach and its computational implementation.Two key issues are discussed for the proposed formulation.First, the definition of the dissipation functional is not given by a Rayleigh function, instead, it is directly incorporated to the action formalism.This procedure opens a possibility for constructing different formulations as well as assures to acquire a minimization problem.Second, the KKT conditions are assumed ab initio, although computationally not problematic, there is no guarantee that the conditions are compatible with the used variational formalism.We briefly explain the formalism and test its feasibility with an implementation utilizing opensource codes for a typical two-dimensional mode I and II crack propagation as well as a three-dimensional a so-called tearing mode (mode III) for gradient metamaterials.

VARIATIONAL FORMULATION
We present a direct formulation leading to the weak form that is necessary for computations.By following [62], we start with the assertion that a Lagrange function, , exists.This function, together with the proper Karush-Kuhn-Tucker conditions (will be defined later), models the system in space-time within a finite domain, Ω, accurately.We will separate space with its infinitesimal element, d, and time, d, composing dΩ = d d.By using the Lagrangean density, , we postulate action: over the domain, Ω, with its first closure, Ω = d d, and second closure, Ω = d d.Subsequently, volume, area, line elements are simply one dimension less quantities giving rise to the second closure of the domain as well.In mechanics, pressure is a typical example of a force applied on an area, analogously, one may realize and model an edge under line forces as well.Herein, the term,   , is neglected, we refer to [63][64][65][66][67] for a more elaborate discussion.We acquire the governing equations as introduced in [68], namely, in the case of a reversible system, the action remains invariant, δ = 0.As the damage is dissipative, we mainly use  for incorporating governing equations related to the damage field and introduce The important distinction is that this formulation is not a Rayleigh functional but more restrictive, hence, a more general approach than herein may be constructed.The reason that we choose the dissipation function analogous to the Lagrange function is visible, as we model the dissipation as a virtual power δ.In this setting, the general formulation, also called the principle of least action, reads In this manner, a minimization problem is acquired and the numerical computation circumvents any convergence problems.We search for adequate  and  in order to obtain the weak form for the phase-field modeling in damage mechanics.
The material system is defined by coordinates, , denoting massive particles.This so-called reference frame is often attached to the initial configuration since the initial placement is known.Furthermore, the reference frame is considered as the ground state that is tantamount to initial conditions with vanishing  and .One of the primitive variables is the displacement, , of particles, .We axiomatically set a displacement function, (, ), with necessary continuity assumptions, i.e., displacement is a primitive variable.For damage mechanics, we introduce two phases and use mixture theory as in [69] in order to define the material as one phase and "void" as another phase.By following [70], we may implement the damage by using a scalar measure, , starting at unity and growing as the damage grows.An analogous approach is to use,  =  −1 , where  starts from 1 and goes to 0 such that it is comprehended as a volume ratio of material in the mixture.We use  = 1 −  as a ratio of the "porosity" such that  = 0 means that the material is fully bonded, no porosity; and  = 1 point out that void "fills in" the unit volume, in other words, the material is broken.Since we use two phases of material and void, the damage variable, , is the phase-field.This phase-field is the second primitive variable with necessary continuity assumptions.Going back to [71][72][73], an auxiliary field is introduced to derive the balance equations in fluid dynamics, we direct the reader to [74] as well as [75] for a brief introduction in this topic.We follow the approach and apply this method to the damage mechanics in order to start from a general formulation, which results in the weak form necessary for computing the primitive variables.

First-order theory
Consider the following densities, , , depending on primitive variables,   , and its first derivatives in space-time,  , , with the index,  = 1, 2, as we incorporate two primitive variables, where a comma notation denotes a derivative (metric is identity) with a Latin index for space (expressed in Cartesian coordinates),  = {, , }, and a Greek index for space-time,  = {, }.We consider this "first" order theory, where the action shall depend on the first derivative.Thus, physically, surfaces, Ω, of the domain, Ω, occupied by the continuum body are closed, Ω = {}.The boundary is smooth and the closure Ω leads to the boundary Ω = Ω∖Ω, where Ω is a Lipschitz domain.Moreover, we assume that the surface energy depends only on the primitive variables; but not on their derivatives,  s =  s (  ( Here and henceforth, we understand Einstein's summation convention over all repeated indices, , , .In time, the discretization will be used in order to calculate in finite time steps.By using a constant time step, time is discretized as the following list: for each of its components, within a time step, the fields are constant.In other words, fields are step-wise functions.Therefore, all the terms with time derivatives are integrated by parts.Since the boundary terms vanish-constant functions within one time step-we obtain the weak form: The continuum body, , possesses the same smoothness as the space-time domain, Ω. Continuum body's boundaries, , are the physical surfaces with the corresponding infinitesimal area element, d.This formulation holds for all systems with the formal assertion that there is a Lagrangean density modeling the system.Defining  and  is arduous and system specific.We give examples in the following section for well-known systems, we emphasize that the proposed  and  functions are models of a physical system.

Damage in brittle materials
We model the first-order theory in damage mechanics, by using the following density functions: with the primitive variables and their derivatives (in space and time), Mass density,  0 , coefficients,  1 ,  2 , and boundary terms, t, q are constants in primitive variables and their rates.They are given functions (in space).We emphasize that displacement, , damage (phase or order) field, , are physically meaningful primitive variables, whereas  is an auxiliary primitive variable with the sole purpose of augmenting the formulation as shown in the following.We refrain ourselves from giving a physical justification for this variable.Even more, we will not solve this variable.The choice of not solving the auxiliary variable is adequate due to the specific choice of  and  functions.From the theoretical point of view,  and  may depend on ,  , , and  • .Nonetheless, we have chosen  and  in Equation ( 8) such that they fail to depend on auxiliary primitive variable or its rates and derivatives.In this way, we skip solving for .We emphasize that the formulation is for one time step, so ,  are the current (unknown) values; their values at the last time step,  0 ,  0 , or even two time steps before,  00 ,  00 , have been calculated-they are known.By inserting density functions (8) into the weak form of the first-order theory in Equation ( 7), we obtain This weak form is applicable to existing formulations for damage elasticity in the literature, as we will demonstrate in the following with an example.By separating the latter as follows: we obtain the first weak form for the displacement, , the second weak form involving the damage field, , and the third weak form depending on the auxilliary variable, .We will make use of the first two for solving  and  in a staggered scheme.We stress that  = (,  , ) and  = (,  , ) such that solving for  is not necessary for obtaining  and .Hence, we discard the third integral form.Until now, the formulation has been general, in other words, universal since material models are not yet incorporated.At the moment of defining  and , we need to model the material response, which is system specific.For moderately large displacements and small strains, we use as usual with a so-called degradation function,  = (), and an adequate strain measure, , called Euler-Lagrange strains.
Especially for the damage mechanism modeled by defining , there are different formulations in the literature [76].We apply a so-called non-standard approach resulting in a diffusive localization band [77,78] by defining  as follows: Fracture toughness,   , is a material parameter, often used as a constant or related to the inner substructure characteristic length in order to model the size effect.The history variable, , is an energy threshold that increases monotonously assuring the irreversibility of the fracture.This variable plays a key role for implementing "non-healing" cracks and it is solely motivated by the phenomenon imposed as in Karush-Kuhn-Tucker conditions.We stress that the history variable fails to be compatible with the variational formulation in the sense that existence of such a threshold may be in contradiction with the "least action" since we synthetically increase its value.However, it is often chosen as a conditional parameter based on the fact that a lost adherence between particles, physically a micro-crack, alters the energy threshold.
We present an alternative derivation by starting with the assertion, where either rate of  vanishes and the damage stops growing; or the stored energy density,  = 1 2       , equals  whenever the damage and hence crack propagates.This definition is based on isotropic hardening in plasticity, where the yield stress (analogous to ) is an increasing function.Hence, in all computations we assume that the crack fails to heal.We circumvent to implement an additional condition for this behavior.For example, again in plasticity, so-called Karush-Kuhn-Tucker conditions are implemented by using a Lagrange multiplier in order to prevent a decrease of accumulated plastic strain.Instead of such a condition, we use a slightly altered formulation,  • | 0 −  0 | = 0, by using the (known) values from the last time step and absolute value of the difference.Now, defining a positive  function by means of absolute value, In computational algebra, the update parameter "∶=" is used to emphasize that the last operation is not a mathematical equation, rather a conditional operation to be implemented into the numerical approach.This non-standard phase-field damage model lacks the boundedness of the phase-field between 0 and 1 such that we will numerically enforce this condition as well.Moreover, the model starts immediately with damage evolution.From the microcracks perspective, this implementation is fine.In some cases, an energy barrier, ē, for the  parameter is introduced to allow damage evolution beyond a threshold value.In other words, damage is prevented for deformation energy below this value, For the stored energy,  0 = 1 2  0     0  , we utilize the known values from the last time step in order to make an implementation possible without additional iterations (as in prediction and correction method often used in the plasticity).The choice of this conditional parameter, , as well as the degradation function, , are of utmost importance in the formulation, as they change the shape and propagation of the crack directly.Especially the specific choice of the degradation function limits the upper bound, which is directly related to the characteristic length, .By using Equation (13), we obtain the following weak forms for  and , On purpose, we circumvent another integration by parts since we search for the weak form.Indeed, it is adequate to perform an integration by parts in space derivatives to address the connection to the boundary conditions, as shown in detail in A. Mechanical loading is modeled by giving a traction t in energy density that is equivalent to force per area.We use q = 0 meaning that the damage field cannot enter or leave the domain across the boundaries, which is reasonable as we know that the crack is not "diffusing" or "migrating" within the material.For time discretization, we utilize the Euler backward time scheme.It is an implicit and thus stable method for problems with real numbers.We acquire the following weak forms for the damage elasticity, ) d . ( This formulation is for a brittle material with cohesive crack propagation derived by ∕ with a crack band described by the degradation function,  = ().

Second-order theory and damage in brittle materials
In order to obtain a second-order theory, we start with the densities, , , depending on primitive variables,   , its first as well as second derivatives in space-time As we consider a "second" order theory, surface and edge energies exist, we assume that they depend on the primitive variables and their first derivatives,  s =  s (  ,  , ),  e =  e (  ,  , ).With the analogous steps, we use the variational formulation and integrate by parts (as many times as necessary) in rate terms leading to the weak form: We stress that the latter weak form is general; however, we need a model of the system, herein we aim at finding a formulation in damage mechanics.Different formulations may be suggested and we propose the simplest model presenting the strength of the approach presented herein.The primitive variables and their derivatives (in space and time) read We model the second-order theory with the densities: where  and  are chosen with the following dependencies:  = (,  , ,  , ),  = (,  , ,  , ).
For the surface energy, we use linear dependency and ignore the edge energy as aforementioned, Based on the last model, we have taken several assumptions: • First, we neglect any additional inertial terms in the Lagrangean density.
• Second, the free energy is extended by the higher gradients of displacement but not of damage.
• Third, only one additional mixed term in  has been added.
• Fourth, as before,  and  do not depend on the auxiliary variable or its space-time gradients.
Varying these assumptions leads to several different formulations, but we have chosen the simplest possible model for generalized damage mechanics.By taking these assumptions, we obtain ) d + ∫  t δ   d = 0 , In the second-order theory, the free energy depends on the first as well as second gradient of displacement.Therefore, the stored energy, , and thus, the free energy, , depend on strain and strain gradients.In the case of an isotropic, centrosymmetric material, we use with additional, in total 7, material parameters, c× , to be determined In the weak form, the following terms arise where we may call them (first order and second order) stresses.For consistency with the first-order theory, we utilize the same relations for ,  1 , and  2 and suggest the simplest possible relation for  3 = −   3 ∕2 such that the formulation for brittle damage strain gradient elasticity becomes We emphasize that we have generalized the formulation to higher gradients by defining  and  independently.Without using such a formulation, the justification of the adequate weak form for the damage variable is obtained phenomenologically.We are not aware of experimental possibilities to directly perform such an analysis.Therefore, the proposed formalism tends to be the feasible approach to generate governing equations for generalized damage mechanics.

COMPUTATIONAL IMPLEMENTATION
In analytical mechanics, all fields belong to  ∞ , i.e., their infinitely many derivatives exist in space.For a numerical implementation, fields need to be discretized in time and space such that their numerical counterparts are finitely differentiable.In time we have already used a stable finite difference method, often called Euler's backward scheme.Hence, the fields are step functions in time, the value during Δ is approximated as constant.In other words, they are  0 in time.
In the case of space discretization as a result of a triangulation, finite element method is used that approximates field functions by   shape functions within every element.We may choose  = 1 for using linear shape functions or  = 2 for utilizing quadratic shape functions within the elements.Across elements, the fields are continuous but not differentiable.This property is of importance for the numerical stability.Since the densities depend on primitive variables and their derivatives, considering the space derivative, we have {  ,  , ,  , , ,  , ,  , } to choose as fields to be approximated.
In the case of the second-order theory, we need to have continuous and differentiable   ,  , , ,  , such that we propose to use a distinct notation for derivatives in order to clearly identify the unknown and the space derivative in the formulation.By introducing ∇  =  , and ∇  =  , , we have the unknowns {  , ∇  , , ∇  }.Hence, in each node, 3 unknowns for   , 9 unknowns for ∇  , 1 unknown for , and 3 unknowns for ∇  are computed.These 3 + 9 + 1 + 3 = 16 degrees of freedom (DOFs) in 3-D or 2 + 4 + 1 + 2 = 9 DOFs are all independent in the formulation.Having derivatives in the formulation, the method is consistent with the usual finite element method, for analogous formulations in this setting, we refer to [25,79,80].Concretely, we choose all unknowns from the same Hilbertian Sobolev space,  1 , spanned by linear shape functions, we refer to [81] for details of the finite element method.Simply stated, we use standard polynomial finite elements constructing a vector space within the computational domain occupied by the continuum body, , with its closure, , and given values on Dirichlet boundaries,  D , and given derivatives on Neumann boundaries,  N .For simplicity, and also as a physical application is rare, we exclude Robin boundaries and define all boundaries either Dirichlet or Neumann type,  D ∩  N = {}.The test functions are chosen from the same space, known as Galerkin procedure, where they vanish on Dirichlet boundaries since the values are given Since we introduce   , ∇  as unknowns and then impose  , = ∇  , a numerical consistency problem may arise by choosing them from the same space.Therefore, we also examine the numerical implementation by using quadratic shape functions for   ,  and linear for ∇  , ∇  with the space In both cases, the weak forms in Equation ( 30) result in ) d . ( In order to enforce the identities, ∇  =  , , ∇  =  , , Lagrange multipliers deliver an action formalism.One possibility is to solve these multipliers, which increases the DOFs greatly.Instead, we use a standard analysis known from the Hu-Washizu principle and by simply comparing the multiplier with the weak from, we obtain that the multipliers are related to the fluxes.For  1  , it is the stress-so we can approximate their variations in {  , ∇  , , ∇  }, as follows: since we set the multiplier by the fluxes, their variations vanish.Especially the computation of the strain and its gradient may be done in different ways.For an increased robustness, we choose from Equation (12), by using the unknowns   and ∇  , as follows: As a consequence, the stored energy is a quadratic function of only first derivatives of unknowns, or in other words, in the weak form in Equation ( 34), the third and fourth terms are quadratic in shape functions, since we use the Galerkin approach and have   and δ   from the same space as well as ∇  and δ∇  .Therefore, the numerical implementation is stable with the ellipticity condition.The complete functional to be minimized reads for each problem of displacement and damage {  , ∇  } = arg min (δΛ 1 + F 1 ), {,∇  } = arg min (δΛ 2 + F 2 ).
The weak forms are given in Equation ( 34) and the multipliers related weak forms in Equation ( 37).This procedure is known as a mixed formulation, where for the classical elasticity it is identical to the Hu-Washizu principle as in [82].By following the computational framework in [83], we solve the problem in Equation (39) by using the open-source FEniCS platform, we refer to [61,[84][85][86] for similar procedures.

SIMULATION
For elasticity and damage problems, different solvers are used from the PETSc library [87] as directly incorporated in FEn-iCS platform [88].For the elasticity problem with the nonlinear weak form in Equation ( 39) 1 , a Newton-Raphson solver is used for the incremental solution of the linearized problem.The linearization is established by the symbolic derivation [89].For each Newton-Raphson iteration, a Krylov solver is utilized based on conjugate gradient (bicgstab) with an algebraic multigrid preconditioner (hypre-amg).For the damage problem with the nonlinear weak form in Equation ( 39) 2 , we use the Scalable Nonlinear Equations Solvers (SNES) component from the PETSc library by utilizing the semi-smooth solver for variational inequalities (vinewtonssls) based on Newton's method [90] using a basic line search.This method allows to use bounds such that the damage parameter is restricted,  ∈ [0; 1].Material parameters for elasticity and strain gradient elasticity need to be chosen.Especially, determination of strain gradient parameters is challenging [91][92][93], there are different homogenization techniques [94][95][96][97][98], specifically for higher order theories [99][100][101][102][103][104] with technical difficulties for strain gradient theory [105,106].An example is the use of an RVE as in [107][108][109], which is problematic because of periodic boundaries with non-constant deformation gradient [110].Therefore, different techniques emerge like applying gamma-convergence [111][112][113] or asymptotic homogenization methods [114][115][116][117][118][119].We apply the granular micromechanical modeling [120][121][122][123] for a porous polymer structure of  Y = 1 GPa Young's mod- for constructing Equation (28).For the damage parameters, we assume the critical energy release rate, as approximated by the mode I fracture toughness   = 300 kPa √ m.This value is realistic for thermosetting polymers such as epoxy.The mobility,  = 2/(Pa s), is chosen based on the study in [124].The geometry and also the diffusive crack length,  = 0.75 mm, are chosen exactly as in [78] for a qualitative comparison of the crack propagation.We emphasize that the choice of the degradation function,  = (), is of importance.One usual choice is Another possible degradation function reads with two additional parameters m and p, we obtain after a straightforward calculation The additional parameters determine the cohesive zone, which models a diffusive crack instead of a sharp jump in the displacement.These parameters are chosen as follows: with a so-called energy barrier, ē.It is challenging to determine such a barrier from experiments, we use a small threshold of ē = 1 kJ/m 3 in the simulation.Both degradation functions are used in simulations and they determine the stiffness change with respect to microcracks.We plot both degradation functions in Equations ( 42), (43) with the chosen parameters in Figure 1.We stress that the parameters,   , , are realistic, m is calculated, and ē, p are simply guessed.Their values affect the numerical convergence as well as the crack propagation.We have chosen all parameters in order to have a consistent set of parameters with the literature, we refer to [125] for a discussion about the sensitivity of these parameters.A two-dimensional plate with symmetric notches as in Figure 2 has been experimentally studied in [126] and numerically proposed to be used with  = 100 mm and  = 100 mm in [127] for higher order theories.The  = 1 mm thick notches are roughly 25 mm long and also rounded with an increased mesh density around it.Bottom of the plate is clamped by using Dirichlet boundary conditions.Equal amount of displacement along  and  on top of the plate is applied by Dirichlet boundary conditions in order to displace the top 45 • tilted and cause a mixed mode crack analogous in the Iosipescu specimen [128].The implementation is done by using Dirichlet boundaries for |  D = (, ) with an amplitude  = 0.01 mm leading to 0.1 mm in  and  directions after 10 s.Therefore, the results are comparable with results in [78] and the asymmetric start of the crack is devoted to the mixed mode loading as demonstrated in Figure 3.
Crack curving and kinking is often devoted to the strain gradient energy; but we stress that the crack propagation depends on the degradation function as well.Not only the speed but also the form and coalescence of the crack is directly altered by the choice of the degradation function as seen in Figure 4. Since the simulation is steered by the position, we point out that the fracture occurs differently in shape and amount.With the degradation function   as in Equation ( 42) the complete fracture is around 60 s meaning 0.6 mm displacement in  and -axis, whereas with   as in Equation ( 43) it is already reached around a third of this displacement.A direct comparison is possible as seen in Figure 5, where we demonstrate color distribution of the displacement magnitude at the same time.Obviously the crack propagation is different with used parameters in Equation (45).Indeed the choice of the degradation function alters the mechanical response as well.
An additional comparison is possible by means of the mechanical response.We refer to Figure 1 for the used   and   .Especially, the shape of fracture is important to notice.In an experiment, force and displacement are measured.The definition of force, , is rather challenging in generalized mechanics.We propose to interpret the free energy from Equation ( 27) integrated over the continuum body, as the work done on the specimen.Hence, the loading on top of the plate has to be equivalent to this energy, which leads to simulation is defined by Dirichlet boundaries, force versus displacement curves are qualitatively identical along  and  axes.Obviously the maximum force depends how quickly the degradation function decreases.As   is steeper for smaller , the strength is smaller.Additionally we compare two different implementations for the degradation function   by using the "linear" space as in Equation ( 31) and the "quadratic" space as in Equation (33).
There is no significant difference in the outcome between these implementations.Also the crack propagation shows no significant differences (not shown herein).Of course, the linear space is computationally more efficient regarding the less amount of degrees of freedom.The proposed approach may be called generalized Castigliano's theorem and constructing such a force-displacement curve is a possible way to calibrate the model by using experimental results.
An analogous geometry of  = 100 mm and  = 100 mm with one notch of  = 1 mm is depicted in Figure 7.The model is a plate with 20 mm thickness.Therefore, the length to thickness ratio of 1 : 5 compels a 3-D model.We utilize the same material parameters as in the 2-D case.As seen from Figure 7, the loading is out of plane and applied on the left and right of the notch in the opposite direction.This loading leads to a tearing (mode III) fracture with a fracture starting around the notch.The boundaries are fully clamped, therefore, also on the boundaries we observe an energy localization that initializes a crack in addition to the notch.The 3-D model has more than 2 million of DOFs for the displacement such that the computational time increase is directly related to the used mesh.We use parallel computing with the mpirun algorithm.The nonlinearity of the problem requires a small time step such that the displacement on the boundary is altered 0.01 mm in each time step.In this manner, the linearized problem is capable of finding the solution by using an iterative solution method as well.We assume that this implementation is novel by using iterative solvers running on several threads in parallel for a damage generalized mechanics problem in 3-D.

CONCLUSION
A variational formalism has been developed for obtaining damage generalized mechanics in a very general form by exploiting auxilliary variables.The suggested formalism entails rate terms with the first derivative in time, which is necessary for the damage problem.Variational formalism involving the second derivative in time is well-established and used for the elasticity problem.Herein, we demonstrate the use of auxilliary variables for damage mechanics.Nevertheless, this method is applicable for multiphysics applications as in thermoelasticity, which is left for further research.Specifically the choice of the Lagrange function allows for a generalization of the damage mechanics for elasticity.This generalization is from the first-order theory to the second-order theory that is equivalent to the strain gradient theory.Such an approach is beneficial, since the dependency of the several terms are difficult to justify, partly because of the ad-hoc implementation of the damage problem.The general formalism is applied for a concrete problem by choosing a Lagrange function leading to a problem formulation with weak forms.They are solved by using the finite element method in space and the finite difference method in time.A mixed formulation has been generated, where displacement and damage variable as well as their space derivatives are computed as independent unknowns.This formulation has the formal benefit of the necessary continuity.Different numerical configurations have been investigated for consistency.No significant difference has been detected in the use of a linear or quadratic solution space.Also two common degradation functions are analyzed and we have observed how the degradation function is altering the crack propagation.For a possible comparison to experimental results, we have proposed a simple method for force versus displacement extraction.All implementation has been done with the aid of open-source packages, for more details of this computational framework, we refer to [83].The codes used for the computation are developed under the GNU Public license [129] made publicly available in [130] for encouraging a transparent and efficient scientific exchange.

A C K N O W L E D G E M E N T S
B. E. Abali's work was partly funded by a grant from the Daimler and Benz Foundation.
Open access funding enabled and organized by Projekt DEAL.

F I G U R E 2
Used geometry and loading condition (left), and mesh (right) with detailed view around the notch end analogous to the Castigliano's theorem in the classical elasticity.The energy calculation is obtained in the simulation and displacement extraction at the top middle point,  top , has been acquired by ParaView Python server.A numerical differentiation by using Numpy packages computes the force.We use the vertical displacement, thus obtain the vertical force component, and show in Figure6acomparative analysis for different configurations.As the F I G U R E 3 Simulation results at 5 s (top left), 10 s (top right), 15 s (bottom left), 20 s (bottom right), deformation is shown (with a scale factor of 10), and the colors denote the phase-field -simulation with the degradation function   as in Equation (43)

F I G U R E 4
Simulation results at 15 s (top left), 30 s (top right), 50 s (bottom left), 60 s (bottom right), deformation is shown (with a scale factor of 5), and the colors denote the phase-field -simulation with the degradation function   as in Equation (42) F I G U R E 5 Simulation results at 20 s, colors denote the displacement magnitude at the initial placement, left with the degradation function   as in Equation (42) and right with the degradation function   as in Equation (43) F I G U R E 6 Vertical force versus vertical displacement on the top end of the double notched specimen in Figure 2 for both degradation functions in Equations (42), (43) as well as different FEM spaces F I G U R E 7 Geometry and loading conditions for 3-D model in top-view (top left), at 30 s the damage field, , in colors and displacement scaled 10 times in perspective view (top right), mesh with decreased ℎ-length around the notch (bottom left), the same damage field in the material frame zoomed in near the notch (bottom right) ).By following a standard variational calculus, we obtain 0  +  00  Δ δ   −     δ  , , },  , = { •  ,  , ,  • ,  , ,  • ,  , ) ,  , = { •  ,  , ,  • , ,  , ,  • ,  , ,  • , ,  , ,  • ,  , ,  • , ,  , ).