Effective elastoplastic properties of the periodic composites

SAVE OFFLINE

Computational Materials Science 22 (2001) 221±239

www.elsevier.com/locate/commatsci

E€ective elastoplastic properties of the periodic composites M. Kami nski a,1, è. Figiel b,* b

a Department of Civil Engineering, Rice University, 6100 Main Street, 77005-1892 Houston, TX, USA Department of Structure and Mechanics, Institute of Polymer Research, Hahe Straûe 6, D-01069 Dresden, Germany

Received 5 June 2000; received in revised form 26 March 2001; accepted 24 May 2001

Abstract The article presented deals with the homogenization of composite materials with elastoplastic constituents. The transformation ®eld analysis (TFA) approach is presented and applied to compute the e€ective nonlinear behavior of multicomponent periodic composite structure. Computational implementation of the method consists in special utilization of the program ABAQUS, which makes it possible to homogenize n-component periodic composites with relatively general con®guration of the periodicity cell. Numerical example of homogenization of a three-component periodic composite shows the comparison between the nonlinear behavior of a real composite and of a homogenized one in a speci®c boundary problem de®ned on its representative volume element (RVE). Ó 2001 Elsevier Science B.V. All rights reserved. Keywords: Periodic composites; Elastoplastic e€ective properties; Transformation ®eld analysis; Finite element method

1. Introduction Computational modeling of composite materials is a very complicated process due to the scale di€erences between the entire composite structure, its constituents and interfaces between them and also because of the randomness of reinforcement shape and location, local lacks of periodicity in composite structures as well as the structural defects appearing within and between the macrohomogeneous components [7,20]. One of the most powerful tools to speed up the modelling process, both the composite discretization and the computer simulation of composites in real conditions, is the homogenization method. The main idea of the method is to ®nd a globally homogeneous medium equivalent to the original composite, where the strain energy stored in both systems is approximately the same. It should be mentioned that this method is only an intermediate tool to eciently analyze composite materials; the next step is to solve the original boundary problem using the homogenized composite. The comparison of these two models veri®es the range of usefulness of a speci®c homogenization theory in computational engineering practice. There are numerous well-established techniques to calculate the e€ective materials characteristics for composite materials. In case of composite components volume fractions speci®ed only, one can use the *

Corresponding author. Present address: Division of Mechanics of Materials, Technical University of è odz Al. Politechniki 6, 93-590 è odz, Poland. Tel.: +48-42-631-3551; fax: +48-42-631-3551. E-mail addresses: ®[email protected] (è. Figiel), [email protected], [email protected] (M. Kami nski). URL: http://kmm.p.lodz.pl/Marcin_Kaminski 1

0927-0256/01/$ - see front matter Ó 2001 Elsevier Science B.V. All rights reserved. PII: S 0 9 2 7 - 0 2 5 6 ( 0 1 ) 0 0 1 9 2 - 6

222

M. Kaminski, è. Figiel / Computational Materials Science 22 (2001) 221±239

closed form algebraic equations on upper and lower bounds or direct estimates for the e€ective material tensors components [4,13,14,18]. On the other hand, for the composite with detailed information on the microgeometry, so-called cell problems are formulated and solved [11,15±18] using their ®nite element method (FEM) or, alternatively, the boundary element method (BEM) numerical implementations, which enable direct computations of the e€ective characteristics; the methods mentioned above deal mainly with the homogenization of a composite in elastic range. Recent advances in the area of computational methods in homogenization of nonlinear e€ective characterization of heterogeneous materials and structures are reported in [4,11,17,21]. In the same time, stochastic analysis is still being developed to estimate or to compute probabilistic moments of homogenized material tensors [16]. The main idea behind the article is to present the homogenization of composite materials with elastoplastic constituents. The transformation ®eld analysis (TFA) approach proposed by Dvorak and co-workers, cf. [8±10], is presented brie¯y and applied to compute the e€ective nonlinear behavior of a three-component periodic composite. The self-consistent model [13] and Mori±Tanaka's theory [5] providing estimation of the overall thermoelastic constants of composites on the basis of constituent properties and volume fractions are used. Computational implementation of the method consists in special utilization of the program ABAQUS [1] to enable automatic homogenization of n-component periodic composites in a general con®guration of the components in the periodicity cell. Numerical example of three-component periodic composite homogenization makes it possible to compare the nonlinear behavior of a composite in its real and homogenized con®guration for some speci®c boundary problem de®ned on the representative volume element (RVE). The next step in development of this approach would be to determine the parameter sensitivity of the homogenized properties of the composite with respect to material characteristics of the constituents as well as to geometrical data de®ning the RVE. Statistical and stochastic simulation of probabilistic moments of the e€ective material tensors would be possible after such a sensitivity determination, using the experimental knowledge on statistical parameters of the composite. 2. Homogenization method The periodic n-component composite in the plane orthogonal to the ®ber direction is considered where perfectly bonded components are assumed as elastoplastic [2]. The symmetric elasticity tensor is denoted by C, while the example of some composite as well as its RVE are shown in Figs. 1 and 2.

