Advances in Materials Science and Engineering
Volume 2016 (2016), Article ID 5289642, 9 pages
http://dx.doi.org/10.1155/2016/5289642

Research Article

Semi-Implicit Algorithm for Elastoplastic Damage Models Involving Energy Integration

Ji Zhang1,2 and Jie Li²

¹School of Civil Engineering and Mechanics, Huazhong University of Science and Technology, Wuhan 430074, China
²School of Civil Engineering, Tongji University, Shanghai 200092, China

Received 11 September 2015; Accepted 28 February 2016

Academic Editor: Claudio Borri

Copyright © 2016 Ji Zhang and Jie Li. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

This study aims to develop a semi-implicit constitutive integration algorithm for a class of elastoplastic damage models where calculation of damage energy release rates involves integration of free energy. The constitutive equations with energy integration are split into the elastic predictor, plastic corrector, and damage corrector. The plastic corrector is solved with an improved format of the semi-implicit spectral return mapping, which is characterized by constant flow direction and plastic moduli calculated at initial yield, enforcement of consistency at the end, and coordinate-independent formulation with an orthogonally similar stress tensor. The tangent stiffness consistent with the updating algorithm is derived. The algorithm is implemented with a recently proposed elastoplastic damage model for concrete, and several typical mechanical tests of reinforced concrete components are simulated. The present semi-implicit algorithm proves to achieve a balance between accuracy, stability, and efficiency compared with the implicit and explicit algorithms and calculate free energy accurately with small time steps.

1. Introduction

Due to its accuracy, fast convergence, and unconditional stability, the (fully) implicit backward Euler return mapping method [1] has become the most favorable one for numerically integrating inelastic constitutive models. In this method, quite large time steps can be used, hence far fewer global (finite element) iterations.

However, sufficiently small time steps are warranted to accurately calculate free energy for a class of elastoplastic damage models. In these constitutive theories, a damage energy release rate is derived from a Helmholtz free energy potential. Unlike the elastic free energy (i.e., strain energy), which is determined from the current state, the plastic free energy (a part of the plastic work [2]) has to be integrated along the plastic loading history [3, 4].

With very small time steps and frequent iterations, the implicit return mapping method is too expensive computationally, especially for quasi-brittle material models, which usually employ complex hardening laws and flow rules. The extremely complex partial derivatives of the plastic moduli and flow direction with respect to the stress and hardening variables are required for the process of return mapping, and the inversion of certain matrices containing these derivatives needs to be carried out to construct the consistent tangent stiffness.

Return mapping can be made much simpler with the so-called semi-implicit method [5], where some quantities, including the plastic moduli and flow direction, are calculated at the beginning of a time step and remained unchanged throughout the process of return mapping, and some, including the plasticity consistency parameter, are determined at the end of the time step. For a single small time step, this approach is much easier to implement and more efficient than the implicit method; at the same time, the stress is return-mapped effectively onto the subsequent yield surface by enforcing the consistency condition, and the problem of drift in the explicit forward Euler method is overcome. Thus, the semi-implicit return mapping method is an ideal choice for the abovementioned plastic free energy integration.

It is the purpose of the present study to develop a semi-implicit constitutive integration algorithm for a class of elastoplastic damage models where calculation of damage energy release rates involves integration of free energy. It is intended that, with small time steps, free energy can be integrated accurately, and the numerical algorithm remains efficient and stable.

2. Constitutive Framework

The class of elastoplastic damage models under discussion are mainly for quasi-brittle materials, which exhibit distinct behavior in tension and compression. The general constitutive framework is presented as follows.

2.1. Overall State Description

Two scalars and , the tensile and compressive damage variables, are used to describe a damage state as

where is the nominal (Cauchy) stress and is the tensile/compressive part of the effective stress :

in which are the principal effective stresses, is the orthonormal basis along the principal directions, and .

Elastoplasticity is formulated in the effective stress space:

where is the effective elastic stiffness, is the strain, and is the plastic strain.

2.2. Damage Evolution: Drive Chain

The tensile/compressive Helmholtz free energy is split into its elastic part and plastic part , mapped into their effective values and , and calculated as

where is the elastic strain, is the tensile/compressive effective plastic energy storage rate [4], and is the tensile/compressive hardening variable.

The tensile/compressive damage energy release rate is defined and calculated as

The tensile/compressive damage threshold is defined as

where is the initial tensile/compressive damage threshold, is any instant throughout the loading history, and is the current instant.

