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 2*a*. 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 $\nu $ be the Poisson's ratio of the uncracked solid, and $E\xaf\alpha $ and $\nu \xaf\alpha $ be, respectively, the corresponding effective Young's modulus and Poisson's ratio of the cracked solid under tension ($\alpha =t$) and under compression ($\alpha =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\xaf\alpha $ and $\nu \xaf\alpha $ 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 $-\pi /2\u2264\phi \u2264\pi /2$. Further, we assume that the solid is subjected to surface traction $p=\sigma \xaf\xb7m$ on its boundary $\u2202A$, where $\sigma \xaf$ is a constant stress tensor and **m** is the outward unit normal of $\u2202A$.

*A*without any crack is subjected to surface traction $p=\sigma \xaf\xb7m$ on its boundary $\u2202A$; and (2) a solid

*A*with

*N*randomly oriented cracks of length 2

*a*, where the surfaces of the cracks are subjected to applied normal and shear stresses $\sigma 0$ and $\tau 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*is the frictional force given in Eq. (2), and

Clearly, $\sigma 0$ is the normal stress to open the crack, $\tau 0$ is the shear stress in the direction of **s** to cause the crack faces to slid over each other, and $\mu \sigma 0H(-\sigma 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.

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-\nu 2)$ for plane strain. The dimensionless functions $gn(\omega ,\phi ,x)$ and $gs(\omega ,\phi ,x)$ are derived in the Appendix. It is noticed here that $gn(\omega ,\phi ,x)$ and $gs(\omega ,\phi ,x)$ are even functions of $\phi $, and are complex-valued. In the static limit of $\omega \u21920$, one has $gn(0,\phi ,x)=gs(0,\phi ,x)=1$. Thus, Eq. (5) reduces to the well-known result for a Griffith crack under static loading, e.g., Refs. [45] and [46].

*A*that contains the line segment 2

*a*, the line Dirac delta function satisfies the following:

## Dilute-Concentration Estimates

*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.,

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

*A*is thus,

**q**applied to the crack faces, i.e.,

*l*, Eq. (12) can be further simplified to

*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

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

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

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

*F*is related to the coefficient of friction through

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 $\mu =0$, and $F\u21920$ when $\mu \u2192\u221e$. 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., $\mu =0$), Eq. (34) yields $F=1$, meaning that all the cracks may participate in sliding. When the crack surfaces are extremely rough (i.e., $\mu \u2192\u221e$), Eq. (34) yields $F\u21920$. 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 $\mu \u2192\u221e$, one can see from Eq. (33) that $h1c\u21920$. Thus, Eqs. (36) and (37) immediately lead to $E\xaf'\u2009c=E'$ and $\nu \xaf'\u2009c=\nu '$, 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., $c\u226a1$. For example, if every square of side 4*a* contains on average one crack of length 2*a*, then $c\u22480.06$. Even if in each square of side 4*a*, there are on average five cracks of length 2*a*, one still has $c\u22480.3$. Five cracks of length 2*a* in a square of side 4*a* is a rather high crack density.

## Self-Consistent Estimates

where an overbar indicates that the quantity is evaluated based on the effective properties $E\xaf\alpha $ and $\nu \xaf\alpha $. For clarity, $g\xafn(\nu \xaf\alpha ,\omega ,\phi ,x)$ and $g\xafs(\nu \xaf\alpha ,\omega ,\phi ,x)$ are explicitly written as functions of $\nu \xaf\alpha $. Dimensional analysis indicates that they are independent of $E\xaf\alpha $.

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

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

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

## Wave Dispersion and Attenuation

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 $\eta 2$. It can then be shown through straightforward asymptotic analyses that the coefficient of attenuation $\alpha (\omega )$ introduced in Eq. (47) is proportional to the third power of frequency, i.e., $\alpha (\omega )\u221d\omega 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 $\alpha (\omega )\u221d\omega n+1$, where *n* is the dimension of the scatterers. Therefore, for two-dimensional cracks, i.e., *n* = 2, Rayleigh scattering leads to $\alpha (\omega )\u221d\omega 3$. In a separate publication [50], we will show that the scattering-induced attenuation by a distribution of penny-shaped cracks is scaled with $\alpha (\omega )\u221d\omega 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 $\alpha (\omega )\u221d\omega 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.

## 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\xd7L$ and it contains *N* microcracks. For numerical expediency, all the cracks have the same length of 2*a*. 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 $\rho =2700\u2003kg/m3$, $E=7\xd71010\u2003Pa$, and $v=0.33$, which yield a longitudinal phase velocity of $cL=6198\u2003m/s$. A linear elastic constitutive law is used in the FEM simulations.

In simulating the dynamic loading, a displacement excitation $u(x,t)=Asin(\omega t)$ is prescribed on the left edge of the $L\xd7L$ 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 $\Lambda =2\pi c/\omega $, and that of the second harmonic is $\Lambda /2$. In this simulation problem, there are four characteristic length parameters, namely, the sample size *L*, the wavelength $\Lambda $ 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, $\Lambda /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, $\omega =5\pi \u2003MHz$ is used, which corresponds to a wavelength of $\Lambda =2.48\u2003mm$. The half crack length used is $a=0.05\u2003mm$. These give us $\Lambda /a\u224850$, 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., $\Lambda /(2h)>20$. So, in our FEM model, the element size is kept at $h<\Lambda /40$, considering that the wavelength of the second harmonic is $\Lambda /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\xd7L$ 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 $\mu =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 $\mu $. As expected, increasing friction leads to higher compressive modulus. In the limit of $\mu \u2192\u221e$, 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 ($\eta =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\pi (tan\u221211/\mu \u2212\mu /\mu 2+1)$ is introduced, which varies between 0 (when $\mu \u2192\u221e$) and 1 (when $\mu \u21920$). 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

*E*and Poisson's ratio $\nu $ through

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

*a*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

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

Since $K\u2227(x,s,\eta )=K\u2227(s,x,\eta )$, the expression for $s<x$ can be obtained from Eq. (A59) by interchanging *s* and *x*.

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

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

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