Open access peer-reviewed chapter

Numerical Simulation of Wave Propagation in 3D Elastic Composites with Rigid Disk-Shaped Inclusions of Variable Mass

By Viktor Mykhas'kiv

Submitted: November 4th 2011Reviewed: May 4th 2012Published: August 22nd 2012

DOI: 10.5772/48113

Downloaded: 1555

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ρ, the shear modulus Gand Poisson’s ratioν, 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 Ox1x2x3coincides with the mass center of the inclusion (see Figure 1), within the described geometrical assumptions the limit values x3=±0correspond 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 tplane longitudinal L-wave or transverse T-wave with the frequencyω, the constant amplitudeU0, the phase velocities cLandcT, and the wave numbers χL=ω/cLandχT=ω/cT, respectively. The displacement vector uin=(u1in,u2in,u3in)of the incident wave is given by the relation

uin(x)=Uexp[iχ(nx)]E1

Here and hereafter the common factor exp(iωt)is omitted, χis the wave number of the incident wave, n=(sinθ0,0,cosθ0)is the direction of propagation of the incident wave, θ0is the angle characterizing the direction of the wave incidence, and Uis the polarization vector with χ=χLand U=U0nfor the L-wave and χ=χT,U=U0eand ne=0for the T-wave.

By using the superposition principle, the total displacement field utotin the solid can be written in the form

utot(x)=uin(x)+u(x),E2

where u=(u1,u2,u3)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 u0=(u10,u20,u30)and the rotation with respect to the coordinate axes with the anglesΩ1, Ω2andΩ3, respectively. Then the displacement components in the domain S can be represented by

uj(x)=ujin(x)+uj0+(1)jΩ3x3j,j=1,2,E3
u3(x)=u3in(x)+u30+Ω1x2Ω2x1,x(x1,x2,±0)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]:

uj(x)=14πGS{Δσj(ξ)exp(iχT|xξ|)|xξ|1χT2xj[Δσ1(ξ)x1+Δσ2(ξ)x2+E5
+Δσ3(ξ)x3][exp(iχL|xξ|)|xξ|exp(iχT|xξ|)|xξ|]}dSξ,
j=1,2,3E6

where the displacement continuity conditions across the inclusion are used, |xξ|is the distance between the field point x=(x1,x2,x3)and integration pointξ=(ξ1,ξ2,0), and Δσj(j=1,2,3)are the jumps of the interfacial stresses across the inclusion, which are defined by

Δσj(x)=σj3(x)σj3+(x),j=1,2,3,xS,σj3±(x)=limx3±0σ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Δσj:

uj0=1ω2MSΔσj(ξ)dSξ,
j=1,2,3E8
Ωj=(1)j+1ω2Mij2Sξ3jΔσ3(ξ)dSξ,
j=1,2E9
Ω3=1ω2Mi32S[ξ1Δσ2(ξ)ξ2Δσ1(ξ)]dSξE10

where ijis the radius of inertia of the inclusion with respect to the xj-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

SΔσ3(ξ)R3(x,ξ)dSξ=4πGχT2u3in(x),
xSE11
S[Δσj(ξ)Rj(x,ξ)+Δσ3j(ξ)Rj(3j)(x,ξ)]dSξ=4πGχT2ujin(x),
j=1,2,
xSE12

In Eq. (7), the kernelsRj, R12and R21have the form

Rj(x,ξ)=L1(|xξ|)(xjξj)2|xξ|2L2(|xξ|)4πρM(1+ξ3jx3ji32),
j=1,2E13
R3(x,ξ)=L1(|xξ|)4πρM(1+ξ1x1i22+ξ2x2i12)E14
Rkj(x,ξ)=(x1ξ1)(x2ξ2)|xξ|2L2(|xξ|)+4πρMξkxji32,k,j=1,2,
kjE15
Lj(d)=lj1(d)exp(iχLd)d3lj2(d)exp(iχTd)d3,
j=1,2E16
l11(d)=1iχLd,
l12(d)=1iχTdχT2d2,E17
l21(d)=33iχLdχL2d2,
l22(d)=33iχTdχT2d2.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Δσ3. After the solution of this equation the displacement u30and the rotations Ω1and Ω2can 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 Δσ1andΔσ2. After these quantities have been computed by solving these equations, the kinematical parameters u10,u20and Ω3can 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

