On reference solutions for ship waves

This paper presents asymptotic solutions to the classic problem of stationary deep water ship waves first formulated by Lord Kelvin.¹ The Kelvin analysis of wave kinematics predicts universal wake patterns in the form of arms of chevron accompanied by a system of transverse arcs inside a gusset of total opening at approximately 39° (twice the Kelvin angle, θ_K ≈ 19.47°). This analysis was later extended^2–5 by studies of amplitude distributions based on asymptotical approaches of the stationary phase⁶ (see also the work of Thomson⁷) and the steepest descent.^8,9 It should be mentioned that these approaches showed their efficiency in many wave problems, mostly one-dimensional. An extension of these methods to two and more dimensions is painstaking or even impossible. Quite often, the solution of the multi-dimensional problem is replaced by a unidirectional one to avoid serious mathematical difficulties.

Wave caustics are the main difficulty of the stationary phase method (as well as of its counterpart, the steepest descent approach). The difficulty lies in an essentially point-wise nature of the method when the evolution of a spatial wave pattern of finite volume is replaced by one of the characteristic points that provides phase stationarity. In this way, the concept of wave group velocity appears^6,10 when the evolution can be adequately described as a propagation of the characteristic point of the wavetrain at the group velocity. The evolution of wave amplitude is accounted for by the dispersion of group velocity (second derivatives of the dispersion relationship). The spatial distribution of the wavetrain remains beyond this method. When wave dispersion vanishes at a caustic (or at its projection when reduced to an omnidirectional problem), this approach faces a mathematical problem: there appears singularity of wave amplitude which can be removed by accounting for higher-order effects of wave dispersion. A physical way of avoiding the “mathematical infinity” “to blur or to smooth it down” was realized by Lord Kelvin¹ (p. 425) quite well by using explicit wavetrain distributions² [Eq. (95)]. This physical idea found its mathematical treatment in many works on asymptotic methods and perturbation theory (Ref. 11 and references therein).

The novel approach recently sketched by the authors,¹² the so-called reference solution method, follows the idea of Lord Kelvin. According to this approach, the reference wavetrains are Gaussian, which allow one to develop an analytic approach to many problems of wave propagation in multidimensional dispersive media. The resulting asymptotic solutions are shown to be finite and continuous in the whole space.

In Sec. II, we outline features of the method of reference solutions for dispersive waves of infinitely small amplitude. The general formulas will be used in Sec. III for the classic problem of stationary ship waves. Solutions will be presented in an explicit form where particular terms have transparent physical meaning. In this way, the effects of amplitude and phase transformation can be analyzed separately. Graphical illustrations unveil the general effects of ship wake formation. In conclusion, we discuss the prospects of the proposed approach.

The new approach to the analysis of linear wave packets in the case of two or more spatial dimensions has been proposed recently by the authors.¹² In this section, we present the basics of the method and general formulas for wavetrains in dispersive multidimensional media.

Following the study of Gnevyshev and Badulin,¹² consider solutions for waves of infinitely small amplitude in the form of linear superposition of harmonics in two-dimensional media (e.g., water waves),

Here, η(x, t) is a wave variable (e.g., surface elevation of water waves), Re is the real part of the corresponding expression, and F(k) is the Fourier image of the function η(x, t) at initial time t = 0. The wavevector k = (k, l) and frequency Ω are assumed to obey the dispersion relation Ω(k).

Consider the initial distribution f(x) = η(x, t = 0) in (1) in the form of a Gaussian wavetrain,

Note that at the limit Δx → 0, Δy → 0, (2) is reduced to a δ-like function,

This property will be used essentially to construct solutions for the ship wake of a finite bandwidth. In the k-space distribution, (2) is a normalized Gaussian pulse,

where the wavevector κ = (κ_x, κ_y) denotes the deviation from the carrier wavevector k₀ = (k₀, l₀), i.e., κ_x = (k − k₀), κ_y = (l − l₀), and the spectral bandwidth is related to its dimension in a coordinate space as follows:

Assuming the pulse in (4) is narrow-banded, i.e., Δk/k₀ ≪ 1, Δl/l₀ ≪ 1, one can approximate the dispersion relation Ω(k) in (1) by its Taylor series in the vicinity of the wave packet carrier harmonic k₀ (⁠$Ω0=Ω∣k=k0$⁠),

