Attenuation and Phase Velocity of Elastic Wave in Textured Polycrystals with Ellipsoidal Grains of Arbitrary Crystal Symmetry

Abstract

This study extends the second-order attenuation (SOA) model for elastic waves in texture-free inhomogeneous cubic polycrystalline materials with equiaxed grains to textured polycrystals with ellipsoidal grains of arbitrary crystal symmetry. In term of this work, one can predict both the scattering-induced attenuation and phase velocity from Rayleigh region (wavelength >> scatter size) to geometric region (wavelength << scatter size) for an arbitrary incident wave mode (quasi-longitudinal, quasi-transverse fast or quasi-transverse slow mode) in a textured polycrystal and examine the impact of crystallographic texture on attenuation and phase velocity dispersion in the whole frequency range. The predicted attenuation results of this work also agree well with the literature on a textured stainless steel polycrystal. Furthermore, an analytical expression for quasi-static phase velocity at an arbitrary wave propagation direction in a textured polycrystal is derived from the SOA model, which can provide an alternative homogenization method for textured polycrystals based on scattering theory. Computational results using triclinic titanium polycrystals with Gaussian orientation distribution function (ODF) are also presented to demonstrate the texture effect on attenuation and phase velocity behaviors and evaluate the applicability and limitation of an existing analytical model based on the Born approximation for textured polycrystals. Finally, quasi-static phase velocities predicted by this work for a textured polycrystalline copper with generalized spherical harmonics form ODF are compared to available velocity bounds in the literature including Hashin–Shtrikman bounds, and a reasonable agreement is found between this work and the literature.

1. Introduction

Polycrystals or polycrystalline materials are solid aggregates of numerous individual crystallites (or grains) of varying morphology, size and crystallographic orientation, where the grains are bonded together by grain boundaries. Polycrystals thus are heterogeneous materials in grain scale. The majority of inorganic solids such as common metals, ceramics and rocks are polycrystals. Due to their significance in various industries, polycrystals are one of the major objects of study in non-destructive evaluation (NDE) and seismology. One significant but challenging topic related to polycrystals in NDE and seismology is the elastic wave scattering and the resulting attenuation in polycrystals. It has been studied for decades [1,2,3,4,5,6,7,8] but elastic wave scattering in complex polycrystals is still not well understood.

Elastic waves propagating in heterogeneous polycrystals are subject to scattering at the grain boundaries due to the misorientation between adjacent grains [9,10,11,12], and this further induces beam attenuation due to energy dissipation from the main ultrasonic beam. Although there are various wave attenuation mechanisms in polycrystalline materials such as viscosity [13] and relaxation [14], in many cases the main source of attenuation is wave scattering [5,11] and it has been experimentally proved by a study [15]. The scattering-induced attenuation in polycrystalline materials has an important implication of microstructure information [10,16,17] and thus has been extensively applied to the nondestructive characterization of grain microstructures [18,19]. It is worthy to mention that many studies [20,21,22,23] recently emerge on finite element modeling (FEM) of elastic wave attenuation in polycrystals. However, FEM is a numerical method and it has limited practical application to the ultrasonic characterization of grain microstructures. Due to the industrial needs of ultrasonic characterization of microstructure characteristics from ultrasound measurements, many studies [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23] have strived to develop analytical attenuation models for polycrystalline materials with a variety of complexities such as the low symmetry constituents [7,24,25,26,27], the duplex microstructure [10] and the columnar grains [28].

Numerous attenuation models have been reported for statistically isotropic polycrystals and they fall into four different catalogues: the Karal and Keller model [29], perturbation theory [3,4,30], Weaver type model [5], and the spectral function method [31,32]. The Karal and Keller model [29] is limited to scalar wave attenuation while the Stanke and Kino model [3], a unified model based on perturbation theory and suitable for the whole frequency range, was constrained to macroscopically isotropic polycrystals with equiaxed cubic grains. It is worth mentioning that the Stanke and Kino model [3] has been extended to cubic polycrystals with ellipsoidal shape grains [28]. Weaver [5] systematically derived another attenuation model for elastic waves in texture-free polycrystals by Dyson equation. The Dyson equation in Ref. [5] is general for all polycrystals but the successive derivations were based on the Born approximation (applicable frequency is below the geometric region where wavelength << scatter size) and its explicit expressions were limited to cubic polycrystals without crystallographic texture. A recent study [33] also targeted on texture-free polycrystals has shown that Weaver’s model [5] can be converged to the Stanke and Kino model [3] using numerical iteration. Cavlet and Margerin [32] proposed a spectral function method to model the wave attenuation and phase velocity in texture-free polycrystalline materials at all frequency regions; their development [32] is limited to orthotropic symmetry crystallites [32]. Cavlet and Margerin [32] also drew a comparison to the Stanke and Kino model for cubic iron polycrystals, and it indicates that their model [32] is comparable to the Stanke and Kino model [3]. An intriguing finding in Ref. [32] is that two branches exist for attenuation and phase velocity while Stanke and Kino [3] only calculated one branch. Recently, Cavlet and Margerin [34] further investigated the grain shape effect on attenuation and phase velocity using the spectral function method for texture-free cubic polycrystals. All these studies [3,5,32,33,34] have developed second-order attenuation models, since they consider some multiple scattering events in the mean Green function [35]. Another second-order attenuation model, a far-field attenuation model [36], following Weaver’s approach [5], has been developed for macroscopically isotropic triclinic polycrystals of ellipsoidal grains [36] and it is reasonably close to the Stanke and Kino model [3] and spectral function method [32]. A recent study [23] has reported a full-wave second-order attenuation (SOA) model for texture-free cubic polycrystals by extending Weaver’s approach [5], and the SOA model has a reasonable agreement with Stanke and Kino model [3] and finite element modeling [23]. All these analytical models [3,5,23,29,32,34,36] above are only applicable for macroscopically isotropic polycrystalline materials; however, crystallographic texture generally exists in natural and man-made materials such as rocks and metals causing macroscopically anisotropic properties [37]. Thus, more general attenuation models are desired to take account of crystallographic texture.

Modeling of elastic wave attenuation in textured polycrystals has also attracted the attention of some researchers. For polycrystals with preferred crystallographic orientation, the Weaver model [5], Hireskorn model [4,30], and the Stanke and Kino model [3] have been extended to textured polycrystalline materials with crystallites of different symmetry classes and shapes by Refs. [6,38,39,40,41]. However, the models in Refs. [6,39,40,41] that presented explicit attenuation coefficients by the Born approximation break at high frequency, and the accuracy of the Born approximation has never been investigated for textured polycrystals. These models in Refs. [6,39,40] also fail to provide any information about the phase velocity in textured polycrystalline media that is of practical importance in ultrasonic nondestructive testing and seismology (analytical expressions for velocities in Ref. [41] are not applicable to high-frequency regions, e.g., stochastic regime). Note that the model developed in Ref. [38] based on the Stanke and Kino model [3] adopted the Green function for isotropic reference media instead of that for anisotropic reference media, making it an unrigorous model although it was a significant development in the 1990s. The latest development of attenuation modeling in Ref. [40] accounted for aggregates of triclinic grains in ellipsoidal shape with arbitrary texture symmetry and obtained explicit equations for attenuation coefficients; however, it was implemented under the frame of the Born approximation [40] and had no attempt to predict phase velocity behavior in textured polycrystal due to the complexity of mass operator (self-energy operator) [5] in anisotropic media. Therefore, the attenuation and phase velocity behaviors of elastic waves in the textured polycrystals have not been well understood.

This work aims to develop a second-order attenuation model for a textured polycrystalline aggregate of triclinic crystallites in ellipsoidal shape by extending an existing study [23] for texture-free polycrystals. It also investigates the accuracy and applicability of the Born approximation, which will define the boundary between single-scattering theory [12] and diffusion theory [12], and the behaviors of phase velocity in different frequency regimes. The work will degenerate to previous analytical models [5,23] when the crystallographic texture of a polycrystal is extremely weak. Although a recent paper has studied grain number impact on attenuation [42] by synthetic grain microstructures, this work still focuses on polycrystals composed of infinite grains or a sufficiently large number of grains (so that grain number impact is negligible), just as most of the previous studies. Several novelties exist in this work compared to previous studies [23,40]: it is the first time integrating the Green function for anisotropic reference media into an attenuation model, implementing numerical calculation of the full mass operator (including Cauchy Principal value in the real part) and solving the dispersion equation for a textured polycrystal of triclinic grains; the SOA model for textured polycrystals is further validated by comparing to available FEM results and experimental data on a textured stainless steel sample; to our best knowledge, it is also the first time to investigate the accuracy and limitation of the Born approximation on textured polycrystals and examine the texture impact both on attenuation and phase velocity at whole frequency range; as another original development of this work, an explicit expression for quasi-static velocities at an arbitrary wave propagation direction in textured polycrystals of triclinic elongated grains is derived in this paper and obtained quasi-static velocities from this work are compared to available bounds in the literature on a polycrystalline copper.

