From ductile damage to unilateral contact via point-wise implicit discontinuity at the infinitesimal element level

Ductile damage models and cohesive laws incorporate the material plasticity entailing the growth of irrecoverable deformations even after complete failure. This unrealistic growth remains concealed until the unilateral effects arising from the crack closure emerge. We address this issue by proposing a new strategy to cope with the entire process of failure, from the very inception in the form of diffuse damage to the final stage, i.e. the emergence of sharp cracks. To this end, we introduce a new strain field, termed discontinuity strain, to the conventional additive strain decomposition to account for discontinuities in a continuous sense so that the standard principle of virtual work applies. We treat this strain field similar to a strong discontinuity, yet without introducing new kinematic variables and nonlinear boundary conditions. In this paper, we demonstrate the effectiveness of this new strategy at a simple ductile damage constitutive model. The model uses a scalar damage index to control the degradation process. The discontinuity strain field is injected into the strain decomposition if this damage index exceeds a certain threshold. The threshold corresponds to the limit at which the induced imperfections merge and form a discrete crack. With three-point bending tests under pure mode I and mixed-mode conditions, we demonstrate that this augmentation does not show the early crack closure artifact which is wrongly predicted by plastic damage formulations at load reversal. We also use the concrete damaged plasticity model provided in Abaqus commercial finite element program for our comparison. Lastly, a high-intensity low-cycle fatigue test demonstrates the unilateral effects resulting from the complete closure of the induced crack.


