Three-dimensional finite element (FE) modelling, with representation of materials at grain scale in realistic sample volumes, is capable of accurately describing elastic wave propagation and scattering within polycrystals. A broader and better future use of this FE method requires several important topics to be fully understood, and this work presents studies addressing this aim. The first topic concerns the determination of effective media parameters, namely, scattering induced attenuation and phase velocity, from measured coherent waves. This work evaluates two determination approaches, through-transmission and fitting, and it is found that these approaches are practically equivalent and can thus be used interchangeably. For the second topic of estimating modelling errors and uncertainties, this work performs thorough analytical and numerical studies to estimate those caused by both FE approximations and statistical considerations. It is demonstrated that the errors and uncertainties can be well suppressed by using a proper combination of modelling parameters. For the last topic of incorporating FE model information into theoretical models, this work presents elaborated investigations and shows that to improve agreement between the FE and theoretical models, the symmetry boundary conditions used in FE models need to be considered in the two-point correlation function, which is required by theoretical models.

## I. INTRODUCTION

A polycrystalline material is composed of elastically anisotropic grains with differently oriented crystallographic axes. The elastic properties within such a medium are thus spatially varied, and elastic waves are scattered as they propagate, leading to the waves being attenuated and their phase velocities being dispersive. These wave behaviors are undesirable to the non-destructive testing for the detection of defects, a major reason being that attenuation can weaken the detected signal that may be masked by the grain scattering noise.^{1} On the other hand, the wave behaviors can also be used advantageously for the characterization of microstructure because both attenuation and dispersion carry bulk information about grain structures.^{2,3} Therefore, understanding these behaviors is practically important.

Extensive studies have been carried out to understand the behaviors of elastic waves in polycrystals. Early works, a review of which can be found in Papadakis,^{4} focused on experimental investigations and the theoretical explanations of experimental results. These early efforts led to a variety of theoretical models that showed good qualitative agreement with experiments,^{4} but each of them applies only to a specific scattering regime. A unified model that is valid across all scattering regimes was first presented by Stanke and Kino^{5} by applying the Keller approximation^{6,7} to polycrystalline media. An equivalent model was later developed by Weaver^{8} based on the Dyson equation with the introduction of the first-order smoothing approximation.^{9} Further theoretical developments^{2,10–19} were mostly extensions of these two models, the Weaver model in particular due to its ease of extension, to various polycrystals of different elastic properties and grain geometries.

Theoretical models are necessarily approximate because it is impossible to carry a rigorous mathematical description through all steps of the model development. Therefore, assumptions and approximations are made in the derivations, leading to expected, but mostly not quantified, approximations in the outputs. The Stanke and Kino^{5} and Weaver^{8} models are accurate up to the second order of material inhomogeneities, offering them possibilities of considering forward multiple scattering. However, except the rare cases^{11,14,17–20} that maintained the second-order accuracy, other developments^{2,10,12,13,15,16} of the Weaver model (including the eventual formulation of Weaver^{8}) were formulated by invoking the Born approximation to facilitate explicit calculations for attenuation. Such treatment limits the validities of the models to even lower-order scattering in frequency ranges below the geometric regime.^{8}

In contrast to theoretical studies, numerical simulations can be performed without the approximations of low order scattering. A wide range of numerical methods can be employed to achieve such numerical experiments and a summary of them can be found in Van Pamel *et al.*^{21} Among these methods, the finite element (FE) method is widely researched for the simulation of elastic waves in complex media.^{22} The use of this method for polycrystals has long been limited by its extreme computational demands to two-dimensional (2D) cases.^{23–28} However, the dependence of scattering on spatial dimensionality^{29} at certain frequencies means that three-dimensional (3D) simulation is always desirable for understanding real-world wave phenomena. Recent advancements of computer hardware and elastodynamic software, especially the development of the GPU-accelerated Pogo program,^{30} have enabled the successes of 3D modelling.^{11,19,21,31,32} These 3D simulation studies, with representations of polycrystals at grain scale in realistic sample volumes, addressed a relatively simple case of plane longitudinal waves propagating in polycrystals with macroscopically isotropic elastic properties and equiaxed grains. They have demonstrated that 3D FE simulations can accurately describe the interactions of elastic waves with grains that involve complex multiple scattering.

Due to the flexibility and accuracy of 3D FE modelling, we can expect that there will be a surge of modelling studies in the future to understand and optimize practical inspection problems, and to evaluate the approximations of theoretical models for model-based applications, such as the inverse characterization of microstructure. However, even for the relatively simple cases addressed so far,^{11,19,21,31,32} some essential topics are not yet understood well enough to allow confidence that the best accuracy of simulation is being achieved. The authors have been investigating this challenge, and the purpose of this article is to report on important topics of modelling methodology for which new knowledge will inform a reliable quality of simulation performance.

The first topic is the determination of effective media parameters from the simulation results, namely, attenuation and phase velocity. An important task of these addressed to 3D FE simulations is to understand how a transmitted plane wave changes with propagation distance. This change is parametrically described by attenuation coefficient and phase velocity, representing, respectively, the amplitude and phase parts of the change. One can perform the determination in a through-transmission (TT) configuration, which is commonly used in actual experiments where access is limited to sample surfaces. This method has seen wide applications in our previous simulations,^{11,19,21,29,31} in which the effective media parameters are determined by comparing the two waves measured on the transmitting surface and its opposite receiving surface. Alternatively, one can also determine the parameters by best fitting the waves acquired on multiple surfaces within the simulation volume that are parallel to the transmitting surface. This fitting method was first used in a recent work^{32} and it is made possible by the intrinsic advantage of FE modelling that waves inside a sample can be easily measured. Each of these two methods has its own advantages: the TT method is consistent with physical experimental configurations, so it is conceptually attractive, and it requires relatively little measured data, while the fitting method can potentially provide more abundant information about wave-microstructure interaction. However, these two methods have not been compared in actual simulations to identify their ranges of application and their possible advantages and drawbacks. It is important to recognize that while the difference between these methods would be trivial for a homogeneous material, this is not the case for a polycrystal: the grain scattering creates high levels of noise in the measurements of the coherent signals as well as local influences of multiple scattering, both of which bring significant uncertainty to the measurements.

The second topic concerns the errors and uncertainties of the determined scattering parameters. It is known that numerical approximations, including the space and time discretization in the FE method and the inexactness in finite-precision computations, can cause numerical errors. In addition, statistical considerations, involving the random distribution of grain structures and the stochastic pattern of wave scattering, can induce both statistical errors and uncertainties. For numerical errors, studies so far were partly limited to homogeneous and isotropic materials^{33,34} and were partly qualitative when they went to polycrystalline materials,^{21,31} whereas for statistical errors and uncertainties, attention has not yet been turned to the specific problem of wave propagation within polycrystals.

The third of our three topics of study involves the evaluation of the performance of theoretical models using the simulation data of 3D FE modelling. The simulations are ideal for the evaluation of the effects of the approximations of theoretical models, as they can be fully controlled, and their errors and uncertainties can to some extent be estimated. The evaluation requires the statistical information of FE material systems be measured and put into theoretical models. Such information is represented by the two-point correlation (TPC) function for widely used theoretical models, like the second-order^{5,8,11,19} and far-field^{20} approximations. Previous studies^{11,19,32} have been successful in incorporating the TPC function of polycrystal models into theoretical models, and have shown good agreement between FE and theoretical results in terms of scattering parameters. However, these studies neglected in their TPC functions the fact that symmetry boundary conditions (SBCs), employed to emulate infinitely wide media that can hold plane waves, actually mirror the boundary grains of FE models and thus make the grains doubled in volume.

For the benefit of future applications of the 3D FE method, this work presents a thorough discussion of the above-mentioned topics, with the aim to solve them appropriately. The discussion focuses on the same simulation case addressed in currently-available works,^{11,19,21,31,32} which is the propagation of plane longitudinal waves in polycrystals with untextured properties and equiaxed grains. The goal of this work, however, is to lay a solid foundation for general cases of arbitrary wave modalities in any polycrystalline material system. For the specific case concerned in this work, its 3D FE model is summarized in Sec. II. The subsequent three sections discuss, respectively, the three topics mentioned above. Section III deals with the determination of effective media parameters. Section IV presents a systematic study of errors and uncertainties, with discussions on how to suppress them. Section V provides a comprehensive discussion covering the incorporation of 3D FE simulation data into theoretical models and the proper consideration of SBCs in this incorporation procedure.

## II. FE METHOD FOR ELASTIC WAVES IN POLYCRYSTALS

The three studies reported here were conducted using a 3D FE model of a volume of a polycrystal, simulating the propagation of plane longitudinal waves in the time domain. These kinds of models and studies using it have been reported previously,^{21,31,32} but for the convenience of subsequent discussions, we provide in the following a summary of the method.

### A. Representation of polycrystals at grain scale

Scattering only occurs when grains are present, so the representation of polycrystals at grain scale is the key to reproducing scattering. Researchers have simulated the spatial properties of polycrystals with several methods, two popular ones being Voronoi tessellation (e.g., Quey *et al.*^{35}) and Dream.3D.^{36} In this work, we have used Voronoi tessellation.^{35}

