Computation of Homogenized Constitutive Tensor of Elastic Solids Containing Evolving Cracks

Abstract

The determination of the equivalent (homogenized) constitutive tensor is one of the most important steps in multiscale models as well as in the classical homogenization theory. In this article, a procedure for determining the homogenized instantaneous (tangent) constitutive tensor of elastic materials containing growing cracks is proposed. The primary purpose of this procedure is its use in two-way coupled multiscale finite element algorithms that can model crack formation and propagation at the local microstructure. The procedure is basically developed by relating the local displacement field to the global strain tensor at each location and using first-order homogenization techniques. The finite element formulation is developed and some example problems are presented in order to verify and demonstrate the model capabilities.

Keywords

homogenization constitutive tensor elastic media evolving cracks multiscale modeling

INTRODUCTION

Numerous materials used in structural applications are heterogeneous, at least at some microscopic level. Examples of heterogeneous materials with structural functionality are Portland cement concrete, asphalt concrete, steel, geological materials, carbon fiber-reinforced composites, and human bone. The use of such materials is mainly due to the fact that they can be suitably designed to take advantage of particular properties of each individual constituent. However, the design process requires a reasonable physical understanding of the interactions among all constituents, as well as the eventually evolving internal boundaries (voids and cracks), and its effect on the overall behavior of the material.

The classical homogenization theory and multiscale models are the most commonly used approaches to modeling and design composite materials. The objective of the classical homogenization theory is the determination of the effective (homogenized) constitutive behavior of heterogeneous materials by solving an initial boundary value problem (IBVP) for a representative volume element (RVE) (Eshelby, 1957; Hill, 1963, 1965a,b; Hashin, 1964; Fish et al., 1997; Allen and Yoon, 1998). Once the homogenized constitutive properties have been determined a priori, modeling of more complex problems can be performed. On the other hand, multiscale models intend to determine the effective constitutive behavior of the RVE simultaneously throughout the analysis (Allen, 1994, 2001; Feyel, 1999; Feyel and Chaboche, 2000; Fish and Shek, 2000; Souza et al., 2008), using concepts of the classical homogenization theory, either using a self-consistent approach (Eshelby, 1957; Hashin and Shtrikman, 1963), or the so-called asymptotic analysis (Bensoussan et al., 1978; Sanchez-Palencia, 1980; Bakhvalov and Panasenko, 1984).

Sometimes, multiscale models are termed two-way coupled because the loads (tractions or deformations) observed at the global-scale object are applied to the boundaries of RVE's (global–local coupling), whereas the homogenized material properties are passed from the local-scale RVE's to the global-scale material coordinates (local–global coupling). Therefore, two-way coupled multiscale models are more appropriate when the problem considers evolving microstructure, such as growing cracks, since the evolution of cracks generally depend on the loading history, especially for inelastic materials. In both classical homogenization theory and multiscale models, however, the determination of the effective constitutive behavior of the material is crucial.

In practical applications, most structural parts experience at least a certain amount of microcracks during their intended life and the existence of microcracks may not adversely affect their structural performance if the amount of microcracks is kept up to a certain level so that they do not coalesce and form a macrocrack that leads to component failure. In order to facilitate the use of two-way coupled multiscale models to design and model such materials with microcracks, this article presents a procedure to determine the effective constitutive tensor of elastic solids containing growing cracks. The same approach has already been used to model similar problems but involving history-dependent viscoelastic materials (Souza, 2009; Souza and Allen, 2010a).

MULTISCALE MODELING

The purpose of multiscale models is to perform simultaneous analyses on both global and all local scales. The effective constitutive behavior of one scale is then determined by the behavior of each consecutively smaller scale in such a manner that only the material properties of the individual constituents considered at the smallest length scale, as well as the fracture properties, need to be known a priori.

In cases where the local structure (RVE) does not change and its behavior does not depend on the loading history so that the effective constitutive properties need to be determined only once, e.g., purely elastic solids, the classical homogenization theory is more cost effective. However, multiscale models are particularly advantageous for problems where the local structure may change due to the formation and growth of internal boundaries because the damage-dependent homogenized constitutive properties need to be recursively determined in each step of the analysis.

Consider the two-scale IBVP schematically shown in Figure 1. The global object, referred to as scale 0, is considered to be statistically homogeneous whereas the local object, referred to as scale 1, is locally heterogeneous containing inclusions and cracks which may be growing or not. In Figure 1, superscripts are used to denote the scale index, while subscripts refer to the rectilinear coordinate system. V^μ, $\partial V_{E}^{μ}$ , and $\partial V_{I}^{μ}$ denote, respectively, the volume, the external, and the internal boundary surfaces of the object at scale μ. The length scale associated with the object and with the cracks at scale μ are denoted by ℓ^μ and $ℓ_{c}^{μ}$ , respectively.

