Computational micro-macro material testing

Vol. 8, 2, 131-228 (2001)

Archives of Computational Methods in Engineering State of the art reviews

Computational Micro-macro Material Testing T.I. Zohdi and P. Wriggers Institut f¨ ur Baumechanik und Numerische Mechanik (IBNM) Universit¨ at Hannover Appelstraße 9A, 30167 Hannover, Germany

Dedicated to Christian Huet and J. Tinsley Oden on the occasion of their 65th birthdays. 1 INTRODUCTION A key to the success of many modern structural components is the tailored behavior of the material. A relatively inexpensive way to obtain macroscopically desired responses is to enhance a base matrix properties by the addition of microscopic matter, i.e. to manipulate the microstructure. Accordingly, in many modern engineering designs, materials with highly complex microstructures are now in use. The macroscopic characteristics of modified base materials are the aggregate response of an assemblage of different “pure” components, for example several particles or fibers suspended in a binding matrix material (Figure 1.1). In the construction of such materials, the basic philosophy is to select material combinations to produce aggregate responses possessing desirable properties from each component. For example, in structural engineering applications, the classical choice is a harder particulate phase that serves as a stiffening agent for the base matrix material. Such inhomogeneities are encountered in metal matrix composites, concrete, etc. A variety of materials are characterized by particulate inhomogeneities as shown in Figures 1.2 and 1.3. BASE MATERIAL MACROSCOPIC STRUCTURE

ADDITIVES

NEW "MATERIAL"

MIXING PROPELLOR Figure 1.1. LEFT: Doping a base material with particulate additives. RIGHT: Typical loose particulate additives

If one were to attempt to perform a direct numerical simulation, for example of the mechanical response of a macroscopic engineering structure composed of a microheterogeneous material, incorporating all of the microscale details, an extremely fine spatial discretization mesh would be needed to capture the effects of the microscale heterogeneities. The resulting system of equations would contain literally billions of numerical unknowns. Such problems c 2001 by CIMNE, Barcelona (Spain).

ISSN: 1134–3060

Received: February 2001

132

T.I. Zohdi and P. Wriggers

Figure 1.2. Left to right and top to bottom: (1) Ductile iron (×100), (2) Titanium-aluminum alloy (×450), (3) Aluminum alloy (×250), (4) Zinc alloy (×1000). For micrographs (1)-(4) Source: The American Society for Metals [60]

Figure 1.3. LEFT: Pearlite ductile iron with matrix etched away to show secondary graphite, and bull’s eye ferrite, around primary graphite nodules (×475). RIGHT: Ductile iron, graphite nodules (spherulites) encased in envelopes of free ferrite, all in a matrix of pearlite (×100). Source: The American Society for Metals [60]

are beyond the capacity of computing machines for the forseeable future. Furthermore, the exact subsurface geometry is virtually impossible to ascertain throughout the structure. In addition, even if one could solve such a system, the amount of information to process would be of such complexity that it would be difficult to extract any useful information on the desired macroscopic behavior. It is important to realize that solutions to partial differential equations, of even linear material models, at infinitesimal strains, of small bodies containing a few heterogeneities are still open problems. In short, complete solutions are virtually impossible. Because of these facts, the use of regularized or homogenized material models (resulting in smooth coefficients in the partial differential equations) are commonplace in

133

Computational Micro-macro Material Testing

virtually all branches of physical sciences. The usual approach is the compute a constitutive “relation between averages”, relating volume averaged field variables. Thereafter, the regularized properties can be used in a macroscopic analysis (Figure 1.4). The volume averaging takes place over a statistically representative sample of material, referred to in the literature as a representative volume element (RVE). The internal fields to be volumetrically averaged must be computed by solving a series boundary value problems with test loadings. Such regularization processes are referred to as “homogenization”, “mean field theories”, “theories of effective properties”, etc. For overviews, we refer the interested reader to Jikov et al. [47] for mathematical reviews or to Aboudi [1], Mura [63] or Nemat-Nasser and Hori [69] for mechanically inclined accounts of the subject. ACTUAL STRUCTURE

STRUCTURE WITH EFFECTIVE PROPERTIES

EFFECTIVE PROPERTIES DETERMINED

Ω Figure 1.4. The use of effective properties in structural engineering

For a sample to be statistically representative it must usually contain a great deal of heterogeneities (Figure 1.5), and therefore the computations over the RVE are still extremely large, but are of reduced computational effort in comparison with a direct attack on the “real” problem. Historically, most classical analytical methods for estimating the macroscopic response of such engineering materials have a strongly phenomenological basis, and are in reality non-predictive of material responses that are unknown a-priori. This is true even in the linearly elastic, infinitesimal strain, range. In plain words, such models require extensive experimental data to “tune” parameters that have little or no physical significance. Criticisms, such as the one stated, have led to computational approaches which require relatively simple descriptions on the microscale, containing parameters that are physically meaningful. In other words, the phenomenological aspects of the material modeling are reduced, with the burden of the work being shifted to high performance computational methods. Stated clearly, the “mission” of computational micro-macro mechanics is to determine relationships between the microstructure and the macroscopic response or “structural property” of a material, using models on the microscale that are as simple as possible.∗ ∗ Appropriately, the following quote, attributed to Albert Einstein: “A model should be as simple as possible, but not simpler”, is a guiding principle.

134


L APPROXIMATELY UNIFORM B.C’s

Ω

Ω b r b κ0 , µ > µ0 , and 0 u) is a minimum it is equal to Π(ˆ u) is an absolute maximum if κ < κ , µ < µ0 , (4) When Π(ˆ the strain energy, (5) When IE 0 is vanishingly small compared to IE, then the variational principle collapses to the PMCPE and (6) When IE 0 is infinitely large compared to IE, then the variational principle collapses to the PMPE. Similar properties hold for a pure traction loaded problem. The relations are then extremized by optimizing with respect to κ0 and 2 ∗ µ0 . Noting that when u|∂Ω = E · x the strain energy is W = 12 (9κ∗ ( trE 3 ) + 2µ E : E ), one can arrive at, for a two-phase microstructure:

κ1 +

v2 1 κ2 −κ1

+

bulk modulus H−S lower

µ1 +

≤ κ∗ ≤

3(1−v2 ) 3κ1 +4µ1

+

6(1−v2 )(κ1 +2µ1 ) 5µ1 (3κ1 +4µ1 )

shear modulus H−S lower

def bound = µ∗,−

1 − v2 2 + 3κ23v +4µ2

1 κ1 −κ2

,

def

bulk modulus H−S upper bound = κ∗,+

def bound = κ∗,−

v2 1 µ2 −µ1

κ2 +

≤

µ∗

≤

µ2 +

(2.15)

(1 − v2 ) 1 µ1 −µ2

+

6v2 (κ2 +2µ2 ) 5µ2 (3κ2 +4µ2 )

shear modulus H−S upper

,

def bound = µ∗,+

where κ1 , µ1 and κ2 , µ2 are the bulk and shear moduli for the phases, while v2 is the phase 2 volume fraction. The proofs, which are algebraically complicated, can be found in Hashin and Shtrikman [29] and [30]. We emphasize that in the derivation of the bounds, the body is assumed to be infinite, the microstructure isotropic, and that the effective responses are isotropic. We remark that while the bounds are the tightest possible, under the previous assumptions, when no geometric information is included. 2.4.2 The concentration tensor

Most classical approximation methods for particulate microstructural analyses start with 1 ( Ω1 dΩ + Ω2 dΩ) = v1 Ω1 + v2 Ω2 and (2) the following definitions (1) Ω = |Ω|

142

σΩ =


1 |Ω| ( Ω1

σ dΩ +

Ω2

σ dΩ) = v1 σΩ1 + v2 σΩ2 . By direct manipulation we obtain

IE ∗ : Ω = = = = =

v1 σΩ1 + v2 σΩ2 v1 IE 1 : Ω1 + v2 IE 2 : Ω2 IE 1 : (Ω − v2 Ω2 ) + v2 IE 2 : Ω2 IE 1 : + v2 (IE 2 − IE 1 ) : Ω2 (IE 1 + v2 (IE 2 − IE 1 ) : C) : Ω ,

=IE

∗

(2.16)

def

where quantity C : Ω = Ω2 , is known as the concentration tensor. One can also recast ˆ these expressions to yield the following form, (IE ∗ )−1 = (IE 1 )−1 +v2 ((IE 2 )−1 −(IE 1 )−1 ) : C), def where Cˆ : σΩ = σΩ2 . In either form there has been no approximation yet. The “burden” in the computations has shifted to the determination of C. Classical methods approximate C. For example, the simplest approximation is C = I, which is the Voigt approximation, IE ∗ = (IE 1 + v2 (IE 2 − IE 1 ) : I) or for the Reuss approximation (IE ∗ )−1 = (IE 1 )−1 + v2 ((IE 2 )−1 − (IE 1 )−1 ) : I), where Cˆ = I. 2.4.3 The Eshelby result

In 1957, Eshelby [25] developed an elegant formalism to determine the solution to the elasticity solution of a single inclusion embedded in an infinite matrix of material with uniform exterior loading. While alone this result is of little practical interest, this solution, which is relatively compact, has been the basis of approximation methods for effective properties for non-interacting and weakly interacting particulate solutions. The essentials of the Eshelby result are as follows. Consider a single linearly elastic particle, at infinitesimal strains, embedded in an infinite matrix: If the shape of the particle is ellipsoidal, under uniform far field stress or strain, then the stress, and hence the strain, in the particle is constant. As we have mentioned, this result is the foundation of most methods of approximation of the effective property, and is used to approximate C for ellipsoidal geometries. Ellipsoidal shapes are qualitatively useful since the geometry can mimic a variety of microstructures (1) platelet behavior when the ellipsoid is oblate with a high aspect ratio (2) needlelike microstructure when it is prolate with high aspect ratio, (3) cracklike behavior as the material in the particle is modeled as extremely soft with very high aspect ratios for the particle geometry and (4) pores when the particle shapes are spheres with very soft material inside. For a single particle in three dimensions (Figure 2.3) we have, ∀x ∈ Ω2 , IE 2 : = IE 1 : (E + ∆ − T ) = IE 1 : (E + P : T − T ) = IE 1 : E + IE 1 : (P − I) : T , (2.17) where E is the far field uniform strain, = E + ∆ is the strain in the particle, T is the def transformation strain and P is the fourth-order Eshelby tensor which satisfies P : T = ∆. The Eshelby result, stated in a compact way, is: P and T are constants in an isolated ellipsoidal particle under uniform far field loading. 2.4.4 Dilute methods

The use of the Eshelby result is straightforward as a method of determining approximate effective properties for regularly shaped particles. When the particle is a sphere, P is isotropic and the relations are relatively easy to write down. Breaking the tensor relations into dilatational and deviatoric pieces yields, ∀x ∈ Ω2 one has tr3 = α tr3T + trE 3 , α =

143


UNIFORM FAR FIELD IMPOSED

Ω2

Ω1

Figure 2.3. The notation for the Eshelby formalism 1+ν1 3κ1 2 4−5ν1 3(1−ν1 ) = 3κ1 +4µ1 , = βT + E and β = 15 1−ν1 . The constitutive relations are, ∀x ∈ Ω2 , 3κ−2µ tr σ tr tr T tr 3 = 3κ2 3 = 3κ1 ( 3 − 3 ) and σ = 2µ2 = 2µ1 ( −T ), where ν = 6κ+2µ . Combining def the previous expressions yields, ∀x ∈ Ω2 , tr = γtrE, = ρE , γ = ακ2 +κκ11(1−α) and def ρ = βµ2 +µµ11(1−β) . The approximations allow the determination of the strains in the particles

as a function of the loading and geometry. Notice that the deviatoric and dilatational strain components in the particle are the far field corresponding strain components multiplied by a factor composed by the corresponding components of the matrix material divided by a convex combination of the matrix and particle materials. The transformation relations are tabulated for various simple particulate shapes in Mura [63]. It is important to realize that the stress fields in the matrix are not constant, and are extremely complicated. Fortunately, in developing approximate expressions for the effective response it unnecessary to know the fields in the matrix. Consequently, the concentration tensor, under the non-interacting particle assumption, is isotropic and can be written as  γ 4ρ 3 + 32 γ  − 2ρ  3 32 def  γ − 2 ρ 3 32 C =   0   0

γ 3 γ 3 γ 3

0

− 23 ρ2 + 43 ρ2 − 23 ρ2 0 0 0

γ 3 γ 3 γ 3

− 23 ρ2 − 23 ρ2 + 43 ρ2 0 0 0

0 0 0 ρ 2

0 0

0 0 0 0 ρ 2

0

0 0 0 0 0

    .   

(2.18)

ρ 2

Therefore, in the case of noninteracting spheres, C is isotropic with deviatoric and dilatational components given by γ and ρ, which are furnished by the Eshelby solution. This is the framework for the Dilute method, which postulates the complete non-interaction between microstructural components. Therefore we have IE ∗ ≈ IE 1 + v2 (IE 2 − IE 1 ) : C implies κ∗ ≈ κ1 + v2 (κ2 − κ1 )γ

and

µ∗ ≈ µ1 + v2 (µ2 − µ1 )ρ.

Dilute approximations are usually only accurate at extremely low volume fractions.

(2.19)

144


2.4.5 The Mori-Tanaka method

Non-interaction of particulates is an unrealistic assumption for materials with randomly dispersed particulate microstructure at even a few percent volume fraction. Weak interaction between particles, via a sensitivity to increased volume fraction, involves slight modification of the dilute method. This involves assuming that there exists another tensor, such that Ω2 = A : Ω , thus implying IE ∗ : Ω = (IE 1 + v2 (IE 2 − IE 1 ) : A) : Ω . The tensor A = A(v2 ) must be approximated. There are two requirements: (1) limv2 →0 A = C and (2) limv2 →1 A = I. The most widely used approach is that of Mori-Tanaka [62]. The Mori-Tanaka method proceeds by constructing a tensor G such that G : Ω1 = Ω2 , thus implying that Ω = v1 Ω1 +v2 Ω2 = v1 G−1 : Ω2 +v2 Ω2 = (v1 G−1 +v2 I) : Ω2 = A−1 : Ω2 . The Mori-Tanaka method approximates G by C. The physical implications are that the particle “sees” a matrix material with an average stress state provided by the average in the matrix. One can consider the Mori-Tanaka as one possible curve fit that satisfies the two mentioned conditions, i.e. A = I when v2 = 1 and A = C when v2 = 0 (Figure 2.4). The result is def

κ∗ = κ1 + v2 (κ2 − κ1 )η def

and

def

µ∗ = µ1 + v2 (µ2 − µ1 )ζ,

(2.20)

def

1 (9κ1 +8µ1 ) , where where η = Θ+(1−v2Θ)(κ2−κ1 ) , ζ = Φ+(1−v2Φ)(µ2−µ1 ) , Θ = κ1 + 43 µ1 , Φ = µ1 + µ6(κ 1 +2µ1 ) ∀x ∈ Ω2 , tr = ηtrE and = ζE . The interaction, via a volume fraction dependency, is explicitly clear in the expressions for η and ζ. It is noted that for spherical particles the Mori-Tanaka matches the well known Hashin-Shtrikman lower bound [30] on effective properties.

