ElasticMatrix: A MATLAB Toolbox for Anisotropic Elastic Wave Propagation in Layered Media

Simulating the propagation of elastic waves in multi-layered media has many applications. A common approach is to use matrix methods where the elastic wave-field within each material layer is represented by a sum of partial-waves along with boundary conditions imposed at each interface. While these methods are well-known, coding the required matrix formation, inversion, and analysis for general multi-layered systems is non-trivial and time-consuming. Here, a new open-source toolbox called ElasticMatrix is described which solves the problem of acoustic and elastic wave propagation in multi-layered media for isotropic and transverse-isotropic materials where the wave propagation occurs in a material plane of symmetry. The toolbox is implemented in MATLAB using an object oriented programming framework and is designed to be easy to use and extend. Methods are provided for calculating and plotting dispersion curves, displacement and stress fields, reflection and transmission coefficients, and slowness profiles.


Motivation and Significance
Matrix models of wave propagation in multi-layered elastic solids have had a significant contribution to research areas such as acoustics, geophysics and electromagnetics. A few examples include: structural health monitoring [1], characterisation of interface bonding [2], detection of debonding in joints [3], measuring material properties [4], designing composite layered structures [5], mode sorting of guided waves [6], the physical interpretation of guided wave structures [7], modelling the directional response of Fabry-Pérot ultrasound sensors [8], reflection and transmission of plane waves [9], elastography of layered soft tissues [10], and ice detection on wind turbines [11].
Matrix methods, in particular the partial-wave and global matrix method, represent the stress and displacement fields as a sum of partial-waves for each material of the layered-structure. Each partial-wave represents an upward or downward travelling (quasi-)compressional or (quasi-)shear wave. By invoking boundary conditions at the interfaces of adjacent layers, the partial-wave amplitudes and field properties of the first layer can be related to the last in the form of a 'global' matrix. The resulting matrix equation can be used in two different ways. Firstly, the roots of the equation can be found which give the modal solutions or dispersion curves. Secondly, a subset of partial-wave amplitudes can be defined and the remaining amplitudes solved for. This can be used to calculate the displacement and stress fields within the multi-layered structure when a plane wave is incident. This method will be discussed further in Section 2.
Despite its usefulness, there are few available implementations of the partial-wave method. The current state-of-theart implementation is Disperse [12]. This software has been in development since 1990 and is primarily focused on calculating the dispersion solutions for multi-layered structures. The Disperse software was originally based the partial-wave method, however, is currently being updated to use the spectral collocation method [12,13,14,15,16]. The main limitation with Disperse is that it is closed-source. For this reason it is not easily adaptable for applications arXiv:1909.12740v1 [physics.app-ph] 27 Sep 2019 Figure 1: Diagram for an n-layered elastic medium. In the 2D plane there are four partial-waves with amplitude B n i , these represent (quasi-)compressional (solid arrows) and (quasi-)shear (dashed arrows) waves travelling upwards and downwards in each layer.
that are not dispersion analysis, for example, extracting reflection coefficients or slowness profiles. There currently is only a single open-source code modelling the partial-wave method, (LAMB [17]), however, this is limited to modelling only an isotropic plate.
In this paper, a new open-source toolbox called ElasticMatrix is introduced which uses the partial-wave method for multi-layered structures with an arbitrary number of isotropic and transverse-isotropic materials. Where possible, it is validated against existing literature and has been implemented so that it is both easy to use and extend. Some potential uses of this software are: 1) plotting the slowness profiles of materials, 2) determining the reflection and transmission coefficients of multi-layered structures, 3) finding the dispersion curves of multi-layered structures, 4) plotting the displacement and stress fields, 5) extending the toolbox for other applications, for example modelling the directional response of Fabry-Perot ultrasound sensors [8,18].  ElasticMatrix uses the partial-wave method to model wave propagation in multi-layered elastic solids. The method describes elastic plane-wave propagation along a plane of symmetry for n-layers of rigidly bonded transverse-isotropic materials. An example of an isotropic material is glass, where the material properties are the same when measured from every direction. An example of a transverse-isotropic material is a bundle of fibres, where the properties have translational symmetry axially along the fibre, and are isotropic in the plane perpendicular to this. ElasticMatrix can model layered transverse-isotropic materials if they are aligned such that they have rotational symmetry about the axis perpendicular to the plane of each layer and the wave-vector of the propagating wave lies in the plane of symmetry. In this case, the multi-layered structure can be modelled in two-dimensions. This is illustrated in Fig. 1.
Each partial-wave represents the superposition of waves that have been multiply reflected or transmitted at the interfaces between each layer in a steady-state. The polarisation vector and wave-vector of each of these partial-waves can be found from the Christoffel Equation which is described in Section 2.2. The degree of reflection and transmission depends on the boundary conditions at the interfaces and material properties of each layer. The coupled equations that arise from the boundary conditions can be combined into a 'global-matrix' which allows them to be solved simultaneously, which is discussed in Section 2.3. This global matrix approach can be used to tackle various problems in elastic wave propagation. For example, the singularities of the global matrix give the dispersion curves, and by specifying an incident wave, the resulting wave-field throughout the structure can be calculated. More detailed descriptions of the partial-wave and global-matrix method can be found in [4,5,19,20,21,22,23,24].

