Large deformation analysis of functionally graded elastoplastic materials via solid tetrahedral finite elements

Computers and Structures 146 (2015) 59–75

Contents lists available at ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

Large deformation analysis of functionally graded elastoplastic materials via solid tetrahedral finite elements João Paulo Pascon a,⇑, Humberto Breves Coda b,1 a b

Materials Engineering Department, Lorena School of Engineering, University of São Paulo, Pólo-Urbo Industrial, Gleba AI-6, s/no, 12602-810 Lorena, SP, Brazil Structural Engineering Department, São Carlos School of Engineering, University of São Paulo, Trabalhador Sãocarlense Avenue, 400, 13566-590 São Carlos, SP, Brazil

a r t i c l e

i n f o

Article history: Received 28 April 2014 Accepted 22 September 2014

Keywords: Functionally graded elastoplastic materials Finite deformations Finite strains Hyperelastoplasticity Isoparametric solid tetrahedral finite element of any-order

a b s t r a c t In this paper, a solid finite element formulation applied to the analysis of functionally graded materials (FGMs) under finite elastoplastic deformation with mixed hardening is presented. The main novelty of this paper is the use of gradually variable material coefficients in finite elastoplastic strain regime. The numerical simulations are performed with the same constitutive models but with different variations of the material coefficients. Full integration and high order tetrahedral elements are used. Results show that high order elements and mesh refinement avoids general locking problems. Finally, the differences (regarding final results) between the homogeneous and the FG cases are depicted. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction The aim of this study is to analyze, via solid finite elements, the behavior of functionally graded (FG) elastoplastic materials under finite deformations and strains. These materials have many engineering applications. At the automotive industry level, one can cite the sheet-metal forming processes, in which some automobile components acquire the desired shape via plastic deformation ([1–3]). Functionally graded materials (FGMs) are special composites in which the composition and, thus, the mechanical properties vary gradually (smooth and continuously) over the volume. This special variation avoids the material mismatch, which may cause delamination failures, common in laminated composites. The main applications of these advanced composites are in thermal environments, such as high-speed aircrafts, nuclear reactors and chemical power plants. Besides, highly deformable materials are very common nowadays, such as metals, which can present large displacements, and rubber-like materials, which can withstand large strains. In order to predict the structural behavior of engineering materials, and to reduce costs and time during the design process, scientists and engineers have been using the Finite Element Method (FEM). In this study, the element adopted is the isoparametric solid tetrahedral finite element of any-order (see, for

⇑ Corresponding author. Tel.: +55 12 3159 9900. 1

E-mail addresses: [email protected], [email protected] (J.P. Pascon). Tel.: +55 16 3373 9455.

http://dx.doi.org/10.1016/j.compstruc.2014.09.005 0045-7949/Ó 2014 Elsevier Ltd. All rights reserved.

instance, [4,5]). Regarding the numerical results, mesh refinement and full numerical integration are performed in order to obtain an accurate and reliable solution. The equilibrium analysis of elastoplastic materials under finite deformations is not a simple task. To analyze flexible structural components under large displacements, the geometrically nonlinear analysis is imperative. In addition, to describe the material behavior in the finite elastoplastic strain regime, a general framework has been recently developed. In this framework, called hyperelastoplascity, the finite elastoplastic deformations are described by means of a multiplicative gradient decomposition ([6–10]). In general, for elastoplastic materials under finite strains, the additive decomposition of the strain tensor, called Green-Naghdi decomposition [11], is not valid (see the work of [12] for further details). In this regime, the multiplicative decomposition, called Kröner-Lee decomposition [13,14], is more suitable and well-accepted. Some concepts from the small strain elastoplasticity can be extended to the finite strain regime, such as yield criterion, plastic flow, hardening, consistency condition and return algorithm. According to [15], FGMs have great potential in many engineering sectors, such as aerospace, automobile and defense industries, as well as electronics and biomedics. Concerning FG elastoplastic materials, some works present in the scientific literature may be cited. Akis and Eraslan [16] have presented plane strain analytical solutions of rotating FG hollow shafts. In [17], spherical pressure vessels composed of a FG elastoplastic material are investigated analytically. A finite element analysis of circular

60

J.P. Pascon, H.B. Coda / Computers and Structures 146 (2015) 59–75

Table 1 The expressions describing the plastic model. Internal dissipation inequality:

d ¼ ðMe  vÞ : Dp 

    @wp @wp    v jP0 @v @j

ð18Þ

Plastic spatial velocity gradient:

L p ¼ Dp þ W p

ð19Þ

Plastic strain rate tensor:

     Ep F1 ¼ symðLp Þ Dp ¼ FT p p

ð20Þ

Plastic spin tensor:

W p ¼ antðLp Þ

ð21Þ

Objective Jaumann rate of the intermediate backstress: 



v ¼ v Wp v þ vWp

ð22Þ

plates made of FG elastoplastic materials under low-velocity impact is performed in [18]. Rotating disks composed of FG elastoplastic materials are analytically and numerically analyzed in [19]. However, most of (if not all) the related works are restricted to small strain theory. Moreover there is no evidence of convergence analysis regarding mesh refinement in all consulted works related to the subject and, according to the research done by the present authors on the scientific literature, there are no published studies regarding finite element analysis of FG hyperelastoplastic materials. Therefore, another aim of this work is to fulfill the lack of studies related to the numerical analysis of FGMs under finite elastoplastic strains. This paper is organized in five sections. The exact kinematics and its numerical approximation are described in section two. In the third section, the constitutive framework, the adopted models, the FG laws and the equilibrium principle are given. Fourth section describes some aspects of the incremental procedure performed in order to achieve force equilibrium and plastic consistency. In the fifth section, the numerical examples used to validate the methodology, as well as the discussion of results, are provided. Finally, the main conclusions are given in the sixth section.

Three-dimensional von-Mises yield criterion:

2. Kinematics rffiffiffi 2 / ¼ /ðMe  v; jÞ ¼ kdev ðMe  vÞk  rj ðjÞ 6 0 3