This paper is structured as follows. The relevant theoretical background is first reviewed in Section 2. It is followed by the SOA model for textured polycrystals and its asymptotes for attenuation and phase velocity at the Rayleigh region. In Section 4, comparisons are drawn between the SOA model developed in this work and other studies on textured polycrystals to justify the SOA model and examine the limitation of the Born approximation. The texture impacts on attenuation and phase velocity are also studied. The quasi-static velocities from the SOA model are further compared on a polycrystalline copper with available velocity bounds in the literature. Finally, conclusions are provided.

2. Theoretical Background about Attenuation Modelling in Textured Polycrystals

This section will introduce the theoretical background of attenuation modeling for textured polycrystalline materials. Firstly, a quantitative texture description, namely orientation distribution function (ODF), is introduced. Furthermore, based on the Weaver type model [5] some fundamentals of attenuation modeling for polycrystalline materials are reviewed. Finally, one key element of attenuation modeling, the two-point correlation function, is addressed.

2.1. Orientation Distribution Function (ODF)

Statistical isotropy is the ideal status of polycrystalline materials and crystallographic texture widely exists in man-made and natural polycrystalline materials such as metals and minerals [7,43,44]. In other words, the crystallographic orientations of crystallites in an actual polycrystals aggregate are not perfectly randomized; instead, grain orientations are aligned in certain directions and the whole aggregate thus presents anisotropic properties.

2.1.1. Generalized Spherical Harmonics ODF

To bridge the gap between the property of individual crystallites and that of the assembly, one needs to know the crystallographic orientation distribution of crystallites that presents the probability density of a given orientation direction in 3D space. Intensive studies were reported to quantify the orientation distribution of crystallites in the past [43,45,46,47]. Well-known terminologies of orientation distribution functions (ODF) for quantitative texture analysis are from Bunge [43] and Roe [46]. Both are in terms of generalized spherical harmonics, although different rotation angles are used.

Here we follow the ODF in Bunge’s terminology [43]. To describe the three-dimensional ODF by generalized spherical harmonics, the representation of rotation from crystallite coordinate to global coordinate should be well defined. In Bunge’s notation [43], the three Eulerian angles describing 3D rotation are φ₁, Φ, φ₂, which are shown in Figure 1. To rotate the global texture coordinate to the crystallite coordinate, one has to follow such a procedure [43]: first, rotate about the T_Z axis through the angle φ₁; then rotate about x’ axis by the angle Φ and finally rotate about z-axis by the angel φ₂. Therefore, the rotation matrix in Bunge’s Euler rotation angles is:

$$R\left({\phi }_1,\Phi ,{\phi }_2\right)=\left[\begin{array}{ccc}\cos {\phi }_1\cos \Phi \cos {\phi }_2-\sin {\phi }_1\sin {\phi }_2& \sin {\phi }_1\cos \Phi \cos {\phi }_2+\cos {\phi }_1\sin {\phi }_2& -\sin \Phi \cos {\phi }_2\\ -\cos {\phi }_1\cos \Phi \sin {\phi }_2-\sin {\phi }_1\cos {\phi }_2& -\sin {\phi }_1\cos \Phi \sin {\phi }_2+\cos {\phi }_1\cos {\phi }_2& \sin \Phi \sin {\phi }_2\\ \cos {\phi }_1\sin {\phi }_2& \sin {\phi }_1\sin \Phi & \cos \Phi \end{array}\right]$$

(1)

The ODF in Bunge’s notation expressed by generalized spherical harmonics from Ref. [43] is:

$$f({\phi }_1,\Phi ,{\phi }_2)={\displaystyle \sum _{l=0}^{\infty }{\displaystyle \sum _{m=-l}^l{\displaystyle \sum _{n=-l}^l{C}_l^{mn}\exp (im{\phi }_2)P_l^{mn}(\Phi )\exp (in{\phi }_1)}}}{,}_{ }$$

(2)

where $C_l^{mn},P_l^{mn}(\Phi )$ are texture coefficients and polynomial function [43], and $l,m,n$ are integers (0 ≤ l ≤ ∞; however, l ≤ 32 in practical texture measurements). The polynomial function $P_l^{mn}(\Phi )$ is defined in Ref. [43] as:

$$\begin{array}{l}{P}_l^{mn}(\Phi )=P_l^{mn}(\cos \Phi )=P_l^{mn}(x)=\\ \frac{{(-1)}^{l-m}{i}^{n-m}}{2^l(l-m)!}{\left[\frac{(l-m)!(l+n)!}{(l+m)!(l-n)!}\right]}^{1/2}(1-x{)}^{-\frac{n-m}{2}}(1+x{)}^{-\frac{n+m}{2}}\frac{d^{l-n}}{d{x}^{l-n}}\left[{(1-x)}^{l-m}{(1+x)}^{l+m}\right],\end{array}$$

(3)

where x is defined as x=cosΦ for simplicity. It represents a series of mutually orthogonal bases since this relation holds [43]:

$${\displaystyle {\int }_0^{\pi }{P}_l^{mn}(\Phi )}{\left(P_{l\prime }^{mn}(\Phi )\right)}^{*}\sin \Phi d\Phi =\frac{2{\delta }_{ll\prime }}{2l+1}.$$

(4)

2.1.2. Gaussian form ODF

Although the generalized spherical harmonics form ODFs [43,46] are widely used in material science and engineering, one should note that there are other forms of ODF existing in the literature [47,48,49,50,51,52]. For example, Gaussian-shape ODF as reported by references [47,48,49,50] has been utilized as a simple tool to investigate the texture effect on piezoelectricity and elasticity. As described in Ref. [50], the 3D Gaussian shape ODF is:

$$F\left(\Phi ,{\phi }_1,{\phi }_2\right)=F_0\exp \left(\frac{\cos \Phi }{2{\sigma }_{\Phi }}+\frac{\cos {\phi }_1}{2{\sigma }_{\phi 1}}+\frac{\cos {\phi }_2}{2{\sigma }_{\phi 2}}\right),$$

(5)

where σ_φ1, σ_Φ, σ_φ2 are texture parameters and the factor F₀=S₁csch(S₁)/[I(0,S₂)I[0,S₃]] satisfies this equality $\frac{1}{8{\pi }^2}{\int }_0^{\pi }{\int }_0^{2\pi }{\int }_0^{2\pi } F\left(\Phi ,{\phi }_1,{\phi }_2\right)sin\Phi d\Phi d{\phi }_{1}d{\phi }_2=1$. Here, S₁=1/(2σ_Φ), S₂= 1/(2σ_φ1), S₃=1/(2σ_φ2) and I(0,S₂), I(0,S₃) are modified Bessel functions of the first kind [50], and φ₁, Φ, φ₂ are three Euler angles defined by Bunge notation as shown in Figure 1. In this representation, the texture parameters σ_φ1, σ_Φ, σ_φ2 in the Gaussian ODF can be directly related to the full width at half maxima (FWHMs) of texture angles φ₁, Φ, φ₂ from experimental orientation imaging microscopy [39,47]. When all three parameters σ_φ1, σ_Φ, σ_φ2 approach infinity the ODF degenerates to the case of statistically isotropic polycrystals. When two of them go to infinity, the ODF would be equivalent to fiber (axisymmetric) texture like Refs. [39,53]. Note that ODFs in different notations are related to each other, for example, the relation between Bunge’s ODF and Roe’s ODF has been investigated in Ref. [54] and the conversion from 1D Gaussian ODF to Roe’s ODF has been proposed in Refs. [49,53].

2.2. Theoretical Background of Attenuation Modelling

Although there are three types of theories available for modeling elastic wave attenuation (perturbation method [3,4,30], spectral function method [32], and Dyson equation [5]), here we only focus on the theoretical background of Dyson equation method [5].

As stated in Refs. [5,6,40], the elastodynamic response at a location $\mathit{X}$ to a point load in inhomogeneous media at $X\prime$ is (all the equations below in the paper are in tensor nation and follow the Einstein summation):

$$\left\{\frac{\partial }{\partial x_i}\delta c_{ijkl}\left(X\right)\frac{\partial }{\partial x_l}+C_{ijkl}^0\frac{\partial }{\partial x_i}\frac{\partial }{\partial x_l}+{\omega }^2{\delta }_{jk}\right\}{G}_{k\alpha }\left(X,X^{\prime };\omega \right)={\delta }_{j\alpha }{\delta }^3\left(X-X^{\prime }\right),$$

