0
Research Papers

Frequency-Dependent Tensile and Compressive Effective Moduli of Elastic Solids With Randomly Distributed Two-Dimensional Microcracks OPEN ACCESS

[+] Author and Article Information
Youxuan Zhao

School of Civil Engineering,
Southwest Jiaotong University,
Chengdu 610031, China
Department of Civil
and Environmental Engineering,
Department of Mechanical Engineering,
Northwestern University,
Evanston, IL 60208

Yanjun Qiu

School of Civil Engineering,
Southwest Jiaotong University,
Chengdu 610031, China

Laurence J. Jacobs

College of Engineering,
Georgia Institute of Technology,
Atlanta, GA 30332-0360

Jianmin Qu

Department of Civil
and Environmental Engineering,
Department of Mechanical Engineering,
Northwestern University,
Evanston, IL 60208
e-mail: j-qu@northwestern.edu

When the crack is closed, the crack opening displacement describes the relative sliding of the crack faces (Mode II).

1Corresponding author.

Contributed by the Applied Mechanics Division of ASME for publication in the JOURNAL OF APPLIED MECHANICS. Manuscript received March 20, 2015; final manuscript received May 3, 2015; published online June 9, 2015. Editor: Yonggang Huang.

J. Appl. Mech 82(8), 081006 (Aug 01, 2015) (13 pages) Paper No: JAM-15-1150; doi: 10.1115/1.4030538 History: Received March 20, 2015; Revised May 03, 2015; Online June 09, 2015

This paper develops micromechanics models to estimate the tensile and compressive elastic moduli of elastic solids containing randomly distributed two-dimensional microcracks. The crack faces are open under tension and closed under compression. When the crack faces are closed, they may slide against one another following the Coulomb's law of dry friction. The micromechanics models provide analytical expressions of the tensile and compressive moduli for both static and dynamic cases. It is found that the tensile and compressive moduli are different. Further, under dynamic loading, the compressive and tensile moduli are both frequency dependent. As a by-product, the micromechanics models also predict wave attenuation in the dynamic case. Numerical simulations using the finite element method (FEM) are conducted to validate the micromechanics models.

FIGURES IN THIS ARTICLE
<>

Microcracks are often observed in quasi-brittle materials such as ceramics, rocks, and concrete. As such, techniques to estimate the effective modulus of such cracked solids have been studied extensively in the literature. A seminal work in this area is the paper by Budiansky and O'Connell [1], in which they developed a self-consistent scheme to estimate the static effective modulus of an elastic solid containing randomly distributed flat elliptical cracks. Explicit expressions of the effective modulus are given in terms of the elastic modulus of the uncracked solid and the crack density. One of the assumptions made in their work is that the cracks remain open all the time, i.e., the crack faces do not come into contact, and thus remain traction-free. This is the case when the solid is subjected to tensile loading. In other words, the effective modulus obtained in Ref. [1] can be viewed as the effective tensile modulus.

When the solid is under compressive loading, the cracks tend to close. Thus, crack faces may come into contact and slide against each other. Due to this tension and compression asymmetry of the cracks, the effective compressive modulus of a cracked solid may differ from its effective tensile modulus [2-4]. Furthermore, sliding of the crack faces may lead to a compressive modulus that is load-dependent, and may even be affected by the sequence of the load applications [3,4]. These early works derived explicit expressions of the effective compressive modulus in terms of the modulus of the uncracked solid, the crack density, and the coefficient of friction between the crack faces.

When a time-harmonic longitudinal wave propagates through a solid, it subjects the solid to alternating tension and compression. In the presence of distributed microcracks, the solid will respond to the tensile and compressive cycles differently, due to its tension and compression elastic asymmetry. Consequently, the waveform will be distorted and higher order harmonic waves are generated [5]. In particular, the magnitude of the second harmonic wave is proportional to the difference between the effective tensile and compressive moduli. Therefore, by measuring the amplitude of the second harmonics, the effective tensile and compressive moduli of the cracked solid can be obtained. Since the effective tensile and compressive moduli are both related to the crack density in the solid, measuring the second harmonic may provide a method to nondestructively characterize material damage due to microcracking. For quantitative nondestructive characterization, however, the relationship between the crack density and the tensile and compressive moduli needs to be developed. This is the objective of the current study.

As discussed earlier, most of the existing results on tensile and compressive moduli of cracked solids are for the static case [1-4]. When a time-harmonic longitudinal wave propagates through a cracked solid, the solid is subjected to alternating tension and compression. Under such dynamic loading, the effective modulus during the tensile cycle is likely frequency-dependent, and thus different from the effective tensile modulus under static tension. Such frequency-dependent effective modulus may lead to wave dispersion. The same can be said about the effective compressive modulus between the dynamic and static cases. In other words, the dynamic effective tensile and compressive moduli may be different from their static counterparts. Further, they may be frequency dependent and complex. Such complex valued elastic moduli give rise to attenuation of the waves [6].

Although wave propagation in cracked solids has been a subject of intensive studies, e.g., Refs. [7-19], and various micromechanics schemes have been developed to estimate the dynamic effective modulus of solids containing distributed inclusions or cracks, e.g., Refs. [20-23], very few studies have considered crack face contact during the compressive cycles of the wave motion. Most of these studies assumed that the solid is under a pretension so that the cyclic wave motion is not large enough to close the cracks.

