How regularization concepts interfere with (quasi-)brittle damage: a comparison based on a unified variational framework

Three regularization concepts are assessed regarding their variational structure and interference with the predicted physics of (quasi-)brittle damage: the fracture energy concept, viscous regularization and micromorphic regularization. They are first introduced in a unified variational framework, depicting how they distinctively evolve from incremental energy minimization. The analysis of a certain time interval of a one-dimensional example is used to show how viscous and micromorphic regularization retains well-posedness within the softening regime. By way of contrast, the fracture energy concept is characterized by ill-posedness—as known from previous non-variational analyses. Numerical examples finally demonstrate the limitations and capabilities of each concept. The ill-posed local fracture energy concept leads by its design to a spatially constant fracture energy—in line with Griffith’s theory. The viscous regularization, in turn, yields a well-posed problem but artificial viscosity can add a bias to unloading and fracture thickness. Furthermore, and even more important, a viscous regularization does not predict a spatially constant fracture energy due to locally heterogeneous loading rates. The well-posed micromorphic regularization is in line with the underlying physics and does not show this undesired dependency. However, it requires the largest numerical efforts, since it is based on a coupled two-field formulation.

Regularization is driven by mathematics and physics applications alike and has been the origin of an impressive amount of studies from various disciplines. We thus find it helpful to first differentiate between possibly distinct and often tacit terminology. Mathematical studies often approach regularity in the sense of singularities and a resulting lack of uniqueness, e.g., by Lipschitz regularization for material softening [6], coercive regularization [7] and for global problems [8]. Tornberg and Engquist [9], for instance, apply regularization to PDEs with singularities. This rather mathematical interpretation of regularization is closely related to regular functions and their properties, e.g., in terms of differentiability. The investigation of mathematical properties can also involve the evaluation of physically sound limits, e.g., for penalty parameters or interface thicknesses [10,11] or the convergence of solution sequences [12,13]. The interested reader is referred to mathematically oriented literature, though, for general conditions on the well-posedness of systems of PDEs, e.g., regarding the domains, compatible initial and boundary conditions, function spaces, and also differentiability and integrability to [10,14].
Engineering research, on the other hand, often introduces regularization to avoid implied problems such as numerical instabilities, branching between multiple solutions and mesh dependence [15][16][17]. Depending on the problem at hand, this is replaced by different criteria. Convexity is one example, cf. [18] for solids undergoing softening due to damage and plasticity. Bourdin et al. [12], for instance, state the strict convexity of the potential energy to be a necessary and sufficient condition for the existence of a unique, smooth crack evolution. The proof is given in the context of Griffith's criterion and smoothness is assumed. Carstensen et al. [19] use the lack of rank-one convexity, implying the lack of quasi-convexity, to demonstrate counterexamples that do not attain a minimum and yield the formation of microstructures in finite-strain plasticity. Ortiz and Repetto [20] study non-convex pseudoelastic energy functions with walls corresponding to single-slip deformations in crystals. Coercivity, growth and monotonicity are often added as requirements to guarantee appropriate behavior at the boundary of physically sound domains of definitions [21,22]. Junker et al. [23] achieve coercivity in viscous regularization via rate limitation. The property of ellipticity-which is equivalent to rank-one convexity under relatively weak assumptions-plays a major role for the well-posedness of systems of partial differential equations [24,25]. The stability of linearized incremental steps or approximative elliptical functionals [12], for example, can thus be easily related to the governing coefficient matrices and eigenvalues thereof. Forest and Lorentz [26] complement the ellipticity condition by appropriate boundary and interface conditions. Material instability is thus prevented as are boundary and interface instabilities. Accordingly, yet other criteria strongly related to well-posedness are stability and branching or bifurcation. Clearly, this is not surprising, since the Legendre-Hadamard ellipticity condition enjoys a long tradition within the analyses of material instabilities in solid mechanics, cf. [27]. Instability can indicate a non-vanishing response for arbitrary small perturbations (and thus a discontinuous dependency on initial or boundary conditions) and branching/bifurcations yield ambiguity (losing uniqueness of the solution) [28,29], violating the last two Hadamard conditions [30,31]. This underlines that uniqueness may not be achievable in the context of well-posedness, let alone physical for many mechanical problems. Pham and Marigo [32] emphasize the influence of geometry and boundary conditions on the stability of gradient-enhanced damage models, while the present focus is on the material stability.