(6)

where $\delta c_{ijkl}\left(\mathit{X}\right)$, a fourth-order tensor, is the variation of normalized elastic constants from the reference medium, δ_jα means Kronecker delta (a second-order tensor) and $C_{ijkl}^0$ stands for the homogenized elastic constants. ω is the angular frequency. The elastic constant at a given location X is equal to the elastic constants of its reference medium plus a spatial fluctuation term, namely:

$$c_{ijkl}\left(X\right)=C_{ijkl}^0+\delta c_{ijkl}\left(X\right),$$

(7)

where the ensemble average of elastic constants is $C_{jkli}^0=\langle c_{ijkl}\left(\mathit{X}\right)\rangle$ and the symbol < > denotes ensemble average over all grains. Note that the reference media is anisotropic for textured polycrystals. It implies that the statistical average of the spatial fluctuation term is zero $\langle \delta c_{ijkl}\left(\mathit{X}\right)\rangle =0$ from Equation (7), as explained in Refs. [5,40].

As stated in Refs. [5,6], the solution of the elastodynamic equation at a certain location is unsolvable, but the mean Green function response $\langle G_{k\alpha }\left(\mathit{X},{\mathit{X}}^{\prime }\right)\rangle$ in the heterogeneous polycrystals can be obtained by the Dyson equation:

$$\langle G_{k\alpha }\left(X,X^{\prime }\right)\rangle =G_{k\alpha }^0\left(X,X^{\prime }\right)+{\displaystyle \iint G_{k\beta }^0\left(X,Y\right)M_{\beta j}\left(Y,Z\right)\langle G_{j\alpha }\left(Z,X^{\prime }\right)\rangle d^3\mathbf{y}{d}^3\mathbf{z}},$$

(8)

where $G_{k\alpha }^0\left(\mathit{X},{\mathit{X}}^{\prime }\right)$ is the Green function of a homogenized medium or reference medium. $M_{\beta j}\left(\mathit{Y},\mathit{Z}\right)$ is the mass operator or self-energy accounting for all possible scattering events [55] between location Y and location Z. The mass operator may be expressed in a diagrammatic form including infinite series [56]. However, the exact equation for the mass operator is difficult to obtain but finite order approximation is employed instead in references [35,56]. Under the weak scattering assumption $\left|\delta c_{ijkl}\right|/C_{ijkl}^0<<1$, the mass operator after first order smoothing approximation (FOSA) is [6]:

$$M_{\beta j}\left(Y,Z\right)\approx \langle \frac{\partial }{\partial Y_{\alpha }}\delta c_{\alpha \beta \gamma \delta }\left(Y\right)\frac{\partial }{\partial Y_{\delta }}{G}_{\gamma q}^0\left(Y,Z\right)\frac{\partial }{\partial Z_i}\delta c_{ijql}\left(Z\right)\frac{\partial }{\partial Z_l}\rangle ,$$

(9)

where $G_{\gamma q}^0$ is the Green function for the homogenized medium and δc_ijql is the spatial variation of elastic constants. The mass operator in Equation (9) includes some second-order scattering events according to Refs. [35,56]. From Equation (9), the mass operator is relevant to the covariance of elastic constants fluctuation at two different points, namely the two-point correlation [5].

After spatial Fourier transform of the double convolution, the Dyson equation can be simplified as [5]:

$$\langle G\left(k\right)\rangle =G_{}^0\left(k\right)+G_{}^0\left(k\right)M\left(k\right)\langle G\left(k\right)\rangle .$$

(10)

Since the Green function $G_{k\alpha }^0\left(\mathit{X},{\mathit{X}}^{\prime }\right)$ in anisotropic reference media satisfies this equation [6,40]:

$$\left\{C_{ijkl}^0\frac{\partial }{\partial x_i}\frac{\partial }{\partial x_l}+{\omega }^2{\delta }_{jk}\right\}{G}_{k\alpha }^0\left(X,X^{\prime }\right)={\delta }_{j\alpha }{\delta }^3\left(X-X^{\prime }\right)$$

(11)

after spatial Fourier transform on both side of the equation above one obtains:

$$\left(-C_{ijql}^0{k}_i{k}_l+{\omega }^2{\delta }_{jq}\right)G_{q\alpha }^0\left(k\right)={\delta }_{j\alpha }.$$

(12)

The tonsorial form (dyadic) Green function in Equation (12) of an anisotropic homogenized medium can be decomposed into three components according to the spectral method (matrix decomposition) [40,57]:

$$G_{}^0\left(k\right)={\displaystyle \sum _{M=1}^3\frac{u_{}^M\otimes u_{}^M}{(\Omega -{\Omega }_M)}}={\displaystyle \sum _{M=1}^3\frac{u_{}^M\otimes u_{}^M}{{\omega }^2-k^2{V}_M^2}}={\displaystyle \sum _{M=1}^3{g}_M^o{u}_{}^M\otimes u_{}^M},$$

(13)

where ${\Omega }_M=k^2{V}_M^2$ are eigenvalues (M is a wave mode in anisotropic media, quasi-longitudinal, quasi-transverse fast or quasi-transverse slow mode), p is the unit wave normal vector, and $\otimes$ denotes dyadic product. The phase velocity V_M and polarization vector ${\mathit{u}}^M$ at a certain wave propagation direction are determined through Christoffel’s Equation [40]:

$$[C_{ijml}^0{p}_j{p}_l-\rho V_M^2{\delta }_{im}^{}]u_m=0,$$

(14)

where p_j is a component of wave normal vector p and ρ denotes the mass density of the medium. The phase velocity V_M is the square root of the eigenvalue of the matrix $\frac{C_{ijml}^0{n}_j{n}_l}{\rho }$ for wave mode M and polarization vector u is the corresponding eigenvector. In a textured polycrystal, there are three wave modes: quasi-longitudinal (L) mode, quasi-transverse fast (T1) mode and quasi-transverse slow (T2) mode. Unlike the isotropic media, the polarization direction of a quasi-longitudinal wave generally deviates from its wave propagation direction while the polarization directions of a quasi-transverse wave are not perpendicular to its wave propagation direction. However, three polarization vectors are still mutually perpendicular to each other.

From Equation (10), the mean Green function $\langle \mathit{G}\left(\mathit{k}\right)\rangle$ in the spatial Fourier domain solved from the Dyson equation is [6,40]:

$$\langle G\left(k\right)\rangle ={\left[G_{}^0{\left(k\right)}^{-1}-M\left(k\right)\right]}^{-1},$$

(15)

where $\mathit{M}\left(\mathit{k}\right)$ is the spatial Fourier transforms of the mass operator. Before we proceed to the expression for the mass operator, two-point correlation function should be introduced first.

2.3. Two-point Statistics for Textured Polycrystals

The mass operator aforementioned is related to the two-point correlation function that is expressed as δc_ijkl(X)δc_αβγδ(X’), a covariance of elastic constants fluctuation from two points X and X’ (see the schematic of a polycrystal in Figure 2). Since we are solely interested in the mean wave response in heterogeneous polycrystals, the ensemble average of this covariance is critical. From statistics of numerous grains, the ensemble average of elastic covariance can be further decomposed into two parts: the volumetric average of elastic constants variation covariance and geometric two-point correlation (GTPC) function [6,40]. The mathematical expression is:

$$\langle \delta c_{ijkl}(X)\delta c_{\alpha \beta \gamma \delta }(X^{\prime })\rangle =\langle \delta c_{ijkl}\delta c_{\alpha \beta \gamma \delta }\rangle w\left(X-X^{\prime }\right),$$

(16)

where <δc_ijkδc_αβγδ> is the ensemble average of elastic constants variation covariance for textured polycrystals by Voigt averaging [58]. It can be numerically calculated through single-crystal elastic constants and ODF by the following equation [40,50,59]:

$$\langle \delta c_{ijkl}\delta c_{\alpha \beta \gamma \delta }\rangle =\langle c_{ijkl}{c}_{\alpha \beta \gamma \delta }\rangle -\langle c_{ijkl}\rangle \langle c_{\alpha \beta \gamma \delta }\rangle ,$$

(17)

where the effective elastic constants by Voigt averaging and the elastic constants covariance for a polycrystals aggregate are [40,50,59]:

