Variational eigenerosion for rate‐dependent plasticity in concrete modeling at small strain

In the context of eigenfracture scheme, the work at hand introduces a variational eigenerosion approach for inelastic materials. The theory seizes situations where the material accumulates large amounts of plastic deformations. For these cases, the surface energy entering the energy balance equation is rescaled to favor fracture, thus energy minimization delivers automatically the crack‐tracking solution also for inelastic cases. The minimization approach is sound and preserves the mathematical properties necessary for the Γ‐limit proof, thus the existence of (local) minimizers is guaranteed by the Γ‐convergence theory. Although it is not possible to demonstrate that the obtained minimizers are global, satisfactory results are obtained with the local minimizers provided by the method. Furthermore, with the goal of addressing the constitutive behavior of concrete, a Drucker‐Prager viscoplastic consistency model is introduced in the microplane setting. The model delivers a rate‐dependent three‐surface smooth yield function that requires hardening and hardening‐rate parameters. The independent evolution of viscoplasticity in different microplanes induces anisotropy in the mechanical response. The sound performance of the model is illustrated via numerical examples for both rate‐independent and rate‐dependent plasticity.


INTRODUCTION
Numerical modeling of fracture as a limit/asymptotic state of a dissipative system is one of the main challenges in material and structural design. Not necessarily starting from the atomic bond breakage, researchers have been able to isolate this complex process and tried to characterize the phenomena considering the underlying physics. Since the failure of the Liberty Ships in World War II, awareness and understanding of the process have rapidly increased toward development of more durable structures. Nowadays, a major challenge in fracture mechanics is to determine the correct crack-driving force when material nonlinearities are considered. A large number of existing models rely on the variational formulation of fracture presented in References 2 and 3 and fully elaborated in Reference 4, where the free discontinuity is regularized through elliptic functions by taking the sought deformation field (including its strong gradients) in the weak functional space of the special bounded variation (SBV), which is introduced in Reference 5. To this end, the developed approximation is proven to converge toward the original formulation in terms of the Γ-convergence concept presented in Reference 6.
Griffith moved the initial steps in building a sound theory for brittle fracture. In this regard, the fracture process is seen as a competition between the stored bulk energy in a body and the surface energy created from a new crack. As a consequence, this energy balance yields a propagation criterion, where a new crack is localized only if the stored energy exceeds the fracture energy. The irreversibility condition is a priori enforced by the assumption of stress-free crack faces. 7 On the other hand, Barenblatt introduced a cohesive energetic approach to overcome the dimensional (geometrical) inconsistencies arising from the Griffith theory between the volume bulk energy and the surface fracture energy. 8 In this theory, the fracture energy is depending not only on the newly created crack surface but also on the crack opening/displacement jump. As it is shown in Reference 9, the former is an asymptotic condition for the latter if sufficiently large crack openings occur. One of the main disadvantages of these approaches, which is the lack of ability to initiate cracks, is overcome through the variational formulation.
In a celebrated work, 10 Rice was able to characterize the fracture process via a single parameter, namely, the J-integral. It was proven that using a line integral (arbitrary contour), which encloses the yield zone at the crack tip, one can calculate the energy release rate (ERR). This concept led later to the computation of configurational forces (material forces), which, in combination with a node-splitting algorithm, gives rise to a discontinuous thermodynamic consistent method to model crack propagation. [11][12][13][14] Nevertheless, discrete approaches lack to describe ductile fracture, where the solution of the problem is not always a sharp crack.
Furthermore, using the extended finite element method, a crack can be represented within finite elements via the partition of unity method. 15,16 In this way, the near tip asymptotic behavior and the displacement jump representing the crack are described by enriching the standard displacement field with additional local functions. Although promising, the method encounters difficulties in the numerical integration over the elements containing a crack in three-dimensional (3D) applications and in extrapolating the history variables to the nodes.
The focus of the present work is the variational approach of fracture, particularly the version introduced in Reference 1. To overcome the difficulties that arise from the free-discontinuity problem, functions of the SBV space depending on a scale parameter are used to approximate the part of the energy functional that contains the displacement jumps. The variational problem reduces to the minimization of this functional with respect to the field variables, namely, displacements and crack. In Reference 17, the formulation is brought to a more engineering application in form of the eigenerosion approach, implemented into the finite element method (FEM) framework. In the limiting case, for vanishing both mesh size and parameter in the right order, the approximation converges to the Griffith description. Very good results have been achieved for the eigenerosion approach in case of static and dynamic applications (see among others Reference 18). The main advantage of using the eigenerosion approach resides in the description of the fracture process without the need of introducing extra degrees of freedom. In the variational eigenfracture, the crack-tracking problem is solved by only considering the energy balance at the crack tip in postprocessing. Although the regularization via the parameter ensures mesh independence, the binary formulation, that constructs the crack by fully eroding the elements one after another, could lead to spurious mesh bias related to the mesh orientation, and hence, the crack has a preferential direction. The effect was studied in Reference 19, where mesh adaptive refinement techniques have been introduced to circumvent the bias.
The eigenerosion approach to fracture provides mesh independent results, 17 as a consequence of considering a surrounding volume of the crack (the -neighborhood). In this respect, it can be considered embedded with a length scale that avoids the annoying mesh-dependency factors. Material or numerical models that include multiple length scales are in general effective in resolving spurious mesh-dependent behavior. For example, to overcome the degeneracy typical for damage mechanics, numerical models that account for the gradient of the deformation have been proposed. However, they are artificial approaches that require a tuning according to the actual problem at the hand. Taking a more physically correct point of view, a few approaches in the modern mechanics of material literature are based on the modeling of the microstructure, that is, they consider the subcomponents of a material to include one or more suitable length scales, that will automatically be able to describe the degenerative behavior of the material under different loading conditions. [20][21][22] An application of the eigenfracture scheme, in the case of a meshfree discretization, is the material point erosion. The basic principles of this technique follow the standard eigenerosion approach, but in a completely different spatial discretization. Using the max-ent shape functions, 23 the derived robust algorithms can model very complex impact simulations, (see among others Reference 24). In spite of very good results in structural examples, that can capture multispalling at high-dynamic impact, the difficulties arising from continuous update of the shape functions and the boundary treatment face vast challenges as far as computation efficiency is concerned.
A very popular approach, derived from the work of Bourdin et al 2 and Francfort and Marigo, 3 is the phase-field method, which models fracture in a diffuse manner. Following the idea of Mumford and Shah 25 for the discontinuous problem in the theory of image segmentation, convex elliptic functions, that smear the sharp gradients at the crack, are used to regularize the energy functional. This functional is well posed by extending the domain to enforce Dirichlet-type boundaries, where the displacement discontinuity occurs. The regularization parameter herein is called length scale. To have a more realistic energetic formulation in terms of the crack-driving force, the bulk energy is split, to relate only a specific part of the deformation to the crack-driving mechanism. In this way, in Reference 26, a spectral operator split, which differs from tension and compression of the bulk energy, is presented. Moreover, in Reference 27, a deviatoric-volumetric split is used, and, newly, in Reference 28, a directional stress decomposition is presented. In the same way, for vanishing length scale, which requires also a vanishing mesh size, the phase-field method converges to the discrete crack approximation, although, in the work of May et al, 29 it was proven via a one-dimensional (1D) example that Γ-convergence for the phase-field approximation cannot be satisfied numerically.
Although not largely diffused, the literature presents some (mostly) mathematical contributions on treating large gradients in variational terms for the case of inelasticity. In Reference 30, a variational formulation for plastic slip analogously to fracture problems is presented. The energy functional to be minimized is composed of the standard bulk energy and two energies defined in a rectifiable Hausdorff space, where the singular terms are addressed, namely, the surface unpinning and surface plastic energies. The regularization of the free-discontinuity problem via elliptic functions (cf, Reference 31) is proven to be convergent in terms of the Γ-limit. Moreover, in Reference 32, a gradient damage formulation coupled with plasticity is introduced in the variational framework.
Standard gradient damage models are regularly used to model failure. [33][34][35] Even though, from the algorithmic and structural point of view, these models are similar to the variational approaches, the physical background is relatively different. 36 A distinctive drawback of gradient damage models is the lack of the Γ-convergence property. In this way, no minimizers are sought via the damage models, but just an extra boundary value problem is solved. Therefore, it is possible to observe during the numerical simulations the evolution of a wide-damaged zone, which fails to converge to a sharp crack. This contradicts the variational models that converge to the Griffithlike crack for vanishing mesh sizes and regularization lengths.
Our goal is the extension of the existing variational fracture framework to rate-dependent and rate-independent plasticity. In a fully variational setting, considering either plasticity or viscoplasticity, this extension requires the existence of minimizers in the SBV space for the energy functional, which also implies the existence of minimizers for the free-discontinuity problem. This rather nontrivial task requires an enormous mathematical effort. We do not claim to be able to prove the existence, but merely we build the framework of the problem using the established foundations of the variational models and we make sure that our model lies therein. To consider the aforementioned effect when considerable amount of plasticity evolves, the surface energy in the energy balance equation is modified via a definite positive function to decrease the energy necessary to create a new crack. Similarly, in a variational setting, a crack-driving force and a degradation function depending on plasticity were introduced in Reference 37. Furthermore, we introduce herein rate-dependent plasticity to model the material and structural characteristics of concrete. A Drucker-Prager yield surface with caps, described in Reference 35, is modified to account for rate effects using the viscoplastic consistency approach.
The article is organized as follows. In Section 2, the full-mathematical and mechanical minimization problem is constructed. Section 3 gives the formulation of the rate-dependent constitutive behavior of concrete in terms of the microplane approach and introduces the coupling of plasticity to the fracture formulation. A series of numerical examples given in Section 4 shows the capabilities of the rate-dependent and rate-independent fracture plasticity model. Finally, Section 5 closes this work with some useful remarks and open discussions related to variational fracture and constitutive behavior of concrete.