Wave-vectors and Polarisation
Firstly, the solution for a plane wave propagating in an unbounded medium is derived. This is needed to calculate the polarisation and wave-vectors for each partial-wave component and the process is repeated independently for every layer. The wave-equation for an anisotropic unbounded medium is where the indices i, j, k, l ∈ {1, 2, 3}, x and t are the spatial and temporal variables and Einstein summation notation is used. The variable u i is the displacement in direction i. The elastic properties of each material are described by the density ρ and the stiffness-tensor C ijkl . The stiffness tensor has 81 components which can be reduced to 21 independent coefficients to describe a fully-anisotropic medium [5]. Here, the analysis is restricted to materials that are either isotropic or transverse-isotropic, which reduces the number of independent coefficients further. As described previously, the wave-vectors of the partial-waves lie in a plane of material symmetry, (x 1 , x 3 ). Here ∂/∂x 2 = 0 and the expanded form of the wave-equation Eq. (1) is where Voigt notation has been used to contract the indices of the stiffness-matrix (where 11 → 1, 22 → 2, 33 → 3, 23 → 4, 13 → 5, 12 → 6). A single-frequency plane wave can be written in the form where i ∈ {1, 2, 3}, ω is the circular frequency, α is the ratio of the vertical and horizontal (ζ) wavenumbers, and A i is the polarisation unit vector which describes the direction of displacement relative to the direction of wave propagation.
where the components of the Christoffel matrix (Γ) are The phase velocity ν along the x 1 axis is calculated from the relation ν = ω/ζ. Solving Eq. (4) admits three solutions for α 2 and therefore 6 solutions for α. From here, the notation α q , where q ∈ {1, 2, ..., 6}, will be used to indicate each solution.
It can be seen from Eq. (4) that the plane wave component A 2 is only dependent on Γ 22 , hence displacement occurring in the (x 1 , x 3 ) plane is independent of displacement in x 2 . Four solutions (q = 1, 2, 3, 4) of α q Eq. (4) can be found when These describe upwards and downward travelling quasi-shear-vertical (qSV) and quasi-longitudinal (qL) waves with the displacement restricted to the (x 1 , x 3 ) plane. The remaining two solutions (q = 5, 6) are found from Γ 22 = 0, and correspond to upward-and downward-travelling quasi-shear-horizontal (qSV) waves. The notation A iq will be used to indicate polarisation vector for each solution q. The displacement field can now be written as where B q is the amplitude of each partial-wave. Additionally, the stress field within the unbounded medium can be found using Hooke's law where i, j, k, l ∈ {1, 2, 3}. For a multi-layered medium, the Christoffel equation Eq. (4) is solved independently for every layer to calculate the polarisation vector and wave-vector of each partial-wave. However, the amplitude B q of each partial-wave is solved by invoking the boundary conditions at the interfaces of adjacent layers. This is discussed in the following section.