Damage evolution is formulated as a function of and :

2.3. Plastic Evolution: Intertwined Components

Two scalars and , the tensile and compressive hardening variables, are used to characterize plasticity, which includes the yield criterion,

and the flow rule

where is the plasticity consistency parameter and is the plastic potential, and the evolution law for is as follows:

where is the tensile/compressive plastic modulus.

Finally, with the Karush–Kuhn-Tucker loading/unloading conditions,

and the consistency condition,

the plastic relation is complete.

3. Operator Split

In accordance with the concept and methods of product formulas and operator split [3, 6], the constitutive equations presented in the previous section are replaced here with the elastic, plastic, and damage parts in Table 1, which can serve as a mathematical foundation for developing first-order accurate algorithms.

Table 1: Split of rate constitutive equations.

Total
=Elastic part
+Plastic part
+Damage part

These three parts are taken as split problems in a series: the first problem provides the initial condition for the second one, the second one for the third one, and the third one for the first one in the next time step; and all the split problems add up to the original elastoplastic damage problem. In this manner, the elastoplastic damage constitutive equations can be integrated incrementally over a sequence of time steps .

3.1. Elastic Predictor

The first problem, usually referred to as the elastic predictor, is simply the update of the strain

and the effective stress

assuming purely elastic behavior.

3.2. Plastic Corrector

The second problem, the plastic corrector, is to relax the plastic strain increment and direct the excursive effective stress back to the updated yield surface:

if violates the yield criterion

Various types of algorithms have been developed for this problem, including the explicit forward Euler method, implicit backward Euler method, and semi-implicit method, which differ from one another in their specific (iterative) updating schemes for (15). In a broader sense, all these methods can be regarded as “return mapping”; many researchers mean only those that accurately enforce consistency by this popular term.

3.3. Damage Corrector

The third problem, the damage corrector, is to evaluate damage evolution:

Like the elastic predictor, it is also straightforward and involves no iteration. However, its formulation plays an important role in the construction of the algorithmic consistent tangent stiffness as demonstrated in the following Section 5.

4. Semi-Implicit Return Mapping

For the plastic corrector problem, the semi-implicit method [5] is adopted and improved in the present study in order to integrate free energy accurately and efficiently with small time steps .

4.1. General Format

Different from the implicit backward Euler method, the flow direction and plastic moduli at the beginning of the elastoplastic loading range are used in the whole elastoplastic range , where is the relative length of the purely elastic range , so that the iterations become much easier and simpler [7]. Different from the explicit forward Euler method, the consistency condition is enforced at the end of the step so that drift no longer occurs.

The set of nonlinear equations (8)–(10) are cast in a semi-implicit format of their residuals:

where

To solve (18), the Newton-Raphson iterations

are used, where the superscript () denotes the th iteration within a time step, leading to

where , , and are obtained by solving (20).

Remark 1. It is noted that when , that is, the stress state is on the current yield surface at the beginning of the time step , the above procedure reduces to the existing semi-implicit format [5].

4.2. Spectral Version

The plastic constitutive equations for quasi-brittle materials often involve terms of principal stresses. As a consequence, eigenvalue calculation needs to be performed for each iteration () in the return mapping process developed in the previous subsection, which substantially increases the computational cost. For these cases, the spectral version [8, 9] of the return mapping methods is much more efficient. It is adopted and formulated rigorously in the present subsection.

The matrix corresponding to the tensor in the structural coordinate system is spectrally decomposed as

where and , which corresponds to in the principal coordinate system.

In this study, a new tensor corresponding to the diagonal matrix in the structural coordinate system is defined, and (22) is expressed in a tensor form:

where is the tensor corresponding to in the structural coordinate system. Because is an orthogonal matrix, is an orthogonal tensor. By definition, the newly introduced tensor and the traditional tensor are orthogonally similar.

It has been proved [8, 9] that, for an isotropic material with coaxial plastic flow, that is, a plastic potential exists such that (9) holds, after return mapping in (15) has the identical eigenspace with . Thus, corresponding to in the structural coordinate system can be spectrally decomposed with the same eigenvector matrix as

where , which corresponds to in the principal coordinate system.

In a similar way, corresponding to in the structural coordinate system is defined, and (24) is expressed in a tensor form:

which indicates that and are orthogonally similar.

The return mapping procedure expressed by (18)–(21) is carried out here on the orthogonally similar tensor as