Mechanical preliminaries
Within the infinitesimal strain theory, consider a solid with volume Ω and boundary Ω as given in Figure 1. The Dirichlet and Neumann boundary conditions are prescribed in Ω D and Ω N by fulfilling the conditions Ω D ∩ Ω N = ∅ and Ω D ∪ Ω N = Ω.
In every point of the medium Ω, the balance of linear momentum writes with t and n denoting the surface tractions and outward normal vector, respectively. Displacements are prescribed on Ω D and represented by u. For the clarity of notation used in the following sections, one deals with hard or soft devices, when the boundary conditions Equations (2a) or (2b) are prescribed. The symmetric part of the displacement gradient is the total strain tensor In a plasticly deformed body at small strain theory, the total strain takes the form where the elastic and plastic part of the strain tensor are identified. A set of internal variables is introduced to account for inelasticity and the energy-dissipation functional is written as follows: where the volume bulk energy W(t, u, ) is defined in the medium Ω without the crack Γ and D( , ) denotes the dissipated energy. Moreover, t represents the time, which primarily could introduce in Equation (5) a rate-dependent process and, second, from continuum mechanics perspective, could also refer to a specific configuration in time where a quantity is evaluated. By minimizing Equation (5) with respect to the field variables (preferably globally), one can obtain the necessary constitutive equations. The total energy functional at hand is composed of the bulk and the surface energy. Whereas the latter one increases monotonically in a deformed body undergoing fracture, the former one decreases monotonically, too.

