Anisotropic Effective Moduli of Microcrack Damaged Media 

by George A. Gazonas 

ARL-SR-191 January 2010 

Approved for public release; distribution is unlimited. 
NOTICES 

Disclaimers 

The findings in this report are not to be construed as an official Department of the Army position unless 

so designated by other authorized documents. 

Citation of manufacturer’s or trade names does not constitute an official endorsement or approval of the 

use thereof. 

Destroy this report when it is no longer needed. Do not return it to the originator. 
Army Research Laboratory 

Aberdeen Proving Ground, MD 21005-5069 

ARL-SR-191 January 2010 

Anisotropic Effective Moduli of Microcrack Damaged Media 

George A. Gazonas 

Weapons and Materials Research Directorate, ARL 

Approved for public release; distribution is unlimited. 
REPORT DOCUMENTATION PAGE Form Approved 

OMB No. 0704-0188

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering 

and maintaining the data needed, and completing and reviewing the collection information. Send comments regarding this burden estimate or any other aspect of this collection of information, 

including suggestions for reducing the burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson 

Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to 

comply with a collection of information if it does not display a currently valid OMB control number. 

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS. 

1. REPORT DATE (DD-MM-YYYY)

January 2010 

2. REPORT TYPE 

Final 

3. DATES COVERED (From - To) 

March 2006–June 2009 

4. TITLE AND SUBTITLE 

Anisotropic Effective Moduli of Microcrack Damaged Media 

5a. CONTRACT NUMBER 

5b. GRANT NUMBER 

5c. PROGRAM ELEMENT NUMBER 

6. AUTHOR(S) 

George A. Gazonas 

5d. PROJECT NUMBER 

AH84 

5e. TASK NUMBER 

5f. WORK UNIT NUMBER 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 

U.S. Army Research Laboratory 

RDRL-WMM-D 

Aberdeen Proving Ground, MD 21005-5069 

8. PERFORMING ORGANIZATION 

 REPORT NUMBER 

ARL-SR-191 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR’S ACRONYM(S) 

11. SPONSOR/MONITOR'S REPORT 

 NUMBER(S) 

12. DISTRIBUTION/AVAILABILITY STATEMENT 

Approved for public release; distribution is unlimited. 

13. SUPPLEMENTARY NOTES 

14. ABSTRACT 

This report combines four recent papers related to using the generalized self-consistent method (GSCM) for determining the 

homogenized constitutive response of microcracked media for use in the development of multiscale constitutive models. In 

“The effect of crack face contact on the anisotropic effective moduli of microcrack damaged media,” the GSCM is used in 

conjunction with a finite-element method to determine the anisotropic effective moduli of a medium containing damage 

consisting of microcracks with an arbitrary degree of alignment. The moduli of the medium subjected to tension, compression, 

and an initially stress-free state are evaluated and shown to be significantly different, affecting the wave speed (illustrated using 

slowness surfaces) in the damaged medium. In “An effective medium model for elastic waves in microcrack damaged media,” 

direct numerical simulations of waves traveling in microcrack damaged media are compared to results using a homogenized 

effective medium calculation. In “Anisotropic effective moduli of microcracked materials under antiplane loading,” the 

anisotropic effective moduli of a cracked solid subjected to antiplane shear deformation are analytically determined. Finally, in 

“On the effective electroelastic properties of microcracked generally anisotropic solids,” concise expressions are derived for the 

effective electroelastic properties of a piezoelectric solid containing insulating, permeable, or conducting microcracks. 

15. SUBJECT TERMS

generalized self-consistent method, multiscale modeling, homogenization, piezoelectric media, finite elements 

16. SECURITY CLASSIFICATION OF: 

17. LIMITATION 

OF ABSTRACT 

UU 

18. NUMBER 

OF PAGES 

68 

19a. NAME OF RESPONSIBLE PERSON

George A. Gazonas 

a. REPORT 

Unclassified 

b. ABSTRACT 

Unclassified 

c. THIS PAGE 

Unclassified 

19b. TELEPHONE NUMBER (Include area code) 

410-306-0863 

Standard Form 298 (Rev. 8/98) 

Prescribed by ANSI Std. Z39.18 
Preface 

This report contains reprints of four papers related to use of the generalized self-consistent 

method (GSCM) for determining the homogenized constitutive response of microcracked media 

for the development of multiscale constitutive models. In “The Effect of Crack Face Contact on 

the Anisotropic Effective Moduli of Microcrack Damaged Media,”1 the GSCM is used in 

conjunction with a finite-element method to determine the anisotropic effective moduli of a 

medium containing damage consisting of microcracks with an arbitrary degree of alignment. 

The moduli of the medium subjected to tension, compression, and an initially stress-free state are 

evaluated and shown to be significantly different, affecting the wave speed (illustrated using 

slowness surfaces) in the damaged medium. In “An Effective Medium Model for Elastic Waves 

in Microcrack Damaged Media,”2 direct numerical simulations of waves traveling in microcrack 

damaged media are compared to results using a homogenized effective medium calculation. For 

incident waves with large wavelengths (1/ka), where k is the wavenumber and a is the half-crack 

length, the scattered elastic energy approaches 0 and the wave does not “see” the obstacle; for 

crack systems modeled using the finite-element approach, this occurs for 1/ka > 60 for media in 

tension and 1/ka > 10 for those in compression. In “Anisotropic Effective Moduli of 

Microcracked Materials Under Antiplane Loading,”3 the anisotropic effective moduli of a 

cracked solid subjected to antiplane shear deformation are analytically determined. When the 

undamaged solid is isotropic, the GSCM can be realized exactly; however, when the undamaged 

solid is anisotropic, coupled nonlinear equations for the unknown effective moduli are 

determined through numerical iteration. Finally, in “On the Effective Electroelastic Properties of 

Microcracked Generally Anisotropic Solids,”4 concise expressions are derived for the effective 

electroelastic properties of a piezoelectric solid containing insulating, permeable, or conducting 

microcracks. 

 1Su, D.; Santare, M. H.; Gazonas, G. A. The Effect of Crack Face Contact on the Anisotropic Effective Moduli of 

Microcrack Damaged Media. Engineering Fracture Mechanics 2007. 

2Su, D.; Santare, M. H.; Gazonas, G. A. An Effective Medium Model for Elastic Waves in Microcrack Damaged Media. 

Engineering Fracture Mechanics 2008. 

3Wang, X.; Santare, M. H.; Gazonas, G. A. Anisotropic Effective Moduli of Microcracked Materials Under Antiplane 

Loading. Engineering Fracture Mechanics 2009. 

4

Wang, X.; Gazonas, G. A.; Santare, M. H. On the Effective Electroelastic Properties of Microcracked Generally Anisotropic 

Solids. Int. J. Fract. 2009. 
The effect of crack face contact on the anisotropic

effective moduli of microcrack damaged media

Dan Su a, Michael H. Santare a,*, George A. Gazonas b

a Department of Mechanical Engineering and Center for Composite Materials, University of Delaware, Newark, DE 19716, USA

b US Army Research Laboratory, Weapons and Materials Research Directorate, Aberdeen Proving Ground, MD 21005, USA

Received 19 March 2006; received in revised form 16 August 2006; accepted 22 August 2006

Available online 18 October 2006

Abstract

The generalized self-consistent method (GSCM) in conjunction with a computational finite element method is used to

calculate the anisotropic effective moduli of a medium containing damage consisting of microcracks with an arbitrary

degree of alignment. Since cracks respond differently under different external loads, the moduli of the medium subjected

to tension, compression and an initially stress-free state are evaluated and shown to be significantly different, which will

further affect the wave speed inside the damaged media. There are four independent material moduli for a 2-D plane stress

orthotropic medium in tension or compression, and seven independent material moduli for a 2-D plane stress orthotropic

cracked medium, which is initially stress free. When friction exists, it further changes the effective moduli. Numerical methods are used to take into account crack face contact and friction. The wave slowness profiles for microcrack damaged

media are plotted using the predicted effective material moduli.

2006 Elsevier Ltd. All rights reserved.

Keywords: Anisotropic damage; Effective moduli; Generalized self-consistent method; Microcracks

1. Introduction

Microcracking is a major cause of damage in brittle material systems, which reduces the material’s effective

moduli by allowing an increase in the large-scale, average deformation resulting from a given external load. As

the number or size of the microcracks increase, there is a commensurate reduction in the effective stiffness,

until the damage reaches a critical value causing the material to fail. A convenient way to quantify the crack

density in 2-D is the parameter g [1].

g 1⁄4 Mc2

A ; ð1Þ

0013-7944/$ - see front matter 2006 Elsevier Ltd. All rights reserved.

doi:10.1016/j.engfracmech.2006.08.015

* Corresponding author. Tel.: +1 302 831 2246; fax: +1 302 831 3619.

E-mail address: santare@me.udel.edu (M.H. Santare).

Engineering Fracture Mechanics 74 (2007) 1436–1455

www.elsevier.com/locate/engfracmech
where M is the number of cracks per unit area A, and c is the average half crack length. In the general case, an

ensemble of microcracks can be characterized by a distribution function, which includes the average orientation of the cracks and a description of how the crack angles are distributed around the average. This will result

in anisotropic effective moduli, even when the undamaged material behaves isotropically. Since microcracks

respond differently under tension and compression, damage growth and the resulting anisotropic effective

moduli of microcracked media will differ depending upon the magnitude and direction of the boundary loads.

The crack-induced changes in the anisotropic effective moduli will have a profound effect on the overall

structural response to loads. The prediction of these changes will help in the understanding of the anisotropic

damage evolution process in brittle materials such as concrete, rock and ceramics. Indeed, continuum damage

finite element modeling of the dynamic response of such brittle materials often requires, at the very least, an

orthotropic description of damage in terms of three orthogonal arrays of interacting cracks [2]. The anisotropic response, is in part, a load-induced effect [3], but can also be attributed to a pre-existing preferred orientation of microcracks. In addition, the mechanical properties of microcracked media have a significant effect

on wave propagation even in cases where the wavelengths are large enough to ignore wave reflection from the

crack faces (see e.g. [4–6]).

Sayers et al. [7,8] used the effective moduli of the damaged media with perfectly aligned or randomly oriented cracks to study the elastic wave propagation. Schubnel and Gue ́guen [9] and Gue ́guen and Schubnel [10]

calculated the elastic wave velocities and permeabilities in microcracked rocks. However, because the interactions of arbitrary orientations of cracks were not fully characterized, these authors only focused on the problems of aligned, perfectly random or specially distributed cracks, and did not consider the difference between

the tensile and compressive behavior of microcracks. Therefore, further work is needed to understand the

complex relationships among microcrack orientation distributions, external loads and the resulting effective

moduli.

Many researchers have proposed methods to calculate the effective moduli for solids with inhomogeneities,

voids or cracks. Without attempting to review the vast literature on the subject, we will focus on the studies

that had the largest direct impact on the current work; the interested reader is referred to a comprehensive

review of the literature provided by Kachanov [11,12]. The first solution that considered the elastic interaction

of elliptically shaped cracks was by Budiansky and O’Connell [13]. In that landmark paper, they used the selfconsistent method (SCM) to calculate the 3-D elastic moduli of a cracked solid, where the distribution of

microcracks was assumed to be isotropic and homogeneous. In the SCM, a crack is embedded directly into

an effective medium and the energy associated with opening the crack under a given load is calculated. By

equating the strain energy associated with the opening of a distribution of such cracks, with the difference

in strain energy between the undamaged and the effective medium, the effective moduli for the cracked medium can be calculated.

Gottesman et al. [14] used variational principles to calculate the effective elastic moduli for materials with

both interacting and non-interacting parallel cracks using the SCM. These two solutions provided upper

and lower bounds for the effective moduli of a medium with aligned cracks, resulting in an orthotropic effective

medium as opposed to the randomly oriented cracks, which results in an isotropic effective medium. Horii and

Nemat-Nasser [3] extended the Budiansky and O’Connell [13] solution, to account for crack closure and crack

face friction effects. They considered load cases in which some cracks opened and some closed, and showed that

when crack closure and friction are considered, the overall material response depends on load history.

Kachanov [11] provided a comprehensive review of earlier results, and compared them to analytical solutions he derived for non-interacting cracks. He noted that there are several drawbacks in the use of the SCM

for determining effective moduli of cracked solids: (1) Because this method inserts the cracks directly into the

effective media, the SCM over estimates the interaction between the cracks. (2) The SCM has a limiting value

for crack density for randomly distributed cracks: it gives zero stiffness when the 2-D crack density g is larger

than 1/p.

In order to overcome these difficulties, several researchers have developed different methods for determining

the effective moduli of microcracked media. Benveniste [15] extended the effective field or Mori–Tanaka

method (MTM) [16] to calculate the effective moduli of 2-D cracked media. This MTM has some similarities

to the SCM, but the energy is calculated based on a model of a crack in an undamaged medium subjected to

an effective far-field stress or strain. The MTM predicts higher moduli than the SCM and does not predict a

D. Su et al. / Engineering Fracture Mechanics 74 (2007) 1436–1455 1437
limiting value for the crack density. Hashin [17] developed a differential scheme method (DSM) for elastic

properties of a cracked material by following a limiting process on the DSM for porous materials. Kachanov

[11] points out that the DSM is an incremental version of the SCM, which does not have a limiting crack density value, but still under-predicts the effective moduli, although not as severely as the SCM. At the same time,

he notes that the DSM is path dependent, since it depends on the order in which the cracks are added to the

formulation.

Aboudi and Benveniste [1] used a generalized self-consistent method (GSCM) to evaluate the 2-D effective

moduli of a material with randomly oriented cracks. This method assumes that each crack resides in a region

of undamaged material, which in turn is surrounded by the effective medium. The general solution for a crack

inside a circular inclusion was used as the basis for this work (see Fig. 1). By using this assumption, the GSCM

gives higher effective moduli than the SCM and gives non-zero stiffness at crack density g = 1/p. However,

because of the use of a circular inclusion, the crack tip touches the inclusion-effective medium interface at

g = 1/p. Therefore, this method can only be used for crack densities less than 1/p and becomes increasingly

less accurate as the crack tip approaches the interface g ! 1/p. Huang et al. [18] provided a crack-elliptical

matrix-composite model for the GSCM; in this improved GSCM, an elliptical inclusion is used instead of a

circular inclusion, which extends the crack density limit seen in Aboudi and Benveniste [1] to a higher limiting

value. The solution, however, is still only presented for a randomly oriented distribution of cracks.

The GSCM accounts for the interactions between cracks in an average sense under the assumption that

each crack is surrounded by a region of intact material. Therefore, in physical situations when cracks are close

together, the strong interactions between them cannot be accounted for within the framework of the GSCM.

For example, Kachanov [11,12] shows that in some pre-arranged crack configurations, the crack density could

become increasingly large, while having very little impact on the effective moduli. For another configuration

with a very low crack density, the effective modulus can approach zero. Kachanov [11,12] points out that the

limitations of the GSCM stem from the fact that this method does not permit overlap of the inclusions containing cracks. This does not preclude modeling of randomly located cracks, so long as the inclusions do not

overlap, yet the maximum crack density is limited to g ’ 0.3. In this statement, Kachanov refers to Aboudi

and Benveniste [1], where the authors use a circular inclusion, which geometrically limits the crack density

E*, μ*

E,μ

2a

Fig. 1. GSCM model from Aboudi and Benveniste [1].

E1

*, E2*, μ12*, K*

θ

E, μ

x

Fig. 2. Model of the crack problem from Santare et al. [19].

1438 D. Su et al. / Engineering Fracture Mechanics 74 (2007) 1436–1455
to g = 1/p as discussed above. At this density, the crack tips touch the inclusion-matrix interface. In Huang

et al. [18] and Santare et al. [19], the authors use an elliptical rather than a circular inclusion. With the elliptical

inclusion this geometric limitation is eliminated and any g value can be used in the calculations. However the

assumption of no overlapping between the inclusions is still present, so the model cannot account for cracks

which nearly touch. Because of this, there is still a practical limitation on the value for g.

In cases where there is strong interaction among multiple cracks, the GSCM is not useful and one would

need other techniques to calculate the effective moduli. Such techniques have been developed. For example,

Greengard and Helsing [20] employ the fast multipole method (FMM) to evaluate the stress field of the elastic

media with different inclusions; Helsing [21] applies FMM in calculating the interactions among large numbers

of cracks; Wang and Chau [22] use a boundary element approach to study the effective stress intensity factor

and the interaction between cracks and holes; Pan [23] derives the single domain BEM formulation to study

the effect of different cracks in different size of domains; Wang and Feng [24] discuss the interaction between

multiple rows of periodical cracks; Dong and Lee [25] apply the boundary integral equation method to evaluate the interaction between the doubly periodic array of cracks and their effect on the mechanical properties.

These methods may be used to account for strong crack interactions resulting from high crack densities and

presumably, from randomly located cracks. But the example solutions in the referenced literature still do not

account for randomly located cracks. In the current paper, we show the results for crack density up to g = 0.4,

which as described in Kachanov [11] ‘‘can be considered as quite high’’ for a 2-D problem, but seems a reasonable upper limit in light of the considerations discussed above.

In most of the above studies, the authors consider either randomly oriented cracks or aligned cracks. For

randomly oriented cracks, the effective medium will be isotropic, since there is no preferred crack orientation.

In the case of aligned cracks, the effective moduli will clearly be anisotropic, since the material behaves differently when the cracks are open than when they are not. In real applications, however, the orientation distribution of the cracks will be neither random nor aligned. For example, the cracks shown in Fig. 3 have an

average orientation parallel to the x-axis and are distributed within ±h0 relative to the x-axis. In this case,

crack opening will influence the moduli in all directions, but not equally. Santare et al. [19] extended Aboudi

and Benveniste’s [1] GSCM to study the 2-D anisotropic case where the cracks have a prescribed orientation

distribution /(h) (Fig. 3). They used the solution for a crack, inside an elliptical inclusion surrounded by an

anisotropic effective medium (Fig. 2). Feltman and Santare [26] extended the Santare et al. [19] work to study

the problem of arbitrarily oriented cracks in a material that is originally anisotropic before damage occurs.

In recent years, the research in this field has continued. Zheng and Du [27] proposed the effective self-consistent method (ESCM) and its simplified explicit version, the interaction direct derivative method (IDDM).

These methods are based on an effective stress/strain field as in the MTM mentioned earlier. Feng et al. [28]

proposed a quasi-micromechanical method (QMM) to calculate the overall constitutive relations for brittle

materials with interacting and growing microcracks. They combined phenomenological observations with

micromechanical models to realize some of the benefits from both approaches. They built the model from

a micromechanical analysis while using the orientation domain of microcrack growth (DMG) from continuum

damage mechanics together with the crack density parameter, to characterize the microcrack damage. In

addition to these quasi-analytical studies, several authors, (see for e.g., [11,29–32]), have provided numerical

solutions to determine the effective moduli in 2-D cracked media.

It is clear that a microcracked medium will behave differently in tension and compression, hence the

effective anisotropic moduli of such a medium should reflect this load-induced anisotropy. Horii and

Nemat-Nasser [3] discuss the response of a medium with randomly distributed cracks under overall

Fig. 3. Cracks distributed within ±h0. h is the angle between the crack and x-axis (label the x-axis in Fig. 2).

D. Su et al. / Engineering Fracture Mechanics 74 (2007) 1436–1455 1439
compressive forces using the SCM. In this solution, some of the pre-existing cracks close, some open, which

lead to the anisotropic response. They show the effective moduli under several different loading conditions.

However, the results are limited to randomly distributed cracks and the limitations of the SCM, such as under

estimating the moduli and predicting zero moduli, are still present in the results.

The problem of crack face overlap has also been recently addressed by Wang and Sun [33] in their work

using boundary element methods. In the present work, we extend the GSCM developed in Santare et al.

[19] and use computational finite element methods to calculate the anisotropic effective moduli for 2-D

cracked media in plane stress, by considering the different behaviors of cracks in tension and compression.

2. Generalized self-consistent method for anisotropic effective media

To begin the discussion, consider the application of the GSCM to the determination of the effective moduli

in a microcrack damaged solid. For a solid body containing discontinuities or voids, the total strain will be the

sum of the deformation of the material itself plus the changes in the void volume. For a specific, external load,

the average deformation of the material is related to the average stress and therefore can easily be found

through the undamaged material’s constitutive relations. To determine the change in void volume, one needs

to know the initial shape and distribution of the voids. In the current model, the voids are assumed to be

homogeneously distributed microcracks and the strain associated with them is therefore related to the crack

opening displacements. Assuming homogeneous loading and elastic material response, the strain relationship

described above can be written in terms of strain energies ([1]; see also Walsh [34] who derives a similar expression for isotropic solids using the reciprocal theorem),

1

2

S

ijklr0

ijr0ij 1⁄4 1

2

Sijklr0

klr0ij þ

1

2V

XM

k1⁄41

Z

Ck

1⁄2uit

0

i dCk : ð2Þ

In this expression, S

ijk‘ is the effective compliance of the damaged material, and Sijk‘ is the compliance of the

undamaged material, r0

ij is the homogeneous, applied Cauchy stress field, and

t

0

i 1⁄4 r0ijnj ð3Þ

is the traction that would be present along the crack faces if the material were not damaged. The term in

brackets [ui] is the crack opening displacement so that the integral in the expression, taken over the crack surface Ck of each of the M cracks in the volume V, gives the total work associated with crack opening. This

integral is evaluated to take into account the crack distribution, by considering all the existing microcrack orientations, appropriately weighted. To determine ‘‘N’’ distinct effective moduli, ‘‘N’’ different homogeneous

stress fields r0

ij are applied to Eq. (2), and the different forms of the equation are solved simultaneously to

determine the ‘‘N’’ effective elastic constants (see for example [19]).

3. Crack face contact

The crack-opening energy integrals used in the previous studies are based on linear elastic fracture mechanics (LEFM), which does not explicitly account for crack face contact. Therefore, the strain energy associated

with compressive loading is the same as that associated with the tensile loading due to the fact that in LEFM

the crack faces overlap in compression. In reality, the crack faces will open or close depending on the applied

loads, resulting in very different strain energies in tension and compression. Consider a two-dimensional plane

stress orthotropic elastic material with cracks, under bi-axial loads P and P* (Fig. 4a) and a small load increment Dq, which is not necessarily a bi-axial load.

We will discuss the following three different cases as shown in Fig. 4b. The r1 and r2 in Fig. 4b represent the

principle stresses in the 2-D plane.

Case 1: P > 0, P* > 0 and jPjjDqj, jP*jjDqj

In this case, the pre-existing bi-axial tensile loads P and P* are much larger than the load increment Dq.

Because of this, all cracks are open and there is no crack face contact. The results in Santare et al. [19] are

valid in this case.

1440 D. Su et al. / Engineering Fracture Mechanics 74 (2007) 1436–1455
Case 2: P < 0,P* < 0 and jPjjDqj, jP*jjDqj

In this case, the pre-existing bi-axial loads P and P* are compressive and much larger than the load increment

Dq. Therefore, all cracks remain closed in mode I (normal). However the cracks may still open in mode II (shear).

Case 3: P = P* = 0

The medium is initially stress free before the load increment is applied. In this case, the cracks may either

open or close in both mode I and II depending on the direction of the load increments Dq.

Fig. 4c shows schematically the cracks’ response to these three cases under a pure shear load increment as

an example of the Dq. In Case 1, all the cracks open; in Case 2, all the cracks close in mode I; and in Case 3,

some cracks open and some close.

For Cases 1 and 2, four independent material constants are required to fully characterize S

ijk‘ in the plane

(plane stress); thus four independent load cases are needed as shown schematically in Fig. 5 (for Case 1) and

Fig. 6 (for Case 2). Taking each load case separately and assuming that the undamaged material is isotropic,

the following equations result [19]:

1

2E

1

p2

1A 1⁄4 1

2E p2

1A þ

1

2

XM

k1⁄41

Z

Ck

1⁄2ui

p1 t

p1

i dCk ; ð4Þ

1

2E

2

p2

2A 1⁄4 1

2E p2

2A þ

1

2

XM

k1⁄41

Z

Ck

1⁄2ui

p2 t

p2

i dCk ; ð5Þ

P* P*

P 

P 

σ1

σ 2

Case 1 

Tension 

Case2 

Compression

Case3 

Initially stress free 

ο 45

Case 2

P 

P 

P* P*

P* P*

P 

Case 1

Case 3 

P 

Fig. 4. (a) A cracked medium under bi-axial load. (b) Three different load cases. (c) Three cases with pure shear load increment.

D. Su et al. / Engineering Fracture Mechanics 74 (2007) 1436–1455 1441
1

2lq2

A 1⁄4 1

2l

q2

A þ

1

2

XM

k1⁄41

Z

Ck

1⁄2ui

q

t

q

i dCk ; ð6Þ

1

2Kh2

A 1⁄4 1

2K

h2

A þ

1

2

XM

k1⁄41

Z

Ck

1⁄2ui

h

t

h

i dCk : ð7Þ

In these equations, p1, p2, q and h are the homogeneous load increments as shown in Figs. 5 and 6. E

1; E2; l

and K* are the effective moduli of the damaged material. E