I

A

C

V_2=0

V_2=1

VOLUME FRACTION Figure 2.4. The Mori-Tanaka “extrapolation”

Remark: Another approach to incorporate particulate interaction is the Self Consistent method. The idea is simply that the particles “see” the effective medium instead of the matrix in the computation of the transformation strain. In other words the “matrix material’ in the Eshelby tensor is simply the effective medium. Therefore one simply replaces the matrix material in the concentration formulae (only in C) with the effective property to obtain ∗

κ∗ = κ1 (1 − v2 γ α−1 ∗ )

and

∗

µ∗ = µ1 (1 − v2 ρ β−1 ∗ ),

(2.21)

145


240

220

EFFECTIVE BULK MODULUS K*

200

180

REUSS VOIGT HASH/SHTRHASH/SHTR+ 160 DILUTE MORI-TANAKA 140

120

100

80

60 0

0.2

0.4

0.6 VOLUME FRACTION

0.8

1

1.2

0.4

0.6 VOLUME FRACTION

0.8

1

1.2

180

EFFECTIVE SHEAR MODULUS MU*

160

REUSS VOIGT HASH/SHTRHASH/SHTR+ DILUTE 140 MORI-TANAKA 120

100

80

60

40

20 0

0.2

Figure 2.5. The predicted effective bulk (κ∗ ) and shear moduli (µ∗ ) for an Aluminum matrix (κ = 77.9 GP a, µ = 24.9 GP a)/Boron particle (κ = 230 GP a, µ = 172 GP a) ∗

∗

∗

1+ν 3κ ∗ = 2 4−5ν . The approach is identical to the where α∗ = 3(1−ν ∗ ) = 3κ∗ +4µ∗ and where β 15 1−ν ∗ dilute approximation, with the result being somewhat more complicated expressions for the approximate effective property, however, which can be computed in closed form in many cases. We refer the reader to the books of Mura [63], Aboudi [1] or Nemat-Nasser and Hori [69]. The self-consistent method produces negative effective bulk and shear responses for voids above volume fractions of 50 %. For rigid inclusions it produces infinite effective bulk responses for any volume fraction and infinite effective shear responses above 40 %. For proofs and discussions on variations of this method, see Aboudi[1]. With these facts in mind, one can safely employ such methods only at very low volume fraction levels.

146


3 A-PRIORI INTERNAL FIELD BOUNDS During material design development, when selecting micro-additives, information about the magnitude of the resulting distortion of the otherwise smooth internal fields corresponding to the matrix material alone, is valuable to characterize a new composite material’s behavior. Such information would be attractive, especially if it involved no intensive numerical computations. Accordingly, in the analysis to follow, a-priori bounds are developed on the difference between solutions to two boundary value problems involving linearly elastic microstructure at infinitesimal strains. Both boundary value problems have the same external geometry, boundary conditions and loading. However, one possesses a uniform microstructure, while the other is “perturbed” with inhomogeneities (Figure 3.1). As will be shown, the bounds developed are solely in terms the solution to the uniform microstructure problem, and the given loading data. There are no unknown constants in the bounds. Furthermore, no calculation of computationally expensive heterogeneous problems are necessary, and thus the bounds can serve as a-priori indicators, via their computed magnitude, of whether the effects of the particles are great enough to warrant an intensive computational analysis. The bounds developed make no assumptions on the character of the microstructure of the heterogeneous problem, other than it be pointwise positive definite, as well as no assumptions on the external loading or geometry. These bounds are a result of work found in Zohdi [100], and are closely related to results also found in Zohdi and Wriggers [93] [94], Zohdi, et al. [91] and Oden and Zohdi [70].

Ω

HOMOGENEOUS PROBLEM (ONLY MATRIX MATERIAL)

Ω

INHOMOGENEOUS PROBLEM (MATRIX AND INHOMOGENEITIES)

Figure 3.1. The two problems considered. LEFT, the homogeneous problem, and RIGHT, the inhomogeneous problem. Both problems have the same external geometry and loading

3.1 Governing Equations We consider a structure which occupies an open bounded domain in Ω ∈ R3 . Its boundary is denoted ∂Ω. The body is in static equilibrium under the action of body forces, f , and surface tractions, t. The boundary ∂Ω = Γu ∪ Γt consists of a part Γu and a part Γt on which displacements and tractions are respectively prescribed. The data are assumed to be such that f ∈ L2 (Ω) and t ∈ L2 (Γt ), but less smooth data can be considered without complications. The mechanical properties of the material are characterized by a (regular) spatially constant inhomogeneity-free elasticity tensor, IR, and a spatially nonconstant

147


(perturbed) elasticity tensor, IE(x), both of which are in R3 ×3 , and whose components satisfy the following symmetry and boundedness conditions: ∀ ∈ R3×3 , = T , α+,R : ≥ : IR : ≥ α−,R : , ∞ > α+,R , α−,R > 0, where Rijkl = Rjikl = Rijlk = Rklij , 1 ≤ i, j, k, l ≤ 3, and α+ : ≥ : IE(x) : ≥ α− : , ∞ > α+ , α− > 0, where Eijkl (x) = Ejikl (x) = Eijlk (x) = Eklij (x), 1 ≤ i, j, k, l ≤ 3. The quantities Rijkl and Eijkl (x) are the Cartesian components of IR and IE(x), respectively. The microstructure is assumed to be perfectly bonded. Following standard notation, we denote H 1 (Ω) as the usual (Sobolev) space of functions with generalized partial derivatives of order ≤ 1 in L2 (Ω). def We define H 1 (Ω) = [H 1 (Ω)]3 as the space of vector-valued functions whose components def have generalized partial derivatives of order ≤ 1 in L2 (Ω) = [L2 (Ω)]3 . 2

2

3.1.1 Homogeneous-coefficient boundary-value problem

The solution to the constant coefficient problem, denoted the regular solution, uR , is characterized by a virtual work formulation:

Find a uR ∈ H 1 (Ω), uR |Γu = d, such that Ω

∇v : σR dΩ =

f · v dΩ +

Ω

def

Γt

t · v dA

∀v ∈ H 1 (Ω), v|Γu = 0,

(3.1)

def

= BR (uR ,v )

= F (v )

where σR = IR : ∇uR . The equivalent complementary form is Find σ R , ∇ · σR + f = 0, σ R · n|Γt = t such that

Ω

τ : IR−1 : σ R dΩ =

def

Γu

τ · n · d dA

∀τ , ∇ · τ = 0, τ · n|Γt = 0.

(3.2)

def

= AR (σ R ,τ )

= T (τ )

For the complementary problem, similar restrictions are placed on the solution and test fields to force the integrals to make sense. In other words, we assume that solutions produce finite global energy. 3.1.2 Inhomogeneous-coefficient boundary-value problem

The solution corresponding to a material with microstructure is u, and is characterized by the following virtual work formulation: Find u ∈ H 1 (Ω), u|Γu = d, such that

Ω

∇v : σ dΩ =

def

= B(u,v )

Ω

f · v dΩ + def

Γt

= F (v )

t · v dA

∀v ∈ H 1 (Ω), v|Γu = 0.

(3.3)

148


In the infinitesimal strain linearly elastic case, then σ = IE : ∇u. The equivalent complementary form is Find σ, ∇ · σ + f = 0, σ · n|Γt = t such that

Ω

τ : IE

−1

: σ dΩ =

def

Γu

τ · n · d dA

∀τ , ∇ · τ = 0, τ · n|Γt = 0.

(3.4)

def

= A(σ ,τ )

= T (τ )

3.2 Homogeneous/inhomogeneous Solution Difference Bounds Our immediate objective is to bound ||u − uR ||E(Ω) and ||σ − σ R ||E −1 (Ω) in terms of IR, IE, uR and σ R . Accordingly, we seek bounds on the energy norm of the solution difference, H−,(u) ≤ ||u − uR ||E(Ω) ≤ H+,(u) , and for the complementary norm H−,(σ) ≤ ||σ − σR ||E −1 (Ω) ≤ H+,(σ) . 3.2.1 Primal and complementary upper bounds From Boxes (3.1) and (3.3), we have Ω ∇v : IE : ∇u dΩ = Ω ∇v : IR : ∇uR dΩ = R dΩ from both sides yields f · v dΩ + t · v dA. Subtracting ∇v : I E : ∇u Ω Γt Ω Ω ∇v : R R R

IE : ∇(u − u ) dΩ = Ω ∇v : (IR : ∇u − IE : ∇u ) dΩ. Since v is an arbitrary virtual displacement, we may set v = u − uR to yield ||u −

uR ||2E(Ω)

∇(u − uR ) : (IR : ∇uR − IE : ∇uR ) dΩ

= Ω

= Ω

(IE 2 : ∇(u − uR )) : IE 2 : IE −1 : ((IR − IE) : ∇uR ) dΩ 1

≤

1 Ω

∇(u − uR ) : IE : ∇(u − uR ) dΩ

2

(3.5)

||u−uR ||E(Ω)

×

1

Ω

1

((IR − IE) : ∇uR ) : IE −1 : ((IR − IE) : ∇uR ) dΩ

2

.

Therefore, a primal solution difference upper bound is ||u − uR ||2E(Ω) ≤

Ω

((IR − IE) : ∇uR ) : IE −1 : ((IR − IE) : ∇uR ) dΩ.

(3.6)

def

= (H+,(u))2

To derive this result we used the fact that a symmetric positive-definite matrix has a unique 1 1 square root, IE = IE 2 : IE 2 and the Cauchy-Schwarz inequality. Note that the H+,(u)bound depends only on IR, IE and uR . The process is repeated in a similar manner to form a complementary upper bound: ||σ −

σR ||2E −1 (Ω)

≤

Ω

((IR−1 − IE −1 ) : σ R ) : IE : ((IR−1 − IE −1 ) : σR ) dΩ .

def

= (H+,(σ) )2

Analogous to the H+,(u) -bound, the H+,(σ) -bound depends only on IR, IE and σR .

(3.7)

149


3.2.2 Exact hybrid bounds

If we choose w = uR , which is a kinematically admissible function, we obtain ||u − uR ||2E(Ω) = 2(J (uR ) − J (u)) which implies J (u) = J (uR ) − 12 ||u − uR ||2E(Ω) . Also, choosing γ = σ R , which is statically admissible, we have ||σ − σR ||2E −1 (Ω) = 2(K(σR ) − K(σ)) which implies K(σ) = K(σ R ) − 12 ||σ − σ R ||2E −1 (Ω) . A so-called “exact” bound on a hybrid norm, which is a combination of the primal and complementary norms, can be obtained by combining the two previous results to yield def

(P(uR , σ R ))2 = 2(J (uR ) + K(σR )) = ||u − uR ||2E(Ω) + ||σ − σR ||2E −1 (Ω)

HYBRID BOUND

+ 2(J (u) + K(σ)),

(3.8)

=0

= ||u − uR ||2E(Ω) + ||σ − σR ||2E −1 (Ω) . This “estimate” is exact for the hybrid homogeneous/inhomogeneous measure of solution difference. 3.2.3 Lower primal and complementary bounds

From Equation (3.8), lower primal and complementary bounds follow, ||u − uR ||2E(Ω) = 2(J (uR ) + K(σ R )) − ||σ − σR ||2E −1 (Ω) , ≥ 2(J (uR ) + K(σ R )) − (H+,(σ) )2 = (H−,(u))2 , def

(3.9) ||σ − σ R ||2E −1 (Ω) = 2(J (uR ) + K(σ R )) − ||u − uR ||2E(Ω) , def

≥ 2(J (uR ) + K(σ R )) − (H+,(u) )2 = (H−,(σ))2 . Therefore, in summary we have (H−,(u))2 ≤ 2(J (uR ) − J (u)) ≤ (H+,(u))2 ,

||u − uR ||2E(Ω)

SOLUTION DIFFERENCE

(H−,(σ) )2 ≤ 2(K(σ R ) − K(σ)) ≤ (H+,(σ) )2 ,

(3.10)

||σ − σ R ||2E −1 (Ω)

SOLUTION DIFFERENCE

(P(uR , σ R ))2 = ||u − uR ||2E(Ω) + ||σ − σ R ||2E −1 (Ω) . These relations hold for any loading, external geometry and pointwise positive-definite microstructure.

150


3.2.4 Observations

If we return to the special case of a microstructure composed of a matrix embedded with particulate matter, the expressions attain very compact forms. Let us write the heterogeneous material as IE = IR + ∆IR. Referring to Figure 3.1, it is clear that ∆IR = IE − IR = 0 in the matrix. Therefore, the upper bounds collapse to †

||u − uR ||2E(Ω) ≤

N i=i

||σ − σR ||2E −1 (Ω) ≤

Ωi

(∆IR : ∇uR ) : IE −1 : (∆IR : ∇uR ) dΩi , (3.11)

N i=i

Ωi

(∆IR−1 : σ R ) : IE : (∆IR−1 : σ R ) dΩi ,

where Ωi are the domains of the (N) inhomogeneities. In other words, the difference can be computed solely in terms of the product of the difference between the elastic properties of the inhomogeneities and the matrix, and the unperturbed homogeneous elastic field. Furthermore, this only needs to be computed over the domains of the inhomogeneities. Similar expressions can be obtained for the hybrid norm, and the lower bounds. 3.2.5 A quantitative example

Let us further consider the case where IE = IR + ∆IR = ξIR, where 1 < ξ < ∞ in the particle domains. Also we consider an applied homogeneous boundary deformation of the form 









u1 |∂Ω E11 E12 E13 x1  u2 |∂Ω  =  E12 E22 E23   x2  , u3 |∂Ω E31 E32 E33 x3

(3.12)

E

where the Eij are interpreted as applied strains, via linear displacements at infinitesimal strains. For this loading case, with f = 0, after some simple algebra, we have

(H

+,(u) 2

) = (H

+,(σ) 2

R

R

2

) = (P(u , σ )) =

(1 − ξ)2 vp E : IR : E|Ω|, ξ

(3.13)

where vp is the particle volume fraction, and where (H−,(u))2 = (H−,(σ))2 = (P(uR , σR ))2 − (H+,(u))2 = 0. We can normalize these expressions by the strain energy of the homogeneous system, ||uR ||2R(Ω) = E : IR : E|Ω|, to obtain ||u − uR ||2E(Ω) (1 − ξ)2 vp = ≤ R 2

||u ||R(Ω)

ξ

(||u − uR ||2E(Ω) + ||σ − σ R ||2E −1 (Ω) ) .

||uR ||2R(Ω)

(3.14)

The results are shown in Figure 3.2 for variation in the volume fraction and mismatch in the material properties. The iso-contour values are also shown. We remark that the lower bounds, which are zero in the special loading case used, could be positive or negative †

Note that IE −1 = IR−1 + ∆IR −1 (∆IR −1 = (∆IR)−1 ).

151