The Voronoi tessellation partitions a cuboid model, which is shown in Fig. 1, into convex polyhedra or grains that have no overlaps and no gaps. Each grain encloses a seed that is specified beforehand and all points in this grain are closer to its own seed than to any other. Grain statistics, like grain size distribution and average grain shape, are determined by how the seeds are specified. The present work uses a relatively simple seeding procedure that utilizes a Poisson point process to drop random seeds into the model space. As would be expected from the central limit theorem, the grain sizes (specifically the cubic roots of grain volumes) of the generated model are normally distributed^{21} and the grain shapes are equiaxed on average. For later discussions, three such models are generated, and their parameters are provided in Table I.

Model . | $dx\xd7dy\xd7dz$ . | $N$ . | $a$ . | $h$ . | d.o.f. . | $fc$ . | |
---|---|---|---|---|---|---|---|

Aluminium . | Inconel . | ||||||

N115200 | 12 × 12×100 | 115 200 | 0.5 | 0.050 | 349 × 10^{6} | 2, 5 | 2 |

N11520 | 12 × 12×10 | 11 520 | 0.5 | 0.025 | 278 × 10^{6} | 10 | 5 |

N16000 | 20 × 20×5 | 16 000 | 0.5 | 0.020 | 755 × 10^{6} | 20 | 10, 20 |

Model . | $dx\xd7dy\xd7dz$ . | $N$ . | $a$ . | $h$ . | d.o.f. . | $fc$ . | |
---|---|---|---|---|---|---|---|

Aluminium . | Inconel . | ||||||

N115200 | 12 × 12×100 | 115 200 | 0.5 | 0.050 | 349 × 10^{6} | 2, 5 | 2 |

N11520 | 12 × 12×10 | 11 520 | 0.5 | 0.025 | 278 × 10^{6} | 10 | 5 |

N16000 | 20 × 20×5 | 16 000 | 0.5 | 0.020 | 755 × 10^{6} | 20 | 10, 20 |

Material properties are then assigned to the generated models. This study considers single-phase polycrystals, and therefore the grains of each model belong to the same crystal material. The two materials considered in this study are given in Table II, both having cubic symmetry and being chosen as examples of weakly and strongly scattering materials, respectively. Although in each model all grains share an identical material, their crystallographic axes are randomly oriented. Therefore, the models are homogeneous and statically isotropic on the macroscopic scale.

Material . | $c11$ . | $c12$ . | $c44$ . | $\rho $ . | $A$ . | $VL$ . | $VT$ . | $QL\u2192L$ . | $QL\u2192T$ . |
---|---|---|---|---|---|---|---|---|---|

aluminum | 103.40 | 57.10 | 28.60 | 2700 | 1.24 | 6317.52 | 3128.13 | 0.08 | 0.33 |

Inconel | 234.60 | 145.40 | 126.20 | 8260 | 2.83 | 6025.37 | 3365.54 | 2.26 | 7.59 |

Material . | $c11$ . | $c12$ . | $c44$ . | $\rho $ . | $A$ . | $VL$ . | $VT$ . | $QL\u2192L$ . | $QL\u2192T$ . |
---|---|---|---|---|---|---|---|---|---|

aluminum | 103.40 | 57.10 | 28.60 | 2700 | 1.24 | 6317.52 | 3128.13 | 0.08 | 0.33 |

Inconel | 234.60 | 145.40 | 126.20 | 8260 | 2.83 | 6025.37 | 3365.54 | 2.26 | 7.59 |

It is important to realize that every single model in Table I has a considerable sample volume and a large number of grains, but still, it is only an incomplete random realization (sample) of the parent polycrystal that one aims to simulate. Using a larger sample size and thus a larger number of grains can undoubtedly better represent the target polycrystal, but the massive sample size is limited by computer memory, and often this is not as large as we would wish because of the intense demands on computer resources of the simulations. An alternative way of achieving a better representation is to use the combination of multiple realizations that are independently generated but follow the same statistics. This is comparable to multiple experimental measurements, performed at various locations of the same sample, that cover different volumes of grains. Multiple realizations can be generated by re-randomizing either grain seeds, crystallographic orientations, or both. It has been demonstrated^{21} that these procedures are statistically equivalent, and this work uses the simplest one that reshuffles crystallographic orientations only.

Thus, the present studies make use of a model volume that is inhomogeneous and anisotropic at grain scale but homogeneous and isotropic in both stiffness and geometry at macro scale. However, the principle established here may be extended to elongated grains and preferred crystallographic texture materials.

### B. Approximation of wave equations by FE

The discretization of a volume of polycrystal material for the simulation of wave propagation has been discussed previously,^{21,31} and we follow closely the setup reported in Cook *et al.*^{37} The spatial representation process is illustrated in Fig. 1, with details in the caption.

While the spatial problem is represented by FE discretization using first-order 8-node regular hexahedral elements, the time-stepping solution is based on the finite difference method, using the well-known equation,

in which $\Delta t$ is the duration of discrete uniform time steps, the superscript $n$ represents the current time step, $F$ is the dynamic load vector describing externally applied forces at individual degrees of freedom, $M$ and $K$ are the mass and stiffness matrices, so $KUn$ is the elastic force vector, with $Un$ representing the global displacement vector at time step $n$. $MU\xa8n$ is the inertia force, and the acceleration vector $U\xa8n$ therein is represented by a standard central difference as $(Un+1\u22122Un+Un\u22121)/\Delta t2$. The equilibrium equation establishes the relation of displacements at three time steps $n\u22121$, $n$, and $n+1$. Re-arranging the equation allows the displacements at step $n+1$ to be determined explicitly from those already known at the two preceding steps. We neglect material damping in our studies.

The above FE formulation involves two major approximations that need to be carefully treated. The first is the discretization of model space into elements, which we choose to be perfect identical cubes, all with an edge size of $h$, and this spatial discretization affects grain representation and numerical accuracy. The design of FE meshes to achieve a satisfactory convergence has been discussed widely in the literature, and it is important to note here that we limit our assessment to our specific present context of the solution of wave propagation in the time domain using the particular spatial and temporal schemes set out here. In other contexts, the convergence can be quite different, for example, Langer *et al.*^{38} show that several hundred elements per wavelength may be needed to achieve satisfactory convergence in the FE calculation of the eigenfrequencies of even a simple structure such as a plate in flexure. Structural vibration problems can incorporate quite complex spatial behavior that places additional demands beyond those of simple plane waves, and indeed this was demonstrated in the same paper by comparison of the plate flexural problem with a wave propagation problem that was much less demanding on spatial discretization. Furthermore, there are differences between the eigenvalue solution for frequency domain calculations and the time stepping for time domain problems.

In previous work using the same methodology as we discuss here, we found that a good representation of grains requires a minimum of ten elements per average grain size *a* (specifically the cubic root of average grain volume) and a satisfactory numerical accuracy requires at least ten elements per wavelength $\lambda $,^{21,31} the latter following on from previous work and prior reports by others. For now, we set out with this mesh refinement in mind, but we will evaluate in detail the accuracy both in wave velocity and wave amplitude in Sec. IV A. We note that increasing either grain or wavelength sampling rate would help achieve better numerical accuracy and convergence but would also increase FE computation intensity, and ten elements per grain dimension or wavelength is a moderate choice within our current computation capability. For the former criterion, the required mesh size of a given polycrystal can be determined without ambiguity from $a/h\u226510$. For the latter, however, wavelength $\lambda $ changes with frequency $f$ and phase velocity $V$. It is convenient to define a spatial sampling number as the ratio of wavelength to mesh size, $S=\lambda /h=V/(fh)$. Thus, the latter discretization criterion can be re-defined as $S\u226510$, which should be met by the smallest wavelength that corresponds to the highest simulation frequency and the slowest phase velocity.

The second approximation is the sampling of time into discrete steps with a spacing of $\Delta t$. This temporal sampling and the spatial discretization collectively determine the stability of FE calculations, which is prescribed by the Courant-Friedrichs-Levy condition. Due to the use of explicit time-marching in Eq. (1), the condition of stability requires the Courant number being no larger than unity, i.e., $C\u22641$, which is to assure that the distance travelled by the wave in one time step must not be greater than the distance between two mesh points. In the present work, the Courant number is given by $C=Vmax\Delta t/h$, where $Vmax$ is the fastest phase velocity in the simulated material, and for cubic materials considered in this work, it is given by^{39} $Vmax=(c11+2c12+4c44)/(3\rho )$.

### C. Excitation of plane longitudinal waves

Boundary and loading conditions are defined to excite the desired plane longitudinal waves, as shown in Fig. 1, with details in the caption. The lateral boundary conditions that are prescribed so that the finite domain represents an infinitely wide one can be achieved^{21} by applying periodic boundary conditions (PBCs) or SBCs to the four lateral boundaries: $x=0$, $x=dx$, $y=0$, and $y=dy$. Both have been shown^{21} to perform equally well for polycrystals, and due to the ease of implementation in the FE method, the SBCs approach is chosen in this work. The SBCs constrain the out-of-plane displacements of the lateral boundaries, i.e., $ux(x=0/dx,y,z)=uy(x,y=0/dy,z)=0$.

Second, loading conditions are imposed to create longitudinal waves of plane wavefronts. This is implemented by applying a uniform dynamic force to the external surface of $z=0$. A $z$-direction force, in the time domain being a three-cycle Hann-windowed toneburst of a center frequency of $fc$, is exerted to all nodes on the surface $z=0$. Due to the SBCs, this simulates a desired plane longitudinal wave travelling in the $z$-direction.

