Numerical Simulation of Wave Propagation in 3D Elastic Composites with Rigid Disk-Shaped Inclusions of Variable 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$$\rho$$, the shear modulus $$G$$ and Poisson’s ratio$$\nu$$, and a rigid disk-shaped inclusion with the mass M, which thickness is much smaller than the characteristic size of its middle-surface S. The center of the Cartesian coordinate system $$O{x}_1{x}_2{x}_3$$ coincides with the mass center of the inclusion (see Figure 1), within the described geometrical assumptions the limit values $$x_3=\pm 0$$ correspond to the opposite interfaces between the matrix and the inclusion, where a welded contact is assumed. The stress-strain state in the solid is induced by harmonic in the time $$\text{t}$$ plane longitudinal L-wave or transverse T-wave with the frequency$$\omega$$, the constant amplitude$$U_0$$, the phase velocities $$c_L$$ and$$c_T$$, and the wave numbers $${\chi }_L=\omega /c_L$$ and$${\chi }_T=\omega /c_T$$, respectively. The displacement vector $$u^{\text{in}}=(u_1^{\text{in}},u_2^{\text{in}},u_3^{\text{in}})$$ of the incident wave is given by the relation

$$u^{\text{in}}(x)=U\exp \left[i\chi \left(n\cdot x\right)\right]$$E1

Here and hereafter the common factor $$\exp \left(-i\omega t\right)$$ is omitted, $$\chi$$is the wave number of the incident wave, $$n=(\sin {\theta }_0,0,\cos {\theta }_0)$$is the direction of propagation of the incident wave, $${\theta }_0$$is the angle characterizing the direction of the wave incidence, and $$U$$ is the polarization vector with $$\chi ={\chi }_L$$ and $$U=U_0\cdot n$$ for the L-wave and $$\chi ={\chi }_T,$$ $$U=U_0\cdot e$$ and $$n\cdot e=0$$ for the T-wave.

By using the superposition principle, the total displacement field $$u^{\text{tot}}$$ in the solid can be written in the form

$$u^{\text{tot}}(x)=u^{\text{in}}(x)+u(x),$$E2

where $$u=(u_1,u_2,u_3)$$ 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 [19]).

Figure 1.

Single disk-shaped inclusion subjected to an incident elastic wave.

The inclusion is regarded as a rigid unit and its motion is described by the translation $$u^0=(u_1^0,u_2^0,u_3^0)$$ and the rotation with respect to the coordinate axes with the angles$${\Omega }_1$$, $${\Omega }_2$$and$${\Omega }_3$$, respectively. Then the displacement components in the domain S can be represented by

$$u_j(x)=-u_j^{\text{in}}(x)+u_j^0+{(-1)}^j{\Omega }_3{x}_{3-j},j=1,2,$$E3

$$u_3(x)=-u_3^{\text{in}}(x)+u_3^0+{\Omega }_1{x}_2-{\Omega }_2{x}_1,x(x_1,x_2,\pm 0)\in S.$$E4

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

$$u_j\left(x\right)=\frac{1}{4\pi G}{\displaystyle {\iint }_S\{\Delta {\sigma }_j\left(\xi \right)\frac{\exp \left(i{\chi }_T\left|x-\xi \right|\right)}{\left|x-\xi \right|}-\frac{1}{{\chi }_T^2}\frac{\partial }{\partial x_j}[\Delta {\sigma }_1\left(\xi \right)\frac{\partial }{\partial x_1}+\Delta {\sigma }_2\left(\xi \right)\frac{\partial }{\partial x_2}+}$$E5

$$+\Delta {\sigma }_3\left(\xi \right)\frac{\partial }{\partial x_3}]\left[\frac{\exp \left(i{\chi }_L\left|x-\xi \right|\right)}{\left|x-\xi \right|}-\frac{\exp \left(i{\chi }_T\left|x-\xi \right|\right)}{\left|x-\xi \right|}\right]\}d{S}_{\xi },$$

$$j=1,2,3$$E6

where the displacement continuity conditions across the inclusion are used, $\left|x-\xi \right|$is the distance between the field point $$x=(x_1,x_2,x_3)$$ and integration point$\xi =({\xi }_1,{\xi }_2,0)$, and $$\Delta {\sigma }_j$$ $$(j=1,2,3)$$ are the jumps of the interfacial stresses across the inclusion, which are defined by

$$\Delta {\sigma }_j(x)={\sigma }_{j3}^{-}(x)-{\sigma }_{j3}^{+}(x),j=1,2,3,x\in S,{\sigma }_{j3}^{\pm }(x)=\underset{x_3\to \pm 0}{\lim }{\sigma }_{j3}(x).$$E7

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$$\Delta {\sigma }_j$$:

$$u_j^0=\frac{1}{{\omega }^{2}M}{\displaystyle {\iint }_S\Delta {\sigma }_j\left(\xi \right)d{S}_{\xi }},$$

$$j=1,2,3$$E8

$${\Omega }_j=\frac{{\left(-1\right)}^{j+1}}{{\omega }^{2}M{i}_j^2}{\displaystyle {\iint }_S{\xi }_{3-j}\Delta {\sigma }_3\left(\xi \right)d{S}_{\xi }},$$

$$j=1,2$$E9

$${\Omega }_3=\frac{1}{{\omega }^{2}M{i}_3^2}{\displaystyle {\iint }_S\left[{\xi }_1\Delta {\sigma }_2\left(\xi \right)-{\xi }_2\Delta {\sigma }_1\left(\xi \right)\right]d{S}_{\xi }}$$E10