Herein, we use the following notations for the second order differentials associated with the bandwidth of the reference function (4):

The resulting approximation for (1) takes a form of a modulated wavetrain,

with the phase of the carrier wave being

and the dimensionless wavevector λ and the coordinate ξ that characterizes formally small deviations from the carrier harmonic k₀ and its trajectory in a coordinate space being

and

The modulation function given by the Fourier integral in (8) is determined by the quadratic form with a dimensionless matrix,

The quadratic form exponent in (8) can be integrated explicitly⁸ as follows:

The real part of the matrix B is strictly positively determined that ensures convergence of the integral (8) and existence of the inverse matrix B⁻¹. The determinant of B,

does not vanish for any dispersion relation Ω(k) and, thus, the problem of wave caustics never appears within this approach as it generally occurs in the lowest orders of the widely used stationary phase and steepest descent methods.

Finally, the reference solution can be written as follows:

We used “mixed” notations introducing “fast” x and “slow” ξ(11) variables in (15) in order to discriminate dependence on the wavetrain carrier parameters (“fast” x) and transformation of the wave packet shape associated with small deviations from the harmonic k₀ (“slow” ξ).

In addition to the carrier wave phase ϑ₀ (9), two phase shifts ϑ₁ and ϑ₂ appear in (15). The first one ϑ₁ shifts phases of all the wavetrain components by the same value,

Generally, ϑ₁(k₀, t) is decaying as t⁻¹ is dependent on signatures of the second differentials μ_x, μ_y, and μ_xy. For omnidirectional propagation and in the special case of quasi-dispersion (⁠$μxy2−μxμy=0$⁠, see Ref. 12), ϑ₁ tends toward the well-known result of the stationary phase method ϑ₁ = ±π/4 (e.g., Ref. 9), thus generalizing this result.

The phase shift ϑ₂(ξ, k₀, t) depends on both the wavetrain dispersion parameters and “slow” coordinate ξ,

It vanishes for the wavetrain carrier harmonic (ξ = 0) and changes peripheral harmonics when ξ ≠ 0. For deep water waves at t → +∞, the lines of constant phases of (17) are hyperbolas, and the phase ϑ₂(ξ, k, t) is decaying as t⁻¹.

The argument of the exponential function in (15) preserves the Gaussian shape of the wavetrain modulation. The isolines of Λ(k₀, ξ, t),

are ellipses whose parameters (orientation and eccentricity) vary with time. Equation (18) can be written in the following form:

that clearly illustrates effects of the wave packet rotation and stretching from the initial state $ξ12+ξ22=const1$ to the asymptotic one (t → +∞) when $ξ1μxy−ξ2μx2+ξ2μxy−ξ1μy2=const2$⁠. It should be stressed that the isotropy of initial distribution and/or of the dispersion relation does not cancel these effects.¹² 

The factor D affects the wavetrain modulation in two ways in (15): the magnitude of perturbation decays, but its width in a coordinate space grows to preserve the total energy. In the general case, when the dispersion factor does not vanish, i.e., $Ωkl″2−Ωkk″Ωll″≠0$ (⁠$μxy2−μxμy≠0$⁠), the wavetrain in two-dimensional media decays as t⁻¹. In the special case of quasi-dispersion, $Ωkl″2−Ωkk″Ωll″=0$⁠, the decay becomes slower and follows the one-dimensional asymptotics t^−1/2.

Thus, the method of this paper reproduces the well-known results of the stationary phase approach. Moreover, it allows for the realization of Lord Kelvin’s idea “to smooth, to smear”¹³ formally infinite growth at wave caustics in a consistent way. At the same time, it is not a general method of solution for arbitrary initial (boundary) conditions. We leave the issues of completeness and uniqueness of the solutions beyond our analysis. An important point to emphasize is that the method is multi-dimensional. Thus, it essentially extends similar results for one-dimensional waves¹⁴ and captures new physical effects like the wave packet rotation.

The general approach presented above is successfully applied here to the problem of ship waves. Formally, it is valid for narrow-banded wavetrains (6), while the ship waves are not. In this regard, additional consideration of wave kinematics is required to substantiate the method for the specific problem. We start this section with an analysis of ship wave kinematics first developed by Lord Kelvin.¹ We reproduce and precise some well-known results on phase distributions treated in previous studies.¹⁵ Our analysis within the approach of the reference functions provides a simple and effective tool for quantifying the amplitude distributions in the ship wake, the problem that has many important implications including the case of shallow and intermediate water depth.¹⁶ In particular, new analytic results of the section explain recent findings for ship wakes at high speeds.^17,18