Figure 1.

Schematic representation of a two-scale IBVP.

The mechanical IBVP for the global object can then be well posed by appropriate initial boundary conditions and a set of equations: conservation of linear and angular momentum (Equations (1) and (2)), strain–displacement relations (Equation (3)) and constitutive equations (Equation (4)).

(1)

where

σ_{ij}^{0}

, ρ⁰,

b_{i}^{0}

, and

u_{i}^{0}

are, respectively, the Cauchy stress tensor, the mass density, the body force vector per unit mass, and the displacement vector at a certain location

x_{i}^{0}

inside the statistically homogeneous object at the global length scale. The right-hand side of Equation (1) can be neglected in the case where the global object is subject to quasi-static loading.

(2)

(3)

where

ɛ_{ij}^{0}

is the infinitesimal strain tensor defined on the global scale.

(4)

where

{\bar{Ω}}_{τ = - \infty}^{τ = t}

is a functional describing the constitutive behavior at each position

x_{i}^{0}

in the object at the global scale, which may account for history-dependent effects, such as viscoelasticity, and is obtained by locally averaging the constitutive behavior of the RVE corresponding to that particular position. In this study, history-dependent effects are not considered.

The IBVP for the local RVE should obey the same set of governing equations, except for constitution, as the global object. The local boundary conditions should be assessed based on the actual traction (or deformation) field determined for the global scale mapped into the local coordinate system $x_{i}^{1}$ . However, this could be very complicated in practice so that some hypotheses are generally assumed in order to simplify the problem. First, it is assumed that the global length scale is much larger than the local length scale:

(5)

This assumption allows the use of homogeneous boundary conditions (uniform boundary tractions or linear boundary displacements) at the local scale, also simplifying the homogenization relationships connecting the global and local scales. The second assumption is that the local length scale is much larger than the length scale associated with cracks at the local scale, thus limiting the model to problems where statistical homogeneity is satisfied:

(6)

Finally, for the case where the global problem cannot be simplified to a quasi-static one, it is assumed that the length of the wave propagating at the global scale, $ℓ_{w}^{0}$ , is much larger than the local length scale:

(7)

In that limiting case, the local boundary conditions can be approximated by homogeneous boundary conditions, mitigating the need to consider waves propagating at the local scale, thus allowing the local IBVP to be approximated as a quasi-static problem:

(8)

(9)

(10)

(11)

(12)

where Equation (12) describes a fracture criterion,

G_{i}^{1}

is the fracture energy release rate at a particular position in the local scale,

G_{ci}^{1}

the critical energy release rate of the material at that particular position, and the index i the mode of fracture. Furthermore, the functional

Ω_{τ = - \infty}^{τ = t}

describing the constitutive behavior of each particle at the local position

x_{i}^{1}

is known a priori for all constituents.

For the cases where the previous assumptions are valid, one can show that the relations connecting the global and local scales are given by the following (mean field) homogenization principles (Allen and Yoon, 1998; Nemat-Nasser and Hori, 1999; Allen, 2001; Souza et al., 2008):

(13)

(14)

(15)

(16)

(17)

where

n_{i}^{1}

is the outward unit normal vector to the external or internal boundary,

\partial V_{E}^{1}

\partial V_{I}^{1}

, of the RVE and it has been assumed that the body forces acting on the local scale,

b_{i}^{1}

, are negligible in order to obtain the second equality in Equation (16). Observe from Equation (14) that the strain tensor in the global scale is equal to the external boundary average of the displacement vector on the local scale, which differs from the volume average of the strain tensor by the internal boundary average of the displacement vector,

α_{ij}^{0}

. Also note that in the case where there is no separation along the internal boundaries of the RVE (no voids or cracks)

α_{ij}^{0}

is zero, so that

α_{ij}^{0}

can be regarded as a measure of damage at the global scale (Allen and Yoon, 1998).

On the other hand, from Equation (16), one can conclude that under the given circumstances, the external boundary average of the traction vector, $t_{i}^{1}$ , is equal to the volume average of the stress tensor. However, this is true only if the traction vector over the internal boundary surfaces is zero (free surfaces, cracks) or self-equilibrating (cohesive zones). Recall from Cauchy's formula that the traction vector on any surface in V^μ is given by:

(18)

The use of Equations (14) and (16) as the state variables at the global scale is commonly termed a mean field theory of homogenization because the behavior of the global object is determined only in terms of the mean stress and strain tensors evaluated at the external boundary of the local RVE. Higher order homogenization theorems can, however, be formulated if one includes gradients of the deformation and moments of the stress tensor in the model formulation (Allen et al., 1987a,b; Phillips et al., 1999; Kouznetsova, 2002; Kouznetsova et al., 2004).

