Kazuho Fujii


In numerical analysis with the finite difference method, the fact is known that a time step can not be more than some values for calculation stability. For example, The condition of the Courant's number and the diffusion number are famous: respectively $u\Delta t/\Delta < 1$ in the advection equation (CFL condition) and $D\Delta t/\Delta^2 < 0.5$ in the diffusion equation, where $\Delta$, $\Delta t$, $u$ and $D$ are respectively the grid size, the time step, the flow speed and the diffusion coefficient.

However, calculation about free surface can be unstable whenever it satisfies these conditions. It requires yet another condition for stability.

In this report, I think about calculation stability with von Neumann’s analysis method in the case of the small amplitude wave.

Analyzed Equations

I think about Euler’s equation and the equation of continuity: $$\begin{align} \frac{\partial \bm{u}}{\partial t} + (\bm{u} \cdot \nabla)\bm{u} &= - \rho^{-1}\nabla P, \\ \nabla \cdot \bm{u} &= 0, \end{align}$$ where $\bm{u}=(u,v,w), \rho, P$ respectively mean velocity of fluid, density of fluid, pressure including gravity potential.

I can not apply the von Neumann’s method to the raw Euler’s equation, which has a nonlinear term. I assume that surface elevation $\eta$ is enough small and depth of water is enough deep. In the case, I can remove the nonlinear term because velocity fluctuations by the wave are negligible small. That is, $$\begin{align} \frac{\partial \bm{u}}{\partial t} &= - \rho^{-1}\nabla P \label{linear-euler}, \\ \nabla \cdot \bm{u} &= 0. \end{align}$$ Boundary conditions in the assumptions are as follows: $$\begin{align} \bm{u} = 0 &{\quad\mathrm{as}\quad}z \to -\infty, \\ - \frac{P}{\rho} = - g\eta + \gamma \left( \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} \right) \eta &{\quad\mathrm{on}\quad}z=0, \\ \frac{\partial \eta}{\partial t} = w &{\quad\mathrm{on}\quad}z=0, \end{align}$$ where $g$ means gravity acceleration and $\gamma$ means surface tension divided by density $\rho$. The boundary at $z = 0$ means the water surface.

Stability Conditions

In this section, I derive stability conditions for calculating free surface flow. I think about the effects of horizontal and time discretization. For vertically the domain is not discretized here. For spacial difference, I analyzed three grid systems: the regular grid system and the collocate grid system and the staggered grid system; for time integration, I analyzed the Euler’s explicit scheme. I write detail just about the case of using the regular grid system and the Euler’s explicit scheme. Other cases can be analyzed in the same way.

I express intervals between grids for $x$ and $y$ directions respectively with $\Delta x$ and $\Delta y$ respectively, express time step with $\Delta t$, express values of $u, v, w, P, \eta$ at a difference point ($x = l \Delta x, y = m \Delta y, t = n \Delta t$) respectively with $u^n_{l,m}(z), v^n_{l,m}(z), w^n_{l,m}(z), P^n_{l,m}(z), \eta^n_{l,m}$. And I express vertical differential $\partial/\partial z, \partial^2/\partial z^2$ with $()', ()''$.

When governing equations are discretized with the regular grid system and the Euler’s explicit scheme, they are written as follows: $$\begin{align} \frac{u^{n+1}_{l,m} - u^{n}_{l,m}}{\Delta t} &= - \rho^{-1}\frac{- P^{n}_{l-1,m} + P^{n}_{l+1,m}}{2\Delta x}, \\ \frac{v^{n+1}_{l,m} - v^{n}_{l,m}}{\Delta t} &= - \rho^{-1}\frac{- P^{n}_{l,m-1} + P^{n}_{l,m+1}}{2\Delta y}, \\ \frac{w^{n+1}_{l,m}-w^{n}_{l,m}}{\Delta t} &= - \rho^{-1}(P^{n}_{l,m})', \end{align}$$ $$\begin{align} 0 = \frac{-u^{n}_{l-1,m}+u^{n}_{l+1,m}}{2\Delta x} + \frac{-v^{n}_{l,m-1}+v^{n}_{l,m+1}}{2\Delta y} + (w^n_{l,m})'. \end{align}$$ Discretized boundary condition are following: $$\begin{align} u^n_{l,m},\,v^n_{l,m},\,w^n_{l,m} = 0 {\quad\mathrm{as}\quad}z \to -\infty, \end{align}$$ $$\begin{align} - \frac{P^n_{l,m}}{\rho} = - g\eta^n_{l,m} + \gamma \left( \frac{\eta^n_{l-1,m} - 2\eta^n_{l,m} + \eta^n_{l+1,m}}{\Delta x^2} + \frac{\eta^n_{l,m-1} - 2\eta^n_{l,m} + \eta^n_{l,m+1}}{\Delta y^2} \right) \nonumber\\ {\quad\mathrm{on}\quad} z=0, \end{align}$$ $$\begin{align} \frac{\eta^{n+1}_{l,m} - \eta^{n}_{l,m}}{\Delta t} = w^{n+1}_{l,m} {\quad\mathrm{on}\quad}z=0. \label{kbc_discritize} \end{align}$$