In the reference frame of the ship moving along the x-axis in a negative direction with a constant speed U, the dispersion relation for deep water waves takes the form of the Doppler-shifted frequency,

Assuming ship waves to be stationary with reference to the ship (Ω = 0), one immediately has the conditions of synchronism that do not depend explicitly on the dispersion relation (type of waves),

where c is the phase speed of deep water waves, c = g/ω with intrinsic wave frequency obeying the dispersion relation $ω=g|k|1/2$⁠, and

where θ is the angle between the opposite direction of axis x and the wavevector direction. Formulas (20)–(22) assume real-valued solutions for the stationarity condition Ω = 0, thus implying the far field approximation of the problem under study.

Our approach defines the carrier wave trajectory as ξ = 0 [(6) and (13)], i.e.,

It is equivalent to the geometric optics approach or the stationary phase method which treats the trajectory as one of the quasi-particles traveling at group speed. For the deep water dispersion relation (20) with conditions of stationarity [(21) and (22)], it gives a parametric form of the trajectory,

Within the problem of stationary ship wake, the time t can be treated as a formal parameter of the stationary phase approach. Indeed, looking for wave phase extrema in (1) with the constraint of wave system stationarity Ω(k) = 0 (20), we obtain the time t as the Lagrangian multiplier and Eq. (21) for the conditional phase extrema.¹⁵ Alternatively, within our approach, the carrier harmonic with k and x is fixed by the reference solution shape and does not imply solution of the extremum problem.

The trajectories of (23) are straight lines, and the phase of the carrier harmonic ϑ₀ of (9) along these trajectories can be written in a trivial way,

It leads to the well-known parametric solutions of Kelvin–Lighthill^1,5 in terms of phase ϑ₀ and angle θ,

Accordingly, lines of constant phases ϑ₁ and ϑ₂ [(16) and (17)] can be parameterized by excluding dependence on parameter t.

Turning points of the curves (26) are determined by stationary points, where dx/dθ = 0, dy/dθ = 0,

i.e., for angle θ₁ ≈ 35.26°. These turning points form a cone with a spread angle

which is the well-known angle of Kelvin.

One can easily show that the stationary point in coordinates (x, y) at the Kelvin angle θ_K (28) corresponds to the inflection point of isofrequency (21) at angle θ = θ₁ (27) in a wavenumber space (k, l), i.e., $∂2⁡k/∂l2|θ=θ1=0$⁠. The latter condition means merging of two wave systems: diverging and transverse at their demarcation line. The inflection point $∂2k/∂l2|θ=θ1=0$ corresponds to maximal group velocity when stationary waves can exist and to the maximal opening angle, the Kelvin angle θ_K, in the coordinate plane (Fig. 2).

Ship wave kinematics is illustrated by two figures. Figure 1 shows wave parameters in a wavevector space, while Fig. 2 describes the phenomenon in the coordinate plane. In Fig. 1(a), two isofrequency curves Ω(k) = 0 [see Eqs. (20) and (21)] are shown for different dimensionless ship speeds and Froude numbers Fr = 0.25; 1.0. The latter can be naturally defined in terms of the bandwidth Δk or perturbation scale Δx,

It should be emphasized that wave kinematics of the carrier wave determined by trajectories (23), i.e., by formally point-wise initial conditions, does not depend on the Froude number (29). The corresponding dependencies appear in dispersive corrections of the wave phases ϑ₁ andϑ₂ [(16) and (17)] and, what is more important, in expressions for wave amplitudes (15) that essentially account for the finite size of an initial perturbation.

