Scale effects in plasticity of random media: status and challenges

International Journal of Plasticity 21 (2005) 1119–1160 www.elsevier.com/locate/ijplas

Scale effects in plasticity of random media: status and challenges Martin Ostoja-Starzewski

*

Department of Mechanical Engineering and McGill Institute for Advanced Materials, McGill University, Montre´al, Quebec, Canada H3A 2K6 Received 31 January 2004; in final revised form 3 June 2004

Abstract When the separation of scales in random media does not hold, the representative volume element (RVE) of classical continuum mechanics does not exist in the conventional sense, and various new approaches are needed. This subject is discussed here in the context of plasticity of random, microheterogeneous media. The first principal topic considered is that of hierarchies of mesoscale bounds, set up over a statistical volume element (SVE), for elastic– plastic-hardening microstructures; these bounds, with growing mesoscale, tend to converge to RVE responses. Following a formulation of the said hierarchies from variational principles and their illustration on two specific examples of power-law hardening materials, we turn to rigid-perfectly-plastic materials. The latter are illustrated by simulations in the setting of a planar random chessboard. The second principal topic is the analysis of spatially non-uniform response patterns of randomly heterogeneous plastic materials. We focus here on the geodesic properties of shear-band patterns, and then on the correlation of strain fields to the underlying microstructures. In the case of perfectly-plastic materials, shear-bands become slip-lines, but their spatial disorder is still present, and is described in ensemble sense by wedges of randomly scattered characteristics.
2004 Elsevier Ltd. All rights reserved.

*

Tel.: +1 514 398 7394; fax: +1 514 398 7365. E-mail address: [email protected].

0749-6419/$ - see front matter
2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijplas.2004.06.008

1120

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

1. Introduction 1.1. Separation of scales Plasticity of heterogeneous media has mostly been studied under the assumption of a periodic microstructure. In reality, of course, responses of heterogeneous materials are not periodic, because the shear bands which facilitate the plastic flow tend to find weakest paths in zones of lower yield around the obstacles of higher resistance, and minute spatial fluctuations in material properties facilitate such weakest zones. This then results in a lower effective yield than what is predicted by the periodic model. One cannot even produce a periodic material (either physical or numerical model) that would respond periodically. Considerations of this type are the motivation for dealing with plasticity of random media and associated scale effects. Deterministic continuum mechanics hinges on the concepts of the ‘‘representative volume element’’ (RVE) and the so-called ‘‘separation of scales.’’ Follow the definition of RVE due to Hill (1963), we have: ‘‘a sample that (a) is structurally entirely typical of the whole mixture on average, and (b) contains a sufficient number of inclusions for the apparent overall moduli to be effectively independent of the surface values of traction and displacement, so long as these values are macroscopically uniform.’’ In other words, there are two conditions involved here: (a) the microstructureÕs statistics must be spatially homogeneous and ergodic; (b) the materialÕs effective constitutive response must be the same under uniform boundary conditions of either traction or displacement type. The RVE is supposed to be a spatial domain of size L, satisfying two inequalities:
d< L
Lmacro : ð1:1Þ d
Here d is the microscale (average size of constituents),while Lmacro is the macroscale (i.e. macroscopic body size); L may be called a mesoscale. On the left in (1.1) we indicate that L needs to be either larger or much larger than d, depending on the particular situation at hand. At the same time, L
Lmacro is always required, because, on that basis, L3 may be treated as a dV element of continuum mechanics. Inequalities (1.1) are the statement of the separation of scales, sometimes called the MMM principle in micromechanics (Hashin, 1983). It will now be convenient to introduce a dimensionless parameter characterizing the mesoscale L ð1:2Þ d¼ : d In essence, the RVE is clearly denned in either one of two situations: (i) a unit cell in a periodic microstructure, with d a rather small number (say, 3); (ii) a domain in a random medium, containing a very large (mathematically infinitely large: d ! 1) set of constituents, having statistically homogeneous and ergodic properties.

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

1121

Case (i) is of interest in studies of random media as a basis for comparison and reference, and henceforth will not be studied here. Note that one may also study large systems with periodic geometry under periodic boundary conditions, and aim at the so-called effective medium response. Henceforth, focusing on case (ii), which admits spatial randomness, we must observe that the conditions of homogeneity and ergodicity as well as the strong inequality on the right in (1.1), may not always be satisfied. In such a case one has to work with a mesoscale random field (Fig. 1(b)) which bridges microscale levels (illustrated by (a)) with the macroscale level (c). Some situations where this happens include:  the material displays inhomogeneous statistical fluctuations on macroscopic/continuum length scales (Lmacro), see. e.g. (Clayton and McDowell, 2003, 2004);  the material displays statistical fluctuations on length scales where continuum assumption is conventionally desired but does not hold, e.g. the finite element size is not larger than the RVE, the wavefront thickness is not larger than the RVE;  the material specimen tested in the lab is smaller than the RVE; such is the situation when trying to determine effective continuum properties of a hydroelectric dam made of boulders several meters in size (Huet, 1999), and many other geophysical problems easily come to our mind;  one or more physical dimensions of the material specimen – as, for example, in a micro- (or nano-) electromechanical system (MEMS or NEMS) – is smaller than the RVE (e.g. Ostoja-Starzewski, 2004);  the material specimen to be tested is too expensive to be manufactured to a ÔlargeenoughÕ RVE size.

1.2. Estimation of the random field The need to deal with material spatial randomness forces us to introduce scalarand vector-valued random fields. In the first case we say that E (e.g. YoungÕs modulus) is a (scalar) random field in D (1, 2, or 3) dimensions if it assigns to an elementary event a specific realization over a physical domain X  RD , that is E : X  X ! R, Eðx,xÞ ¼ e,

ð1:3Þ

where x 2 X. In the second case we say that C (e.g. stiffness tensor) is an n-component vector random field in D parameter dimensions if it assigns to an elementary event a realization over X, that is C : X  X ! Rn :

ð1:4Þ

In principle, in order to set up a model of a spatially random medium over a domain X, we need to estimate the random field of a specific material property, e.g., the stiffness tensor C (”Cijkl). This estimation is usually done from one realization C(x) available to us (in a laboratory or field observation), that is

1122

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

Fig. 1. Scales and randomness: (a) sequential alternate random function (modeling a random microstructure); (b) diffusion random function (modeling smoothly inhomogeneous media on mesoscale); (c) macroscopic body.

CðxÞ  lim

L!1

1 L

Z

L

Cðx,xÞdx ¼ 0

Z X

Cðx,xÞdP ðxÞ  hCðxÞi:

ð1:5Þ

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

1123

Here the overbar indicates the volume average, while Æ æ stands for the ensemble average, whereby P(x) is the probability measure (distribution) assigned to the ensemble {C(x,x); x 2 X, x 2 X} The equality (1.5) is the statement of the aforementioned ergodic property, without which one just cannot infer any tangible information from one specimen about the other ones in the ensemble. The lack of ergodicity may, of course, be the case in applications, and one needs to be careful about this fact. In general, the degree to which ergodicity holds depends on the length scale involved: for example, it may hold on moderate length scales, but not on large ones where one may have long range inhomogeneities. Even if the ergodicity holds in a particular situation, the limiting process L ! 1 in (1.5) cannot truly be carried out! Therefore, in practice, the left-hand side in the above is replaced by a spatial average from a finite number M of sampling points m over one x M 1 X CðxÞ  Cðx,xm Þ, ð1:6Þ M m¼1 and the right-hand side by an ensemble average from a finite number N of realizations xn at one sampling point x hCðxÞi 

N 1 X Cðxn ,xÞ: N n¼1

ð1:7Þ

On the mathematical side, the ergodicity of the estimator – i.e., CðxÞ ¼ hCðxÞi – is assured, for sufficiently large M and N, by these properties of the autocorrelation function RC(x)of {C(x,x); x 2 X, x 2 X}: 2

lim RC ðxÞ ¼ hCi ,

jxj!1

ð1:8Þ

Z 1 L RC ðxÞ  limL!1 Cðx,x1 þ xÞCðx,x1 Þdx1 L 0 Z ¼ Cðx,x1 þ xÞCðx,x1 ÞdP ðxÞ  hCðx,x1 þ xÞCðx,x1 Þi X

¼ RC ðxÞ:

ð1:9Þ

1.3. Outline of the paper In this paper we review results and discuss future challenges in plasticity of random heterogeneous media from the perspective of our ongoing studies. Our starting point is the recognition that, when the separation of scales does not hold, the representative volume element (RVE) of classical continuum mechanics does not exist in the sense of Hill condition, and various new techniques are needed. First, we discuss hierarchies of mesoscale bounds for linear elastic microstructures, which, with growing mesoscale, tend to converge to RVE responses such as the macroscopic HookeÕs law. The same type of behavior is then displayed for physically nonlinear elastic and

1124

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

elastic-hardening plastic microstructures, whereby it is pointed out how to modify the Hill condition when dealing the high strain rate phenomena. In particular, we outline the derivation of the said hierarchies from variational principles, and give two particular examples for power-law hardening materials: two-phase composite, and paper. This sets the stage for the original results of this article: hierarchies of mesoscale bounds for rigid-perfectly-plastic materials. These are also illustrated by finite element simulations in the setting of a planar random chessboard. The direct computational accessibility of spatially non-uniform responses of randomly heterogeneous plastic materials prompts us to analyze and exploit their various features. Thus, we look here at the geodesic properties of shear-band patterns and then at the correlation of strain fields to the underlying microstructures. In the latter case, one can infer microscale material properties. In the case of perfectly-plastic materials, shear-bands become slip-lines, but their spatial disorder is still present. Given the stochastic hyperbolic nature of the field equations, this disorder is described in ensemble sense by wedges of randomly scattered characteristics. In particular, as an example of solution of stochastic boundary value problems, we focus here on a plane-strain formulation for materials with the Huber–Mises criterion; our approach is illustrated by a limit analysis of a pipe under internal loading.

