How to calculate dispersion relation from a Finite Difference (FD) wave simulation

Question

I have a python code that calculates the solution of the inhomogeneous acoustic wave equation for a 2D medium with any velocity and source configuration. It was implemented using Finite Differences from the following set of equations:

$$ \rho( \mathbf{r}) \mathbf{\nabla} \cdot \left( \frac{1}{\rho (\mathbf{r})} \mathbf{\nabla} p(\mathbf{r},t) \right) - \frac{1}{c^2(\mathbf{r})} \frac{\partial^2 p(\mathbf{r},t)}{\partial^2 t} = -s(\mathbf{r},t) $$

where,

$$ s(\mathbf{r},t) = \rho( \mathbf{r}) \frac{\partial^2 i_v(\mathbf{r},t)}{\partial^2 t} $$

and with the propagation velocity

$$ c(\mathbf{r}) = \sqrt{\frac{\kappa(\mathbf{r})}{\rho(\mathbf{r})}} $$

Where $\mathbf{\kappa}$ is the adiabatic compression modulus of the medium, $i_v$ is the source and $\rho$ is the density of the medium.

The FD schema uses 2nd order discretization in time and 4th in space and is implemented for 2D space.

My question is how do I calculate from my numerical simulation the dispersion relation of my code? In fact i want to calculate the phase velocity $v(k) = \frac{\omega(k)}{k}$ dispersion.

I want to do that to compare with the expected relation from literature. I expect it varying with different grid angles, source frequencies and stability parameters. I know how to input those in my simulation but I dont know how to use the results (2D time panels) to establish the dispersion relation.

Answer 1

I'm trying to answer it with physical intuition and arguments.

I assume you'll get $p(\mathbf{r},t)$ as solution. Take only one dimension for example i.e. writing as $p(x,t)=\Sigma_n \Sigma_m a_{mn} Sin(k_mx-\omega_nt)$. If you Fourier decompose it in space first, you'll $k_m$s. Then for each $k_m$, do Fourier transform in time to get $\omega_n$. Now plot of $\omega_n$ vs $a_{mn}$ for one n at a time. The peak will correspond to the leading $k_m$ for that $\omega_n$ (Here I'll argue that only leading order matters; this can be refined to find second and later order dependence on $k_m$). Think about this way-- if some $a_{mn}$ is numerically zero for some m (call $m_0$) for one particular n (call $n_o$), this means that the component of wave, $Sin(k_{m_0}x-\omega_{n_{0}}t)$ doen't exist in the full wave. Similarly, for other n. This way you'll end up with the pairs $(n,m)$. In general, you'll get different m for each n. Finally, find a fit to this curve or use the plot to define your dispersion relation.

Please look into care and considerations of discrete fourier transform in space and time which has some limits on your choice of wavelengths ($k_m$) and time period ($\omega_n$) by cell size and time step.

Sweep long $x, y$ and $z$ separately to get $(k_x,k_y,k_z)$ in each cell to expand this to three dimensions.