Various types of density distribution functions of inclusion mass (

## 1. Introduction

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 [1], the quasi-crystalline approximation of Lax [2] 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 [15]. 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 [1] 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* M*, which thickness is much smaller than the characteristic size of its middle-surface

*. The center of the Cartesian coordinate system*S

*-wave or transverse*L

*-wave with the frequency*T

Here and hereafter the common factor * L*-wave and

*-wave.*T

By using the superposition principle, the total displacement field

where

The inclusion is regarded as a rigid unit and its motion is described by the translation * S*can be represented by

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 [18]:

where the displacement continuity conditions across the inclusion are used,

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

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

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

The kernels of the BIEs (7) contain weakly singular integrals only. To isolate these singularities explicitly we rewrite the BIEs (7) as

where

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

The singularities of the BIEs (9) are identical to those of the corresponding BIEs for the static inclusion problems, which have been investigated in [20] 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

Substitution of Eq. (10) into Eq. (9) results in a system of BIEs for the functions

Here the special integral identities, taken from [20], are used, namely:

Next we perform the following transformation of the variables:

where

By applying Eqs. (11)-(13) to the BIEs (9) we obtain their regularized version as

where

In Eq. (15), * T*-wave,

For the discretization of the domain

The far-field quantities of the scattered elastic waves can be computed from the stress jumps

the asymptotic expressions for the scattered radial

Here,

where

The forward scattering amplitudes are defined as the values of

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 * P*near the inclusion in the plane

Here * P*,

By using the Eq. (20) the

where the dependence of

## 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

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

where

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

In Eq. (24), * T*-wave, we assume that the normal to the inclusions lie in the plane of incidence of the incoming

*- or*TV

*-wave). Hence, in the case of parallel inclusions the expressions for the average forward scattering amplitudes of corresponding wave mode are*TH

and in the case of randomly oriented inclusion they become the form

Here

A suitable set of inclusion mass variations, which corresponds to aligned, normal and uniform distributions, is defined in Table 1.

Distribution type | |

Aligned | |

Normal | |

Uniform |

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

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 * L*- and

*-waves, and thus a change in the effective polarization vector. However, it is reasonable to apply Foldy’s theory, when the wave propagation is along the principal axes because of the decoupling of the*T

*- and*L

*-waves. In this special case, wave propagation can be treated like in the isotropic case.*T

## 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 * L*- or

*-waves. For numerical discretization of the inclusion surface, the domain*TV

For comparison purpose, normalized effective wave velocities and normalized attenuation coefficients are introduced as* L*- and

*-waves, respectively.*TV

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 * L*- and

*-waves and grazing incidence*T

*- wave. This choice provides the vanishing of the dynamic torque on the inclusions and, therefore, their zero-rotations.*L

For normal incidence of a plane * L*-wave, the normalized attenuation coefficient

In the low frequency range, the normalized attenuation coefficient

The corresponding numerical results for grazing incidence of a plane * L*-wave are presented in Figure 4. As followed from Figure 4(a), the normalized attenuation coefficient

*-wave is also very similar to Figure 3(b). Compared to the normalized effective wave velocity*L

*-wave, the maximum effective wave velocity*L

For normal incidence of a plane * T*-wave (then

*- and*TV

*-wave are the same), the numerical results for the normalized attenuation coefficient*TH

*-wave as presented in Figure 4. In comparison to the peak values of*L

*-wave, the peak values of*L

*-wave are increased, what follows from Figures 3(a), 4(a) and 5(a). Comparison of Figures 3(b), 4(b) and 5(b) shows, that the minimum values of the normalized effective wave velocity*T

*-wave for the same inclusion mass distribution, while the maximum values of*L

*-wave.*L

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 * L*-wave, the normalized attenuation coefficient

Figure 7 demonstrates the corresponding results for the normalized attenuation coefficient * TV*-wave. In contrast to parallel inclusions (see Figure 5), now the variations of

## 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 * L*-wave with constant amplitude

At normal * L-*wave incidence on the inclusion (antisymmetric problem)

At grazing * L-*wave incidence on the inclusion (symmetric problem)

## 6. Conclusion

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.

## Acknowledgement

This work is sponsored by the State Foundation for Fundamental Researches of Ukraine (Project No. 40.1/018), which is gratefully acknowledged.