If the global-scale problem is dynamic and one adopts an explicit time integration scheme to solve the global IBVP, the above equations are sufficient to perform the multiscale analysis, because there is no need to compute a stiffness matrix. Therefore, for each time step, one can simply apply displacements on the boundary of the RVE that satisfy Equation (14), solve the step-wise local IBVP, and then compute the global stress tensor, $σ_{ij}^{0}$ using Equation (16), even though it is always desirable to compute the homogenized constitutive tensor.

However, if the global IBVP is approximated by a quasi-static problem, the computation of the homogenized constitutive tensor is required in order to correctly calculate the tangent stiffness matrix. If one uses an incremental solution scheme, it is then preferable to seek an incremental constitutive equation of the form:

(19)

COMPUTATION OF INCREMENTAL (TANGENT) HOMOGENIZED CONSTITUTIVE TENSOR

Since incremental algorithms are most commonly used for solving IBVP in general purpose Finite Element (FE) codes, this study seeks to develop a procedure to compute the incremental (tangent) homogenized constitutive tensor of elastic solids containing cracks. The model to be demonstrated in this article does not consider any kind of history-dependent constitutive behavior, neither for the bulk material nor for the cohesive zones ahead of the crack tips that may eventually develop. Readers interested in the use of such approach for viscoelastic media are referred to Souza (2009) and Souza and Allen (2010a).

At least to the authors' knowledge, only purely numerical techniques have been proposed in the literature so far to compute the homogenized tangent constitutive tensor. For example, Kouznetsova (2002) employs the technique of condensation of the total RVE stiffness, which has been presented in Kouznetsova et al. (2001). Feyel and Chaboche (2000) use the approximate numerical differentiation technique, which requires the solution of four (2D) or six (3D) IBVP's for every solution step. The approach to be proposed in this article, on the other hand, is based on continuum mechanics homogenization principles and allows the computation of the full homogenized anisotropic tangent constitutive tensor of elastic solids with evolving cracks by solving the IBVP only once, which greatly reduces the required amount of computational time.

First, assume that the stress tensor at both global and local scales can be written in incremental form:

(20)

(21)

Then, substituting the above equations into Equation (16), gives:

(22)

(23)

Assume that the constitutive behavior of the global object can be approximated by the following incremental form:

(24)

where

C_{ijkl}^{0} (t)

is the instantaneous (tangent) constitutive tensor evaluated at the previous time step t which is a function of time through its dependence on the amount of damage accumulated at the local RVE, thus producing a nonlinear behavior at the global scale. Also assume that the bulk materials at the local scale follow a simple linear elastic model which can be written in incremental form as:

(25)

Substituting Equations (24) and (25) into Equation (23), gives:

(26)

where

C_{ijkl}^{0} (t)

is now written as

C_{ijkl}^{0}

for the shortness of notation. Let us now relate the field variables of the local scale problem to the external boundary average quantities corresponding to field variables of the global scale. Allen and Yoon (1998) have used a fourth-order strain localization tensor which relates

ɛ_{kl}^{1}

ɛ_{kl}^{0}

. This approach may be sometimes impractical or inconvenient, especially when using the FE method to solve the local IBVP. Since the displacement vector is generally used as the primary variable in FE formulations, it is therefore preferable to relate the local displacement field, instead of the local strain field, to the global strain tensor (external boundary average of the local displacement field) according to:

(27)

or in an equivalent incremental form (for elastic materials):

(28)

where

λ_{ijk}^{1}

and

Λ_{ijk}^{1}

are herein called the localization and incremental localization tensors (Hill, 1967; Allen and Yoon, 1998). Note that the first term in Equations (27) and (28) corresponds to the homogeneous displacement field that would be experienced by the RVE if it was a homogeneous material. The second term, on the other hand, represents the contribution to the local displacement field due to the existence of local asperities and internal boundaries, including cracks and cohesive zones.

Oskay and Fish (2007) have used a similar decomposition of the local displacement field in order to analyze failure of heterogeneous materials using continuous damage mechanics to model both bulk degradation and fiber-matrix debonding.

Substituting (28) into (10), then into (26) yields:

(29)

where symmetry of the local constitutive tensor

C_{ijkl}^{1}

and the fact that

Δ ɛ_{kl}^{0}

is constant in V¹ have been used. Since

Δ ɛ_{kl}^{0}

is arbitrary, one obtains:

(30)

Note that the effect of heterogeneity, internal boundaries and cracks at the local scale on the global constitutive behavior is accounted for by the terms involving the incremental localization tensor, $Λ_{ijk}^{1}$ , even though it may not be instantly obvious. It will, however, become more clear when the FE formulation is developed in the following section.