2. Hierarchies of mesoscale bounds for elastic and elastic–plastic materials 2.1. Elastic materials 2.1.1. Linear elastic materials In situations lacking the separation of scales, we propose setting up a statistical volume element (SVE), a material domain necessarily smaller than the RVE and definitely finite relative to the microscale. The SVE is then the mesoscale domain alluded to earlier, and its response depends on the boundary conditions. Here the starting point is Z r  e ¼ r  e () ðt  r  nÞ  ðu  e  xÞdS ¼ 0: ð2:1Þ oBd

The equality on the right, called the Hill condition, states that the energetic ðe  rÞ and mechanical ðe  rÞ definitions of response are identical. Clearly, this happens providing the body is loaded according to the equality on the right in (2.1). The latter is satisfied by either one of three types of boundary conditions: (i) uniform displacement (also called kinematic, essential, or Dirichlet) boundary condition (d) uðxÞ ¼ e0  x

8x 2 oBd ;

ð2:2Þ

(ii) uniform traction (also called static, natural, or Neumann) boundary condition (t) tðxÞ ¼ r0  nðxÞ

8x 2 oBd ;

ð2:3Þ

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

1125

(iii) uniform displacement-traction (also called mixed orthogonal) boundary condition (dt) ½uðxÞ  e0  x  ½tðxÞ  r0  nðxÞ ¼ 0 8x 2 oBd :

ð2:4Þ

Suppose for the moment that we deal with a linear elastic microstructure. Then, each of these boundary conditions results in a different stiffness, or compliance tensor. The late Christian Huet (1995, 1999) used the term ÔapparentÕ to distinguish the mesoscale properties from the ÔeffectiveÕ (macroscopic, global, or overall) ones that are typically denoted by the superscript ÔeffÕ. Thus, condition (2.2) yields a mesoscale random stiffness Cdd ðxÞ. On the other hand, condition (2.3) results in a mesoscale random compliance Std ðxÞ. Finally, condition (2.4) results in a mesoscale random stiffness Cdtd ðxÞ. The argument x points to the explicit dependence of mesoscale response on the actual microstructure (recall Fig. 1) – that is, the random nature in the ensemble sense – of the resulting average stress field and of the apparent stiffness tensor, with the fluctuations disappearing in the RVE limit: d ! 1. Hill (1967) and Mandel (1972) gave a qualitative estimate of the error between Std ðxÞ and C dd ðxÞ 3

Std ðxÞ  Cdd ðxÞ ¼ I þ Oð1=dÞ :

ð2:5Þ

Many quantitative estimates of d-dependence – or, what is called finite-size scalingin condensed matter physics – were computed for many different materials by Huet and co-workers, this author and co-workers; see also Terada et al. (2000), and Kanit et al. (2001). Let us now recall hierarchies of mesoscale bounds on macroscopically effective stiffness tensor C1 for linear elastic materials hSt1 i1 6    6 hStd0 i1 6 hStd i1 6    6 C1 6    6 hCdd i 6 hCdd0 i 6    6 hCd1 i Cdd ðxÞ

8d0 ¼ d=2,

ð2:6Þ

Std ðxÞ

where and are defined as explained above. The actual determination of those two tensors necessarily involves some computational method – such as finite elements or spring network discretization (e.g., Ostoja-Starzewski, 2001a,b) – carried out for a single, and finite scale, realization of the microstructure. Thus, both tensors are parametrized by x 2 X, which says there is a random scatter on mesoscale d < 1. As d ! 1, this scatter goes to zero, and we approach the RYE. The brackets Æ æ are dropped for C1, because of the assumed ergodic property of the material. From a stochastic standpoint, given a random two-phase (generally multiphase) composite, we start with a random indicator (or characteristic) field with discretevalued realizations: fvm ðx,xÞ; x 2 X, x 2 R2 g,

ð2:7Þ

and we smooth it by two mesoscale random fields: C : X  R2 ! fIC ðiÞ ,IC ðmÞ g

or

fCðx,xÞ; x 2 X,x 2 R2 g:

ð2:8Þ

1126

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

The point is that we can set up mesoscale domains at any location x, and so we have two approximating continua of bounding character (Ostoja-Starzewski, 1998). Note: the so-called local averaging, such as employed in some stochastic finite element (SFE) models, is not correct. The inequalities on the right (left) of the hierarchy above follow from the principle of minimum potential (respectively, complementary) energy, (Huet, 1990). Thus, in a nutshell, (2.6) is qualitative, and requires some computational mechanics method to be made quantitative. The entire idea of hierarchies of mesoscale bounds can subsequently be – and, indeed, has already been – extended to other than linear elastic materials. 2.1.2. Physically nonlinear elastic materials An elastic–plastic hardening material under monotonically increasing loading may be modeled as a physically nonlinear elastic one. Following up on that well known idea of micromechanics, consider a random medium with the relations (Jiang et al., 2001) r ¼ rðeÞ ¼

owðeÞ ow ðrÞ e ¼ eðrÞ ¼ oe or

ð2:9Þ

holding locally, i.e. pointwise. From this, using variational principles on the mesoscale level, we arrive at a hierarchy of bounds on the strain energy w under displacement (d) boundary conditions wd1 6    6 hwdd i 6 hwdd0 i 6    6 hwd1 i 8d0 ¼ d=2:

ð2:10Þ

On the other hand, we can also obtain a hierarchy of upper bounds on the complementary energy w* t t t w t 8d0 ¼ d=2: 1 6    6 hwd i 6 hwd0 i 6    6 hw1 i

ð2:11Þ

With w and w* being related through the Legendre transformation, the hierarchy (2.11) provides, in effect, a hierarchy of lower bounds on the RVE response modeled by wd1 . In (2.8) wV is the Voigt-type bound for the energy on the smallest scale, while w*R in (2.9) is the Reuss-type bound for the complementary energy. These hierarchies were proved in Jiang et al. (2001) and also illustrated in OstojaStarzewski and Castro (2003) for paper, see Figs. 2 and 3. The motivation for the latter application was twofold: (i) paper displays random formation, which is seen as greyscale spatial fluctuations when it is held against light, and (ii) paperÕs stress–strain curve displays a smooth transition from a linear to physically nonlinear elastic response in loading (the latter displaying plastic strain upon unloading). The random formation involves small floes of cellulose fibers randomly deposited during the paper manufacturing process, and, given the finite element mesh being five times smaller than the floe diameter, the smallest mesoscale d was set equal to 1/5, rather than 1 as in (2.8–9). Evidently, the larger is the test domain, the smaller is the strain or complementary energy, stored under either (2.2) or (2.3), respectively, and this gives two energy bounds on the effective medium on the infinite scale. Evident from this is the trend to approach the RVE with mesoscale d increasing, and this occurs here approximately on scales about ten times larger than the floe size. The actual choice of

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

1127

Fig. 2. Comparison of shear stress fields according to applied boundary conditions (2.3) and (2.4) and the mesoscale d = 4 or 10; after Ostoja-Starzewski and Castro (2003).

Fig. 3. Stress–strain curves of boundary condition- and scale-dependent shear responses. Upper and lower bounds resulting from uniform displacement (resp. traction) boundary conditions are shown; after OstojaStarzewski and Castro (2003).

1128

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

RVE size depends on oneÕs preferred choice of the acceptable error between both responses. To be exact, paper is an elastic–viscoplastic material. However, techniques presented here can be extended to handle its scale effects as well. This would be a more complicated development, which we leave as one of future ÔchallengesÕ according to the title of this article. 2.2. Power-law hardening materials In the case we have a random, plastic material, the result of the previous subsection can be applied to the loading, but not unloading, response regime. Let us consider a multiphase (p = 1, . . ., ptot) elastic–plastic-hardening material described by an associated flow rule (e.g., Hill, 1950) de0ij ¼

dr0ij of þk dfp orij 2Gp

de0ij ¼

dr0ij 2Gp

de ¼

whenever f p ¼ cp and df P 0,

whenever f p < cp ,

1  2mp dr 2Gp ð1 þ mp Þ

ð2:12Þ

everywhere ðde ¼ deii =3 dr ¼ drii =3Þ

where Gp (shear modulus), mp (PoissonÕs ratio), and cp (yield limit) form a vector, whose each component (described by its indicator, or characteristic, function vp recall (2.7)) gives rise to a scalar random field, such as Gðx,xÞ ¼

ptot X

Gp vp ðx,xÞ 8x 2 X:

ð2:13Þ

p¼1

The entire body B ¼ fBðxÞ; x 2 Xg is described by a random vector field H = {G,m,c}. On the microscale we have tangent moduli – CTd or STd – of the body Bd(x), which connect stress increments with strain increments applied to it dr ¼ CTd de de ¼ STd dr:

ð2:14Þ

The Hill condition, and its implication for the type of admissible boundary conditions, is Z dr de ¼ dr de () ðdt  drnÞðdu  dexÞdS ¼ 0, ð2:15Þ oBd

where dr ¼ dr0 and de ¼ de0 . Next, we recall (e.g., Hill, 1950) two extremum principles: one for kinematically admissible fields Z Z Z Z 1 1 dt d~ u dS  d~ r d~e dV 6 dt du dS  dr de dV ð2:16Þ 2 Bd 2 Bd oBtd oBtd

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

and another for statically admissible fields Z Z Z Z 1 1 dr de dV  dt du dS 6 d~ r d~e dV  d~t du dS: 2 Bd 2 Bd oBud oBud

1129

ð2:17Þ

From this, we find Td Tt Tt 8d0 ¼ d=2: CTd d 6 Cd0 ,Sd 6 Sd0 ,

ð2:18Þ

Upon ensemble averaging, and applying this to ever larger windows ad infinitum we get a hierarchy of bounds on the macroscopically effective tangent modulus CT1 ¼ ðST1 Þ1 1

hSTt 1 i

6    6 hSTt d0 i

1

6 hSTt d i

1

Td 6    6 CT1 6    6 hCTd d i 6 hCd0 i

6    6 hCTd 8d0 ¼ d=2, 1 i, 1 hSTt 1 i

ð2:19Þ

hCTd 1 i