1 and E2 are the elastic moduli in the principal directions, l*, is the shear modulus, and K* is the effective 2-D plane stress bulk modulus. These moduli are related

to Poisson’s ratio through the equation [19],

K1⁄4 E

1E2

E

1 þ E2ð1  2v

1Þ

; ð8Þ

where v

1 is one of the effective Poisson’s ratios of the damaged material. E, l and K are the undamaged material moduli where K is analogous to K* for the case of an isotropic material.

In Case 3, the cracked medium responds differently under tensile and compressive loads, hence there will be

different Young’s moduli, and bulk moduli under tension and compression (Fig. 7). We assume the crack distribution is symmetric with respect to the x-axis. Therefore, when we shear the medium relative to this axis,

half of the cracks will open and half will close in mode I (as depicted in Fig. 4c), so there is only one effective

shear modulus l* in Case 3. Seven independent material constants are thus required to describe this 2-D orthotropic plane stress model

P 

p1

p2

P 

q 

P 

h

P 

P* P*

P* P*

Fig. 5. Four load cases for medium in tension.

p1

P 

p2

P 

q 

P 

h

P 

P* P* 

P* P* 

Fig. 6. Four load cases for medium in compression.

1442 D. Su et al. / Engineering Fracture Mechanics 74 (2007) 1436–1455
1

2E

1T

p2

1TA 1⁄4 1

2E p2

1TA þ

1

2

XM

k1⁄41

Z

Ck

1⁄2ui

p1T t

p1T

i dCk ; ð9Þ

1

2E

1C

p2

1CA 1⁄4 1

2E p2

1CA þ

1

2

XM

k1⁄41

Z

Ck

1⁄2ui

p1C t

p1C

i dCk ; ð10Þ

1

2E

2T

p2

2TA 1⁄4 1

2E p2

2TA þ

1

2

XM

k1⁄41

Z

Ck

1⁄2ui

p2T t

p2T

i dCk ; ð11Þ

1

2E

2C

p2

2CA 1⁄4 1

2E p2

2CA þ

1

2

XM

k1⁄41

Z

Ck

1⁄2ui

p2C t

p2C

i dCk ; ð12Þ

1

2lq2

A 1⁄4 1

2l

q2

A þ

1

2

XM

k1⁄41

Z

Ck

1⁄2ui

q

t

q

i dCk ; ð13Þ

1

2K

T

h2

TA 1⁄4 1

2K

h2

TA þ

1

2

XM

k1⁄41

Z

Ck

1⁄2ui

hT t

hT

i dCk ; ð14Þ

1

2K

C

h2

CA 1⁄4 1

2K

h2

CA þ

1

2

XM

k1⁄41

Z

Ck

1⁄2ui

hC t

hC

i dCk : ð15Þ

In these equations, p1T, p1C, p2T, p2C, q, hT, hC represent the homogeneous load increments shown in Fig. 7.

E

1T and E2T are the effective Young’s moduli under tensile load, and E1C; E2C are the effective Young’s moduli

under compressive load. K

T; KC are the effective 2-D plane stress bulk moduli under tensile and compressive

loads. By generalizing the results of Taya and Mura [35] and Feltman and Santare [26], they can be written as

K

T 1⁄4 E

1TE2T

E

1T þ E2Tð1  2v

12TÞ

; ð16Þ

K

C 1⁄4 E

1CE2C

E

1C þ E2Cð1  2v

12CÞ

; ð17Þ

where v

12T represents the effective Poisson’s ratios of the damaged material under tensile loads, and v12C is for

compressive loads.

4. Finite element model

The results in Santare et al. [19] are valid only for Case 1, as discussed in the previous section. In that case,

the crack surfaces will never overlap due to the applied bi-axial tensile loads. In order to prevent the unphysical crack face overlap under compression loads in LEFM, and discuss other cases mentioned in previous section, we have developed a finite element model to calculate the energy associated with the crack opening

p1T

p2T

q hT

p1C

p2C

hC 

Fig. 7. Seven load cases for an initially stress-free medium.

D. Su et al. / Engineering Fracture Mechanics 74 (2007) 1436–1455 1443
displacements. Nonlinear contact elements are used along the crack surfaces, so that overlap is eliminated, and

crack face friction can be added easily.

The finite element model was designed to mimic as closely as possible the analytical boundary value problem used in the GSCM. This model is shown schematically in Fig. 8. The use of finite elements introduces several new issues. For example, the analytical solution domain used in GSCM is an infinite medium. To simulate

this using an FE model we need a large enough domain so that the embedded crack does not interact with the

external boundary. However, the larger the model, the more calculation time needed to converge to a set of

effective moduli. A second issue is the interface between the ellipse and the surrounding medium. In the FE

model, the mesh needs to be sufficiently refined to simulate the smooth elliptical interface in the analytical

solution. A comprehensive meshing study was undertaken and the numerical results were compared to known

analytical solutions [36] to resolve these domain size and finite element mesh refinement issues.

The resulting FE model is constructed with 2-D plane stress, isoparametric elements. As in the analytical

solution, the material inside the ellipse is isotropic and homogeneous and the material outside is a homogeneous

orthotropic effective material with a known principal axis orientation. Note that in Fig. 8, the crack shown is

horizontal. In order to calculate the crack opening displacement for cracks at other angles, we simply modify

the externally applied loads as calculated from Mohr’s circle. The external boundary length of the FE model is

L = 500 and the focal length of the ellipse equals the half crack length c = 30. In order to model a specific crack

density g, the crack length is kept constant, and the size of the ellipse is varied. By using Eq. (1), together with

the geometrical relations for an ellipse, A = pab and c

2 = a2  b2

, we can determine a unique a and b, the semimajor and semi-minor axes of the ellipse for a given g. For the results that follow, the FE model has 31,882

elements in the entire 2-D plane stress model, with 4054 elements inside the ellipse when g = 0.4.

The effective moduli for the cracked media are obtained using the FE solutions for the crack energies in the

method described in Section 2 and 3. We use ABAQUS 6.4.4 [37] as the FE solver and the program is written

in the ABAQUS script language Python. For the initially stress-free case (Case 3), we wrote a user subroutine

to define a new 2-D material with seven material constants. On a PC with a Pentium 4 2.0 GHz cpu and 1 GB

RAM, it took from 30 min to 3 h for a single iteration with a specific g and h0, depending on the number of

unknown material properties (4 for Cases 1 and 2 and 7 for Case 3), whether the friction existed and whether

we used the user defined material. There are about 5–15 iterations in a ‘‘typical’’ solution, when we use the

virgin material properties as the initial guess.

To describe the algorithm of the FE script, consider the medium in tension (Case 1) as an example. Let us

assume that the crack distribution function is /(h) and the individual crack angles are distributed in the interval [h0,h0] relative to the x-axis. With this, we can replace the summations in Eqs. (4)–(7) with the following

integrals (adapted from [19]):

1

2E

1

p2

1 1⁄4 1

2E p2

1 þ

1

4

g

c2h0

Z h0

h0

/ðhÞ

Z

Ck

1⁄2ui

p1 t

p1

i dCk dh; ð18Þ

1

2E

2

p2

2 1⁄4 1

2E p2

2 þ

1

4

g

c2h0

Z h0

h0

/ðhÞ

Z

Ck

1⁄2ui

p2 t

p2

i dCk dh; ð19Þ

1

2lq2 1⁄4 12l

q2 þ

1

4

g

c2h0

Z h0

h0

/ðhÞ

Z

Ck

1⁄2ui

q

t

q

i dCk dh; ð20Þ

L 

L 

2c

2a

Crack

Fig. 8. FEM model for single crack inside an elliptical inclusion.

1444 D. Su et al. / Engineering Fracture Mechanics 74 (2007) 1436–1455
1

2Kh2 1⁄4 12K

h2 þ

1

4

g

c2h0

Z h0

h0

/ðhÞ

Z

Ck

1⁄2ui

h

t

h

i dCk dh: ð21Þ

Fig. 9 is a flow chart of the algorithm used for calculating the effective moduli of the medium in tension. The

quantities assumed known, are the crack density, g (hence the elliptical geometry and crack length) the undamaged material moduli E, l and K and the crack orientation distribution, /(h), and range [h0,h0]. An initial

guess is made for the four effective moduli in the surrounding matrix. By applying these and the traction

boundary conditions to the model as shown in Fig. 5, we calculate the crack opening displacement [ui] in each

load case for the full range of crack orientations considered. These displacements allow us to approximate the

double integrals in Eqs. (18)–(21). We can then use these four equations to calculate a new set of the four effective moduli E

1; E2; land K*. These new moduli are then compared to the initial estimates to determine if

another iteration is necessary to converge within a user-defined level of accuracy.

5. Computational results

As discussed before, different crack distribution functions /(h) will result in different orthotropic effective

moduli. In the following solutions, we choose /(h) to be a constant distribution between [h0,h0] as shown in

Fig. 9. Flow chart for calculating effective moduli of the medium in tension.

D. Su et al. / Engineering Fracture Mechanics 74 (2007) 1436–1455 1445
Fig. 10. This is chosen so that when h0 = 90, the cracks are randomly distributed as is the case in most of the

literature. Other distributions could be used simply by inputting the proper function into Eqs. (18)–(21). We

assume that the intact Poisson’s ratio is 1/3 and that the crack tips are placed at the foci of the ellipse in Fig. 8.

In the tables, we list the results for the effective moduli for all three load cases described in Section 3. The

effective moduli for the medium in tension (Case 1) are all within 4% relative error when compared to the

results in Santare et al. [19]. This verifies that the FE model developed approximates the analytical solution

to within an acceptable level of accuracy when all the cracks are open. In Figs. 11 and 12, we compare our

current results with several classic results for the effective moduli of cracked media in 2-D plane stress models.

Fig. 11 is the effective Young’s modulus for the medium with randomly distributed cracks (h0 = 90) and

Fig. 12 shows the effective Young’s modulus in the direction perpendicular to an aligned array of cracks

(h0 = 0).

As discussed above, in a 2-D orthotropic plane stress model, four material moduli (Cases 1 and 2) or seven

material moduli (Case 3) are required to fully define the material compliance. In plane stress, the mechanical

properties in the out-of-plane direction are uncoupled from those in the in-plane directions and are therefore

completely arbitrary.

However, in the literature, many authors predict the effective moduli for the cracked medium using a 2-D

plane strain model. In plane strain, the out-of-plane compliances affect the in-plane moduli. Yet, many

authors do not report the mechanical properties in the out-of-plane direction, and this makes it difficult to

compare our result to theirs. Alternatively, it is possible to use the current procedure, along with an assumed

set of out-of-plane properties, to generate an analogous set of plane strain results.

-θ0 +θ0

φ(θ)

Fig. 10. Constant distribution of cracks.

0.0 0.1 0.2 0.3 0.4

0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

E*/E

η

 Non-interacting Kachanov [11]

 Current result for a medium in tension

 Current result for a medium in compression

 SCM Kachanov [11]

Fig. 11. Normalized effective Young’s modulus for randomly distributed cracks vs. crack density.

1446 D. Su et al. / Engineering Fracture Mechanics 74 (2007) 1436–1455
The data in Tables 1–4 are plotted in a series of contour plots (Figs. 13–15) to facilitate visualization and

comparison. In each contour plot, the horizontal axes are the limits of the orientation distribution, h0 from 0

to 90 and the crack density, g from 0 to 0.4. The vertical axis is the effective modulus normalized by the corresponding undamaged modulus. Figs. 13 and 14 show the moduli for the medium in tension (Case 1) and

compression (Case 2). Comparing Fig. 13 with Fig. 14, we see that the medium is generally stiffer under compressive loads than under tensile loads. This is because the cracks can only experience mode II (shear) opening

under compressive loads and therefore the energy dissipated through crack opening is smaller than that under

tensile loads where modes I and II are both possible. As shown in Fig. 14, the 2-D effective bulk modulus is

unchanged under compressive loads. This is due to the fact that when the medium is loaded in bi-axial compression, all the cracks will be closed (in both normal and shear modes).

Fig. 15 shows how the effective moduli of the initially stress-free medium (Case 3) vary with h0 and g. In this

figure, when h0 = 90, the crack orientations are random and therefore the medium is isotropic. The classical

relationship among isotropic material moduli is shown as

l 1⁄4 E

2ð1 þ vÞ

: ð22Þ

However, the effective properties in Fig. 15 do not satisfy Eq. (22). In an isotropic medium, where the tensile

properties are different from the compressive properties, the following relationship can be easily derived:

l 1⁄4 ECET

ECð1 þ vTÞ þ ETð1 þ vCÞ

; ð23Þ

where EC, ET are Young’s moduli under tensile and compressive loads and vT, vC are Poisson’s ratios in these

two load cases. In this case, the Poisson’s ratios in Eq. (23) can be calculated from the 2-D bulk moduli and

Young’s moduli through the following relations:

KT 1⁄4 ET

2ð1  vTÞ

; ð24Þ

KC 1⁄4 EC

2ð1  vCÞ

; ð25Þ

Eqs. (24) and (25) are the isotropic versions of Eqs. (16) and (17).

0.0 0.1 0.2 0.3 0.4

0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

E2*/E

η

 Non-interacting Kachanov [11]

 Current result for a medium in tension

 Current result for a medium in compression

 SCM Gottesman et al [14]

Fig. 12. Normalized effective Young’s modulus, in the direction perpendicular to the cracks, for aligned cracks vs. crack density.

D. Su et al. / Engineering Fracture Mechanics 74 (2007) 1436–1455 1447
When friction is present on the crack surfaces, the moduli become path and load dependent. This is because

the friction force is a non-conservative force, and therefore the work of friction depends on the sequence of

Table 1

Effective moduli of a medium in tension

h0 E

1=E E2=E l

12=l K*/K h0 E

1=E E2=E l

12=l K*/K

g = 0.1 g = 0.2

0 1.000 0.553 0.796 0.623 0 1.000 0.316 0.633 0.381

10 0.993 0.556 0.795 0.623 10 0.985 0.319 0.632 0.382

20 0.974 0.566 0.795 0.625 20 0.944 0.328 0.630 0.384

30 0.944 0.581 0.794 0.627 30 0.883 0.344 0.628 0.388

40 0.907 0.601 0.794 0.629 40 0.812 0.366 0.626 0.391

50 0.866 0.626 0.793 0.631 50 0.739 0.393 0.624 0.395

60 0.826 0.653 0.793 0.633 60 0.670 0.424 0.624 0.397

70 0.789 0.681 0.793 0.634 70 0.611 0.459 0.624 0.399

80 0.758 0.708 0.794 0.634 80 0.562 0.494 0.625 0.400

90 0.732 0.732 0.794 0.634 90 0.526 0.526 0.626 0.400

g = 0.3 g = 0.4

0 1.000 0.182 0.502 0.228 0 1.000 0.106 0.398 0.137

10 0.976 0.184 0.500 0.229 10 0.965 0.108 0.395 0.138

20 0.910 0.191 0.496 0.231 20 0.870 0.113 0.388 0.139

30 0.816 0.203 0.491 0.234 30 0.740 0.120 0.380 0.140

40 0.711 0.220 0.486 0.237 40 0.611 0.131 0.372 0.140

50 0.611 0.242 0.482 0.239 50 0.484 0.145 0.366 0.139

60 0.523 0.268 0.481 0.240 60 0.387 0.163 0.364 0.138

70 0.452 0.298 0.481 0.240 70 0.316 0.184 0.364 0.137

80 0.399 0.330 0.483 0.241 80 0.266 0.208 0.366 0.136

90 0.361 0.361 0.485 0.241 90 0.234 0.234 0.368 0.136

Table 2

Effective moduli of a medium in compression

h0 E

1=E E2=E l

12=l K*/K h0 E

1=E E2=E l

12=l K*/K

g = 0.1 g = 0.2

0 1.000 1.000 0.797 1.000 0 1.000 1.000 0.639 1.000

10 0.994 0.994 0.804 1.000 10 0.987 0.987 0.650 1.000

20 0.977 0.977 0.823 1.000 20 0.953 0.953 0.680 1.000

30 0.955 0.955 0.850 1.000 30 0.909 0.909 0.724 1.000

40 0.933 0.933 0.878 1.000 40 0.868 0.868 0.772 1.000

50 0.917 0.917 0.900 1.000 50 0.838 0.838 0.810 1.000

60 0.909 0.909 0.911 1.000 60 0.824 0.824 0.830 1.000

70 0.910 0.910 0.911 1.000 70 0.825 0.825 0.828 1.000

80 0.916 0.916 0.902 1.000 80 0.836 0.836 0.813 1.000

90 0.924 0.924 0.890 1.000 90 0.851 0.851 0.793 1.000

g = 0.3 g = 0.4

0 1.000 1.000 0.516 1.000 0 1.000 1.000 0.418 1.000

10 0.980 0.980 0.528 1.000 10 0.973 0.973 0.432 1.000

20 0.929 0.929 0.564 1.000 20 0.903 0.903 0.470 1.000

30 0.864 0.864 0.618 1.000 30 0.819 0.819 0.528 1.000

40 0.805 0.805 0.677 1.000 40 0.745 0.745 0.595 1.000

50 0.764 0.764 0.728 1.000 50 0.693 0.693 0.653 1.000

60 0.744 0.744 0.754 1.000 60 0.670 0.670 0.684 1.000

70 0.745 0.745 0.752 1.000 70 0.671 0.671 0.682 1.000

80 0.761 0.761 0.731 1.000 80 0.690 0.690 0.657 1.000

90 0.782 0.782 0.705 1.000 90 0.715 0.715 0.626 1.000

1448 D. Su et al. / Engineering Fracture Mechanics 74 (2007) 1436–1455
loading. To simplify the current discussion, here we will only consider the initially stress-free medium with the

external loads applied proportionally.

Table 3

Effective moduli of initially stress-free media

h0 E

1T=E E1C=E E2T=E E2C=E l

12=l K

t =K Kc=K h0 E1T=E E1C=E E2T=E E2C=E l

12=l K

t =K Kc=K

g = 0.1 g = 0.2

0 1.000 1.000 0.553 1.000 0.796 0.623 1.000 0 1.000 1.000 0.316 1.000 0.634 0.381 1.000

10 0.993 0.994 0.556 0.994 0.799 0.623 1.000 10 0.986 0.987 0.319 0.987 0.640 0.382 1.000

20 0.975 0.977 0.565 0.977 0.809 0.625 1.000 20 0.947 0.953 0.328 0.953 0.654 0.385 1.000

30 0.946 0.954 0.580 0.954 0.821 0.627 1.000 30 0.889 0.909 0.344 0.909 0.673 0.389 1.000

40 0.910 0.933 0.600 0.933 0.834 0.629 1.000 40 0.821 0.867 0.366 0.867 0.693 0.394 1.000

50 0.871 0.917 0.624 0.917 0.844 0.632 1.000 50 0.749 0.837 0.393 0.837 0.708 0.398 1.000

60 0.831 0.909 0.651 0.909 0.849 0.633 1.000 60 0.681 0.823 0.424 0.823 0.716 0.401 1.000

70 0.793 0.910 0.679 0.910 0.849 0.635 1.000 70 0.620 0.823 0.459 0.823 0.715 0.403 1.000

80 0.760 0.916 0.707 0.916 0.845 0.635 1.000 80 0.570 0.835 0.496 0.835 0.709 0.404 1.000

90 0.733 0.924 0.733 0.924 0.840 0.635 1.000 90 0.530 0.850 0.530 0.850 0.701 0.404 1.000

g = 0.3 g = 0.4

0 1.000 1.000 0.182 1.000 0.505 0.228 1.000 0 1.000 1.000 0.107 1.000 0.402 0.137 1.000

10 0.978 0.980 0.184 0.980 0.511 0.230 1.000 10 0.968 0.973 0.108 0.973 0.408 0.138 1.000

20 0.916 0.928 0.192 0.928 0.526 0.233 1.000 20 0.880 0.902 0.114 0.902 0.423 0.141 1.000

30 0.828 0.863 0.205 0.863 0.547 0.237 1.000 30 0.761 0.816 0.122 0.816 0.443 0.144 1.000

40 0.728 0.803 0.222 0.803 0.569 0.241 1.000 40 0.632 0.739 0.134 0.739 0.463 0.146 1.000

50 0.630 0.760 0.244 0.760 0.586 0.245 1.000 50 0.514 0.686 0.149 0.686 0.478 0.147 1.000

60 0.542 0.740 0.271 0.740 0.595 0.247 1.000 60 0.415 0.661 0.168 0.661 0.486 0.146 1.000

70 0.469 0.741 0.302 0.741 0.594 0.248 1.000 70 0.339 0.662 0.191 0.662 0.485 0.145 1.000

80 0.412 0.756 0.336 0.756 0.587 0.248 1.000 80 0.283 0.680 0.217 0.680 0.479 0.144 1.000

90 0.370 0.777 0.370 0.777 0.578 0.248 1.000 90 0.246 0.706 0.245 0.706 0.470 0.143 1.000

Table 4

Effective moduli of initially stress-free media with friction coefficient equal to 0.2

h0 E

1T=E E1C=E E2T=E E2C=E l

12=l K

t =K Kc=K h0 E1T=E E1C=E E2T=E E2C=E l

12=l K

t =K Kc=K

g = 0.1 g = 0.2

0 1.000 1.000 0.553 1.000 0.796 0.623 1.000 0 1.000 1.000 0.316 1.000 0.634 0.381 1.000

10 0.993 0.994 0.556 1.000 0.802 0.623 1.000 10 0.986 0.987 0.319 1.000 0.644 0.382 1.000

20 0.975 0.978 0.565 0.994 0.814 0.625 1.000 20 0.947 0.955 0.329 0.988 0.661 0.385 1.000

30 0.946 0.958 0.580 0.978 0.828 0.627 1.000 30 0.890 0.915 0.344 0.955 0.683 0.390 1.000

40 0.910 0.940 0.600 0.959 0.841 0.629 1.000 40 0.821 0.880 0.366 0.917 0.703 0.394 1.000

50 0.871 0.928 0.624 0.943 0.850 0.632 1.000 50 0.750 0.857 0.393 0.886 0.717 0.399 1.000

60 0.831 0.924 0.651 0.934 0.855 0.634 1.000 60 0.681 0.850 0.425 0.868 0.725 0.402 1.000

70 0.793 0.928 0.680 0.932 0.855 0.635 1.000 70 0.621 0.856 0.460 0.864 0.725 0.404 1.000

80 0.760 0.935 0.708 0.936 0.851 0.635 1.000 80 0.570 0.869 0.496 0.871 0.720 0.405 1.000

90 0.733 0.942 0.733 0.942 0.846 0.635 1.000 90 0.531 0.883 0.531 0.883 0.711 0.405 1.000

g = 0.3 g = 0.4

0 1.000 1.000 0.182 1.000 0.505 0.229 1.000 0 1.000 1.000 0.107 1.000 0.413 0.138 1.000

10 0.978 0.980 0.184 1.000 0.515 0.230 1.000 10 0.968 0.972 0.109 1.000 0.418 0.139 1.000

20 0.916 0.930 0.192 0.981 0.535 0.233 1.000 20 0.881 0.904 0.114 0.974 0.436 0.141 1.000

30 0.828 0.871 0.205 0.931 0.559 0.238 1.000 30 0.762 0.824 0.123 0.905 0.458 0.144 1.000

40 0.729 0.819 0.222 0.874 0.581 0.242 1.000 40 0.634 0.758 0.135 0.828 0.479 0.147 1.000

50 0.631 0.786 0.245 0.827 0.597 0.245 1.000 50 0.515 0.716 0.150 0.768 0.493 0.147 1.000

60 0.543 0.775 0.272 0.801 0.606 0.248 1.000 60 0.417 0.703 0.169 0.735 0.501 0.147 1.000

70 0.470 0.784 0.303 0.796 0.606 0.249 1.000 70 0.340 0.713 0.192 0.728 0.501 0.146 1.000

80 0.413 0.804 0.337 0.806 0.600 0.249 1.000 80 0.285 0.738 0.219 0.742 0.495 0.145 1.000

90 0.371 0.824 0.371 0.824 0.589 0.248 1.000 90 0.247 0.764 0.247 0.764 0.486 0.144 1.000

D. Su et al. / Engineering Fracture Mechanics 74 (2007) 1436–1455 1449
Fig. 16 shows the comparison of normalized effective Young’s moduli for the initially stress-free media with

and without friction under compressive load. The crack density g = 0.4, and in the friction case, the friction

Fig. 13. Effective moduli of a medium in tension.

Fig. 14. Effective moduli of a medium in compression.

1450 D. Su et al. / Engineering Fracture Mechanics 74 (2007) 1436–1455
coefficient is equal to 0.2. We can see that the friction affects to the Young’s moduli under compressive loads.

When g = 0.4, E

1C and E2C increase about 8% when the friction exists. Also, we see that there is a clear ‘‘flat

zone’’ ðE

2C=E 1⁄4 1Þ for E2C when h0 < 11. This is because the friction prevents the cracks from opening in any

mode under this situation. The maximum extent of the ‘‘flat zone’’ is determined by the friction coefficient.

Fig. 15. Effective moduli of an initially stress-free medium.

