GenCade:Numerical Stability

The GenCade model utilizes an explicit solution scheme. The main advantages are: easy programming, simple (and sometimes the only possible) expressions of boundary conditions, and shorter computer run-time as compared to an implicit scheme (for a single time increment). A major disadvantage is, however, the stability of the solution. This means that smaller time steps are often needed, and thus a longer simulation time as compared to an implicit scheme. To minimize the computational effort, the longest time step that may be used for a specific calculation must be determined. Under certain idealized conditions, the CERC equation can be reduced to a simpler form to examine the dependence of the solution on the time and space steps. First, rewrite the CERC equation in the form:

$Q=Q_{O}[\alpha _{1}sin2\alpha _{bs}-a_{2}cos(\alpha _{bs}){\frac {\partial H_{b}}{\partial x}}]$

where $Q_{O}=H_{b}^{2}C_{gb}$ (cubic meters/second). A useful approximate stability criterion can be obtained by linearizing the governing equation with respect to y. The linearization is made by assuming small breaking wave and shoreline angles, which leads to:

$sin2\alpha _{bs}\approx 2\alpha _{bs}$ and $\alpha _{bs}=\alpha _{b}-atan{\frac {\partial y}{\partial x}}\approx \alpha _{b}-{\frac {\partial y}{\partial x}}$

where $\alpha _{bs}$ is the angle of breaking waves to the local shoreline orientation, $\alpha _{b}$ is the angle of breaking waves to the x-axis, and ${\frac {\partial y}{\partial x}}$ is the local shoreline orientation. Assuming that ${\frac {\partial q}{\partial x}}$ is zero, the governing equation becomes (Kraus and Harikai 1983):

${\frac {\partial y}{\partial t}}=(\varepsilon _{1}+\varepsilon _{2}){\frac {\partial ^{2}y}{\partial x^{2}}}$

where $\varepsilon _{1}={\frac {2Q_{O}a_{1}}{D_{B}+D_{C}}}$ and $\varepsilon _{2}={\frac {Q_{O}a_{2}sin\alpha _{b}}{D_{B}+D_{C}}}{\frac {\partial H_{b}}{\partial x}}$

In the presence of an external current not generated by breaking waves, the second diffusion coefficient in the above partial differential will change to:

$\varepsilon _{2}={\frac {Q_{O}sin\alpha _{b}}{D_{B}+D_{C}}}(a_{2}{\frac {\partial H_{b}}{\partial x}}-a_{3}{\frac {{\bar {v}}_{t}}{u_{m}}})$

As the above partial differential is a diffusion-type equation, its stability properties are well known. The numerical stability of the calculation scheme is governed by:

$R_{s}={\frac {\Delta t(\varepsilon _{1}+\varepsilon _{2})}{(\Delta x)^{2}}}$

where the quantity $R_{s}$ is known as the Courant number in numerical methods; here it is called the stability parameter. The finite difference form of the above governing partial differential equation shows that $\Delta y\approx {\frac {\Delta t}{(\Delta x)^{2}}}$ which means that if $\Delta x$ is reduced by a factor of two, the time step will need to be reduced by a factor of four to maintain the same stability of the calculation scheme. If an explicit solution scheme is used to solve the diffusion equation, the following condition must be satisfied (Crank 1975):

$R_{s}\leq 0.5$

If the value of $R_{s}$ exceeds 0.5 at any point on the grid, the calculated shoreline will start to become unstable and show unrealistic oscillations in time. The parameters $\varepsilon _{1}$ and $\varepsilon _{2}$ can change substantially alongshore since they depend on the local wave conditions. Assuming that the grid cell spacing is fixed by engineering requirements, a large wave height would necessitate a small value of $\Delta t$. The GenCade model will issue a warning if the stability requirement is violated at any point in the domain. If such a warning is issued, either the time step $\Delta t$ or the spatial resolution ${\frac {1}{\Delta x}}$ (or both) will need to be reduced. Thus, it is necessary to increase $\Delta t$ and/or decrease $\Delta x$.

Useful Links

GenCade References

GenCade Home Page

GenCade Documentation Portal