$$\begin{array}{l}{C}_{ijkl}^0=\langle c_{ijkl}\rangle =\frac{1}{8{\pi }^2}{\displaystyle \underset{0}{\overset{\pi }{\int }}{\displaystyle \underset{0}{\overset{2\pi }{\int }}{\displaystyle \underset{0}{\overset{2\pi }{\int }}{c}_{ijkl}({\phi }_1,\Phi ,{\phi }_2)f({\phi }_1,\Phi ,{\phi }_2)\sin (\Phi }}})d{\phi }_{1}d{\phi }_{2}d\Phi ;\\ \langle c_{ijkl}{c}_{\alpha \beta \gamma \delta }\rangle =\frac{1}{8{\pi }^2}{\displaystyle {\int }_0^{\pi }{\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_0^{2\pi } }}}{c}_{ijkl}({\phi }_1,\Phi ,{\phi }_2)c_{\alpha \beta \gamma \delta }({\phi }_1,\Phi ,{\phi }_2)f({\phi }_1,\Phi ,{\phi }_2)\sin \Phi d{\phi }_{1}d{\phi }_{2}d\Phi .\end{array}$$

(18)

Here c_ijkl(φ₁, Φ, φ₂) and f(φ₁, Φ, φ₂) are the rotated single-crystal elastic tensor and ODF (either Gaussian ODF or generalized spherical harmonics ODF), respectively. The elastic constants rotation matrix and details about elastic constants covariance calculation can be found in Refs. [39,40,50,53,59]. It is worthy to mention that Voigt effective elastic constants (or called Voigt reference media) are used in this whole paper because it is more suitable for wave scattering modeling [23,60], although different homogenization methods like self-consistent reference medium are available for polycrystalline media [58,61]. Ref. [23] has compared the SOA models (for texture-free polycrystals) with Voigt reference medium and self-consistent reference medium to finite element modeling and found that the SOA model with Voigt reference medium agrees much better with FEM results.

The spatial two-point correlation function $w\left(\mathit{X}-{\mathit{X}}^{\prime }\right)$ from Equation (16) only depends on the geometric characteristics of grain microstructure. Physically, it addresses the probability that two points $X$ and $X^{\prime }$ fall into the same grain. In a polycrystal aggregate shown in Figure 2, the probability of two points with a shorter distance falling into the same grain is always higher than that of two points departing by a larger distance. For a realistic polycrystalline material, its spatial two-point correlation function can be determined by the line intercept method from numerous microstructure cross-sections [62]. Under the assumption that line intercepts on a cross-section follow a Poisson distribution, an exponential form GTPC function was initially proposed by Stanke [62] and it has been adopted by many other publications [6,7,9,24,40,48,50,59,63]. The mere parameter governing the exponential GTPC is the mean line intercept (or called cord length) [62].

In this paper, it is still assumed that the geometric two-point correlation function is in exponential form but the polycrystals contain ellipsoidal shape grains with three (statistically mean) semi-principal axes $a_X,a_Y$ and $a_Z$ (see Figure 3). According to Refs. [9,24,40,50,59,64,65], the spatial two-point correlation function W(r) for general ellipsoidal grains is:

$$W\left(r\right)=\exp (-\sqrt{\frac{x^2}{a_X^2}+\frac{y^2}{a_Y^2}+\frac{z^2}{a_Z^2}}).$$

(19)

Considering an incident wave propagating in direction $\mathit{p}$ and a scattered wave in direction $\mathit{s}$ (shown in Figure 3), according to Refs. [9,24,40,50,59,64], the two-point correlation function in spatial Fourier domain W (q) is:

$$W\left(q\right)=\frac{a_X{a}_Y{a}_Z}{{\pi }^2{\left[1+\left(q_X^2{a}_X^2+q_Y^2{a}_Y^2+q_Z^2{a}_Z^2\right)\right]}^2},$$

(20)

where $\mathit{q}=\mathit{k}-{\mathit{k}}_{}^{S,N}$ related to the incident wavenumber k and the scattered wavenumber ${\mathit{k}}^{S,N}$. Here N stands for a scattering wave mode in an anisotropic medium, namely a quasi-longitudinal mode, a quasi-transverse fast mode or a quasi-transverse slow mode.

3. Second-order Scattering Model for Textured Polycrystals with Ellipsoidal Triclinic Grains

This section will develop the SOA model for textured polycrystals with triclinic grains of ellipsoidal shape. The asymptotes at Rayleigh limit for attenuation and phase velocity will also be derived.

3.1. Dispersion Equation

As stated in Ref. [40], the mass operator after spatial Fourier transform is:

$$M_{\beta j}\left(k\right)=k_{\alpha }{k}_l{\displaystyle \int G_{\gamma q}^0\left(k_{}^S\right)k_{\delta }^S{k}_i^S\langle \delta c_{ijql}\delta c_{\alpha \beta \gamma \delta }\rangle W\left(k-k^S\right)d^3{k}_{}^S},$$

(21)

where $\mathit{k}=k\mathit{p}$ means incident wavenumber and ${\mathit{k}}^S=k^S\mathit{s}$ denotes scattered wavenumber. $\langle \delta c_{ijql}\left(\sigma \right)\delta c_{\alpha \beta \gamma \delta }\left(\sigma \right)\rangle$ is the ensemble average of elastic constant covariance and $W\left(\mathit{k}-{\mathit{k}}^S\right)$ stands for the GTPC function in the spatial Fourier domain. Similar to the Green function in a homogenized medium, the mass operator and the mean Green function can be rewritten into dyadic form via spectral method [36,40]:

$$M\left(k,\omega \right)={\displaystyle \sum _{M=1}^3{m}_M{u}_{}^M\otimes u_{}^M},\langle G\left(k,\omega \right)\rangle ={\displaystyle \sum _{M=1}^3{g}_M{u}_{}^M\otimes u_{}^M},$$

(22)

where $m_M$ means the eigenvalue of ${\mathit{M}}_{\beta j}\left(\mathit{k}\right)$ at a wave mode M and the mean response of an incident wave mode M is [6,40]:

$$g_M(k)={\left[g_M^o{(k)}^{-1}-m_M(k)\right]}^{-1}={\left[{\omega }^2-k^2{V}_M^2-m_M(k)\right]}^{-1}.$$

(23)

Since the denominator of the mean Green function should equal to zero, the dispersion equation is expressed as [6,40]:

$${\omega }^2-k^2(\mathbf{p})V_M^2(\mathbf{p})-m_M(k)=0,$$

(24)

where k=Re(k)+iα is a complex wavenumber in the perturbed medium. The real part provides frequency-dependent velocity while the imaginary part corresponds to the beam attenuation. One should note that $k\left(\mathit{p}\right),V_M\left(\mathit{p}\right),m_M\left(\mathit{k}\right)$ are direction-dependent.

After substituting Equation (22) into Equation (21) and rewriting the tonsorial mass operator via Equation (22), one can obtain the expression for $m_M\left(\mathit{k}\right)$. For a certain incidence wave M, the mass operator is the sum of three scattering terms ( L-L scattering, L-T1 scattering and L-T2 scattering [40]):

$$m_M(k)={\displaystyle \sum _{N=1}^3{m}_{M\to N}}=\frac{k^2(\mathbf{p})}{{\rho }^2}{\displaystyle \sum _{N=1}^3{\displaystyle \int g_0^N(k_{}^{S,N})I{P}_{M\to N}^{}(\mathbf{p},\mathbf{s})k_s^{2}W(k-k_{}^{S,N})}{d}^3{\mathit{k}}_{}^{S,N}},$$

(25)

where the inner product is expressed as [40,50]:

$$I{P}_{M\to N}=\langle \delta c_{ijql}\delta c_{\alpha \beta \gamma \delta }\rangle u_{\beta }{u}_q{p}_{\alpha }{p}_l{s}_i{s}_{\delta }{v}_{\gamma }{v}_j,$$

(26)

where u and v are the polarization vectors corresponding to the incident wave and scattered wave, respectively.

To calculate the integration in mass operator, an explicit branch (upper half-space or lower half-space) of the spectral Green function for anisotropic reference media should be chosen. Here we introduce an infinite small imaginary value ε in Green function to choose the upper half-space, and according to the Sokhotskyi–Plemelj theorem [32,36] the Green function is rewritten as:

$$g_0^N(k^{S,N})=\frac{1}{{(\omega +i\epsilon )}^2-{(k^{S,N}{V}_N^{})}^2}=P.V.\left(\frac{1}{{\omega }^2-{(k^{S,N}{V}_N^{})}^2}\right)-i\pi \delta [{\omega }^2-{(k^{S,N}{V}_N^{})}^2],$$

(27)

where P.V. stands for Cauchy principal value.

After substituting Equation (27) into Equation (25), the components of the mass operator can be rewritten as:

$$\begin{array}{l}{m}_{M\to N}=\frac{k^2(\mathbf{p})}{{\rho }^2}\{P.V.{\displaystyle {\int }_0^{\infty }\frac{k_S^4}{k_N^2(\mathbf{s})-k_S^2}}d{k}_S{\displaystyle {\int }_0^{\pi }{\displaystyle {\int }_0^{2\pi }\frac{I{P}_{M\to N}^{}(\mathbf{p},\mathbf{s})W(k-k_{}^S)}{V_N^2(\mathbf{s})}sin\theta d\theta d\phi }}\\ -\frac{i\pi }{2}{\displaystyle {\int }_0^{\pi }{\displaystyle {\int }_0^{2\pi }\frac{k_N^3(\mathbf{s})I{P}_{M\to N}^{}(\mathbf{p},\mathbf{s})W(k-k_N^S)}{V_N^2(\mathbf{s})}}}sin\theta d\theta d\phi \},\end{array}$$

(28)

where the GTPC function spatial Fourier domain is [40,50,59]:

$$W\left(q=k-k_{}^S\right)=\frac{R_0{R}_1{a}_X^3}{{\pi }^2{\left[1+\left(q_X^2{a}_X^2+q_Y^2{a}_Y^2+q_Z^2{a}_Z^2\right)\right]}^2},$$

(29)

where $R_0=a_Z/a_X$ and $R_1=a_Y/a_Z$.

Equation (28) is composed of nine scattering components, namely nine M-N scattering scenarios. For example, three scattering components of the mass operator corresponding to a quasi-longitudinal incidence wave are:

$$\begin{array}{l}{m}_{L\to L}=\frac{k^2(\mathbf{p})}{{\rho }^2}\{P.V.{\displaystyle {\int }_0^{\infty }\frac{k_S^4}{k_L^2(\mathbf{s})-k_S^2}}d{k}_S{\displaystyle {\int }_0^{\pi }{\displaystyle {\int }_0^{2\pi }\frac{I{P}_{L\to L}^{}(\mathbf{p},\mathbf{s})W(k-k_{}^S)}{V_L^2(\mathbf{s})}sin\theta d\theta d\phi }}\\ -\frac{i\pi }{2}{\displaystyle {\int }_0^{\pi }{\displaystyle {\int }_0^{2\pi }\frac{k_L^3(\mathbf{s})I{P}_{M\to N}^{}(\mathbf{p},\mathbf{s})W(k-k_L^S)}{V_L^2(\mathbf{s})}}}sin\theta d\theta d\phi \}\end{array}$$

(30)

$$\begin{array}{l}{m}_{L\to T1}=\frac{k^2(\mathbf{p})}{{\rho }^2}\{P.V.{\displaystyle {\int }_0^{\infty }\frac{k_S^4}{k_{T1}^2(\mathbf{s})-{(k_S)}^2}d{k}_S{\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_0^{\pi }\frac{I{P}_{L\to T1}^{}(\mathbf{p},\mathbf{s})W(k-k_{}^S)}{V_{T1}^2(\mathbf{s})}}}}\sin \theta d\theta d\varphi \\ -\frac{i\pi }{2}{\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_0^{\pi }I{P}_{L\to T1}^{}(\mathbf{p},\mathbf{s})W(k-k_{T1}^S)}}\frac{k_{T1}^3(\mathbf{s})}{V_{T1}^2(\mathbf{s})}\sin \theta d\theta d\varphi \}\end{array}$$

(31)

$$\begin{array}{l}{m}_{L\to T2}=\frac{k^2(\mathbf{p})}{{\rho }^2}\{P.V.{\displaystyle {\int }_0^{\infty }\frac{k_S^4}{k_{T2}^2(\mathbf{s})-{(k_S)}^2}d{k}_S{\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_0^{\pi }\frac{I{P}_{L\to T2}^{}(\mathbf{p},\mathbf{s})W(k-k_{}^S)}{V_{T2}^2(\mathbf{s})}}}}\sin \theta d\theta d\varphi \\ -\frac{i\pi }{2}{\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_0^{\pi }I{P}_{L\to T2}^{}(\mathbf{p},\mathbf{s})W(k-k_{T2}^S)\frac{k_{T2}^3(\mathbf{s})}{V_{T2}^2(\mathbf{s})}}}\sin \theta d\theta d\varphi \}\end{array}$$

(32)

Similarly, one can find the mass operator components for a quasi-transverse fast (T1) or quasi-transverse slow (T2) incidence wave.

Since analytical solutions to complex wavenumber k are difficult to obtain, the mass operator $m_M\left(\mathit{k}\right)$ including the inner product and the dispersion Equation (24) have to be solved numerically. In this study, all calculations are implemented by FORTRAN code since its IMSL numerical libraries provide subroutines that can compute the Cauchy principal values in the mass operator elegantly. FORTRAN also has relatively high computational efficiency. Compared with the SOA model in Ref. [23] for the texture-free polycrystals, computation of the Cauchy principal value and solving dispersion equation are much more time-consuming for textured polycrystal cases even on a workstation. Empirically, obtaining one data point of perturbed wavenumber may take tens of hours.

3.2. Analytical Equations for Attenuation Coefficient and Phase Velocity

3.2.1. Attenuation Coefficient and Phase Velocity by the Born Approximation

In practical applications such as nondestructive evaluation and seismology, the Born approximation is widely used since it simplifies the interpretation of the acoustic signals [9,10]. Physically, the Born approximation assumes very weak scattering, with a result that the perturbed wavenumber in grains is equal to the unperturbed wavenumber. After applying the Born approximation to the dispersion equation (24), the perturbed complex wavenumber can be explicitly expressed as:

$$k(\mathbf{p})=k_M^{}(\mathbf{p})-\frac{m_M(k_M(\mathbf{p}))}{2k_M^{}(\mathbf{p})V_M^2(\mathbf{p})},$$

(33)

where the mass operator after applying the Born approximation is:

$$\begin{array}{l}{m}_M(k_M(\mathbf{p}))={\displaystyle \sum _{N=1}^3\frac{k_M^2(\mathbf{p})}{{\rho }^2}\{P.V.{\displaystyle {\int }_0^{\infty }\frac{k_S^4}{k_N^2(\mathbf{s})-k_S^2}}d{k}_S{\displaystyle {\int }_0^{\pi }{\displaystyle {\int }_0^{2\pi }\frac{I{P}_{M\to N}^{}(\mathbf{p},\mathbf{s})W(k_M-k_{}^S)}{V_N^2(\mathbf{s})}sin\theta d\theta d\phi }}}\\ -\frac{i\pi }{2}{\displaystyle {\int }_0^{\pi }{\displaystyle {\int }_0^{2\pi }\frac{k_N^3(\mathbf{s})I{P}_{M\to N}^{}(\mathbf{p},\mathbf{s})W(k_M-k_N^S)}{V_N^2(\mathbf{s})}}}sin\theta d\theta d\phi \}.\end{array}$$

(34)

Note that the difference between Equation (34) and Equation (28) is that the complex wavenumber $k\left(\mathit{p}\right)$ in Equation (28) is replaced by the unperturbed wavenumber $k_M\left(\mathit{p}\right)$.

Based on the Born approximation, one can retrieve the attenuation coefficient from the imaginary part of the perturbed wavenumber by Equations (33) and (34), and the result is:

$$\begin{array}{l}{\alpha }_M(\mathbf{p})=-\mathrm{Im}\left[\frac{m_M(k_M)}{2k_M^{}(\mathbf{p})V_M^2(\mathbf{p})}\right]=\\ {\displaystyle \sum _{N=1}^3\frac{{\pi \omega }^4}{4{\rho }^2{V}_M^3(\mathbf{p})}{\displaystyle {\int }_0^{\pi }{\displaystyle {\int }_0^{2\pi }\frac{I{P}_{M\to N}^{}(\mathbf{p},\mathbf{s})W(k_M-k_N^S)}{V_N^5(\mathbf{s})}}}sin\theta d\theta d\phi },\end{array}$$

(35)

which has been reported in Ref. [40].

Now we examine the real part of the perturbed wavenumber that governs the phase velocity in textured polycrystalline materials, and it has never been studied in previous studies [6,39,40,59]. From Equations (33) and (34), the real part of the complex wavenumber $k\left(\mathit{p}\right)$ is explicitly rewritten as:

$$\begin{array}{l}\mathrm{Re}[k(\mathbf{p})]=k_M^{}(\mathbf{p})-\mathrm{Re}\left[\frac{m_M(k_M)}{2k_M^{}(\mathbf{p})V_M^2(\mathbf{p})}\right]=k_M^{}(\mathbf{p})-\\ {\displaystyle \sum _{N=1}^3\frac{k_M(\mathbf{p})}{2{\rho }^2{V}_M^2(\mathbf{p})}P.V.{\displaystyle {\int }_0^{\infty }\frac{k_S^4}{k_N^2(\mathbf{s})-k_S^2}}d{k}_S{\displaystyle {\int }_0^{\pi }{\displaystyle {\int }_0^{2\pi }\frac{I{P}_{M\to N}^{}(\mathbf{p},\mathbf{s})W(k_M-k_{}^S)}{V_N^2(\mathbf{s})}sin\theta d\theta d\phi }}}.\end{array}$$

(36)

From Equation (36), the frequency-dependent phase velocity $V\left(\mathit{p}\right)$ can be finally obtained as:

$$V(\mathbf{p})=V_M^{}(\mathbf{p}){\left\{1-{\displaystyle \sum _{N=1}^3\frac{1}{2{\rho }^2{V}_M^2(\mathbf{p})}P.V.{\displaystyle {\int }_0^{\infty }\frac{k_S^4}{k_N^2(\mathbf{s})-k_S^2}}d{k}_S{\displaystyle {\int }_0^{\pi }{\displaystyle {\int }_0^{2\pi }\frac{I{P}_{M\to N}^{}(\mathbf{p},\mathbf{s})W(k_M-k_{}^S)}{V_N^2(\mathbf{s})}sin\theta d\theta d\phi }}}\right\}}^{-1}.$$

(37)

Explicit expressions for the attenuation and phase velocity thus have been obtained in this section. However, due to the intrinsic limitation of the Born approximation [66], the explicit expressions, Equations (35) and (37) for attenuation and phase velocity solely validate for the frequency below the geometric region.

3.2.2. Attenuation at Rayleigh Limit

The attenuation and phase velocity obtained from the Born approximation also can predict the asymptotic behavior of the attenuation and phase velocity. At Rayleigh limit, the wavenumber k_M approaches zero and $W\left(\mathit{k}-{\mathit{k}}_{}^S\right)=\frac{R_0{R}_1{a}_X^3}{{\pi }^2}$, thus the attenuation coefficient at Rayleigh limit ${\alpha }_M^R\left(\mathit{p}\right)$ by Equation (35) becomes:

$${\alpha }_M^R(\mathbf{p})={\displaystyle \sum _{N=1}^3\frac{{\omega }^{4}V}{32{\pi }^2{\rho }^2{V}_M^3(\mathbf{p})}{\displaystyle {\int }_0^{\pi }{\displaystyle {\int }_0^{2\pi }\frac{I{P}_{M\to N}^{}(\mathbf{p},\mathbf{s})}{V_N^5(\mathbf{s})}}}sin\theta d\theta d\phi },$$

(38)

where $V=8\pi a_X{a}_Y{a}_Z$ is the effective grain volume. Thus, the attenuation is proportional to the grain volume and fourth power of frequency at the Rayleigh region, and Equation (38) is consistent with the Rayleigh limit derived in Ref. [40].

3.2.3. Quasi-Static Phase Velocity

One remarkable aspect of the asymptotic complex wavenumber at Rayleigh limit is that it can provide quasi-static velocity for textured polycrystalline materials at an arbitrary wave propagation direction, from which one can retrieve the effective elastic constants of textured polycrystals (this inverse problem will not be addressed here). Homogenization of polycrystals is of great importance for predicting the mechanical properties of polycrystals and has drawn the attention of many researchers [58,61,67,68,69]. For statistically isotropic polycrystals, different velocity bounds, including Voigt–Reuss and Hashin–Shtrikman bounds [67], and a rather accurate prediction named self-consistent approach [61] have been developed to predict the effective elastic constants of grains ensemble. For textured polycrystals, Voigt–Reuss bounds [70,71] Hashin–Shtrikman bounds [72,73,74] and self-consistent approach [73,75] have also been reported, from which one can calculate the phase velocity of an incident elastic wave. Although the self-consistent approach in Refs. [73,75] is still limited to high-symmetry crystallites and certain texture symmetry, a recent study [74] developed Hashin–Shtrikman bounds for polycrystals with arbitrary crystallographic texture and triclinic crystallites. The quasi-static velocities also can be obtained in this work for aggregates of triclinic grains with arbitrary texture symmetry and it counts grain shape impact as well. The detailed derivation is given below.

When k_M reaches zero, the Cauchy principal value can be analytically obtained by contour integral. After much simplification, the real part of the perturbed wavenumber at Rayleigh limit is obtained as:

$$\begin{array}{l}\mathrm{Re}[k(\mathbf{p})]=k_M^{}(\mathbf{p})+\\ {\displaystyle \sum _{N=1}^3\frac{k_M^{}(\mathbf{p})R_0{R}_1}{8\pi {\rho }^2{V}_M^2(\mathbf{p})}}{\displaystyle {\int }_0^{\pi }{\displaystyle {\int }_0^{2\pi }\frac{I{P}_{M\to N}^{}(\mathbf{p},\mathbf{s})}{V_N^2(\mathbf{s}){[1+(R_0^2-1){\cos }^2\theta +(R_1^2-1){\sin }^2\theta {\sin }^2\varphi ]}^{3/2}}\sin \theta d\theta d\varphi }}.\end{array}$$

(39)

From the equation above, the quasi-static velocity at a certain wave propagation direction p also can be obtained as:

$$\begin{array}{l}V(\mathbf{p})=V_M^{}(\mathbf{p})\{1+\\ {{\displaystyle \sum _{N=1}^3\frac{R_0{R}_1}{8\pi {\rho }^2{V}_M^2(\mathbf{p})}}{\displaystyle {\int }_0^{\pi }{\displaystyle {\int }_0^{2\pi }\frac{I{P}_{M\to N}^{}(\mathbf{p},\mathbf{s})}{V_N^2(\mathbf{s}){[1+(R_0^2-1){\cos }^2\theta +(R_1^2-1){\sin }^2\theta {\sin }^2\varphi ]}^{3/2}}\sin \theta d\theta d\varphi }}\}}^{-1}.\end{array}$$

(40)

It indicates that the quasi-static velocity is determined collectively by grain shape (related to R₀ and R₁), texture (affecting wave velocity and inner product) and wave scattering. For a polycrystal with known texture coefficient, single-crystal elastic constants, and grain sizes, the quasi-static velocity of an incident wave (QL, QT1 or QT2) at a certain direction can be predicted by Equation (40) and the elastic constants of the homogenized polycrystal also can be obtained inversely if sufficient phase velocities at different angles are available. For textured polycrystals with equiaxed grains, the expression for the quasi-static velocity can be further simplified as:

$$V(\mathbf{p})=V_M^{}(\mathbf{p}){\left\{1+{\displaystyle \sum _{N=1}^3\frac{1}{8\pi {\rho }^2{V}_M^2(\mathbf{p})}}{\displaystyle {\int }_0^{\pi }{\displaystyle {\int }_0^{2\pi }\frac{I{P}_{M\to N}^{}(\mathbf{p},\mathbf{s})}{V_N^2(\mathbf{s})}\sin \theta d\theta d\varphi }}\right\}}^{-1}.$$

(41)

From Equation (41), the effective phase velocity is governed by texture and wave scattering. A similar homogenization method based on scattering theory for texture-free polycrystalline materials has been reported in Ref. [76], however, this work not only accounts for crystallographic texture impact on phase velocity but also provides explicit expressions for phase velocities (limited to second-order scattering). Therefore, Equation (40) provides a general formula for quasi-static velocities in polycrystals, from which one can estimate the effective elastic constants of grains ensemble.

4. Computational Results and Discussion

The focus of this section is on the computational results from the SOA model and comparisons to other models and available experimental results in the literature. First, the SOA model is validated by comparison to reported experimental attenuation data and finite element modeling results for a 316L stainless steel sample in references [22,77]. Furthermore, the applicability and limitation of the Born approximation are investigated and the texture impacts on attenuation and velocity are studied. Finally, the quasi-static phase velocities obtained in this work by Equation (40) are compared with available bounds on a textured polycrystalline copper.

4.1. Comparison of the SOA Model to Available FEM and Experimental Results

Since two recent studies [22,77] reported experimental attenuation results and FEM results on textured stainless steel and nickel alloy, it will be interesting to see the comparison between the SOA model in this work and references [22,77] considering these three are independent approaches. Here a textured stainless steel polycrystal with perfect transverse isotropy texture, mean grain diameter 0.25 mm and mean grain length 5 mm is chosen because its attenuation was experimental measured in Ref. [77] and numerically modeled by FEM in Ref. [22]. The single-crystal elastic constants, texture coefficients, and material density are also given in references [22,77]. Thus, attenuation coefficients can be calculated using the FORTRAN code by the SOA model for a quantitative comparison with references [22,77].