I transform $u^n_{l,m}, v^n_{l,m}, w^n_{l,m}, P^n_{l,m}, \eta^n_{l,m}$ to Fourier series as below: $$\begin{align} u^n_{l,m} &= \sum_{k_x}\sum_{k_y} U^n \exp(ik_x l\Delta x + ik_y m\Delta y), \\ v^n_{l,m} &= \sum_{k_x}\sum_{k_y} V^n \exp(ik_x l\Delta x + ik_y m\Delta y), \\ w^n_{l,m} &= \sum_{k_x}\sum_{k_y} W^n \exp(ik_x l\Delta x + ik_y m\Delta y), \\ P^n_{l,m} &= \sum_{k_x}\sum_{k_y} \Phi^n \exp(ik_x l\Delta x + ik_y m\Delta y), \\ \eta^n_{l,m} &= \sum_{k_x}\sum_{k_y} H^n \exp(ik_x l\Delta x + ik_y m\Delta y). \end{align}$$

The equations about Fourier series at a wave number $(k_x, k_y)$ are $$\begin{align} \frac{U^{n+1}-U^{n}}{\Delta t} &= - \rho^{-1}\frac{- \exp(-i k_x \Delta x) + \exp(i k_x \Delta x)}{2\Delta x} \Phi^n \label{euler_fourier1}, \\ \frac{V^{n+1}-V^{n}}{\Delta t} &= -\rho^{-1}\frac{-\exp(-i k_y \Delta y) + \exp(i k_y \Delta y)}{2\Delta y} \Phi^n \label{euler_fourier2}, \\ \frac{W^{n+1}-W^{n}}{\Delta t} &= -\rho^{-1}(\Phi^n)', \label{euler_fourier3} \end{align}$$ $$\begin{align} \frac{- \exp(- i k_x \Delta x) + \exp(i k_x \Delta x)}{2\Delta x}U^n + \frac{-\exp(- i k_y \Delta y) + \exp(i k_y \Delta y)}{2\Delta y}V^n \nonumber\\ +(W^n)' = 0, \label{conti_fourier} \end{align}$$ $$\begin{align} U^n,V^n,W^n = 0 &{\quad\mathrm{as}\quad}z \to -\infty, \label{bc_bottom_fourier}\\ - \rho^{-1}\Phi^n = - (g + K^2 \gamma) H^n &{\quad\mathrm{on}\quad}z=0, \label{dbc_bottom_fourier} \\ \frac{H^{n+1} - H^n}{\Delta t} = W^{n+1} &{\quad\mathrm{on}\quad}z=0, \label{kbc_bottom_fourier} \end{align}$$ where $$\begin{align} K := \sqrt{\frac{2(1-\cos (k_x \Delta x))}{\Delta x^2} + \frac{2(1-\cos (k_y \Delta y))}{\Delta y^2}}. \end{align}$$

With equations $\eqref{euler_fourier1}\eqref{euler_fourier2}\eqref{euler_fourier3}\eqref{conti_fourier}$, I get $$\begin{align} -\kappa^2 \Phi^n + (\Phi^n)'' = 0, \end{align}$$ where $$\begin{align} \kappa := \sqrt{\frac{2(1-\cos (2k_x \Delta x))}{(2\Delta x)^2} + \frac{2(1-\cos (2k_y \Delta y))}{(2\Delta y)^2}}. \end{align}$$

With boundary conditions $\eqref{bc_bottom_fourier}\eqref{dbc_bottom_fourier}$, I get $$\begin{align} \Phi^n = \rho(g+\gamma K^2)H^n \exp (\kappa z). \end{align}$$ With $\eqref{euler_fourier3} \eqref{kbc_bottom_fourier}$, I get $$\begin{align} \frac{1}{\Delta t}\left(\frac{H^{n+2} - H^{n+1}}{\Delta t} - \frac{H^{n+1} - H^{n}}{\Delta t}\right) = -(g+\gamma K^2)\kappa H^{n+1}. \end{align}$$ Organizing this, I get $$\begin{align} H^{n+2} - 2BH^{n+1} + H^{n} = 0 \label{recurrence2}, \end{align}$$ where $$\begin{align} B := 1 - \frac{\Delta t^2}{2}\left(g + \gamma K^2 \right)\kappa. \end{align}$$

Defining $\alpha$ and $\beta$ as $$\begin{align} \alpha := B + \sqrt{B^2-1} \label{solution_alpha},\\ \beta := B - \sqrt{B^2-1} \label{solution_beta}, \end{align}$$ the recursion about $H_n$ $\eqref{recurrence2}$ can be transformed as $$\begin{align} H^{n+2} - \alpha H^{n+1} = \beta (H^{n+1} - \alpha H^{n}), \\ H^{n+2} - \beta H^{n+1} = \alpha (H^{n+1} - \beta H^{n}). \end{align}$$