ð23Þ

Associative plastic flow rule:





Dp ¼ k R P ¼ k

 dev ðM  vÞ @/ e ¼k @Me kdev ðMe  vÞk

ð24Þ

2.1. Elastoplastic decomposition

Swift isotropic hardening law:

rj ¼ Y 0  ðe0 þ jÞn

By kinematics one understands the relation among positions (or displacements) and the strain measure, including its elastic and plastic parts. In this section, the exact kinematics of hyperelastoplasticity and the adopted finite element approximation are described.

ð25Þ

In this study, finite elastoplastic strains are considered. So, the decomposition of the deformation gradient (or change of configuration function) into its elastic and plastic parts is performed here by means of the Kröner-Lee decomposition:

Isotropic hardening parameter:

F ¼ Fe Fp rffiffiffi rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rffiffiffi     2 2 2 2 kDp k ¼ j ¼ k rj ¼ ðDp : Dp Þ ¼ k ðRP : RP Þ ¼ k 3 3 3 3

ð26Þ

Nonlinear Armstrong-Frederick hardening law:







v ¼ k Rv ¼ cDp  k bv

ð27Þ

where F is the deformation gradient; and the subscripts ()e and ()p denote, respectively, its elastic and the plastic parts. Moreover, it is assumed the existence of an intermediate configuration, which is stress-free and locally defined. From expression (1), one can write:

J ¼ detðFÞ ¼ J e J p

Evolution of the initial backstress tensor:

E¼ 



X ¼ k RX ¼

 F1 p

v

FT p



symð2C1 p

ð1Þ

ð2Þ

1 T ðF F  IÞ ¼ FTp Ee Fp þ Ep 2

ð3Þ



Ep XÞ

ð28Þ

where J is the Jacobian; E is the Green–Lagrange strain tensor; and I is the identity matrix.

Null plastic spin tensor:

2.2. Tetrahedral finite element 



Wp ¼ O ) Fp ¼ k RP FP ;





v¼v

ð29Þ

Consistency condition:  

k/¼0

ð30Þ



 @/=@E :E ½ð@/=@Fp Þ : ðRP Fp Þ þ ð@/=@XÞ : RX þ ð@/=@ jÞrj 

x ¼ Xk /k ðnÞ or xi ¼ ðxi Þk /k ðnm Þ k

k

Plastic multiplier:

k¼

In the present study, the positional version of FEM is used (see, for instance, [4,5,20–22]). In this case, both initial and current positions, represented here by x and y, are mapped from a non-dimensional space (n):

y ¼ Y /k ðnÞ or yi ¼ ðyi Þ /k ðnm Þ k

ð31Þ

k

ð4Þ ð5Þ

where X and Y denote, in this order, the three initial and final coordinates of node k; and /k is the shape function associated with node k. The deformation gradient (1) is multiplicatively decomposed into two auxiliary gradients:

J.P. Pascon, H.B. Coda / Computers and Structures 146 (2015) 59–75

61

Fig. 1. Uniaxial extension of a block (deformation gradient, material coefficients and boundary conditions). The symbol u⁄ denotes the displacement along the x1-direction of the face with initial coordinate x1 = L0.

F ¼ F1 F1 0 @/k ðnÞ @/ ðnÞ F0 ¼ X or ðF 0 Þij ¼ ðxi Þk k @n @nj k @/k ðnÞ k @/k ðnÞ F1 ¼ Y or ðF 1 Þij ¼ ðyi Þ @n @nj k

ð6Þ

3. Constitutive law

ð7Þ

This section deals with the constitutive modeling of FGMs under finite elastoplastic strains. First, the general framework for homogeneous elastoplastic materials is briefly explained. After that, the specific constitutive models adopted are provided. The expressions, employed here to describe the gradual variation of the elastoplastic coefficients, are given in the sequence. Finally, the equilibrium principle adopted, as well as the numerical integration of the expressions involved, are given.

ð8Þ

As the finite element adopted here is the isoparametric solid tetrahedron, the non-dimensional space is defined here by the set:

n ¼ fðn1 ; n2 ; n3 Þ 2 R3 =0 6 n1 ; n2 ; n3 ; n1 þ n2 þ n3 6 1g

ð9Þ

In this study, shape functions are Lagrangian polynomials and the approximation order of elements is generic, i.e., the element degree may be any one. In order to employ this generalized approximation order, the numerical strategy presented in [4] is adopted here. The degrees of freedom of the present finite element are the current nodal positions (given by the vector Y).

3.1. Hyperelastoplasticity For hyperelastoplastic materials, the decomposition (1) is employed, and both elastic and plastic responses are derived from a scalar potential energy, which is – in the present study – the Helmholtz free energy (w). Considering elastic material isomorphism [7] and rate-independent elastoplasticity [23], one can write:

w ¼ wðE; Ep ; qÞ ¼ we ðEe Þ þ wp ðEp ; qÞ

ð10Þ

where q denotes hardening parameters. Following [7], the relation between stresses is: T 1 S ¼ F1 p Se Fp ¼ Fp

@we T F @Ee p

ð11Þ

where S is the second Piola–Kirchhoff stress tensor; and Se is the corresponding stress tensor defined at the intermediate configuration. The kinematic hardening parameter is the backstress tensor defined at the intermediate configuration (v), and the isotropic hardening is described by the parameter j. So, the yield criterion can be written as [7]:

/ ¼ /ðMe ; v; jÞ ¼ f MX ðMe  vÞ  f j ðjÞ 6 0 ðFTe Fe ÞSe

Me ¼ Ce Se ¼   @wp T F ¼ Fp XFTp v ¼ Fp @Ep p   v q¼

j

Fig. 2. Meshes of tetrahedral finite elements and histories of stretching employed.

ð12Þ ð13Þ ð14Þ ð15Þ

where / is the yield function; fMX and fj are auxiliary scalar functions; Me is the elastic Mandel stress tensor, defined at the intermediate configuration; Ce is the elastic right Cauchy-Green stretch tensor; and X is the backstress tensor defined at the initial configuration.

