An efficient surface-integral algorithm applied to unsteady gravity waves

A computationally fast method for calculating the unsteady motion of a surface on a two-dimensional fluid is described. Cauchy’s integral theorem is used iteratively to solve Laplace’s equation for successive time derivatives of the surface motion and time-stepping is performed using truncated Taylor series. This allows fairly large time-steps to be made for a given accuracy while the required number of spatial points is minimised by using high order differencing formulae. This reduces the overall number of required calculations. The numerical implementation of the method is found to be accurate and efficient. A fairly thorough examination of this implementation is carried out, revealing that high accuracies are often achievable using surprisingly few numerical surface points. Extensive calculations are also performed using modest computing resources. Some numerical instabilities are identified, although these would not usually be significant in practical calculations. A model analysis reveals that two of these instabilities can be eliminated by using suitable methods of time-stepping. Should the third “steep-wave instability” become significant, it is shown that it can be completely controlled by using high-order smoothing techniques, at little cost to accuracy. Using a routine to ensure asymptotic conservation of energy, this is confirmed by time-stepping a very steep (but stable) wave over thousands of wave-periods. 0 1992 Academx Press. Inc. Since all of the interior properties of a body of fluid undergoing inviscid, incompressible irrotational motion are fully determined by properties at its boundary, it is possible to reduce the calculation of the motion of such a fluid to the evaluation of the motion of its surface alone. In a numerical scheme, the entire motion can thus be modelled using only a point discretisation of the surface. Such an approach was first suggested by Svendsen in 1971 as a means of calculating the motion of two-dimensional gravity waves on water [ 13. Since then the idea has been implemented using a variety of boundary-integral or conformal mapping techniques, by Longuet-Higgins and Cokelet [2], Vinje and Brevig [3], Baker, Meiron, and Orszag [4], Roberts [S], Fornberg [6] and (based on his original formulation) by Svendsen and co-workers [7]. Although this paper presents yet another numerical method, the approach is aimed specifically at producing a computationally fast and efficient numerical scheme in order to be able to study problems of increased complexity. In numerically representing complicated surface motions, ranging from the breaking of waves to the evolution of instabilities in a train of travelling waves, it typically becomes necessary to use a large number of surface points. If running times tend to increase too rapidly with the number of points, then this can destroy any apparent advantage in being able to treat the whole fluid only in terms of its moving surface. The penalty is particularly severe for algorithms that require matrix inversion or factorisation techniques-yielding running times per time step that increase like N3 (the cube of the number of surface points). By using iterative techniques, the penalty can be reduced to running times per time-step that increase in proportion to N2. While being less severe, this still has the effect that increased numbers of surface points become increasingly expensive to calculate. The use of a conformal mapping technique [6] can lead to running times per timestep that vary like N log N. However, this method is severely limited in its ability to resolve steep wave surfaces (as are found in breaking waves), requiring disproportionately many surface points to achieve a given accuracy. It is therefore vital to ensure computational elliciency if it is to be possible accurately to describe wave processes for which large numbers of surface points may be required. The most direct way of achieving this is to employ the lowest number of points that can adequately describe a fluid surface. In this paper, a further improvement is introduced by showing that higher time-derivatives of the motion can be calculated directly using exactly the same technique that is used to determine the velocity of the surface at each timestep. This yields significantly more information about the motion for relatively little extra computation. In turn, this allows larger time-steps to be performed for a given accuracy so that the overall number of calculations can be reduced. The method presented here combines these two ideas into producing a practical numerical algorithm, the properties of which are examined in the following pages. In Appendix A, the formulae are derived that make it