The mathematical problem
A major obstacle in a crack-tracking problem is the disturbance of the continuity due to the localization of cracks in the considered domain. Difficulties arising from the mathematical problem and the need of a predefined crack path require a more general representation of the phenomenon. Thus, we start by introducing the mathematical formulation of the original problem in terms of the variational framework and build the approximative solution using the tools of Γ-convergence theory. Our main focus is an energy functional, which develops discontinuities, at a certain extent similar to the Mumford-Shah function (cf, Reference 25), and its approximation via an eigendeformation field following Reference 1. In Reference 31, the same energy functional is approximated via bounded elliptic functions defined in a specific topological space. This space is introduced to account for weak derivatives in the solution of partial differential equations that lack the differentiability of functions that develop jumps (discontinuities). For the sake of simplicity, in the sequel, we omit any tensorial representation for the considered quantities. In mathematical terms, our medium Ω is bounded by a Lipschitz (regular) boundary Ω. Let n ∈ N be the dimension of the space, Ω ⊂ R n be an open set and Γ ⊂ Ω closed in Ω. Considering the brittle medium given in Figure 1, similar to the representation depicted in Reference 1, one can write where G, a positive scalar parameter in mathematical sense, denotes the ERR, and ℋ n−1 is the n − 1 dimensional finite Hausdorff measure for a set Γ representing for n = 1 the set of discontinuous points, for n = 2 a set of lines, and for n = 3 a set of surfaces in Ω. Functions of the Hausdorff topological space have the property that their value approaches zero if moving away to infinity (far field) from a given point, which greatly helps to smear the local effect of a singularity/crack. Equation (6) is often known as the strong variational evolution, that defines the task of finding the pair (u, Γ), which minimizes the problem The displacement field u is defined in a Sobolev space and, as mentioned before, Γ is a closed set in Ω. Minimization of the total potential energy with respect to the field variables, namely, displacement and a previously assigned fracture parameter, automatically yields the crack path and also crack initiation. In this way, the necessity of a predefined crack to solve the fracture problem is omitted. The minimization in Equation (7) suffers from finding a suitable topology for the pair (u, Γ). To overcome this issue, a new form of the energy functional is represented and defined as the weak form. 25 Minimizers of the weak variational problem are defined in a specific topology, namely, the functional space of SBV, which is part of the space of bounded variations (BV). By exploiting the SBV space characteristics, it is possible to represent the gradient terms of u in a bounded form. In this way, the distributional derivative of u denoted as Du writes which is a finite measure in BV. In Equation (8), D a u represents the absolutely continuous part of the derivative (∇u) with respect to an n dimensional Lebesgue measure, D s u is the singular part of the derivative, where the jumps of the displacement field are considered, and D c u is the diffuse part of the derivative (Cantor). For the SBV space, the Cantor part of the derivative vanishes. This will not be the case for the Barenblattlike functional, where a cohesive energy representation is required. Without further ado, a weak variational form of Equation (6) is given by where S u is a jump set (cracks) for the field u. Representation via Equation (9) keeps the generality with respect to Equation (6) and introduces a well-posed one-field minimization problem, once the Dirichlet boundary is preserved in the SBV space if discontinuities appear in Ω D . In soft devices, where traction boundary conditions (or body forces) are prescribed, the resulting unbounded energy functional fails to ensure global minimality. 4 Bourdin et al 4 enlarge the domain to Ω D ⊃ Ω such that Ω D = Ω D ∩ Ω, where the displacements are prescribed and, hence, enforce a solution only for hard devices. The restrictions to hard-device conditions, where only displacement type loads are permitted, need a special treatment, limiting often the applicability of variational methods in fracture. In Reference 1, the idea of considering the enlarged domain Ω D is enriched by avoiding the evolution of discontinuities in a neighborhood of Ω N , which ensures the existence of minimizers also for soft devices. In the introduced SBV space, using the compactness and lower semicontinuity theorems, Reference 38 prove the existence of minimizers for Mumford-Shah-like functions generally stated in Equation (9). Nevertheless, minimization of Equation (9) is difficult since the position of the singularities S u is unknown. 30 To this end, the regularization of Equation (9) and the solution of the minimization problem in terms of so-called Γ-convergence produces a very powerful tool in fracture mechanics to capture crack evolution.

-convergence
In this section, the application of the Γ-convergence concept is introduced in context of the eigenfracture scheme. For general application of the technique in the field of variational calculus and partial differential equations, the reader is referred to Reference 39. This type of convergence is the most used approach to prove convergence of the variational problems. In this regard, the solution of the minimization problem for the weak variation Equation (9) is carried out by means of Γ-convergence theory. Treating the jumps S u of the discretized form of Equation (9) via the Γ-limit, provides a suitable fracture approach for numerical implementations in a standard FEM framework. The basic idea of this application is that it allows studying limit or asymptotic behavior for very complex problems, for example, Equation (7), instead of finding the exact solution, which faces numerous difficulties due to the presence of discontinuous fields. To this end, as initially introduced in Reference 1, Equation (9) is approximated using the eigendeformation field (ie, Γ) and a scalar parameter , which regularizes the functional with the variables u and Γ in a standard Sobolev space equipped, respectively, with a specific norm to prove the convergence. The term |(S u ) | scaled by 1∕ in Equation (10) is considered the neighborhood volume of the jump S u and, henceforth, will be denoted the -neighborhood. As it will be seen in the following sections, this regularization ensures the calculation of the right amount of the surface energy that will be dissipated after the creation of a new crack (S u ). On the other hand, a limiting case of Equation (10) is written as In Equation (11), one can observe that the eigendeformation field Γ is nothing but the singular part of the distributional derivative Equation (8). Noteworthy to mention, a main characteristic of functions in the SBV space is that they are approximately differentiable almost everywhere, which makes it possible to describe mathematically problems that contain strong discontinuities, see Reference 31 for this notion and its applicability to the Γ-limit. This property is of major interest for the discretized solution of the minimization problem as it lays ground for the existence of these minimizers. Given a family of functions F (u ), the main result of Γ-convergence can be summarized as Equation (12a) is known as the lower bound and compactness property, meanwhile Equation (12b) as the recovery sequence. Equation (12) imply that for every sequence u of minimizers of F , a subsequence u exists that converges to the minimizers u of F. Schmidt et al 1 have proven the Γ-convergence results of Equation (12) for the general case by use of the slicing properties, compactness theorem, and density results of the SBV space (note the analogy between the SBV and the space of special bounded deformations). In general, for these minimizers to exist, the bulk energy (also its sequences when time discretization takes place) should be quadratic definite positive and the surface fracture energy should be lower semicontinuous. The details of the Γ-limit proofs go beyond the scope of this work but we refer the interested reader to References 1,4 and the references therein.
The minimizers of the approximation functions have the characteristics that they take a specific value (bounded, usually equal to 1) close to localization regions, and they vanish in the far field. Consequently, it is possible to represent the strong discontinuities in a mathematical sense. Existence of the Γ-limits ensures the uniqueness of the solution, although it is still questionable if these minimizers are global or local in a time discretized problem. Although global minimizers are favored more than local minimizers, in many applications, it is almost impossible to find them and one should be satisfied with the local minimizers, or in the worst case, with stationary points. In Reference 4, a "backtracking" algorithm is introduced, which corrects itself at each time step by checking if the energy is monotonically increasing. If not, it returns back to the step where it was monotonically increasing, it decreases the time step and searches for the minimizers again. In this way, the authors argue that the global minima are reached with the drawback of being computationally very costly.