where $$i_j$$ is the radius of inertia of the inclusion with respect to the $$x_j$$-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

$${\displaystyle {\iint }_S\Delta {\sigma }_3\left(\xi \right)R_3\left(x,\xi \right)d{S}_{\xi }}=-4\pi G{\chi }_T^2{u}_3^{\text{in}}\left(x\right),$$

$$x\in S$$E11

$${\displaystyle {\iint }_S[\Delta {\sigma }_j\left(\xi \right)R_j\left(x,\xi \right)+\Delta {\sigma }_{3-j}\left(\xi \right)R_{j\left(3-j\right)}\left(x,\xi \right)]d{S}_{\xi }}=-4\pi G{\chi }_T^2{u}_j^{\text{in}}\left(x\right),$$

$$j=1,2,$$

$$x\in S$$E12

In Eq. (7), the kernels$R_j$, $R_{12}$and $R_{21}$ have the form

$$R_j\left(x,\xi \right)=L_1\left(\left|x-\xi \right|\right)-\frac{{\left(x_j-{\xi }_j\right)}^2}{{\left|x-\xi \right|}^2}{L}_2\left(\left|x-\xi \right|\right)-\frac{4\pi \rho }{M}\left(1+\frac{{\xi }_{3-j}{x}_{3-j}}{i_3^2}\right),$$

$$j=1,2$$E13

$$R_3\left(x,\xi \right)=L_1\left(\left|x-\xi \right|\right)-\frac{4\pi \rho }{M}\left(1+\frac{{\xi }_1{x}_1}{i_2^2}+\frac{{\xi }_2{x}_2}{i_1^2}\right)$$E14

$$R_{kj}\left(x,\xi \right)=-\frac{\left(x_1-{\xi }_1\right)\left(x_2-{\xi }_2\right)}{{\left|x-\xi \right|}^2}{L}_2\left(\left|x-\xi \right|\right)+\frac{4\pi \rho }{M}\frac{{\xi }_k{x}_j}{i_3^2},$$$k,j=1,2,$

$$k\ne j$$E15

$$L_j\left(d\right)=l_{j1}\left(d\right)\frac{\exp \left(i{\chi }_{L}d\right)}{d^3}-l_{j2}\left(d\right)\frac{\exp \left(i{\chi }_{T}d\right)}{d^3},$$

$$j=1,2$$E16

$$l_{11}\left(d\right)=1-i{\chi }_{L}d,$$

$$l_{12}\left(d\right)=1-i{\chi }_{T}d-{\chi }_T^2{d}^2,$$E17

$$l_{21}\left(d\right)=3-3i{\chi }_{L}d-{\chi }_L^2{d}^2,$$

$$l_{22}\left(d\right)=3-3i{\chi }_{T}d-{\chi }_T^2{d}^2.$$E18

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$$\Delta {\sigma }_3$$. After the solution of this equation the displacement $$u_3^0$$ and the rotations $${\Omega }_1$$ and $${\Omega }_2$$ 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 $$\Delta {\sigma }_1$$ and$$\Delta {\sigma }_2$$. After these quantities have been computed by solving these equations, the kinematical parameters $$u_1^0,u_2^0$$ and $${\Omega }_3$$ 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

$$A{\displaystyle {\iint }_S\frac{\Delta {\sigma }_3\left(\xi \right)}{\left|x-\xi \right|}d{S}_{\xi }+}{\displaystyle {\iint }_S\Delta {\sigma }_3\left(\xi \right)[\frac{1}{{\chi }_T^2}{R}_3\left(x,\xi \right)-\frac{A}{\left|x-\xi \right|}]d{S}_{\xi }}=-4\pi G{u}_3^{\text{in}}\left(x\right),$$

$$x\in S$$E19

$${\displaystyle {\iint }_S\frac{\Delta {\sigma }_j\left(\xi \right)}{\left|x-\xi \right|}\left[A+B\frac{{\left(x_j-{\xi }_j\right)}^2}{{\left|x-\xi \right|}^2}\right]d{S}_{\xi }+}B{\displaystyle {\iint }_S\Delta {\sigma }_{3-j}\left(\xi \right)\frac{\left(x_1-{\xi }_1\right)\left(x_2-{\xi }_2\right)}{{\left|x-\xi \right|}^3}d{S}_{\xi }}+$$E20

$$+{\displaystyle {\iint }_S\{\Delta {\sigma }_j\left(\xi \right)\left[\frac{1}{{\chi }_T^2}{R}_j\left(x,\xi \right)-\frac{A}{\left|x-\xi \right|}-B\frac{{\left(x_j-{\xi }_j\right)}^2}{{\left|x-\xi \right|}^3}\right]+}\Delta {\sigma }_{3-j}\left(\xi \right)[\frac{1}{{\chi }_T^2}{R}_{j\left(3-j\right)}\left(x,\xi \right)\begin{array}{c}\\ -\\ \end{array}$$E21

$$-B\frac{\left(x_1-{\xi }_1\right)\left(x_2-{\xi }_2\right)}{{\left|x-\xi \right|}^3}]\}d{S}_{\xi }=-4\pi G{u}_j^{\text{in}}\left(x\right),$$

$$j=1,2,$$

$$x\in S$$E22

where

$$A=\left(3-4\nu \right)/\left[4\left(1-\nu \right)\right],$$

$$B=1/\left[4\left(1-\nu \right)\right]$$E23

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$$\xi \to x$$. Therefore, in the numerical evaluation of these integrals it is sufficient to perform the integration over $S_x^0$ by excluding a small region (the neighborhood of the $$x$$-point) around $$x$$ from 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