Fig. 1. The periodic composite domain.

M. Kaminski, è. Figiel / Computational Materials Science 22 (2001) 221±239

223

Fig. 2. The RVE of the composite.

Mechanical behavior of the composite constituents is represented by time and temperature-dependent constitutive relations under the assumption that for any time s the total strains and stresses can be decomposed as er …y; s† ˆ eelr …y; s† ‡ er …y; s†;

…1†

rr …y; s† ˆ relr …y; s† ‡ rr …y; s†;

…2†

where eelr , er denote the elastic strain resulting from a given displacement boundary conditions applied on the region Vr and the eigenstrain in the same subvolume, respectively; relr , rr denote the elastic stress and eigenstress tensor components in Vr . The eigenstrain and eigenstress ®elds considered here as transformation ®eld may be decomposed in case of thermal and inelastic e€ects as er …y; s† ˆ mr h…s† ‡ einel r …y; s†;

…3†

rr …y; s† ˆ lr h…s† ‡ rrel r …y; s†;

…4†

where mr denotes the thermal strain tensor and lr is the thermal stress tensor. Since einel is the inelastic strain r and rrel is the relaxation stress, Eqs. (3) and (4) can be written that r er …y; s† ˆ Mr rr …y; s† ‡ mr h…s† ‡ einel r …y; s†; rr …y; s† ˆ Cr er …y; s† ‡ lr h…s† ‡

rrel r …y; s†;

…5† …6†

where Cr and Mr are the elastic and compliance tensors components for the subregion Vr . Hence, it is possible to write the following relations between mr ; lr ; Mr ; Cr ; einel and rrel r r : Mr ˆ C r 1 ;

…7†

mr ˆ

Mr l r ;

…8†

Cr mr ;

…9†

lr ˆ

einel r …y; s† ˆ

Mr rrel r …y; s†;

…10†

rrel r …y; s† ˆ

Cr einel r …y; s†;

…11†

224

M. Kaminski, è. Figiel / Computational Materials Science 22 (2001) 221±239

while the mechanical and thermal elastic in¯uence functions are given by the following relations: er …y; s† ˆ Ar …y†e…s† ‡ Drs …y; y0 †es …y; s†;

…12†

rr …y; s† ˆ Br …y†e…s† ‡ Frs …y; y0 †rs …y; s†:

…13†

Matrices Br …y† and Ar …y† in Eqs. (12) and (13) denote stress and strain concentration factor tensors representing the volume averages of corresponding functions over the periodicity cell [12], as it is proposed in Eqs. (16)±(19). To describe the overall homogenized response of volume V , the resulting strains and stresses are combined with their corresponding local components described by Eqs. (3)±(6) as Z Z  el  1 1 e…s† ˆ er …y; s† dV ˆ er …y; s† ‡ er …y; s† dV ˆ eel …s† ‡ e …s†; …14† jV j V jV j V Z Z  el  1 1 r…s† ˆ rr …y; s† dV ˆ rr …y; s† ‡ rr …y; s† dV ˆ rel …s† ‡ r …s†: …15† jV j V jV j V The local elastic ®elds related to their counterparts may be written in the form of Z 1 eel …s† ˆ ‰Ar …y†e…s† ‡ ar …y†aT …s†Š dV ; jV j V Z 1 el r …s† ˆ ‰Br …y†r…s† ‡ br …y†aT …s†Š dV ; jV j V

…16† …17†

where ar …y† and br …y† are the thermoelastic strain and stress concentration factors [12]. The strain transformation ®eld e …y; s† de®ned in V results in the displacements on the unconstrained part of surface oV , while transformation stress r …y; s† generates surface tractions on V being constrained. The relation between the local and global transformation ®elds was proposed as [8±10] Z Z 1 1 e …s† ˆ ‰er …y; s† Ar …y†e…s†Š dV ˆ Ar …y†e …y; s† dV ; …18† jV j V jV j V Z Z 1 1  r …s† ˆ ‰rr …y; s† Br …y†r…s†Š dV ˆ BT …y†r …y; s† dV : …19† jV j V jV j V r The elastic local ®elds of strain er …y; s† and stress rr …y; s† are found from superposition of the elastic local ®elds eelr …y; s† and relr …y; s† with local eigenstrains er …y; s† and eigenstresses rr …y; s†, respectively [cf. Eqs. (3) and (4)]; the same model in the context of global scale is introduced analogously to the equations posed above. These two di€erent scales are linked using formulation [8±10] for local strain or/and stress ®elds in the form of er …y0 † ˆ Ar …y0 †e ‡ Drs …y; y0 †es …y0 †; 0

0

rr …y † ˆ Br …y †r ‡ Frs …y; y

0

†rs …y0 †;

…20† …21†