## III. EFFECTIVE MEDIA PARAMETERS FOR COHERENT WAVE: ATTENUATION AND PHASE VELOCITY

Although the initiation of the wave is coherent across the wavefront, the coherence is lost as the wave propagates through the polycrystal. The wavefront degenerates, scattering energy locally into incoherent diffuse components. Consequently, coherent waves are attenuated, and their phase velocities vary with frequency. These two scattering behaviors are parametrically described by attenuation coefficient $\alpha (f)$ and dispersive phase velocity $V(f)$, both as functions of frequency $f$.

### A. Coherent waves

To calculate these two effective media parameters, displacements are numerically measured on planes normal to the propagation direction. In the case of plane longitudinal waves travelling in the $z$-direction as addressed in this work, $z$-displacements are acquired on cross-sections normal to the $z$-direction. Note that the displacement needs to be divided by 2 if the cross-section is the stress-free end at $z=dz$, as the displacement is the sum of the incident and (equally) reflected waves. For an arbitrary cross-section $z=z0$, the obtained displacement field is denoted as $u(x,y,z0;t)$. The model length $dz$ for the $z$-direction will be used repeatedly hereafter, and for brevity, it will be simplified as $d$.

As a result of grain scattering, the phase of the wave $u(x,y,z0;t)$ varies across points on the cross-section. The coherent part of this wave can be obtained by averaging this wave across the cross-section;^{40} i.e., $U(z0;t)=\u3008u(x,y,z0;t)\u3009x,y$, where $\u3008\u22c5\u3009x,y$ represents the spatial average over all points on the surface $z=z0$. In the case of FE modelling, a significant number of nodes across the wavefront $z=z0$ contribute to the average and thus a significant wave information is included in the coherent wave. Fourier transforming this coherent wave with respect to time $t$ leads to its counterpart $U(z0;f)$ in the frequency domain. Upon obtaining the coherent waves for two or more cross-sections, the two effective media parameters can be conveniently calculated, which will be discussed in the following subsection.

In addition to the coherent part, there is the incoherent, fluctuation part in the total field $u(x,y,z0;t)$ which scatters with a random phase across the wavefront and therefore its average over the cross-section is nearly zero.^{40} Although this fluctuation part is discarded in the calculation of effective media parameters, it is used below to demonstrate the level of scattering over a cross-section. For a clearer demonstration, the fluctuation part is normalized by the coherent part and denoted as $uf(x,y,z0;t)=[u(x,y,z0;t)\u2212U(z0;t)]/U(z0;t)$.

Here we use an example to illustrate coherent waves and fluctuations on wavefronts. The example simulates a plane longitudinal wave travelling in the model N11520 (Table I) using the material properties of Inconel (Table II). The wave is excited with a center frequency of 5 MHz. The $z$-displacement fields are recorded and denoted as $u(x,y,0;t)$ and $u(x,y,d;t)$ for the transmitting $z=0$ and receiving $z=d$ surfaces, respectively.

In the time domain, the spatial averages of the displacement fields, namely, $U(0;t)=\u3008u(x,y,0;t)\u3009x,y$ and $U(d;t)=\u3008u(x,y,d;t)\u3009x,y$, are calculated and plotted in Fig. 2(a). For these two coherent signals, Figs. 2(b) and 2(c) provide their corresponding fluctuations $uf(x,y,0;t1)$ and $uf(x,y,d;t2)$ at the time instances of $t1=0.30$ µs and $t2=1.98$ µs, respectively. As shown in Fig. 2(a), the time instances $t1$ and $t2$ are chosen such that the coherent signals are at their maxima.

To obtain the spectral fields $u(x,y,0;f)$ and $u(x,y,d;f)$, the time-domain fields are cropped by using the time windows that span from $t=0$ to the vertical marks shown in Fig. 2(a) and then transformed to the frequency domain by utilizing the Fourier transform. Spatial averaging these spectral fields leads to the coherent parts $U(0;f)=\u3008u(x,y,0;f)\u3009x,y$ and $U(d;f)=\u3008u(x,y,d;f)\u3009x,y$, whose amplitude spectra are plotted in Fig. 2(d). At the frequency of $f0$ where the transmitting amplitude is at its maximum, the amplitude fluctuations $uf(x,y,0;f0)$ and $uf(x,y,d;f0)$ are illustrated in Figs. 2(e) and 2(f).

It is clear from Figs. 2(b) and 2(e) that the transmitting wavefront at $z=0$ is not planar. This is because we apply forces on this plane and then measure displacements at the same location where random elastic deformations occur immediately upon application of forces. As a result of multiple scattering, the field fluctuation on the receiving wavefront [Figs. 2(c) and 2(f)] is substantially larger than that of the source, leading to a smaller coherent wave (spatial average) on the receiving surface, as was seen in Figs. 2(a) and 2(d). As previously explained, the change of a coherent wave with propagation distance is represented by attenuation and dispersive phase velocity. The calculation of these parameters from the evolution of coherent waves will be addressed in the following subsection.

### B. Determination of attenuation and phase velocity

Before determining the scattering parameters, we demonstrate how coherent waves change as they propagate through polycrystal media. In Fig. 3, example coherent waves in polycrystalline aluminum and Inconel are given as functions of propagation distances $z$. In the figure, a coherent wave $U(z;f)$ is normalized by its corresponding initiating wave $U(0;f)$, and the amplitude and phase parts of the normalized coherent wave, $U(z;f)/U(0;f)$, are denoted as $A(z;f)$ and $\phi (z;f)$, respectively. These two parts are shown separately in Figs. 3(a) and 3(b).

Despite the large wavefront fluctuations (illustrated in Fig. 2), Fig. 3 shows that the amplitude and phase parts of the coherent wave change exponentially and linearly with propagation distance, respectively. It is important to emphasize that these exponential and linear relations are nearly perfectly described by their corresponding fitted straight lines (the goodness of fit is characterized by the annotated $R2$ measure, and $R2=1$ corresponds to a perfect fit). This is essential evidence that the change of the coherent wave with distance can be very well represented by the effective media parameters, namely, attenuation coefficient $\alpha (f)$ and phase velocity $V(f)$, by the following relations:

It is clear from Eq. (2) that attenuation coefficient and phase velocity can be determined in two different ways; these will be described in the following subsections.

#### 1. Measurement by through transmission method

The first way to measure attenuation and phase velocity is to use only $U(0;f)$ and $U(d;f)$, which are collected on the transmitting $z=0$ and receiving $z=d$ surfaces, see Fig. 1. This method is consistent with experiments where measurements can only be made on free boundaries. This method is named the through transmission (TT) method in this work and the two scattering parameters can be given explicitly as

where $A(d;f)$ and $\phi (d;f)$ are the amplitude and phase parts of $U(d;f)/U(0;f)$, respectively. This method is graphically shown in Fig. 3 as the solid (for aluminum) and dash (for Inconel) lines connecting the transmitting points at $z=0$ to the receiving ones at $z=d=10$ mm.

#### 2. Measurement by fitting method

The second way is to use all the coherent waves, measured on free boundaries as well as inside models. This can be fairly easily achieved by fitting the Eq. (3) to all the coherent waves, using attenuation coefficient $\alpha (f)$ and phase velocity $V(f)$ as fitting parameters. Similarly, this method is shown in Fig. 3 as the dashed-dotted (for aluminum) and dotted (for Inconel) lines that best fit all the data points in a least-squares sense.

#### 3. Comparison

Previously, the authors^{11,19,21,29,31} used the through transmission method, while Ryzy *et al.*^{32} employed the fitting method. In comparison, the results in Fig. 3 show a very small but visible systematic difference for Inconel. The methods are further compared in Fig. 4 for aluminum and Inconel over the frequency range of 3–6 MHz; details of the models are given in the caption. Subtle differences can be observed between the two methods: for the attenuation coefficient, the relative difference is smaller than 2%, and for phase velocity, it is smaller than 0.04%. Although the differences are small, it seems that the TT method gives systematically larger attenuation coefficients and smaller phase velocities than the fitting method. Thus, it is important to understand the reason for this seemingly systematic difference.

There are two key differences between these methods. (1) The TT method uses waves collected on traction-free boundaries, while the majority of the waves of its counterpart are acquired inside models; this difference is illustrated in Fig. 5. Examples are provided in Fig. 6 to illustrate the difference between measurements made on boundaries and inside models. In the figure, circular data sets are acquired inside the model N11520 at the frequency of 5 MHz, with each point and error bar representing the mean and standard deviation of 15 realizations. Square data sets are collected on the same surfaces and under the same modelling conditions, but the collection surfaces are free boundaries that are formed by truncating the model N11520 at these measuring surfaces. Thus, the grains are identical for each of the two models under comparison at each propagation distance. The figure shows that the signal amplitudes on free boundaries are mostly smaller than those inside models, hence the resulting attenuation coefficient from TT method is inevitably larger than that from the fitting.