$$\Delta {\sigma }_j(x)=\frac{f_j(x)}{\sqrt{a^2-x_1^2-x_2^2}},j=1,2,3,x\in S,$$E24

where $$f_j(x)$$ 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$$f_j(x)$$. These BIEs have a weak singularity $$1/|x-\xi |$$ at the source point $$\xi =x$$ 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$$x\in \text{S}$$:

$${\displaystyle {\iint }_S\frac{f(\xi )}{\sqrt{a^2-{\xi }_1^2-{\xi }_2^2}\cdot |x-\xi |}d{S}_{\xi }}={\pi }^{2}f(x)+{\displaystyle {\iint }_{S_x^0}\frac{f(\xi )-f(x)}{\sqrt{a^2-{\xi }_1^2-{\xi }_2^2}\cdot |x-\xi |}d{S}_{\xi }},$$E25

$${\displaystyle {\iint }_S\frac{f(\xi )}{\sqrt{a^2-{\xi }_1^2-{\xi }_2^2}}\cdot \frac{{(x_j-{\xi }_j)}^2}{|x-\xi {|}^3}d{S}_{\xi }}=\frac{{\pi }^2}{2}f(x)+{\displaystyle {\iint }_{S_x^0}\frac{\left[f(\xi )-f(x)\right]}{\sqrt{a^2-{\xi }_1^2-{\xi }_2^2}}\cdot \frac{{(x_j-{\xi }_j)}^2}{|x-\xi {|}^3}d{S}_{\xi }},$$E26

$${\displaystyle {\iint }_S\frac{f(\xi )}{\sqrt{a^2-{\xi }_1^2-{\xi }_2^2}}\cdot \frac{(x_1-\xi )(x_2-{\xi }_2)}{|x-\xi {|}^3}d{S}_{\xi }}={\displaystyle {\iint }_{S_x^0}\frac{\left[f(\xi )-f(x)\right]}{\sqrt{a^2-{\xi }_1^2-{\xi }_2^2}}\cdot \frac{(x_1-{\xi }_1)(x_2-{\xi }_2)}{|x-\xi {|}^3}d{S}_{\xi }.}$$E27

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

$${\displaystyle {\iint }_S\frac{d{S}_{\xi }}{\sqrt{a^2-{\xi }_1^2-{\xi }_2^2}|x-\xi |}}={\pi }^2,$$E28