FE FORMULATION OF THE LOCAL-SCALE IBVP

The FE formulation of the local IBVP is developed in this section. The reader should note that the global IBVP remains the usual problem which can be solved by commonly used FE algorithms. The local-scale IBVP, however, has been changed due to the introduction of the incremental localization tensor, $Λ_{ijk}^{1}$ . Consider the variational form of conservation of linear momentum at the local scale, Equation (8), neglecting the body force vector:

(31)

From the chain rule of differentiation and conservation of angular momentum, Equation (9), the above equation results in:

(32)

Applying the divergence theorem to the above equation at time t + Δt and using Cauchy's formula (Equation (18)) results in:

(33)

Besides the incremental Equation (21), consider that the displacement, strain, and traction fields can also be approximated by the following incremental forms:

(34)

(35)

(36)

Thus, substituting into (33), but noting that the variation of $ɛ_{ij}^{1} (t)$ and $u_{i}^{1} (t)$ are null because these are known quantities (determined in previous step), yields:

(37)

Also noting that,

(38)

results in:

(39)

Now consider only the first integral on the right-hand side of the equation above. This integral concerns with the loading applied on the external boundary of the local RVE. If displacements are applied, which is the case, then, the first integral on the right-hand side of the equation above is zero because the variation of the displacement field on the external boundary is null. Therefore,

(40)

Noticing that the virtual work done by tractions over the internal boundary (right-hand side of equation above) is zero unless there is separation along that boundary, the integral over the internal boundaries can be reduced to an integral over the cohesive zones (note that tractions are null over crack surfaces):

(41)

Substituting (28) into (10), then into (25), gives:

(42)

Also, assume that the following incremental traction–displacement relation holds for each cohesive zone at the local scale:

(43)

where the brackets in

Δ [u_{j}^{1}]

represent the displacement jump across the cohesive zone, which, from Equation (28), can be written as:

(44)

(45)

where one should note that the initial coordinates,

x_{i}^{1}

, of the cohesive zone surfaces coincide.

The incremental localization field can now be discretized in space by expressing it as a function of its nodal values:

(46)

where the above equation has been written in matrix form. [N¹] is the matrix of nodal shape functions and the over-tilde indicates nodal value. It is important to note that the shape functions used to interpolate Λ¹ are the same as for the displacement field. Therefore, one can also write:

(47)

(48)

(49)

where [B¹] is the matrix of derivatives of the shape functions with respect to the spatial coordinates. Substituting Equations (42) – (49) into Equation (41) and noting that

δ {\tilde{Λ}}_{ijk}^{1}

and

Δ ɛ_{ij}^{0}

are arbitrary quantities, results in:

(50)

where

(51)

(52)

Once $Λ_{ijk}^{1}$ has been obtained using Equation (50), one can compute the global homogenized incremental (tangent) constitutive tensor, $C_{ijkl}^{0}$ , from Equation (30) and calculate the incremental displacements using Equation (28). Then, one can use standard FE procedures to compute the incremental strain and stress tensors from the incremental displacements, thus minimizing the code implementation effort.

Now, one can clearly see that the incremental localization tensor, $Λ_{ijk}^{1}$ , is affected by cracks and cohesive zones through the stiffness matrix, [K¹]. Therefore, one can also observe how the global homogenized incremental constitutive tensor, $C_{ijkl}^{0}$ , is affected by cracks and cohesive zones through $Λ_{ijk}^{1}$ (Equation (30)). It is important to note that the stiffness matrix, [K¹], is the same stiffness matrix usually assembled in common FE algorithms, what greatly simplifies the implementation of the proposed algorithm.

CODE VERIFICATION

In order to verify the code implementation, consider the simple case where the local IBVP is approximated by two 4-noded quadrilaterals with one cohesive zone element in between, as shown in Figure 2. The cases of pure normal and pure shear loading are considered. In the first case, boundary conditions are set such that $Δ ɛ_{22}^{0}$ is the only nonzero component of the incremental external boundary average strain tensor. In the second case, boundary conditions are set such that $Δ ɛ_{12}^{0}$ is the only nonzero component of the incremental external boundary average strain tensor. In the first case, the component $C_{2222}^{0}$ of the homogenized (tangent) constitutive tensor is evaluated, while $C_{1212}^{0}$ is evaluated in the second case.

Figure 2.

Verification problem: two 4-noded quadrilateral elements with one cohesive zone element.

Even though a variety of cohesive zone model may be used, in this section, the cohesive zone is assumed to follow the Tvergaard model (Tvergaard, 1990) which is neither history nor rate dependent:

(53)

where

σ_{i}^{max}

is a material property of the cohesive zone, δ_i material length scales associated with debonding, the index i normal or tangential directions, and λ a normalized quantity coupling normal and tangential responses:

(54)

where n indicates normal direction and r and s tangential directions. It is further assumed that the cohesive zone is fully debonded when λ ≥ 1, in which case t_n = t_r = t_s = 0. Figure 3 depicts the normal traction as a function of the normal opening displacement, where both

σ_{i}^{max}

and δ_i have been set to unity. Moreover, as described in Tvergaard (1990), if λ > 1/3 and

\overset{\cdot}{λ} <0

, the following traction–displacement relationships are used to avoid physically unreasonable behavior:

(55)

where λ_max is the maximum value of λ prior to unloading and the above equations remain valid until reloading, λ > λ_max.

Figure 3.

Tvergaard traction–displacement relation for pure normal separation.

Now, the terms k_ij in Equation (43) can be obtained by differentiating Equation (53) with respect to the displacement jumps, yielding zero off-diagonal terms and the following for the diagonal terms:

(56)

where δ_ij is the Kronecker delta.

First, let the bulk material be infinitely rigid, E → ∞, so that the response of the body shown in Figure 2 will be governed by the response of the cohesive zone, which analytical solution for the incremental (tangent) constitutive behavior is given in Equation (56). The actual numerical values used for the Young's modulus and Poisson's ratio are 10¹² and 0.3, respectively. The cohesive zone material properties used are δ_n = δ_r = δ_s = 1.0 and $σ_{n}^{max} = σ_{r}^{max} = σ_{s}^{max} = 10^{3}$ . Plane stress conditions are assumed.

Figure 4 shows the computed homogenized incremental (tangent) moduli, $C_{2222}^{0}$ and $C_{1212}^{0}$ , for the case where the bulk material is infinitely rigid compared with the expected analytical solution for the cohesive zone itself. It is important to note that $C_{2222}^{0}$ and $C_{1212}^{0}$ assume negative values, because it is a tangent modulus, as opposed to a secant modulus. Figure 4 also shows the results for E = 10⁴, in which case the bulk modulus is of same order of magnitude as the cohesive zone stiffness, so that the homogenized behavior is now governed by the interaction between the bulk and cohesive zone materials. The effect of different bulk Young's modulus on the homogenized stress–strain behavior is presented in Figure 5. It is important to mention that the results obtained for the normal opening mode have been extracted from Souza and Allen (2010a).

Figure 4.

Comparison of analytical and numeric results for normal and shear modes. Results shown in part have been extracted from Souza and Allen (2010a).

Figure 5.

Normal and shear homogenized constitutive behavior for two values of Young's modulus. Results shown in part have been extracted from Souza and Allen (2010a).

One can also observe from Figure 6 that a damage-induced anisotropy is introduced in the material due to the cohesive zone debonding. This is an important phenomenon that occurs in certain applications where materials undergo crack propagation in preferable directions. Therefore, models that can accurately capture this behavior are desirable when modeling and designing the materials used in these applications.

Figure 6.

Damage induced anisotropy for E = 10⁴.

EXAMPLE PROBLEMS

In this section, two example problems are presented in order to demonstrate the model. The first example consists of a multiscale tapered bar problem using a simple local structure while the second problem considers the homogenization of fiber-reinforced composites containing cracks. These are illustrative examples using arbitrary geometries and material properties which are shown here in an attempt to provide a better understanding of the intended goal of the presented procedure for computing the homogenized incremental (tangent) constitutive tensor.

Multiscale Tapered Bar

Consider the tapered bar shown in Figure 7 as the global-scale object of interest. As can be observed, one has taken advantage of the geometric symmetry in order to minimize computational effort. Horizontal displacements are applied at the right end of the bar in a quasi-static fashion. The local structure is that shown in Figure 2 rotated by 90° so that the cohesive zone is perpendicular to the loading direction. A local structure of same initial geometry is used at all global scale integration points. Same material properties have been used for the cohesive zone and the bulk Young's modulus and Poisson's ratio are 10⁴ and 0.3, respectively. Plane stress condition is assumed for both global and local scales.

Figure 7.

Geometry and FE mesh of tapered bar problem: (a) global scale, (b) local scale.

The reader should notice that this is actually a one-dimensional problem, even though the problem has been simulated using two-dimensional meshes. The field variables do not change with respect to the vertical spatial coordinate. Therefore, the microcracks do not interact with the boundaries.

Figure 8 shows the contour of the component $C_{1111}^{0}$ of the computed homogenized incremental (tangent) constitutive tensor as well as selected local structures at the end of the simulation when the total applied tip displacement is of 2.35. The bulk elements at the local structures are colored in gray, while the cohesive zone is colored in dark gray in order to facilitate visualization and the used colors are not related to any contour color. Figure 9 presents how the homogenized modulus $C_{1111}^{0}$ varies with the horizontal position x along the symmetry line of the tapered bar. Values shown in the graph also correspond to the end of the simulation.