where Drs …y; y0 †, Frs …y; y0 † are transformation strain and stress in¯uence functions, that enable to relate the strain and stress tensor components on the macroscale de®ned by y and the microscale identi®ed by y0 . Solve the following boundary value problem on the RVE div r…y† ˆ

or…y† ˆ 0; oy

er …y† ˆ Mr rr …y† ‡ er …y†;

…22† …23†

M. Kaminski, è. Figiel / Computational Materials Science 22 (2001) 221±239

1 jV j

225

Z V

er …y† dV ˆ heiV ;

…24†

u…y† ˆ ey ‡ u …y†;

…25†

where the local uniform strain ®eld er is found using the matrices Ar …y†, Drs …y; y0 †. Further, it is possible to determine the approximated expression of the averaged strain in the subvolume Vr given as er ˆ A r e ‡

N X r;sˆ1

Drs er :

…26†

Analogously to Eq. (26), the averaged stresses in the subregion Vr can be written in the form of rr ˆ Br r ‡

N X r;sˆ1

Frs rr :

…27†

It is observed that Fr …y; y0 † and Dr …y; y0 † are eigenstress and eigenstrain in¯uence functions, that re¯ect the e€ect on the scale y resulting from a transformation on the scale y0 under overall uniform applied stress or strain. Respective in¯uence functions are determined in terms of averages and piecewise constant approximations in introduced subvolumes of the RVE [9,10]. Therefore, under overall strain e…t† ˆ 0, the transformation concentration factor tensor Drs gives the strain induced in subvolume Vr by a unit uniform eigenstrain located in Vs . Under overall stress r…t† ˆ 0, the concentration factor tensor Frs de®nes the stress in Vr resulting from the unit eigenstrain located in Vs . It can be shown that these tensors can be determined in case of multiphase media as
Dsr ˆ …I Ar †…Cr C† 1 drs I cs ATs Cs ; …28†
…29† Fsr ˆ …I Br †…Mr M† 1 drs I cs BTs Ms …r; s ˆ 1; . . . ; N , without summation over repeated indices) which for two-component composite gives Dpa ˆ …I

Ap †…Ca

Cb † 1 Ca

Fpa ˆ …I

Bp †…Ma

Mb † M a

1

and and

Dpb ˆ Fpb ˆ

…I …I

Ap †…Ca Bp †…Ma

Cb † 1 Cb ; 1

Mb † Mb

…30† …31†

for p ˆ a; b, which completes description of the homogenization method for a composite with elastoplastic coecients by the use of the TFA. It should be underlined that, in comparison to homogenization model valid for linear elastic range, the necessity of transformation matrices computations is crucial for the proposed nonlinear FEM analysis. 2.1. Governing equations of elastoplasticity problem Therefore, the following boundary value problem is to be solved: Drkl;l ˆ 0;

x 2 X;

D~ rkl ˆ Cklmn Demn ;

…32† x 2 X;

1 Demn ˆ ‰Duk;l ‡ Dul;k ‡ ui;k Dui;l ‡ Dui;k ui;l ‡ Dui;k Dui;l Š; 2 with the following boundary conditions:

…33† x 2 X;

…34†

226

M. Kaminski, è. Figiel / Computational Materials Science 22 (2001) 221±239

Drkl  nl ˆ Dtk; uk^; Duk^ ˆ D^

k ˆ 1; 2; 3;

x 2 oXr ;

…35†

k^ ˆ 1; 2; 3:

x 2 oXu ;

…36†

The problem is solved for displacements uk …x†, strain ekl …x† and stress rkl …x† tensor components ful®lling equilibrium equations (32)±(36). Let us note that the stress tensor increments Drkl …x†, D~ rkl …x† denote here the ®rst and the second Piola±Kirchho€ tensors Drkl ˆ DFkm D~ rml ‡ Fkm D~ rml ‡ DFkm r~ml ;

x 2 X;

…37†

where DFkm ˆ Duk;m ;

x 2 X:

…38†

To get the solution, the following functional de®ned on the displacements increments as Duk is introduced:  Z  Z 1 1 ~ Cklmn Dekl Demn ‡ rkl Dui;k Dui;l dX D^tk Duk d…oX†: …39† J …Duk † ˆ 2 2 X oX The discretization of the problem according to the FEM technique is presented in the following section. 3. Numerical analysis 3.1. Finite element analysis Let us introduce the displacement increments function Duk …x† being continuous and di€erentiable on X and, consequently, including all geometrically continuous and coherent subsets (®nite elements) Xe , e ˆ 1; . . . ; E, discretizing the whole X. It is not assumed that Duk …x† is di€erentiable on the interelements surfaces and boundaries oXef (for e; f ˆ 1; . . . ; E, e 6ˆ f ). Next, let us consider the following approximation of Duk …x† for x 2 X: Duk …x† ˆ

Ne X fˆ1

…N †

ufk …x†Duf ;

…40† …N †