62

J.P. Pascon, H.B. Coda / Computers and Structures 146 (2015) 59–75

Fig. 3. Curves stress–strain for the block under uniaxial extension: (a) homogeneous case p = 0; (b) FGM with p = 0.25; (c) FGM with p = 4.0; (d) homogeneous case p ? 1.

3.2. Elastic models Regarding the elastic material response, two models for isotropic materials have been adopted: the linear Saint–Venant Kirchhoff hyperelastic model (SVK), and the nonlinear and compressible neoHookean hyperelastic model (nH). For these models, the specific expressions for the elastic Helmholtz free energy are:

h  i k ½trðEe Þ2 þ l tr E2e 2 K l 2 ¼ ½lnðJ e Þ þ ½trðCe Þ  3  2 lnðJ e Þ 2 2

ðwe ÞSVK ¼

ð16Þ

ðwe ÞnH

ð17Þ

In Table 1, the symbols sym(), ant() and dev() denote the symmetric part, the anti-symmetric part and the deviator operator, respectively. The coefficients Y0, e0 and n are the material coefficients of the Swift isotropic hardening law (25), for which the ini-

where k, l and K are material coefficients.

3.3. Plastic models The plastic response is usually defined by means of the evolution equations. In the context of Section 3.1, one must set the yield function (/), the rate of the plastic strain ðEp Þ, the rate of the back stress tensor ðvÞ and the rate of the isotropic hardening parameter  ðjÞ. In order to save space in the text, the expressions that define the plastic model adopted here are provided in Table 1. This plastic model is based on the works of Svendsen [7], Dettmer and Reese [9] and Vladimirov, Pietryga and Reese [10], in which further details regarding homogeneous hyperelastoplastic materials are given.

Fig. 4. Comparison of the graphs stress–strain for the block under uniaxial extension. The longitudinal engineering stress is the ratio between the applied force and initial cross-sectional area, and the linearized strain is the ratio between the length variation and the initial length.

J.P. Pascon, H.B. Coda / Computers and Structures 146 (2015) 59–75

63

Fig. 5. Distributions of the isotropic hardening parameter j along the x2-direction for a longitudinal stretch of: (a) k1 = 1.80; (b) k1 = 2.20; (c) k1 = 2.60; and (d) k1 = 3.0.

Fig. 6. Distributions of yield function / along the x2-direction for a longitudinal stretch of: (a) k1 = 1.80; and (b) k1 = 2.20.

64

J.P. Pascon, H.B. Coda / Computers and Structures 146 (2015) 59–75 Table 3 Final results for the cantilever beam composed of homogeneous high strength steel HSLA340.

Fig. 7. Cantilever metallic beam.

Table 2 Final results for the cantilever beam composed of homogeneous mild steel DC06. OPA

NDOF

NTE

NIP

F3

NI

PT

1

48 84 120 156 300 588 4365

18 36 54 72 144 288 4608

1

31.29 24.04 35.57 28.53 22.57 32.86 17.29

10,235 10,802 10,578 10,612 11,725 12,988 10,985

24 76 58 132 202 590 18,430

2

189 351 513 675

18 36 54 72

11

23.15 14.52 12.58 11.46

10,840 9950 10,137 10,263

408 749 1184 1312

3

480 912 1344 1776

18 36 54 72

15

10.20 9.36 9.01 8.58

10,273 10,084 10,021 10,444

975 1875 2775 3675

18 36 54 72

24

8.97 8.98 8.70 8.96

9725 9689 10,791 10,845

4

NTE

NIP

F3

NI

PT

48 84 120 156 300 588 4365

18 36 54 72 144 288 4608

1

34.39 28.52 54.73 44.08 30.31 43.65 30.04

9757 10,302 10,286 10,662 11,161 12,127 10,827

23 67 60 132 182 542 17,990

2

189 351 513 675

18 36 54 72

11

38.64 24.69 22.19 20.66

10,210 9154 9809 10,099

367 650 1078 1314

3

480 912 1344 1776

18 36 54 72

15

17.32 16.84 16.72 16.18

9420 10,060 10,087 10,732

1457 2125 3740 5481

4

975 1875 2775 3675

18 36 54 72

24

16.08 15.69 15.74 14.98

9770 9773 10,530 10,319

5576 9135 17,225 25,905

OPA = order of polynomial approximation. NDOF = number of degrees of freedom. NTE = number of tetrahedral finite elements. NIP = number of integration points (per element). F3 = final applied force along the x3-direction. NI = total number of iterations. PT = processing time (in seconds).

Table 4 Final results for the cantilever beam composed of a FG steel (DC06 and HSLA340). NDOF

NTE

NIP

F3

NI

PT

1

36 72 108 144 288 576

4

1662 2137 3724 5285

72 126 180 234 450 882

33.69 63.93 42.88 33.53 25.78 39.29

11,851 11,664 11,650 11,984 13,502 14,051

287 473 642 860 1870 3884

2

36 72 108 144

15

7030 9042 17,870 26,596

315 585 855 1125

23.60 18.89 16.14 14.70

11,438 10,920 10,777 10,659

1033 1937 3014 4230

3

840 1596 2352 3108

36 72 108 144

24

13.86 12.28 11.59 11.35

10,210 9852 10,281 11,060

3078 6859 12,533 20,503

4

1755 3375 4995 6615

36 72 108 144

45

11.96 11.44 11.31 11.18

10,741 9940 19,643 11,002

15,267 87,449 107,240 122,400

