A direct method for wave propagation in elastic fiber composites

Abstract

A direct method is proposed to study wave dispersion in stiff unidirectional fiber-reinforced elastic random composites. The present model is featured by a modified form of classical elastodynamic equations supplied with a simple differential relation between the displacement field of the composite and the displacement field of the mass center of representative cell. Explicit formulas are derived for the attenuation coefficient and effective phase velocity of in-plane P-waves, SV-waves, and anti-plane SH-waves, respectively. The efficiency and accuracy of the model are demonstrated by comparing predicted results (on the frequency for the maximum attenuation and effective phase velocity at frequencies below the bandgap) with a large group of known data reported in past decades. The proposed model offers a novel general method for dynamics of stiff unidirectional fiber-reinforced elastic random composites.

1. Introduction

Unidirectional fiber-reinforced composites find important application in numerous industrial areas [1–4], and the study of their mechanical behavior is a topic of long history [5–9]. In the past decades, dynamic behavior, particularly scattering-caused wave attenuation, of unidirectional fiber-reinforced elastic composites has been the topic of numerous researches [10–24]. Existing works studied wave scattering of unidirectional fiber-reinforced composites and developed numerical models for wave dispersion based on various effective medium theories. Extensive numerical calculation is usually requested for wave dispersion analysis, and different numerical models often give quite different results for attenuation and phase velocity, typically with 20%–30% or even larger relative errors between different models. To mention a few, for example, Wang and Gan [19] studied phase velocity and attenuation of anti-plane SH-waves of fiber-reinforced composites with a self-consistent scheme; Sato et al. [12] developed a plane-strain boundary-element analysis of P/SV-waves for stiff fiber-reinforced concretes; Kim [21] conducted a comparative study on anti-plane SH-waves in fiber-reinforced composites; and more recently, Sumiya et al. [16] developed a computational analysis of P/SV-waves in fiber-reinforced composites, Zhang et al. [24] conducted a numerical simulation on effective phase velocity and attenuation of SH-waves in fiber composites, and Yu et al. [13] studied P-waves in aluminum rod-reinforced polymers with a plane-strain two-dimensional numerical model.

In spite of extensive research, however, most existing models usually request rather complicated mathematical formulation and extensive numerical calculations, and different models often give quite different, or even qualitatively different, predictions particularly on wave attenuation, see, for example, comparative studies [13,15,18,19,21]. For example, it remains a relevant question how to predict the frequency for the maximum wave attenuation. Therefore, it is of theoretical and practical interest to develop new analytical models for wave dispersion in unidirectional fiber-reinforced elastic random composites.

The present work aims to develop a new general model to analyze elastic wave propagation in stiff unidirectional fiber-reinforced random composites. Unlike existing models, the present model is based on the concept that the derivation of the displacement field of embedded stiff fibers from the displacement field of the composite dominates wave dispersion of stiff unidirectional fiber-reinforced elastic composites. As shown in section 2, the novelty of the present model is featured by the idea to modify the classical elastodynamic equations by substituting the inertia term by the acceleration field of the mass center of the representative unit cell, which is supplied with a simple differential relation between the mass center’s displacement field and the composite’s displacement field. Explicit dispersion relations are derived in section 3 for P-waves, SV-waves, and SH-waves. The efficiency and accuracy of the proposed model are demonstrated in section 4 by comparing its predicted results to known data reported in past decades. Finally, the main conclusions are summarized in section 5.

2. Basic equations of the present model

Since transversely isotropic (such as carbon) fibers are widely used in industries, let us consider a transversely isotropic elastic composite comprised of an isotropic elastic matrix (of the Young’s and shear moduli E_m and G_m) reinforced by transversely isotropic cylindrical fibers (defined by five independent elastic constants: E₁ ^f , E₂ ^f , G₁₂ ^f , G₂₃ ^f , ν₁₂ ^f ) aligned along the x₁-axis and randomly distributed with a constant volume fraction δ, as shown in Figure 1. Here, it is assumed that the identical fibers are of circular cross section of radius a, and the elastic moduli of fibers are much (say, at least 4–5 times) higher than that of the compliant matrix so that the elastic solutions given in Appendixes 1 and 2 for a rigid circular fiber can be used, and the characterized wavelengths are much larger than the diameter of fibers and the average distance between the centers of adjacent fibers.

Figure 1.

Unidirectional fiber-reinforced elastic random composites.