where eigenvalue calculation needs to be performed only for the first iteration (), and then are iterated upon with fixed. The existence of the functions , and in (26) and (27) is based on the assumption that the material is isotropic [10].

Remark 2. It can be seen from (28) and (29) that the introduction of the orthogonally similar tensor is necessary for the coordinate-independent implementation of the double dot operation between higher-order tensors.

5. Consistent Tangent Stiffness

The differentiation with respect to the strain increment of the resulting and gives, respectively, the tangent stiffness consistent with the updating algorithms in elastic-plastic parts and elastic-plastic-damage parts:

which lead to

On the basis of the continuum elastoplastic tangent stiffness expression

replacing the flow direction and plastic moduli with their values at and the other quantities at , one obtains

in (31) needs to be derived in accordance with specific damage formulation. From (1), one can obtain

where

as derived in [11], and

In (36), depend on damage evolution functions and on definitions of . If the formulation of in [4] is adopted, the specific expressions read

where

6. Implementation

The complete algorithm of the present study is summarized as follows.

Summary of integration algorithm is as follows:

As numerical examples, this algorithm is implemented with the specific elastoplastic damage model for concrete presented in [4]. Several typical mechanical tests of reinforced concrete components are simulated with the commercial finite element program Abaqus and its user material subroutine UMAT. Solid elements are used for concrete and truss elements for steel rebars, between which perfect bonding is assumed.

The uniaxial mechanical properties for concrete are first calibrated to fit the stress-strain curves of uniaxial tension and compression recommended in the Chinese code [12]. The tension stiffening effects due to the presence of reinforcement are accounted for by using the descending stress-strain curve by Stevens et al. [13] and the confining effects resulting from closely spaced stirrups by fitting the stress-strain curve for core concrete by Kent and Park [14] and Scott et al. [15]. Besides, the common values for the equibiaxial strength ratio , Poisson’s ratio , and dilation parameter are set for multiaxial behavior. For reinforcing steel, uniaxial elastoplastic models are employed.

6.1. Slab

For application in reinforced concrete slabs, the test by McNeice [16] (Figure 1(a)) is simulated. The two-way slab has a bottom reinforcement ratio of 0.85% in each direction. In the simulation, Young’s modulus and yield point are set for steel, and effective Young’s modulus , initial damage stresses in uniaxial tension and compression , , and initial effective yield stresses in uniaxial tension and compression , for concrete. The calculated load-deflection curve is shown along with the test data in Figure 1(b).

(a)

(b)

Figure 1: Reinforced concrete slab under out-of-plane loading [16]. (a) Setup (dimensions in mm). (b) Results.

For comparison, the implicit backward Euler return mapping algorithm is also implemented, and the calculated load-deflection curve is precisely the same as in Figure 1(b). The times of calculation consumed with the semi-implicit and implicit methods are 179 s and 804 s, respectively, at 2000 displacement steps.

6.2. Wall

For application in reinforced concrete walls, the tests by Lefas et al. [17] (Figure 2(a)) are simulated. In one of the two selected tests (SW22 and SW24), a constant vertical load is held while the horizontal load increases. These loads are transferred through the top spreader beams, which are cast integrally with the walls. The vertical and horizontal reinforcement ratios in the walls are 2.49% and 0.82%, respectively, and the ratio of closed stirrups to confined concrete of 140 mm-wide concealed columns at both edges is 0.52%. In the simulations, , (vertical rebars), (horizontal rebars), and (stirrups) are set for steel, , , , , and for concrete of SW24, and , , , , and for concrete of SW22. The calculated load-deflection curves are displayed along with the test data in Figures 2(b) and 2(c), respectively, where is the cubic strength of concrete.

(a) 

(b)

(c)

Figure 2: Reinforced concrete walls under in-plane loading [17]. (a) Setup (dimensions in mm). (b) Results for SW24 (, N/mm²). (c) Results for SW22 ( kN,  N/mm²).

For comparison, the implicit return mapping algorithm is also implemented, and the calculated load-deflection curves are precisely the same as in Figures 2(b) and 2(c). The times of calculation consumed with the semi-implicit and implicit methods are 495 s/508 s and 1870 s/1922 s, respectively, at 2000 displacement steps.

7. Conclusions

History-dependent energy calculation makes the quality of stability in large time steps of the popular implicit backward Euler return mapping algorithm no longer attractive for some inelastic constitutive models. More flexible and efficient integrating schemes are favored.