Free-discontinuity problem regularization in terms of eigenerosion
The mechanical problem at hand lays within the infinitesimal strain theory for rate-dependent plasticity. In many materials, plastic deformations develop in the body and in singularity regions, they localize to create a crack. To model these phenomena, we couple plasticity, represented by an internal variable, with the fracture process. For this, the surface energy is multiplied by a continuous definite positive function f( ). The motivation behind this rescaling is that the more plastic deformations evolve at a specific point, the more the point is prone to create a crack, and hence, it should facilitate to bypass the energy barrier for crack formation. Under these circumstances, we introduce the regularized energy functional as As mentioned in Reference 39, a remarkable property of Γ-convergence is that, if continuous terms are added to the original problem, its limit convergence (and of its minimizers) remains valid, defining a general framework for different applications without restricting to a specific case. In this way, any new fracture terms added to the energy functional would maintain the generality of the problem. By rescaling of the fracture surface energy with a continuous strictly positive function f( ) would preserve the convexity of the problem with respect to the bulk energy and, at the same time, keep the fracture energy as lower semicontinuous. Under these circumstances, it is possible to find the minimizers and the Γ-limit convergence of these minimizers is ensured.
We adopt here the eigenerosion scheme by assigning two values to the eigendeformation field Γ which correspond either to complete failure or undamaged material. Thus, in the spirit of Griffith and analogously to the derivations in Reference 17, Equation (13) is reduced to a compact form suitable for FEM applications, reading where −ΔE K is the difference of the global energies between the medium with an intact element K and the medium with an eroded element K. In different words, it is the energy that will be released if the element K is to be eroded. ΔA K is the volume of the neighborhood support, which is analogous to the surface of a crack extension. The neighborhood of each element consists of a list of elements whose barycenter is within the radius from the barycenter of element K. If a crack set is denoted as C, see Figure 2, then, its neighborhood is denoted as C . If element K is eroded, the new crack neighborhood becomes (C ∪ K) . In this regard, ΔA K is calculated as: where "∖" denotes the set difference and || • || is the volume of the set (length for 1D, surface for 2D, and volume for 3D). The length of the neighborhood is usually chosen in dependency of the mesh size h as = c 1 h, with c 1 suggested between 3 and 6. −ΔF K in Equation (14) is considered the element energy gain. If this quantity is positive, the element K is set to a list called "priority queue", where the elements are arranged from the highest to the lowest energy gain. Within a user-defined tolerance, the first element(s) is (are) eroded. In a thermodynamic consistent setting, the erosion is modeled by reducing the stresses and stiffness of a specific element, as it will be explained in the next section. The degradation has F I G U R E 2 Discretized representation of eigenfracture for eroding element K a binary form in a sense that an element is either intact, and still carries load for −ΔF K < 0, or, if the fracture criterion is fulfilled, it is eroded and has no more load-bearing capacity. Moreover, the irreversibility constraint is ensured by the second law of thermodynamics which prevents the crack healing. 17

Microplane model
The microplane model is a common approach that is used in FEM to describe the behavior of concretelike materials. The idea follows the fact that a macroscopic quantity (eg, stress, strain) at the material point is decomposed into quantities of lower-order tensor onto some microplanes that surround the integration point in the form of a sphere. Each of these microplanes has different normal directions, which is used to project the aforementioned quantities. The motivation for using a microplane approach to characterize the behavior of concrete is inspired from the fact that, in geomaterials, the evolution of plasticity or cracks is a direction-dependent phenomenon. For that matter, all inelastic constitutive laws can be described on each microplane independently followed by a homogenization on the material point level. As a result, an anisotropic homogenized consistent tangent modulus is obtained. Using the microplane approach, it is possible to decompose the complexity of a macroscopic three-dimensional stress state into microscopic vectorial quantities, where simple constitutive laws can be adopted, see Reference 40 for a thorough description of the approach. Figure 3 shows a conceptual representation of this model. Considering a material point as a sphere, its outer surface can be approximated via different planes, where the macroscopic strain/stress tensor can be projected into vectorial quantities.
We utilize in the current research the microplane model introduced in Reference 41, where the kinematic constraint is used to project the macroscopic strain tensor into its volumetric and deviatoric parts as defined for each microplane and represented by the microscopic scalar and vectorial strains V and D , respectively. The projection tensors used in this operation read with the second-order identity tensor 1 and the symmetric part of the fourth-order identity tensor I sym . Moreover, n is the outward normal unit vector of each microplane. Given a constitutive relation within the kinematic constraint, from the strain tensor, one can derive the stresses and tangent moduli counterparts. As a result, they are splitted into volumetric and deviatoric components. Practically, the normal and shear parts of the stress tensor are used to characterize physical phenomena such as friction, strength, or yield limit, 42 which makes the application to concrete very convenient.
Obviously, the number of microplanes can vary depending on the integration rule chosen, but in the framework used here, a total number of 42 microplanes is considered to be adequate for the homogenization. Exploiting the symmetry properties, only 21 microplanes are sufficient to calculate a homogenized quantity (•) by numerically integrating over these planes using with the weight w mic for each microplane. Note that in the microplane framework, the symbol Ω representing the domain of an integration point is used for convenient notation and it should not be confused with the domain Ω introduced in Section 2.