Regularization options and reviewing them within a variational formulation
The variations of regularization methods are numerous and we shall first identify a unifying framework before introducing different regularization methods. The framework of choice in this work is the variational formulation. The regularization methods then follow either by slight modifications of its underlying Helmholtz energy or by the dissipation function. This will allow us to rigorously highlight the differences of the governing equations. The fundamentals of the variational approach go back to Petryk [33], Comi and Perego [34], Ortiz and Repetto [20], Carstensen et al. [19] and also to Miehe et al. [35]. Nowadays, a broad class of dissipative solids are based on the variational formulation [35][36][37]. Even existing models can often be reformulated in terms of an incremental energy, from which the balance laws and the evolution equations follow through stationarity conditions [19,20]. Due to the variational approach, the analysis of the Hadamard conditions can be replaced by the analysis of the convexity properties of the underlying potential. Therefore, the variational approach is a suitable framework to compare the structure of the regularizations on the one hand and to investigate the well-posedness on the other hand [35].
Based on this framework, we will apply known regularization methods to an underlying ill-posed material model for (quasi)-brittle damage. We thus want to sketch three representative concepts. The first one is the so-called fracture energy concept (also known as crack-band theory), cf. [38,39] among others. The material parameters are chosen in dependence of the characteristic length of the underlying mesh, such that each element has the same fracture energy. Clearly, this length also depends on the orientation of the element to the cracking path which lead to the development of advanced calculation methods of the characteristic length [40,41]. By doing so, the loading behavior and the fracture energy are unaffected by mesh refinements; the distribution of the damage variable will still localize in only one element row and therefore remains dependent on the underlying mesh [42].
Viscous regularization constitutes an alternative concept [43]. It is applied to a broad range of materials, e.g., to the modeling of quasi-brittle damage [44] or softening elasto-plasticity [45]. The key idea is the incorporation of a time dependence into the constitutive equations. By choosing that influence to be sufficiently large, the softening behavior becomes mainly dependent on the time step and no longer on the mechanical behavior, i.e., the strains and stresses. Hence, the mechanical tangent remains positive definite as pointed out in [46]. On the one hand, this leads to mesh-independent loading behavior as well as to mesh independent distributions of the damage variable. On the other hand, the results are now additionally influenced by the loading rate, which can be heterogeneous within the structure and which can interfere with additional time-dependent processes.
A third option for regularization are non-local theories [2,3]. The relevant variables are averaged over a certain spatial domain, which exceeds the dimension of the one individual finite element. As a result, a length scale defining the crack width is introduced into the constitutive model. This leads to mesh-independent loading behavior as well as to mesh-independent distributions of local softening variables. Based on this, the gradient regularization has been developed [47]. It results from the non-local theory through a Taylor series expansion up to the second derivative [2,48]. The key idea is to include the gradient of the softening variable into the constitutive framework. Penalization of these gradients hinders localization into one element row. The neighboring elements will also show softening behavior in order to flatten the gradient. The gradient regularization hence introduces a spatial length scale into the constitutive framework. Alternatively, the gradient regularization might be approximated in a micromorphic manner [49,50]. This is often beneficial from an implementation point of view, since it preserves the structure of the underlying local model. It is based on the introduction of an additional global field, which is coupled to the softening variable. Instead of penalizing the gradient of the damage variable directly, the gradient of the global additional field is considered. Similar to the viscous regularization, the micromorphic regularization is applied to a broad range of materials being isotropic brittle damage [50], anisotropic brittle damage [51] crystal plasticity and phase transformation [52], isotropic ductile damage [53,54] and anisotropic ductile damage [55] to mention a few.
There are further concepts for well-posed modeling of damage-induced material softening such as cohesive zone models [56][57][58][59][60] based on a traction-separation law. These models fall into the range of so-called sharp interface frameworks where the fracture zone shows a zero thickness-in contrast to the diffuse or smeared approximations considered within this paper, where the damage zone shows a finite thickness. Smeared approximations of fracture consistent with an underlying sharp interface theory are phase-field models [36,61]. The present focus, nevertheless, is on the explicit regularization of local models as a tool for achieving a well-posed problem description.
While the number of individual regularization methods is quite impressive, a synoptic review is missing. This is probably due to the novelty of recent developments and the challenge of a unified assessment. The present work aims at contributing to this comparison by an analysis of the three aforementioned regularization methods: the fracture energy concept, the viscous regularization and the micromorphic regularization. The comparison will focus on the following key aspects: the variational structure, the well-posedness and the characteristic response (loading behavior, damage field, fracture energy). The approaches are then compared in terms of their numerical feasibility and their mutual compatibility and distinctive limitations. For that reason, Sect. 2 recapitulates the variational formulation in its general form. Based on this method, a local prototype model suitable for modeling isotropic quasi-brittle damage is presented in Sect. 3. In Sect. 4 the prototype model is regularized with respect to the fracture energy concept, the viscous regularization and the micromorphic approximation of the gradient regularization. It concludes with the first key result • a juxtaposition of the mathematical formulations in a unified variational framework.
Secion 5 adds the examination of well-posedness. It closes by the second key result: • highlighting how well-posedness is achieved via the respective coefficient tensors of the incremental step.
The numerical simulations in Sect. 6 finally extend the analysis to complex boundary value problems. The numerical evaluation yields the third key result in the form of the regularization methods' • physics predicted by the regularization and • numerical feasibility The main findings and implications are summarized in Sect. 7.

Variational formulation
The variational derivation of material models will be used as a unifying framework for the analytical comparison of the regularization methods. This method can be traced back, at least, to Petryk [33], Comi and Perego [34], Ortiz and Repetto [20] and also Carstensen et al. [19]. Since then, the topic has been strongly promoted by Francfort and Marigo [36], Ortiz and Stainier [37], Miehe and Lambrecht [35], Petryk [62] and Mosler [63], among others. The key idea is the introduction of potential I and its rateİ. The superposed dot (• ) denotes the material time derivative. After adding powers due to external forces, balance and evolution equations then follow jointly from stationarity. Due to the energetic approach, the convexity properties of the incremental potential I are sufficient to determine the well-posedness of the resulting constitutive equations. In order to recall the nomenclature and structure as briefly as necessary, the rate potential is chosen aṡ Helmholtz energy ψ depends on the linearized strain tensor ε and scalar internal variable κ. It splits into two parts: ψ loc defining solely the local model and ψ mic accounting for the micromorphic gradient regularization. The latter employs additional global field ϕ, cf. [49,50]. It will enter neither the local model nor the viscous regularization. It bears emphasis that the choice of ϕ and ∇ϕ being energetic is not mandatory and actually excludes certain material models, e.g., the phase-field model for fracture. However, these assumptions cover all regularization methods addressed within the manuscript. The state dependent (dependence on κ andκ) dissipation function is split in a similar manner. D loc defines the local model and is homogeneous of order one inκ. D vis is an extension in order to include viscosity as a regularization method. This set of equations allows the representation of all three regularization methods. In what follows, the time-continuous case will first be specified further, followed by the time-discrete counterpart.