The comparison of quasi-longitudinal attenuation coefficients between the SOA model in this work and Refs. [22,77] for stainless steel is shown in Figure 4, where the data points for experiment and FEM are directly from Ref. [22,77] by an image digitizer. Unlike Refs. [22,77], where wave propagation direction is fixed while the grain orientation is treated as a variable, this work assumes that grain elongation direction is fixed (namely aligned with Z-direction) but the wave propagation direction τ changes from 0° to 90°. Obviously, these two computation methods are equivalent. From Figure 4, one can see that attenuation results from this work have a reasonable agreement with FEM [22] since the attenuation coefficients from this work fall within the error bar ranges of FEM. However, the SOA model somehow underestimates the attenuation coefficient (the predicted attenuation form the SOA model is lower than the corresponding mean attenuation of FEM at different angles). Besides, the SOA model agrees well with experiment measurements [77] except for angles 0°, 10° and 45°. Accurate attenuation measurement is difficult when the attenuation is too small (at angles 0 and 10 degrees) because the real attenuation is comparable to the measurement system error [78]. Normally, the experiment may deviate from the models by a fact of two [79] even when the attenuation is reasonably high. Therefore, a good agreement between the SOA model and the experiment is justified. Compared with experimental results, the SOA model also underestimates the attenuation. Such a systematic discrepancy may be caused by a finite number of grains in the stainless steel samples, since a recent study [42] reveals that grain number can impose an obvious impact on attenuation coefficient and attenuation shows more fluctuation on a finite number of grains. In other words, the FEM modeling [22] and experiment measurements [77] include grain number effect on attenuation coefficients while the SOA model in this work excludes this effect by assuming infinite grains. Another reason may be the assumption using perfect ellipsoidal grains in the SOA model to represent columnar 316L grains in Refs. [22,77]. Further investigation would require additional modeling and it is out of the scope of this work.

4.2. Applicability and Limitation of the Born Approximation on Textured Polycrystals

In this section, a triclinic Ti polycrystal from Refs. [40,50] is employed as an example to calculate attenuation and phase velocity via the SOA model. The single-crystal elastic constants and density of this triclinic Ti crystal from Refs. [40,50] are listed in

Table 1. The elastic constants $c_{ij}$ (in GPa) and the mass density ρ(g/cm³) of a triclinic Ti crystallite reported in Ref. [40,50].

	$c_{11}$	$c_{12}$	$c_{13}$	$c_{14}$	$c_{15}$	$c_{16}$	$c_{22}$	$c_{23}$	$c_{24}$
Triclinic Ti	161.6	67.8	88.2	1.8	1.1	−0.64	184.4	68.0	−2.5
	$c_{25}$	$c_{26}$	$c_{33}$	$c_{34}$	$c_{35}$	$c_{36}$	$c_{44}$	$c_{45}$	$c_{46}$
Triclinic Ti	−1.9	−2.6	161.3	−0.18	0.83	2.4	50.9	−5.1	0.20
	$c_{55}$	$c_{56}$	$c_{66}$	ρ
Triclinic Ti	39.0	−5.2	50.9	4.54

. Assume that the polycrystal is composed of triclinic Ti crystals of ellipsoidal shape, whose main axes coincide with the global coordinate axes as shown in Figure 3. The same example in Ref. [40] is used here to facilitate the comparison between the SOA model in this work and the Born approximation in Ref. [40]: the mean grain radii are $a_X=0.1\text{mm}$, $a_Y=0.3 \text{mm}$ and $a_Z=0.5 \text{mm}$ and the macrotexture parameters are ${\sigma }_{\theta }=0.03,{\sigma }_{\varphi }={\sigma }_{\theta }$ and ${\sigma }_{\zeta }=3{\sigma }_{\theta }$. In the SOA model, the calculation procedures of effective elastic constants and elastic constants covariances using 3D Gaussian ODF are the same as references [40,50], thus they will not be repeated here. The mass operator, Cauchy principal value, and dispersion equation are solved numerically by FORTRAN code, and only three wave propagation directions (X, Y and Z) are considered in this work. The computed attenuation results from the SOA model in this work will be compared with those from Ref. [40] by the Born approximation.

Figure 5 shows the longitudinal attenuation comparison between the SOA model and the Born approximation [40] for the Ti polycrystal with ellipsoidal grains ($a_X=0.1 \text{mm},a_Y=0.3 \text{mm}$ and $a_Z=0.5 \text{mm}$). The attenuation and wavenumber are normalized by a_X so that the horizontal and vertical axes are dimensionless. It can be seen from Figure 5 that the Born approximation is in reasonable agreement with our SOA model below the geometric region (also refer to the mean deviation in

Table 2. The mean deviation between the SOA model and the Born approximation below break frequency (SOA model is treated as the reference) and the break frequency $k_L{a}_X$ of the Born approximation for a textured polycrystal with ellipsoidal Ti grains.

Direction	X	Y	Z
Mean deviation	−13%	−9.7%	−14%
Departure frequency	31.6	12.7	25.1

) and the transition frequency to geometric region varies in different directions. The deviation is defined as $\left({\alpha }_{Born}-{\alpha }_{SOA}\right)/{\alpha }_{SOA}$. The signs for all mean deviations at different wave propagation directions in Table 2 are negative, which indicates that the Born approximation underestimates the wave attenuation because it accounts for fewer scattering events than the SOA model. In the geometric region, the Born approximation [40] breaks due to its intrinsic limitation, and the break frequency $k_L{a}_X$ of the Born approximation for each direction is also listed in Table 2. Since the departure frequency $k_L{a}_X$ of the Born approximation from the SOA model is comparable to the stochastic-to-geometric transition frequency in the SOA model (sharp corner shown in Figure 5 from our SOA model), this departure frequency is also governed by the grain elongation and longitudinal wave scattering strength at high frequency. The attenuation at the geometric region from this work for textured polycrystals is also inversely proportional to grain size in the wave propagation direction as texture-free polycrystal case in Ref. [64,80], thus attenuation coefficients in geometric region predicated by the SOA model obey this relation: ${\alpha }_X>{\alpha }_Y>{\alpha }_Z.$ In the geometric region, the attenuation is so high that the penetration depth of an elastic wave is limited to a couple of grains [81].

Figure 6 indicates the longitudinal attenuation comparison between SOA model and the Born approximation [40] for a Ti polycrystal with equiaxed grains ($a_X=a_Y=a_Z=0.1 \text{mm}$), where the texture parameters are the same as Figure 5. From Figure 6, it can be seen that the Born approximation still reasonably agrees with the SOA model except for the geometric region. As shown in

Table 3. The mean deviation between SOA and the Born approximation below break frequency (SOA model is treated as reference) and the break frequency $k_L{a}_X$ of the Born approximation for a textured polycrystal with equiaxed Ti grains.

Direction	X	Y	Z
Mean deviation	−10.9%	−7.5%	−14.4%
Departure frequency	31.0	34.4	102.7

, the mean deviations (below departure frequency) of the Born approximation from the SOA model are −10.9% for X direction, −7.5% for Y direction and −14.4% for the Z direction, respectively. The model in Ref. [40] underestimates the attenuation coefficient on this case as well. The break frequency $k_L{a}_X$ of the Born approximation [40] is also direction-dependent (see Table 3), and according to Table 3 the crystallographic texture affects the departure frequency most in the Z direction. In the geometric region, attenuation coefficients for X, Y and Z directions are comparable due to equiaxed grains.

The comparisons between the SOA model and the Born approximation [40] on these two cases (shown in Figure 5 and Figure 6, respectively) indicate that the Born approximation may underestimate the attenuation coefficient up to 15% below the geometric region. However, this discrepancy is still acceptable, considering that experimental attenuation can deviate from the theoretical prediction by a factor of two [82]. As for the break frequency of the Born approximation, it is affected collectively by texture extent, wave scattering extent, and the grain shape.