The inflection point of the dispersion relation (20) [black dotted lines in Fig. 1(a)] demarcates two systems of ship waves: transverse (shown in blue) and diverging (red) ones. At any coordinate point inside the Kelvin cone [$arctany/x<θK$], two solutions of the condition of synchronism (21) exist that correspond to one transverse (θ < θ₁, blue) and one diverging wave (θ > θ₁, red) with corresponding wavevectors (shown in the same color) and different phase shifts ϑ₀, ϑ₁, andϑ₂ [see general expressions (16) and (17)]. Evidently, wavevectors of each family are the same at any point of a wave trajectory (24), shown by a green line in Fig. 2. These wavevectors (shown by colored lines in Fig. 2) are perpendicular to isophase curves ϑ₀ = const of transverse and diverging waves. The corresponding identical direction of group velocities of both waves (α_T = α_D) for the particular trajectory (green line) is shown by a black arrow in Fig. 2. Note that absolute values of group velocities are different and, hence, phases ϑ₀, ϑ₁, and ϑ₂ of diverging and transverse waves in a point are, generally, different.

Two points of the analysis of Sec. III A should be emphasized. First, all the wave trajectories x(t), y(t) as functions of t are straight lines, and there are no intersections between them (caustics). Second, at every point of these trajectories, two harmonics (transverse and diverging) obeying the synchronism condition (21) coexist. Thus, wave interference at any point of the stationary ship wake is represented by two harmonics only. The stationarity condition (21) wipes out all the satellites of the carrier wave making the wave packet narrow. In this way, the inherently wide-banded problem of ship waves is reduced to a problem of interference of just two narrow-banded wavetrains. Thus, the formulas for the evolution of narrow-banded perturbation can be modified by formal passage to δ-pulse reference wavetrains (3).

In fact, key expressions for wavetrain evolution (15)–(19) contain the bandwidth parameters Δk, Δl explicitly [(6) and (8)],

These small but still finite quantities (Δk and Δl) should be taken into account as they essentially determine wave amplitudes via the determinant D(14).

The bandwidth parameter Δk (Δx) can be associated with the dimension of the traveling ship and should satisfy, formally, the condition of the narrow-banded perturbation that matches the stationarity condition (21). The latter can be re-written for the wavevector modulo æ $=(k2+l2)1/2$ and its maximal perturbation Δæ,

Simple algebra for determinant D in (14) gives

where the parameter

is a length–beam ratio of the ship.

Bandwidth narrowness is provided by Gaussian distribution (4) that can be re-written taking into account the stationarity condition [(21) and (22)] as follows:

The representation (35) illustrates an effect of amplitude (energy) partitioning between different wave systems that depends on the angle parameter θ and the wavevector direction (22). Figure 1(b) shows the correspondence between the ranges of wavenumbers (frequencies) and the wave amplitude content of different wave systems. Waves with relatively low transversal wavenumbers (transversal waves) gain energy from the core of distribution F(k, l) (4), while diverging waves with high transversal components of wavevectors are controlled by low energy tails of F(k, l). Thus, at low Froude numbers Fr, transverse waves (relatively small θ) predominate, while diverging waves with large angles θ (1/cos θ ≫ 1) are suppressed.

At the Kelvin cone cos²θ = 2/3 that demarcates diverging and transverse waves, the “weight” of the corresponding harmonics is written as follows (35):

which means an effective suppression of diverging waves at small Froude numbers. For Fr = 0.3 and β = 5 (length–beam ratio), wave amplitude at the edge of the Kelvin cone appears to be less than 1% of the maximal amplitude of the transverse waves at the wake axis. For Fr = 0.5, this fraction becomes about 53% and decays slower inside the Kelvin cone. This can explain the observed growth of the diverging fraction of the ship wake at high Froude numbers and the visible tendency of the ship wake to shrink at high Froude numbers Fr ≫ 1. The latter effect is interpreted as shrinking of the Kelvin cone,¹⁹ while in terms of our approach is naturally treated as energy redistribution between transverse and diverging wave fractions inside the permanent Kelvin cone.

The limit of small Froude numbers (32) does not simplify analysis. Moreover, the Froude number Fr enters D in (33) in combination with parameters β and ϑ₀ that could be and are indeed large values. Thus, asymptotics

make sense and can be taken into consideration as well.

An additional phase shift ϑ₁ associated with the wave dispersion (16) depends on the bandwidth Δk (Δx or Fr),

At small Froude numbers (32) and finite β, ϑ₀, it gives the result by Lighthill,⁵ 

For large values of the length–beam ratio β and ϑ₀, i.e., at the limit (37), a new characteristic value of angle cos²θ₂ = 1/3 appears in addition to cos²θ₁ = 2/3. This new angle is implied by symmetry of elliptical initial perturbation in our problem [cf. (28)],