Based on a widely accepted idea that the derivation of the displacement field of embedded stiff fibers from the displacement field of the composite is largely responsible for novel dynamic behavior of stiff unidirectional fiber-reinforced elastic metacomposites, let us consider two x₁-independent (plane or anti-plane) displacement fields: the displacement field u _f (x₂, x₃, t) = (w_f, u_f, v_f) of the embedded fibers (defined by the displacement of the central axis of a fiber), and the displacement field u(x₂, x₃, t) = (w, u, v) of the composite (which is considered as the displacement field of the geometrical center of the representative unit cell). Thus, strain ε of the composite is calculated based on the composite’s displacement u(x₂, x₃, t), and the three-dimensional (3D) transversely isotropic strain ε/stress σ relations of the composite can be obtained from the effective composite model given in Appendix 1. The mass density ρ of the composite (per unit volume of the composite) is defined by,

$$\rho =(1-\delta ){\rho }_m+{\rho }_f\delta ,$$

(1)

where δ is the volume fraction of the embedded fibers, and ρ_m and ρ_f are the mass densities of the matrix phase and the embedded fibers, respectively. For any deformable material element such as a representative unit cell, the resultant external force acting on the element equates to the mass of the element multiplied by the acceleration of its mass center. Thus, in the absence of body forces, the following modified form of classical (plane or anti-plane) elastodynamic equations is proposed where the inertia term on LHS is substituted by the acceleration field of the mass center of the representative unit cell:

$$\begin{array}{c}\rho \frac{{\partial }^2{w}_{\mathrm{mass}}}{\partial t^2}=\frac{\partial {\mathit{\sigma }}_{12}}{\partial x_2}+\frac{\partial {\sigma }_{13}}{\partial x_3},\hfill \\ \rho \frac{{\partial }^2{u}_{\mathrm{mass}}}{\partial t^2}=\frac{\partial {\mathit{\sigma }}_{22}}{\partial x_2}+\frac{\partial {\mathit{\sigma }}_{23}}{\partial x_3},\hfill \\ \rho \frac{{\partial }^2{v}_{\mathrm{mass}}}{\partial t^2}=\frac{\partial {\mathit{\sigma }}_{23}}{\partial x_2}+\frac{\partial {\mathit{\sigma }}_{33}}{\partial x_3},\hfill \end{array}$$

(2)

where σ_ij (i, j = 1, 2, 3) are the (symmetric) stress components which can be given in terms of the displacement field u of the composite, the RHS of equation (2) represents the resultant external force acting on the representative unit cell, and u_mass = (w_mass, u_mass, v_mass) is the displacement field of the mass center of representative unit cell given by,

$$\rho {\mathbf{u}}_{\mathrm{mass}}={\rho }_m(1-\delta )\mathbf{u}+\delta {\rho }_f{\mathbf{u}}_f.$$

(3)

For stiff unidirectional fiber-reinforced composites characterized by stiff (such as steel/glass/carbon) fibers and soft matrix phase (such as polymer or rubber), the deviation (u _f –u) of the displacement field of embedded fibers from the displacement field of the composite is responsible for novel dynamic behavior of the composite, which implies that u_mass ≠ u on the two sides of equation (2).

Based on the spring model for the x₁-independent displacement field of a single rigid fiber surrounded by an infinite effective medium (see Appendix 2), we have a simple differential relationship between u_mass(w_mass, u_mass, v_mass) and u(w, u, v) (of two spatial coordinates (x₂, x₃)) as follows:

$${\mathbf{u}}_{\mathrm{mass}}=\mathbf{u}-\frac{\delta {\rho }_f}{\rho \beta }[\rho \frac{{\partial }^2{\mathbf{u}}_{\mathrm{mass}}}{\partial t^2}-{\rho }_m(1-\delta )\frac{{\partial }^2\mathbf{u}}{\partial t^2}],$$

(4)

where the spring constant β (per unit volume of the composite) is given in Appendix 2 for the decoupled plane strain–stress case (between (u_mass, v_mass) and (u, v)) and anti-plane case (between w_mass and w), respectively, also see equations (9) and (23). Thus, with the present direct method, the x₁-independent (plane or anti-plane) dynamics of stiff unidirectional fiber-reinforced composites is governed by the modified form of classical equation (2) coupled with the differential relation equation (4) between the mass-center displacement u_mass of the representative unit cell and the displacement u of the composite.

3. Dispersion relations for plane waves

As compared to existing methods [10–24], the present model based on governing equations (2) and (4) enjoys conceptual and mathematical simplicity and the generality although it is restricted to stiff fibers whose elastic moduli are much higher than those of the compliant matrix. In this section, explicit dispersion relations will be derived from equations (2) and (4), for various plane waves in the stiff unidirectional fiber-reinforced elastic random composites.

