Variational regularization of damage models based on the emulated RVE

Material models exhibiting softening effects due to damage or localization share the problem of leading to ill-posed boundary value problems that lead to physically meaningless, mesh-dependent finite element results. It is thus necessary to apply regularization techniques that couple local behavior, described, e.g., by internal variables, at a spatial level. The common way to do this is to take into account higher gradients of the field variables, thus introducing an internal length scale. In this paper, we suggest a different approach to regularization that does not make use of any nonlocal enhancement like the inclusion of higher gradients or integration over local sub-domains nor of any classical viscous effects. Instead we perform an appropriate relaxation of the (condensed) free energy in a time-incremental setting which leads to a modified energy that is coercive and satisfies quasiconvexity in an approximate way. Thus, in every time increment a regular boundary value problem is solved. The proposed approach holds the same advantage as other methods, but with less numerical effort. We start with the theoretical derivation, discuss a rate-independent version of the proposed model and present details of the numerical treatment. Finally, we give finite element results that demonstrate the efficiency of this new approach.

We give an overview of continuum damage models and associated regularization techniques in Sect. 2. In Sect. 3, we will introduce the variational point of view of damage modeling which will allow a specific approach to the ill-posedness of the associated boundary value problems. In Sects. 4 and 5, we will address the lack of coercivity of corresponding potentials to be derived, and in Sect. 6 the failing of those to be quasiconvex. By mending these two deficiencies, we will create our new modeling approach in Sects. 5 and 6, where we will present the novel concept of the "emulated representative volume element." In Sects. 7, 8 and 9 we will demonstrate the behavior of the novel approach via various numerical examples. In Sect. 10, we will introduce a global parameter related to the energy release rate of the system which allows to render the approach rateindependent. Moreover, this will allow to establish a distinct relation to the parameters of Griffith theory and perform a comparison with a model employing gradient-regularization. Finally, and in Sect. 12 we will conclude this paper.
Let us mention that the models derived in this paper are symmetric with respect to tension and compression which is not physical in many situations. However, we felt that addressing this topic would constitute a distraction from the essential message we are trying to convey. So we refrained from it.

Softening phenomena and regularization methods
The perhaps most well-known material behavior that includes softening effects are damage processes: Increasing the external load to a certain threshold value results in local loss of elasticity due to the formation of voids, microcracks and cracks (e.g., [2]). In a global setting, e.g., looking at force-displacement diagrams, this local microstructural evolution leads to a decreasing resultant force as the displacement increases. Consequently, as the strain increases, the stress starts to decrease after a critical threshold value has been reached. Due to the relation ψ = σ dε (one-dimensional case), the Helmholtz free energy ψ initially increases quadratically until it changes its curvature to become a nonconvex function when damage is evolving, and the energy is no longer a coercive function.
The fundamental balance of momentum in continuum mechanics can be derived from different standpoints. One is the minimum condition of the total energy E where denotes the potential of external forces, is the body's volume, b are internal (body) forces (e.g., gravity) and ∂ is the surface of the body on which external tractions t act. The displacement field u that minimizes E is, according to fundamental physical axioms, the one which is realized by the material. The balance of momentum thus equals the minimization of the total energy, subject to appropriate boundary conditions. In case of softening effects, when the Helmholtz free energy ψ is no longer a coercive function, the functional E loses its lower semicontinuity, i.e., it does not possess a minimizer anymore, [6,10,12,13].
Any finite-dimensional discretization of E, e.g., by the finite element method, will of course posses a minimum. This one, however, will only reflect the specifics of the discretization. This is the deeper reason behind the observed mesh-dependence. As a consequence, any method of regularization applied has to replace E by a lower semicontinuous functional that in some sense is sufficiently close to E. There are various strategies that can be used to overcome mesh-dependence. They can usually be grouped into gradient-based, integration-based and viscous ones. Viscous regularization, introduced first by Needleman in [30], limits the evolution of damage in a specific element causing the surrounding elements to also evolve damage in order to maintain the same amount of dissipated energy than in the case of no regularization. Further publications regarding viscous strategies can be found in, e.g., [16,34].
Integral-based approaches, such as in [7,11], postulate a nonlocal measure, which is obtained by averaging over neighborhoods of material points. Information relating undamaged regions to nearby damaged ones is taken into account to achieve mesh-independence. An approach to multiple-scale mesh-free methods that also use an integral-based strategy can be found in [25], whereas an integral-based formulation in an energetic framework was presented in [26].
Gradient-based formulations generally employ a Helmholtz-type equation for the local variable (e.g., see [14,15,32]). Therein, the gradient of the internal variable is explicitly present such that the material behavior is regularized. The drawback is an increased numerical effort. A comparison of nonlocal and gradient-enhanced softening continua was given in [33]. Numerical aspects were discussed in [1] and [27].
For a more elaborated description of the problem of localization and strategies for regularization, we refer to Forest and Lorentz in [17], which gives a very detailed overview. Another classification of regularization techniques is presented in the paper of Yang and Misra [36].
Our aim is to contribute to regularization techniques by proposing a novel approach involving strategies for energy relaxation, which possesses one specific characteristic: There is no additional partial differential equation that has to be solved along with the balance of momentum. This yields a substantial numerical advantage. Only the displacement field is calculated using the finite element method. In comparison with other strategies, in which field functions are used, the number of unknowns at the nodes does not increase. In our strategy, a relatively small adaption for relaxation of the condensed energy yields an approximately quasiconvex energy, which renders the underlying mathematical problem well-posed, although exhibiting stress decrease and localization of the damage parameter. Consequently, our new approach can be demonstrated to give meshindependent results. We present the theoretical derivation, discuss the numerical implementation and give numerical examples.

