Numerical scheme
For convenience of numerical realisation, we modify the original problem (1, 3 ,4) and write it in the following form
where the operator At denotes the change between time steps. Equations (7) mean that values in brackets at two time steps are equal to one another. This formulation does not require explicit specification of functions J(a,c) and ff(a,c), which are specified implicitly by initial conditions. We also modify the dynamic free surface condition by adding various supplementary artificial and physical terms on the right-hand side
We use the following set of additional terms
where к is the surface curvature. The first term in (9) introduces damping of displacement of surface particles with the damping coefficient k(a). This term is used for absorbing waves approaching a boundary of a numerical wave tank opposite to a wavemaker and minimising reflections. The second term represents damping of surface curvature and a(a) is the corresponding damping coefficient. As described later in this chapter, it is used to simulate the dissipative effects of wave breaking. The third term represents surface tension. The coefficient 7 is the ratio of the surface tension over density. In calculations presented in this work we use the value of 7 = 0.00073m3/s2 corresponding to fresh water at 20°C. The last term is the prescribed surface pressure gradient, which can be used to create a pneumatic wave generator. The set of Equations (7) with boundary conditions (8, 5, 6) is solved numerically using a finite-difference technique described below.
Since Equations (7) for internal points of a computational domain include only first order spatial derivatives, a compact four-point Keller box scheme (Keller, 1971) can be used for finite-difference approximation of these equations. For our selection of the Lagrangian computational domain the stencil box can be chosen with sides parallel to the axes of the Lagrangian coordinate system, which significantly simplifies the final numerical scheme. The values of the unknown functions x and z on the sides of the stencil box are calculated as averages of the values at adjacent points and then used to approximate the derivatives across the box by first-order differences. The scheme provides the second-order approximation for the box central point. Time derivatives in the second Equation in (7) are approximated by second-order backward differences. It should be noted that the same scheme can be constructed by applying conservation of volume and circulation to elementary rectangular contours with linear approximation of unknown functions on boundaries of elementary volumes. Spatial derivatives in the free-surface boundary condition (8) are approximated by second-order central differences and second-order backward differences are used to approximate time derivatives. As demonstrated below, this leads to a stable numerical scheme with weak dissipation. The overall numerical scheme is of second order accuracy in both time and space.
For finite-difference approximation of a boundary-value problem the number of algebraic mesh equations must be equal to the number of unknown values of functions at grid nodes. Let us consider a finite-difference approximation of the problem (7, 8, 5, 6) on a rectangular N x M mesh, where N and M are the numbers of mesh points in a and c directions. We are required to calculate values of unknown functions x and z at each mesh point. This gives 2 x N x M unknowns which should be determined by solving the same number of finite difference equations. Keller-box approximation of (7) gives two equations for every mesh cell or 2 x (N — 1) x (M — 1) equations. Approximation of boundary conditions (8, 5, 6) at internal points of corresponding boundaries results in 2 x (N — 2) + 2 x (M — 2) equations. Numerical tests demonstrated that both bottom condition (5) and vertical wall conditions (6) should be satisfied at lower corner points to provide their stability (4 equations). The last two equations are provided by vertical boundary conditions (6) at upper corner points. Altogether, this adds up to the required 2 x N x M equations.
A fully-implicit time marching is applied and the Newton-Raphson method is used on each time step to solve the nonlinear algebraic difference equations. It is important to note that the scheme uses only 4 mesh points in the corners of the box for internal points of the fluid domain. Therefore, the resulting Jacobi matrix used by nonlinear Newton iterations has a sparse 4-diagonal structure and can be effectively inverted using algorithms which are faster and less demanding for computational memory than general algorithms of matrix inversion. The current version of the solver uses a standard routine for inversion of general sparse matrices (NAG, 2016). To reduce calculation time, the inversion of a Jacobi matrix is performed at a first step of Newton iterations and if iterations start to diverge. Otherwise, a previously calculated inverse Jacobi matrix is used. Usually only one matrix inversion per time step is required. To start time marching, positions of fluid particles at three initial time steps should be provided, which specifies initial conditions for both particle positions and velocities. An adaptive mesh is used in the horizontal direction with an algorithm based on the shape of the free surface in Lagrangian coordinates z(a, 0, t) to refine the mesh at each time step in regions of high surface gradients and curvatures. Constant mesh refinement near the free surface is used in the vertical direction.
Since the finite-difference approximation presented above uses quadrangular mesh cells, it may be subjected to the so called alternating errors caused by non-physical deformations of cells (Hilt et ah, 1970; Chan, 1975). This deformation preserves the overall volume (area) of a cell, as prescribed by finite difference approximation of the Lagrangian continuity equation. However, volumes of triangles built on cell vertexes can change. The area of one of the triangles increases while the area of the second triangle decreases by the same value. For a cell which is originally rectangular, this corresponds to trapezoidal deformation. This effect finally leads to alternating trapezoidal distortions of neighbouring cells occupying the whole domain, resulting in instability of short-wave disturbances with a wavelength equal to two cell spacing. In earlier models it was usual to implement artificial smoothing to suppress this instability (e.g., Chan, 1975). For the current model, efficient suppression of the instability caused by alternating errors is provided by an adaptive mesh. After several time steps the solution is transferred to a new mesh using quadratic interpolation. The effect of this procedure is similar to regridcling used by Dommermuth and Yue (1987b). It reduces short-wave disturbances and does not produce unwanted damping. Alternating errors still remain the main reason for a breakdown in calculations. However, this happens for a large deformation of the computational domain associated with break in continuity of a physical domain caused, for example, by wave breaking. Otherwise, the numerical scheme proves to be stable with respect to this type of instability.