3.1. Longitudinal P-waves

First, let us consider the plane P-waves under the conditions of plane strain with u(x₂, x₃, t) and v(x₂, x₃, t) and w ≡ 0 (see Figure 1). In terms of the two effective elastic coefficients (c₂₂, c₂₃), the stress–strain relations in plane-strain conditions are given by (see Appendix 1),

$$\left(\begin{array}{c}{\sigma }_{22}\hfill \\ \multicolumn{1}{c}{{\sigma }_{33}}\end{array}\right)=\left(\begin{array}{cc}{c}_{22}\hfill & c_{23}\hfill \\ \multicolumn{1}{c}{c_{23}}& c_{22}\hfill \end{array}\right)\left(\begin{array}{c}{\epsilon }_{22}\hfill \\ \multicolumn{1}{c}{{\epsilon }_{33}}\end{array}\right),{\sigma }_{23}=2G_{23}{\epsilon }_{23},G_{23}=\frac{1}{2}(c_{22}-c_{23})=\frac{E_2}{2(1+{\nu }_{23})},$$

(5)

where,

$$c_{22}=\frac{E_2(E_1-E_2{\nu }_{12}^2)}{(1+{\nu }_{23})[E_1(1-{\nu }_{23})-2E_2{\nu }_{12}^2]},c_{23}=\frac{E_2(E_1{\nu }_{23}+E_2{\nu }_{12}^2)}{(1+{\nu }_{23})[E_1(1-{\nu }_{23})-2E_2{\nu }_{12}^2]}.$$

(6)

For plane P-waves propagating in x₂-direction, we have u(x₂, t) with v = 0, thus ${\sigma }_{23}=0;{\sigma }_{22}=c_{22}{u}_{,2},{\sigma }_{33}=c_{23}{u}_{,2}$ , and $\rho ({\partial }^2{u}_{\mathrm{mass}}/\partial t^2)=c_{22}{u}_{,22}$ . Substituting it into the above relationship equation (4), we have,

$$c_{22}{u}_{,22}=\rho u_{,\mathrm{tt}}-\frac{\delta {\rho }_f}{\beta }[c_{22}{u}_{,22\mathrm{tt}}-{\rho }_m(1-\delta )u_{,\mathrm{tttt}}].$$

(7)

Thus, for harmonic P-waves u(x₂, t) =u₀e^-i(kx₂^-ωt), where u₀ is a constant, ω is circular frequency, and k is the wave number, then,

$$k^2=\rho {\omega }^2\frac{(1-\frac{\delta {\rho }_f}{\beta \rho }{\rho }_m(1-\delta ){\omega }^2)}{c_{22}(1-\frac{\delta {\rho }_f}{\beta }{\omega }^2)}.$$

(8)

Within the following range of frequencies (called the “bandgap” for P-waves and SV-waves under plane-strain, where κ = 3-4ν₂₃, see Appendix 2),

$${\omega }_0^2\le {\omega }^2\le {\omega }_0^2\left(\frac{\rho }{{\rho }_m(1-\delta )}\right),{\omega }_0^2=\frac{\beta }{\delta {\rho }_f},\beta =\frac{\delta {\beta }_0}{\pi a^2}=\frac{4G_{23}(1+\kappa )\delta }{a^2[\kappa \ln \left(\frac{\pi }{\delta }\right)-(1-\frac{\delta }{\pi })]}.$$

(9)

We have,

$$k=-i\omega \sqrt{\frac{\rho (1-\frac{\delta {\rho }_f}{\beta \rho }{\rho }_m(1-\delta ){\omega }^2)}{c_{22}(\frac{\delta {\rho }_f}{\beta }{\omega }^2-1)}},c_P=\infty ,\alpha =\omega \sqrt{\frac{\rho (1-\frac{\delta {\rho }_f}{\beta \rho }{\rho }_m(1-\delta ){\omega }^2)}{c_{22}(\frac{\delta {\rho }_f}{\beta }{\omega }^2-1)}}.$$

(10)