D. Su et al. / Engineering Fracture Mechanics 74 (2007) 1436–1455 1451
0 10 20 30 40 50 60 70 80 90

0.60

0.65

0.70

0.75

0.80

0.85

0.90

0.95

1.00

Normalized Effective Young's moduli

θ0

 E1*/E and E2*/E without friction

 E1*/E with friction 

 E2*/E with friction 

Fig. 16. Normalized E

1C; E2C of initially stress-free media with and without friction. g = 0.4 and the friction coefficient equal to 0.2.

Fig. 17. Normalized slowness of cracked media under tension and compression. Crack density g = 0.4. Outside lines: quasi-transverse

wave; inside lines: quasi-longitudinal wave. All the lines are normalized by the slowness of transverse wave in virgin material.

1452 D. Su et al. / Engineering Fracture Mechanics 74 (2007) 1436–1455
6. Wave propagation in microcrack damaged media

When applying the predicted effective moduli to the wave equation, we can predict the elastic wave speed

within the cracked media. To find out the wave speed of a plane harmonic wave in the L direction, we need to

solve the following eigenvalue equation [38]:

detfAij  dijqv2

g 1⁄4 0; ð26Þ

in which, A = LTCL is the acoustic tensor, C is the elastic tensor and L is a 3 · 2 direction cosine matrix.

1⁄2L 1⁄4

l1 0

0 l2

l2 l1

2

6

4

3

7

5: ð27Þ

The slowness is the inverse of wave speed

SðiÞ

ðlÞ 1⁄4 l

vðiÞ : ð28Þ

Fig. 17 shows the normalized slowness (inverse of the wave speed) in the media with cracks in tension and

compression. The outside lines depict curves of maximum slowness (quasi-transverse wave) and the inside lines

depict curves of minimum slowness (quasi-longitudinal wave). The figures on the first line are the media under

tensile loads (Case 1 in Fig. 4b) and the figures on the second line are the media under compressive loads (Case

2 in Fig. 4b). Comparing the figures in the same column, we can see that for identical crack distributions and

densities, the wave speed (slowness) profiles are very different under different loads.

7. Conclusions

In this study, a finite element model was used to solve the analytical crack-inclusion anisotropic boundary

value problem used in the GSCM. The analytical results of Santare et al. [19] were used to validate the finite

element model within a 4% relative error (Case 1), which confirms the accuracy of the model. Using the same

basic FE model, contact elements were introduced to eliminate the crack surface overlap, which occurs under

compression in linear elastic fracture mechanics. With this model, we calculated effective moduli for a damaged medium under tensile, compressive, and initially stress-free conditions. Base on the effective moduli,

we also predict the elastic wave slowness within the 2-D cracked media.

Generally speaking, the results show that the existence of microcracks decreases the effective stiffness of the

medium and the effective moduli further decrease as the microcrack density is increased. The crack orientation

distribution also has a significant effect on the anisotropy of the effective moduli which can vary from highly

anisotropic when the cracks are unidirectional (h0 = 0) to isotropic when they are randomly distributed

(h0 = 90). Furthermore, we see that the damaged medium responds quite differently depending on whether

the external loads are tension, compression or initially zero. And we see that the medium is much stiffer when

the overall external loads are compressive.

We have described how, in an initially stress-free medium under 2-D plane stress loading, there are seven

independent orthotropic effective material properties for the microcracked medium (Case 3); for randomly distributed microcracks, the number of independent moduli are reduced to four, and satisfy Eq. (23). We have

also shown how friction can be introduced into the model, and that in general, friction makes the medium

stiffer particularly under compression loads. From the literature [3] we know that the moduli become path

and load dependent when crack face friction is included in the analysis, and to adequately account for these

effects, further study is required.

Finally, we plot wave slowness profiles for microcrack damaged media and show that for media with identical microcrack densities and distributions the wave speeds are very different under tensile and compressive

loads, because of the different effective moduli for these two load cases.

D. Su et al. / Engineering Fracture Mechanics 74 (2007) 1436–1455 1453
Acknowledgement

This research was supported by the United States Army Research Laboratory through the Composite

Materials Technology cooperative in agreement with the Center for Composite Materials at the University

of Delaware.

References

[1] Aboudi J, Benveniste Y. The effective moduli of cracked bodies in plane deformations. Engng Fract Mech 1987;26(2):171–84.

[2] Fahrenthold EP. A continuum damage model for fracture of brittle solids under dynamic loading. Trans ASME 1991;58:

904–9.

[3] Horii H, Nemat-Nasser S. Overall moduli of solid with microcracks: load-induced anisotropy. J Mech Phys Solids 1983;31(2):

155–71.

[4] Bayuk IO, Kalinin VA. Elastic anisotropy of rocks 2. An oriented system of cracks of arbitrary concentration. Phys Solid Earth

1995;31(3):194–200.

[5] Monsen K. Acoustic velocity in fractured rock. J Geophys Res 2001;106(B7):13261–7.

[6] Liu E, Hudson JA, Pointer T. Equivalent medium representation of fractured rock. J Geophys Res 2000;105(B2):2981–3000.

[7] Sayers CM, van Munster JG. Microcrack-induced seismic anisotropy of sedimentary rocks. J Geophys Res 1991;96(B10):529–33.

[8] Sayers CM, Kachanov M. Microcrack-induced elastic wave anisotropy of brittle rocks. J Geophys Res 1995;100(B3):4149–56.

[9] Schubnel A, Gue ́guen Y. Dispersion and anisotropy of elastic waves in cracked rocks. J Geophys Res 2003;108(B2):2101.

[10] Gue ́guen Y, Schubnel A. Elastic wave velocities and permeability of cracked rocks. Tectonophysics 2003;370:163–76.

[11] Kachanov M. Effective elastic properties of cracked solids: critical review of some basic concepts. Appl Mech Rev 1992;45(8):

304–35.

[12] Kachanov M. Elastic solids with many cracks and related problems. In: Hutchinson JW, Wu T, editors. Advances in applied

mechanics, vol. 30. Academic Press; 1993. p. 259–445.

[13] Budiansky B, O’Connell RJ. Elastic moduli of a cracked solid. Int J Solids Struct 1976;12:81–97.

[14] Gottesman T, Hanshin Z, Brull MA. Effective elastic moduli of cracked fiber composites. In: Bussell AR, editor. Adv Compos Mater

1980:749–58.

[15] Benveniste Y. On the Mori–Tanaka’s method in cracked bodies. Mech Res Commun 1986;13:193–201.

[16] Mori T, Tanaka K. Average stress in matrix and average elastic energy of materials with misfitting inclusions. Acta Metall

1973;21:571–4.

[17] Hashin Z. The differential scheme and its application to cracked materials. J Mech Phys Solids 1988;36(6):719–34.

[18] Huang Y, Hu KX, Chandra A. A generalized self-consistent mechanics method for microcracked solids. J Mech Phys Solids

1994;42(8):1273–91.

[19] Santare MH, Crocombe AD, Anlas G. Anisotropic effective moduli of materials with microcracks. Engng Fract Mech

1995;52(5):833–42.

[20] Greengard L, Helsing J. On the numerical evaluation of elastostatic fields in locally isotropic two-dimensional composites. J Mech

Phys Solids 1998;46(8):1441–62.

[21] Helsing Johan. Fast and accurate numerical solution to an elastostatic problem involving ten thousand randomly oriented cracks. Int

J Fracture 1999;100:321–7.

[22] Wang YB, Chau KT. A new boundary element for plane elastic problems involving cracks and holes. Int J Fracture 1997;87(1):

1–20.

[23] Pan EN. A general boundary element analysis of 2-D linear elastic fracture mechanics. Int J Fracture 1997;88(1):41–59.

[24] Wang GS, Feng XT. The interaction of multiple rows of periodical cracks. Int J Fracture 2001;110(1):73–100.

[25] Dong CY, Lee Kangyong. Numerical analysis of doubly periodic array of cracks/rigid-line inclusions in an infinite isotropic medium

using the boundary integral equation method. Int J Fracture 2005;133:389–405.

[26] Feltman RS, Santare MH. Anisotropic effective moduli of cracked short-fiber reinforced composites. J Appl Mech 1999;66:1–5.

[27] Zheng Q-S, Du D-X. An explicit and universally applicable estimate for the effective properties of multiphase composites which

accounts for inclusion distribution. J Mech Phys Solids 2001;49:2765–88.

[28] Feng X-Q, Qin Q-H, Yu S-W. Quasi-micromechanical damage model for brittle solids with interacting microcracks. Mech Mater

2004;36:261–73.

[29] Kachanov M. Elastic solids with many cracks: a simple method of analysis. Int J Solids Struct 1987;23(1):23–43.

[30] Huang Y, Chandra A, Jiang ZQ, Wei X, Hu KX. The Numerical calculation of two-dimensional effective moduli for microcracked

solids. Int J Solids Struct 1996;33(11):1575–86.

[31] Shen L, Yi S. An effective inclusion model for effective moduli of heterogeneous materials with ellipsoidal inhomogeneities. Int J

Solids Struct 2001;38:5789–805.

[32] Shen L, Li J. A numerical simulation for effective elastic moduli of plates with various distributions and sizes of cracks. Int J Solids

Struct 2004;41:7471–92.

[33] Wang YB, Sun YZ. A new boundary integral equation method for cracked 2-D anisotropic bodies. Engng Fract Mech 2005;72:

2128–43.

1454 D. Su et al. / Engineering Fracture Mechanics 74 (2007) 1436–1455
[34] Walsh JB. The effect of cracks on the uniaxial elastic compression of rocks. J Geophys Res 1965;70:399–411.

[35] Taya M, Mura T. On the stiffness and strength of aligned short fiber reinforced composites containing fiber end cracks under uniaxial

applied stress. ASME J Appl Mech 1981;48:361–7.

[36] Anlas G, Santare MH. An arbitrarily oriented crack inside an elliptical inclusion. J Appl Mech 1993;60:589–94.

[37] ABAQUS analysis user’s manual, 2003. ABAQUS, INC.

[38] Auld BA. Acoustic fields and waves in solids, vol. I. Malabar, FL: Krieger Publ. Co.; 1990. 435 pp.

D. Su et al. / Engineering Fracture Mechanics 74 (2007) 1436–1455 1455
An effective medium model for elastic waves in microcrack damaged media

Dan Su a, Michael H. Santare a,*, George A. Gazonas b

aDepartment of Mechanical Engineering and Center for Composite Materials, University of Delaware, Newark, DE 19716, USA

bUS Army Research Laboratory, Weapons and Materials Research Directorate, Aberdeen Proving Ground, MD 21005, USA

article info

Article history:

Received 29 June 2007

Received in revised form 3 April 2008

Accepted 7 April 2008

Available online 14 April 2008

Keywords:

Waves

Anisotropic damage

Numerical experiment

Effective moduli

Microcracks

abstract

Direct numerical simulations of waves traveling through microcrack-damaged media are

conducted and the results are compared to effective medium calculations to determine

the applicability of the latter for studying wave propagation. Both tensile and compressive

waves and various angular distributions of randomly-located cracks are considered. The

relationships between the input wavelength and the output wave speed and output signal

strength are studied. The numerical simulations show that the wave speed is nearly constant when 1/ka > 60 for tensile waves and 1/ka > 10 for compressive waves, where k is

the wave number and a is the average half-crack length. The direct simulations also show

that when the input wavelength is much longer than the crack length, 1/ka > 60, the wave

can pass through the damaged medium relatively unattenuated. On the other hand, when

the input wavelength is shorter than a ‘‘cut off” wave length, the output wave magnitude

decreases linearly with the input wavelength. The effective medium wave speed and magnitude calculations are not dependent on the input wavelength and therefore the results

correspond well with the numerical simulations for large 1/ka. This suggests a minimum

wavelength for which the homogenized methods can be used for studying these problems.

2008 Elsevier Ltd. All rights reserved.

1. Introduction

For many brittle materials, the most common damage mechanism is microcracking. As microcracks develop, they change

the local mechanical response to input loads. One way, commonly used to evaluate the quasi-static response of the damaged

medium is by directly simulating an actual array of cracks. Solutions of this type include the fast multipole method (FMM,

[1,2]), boundary element method (BEM [3,4]) and many other numerical solutions, such as described in [5–9]. The major

advantage of these methods is that they take into account the actual crack–crack interactions. There is no theoretical limitation to the crack density or crack distribution and these direct methods can give highly accurate results if the model is

an accurate representation of the damaged medium. However, the detailed calculation of the interaction between the cracks

is also the major disadvantage of the direct methods. As the crack density increases, the modeling time and computational

demands can increase dramatically.

To address this increase in computational demand, effective medium methods can be used. In these methods, one models

the quasi-static response of a medium with diffuse microcrack damage by replacing the damaged material with an effective

medium, possessing the same local mechanical properties. Many such effective medium models have been developed, such

as the self consistent method (SCM) [10–12], Generalized self consistent method (GSCM) [13–17], Mori–Tanaka method

(MTM) [18,19], and the differential scheme method (DSM) [20]. For a comprehensive review of the literature in effective

mechanical properties prediction, readers can refer to [6,21].

0013-7944/$ - see front matter 2008 Elsevier Ltd. All rights reserved.

doi:10.1016/j.engfracmech.2008.04.003

* Corresponding author. Tel.: +1 302 831 2246; fax: +1 302 831 3619.

E-mail address: santare@me.udel.edu (M.H. Santare).

Engineering Fracture Mechanics 75 (2008) 4104–4116

Contents lists available at ScienceDirect

Engineering Fracture Mechanics

journal homepage: www.elsevier.com/locate/engfracmech
While the above-cited literature primarily discusses the quasi-static response of microcrack-damaged media, several

researchers have addressed elastic wave propagation behavior in these materials. In 1995 Sayers and Kachanov [22] used

the effective moduli of a cracked medium to study microcrack-induced elastic wave anisotropy. In that paper, the authors

used second-rank and fourth-rank crack density tensors to evaluate the effective elastic moduli. They calculated the effective

moduli for the (randomly-oriented) isotropic distribution of cracks and the perfectly aligned cracks. For the isotropic cracked

medium, they also back-calculated the crack density from experimental data of the ultrasonic wave speed. Following Sayers

and Kachanov’s results, Schubnel and Gueguen [23,24] used micromechanical and statistical physics to evaluate the elastic

wave velocities and permeability of cracked rocks. They calculated the velocity for high and low frequency waves for both

wet and dry cracks for an aligned array of cracks and an array composed of ±60 cracks under an external confining pressure.

Maurel et al. [25] studied the elastic wave propagation through a random array of dislocations. Markov [26] used the Frenkel–

Biot theory to study the longitudinal harmonic wave propagation in an isotropic porous matrix with inclusions. Levin and

Markov [27] used an effective medium approximation (EMA), based on the self-consistent method, to determine the effective

elastic moduli and elastic wave propagation velocities in a transversely isotropic solid containing aligned inclusions.

Yet, with all this literature, there are several questions which remain. For example, under what conditions can we model

an actual damaged material with an equivalent effective medium? How accurate is the effective medium representation?

What factors affect the accuracy?

Several researchers have performed studies to address some of these topics. Zhang and Gross [28] studied the propagation

of an elastic wave through a medium with randomly distributed cracks. The cracks are treated as finite length lines with displacement discontinuity and the crack surfaces are assumed to be frictionless. They studied the scattering function for a single crack and expanded this using a numerical approach to study the response of multiple cracks. They used both the theory

of Foldy and the causal approach based on K–K relations to compute the attenuation coefficient and the effective wave velocity. Littles et al. [29] performed an experimental and theoretical investigation to study how longitudinal waves are scattered

by a distribution of cracks. The results show that the transmission coefficients are a function of the incident wave number,

the crack size and crack density. Kawahara and Yamashita [30] studied the scattering of P, SV and SH waves by a zonal distribution of uniaxial cracks using a theoretical analysis. They showed that the attenuation coefficient peaks around ka 1,

the phase velocity is almost independent of k in the range of ka < 1 and increases monotonically when k is in the range of

ka > 1. They also found that the effect of crack-face friction is to shift the peak of the attenuation coefficient toward the lower

wave number range. Kelner et al. [31] used the boundary element method to simulate P wave propagation in media with

uniaxial cracks. The authors conclude that when the wavelength of the incident wave is close to, or shorter than, the crack

length, the scattered waves are efficiently excited. They also observed that the attenuation factor of the direct P wave peaks

at around ka = 2, where k is the incident wave number and a is the crack length, and decreases proportionally with (ka)

1 in

the high wave number range.

Most of the above papers focus on low crack densities and high frequency input waves, a combination which corresponds

to high attenuation. Also, the crack orientations in most of the papers are limited to aligned or randomly-oriented distributions. Additionally, in most of the literature cited, the authors assume the crack surfaces are stress free, which means the

cracks are opened under all loads. In our previous paper [16], we use the generalized self consistent method (GSCM) to predict the anisotropic effective moduli of a medium containing an arbitrarily-oriented distribution of cracks. In that paper, we

show that because the cracks open or close under different external loads, the effective moduli under different loading conditions are quite different. In that paper, we discussed three different loading conditions: (1) overall applied tensile loads; (2)

overall compressive loads; (3) initially stress free media. The effective moduli under these three different loading conditions

are evaluated taking into account the effects of crack-face contact. In that paper we also mentioned that tensile and compressive wave propagation in the same medium behave very differently. In this paper, we study under what conditions

one may use the calculated effective moduli to simulate wave propagation in the microcrack damaged medium. We conduct

direct numerical simulations to relate input wavelength to propagating wave speed and output signal strength, and compare

these results to those obtained via the effective medium models. Once these comparisons are established, we can establish

guidelines for the use of effective medium calculations to analyze ultrasound data for determining the internal state of

microcrack damage.

2. FEM model

In the following numerical experiments, we use the definition of crack density and crack orientation distribution function

/(h) described in [15,16]. We define the crack density as

g 1⁄4 Ma2

A ; ð1Þ

where M is the number of cracks per unit area A, and a is the average half-crack length. We assume the crack orientation can

be described by an orientation distribution function /(h), where h is the angle of the individual crack with respect to the

average crack orientation, assigned to be the x-axis. In the following analysis, we assume that the cracks are evenly

distributed between the angles h0 and +h0, as shown in Fig. 1. A more detailed description of the crack distribution function

/(h), and crack density function g, can be found in [16].

D. Su et al. / Engineering Fracture Mechanics 75 (2008) 4104–4116 4105
To simulate uni-directional plane wave propagation, we assume an infinite half plate, and apply a time-dependent load

uniformly along the top surface (Fig. 2a). Because the plate is infinitely extended in the x-direction, we can pick a strip from

anywhere in the half plane as a representative sample (Fig. 2b) by applying symmetry boundary conditions to the left and

right edges of the strip. The location of each crack within the strip was generated by a random number generator as was the

crack angle h chosen to be in the interval +h0 and h0. If the randomly generated location of a crack caused it to overlap with

an already existing crack, that crack was eliminated from the model and the process was continued until the desired crack

density was reached. Cracks that contacted the left or right boundaries were kept in the model, but because of the symmetry

conditions imposed on these boundaries, their lengths could be quite different from the average. Building the model in this

manner, the array is not completely random, but represents a relatively homogeneous approximation to randomly distributed cracks. Contact elements are placed along the cracks to prevent crack-face interpenetration. Fig. 2b is a representation

of the models used in the numerical experiments. A number of trial models were run to determine an appropriate length and

width for the strips and the chosen dimensions represent a compromise between a model large enough to eliminate edge

effects and one small enough to allow for efficient computations. All lengths are normalized with respect to the average crack

length (2a) and the resulting width of each model is 4, and the length is 80. In all the results that follow, the crack density

parameter is chosen to be g = 0.4. (Most of the previously reviewed literature studies crack densities in the range of 0.1–0.2

when defined by the term g, in Eq. (1)) Larger densities proved difficult to model, due to interference between individual

cracks and smaller densities show the same overall trends, but to a lesser extent. Since the cracks are randomly located

and oriented, we constructed four or five models for each orientation distribution to show the variability introduced by

the randomness.

In the simulation, we use ABAQUS 6.4.4 as the finite element solver. In Figs. 3–6, we show the finite element models for

samples 1–14. As shown in Fig. 3d, each crack is represented by two crack surfaces with zero distance between them. Therefore, in the stress free state, all the cracks are closed.

In order to study the relationship between the input wavelength and the output signal, we apply a single frequency timedependent load, uniformly along the top of the strip. Fig. 7 shows a sample input load with the function

PðtÞ 1⁄4 A 1  cosðxtÞ

2

; for 0 6 xt

2p 6 1; ð2Þ

where, P(t) is the time-dependent load, A is the maximum load amplitude, x is the circular frequency of the input and t is the

time. We can apply positive or negative pressure in order to generate a compressive or tensile wave inside the sample.

In the numerical experiments, we use the dimensionless number 1/ka as the effective wave length parameter, where a is

the average half-crack length, k = x/cp, x is the input wave frequency and cpis the plane-stress longitudinal wave velocity in

the un-damaged medium [32], which is defined as

−θ 0 +θ 0

φ(θ )

x

θ

crack

Fig. 1. Even distribution of crack orientation.

b

Distance

from the 

top

a

y

x

Fig. 2. Model schematics. (a) Time-dependent load is applied along the top boundary of the half plane. (b) A representative strip from the plane is used in

the finite element simulation.

4106 D. Su et al. / Engineering Fracture Mechanics 75 (2008) 4104–4116
cp 1⁄4

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

E

ð1  m2Þq

s

: ð3Þ

This effective wave number 1/ka, therefore represents the ratio between the input wave length in the undamaged media and

the average half crack length.

Fig. 3. Finite element model of sample 06 (h0 = 30). (a) Whole model without mesh detail. (b1, b2) Zoom in for the top part with and without mesh detail.

(c1, c2) Further zoom in with and without mesh detail. (d1, d2) Zoom in for a single crack with and without mesh detail.

Fig. 4. The top part of finite element model of samples 01–05 (h0 = 0).

D. Su et al. / Engineering Fracture Mechanics 75 (2008) 4104–4116 4107
Fig. 5. The top part of finite element model of samples 06–09 (h0 = 30).

Fig. 6. The top part of finite element model of samples 10–14 (h0 = 60).

4108 D. Su et al. / Engineering Fracture Mechanics 75 (2008) 4104–4116
3. Wave speed calculation

In this section we determine the wave speeds in the direct numerical simulations and compare these to the speeds in the

effective media models. The properties of the effective media, and the generalized-self-consistent method used to calculate

them, are shown in [16]. All the wave speeds are normalized by the plane-stress wave speed in the intact material cp, shown

in Eq. (3). We trigger the wave from the top of the models with time dependent pressure loads and trace the stress component ry. We chose 1% of the maximum input as the threshold value for the stress at the wave front. This means that we assume that the wave arrives at a particular cross section when the average stress ry in that cross section reaches 1% of the

maximum input stress.

In Figs. 8–10, we show the wave speed in each of the direct numerical simulation models as a function of the input wave

length number 1/ka, as well as the effective medium predictions and an analytical solution. The analytical solution is the

uniaxial longitudinal wave speed, which can be derived from [33], for an anisotropic plane stress medium calculated from

cp2 1⁄4

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

E2

ð1  m12m21Þq

s

; ð4Þ

where E2, m12, m21 are the corresponding effective Young’s moduli in y direction, and the effective Poisson’s ratios, which are

shown in [16] and q is the corresponding mass density. Fig. 8 shows results for the five different models with an aligned

array of cracks h0 = 0. Fig. 9 contains the models with crack orientations evenly distributed between the values h0 = ±30

and Fig. 10, the models with angles between h0 = ±60. The left column in each figure shows the tensile wave speed, while

the right column shows the speed of compressive waves.

0

0.2

0.4

0.6

0.8

1

0 0.2 0.4 0.6 0.8 1

ωt/(2π)

P/A

Fig. 7. Sample input load. The function of the input load is P 1⁄4 A 1cosðxtÞ

2 for 0 6 xt

2p 6 1.

ο θ0 = 0

0 20 40 60 80 100 120 140

0.0

0.2

0.4

0.6

0.8

1.0

Normalized Speed

1/ka 1/ka

 Average Speed (Experiments)

 Effective Media (FEM)

 Analytical Solution

 Sample 01

 Sample 02

 Sample 03

 Sample 04

 Sample 05

0 20 40 60 80 100 120 140

0.0

0.1

0.2

0.3

0.4

0.5 

Normalized Speed

 Average Speed (Experiments)

 Effective Media (FEM)

 Analytical Solution

 Sample 01

 Sample 02

 Sample 03

 Sample 04

 Sample 05

Fig. 8. Normalized wave speed vs. 1/ka for h0 = 0. Left : speed of tensile waves. Right : speed of compressive wave.

D. Su et al. / Engineering Fracture Mechanics 75 (2008) 4104–4116 4109
From the results, we can see that in all cases, the compressive waves travel much faster than the tensile waves. This is

because under compression, the cracks can only open in mode II (shear opening) (As we discuss in Section 2, in the

stress-free state, the distance between the cracks surfaces is zero. Therefore, the crack surfaces cannot move closer to each

other under compression.) while in tension, the cracks can open in both mode I (normal opening) and mode II. Therefore for a