tial yield stress is r0 = Y0(e0)n. The letters b and c, in expression (27), denote the material coefficients regarding the nonlinear Armstrong-Frederick kinematic hardening law. The variables RP, Rv, RX and rj define, in this order, the plastic flow direction, the evolution of the initial and intermediate backstresses, and the evolution of the isotropic hardening parameter. 3.4. Functionally graded materials It is usual to apply FG laws in order to define the gradual variation of material mechanical properties. In this work three FG laws, written in terms of the initial coordinates, have been chosen: power law (P-FGM), sigmoidal law (S-FGM) and logarithmic law (L-FGM). Assuming that some material coefficient, represented by u, varies gradually along a general initial coordinate x, from u2  ¼ 0 to u1 at x  ¼ h, then the mathematical expressions of the at x adopted FG laws are:

 p  p

  x x þ u2  1  h h

NDOF

1

OPA

OPA = order of polynomial approximation. NDOF = number of degrees of freedom. NTE = number of tetrahedral finite elements. NIP = number of integration points (per element). F3 = final applied force along the x3-direction. NI = total number of iterations. PT = processing time (in seconds).

uPFGM ¼ u1 

OPA

ð32Þ

OPA = order of polynomial approximation. NDOF = number of degrees of freedom. NTE = number of tetrahedral finite elements. NIP = number of integration points (per element). F3 = final applied force along the x3-direction. NI = total number of iterations. PT = processing time (in seconds).

  

 

  s  s 1 x 1 1 x x 22 22 6 6 1 ¼ u1  1  þ u2  2 h 2 2 h h    s

 s

 1   1 1 x x x 2 2 uSFGM 0 6 6 ¼ u1  þ u2  1  2 2 h 2 h h

uSFGM







 x h

uLFGM ¼ u1  log10 9 þ 1



ð33:aÞ ð33:bÞ

     x þ u2  1  log10 9 þ 1 h ð34Þ

where p and s are the FG coefficients for the power and sigmoidal laws, respectively. The logarithmic law does not depend on FG coefficients. The first two FG laws are the most common in the scientific literature (see, for instance, [24]). One should note that we can have one or more material coefficients with gradual variation, we can have one material coefficient with gradual variation along one or more directions, and we can have different FG laws applied to differ-

J.P. Pascon, H.B. Coda / Computers and Structures 146 (2015) 59–75

65

evolution and, hence, the derivatives of these plastic variables regarding the strain tensor E vanish. In this case, the derivatives of the free energy w with respect to the Green–Lagrange strain tensor E are:

@w ¼S @E i @2w @S @ h ¼ ¼ ðFp Þ1  Se  ðFp ÞT @E@E @E @E n h io ¼ ðFp Þ1  }e : ðFp ÞT  II  ðFp Þ1  ðFp ÞT or ð}ep Þijmn   h i @Sij T 1 ðFp Þ1 ¼ ¼ ðFp Þ1 ð} Þ ðF Þ IIðF Þ p p e klop ik jl opmn @Emn

ð35Þ

}ep ¼

Fig. 8. Graphs force–displacement for the cantilever beam. For the homogeneous cases, the mesh has 72 fourth-order tetrahedral elements, and for the FG case, the mesh has 144 fourth-order tetrahedral elements.

ent material parameters. It is possible to propose other FG laws, but according to the definition, these new laws must be smooth and continuous. The FG laws adopted here can be used to describe the gradual variation of the elastic material coefficients k, l and K (see expressions (16) and (17)), and the plastic material coefficients Y0, e0, n, b and c (see expressions (25) and (27)). 3.5. Equilibrium In this study, to describe the equilibrium of forces at the current configuration, the Principle of the Stationary Total Potential Energy has been employed (see, for instance, the works of [4,5,20–22]). The resultant system of equilibrium equations is nonlinear regarding the current position y. In order to solve this system, the Newton–Raphson iterative technique is used here. The integrals that appear in the calculation of the internal force vector and Hessian matrix are numerically evaluated via tetrahedral Gaussian points. For FGMs, as the material parameters vary in space, more integration points are needed to evaluate these integrals when compared to the homogeneous cases. The position of these points and their corresponding weights can be found in [25,26]. The numerical formulation related to the aspects discussed in the present paragraph can be found in [4,5]. 4. Incremental procedure As the numerical strategy adopted here is iterative, the solution may diverge. Moreover, the integration of the plastic evolution equations cannot be performed analytically. Thus, this integration is done numerically. In this study, to improve the convergence of the solution and to reduce the errors related to the plastic evolution integration, the prescribed forces (or displacements) are divided into a finite number of incremental forces (or displacements), called steps. Besides at the equilibrium position, the current state of variables (F, Fp, j, X) must satisfy the yield criterion (23) and the consistency condition (30). In order to satisfy these conditions, the elastic predictor/plastic corrector algorithm is employed (for further operational details, see [27,28]). 4.1. Elastic prediction In general, the evolution of the free energy w depends on the evolution of the plastic variables (Fp, j, X) which, in turn, depends on the evolution of the total strain E. Thus, the derivatives of w in respect to E depends on the derivatives of the plastic variables in terms of E. At the beginning of each step, it is assumed that the plastic variables (Fp, j, X) remain constant, i.e., there is no plastic

}e ¼

@Se @ 2 we ¼ @Ee @Ee @Ee

ð36Þ

ð37Þ

where the fourth-order tensors }ep, }e and II are, in this order, the elastoplastic consistent tangent operator, the elastic consistent tangent operator and the unit tensor. 4.2. Plastic correction (return algorithm) Once the trial position is achieved, the elastic prevision of strain and stress prediction is performed. If yield criterion (23) is satisfied at all the points, the step is elastic and the numerical simulation can proceed to the next step. Otherwise, the plastic variables (Fp, j, X) must be updated in order to obtain / = / (F, Fp, j, X) = 0 respecting the evolution laws (see expressions 26–29). To perform this update, called return algorithm, the fully implicit backward Euler method is employed. In this study, the following residual vector is used:

rret

8 9 ðFp Þ  ðFp Þð0Þ  DkðRP ÞðFp Þ > > > > > > < = j  jð0Þ  Dkðrj Þ ¼ > > > X  Xð0Þ  DkðRX Þ > > > : ; /

ð38Þ