The phase velocities from the SOA model in this work for these two cases (textured Ti polycrystal with equiaxed grains and textured Ti polycrystal with elongated grains) are plotted in Figure 7 to demonstrate the dispersion behaviors of phase velocity at different frequency regimes, where normalized phase velocities (deviation from Voigt velocity) are shown. The phase velocity results in different directions for the Ti polycrystal with ellipsoidal grains are shown in Figure 7a, while the phase velocity results for the equiaxed grain case are depicted in Figure 7b. In the Rayleigh region, the phase velocities (quasi-static velocities) at X, Y and Z directions are non-dispersive for both cases and velocities in different directions follow this relation: Y > X > Z because of texture impact. Note that quasi-static phase velocities from these two cases for X, Y and Z directions are slightly different (less than 0.06%), which manifests the grain shape impact on quasi-static phase velocity. From Figure 7, the phase velocities for both cases increase from the Rayleigh region to the stochastic region and they reach a plateau before they become highly dispersive at the beginning of the geometric region. Although two branches should exist in wave propagation [32], the dominant branch is exclusively shown in Figure 7 because the subdominant branch is even harder to compute for textured polycrystals (like Ref. [36], a special approach has to be used for calculating the subdominant branch). In Figure 7a the phase velocity from X direction is above its Voigt velocity at geometric region while these for Y and Z directions in geometric region are below their corresponding Voigt velocities; intriguingly, in Figure 7b the phase velocity in the Z direction at geometric region is above its Voigt velocity while the rest are below their corresponding Voigt velocities. This interesting observation was also reported on texture-free polycrystals [3,34], but the physical meaning has not been well-interpreted. Literally, the wave chooses to propagate along faster phase velocity direction of grains at some scenarios but slower velocity direction at other scenarios. A highly dispersive phase velocity at the beginning of the geometric region also means that the energy propagation direction deviates from the propagation direction of the wavefront. Further investigation on the phase velocity in the geometric region is beyond the scope of this work.

From the phase velocity results, we can conclude that the texture impact on phase velocity is obvious in the Rayleigh region (the quasi-static phase velocities are changed) and geometric region (at some direction the phase velocity is larger than the Voigt velocity while at another direction the phase velocity becomes lower than its corresponding Voigt velocity). Additionally, the texture impact on phase velocity decreases with frequency in the transition region (from the Rayleigh region to the stochastic region) and its influence is minimal at the stochastic region. The texture also affects the transition frequency to the geometric region significantly. Last but not least, a competing effect between texture and grain shape is seen on phase velocity for textured polycrystals with elongated grain, but the grain shape impact is relatively weaker compared with crystallographic texture.

4.3. Comparison of Quasi-Static Velocity from the SOA Model with Other Velocity Bounds

An example is given here to validate the quasi-static velocity in Equation (40). The Hashin–Shtrikman bounds for anisotropic aggregates of cubic grains with arbitrary texture symmetry have been derived in Ref. [73] using variational principals while a recent study [74] explicitly represented the Hashin–Shtrikman bounds by second and fourth-order ODCs for polycrystals with arbitrary texture symmetry and crystal symmetry. Since the elastic modulus bounds of a textured polycrystalline copper (orthotropic texture symmetry) with equiaxed grains are available in Refs. [73,74], it makes possible a quantitative comparison between this work and Refs. [73,74]. Note that the derivation of Hashin–Shtrikman bounds in Ref. [74] is more rigorous than Ref. [73] because of a more generous assumption for stress polarization in Ref. [74] and the Hashin–Shtrikman bounds in Ref. [74] apply to more general texture and crystallite system cases; however, these two methods in Refs. [73,74] are comparable as demonstrated by Ref. [74]. For this textured polycrystalline copper, the orientation distribution coefficients (ODCs) of the Wigner-D function ODF from Ref. [73] are $c_{00}^4=-0.0055154$, $c_{20}^4=0.006364$, $c_{40}^4=-0.0165463$ and all other $c_{m0}^4=0$ (m stands for an integer). These ODCs in Wigner-D function ODF are converted to ODCs in Roe’s notation before they are applied to SOA model and the conversion formula is $W_{lmn}={(-1)}^{m-n}\sqrt{\frac{2}{2l+1}}{c}_{mn}^l$ ($W_{lmn}$ are ODCs in Roe’s notation). Single-crystal elastic constants of the copper grain are also from Ref. [73]. The Voigt and Reuss bounds for this textured copper reported in Refs. [73,74] were first repeated using Roe’s ODF as a cross-check of author’s FORTRAN code before calculating the quasi-static phase velocities.

The quasi-static phase velocities for this textured copper at three different wave propagation directions (X, Y and Z) are obtained based on Equation (40), and compared with velocities calculated from corresponding Hashin–Shtrikman (H-S) and self-consistent (SC) effective moduli in Refs. [73,74]. The comparison is shown in

Table 4. Comparison of the longitudinal velocities (in km/s) at different propagation directions for a macroscopically orthotropic copper with equiaxed grains [73,74]. Hashin–Shtrikman (H-S) bounds and self-consistent (SC) velocity are directly from Ref. [73,74]. Assume copper density is 8.935g/cm³.

Directions	Reuss	H-S Lower [73]	H-S Lower [74]	SC [73]	This work, Equation (40)	H-S Upper [74]	H-S Upper [73]	Voigt
X	5.038	4.919	4.935	4.969	4.979	4.983	4.987	4.841
Y	4.965	4.834	4.842	4.883	4.892	4.897	4.900	4.743
Z	4.902	4.757	4.759	4.805	4.816	4.820	4.823	4.656

, where the Voigt and Reuss velocities are also listed. A reasonable homogenization model should follow Hashin–Shtrikman bounds and provide results close to the SC approach. From Table 4, the quasi-static velocities from this work are much better than Voigt–Reuss bounds and they obey the H-S upper and lower bounds provided in Refs. [73,74] (H-S upper and lower bounds from Refs. [73,74] are comparable), and they are also reasonably close to SC velocity given in Ref. [73]. Therefore, the explicit expression Equation (40) for quasi-static velocities obtained in this work is justified.

Although sophisticated bounds including Hashin–Shtrikman and self-consistent approaches have been developed for polycrystalline aggregates with different texture symmetry and grain symmetry like Refs. [72,74,75], no further attempt is made in this paper to systematically compare with these studies [72,74,75] because of certain constraints: these models in Refs. [72,74,75] are not straightforward to apply; no additional computational example is available in Refs. [72,74]; ODCs are not provided in Ref. [75]. More will be exploited on this aspect in further work.

To inversely determine elastic moduli from phase velocity, one has to collect sufficient phase velocity data at multiple angles and fit the elastic moduli by resorting to Christoffel’s equation. This will not be demonstrated here since details about the determination of effective elastic moduli from ultrasonic phase velocities can be found in Ref. [15].

5. Conclusions

A second-order attenuation model for textured polycrystals, which is applicable to triclinic crystallites of ellipsoidal shape and valid for the whole frequency range, has been derived in terms of the Dyson equation in this study. The SOA model can predict both attenuation and phase velocity including the quasi-static velocity. Moreover, the quasi-longitudinal attenuation coefficient predicted by the SOA model is compared to available FEM results and experimental attenuation data on a textured stainless steel polycrystal. The SOA model is further compared with an existing model [40] based on the Born approximation on two textured polycrystals with triclinic titanium grains but different grain morphology. The phase velocity results predicted by the SOA model for these two textured Ti polycrystals are also presented to demonstrate the texture impact on phase velocity. Finally, quasi-static velocities from this work on a textured polycrystalline copper are compared with other velocity bounds reported in the references [73,74] including self-consistent velocity. Several conclusions can be drawn in this study.

The attenuation comparison of the SOA model to available FEM results [22] and experimental attenuation data [77] on the textured stainless steel polycrystal indicates the SOA model is in reasonable agreement with FEM and experimental measurement. This justifies the SOA model developed in this work.

From comparisons between the SOA model and the existing model [40] based on the Born approximation on two textured polycrystalline aggregates of triclinic Ti grains (one ellipsoidal grain case and one equiaxed grain case), the SOA model validates in the whole frequency range and quantitatively shows the applicable frequency range of the Born approximation [40]. It is also concluded the accuracy of the existing model in Ref. [40] is still acceptable below the geometric region.

The computational results for these two textured Ti polycrystals show that the phase velocity predicted by the SOA model in the whole frequency range is dispersive (frequency-dependent), which is affected both by crystallographic texture and the grain shape. The texture impact on phase velocity is most phenomenal in the Rayleigh region and geometric region. Therefore, one big advantage of the SOA model is that it simultaneously provides attenuation and phase velocity.

The quasi-static velocity obtained from the SOA model applies to polycrystals with arbitrary texture symmetry and triclinic grains and also accounts for different grain shapes, which had never been achieved before. The quasi-static velocities from the SOA model obey the Hashin–Shtrikman bounds and agree reasonably with self-consistent velocity. Therefore, it may provide another homogenization approach to determine effective elastic constants of texture polycrystals.

This work can be further extended to polycrystals with a generalized two-point correlation function. Recent studies [23,83] have shown that the GTPC function of a Voronoi grain microstructure may not follow an exponential two-point correlation function. However, the numerically determined two-point correlation function can be incorporated into this work by following the approach in Ref. [23] using a summation of multiple exponential terms. Further development will be addressed in the future.