given medium, the material is much less stiff when tensile waves pass through, than when compressive waves pass through.

For a detailed discussion regarding the different effective material properties under different loads, the reader is referred to

[16]. The figures also show that the effective medium models correspond well with the numerical experiments in the cases

considered. For the case of tensile wave speed when h0 = 0, the analytical equations and the effective media model give

higher wave speed predictions than the direct numerical simulations. This is because for tensile waves, when h0 = 0, the

cracks in the numerical simulations cause significant energy dispersion as the wave travels along the strip. Therefore, for

the same input magnitude, it takes longer for the average stress across the strip to reach the 1% threshold value, resulting

in a slightly lower calculated speed for the direct simulations.

Kawahara and Yamashita [30] found analytical solutions for a similar problem containing randomly distributed, aligned

cracks under confining pressure. They found that when the crack surfaces are friction free, the phase velocity is nearly independent of k when 1/ka > 1 and increases monotonically with k when 1/ka < 1. From the right column of Figs. 8–10 we can

see that the wave speed in our numerical experiments follows this same trend and are completely consistent with the anaο θ0 = 30

0 20 40 60 80 100 120 140

0.0

0.1

0.2

0.3

0.4

0.5

Normalized Speed

1/ka

 Average Speed (Experiments)

 Effective Media (FEM)

 Analytical Solution

 Sample 06

 Sample 07

 Sample 08

 Sample 09

0 20 40 60 80 100 120 140

0.0

0.2

0.4

0.6

0.8

1.0

Normalized Speed

1/ka

 Average Speed (Experiments)

 Effective Media (FEM)

 Analytical Solution

 Sample 06

 Sample 07

 Sample 08

 Sample 09

Fig. 9. Normalized wave speed vs. 1/ka for h0 = 30. Left figure: speed of tensile waves. Right figure: speed of compressive wave.

ο 60 θ0 =

0 20 40 60 80 100 120 140

0.0

0.1

0.2

0.3

0.4

0.5

Normalized Speed

1/ka

0 20 40 60 80 100 120 140

1/ka

 Average Speed (Experiments)

 Effective Media (FEM)

 Analytical Solution

 Sample 10

 Sample 11

 Sample 12

 Sample 13

 Sample 14

0.0

0.2

0.4

0.6

0.8

1.0

Average Speed

 Average Speed (Experiments)

 Effective Media (FEM)

 Analytical Solution

 Sample 10

 Sample 11

 Sample 12

 Sample 13

 Sample 14

Fig. 10. Normalized wave speed vs. 1/ka for h0 = 60. Left : speed of tensile waves. Right : speed of compressive wave.

4110 D. Su et al. / Engineering Fracture Mechanics 75 (2008) 4104–4116
Fig. 11. Tensile waves with different effective wave number traveling in sample 05, g = 0.4, h0 = 0. Normalized ry is the average stress ry in the specific

cross session normalized by the maximum average stress ry at distance equal to zero (the top boundary). The distance from the top is shown in Fig. 2b.

Fig. 12. Compressive waves with different effective wave number traveling in sample 05, g = 0.4,h0 = 0. Normalized ry is the average stress ry in the specific

cross session normalized by the maximum average stress ry at distance equal to zero (the top boundary). The distance from the top is shown in Fig. 2b.

D. Su et al. / Engineering Fracture Mechanics 75 (2008) 4104–4116 4111
lytical results in Kawahara and Yamasita’s paper for cracks under compression. At the same time, from the left column in

Figs. 8–10 we see that the tensile wave speed is nearly constant when 1/ka > 60 and increases monotonically when 1/

ka < 60. Therefore, we can see that the Kawahara and Yamasita’s general conclusion still applies to tensile wave propagation

in our current study, although the cut off ka (or 1/ka) value is much different for tensile waves than compressive waves.

4. Signal strength

In the direct numerical simulations, we applied a uniform, time varying load along the top boundary of the finite element

models, as shown in Fig. 2b. For a specific time and a specific cross section in the model, we record the average stress component ry through this cross section (Fig. 2b) then normalize this value by the maximum average stress ry calculated at the

top boundary of the model. (Recall that the crack length 2a equals 1, except for the cracks that touch the boundary.) The

properties of the intact material are E/l = 8/3, m = 1/3, where E is Young’s modulus, l is shear modulus and m is Poisson’s ratio.

The normalized time is defined as t

_

1⁄4 t

2a=cp 1⁄4 tcp=ð2aÞ.

Fig. 11 shows the results for the tensile waves with different effective wave numbers traveling in sample 05. We can see

that for tensile waves, when the input wavelength is much larger than the crack length, (1/ka > 10), the input signal passes

through the damaged medium relatively unattenuated (the first row of the figures). When the input wavelength is smaller,

(1/ka = 8.24 and 3.53) the cracks show a larger ‘‘blocking” effect and the transmitted wave shows smaller stress amplitude.

When the input wavelength is about the same as the crack length, (1/ka = 0.706) the signal is almost completed blocked by

the damage.

0.0

0.2

0.4

0.6

0.8

1.0

1.2

A10/A0

 Sample 01

 Sample 02

 Sample 03

 Sample 04

 Sample 05

0 10 20 30 40 50 60 70

0.0

0.2

0.4

0.6

0.8

1.0

1.2

A20/A0

1/ka

 Sample 01

 Sample 02

 Sample 03

 Sample 04

 Sample 05

0.0

0.2

0.4

0.6

0.8

1.0

1.2

A10/A0

 Sample 01

 Sample 02

 Sample 03

 Sample 04

 Sample 05

0 20 40 60 80 100 120

0.0

0.2

0.4

0.6

0.8

1.0

1.2

A20/A0

1/ka

0 10 20 30 40 50 60 70

1/ka

0 20 40 60 80 100 120

1/ka

 Sample 01

 Sample 02

 Sample 03 

 Sample 04

 Sample 05

ο 0 θ0 =

Fig. 13. Output signal strength vs. input effective wave number (1/ka) for the models with g = 0.4, h0 = 0 and 2a = 1. Left column: Tensile waves output.

Right column: Compressive waves output. First row: The normalized peak averaged stress measured at distance equal to 10. Second row: The normalized

peak averaged stress measured at distance equal to 20.

4112 D. Su et al. / Engineering Fracture Mechanics 75 (2008) 4104–4116
Fig. 12 shows the compressive wave propagation through the same model sample 05. Comparing Fig. 11 with Fig. 12, we

can see that the compressive wave signal does not decrease appreciably, even when the input wavelength (1/ka) is relatively

small. This is because when tensile waves pass through the sample, the cracks are opened and parts of the waves are reflected by the open crack surfaces, causing the output signal to decrease. Yet, when compressive waves travel in the model,

all the cracks remain closed and therefore have almost no effect on the wave propagation in this situation.

5. Output signal strength

5.1. Uniform crack length

In this section, we show how the peak stress values at different cross sections change with varying input wavelength. In

Figs. 13–15, A0, A10 and A20 represent the peak average stress ry at cross sections with normalized distances of 0, 10 and 20,

respectively from the top of the model. Fig. 13 shows results for the models with an aligned array of cracks h0 = 0 and Figs.

14 and 15 contain the models with cracks oriented between the values h0 = ± 30 and h0 = ± 60 respectively. As before, the

left column in each figure shows the results for the tensile waves, while the right column shows the results for the compressive waves.

In Figs. 13–15, the individual lines represent the different model samples. In all of these models g = 0.4 and 2a = 1 (except

for the cracks that touch the boundaries). All of these models show the same phenomenon observed in the previous section

ο θ0 = 30

 Sample 06

 Sample 07

 Sample 08

 Sample 09

 Sample 07

 Sample 08

 Sample 09

 Sample 10

 Sample 07

 Sample 08

 Sample 09

 Sample 10

 Sample 06

 Sample 07

 Sample 08

 Sample 09

0.0

0.2

0.4

0.6

0.8

1.0

1.2

A10/A0

0 10 20 30 40 50 60 70

0.0

0.2

0.4

0.6

0.8

1.0

1.2

A20/A0

1/ka

0.0

0.2

0.4

0.6

0.8

1.0

1.2

A10/A0

0 20 40 60 80 100 120

0.0

0.2

0.4

0.6

0.8

1.0

1.2

A20/A0

1/ka

0 10 20 30 40 50 60 70

1/ka

0 20 40 60 80 100 120

1/ka

Fig. 14. Output signal strength vs. input wave length number (1/ka) for the models with g = 0.4, h0 = 30 and 2a = 1. Left column: Tensile waves output.

Right column: Compressive waves output. First row: The normalized peak averaged stress measured at distance equal to 10. Second row: The normalized

peak averaged stress measured at distance equal to 20.

D. Su et al. / Engineering Fracture Mechanics 75 (2008) 4104–4116 4113
when we discussed the results of sample 05. What these figures bring out, is the variability from the different arrays of cracks

(different experimental samples) as well as more clearly delineating the signal strength–wavelength relationship.

For tensile waves, when the input wave length is much larger than the crack length (1/ka > 20), the input signal passes

through the microcrack damaged medium relatively easily. The output signal strength remains nearly constant for all input

wavelengths above (1/ka > 20), and the output signal strengths are usually more than 80% of the input signal strength. On the

other hand, the cracks show larger ‘‘blocking” effect when the input wave length decreases. In the area where 1/ka < 15, the

output signal strength decreases almost linearly with the input 1/ka decreases. Therefore, we can think of 1/ka = 15 as a ‘‘cut

off” wave length below which tensile waves will not travel through a microcrack damaged medium. When the input wavelength is close to or smaller than the crack length (1/ka 1), almost no signal can pass through the medium.

For compressive waves, we can see the similar phenomenon but there are two differences that distinguish the compressive wave output from the tensile wave output. First, when the 1/ka number is large, the output signal of the compressive

waves is much more uniform than that of the tensile waves. Secondly, the ‘‘cut off” wave length for the compressive waves is

much smaller than that for the tensile waves. In this study we see that for compressive waves, the output signal decreases

only when 1/ka < 5.

These results suggest that there is a ‘‘safe range” for modeling microcrack-damaged media with homogenized effective

medium approximations. When the input wavelength is larger than the ‘‘cut off” wavelength, the signal can travel relatively

unattenuated through the medium. In this range, the effective medium FE models give a very good approximation to the

actual medium with damage. On the other hand, applying the homogenized methods in the range with input wavelength

ο θ 0 = 60

 Sample 10

 Sample 11

 Sample 12

 Sample 13

 Sample 14

 Sample 11

 Sample 12

 Sample 13

 Sample 14

 Sample 15

 Sample 11

 Sample 12

 Sample 13

 Sample 14

 Sample 15

 Sample 10

 Sample 11

 Sample 12

 Sample 13

 Sample 14

0.0

0.2

0.4

0.6

0.8

1.0

1.2

A10/A0

0 10 20 30 40 50 60 70

0.0

0.2

0.4

0.6

0.8

1.0

1.2

A20/A0

1/ka

0.0

0.2

0.4

0.6

0.8

1.0

1.2

A10/A0

0 20 40 60 80 100 120

0.0

0.2

0.4

0.6

0.8

1.0

1.2

A20/A0

1/ka

0 10 20 30 40 50 60 70

1/ka

0 20 40 60 80 100 120

1/ka

Fig. 15. Output signal strength vs. input wave length number (1/ka) for the models with g = 0.4, h0 = 60 and 2a = 1. Left column: Tensile waves output.

Right column: Compressive waves output. First row: The normalized peak averaged stress measured at distance equal to 10. Second row: The normalized

peak averaged stress measured at distance equal to 20.

4114 D. Su et al. / Engineering Fracture Mechanics 75 (2008) 4104–4116
below the ‘‘cut off” wavelength is inappropriate, unless we modify the homogenized methods to take the wave length dependent energy dissipation into account.

We also expect the ‘‘cut off” 1/ka to decrease when the crack density decreases. This is because when the crack density

decreases, there are fewer cracks per unit area. Therefore, less energy will be dissipated by the cracks when the wave passes

through the media.

On a PC with a Pentium 4, 2.0 GHz CPU and 1GB RAM, it took from 3 to 30 h for a single point in the tensile Figs. 13–16

depending on the crack distribution and the time frame of simulation for the wave propagation. For the compressive cases,

the time varied from 2 to 5 h.

5.2. Distributed crack lengths

Fig. 16 shows the output signal strength as a function of 1/ka for a set of models with h0 = 0 and average crack length 2a.

These are different than the previous models however, because the individual crack length is variable. For simplicity, we

choose an even distribution of crack lengths between 0.6 and 1.4 with the average of 1.0, (again, except for the cracks that

touch the boundaries). Comparing Fig. 13 with Fig. 16, we can see that for the models with distributed crack lengths, the

signal strength values have a wider band for the range of 10 <1/ka< 60. Outside of this region, the results in Figs. 13 and

16 are quite similar to each other. We can conclude from this result, that a distributed crack length will not have a large

effect on the ‘‘cut off” wavelength and the ‘‘safe range” for using the homogenized effective medium calculations.

From the above figures we can see that in most of the cases the output signal strength is weaker than the input signal

strength. The reason is obvious: the reflection and diffusion between the cracks reduce the average signal strength. However,

in some regions of several samples, we see that the output signal is actually stronger than the input signal. This can be explained by the random distribution of the cracks. When the wave travels through a field of cracks in the direct numerical

simulations, part of the wave is reflected by the crack surfaces and interferes with other components of the wave which follow. Therefore, the local distributions of the cracks in the models can actually increase the local stress wave amplitude in

some regions.

6. Conclusions and discussion

In this paper, we performed direct numerical simulations of waves traveling in micro-crack damaged media and compared these to results using a homogenized effective medium calculation. For the numerical simulations, we built several

models with different crack distributions for each of three angular distributions; h0 = 0, h0 = 30 and h0 = 60. We studied

the wave speed variation as a function of the input wavelength in each model for both tensile waves and compressive waves.

The results for wave speed in the effective medium models are close to the numerical simulation results.

When a plane wave approaches an obstacle, the obstacle will reflect (or block) a portion of the advancing wave energy. If

the obstacle is small with respect to the wavelength of the advancing wave, there will be a small disturbance in the wave

front, but the majority of the energy will pass. Essentially, the wave travels around the obstacle. However, if the obstacle

is large with respect to the wavelength, it will reflect a significant portion of the wave, detracting from the energy remaining

in the transmitted wave. These observations are supported by Eringen’s [34] analytical solutions of scattering of a plane

0 20 40 60 80 100 120

0.0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

A10/A0

1/ka

 Sample 15

 Sample 16

 Sample 17

 Sample 18

 Sample 19

0 20 40 60 80 100 120

0.0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

A20/A0

1/ka

 Sample 15

 Sample 16

 Sample 17

 Sample 18

 Sample 19

Fig. 16. Output signal strength v.s. 1/ka for tensile waves in the models with h0 = 0. The crack length 2a is evenly distributed between 0.6 and 1.4, with

average of 1.0. Left figure: The normalized peak averaged stress measured at distance equal to 10. Right figure: The normalized peak averaged stress

measured at distance equal to 20.

D. Su et al. / Engineering Fracture Mechanics 75 (2008) 4104–4116 4115
wave with a single cylindrical obstacle. Eringen found that the total scattered power increases with decreasing wavelength

(ka = 2pa/k) of the incident plane wave. Here, a is the radius of the cylinder and k is the wavelength of the incident plane

wave. For incident waves with large wavelengths (large 1/ka) the scattered elastic energy approaches zero and the wave does

not ‘‘see” the obstacle; we have determined that for crack systems modeled using our finite element approach this occurs for

1/ka > 60 for media in tension, and 1/ka > 10 for those in compression. In the case of distributed obstacles (such as the microcracks in the present study), the effect is more evident since the plane wave energy that passes by the obstacles (cracks), will

be further blocked by other obstacles behind the first ones.

We also studied the relationships between the input wavelength and output signal strength. These results show that for

both tensile and compressive waves, there is a ‘‘cut off” wavelength above which the wave can easily pass through the damaged medium and that this ‘‘cut off” is different in tension and compression. When the input wavelength is shorter than the

‘‘cut off” wave length, the output wave strength decreases linearly with the input wavelength. When the input wavelength is

close to average crack length, (1/ka = 1) the wave is highly diffused. This suggests a ‘‘safe range” for using the homogenized

effective media methods in studying these problems without considering the wave length dependent energy diffusion.

Acknowledgement

This research was supported by the United States Army Research Laboratory through the Composite Materials Technology cooperative in agreement with the Center for Composite Materials at the University of Delaware.

References

[1] Greengard L, Helsing J. On the numerical evaluation of elastostatic fields in locally isotropic two-dimensional composites. J Mech Phys Sol

1998;46(8):1441–62.

[2] Helsing J. Fast and accurate numerical solution to an elastostatic problem involving ten thousand randomly oriented cracks. Int J Fract

1999;100(4):321–7.

[3] Wang YB, Chau KT. A new boundary element for plane elastic problems involving cracks and holes. Int J Fract 1997;87(1):1–20.

[4] Pan EN. A general boundary element analysis of 2-D linear elastic fracture mechanics. Int J Fract 1997;88(1):41–59.

[5] Kachanov M. Elastic solids with many cracks – a simple method of analysis. Int J Sol Struct 1987;23(1):23–43.

[6] Kachanov M. Effective elastic properties of cracked solids: critical review of some basic concepts. Appl Mech Rev 1992;45(8):304–35.

[7] Huang Y, Chandra A, Jiang ZQ, Wei X, Hu KX. The numerical calculation of two-dimensional effective moduli for microcracked solids. Int J Sol Struct

1996;33(11):1575–86.

[8] Shen LX, Li J. A numerical simulation for effective elastic moduli of plates with various distributions and sizes of cracks. Int J Sol Struct

2004;41(26):7471–92.

[9] Shen LX, Yi S. An effective inclusion model for effective moduli of heterogeneous materials with ellipsoidal inhomogeneities. Int J Sol Struct

2001;38(32–33):5789–805.

[10] Budiansky B, O’Connell RJ. Elastic moduli of a cracked solid. Int J Sol Struct 1976;12(2):81–97.

[11] Gottesman T, Hanshin Z, Brull MA. Effective elastic moduli of cracked fiber composites. In: Bussell AR, editor. Advances in composite materials.

Proceedings of the ICCM 3, Pergamon Press, Oxford; 1980. p. 749–58.

[12] Horii H, Nemat-Nasser S. Overall moduli of solids with microcracks: load-induced anisotropy. J Mech Phys Sol 1983;31(2):155–71.

[13] Aboudi J, Benveniste Y. The effective moduli of cracked bodies in plane deformations. Engng Fract Mech 1987;26(2):171–84.

[14] Huang Y, Hu KX, Chandra A. A generalized self-consistent mechanics method for microcracked solids. J Mech Phys Sol 1994;42(8):1273–91.

[15] Santare MH, Crocombe AD, Anlas G. Anisotropic effective moduli of materials with microcracks. Engng Fract Mech 1995;52(5):833–42.

[16] Su D, Santare MH, Gazonas GA. The effect of crack face contact on the anisotropic effective moduli of microcrack damaged media. Engng Fract Mech

2007;74(9):1436–55.

[17] Feltman RS, Santare MH. Anisotropic effective moduli of cracked short-fiber reinforced composites. J Appl Mech–Trans ASME 1999;66(3):709–13.

[18] Mori T, Tanaka K. Average stress in matrix and average elastic energy of materials with misfitting inclusions. Acta Metall 1973;21(5):571–4.

[19] Benveniste Y. On the Mori-Tanaka’s method in cracked bodies. Mech Res Commun 1986;13(4):193–201.

[20] Hashin Z. The differential scheme and its application to cracked materials. J Mech Phys Sol 1988;36(6):719–34.

[21] Kachanov M. Elastic solids with many cracks and related problems. In: Hutchinson JW, Te Wu, editors. Advances in applied mechanics. Academic

Press; 1993. p. 259–445.

[22] Sayers CM, Kachanov M. Microcrack-induced elastic-wave anisotropy of brittle rocks. J Geophys Res–Sol Earth 1995;100(B3):4149–56.

[23] Schubnel A, Gueguen Y. Dispersion and anisotropy of elastic waves in cracked rocks. J Geophys Res–Sol Earth 2003;108(B2).

[24] Gueguen Y, Schubnel A. Elastic wave velocities and permeability of cracked rocks. Tectonophysics 2003;370(1–4):163–76.

[25] Maurel A, Mercier JF, Lund F. Elastic wave propagation through a random array of dislocations. Phys Rev B 2004;70(2).

[26] Markov MG. Propagation of longitudinal elastic waves in a fluid-saturated porous medium with spherical inclusions. Acoust Phys 2005;51:S115–21.

[27] Levin VM, Markov MG. Elastic properties of inhomogeneous transversely isotropic rocks. Int J Sol Struct 2005;42(2):393–408.

[28] Zhang C, Gross D. Wave attenuation and dispersion in randomly cracked solids. 1. Slit cracks. Int J Engng Sci 1993;31(6):841–58.

[29] Littles JW, Jacobs LJ, Qu JM. Experimental and theoretical investigation of scattering from a distribution of cracks. Ultrasonics 1995;33(1):37–43.

[30] Kawahara J, Yamashita T. Scattering of elastic-waves by a fracture-zone containing randomly distributed cracks. Pure Appl Geophys

1992;139(1):121–44.

[31] Kelner S, Bouchon M, Coutant O. Numerical simulation of the propagation of P waves in fractured media. Geophys J Int 1999;137(1):197–206.

[32] Cemal EA, Suhubi ES. Elastodynamics. Linear theory, vol. II. New York, San Francisco, London: Academic press, INC; 1975.

[33] B.A. Auld, Acoustic fields and waves in solids vol. I. New York, London, Sydney, Toronto;1990.

[34] Eringen AC. Mechanics of continua. R.E. Krieger Publishing Co.; 1989.

4116 D. Su et al. / Engineering Fracture Mechanics 75 (2008) 4104–4116
Anisotropic effective moduli of microcracked materials under

antiplane loading

Xu Wang a, Michael H. Santare a,*, George A. Gazonas b

aDepartment of Mechanical Engineering and Center for Composite Materials, University of Delaware, Newark, DE 19716, USA

bUS Army Research Laboratory, Weapons and Materials Research Directorate, Aberdeen Proving Ground, MD 21005, USA

article info

Article history:

Received 27 January 2009

Received in revised form 9 April 2009

Accepted 27 April 2009

Available online 5 May 2009

Keywords:

Anisotropic damage

Crack orientation distribution function

Effective moduli

Generalized self-consistent method

abstract

This study focuses on the prediction of the anisotropic effective elastic moduli of a solid

containing microcracks with an arbitrary degree of alignment by using the generalized

self-consistent method (GSCM). The effective elastic moduli pertaining to anti-plane shear

deformation are discussed in detail. The undamaged solid can be isotropic as well as anisotropic. When the undamaged solid is isotropic, the GSCM can be realized exactly. When the

undamaged solid is anisotropic it is difficult to provide an analytical solution for the crack

opening displacement to be used in the GSCM, thus an approximation of the GSCM is pursued in this case. The explicit expressions of coupled nonlinear equations for the unknown

effective moduli are obtained. The coupled nonlinear equations are easily solved through

iteration.

2009 Elsevier Ltd. All rights reserved.

1. Introduction

Microcracks are common defects in solids and multiple microcracks usually coexist in a single solid. The prediction of the

effective elastic properties of a microcracked solid is technically challenging and can find many practical applications.

The major methods developed so far to predict the effective elastic moduli of a microcracked solid include the following:

the self-consistent method (SCM) in which a crack is embedded directly into an effective medium [1]; the generalized

self-consistent method (GSCM) in which a crack is surrounded by an undamaged matrix region, and then embedded in

the effective medium [2,3]; the Mori-Tanaka method (MTM) [4]; the differential scheme method (DS) [5–7]; and the modified differential scheme (MDS) [8]. Recently the GSCM in conjunction with a finite element method (FEM) was developed to

take into account crack face contact and friction [9]. Most recently the representative unit cell approach was proposed by

Kushch et al. [10] to calculate effective elastic properties of a microcracked solid. Here it shall be noted that the effect of crack

orientation statistics on the anisotropic effective moduli of a microcracked solid was in fact first investigated by Santare et al.

[3] through the introduction of the crack orientation distribution function /(h) which was later adopted in [11] within the

framework of the SCM to study the problem of cracks with an arbitrarily degree of alignment in a material that is originally

anisotropic before the damage occurs.

In this research we analytically study the anisotropic effective elastic moduli of a solid containing microcracks with an

arbitrary degree of alignment by using the GSCM. Our model is in principle based on the GSCM developed by Santare

et al. [3]. In [3] the GSCM was used to predict the anisotropic moduli under plane stress loading. Here we are concerned with

the effective elastic moduli pertaining to anti-plane shear deformation. In our model the undamaged material can be isotropic as well as anisotropic. An exact solution to the cracked elliptical inclusion problem, which is essential in the realization of

0013-7944/$ - see front matter 2009 Elsevier Ltd. All rights reserved.

doi:10.1016/j.engfracmech.2009.04.010

