TRANSIENT ANALYSIS OF GEOMETRICALLY NON-LINEAR TRUSSES CONSIDERING COUPLED PLASTICITY AND DAMAGE

The fracture of ductile materials is usually preceded by considerable levels of plastic strain. After a first stage in which plastic strain strengthens the material due to the introduction and increase of dislocations (hardening), a degradation phenomenon begins to take place due to the nucleation of microcracks and microvoids. The nucleation, growth and coalescence of these defects can be modeled using the concepts of Continuum Damage Mechanics. In this theory, a continuous damage variable is defined, which evolves coupled to plastic strain until attaining a critical value associated to rupture. Several damage models have been proposed for ductile materials, two of the most important being attributed to Gurson and to Lemaitre. This work evaluates the effect of Lemaitre’s damage model when applied to 3D trusses subjected to geometrical nonlinearities including inertial forces. To this end, two different damage evolution laws found in the literature are studied and compared. Special attention is given to unstable problems such as snap-through, in which results show that the effect of inertial forces is predominant. Furthermore, it becomes evident that for a realistic description of damage evolution and failure prediction, a different treatment must be given to tensile and compressive states. The work is closed by the discussion of damping effects on the damaged dynamic problem. It should be remarked that the evaluation of coupled plasticity and damage including geometrical nonlinearities, inertial forces and damping is a complex problem. Hence, setting these phenomena in a simple 3D truss framework makes it possible to focus on the description of physical behavior rather than on element technology complexities. This allows a clear understanding of some effects of Lemaitre’s damage, paving the way for implementations using 2D and 3D continuum finite elements.


INTRODUCTION
The analysis of three-dimensional trusses subjected to large displacements is geometrically nonlinear.In addition, often the consideration of finite strains no longer resides in the linear range of material response, because of plasticity and hardening effects.Yielding affects the propagation and growth of microcracks and microvoids, characterizing an associated damage phenomenon.
Several studies on Continuum Damage Mechanics have been conducted in the last decades [11,7], until the development of the principle of equivalent deformation [9], which has led to a wide spectrum of applications in engineering.

ONE-DIMENSIONAL ELASTOPLASTICITY
The elastoplastic rheological model adopted in this work is shown in Fig. 1 [13].This is a frictional device with unit length and area in which the yield stress σ y > 0.
Figure 1.One dimensional friccional model, rate independent plasticity.When applying external forces, there is an associated strain, which can be split into two parts: one purely elastic and the other one purely plastic.The strain measures are work conjugate to associated stress measures, and a list of Hill's tensors can be found in the work developed by Driemeier et al. [5].A stretch ratio λ is defined by Eq. 1, where l and L are the current and initial lengths of the bar, respectively.When l = L, λ = 1 and the bar is unstrained.
A multiplicative decomposition can be defined, for elastic and plastic stretch contributions (Eq.2).The application of logarithmic strain leads to an additive decomposition, as shown by Eq. 3), similar to the additive relation used for small strains.
lnλ = lnλ e + lnλ p = e + p (3) The stress-strain relation for the elastic behaviour is given by Eq. 5.
Plastic strains will only occur if the associated stress is larger than the yield stress.A yield function defined by Eq. 6 must be satisfied.Or else, the material is being subjected to an irreversible state that usually requires the calculation of a new yield stress.
A plastic strain rate is defined, which is function of a plastic multiplier and of the sign of the stress (Eq.7).The Kuhn-Tucker conditions are then defined by Eq. 8.A consistency condition is also defined, which states that when an increment in plastic deformation takes place, there is a physical requirement of nonzero plastic strain rate, while the stress state lies within the yield surface defined by Eq. 6 [13].