Rate-dependent Drucker-Prager plasticity
We adopt in this work the plasticity formulation represented by a Drucker-Prager yield surface equipped with a tension and a compression cap formulated in Reference 35 and extend it to rate dependency to account for viscous effects that are observed experimentally in concrete behavior. The stress-strain relation in a deformed body, which undergoes plastic deformation, can be given as follows: where is the homogenized macroscopic stress tensor, vp V and vp D are the microscopic volumetric and deviatoric viscoplastic strains. The transposed deviatoric projection tensor is calculated from Dev T = I dev ⋅ n. The elasticity parameters K mic and G mic on the microplane level are calculated from the macroscopic bulk and shear modulus as K mic = 3K and G mic = G.
The most common viscoplastic formulations are based on the Perzyna and Duvaut-Lions type models, which allow the evolution of stresses also outside the yield surfaces, although it is possible to additionally introduce the consistency condition to the latter one, see Reference 43 for a full treatment of both models. Moreover, introducing rate dependency in the material behavior regularizes the localization of plastic strains, thus leading to mesh-independent formulations. 44,45 A very straightforward and simple to implement approach is the viscoplastic consistency model, 46 which requires that the stress state should be within the yield surface. An application of the consistency model to the Hoffman yield surface is formulated in Reference 47. The crux of this implementation is that the rate effect is introduced into the yield function itself. To this end, the yield surface in the microplane model is constructed as  (20) is a function of the hardening variable and its ratė, the volumetric effective stresses e V and the deviatoric effective stresses e D . This functional dependency makes it possible for the yield surface to expand from both anḋin a so-called "breathable" yield function. Figure 4 shows the graphical representation of the yield surface from a section cut in the space of stress invariants.
In the smooth three-surface yield criterion Equation (20), we make use of the standard Drucker-Prager function f 1 enclosed by a tension cap f t and a compression cap f c , that, additionally to , depend also oṅ. These three surfaces are constructed as In the standard Drucker-Prager yield function f 1 , 0 is the initial yield stress, is the friction coefficient, H and Y are hardening moduli of the linear isotropic hardening laws. The second hardening modulus Y is given as in terms of the initial yield stress and the viscosity material parameter . On the other hand, for the compression and tension caps, the auxiliary quantities X and T are calculated as and the Heaviside functions H c and H t are generally defined as The other material parameters, that appear in the caps, are explained in the following. C V is the volumetric stress of the intersection between the compression cap f c and the Drucker-Prager function f 1 . The same analogy holds for T V . T 0 is the initial intersection of the tensile cap with the horizontal axis. R and R t control the increase of the intersections' absolute values of the two caps with the abscissa axis due to hardening. In Reference 35, a physical interpretation of most of these parameters is explained with respect to standard concrete. To quantify the accumulation of plasticity on the material point level, the following quantity is calculated Considering the multiplication of f 1 with f c and f t , where the activation of the latter ones is controlled by the Heaviside functions, and the construction of the caps, make the overall yield surface a convex continuous function with continuous derivatives. This ensures a unique solution of the stress return-mapping algorithm, given that the convexity is not disturbed due to nonphysical values of the trial state. To this end, in terms of a plastic multiplier , the flow rules are given aṡv From the algorithmic point of view, the treatment of the consistency model does not differ at all from standard elastoplastic models. In this way, for a given stress state, if f mic > 0, the trial state represented by the stresses e,tr V and e,tr D should be projected onto the yield surface via a return mapping algorithm. The solution of the problem is sought by making use of a standard backward Euler scheme to calculate the stresses and internal variables at current time t n+1 as The implicit scheme calculates the stresses on the yield surface depending on the inelastic strains vp,n+1 V and vp,n+1 D at the current time step and updates the hardening variable n+1 as accumulation of the plastic multipliers Δ n+1 for each time step. Combining Equation (27) with the condition f mic = 0 forms a nonlinear system of equations for the unknowns Δ n+1 , e,n+1 V , and e,n+1 D . As for the last remaining constituent, the evolution oḟtakes the simple discrete forṁn The generality of the consistency model is very prominent as it fully recovers the rate-independent plasticity model by only setting the viscosity parameter = 0.0 s. Moreover, due to the fact that, for f mic > 0, the return mapping algorithm projects the trial state onto the yield surface and that Δ < 0 is not admissible, loading/unloading, Kuhn-Tucker, and consistency conditions of the plasticity theory are all automatically fulfilled.

