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.

Introduction

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.

Problem Statement

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

 
n=cosφe1+sinφe2
(1)

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. 
f=μσ0H(-σ0)
(2)
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 by 
q=σ0n+(τ0+f)H(τ0+f)s=σ0n+[τ0+μσ0H(-σ0)]H[τ0+μσ0H(-σ0)]s
(3)
where f is the frictional force given in Eq. (2), and 
σ0=n·σ¯·n,   τ0=|σ·n-σ0n|,   s=σ·n-σ0n|σ·n-σ0n|
(4)

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 as 
b(φ)u+-u-=4E'[gn(ω,φ,x)σ0H(σ0)n   +gs(ω,φ,x)(τ0+f)H(τ0+f)s]a2-x2
(5)

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] 
ɛ*(φ,x)=12[b(n)n+nb(n)]δl(2a-x)
(6)
where the symbol means cross product between two vectors and δl(2a-x) is the line Dirac delta function defined by 
δl(2a-x)=2aδ(z-x)dl(z)
(7)
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: 
Af(x)δl(2a-x)dA(x)=2af(x)dl(x)
(8)
Substituting Eq. (5) into Eq. (6) and making use of Eq. (8) yields 
ɛ*(φ,x)=2a2-x2E'[2σ0gn(ω,φ,x)H(σ0)nn   +gs(ω,φ,x)(τ0+f)H(τ0+f)(sn+ns)]δl(2a-x)
(9)

Dilute-Concentration Estimates

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., 
ɛ*(x)-π/2π/2ɛ*(n,x)ψ(n)sinφdφ=Nπ-π/2π/2ɛ*(n,x)dφ
(10)

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, 
ɛ¯=ɛ¯0+1AAɛ*(x)dA(x)=ɛ¯0+ɛ¯*
(11)
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., 
ɛ¯*1AAɛ*(x)dA(x)=N2πA-π/2π/2{l[b(φ)n+nb(φ)]dl(x)}dφ
(12)
By carrying out the integral over l, Eq. (12) can be further simplified to 
ɛ¯*=c2πa2-π/2π/2[v(φ)n+nv(φ)]dφ
(13)
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 by 
v(φ)=lb(φ)dl=2πa2E'[gn(ω,φ)σ0H(σ0)n   +gs(ω,φ)(τ0+f)H(τ0+f)s]
(14)
In Eq. (14) 
gβ(ω,φ)=2πa2-aagβ(ω,φ,x)a2-x2dx,   β=n,s
(15)
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) produces 
ɛ¯*=cE'-π/2π/2[2gn(ω,φ)σ0H(σ0)nn   +gs(ω,φ)(τ0+f)H(τ0+f)(sn+ns)]dφ
(16)
We now consider uniaxial loading in the e1-direction, i.e., σ¯=pe1e1. This leads to 
ɛ¯=pE¯'α(e1e1-ν¯'e2e2),   ɛ¯0=pE'[e1e1-ν'(e2e2)]
(17)
 
σ0=pcos2φ,   τ0=|pcosφsinφ|,   f=pμcos2φH(-p)
(18)
 
v(φ)=2πa2pE'[gn(ω,φ)cos2φH(p)n   +gs(ω,φ)[|p|pcosφsinφ+μcos2φH(-p)]H(τ0+f)s]
(19)
 
s=sgn(p)sgn(φ)(sinφe1-cosφe2)
(20)
In the case of uniaxial tension, i.e., p>0, one has f=0 and τ0+f>0. Thus, Eq. (16) is simplified to 
ɛ¯*=cpE'-π/2π/2[2gn(ω,φ)cos2φnn   +gs(ω,φ)|cosφsinφ|(sn+ns)]dφ
(21)
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 have 
ɛ¯*=πcpE'(h1te1e1+h2te2e2)
(22)
where 
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
(23)
 
h2t=2π-π/2π/2[gn(ω,φ)-gs(ω,φ)]cos2φsin2φdφ=-Bη2lnη+[1+2γ-2γκ2-2log232κ2   +B(log2-logκ)]η2+iπB2η2
(24)
In the above, γ is the Euler–Mascheroni constant, η=kTa=ωa/cT is a dimensionless wave number and 
A=5-6κ2+5κ416κ2(κ2-1),   B=κ2-116κ2
(25)

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 ν¯t 
1E¯'t=1E'+πch1tE',   ν¯'tE¯'t=ν'E'-πch2tE'
(26)
Solving these equations leads to 
E¯'tE'=11+πch1t,   ν¯'t=ν'-πch2t1+πch1t
(27)

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, 
E¯'tE'=11+πc,   ν¯t=ν'1+πc
(28)

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, 
σ0=pcos2φ,   τ0=|pcosφsinφ|,   f=pμcos2φ,   s=-sgn(φ)(sinφe1-cosφe2)
(29)
 
