- Email: [email protected]
Three-dimensional finite element analysis of finite deformation micromorphic linear isotropic elasticity
International Journal of Engineering Science 49 (2011) 1326–1336
Contents lists available at ScienceDirect
International Journal of Engineering Science journal homepage: www.elsevier.com/locate/ijengsci
Three-dimensional finite element analysis of finite deformation micromorphic linear isotropic elasticity Volkan Isbuga, Richard A. Regueiro ⇑ Department of Civil, Environmental, and Architectural Engineering, University of Colorado at Boulder, United States
a r t i c l e
i n f o
Article history: Received 11 October 2010 Accepted 19 April 2011 Available online 12 May 2011 Keywords: Micromorphic elasticity Finite deformation Finite element implementation Three dimensional
a b s t r a c t A finite deformation micromorphic materially linear isotropic elastic model is formulated and implemented for three dimensional finite element analysis. The model is based on the kinematics, balance equations and thermodynamic equations proposed by Eringen and Suhubi (1964). The constitutive equations are calculated in the reference configuration, and the resulting stresses are mapped to the current configuration. The balance of linear momentum and the balance of first moment of momentum are linearized to construct the consistent tangent for three dimensional finite element implementation for solution by the Newton–Raphson method. Three dimensional numerical examples are analyzed to demonstrate preliminarily the implementation. Ó 2011 Elsevier Ltd. All rights reserved.
1. Introduction The nonlinear theory of a micromorphic continuum has been first proposed by Eringen and Suhubi (1964) and Suhubi and Eringen (1964). The micromorphic continuum in which each material point is endowed with nine additional kinematic degrees of freedom (dof) represents the most general case of higher order continua of ‘‘grade one’’ (Eringen, 1999) (i.e., condef sidering only the first gradient of micro-structure deformation v ¼ @n=@N) (Eringen & Suhubi, 1964; Green & Naghdi, 1995; Kröner, 1968). These additional dof allow modeling material response when classical continuum theories may not be sufficient, such as when size effects of the micro-structure become important under large deformation gradient (Chambon, Caillerie, & Tamagnini, 2004). For example, in particulate materials, strain localization has a size effect related to the grain-size, and thus it is desirable to capture this size effect using multiscale models. Most of the micromorphic continuum based models follow two main approaches for the kinematic assumptions and balance equations which were proposed by Mindlin (1964) and Eringen and Suhubi (1964). Germain (1973) also provided an approach in determining the balance equations of generalized continua by implementing the principle of virtual work. Neff and Forest (2007) applied a geometrically exact micromorphic model with a homogenization approach for elastic metallic foams. Forest and Sievert (2006) provided some guidelines based on previous literature to choose the proper higher order continuum or higher grade continuum in modeling length-scale-dependent material response. Sansour, Skatulla, and Zbib (2010) developed a micromorphic model for finite strain inelasticity, and they applied the model to metals. Numerous works on micromorphic continuum based modeling of material response with elastic and plastic models are still on-going, and it is not possible to mention all such works here. However, a summary and development of the micro-continuum theories may be found in Eringen (1999). In this work, we implemented a finite deformation micromorphic materially linear isotropic elastic model to simulate three dimensional problems using the finite element method, which establishes the basis for future elastoplastic model ⇑ Corresponding author. Fax: +1 303 492 7317. E-mail address: [email protected] (R.A. Regueiro). 0020-7225/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijengsci.2011.04.006
V. Isbuga, R.A. Regueiro / International Journal of Engineering Science 49 (2011) 1326–1336
1327
implementation (Regueiro, 2009, 2010). The model follows the kinematics and balance equations proposed by Eringen and Suhubi (1964), Suhubi and Eringen (1964) and Eringen (1999). One of the main reasons that we follow this approach is that all macro-scale stress tensors and higher order couple stress tensor together with body force vectors and body force couple tensors are expressed in terms of micro-scale tensors and parameters (Eringen, 1999). This provides an opportunity to use the model in a multiscale modeling approach where a bridging/coupling between scales may be needed. Also, the nonlinear continuum mechanics of the formulation are developed rigorously by Eringen and co-workers. We assume Catersian coordinates for continuum finite element implementation. Index notation will be used throughout so as to be as clear as possible with regard to details of the formulation. Some symbolic/direct notation is also given, such as (ab)ik = aijbjk. Boldface denotes a tensor or vector, where its index notation has been given uniformly throughout the paper. Generally, variables in uppercase letters live in the reference configuration B0 (such as the reference differential volume dV), and variables in lowercase live in the current configuration B (such as the current differential volume dv). The same applies to their indices, such that a differential line segment in the current configuration dxi is related to a differential line segment in the reference configuration dXI through the deformation gradient: dxi = FiIdXI (Einstein’s summation convention assumed (see Eringen, 1962; Holzapfel, 2000)). Subscripts (),i and (),I imply partial differentiation in the current and reference configurations, respectively. A superscript prime symbol ()0 denotes a variable associated with the micro-element (see Section 2 _ ¼ DðÞ=Dt denotes material time derivative. on kinematics). Superposed dot ðÞ An outline of the remainder of the paper is as follows: Section 2 summarizes the kinematics of a micromorphic continuum; Section 3 presents the momentum balance equations as well as Clausius–Duhem inequality in the current configuration; and Section 4 shows the mapping of Clausius–Duhem inequality to reference configuration in which constitutive equations are derived. Constitutive equations are given in Section 5. Weak form and finite element equations are given in Section 6. Choice of elastic parameters and numerical results are presented in Section 7, and conclusions in Section 8. 2. Micromorphic kinematics Fig. 1 illustrates the mapping of the macro-element and micro-element in the reference configuration B0 to the current configuration B through the deformation gradient F and micro-deformation tensor v. The macro-element continuum point is denoted by P(X, N) and p(x, n, t) in the reference and current configurations, respectively, with centroid C and c. The microelement continuum point centroid is denoted by C0 and c0 in the reference and current configurations, respectively. The micro-element in general represents a grain/particle/fiber micro-structural sub-volume (cluster) of the heterogeneous material. The relative position vector of the micro-element centroid with respect to the macro-element centroid is denoted by N and n(X, N, t) in the reference and current configurations, respectively, such that the micro-element centroid position vectors are written as (Fig. 1)
X 0K ¼ X K þ NK ;
x0k ¼ xk ðX; tÞ þ nk ðX; N; tÞ
ð1Þ
Eringen and Suhubi (1964) assumed that for sufficiently small lengths kNk 1 (kk is the L2 norm), n is linearly related to N through the micro-deformation tensor v, such that
nk ðX; N; tÞ ¼ vkK ðX; tÞNK
ð2Þ
where then the spatial position vector of the micro-element centroid is written as
x0k ¼ xk ðX; tÞ þ vkK ðX; tÞNK
Fig. 1. Deformation of macro differential volume dV and micro differential volume dV0 (Eringen, 1968).
ð3Þ
1328
V. Isbuga, R.A. Regueiro / International Journal of Engineering Science 49 (2011) 1326–1336
This is called a micromorphic continuum theory of ‘‘grade one’’ (Eringen, 1999). This is equivalent to assuming an affine, or homogeneous, deformation of the macro-element differential volume dV (but not necessarily the body B; i.e., the continuum body B may experience heterogeneous deformation because of v, even if boundary conditions (BCs) are uniform). It also simplifies considerably the formulation of the micromorphic continuum balance equations as presented in Eringen and Suhubi (1964) and Eringen (1999). This micro-deformation tensor v is analogous to the small strain ‘‘micro-deformation’’ tensor w in Mindlin (1964), physically described in Mindlin’s Fig. 1. Eringen (1968) also provides a physical interpretation of v generally, but then simplifies for the micropolar case. For example, v can be interpreted as calculated from a micro-displacement tensor U as v = 1 + U, where U is not calculated from a micro-displacement vector u0 , because in general U, and in turn v, are incompatible. For finite element implementation, U will be interpolated at the nodes, providing an additional nine dofs because it is unsymmetric. 3. Micromorphic balance equations and Clausius–Duhem inequality In this section, the balance of linear momentum, first moment of momentum, and Clausius–Duhem inequality are presented in the current configuration B. The derivation of the equations and the detailed integral definitions of all the vectors and the stress tensors as well as the higher order stress and couple tensors can be found in Eringen and Suhubi (1964) and Eringen (1999), and therefore will not be given here. Eq. (4) is the balance of linear momentum in a classical Cauchy continuum theory, however, the Cauchy stress tensor ris now unsymmetric. Eq. (5) was first introduced by Eringen and Suhubi (1964), which is called the balance of first moment of momentum that is more general than the balance of angular momentum. It contains six additional equations, for a total of nine equations for nine unknowns U. Eq. (6) is the Clausius–Duhem inequality in the current configuration B. The equations are written as
rlk;l þ qðfk ak Þ ¼ 0
ð4Þ
rml sml þ mklm;k þ qðklm xlm Þ ¼ 0
ð5Þ
1 q w_ þ gh_ þ qk h;k þ rkl v l;k þ ðskl rkl Þmlk þ mklm mlm;k P 0 h
ð6Þ
where f and a are the body force vector per unit mass and the acceleration vector respectively, s the micro-stress tensor, m the higher order couple stress tensor, k the body force couple tensor per unit mass, x the micro-spin inertia tensor per unit _ 1 the micro-gyration tensor, $v the velocity gradient, q the heat flux vector per unit mass, w the Helmholtz free mass, m ¼ vv energy function per unit mass, q the mass density in the current configuration, and h the temperature. In the derivation of constitutive equations, we map these equations into the reference configuration in which constitutive equations are more naturally introduced, such as to incorporate anisotropy in a finite strain setting and/or to introduce plasticity in a naturally elastically-relaxed intermediate configuration (Regueiro, 2010). 4. Clausius–Duhem inequality in B0 The Clausius–Duhem inequality given in (6) expressed in current configuration B can be rewritten in the reference configuration B0 by using the relations q = q0/J and dv = J dV, as well as Nanson’s formula and the Piola transform in local form as
1 q0 w_ þ gh_ þ J rkl ðv l;k mlk Þ þ Jskl mlk þ J mlm;k mklm þ Q K h;K P 0 h
ð7Þ
where the Piola transforms for the unsymmetric Cauchy stress tensor rkl, symmetric micro-stress tensor skl, and higher order couple stress tensor mklm are written as
1 J
1 J
1 J
rkl ¼ F kK SKL F lL ; skl ¼ F kK RKL F lL ; mklm ¼ F kK F lL MKLM vmM
ð8Þ
Inserting (8) into (7) gives us the Clausius–Duhem inequality in the reference configuration as
1 q0 w_ þ gh_ þ F kK SKL F lL ðv l;k mlk Þ þ F kK RKL F lL mlk þ mlm;k F kK F lL M KLM vmM þ Q K h;K P 0 h
ð9Þ
The Clausius–Duhem inequality in conjunction with the Helmholtz free energy function, which is presented in the next section, is applied to derive the elastic constitutive equations in the reference configuration B0 . The rest of the paper focuses on the Helmholtz free energy function assumption, the derivation of constitutive equations, weak form, finite element formulation and linearization, and numerical examples.
V. Isbuga, R.A. Regueiro / International Journal of Engineering Science 49 (2011) 1326–1336
1329
5. Helmholtz free energy function and constitutive equations Eringen and Suhubi (1964) proposed unique sets of elastic deformation measures which are invariant with respect to rigid body motion in the current configuration. For the Helmholtz free energy function per unit mass in B0 , we use the set consisting of C, W, and C defined as
C KL ¼ F iK F iL ;
WKL ¼ F iK viL ;
CKLM ¼ F iK viL;M
ð10Þ
These deformation tensors can be written in terms of the displacement vector U and micro-displacement tensor U as follows (Eringen, 1968; Suhubi & Eringen, 1964)
C KL ¼ dKL þ U K;L þ U L;K þ U M;K U M;L
ð11Þ
WKL ¼ dKL þ UKL þ U L;K þ U M;K UML
ð12Þ
CKLM ¼ UKL;M þ U N;K UNL;M
ð13Þ
Then, the Helmholtz free energy density in terms of these deformation measures is expressed as
q0 wðC; W; CÞ; q0 wðC KL ; WKL ; CKLM Þ
ð14Þ
A quadratic form assumption for the Helmholtz free energy function in the reference configuration will lead to materially linear elasticity at finite strain as
1 2
1 2
1 2
q0 w ¼ EKL AKLMN EMN þ E KL BKLMN E MN þ CKLM C KLMNPQ CNPQ þ EKL DKLMN E MN
ð15Þ
where elastic strain tensors are defined in Suhubi and Eringen (1964) as
EKL ¼
1 ðC KL dKL Þ; 2
E KL ¼ WKL dKL
ð16Þ
AKLMN, BKLMN, DKLMN, and CKLMNPQ are the elastic material moduli tensors (Regueiro, 2010) defined for isotropic linear elasticity and obtained from back analysis of linear constitutive equations and energy function given in Suhubi and Eringen (1964). They are written as
AKLMN ¼ kdKL dMN þ lðdKM dLN þ dKN dLM Þ BKLMN ¼ ðg sÞdKL dMN þ ðj rÞdKM dLN þ ðm rÞdKN dLM C KLMNPQ ¼ s1 ðdKL dMN dPQ þ dKQ dLM dNP Þ þ s2 ðdKL dMP dNQ þ dKM dLQ dNP Þ þ s3 dKL dMQ dNP þ s4 dKN dLM dPQ þ s5 ðdKM dLN dPQ þ dKP dLM dNQ Þ
ð17Þ
þ s6 dKM dLP dNQ þ s7 dKN dLP dMQ þ s8 ðdKP dLQ dMN þ dKQ dLN dMP Þ þ s9 dKN dLQ dMP þ s10 dKP dLN dMQ þ s11 dKQ dLP dMN DKLMN ¼ sdKL dMN þ rðdKN dLM þ dLN dKM Þ where major and minor symmetry hold for AKLMN and DKLMN, whereas only major symmetry holds for the tensors BKLMN and CKLMNPQ, k, l, g, s, j, m, r, s1, . . . , s11 are the eighteen micromorphic isotropic elastic parameters. Note that the units for s1. . .s11 are stress length2 (Pa m2), thus there is a built in length scale to these elastic parameters for the higher order stress. These elastic material moduli tensors are different than those of elastic material moduli tensors given in Eringen (1999) due to different strain measures in our formulation, except the higher order elastic modulus tensor CKLMNPQwhich has the same strain measures for the higher order couple stress tensor in Eringen (1999). Following the classical argument of Coleman and Noll (1963) and Coleman and Gurtin (1967) for independent rate pro_ Wand _ _ the constitutive equations in terms of partial derivatives of Helmholtz free energy function with recesses for C; C, spect to the deformation measures are expressed as
@ðq0 wÞ @ðq0 wÞ 1 @ðq0 wÞ 1 þ C WAB þ C CABC @C KL @ WKB LA @ CKBC LA
ð18Þ
RKL ¼ 2
@ðq0 wÞ @ðq0 wÞ 1 @ðq0 wÞ 1 þ 2sym C LA WAB þ 2sym C LA CABC @C KL @ WKB @ CKBC
ð19Þ
MKLM ¼
@ ðq0 wÞ @ CLMK
ð20Þ
SKL ¼ 2
Inserting the Helmholtz free energy function into (18)–(20) and also assuming elastic strains are small (but rotations could be large), we ignore quadratic terms in (18) and (19), leading to simplified stress constitutive equations for SKL and RKL, as well as the equation for MKLM (Regueiro, 2010; Suhubi & Eringen, 1964) as
1330
V. Isbuga, R.A. Regueiro / International Journal of Engineering Science 49 (2011) 1326–1336
SKL ¼ ðk þ sÞEMM dKL þ 2ðl þ rÞEKL þ gE MM dKL þ jE KL þ mE LK
ð21Þ
RKL ¼ ðk þ 2sÞEMM dKL þ 2ðl þ 2rÞEKL þ ð2g sÞE MM dKL þ ðm þ j rÞðE KL þ E LK Þ
ð22Þ
MKLM ¼ s1 ðCKRR dLM þ CRRL dKM Þ þ s2 ðCRKR dLM þ CRRM dKL Þ þ s3 CRRK dLM þ s4 CLRR dKM þ s5 ðCRLR dKM þ CMRR dKL Þ þ s6 CRMR dKL þ s7 CLMK þ s8 ðCMKL þ CKLM Þ þ s9 CLKM þ s10 CMLK þ s11 CKML
ð23Þ
where SKL is the unsymmetric second Piola–Kirchhoff stress tensor, RKL is the symmetric second Piola–Kirchhoff micro-stress tensor, and MKLM is the higher order couple stress tensor in B0 . Note that the stress tensors SKL, RKL, and MKLM are mapped to the current configuration by (8). We can see that if all micromorphic isotropic elastic parameters are zero except k and l, then k⁄ = k and l⁄ = l, where k⁄ and l⁄ are the Lamé parameters. Otherwise, k and l are just two of the eighteen isotropic micromorphic linear elastic parameters. 5.1. Choice of elastic parameters In the Cauchy continuum, the constitutive equations involve two Lamé parameters for linear isotropic elasticity, k⁄ and l⁄. The isotropic micromorphic elasticity approach introduces five more elastic moduli in the unsymmetric Cauchy stress tensor r, as well as in microstress tensor s, and eleven length-scale elastic constants in the higher order couple stress tensor m. Smith (1968) proposed the restrictions among those material parameters in the form of inequalities to achieve the positive definiteness of a quadratic strain energy function. Smith (1968) used the same strain tensors as Suhubi and Eringen (1964) in constructing the free energy function for small strains. This allows us to indirectly use these inequalities since our model is also based on the same strain measures applied by Suhubi and Eringen (1964), but in the reference configuration for finite deformation assuming linear elasticity. In this section, we present the results obtained by Smith (1968). A linearized form of Eulerian strain measures reduces to the small strain tensors
ekl ¼
1 ðuk;l þ ul;k Þ; 2
ekl ¼ /kl þ ul;k ; cklm ¼ /kl;m
ð24Þ
where uk is the displacement vector, and /kl is the micro-displacement tensor in B. We can see the similarity to a linearization of the Lagrangian strain measures in (16) and (11)–(13). Smith (1968) defined sets of variables consisting of the components of the strain measures given in (24), and took the second partial derivative of a strain energy function with respect to those components to obtain four different coefficient matrices each of which must have positive eigenvalues to satisfy the positive definiteness of the strain energy function. After some algebra, the positive eigenvalue condition of those matrices yielded these sets of inequalities among the elastic parameters that must hold (Smith, 1968) as
l > 0; j þ m > 2r; j m > 0; 3k þ 2l > 0; ðj þ m 2rÞl > 2r2 ; j þ m þ 3g > 3s þ 2r ðj þ m þ 3g 3s 2rÞð3k þ 2lÞ > ð3s þ 2rÞ2
ð25Þ
The following set of inequalities specify the constraints on the elastic parameters si as
1
1=2
s7 þ 2s8 > js9 þ s10 þ s11 j; s7 s8 > pffiffiffi ðs9 s10 Þ2 þ ðs10 s11 Þ2 þ ðs11 s9 Þ2
ð26Þ
2
tr ðTÞ > 0;
tr ðCOTÞ > 0;
det ðTÞ > 0
where tr is the trace operator, COT is the cofactor of T where T is defined as
2
3
s1 þ s2 þ 3s3 þ s7 þ s10 3s1 þ s4 þ s5 þ s8 þ s11 3s2 þ s5 þ s6 þ s8 þ s9 6 7 T ¼ 4 3s1 þ s2 þ s3 þ s8 þ s11 s1 þ 3s4 þ s5 þ s7 þ s9 s2 þ 3s5 þ s6 þ s8 þ s10 5 s1 þ 3s2 þ s3 þ s8 þ s9 s1 þ s4 þ 3s5 þ s8 þ s10 s2 þ s5 þ 3s6 þ s7 þ s11 Note that the inequalities in (25) are different than those given in Eringen (1999, p. 278), except the constraints on si, because of the different strain measures and constitutive model for rij and sij used by Eringen (1999) as compared to Suhubi and Eringen (1964), and what we use here. Although we have restrictions on micromorphic elastic material parameters by these inequalities based on the positive definiteness of the strain energy function, determination of the numerical values of these parameters is necessary to apply to physical examples. There are some recent works on determination of micromorphic elastic parameters. Chen and Lee (2003) proposed a method to determine the micromorphic elastic parameters for crystals by using atomic dispersion curves. Vernerey, Liu, and Moran (2007) applied a volume averaging method starting from the micro-scale to determine the material constants in upper-scales for hierarchical materials. Neff and Forest (2007) related the micro-scale parameters to the macroscale Lamé parameters in conjunction with a homogenization approach. They multiplied the macro-scale parameters with a factor to ensure larger values for the micro-scale parameters and applied homogenization to obtain elastic parameters appearing in the constitutive equations. The reasoning behind this approach is based on the physical interpretation of behavior of small samples approaching the microstructural representative volume element (RVE) that are stiffer than the behavior
1331
V. Isbuga, R.A. Regueiro / International Journal of Engineering Science 49 (2011) 1326–1336
of the macro-scale samples. They related the additional elastic parameters associated with the volumetric strain to kRVE, and the additional elastic parameters associated with the shear strain components to lRVE. Zervos, Papanicolopulos, and Vardoulakis (2009) followed a similar approach to determine the additional elastic parameters appearing in the constitutive equations which are based on the gradient elasticity approach. In our model, we follow a similar approach because of its simplicity to relate the micromorphic elastic parameters associated with the volumetric strain to Lamé parameter k⁄, and the micromorphic elastic parameters associated with the full strain tensors to Lamé parameter l⁄, together with considering the inequalities proposed by Smith (1968). The relations between the micromorphic elastic parameters and Lamé parameters k⁄ and l⁄ will be given with the numerical examples in Section 7. 6. Weak form and finite element equations In this section, we summarize the weak form and finite element formulation which is used to implement the model for three dimensional finite element problems. The strong form of the balance equations were given in (4) and (5), for the balance of linear momentum and the balance of first moment of momentum, respectively. Applying the method of weighted residuals and integration by parts to the balance equations (Hughes, 1987) with Nanson’s formula and Piola transform, we get the following equations in the reference configuration. First, we multiply (4) in residual form with the weight function wk and apply integration by parts to obtain
Z
wk ½rlk;l þ qðfk ak Þdv ¼
V
Z
A
wk rlk nl da
Z
wk;l rlk þ wk qðfk ak Þ dv ¼ 0
If we also apply Piola transforms of the Cauchy stress tensor P lL ¼ J rlk F 1 Lk ; nl da ¼ JF 1 Kl N K dA for area change, we obtain
Z
wk PkK NK dA
A
Z
V
ð27Þ
V
rlk ¼ F lL SLK F kK =J, and Nanson’s formula
½wk;l F lL SLK F kK þ wk q0 ðfk ak ÞdV ¼ 0
ð28Þ
where PkK is the first Piola–Kirchhoff stress tensor. If we define a traction Tk = PkLNL, Eq. (28) can be written as
Z
wk T k dA
A
Z
V
½wk;l F lL SLK F kK þ wk q0 ðfk ak ÞdV ¼ 0
ð29Þ
Application of the Newton–Raphson method requires linearization of the balance equations to construct a consistent tangent, where linearization of the balance of linear momentum may be expressed as
Z
wk T k dA
A
Z V
wk;l F lL SLK F kK þ wk q0 ðfk ak Þ dV þ d
Z
wk T k dA
A
Z
V
wk;l F lL SLK F kK þ wk q0 ðfk ak Þ dV
¼0
ð30Þ
where d() is the increment operator within a linearization procedure. Similarly, we follow the same method for the balance of first moment of momentum. Multiplying the residual form of (5) with the weight function gml and applying integration by parts together with Piola transforms and area change relations expressed above, we obtain
Z V
gml ½F mM SML F lL F mM RML F lL þ q0 ðklm xlm ÞdV
Z V
gml;k F kK F lL MKLM vmM dV þ
Z A
gml Mlm dA ¼ 0
ð31Þ
where Mlm is the traction couple tensor, defined as Mlm ¼ Jmklm F 1 Kk N K . Then, linearization of (31) is expressed as
Z
V
gml ½F mM SML F lL F mM RML F lL þ q0 ðklm xlm ÞdV
Z
gml;k F kK F lL MKLM vmM dV þ
Z
gml Mlm dA V A
Z Z Z þd gml ½F mM SML F lL F mM RML F lL þ q0 ðklm xlm ÞdV gml;k F kK F lL MKLM vmM dV þ gml Mlm dA V
¼0
V
A
ð32Þ
Eqs. (30) and (32) give the coupled equations to solve using a Newton–Raphson method after finite element interpolation functions are introduced. For elastostatics, we ignore the acceleration a and micro-spin inertia x, as well as the body force f, body couple k, traction T, and traction couple M. We use a mixed 27 node hexahedral element for finite element implementation (Fig. 2(a)): 27 nodes for uh and 8 vertex nodes for Uh (Q27P8), where h is the discretization parameter when expressing the linearized weak form in Galerkin form (Hughes, 1987); this element is used because its 2D analog has been shown to be convergent for the small strain case, Hughes (1987, p. 215). The mixed formulation is used in recognition that Uh is a micro-displacement tensor (even though no direct gradient is calculated on a micro-displacement vector u0 ), where the gradient of a quadratic interpolation is linear, such that $uh and Uh would be of approximate same order interpolation. Such mixed methods are shown to be convergent (Hughes, 1987) for the small strain case, but no formal proof of convergence is presented here; i.e., we do not show that u = limh?0 uh and U = limh?0 Uh. The global nodal displacement vector is d and the nodal micro-displacement tensor in vector form is /. Displacement vector uhk is interpolated at each node shown by solid dots, and the micro-displacement tensor UhkK is interpolated only at the vertices shown by open circles in Fig. 2(a).
1332
V. Isbuga, R.A. Regueiro / International Journal of Engineering Science 49 (2011) 1326–1336
(b) 5
(a)
4.5 4 3.5
MPa
X3
−σ 33
u3 = −0.5 m
X2
3 2.5 2
*
*
Finite strain elasticity (FSE) with λ ,μ Micromorphic FSE with λ,μ,η, BC1 Micromorphic FSE with λ,μ,η,κ, and ν, BC1 Micromorphic FSE with λ,μ,η,κ,ν, and τ7, BC1
1.5
Micromorphic FSE with λ,μ,η,κ,ν,τ,σ, and τ7, BC1
1
Micromorphic FSE with λ,μ,η, BC2 Micromorphic FSE with λ,μ,η,κ,and ν, BC2 Micromorphic FSE with λ,μ,η,κ,ν, and τ7, BC2
0.5
Micromorphic FSE with λ,μ,η,κ,ν,τ,σ, and τ , BC2 7
0 0
X1
0.05
0.1
0.15
0.2
−e33
0.25
0.3
0.35
0.4
Fig. 2. (a)2 2 2 m cube compression example with a one element mesh. (b) Comparison of Cauchy stress tensor component r33 obtained by standard finite strain elasticity (FSE) and finite strain micromorphic elasticity with various combinations of micromorphic elastic material parameters as well as different BCs.
In the finite element implementation, the terms in (30) and (32) that are not linearized will yield the residual vectors Rd and R/, respectively. The terms involving du and dU in (30) result in the associated consistent tangent matrices Kdd and Kd/, respectively. Similarly, in (32) the associated consistent tangent matrices are, respectively, K/d and K//. We have the following coupled finite element system of equations to solve for nodal displacement increments dd and nodal micro-displacement tensor increments d/ at each iteration in a Newton–Raphson algorithm:
K dd K /d
K d/ K //
dd d/
¼
Rd
ð33Þ
R/
The next section focuses on numerical examples to demonstrate the implementation of the three dimensional finite deformation micromorphic isotropic linear elastic model. 7. Numerical examples The first example is a one element 2 m cube uniaxial stress in compression analysis, allowing one free micro-displacement tensor degree of freedom Uh33 with various combinations of micromorphic isotropic elastic material parameters. As mentioned in Section 5.1, the numerical values may be chosen by a set of relations between Lamé parameters k⁄ and l⁄ and the other micromorphic elastic parameters. The following set of relations satisfy the positive definiteness of strain energy for specific values of k⁄ = 39 103 kPa and l⁄ = 12 103 kPa for a geomaterial as:
s1 s2 ¼ s3 s4 s5 s6
0:111k L2c
0:185k L2c 0:204k L2c 0:1k L2c 0:256k L2c
s7 s8 ¼ s9 s10 s11
0:670l L2c 2 Lc 2 Lc 2 Lc
0:495l 0:408l
0:495l
k
0:7435k
s 0:256k g 1:53k
l m r j
0:583l 0:667l 0:4167l
ð34Þ
0:833l
where Lc is an elastic length scale parameter. Numerical values for the micromorphic elastic material parameters are then calculated as
k ¼ 29 103 kPa;
l ¼ 7 103 kPa; g ¼ 60 103 kPa m ¼ 8 103 kPa; j ¼ 10 103 kPa; s ¼ 10 103 kPa r ¼ 5 103 kPa
ð35Þ
Note that to simplify the computations from the set of relations which satisfy the positive definiteness of strain energy expressed in (34), we take
s1 ¼ s2 ¼ s3 ¼ s4 ¼ s5 ¼ s6 ¼ s8 ¼ s9 ¼ s10 ¼ s11 ¼ 0; s7 ¼ 8:103 kPa m2
ð36Þ
where the characteristic length is assumed to be Lc = 1 m. Fig. 2(a) shows the compression example. The displacement boundary condition u3 = 0.5 m is applied to the top surface in the negative X3 direction. The displacement boundary conditions were chosen to constrain the rigid body motions as:
1333
V. Isbuga, R.A. Regueiro / International Journal of Engineering Science 49 (2011) 1326–1336
u1 = 0 on X1 face, u2 = 0 on X2 face, and u3 = 0 on X3 face. We considered two different sets of boundary conditions for Uh: (1) all the micro-displacement tensor components UhiI are set = 0 except the component Uh33 in the X3 direction as boundary condition type 1 (BC1); (2) Uh33 is free on the top surface of the cube (X3 = 2 m), and Uh33 ¼ 0 at X3 = 0 for boundary condition type 2 (BC2). The different BCs will show the effect of the gradient of Uh33 in the X3 direction. Another example in Fig. 3 illustrates a boundary layer effect. Fig. 2(b) is the comparison of Cauchy stress tensor component r33 obtained by standard finite strain elasticity (FSE) with the same stress component in micromorphic finite strain elasticity calculated by Eqs. (8)1 and (21). Fig. 2(b) shows that the micromorphic approach may result in lower stress for BC1 as compared to BC2. In Fig. 2(b), for BC1, all the Uh33 are the same which, in turn, gives zero gradient of micro-displacement tensor $U ¼ 0 between the top and the bottom surfaces of the cube. For BC2, there is a gradient of Uh33 between the top and bottom surfaces of the cube element that generates a non-zero higher order stress tensor m which results in different relative stress s r, and different Cauchy stress r in the first moment of momentum balance due to the coupling between the two balance equations. Fig. 2(b) shows an increase in r33 by including additional parameters g, j, m, s, and r where for BC1, no contribution from the parameter s7 is expected due to zero gradient of U. The effect of boundary conditions on the micro-displacement tensor U may be presented more clearly in the next example in which a 10 m column with a 1 1 m2 cross section is loaded on the +X3 surface (top surface) with a prescribed displacement, u3 = 1 m. Similar boundary conditions to the one element cube example (BC1 and BC2) are employed, except displacements in X1 and X2 directions are zero. Fig. 3(a) shows the boundary effect on micro-displacement tensor where we have slightly lower values of Uh33 near the bottom surface that is causing a gradient in micro-displacement tensor which results in different values of the higher order stress tensor m at different heights of the column as shown in Fig. 3(c) and (d). We see higher values for invariants of the higher order stress tensor at locations where we have a higher gradient of Uh33 near X3 = 0 (see Fig. 3(b)) which is expected because of the definition of the higher order stress tensor in the reference configuration, for this specific example, M = s7C where C ¼ F T $v. We do not see a similar behavior for BC1 in which $U ¼ 0, that makes the higher order stress tensor zero; therefore, for BC1, stress invariants given in Fig. 3(c)–(f) are the same at different heights along the column (i.e., the stress is uniform along the height of the column for BC1).
(c) 0.9
(e) z=0.1125 m, BC2 z=1.1125 m, BC2 z=9.8875 m, BC2 z=0.1125 m, BC1 z=1.1125 m, BC1 z=9.8875 m, BC1
0.8
||dev(m)|| MPa
0.7 0.6
z=9.8875 m, BC2 z=9.8875 m, BC1 z=1.1125 m, BC1 z=0.1125 m, BC1 z=1.1125 m, BC2 z=0.1125 m, BC2
0.8 0.7 0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1 0 0
0.1 0.1
0.2
0.3
0.4
0.5
time
0.6
0.7
0.8
0.9
0 0
1
(d) 0.9
(f) z=0.1125 m, BC2 z=1.1125 m, BC2 z=9.8875 m, BC2 z=0.1125 m, BC1 z=1.1125 m, BC1 z=9.8875 m, BC1
0.8
||tr(m)|| MPa
0.7 0.6
0.1
0.5 0.4
0.4
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0.5
0.6
0.7
0.8
0.9
1
time
z=9.8875 m, BC2 z=9.8875 m, BC1 z=1.1125 m, BC1 z=0.1125 m, BC1 z=1.1125 m, BC2 z=0.1125 m, BC2
0.6
0.5
0.2
0.7
||dev(s − σ )|| MPa
(b)
1 0.9
||dev(σ )|| MPa
(a)
0.3
0.3
0.2
0.2
0.1
0.1 0 0
0.1
0.2
0.3
0.4
0.5
time
0.6
0.7
0.8
0.9
1
0 0
0.1
0.2
0.3
0.4
time
Fig. 3. (a) Variation of norm of micro-displacement tensor, kUk, through the column height that shows the boundary effect on U, (b) norm of gradient of micro-displacement tensor, k$Uk, that generates non-zero m, (c) norm of deviatoric higher order stress kdev (m)k, which is largest where the highest gradient is observed, (d) norm of trace of m, ktr (m)k, (e) norm of deviatoric Cauchy stress, kdev (r)k, (f) norm of deviatoric relative stress, kdev (s r)k.
1334
V. Isbuga, R.A. Regueiro / International Journal of Engineering Science 49 (2011) 1326–1336
The third example is a three dimensional square corner punch problem to demonstrate the three dimensionality of the finite element implementation. It is a 10 m cube with different li lim2 square areas loaded with downward prescribed displacement u3 = 1 m (see Fig. 4). Note that all simulations were run to large strains where the assumption that elastic strains are small in (21) and (22) (but rotations can still be large) would be invalid. This implementation, however, will eventually be used as a precursor to a finite strain elasto-plastic model implementation, where plastic strains can be large for such materials like soil, rock, concrete, etc., where elastic strains are typically small. We assumed that only three micro-displacement tensor components Uh11 ; Uh22 , and Uh33 are free, and all the other shear components UhiI ¼ 0 ði – IÞ. The boundary conditions on micro-displacement tensor are chosen as: Uh11 ¼ 0 on ±X1 faces, Uh22 ¼ 0 on ±X2 faces, and Uh33 ¼ 0 on X3 face according to Fig. 4. Elastic parameters given in (35) and (36) were used in the analysis. Fig. 5 shows the variation over loading of the norms of deviatoric stress kdev rk, kdev (s r)k, and kdev mk and traces tr r, tr (s r), and ktr mk for the unsymmetric Cauchy stress tensor r, the relative stress tensor (s r), and the higher order couple stress tensor m, respectively, at the Gauss point near the node A underneath the punch area. Definitions of these invariant stress measures are
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðdev rij Þðdev rij Þ tr r ¼ rkk pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kdev ðs rÞk ¼ ½dev ðsij rij Þ½dev ðsij rij Þ tr ðs rÞ ¼ ðskk rkk Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kdev mk ¼ ðdev mijk Þðdev mijk Þ ktr mk ¼ maak mbbk kdev rk ¼
ð37Þ
where
dev rij ¼ rij
1 1 1 rkk dij ; dev sij rij ¼ ðsij rij Þ ðskk rkk Þ dij ; dev mijk ¼ mijk dij maak 3 3 3
ð38Þ
Fig. 5(a) and (b) show the deviatoric stress norms for various combinations of elastic parameters k, l, g, j, m, s, r, s7. Upon including terms with s and r, kdev rk and k dev (s r)kdecrease until the end of loading, when they increase. The stress state and deformation is three dimensional, so unlike the previous simple examples of uniaxial compression, the combination of deformation and micromorphic elastic parameters can lead to unexpected trends in stress. Upon introducing s7 to include the effect of higher order stress tensor m, the relative deviatoric stress norm kdev (s r)k tends to decrease, whereas kdev rkincreases. For k dev (s r)k, this behavior is also observed in the column compression example: kdev (s r)k is smaller where higher gradient of Uh33 is seen that causes higher values of higher order stress tensor. Fig. 5(c) and (f) agree with results shown in Fig. 3(c) and (d) regarding the effect of the gradient of micro-displacement tensor $U which leads to an increase in magnitude of higher order stress tensor m that in turn induces a change in relative and Cauchy stress tensors, s r and r, respectively. Fig. 5(d)–(f) show the first invariants of r, (s r), and m. Fig. 5(d) shows decreased first invariant (increasing compressive stress, which is negative) for tr r, while tr (s r) in Fig. 5(e) goes into tension when s and r are introduced. Because of the definition ktr mk, it is always positive, a measure that must be changed to reflect the sign of tr m(tension or compression) when used in a pressure-sensitive plasticity model (Regueiro, 2009). Recall that onlys7 is nonzero in (36), so that the higher order stress is mklm = s7FkKFlLCLMKvmM/J. The stress measures kdev rk,kdev (s r)k, kdev mk, tr r, tr (s r), and ktr mk used in this work will provide a basis for a future implementation of the three dimensional micromorphic finite strain elasto-plastic pressure sensitive model in which three yield functions are defined for macro-scale, micro-scale, and micro-scale gradient plasticity concept (Regueiro, 2009, 2010). Fig. 5(a) and (d) also include the comparison of kdev rk and tr r for the standard finite strain symmetric Cauchy stress tensor with those of micromorphic finite strain unsymmetric Cauchy stress tensor with the various combinations of micromorphic isotropic elastic parameters. The micromorphic results, depending on the deformation history, show either stiffer or softer response than standard elasticity. These
u3 = −1 m
A x
10 m
li
li X3 X2 X1 10 m 10 m Fig. 4. (left) 10 10 10 m cube with various square punch areas, si (top right), showing Gauss point as X near nodal point A where stresses are plotted in Fig. 5. Convergence profile obtained by Newton–Raphson algorithm at the first time step for the largest punch area (bottom right). The time step is Dt = 0.025 and total time = 1. There are 10 10 10 = 1000 mixed Q27P8 hexahedral elements in the mesh.
1335
V. Isbuga, R.A. Regueiro / International Journal of Engineering Science 49 (2011) 1326–1336
1.4
*
1
7
0.8
0.4
0
*
−10
0.25
0.4
time
0.6
0.8
−14 0
1
(e) Micromorphic FSE with Micromorphic FSE with Micromorphic FSE with Micromorphic FSE with
λ, μ, κ,ν λ, μ, η,κ,ν λ, μ, η,κ,ν,τ,σ λ, μ, η,κ,ν,τ,σ and τ7
(MPa)
0.3
0.2
0.15 0.1
2
0.2
0.4
time
0.6
0.8
0.2
0.3
Micromorphic FSE with Micromorphic FSE with Micromorphic FSE with Micromorphic FSE with
0.4
0.5
time
0.6
0.7
0.8
0.9
1
λ, μ, κ,ν λ, μ, η,κ,ν λ, μ, η,κ,ν,τ,σ λ, μ, η,κ,ν,τ,σ and τ
7
1.5
(f)
Micromorphic FSE with λ, μ, η,κ,ν,τ,σ,τ7, punch area s5
0.5
7
0.1
0.2
0.3
0.4
0.5
time
0.6
0.7
0.8
0.9
1
Micromorphic FSE with λ, μ, η,κ,ν,τ,σ,τ7, punch area s5 Micromorphic FSE with λ, μ, η,κ,ν,τ,σ,τ , punch area s
1
Micromorphic FSE with λ, μ, η,κ,ν,τ,σ,τ7, punch area s2 Micromorphic FSE with λ, μ, η,κ,ν,τ,σ,τ , punch area s 7
2.5
2
Micromorphic FSE with λ, μ, η,κ,ν,τ,σ,τ , punch area s
(MPa)
1.5
1
−0.5 0
1
Micromorphic FSE with λ, μ, η,κ,ν,τ,σ,τ7, punch area s4
3
7
4
7
1
7
2
Micromorphic FSE with λ, μ, η,κ,ν,τ,σ,τ , punch area s Micromorphic FSE with λ, μ, η,κ,ν,τ,σ,τ , punch area s
1.5
Micromorphic FSE with λ, μ, η,κ,ν,τ,σ,τ7, punch area s3
1
||trm||
(MPa)
2
0.1
0
0.05 0 0
||devm||
2.5
*
Finite strain elasticity (FSE) with λ , μ Micromorphic FSE with λ, μ, κ,ν Micromorphic FSE with λ, μ, η,κ,ν Micromorphic FSE with λ, μ, η,κ,ν,τ,σ Micromorphic FSE with λ, μ, η,κ,ν,τ,σ and τ7
0.2
tr (s − σ)
(MPa)
(b) 0.35
||dev (s − σ) ||
−6
−12
0.2
(c)
−4
−8
0.6
0
0 −2
trσ
(MPa)
1.2
||devσ ||
(d)
*
Finite strain elasticity (FSE) with λ , μ Micromorphic FSE with λ, μ, κ,ν Micromorphic FSE with λ, μ, η,κ,ν Micromorphic FSE with λ, μ, η,κ,ν,τ,σ Micromorphic FSE with λ, μ, η,κ,ν,τ,σ and τ
(MPa)
(a)
0.5
0 0
1
0.5
0.1
0.2
0.3
0.4
0.5
time
0.6
0.7
0.8
0.9
1
0 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
time
Fig. 5. Deviatoric stress norms (a) kdev rk, (b) kdev (s r)k, and (c) kdev mk, and traces (d) tr r, (e) tr (s r), and norm of trace of m (f) k trmk at the Gauss point closest to the node A under the punch area. Plots (c) and (f) show the results under five different punch areas.
results signify the importance of the selection of micromorphic elastic parameters in addition to micro-displacement tensor Uh boundary conditions. The table in Fig. 4 illustrates the convergence profile for the large deformation punch problem, with parameters k, l, g, j, m, s, r and s7. 8. Conclusions We have implemented a three dimensional finite element model for finite deformation micromorphic materially linear isotropic elasticity into an open source finite element code Tahoe. The paper presented preliminary results of the model applied to three three-dimensional examples. These results motivate further study of the meaning of certain micromorphic isotropic elastic parameters and boundary conditions on micro-displacement tensor U, as well as their influence on the simulations. In the numerical examples, micromorphic elastic parameters have been chosen to satisfy the positive definiteness of strain energy that obey the inequalities proposed by Smith (1968). Numerical values have been chosen in a similar way followed by Neff and Forest (2007) and Zervos et al. (2009). The boundary conditions on U may not be explicitly anticipated from the physical problem as discussed by Eringen (1968). In addition, we noticed from several analyses that U boundary conditions together with the associated micromorphic elastic parameters play an important role in convergence of the nonlinear Newton–Raphson solution algorithm. Even for the comparatively simple example of uniaxial compression, specification of the boundary conditions on U and the choice of micromorphic elastic parameters may affect the numerical results. A square corner punch problem demonstrated the full three dimensionality of the finite element implementation and some interesting trends in stress response at a Gauss point beneath the punch with various combinations of micromorphic elastic parameters, and the length scale effect for various punch areas.
1336
V. Isbuga, R.A. Regueiro / International Journal of Engineering Science 49 (2011) 1326–1336
Acknowledgements The support of National Science Foundation grant CMMI-0700648, Army Research Office Grant W911NF-09-1-0111, and the Army Research Laboratory is gratefully acknowledged. The authors would like to recognize and dedicate this paper to Dr. A.C. Eringen, who even in his later years continued to inspire young researchers in materials mechanics and continuum physics. References Chambon, R., Caillerie, D., & Tamagnini, C. (2004). A strain space gradient plasticity theory for finite strain. Computer Methods in Applied Mechanics and Engineering, 193(27–29), 2797–2826. Chen, Y., & Lee, J. D. (2003). Determining material constants in micromorphic theory through phonon dispersion relations. International Journal of Engineering Science, 41(8), 871–886. Coleman, B. D., & Gurtin, M. E. (1967). Thermodynamics with internal state variables. Journal of Chemical Physics, 47, 597–613. Coleman, B. D., & Noll, W. (1963). The thermodynamics of elastic materials with heat conduction and viscosity. Archive for Rational Mechanics and Analysis, 13, 167–178. Kröner, E. (Ed.). (1968). Mechanics of Generalized Continua. Springer-Verlag. Eringen, A. (1962). Nonlinear theory of continuous media (1st ed.). McGraw-Hill. Eringen, A. (1968). Theory of micropolar elasticity. In H. Liebowitz (Ed.), Fracture. An Advanced Treatise (Vol. 2, pp. 622–729). Academic Press. Eringen, A. (1999). Microcontinuum field theories. I: Foundations and solids. Springer-Verlag. Eringen, A. C., & Suhubi, E. S. (1964). Nonlinear theory of simple micro-elastic solids – 1. International Journal of Engineering Science, 2(2), 189–203. Forest, S., & Sievert, R. (2006). Nonlinear microstrain theories. International Journal of Solids and Structures, 43(24), 7224–7245. Germain, P. (1973). The method of virtual power in continuum mechanics. Part 2: Microstructure. SIAM Journal on Applied Mathematics, 25(3), 556–575. Green, A. E., & Naghdi, P. M. (1995). A unified procedure for construction of theories of deformable media. II: Generalized continua. Proceedings of the Royal Society of London, Series A: Mathematical and Physical Sciences, 448(1934), 357–377. Holzapfel, G. A. (2000). Nonlinear solid mechanics: A continuum approach for engineering. John Wiley & Sons. Hughes, T. J. R. (1987). The finite element method. New Jersey: Prentice-Hall. Mindlin, R. (1964). Micro-structure in linear elasticity. Archive for Rational Mechanics and Analysis, 16, 51–78. Neff, P., & Forest, S. (2007). A geometrically exact micromorphic model for elastic metallic foams accounting for affine microstructure. modelling, existence of minimizers, identification of moduli and computational results. Journal of Elasticity, 87(2), 239–276. Regueiro, R. (2009). Finite strain micromorphic pressure-sensitive plasticity. Journal of Engineering Mechanics, 135, 178–191. Regueiro, R. (2010). On finite strain micromorphic elastoplasticity. International Journal of Solids and Structures, 47, 786–800. Sansour, C., Skatulla, S., & Zbib, H. (2010). A formulation for the micromorphic continuum at finite inelastic strains. International Journal of Solids and Structures, 47(11–12), 1546–1554. Smith, A. C. (1968). Inequalities between the constants of a linear micro-elastic solid. International Journal of Engineering Science, 6(2), 65–74. Suhubi, E., & Eringen, A. (1964). Nonlinear theory of simple micro-elastic solids – II. International Journal of Engineering Science, 2(2), 389–404. Vernerey, F., Liu, W., & Moran, B. (2007). Multi-scale micromorphic theory for hierarchical materials. Journal of the Mechanics and Physics of Solids, 55(12), 2603–2651. Zervos, A., Papanicolopulos, S.-A., & Vardoulakis, I. (2009). Two finite-element discretizations for gradient elasticity. Journal of Engineering Mechanics, 135(3), 203–213.
Contents lists available at ScienceDirect
International Journal of Engineering Science journal homepage: www.elsevier.com/locate/ijengsci
Three-dimensional finite element analysis of finite deformation micromorphic linear isotropic elasticity Volkan Isbuga, Richard A. Regueiro ⇑ Department of Civil, Environmental, and Architectural Engineering, University of Colorado at Boulder, United States
a r t i c l e
i n f o
Article history: Received 11 October 2010 Accepted 19 April 2011 Available online 12 May 2011 Keywords: Micromorphic elasticity Finite deformation Finite element implementation Three dimensional
a b s t r a c t A finite deformation micromorphic materially linear isotropic elastic model is formulated and implemented for three dimensional finite element analysis. The model is based on the kinematics, balance equations and thermodynamic equations proposed by Eringen and Suhubi (1964). The constitutive equations are calculated in the reference configuration, and the resulting stresses are mapped to the current configuration. The balance of linear momentum and the balance of first moment of momentum are linearized to construct the consistent tangent for three dimensional finite element implementation for solution by the Newton–Raphson method. Three dimensional numerical examples are analyzed to demonstrate preliminarily the implementation. Ó 2011 Elsevier Ltd. All rights reserved.
1. Introduction The nonlinear theory of a micromorphic continuum has been first proposed by Eringen and Suhubi (1964) and Suhubi and Eringen (1964). The micromorphic continuum in which each material point is endowed with nine additional kinematic degrees of freedom (dof) represents the most general case of higher order continua of ‘‘grade one’’ (Eringen, 1999) (i.e., condef sidering only the first gradient of micro-structure deformation v ¼ @n=@N) (Eringen & Suhubi, 1964; Green & Naghdi, 1995; Kröner, 1968). These additional dof allow modeling material response when classical continuum theories may not be sufficient, such as when size effects of the micro-structure become important under large deformation gradient (Chambon, Caillerie, & Tamagnini, 2004). For example, in particulate materials, strain localization has a size effect related to the grain-size, and thus it is desirable to capture this size effect using multiscale models. Most of the micromorphic continuum based models follow two main approaches for the kinematic assumptions and balance equations which were proposed by Mindlin (1964) and Eringen and Suhubi (1964). Germain (1973) also provided an approach in determining the balance equations of generalized continua by implementing the principle of virtual work. Neff and Forest (2007) applied a geometrically exact micromorphic model with a homogenization approach for elastic metallic foams. Forest and Sievert (2006) provided some guidelines based on previous literature to choose the proper higher order continuum or higher grade continuum in modeling length-scale-dependent material response. Sansour, Skatulla, and Zbib (2010) developed a micromorphic model for finite strain inelasticity, and they applied the model to metals. Numerous works on micromorphic continuum based modeling of material response with elastic and plastic models are still on-going, and it is not possible to mention all such works here. However, a summary and development of the micro-continuum theories may be found in Eringen (1999). In this work, we implemented a finite deformation micromorphic materially linear isotropic elastic model to simulate three dimensional problems using the finite element method, which establishes the basis for future elastoplastic model ⇑ Corresponding author. Fax: +1 303 492 7317. E-mail address: [email protected] (R.A. Regueiro). 0020-7225/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijengsci.2011.04.006
V. Isbuga, R.A. Regueiro / International Journal of Engineering Science 49 (2011) 1326–1336
1327
implementation (Regueiro, 2009, 2010). The model follows the kinematics and balance equations proposed by Eringen and Suhubi (1964), Suhubi and Eringen (1964) and Eringen (1999). One of the main reasons that we follow this approach is that all macro-scale stress tensors and higher order couple stress tensor together with body force vectors and body force couple tensors are expressed in terms of micro-scale tensors and parameters (Eringen, 1999). This provides an opportunity to use the model in a multiscale modeling approach where a bridging/coupling between scales may be needed. Also, the nonlinear continuum mechanics of the formulation are developed rigorously by Eringen and co-workers. We assume Catersian coordinates for continuum finite element implementation. Index notation will be used throughout so as to be as clear as possible with regard to details of the formulation. Some symbolic/direct notation is also given, such as (ab)ik = aijbjk. Boldface denotes a tensor or vector, where its index notation has been given uniformly throughout the paper. Generally, variables in uppercase letters live in the reference configuration B0 (such as the reference differential volume dV), and variables in lowercase live in the current configuration B (such as the current differential volume dv). The same applies to their indices, such that a differential line segment in the current configuration dxi is related to a differential line segment in the reference configuration dXI through the deformation gradient: dxi = FiIdXI (Einstein’s summation convention assumed (see Eringen, 1962; Holzapfel, 2000)). Subscripts (),i and (),I imply partial differentiation in the current and reference configurations, respectively. A superscript prime symbol ()0 denotes a variable associated with the micro-element (see Section 2 _ ¼ DðÞ=Dt denotes material time derivative. on kinematics). Superposed dot ðÞ An outline of the remainder of the paper is as follows: Section 2 summarizes the kinematics of a micromorphic continuum; Section 3 presents the momentum balance equations as well as Clausius–Duhem inequality in the current configuration; and Section 4 shows the mapping of Clausius–Duhem inequality to reference configuration in which constitutive equations are derived. Constitutive equations are given in Section 5. Weak form and finite element equations are given in Section 6. Choice of elastic parameters and numerical results are presented in Section 7, and conclusions in Section 8. 2. Micromorphic kinematics Fig. 1 illustrates the mapping of the macro-element and micro-element in the reference configuration B0 to the current configuration B through the deformation gradient F and micro-deformation tensor v. The macro-element continuum point is denoted by P(X, N) and p(x, n, t) in the reference and current configurations, respectively, with centroid C and c. The microelement continuum point centroid is denoted by C0 and c0 in the reference and current configurations, respectively. The micro-element in general represents a grain/particle/fiber micro-structural sub-volume (cluster) of the heterogeneous material. The relative position vector of the micro-element centroid with respect to the macro-element centroid is denoted by N and n(X, N, t) in the reference and current configurations, respectively, such that the micro-element centroid position vectors are written as (Fig. 1)
X 0K ¼ X K þ NK ;
x0k ¼ xk ðX; tÞ þ nk ðX; N; tÞ
ð1Þ
Eringen and Suhubi (1964) assumed that for sufficiently small lengths kNk 1 (kk is the L2 norm), n is linearly related to N through the micro-deformation tensor v, such that
nk ðX; N; tÞ ¼ vkK ðX; tÞNK
ð2Þ
where then the spatial position vector of the micro-element centroid is written as
x0k ¼ xk ðX; tÞ þ vkK ðX; tÞNK
Fig. 1. Deformation of macro differential volume dV and micro differential volume dV0 (Eringen, 1968).
ð3Þ
1328
V. Isbuga, R.A. Regueiro / International Journal of Engineering Science 49 (2011) 1326–1336
This is called a micromorphic continuum theory of ‘‘grade one’’ (Eringen, 1999). This is equivalent to assuming an affine, or homogeneous, deformation of the macro-element differential volume dV (but not necessarily the body B; i.e., the continuum body B may experience heterogeneous deformation because of v, even if boundary conditions (BCs) are uniform). It also simplifies considerably the formulation of the micromorphic continuum balance equations as presented in Eringen and Suhubi (1964) and Eringen (1999). This micro-deformation tensor v is analogous to the small strain ‘‘micro-deformation’’ tensor w in Mindlin (1964), physically described in Mindlin’s Fig. 1. Eringen (1968) also provides a physical interpretation of v generally, but then simplifies for the micropolar case. For example, v can be interpreted as calculated from a micro-displacement tensor U as v = 1 + U, where U is not calculated from a micro-displacement vector u0 , because in general U, and in turn v, are incompatible. For finite element implementation, U will be interpolated at the nodes, providing an additional nine dofs because it is unsymmetric. 3. Micromorphic balance equations and Clausius–Duhem inequality In this section, the balance of linear momentum, first moment of momentum, and Clausius–Duhem inequality are presented in the current configuration B. The derivation of the equations and the detailed integral definitions of all the vectors and the stress tensors as well as the higher order stress and couple tensors can be found in Eringen and Suhubi (1964) and Eringen (1999), and therefore will not be given here. Eq. (4) is the balance of linear momentum in a classical Cauchy continuum theory, however, the Cauchy stress tensor ris now unsymmetric. Eq. (5) was first introduced by Eringen and Suhubi (1964), which is called the balance of first moment of momentum that is more general than the balance of angular momentum. It contains six additional equations, for a total of nine equations for nine unknowns U. Eq. (6) is the Clausius–Duhem inequality in the current configuration B. The equations are written as
rlk;l þ qðfk ak Þ ¼ 0
ð4Þ
rml sml þ mklm;k þ qðklm xlm Þ ¼ 0
ð5Þ
1 q w_ þ gh_ þ qk h;k þ rkl v l;k þ ðskl rkl Þmlk þ mklm mlm;k P 0 h
ð6Þ
where f and a are the body force vector per unit mass and the acceleration vector respectively, s the micro-stress tensor, m the higher order couple stress tensor, k the body force couple tensor per unit mass, x the micro-spin inertia tensor per unit _ 1 the micro-gyration tensor, $v the velocity gradient, q the heat flux vector per unit mass, w the Helmholtz free mass, m ¼ vv energy function per unit mass, q the mass density in the current configuration, and h the temperature. In the derivation of constitutive equations, we map these equations into the reference configuration in which constitutive equations are more naturally introduced, such as to incorporate anisotropy in a finite strain setting and/or to introduce plasticity in a naturally elastically-relaxed intermediate configuration (Regueiro, 2010). 4. Clausius–Duhem inequality in B0 The Clausius–Duhem inequality given in (6) expressed in current configuration B can be rewritten in the reference configuration B0 by using the relations q = q0/J and dv = J dV, as well as Nanson’s formula and the Piola transform in local form as
1 q0 w_ þ gh_ þ J rkl ðv l;k mlk Þ þ Jskl mlk þ J mlm;k mklm þ Q K h;K P 0 h
ð7Þ
where the Piola transforms for the unsymmetric Cauchy stress tensor rkl, symmetric micro-stress tensor skl, and higher order couple stress tensor mklm are written as
1 J
1 J
1 J
rkl ¼ F kK SKL F lL ; skl ¼ F kK RKL F lL ; mklm ¼ F kK F lL MKLM vmM
ð8Þ
Inserting (8) into (7) gives us the Clausius–Duhem inequality in the reference configuration as
1 q0 w_ þ gh_ þ F kK SKL F lL ðv l;k mlk Þ þ F kK RKL F lL mlk þ mlm;k F kK F lL M KLM vmM þ Q K h;K P 0 h
ð9Þ
The Clausius–Duhem inequality in conjunction with the Helmholtz free energy function, which is presented in the next section, is applied to derive the elastic constitutive equations in the reference configuration B0 . The rest of the paper focuses on the Helmholtz free energy function assumption, the derivation of constitutive equations, weak form, finite element formulation and linearization, and numerical examples.
V. Isbuga, R.A. Regueiro / International Journal of Engineering Science 49 (2011) 1326–1336
1329
5. Helmholtz free energy function and constitutive equations Eringen and Suhubi (1964) proposed unique sets of elastic deformation measures which are invariant with respect to rigid body motion in the current configuration. For the Helmholtz free energy function per unit mass in B0 , we use the set consisting of C, W, and C defined as
C KL ¼ F iK F iL ;
WKL ¼ F iK viL ;
CKLM ¼ F iK viL;M
ð10Þ
These deformation tensors can be written in terms of the displacement vector U and micro-displacement tensor U as follows (Eringen, 1968; Suhubi & Eringen, 1964)
C KL ¼ dKL þ U K;L þ U L;K þ U M;K U M;L
ð11Þ
WKL ¼ dKL þ UKL þ U L;K þ U M;K UML
ð12Þ
CKLM ¼ UKL;M þ U N;K UNL;M
ð13Þ
Then, the Helmholtz free energy density in terms of these deformation measures is expressed as
q0 wðC; W; CÞ; q0 wðC KL ; WKL ; CKLM Þ
ð14Þ
A quadratic form assumption for the Helmholtz free energy function in the reference configuration will lead to materially linear elasticity at finite strain as
1 2
1 2
1 2
q0 w ¼ EKL AKLMN EMN þ E KL BKLMN E MN þ CKLM C KLMNPQ CNPQ þ EKL DKLMN E MN
ð15Þ
where elastic strain tensors are defined in Suhubi and Eringen (1964) as
EKL ¼
1 ðC KL dKL Þ; 2
E KL ¼ WKL dKL
ð16Þ
AKLMN, BKLMN, DKLMN, and CKLMNPQ are the elastic material moduli tensors (Regueiro, 2010) defined for isotropic linear elasticity and obtained from back analysis of linear constitutive equations and energy function given in Suhubi and Eringen (1964). They are written as
AKLMN ¼ kdKL dMN þ lðdKM dLN þ dKN dLM Þ BKLMN ¼ ðg sÞdKL dMN þ ðj rÞdKM dLN þ ðm rÞdKN dLM C KLMNPQ ¼ s1 ðdKL dMN dPQ þ dKQ dLM dNP Þ þ s2 ðdKL dMP dNQ þ dKM dLQ dNP Þ þ s3 dKL dMQ dNP þ s4 dKN dLM dPQ þ s5 ðdKM dLN dPQ þ dKP dLM dNQ Þ
ð17Þ
þ s6 dKM dLP dNQ þ s7 dKN dLP dMQ þ s8 ðdKP dLQ dMN þ dKQ dLN dMP Þ þ s9 dKN dLQ dMP þ s10 dKP dLN dMQ þ s11 dKQ dLP dMN DKLMN ¼ sdKL dMN þ rðdKN dLM þ dLN dKM Þ where major and minor symmetry hold for AKLMN and DKLMN, whereas only major symmetry holds for the tensors BKLMN and CKLMNPQ, k, l, g, s, j, m, r, s1, . . . , s11 are the eighteen micromorphic isotropic elastic parameters. Note that the units for s1. . .s11 are stress length2 (Pa m2), thus there is a built in length scale to these elastic parameters for the higher order stress. These elastic material moduli tensors are different than those of elastic material moduli tensors given in Eringen (1999) due to different strain measures in our formulation, except the higher order elastic modulus tensor CKLMNPQwhich has the same strain measures for the higher order couple stress tensor in Eringen (1999). Following the classical argument of Coleman and Noll (1963) and Coleman and Gurtin (1967) for independent rate pro_ Wand _ _ the constitutive equations in terms of partial derivatives of Helmholtz free energy function with recesses for C; C, spect to the deformation measures are expressed as
@ðq0 wÞ @ðq0 wÞ 1 @ðq0 wÞ 1 þ C WAB þ C CABC @C KL @ WKB LA @ CKBC LA
ð18Þ
RKL ¼ 2
@ðq0 wÞ @ðq0 wÞ 1 @ðq0 wÞ 1 þ 2sym C LA WAB þ 2sym C LA CABC @C KL @ WKB @ CKBC
ð19Þ
MKLM ¼
@ ðq0 wÞ @ CLMK
ð20Þ
SKL ¼ 2
Inserting the Helmholtz free energy function into (18)–(20) and also assuming elastic strains are small (but rotations could be large), we ignore quadratic terms in (18) and (19), leading to simplified stress constitutive equations for SKL and RKL, as well as the equation for MKLM (Regueiro, 2010; Suhubi & Eringen, 1964) as
1330
V. Isbuga, R.A. Regueiro / International Journal of Engineering Science 49 (2011) 1326–1336
SKL ¼ ðk þ sÞEMM dKL þ 2ðl þ rÞEKL þ gE MM dKL þ jE KL þ mE LK
ð21Þ
RKL ¼ ðk þ 2sÞEMM dKL þ 2ðl þ 2rÞEKL þ ð2g sÞE MM dKL þ ðm þ j rÞðE KL þ E LK Þ
ð22Þ
MKLM ¼ s1 ðCKRR dLM þ CRRL dKM Þ þ s2 ðCRKR dLM þ CRRM dKL Þ þ s3 CRRK dLM þ s4 CLRR dKM þ s5 ðCRLR dKM þ CMRR dKL Þ þ s6 CRMR dKL þ s7 CLMK þ s8 ðCMKL þ CKLM Þ þ s9 CLKM þ s10 CMLK þ s11 CKML
ð23Þ
where SKL is the unsymmetric second Piola–Kirchhoff stress tensor, RKL is the symmetric second Piola–Kirchhoff micro-stress tensor, and MKLM is the higher order couple stress tensor in B0 . Note that the stress tensors SKL, RKL, and MKLM are mapped to the current configuration by (8). We can see that if all micromorphic isotropic elastic parameters are zero except k and l, then k⁄ = k and l⁄ = l, where k⁄ and l⁄ are the Lamé parameters. Otherwise, k and l are just two of the eighteen isotropic micromorphic linear elastic parameters. 5.1. Choice of elastic parameters In the Cauchy continuum, the constitutive equations involve two Lamé parameters for linear isotropic elasticity, k⁄ and l⁄. The isotropic micromorphic elasticity approach introduces five more elastic moduli in the unsymmetric Cauchy stress tensor r, as well as in microstress tensor s, and eleven length-scale elastic constants in the higher order couple stress tensor m. Smith (1968) proposed the restrictions among those material parameters in the form of inequalities to achieve the positive definiteness of a quadratic strain energy function. Smith (1968) used the same strain tensors as Suhubi and Eringen (1964) in constructing the free energy function for small strains. This allows us to indirectly use these inequalities since our model is also based on the same strain measures applied by Suhubi and Eringen (1964), but in the reference configuration for finite deformation assuming linear elasticity. In this section, we present the results obtained by Smith (1968). A linearized form of Eulerian strain measures reduces to the small strain tensors
ekl ¼
1 ðuk;l þ ul;k Þ; 2
ekl ¼ /kl þ ul;k ; cklm ¼ /kl;m
ð24Þ
where uk is the displacement vector, and /kl is the micro-displacement tensor in B. We can see the similarity to a linearization of the Lagrangian strain measures in (16) and (11)–(13). Smith (1968) defined sets of variables consisting of the components of the strain measures given in (24), and took the second partial derivative of a strain energy function with respect to those components to obtain four different coefficient matrices each of which must have positive eigenvalues to satisfy the positive definiteness of the strain energy function. After some algebra, the positive eigenvalue condition of those matrices yielded these sets of inequalities among the elastic parameters that must hold (Smith, 1968) as
l > 0; j þ m > 2r; j m > 0; 3k þ 2l > 0; ðj þ m 2rÞl > 2r2 ; j þ m þ 3g > 3s þ 2r ðj þ m þ 3g 3s 2rÞð3k þ 2lÞ > ð3s þ 2rÞ2
ð25Þ
The following set of inequalities specify the constraints on the elastic parameters si as
1
1=2
s7 þ 2s8 > js9 þ s10 þ s11 j; s7 s8 > pffiffiffi ðs9 s10 Þ2 þ ðs10 s11 Þ2 þ ðs11 s9 Þ2
ð26Þ
2
tr ðTÞ > 0;
tr ðCOTÞ > 0;
det ðTÞ > 0
where tr is the trace operator, COT is the cofactor of T where T is defined as
2
3
s1 þ s2 þ 3s3 þ s7 þ s10 3s1 þ s4 þ s5 þ s8 þ s11 3s2 þ s5 þ s6 þ s8 þ s9 6 7 T ¼ 4 3s1 þ s2 þ s3 þ s8 þ s11 s1 þ 3s4 þ s5 þ s7 þ s9 s2 þ 3s5 þ s6 þ s8 þ s10 5 s1 þ 3s2 þ s3 þ s8 þ s9 s1 þ s4 þ 3s5 þ s8 þ s10 s2 þ s5 þ 3s6 þ s7 þ s11 Note that the inequalities in (25) are different than those given in Eringen (1999, p. 278), except the constraints on si, because of the different strain measures and constitutive model for rij and sij used by Eringen (1999) as compared to Suhubi and Eringen (1964), and what we use here. Although we have restrictions on micromorphic elastic material parameters by these inequalities based on the positive definiteness of the strain energy function, determination of the numerical values of these parameters is necessary to apply to physical examples. There are some recent works on determination of micromorphic elastic parameters. Chen and Lee (2003) proposed a method to determine the micromorphic elastic parameters for crystals by using atomic dispersion curves. Vernerey, Liu, and Moran (2007) applied a volume averaging method starting from the micro-scale to determine the material constants in upper-scales for hierarchical materials. Neff and Forest (2007) related the micro-scale parameters to the macroscale Lamé parameters in conjunction with a homogenization approach. They multiplied the macro-scale parameters with a factor to ensure larger values for the micro-scale parameters and applied homogenization to obtain elastic parameters appearing in the constitutive equations. The reasoning behind this approach is based on the physical interpretation of behavior of small samples approaching the microstructural representative volume element (RVE) that are stiffer than the behavior
1331
V. Isbuga, R.A. Regueiro / International Journal of Engineering Science 49 (2011) 1326–1336
of the macro-scale samples. They related the additional elastic parameters associated with the volumetric strain to kRVE, and the additional elastic parameters associated with the shear strain components to lRVE. Zervos, Papanicolopulos, and Vardoulakis (2009) followed a similar approach to determine the additional elastic parameters appearing in the constitutive equations which are based on the gradient elasticity approach. In our model, we follow a similar approach because of its simplicity to relate the micromorphic elastic parameters associated with the volumetric strain to Lamé parameter k⁄, and the micromorphic elastic parameters associated with the full strain tensors to Lamé parameter l⁄, together with considering the inequalities proposed by Smith (1968). The relations between the micromorphic elastic parameters and Lamé parameters k⁄ and l⁄ will be given with the numerical examples in Section 7. 6. Weak form and finite element equations In this section, we summarize the weak form and finite element formulation which is used to implement the model for three dimensional finite element problems. The strong form of the balance equations were given in (4) and (5), for the balance of linear momentum and the balance of first moment of momentum, respectively. Applying the method of weighted residuals and integration by parts to the balance equations (Hughes, 1987) with Nanson’s formula and Piola transform, we get the following equations in the reference configuration. First, we multiply (4) in residual form with the weight function wk and apply integration by parts to obtain
Z
wk ½rlk;l þ qðfk ak Þdv ¼
V
Z
A
wk rlk nl da
Z
wk;l rlk þ wk qðfk ak Þ dv ¼ 0
If we also apply Piola transforms of the Cauchy stress tensor P lL ¼ J rlk F 1 Lk ; nl da ¼ JF 1 Kl N K dA for area change, we obtain
Z
wk PkK NK dA
A
Z
V
ð27Þ
V
rlk ¼ F lL SLK F kK =J, and Nanson’s formula
½wk;l F lL SLK F kK þ wk q0 ðfk ak ÞdV ¼ 0
ð28Þ
where PkK is the first Piola–Kirchhoff stress tensor. If we define a traction Tk = PkLNL, Eq. (28) can be written as
Z
wk T k dA
A
Z
V
½wk;l F lL SLK F kK þ wk q0 ðfk ak ÞdV ¼ 0
ð29Þ
Application of the Newton–Raphson method requires linearization of the balance equations to construct a consistent tangent, where linearization of the balance of linear momentum may be expressed as
Z
wk T k dA
A
Z V
wk;l F lL SLK F kK þ wk q0 ðfk ak Þ dV þ d
Z
wk T k dA
A
Z
V
wk;l F lL SLK F kK þ wk q0 ðfk ak Þ dV
¼0
ð30Þ
where d() is the increment operator within a linearization procedure. Similarly, we follow the same method for the balance of first moment of momentum. Multiplying the residual form of (5) with the weight function gml and applying integration by parts together with Piola transforms and area change relations expressed above, we obtain
Z V
gml ½F mM SML F lL F mM RML F lL þ q0 ðklm xlm ÞdV
Z V
gml;k F kK F lL MKLM vmM dV þ
Z A
gml Mlm dA ¼ 0
ð31Þ
where Mlm is the traction couple tensor, defined as Mlm ¼ Jmklm F 1 Kk N K . Then, linearization of (31) is expressed as
Z
V
gml ½F mM SML F lL F mM RML F lL þ q0 ðklm xlm ÞdV
Z
gml;k F kK F lL MKLM vmM dV þ
Z
gml Mlm dA V A
Z Z Z þd gml ½F mM SML F lL F mM RML F lL þ q0 ðklm xlm ÞdV gml;k F kK F lL MKLM vmM dV þ gml Mlm dA V
¼0
V
A
ð32Þ
Eqs. (30) and (32) give the coupled equations to solve using a Newton–Raphson method after finite element interpolation functions are introduced. For elastostatics, we ignore the acceleration a and micro-spin inertia x, as well as the body force f, body couple k, traction T, and traction couple M. We use a mixed 27 node hexahedral element for finite element implementation (Fig. 2(a)): 27 nodes for uh and 8 vertex nodes for Uh (Q27P8), where h is the discretization parameter when expressing the linearized weak form in Galerkin form (Hughes, 1987); this element is used because its 2D analog has been shown to be convergent for the small strain case, Hughes (1987, p. 215). The mixed formulation is used in recognition that Uh is a micro-displacement tensor (even though no direct gradient is calculated on a micro-displacement vector u0 ), where the gradient of a quadratic interpolation is linear, such that $uh and Uh would be of approximate same order interpolation. Such mixed methods are shown to be convergent (Hughes, 1987) for the small strain case, but no formal proof of convergence is presented here; i.e., we do not show that u = limh?0 uh and U = limh?0 Uh. The global nodal displacement vector is d and the nodal micro-displacement tensor in vector form is /. Displacement vector uhk is interpolated at each node shown by solid dots, and the micro-displacement tensor UhkK is interpolated only at the vertices shown by open circles in Fig. 2(a).
1332
V. Isbuga, R.A. Regueiro / International Journal of Engineering Science 49 (2011) 1326–1336
(b) 5
(a)
4.5 4 3.5
MPa
X3
−σ 33
u3 = −0.5 m
X2
3 2.5 2
*
*
Finite strain elasticity (FSE) with λ ,μ Micromorphic FSE with λ,μ,η, BC1 Micromorphic FSE with λ,μ,η,κ, and ν, BC1 Micromorphic FSE with λ,μ,η,κ,ν, and τ7, BC1
1.5
Micromorphic FSE with λ,μ,η,κ,ν,τ,σ, and τ7, BC1
1
Micromorphic FSE with λ,μ,η, BC2 Micromorphic FSE with λ,μ,η,κ,and ν, BC2 Micromorphic FSE with λ,μ,η,κ,ν, and τ7, BC2
0.5
Micromorphic FSE with λ,μ,η,κ,ν,τ,σ, and τ , BC2 7
0 0
X1
0.05
0.1
0.15
0.2
−e33
0.25
0.3
0.35
0.4
Fig. 2. (a)2 2 2 m cube compression example with a one element mesh. (b) Comparison of Cauchy stress tensor component r33 obtained by standard finite strain elasticity (FSE) and finite strain micromorphic elasticity with various combinations of micromorphic elastic material parameters as well as different BCs.
In the finite element implementation, the terms in (30) and (32) that are not linearized will yield the residual vectors Rd and R/, respectively. The terms involving du and dU in (30) result in the associated consistent tangent matrices Kdd and Kd/, respectively. Similarly, in (32) the associated consistent tangent matrices are, respectively, K/d and K//. We have the following coupled finite element system of equations to solve for nodal displacement increments dd and nodal micro-displacement tensor increments d/ at each iteration in a Newton–Raphson algorithm:
K dd K /d
K d/ K //
dd d/
¼
Rd
ð33Þ
R/
The next section focuses on numerical examples to demonstrate the implementation of the three dimensional finite deformation micromorphic isotropic linear elastic model. 7. Numerical examples The first example is a one element 2 m cube uniaxial stress in compression analysis, allowing one free micro-displacement tensor degree of freedom Uh33 with various combinations of micromorphic isotropic elastic material parameters. As mentioned in Section 5.1, the numerical values may be chosen by a set of relations between Lamé parameters k⁄ and l⁄ and the other micromorphic elastic parameters. The following set of relations satisfy the positive definiteness of strain energy for specific values of k⁄ = 39 103 kPa and l⁄ = 12 103 kPa for a geomaterial as:
s1 s2 ¼ s3 s4 s5 s6
0:111k L2c
0:185k L2c 0:204k L2c 0:1k L2c 0:256k L2c
s7 s8 ¼ s9 s10 s11
0:670l L2c 2 Lc 2 Lc 2 Lc
0:495l 0:408l
0:495l
k
0:7435k
s 0:256k g 1:53k
l m r j
0:583l 0:667l 0:4167l
ð34Þ
0:833l
where Lc is an elastic length scale parameter. Numerical values for the micromorphic elastic material parameters are then calculated as
k ¼ 29 103 kPa;
l ¼ 7 103 kPa; g ¼ 60 103 kPa m ¼ 8 103 kPa; j ¼ 10 103 kPa; s ¼ 10 103 kPa r ¼ 5 103 kPa
ð35Þ
Note that to simplify the computations from the set of relations which satisfy the positive definiteness of strain energy expressed in (34), we take
s1 ¼ s2 ¼ s3 ¼ s4 ¼ s5 ¼ s6 ¼ s8 ¼ s9 ¼ s10 ¼ s11 ¼ 0; s7 ¼ 8:103 kPa m2
ð36Þ
where the characteristic length is assumed to be Lc = 1 m. Fig. 2(a) shows the compression example. The displacement boundary condition u3 = 0.5 m is applied to the top surface in the negative X3 direction. The displacement boundary conditions were chosen to constrain the rigid body motions as:
1333
V. Isbuga, R.A. Regueiro / International Journal of Engineering Science 49 (2011) 1326–1336
u1 = 0 on X1 face, u2 = 0 on X2 face, and u3 = 0 on X3 face. We considered two different sets of boundary conditions for Uh: (1) all the micro-displacement tensor components UhiI are set = 0 except the component Uh33 in the X3 direction as boundary condition type 1 (BC1); (2) Uh33 is free on the top surface of the cube (X3 = 2 m), and Uh33 ¼ 0 at X3 = 0 for boundary condition type 2 (BC2). The different BCs will show the effect of the gradient of Uh33 in the X3 direction. Another example in Fig. 3 illustrates a boundary layer effect. Fig. 2(b) is the comparison of Cauchy stress tensor component r33 obtained by standard finite strain elasticity (FSE) with the same stress component in micromorphic finite strain elasticity calculated by Eqs. (8)1 and (21). Fig. 2(b) shows that the micromorphic approach may result in lower stress for BC1 as compared to BC2. In Fig. 2(b), for BC1, all the Uh33 are the same which, in turn, gives zero gradient of micro-displacement tensor $U ¼ 0 between the top and the bottom surfaces of the cube. For BC2, there is a gradient of Uh33 between the top and bottom surfaces of the cube element that generates a non-zero higher order stress tensor m which results in different relative stress s r, and different Cauchy stress r in the first moment of momentum balance due to the coupling between the two balance equations. Fig. 2(b) shows an increase in r33 by including additional parameters g, j, m, s, and r where for BC1, no contribution from the parameter s7 is expected due to zero gradient of U. The effect of boundary conditions on the micro-displacement tensor U may be presented more clearly in the next example in which a 10 m column with a 1 1 m2 cross section is loaded on the +X3 surface (top surface) with a prescribed displacement, u3 = 1 m. Similar boundary conditions to the one element cube example (BC1 and BC2) are employed, except displacements in X1 and X2 directions are zero. Fig. 3(a) shows the boundary effect on micro-displacement tensor where we have slightly lower values of Uh33 near the bottom surface that is causing a gradient in micro-displacement tensor which results in different values of the higher order stress tensor m at different heights of the column as shown in Fig. 3(c) and (d). We see higher values for invariants of the higher order stress tensor at locations where we have a higher gradient of Uh33 near X3 = 0 (see Fig. 3(b)) which is expected because of the definition of the higher order stress tensor in the reference configuration, for this specific example, M = s7C where C ¼ F T $v. We do not see a similar behavior for BC1 in which $U ¼ 0, that makes the higher order stress tensor zero; therefore, for BC1, stress invariants given in Fig. 3(c)–(f) are the same at different heights along the column (i.e., the stress is uniform along the height of the column for BC1).
(c) 0.9
(e) z=0.1125 m, BC2 z=1.1125 m, BC2 z=9.8875 m, BC2 z=0.1125 m, BC1 z=1.1125 m, BC1 z=9.8875 m, BC1
0.8
||dev(m)|| MPa
0.7 0.6
z=9.8875 m, BC2 z=9.8875 m, BC1 z=1.1125 m, BC1 z=0.1125 m, BC1 z=1.1125 m, BC2 z=0.1125 m, BC2
0.8 0.7 0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1 0 0
0.1 0.1
0.2
0.3
0.4
0.5
time
0.6
0.7
0.8
0.9
0 0
1
(d) 0.9
(f) z=0.1125 m, BC2 z=1.1125 m, BC2 z=9.8875 m, BC2 z=0.1125 m, BC1 z=1.1125 m, BC1 z=9.8875 m, BC1
0.8
||tr(m)|| MPa
0.7 0.6
0.1
0.5 0.4
0.4
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0.5
0.6
0.7
0.8
0.9
1
time
z=9.8875 m, BC2 z=9.8875 m, BC1 z=1.1125 m, BC1 z=0.1125 m, BC1 z=1.1125 m, BC2 z=0.1125 m, BC2
0.6
0.5
0.2
0.7
||dev(s − σ )|| MPa
(b)
1 0.9
||dev(σ )|| MPa
(a)
0.3
0.3
0.2
0.2
0.1
0.1 0 0
0.1
0.2
0.3
0.4
0.5
time
0.6
0.7
0.8
0.9
1
0 0
0.1
0.2
0.3
0.4
time
Fig. 3. (a) Variation of norm of micro-displacement tensor, kUk, through the column height that shows the boundary effect on U, (b) norm of gradient of micro-displacement tensor, k$Uk, that generates non-zero m, (c) norm of deviatoric higher order stress kdev (m)k, which is largest where the highest gradient is observed, (d) norm of trace of m, ktr (m)k, (e) norm of deviatoric Cauchy stress, kdev (r)k, (f) norm of deviatoric relative stress, kdev (s r)k.
1334
V. Isbuga, R.A. Regueiro / International Journal of Engineering Science 49 (2011) 1326–1336
The third example is a three dimensional square corner punch problem to demonstrate the three dimensionality of the finite element implementation. It is a 10 m cube with different li lim2 square areas loaded with downward prescribed displacement u3 = 1 m (see Fig. 4). Note that all simulations were run to large strains where the assumption that elastic strains are small in (21) and (22) (but rotations can still be large) would be invalid. This implementation, however, will eventually be used as a precursor to a finite strain elasto-plastic model implementation, where plastic strains can be large for such materials like soil, rock, concrete, etc., where elastic strains are typically small. We assumed that only three micro-displacement tensor components Uh11 ; Uh22 , and Uh33 are free, and all the other shear components UhiI ¼ 0 ði – IÞ. The boundary conditions on micro-displacement tensor are chosen as: Uh11 ¼ 0 on ±X1 faces, Uh22 ¼ 0 on ±X2 faces, and Uh33 ¼ 0 on X3 face according to Fig. 4. Elastic parameters given in (35) and (36) were used in the analysis. Fig. 5 shows the variation over loading of the norms of deviatoric stress kdev rk, kdev (s r)k, and kdev mk and traces tr r, tr (s r), and ktr mk for the unsymmetric Cauchy stress tensor r, the relative stress tensor (s r), and the higher order couple stress tensor m, respectively, at the Gauss point near the node A underneath the punch area. Definitions of these invariant stress measures are
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðdev rij Þðdev rij Þ tr r ¼ rkk pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kdev ðs rÞk ¼ ½dev ðsij rij Þ½dev ðsij rij Þ tr ðs rÞ ¼ ðskk rkk Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kdev mk ¼ ðdev mijk Þðdev mijk Þ ktr mk ¼ maak mbbk kdev rk ¼
ð37Þ
where
dev rij ¼ rij
1 1 1 rkk dij ; dev sij rij ¼ ðsij rij Þ ðskk rkk Þ dij ; dev mijk ¼ mijk dij maak 3 3 3
ð38Þ
Fig. 5(a) and (b) show the deviatoric stress norms for various combinations of elastic parameters k, l, g, j, m, s, r, s7. Upon including terms with s and r, kdev rk and k dev (s r)kdecrease until the end of loading, when they increase. The stress state and deformation is three dimensional, so unlike the previous simple examples of uniaxial compression, the combination of deformation and micromorphic elastic parameters can lead to unexpected trends in stress. Upon introducing s7 to include the effect of higher order stress tensor m, the relative deviatoric stress norm kdev (s r)k tends to decrease, whereas kdev rkincreases. For k dev (s r)k, this behavior is also observed in the column compression example: kdev (s r)k is smaller where higher gradient of Uh33 is seen that causes higher values of higher order stress tensor. Fig. 5(c) and (f) agree with results shown in Fig. 3(c) and (d) regarding the effect of the gradient of micro-displacement tensor $U which leads to an increase in magnitude of higher order stress tensor m that in turn induces a change in relative and Cauchy stress tensors, s r and r, respectively. Fig. 5(d)–(f) show the first invariants of r, (s r), and m. Fig. 5(d) shows decreased first invariant (increasing compressive stress, which is negative) for tr r, while tr (s r) in Fig. 5(e) goes into tension when s and r are introduced. Because of the definition ktr mk, it is always positive, a measure that must be changed to reflect the sign of tr m(tension or compression) when used in a pressure-sensitive plasticity model (Regueiro, 2009). Recall that onlys7 is nonzero in (36), so that the higher order stress is mklm = s7FkKFlLCLMKvmM/J. The stress measures kdev rk,kdev (s r)k, kdev mk, tr r, tr (s r), and ktr mk used in this work will provide a basis for a future implementation of the three dimensional micromorphic finite strain elasto-plastic pressure sensitive model in which three yield functions are defined for macro-scale, micro-scale, and micro-scale gradient plasticity concept (Regueiro, 2009, 2010). Fig. 5(a) and (d) also include the comparison of kdev rk and tr r for the standard finite strain symmetric Cauchy stress tensor with those of micromorphic finite strain unsymmetric Cauchy stress tensor with the various combinations of micromorphic isotropic elastic parameters. The micromorphic results, depending on the deformation history, show either stiffer or softer response than standard elasticity. These
u3 = −1 m
A x
10 m
li
li X3 X2 X1 10 m 10 m Fig. 4. (left) 10 10 10 m cube with various square punch areas, si (top right), showing Gauss point as X near nodal point A where stresses are plotted in Fig. 5. Convergence profile obtained by Newton–Raphson algorithm at the first time step for the largest punch area (bottom right). The time step is Dt = 0.025 and total time = 1. There are 10 10 10 = 1000 mixed Q27P8 hexahedral elements in the mesh.
1335
V. Isbuga, R.A. Regueiro / International Journal of Engineering Science 49 (2011) 1326–1336
1.4
*
1
7
0.8
0.4
0
*
−10
0.25
0.4
time
0.6
0.8
−14 0
1
(e) Micromorphic FSE with Micromorphic FSE with Micromorphic FSE with Micromorphic FSE with
λ, μ, κ,ν λ, μ, η,κ,ν λ, μ, η,κ,ν,τ,σ λ, μ, η,κ,ν,τ,σ and τ7
(MPa)
0.3
0.2
0.15 0.1
2
0.2
0.4
time
0.6
0.8
0.2
0.3
Micromorphic FSE with Micromorphic FSE with Micromorphic FSE with Micromorphic FSE with
0.4
0.5
time
0.6
0.7
0.8
0.9
1
λ, μ, κ,ν λ, μ, η,κ,ν λ, μ, η,κ,ν,τ,σ λ, μ, η,κ,ν,τ,σ and τ
7
1.5
(f)
Micromorphic FSE with λ, μ, η,κ,ν,τ,σ,τ7, punch area s5
0.5
7
0.1
0.2
0.3
0.4
0.5
time
0.6
0.7
0.8
0.9
1
Micromorphic FSE with λ, μ, η,κ,ν,τ,σ,τ7, punch area s5 Micromorphic FSE with λ, μ, η,κ,ν,τ,σ,τ , punch area s
1
Micromorphic FSE with λ, μ, η,κ,ν,τ,σ,τ7, punch area s2 Micromorphic FSE with λ, μ, η,κ,ν,τ,σ,τ , punch area s 7
2.5
2
Micromorphic FSE with λ, μ, η,κ,ν,τ,σ,τ , punch area s
(MPa)
1.5
1
−0.5 0
1
Micromorphic FSE with λ, μ, η,κ,ν,τ,σ,τ7, punch area s4
3
7
4
7
1
7
2
Micromorphic FSE with λ, μ, η,κ,ν,τ,σ,τ , punch area s Micromorphic FSE with λ, μ, η,κ,ν,τ,σ,τ , punch area s
1.5
Micromorphic FSE with λ, μ, η,κ,ν,τ,σ,τ7, punch area s3
1
||trm||
(MPa)
2
0.1
0
0.05 0 0
||devm||
2.5
*
Finite strain elasticity (FSE) with λ , μ Micromorphic FSE with λ, μ, κ,ν Micromorphic FSE with λ, μ, η,κ,ν Micromorphic FSE with λ, μ, η,κ,ν,τ,σ Micromorphic FSE with λ, μ, η,κ,ν,τ,σ and τ7
0.2
tr (s − σ)
(MPa)
(b) 0.35
||dev (s − σ) ||
−6
−12
0.2
(c)
−4
−8
0.6
0
0 −2
trσ
(MPa)
1.2
||devσ ||
(d)
*
Finite strain elasticity (FSE) with λ , μ Micromorphic FSE with λ, μ, κ,ν Micromorphic FSE with λ, μ, η,κ,ν Micromorphic FSE with λ, μ, η,κ,ν,τ,σ Micromorphic FSE with λ, μ, η,κ,ν,τ,σ and τ
(MPa)
(a)
0.5
0 0
1
0.5
0.1
0.2
0.3
0.4
0.5
time
0.6
0.7
0.8
0.9
1
0 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
time
Fig. 5. Deviatoric stress norms (a) kdev rk, (b) kdev (s r)k, and (c) kdev mk, and traces (d) tr r, (e) tr (s r), and norm of trace of m (f) k trmk at the Gauss point closest to the node A under the punch area. Plots (c) and (f) show the results under five different punch areas.
results signify the importance of the selection of micromorphic elastic parameters in addition to micro-displacement tensor Uh boundary conditions. The table in Fig. 4 illustrates the convergence profile for the large deformation punch problem, with parameters k, l, g, j, m, s, r and s7. 8. Conclusions We have implemented a three dimensional finite element model for finite deformation micromorphic materially linear isotropic elasticity into an open source finite element code Tahoe. The paper presented preliminary results of the model applied to three three-dimensional examples. These results motivate further study of the meaning of certain micromorphic isotropic elastic parameters and boundary conditions on micro-displacement tensor U, as well as their influence on the simulations. In the numerical examples, micromorphic elastic parameters have been chosen to satisfy the positive definiteness of strain energy that obey the inequalities proposed by Smith (1968). Numerical values have been chosen in a similar way followed by Neff and Forest (2007) and Zervos et al. (2009). The boundary conditions on U may not be explicitly anticipated from the physical problem as discussed by Eringen (1968). In addition, we noticed from several analyses that U boundary conditions together with the associated micromorphic elastic parameters play an important role in convergence of the nonlinear Newton–Raphson solution algorithm. Even for the comparatively simple example of uniaxial compression, specification of the boundary conditions on U and the choice of micromorphic elastic parameters may affect the numerical results. A square corner punch problem demonstrated the full three dimensionality of the finite element implementation and some interesting trends in stress response at a Gauss point beneath the punch with various combinations of micromorphic elastic parameters, and the length scale effect for various punch areas.
1336
V. Isbuga, R.A. Regueiro / International Journal of Engineering Science 49 (2011) 1326–1336
Acknowledgements The support of National Science Foundation grant CMMI-0700648, Army Research Office Grant W911NF-09-1-0111, and the Army Research Laboratory is gratefully acknowledged. The authors would like to recognize and dedicate this paper to Dr. A.C. Eringen, who even in his later years continued to inspire young researchers in materials mechanics and continuum physics. References Chambon, R., Caillerie, D., & Tamagnini, C. (2004). A strain space gradient plasticity theory for finite strain. Computer Methods in Applied Mechanics and Engineering, 193(27–29), 2797–2826. Chen, Y., & Lee, J. D. (2003). Determining material constants in micromorphic theory through phonon dispersion relations. International Journal of Engineering Science, 41(8), 871–886. Coleman, B. D., & Gurtin, M. E. (1967). Thermodynamics with internal state variables. Journal of Chemical Physics, 47, 597–613. Coleman, B. D., & Noll, W. (1963). The thermodynamics of elastic materials with heat conduction and viscosity. Archive for Rational Mechanics and Analysis, 13, 167–178. Kröner, E. (Ed.). (1968). Mechanics of Generalized Continua. Springer-Verlag. Eringen, A. (1962). Nonlinear theory of continuous media (1st ed.). McGraw-Hill. Eringen, A. (1968). Theory of micropolar elasticity. In H. Liebowitz (Ed.), Fracture. An Advanced Treatise (Vol. 2, pp. 622–729). Academic Press. Eringen, A. (1999). Microcontinuum field theories. I: Foundations and solids. Springer-Verlag. Eringen, A. C., & Suhubi, E. S. (1964). Nonlinear theory of simple micro-elastic solids – 1. International Journal of Engineering Science, 2(2), 189–203. Forest, S., & Sievert, R. (2006). Nonlinear microstrain theories. International Journal of Solids and Structures, 43(24), 7224–7245. Germain, P. (1973). The method of virtual power in continuum mechanics. Part 2: Microstructure. SIAM Journal on Applied Mathematics, 25(3), 556–575. Green, A. E., & Naghdi, P. M. (1995). A unified procedure for construction of theories of deformable media. II: Generalized continua. Proceedings of the Royal Society of London, Series A: Mathematical and Physical Sciences, 448(1934), 357–377. Holzapfel, G. A. (2000). Nonlinear solid mechanics: A continuum approach for engineering. John Wiley & Sons. Hughes, T. J. R. (1987). The finite element method. New Jersey: Prentice-Hall. Mindlin, R. (1964). Micro-structure in linear elasticity. Archive for Rational Mechanics and Analysis, 16, 51–78. Neff, P., & Forest, S. (2007). A geometrically exact micromorphic model for elastic metallic foams accounting for affine microstructure. modelling, existence of minimizers, identification of moduli and computational results. Journal of Elasticity, 87(2), 239–276. Regueiro, R. (2009). Finite strain micromorphic pressure-sensitive plasticity. Journal of Engineering Mechanics, 135, 178–191. Regueiro, R. (2010). On finite strain micromorphic elastoplasticity. International Journal of Solids and Structures, 47, 786–800. Sansour, C., Skatulla, S., & Zbib, H. (2010). A formulation for the micromorphic continuum at finite inelastic strains. International Journal of Solids and Structures, 47(11–12), 1546–1554. Smith, A. C. (1968). Inequalities between the constants of a linear micro-elastic solid. International Journal of Engineering Science, 6(2), 65–74. Suhubi, E., & Eringen, A. (1964). Nonlinear theory of simple micro-elastic solids – II. International Journal of Engineering Science, 2(2), 389–404. Vernerey, F., Liu, W., & Moran, B. (2007). Multi-scale micromorphic theory for hierarchical materials. Journal of the Mechanics and Physics of Solids, 55(12), 2603–2651. Zervos, A., Papanicolopulos, S.-A., & Vardoulakis, I. (2009). Two finite-element discretizations for gradient elasticity. Journal of Engineering Mechanics, 135(3), 203–213.