This difference is due to the fact that on boundaries and inside models, incoherent waves contribute to the coherent signals differently (see Fig. 5). It is not hard to imagine that on a measuring surface the incoherent waves induced by early arrivals enter the coherent part of the arrival signals. Inside models, coherent waves harvest incoherent ones in both forward- and backward-propagating directions, while on free boundaries, incoherent backward-scattered waves are the main contributors. As a result, coherent waves collected inside models have slightly larger amplitudes than those acquired on the corresponding surfaces that are traction-free. Also, additional mode conversion of scattered waves on a free boundary may contribute to the difference. The incidence of random waves on complex free boundaries with spatially-varied elastic properties produces complicated reflections that involve mode conversions. Such mode conversions occur differently on a free boundary compared to the equivalent plane inside the material.

(2) The second aspect is that the TT method reflects the overall scattering behavior of the whole cuboid of a polycrystal medium, while its counterpart is largely influenced by the smaller subsets of the medium. It suggests that these two methods may have different statistical significances in representing the target polycrystal. This volume effect, however, is not the main cause of the difference between the two methods. The reason can be clearly seen from Fig. 6, which illustrates the use of identical volumes at each propagation distance still leads to statistically different results between the two methods.

However, it is important to realize that the difference between the two methods is generally within the range of error bars as shown in Fig. 6. This suggests that the TT and fitting methods are practically equivalent and both can be used. In this work, the TT method is used in later discussions due to its convenience in the collection of waves and in the calculation of scattering parameters.

## IV. NUMERICAL ERRORS AND STATISTICAL UNCERTAINTIES

It is certainly true that disagreement exists between results calculated by the FE method and the exact solutions. This disagreement should be estimated to establish the validity of the FE model and to aid improvement of its accuracy and precision.

In the case of modelling elastic wave propagation in polycrystals, the disagreement can be divided into two parts. One part is the systematic difference between an FE result and the exact solution, and another is the random variation around the FE result. These two parts are denoted as error and uncertainty, respectively. As mentioned in Sec. II A, we obtain an FE result by averaging those of multiple FE realizations that are independently generated but follow the same microstructure statistic. Error in this case is the difference between the average result and the exact solution, while uncertainty is the standard deviation of the multiple results relative to the average. Error and uncertainty arise in the FE model mainly as results of numerical approximations and statistical fluctuations.

Numerical approximations include two major facets: one is the approximation of an exact wave propagation problem by an FE model that is discretized in space and time, and another is the truncation or rounding of an exact number to a computer number that has a finite precision. These approximations cause numerical errors, and the estimation of these errors will be provided in Sec. IV A.

Statistical fluctuations relate to the random nature of polycrystalline microstructure and the stochastic scattering it causes. A combination of multiple realizations is used in the FE model to average out the fluctuations and thus to obtain statistically meaningful FE results. Such statistical considerations are accompanied by statistical errors and uncertainties, which will be estimated in Sec. IV B.

Following the estimations of errors and uncertainties in Secs. IV A and IV B, we examine in Sec. IV C the propagation of errors and uncertainties through coherent waves to effective media parameters, namely, attenuation coefficient and phase velocity. Finally, we provide examples in Sec. IV D to illustrate some practical approaches to suppress errors and uncertainties. In this work, absolute and relative errors are denoted, respectively, as $\Delta $ and $\delta $, while absolute and relative uncertainties are denoted as $\sigma $ and $\epsilon $.

### A. Errors in numerical approximations

The spatial and temporal discretization approximations in the FE method introduce numerical errors, and here we look at both space and time. Our interest is to establish knowledge of the underlying errors in amplitude and in wave velocity that come from the FE discretization of the propagation behavior; key to this is whether these errors are significant in comparison with, respectively, the changes of amplitudes due to scattering induced attenuation and the changes of phase velocity due to dispersion by the polycrystal. Such errors, if sufficiently understood, may also be used in practice in order to make corrections to simulated results. Therefore, we choose to assess the errors in a homogeneous isotropic domain, being the background on which the polycrystal behavior will be built, keeping with the same meshes and time stepping schemes of the latter. We recognize that this omits possible errors that might be specific to the discretization of the spatial changes in the polycrystal, but we suggest that the target of achieving good recommendations for choices of mesh and timestep is addressed most usefully by studying the non-specific homogeneous background case. Furthermore, we believe that a good discretization for the homogeneous background will automatically ensure good representation in the polycrystal case because it is already well established that the wave behavior is not sensitive to specific features at the space scale of the crystal boundaries, and indeed the analytical models do not include information about these boundaries at all.^{11}

The evaluation of the accuracy of FE models for wave propagation is a highly developed subject with a long history. Considerations of phase velocity dispersion in particular have motivated many different spatial and temporal schemes; dispersion is strongest for coarse meshes, so the challenge is to achieve best efficiency for a desired performance in modeling short wavelength waves. The monograph by Bathe^{41} provides an excellent description and evaluation of the basic methodologies of the main schemes in use for time domain simulation. A variety of studies of the performance of these and their variants, limiting to relevance to the method used in the present work, can be found in, e.g., Refs. 34 and 42–50.

We limit our reporting here to our specific choice of scheme and include results in brief for the benefit of modelers in the applications context of this article, leaving detailed evaluation of the mathematics and performance of the scheme to the cited literature. Our specific choice is the central difference explicit scheme, using lumped mass matrices, no material damping, and a structured mesh of identically-sized cube-shaped linear hexahedral elements. The basis of this choice has been optimized and explained in our previous work.^{21,31}

#### 1. Analytical prediction of numerical errors

The exact solution to the propagation of elastic plane waves within an isotropic non-damping material is given by $u(x,t)=p0exp[i(k0x\u22122\pi ft)]$, where $p0$ and $k0=k0n$ are the exact polarization and wave vectors, respectively. It is implied in the solution that the displacement amplitude is unity. There is no dissipation and the wave number $k0=2\pi f/V0$ is thus a real number where $V0$ is the exact phase velocity.

Due to numerical approximations, the FE solution to the same problem is approximate. At a given node $x$ and time instance $t=n\Delta t$, the approximate solution $uxn$ takes a similar form to the exact one

where $p$ and $k=kn$ are the approximate polarization and wave vectors. Numerical attenuation $\alpha $, if any, and phase velocity $V$, are included in the complex wave number $k$, i.e., $k=2\pi f/V+i\alpha $. By using Eq. (3), the wave number can be further related to amplitude and phase changes via $k=\phi /d\u2212i\u2009ln\u2009A/d$, where $d$ is the wave propagation distance. Thus, the estimation of amplitude and phase errors depends on the solution of wave number $k$.

The central difference explicit scheme is known to be non-dissipative, correctly conserving the energy of the propagating waves,^{49,50} although the literature devotes little attention to the details of this, concentrating in preference to the real concern, which is the velocity dispersion. Here, we suggest it is useful to think a bit more about the non-dissipative property before we look at the dispersion.

The non-dissipative property can be examined by looking at the numerical wave propagation solution for a monochromatic plane wave in an infinite domain (see also, e.g., Refs. 43 and 50). The elements are all identical, so we arbitrarily consider a single node, with waves propagating in the $z$ direction and no external forces. The scheme requires only local information at each time step, so we need only consider information from the eight elements that adjoin this node. This local assembly has mass $M$ and stiffness $K$ matrices with 81 × 81 coefficients, as defined in Sec. II B. Substituting Eq. (4) into Eq. (1) gives a system of 81 equations. However, the eight-element assembly only carries complete information about nodal connectivity for the target node. As a result, only the three equations which correspond to the degrees of freedom of this target node are meaningful in terms of describing wave motion, and, considering compression waves, only the displacement in the chosen coordinate $z$ relates to the example wave propagating along $z$. Simplifying this equation results in

and solving this equation leads to

It is clear from Eq. (6) that $k$ is a real number. This confirms the non-dissipative nature of the solution, which holds for plane waves that are continuous in space and time. A similar analysis would show the same for shear waves.

Equation (6) also shows the relative phase error, which is given by $\delta \phi =(k\u2212k0)/k0$. We define the spatial sampling and Courant numbers with respect to the longitudinal phase velocity $V0=VL$, i.e., $S=VL/(fh)$ and $C=\Delta tVL/h$. Considering these relations, relative phase error can be expressed as

The equation shows that phase error is directly related to spatial sampling number $S$ and Courant number $C$, which correspond to spatial discretization and time stepping, respectively. It follows now that we also have anisotropic behavior because, however well the element behaves for waves at different angles, the spatial sampling is certainly different in different directions; this is a well-known property of these solutions.^{34,43,47} Thus, we see that we should expect perfect preservation of the amplitude of monochromatic plane waves, but a frequency-dependent error of phase velocity. The latter will imply that the amplitude of the envelope of a wave packet in the time domain will in general exhibit distortion as it travels. We will evaluate this performance in simulation studies in the following two sub-sections.

#### 2. Evaluation of numerical phase errors

A simple cuboid model, with the dimensions of $dx\xd7dy\xd7d=12\xd712\xd710$ mm, is established using isotropic steel, with Young's modulus $E=210$ GPa, Poisson's ratio $v=0.3$, and density $\rho =8000$ kg/m^{3}. The longitudinal wave velocity of this material is $VL=5944.45$ m/s. Similar to the FE model given in Sec. II, the model is discretized with uniform eight-node cube elements and configured to accommodate a plane longitudinal wave that has the center frequency of 5 MHz. The transmitting and receiving waves, $U(0;f)$ and $U(d;f)$, are collected and the relative phase error is thus obtained as $\delta \phi =\phi /(2\pi fd/VL)\u22121$, where $\phi $ is the phase part of $U(d;f)/U(0;f)$.