H(τ0+f)=H[|cosφsinφ|-μcos2φ)]=H(|tanφ|-μ)
(30)
The eigenstrain then follows from substituting Eq. (30) and the first three equations in Eq. (29) into Eq. (16) 
ɛ¯*=pcE'-π/2π/2gs(ω,φ)(-cosφ|sinφ|+μcos2φ)H(|tanφ|-μ)   ×(sn+ns)dφ                                          (31)
(31)
Making use of Eq. (1) and the last equation in Eq. (29) into Eq. (31) yields 
ɛ¯*=πcpE'h1c(e1e1-e2e2)
(32)
where by making use of Eq. (A101) 
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)
(33)
In the above, the parameter F is related to the coefficient of friction through 
F=2π(tan-11μ-μμ2+1)
(34)
By substituting Eqs. (17) and (32) into Eq. (11), we obtain the equations for the compressive effective Young's modulus and Poisson's ratio 
1E¯'c=1E'+πch1cE',   ν¯'cE¯'c=ν'E'+cπh1cE'
(35)
These lead to 
E¯'cE'=11+πch1c,   ν¯'c=ν'+πch1c1+πch1c
(36)
In the static limit (η0) 
E¯'cE'=44+πFc,   ν¯'c=4ν'+πFc4+πFc
(37)

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.