When the crack faces are allowed to come into contact, the tension and compression asymmetry introduces nonlinearity into the wave field. The first study on such contact-induced acoustic nonlinearity is probably the paper by Richardson [24] who considered the contact interface between two semi-infinite half spaces. Since then, the topic has been studied in considerable details [25-36]. However, most of these studies focused on either the scattering of elastic waves by a single or an array of cracks, or on the propagation of elastic waves in a cracked medium. Such studies, although in principle, may lead to the dynamic effective modulus of cracked solids, the actual practice is not straightforward [35].

In the present work, attention is focused on obtaining the dynamic effective tensile and compressive moduli of elastic solids containing randomly distributed two-dimensional microcracks. Our results show that (1) such dynamic moduli are complex and frequency-dependent. The imaginary part of the complex moduli gives rise to wave attenuation, and (2) the compressive and tensile moduli are different, and the compressive modulus depends on the friction between the crack faces in contact. In a separate publication [5], we show that such tension and compression asymmetry gives rise to the acoustic nonlinearity, which forms the basis for nonlinear ultrasonic nondestructive evaluation techniques, e.g., Refs. [37-42].

As a longitudinal wave propagates through a cracked medium, the cracks are subjected to cyclic tension and compression. Therefore, it is somewhat ambiguous to distinguish effective tensile and compressive moduli during a wave motion. To overcome this difficulty, we propose a two-step approach in this paper. In step 1, we treat the tensile and compressive cycles in a wave motion separately as static loadings. In doing so, we obtain the “static” or instantaneous effective tensile and compressive moduli as functions of the crack opening displacement (COD).2 In this step, we allow the cracks to open under tension and close under compression. Crack faces are also allowed to slide against one another with friction. In step 2, we calculate the COD by treating the crack as a scatterer insonified by a harmonic plane wave. To be precise, crack face contact should be considered in calculating the COD. However, doing so will make the scattering problem nonlinear and the solution can only be obtained numerically [43]. To reduce the complexity, the COD is calculated in this step by assuming that the crack remains open and its faces do not come into contact. Such an assumption introduces negligible error at the end, because crack face contact has already been accounted for in step 1, when the effective compressive modulus was derived as a function of the COD. To validate this assertion, results based on this assumption will be compared with numerical simulations using the FEM, which provides the numerically exact solutions.

The paper is arranged as follows. The problem to be solved is described in Sec. 2, where the governing equations for the field quantities are presented and the eigenstrain induced by the cracks is written as a function of the applied load and the COD. Micromechanics models for the effective tensile and compressive moduli are derived based on a dilute-concentration approach in Sec. 3 and on a self-consistent approach in Sec. 4. The dilute-concentration approach assumes that the cracks are far apart from each other, so that they do not interact, while the interactions among the cracks are partially taken into account in the self-consistent approach [44]. The effective tensile and compressive moduli obtained from Secs. 3 and 4 are both functions of the COD, which is obtained in the Appendix. Both micromechanics models yield complex and frequency-dependent effective moduli. They are used in Sec. 5 to obtain the dispersion relationship and scattering-induced attenuation. To validate the micromechanics models, numerical results from the FEM simulations are presented in Sec. 6. Comparison between the micromechanics model predictions and the FEM simulation results shows good agreement. Finally, a summary and some remarks are provided in Sec. 7.

Consider an isotropic elastic solid of area A containing N randomly distributed and randomly oriented two-dimensional microcracks (Griffith cracks) of length 2a. By microcracks, we mean that the crack length is much smaller than the other characteristic lengths of interest, such as the size of the representative area A, or the wavelength in case of dynamic loading. Let E be the Young's modulus and ν be the Poisson's ratio of the uncracked solid, and E¯α and ν¯α be, respectively, the corresponding effective Young's modulus and Poisson's ratio of the cracked solid under tension (α=t) and under compression (α=c). We note that the elastic moduli of cracked solids might be anisotropic under general loading condition [2-4]. However, for the uniaxial loading considered here, E¯α and ν¯α are sufficient to describe the uniaxial deformation as will be seen later.

Let n be a unit vector normal to the faces of a typical crack, see Fig. 1. Without loss of generality, we can assume that

where -π/2φπ/2. Further, we assume that the solid is subjected to surface traction p=σ¯·m on its boundary A, where σ¯ is a constant stress tensor and m is the outward unit normal of A.

By a simple superposition argument, the original problem can be decomposed into two problems: (1) a solid A without any crack is subjected to surface traction p=σ¯·m on its boundary A; and (2) a solid A with N randomly oriented cracks of length 2a, where the surfaces of the cracks are subjected to applied normal and shear stresses σ0 and τ0, respectively. Further, we assume that the crack faces, when sliding against each other, follow the Coulomb dry friction law, i.e., the frictional force against sliding of the crack faces is proportional to the normal pressure.Display Formula

(2)f=μσ0H(-σ0)

where μ is the coefficient of friction and H() is the Heaviside step function. Therefore, by accounting for the frictional forces, the actual traction applied to the crack faces is thus given byDisplay Formula

