Numerical dispersion relation and dispersion correction

Special attent ion must be paid to approximation of second time derivatives in the free surface condition (4) since it defines the form of the numerical dispersion relation and is crucial for the overall stability of the scheme. For simplicity let us first consider a case of continuous spatial field in (1, 3, 4) combined with implicit discrete time approximation in (4). Let us approximate second derivatives by 3-point backward differences and expand this approximation to Taylor series with respect to a small time step r. We get

As can be seen, the approximation is of the first order with the leading term of the error proportional to the third derivative of a function, which gives the main contribution to the error of the dispersion relation. Under an assumption of small perturbations of original particle positions we represent unknown functions in the form

and keep only linear terms of expansions with respect to the small displacement amplitude £ —> 0. Introducing a displacement potential ф

we satisfy the vorticity conservation (3) to the first order as £ —¥ 0 and the corresponding approximation of the continuity Equation (1) is the Laplace equation for ф. The dynamic surface condition (4) becomes

where primes denote time derivatives and only the leading term of the approximation error from (10) is taken into account. To derive the numerical dispersion relation we are looking for a solution in the form of a regular wave in deep water:

which satisfies the Laplace equation. The dynamic condition (11) is satisfied when ш and к are related by a dispersion relation. Similar analysis can be performed for higher orders of approximation of the derivatives. Below is the summary of dispersion relations obtained for orders n = 1 to 4:

We use a nondimensional expansion parameter f = fgkr, which is the measure of the problem discretisation representing the ratio of the time step to a typical problem period. As can be seen, the first-order scheme (12a) introduces numerical viscosity proportional to f which leads to fast non-physical decay of perturbations. The higher-order schemes (12c,12d) include terms proportional to —г, leading to growth of perturbations, making the numerical scheme unstable. We therefore use the second-order scheme (12b), which incorporates a numerical error to dispersion at the second order f² and a weak third order (t³) dissipation term.

Let us now include spatial discretisation according to the numerical scheme described above with discretisation steps (8a: Sc). As before, we consider a linearised approximation for a regular travelling wave in deep water. However, differential approximation of the discretised field Equations (1,3) includes higher spatial derivatives and the solution can not be represented in the form of the displacement potential. A suitable form of the solution is

where the constant for the exponential decay of displacement with depth я is not equal to the wavenumber. The expression for я to satisfy the discrete versions of (1) and (3) can be found as an expansion by discretisation parameters. As before, the expansion for ui defines the numerical dispersion relation and is used to satisfy the free-surface condition (4). The corresponding expansions are found to be

and

where the nondimensional discretisation steps f = s/дкт, 8 a = к 8a and 8c = к 8c are used. It is interesting to note that the dispersive error is affected only by the horizontal discretisation step, while the vertical discretisation affects wave kinematics. Therefore, if we are interested in the evolution of the waveform alone, we can use relatively few vertical mesh points.

Validation tests presented later in this section show that the dispersion error is crucial for travelling waves, and the achieved convergence rate for the second-order dispersion approximation is not sufficient. To increase the approximation order for the dispersion relation we introduce dispersion correction terms to the free surface boundary condition. These terms should satisfy the following conditions: (i) to have the order of ()(t²: 8a²): (ii) to be linear; (iii) to not include high derivatives; (iv) to use the same stencil as the original scheme and (v) to reduce the order of the dispersion error. It has been found that the free surface boundary condition (4) with the dispersion correction term satisfying these conditions can be written as follows

where the dispersion correction term is given in parentheses. The term x_aa with the second- order spatial derivative leads to high-wavenumber nonlinear instability for large wave amplitudes. To suppress this instability, we apply 5-point quadratic smoothing to the function x(a) before applying the finite-difference operator. The smoothing is applied to this term in (14) only, and at the future time layer to ensure that the scheme remains fully implicit. It can be shown that the numerical dispersion relation becomes

We now have weak numerical dissipation at 3-rcl order and the dispersion error at 4-th order. A term with weak negative dissipation at 5-th order should also be noted. This term can potentially lead to solution instability for large time steps.