where ufk …x† are the shape functions for node k, Duf represents the degrees of freedom (DOF) vector, while Ne is the total number of the DOF in this node. Considering the above, the displacements and strains gradients are rewritten as follows: …N †

Duk;l …x† ˆ ufk;l …x†Duf ; h i …1†f …2†f …N † …N † Dekl …x† ˆ Bkl ‡ Bkl Duf ˆ Bfkl Duf ;

…42†

fn …N † …N † Dekl …x† ˆ Bkl ‡ Duf Dun ;

…43†

…41†

and ®nally Dekl …x† ˆ Dekl …x† ‡ ekl …x†:

…44†

The following notation is applied in Eqs. (42) and (43): …1†f Bkl …x† ˆ ufk;l …x†;

…45†

M. Kaminski, è. Figiel / Computational Materials Science 22 (2001) 221±239

227

…2†f …N † Bkl …x† ˆ ufi;k …x†uni;l …x†un ;

…46†

1 fn Bkl …x† ˆ ufi;k …x†uni;l …x†: 2

…47†

All these equations are substituted into the variational formulation of the problem [cf. Eq. (39)]. There holds 1 1 Cklmn Dekl Demn ˆ Cklmn …Dekl ‡ Dekl †…Demn ‡ Demn † 2 2 1 ˆ Cklmn …Dekl Demn ‡ Dekl Demn ‡ Dekl Demn ‡ Dekl Demn † 2  1 fn …N † …N † …N † lm …N † …N † † …N † ˆ Cklmn Bfkl Duf Bnmn Dun ‡ Bfkl Duf Bmn Du…N ‡ Bkl Duf Dun Blmn Dul…N † l Dum 2  fn …N † …N † …N † lm ‡ Bkl Duf Duf Duf Bmn Dul…N † Dum…N † …48† and 1 1 …N † …N † r~kl Dui;k Dui;l ˆ r~kl ufi;k …x†Duf uni;l …x†Dun : 2 2

…49†

Next, the following notation is applied: Z …r†e …N † kfn ˆ r~kl ufi;k …x†uf uni;l …x† dX;

…50†

Xe

…u†e

kfn

Z

1 …1†f …1†n Cklmn Bkl Bmn dX; Xe 2 Z   1 …1†f …2†n …2†f …1†n …2†f …2†n Cklmn Bkl Bmn ˆ ‡ Bkl Bmn ‡ Bkl Bmn dX; Xe 2

…con†e

kfn

ˆ

…51† …52†

where …1†e

…r†e

…con†e

kfn ˆ kfn ‡ kfn

…u†e

‡ kfn

…53†

and for the second- and third-order terms Z   3 …2†e …N † lm † …N † fn Du…N † Du…N † Bl Du…N † dX; Cklmn Bfkl uf Bmn Du…N kfn ˆ B Du ‡ f n kl l m mn l Xe 2 Z   fn …3†e …N † …N † lm kfn ˆ 2Cklmn Bkl Duf Dun Bmn Dul…N † Dum…N † dX: Xe

…54† …55†

…i†

Introducing kfn for i ˆ 1; 2; 3 to the functional J …Duk † in Eq. (16) and applying transformation from the local system to the global one by the use of the following formula …N †

Duf

ˆ ana Dqa ;

…56†

it is obtained that 1 …1† 1 …2† 1 …3† J …Dqa † ˆ Kab Dqa Dqb ‡ Kabc Dqa Dqb Dqc ‡ Kabcd Dqa Dqb Dqc qd 2 3 4

DQa Dqa :

The stationarity of the functional J …Dqa † leads to the following algebraic equation:

…57†

228

M. Kaminski, è. Figiel / Computational Materials Science 22 (2001) 221±239 …1†

…2†

…3†

Kab Dqb ‡ Kabc Dqb Dqc ‡ Kabcd Dqb Dqc Dqd ˆ DQa

…58†

being ful®lled for any con®guration of X. Iterative numerical solution of this equation makes it possible, according to the speci®ed boundary conditions, to compute the e€ective constitutive tensor components of the homogenized composite. It should be stressed that the ®rst two components of the sti€ness matrix are considered in further numerical analysis only (geometrical nonlinearity is omitted in the homogenization process); detailed description of the numerical integration and solution of Eq. (58) can be found in [3,6], for instance. 3.2. Transformation ®eld analysis As it was mentioned above, the main goal of the transformation-based approach is to compute the transformation matrices Ar , Drs that are computed only once for original geometry of the composite and with the assumption about elastic characteristics of the constituents. There hold eff hdrr iVr ˆ Ceff r hder iVr ˆ Cr

el 

 drr Vr ˆ Celr deelr Vr



deelr ‡ deinel r

and

drinel r

 Vr




; Vr

hder iVr ˆ Ar de ‡

 ˆ Cr deinel ; r Vr

N X r;sˆ1

 Drs deinel ; s Vs

…59† …60†

el Ceff r ˆ Cr ‡ Cr :

…61†