where the quantities with index (0) remain constant along the plastic correction and correspond to the values at the beginning of this correction; and the quantities without index correspond to the values at the end of the plastic correction. In this expression, one should note that Fp and X have nine components. The variables Fp, j, X and Dk are updated, via a Newton–Raphson iterative scheme, until the residual vector rret is sufficiently small. As the von Mises yield criterion is adopted here, the corresponding return algorithm is called radial return algorithm (see [27] for further details). 4.3. Elastoplastic prediction After performing the return algorithm, for the Gaussian points in which the plastic variables had to be updated, the derivatives of the plastic variables (Fp, j, X) with respect to the strain tensor E does not vanish. So, expression 1(36) is replaced by:

    1 @ FT p @S @ Fp T 1 @Se T 1 Se Fp þ Fp }hep ¼ ¼ or ð}hep Þijkl Fp þ Fp Se @E @E @E @E   1   @ðS Þ @Sij @ F p im e mn 1 ¼ ¼ ðSe Þmn ðF 1 ðF 1 p Þjn þ F p p Þjn im @Ekl @Ekl @Ekl   @ F T   p nj 1 þ Fp ðSe Þmn ð39Þ im @Ekl

66

J.P. Pascon, H.B. Coda / Computers and Structures 146 (2015) 59–75

Fig. 9. Final distribution of the isotropic hardening parameter for the cantilever beam composed of homogeneous mild steel DC06: (a) 36 elements of linear order; (b) 72 elements of linear order; (c) 36 elements of quadratic order; (d) 72 elements of quadratic order; (e) 36 elements of cubic order; (f) 72 elements of cubic order; (g) 36 elements of fourth- order; (h) 72 elements of fourth-order.

The partial derivative @(Fp)/@E can be found by means of the evolution laws and the consistency condition (see Table 1). 5. Numerical examples To validate the present formulation, some numerical examples have been analyzed. In order to avoid general locking problems and, thus, provide reliable results, convergence analysis regarding mesh refinement is carried out. Following a usual procedure, the same problem is analyzed with different constitutive models:

two homogeneous materials; and some intermediate FGMs. For the FG cases, the limiting values of the variable material coefficients correspond to the two homogeneous materials. So, in these cases, one needs to specify the spatial coordinate  x in which the material coefficients vary, as well as their limiting values (at  x¼ xmin and at  x¼ xmax ) and the used FG law (one of the Eqs. (32)–(35)). Regarding the numerical integration scheme, full integration is performed here, i.e., we do not use alternative techniques to avoid locking, such as reduced/selective integration, enhanced assumed strain method and hourglass modes stabilization.

J.P. Pascon, H.B. Coda / Computers and Structures 146 (2015) 59–75

67

Fig. 10. Final distribution of the isotropic hardening parameter for the cantilever beam composed of homogeneous high strength steel HSLA340: (a) 36 elements of linear order; (b) 72 elements of linear order; (c) 36 elements of quadratic order; (d) 72 elements of quadratic order; (e) 36 elements of cubic order; (f) 72 elements of cubic order; (g) 36 elements of fourth-order; (h) 72 elements of fourth-order.

5.1. Uniaxial extension In the first example, a block with dimensions L0, b0 and h0 is converted into a block with dimensions Lf = 3L0, bf and hf. The corresponding deformation gradient, the material coefficients and the boundary conditions are shown in Fig. 1. In this case, only isotropic hardening is considered, and the elastic Young modulus E and the Swift coefficient Y0 vary according to the FG power law (32) along the x2-direction. Four values of the FG power coefficient p have been selected: p = 0 (which corresponds to the homogeneous case E = E1 and Y0 = (Y0)1); p = 0.25; p = 4.0; and p ? 1 (which

corresponds to the homogeneous case E = E2 and Y0 = (Y0)2). The simulation was carried out via displacement control. Two meshes of tetrahedral finite elements and two histories of stretching are used (see Fig. 2). As expected, for each of the four FG power coefficients, both stretching histories provide the same final results (see Fig. 3). This fact is due to the material isomorphism considered here, i.e., the elastic response does not depend on the plastic strains [7]. Regarding the second history, the comparison among the four cases analyzed is shown in Fig. 4, in which one can see that increasing the power coefficient p decreases the longitudinal stresses during the

68

J.P. Pascon, H.B. Coda / Computers and Structures 146 (2015) 59–75

Fig. 11. Final distribution of the isotropic hardening parameter for the cantilever beam composed of FG steel (DC06 and HSLA340): (a) 72 elements of linear order; (b) 144 elements of linear order; (c) 72 elements of quadratic order; (d) 144 elements of quadratic order; (e) 72 elements of cubic order; (f) 144 elements of cubic order; (g) 72 elements of fourth- order; (h) 144 elements of fourth-order.

elastic regime and increases these stresses during the elastoplastic loading. The reason for these results is that, in this case, when we increase the FG coefficient p, we increase the volume fraction of the material whose elastic rigidity is smaller (E2) and initial yield stress is larger ((Y0)2). Some specific distributions of the isotropic hardening parameter and the yield function along the x2-direction are depicted in Figs. 5 and 6, in which one can see that, as expected, increasing the FG power coefficient p provides smaller values for the isotropic hardening parameter. In these figures, the points with U < 0 and j = 0 remain in elastic regime, and the points with U = 0 and j > 0 are in elastoplastic regime. Thus, for the functionally graded cases, once the material reaches the plastic limit, there

are an elastic portion and an elastoplastic portion along the cross-section. For these cases, as the x2 coordinate increases, the plastic limit decreases and the isotropic hardening parameter increases. 5.2. Prismatic cantilever beam In order to analyze the performance of the present finite element in bending-dominated problems with large displacements, a prismatic cantilever beam composed of a metallic material has been simulated (see Fig. 7). The simulation was carried out via displacement control and the necessary forces on the edge with pre-

69

J.P. Pascon, H.B. Coda / Computers and Structures 146 (2015) 59–75

Fig. 12. Partially loaded elastomeric block (geometry, loading and material coefficients).