* Corresponding author. Tel.: +1 302 831 2246; fax: +1 302 831 3619.

E-mail address: santare@udel.edu (M.H. Santare).

Engineering Fracture Mechanics 76 (2009) 1910–1919

Contents lists available at ScienceDirect

Engineering Fracture Mechanics

journal homepage: www.elsevier.com/locate/engfracmech
GSCM, can be derived when the undamaged material is isotropic. On the other hand an approximate analytical solution can

still be derived when the undamaged material is anisotropic and when the crack density is not very high.

2. The effective moduli of a microcracked solid

For a microcracked solid, the strain energy relationship between the effective medium and the actual microcracked solid

can be expressed as [2,3]

1

2 S

ijklr0

klr0ij 1⁄4 1

2 Sijklr0

klr0ij þ

1

2V

XM

k1⁄41

Z

Ck

1⁄2uit

0

i dCk; ð1Þ

where S

ijkl is the effective compliance of the damaged material, Sijkl is the compliance of the undamaged material and V is the

sample volume. The above integral is the energy dissipated through the opening of each microcrack, summed over all M

cracks, and r0

ij is the applied homogeneous stress field while t

0

i 1⁄4 r0

ijnj is the traction along the crack face if the crack did

not exist and [ui] is the crack opening displacement. In this research we focus on the two dimensional case in which all

the cracks penetrate the solid through the x3-direction. In addition we only discuss the effective elastic moduli pertaining

to anti-plane shear deformation. In the following we will address two cases: (i) the undamaged material is isotropic; (ii)

the undamaged material is anisotropic.

2.1. Isotropic undamaged material

The degree of crack alignment can be described by the crack orientation distribution function /(h) with h, (|h| < p/2) being

the angle between an individual crack and the positive x1-axis [3]. Without losing generality, /(h) can be taken as an even

function of h since we have assumed that the undamaged solid is isotropic. For simplicity, /(h) is specifically given by [3,9]

/ðhÞ 1⁄4

1

2h0 ; jhj  h0

0; jhj > h0