Further, using spatial averaging de®nitions, the averaged stress tensor components are calculated as follows: her iV ˆ e and

hrr iV ˆ r:

…62†

Hence, the e€ective elasticity tensor components Ceff are derived for given increment as dr ˆ Ceff de; Ceff ˆ

N X rˆ1

…63†

cr Celr ‡

N X r;sˆ1

cr Drs einel s




1

: rinel r :

…64†

In the particular case of a two-component composite, the transformation concentration matrices are obtained as 1

…65†

1

…66†

D11 ˆ …I

A1 †…C1

C2 † C1 ;

D21 ˆ …I

A2 †…C1

C2 † C1 ; 1

D12 ˆ

…I

A1 †…C1

C2 † C2 ;

…67†

D22 ˆ

…I

A2 †…C1

C2 † 1 C2 :

…68†

C1 , C2 denote here the components corresponding to elastic properties, while A1 ; A2 are mechanical concentration matrices. Finally, using Eq. (64) it is obtained that Ceff ˆ

N X rˆ2

‡

c1 Cel1 ‡ c2 Cel2 ‡

X

c1 D21 einel 1




1

X

c1 D11 einel 1




1

inel : rinel 1 ‡ c2 D12 e2

inel : rinel 1 ‡ c2 D22 e2




1

: rinel 2 :




1

: rinel 2 …69†

M. Kaminski, è. Figiel / Computational Materials Science 22 (2001) 221±239

229

The FEM aspects of the TFA computational implementation are discussed in detail in Section 4. Further, it should be noticed that there were some approaches in the elastoplastic approach to composites, where, analogously to the linear elasticity homogenization method, the approximation of the e€ective yield limit stresses of a composite was proposed as a quite simple closed form function [17] p …70† r…eff† ˆ r1 r2 or, in terms of e€ective yield surface, in the following form: 9 82 !2 3 1  <  2 = ry…1† U… r† ˆ max 4c2 ‡ c1 ; y 5 r m…y†r r…2† y …2† y P0 : ; ry

…71†

where m…y ˆ l1 =l2 † ˆ 3l2 V …l1 ; l2 † and V is any estimate of the viscosity compliance tensor de®ned starting from the viscosities l1 and l2 [21]. The review of the most recent theories in this ®eld can be found in [17], for instance.

4. Computational illustration The main purpose of the computational experiment presented is to determine the global nonlinear homogenized constitutive law for two component composites with elastoplastic components. The FEMbased program ABAQUS [1] is used in all computational procedures, however the method presented can be implemented in any nonlinear FEM plane strain/stress code as in [19], for instance. The numerical experiments are carried out in the microstructural (RVE) level, and that is why the global behavior of the composite is predicted starting from the behavior of the periodicity cell. The numerical micromechanical model consists of a three-component periodicity cell with horizontal and vertical symmetry axes and the dimensions 3:0 cm …horizontal†  2:13 cm …vertical† (cf. Figs. 3 and 4). The composite is made of epoxy matrix and metal reinforcement with material properties of the components collected in Table 1. The void embedded into the steel casting simulates the lack of any matrix in the periodicity cell. Some nonzero material data are introduced to avoid numerical singularities during homogenization problem solution.

Fig. 3. Cross-section of a superconducting coil.

230

M. Kaminski, è. Figiel / Computational Materials Science 22 (2001) 221±239

Fig. 4. 3D view of the superconducting coil part.

Table 1 Material characteristics of composite constituents No.

Material

Young's modulus

Poisson's ratio

Yield stress

1 2 3

Epoxy resin Metal Void

7000.0 42000.0 70.0

0.3 0.2 0.1

10.0 22.0 0.1

To calculate the e€ective tensor components, the boundary value problem given by Eqs. (22)±(25) is solved ®rst. The periodicity cell is discretized with 25 CGPE10R ®nite elements (10-node biquadratic, quadrilateral hybrid linear pressure reduced integration plane strain element with four integration Gaussian points). Periodic boundary conditions are imposed to assure periodic character of the entire structure behavior. Suitable formulation of displacement boundary conditions has the following form: ui ˆ eij …y…P2 †

y…P1 ††;

…72†

where ui ˆ fu1 ; u2 g represents the displacement function components, eij is the global total strain imposed to the periodicity cell, while y…P1 † and y…P2 † are the coordinates of the points lying on the opposite sides of the RVE. The displacement boundary conditions are introduced at the edges of RVE quarter as it is shown in Figs. 5 and 6. Further, since the generalized plane strain is considered, the matrices computed have a rank a ˆ 4 and the total dimensions of the matrices Ar and Drs are ‰4  4Š. The ®rst step in the numerical analysis is to compute mechanical and transformation concentration matrices Ar and Drs , which is carried out according to the implementation in computer system ABAQUS (cf. Fig. 7). The transformation matrices Ar and Drs are evaluated as follows: T 1. Matrix Ar by means of the overall strain loading case e ˆ feij g ˆ ‰e11 ; e22 ; 2e12 ; e33 Š introduced as displacements ui ˆ eij yj :