Table 5 Vertical displacements of point A for the partially loaded block composed of homogeneous polymer (s = 0.0).

Table 7 Vertical displacements of point A for the partially loaded block composed of a FG polymer (s = 2.0).

OPA

NDOF

NE

NIP

u1 (A) [100]

u1 (A) [300]

NI

PT

OPA

NDOF

NE

NIP

u1 (A) [100]

u1 (A) [300]

NI

PT

1

54 225 588 1215

24 192 648 1536

1

0.67 0.99 1.14 1.22

0.73 1.06 1.20 1.27

2400 2400 2400 2400

7 50 171 468

1

81 375 1029 2187

48 384 1296 3072

1

1.12 1.54 1.72 1.81

1.30 1.72 1.87 1.93

2426 2750 2897 2958

38 376 1583 4864

2

225 1215 3549

24 192 648

11

1.33 1.36 1.37

1.37 1.39 1.39

2462 2491 2412

77 747 4726

2

375 2187 6591

48 384 1296

11

1.98 1.99 1.99

2.06 2.04 2.04

2999 3069 3036

267 3642 27,637

3

588 3549

24 192

15

1.34 1.36

1.41 1.39

2927 2848

325 5129

3

1029 6591

48 384

15

1.96 1.97

2.02 2.03

3841 3322

38,306 30,692

4

1215 7803

24 192

24

1.37 1.37

1.39 1.39

2400 2405

1204 29,051

4

2187 14,739

48 384

24

1.99 1.99

2.04 2.04

3086 3111

5211 216,000

OPA = order of polynomial approximation. NDOF = number of degrees of freedom. NE = number of elements. NIP = number of integration points. u1 (A) [X] = displacement along the x1-direction of point A at the end of the step X. NI = number of iterations. PT = processing time (in seconds).

Table 6 Vertical displacements of point A for the partially loaded block composed of a FG polymer (s = 0.50).

OPA = order of polynomial approximation. NDOF = number of degrees of freedom. NE = number of elements. NIP = number of integration points. u2 (A) [X] = displacement along the x1-direction of point A at the end of the step X. NI = number of iterations. PT = processing time (in seconds).

Table 8 Vertical displacements of point A for the partially loaded block composed of homogeneous glassy polymer OPET (s ? 1). OPA

NDOF

NE

NIP

u1 (A) [100]

u1 (A) [300]

NI

PT

43 337 1448 4427

1

54 225 588 1215

24 192 648 1536

1

1.29 2.01 2.49 2.79

1.55 2.57 3.15 3.44

2463 2932 3152 3282

7 73 266 787

2926 2918 2908

254 3398 26,247

2

225 1215 3549

24 192 648

11

3.18 3.48 3.56

3.72 3.90 3.93

3386 3571 3618

132 1323 8186

1.98 1.98

3945 3418

39,671 31,535

3

588 3549

24 192

15

3.53 3.58

4.56 4.14

4224 4228

529 8518

1.99 1.99

2943 2988

4895 216,000

4

1215 7803

24 192

24

3.58 3.59

3.95 3.91

3571 3809

1966 48,493

OPA

NDOF

NE

NIP

u1 (A) [100]

u1 (A) [300]

NI

PT

1

81 375 1029 2187

48 384 1296 3072

1

1.11 1.50 1.67 1.76

1.29 1.67 1.81 1.88

2400 2578 2743 2814

2

375 2187 6591

48 384 1296

11

1.93 1.94 1.94

2.00 1.99 1.99

3

1029 6591

48 384

15

1.92 1.92

4

2187 14,739

48 384

24

1.94 1.94

OPA = order of polynomial approximation. NDOF = number of degrees of freedom. NE = number of elements. NIP = number of integration points. u2 (A) [X] = displacement along the x1-direction of point A at the end of the step X. NI = number of iterations. PT = processing time (in seconds).

OPA = order of polynomial approximation. NDOF = number of degrees of freedom. NE = number of elements. NIP = number of integration points. u1 (A) [X] = displacement along the x1-direction of point A at the end of the step X. NI = number of iterations. PT = processing time (in seconds).

70

J.P. Pascon, H.B. Coda / Computers and Structures 146 (2015) 59–75

Fig. 13. Force versus displacement for the partially loaded block. The coefficient s is the sigmoidal FG coefficient (see Eq. 33). The mesh employed here for the homogeneous cases has 2601 nodes (7803 degrees of freedom) and 192 fourthorder elements, and the mesh for the FG cases has 4913 nodes (14,739 degrees of freedom) and 384 fourth-order elements.

scribed displacement were computed at the end of some steps. For this numerical example, three material models have been selected: the homogeneous mild steel DC06; the homogeneous high strength steel HSLA340; and the FGM whose coefficients vary between the homogeneous limits according to the logarithmic

law (34) along the x3-direction. For the homogeneous cases, the tetrahedral meshes employed have one layer of elements along the x3-direction and only one half of the beam was discretized due to the symmetry regarding the plane x3 = 0.0. To simulate the beam with the FGM, two layers of elements along the x3-direction have been employed. In order to provide reliable results, convergence analysis was carried out regarding the approximation order and the number of tetrahedral elements in the longitudinal direction. For the three material models employed, the results in terms of the applied forces converge with mesh refinement (see Tables 2–4). Some additional information – as the number of degrees of freedom, processing time and number of iterations – is also provided, so future finite element formulations can be compared to the present one. According to Tables 2–4, considering the meshes employed, as the number of degrees of freedom NDOF increases, the applied force F3 decreases, except for the linear tetrahedron elements, which present an oscillatory behavior regarding the values of F3 and NDOF. The resultant processing time PT becomes very high for the most refined meshes. The total number of iterations NI is not significantly affected by the mesh size. For the most refined meshes (the most refined for the homogeneous material and the most refined for the FGM), the comparison in terms of the graph force–displacement is depicted in Fig. 8. Finally, the distribution of the isotropic hardening parameter at the last step, for some specific meshes, is illustrated in Figs. 9–11. As this parameter is obtained at each Gauss point during the simulation, the corre-