The condition about calculation stability is $$\begin{align} |\alpha| < 1, |\beta| < 1 \quad\mathrm{for}\ \mathrm{all}\ k_x, k_y. \end{align}$$

In the case of $B^2 - 1 \leq 0$, the calculation is stable because $\alpha$ and $\beta$ are complex numbers and $ |\alpha|^2 = |\beta|^2 = 2B^2 - 1 \leq 1 $ . In the case of $B^2 - 1 > 0$, the calculation is unstable because $B > 1$ if $\alpha > 1$ and $B < -1$ if $\beta < -1$. Then the stability condition is $$\begin{align} B^2 - 1 \leq 0 \quad\mathrm{for}\ \mathrm{all}\ k_x, k_y. \label{stability_cond} \end{align}$$

I define $\Delta$ as $$\begin{align} \frac{1}{\Delta^2} := \frac{1}{\Delta x^2} + \frac{1}{\Delta y^2}, \, \Delta > 0. \end{align}$$ $B$ takes the maximum value $$\begin{align} B = 1 \quad\mathrm{as}\quad K = \kappa = 0. \end{align}$$ It always satisfies $\eqref{stability_cond}$. $B$ takes the minimum value $$\begin{align} B= 1 - \frac{\Delta t^2}{4 \Delta} \left(g + \frac{4}{\Delta^2}\gamma\right) {\quad\mathrm{as}\quad}K = \frac{2}{\Delta}, \kappa = \frac{1}{2\Delta}. \end{align}$$ Therefore, in order to satisfy $\eqref{stability_cond}$ the minimum value of $B$ satisfies it.

In conclusion, the condition of calculation stability is $$\begin{align} \frac{\Delta t}{2\Delta}\sqrt{\frac{\Delta}{2}g + \frac{2}{\Delta}\gamma} < 1. \label{waveCFL} \end{align}$$ In the cases of the collocate grid system, the same result is gotten. In the case of the staggered grid system, I can get the stability condition with the same way: $$\begin{align} \frac{\Delta t}{\Delta}\sqrt{\frac{\Delta}{2}g + \frac{2}{\Delta}\gamma} < 1. \label{waveCFL2} \end{align}$$

The condition is more strict with the staggered grid system because the resolution with the staggered grid system is higher with the regular or collocate system. We take difference between points having $2\Delta x$ or $2\Delta y$ distance with the regular or collocate grid system. However, we take difference between points having $\Delta x$ or $\Delta y$ distance in the staggered grid system.


As $\sqrt{(\Delta /2)g + (2/\Delta)\gamma}$ means the velocity of a wave having $\pi \Delta$ length, the condition $\eqref{waveCFL}$ or $\eqref{waveCFL2}$ can be explained as CFL condition about the $\pi \Delta$ length wave, which is a capillary wave of the grid size. This condition is usually more strict than CFL condition about flow because $\Delta^{−3/2} > \Delta^{−1}$.

Although viscosity and water depth are ignored, they don’t affect the result because they reduce wave speed.

When we don’t resolute the transport of the minimum wave which length is $2 \Delta$, we can conduct a stable calculation. This is because velocity of wave is reduce by discretization as below.

When express $\alpha$ defined at $\eqref{solution_alpha}$ as $$\begin{align} \alpha = |\alpha|\exp{i\theta_\alpha}, \end{align}$$ $\theta_\alpha$ is $$\begin{align} \theta_\alpha = \tan^{-1} \sqrt{B^{-2}-1}. \end{align}$$ The equation $\eqref{waveCFL}$ also means the condition where $\theta_\alpha$ exists for all wave numbers. $\theta_\alpha$ means the phase shift by a time step calculation in a wave which wave number is $k$. Because an angular frequency $\sigma$ is a phase shift per a unit time, $$\begin{align} \sigma &= \frac{\theta_\alpha}{\Delta t} \\ &= \frac{1}{\Delta t}\tan^{-1} \sqrt{B^{-2}-1} \\ &= \frac{1}{\Delta t}\tan^{-1} \sqrt{\left\{1 - \frac{\Delta t^2}{2}\left(g +\gamma K^2\right)\kappa\right\}^{-2}-1}. \label{DescritizeAngularFrequency} \end{align}$$ The case about $\beta$ is the same.

Compared with the angular frequency with the non-discretized case $$\begin{align} \sigma = \sqrt{gk+\gamma k^3}, \end{align}$$ they are matched with each other in low wave numbers. However, the discretized case is later than the non-discretized case in high wave numbers.


In this report I derived another stability condition of a numerical calculation about free surface flow. I note that we must resolute transportation of the grid size capitally wave. This condition is not applied only in calculation of surface wave but also applied in calculation of maltiphase flow with level-set method or volume-of-fluid method or lattice Boltzmann method.

Additionally, I comment that the mean flow is ignored here. In the case where the mean flow or fluctuation of large waves affect, small waves are transported on them. It means that the stability condition is to be more strict.