Boundary Conditions and Partial-wave Amplitudes
As mentioned previously, the wave-vector of the plane waves are constrained to a plane of symmetry of the transverseisotropic material reducing the analysis to two dimensions, (x 1 , x 3 ). The normal and transverse displacement and stress describing (quasi-)longitudinal and (quasi-)shear-vertical waves for a single layer is written in the form where Only the first four solutions of α q are needed as the motion is restricted to two dimensions (x 1 , x 3 ). The left hand vector of Eqs. (8) contains the components of the displacement and stress, and the right hand vector contains the amplitude of the partial-wave components. The product of the matrices in Eq. (8) will be written as a field matrix F. At an interface x 3 = d between material layers in welded contact, the normal transverse stress and displacement must be continuous across the interface. Therefore, the product of the field matrix and wave amplitudes at the interface of one layer is set equal to the field matrix and wave amplitudes of the adjacent layer. This process is repeated for every interface of the layered medium. For n-layers, there are 4(n − 1) boundary conditions and 4n wave amplitudes which can be arranged into a global matrix. For example, for a medium consisting of 4 layers, the global matrix equation is written  Here, F N n is the 4 × 4 field matrix and B n a 4 × 1 vector of partial-wave amplitudes of layer n at interface N . By assigning values to four of the partial-wave amplitudes, Eq. (9), can be rearranged and solved for the remaining partialwave amplitudes. For example, if a compressional wave in the first medium is incident on the layered-structure, the downward-travelling partial-wave amplitude relating to shear (B 1 2 , Figure 1) in the first layer and upward-travelling partial-wave amplitudes relating to compressional (B n 3 ) and shear (B n 4 ) waves in the last (nth) layer are set to zero. Finally, the downward-travelling partial-wave amplitude relating to a compressional wave in the first layer is set to an arbitrary value (B 1 1 = 1). In this case, the solved amplitudes describe the solution for an incident single-frequency plane-wave at an angle θ or wavenumber ζ and frequency f . For a four-layered medium this would be Here, + and − superscripts indicate the upwards and downwards travelling partial-wave amplitudes and their respective columns in the field matrices. For example, F 1+ 1 would be the third and fourth columns of F 1 1 , and B + 1 would be the third and fourth elements of B 1 . In the example described above, the first element of B − 1 is 1 and all the elements of B + 4 and the second element of B − 1 are 0. Additionally, once the remaining wave-amplitudes for each layer are found, Eq. (8) can be used to find the displacement and stress anywhere in the layered structure. Alternatively, the dispersion curves can be extracted from the model by setting the incident wave-amplitudes of the layered structure to zero and finding the frequency-wavenumber pairs in which the resulting left-hand-side matrix of Eq. (10) becomes singular. The algorithm used is described further in Section 2.5

Shear-Horizontal Waves
Shear-horizontal waves propagate independently of (quasi-)shear-vertical and (quasi-)compressional waves and their propagation is analogous to compressional waves in a liquid [25]. Solutions 5 and 6 of α q correspond to shearhorizontal waves and the displacement and shear stress can be written in the form where D 3q = C 44 α q A 2q (iζ).
By setting the properties of a new medium C so that C 11 , C 22 , C 33 = C 44 and setting the remaining coefficients to zero, Eq. (8) reduces to Eq. (11). Hence, the propogation of shear-horizontal waves can be modelled with Eq. (8) by reassigning the relevant components of the stiffness matrix.