in other cases. As Figure 3.2 indicates, the normed solution difference can be quite high for moderate volume fractions and material mismatches. For example, for a mismatch of ξ = 2.5 and volume fraction of vp = 0.125, the normed solution difference is approximately 25.0 %, while for ξ = 5 and volume fraction of vp = 0.25, the normed solution difference is approximately 89.4 %, and for a mismatch of ξ = 10 and volume fraction of vp = 0.5, the normed solution difference is approximately 201.5 %.

2.5 2 1.5 1 0.5 0

1

2

3

4 5 6 7 ADDITIVE MISMATCH

8

9

10

0

0.5 0.45 0.4 0.35 0.3 0.25 0.2 VOLUME FRACTION 0.15 0.1 0.05

Figure 3.2. The upper and hybrid (exact) bound behavior. The plot is the same for both in this case

Remarks: Throughout the analysis, IR was assumed to be the homogeneous matrix material, however, as illustrated by the proofs, this is unnecessary. The tensor IR can also be a spatially-variable function, with no alteration in the analysis. Other choices for IR are exploited later in the monograph. The main point to be gathered from this chapter is that at relatively moderate volume fractions and phase contrasts, the internal fields are quite perturbed throughout the body due to the presence of the microheterogeneities. Consequently, for a reliable evaluation of a material’s behavior, particularly in the evaluation of progressive material failure, numerical simulations are virtually unavoidable. 4 COMPUTATIONAL/STATISTICAL TESTING METHODS The amount of uncertainty in the effective property bounds, as well the magnitude of the internal solution difference bounds, leads one to consider direct numerical simulation, whereby the effective responses can be obtained by volumetrically averaging (post processing) numerical solutions of boundary value problems representing the response of aggregate assemblages of heterogeneous material. Accordingly, in this chapter we investigate topics related to the numerical simulation of the testing of mechanical responses of samples of microheterogeneous solid materials formed by aggregates of particulates suspended in a binding matrix. Consistent with what is produced in associated manufacturing methods, the microstructures considered are generated by randomly distributing aggregates of particulate material throughout an otherwise homogeneous matrix material. Therefore, the resulting microstructures are irregular and nonperiodic. A primary issue in testing such materials is the fact that only finite sized samples can be tested, leading to no single response, but a distribution of responses. A technique employing potential energy principles

152


is developed to interpret the results of testing groups of samples. Three dimensional numerical examples employing the finite element method are given to illustrate the overall analysis and computational testing process. Afterwards, in the chapters to follow, a variety of nonlinear problems treated, building on the techniques developed in this chapter. The basis for this chapter follows from work found in Zohdi and Wriggers [97]. 4.1 A Boundary Value Formulation We consider a sample of heterogeneous material (Figure 4.1), with domain Ω, under a given set of specified boundary loadings. The variational boundary value problem is Find u ∈ H 1 (Ω), u|Γu = d, such that Ω

∇v : IE : ∇u dΩ =

Ω

f · v dΩ +

Γt

t · v dA

∀v ∈ H 1 (Ω), v|Γu = 0.

(4.1)

When we perform material tests satisfying Hill’s condition, we have, in the case of displacement controlled tests (loading case (1)) Γu = ∂Ω and u|∂Ω = E · x or for traction controlled tests (case (2)) Γt = ∂Ω and t|∂Ω = L · n. In either case we consider f = 0 and that the material is perfectly bonded. However, we note that (1) Hill’s condition is satisfied with f = 0 (and no debonding) and in case (2) it is satisfied even with debonding, however only if f = 0. The boundary value problem in Box 4.1 must be solved for each new sample, each possessing a different random microstructure (IE). The solution is then post processed for the effective quantities. It is convenient to consider the RVE domain Ω as a cube, and we shall do so for the remainder of the work.

X2

X 1 X3 Figure 4.1. A cubical sample of microheterogeneous material

4.2 Numerical Discretization In order to computationally simulate effective responses, our choice for spatial discretization is the finite element method. There are essentially two choices to mesh the microstructure with the finite element method, an unaligned or an aligned approach. We refer to an


153

unaligned approach as one which does not require the finite element boundaries to coincide with material interfaces when meshing the internal geometry (Figure 4.2). This leads to material discontinuities within the finite elements. An aligned approach would impose that the element boundaries coincide with material interfaces and therefore the elements have no material discontinuities within them. There are advantages and disadvantages to both approaches. Unaligned meshing has the advantages of rapid generation of structured internal meshes and consequently no finite element distortion arising from the microstructure. This is critical to computational performance if iterative solvers are to be used. The aligned meshing usually will require less finite elements than the unaligned approach for the same pointwise accuracy. However, the disadvantages are the mesh generation for irregular microstructures in three dimensions. Even if such microstructures can be meshed in an aligned manner, the finite element distortion leads to stiffness matrix ill conditioning and possible element instability (shape nonconvexity). For numerical studies comparing the meshing approaches, see Zohdi et at. [92]. Here our emphasis is on admitting general irregular microstructures, and rapidly evaluating them during the testing process, thus we have taken the unaligned approach, which we discuss in more detail next.

Figure 4.2. a) Unaligned meshing with material discontinuities within an element

4.2.1 Topological error

Inherent in the unaligned approach is the integration of discontinuous integrands (Figure 4.3). The topology is not embedded into the finite element a-priori, as it would be in an aligned approach, via isoparametric maps onto material interfaces. To some extent if the elements are much smaller than the particle length scales, the topology will be approximately captured. However, one can improve this representation (Figure 4.2). Since the finite element method is an integral-based method, the quadrature rules can be increased in an element by element fashion to better capture the geometry in elements with material discontinuities. A primary question is: If there is an integrand discontinuity in a finite element, how high should the quadrature rules be to perform accurate integration?

154


Topological error bounds for unaligned meshing We first consider a one dimensional reference element as shown in Figure 4.3. We assume that the elements are small compared to the length scales of the particulate matter (this is automatically dictated by a subspatial approximation error estimator presented in the next section) and as a consequence we assume that there is at most one discontinuity within the element. However, this assumption makes no difference in practice. We consider integrands, F (ζ), defined on the reference element, with jump discontinuity at δ, which admit the following decomposition F (ζ) =

C(ζ) +|]F (δ)[| H(ζ − δ),

continuous

(4.2)

jump

def

where ζ is the local coordinate, and |](·)[| = (·)|+ − (·)|− is the jump operator. Integrating 1 1 over the reference element, we have −1 F (ζ)dζ = −1 (C(ζ) + |]F (δ)[|H(ζ − δ)) dζ. We perform a straightforward Gauss-Legendre quadrature, with G quadrature points, where m points lie before the discontinuity at δ (Figure 4.3). We assume that the continuous function, possibly a polynomial, can be integrated exactly, or nearly exactly, with standard 1 G quadrature −1 C(ζ)dζ ≈ i=1 C(ζi )wi , where wi are the Gauss weights and where ζi 1 are the Gauss point locations. The remainder is the discontinuous part |]F (δ)[| −1 H(ζ − G m δ)dζ = |]F (δ)[|(1 − δ). We may write |]F (δ)[| i=1 H(ζi − δ)wi = |]F (δ)[| i=1 0 × wi + |]F (δ)[| G i=m+1 1 × wi . As a consequence, the integration error is

G G −1 F (ζ) dζ − (C(ζi ) + |]F (δ)[|H(ζi − δ))wi = ||]F (δ)[|| (2 − δ) − (2 − wi ) , −1 i=1 i=m+1 G = ||]F (δ)[|| ( wi − δ) , (4.3) i=m+1 G G wi − wi ) , ≤ ||]F (δ)[|| ( i=m+1 i=m

≤ ||]F (δ)[|| max wi max |ζj − ζj±1 | . i

j

In order to bound the dependence of the largest quadrature weight, wi , and the maximum possible distance |ζi − ζi±1 | in the interval (−1, 1), with the Gauss-Legendre rule, we least squares curve-fit (Gauss rules of 1 ≤ G ≤ 10) with the following, maxi≤g wi ≈ 1.93G−0.795 and maxξm ∈(−1,1) |ζi −ζi±1 | ≈ 2.6G−1.02 , with R2 = 0.99, where R2 = 1.0 indicates a perfect regression value of the curve fit. The product is maxi wi maxi |ζi − ζi±1 | ≈ 5.07G−1.82 .

(4.4)

Simple three dimensional estimates can be made by applying this procedure in all three directions on a reference element. It is clear that the error is directly related to the distance between Gauss point sampling locations. Table 4.1 illustrates the behavior. The “oversampling” only takes place in elements where there is a discontinuity (Figure 4.2). In elements with no material discontinuities the canonical rules should be used. This reduces computation time in forming the stiffness matrix. Consequently, we apply a “2/5” rule, i.e. a 2 × 2 × 2 Gauss rule if there is no discontinuity in the element, and a 5 × 5 × 5 rule if there is a discontinuity. We emphasize that this procedure is used simply to integrate elemental quantities with discontinuities accurately.

155


ORIGINAL=F(ζ ) 2 - δ ζ

CONTINUOUS=C( ζ ) ζ

M POINTS BEFORE

JUMP= [ F( ζ )] H( ζ - δ ) 2- δ

ζm ζm +1

ζ

TOTAL LENGTH=

a)

2

b) Figure 4.3. a) Decomposition of the discontinuous integrand on the reference element level. b) Gaussian quadrature points

Gauss rule 2×2×2 3×3×3 4×4×4 5×5×5

−1.82 )3