Introduction
The last decades have witnessed extensive studies on computational failure mechanics.With ever growing computational resources at our disposal, various advanced numerical tools representing different aspects of failure problems have emerged, enabling us to predict complex failure scenarios.The complexity arises from the variety of stages involved, comprising the presence of diffuse imperfections, localization of intense deformations, and the formation of distinct sharp cracks.There is a spectrum of options to cope with the failure process, ranging from purely continuous-based approaches to purely discontinuous ones [1].On the one hand, continuous approaches span the damaged area over a finite width of the medium from the emergence of induced defects, up to the formation of discontinuity surfaces.Therefore, they recast discrete cracks into narrow bands of highly localized deformation with (almost) zero stiffness.Continuum damage models [2][3][4] and phase field methods [5][6][7][8][9][10] belong to this group.On the other hand, discontinuous approaches represent the fracture process zone in the form of displacement jumps, no matter if there truly exists an abrupt change in the displacement field or not.Some prominent tools belonging to this group are the embedded discontinuity model [11][12][13] and the extended finite element method [14][15][16].
Neither of the aforementioned strategies is self-sufficient in dealing with the entire fracture process.Continuous approaches smear sharp cracks over a finite width of the medium, while discontinuous approaches sample the damaged region on the fictitious crack faces.Both represent the two extreme ends of the spectrum of physical effects well and their strengths and deficiencies complement each other to some extent.Therefore, numerous studies have been devoted to leveraging both approaches by establishing a transition from distributed damage to strong discontinuities.Simone et al. [17] alleviated the unrealistic damage growth resembled by the class of continuous models through combining a regularised continuum framework with a partition of unity finite element method.Introducing the thick level set approach, Moës et al. [18] allowed a straightforward transition from damaged zone to complete fracture once the material is totally damaged.Han et al. [19] proposed a continuous-discontinuous framework considering large deformation kinematics to recast diffuse degraded topologies into sharp crack paths.Dynamic crack propagation induced by gradient-enhanced damage growth was presented by Sun and Löhnert [20].They coupled transient thermal and dynamic mechanical fields and, used the extended finite element method in three-dimensional settings for representing discrete cracks.Pandey et al. [21] introduced a hybrid methodology by incorporating the extended finite element method and continuum damage mechanics to represent creep-fatigue crack propagation.A continuous-discontinuous strategy within which a thin layer represents the discontinuity emerging due to localized failure is presented by Puccia et al. [22].This layer behaves similar to the bulk material so that no additional constitutive model associated with cohesive-like approaches is required.Regarding the conservation of energy during the transition, some continuous-discontinuous approaches inject the discontinuity at the final stage of the fracture process when the material is almost fully degraded to avoid spurious energy release.The combined model of Seabra et al. [23] for damagedriven crack propagation in ductile metals, the fracture-based continuous-discontinuous approach of Sarkar et al. [24], and the enriched continuum model of Negi and Kumar [25] fall into this category.Others use the concept of energy equivalence to define an intermediate state so that the transition can be triggered at any stage in a smooth manner.The hybrid model of Cuvilliez et al. [26], the continuum to discontinuum transition strategy of Roth et al. [27], and the thermodynamically consistent model of Wang and Waisman [28] are some examples of this type.
All of the above and numerous other can accurately address the failure process over moderately complex loading scenarios.Yet they are not necessarily well-equipped to deal with complex loading paths.The most overlooked aspect of the fracture process in the strategies mentioned is the material response during load reversal.Elegant constitutive models and cohesive laws incorporate plasticity by which irreversible strains and displacement jumps are coupled with reversible ones [29].As a result, all continuous, discontinuous or hybrid models exhibit excessive unrealistic permanent deformations by construction which are accumulated within the smear or on the discrete fracture process zone.Figure 1 depicts a generic stress-strain response and traction-separation law of such models.The crack closure effects lead to stiffness recovery upon transition from tension to compression.This unilateral behavior is essential in general loading scenarios, such as cyclic loading [30][31][32][33], fatigue tests [34][35][36][37], and seismic modeling [38][39][40], that entail complex strain paths.Obeying the responses shown, fully degraded regions develop excessive permanent deformations that shift the transition point to the right of the curve.As a result, the unilateral effects may get triggered in a faulty state and lead to artificial stiffening.This is elaborated further in Figure 2, which shows the normal stress distribution along the crack axis.According to the figure, the fracture process zone nucleated in front of the crack tip spans over the domain that experiences a nonlinear response.As a consequence of balancing the internal forces, part of this zone experiences compressive stresses upon external load removal (see Figure 2 (b)).However, the crack faces spanning from the notch tip to the crack tips remain traction-free in this relaxed state.This is not the case for the relaxed numerical model as the fully degraded section of the crack develops compressive stresses that affect the global response of the numerical model.This artificial stiffening remains untraced when a constantly increasing load is applied and can only be detected in scenarios that involve unloading.A remedy is to limit the growth of permanent deformations to an ultimate state corresponding to the limit at which discrete crack faces emerge.This treatment is completely impractical in the continuous models, yet can be applied in the discontinuous ones by permanently detaching the opposing crack faces and removing the stress continuity condition.In this case, resembling the unilateral effects still invokes the use of computational contact mechanics, which renders the model even more complex.
By contrast, the paper at hand introduces the discontinuity in the strain field to mimic discrete crack faces.As a result, no nonlinear boundary conditions or kinematic variables are involved, preserving the standard form of the virtual work principle.To this end, we organize the body of this text as follows.After this introduction, the additive strain decomposition including the discontinuity strain field is introduced in Section 2, followed by the underlying theory of a ductile damage model involving the unilateral effects.It is worth mentioning that a simple damage model is chosen to avoid unnecessary complications arising from the constitutive modeling.However, the discontinuity strain field can be incorporated into other material models regardless of the hypothesis involved.Section 3 is devoted to numerical examples.A three-point bending specimen under pure opening mode is analyzed first.Then, a similar test is conducted under mixed-mode conditions.The complete closure of the crack is also resembled by applying a high-intensity low-cycle load on a double-edge notched specimen before conclusions drawn in Section 4 close the paper.