Figure 7 compares the actual simulations with analytically predicted results from Eq. (7). Figure 7(a) assesses the influence of spatial discretization, as characterized by sampling number $S$, while Fig. 7(b) evaluates the impact of time stepping, as characterized by the Courant number $C$. Model details are given in the figure caption. The figure shows that actual phase errors can be well represented by Eq. (7).

A key observation from the figure is the high degree of accuracy of the calculated phase. The lowest sampling number, $S$, in the results in Fig. 7(a) of the figure is ten elements per wavelength, for which the phase error (and thus velocity error) is about 0.5%, and this error is reduced to 0.1% by increasing $S$ to 20 elements per wavelength. Thus, the policy of using at least ten elements per wavelength is justified for all routine analysis by this method unless particularly high precision is to be pursued.

The figure also defines the two approaches to reduce phase errors, pertaining to spatial discretization and time sampling, respectively. The first is to use finer meshes, which can be discovered from Fig. 7(a) that phase error decreases exponentially $[\delta \phi \u221dexp(\u22122S)]$ as spatial sampling number increases. The second approach is to use larger Courant numbers, which can be found from Fig. 7(b). When this number is at its limit of unity, phase error can nearly be removed. However, for stability reasons it is not recommended to use $C=1$ because the actual $C$ in a simulation could be greater than 1 (thus violate the stability condition) as a result of accelerated wave velocity caused by numerical dispersion, and also this condition cannot be achieved for all waves when there are spatial variations of material properties or multiple wave modes. An important implication of the second approach is that a structured mesh is likely to perform better than a free mesh in this analysis. This is because the free mesh has a range of element sizes and the range can be very wide for sharp spatial features, and therefore, it is impossible to achieve the same Courant number everywhere, even if there is only one wave mode and material present. In Sec. IV D 1, examples will be given to show actual errors in phase velocity calculations and to illustrate the corrections of these errors by using the two approaches discussed here.

#### 3. Evaluation of numerical amplitude errors

We use the same models as for Fig. 7(a) and calculate their relative amplitude errors by $\delta A=A\u22121$, where $A$ is the amplitude of $U(d;f)/U(0;f)$ and its exact solution is unity in the absence of amplitude error. In the frequency range of $f=$ 3–6 MHz, the errors are plotted against frequency in Fig. 8.

Figure 8 shows that the amplitude errors of the FE simulations are small and that they reduce as the mesh is refined, but they are not zero. We recall here that we expect zero loss for monochromatic plane waves, but that analysis did not include the FE calculations for the transient behavior of a wave packet. Critically, we need to distinguish here between energy loss, which would violate the non-dissipative nature, and signal distortion without loss of energy. Therefore, we have calculated the integral of the energy content for the whole wave packet in each case and this revealed that the energy loss is negligible, at around 10^{−4}% for the coarsest mesh. Thus, what we see here is that the energy is correctly retained, so the solution is convincingly non-dissipative, but that there is some small distortion of the amplitude spectrum, resulting in time-shifting of some spectral amplitudes in the wave packet.

For interest, Fig. 8 also shows the normalized amplitude spectrum of the transmitting wave for the model with the mesh size of $h=0.100$ mm, which can be seen, when scaled appropriately (right *y*-axis), to match the same shape as the spectrum of the distortion. Further work on this observation could be academically interesting, but this would be beyond the present scope. For now, we note that the wave packet does not dissipate energy, but that there is a small distortion across the amplitude spectrum; furthermore, the latter is very small (0.2%), even for the coarsest mesh ($S=12$ elements per wavelength at center frequency), while for the preferred finer meshes it is below 0.003% across the bandwidth ($S=60$ elements per wavelength at center frequency). These are very small numbers compared to the scattering induced attenuation in the polycrystal.

Thus, the analysis and practice both show that amplitude errors are insignificant in comparison to phase ones. Therefore, we can expect that limitations of the FE modelling will impact the predicted phase velocity but have a negligible influence on the wave amplitude when reported in the frequency domain. The key result to achieve accurate simulation is the phase error in Fig. 7(a), showing dramatic improvement as the number of elements per wavelength, $S$, is increased. However, this is in direct conflict to the desire to minimize $S$ for large models and high frequency simulations.

### B. Errors and uncertainties in statistical considerations

Generally, we aim to simulate polycrystals of prescribed statistical properties. A single FE model is not capable of fully describing one of such polycrystals, due to its relatively small sample volume and finite number of grains. Instead, the statistical combination of the results of multiple realizations that are randomly generated following the same statistical property is desirable for drawing meaningful conclusions about the given polycrystal.

The effectiveness of using multiple realizations to achieve statistical convergence is exemplified here. The model N11520 is used to simulate polycrystalline aluminum and Inconel that have statistically isotropic and macroscopically homogeneous properties. Due to random scattering, the wave field $u(x,y,d;t)$ on the receiving plane $z=d$ should be randomly distributed. This means that at a given time instance of $t=t0$ the fluctuation $uf(x,y,d;t0)$, which is defined in Sec. III A, should follow the Gaussian distribution. However, the actual distribution of a single realization differs from the Gaussian for both aluminum and Inconel, which are, respectively, shown in Figs. 9(a) and 9(b) as probability density histograms. Instead, the combination of multiple realizations (15 in this example) delivers satisfactory distributions that agree very well with the Gaussian. Thus, a combination of multiple realizations can be regarded as a confident description of the target polycrystal.

As mentioned in Sec. III A, we use the coherent wave of a model to calculate effective media parameters. In the case of multiple realizations, the calculation needs to use the average of the multiple coherent waves that can be denoted as $\mu (z;f)=\u3008Ui(z;f)\u3009i$, where $Ui(z;f)$ is the coherent wave of the $i$ th realization. As prescribed by the central limit theorem, the multiple coherent waves follow the Gaussian distribution because the multiple realizations are independently randomly generated. Also, the average $\mu (z;f)$ is an appropriate estimate of the true coherent wave, and this estimation includes statistical uncertainty and error as discussed below.

Uncertainty is characterized by the standard deviation $\sigma (z;f)$ of the coherent waves. Here, we use the 15-realization example of Fig. 9 to investigate the uncertainties associated with our simulations. For the received coherent waves $Ui(d;f)\u2009(i=1,2,\u2026,15)$, their uncertainties are calculated and split into amplitude $\sigma A$ and phase $\sigma \phi $ parts, respectively; normalization by their respective mean values results in relative uncertainties $\epsilon A=\sigma A/\mu A$ and $\epsilon \phi =\sigma \phi /\mu \phi $. The results show that Inconel has a larger uncertainty than aluminum, which is consistent with Fig. 9. This means that uncertainty is positively correlated with material anisotropy. In addition, the amplitude uncertainty is around ten times larger than that of phase. This is the opposite of what we found for numerical errors, which was that amplitude error is in general smaller than phase error; this was shown in Fig. 8(b).

Error in this case is termed the standard error of the mean. It is a statistical error between the mean and the true value, and it can be probabilistically estimated by $\sigma \mu =\sigma /N$, where $N$ is the number of realizations.

It is important to realize that uncertainty cannot be suppressed, but rather it can be better determined by using more realizations for example. For standard error, however, its expression indicates that it can be reduced by increasing the number of realizations $N$. Practically, a simple trial-and-error process can be employed to suppress standard error to a satisfactory extent. This process increases the number $N$ until the relative standard error of either amplitude or phase $(\delta \mu =\epsilon A,\phi /N)$ is below a set threshold (e.g., 0.1%).

### C. Propagation of errors and uncertainties

It is known that errors and uncertainties propagate through independent variables to dependent ones. Thus, those of coherent waves discussed in Secs. IV A and IV B will be brought into scattering parameters, namely, attenuation coefficient and phase velocity, via Eq. (3).

Given that propagation distance and frequency are accurate, the attenuation coefficient in Eq. (3) is a univariate function of amplitude. We denote amplitude error as $\Delta A=A\u2212AT$, where $AT$ is the true amplitude and $A$ is its estimation given by the mean of multiple realizations. Replacing $A$ by $AT$ in Eq. (3) and expanding this equation about the point $A$ in a Taylor series, we can then obtain the relative error of attenuation coefficient

where $\alpha $ is the true value of attenuation coefficient; frequency $f$ is implicit. In general, the true attenuation coefficient $\alpha $ is unknown, so instead, it is approximately replaced by the mean value of multiple realizations. Note that only the first-order term of the Taylor expansion is used in the equation by assuming that amplitude error is small. Higher-order terms should be included if this assumption does not hold. For uncertainty, we need to express the attenuation coefficient of each realization as a Taylor series. Substituting the expansion into the formula of standard deviation, we can then obtain the relative uncertainty

where $\sigma A$ is amplitude uncertainty. Similarly, assuming that phase error and uncertainty are small, then their propagation from phase to velocity can be easily obtained from Eq. (3),

where $V$ is the true phase velocity that is approximated to the mean value of multiple realizations.