$${\displaystyle {\iint }_S\frac{{(x_1-{\xi }_1)}^k{(x_2-{\xi }_2)}^j}{\sqrt{a^2-{\xi }_1^2-{\xi }_2^2}|x-\xi {|}^3}}d{S}_{\xi }=\{\begin{array}{cc}{\pi }^2/2,& \text{when}k=2,j=0\text{or}k=0,j=2,\\ 0,& \text{when}k=1,j=1.\end{array}$$E29

Next we perform the following transformation of the variables:

$$\{\begin{array}{c}{x}_1=a\sin y_1\cos y_2,\\ x_2=a\sin y_1\sin y_2,\end{array}$$

$$\{\begin{array}{l}{\xi }_1=a\sin {\eta }_1\cos {\eta }_2,\\ {\xi }_2=a\sin {\eta }_1\sin {\eta }_2,\end{array}$$E30

where $$y=(y_1,y_2)$$ and $$\eta =({\eta }_1,{\eta }_2)$$ are new variables in the rectangular domain$$\tilde{S}:\left\{0\le y_1,{\eta }_1\le \pi /2;0\le y_2,{\eta }_2\le 2\pi \right\}$$. 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$${\eta }_1=\pi /2$$.

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

$$A{\tilde{f}}_3\left(y\right)\left[{\pi }^2-{\displaystyle {\iint }_{{\tilde{S}}_y^0}\frac{\sin {\eta }_1}{R\left(y,\eta \right)}d{S}_{\eta }}\right]+\frac{1}{{\vartheta }^2}{\displaystyle {\iint }_{{\tilde{S}}_y^0}{\tilde{f}}_3\left(\eta \right){\tilde{R}}_3\left(y,\eta \right)\sin {\eta }_{1}d{S}_{\eta }}=-4\pi G{\tilde{u}}_3^{\text{in}}\left(y\right)$$

$$y\left(y_1,y_2\right)\in \tilde{S}$$E31

,

$${\tilde{f}}_j\left(y\right)\{A\left[{\pi }^2-{\displaystyle {\iint }_{{\tilde{S}}_y^0}\frac{\sin {\eta }_1}{R\left(y,\eta \right)}d{S}_{\eta }}\right]+B\left[\frac{{\pi }^2}{2}-{\displaystyle {\iint }_{{\tilde{S}}_y^0}{\Phi }_j\left(y,\eta \right)\sin {\eta }_{1}d{S}_{\eta }}\right]\}-$$E32

$$-B{\tilde{f}}_{3-j}\left(y\right){\displaystyle {\iint }_{{\tilde{S}}_y^0}\Psi \left(y,\eta \right)\sin {\eta }_1}d{S}_{\eta }+\frac{1}{{\vartheta }^2}{\displaystyle {\iint }_{{\tilde{S}}_y^0}[{\tilde{f}}_j\left(\eta \right){\tilde{R}}_j\left(y,\eta \right)+}$$E33

$$+\tilde{f}{}_{3-j}\left(\eta \right){\tilde{R}}_{j\left(3-j\right)}\left(y,\eta \right)]\sin {\eta }_{1}d{S}_{\eta }=-4\pi G{\tilde{u}}_j^{\text{in}}\left(y\right),$$

$$j=1,2,$$

$$y\left(y_1,y_2\right)\in \tilde{S}$$E34

,

where

$${\tilde{f}}_j\left(y\right)=f_j\left(x\right),$$

$${\tilde{u}}_j^{\text{in}}\left(y\right)=u_j^{\text{in}}\left(x\right),$$

$${\tilde{R}}_j\left(y,\eta \right)=a^3{R}_j\left(x,\xi \right),j=1,2,3,$$E35

$${\tilde{R}}_{kj}\left(y,\eta \right)=a^3{R}_{kj}\left(x,\xi \right)$$

$$k,j=1,2,$$

$$k\ne j$$E36

.

In Eq. (15), $x$and $\xi$ are defined by Eq. (13), $\vartheta ={\chi }_{T}a$is the normalized wave number of the T-wave, $${\tilde{S}}_y^0$$is the mapping of the domain $S_x^0$ due to the transformation (13) (in the domain $${\tilde{S}}_y^0$$ the points $y$ and $$\eta$$ do not coincide), and

$$R\left(y,\eta \right)={\left[{\sin }^2{y}_1+{\sin }^2{\eta }_1-2\sin y_1\sin {\eta }_1\cos \left(y_2-{\eta }_2\right)\right]}^{1/2}$$E37

$${\Phi }_1\left(y,\eta \right)=(\sin y_1\cos y_2-{\sin {\eta }_1\cos {\eta }_2)}^2{\left[R\left(y,\eta \right)\right]}^{-3}$$E38

$${\Phi }_2\left(y,\eta \right)=(\sin y_1\sin y_2-{\sin {\eta }_1\sin {\eta }_2)}^2{\left[R\left(y,\eta \right)\right]}^{-3}$$E39

$$\Psi \left(y,\eta \right)=(\sin y_1\cos y_2-\sin {\eta }_1\cos {\eta }_2)(\sin y_1\sin y_2-\sin {\eta }_1\sin {\eta }_2){\left[R\left(y,\eta \right)\right]}^{-3}$$E40

For the discretization of the domain$\tilde{S}$, 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 $${\tilde{f}}_j$$ is obtained. After solving the system of linear algebraic equations numerically, the stress jumps $$\Delta {\sigma }_j$$ 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$$\Delta {\sigma }_j$$. For this purpose we use the asymptotic relations for an observation point far away from the inclusion, namely $\left|x-\xi \right|\approx \left|x\right|-\left(x\cdot \xi \right)/\left|x\right|$ and${\left|x-\xi \right|}^{-1}\approx {\left|x\right|}^{-1}$, when$\left|x\right|\to \infty$. 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

$$x_1=R\sin \theta \cos \phi ,$$

$$x_2=R\sin \theta \sin \phi ,$$

$$x_3=R\cos \theta$$

$$0\le \theta ,\phi \le 2\pi$$E41

the asymptotic expressions for the scattered radial $u_R$ and tangential $u_{\theta }^{},u_{\phi }^{}$ displacements in the far-field are obtained in the form

$$u_R\left(R,\theta ,\phi \right)=\frac{\text{exp}\left(i{\chi }_{L}R\right)}{4\pi R}{F}_L\left(\theta ,\phi ,{\gamma }_0\right),$$

$$R\to \infty$$E42

$$u_{\theta }^{}\left(R,\theta ,\phi \right)=\frac{\text{exp}\left(i{\chi }_{T}R\right)}{4\pi R}{F}_{TV}\left(\theta ,\phi ,{\gamma }_0\right),$$

$$R\to \infty$$E43

$$u_{\phi }^{}\left(R,\theta ,\phi \right)=\frac{\text{exp}\left(i{\chi }_{T}R\right)}{4\pi R}{F}_{TH}\left(\theta ,\phi ,{\gamma }_0\right),$$

$$R\to \infty$$E44

Here, $F_L$, $F_{TV}$, and $F_{TH}$ are the longitudinal, vertically polarized transverse, and horizontally polarized transverse wave scattering amplitudes, respectively, which are related to the inclusion of normalized mass${\gamma }_0=M/(\rho a^3)$. They are given by

$$F_L\left(\theta ,\phi ,{\gamma }_0\right)=\frac{1-2\nu }{\text{2(1-}\nu \text{)G}}{\displaystyle \sum _{j=1}^3{p}_j{\displaystyle {\iint }_S\Delta {\sigma }_j\left(\xi \right)\text{exp}\left[-i{\chi }_L\left(p\cdot \xi \right)\right]d{S}_{\xi }}}$$E45

$$F_{TV}\left(\theta ,\phi ,{\gamma }_0\right)=\frac{1}{G}{\displaystyle \sum _{j=1}^3{r}_j{\displaystyle {\iint }_S\Delta {\sigma }_j\left(\xi \right)\text{exp}\left[-i{\chi }_T\left(p\cdot \xi \right)\right]d{S}_{\xi }}}$$E46

$$F_{TH}\left(\theta ,\phi ,{\gamma }_0\right)=\frac{1}{G}{\displaystyle \sum _{j=1}^3{\tau }_j{\displaystyle {\iint }_S\Delta {\sigma }_j\left(\xi \right)\text{exp}\left[-i{\chi }_T\left(p\cdot \xi \right)\right]d{S}_{\xi }}}$$E47

where $$p_j,r_j,{\tau }_j(j=1,2,3)$$ are the coordinates of the spherical unit vectors$p=(\text{sin}\theta \text{cos}\phi ,\text{sin}\theta \text{sin}\phi ,\text{cos}\theta )$, $r=(\text{cos}\theta \text{cos}\phi ,\text{cos}\theta \text{sin}\phi ,\text{-sin}\theta )$and$\tau =(\text{-sin}\phi ,\text{cos}\phi ,0)$.

The forward scattering amplitudes are defined as the values of $F_Z\left(\theta ,\phi ,{\gamma }_0\right)(Z=L,TV,TH)$ in the direction of the wave incidence, i.e.,$F_Z\left({\theta }_0,0,{\gamma }_0\right)$.

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 $$Otqz$$ with the center in the inclusion contour point, so that the value $$z=0$$ corresponds to the inclusion plane, the axes $$Ot$$ and $$Oq$$ 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 P near the inclusion in the plane $$q=0$$ can be approximated as [21]:

$$u_z(r,\psi ,\varphi )=\frac{\sqrt{2r}}{2G}\sin \frac{\psi }{2}\left(1-{\sin }^2\frac{\psi }{2}\right)K_{\text{I}}(\varphi )+\frac{\sqrt{2r}}{2G}\cos \frac{\psi }{2}\left(2-4\nu +{\cos }^2\frac{\psi }{2}\right)K_{\text{II}}(\varphi )+O(r^{\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.})$$E48

$$u_t(r,\psi ,\varphi )=-\frac{\sqrt{2r}}{2G}\cos \frac{\psi }{2}\left(4-4\nu -{\cos }^2\frac{\psi }{2}\right)K_{\text{I}}(\varphi )-\frac{\sqrt{2r}}{2G}\sin \frac{\psi }{2}\left(1-{\sin }^2\frac{\psi }{2}\right)K_{\text{II}}(\varphi )+O(r^{\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.})$$E49

$$u_q(r,\psi ,\varphi )=\frac{\sqrt{2r}}{2G}\cos \frac{\psi }{2}{K}_{\text{III}}(\varphi )+O(r^{\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.})$$E50

$${\sigma }_{zz}(r,\psi ,\varphi )=-\frac{1}{\sqrt{2r}}\cos \frac{\psi }{2}\left(1-2\nu -\sin \frac{\psi }{2}\sin \frac{3\psi }{2}\right)K_{\text{I}}(\varphi )+$$E51

$$+\frac{1}{\sqrt{2r}}\sin \frac{\psi }{2}\left(2-2\nu -\cos \frac{\psi }{2}\cos \frac{3\psi }{2}\right)K_{\text{II}}(\varphi )+O(1)$$E52

$${\sigma }_{tt}(r,\psi ,\varphi )=\frac{1}{\sqrt{2r}}\cos \frac{\psi }{2}\left(3-2\nu -\sin \frac{\psi }{2}\sin \frac{3\psi }{2}\right)K_{\text{I}}(\varphi )+\frac{1}{\sqrt{2r}}\sin \frac{\psi }{2}\left(2\nu +\cos \frac{\psi }{2}\cos \frac{3\psi }{2}\right)K_{\text{II}}(\varphi )+O(1)$$E53

$${\sigma }_{qq}(r,\psi ,\varphi )=\frac{2\nu }{\sqrt{2r}}\cos \frac{\psi }{2}{K}_{\text{I}}(\varphi )+\frac{2\nu }{\sqrt{2r}}\sin \frac{\psi }{2}{K}_{\text{II}}(\varphi )+O(1)$$E54

$${\sigma }_{\text{zq}}(r,\psi ,\varphi )=\frac{1}{\sqrt{2r}}\sin \frac{\psi }{2}{K}_{\text{III}}(\varphi )+O(1)$$E55

$${\sigma }_{tq}(r,\psi ,\varphi )=-\frac{1}{\sqrt{2r}}\cos \frac{\psi }{2}{K}_{\text{III}}(\varphi )+O(1)$$E56

$${\sigma }_{zt}(r,\psi ,\varphi )=\frac{1}{\sqrt{2r}}\sin \frac{\psi }{2}\left(2-2\nu +\cos \frac{\psi }{2}\cos \frac{3\psi }{2}\right)K_{\text{I}}(\varphi )+$$E57

$$+\frac{1}{\sqrt{2r}}\cos \frac{\psi }{2}\left(1-2\nu +\sin \frac{\psi }{2}\sin \frac{3\psi }{2}\right)K_{\text{II}}(\varphi )+O(1)$$E58

Here $$r$$ and $$\psi$$ are the polar coordinates of the point P, $$\varphi$$is the angular coordinate of the inclusion contour point (see Figure 1), $$K_{\text{I}},K_{\text{II}}$$and $$K_{\text{III}}$$ are the mode-I, II, and III dynamic stress intensity factors in the inclusion vicinity.

By using the Eq. (20) the ${\text{K}}_{\text{I,II,III}}$-factors can be defined directly from the stress jumps $$\Delta {\sigma }_j$$ or the solutions of BIEs (7) by the following relations:

$$K_{\text{I}}(\varphi ,{\gamma }_0)=-\frac{1}{4(1-\nu )\sqrt{a}}{\left\{\sqrt{a^2-x_1^2-x_2^2}\left[\Delta {\sigma }_1(x)\cos \varphi +\Delta {\sigma }_2(x)\sin \varphi \right]\right\}|}_{\begin{array}{l}{x}_1=a\cos \varphi ;\\ x_2=a\sin \varphi \end{array}},$$E59

$$K_{\text{II}}(\varphi ,{\gamma }_0)=-\frac{1}{4(1-\nu )\sqrt{a}}{\left[\sqrt{a^2-x_1^2-x_2^2}\Delta {\sigma }_3(x)\right]|}_{\begin{array}{l}{x}_1=a\cos \varphi ;\\ x_2=a\sin \varphi \end{array}},$$E60

$$K_{\text{III}}(\varphi ,{\gamma }_0)=\frac{1}{4\sqrt{a}}{\left\{\sqrt{a^2-x_1^2-x_2^2}\left[\Delta {\sigma }_1(x)\sin \varphi -\Delta {\sigma }_2(x)\cos \varphi \right]\right\}|}_{\begin{array}{l}{x}_1=a\cos \varphi ;\\ x_2=a\sin \varphi \end{array}}$$E61

where the dependence of ${\text{K}}_{\text{I,II,III}}$-factors on the inclusion mass also is fixed by the variable$${\gamma }_0$$.

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 $$x_1-x_2$$-plane. The radius $$a$$ 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.

Figure 2.

Multiple disk-shaped inclusions in the elastic wave field: parallel inclusions with the angle of wave incidence (left); randomly oriented inclusions (right).

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 $$K\left(\omega \right)$$ as

$$K\left(\omega \right)=\frac{\omega }{c^{\ast }\left(\omega \right)}+i\alpha \left(\omega \right)$$E62

where $$c^{\ast }\left(\omega \right)$$ is the effective phase velocity and $$\alpha \left(\omega \right)$$ 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

$$u(x,\omega )=U\exp \left[iK\left(n\cdot x\right)\right]=U\exp \left[-\alpha (\omega )\left(n\cdot x\right)\right]\exp \left[i\omega \left(n\cdot x\right)/c^{*}(\omega )\right]$$E63

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 $$K_Z(Z=L,T)$$ of plane $$L$$- and $$T$$-waves can be calculated by using the Foldy-type dispersion relation, which was extended to elastic waves in [3] and may be stated as

$$K_Z^2={\chi }_Z^2+\epsilon a^{-3}{F}_Z^{*},Z=L,T$$E64

In Eq. (24), $$\epsilon$$is the density parameter of the inclusions, $$\epsilon a^{-3}$$corresponds to the number density of inclusions of the same radius, i.e. the number of inclusions per unit volume, $F_Z^{*}$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 T-wave, we assume that the normal to the inclusions lie in the plane of incidence of the incoming TV- or TH-wave). Hence, in the case of parallel inclusions the expressions for the average forward scattering amplitudes of corresponding wave mode are

$$F_Z^{*}\left({\theta }_0\right)={\displaystyle {\int }_{\tau }^{\beta }{F}_Z\left({\theta }_0,0,\gamma \right)g(\gamma )d\gamma },Z=L,TV,TH,$$E65

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

$$F_Z^{*}=\frac{1}{2}{\displaystyle {\int }_0^{\pi }{\displaystyle {\int }_{\tau }^{\beta }{F}_Z\left({\theta }_0,0,\gamma \right)g(\gamma )\sin {\theta }_{0}d\gamma d{\theta }_0}},Z=L,TV,TH.$$E66

Here $${\theta }_0$$ is the angle characterizing the direction of the wave incidence, the parameters $$\tau$$ and $$\beta$$ characterize the minimal and maximal masses, respectively, in the system of distributed inclusions with the density function $$g$$ of inclusion mass, and $F_Z({\theta }_0,0,\gamma )$ are the forward scattering amplitudes given by Eq. (19). The density distribution function of inclusion mass should satisfy the normalization condition

$${\displaystyle {\int }_{\tau }^{\beta }g(\gamma )d\gamma }=1$$E67

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

Distribution type $$g(\gamma )$$
Aligned	$$\delta (\gamma -{\gamma }_0)$$
Normal	$$3{\gamma }^2/({\beta }^3-{\tau }^3)$$
Uniform	$$1/(\beta -\tau )$$

Table 1.

Various types of density distribution functions of inclusion mass ($$\delta$$is the Dirac delta function).

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 $$K_Z$$ have been determined via Eq. (24), the effective wave velocities $$c_Z^{\ast }$$ and the attenuation coefficients $${\alpha }_Z$$ of the plane $$L$$- and $$T$$-waves can be obtained by considering the definition (22). This results in

$$c_Z^{\ast }\left(\omega \right)=\frac{\omega }{\mathrm{Re}\left[K_Z\left(\omega \right)\right]},$$

$${\alpha }_Z\left(\omega \right)=\mathrm{Im}\left[K_Z\left(\omega \right)\right]$$

$$Z=L,TV,TH$$E68

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 T-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 L- and T-waves. In this special case, wave propagation can be treated like in the isotropic case.

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 TV-waves. For numerical discretization of the inclusion surface, the domain $$\tilde{S}$$ is divided into 264 rectangular elements of length $$\Delta y_1=\pi /22$$ and $$\Delta y_2=\pi /12$$ in the $$y_1$$- and $$y_2$$-directions, the domain ${\tilde{S}}_y^0$ is chosen as$${\tilde{S}}_y^0=\tilde{S}\backslash \Delta {\tilde{S}}_y$$, where $$\Delta {\tilde{S}}_y$$ is the boundary element with the point y in the center of the element. Poisson’s ratio is selected as 0.3, the radii of inertia of the circular inclusion are defined as$$i_1=i_2=a/2,i_3=a/\sqrt{2}$$. In the numerical examples, the radius $$a$$ of the inclusions is assumed to be equal, while their masses are varied by the distribution laws defined in Table 1 from the minimal value $$\tau =5$$ to the maximal value$$\beta =20$$. The inclusion density parameter is fixed as $$\epsilon =$$0.01 in accordance to the assumption of low concentration of inclusions in a matrix.

For comparison purpose, normalized effective wave velocities and normalized attenuation coefficients are introduced as$${\overline{c}}_Z=c_Z^{\ast }/c_Z$$, $${\overline{\alpha }}_Z=2a\cdot {\alpha }_Z/\left(\pi \epsilon \right),$$where the subscript $$Z=L,TV$$ stands for L- and TV-waves, respectively.

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 $$({\theta }_0=0^{\text{o}})$$ of L- and T-waves and grazing incidence $$({\theta }_0={90}^{\text{o}})$$ of L- wave. This choice provides the vanishing of the dynamic torque on the inclusions and, therefore, their zero-rotations.

For normal incidence of a plane L-wave, the normalized attenuation coefficient $${\overline{\alpha }}_L$$ and the normalized effective wave velocity $${\overline{c}}_L$$ versus the dimensionless wave number $$\vartheta$$ are presented in Figure 3. Three separate parametrical studies are performed to show the effects of different inclusion mass distributions: aligned with$${\gamma }_0=\beta =20$$, normal and uniform. To check the accuracy of the implemented BEM, the present numerical results are compared with the analytical low-frequency solutions given in [9]. Here and hereafter a very good agreement between both results is observed in the frequency range$$0\le \vartheta \le 0.3$$.

Figure 3.

Effects of the inclusion mass distribution on the normalized attenuation coefficient (a) and the normalized effective wave velocity (b) as the functions of dimensionless wave number for parallel inclusions and normal L-wave incidence.

In the low frequency range, the normalized attenuation coefficient $${\overline{\alpha }}_L$$ increases rapidly with increasing$$\vartheta$$, after reaching a maximum it then decreases and approaches its high frequency limit (see Figure 3(a)). The peak value of $${\overline{\alpha }}_L$$ 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 $${\overline{\alpha }}_L$$ 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 $${\overline{\alpha }}_L$$ can be obtained by using the Kirchhoff approximation for short waves [19]. At low frequencies, the normalized effective wave velocity $${\overline{c}}_L$$ 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 $${\overline{c}}_L$$ 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 $${\overline{c}}_L$$ in the composite can be bigger than that in the homogeneous matrix material, for instance for the considered inclusions of aligned mass and$$\vartheta >1$$. The high-frequency limit $${\overline{c}}_L\to 1$$ at $$\vartheta \to \infty$$ means that the velocity of the $$L$$-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 L-wave are presented in Figure 4. As followed from Figure 4(a), the normalized attenuation coefficient $${\overline{\alpha }}_L$$ shows a similar dependence on the dimensionless wave number $$\vartheta$$ and the inclusion mass distribution. For the same density function $$g$$ of inclusion mass, the peaks of $${\overline{\alpha }}_L$$ are larger than that for normal incidence as depicted in Figure 3(a), but the high-frequency limit of $${\overline{\alpha }}_L$$ in the case with $${\theta }_0={90}^{\text{o}}$$ is smaller than that in the case with$${\theta }_0=0^{\text{o}}$$. The normalized effective wave velocity $${\overline{c}}_L$$ in Figure 4(b) for grazing incidence of a plane L-wave is also very similar to Figure 3(b). Compared to the normalized effective wave velocity $${\overline{c}}_L$$ for a normal incidence of plane L-wave, the maximum effective wave velocity $${\overline{c}}_L$$ is increased while its minimum value is decreased. Also in contrast to Figure 3(b), the normalized effective wave velocity $${\overline{c}}_L$$ for all considered functions $$g$$ is larger than that of the matrix material at high frequencies. In general, the quite tangled effects of the inclusions’ stiffness and mass on the average phase velocity are observed at different frequencies.

Figure 4.

Effects of the inclusion mass distribution on the normalized attenuation coefficient (a) and the normalized effective wave velocity (b) as the functions of dimensionless wave number for parallel inclusions and grazing L-wave incidence.

For normal incidence of a plane T-wave (then TV- and TH-wave are the same), the numerical results for the normalized attenuation coefficient $${\overline{\alpha }}_T$$ and the normalized effective wave velocity $${\overline{c}}_T$$ are presented in Figure 5 versus the dimensionless wave number$$\vartheta$$. The global behavior of the normalized attenuation coefficient $${\overline{\alpha }}_T$$ and the normalized effective wave velocity $${\overline{c}}_T$$ is very similar to that for grazing incidence (i.e.$${\theta }_0={90}^{\text{o}}$$) of plane L-wave as presented in Figure 4. In comparison to the peak values of $${\overline{\alpha }}_L$$ for normal and grazing incidence of a plane L-wave, the peak values of $${\overline{\alpha }}_T$$ for an normal incidence of a plane T-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 $${\overline{c}}_T$$ are reduced compared to $${\overline{c}}_L$$ for normal and grazing incidence of a plane L-wave for the same inclusion mass distribution, while the maximum values of $${\overline{c}}_T$$ lie between the $${\overline{c}}_L$$-values for normal and grazing incidence of a plane L-wave.

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.

Figure 5.

Effects of the inclusion mass distribution on the normalized attenuation coefficient (a) and the normalized effective wave velocity (b) as the functions of dimensionless wave number for parallel inclusions and normal T-wave incidence.

For an incident plane L-wave, the normalized attenuation coefficient $${\overline{\alpha }}_L$$ and the normalized effective wave velocity $${\overline{c}}_L$$ are presented in Figure 6 versus the dimensionless wave number$$\vartheta$$. A comparison of Figure 6 with Figures 3 and 4 for parallel inclusions shows a similar dependence of the $${\overline{\alpha }}_L$$ and $${\overline{c}}_L$$ on the dimensionless wave number $$\vartheta$$ and the inclusion mass distribution. As expected, the peak values of $${\overline{\alpha }}_L$$ and $${\overline{c}}_L$$ for randomly oriented penny-shaped inclusions are larger than that for parallel inclusions with $${\theta }_0=0^{\text{o}}$$ but smaller than that for parallel inclusions with$${\theta }_0={90}^{\text{o}}$$.

Figure 6.

Effects of the inclusion mass distribution on the normalized attenuation coefficient (a) and the normalized effective wave velocity (b) as the functions of dimensionless wave number for randomly oriented inclusions and L-wave incidence.

Figure 7 demonstrates the corresponding results for the normalized attenuation coefficient $${\overline{\alpha }}_T$$ and the normalized effective wave velocity $${\overline{c}}_T$$ for an incident plane TV-wave. In contrast to parallel inclusions (see Figure 5), now the variations of $${\overline{\alpha }}_T$$ and $${\overline{c}}_T$$ with the dimensionless wave number $$\vartheta$$ are rather complicated at low frequencies. For instance, the normalized attenuation coefficient $${\overline{\alpha }}_T$$ for aligned distribution of inclusion mass in Figure 7(a) shows two distinct peaks, which are not observed in the case of parallel inclusions. The normalized attenuation coefficient $${\overline{\alpha }}_T$$ for randomly oriented inclusions is smaller than that for parallel inclusions with$${\theta }_0=0^{\text{o}}$$, while the normalized effective wave velocity $${\overline{c}}_T$$ could be larger or smaller than that for parallel inclusions depending on the dimensionless wave number $$\vartheta$$ and the density function $$g$$ of inclusion mass (see Figure 7(b)).

Figure 7.

Effects of the inclusion mass distribution on the normalized attenuation coefficient (a) and the normalized effective wave velocity (b) as the functions of dimensionless wave number for randomly oriented inclusions and TV-wave incidence.

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 $$K_{\text{I}},K_{\text{II}}$$ and $$K_{\text{III}}$$ in the vicinity of separate inclusion. We suppose the impinge on the inclusion of plane longitudinal L-wave with constant amplitude$$U{}_0$$, which should be corrected for each inclusion in accordance to its localization in a composite with the considering of attenuation and dispersion law (23). The value $$K_{*}=U_{0}G/\left[4(1-\nu )\sqrt{a}\right]$$ is chosen as the normalizing factor for the amplitudes of dynamic stress intensity factors, so that$${\overline{K}}_{\text{I,II,III}}=\left|K_{\text{I,II,III}}\right|/K_{*}$$. All material and discretization parameters are the same as it is fixed in the previous Section. Different inclusion masses are involved into analysis by the changing of parameter$${\gamma }_0$$.

At normal L-wave incidence on the inclusion (antisymmetric problem) $$K_{\text{I}}=K_{\text{III}}=0,$$and the mode-II dynamic stress intensity factor $$K_{\text{II}}$$ does not vary along the inclusion contour (see Figure 8). It follows from the Figure 8 that in the initial range of wave numbers $${\overline{K}}_{\text{II}}$$-factor rapidly increases from a zero value, what is more pronounced for the inclusions of large mass. A further increase in $$\vartheta$$ in the case of an inclusion with the mass characteristic $${\gamma }_0=20$$ leads to a local maximum of$${\overline{K}}_{\text{II}}$$. For higher wave numbers, a regularity is the approaching of $${\overline{K}}_{\text{II}}$$-factors for the inclusions of different mass, in addition, this approaching is from above as the mass of inclusion increases. Subsequently, a linear relationship between $${\overline{K}}_{\text{II}}$$ and $$\vartheta$$ is reached in accordance to the increasing order of stresses in an incident wave.

Figure 8.

Effects of the inclusion mass on the normalized amplitude of the mode-II dynamic stress intensity factor as the function of dimensionless wave number for normal L-wave incidence on an inclusion.

Figure 9.

Effects of the inclusion mass on the normalized amplitude of the mode-I dynamic stress intensity factor as the function of dimensionless wave number for grazing L-wave incidence on an inclusion: solid lines correspond to the contour point where the wave runs on the inclusion; marked lines correspond to the contour point where the wave runs down from the inclusion.

At grazing L-wave incidence on the inclusion (symmetric problem) $$K_{\text{II}}=0,$$and the mode-I and mode-III dynamic stress intensity factor $$K_{\text{I}}$$ and $$K_{\text{III}}$$ depend on the angular coordinate $$\varphi$$ of the inclusion contour point. Hence, in Figure 9 $${\overline{K}}_{\text{I}}$$-factor at the most representative points of the inclusion contour is showed, where $$K_{\text{III}}=0$$ due to the symmetry conditions. A similar frequency behaviour (with a difference in the quantitative values) as in the previous antisymmetric problem is observed here for $${\overline{K}}_{\text{I}}$$-factor at the point, where the wave runs on the inclusion. It should be mentioned that$${\overline{K}}_{\text{I}}$$-factor at this point exceed one at the point, where the wave runs down from the inclusion. In the range of high wave numbers $${\overline{K}}_{\text{I}}$$-factors for the inclusions of different mass approach also, but now from below as the mass of inclusion increases.

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.