For a visual illustration of general wave behavior predicted by the present model without any restriction to specific material and geometrical parameters, phase velocity and attenuation predicted by the present model are shown graphically in Figures 2(b) and 3(b) without specific parameters, with a comparison to actual wave behavior of unidirectional fiber-reinforced elastic composites shown in Figures 2(a) and 3(a). Thus, it is seen from Figures 2(b) and 3(b) that the phase velocity given by the present model is infinite within the bandgap equation (9) while the attenuation is infinite at the lower cut-off frequency ω₀ and decreases to zero at the upper cut-off frequency of the bandgap, while actual phase velocity and attenuation of unidirectional fiber-reinforced composites observed in experiments or predicted by more complicated numerical models (see, e.g., several studies [10–13,17–21]) are illustrated in Figures 2(a) and 3(a), respectively. Clearly, the actual attenuation at the local resonance frequency ω = ω₀ cannot be infinite, and the actual phase velocity within the bandgap cannot be infinite. It is not surprising that, just similar to all linear undamped theories on resonance, the present (linear and undamped) model cannot give quantitative details of resonant behavior at and around the local resonance frequency ω = ω₀. In spite of this, however, the results given by the present model (shown in Figures 2(b) and 3(b)) are qualitatively consistent with actual wave behavior (shown in Figures 2(a) and 3(a)), because both show that the phase velocity rises up to the maximum around the bandgap while the attenuation rises up abruptly to the highest peak around the local resonance frequency ω = ω₀.

Figure 2.

A graphical illustration of phase velocity given by the present model with comparison to actual phase velocity (e.g., see Sato and his colleagues [11,12] for P/SV-waves, and several studies [17–19,21] for SH-waves): (a) Actual phase velocity and (b) the present model.

Figure 3.

A graphical illustration of attenuation given by the present model with comparison to actual attenuation (e.g., see several studies [10–13] for P/SV-waves, and several studies [17–21] for SH-waves): (a) Actual attenuation and (b) the present model.

On the other hand, outside the bandgap equation (9), the present model gives that the attenuation α = 0 and the longitudinal phase velocity c _P for P-waves decreases with increasing ω in two separate ranges, below or above the bandgap respectively, and given by,

$$k=\omega \sqrt{\frac{\rho (1-\frac{\delta {\rho }_f}{\beta \rho }{\rho }_m(1-\delta ){\omega }^2)}{c_{22}(1-\frac{\delta {\rho }_f}{\beta }{\omega }^2)}},\alpha =0,c_P=\frac{\omega }{k}=\sqrt{\frac{c_{22}(1-\frac{\delta {\rho }_f}{\beta }{\omega }^2)}{\rho (1-\frac{\delta {\rho }_f}{\beta \rho }{\rho }_m(1-\delta ){\omega }^2)}}.$$

(11)

In particular, for low frequencies well below the bandgap, the phase velocity is given approximately by,

$$c_P\approx \sqrt{\frac{c_{22}}{\rho }}(1-\frac{\delta {\rho }_f}{2\rho }\left(\frac{\delta {\rho }_f}{\beta }{\omega }^2\right)),\left(\frac{\delta {\rho }_f}{\beta }{\omega }^2\right)<<1.$$

(12)

For high frequencies well above the bandgap, we have,

$$c_P\approx \sqrt{\frac{c_{22}}{{\rho }_m(1-\delta )}}>\sqrt{\frac{c_{22}}{\rho }}.$$

(13)

It is seen from equation (12) that the phase velocity c _P for P-waves predicted by the present model decreases very slowly with increasing frequency well below the bandgap and drops drastically beyond the bandgap, as shown in Figure 2(b), qualitatively consistent with actual phase velocity observed in experiments or predicted by more complicated numerical models shown in Figure 2(a) which show that the phase velocity decreases slowly with increasing frequency well below the bandgap and drops from the peak value beyond the bandgap, while the actual attenuation remains vanishingly small (more than 2–3 orders of magnitude lower than the maximum attenuation) outside the bandgap, as shown in Figure 3(a), qualitatively consistent with the attenuation shown in Figure 3(b) predicted by the present model. Therefore, the results given by the present model are qualitatively consistent with actual wave behavior, because both show that the phase velocity drops considerably when the bandgap is approached and rises up within the bandgap, while the attenuation remains vanishingly low for frequencies well below or above the bandgap and rises up abruptly to the highest peak around the local resonance frequency ω = ω₀. A more detailed quantitative comparison of the present model with known data on P-waves will be made in section 4.

3.2. Transverse SV-waves

Similarly, for transverse SV-waves propagating in x₂-direction, we have v(x₂, t) and u = 0, thus ${\sigma }_{23}=G_{23}{v}_{,2};{\sigma }_{22}={\sigma }_{33}=0$ and $\rho ({\partial }^2{v}_{\mathrm{mass}}/\partial t^2)=G_{23}{v}_{,22}$ . For harmonic SV-waves, v(x₂, t) =v₀e^-i(kx₂^-ωt) with a constant v₀, substituting it into the relationship equation (4) gives,