Implementation Details
To construct the global matrix, Eq. (9), a field matrix F i must be calculated for each layer. To improve the conditioning of the matrix, rows relating to displacement equations are scaled by ζ the horizontal wavenumber, and rows relating to stress are scaled by ρω 2 . The system matrix is constructed by looping over each interface, calculating the 4 × 4 field matrices above and below each interface and arranging them into a single matrix. This leads to a rectangular matrix which has 4n columns and 4n − 1 rows. Additionally there are 4n partial-wave amplitudes. Four partial-wave amplitudes are defined and the global-matrix in Eq. (9) is rearranged to be square Eq. (10). The resulting equation is solved using the mldivide function in MATLAB. This function solves a system of linear equations using the fastest algorithm based on the matrix structure. However, the global-matrix becomes singular at values of ζ and ω on or close-to dispersion curve solutions.
The computation of dispersion curves follows the algorithm described in [4]. Firstly, the wavenumber parameter is fixed and the determinant of the global-matrix is found over a range of frequencies. Close to dispersion solutions the determinant of the global-matrix tends to zero. Using these as starting points and taking a limit either side, the exact frequency and wavenumber of the dispersive solution is found using a bisection algorithm. ElasticMatrix makes use of MATLAB's fmincon() function for this. These solutions are the starting points for each dispersion curve. To find the second point on each dispersion curve, the fixed value of wavenumber is increased and the search is performed again. The algorithm then uses linear interpolation to estimate the location of the third, fourth and fifth points on the dispersion curve, similarly using a bisection algorithm to find the exact frequency-wavenumber pairs. After five points have been found, a higher-order polynomial interpolation scheme is used to more accurately predict points on the dispersion curve. The algorithm implemented in ElasticMatrix only searches in the real domain of ζ which is a good-estimate for simple plate structures in a vacuum, however, it may be inaccurate for leaky solutions, for example an plate embedded in soil.

Overview
The ElasticMatrix toolbox implements the partial-wave method using an object-oriented framework in MATLAB. This allows the toolbox to be used with either a simple scripting or command line interface, and makes it easy to use and expand. The software is divided into three classes. The first class, Medium, defines the multi-layered geometry and material properties of each layer. The second class ElasticMatrix is initialised by a Medium object. This class contains the partial-wave method implementation and methods for extracting additional details such as dispersion curves and reflection coefficients. By default, all the calculations use 64 bit precision. The final class, FabryPerotSensor, is an example of how numerical models can be built from the ElasticMatrix and Medium objects. This class inherits ElasticMatrix and can be used to model the directional response of a Fabry-Pérot ultrasound sensor. More details can be found in [8,18]. Each class in the toolbox inherits the MATLAB handle class. Consequently, the object does not need to be reassigned when a method is called. The classes and their respective attributes and methods can be seen in Fig. 2. The toolbox is self contained and has been tested with MATLAB 2016a and above. There are three steps to using the toolbox. Firstly, the geometry of the layered medium must be defined. Secondly, the input parameters to the model should be defined, which are generally a range of angles, frequencies or wavenumbers. Finally, the model can be solved and details such as the reflection coefficients and dispersion curves can be extracted. Note, for clarity in the code implementation, the x 1 and x 3 coordinates are referred to as x and z, respectively.

Medium
The Medium class is used to define the material properties and thickness of each layer. The class is initialised by calling the class constructor with input arguments of the material name followed by its thickness. However, the thickness of the first and last layers are semi-infinite and their values should be set with the Inf keyword. The Medium class will automatically set the thickness of the first and last layer to Inf if another value is used. An example is given below. my_medium = Medium('water', Inf, 'blank', 3e-3, 'PVDF', 1e-3, 'glass', Inf); Here, my medium is an object array and every index in the object array corresponds to a different layer in the medium. In the current example, my medium(3) will return a object with the material properties and thickness associated with PVDF. The 'blank' keyword can used for a material which is not predefined. The material properties and names can be set using their respective set functions. User defined materials can be added to the script materialList.m.

