Contents

- JE5: Hyperbolic balance laws with dispersive source terms
- Linearization and hyperbolic balance laws with dispersive sources
- The Dispersive Euler Equations
- Exact solution for initial step function perturbation
- A note on solving dispersive Euler equations with Lucee
- Problem 1: Case \(\omega_c = 10\)
- Problem 2: Case \(\omega_c = 100\)
- Problem 3: Sod-shock problem
- Conclusions

While solving the *two-fluid Riemann problems* small oscillations are observed in the
solution. On first sight it might be thought that these oscillations
are numerical and do not reflect a correct solution to the ideal
two-fluid equation system. In this note I show that such solutions can
(and do) occur when certain types of source terms are present and the
homogeneous system is hyperbolic.

Linearizing a hyperbolic balance law can lead to important insights into the structure of solutions, although linearization can not reveal the rich non-linear phenomena (shocks, rarefactions and contact discontinuities) that are described by these equations. We start from the non-conservative form of the equations

\[\frac{\partial \mathbf{v}}{\partial t}
+ \mathbf{A}_p\frac{\partial \mathbf{v}}{\partial x} = \mathbf{s}_p\]

where \(\mathbf{v}\) are the primitive variables, \(\mathbf{A}_p\) is a Jacobian matrix and \(\mathbf{s}_p\) are source terms. Linearize about a uniform equilibrium \(\mathbf{v}_0\) to write \(\mathbf{v} = \mathbf{v}_0 + \mathbf{v}_1\), where \(\mathbf{v}_1\) is a small perturbation. Using a Taylor series expansion to first-order to write \(\mathbf{s}_p(\mathbf{v}) = \mathbf{s}_p(\mathbf{v}_0) + \left( {\partial \mathbf{s}_p}/{\partial \mathbf{v}} \right)_{\mathbf{v}_0} \mathbf{v}_1\) and \(\mathbf{A}_p(\mathbf{v}) = \mathbf{A}_p(\mathbf{v}_0) + \left( {\partial \mathbf{A}_p}/{\partial \mathbf{v}} \right)_{\mathbf{v}_0} \mathbf{v}_1\) and letting \(\mathbf{M}_p \equiv {\partial \mathbf{s}_p}/{\partial \mathbf{v}}\), the linear form of the non-conservative equation becomes

\[\frac{\partial \mathbf{v}_1}{\partial t}
+ \mathbf{A}_p(\mathbf{v}_0)\frac{\partial \mathbf{v}_1}{\partial x}
= \mathbf{M}_p(\mathbf{v}_0)\mathbf{v}_1.\]

To understand the mathematical structure of the linear equation a Fourier representation of the solution is assumed and each mode is represented as \(\mathbf{v}_1 = \mathbf{\hat{v}}_1 e^{i\omega t} e^{i k x}\), where \(\omega\) is the frequency and \(k\) is the wavenumber. Using this we obtain

\[\left[
i\omega\mathbf{I} + ik\mathbf{A}_p(\mathbf{v}_0) - \mathbf{M}_p(\mathbf{v}_0)
\right] \mathbf{v}_1 = 0,\]

where \(\mathbf{I}\) is a unit matrix. For non-trivial solutions
the determinant of the matrix in the square brackets must
vanish. Another way to state this condition is that the frequency and
wavenumbers must be related by the *dispersion relations*
\(\omega = \omega(k)\) which are the eigenvalues of the matrix
\(-k\mathbf{A}_p(\mathbf{v}_0) -
i\mathbf{M}_p(\mathbf{v}_0)\). Thus if
\(\lambda^p(k,\mathbf{v}_0)\) is the \(p^{\textrm{th}}\)
eigenvalue of this matrix then the \(p^{\textrm{th}}\) branch of
the dispersion relation is \(\omega = \lambda^p(k,\mathbf{v}_0)\).

If the dispersion relation is purely real then the equation system
will support non-decaying waves. Further, if the dispersion relation
is *linear*, then a wave packet will propagate without distortion. If,
however, if the dispersion relation is non-linear (but still real),
the wave packet will suffer dispersion, i.e. waves with different
wave-numbers will propagate with different group and phase speeds.