$$G_{23}{v}_{,22}=\rho v_{,\mathrm{tt}}-\frac{\delta {\rho }_f}{\beta }[G_{23}{v}_{,22\mathrm{tt}}-{\rho }_m(1-\delta )v_{,\mathrm{tttt}}].$$

(14)

Thus, we have,

$$k^2=\rho {\omega }^2\frac{(1-\frac{\delta {\rho }_f}{\beta \rho }{\rho }_m(1-\delta ){\omega }^2)}{G_{23}(1-\frac{\delta {\rho }_f}{\beta }{\omega }^2)}.$$

(15)

It is noticed from equations (8) and (15) that, because the present model does not consider the effect of rotation of embedded fibers on local resonance, the predicted local resonance frequency for SV-waves (obtained by setting the denominator in equation (15) to be zero) is identical to that given by equation (9) for P-waves. Thus, for SV-waves within the bandgap defined by equation (9), we have,

$$c_{\mathrm{SV}}=\infty ,\alpha =\omega \sqrt{\frac{\rho (1-\frac{\delta {\rho }_f}{\beta \rho }{\rho }_m(1-\delta ){\omega }^2)}{G_{23}(\frac{\delta {\rho }_f}{\beta }{\omega }^2-1)}}.$$

(16)

On the other hand, outside the bandgap, we have α = 0, and the transverse phase velocity c_SV for SV-waves decreases with increasing ω in two separate ranges, below or above the bandgap, respectively,

$$\alpha =0,c_{\mathrm{SV}}=\frac{\omega }{k}=\sqrt{\frac{G_{23}(1-\frac{\delta {\rho }_f}{\beta }{\omega }^2)}{\rho (1-\frac{\delta {\rho }_f}{\beta \rho }{\rho }_m(1-\delta ){\omega }^2)}}.$$

(17)

In particular, for low frequencies well below the bandgap, the phase velocity c_SV is given approximately by,

$$c_{\mathrm{SV}}\approx \sqrt{\frac{G_{23}}{\rho }}(1-\frac{\delta {\rho }_f}{2\rho }\left(\frac{\delta {\rho }_f}{\beta }{\omega }^2\right)),\left(\frac{\delta {\rho }_f}{\beta }{\omega }^2\right)<<1.$$

(18)

For high frequencies well above the bandgap, we have,

$$c_{\mathrm{SV}}\approx \sqrt{\frac{G_{23}}{{\rho }_m(1-\delta )}}>\sqrt{\frac{G_{23}}{\rho }}.$$

(19)

The phase velocity c_SV and the attenuation of plane SV-waves predicted by the present model are qualitatively similar to those of plane P-waves, as shown in Figures 2(b) and 3(b). A more detailed quantitative comparison of the present model with known data on SV-waves will be made in section 4.

3.3. Anti-plane SH-waves

For anti-plane deformation w(x₂, x₃, t) with u = v = 0, the stress–strain Hookean law gives (see Appendix 1),

$${\sigma }_{12}=2G_{12}{\epsilon }_{12},{\sigma }_{13}=2G_{12}{\epsilon }_{13}.$$

(20)

For SH-wave propagating along x₂-direction, we have w(x₂, t) =w₀e^-i(kx₂^-ωt) with a constant w₀, thus ${\sigma }_{12}=G_{12}{w}_{,2};{\sigma }_{13}=0$ , and $\rho ({\partial }^2{w}_{\mathrm{mass}}/\partial t^2)=G_{12}{w}_{,22}$ . Substituting it into the relationship equation (9) gives,

$$G_{12}{w}_{,22}=\rho w_{,\mathrm{tt}}-\frac{\delta {\rho }_f}{\beta }[G_{12}{w}_{,22\mathrm{tt}}-{\rho }_m(1-\delta )w_{,\mathrm{tttt}}].$$

(21)

Thus, we have,

$$k^2=\rho {\omega }^2\frac{(1-\frac{\delta {\rho }_f}{\beta \rho }{\rho }_m(1-\delta ){\omega }^2)}{G_{12}(1-\frac{\delta {\rho }_f}{\beta }{\omega }^2)}.$$

(22)

Within the following range of frequencies (called the “bandgap” for SH-waves under anti-plane shearing, see Appendix 2) defined by,