Fig. 14. Vertical displacements of the partially loaded block at the end of step 100: (a) homogeneous polymer (s = 0.0); (b) FG polymer (s = 0.50); (c) FG polymer (s = 2.0); (d) homogeneous glassy polymer OPET (s ? 1).

J.P. Pascon, H.B. Coda / Computers and Structures 146 (2015) 59–75

71

Fig. 15. Deviatoric backstress norm (40) for the partially loaded block composed of glassy polymer (s ? 1) at the end of step 100: (a) 648 elements of linear order; (b) 1536 elements of linear order; (c) 192 elements of quadratic order; (d) 648 elements of quadratic order; (e) 24 elements of cubic order; (f) 192 elements of cubic order; (g) 24 elements of fourth-order; (h) 192 elements of fourth-order.

72

J.P. Pascon, H.B. Coda / Computers and Structures 146 (2015) 59–75

Fig. 16. Deviatoric backstress norm (47) for the partially loaded block composed of FG polymer (s = 0.50) at the end of step 100: (a) 1296 elements of linear order; (b) 3072 elements of linear order; (c) 384 elements of quadratic order; (d) 1296 elements of quadratic order; (e) 48 elements of cubic order; (f) 384 elements of cubic order; (g) 48 elements of fourth-order; (h) 384 elements of fourth-order.

J.P. Pascon, H.B. Coda / Computers and Structures 146 (2015) 59–75

sponding values at the nodes, used in the post-processing program, are determined via a least squares method, in which a linear extrapolation is employed. As one can see in Figs. 8–11, the beam composed of homogeneous high strength steel (HSLA340) has the largest plastic limit and provides the largest values of the isotropic hardening parameter. The reason for these facts is that the material coefficients adopted (see Fig. 7) provide a larger initial yield stress and a smaller derivative d//dj for the high strength steel. Finally, as expected, the FG beam has an intermediate behavior, and the plastic region is around the hinge (x1 = 0) due to the stress concentration in this area. 5.3. Partially loaded elastomeric block The last example is a full three-dimensional problem involving large compressive strains, large mesh distortions and an elastomeric FGM under elastoplastic regime (see Fig. 12). This problem is adapted from the reference [9], which employed a finite-strain elastoplastic framework for homogeneous materials, and a linear solid finite element based on reduced integration with hourglass stabilization. To describe the mechanical behavior, four material models have been selected: two homogeneous and two FGMs. For the FG cases, the sigmoidal law (33) was employed, and the values of the coefficient s are: 0.50 and 2.0. The limiting material property values at face x1 = 0.0 are the adopted values in the homogeneous analysis done by Dettmer and Reese [9]. Only kinematic hardening is considered in this example. Due to the double symmetry regarding the planes x2 = 10.0 and x3 = 10.0, one quarter of

73

the block is discretized. In this case, the simulation was performed via a load control. Regarding the mesh division, the number of layers along the x3-direction for the FG case is always the double of the number of layers for the homogeneous case. For example, for the least refined meshes, we have two layers along the x1- and x2-directions, and along the x3-direction, we have one layer for the homogeneous case and two layers for the FG case. This was done because, for the FG cases, the material coefficients vary along the x3-direction. The convergence analysis regarding the vertical displacement of point A is provided in Tables 5–8 for all the cases analyzed. One can note that a large number of elements of linear order is necessary to provide displacement results close to the solution. As in the last example, considering the meshes employed, as the number of degrees of freedom NDOF increase, the absolute values of the vertical displacement of point A increase, the processing time PT becomes very high and the total number of iterations NI does not change significantly (see Tables 5–8). For the most refined meshes employed here (the most refined for the homogeneous material and the most refined for the FGM), the force–displacement solutions, as well as the agreement with the reference solution, are depicted in Fig. 13. The differences among the four materials analyzed regarding the vertical displacements are depicted in Fig. 14, in which one can see that the homogeneous glassy polymer (OPET) is the most flexible. As the kinematic hardening is included, the initial backstress tensor (X) has been computed. The convergence analysis with respect to this tensor is performed via the following norm:

Fig. 17. Elastic Green–Lagrange strain (Ee)11 of the partially loaded block at the end of step 100: (a) homogeneous polymer (s = 0.0); (b) FG polymer (s = 0.50); (c) FG polymer (s = 2.0); (d) homogeneous glassy polymer OPET (s ? 1).

74

J.P. Pascon, H.B. Coda / Computers and Structures 146 (2015) 59–75

Fig. 18. Plastic Green–Lagrange strain (Ep)11 of the partially loaded block at the end of step 100: (a) homogeneous polymer (s = 0.0); (b) FG polymer (s = 0.50); (c) FG polymer (s = 2.0); (d) homogeneous glassy polymer OPET (s ? 1).

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi



ffi 1 1 X  ðtrXÞI : X  ðtrXÞI 3 3

The distribution of this norm at the end of step 100 is shown in Figs. 15 and 16, for one homogeneous case and one FG case. These figures reveal that, in terms of the norm (40), the meshes with linear order present severe locking even for a large number of elements. As one can see, in all the cases the backstresses are concentrated at the region near the applied pressure. Regarding the most refined meshes, to illustrate the differences among the four materials analyzed in respect to the strain field, the vertical normal strains at the end of the step 100 are showed in Figs. 17 and 18, in which one can see that increasing the sigmoidal coefficient s increases the strains. Moreover, the results presented in this section show that the numerical solutions for s = 0.50 and for s = 2.0 are very close. This fact may be due to the symmetry of the sigmoidal law (33). Again, the FG cases provide an intermediate behavior. Finally, for the most rigid case (s = 0.0), the block remains in the elastic regime (see Figs. 13 and 18a).