…73†

2. Matrix Drs imposing the uniform eigenstrain in the subvolume Vr or Vs as the uniform stress; since it is not possible to introduce the eigenstrain directly in each subvolume in the program ABAQUS, the stress tensor components are calculated as rr ˆ

Cr : er

and

r ˆ 1; . . . ; N

…74†

M. Kaminski, è. Figiel / Computational Materials Science 22 (2001) 221±239

231

Fig. 5. Boundary conditions for E11 6ˆ 0.

Fig. 6. Boundary conditions for E22 6ˆ 0.

and applied to each of N subvolumes, where the elasticity tensor Cr is given as 2 3 vr vr 0 1 vr 6 vr 1 vr vr 0 7 Er 6 7 Cr ˆ 6 7: …1 vr †…1 2vr † 4 vr vr 1 vr 0 5 0

0

0

…75†

1 2vr 2

The accuracy of homogenization method applied for a given material model is veri®ed by the comparison with the results obtained for the real heterogeneous composite under the same boundary conditions. For this purpose, the same boundary value problem is solved with four di€erent loading cases. The elastoplastic static analysis consists of 25 incremental load steps (with constant increment in each step) and is performed using the Radial Return algorithm for the perfect J2 elastoplastic material [6,19]. The results in the form of stress±strain relation are shown in Figs. 8±11, while the stresses distribution in the periodicity cell can be compared in Figs. 12±15. Generally, it is observed that the elastic range is very well approximated by TF model results. However, the homogenized material seems to be a little sti€er than the heterogeneous

232

M. Kaminski, è. Figiel / Computational Materials Science 22 (2001) 221±239

Fig. 7. Flowchart of the numerical analysis.

one, especially in nonlinear range in the direction y1 of the RVE (cf. Fig. 3). In the same time, for the interrelation of shear strain and stress, the last incremental steps show almost linear behavior and that is why there are no di€erences between heterogeneous and homogeneous material. To obtain e€ective elastoplastic properties more eciently, the homogenization method presented above should be corrected to include the increments of transformation matrices during loading process. Further, the e€ective properties of the homogenized material are computed starting from the properties of composite constituents and the constitutive relation veri®ed for all strain increments during the computational incremental analysis. For this purpose the relation (69) is used and therefore

M. Kaminski, è. Figiel / Computational Materials Science 22 (2001) 221±239

Fig. 8. Constitutive r11 ±e11 relation for homogenized and real composites.

Fig. 9. Constitutive r22 ±e22 relation for homogenized and real composites.

Fig. 10. Constitutive r33 ±e33 relation for homogenized and real composites.

233

234

M. Kaminski, è. Figiel / Computational Materials Science 22 (2001) 221±239

Fig. 11. Constitutive r12 ±e12 relation for homogenized and real composites.

Fig. 12. The equivalent stress req 11 distribution in RVE.

2

k‡l 6 k l Cel ˆ 6 4k l 0

k l k‡l k l 0

k l k l k‡l 0

3 2 0 15776:83 2601:7 6 2601:7 15776:83 07 7ˆ6 0 5 4 2601:7 2601:7 l 0 0

while the inelastic part of the e€ective constitutive 2 3393:88 1505:92 6 749:44 2653:36 inel inel 1 6 C ˆr:e !C ˆ4 1631:4 1546:4 0 0

2601:7 2601:7 15776:83 0

tensor can be written as 3 2080:0 0 7 1576:12 0 7: 5 5017:2 0 0 499:600

3 0 7 0 7; 5 0 6587:57

M. Kaminski, è. Figiel / Computational Materials Science 22 (2001) 221±239

Fig. 13. The equivalent stress req 22 distribution in RVE.

Fig. 14. The equivalent stress req 12 distribution in RVE.

Finally, the overall e€ective constitutive tensor Ceff is obtained as 2 3 19170:71 4107:62 4681:7 0 6 3351:14 18430:2 4177:82 7 0 7; Ceff ˆ 6 4 4233:1 5 4148:1 20794:03 0 0 0 0 7087:17 what completes the calculations of e€ective elastoplastic characteristics of the composite considered.

235

236

M. Kaminski, è. Figiel / Computational Materials Science 22 (2001) 221±239

Fig. 15. The equivalent stress req 33 distribution in RVE.