Isotropic hardening
The effects of isotropic hardening are introduced by means of a plastic modulus K and a dissipative variable α, resulting in linear hardening.The following considerations are stated for the isotropic hardening elastoplastic model [13]: • The hardening is isotropic, in the sense that there is no displacement of the yield surface center; • The linear behaviour is independent of the load direction.
The yield function is then defined by Eq. 10, where σ Y > 0, K ≥ 0 is the isotropic hardening parameter (plastic modulus), and α : [0, T ] → is a non-negative function of the plastic strain rate, which is called internal hardening variable (Eq.11), According to previous definitions by Eq. 7, one has Figure 2 shows the behaviour of the isotropic hardening in a traction-compression cycle.Considering points 1, 2, 3, and 4, it can be noticed that it is a closed cycle and the final plastic strain is zero.However, the accumulated plastic strain α increases its value along with the plastic strain, independent of the loading direction.Stress strain relationship for the elastic state: σ = E( − p ) Isotropic hardening evolution law: Box 1: One dimensional elastoplasticity with isotropic hardening.
The consideration of isotropic hardening results in a nonzero plastic tangent modulus.Eq. 14 is obtained by Eq. 13 applying partial derivatives of the yield function with respect to the stress and the hardening parameter α (where R = Kα), and describes the relation between strain and stress rates.The detailed deduction of these equations can be found in several sources [1,10,13].

Return-mapping algorithm
The incremental implicit form of the equations for the elastoplastic behaviour is presented by Eqs 16-22.The indexes n and n + 1 are the previous and current load steps, respectively.
1. Given p n , α n (from initial state or last load step); 2. Compute the strain field n+1 = n + ∆ n ; 3. Elastic predictor and yield function axiom verification: Algorithm 1 shows the return-mapping procedure.It calculates an elastic predictor for the stress state from the strain field, giving information if the predicted state is elastic or plastic.If the state is within the elastic field, no hardening is associated, and there is no need to modify the hardening variables.On the contrary, if the predicted stress exceeds the yield surface, it means that the material is in the plastic domain, and the return-mapping procedure must be performed to calculate updated values for the hardening parameters.Note that for a linear hardening model, the algorithm is not iterative, that is, expressions in step 4 have a closed form.

CONTINUUM DAMAGE MECHANICS
The damage approach of this work is based on Lemaitre's formulation [8] for ductile materials.It defines a damage variable based on the influence of microcracks and microvoids on the macroscopic material properties such as the elasticity modulus and the yield stress.
The damage variable is defined by a relation between the original and effective areas on a representative volume element (RVE) [6], shown in Fig. 3.The size of this element varies depending on the homogeneity of the material, and must be small enough to be considered a point in continuum mechanics, and large enough to correctly represent the material properties on a mesoscale.The damage relation is given by Eq. 23.
where Ā(n) is the effective cross-section area of the volume element, and A (n) is the original area (undamaged state) both measured in a plane normal to the unit vector n.D (n) is the damage with dependency on the normal unit area direction.For isotropic damage, there is no such dependency and the damage variable is defined by Eq. 24, where A D is the damaged area.The damage domain then must be within the range 0 ≤ D ≤ 1, so that the equivalent stress is given by Eq. 25. Figure 4 shows a damaged bar, and its equivalent model after removing the cracks and voids, and substituting the equivalent stress.The damage variable can be expressed in terms of the elasticity modulus, using the deformation equivalence principle, as show by Eq. 26.This is called phenomenological damage, defining the influence of the damage variable on the macroscopical properties of the material.
The law of damage evolution is given by Eq. 27.Ḋ is the damage rate, S is a plasticdependent parameter, ṗ is the effective plastic deformation rate (Eq.28), and Y (Eq.29) is the damage energy release rate, which represents the variation of internal energy density due to damage growth at constant stress.
Taking the partial derivatives of the yield function (Eq.30), now with respect to the stress, hardening variable and the damage variable, Eq. 31 is obtained for the tangent elastoplastic modulus with damage.It describes the constitutive behaviour for the elastic and plastic damaged states. where . It can be noticed that Eq. 27 is recursive, and the damage evolution behaviour is nonlinear with respect to the plastic deformation.A simple relation can be stated between the rupture stress σ R and the ultimate tensile strength σ u to express the critical damage associated to rupture [6], as seen in Eq. 33.
The following assumption by Eq. 34 is made based on hardening principles, for the ultimate tensile strength.With the increase of the hardening effects for very large deformations, the tendency is to have a saturation of that hardening, with an assymptotic value for the stress.Thus, Eq. 34 just states that this limit point is reached at the ultimate tensile strength, with the values for stress very close to the effective rupture stress σ inf , as illustrated in Fig. 5 σ Figure 5. Linear damage evolution assumption with saturation on hardening.
Substituting Eq. 34 in Eqs 27 and 29, one has, Stating initial conditions with the damage threshold, which is the value of strain when the damage is initiated, Hence, the damage evolution is now dependent on a damage threshold.For the critical damage (D c ) associated to rupture (R), when p = R then D = D c on Eq. 36.
Separating some terms in Eq. 37, Equation 35 can be expressed as, Hence, Eq. 39 characterizes a linear damage evolution law [6].However, it can be noticed that the use of this equation can lead to very different results than Eq.35, specially if the damage threshold holds a small value.Some analyses concerning Eqs.35 and 39 are presented in Section 5.