(

ð2Þ

where h0 6 p/2. The two special cases of perfectly aligned cracks and randomly oriented cracks correspond to h0 = 0 and

h0 = p/2, respectively, in Eq. (2).

Once we have introduced /(h), the summation in the energy relationship Eq. (1) can be written as an integral over orientation angle h,

1

2 S

ijklr0

klr0ij 1⁄4 1

2 Sijklr0

klr0ij þ

M

2A

Z h0

h0

/ðhÞ

Z

Ck

1⁄2uit

0

i dCkdh; ð3Þ

where A is a representative area of the sample.

In this study, we assume that the material is under anti-plane shear deformation. As a result the above energy relationship Eq. (3) can be simplified as

ðr0

31Þ

2

2C

55

þ ðr0

32Þ

2

2C

44

1⁄4 ðr0

31Þ

2 þ ðr0

32Þ

2

2l

þ g

2c2

Z h0

h0

/ðhÞ

Z

Ck

1⁄2u3t

0

3dCkdh; ð4Þ

where r0

31 and r032 are the anti-plane, far-field stresses, C

44 and C55 are the two relevant effective moduli of the damaged

material, l is the shear modulus of the undamaged material, c is the average half crack length and g = Mc2

/A is the crack

density parameter.

If the crack did not exist, the uniform traction due to the far-field stresses, t

0

3, along the line of the crack face is given by

t

0

3 1⁄4 cos hr032  sin hr031: ð5Þ

Next we introduce the GSCM [3] to approximately take into account the interaction between the cracks, as shown in

Fig. 1. The inclusion is assumed to have the properties of the undamaged material with known elastic moduli. The surrounding matrix is composed of the effective orthotropic, damaged material with the principal directions along the x1 and x2 axes.

The half-length of the crack is c, the semi-major and semi-minor axes of the ellipse are a and b, respectively. The crack density parameter g relates average crack length to the area of the ellipse, but in general, this leaves one of the three parameters

a, b and c, unspecified. Therefore, as an additional condition, we will require a2  b2 = c2 to be satisfied. This is the same relationship that was used in [3] for convenience, but here it is necessary in order to make analytical solutions possible. The twophase composite is subjected to uniform anti-plane shearing r0

31 and r032 at infinity.

By using the complex variable method [12,13], the crack opening displacement for the elliptical domain, depicted in

Fig. 1, [u3] can be obtained exactly as

1⁄2u3 1⁄4 4 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

c2  x2 p

l1⁄21 þ C  R2ð1  CÞ

alþ bC

55

ða þ bÞlcos hr0

32  alþ bC

44

ða þ bÞlsin hr0

31   ðjxj  cÞ ð6Þ

X. Wang et al. / Engineering Fracture Mechanics 76 (2009) 1910–1919 1911
where a 1⁄4 c

2 ðR þ R1

Þ; b 1⁄4 c

2 ðR  R1

Þ, (R > 1), C = l*/l with l1⁄4 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiC

44C55 p . Here the parameter R, which is introduced during

the analysis by using the complex variable method and conformal mapping, can be expressed in terms of a and b as R 1⁄4

ffiffiffiffiffiffi

aþb

ab

q

.

A detailed derivation can be found in the Appendix.

Consequently the integral in Eq. (4) can also be carried out exactly to arrive at the following expression

ðr0

31Þ

2

2C

55

þ ðr0

32Þ

2

2C

44

1⁄4 ðr0

31Þ

2 þ ðr0

32Þ

2

2l

þ g

4l

alþ bC

44

ða þ bÞl

p1⁄22h0  sinð2h0Þðr0

31Þ

2

h01⁄21 þ C  R2ð1  CÞ

þ

alþ bC

55

ða þ bÞl

p1⁄22h0 þ sinð2h0Þðr0

32Þ

2

h01⁄21 þ C  R2ð1  CÞ " #

;

ð7Þ

which immediately leads to the following two coupled nonlinear equations for C

44 and C55 by observing the fact that r0

31 and

r0

32 are arbitrary

l

C

55

1⁄4 1 þ

alþ bC

44

ða þ bÞl

pg1⁄22h0  sinð2h0Þ

2h01⁄21 þ C  R2

ð1  CÞ

;

l

C

44

1⁄4 1 þ

alþ bC

55

ða þ bÞl

pg1⁄22h0 þ sinð2h0Þ

2h01⁄21 þ C  R2

ð1  CÞ

;

ð8Þ

where g 1⁄4 c2

pab 1⁄4 4

pðR2R2Þ

. The above set of nonlinear equations can be solved easily through iteration.

In the following we discuss in more detail some special cases for the above solution.

2.1.1. Aligned cracks (h0 = 0)

In the case of aligned cracks, it follows from Eq. (8) that

C

55 1⁄4 l;

1

C2 1⁄4 1 þ

aC þ b

ða þ bÞC

2pg

1 þ C  R2ð1  CÞ

; ð9Þ

which means that there is no reduction in stiffness along the x1-direction for aligned cracks.

When the crack density parameter is extremely low, i.e., g  1, Eq. (9)2 further reduces to

C

44 1⁄4 l=ð1 þ pgÞ; ð10Þ

which is exactly the non-interaction approximation (NIA) for aligned cracks described by Kachanov [14].

Fig. 2 illustrates C

44 calculated by using Eq. (9). The dashed line is the result of NIA. We see that the predictions of the

GSCM are lower than that of NIA, especially when g is large.

2.1.2. Randomly oriented cracks (h0 = p/2)

In the case of randomly oriented cracks, it follows from Eq. (8) that

C

44 1⁄4 C55 1⁄4 l

;

1

C 1⁄4 1 þ

pg

1 þ C  R2ð1  CÞ

; ð11Þ

x1

x2

*

55

*

44 ,CC

μ

x

y

a b

c

Fig. 1. The generalized self-consistent method.

1912 X. Wang et al. / Engineering Fracture Mechanics 76 (2009) 1910–1919
through which C can be exactly determined as

C 1⁄4 ðpg  2R2

Þ þ

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

p2g2 þ 4ð1  pgR2Þ

q

2ð1 þ R2

Þ : ð12Þ

Eq. (11)1 shows that the effective material properties are still isotropic when the cracks are randomly oriented. In addition

when the crack density parameter is extremely low, i.e., g  1, Eq. (11)2 further reduces to

l1⁄4 l=ð1 þ 0:5pgÞ; ð13Þ

which is again, the NIA described by Kachanov, this time for randomly oriented cracks [14].

Fig. 3 shows l* calculated by using Eq. (12). The dashed line is the result of NIA. We observe a similar behavior here to that

of C

44 for aligned cracks. We also observe that the values of l* for a solid with randomly oriented cracks are higher than those

of C

44 for a solid with aligned cracks.

2.1.3. An arbitrary degree of alignment (0 < h0 < p/2)

The present model can be used conveniently to predict the anisotropic effective moduli of a solid containing microcracks

with an arbitrary degree of alignment through the introduction of the crack orientation distribution function /(h).

We demonstrate in Figs. 4 and 5 the predicted values of C

44 and C55 as functions of g for five different values of h0 = p/16,

p/6, p/4, p/3, p/2. It is clear from the two figures that a decrease in h0 will result in a decrease in C

44 and an increase in C55.

Fig. 2. The calculated C

44 for an isotropic solid with perfectly aligned cracks.

Fig. 3. The calculated l1⁄4 C

44 1⁄4 C55 for an isotropic solid with randomly oriented cracks.

X. Wang et al. / Engineering Fracture Mechanics 76 (2009) 1910–1919 1913
The anisotropy of the effective medium monotonically increases as h0 decreases from p/2 (for randomly oriented cracks) to

zero (for perfectly aligned cracks).

Due to the fact that the parameter l1⁄4 ffiffiffiffiffiffiffiffiffiffiffiffiffiffi

C

44C55 p is an invariant under coordinate rotation, it can be considered as a measure of the overall stiffness of the cracked material. Therefore, it is of interest to check how l* is influenced by h0 and g. We

demonstrate in Fig. 6 the predicted values of l* as a function of g for four different values of h0. It is observed that an increase

in h0 will cause a decrease in the overall stiffness characterized by l*

, especially when the crack density g is large. This fact

can be more clearly observed in Fig. 7 in which l* is plotted as a function of h0 for six different values of g. Figs. 6 and 7 demonstrate that for a microcracked solid, randomly oriented cracks (h0 = p/2) gives the lowest stiffness among all possible crack

orientation distributions, while aligned cracks (h0 = 0) gives the greatest stiffness.

2.2. Anisotropic undamaged material

In the previous discussion we have assumed that the undamaged material is isotropic. Next we consider the more complex situation in which the undamaged material is anisotropic. Without losing generality we can assume that the intact

material is orthotropic with its principal directions along the x1 and x2 axes, respectively. When the intact material is orthotropic, we cannot necessarily take the crack orientation distribution function /(h) as an even function of h. Here /(h) takes

the following form

/ðhÞ 1⁄4

1

h2h1 ; h1  h  h2

0; h > h2 or h < h1

(

ð14Þ

Fig. 4. The predicted values of C

44 for a cracked isotropic solid as a function of g for five different values of h0.

Fig. 5. The predicted values of C

55 for a cracked isotropic solid as a function of g for five different values of h0.

1914 X. Wang et al. / Engineering Fracture Mechanics 76 (2009) 1910–1919
where |h1|, |h2| 6 p and h2  h1 6 p. The definition region of h in Eq. (14) has been enlarged relative to Eq. (2), to |h| 6 p to

incorporate the complex situation in which the average orientation of the cracks h = (h1 + h2)/2 may be arbitrary with respect

to the x1, x2 axes. In this case, perfectly aligned cracks and randomly oriented cracks correspond to h1 = h2 and h2 = h1 = p/2

(or more generally but equivalently h2h1 = p), respectively, in Eq. (14).

Now we can write the energy relationship in Eq. (3) as follows

C

44ðr0

31Þ

2

2l2  C

45r0

31r032

l2 þ

C

55ðr0

32Þ

2

2l2 1⁄4 ðr0

31Þ

2

2C55

þ ðr0

32Þ

2

2C44

þ g

2c2

Z h2

h1

/ðhÞ

Z

Ck

1⁄2u3t

0

3dCkdh; ð15Þ

where r0

31 and r032 are the anti-plane, far-field stresses, C

44, C55 and C45 are the three relevant effective moduli of the damaged material, l1⁄4

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

C

44C55  C245 q

> 0 and C44 and C55 are the two material moduli of the undamaged material. When the

undamaged material is anisotropic, it may very well be impossible to calculate exact values for [u3] using the boundary-value

problem shown in Fig. 1. Therefore, in the following we use an approximation to the BVP to calculate [u3]. When the crack

density parameter is relatively low, we have R  1 or equivalently a  b  c. As a result, [u3] can be obtained approximately

through the following method: (i) first calculate the uniform stress field within an uncracked circular inclusion with material

moduli C44 and C55 surrounded by the effective medium with material moduli C

44, C55 and C45 [15,16]; (ii) then solve the

problem of a crack in an infinite homogeneous material with material moduli C44 and C55 subjected to the remote loading

which is equal to the internal uniform stress field obtained in (i). By using this method we obtain the following approximate

expression of [u3]

Fig. 6. The predicted values of l* for a cracked isotropic solid as a function of g for four different values of h0.

Fig. 7. The predicted values of l* for a cracked isotropic solid as a function of h0 for six different values of g.

X. Wang et al. / Engineering Fracture Mechanics 76 (2009) 1910–1919 1915
1⁄2u3 1⁄4 2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

c2  x2 p

l

2C44

l þ C44

ðlþ C

55Þr0

32  C

45r0

31

l1⁄21 þ C þ qð1  CÞ cos h  2l

l þ C44

ðlþ C

44Þr0

31  C

45r0

32

l1⁄21 þ C  qð1  CÞ sin h

 

ðjxj  cÞ ð16Þ

where q 1⁄4

ffiffiffiffiffi

C44 p 

ffiffiffiffiffi

C55 p

ffiffiffiffiffi

C44 p þ

ffiffiffiffiffi

C55 p , l 1⁄4 ffiffiffiffiffiffiffiffiffiffiffiffiffiffi

C44C55 p and C = l*/l.

The integral in Eq. (15) can now be carried out to arrive at the following relationship

C

44ðr0

31Þ

2

2l2  C

45r0

31r032

l2 þ

C

55ðr0

32Þ

2

2l2 1⁄4 ðr0

31Þ

2

2C55

þ ðr0

32Þ

2

2C44

þ

pg

2lðh2  h1Þ

	 C44

2ðl þ C44Þ

ðlþ C

55Þðr0

32Þ

2  C

45r0

31r032

l1⁄21 þ C þ qð1  CÞ 1⁄22ðh2  h1Þ þ sinð2h2Þ  sinð2h1Þ "

þ l

2ðl þ C44Þ

ðlþ C

44Þðr0

31Þ

2  C

45r0

31r032

l1⁄21 þ C  qð1  CÞ 1⁄22ðh2  h1Þ  sinð2h2Þ þ sinð2h1Þ

 l

l þ C44

ðlþ C

44Þr0

31r032  C

45ðr0

32Þ

2

l1⁄21 þ C  qð1  CÞ þ

C44

l þ C44

ðlþ C

55Þr0

31r032  C

45ðr0

31Þ

2

l1⁄21 þ C þ qð1  CÞ " #

1⁄2sin2 h2  sin2 h1

#

: ð17Þ

In view of the fact that r0

31 and r032 are arbitrary, we arrive at the following three coupled, nonlinear equations for C

44, C55

and C

45

C

44

l2 1⁄4 1

C55

þ

pg

lðh2  h1Þðl þ C44Þ

ðlþ C

44Þ1⁄22ðh2  h1Þ  sinð2h2Þ þ sinð2h1Þ

21⁄21 þ C  qð1  CÞ þ

C44C

451⁄2sin2 h2  sin2 h1

l1⁄21 þ C þ qð1  CÞ " #;

C

55

l2 1⁄4 1

C44

þ

pg

lðh2  h1Þðl þ C44Þ

C44ðlþ C

55Þ1⁄22ðh2  h1Þ þ sinð2h2Þ  sinð2h1Þ

2l1⁄21 þ C þ qð1  CÞ þ

C

451⁄2sin2 h2  sin2 h1

1 þ C  qð1  CÞ

" #;

C

45

l2 1⁄4 pg2lðh2  h1Þðl þ C44Þ

C44C

451⁄22ðh2  h1Þ þ sinð2h2Þ  sinð2h1Þ

2l1⁄21 þ C þ qð1  CÞ þ

C

451⁄22ðh2  h1Þ  sinð2h2Þ þ sinð2h1Þ

21⁄21 þ C  qð1  CÞ 

þ lþ C

44

1 þ C  qð1  CÞ

þ

C44ðlþ C

55Þ

l1⁄21 þ C þ qð1  CÞ  1⁄2sin2 h2  sin2 h1



;

ð18Þ

which can also be solved easily through iteration. We see from Eq. (18) that when /(h), defined by Eq. (14), is an even function of h, i.e., h1 = h2, then the effective medium will be orthotropic with C

45 1⁄4 0. In other words, when the average orientation of the cracks is zero the undamaged material and the effective medium will possess the same principal axes. In the

following we concentrate our discussion on the effective elastic moduli of a solid containing perfectly aligned cracks and

one containing randomly oriented cracks, respectively.

2.2.1. Aligned cracks (h1 = h2)

When the cracks are perfectly aligned, we have h1 = h2 and both the numerator and denominator in the right hand side of

Eq. (18) vanish. In this case applying the L’Hospital’s rule to Eq. (18) when h2 ? h1 yields the following:

C

44

l2 1⁄4 1

C55

þ

pg

lðl þ C44Þ

ðlþ C

44Þ1⁄21  cosð2h1Þ

1 þ C  qð1  CÞ þ

C44C

45 sinð2h1Þ

l1⁄21 þ C þ qð1  CÞ  ;

C

55

l2 1⁄4 1

C44

þ

pg

lðl þ C44Þ

C44ðlþ C

55Þ1⁄21 þ cosð2h1Þ

l1⁄21 þ C þ qð1  CÞ þ

C

45 sinð2h1Þ

1 þ C  qð1  CÞ

 ;

C

45

l2 1⁄4 pg2lðl þ C44Þ

C44C

451⁄21 þ cosð2h1Þ

l1⁄21 þ C þ qð1  CÞ þ

C

451⁄21  cosð2h1Þ

1 þ C  qð1  CÞ

þ lþ C

44

1 þ C  qð1  CÞ

þ

C44ðlþ C

55Þ

l1⁄21 þ C þ qð1  CÞ   sinð2h1Þ

 

:

ð19Þ

Now, we can further look into two special cases of the above formulation: the cracks are aligned horizontally (h1 = h2 = 0)

and the cracks are aligned vertically (h1 = h2 = p/2). When h1 = h2 = 0, it follows from Eq. (19) that

C

45 1⁄4 0; C55 1⁄4 C55;

1

C

44

1⁄4 1

C44

þ

2pgC44ðlþ C55Þ

llðl þ C44Þ1⁄21 þ C þ qð1  CÞ ; ð20Þ

which means that there is no reduction in stiffness along the x1-direction for horizontally aligned cracks.

On the other hand when h1 = h2 = p/2, it follows from Eq. (19) that

C

45 1⁄4 0; C44 1⁄4 C44;

1

C

55

1⁄4 1

C55

þ

2pgðlþ C44Þ

lðl þ C44Þ1⁄21 þ C  qð1  CÞ ; ð21Þ

which means that there is no reduction in stiffness along the x2-direction for vertically aligned cracks.

1916 X. Wang et al. / Engineering Fracture Mechanics 76 (2009) 1910–1919
Fig. 8 illustrates the variations of C

44 for h1 = h2 = 0 and C55 for h1 = h2 = p/2 as functions of g. To plot the results on one

graph, the undamaged anisotropic moduli are specifically chosen as C44 = 1.2l and C55 = 0.8333l. It is interesting to observe

that C

44 for h1 = h2 = 0 and C55 for h1 = h2 = p/2, both of which are decreasing functions of g, converge to basically the same

value as g ? 1. We also observe that the effective material is isotropic ðC

44 1⁄4 C55 1⁄4 0:8333lÞ at low crack density g  0.10 for

horizontally aligned cracks (h1 = h2 = 0).

We can also see from Eq. (19) that when the cracks are not aligned horizontally or vertically, the principal directions of

the anisotropic effective medium are no longer the same as those of the undamaged material. We illustrate in Fig. 9, the variations of C

44, C55 and C45 as functions of h1 (=h2) for an orthotropic solid with undamaged moduli, C44 = 1.2l and

C55 = 0.8333l containing aligned cracks with g = 1/p. We observe from Fig. 9 that: (i) C

45 attains a maximum value of

C

45 1⁄4 0:3073l at h1 = 48 (the reason why C45 does not attain its maximum value at 45 is due to the anisotropy of the

undamaged material); (ii) C

44 is an increasing function and C55 is a decreasing function of h1. As expected

C

44 1⁄4 C44 1⁄4 1:2l for vertically aligned cracks (h1 = 90) and C55 1⁄4 C55 1⁄4 0:8l for horizontally aligned cracks (h1 = 0). More

interestingly, our numerical results also verify that the effective compliance constants eS

45 1⁄4 e

C

45=l2 and eS55 1⁄4 e

C

44=l2

in the Cartesian coordinate system ðx0

1; x02Þ, which is rotated by the crack alignment angle h1 with respect to the principal

coordinate system (x1, x2), are exactly the same as the compliance constants eS45 1⁄4 e

C 45=l2 and eS55 1⁄4 eC 44=l2 for the undamaged material in the same coordinate system ðx0

1; x02Þ .

Fig. 8. The variations of C

44 (h1 = h2 = 0) and C55 (h1 = h2 = p/2) for an orthotropic solid C44 = 1.2l and C55 = 0.8333l containing perfectly aligned cracks as

functions of g.

Fig. 9. The variations of C

44, C55 and C45 as functions of h1 (=h2) for an orthotropic solid C44 = 1.2l and C55 = 0.8333l containing perfectly aligned cracks with

g=1/p.

X. Wang et al. / Engineering Fracture Mechanics 76 (2009) 1910–1919 1917
2.2.2. Randomly oriented cracks (h2 = h1 = p/2)

When the cracks are randomly oriented, Eq. (18) reduces to

C

45 1⁄4 0;

1

C

55

1⁄4 1

C55

þ

pgðlþ C

44Þ

lðl þ C44Þ1⁄21 þ C  qð1  CÞ ;

1

C

44

1⁄4 1

C44

þ

pgC44ðlþ C

55Þ

llðl þ C44Þ1⁄21 þ C þ qð1  CÞ :

ð22Þ

From the above expressions, we can see that C

44 – C55 when C44 – C55 even when the cracks are randomly oriented. To

illustrate this fact more clearly we present in Fig. 10, the variations of C

44 and C55 for an orthotropic solid containing randomly oriented cracks as functions of g. The undamaged anisotropic moduli are chosen as C44 = 1.2l and C55 = 0.8333l. It

is clear that C

44 > C55, and in addition the difference in C44 and C55 decreases as the crack density increases.

Finally it is of interest to point out that when the undamaged material is isotropic, our calculations show that the result by

using the approximation of the GSCM always lies between that by using the exact GSCM and that by using NIA. In view of the

fact that the GSCM is only an approximate scheme to account for crack interactions, the predictions derived by using the

GSCM are not necessarily more accurate than those derived by using an approximation to the GSCM.

3. Conclusions

By means of the GSCM, we analytically investigated the anisotropic effective moduli of a cracked solid subjected to antiplane shear deformation. When the undamaged solid is isotropic, the two coupled, nonlinear equations for the two effective

moduli C

44 and C55 were obtained and shown in Eq. (8). When the undamaged solid is anisotropic, the three coupled nonlinear equations for the three effective moduli C

44, C55 and C45 were derived and shown in Eq. (18). Detailed numerical results were presented to illustrate how the anisotropic effective moduli are influenced by various degrees of crack alignment

(characterized by h0 for the isotropic undamaged material or h1 and h2 for the anisotropic undamaged material) and crack

density (characterized by g).

Acknowledgements

The three referees’ comments and suggestions on revising the initial manuscript are highly appreciated. This research was

supported by the United States Army Research Laboratory through the Composite Materials Technology cooperative agreement with the Center for Composite Materials at the University of Delaware.

Appendix A. Detailed derivation of Eq. (6)

It is convenient to solve the boundary-value problem in the Cartesian coordinate system (x, y, z) with z = x3, as shown in

Fig. 1. First the remote uniform loading r0

zx and r0zy can be written as

Fig. 10. The variations of C

44 and C55 for an orthotropic solid C44 = 1.2l and C55 = 0.8333l containing randomly oriented cracks as functions of g.

1918 X. Wang et al. / Engineering Fracture Mechanics 76 (2009) 1910–1919
r0

zx 1⁄4 r031 cos h þ r032 sin h;

r0

zy 1⁄4 r032 cos h  r031 sin h: ðA1Þ

Secondly the effective moduli of the damaged material in the (x, y, z) coordinate system, eC

44; e

C

55 and e

C

45 can be derived

as

e

C

55 1⁄4 C

55 cos2 h þ C44 sin2 h;

e

C

44 1⁄4 C

44 cos2 h þ C55 sin2 h;

e

C

45 1⁄4 ðC

44  C55Þ cos h sin h:

ðA2Þ

Due to the fact that the cracked elliptical inclusion is isotropic, if one is only concerned with the field within the inclusion,

it is sufficient to treat the matrix as isotropic with shear modulus l* subject to virtual remote uniform shear stresses r~0

zx and

r~ 0

zy such that [13]

r~0

zx  ir~0zy 1⁄4 ða  ipbÞðr0

zx þ pr0zyÞ

ða þ bÞImfpg ; ðA3Þ

where the complex constant p is given by

p 1⁄4 e

C

45 þ il

e

C

44

1⁄4 ðC

55  C44Þ cos h sin h þ i ffiffiffiffiffiffiffiffiffiffiffiffiffiffi

C

44C55 p

C

44 cos2 h þ C55 sin2 h : ðA4Þ

Here we are only interested in the expression for r~0

zy since the loading r~0zx will not induce any crack opening displacement.

It follows from Eqs. (A1), (A3), and (A4) that

r~0

zy 1⁄4 bRefpgr0

31 þ 1⁄2aImfpg þ bjpj

2

r0

32

ða þ bÞImfpg cos h þ

bRefpgr0

32  1⁄2aImfpg þ bjpj

2

r0

31

ða þ bÞImfpg sin h

1⁄4 alþ bC

55

ða þ bÞlcos hr0

32  alþ bC

44

ða þ bÞlsin hr0

31: ðA5Þ

On the other hand it can be deduced from Ref. [12] that the crack opening displacement [u3] due to remote uniform loading r~0

zy can be expressed as

1⁄2u3 1⁄4

4r~0

zy

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

c2  x2 p

l1⁄21 þ C  R2ð1  CÞ ðjxj  cÞ: ðA6Þ

Starting with Eqs. (A5) and (A6), we can easily arrive at Eq. (6).

References

[1] Budiansky B, O’Connel RJ. Elastic moduli of a cracked solid. Int J Solid Struct 1976;12:81–97.

[2] Aboudi J, Benveniste Y. The effective moduli of cracked bodies in plane deformations. Engng Fract Mech 1987;26(2):171–84.

[3] Santare MH, Crocombe AD, Anlas G. Anisotropic effective moduli of materials with microcracks. Engng Fract Mech 1995;52(5):833–42.

[4] Benveniste Y. On the Mori-Tanaka’s method in cracked bodies. Mech Res Commun 1986;13:193–201.

[5] Salganik RL. Mechanics of bodies with many cracks. Mech Solid 1973;8:135–43 [English translation].

[6] Zimmerman RW. The effect of microcracks on the elastic moduli of brittle materials. J Mater Sci-Lett 1985;4:1457–60.

[7] Hashin Z. The differential scheme and its application to cracked materials. J Mech Phys Solid 1988;13:193–201.

[8] Sayers C, Kachanov M. A simple technique for finding effective elastic constants of cracked solids for arbitrary crack orientation statistics. Int J Solid

Struct 1991;7:671–80.

[9] Su D, Santare MH, Gazonas GA. The effect of crack face contact on the anisotropic effective moduli of microcrack damaged media. Engng Fract Mech

2007;74:1436–55.

[10] Kushch VI, Sevostianov I, Mishnaevsky L. Effect of crack orientation statistics on effective stiffness of microcracked solid. Int J Solid Struct

2009;46:1574–88.

[11] Feltman RS, Santare MH. Anisotropic effective moduli of cracked short fiber reinforced composites. ASME J Appl Mech 1999;66:709–13.

[12] Wu CH, Chen CH. A crack in a confocal elliptic inhomogeneity embedded in an infinite medium. ASME J Appl Mech 1990;57:91–6.

[13] Wang X, Pan E. Antiplane shear deformations of an anisotropic elliptical inhomogeneity with imperfect or viscous interface. Z Angew Math Mech

(ZAMM) 2008;88:142–50.

[14] Kachanov M. Effective elastic properties of cracked solids: critical review of some basic concept. Appl Mech Rev 1992;45(8):304–35.

[15] Ting TCT. Anisotropic elasticity—theory and applications. New York: Oxford University Press; 1996.

[16] Kattis MA, Providas E. Two-phase potentials in anisotropic elasticity: antiplane deformation. Int J Engng Sci 1998;36:801–11.

X. Wang et al. / Engineering Fracture Mechanics 76 (2009) 1910–1919 1919
Int J Fract (2009) 158:27–40

DOI 10.1007/s10704-009-9363-2

ORIGINAL PAPER

On the effective electroelastic properties of microcracked

generally anisotropic solids

Xu Wang · George A. Gazonas ·

Michael H. Santare

Received: 9 March 2009 / Accepted: 20 May 2009 / Published online: 9 June 2009

© Springer Science+Business Media B.V. 2009

Abstract In this study we first obtain the explicit

expressions for the 15 effective reduced elastic compliances of an elastically anisotropic solid containing multiple microcracks with an arbitrary degree of alignment

under two-dimensional deformations within the framework of the non-interaction approximation (NIA).

Under special situations, our results can reduce to the

classical ones derived by Bristow (J Appl Phys 11:

81–85, 1960), and Mauge and Kachanov (J Mech Phys

Solids 42(4):561–584, 1994). Some interesting phenomena are also observed. For example, when the

undamaged solid is orthotropic, the effective in-plane

shear modulus is dependent on the degree of the crack

alignment. The NIA method is then extended to obtain

the effective electroelastic properties of an anisotropic

piezoelectric solid containing two-dimensional insulating, permeable or conducting microcracks with an

arbitrary degree of alignment. We also derive a set of

fifteen coupled nonlinear equations for the unknown

effective reduced elastic compliances of a microcracked, anisotropic, elastic solid by using the generalized

X. Wang · M. H. Santare

Department of Mechanical Engineering and Center for

Composite Materials, University of Delaware, Newark,

DE 19716, USA

G. A. Gazonas (B)

US Army Research Laboratory, Weapons and Materials

Research Directorate, Aberdeen Proving Ground,

Aberdeen, MD, 21005, USA

e-mail: gazonas@arl.army.mil

self-consistent method (GSCM). The set of coupled

nonlinear equations can be solved through iteration.

Keywords Non-interaction approximation ·

Effective properties · Microcracks · Elastic material ·

Piezoelectric material

1 Introduction

The prediction of the effective elastic properties of

anisotropic materials containing microcracks has been

a topic of micromechanics for nearly 30 years (see for

example, Gottesman et al. 1980; Hashin 1988; Mauge

and Kachanov 1994; Feltman and Santare 1999; among

others). The limitations of the previous works lie in the

following: (i) the two-dimensional, anisotropic matrix

studied so far is confined to the special case in which

the in-plane deformations and the out-of-plane deformations are decoupled (see for example, Mauge and

Kachanov 1994); (ii) in most of the studies the microcracks were assumed to be either perfectly aligned or

randomly oriented (Mauge and Kachanov 1994).

In this study we first address the problem of the

effective elastic properties of a generally anisotropic

elastic matrix containing microcracks with an arbitrary

degree of alignment under two-dimensional deformations. Here we consider any degree of crack alignment,

from perfectly aligned to completely random, through

the introduction of a crack orientation distribution function (Santare et al. 1995). Very concise expressions for

123
28 X. Wang et al.

the effective reduced elastic compliances of the cracked

material are obtained by making use of energy considerations and by ignoring crack interactions. During

the development, we adopt some results (for example,

the crack opening displacement) based on the Stroh

formalism for two-dimensional anisotropic elasticity

(Ting 1996). The current result based on a non-interaction approximation (NIA) still retains accuracy at substantially higher crack densities due to the cancellation of the competing interaction effects of shielding

and antishielding (Kachanov 1992; Mauge and Kachanov 1994; Kachanov 2007). The NIA method is then

extended to study the effective electroelastic properties of a microcracked, anisotropic, piezoelectric solid.

Here the microcracks can be insulating (Pak 1990; Suo

et al. 1992), permeable (Suo et al. 1992) or conducting

(McMeeking 1987; Suo 1993). It is of interest to point

out that the non-convex electric enthalpy used for insulating cracks is different than the convex energy density function used for conducting cracks. The effective

elastic moduli of a microcracked, generally anisotropic

elastic solid are also discussed within the framework of

the generalized self-consistent method (GSCM), which

approximately takes into account crack interaction

(Aboudi and Benveniste 1987; Santare et al. 1995).

It is found that some observations made by using the

NIA assumption are no longer valid when using the

GSCM.

2 The theory within the framework of NIA

The strain energy balance equation for a microcracked,

elastically anisotropic solid can be expressed as

(Eshelby 1956; Benveniste 1985; Aboudi and Benveniste 1987)

1

2

S∗

ijklσ0klσ0i j=1

2

Sijklσ0

klσ0i j+

1

2V



M

k=1



Ck

[ui]t

0

i dCk ,

(1)

where S∗

ijkl is the effective compliance of the damaged material, Sijkl is the compliance of the undamaged

material, V is the sample volume, M is the number of

microcracks within the solid, σ0

i j is the applied homogeneous stress field while t0

i = σ0i jn j is the traction

along the crack face if the crack did not exist and [ui]

is the crack opening displacement.

x1

x2

θ

x

y

a

−a

Fig. 1 A crack of half crack length a oriented at an angle θ to

the positive x1-axis

In this study we assume that all the cracks penetrate the solid through the x3-axis. The degree of crack

alignment can be described by the crack orientation

distribution function φ(θ) with θ,(|θ| < π/2) being

the angle between an individual crack and the positive

x1-axis as shown in Fig. 1 (Santare et al. 1995). Without

losing generality, φ(θ) can be taken as an even function of θ such that the average orientation is parallel

to the x1-axis. For simplicity, φ(θ) is specifically given

by (Santare et al. 1995; Wang et al. 2009)

φ(θ ) =

 1

2θ0 , |θ| ≤ θ0

0, |θ| > θ0

(2)

where θ0 ≤ π/2. The two special cases of perfectly

aligned cracks and randomly oriented cracks correspond to θ0 = 0 and θ0 = π/2, respectively in Eq. 2.

As a result, under two-dimensional deformations in

which the displacements depend on x1 and x2 only,

Eq. 1 can be further expressed as (Santare et al. 1995;

Ting 1996)

1

2

(σ0)

T S∗

σ0 = 1

2

(σ0)

T S

σ0

+

M

2A



θ0

−θ0

φ(θ ) 

Ck

(t

0)T [u]dCkdθ ,

(3)

where A is the sample area, and

σ0 =



σ0

11 σ022 σ023 σ013 σ012T

, (4)

123
On the effective electroelastic properties 29

S =

⎡

⎢

⎢

⎢

⎢

⎣

S

11 S12 S14 S

15 S

16

S

12 S22 S24 S

25 S

26

S

14 S24 S44 S

45 S

46

S

15 S25 S45 S55 S56

S

16 S26 S46 S

56 S

66

⎤

⎥

⎥

⎥

⎥

⎦

, (5)

S∗ =

⎡

⎢

⎢

⎢

⎢

⎣

S

∗

11 S

∗

12 S

∗

14 S

∗

15 S

∗

16

S

∗

12 S

∗

22 S

∗

24 S

∗

25 S

∗

26

S

∗

14 S

∗

24 S

∗

44 S

∗

45 S

∗

46

S

∗

15 S

∗

25 S

∗

45 S

∗

55 S

∗

56

S

∗

16 S

∗

26 S

∗

46 S

∗

56 S

∗

66

⎤

⎥

⎥

⎥

⎥

⎦

, (6)

t

0 =



t

0

1 t

0

2 t

0

3

T

, [u] = [[u1] [u2] [u3]]

T . (7)

We add that S

αβ = Sαβ−Sα3 S3β/S33 in Eq. 5 are the

reduced elastic compliances of the undamaged solid,

and S

∗

αβ = S∗

αβ − S∗

α3 S∗3β/S∗33 in Eq. 6 are the effective

reduced elastic compliances of the cracked solid (Ting

1996).

The traction vector t

0 is given by

t

0 =

⎡

⎢

⎣

σ0

12 cos θ − σ011 sin θ

σ0

22 cos θ − σ012 sin θ

σ0

23 cos θ − σ013 sin θ

⎤

⎥

⎦ . (8)

If we ignore the crack interactions, i.e., η 1, then the

crack opening displacement vector [u] = [[u1] [u2]

[u3]]

T can be simply given by (Ting 1996)

[u] = 2



a2 − x2L−1t

0, (|x| ≤ a) (9)

where a is the half crack length, the 3 × 3 real and

symmetric Barnett–Lothe tensor L for the undamaged

solid is positive definite and L−1 is further written in

the following form for the convenience of the following

theoretical development

L−1 =

⎡

⎣

Y11 Y12 Y13

Y12 Y22 Y23

Y13 Y23 Y33

⎤

⎦ . (10)

We stress here that in writing Eq. 9 we have employed

the property of the second-order tensor L under coordinate rotation (also see Fig. 114 on page 421 in Ting

1996).

Substituting Eqs. 2, 8 and 9 into Eq. 3, we can arrive

at

1

2

(σ0)

T S∗

σ0 = 1

2

(σ0)

T S

σ

0 + η

2a2θ0

⎡

⎣

a

−a



a2 − x2dx

⎤

⎦

⎡

⎢

⎣



θ0

−θ0

(t

0)TL−1t0dθ

⎤

⎥

⎦ , (11)

where η = Ma2/A is the crack density parameter.

We can clearly observe from the above expression

that S∗ is also positive definite due to the fact that

the value of the second term on the right hand side of

Eq. 11 is always positive (i.e., the elastic energy density

is always increased through the introduction of the microcracks). The two integrals in Eq. 11 can be exactly

carried out as

a

−a



a2 − x2dx = πa2

2 , (12)

and



θ0

−θ0

(t

0)TL−1t0dθ = Y11

2 [2θ0 − sin(2θ0)](σ0

11)

2

+

Y22

2 [2θ0 + sin(2θ0)](σ0

22)

2

+

Y33

2 [2θ0 + sin(2θ0)](σ0

23)

2

+

Y33

2 [2θ0 − sin(2θ0)](σ0

13)

2

+

1

2 {Y11[2θ0 + sin(2θ0)]

+Y22[2θ0 − sin(2θ0)]}(σ0

12)

2

+Y13[2θ0 − sin(2θ0)]σ0

11σ013

+Y12[2θ0 − sin(2θ0)]σ0

11σ012

+Y23[2θ0 + sin(2θ0)]σ0

22σ023

+Y12[2θ0 + sin(2θ0)]σ0

22σ012

+Y13[2θ0 + sin(2θ0)]σ0

23σ012

+Y23[2θ0 − sin(2θ0)]σ0

13σ012.

(13)

In view of the fact that the applied homogeneous

stresses σ0

11, σ022, σ023, σ013 and σ012 can be arbitrary,

then we obtain the following expressions for the effective reduced elastic compliances

S

∗

11 = S

11 + πηY11

4θ0

[2θ0 − sin(2θ0)],

S

∗

22 = S

22 + πηY22

4θ0

[2θ0 + sin(2θ0)],

S

∗

44 = S

44 + πηY33

4θ0

[2θ0 + sin(2θ0)],

S

∗

55 = S

55 + πηY33

4θ0

[2θ0 − sin(2θ0)],

S

∗

66 = S

66 + πη

4θ0

{Y11[2θ0 + sin(2θ0)]

+Y22[2θ0 − sin(2θ0)]},

123
30 X. Wang et al.

S

∗

12 = S

12,

S

∗

14 = S

14,

S

∗

15 = S

15 + πηY13

4θ0

[2θ0 − sin(2θ0)],

S

∗

16 = S

16 + πηY12

4θ0

[2θ0 − sin(2θ0)],

S

∗

24 = S

24 + πηY23

4θ0

[2θ0 + sin(2θ0)],

S

∗

25 = S

25,

S

∗

26 = S

26 + πηY12

4θ0

[2θ0 + sin(2θ0)],

S

∗

45 = S

45,

S

∗

46 = S

46 + πηY13

4θ0

[2θ0 + sin(2θ0)],

S

∗

56 = S

56 + πηY23

4θ0

[2θ0 − sin(2θ0)]. (14)

Equation 14 indicates that the four compliance components S

∗

12, S

∗

14, S

∗

25 and S

∗

45 are equal to the corresponding ones for the undamaged material and that S

∗

ii ≥

S

ii, (i = 1, 2, 4, 5, 6) since Y11, Y22, Y33 > 0. It is of

interest to notice that the present result of S

∗

12 = S

12 is

consistent with the observation by Mauge and Kachanov (1994) that ν∗

12/E∗1 = ν12/E1 for an orthotropic

solid containing perfectly aligned or randomly oriented

cracks under plane stress deformation. The result we

obtain is valid for generally anisotropic materials containing cracks with an arbitrary degree of alignment.

When the undamaged material is orthotropic, it

follows from the above expression that

1

μ∗

12

= 1

μ12

+ πη

4θ0

×




2θ0 + sin (2θ0)

L11

+

2θ0 − sin (2θ0)

L22 

,

(15)

which clearly demonstrates that the effective in-plane

shear modulus μ∗

12 depends on the degree of crack

alignment characterized by θ0 when L11 = L22. Furthermore when the undamaged material is isotropic

such that L11 = L22 = μ

1−ν , with μ and ν being

the shear modulus and the Poisson’s ratio, respectively,

then Eq. 15 reduces to

μ∗

12 = μ

1 + πη(1 − ν), (16)

which is the classical NIA result obtained by Bristow

(1960), and is independent of the degree of crack align0 10 20 30 40 50 60 70 80 90 0.76

0.77

0.78

0.79

0.8

0.81

0.82

0.83

0.84

0.85

θ (Deg) 0

πημ /L =0.25 12 11

L /L =1.5

L /L =1.0

L11 /L22 =0.5

11 22

11 22

μ* /μ 12 12

Fig. 2 The variations of μ∗

12 of a cracked orthotropic elastic

material as a function of θ0 by using Eq. 15 for three different

values of L11/L22 = 0.5, 1.0, 1.5 with πημ12/L11 = 0.25

ment (Kachanov 1992). We illustrate in Fig. 2 the variation of μ∗

12 as a function of θ0 by using Eq. 15 for

three different values of L11/L22 = 0.5, 1.0, 1.5 with

πημ12/L11 = 0.25. We observe from Fig. 2 that μ∗

12 is

an increasing function of θ0 when L11 < L22; whereas

it is a decreasing function of θ0 when L11 > L22.

In the following we discuss two typical cases for the

crack orientations:

2.1 Perfectly aligned cracks

When the microcracks are perfectly aligned, we have

θ0 = 0. Consequently it follows from Eq. 14 that

S

∗

11 = S

11,

S

∗

22 = S

22 + πηY22,

S

∗

44 = S

44 + πηY33,

S

∗

55 = S

55,

S

∗

66 = S

66 + πηY11,

S

∗

12 = S

12,

S

∗

14 = S

14,

S

∗

15 = S

15,

S

∗

16 = S

16,

S

∗

24 = S

24 + πηY23,

S

∗

25 = S

25,

S

∗

26 = S

26 + πηY12,

123
On the effective electroelastic properties 31

S

∗

45 = S

45,

S

∗

46 = S

46 + πηY13,

S

∗

56 = S

56. (17)

In addition, if the undamaged material is orthotropic,

Eq. 17 further reduces to Eq. 27 derived by Mauge and

Kachanov (1994) in view of the fact that Y11 and Y22

for orthotropic materials are explicitly given by (Suo

1990)

Y11 =



S

11 


2



S

11 S22 + 2S

12 + S661

2

,

(18)

Y22 =



S

22 


2



S

11 S22 + 2S

12 + S661

2

.

Here it should be noticed that Mauge and Kachanov

(1994) discussed the plane stress problem.

2.2 Randomly oriented cracks

When the microcracks are oriented randomly, we have

θ0 = π/2. Consequently it follows from Eq. 14 that

S

∗

11 = S

11 +

1

2

πηY11,

S

∗

22 = S

22 +

1

2

πηY22,

S

∗

44 = S

44 +

1

2

πηY33,

S

∗

55 = S

55 +

1

2

πηY33,

S

∗

66 = S

66 +

1

2

πη(Y11 + Y22),

S

∗

12 = S

12,

S

∗

14 = S

14,

S

∗

15 = S

15 +

1

2

πηY13,

S

∗

16 = S

16 +

1

2

πηY12,

S

∗

24 = S

24 +

1

2

πηY23,

S

∗

25 = S

25,

S

∗

26 = S

26 +

1

2

πηY12,

S

∗

45 = S

45,

S

∗

46 = S

46 +

1

2

πηY13,

S

∗

56 = S

56 +

1

2

πηY23. (19)

In addition if the undamaged material is orthotropic,

Eq. 19 reduces to Eq. 32 derived by Mauge and Kachanov (1994). It seems that the factor ‘π/4’ appearing in

the third expression of Eq. 32 in Mauge and Kachanov

(1994) should read ‘π/2’. Furthermore it follows from

Eq. 19 that the overall effective elastic properties will be

transversely isotropic for an originally isotropic material containing randomly oriented cracks (θ0 = π/2).

Once the effective reduced elastic compliances have

been obtained, the effective stiffnesses can be conveniently obtained as

C0∗ = (S∗

)

−1, (20)

where

C0∗ =

⎡

⎢

⎢

⎢

⎢

⎣

C∗

11 C∗12 C∗14 C∗

15 C∗

16

C∗

12 C∗22 C∗24 C∗

25 C∗

26

C∗

14 C∗24 C∗44 C∗

45 C∗

46

C∗

15 C∗25 C∗45 C∗55 C∗56

C∗

16 C∗26 C∗46 C∗

56 C∗

66

⎤

⎥

⎥

⎥

⎥

⎦

. (21)

3 Extension to piezoelectric solids

In the previous section we have derived the effective

reduced elastic compliances of a purely elastic solid

containing microcracks. In fact, the energy method can

also be conveniently adopted to derive the effective

electroelastic properties of a piezoelectric solid containing insulating, permeable or conducting microcracks.

3.1 A piezoelectric solid containing insulating

microcracks

In this case, the energy balance relationship for twodimensional problems, in which the displacements ui

and the electric potential φ depend on x1 and x2 only,

can be expressed in the following form when the piezoelectric solid containing insulating microcracks is subjected to the homogeneous stresses σ0

i j and the homogeneous electric displacements D0

i (Suo et al. 1992)

1

2

(σ ̃ 0)

T (C ̃ 0∗)−1σ ̃ 0 = 1

2

(σ ̃ 0)

T (C ̃ 0)−1σ ̃ 0

+

M

2A



θ0

−θ0

φ(θ )

Ck

( ̃t

0)T [u ̃]dCkdθ ,

(22)

123
32 X. Wang et al.

where

σ ̃ 0 = [σ0

11 σ022 σ023 σ013 σ012 D01 D02]

T , (23)

[u ̃] = [[u1] [u2] [u3] [φ]]

T , (24)

C ̃ 0 =

⎡

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎣

C11 C12 C14 C15 C16 e11 e21

C12 C22 C24 C25 C26 e12 e22

C14 C24 C44 C45 C46 e14 e24

C15 C25 C45 C55 C56 e15 e25

C16 C26 C46 C56 C66 e16 e26

e11 e12 e14 e15 e16 −
11 −
12

e21 e22 e24 e25 e26 −
12 −
22

⎤

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎦

,

(25)

C ̃ 0∗ =

⎡

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎣

C∗

11 C∗12 C∗14 C∗

15 C∗

16 e∗11 e∗21

C∗

12 C∗22 C∗24 C∗

25 C∗

26 e∗12 e∗22

C∗

14 C∗24 C∗44 C∗

45 C∗

46 e∗14 e∗24

C∗

15 C∗25 C∗45 C∗55 C∗56 e∗15 e∗25

C∗

16 C∗26 C∗46 C∗

56 C∗

66 e∗16 e∗26

e∗

11 e∗12 e∗14 e∗

15 e∗

16 −
∗11 −
∗12

e∗

21 e∗22 e∗24 e∗

25 e∗

26 −
∗12 −
∗22

⎤

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎦

,

(26)

 ̃t

0 =

⎡

⎢

⎢

⎢

⎣

σ0

12 cos θ − σ011 sin θ

σ0

22 cos θ − σ012 sin θ

σ0

23 cos θ − σ013 sin θ

D0

2 cos θ − D01 sin θ

⎤

⎥

⎥

⎥

⎦ . (27)

In Eqs. 25 and 26, Ci j, ei j and 
i j are the elastic, piezoelectric and dielectric constants of the undamaged

solid, while C∗

i j, e∗i j and 
∗i j are the corresponding

(unknown) constants for the damaged solid. The energy

expression in Eq. 22 is in fact the electric enthalpy

(Suo et al. 1992). When ignoring crack interactions

and assuming that the crack surfaces are traction-free

and charge-free (Pak 1990; Suo et al. 1992), the jump

in displacement and electric potential across the crack

surface can be simply given by

[u ̃] = 2



a2 − x2L ̃ −1 ̃t

0, (|x| ≤ a) (28)

where L ̃ is the 4 × 4 real and symmetric Barnett–

Lothe tensor for the undamaged piezoelectric solid

(Ting 1996). We add that L ̃ is not positive definite and

can be further expressed as

L ̃ =




L11 L14

LT

14 −L44 

, (29)

where the 3 × 3 symmetric matrix L11 is positive definite and L44 > 0.

Substituting the above results into Eq. 22, we finally

arrive at a concise expression for the effective electroelastic properties of the microcracked piezoelectric

solid

(C ̃ 0∗)

−1 = (C ̃ 0)−1 + πη[2θ0 − sin(2θ0)]

4θ0

J ̃T

1 L ̃ −1J ̃1

+πη[2θ0 + sin(2θ0)]

4θ0

J ̃T

2 L ̃ −1J ̃2, (30)

where

J ̃1 =

⎡

⎢

⎢

⎣

1000000

0000100

0001000

0000010

⎤

⎥

⎥

⎦ ,

J ̃2 =

⎡

⎢

⎢

⎣

0000100

0100000

0010000

0000001

⎤

⎥

⎥

⎦ . (31)

It can be easily verified from Eqs. 29 and 30 that both

C0∗ defined by Eq. 21 and 



∗

11 
∗12


∗

12 
∗22 

are positive definite. It can also be observed from Eq. 30 that the eighteen components (12)=(21), (13)=(31), (24)=(42),

(34)=(43), (26)=(62), (36)=(63), (17)=(71), (47)=

(74) and (67)=(76) of (C ̃ 0∗)−1 do not change due to

the introduction of the insulating microcracks within

the framework of NIA.

For a hexagonal piezoelectric solid with its poling

direction along the x3-axis, it follows from Eq. 30 that

the effective electroelastic properties of the microcracked, piezoelectric solid pertaining to the antiplane

deformation and in-plane electric fields can be explicitly given by

C∗

55

C44

= e∗

15

e15

= 
∗

11


11

= 1

1 + πη[2θ0−sin(2θ0)]

4θ0

,

C∗

44

C44

= e∗

24

e15

= 
∗

22


11

= 1

1 + πη[2θ0+sin(2θ0)]

4θ0

, (32)

C∗

45 = e∗

14 = e∗

25 = 
∗

12 = 0.

For a hexagonal piezoelectric solid with its poling

direction along the x2-axis, it follows from Eq. 30 that

the effective electroelastic properties of the microcracked, piezoelectric solid pertaining to the in-plane

deformation and in-plane electric fields can be explicitly given by

⎡

⎣

C∗

11 C∗12 e∗21

C∗

12 C∗22 e∗22

e∗

21 e∗22 −
∗22

⎤

⎦

−1

=

⎡

⎣

C11 C12 e21

C12 C22 e22

e21 e22 −
22

⎤

⎦

−1

123
On the effective electroelastic properties 33

+πη [2θ0 − sin(2θ0)]

4θ0

⎡

⎣

1

CL 0 0

0 00

0 00

⎤

⎦

+πη [2θ0 + sin(2θ0)]

4θ0

⎡

⎣

00 0

0 1

CT

1

e

0 1

e −1


⎤

⎦ ,




C∗

66 e∗16

e∗

16 −
∗11 −1

=




C66 e16

e16 −
11 −1

+πη [2θ0 − sin(2θ0)]

4θ0

 1

CT

1

e 1

e −1




+πη [2θ0 + sin(2θ0)]

4θ0


 1

CL 0

0 0 

,

C∗

16 = C∗26 = e∗11 = e∗12 = e∗26 = 
∗12 = 0, (33)

where

Re {Y} =

⎡

⎢

⎣

1

CL 0 0

0 1

CT

1

e

0 1

e −1


⎤

⎥

⎦ , (34)

has been defined in Eq. C12 by Suo et al. (1992). The

elements in Eq. 34 have to be determined numerically

for the specific piezoelectric material. We illustrate in

Figs. 3, 4, 5 the variations of the effective electroelastic

moduli as functions of θ0 with η = 0.2 for a microcracked piezoelectric material BaTiO3 with its material

properties given by (Pan 2001)

C11 = C33 = 166 × 109 N/m2,

C13 = 77 × 109 N/m2,

C12 = C23 = 78 × 109 N/m2,

C22 = 162 × 109 N/m2,

C44 = C66 = 43 × 109 N/m2,

C55 = 44.5 × 109 N/m2,

e21 = e23 = −4.4 C/m2,

e22 = 18.6 C/m2,

e34 = e16 = 11.6 C/m2,


11 = 
33 = 11.2 × 10−9 C2/(Nm2),


22 = 12.6 × 10−9 C2/(Nm2).

It can be observed from the three figures that the

degree of the crack alignment characterized by θ0 has

a significant influence on most of the effective electroelastic moduli except C∗

66 due to the fact that the

ratio CT /CL = 1.0736 for BaTiO3 is very close to

unity. We also observe that all the effective electroelastic moduli are monotonically increasing or decreasing

0 10 20 30 40 50 60 70 80 90

2

4

6

8

10

12

14

16

x 1010

θ0 (Deg)

η=0.2

(N/m2

)

11

22

66 C12

* C*

C*

C*

Fig. 3 The variations of the four effective elastic constants

C∗

11,C∗22,C∗12 and C∗66 as functions of θ0 with η = 0.2 for a

piezoelectric material BaTiO3 containing two-dimensional insulating microcracks. The poling direction of the piezoelectric solid

is along the x2-axis

0 10 20 30 40 50 60 70 80 90

-6

-4

-2

0

2

4

6

8

10

12

14

θ0 (Deg)

e*

e*

e*

22

16

21

η=0.2

(C/m2

)

Fig. 4 The variations of the three effective piezoelectric constants e∗

22, e∗21 and e∗16 as functions of θ0 with η = 0.2 for a

piezoelectric material BaTiO3 containing two-dimensional insulating microcracks. The poling direction of the piezoelectric solid

is along the x2-axis

functions of θ0 at a fixed crack density η. It can be

clearly seen from Fig. 3 that the effective elastic anisotropy decreases as θ0 increases from zero for perfectly

123
34 X. Wang et al.

0 10 20 30 40 50 60 70 80 90

0.8

0.85

0.9

0.95

1

1.05

1.1

1.15

1.2

1.25

x 10-8

θ0 (Deg)

(C2/(Nm2))

∈*

∈*

11

22

η=0.2

Fig. 5 The variations of the two effective dielectric constants


∗

11 and 
∗22 as functions of θ0 with η = 0.2 for a piezoelectric

material BaTiO3 containing two-dimensional insulating microcracks. The poling direction of the piezoelectric solid is along

the x2-axis

aligned microcracks to 90◦ for randomly oriented microcracks. When θ0 = 90◦,C∗

11 ≈ C∗22 and C∗66 ≈ 1

2



C∗

11 − C∗12

. We observe from Fig. 5 that 
∗

11 > 
11

when θ0 < 40◦. This observation indicates that the

piezoelectric effect plays an important role in the effective dielectric properties of a microcracked piezoelectric solid. On the other hand our results show that some

of the effective piezoelectric and dielectric moduli at a

certain fixed θ0 can be rather complex functions of η.

This complex phenomenon is demonstrated in Figs. 6

and 7 for e∗

21 and 
∗11 as functions of η for four different values of θ0 = 0,π/8,π/4,π/2. We can see from

Fig. 6 that when θ0 is below a certain value (the detailed

results show that this value is θ0 = π/3.03), e∗

21 can

attain its minimum value at an extremely low crack

density η (for example when θ0 = π/8, e∗

21 attains its

minimum value of −4.5851 C/m2 at η = 0.04); on the

other hand when π/3.03 < θ0 ≤ π/2, e∗

21 is a monotonically increasing function of η. Our detailed results

also show that depending on the value of θ0, there exist

three possible distributions of 
∗

11: (i) a monotonically

increasing function of η when θ0 = 0; (ii) a function

with its maximum value at a certain value of η when

0 < θ0 ≤ π/4.09 (for example as shown in Fig. 7 when

θ0 = π/8, 
∗

11 attains its maximum value of 11.86 ×

10−9 C2/(Nm2) at η = 0.436); (iii) a monotonically

decreasing function of η when θ0 > π/4.09. We add

0 0.2 0.4 0.6 0.8 1

-5

-4.5

-4

-3.5

-3

-2.5

-2

-1.5

-1

-0.5

η

e*21

(C/m2

)

θ =π/2 θ =π/4 θ =π/8 θ =0 0 0 0 0

Fig. 6 The variations of e∗

21 as a function of η at four different values of θ0 = 0,π/8,π/4,π/2 for a piezoelectric material

BaTiO3 containing two-dimensional insulating microcracks. The

poling direction of the piezoelectric solid is along the x2-axis

0 0.2 0.4 0.6 0.8 1

5

6

7

8

9

10

11

12

13

14 x 10-9

η

(C2/(Nm2))

∈*

θ =π/8

θ =π/4

θ =π/2

θ =0 0

0

0

0

11

Fig. 7 The variations of 
∗

11 as a function of η at four different values of θ0 = 0,π/8,π/4,π/2 for a piezoelectric material

BaTiO3 containing two-dimensional insulating microcracks. The

poling direction of the piezoelectric solid is along the x2-axis

that this unexpected phenomenon is caused by the coupling effect between the mechanical field and the electric field (or the piezoelectric effect). For example when

there exists no piezoelectric effect, 
∗

11 ≡ 
11 for perfectly aligned microcracks (θ0 = 0).

3.2 A piezoelectric solid containing permeable

microcracks

The effective electroelastic properties of a piezoelectric solid containing permeable microcracks (i.e., φ+ =

123
On the effective electroelastic properties 35

φ−, D+

n = D−n across the crack surface, see Suo et al.

1992) can be simply derived as

(C ̃ 0∗)

−1 = (C ̃ 0)−1

+πη [2θ0 − sin(2θ0)]

4θ0

JT

3 L−1

11 J3

+πη [2θ0 + sin(2θ0)]

4θ0

JT

4 L−1

11 J4, (35)

where the 3×3 symmetric matrix L11 has been defined

in Eq. 29 and

J3 =

⎡

⎣

1000000

0000100

0001000

⎤

⎦ ,

J4 =

⎡

⎣

0000100

0100000

0010000

⎤

⎦ . (36)

We see from Eq. 35 that the thirty-two components

(12)=(21), (13)=(31), (24)=(42), (34)=(43), (ij)=

(ji), (i=1-7, j=6,7) of (C ̃ 0∗)−1 do not change due to

the introduction of the permeable microcracks within

the framework of NIA. For a hexagonal piezoelectric

solid with its poling direction along the x3-axis, it follows from Eq. 35 that the effective electroelastic properties of the microcracked, piezoelectric solid pertaining to the antiplane deformation and in-plane electric

fields can be explicitly given by

C∗

55

C44

= e∗

15

e15

= 1

1 + πη[2θ0−sin(2θ0)]

4θ0

,


∗

11


11

= 1 +

πk1η [2θ0 − sin(2θ0)]

4θ0 + πη [2θ0 − sin(2θ0)] ≥ 1,

C∗

44

C44

= e∗

24

e15

= 1

1 + πη[2θ0+sin(2θ0)]

4θ0

,


∗

22


11

= 1 +

πk1η [2θ0 + sin(2θ0)]

4θ0 + πη [2θ0 + sin(2θ0)] ≥ 1,

C∗

45 = e∗

14 = e∗

25 = 
∗

12 = 0, (37)

where k1 = e2

15/(
11C44) < 1 is the electromechanical

coupling factor.

We observe from the above expression that only the

expressions for 
∗

11 and 
∗22 are different from the corresponding ones in Eq. 32 for insulating microcracks.

When θ0 is fixed, 
∗

11 and 
∗22 for permeable microcracks are increasing functions of η; whereas those for

insulating microcracks are decreasing functions of η.

Similarly, for a hexagonal piezoelectric solid with

its poling direction along the x2-axis, it follows from

Eq. 35 that the effective electroelastic properties of the

microcracked piezoelectric solid pertaining to the inplane deformation and in-plane electric fields can be

explicitly given by

⎡

⎣

C∗

11 C∗12 e∗21

C∗

12 C∗22 e∗22

e∗

21 e∗22 −
∗22

⎤

⎦

−1

=

⎡

⎣

C11 C12 e21

C12 C22 e22

e21 e22 −
22

⎤

⎦

−1

+ πη

4θ0

⎡

⎢

⎣

2θ0−sin(2θ0)

CL 0 0

0 (1+k2)[2θ0+sin(2θ0)]

k2CT 0

0 00

⎤

⎥

⎦ ,




C∗

66 e∗16

e∗

16 −
∗11 −1

=




C66 e16

e16 −
11 −1

+ πη

4θ0




(1 + k2)[2θ0 − sin(2θ0)]

k2CT

+

2θ0 + sin(2θ0)

CL

 
 1 0

0 0 

,

C∗

16 = C∗26 = e∗11 = e∗12 = e∗26 = 
∗12 = 0, (38)

where k2 = e2/(
CT ).

3.3 A piezoelectric solid containing conducting

microcracks

The effective electroelastic properties of a piezoelectric

solid containing conducting microcracks (McMeeking

1987; Suo 1993) can also be similarly derived. The

energy balance relationship for two-dimensional problems can be expressed in the following form when the

piezoelectric solid containing conducting microcracks

is subjected to the homogeneous stresses σ0

i j and the

homogeneous electric fields E0

i (Suo 1993)

1

2

(σˆ 0

)

T (Cˆ 0∗)−1σˆ 0 = 1

2

(σˆ 0

)

T (Cˆ 0)−1σˆ 0

+

M

2A



θ0

−θ0

φ(θ )

Ck

(ˆt

0)T [uˆ]dCkdθ ,

(39)

where

σˆ 0 =



σ0

11 σ022 σ023 σ013 σ012 E01 E02

T

, (40)

[uˆ] = [[u1] [u2] [u3] [ξ ]]

T , (41)

123
36 X. Wang et al.

Cˆ 0 =

⎡

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎣

C11 C12 C14 C15 C16 h11 h21

C12 C22 C24 C25 C26 h12 h22

C14 C24 C44 C45 C46 h14 h24

C15 C25 C45 C55 C56 h15 h25

C16 C26 C46 C56 C66 h16 h26

h11 h12 h14 h15 h16 β11 β12

h21 h22 h24 h25 h26 β12 β22

⎤

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎦

, (42)

Cˆ 0∗ =

⎡

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎣

C∗

11 C∗12 C∗14 C∗

15 C∗

16 h∗11 h∗21

C∗

12 C∗22 C∗24 C∗

25 C∗

26 h∗12 h∗22

C∗

14 C∗24 C∗44 C∗

45 C∗

46 h∗14 h∗24

C∗

15 C∗25 C∗45 C∗55 C∗56 h∗15 h∗25

C∗

16 C∗26 C∗46 C∗

56 C∗

66 h∗16 h∗26

h∗

11 h∗12 h∗14 h∗

15 h∗

16 β∗11 β∗12

h∗

21 h∗22 h∗24 h∗

25 h∗

26 β∗12 β∗22

⎤

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎦

, (43)

ˆt

0 =

⎡

⎢

⎢

⎢

⎣

σ0

12 cos θ − σ011 sin θ

σ0

22 cos θ − σ012 sin θ

σ0

23 cos θ − σ013 sin θ

E0

1 cos θ + E02 sin θ

⎤

⎥

⎥

⎥

⎦ . (44)

In Eq. 41, ξ is defined such that D1 = ξ,2, D2 =

−ξ,1. In Eqs. 42 and 43, Ci j, hi j and βi j are the elastic, piezoelectric and dielectric constants of the undamaged solid, while C∗

i j, h∗i j and β∗i j are the corresponding

(unknown) constants for the damaged solid. Here we

have adopted the material constant notations in Suo

(1993). When ignoring crack interactions, the jump in

displacement and ξ across the crack surfaces can be

simply given by (Suo 1993)

[uˆ] = 2



a2 − x2Lˆ −1ˆt

0, (|x| ≤ a) (45)

where Lˆ is a 4×4 real and symmetric, positive definite

matrix.

Substituting the above results into Eq. 22, we finally

arrive at a concise expression for the effective electroelastic properties of the microcracked piezoelectric

solid

(Cˆ 0∗)

−1 = (Cˆ 0)−1 + πη[2θ0 − sin(2θ0)]

4θ0

JˆT

1 Lˆ −1Jˆ1

+πη[2θ0 + sin(2θ0)]

4θ0

JˆT

2 Lˆ −1Jˆ2, (46)

where

Jˆ1 =

⎡

⎢

⎢

⎣

100000 0

000010 0

000100 0

000000 −1

⎤

⎥

⎥

⎦,

Jˆ2 =

⎡

⎢

⎢

⎣

0000100

0100000

0010000

0000010

⎤

⎥

⎥

⎦. (47)

It is apparent that the energy density function always

increases through the introduction of conducting microcracks, and that Cˆ 0∗ is positive definite in view of

the fact that both Cˆ 0 and Lˆ are positive definite. Interestingly, this result is very similar to the situation when

the matrix is purely elastic. We stress that the convex

energy density function used here for conducting cracks

is different than the non-convex electric enthalpy used

for insulating cracks.

For a hexagonal piezoelectric solid with its poling

direction along the x3-axis, it follows from Eq. 46 that

the effective electroelastic properties of the microcracked piezoelectric solid pertaining to the antiplane

deformation and in-plane electric fields can be explicitly given by

C∗

55

C44 

1 + πη[2θ0−sin(2θ0)]

4θ0

 = β∗

11

β11 

1 + πη[2θ0−sin(2θ0)]

4θ0



= h∗

15

h15

= C44β11 − h2

15

C44β11 

1 + πη[2θ0−sin(2θ0)]

4θ0

2

− h2

15

,

C∗

44

C44 

1 + πη[2θ0+sin(2θ0)]

4θ0

 = β∗

22

β11 

1 + πη[2θ0+sin(2θ0)]

4θ0



= h∗

24

h15

= C44β11 − h2

15

C44β11 

1 + πη[2θ0+sin(2θ0)]

4θ0

2

− h2

15

,

C∗

45 = h∗

14 = h∗

25 = β∗

12 = 0, (48)

where C44β11 > h2

15. Apparently C∗

44,C∗

55, β∗

11, β∗22 >

0 and C∗

55β∗

11 > h∗2

15,C∗

44β∗22 > h∗224.

For a hexagonal piezoelectric solid with its poling

direction along the x2-axis, it follows from Eq. 46 that

the effective electroelastic properties of the microcracked, piezoelectric solid pertaining to the in-plane

deformation and in-plane electric fields can be explicitly given by

⎡

⎣

C∗

11 C∗12 h∗21

C∗

12 C∗22 h∗22

h∗

21 h∗22 −β∗22

⎤

⎦

−1

=

⎡

⎣

C11 C12 h21

C12 C22 h22

h21 h22 −β22

⎤

⎦

−1

+πη [2θ0 − sin(2θ0)]

4θ0

⎡

⎣

S1 0 −d

000

−d 0 


⎤

⎦

+πη [2θ0 + sin(2θ0)]

4θ0

⎡

⎣

000

0 S3 0

000

⎤

⎦ ,




C∗

66 h∗16

h∗

16 −β∗11 −1

=




C66 h16

h16 −β11 −1

123
On the effective electroelastic properties 37

+πη [2θ0 − sin(2θ0)]

4θ0




S3 0

0 0 

+πη [2θ0 + sin(2θ0)]

4θ0




S1 d

d 




,

C∗

16 = C∗26 = h∗11 = h∗12 = h∗26 = β∗12 = 0, (49)

where

Re {Y} =

⎡

⎣

S1 0 d

0 S3 0

d 0 


⎤

⎦ , (50)

has been defined by Suo (1993). Similar to the

case of insulating cracks, the elements in Eq. 50

have to be computed numerically. Apparently both

⎡

⎣

C∗

11 C∗12 h∗21

C∗

12 C∗22 h∗22

h∗

21 h∗22 β∗22

⎤

⎦ and 


C∗

66 h∗16

h∗

16 β∗11 

are positive definite.

3.4 A special distribution of the insulating

microcracks

Next we discuss a special distribution of the microcracks: all the insulating microcracks with common

half-length a are distributed along the x1-axis in a

homogeneous piezoelectric solid subjected to the

homogeneous stresses σ0

i j and the homogeneous

electric displacements D0

i . As pointed out by Fan and

Sze (2001), we may merge the microcracks along the

x1-axis into a continuously damaged interface such that


 ̃t2 = k ̃
[u ̃], on x2 = 0 (51)

where 
∗ stands for the average over an area of a scale

much greater than the dimension of the microcracks,

 ̃t2 = [σ12 σ22 σ32 D2]

T and [u ̃] being defined by

Eq. 24. In addition 
 ̃t2 =  ̃t

0

2 = 

σ0

12 σ022 σ032 D02

T .

Equation 51 can be termed the generalized, spring-type,

imperfect interface model with the 4 × 4 symmetric

matrix k ̃ being the generalized stiffness matrix to be

determined. When ignoring the crack interactions and

noticing Eqs. 12 and 28 with θ = 0, we can easily

obtain the following


[u ̃] =

πaρ

2

L ̃ −1 ̃t

0

2, (52)

where ρ = Ma/L with M being the number of microcracks within the sample of half-length L taken out

of the x1-axis. Comparison of Eq. 52 with Eq. 51 will

immediately lead to the following closed-form expression of k ̃ as

k ̃ = 2

πaρ

L ̃ , (53)

which remains valid when ρ  1. We see from the

above expression that the structure of k ̃ is exactly the

same as that of L ̃ described in Eq. 29.

3.4.1 Remark

Fan and Sze (2001) derived a closed-form expression in

(2.9) for the spring constant k for two-dimensional noninteracting cracks. In fact, the expression for k in threedimensions for penny-shaped, non-interacting cracks

can also be similarly derived as

k = 3πC

8a

1

ρ2 , (54)

which matches quite well the finite element results in

Table 2 in Fan and Sze (2001) up to ρ = 0.5.

4 An approximation of the GSCM

We point out that the problem in Sect. 2 can also be

discussed within the framework of the generalized selfconsistent method (GSCM), which approximately takes

into account crack interaction (Aboudi and Benveniste

1987; Santare et al. 1995). Here we can adopt a simplified version of the GSCM for a microcracked, elastically anisotropic solid: first calculate the uniform stress

field σ ̃ 0

i j within an intact circularinclusion with undamaged material properties surrounded by the effective

medium with (unknown) material properties subjected

to remote uniform stresses σ0

i j (Ting 1996); second

solve the problem of a crack in an infinite homogeneous

material with undamaged material properties subjected

to remote uniform loading σ ̃ 0

i j obtained in step 1. After

adopting the above, simplified version of the GSCM,

the crack opening displacement can be finally

obtained as

[u] = 2



a2 − x2L−1 ̃t, (|x| ≤ a) (55)

where

 ̃t = L(E1t

0

1 + E2t

0

2) cos θ − L(F1t

0

1 + F2t

0

2) sin θ ,

(56)

with

t

0

1 = [σ011 σ021 σ031]

T , (57)

t

0

2 = [σ012 σ022 σ032]

T ,

[E1 E2] = L−1[03×3 I3×3](N ̃ ∗ + N)

−1

×(N ̃ ∗ + N∗)


 q2 q2N∗T

1

03×3 I3×3



,

123
38 X. Wang et al.

[F1 F2]=−L−1[N3 NT

1 ](N ̃ ∗ + N)

−1

×(N ̃ ∗ + N∗)


 q2 q2N∗T

1

03×3 I3×3



, (58)

and

N =




N1 N2

N3 NT

1



, N ̃ =


 S H

−L ST



, (59)

q2 =

⎡

⎣

S

∗

11 0 S

∗

15

0 00

S

∗

15 0 S

∗

55

⎤

⎦ . (60)

In the above expressions the quantities attached with

the superscript * are those pertaining to the effective

medium. Here we point out that the above derivations

are based on the elliptic inclusion and crack solutions in

Ting (1996) and more detailed structures and identities

of N1, N2, N3, S, H andL, all of which can be explicitly

expressed in terms of the reduced elastic compliances,

can be found in Chapter 6 of Ting (1996). The Barnett–Lothe tensor, S, should not be confused with the

elastic compliance. It can be easily verified from the

above expressions that  ̃t = t

0

2 cos θ − t

0

1 sin θ when

the circular inclusion and the surrounding matrix possess exactly the same material properties. Substituting

the above results into the energy balance expression,

Eq. 3, we finally arrive at the following set of fifteen

coupled, nonlinear equations for the fifteen unknown

effective reduced elastic compliances:

S

∗

11 = S

11 + πηF(1)

11

4θ0

[2θ0 − sin(2θ0)] ,

S

∗

22 = S

22 + πηE(2)

22

4θ0

[2θ0 + sin(2θ0)] ,

S

∗

44 = S

44 + πηE(2)

33

4θ0

[2θ0 + sin(2θ0)] ,

S

∗

55 = S

55 + πηF(1)

33

4θ0

[2θ0 − sin(2θ0)] ,

S

∗

66 = S

66 + πη

4θ0

E(1)

12 + E(2)11 

[2θ0 + sin(2θ0)]

+



F(1)

22 + F(2)21 

[2θ0 − sin(2θ0)]



, (61)

S

∗

12 = S

12 + πη

8θ0

×



E(1)

21 [2θ0 + sin(2θ0)]+F(2)12 [2θ0 − sin(2θ0)]



,

S

∗

14 = S

14 + πη

8θ0

×



E(1)

31 [2θ0 + sin(2θ0)]+F(2)13 [2θ0 − sin(2θ0)]



,

S

∗

15 = S

15 + πη(F(1)

13 + F(1)31 )

8θ0

[2θ0 − sin(2θ0)] ,

S

∗

16 = S

16 + πη

8θ0



E(1)

11 [2θ0 + sin(2θ0)]

+



F(1)

12 + F(1)21 + F(2)11 

[2θ0 − sin(2θ0)]



,

(62)

S

∗

24 = S

24 +

πη E(2)

23 + E(2)32 

8θ0

[2θ0 + sin(2θ0)] ,

S

∗

25 = S

25 + πη

8θ0



E(1)

23 [2θ0 + sin(2θ0)]

+F(2)

32 [2θ0 − sin(2θ0)]



,

S

∗

26 = S

26 + πη

8θ0

E(1)

22 + E(2)12 + E(2)21 

[2θ0 + sin(2θ0)] + F(2)

22 [2θ0 − sin(2θ0)]



,

S

∗

45 = S

45 + πη

8θ0



E(1)

33 [2θ0 + sin(2θ0)]

+F(2)

33 [2θ0 − sin(2θ0)]



,

S

∗

46 = S

46 + πη

8θ0

E(1)

32 + E(2)13 + E(2)31 

[2θ0 + sin(2θ0)] + F(2)

23 [2θ0 − sin(2θ0)]



,

S

∗

56 = S

56 + πη

8θ0



E(1)

13 [2θ0 + sin(2θ0)]

+



F(1)

23 + F(1)32 + F(2)31 

[2θ0 − sin(2θ0)]



,

(63)

where E(1)

i j , E(2)i j , F(1)i j and F(2)i j , all of which are functions of the unknown effective reduced elastic compliances, are the components of four 3 × 3 matrices

E1,E2, F1 and F2 respectively.

Equations 61–63 can be solved through iteration

(Santare et al. 1995; Wang et al. 2009). During specific iteration for the unknowns, it is more convenient

to adopt the following equivalent matrix expression

S∗ = S + πη [2θ0 + sin(2θ0)]

8θ0

×




JT

2 [E1 E2]




J1

J2



+



JT

1 JT2

 


ET

1

ET

2



J2



+πη [2θ0 − sin(2θ0)]

8θ0

×




JT

1 [F1 F2]




J1

J2



+



JT

1 JT2

 


FT

1

FT

2



J1



,

(64)

123
On the effective electroelastic properties 39

where

J1=

⎡

⎣

10000

00001

00010

⎤

⎦, J2=

⎡

⎣

00001

01000

00100

⎤

⎦. (65)

We observe that when using the GSCM, the fifteen

effective, reduced elastic compliances will, in general,

be different from the corresponding compliances for

the undamaged material. Thus the result that S

∗

12 =

S

12, S

∗

14 = S

14, S

∗

25 = S

25, S

∗

45 = S

45 (see Eq. 14) is

only valid for NIA.

5 Effective thermal conductivities of a solid

containing insulated penny-shaped cracks

In the above, the crack orientation distribution function

φ(θ) has been introduced for two-dimensional cracks.

In fact, the crack orientation distribution function can

also be adopted to describe three-dimensional cracks.

As an illustration, we consider the effective thermal

conductivities of a solid containing insulated pennyshaped cracks with an arbitrary degree of alignment.

Within the framework of steady-state heat conduction,

we consider insulated penny-shaped cracks of radius a

with an arbitrary degree of alignment, distributed in a

homogeneous and isotropic medium with thermal conductivity κ. The orientation of the penny-shaped crack

can be described by the azimuth angle θ (|θ| ≤ π/2)

and the zenith angle φ (0 ≤ φ ≤ π) which define

the direction of the unit normal of the crack surface as

shown in Fig. 8 (Yang and Turner 2003). The crack orientation distribution function φ (θ,φ) on a unit sphere

can be given by

φ(θ , φ)

=

⎧

⎨

⎩

1

4θ0 sin φ0

, |θ| ≤ θ0 and π

2 − φ0 ≤ φ ≤ π2 + φ0

0, else

(66)

where θ0, φ0 ≤ π/2. By adopting the energy method

and ignoring the crack interactions, closed-form

expressions for the effective thermal conductivities can

be derived as

κ∗

11 = κ

1 + 2η[2θ0+sin(2θ0)]

3θ0



1 − sin2 φ0

3

,

κ∗

22 = κ

1 + 2η[2θ0−sin(2θ0)]

3θ0



1 − sin2 φ0

3

,

x2

x3

x θ 1

φ

o

Unit sphere

Fig. 8 Geometry of a penny-shaped crack

κ∗

33 = κ

1 + 8η sin2 φ0

9

, (67)

where η = Ma3/V is the crack density parameter and

κ∗

i j = 0, (i = j), which implies that the three principal directions of the effective medium are along the

x1, x2 and x3 axes. The above derivation uses the following expression for the temperature jump [T ] across

the insulated crack surfaces (Barber 1975)

[T ]=−

4q0

√

a2 − r 2

πκ

, (0 ≤ r ≤ a) (68)

where q0 is the uniform heat flow prescribed on the

crack face.

If all the microcracks are randomly oriented, i.e.,

φ0 = θ0 = π/2, then we have

κ∗

11 = κ∗22 = κ∗33 = κ

1 + 8

9η , (69)

which means that the effective property is still isotropic.

If the unit normals of all the microcracks lie in the

x1 − x2 plane, i.e., φ0 = 0, then we have

κ∗

11 = κ

1 + 2η[2θ0+sin(2θ0)]

3θ0

,

κ∗

22 = κ

1 + 2η[2θ0−sin(2θ0)]

3θ0

,

κ∗

33 = κ, (70)

which indicates that the effective thermal conductivity κ∗

33 is not influenced by the existence of the microcracks. Furthermore, if the unit normals, which lie

in the x1 − x2 plane, have a random distribution, i.e.,

θ0 = π/2 (or the so called uniaxially aligned cracks,

see Yang and Turner 2003), then it follows from Eq. 70

that

κ∗

11 = κ∗22 = κ

1 + 4

3η

, κ∗

33 = κ, (71)

123
40 X. Wang et al.

which indicates that the effective property is transversely isotropic with the x3-axis taken as the uniaxial

symmetry axis. Alternatively, if the unit normals lying

in the x1 − x2 plane are along the x1-axis, i.e., θ0 = 0

(i.e, the so-called perfectly aligned cracks, see Yang

and Turner 2005), then it follows from Eq. 70 that

κ∗

11 = κ

1 + 8

3η

, κ∗

22 = κ∗33 = κ, (72)

which indicates that the effective property is also transversely isotropic, but now with the x1-axis as the

uniaxial symmetry axis.

6 Conclusions

We first derived in Eq. 14 the closed-form expressions

for the fifteen effective reduced elastic compliances of

a microcracked, elastically anisotropic solid by using

NIA. Then we further obtained in Eqs. 30, 35 and 46

the concise expressions for the effective electroelastic

properties of a piezoelectric solid containing insulating, permeable or conducting microcracks by extending the NIA for purely elastic materials. We also derived

a set of coupled nonlinear Eq. 64 for the unknown effective reduced elastic compliances by using the GSCM.

Finally the crack orientation distribution function was

introduced to investigate the effective thermal conductivities of a solid containing insulated penny-shaped

cracks with an arbitrary degree of alignment, and the

effective thermal conductivities were derived in Eq. 67.

Acknowledgement The authors wish to thank the two referees

for their comments and suggestions on revising the initial manuscript. This research was supported by the United States Army

Research Laboratory through the Composite Materials Technology cooperative agreement with the Center for Composite Materials at the University of Delaware.

References

Aboudi J, Benveniste Y (1987) The effective moduli of cracked

bodies in plane deformations. Eng Fract Mech 26(2):171–

184. doi:10.1016/0013-7944(87)90195-0

Barber JR (1975) Steady-state thermal stresses in an elastic solid

containing an insulated penny-shaped crack. J Strain Anal

10(1):19–24. doi:10.1243/03093247V101019

Benveniste Y (1985) The effective mechanical behavior of

composite materials with imperfect contact between the

constituents. Mech Mater 4(2):197–208. doi:10.1016/

0167-6636(85)90016-X

Bristow JR (1960) Microcracks, and the static and dynamic elastic constants of annealed heavily cold-worked metals. Br J

Appl Phys 11:81–85. doi:10.1088/0508-3443/11/2/309

Eshelby JD (1956) The continuum theory of lattice defects. In:

Seitz F, Turnbull D (eds) Progress in solid state physics.

Academic Press, New York p 79

Fan H, Sze KY (2001) A micro-mechanics model for imperfect

interface in dielectric materials. Mech Mater 33:363–370.

doi:10.1016/S0167-6636(01)00053-9

Feltman RS, Santare MH (1999) Anisotropic effective moduli

of cracked short fiber reinforced composites. ASME J Appl

Mech 66:709–713

Gottesman T, Hashin Z, Brull MA (1980) Effective elastic

moduli of cracked fiber composites. In: Bunsell et al (ed)

Advance in composites materials. Pergamon Press, Oxford

pp 749–758

Hashin Z (1988) The differential scheme and its application to

cracked materials. J Mech Phys Solids 36:719–734. doi:10.

1016/0022-5096(88)90005-1

Kachanov M (1992) Effective elastic properties of cracked solids: critical review of some basic concepts. Appl Mech Rev

45(8):305–336. doi:10.1115/1.3119761

Kachanov M (2007) On the effective elastic properties of cracked

solids—editor’s comments. Int J Fract 146:295–299. doi:10.

1007/s10704-007-9170-6

Mauge C, Kachanov M (1994) Effective elastic properties of

an anisotropic material with arbitrarily oriented interacting

cracks. J Mech Phys Solids 42(4):561–584. doi:10.1016/

0022-5096(94)90052-3

McMeeking RM (1987) On mechanical stresses at cracks in

dielectrics with application to dielectric-breakdown. J Appl

Phys 62(8):3116–3122. doi:10.1063/1.339361

Pak YE (1990) Crack extension force in a piezoelectric material.

ASME J Appl Mech 57:647–653

Pan E (2001) Exact solution for simply supported and multilayered magneto-electro-elastic plates. ASME J Appl Mech

68:608–618

Santare MH, Crocombe AD, Anlas G (1995) Anisotropic effective moduli of materials with microcracks. Eng Fract Mech

52(5):833–842. doi:10.1016/0013-7944(95)00059-5

Suo Z (1990) Singularities, interfaces and cracks in dissimilar

anisotropic materials. Proc R Soc Lond A Math Phys Sci

427:331–358. doi:10.1098/rspa.1990.0016

Suo Z (1993) Models for breakdown-resistant dielectric and ferroelectric ceramics. J Mech Phys Solids 41:1155–1176.

doi:10.1016/0022-5096(93)90088-W

Suo Z, Kuo CM, Barnett DM, Willis JR (1992) Fracture mechanics for piezoelectric ceramics. J Mech Phys Solids 40:739–

765. doi:10.1016/0022-5096(92)90002-J

Ting TCT (1996) Anisotropic elasticity-theory and applications.

Oxford University Press, New York

Wang X, Santare MH, Gazonas GA (2009) Anisotropic effective

moduli of microcracked materials under antiplane loading.

Eng Fract Mech (in press)

Yang LY, Turner JA (2003) Elastic wave propagation and scattering in solids with uniaxially aligned cracks. J Acoust Soc

Am 114(2):591–600. doi:10.1121/1.1592158

Yang LY, Turner JA (2005) Wave attenuations in solids with perfectly aligned cracks. Acoust Res Lett Online 6(2):99–105.

doi:10.1121/1.1861092

123
NO. OF 

COPIES ORGANIZATION 

 1 DEFENSE TECHNICAL 

 (PDF INFORMATION CTR 

 only) DTIC OCA 

 8725 JOHN J KINGMAN RD 

 STE 0944 

 FORT BELVOIR VA 22060-6218 

 1 DIRECTOR 

 US ARMY RESEARCH LAB 

 IMNE ALC HRR 

 2800 POWDER MILL RD 

 ADELPHI MD 20783-1197 

 1 DIRECTOR 

 US ARMY RESEARCH LAB 

 RDRL CIM L 

 2800 POWDER MILL RD 

 ADELPHI MD 20783-1197 

 1 DIRECTOR 

 US ARMY RESEARCH LAB 

 RDRL CIM P 

 2800 POWDER MILL RD 

 ADELPHI MD 20783-1197 

ABERDEEN PROVING GROUND

 1 DIR USARL 

 RDRL CIM G (BLDG 4600) 
NO. OF NO. OF 

COPIES ORGANIZATION COPIES ORGANIZATION

 1 DPTY ASSIST SCT FOR R&T 

 (CD SARD TT 

 only) ASA (ACT) 

 J PARMENTOLA 

 THE PENTAGON RM 3E479 

 WASHINGTON DC 20310-1714 

 1 PRIN DPTY FOR TCHNLGY HQ 

 US ARMY MATCOM 

 AMCDCGR 

 R PRICE 

 9301 CHAPEK RD 

 FT BELVOIR VA 22060-5527 

 3 AIR FORCE ARMAMENT LAB 

 AFATL DLJW 

 W COOK 

 D BELK 

 J FOSTER 

 EGLIN AFB FL 32542 

 2 NSF 

 S MCKNIGHT 

 G PAULINO 

 4201 WILSON BLVD STE 545 

 ARLINGTON VA 22230-0002 

 1 DARPA 

 W COBLENZ 

 3701 N FAIRFAX DR 

 ARLINGTON VA 22203-1714 

 1 DIRECTOR 

 US ARMY ARDEC 

 AMSRD AAR AEE W 

 E BAKER 

 BLDG 3022 

 PICATINNY ARSENAL NJ 

 07806-5000 

 2 US ARMY TARDEC 

 AMSTRA TR R MS 263 

 K BISHNOI 

 D TEMPLETON MS 263 

 WARREN MI 48397-5000 

 1 COMMANDER 

 US ARMY RSRCH OFC 

 AMSRD ARL RO EN 

 B LAMATTINA 

 PO BOX 12211 

 RESEARCH TRIANGLE PARK NC 

 27709-2211 

 1 COMMANDER 

 US ARMY RSRCH OFC 

 AMSRD ARL RO MM 

 J LAVERY 

 PO BOX 12211 

 RESEARCH TRIANGLE PARK NC 

 27709-2211 

 1 COMMANDER 

 US ARMY RSRCH OFC 

 AMSRD ARL RO EM 

 D STEPP 

 PO BOX 12211 

 RESEARCH TRIANGLE PARK NC 

 27709-2211 

 5 NAVAL RESEARCH LAB 

 E R FRANCHI CODE 7100 

 M H ORR CODE 7120 

 J A BUCARO CODE 7130 

 J S PERKINS CODE 7140 

 S A CHIN BING CODE 7180 

 4555 OVERLOOK AVE SW 

 WASHINGTON DC 20375 

 1 DTRA 

 M GILTRUD 

 8725 JOHN J KINGMAN RD 

 FORT BELVOIR VA 22060 

 1 ERDC 

 US ARMY CORPS OF ENGRS 

 USACEGSL 

 P PAPADOS 

 7701 TELEGRAPH RD 

 ALEXANDRIA VA 22315 

 1 AFOSR NL 

 F FAHROO 

 875 N RANDOLPH ST 

 STE 325 RM 3112 

 ARLINGTON VA 22203 

 1 CLEMSON UNIV 

 DEPT MECH ENGR 

 M GRUJICIC 

 241 ENGRG INNOVATION BLDG 

 CLEMSON SC 29634-0921 

 1 UNIV OF CALIFORNIA 

 CTR OF EXCELLENCE FOR 

 ADV MATLS 

 S NEMAT NASSER 

 SAN DIEGO CA 92093-0416 
NO. OF NO. OF 

COPIES ORGANIZATION COPIES ORGANIZATION

 5 DIRECTOR 

 LANL 

 P MAUDLIN 

 R GRAY 

 W R THISSELL 

 A ZUREK 

 F ADDESSIO 

 PO BOX 1663 

 LOS ALAMOS NM 87545 

 7 DIRECTOR 

 SANDIA NATL LABS 

 J BISHOP MS 0346 

 E S HERTEL JR MS 0382 

 W REINHART MS 1181 

 T VOGLER MS 1181 

 L CHHABILDAS MS 1811 

 M FURNISH MS 1168 

 M KIPP MS 0378 

 PO BOX 5800 

 ALBUQUERQUE NM 87185-0307 

 1 DIRECTOR 

 LLNL 

 M J MURPHY 

 PO BOX 808 

 LIVERMORE CA 94550 

 3 CALTECH 

 M ORTIZ MS 105 50 

 G RAVICHANDRAN 

 T J AHRENS MS 252 21 

 1201 E CALIFORNIA BLVD 

 PASADENA CA 91125 

 5 SOUTHWEST RSRCH INST 

 C ANDERSON 

 K DANNEMANN 

 T HOLMQUIST 

 G JOHNSON 

 J WALKER 

 PO DRAWER 28510 

 SAN ANTONIO TX 78284 

 1 TEXAS A&M UNIV 

 DEPT OF MATH 

 J WALTON 

 COLLEGE STATION TX 77843 

 1 WORCESTER POLYTECHNIC INST 

 MATHEMATICAL SCI 

 K LURIE 

 WORCESTER MA 01609 

 2 SRI INTERNATIONAL 

 D CURRAN 

 D SHOCKEY 

 333 RAVENSWOOD AVE 

 MENLO PARK CA 94025 

 1 VIRGINIA POLYTECHNIC INST 

 COLLEGE OF ENGRG 

 R BATRA 

 BLACKSBURG VA 24061-0219 

 8 UNIV OF NEBRASKA 

 DEPT OF ENGRG MECH 

 D ALLEN 

 F BOBARU 

 Y DZENIS 

 G GOGOS 

 M NEGAHBAN 

 R FENG 

 J TURNER 

 Z ZHANG 

 LINCOLN NE 68588 

 1 JOHNS HOPKINS UNIV 

 DEPT OF MECH ENGRG 

 K T RAMESH 

 LATROBE 122 

 BALTIMORE MD 21218 

 4 UNIV OF UTAH 

 DEPT OF MATH 

 A CHERKAEV 

 E CHERKAEV 

 E S FOLIAS 

 R BRANNON 

 SALT LAKE CITY UT 84112 

 1 PENN STATE UNIV 

 DEPT OF ENGRG SCI & MECH 

 F COSTANZO 

 UNIVERSITY PARK PA 168023 

 3 UNIV OF DELAWARE 

 DEPT OF MECH ENGRG 

 T BUCHANAN 

 T W CHOU 

 M SANTARE 

 126 SPENCER LAB 

 NEWARK DE 19716 

 1 UNIV OF DELAWARE 

 CTR FOR COMPST MATRLS 

 J GILLESPIE 

 NEWARK DE 19716
NO. OF NO. OF 

COPIES ORGANIZATION COPIES ORGANIZATION

 1 COMPUTATIONAL MECH 

 CONSULTANTS 

 J A ZUKAS 

 PO BOX 11314 

 BALTIMORE MD 21239-0314 

 1 LOUISIANA STATE UNIV 

 R LIPTON 

 304 LOCKETT HALL 

 BATON ROUGE LA 70803-4918 

 1 INST OF ADVANCED TECH 

 UNIV OF TX AUSTIN 

 S BLESS 

 3925 W BRAKER LN STE 400 

 AUSTIN TX 78759-5316 

 1 APPLIED RSCH ASSOC 

 D E GRADY 

 4300 SAN MATEO BLVD NE 

 STE A220 

 ALBUQUERQUE NM 87110 

 1 INTERNATIONAL RSRCH 

 ASSOC INC 

 D L ORPHAL 

 4450 BLACK AVE 

 PLEASANTON CA 94566 

 3 ORNL 

 ENVIRONMENTAL SCI DIV 

 W DOLL 

 T GAMEY 

 L BEARD 

 PO BOX 2008 

 OAK RIDGE TN 37831 

 1 UNIV OF ILLINOIS 

 DEPT OF MECHL SCI & ENGRG 

 A F VAKAKIS 

 1206 W GREEN ST MC 244 

 URBANA CHAMPAIGN IL 61801 

 1 UNIV OF ILLINOIS 

 ARSPC ENGRG 

 J LAMBROS 

 104 S WRIGHT ST MC 236 

 URBANA CHAMPAIGN IL 61801 

 2 WASHINGTON ST UNIV 

 INST OF SHOCK PHYSICS 

 Y M GUPTA 

 J ASAY 

 PULLMAN WA 99164-2814

 1 NORTHWESTERN UNIV 

 DEPT OF CIVIL & ENVIRON ENGRG 

 Z BAZANT 

 2145 SHERIDAN RD A135 

 EVANSTON IL 60208-3109 

 1 UNIV OF DAYTON 

 RSRCH INST 

 N S BRAR 

 300 COLLEGE PARK 

 MS SPC 1911 

 DAYTON OH 45469 

 2 TEXAS A&M UNIV 

 DEPT OF GEOPHYSICS MS 3115 

 F CHESTER 

 T GANGI 

 COLLEGE STATION TX 778431 

1 UNIV OF SAN DIEGO 

 DEPT OF MATH & CMPTR SCI 

 A VELO 

 5998 ALCALA PARK 

 SAN DIEGO CA 92110 

 1 NATIONAL INST OF 

 STANDARDS & TECHLGY 

 BLDG & FIRE RSRCH LAB 

 J MAIN 

 100 BUREAU DR MS 8611 

 GAITHERSBURG MD 20899-8611 

 1 MIT 

 DEPT ARNTCS ASTRNTCS 

 R RADOVITZKY 

 77 MASSACHUSETTS AVE 

 CAMBRIDGE MA 02139 

 1 UNIV OF DELAWARE 

 DEPT ELECTRICAL & CMPTR ENGRG 

 D WEILE 

 NEWARK DE 19716 

 3 DIR USARL 

 RDRL D 

 V WEISS 

 RDRL SE 

 J PELLEGRINO 

 RDRL SES P 

 A EDELSTEIN 

 2800 POWDER MILL RD 

 ADELPHI MD 20783-1197 
NO. OF NO. OF 

COPIES ORGANIZATION COPIES ORGANIZATION

 2 MATERIALS SCI CORP 

 A CAIAZZO 

 R LAVERTY 

 181 GIBRALTAR RD 

 HORSHAM PA 19044 

 1 UNIV OF MISSISSIPPI 

 DEPT OF MECH ENGRG 

 A M RAJENDRAN 

 201 B CARRIER HALL 

 UNIVERSITY MS 38677 

ABERDEEN PROVING GROUND

 67 DIR USARL 

 RDRL WM 

 B FORCH 

 S KARNA 

 J MCCAULEY 

 P PLOSTINS 

 J SMITH 

 RDRL WML 

 J NEWILL 

 M ZOLTOSKI 

 RDRL WML B 

 I BATYREV 

 B RICE 

 N WEINGARTEN 

 RDRL WML D 

 P CONROY 

 M NUSCA 

 RDRL WML G 

 J ANDZELM 

 M BERMAN 

 W DRYSDALE 

 A RAWLETT 

 RDRL WML H 

 D SCHEFFLER 

 S SCHRAML 

 B SCHUSTER 

 RDRL WMM 

 R DOWDING 

 J ZABINSKI 

 RDRL WMM A 

 M MAHER 

 J TZENG 

 E WETZEL 

 RDRL WMM B 

 T BOGETTI 

 B CHEESEMAN 

 C FOUNTZOULAS 

 D HOPKINS 

 B POWERS 

 C RANDOW 

 M VANLANDINGHAM 

 R WILDMAN 

 C F YEN 

 RDRL WMM D 

 E CHIN 

 K CHO 

 R HOWELL 

 RDRL WMM E 

 M COLE 

 T JESSEN 

 J LASALVIA 

 J SANDS 

 RDRL WMM F 

 L KECSKES 

 S MATHAUDHU 

 RDRL WMP 

 P BAKER 

 S SCHOENFELD 

 RDRL WMP B 

 R BECKER 

 S BILYK 

 D CASEM 

 J CLAYTON 

 C HOPPEL 

 M GREENFIELD 

 R KRAFT 

 B LEAVY 

 B LOVE 

 M SCHEIDLER 

 T WEERASOORIYA 

 RDRL WMP C 

 T BJERKE 

 S SEGLETES 

 RDRL WMP D 

 R DONEY 

 J RUNYEON 

 B SCOTT 

 RDRL WMP E 

 M BURKINS 

 RDRL WMP F 

 M CHOWDHURY 

 A FRYDMAN 

 N GNIAZDOWSKI 

 R GUPTA 

 RDRL WMP G 

 D KOOKER 

 G R PEHRSON