ASΔσ3(ξ)|xξ|dSξ+SΔσ3(ξ)[1χT2R3(x,ξ)A|xξ|]dSξ=4πGu3in(x),
xSE19
SΔσj(ξ)|xξ|[A+B(xjξj)2|xξ|2]dSξ+BSΔσ3j(ξ)(x1ξ1)(x2ξ2)|xξ|3dSξ+E20
+S{Δσj(ξ)[1χT2Rj(x,ξ)A|xξ|B(xjξj)2|xξ|3]+Δσ3j(ξ)[1χT2Rj(3j)(x,ξ)E21
B(x1ξ1)(x2ξ2)|xξ|3]}dSξ=4πGujin(x),
j=1,2,
xSE22

where

A=(34ν)/[4(1ν)],
B=1/[4(1ν)]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ξx. Therefore, in the numerical evaluation of these integrals it is sufficient to perform the integration over Sx0by excluding a small region (the neighborhood of the x-point) around xfrom 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

Δσj(x)=fj(x)a2x12x22,j=1,2,3,xS,E24

where fj(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 functionsfj(x). These BIEs have a weak singularity 1/|xξ|at the source point ξ=xand 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 whenxS:

Sf(ξ)a2ξ12ξ22|xξ|dSξ=π2f(x)+Sx0f(ξ)f(x)a2ξ12ξ22|xξ|dSξ,E25
Sf(ξ)a2ξ12ξ22(xjξj)2|xξ|3dSξ=π22f(x)+Sx0[f(ξ)f(x)]a2ξ12ξ22(xjξj)2|xξ|3dSξ,E26
Sf(ξ)a2ξ12ξ22(x1ξ)(x2ξ2)|xξ|3dSξ=Sx0[f(ξ)f(x)]a2ξ12ξ22(x1ξ1)(x2ξ2)|xξ|3dSξ.E27

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

SdSξa2ξ12ξ22|xξ|=π2,E28
S(x1ξ1)k(x2ξ2)ja2ξ12ξ22|xξ|3dSξ={π2/2,whenk=2,j=0ork=0,j=2,0,whenk=1,j=1.E29

Next we perform the following transformation of the variables:

{x1=asiny1cosy2,x2=asiny1siny2,
{ξ1=asinη1cosη2,ξ2=asinη1sinη2,E30

where y=(y1,y2)and η=(η1,η2)are new variables in the rectangular domainS˜:{0y1,η1π/2;0y2,η22π}. 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η1=π/2.

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

Af˜3(y)[π2S˜y0sinη1R(y,η)dSη]+1ϑ2S˜y0f˜3(η)R˜3(y,η)sinη1dSη=4πGu˜3in(y)
y(y1,y2)S˜E31
,
f˜j(y){A[π2S˜y0sinη1R(y,η)dSη]+B[π22S˜y0Φj(y,η)sinη1dSη]}E32
Bf˜3j(y)S˜y0Ψ(y,η)sinη1dSη+1ϑ2S˜y0[f˜j(η)R˜j(y,η)+E33
+f˜3j(η)R˜j(3j)(y,η)]sinη1dSη=4πGu˜jin(y),
j=1,2,
y(y1,y2)S˜E34
,

where

f˜j(y)=fj(x),
u˜jin(y)=ujin(x),
R˜j(y,η)=a3Rj(x,ξ),     j=1,2,3,E35
R˜kj(y,η)=a3Rkj(x,ξ)
k,j=1,2,
kjE36
.

In Eq. (15), xand ξare defined by Eq. (13), ϑ=χTais the normalized wave number of the T-wave, S˜y0is the mapping of the domain Sx0due to the transformation (13) (in the domain S˜y0the points yand ηdo not coincide), and

R(y,η)=[sin2y1+sin2η12siny1sinη1cos(y2η2)]1/2E37
Φ1(y,η)=(siny1cosy2sinη1cosη2)2[R(y,η)]3E38
Φ2(y,η)=(siny1siny2sinη1sinη2)2[R(y,η)]3E39
Ψ(y,η)=(siny1cosy2sinη1cosη2)(siny1siny2sinη1sinη2)[R(y,η)]3E40

For the discretization of the domainS˜, 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 f˜jis obtained. After solving the system of linear algebraic equations numerically, the stress jumps Δσjacross 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Δσj. For this purpose we use the asymptotic relations for an observation point far away from the inclusion, namely |xξ||x|(xξ)/|x|and|xξ|1|x|1, when|x|. 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

x1=Rsinθcosφ,
x2=Rsinθsinφ,
x3=Rcosθ
0θ,φ2πE41

the asymptotic expressions for the scattered radial uRand tangential uθ,uφdisplacements in the far-field are obtained in the form

uR(R,θ,φ)=exp(iχLR)4πRFL(θ,φ,γ0),
RE42
uθ(R,θ,φ)=exp(iχTR)4πRFTV(θ,φ,γ0),
RE43
uφ(R,θ,φ)=exp(iχTR)4πRFTH(θ,φ,γ0),
RE44

Here, FL, FTV, and FTHare the longitudinal, vertically polarized transverse, and horizontally polarized transverse wave scattering amplitudes, respectively, which are related to the inclusion of normalized massγ0=M/(ρa3). They are given by

FL(θ,φ,γ0)=12ν2(1-ν)Gj=13pjSΔσj(ξ)exp[iχL(pξ)]dSξE45
FTV(θ,φ,γ0)=1Gj=13rjSΔσj(ξ)exp[iχT(pξ)]dSξE46
FTH(θ,φ,γ0)=1Gj=13τjSΔσj(ξ)exp[iχT(pξ)]dSξE47

where pj,rj,τj(j=1,2,3)are the coordinates of the spherical unit vectorsp=(sinθcosφ,sinθsinφ,cosθ), r=(cosθcosφ,cosθsinφ,-sinθ)andτ=(-sinφ,cosφ,0).

The forward scattering amplitudes are defined as the values of FZ(θ,φ,γ0)(Z=L,TV,TH)in the direction of the wave incidence, i.e.,FZ(θ0,0,γ0).

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 Otqzwith the center in the inclusion contour point, so that the value z=0corresponds to the inclusion plane, the axes Otand Oqlie 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=0can be approximated as [21]:

uz(r,ψ,ϕ)=2r2Gsinψ2(1sin2ψ2)KI(ϕ)+2r2Gcosψ2(24ν+cos2ψ2)KII(ϕ)+O(r32)E48
ut(r,ψ,ϕ)=2r2Gcosψ2(44νcos2ψ2)KI(ϕ)2r2Gsinψ2(1sin2ψ2)KII(ϕ)+O(r32)E49
uq(r,ψ,ϕ)=2r2Gcosψ2KIII(ϕ)+O(r32)E50
σzz(r,ψ,ϕ)=12rcosψ2(12νsinψ2sin3ψ2)KI(ϕ)+E51
+12rsinψ2(22νcosψ2cos3ψ2)KII(ϕ)+O(1)E52
σtt(r,ψ,ϕ)=12rcosψ2(32νsinψ2sin3ψ2)KI(ϕ)+12rsinψ2(2ν+cosψ2cos3ψ2)KII(ϕ)+O(1)E53
σqq(r,ψ,ϕ)=2ν2rcosψ2KI(ϕ)+2ν2rsinψ2KII(ϕ)+O(1)E54
σzq(r,ψ,ϕ)=12rsinψ2KIII(ϕ)+O(1)E55
σtq(r,ψ,ϕ)=12rcosψ2KIII(ϕ)+O(1)E56
σzt(r,ψ,ϕ)=12rsinψ2(22ν+cosψ2cos3ψ2)KI(ϕ)+E57
+12rcosψ2(12ν+sinψ2sin3ψ2)KII(ϕ)+O(1)E58

Here rand ψare the polar coordinates of the point P, ϕis the angular coordinate of the inclusion contour point (see Figure 1), KI,KIIand KIIIare the mode-I, II, and III dynamic stress intensity factors in the inclusion vicinity.

By using the Eq. (20) the KI,II,III-factors can be defined directly from the stress jumps Δσjor the solutions of BIEs (7) by the following relations:

KI(ϕ,γ0)=14(1ν)a{a2x12x22[Δσ1(x)cosϕ+Δσ2(x)sinϕ]}|x1=acosϕ;x2=asinϕ,E59
KII(ϕ,γ0)=14(1ν)a[a2x12x22Δσ3(x)]|x1=acosϕ;x2=asinϕ,E60
KIII(ϕ,γ0)=14a{a2x12x22[Δσ1(x)sinϕΔσ2(x)cosϕ]}|x1=acosϕ;x2=asinϕE61

where the dependence of KI,II,III-factors on the inclusion mass also is fixed by the variableγ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 x1x2-plane. The radius aof 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(ω)as

K(ω)=ωc(ω)+iα(ω)E62

where c(ω)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

u(x,ω)=Uexp[iK(nx)]=Uexp[α(ω)(nx)]exp[iω(nx)/c*(ω)]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 KZ(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

KZ2=χZ2+εa3FZ*,Z=L,TE64

In Eq. (24), εis the density parameter of the inclusions, εa3corresponds to the number density of inclusions of the same radius, i.e. the number of inclusions per unit volume, FZ*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

FZ*(θ0)=τβFZ(θ0,0,γ)g(γ)dγ,Z=L,TV,TH,E65

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

FZ*=120πτβFZ(θ0,0,γ)g(γ)sinθ0dγdθ0,Z=L,TV,TH.E66

Here θ0is 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 gof inclusion mass, and FZ(θ0,0,γ)are the forward scattering amplitudes given by Eq. (19). The density distribution function of inclusion mass should satisfy the normalization condition

τβg(γ)dγ=1E67

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

Distribution typeg(γ)
Alignedδ(γγ0)
Normal3γ2/(β3τ3)
Uniform1/(βτ)

Table 1.

Various types of density distribution functions of inclusion mass (δ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 KZhave been determined via Eq. (24), the effective wave velocities cZand the attenuation coefficients αZof the plane L- and T-waves can be obtained by considering the definition (22). This results in

cZ(ω)=ωRe[KZ(ω)],
αZ(ω)=Im[KZ(ω)]
Z=L,TV,THE68

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 S˜is divided into 264 rectangular elements of length Δy1=π/22and Δy2=π/12in the y1- and y2-directions, the domain S˜y0is chosen asS˜y0=S˜\ΔS˜y, where ΔS˜yis 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 asi1=i2=a/2,i3=a/2. In the numerical examples, the radius aof the inclusions is assumed to be equal, while their masses are varied by the distribution laws defined in Table 1 from the minimal value τ=5to the maximal valueβ=20. The inclusion density parameter is fixed as ε=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 asc¯Z=cZ/cZ, α¯Z=2aαZ/(πε),where the subscript Z=L,TVstands 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 (θ0=0o)of L- and T-waves and grazing incidence (θ0=90o)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 α¯Land the normalized effective wave velocity c¯Lversus the dimensionless wave number ϑare presented in Figure 3. Three separate parametrical studies are performed to show the effects of different inclusion mass distributions: aligned withγ0=β=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 range0ϑ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 α¯Lincreases rapidly with increasingϑ, after reaching a maximum it then decreases and approaches its high frequency limit (see Figure 3(a)). The peak value of α¯Lincreases 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 α¯Lin 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 α¯Lcan be obtained by using the Kirchhoff approximation for short waves [19]. At low frequencies, the normalized effective wave velocity c¯Lin 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 c¯Lin comparison with the aligned (uniform) situation. An opposite tendency is observed in the range of higher frequencies. In addition, the normalized effective wave velocity c¯Lin the composite can be bigger than that in the homogeneous matrix material, for instance for the considered inclusions of aligned mass andϑ>1. The high-frequency limit c¯L1at ϑ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 α¯Lshows a similar dependence on the dimensionless wave number ϑand the inclusion mass distribution. For the same density function gof inclusion mass, the peaks of α¯Lare larger than that for normal incidence as depicted in Figure 3(a), but the high-frequency limit of α¯Lin the case with θ0=90ois smaller than that in the case withθ0=0o. The normalized effective wave velocity c¯Lin 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 c¯Lfor a normal incidence of plane L-wave, the maximum effective wave velocity c¯Lis increased while its minimum value is decreased. Also in contrast to Figure 3(b), the normalized effective wave velocity c¯Lfor all considered functions gis 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 α¯Tand the normalized effective wave velocity c¯Tare presented in Figure 5 versus the dimensionless wave numberϑ. The global behavior of the normalized attenuation coefficient α¯Tand the normalized effective wave velocity c¯Tis very similar to that for grazing incidence (i.e.θ0=90o) of plane L-wave as presented in Figure 4. In comparison to the peak values of α¯Lfor normal and grazing incidence of a plane L-wave, the peak values of α¯Tfor 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 c¯Tare reduced compared to c¯Lfor normal and grazing incidence of a plane L-wave for the same inclusion mass distribution, while the maximum values of c¯Tlie between the 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 α¯Land the normalized effective wave velocity c¯Lare presented in Figure 6 versus the dimensionless wave numberϑ. A comparison of Figure 6 with Figures 3 and 4 for parallel inclusions shows a similar dependence of the α¯Land c¯Lon the dimensionless wave number ϑand the inclusion mass distribution. As expected, the peak values of α¯Land c¯Lfor randomly oriented penny-shaped inclusions are larger than that for parallel inclusions with θ0=0obut smaller than that for parallel inclusions withθ0=90o.

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 α¯Tand the normalized effective wave velocity c¯Tfor an incident plane TV-wave. In contrast to parallel inclusions (see Figure 5), now the variations of α¯Tand c¯Twith the dimensionless wave number ϑare rather complicated at low frequencies. For instance, the normalized attenuation coefficient α¯Tfor 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 α¯Tfor randomly oriented inclusions is smaller than that for parallel inclusions withθ0=0o, while the normalized effective wave velocity c¯Tcould be larger or smaller than that for parallel inclusions depending on the dimensionless wave number ϑand the density function gof 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 KI,KIIand KIIIin the vicinity of separate inclusion. We suppose the impinge on the inclusion of plane longitudinal L-wave with constant amplitudeU0, 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*=U0G/[4(1ν)a]is chosen as the normalizing factor for the amplitudes of dynamic stress intensity factors, so thatK¯I,II,III=|KI,II,III|/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γ0.

At normal L-wave incidence on the inclusion (antisymmetric problem) KI=KIII=0,and the mode-II dynamic stress intensity factor KIIdoes not vary along the inclusion contour (see Figure 8). It follows from the Figure 8 that in the initial range of wave numbers K¯II-factor rapidly increases from a zero value, what is more pronounced for the inclusions of large mass. A further increase in ϑin the case of an inclusion with the mass characteristic γ0=20leads to a local maximum ofK¯II. For higher wave numbers, a regularity is the approaching of K¯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 K¯IIand ϑ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) KII=0,and the mode-I and mode-III dynamic stress intensity factor KIand KIIIdepend on the angular coordinate ϕof the inclusion contour point. Hence, in Figure 9 K¯I-factor at the most representative points of the inclusion contour is showed, where KIII=0due 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 K¯I-factor at the point, where the wave runs on the inclusion. It should be mentioned thatK¯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 K¯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.

© 2012 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Viktor Mykhas'kiv (August 22nd 2012). Numerical Simulation of Wave Propagation in 3D Elastic Composites with Rigid Disk-Shaped Inclusions of Variable Mass, Composites and Their Applications, Ning Hu, IntechOpen, DOI: 10.5772/48113. Available from:

chapter statistics

1555total chapter downloads

1Crossref citations

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

Structural Health Monitoring for Composite Materials

By Jian Cai, Lei Qiu, Shenfang Yuan, Lihua Shi, PeiPei Liu and Dong Liang

Related Book

First chapter

Graphene-Boron Nitride Composite: A Material with Advanced Functionalities

By Sumanta Bhandary and Biplab Sanyal

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us