(3)q=σ0n+(τ0+f)H(τ0+f)s=σ0n+[τ0+μσ0H(-σ0)]H[τ0+μσ0H(-σ0)]s

where f is the frictional force given in Eq. (2), andDisplay Formula

(4)σ0=n·σ¯·n,   τ0=|σ·n-σ0n|,   s=σ·n-σ0n|σ·n-σ0n|

Clearly, σ0 is the normal stress to open the crack, τ0 is the shear stress in the direction of s to cause the crack faces to slid over each other, and μσ0H(-σ0) is the frictional force that resists the sliding of crack faces when they are pressed against each other. Note that s is in the plane of the crack.

The COD for a Griffith crack under unit normal and shear tractions can be written asDisplay Formula

(5)b(φ)u+-u-=4E'[gn(ω,φ,x)σ0H(σ0)n   +gs(ω,φ,x)(τ0+f)H(τ0+f)s]a2-x2

where x is the distance from the crack center in the plane of the crack, and E'=E for plane stress, and E'=E/(1-ν2) for plane strain. The dimensionless functions gn(ω,φ,x) and gs(ω,φ,x) are derived in the Appendix. It is noticed here that gn(ω,φ,x) and gs(ω,φ,x) are even functions of φ, and are complex-valued. In the static limit of ω0, one has gn(0,φ,x)=gs(0,φ,x)=1. Thus, Eq. (5) reduces to the well-known result for a Griffith crack under static loading, e.g., Refs. [45] and [46].

Once the COD is known, the corresponding strain field can be written as Ref. [44]Display Formula

(6)ɛ*(φ,x)=12[b(n)n+nb(n)]δl(2a-x)

where the symbol means cross product between two vectors and δl(2a-x) is the line Dirac delta function defined byDisplay Formula

(7)δl(2a-x)=2aδ(z-x)dl(z)

and the integration is along the length of the crack. It can be easily shown that for any continuous function defined in an area A that contains the line segment 2a, the line Dirac delta function satisfies the following:Display Formula

(8)Af(x)δl(2a-x)dA(x)=2af(x)dl(x)

Substituting Eq. (5) into Eq. (6) and making use of Eq. (8) yieldsDisplay Formula

(9)ɛ*(φ,x)=2a2-x2E'[2σ0gn(ω,φ,x)H(σ0)nn   +gs(ω,φ,x)(τ0+f)H(τ0+f)(sn+ns)]δl(2a-x)

By dilute concentration, we mean that the crack density is low enough that the interactions among the cracks can be neglected. Under such dilute-concentration assumption, the strain field induced by all N cracks is simply the sum of all the strain fields induced by each individual crack as if it were all by itself, i.e.,Display Formula

(10)ɛ*(x)-π/2π/2ɛ*(n,x)ψ(n)sinφdφ=Nπ-π/2π/2ɛ*(n,x)dφ

where ψ(n) is the crack orientation distribution. The second equality in Eq. (10) follows from the assumption that the cracks are randomly oriented, i.e., ψ(n)=N/π.

The average strain in the area A is thus,Display Formula

(11)ɛ¯=ɛ¯0+1AAɛ*(x)dA(x)=ɛ¯0+ɛ¯*

where ɛ¯0 is the strain field induced by σ¯ in the absence of the cracks and ɛ¯* is the average strain induced by the traction q applied to the crack faces, i.e.,Display Formula

(12)ɛ¯*1AAɛ*(x)dA(x)=N2πA-π/2π/2{l[b(φ)n+nb(φ)]dl(x)}dφ

By carrying out the integral over l, Eq. (12) can be further simplified toDisplay Formula

(13)ɛ¯*=c2πa2-π/2π/2[v(φ)n+nv(φ)]dφ

where a is the average half-length of the cracks, c=Na2/A is the crack density, and the crack opening area v(n) is given byDisplay Formula

(14)v(φ)=lb(φ)dl=2πa2E'[gn(ω,φ)σ0H(σ0)n   +gs(ω,φ)(τ0+f)H(τ0+f)s]

In Eq. (14)Display Formula

(15)gβ(ω,φ)=2πa2-aagβ(ω,φ,x)a2-x2dx,   β=n,s

where gβ(ω,φ,x) for β=n,s are derived in the Appendix, and the integral in Eq. (15) has been carried out to yield explicit expressions for gβ(ω,φ), see Eq. (A101). Combining Eqs. (13) and (14) producesDisplay Formula

(16)ɛ¯*=cE'-π/2π/2[2gn(ω,φ)σ0H(σ0)nn   +gs(ω,φ)(τ0+f)H(τ0+f)(sn+ns)]dφ

We now consider uniaxial loading in the e1-direction, i.e., σ¯=pe1e1. This leads toDisplay Formula

(17)ɛ¯=pE¯'α(e1e1-ν¯'e2e2),   ɛ¯0=pE'[e1e1-ν'(e2e2)]
Display Formula
(18)σ0=pcos2φ,   τ0=|pcosφsinφ|,   f=pμcos2φH(-p)
Display Formula
(19)v(φ)=2πa2pE'[gn(ω,φ)cos2φH(p)n   +gs(ω,φ)[|p|pcosφsinφ+μcos2φH(-p)]H(τ0+f)s]
Display Formula
(20)s=sgn(p)sgn(φ)(sinφe1-cosφe2)

