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 rulewhere 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.
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)
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)
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.