Return mapping algorithm including damage
The inclusion of damage to the elastoplastic material requires now an iterative procedure, different from the undamaged case, where the solution could be obtained in a straightforward manner.The procedure is described below [5].
Let the initial conditions be λ p i+1 = λ p i , α i+1 = α i and D i+1 = D i , where i+1 representing the current equilibrium iteration, i stands for the last converged equilibrium state, λ denotes the stretch, α the isotropic hardening, and D the damage variable.Algorithm 2 in the box shows the procedure for calculating a trial state, where the dissipative variables (α, p, D) are not considered, and only an elastic prediction is made.
Algorithm 2 Equations for the elastic prediction phase.
Estimating a trial stretch, and then the appropriate tensor (logarithmic strains in this case): Calculating a trial stress.
The condition for the yield surface on Eq. 40 is then tested, to verify if the material is in the elastic or plastic domain For elastoplasticity with isotropic hardening, If the axiom in Eq. 40 is satisfied, the material is in the elastic domain, and the dissipative variables are updated with their values "frozen", and the next equilibrium step is achieved, that is, Otherwise, a correction phase is needed.

Correction phase:
Before showing the algorithmic step of this phase, the need of an iterative procedure is demonstrated.The incremental form is expressed by the following equations, with, Notice by Eqs 45 -48 that it is not possible to represent the damage variable D i+1 in an explicit form for the plastic multiplier ∆γ, so the elastoplastic procedure for solving both for the plastic multiplier and the damage variable must be solved iteratively.Following Driemeier [5], the iterative procedure is presented by Algorithm 3. Different indexes refer to the equilibrium step (i) and the current iteration for the correction phase (j).

DEVELOPED SOFTWARE
A software intended for academic support was developed based on the theory presented in the previous sections.The resulting code is a continuation of software ATENAS, which has the first version in 2001 [10].The code already had implementations for large displacements, elastoplasticity and viscoplasticity [1] for static and dynamic analyses [2].In the present version -ATENAS 5.0 -the software gives a first insight into elastoplastic damage in bars and into a simplified bending stiffness for beam elements [12].The equilibrium nonlinear equations are solved using either the Newton-Raphson Method or the Displacement Control Method.In the case of nonlinear transient problems, two methods are implemented: an implicit one, using a Newmark algorithm and an explicit one which employs central differences (average accelerations) [3,4].Inspired in the work of Driemeier et al. [5], some examples have been studied and the results are discussed in Section 6.All the transient cases in this paper were solved using the Newmark algorithm.
Algorithm 3 Iterative procedure for the correction phase.
1.At the first iteration (j = 0), assuming D j = D i : With: Then the third step is initiated: and,