$${\omega }_w^2\le {\omega }^2\le {\omega }_w^2\left(\frac{\rho }{{\rho }_m(1-\delta )}\right),{\omega }_w^2=\frac{\beta }{\delta {\rho }_f},\beta =\frac{\delta {\beta }_w}{\pi a^2}=\frac{4\delta G_{12}}{a^2\ln (\pi /\delta )}.$$

(23)

We have,

$$c_{\mathrm{SH}}=\infty ,\alpha =\omega \sqrt{\frac{\rho (1-\frac{\delta {\rho }_f}{\beta \rho }{\rho }_m(1-\delta ){\omega }^2)}{G_{12}(\frac{\delta {\rho }_f}{\beta }{\omega }^2-1)}}.$$

(24)

On the other hand, outside the bandgap, we have α = 0, and the phase velocity c_SH for SH-waves decreases with increasing ω in two separate ranges, below or above the bandgap, respectively,

$$\alpha =0,c_{\mathrm{SH}}=\frac{\omega }{k}=\sqrt{\frac{G_{12}(1-\frac{\delta {\rho }_f}{\beta }{\omega }^2)}{\rho (1-\frac{\delta {\rho }_f}{\beta \rho }{\rho }_m(1-\delta ){\omega }^2)}}.$$

(25)

In particular, for low frequencies well below the bandgap, the phase velocity c_SH is given approximately by,

$$c_{\mathrm{SH}}\approx \sqrt{\frac{G_{12}}{\rho }}(1-\frac{\delta {\rho }_f}{2\rho }\left(\frac{\delta {\rho }_f}{\beta }{\omega }^2\right)),\left(\frac{\delta {\rho }_f}{\beta }{\omega }^2\right)<<1.$$

(26)

For high frequencies well above the bandgap, we have,

$$c_{\mathrm{SH}}\approx \sqrt{\frac{G_{12}}{{\rho }_m(1-\delta )}}>\sqrt{\frac{G_{12}}{\rho }}.$$

(27)

Thus, the phase velocity c_SH and the attenuation of plane SH-waves predicted by the present model are qualitatively similar to those of plane P-waves and SV-waves, as shown in Figures 2(b) and 3(b). A more detailed quantitative comparison of the present model with known data on SH-waves will be made in section 4.

4. Comparison with known data

To demonstrate the efficiency and reasonable accuracy of the present model, let us make a detailed quantitative comparison of the results predicted by the present model with some known data published in past decades. Since most experimental data and numerical simulations show that the variation of phase velocity with frequency is relatively minor (usually not more than 10%–20%), the specific frequency at which the maximum attenuation appears is of major interest. Therefore, the local resonance frequency ω₀ predicted by equation (9) for P/SV-waves and the local resonance frequency ω_w predicted by equation (23) for SH-waves will be compared to the known data on the frequency for the maximum attenuation.

4.1. The frequency for the maximum attenuation (α)_max

The local resonance frequency ω₀ predicted by equation (9) is compared to the frequency for the maximum attenuation in P/SV-waves reported by more accurate numerical models and experiments, as shown in Table 1. The local resonance frequency ω_w predicted by equation (23) is compared to the frequency for the maximum attenuation in SH-waves reported by more accurate numerical models, as shown in Table 2.

Table 1.

The local resonance frequency ω₀ predicted by the present model equation (9) with comparison to known simulation/experimental data on the frequency for the maximum attenuation coefficient (α)_max in P/SV-waves reported in several studies [10–13] (where all fibers are isotropic, with the effective Poisson ratio ν of the composite in [0.3, 0.4], c_Pm and c_Sm are the speeds of longitudinal/shear waves in the isotropic matrix).

Reference	Yu et al. [13]P-wave experiment, Al-alloy/resin	Sato et al. [12], P/SV-wave,steel/mortar	Sato and Shindo [11], P/SV-wave,Sic/aluminum	Varadan et al. [10], P-wave, boron/aluminum
Parameters cited in reference	ρf ≈ 1.7ρm, Gf ≈ 7Gm, δ = 0.2,cSm = 1580 m/s, G23 ≈ 1.6Gm, a = 0.8 mm	ρf ≈ 5.6ρm, Gff ≈ 7.1Gm, δ = 0.3,G23 ≈ 1.9Gm	ρf ≈ 1.2ρm, Gf ≈ 7Gm, δ = 0.3, G23 ≈ 1.9Gm	ρf ≈ 0.93ρm, Gf ≈ 6.5Gm, δ = 0.56,G23 ≈ 2.7Gm, (cPm/cSm) ≈ 1.94
ω 0 given by present model equation (9)	(aω0/cSm) ≈ 1.78 f = ω0/2π≈ 560 (kHz)	(aω0/cSm) ≈ 1.35	(aω0/cSm) ≈ 2.9	(aω0/cPm) ≈ 2.2
Reported frequency for (α)max	f ≈ 380 (kHz) Figure 8(b) in Yu et al. [13]	(aω/cSm) ≈ 1.3/1.4 Figures 6 and 7 in Sato et al. [12]	(aω/cSm) ≈ 2.5/2.2 Figures 8 and 9 in Sato and Shindo [11]	(aω/cPm) ≈ 2.3 Figure 6in Varadan et al. [10]

