An Efficient Treatment of the Laplacian in a Gradient‐Enhanced Damage Model

As is well‐known, softening effects that are characteristic for damage models are accompanied by ill‐posed boundary value problems. This arises from non‐convex and non‐coercive energies and results in mesh‐dependent finite element results. For that reason, regularization strategies, which somehow take into account the non‐local behavior, have to be applied in order to prevent ill‐posedness and to achieve mesh‐independence. Hereto, most commonly gradient‐enhanced formulations [1, 2] are considered, but also integral‐type [3,4] and viscous [5,6] regularization are well‐known.


Material Model
The underlying Helmholtz free energy Ψ, describing the isotropically damaging material, is given by The individual evaluation of δ u H provides the stationarity condition for the displacement field u as follows Moreover, the individual evaluation of δ d H provides the stationarity condition for the internal variable d including a set-valued sub-differential ∂D handling the non-differentiability, thus Integration by parts of the second term yields the Laplace operator ∇ 2 according to Herewith, the strong form of the stationarity condition of d including the Neumann boundary condition is obtained by In order to achieve an indicator function Φ deciding on the damage evolutionḋ, a Legendre transformation is performed yielding Φ := p − r ≤ 0 on the assumption of no healing by sgnḋ = {0, 1} and including the driving force p derived by Finally, the damage evolution is completely described bẏ

Algorithmic Procedure for Displacement Field
The displacement field u can easily be solved by means of the respective stationarity condition for δ u and application of the finite element method. According to typical procedure, shape functions S and the B-operator matrix B are applied as follows A first operator split is introduced by using the damage functions f n from the previous time step n, whereby their current values are uncoupled from the calculation of the displacement field. Thus, the stresses σ are only depending on the current displacements and the tangent K becomes constant for each time step. A simple and fast numerical handling follows by but a minimum number of loading steps is required to overcome the influences due to the operator split.

Algorithmic Procedure for Damage Function
The damage function f i is defined at the center of mass of each element i and thus the indicator function Φ i is element-wise formulated including the respective element-wise averaged elastic energyΨ 0,i with volume Ω i . Therefore A second operator split is introduced by using the neighbored damage functions f k i from the previous iteration step k within the Newton scheme. This simplifies the system of coupled inequalities with a coefficient matrix n e × n e into uncoupled inequalities of number n e . Again, a minimum number of loading steps is required. The indicator function Φ i becomes At this point, the Taylor series expansion for a central element x i and nine neighbors x j (in 3D, five in 2D) is presented by considering the spatial increments ∆x j := x i − x j between the central and neighbored elements. The linear system of equations by ∆f i = D · ∂f i includes the 9×1 vector of known increments of the damage function ∆f i = (∆f i,j )e j , the 9×9 matrix of known spatial increments D with the following entries for each neighbor j D j,1 = ∆x j,1 , D j,4 = ∆x j,1 ∆x j,2 , D j,7 = (∆x j,1 ) 2 /2 , D j,2 = ∆x j,2 , D j,5 = ∆x j,1 ∆x j,3 , D j,8 = (∆x j,2 ) 2 /2 , D j,3 = ∆x j,3 , D j,6 = ∆x j,2 ∆x j,3 , D j,9 = (∆x j,3 ) 2 /2 , and the 9×1 vector of unknown partial derivatives Solving for the unknown partial derivatives ∂f i = D −1 · ∆f i and introducing the constant 3×9 matrix B ∇ 2 allows to capture all necessary unmixed second derivatives for the Laplace operator Λ f and its derivative DΛ f , consequently Finally, the damage evolution for

Numerical Results
Finite element calculations are performed in order to investigate the behavior of the model with respect to different loading rates as well as discretizations. Hereto, the double-notched plate (see right-hand side) is loaded by a displacement u = 0.015 mm and two sets of regularization are considered: a slight one with β = 0.002 MPa mm 2 for a localized and a stronger one with β = 0.015 MPa mm 2 for a diffusive characteristic. Further parameters are E = 210 GPa, ν = 0.33 and r = 1 MPa. The results are provided by Figure 2 presenting the zoomed damage distribution by means of the damage function f (d) and by Figure 3 for the corresponding global structural response by means of force-displacement diagrams.
The variation of the loading rate is realized by different number of time steps. Both the damage distribution as well as the global structural response demonstrate a converging behavior of the results with an increasing number of time steps. Obviously, in case of 1000 time steps hardly any effects due to the operator splits are included anymore. Consequently, this loading rate is considered to be suitable for further investigations. Accordingly, the variation of the discretization is realized by different number of elements. The results certify a mesh-independent behavior of the calculations by very similar damage distributions and wellagreeing force-displacement curves. Besides, these curves are characterized by very steep drops which fall to zero. This fall occurs more quickly in case of the stronger regularization described by the diffusive characteristic.

Conclusion
A novel approach for gradient-enhanced damage modeling based on a very efficient evaluation of the Laplace operator is presented. For the purpose of efficiency, two operator splits are applied and possible influences on the results are prevented by application of a suitable loading rates as investigated. Numerical examples for localized as well as diffusive damage characteristics provided mesh-independent results by very similar distributions of the damage function and well-agreeing global structural responses. Also, the computation times are very fast and for each load step close to purely elastic calculations.