where and are recognized as the Sachs (1928) and Taylor (1938) bounds, respectively. See Suquet (1997) and Ponte Castan˜eda and Suquet (1998) for comprehensive reviews of effective RVE level properties of nonlinear composites. Indeed, one particular problem we studied involved a matrix-inclusion composite material (Jiang et al., 2001), with the matrix phase being elastic-hardening plastic, and the inclusions being elastic and, indeed, of the same modulus as that of the matrix in the elastic range. Clearly, we see the bounding character of responses computed under (2.2) and (2.3), respectively. Patterns of plastic shear bands under these two boundary conditions are displayed in Fig. 4, while the corresponding hierarchy of bounds, computed here on two scales, is shown in Fig. 5. Note: (i) at a smaller scale (d = 6), the response under uniform kinematic b.c. (2.2) is much more uniform than that under uniform stress b.c. (2.3); (ii) at a larger scale (d = 20), the discrepancy due to different types of loading is much smaller, which fact shows a tendency to homogenize with mesoscale tending to 1; (iii) depending on the desired approximation of RVE, we have either or on the lefthand side of (1.1); (iv) the macroscopic response is quite well approximated by the uniform stress (but not the uniform strain) assumption. This situation would be reversed for an elastic–plastic matrix with soft, rather than hard, inclusions, and helps explain why Taylor bound works well for polycrystals with dislocations. Another setting where mesoscale tangent moduli have recently been employed is mechanics of paper, see Section 4.2 below. 2.3. Homogenization in dynamic response In the case of dynamics – that is, when inertial effects cannot be neglected – the Hill condition (2.1) needs to be modified. Starting from the local equations of motion rij,j ¼ q€ ui :

ð2:20Þ

1130

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

Fig. 4. Two samples (realizations Bd(x)) of a random matrix-inclusion composite at d = 6 (left column, top row) and d = 20 (right column, top row); corresponding von Mises strain patterns shown under displacement (middle row) and traction boundary conditions (bottom row).

Wang (1997) found eij rij  eij rij þ ui rij,j  eij rki,k xj ¼

1 jBd j

Z

½ui ðxÞ  ui,j xj   ½ti ðxÞ  r0ki nk ðxÞdS:

oBd

ð2:21Þ

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

1131

Fig. 5. Stress–strain curves of a random matrix-inclusion composite, shown in Fig. 4, with a linear elastic inclusion (i) phase and a power-law hardening matrix (m) phase: ensemble averages of eight sample windows at mesoscale d = 6, and four sample windows at mesoscale d = 20, under displacement, mixed, and traction boundary conditions.

When dealing with strain rates and velocities instead of strains and displacements, the left hand side in the above may also be written as e_ ij rij  e_ ij rij þ 1 q d ðvi vi Þ  e_ ij q_vi xj : 2 dt In the case ðq€ui ¼ 0Þ and when either one of three boundary conditions (2.2)–(2.4) is prescribed, the Hill condition ðeij rij ¼ eij rij Þ is satisfied. However, in the dynamic case, when velocities prescribed on oBd are affine, the right hand side of (2.20) equals zero, and that equation reduces to 1 d e_ ij rij ¼ e_ ij rij þ vi rij,j  e_ ij rki,k xj ¼ e_ ij rij þ q ðvi vi Þ  e_ ij q_vi xj : 2 dt

ð2:22Þ

On the other hand, in the dynamic case with tractions prescribed by (2.3), one finds a simpler relation involving only the rate of strain energy computed from volume averages ð_eij rij Þ, the volume average of rate of strain energy, and the rate of kinetic energy 1 d e_ ij rij ¼ e_ ij rij þ vi rij,j ¼ e_ ij rij þ q ðvi vi Þ: 2 dt

ð2:23Þ

Clearly, we see a generalization of the Hill condition given in (2.1). Wang (1997) has developed dynamic plasticity models of porous materials on the RVE level with this formulation. An outstanding challenge is the development of such models to randomly heterogeneous materials on the SVE (mesoscale) level, with explicit inclusion of scale effects.

1132

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

It is of interest to note the effective equation of motion (for the RVE), derived by Wang and Sun (2003) rij,j þ F i ¼ qu€i ,

ð2:24Þ

where Fi may be regarded as an effective body force resulting from the microinertia.

3. Hierarchies of mesoscale bounds for rigid-perfectly-plastic materials 3.1. Background Here a random rigid-plastic material B ¼ fBðxÞ; x 2 Xg,8x 2 X, is defined by stating that, for any grain (i.e. on microscale), d ij ¼ k

oF ðrij ,xÞ oF ðrij ,xÞ , rij d ij ¼ rij k P 0: orij orij

ð3:1Þ

Thus, in each grain the admissible stress states lie within the set P 1 ðxÞ ¼ frjF ðr,xÞ 6 0g  frjreq 6 rY g, where req is the equivalent stress; rY is the yield stress. Next, on any mesoscale d = L/d, we define volume averages: Z Z 1 1 rðxÞdV dd ¼ dðxÞdV : rd ¼ Bd Bd

ð3:2Þ

ð3:3Þ

For macroscale, Bd!1(x) in the notation of homogenization theory (Suquet, 1986, 1993): R ¼ rd!1 ðxÞ D ¼
dd!1 ðxÞ: ð3:4Þ For the RVE (d ! 1), the yield (or external) surface of the composite B1 delimits the set   P hom ¼ R 2 R33 j9rðxÞ with r ¼ R,div r ¼ 0,F ðr,xÞ 6 0, 8x 2 B ð3:5Þ and the elementary Taylor and Sachs bounds on Phom are expressed via a hierarchy of inclusion relations 
Sachs  Taylor  P hom  RjReq 6 rY : ð3:6Þ RjReq 6 inf rY ðxÞ x2Bd

On any finite mesoscale – a scale below RVE – there holds some form of an apparent yield surface F d ðr,xÞ, bounding some set P d ðxÞ ¼ frd j9rðxÞ,div r ¼ 0,F ðr,xÞ 6 0,8x 2 Bd g: hom

ð3:7Þ

In contradistinction to p , Pd(x) depends on the configuration x and the mesoscale d, and these are two issues we investigate in the following.

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

1133

3.2. Bounding on mesoscales via kinematic and traction boundary conditions Let us start by recalling the upper bound theorem allowing for discontinuities in the velocity field, according to Kachanov (1971), but in a form more suited for our purposes. First, we consider an arbitrary kinematically admissible velocity field v* of the body Bd. If r ij is the stress field associated with v i by (3.1)1 and also satisfying (3.2)2, the basic energy balance equation can now be written as Z Z Z ti v i dS ¼ rij d ij dV þ sY j½v jdS, ð3:8Þ oBd ½v 

S ½v

Bd ½v 