Theory
Material inelasticity involves both reversible and non-reversible changes of shape in the consequence of applied loads.The classical models of infinitesimal plasticity characterize these two by means of an additive split with which the total strain is decomposed into an elastic and a plastic part.We augment this additive strain decomposition by introducing the discontinuity strain field ε d so that the total strain tensor reads where ε e and ε p are the conventional elastic and plastic parts, respectively.Here, the strain field ε d represents the induced crack due to the excessive damage growth, but in general, it can mimic any kind of discontinuity.In contrast to the plastic strain ε p , the discontinuity strain field is reversible, allowing it to resemble the opening and closing of the hypothetical crack.Figure 3 shows an abstraction of this additive decomposition.By activating the discontinuity strain field at a certain damage threshold, which corresponds to the limit of a sharp crack forming due to the complete detachment of material, further plastic straining is prevented.As a consequence, crack closure can occur at this detachment limit in the case of load reversal.Here, we incorporate this strain split in a simple isotropic damage model to keep the constitutive modeling minimal.However, it is worth mentioning once again that the discontinuity strain field can be used in any material model and also for representing any kind of discontinuity.Following the concept of effective quantities originally introduced by Kachanov [41], we characterize the material elastoplastic regime within a fictitious intact configuration at which the mechanical stress, known as the effective stress σ, is the measure of average force acting on the undamaged cross-sectional area of a body.Hence, whether the infinitesimal volume element of material is influenced by damage or not, the following general constitutive law of isotropic material applies or, equivalently, σ = D : (ε where D is the tangent operator of intact material.Note that we preserve this naming convention throughout the text so that the tilde accent refers to the effective configuration.
According to the work of Feenstra and de Borst [42], the yield locus is given in the effective configuration by means of the Rankine maximum principal stress criterion through where σy is the yield strength.Hence, the admissibility of the stress state must be preserved by enforcing it to remain inside or on this locus.To this end, the evolution of the plastic strain is given by the associated flow rule wherein γ is a Lagrange multiplier.The pair f and γ is subjected to the Kuhn-Tucker optimality conditions which enforce the stress state admissibility determining whether yielding occurs or not.Now by defining the stress tensor in the effective configuration can be mapped to its macroscopically observed counterpart, the Cauchy stress tensor σ.In the above mapping, the damage index d represents the density of imperfections reducing the load-bearing area of the infinitesimal volume element of the material.However, in order to establish the unilateral conditions, the above relation must be modified such that the damage index d only applies if the stress state is in tension.  is defined wherein the subscripts t and c denote the tensile and compressive parts, respectively [43].The tensile part can be given by [44] σt = wherein is the Macaulay bracket, and σi and ei are the i th principal value and principal direction, respectively.Employing the additive split (8), the compressive part reads With the tensile and compressive parts at hand, the stress mapping including the unilateral effects is given by As a result, due to the micro-cracks closure upon transition from tension to compression, the initial stiffness of the material is recovered.The damage evolution is closely linked to the growth of plastic strain [45].However, depending on the material being studied, its growth may start as the plastic strain is mobilized or when the yielding process reaches a certain stage.Inspired by the work of Grassl and Jirásek [46], we use the exponential law where a is considered as a material constant identifying the magnitude of damage growth, and is the accumulated equivalent plastic strain.It is worth mentioning that, for a model aiming at realistic constitutive modeling, a more complicated relation can be opted for.Now, at the onset of cracking we have where dcr is the critical damage denoting the threshold at which the induced imperfections merge and form a sharp crack.
Utilizing the above expression, we can define the critical accumulated equivalent plastic strain at the inception of cracking as Thus, the discontinuity strain field mobilizes if the accumulated equivalent plastic strain reaches the above limit.Referring to (14), for two arbitrary successive (pseudo-)time steps t and t + dt we have However, the accumulated equivalent plastic strain cannot exceed the limit εp cr .Hence, if the standard return mapping causes an overshoot, there exists an intermediate state between these two successive steps, identified by t + αdt, at which and the maximum principal stress hits the yield locus at the same time.Defining λ that denotes the evolution of εp from the time step t to the time step t + αdt, we arrive at By plugging the magnitude of εp t+αdt from (18) into the above expression, λ is obtained as As a result, the growth of the equivalent plastic strain is known beforehand and the evolution of the plastic strain from t to t + αdt can be given explicitly by εp = λ∂σf.
Hence, instead of solving the standard nonlinear return mapping, we must find α such that From the intermediate state t + αdt onward, further strains are accumulated on the discontinuity strain field.As a result, the discontinuity strain ε d emerges at t + αdt and receives the remaining part of the total strain rate during t + αdt to t + dt so that In addition, to comply with the Rankine maximum principal stress criterion, the discontinuity plane is set to the plane of maximum tensile stress at the intermediate state t + αdt.The normal vector to this plane, denoted by n, will be used in the subsequent stages of the analysis to check whether the crack is open or not.Crack closure is detected if meaning that the opposing crack faces penetrate each other.Once again, we must find an intermediate step, denoted by t + βdt, at which the opposing crack faces meet.By writing the incremental evolution of the discontinuity strain field from t to t + βdt, we arrive at Solving for β, the exact time of crack closure is obtained.After crack closure, the discontinuity strain field vanishes and the standard elastoplasticity applies until the inception of another crack opening regime.It should be noted that, by employing this formulation, the maximum damage that can be reached is limited to its critical value since no plastic straining occurs afterward.As a result, the degradation process freezes once the damage index reaches its upper limit, which is only acceptable if dcr is chosen extremely close to unity.Otherwise the material preserves its load-bearing capacity and never degrades completely.To allow further damage growth, the exponential function ( 13) is redefined as where b is a dimensionless constant used to adjust the damage growth after cracking, and εd = max is the maximum strain jump that is experienced during the loading history.Needless to say that a proper strategy is to establish a relation between the damage growth function (28) and the fracture energy of material so that the constants a and b get linked to a physical property.In addition, it can provide the objectivity of global responses with respect to numerical discretization.However, since this paper is not aimed at the constitutive modeling of material, we will pursue this target in a separate paper in the future.

Numerical results
This section is devoted to showing the contribution of the discontinuity strain field in the global responses of specimens undergoing the failure process.To this end, three laboratory tests are chosen, including a pure mode I three-point bending test, a mixed-mode test, and a double-edge notched specimen subjected to cyclic loading.The standard finite element method is used to discretize the domain of problems.Note that the solutions are not unique since the softening response during the failure process induces a local material instability that leads to the ill-posedness of the boundary value problem.This can be cured by injecting a measure of length, known as the material length scale, representing the width of the fracture process zone through non-local or gradient-enhanced formulations or by using the concept of fracture energy equivalence [47].In the paper at hand, this ill-posedness is left untreated to focus only on the introduced discontinuity strain field and keep the message of the paper as clear as possible.As a result, due to the unregularized nature of the field equations, the damage constants a and b are calibrated in accordance with the finite element meshes.This mesh-dependency vanishes once a material length scale is introduced.We analyze each test using different assumptions of the critical damage dcr to show its effects on the unloading branches.To consider the cases in which the conventional strain decomposition is used, the critical damage dcr is set to unity.As a result, since d is defined by means of an exponential function that never reaches this ultimate value (see the damage growth function in (28)), the discontinuity strain field is simply neglected in the simulations.In addition, to further illustrate the shortcoming of the conventional formulations, the concrete damaged plasticity model of Abaqus finite element program is also used to simulate the tests.This model is capable of reproducing the unilateral effects, which is the key ingredient for revealing the artificial stiffening upon load reversal.Note that one-dimensional stress-strain and damage-strain curves are required to define the post-linear behavior in the Abaqus model.These curves are provided for each example using the calibrated values of E, σy, and a in conjunction with the one-dimensional Hooke's law and the exponential damage growth function in (13).Other parameters of the model are set to their default values.

Opening-mode test
This numerical example is a replication of the laboratory test of Perdikaris and Romeo [48].A series of concrete beams with similar dimensions, yet notched at different offsets were tested.We use the mid-span notched beam for testing the model under a pure mode I loading.Figure 4 shows the finite element model of the beam on the test setup.The dimensions s, h, and e are 304.8,76.2, and 25.4 millimeters, respectively.In addition, the notch depth is one-third of the height of the beam.The test is simulated under plane stress conditions with a thickness of 28.6 millimeters.The material properties of the beam are reported in Table 1.Four cases of dcr are considered, including 1.00, 0.85, 0.60, and 0.45.The former, dcr = 1.0, denotes the condition in which the discontinuity strain field is disregarded.Decreasing the value of dcr, the contribution of the discontinuity strain field becomes more pronounced.The test is also modeled by means of the concrete damaged plasticity model of Abaqus.
The global responses of the numerical models are plotted against the experimental data in Figure 5.The horizontal axis, i.e. the crack mouth displacement (CMD), represents the absolute relative displacement of the crack mouth.Figure 5 (a) clearly shows that the conventional additive strain decomposition, which is the case in the model with dcr = 1.0 and the constitutive mode of Abaqus, causes the unloading branches of the curves to have a slope similar to the initial elastic  branch.On the other hand, by decreasing the value of dcr, which is accompanied by earlier activation of the discontinuity strain field, the unloading slopes gradually become less steep.The cause of this behavior is elaborated further by comparing the stress distributions before and after the load removal.To this end, consider the loaded and the corresponding relaxed states of the Abaqus model and the dcr = 0.45 case depicted in Figure 6.The stress contours of these cases are shown in Figure 7.Note that the absolute values of the principal stresses are used so that both tensile and compressive distributions are shown.Figure 7 (a) depicts the stress distributions at the beginning of the unloading phase for the Abaqus model and Figure 7 (c) for the case of including the discontinuity strain field with the critical damage dcr = 0.45.The domain occupied by the fully formed crack is almost stress-free in both cases and the crack front experiences tensile stress with a peak value identical to the yield strength.However, in the relaxed state of the Abaqus model (Figure 7 (b)), the notch tip experiences a considerable amount of compressive stress.Hence, a transition from compression to tension exists by moving along the crack axis from the notch tip to the crack front.However, the material must be completely detached and remain stress-free along some part of this axis.This faulty behavior is caused due to the unrealistic growth of plastic strain within the fully degraded regions.This excessive plastic strain causes an early activation of the unilateral effects upon a slight reduction of the total strain, leading to an unrealistically early transition from tension to compression.On the other hand, according to Figure 7 (d), by including the discontinuity stain field, some part of the fully damaged region remains stress-free after unloading, and the core of the compressive region is shifted from the notch tip to the middle of the crack.The damaged region located between the tensile and compressive cores can be interpreted as the active fracture process zone, and these two cores indicate the head and tail of this zone.Referring to the schematic stress distributions shown previously in Figure 2, we can deduce that by considering the discontinuity strain filed, the fully detached section of the crack is reproduced within the model, which is not the case when the discontinuity strain field is not included.

Mixed-mode test
The second example is also selected from the work of Perdikaris and Romeo [48].The test setup is similar to the previous one, with the exception that the notch is offset by 75.6 millimeters from the mid-span.The properties presented in Table 2 are used for the numerical model.Similar to the previous test, in addition to the Abaqus model, four cases of dcr = 1.00, dcr = 0.85, dcr = 0.60, and dcr = 0.45 are considered.The applied force versus the crack mouth displacement curves, plotted on the experimental data in Figure 9, conspicuously demonstrate the contribution of the discontinuity strain field.Referring to the stages marked in Figure 10 for the Abaqus model and the dcr = 0.45 case, a deduction similar to that of the previous test regarding the stress contours shown in Figure 11 can be drawn.Upon the external load removal, the fully degraded part of the crack develops compressive stress in the Abaqus model, while this part remains stress-free in the dcr = 0.45 case.

Full-cycle test
The final example is dedicated to showing the crack closing and reopening behavior occurring in full-cycle tests.For this purpose, the high-intensity low-cycle fatigue test of Reinhardt [49], conducted on double-edge notched specimens, is chosen.The test setup and finite element mesh of its numerical model are depicted in Figure 12.As shown in the figure, the length and width of the beam are 250 and 60 millimeters, respectively.The thickness of the beam is also 50 millimeters, included in the numerical model by means of the plane stress assumption.The test is simulated by imposing prescribed opposing displacements on the left and right edges of the specimen.The recorded reaction forces are divided by the cross-section area of the specimen to represent the average stress.The average relative displacement of two vertical lines that are offset by 17.5 millimeters from the mid-section of the specimen is also recorded to generate the stress-deflection curves.
The contribution of the discontinuity strain field is assessed using four different values of d cd , including 1.00, 0.60, 0.40, and 0.20.The concrete damaged plasticity model of Abaqus is also utilized to simulate the test.The material properties employed in the numerical model are presented in Table 3.The stress versus deflection responses of the numerical models are plotted on the experimental data in Figure 13.No sign of stiffness degradation can be observed in the unloading branches of the Abaqus model and the d cd = 1.00 case.On the other hand, by decreasing the value of dcr, the contribution of the discontinuity strain field, which reveals itself in the   activation of the unilateral effects under less permanent deflection, becomes more pronounced.

Conclusion
This paper introduces a new term to the conventional additive strain decomposition to prevent the excessive unrealistic growth of irrecoverable deformations to reproduce the entire failure process of solids.The issue of excessive growth of irrecoverable deformations is rarely addressed (see for example [50,51]) since it only manifests itself if the unilateral effects resulting from to the closing and reopening of fully formed cracks are of interest.However, this is the case in the area of low cycle fracture/fatigue as caused by earthquakes in civil engineering structures, for example.We demonstrated the problem arising from the mentioned excessive strains which are inherent to other plastic damage models and showed how these strains cause a faulty stress redistribution upon load reversal.
We called the newly introduced term discontinuity strain field since it mimics a discontinuity, yet it is defined at the infinitesimal element level.We demonstrated that introduction of the discontinuity stain field to the conventional additive strain decomposition of the plasticity theory avoids the built up of excessive strains by means of three numerical examples: a pure mode I cracking, a mixed-mode cracking, and a high-intensity low-cycle fatigue test.It is worth mentioning that, although we used the discontinuity strain field to resemble the unilateral effects arising from the collision of the opposing crack faces, it can be used to represent any kind of discontinuity.

Figure 1 :
Figure 1: Constitutive response of generic ductile damage models and cohesive laws: (a) stress-strain response, and (b) traction-separation law.

Figure 3 :
Figure 3: An abstraction of the additive strain split including the discontinuity strain field.

FFigure 4 :
Figure 4: Geometry, boundary conditions, and finite element mesh of the opening-mode test.

Figure 5 :
Figure 5: Comparison of the global responses with the experimental data for the opening-mode test: (a) Abaqus model and d cr = 1.00,(b) d cr = 0.85, (c) cr = 0.60, and (d) d cr = 0.45.

Figure 6 :
Figure 6: Loaded and relaxed states at the beginning and the end of the second unloading phase for the openingmode test: (a) Abaqus model, and (b) d cr = 0.45.

Figure 7 :
Figure 7: Stress contours (MPa) at the loaded and relaxed states for the opening-mode test: (a) loaded state in the Abaqus model, (b) relaxed state in the Abaqus model, (c) loaded state in the d cr = 0.45 case, and (d) relaxed state in the d cr = 0.45 case.

Figure 11 :
Figure 11: Stress contours (MPa) at the loaded and relaxed states for the mixed-mode test: (a) loaded state in the Abaqus model, (b) relaxed state in the Abaqus model, (c) loaded state in the d cr = 0.45 case, and (d) relaxed state in the d cr = 0.45 case.

Figure 12 :Figure 13 :
Figure 12: Geometry, boundary conditions, and finite element mesh of the full-cycle test.

Table 1 :
Parameters of the opening-mode test.

Table 2 :
Parameters of the mixed-mode test.

Table 3 :
Parameters of the full-cycle test.