In the case of uniaxial tension, i.e., p>0, one has f=0 and τ0+f>0. Thus, Eq. (16) is simplified toDisplay Formula

(21)ɛ¯*=cpE'-π/2π/2[2gn(ω,φ)cos2φnn   +gs(ω,φ)|cosφsinφ|(sn+ns)]dφ

Making use of Eqs. (1) and (20) in Eq. (21), and noting that the gn(ω,φ) and gs(ω,φ) given in Eq. (A101) are even functions of φ, we haveDisplay Formula

(22)ɛ¯*=πcpE'(h1te1e1+h2te2e2)

whereDisplay Formula

(23)h1t=2π-π/2π/2[gn(ω,φ)cos4φ+gs(ω,φ)sin2φcos2φ]dφ=1-Aη2lnη+[κ2-1+κ4(1-2logκ)16κ2(κ2-1)   +A(log4+logκ-γ)]η2+iπA2η2
Display Formula
(24)h2t=2π-π/2π/2[gn(ω,φ)-gs(ω,φ)]cos2φsin2φdφ=-Bη2lnη+[1+2γ-2γκ2-2log232κ2   +B(log2-logκ)]η2+iπB2η2

In the above, γ is the Euler–Mascheroni constant, η=kTa=ωa/cT is a dimensionless wave number andDisplay Formula

(25)A=5-6κ2+5κ416κ2(κ2-1),   B=κ2-116κ2

We note that A and B are dimensionless parameters that depend on the Poisson's ratio only. They are independent of the frequency.

Making use of Eqs. (22) and (17) in Eq. (11), we obtain the equations for the effective tensile Young's modulus E¯t and effective tensile Poisson's ratio ν¯tDisplay Formula

(26)1E¯'t=1E'+πch1tE',   ν¯'tE¯'t=ν'E'-πch2tE'

Solving these equations leads toDisplay Formula

(27)E¯'tE'=11+πch1t,   ν¯'t=ν'-πch2t1+πch1t

when h1t and h2t are complex, the effective modulus becomes complex, which gives rise to wave attenuation.

In the static limit (η0), it can be easily shown from Eqs. (23) and (24) that h1t=1, h2t=0. Thus,Display Formula

(28)E¯'tE'=11+πc,   ν¯t=ν'1+πc

These are the same static properties derived in Ref. [47] under the dilute-concentration assumption.

In the case of uniaxial compression, one has p<0. Thus,Display Formula

(29)σ0=pcos2φ,   τ0=|pcosφsinφ|,   f=pμcos2φ,   s=-sgn(φ)(sinφe1-cosφe2)
Display Formula
(30)H(τ0+f)=H[|cosφsinφ|-μcos2φ)]=H(|tanφ|-μ)

The eigenstrain then follows from substituting Eq. (30) and the first three equations in Eq. (29) into Eq. (16)Display Formula

(31)ɛ¯*=pcE'-π/2π/2gs(ω,φ)(-cosφ|sinφ|+μcos2φ)H(|tanφ|-μ)   ×(sn+ns)dφ                                          (31)

Making use of Eq. (1) and the last equation in Eq. (29) into Eq. (31) yieldsDisplay Formula

(32)ɛ¯*=πcpE'h1c(e1e1-e2e2)

where by making use of Eq. (A101)Display Formula

(33)h1c=2π-π/2π/2gs(ω,φ)(|tanφ|-μ)cos2φ|sinφ|cosφH(|tanφ|-μ)dφ=F4-F(κ4+1)32κ2(κ2-1)η2lnη+F[κ4+2(1+κ4)(log4-γ)+2logκ]64κ2(κ2-1)η2   -132πκ2[tan-11μ-μ(3μ2+1)3(μ2+1)2]η2+iπF(κ4+1)64κ2(κ2-1)η2(33)

In the above, the parameter F is related to the coefficient of friction throughDisplay Formula

(34)F=2π(tan-11μ-μμ2+1)

By substituting Eqs. (17) and (32) into Eq. (11), we obtain the equations for the compressive effective Young's modulus and Poisson's ratioDisplay Formula

(35)1E¯'c=1E'+πch1cE',   ν¯'cE¯'c=ν'E'+cπh1cE'

These lead toDisplay Formula

(36)E¯'cE'=11+πch1c,   ν¯'c=ν'+πch1c1+πch1c

In the static limit (η0)Display Formula

(37)E¯'cE'=44+πFc,   ν¯'c=4ν'+πFc4+πFc

The compressive effective modulus and Poisson's ratio given in Eqs. (36) and (37) appear to be new to the literature.

Note that F=1 when μ=0, and F0 when μ. It is seen that, in the static limit, the product Fc stays together all the time as an effective crack density. In other words, F determines the amount of cracks that may have sliding faces. When the crack faces are smooth (i.e., μ=0), Eq. (34) yields F=1, meaning that all the cracks may participate in sliding. When the crack surfaces are extremely rough (i.e., μ), Eq. (34) yields F0. Thus, none of the cracks can slide. Under compression, it is the crack face sliding that produces the additional compressive deformation. So, in the limit of μ, one can see from Eq. (33) that h1c0. Thus, Eqs. (36) and (37) immediately lead to E¯'c=E' and ν¯'c=ν', as if the cracks were absent.