Coupling of eigenerosion and plasticity
It has been observed experimentally for concrete under high-speed loading conditions, for example, impact, that the ERR increases. Consequently, this increase would yield a higher load-bearing capacity for higher rates, which in terms of concrete strength is denoted as the dynamic increase factor (DIF). Nevertheless, although attempts have been made to characterize this phenomenon, a clear relation between the static and dynamic ERR still needs to be formulated. Difficulties arise not only from the understanding of material softening at these loading rates but also from the experimental measurements of the ERR. In Reference 48, the increase of the ERR due to high strain rates is motivated by an extensive increase of microcracking, which might postpone the localization for the same magnitude of loading. The authors therein propose that the dynamic ERR increases proportional to the rate effect of the tensile strength up to a maximum factor of 2.5. To account for this fundamental characteristic, we scale the fracture surface energy with the function f( ) = f( ) in the framework of the constitutive behavior introduced in the previous section (see Equation 14). The function f( ) is defined as where p is a positive scalar depending on the maximum value of plasticity accumulated in all the integration points of an element denoted bŷh om and an algorithmic parameter hom crit such that The critical value hom crit is material dependent and it will define the level where plasticity favors fracture. Recalling again Equation (14), the accumulated elastic energy for each finite element is approximated by the volume integral with the homogenized macroscopic free energy In Equation (32), d and d ′ are given as Using the definition of d ′ , erosion is activated also for the volumetric part of the free energy if volumetric expansion accrues. In a fully thermodynamic consistent setting, erosion is applied to the macroscopic stress tensor and similarly to the effective tangent moduli for all the integration points of an eroded element. We recall the compactness results for the Γ-convergence problem, Equation (12), and note that the energy terms of the proposed model are bounded. The representation in Equation (32) ensures (also in the discretized form) that the energy sequences are quadratic definite positive, which is one of the main components for the proof of lim-inf inequality (Equation 12a). Moreover, during the numerical simulation, it is assumed that these energy sequences are always continuous, and this continuity is not disrupted by local divergence of the return mapping algorithm.
In case of isotropic degradation of the elastic energy, a small artificial stiffness is introduced for the eroded elements, which, from a mathematical standpoint, prevents the problem from becoming ill-posed. This issue is overcome if an energy split leading to an anisotropic energy degradation is considered because it gives the possibility to partially erode the stiffness matrix. Note that the term anisotropic used here refers exclusively to a partial degradation of the energy. The anisotropic fracture mechanism chosen here relates the energy of the tensile volumetric part, the energy of the deviatoric part, and the hardening term given byf (H, ) to the fracture-driving force. Considering the fact that the global effect of the hardening part is relatively small compared with the others, its value is approximated using the standard Drucker-Prager formulation asf (H, ) = 1∕2h(H) 2

