It would be interesting to use the AFM formulation to quickly derive not only the dispersion relationship for N-layers of linear shear, but also develop numerical simulations of traveling waves and relationships between the pressure and other internal wave properties.

## Setting up the equations of motion

Let [mathjax] $$\psi_j$$ represent the stream-function within each fluid layer, and let $$\varphi$$ be some harmonic function that is periodic in $$x$$.

In each fluid layer, we can write the integral identity $\oint_{\mathscr{D}_j}\left(\varphi\frac{d\psi_j}{dn} + \omega \varphi\right)\,ds = 0$

This allows us to relate develop a series of integral relationships in the following form: $\int_{\mathscr{S}_{j+1}}\left(\varphi_z\frac{d\psi_j}{dn_{j+1}} + \omega\varphi\bigg\vert_{\eta_{j+1}}\right)\,dx = \int_{\mathscr{S}_{j+1}}\left(\varphi_z\frac{d\psi_j}{dn_{j}} + \omega\varphi\bigg\vert_{\eta_{j}}\right)\,dx$

Here, since both $$\varphi$$ and the streamfunction $$\psi_j$$ are periodic in $$x$$, we are able to eliminate the contributions from the side boundaries, leaving only contributions from the interface above and below the fluid domain.

Using the pressure condition within each fluid layer and evaluating it at the interface $$z = \eta_{j+1}-d_{j+1}$$, we find equations of the form $\frac{1}{2}\left(1 + \eta_{x,j+1}^2\right)\psi_{z,j}^2 + g\left(\eta_{j+1} – d_{j+1}\right) + \omega_j\psi_j = -p\bigg\vert_{\eta_{j+1} – d_{j+1}} + \frac{1}{2}Q_{j}^2 + \frac{1}{2}c^2$and$\frac{1}{2}\left(1 + \eta_{x,j+1}^2\right)\psi_{z,j+1}^2 + g\left(\eta_{j+1} – d_{j+1}\right) + \omega_{j+1}\psi_{j+1} = -p\bigg\vert_{\eta_{j+1} – d_{j+1}} + \frac{1}{2}Q_{j+1} + \frac{1}{2}c^2$again, both evaluated at $$z = \eta_{j+1}$$_. Since we know that the pressure must be the same across the interface, the above equations reduce to$\frac{1}{2}\left(1 + \eta_{x,j+1}^2\right)\left(\psi_{z,j+1}^2 – \psi_{z,j}^2\right) + \omega_{j+1}\psi_{j+1} – \omega_j\psi_j = \frac{1}{2}\left(Q_{j+1}-Q_j\right)$

In trying to reduce the number of variables in the problem, we would like a way to simplify the above expression. Specifically, we would like to address the portion $\omega_{j+1}\psi_{j+1} – \omega_j\psi_j – \frac{1}{2}\left(Q_{j+1}-Q_j\right)$

We know that the value of $$\psi$$ evaluated at an interface is some arbitrary constant. Typically, we normalize $$\psi_1(x,\eta_1) = 0$$ and let $$\psi_1(x,\eta_2-d_2) = -M_1$$, which is defined as the mass-flux. There are several options for normalizing the remaining values of $$\psi$$ along each interface. They are

1. Zero at top. We could normalize such that $$\psi_j(x,\eta_j-d_j) = 0$$, and $$\psi_j(x,\eta_{j+1}-d_{j+1}) = -M_j$$. This would reduce the above difference to be $$\omega_{j+1}M_j – \frac{1}{2}\left(Q_{j+1} – Q_j\right)$$.
2. Continuity of the Stream-functions. We could simply normalize such that $$\psi_j(x,\eta_{j+1}-d_{j+1}) = \psi_{j+1}(x,\eta_{j+1}-d_{j+1})$$, and still define $$\psi_j(x,\eta_{j+1}-d_{j+1}) = -M_j$$. This would reduce the above different to be $$(\omega_{j}-\omega_{j+1})M_j – \frac{1}{2}\left(Q_{j+1}-Q_j\right)$$.
3. Canceling the Constants. Perhaps there is a better normalization. If we could say that there was a normalization such that $$\left(\omega_{j+1}\psi_{j+1} – \omega_j\psi_j\right)\vert_{z = \eta_{j+1}} = \frac{1}{2}\left(Q_{j+1}-Q_j\right)$$ then we could dramatically simplify the above calculation.
4. Alternative. Just define $$Q_j$$ with $$Q_1 = 2\omega_1\psi_1^+$$ such that
\begin{align*}
Q_{j+1} &= Q_j + 2\left(\omega_{j+1}\psi_{j+1}^+ – \omega_j\psi_j^-\right)\\
&=Q_{j-1} + 2\left(\omega_{j+1}\psi_{j+1}^- – \omega_{j-1}\psi_{j-1}j^-\right) + 2\left(\omega_{j+1}\psi_{j+1}^+ – \omega_j\psi_j^-\right)\\
&= Q_{j-1} + 2\omega_{j+1}\psi_{j+1}^+ – 2\omega M_j – 2\omega_{j-1}\psi_{j-1}^-\\
&=2\omega_{j+1}\psi_{j+1}^+ – 2\sum_{k = 1}^j\omega_kM_k
\end{align*}
If we proceeded with the third or forth option, we would then find that $$\psi_{z,j+1} = \pm \psi_{z,j}$$ when evaluated at $$z = \eta_{j+1}-d_{j+1}$$. For simplicity, let’s define that $$\psi_{z,j+1} = \sigma_j\psi_{z,j}$$ when evaluated at $$z = \eta_{j+1}-d_{j+1}$$.

### Finding the Mass Flux

The above option relies on determining the mass-flux $$M$$.

## Reducing the Equations

Introduce the variables $$q^+_j = \psi_{z,j}(x,\eta_{j}-d_j)$$ (top of the layer) and $$q^-_j = \psi_{z,j}(x,\eta_{j+1}-d_{j+1})$$ (bottom of the layer), then we have $$q^+_{j+1} + \sigma_j q^-_j$$.

Noting that we can re-write the normal derivatives as

then the integral relationships become

We can then write the equations solely in terms of $$\eta$$ and $$q_j^+$$. Thus,

## Determining the Dispersion Relationship

Now that we have the system of integral equations, we should note some special cases of $$q_j^+$$.

At the top surface: When $$j = 1$$, we know that $$q_1^-$$ can be expressed in terms of $$\eta_1$$ as

At the bottom surface: When $$j = N$$, we know that we can express the term $$q_{N+1}^+$$ as

### Putting it all together

#### Middle layers $$j = 2, .. N-1$$

We can also collapse the above! Note that there is a recursive structure that will allow us to collapse and eliminate the $$q_j^+$$’s.

To determine the value of the mass $$M$$, need to do more work.

## Using the above to find the dispersion relationship

Still a work in progress… At this point, the idea is as follows
formally eliminate the $$q_j^+$$… To be continued.