Time-continuous setting
According to Petryk [62], Ortiz and Repetto [20], Carstensen et al. [19] and Miehe and Lambrecht [35], among others, a rate potentialĖ can be defined aṡ whereİ is the potential presented in Eq. (1), u is the displacement field, ρ is the density of the material and b and t are forces acting in body B and on its surface ∂B, respectively. The choice of potential (4) is not mandatory. It is also admissible to apply external loads with respect to the micromorphic field ϕ, as it is done for instance in [64]. The total variation of potential (4), given as δĖ = δuĖ · δu + δφĖ δφ + ∂κĖ δκ (5) and forcing stationary, δĖ = 0, leads to the balance laws and the evolution equation as Here, operator ∂ κ E denotes the subdifferential with respect toκ. Starting with balance of linear momentum, i.e., the computation of δuĖ · δu = 0, yields where σ is the Cauchy stress tensor in the framework of linearized kinematics following from the Helmholtz energy as σ = ∂ ε ψ. Through the application of the localization theorem to the weak form of balance of linear momentum (7), it is transformed into its strong form known as By following the same idea, the balance law for the additional global field ϕ is derived as which is also known as the balance of microforces, cf. [65]. Here, ω = ∂ ϕ ψ and = ∂ ∇ϕ ψ follow as derivatives from Helmholtz energy (2). By, once again, applying the localization theorem to weak form (9), the strong form of the balance of micro forces results in The evolution equation of internal variable κ is left to be determined. It follows from stationarity of potential (4) with respect toκ, resulting in Biot's equation, cf. [66], as or in its local format as Here, Y is the energetic dual quantity to internal variable κ and follows from Helmholtz energy (2) as Y = −∂ κ ψ. Evolution equation (12) can be recast in terms of an elastic domain. Its boundary is characterized by δκ I = 0 and internal variable κ evolves, if the boundary of the elastic domain is reached (or exceeded within the trial step). A further restriction to the evolution equation follows from the second law of thermodynamics. The stress power associated with the model is thus leading to dissipation inequality 2.2 Time-discrete setting The time discretized formulation is not only required for the numerical implementation but will also allow the assessment of the well-posedness of the incremental step. The variational framework turns thus into a minimization problem. The time discretization ofu,φ andκ are chosen aṡ The time index n + 1 is omitted for readability. By inserting approximations (15) into potentials (1) and (4) and by subsequently performing a time integration leads to incremental potentials where ψ n = ψ(ε n , ϕ n , κ n ) denotes the Helmholtz energy of the previous time step and where ψ = ψ(ε, ϕ, κ) denotes the Helmholtz energy of the current time step n + 1. In contrast to rate potential (1), incremental potential (16) depends only on u, ϕ, κ and no longer on the ratesu,φ andκ. Similar to Sect. 2.1, the balance and evolution laws follow jointly through stationarity of potential (17), i.e., has to hold. The first term results in balance of linear momentum that is equivalent to time-continuous variation (7). Consequently, the local format of this balance law is equivalent to the strong form of balance of linear momentum (8) at time step t n+1 . A straightforward calculation of the second expression in total variation (18) yields which locally yields balance of microforces (10) at time t n+1 . The remaining variation results in Biot's equation for internal variable κ. It is given in its integral format as and in its local format as The additional constraint due to the second law of thermodynamics is transformed in an analogous manner into a constraint for the evolution equation. It results in

Determining well-posedness for local models
Due to the time-discrete variational framework, an analysis of the convexity properties of the incremental energy E is sufficient to characterize the boundary value problem. Since potential E differs only in terms of linear expressions from potential I, the convexity properties of potential E are identical to the convexity properties of potential I, which allows to shift the analysis of the convexity properties to potential I. Following Miehe and Lambrecht [35], among others, the thermodynamic driving forces σ , ω and are given as (26) and hence, follow also as total derivatives of the incremental energy. The product ∂ κ I d ε κ vanishes, since ∂ κ I = 0 if internal variable κ evolves and otherwise d ε κ = 0. By focusing on local models only, i.e., by neglecting ψ mic and, thus, ϕ, the analysis of well-posedness is reduced to the analysis of the convexity properties of I(u) or I(ε), respectively. In this case, the Hessian follows as Thus, the condition for well-posedness (neglecting degree of freedom ϕ) reduces to

Localized incremental step
Since the analysis in Sect. 2.3 is only applicable to local models, this subsection shows a unified framework in order to determine the well-posedness of constitutive models at specific time steps. The approximation of the evolution equations by an incremental step will later allow an analytical investigation of well-posedness. For that reason, we simplify it further by assuming a localized variation for a spatially homogeneous problem, e.g., for a homogeneous bar under tension and focus solely on the critical time step at which the softening process starts, i.e., when ∂ κ I = 0 is reached for the first time. The incremental step for a quadratic approximation of δI = 0 then reads where L ext accounts for the influence of external loads. The first linearization with respect to time can be written as with symbolξ = (ε,κ,φ) T representing the summary of all variables. Increments are denoted as•, the starting equilibrium configuration is denoted as• and Kelvin notation is assumed for interpretation of the products. Matrix A follows as a partial derivative of the incremental energy as Similar, matrix G contains the partial derivatives with respect to Laplacian ξ = ( ε, κ, φ) T . Since the incremental step relies on a quadratic approximation of I, matrix G follows as Lastly, vector b total is left to be determined. It contains the external loads z ext , the variation evaluated at the last converged time step and boundary terms of ∇ϕ that are typically assumed to vanish at the body's boundaries.
(30) shall be used to analyze whether and how regularization methods provide well-posedness. According to the underlying variational structure, all matrices in 30 are defined by means of derivatives of the potential I. This analysis is limited to the one-dimensional, incremental approximation of the respective models for the sake of analytical accessibility. The simulations in Sect. 6 will then extend the findings to the numerical analysis of complex problems.