Slowness profiles
Slowness profiles are a plot of the inverse-phase velocity of each bulk wave component. They can be used to determine the angles of reflection and transmission between multi-layered media as well as the direction of energy propagation and skew angle [22]. Slowness profiles are found by solving the Christoffel equation, Eq. 4, and only depend on the material properties of each material. The method .calculateSlowness is part of the Medium class, and calls the function which is an implementation of Eqs. (4) and (5). This takes input arguments of the material properties and phasevelocity and returns the polarisation and wave-vectors. The slowness profiles given by this function are plotted in terms of k x /ω vs k z /ω. For an isotropic material, the slowness profiles for each bulk wave are spherical, however, this is not true for an anisotropic material. An example of the slowness profiles for isotropic-glass and transverse-isotropic beryl is shown in Fig. 3. This figure has been reproduced from [26]. The slowness profiles of the (quasi-)longitudinal, (quasi-)shear-vertical and (quasi-)shear-horizontal bulk waves are shown. As glass is an isotropic material, the slowness profiles are spherical and the magnitudes of L,SV and SH when k x /ω = 0 or k z /ω = 0 are equal to the reciprocal of the compressional-and shear-speeds of glass. For the transverse-isotropic case, when k x /ω = 0, the value of qL is equal to ρ/C 33 and qSV is equal to ρ/C 55 . When k z /ω = 0, the value of qL is equal to ρ/C 11 and the value of qSV is equal to ρ/C 55 . These have been checked in the toolbox example script and all are within numerical precision. my_medium = Medium('glass', Inf, 'beryl', Inf); my_medium.calculateSlowness; my_medium.plotSlowness;

ElasticMatrix
The medium class is used to initialise the ElasticMatrix class which runs the partial-wave method over a range of frequencies, wavenumbers, phasespeeds and angles. Two of these must be defined using the .set functions.
The .calculate method is then used to run the partial-wave procedure. If the properties are not set before the .calculate method is called, the model will run over a predefined range of frequencies and angles. The .calculate method constructs, rearranges and solves the global-matrix, Eq. (10), using the function This function takes input arguments of the material properties and the parameters to calculate over (angles, frequencies, wavenumbers). It returns the determinant of the system matrix and the stresses and displacements at the layer interfaces. Each individual field-matrix is calculated using the function which is an implementation of Eq. (8). This takes input arguments of the material properties, the wave-vector components, polarisation components and the phase velocity and returns the field-matrix. The default calculation is to find the partial-wave amplitudes and interface stresses and displacements when there is a single-frequency compressional wave incident on the structure from the first layer. An example is given below for a titanium plate. my_medium = Medium('water', Inf, 'titanium', 1e-3, 'water', Inf); my_model = ElasticMatrix(my_medium); % initialise class my_model.setFrequency(linspace(1e6, 5e6, 100)); my_model.setAngle(linspace(0, 45, 100)); my_model.calculate;