Before closing this section, it is worthwhile to point out that the crack density c=Na2/A is the only parameter introduced here to characterize the cracks, where N is the number of cracks within area A, and a is the half-crack length. For a fixed N, increasing c means increasing crack size. On the other hand, for a fixed a, increasing c means increasing the number of cracks per unit area. In other words, the crack density c represents a combined effect of crack size and number of cracks per unit area. It turns out that the value of c is usually much smaller than one, i.e., c1. For example, if every square of side 4a contains on average one crack of length 2a, then c0.06. Even if in each square of side 4a, there are on average five cracks of length 2a, one still has c0.3. Five cracks of length 2a in a square of side 4a is a rather high crack density.

At higher crack densities, crack interactions must be considered. One way to do so is by using the self-consistent scheme [44]. The self-consistent scheme assumes that each crack is embedded in an effective medium with E¯α and ν¯α, which are yet to be determined. Consequently, the COD given in Eq. (5) should be written asDisplay Formula

(38)b(n)=4E¯'[g¯n(ν¯α,ω,φ,x)σ0H(σ0)n   +g¯s(ν¯α,ω,φ,x)(τ0+f)H(τ0+f)s]a2-x2

where an overbar indicates that the quantity is evaluated based on the effective properties E¯α and ν¯α. For clarity, g¯n(ν¯α,ω,φ,x) and g¯s(ν¯α,ω,φ,x) are explicitly written as functions of ν¯α. Dimensional analysis indicates that they are independent of E¯α.

Therefore, under uniaxial tension, Eq. (26) becomesDisplay Formula

