A direct method for wave propagation in elastic fiber composites

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.


Introduction
Unidirectional fiber-reinforced composites find important application in numerous industrial areas [1][2][3][4], and the study of their mechanical behavior is a topic of long history [5][6][7][8][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][11][12][13][14][15][16][17][18][19][20][21][22][23][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 fiberreinforced concretes; Kim [21] conducted a comparative study on anti-plane SH-waves in fiberreinforced 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.

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: ) aligned along the x 1 -axis and randomly distributed with a constant volume fraction d, 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.
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 1 -independent (plane or anti-plane) displacement fields: the displacement field u f (x 2 , x 3 , 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 2 , x 3 , 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 e of the composite is calculated based on the composite's displacement u(x 2 , x 3 , t), and the three-dimensional (3D) transversely isotropic strain e/stress s relations of the composite can be obtained from the effective composite model given in Appendix 1.The mass density r of the composite (per unit volume of the composite) is defined by, where d is the volume fraction of the embedded fibers, and r m and r 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: where s 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, 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 fu) 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 6 ¼ u on the two sides of equation (2).
Based on the spring model for the x 1 -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 2 , x 3 )) as follows: where the spring constant b (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 1 -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.

Longitudinal P-waves
First, let us consider the plane P-waves under the conditions of plane strain with u(x 2 , x 3 , t) and v(x 2 , x 3 , t) and w [ 0 (see Figure 1).In terms of the two effective elastic coefficients (c 22 , c 23 ), the stress-strain relations in plane-strain conditions are given by (see Appendix 1), where, For plane P-waves propagating in x 2 -direction, we have u(x 2 , t) with v = 0, thus . Substituting it into the above relationship equation ( 4), we have, Thus, for harmonic P-waves u(x 2 , t) =u 0 e -i(kx , where u 0 is a constant, v is circular frequency, and k is the wave number, then, Within the following range of frequencies (called the ''bandgap'' for P-waves and SV-waves under plane-strain, where k = 3-4n 23 , see Appendix 2), We have, 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 v 0 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][11][12][13][17][18][19][20][21]) are illustrated in Figures 2(a) and 3(a), respectively.Clearly, the actual attenuation at the local resonance frequency v = v 0 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 v = v 0 .In spite of this, however, the results given by the present model (shown in Figures 2(b On the other hand, outside the bandgap equation ( 9), the present model gives that the attenuation a = 0 and the longitudinal phase velocity c P for P-waves decreases with increasing v in two separate ranges, below or above the bandgap respectively, and given by, In particular, for low frequencies well below the bandgap, the phase velocity is given approximately by, For high frequencies well above the bandgap, we have, 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 v = v 0 .A more detailed quantitative comparison of the present model with known data on P-waves will be made in section 4.

Transverse SV-waves
Similarly, for transverse SV-waves propagating in x 2 -direction, we have v(x 2 , t) and u = 0, thus with a constant v 0 , substituting it into the relationship equation ( 4) gives, Thus, we have, 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, On the other hand, outside the bandgap, we have a = 0, and the transverse phase velocity c SV for SVwaves decreases with increasing v in two separate ranges, below or above the bandgap, respectively, In particular, for low frequencies well below the bandgap, the phase velocity c SV is given approximately by, For high frequencies well above the bandgap, we have, 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.

Anti-plane SH-waves
For anti-plane deformation w(x 2 , x 3 , t) with u = v = 0, the stress-strain Hookean law gives (see Appendix 1), For SH-wave propagating along x 2 -direction, we have w(x 2 , t) =w 0 e -i(kx 2 -vt) with a constant w 0 , thus s 12 = G 12 w , 2 ; s 13 = 0, and r(∂ 2 w mass ∂t 2 ) = G 12 w , 22 .Substituting it into the relationship equation ( 9) gives, Thus, we have, Within the following range of frequencies (called the ''bandgap'' for SH-waves under anti-plane shearing, see Appendix 2) defined by, We have, On the other hand, outside the bandgap, we have a = 0, and the phase velocity c SH for SH-waves decreases with increasing v in two separate ranges, below or above the bandgap, respectively, In particular, for low frequencies well below the bandgap, the phase velocity c SH is given approximately by, For high frequencies well above the bandgap, we have, 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.

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 v 0 predicted by equation ( 9) for P/SV-waves and the local resonance frequency v w predicted by equation ( 23) for SH-waves will be compared to the known data on the frequency for the maximum attenuation.

The frequency for the maximum attenuation (a) max
The local resonance frequency v 0 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 v 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.
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 v 0 or v 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 (a) 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.

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 (a = 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 (v = 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.  in Sato and Shindo [11] (av/c Pm ) '

2.3
Figure 6 in Varadan et al. [10] Table 2.The local resonance frequency v w predicted by the present model equation ( 23) with comparison to known simulation data on the frequencies for the maximum attenuation coefficient (a) max for SH-waves reported in several studies [17][18][19][20][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 cited in reference  Table 4. Effective phase velocity of SH-waves at the long-wavelength limit (v = 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 3.8G The value given by present model Table 3. Effective phase velocity of P/SV-waves at the long-wavelength limit (v = 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]  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 v = 0 to v = v 0 /2 (for P/SV-waves) or v = v w /2 (for SHwaves), the ratio of the decrease in phase velocity to the phase velocity at v = 0 is given by dr f /(8r).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.

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.
A unidirectional fiber-reinforced transversely isotropic composite shown in Figure 1 where, and, Since transversely isotropic (such as carbon) fibers are widely used in industries, let us consider an isotropic matrix (of the moduli E m and G m , and Poisson ratio n m with G m = E m /[2(1 + n m )]) reinforced by transversely isotropic fibers defined by five independent elastic constants ( ).In existing literature, the two effective elastic constants (E 1 , n 12 ) of the composite can be calculated with reasonable accuracy by the simple rule of mixture [5][6][7][8][9]: The other three effective elastic constants (E 2 , G 12 , G 23 ) can be calculated with reasonable accuracy (as confirmed in Tables 3 and 4) by the ''Chamis formulas'' [3,[6][7][8][9]: and n 23 can be calculated based on the mentioned relation ), E 1 of the composite is given by equation ( 31) with E 1 f = E f while E 2 is given by equation ( 32) with E 2 f = E f , and G 12 = G 23 of the composite is given by equation (32) with

Appendix 2
The spring constant b The spring constant b for the u mass -u relation equation ( 4) is given here for plane strain-stress case and the anti-plane case, respectively.
Plane strain-stress.Let us consider a single rigid fiber surrounded by the infinite composite of effective shear modulus G 23 and Poisson ratio n 23 (see Figure 1 and Appendix 1).The interaction force (per unit length along the x 1 direction) between the fiber and the composite is linearly proportional to the relative displacement (uu f ) in the x 2x 3 plane with a spring constant b 0 (of the dimension of force/ length 2 , per unit length in the x 1 direction) [25]; thus, we have, Eliminating u f from equations (33) and (3), we have the simple differential relationship equation ( 4) between the displacement of the mass center (u mass , v mass ) and the composite's displacement (u, v) in the x 2x 3 plane, with b = b 0 d=a f .An expression for the spring constant b 0 for plane strain-stress is given in Ru [26] as, b 0 ' where, k = 3 À 4n 23 , plane À strain; Anti-plane.Consider (axisymmetric) anti-plane deformation of the infinite composite (of effective shear modulus G 12 ) cause by an x 1 -directional displacement w 0 of a rigid circular fiber (of radius a) perfectly bounded to the composite (see Figure 1 and Appendix 1).The x 1 -directional anti-plane displacement w(r) (r ø a) outside the fiber meets the Laplace equation and the boundary condition w(a) = w 0 and is given by (where C is a constant) [27], w r ð Þ = w 0 + Cln r a , r ø a: Thus the x 1 -directional force (per unit length in x 1 -direction) requested for such a displacement of the fiber is, For the present problem of fiber composites with at least moderately densely packed fibers (say, the average distance d between the centers of adjacent fibers is not larger than five times the diameter of fibers, for example d ø 0.05), the spring constant can be defined based on the relative displacement of the fiber with respect to a reference point of the distance r = d [26], therefore the spring constant is given by, Thus, the single fiber is governed by the spring model with the spring constant b w (of the dimension of force/length 2 , per unit length in the x 1 direction) as, Eliminating w f from equations (39) and (3), we have the simple differential relationship equation (4) between the displacement of the mass center w mass and the composite's displacement w, with b = b w d=a f .
) and 3(b)) are qualitatively consistent with actual wave behavior (shown in Figures2(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 v = v 0 .

Figure 2 .
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 .
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.

Table 1 .
[10][11][12][13]nce frequency v 0 predicted by the present model equation(9)with comparison to known simulation/experimental data on the frequency for the maximum attenuation coefficient (a) max in P/SV-waves reported in several studies[10][11][12][13](where all fibers are isotropic, with the effective Poisson ratio n Pm and cSm are the speeds of longitudinal/shear waves in the isotropic matrix). 4G