Reflection and Transmission Coefficients
For a plane wave incident at an oblique angle on a multi-layered structure, the reflection and transmission coefficients determine the amplitude of the wave that is reflected and transmitted at each interface. Knowing these coefficients is useful for a number of applications. For example, selecting the appropriate launch angle when coupling energy into particular modes in a wave-guide, or determining the thickness and material properties of matching layers for ultrasonic transducers [22,27].
The angles of refraction at the interfaces between multi-layered media can be found by studying the slowness profiles. However, slowness profiles do not take into account the boundary conditions at the interfaces. Consequently, the magnitude of each of the refracted waves cannot be calculated directly. For a plane wave incident on a multi-layered structure, the magnitude of the reflection and transmission coefficients are found by normalising the partial-wave amplitudes B n i by the incident plane wave amplitude B 1 1 . This is automatically calculated when using the .calculate method.
An example of the reflection and transmission coefficients at a PVDF-aluminium interface is given below and shown in Fig. 4. For a plane compressional wave incident on a PVDF-aluminium interface, there are four resulting refracted waves. These are a reflected R and transmitted T compressional L and shear S wave. The reflection and transmission coefficients have been compared to the analytic solutions for a two-layered elastic-medium from [22] and have an average error less than 1e −15 which is within numerical precision for a 64 bit floating point number. For clarity the analytical solutions have not been plotted but can be seen in the toolbox examples folder. my_medium = Medium('PVDF', Inf, 'aluminium', Inf); my_model = ElasticMatrix(my_medium); my_model.setFrequency(1e6); my_model.setAngle(linspace(0, 90, 90)); my_model.calculate; my_model.plotRTCoefficients; Figure 4: Longitudinal L and shear S reflection R and transmission T coefficients for a PVDF-Aluminium interface.

Dispersion Curves
Dispersion curves describe the modal solutions of the multilayer structure and describe a wave-mode which propagates parallel to the layer-interfaces independently of a bulk wave. As one example, exciting these modes is essential in ultrasonic inspection. Knowledge of the dispersion curves is useful for determining the most appropriate modes to excite and for optimising the inspection process.
As mentioned in Section 2.5, the modal solutions are found when the global matrix becomes singular. The ElasticMatrix software can calculate dispersion curves for simple layered structures (i.e., a plate in a vacuum or water). However, it is not robust for very-leaky cases, for example a plate embedded in soil. For these types of cases either Disperse, or other techniques based on the spectral-collocation method or semi-analytic finite element method are more appropriate [12,13,14,15,16,28]. An example of the dispersion curves for a 1 mm titanium plate in a vacuum is shown in Fig. 5(a). The dispersion curves are plotted on a graph of frequency vs wavenumber and show the first three symmetric S and anti-symmetric A Lamb modes. The results from Disperse are also plotted and have excellent agreement. my_medium = Medium('vacuum', Inf, 'titanium', 0.001, 'vacuum', Inf); my_model = ElasticMatrix(my_medium); my_model.setFrequency(linspace(0.5e6, 5e6, 100)); my_model.calculateDispersionCurves; my_model.plotDispersionCurves; Figure 5: (a) Dispersion curves for a titanium plate in a vacuum. The solid lines are from ElasticMatrix and the points are generated using Disperse [12]. The first three symmetric (S, black) and anti-symmetric (A, blue) have been plotted. (b) Displacement field for an anti-symmetric and symmetric mode shape.

Displacement and Stress Fields
More information about the wave-physics and guided wave structures can be taken from dispersion curves by plotting the displacement and stress fields at different points. In the ElasticMatrix software implementation, the x and z ranges over which to plot the displacement or stress fields must be specified. The .calculateField(...) method returns a structure with the input ranges and field values at each point of the resulting grid. The values of the structure can be plotted independently or given as an argument to the .plotField method. An example is given below for the displacement field within an titanium plate for a symmetric and anti-symmetric mode. The resulting plot can be seen in Fig. 5(b). field_values = myModel.calculateField(freq, angle, {x_range, z_range}); myModel.plotField(field_values, plot_style);

Impact and Conclusions
This paper introduces a new open-source toolbox called ElasticMatrix which models elastic wave propagation in multi-layered media with anisotropic materials with isotropic or transverse-isotropic symmetry. The toolbox uses the partial-wave method which allows the calculation of slowness profiles, reflection and transmission coefficients, dispersion curves and stress and displacement fields. The software has been implemented using the object-oriented capabilities of MATLAB allowing for a simple command line or scripting interface. The implementation allows researchers to add functionality and integrate the software into other projects. For example, this toolbox has already been used to model the directional response of Fabry-Pérot ultrasound sensors [8]. It is anticipated the research userbase will actively contribute to ElasticMatrix and add to the functionality.