Figure 8.

Contour of $C_{1111}^{0}$ and local structures at different positions.

Figure 9.

Component $C_{1111}^{0}$ of homogenized incremental (tangent) constitutive tensor as a function of $x_{1}^{0}$ coordinate.

As expected, one can observe from Figure 8 that the displacement jump across the cohesive zone is larger at regions of smaller cross-section, where stresses are concentrated. Consequently, $C_{1111}^{0}$ is lower at smaller cross-sections (Figures 8 and 9), once more damage has been induced at the corresponding local structures. It is important to note that the local structures experience a homogeneous state of stress so that statistical homogeneity is not lost even in the presence of the cohesive zone. Also no traction-free cracks have developed in this example. The reader should notice that this example problem has been oversimplified for demonstration purpose.

Figure 10 presents the global stress component $σ_{11}^{0}$ at the right-hand side tip of the bar as a function of the global applied tip displacement. Results for both damaged and undamaged cases are shown, where in the undamaged case the tangent constitutive tensor, $C_{ijkl}^{0}$ , is not corrected according to the proposed procedure, thus producing results as if no cracks are allowed to initiate and grow at the local scale. As expected, the existence of microcracks considerably affects the behavior of the global object. One can also observe the nonlinear behavior of the global stress when microcracks are considered, as opposed to the linear behavior obtained when there is no microcrack. The model herein presented has been developed for use in such cases where microcracks, voids, and/or defects are important design variables.

Figure 10.

Comparison between predicted results with and without damage corrected tangent constitutive tensor.

Homogenization of Fiber-reinforced Composites with Cracks

Fiber-reinforced composites have a very broad range of application in real structures. These materials are, in general, filled with microcracks initiated by service loading or small defects due to the manufacturing process. It is, therefore, interesting to evaluate how the material properties change as a function of crack density.

Assume that the microstructure shown in Figure 11 corresponds to the RVE of an arbitrary fiber-reinforced composite. The fibers diameter has been arbitrarily defined as 1.0 µm. The fibers Young's modulus and Poisson's ratio are 100 GPa and 0.3, respectively, while 10 GPa and 0.3 are used as the matrix material properties. It is important to observe that all material properties herein used may not represent actual properties of fibers and matrix materials commonly used in industry.

Figure 11.

Initial geometry and FE mesh of RVE.

The FE mesh contains initially no cohesive zone elements, which are automatically inserted into the mesh once the crack initiation criterion is satisfied; in this case, if the traction vector along the edges between elements exceeds the maximum cohesive stress, $t_{i} \geq σ_{i}^{c}$ . However, for this example, cracks are not allowed to initiate inside the fibers. The cohesive zone insertion algorithm used has been implemented based on that introduced by Camacho and Ortiz (1996) and Pandolfi and Ortiz (2002). However, no coupling between normal and shear tractions is considered, as opposed to Camacho and Ortiz (1996).

This approach is not based on mesh refinement; instead, cohesive zone elements are inserted on the fly along the existing edges of the original FE mesh, by simply doubling the nodes that violated the crack initiation criterion. More details can be found in Souza and Allen (2010b).

The cohesive zone constitutive behavior is modeled by an irreversible linear decreasing law, depicted in Figure 12 and described in Camacho and Ortiz (1996) and Zhou and Molinari (2004), for which the fracture energy release rate for each mode of fracture, $G_{i}^{c}$ , is given by:

(57)

Figure 12.

Linear decreasing cohesive zone constitutive behavior.

The arbitrarily selected cohesive zone material properties are $σ_{n}^{c} = σ_{r}^{c} = σ_{s}^{c} =$ 10 MPa and $δ_{n}^{c} = δ_{r}^{c} = δ_{s}^{c} =$ 0.1 µm. Biaxial displacements are applied to the boundaries of the RVE resulting in external boundary average strain rates of ${\overset{\cdot}{ɛ}}_{22}^{0} = 1.0 \times 10^{- 4}$ and ${\overset{\cdot}{ɛ}}_{11}^{0} = 0.5 \times 10^{- 4}$ 1/s in the vertical and horizontal directions, respectively.

Figure 13 presents how the components of the homogenized incremental (tangent) constitutive tensor change as a function of the applied external boundary average strain rate ${\overset{\cdot}{ɛ}}_{22}^{0}$ . Snapshots of the RVE for selected states of strain are shown in Figure 14 from which one can observed the amount of cracking occurring at each stage. These selected states of strain are marked in Figure 13 with vertical dashed lines so that one can correlate how the components of the homogenized incremental (tangent) constitutive tensor change as a function of crack density by comparing Figures 13 and 14.