Variational damage modeling and mesh-dependence
Modeling of the decrease in stress associated with material deterioration can be achieved quite easily by postulating a specific free energy in the form where denotes the free energy of the undamaged material, d is the damage parameter and f (d) is a constitutive function describing the reduction in material stiffness. There exist several opportunities to define f (d), e.g., f (d) = e −d with d ∈ [0, ∞[ and d = 0 as the undamaged material (see [14,15]). A more popular definition is given as f (d) = (1 − d) 2 with d ∈ [0, 1] implying d = 1 for the fully damaged material (see [8,28]). The choice f (d) = 1 − d (see [23,24]), however, results in a thermodynamic driving force which is independent of d yielding an ill-posed evolution equation. The evolution equation for d can be obtained in various ways. A variationally consistent method is the usage of the so-called principle of the minimum of the dissipation function (see [9,20,31]), stating that the rate of d can be determined asḋ for which an appropriate dissipation function (ḋ) has to be defined. According to the assumed material behavior, elasto-viscoplastic, viscous or rate-independent, various approaches for (ḋ) can be made. Here we will focus on a rate-independent formulation of the form taking account of the irreversibility of accumulated damage. Here r is a constitutive parameter that can be interpreted as a yield limit for the initiation of damage. Application of this procedure can be easily performed and the decreasing stress-strain behavior can be modeled without any (mathematical) obstacles for a material point analysis, when the strains are given as a direct function of time. However, things become much more complicated if such a material model is employed in context of a boundary value problem, and implemented, for example, within the context of the finite element method. As discussed above, mesh-dependence is a direct result of ill-posedness for such models. This effect becomes apparent in both the spatial distribution of the damage parameter and the global behavior of the system: The force/displacement curves depart for different finite element discretizations, and even worse, they do not converge to a final result when the mesh size is decreased.

Condensed energy and lack of coercivity
The ill-posedness of continuum damage problems can be best understood and consequently remedied in a time-incremental setting. For this purpose, let us introduce the so-called dissipation distance as where t = t 1 − t 0 ; see [21,29]. This may be interpreted as the energy which has to be dissipated for changing the damage state from d 0 to d 1 within the time interval t. For the dissipation function given in Eq. (7), we obtain Hence, the result is independent of t in this case. Let us consider the potential and the corresponding minimization problem While minimization with respect to the variables u, subject to appropriate boundary conditions, gives elastic equilibrium as defined in Eq. (3), minimization with respect to d returns a time-incremental version of the evolution equation Eq. (6). This second minimization is a point-wise operation giving rise to the so-called condensed energy Subsequently, a purely elastic problem has to be solved for the displacement field: The present formulation has the distinct advantage that the well-posedness of the time-incremental problem can be studied solely in terms of problem defined by Eq. (13), which is well understood (see e.g., [6]). Two conditions are required for the existence of a solution: coercivity, i.e., superlinear growth of the energy as a function of strain, and quasiconvexity, i.e., stability of the energy under superimposed fluctuations of the displacement field. Let us discuss the coercivity of ψ cond d 0 . The stress of the damaged material is given as A fundamental requirement for damage models is the complete deterioration of the material for large strains, i.e., lim where d and ε are related via the stationarity condition where the constant A > 0 depends on d 0 and E 0 only. This, however, cannot be in agreement with the notion of coercivity. A plot of the material behavior due to damage processes is shown in Fig. 1: The stress/strain diagram is given on the left-hand side and the corresponding energy/strain diagram is plotted on the right-hand side.

Coercivity via rate limitation
A straightforward way to achieve coercivity of the condensed energy is to limit the rate at which damage can progress. Rate limitation has been introduced before by Allix, [5]. This can be done by introducing a dissipation function of the form where k is a positive parameter. This yields a dissipation distance Note that this renders the dissipation distance explicitly dependent on t. Hence, our model is rate-dependent now. However, this limitation can be easily justified in the context of high-cycle fatigue. If time is replaced by number of cycles, then our approach limits the amount of damage occurring within an individual cycle, which is a reasonable assumption. Moreover, constant stress relates to constant amplitude of stress in this case. Thus damage evolution has to be expected under these circumstances. The regularizing effect of rate-dependence has been studied in great detail in [35]. The approach to limit damage to two states, here d 0 and d 1 , has been pursued already in the mathematical literature (see [3,4,18,19]). There however, this restriction has been imposed globally, not limited to individual time increments. This resulted in a violation of Eq. (15) giving a nonvanishing limiting stress, i.e., a behavior more remnant of an elastoplastic material.
Let us denote byd(ψ 0 /r ) the unique positive solution of f (d) ψ 0 + r = 0. Then the condensed energy can be written as It can be seen that the graph of ψ cond d 0 consists of two convex parts given by the first (marked green) and third (marked red) alternative in Eq. (20), and a concave part described by the second one (marked blue), in accordance with Figs. 2, 3, 4, and 5. The green line always corresponds to the energy landscape of the previous time step (green part in Eq. (20)), whereas the red curve always corresponds to that of the current time step (red part in Eq. (20)). The mixture energy is indicated by the blue relaxation curve being tangential to both the energy of the previous and current time step (blue part in Eq. (20)). During loading, damage evolves, which, in turn, reduces the stiffness and thus the slope of the respective green and red curves (see Figs. 2, 3, and 4). Consequently, the mixture energy also evolves in time. The damaged microstructure that is realized in the material depends on both the discretized current load (given in terms of the external strain ε(t) and the (fixed) damage state of the previous time step. Since the damage increment is limited, the red curve is also given which allows for a direct computation of ψ cond d 0 . The set of relaxed energy states for all discretized time steps represents the evolving energy landscape and corresponds to the microstructure realized in the material. The condensed energy over time is indicated by the black curve in Fig. 5.
Obviously, ψ cond d 0 still is not quasiconvex. Even after having remedied the lacking coercivity, our damage model is still ill-posed.

Quasiconvexity via relaxation: the emulated representative volume element
The ultimate way to mend lacking quasiconvexity of a potential is to compute its quasiconvex envelope, in our case i.e., the minimum with respect to all superimposed fluctuation fields u = u(x) defined on a representative volume element rep of unit volume (see [12]). While such envelopes have been calculated in [3,4,18], this has been possible for two global damage state only, not for an evolution problem involving many time increments, as it is the case here. Indeed, the initial damage within a time increment is now a field as well, d 0 = d 0 (x). This follows immediately from the stationarity condition Eq. (16). Here, a distributed strain field ε = ε(x) in rep gives rise to a distributed damage field becoming the initial condition for the subsequent time increment. At this point, it turns out to be advantageous not to look at the condensed energy, but to take one step back and consider the cross-quasiconvex envelope of the original energy for fixed damage distribution Now, a corresponding dissipation functional involvingḋ =ḋ(x) is defined aŝ and the evolution of d is given byḋ Calculating Qψ, however, is a formidable task that can only be carried out by numerically solving a high dimensional global optimization problem. Hence, Eq. (24) is not suitable for modeling damage evolution. Therefore, we will resort to calculating the convex envelope, disregarding compatibility of the fluctuation field. Note, however, that the convex envelope is quasiconvex, thus guaranteeing well-posedness of the resulting damage model. Moreover, we will assume a specific damage distribution inside rep . Effectively, we will not calculate what's happening inside rep , but rather model it. We will call the approach the emulated representative volume element, ERVE. In this context, Fig. 6 actually describes what's happening inside rep by the interpretation as a microstructured damage variable and the numerically simple handling by using a fixed number of sub-domains. The development of a microstructured damage variable can be described as follows: When damage, whose increase . This can be interpreted as a microstructured damage variable. However, from a numerical point of view it is more convenient to fix the number of sub-domains to a constant value n for all time steps. This is the condition for the mentioned ERVE (right-hand side). A more detailed view on the ERVE and an exemplary application is given in Fig. 7. On the left-hand side, a simplified case of the microstructured damage distribution inside rep is presented and compared with a corresponding ERVE shown on the right-hand side. In this schematic illustration, the damage structures (left-hand side) are emulated by the sub-domains in the ERVE (right-hand side). The total amount of damage correlates in both representations; consequently, the purple and green structures are related to the damage distribution by the purple and green sub-domains, whereby the increase is again prescribed by exemplary k t = 0.1.
As described, let us assume rep to be divided into n sub-domains of equal volume, each one possessing a constant strain ε i as well as constant damage d i . Instead of minimizing with respect to u, we will do this with respect to all ε i yielding a given overall strain, resulting in Note, that this approach ignores the compatibility of the strain field, thus giving not the true cross-quasiconvex envelope but a lower bound to it. Nevertheless, this lower bound is cross-quasiconvex by itself and ensures the lower semicontinuity of the corresponding functional.
Solving the minimization problem in Eq. (25) gives and The corresponding dissipation function is again defined as and the evolution ofḋ i is given via The time-incremental setting is now canonical. The dissipation distance is given as and the corresponding variational formulation reads Let us introduce driving forces q i defined as Then the minimizers in Eq. (32) satisfy However, the second alternative in Eq. (34) is in our experience never realized during numerical simulations. This is consistent with the fact that these values of d i correspond to concave regions of the energy landscape, compare Fig. 1, which are removed during the relaxation process. A mathematical proof of this would be desirable, of course.
This observation results in the following fairly simple Algorithm 1 on material point level:

Numerical investigation at the material point level
Before switching to more complex simulations, some material point calculations are performed in order to provide a better insight into the model's functioning and behavior. The calculations on material point level are displacement-controlled, three-dimensional, and the chosen material parameters are given in Table 1.
The damage function f (d) = e −d was chosen. The maximum strain ε max = 0.01 is applied in 100 loading steps, and the results are presented in Figs. 8, 9, 10, 11, 12. In the first mentioned result (Fig. 8), the damage variable of each sub-domain is plotted for the last time step of the material point calculation. Hereby, the sub-domains are mirrored so that the damage development can be interpreted as a crack in the microstructure of the ERVE. The crack-like appearance is of course made possible due to the mirroring but also due to behavior of the damage evolution in the sub-domains: Damage always starts to evolve in the first sub-domain, and further sub-domains are filled successively until the total amount of damage over all sub-domains complies with the energetic equilibrium. Therefore, the biggest amount of damage is always located in the middle of the mirrored x-axis.
The just mentioned figure can directly be linked to the next result in Fig. 9 in which the damage variable of each sub-domain is now plotted over the strains. In the first part up to approximately ε = 0.001, no damage evolves due to the elastic behavior. After this point, when the driving force of one or several sub-domains overcome the energetic threshold limit r , damage starts to evolve and increases as long as the energetic equilibrium is not fulfilled. Thereby, as already described before, the sub-domains are filled successively, and through this, the branching structure of the evolving damage arises. At maximum strain, a final damage structure has arisen that complies with the damage development over the sub-domains in the last Fig. 8. Hereby, the top line with the highest amount of damage corresponds to the first sub-domain (see "1"), the last damage-filled sub-domain is the eighteenth one represented by the lowest branched line (see "18"), and finally the nineteenth and twentieth sub-domains that remain unaffected by damage correspond to the lines at the bottom (see "19" and "20").
Applying the chosen damage function f (d) = e −d to the damage variables of each sub-domain provides a plot of these damage functions over the strains (see Fig. 10). According to the last figure, the same elastic and damage parts can be observed. The damage function directly describes the remaining material stiffness and consequently decreases from 1 to 0. Since this plot is only a functional description of the damage plot before, the same branching damage structures are evolving. Therefore, the mentioned sub-domains ("1," "18," "19" and "20") again are illustrated but accordingly in reverse order. Figure 11 combines the remaining material stiffness of all sub-domains into an effective damage function f that also decreases from 1 to 0 and is also plotted over the strains. This effective value is directly used for the calculation of the stresses over strains, pictured in Fig. 12.

Numerical investigation of finite boundary value problems
After providing some material point simulations before, we also want to perform some more complex finite element calculations in order to present the model's global behavior, efficiency and stability. Hereto, the applied stresses and material tangent for the finite element scheme are provided in Appendix A. Three different boundary value problems with plane stress conditions and with, respectively, three different meshes are chosen: the cracked plate (see Fig. 13), the plate with centered hole (see Fig. 14) and the double-notched tensile specimen (see Fig. 15).
The simulations are again displacement-controlled, three-dimensional, and the chosen material parameters are now given in Table 2     The first finite element example is the plate with centered hole with dimensions 1 × 1 × 0.06 mm and three different meshes: 920, 4136 and 9656 structured elements (see Fig. 14). The left-hand side is fixed, and a displacement is applied on the right-hand side in 200 loading steps with an increment of u = 0.05 × 10 −5 mm. The force-displacement diagram in Fig. 19 demonstrates the mesh-independent behavior with an even extraordinary good agreement of the different meshes. It is apparently not possible to differentiate the three curves. Figure 16 provides the plots of the damage variable d that show the fine but smooth Y-shaped evolution of the crack-like evolving damage zones. All meshes provide the same damage distributions and thicknesses of the cracks.
The second finite element example is the double-notched tensile specimen with dimensions 18 × 50 × 3 mm and also three different meshes: 896, 4046 and 10,206 structured elements (see Fig. 15). The bottom is fixed, and a displacement is applied on the top in 200 loading steps with an increment of u = 2.5 × 10 −5 mm. Again, the mesh-independent behavior is demonstrated by the force-displacement diagram in Fig. 20 by extraordinary good agreeing results of the different meshes and it is apparently not possible to differentiate the three curves. As before, Fig. 17 underlines the mesh-independent behavior by equal damage distributions and thicknesses for all three meshes by fine but smooth crack-like evolving damage zones.
The third finite element example is the crack boundary value problem with dimensions 2 × 2 × 0.01 mm and also three different meshes: 935, 4227 and 9727 unstructured elements (see Fig. 13). The bottom is fixed and a displacement is applied on the top in 500 loading steps with an increment of u = 4 × 10 −7 mm. Again, the force-displacement diagram in Fig. 21 demonstrates the mesh-independent behavior of the model since all three curves of the different meshes have a good agreement, even though the coarse mesh with 935 elements differs a little bit. Obviously, the crack singularity has greater demands on the mesh size than the other boundary value problems; the elements at the crack tip are too large in case of the coarse mesh. As before, the results in Fig. 18 underline the mesh-independence of the distributions of the damage variable d that fairly correspond to each other. The mentioned small difference of the coarse mesh is also visible here: At the crack tip, the elements are too large to describe the same sharp damage singularity as illustrated by the finer meshes.
All finite element simulations provide mesh-independent results with extraordinary good agreeing forcedisplacement diagrams, even for coarse meshes. Only the crack boundary value problem needs a specific number of elements at the crack tip in order to capture the crack singularity sufficiently. The plotted damage distributions confirmed these observations: mesh-independent damage zones with equal shape and thickness evolved for all numerical examples.

Numerical investigation of rate-dependence
In this section, the effects of rate-dependence are object of investigation. Hereto, several loading velocities were applied to the boundary value problem of the double-notched tensile specimen (see geometry in Fig.  15). As mentioned before, the bottom of the specimen is fixed and a displacement is applied on the top. The chosen loading velocities areu = {12.50, 8.33, 6.25, 5.00, 2.50, 1.25, 0.25} × 10 −6 mm/s. The simulations were performed with all three given meshes (896, 4046 and 10,206 structured elements, again (see Fig. 15)).
In Fig. 22, the force-displacement diagrams for the simulations with 896, 4046 and 10,206 elements are provided. With decreasing loading velocity, the peak forces are decreasing as well while the drop in the forcedisplacement diagram is getting steeper and more abrupt. This can be interpreted by a more brittle behavior for slow loading velocities in contrast to a more ductile behavior for quick ones, a result of the rate-dependency. The diagrams appear identically for all meshes due to the mesh-independence. Only in case of very slow loading velocities, a slight waviness of the curves can be observed particularly and the coarse mesh with 896 elements slightly differs from the other ones. This indicates that with decreasing loading velocity finer meshes  are required in order to achieve a converged solution. This is a clear indication that mesh-independence is not due to rate effects. Indeed, the required regularization is achieved by relaxation. Additionally, the damage distribution for the loading velocitiesu = {12.50, 8.33, 5.00} × 10 −6 mm/s are shown in Fig. 23, whereby only the finest mesh with 10,206 elements is considered at this point. Interestingly, the three simulations produce quite different results: The faster the loading velocity (left-hand side), the bigger is the reciprocal influence of the propagating damage zones starting from the two notches. The result is one crack connecting both damage origins. In contrast, the slower the loading velocity (right-hand side), the smaller is the reciprocal influence and the two cracks are propagating independently from each other. Between both mentioned cases, there is a crossover (central) that combines both behavior. Moreover, the transition zones between damaged and undamaged material are reducing and are becoming sharper with slower loading velocity. This observation meets expectations because the viscous effects are decreasing with slower loading velocity and the same applies for the transition zones, accordingly.

Rate-independent formulation
A straightforward way to render the present formulation rate-independent is to relate the maximum damage rateḋ max = k to the power of external forces Let us replace the condition 0 ≤ḋ ≤ k in Eq. (18) by In order to check rate-independence, let us for simplicity assume Dirichlet conditions in form of (at least locally in time) linear displacement of the form u = u 0 + vt prescribed on a subset ∂ u . The quantity |v|   corresponds to the loading rate. Finally, it is easy to see that all model equations become independent of |v| by introducing the rescaled time variablet = t/|v|. Please note that this procedure does not work if, instead of restricting the rateḋ, we introduce standard viscosity via adding a quadratic term inḋ to the dissipation potential. Then (ḋ) would not be a homogeneous function inḋ anymore, rendering the rescaling argument above invalid. Removing the first-order term in (ḋ) on the other hand would result in eliminating the damage threshold.
Moreover, ext (u) constitutes a global quantity not related to a material property. This makes the regularizing effect and especially the width of the resulting damage zone dependent on the specific boundary value problem investigated.
The numerical results given below indeed demonstrate the elimination of the rate-dependence. However, let us first examine whether the procedure introduced has an impact on mesh-independence. Figure 24 presents the force-displacement curves and Fig. 25 the related damage distributions for the three different meshes with 896, 4046 and 10,206 elements of the double-notched tensile specimen. The rate-limiting parameter α, which from now serves as the regularization parameter, is set to α = 5 1/(N mm) for all calculations. Obviously, the global structural responses for all three meshes match perfectly.
Since the intended goal at this point is investigating the elimination of rate-dependence, five different loading velocitiesu = {12.50, 8.33, 6.25, 5.00, 2.50} × 10 −6 mm/s are applied for the same problem. The obtained force-displacement curves are compared with each other in Fig. 26 with respect to the mesh with 10,206 elements. The force-displacement curves perfectly match each other, and therefore, the rate-dependence can be assumed to be successfully eliminated. Further observations can be done by comparing the related damage distributions for the three selected loading velocitiesu = {12.50, 5.00, 2.50}×10 −6 mm/s. This is done in Fig.  27. As expected from the global structural responses, also the damage distributions perfectly match for different loading velocities and confirm rate-independence. Once again, let us stress that this result demonstrates the regularizing effect of relaxation, independent of viscous effects or higher gradients.
Lastly, also different strength of regularization can be considered by varying the parameter α. This has been done for different meshes with 896, 4046 and 10,206 elements in Fig. 28, presenting the force-displacement curves for α = {5, 10, 400} 1/(N mm). The related damage distributions are provided by Fig. 29. Obviously, the regularization parameter allows to adjust the damage characteristic for different applications. Thereby, the mesh-independence is preserved, only in case of a very tenuous regularization with α = 5 1/(N mm), resulting in a very local damage behavior, the coarsest mesh slightly differs. However, the medium and fine mesh still have a good agreement.
Besides the aforementioned results showing the double-notched tensile specimen, the rate-independent formulation of the model is also applied to the crack boundary value problem. Since the dimensions of the geometry are different, the parameters are different as well: the loading velocity is now set tou = 2.00 × 10 −7 mm/s and the regularization parameter to α = 4,000,000 1/(N mm). The force-displacement curves for the three different meshes with 935, 4227 and 9727 elements are presented in Fig. 30 and the related damage distributions in Fig. 31. Obviously, the slight regularization allows for a finer crack width, although the coarse mesh becomes mesh independent. However, the medium and fine mesh still have a good agreement. Besides that, the splitting behavior at the end of the crack cannot completely be removed by a slight regularization.
For the purpose of a straight crack behavior, the boundary conditions have to be changed as follows: The bottom and top of the crack boundary value problem are now fixed to all directions, and to illustrate the straight crack behavior, the size of the prescribed crack is reduced. The respective results of the new crack boundary   Fig. 35. Obviously, the curves and damage distributions perfectly match each other and proof the convergence of the numerical calculations.

Comparison with gradient regularization
Let us finally attempt a comparison to a different damage model which is regularized by introducing a gradient of the damage variable, by relating the respective internal lengths with each other. For this purpose, a planar body of thickness w is considered. Then, we can identify an internal length scale h of the rate-independent model given by denotes the energy required for crack initiation, following from Eq. (16). As comparison, we employ the gradient-enhanced damage model presented in [22]. The internal length of this comparative ("comp") damage model is given by In the following, we have adapted the model parameters in order to obtain h = h comp . The following numerical results include two different cases: Firstly, Fig. 36 compares both models at a smaller internal length of h = h comp = 1.83 mm and secondly, Fig. 37 compares both models at a larger Although the characteristic of the damage distribution is quite different, both damage models provide process zones at the crack tip which are quite similar and in the same order of size. For that reason, the internal length scale relates different regularization strategies and damage models and can be set and adjusted in accordance with respective requirements. Moreover, the internal length, or crack width, respectively, scales in accordance with Eq. (39)

Conclusions and outlook
We have presented a novel approach to continuum damage modeling based on a relaxed envelope of the free energy. As expected, this results in a well-posed formulation which gives rise to a robust numerical implementation with mesh-independent behavior. Differently from existing approaches, the presented one does not resort to higher gradients of the field variables and hence does not a priory include an internal length scale. This results in converged force-displacement diagrams even for rather coarse meshes. The presented algorithm is fast and easy to implement.
We derived a rate-dependent and a rate-independent version of the proposed model. The former is well suited to model high-cycle fatigue where time is replaced by the number of cycles which formally leads to ratedependence. The latter requires the calculation of the power of external forces, a global quantity. Employing As a critical point to be mentioned, the results are not very localized in some examples. We suspect this effect is due to the usage of Reuss-type envelopes of the energies. To improve this situation by using Hashin-Shtrikman-type envelopes will be the subject of future work. Another issue to be addressed in the future is the physically required tension-compression asymmetry of damage models. Moreover, the regularization proposed in this work is only shown to be valid for individual time increments. In the future, the behavior over long time periods should be investigated.
Another deficiency of the proposed model is that the rate-independent formulation does not correspond to a true material model, but is dependent on the specific boundary value problem investigated. A possible solution for this problem will be presented in a subsequent publication.
and for the derivative of the stress vector with respect to all individual damage variables in the ERVE as From Eq. (41), we compute the semi-implicit consistent material tangent ∂σ m+1 /∂ε m+1 is given by where (∂ d/∂ε) | m+1 remains to be determined. In the case of, e.g., viscous damage modeling, the missing derivative ∂ d/∂ε can be computed analytically. In our model, however, the update is described by an iterative algebraic scheme (see Algorithm 1). This renders an analytic computation of the derivative intricate. It is therefore convenient to perform a numerical strategy to compute the derivative. To this end, we make use of d = d(ψ 0 (ε)). Application of the chain rule thus yields Consequently, the (numerical) computation of the unknown derivative ∂ d/∂ε transforms to the (numerical) computation of ∂ d/∂ψ 0 . Since the derivative of the vector of damage variables in the ERVE in the current integration point d with respect to the strains reduces to the derivative with respect to the (scalar-valued) free energy ψ 0 , a numerical evaluation is computationally very cheap and is performed by where the tolerance tol = 10 −8 is applied. The derivative of the damage variable d with respect to the strain vector is finally given by In case of no evolving damage, the last term in Equation (44) is equal to zero.