As shown in Eqs. (8) and (9), the error and uncertainty of the attenuation coefficient are related to those of amplitude by the same factor, $F\alpha =exp(\alpha d)/(\alpha d)$. Due to the lack of true values, such propagation of error cannot be evaluated. However, it is easy to do so for uncertainty. For this purpose, $\sigma A$ and $\epsilon \alpha $ are calculated from 15 realizations and the actual propagation factor is calculated from their division. The resulting factor is plotted against $\alpha d$ in Fig. 10(a) for aluminum and Inconel, and each material covers a wide frequency range by using three models: N115200, N11520, and N16000. The actual uncertainty propagation can be seen to be very well represented by the analytical factor $F\alpha $, and this should always be valid as long as amplitude uncertainty is small. Likewise, this factor is also applicable to error propagation on the same condition that amplitude error is small. The figure conveys important information that the propagation factor is at its minimum of $e=2.72$ when $\alpha d=1$ neper. Thus, in order to suppress the magnification of error and uncertainty, the total attenuation $\alpha d$ should be kept around 1 neper^{51} Practically, we have observed that moderate magnifications can be achieved by maintaining $\alpha d$ between 0.01 and 6 neper.

For phase velocity, Eq. (10) shows that the same factor $FV=V/2\pi fd=1/\phi $ governs both error and uncertainty propagation. Similarly, uncertainty propagation is evaluated in Fig. 10(b) by using the same models as its left panel. The figure shows that the factor $FV$ agrees very well with actual simulations, and it is no doubt that it is also valid to error propagation. An elaborate analysis would suggest that error and uncertainty are naturally suppressed during their propagation. This is based on the fact that achieving a good statistical significance requires waves to propagate at least a quarter of a wavelength, which means $\phi \u22652\pi /4=\pi /2$ and further leads to $FV=1/\phi \u22642/\pi <1$. We note that in Fig. 10(b), the phase $\phi $ on the horizontal axis is in the range from 10 to 1000 and the propagation factor $FV$ on the vertical axis is not greater than 0.1. It is also important to emphasize that better suppression can be achieved by using a longer model length $d$. However, increasing $d$ can make it difficult to meet the condition of $\alpha d$ being close to 1 neper, which is necessary to suppress the error and uncertainty of attenuation. It is thus advised to only consider the suppression condition exerted to attenuation and that of phase velocity should be fulfilled accordingly.

### D. Example suppression of errors and uncertainties

The preceding studies have provided a variety of ways to understand and then suppress errors and uncertainties. Following these studies, this subsection presents some typical examples to demonstrate the practical use of the suppression methods in actual simulations.

Suppressing statistical error and uncertainty is relatively simple, and it is thus not exemplified in this subsection. As illustrated in Sec. IV B, statistical error can be reduced to a satisfactory extent by increasing the number of realizations $N$, while uncertainty cannot be suppressed and can only be better determined through the use of large $N$. For all simulations considered in this work, we use $N=15$ to achieve a satisfactory statistical error of at most 0.1% and a sufficiently good determination of uncertainty. Section IV C further suggests that the propagations of statistical error and uncertainty can be effectively suppressed by minimizing the factors $F\alpha $ and $FV$, respectively. However, because such minimization equally suppresses numerical errors, we leave its discussion to the following example suppressions of numerical errors.

Suppressing numerical errors requires a good combination of modelling parameters, and the suppression approach differs from phase velocity to attenuation. Thus, the following gives suppression examples separately for velocity and attenuation.

#### 1. Suppression of numerical velocity errors

As shown in Eq. (10), the relative error of phase velocity $\delta V$ is the product of phase error $\Delta \phi $ and propagation factor $FV$. Thus, reducing either of these two parts leads to the suppression of velocity error. For $FV$, Sec. IV C infers that it naturally reduces errors as it is smaller than unity. This factor can be further reduced by using a longer model length $d$. However, we do not aim to do so because $d$ is a critical parameter for suppressing attenuation error.

For $\Delta \phi $, we have obtained its analytical expression for homogeneous materials as given in Eq. (7). Although we have not seen a means yet to extend it analytically to polycrystals, we shall see that the conclusions drawn from it can be quantitatively applied to polycrystals for error suppression. To use this equation for polycrystals, we define its spatial sampling number $S$ with respect to the Voigt velocity $VL$, i.e., $S=VL/(fh)$, and its Courant number $C$ with respect to the peak velocity $Vmax$, namely, $C=\Delta tVmax/h$. For cubic materials considered in this work, the Voigt and peak velocities are given, respectively, as $VL=(3c11+2c12+4c44)/(5\rho )$^{11} and $Vmax=(c11+2c12+4c44)/(3\rho )$.^{39}

Reducing phase error $\Delta \phi $ can be achieved by using the two approaches summarized in Sec. IV A 2. The first approach is to use a larger spatial sampling number $S$. However, this approach is practically difficult for the present example. This is because the three material models (Table I) used in this work are already very big in terms of their required computational resources, and further increasing their $S$ requires the use of even finer meshes that would easily exceed our current computer capability. Instead, this first approach is substituted by a correction approach that uses the numerical phase error predicted from Eq. (7) to correct the error included in simulation result. The second approach is to control phase error $\Delta \phi $ by using a Courant number that is as close to unity as possible.

Examples are provided in Fig. 11 to illustrate these two approaches. The examples use the same FE models and configurations as Fig. 10. In the figure, the phase velocities obtained from FE simulations are provided as discrete points and theoretically predicted results are given as lines. The accurate statistical information of the FE material models is incorporated into the theoretical curves (see Sec. V for details), and these curves are of good accuracy for the studied materials and frequency ranges. In particular, the curve of aluminum in Fig. 11(a) is very accurate because the underlying theory can take almost all its weak scattering into consideration. This curve can thus be used as a reference to evaluate the error suppression approaches.

For the first approach of error correction, phase velocities before correction are provided as hollow circles and corrected ones as solid circles. In this case, small Courant numbers (0.81 for aluminum and 0.86 for Inconel) are used in order to get large phase errors and thus to show clear effects of correction. The correction is achieved by subtracting the numerical phase error $\phi \delta \phi $ [$\delta \phi $ is the relative phase error calculated by Eq. (7)] from the simulation result $\phi $ and then calculating the corrected phase velocity by $V=2\pi fd/(\phi \u2212\phi \delta \phi )$. It is obvious that before correction, the three models (N115200, N11520, N16000) do not give consistent velocities as expected at the transitions between individual models. Note that the transitional inconsistency is less clear for Inconel [Fig. 11(b)] because its *y*-axis scale is too large to show this. This inconsistency is nicely averted in the corrected results, implying the effectiveness of the correction approach. Also, by comparing the hollow and solid circles with the theoretical curve of aluminum [Fig. 11(a)], we can observe that the use of the correction approach significantly suppresses numerical phase velocity errors.

For the second approach of error control, both aluminum and Inconel use a large Courant number of $C=0.99$. This is the largest $C$ number attempted in this study that does not encounter stability issues; the Courant numbers above $C=0.99$ were tested but led to unstable FE computations for both aluminum and Inconel. The velocity results are provided in the figure as solid squares. The results show good consistencies at model transitions, but not as good as the corrected results of the first approach. The theoretical curve of aluminum in Fig. 11(a) is also closer to the corrected points, suggesting that the correction approach performs better than the control approach. This error control approach, however, may be more advantageous for practical simulations as it requires less post-simulation calculations compared to the correction approach.

#### 2. Suppression of numerical attenuation errors

Equation (8) indicates that the relative error of attenuation $\delta \alpha $ is the multiplication of amplitude error $\Delta A$ and propagation factor $F\alpha $, and thus reducing either of these two parts leads to the suppression of attenuation error. It is suggested in Sec. IV A 3 that amplitude error $\Delta A$ can be reduced by using either a larger spatial sampling rate $S$ (thus, a finer mesh size $h$) or an appropriate spectrum range that is in close proximity to its maximum. Also, Sec. IV C indicates that propagation factor $F\alpha $ can be suppressed by limiting $\alpha d$ to a certain range in the vicinity of 1 neper.

Thus, three opportunities are available to suppress numerical attenuation errors. For the first opportunity of using a larger spatial sampling rate $S$, we have mentioned above in Sec. IV D 1 that the three models utilized in this work are already very computationally intensive, and therefore we do not aim to further increase their sampling rates. For the third opportunity of using a $\alpha d$ being close to 1 neper, the creation of the currently used FE models has ensured that this condition can be relatively easily satisfied. It is certainly true that we can create more models to accommodate different modelling frequencies and polycrystalline media, such that we always have an optimal suppression condition of $\alpha d\u22481$ neper. However, such treatment would introduce huge computation and time costs, and therefore it is not pursued in this example. Instead, this example uses a combination of all three opportunities, rather than any single one, to reduce attenuation errors.

The principle of combining the suppression opportunities is relatively simple. For a given material, trial simulations (or theoretical predictions) are first conducted to determine the frequency range of each FE model in which $\alpha d$ is within a bound of 1 neper, e.g., $0.01\u2264\alpha d\u22646$ neper. If there is an overlapped frequency range between two models, the overlapped range is designated to the model with either a larger sampling rate $S$, a closer attenuation $\alpha d$ to 1 neper, or both. Finally, within the chosen frequency range of each model, transmitting center-frequencies are tailored to ensure that the whole range uses appropriate spectra of the waves.