For the simple case of vanishing sources (\(\mathbf{s}_p=0\)) the dispersion relation reduces to \(\omega = \lambda^p k\), where \(\lambda^p\) are the eigenvalues of the Jacobian matrix \(\mathbf{A}_p\). As the homogeneous equation is assumed to be hyperbolic the eigenvalues are all real, indicating that for the homogeneous case waves will propagate without dispersion or decay.

This simple analysis indicates that the linear solution will depend on
the nature of the eigenvalues of the source Jacobian matrix
\(\mathbf{M}_p\). In case the eigenvalues of this matrix are
*purely imaginary*, the dispersion relation will be real but waves
will suffer dispersion as they propagate in the background uniform
equilibrium. In this case the system of equations is called
*hyperbolic balance laws with dispersive source terms*. (This is my
own terminology and I have not seen such equations discussed in the
literature).

Several simple hyperbolic balance laws with dispersive source terms can
be constructed. However, a particularly useful system that has
properties similar to the ideal two-fluid equations is the following
*dispersive Euler* system

\[\begin{split}&\frac{\partial \rho}{\partial t} + \nabla\cdot(\rho\mathbf{u}) = 0 \\
&\frac{\partial \mathbf{u}}{\partial t} +
\mathbf{u}\cdot\nabla\mathbf{u} =
-\nabla p/\rho + \lambda\mathbf{u}\times\mathbf{b} \\
&\frac{\partial p}{\partial t} + \mathbf{u}\cdot\nabla p =
-\gamma p \nabla\cdot\mathbf{u}\end{split}\]

where \(\mathbf{b}(\mathbf{x})\) is a time-independent vector field and \(\lambda\) is a constant.

Consider solving the linearized system in one-dimension. Linearizing around stationary uniform initial state \(\rho = \rho_0\), \(p = p_0\) and assuming \(\mathbf{b} = (0,0,b_z)\), where \(b_z\) is constant leads to the linear system

\[\begin{split}\frac{\partial \rho_1}{\partial t}
&= -\rho_0\frac{\partial u_1}{\partial x} \\
\rho_0\frac{\partial u_1}{\partial t} &=
-\frac{\partial p_1}{\partial x} + \rho_0 \lambda v_1 b_z \\
\rho_0\frac{\partial v_1}{\partial t} &= -\rho_0 \lambda u_1 b_z \\
\frac{\partial p_1}{\partial t} &=
-\gamma p_0 \frac{\partial u_1}{\partial x}\end{split}\]

Assuming solutions of the form

\[f(x,t) = \sum_{n=0}^\infty f_n e^{i(k_n x + w_n t)}\]

for \(f\in \{\rho_1,u_1,v_1,p_1\}\) and \(k_n\) is the wave number and \(\omega_n\) is frequency we get the algebraic equations

\[\begin{split}i \omega_n \rho_1 &= - i k_n u_1 \rho_0 \\
i \omega_n u_1 \rho_0 &= -i k_n p_1 + \lambda v_1 b_z \\
i \omega_n v_1 \rho_0 &= - \rho_0 \lambda u_1 b_z \\
i \omega_n p_1 &= i \gamma p_0 k_n u_1.\end{split}\]

From this the dispersion relation can be computed as

\[\omega_n = \pm ( k_n^2 c_{s0}^2 + \omega_c^2 )^{1/2}\]

Here \(c_{s0} \equiv \sqrt{\gamma p_0/\rho_o}\) is the speed of sound and \(\omega_c \equiv \lambda b_z\) is the eigenvalue of the source Jacobian.

Consider a initial perturbation of the form \(u(x,0)\) where

\[u_1(x,t) = U_0 \sum_{n=0}^N
\frac{i}{2n+1} e^{i k_nx} e^{i \omega_n t}\]