Table 2.

The local resonance frequency ω_w predicted by the present model equation (23) with comparison to known simulation data on the frequencies for the maximum attenuation coefficient (α)_max for SH-waves reported in several studies [17–21] (where all fibers are isotropic, and c_Sm is the speed of shear waves in the isotropic matrix).

Reference	Kim [21],Sic/Titanium	Wang and Gan [19]	Kim [18], steel/plexiglas	Beltzer and Brauber [17],stiff fibers	Kanaun et al. [20],glass/epoxy
Parameters citedin reference	ρf≈ 0.7ρm, Gf≈ 4.4Gm,δ = 0.33, G12≈ 1.8Gm,	ρf≈ 9.4ρm, Gf≈ 4.8Gm,δ = 0.20, G12≈ 1.55Gm	ρ f ≈ 6.8ρm, Gf≈ 51Gm, δ = 0.27, G12≈ 2.1Gm	ρf≈ 5ρm, Gf≈ 5Gm, δ = 0.3, G12≈ 1.8Gm	ρf≈ 2.1ρm, Gf≈ 18.8Gm, δ = 0.5, G12≈ 3.4Gm
ωw by equation (23)	(aωw/cSm) ≈ 2.15	(aω w /c Sm ) ≈ 0.49	(aω w /c Sm ) ≈ 0.71	(aω w /c Sm ) ≈ 0.78	(aω w /c Sm ) ≈ 1.88
Frequency for (α)max	(aω/cSm) ≈ 2.0–2.7 Figure 6 in Kim [21]	(aω/cSm) ≈ 0.67 Figure 5 in Wang andGan [19]	(aω/cSm) ≈ 1.0 Figure 4 in Kim [18]	(aω/cSm) ≈ 0.6 Figure 6(b) in Beltzer and Brauber [17]	(aω/cSm) ≈ 2.0 Figure 3 in Kanaun et al. [20]

In view of the fact that the existing numerical models often give quite scattering data (typically, with 20%–30% or even larger relative errors between different models) for the frequency of maximum attenuation, it is reasonable to conclude from Tables 1 and 2 that the local resonance frequency ω₀ or ω_w predicted by the present model equations (9) or (23), for P/SV-waves or SH-waves, respectively, is in reasonable agreement with the known data on the frequency for the maximum attenuation (α)_max, consistent with the concept that the local resonance of embedded fibers dominates wave attenuation of stiff unidirectional fiber-filled elastic random composites. In particular, the present model gives the same local resonance frequency equation (9) for P-waves and SV-waves, reasonably consistent with most known data, for instance, the simulation results of Sato and his colleagues [11,12] cited in Table 1 which show that the two frequencies for the maximum attenuation of P-waves and SV-waves are close to each other.

4.2. Phase velocity at low frequencies

In addition, some previous works have focused on lower frequencies well below the bandgap. In such a case, the present model predicts that the attenuation vanishes (α = 0), in agreement with known numerical simulations which showed that the attenuation at frequencies well below the bandgap is at least 2–3 orders of magnitude smaller than the maximum attenuation. Therefore, for frequencies well below the bandgap, let us focus on the effective phase velocity of the composites.

For example, the effective phase velocities at the long-wavelength limit (ω = 0) predicted by the present model, equation (12) for P-waves, equation (18) for SV-waves, and equation (26) for SH-waves, are compared with known data, as shown in Tables 3 and 4, respectively. Again, it is concluded from Tables 3 and 4 that the effective phase velocity predicted by the present model for frequencies well below the bandgap is generally in reasonable agreement with known data.

Table 3.

Effective phase velocity of P/SV-waves at the long-wavelength limit (ω = 0) predicted by the present model equations (12) and (18) with comparison to known simulation data (where the fibers in Verbis et al. [15] and Sumiya et al. [16] are isotropic while the graphite fibers in Yang and Mal [14] are transversely isotropic).