(39)1E¯'t=1E'+πch¯1t(ν¯'t)E¯'t,   ν¯'tE¯'t=ν'E'-πch¯2t(ν¯'t)E¯'t

where h¯1t(ν¯'t) and h¯2t(ν¯'t) are obtained from Eqs. (23) and (24), respectively, by replacing the Poisson's ratio ν' with the effective Poisson's ratio ν¯'t. Combining the two equations in Eq. (39) provides a pair of nonlinear algebraic equations for E¯'t and ν¯'tDisplay Formula

(40)E¯'tE'=1-πch¯1t(ν¯'t),   E¯'tE'=ν¯'t+πch¯2t(ν¯'t)ν'

Further, combining the two equations in Eq. (40) yields a nonlinear algebraic equation for ν¯'tDisplay Formula

(41)c=ν'-ν¯'tπh¯2t(ν¯'t)+πν'h¯1t(ν¯'t)

Once ν¯'t is solved from Eq. (41), it can be substituted into either of the two equations in Eq. (40) to find E¯'t. The ν¯'t and E¯'t so solved are the self-consistent estimates of the effective tensile properties of the cracked solid.

In the static limit, Eq. (41) and the second of Eq. (40) reduce toDisplay Formula

(42)ν¯'tν'=1-πc,   E¯'tE'=1-πc

These are the same static properties given in Ref. [47].

Analogously, under uniaxial compression, Eq. (35) becomesDisplay Formula

(43)1E¯'c=1E'+πch¯1c(ν¯'c)E¯'c,   ν¯'cE¯'c=ν'E'+πch¯1c(ν¯'c)E¯'c

where h¯1c(ν¯'c) is calculated from Eq. (33) by replacing the Poisson's ratio ν with the effective Poisson's ratio ν¯'c. Equation (43) can be recast intoDisplay Formula

(44)E¯'cE'=1-πch¯1c(ν¯'c),   E¯'cE'=ν¯'c-πch¯1c(ν¯'c)ν'

Combination of the above renders the following nonlinear algebraic equation for ν¯'c:Display Formula

(45)c=ν'-ν¯'cπν'h¯1c(ν¯'c)-πh¯1c(ν¯'c)

Once ν¯'c is solved (45), it can be substituted into either Eq. (43) or (44) to obtain E¯'c. This completes the self-consistent estimate of the effective properties of a cracked solid under uniaxial compression.

In the static limit, h¯1c(ν¯'c)=1. Thus,Display Formula

(46)ν¯'cν'=1+π(1-ν')Fc4ν',   E¯'cE'=1-πFc4

Since the elastic moduli are complex numbers, wave propagation is typically dispersive and attenuated. Thus, we introduce a complex wavenumber byDisplay Formula

(47)kL=k(ω)+iα(ω)=ωρλ+2μ

whereDisplay Formula

(48)λ=12[E¯'tν¯'t(1+ν¯'t)(1-2ν¯'t)+E¯'cν¯'c(1+ν¯'c)(1-2ν¯'c)]μ=12[E¯'t2(1+ν¯'t)+E¯'c2(1+ν¯'c)]

are the average (between tension and compression) Lamé constants of the cracked solid. Clearly, k(ω) gives the dispersion relationship and α(ω) is the coefficient of attenuation. The phase and group velocities are then given byDisplay Formula

(49)cp=ωk(ω),   cg=dωdk=cp+kdcpdk

Figure 2 shows the phase velocity (normalized by the phase velocity of the uncracked solid) versus crack density. It is seen that the dilute-concentration and self-consistent estimates yield very similar results. Both predict decreasing phase velocity with increasing crack density. Figure 3 shows the normalized phase velocity versus frequency. It seems that the self-consistent estimate predicts slightly higher phase velocity than does the dilute-concentration estimate. Both micromechanics models predict that plane waves in cracked medium are dispersive in that their phase velocity depends on the frequency of the wave, albeit only slightly.

Next, consider the scattering-induced attenuation. Notice from Eqs. (23), (24), and (33) that the imaginary parts of h1t, h2t, and h1c are all on the order of η2. It can then be shown through straightforward asymptotic analyses that the coefficient of attenuation α(ω) introduced in Eq. (47) is proportional to the third power of frequency, i.e., α(ω)ω3. This is consistent with the asymptotic solution obtained in Ref. [48]. In fact, as demonstrated by Nagy [49], scattering-induced attenuation by finite size scatterers in the low frequency regime (Rayleigh scattering) is scaled with α(ω)ωn+1, where n is the dimension of the scatterers. Therefore, for two-dimensional cracks, i.e., n = 2, Rayleigh scattering leads to α(ω)ω3. In a separate publication [50], we will show that the scattering-induced attenuation by a distribution of penny-shaped cracks is scaled with α(ω)ω4.

Smyshlyaev and Willis [43] showed that when the crack faces are allowed to be in contact with frictional sliding, the scattering cross section is scaled with frequency square. One may then conclude that the attenuation should be α(ω)ω2. They attribute such behavior to shock waves radiated at the moments of opening and closure of the cracks faces. It appears that the effects of such shock waves are not captured by the current micromechanics models.

The solid and dashed lines in Fig. 4 show the model-predicted attenuation versus frequency. It is seen that the dilute-concentration and self-consistent estimates yield very similar results. The symbols are from the FEM simulations to be discussed in Sec. 6.

In Secs. 3 and 4, we have developed micromechanics models for the frequency-dependent tensile and compressive moduli of elastic solids with randomly distributed and oriented microcracks. These micromechanics models are now validated by comparing their predictions with the numerical simulation results in this section. The numerical simulations are performed using the FEM. The commercial FEM software abaqus is used for this purpose.

A two-dimensional FEM model is constructed using the four-node plane strain (CPE4R) elements. The size of the model is L×L and it contains N microcracks. For numerical expediency, all the cracks have the same length of 2a. They are randomly distributed spatially, and their orientations are also random. Crack faces are modeled as contact surfaces in that crack faces can separate from, but cannot interpenetrate into each other. On contact, the dry friction law is used to model the relative sliding of the crack faces as discussed earlier. Figure 5 shows a typical FEM model with random cracks.

The material parameters used in the FEM analysis are taken from a typical polycrystalline, aluminum. They are ρ=2700kg/m3, E=7×1010Pa, and v=0.33, which yield a longitudinal phase velocity of cL=6198m/s. A linear elastic constitutive law is used in the FEM simulations.

In simulating the dynamic loading, a displacement excitation u(x,t)=Asin(ωt) is prescribed on the left edge of the L×L FEM model. As a result, a plane longitudinal wave is generated to propagate in the x-direction (from the left to the right of the sample). The wavelength of the fundamental frequency is Λ=2πc/ω, and that of the second harmonic is Λ/2. In this simulation problem, there are four characteristic length parameters, namely, the sample size L, the wavelength Λ of interest, the half crack length a, and the FEM mesh size h. These length parameters need to be chosen carefully.

First, L/a needs to be large enough so that the FEM model is representative of a solid with a random distribution of microcracks of random orientations. To this end, we constructed several FEM models with various values of L/a. We found that L/a=120 yields a good compromise between computational accuracy and efficiency.

Second, Λ/a needs to be large enough so the low frequency assumption is valid and the cracked solid can be treated as an effective medium for wave propagation. In this work, ω=5πMHz is used, which corresponds to a wavelength of Λ=2.48mm. The half crack length used is a=0.05mm. These give us Λ/a50, which is well within the low frequency regime.

Finally, in order to capture wave motion in the FEM simulations, there should be no less than 20 elements per wavelength of the highest frequencies [51,52], i.e., Λ/(2h)>20. So, in our FEM model, the element size is kept at h<Λ/40, considering that the wavelength of the second harmonic is Λ/2.

An additional consideration is the randomness of the crack distribution. To capture the stochastic nature of the random distribution, multiple FEM models are constructed with different realizations of the random distribution. Results reported here are the averages of 100 FEM models. In Figs. 4 and 6–8, the FEM results are indicated by symbols with error bars. The size of the error bars indicates the extreme values among the different realizations.

In the FEM model described above, meshing becomes extremely difficult for higher (c>0.04) crack densities. Therefore, for c>0.04, a different approach is used, where, instead of constructing an FEM model with many cracks, we include only one crack in each L×L square. Periodic boundary conditions are prescribed on all four sides of the square. Thus, this square can be viewed as a “unit cell” in an infinite solid containing multiple cracks. To account for the randomness of the crack distribution, we constructed many such squares, each containing a single crack of a random orientation. Results from each of such single squares are then averaged. The number of squares used is large enough so that their average converges (i.e., adding more squares will not change the average anymore).

The above FEM models are also used to simulate the static loading. The static effective modulus is obtained by applying a static uniaxial stress to the FEM model and then computing the corresponding strain in the direction of the load. The compressive and tensile moduli are both obtained in this way. Figure 6 shows the static tensile modulus as a function of the crack density for the plane strain case. Solid and dashed lines are model predictions and symbols are from the FEM models. The error bars indicate the extreme values from all the realizations. It is seen that both the dilute-concentration and self-consistent estimates give good predictions of the static modulus at lower crack densities, although the dilute estimate seems to be closer to the FEM simulations. At higher crack densities, the micromechanics model predictions start to deviate from the FEM simulation results.

The static compressive modulus is shown in Fig. 7 for μ=0.5. Again, it is seen that both the dilute-concentration and self-consistent estimates yield good agreement with the FEM simulations. Although not shown, we carried out FEM simulations for various values of μ. As expected, increasing friction leads to higher compressive modulus. In the limit of μ, the compressive modulus approaches that of the uncracked solid.

In the dynamic simulations, wave attenuation is also calculated from the FEM models and compared with the micromechanics model predictions from Eq. (47). Figure 8 shows the attenuation versus crack density for a fixed frequency (η=0.25). It is seen that, at very low crack densities, the micromechanics model predictions show a linear increase in the attenuation with increasing crack density, and are very close to the FEM results. However, the FEM results seem to show a nonlinear increase in the attenuation with increasing crack density.

This paper develops micromechanics models to estimate the effective tensile and compressive elastic moduli of elastic solids containing randomly distributed two-dimensional microcracks. Two approaches, dilute-concentration and self-consistent, are employed to obtain explicit expressions of the effective tensile and compressive moduli. The solutions are valid for both static and dynamic cases. As a by-product, the micromechanics models also predict wave dispersion and attenuation in the dynamic case. These micromechanics models are validated by comparing their predictions with the numerical simulations using the FEM. Based on our micromechanics model predictions, the following observations can be made.

First, effective tensile and compressive moduli of a cracked solid are different. The effective tensile modulus depends on the crack density in a nonlinear fashion. Predictions from the dilute-concentration and self-consistent estimates agree well at the lower crack densities, and start to deviate from each other at higher crack densities. Somewhat counter-intuitively, the dilute-concentration model seems to yield better agreement with the FEM simulations. This could be due to the fact that the self-consistent approaches generally perform poorly when the inhomogeneities, such as cracks, voids or rigid inclusions, have a strong contrast with the matrix material [44].

Second, the effective compressive modulus is dependent on the friction of the crack faces in contact. The effective compressive modulus is the lowest when the crack faces are smooth, and it increases with increasing friction coefficient. In the limit of extreme roughness crack faces, the effective compressive modulus reaches that of the uncracked solid. The relationship between the effective compressive modulus and the coefficient of friction is rather nonlinear. However, a dimensionless parameter F=2π(tan11/μμ/μ2+1) is introduced, which varies between 0 (when μ) and 1 (when μ0). This parameter gives a measure of the percentage of the cracks that can be activated to slide under uniaxial compression.

Third, under dynamic loading, both the effective compressive and tensile moduli are complex and frequency-dependent. Consequently, a plane wave propagating in a cracked medium will be attenuated and dispersive. In the Rayleigh scattering regime considered here, our micromechanics models predict that the scattering-induced attenuation is scaled with the third power of frequency. Further, our model predictions indicate that the dispersion is rather weak. For example, the phase velocity decreases only by a few percent, when the dimensionless wavenumber increases from 0 to 0.5. However, both phase velocity and attenuation are strongly affected by crack density. The phase velocity is reduced by more than 20%, when the dimensionless crack density is increased to 0.2.

Finally, the explicit nature of the solutions derived in this paper enables the development of a relationship between the microcrack density and the acoustic nonlinearity parameter, which provides a useful tool for nondestructive evaluation of material damage [5].

The work was supported in part by the U.S. National Science Foundation through CMMI-1363221 and in part by the U.S. Department of Energy's Nuclear Energy University Program through Standard Research Contracts Nos. 00126931 and 00127346.

Appendix—Low Frequency Scattering by a Griffith Crack

Consider an isotropic and linearly elastic solid of infinite extent. Let the elastic stiffness tensor of the solid be givenDisplay Formula

(A1)Cijkl=λδijδkl+μ(δikδjl+δilδjk)

where λ and μ (not to be confused with the coefficient of friction introduced in the main body of the paper) are the Lamé constants. They are related to the Young's modulus E and Poisson's ratio ν throughDisplay Formula

(A2)λ=E'ν'(1+ν')(1-2ν'),   μ=E'2(1+ν')

where E'=E, ν'=ν for plane stress, and E'=E/(1-ν2), ν'=ν/(1-ν) for plane strain.

We assume that a two-dimensional Griffith crack of length 2a is contained in this infinite solid. For convenience, we attach a Cartesian coordinate system xi to the crack such that the crack occupies the line segment

L{x;x2=0,|x1|1}

Let an incident longitudinal time-harmonic plane wave propagate in the direction of pjDisplay Formula

(A3)vj=apjexp(ikLaxmpm)=apjexp(iηxmpm/κ)

where kL=ω/cL and kT=ω/cT are, respectively, the longitudinal and transverse wavenumbers, κ=cL/cT=2(1-ν)/(1-2ν) is a dimensionless parameter that depends only on the Poisson's ratio, and η=kTa is a dimensionless wavenumber. The time factor e-iωt has been omitted. The pertinent corresponding stresses areDisplay Formula

(A4)qj(x1,x2)=Cj2klvk,l=iη(τ0δj1+σ0δj2)exp(iηxmpm/κ)/κ

where σ0=C22klpkpl, τ0=C12klpkpl. Further, using the angle φ introduced in Fig. 1, we haveDisplay Formula

(A5)p1=sinφ,   p2=cosφ

Interaction between the incident wave and the crack produces a scattered wave field so that the total wave field ujtotal is the sum of the incident and the scattered fields, i.e., ujtotal=vj+uj. If the crack faces remain traction-free, the boundary value problem for governing the scattered field uj is given byDisplay Formula

(A6)Cijkluk,lj(x1,x2)+(λ+2μ)η2ui(x1,x2)=0
Display Formula
(A7)Δuj(x1)uj(x1,0+)-uj(x1,0-)=0for|x1|1
Display Formula
(A8)tj(x1)Cj2kluk,l(x1,0)=-qj(x1,0)=-iη(τ0δj1+σ0δj2)exp(iηx1p1/κ)/κfor|x1|<1

This boundary value problem can be solved asymptotically for small η (low frequency approximation) by a perturbation technique developed in Refs. [12] and [53]. Here, we briefly outline the procedure.

To this end, we express the scattered displacement field in an integral formDisplay Formula

(A9)u=-A(ξ,x2)E(ξ,x2)B-1(ξ,x2)h(ξ)exp(-iξx1η)dξ

It can be shown by direct substitution that the above satisfies equations of motion, ifDisplay Formula

(A10)A(ξ)=[-iξisgn(x2)γTisgn(x2)γLiξ],B(ξ)=μ[2sgn(x2)ξγL2ξ2-12ξ2-1-2sgn(x2)ξγT]
Display Formula
(A11)E(ξ,x2)=[exp[isgn(x2)γLηx2]00exp[isgn(x2)γTηx2]]
Display Formula
(A12)γL={1/κ2-ξ2for|ξ|<1/κiξ2-1/κ2for|ξ|1/κγT={1-ξ2for|ξ|<1iξ2-1for|ξ|1

The corresponding traction vector is thenDisplay Formula

(A13)t=η-B(ξ,x2)E(ξ,x2)B-1(ξ,x2)h(ξ)exp(-iξx1η)dξ

The boundary conditions (A8) and (A7) at x2=0 lead toDisplay Formula

(A14)Δu(x1)=2-L(ξ)h(ξ)exp(-iξx1η)dξ=0for|x1|>1
Display Formula
(A15)t(x1,0)=η-h(ξ)exp(-iξx1η)dξ=q(x1)for|x1|<1

whereDisplay Formula

(A16)L(ξ)=-iμ[4ξ2γLγT+(2ξ2-1)2][γT00γL]

For convenience, let f(ξ)=L(ξ)h(ξ), then Eqs. (A14) and (A15) can be recast into a pair of dual integral equations.Display Formula

(A17)2-f(ξ)exp(-iξx1η)dξ=0for|x1|>1
Display Formula
(A18)η-L-1(ξ)f(ξ)exp(-iξx1η)dξ=q(x1,0)for|x1|<1

In component formDisplay Formula

(A19)2-fi(ξ)exp(-iξx1η)dξ=0for|x1|>1
Display Formula
(A20)ημ-Kik(ξ)fk(ξ)exp(-iξx1η)dξ=qi(x1,0)for|x1|<1

where K21(ξ)=K12(ξ)=0 andDisplay Formula

(A21)K11(ξ)=iγT[4ξ2γLγT+(2ξ2-1)2],K22(ξ)=iγL[4ξ2γLγT+(2ξ2-1)2]

It is easy to verify thatDisplay Formula

(A22)limξ[K11(ξ)/ξ]=limξ[K22(ξ)/ξ]=-2(1-1κ2)-m

We can now divide the problem into mode IDisplay Formula

(A23)2-f1I(ξ)exp(-iξx1η)dξ=0for|x1|>1
Display Formula
(A24)η-K11(ξ)f1I(ξ)exp(-iξx1η)dξ=0for|x1|<1
Display Formula
(A25)2-f2I(ξ)exp(-iξx1η)dξ=0for|x1|>1
Display Formula
(A26)ημ-K22(ξ)f2I(ξ)exp(-iξx1η)dξ=-iησ0exp(iηx1p1/κ)/κfor|x1|<1

and mode IIDisplay Formula

(A27)2-f1II(ξ)exp(-iξx1η)dξ=0for|x1|>1
Display Formula
(A28)η-K11(ξ)f1II(ξ)exp(-iξx1η)dξ=-iητ0exp(iηx1p1/κ)/κfor|x1|<1
Display Formula
(A29)2-f2II(ξ)exp(-iξx1η)dξ=0for|x1|>1
Display Formula
(A30)η-K22(ξ)f2II(ξ)exp(-iξx1η)dξ=0for|x1|<1