Elastic Wave Propagation in a Stainless-Steel Standard and Veriﬁcation of a COMSOL Multiphysics Numerical Elastic Wave Toolbox

: Laboratory-based elastic wave measurements are commonly used to quantify the seismic properties of Earth’s crust and upper mantle. Different types of laboratory apparatuses are available for such measurements, simulating seismic properties at different pressure and temperature. To complement such laboratory measurements, we present a numerical toolbox to investigate the seismic properties of rock samples. The numerical model is benchmarked against experimental results from a multi-anvil apparatus, using measurements of a stainless steel calibration standard. Measured values of the mean compressional- and shear-wave velocities at room conditions of the steel block were 6.03 km/s and 3.26 km/s, respectively. Calculated numerical results predicted 6.12 km/s and 3.30 km/s for compressional and shear-wave velocities. Subsequently, we measured Vp and Vs up to 600 MPa hydrostatic conﬁning pressure and 600 ◦ C. These measurements, at pressure and temperature, were then used as the basis to predict numerical wave speeds. There is, in general, good agreement between measurement and predicted numerical results. The numerical method presented in this study serves as a ﬂexible toolbox, allowing for the easy setup of different model geometries and composite materials.

As such, laboratory-derived elastic properties are invaluable for interpreting seismic data and provide an important method for investigating tectonic and geodynamic processes and understanding the inner Earth structure [7,24,25]. A complement to laboratory measurements are provided through predictive (e.g., effective medium) and numerical models e.g., [26][27][28][29][30][31][32][33][34]. Predictive and numerical experiments can help understand elastic wave propagation in a material by providing additional information to laboratory measurements. This is because models provide the possibility to investigate the effects of composition, structure, and geometry using inferences of the compositional e.g., [31,35] and structural characterization of materials e.g., [18][19][20][21][32][33][34]36]. We note that laboratory and numerical experiments are complementary and best used in combination.
In this study, the compressional (Vp) and shear wave (Vs) velocities of a near isotropic standard stainless steel were measured up to 600 MPa confining pressure and temperatures up to 600 • C using an ultrasonic frequency measurement setup in a multi-anvil apparatus. A numerical model created with the COMSOL Multiphysics ® software was set up and benchmarked against the laboratory measurement data. COMSOL Multiphysics is a general-purpose commercial finite element software. Models were set up and tested in 2D and 3D to simulate elastic wave propagation in the multi-anvil apparatus for the same stainless steel standard in order to verify the model based on laboratory data. This numerical setup can be used as a toolbox to obtain a better understanding of the role of different features, ranging, for example, from microstructural effects of minerals in rocks to compositional layering at the macroscopic scale. In this study, the numerical benchmark models consider an elastic isotropic material with various boundary conditions.
In this study, the compressional (Vp) and shear wave (Vs) velocities of a near isotropic standard stainless steel were measured up to 600 MPa confining pressure and temperatures up to 600 °C using an ultrasonic frequency measurement setup in a multi-anvil apparatus. A numerical model created with the COMSOL Multiphysics ® software was set up and benchmarked against the laboratory measurement data. COMSOL Multiphysics is a general-purpose commercial finite element software. Models were set up and tested in 2D and 3D to simulate elastic wave propagation in the multi-anvil apparatus for the same stainless steel standard in order to verify the model based on laboratory data. This numerical setup can be used as a toolbox to obtain a better understanding of the role of different features, ranging, for example, from microstructural effects of minerals in rocks to compositional layering at the macroscopic scale. In this study, the numerical benchmark models consider an elastic isotropic material with various boundary conditions.

Experimental Setup
The laboratory measurements were performed with a triaxial multi-anvil apparatus ( Figure 2) at the Institute for Geosciences at Kiel University, Germany. The apparatus has six pistons that are used to apply stresses. Ultrasonic velocities (P-and S-waves) were measured using the ultrasonic pulse transmission technique e.g., [1,2,6,[8][9][10][11][12][13][14][15]36]. Three ultrasonic transducers were mounted at the end of each piston in a configuration that allows for the transmission of one compressional wave velocity, Vp, and two orthogonally polarized shear wave velocities, Vs 1 and Vs 2 , along the three sample axes. The dominant frequencies of the transducers were 2 MHz and 1 MHz for P-and S-waves, respectively.