Local prototype model
This section introduces and analyzes a local prototype model for isotropic quasi-brittle damage. This model demonstrates the limitations of not regularizing and constitutes the starting point for doing so.

Constitutive equations
Following the variational formulation presented in Sect. 2, the constitutive equations are derived from incremental potential I, which consists of Helmholtz energy ψ and dissipation function D. They are defined as neglecting parts ψ mic and D vis . ψ 0 represents the Helmholtz energy of an undamaged elastic solid and f is a degradation function. f is assumed to be strictly monotonous and bounded between 1 and 0, i.e., f = 1 denotes virgin material points while f → 0 denotes fully damaged material points. A frequently made choice is where λ and μ are the Lamé parameters, κ 0 is the threshold value for damage initiation and κ u and p are parameters defining the fracture energy and the softening response of the material. Insertion of both the Helmholtz energy and the dissipation function into incremental potential (16) results in where subscript n indicates quantities belonging to the previous time step. By way of contrast, variables without this subscript are associated with current time step t n+1 . Since balance laws (19) and (20) are independent of the specific choice of the dissipation function, only the derivation of the evolution equation remains of interest. A straightforward computation of Biot's equation (22) leads to By considering the second law of thermodynamics in incremental format, i.e., the additional constraint κ ≥ κ n has to be complied with. Due to the time-discrete setting, both equations (39) and (40) can be combined into as the final evolution equation. A further analysis of the prototype model is given in "Appendix A" in the form of intermediate calculations and explicit expressions for characteristic values required for parameter identification.

Ill-posedness
Following Sect. 2.4, well-posedness (or the lack thereof) is illustrated for the incremental approximation of an initially homogeneous bar under tension. The stationary conditions for the local prototype model are The solution to this equation is obtained by inverting matrix A pro as long as it is regular. We can check regularity via the eigenvalues pro 1, 2 = fˆ At the beginning, both eigenvalues are real and positive for small ε and κ. A pro thus starts to be regular and the problem starts to be well-posed. A change in sign occurs for one eigenvalue, though, at κ = c E ψ 0 = [κ u /[2 p]] 1/ p and defines the onset of softening-rendering the model ill-posed. The other eigenvalue remains always positive and corresponds to the evolution of κ that remains well defined itself. The aim of the subsequent regularization techniques is to recover well-posedness, including secondary effects such as a mesh-independent solution. Alternatively, the incremental energy can be analyzed according to Sect. 2.3. For the inelastic region, the tangent follows as Thus, as soon as κ > [κ u /[2 p]] 1/ p the model becomes ill-posed. The ill-posedness is highlighted numerically by a one-dimensional example in "Appendix B" with more details on the evolution of eigenvalues and corresponding eigenmodes.

Regularization concepts and selected example models
In view of the ill-posed local prototype model for (quasi-)brittle damage, three regularization concepts are presented and examined analytically: the fracture energy concept, a viscous regularization and a micromorphic gradient regularization.

Fracture energy concept
The so-called fracture energy concept is not a regularization method from a mathematical point of view, since the underlying equations remain ill-posed. In contrast to modifying the governing equations, the fracture energy concepts adapts the solution scheme to fit the physical response. More precisely, the material parameters are chosen in dependence on the underlying finite element discretization. As one major result, the integrated physical quantities, e.g., the load-displacement diagram and the (constant) fracture energy, are matched to experimental observations. The fracture energy concept relies on an immediate localization after leaving the elastic region. Localized inelastic deformation is observed in one element row as far as the finite element method is concerned. A brief analysis of the properties of the prototype model is given in "Appendix A". According to Eqs. (A.8) and (A.9), immediate localization is guaranteed if the material parameters are chosen such that the inequality κ u ≤ 2 κ 0 p holds. Additional constraints are associated with the maximum stress (A.9) and the fracture energy (A.1) and have to be considered with respect to the element's characteristic diameter d e . While the maximum load turns out to be independent of the underlying mesh, the fracture energy G f is reformulated into where g f is the volume specific energy. By means of a structured mesh, the element's characteristic diameter can be approximated well by the element length orthogonal to the crack surface, cf., e.g., [67]. Since, fracture energy G f is given as a combination of material parameters and element characteristic diameters, the material parameters have to be chosen for each element separately. It ensures that each element has the same dissipation per crack surface. This, along with the immediate localization of damage, removes the effect of the numerical discretization on the simulation results.
In summary, the material parameters have to be chosen for each element separately such that (i) the maximum stress (A.9) as well as (ii) the fracture energy (A.1) show identical values for each element and match the experimental observations. The problem description remains ill-posed but secondary effects such as unphysical behavior or mesh dependence shall be avoided by a calibration between material and numerical parameters.