the standard isoparametric solid tetrahedron of any-order. The modeling of the FG variation of the mechanical properties is done via smooth and continuous FG laws. The solution of the incremental elastoplastic problem is based on the elastic prediction/plastic correction method. In order to validate the present formulation, three numerical examples have been analyzed. A simple problem is simulated to check the influence of the material model on the final results, regarding displacements, stresses and plastic parameters. The performance of the finite element (up to fourth-order) is tested in two more elaborated problems: a prismatic cantilever beam made of a FG metal; and a partially loaded block composed of a FG polymer. High order elements are proved to be free of locking even for FGMs. Based on the realized bibliographic research there is no other finite elements designed for this kind of analysis; therefore, this study can be used as a basis for future comparisons with other finite element formulations, and for future implementations of other elastoplastic models. Future studies related to the behavior of the proposed formulation regarding mesh distortion are the next objectives of the authors.

6. Conclusions

Acknowledgements

A finite element formulation to solve general solids composed of functionally graded (FG) elastoplastic materials under finite strains is presented. The finite-strain elastoplastic framework, called hyperelastoplastic, is based on the multiplicative gradient decomposition and on the Helmholtz free-energy function, and accounts for finite elastic and plastic strains. The finite element is

The authors appreciate the essential help of the Materials Engineering Department (DEMAR/EEL/USP) and the Structural Engineering Department (SET/EESC/USP), and the Brazilian agencies Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP) and Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES).

kdev ðXÞk ¼

ð40Þ

J.P. Pascon, H.B. Coda / Computers and Structures 146 (2015) 59–75

References [1] Hol J, Alfaro MVC, Rooij MB, Meinders T. Advanced friction modeling for sheet metal forming. Wear 2012;286–287:66–78. [2] Parsa MH, Ahkami SN, Pishbin H, Kazemi M. Investigating spring back phenomena in double curved sheet metals forming. Mater Des 2012;41:326–37. [3] Neto DM, Oliveira MC, Menezes LF, Alves JL. Improving Nagata patch interpolation applied for tool surface description in sheet metal forming simulation. Comput-Aid Des 2013;45:639–56. [4] Pascon JP, Coda HB. Analysis of elastic functionally graded materials under large displacements via high-order tetrahedral elements. Finite Elements Anal Des 2012;50:33–47. [5] Pascon JP, Coda HB. High-order tetrahedral finite elements applied to large deformation analysis of functionally graded rubber-like materials. Appl Math Modell 2013;37:8757–75. [6] Dogui A, Sidoroff F. Kinematic hardening in large elastoplastic strain. Eng Fract Mech 1985;21:685–95. [7] Svendsen B. A thermodynamic formulation of finite-deformation elastoplasticity with hardening based on the concept of material isomorphism. Int J Plast 1998;16:473–88. [8] Ekh M, Runesson K. Modeling and numerical issues in hyperelasto-plasticity with anisotropy. Int J Solids Struct 2001;38:9461–78. [9] Dettmer W, Reese S. On the theoretical and numerical modelling of ArmstrongFrederick kinematic hardening in the finite strain regime. Comp Meth Appl Mech Eng 2004;193:87–116. [10] Vladimirov IN, Pietryga MP, Reese S. On the modelling of non-linear kinematic hardening at finite strains with application to springback – comparison of time integration algorithms. Int J Numer Meth Eng 2008;75:1–28. [11] Green AE, Naghdi PM. A general theory of an elastic-plastic continuum. Arch Rational Mech Anal 1965;18:251–81. [12] Eterovic AL, Bathe KJ. A note on the use of the additive decomposition of the strain tensor in finite deformation inelasticity. Comp Meth Appl Mech Eng 1991;93:31–8. [13] Kröner E. Allgemeine kontinuumstheorie der versetzungen und eigenspannungen. Arch Rational Mech Anal 1960;4:273–334.

75

[14] Lee EH. Elastic-plastic deformation at finite strain. ASME J Appl Mech 1969;36:1–6. [15] Ben-Oumrane S, Abedlouahed T, Ismail M, Mohamed B, Mustapha M, El Abbas AB. A theoretical analysis of flexional bending of Al/Al2O3 S-FGM thick beams. Comput Mater Sci 2009;44:1344–50. [16] Akis T, Eraslan AN. Exact solution of rotating FGM shaft problem in the elastoplastic state of stress. Arch Appl Mech 2007;77:745–65. [17] Akis T. Elastoplastic analysis of functionally graded spherical pressure vessels. Comput Mater Sci 2009;46:545–54. [18] Gunes R, Aydin M, Apalak MK, Reddy JN. The elasto-plastic impact analysis of functionally graded circular plates under low-velocities. Compos Struct 2011;93:860–9. [19] Jahromi BH, Hashemi HN, Vaziri A. Elasto-plastic stresses in a functionally graded rotating disk. J Eng Mater Tech 2012;134:1–11. [20] Coda HB, Greco MA. A simple FEM formulation for large deflection 2D frame analysis based on position description. Comp Meth Appl Mech Eng 2004;193:3541–57. [21] Coda HB, Paccola RR. An alternative positional FEM formulation for geometrically non-linear analysis of shells: curved triangular isoparametric elements. Comput Mech 2007;40:185–200. [22] Coda HB. A solid-like FEM for geometrically non-linear 3D frames. Comput Methods Appl Mech Eng 2009;198:3712–22. [23] Simo JC, Taylor RL. Consistent tangent operators for rate-independent elastoplasticity. Comp Meth Appl Mech Eng 1985;48:101–18. [24] Chi SH, Chung YL. Mechanical behavior of functionally graded material plates under transverse load – Part II: Numerical results. Int J Solids Struct 2006;43:3675–91. [25] Jinyun Y. Symmetric Gaussian quadrature formulae for tetrahedral regions. Comp Meth Appl Mech Eng 1984;43:349–53. [26] Keast P. Moderate-degree tetrahedral quadrature formulas. Comp Meth Appl Mech Eng 1986;55:339–48. [27] Simo JC, Hughes TJR. Computational inelasticity. New York: Springer; 1998. [28] Crisfield MA. Non-linear finite element analysis of solids and structures, vol. 1. Chichester, England: John Wiley & Sons; 2000.