Here, we provide an example in Fig. 12 to demonstrate an incorrect combination of suppression opportunities and the corrected combination that appropriately reduces numerical attenuation errors. The example is given for the material of aluminum. In the frequency range of 6.5–13.5 MHz, Fig. 12(b) shows that the models N11520 (solid triangles) and N16000 (hollow circles) both satisfy the condition of $0.01\u2264\alpha d\u22646$ neper $(2.72\u2264F\alpha \u2264101.00)$. The incorrect combination chooses the model N16000 to cover this frequency range because it has a larger sampling rate $S$ than N11520, which is shown in Fig. 12(c). To save modelling costs, the incorrect combination excites a transmitting wave in the form of a three-cycle Hann-windowed toneburst with a center frequency of 20 MHz to provide a wide frequency coverage. In this case, the range of 6.5–13.5 MHz is far away from the maximum of the transmitted amplitude spectrum as shown in Fig. 12(d). This deviation from spectral maximum may cause a relatively small amplitude error $\Delta A$, but the error is significantly magnified by the large propagation factor $F\alpha $. This leads to incorrect attenuation results (hollow circles) as shown in Fig. 12(a). Instead, the corrected combination uses the model N11520 having smaller propagation factors in the frequency range of 6.5–13.5 MHz as shown in Fig. 12(b) and performs the FE simulations in an appropriate spectral range as shown in Fig. 12(d). The corrected results [solid triangles in Fig. 12(a)] agree very well with the theoretical curve, which is very accurate in the range of concern as explained in Sec. IV D 1, while the incorrect results deviate substantially from both corrected and theoretical results. It is important to realize that the peak difference between the incorrect and corrected results, with the latter as reference, is as large as −54.78% occurring at the frequency of 6.5 MHz, see Fig. 12(a).

## V. INCORPORATION OF STATISTICAL INFORMATION INTO THEORETICAL MODELS

The preceding FE method is capable of performing numerical experiments with fully-controlled conditions, and from the experiments, effective media parameters can be ideally measured with errors and uncertainties being appropriately estimated. Thus, this method is well suited to evaluating theoretical models that are generally approximate due to the use of assumptions and simplifications. However, there remains one more challenge for accurate simulations, which is to determine the set of characteristics that are defined by a chosen FE model. Clearly, the choices of most parameters, such as density, dimensions, or crystal properties are simply recorded. However, the spatial characteristics of the grains are a product of the creation of the polycrystal domain, and these need to be measured from the FE model. This is a vital step for any comparative studies with analytical models to ensure that these properties are represented properly in parallel in the FE and analytical models. In this third study, we will investigate this challenge.

### A. Statistical information required: TPC function

The evaluation of a theoretical model requires its input statistical characteristics of the material system being measured from FE material models and its output being assessed against modelling results. For this reason, we aim to describe FE material models in a way that closely matches the need for theoretical models. As we have identified in our previous studies,^{11,19,21} the most significant element is the statistical description of the spatial variation of the material properties, namely, the geometric TPC function that is generally required for the advanced theoretical models.

The geometric TPC function, $w(r)$, is the probability of two points separated by a distance of $r$ being in the same grain. A single-term exponential function, $w(r)=exp(\u2212r/a)$ ($a$ is the correlation length), was traditionally used to describe polycrystals with untextured, equiaxed grains.^{5,8} However, this simplified function is substantially different from the TPC curve measured from a typical polycrystal model,^{19,32} and therefore it is essential to incorporate measured TPC data into theoretical models for more accurate evaluations.

To measure the TPC $w(r)$, previous procedures^{19,32} drop a large number of random pairs of points, each pair being separated by a distance of $r$ [see Fig. 13(a)], into a polycrystal model. Then they count the number of pairs whose two points occur in the same grain and compare this number to the total number of pairs to obtain the probability $w(r)$. It is important to emphasize that a substantial number of tests, usually 1 × 10^{6} tests per discrete $r$, are required to obtain statistically converged TPC data.

To incorporate measured TPC data into a theoretical model, a mathematical function is constructed by best fitting the measured data. An optimal function of choice is the multi-term exponential function,^{19} $w(r)=\u2211iAiexp(\u2212r/ai)$, that has the advantage of delivering a simpler formulation for the theoretical model. The physical constraints on the coefficients $Ai$ and $ai$ have been summarized in Van Pamel *et al.*^{19} and the fitting can be achieved by using any nonlinear curve fitting algorithms.

In summary, the statistical information required for theoretical evaluations is the TPC function, which needs to be measured from FE models and fitted into a mathematical form.

### B. Effect of SBCs

As shown in Fig. 13(a), previous procedures^{19,32} to measure TPC data limit both points of a random pair to the inside of a model, and thus consider only the statistical information within the cuboid of the model. However, it is mentioned in Sec. II C that we use SBCs to emulate polycrystal media of infinite size in the transverse direction to that of wave propagation. Such conditions mirror the boundary grains that are associated with SBCs, which is demonstrated in Fig. 13(b), and these grains are therefore likely to be larger on average than the grains located within the body of the model. Also, in the FE computations, this results in stronger scattering of propagating waves when they encounter enlarged boundary grains. Thus, it might be desirable to consider this mirroring effect in the measurement of TPC data with a goal to evaluate, by comparison with the FE results, which TPC is more appropriate for evaluating theoretical results.

#### 1. Effect on TPC function

It can be inferred from the illustration in Fig. 13(b) that the statistical significance and impact of the boundary grain mirroring increases as the proportion of boundary grains grows. Based on this inference, now we evaluate the actual effect of SBCs on the TPC function. Achieving this evaluation requires the grains associated with SBCs being parametrically described and the TPC functions with consideration of SBCs being numerically measured.

To describe the grains associated with SBCs, we use three parameters: the number fraction, NF, and volume fraction, VF, describing, respectively, the relative number and volume of boundary grains to all grains within the model, and the surface area to volume ratio SA:V representing the amount of symmetry boundary area per unit model volume. We calculate these parameters for the models given in Table I and provide the results in Table III with each of the coordinate axes being wave propagation direction. The table shows that for the provided models having significant numbers of grains, these three parameters are actually equivalent, with NF and VF being approximately equal to each other and being proportional to SA:V. It is important to realize that NF is slightly smaller than VF because the grains generated by the Voronoi tessellation do not have uniform volume, but rather boundary grains are to some degree larger than internal grains. Since these three parameters are equivalent, the following evaluation only uses the surface area to volume ratio SA:V to represent SBCs because this parameter can be easily calculated without the knowledge of actual grains. Also, the studies reported below just consider, for the sake of clarity, the five cases in Table III with unique SA:V values: N115200-$x$, N11520-$x$, N11520-$z$, N16000-$x$, and N16000-$z$. Here, the label of each case is composed of its model name and propagation direction.

Model . | $x$-direction propagation
. | $y$-direction propagation
. | $z$-direction propagation
. | ||||||
---|---|---|---|---|---|---|---|---|---|

NF . | VF . | SA:V . | NF . | VF . | SA:V . | NF . | VF . | SA:V . | |

N115200 | 0.0860 | 0.0900 | 0.1867 | 0.0858 | 0.0896 | 0.1867 | 0.1482 | 0.1541 | 0.3333 |

N11520 | 0.1582 | 0.1714 | 0.3667 | 0.1614 | 0.1680 | 0.3667 | 0.1461 | 0.1555 | 0.3333 |

N16000 | 0.2244 | 0.2323 | 0.5000 | 0.2259 | 0.2329 | 0.5000 | 0.0908 | 0.0958 | 0.2000 |

Model . | $x$-direction propagation
. | $y$-direction propagation
. | $z$-direction propagation
. | ||||||
---|---|---|---|---|---|---|---|---|---|

NF . | VF . | SA:V . | NF . | VF . | SA:V . | NF . | VF . | SA:V . | |

N115200 | 0.0860 | 0.0900 | 0.1867 | 0.0858 | 0.0896 | 0.1867 | 0.1482 | 0.1541 | 0.3333 |

N11520 | 0.1582 | 0.1714 | 0.3667 | 0.1614 | 0.1680 | 0.3667 | 0.1461 | 0.1555 | 0.3333 |

N16000 | 0.2244 | 0.2323 | 0.5000 | 0.2259 | 0.2329 | 0.5000 | 0.0908 | 0.0958 | 0.2000 |

To measure the TPC functions with consideration of SBCs, we use a procedure that is mostly the same as the one summarized in Sec. V A except for the two treatments adopted to consider SBCs. The first treatment is that for a pair of random points, only one of them is confined within the model, while the other is completely random as long as their distance is $r$, meaning that the latter point can be outside of the model. Given a sufficiently large number of pairs, this allows the ergodicity to be satisfied, as the pairs can cover the entire space of infinite lateral extents that are prescribed by SBCs. The second treatment involves the counting of correlated pairs. As illustrated in Fig. 13(b), an extra situation of two points being correlated may occur when the second point is outside of the model. In this case, a pair is considered to be correlated if one point occurs in a symmetry boundary grain and the other falls into its mirrored counterpart. By using this measurement procedure, the TPC data of the five cases given above are obtained and plotted in Fig. 14(a) as hollow points (the data points for the case N115200-$z$ are also given in the figure for later use). To evaluate the influence of SBCs, the TPC statistics neglecting SBCs are also provided as solid points in the figure.

A clear difference can be seen in Fig. 14(a) between the two sets of TPC points with and without consideration of SBCs, indicating that the mirroring of symmetry boundary grains has a substantial influence on TPC statistics. When SBCs are not considered, the data points of different models are almost overlapped. This suggests that the models adopted in the present work are essentially statistically equivalent and are of satisfactory statistical significance. When SBCs are taken into consideration, the measured statistics show evident discrepancies, and these statistics are positively correlated with SA:V. The correlation is more pronounced at TPC tails, which reflects the increased grain sizes (or correlation lengths) caused by grain mirroring.