Analysis of the damage evolution laws
The main goal of this analysis is to determine when it would be convenient to use Eq.39 (linear damage evolution) instead of Eq. 35 (nonlinear damage evolution), and to verify the behaviour of the nonlinear evolution law at high values for the damage threshold.This can be made by changing the parameter S in Eq. 35.All the tests in this section were performed with D c = 0.40 and R = 0.20.Figure 7 shows the comparison for a damage threshold = 0.07 and S = 0.75, 1.00, 1.25, 1.50 M P a.Even with the proximity between S = 1.50 M P a and the linear evolution law, for small values of damage, there is a disagreement for greater values, with smaller values of S achieving prior rupture.Figure 6 presents results for S = 1.50, 1.75, 2.00, and 2.25 M P a.Even with a smaller damage rate, the behaviour is still distinct.A higher damage threshold of D = 0.165 was tested, seeking more proximity between solutions, as predicted by the linear approximation from Eq. 39. Figure 8 shows the results for S = 0.75, 1.00, 1.25, and 1.50 M P a.It can be noticed that for S = 1.25 M P a, the correspondent curve almost fits the curve from the linear damage evolution.

Snap-through problem
A classical non-linear problem is shown in Fig. 9.It is characterized by two limit points, in which the force x displacement curve changes the slope signal (see Fig. 10).The material and geometrical properties are E = 210 GP a, K = 100 GP a, σ y = 150 M P a, S = 10 M P a and A = 7 mm 2 .Static and dynamic analyses were performed, using logarithmic strain measure and no damage threshold.

Static analysis
With the purpose of comparing results to those reported by [5], this analysis was performed using a Newton Raphson procedure with prescribed displacements applied on node 2, from 0 mm to -300 mm.
Figure 10 shows the force x displacement relation for elasticity, elastoplasticity and elastoplasticity coupled to damage.The damage consideration led to failure briefly after the limit point in traction.This failure occurred with a displacement of approximately -231 mm and an accumulated plastic strain of ≈ 0.105.Although the results are quite similar to those obtained by the reference work [5], the latter reports different values for S, K and E (6 M P a, 10 GP a and 21 GP a, respectively).
It should be remarked that a more realistic formulation would consider different rates of damage evolution for compressive and tensile stress states.Compressive stress states lead to a much lower damage evolution rate and hence, the truss would not fail until attaining a larger final displacement.

Dynamic analysis
The following parameters were used for the dynamic analysis: E = 210 GP a, K = 100 GP a, σ y = 700 M P a, A = 7 mm 2 , ρ = 7.850 × 10 −6 kg mm 3 .A load of 1100 kN was applied in a time interval of 0.250 s, using 3000 increments (hence ∆t = 8.33 × 10 −5 s for each).
The first test was carried out without damping effects.Figure 11 shows the results for the displacement.The snap-through phenomenon occurs at t ≈ 0.22 s for elasticity and at t ≈ 0.09 s for isotropic elastoplasticity coupled to damage.The latter consideration led to 8 cycles before achieving the critical damage (D = 1).When compared to the case of isotropic elastoplasticity, it is possible to notice a reduction of the system natural frequency, given the increasing relative distance between the displacement peaks after the occurrence of the snapthrough phenomenon.Figure 12 presents the data obtained for the stress x total deformation relationship.In each loading cycle, traction follows compression, both for elastoplasticity and for elastoplasticity coupled to damage.Also, in each loading step, the elastoplastic material suffers additional yielding due to the cyclic pattern shown in Fig. 11.The dashed curve, which includes damage effects, shows considerable and increasing reduction of the elasticity modulus.
The second test was performed including proportional mass damping, arbitrarily setting C t = 10M, where C t and M are the damping matrix and the lumped mass matrix respectively.Figure 13 shows the displacement over time.The damping effect supressed considerably the stresses caused by inertial effects, and did not lead to material failure after the load application.Figure 14 shows a stress behaviour similar to the undamped analysis.The material was still severely damaged, although with a slower damage evolution, due to lower yielding.

Star truss
Another problem which is also presented by Driemeier et al. [5] is the star truss shown in Fig. 15.Node 1 has a z-coordinate of 72.16 mm, white nodes have z = 42.16mm, and the dark nodes have a z-coordinate of 0.00 mm.Material and geometrical properties are E = 200 GP a, K = 200 M P a, σ y = 250 M P a, S = 0.003 M P a, A = 20 mm 2 .There are three geometrically distinct elements, namely 1, 2 and 3.As in the 2-bar snap through analysis, this test used logarithmic deformations and zero damage threshold.