Bound = (5.07G8 0.37 0.04 0.008 0.002

Table 4.1. Integration error bound behavior for a unit jump. Domain: (−1, 1) × (−1, 1) × (−1, 1), for a 3-D step function discontinuity (8=2 × 2 × 2=volume of the reference element)

4.2.2 A Galerkin error monitor

In order to judge the quality of the numerical results, we need an indication of the Galerkin “subspatial” approximation error. As stated earlier, under standard assumptions the classical a-priori error estimate for the finite element method is ||u − uh ||E(Ω) ≤ C(u, p) def

hmin(s−1,p) = r , where p is the (complete) polynomial order of the finite element method used, s is the regularity of the exact solution and C is a global constant dependent on the exact solution and the polynomial approximation, but independent of h, the maximum element diameter. Related forms of this estimate hold locally, but with constants that are element dependent. The smoothness of the solution is no better than s = 1.5 (r = 0.5) for elements containing material discontinuities. However, the solution can be quite smooth, i.e. s > 1.5 (r > 0.5), for elements with no material discontinuities. For simplicity we employ the h-version finite element method, i.e. successive mesh subdivision using trilinear hexahedra (p = 1), and therefore r = 1 in elements with no discontinuities. We recall the well known property (PMPE), introduced earlier, ||u − uh ||2E(Ω) = 2(J (uh ) − J (u)). By solving the boundary value problem associated for a given sample for three successively finer meshes

156


MESH=h 1

~ 1 J(u ) h

~ J(u ) h

MESH=h 2

h

J(u~ 2) h

MESH=h 3

~ 3 J(u )

1/h h1

h

~ J(u )

~ 2 J(u )

h

~ 3 J(u )

Figure 4.4. Successively refined (embedded) meshes used to estimate the error

(Figure 4.4), h1 > h2 > h3 , with the following property: J (uh1 ) ≥ J (uh2 ) ≥ J (uh3 ) ≥ J (uh=0 ), we can set up the following system of equations for unknown constants a, b, and c, ahα1 + bhβ1 + c = J (uh1 ), ahα2 + bhβ2 + c = J (uh2 ),

(4.5)

ahα3 + bhβ3 + c = J (uh3 ), where the exponents α and β are selected to reflect the estimated regularity of the solution in the elements. In the case of using trilinear hexahedra we are led to the logical choice of exponents: α = 2r = 2, representing elements with no material discontinuities, and β = 2r = 1, representing elements with material discontinuities. It is clear that c ≈ J (uh=0 ), and we can thus define the following quantity to monitor the Galerkin error def

ξ =

2(J (uh3 ) − c) . E : IEΩ : E|Ω|

(4.6)

157


In summary, to monitor the discretization error, we apply the following algorithm:

STEP 1 : COARSE MESH h1 : uh1 ⇒ ||uh1 ||2E(Ω) STEP 2 : FINER MESH h2 = K1 × h1 : uh2 ⇒ ||uh2 ||2E(Ω)

(0 < K1 < 1)

STEP 3 : FINER MESH h3 = K2 × h2 : uh3 ⇒ ||uh3 ||2E(Ω)

(0 < K2 < 1)

def

STEP 4 : COMPUTE a, b, AND c ⇒ ξ =

(4.7)

2(J (uh3 ) − c) . E : IEΩ : E|Ω|

These operations are performed for each sample, since they each have a different random microstructure. 4.2.3 Remarks on local solution fields

As we have mentioned, for elements with no material discontinuities, the error in the energy 1 norm is bounded from above by Ch 2 . However, if we ask only to compare the average of the solutions over each element we have

(Ω − Ω ) : (σΩ − σ Ω ) dΩ ≤ h

h

Ω

Ω

12

(Ω − Ω ) : (Ω − Ω ) dΩ

h

h

≤C1 ha

×

12

(σΩ − σ Ω ) : (σΩ − σ Ω ) dΩ Ω

h

h

(4.8)

≤C2 hb

≤ C3 ha+b ,

where a, b ≥ 1. In short, the rates of convergence between the averaged quantities are much faster than between the unaveraged quantities. For example consider a one-dimensional example. Suppose the true solution contained in an element domain containing material du $ $ discontinuities is u = Ax ⇒ du dx = A, if x ≤ x and u = Bx + (A − B)x ⇒ dx = B if $ x ≥ x , as shown in Figure 4.5. Now consider a field produced by a trilinear finite element approximation The field will be a purely linear function, which we denote uh = Dx. The difference between the true and the approximate solution is  2  h  ∗ ∗ 1 h 1 h h    dx = h|(A x + B h − x − D)| ≤ Ch. dx − dx   0 h h  h 0 h 0

 %Ωe dx

%h Ωe dx

(4.9)

158


For example, take the special case when x∗ =

h 2

h h A+B 1 1 h h ( dx − dx)2 dx = h| − D|.

h

0

h

0

2

0

(4.10)

Therefore, we see that the rate of convergence is at least linear. However, this can be increased further, for example via more sophisticated post processing “extraction” techniques (see Szabo and Babuska [82]). Essentially the error in the element averaged stress and strain fields scale at least linearly with element size.

u

h u =Dx u=Bx+(A-B)x* u=Ax x

x* MATERIAL 1

MATERIAL 2

h Figure 4.5. An intuitive characterization of the rates of convergence for element averaged quantities

4.3 Iterative Krylov Solvers/microstructural “Correctors” For micro-macro mechanics problems, one important feature of iterative solvers, such as the CG method, is the fact that since they are based on successive updates of a starting guess solution vector. Consequently, they can be given a head start by a good solution guess, for example, provided by an inexpensive analytical approximation solution. For example, in a micro-macro mechanical context, the Conjugate Gradient (CG) method can be interpreted as making fine scale (high frequency) corrections with each iteration. Importantly, one should use inexpensive approximate solutions as initial guesses for the iterative solver. This further enhances the effectiveness of the CG searches for this class of problems. A well known fact, with regard to iterative solvers is the fact that they are quite adept at capturing local high-frequency responses. In other words, localized effects are resolved quickly during

159


CG iterations. However, Krylov methods typically may require many iterations to capture “long-wave” modes. For a discussion see Briggs [12]. In other words, a starting vector that captures a-priori the “long wave” components of the solution is advantageous. For example, if we assume that either of the the two classical boundary conditions that satisfy Hill’s conditions hold: (1) u|∂Ω = E · x ⇒ Ω = E and (2) t|∂Ω = L · n ⇒ σΩ = L. The classical assumptions of Voigt [87] (uniform strain) and Reuss [79] (uniform stress) produce the following: • u|∂Ω = E · x, which implies VOIGT: σ = IE : E, = E REUSS: σ = IE −1 Ω : E, = IE −1 : IE −1 −1 Ω :E • t|∂Ω = L · n, which implies VOIGT: = IE−1 : L, σ = IE : IE −1 Ω : L REUSS: σ = L, = IE −1 : L. The Reuss solution, while statically admissible, is highly oscillatory in the displacement field variable and kinematically inadmissible. The Voigt assumption, which is statically inadmissible and kinematically admissible, is not oscillatory in the displacement field variable. Therefore, from an iterative solver point of view, the Reuss field assumption is not of much value as a starting vector for a displacement based variational method. Furthermore, such guess cannot be projected onto the initial nodal values because it is double valued at material discontinuities (not kinematically admissible). However, it is important to note that if a complementary based numerical method were to be used, where the primary variable is the stress, the Reuss guess would be advantageous due to the smooth nature of stress field assumed. The Voigt (constant strain) assumption produces a low frequency initial guess. Therefore, from an iterative solver point of view, it is of value when using displacement based variational principles such as the standard FEM. Furthermore, such a guess can easily be projected onto initial nodal values, because it is single valued. The central point is that each Conjugate Gradient iteration can be viewed as a microstructural “corrector” to an initially statically inadmissible initial guess (Figure 4.6), and that u|Ω = E · x should be used for the displacement tests.

u=

ε.x

CORRECTIONS TO u=

ε.x

Figure 4.6. Krylov “corrections” to a smooth “Voigt” starting vector

160


4.4 Overall Testing Process: Numerical Examples As considered before, a typical example of a composite material combination is that of an aluminum matrix (77.9, 24.9 GPa) embedded with (stiffening) boron particles (230, 172 GPa). We chose AL/B as a material combination which exhibits significant enough mismatch in the mechanical properties to be representative of a wide range of cases. All tests were run on a single IBM RISC 6000 workstation. Comparable hardware is available in most academic and industrial work places, therefore such simulations are easily reproducible elsewhere for other parameter selections. 4.4.1 Successive sample enlargement

In a first set of tests, the number of particles contained in a sample were increased holding the volume fraction constant. During the tests, we repeatedly refined the mesh to obtain invariant macroscopic responses. A sample/particle size ratio was used as a microstrucdef , where tural control parameter. This was done by defining a subvolume size V = L×L×L N N is the number of particles in the entire sample and where L is the length of the (cubical L × L × L) sample. A generalized diameter was defined, d, which was the diameter of the smallest sphere that can enclose a single particle, of possibly non-spherical shape if desired. def The ratio between the generalized diameter and the subvolume was defined by ζ = d1 . V

3

For a variety of numerical tests, discussed momentarily, the typical mesh density to deliver invariant volumetrically averaged responses was 9 × 9 × 9 trilinear finite element hexahedra (approximately 2200-3000 degrees of freedom) per particle. We used ζ = 0.75, which resulted in a (fixed) volume fraction of approximately 22 %. The following particle per sample sequence (Figure 4.7) was used to study the dependence of the effective responses on the sample size: 2 (5184 DOF), 4 (10125 DOF), 8 (20577 DOF), 16 (41720 DOF), 32 (81000 DOF) and 64 (151959 DOF) particles. In order to get more reliable response data for each particle number set, the tests were performed five times (each time with a different particulate distribution) and the responses averaged. Throughout the tests, we considered a single combined boundary loading satisfying Hill’s condition, u|∂Ω = E·x, Eij = 0.001, i, j = 1, 2, 3. We tracked the strain energy, as well as κ∗ and µ∗ , as defined in Equation (2.1). Table 4.2 and Figure 4.8 depict the dependency of the responses with growth in particle number. Justified by the somewhat ad-hoc fact that for three successive enlargements of

2

4

8

16

32

64

SAME VOLUME FRACTION INCREASING NUMBER OF PARTICLES Figure 4.7. A series of test samples with increasingly more particles, but with the volume fraction fixed

161


98.5 KSTAR

EFFECTIVE BULK MODULI (GPa)

98

97.5

97

96.5

96

95.5 0

10

20

30 40 NUMBER OF PARTICLES

50

60

70

47 MUSTAR

EFFECTIVE SHEAR MODULI (GPa)

46

45

44

43

42

41 0

10

20


50

60

70

0.00145 ENERGY 0.00144 0.00143

SYSTEM ENERGY (GPa)

0.00142 0.00141 0.0014 0.00139 0.00138 0.00137 0.00136 0.00135 0

10

20


50

60

70

Figure 4.8. The values of the effective bulk and shear responses of the system for samples containing various numbers of particles. Five tests were performed per particle/sample size ratio and the results averaged

162


the number of particles, i.e. 16, 32 and 64 particle samples, the responses differed from one another, on average, by less than 1 %, we selected the 20-particle microstructures for further tests. We remark that we applied a “2/5” rule, i.e. a 2 × 2 × 2 Gauss rule if there is no discontinuity in the element, and a 5 × 5 × 5 rule if there is a discontinuity, which is consistent with the earlier derivation in the work (Table 4.1). The microstructure, as seen by this mesh density, is shown in Figure 4.9.

PART 2 4 8 16 32 64

d L

0.595 0.472 0.375 0.298 0.236 0.188

DOF 5184 10125 20577 41720 81000 151959

W (GPa) 0.001444 0.001408 0.001386 0.001375 0.001365 0.001358

κ∗ (GPa) 98.2 97.3 96.5 96.2 95.9 95.7

µ∗ (GPa) 46.7 44.3 43.2 42.5 41.6 41.4

Table 4.2. Results of successive sample “enlargements”. Five tests, each with a different random distribution, were performed at each sample/particulate size ratio level to obtain somewhat representative data

Figure 4.9. LEFT: a random microstructure consisting of 20 non-intersecting boron spheres, occupying approximately 22 % of the volume in an aluminum matrix, as seen by the algorithm with a 24 × 24 × 24 trilinear hexahedra mesh density for a total of 46875 degrees of freedom (approximately 9 × 9 × 9 hexahedra or 2344 degrees of freedom per element). A “2/5” rule, i.e. a 2 × 2 × 2 Gauss rule if there is no discontinuity in the element, and a 5 × 5 × 5 rule if there is a discontinuity, was used. RIGHT: a zoom on one particle

4.4.2 Multiple sample tests

For further tests, we simulated 100 different samples, each time with a different random distribution of 20 nonintersecting particles occupying 22 % (ζ = 0.75). Consistent with the previous test’s mesh densities per particle, we used a 24 × 24 × 24 mesh ( 9 × 9 × 9 trilinear hexahedra or 2344 DOF per particle, 46875 DOF per test sample). The numerical error

163


monitor, ξ, which measures microscopic differences (described in Box 4.7), was constructed for the each sample with a 24 × 24 × 24 mesh density by using a mesh sequence of embedded finite element meshes: 6 × 6 × 6, 12 × 12 × 12 and 24 × 24 × 24 (Figure 4.4). The 24 × 24 × 24 mesh was used in all final computations for each sample, with an accompanying ξ computed each time.

QUANTITY Ω ∇u : IE : ∇u dΩ(GPa) κ∗ (GPa) µ∗ (GPa) ξ CG − IT ERAT ION S

AV ERAGE 0.001373 96.171 42.350 0.0376 49.079

ST AN DARD DEV. 7.4 × 10−6 0.2025 0.4798 0.0045 1.5790

M AX − M IN 3.6 × 10−5 0.950 2.250 0.024 9

Table 4.3. Results of 100 material tests for randomly distributed particulate microstructures (20 spheres). We note that for a zero starting guess for the iterative solver, the average number of CG-iterations was 55, therefore using the initial Voigt (u = E · x) guess saved approximately 12% of the computational solution effort

The plots of the behavior of the various quantities of interest are shown in Figures 4.104.14. The averages, standard deviations and maximum/minimum of these quantities are tabulated in Table 4.3. For the 100 sample tests, with 20 particles per sample, the results for the effective responses were 91.37 = κ−1 −1 ˜ ∗ = 96.17 ≤ κΩ = 111.79, Ω ≤κ 30.76 =

µ−1 −1 Ω

≤

µ ˜∗

(4.11) = 42.35 ≤ µΩ = 57.68,

where κ ˜∗ and µ ˜∗ are the averaged effective responses from the 100 tests, and where the lower and upper bounds are, respectively, the classical Reuss [79] and Voigt [87] bounds. We also compared the computed results to the well-known Hashin-Shtrikman bounds ([29], [30]) which are, strictly speaking, only applicable to asymptotic cases of an infinite (sample length)/(particulate length) ratio and purely isotropic macroscopic responses. The “bounds” were as follows:

94.32 = κ1 +

v2 1 κ2 −κ1

+

3(1−v2 ) 3κ1 +4µ1

≤κ ˜∗ = 96.17 ≤

1 µ2 −µ1

+

v2 6(1−v2 )(κ1 +2µ1 ) 5µ1 (3κ1 +4µ1 )

shear modulus H/S lower bound

1 − v2 2 + 3κ23v +4µ2

1 κ1 −κ2

= 102.38,

bulk modulus H/S upper bound

bulk modulus H/S lower bound

35.43 = µ1 +

κ2 +

≤µ ˜∗ = 42.35 ≤ µ2 +

(1 − v2 ) 1 µ1 −µ2

+

6v2 (κ2 +2µ2 ) 5µ2 (3κ2 +4µ2 )

= 45.64,

shear modulus H/S upper bound

(4.12) where κ1 , µ1 and κ2 , µ2 are the bulk and shear moduli for the matrix and particle phases. Despite the fact that the bounds are technically inapplicable for finite sized samples, the

164


computed results did fall within them. The time to preprocess, solve and postprocess each 20 particle finite element test took no more than one minute on a single RISC 6000 workstation. Therefore, as before, 100 of such tests lasted approximately 1.5 hours.

25

SAMPLES

20

15

10

5

0.00136

0.00137

0.00138

0.00139

ENERGY (GPa)

Figure 4.10. 100 SAMPLES: The corresponding histogram for the variations in en1 ∇u : IE : ∇u dΩ(GPa) ergy, W = |Ω| Ω

Remarks: A primary question is: “What other information can we extract from the average of the responses of many samples? This is discussed next. 4.5 A Minimum Principle Interpretation Consider the following process for a sample of material with u|∂Ω = E · x: 1. Step 1: Take the sample, and cut it into N pieces, Ω = ∪N K=1 ΩK . The pieces do not have to be the same size or shape, although for illustration purposes it is convenient to take a uniform (regular) partitioning (Figure 4.15) 2. Step 2: Test each piece (solve the subdomain BVP) with the loading: u|∂ΩK = E · x. ˜ K is the solution to the BVP posed over subsample ΩK The function u 3. Step 3: One is guaranteed the following properties:

def ˜ ∗ ˜ ΩK = IE σ ΩK , K : ˜

˜ ∗ def ˜ ∗ |ΩK | IE = N K=1 IE K |Ω| ,

∗

∗

˜ − IE ∗ ) : E|Ω| ≤ E : (IE ˜ − IE −1 −1 ) : E|Ω|, ˜ ||2E(Ω) = E : (IE ||u − u Ω ∗

∗ ˜ IE −1 −1 Ω ≤ IE ≤ IE ≤ IEΩ , def

˜ 1 |Ω1 + u ˜ 2 |Ω2 ...˜ ˜ = u uN |ΩN . u

(4.13)

165


SAMPLES

20

15

10

5

95.75

96

96.25

96.5

EFFECTIVE BULK MODULUS (GPa)

Figure 4.11. 100 SAMPLES: The corresponding histogram for the variations in the effective bulk responses, κ∗

20

SAMPLES

15

10

5

0 41

41.5

42

42.5

43

EFFECTIVE SHEAR MODULUS (GPa)

Figure 4.12. 100 SAMPLES: The corresponding histogram for the variations in the effective shear responses, µ∗

166


SAMPLES

20

15

10

5

0.025

0.03

0.035

0.04

0.045

NORMALIZED NUMERICAL ERROR

Figure 4.13. 100 SAMPLES: The corresponding histogram for the variations in the numerical error monitor ξ

25

SAMPLES

20

15

10

5

0

46

48

50

52

CONJUGATE GRADIENT ITERATIONS

Figure 4.14. 100 SAMPLES: The corresponding histogram for the variations in Conjugate Gradient iterations needed for a solution

167


X2

Ω LARGE PROBLEM

Ωk

X1

PARTITIONED PROBLEM

STATISTICALLY REPRESENTATIVE SAMPLE

"SUBSAMPLES"

Figure 4.15. The idea of partitioning a sample into smaller samples or equivalently combining smaller samples into a larger sample

The same process can be done for traction test loading cases: t|∂ΩK = L·n. The effective material ordering, line three in Box 4.13, has been derived by Huet [40]. The second line of Box 4.13, and generalizations to nonuniform loading, were developed in Zohdi and Wriggers [93]. The proofs, which are constructive, and illustrative of the main ideas, are provided next. They result from a direct manipulation of classical energy minimization principles. 4.5.1 A proof based on partitioning results

Consider a large sample whose domain consists of the union of many smaller domains, ∪N K=1 ΩK = Ω, as depicted in Figure 4.15. The corresponding solution u, denoted the globally exact solution, is characterized by the following variational boundary value problem Find u ∈ H 1 (Ω), u|Γu = d, such that Ω

∇v : IE : ∇u dΩ =

Ω

f · v dΩ +

Γt

t · v dA

∀v ∈ H 1 (Ω), v|Γu = 0.

(4.14)

Specific choices for d, t and f will be given momentarily. Now consider the boundary of an individual subdomain ΩK , as ∂ΩK , 1 ≤ K ≤ N (N is the total number of subdomains) ˜ K ∈ H 1 (ΩK ), u ˜ K |∂ΩK ∩(Ω∪Γu ) = U ∈ H 1 (Ω), such that Find u ΩK

˜ K dΩ = ∇v K : σ

ΩK

f · v K dΩ +

∂ΩK ∩Γt

t · v K dA

(4.15)

∀v K ∈ H 1 (ΩK ), v K |∂ΩK ∩(Ω∪Γu ) = 0. A specific choice for U will also be given momentarily. The individual subdomain solutions form an approximate solution to the globally exact problem, u. This approximate solution is constructed by a direct assembly process def

˜ = U + (˜ u u1 − U )|Ω1 + (˜ u2 − U )|Ω2 + ... + (˜ uN − U )|ΩN .

(4.16)

168


The approximate displacement field is continuous, however, the approximate traction field is usually discontinuous. Clearly, if there are no jumps in the tractions, the solution is exact. 4.5.2 Relation to the material tests

We have for any kinematically admissible function w, using the (bilinear operator) notation, def 0 ≤ ||u − w||2E(Ω) = B(u − w, u − w) = 2J (w) − 2J (u) ⇒ J (u) ≤ J (w), where J (w) = 1 2 B(w, w) − F(w). In other words, the true solution has a minimum potential (PMPE). By applying the PMPE, for the test loading u|∂Ω = E · x = d, f = 0, with the specific (kinematically admissible) choice w = U , we obtain ||u − U ||2E(Ω) = 2(J (U ) − J (u)) ⇒ def

J (u) = J (U ) − 12 ||u − U ||2E(Ω) . The critical observation is that when choosing w = U = ˜ , as defined in Box 4.16, is also kinematically admissible. Therefore, E ·x in Box 4.15, then u ˜ ||2E(Ω) = 2(J (˜ we have by direct expansion ||u − u u) − J (u)) = 2(J (˜ u) − J (U )) + ||u − ˜ K is a solution to a subdomain boundary value problem posed over ΩK , U ||2E(Ω) . Since u it minimizes the corresponding subdomain potential energy function JK (·). Therefore, uK ) = J (˜ u). Consequently, J (U ) = SK=1 JK (U ) ≥ SK=1 JK (˜ ˜ ||2E(Ω) = 2(J (˜ ||u − u u) − J (U )) +||u − U ||2E(Ω) .

(4.17)

negative

By direct expansion we have ||u−U||2E ( Ω) = E : (IE−IE ∗ ) : E|Ω| ≤ E : (IE−IE −1 −1 ) : E|Ω|. We have by definition, 2J (u) = E : IE ∗ : E|Ω|, 2J (U ) = E : IE : E|Ω|, and ˜ ∗ : E|ΩK |, which implies 2J (˜ ˜ ∗ : E|Ω|. By direct substitution u) = E : IE u) = E : IE 2JK (˜ K this completes the first of the assertions in Box 4.13. We can directly interpret the testing of many smaller samples as simply the partitioning of a very large one which we cannot easily test (Figure 4.15). Therefore, for displacement tests, the averaged effective responses generated will always bound the response of the very large sample from above. Therefore, the average of the 100 sample tests provide us with tighter upper bounds on the response of a very large sample. Therefore, by the previously proven results in Box 4.13, the Reuss-Voigt bounds in Box 4.11 are tightened by the following factors 96.17 − 91.37 = 0.2747 111.79 − 91.37

and

42.35 − 30.76 = 0.4305. 57.68 − 30.76

(4.18)

4.6 Partitioning and Traction Test Cases Suppose we repeat the partitioning process for an applied internal traction set of tests. As will be shown traction tests will bound the response of the very large sample from below. Employing the Principle of Complementary Potential Energy, we obtain for the test loading t|∂ΩK = L · n, with γ = L (a statically admissible trial field), ˜ − K(σ)) = 2(K(σ) ˜ − K(γ)) +||σ − γ||2E −1 (Ω) , ˜ 2E −1 (Ω) = 2(K(σ) ||σ − σ||

(4.19)

negative

˜ being the stress field produced by where K(σ) = 12 Ω σ : IE −1 : σ dΩ − Γu σ · n · d dA, σ solving each subsample problem with the uniform stress boundary conditions and forming def ˜ = σ ˜ 1 |Ω1 + σ ˜ 2 |Ω2 ...σ ˜ N |ΩN . By direct expansion we have ||σ − γ||2E −1 (Ω) = L : (IE −1 − σ def


169

IE −1∗ ) : L|Ω| ≤ L : (IE −1 − IE−1 ) : L|Ω|. The complementary forms collapse to ˜ −1∗ : L|ΩK | ⇒ ˜ = L : IE 2K(σ) = L : IE −1∗ : L|Ω|, 2K(γ) = L : IE −1 : L|Ω| and 2KK (σ) K ˜ −1∗ : L|Ω|, where KK (·) is the subdomain complementary potential energy ˜ = L : IE 2K(σ) function for subdomain K. By direct substitution this yields ˜ −1∗ : σ ˜ −1∗ def ˜ −1∗ |ΩK | , ˜ ΩK , ˜ΩK = IE IE = N K=1 IE K K |Ω| ˜ −1∗ − IE ∗−1 ) : L|Ω|, ˜ 2 −1 = L : (IE ||σ − σ|| E

(Ω)

˜ −1∗ − IE−1 ) : L|Ω|, ≤ L : (IE Ω −1∗ −1 −1 −1∗ ˜ ≥ IE ≥ IEΩ . IE Ω ≥ IE

(4.20)

As for the primal (displacement) case, the second line of Box 4.20, and generalizations to nonuniform loading, were developed in Zohdi and Wriggers [93]. The fourth line has been derived in Huet [40] by other means. Consequences/difficulties with traction tests If the sample were an RVE, we have IE −1∗ = IE ∗−1 , then the preceding analysis yields the following two sided ordering of approximate effective material responses, ˜ −1∗ )−1 ≤ IE ∗ ≤ IE ˜ ∗ ≤ IEΩ . IE −1 −1 Ω ≤ (IE

(4.21)

We emphasize that IE −1∗ = IE ∗−1 is an assumption, which may not be true for a finite sized sample. Therefore, in theory, under the RVE assumption, traction tests form lower bounds on the effective responses. However, traction tests pose great difficulties, which are as follows: 1. Numerically pure traction boundary data cause rigid motions (singular FEM stiffness matrices), ˜ 2. The FEM is a method based upon generating kinematically admissible solutions σ. The traction tests result is based upon the assumption that statically admissible trial field are generated. Statically admissible fields cannot be achieved by a standard FEM approach, 3. True laboratory tests specifying the force on a sample are far more difficult than specifying the displacements. Isolating the subsampling error/numerical error orthogonality According to the results in Box 4.13, we have the following normalized estimate: ˜ ∗ − IE ∗ ) : E|Ω| ˜ ∗ − IE −1 −1 ) : E|Ω| E : (IE E : (IE Ω = ≤ , E : IEΩ : E|Ω| E : IEΩ : E|Ω| E : IEΩ : E|Ω| ˜ ||2E(Ω) ||u − u

(4.22)

where, since we have used displacement controlled tests, we have used the Voigt material’s energy, which corresponds to assuming a constant strain throughout the material equal to E, to normalize the results.

170


It is important to note that the numerical error, which is implicitly included in the estimate in Equation 4.22, is orthogonal to the “partitioning” or “subsampling” error. In other words, while the numerical error is known, it can be directly filtered out of the sampling error estimates. To see this, consider that the boundary value formulations in ˜ ) dΩ = 0, Boxes 4.14 and 4.15 directly imply, for any type of loading Ω ∇v : IE : ∇(u − u 1 ˜ h ) : IE : ∀v ∈ H (Ω), v|Ω∩∂ΩK = 0, K = 1, 2, 3...N . This directly implies that Ω ∇(˜ u−u ˜ ) dΩ = 0. Therefore, ∇(u − u ˜ h ||2E(Ω) ||u − u

= = −

˜+u ˜−u ˜ h ||2E(Ω) ||u − u ˜ ) : IE : ∇(u − u ˜ ) dΩ ∇(u − u Ω ˜ h ) : IE : ∇(u − u ˜ ) dΩ + ˜ h ) : IE : ∇(˜ ˜ h ) dΩ 2 ∇(˜ u−u ∇(˜ u−u u−u Ω Ω

=0 by the above

=

||u −

˜ ||2E(Ω) u

˜ h ||2E(Ω) . + ||˜ u−u

(4.23)

def

Under the special case that U = E · x, we have the following orthogonal decomposition of sampling and numerical error:

˜ E : (IE

∗,h

∗

˜ − IE ∗ ) : E|Ω| + − IE ∗ ) : E|Ω| = E : (IE h

˜ ||2 ||u−u E(Ω)

˜ ||2 ||u−u E(Ω)

N K=1

˜ K ||2E(Ω) , ||˜ uhK − u

˜ ∗,h ˜ ∗ E:(IE −IE ):E|Ω|

(4.24)

def ˜ ∗,h ˜ ∗,h |ΩK | ˜ ∗,h def ˜ h ΩK = IE h ΩK and IE = N where σ K : ˜ K=1 IE K |Ω| . Therefore, the expression in Box 4.24 indicates that the estimates made on the effective response is an orthogonal sum of both the numerical and sampling error. Therefore, one can isolate either the numerical or sampling error, if the other is known, or estimated, by simply subtracting it from the the total expression (the lefthand side of Box 4.24). Therefore, for example, approximating the numerical error by the error monitor, ξ, we have

˜ ∗ − IE −1 −1 ) : E|Ω| E : (IE Ω −ξ SU BSAM P LIN G ERROR = E : IEΩ : E|Ω|

(4.25)

= 0.1073 − 0.0376 = 0.0697. This quantity represents the isolated error in the internal fields induced by the partitioning or subsampling only. More general theoretical relationships between sampling and numerical error, which are beyond the scope of this monograph, can be found in Zohdi and Wriggers [93]. 4.7 Dependency on Volume Fraction We repeated the 100 sample tests procedure for samples containing significantly higher volume fraction, approximately 32 %. Unlike fiber-reinforced composites, which can contain well over 50 % volume fraction of fibers, particulate composites usually contain no more

171


than approximately 25 % volume fraction of particulates. Particulate volume fractions over 15 % are already high in three dimensions. For the relevant examples, 22 % particulate volume fraction in three dimensions corresponds roughly to 44 % in two dimensions, while 32 % in three dimension corresponds roughly to 57 % in two dimensions. Cross-sections of these volume fractions are depicted in Figure 4.16. The results for the 32 % case are tabulated in Table 4.4. The averages for the effective properties of the samples were 99.04 = κ−1 −1 ˜ ∗ = 106.83 ≤ κΩ = 126.99, Ω ≤κ (4.26)

34.39 = µ−1 −1 ˜∗ = 52.63 ≤ µΩ = 72.38, Ω ≤µ

where κ ˜ ∗ and µ ˜∗ are the averaged quantities from the 100 tests, and where the lower and upper bounds are, respectively, the classical Reuss [79] and Voigt [87] bounds. The ReussVoigt bounds in Box 4.26 are tightened by the following factors 106.831 − 99.043 = 0.2785 126.999 − 99.043

and

52.635 − 34.395 = 0.4801. 72.385 − 34.395

(4.27)

The error estimate in this case is ∗

SU BSAM P LIN G ERROR ≈

˜ − IE −1 −1 ) : E|Ω| E : (IE Ω −ξ E : IEΩ : E|Ω|

(4.28)

= 0.1436 − 0.0523 = 0.0913, which is approximately 2 % higher error estimate than the case of 22 % particulate volume fraction. Intuitively one would expect such a result, since the particles interact more at such a high volume fraction (Figure 4.16). In other words, they “feel” the presence of the other particles more at higher volume fractions than at lower volume fractions. As before, we also compared the computed results to the Hashin-Shtrikman bounds [29] and [30] which are, strictly speaking, only applicable to asymptotic cases of an infinite (sample length)/(particulate length) ratio, and purely isotropic macroscopic responses

103.378 = κ1 +

v2 1 κ2 −κ1

+

3(1−v2 ) 3κ1 +4µ1

≤κ ˜ ∗ = 106.831 ≤

1 µ2 −µ1

v2 6(1−v2 )(κ1 +2µ1 ) 5µ1 (3κ1 +4µ1 )

+

shear modulus H/S lower bound

1 − v2 2 + 3κ23v +4µ2

1 κ1 −κ2

= 114.671,

bulk modulus H/S upper bound

bulk modulus H/S lower bound

41.624 = µ1 +

κ2 +

≤µ ˜∗ = 52.635 ≤ µ2 +

(1 − v2 ) 1 µ1 −µ2

6v2 (κ2 +2µ2 ) + 5µ (3κ +4µ2 ) 2 2

= 56.437.

shear modulus H/S upper bound

(4.29)

172


22 % VOLUME FRACTION SIZE 32 % VOLUME FRACTION SIZE Figure 4.16. A comparison of relative sizes of particles occupying approximately 22 % (ζ = 0.75) and 32 % (ζ = 0.85) volume fraction

SUBDOMAIN QUANTITY W (GPa) κ∗ (GPa) µ∗ (GPa) ξ CG − IT ERAT ION S

AV ER. 0.001593 106.831 52.635 0.0523 48.950

ST AN. DEV. 9.11 × 10−6 0.2505 0.5876 0.0044 1.716

M AX − M IN 5 × 10−5 1.350 3.095 0.0255 8

Table 4.4. Results of 100 material tests for randomly distributed particulate microstructures (20 spheres, ζ=0.85, or approximately 32 % volume fraction)

4.8 Increasing the Number of Samples Suppose we increase the number of samples tested. For purposes illustration purposes, let us test 512 samples of material. The number 512 is not accidental, since it is a common number of independent processors in modern parallel processing machines. Table 4.5 illustrates that the averaged results are virtually identical to the 100 sample tests for all the quantities. Testing more and more samples will not help obtain better average results. However, despite practically the same average values, one can observe from the Figures 4.17-4.20 that the 512 sample tests have a more Gaussian distribution, relative to the 100 sample tests, for the responses. However, for even more accurate average responses, we must test larger samples of material. This is explored further in the next section. 4.9 Increasing Sample Size Beyond a certain threshold, it is simply impossible to obtain any more information by testing samples of a certain size. As we have shown the reason for the inability for obtaining further information is the fact that the testing conditions are uniform on the “subsamples”. This idealization is valid only for an infinitely large sample. However, suppose we increase the number of particles per subsample even further, from 20 to say 40 then 60, each time performing the 100 tests procedure. With this information one could then possibly extrapolate to a (giant) sample limit. The results for the 40 and 60 particle cases are shown in

173


SUBDOMAIN QUANTITY def

W =

1 |Ω|

AV ERAGE

ST AN DARD DEV.

M AX − M IN

0.001373

7.2 × 10−6

7.5 × 10−5

96.169 42.353 48.8984

0.1967 0.4647 1.7737

1.203 3.207 11

∇u : IE : ∇u dΩ(GPa) Ω

κ∗ (GPa) µ∗ (GPa) CG − IT ERAT ION S

Table 4.5. Results of 512 material subdomains, each containing 20 randomly distributed spheres for a total of 10240

140 130 120 110

SAMPLES

100 90 80 70 60 50 40 30 20 10 0.00135

0.00136

0.00137

0.00138

0.00139

ENERGY (GPa)

Figure 4.17. 512 SAMPLES: The corresponding histogram for the variations in endef 1 ∇u : IE : ∇u dΩ(GPa) ergy, W = |Ω| Ω

Table 4.6 for 22 % boron volume fraction. Using these results, along with the 20 particle per sample tests, we have the following curve fits

W = 0.0013205 + 0.0001895 Ld , κ∗ = 94.527 + 5.909 Ld , µ∗ = 39.345 + 10.775 Ld ,

R2 = 0.9997,

R2 = 0.986, R2 = 0.986,

(4.30)

174


120 110 100

SAMPLES

90 80 70 60 50 40 30 20 10 96

96.5

EFFECTIVE BULK MODULUS (GPa)

Figure 4.18. 512 SAMPLES: The corresponding histogram for the variations in the effective bulk responses, κ∗

160 140 120

SAMPLES

100 80 60 40 20

41

42

43

EFFECTIVE SHEAR MODULUS (GPa)

Figure 4.19. 512 SAMPLES: The corresponding histogram for the variations in the effective shear responses, µ∗

175


110 100 90

SAMPLES

80 70 60 50 40 30 20 10 44

46

48

50

52

CONJUGATE GRADIENT ITERATIONS

Figure 4.20. 512 SAMPLES: The corresponding histogram for the variations in Conjugate Gradient iterations needed for a solution

where L is the sample size, d is the diameter of the particles. Thus as Ld → 0, we obtain estimates of W = 0.0013205 GPa, κ∗ = 94.527 GPa and µ∗ = 39.345 GPa as the asymptotic energy, effective bulk modulus, and effective shear modulus, respectively. Indeed, judging from the degree of accuracy of the curve-fit, (R2 = 1.0 is perfect) the curve fit appears to be perfect for W . For κ∗ and µ∗ , the slightly less accurate reliability (regression values of R2 = 0.98) is attributed to the fact that absolute perfect isotropy is impossible to achieve with finite sized samples. In other words the extrapolations using various samples exhibit slight isotropic inconsistencies. However, the energy W , has no built-in assumptions whatsoever, thus leading to the nearly perfect curve-fit. The monotonicity of the testing curves is to be expected, and is explained further in the next chapter.

PART 40 40 40 40 60 60 60 60

d L

0.2193 0.2193 0.2193 0.2193 0.1916 0.1916 0.1916 0.1916

DOF 98304 98304 98304 98304 139968 139968 139968 139968

QUANTITY W (GPa) κ∗ (GPa) µ∗ (GPa) CG − IT ER W (GPa) κ∗ (GPa) µ∗ (GPa) CG − IT ER

AV ER. 0.0013617 95.7900 41.6407 60.3200 0.0013586 95.6820 41.4621 66.500

ST AN. DEV. 5.0998 × 10−6 0.1413 0.3245 1.6302 4.3963 × 10−6 0.1197 0.2801 1.9974

M AX − M IN 2.35 × 10−5 0.6600 1.590 7 2.25 × 10−5 0.6214 1.503 10

Table 4.6. Results of material tests for randomly distributed particulate microstructures for 100 (ζ = 0.75, approximately 22 %) samples of 40 and 60 particles per sample

176


4.10 An Ergodic Interpretation of the Results It is possible to interpret the results in ergodic manner, which we now briefly discuss. In statistical mechanics a frequently invoked assumption is that of statistical ergodicity. Importantly, in the context of material testing, Hill’s condition can be interpreted as such an assumption. To see this, consider several bodies with the same external geometry composed of the same material volume fractions but with different random microstructure. Let us focus on the same spot on each body (indexed by i = 1, 2, ...N ) and compute 1 N (i) def (i) = IEI ∗(i) : ( N1 N i=1 σ i=1 ), where N is the number of macroscopic structures. N ∗ The tensorial quantity, IEI , is called an ensemble average (Figure 4.21). Now we perform another experiment. We have a single body, and focus on one spot, but repeatedly compute and enlarge the sample size (indexed by j = 1, 2, ...N ), σ(j) Ωj = IE ∗(j) : (j) Ωj . This is our familiar (volumetric averaged) effective property. A classical ergodicity assumption is that as i, j → ∞ that IEI ∗ = IE ∗ . Therefore, volumetric averaging and ensemble averaging yield the same result. In other words, an infinitely large sample’s volumetric average must equal the ensemble average of infinitely many finite samples at a point, from different bodies. ˜ ˜) and If we split the stress and strain fields into a purely fluctuating (zero mean) (σ, average parts (σΩ , Ω ) we see by direct expansion:

˜ ) : (Ω + ˜ ˜ Ω : Ω + σΩ : ˜Ω +σ ˜ : ˜Ω (σΩ + σ )Ω = σΩ : Ω + σ

=0

=0

˜ : ˜Ω . = σΩ : Ω + σ

(4.31)

˜ : ˜Ω = 0. In other words, the product of two purely The ergodicity assumption is that σ fluctuating random fields is again purely fluctuating (Figure 4.22). The implication is that ˜ : ˜ is purely fluctuating. This is exactly the as the sample becomes infinitely large that σ implication of the Hill condition. This motivates the use of exterior uniform loading for an infinitely large sample. Clearly, the partitioning, and the error estimates, can be interpreted testing an ergodic assumption. However, the loadings are still not “rough” enough, since it is clear that uniform loading is an idealization and thus will only be present within the microstructure of a finite sized engineering (macro)structure to an error. For more on ergodic hypotheses see the classical work of Kröner [52].

5 MATHEMATICAL GENERALIZATIONS OF MULTIPLE TESTING One can consider the multiple sample testing introduced in the last chapter as a type of domain decomposition. In order to see this, we consider a “partitioning” approach whereby a global domain, under arbitrary loading, is divided into nonoverlapping subdomains. On the interior subdomain partitions an approximate globally kinematically admissible solution is projected. This allows the subdomains to be mutually decoupled, and therefore separately solvable. The subdomain boundary value problems are solved with the exact microstructural representation contained within their respective boundaries, but with approximate displacement boundary data. The resulting microstructural solution is the assembly of the subdomain solutions, each restricted to its corresponding subdomain. The approximate solution is far more inexpensive to compute than the direct problem. The work follows from results found in Zohdi et al. [95].

177


ENSEMBLE AVERAGING

MORE BODIES VOLUMETRIC AVERAGING

LARGER SAMPLES

Figure 4.21. Ensemble and volumetric averaging processes.

σ

SAMPLE LENGTH ε

SAMPLE LENGTH

σ:ε

SAMPLE LENGTH Figure 4.22. An illustration of an ergodicity assumption. Essentially it states that the product of two random purely fluctuating fields is also purely fluctuating and random

178


5.1 Boundary Value Problem Formulations The globally exact solution, u, is characterized by the following virtual work formulation: Find u ∈ H 1 (Ω), u|Γu = d, such that ∇v : σ dΩ = f · v dΩ + t · v dA Ω

Ω

Γt

∀v ∈ H 1 (Ω), v|Γu = 0,

(5.1)

where σ is the Cauchy stress. In the infinitesimal strain, linearly elastic, case σ = IE : ∇u. In order to construct approximate solutions, next we consider the subdomain boundary value problems, which have the exact microstructural representation contained within the domain, but approximate displacement boundary data on the interior subdomain boundaries. To construct the approximate microstructural solutions we first partition the domain, Ω, into N nonintersecting open subdomains ∪N K=1 ΩK = Ω. We define the boundary of an individual subdomain ΩK , as ∂ΩK . When employing the applied internal displacement approach, a kinematically admissible function, U ∈ H 1 (Ω) and U |Γu = d, is projected onto the internal boundaries of the subdomain partitions. Any subdomain boundaries coinciding with the exterior surface retain their original boundary conditions (Figure 5.1). Accordingly, we have the following virtual work formulation, for each subdomain, 1 ≤ K ≤ N : ˜ K ∈ H 1 (ΩK ), u ˜ K |∂ΩK ∩(Ω∪Γu ) = U ∈ H 1 (Ω), such that Find u ˜ K dΩ = ∇v K : σ f · v K dΩ + t · v K dA ΩK

∂ΩK ∩Γt

ΩK

∀v K ∈ H 1 (ΩK ), v K |∂ΩK ∩(Ω∪Γu ) = 0.

PROJECT A KINEMATICALLY ADMISSIBLE FUNCTION

Figure 5.1. Construction of the approximate solution

(5.2)

179


The constitutive law and microstructure are identical to that of the globally exact problem. ˜ = IE : ∇˜ In the infinitesimal strain, linearly elastic, case σ uK . The individual subdomain ˜ K , are zero outside of the corresponding subdomain ΩK . In this case the solutions, u approximate solution is constructed by a direct assembly process def

˜ = U + (˜ u1 − U )|Ω1 + (˜ u2 − U )|Ω2 + ... + (˜ uN − U )|ΩN . u

(5.3)

The approximate displacement field is in H 1 (Ω), however, the approximate traction field is possibly discontinuous. Logical choices of U will be given later in the work. It should be clear that if U = u on the internal partition boundaries, then the approximate solution is exact. Remark I: In theory, an approach of projecting the tractions could be developed. However, such an approach has a number of difficulties, in particular the pure traction subdomain problems must have equilibrated tractions to be well posed. To impose this on the numerical implementation level is not a trivial task. Primarily for this reason the projected displacement approach is preferable. Remark I: For example, if one considers a cubical domain, whose response is simulated using n numerical unknowns, for example employing a finite element discretization, approximately O(nγ ) floating point operations are needed to solve the system. The value of γ, typically 2 ≤ γ ≤ 3, depends on the type of algebraic system solver used. If the cube was divided into N equal subdomains (subcubes), n γthe number of floating point operations needed would N , and thus DIRECT COST S/DECOM P OSED be approximately on the order of N γ−1 . For example, if one had 1000 subdomains, the “broken” solution costs COST S ≈ N between 1000 and 1000000 times less to compute than the globally exact solution. The advantages are not limited to the possible reduction of operation counts, since the subdomain problems can be solved separately, and trivially in parallel. If parallel processing is used for the decomposed problem, which by construction has a perfect speedup, then the ratio of costs becomes P N γ−1 , where P denotes the number of processors. 5.2 Error in the “Broken” Problems Since we employ energy type variational principles to generate approximate solutions, we use energy type measures for the error. Directly by use of the divergence theorem one has ˜ ) : (σ − σ) ˜ dΩ = ∇(u − u Ω

N

˜ K ) : (σ − σ ˜ K ) dΩK , ∇(u − u

K=1 ΩK N

˜ K ) · (σ − σ ˜ K )) dΩK , ∇ · ((u − u

=

−

K=1 ΩK N K=1

=

˜ K ) · ∇ · (σ − σ ˜ K ) dΩK , (u − u

N

K=1

=

ΩK

∂ΩK ∩Ω

NI

I=1

ΓI

=0

˜ K · nK ) · (˜ (σ uK − u) dAI ,

not continuous

error

(traction jump) · (displacement error) dAI ,

(5.4)

180


DISCONTINUOUS TRACTIONS

Figure 5.2. The consequences of the projection approach

where ΓI is an interior subdomain interface, and I = 1, 2, ...NI =number of interior subdomain interfaces. For the applied internal displacement case, the tractions may suffer discontinuities at the interior subdomain boundaries (see Figure 5.2). In the case of infinitesimal strain, linear elasticity, we have

||u −

˜ ||2E(Ω) u

=

N K=1

def

˜ ||2E(Ω) = ||u − u

∂ΩK ∩Ω

(IE : ∇˜ uK · nK ) · (˜ uK − u) dAI ,

not continuous

error

(5.5)

˜ ) : IE : ∇(u − u ˜ ) dΩ. ∇(u − u Ω

It is convenient to cast the error in terms of the potential energy for the case of lin def ear elasticity, J (w) = 12 Ω ∇w : IE : ∇w dΩ − Ω f · w dΩ − Γt t · w dA, where w is any kinematically admissible function. The well known relationship for any kinematically admissible function w is ||u − w||2E(Ω) = 2(J (w) − J (u)) or J (u) ≤ J (w),

(5.6)

which is the Principle of Minimum Potential Energy (PMPE). In other words, the true ˜ is kinematically admissible, we solution possesses a minimum potential. Therefore, since u ˜ ||2E(Ω) = 2(J (˜ u) − J (u)). The critical observation is that if we immediately have ||u − u can bound J (u) from below, then we can bound the error from above. In other words, what we seek is J − ≤ J (u). To bound J (u) in terms of easily accessible quantities, we first calculate a computationally inexpensive, kinematically admissible, regularized test solution. The regularized test solution, uR , is characterized by a virtual work formulation: Find a uR ∈ H 1 (Ω), uR |Γu = d, such that ∇v : σR dΩ = f · v dΩ + t · v dA Ω

Ω

Γt

∀v ∈ H 1 (Ω), v|Γu = 0.

(5.7)

181


Here the constitutive law for σR should be taken to be as simple as a possible, i.e. a constant (regularized) fourth order symmetric positive definite linear elasticity tensor IR, defined via σR = IR : ∇uR . In general, the regularized solution (σ R , uR ) does not coincide with the field associated with the decoupling solution U . If we choose w = uR in Equation 5.6, which is an admissible function, we obtain ||u−uR ||2E(Ω) = 2(J (uR )−J (u)), which implies J (u) = J (uR )− 12 ||u−uR ||2E(Ω) . Our objective is to form an upper bound on ||u−uR ||E(Ω) in terms of IR, IE and ∇uR , to obtain a lower bound on J (u). For the bound to be useful, it should contain no unknown constants, and should be solely in terms of the regularized solution and the material data. In other words we seek ||u − uR ||E(Ω) ≤ H+,(u), which will lead to J − = J (uR ) − 12 (H+,(u) )2 ≤ J (u), thus supplying an upper bound for the quantities in Box 5.5. By the proof and results in Box 3.6, one immediately has def

||u −

uR ||2E(Ω)

≤

2

((IE − IR) : ∇uR ) : IE −1 : ((IE − IR) : ∇uR ) dΩ = H+,(u) . def

(5.8)

Ω

To derive this result we used the fact that a symmetric positive definite matrix has a 1 1 unique positive square root, IE = IE 2 : IE 2 and the Cauchy-Schwarz inequality. Note that the H+,(u) -bound depends only on IR, IE and ∇uR . Therefore, using the Box 5.8 we have def J − = J (uR ) − 12 (H+,(u) )2 ≤ J (u), and thus ˜ ||2E(Ω) ||u − u ≤ 2(J (˜ u) − J − ) .

U P P ER BOU N D ˜ )−J (u)) ERROR=2(J (u

(5.9)

Remark I: Bounds of such type have been studied in Zohdi et. al [91], Oden and Zohdi [70] and Zohdi and Wriggers [92], however, they required that the decoupling solution U be identical to uR . Clearly this is unnecessary. In the preceding bound, the tensor IR did not have to be constant, and this has been exploited to generate hierarchical micromechanical models in Zohdi et. al [91]. More general theoretical aspects of such bounds, including their complementary (dual) formulation counterparts, can be found in Zohdi and Wriggers [92] and Zohdi [100]. Under certain conditions, such bounds coincide with results found in Huet [37], [38]. Remark II: It is important to realize that the only approximation in this entire estimate +,(u) def is an extremely good indicator is ||u − uR ||E(Ω) ≤ H+,(u) . Thus the ratio > = ||u−HuR || E(Ω) of the quality of the estimate. Numerical and theoretical studies of the proximity of > to unity has been studied in Zohdi et. al [91], Oden and Zohdi [70], Zohdi and Wriggers [92][97] and Zohdi [100]. It is clear that the tensorial parameter IR is free to choose, provided it is symmetric and positive definite. Thus, ideally, one would want to minimize 2 = 0, ∂ (,−1) > 0 ⇒ IR$ . Unfortunately, minimization of > − 1, in other words ∂(,−1) 2 ∂ IR ∂ IR > − 1 is usually impossible, since ||u − uR ||E(Ω) is unknown. However, another avenue to optimize the bound is possible. Formally, the process to obtain the optimal choice, IR$ , − $ ∂2J − is to maximize the lower bound J − , in other words ∂∂J IR = 0, ∂ IR2 < 0 ⇒ IR . For example, if we consider a one-dimensional displacement controlled structure d 0 < x < L : dx (E(x) du )=0 dx

u(0) = ∆0

u(L) = ∆L ,

(5.10)

182


it is easy to show that

H ||u − uR ||E(Ω) +,(u)

>=

1 1 ( − 1 )2 E dx = 10 1R 1E ( − E β)2E dx 0 R

where the estimate is exact (> = 1) for R = E1 −1 L =

def

def

β =

1 R L , E1 L

1 L 1 L 0 E

!−1

dx

(5.11)

. This value of R

J −.

We remark that for traction controlled loading u(0) = 0, E du also maximizes dx (L) = G, for a given G, the error estimate is exact for any R. For three dimensional problems, the optimal tensor IR$ , which is problem dependent, is easily computable, since it requires no microscale computations. 5.3 The Connection to Material Testing As an example, we apply the decomposition approach to a sample of material containing a significant amount of microstructure (Figure 5.3). We simulated a rectangular block composed of an Aluminum matrix, containing 10240 Boron particles. The partitioning consisted of 512 subblocks (subdomains), each containing a random distribution of 20 nonintersecting particles occupying approximately 22 % of volume. We considered external boundary loading of the form Eij = 0.001, i, j = 1, 2, 3: 









u1 |∂Ω E11 E12 E13 x1  u2 |∂Ω  =  E12 E22 E23   x2  , u3 |∂Ω E31 E32 E33 x3

E

(5.12)

which is relatively standard in micromechanical analyses. The usual motivation for such a loading is that the sample is considered so large that it experiences nearly uniform strain (linear displacements) on its boundary. For this loading, a straightforward kinematically admissible choice for U (in Box 5.2) is U = E ·x. Before we ran the full tests, we repeatedly refined the finite element mesh for several subdomains (subblocks), each containing the same number of spheres, but with different random distributions, until we obtained invariant total strain energy responses. The result was that beyond typically a 24 × 24 × 24 mesh (46875 dof per test) per subdomain, roughly 2344 dof per particle, delivered macroscopically (volumetrically averaged) mesh invariant responses. The resulting number of total degrees of freedom (n), if the entire domain was directly simulated with this mesh density, would be n=21,567,171. Important: Under these special loading conditions, this decomposition process can be considered as equivalent to the statistical testing of 512 samples of heterogeneous materials in the previous chapter, since the loading on each subdomain will be uniform. It is important to realize that since all subdomains are the same size, that the projection of U = E · x on the interfacial boundaries would produce the same strain energy in each subdomain if the material had been uniform.‡ Therefore, it is logical to (statistically) compare the solution behavior between the subdomains. Throughout the tests a standard preconditioned Conjugate Gradient (CG) solver was used. The primary premise for its (wide) use is that a solution, to tolerable accuracy, can be achieved in much less than O(n3 ) operations, which is required with most Gaussian-type techniques. For overviews of a variety of such methods see Axelsson [8]. The CG method, is guaranteed to converge in n iterations, provided the algebra is performed exactly. At each iteration the ‡

Clearly, it is a homogeneous deformation in this case.

183


512 SUBDOMAINS

20 PARTICLES PER SUBDOMAIN

MESH DENSITY VARIABLE GAUSS RULES

Figure 5.3. A “large” sample of Aluminum with 10240 embedded Boron particles. In reality the sample is approximately 0.1mm × 0.1mm × 0.1mm, while the particles are approximately 0.0035 mm in diameter

computational cost is O(n2 ) operations, and if the number of CG iterations is C, then the total cost is Cn2 operations, where typically C mloc . In this case the approximate cost savings are approximately N Mglob ≥ 512. The time to preprocess, DIRECT COST S/DECOM P OSED COST S ≈ mloc solve and postprocess each 20 particle subdomain took no more than one minute on a single RISC 6000 workstation. Clearly, the overall solution process is trivially parallelizable. We computed the following upper bound, which coincided with the previous chapter’s results, normalized by the energy of U = E · x, on the error

0≤

˜ ||2E(Ω) ||u − u ||U||2E(Ω)

≤

2(J (˜ u) − J − ) = 0.1073, ||U||2E(Ω)

(5.13)

!−1 −1 1 dΩ , |Ω| Ω IE 2 − ∂ J $ 2 < 0 ⇒ IR = ∂ IR

$ −1 where the maximized J − = 12 E : IE −1 −1 Ω : E|Ω|, IR = IE Ω =

has been used. The optimal IR$ stems from computing E −1 −1 Ω .

∂J − ∂ IR

= 0,

184


5.4 A “Total” Orthogonal Sum Clearly, for the class of problems under consideration, solutions must be generated numerically, for example, as in the last section, by the finite element method. We now consider the effect of the use of the finite element method to generate a subspatial approximation ˜ , denoted u ˜ h , governed per subdomain by to u ˜ hK ∈ H hu˜ (ΩK ) ⊂ H 1 (Ω), u ˜ hK |∂ΩK ∩(Ω∪Γu ) = U ∈ H 1 (Ω), such that Find u ˜ hK dΩ = ∇v K : σ f · v K dΩ + t · v K dA ΩK

∂ΩK ∩Γt

ΩK

(5.14)

∀v hK ∈ H hv (ΩK ) ⊂ H 1 (ΩK ), v hK |∂ΩK ∩(Ω∪Γu ) = 0.

A critical point is that H hu˜ (Ω), H hv (Ω) ⊂ H 1 (Ω). This “inner” approximation allows the development of a straightforward decomposition. Using a subspatial approximation, by direct expansion we have

J (˜ uh ) =

N K=1

=

N

JK (˜ uK ) +

+

N K=1

=

N

N

K=1

N

K=1

˜ K )), JK (˜ uK + (˜ uhK − u

K=1

K=1

−

N

JK (˜ uhK ) =

f·

(˜ uhK

˜ K ) : IE : ∇˜ ∇(˜ uhK − u uK dΩ,

ΩK

˜ K ) dΩ − −u

ΩK

1 2

N K=1

∂ΩK ∩Γt

˜ K ) dA, t · (˜ uhK − u

˜ K ) : IE : ∇(˜ ˜ K ) dΩ, ∇(˜ uhK − u uhK − u ΩK

J (˜ uK ) +

K=1

N 1 ˜ K ) : IE : ∇(˜ ˜ K ) dΩ, ∇(˜ uhK − u uhK − u 2 Ω K K=1

(5.15)

where, by Box 5.14, for K = 1, 2, ..., N ΩK

˜ K ) : IE : ∇˜ ∇(˜ uhK − u uK dΩ −

ΩK

˜ K ) dΩ − f · (˜ uhK − u

∂ΩK ∩Γt

˜ K ) dA = 0. t · (˜ uhK − u

Subtracting J (u) from both sides of Equation (5.15) leads to u) − J (u) + J (˜ uh ) − J (u) = J (˜

1 2

Ω

˜ h ) : IE : ∇(˜ ˜ h ) dΩ, ∇(˜ u−u u−u

(5.16)

(5.17)

which immediately implies ˜ h ||2E(Ω) = ||u − u ˜ ||2E(Ω) + ||˜ ˜ h ||2E(Ω) . ||u − u u−u

(5.18)

Therefore the estimate of the error, for example in Boxes 5.9 and 5.13 implicitly contains the numerical discretization error. Clearly, the total error is composed of two mutually


185

orthogonal parts: (1) a partitioning error introduced by subdividing the body and applying inexact boundary data on the interfaces and (2) a numerical error, introduced by spatial discretization. Therefore, if one has a numerical error estimator, one can isolate the pure decomposition error. Such “filtering” processes have been performed in Zohdi and Wriggers [97]. For details on elementary a priori numerical error estimates, see the text of Hughes [46], while for a posteriori estimates see Ainsworth and Oden [2]. 5.5 Extensions and Future Work In order to achieve further reduction of the global/local error, the subdomain interfacial boundary conditions must be improved. Clearly, the error bound in Box 5.9, and the orthogonality expression in Box 5.18 still holds, provided that the modified, hopefully improved, local boundary data deliver a globally kinematically admissible solution. Thus, we now discuss a few methods which consist of updating the numerical degrees of freedom on the interface of a finite element discretization. Without loss of generality, and for simplicity of presentation, we consider that we have meshed each subdomain in the substructuring process so that a global nodally-conforming trilinear hexahedral finite element mesh results. 5.5.1 Method I: global/local CG iterations

One can consider a process of updating interfacial nodal values as having as its goal the minimization of the discrete analog of J (˜ u), namely J (˜ uh ), which is convex in terms of the finite element nodal displacements. There exist a variety of classical techniques which could drive a process to minimize potential energies, the most applicable in the linearly elastic case being the CG method, which, for the sake of completeness, we briefly discuss in the context of a global/local strategy. Suppose we wish to solve the discrete system [K]{uh } = [B]. [K] is a symmetric positive definite n × n matrix, {uh } is the n × 1 solution vector, [B] is the n × 1 righthand side, and n is the number of discrete unknowns. We define a discrete potendef tial Π = 12 {uh }T [K]{uh } − {uh }T [B]. Correspondingly, from basic calculus we have def

∂Π ∂Π ∂Π T , , .... ∂u } = 0 ⇒ [K]{uh } − [B] = 0. Therefore the minimizer of the ∇Π = { ∂u 1 ∂u2 n potential Π is also the solution to the discrete system. The CG method is based upon minimizing Π by successively updating a starting vector. Therefore, the substructuring method described in this work can be considered as a construction of an advanced initial starting ˜ h , to the solution of the global problem formed by directly discretizing the entire guess, u body with no decomposition. If we are to directly apply the CG method, one will reduce J (˜ uh ), to the minimum value for the corresponding global mesh discretization. Since a CG approach updates all values in the entire domain simultaneously, via a global matrix-vector multiply, each iteration can be thought of as providing global transfer of information, i.e. ˜ h , and thus the interfacial boundary microstructural corrections to the initial solution u conditions. The minimum of the discrete potential is the global solution for the given finite element discretization. However, instead of a direct application of the CG method, one can perform a set number of global CG iterations, Mglob, then reevaluate the local boundn 2 ary value problems with the new local boundary data. Per subdomain we have mloc ( N ) operation counts, where mloc is the number of local CG iterations for subdomain solution convergence. For simplicity, we have assumed that one would use the CG method to solve each subdomain problem, and that mloc is approximately the same for each subdomain. n 2 Therefore, for N total subdomains we have mloc ( N ) N operation counts. Let us denote the number of times that we globally iterate and resolve the local problems as, S. Thus,

186


n 2 in total, we have S(Mglobn2 + N mloc ( N ) ) = S(Mglob + mNloc )n2 . Therefore, a direct CG approach and a projecting-global/local resolve approach compare as follows

mloc 2 ˆ 2 def Cn = S(Mglob + )n N

2 Cn

and

Conjugate Gradient

.

(5.19)

projecting and resolving

Clearly the magnitudes of Cˆ and C determine the relative performance of the two methods. 5.5.2 Method II: iterative equilibration

Clearly, when decomposing the structure by a projection of a kinematically admissible function onto the partitioning interfaces, regardless of the constitutive law, the error is due to the jumps in tractions at the interfaces (Box 5.4). If the interfaces would be in equilibrium, then there would be no traction jumps. Therefore, an approach which attempts to eliminate the jumps in tractions is sought next. Global equilibration For global equilibration one treats the entire interface as a single subdomain, denoted ω ⊂ Ω, and solves the following (Figure 5.4) ˜ ∈ H 1 (Ω), such that Find s ∈ H 1 (ω), s| ∂ω∩(ω∪Γu ) = u ∇v : σ(s) dω = f · v dω + t · v da ω

ω

∀v ∈ H 1 (ω), v|∂ω∩(ω∪Γu ) = 0.

(5.20)

∂ω∩Γt

The new interfacial displacements are then computed, and the subdomain problems are re-solved. The process is repeated until convergence. The entire interface problem can be quite large. A less expensive approach, amenable to trivial parallel processing, follows. Local (nodal) equilibration For local equilibration, one prescribes small non-intersecting “nodal domains” , ωK (Figure 5.4), surrounding each interfacial node 1 ≤ K ≤ M : (M = number of interfacial nodes) ˜ ∈ H 1 (Ω), such that Find sK ∈ H 1 (ωK ), sK |∂ω K ∩(ωK ∪Γu ) = u ∇vK : σK (s) dωK = f · v K dωK + t · v K daK ωK

ωK

∀v K ∈ H (ωK ), vK |∂ωK ∩(ωK ∪Γu ) = 0. 1

∂ωK ∩Γt

(5.21)

The new interfacial displacements are then computed, and the subdomain problems are resolved. The process is repeated until convergence. The advantage of such an approach is that the “nodal domains” can be processed independently of one another. The disadvantage is that global equilibrium is not satisfied throughout the entire interface. Either of the introduced approaches can be considered as a fixed-point iteration, which we discuss next.

187


GLOBAL INTERFACE PROBLEM IS SOLVED

LOCAL INTERFACE (NODAL) PROBLEMS ARE SOLVED

Figure 5.4. Global equilibration and local equilibration

General fixed-point iterations Consider a general system of coupled partial differential equations given by A(s) = F, where s is a solution, and where it is assumed that the operator, A, is invertible. One desires that the sequence of iterated solutions, sI , I = 1, 2, ..., converge to A−1 (F) as I → ∞. If sI is a function of A, F, sI , ...sI−K one says that K is the order of iteration. It is assumed that the I−th iterate can be represented by some arbitrary function sI = T I (A, F, sI−1 ). One makes the following split sI = G I (sI−1 )+rI . For this method to be useful the exact solution should be reproduced. In other words, when s = A−1 (F), then s = A−1 (F) = G I (A−1 (F)) + rI . Therefore, one has the following consistency condition rI = A−1 (F) − G I (A−1 (F)), and as a consequence, sI = G I (sI−1 ) + A−1 (F) − G I (A−1 (F)). Convergence of the iteration can be studied by defining the error vector: eI

=

sI − s = sI − A−1 (F) = G I (sI−1 ) + A−1 (F ) − G I (A−1 (F )) − A−1 (F )

=

G I (sI−1 ) − G I (A−1 (F)).

(5.22)

GI

eI

G I (sI−1 −A−1 (F))

One sees that, if is linear and invertible, the above reduces to = = G I (eI−1 ). Therefore, if the spectral radius of G I , i.e. the magnitude of its largest eigenvalue, is less than unity for each iteration I, then eI → 0 for any arbitrary starting solution sI=0 as I → ∞. An example d (E(x) du Again consider the simple one dimensional structure 0 < x < L : dx dx ) = 0, u(0) = ∆0 , u(L) = ∆L in Figure 5.5. Consider two subdomains, and one interfacial problem, which are solved in an iterative manner i = 1, 2, ...N , until |eI ( L2 )| = |sI ( L2 ) − u( L2 )| ≤ T OL

SUBD#1 (0 < x
α > − 12 . The quantities of interest in Box 6.11 are shown in Figure 6.6. 1

PERTURBATIONS

0.8

0.6

ENERGY-BULK PERT ENERGY-SHEAR PERT STRESS-BULK PERT STRESS-SHEAR PERT

0.4

0.2

0 -0.6

-0.4

-0.2

0

0.2 0.4 LOAD VARIABLE

0.6

0.8

Figure 6.6. Shear loading (φ = 0.5): perturbation (ηW ∗ −κ∗ , ηW ∗ −µ∗ , ησ−κ∗ , ησ−µ∗ ) with multiaxial loading

1

behavior

Remark 1: As the theoretical results assert, in all three loading cases the perturbations in the finite deformation range are bounded by the perturbation ratio in the infinitesimal range. Remark 2: An analyses for more sophisticated models, for example Ogden response functions, is quite more involved. It is being currently investigated whether the analysis presented here can be extended to obtain meaningful information for such hyperelastic functions, since they too are influenced by uncertainties in the linear elastic constants. Remark 3: Problems of geometric nonlinearities will not be considered further in this work, except in the appendix. For the rest of the monograph, emphasis is placed on material nonlinearities at infinitesimal strains, in particular the effects of damage. 7 MECHANICAL LOSSES IN MATERIAL RESPONSES A drawback in the use of solid structures composed of heterogeneous media is the existence of highly distorted and intense microfields, which may cause a significant amount of microscale failure in the form of subcontinuum scale microfissures, interface separation, void nucleation, etc. As the loading progresses, these local failures increase in size and number and eventually merge to produce observed macroscopic failure. Therefore, the desirable aggregate response of the heterogeneous material may be lost. In particular, of interest to a structural engineer is knowledge of the deterioration of IE ∗ with progressive loading, or overloading. The typical approach to describe the deterioration of the macroscopic response is to employ phenomenologically based macroscopic damage parameters, and corresponding evolution laws. These are applied directly on the macroscopic level in order to mimic the observed weakening (softening) responses of material samples with progressive loading. There are a variety of such approaches, loosely referred to as Continuum Damage

199


Mechanics (CDM) models. For an overview of the field see Krajcinovic [51]. Common to such approaches is that the mentioned evolution law is based on, or operates solely upon, macroscopic quantities such as IE ∗ . Typically, the evolution law is a curve fit of laboratory observed macroscopic behavior relating σΩ to Ω . Consequently, it is fair to state that such phenomenological models are not really predictive of responses that are unknown a-priori. The goal of this chapter is to develop a microscale model, amenable to direct numerical simulation, in order to describe the deterioration of aggregate responses of assemblages of heterogeneous material. Afterwards, if relations between volumetrically averaged quantities are desired, they can be post-processed from the microfield solution. The extent of deterioration, which is induced by reducing the eigenvalues of the elasticity tensor at a point in the heterogeneous body, is dictated by the condition that the solution must satisfy the equations of equilibrium, and simultaneously the constraints at that point. In a sense, the model is quite similar to plasticity formulations. However, the main difference is that the material properties change, as opposed to the generation of plastic strain. Here we concentrate only on the phenomena of microscopic deterioration, and exclude inelastic strains. However, the model is to be viewed as complementary to plasticity models, when they are applied on the microscale, and can be combined with them, although this will not be carried out in this work. For macroscopic approaches combining and comparing damage and plasticity models, we refer the reader to the representative articles of Carol et. al, [14], Ju [48], Lemaitre [55] and Yazdani and Schreyer [90]. Throughout the analysis the microscopic material is assumed to remain a continuum. Therefore, to an extent, this model is also phenomenological, but on the microscale, since no attempts are made to geometrically describe the true subcontinuum fissuring mechanisms that cause the changes in microscopic properties. The work follows from results found in Zohdi et al. [98] and Zohdi and Wriggers [94] and [97]. 7.1 A Micro-macro Mechanical Degradation Model We denote u as the undeteriorated “undegraded” solution field of the following variational boundary value problem: Find u ∈ H 1 (Ω), u|Γu = d, such that ∇v : IE : ∇u dΩ = f · v dΩ + t · v dA Ω

Ω

Γt

∀v ∈ H 1 (Ω), v|Γu = 0.

(7.1)

This formulation contains no description of microscopic deterioration, which we refer to in this work as “degradation”. We now construct a model for a degrading microstructure, by enforcing two physically motivated requirements for a weaking material undergoing irreversible changes: (1) under pure displacement boundary loading control, Γu = ∂Ω, the body should exhibit energy dissipation and (2) for pure traction loading control, Γt = ∂Ω, an energy increase. Energetic restrictions on microstructural degradation Consider two symmetric positive definite material property (elasticity tensor) distributions, the undegraded material, IE, and a deteriorated or degraded material IE deg . Both are separately used to generate solutions to two boundary value problems, with different material coefficients, but having the same exterior geometry and loading. The corresponding displacement (u), stress (σ) and strain ( = 12 (∇u+(∇u)T )) states produced when using these

200


materials are denoted (u, σ, ) and (udeg , σ deg , deg ) respectively. Consider their respective strain energies: def

Π =

Ω

: IE : dΩ

def

and

Πdeg =

deg : IE deg : deg dΩ .

Ω

undegraded energy

(7.2)

degraded energy

We have the following result for Γu = ∂Ω: (IE − IE deg ) ≥ 0 ⇒

Ω

: IE : dΩ −

Ω

deg : IE deg : deg dΩ ≥ 0.

Π

(7.3)

Πdeg

Similarly, defining def

Π =

Ω

σ : IE

−1

: σ dΩ

def

and

Πdeg =

undegraded energy

Ω

σdeg : IE −1 deg : σ deg dΩ,

(7.4)

degraded energy

we have δσ = σdeg − σ, if Γt = ∂Ω (IE −1 − IE −1 deg ) < 0 ⇒

σ : IE −1 : σ dΩ − σ deg : IE −1 deg : σ deg dΩ ≤ 0. Ω Ω

Π

(7.5)

Πdeg

The relations in Boxes 7.3 and 7.5 provide sufficient restrictions on the properties of a degraded microstructure, such that physically realistic macroscopic responses can be insured. For the sake of completeness, we prove such results. The proofs are constructive, and we believe add extra insight the analysis to follow later. 7.1.1 Proof of energetic degradation restrictions

Consider two symmetric positive definite material property (elasticity tensor) distributions, IE (I) and IE (II) used in same boundary value problem, for example Box 7.1. The corresponding stress and strain states using these materials are denoted ((I), σ (I)) and ((II), σ (II)) respectively. Consider their respective strain energies: def

ΠI =

(I) : IE (I) : (I) dΩ

Ω

def

and

ΠII =

I s energy

((I) + δ) : IE (II) : ((I) + δ) dΩ .

Ω

II s energy

We have, denoting δ = (II) − (I), ΠI − ΠII

(I)

=

(I)

: IE (I) :

dΩ −

Ω

−

(I) : IE (II) : ((I) + δ) dΩ Ω

δ : IE (II) : ((I) + δ) dΩ

Ω

=0 if Γu =∂Ω

(7.6)

201


(I)

=

: (IE (I) − IE (II) ) :

dΩ −

(I)

Ω

(I) : (IE (I) − IE (II) ) : (I) dΩ −

=

(I) : IE (II) : δ dΩ

Ω

Ω

((I) + δ) : IE (II) : δ dΩ

Ω

=0 if Γu =∂Ω

+

δ : IE (II) : δ dΩ

Ω

: (IE (I) − IE (II) ) : dΩ + δ : IE (II) : δ dΩ .

Ω

≥0 pos. def. for (IE (I) − IE (II) ) > 0 ∀ x (I)

(I)

=

(7.7)

Ω

The vanishing terms in the preceding analysis terms are due to a direct application of the principle of virtual work. The end result is that (IE (I) − IE (II)) ≥ 0 implies (I) : IE (I) : (I) dΩ − ((I) + δ) : IE (II) : ((I) + δ) dΩ ≥ 0.

Ω

Ω

I s

II s

energy

(7.8)

energy

Similarly, defining δσ = σ(II) − σ(I), if Γt = ∂Ω def

ΠI =

and def

ΠII =

Ω

Ω

(I) σ (I) : IE −1 dΩ, (I) : σ

(7.9)

I s energy

(I) (σ(I) + δσ) : IE −1 + δσ) dΩ, (II) : (σ

(7.10)

II s energy

and we have ΠI − ΠII

= Ω

(I) σ (I) : IE −1 dΩ − (I) : σ

δσ :

+

Ω

IE −1 (II)

σ

(I)

σ

(I)

:

(IE −1 (I)

:

(IE −1 (I)

Ω = Ω

−

Ω

=

Ω

(I) σ (I) : IE −1 + δσ) dΩ (II) : (σ

− δσ) dΩ

=0 if Γt =∂Ω

=

: (σ

(I)

−

IE −1 (II) )

−

IE −1 (II) )

:σ

(I)

:σ

(I)

dΩ − Ω dΩ +

σ (I) : IE −1 (II) : δσ dΩ

(σ(I) + δσ) : IE −1 (II) : δσ dΩ Ω

=0 if Γt =∂Ω

δσ : IE −1 (II) : δσ dΩ (IE −1 (I)

IE −1 (II) )

− : σ dΩ −

Ω

−1 −1 neg. def. for (IE (I) − IE (II) ) < 0 ∀ x σ

(I)

:

(I)

δσ : IE −1 (II) : δσ dΩ . Ω

≤0

(7.11)

202


−1 Therefore (IE −1 (I) − IE (II) ) < 0 implies

(I) (I) σ(I) : IE −1 : σ dΩ − (σ(I) + δσ) : IE −1 + δσ) dΩ ≤ 0. (I) (II) : (σ Ω Ω

I s energy

(7.12)

II s energy

Similar results, however for more restrictive (uniform) loading cases, and not in the context of a microfailure analysis, can be found in Huet, Navi, and Roelfstra [38]. The above results assert that a sufficient condition to meet our physical restrictions is to force the eigenvalues of IE deg to decrease at a location where failure occurs, and to remain unchanged otherwise. It is therefore convenient to define a spatially varying tensorial degradation function, α, which operates on the eigenvalues of an otherwise undegraded elasticity tensor: def

IE deg = α ) IE ⇒ EIG(IE deg ) = {Λdeg/1, Λdeg/2 , Λdeg/3 , Λdeg/4 , Λdeg/5 , Λdeg/6 },

(7.13)

= {α1 Λ1 , α2 Λ2 , α3 Λ3 , α4 Λ4 , α5 Λ5 , α6 Λ6 }, where 0 < α1→6 ≤ 1. The simplest representation, which is the one we employ in this work is IE deg = αIE

0 < α ≤ 1.

with

(7.14)

The scalar degradation function α takes on different values throughout the body, which are dictated by the solution to a so-called degraded boundary value problem introduced next. Remark: When the material is isotropic, the eigenvalues are (3κ, 2µ, 2µ, µ, µ, µ). Therefore, we must have κ > 0 and µ > 0 to retain positive definiteness of IE, and consequently ακ > 0 and αµ > 0 for IE deg . 7.1.2 Degraded formulation

Consistent with the previous results, to describe losses of stiffness in the material on the microscopic level, we construct a degraded solution, udeg , that is generated by a weakened material modulus, IE deg ≤ IE. The corresponding variational boundary value problem governing the degraded solution is: Find udeg ∈ H 1 (Ω), udeg |Γu = d, such that ∀v ∈ H 1 (Ω), v|Γu = 0 ∇v : IE deg : ∇udeg dΩ = f · v dΩ + t · v dA Ω

If M < K then IE deg = IE

Ω

Γt

(7.15)

(no degradation)

If M = K then 0 < IE deg = αIE ≤ IE

(degradation).

The functions M = M(IE deg , udeg , α) and K = K(IE deg , udeg , α) respectively serve as a measure of, and constraint on, selected internal fields of interest to the analyst. The functions M and K could be scalar, vector, or tensor valued. In this work we consider only scalar constraints, of which specific forms are introduced later. The values of α are dictated by the fact that the solution udeg must satisfy the equations of equilibrium, and simultaneously the constraints.

203


7.1.3 One dimensional examples

Consider a simple one-dimensional bar composed of a material that is initially linearly elastic, at infinitesimal strains. For a displacement tension test load: udeg (0) = 0 and udeg (x = L) = E × L, where E is an arbitrary constant that ! represents “applied strain”, the dudeg dudeg d corresponding boundary value problem is dx Edeg dx = 0, which implies dx =E ⇒ σdeg = Edeg E, where Edeg = αE is a material which may degrade (Figure 7.3). To satisfy a fixed stress constraint, such as σdeg ≤ σ crit , the material must obey (Figure 7.1): (1) if crit σdeg > σ crit , then Edeg E = σ crit and Edeg = σ E or (2) if σdeg ≤ σ crit then Edeg = E. In this def σ crit σdeg

case, the degradation function (α) obeys (1) if σdeg > σ crit then α =

or (2) if σdeg ≤

def

σ crit then α = 1. A relatively general constraint equation, for the case when the constraint surface is dependent on the deformation, is σ crit = σ lim +(σ0crit −σ lim )αP (Figure 7.2), where σ0crit represents an initial critical stress which must be met to initiate degradation, and where σ lim is a limit (asymptotic) stress. It is clear, in an algorithmic sense, if σdeg − σ crit > 0, then we must enforce σdeg = σ crit , which implies αE deg − σ lim − (σ0crit − σ lim )αP = 0. For simplicity, consider the case when σ lim = 0, and therefore αEE − σ0crit αP = 0, which

yields α =

1 P −1

EE σ0crit − σ0crit )

=

σ0crit EE

1 1−P

. Clearly, a critical value is P = 1. For this case we

have α(EE = 0. Therefore, since for a fixed σ0crit , in general (EE − σ0crit ) += 0, therefore α = 0. This violates the restriction that α > 0, and thus the constraints cannot be satisfied, implying that the boundary value problem has no solution. For P > 1, for example P = 2 ⇒ |α| > 1, which is physically absurd behavior. No solution exists for P ≥ 1, that satisfies the constraint conditions and equilibrium. For these reasons, we consider ranges of P < 1 in this work, even though in some cases, for nonzero values of σ lim , it may be possible to attain a solution for P > 1. To give an idea of the solution behavior, for the one dimensional bar we plotted the stress strain curve for various values of P < 1, corresponding to σlim