Figure 13.

Components $C_{1111}^{0}$ , $C_{2222}^{0}$ , $C_{1122}^{0}$ , and $C_{1212}^{0}$ of homogenized incremental (tangent) constitutive tensor as a function of the strain rate ${\overset{\cdot}{ɛ}}_{22}^{0}$ .

Figure 14.

Snapshots of the damage accumulated in the RVE.

It should be emphasized that the components of the homogenized tangent constitutive tensor $C_{ijkl}^{0}$ obtained in this example may be different for different loading conditions since the crack pattern may likely change.

It is important to say that this multiscale approach is computationally intense due to the complexity of the problem: multiple length scales and multiple evolving cracks. For example, a single processor at 2.16 GHz running under Fedora Linux 8 took 33 min to solve 76 time steps of the single RVE problem considered in this section. On the other hand, an alternative direct approach, wherein a highly refined mesh is used at a single global scale, is well beyond the capability of currently available computing platforms.

In order to minimize computational effort, one may, however, use multiscaling only at selected regions of the global structure where damage evolution is prone, thus reducing the number of local meshes. Another very efficient strategy is to use parallel computing in such a way that multiple processors can handle different local meshes simultaneously (Feyel and Chaboche, 2000; Souza and Allen, 2010a). For example, a plate impact problem of 600,000 time steps and using 48 local meshes (each with 1660 elements and 871 nodes) took 101 h (4.2 days) to be solved in parallel by eight processors at 3.0 GHz running under Red Hat Linux 5.0. The same problem took more than 4 weeks to be solved by a single processor. The authors believe that the advent of faster computers, the advance of parallel computing and the development of faster algorithms will improve feasibility of multiscale simulations of real structures soon.

CONCLUSIONS

This article has presented a procedure to compute the homogenized incremental (tangent) constitutive tensor of elastic solids containing cracks using the FE method. The case of viscoelastic solids has been reported in Souza and Allen (2010a). This procedure can be useful for both homogenization and multiscale models which account for cracking at the local levels. One can find applications for this procedure not only for composite materials, but also for homogeneous materials containing microcracks.

Some example problems have been presented in order to verify the correctness of implementation and to demonstrate the capabilities of the model, such as accounting for anisotropy and stiffness degradation induced by crack initiation and propagation.

These features adjoined to the ability of explicitly include important microstructural design variables, such as volume fraction of inclusion, combination of individual materials properties, spatial and orientation distribution of inclusions, defects and microcracks, are important reasons for the potential use of two-way coupled multiscale models in the analysis and design of heterogeneous materials.

Footnotes

FUNDING

This study was supported by the US Army Research Laboratory [grant number W911NF-07-D-0001] and the Brazilian National Council of Scientific and Technological Development – CNPq/Brazil [grant number 200372/2006-8].

References

Allen

D.H.

(1994). Damage Evolution in Laminates. In: Talreja

(ed.), Damage Mechanics of Composite Materials, Composite Materials Series, 9, Amsterdam, Elsevier, pp. 79–114.

Allen

D.H.

(2001). Homogenization Principles and Their Application to Continuum Damage Mechanics, Composites Science and Technology, 61: 2223–2230.

Allen

D.H.

Harris

C.E.

Groves

S.E.

(1987a). A Thermomechanical Constitutive Theory for Elastic Composites with Distributed Damage–I. Theoretical Development, International Journal of Solids and Structures, 23: 1301–1318.

Allen

D.H.

Harris

C.E.

Groves

S.E.

(1987b). A Thermomechanical Constitutive Theory for Elastic Composites with Distributed Damage–II. Application to Matrix Cracking in Laminated Composites, International Journal of Solids and Structures, 23: 1319–1338.

Allen

D.H.

Yoon

(1998). Homogenization Techniques for Thermoviscoelastic Solids Containing Cracks, International Journal of Solids and Structures, 35: 4035–4053.

Bakhvalov

N.S.

Panasenko

G.P.

(1984). Homogenization in Periodic Media. Mathematical Problems of the Mechanics of Composite Materials. Moscow, Nauka.

Bensoussan

Lions

J.L.

Papanicolaou

(1978). Asymptotic Analysis for Periodic Structures. Studies in Mathematics and its Applications, 5, Amsterdam, North-Holland Pub. Co.

Camacho

G.T.

Ortiz

(1996). Computational Modelling of Impact Damage in Brittle Materials, International Journal of Solids and Structures, 33: 2899–2938.

Eshelby

J.D.

(1957). The Determination of the Elastic Field of an Ellipsoidal Inclusion, and Related Problems, Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences (1934–1990), 241: 376–396.

10.

Feyel

(1999). Multiscale FE² Elastoviscoplastic Analysis of Composite Structures, Computational Materials Science, 16: 344–354.