5. Conclusions 1. As it is shown in the computational experiments, the homogenized material obtained as a result of the transformation ®eld analysis (TFA) is sti€er in nonlinear range than the original composite. It is caused by the fact that the constitutive relations during the whole iteration procedure are based on constant constitutive tensors Cr with elastic properties. To obtain better e€ective approximation of composite behavior, these matrices should be divided into elastic and plastic part, by means of the consistent tangent matrices. 2. Since the TFA makes it possible to characterize explicitly the e€ective elastoplastic behavior starting from composite components material properties, it is possible to carry out the numerical sensitivity studies of homogenized composite properties with respect to its original material characteristics. Such computational studies make it possible to determine the most crucial material parameter for the overall elastoplastic behavior of the composite, which may be important in the context of optimization techniques applied in composite engineering design studies. 3. Due to the fact that most of composite components material characteristics are obtained experimentally as statistical estimators, the next step to utilize the presented approach is probabilistic implementation of the homogenization problem. It will generally enable to compute respective probabilistic moments and coecients of e€ective properties, starting from the expected values and standard deviations of composite components elastoplastic characteristics. As it is well known, it can be done using the Monte Carlo simulation technique [16], for instance. Further, it should be mentioned that such an implementation makes it possible to determine the stochastic sensitivity of composite e€ective characteristics to randomness of component materials elastoplastic behavior.

Acknowledgements The ®rst author would like to acknowledge the ®nancial support of the Foundation for Polish Science in the framework of postdoctoral fellowship at Rice University in Houston, TX, USA.

M. Kaminski, è. Figiel / Computational Materials Science 22 (2001) 221±239

237

Appendix The two-component transversely isotropic RVE of volume V is subjected to a uniform overall strain increment DE or stress DR. A possible description of the local uniform and stress increment ®eld is suggested as Der ˆ Ar DE;

r ˆ 1; 2;

…A:1†

Drr ˆ Br DR;

r ˆ 1; 2;

…A:2†

with further relations c1 A1 ‡ c2 A2 ˆ I;

…A:3†

c1 B1 ‡ c2 B2 ˆ I:

…A:4†

Further, the local and overall increments are expressed as DR ˆ c1 Dr1 ‡ c2 Dr2

and DE ˆ c1 De1 ‡ c2 De2 ;

…A:5†

with composite constituents volume fractions c1 ˆ

V1 ; V

c2 ˆ

V2 V

and

V1 ‡ V2 ˆ V :

…A:6†

The constituent properties in the elastic and plastic range are de®ned as Drr ˆ Cr Der

Der ˆ Mr Drr and Mr ˆ Cr 1 ;

…A:7†

while the overall properties as DR ˆ CDE

and

DE ˆ MDR; M ˆ C 1 :

…A:8†

The constitutive and compliance matrices are given as spatial averages over the RVE Cˆ

2 X

cr C r A r

and

Mˆ

rˆ1

2 X

cr Mr Br :

…A:9†

rˆ1

The individual components of Br and Ar may be found as the solutions of Eqs. (A.3) and (A.4). There holds 2 3 0; 5…k=kr ‡ l=lr † 0; 5…k=kr l=lr † 0; 5…l lr †=kr 0 6 0; 5…k=k l=lr † 0; 5…k=kr ‡ l=lr † 0; 5…l lr †=kr 0 7 r 6 7 …A:10† Ar ˆ 6 7; r ˆ 1; 2; 4 0 0 1 0 5 0 as well as

2

EC 6 0 1 6 Br ˆ 6 EC 4 …1 cr †ar 0

0

0 EC …1 cr †ar 0

0 0 ErL 0

0 3 0 0 7 7 7; 0 5 EC

l=lr

r ˆ 1; 2;

…A:11†

where EC ˆ c1 E…1†L ‡ c2 E…2†L

and

a1 ˆ …m…1†L E…2†L

m…2†L E…1†L † ˆ

a2 :

…A:12†

238

M. Kaminski, è. Figiel / Computational Materials Science 22 (2001) 221±239

Next, the transformation and concentration matrices Drs , Frs are calculated as Dr1 ˆ …I

Ar †…C1

C2 † 1 C1 ;

Fr1 ˆ …I

Br †…M1

M2 † 1 M1 ;

Dr2 ˆ …I

Ar †…C1

Fr2 ˆ …I

Ar †…M1

C2 † 1 C2 ;

r ˆ 1; 2;

M2 † 1 M 2 ;

…A:13†

r ˆ 1; 2:

…A:14†

where D11 ; D12 ; D21 ; D22 ; F11 ; F12 ; F21 ; F22 are to be calculated. The components of the matrix D11 are for instance obtained as 2 3 D1111 D1112 D1113 D1114 6 D1121 D1122 D1123 D1124 7 7 …A:15† D11 ˆ 6 4 D1131 D1132 D1133 D1134 5; D1141 D1142 D1143 D1144 then 1

k 2k1

2  l 2l1 …k1 …2k1 ‡ l1 †4 2l1

1

k 2k1

2  l 2l1 …k1 …2k1 ‡ l1 †4 2l1

 D1111 ˆ  D1112 ˆ  D1113 ˆ 1  D1121 ˆ

1 

D1122 ˆ

1

2

k2 † ‡ 2l2 …k2 k1 † ‡ …k1 k2 † ‡ …l1   2 2 …k1 k2 † ‡ …l1 l2 † …l1 l2 † 2

