Efficient and robust numerical treatment of a gradient-enhanced damage model at large deformations
EEfficient and robust numerical treatment of agradient-enhanced damage model at large deformations
Philipp Junker, Johannes Riesselmann, Daniel Balzani
Chair of Continuum Mechanics, Ruhr University Bochum, Bochum, GermanyCorresponding author:Philipp Junker, (cid:0) [email protected]
Abstract
The modeling of damage processes in materials constitutes an ill-posed mathematical problemwhich manifests in mesh-dependent finite element results. The loss of ellipticity of the discretesystem of equations is counteracted by regularization schemes of which the gradient enhancementof the strain energy density is often used. In this contribution, we present an extension of theefficient numerical treatment, which has been proposed in [1], to materials that are subjectedto large deformations. Along with the model derivation, we present a technique for elementerosion in the case of severely damaged materials. Efficiency and robustness of our approach isdemonstrated by two numerical examples.
The modeling of damage at finite strains roots back to the eighties of the 20 th century when firstphenomenological models related to material degradation were developed. These models weremainly constructed based on the original ideas for small strains where damage was considered asa reduction of the effective cross section area on which the load acts resulting from the nucleationand evolution of cavities. The reduced cross section thus softens the material which gives riseto the stress- and strain-softening effects. On a material point level, this means that the stressreduces with increasing strains after a critical threshold value is reached. This phenomenologicalapproach of describing the complex microstructural processes during material failure, that mightbe accompanied by the evolution of dislocations, micro-cracks and delamination effects, providesa convincing method for material modeling. Unfortunately, softening may also be accompaniedby a loss of ellipticity of the system of governing equations once a certain amount of damage hasevolved which renders boundary value problems ill-posed. An obvious indicator for such problemsare severe dependencies of the finite element results on the individual mesh discretization: thesmaller the finite element size is chosen, the more localized resulting damage bands will be.Furthermore, and even worse, the numerical solution of such problems does not converge andmight be unstable regarding small perturbations of the boundary conditions, i.e., they lacknumerical robustness. These effects are not only observed at the level of spatial distribution of thedamage variable but also for the global behavior in terms of reaction forces [2, 3]. Consequently, a1 a r X i v : . [ c s . C E ] F e b athematical correction of the problem is inevitable. This process is referred to as regularizationwhich might be performed in time and space. Overview articles may be found, e.g., in [4] and [5].Alternatively, a relaxed variational formulation might be employed which relies on a carefulanalysis of the convexity properties of the underlying energy. This energetic analysis demands adirect formulation of the damage variable as pure function of the deformation state such that thecondensed energy (energy plus time-integrated dissipation) may be relaxed. It is thus convenientto formulate the model in a time-incremental way, cf. [6, 7, 8, 9]. Thereby, a dissipative materialbehavior may be described by an incrementally hyperelastic representation, which is thus open foranalysis with respect to generalized convexity conditions. This enables the construction of convexhulls (see [10] for small strains and [11] and [12] for finite strains) which on the one hand ensuresthe recovery of ellipticity and on the other hand guarantees the existence of minimizers. Incontrast to spacial regularization approaches, no internal length scale parameter has to be chosensince the micromechanical parameters are obtained as result of the convexification procedure.On the other hand, the effort at the material point level may be increased if a nonconvexminimization problem is solved as part of the convexification at each integration point, whichyields an increased expense in the assembling procedure. Another drawback may be seen inthe fact that although the stress softening hysteresis as explained above may in principle becaptured, a softening in the sense of negative tangents in stress-strain curves can not yet bemodeled at the material point. One exception so far is the approach of Schwarz et al. in [13] forsmall strains, where a simplified microstructure consisting of damaged, linear elastic materialsand whose resulting convexified energy indeed enables the description of strain-softening.In the context of time-wise regularization, viscous effects are included into the model. Ex-amples for such models may be found in [14, 15, 16, 4, 17, 18]. These models introduce arate-dependence of the damage evolution which damps the local reduction of stiffness. Conse-quently, damage evolves also at neighbored material points to keep the amount of dissipatedenergy constant. A huge benefit of such models (as well as of relaxed incremental formulations)is their implementability in any finite element framework using solely the material law inter-face. However, viscous damage models possess a remarkable drawback: the viscosity associatedto damage evolution has to be chosen quite high, reaching values that are unphysical, in orderto obtain the desired regularization. Furthermore, the effectiveness of viscous regularizationstrongly depends on the applied loading velocity which does not necessarily need to agree withexperimental velocities. Finally, the interplay of viscosity and loading velocity has to be set andthus well-posedness has to be proven for each individual boundary value problem rendering suchapproaches rather impractical for application.Spacial regularization “homogenizes” the damage evolution on a local level. To this end,either integral formulations are introduced or a gradient enhancement of the strain energy isemployed. A major benefit of such models as compared to viscous regularization is that itcan be shown that they are inherently well-posed. Examples for spatially regularized damagemodels may be found, e.g., in [2, 3, 19, 20, 21, 22, 23, 24]. Also for the finite-strain regime, theyhave been used, cf. [25, 26, 27, 28, 29, 30, 31, 32]. Such models may be applied to arbitraryloading velocities since any dependence of the material behavior on this boundary condition isincluded in a physically sound manner. Furthermore, they might be coupled to other physicalprocesses as for instance plasticity and temperature dependence without unphysical parametersbeing present in the complete system of model equations. For gradient regularization, we referto the prototype approach of [20, 21] and similarly [22], and its finite-strain extension in [31],where a nonlocal function carries information on the damage state. This procedure transformsthe associated equation for the damage state from a local evolution equation at the materialpoint (as is the case for viscous regularization) to a field equation (integral equation for integralregularization or partial differential equation for gradient regularization) which has to be solved2long with the balance of linear momentum for the displacements. Consequently, the numberof unknowns increases during numerical treatment, e.g., the number of nodal variables increasesfrom three (displacements) to four (displacements + damage state variable) in a mixed finiteelement implementation. While thermodynamic consistency and mesh-independent simulationscan be ensured with these approaches, unfortunately, the increased number of degrees of freedomin the discrete global system of equations significantly increases computational costs.To limit the usual numerical deficiency of gradient regularization, we present in this contri-bution a novel finite strain gradient damage formulation. A similar approach for the small strainregime has been introduced in [1] (see also [33] for a detailed numerical analysis including meshadaptivity). We begin with the variational derivation of the damage model that is formulatedin terms of a partial differential inequality which complements the partial differential equationdescribing the balance of linear momentum. Then, we apply the neighbored element methodwhich has proven to be beneficial for a robust and time-efficient numerical evaluation of the setof governing equations given by the balance of linear momentum and an additional field equa-tion. Examples have been given in the context of small-strain damage [1, 33] and small-strain[34] and finite-strain topology optimization [35]. The neighbored element method combines afinite element approach for the balance of linear momentum and a finite difference approach forthe additional field equation along with operator splitting techniques. Thereby, the neighboredelement method conserves the size of the global system of equations resulting from the finiteelement formulation for the displacements: to be more precise, the number of nodal unknownsfor the finite elements remains unchanged from a corresponding purely elastic formulation. Thepartial differential inequality for the damage state, however, is solved through a correspondinglymodified internal update procedure. Furthermore, a novel stabilization technique based on el-ement erosion methods enables us to obtain numerical results even for substantial deformationstates including regions of highly damaged material for which the model is close to the lossof its validity. The numerical results show mesh-independence while requiring minimal extracomputation time compared to purely hyperelastic simulations. There exist various strategies for the derivation of fundamental material models of which anextended Hamilton principle offers a variational approach. Its stationarity conditions agree withthe 2 nd Law of Thermodynamics and Onsager’s principle by construction if physically reasonableansatzes for the strain energy density and the dissipation function are chosen. For details on theextended Hamilton principle and its relation to thermodynamics and other modeling strategies,such as the principle of virtual work or the principle of the minimum of the dissipation potential,we refer to [36].The extended Hamilton principle is related to the stationarity of the action function. It canbe shown, cf. [36], that it constitutes as the following condition: for the quasi-static case, thesum of the Gibbs energy G , which is also referred to as total potential, and the work due todissipative processes D tends to be stationary:(1) H [ u , α ] := G [ u , α ] + D [ α ] → stat u ,α . A detailed investigation of the extended Hamilton principle as unifying theory for coupled prob-lems, i.e. thermo-mechanical processes, and dissipative microstructure evolution has been pre-sented in [36].In (1), the displacements are denoted by u and the state of microstructure is expressed in3erms of the internal variable α . The Gibbs energy reads(2) G [ u , α ] = (cid:90) Ω Ψ( C , α ) d V − (cid:90) Ω b (cid:63) · u d V − (cid:90) ∂ Ω t (cid:63) · u d A with the strain energy density Ψ, the prescribed body forces b (cid:63) , and the traction vector t (cid:63) . Theintegrals are evaluated for the body’s volume Ω and its surface ∂ Ω in its reference configurationwith the position vector X . The deformation is measured by the right Cauchy Green tensor C := F T F with the deformation gradient F = ∂ x /∂ X = I + u ⊗ ∇ and the spatial coordinate x = X + u in the current configuration. Throughout this contribution, the nabla operator iscomputed with respect to the reference configuration, i.e., ∇ ≡ ∇ X .The work due to dissipative processes is given by(3) D [ α ] := (cid:90) Ω p diss α d V when the non-conservative force p diss performs work “along” the microstructure state which isdescribed in terms of the internal variable α . The stationarity condition of (1) thus reads(4) δ H [ u , α ]( δ u , δα ) = δ G [ u , α ]( δ u , δα ) + δ D [ α ]( δα ) = 0 ∀ δ u , δα with(5) δ G [ u , α ]( δ u , δα ) = (cid:90) Ω (cid:16) δ u Ψ( C , α ) d V + δ α Ψ( C , α ) (cid:17) d V − (cid:90) Ω b (cid:63) · δ u d V − (cid:90) ∂ Ω t (cid:63) · δ u d A .
The non-conservative force is assumed to be a material-specific quantity which is determinedonce the process conditions are set, similarly as the external forces b (cid:63) and t (cid:63) . Consequently, itdoes not follow from a stationarity condition but it needs to be modeled. This implies δ α p diss = 0such that the variation of the dissipated energy reads(6) δ D [ α ]( δα ) = (cid:90) Ω p diss δα d V .
Since the variations for displacements and microstructural state are independent, we obtain(7) (cid:90) Ω δ u Ψ( C , α ) d V − (cid:90) Ω b (cid:63) · δ u d V − (cid:90) ∂ Ω t (cid:63) · δ u d A = 0 ∀ δ u (cid:90) Ω δ α Ψ( C , α ) d V + (cid:90) Ω p diss δα d V = 0 ∀ δα . Specifications of the strain energy density Ψ and the non-conservative force p diss allow the appli-cation of the latter equations for the derivation of models for various materials. Based thereon,let us provide specific formulas for a gradient-enhanced damage model by first setting the internalvariable α as damage variable. Then, we define the strain energy density by(8) Ψ( C , α ) = (cid:0) − D ( α ) (cid:1) Ψ ( C ) + 12 β ||∇ f || with the damage function D ( α ) := 1 − f ( α ), the damage variable α and f ( α ) = exp( − α ). Theparameter β serves as regularization parameter and thus controls the thickness of the damagezone. The formulation of the strain energy density is the same approach as in [1]; however, incontrast to the original work which was restricted to small deformations, we now make use of ahyperelastic strain energy density of the undamaged material Ψ = Ψ ( C ).4he non-conservative force is usually derived from a dissipation function Π diss such that(9) p diss = ∂ Π diss ∂ ˙ α holds. For rate-independent microstructural evolution, which is the case for the present modelfor damage evolution, functions have to be used that are homogeneous of order one. Hence, thedissipation function(10) Π diss := r | ˙ α | is used [1] from which(11) p diss = ∂ Π diss := r ˙ α | ˙ α | for ˙ α (cid:54) = 0 {− r, r } for ˙ α = 0follows. The parameter r is referred to as dissipation parameter and it will be become obviousthat it represents an energetic threshold value for the onset and evolution of damage. It is worthmentioning that the partial derivative in (9) transforms to a subdifferential in our case since ∂ Π diss /∂ ˙ α is not unique at ˙ α = 0 for our rate-independent ansatz for Π diss in (10). Then, thestationarity condition in (7) result in the weak form of the balance of linear momentum(12) (cid:90) Ω S : 12 δ C d V − (cid:90) Ω b (cid:63) · δ u d V − (cid:90) ∂ Ω t (cid:63) · δ u d A = 0 ∀ δ u , where S = 2 ∂ Ψ( C , α ) /∂ C denotes the 2 nd Piola Kirchhoff stress tensor. The stationarity con-dition in (7) reads(13) (cid:90) Ω f (cid:48) Ψ δα d V + (cid:90) Ω β ∇ f · ∇ ( f (cid:48) δα ) d V + (cid:90) Ω ∂ Π diss δα d V = 0where we considered δ α f = f (cid:48) δα . Integration by parts of the second term results in(14) (cid:90) Ω βf (cid:48) ∇ f · ∇ δα d V = (cid:90) ∂ Ω βf (cid:48) n · ∇ f δα d A − (cid:90) Ω βf (cid:48) (cid:52) f δα d V and thus, (13) transforms to(15) − (cid:90) Ω (cid:0) f Ψ − βf (cid:52) f − ∂ Π diss (cid:1) δα d V − (cid:90) ∂ Ω βf n · ∇ f δα d A = 0due to f (cid:48) = − f . The local evaluation of the first integral in (15) results in the differentialinclusion(16) f Ψ − βf (cid:52) f − ∂ Π diss (cid:51) ∀ X ∈ Ωdue to the set-valued character of the subdifferential ∂ Π diss . Therefore, it is convenient to performa Legendre transformation such that we transform the dissipation function, which is formulatedin terms of the thermodynamic flux ˙ α , into a function ˜Π diss which depends on the thermodynamicforce p := f Ψ − βf (cid:52) f . Consequently, we obtain(17) ˜Π diss = sup ˙ α (cid:110) p ˙ α − Π diss (cid:111) = sup ˙ α (cid:110) | ˙ α | ( p sgn ˙ α − r ) (cid:111) . α = { , } and the Legendre transform reads(18) ˜Π diss = (cid:40) p − r ≤ ∞ elseFrom the first case, we identify(19) (cid:40) p < r : ˙ α = 0 p = r : ˙ α > . It is thus convenient to introduce the indicator function Φ := p − r which separates stationarity ofthe damage state for Φ < r hasa similar meaning as a scalar-valued yield stress. However, it is an energetic threshold value forthe present case of damage. Finally, we can collect the governing equations for damage evolutionin the following form(20) ˙ α ≥ , Φ := f Ψ − βf (cid:52) f − r ≤ , Φ ˙ α = 0 ∀ X ∈ Ωwhich are identified as Karush Kuhn Tucker conditions. The surface integral constitutes asNeumann condition(21) n · ∇ f = 0 ∀ X ∈ ∂ Ωfor f . Here, n denotes the normal vector on the surface ∂ Ω in the reference configuration. Moredetails on the fundamentals of the damage model can be found in the original publication [1].
The model above consists of two unknowns: the displacement field u and the damage variable α . They can be determined by solving(22) (cid:90) Ω S : 12 δ C d V − (cid:90) Ω b (cid:63) · δ u d V − (cid:90) ∂ Ω t (cid:63) · δ u d A = 0 ∀ δ u f Ψ − βf (cid:52) f − r ≤ ∀ x ∈ Ωwith the Dirichlet and Neumann boundary conditions u = u (cid:63) ∀ X ∈ ∂ Ω u and F Sn = t (cid:63) ∀ X ∈ ∂ Ω σ with ∂ Ω = ∂ Ω u ∪ ∂ Ω σ and ∂ Ω ∩ ∂ Ω σ = ∅ for (22) . For (22) , the initial condition f ( t = 0) = 1 ∀ X ∈ Ω and the Neumann condition n · ∇ f = 0 ∀ X ∈ ∂ Ω apply.From (22) , we recognize that a partial differential inequality has to be solved for the primalvariable f . Once f is determined, the value of the damage variable could be computed subse-quently for each material point X by α = − log[ f ( X )]; however, this information is of minorinterest. More importantly, the damage function D = 1 − f can be computed which governs thestresses by(23) S = (1 − D ) S with S := ∂ Ψ ∂ C . The partial differential inequality in (22) is rather untypical in the context here and a standardfinite element treatment would be complex due to the subdifferential in its weak form, cf. (13).6urthermore, additional degrees of freedom at the nodes would be present in the numericalsolution scheme. To avoid both drawbacks, we follow [1] and make use of the neighbored elementmethod: standard finite element approaches are applied to the weak form of the balance of linearmomentum whereas a finite difference approach for unstructured grids is used for discretizingthe strong form of the evolution equation for f . The neighbored element method is completed byutilization of an operator split, i.e., the discretized equations are solved in a staggered manner.This procedure of combining staggered FEM and FDM has also been proven advantageous forthe small-strain regime in [1]. In [33], it has been demonstrated that convergence is achievedwhich justifies the operator split. The finite element treatment of the nonlinear weak form ofthe balance of linear momentum is standard such that we dispense with a detailed presentation.More details on the finite element method can be found in standard textbooks as in, e.g. [37].The only modification is, of course, that both the stress and the stiffness are scaled by thedamage state, cf. (23).For the derivation of an appropriate finite difference method that also operates on unstruc-tured grids, we employ a Taylor series expansion up to order two in order to approximate f .Usually, the value of the function, i.e., its zeroth derivatives, along with the derivatives of higherorder are known when a Taylor series expansion is used. Then, due to the chosen spatial in-crement, the value of the function can be approximated at a neighbored spatial point. Thisoperation is inverted here: the value of the damage function at specific points are all known.These points are the centers of gravity of all finite elements. Then, the spatial increments inthe Taylor series expansion are given by the spatial increments of the centers of gravity; theyremain fixed if no mesh adaptation is employed during the computation. Consequently, the onlyunknown in the Taylor series expansion are the (mixed) derivatives of order ≥
1. Collecting anappropriate set of neighbored elements, a linear system of algebraic equations can be constructedto compute the individual (partial) derivatives from which the local value of the Laplace operatorfollows by simple summation of the unmixed derivatives of order two.Our finite differences approach at unstructured meshes can be recast in mathematical formu-las by first introducing the Taylor series expansion as(24) f ( N ( e ) k ) = f ( e ) + (cid:88) o A ( i,k ) o f ( e ) ∂,o where f ( e ) ∂,o stores the value of the partial derivatives of varying order o ≥
1, such that a matrixincluding all derivatives can be defined as(25) f ( e ) ∂ := (cid:18) ∂f ( e ) ∂X ∂f ( e ) ∂Y ∂f ( e ) ∂Z ∂ f ( e ) ∂X∂Y ∂ f ( e ) ∂Y ∂Z ∂ f ( e ) ∂X∂Z ∂ f ( e ) ∂X ∂ f ( e ) ∂Y ∂ f ( e ) ∂Z (cid:19) . The quantity A ( i,k ) o comprises all spatial increments with varying power such that again theassociated column matrix is defined as A ( k ) := (cid:32) X ( N ( e ) k )∆ Y ( N ( e ) k )∆ Z ( N ( e ) k )∆ X ( N ( e ) k )∆ Y ( N ( e ) k )∆ Y ( N ( e ) k )∆ Z ( N ( e ) k )∆ X ( N ( e ) k )∆ Z ( N ( e ) k )∆ (26) 12 (cid:16) X ( N ( e ) k )∆ (cid:17) (cid:16) Y ( N ( e ) k )∆ (cid:17) (cid:16) Z ( N ( e ) k )∆ (cid:17) (cid:33) with X ( N ( e ) k )∆ := X ( N ( e ) k ) − X ( e ) , Y ( N ( e ) k )∆ := Y ( N ( e ) k ) − Y ( e ) , Z ( N ( e ) k )∆ := Z ( N ( e ) k ) − Z ( e ) .
7o ensure that only the required number of elements and thus extrapolation points is used, thequantity N ( e ) k has been introduced: N ( e ) is the set of neighbored elements around element ( e )and, hence, N ( e ) k returns the global element number for k th neighbor. All column matrices A ( k ) in the neighborhood N ( e ) k of the element of interest ( e ) are collected in the matrix A ( e ) . Then,by defining(27) f ( e )∆ ,k := f ( N ( e ) k ) − f ( e ) , we can shortly write for the Taylor series expansion (24)(28) f ( e )∆ = A ( e ) f ( e ) ∂ such that the unknown vector of derivatives follows from(29) f ( e ) ∂ = A ( e ) − f ( e )∆ . The Laplace operator can be computed by aid of f ( e ) ∂ according to(30) (cid:52) f ( e ) = ∂ f ( e ) ∂X + ∂ f ( e ) ∂Y + ∂ f ( e ) ∂Z = f ( e ) ∂, + f ( e ) ∂, + f ( e ) ∂, , cf. (25). Consequently, only some components of f ( e ) ∂ are of interest. Introducing(31) a := (cid:0) (cid:1) , the Laplace operator is simply given by(32) (cid:52) f ( e ) = l ( e ) · f ( e )∆ where the vector(33) l ( e ) := a T A ( e ) − can be computed once for each element ( e ) in advance of the actual solution of the boundaryvalue problem. In contrast, the evaluation of (32) has to be repeated during computationalevaluation of (22) to find the (current) field f ; this evaluation, however, is computationallyvery cheap. It is worth mentioning that the dimension of A ( e ) depends on cardinality of the setof neighbored elements N ( e ) : to close the system of equations in (28), i.e., to ensure that A ( e ) is aregular matrix, the cardinality of the set of neighbored elements has at least to equal the lengthof f ( e ) ∂ = 9 in the considered three-dimensional case. For boundary elements, however, we usuallyfind a smaller cardinality of N ( e ) , i.e. less than nine neighbored elements can be identified. Tocircumvent this problem, we introduce ghost elements at the boundary which mirror the valueof the damage field inside of Ω to the outer vicinity of the boundary ∂ Ω which equals the usualtreatment in the finite difference method. It is worth mentioning that (21) is fulfilled identicallyby usage of ghost elements. For the elements inside of Ω, the opposite case is present: consideringall neighbored elements around ( e ), i.e., the stencil of six elements plus eight diagonal elementsplus twelve elements in each “plane”, results in an overdetermined system of equations for (28).Here, two possibilities can be employed: the first option is to take into account element ( e ) andall 26 neighbored elements and compute the inverse of A ( e ) by use of the right-inverse as(34) A ( e ) − = A ( e ) T (cid:16) A ( e ) A ( e ) T (cid:17) − , e ) as neighborhood. For sufficiently regular meshes, A ( e ) contains full rank and is regular, cf. [33]. In all computational results we present later, we madeuse of the latter option. We refer to [1, 34, 33] for more details on the neighbored elementmethod.Finally, it is mandatory to chose a strategy for solving the discretized form of (22) given by(35) Φ ( e ) = f ( e ) ¯Ψ ( e )0 − βf ( e ) l ( e ) · f ( e )∆ − r ≤ ( e )0 := 1Ω ( e ) (cid:90) Ω ( e ) Ψ d V .
The Laplace operator turns (35) into a system of coupled inequalities which could be solved ina monolithic manner. Although being the obvious strategy, the coupling demands the usage of(external) solver routines. In combination with the inequalities, this renders such a monolithicstrategy inappropriate to be implemented into an arbitrary finite element program. Therefore,a simple yet successful numerical solution procedure has been proposed in [1]: the system ofcoupled inequalities is solved by means of a Jacobi method. To be more precise, the couplingis ignored at first glance and each inequality is checked individually such that an update of f isperformed for all elements with Φ ( e ) >
0. Of course, a modification of f ( e ) may have an impact onthe indicator function in a different element ( k ) due to the Laplace operator which “transports”the information. Therefore, the update has to be repeated until Φ ( e ) ≤ ∀ ( e ). Usually, such anupdate scheme is known to be rather time consuming since the numbers of repetitions needed forconvergence depends on the number of unknowns which equals the number of finite elements inour case. This should turn this strategy to be more expensive than monolithic solution schemes.However, the damage evolution is limited to only a small number of elements as compared totheir total number and thus, an update is mandatory only for a negligible fraction of inequalities.As will be shown by our numerical experiments, the staggered solution scheme is still very fastand simulations involving damage evolution consume hardly more computation time. A similarbehavior has already been demonstrated for the small deformation setting of damage modelingin [1] and even for an evolutionary method for topology optimization both at small strains in[34] and at finite strains in [35]. After the damage function has been updated, the next time stepis investigated without checking the impact of the updated damage field on the displacementsfor the current time step. This might be interpreted as (additional) operator step and may seemquestionable whether a reliable result is obtained by this method. However, time steps withvanishing increment are mandatory for an appropriate resolution e.g. of the force/displacementdiagram. For this case, the method of operator splits converges. It has been numerically shownin [1] that the artificially included pseudo-viscous behavior due to the operator step does notinfluence the final result on a relevant scale. Furthermore, it is worth mentioning that a detailednumerical study on the numerical behavior of the damage model for the small deformation settingin [33] revealed that almost identical results are obtained when both fields (displacements anddamage) are repeatedly updated until a common convergence is achieved. The update of eachinequality is performed by a one-step Newton procedure as(37) f ( e ) ← f ( e ) − Φ ( e ) dΦ ( e ) with(38) dΦ ( e ) = Ψ ( e )0 − β l ( e ) · f ( e )∆ + βf ( e ) l ( e ) · equals the cardinality of the set N ( e ) (9 for our chosen method fordefining the neighborhood) and all elements in are set to one. The one-step Newton method hasbeen proven beneficial due to its smooth convergence behavior for neighbored elements causedby the undershooting. Additionally, this ensures a monotonic decrease of the damage functionand thus agrees to (20) : the indicator function Φ ( e ) is a convex quadratic function in f ( e ) , andsince the update for f ( e ) is only performed for Φ ( e ) >
0, it holds for the derivative dΦ ( e ) > ← Φ ( e ) , and thus f ( e ) ismonotonically decreasing which implies ˙ α ≥ Algorithm 1
Element erosion strategy for , . . . , n e do (cid:46) apply to all finite elements if D ( e ) = D crit then (cid:46) check damage criterionset r ( e ) = (cid:46) eliminate element residualset d r ( e ) d ˆ u ( e ) = (cid:32) d r ( e ) i dˆ u ( e ) j (cid:33) = s crit ( δ ij ) (cid:46) erode element ( e ) end ifend for The above mentioned numerical procedure is analogous to the one used for the small deforma-tion setting in [1]. However, a remarkable difference is accompanied by the large deformation set-ting: here, severe deformations that occur during damage evolution are correctly described. This,of course, is the primal intention of using the large deformation setting. Unfortunately, this sen-sitivity might result in numerical instabilities causing det ¯ F ( e ) < F ( e ) := (cid:82) Ω ( e ) F d V / Ω ( e ) is the averaged deformation gradient in element ( e ). To be more precise, deformation states arecomputed that do not correspond to a physically accurate material behavior: once the deforma-tion state has reached some critical threshold, cracks will be present in a real material such thata consideration of the specimen as one continuous body cannot be justified and a description ina continuum mechanical sense is not justified anymore. To circumvent this numerical artifactand also to improve physical accuracy, a stabilization technique in terms of element erosion isapplied. To this end, we define a critical value D crit when damage modeling is not suitable any-more but rather cracking needs to be considered. To correctly account for this physical behaviorin our finite element simulations, we thus modify the element stiffness by setting it to a diagonalmatrix with a virtual remaining stiffness value of s crit representing an eroded element throughwhich rank deficiency of the affected elements is avoided. Finally, we set the residual forces ofthis element to zero and no further evolution of the damage value is considered anymore at thispoint.Summarizing, we follow the element erosion technique in Alg. 1 and employ Alg. 2 for thedamage update. The flowchart in Fig. 1 visualizes the neighbored element approach includingan interplay between a finite element and finite difference update for the displacements and f ,respectively. In this section, the introduced approach is numerically tested in a 3D finite strain setting. To thisend, corresponding 8-node hexahedral finite elements were implemented and evaluated within10 lgorithm 2
Damage updateinput ¯Ψ ( e )0 ,n +1 = (cid:82) Ω ( e ) Ψ ,n +1 d V / Ω ( e ) (cid:46) averaged energy for each finite element ( e )¯ F ( e ) = (cid:82) Ω ( e ) F d V / Ω ( e ) (cid:46) averaged deformation gradient for each element ( e ) for , . . . , n loop do (cid:46) Jacobi method: repeat n loop timescompute (cid:52) f ( e ) according to (32) (cid:46) Laplace operatorcompute Φ ( e ) according to (35) (cid:46) indicator function if Φ ( e ) > and det ¯ F ( e ) > then update f ( e ) ← f ( e ) − f ( e ) dΦ according to (38) (cid:46) new damage value for inelastic Φ ( e ) if D ( e ) = 1 − f ( e ) > D crit then set D ( e ) = D crit and f ( e ) = 1 − D crit (cid:46) activation of element erosion end ifelse set Φ ( e ) = 0 (cid:46) elimination of elastic Φ ( e ) end ifif max { Φ ( e ) } < − then exit (cid:46) terminate Jacobi method end ifend for the finite element program FEAP, cf. [38]. For the strain energy corresponding to the fictivelyundamaged state, the Neo-Hooke energy given asΨ = µ I −
3) + g ( J )(39)with g ( J ) = λ ( J − − (cid:0) λ + µ (cid:1) ln J and I = tr[ C ], J = det F is employed (see also [37]).The Lam´e constants λ = Eν/ ((1 + ν )(1 − ν )) and µ = E/ (2(1 + ν )) are computed fromthe Young’s modulus E and the Poisson ratio ν . Throughout the numerical tests, an elementerosion approach is incorporated as discussed in Sec. 3 through which, once the damage valueof an element succeeds a critical value D ( e ) > D crit , the element stiffness is replaced by a virtualremaining stiffness s crit and residual forces are eliminated, cf. Alg. 1. Moreover, to avoid damageevolution under compression, we add the condition det ¯ F ( e ) > As first example, we consider a plate with a circular hole which is subjected to the prescribeddisplacement u (cid:63) = (0 , u (cid:63) , T with u (cid:63) = 25 mm at its upper boundary ( Y = L , with L = 100 mm).The geometry of the considered domain is depicted in Fig. 2(a). Due to the symmetry of theproblem, the computational domain amounts only to the upper right quarter of the geometry.11 nitialization • damage function: f = 1 ∀ x ∈ Ω • boundary conditions: u (cid:63) , b (cid:63) , t (cid:63) finite element analysis • update boundary conditions: u (cid:63)n +1 , b (cid:63)n +1 , t (cid:63)n +1 • employ element erosion technique for severe damage, see Alg. 1 • solve nonlinear equilibrium equation for ˆ u with fixated field f ( e ) converged?abort &error report • collect ¯Ψ ( e )0 ,n +1 = ( e ) (cid:82) Ω ( e ) Ψ ( C n +1 ) d V , Ω ( e ) ∀ e = 1 , . . . , n e update of the damage function • update damage function according to Alg. 2 • initialize next time step: t n ← t n +1 t n +1 > t max ?exitno yes yesFigure 1: Flowchart for the proposed damage model.Moreover, as a consequence of the symmetry, the displacement is fixed in X-direction with u X = 0at the left edge X = 0 and in Y-direction with u Y = 0 at the lower edge Y = 0. The radius ofthe circular hole as shown in Fig. 2(a) is given with R = 50 mm and the thickness is given with H = 10 mm. An overview of the used material and load parameters is collected in Tab. 1. In Fig. 3(a), force/displacement curves corresponding to various mesh refinements (see eg.Fig. 2(b)) are shown. The value of the prescribed displacement u (cid:63) is depicted on the abscissawhereas the value of the reaction force F in Y-direction, which is recovered from the nodes atthe upper surface Y = L , is plotted on the ordinate. Here, the nonlocal parameter has thevalue β = 100 N and the prescribed displacement is applied in increments of 0 .
025 mm, whichis equivalent to 1000 load steps. For the given plot resolution in Fig. 3(a), it can be observedthat the force/displacement curve of the lowest refinement stage (400 elements) almost coincideswith the curves corresponding to the higher refinement stages. To illustrate the convergenceof the curves as the number of elements increases, an enlarged resolution of the plot is givenin Fig. 3(b). The area of the enlargement is the transition regime at which material softeningoccurs (marked with gray dashed lines in Fig. 3(a)). The evolution of the material softening isalso illustrated through the contourplots of the damage function D in Fig. 3(c). The contour-plots correspond to the computation with the 6250 element mesh and the depicted stages ofdamage evolution are marked with bullets in Fig. 3(a) and (b). It can be observed that alreadyat low damage propagation states, as depicted in the left contourplot, the damage value D ofseveral elements reaches the critical damage value D crit . For all computations, convergence of12able 1: Material and boundary parameters used throughout the numerical tests.Description Symbol Value UnitPlate with hole U-shapeE-modulus E
500 1000 MPapoisson ratio ν . . r . . β ∈ { , , }
100 Ncritical damage value D crit .
95 0 .
995 -critical stiffness s crit − − N/mmload increment - { . , . , . , . } . XZY R ˆ uL LH (b)Figure 2: (a) Description of the geometry of the plate with hole benchmark problem. (b) Examplefinite element meshes (400 and 6250 elements).13
00 elements1350 elements3200 elements6250 elements u * [ mm ] F [ k N ] β =
100 N (a)
400 elements1350 elements3200 elements6250 elements u * [ mm ] F [ k N ] β =
100 N (b) D . . .
95 = D crit (c)Figure 3: (a) Force/displacement curves (higher resolution in (b)) and (c) corresponding con-tourplots depicting the damage evolution for the 6250 element mesh in three different load stepswhich are depicted as bullets in (a). 14he iterative solution procedure was obtained. In this subsection, the influence of the value of the nonlocal parameter β on the convergencebehavior and the damage evolution is investigated. Therefore, in Fig. 4 corresponding force/dis-placement curves are shown for values β = 10 ,
100 and 1000 N, respectively. The same loadincrement of 0 .
025 mm as in the previous subsection is used. From the plots of Fig. 4, it canbe observed that as β increases the convergence of the curves appears more rapidly. However,as β increases also the starting point of damage initialization is shifted to higher load values.This goes along with the increased nonlocality of the distribution of the damage function D asdepicted in the contourplots of Fig. 5(c): while in the left contourplot ( β = 10 N) the damagedistribution is mostly concentrated to the elements adjacent to the lower edge, in the rightcontourplot ( β = 1000 N), a more smeared damage distribution can be observed. To investigate the influence of the load increment on the solution, we show the convergence of theforce/displacement curves for the specific displacement increments { . , . , . , . } mm,while the mesh refinement stage (400 elements) and nonlocal parameter ( β = 1000 N) remainfixed. For the chosen parameters, the corresponding curves are depicted in Fig. 6(a) and showthe tendency to converge as the value of the load increment decreases and, consequently, thenumber of steps increases. Furthermore, Fig. 6(b) shows the computing time for the gradientdamage model (labeled with GD) and a comparative purely elastic simulation (labeled withEL). The time is measured until the crack has completely evolved which corresponds to a loadof u (cid:63) = 10 .
85 mm. Considering also later loading steps would result in an overestimation of thecomputation time of the elastic simulation: whereas for the elastic material behavior a non-linearfinite element simulation has to be performed which requires more than two iteration steps forconvergence, the damage model has to solve the rigid body movement due to the completelyevolved horizontal crack that separated the plate. Consequently, only two iteration steps arerequired for convergence yielding a total computing time that is larger for the elastic problemthan for the damage model. The elastic problem has been solved by deactivating the updatesubroutine for the damage field. This implies that the same effort for writing the results to textfiles was needed and no unfair acceleration of either simulation has been employed. For reachingthe critical load of u (cid:63) = 10 .
85 mm, the elastic simulation required t = 154 . t = 159 . . As second numerical test, a u-shaped geometry is analyzed in which, caused by its geometry andboundary conditions, large deformations are present (cf. Fig. 7(a)) although the strains may becomparatively moderate. Due to its symmetry, only the left half of the domain is considered andthe displacement u X = 0 is fixed at the symmetry plane X = 0. The dimensions (cf. Fig. 7(a))15a)
400 elements1350 elements3200 elements6250 elements u * [ mm ] F [ k N ] β =
10 N
400 elements1350 elements3200 elements6250 elements u * [ mm ] F [ k N ] β =
10 N (b)
400 elements1350 elements3200 elements6250 elements u * [ mm ] F [ k N ] β =
100 N
400 elements1350 elements3200 elements6250 elements u * [ mm ] F [ k N ] β =
100 N (c)
400 elements1350 elements3200 elements6250 elements u * [ mm ] F [ k N ] β =
400 elements1350 elements3200 elements6250 elements u * [ mm ] F [ k N ] β = Figure 4: Influence of the nonlocal parameter: force/displacement curves for nonlocal parameterswith enlarged resolution on the right hand side. Figures (a), (b) and (c) correspond to parameters β = { , , } N, respectively. 16 . . .
95 = D crit Figure 5: Influence of the nonlocal parameter: contourplots for varying values of β (from left toright) corresponding to the bullet marks in Fig. 4(a)-(c). Left: β = 10 N, middle: β = 100 N,right: β = 1000 N. load increment0.25 mm0.1 mm0.025 mm0.01 mm u * [ mm ] F [ k N ] (a) [ s ] (b)Figure 6: (a) Influence of the load increment on the solution: force/displacement curves forvarying load increments for a 400 element mesh and β = 1000 N. (b) GD: computing time untilreaching the marked point at u (cid:63) = 10 .
85 mm in (a) with the load increment 0 .
025 mm. EL:computing time of a purely elastic 8-node hexahedral reference element until the same load of u (cid:63) = 10 .
85 mm. 17 uYZ XRS S ˆ u (a)
960 elements3240 elements4680 elements u * [ mm ] F [ k N ]
960 elements3240 elements4680 elements u * [ mm ] F [ k N ] (b) D . . .
995 = D crit (c)Figure 7: (a) Description of the geometry of the u-shaped boundary value problem. (b)Force/displacement curves for various mesh refinement stages (enlarged resolution in the bot-tom). (c) Contourplots of the u-shaped problem visualizing the damage evolution correspondingto the marked load stages depicting the damage evolution.18re given with R = 50 mm and S = 10 mm. The prescribed displacement u (cid:63) = ( u (cid:63) , ,
0) with u (cid:63) = 100 mm is imposed on the upper right edge ( X, Y ) = ( − ,
50) mm. The used material andboundary parameters are shown in Tab. 1. From the force/displacement curves shown in theupper plot in Fig. 7(b) and with refined resolution in the lower figure, the tendency to convergecan be concluded. Due to the relatively large displacement, the additional material nonlinearitybecomes visible in the force/displacement curves (pile-up of the reaction force close to failure).The contourplots in Fig. 7(c) illustrate the damage propagation at the load stages marked withbullets in the force/displacement plots. Here, to visualize the crack formation in the entirelydamaged area of the domain, the elements in which the damage value has reached the upperbound D ( e ) = D crit = 0 .
995 were made invisible. A remarkable elastic snapback is observed oncethe geometry fails, cf. Fig. 7(c) 7). Despite the relatively coarse time discretization, convergenceof the iterative solution was obtained at all refinement stages. It is worth mentioning that thesnapback phenomenon constitutes a non-trivial process in a numerical treatment. Since we alsoobserve convergence for this regime, the u-shaped boundary value problem empirically shows thenumerical robustness of our approach to damage modeling for large deformations.
We presented a novel gradient damage model for hyperelastic structures undergoing large defor-mation along with an efficient and stable numerical update scheme. Numerical results showedconvergence for varying mesh size and increasing regularization parameter β . For a smoothiterative solution scheme, it was necessary to develop a stabilization technique that allows forerosion of finite elements with damage exceeding a critical damage state without the need ofremeshing. Concluding, we obtained an efficient numerical approach for damage processes to bedescribed in the large deformation setting. In future investigations, important physical aspectslike plasticity and hardening will be included. Acknowledgment
The authors J. Riesselmann and D. Balzani greatly appreciate funding by the German ScienceFoundation (Deutsche Forschungsgemeinschaft, DFG) as part of the Priority Program German1748 “Reliable simulation techniques in solid mechanics. Development of non-standard dis-cretization methods, mechanical and mathematical analysis”, project ID BA2823/15-1.
References [1] P. Junker, S. Schwarz, D.R. Jantos, and K. Hackl. A fast and robust numerical treat-ment of a gradient-enhanced model for brittle damage.
International Journal for MultiscaleComputational Engineering , 17(2), 2019.[2] Z.P. Baˇzant and M. Jir´asek. Nonlocal integral formulations of plasticity and damage: surveyof progress.
Journal of engineering mechanics , 128(11):1119–1149, 2002.[3] R.H.J. Peerlings, M.G.D. Geers, R. de Borst, and W.A.M. Brekelmans. A critical comparisonof nonlocal and gradient-enhanced softening continua.
International Journal of solids andStructures , 38(44):7723–7746, 2001.[4] S. Forest and E. Lorentz. Localization phenomena and regularization methods.
Local ap-proach to fracture , 1:311–371, 2004. 195] Y. Yang and A. Misra. Micromechanics based second gradient continuum theory for shearband modeling in cohesive granular materials following damage elasticity.
InternationalJournal of Solids and Structures , 49(18):2500–2514, 2012.[6] M. Ortiz and L. Stainier. The variational formulation of viscoplastic constitutive updates.
Computer methods in applied mechanics and engineering , 171(3-4):419–444, 1999.[7] A. Mielke. Energetic formulation of multiplicative elasto-plasticity using dissipation dis-tances.
Continuum Mechanics and Thermodynamics , 15(4):351–382, 2003.[8] K. Hackl, A. Mielke, and D. Mittenhuber. Dissipation distances in multiplicative elasto-plasticity. In W. Wendland and M. Efendiev, editors,
Analysis and Simulation of Multifieldproblems. Lecture Notes in Applied and Computational Mechanics , volume 12, pages 87–100.Springer, 2003.[9] C. Miehe. Strain-driven homogenization of inelastic microstructures and composites basedon an incremental variational formulation.
International Journal for numerical methods inengineering , 55(11):1285–1322, 2002.[10] E. G¨urses and C. Miehe. On evolving deformation microstructures in non-convex partiallydamaged solids.
Journal of the Mechanics and Physics of Solids , 59(6):1268–1290, 2011.[11] D. Balzani and M. Ortiz. Relaxed incremental variational formulation for damage atlarge strains with application to fiber-reinforced materials and materials with truss-like mi-crostructures.
International journal for numerical methods in engineering , 92(6):551–570,2012.[12] T. Schmidt and D. Balzani. Relaxed incremental variational approach for the modeling ofdamage-induced stress hysteresis in arterial walls.
Journal of the mechanical behavior ofbiomedical materials , 58:149–162, 2016.[13] S. Schwarz, P. Junker, and K. Hackl. Variational regularization of damage models based onthe emulated RVE.
Continuum Mechanics and Thermodynamics , pages 1–27, 2020.[14] R. Faria, J. Oliver, and M. Cervera. A strain-based plastic viscous-damage model for massiveconcrete structures.
International journal of solids and structures , 35(14):1533–1558, 1998.[15] A. Needleman. Material rate dependence and mesh sensitivity in localization problems.
Computer Methods in Applied Mechanics and Engineering , 67:69–85, 1988.[16] A. Suffis, Ton A. A. Lubrecht, and A. Combescure. Damage model with delay effect: Analyt-ical and numerical studies of the evolution of the characteristic damage length.
InternationalJournal of Solids and Structures , 40(13-14):3463–3476, 2003.[17] P. Junker, S. Schwarz, J. Makowski, and K. Hackl. A relaxation-based approach to damagemodeling.
Continuum Mechanics and Thermodynamics , 29:291–310, 2017.[18] K. Langenfeld, P. Junker, and J. Mosler. Quasi-brittle damage modeling based on incremen-tal energy relaxation combined with a viscous-type regularization.
Continuum Mechanicsand Thermodynamics , 30(5):1125–1144, 2018.[19] R. de Borst. Some recent issues in computational failure mechanics.
International Journalfor Numerical Methods in Engineering , 52(1-2):63–95, 2001.2020] R.H.J. Peerlings, R. de Borst, A.M. Breckelmans, and J.H.P. de Vree. Gradient enhanceddamage for quasi-brittle materials.
International Journal for Numerical Methods in Engin-geering , 39:3391–3403, 1996.[21] R. Peerlings, T.J. Massart, and M.G.D Geers. A thermodynamically motivated implicitgradient damage framework and its application to brick masonry cracking.
Comp. MethodsAppl. Mech. Engrg. , 193:3403–3417, 2003.[22] B. Dimitrijevic and K. Hackl. A method for gradient enhancement of continuum damagemodels.
Tech. Mech , 28(1):43–52, 2008.[23] D. Gross and T. Seelig.
Fracture mechanics: with an introduction to micromechanics .Springer, 2017.[24] B. Kiefer, T. Waffenschmidt, L. Sprave, and A. Menzel. A gradient-enhanced damage modelcoupled to plasticity—multi-surface formulation and algorithmic concepts.
InternationalJournal of Damage Mechanics , 27(2):253–295, 2018.[25] C. Miehe, F. Aldakheel, and A. Raina. Phase field modeling of ductile fracture at finitestrains: A variational gradient-extended plasticity-damage theory.
International Journal ofPlasticity , 84:1–32, 2016.[26] M. Ambati, R. Kruse, and L. De Lorenzis. A phase-field model for ductile fracture at finitestrains and its experimental verification.
Computational Mechanics , 57(1):149–167, 2016.[27] V. Carollo, J. Reinoso, and M. Paggi. A 3d finite strain model for intralayer and interlayercrack simulation coupling the phase field approach and cohesive zone model.
CompositeStructures , 182:636–651, 2017.[28] M.J. Borden, T.J.R. Hughes, C.M. Landis, A. Anvari, and I.J. Lee. A phase-field formu-lation for fracture in ductile materials: Finite deformation balance law derivation, plasticdegradation, and stress triaxiality effects.
Computer Methods in Applied Mechanics andEngineering , 312:130–166, 2016.[29] O. G¨ultekin, H. Dal, and G.A. Holzapfel. A phase-field approach to model fracture ofarterial walls: theory and finite element analysis.
Computer Methods in Applied Mechanicsand Engineering , 312:542–566, 2016.[30] F. Fathi, S.H. Ardakani, P.F. Dehaghani, and S. Mohammadi. A finite strain integral-type anisotropic damage model for fiber-reinforced materials: Application in soft biologicaltissues.
Computer Methods in Applied Mechanics and Engineering , 322:262–295, 2017.[31] T. Waffenschmidt, C. Polindara, A. Menzel, and S. Blanco. A gradient-enhanced large-deformation continuum damage model for fibre-reinforced materials.
Comp. Methods Appl.Mech. Engrg. , 268:801–842, 2013.[32] T. Brepols, S. Wulfinghoff, and S. Reese. Gradient-extended two-surface damage-plasticity:micromorphic formulation and numerical aspects.
International Journal of Plasticity , 97:64–106, 2017.[33] A. Vogel and P. Junker. Adaptive and highly accurate numerical treatment for a gradient-enhanced brittle damage model.
International Journal for Numerical Methods in Engineer-ing , 121(14):3108–3131, 2020. 2134] D.R. Jantos, K. Hackl, and P. Junker. An accurate and fast regularization approach tothermodynamic topology optimization.
International Journal for Numerical Methods inEngineering , 117(9):991–1017, 2019.[35] P. Junker and D. Balzani. A new variational approach for the thermodynamic topologyoptimization of hyperelastic structures.
Computational Mechanics , pages 1–26, 2020.[36] P. Junker and D. Balzani. An extended hamilton principle as unifying theory for coupledproblems and dissipative microstructure evolution. doi: arxiv.org/abs/2010.09440, 2020 .[37] P. Wriggers.
Nonlinear finite element methods . Springer Science & Business Media, 2008.[38] FEAP. A finite element analysis program. http://projects.ce.berkeley.edu/feap/http://projects.ce.berkeley.edu/feap/