11.

Feyel

Chaboche

J.-L.

(2000). FE² Multiscale Approach for Modelling the Elastoviscoplastic Behaviour of Long Fibre SiC/Ti Composite Materials, Computer Methods in Applied Mechanics and Engineering, 183: 309–330.

12.

Fish

Shek

(2000). Multiscale Analysis of Composite Materials and Structures, Composites Science and Technology, 60: 2547–2556.

13.

Fish

Shek

Pandheeradi

Shephard

M.S.

(1997). Computational Plasticity for Composite Structures Based on Mathematical Homogenization: Theory and Practice, Computer Methods in Applied Mechanics and Engineering, 148: 53–73.

14.

Hashin

(1964). Theory of Mechanical Behavior of Heterogeneous Media, Applied Mechanics Review, 17: 1–9.

15.

Hashin

Shtrikman

(1963). A Variational Approach to the Theory of the Elastic Behaviour of Multiphase Materials, Journal of the Mechanics and Physics of Solids, 11: 127–140.

16.

Hill

(1963). Elastic Properties of Reinforced Solids: Some Theoretical Principles, Journal of the Mechanics and Physics of Solids, 11: 357–372.

17.

Hill

(1965a). Continuum Micro-mechanics of Elastoplastic Polycrystals, Journal of the Mechanics and Physics of Solids, 13: 89–101.

18.

Hill

(1965b). A Self-consistent Mechanics of Composite Materials, Journal of the Mechanics and Physics of Solids, 13: 213–222.

19.

Hill

(1967). The Essential Structure of Constitutive Laws for Metal Composites and Polycrystals, Journal of the Mechanics and Physics of Solids, 15: 79–95.

20.

Kouznetsova

V.G.

(2002). Computational Homogenization for the Multi-scale Analysis of Multi-phase Materials, PhD Dissertation , Technische Universiteit Eindhoven.

21.

Kouznetsova

V.G.

Brekelmans

W.A.M.

Baaijens

F.P.T.

(2001). An Approach to Micro-macro Modeling of Heterogeneous Materials, Computational Mechanics, 27: 37–48.

22.

Kouznetsova

V.G.

Geers

M.G.D.

Brekelmans

W.A.M.

(2004). Multi-scale Second-order Computational Homogenization of Multi-phase Materials: A Nested Finite Element Solution Strategy, Computer Methods in Applied Mechanics and Engineering, 193: 5525–5550.

23.

Nemat-Nasser

Hori

(1999). Micromechanics: Overall Properties of Heterogeneous Materials, 2nd edn, Amsterdam, North-Holland.

24.

Oskay

Fish

(2007). Eigendeformation-based Reduced Order Homogenization for Failure Analysis of Heterogeneous Materials, Computer Methods in Applied Mechanics and Engineering, 196: 1216–1243.

25.

Pandolfi

Ortiz

(2002). An Efficient Adaptive Procedure for Three-dimensional Fragmentation Simulations, Engineering with Computers, 18: 148–159.

26.

Phillips

M.L.

Yoon

Allen

D.H.

(1999). A Computational Model for Predicting Damage Evolution in Laminated Composite Plates, Journal of Engineering Materials and Technology, 121: 436–444.

27.

Sanchez-Palencia

(1980). Non-homogeneous Media and Vibration Theory, Berlin, Springer-Verlag.

28.

Souza

F.V.

(2009). Multiscale Modeling of Impact on Heterogeneous Viscoelastic Solids with Evolving Microcracks, PhD Dissertation , University of Nebraska–Lincoln.

29.

Souza

F.V.

Allen

D.H.

(2010a). Multiscale Modeling of Impact on Heterogeneous Viscoelastic Solids Containing Evolving Microcracks, International Journal for Numerical Methods in Engineering, 82: 464–504.

30.

Souza

F.V.

Allen

D.H.

(2010b). Modeling Failure of Heterogeneous Viscoelastic Solids Under Dynamic/Impact Loading due to Multiple Evolving Cracks Using a Two-way Coupled Multiscale Model, Mechanics of Time-Dependent Materials, 14: 125–151.

31.

Souza

F.V.

Allen

D.H.

Kim

Y.-R.

(2008). Multiscale Model for Predicting Damage Evolution in Composites due to Impact Loading, Composites Science and Technology, 68: 2624–2634.

32.

Tvergaard

(1990). Effect of Fibre Debonding in a Whisker-reinforced Metal, Materials Science and Engineering: A, 125: 203–213.

33.

Zhou

Molinari

J.F.

(2004). Dynamic Crack Propagation with Cohesive Elements: A Methodology to Address Mesh Dependency, International Journal for Numerical Methods in Engineering, 59: 1–24.