Viscous regularization
The key idea of the viscous regularization is the introduction of a time dependence into the evolution of the internal variable responsible for ill-posedness [43]. A straightforward option is to modify the dissipation function for models based on an incremental energy. Helmholtz energy (2) and dissipation function (3) are then given as While ψ loc and D loc are adopted from the local prototype model, cf. Eqs. (34) and (35), additional dissipation function is the only difference compared to local prototype model 3.1. The quadratic rate term corresponds to the viscosity, whereas a term linear depending on the rate results in a rate-independent model. By choosing an exponent larger than two, the viscous model is extended to power law viscosity, cf. [37]. Since the viscosity is employed here only to regularize the local model, the exponent two has been chosen for the sake of simplicity. Insertion of Eqs. (46) and (47) into incremental energy (16) leads to Since the second integral cannot easily be transformed into an integral over κ, the integral is approximated numerically. By means of an implicit approximation, the incremental energy takes the form The additional time dependence in (50) affects the evolution equation of internal variable κ as where the first part is identical to the time-continuous variation ∂κİ 0. Due to the continuity of κ, the second part converges towards zero when time step t is reduced, cf. [35]. Note that viscosity still affects the evolution of κ under this approximation. A numerical demonstration is given in "Appendix" C. Assuming sufficiently small time steps, the quadratic part of Eq. (51) is thus omitted in the following, so that the evolution equation is approximated by By combining this with constraint (40), the final form of the evolution equation is given as

Micromorphic gradient regularization
The by now classic form of a gradient-enhanced continuum is where the gradient of the internal variable to be regularized enters the constitutive equations. This changes the structure of the underlying material model by shifting the evolution equations to the global/non-local level. Alternatively, a gradient regularization can be approximated in a micromorphic manner. By following Forest [49] and also Dimitrijevic and Hackl [50], the micromorphic part ψ mic of Helmholtz energy (2) is given as κ is coupled to the additional, global field ϕ and the latter gradient is penalized instead. The evolution equation for κ remains at the local level. Starting from the local prototype model, i.e., D = D loc , the micromorphic extension yields an incremental energy of the form Material parameter stress power This extension does not alter the balance equations (Sect. 2) but only the evolution equation. It is given through Biot's Eq. (22) as Once again, by considering the second law of thermodynamics with an associated stress power of the type P = σ :ε + ωφ + · ∇φ, the dissipation inequality follows as Choosing a sufficiently large material parameter c κ enforces the identity κ ≈ ϕ and the dissipation inequality can thus be approximated by means of The dissipation inequality is formally identical to the dissipation of the underlying local prototype model and thus constraint κ ≥ κ n has also to be enforced for the gradient-enhanced model. Evolution Eq. (57) and constraint (58) are combined via Karush-Kuhn-Tucker conditions with indicator function as

Comparison of the regularization approaches
The variational approach allows the comparison of the three regularization concepts in a unified mathematical structure. This is the first intermediate result, see Table 1 for a summarized juxtaposition. By applying a general Helmholtz energy and a general dissipation function, each regularization concept can be derived from a specific combination of material parameters η and c κ . The local prototype model is derived by choosing viscosity η = 0 and penalty parameter c κ = 0, cf. Sect. 3-as well as the model regularized with respect to the fracture energy. The model becomes rate dependent by choosing penalty parameter c κ = 0 and a viscosity η > 0. The internal variable is coupled to additional field ϕ by choosing viscosity η = 0 and penalty parameter c κ sufficiently large. Since the gradient of ϕ enters the constitutive behavior, the gradient of κ also implicitly enters the constitutive behavior.

Well-posedness of the regularization concepts
The regularization concepts introduced in the previous chapter are supposed to restore well-posedness or, at least, diminish effects deriving from ill-posedness. In this chapter, a one-dimensional bar under tension is used to examine their well-posedness. It continues the analysis of the local prototype model from Sect. 3.2. Well-posedness is thus only achieved in the sense of mesh-independent simulation results.

Adaption of the fracture energy concept to ill-posedness
As mentioned previously in Sect. 3.2, the fracture energy concept is still characterized as an ill-posed boundary value problem from a mathematical point of view. The material parameters are chosen in accordance to the element's characteristic diameter, see Sect. 4.1, such that the characteristic observations (here, the fracture energy and the peak load) are identical for all elements.

Viscous regularization
Focusing on an incremental step (Sect. 2.4), the analysis of matrices A and G is sufficient to determine the well-posedness of the constitutive equations. Since the viscous regularization does not introduce a gradient dependence, matrix G vanishes and only stationary condition (30) has to be solved. Clearly, the existence of a unique solution depends on the regularity of matrix A, given for the one-dimensional case as Here, the difference between matrix A pro and A vis is highlighted in gray color rendering the evolution equation The eigenvalue belonging to the evolution equation of κ remains unconditionally positive. The sign of the eigenvalue belonging to ε depends on the time step t. It remains positive as long as By the estimate of 2κ −κ n >κ > κ 1/ p 0 a lower bound for the critical time step can thus be specified as since effective energyψ 0 is limited due to the boundary value problem. As a consequence, the viscous regularization guarantees well-posedness, if the time step is smaller than the critical one. Alternatively, the well-posedness is investigated according to Sect. 2.3. Considering inelastic loading, the tangent is given as which converges for t → 0 to Since a lower bound for the time step is specified in Eq. (65), the viscous example model is thus well-posed. Again, a detailed numerical analysis of the viscous regularization on the basis of the one-dimensional bar is given in "Appendix B". It highlights the numerical evolution of eigenvalues and eigenmodes.