with \(k_n = 2\pi(2n+1)\). For \(N\rightarrow \infty\) this represents the propagation of a step function perturbation. Letting \(u_i^{(n)} \equiv i U_0 /(2n+1) e^{i(k_nx+\omega_nt)}\) the Fourier components of the other flow variable perturbations are given by

\[\begin{split}\rho_1^{(n)} &= -\frac{k_n\rho_0}{\omega_n} u_1^{(n)} \\
v_1^{(n)} &= -i\frac{\lambda b_z}{\omega_n} u_1^{(n)} \\
p_1^{(n)} &= -\frac{\gamma k_n p_0}{\omega_n} u_1^{(n)},\end{split}\]

summing which over \(n=0,\ldots,N\) gives the exact solution to the linear problem. The following figure shows the exact solution for \(N=5000\), \(\omega_c = 10\) and \(c_s = \sqrt{2}\) at time 1000.

The dispersive Euler equations can be solved by adding a source term
to the Euler equations. The source terms can be implemented using a
Lorentz force object. This object needs an electric and magnetic field
as input. Hence, we need to allocate memory for all the field
components and set the electric field to zero. Due to the peculiarity
of the point ODE integrator, this memory needs to be part of the fluid
fields. Hence, in the simulations shown below (see, for example,
*s40*) the fields have 11
components (5 for fluids and 3 for electric field and 3 for magnetic
field).

A series of simulations was performed for the case of \(\omega_c = 10\) and \(c_s = \sqrt{2}\). To avoid exciting all the Fourier modes in the step function, the expansion was carried out to only \(N=9\) modes. The solution was computed on grids of 100, 200, 300 and 400 cells. The results of velocity \(u(x,t)\) are shown below at \(t=3\). The wave-propagation scheme has intrinsic diffusion due to which the small wavelength features are poorly resolved when the grid is relatively coarse.

In these simulations, the influence from sources was increased by setting \(\omega_c = 100\). The simulation is run on a grid with 200 cells. The time-step for this case is constrained by the need to resolve the oscillations from the source terms. Taking \(k \approx 1/\Delta x = 1/200\) we get the largest frequency as approximately 283. To resolve this the time-step needs to much smaller than \(1/238 \approx 0.0035\). This forces a more restrictive CFL number (0.5) than allowed by stability of just the hyperbolic part. If the oscillations are not resolved significant phase errors are seen in the solution.

The simulation was next run with 400 cells. This significantly improves the numerical solution even though some small-scale features are still not resolved correctly.

The previous simulations show the effect of dispersive source terms on linear problems. In this simulation I solve the sod-shock problem for the dispersive Euler equations. This is a highly non-linear problem and shows complex shock structure. The problem is initialized with a discontinuity at \(x=0.5\) and with left and right states

\[\begin{split}\left[
\begin{matrix}
\rho_l \\
p_l
\end{matrix}
\right]
=
\left[
\begin{matrix}
3.0 \\
3.0
\end{matrix}
\right],
\qquad
\left[
\begin{matrix}
\rho_r \\
p_r
\end{matrix}
\right]
=
\left[
\begin{matrix}
1.0 \\
1.0
\end{matrix}
\right].\end{split}\]

and is run to \(t=0.1\) on a grid of 800 cells with \(\mathbf{b} = (0.75, 0.0, 1.0)\), \(\lambda=100\) and \(\gamma = 5/3\).

The results are shown below. These show significant differences between the zero-source case and the one with the dispersive sources. Note that the solution looks like the two-fluid solutions to the Riemann problem.

One conclusion from these series of simulations is that dispersive
source terms can cause small-scale features in the solution. To
resolve these features sufficient spatial *and* temporal resolution is
needed. Poor spatial resolution can diffuse the oscillations while
poor temporal resolution can lead to phase errors. In physical
problems (for example multi-fluid plasmas) there is usually some
diffusion that is active at small scales and can be important when the
gradients change rapidly over a few grid cells. This physical
diffusion will wipe out the oscillations from the dispersive
sources. Hence, in such cases the resolution of the oscillations might
not be so important. However, from a mathematical view-point the
numerical schemes need to be accurate enough to resolve such features.