Homogeneous uniaxial loading of a unit cube
To show the basic characteristics of the constitutive material behavior presented in the previous section, a unit cube is tested under uniaxial tensile and compression loading in a homogeneous test. Initially, cyclic loading is applied to the unit cube in form of a displacement-controlled simulation in tension and compression for a number of two cycles. The geometry and one loading cycle are shown in Figure 5. The applied load takes the form û = nf (t)ū, with the number of cycles n = 2, time-dependent loading function f(t) given in Figure 5B, and the prescribed displacement û = 1 mm. The duration of one cycle is t 0 = 1.33 × 10 −4 s. In Reference 35, a material parameter study for concrete was presented. These parameters are adopted also in this work to highlight the rate effects of the overall response. Based on the characteristics of normal concrete, from the uniaxial compression strength f uc , several other parameters can be derived such as biaxial compression strength f bc = 1.15f uc , uniaxial tensile strength f ut = 1.4(f uc ∕10) 2/3 , friction angle = ( Clearly, these suggestions should not always be considered a priori, but consulting also with any test data, their choice, hence the shape of the yield function, should deliver stable numerical simulations and reliable results. Table 1 shows the chosen parameters for this example. Four different cases are considered, where it is started with the elastic case, continued with the elastoplastic case by neglecting rate effects simply with = 0.0 s, and finally, two rate-dependent viscoplastic results with an arbitrary = 10.0 s and strain rateṡ1 = 100∕s, 2 = 300∕s. The stress-strain relations in loading direction for the homogeneous tests are shown in Figure 6. Note that for the current example, no fracture criterion is activated and the specimen cannot fail. Observing Figure 6, no inelastic strains are present for the elastic case and symmetric behavior in tension and compression is achieved. On the other hand, the effects of the two caps in the yield function (Equation 20) are reflected in the elastoplastic behavior with hardening, represented by the dashed curve. It can be seen that the amount of evolved compressive stresses is higher than the tensile ones. This follows due to the larger size of the compression cap of the yield surface, shown in Figure 4. Moreover, when the rate effects are considered, more hardening occurs and less inelastic strains develop. Similar to elastoplasticity, the asymmetry between compression and tension is evident also in this case. One can also see the response at different rates, where foṙ2 >̇1, more hardening develops in the material. In general, it can be stated that at specific loading conditions, the rate-dependent behavior always lies between the elastic and rate-independent elastoplastic response.

Four-point bending test of a beam with notch
To check the behavior of the model for a structural example, a concrete beam with a notch at its midspan is investigated at standard four-point-bending conditions. Its geometry and boundary conditions are illustrated in Figure 7A, where the shown dimensions are given in millimeters. The finite element discretization is depicted in Figure 7B. The beam has a depth of 50 mm. A total number of 3680 low-order eight-node solid elements are used for meshing, with finer elements in the region where the crack is expected to propagate. The experimental details for this setup can be found in Reference 49. The simulation is run under a displacement-controlled regime to avoid abrupt failure of the specimen. Table 2 shows the chosen material parameters. The critical ERR takes the value G c = 5.5 × 10 −4 kN/mm and the parameter controlling the effect of plasticity on the evolution of eigenerosion is chosen to hom crit = 5.0 × 10 −5 1/MPa. Note that the used G c should be seen as an effective value during the simulations, which might differ from the experimental measurements due to the nature of how the numerical simulation is constructed and the formulation of the eigenerosion approach. Moreover, due to the scaling function used in the current implementation, the initial value of G c is decreased with the evolution of plasticity, which favors fracture over indefinite plastic deformation of concrete. As quasi-static-loading is applied, rate effects are not considered. As it was mentioned in the previous section, the rate-independent model is recovered by assigning = 0.0 s. A certain amount of elements, which coincide with the notch geometry, are preeroded (red elements in Figure 7B. In this way, a numerical crack, that is, eigenerosion, is introduced into the setup, and consequently, a fracture problem with a preexisting crack is solved and no special assumption for the crack initiation criterion is made. Elements located at singularity regions, which develop at load applications and boundary conditions, are prone to initial evolution of erosion. The stiffness loss of these elements would compromise the numerical solution of the boundary value problem. To avoid this, the elements situated at these regions are excluded from the eigenerosion scheme. Figure 8

F I G U R E 8 Simulation and experimental
results for the vertical displacement of the midpoint A as function of the reaction forces (see Figure 7). As it can be seen in Figure 8, the simulation results show a very good agreement with the experimental ones. The "zig-zag"-like shape of the numerical results graph occurs due to the binary approach of the eigenerosion, where an element can either be intact or eroded avoiding any other intermediate configuration.
Obviously, the choice of crit hom can change the results dramatically by making the erosion less or highly depending on plasticity. To identify this parameter, the evolution of plasticity is first determined when tensile strains of around 0.15‰ are achieved at the crack tip. This is actually the maximum value of tensile strains at the fracture zone proposed in fib Model Code for Concrete Structures 2010. At this state, the evolution of the homogenized plasticity parameter around the crack tip ranges between 2.0 × 10 −5 < hom < 5.0 × 10 −5 1/MPa. Thus, crit hom is chosen in the same range. Based on the original variational formulation of fracture and on Griffith criterion, new cracks form when stored elastic energy is sufficiently high to overcome the energy necessary for the advance of the crack. If an elastic perfectly plastic medium with low-yield limit is considered, plastic deformations cumulate at the crack tip reducing the increment of the elastic energy and precluding the attainment of the energy criterion. This concept is tackled with the present formulation where the evolution of plastic deformations at the crack tip, at the same time, is a precursor of crack propagation and prevents localization.
The eigenerosion representation and the evolution of plasticity for this example are shown in Figure 9. The crack starts as purely mode I failure until it reaches the upper part of the beam. At this point, due to bending, compressive stresses increase rapidly through permanent activation of the compression cap (see Equations 20 and 21). Hence, plasticity due to compression dominates the uncracked elements ahead of the crack. Therefore, as crack propagation under compression requires more energy (only the deviatoric part of the free-energy function is the crack-driving force in this case), remaining reactions are shown in Figure 8. Clearly, one can safely conclude that the minimizers of this problem are global ones due to the mesh's homogeneity and complexity of the crack path.

Compact tension specimen
The last example is a standard compact tension specimen tested at different loading rates. Figure 10A shows the geometry and boundary conditions of the specimen (dimensions in millimeters). The sample has a thickness of 25 mm. To verify the mesh independence of the proposed approach, three different mesh sizes are introduced for the finite element (FE) discretization with the respective minimum element lengths h min,1 = 3.05 mm, h min,2 = 1.52 mm, and h min,3 = 0.76 mm. The second FE discretization is depicted in Figure 10B. The specimen has a depth of 25 mm. This test is documented experimentally in Reference 50, where also some numerical analyses in terms of rate sensitive dynamic fracture occur.
For the simulation at hand, the load application frames are not modeled but rather simplified by applying the load and boundaries directly at the crack faces. A displacement-controlled simulation is conducted by fixing the upper part of the notch and applying different loading rates at its bottom part. Four different displacement rates are depicted here, namely, 35, 1300, 3300, and 4500 mm/s. The importance of this test is two-fold. First, it is possible to observe the rate effects from the material point of view, and second, in terms of fracture mechanics, one can observe branching for high-loading rates. The used material parameters are presented in Table 3. The critical ERR is chosen as G c = 5.5 × 10 −4 kN/mm and the parameter controlling the effect of plasticity on the evolution of eigenerosion is chosen to hom crit = 1.0 × 10 −6 1/MPa. The material parameters are identified as follows. First, depending on the type of concrete, we chose an initial set of parameters. Then, we fit the almost quasi-static case (35 mm/s) by changing the parameters for the numerical simulation. These parameters are then kept constant and only the loading rate is changed.
The implicit Newmark time integration scheme is used to compute the solution of the dynamic process. To avoid spurious oscillations, the Newmark parameters are chosen as = 0.5 and = 0.7. As the iterative nature of the implicit scheme is unconditionally stable, theoretically, there is no restriction in the choice of the time step Δt. However, to capture wave propagation in the discretized medium, Δt is often calculated depending on the wave speed in the elastic solid and  the minimum element length to correctly pass the wave through the mesh. Caution should be paid as well not to take very small values as it increases the numerical instability due to the division by Δt 2 during the solution. For the different mesh sizes considered here and referring to the speed of the shear wave in concrete, the time step Δt is calculated accordingly. As the deviatoric part of the energy yields very high values at the load application location and very close to it, the elements at this region would instantaneously be eroded once the load is applied. To circumvent this situation, the evolution of erosion is artificially avoided at this region to have a realistic modeling of the experimental setup. The crack patterns for different rates of the experimental results are shown in Figure 11, meanwhile, the simulation results are depicted in Figure 12. For convenience, only three different rates are visualized in these figures. Comparing Figures 11 and 12, one can realize the very good agreement between simulation and experiments with respect to the crack patterns. As expected with an eigenerosion-based method, results are independent of the mesh size (see Figure 12). For the loading rate of 35 mm/s, a horizontal crack is achieved. While the rate increases, the crack inclines not only toward the load application region but also toward branches. For high-loading rates, the amount of accumulated energy becomes so large that it requires more than one element to be dissipated by the fracture mechanism. In this way, branching occurs, and henceforth, the energy is dissipated by several cracks. Spurious erosion of elements perpendicular to the crack direction shows that the minimizers in this example cannot be proven to be global. Nevertheless, the local minimizers can be accepted as a satisfactory condition for the solution of the crack-tracking problem instead of a necessary condition. Despite the well satisfying results depicted in Figures 11 and 12, the current approximation is far from being perfect. To have a realistic physical model of wave reflection and transmission that occur between the crack faces, appropriate contact modeling is required. Moreover, as a heterogeneous material, concrete constituents effect the global structural response, and a more accurate modeling at smaller scales is often necessary (see eg, application to metaconcrete in Reference 51). The measured reaction forces as a function of the vertical displacements of the notch at the lower left corner for the element size h min,2 are shown in Figure 13. One can observe in Figure 13, the response at different rates and the peaks (maximum reactions) that are achieved. The proposed formulation is able to model a raise in strength for increasing loading rates. As it is often documented in literature, with the increase of the loading rate, an increase in concrete strength is observed (DIF). For this simulation, at static loading condition, we achieved a maximum reaction of R stat = 2.15 kN. Considering the maximum load for the other rates, the well-known DIF-curve for reaction forces is shown in Figure 14, where the results documented in literature and those from the simulations are illustrated. As documented in Reference 52, the literature results are taken from various authors that tested concrete structures at different loading rates using setups such as explosion tests, spallation tests, direct tensile tests, Split Hopkinson tests, and so on. It is clear from the graph that at high-loading rates, the strength increases rapidly. The reason for this kind of behavior is still under investigation in literature, but it is accepted widely that the structural inertia, which delays the fracture process, plays a major role. Another reason for this effect is also the change of stress distribution in dynamic loads that is introduced by the wave propagation. Water content also plays a role in the increase of DIF. At the current work, the increase in strength is directed from two mechanisms, namely, the material effect that delays the damage initiation due to viscous effects in F I G U R E 13 Reaction forces at the boundaries as a function of vertical displacements at the load application F I G U R E 14 DIF for reactions. DIF, dynamic increase factor the evolution of plasticity and the structural effect due to inertia of the system. Nevertheless, the reader should be very critical when examining this graph, as the measurements, evaluations, and type of tests used by different authors vary significantly. Moreover, caution should be paid to the effect of cracking mechanism while increasing the strength . From the mechanics point of view, an increase in the crack length implies larger energy dissipation, which consequently is reflected in decreased reactions. The same conclusion can be drawn when branching occurs. Thus, building a general representation of the graph shown in Figure 14 should always be conducted in an objective manner.

CONCLUSION AND DISCUSSION
In this work, two main developments have been introduced. First, we proposed an approach to describe cracks for nonlinear material behavior within the variational framework of the eigenfracture. Second, in the framework of the microplane approach, a rate-dependent consistency viscoplastic model was formulated based on the Drucker-Prager yield function equipped with a tension and compression cap. Both developments are combined in a solution strategy to model failure of quasi-brittle material such as concrete. Through three representative numerical examples, we have shown the potential of the proposed approach and established a solid ground for further developments in nonlinear variational fracture.
The original energetic description has been modified by scaling the surface fracture energy to account for plasticity evolution in the material. Demonstrating that the resulting energy term is lower semicontinuous and the bulk energy is quadratic positive definite, we built a regularized variational problem whose convergence is ensured in the spirit of Γ-limit analysis. Nevertheless, the used eigenerosion scheme does not guarantee the solution for global minimizers. For the problems at hand, local minimizers are accepted as satisfactory condition of the variational formulation. Moreover, mesh dependency of the surface energy can also produce local minimizers leading to preferential directions of the crack path. While remeshing can be beneficial in finding the global minimizers, mapping the state variables is not always a trivial procedure.
Another noteworthy remark is that in the absence of singularities (eg, notches), it has been proven in variational brittle fracture that a crack initiates abruptly and can propagate either frantically or monotonically depending on the convexity of the actual surface energy function. Consequently, abrupt failure sequences (or even instantaneous failure) occur (see the theoretical proofs in Reference 3 and also numerical results in Reference 2). This is observed in numerous simulation results obtained by different methods, which stem from the variational formulation of brittle fracture. The coupling to the inelastic internal variable, introduced in this work, is a numerical tool that bypasses the brittleness of the classical eigenfracture scheme, as the evolution of plasticity and the rate effects adequately deliver a more ductile solution.
To model the constitutive behavior of concrete at dynamic loads, we introduced a rate-dependent viscoplastic consistency model by constructing a smooth three-surface Drucker-Prager yield function with caps as a dependency also on the rate of the internal variable. Thus, isotropic hardening can occur not only due to the hardening variable but also due to its rate. The generality of the formulation ensures switching between the rate-dependent and rate-independent description without any extra effort. From the material point of view, considering viscous effects in concrete modeling are motivated by the presence of water in small quantities in its microstructure. 50 The increase in strength (DIF) observed in the dynamic example is attributed to two different mechanisms from the numerics point of view. On the one hand, erosion is delayed due to viscous effects. On the other hand, the inertia of the crack faces plays also a major role in delaying the failure. The introduced rate effect is independently applied on each microplane. Thus, an induced anisotropy is also theoretically feasible for the rate effect variable. In previous works done on the microplane approach, the rate effect is dependent on the macroscopic strain tensor, making it an invariant quantity on different microplanes. Based on our experience with the present implementation, no arguments that conflict with the physics behind this formulation have been encountered. Nevertheless, the focus of the present research is on the differences of the global material response coming from the two approaches that account for rate dependency.
The straightforward description of the eigenerosion algorithm, its easy implementation strategy, and the low computational effort make the proposed approach very attractive for failure modeling. Moreover, the introduced constitutive relation that models concrete behavior is able to accurately capture the characteristics of the material realistically. Certainly, the crack-tracking task at hand faces numerous difficulties due to the number of material parameters used in the constitutive formulation and the challenges regarding the appropriate crack-driving force for nonlinear fracture problems. Combining these obstacles in a dynamic solution scheme requires further intensive research.