where S ¼ [M m¼1 S m is the set of internal surfaces of discontinuity in v*. Recalling the inequality (e.g., Hill, 1950) rij d ij 6 r ij d ij , we find Z

ti v i dS 6

oBd

ð3:9Þ Z

r ij d ij dV þ sY Bd

Z S ½v 

j½v jdS:

ð3:10Þ

Upon subjecting the body Bd(x) on its entire boundary oBd ¼ oBvd to a uniform kinematic boundary condition (d ij is a constant and equals the volume average of dij). vi ðxÞ ¼ d ij xj

8x 2 oBd ,

ð3:11Þ

noting (3.8), we write this upper bound theorem as Z Z Z Z Z rij d ij dV þ jsY j ½vdS ¼ ti vi dS 6 rij d ij dV þ jsY j Bd

S ½v

oBd

j½v jdS:

S ½ v 

Bd

ð3:12Þ This states that the upper bound on the actual external forces at which plastic deformations begin may be found by assuming an arbitrary kinematically admissible deformation mechanism of the body under consideration. Next, as v i let us take a ÔrestrictedÕ uniform kinematic boundary condition v i ðxÞ ¼ d ij xj

8x 2 [Ss¼1 oBd0s s ¼ 1, . . . ,S,

ð3:13Þ

which acts on all the boundaries of the partition. In the context of Fig. 6(a), (3.13) applies to the external square-shaped boundary oBd as well as to the internal cross. Because the solution under the condition (3.13) on the partition [Ss¼1 Bd0s is a kinematically admissible distribution under the condition (3.11), but, not vice versa, from (3.12) we have Z Z Z Z rij d ij dV þ jsY j ½vdS 6 r ij d ij dV þ jsY j j½v jdS: ð3:14Þ Bd

S ½v

[Ss¼1 Bd0

s

S ½v



This says that the actual external forces at which plastic deformations begin in Bd under (3.11) are bounded from above by forces at which plastic deformations begin in the partition [Ss¼1 Bd0s under (3.13). If by P vd and P vd0s , we denote domains in stress space bounded by volume average yield stresses rdY and rd0 Y corresponding to Bd and [Ss¼1 Bd0s respectively, we can write

1134

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

Fig. 6. Two realizations of a random, two-phase chessboard material, of nominal volume fractions 50%, in mesoscale domains 10 · 10 and 50 · 50. The corresponding apparent rigid-plastic stress–strain responses, at contrast 5, under uniform kinematic and uniform traction boundary conditions of the 10 · 10 and 50 · 50 domains. A partition is shown in the top left figure.

P vd ðxÞ 6 P vd0s ðxÞ:

ð3:15Þ

These two domains are bounded by the surfaces F vd ðr,xÞ and F vd0s ðr,xÞ, respectively. The randomness above – i.e., the explicit dependence on a particular configuration x – can now be removed thanks to the assumption of statistical homogeneity

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

1135

and ergodicity of the material, providing we take ensemble averages of rdY and rd0 Y resulting in average yield surfaces hF vd i and hF vd0 i, respectively. Thus, if c P vd and Pcvd0 v v denote two domains bounded, respectively, by hF d i and hF d0 i, we can write c P vd P vd c

8d0 ¼ d=2:

ð3:16Þ

0

d0s ,

because the information about a particular Here we simply write d in place of part s of the partition is lost. Finally, using induction, we obtain a hierarchy of inclusions P v1    c P vd Pcvd0    c P v1 P v1

8d0 ¼ d=2,

ð3:17Þ

hom

where is the yield locus P of RVE, for which the statistical average may be dropped. On the other hand, c P v1 is the domain bounded by the average yield locus of a single grain, or the Taylor (1938) bound of (3.6), which results from prescribing a uniform deformation rate field everywhere. At present, we do not have a proof based on energy principles, but we do have computational results which indicate that bounding on mesoscales from below occurs via traction boundary condition rij nj ðxÞ ¼ ti ðxÞ ¼ md t0i ðxÞ 8x 2 oBtd :

ð3:18Þ

Here we define P td ðxÞ  P td0s ðxÞ

ð3:19Þ

and our results of Section 3.3 suggest this hierarchy of inclusions P t1      c P td  c P td      c P t1

8d0 ¼ d=2,

ð3:20Þ

where c P td and Pctd0 , denotes a domain in stress space bounded by an ensemble averaged yield surfaces hF td i and hF td0 i. 3.3. Plane deformation rate and plane stress responses of a random chessboard Huber– Mises material 3.3.1. Model formulation We consider a two phase (r = 1,2) material with the microgeometry of a random chessboard in two dimensions (x1 and x2), see Figs. 6 and 7. This material is a system of phase 1 (light red) and phase 2 (dark red) square-shaped grains d · d. More precisely, the microgeometry of phase 1 is taken as a Bernoulli lattice process Pp,d (p is a site probability) on a Cartesian lattice of spacing d in R2 (e.g. Stoyan et al., 1987) Ld ¼ fx ¼ ðm1 d,m2 dÞg

m1 ,m2 2 N:

ð3:21Þ

Any configuration B(x) of this piecewise-constant material is described by an indicator function v(x,x) taking values 1 and 0 according as x falls in phase r = 1 and outside (i.e., phase 2), respectively. We take grains to be incompressible, rigid-plastic, isotropic of the Huber–Mises type (HM), that is, having a yield condition

1136

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

Fig. 7. Apparent plastic responses (shown in terms of the von Mises strain) of the two-phase 10 · 10 chessboard of Fig. 5, in plane deformation rate under kinematic (left column) and traction (right column) boundary conditions.

F ðrij ,xÞ ¼ r0ð2Þ  k 2r ðxÞ ¼ 0,

r ¼ 1,2:

ð3:22Þ

Here r0ð2Þ is the second invariant of the stress deviator r0ij . The above shows that randomness appears at the level of yield stress kr(x) of a homogeneous grain; prime denotes a deviator. As is well known, in this case, the flow rule becomes: d ij ¼ kr0ij whenever F(rij,x) = 0 and dF P 0. Given k1 and k2 of both phases, we also introduce a mismatch (or contrast) a = k2/k1. We set a and the nominal volume fraction of light red and dark red phases to 50%, so that we have P12,d . On finite scales (d < 1) there are, of course, local fluctuations from one Bd(x) to another, and this is reflected in Figs. 6–8, 10 and 11. More specifically, we consider two domain sizes d = 10 and 50. The response of any given

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

1137

Fig. 8. Apparent plastic responses (shown in terms of the von Mises strain) of the two-phase 50 · 50 chessboard of Fig. 5, in plane deformation rate under kinematic (left column) and traction (right column) boundary conditions.

specimen Bd(x), generated according to the Bernoulli lattice process, and subjected to in-plane loadings is obtained by a computational mechanics method described below. We study mechanical responses of 50 specimens at mesoscale d = 10, and 20 specimens at d = 50. These numbers reflect the fact that the larger is the scale, the smaller is the scatter, and hence, the number of realizations required. The computational mechanics method we use is an explicit-in-time finite difference/finite element-in-space program, called FLAC (1999) – the acronym stands for a fast Lagrangian analysis of continua. This method has been extensively used in geomechanics (see Cundall, 1976, 1990), and its major advantage over the conventional implicit finite element method is that the inversion of a large global stiffness matrix is avoided. Another benefit of FLAC is that it can model either quasi-static

1138

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

(as in our study) or dynamic phenomena. Also, FLAC, just like other finite element methods, allows spatially nonuniform assignment of constitutive properties. The general calculation sequence embodied in FLAC first invokes the equations of motion to derive new velocities and displacements from stresses and forces. Then, strain rates are derived from velocities for all of the grid points, and new stresses from those strain rates. One time step is taken for every cycle around the loop. The central concept in FLAG is that the computational wave speed always keeps ahead of the physical wave speed, so that the equations always operate on known values that are fixed for the duration of the calculation. There are several distinct advantages to this approach, but, most importantly, no iteration process is necessary when computing stresses from strains in an element, even if the constitutive law is very nonlinear. 3.3.2. Plane deformation rate It is well known that the yield condition (3.22) now simplifies to F r ðrÞ ¼ ðr11  r22 Þ2 þ 4r212  4k 2r , r ¼ 1,2

ð3:23Þ

and the flow rule (2.12) leads to d 11 d 22 d 12 ¼ ¼ , r11  rm r22  rm r12

ð3:24Þ

where rm = (r11 + r22)/2 is the mean stress. In Figs. 7 and 8 (top pictures) we show apparent plastic responses – in terms of the von Mises strain – of the HM two-phase chessboards at contrast 5, of respective Fig. 6(a) and (b) on mesoscales d = 10 (a) and 50 (b). The loading is of shear type d 011 ¼ d 022 , d 012 ¼ 0

ð3:25Þ

for the kinematic boundary conditions (left column), and r011 ¼ r022 , r012 ¼ 0

ð3:26Þ

for the traction boundary conditions (right column). In the latter case, the loading was increased until the deformation started. In each case, bottom pictures show the resulting plastic zones superposed onto the original mesh. The inhomogeneity of strain fields, especially under the traction conditions, is evident. As the scale d increases, there is a clear tendency for responses under both types of loading to attain spatial uniformity and similarity. While this tendency to homogenization was noted in a similar study in (Jiang et al., 2001), the key difference of the present model from the matrix-inclusion composite studied there is that the random two-phase chessboard at nominal 50% volume fraction is below the percolation of either phase, because the critical probability in a site-percolation problem is about 59% (Ziman, 1979). [Of course, the actual volume fraction varies from one x to another.] Thus, in general, we do not see shear bands spanning the entire material domain through zones of soft material only, because the soft phase does not percolate (or interconnect the entire material domain). To enable a global plastic flow, portions of hard material need to yield in the most critical spots

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

1139

separating the soft phases. Strictly speaking, the shear bands probably involve infinitesimally thin slip-lines, which, however, cannot be well resolved given our mesh with one finite element per one chessboard square. The trend of theoretically derived bounding yield surfaces to converge with increasing mesoscale is displayed in Fig. 9. In particular, for two mesoscales d = 10 and 50 we show apparent yield surfaces under kinematic and traction loadings. Interestingly, the ensemble averaged yield loci under the latter loading do show the scale effect: larger scale leads to a higher yield limit. The corresponding scatter, only shown for kinematic loading at d = 10, decreases with d going up. The vectors of deformation rate d (not drawn here) do satisfy the normality: flow rule holds in the ensemble average sense, but generally not per single realization x 2 X. Before we leave this section, we notice the similarity of our patterns to those in a study of stress and inelastic strain fields in heterogeneous bond coat microstructure thermal barrier coatings (TBCs) due to Pindera et al. (2004); note the random-microstructure strain localization effects evident in Figs. 7 and 8 there.

Fig. 9. An ensemble view of apparent yield surfaces of random Huber–Mises two-phase chessboards in plane deformation rate on mesoscales d = 10 and 50. For reference, the hard and soft phase yield loci are shown too.

1140

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

3.3.3. Plane stress The yield condition (3.22) now simplifies to F r ðrÞ ¼ r211  r11 r22 þ r222 þ 3r212  3k 2r ,

ð3:27Þ

and the flow rule becomes d 11 d 22 d 12 ¼ ¼ : 2r11  r22 2r22  r11 6r12

ð3:28Þ

It will turn out below that an anisotropic form of apparent yield surfaces Fd(r) will be relevant; superscript v or t will apply accordingly. A convenient form of Fd(r) in the space of principal stresses r1 and r2, actually employed for homogeneous transversely isotropic materials (e.g. Mellor, 1982), could be

Fig. 10. Apparent plastic responses (shown in terms of the von Mises strain) of the two-phase 10 · 10 chessboard of Fig. 5, in plane stress under kinematic (left column) and traction (right column) boundary conditions.

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160 2

2

F d ðrÞ ¼ ðr1 þ r2 Þ þ ð1 þ 2Rd Þðr1  r2 Þ  2ð1 þ Rd Þr2dY :

1141

ð3:29Þ

Here rdY is the yield stress in uniaxial tension, with d denoting its scale dependence; in general, Rd is scale-dependent too. Figs. 10 and 11 display apparent plastic responses – in terms of the von Mises strain – of the Huber–Mises two-phase chessboards at contrast 5, of respective Fig. 6(a) and (b) on mesoscales d = 10 (a) and 50 (b). The loading is again of shear type: in the left column for the kinematic boundary conditions (3.4), and in right column for the traction boundary conditions (3.5). As before, in each case, bottom pictures show the resulting plastic zones superposed onto the original mesh.

Fig. 11. Apparent plastic responses (shown in terms of the von Mises strain) of the two-phase 50 · 50 chessboard of Fig. 5, in plane stress under kinematic (left column) and traction (right column) boundary conditions.

1142

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

The inhomogeneity of strain fields, especially under the kinematic condition, is evident. However, as the scale d increases from 10 to 50, there is a clear tendency for responses under both types of loading to attain spatial uniformity and similarity. It is interesting to compare the apparent plastic responses of Figs. 10 and 7 for the two-phase 10 · 10 chessboard: there is a strong difference between the kinematiccontrolled responses (left columns), but not between the traction-controlled (right columns), as we switch from the plane deformation rate to the plane stress. This comparison may be continued for the two-phase 50 · 50 chessboard by considering apparent plastic responses (always shown in terms of the von Mises strain) in Fig. 11 versus Fig. 8. In Fig. 11 note the formation of shear band features under the kinematic condition, which was brought out at 0.01 of the maximum von Mises strain.

Fig. 12. An ensemble view of apparent yield surfaces under kinematic conditions in plane stress on mesoscales d = 10 and 50. Hard and soft phase yield loci, as well as an average kinematic-controlled locus and two most extreme ones on d = 50, are shown.

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

1143

Fig. 12 gives an ensemble view of apparent yield surfaces in the (r1,r2) plane on mesoscales d = 10 and 50, whereby, respectively, fifty and twenty realizations were employed. Hard and soft phase yield loci at contrast 5, having the classical HM

Fig. 13. A close-up from Fig. 11 of the ensemble view of apparent yield surfaces in plane stress, focusing on the first quarter in the space of principal stresses, on mesoscales d = 10 and 50. The
d vectors, shown for ensemble average yield surfaces at d = 10 and 50 under kinematic and traction conditions, display the departure from normality.

1144

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

elliptical shape, are shown for reference. Also drawn, is the ensemble average kinematic-controlled yield locus and two most extreme ones for scale d = 50; their symmetry with respect to the origin of (r1,r2)-plane was noted. Two Huber–Mises ellipses are drawn in broken lines for reference so as to display the departure of these apparent (mesoscale) responses from the HM response. This departure is made more visible in Fig. 13 representing a close-up view of a part of Fig. 12, focusing on the quarter of apparent yield surfaces ranging from pure shear up to uniaxial (ex)tension, on mesoscales d = 10 and 50. It is evident that the form (3.29) provides a convenient model, whereby the Rd value is about 0.7 albeit without betraying any clear trend with d. Larger mesoscale domains will need to be studied in the future to uncover such a trend. It follows that, to set up a simplest stochastic model, one can take Rd = 1 and deal only with scale dependence in rYd in (3.29), plus its scatter about this mean value. That model could describe the ensemble of yield surfaces Fd ¼ fF d ðxÞ; x 2 Xg,

ð3:30Þ

where Fd is of an elliptical (finite thickness) shell shape. The shell converges to an infinitesimal thickness in the d ! 1 limit, but bifurcates into two distinct HM ellipses in the d ! 1 limit. Fig. 13 also displays the hdi vectors, ensemble averaged for each one of four yield surfaces here (at d = 10 and 50 under both kinematic and traction conditions). Noteworthy is the departure of apparent (mesoscale) flow rule, for each case, from the normality. 3.3.4. Effects of mesh refinement To study the effect of increasing mesh resolution, some of the realizations on mesoscale d = 10 were solved for their apparent yield surface and flow rule using a grid of 50 · 50 instead of 10 · 10 elements. In practice, however, a 50 · 50 mesh thus introduced can be thought of as a particular realization on a mesoscale d = 50, where the hard and soft phases appear to be lumped in 5 · 5 squared clusters. Therefore, the issue of mesh dependence of a computational mechanics solution cannot be easily divorced from the issue of scale dependence. Our results on the effects of such a five-fold mesh refinement -and beyond – are summarized as follows. 3.3.4.1. Plane deformation rate. As the mesh resolution increases five-fold under the traction boundary condition, the apparent yield surface of each realization, as well as its ensemble average, moves in the space of principal stresses towards the soft phaseÕs yield surface. As it is refined further, there is almost no scatter of the apparent yield limit among the realizations, and normality is recovered. In the case of the five-fold mesh refinement under the kinematic boundary condition, the apparent yield surface of each realization Bd = 10(x) moves towards the ensemble of yield surfaces of Bd = 50(x) materials. The ensemble average yield surfaces coincide then. As in the case of a coarse mesh, the vectors of deformation rate satisfying normality of flow rule hold in the ensemble average sense, but not per single realization. Finally, the scatter of yield points around the ensemble average is the same as with the coarse mesh.

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

1145

3.3.4.2. Plane stress. As the mesh resolution increases five-fold under the traction boundary condition, the apparent yield surface of each realization, as well as its ensemble average, moves in the space of principal stresses towards the soft phaseÕs yield surface. As it is refined further, the apparent yield occurs at the yield limit of the soft phase and normality is recovered. On the other hand, as the mesh resolution increases five-fold under the kinematic boundary condition, the apparent yield surface of each realization Bd = 10(x) moves towards the ensemble of yield surfaces of Bd = 50(x) materials. The ensemble average yield surfaces then coincide. 3.4. Observations  Neither in plane deformation rate, nor in plane stress, is the flow rule normal for any given realization of the random medium.  In plane deformation rate, the flow rule is normal upon ensemble averaging, but not so in plane stress.  In plane deformation rate the yield surface, for any d and x of random medium, consists of two parallel lines. In plane stress, the yield surface is ellipse-shaped, but, for both bcÕs (kinematic and traction), of a smaller aspect ratio than that of the HM material. This suggests a convenient form of Fd(r) for homogeneous transversely isotropic materials (e.g. Mellor, 1982) 2

2

F d ðrÞ ¼ ðr1 þ r2 Þ þ ð1 þ 2Rd Þðr1  r2 Þ  2ð1 þ 2Rd Þr2dY ,

ð3:31Þ

where rdY is the yield stress in uniaxial tension and Rd . 0.7.  It is not yet clear what scale is necessary to approximately attain the classical HM material for the RVE. All we can say is that much larger domains need to be simulated, and these involve much more computing power.  Scatter of yield loci is decreasing with mesoscale increasing, under either type of boundary conditions, be it plane strain or plane stress.

4. Patterns of plastic responses in random materials 4.1. Geodesic properties of shear-band patterns It is clear from Fig. 4 that shear bands display geometric shapes which conform to the actual spatial distribution of the material microstructure. With the latter being spatially irregular, the shear bands are irregular too, but have a certain overriding characteristic: they seem to minimize the path through the specimen at angles roughly 45 from the principal axes of tensile loading. Conceptually, therefore, shear bands may well behave as geodesics – curves of shortest path joining two specific points in space, say, P1 and P2. Let us recall here that the geodesic distance between P1 and P2 is defined as

1146

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

dðP 1 ,P 2 Þ ¼ inffall paths from P 1 to P 2 in a set in R2 g:

ð4:1Þ

In accordance with this hypothesis, the geodesics should join the opposite faces of a specimen Bd(x), such as that in Fig. 14(a), to which pure shear loading is applied. Taking the exactly same material model as in Section 2.2, Fig. 14(b) displays patterns of equivalent plastic strain e, found by computational micromechanics, under uniform kinematic loading (2.2) on the entire boundary oB with e011 ¼ e022 ¼ 0, e012 ¼ e021 . While in a homogeneous material perfectly straight (horizontal and vertical) shear bands would form, in the heterogeneous material of Fig. 14(a) a distortion of these bands, due to the presence of elastic inclusions, occurs. We postulate that these distorted shear bands follow the shortest paths in the matrix, and proceed to estimate, by geodesic propagations, shortest paths on the two-dimensional configuration of the composite of domain Bd(x). As potential shear bands, we therefore consider two families of shortest paths obtained by geodesic propagations for each and every point in the material domain in two orthogonal directions that avoid (black) inclusions. More precisely, for a horizontal geodesic propagation, for each and every point in the material domain, as source S and destination D we first take the left and right faces, and then invert their roles thus obtaining the set G – hor. For a vertical geodesic propagation as source S and destination D we first take the top and bottom faces, and then invert their roles obtaining the set G – ver. We consider propagations obtained on hexagonal lattices, so that we deal with the so-called hexagonal geodesic distances. For basic concepts of this methodology, in the context of fracture mechanics of random media, see (Jeulin, 1993). Final results, following Jeulin and Ostoja-Starzewski (2000), are synthesized as follows:

Fig. 14. (a) A matrix-inclusion composite; (b) shear bands (bright points) obtained by computational micromechanics under pure shear loading e011 ¼ e022 ¼ 0, e012 ¼ e021 .

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

1147

(1) Addition of G – hor with G – ver to obtain the set G – add, Fig. 15(a). (2) Supremum of G – hor with G – ver to obtain the set G – sup, Fig. 15(b). The grayscale at any given point in Fig. 15(a) as well as 15(b) indicates the (appropriately normalized) distance it takes to connect that point by a geodesic to the opposite edges of the square shaped domain Bd(x): the darker is the point, the shorter is this distance. Given two fields such as {g(i); i = 1, . . ., I} and {e(i); i = 1, . . ., I} both defined on a square lattice of the same size (256 · 256), so that I = 2562, we can compare them just as we would compare two random variables, say, g and e by using the normalized covariance (or cross-correlation coefficient) qge ¼

h½g  hgi½e  heii , rg re

ð4:2Þ

wherein rg and re are standard deviations. Noting that the numerator equals Ægeæ  Ægæ Æeæ, and invoking the ergodic assumption, we therefore compute qge from PI PI PI 1 1 1 i¼1 gðiÞeðiÞ  ½ I i¼1 gðiÞ½ I i¼1 eðiÞ I ffi: qge ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð4:3Þ P P P PI I I 2 I 2 f1I i¼1 g2 ðiÞ  ½1I i¼1 gðiÞ gf1I i¼1 e2 ðiÞ  ½1I i¼1 eðiÞ g It turns out that the cross-correlation of G – sup (Fig. 15(b)) with the true solution, e, obtained by finite elements (Fig. 14(b)) is only 0.2, and a lower number ( 0.1) is obtained for the G – add geodesic propagation (Fig. 15(a)). The advantage of this purely geometric method is an extremely fast computation of the pattern of

Fig. 15. (a) Shortest paths (in dark) that avoid (black) inclusions obtained by addition – or supremum – of horizontal with vertical geodesics. Note that the darkest (brightest) points here correspond to the brightest (respectively, darkest) ones in (b) of the previous figure.

1148

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

plastic deformation as opposed to a full computational mechanics approach. This may then offer a very rapid way of a purely geometric assessment of zones of plastic flow in disordered heterogeneous materials, although much further research needs to be done to assess the appropriate ranges of parameters. One situation where the geodesics may be very useful occurs in wintry ship operations, when a highly heterogeneous ice field is covering a large body of water. Given only a finer scale image (Fig. 1(a)) or, say, a coarser scale image (Fig. 1(b)) – with gray scale roughly indicative of the ice thickness, and hence of its strength – an ice breakerÕs captain wants to find the shortest path, yet one of the least resistance (!), to go from point A to B. 4.2. Random formation versus inelastic response of paper One of the key problems in paper physics concerns the nature of constitutive law of paper and its relation to a spatially inhomogeneous paper structure, the so-called formation. The latter is commonly seen as grayscale fluctuations when a sheet of paper is held against light, and is a reflection of the fact that mass per planar area of the paper (the so-called basis weight, b, or weight per unit area) is non-uniform. While all the paper industries worldwide strive for the most uniformly manufactured paper, some random fluctuations of formation always remain due to the unavoidable presence of flocs, streaks, and, last but not least, the random shapes and sizes of fibers which make it impossible to construct a uniform fiber network. Recently, formation has been found as a possible cause and explanation of the special elastic orthotropy of paper (Ostoja-Starzewski and Stahl, 2000), which says that the in-plane shear modulus of paper is invariant with respect to in-plane rotations. Now, on conventional time scales where viscous effects may be disregarded, and beyond the elastic response regime, paper may be taken as a planar, elastic–plastichardening material. The constitutive coefficients (especially the shear moduli), however, are difficult to measure precisely and are highly dependent on the formation. The measurement of formation is conventionally done by either radiography with resolution down to one millimeter, or, as proposed in Castro and Ostoja-Starzewski (2003), with a scanner with an order of magnitude better quality. An image obtained by a scanner is shown in Fig. 16(a). These considerations offer a basis for carrying out a cross-correlation of two planar fields: the formation and the strain (or stress) field resulting from a computational mechanics simulation of a sheet of paper subjected to some in-plane loading. The simulation requires some assumption of constitutive coefficients in function of the formation, and the higher the cross-correlation comes out to be, the better is the postulated constitutive law. This is illustrated in Fig. 16, which on the left shows the scanned image and the von Mises strain field (e) resulting from a computational mechanics simulation of a uniaxial tension test in the x-direction. The apparent spatial correlation of b with, say, e is quantified as follows: we have two number series: b(i) and e(i), i = 1, . . ., 482 (because the images are defined on a square 48 · 48) each taking values from 0 to 255. Both series being of the same

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

1149

Fig. 16. Cross-correlation between formation field (a) and strain patterns (b–d) for a paperboard (50 mm · 50 mm); after Ostoja-Starzewski and Castro (2003).

length I, the normalized cross-correlation coefficient of b(i) and e(i) is computed by a formula entirely analogous to (4.1), that is PI PI PI 1 1 1 i¼1 bðiÞeðiÞ  ½ I i¼1 bðiÞ½ I i¼1 eðiÞ I ffi: qbe ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð4:4Þ PI 2 PI PI 2 PI 2 2 1 1 1 fI i¼1 b ðiÞ  ½ I i¼1 bðiÞ gf I i¼1 e ðiÞ  ½1I i¼1 eðiÞ g The cross-correlation is 0.77; and very similar numbers are obtained for crosscorrelation of formation with the e11-strain (Fig. 16(c)), as well as for cross-correlation of formation with the volumetric strain increment (Fig. 16(d)). Among the paper physicists there is a simple idea that the higher is the basisweight, the stiffer and stronger is the paper locally. But over the years this has been countered by this consideration of the manufacturing process: paper undergoes pressing by series of large steel rolls, so that areas of high concentration of cellulose fibers are squashed more than those of low to medium concentration. It is therefore speculated by some researchers that the relation between basis-weight and paper properties is not proportional. Our results are in support of the proportional relation for paper with an allowance for some added nonlinearity.

1150

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

4.3. When shear-bands become slip-lines. . . 4.3.1. Stochastic slip-line systems Materials studied in the previous two subsections were of hardening type, and plasticity occurred in them through shear-bands. When hardening vanishes and the material becomes perfectly-plastic, shear-bands become infinitesimally thin and turn into slip-lines. One such system studied by us involved a planar rigidperfectly-plastic material described by a random field of plastic limit in shear, k, with continuum realizations (Ostoja-Starzewski, 1992; Ostoja-Starzewski and Ilies, 1996). Thus, in the present review, we are looking at a mesoscale model, such as depicted in Fig. 1(b), governed pointwise by an isotropic form of Huber–Mises or Tresca type in plane strain: F d ðrÞ ¼ ðr11  r22 Þ2 þ 4r212  4k 2d :

ð4:5Þ

Here kd is the random field fk d ðx,xÞ; x 2 X,x 2 R2 g, which depends on the mesoscale d. A particular choice of the latter defines the random medium B ¼ fBðxÞ; x 2 Xg. Just like in many other models of random media with continuous realizations, the plastic limit is taken as the mean Ækdæ plus a zero-mean fluctuation k 0d k d ðx,xÞ ¼ k 0 þ k 0d ðx,xÞ, hk d i ¼ k 0 :

ð4:6Þ

Considering the equilibrium equations and the yield condition (4.5) of any realization B(x), we have a quasi-linear hyperbolic problem, which, just like conventional homogeneous media problems, lends itself to a solution by the method of characteristics (slip-lines). [Note: in the ensemble sense we have a stochastic hyperbolic problem.] Thus, upon introduction of two functions p and u according to r11 ¼ p þ k cos 2u, r22 ¼ p  k cos 2u, r12 ¼ k sin 2u,

ð4:7Þ

a subsequent substitution of the above into the equilibrium equations, and setting u = p/4 on differentiation, we get op ou ok þ 2k ¼ , ox ox oy

op ou ok  2k ¼ : oy oy ox

ð4:8Þ

Here the rectangular axes are now along the local slip-line directions. Eq. (4.8) can be made independent of the orientation of axes by replacing o/ox1 and o/ox2 by the tangential derivatives o/osa and o/osb along the sa and sb characteristics, respectively, to get dp þ 2k du ¼

ok dsa , osb

dp  2k du ¼

ok dsb , osa

ð4:9Þ

with k 0 = 0, the right-hand sides become zero: we recover the relations of classical plasticity. In any case, the characteristic directions are given by tangents dx2 ¼ tanðu þ p=4Þ, dx1

dx2 ¼ tanðu  p=4Þ: dx1

ð4:10Þ

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

1151

Eqs. (4.9) and (4.10) form the basis for determination of the Hencky–Prandtl network of slip-lines in a given boundary value problem. This relies on the method of finite differences for finding x1, x2, p, and u at a new point N, given the data at the two preceding points i = 1,2, Fig. 17. Now, due to the randomness in k1, k2, and kN, as well as the possible randomness in the initial data (p,u)1 and (p,u)2 at the points 1 and 2, two characteristics of the deterministic problem are replaced by two forward causality wedges, Fig. 18. This is seen in parts (b) and (c), where the plastic limit is taken according to (4.6), while part (a) depicts the solution of a reference Cauchy problem for a homogeneous medium where k d ðx,xÞ ¼ k 0 ,

k 0d ðx,xÞ ¼ 0

8x,x:

ð4:11Þ

The Cauchy problem is defined as one in which the normal stress and shear stress are specified on a segment [0.0,9.0] of the x1-axis; this segment intersects each of the characteristics only once, and on it we assume r22 ¼ 1:5, r12 ¼ 0:6 þ x1 =20:

ð4:12Þ

Adopting k0 = 1, we consider two cases of fluctuation: (i) ‘‘very small noise’’ with k 0d taken as a uniform random number in the interval k 0d 2 ½0:0025,0:0025, (ii) ‘‘small noise’’ with

k 0d

k 0d 2 ½0:025,0:025:

ð4:13Þ

taken as a uniform random number in the interval ð4:14Þ

Fig. 17. The setup of two mutually orthogonal characteristics a and b (modeling two slip-lines) superposed on a disordered, granular microstructure, showing that the choice of spacing of points 1 and 2 introduces a mesoscale in the problem. Three mesoscale windows are shown, centered at points 1, 2, and 3, respectively; after Ostoja-Starzewski and Ilies (1996).

1152

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

Fig. 18. Solution of a Cauchy problem with inhomogeneous data on y = 0, in function of increasing random noise in the material: (a) no noise, (b) very small noise, and (c) small noise; after Ostoja-Starzewski and Ilies (1996).

With this, in Fig. 18 we study the effects of increasing randomness in the plastic limit on slip-line patterns. Overall, we find that in the Cauchy as well as the characteristic problems, those effects are strongly amplified by any in-homogeneity in the boundary data. Furthermore, the ensemble average of the solution is very close to the solution of the homogeneous medium problem for very small noise, but the type of randomness of kd (uniform versus Gaussian) has but a weak effect. 4.3.2. Limit analysis of a pipe under internal loading As a more practical application of this method we consider limit analysis of a cylindrical tube loaded by a uniform traction on its internal boundary. Approach to this classical problem (e.g., Kachanov, 1971) in the case of a homogeneous mate-

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

1153

rial is based on the solution of a Cauchy problem for an axisymmetric stress field (rr, rh, rrh) in the (r,h) coordinate system in a plate with a circular hole, under a pressure boundary condition rr ¼ p < 0

rrh ¼ 0

at r ¼ a:

ð4:15Þ

With the yield condition rh  rr = 2k the stress field is determined explicitly by the formula r rh ¼ rr þ 2k, rr ¼ p þ 2k ln ð4:16Þ a and the slip-line field is now given by a system of two orthogonal, logarithmic spirals r r ð4:17Þ h þ ln ¼ a h  ln ¼ b: a a It follows from (4.12)1 that, in case of a tube of external radius b, the condition rr = 0 at r = b gives the limit pressure b p ¼ 2k ln : a

ð4:18Þ

Having this formula, we may regard p* as given and ask for the external radius b as a function of p*. To make the example concrete, we set k = 1.5, p* = 3.0 and a = 1, which then results in bðp Þ ¼ ep =2k ¼ 2:723. In the case of a tube made of an inhomogeneous material, however, the net of sliplines becomes distorted from such perfect pattern and the above ÔniceÕ relations no longer apply. The finite differencing outlined above has now to be used to determine the slip-lines as well as the stress field in any particular realization of B. Now, the evolution of rr along any slip-line emanating from r = a is a certain random walk in the body domain D starting from the value p* there and stopping at 0 at r = b. Thus, the condition rr = 0 plays the key role in the definition of an excursion set of the random field rr(r,h) = {rr(x); x 2 X} (e.g. Adler, 1981) A0 ðrr ,DÞ ¼ fðr,hÞ 2 Djrr ðr,hÞ P 0g:

ð4:19Þ

This then leads to a set of level crossings oA0 ðrr ,DÞ ¼ fðr,hÞ 2 Djrr ðr,hÞ ¼ 0g,

ð4:20Þ

where oA0(rr,D) is a set of closed contours of plastic zone, which in the case of a homogeneous material Bdet under the pressure boundary condition would be a circle of radius b(p*) = 2.723. Thus, in any given random material the set of level crossings is a random set in plane, which, assuming spatial homogeneity and isotropy of field {kd(x,r,h); x 2 X,(r,h) 2 D}, is a circle containing all the possibilities. Now, let bmax and bmin be a maximum and minimum distance for a given realization, respectively, from the origin to the contour. Since bmax determines the minimal amount of material needed for a tube to withstand the internal pressure p*, we ask: Is bmax(p*) smaller, equal to, or larger than b(p*) of the homogeneous medium problem? In Fig. 19(a) we plot patterns of slip-lines – actually as wedges of forward dependence – under condition (4.15) corresponding to four hundred realizations of B with

1154

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

Ækdæ = 1.5 and k 0d 2 ½0:025,0:025. The set of level crossings is shown as a ring containing all four hundred piecewise-constant, non-circular closed curves. Clearly, bmax(p*) is always larger than b(p*),whereby bmax increases as the noise level in k increases. Thus, the principal conclusion is that the presence of material inhomogeneities requires a thicker tube than what is predicted by the homogeneous medium theory. Let us also note that the higher the pressure p, the greater is the external radius b, and hence, the greater is the spread of the forward evolution cones implying

Fig. 19. Slip-line patterns in a randomly inhomogeneous material, with k 0 following (4.14), under: (a) boundary condition (4.15) and (b) boundary condition (4.21). In each case, 400 realizations of the random medium were used to generate wedges of characteristics; after Ostoja-Starzewski and Ilies (1996).

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

1155

an increase in the scatter of bmax and bmin. Interestingly, both random variables bmax and bmin are symmetrically distributed about the deterministic radius b(p*) of the homogeneous medium problem. The same qualitative conclusions carry over to the case of a combined pressure and shear boundary condition rr ¼ p < 0

rrh ¼ q 6¼ 0

at r ¼ a:

ð4:21Þ

The slip-lines are still spirals, albeit no longer of logarithmic type. The stress field is now given by Kachanov (1971) 2 3 ffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi q q r 2 r 2  r 4 q 2 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2ffi ð þ Þ  þ a q k k a a 6 7  rr ¼ p k 42 ln pffiffiffiffiffiffiffiffiffiffiq pffiffiffiffiffiffiffiffiffiffiq  þ 1 5, r a k k 1kþ 1þk rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a 4ffi 2 2 rh  rr ¼ k  q : r

ð4:22Þ

The general character of these spirals is displayed in Fig. 19(b), and they would be perfectly recovered had the noise k 0d been zero. Thus, we have here a random slip-line network, made of an ensemble of four hundred realizations, corresponding to q* = 1.3, with all the other parameters being the same as before. These plots result in the external radius b(p*,q*) = 3.028. Note: the addition of shear traction has a strongly amplifying effect on the scatter of dependent field quantities, and most notably on the spread of slip-lines. 4.3.3. Related work on stochastic plasticity To the best of our knowledge, the need to study plasticity of randomly inhomogeneous media was pointed out for the first time by Olszak et al. (1962). While randomness per se was not treated, this important reference provided, among others, a very good discussion of methods for solution of boundary value problems of plasticity of inhomogeneous media described by deterministic fields, which, in principle, form a starting point for stochastic problems. Discussed therein were, four principal kinds of solution methods: analytical, approximate (perturbation-type), inverse and semi-inverse, and numerical. Since a random medium is understood as a family of deterministic specimens – i.e. as a random field – the difficulty in finding an explicit analytical solution in the deterministic problem definitely carries over to the stochastic one. We observe here that the analytical solution can be found in a limited number of specific problems described by deterministic fields only. Next, the perturbation method solutions are not placed on a firm theoretical foundation even in the deterministic case. Finally, the inverse-type methods would require the knowledge of an explicit solution in a corresponding stochastic problem of, say, elastic-type, which is a formidable obstacle itself. Plasticity of randomly inhomogeneous media has recently been studied by Nordgren (1992) and Ku and Nordgren (2001) from a different standpoint. The focus there has been on the formulation of upper-bound and lower-bound theorems for random plastic media in combination with methods of reliability theory, and

1156

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

application thereof to plastic collapse of a wedge, a half-plane, and a half-space. The collapse load calculation based on the upper- (respectively, lower-) bound theorem involves a nonlinear (respectively, linear) programming method. Material yield limit has been assumed as kðx,xÞ ¼ k 0 þ k 0 ðx,xÞ,

k 0 ¼ hki,

ð4:23Þ

0

with k being a zero-mean Gaussian random field with an exponentially decaying covariance function ( "

2

2 #) x1  x2 y1  y2 2 Cðx1 ,y 1 ; x2 ,y 2 Þ ¼ r exp p þ , ð4:24Þ rx ry where r2 is the variance, while rx and ry are the correlation lengths in the x and y directions. A different approach was recently advanced by Anders and Hori (1999, 2001), who employed a stochastic finite element approach based on the Karhunen–Loe`ve polynomial chaos expansion (Ghanem and Spanos, 1991) to study 3-D softening elasto-plastic bodies. The material spatial randomness was assumed in the stochastic plastic multiplier, the stiffness tensor, and the yield surface (with an associated flow rule). The presence of random elastic modulus led the authors to compute the mean behavior through a bounding media analysis; the latter involved two materials: one defined by the arithmetic average of E, i.e. ÆEæ, and the other by its harmonic average ÆE1æ1. Stochastic methods have recently been introduced in research on metal forming. These methods combine an analytical/numerical analysis with a computational methodology; see (Sluzalec, 2004). 5. Conclusions This paper was written with an ambitious aim in mind: to outline the status of plasticity theory of random materials, and outstanding challenges. While the results are incomplete, we hope we have met at least some expectations that a reader of such a review paper might have, whereby we point out that the presentation was done from the particular perspective of our own work. Thus, we have said was basically focused on the mechanics of plastic materials where the separation of scales does not hold, and, in particular, on solving the resulting boundary value problems either to bridge the micro-mesoscale gap, or the meso-macroscale gap. The first class of these problems pertains to the setup of SVE, rather than the conventional RVE, while the second one concerns the solution of a stochastic – rather than a deterministic – boundary value problem of plasticity. We make these observations to summarize the current status of plasticity of random media: 1. At the mesoscale level, the yield condition is generalized to bound from above (and below) the macroscopic homogenized response when kinematic (respectively, traction) boundary condition is applied. In essence, the larger is the mesoscale, the

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

1157

tighter are the bounds, albeit the kinematic condition displays a much stronger scale-dependence than the traction one. Just like in the case of elastic and viscoelastic materials studied earlier in the nineties, one needs a combination of extremum principles with material ergodicity and homogeneity. At present, though, the proof of bounding character of traction boundary condition for a rigid-perfectlyplastic material is still lacking. Some more specific issues are: (a) Even though at the microscale level the material possesses normality, neither in plane deformation rate, nor in plane stress, is the mesoscale flow rule normal for any given realization of the random medium. (b) In the case of plane deformation rate, the flow rule is normal upon ensemble averaging. This property does not hold in the plane stress case. (c) In the case of plane deformation rate the yield surface, for any given mesoscale and realization of the random medium, consists of two parallel lines. In the plane stress case, the yield surface is ellipse-shaped, albeit, for both boundary conditions (kinematic and traction), of a smaller aspect ratio than that of the HM material. It is not yet clear what scale effects, beyond those found here, apply. In other words, what mesoscale is necessary to approximately attain the classical HM material for the RVE? An answer to this question, will very likely require much larger domains, which, in turn, will necessitate more computational work. (d) The scatter of yield loci is decreasing with mesoscale increasing, under either type of boundary conditions. It is yet to be seen what behaviors will be displayed by other than HM materials. (e) All the numerical results of Section 3.3 have, strictly speaking, been obtained for one particular case of volume fraction of a specific random material – a planar, two-phase chessboard at 50% – and it still remains an open issue to see to what extent do the particular behaviors found here carry over to other models. Also, other microgeometric features – such as shapes of inclusions may have specific effects on the mesoscale and macroscale yield surfaces. 2. Three proposals on the use of geometric patterns of plastic zones have been outlined here. (a) One proposal involved a very inexpensive (i.e. computationally fast) determination of geodesies in place of the precise but costly solution of, say, an elastic– plastic heterogeneous material by computational mechanics. However, the cross-correlation was rather low, and more work needs to be done to find what material models (other hardening laws and random geometries) would display a better match. (b) It was proposed to use cross-correlations of computationally found strain fields with the actual material microstructure in order to determine the very-hard-to-measure microscale material properties. New work needs to be done on materials other than paper, and also employing experimentally found strain fields. (c) Solution of hyperbolic boundary value problems of rigid-perfectly-plastic materials on meso-to-macroscales involves scattered slip-lines. This technique is rapid, but locally anisotropic, rather isotropic yield functions should be

1158

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

used, and these must, in principle, come from studies of the SVE discussed above. Finally, finite elements offer a more general framework than finite differences for solution of stochastic boundary value problems, and these are appropriately called Ôstochastic finite elements.Õ 3. It is well known that the thermomechanics with internal variables (so- called T.I.V.) offers a general framework for inelastic material behaviors. Basically, the RVE response in T.I.V. is described by the free energy (W) and dissipation (U) functionals, e.g. (Collins and Houlsby, 1997; Germain et al., 1983; Maugin, 1999; Ziegler, 1983). Below the RVE level W and U are random functionals. Scale dependent hierarchies of bounds on them, similar in vein to what we have discussed in this paper, may be developed along the lines suggested in (Ostoja-Starzewski, 2002a,b, 2004). See also Spearot et al. (2004) for a conceptual framework involving internal variables at the SVE level.

Acknowledgement We have benefitted from various constructive comments of M.-J. Pindera and the anonymous reviewers. This material is based upon work supported by the CRC, NSERC, and CFI.

References Adler, R.J., 1981. The Geometry of Random Fields. John Wiley, New York. Anders, M., Hori, M., 1999. Stochastic finite element method for elasto-plastic body. Int. J. Numer. Meth. Eng. 46, 1897–1916. Anders, M., Hori, M., 2001. Three-dimensional stochastic finite element method for elasto-plastic bodies. Int. J. Numer. Meth. Eng. 51, 449–478. Castro, J., Ostoja-Starzewski, M., 2003. Elasto-plasticity of paper. Int. J. Plast. 19 (10), 2083–2098, 2003. Clayton, J.D., McDowell, D.L., 2003. A multiscale multiplicative decomposition for elastoplasticity of polycrystals. Int. J. Plast. 19 (9), 1401–1444. Clayton, J.D., McDowell, D.L., 2004. Homogenized finite elastoplasticity and damage: theory and computations. Mech. Mater. 36, 799–824. Collins, I.F., Houlsby, G.T., 1997. Application of thermomechanical principles to the modeling of geotechnical materials. Proc. Roy. Soc. Lond. A 453, 1975–2001. Cundall, P., 1976. Explicit finite-difference methods in geomechanics. Numerical Methods in Geomechanics 1, 132–150 (Proceedings of the EF conference on Numerical Methods in Geomechanics, Blacks burg, VA, 1976). Cundall, P.A., 1990. Numerical modelling of jointed and faulted rock. In: Mechanics of Jointed and Faulted Rock. A.A. Balkema, Rotterdam, pp. 11–18. FLAC, 1999. Itasca Consulting Group, Inc., Minneapolis, MN 55415, USA. Germain, P., Nguyen, Q.S., Suquet, P., 1983. Continuum thermodynamics. ASME J. Appl. Mech. 50, 1010–1020. Ghanem, R., Spanos, P.D., 1991. Stochastic Finite Elements: A Spectral Approach. Springer-Verlag, Berlin. Hashin, Z., 1983. Analysis of composite materials. ASME J. Appl. Mech. 50, 481–505.

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

1159

Hill, R., 1950. The Mathematical Theory of Plasticity. Clarendon Press, Oxford. Hill, R., 1963. Elastic properties of reinforced solids: some theoretical principles. J. Mech. Phys. Solids 11, 357–372. Hill, R., 1967. Essential structure of constitutive laws for metal composites and polycrystals. J. Mech. Phys. Solids 15, 79–95. Huet, C., 1990. Application of variational concepts to size effects in elastic heterogeneous bodies. J. Mech. Phys. Solids 38, 813–841. Huet, C., 1995. Bounds and hierarchies for the overall properties of viscoelastic heterogeneous and composite materials. Arch. Mech. 47 (6), 1125–1155. Huet, C., 1999. Coupled size and boundary-condition effects in viscoelastic heterogeneous and composite bodies. Mech. Mater. 31 (12), 787–829. Jeulin, D., 1993. Damage simulation in heterogeneous materials from geodesic propagations. Eng. Ccomput. 10, 81–91. Jeulin, D., Ostoja-Starzewski, M., 2000. Shear bands in elasto-plastic response of random composites by geodesics, 8th ASCE Specialty Conference on Probabilistic Mechanics and Structural Reliability, University of Notre Dame, 2000. Jiang, M., Ostoja-Starzewski, M., Jasiuk, I., 2001. Scale-dependent bounds on effective elastoplastic response of random composites. J. Mech. Phys. Solids 49 (3), 655–673. Kachanov, L.M., 1971. Foundations of the Theory of Plasticity. North-Holland, Amsterdam. Kanit, T., Forest, S., Galliet, I., Monoury, V., Jeulin, D., 2001. Determination of the size of the representative volume element for random composites: statistical and numerical approach. Int. J. Solids Struct. 40, 3647–3679. Ku, A.P.-D., Nordgren, R.P., 2001. On plastic collapse of media with random yield strength. ASME J. Appl. Mech. 68, 715–724. Mandel, J., 1972. Plasticite´ Classique et Viscoplasticite´, CISM Courses and Lectures 97. Springer, WienNew York. Maugin, G.A., 1999. The Thermomechanics of Nonlinear Irreversible Behaviors – an Introduction. World Scientific, Singapore. Mellor, P.B., 1982. Experimental studies of plastic anisotropy in sheet metal. In: Hopkins, H.G., Sewell, M.J. (Eds.), Mechanics of Solids – The Rodney Hill 60th Anniversary Volume. Pergamon Press, pp. 383–415. Nordgren, R.P., 1992. Limit analysis of a stochastically inhomogeneous plastic medium with application to plane contact problems. ASME J. Appl. Mech. 59, 477–484. Olszak, W., Rychlewski, J., Urbanowski, W., 1962. Plasticity under inhomogeneous conditions. Adv. Appl. Mech. 7, 132–214. Ostoja-Starzewski, M., 1992. Plastic flow of random media: microme-chanics, Markov property and sliplines, Appl. Mech. Rev. (Special Issue: Material Instabilities), Zbib, H.M., Shawki, T.G., Batra, R.C., (Eds.) 45(3, Part 2), S75–S82. Ostoja-Starzewski, M., 1998. Random field models of heterogeneous materials. Int. J. Solids Struct. 35 (19), 2429–2455. Ostoja-Starzewski, M., 2001a. On geometric acoustics in random, locally anisotropic media. Cont. Mech. Thermodyn. 13, 131–134. Ostoja-Starzewski, M., 2001. Mechanics of random materials: Stochas-tics, scale effects, and computation. In: Jeulin, D., Ostoja-Starzewski, M., (Eds.), Mechanics of Random and Multiscale Micro structures. 93-161, CISM Courses and Lectures 430, Springer, Wien-New York. Ostoja-Starzewski, M., 2002a. Microstructural randomness versus representative volume element in thermomechanics. ASME J. Appl. Mech. 69, 25–35. Ostoja-Starzewski, M., 2002b. Towards stochastic continuum thermodynamics. J. Non-Equilib. Thermodyn. 27 (4), 335–348. Ostoja-Starzewski, M., 2004. Fracture of brittle micro-beams. ASME J. Appl. Mech. 71, 424–427. Ostoja-Starzewski, M., Castro, J., 2003. Random formation, inelastic response, and scale effects in paper. Phil. Trans. Roy. Soc. Lond. A 361 (1806), 965–986. Ostoja-Starzewski, M., Ilies, H., 1996. The Cauchy and characteristic boundary value problems for weakly random rigid-perfectly plastic media. Int. J. Solids Struct. 33 (8), 1119–1136.

1160

M. Ostoja-Starzewski / International Journal of Plasticity 21 (2005) 1119–1160

Ostoja-Starzewski, M., Stahl, D.C., 2000. Random fiber networks and special elastic orthotropy of paper. J. Elast. 60 (2), 131–149. Pindera, M.-J., Aboudi, J., Arnold, S.M., 2004. Analysis of the spallation mechanism suppression in plasma-sprayed TBCs through the use of heterogeneous coat architectures, Int. J. Plasticity, this issue. Ponte Castan˜eda, P., Suquet, P.M., 1998. Nonlinear composites. Adv. Appl. Mech. 34, 171–302. Sachs, G., 1928. Zur Ableitung einer Fliessbedingung. Z. Ver. Deutsch. Ing. 72, 734–736. Sluzalec, A., 2004. Theory of Metal Forming Plasticity: Classical and Advanced Topics. Springer Verlag, Berlin. Spearot, D.E., Jacob, K.I., McDowell, D.L., 2004. Non-local separation constitutive laws for interfaces and their relation to nanoscale simulations. Mech. Mater. 36 (9), 825–847. Stoyan, D., Kendall, W.S., Mecke, J., 1987. Stochastic Geometry and its Applications. J. Wiley & Sons, New York. Suquet, P., 1986. Elements of homogenization for solid mechanics. In: Sanchez-Palencia, E., Zaoui, A. (Eds.), Homogenization Techniques for Composite Media, Lecture Notes in Physics 272, SpringerVerlag, Berlin, pp. 193–278. Suquet, P.M., 1993. Overall potentials and extremal surfaces of power law or ideally plastic composites. J. Mech. Phys. Solids 41 (6), 981–1002. Suquet, P., 1997. Effective properties of nonlinear composites. In: Suquet, P.M. (Ed.), Continuum Micromechanics, CISM Courses and Lectures 377, Springer-Verlag, Berlin, pp. 197–264. Taylor, G.I., 1938. Plastic strain in metals. J. Inst. Met. 62, 307–324. Terada, K., Hori, M., Kyoya, T., Kikuchi, N., 2000. Simulation of the multi-scale convergence in computational homogenization approaches. Intl. J. Solids Struct. 37, 2285–2311. Wang, Z.-P., 1997. Void-containing nonlinear materials subject to high-rate loading. J. Appl. Phys. 81 (11), 7213–7227. Wang, Z.-P., Sun, C.T., 2003. Modeling micro-inertia in heterogeneous materials under dynamic loading. Wave Motion 36, 473–485. Ziegler, H., 1983. An Introduction to Thermomechanics. North-Holland, Amsterdam. Ziman, J.M., 1979. Models of Disorder. Cambridge University Press, Cambridge.