Various types of density distribution functions of inclusion mass (is the Dirac delta function).
In many practical situations elastic composites are subjected to dynamic loadings of different physical nature, which origin the wave propagation in such structures. Then overall dynamic response of composite materials is characterized by the wave attenuation and dispersion due to the multiple wave scattering, in the local sense these materials are exhibited also by the dynamic stress intensification due to the wave interaction with the composite fillers. Essential influences on the mentioned phenomena have the shapes and the space distributions of inclusions, i.e. composite architecture, as well as matrix-inclusion materials characteristics. In this respect the numerical investigation of elastic wave propagation in the composite materials with inclusions of non-classical shape and contrast rigidity in comparison with the matrix material is highly demanded. A deep insight into their dynamic behavior, especially on the microscale, is extremely helpful to the design, optimization and manufacturing of composite materials with desired mechanical qualities, fracture and damage analysis, ultrasonic non-destructive testing of composites, and modeling of seismic processes in complex geological media.
The macroscopic dynamic properties of particulate elastic composites can be described by effective dynamic parameters of the equivalent homogeneous effective medium via a suitable homogenization procedure. Generally speaking, the homogenization procedure to determine the effective dynamic properties of particulate elastic composites is much more complicated than its static counterpart because of the inclusion interactions and multiple wave scattering effects. For small inclusion concentration or dilute inclusion distribution, their mutual interactions and the multiple wave scattering effects can be neglected approximately. In this case, the theory of Foldy , the quasi-crystalline approximation of Lax  and their generalizations to the elastic wave propagation [3-5] can be applied to determine the effective wave (phase) velocities and the attenuation coefficients in the composite materials with randomly distributed inclusions. In these models, wave scattering by a single inclusion has to be considered in the first step. Most previous publications on the subject have been focused on 3D elastic wave propagation analysis in composite materials consisting of an elastic matrix and spherical elastic inclusions (for example, see [6,7]). Aligned and randomly oriented ellipsoidal elastic inclusions have been considered in [8-10] under the assumption that the wavelength is sufficiently long compared to the dimensions of the individual inclusions (quasi-static limit). As special cases, the results for a random distribution of cracks and penny-shaped inclusions can be derived from those for ellipsoidal inclusions. In the long wavelength approximation, analytical solutions for a single inclusion as a series of the wave number have been presented in these works. However, this approach is applicable only for low frequencies or small wave numbers. For moderate and high frequencies, numerical methods such as the finite element method or the boundary element method can be applied. By using the boundary integral equation method (BIEM) or the boundary element method (BEM) in conjunction with Foldy’s theory the effective wave velocities and the wave attenuations in linear elastic materials with open and fluid-filled penny-shaped cracks as well as soft thin-walled circular inclusions have been calculated in [11,12]. Both aligned and randomly oriented defect configurations have been studied, where a macroscopic anisotropy for aligned cracks and non-spherical inclusions appears. Previous results have shown that distributed crack-like defects may cause a decrease in the phase velocity and an increase in the wave attenuation. The efficiency and the applicability ranges of 2D homogenization analysis of elastic wave propagation through a random array of scatters of different shapes and dilute concentrations based on the BEM and Foldy-type dispersion relations were demonstrated also by many authors, for instance, in the papers [13,14]. In 3D case this approach was applied for the numerical simulation of the average dynamic response of composite material containing rigid disk-shaped inclusions of equal mass only . Dynamic stresses near single inclusion of such type under time-harmonic and impulse elastic waves incidence where also investigated [16-18].
In this Chapter the effective medium concept is extended to the time-harmonic plane elastic wave propagation in an infinite linear elastic matrix with rigid disk-shaped movable inclusions of variable mass. Both time-harmonic plane longitudinal and transverse waves are considered in the analysis. The solution procedure consists of three steps. In the first step, the wave scattering problem is formulated as a system of boundary integral equations (BIEs) for the stress jumps across the inclusion surfaces. A BEM is developed to solve the BIEs numerically, where the kinetics of the inclusion and the “square-root” singularity of the stress jumps at the inclusion edge are taken into account properly. The improved regularization procedure for the obtained BIEs involving the analytical evaluation of regularizing integrals and results of mapping theory is elaborated to ensure the stable and correct numerical solution of the BIEs. The far-field scattering amplitudes of elastic waves induced by a single inclusion are calculated from the numerically computed stress jumps. In the second step, the simple Foldy-type approximation  is utilized to calculate the complex effective wave numbers for a dilute concentration of inclusions, where their interactions and multiple wave scattering can be neglected. The averages of the forward scattering amplitudes over 3D inclusion orientations or directions of the wave incidence and over inclusions masses are included into the resulting homogenization formula (dispersion relations). Finally, the effective wave velocity and the attenuation coefficient are obtained by taking the real and the imaginary parts of the effective wave numbers. To investigate the influence of the wave frequency on the effective dynamic parameters, representative numerical examples for longitudinal and transverse elastic waves in infinite elastic composite materials containing rigid disk-shaped inclusions with aligned and random orientation, as well as aligned, normal and uniform mass distribution are presented and discussed. Besides the global dynamic parameters, the mixed-mode dynamic stress intensity factors in the inclusion vicinities are calculated. They can be used for the fracture or cracking analysis of a composite.
2. Boundary integral formulation of 3D wave scattering problem for a single massive inclusion
Let us consider an elastic solid consisting of an infinite, homogeneous, isotropic and linearly elastic matrix specified by the mass density, the shear modulus and Poisson’s ratio, and a rigid disk-shaped inclusion with the mass
Here and hereafter the common factor is omitted, is the wave number of the incident wave, is the direction of propagation of the incident wave, is the angle characterizing the direction of the wave incidence, and is the polarization vector with and for the
By using the superposition principle, the total displacement field in the solid can be written in the form
where is the unknown displacement vector of the scattered wave, which satisfies the equations of motion and the radiation conditions at infinity (these well-known governing relations of elastodynamic theory can be found in ).
The inclusion is regarded as a rigid unit and its motion is described by the translation and the rotation with respect to the coordinate axes with the angles, and, respectively. Then the displacement components in the domain
In order to obtain the integral representations for the displacement components we apply the Betty-Rayleigh reciprocity theorem in conjunction with the properties of the elastodynamic fundamental solutions. As a result, the displacement components of the scattered waves can be written in the form :
where the displacement continuity conditions across the inclusion are used, is the distance between the field point and integration point, and are the jumps of the interfacial stresses across the inclusion, which are defined by
Eqs. (5) together with the equations of motion of the inclusion as a rigid unit yields the following relations between the translations and the rotations of the inclusion and the stress jumps:
where is the radius of inertia of the inclusion with respect to the -axis.
The displacement components in the matrix and the kinematical parameters of the inclusion are related to the stress jumps across the inclusion by the relations (4) and (6). Substitution of Eqs. (4) and (6) into Eqs. (3) results in three boundary integral equations (BIEs) for the stress jumps as
In Eq. (7), the kernels, and have the form
The problem governed by the BIEs (7) can be divided into an antisymmetric problem and a symmetric problem. The antisymmetric problem corresponding to the transverse motion of the inclusion is described by first equation of the BIEs (7) for the stress jump. After the solution of this equation the displacement and the rotations and can be obtained by using the relations (6). The symmetric problem corresponds to the motion of the inclusion in its own plane, which is governed by the last two equations of the BIEs (7) for the stress jumps and. After these quantities have been computed by solving these equations, the kinematical parameters and can be obtained by using the relations (6).
The kernels of the BIEs (7) contain weakly singular integrals only. To isolate these singularities explicitly we rewrite the BIEs (7) as
In Eq. (9), the last integrals on the left-hand sides exist in the ordinary sense. This fact follows from an analysis of the integrand in the limit. Therefore, in the numerical evaluation of these integrals it is sufficient to perform the integration over by excluding a small region (the neighborhood of the -point) around from
The singularities of the BIEs (9) are identical to those of the corresponding BIEs for the static inclusion problems, which have been investigated in  both for the antisymmetric and symmetric cases. The local behavior of the stress jumps at the front of the inclusion is also the same as in the static case. For a circular disk-shaped inclusion, the stress jumps have a “square-root” singularity, which can be expressed as
where are unknown smooth functions, and a is the radius of the inclusion.
Substitution of Eq. (10) into Eq. (9) results in a system of BIEs for the functions. These BIEs have a weak singularity at the source point and a “square-root” singularity at the edge of the inclusion. To regularize the singular BIEs, the following integral relations for the elastostatic kernels are utilized when:
Here the special integral identities, taken from , are used, namely:
Next we perform the following transformation of the variables:
where and are new variables in the rectangular domain. Equation (13) transforms the circular integration domain to a rectangular integration domain and eliminates the “square-root” singularity at the front of the inclusion corresponding to.
By applying Eqs. (11)-(13) to the BIEs (9) we obtain their regularized version as
In Eq. (15), and are defined by Eq. (13), is the normalized wave number of the
For the discretization of the domain, a boundary element mesh with equal-sized rectangular elements is used. For simplicity, constant elements are adopted in this analysis. By collocating the BIEs (14) at discrete points coinciding with the centroids of each element, a system of linear algebraic equations for discrete values of is obtained. After solving the system of linear algebraic equations numerically, the stress jumps across the inclusion can be obtained by the relations (10) and (15).
The far-field quantities of the scattered elastic waves can be computed from the stress jumps. For this purpose we use the asymptotic relations for an observation point far away from the inclusion, namely and, when. By substituting of these relations into the integral representation formula (4) and introducing the spherical coordinate system with the origin at the center of the inclusion as
the asymptotic expressions for the scattered radial and tangential displacements in the far-field are obtained in the form
Here, , , and are the longitudinal, vertically polarized transverse, and horizontally polarized transverse wave scattering amplitudes, respectively, which are related to the inclusion of normalized mass. They are given by
where are the coordinates of the spherical unit vectors, and.
The forward scattering amplitudes are defined as the values of in the direction of the wave incidence, i.e.,.
Thus, the scattering problem in the far-field is reduced to the numerical solution of the BIEs (14) and the subsequent computation of the scattering amplitudes by using Eq. (19), where the transformation or mapping relations (13) have to be considered.
For the convenient description of the wave parameters in the inclusion vicinity let us introduce the local coordinate system with the center in the inclusion contour point, so that the value corresponds to the inclusion plane, the axes and lie in the normal and tangential directions relative to the inclusion contour line, respectively, as depicted in Figure 1. Then the corresponding displacement and stress components at the arbitrary point
Here and are the polar coordinates of the point
By using the Eq. (20) the -factors can be defined directly from the stress jumps or the solutions of BIEs (7) by the following relations:
where the dependence of -factors on the inclusion mass also is fixed by the variable.
3. Dispersion relations for distributed inclusions of variable mass
We consider now a statistical distribution of rigid disk-shaped micro-inclusions in the matrix. The location of the micro-inclusions is assumed to be random, while their orientation is either completely random or aligned, see Figure 2. In the case of aligned inclusions, it is postulated that the inclusions are parallel to the -plane. The radius of the inclusions is assumed to be equal, while their masses can be variable due to the different material properties of inclusions and their geometric aspect ratios.
The average response of the composite materials to the wave propagation is characterized by the geometrical dispersion and attenuation of waves due to the wave scattering process. To describe these phenomena within the coherent wave field, the dynamic properties of the composite can be modeled by a complex and frequency-dependent wave number as
where is the effective phase velocity and is the attenuation coefficient for the wave of corresponding mode. With Eq. (22), the amplitude of a plane time-harmonic elastic wave propagating in the n-direction can be expressed as
For low concentration of inclusions or small number density, the interaction or multiple scattering effects among the inclusions can be neglected. Under these assumptions the complex effective wave numbers of plane - and -waves can be calculated by using the Foldy-type dispersion relation, which was extended to elastic waves in  and may be stated as
In Eq. (24), is the density parameter of the inclusions, corresponds to the number density of inclusions of the same radius, i.e. the number of inclusions per unit volume, is the average forward scattering amplitude of the corresponding wave mode by a single inclusion. For randomly oriented inclusions of variable mass the averages should be taken both over all possible inclusion orientations and masses. It should be noted here that the average over all inclusion orientations is the same as the average over all directions of the wave incidence (to avoid the additional average over all wave polarizations for an incoming
and in the case of randomly oriented inclusion they become the form
Here is the angle characterizing the direction of the wave incidence, the parameters and characterize the minimal and maximal masses, respectively, in the system of distributed inclusions with the density function of inclusion mass, and are the forward scattering amplitudes given by Eq. (19). The density distribution function of inclusion mass should satisfy the normalization condition
A suitable set of inclusion mass variations, which corresponds to aligned, normal and uniform distributions, is defined in Table 1.
The approximation for the complex wave number (24) can be considered as a special case of the multiple wave scattering models of higher orders [4,5], and it involves only the first order in the inclusion density and is thus only valid for a dilute or small inclusion density. In the case of a large density or high concentration of inclusions, more sophisticated models such as the self-consistent approach or the multiple scattering models should be applied, to take the mutual dynamic interactions between individual inclusions into account.
Once the complex effective wave numbers have been determined via Eq. (24), the effective wave velocities and the attenuation coefficients of the plane - and -waves can be obtained by considering the definition (22). This results in
It should be remarked here that Foldy’s theory was derived for isotropic wave scattering, which is appropriate macroscopically for the configuration of randomly oriented inclusions. A composite solid with aligned (parallel) disk-shaped inclusions exhibits a macroscopic anisotropy, namely a transversal anisotropy. When an incident plane wave propagates in an arbitrary direction, this gives rise to a coupling between the
4. Numerical analysis of global dynamic parameters of a composite
The method presented in the previous sections is used to calculate the effective dynamic parameters of a composite elastic solid with both parallel and randomly oriented rigid disk-shaped inclusions of variable mass for the propagation of time-harmonic plane
For comparison purpose, normalized effective wave velocities and normalized attenuation coefficients are introduced as, where the subscript stands for
For parallel or aligned disk-shaped inclusions, the macroscopic dynamic behavior of the composite materials is transversely isotropic. Thus, the effective wave velocities and the attenuation coefficient are dependent on the direction of the wave incidence. In this analysis, only two wave incidence directions are considered, namely normal incidence of
For normal incidence of a plane
In the low frequency range, the normalized attenuation coefficient increases rapidly with increasing, after reaching a maximum it then decreases and approaches its high frequency limit (see Figure 3(a)). The peak value of increases and is shifted to a smaller value of the dimensionless frequency with changing of inclusion mass distribution from uniform to normal and then to aligned. The normalized attenuation coefficient in the high-frequency or short-wave limit does not depend on the frequency and the inclusion mass. Also, the inclusions are unmovable in the high-frequency limit and the normalized attenuation coefficient can be obtained by using the Kirchhoff approximation for short waves . At low frequencies, the normalized effective wave velocity in the composite is smaller than that in the homogeneous matrix material (see Figure 3(b)). Then the normal inclusion mass distribution is characterized by the bigger (smaller) value of in comparison with the aligned (uniform) situation. An opposite tendency is observed in the range of higher frequencies. In addition, the normalized effective wave velocity in the composite can be bigger than that in the homogeneous matrix material, for instance for the considered inclusions of aligned mass and. The high-frequency limit at means that the velocity of the -wave in the short-wave limit coincides with that in the matrix. The explanation of this high-frequency limit follows from the geometrical optical interpretation of the wave field. The wave field at high frequencies may be considered as a set of independent beams propagating through the medium. Because of the existing continuous matrix material, the effective wave velocity should coincide with the wave velocity in the matrix in the high-frequency limit.
The corresponding numerical results for grazing incidence of a plane
For normal incidence of a plane
Next numerical examples concern the randomly oriented micro-inclusions, when the macroscopic dynamic behavior of the composite material is isotropic. It means that the effective wave velocity and the attenuation coefficient do not depend on the direction of the wave incidence. Both the translations and the rotations of the inclusions are exhibited in this case.
For an incident plane
Figure 7 demonstrates the corresponding results for the normalized attenuation coefficient and the normalized effective wave velocity for an incident plane
5. Numerical analysis of local dynamic parameters of a composite
Description of macroscopic dynamic response of a composite to elastic wave propagation by Eqs. (23) and (28) allows us the extension of analysis on the near-field quantities connected with each inclusion. Special attention should be paid to the dynamic stress intensity factors as the most important fracture parameters. Taking in mind the assumptions of neglecting the inclusions interaction, the relations (21) are applied for the estimation of the mode-I, II, and III dynamic stress intensity factors and in the vicinity of separate inclusion. We suppose the impinge on the inclusion of plane longitudinal
Attenuation and dispersion of time-harmonic elastic waves, as well as dynamic stress concentration in 3D composite materials consisting of a linear elastic matrix and rigid disk-shaped inclusions of variable mass is simulated numerically. Translations and rotations of the inclusions in the matrix are taken into account in the analysis. Wave scattering by a single disk-shaped inclusion is investigated by a boundary element method to obtain the stress jumps across the inclusion surfaces. Then, far-field scattering amplitudes of elastic waves are computed by using the stress jumps. To describe the average macroscopic dynamic properties of the composite materials with a random distribution of disk-shaped micro-inclusions, complex wave numbers are computed by the Foldy-type dispersion relations, from which the effective wave velocities and the wave attenuation can be obtained. The present analysis concerns a dilute distribution of micro-inclusions, when the mutual inclusion interactions and the multiple scattering effects are approximately neglected. Numerical examples involve:
both longitudinal and transversal waves propagation in a composite material;
parallel and randomly oriented rigid disk-shaped inclusions;
aligned, normal and uniform distributions of inclusion mass;
frequency-domain analysis of global dynamic parameters, such as the wave attenuation coefficients and effective wave velocities;
frequency-domain analysis of local dynamic parameters, such as the dynamic stress intensity factors in the inclusion vicinities.
As shown, particular dynamic properties of composite materials can be varied by controlled changes in the microstructure.
This work is sponsored by the State Foundation for Fundamental Researches of Ukraine (Project No. 40.1/018), which is gratefully acknowledged.
Foldy LL 1945Multiple Scattering Theory of Waves. Physical Review 67 107 119
Lax M. 1952Multiple Scattering of Waves. Physical Review 85 621 629
Gubernatis J. E. Domany E. 1984Effects of Microstructure on the Speed and Attenuation of Elastic Waves in Porous Materials. Wave Motion 6 579 589
Martin PA 2006Multiple Scattering: Interaction of Time-harmonic Waves with N Obstacles. Cambridge: Cambridge University Press 437 p.
Conoir JM, Norris AN 2010Effective Wavenumbers and Reflection Coefficients for an Elastic Medium Containing Random Configurations of Cylindrical Scatterers. Wave Motion 47 183 197
Kerr FH 1992The Scattering of a Plane Elastic Wave by Spherical Elastic Inclusions. International Journal of Engineering Science 30 169 186
Kanaun SK, Levin VM 2007Propagation of Longitudinal Elastic waves in Composites with a Random Set of Spherical Inclusions. Archive of Applied Mechanics 77 627 651
Sabina FJ, Smyshlaev VP, Willis JR 1993Self Consistent Analysis of Waves in a Matrix-Inclusion Composite I. Randomly Oriented Spheroidal Inclusions. Journal of the Mechanics and Physics of Solids 41 1573 1588
Kanaun SK, Levin VM 2008Self-Consistent Methods for Composites. 2-WavePropagation in Heterogeneous Materials. Heidelberg: Springer 294 p.
Levin V. Markov M. Kanaun S. 2008Propagation of Long Elastic Waves in Porous Rocks with Crack-Like Inclusions. International Journal of Engineering Science 46 620 638
AS Eriksson Boström. A. Datta S. K. 1995Ultrasonic Wave Propagation through a Cracked Solid. Wave Motion 22 297 310
Zhang Ch. Gross D. 1998On Wave Propagation in Elastic Solids with Cracks. Southampton: Computational Mechanics Publications 248 p.
Sato H. Shindo Y. 2002Influence of Microstructure on Scattering of Plane Elastic Waves by a Distribution of Partially Debonded Elliptical Inclusions. Mechanics of Materials 34 401 409
Maurel A. Mercier J. F. Lund F. 2004Elastic Wave Propagation through a Random Array of Dislocations. Physical Review B 70: 024303 EOF
Mykhas’kiv V. V. Khay O. M. Zhang Ch. Boström A. 2010Effective Dynamic Properties of 3D Composite Materials Containing Rigid Penny-Shaped Inclusions. Waves in Random and Complex Media 20 491 510
Kit HS, Mykhas’skiv VV, Khay OM 2002Analysis of the Steady Oscillations of a Plane Absolutely Rigid Inclusion in a Three-Dimensional Elastic Body by the Boundary Element Method. Journal of Applied Mathematics and Mechanics 66 817 824
Mykhas’kiv VV 2005Transient Response of a Plane Rigid Inclusion to an Incident Wave in an Elastic Solid. Wave Motion 41 133 144
Mykhas’kiv VV, Khay OM 2009Interaction between Rigid-Disc Inclusion and Penny- Shaped Crack under Elastic Time-Harmonic Wave Incidence. International Journal of Solids and Structures 46 602 616
Achenbach JD 1973Wave Propagation in Elastic Solids. Amsterdam: North-Holland Publishing Company 425 p.
Khaj MV 1993Two-Dimensional Integral Equations of Newton Potential Type and Their Applications. Kyiv: Naukova Dumka 253 p. (in Russian).
Kassir MK, Sih GC 1968Some Three-Dimensional Inclusion Problems in Elasticity. International Journal of Solids and Structures 4 225 241