The simulated waves travel in the $z$-direction, meaning that the consideration of SBCs involves only the cases associated with $z$-direction: N115200-$z$, N11520-$z$ and N16000-$z$. It is shown in Fig. 14(a) that these three cases are quite close to each other, but N11520-$z$ has the largest tail and thus the strongest SBCs effect. Therefore, this case is used in the following as an example for the consideration of SBCs, and it is further plotted in Fig. 14(b) as hollow circles. For the case of neglecting SBCs, the TPC data of the three models [Fig. 14(a)] are nearly identical. Thus, their average, provided in Fig. 14(b) as solid circles, is employed to describe the case of neglecting SBCs.

It would be practically meaningful to reduce the impact of SBCs such that measuring TPC data can be performed in a simpler way by neglecting SBCs. This can be achieved by using smaller SA:Vs due to the positive correlation between TPC data and SA:V. It is shown in Table III that the way to get a smaller SA:V is through widening the lateral dimensions of polycrystal models, and in extreme cases, the effect of SBCs can be negligible if the lateral dimensions are sufficiently wide. However, the models adopted in the present work are already very big, and increasing their lateral dimensions will lead to added challenges to our computation resources. In addition, an alternative possibility of minimizing SBCs effect is to reduce the average size of the symmetry boundary grains such that the mirrored grains would have similar sizes to internal grains. It is practically difficult to generate grains with such a distribution but implementing a method for this aim would be very useful.

#### 2. Effect on attenuation and phase velocity

Having demonstrated in the preceding examination that SBCs have a significant effect on TPC functions, it is now essential to identify if SBCs would equally affect effective media parameters in terms of attenuation coefficient and phase velocity. For this purpose, we input the TPC functions, with and without consideration of SBCs, into one of the advanced theoretical models, the second-order approximation (SOA),^{19} to calculate the parameters. The degrees of agreement between the resulting parameters and FE results are then estimated to assess the improvement brought in by the consideration of SBCs.

Incorporating TPC functions into the SOA model requires the two sets of data in Fig. 14(b) to be fitted to mathematical forms. The process discussed in Sec. V A is used here and a multi-term exponential function is chosen as the fitting function. The fitting is performed in a least squares sense, and the number of exponential terms is selected as eight because the fitting can be unstable if a larger number is used. The coefficients of the fitted functions are listed in Table IV. The functions are further plotted in Fig. 14(b), showing reasonably good fittings.

Terms . | Without SBCs . | With SBCs . | ||
---|---|---|---|---|

$Ai$ . | $ai$ . | $Ai$ . | $ai$ . | |

1 | −2922.66 | 0.115726 | −129.485 | 0.1280002 |

2 | −10.6181 | 0.172221 | −388.413 | 0.1279995 |

3 | 3914.370 | 0.110790 | 797.091 | 0.1517330 |

4 | −3305.68 | 0.103201 | −1372.040 | 0.1898924 |

5 | 696.919 | 0.0910265 | 1982.690 | 0.1867733 |

6 | 54.3582 | 0.152263 | 694.300 | 0.1242404 |

7 | −42.4732 | 0.0723540 | −469.525 | 0.1279995 |

8 | 1616.78 | 0.110790 | −1113.610 | 0.1709402 |

Terms . | Without SBCs . | With SBCs . | ||
---|---|---|---|---|

$Ai$ . | $ai$ . | $Ai$ . | $ai$ . | |

1 | −2922.66 | 0.115726 | −129.485 | 0.1280002 |

2 | −10.6181 | 0.172221 | −388.413 | 0.1279995 |

3 | 3914.370 | 0.110790 | 797.091 | 0.1517330 |

4 | −3305.68 | 0.103201 | −1372.040 | 0.1898924 |

5 | 696.919 | 0.0910265 | 1982.690 | 0.1867733 |

6 | 54.3582 | 0.152263 | 694.300 | 0.1242404 |

7 | −42.4732 | 0.0723540 | −469.525 | 0.1279995 |

8 | 1616.78 | 0.110790 | −1113.610 | 0.1709402 |

Substituting the two fitted functions into the SOA model and solving the equation results in the theoretical solutions to attenuation and phase velocity. Here we use Inconel, the properties of which are given in Table II, as the material input to calculate theoretical solutions. The resulting attenuation and velocity curves are plotted in Figs. 15(a) and 15(b), respectively. In the figure, the corresponding finite element model (FEM) results are also provided as baseline data to evaluate the approximations of these solutions.

It is shown in Figs. 15(a) and 15(b) that the overall agreement between the theoretical and FEM results is very good, no matter whether SBCs are considered in the TPC functions or not. However, there are clear differences between the theoretical results, which demonstrate that the SBCs effect is actually prominent. In most frequency regions, the curves with consideration of SBCs agree better with FEM points, conveying an essential message that the FEM results are influenced by the doubling of grain size at symmetry boundaries. It is also important to notice that the theoretical and FEM results are not perfectly matched even when boundary grains are considered. This is because in the FEM results, the SBCs affect scattering only at boundaries, while in the theoretical model, the effect is not reproduced exactly as the TPC function assumes random, homogeneous distribution of those grains.

There is no doubt that the TPC function with SBCs should be selected for our material systems when evaluating the theoretical model. To quantify the improvement of theoretical prediction brought in by the consideration of SBCs, we calculate the relative differences of theoretical curves to the FEM points and plot them in Figs. 15(c) and 15(d). The agreement is clearly better at all frequencies for the TPC function with SBCs.

The conclusion of this investigation is that it is best to consider the SBCs effect in the TPC function. In practice, however, it is always desirable to reduce the effect of the boundaries to a negligible level by making a model as wide as possible. It is important to do so because one of the key interests of FE modelling at the moment is to check the significance of assumptions in analytical models, and thus there is a strong motivation to eliminate this TPC source of differences from the modelling.

## VI. SUMMARY AND CONCLUSION

This work presents a thorough study of the FE method for the propagation and scattering of elastic waves in polycrystals. Specifically, this work addresses the following three topics that have not been covered by the modelling studies so far but are essential to future FE simulations.

The determination of effective media parameters, namely, attenuation and phase velocity, is first discussed. These parameters describe the amplitude and phase changes of a plane wave as it travels through a polycrystalline medium. It is demonstrated that the travelling wave essentially does not have a planar wavefront, which is caused by random scattering. However, it is demonstrated that the coherent part, belonging to the plane wave, can be conveniently calculated as the spatial average of the wavefront. From a set of coherent waves, attenuation and phase velocity can be determined in two ways. The TT approach uses the coherent waves measured on two external surfaces: the transmitting surface and its opposite receiving surface. In addition to the two waves used in the TT method, the fitting approach also utilizes extra coherent waves collected inside a FE model. These two approaches are found to be practically equivalent and can thus be used interchangeably. The TT method might be preferable in actual simulations due to its ease of wave collection and calculation.

The estimation and suppression of the errors and uncertainties of effective media parameters are then attempted in this study. It is found that two factors can cause errors and uncertainties: numerical approximations and statistical considerations. Numerical approximations only cause errors, while statistical considerations involve both errors and uncertainties. For numerical errors, both attenuation and velocity errors can be suppressed by the use of a larger spatial sampling number $S$ (thus a finer mesh size $h$) if computational capability allows. In addition to this, numerical attenuation and velocity errors can also be reduced in different ways: (1) Numerical attenuation error can be suppressed by using an appropriate spectrum range that is in close proximity to its maximum, or by limiting attenuation $\alpha d$ to a certain range in the vicinity of 1 neper. (2) Numerical velocity error can be controlled by using a large Courant number $C$ that is close to unity and can also be corrected by utilizing the analytically predicted error from Eq. (7). For statistical errors and uncertainties, the suppression approaches developed for numerical errors also ensure that they can be satisfactorily controlled to a certain extent. The only requirement is that a sufficiently large number $N$ of realizations have to be used to obtain a reasonably small statistical error and a fairly converged statistical uncertainty.

The incorporation of FE model information into theoretical models is finally discussed. It is indicated that the TPC function is the statistic required to put into theories. The SBCs, employed to emulate infinitely wide media that can hold plane waves, have a considerable effect on the TPC function and the theoretically predicted attenuation and velocity. Thus, this SBCs effect needs to be considered in the TPC function for a more accurate incorporation of statistical information.

The FE method discussed in this work focuses on the specific case of plane longitudinal waves travelling in untextured polycrystals with equiaxed grains. However, it is expected that the principle established in this work may be extended reliably to more general cases, for example to elongated grains and preferred crystallographic texture materials.

## ACKNOWLEDGMENTS

The authors acknowledge Dr. A. Van Pamel for his valuable suggestions on studying SBCs. M.H. was supported by the Chinese Scholarship Council and the Beijing Institute of Aeronautical Materials. P.H. was funded by the UK Engineering and Physical Sciences Research Council (EPSRC) Fellowship No. EP/M020207/1. M.J.S.L. was partially sponsored by the EPSRC. G.S. and S.I.R. were partially supported by the AFRL (USA) under Prime Contract No. FA8650-10-D-5210.