Reference	Sumiya et al. [16],Sic/Ti-alloy	Verbis et al. [15],Boron/epoxy	Yang and Mal [14]Graphite/epoxy
Related parameters	ρf≈0.6ρm, δ = 0.25, Gf≈177 GPa≈4Gm, ρ≈0.9ρm, G23≈1.6Gm	ρf≈2.1ρ m , δ = 0.54, Gf≈161 GPa≈89Gm,ρ≈1.6ρm, G23≈3.7Gm,a = 0.05 mm, cSm = 1190 m/s	ρf≈1.4ρm, δ = 0.65, (G23)f≈5 GPa≈2.5Gm, ρ≈1.26ρm, G23≈3.8 GPa≈1.9Gm
The values reportedin reference	(cP/c Sm )≈2.37, (cSV/cSm)≈1.2, Figure 9 in Sumiya et al. [16]	c P ≈3.8 km/s, Figure 5 in Verbis et al. [15] cSV≈1.8 km/s, Figure 6 in Verbis et al. [15]	(cP/cSm)≈2.36, (cSV/cSm)≈1.17, Figure 6 in Yang and Mal [14]
The values given bypresent model	c 22/G23≈3.6, (cP/cSm)≈2.5, (cSV/cSm)≈1.3	c 22/G23≈3.3, cP≈3.3 km/s, cSV≈1.82 km/s	c 22/G23≈4.0, (cP/cSm)≈2.4, (cSV/cSm)≈1.2

Table 4.

Effective phase velocity of SH-waves at the long-wavelength limit (ω = 0) predicted by the present model equation (26) with comparison to known simulation data (where the fibers in Verbis et al. [15] are isotropic while the graphite fibers in Yang and Mal [14] and Kanaun et al. [21] are transversely isotropic).

Reference	Kim [21],graphite/epoxy	Verbis et al. [15],boron/epoxy	Yang and Mal [14]graphite/epoxy
Related parameters	ρf≈1.04ρm, cSm = 1254 m/s, (G12) f ≈21 GPa = 10.6 Gm, δ = 0.52, G12≈2.86Gm	ρf≈2.1ρm, δ = 0.54, Gf≈161 GPa≈89Gm, ρ≈1.6ρm, G23≈3.7Gm, a = 0.05 mm, cSm = 1190 m/s	ρf≈1.4ρ m , δ = 0.65, (G12) f ≈24 GPa≈12Gm, ρ≈1.26ρm, G12≈3.8Gm
The value reported in reference	c SH≈1980 m/s Figure 4 in Kim [21]	c SH≈1.75 km/s Figure 4(a) in Verbis et al. [15]	(cSH/cSm)≈1.65 Figure 6 in Yang and Mal [14]
The value given by present model	c SH≈2080 m/s	c SH≈1.80 km/s	(cSH/cSm)≈1.70

In addition, for frequencies well below the bandgap, phase velocity predicted by the present model decreases very slowly with increasing frequency. For instance, it is seen from equations (12), (18), and (26) that when frequency increases from ω = 0 to ω = ω₀/2 (for P/SV-waves) or ω = ω_w/2 (for SH-waves), the ratio of the decrease in phase velocity to the phase velocity at ω  = 0 is given by δρ_f/(8ρ). For example, for all nine cases shown in Tables 1 and 2, it can be verified that phase velocity predicted by the present model decreases by 3%–9% when frequency increases from zero to a half of the local resonance frequency. This prediction is comparable to some known simulation results, see the examples in several studies [10,11,19] shown in Tables 1 and 2 where the phase velocity decreases by about 4% to 5% when frequency increases from zero to a half of the local resonance frequency.

5. Conclusion

A direct method is proposed to study wave propagation in stiff unidirectional fiber-reinforced elastic random composites. Compared to existing analytical or numerical models, the present model highlights the role of local resonance of embedded stiff fibers and enjoys conceptual and mathematical simplicity and the generality. The governing equations of the present model are given by a modified form of classical equations of elastodynamics (in which the inertia term is replaced by the acceleration field of the mass center of the representative unit cell) coupled with a simple differential relation between the displacement field of the mass center of representative unit cell and the displacement field of the composite. The efficiency and reasonable accuracy of the present model are demonstrated by good agreement between the predicted results and a large group of known data on the frequency for the maximum attenuation and the effective phase velocity at lower frequencies. The proposed model could offer a new general method to study various dynamic problems of stiff unidirectional fiber-reinforced elastic random composites, and an extension of the present method to hard sphere-reinforced elastic random composites will be studied in future work.