Micromorphic gradient regularization
The coefficient matrices of the quadratically approximated incremental potential (30) read, until reaching the point of ill-posedness, Overbars (•) indicate values of the starting point of the incremental step while a tilde (•) indicates the incremental difference. The differences compared to the underlying prototype model are highlighted in gray color. A trivial solution to the micromorphic problem is exactly the solution of the local prototype model, as long as A mic is regular. More precisely speaking, ξ mic = −A mic −1 · b total with ϕ = κ = c E ψ 0 . For the sake of simplicity, we assume the exponent p = 1 for the local prototype model. Also well-posedness of the local prototype model and its micromorphic regularization are equivalent in this case (det A pro = 0 ⇔ det A mic = 0).
The regularization now becomes effective at the point of ill-posedness. It appears for ε = √ κ u /[E c E ] ⇔ κ = κ u /2. We now focus on how the micromorphic regularization affects the evolution of the solution. We started from an initially homogeneous, uni-axial problem for the sake of illustration. At the onset of illposedness, the system attains ε = σ * /[ f E], κ = 1/2 c E E ε 2 and ϕ = κ. Neither A mic nor G have full rank to allow for a classic solution, though. To find a solution, we have to split the problem. First, it can be shown that the upper part of the equations system applying to ε and κ can be stated as Moreover, the parameter c κ > 0 of the gradient extension causes rank A 2×3 = dim(b 2×1 total ) = dim(ε) + dim(κ) = 2. Given this equation system, there is thus a linear reformulation such thatε = [. . .]φ and κ = [. . .]φ. The latter can be used to reformulate the last line of the first variation of the equilibrium condition as where b total,ϕ denotes the third component of b total . A one-dimensional stress (σ * ) results in the boundary for a homogeneous state. As the state becomes heterogeneous once damage localization appears, extra terms add to γ due to the gradients at the damage region. These extra terms will also be necessary for well-posedness. The role of γ will thus be discussed after we clarified the requirements for the respective solution structure.
As an intermediate result regarding well-posedness, we found that the incremental micromorphic problem can be reduced to a second-order PDE of the auxiliary variableφ. According to the previous linear relationships,ε andκ can be back calculated. Finding a solution, nevertheless, still requires further reformulations and considerations. Notably, the problem description so far is still ill-posed as a substantial condition is yet missing, namely, localization. This becomes obvious by two properties of the mathematical structure. First of all, simply inserting the parameters of the softening state into the pure reformulation, i.e., performing the quadratic incremental approximation just around the softening point, yields γ = 0 ⇒ φ = const. Assuming homogeneous Neumann boundary conditions ∇ ϕ = 0, the solution only reaches ϕ = const. Damage localization as observed in experiments is clearly missing, rendering the solution unphysical. A second indicator for an unphysical solution is the lack of a unique position for damage localization.
To achieve a regularized localization we must further assume two things, in addition to the micromorphic regularization scheme. Firstly, an imperfection to uniquely locate the onset of damage localization and, secondly, a function space large enough to describe its shape. The imperfection can be artificially generated by a heterogeneous initial damage field or geometry, for instance. The description by piecewise solution functions is required mathematically, because a single function cannot solve the governing equation in this example. γ = 0 and ∇ϕ = 0 enforce a single function to remain spatially constant. Only a split into piecewise functions allows a localized damage progression. To be more precise, a split allows for a combination of trigonometric wave solutions with decaying exponential functions and for non-zero gradients at their transition points. This is well explained in the example found in [68], which provides an accessible (yet also extensive) analytical solution of a simpler problem. This split into piecewise functions can also be seen in phase-field problems, e.g., by a discontinuity in the non-zero gradient at the position of damage localization, cf. [61]. Non-academic examples typically do not show this particular need for triggering localization. Applying the regularization scheme suffices, because heterogeneity is usually given in more realistic problems.
Eventually, the piecewise solution for the increment of ϕ can be split into a particular solution and a harmonic one. The particularφ p solution can simply be a constant in this example. The harmonic solutionφ h adds the regularization property. By exploiting its harmonic nature with φ h = λ 2φ h , the incremental problem statement can be written as ⎡ which is well-posed. This can indeed be seen as a new, convex(ified) variational formulation of the harmonic micromorphic problem. In view of the elaborated derivations in [68] for an even simpler study, we omit a further detailed analytic discussion and refer to the illustrative example given therein. The two important conditions to make the problem finally well-posed-in addition to the micromorphic extension-are the introduction of an imperfection and giving up the previous assumption of a perfectly homogeneous problem. A detailed numerical analysis of the (harmonic) evolution of the eigenvalues and eigenmodes associated with the one-dimensional bar is given in "Appendix B".

Direct comparison of the regularization concepts by the localized incremental steps
The following comparison of the governing matrices summarizes the major findings of this section in the form of the adaptions of each regularization concept It is clear that the example models contain the equations of the local models-only the part highlighted in gray is different. This belongs to derivative ∂ κκ I, responsible for the evolution of internal variable κ.

Numerical results
The three regularization concepts are now examined numerically by two illustrative boundary value problems, a pre-cracked plate and an L-shape. The properties of interest cover, among others, mesh objectivity, fracture energy, unloading behavior, crack shape and transferability of model parameters to other boundary value problems.

Boundary value problem and numerical setup
The pre-cracked plate has a height and a width of 500 mm, see Fig. 1. The initial crack is located at the left side with a length of 125 mm. A linearly distributed prescribed displacement (ū ∈ [0, 0.01] mm,u = 0.01 mm/s) is applied to the top and bottom. The assumption of plane stress is made for each calculation. Starting from an initial mesh consisting of 16 × 16 elements, the area of interest is refined recursively. The two meshes after three and four refinements will be used for evaluation, Fig. 7. The model parameters are given in Table 2. The area-specific fracture energy is computed in the post-processing for each element according to Eq. A.1 and summed up orthogonal to the emerging crack.