Experimental Setup
The laboratory measurements were performed with a triaxial multi-anvil apparatus ( Figure 2) at the Institute for Geosciences at Kiel University, Germany. The apparatus has six pistons that are used to apply stresses. Ultrasonic velocities (P-and S-waves) were measured using the ultrasonic pulse transmission technique e.g., [1,2,6,[8][9][10][11][12][13][14][15]36]. Three ultrasonic transducers were mounted at the end of each piston in a configuration that allows for the transmission of one compressional wave velocity, Vp, and two orthogonally polarized shear wave velocities, Vs1 and Vs2, along the three sample axes. The dominant frequencies of the transducers were 2 MHz and 1 MHz for P-and S-waves, respectively. The apparatus is capable of operating at temperatures up to 750 ℃; however, in order to prevent the initiation of thermal cracks or the propagation of existing cracks with thermal-induced stresses, as a laboratory routine, the samples are generally tested to a maximum confining pressure of 600 MPa and a temperature up to 600 °C. Changes in sample dimensions, as a function of pressure and temperature, are monitored by linear strain sensors, which are used to correct Vp and vs. measurements. In the current study, the experiments began at room pressure and temperature conditions, followed by an incremental increase in pressure to 600 MPa. In the second phase, the confining pressure was kept constant, and the temperature was increased to 600 °C, in steps of 100 °C. During these two phases, when the apparatus reached the desired condition, ultrasonic wave pulses were applied to the sample along the three principal axes. V and V of the corresponding nine ultrasonic waves were measured along the three sample axes, three were P, and six were S waves.
The elastic wave speeds were calculated by the following equations: The apparatus is capable of operating at temperatures up to 750°C; however, in order to prevent the initiation of thermal cracks or the propagation of existing cracks with thermal-induced stresses, as a laboratory routine, the samples are generally tested to a maximum confining pressure of 600 MPa and a temperature up to 600 • C. Changes in sample dimensions, as a function of pressure and temperature, are monitored by linear strain sensors, which are used to correct Vp and Vs measurements. In the current study, the experiments began at room pressure and temperature conditions, followed by an incremental increase in pressure to 600 MPa. In the second phase, the confining pressure was kept constant, and the temperature was increased to 600 • C, in steps of 100 • C. During these two phases, when the apparatus reached the desired condition, ultrasonic wave pulses were applied to the sample along the three principal axes. V p and V s of the corresponding nine ultrasonic waves were measured along the three sample axes, three were P, and six were S waves.
The elastic wave speeds were calculated by the following equations: where L is the sample length (in m), ∆t is the time of flight of the ultrasonic wave, measured in seconds. To compensate for the change in length at higher pressure and temperature conditions, the linear strain was used to recalculate the velocity, using: where ∆L represents the change in length of the sample, measured based on the piston displacement at elevated pressures and temperatures. Similarly, the density (ρ) of the sample is calculated based on the experimental condition, which at room conditions is: where m is the sample mass (kg) and V is the sample volume (m 3 ). Mass and volume were determined by measurements of weight and sample length, and since the sample is small, any influence due to gravity is neglected. At elevated pressure and temperature conditions, the change in volume (∆V) is considered, and: The first and second dynamic elastic moduli (Lamé's constants, λ and µ) are calculated from the measured Vp and Vs, such that: from which Young's modulus (E) and the Poisson ratio (υ) can be calculated, according to: Finally, we also calculate the mean Vp and Vs, as well as the anisotropy of Vp (%AVp), using: where Vp max , Vp min, and Vp int are the maximum, minimum, and intermediate P-wave velocities measured along the principal sample axes ( Figure 2c). In addition, we calculated the difference in shear wave speed of the two polarized shear waves in terms of percentage (dVs%), using dVs (%) = (Vs1 + Vs2)/(Vs1 − Vs2) × 100 (%). Measurements were carried out in two phases for increasing pressure and temperature. In the first phase of measurements, the ultrasonic wave velocities (V p , V s ) were measured while the pressure was increased step-wise from the room conditions up to 600 MPa (Figure 3a,c,e,g). For each pressure increment, the sample was allowed to equilibrate at the new pressure for 5 min before the measurement was conducted. In the second phase, the confining pressure was kept constant at 600 MPa, and the temperature was increased step-wise to 600 • C (Figure 3b). A period of 30 min was used to let the sample equilibrate at the new temperature. The confining pressure was kept near-constant, at 600 MPa, to limit thermally induced stresses due to expansion and re-opening of pores, flaws, and micro-cracks in the sample while measuring the ultrasonic wave speeds (Figure 3b,d,f,h). Linear strains from all 3 sample axes were measured during confining pressure conditions from room pressure to 600 MPa and temperatures from room conditions up to 600 • C (Figure 4a,b). Linear strain measurements reflect sample volume changes because of increasing pressure and temperature. This volumetric change, calculated from the linear strain, indicates how the sample changes dimensions during pressurization and increase in temperature. The linear strain is small for the pressure increase and relatively much larger for an increase in temperature. In other words, the effect of temperature is significantly greater than the effect of pressurization in influencing the sample volume (Figure 4a,b). Figure 4c,d show the calculated Young's modulus and Poisson ratio, from the Vp and Vs measurements, as a function of pressure and temperature. For further details on the laboratory measurement technique, we refer to [6,36].

Governing Equations and Constitutive Laws
For a general solid in two and three dimensions, with boundary conditions in one to two dimensions and neglecting body force, the equation of motion can be expressed in Einstein notation as: where σij is the stress tensor, ρ is density, and u is the displacement vector. Dots above variables specify time derivatives. This equation is supplemented with suitable sets of boundary and initial conditions. For rheology applications, a pure linear elastic constitutive law (Hooke's law) is used [37]: In Equation (13), the isotropic form of Hooke's law is given in terms of E and v; εij is the strain tensor, εαα is the volumetric strain, and δij is the Kronecker delta. The ε are also the infinitesimal strain tensor components which are kinematically related to the displacement field u , as follows:

Numerical Methodology Governing Equations and Constitutive Laws
For a general solid in two and three dimensions, with boundary conditions in one to two dimensions and neglecting body force, the equation of motion can be expressed in Einstein notation as: where σ ij is the stress tensor, ρ is density, and u i is the displacement vector. Dots above variables specify time derivatives. This equation is supplemented with suitable sets of boundary and initial conditions. For rheology applications, a pure linear elastic constitutive law (Hooke's law) is used [37]: In Equation (13), the isotropic form of Hooke's law is given in terms of E and v; ε ij is the strain tensor, ε αα is the volumetric strain, and δ ij is the Kronecker delta. The ε ij are also the infinitesimal strain tensor components which are kinematically related to the displacement field u i , as follows: In COMSOL Multiphysics software, version 6.0, the Solid Mechanics module part of the Structural Mechanics Module has implemented the weak form of the above sets of equations in the framework of the finite element method and can be efficiently used to describe the seismic wave propagation in a solid media that is modeled as an isotropic linear elastic media.
The 3D FE mesh used was tetrahedral. The model steps are given in microseconds due to the very high frequency (MHz) of the incident wave, and the size of the samples is given in meters ( Figure 5). The edge boundary condition for the compressional wave is zero displacement in the normal direction and free slip in the direction parallel to the boundary. Symmetry and roller boundary conditions are both the same in two-dimensional. For the shear wave, the displacements perpendicular to the polarization direction are set to zero, whereas in the polarization direction, the corresponding stress is set to zero (that is, on the plane with normal motions, the normal stress is set to zero, and on the plane with shear motions, the associated shear stress is set to zero). The buffer zone (piston setup in the current model) is formed out of the same material as the steel samples tested for this bench-mark study. Table 1. Young's modulus, Poisson's ratio, and density for the laboratory measurements (used as input for the numerical model) and the output results in the form of the P wave and a shear wave (note that no polarization is reported in the numerical model).   Table 1.   Table 1.

Laboratory Measurements of Ultrasonic Wave Velocities and Dynamic Elastic Properties
From the Vp and Vs measurements along sample X, Y and Z-axes, it was possible to obtain Vp mean and Vs mean (Equation (9)), as well as AVp, AVs (Equation (4)), and linear strain. Vp mean increased from 5.97 km/s to 6.08 km/s, from room pressure conditions to 600 MPa. The values of Vp mean decreased from 6.08 to 5.80 km/s during increasing temperature to 600 • C (Figure 3b). Anisotropy of Vp is below 1%, indicating that the steel can be considered effectively isotropic. The anisotropy increases slightly from temperatures of 500 • C to 600 • C but is still below 1%. The intrinsic Vp and Vs presented in Figure 3 represent linear fits to the mean Vp above 200 MPa for Vp and Vs as a function of pressure, whereas the intrinsic Vp and Vs as a function of temperature take into account all Vp and Vs measurements as a function of temperature.
Linear strain along the 3 sample axes (ε1, ε2, and ε3) was determined from the piston displacement and shows how the sample deforms as a function of pressure and temperature. Negative linear strain is considered as an expansion of the sample in the laboratory measurement, whereas positive strain is considered compaction of the sample. Importantly, these data provide changes in sample length, which are subsequently used to calculate Vp and Vs and recalculate density. Linear strain appears similar along all sample axes during pressurization, with somewhat higher positive values of ε2 (sample shortening along the y-axis), compared to ε1 and ε3. During increasing temperature, linear strain along all sample axes is very similar, with negative values (expansion) but higher relative values compared to linear strain during pressurization. This indicates that heating the steel has a larger influence on sample volume compared to the pressurization, and hence also on Vp and Vs. From the Vp, Vs measurements, the material properties, including Young's modulus, Poisson's ratio, and density, were calculated (Figure 4c,d). The elastic constants and density served as input parameters for the numerical modeling of Vp and Vs and benchmarking of the model.

Numerical Modelling Results
The numerical models (used for benchmarking) were developed in two and three dimensions with the same material properties (Table 1) and similar geometry and boundary/initial conditions. The three-dimensional model is drawn to simulate a realistic setup of pistons and sample with similar dimensions used as in the laboratory ( Figure 5). Since the geometry is symmetric, the extension from the two-dimensional case to the threedimensional case is considered trivial, although it could be affected if different boundary conditions were used in the two different model setups. The geometry of the model, including the steel sample (in orange) in the middle and the pistons of the multi-anvil apparatus (in grey), are shown in Figure 5. Two buffer zones, representing the apparatus pistons, were placed on the left and right sides of the specimen to resemble the laboratory conditions as close as possible and eliminate model artifacts and edge effects as much as possible. Comparing the first peak arrival time of the incoming and outgoing waves of the constant geometry (in the sample), it is possible to calculate the Vp and Vs using Equation (1) from ( Figure 6). The numerical test specimen was set up with material properties from Table 1. The modeling results for two-and three-dimensional cases were essentially identical and differed by less than 2% in general compared to the laboratory measurements under all experimental conditions (Table 2). Even though it is expected that the measured and modeled Vp and vs. should perfectly coincide (given that the measured values are used as input for the models), the small difference likely arises from uncertainties in the picking of the modeled velocities and effects related to the setup of the model (e.g., size of the model mesh). In this model, the polarization is imposed as in Figure 5c,d. The polarized shear waves are equal for all two-dimensional cases, and no difference is observed in the S-waves for isotropic three-dimensional cases. Table 2. Laboratory measurements, numerical modeling, and their Vp data. Abbreviations: Temptemperature (°C); Pc-Confining pressure; VpL-laboratory-measured Vp; VpN-numerically modeled Vp. Vpmean are calculated from laboratory measurements using Equation (9).
Pc  The modeling results for two-and three-dimensional cases were essentially identical and differed by less than 2% in general compared to the laboratory measurements under all experimental conditions (Table 2). Even though it is expected that the measured and modeled Vp and Vs should perfectly coincide (given that the measured values are used as input for the models), the small difference likely arises from uncertainties in the picking of the modeled velocities and effects related to the setup of the model (e.g., size of the model mesh). In this model, the polarization is imposed as in Figure 5c,d. The polarized shear waves are equal for all two-dimensional cases, and no difference is observed in the S-waves for isotropic three-dimensional cases.

km/s) (km/s) (km/s) (km/s) (km/s) (km/s) (km/s) (km/s) (km/s) (km/s)
The density and elastic moduli of steel measurements make it possible to consider the influence of pressure and temperature effects on Vp and Vs in the numerical model (Figures 3 and 4). The changes in density and elastic moduli during applied confining pressure were small, with~2% change in density and~5% change in elastic moduli from room pressure to 600 MPa. However, these changes were much more significant for increasing temperature (at a fixed confining pressure of 600 MPa), where density and Young's modulus changed from 8256 kg/m 3 and 228 GPa at room temperature to 8053 kg/m 3 and 188 GPa at 600 • C (Figure 4c,d). Poisson's ratio changed from 0.298 to 0.320 for the same increase in temperature. In the model, we use the recalculated density and elastic moduli to model Vp and Vs in order to incorporate the effect of pressure and temperature on elastic wave velocities. The numerical modeling results, presented in Figure 7 and in Table 2, generally illustrate no difference compared to the measured Vp and Vs. The models were ran as an isotropic elastic medium and therefore the models currently run very fast. However, more complicated models in three-dimensional geometries will have much longer run times than two-dimensional ones.

Discussion and Conclusions
This study presents a new COMSOL Multiphysics-based numerical model for elastic wave propagation in materials in two-and three-dimensions, which is verified (benchmarked) using the results from laboratory measurements of steel up to 600 MPa confining pressure and temperatures up to 600 • C. The numerical modeling results show less than a 2% difference with measured experimental values; there was no difference between results from two-and three-dimensional models.
Although the model presented in this study presents a relatively simple application, in which the velocities obtained are based on known laboratory-measured Vp and Vs as a function of pressure and temperature, there is considerable flexibility in setting up different kinds of models. For example, geometries with anisotropic minerals due to lattice preferred orientation, stratified geometry, grain size and shape and structural effects, such as inclusions of partial melting, the existence of microcracks, the effect of grain boundaries, and grain sliding. The advantage of the numerical model, developed in COMSOL Multiphysics, is the possibility to investigate elastic wave propagation in a variety of model setups, considering numerical experiments both on a laboratory scale and field scale (as well as upscaling of laboratory results). This numerical tool can include a variety of mineral arrangements (geometry) and can propagate P-and S-waves from different directions to determine seismic anisotropy; it is simple to extend investigations to study the directional dependency of waves and anisotropy into three dimensions. Three-dimensional models are valuable since all laboratory experiments are performed in three dimensions.
A notable effect in laboratory measurements is the closure of pores and microcracks during sample pressurization, which is illustrated by an exponential increase in elastic wave speed at low confining pressure, which generally occurs below 200 MPa e.g., [1,2,[38][39][40]. In geological samples, this is known as the crack-closure domain, in which no notable contribution to sample volume change occurs (i.e., no notable linear strain changes), but with a pronounced effect on Vp and Vs. To incorporate such crack closure requires a modification of the numerical model but remains possible with the COMSOL Multiphysics software. In contrast, the effect of temperature on samples, which results in thermal expansion, is evident in linear strain measurements and associated sample volume change. An extension of the model can be set up to investigate the effects of pressure and temperature changes on Vp and Vs by incorporating thermal expansion and pressure/temperature-dependent derivatives of Vp and Vs.
Current numerical results allow studying the relations between the structural framework of the rocks (foliation, lineation), velocity anisotropy, shear wave splitting, and shear wave polarization. Anisotropic materials can be considered in the model using the complete elastic stiffness tensor (as defined by Hooke's law in anisotropic form) with up to 21 independent elastic constants e.g., [29] (see appendix/online Supplementary Materials). Such a model setup can play an essential role when studying the influence of crystallographic preferred orientation (CPO) developed by viscoplastic deformation in ductile materials e.g., [16,17,[24][25][26][27][28]. Data sets for CPOs can be obtained from measurements using Scanning Electron Microscope Electron Backscatter Diffraction or, for example, from viscoplastic self-consistent models e.g., [41,42]. Changes in the shape and sizes of geometrical inclusions can, for example, represent grain shape and size e.g., [6,36,43,44]), as such numerical models are used to further investigate the effect of grain sizes and grain boundary on ultrasonic waves. The numerical COMSOL Multiphysics model introduced in this study provides a relatively simple option to study the elastic properties of complex geological materials and serves as a useful addition in the fields of digital rock physics and materials science. Further information and examples of applications with the COMSOL Multiphysics model presented in this study are available in Supplementary Materials and the YouTube channel: https://youtu.be/XxCVrix54V4 (accessed on 20 February 2022).

Supplementary Materials:
The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/resources11050049/s1, S1: Numerical model setup for anisotropic medium; The module recipe used with this COMSOL model is freely available in Supplementary Materials and the YouTube channel https://youtu.be/XxCVrix54V4 (accessed on 20 February 2022)