The semi-implicit constitutive integration algorithm developed for a class of elastoplastic damage models that involve integration of free energy in the present study achieves a balance between accuracy, stability, and efficiency. In comparison with the implicit backward Euler algorithm, it is by far less challenging mathematically and less time-consuming computationally. While it is less stable, the problem can be fixed with small time steps, which are necessary for free energy integration. Another drawback is the need to find the stress where the yield surface is first reached within a time step. In comparison with the explicit forward Euler algorithm, it is much more accurate and stable at the cost of iterations for consistency.

The splitting of operators and the elastic-plastic-damage updating algorithm on its basis presented in this paper are particularly formulated for damage energy release rates that are derived from elastoplastic free energy potentials. Nevertheless, the basic idea and method are also applicable to other types of constitutive models that include history-dependent calculation of certain quantities.

Competing Interests

The authors declare that they have no competing interests.

Acknowledgments

The financial support from the National Science Foundation of China-National Science Foundation (NSFC–NSF) Joint Project (Grant no. 51261120374) and NSFC Projects (Grant nos. 51108336 and 51378377) is greatly appreciated.

References

J. C. Simo and T. J. R. Hughes, Computational Inelasticity, Springer, New York, NY, USA, 1998.

J. Zhang and J. Li, “Microelement formulation of free energy for quasi-brittle materials,” Journal of Engineering Mechanics, vol. 140, no. 8, Article ID 06014008, 2014.

J. W. Ju, “On energy-based coupled elastoplastic damage theories: constitutive modeling and computational aspects,” International Journal of Solids and Structures, vol. 25, no. 7, pp. 803–833, 1989.

J. Zhang and J. Li, “Elastoplastic damage model for concrete based on consistent free energy potential,” Science China Technological Sciences, vol. 57, no. 11, pp. 2278–2286, 2014.

F. Dunne and N. Petrinic, Introduction to Computational Plasticity, Oxford University Press, New York, NY, USA, 2005.

A. J. Chorin, T. J. R. Hughes, M. F. McCracken, and J. E. Marsden, “Product formulas and numerical algorithms,” Communications on Pure and Applied Mathematics, vol. 31, no. 2, pp. 205–256, 1978.

J. Zhang and J. Li, “Construction of homogeneous loading functions for elastoplastic damage models for concrete,” Science China: Physics, Mechanics and Astronomy, vol. 57, no. 3, pp. 490–500, 2014.

J. C. Simo, “Algorithms for static and dynamic multiplicative plasticity that preserve the classical return mapping schemes of the infinitesimal theory,” Computer Methods in Applied Mechanics and Engineering, vol. 99, no. 1, pp. 61–112, 1992.

J. Lee and G. L. Fenves, “Return-mapping algorithm for plastic-damage models: 3-D and plane stress formulation,” International Journal for Numerical Methods in Engineering, vol. 50, no. 2, pp. 487–506, 2001.

I. H. Shames and F. A. Cozzarelli, Elastic and Inelastic Stress Analysis, CRC Press, Boca Raton, Fla, USA, 1997.

R. Faria, J. Oliver, and M. Cervera, On Isotropic Scalar Damage Models for the Numerical Analysis of Concrete Structures, CIMNE Monograph no. 198, Centro Internacional de Métodos Numéricos en Ingeniería, Barcelona, Spain, 2000.

Chinese National Standard, “Code for design of concrete structures,” Tech. Rep. GB50010, China Architecture & Building Press, Beijing, China, 2002.

N. J. Stevens, S. M. Uzumeri, M. P. Collins, and G. T. Will, “Constitutive model for reinforced concrete finite element analysis,” ACI Structural Journal, vol. 88, no. 1, pp. 49–59, 1991.

C. D. Kent and R. Park, “Flexural members with confined concrete,” Journal of Structural Engineering, vol. 97, no. 7, pp. 1969–1990, 1971.

B. D. Scott, R. Park, and M. J. N. Priestley, “Stress-strain behavior of concrete confined by overlapping hoops at low and high strain rates,” Journal of the American Concrete Institute, vol. 79, no. 1, pp. 13–27, 1982.

G. M. McNeice, Elastic-Plastic Bending Analysis of Plates and Slabs by the Finite Element Method, University of London, London, UK, 1967.

I. D. Lefas, M. D. Kotsovos, and N. N. Ambraseys, “Behavior of reinforced concrete structural walls: strength, deformation characteristics, and failure mechanism,” ACI Structural Journal, vol. 87, no. 1, pp. 23–31, 1990.