Static analysis
In this case, node 1 in Fig. 15 is displaced from 0 mm to 100 mm in the negative direction of z. Figure 16 shows the results.The Newton-Raphson procedure could not achieve convergence beyond u = −25 mm for elastoplasticity due to numerical difficulties.A bifurcation phenomenon occurs close to u = −100 mm, with the sudden reduction of forces.For elastoplasticity coupled to damage, the results are close to those obtained by Driemeier et.al [5] until u ≈ −60 mm.Beyond that point, results for elastoplasticity coupled to damage are different, possibly because of a bifurcation.

Dynamic Analysis
For this analysis, a load of −40 kN in the z direction is applied on node 1 during 0.3 s.Then, the load is removed, and the behavior of the structure is analyzed in the following 0.3 s.Thus, the total time is 0.6 s with 2000 increments, resulting in ∆t = 3, 00 × 10 −4 s.The properties are E = 200 GP a, K = 100 GP a, σ y = 250 M P S = 0.1 M P a and A = 20 mm 2 .
Figure 17 shows the displacement of node 1 over time, without damping.In comparison to the 2-bar snap-through problem (Fig. 11), the displacement pattern is more complex due to the geometry and to interaction between elements of the star truss.The snap-through phenomenon (node 1 with the first intersection plane) occurs around t = 0.075 s.Larger displacements are evident when elastoplasticity is considered, and the damage variable causes loss of stiffness, achieving critical damage and leading to failure at t = 0.18 s.At failure, critical damage occurs in elements 3, as shown in Fig. 18.The stress-strain relationship is presented in Fig. 19.Different from the elastoplasticity case, the inclusion of damage did not result in compression cycles after the snap-through occurence.Once more, it must be noted that the present formulation overestimates the damage effect because the same treatment is given to tractive and compressive stress states.The second dynamic analysis for the star truss uses the same material parameters as the previous one, but considers a proportional mass damping C t = 10M, where M is the lumped mass matrix.The number of time increments is arbitrarily set to 4000 during 0.6 s of analysis, and hence ∆t = 1.50 × 10 −4 s. Figure 20 shows the displacement over time.The effects of damping are clear, decreasing the oscillations and taking more time to reach the critical damage (t = 0.23 s) for the damaged element.The results for damage evolution over time and stress-strain are quite similar to those of the undamped test.

CONCLUSIONS
Regarding the damage evolution law comparison, there is a proximity between both linear and nonlinear laws for sufficiently high values of damage threshold, and specific values of the S parameter.A better prediction of a convenient damage threshold could be performed with the implementation of a nonlinear hardening law, which would identify an assymptotic line correspondent to a saturated hardening.
For the structural examples, it is possible to notice the severe degradation of the mechanical properties when damage effects are considered.Also, the damage parameter S plays a crucial role on the evolution, along with the critical damage and damage threshold.Damage effects could be identified on the analyzed structures, by the reduction of the natural frequency as shown in displacement over time plots.The dynamic analysis showed the influence of inertia forces and damping on damage evolution.In particular, inertia forces cannot be neglected in unstable problems, since accelerations are typically high and stresses are affected drastically.Therefore, yielding and damage are much higher than in the quasi-static counterpart.
The formulation presented in this work constitutes an interesting academic development, giving insight to the study of damage inclusion in transient elastoplastic problems.However, unrealistic results are still obtained because damage evolution in traction and compression is treated in the same manner.For better results and further analyses, it is important to introduce variables h+ and h-to account for the difference of damage evolution in tensile (crack opening) and compressive (crack closure) stress states, where the latter causes less damage to the material [14].

Figure 4 .
Figure 4. Bar subjected by a traction force, after the removal of cracks and voids.

Figure 10 .
Figure 10.Force x displacement for the static analysis.

Figure 11 .
Figure 11.Displacement x time for the dynamic analysis, no damping.

Figure 12 .
Figure 12.Stress x total deformation for the dynamic analysis, no damping.

Figure 13 .Figure 14 .
Figure 13.Displacement x time for the dynamic analysis, with damping.