Self-Consistent Estimates

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 as 
b(n)=4E¯'[g¯n(ν¯α,ω,φ,x)σ0H(σ0)n   +g¯s(ν¯α,ω,φ,x)(τ0+f)H(τ0+f)s]a2-x2
(38)

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) becomes 
1E¯'t=1E'+πch¯1t(ν¯'t)E¯'t,   ν¯'tE¯'t=ν'E'-πch¯2t(ν¯'t)E¯'t
(39)
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 ν¯'t 
E¯'tE'=1-πch¯1t(ν¯'t),   E¯'tE'=ν¯'t+πch¯2t(ν¯'t)ν'
(40)
Further, combining the two equations in Eq. (40) yields a nonlinear algebraic equation for ν¯'t 
c=ν'-ν¯'tπh¯2t(ν¯'t)+πν'h¯1t(ν¯'t)
(41)

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 to 
ν¯'tν'=1-πc,   E¯'tE'=1-πc
(42)

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

Analogously, under uniaxial compression, Eq. (35) becomes 
1E¯'c=1E'+πch¯1c(ν¯'c)E¯'c,   ν¯'cE¯'c=ν'E'+πch¯1c(ν¯'c)E¯'c
(43)
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 into 
E¯'cE'=1-πch¯1c(ν¯'c),   E¯'cE'=ν¯'c-πch¯1c(ν¯'c)ν'
(44)
Combination of the above renders the following nonlinear algebraic equation for ν¯'c: 
c=ν'-ν¯'cπν'h¯1c(ν¯'c)-πh¯1c(ν¯'c)
(45)

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, 
ν¯'cν'=1+π(1-ν')Fc4ν',   E¯'cE'=1-πFc4
(46)

Wave Dispersion and Attenuation

Since the elastic moduli are complex numbers, wave propagation is typically dispersive and attenuated. Thus, we introduce a complex wavenumber by 
kL=k(ω)+iα(ω)=ωρλ+2μ
(47)
where 
λ=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)]
(48)
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 by 
cp=ωk(ω),   cg=dωdk=cp+kdcpdk
(49)

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.

Numerical Validation of the Micromechanics Models

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.

Summary and Conclusions

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].

Acknowledgment

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 given 
Cijkl=λδijδkl+μ(δikδjl+δilδjk)
(A1)
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 ν through 
λ=E'ν'(1+ν')(1-2ν'),   μ=E'2(1+ν')
(A2)

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 pj 
vj=apjexp(ikLaxmpm)=apjexp(iηxmpm/κ)
(A3)
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 are 
qj(x1,x2)=Cj2klvk,l=iη(τ0δj1+σ0δj2)exp(iηxmpm/κ)/κ
(A4)
where σ0=C22klpkpl, τ0=C12klpkpl. Further, using the angle φ introduced in Fig. 1, we have 
p1=sinφ,   p2=cosφ
(A5)
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 by 
Cijkluk,lj(x1,x2)+(λ+2μ)η2ui(x1,x2)=0
(A6)
 
Δuj(x1)uj(x1,0+)-uj(x1,0-)=0for|x1|1
(A7)
 
tj(x1)Cj2kluk,l(x1,0)=-qj(x1,0)=-iη(τ0δj1+σ0δj2)exp(iηx1p1/κ)/κfor|x1|<1
(A8)

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 form 
u=-A(ξ,x2)E(ξ,x2)B-1(ξ,x2)h(ξ)exp(-iξx1η)dξ
(A9)
It can be shown by direct substitution that the above satisfies equations of motion, if 
A(ξ)=[-iξisgn(x2)γTisgn(x2)γLiξ],B(ξ)=μ[2sgn(x2)ξγL2ξ2-12ξ2-1-2sgn(x2)ξγT]
(A10)
 
E(ξ,x2)=[exp[isgn(x2)γLηx2]00exp[isgn(x2)γTηx2]]
(A11)
 
γL={1/κ2-ξ2for|ξ|<1/κiξ2-1/κ2for|ξ|1/κγT={1-ξ2for|ξ|<1iξ2-1for|ξ|1
(A12)
The corresponding traction vector is then 
t=η-B(ξ,x2)E(ξ,x2)B-1(ξ,x2)h(ξ)exp(-iξx1η)dξ
(A13)
The boundary conditions (A8) and (A7) at x2=0 lead to 
Δu(x1)=2-L(ξ)h(ξ)exp(-iξx1η)dξ=0for|x1|>1
(A14)
 
t(x1,0)=η-h(ξ)exp(-iξx1η)dξ=q(x1)for|x1|<1
(A15)
where 
L(ξ)=-iμ[4ξ2γLγT+(2ξ2-1)2][γT00γL]
(A16)
For convenience, let f(ξ)=L(ξ)h(ξ), then Eqs. (A14) and (A15) can be recast into a pair of dual integral equations. 
2-f(ξ)exp(-iξx1η)dξ=0for|x1|>1
(A17)
 
η-L-1(ξ)f(ξ)exp(-iξx1η)dξ=q(x1,0)for|x1|<1
(A18)
In component form 
2-fi(ξ)exp(-iξx1η)dξ=0for|x1|>1
(A19)
 
ημ-Kik(ξ)fk(ξ)exp(-iξx1η)dξ=qi(x1,0)for|x1|<1
(A20)
where K21(ξ)=K12(ξ)=0 and 
K11(ξ)=iγT[4ξ2γLγT+(2ξ2-1)2],K22(ξ)=iγL[4ξ2γLγT+(2ξ2-1)2]
(A21)
It is easy to verify that 
limξ[K11(ξ)/ξ]=limξ[K22(ξ)/ξ]=-2(1-1κ2)-m
(A22)
We can now divide the problem into mode I 
2-f1I(ξ)exp(-iξx1η)dξ=0for|x1|>1
(A23)
 
η-K11(ξ)f1I(ξ)exp(-iξx1η)dξ=0for|x1|<1
(A24)
 
2-f2I(ξ)exp(-iξx1η)dξ=0for|x1|>1
(A25)
 
ημ-K22(ξ)f2I(ξ)exp(-iξx1η)dξ=-iησ0exp(iηx1p1/κ)/κfor|x1|<1
(A26)
and mode II 
2-f1II(ξ)exp(-iξx1η)dξ=0for|x1|>1
(A27)
 
η-K11(ξ)f1II(ξ)exp(-iξx1η)dξ=-iητ0exp(iηx1p1/κ)/κfor|x1|<1
(A28)
 
2-f2II(ξ)exp(-iξx1η)dξ=0for|x1|>1
(A29)
 
η-K22(ξ)f2II(ξ)exp(-iξx1η)dξ=0for|x1|<1
(A30)
It is easily seen that f1I(ξ)=f2II(ξ)=0. Thus, we only need to solve for f1II(ξ) and f2I(ξ). Each of these two problems can be further split into parts that are even and odd in x1. For example, the mode I problem can be split into f2I(ξ)=f2Ie(ξ)+f2Io(ξ) so that f2Ie(ξ)=f2Ie(-ξ) and f2Io(ξ)=-f2Io(-ξ). Thus, for f2Ie(ξ), we have 
0f2Ie(ξ)cos(ξx1η)dξ=0for|x1|>1
(A31)
 
0K22(ξ)f2Ie(ξ)cos(ξx1η)dξ=-iσ02κμcos(ηx1p1/κ)for|x1|<1
(A32)
Similarly, for f2Io(ξ) 
0f2Io(ξ)sin(ξx1η)dξ=0for|x1|>1
(A33)
 
0K22(ξ)f2Io(ξ)sin(ξx1η)dξ=-iσ02κμsin(ηx1p1/κ)for|x1|<1
(A34)
For mode II, the equations are 
0f1IIe(ξ)cos(ξx1η)dξ=0for|x1|>1
(A35)
 
0K11(ξ)f1IIe(ξ)cos(ξx1η)dξ=-iτ02κcos(ηx1p1/κ)for|x1|<1
(A36)
 
0f1IIo(ξ)sin(ξx1η)dξ=0for|x1|>1
(A37)
 
0K11(ξ)f1IIo(ξ)sin(ξx1η)dξ=-iτ02κsin(ηx1p1/κ)for|x1|<1
(A38)
Since, 
cosx=πx2J-1/2(x),   sinx=πx2J1/2(x)
(A39)
It is not difficult to see that all the four pairs of dual integral equations are similar type. They are all special cases of the following more general dual integral equations: 
0G(ξ)Jν(xξ)dξ=0forx>d
(A40)
 
0ξαG(ξ)[1+R(ξ)]Jν(xξ)dt=p(x)for0x<d
(A41)
where limξR(ξ)=0. Equations (A40) and (A41) can be converted into a Fredholm integral equation of the second kind [54]. For 0<α<2, the corresponding Fredholm integral equation is 
g(x)+1π0dK(x,s)g(s)ds=x12-ν-α2F(x)for0x<d
(A42)
where 
K(x,s)=πxs0tR(t)Jν+α/2(xt)Jν+α/2(st)dtF(x)=0xp(t)tν+1(x2-t2)-1+α/2dt
(A43)
and G(x) is related to g(x) through 
G(ξ)=ξ1-α/22-1+α/2Γ(α/2)0dsg(s)Jν+α/2(sξ)ds
(A44)
Let us first consider Eqs. (A31) and (A32). By using Eq. (A39), we can rewrite them as 
0ξf2Ie(ξ)J-1/2(ξx1η)dξ=0for|x1|>1
(A45)
 
0ξξf2Ie(ξ)[1+R(ξ,η)]J-1/2(ξx1η)dξ=h(x1,η)for|x1|<1
(A46)
where 
R(ξ,η)=-K22(ξ)ξm-1,   m=2(1-1κ2),h(x1,η)=iσ0μκm2πx1ηcos(ηx1p1/κ)
(A47)
If we identify G(ξ)=ξf2Ie(ξ), α=1, ν=-1/2, and d=η in Eqs. (A45) and (A46), the corresponding Fredholm integral equation follows directly from Eq. (A42): 
g(x)+1π0ηK(x,s)g(s)ds=x12F(x)for0x<d
(A48)
where following Eqs. (A43) and (A44) 
K(x,s)=πxs0tR(t)J0(xt)J0(st)dt,G(x)=x1/22-1/2Γ(1/2)0dsg(s)J0(sx)ds
(A49)
 
F(x)=0xp(t)t1/2(x2-t2)-1/2dt
(A50)
The variables can be scaled by η, 
g(ηx)+1πη01K(ηx,ηs)g(ηs)ds=(ηx)12F(ηx)for0x<1
(A51)
where 
K(ηx,ηs)=πηxsK(x,s,η)F(ηx)=0xp(ηt)(ηt)1/2[(x)2-(t)2]-1/2dt
(A52)
with 
K(x,s,η)=0tR(t,η)J0(ηxt)J0(ηst)dt
(A53)
Once g(ηs) is solved from Eq. (A51), the unknown f2Ie(ξ) can be obtained through 
ξf2Ie(ξ)=G(ξ)=ξ1/2η2-1/2Γ(1/2)01ηsg(ηs)J0(ηsξ)ds
(A54)
Equation (A51) can be further simplified by introducing a new variable g(x)=g(ηx)/ηx 
g(x)+η201sK(x,s,η)g(s)ds=F(x,η)for0x<1
(A55)
where using p(ηt)=h(t), we obtain 
F(x,η)F(ηx)=0xp(ηt)(ηt)1/2[(x)2-(t)2]-1/2dt=iσ0μκmπ8J0(xηp1/κ)
(A56)
Analogously, Eq. (A54) can be converted to 
ξf2Ie(ξ)=η22ξπ01sg(s)J0(ηsξ)ds
(A57)
Finally, the COD is obtained from Eq. (A14) 
Δu1Ie(x1)=0
 
Δu2Ie(x1)=40ξf2Ie(ξ)ξ-1/2cos(ξx1η)dξ=4η22π01sg(s)ds0J0(ηsξ)cos(ξx1η)dξ=4η2πx11sg(s)s2-x12dsfor|x1|<1
(A58)
To obtain an asymptotic solution for g(x), we consider the asymptotic expansion of the kernel K(x,s,η) for small η. To this end, K(x,s,η) needs to be converted into integrals over finite intervals. This has been done in Refs. [11] and [53]. For s>x, 
K(x,s,η)=-im01[(2t2-κ2)2κ41-t2J0(ηxt/κ)H0(1)(ηst/κ)   +4t21-t2J0(ηxt)H0(1)(ηst)]dt
(A59)

Since K(x,s,η)=K(s,x,η), the expression for s<x can be obtained from Eq. (A59) by interchanging s and x.

Next, by using the asymptotic expansion of the Bessel functions for small argument 
J0(x)=n=0a2nx2n,   H0(1)(x)=(1+2iπlnx2)J0(x)+n=0b2nx2n
(A60)
One can expand K(x,s,η) into 
K(x,s,η)=lnηK1(x,s)+K2(x,s)+η2lnηK3(x,s)+O(η2)
(A61)
where for s>x 
K1(x,s)=m02mκ4
(A62)
 
K2(x,s)=m02mκ4lns-iπm04mκ4(1+2iγπ)+12m(1-m0κ4)lnκ   +18mκ4(7-8κ2+κ4)-m0mκ4ln2
(A63)
 
K3(x,s)=-(s2+x2)(5-6κ2+2κ4+κ6)16mκ6
(A64)

with m0=3-4κ2+3κ4 and γ the Euler–Mascheroni constant. For s<x, one only needs to interchange s and x in Eq. (A63).

Further, g(x) and F(x,η) can be expanded in like orders of η. 
g(x)=g0(x)+η2g2(x)+η2lnηg4(x)+O(η2)
(A65)
 
F(x,η)=iσ0μκmπ8[1+a2(xηp1/κ)2+a4(xηp1/κ)4+O(η4)]
(A66)
Finally, one can substitute all the expansions into the integral equation (A55). By grouping terms with the same order of η, a series of equations arrive 
g0(x)=iσ0μκmπ8,g2(x)=iσ0μκmπ8a2(xp1/κ)2-01sK2(x,s)g0(s)ds
(A67)
 
g4(x)=-01sK1(x,s,η)g0(s)ds
(A68)
Carrying out the integrals, we obtain 
g0(x)=iσ0μκmπ8
(A69)
 
g2(x)=-iσ04μκmπ8(xp1/κ)2-12iσ0μκmπ8{m04mκ4(x2-1)   -iπm04mκ4(1+2iγπ)+12m(1-m0κ4)lnκ   +18mκ4(7-8κ2+κ4)-m0mκ4ln2}
(A70)
 
g4(x)=-iσ0m04μm2κ5π8
(A71)
The asymptotic expansion of the corresponding COD is obtained by substituting the above expansions into Eq. (A58) 
Δu2Ie(x1)=Δu20Ie(x1)+η2Δu22Ie(x1)+η2lnηΔu24Ie(x1)   +O(η2)   for   |x1|<1
(A72)
where 
Δu20Ie(x1)=2iησ0μκm1-x12
(A73)
 
Δu22Ie(x1)=iησ0μκm1x12[p126κ2(1+2x12)+m06mκ4(1x12)+iπm04mκ4(1+2iγπ)12m(1m0κ4)lnκ18mκ4(78κ2+κ4)+m0mκ4ln2]
(A74)
 
Δu24Ie(x1)=-iησ0m02μm2κ51-x12
(A75)
Equation (A72) can be recast into 
Δu2Ie(x1)=Δu20Ie(x1)gn(ω,φ,x1)
(A76)
The function gn(ω,φ,x1) is given by 
gn(ω,φ,x1)=1-Anη2lnη+12[-p126κ2-(p123κ2+2A13)x12+iπAn+Bn]η2
(A77)
where 
An=m04mκ4,   Bn=1mκ4[m06-m0γ2-12(κ4-m0)lnκ   -18(7-8κ2+κ4)+m0ln2]
(A78)
It then follows from Eq. (15) that 
gn(ω,φ)=1-Anη2lnη-12(p124κ2+A16-iπAn-Bn)η2
(A79)
Going through a similar procedure, we can also obtain ΔujIo(x1), as well as ΔujIIe(x1) and ΔujIIo(x1). Without presenting the details, we list below the results. 
Δu2Io(x1)=ηΔu21Io(x1)+η3Δu23Io(x1)+O(η3)for|x1|<1
(A80)
where 
Δu21Io(x1)=iησ0p1x1μmκ21-x12
(A81)
 
Δu23Io(x1)=iησ0p148μmκ4[-2p12(2x13+x1)+m0mκ2(5x1-2x13)]1-x12
(A82)
Alternatively, one can write 
Δu2Io(x1)=Δu20Ie(x1)gn(ω,φ,x1)
(A83)
where 
gn(ω,φ,x1)=ηp1x12κ+η3p196κ3[-2p12(2x13+x1)   +m0mκ2(5x1-2x13)]
(A84)
Since the right hand side of Eq. (A84) is an odd function of x1, it follows from Eq. (15) that: 
gn(ω,φ)=0
(A85)
For the even part of mode II 
Δu2IIe(x1)=0
(A86)
 
Δu1IIe(x1)=Δu10IIe(x1)+η2Δu12IIe(x1)+η2lnηΔu14IIe(x1)   +O(η2)for|x1|<1
(A87)
where 
Δu10IIe(x1)=2iητ0κm1-x12
(A88)
 
Δu12IIe(x1)=iητ0κm1-x12[-p126κ2(1+2x12)+κ4+16mκ4(1-x12)   +iπ(κ4+1)4mκ4(1+2iγπ)+12mκ4lnκ+κ4+1mκ4ln2-1-κ48mκ4]
(A89)
 
Δu14IIe(x1)=-iητ0(κ4+1)2m2κ51-x12
(A90)
Thus, 
Δu1IIe(x1)=Δu10IIe(x1)gs(ω,φ,x1)
(A91)
where 
gs(ω,φ,x1)=1-Asη2lnη   +12[-p126κ2-(p123κ2+2As3)x12+iπAs+Bs]η2
(A92)
with 
As=κ4+14mκ4,   Bs=1mκ4[κ4+16-γ(κ4+1)2+12lnκ   +(κ4+1)ln2-18(1-κ4)]
(A93)
Consequently, it follows from Eq. (15) that 
gs(ω,φ)=1-Asη2lnη-12(p124κ2+As6-iπAs-Bs)η2
(A94)
Similarly, for the odd part of mode II 
Δu1IIo(x1)=ηΔu11IIo(x1)+η3Δu13IIo(x1)+O(η3)for|x1|<1
(A95)
where 
Δu11IIo(x1)=iητ0p1x1mκ21-x12
(A96)
 
Δu13IIo(x1)=iητ0p148mκ4[-2p12(2x13+x1)+κ4+1mκ2(5x1-2x13)]1-x12
(A97)
Therefore, 
Δu1IIo(x1)=Δu10IIe(x1)gs(ω,φ,x1)
(A98)
 
gs(ω,φ,x1)=ηp1x12κ+η3p196κ3[-2p12(2x13+x1)   +κ4+1mκ2(5x1-2x13)]
(A99)
Since the right hand side of Eq. (A99) is an odd function of x1, we immediately have 
gs(ω,φ)=0
(A100)
Summarizing all the four cases together, we arrive at 
gα(ω,φ)=1-Aαη2lnη-12(sin2φ4κ2+Aα6-iπAα-Bα)η2forα=n,s
(A101)

where Aα and Bα are constants given in Eq. (A78) for α=n and Eq. (A93) for α=s, respectively.

2

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

References

References
1.
Budiansky
,
B.
, and
O'connell
,
R. J.
,
1976
, “
Elastic Moduli of a Cracked Solid
,”
Int. J. Solids Struct.
,
12
(
2
), pp.
81
97
.10.1016/0020-7683(76)90044-5
2.
Horii
,
H.
, and
Nematnasser
,
S.
,
1983
, “
Overall Moduli of Solids With Microcracks—Load-Induced Anisotropy
,”
J. Mech. Phys. Solids
,
31
(
2
), pp.
155
171
.10.1016/0022-5096(83)90048-0
3.
Kachanov
,
M.
,
1980
, “
Continuum Model of Medium With Cracks
,”
ASCE J. Eng. Mech. Div.
,
106
(
5
), pp.
1039
1051
.
4.
Kachanov
,
M. L.
,
1982
, “
A Microcrack Model of Rock Inelasticity Part I: Frictional Sliding on Microcracks
,”
Mech. Mater.
,
1
(
1
), pp.
19
27
.10.1016/0167-6636(82)90021-7
5.
Zhao
,
Y.
,
Qiu
,
L.
,
Jacobs
,
J.
, and
Qu
,
J.
,
2015
, “
A Micromechanics Model for the Acoustic Nonlinearity Parameter in Solids With Distributed Microcracks
,”
Ultrasonics
(submitted).
6.
Walsh
,
J. B.
,
1969
, “
New Analysis of Attenuation in Partially Melted Rock
,”
J. Geophys. Res.
,
74
(
17
), pp.
4333
4337
.10.1029/JB074i017p04333
7.
Achenbach
,
J. D.
,
1982
,
Ray Methods for Waves in Elastic Solids: With Applications to Scattering by Cracks
,
A. K.
Gautesen
, and
H.
McMaken
, eds.,
Pitman Advanced Publishing
,
Boston
.
8.
Angel
,
Y.
, and
Achenbach
,
J.
,
1985
, “
Reflection and Transmission of Elastic Waves by a Periodic Array of Cracks
,”
ASME J. Appl. Mech.
,
52
(
1
), pp.
33
41
.10.1115/1.3169023
9.
Angel
,
Y.
, and
Achenbach
,
J.
,
1985
, “
Reflection and Transmission of Elastic Waves by a Periodic Array of Cracks: Oblique Incidence
,”
Wave Motion
,
7
(
4
), pp.
375
397
.10.1016/0165-2125(85)90006-X
10.
Mal
,
A.
,
1970
, “
Interaction of Elastic Waves With a Penny-Shaped Crack
,”
Int. J. Eng. Sci.
,
8
(
5
), pp.
381
388
.10.1016/0020-7225(70)90075-3
11.
Mal
,
A. K.
,
1970
, “
Interaction of Elastic Waves With a Griffith Crack
,”
Int. J. Eng. Sci.
,
8
(
9
), pp.
763
776
.10.1016/0020-7225(70)90003-0
12.
Qu
,
J.
,
1990
, “
Low Frequency Scattering by a Planer Crack
,”
Rev. Prog. QNDE
,
9
, pp.
61
68
.
13.
Eriksson
,
A. S.
,
Boström
,
A.
, and
Datta
,
S. K.
,
1995
, “
Ultrasonic Wave Propagation Through a Cracked Solid
,”
Wave Motion
,
22
(
3
), pp.
297
310
.10.1016/0165-2125(95)00036-I
14.
Gross
,
D.
, and
Zhang
,
C.
,
1992
, “
Wave Propagation in Damaged Solids
,”
Int. J. Solids Struct.
,
29
(
14
), pp.
1763
1779
.10.1016/0020-7683(92)90169-T
15.
Itou
,
S.
,
1978
, “
Three-Dimensional Wave Propagation in a Cracked Elastic Solid
,”
ASME J. Appl. Mech.
,
45
(
4
), pp.
807
811
.10.1115/1.3424423
16.
Kikuchi
,
M.
,
1981
, “
Dispersion and Attenuation of Elastic Waves Due to Multiple Scattering From Cracks
,”
Phys. Earth Planet. Inter.
,
27
(
2
), pp.
100
105
.10.1016/0031-9201(81)90037-6
17.
Kikuchi
,
M.
,
1981
, “
Dispersion and Attenuation of Elastic Waves Due To Multiple Scattering From Inclusions
,”
Phys. Earth Planet. Inter.
,
25
(
2
), pp.
159
162
.10.1016/0031-9201(81)90148-5
18.
Zhang
,
C.
, and
Gross
,
D.
,
1993
, “
Wave Attenuation and Dispersion in Randomly Cracked Solids—I. Slit Cracks
,”
Int. J. Eng. Sci.
,
31
(
6
), pp.
841
858
.10.1016/0020-7225(93)90097-E
19.
Zhang
,
C.
, and
Gross
,
D.
,
1993
, “
Wave Attenuation and Dispersion in Randomly Cracked Solids—II. Penny-Shaped Cracks
,”
Int. J. Eng. Sci.
,
31
(
6
), pp.
859
872
.10.1016/0020-7225(93)90098-F
20.
Smyshlyaev
,
V. P.
,
Willis
,
J. R.
, and
Sabina
,
F. J.
,
1993
, “
Self-Consistent Analysis of Waves in a Matrix-Inclusion Composite—III. A Matrix Containing Cracks
,”
J. Mech. Phys. Solids
,
41
(
12
), pp.
1809
1824
.10.1016/0022-5096(93)90071-M
21.
Smyshlyaev
,
V. P.
,
Willis
,
J. R.
, and
Sabina
,
F. J.
,
1993
, “
Self-Consistent Analysis of Waves in a Matrix-Inclusion Composite—II. Randomly Oriented Spheroidal Inclusions
,”
ASME J. Mech. Phys. Solids
,
41
(
10
), pp.
1589
1598
.10.1016/0022-5096(93)90015-8
22.
Sabina
,
F. J.
, and
Willis
,
J. R.
,
1993
, “
Self-Consistent Analysis of Waves in a Polycrystalline Medium
,”
Eur. J. Mech. A-Solids
,
12
(
2
), pp.
265
275
.
23.
Kim
,
J. Y.
,
1996
, “
Dynamic Self-Consistent Analysis for Elastic Wave Propagation in Fiber Reinforced Composites
,”
J. Acoust. Soc. Am.
,
100
(
4
), pp.
2002
2010
.10.1121/1.417910
24.
Richardson
,
J. M.
,
1979
, “
Harmonic Generation at an Unbonded Interface—I. Planar Interface Between Semi-Infinite Elastic Media
,”
Int. J. Eng. Sci.
,
17
(
1
), pp.
73
85
.10.1016/0020-7225(79)90008-9
25.
Biwa
,
S.
,
Hiraiwa
,
S.
, and
Matsumoto
,
E.
,
2006
, “
Experimental and Theoretical Study of Harmonic Generation at Contacting Interface
,”
Ultrasonics
,
44
(
S
), pp.
e1319
e1322
.10.1016/j.ultras.2006.05.010
26.
Biwa
,
S.
,
Nakajima
,
S.
, and
Ohno
,
N.
,
2004
, “
On the Acoustic Nonlinearity of Solid–Solid Contact With Pressure-Dependent Interface Stiffness
,”
ASME J. Appl. Mech.
,
71
(
4
), pp.
508
515
.10.1115/1.1767169
27.
Blanloeuil
,
P.
,
Meziane
,
A.
, and
Bacon
,
C.
,
2014
, “
Numerical Study of Nonlinear Interaction Between a Crack and Elastic Waves Under an Oblique Incidence
,”
Wave Motion
,
51
(
3
), pp.
425
437
.10.1016/j.wavemoti.2013.10.002
28.
Hirose
,
S.
, and
Achenbach
,
J. D.
,
1993
, “
Higher Harmonics in the Far Field Due to Dynamic Crack‐Face Contacting
,”
J. Acoust. Soc. Am.
,
93
(
1
), pp.
142
147
.10.1121/1.405651
29.
Kimoto
,
K.
, and
Ichikawa
,
Y.
,
2015
, “
A Finite Difference Method for Elastic Wave Scattering by a Planar Crack With Contacting Faces
,”
Wave Motion
,
52
, pp.
120
137
.10.1016/j.wavemoti.2014.09.007
30.
Len
,
K. S.
,
Severin
,
F.
, and
Solodov
,
I.
,
1991
, “
Experimental Observation of the Influence of Contact Nonlinearity on the Reflection of Bulk Acoustic Waves and the Propagation of Surface Acoustic Waves
,”
Sov. Phys. Acoust.
,
37
(
6
), pp.
610
612
.
31.
Nazarov
,
V. E.
, and
Sutin
,
A. M.
,
1997
, “
Nonlinear Elastic Constants of Solids With Cracks
,”
J. Acoust. Soc. Am.
,
102
(
6
), pp.
3349
3354
.10.1121/1.419577
32.
Pecorari
,
C.
,
2003
, “
Nonlinear Interaction of Plane Ultrasonic Waves With an Interface Between Rough Surfaces in Contact
,”
J. Acoust. Soc. Am.
,
113
(
6
), pp.
3065
3072
.10.1121/1.1570437
33.
Rudenko
,
O.
, and
Uv
,
C. A.
,
1994
, “
Nonlinear Acoustic Properties of a Rough Surface Contact and Acoustiodiagnostics of a Roughness Height Distribution
,”
Acoust. Phys.
,
40
(
4
), pp.
593
596
.
34.
Smyshlyaev
,
V.
, and
Willis
,
J.
,
1994
, “
Linear and Nonlinear Scattering of Elastic Waves by Microcracks
,”
J. Mech. Phys. Solids
,
42
(
4
), pp.
585
610
.10.1016/0022-5096(94)90053-1
35.
Smyshlyaev
,
V.
, and
Willis
,
J.
,
1996
, “
Effective Relations for Nonlinear Dynamics of Cracked Solids
,”
J. Mech. Phys. Solids
,
44
(
1
), pp.
49
75
.10.1016/0022-5096(95)00059-3
36.
Solodov
,
I. Y.
,
1988
, “
Acoustic Nonlinearity of a Piezoelectric Crystal Surface
,”
J. Appl. Phys.
,
64
(
6
), pp.
2901
2906
.10.1063/1.341574
37.
Chen
,
Z.
,
Tang
,
G.
,
Zhao
,
Y.
,
Jacobs
,
L. J.
, and
Qu
,
J.
,
2014
, “
Mixing of Collinear Plane Wave Pulses in Elastic Solids With Quadratic Nonlinearity
,”
J. Acoust. Soc. Am.
,
136
(
5
), pp.
2389
2404
.10.1121/1.4896567
38.
Jacobs
,
L.
,
Matlack
,
K.
,
Kim
,
J.-Y.
,
Qu
,
J.
, and
Joe
,
W. J.
,
2014
, “
Using Nonlinear Ultrasound to Measure Microstructural Changes Due to Radiation Damage in Steel
,”
J. Acoust. Soc. Am.
,
135
(
4
), pp.
2274
2274
.10.1121/1.4877452
39.
Matlack
,
K. H.
,
Kim
,
J. Y.
,
Wall
,
J. J.
,
Qu
,
J.
,
Jacobs
,
L. J.
, and
Sokolov
,
M. A.
,
2014
, “
Sensitivity of Ultrasonic Nonlinearity to Irradiated, Annealed, and Re-Irradiated Microstructure Changes in RPV Steels
,”
J. Nucl. Mater.
,
448
(
1–3
), pp.
26
32
.10.1016/j.jnucmat.2014.01.038
40.
Morlock
,
M. B.
,
Kim
,
J.-Y.
,
Jacobs
,
L. J.
, and
Qu
,
J.
,
2015
, “
Mixing of Two Co-Directional Rayleigh Surface Waves in a Nonlinear Elastic Material
,”
J. Acoust. Soc. Am.
,
137
(
1
), pp.
281
292
.10.1121/1.4904535
41.
Tang
,
G.
,
Liu
,
M.
,
Jacobs
,
L. J.
, and
Qu
,
J.
,
2014
, “
Detecting Localized Plastic Strain by a Scanning Collinear Wave Mixing Method
,”
J. Nondestr. Eval.
,
33
(
2
), pp.
196
204
.10.1007/s10921-014-0224-1
42.
Zeitvogel
,
D. T.
,
Matlack
,
K. H.
,
Kim
,
J.-Y.
,
Jacobs
,
L. J.
,
Singh
,
P. M.
, and
Qu
,
J.
,
2014
, “
Characterization of Stress Corrosion Cracking in Carbon Steel Using Nonlinear Rayleigh Surface Waves
,”
NDT & E Int.
,
62
, pp.
144
152
.10.1016/j.ndteint.2013.12.005
43.
Smyshlyaev
,
V. P.
, and
Willis
,
J. R.
,
1994
, “
Linear and Nonlinear Scattering of Elastic-Waves by Microcracks
,”
J. Mech. Phys. Solids
,
42
(
4
), pp.
585
610
.10.1016/0022-5096(94)90053-1
44.
Qu
,
J.
, and
Cherkaoui
,
M.
,
2006
,
Fundamentals of Micromechanics of Solids
,
Wiley
,
Hoboken, NJ
.
45.
Qu
,
J.
, and
Li
,
Q. Q.
,
1991
, ‘
Interfacial Dislocation and Its Applications to Interface Cracks in Anisotropic Bimaterials
,”
J. Elasticity
,
26
(
2
), pp.
169
195
.10.1007/BF00041220
46.
Qu
,
J.
, and
Bassani
,
J. L.
,
1993
, “
Interfacial Fracture Mechanics for Anisotropic Bimaterials
,”
ASME J. Appl. Mech.
,
60
(
2
), pp.
422
431
.10.1115/1.2900810
47.
Krajcinovic
,
D.
,
1996
,
Damage Mechanics
,
Elsevier
,
Dordrecht
.
48.
Angel
,
Y. C.
, and
Koba
,
Y. K.
,
1998
, “
Complex-Valued Wave Number, Reflection and Transmission in an Elastic Solid Containing a Cracked Slab Region
,”
Int. J. Solids Struct.
,
35
(
7–8
), pp.
573
592
.10.1016/S0020-7683(97)00068-1
49.
Nagy
,
P.
,
2015
, private communication.
50.
Zhao
,
Y.
,
Qiu
,
Y.
,
Jacobs
,
L. J.
, and
Qu
,
J.
,
2015
, “
Frequency-Dependent Tensile and Compressive Effective Moduli of Elastic Solids With Distributed Penny-Shaped Microcracks
,”
Acta Mech.
, (in press).
51.
Valle
,
C.
,
Qu
,
J. M.
, and
Jacobs
,
L. J.
,
1999
, “
Guided Circumferential Waves in Layered Cylinders
,”
Int. J. Eng. Sci.
,
37
(
11
), pp.
1369
1387
.10.1016/S0020-7225(98)00133-5
52.
Valle
,
C.
,
Qu
,
J. M.
, and
Jacobs
,
L. J.
,
2000
,
Scattering of Circumferential Waves in a Cracked Annulus, in Review of Progress in Quantitative Nondestructive Evaluation
, Vol.
19a and 19b
,
D. O.
Thompson
, and
D. E.
Chimenti
, eds.,
American Institute Physics
,
Melville, NY
, pp.
217
224
.
53.
Qu
,
J.
,
1987
,
Backscattering From Flaws in Fiber-Reinforced Composites, in Theoretical and Applied Mechanics
,
Northwestern University
,
Evanston, IL
.
54.
Noble
,
B.
,
1963
, “
The Solution of Bessel Function Dual Integral Equations by a Multiplying-Factor Method
,”
Proc. Cambridge Philos. Soc.
,
59
(
2
), pp.
351
362
.10.1017/S0305004100036987