For the length–beam ratio β² ≫ 1, the second term in (33) has a minimum near θ_N ≈ 15.8°, i.e., in the Kelvin cone interior. It can be observed in the ship wave patterns irrespective of the value of Froude number.^17,18

At large phases ϑ₀, the ship wake amplitude in (15) tends to simple power-like dependence,

which is also consistent with results of Rabaud and Moisy.^17,18

With explicit formulas at hand, one can easily present the role of different mechanisms in stationary ship wake formation. Formally, the approach implies an essential restriction: the bandwidth of an initial perturbation (boundary conditions in the case of stationary ship waves) should be narrow which leads to a formal requirement Fr cos θ ≪ 1 (32). In fact, this limitation is not essential as long as the stationarity condition (21) guarantees the bandwidth narrowness by wiping out satellite harmonics, thus making the wave packet spectrum narrow. In contrast to the general reference solution (15), which depends on both time t and coordinate x, the stationarity condition (21) reduces the number of arguments of the particular solution (33) and, additionally, prunes wave interference in a wavevector space to a superposition of just two harmonics of diverging and transverse waves.

The resulting modulation of the ship wake spectrum by the Fourier image F(k, l) (4) can be formally used in formulas of Sec. II. Figure 3 presents patterns of ship wakes at different Froude numbers Fr = 0.5, 1.0, and 2.0 (from the top to the bottom). The length–beam ratio β = 10 is fixed. The scaling of Rabaud and Moisy¹⁸ is used to relate the results of different approaches. The amplitudes are normalized by maximal values for x/L_x > 2. The “near-field” range x/L_x < 2 is left blank.

The left column of Fig. 3 shows diverging waves. The central column shows the transverse wave system, and the right one shows the resulting superposition of transverse and diverging waves. The figure illustrates the mentioned amplitude effects. The diverging waves are quite weak at Fr = 0.5, as expected from formulas (35) and (36). At high Froude numbers (middle and bottom rows in Fig. 3), they become dominating as compared to the transverse wave system. A visible collapse of the ship wake cone at high Froude numbers is associated with the mechanism of energy partitioning between different wave systems, as described and illustrated in Fig. 1(b). This transparent physical treatment of the effect of transformation of the ship wake at high Froude numbers does not require an introduction of long scale cutoff of initial perturbation spectra as proposed by Rabaud and Moisy¹⁸ and Darmon et al.¹⁹ 

In this paper, the reference solution approach for linear water waves is presented. It allows for transparent physical analysis of the classic problem of stationary ship waves. Explicit formulas are obtained for wave variables (free surface elevation in our case) without explicit reference to dynamical equations. All the required information is contained in the fact of wave linearity, i.e., the possibility of superposition of partial solutions and in the specific dispersion relation Ω(k).

The approach allows for separate analysis of amplitude and dispersion effects for independent systems of diverging and transverse waves. A transparent physical treatment of the effect of amplitude partitioning (Fig. 1) is presented and illustrated by wave patterns at different Froude numbers. Snapshots of ship wakes in Fig. 3 obtained with explicit analytical formulas show consistency with the results of previous studies based on conventional asymptotic approaches.^17,18 Within the novel analytical approach, the important effect of visible ship wake shrinking at high Froude numbers found its natural physical explanation as a phenomenon of energy partitioning between different families of stationary ship waves (diverging and transverse).

In contrast to the conventional stationary phase and the fastest descent approaches, the two-dimensionality of the problem is not an essential difficulty. Moreover, the resulting formulas allow for easy discrimination between different physical effects associated with wave dispersion. The effect of phase ϑ₁ (38), which is a straightforward counterpart of the phase shift ±π/4 of the method of the stationary phase and the dependencies of wave parameters on the Froude number and length–beam ratio β have not been detailed here.

The novel approach of reference solutions has good perspectives. In particular, the problem of a ship wake in the presence of shear current²⁰ can be analyzed easily within this approach.

The presented method of reference solutions by Gnevyshev and Badulin¹² is associated with a specific class of initial (boundary) conditions and cannot be used in general cases. At the same time, the linearity of the problem essentially extends the approach applicability to qualitative and quantitative analysis as a testbed for numerical simulations.

This research was supported by the Russian Science Foundation (Grant No. 19-72-30028) with the contribution of MIGO GROUP (http://migogroup.ru).