k2 † ‡ 2l2 …k2 k1 † ‡ …k1 k2 † ‡ …l2   …k1 k2 †2 ‡ …l2 l1 †2 …l1 l2 †    k1 l1 l2 ‡ 2…n1 n2 † …l l1 † …2l1 ‡ n1 † ; …l1 l2 †…n1 n2 † 2 2  2 k l 2l1 …k1 k2 † ‡ 2l2 …k2 k1 † ‡ …k1 k2 † ‡ …l2   ‡ …2k1 ‡ l1 †4 2k1 2l1 …k1 k2 †2 ‡ …l2 l1 †2 …l1 l2 † 2  2 k l 2l1 …k1 k2 † ‡ 2l2 …k2 k1 † ‡ …k1 k2 † ‡ …l1   …2k1 ‡ l1 †4 2 2 2k1 2l1 …k1 k2 † ‡ …l l † …l1 l2 † 1

 D1123 ˆ 1 D1144

 ˆ l1 1

   l1 l2 ‡ 2…n1 n2 † l1 † …2l1 ‡ n1 † ; …l1 l2 †…n1 n2 †   l 1 : l1 l1 l2

k1 …l 2

3 2 l2 † 5 ;

…A:16†

3 2 l1 † 5 ;

…A:17†

…A:18† 3 2 l1 † 5 ;

…A:19†

3 2 l2 † 5 ;

…A:20†

2

…A:21† …A:22†

The remaining coecients of D11 are equal to 0. References [1] ABAQUS User's Manual, v. 5.8, Hibbitt, Karlsson & Sorensen Inc., Pawtucket, 1999. [2] J. Aboudi, Elastoplasticity theory for composite materials, Solid Mech. Arch. 11 (1986) 141±183. [3] N. Aravas, On the numerical integration of a class of pressure-dependent plasticity models, Int. J. Num. Meth. Eng. 24 (1987) 1395±1416. [4] Y.A. Bahei-El-Din, Finite element analysis of viscoplastic composite materials and structures, Int. J. Comp. Mat. Struct. 3 (1) (1996) 1±28. [5] Y. Benveniste, A new approach to the Mori±Tanaka's theory in composite materials, Mech. Mat. 6 (1987) 147±157. [6] M.A. Cris®eld, Non-linear Finite Element Analysis of Solids and Structures ± Essentials, Wiley, New York, 1991.

M. Kaminski, è. Figiel / Computational Materials Science 22 (2001) 221±239

239

[7] L.N. Philips (Ed.), Design with Advanced Composite Materials, Springer, Berlin, 1989. [8] G.J. Dvorak, Y.A. Bahei-El-Din, A.M. Wafa, Implementation of the transformation ®eld analysis for inelastic composites materials, Comput. Mech. 14 (1994) 201±228. [9] G.J. Dvorak, Y. Benveniste, On transformation strains and uniform ®elds in multiphase elastic media, Proc. Roy. Soc. London A 437 (1992) 291±310. [10] G.J. Dvorak, Transformation ®eld analysis of inelastic composite materials, Proc. Roy. Soc. London 437 (1992) 311±327. [11] U. Galvanetto, C. Pellegrino, B.A. Schre¯er, Plane stress plasticity in periodic composites, Comput. Mat. Sci. 624 (1998) 1±11. [12] R. Hill, Elastic properties of reinforced solids ± some theoretical principles, J. Mech. Phys. Sol. 11 (1963) 55±372. [13] R. Hill, A self-consistent mechanics of composite materials, J. Mech. Phys. Sol. 13 (1965) 213±222. [14] R. Hill, The essential structure of constitutive laws for metal composites and polycrystals, J. Mech. Phys. Sol. 15 (1967) 79±95. [15] M. Kami nski, Homogenized properties of the n-component composites, Int. J. Eng. Sci. 38 (4) (2000) 405±427. [16] M. Kami nski, M. Kleiber, Numerical homogenization of n-component composites with stochastic interface defects, Int. J. Num. Meth. Eng. 47 (2000) 1001±1025. [17] J.C. Michel et al., E€ective properties of composite materials with periodic microstructures: a computational approach, Comput. Meth. Appl. Mech. Eng. 172 (1-4) (1999) 109±144. [18] S. Nemat-Nasser, M. Hori, Micromechanics: Overall Properties of Heterogeneous Materials, North-Holland, Amsterdam, 1993. [19] D.R.J. Owen, E. Hinton, Finite Elements in Plasticity Theory and Practice, Pineridge, Swansea, 1980. [20] J.C.J. Schellekens, Computational Strategies for Composite Structures, TU Delft, 1990. [21] M. Zaidman, P. Ponte-Castaneda, E€ective yield surfaces for anisotropic composite materials, in: D.F. Parker, A.H. England (Eds.), IUTAM Symposium on Anisotropy, Inhomogeneity and Nonlinearity in Solid Mechanics, Kluwer, Dordrecht, 1995, pp. 415±422.