2014, Vol. 35
The Chinese Meteorological Society
Article Information
- Chao-jun LI, Ji-li FENG 2014.
- Algorithmic tangent modulus at finite strains based on multiplicative decomposition
- Appl. Math. Mech. -Engl. Ed., 35 (3) : 345–358
- http: //dx. doi. org/10.1007/s10483-014-1795-6
Article History
- Received 2013-01-09;
- in final form 2013-05-31
2 School of Mechanics and Civil Engineering, China University of Mining and Technology (Beijing), Beijing 100083, P. R. China
1 Introduction
The initial developments of elastoplastic constitutive relation at finite strains depend on the use of hypoelastic-based constitutive formulations which are characterized in terms of suitably chosen objective stress rates and extended from standard infinitesimal elastoplasticity models[1, 2]. The hypoelastic-based descriptions have two controversial issues which include fundamental drawbacks such as the possible lack of the objectivity of (algorithmic) incremental constitutive laws[3] and dissipative behaviors within the “elastic” range[4]. In the abovementioned context,the hyperelastic-based formulations of finite plasticity by the local multiplicative decomposition have emerged[5],which naturally bypass the inherent drawbacks of hypoelastic-based approaches. As a consequence,the dissipative response becomes impossible within the elastic range. Another merit is that the requirement of the incremental objectivity (the frame invariance of the algorithmic constitutive rule) is trivially satisfied. Generally speaking,the closest point projection algorithm is applied for the analysis of the elastoplastic constitutive equation. The algorithm is based on the elastic predictor-plastic corrector method which deals with two steps. The first step is to compute the trial stress with no plastic strain, i.e.,all strains are considered as elastic strains,and the second step is to analyze the yield function value to check the yield condition according to the trial stress. If the elastic state is violated,the trial stress must be returned to the yield surface according to the discrete Kuhn-Tucker condition. Otherwise,the trial stress is the real stress and exits. In the nonlinear finite element method,the algorithmic modulus must be computed and the chosen algorithmic tangent modulus can gain the quadratic convergence rate when using the Newton-Raphson method. Consequently,the analysis of the algorithmic tangent modulus in the current configuration plays an important role in the nonlinear finite element method. Though the algorithmic tangent modulus is well understood under small deformations in the literature,there are some ambiguities when considering finite strains. The algorithmic tangent modulus in the principal space is one valid method to investigate the finite deformation of solid. The derivations of the algorithmic tangent modulus at finite strains in the principal space have been presented for the case λ1 ≠ λ2 ≠ λ3 and the concise forms have been obtained[6]. However,the matrix expressions what we expect in computational mechanics are not obtained. In particular,the other two cases that include λ1 = λ2 ≠ λ3 and λ1 = λ2 = λ3 are not clear,and their matrix expressions are also not provided. Moreover,the valid matrix expressions of the fourth-order tensor are controversial issues[7].
This paper is organized as follows. In Section 2,the exact tensorial forms of the algorithmic tangent modulus at finite strains are derived,and their corresponding matrix expressions
are given in Section 3. The elastoplastic matrix in the principal space associated with the
specific yield surface is derived by the logarithmic strain in terms of the local multiplicative
decomposition in Section 4. The Drucker-Prager yield function of elastoplastic material is used
as a numerical example to verify the present algorithmic tangent modulus at finite strains in
Section 5. In the present work,the number of the underlines of a notation corresponds to its
order,e.g.,
represents a first-order tensor,
is a second-order tensor,and 
indicates a
fourth-order tensor; and the symbol “:” in equations denotes as a double contract operator of
two tensors,e.g.,
means a scalar,and 
represents a second-order tensor.
The local multiplicative decomposition of the deformation gradient 
is
where 
and 
are the elastic and plastic parts of the deformation gradient,respectively[11].
The right and left Cauchy-Green deformation tensors are defined,respectively,as
The logarithmic strain tensor is defined as
In the principal space,the Kirchhoff stress tensor 
and the left Cauchy-Green deformation
tensor 
are coaxial because of the restriction of isotropy,and their respective spectral decompositions take the forms[8]
where λA2 is the eigenvalue of the left Cauchy-Green deformation tensor,βA indicates the
eigenvalue of the Kirchhoff stress,and 
(A = 1,2,3) represent the eigenvectors in the current
configuration.
The eigenvalues of the right Cauchy-Green deformation tensors are λA2 (A = 1,2,3),
while their corresponding eigenvectors are 
in the reference configuration,expressed as
. Consequently,its spectral decomposition is
Let us define
For the case λ1 ≠ λ2 ≠ λ3,
can also be denoted as[6]
where 
,and 
is the second-order identity tensor.
In the reference configuration,the algorithmic tangent modulus is defined as
where 
is the second Piola-Kirchhoff stress tensor.
From Eq. (8),the algorithmic tangent modulus expressed in the current configuration can be obtained as follows:
The following three basic relations are given:
which are proved in Appendix A.
The exact derivations of the algorithmic tangent modulus in the current configuration are given for the three different cases as follows.
Case I λ1 ≠ λ2 ≠ λ3
Since
where I1 and I3 are the first and third invariants of ,respectively.
Based on Eqs. (14) and (15),Eq. (13) is exactly expressed as
where
According to Eq. (9),we finally obtain
where aep derived in Section 4 is the elastoplastic modulus associated with the special yield
surface,and 
is independent of the specific model of plasticity since the principal directions
are functions of the total deformations expressed by
Case II λ = λ1 = λ2 ≠ λ3
After some operations,the exact expression of Eq. (18) is
The algorithmic tangent modulus in the current configuration reads
where
and the fourth-order tensors of the last two terms in the parentheses are independent of the specific yield surface.
Case III λ = λ1 = λ2 = λ3
The algorithmic tangent modulus in the current configuration is
where 
,and 
is independent of the specific yield surface.
3 Matrix expressions of algorithmic tangent modulus in current configuration
As done in Ref. [7],let the second-order base-tensors be
Therefore,the following matrix expressions are obtained:
From the definition of Eq. (21),the matrix expressions of Eqs. (17),(19),and (20) are summarized in Box 1 after some operations.
Box 1 Matrix expressions of algorithmic tangent modulus
Case I λ1 ≠ λ2 ≠ λ3
Case II λ = λ1 = λ2 ≠ λ3
where
Case III λ = λ1 = λ2 = λ3
4 Modulus aep associated with specific yield surface
For a mixed hardening model,the general yield surface can be written as
where 
and H are the Cauchy stress tensor,the determinant of the deformation
gradient,the back stress tensor associated with the kinematic hardening,and the isotropic
hardening variable,respectively.
From the consistent condition,we obtain 
= 0,i.e.,
The incremental Cauchy stress tensor can be denoted by the incremental logarithmic strains
where 
is the modified elastic modulus in the principal space.
If the logarithmic strain is employed,the additive decomposition of the strain is expressed as
Proof Considering the isotropic hypothesis,from Eqs. (1) and (3),we have
After expanding Eq. (26) and recombining it,we obtain
where 
,and 
.
In the principal space,the principal directions of the logarithmic strain identify with the left Cauchy-Green tensor,and the corresponding respective spectral decomposition is
where εA = lnλA.
Due to aABep =(∂βA/∂λB)λB,we can rewrite aABep in the form
At finite strains,the non-associated flow rule can be,from the return mapping method,approximately expressed as
where gp is the plastic potential function,which is identified by the yield surface function under the consideration of the associated flow rule.
Combining Eqs. (23),(24),(25),and (30),after some operations,we obtain
where 
is associated with the hardening modulus.
Substituting Eqs. (25) and (31) into Eq. (24),we finally obtain the exact form of the stressstrain formulation
The corresponding matrix expression of the elastoplastic tangent modulus 
is
In Eq. (32),the term 
is written as[9]
where
5 Numerical example
Large deformation of geomaterial is one of the mostly concerned issues in the mechanical research field[6, 12, 13, 14]. In order to examine the present algorithmic scheme,the Drucker-Prager model of elasto-plasticity most commonly used in the geotechnical engineering analysis with the associated flow rule is employed in the numerical example. The model formulation with respect to three principal stresses is first presented. The algorithmic tangent modulus is then coded by MATLAB. As a comparison,the elastoplastic tangent modulus at infinitesimal strains is also presented.
The Drucker-Prager yield surface or function reads
where 
is the second invariant of the deviatoric
stress tensor,I1 = σ1 + σ2 + σ3 is the first invariant of the stress tensor,and α and k are
material parameters,respectively.
The first-order derivatives of the yield function with respect to the three principal stresses are
Under the conditions of the no hardening and associated flow rule,Eq. (33) is rewritten as
The corresponding elastoplastic tangent modulus at infinitesimal strains is simplified as follows:
As a typical example of clay,Poisson’s ratio and Young’s modulus are µ = 0.28 and E = 10.3 MPa,respectively. The elastic modulus is expressed as
Furthermore,c = 0.11 MPa and φ = 0.14π are used. It is provided that the yield surface of the Drucker-Prager and Mohr-Coulomb criteria coincide at the outer edges of the Mohr-Coulomb surface. Consequently,we can obtain
The following procedures with one main program are developed by MATLAB to compute the elastoplastic tangent modulus of the Drucker-Prager model. It should be noted that,in Boxes 2 and 3,[·]n+1k and (·)n+1k are a matrix [·] and a variable at the kth iteration of the load step during the time [tn,tn+1],respectively.
i) Compute Eq. (34) and return a 1 × 3 matrix,and its implementation is shown in Box 2.
ii) Compute the elastoplastic tangent modulus associated with the Drucker-Prager yield surface and update the stresses and strains by means of the closest point projection method[10], which is summarized in Box 3.
iii) Compute the algorithmic tangent modulus.
Box 2 Compute the derivative of J with respect to the logarithmic strains ε in the principal space
(i) Give 
. Set the initial values and the numerical tolerance ρtol: [Dea]0= 0,n0= 1,and
i := 1.
(ii) Set [a]i= 0. Then,compute ni= ni−1· i.
(iii) For j = 1 : i,compute 
.
(iv) Compute 
and Ttol = det[a]i,where Dea is a matrix corresponding
to the tensor form 
mentioned in Eq. (35),and Ttol is a summary of residuals at every
computational step for [Dea].
Check the condition. If Ttol < ρtol,exit; else,set i := i + 1 and go to Step (ii).
Box 3 Closest point projection method with no hardening and associated flow rule
(i) k := 0. Assume that the total strains are elastic strains. Set the initial value of the
Lagrange multiplier as 
.
(ii) Compute the values of the yield function,the trial stress,and the residual as follows:
Check the yield condition. If fn+1k < ρtol,then set (*)n+1 = (*)n+1k and exit; else,do the following calculations.
(iii) Compute the matrix [B],
(iv) Compute the incremental value of the Lagrange multiplier,
(v) Compute the incremental value of the Cauchy stress and update it,
Set k := k + 1,and go to Step (ii).
The present example is concerned with the deformation that the changes of the three principal values of the Cauchy-Green deformation are first calculated by the elastic constitutive
equation when the compressive stress σ1 increases and σ2 = σ3 = 0. The yield condition is
then checked by the trial elastic stress. If the value of the yield function is less than zero,the
corresponding stresses are at the elastic state; else,compute the elastoplastic tangent modulus
associated with the Drucker-Prager yield surface and update the stresses and strains by means
of the closest point projection method[10] and obtain the real stresses. The material parameters
of the model are summarized in Table 1. The three trial principal values of the left CauchyGreen deformation tensor are calculated for 50 deformation incremental steps from the elastic
constitutive equation,as shown in Table 2.
The typical computational results of the algorithmic tangent modulus at finite strains are given in Table 3,where aep(i,j) is the element of the ith row and jth column of the matrix [aep] expressed in Eq. (38). Before Step 47 of the Cauchy-Green deformation incremental,the stress state is elastic. However,the algorithmic tangent modulus decreases slightly due to not only the change of the principal directions which are functions of the total deformation tensor but also the change of the elastic modulus shown in Eq. (38) at finite strains,excluding the terms aep(2,3) and aep(3,2) which slightly increase. The stress state is plastic after Step 47. There are significant decreases in the diagonal terms when the stress state is changed from elasticity into plasticity between Step 46 and Step 47,and then such diagonals slightly stably decrease. However,the off-diagonal terms slightly increase,and then slowly decrease,excluding the terms aep(2,3) and aep(3,2). In the whole stage,the algorithmic tangent modulus is not symmetric, which coincides with the expression shown in Box 1.
Similarly,the corresponding typical computational results of the elastoplastic tangent modulus at finite strains are given in Table 4,where aep(i,j) is the element of the ith row and
the j column of the matrix [aep] expressed in Eq. (38). In elasticity,the elastoplastic tangent
modulus decreases slightly due to the term 
,which provokes to the reduction of the
algorithmic tangent modulus. When the stress state is changed from elasticity to plasticity,the diagonal terms first significantly decrease,and then decrease in a stabilized manner. In contrast,the off-diagonal terms first slightly increase,and then slightly decrease,excluding
aep(2,3) and aep(3,2). In the whole deformation process,the elastoplastic tangent modulus is
asymmetrical due to the term .
Table 5 lists the typical computational results of the elastoplastic tangent modulus at infinitesimal strains,where cep(i,j) is the element of the ith row and the jth column of the matrix [cep] expressed in Eq. (39). Before the Cauchy-Green deformation incremental Step 46, the stress state of material is in elasticity. The values of the elastoplastic tangent modulus are precisely the same as those of the elastic modulus,which coincide with Eq. (39). The stress state of material is in plasticity after Step 46. There are also significant decreases in the diagonal terms when the stress state of material is changed from elasticity into plasticity,and then keeps gradually reduction excluding term cep(1,1). However,the off-diagonal terms also first increase,and then slightly decrease,excluding the terms cep(2,3) and cep(3,2) as well. In the whole process,the elastoplastic tangent modulus is symmetrical.
Comparing with the computational results at finite strains in Table 4 with those at infinitesimal strains in Table 5,we can find that they are approximately identified. It is also interestingly shown from the numerical results that the tendency of the general change of the elastoplastic tangent modulus at finite and infinitesimal strains is similar during the deformation process. 6 Conclusions
The exact tensorial forms of the algorithmic tangent modulus at finite strains are derived in the pricipal space for three different eigenvalues. Their corresponding specific matrix expressions are also obtained. The algorithmic tangent modulus consists of two terms,one is independent of the specific model of plasticity,and the other associates with the chosen yield surface. The matrix expressions are detailed by establishing proper second-order tensors,as shown in Box 1. The explicit expression of the elastoplastic tangent modulus at the logarithmic strain by the local multiplicative decomposition is obtained by considering a general yield surface and a non-associated flow rule. The present algorithmic tangent modulus at finite strains is verified by a numerical example,in which the Drucker-Prager model of elastoplastic material and the associated flow rule most commonly used in engineering analysis under the conditions of no hardening are employed. Appendix A
(i) λ1≠λ2≠λ3
Differentiating the spectral decomposition of , we obtain
Contracting the relation with 
and using 
= 1 and d = 0 for A = 1, 2, 3,
we obtain
Similarly,
(ii) λ = λ1 = λ2≠λ3
(ii) λ = λ1 = λ2≠λ3 Similar to (i), we obtain
Since
we obtain
(iii) λ = λ1 = λ2 = λ3
Similar to (ii), we obtain
| [1] | Hibbitt, H. D., Marcal, P. V., and Rice, J. R. A finite element formulation for problems of large strain and large displacement. International Journal of Solids and Structures, 6(8), 1069–1086(1970) |
| [2] | McMeeking, R. M. and Rice, J. R. Finite-element formulations for problems of large elasticplastic deformation. International Journal of Solids and Structures, 11(5), 601–616 (1975) |
| [3] | Hughes, T. J. R. Generalization of selective integration procedures to anisotropic and nonlinearmedia. International Journal of Numerical Methods Engineering, 15(9), 1413–1418 (1980) |
| [4] | Simo, J. C. and Pister, K. S. Remarks on the rate constitutive equations for finite deformationproblems: computational implications. Computer Methods in Applied Mechanics and Engineering,46(2), 201–215 (1984) |
| [5] | Rouainia, M. and Wood, D. M. Computational aspects in finite strain plasticity analysis of geotechnical materials. Mechanics Research Communications, 33(2),123–133 (2006) |
| [6] | Simo, J. C. and Meschke, G. A new class of algorithms for classical plasticity extended to finitestrains, application to geomaterials. Computational Mechanics, 11(4), 253–278 (1993) |
| [7] | Helnwein, P. Some remarks on the compressed matrix representation of symmetric second-orderand fourth-order tensors. Computer Methods in Applied Mechanics and Engineering, 190(22-23),2753–2770 (2001) |
| [8] | Simo, J. C. and Taylor, R. L. Quasi-incompressible finite elasticity in principal stretches: continuum basis and numerical algorithms. Computer Methods in Applied Mechanics and Engineering,85(3), 273–310 (1991) |
| [9] | Neto, S., Peric, D., and Owen, D. R. J. Computational Methods for Plasticity, John Wiley and Sons, New York, 590–599 (2008) |
| [10] | Simo, J. C. and Hughes, T. J. R. Computational Inelasticity, Springer-Verlag, Berlin, 143–147(1998) |
| [11] | Simo, J. C. 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, 99(1), 61–112 (1992) |
| [12] | Borja, R. I., Sama, K. M., and Sanz, P. F. On the numerical integration of three-invariant elastoplastic constitutive models. Computer Methods in Applied Mechanics and Engineering, 192(9-10),1227–1258 (2002) |
| [13] | Simo, J. C. and Taylor, R. L. Consistent tangent operators for rate-independent elastoplasticity.Computer Methods in Applied Mechanics and Engineering, 48(1), 101–118 (1985) |
| [14] | He, M. C., Xie, H. P., Peng, S. P., and Jiang, Y. D. Study on rock mechanics in deep mining engineering. Chinese Journal of Rock Mechanics and Engineering, 24(16), 2803–2813 (2005) |