Results and discussion
The material parameters of the regularized models are chosen so that they match the same peak load and integrated fracture energy, to make comparability as fair as possible. Consequently, the parameters differ from model to model except for the elastic parameters. As a first result, a simple transferability between the models is not available and the regularization of a model involves a subsequent parameter adjustment. Moreover, the fracture energy concept requires mesh-dependent material parameters. Here, only κ u needs to be adjusted for the selected prototype models, cf. Table 2. This is due to the fact that the crack emerges exclusively in the refined area. Apart from that, the loading behavior of the micromorphic model (Fig. 3c) and of the prototype model (Fig. 3a) agrees qualitatively and quantitatively.  Both models show a fast softening behavior up to a load below 2 kN, where the load-displacement diagram shows a kink and slowly converges towards zero, as expected for degradation functions of exponential type. In contrast to the prototype model and the micromorphic example model, the viscous model (Fig. 3b) does not show any change in the softening behavior caused by the exponential function. All three regularization methods eventually yield mesh-objective load-displacement diagrams.
The crack orientation also coincides in all simulations considering the distribution of degradation function f , cf. Fig. 4. However, the damage distribution for the fracture energy concept localizes into one element row. This behavior is neither observed for the viscous example model (Fig. 3b) nor for the micromorphic example model (Fig. 3c). Both show a mesh-independent crack width. However, the viscous example model shows a spatially varying crack width. This is due to the dependence of the model on the strain rate, which varies within the structure.
The fracture energy concept shows a spatially constant fracture energy (Fig. 5). This is expected as the material and numerical parameters have been adjusted accordingly. Also the micromorphic prototype model shows an (almost) constant fracture energy with deviations below 9%. The fracture energy predicted by the viscous prototype model, on the contrary, is far from constant. Also this unphysical influence is caused by the heterogeneous strain rates within the structure, which already influenced the width of the emerging crack.   Fig. 7. Compared to the previous example, two further mesh refinements have been studied due to higher sensitivity regarding mesh dependence. The finest mesh consists of 54532 elements after six refinements. The material parameters have been adopted from Sect. 6.1, cf. Table 2.

Results and discussion
A major difference is the unloading behavior of the viscous model (Fig. 8b). While both the micromorphic example model (Fig. 8c) and its local counterpart (Fig. 8a) show purely linear elastic unloading, this is not the case for the viscous model. This is due to the introduced time dependence of the evolution equation, which leads to further evolution of damage even when the Helmholtz energy-the mechanical driving force of internal variable κ, cf. (53)-decreases. Apart from this, the findings regarding the load-displacement diagram agree with those observed in the previous example. The micromorphic regularization (Fig. 8c) agrees qualitatively with its local counterpart (Fig. 8a). However, the peak loads deviate by approximately 18%. A larger deviation of approximately 20% of the peak load is observed between the viscous and the local example model. The peak load of the viscous model is thereby strongly influenced by external loading rateu. Although being a modest deviation compared to other influences, it indicates the limited transferability of calibrated model parameters from one boundary value problem to the other. A possible explanation lies in the shape of the emerging cracks that may induce additional effects due to their curvature, see Fig. 9.  The regularization methods prove again that they can reliably provide mesh-objective results. Similar to the pre-cracked plate, the fracture energy concept results in a crack band width of one element, highlighting the ill-posed equations. Furthermore, the fracture energy concept yields an (almost) horizontal crack orientation. This is in significant contrast to the viscous regularization and the micromorphic regularization. Both show a curved crack path as observed in experiments, cf. [69]. The viscous model, however, causes a varying crack thickness, indicating again the unphysical dependence on the heterogeneous strain rate within the structure.

Implications for predictions of (quasi-)brittle damage
In summary, all three regularization techniques provide mesh-objective load-displacement diagrams. Deviations for the presented examples can reach moderate amounts, though, such as 20% deviation regarding the peak load.
The damage field is only mesh-objective for the viscous and the micromorphic regularization. Considering the fracture energy concept, the damage width is related to the underlying discretization and requires an adaption of material parameter κ u . Furthermore, the fracture energy concept also shows the least reproducibility of the crack path geometry as it does not accurately capture the curved crack for the L-shape problem, see also [67].
Further influences can be seen in the physical properties, e.g., the fracture energy, which is supposed to be constant for (quasi-)brittle materials. This only holds for the fracture energy concept, since the parameters are specifically chosen for that purpose. It is also almost constant for the micromorphic example model. If fracture energy is a key observation to be made, the viscous regularization shows the least accuracy. This is explained by the spatially heterogeneous strain rate distribution. The artificial viscosity moreover causes an unphysical crack width and inelastic unloading behavior. This may be of increased relevance when modeling fatigue.
Aiming for the least interference with the actual damage physics, the micromorphic regularization seems to be most promising. The numerical effort of the micromorphic example model, however, by far exceeds that of the other two example models. It is based on an additional global field and increases the system of equations by the number of nodes. In terms of the computation time the L-shaped specimen with four refinements times took ∼ 25 minutes for the fracture energy concept and the viscous regularization and ∼ 92 minutes for the micromorphic model. The fracture energy concept and the gradient model, moreover, would allow for further improvement by choosing much larger time steps during unloading. Note that the implementational effort of the micromorphic approach is increased due to the implicit form of its evolution equation. The restrictions of the fracture energy concept and the viscous regularization, on the other hand, require a substantial amount of cumbersome corrections to achieve the reliability of the micromorphic approach.

