- Email: [email protected]
Small deformation multi-scale analysis of heterogeneous materials with the Voronoi cell finite element model and homogenization theory
COMPUTATIONAL MATERIALS SCIENCE ELSEVIER
Computational
Materials Science 7 (1996) 13 I-
146
Small deformation multi-scale analysis of heterogeneous materials with the Voronoi cell finite element model and homogenization theory Kyunghoon Lee Applied
Mechanics
Section, Department
ofAerospace Engineering,
‘,
Somnath Ghosh Applied
Mechanics
* ,*
and Auiution,
The Ohio State University,
Columbus,
OH 43210, USA
Abstract In this paper, a multiple scale finite element model (VCFEM-HOMO) has been developed for elastic-plastic analysis of heterogeneous (porous and composite) materials by combining asymptotic homogenization theory with the Voronoi cell finite element model (VCFEM). VCFEM for microstructural modeling originates from representative material elements at sampling points in the structure. Structural modeling is done by the general purpose finite element code ABAQUS, and interfacing with the microscale VCFBM analysis is done through the user subroutine in ABAQUS for material constitutive relation, UMAT. Asymptotic homogenization in UMAT generates macroscopic material parameters for ABAQUS. Keywords:
Voronoi cell finite element method; Heterogeneous materials; Asymptotic homogenization
1. Introduction The last three decades we have experienced a surge in the advancement of science and technology for heterogeneous materials that contain dispersions of multiple phases in the microstructure. Examples are metal/alloy systems with second phase in the form of precipitates and pores, and composite materials containing a dispersion of fibers, whiskers or particulates in the matrix. Rigorous fundamental studies are needed for understanding deformation and failure mechanisms prior to design of advanced materials in high performance applications. These models establish the effect of second phase shapes,
* Corresponding
author. Tel.: + 1-614-2922599; e-mail: [email protected]. ’ Graduate Research Associate. * Associate Professor.
2927369;
0927-0256/96/$15.00 PII
Copyright
SO927-0256(96)00072-9
fax: + I -614-
sizes and distributions on evolving state variables like stresses, plastic strains, void initiation and growth, and evolving material variables like strain hardening and flow stress. Ghosh and coworkers have developed a material based Voronoi cell finite element model (VCFEM) [l-5] in an attempt to overcome difficulties in modeling arbitrary microstructures by conventional finite element methods. The VCFEM mesh evolves naturally from a heterogeneous microstructural material element (RME) by Dirichlet tessellation into a network of multi-sided ‘Voronoi’ polygons containing one second phase heterogeneity each, at most. The multi-phase Voronoi polygons identified with the material’s basic structural elements, constitute elements in VCFEM. VCFEM with asymptotic homogenization for elastic problems have been presented in Ref. [3]. This paper presents a coupled multiple scale computational model for heterogeneous elastic-plas-
0 1996 Elsevier Science B.V. All rights reserved.
132
K. Lee. S. Ghosh / Computationul
tic structures. Only two dimensional problems are considered. Microstructural analysis for various different representative material element arrangements, is done with the Voronoi cell finite element model. The commercial general purpose code ABAQUS [6] is used for global analysis at the level of overall structural geometry and applied loads. The interfacing between macro- and micro-calculations are done through the user subroutine window UMAT in ABAQUS. Numerical examples are conducted to investigate the advantages of coupled multiple scale analysis over other unit cell and continuum based theories.
Muterids
Science 7 (1996) 131-146
element boundary. Also let Au correspond to an equilibrated stress increment in Y,, Au to a compatible displacement increment on aY,, and At to a traction increment on the traction boundary &,. The incremental problem is solved by using a two field assumed stress hybrid variational principle, derived from an element energy functional as: &(Ao, =-
AB(cr,Aa)dY-je:AodY / Y, r, +
/ AY(
a+Au).n’.(u+Au)aY
-,r’(i+ai)(u+Au)dl In?
2. Microstructural analysis with VCFEM Voronoi cells, resulting from Dirichlet tessellation of a heterogeneous microstructure, make rather unconventional elements, due to the arbitrariness in the number of edges. The Voronoi cell finite element model developed by Ghosh and coworkers [l-4,7] avoids difficulties of conventional displacement based FEM formulations by invoking the assumed stress hybrid method, introduced by Pian [8]. In this formulation, independent assumptions are made on an equilibrated stress field in the interior of each element and a compatible displacement field on the element boundary. Small deformation elastic-plastic analysis of materials with embedded second phase has shown significant promise with respect to efficiency and accuracy. It is based on a hybrid formulation originally proposed by Atluri and coworkers [9,10]). Details of this development are presented in Refs. [ 1,2]. Consider a typical representative material element (RME) Y tessellated into N Voronoi cells, as shown in Fig. lb. This is based on the location, shape and size of N heterogeneities as explained in [5]. The matrix phase in each Voronoi cell Y, is denoted by Y, and the heterogeneity (void or inclusion) is denoted by Y,. The matrix-heterogeneity interface aY, has an outward normal nc, while ne is the outward normal to the element boundary JY,. An incremental finite element formulation is invoked to account for rate independent plasticity. At the beginning of the p-th increment, let CTbe an equilibrated stress field with a strain field e( 0, loud history), and u be a compatible displacement field on the
Au)
- / av,(
urn+,,“‘-u’-Au’)
nc . (u’ + Au’) c?Y,
(1)
where Au’ is the displacement of the interface and AB is the increment in element complimentary energy. Superscripts m and c represent respectively the matrix and second phase parts of the Voronoi cell element. The energy functional for the entire domain is obtained by adding each element functional as fl= Cr= , Lr,. The first variation of fl, with respect to the stress increments A u, results in the kinematic relations as the Euler equation, while the first variation of n with respect to boundary displacement increments Au yields traction conditions as Euler equations. Independent assumptions on stress increments A u are made in the matrix and heterogeneity phases to accommodate stress jumps across the interface. This can be done using stress functions @(x, y) for the matrix and inclusion phases, resulting in stress expressions as: {Au”)
= [P”‘(L
y)]IApm}
{Au’] = [Pc( x, Y)]{W)
t
(2)
where (Ap)‘s correspond to a set of yet undetermined stress coefficients and [PI is a matrix of interpolation functions. Compatible displacement increments on the element boundary aY, as well as on the interface aY,, are generated by interpolation in terms of generalized nodal values. The displacement
K. Lee, S. Ghosh / Computatiod
Marerials Science 7 (1996) I3/-146
(b)
133
(cl
Fig. I. A heterogeneous structure with various levels. (a) The global structure and computational geometry, (b) representative material elements (RME) at a point and the corresponding VCFJ3 model and, (cl a basic structural element (BSE) represented by a Voronoi cell.
increments on the element boundary and interface may then be written as, {Au] = [L’](Aq] 9
(Au’}
= [L’](Aq’},
(3)
where {Aq} and (Aq’) are generalized displacement increment vectors and [L] is an interpolation matrix. Substituting element approximations for stresses Eq. (2) and displacements Eq. (3), in the energy functional Eq. (l), and setting the first variations with respect to the stress coefficients A p m and A /3’ respectively to zero, results in the following two weak forms of the kinematic relations, /,
PITI
:,
[Pm]T[n’][L’]
0
I? c= I
I,
[L’f[nc]T[Pc]dY c
/3” +Aj3”
,f3’+ Ap’
For an elastic-plastic material, the strain increments be in Eq. (4) are non-linear functions of the current state of stress u as well as of their increments A u. The non-linear finite element Eq. (4) and Eq. (5) are solved for the stress parameters (A p m, A /3“) and the nodal displacement increments (Aq, Aq’) in the p-th increment.
e+Ae}dY
m
= /
the weak form of the traction reciprocity conditions as
dY{Aq}
[Pm]T[nC][LC]dY{ Aq’), arc
2.1. Shape based stress representations /r,[PclT{e + Ae} dY = i, <
[Pc]T[nC][L’]dY{Aq’). c (4)
Setting the first variation of the total energy functional with respect to Aq and Aq’ to zero, results in
An important criterion, affecting the convergence of multiple phase Voronoi cell elements, is the proper choice of stress fields in each of the constituent phases. Choosing stress functions from micromechanics considerations, adds considerably to the ele-
K. Lee, S. Ghosh/ Computational Materids
134
Science 7 (19%)
131-146
Initial Global &Local Information
p=l
UMAT (element #l)
sm
/ Minocopic VCFEM-..
Microcopic VCFEM-..
IPI
El II
~
L-
p=p+l
Fig.2. Flowthan
ment efficiency. Three different conditions that are indispensable in this regard are: 1. Stress functions should, in some way, account for the shape of the heterogeneity. 2. Effects of the heterogeneity shape should vanish
of
microscopic analysis.
at large distances from the interface for matrix stress functions. 3. Shape effects in matrix stress functions should facilitate traction reciprocity at the interface. Pure polynomial forms Airy’s stress functions do
‘7 (4
Fig. 3. (a) ANSYS and (b) VCFEM meshes for randomly
r*
packed RME with circular voids CV, = 20%).
K. Lee, S. Ghosh/ComputationulMaterialu
not explicitly account for the shape of the heterogeneity and require very high order terms for convergence. To deal with arbitrary shaped interfaces in VCFEM, a special procedure is pursued. Suppose that equation of the interface aY, in Fig. lc can be expressed in polar coordinates as ,g(r, 0) = 0, where the r coordinate is measured from the centroid of the heterogeneity. A Fourier series expansion for r in terms of the polar angle 13may be expressed as: r=a,+
Ca,cos(nB) n
+ Cb,sin(&) n
Science 7 (1996) 131-146
135
6.0934
R,, = 1.0
on JY,, (6)
where a, and b, are the Fourier coefficients. The interface equation may then be expressed from Eq. (6) as, g( r, 0) =f-
;
- c %os( n a0
- F :sin(no)
ne)
=O.
(4
(7)
0
Here f corresponds to a function that transforms any arbitrary shaped interface to an approximate unit circle, since fl X, y) = 1 on JY,, with the property that l/f+ 0 as (x, y) + x. This function is now used to construct stress functions in the matrix region Y, from the two part expression Grn = @pomly + @e”,,
-
Fig. 5. (a) ABAQUS model for thick various RME’s at 40% volume fraction, 6) square diagonal packing, (d) random heterogeneities and (e) random packing geneities.
(8)
expressed in terms of a polynomial expanwith @pOmly (a) 0.050
plate with circular holes, (b) square edge packing, packing with 15 circular with 15 elliptical hetero-
ANSYS
.._.._._... ,,CFEM 8 =
0.040
-Ti -”
v) 2 e 2 .Y
0.030
I
0.0x
,
-
0.01c
,
-
WVCFEM D-OANSYS
a 2 H0
O.oOC1 0.0
Fig. 4. (a) Macroscopic x/L= 0.5.
0.2
Macroscopic
stress-strain
0.4
0.6
0.8
-0.02
’ 0.0
1 0.2
tensile strain(%)
response
and (b) Microscopic
0.4
0.6
0.8
1.0
Y/L
stress distribution
at 0.8% strain for random
packed
RME along
K. Lee, S. Ghosh / Cmtpututiortd
136
Matrriuk Scxence 7 (1996) 131-146
(b)
(4 Fig. 6. Von Mises stress distribution
for porous material with square edge packing (a) macroscopic
sion of coordinates, i.e. GP’& = ~,,,y~py~pyy. For each polynomial term in @& there exist shape based reciprocal terms cP~~, that give rise to stresses equilibrating the traction field on a Y,. This is written as:
=
c xy$fyi, . p,q,i
The inclusion is modeled using polynomial stress functions only. A detailed treatment of stress func-
stress, (b) microscopic
tion selection and convergence Ref. [2]. 2.2. Constitutive
stress at point A.
issues is presented in
relations
A rate independent small deformation elastic-plastic constitutive relation, following J, flow theory with isotropic hardening is considered in this paper. An additive decomposition of the strain increments Ae into an elastic part Ae’ and a plastic part her’ is assumed. The yield surface in stress-space at the beginning of the p-th increment is expressed as, Y”(e,
Ae) = \/$v’:u’,
(9)
2.472EG
hln.J+
(b)
(a) Fig. 7. Von Mises stress distribution A.
for porous material with square diagonal packing (a) macroscopic
strain, (b) microscopic
strain at point
K. Lee, S. Ghosh / Compurutional Muterids
where g’ is the deviatoric stress and Y(e, Ae) is the radius of the flow surface. The plastic strain increment Aep’ is obtained by numerically integrating the flow rule by the backward Euler method to yield: AeP’=Ah(a+Aa)‘,
(10) where AA is a non-negative incremental flow parameter. Since Ae is in general a function of AA, u and A g, the flow surface radius can be expressed in the form, YP+’ = Yp+‘(AA, u, Au) and the flow parameter AA can be evaluated from the following relations: A A = 0 if u$+ I < YP (elastic unloading) P+ ffeff
I =
yP+
I
if
(+P+ eff
I > -
,
(11)
YP
(neutral and plastic loading) .
(12)
Details on numerical implementation are provided in
[1,21. 2.3. A numerical example with VCFEM The effectiveness of the elastic-plastic Voronoi cell finite element model in analyzing heterogeneous microstructures is established by several numerical examples in [1,2]. In this example, the representative material element (RME) consists of 29 randomly located circular voids with volume fraction V, = 20%. Results generated by VCFEM are compared with those of a displacement based commercial finite
Fig. 8. Von Mises stress distribution A.
Science 7 (1996) 131-146
137
element code ANSYS with a very high resolution mesh necessary for convergence. Fig. 3b and a show the VCFEM mesh with 29 elements and the ANSYS mesh with 5282 QUAD4 elements, respectively. Plane strain uniaxial tension loading is considered to a maximum of 0.8% macroscopic strain. Periodic strain uniaxial tension loading is considered to a maximum of 0.8% macroscopic strain. Periodic construct of the RME is achieved through repeatability conditions on the traction free face. The matrix material is Aluminum with the following properties: E = 69 GPa (Young’s modulus); v= 0.33 (Poisson’s ratio); Y, = 43 MPa (initial yield stress); o& = Y + e”.‘25 (post yield hardening law). Compzson of macroscopic response in Fig. 4a clearly establishes VCFEM as an accurate method for modeling overall behavior of random microstructures at a significant advantage in computing efficiency. The microscopic stress distribution through the section A-A in Fig. 3b at 0.8% strain is depicted in Fig. 4b. Once again the comparison is very satisfactory, with VCFEM producing similar patterns and peak stresses with the highly refined ANSYS model. 3. Asymptotic
homogenization
Consider a heterogeneous body occupying a rein Fig. la, for which the microstrucgion oSrlUCtUTe ture constitutes of spatially periodic representative material elements @ME’s) as shown in Fig. l(b). In real heterogeneous materials, dimensions of the RME
for porous material with random circular packing (a) macroscopic
stress, (b) microscopic
stress at point
K. Lee, S. Ghosh/
138
Computdod
Mutcriu1.v
of characteristic length 1 are typically very small in comparison with the dimensions of the body of characteristic length L. The ratio of these microscopic and macroscopic scales I/L is represented by a very small positive number E. A high level of heterogeneity in the microstructure causes a rapid variation of evolutionary variables e.g. deformation and stresses in a small neighborhood E of the macroscopic point x. This corresponds to a microscopic scale X/E and consequently, all variables are assumed to exhibit dependence on both length scales i.e. @‘= @(x, X/E). The superscript E denotes association of the function with the two length scales. In this notation, 0’ corresponds to a connected domain that extends the structural domain to its microstructure. In most of the work on homogenization theory [l l- 161, a periodic repetition of the microstructure about a macroscopic point x has been assumed, thereby making the dependence of the function on y (= X/E), periodic. This characteristic is often termed as Y-periodicity. where Y corresponds to a RME. For small deformation elasto-plasticity, the rate or incremental forms of the equilibrium equation, kinematic relation, and constitutive relations are given as:
Science
7 (19961
131-146
where tif = tif(x, y) is a Y-periodic rate of displacement field in y. Furthermore the boundary conditions are assumed to satisfy the following equations on the prescribed traction and displacement boundaries, respectively. ci,fnj = ii on r,,
(‘6)
ST=Gi on c,,
(‘7) where n is the unit normal to the boundary. In homogenization theory, the Y-periodic displacement rate or increment field is approximated by an asymptotic expansion with respect to parameter E: C’(x) = tiO(x, y) + Eil’(X, y) + EW(X, y) + . . . . y=;
(18)
The stress rate tensor c+~;.S can then be expressed as I (i,; = -&;;.9 + &;.f + E&j + E*&;; . . . ,
(19)
E
where
(13) in a’, ki;.5= ETjk, S;,
in L?‘,
A Fig. 9. Van point A.
Mises stress distribution for
(‘4)
ti,;: = EFjk,
(15)
Putting the expansion of &;if Eq. (19) in the rate
(b)
(8) porous
material with
(20)
random elliptical packing (a) macroscopic strain, (b)
microscopicswainat
K. Lee, S. Ghosh/C.ompututionalMaterials
form of equilibrium Eq. (13), and setting each coefficientofei,(i=--1,0,1,2 ,... )tozero,resultsin the following set of equations:
Science 7 (19%) 131-146
a Y-periodic function representing characteristic modes of the RME. Substituting Eq. (26) in Eq. (20) yields the microscopic constitutive relations as:
adi; =o, aYj
t+;‘(y)
a&.‘-.‘. ac$ I’+ =o, axj aYj
ax,"' = E,;pm T;; + aY, [
1
(28)
’
where Tit’ is a fourth order identity tensor. The mean of Eq. i 28) yields the homogenized elastic-plastic
ah.2 ski)
1’+ aYj
139
x+f,=o.
(21)
I
Eq. (21) and Eq. (20) lead to the trivial value for &$, and therefore establishes that rio is only a function of x as shown in [ 161, i.e. ki;.s= 0 and tip = L$‘(x) .
(4 0.10
(22)
Substituting Eq. (22) into the second of Eq. (21) leads to the Y-domain equilibrium equation
a&.'.
--‘/ -=
(23)
0.
aYj
From Eq. (15) and Eq. (191, by neglecting the terms associated with E or higher, the constitutive relation in Y is expressed as
random elllpse(40% -
square edge@o%
void) void)
0 - - * square dlagona&?C% e--a
(24)
Macroscopic
random
circle(20%
void) void)
strain(%)
where
ati: f-. JY,
(b) 0.25
(25) 0
Here t;, is the local or microstructural strain rate tensor, for which Sk, = &i/ax,) is an averaged macroscopic part, and k& = &$/ay, is denoted as a fluctuating strain rate tensor (Suquet [17]). Due to linearity of the rate (incremental) problem, di;.f and I$ can be expressed in the forms
b
0.20
1 B .;
0.15
f s
0.10
8 P 0.05
L I
where
0.00
0.0
aeig'(y)
-
= 0 (microscopic equilibrium).
0.1
0.2 Macroscopic
0.3
0.4
0.5
strain(%)
(27)
aYj
In Eq. (26), eiijk’is a Y-periodic function and x/’ is
Fig. 10. Evolution of effective stress with evolving strain at point A in the heterogeneous plate with 40% second phase volume fraction for (a) porous material (b) composite material.
K. Lee. S. Ghosh / Compututional
140
Fig. 11. Von Mises stress distribution point A.
for composite
material
Muterids
Science 7 (19961 131-146
with square edge packing
tangent modulus, for use in the macroscopic analysis, in the form 1 E,& = ( $;‘> = IyI ,c?;;’ dY /
aXj
3.1. Macroscopic
stress, (b) microscopic
stress at
tions for an elastic-plastic problem heterogeneous material as: aiij -=
dY.
(a) macroscopic
-.i
in
%ructure
T
(29)
equations
Sijnj = ii on r,,
The mean of the third equation in Eq. (21) on Y yields the macroscopic form of the governing equa-
lip=Gi
onr,,
(30)
2 OWE41
172lE-01
1,UlEO)
b.UMWQ
(4 Fig. 12. Von Mises stress distribution point A.
for composite
material with square diagonal
packing (a) macroscopic
strain, (b) microscopic
strain at
K. Lee, S. Ghosh/Computationd
Materials Science 7 (1996) 131-146
141
Fig. 13. Von Mises stress distribution for composite material with random circular packing (a) macroscopic stress, (b) microscopic stress at point A.
corresponds to the global domain where%mCtUre with f, and r, as traction and displacement boundaries respectively, % = (k ‘) and ire are the averaged macroscopic stress rate tensor and the displacement rate respectively. Small deformation analysis is pursued at the macroscopic level with the Eq. (301, and with the Voronoi cell finite element model (VCFEM) for the microscopic level.
3.2. Elastic-plastic homogenization with VCFEM In an incremental formulation, the equilibrated microscopic stress increment corresponds to A u ‘( = A u “1 in Eq. (26) and the microstructural strain increments are designated as Ae’ in Eq. (25). Similarly, the increments in microscopic displacements on the cell boundaries aY, are identified with
Fig. 14. Von Mises stress distribution for composite material with random elliptical packing (a) macroscopic strain, (b) microscopic strain at point A.
K. Lee, S. Ghosh/Compututiod
142
Au’
in Eq. (26) and those on the interface are denoted by Au”. In the absence of traction boundaries due to periodicity conditions, the incremental energy functional for each Voronoi cell element in Eq. (1) is modified for the homogenization process as:
-
j( a*‘
uifm + Api;” - c_$’ - Au;;~)
x(u;‘+Auj’)n;dc?Y
Materids
Science 7 (1996) 131-146
scopic module. The first is to evaluate the microscopic stress increments A u ’ from given values of the macroscopic strain G at the beginning of the step, and its increment AZ. This step involves the iterative solution of modified forms of the kinematic relations Eq. (4) and the traction reciprocity conditions Eq. (5) to yield the incremental stress parameters A/3 and the nodal displacement increments Aq and Aq’. The second is to calculate the instantaneous homogenized tangent modulus E& at the end of the increment in the macroscopic module. From the basic definition in Eq. (291, components of E&, are obtained by averaging the true stress increments in response to unit increments in components of the strain tensor G. Details of this procedure are given in [7].
(31)
3.3. Incorporation in the macroscopic analysis mod-
where S,7klis an instantaneous elastic-plastic compliance tensor. The last term in Eq. (31) incorporates the effect of macroscopic strains in the microstructure. The microscopic VCFEM module is executed for two purposes in each increment of the macro-
The Voronoi cell finite element module is incorporated in a macroscopic analysis module with the interface being created by the homogenization procedure. The general purpose commercial code
+ j ( Zij + AZij)Aui; Y,
dY,
ule
(b) i
Fig. 15.(a) Macroscopic ABAQUS model for connecting rod having regions with three different microstructures, microstructure in region A, (c) hexagonal packing microstructure in region B.
(b) square edge packing
K. Lee, S. Ghosh/Computationnl
ABAQUS is chosen to serve as the macroscopic analysis program. The material constitutive relation at each integration point of ABAQUS elements is input through the homogenization process by using results from the microscopic VCEEM analysis. This
Materials Science 7 (1996) 131-146
143
interface between ABAQUS and VCFJZM is created through UMAT in ABAQUS. The analysis code resulting from this macro-micro coupling is termed as VCFEM-HOMO. Within each iteration loop of macroscopic analysis in the p-th increment, micro-
1.438E-01
1.171E-01
9.038E-02
(4 8.382E-02
3.887E-02
1.013E-O2
(d) Fig. 16. Von Mises stress distribution (a) macroscopic stress for homogeneous microstructure, (b) macroscopic stress for heterogeneous microstructure, (c) microscopic stress within a RME at global point A, (d) microscopic stress within a RME at global point B.
144
K. Lee, S. Ghosh/Compurutioncrl
scopic state variables are computed in VCFEM using given values of the macroscopic strains at the beginning of the increment, as shown in the flow chart of Fig. 2.
Moterials
Scrence 7 (19961 131-146
4. Multiple scale analysis examples The homogenization module VCFEM-HOMO is used to solve two multiple scale problems with
1.136E-02
9.066Eo3
6.616E-03
(a) 4.544E-03
2.272E-03
(c)
(4
Fig. 17. Effective plastic strain distribution (a) macroscopic strain for homogeneous microstructure, (b) macroscopic strain for heterogeneous microsuuchxe. (c) microscopic strain within a RME at a global point A, Cd) microscopic strain within a RME at global point B.
K. Lee, S. Ghosh/CompututionalMaterial.s
various microstructural morphologies. In the first problem, the macroscopic structure is a thick plate with an uniform array of large circular holes in plane strain uniaxial tension. Only a portion of the plate, shown in Fig. 5a, is analyzed from symmetry considerations. Effects of various microstructures on the overall behavior of the heterogeneous material is investigated. In particular, the RME’s considered are (a> square edge packing with a circular heterogeneity, (b) square diagonal packing with a circular heterogeneity, cc> random packing with 15 circular heterogeneities and (d) random packing with 15 elliptical heterogeneities, shown in Fig. 5b-e. The volume fraction of heterogeneities is assumed to be 40%. The material is assumed to be a IT/AL composite and aluminum with voids. The boron inclusion is assumed to be elastic, while the aluminum matrix is an elastic-plastic material with the following properties. FT Jiber: Young’s modulus (EC): 344.5 GPa; Poisson ratio (v,): 0.26; Aluminum matrix: Young’s modulus (E,): 68.9 GPa; Poisson ratio (v,,,): 0.32; initial yield stress (Y,): 94 MPa; post yield flow rule: Eequ= W&n~a,,“/W The macroscopic ABAQUS model consists of 128 QUAD4 elements. The left and right edges are constrained to move in vertical straight lines, and the top and bottom edges are pulled to an overall strain of 0.5%. Figs. 6-9 show the effective stress contours plots at the structural level and also in a microstructural RME at the comer point A. It is observed that the maximum effective stress generally occurs at this point. There is a narrow ligament between two large holes along which strains localize as the deformation intensifies. Similarity in the patterns of macroscopic stress distribution is explained from an observation that the macroscopic state of stress at each point is essentially uniaxial. Thus anisotropy emanating from the microstructural morphology does not play an important role. Magnitude of stresses are however considerably different depending on the microstructural arrangement. This is further evidenced in a macroscopic stress-strain plot (Fig. lOa at the point A. The square edge and square diagonal distributions yield nearly identical response. The random void distributions exhibit significantly more ductile behavior compared to the regular distributions. True
Science 7 (1996) 131-146
145
stress in the microstructure is significantly higher than the macroscopic stresses. For example, at point A the maximum microscopic effective stress is 160% higher for square edge, 175% higher for square diagonal, 265% higher for random circular and 340% higher for random elliptical packing. The same problem is considered again with the IT/AL composite microstructure. The same microstructural morphologies are used with 40% fiber volume fraction. Figs. 1l-14 show the effective stress contour plots at the end of loading. The macroscopic plots show that the localization ligament between two macroscopic holes is much less severe and more diffused compared to the porous material. For the same level of macroscopic strain, the square edge packing yields the highest values of effective stress while the square diagonal packing yields the lowest, with the random packings in between. This is also evidenced in the macroscopic stress-strain plot at point A in Fig. lob. The microscopic contour plots reveal that the true stress in the microstructure is significantly higher than the macroscopic counterparts. For example, at point A the maximum microscopic effective stress is 91% higher for square edge, 50% higher for square diagonal, 180% higher for random circular and 177% higher for random elliptical packing. This shows that high stress concentration occurs in certain fibers for random microstructures, and are consequently more susceptible to microstructural damage by fiber cracking or debonding for identical levels of macroscopic stresses. The second problem involves a connecting rod with different microstructures at different regions as shown in Fig. 15. Region A consists of a square edge packing representative material element (RME), region B has a hexagonal packing RME, while the remainder of the connecting rod has a homogeneous aluminum microstructure. The material composition for the heterogeneous regions A and B is 20% IT/AL composite, with same properties as in the first set of examples. Boundary conditions on the top half of the connecting rod being modeled are shown in Fig. 15a. The rod is assumed to be loaded by a linearly varying internal pressure, applied on one half of the inner race as shown in Fig. 15a, representing contact with the crank shaft. At the end of 10 uniform increments, the pressure variation stands at zero at the top to a maximum of 0.075 GPa at the
146
K. Lre. S. Ghosh / Compurcrtiond
center. For comparison, the same problem is also executed with a purely homogeneous, aluminum microstructure. 146 QUAD4 elements are used in ABAQUS for macroscopic model, and only one Voronoi cell element is required for microscopic evolution. Fig. 16a and b show the distribution of effective macroscopic Von Mises stress for the homogeneous and heterogeneous rod respectively, and Fig. 16c and d depict the corresponding microstructural stresses at macroscopic points A and B in Fig. 16b. The corresponding effective plastic strains are depicted in Fig. 17a-d. It is observed that the effective macroscopic stress is relatively unaltered by the change from homogeneous to heterogeneous microstructure, though the microstructural stress is significantly higher especially for the square edge packing. However, the microstructure alteration has a marked effect on the effective macroscopic plastic strain. The maximum value of plastic strain reduces from 0.011361 in the homogeneous problem to 0.008217 at the same point A, for the heterogeneous case. However, very high values of plastic strain are observed at the matrix-inclusion interface in the microstructure, signalling potential sites for damage initiation.
5. Conclusions In this paper, a multiple-scale computational tool (VCFEM-HOMO) is devised for performing small deformation elastic-plastic analysis of heterogeneous materials with inclusions and voids in the microstructure. The microscopic analysis is conducted with the Voronoi cell finite element model while a conventional displacement based FEM code (ABAQUS in this paper) executes the macroscopic analysis. Coupling between the scales is accomplished through the user based UMAT window in ABAQUS by asymptotic homogenization. The accuracy and efficiency of VCFEM are established by comparing with conventional FEM commercial packages. For a wide range of problems VCFEM delivers very high accuracy at a considerably low computational effort. The effect of various microscopic arrangements on the mechanical response at the two scales are investigated. For complex microstructures, the efficiency of microstructural VCFEM makes it possible to realize exhaustive
Muterids
Scirnce 7 (1996) 131- 146
multi-scale computations for these complex microstructures in comparison with conventional FEM. In conclusion, the Voronoi cell finite element model with asymptotic homogenization (VCFEM-HOMO) emerges as an important tool for analyzing arbitrary microstructures in many materials. It is easily adapted with commercial packages at the structural scale which makes it very attractive.
Acknowledgements Support of this work by the Unites States Army Research Office through grant No. DAAL03-91-G0168 (Program Director: Dr. K.R. Iyer), and by the National Science Foundation through grant No. MSS-9301807 (Mechanics and Materials) is gratefully acknowledged. Computer support by the Ohio Supercomputer Center through grant No. PAS813-2 is also gratefully acknowledged.
References [I] S. Ghosh and S. Moorthy, Comput. Methods Appl. Mech. Eng. 121 (1995) 373-409. [2] S. Moorthy and S. Ghosh, Int. J. Numer. Methods Eng. 39 (1996) 2363-2398. [3] S. Ghosh, K. Lee and S. Moorthy, Int. J. Solids Struct. 32(l) ( 1995) 27-62. [4] S. Ghosh and Y. Liu, Int. J. Numer. Methods Eng. 38(8) (1995) 1361-1398. [5] S. Ghosh and S.N. Mukhopadhyay, Comput. Struct. 41 (1991) 245-256. [61 ABAQUS/Standard User’s Manual, Version 5.4 (Hibbit Karlsson and Sorensen Inc., 1994). [71 S. Ghosh, K. Lee and S. Moorthy, Comput. Methods Appl. Mech. Eng. I32 (1996) 63- 116. [8] T.H.H. Pian. AIAA J. 2 (1964) 1333-1336. [9] S.N. Athni, Int. J. Solids Struct. 9 (1973) 1177-l 191. [IO] S.N. Atluri, J. Stmct. Mech. 8(l) (1980) 61-92. 11 I] J.M. Guedes and N. Kikuchi, Comput. Methods Appl. Mech. Eng. 83 (1991) 143-198. 1121 J.M. Guedes, Ph.D. Dissertation (University of Michigan, MI, 1990). (131 C-H. Cheng. Ph.D. Dissertation (University of Michigan, MI, 1992X [I41 V.Z. Parton and B.A. Kudryavtsev, Engineering Mechanics of Composite Structure (CRC Press, 1993). [I51 S.J. Hollister and N. Kikuchi, Comput. Mech. 10 (1992) 73-95. [I61 F. Devries, H. Dumontet, G. Duvant and F. Lene. Int. J. Numer. Methods Eng. 27 (1989) 285-298. (171 P. Suquet, in: Homogenization Techniques for Composite Media, eds. E. Sanchez-Palencia and A. Zaoui (SpringerVerlag, 1987) pp. 194-278.
Computational
Materials Science 7 (1996) 13 I-
146
Small deformation multi-scale analysis of heterogeneous materials with the Voronoi cell finite element model and homogenization theory Kyunghoon Lee Applied
Mechanics
Section, Department
ofAerospace Engineering,
‘,
Somnath Ghosh Applied
Mechanics
* ,*
and Auiution,
The Ohio State University,
Columbus,
OH 43210, USA
Abstract In this paper, a multiple scale finite element model (VCFEM-HOMO) has been developed for elastic-plastic analysis of heterogeneous (porous and composite) materials by combining asymptotic homogenization theory with the Voronoi cell finite element model (VCFEM). VCFEM for microstructural modeling originates from representative material elements at sampling points in the structure. Structural modeling is done by the general purpose finite element code ABAQUS, and interfacing with the microscale VCFBM analysis is done through the user subroutine in ABAQUS for material constitutive relation, UMAT. Asymptotic homogenization in UMAT generates macroscopic material parameters for ABAQUS. Keywords:
Voronoi cell finite element method; Heterogeneous materials; Asymptotic homogenization
1. Introduction The last three decades we have experienced a surge in the advancement of science and technology for heterogeneous materials that contain dispersions of multiple phases in the microstructure. Examples are metal/alloy systems with second phase in the form of precipitates and pores, and composite materials containing a dispersion of fibers, whiskers or particulates in the matrix. Rigorous fundamental studies are needed for understanding deformation and failure mechanisms prior to design of advanced materials in high performance applications. These models establish the effect of second phase shapes,
* Corresponding
author. Tel.: + 1-614-2922599; e-mail: [email protected]. ’ Graduate Research Associate. * Associate Professor.
2927369;
0927-0256/96/$15.00 PII
Copyright
SO927-0256(96)00072-9
fax: + I -614-
sizes and distributions on evolving state variables like stresses, plastic strains, void initiation and growth, and evolving material variables like strain hardening and flow stress. Ghosh and coworkers have developed a material based Voronoi cell finite element model (VCFEM) [l-5] in an attempt to overcome difficulties in modeling arbitrary microstructures by conventional finite element methods. The VCFEM mesh evolves naturally from a heterogeneous microstructural material element (RME) by Dirichlet tessellation into a network of multi-sided ‘Voronoi’ polygons containing one second phase heterogeneity each, at most. The multi-phase Voronoi polygons identified with the material’s basic structural elements, constitute elements in VCFEM. VCFEM with asymptotic homogenization for elastic problems have been presented in Ref. [3]. This paper presents a coupled multiple scale computational model for heterogeneous elastic-plas-
0 1996 Elsevier Science B.V. All rights reserved.
132
K. Lee. S. Ghosh / Computationul
tic structures. Only two dimensional problems are considered. Microstructural analysis for various different representative material element arrangements, is done with the Voronoi cell finite element model. The commercial general purpose code ABAQUS [6] is used for global analysis at the level of overall structural geometry and applied loads. The interfacing between macro- and micro-calculations are done through the user subroutine window UMAT in ABAQUS. Numerical examples are conducted to investigate the advantages of coupled multiple scale analysis over other unit cell and continuum based theories.
Muterids
Science 7 (1996) 131-146
element boundary. Also let Au correspond to an equilibrated stress increment in Y,, Au to a compatible displacement increment on aY,, and At to a traction increment on the traction boundary &,. The incremental problem is solved by using a two field assumed stress hybrid variational principle, derived from an element energy functional as: &(Ao, =-
AB(cr,Aa)dY-je:AodY / Y, r, +
/ AY(
a+Au).n’.(u+Au)aY
-,r’(i+ai)(u+Au)dl In?
2. Microstructural analysis with VCFEM Voronoi cells, resulting from Dirichlet tessellation of a heterogeneous microstructure, make rather unconventional elements, due to the arbitrariness in the number of edges. The Voronoi cell finite element model developed by Ghosh and coworkers [l-4,7] avoids difficulties of conventional displacement based FEM formulations by invoking the assumed stress hybrid method, introduced by Pian [8]. In this formulation, independent assumptions are made on an equilibrated stress field in the interior of each element and a compatible displacement field on the element boundary. Small deformation elastic-plastic analysis of materials with embedded second phase has shown significant promise with respect to efficiency and accuracy. It is based on a hybrid formulation originally proposed by Atluri and coworkers [9,10]). Details of this development are presented in Refs. [ 1,2]. Consider a typical representative material element (RME) Y tessellated into N Voronoi cells, as shown in Fig. lb. This is based on the location, shape and size of N heterogeneities as explained in [5]. The matrix phase in each Voronoi cell Y, is denoted by Y, and the heterogeneity (void or inclusion) is denoted by Y,. The matrix-heterogeneity interface aY, has an outward normal nc, while ne is the outward normal to the element boundary JY,. An incremental finite element formulation is invoked to account for rate independent plasticity. At the beginning of the p-th increment, let CTbe an equilibrated stress field with a strain field e( 0, loud history), and u be a compatible displacement field on the
Au)
- / av,(
urn+,,“‘-u’-Au’)
nc . (u’ + Au’) c?Y,
(1)
where Au’ is the displacement of the interface and AB is the increment in element complimentary energy. Superscripts m and c represent respectively the matrix and second phase parts of the Voronoi cell element. The energy functional for the entire domain is obtained by adding each element functional as fl= Cr= , Lr,. The first variation of fl, with respect to the stress increments A u, results in the kinematic relations as the Euler equation, while the first variation of n with respect to boundary displacement increments Au yields traction conditions as Euler equations. Independent assumptions on stress increments A u are made in the matrix and heterogeneity phases to accommodate stress jumps across the interface. This can be done using stress functions @(x, y) for the matrix and inclusion phases, resulting in stress expressions as: {Au”)
= [P”‘(L
y)]IApm}
{Au’] = [Pc( x, Y)]{W)
t
(2)
where (Ap)‘s correspond to a set of yet undetermined stress coefficients and [PI is a matrix of interpolation functions. Compatible displacement increments on the element boundary aY, as well as on the interface aY,, are generated by interpolation in terms of generalized nodal values. The displacement
K. Lee, S. Ghosh / Computatiod
Marerials Science 7 (1996) I3/-146
(b)
133
(cl
Fig. I. A heterogeneous structure with various levels. (a) The global structure and computational geometry, (b) representative material elements (RME) at a point and the corresponding VCFJ3 model and, (cl a basic structural element (BSE) represented by a Voronoi cell.
increments on the element boundary and interface may then be written as, {Au] = [L’](Aq] 9
(Au’}
= [L’](Aq’},
(3)
where {Aq} and (Aq’) are generalized displacement increment vectors and [L] is an interpolation matrix. Substituting element approximations for stresses Eq. (2) and displacements Eq. (3), in the energy functional Eq. (l), and setting the first variations with respect to the stress coefficients A p m and A /3’ respectively to zero, results in the following two weak forms of the kinematic relations, /,
PITI
:,
[Pm]T[n’][L’]
0
I? c= I
I,
[L’f[nc]T[Pc]dY c
/3” +Aj3”
,f3’+ Ap’
For an elastic-plastic material, the strain increments be in Eq. (4) are non-linear functions of the current state of stress u as well as of their increments A u. The non-linear finite element Eq. (4) and Eq. (5) are solved for the stress parameters (A p m, A /3“) and the nodal displacement increments (Aq, Aq’) in the p-th increment.
e+Ae}dY
m
= /
the weak form of the traction reciprocity conditions as
dY{Aq}
[Pm]T[nC][LC]dY{ Aq’), arc
2.1. Shape based stress representations /r,[PclT{e + Ae} dY = i, <
[Pc]T[nC][L’]dY{Aq’). c (4)
Setting the first variation of the total energy functional with respect to Aq and Aq’ to zero, results in
An important criterion, affecting the convergence of multiple phase Voronoi cell elements, is the proper choice of stress fields in each of the constituent phases. Choosing stress functions from micromechanics considerations, adds considerably to the ele-
K. Lee, S. Ghosh/ Computational Materids
134
Science 7 (19%)
131-146
Initial Global &Local Information
p=l
UMAT (element #l)
sm
/ Minocopic VCFEM-..
Microcopic VCFEM-..
IPI
El II
~
L-
p=p+l
Fig.2. Flowthan
ment efficiency. Three different conditions that are indispensable in this regard are: 1. Stress functions should, in some way, account for the shape of the heterogeneity. 2. Effects of the heterogeneity shape should vanish
of
microscopic analysis.
at large distances from the interface for matrix stress functions. 3. Shape effects in matrix stress functions should facilitate traction reciprocity at the interface. Pure polynomial forms Airy’s stress functions do
‘7 (4
Fig. 3. (a) ANSYS and (b) VCFEM meshes for randomly
r*
packed RME with circular voids CV, = 20%).
K. Lee, S. Ghosh/ComputationulMaterialu
not explicitly account for the shape of the heterogeneity and require very high order terms for convergence. To deal with arbitrary shaped interfaces in VCFEM, a special procedure is pursued. Suppose that equation of the interface aY, in Fig. lc can be expressed in polar coordinates as ,g(r, 0) = 0, where the r coordinate is measured from the centroid of the heterogeneity. A Fourier series expansion for r in terms of the polar angle 13may be expressed as: r=a,+
Ca,cos(nB) n
+ Cb,sin(&) n
Science 7 (1996) 131-146
135
6.0934
R,, = 1.0
on JY,, (6)
where a, and b, are the Fourier coefficients. The interface equation may then be expressed from Eq. (6) as, g( r, 0) =f-
;
- c %os( n a0
- F :sin(no)
ne)
=O.
(4
(7)
0
Here f corresponds to a function that transforms any arbitrary shaped interface to an approximate unit circle, since fl X, y) = 1 on JY,, with the property that l/f+ 0 as (x, y) + x. This function is now used to construct stress functions in the matrix region Y, from the two part expression Grn = @pomly + @e”,,
-
Fig. 5. (a) ABAQUS model for thick various RME’s at 40% volume fraction, 6) square diagonal packing, (d) random heterogeneities and (e) random packing geneities.
(8)
expressed in terms of a polynomial expanwith @pOmly (a) 0.050
plate with circular holes, (b) square edge packing, packing with 15 circular with 15 elliptical hetero-
ANSYS
.._.._._... ,,CFEM 8 =
0.040
-Ti -”
v) 2 e 2 .Y
0.030
I
0.0x
,
-
0.01c
,
-
WVCFEM D-OANSYS
a 2 H0
O.oOC1 0.0
Fig. 4. (a) Macroscopic x/L= 0.5.
0.2
Macroscopic
stress-strain
0.4
0.6
0.8
-0.02
’ 0.0
1 0.2
tensile strain(%)
response
and (b) Microscopic
0.4
0.6
0.8
1.0
Y/L
stress distribution
at 0.8% strain for random
packed
RME along
K. Lee, S. Ghosh / Cmtpututiortd
136
Matrriuk Scxence 7 (1996) 131-146
(b)
(4 Fig. 6. Von Mises stress distribution
for porous material with square edge packing (a) macroscopic
sion of coordinates, i.e. GP’& = ~,,,y~py~pyy. For each polynomial term in @& there exist shape based reciprocal terms cP~~, that give rise to stresses equilibrating the traction field on a Y,. This is written as:
=
c xy$fyi, . p,q,i
The inclusion is modeled using polynomial stress functions only. A detailed treatment of stress func-
stress, (b) microscopic
tion selection and convergence Ref. [2]. 2.2. Constitutive
stress at point A.
issues is presented in
relations
A rate independent small deformation elastic-plastic constitutive relation, following J, flow theory with isotropic hardening is considered in this paper. An additive decomposition of the strain increments Ae into an elastic part Ae’ and a plastic part her’ is assumed. The yield surface in stress-space at the beginning of the p-th increment is expressed as, Y”(e,
Ae) = \/$v’:u’,
(9)
2.472EG
hln.J+
(b)
(a) Fig. 7. Von Mises stress distribution A.
for porous material with square diagonal packing (a) macroscopic
strain, (b) microscopic
strain at point
K. Lee, S. Ghosh / Compurutional Muterids
where g’ is the deviatoric stress and Y(e, Ae) is the radius of the flow surface. The plastic strain increment Aep’ is obtained by numerically integrating the flow rule by the backward Euler method to yield: AeP’=Ah(a+Aa)‘,
(10) where AA is a non-negative incremental flow parameter. Since Ae is in general a function of AA, u and A g, the flow surface radius can be expressed in the form, YP+’ = Yp+‘(AA, u, Au) and the flow parameter AA can be evaluated from the following relations: A A = 0 if u$+ I < YP (elastic unloading) P+ ffeff
I =
yP+
I
if
(+P+ eff
I > -
,
(11)
YP
(neutral and plastic loading) .
(12)
Details on numerical implementation are provided in
[1,21. 2.3. A numerical example with VCFEM The effectiveness of the elastic-plastic Voronoi cell finite element model in analyzing heterogeneous microstructures is established by several numerical examples in [1,2]. In this example, the representative material element (RME) consists of 29 randomly located circular voids with volume fraction V, = 20%. Results generated by VCFEM are compared with those of a displacement based commercial finite
Fig. 8. Von Mises stress distribution A.
Science 7 (1996) 131-146
137
element code ANSYS with a very high resolution mesh necessary for convergence. Fig. 3b and a show the VCFEM mesh with 29 elements and the ANSYS mesh with 5282 QUAD4 elements, respectively. Plane strain uniaxial tension loading is considered to a maximum of 0.8% macroscopic strain. Periodic strain uniaxial tension loading is considered to a maximum of 0.8% macroscopic strain. Periodic construct of the RME is achieved through repeatability conditions on the traction free face. The matrix material is Aluminum with the following properties: E = 69 GPa (Young’s modulus); v= 0.33 (Poisson’s ratio); Y, = 43 MPa (initial yield stress); o& = Y + e”.‘25 (post yield hardening law). Compzson of macroscopic response in Fig. 4a clearly establishes VCFEM as an accurate method for modeling overall behavior of random microstructures at a significant advantage in computing efficiency. The microscopic stress distribution through the section A-A in Fig. 3b at 0.8% strain is depicted in Fig. 4b. Once again the comparison is very satisfactory, with VCFEM producing similar patterns and peak stresses with the highly refined ANSYS model. 3. Asymptotic
homogenization
Consider a heterogeneous body occupying a rein Fig. la, for which the microstrucgion oSrlUCtUTe ture constitutes of spatially periodic representative material elements @ME’s) as shown in Fig. l(b). In real heterogeneous materials, dimensions of the RME
for porous material with random circular packing (a) macroscopic
stress, (b) microscopic
stress at point
K. Lee, S. Ghosh/
138
Computdod
Mutcriu1.v
of characteristic length 1 are typically very small in comparison with the dimensions of the body of characteristic length L. The ratio of these microscopic and macroscopic scales I/L is represented by a very small positive number E. A high level of heterogeneity in the microstructure causes a rapid variation of evolutionary variables e.g. deformation and stresses in a small neighborhood E of the macroscopic point x. This corresponds to a microscopic scale X/E and consequently, all variables are assumed to exhibit dependence on both length scales i.e. @‘= @(x, X/E). The superscript E denotes association of the function with the two length scales. In this notation, 0’ corresponds to a connected domain that extends the structural domain to its microstructure. In most of the work on homogenization theory [l l- 161, a periodic repetition of the microstructure about a macroscopic point x has been assumed, thereby making the dependence of the function on y (= X/E), periodic. This characteristic is often termed as Y-periodicity. where Y corresponds to a RME. For small deformation elasto-plasticity, the rate or incremental forms of the equilibrium equation, kinematic relation, and constitutive relations are given as:
Science
7 (19961
131-146
where tif = tif(x, y) is a Y-periodic rate of displacement field in y. Furthermore the boundary conditions are assumed to satisfy the following equations on the prescribed traction and displacement boundaries, respectively. ci,fnj = ii on r,,
(‘6)
ST=Gi on c,,
(‘7) where n is the unit normal to the boundary. In homogenization theory, the Y-periodic displacement rate or increment field is approximated by an asymptotic expansion with respect to parameter E: C’(x) = tiO(x, y) + Eil’(X, y) + EW(X, y) + . . . . y=;
(18)
The stress rate tensor c+~;.S can then be expressed as I (i,; = -&;;.9 + &;.f + E&j + E*&;; . . . ,
(19)
E
where
(13) in a’, ki;.5= ETjk, S;,
in L?‘,
A Fig. 9. Van point A.
Mises stress distribution for
(‘4)
ti,;: = EFjk,
(15)
Putting the expansion of &;if Eq. (19) in the rate
(b)
(8) porous
material with
(20)
random elliptical packing (a) macroscopic strain, (b)
microscopicswainat
K. Lee, S. Ghosh/C.ompututionalMaterials
form of equilibrium Eq. (13), and setting each coefficientofei,(i=--1,0,1,2 ,... )tozero,resultsin the following set of equations:
Science 7 (19%) 131-146
a Y-periodic function representing characteristic modes of the RME. Substituting Eq. (26) in Eq. (20) yields the microscopic constitutive relations as:
adi; =o, aYj
t+;‘(y)
a&.‘-.‘. ac$ I’+ =o, axj aYj
ax,"' = E,;pm T;; + aY, [
1
(28)
’
where Tit’ is a fourth order identity tensor. The mean of Eq. i 28) yields the homogenized elastic-plastic
ah.2 ski)
1’+ aYj
139
x+f,=o.
(21)
I
Eq. (21) and Eq. (20) lead to the trivial value for &$, and therefore establishes that rio is only a function of x as shown in [ 161, i.e. ki;.s= 0 and tip = L$‘(x) .
(4 0.10
(22)
Substituting Eq. (22) into the second of Eq. (21) leads to the Y-domain equilibrium equation
a&.'.
--‘/ -=
(23)
0.
aYj
From Eq. (15) and Eq. (191, by neglecting the terms associated with E or higher, the constitutive relation in Y is expressed as
random elllpse(40% -
square edge@o%
void) void)
0 - - * square dlagona&?C% e--a
(24)
Macroscopic
random
circle(20%
void) void)
strain(%)
where
ati: f-. JY,
(b) 0.25
(25) 0
Here t;, is the local or microstructural strain rate tensor, for which Sk, = &i/ax,) is an averaged macroscopic part, and k& = &$/ay, is denoted as a fluctuating strain rate tensor (Suquet [17]). Due to linearity of the rate (incremental) problem, di;.f and I$ can be expressed in the forms
b
0.20
1 B .;
0.15
f s
0.10
8 P 0.05
L I
where
0.00
0.0
aeig'(y)
-
= 0 (microscopic equilibrium).
0.1
0.2 Macroscopic
0.3
0.4
0.5
strain(%)
(27)
aYj
In Eq. (26), eiijk’is a Y-periodic function and x/’ is
Fig. 10. Evolution of effective stress with evolving strain at point A in the heterogeneous plate with 40% second phase volume fraction for (a) porous material (b) composite material.
K. Lee. S. Ghosh / Compututional
140
Fig. 11. Von Mises stress distribution point A.
for composite
material
Muterids
Science 7 (19961 131-146
with square edge packing
tangent modulus, for use in the macroscopic analysis, in the form 1 E,& = ( $;‘> = IyI ,c?;;’ dY /
aXj
3.1. Macroscopic
stress, (b) microscopic
stress at
tions for an elastic-plastic problem heterogeneous material as: aiij -=
dY.
(a) macroscopic
-.i
in
%ructure
T
(29)
equations
Sijnj = ii on r,,
The mean of the third equation in Eq. (21) on Y yields the macroscopic form of the governing equa-
lip=Gi
onr,,
(30)
2 OWE41
172lE-01
1,UlEO)
b.UMWQ
(4 Fig. 12. Von Mises stress distribution point A.
for composite
material with square diagonal
packing (a) macroscopic
strain, (b) microscopic
strain at
K. Lee, S. Ghosh/Computationd
Materials Science 7 (1996) 131-146
141
Fig. 13. Von Mises stress distribution for composite material with random circular packing (a) macroscopic stress, (b) microscopic stress at point A.
corresponds to the global domain where%mCtUre with f, and r, as traction and displacement boundaries respectively, % = (k ‘) and ire are the averaged macroscopic stress rate tensor and the displacement rate respectively. Small deformation analysis is pursued at the macroscopic level with the Eq. (301, and with the Voronoi cell finite element model (VCFEM) for the microscopic level.
3.2. Elastic-plastic homogenization with VCFEM In an incremental formulation, the equilibrated microscopic stress increment corresponds to A u ‘( = A u “1 in Eq. (26) and the microstructural strain increments are designated as Ae’ in Eq. (25). Similarly, the increments in microscopic displacements on the cell boundaries aY, are identified with
Fig. 14. Von Mises stress distribution for composite material with random elliptical packing (a) macroscopic strain, (b) microscopic strain at point A.
K. Lee, S. Ghosh/Compututiod
142
Au’
in Eq. (26) and those on the interface are denoted by Au”. In the absence of traction boundaries due to periodicity conditions, the incremental energy functional for each Voronoi cell element in Eq. (1) is modified for the homogenization process as:
-
j( a*‘
uifm + Api;” - c_$’ - Au;;~)
x(u;‘+Auj’)n;dc?Y
Materids
Science 7 (1996) 131-146
scopic module. The first is to evaluate the microscopic stress increments A u ’ from given values of the macroscopic strain G at the beginning of the step, and its increment AZ. This step involves the iterative solution of modified forms of the kinematic relations Eq. (4) and the traction reciprocity conditions Eq. (5) to yield the incremental stress parameters A/3 and the nodal displacement increments Aq and Aq’. The second is to calculate the instantaneous homogenized tangent modulus E& at the end of the increment in the macroscopic module. From the basic definition in Eq. (291, components of E&, are obtained by averaging the true stress increments in response to unit increments in components of the strain tensor G. Details of this procedure are given in [7].
(31)
3.3. Incorporation in the macroscopic analysis mod-
where S,7klis an instantaneous elastic-plastic compliance tensor. The last term in Eq. (31) incorporates the effect of macroscopic strains in the microstructure. The microscopic VCFEM module is executed for two purposes in each increment of the macro-
The Voronoi cell finite element module is incorporated in a macroscopic analysis module with the interface being created by the homogenization procedure. The general purpose commercial code
+ j ( Zij + AZij)Aui; Y,
dY,
ule
(b) i
Fig. 15.(a) Macroscopic ABAQUS model for connecting rod having regions with three different microstructures, microstructure in region A, (c) hexagonal packing microstructure in region B.
(b) square edge packing
K. Lee, S. Ghosh/Computationnl
ABAQUS is chosen to serve as the macroscopic analysis program. The material constitutive relation at each integration point of ABAQUS elements is input through the homogenization process by using results from the microscopic VCEEM analysis. This
Materials Science 7 (1996) 131-146
143
interface between ABAQUS and VCFJZM is created through UMAT in ABAQUS. The analysis code resulting from this macro-micro coupling is termed as VCFEM-HOMO. Within each iteration loop of macroscopic analysis in the p-th increment, micro-
1.438E-01
1.171E-01
9.038E-02
(4 8.382E-02
3.887E-02
1.013E-O2
(d) Fig. 16. Von Mises stress distribution (a) macroscopic stress for homogeneous microstructure, (b) macroscopic stress for heterogeneous microstructure, (c) microscopic stress within a RME at global point A, (d) microscopic stress within a RME at global point B.
144
K. Lee, S. Ghosh/Compurutioncrl
scopic state variables are computed in VCFEM using given values of the macroscopic strains at the beginning of the increment, as shown in the flow chart of Fig. 2.
Moterials
Scrence 7 (19961 131-146
4. Multiple scale analysis examples The homogenization module VCFEM-HOMO is used to solve two multiple scale problems with
1.136E-02
9.066Eo3
6.616E-03
(a) 4.544E-03
2.272E-03
(c)
(4
Fig. 17. Effective plastic strain distribution (a) macroscopic strain for homogeneous microstructure, (b) macroscopic strain for heterogeneous microsuuchxe. (c) microscopic strain within a RME at a global point A, Cd) microscopic strain within a RME at global point B.
K. Lee, S. Ghosh/CompututionalMaterial.s
various microstructural morphologies. In the first problem, the macroscopic structure is a thick plate with an uniform array of large circular holes in plane strain uniaxial tension. Only a portion of the plate, shown in Fig. 5a, is analyzed from symmetry considerations. Effects of various microstructures on the overall behavior of the heterogeneous material is investigated. In particular, the RME’s considered are (a> square edge packing with a circular heterogeneity, (b) square diagonal packing with a circular heterogeneity, cc> random packing with 15 circular heterogeneities and (d) random packing with 15 elliptical heterogeneities, shown in Fig. 5b-e. The volume fraction of heterogeneities is assumed to be 40%. The material is assumed to be a IT/AL composite and aluminum with voids. The boron inclusion is assumed to be elastic, while the aluminum matrix is an elastic-plastic material with the following properties. FT Jiber: Young’s modulus (EC): 344.5 GPa; Poisson ratio (v,): 0.26; Aluminum matrix: Young’s modulus (E,): 68.9 GPa; Poisson ratio (v,,,): 0.32; initial yield stress (Y,): 94 MPa; post yield flow rule: Eequ= W&n~a,,“/W The macroscopic ABAQUS model consists of 128 QUAD4 elements. The left and right edges are constrained to move in vertical straight lines, and the top and bottom edges are pulled to an overall strain of 0.5%. Figs. 6-9 show the effective stress contours plots at the structural level and also in a microstructural RME at the comer point A. It is observed that the maximum effective stress generally occurs at this point. There is a narrow ligament between two large holes along which strains localize as the deformation intensifies. Similarity in the patterns of macroscopic stress distribution is explained from an observation that the macroscopic state of stress at each point is essentially uniaxial. Thus anisotropy emanating from the microstructural morphology does not play an important role. Magnitude of stresses are however considerably different depending on the microstructural arrangement. This is further evidenced in a macroscopic stress-strain plot (Fig. lOa at the point A. The square edge and square diagonal distributions yield nearly identical response. The random void distributions exhibit significantly more ductile behavior compared to the regular distributions. True
Science 7 (1996) 131-146
145
stress in the microstructure is significantly higher than the macroscopic stresses. For example, at point A the maximum microscopic effective stress is 160% higher for square edge, 175% higher for square diagonal, 265% higher for random circular and 340% higher for random elliptical packing. The same problem is considered again with the IT/AL composite microstructure. The same microstructural morphologies are used with 40% fiber volume fraction. Figs. 1l-14 show the effective stress contour plots at the end of loading. The macroscopic plots show that the localization ligament between two macroscopic holes is much less severe and more diffused compared to the porous material. For the same level of macroscopic strain, the square edge packing yields the highest values of effective stress while the square diagonal packing yields the lowest, with the random packings in between. This is also evidenced in the macroscopic stress-strain plot at point A in Fig. lob. The microscopic contour plots reveal that the true stress in the microstructure is significantly higher than the macroscopic counterparts. For example, at point A the maximum microscopic effective stress is 91% higher for square edge, 50% higher for square diagonal, 180% higher for random circular and 177% higher for random elliptical packing. This shows that high stress concentration occurs in certain fibers for random microstructures, and are consequently more susceptible to microstructural damage by fiber cracking or debonding for identical levels of macroscopic stresses. The second problem involves a connecting rod with different microstructures at different regions as shown in Fig. 15. Region A consists of a square edge packing representative material element (RME), region B has a hexagonal packing RME, while the remainder of the connecting rod has a homogeneous aluminum microstructure. The material composition for the heterogeneous regions A and B is 20% IT/AL composite, with same properties as in the first set of examples. Boundary conditions on the top half of the connecting rod being modeled are shown in Fig. 15a. The rod is assumed to be loaded by a linearly varying internal pressure, applied on one half of the inner race as shown in Fig. 15a, representing contact with the crank shaft. At the end of 10 uniform increments, the pressure variation stands at zero at the top to a maximum of 0.075 GPa at the
146
K. Lre. S. Ghosh / Compurcrtiond
center. For comparison, the same problem is also executed with a purely homogeneous, aluminum microstructure. 146 QUAD4 elements are used in ABAQUS for macroscopic model, and only one Voronoi cell element is required for microscopic evolution. Fig. 16a and b show the distribution of effective macroscopic Von Mises stress for the homogeneous and heterogeneous rod respectively, and Fig. 16c and d depict the corresponding microstructural stresses at macroscopic points A and B in Fig. 16b. The corresponding effective plastic strains are depicted in Fig. 17a-d. It is observed that the effective macroscopic stress is relatively unaltered by the change from homogeneous to heterogeneous microstructure, though the microstructural stress is significantly higher especially for the square edge packing. However, the microstructure alteration has a marked effect on the effective macroscopic plastic strain. The maximum value of plastic strain reduces from 0.011361 in the homogeneous problem to 0.008217 at the same point A, for the heterogeneous case. However, very high values of plastic strain are observed at the matrix-inclusion interface in the microstructure, signalling potential sites for damage initiation.
5. Conclusions In this paper, a multiple-scale computational tool (VCFEM-HOMO) is devised for performing small deformation elastic-plastic analysis of heterogeneous materials with inclusions and voids in the microstructure. The microscopic analysis is conducted with the Voronoi cell finite element model while a conventional displacement based FEM code (ABAQUS in this paper) executes the macroscopic analysis. Coupling between the scales is accomplished through the user based UMAT window in ABAQUS by asymptotic homogenization. The accuracy and efficiency of VCFEM are established by comparing with conventional FEM commercial packages. For a wide range of problems VCFEM delivers very high accuracy at a considerably low computational effort. The effect of various microscopic arrangements on the mechanical response at the two scales are investigated. For complex microstructures, the efficiency of microstructural VCFEM makes it possible to realize exhaustive
Muterids
Scirnce 7 (1996) 131- 146
multi-scale computations for these complex microstructures in comparison with conventional FEM. In conclusion, the Voronoi cell finite element model with asymptotic homogenization (VCFEM-HOMO) emerges as an important tool for analyzing arbitrary microstructures in many materials. It is easily adapted with commercial packages at the structural scale which makes it very attractive.
Acknowledgements Support of this work by the Unites States Army Research Office through grant No. DAAL03-91-G0168 (Program Director: Dr. K.R. Iyer), and by the National Science Foundation through grant No. MSS-9301807 (Mechanics and Materials) is gratefully acknowledged. Computer support by the Ohio Supercomputer Center through grant No. PAS813-2 is also gratefully acknowledged.
References [I] S. Ghosh and S. Moorthy, Comput. Methods Appl. Mech. Eng. 121 (1995) 373-409. [2] S. Moorthy and S. Ghosh, Int. J. Numer. Methods Eng. 39 (1996) 2363-2398. [3] S. Ghosh, K. Lee and S. Moorthy, Int. J. Solids Struct. 32(l) ( 1995) 27-62. [4] S. Ghosh and Y. Liu, Int. J. Numer. Methods Eng. 38(8) (1995) 1361-1398. [5] S. Ghosh and S.N. Mukhopadhyay, Comput. Struct. 41 (1991) 245-256. [61 ABAQUS/Standard User’s Manual, Version 5.4 (Hibbit Karlsson and Sorensen Inc., 1994). [71 S. Ghosh, K. Lee and S. Moorthy, Comput. Methods Appl. Mech. Eng. I32 (1996) 63- 116. [8] T.H.H. Pian. AIAA J. 2 (1964) 1333-1336. [9] S.N. Athni, Int. J. Solids Struct. 9 (1973) 1177-l 191. [IO] S.N. Atluri, J. Stmct. Mech. 8(l) (1980) 61-92. 11 I] J.M. Guedes and N. Kikuchi, Comput. Methods Appl. Mech. Eng. 83 (1991) 143-198. 1121 J.M. Guedes, Ph.D. Dissertation (University of Michigan, MI, 1990). (131 C-H. Cheng. Ph.D. Dissertation (University of Michigan, MI, 1992X [I41 V.Z. Parton and B.A. Kudryavtsev, Engineering Mechanics of Composite Structure (CRC Press, 1993). [I51 S.J. Hollister and N. Kikuchi, Comput. Mech. 10 (1992) 73-95. [I61 F. Devries, H. Dumontet, G. Duvant and F. Lene. Int. J. Numer. Methods Eng. 27 (1989) 285-298. (171 P. Suquet, in: Homogenization Techniques for Composite Media, eds. E. Sanchez-Palencia and A. Zaoui (SpringerVerlag, 1987) pp. 194-278.