%matplotlib inline
%config InlineBackend.figure_format = 'svg'
import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import widgets, interact
from exact_solvers import acoustics_heterogeneous, acoustics_heterogeneous_demos
from utils import riemann_tools
import seaborn as sns
sns.set_style('white',{'legend.frameon':'True'});
We would like to model acoustic waves propagation through heterogeneous media, like non-homogeneous materials, layered media or any kind of interface, like walls. The materials in acoustic equations are modeled by the density and bulk modulus coefficients. The natural generalization of constant-coefficient acoustic equations is to have spatially dependent coefficients, so the equations take the form $q_t+A(x)q_x=0$. The explicit form is given by
\begin{align*} \left[ \begin{array}{c} p \\ u \end{array} \right]_t
with $p$ and $u$ the pressure and velocity and $\rho(x)$ and $K(x)$ the spatially dependent density and bulk modulus of compressibility. Note this equation is in non-conservative form. An acoustics equation for heterogeneous media in conservative form can be derived in terms of the momentum and strain. However, that is a particular case of the elasticity equations, which are explored in detail in other sections of this book. Furthermore, it is convenient to use the pressure and velocity as variables since they are continuous at interfaces between materials, and they are also more physically intuitively. It is also important to recognize that more complicated systems emerging in applications might not be written in conservation form. Therefore, studying these problems might provide insight and algorithms on how to solve more complicated cases.
We proceed to do the usual analysis. The eigenvalues of the coefficient matrix $A$ are $\pm c(x)$, and the matrix of column eigenvalues is
\begin{align*} R(x) = \left[ \begin{array}{ccccc} -Z(x) & Z(x) \\ 1 & 1 \\ \end{array} \right]. \end{align*}
where
\begin{align*} c(x) &= \sqrt{\frac{K(x)}{\rho(x)}}, \\ Z(x) &= \rho(x) c(x) = \sqrt{K(x)\rho(x)} \end{align*}
are the spatially dependent sound speed and impedance.
In order to solve these equations, we need to do a numerical discretization. Following the spirit of finite volume methods, we approximate $\rho(x)$ and $K(x)$ by piecewise constant functions that are constant in any given grid cell. Therefore, if we can solve the Riemann problem across two given cells, we can extrapolate the solution to the whole grid using standard finite volume method techniques, see (LeVeque 2002). The Riemann problem to solve consists of the acoustic equations Riemann problem with discontinuous initial data and coefficients (discontinuous $\rho$ and $K$). Once this problem is solved it can be used to approximate a continuous density varying material or other similar examples. As the value of $\rho$ and $K$ is different on the left side than on the right side, the eigenvalues and eigenvectors are as follows. The eigenvalues will be given by the sound speed in each of the two mediums,
\begin{align} s_l = -c_l \ \ \ \ \ s_r = c_r \ \ \ \ \ \mathrm{with:} \ \ \ \ \ c_i = \sqrt{\frac{K_{i}}{\rho_{i}}}, \label{eq:achetero} \end{align}
and the eigenvalues by the impedances of each medium as well, so we can write the matrix of column eigenvectors $R=[r_1, r_2]$ as,
\begin{align*} R = \left[ \begin{array}{ccccc} -Z_{l} & Z_{r} \\ 1 & 1 \\ \end{array} \right]. \end{align*}
Once again, we only need to solve $\mathbf{R} \bar{\alpha} = \Delta \bar{q}$, which yields the values of $\alpha$
\begin{align*} \alpha_1 = \frac{-\Delta p + Z_r\Delta u}{Z_l + Z_r}, \ \ \ \ \ \ \alpha_2 = \frac{\Delta p + Z_l\Delta u}{Z_l + Z_r}. \end{align*}
The middle state is again given simply by $q_m = q_\ell + \alpha_1 r_1 = q_r - \alpha_2 r_2$.
This interactive plots allows you to change all of the parameters, as well as the left and right density and bulk modulus so their influence in the phase plane solution can be obersved.
# Initial states [pressure, velocity]
ql = [25.0, 15.0]
qr = [10.0,-15.0]
# Acoustic eq. parameters [rho, bulk(K)]
paramsl = [1.0, 0.5]
paramsr = [5.0, 3.0]
acoustics_heterogeneous_demos.interactive_phase_plane(ql,qr,paramsl,paramsr)
We will show some examples of where this Riemann problem becomes relevant. As in the previous case, we will begin by defining a function to do the interactive plotting for the different cases.
We repeat the shock tube problem for acoustics but now with two materials. The material properties in the acoustic equations are defined by the density and bulk modulus. Therefore, we can solve the acoustics Riemann problem for two materials by simply choosing different densities and bulk modulus and on the left and on the right. This Riemann problem is fundamental to model acoustic wave propagation across interfaces or heterogenous materials. Note the symmetry of the wave speeds is lost since the eigenvalues are the sound speed and the sound speed depends on the material, i.e. on the density and bulk modulus, like shown in equation (\ref{eq:achetero}). Also note the characteristics bend when crossing the origin. This is a consequence of having differente materials on the left and right sides since different materials will yield different sound speeds.
ql = np.array([5,0])
qr = np.array([1,0])
rhol, rhor = 1.0, 20.0 # left and right density
bulkl, bulkr = 4.0, 15.0 # left and right bulk modulus
auxl = [rhol, bulkl]
auxr = [rhor, bulkr]
interact(acoustics_heterogeneous.riemann_plot_func(ql,qr,auxl, auxr), t=widgets.FloatSlider(value=0.0,min=0,max=1.0),
which_char=widgets.Dropdown(options=[None,1,2],description='Show characteristics'));
The solution in the phase plane is
acoustics_heterogeneous_demos.phase_plane_plot()(ql[0],ql[1],qr[0],qr[1],rhol,rhor,bulkl,bulkr,ymin=-2,ymax=2)
In a previous example, we showed the flow into a wall, which basically models the wall as a completely reflective surface. In most cases, this is a good approximation for the reflected waves; however, we could also ask what is the propagated acoustic wave through the wall. We can answer this question by using the air's bulk modulus and density in the right side and the wall's density and bulk modulus on the right. Air actually has density of $\rho \approx 1 kg/m^3$ and $K\approx 100 kPa$, steel on the other hand has $\rho\approx 8000 kg/m^3$ and $K=160 GPa$. Considering the atmospheric pressure to be $p_{atm} = 101325 Pa$, and an acoustic wave hitting the steel at $340 m/s$, we have all the parameters. As expected you will notice the acoustic wave on the steel propagates extremely faster than in the air, which is around $5000 m/s$, around 14 times faster than in air.
patm = 101325.0
ql = np.array([patm,340])
qr = np.array([patm,0])
rhol, rhor = 1.0, 8000.0 # left and right density
bulkl, bulkr = 100000.0, 160000000000.0 # left and right bulk modulus
auxl = [rhol, bulkl]
auxr = [rhor, bulkr]
interact(acoustics_heterogeneous.riemann_plot_func(ql,qr,auxl, auxr), t=widgets.FloatSlider(value=0.0,min=0,max=1.0),
which_char=widgets.Dropdown(options=[None,1,2],description='Show characteristics'));