Conclusions
Three regularization concepts suitable for (quasi-)brittle damage evolution have been presented and compared to each other: the fracture energy concept, the viscous regularization and the micromorphic gradient regularization. In order to make the comparison as fair as possible, they were derived by the method of incremental energy minimization by only extending the incremental energy of an underlying ill-posed prototype model.
By doing so, the fracture energy concept does not modify the underlying prototype model and thus remains ill-posed. Side-effects such as mesh dependence of the integrated structural response such as the load-displacement diagram are avoided by adapting numerical and material parameters. The crack width and path are nevertheless still affected by the mesh size.
The viscous and the micromorphic example models are well-posed. However, they can affect the physics of the prototype model, which has been negligible for the micromorphic approach in the presented examples. The viscous regularization most notably caused an unphysical crack width variation and inelastic unloading behavior. Furthermore, the viscous regularization does not result in a spatially constant fracture energy-as required for (quasi-)brittle material models.
A practical aspect to be considered is the computing time. For the L-shaped specimen discretized with 4558 elements, the computing time increases by a factor of 4 in the micromorphic case. This relates to the additional global degrees of freedom and hence scales by the number of elements. In the viscous case, the computation time is approximately identical to the fracture energy concept, but in general more stable.
In conclusion, the fracture energy concept can only be a valid option if new boundary value problems do not deviate much from successfully calibrated settings and if computation time is relevant. In the ideal case, remeshing of the damage domain is prevented. However, if complex geometries of the crack path are important, the fracture energy concept is not promising, since a mesh bias is observed within the numerical analyses. The viscous regularization may be a viable option for spatially constant or controllable strain rates or if viscosity is naturally given due to the material. In other cases, the micromorphic model remains with the most reliable prediction. Its advantages due to unaltered material response then outweigh the higher implementation and computational effort. The presented analysis aims at supporting such estimates and future developments.
where A c is the area of the emerging crack and g f the volume specific energy. By using identity due to dissipation function D being homogeneous of degree one inκ, specific energy g f follows as A subsequent calculation leads to The expression g 1 f follows as and g 2 f in dependence on material parameter p given in Table 3. For further analysis of the prototype model, the one-dimensional bar is considered. The stress strain relation is given as We further restrict the equations to the purely monotonously increasing case since load reversal leads to purely elastic behavior. Hence, the stress-strain relation can be reformulated as The maximum stress occurs at with its maximum value being (A.9)

Appendix B: Numerical analysis-One-dimensional bar
The present numerical framework highlights the damage evolution by monitoring the eigenvalues and the eigenmodes of the one-dimensional problem, respectively. The one-dimensional bar has a length of 20 mm, is fixed on the left side and loaded on the right side. The material parameters are taken from Table 2 and an imperfection in form of a lowered threshold value κ 0 is applied at coordinate x = 10 mm. The local model shows a change of sign in the lowest eigenvalue at the onset of damage-induced softening, see Fig. 10a, where the eigenmode corresponding to the negative eigenvalue changes from a continuous mode (associated to the previous elastic domain) to a discontinuous mode (Fig. 10b). This mode leads to a change of sign of the displacement field at the imperfection. Hence, it leads to a discontinuous displacement field and to a damage field, which localizes solely in the imperfect element. The viscous example model shows a different behavior. The evolution of the eigenvalues and the corresponding eigenmodes is shown in Fig. 11. In contrast to the local model, no instantaneous changes of the eigenvalues and the corresponding eigenmodes is observed for the viscous example model. Furthermore, the eigenvalues decrease slowly towards zero (but remain positive) (Fig. 11a), such that even at the onset of softening the eigenmodes (Fig. 11b) are identical to those associated to the elastic region. Unfortunately, this is associated with a degradation of all elements of the one-dimensional bar. The eigenmodes develop a discontinuity at the imperfection with further decreasing eigenvalues. These discontinuities do not appear instantaneously (as for the local model), but continuously until they obtain the final form associated with a completely softened state, cf. Fig. 11c. At this state, almost all elements show an identical material degradation. Only the imperfect element is softened a little further, where the difference between these values is dependent on the size of the applied imperfection. The micromorphic example model shows again a different behavior, see Fig. 12. The eigenvalues decrease instantaneously at the onset of softening towards zero (but remain positive), cf. Fig. 12a. The eigenmodes associated with the displacement field at the onset of softening qualitatively agree with those for the elastic region and (more important) do not show any discontinuities (Fig. 12b). Even at a fully softened state Fig. 12c does not show any discontinuities but a consistent localization width. Furthermore, the transition of the eigenmodes from those in Fig. 12b to the ones in Fig. 12c takes place in a continuous manner.
According to the analyses in Sect. 5, well-posedness was investigated in this paper by means of the governing coefficient matrices and their eigenvalues. Within a finite element framework based on a Bubnov-Galerkin ansatz, positive eigenvalues of theses matrices result in positive eigenvalues of the stiffness matrix. Hence, the transition of the smallest eigenvalue of the local model from a positive to a negative sign highlights again the ill-posedness of this model-in contrast to the viscous and the micromorphic relaxed models for which the eigenvalues remain positive.

Appendix C: Difference between the time-continuous and the time-discrete variation
As shown in Sect. 4.2, the evolution equation for the time-continuous case has been further approximated. In order to quantify the difference, both evolution equations-the continuous one ∂κ D ∂κ 0 and the time discrete one ∂ κ D ∂κ 0-have been applied to a one-dimensional example. The error intrinsic to this approximation is numerically analyzed here. The material parameters have been chosen according to Table 2 and the strain rate is set toε = 10 −3 1/s. The results are shown in terms of degradation function f (Fig. 13a) and in terms of the relative difference with respect to the time-continuous variation (Fig. 13b). The black graph provides the reference and corresponds to the time-continuous variation with a chosen time step of t = 10 −4 s. The blue, green and red graphs belong to the time-discrete counterpart and are calculated with different time steps. The evolution of degradation function f (Fig. 13a) only differs marginally. According to Fig. 13b the relative difference decreases with a decreasing time step. For a time step of t = 10 −3 the relative difference is below 1 % and for a time step of t = 10 −3 even below 1 . Thus, the time-discrete evolution equation indeed converges to the time-continuous one.