Stability of hydroelastic waves in deep water

Two-dimensional periodic travelling hydroelastic waves on water of infinite depth are investigated. A bifurcation branch is tracked that delineates a family of such solutions connecting small amplitude periodic waves to the large amplitude static state for which the wave is at rest and there is no fluid motion. The stability of these periodic waves is then examined using a surface-variable formulation in which a linearised eigenproblem is stated on the basis of Floquet theory and solved numerically. The eigenspectrum is discussed encompassing both superharmonic and subharmonic perturbations. In the former case the onset of instability via a Tanaka-type collision of eigenvalues at zero is identified. The structure of the eigenvalue spectrum is elucidated as the travelling-wave branch is followed revealing a highly intricate structure. This is dedicated to Professor Jean-Marc Vanden-Broeck on the occasion of his 70th birthday


Introduction
The study of two-dimensional hydroelastic waves has received an increasing amount of attention in recent years.Aside from pure scientific interest, the motivation comes in part from various applications in engineering, such as off-shore floating structures, as well as from the desire to understand certain natural phenomena, for example the dynamics of ice covers on the ocean.For a review of hydroelastic wave dynamics, see [1].Steadily propagating periodic waves and solitary waves have been found, for example by [2], [3], [4], [5] and [6].Questions of existence and well-posedness have also been addressed, for example by [7] and [8].The aforementioned studies included the effect of gravity.Studies that focus purely on hydroelastic effects include [9] and [10].
The stability of periodic travelling waves has also been examined building on the voluminous literature dealing with the stability of pure gravity waves initiated by Benjamin & Feir [11].Relevant studies include those by [12] and by [13], who analysed the stability of two-dimensional flexural-gravity waves using the AFM method.
In the present article we consider the stability of hydroelastic waves on fluid of infinite depth in the absence of gravity.We compute two-dimensional periodic waves which represent the infinite-depth counterparts of the finite-depth waves computed by [9].The linear stability of these waves is examined by linearising about the steady states and using the Floquet-Fourier-Hill method [14] to formulate an eigenvalue problem which is solved numerically.We identify a single branch of periodic travelling-wave solutions that bifurcates from the state in which the free surface is flat.The wave speed decreases as this branch is followed to large amplitude waves, and ultimately the branch terminates when the wave speed reaches zero and the wave profile adopts that for a static elastic sheet.
Notably the wave energy along the bifurcation branch exhibits a local maximum before the static state is reached.At this point we observe the onset of a Tanaka [9] type superharmonic instability via a collision of eigenvalues at zero.We also investigate subharmonic perturbations and demonstrate a highly intricate structure for the eigenspectrum.This includes figure-ofeight type patterns such as those seen in traditional gravity wave stability computations [15,16,17,18,19].
The outline of the paper is as follows.In section 2 we present the governing equations which are formulated following the approach of Dyachenko [20] to express the problem in terms of surface variables only.In section 3 we discuss periodic travelling-wave solutions.In section 4 we formulate the linear stability problem, which is then solved numerically in section 5. Finally we summarise our findings in section 6.

Problem formulation
We consider the flow of water of infinite depth and of density ρ beneath an elastic sheet with the particular goal of examining periodic waves on the surface.The elastic sheet has bending modulus EB which, according to thin-shell theory, is related to the Young's modulus where d is the thickness of the sheet and ν is the Poisson ratio.We use a Cartesian coordinate system (x, y) such that in the simplest state when the surface is flat it is described by y = 0.More generally the surface is deformed and our interest lies in calculating steadily propagating waves and determining their stability.In this case we analyse the dynamics in a frame of reference that is moving at the wave speed c > 0.
Assuming inviscid, irrotational and incompressible flow, the governing equation for the velocity potential in the fluid is Laplace's equation, namely The velocity potential is subject to the far-field condition The kinematic condition at the surface requires that where n is the unit normal to the surface and x describes the location of a material point in the surface at time t.Neglecting gravity, the dynamic condition at the surface requires where s denotes arclength along the surface.We note that the arbitrary function B(t) can be subsumed into ϕ.Condition (2.4) includes the jump in pressure across the surface that is balanced by the elastic force prevailing in the flexible sheet.The latter is represented by the final term on the left hand side and has been derived by considering an equilibrium balance of forces in the flexible sheet (see, for example, [9]).The same form is formally derived by [21].
In the sequel we will need to compute the kinetic energy in the stationary lab frame.This can be done using the formula [22] where E k is the kinetic energy per unit length in the transverse direction, and λ is the wave length.Here y = Y corresponds to the free surface, and φ = ϕ + cx is the velocity potential in the stationary lab frame.We shall also want to compute the elastic energy per unit length in the transverse direction, defined to be (see, for example, [23,5]) where L is the total arclength of the deformed plate in one wave period, and s is arclength along the plate.
In the limit of small amplitude travelling waves, the dispersion relation relating the wave speed c to the wave length λ is (see [24]) We make all quantities dimensionless using λ/2π as the length scale and (ρ/EB) 1/2 (λ/2π) 5/2 as the time scale.In dimensionless units the wavelength is 2π and the dispersion relation (2.7) reduces to c = 1.Henceforth all quantities will be assumed to have been non-dimensionalised in this way.We reformulate the governing equations in terms of surface variables using the conformal mapping approach developed by [25], [20] and [26] for water waves.Briefly the potential flow problem in the physical plane is mapped to the lower half image plane and the governing equations for the flow are manipulated to yield a coupled set of partial differential equations in terms of surface variables only.If the location of the surface is described parametrically by (x(ξ, t), y(ξ, t)), where t is time and −∞ < ξ < ∞, and if ϕ(ξ, t) and ψ(ξ, t) are, respectively, the velocity potential and stream function on the surface, then the non-local system of governing equations is and where B(t) is the Bernoulli constant, and Here J = x 2 ξ + y 2 ξ and the curvature (2.12) The Hilbert transform is defined to be (2.13)

Travelling waves
In this section we construct steadily propagating periodic travelling-wave solutions.To this end we remove the time derivatives from (2.8) and (2.9) and write Representing the surface profile as x = X(ξ) and y = Y (ξ), where X(ξ + 2π) = 2π + X(ξ) and Y (ξ + 2π) = Y (ξ), the dynamic boundary condition (2.9) requires that where where a prime indicates a derivative with respect to ξ.
We can compute a solution numerically using a suitably truncated form of the Fourier representation for the surface profile, taking α0 = 0 without loss of generality, that is Then X(ξ) is obtained by integrating (2.11) and fixing the integration constant so that X = ξ − H (Y ).Note that We use a collocation point method to compute the Fourier coefficients in (3.4).Fixing the collocation points at for j = 1, . . ., 2N + 1 for some integer N , we truncate the Fourier series at |n| = N and evaluate (3.2) at each collocation point to yield the required number of nonlinear algebraic equations for the unknowns comprising the 2N Fourier coefficients αn (n ̸ = 0) and c.The nonlinear set of equations is solved numerically using Newton's method with finite differences employed to approximate the necessary derivatives.The branch of nonlinear travelling wave solutions bifurcates from the flat free surface solution at the wave speed corresponding to the linear dispersion relation noted above, namely c = 1.This is evident from the bifurcation curve shown in the left upper panel of Figure 1.This plot shows how the wave height H ≡ max(Y ) − min(Y ) varies with the wave speed.As can be seen the bifurcation curve approaches c = 0 in the limit.Wave profiles at various points along the bifurcation curve are shown in the left lower panel of Figure 1, together with the static profile.The latter was computed using the method described by [9].The total dimensionless wave energy E varies with the wave speed.Here the total dimensionless wave energy is defined to be where is the dimensionless elastic energy.The kinetic energy defined in (2.5) is most conveniently computed using the equivalent (and dimensionless) formula The computed energy will be discussed later when we come to the linear stability results.
The numerical results can be corroborated by performing a classical Stokes expansion near to the point of bifurcation.The details are given in Appendix C, where the first approximation to the elastic energy is found to be (see (C.38)) This formula is compared with the numerical result in the rightmost panel of Figure 1 demonstrating excellent agreement between the two.

Linear stability
To examine the linear stability of the travelling wave solutions presented in section 3, we perturb about the basic state by writing where the tilde decorated terms are assumed to be small.Substituting into (2.8)-(2.11)and neglecting products of small terms, we obtain Here where Following Tiron & Choi [27] we express the perturbations in the form of Fourier expansions, writing where p ∈ [0, 1) and 'c.c.' denotes the complex conjugate.Using (4.4) we find if (j + p) ̸ = 0 that aj = −i sgn(j + p)bj, dj = i sgn(j + p)cj.
Note that although these relations do not specify a0 and d0 when j = 0 and p = 0, neither of these constants will appear in the perturbation equations (4.2)-(4.4)once we have made the substitutions (4.6), due to the differentiations involved.Equations (4.2) and (4.3) become, respectively, where ) ) and where Pe(ξ) is given in Appendix A. We note that Pe(ξ) = 0 if j + p = 0.

Numerical method for arbitrary wave amplitude
For arbitrary amplitude travelling waves, the linear stability problem must be solved numerically.We achieve this by truncating the infinite sums at N terms, for a suitable choice of integer N , and we use the 2N + 1 equally-spaced collocation points defined in (3.7).Applying (4.9) and (4.10) at the grid points we obtain a set of algebraic equations which we form into the matrix system where x = (b−N , . . ., bN , c−N , . . ., cN ) T , and A and B are (4N + 2) × (4N + 2) coefficient matrices.We solve (4.24) numerically to calculate the eigenvalues σ using the in-built Matlab routine eigs.As has been noted by previous authors, in the superharmonic case of p = 0 the system (4.24) has a zero eigenvalue with algebraic multiplicity four.Two of these stem from the fact that, as was noted by [27], A has its (N + 1)th and (3N + 2)th columns both zero.Thus we may constructed two eigenvectors corresponding to the eigenvalue σ = 0 which are filled with zeros except in the (N + 1)th and (3N + 2)th entries, respectively.These correspond to the eigenvectors (b0, c0) = (1, 0) and (0, 1) discussed in section 4.1.Two further null eigenvalues appear due to the freedom to introduce a pure phase shift in the horizontal direction [28].

Numerical results for linear stability
In this section we present numerical results using the method described in section 4.2.In all of the results shown the truncation level N was increased until converged results were obtained.For the most challenging computations we took N = 2048 but typically fewer modes than this are necessary.We discuss first results for superharmonic perturbations with p = 0.As was indicated in section 4.1 the growth rates are all purely imaginary in the limit H → 0. In Figure 2 given by (4.22).The eigenvalues on the branches labelled σ (0,1) − and σ (0,−1) + collide at zero when c = 0.415.This gives rise to what [29] have referred to as a Tanaka-type instability.This is characterised by the vanishing of an eigenvalue at the point where the energy E attains a maximum as a function of c; and, critically, by the fact that the eigenvectors corresponding to the colliding eigenvalues become linearly dependent at the point of collision.To demonstrate such an occurrence here, we must discuss carefully all of the possible zero eigenvalues and how their corresponding eigenvectors behave as the point of collision is approached.
We denote by v1 and v2 the two linearly independent eigenvectors corresponding to the branches labelled σ (0,1) − and σ (0,−1) + in Figure 2: these two eigenvectors are in general linearly independent.Also of interest are the two phase shift modes, whose eigenvectors we denote by v3 and v ′ 3 , which connect to the two flat state modes (4.23) at zero amplitude (however, v3 and v ′ 3 are linearly dependent, and hence we ignore v ′ 3 in the following discussion).We must also consider the two eigenvectors which are present due to the zero columns of the matrix A; we call these v4 and v5, and we note that they are obviously mutually orthogonal.Also v4 and v5 are both orthogonal to v3 since a pure phase shift mode has zero ξ-average.
Consider the set S = {v1, v2, v3} of generally mutually linearly independent vectors.We find that this set collapses to a single eigenvector at the point c = 0.415.To see this we follow [28] and construct the Gram determinant where (x1, x2) = x T 1 x * 2 with the asterisk denoting the complex conjugate.Figure 3(b) shows that G(v1, v3) is non-zero except at c = 0.415 where it vanishes, meaning that v1 and v3 become linearly dependent.Similar graphs can be produced for G(v1, v2) and G(v2, v3).Therefore at c = 0.415 we have a zero eigenvalue of algebraic multiplicity six but geometric multiplicity three, comprising the any one member of S (the other two being linearly dependent) and v4, v5.We have also investigated higher frequency superharmonic modes, and have traced the branches in a graph similar to that shown in Figure 2 for Im(σ) ∈ [−150, 150].We do not find any further collisions of eigenvalues, either at zero or away from zero.This suggests that there are no other superharmonic instabilities that occur before the one we have described above.Turning now to subharmonic modes we consider values of p ∈ (0, 1). Figure 4 shows the real and imaginary parts of the complex growth rate σ for several values of p.The two intermediate values p = 0.175 and p = 0.392 correspond to cases where there is a collision of eigenvalues at zero wave amplitude.Specifically, these occur where, with reference to (4. respectively.Many other similar examples can be readily constructed.These double, purely imaginary, eigenvalues split apart for c < 1 as the wave amplitude increases.Since there is no corresponding instability (the real part of σ remains zero as they split) we surmise that these two eigenvalues have the same Krein signature.We refer here to the result due to [29] that a necessary condition for instability is that colliding eigenvalues should have opposite Krein signature.The Krein signature, sK , is confirmed by direct computation: where s represents ±. (See [27] and [30] for similar calculations of the Krein signature for capillary waves.) Figure 5  , which have negative signature.Notably the double eigenvalues seen in Figure 4 for p = 0.175 and p = 0.392 have the same Krein signature and so cannot produce instability.The instability that occurs in the various panels in Figure 4 arise from the collision at finite amplitude of two modes of opposite Krein signature.For example, the instability that sets in at c = 0.974 for p = 0.175 occurs via a collision of the modes σ positive signature in each panel in Figure 4, we know that there cannot be any other collisions at higher frequencies than those shown in the panels which produce instability.Figures 6 and 7 show the eigenvalue spectrum in the complex plane for a number of different values of c. Evidently subharmonic instability occurs at any value of the wave speed.The structure of the spectrum changes in a highly intricate manner as the travelling-wave branch is followed.Just after the bifurcation point, at c = 0.99, we see two figure-of-eight type structures superposed on each other.This may be contrasted with what is observed for pure gravity waves in infinite depth, where a solitary figure-of-eight is seen for small wave amplitude [19].As c decreases the spectrum evidently undergoes a sequence of changes demonstrating a complex and intricate structure.Note that the computations for the various panels in Figures 6 and 7 required an increasing number of collocation points up to N = 2048 for c = 0.035.

Summary
We have constructed a branch of periodic travelling-wave solutions for hydroelastic waves on water of infinite depth and discussed the stability of these waves to small amplitude perturbations.For superharmonic perturbations we demonstrated that there occurs a collision of eigenvalues that corresponds to the onset of superharmonic instability, but not to the appearance of a bifurcating solution branch.We also examined higher frequency superharmonic modes and tracked the various branches.The absence of further collisions of eigenvalues, either at zero or away from zero indicates that no further superharmonic instabilities occur prior to the aforementioned one.This contrasts with the findings of [29] for gravity waves on water of infinite depth.In the latter case 'high' frequency collisions occur before the one at zero.[27] also found eigenvalue collisions for capillary waves on water of infinite depth but in this case the relevant eigenvalues have the same Krein signature and hence, according to the theory of [29], the collisions do not produce instability .
Additionally we have allowed for the presence of subharmonic perturbations and found that the waves are always subharmonically unstable for any amplitude along the travellingwave branch.Moreover we have described the eigenspectrum for different wave speeds and demonstrated that its structure becomes highly intricate, including the presence of figure-ofeight type patterns as are seen in traditional gravity-wave calculations.Appendix A. The form of Pe (ξ).
We note that Pe where k = 2π/λ is the fundamental wavenumber, and c1 will be determined at the next order.At O(ϵ 2 ) we find ϕ2xx + ϕ2yy = 0 , for y < 0 , (C.

Figure 1 :
Figure 1: Travelling-wave solutions: (left, upper) the bifurcation curve showing wave height H ≡ max(Y ) − min(Y ) versus the wave speed c; (left, lower) wave profiles corresponding to the filled circle symbols in the upper panel.The red dots show the limiting static profile corresponding to c = 0; (right) the solution curve for the elastic energy Ee, shown with a solid line, together with the Stokes wave expansion approximation (3.11), shown with a broken line.

Figure 2 :
Figure 2: Superharmonic case, p = 0. (a) The imaginary part of the growth rate Im(σ) plotted against the wave speed c.(b) Close-up near to c = 0.415 showing Re(σ) and Im(σ) with filled circles and empty diamonds, respectively.

Figure 3 :
Figure 3: (a) The wave energy E shown against wave speed c, shown with a solid line, with the limiting static wave (c = 0) energy E = 1.827 shown with a dotted line.The local maximum occurs at (c, E) = (0.415, 2.209) shown with a solid circle.(b) The Gram determinant G as defined in (5.1) shown against c.
(a) we show how the imaginary part of σ varies as the wave speed c is decreased from unity.The labels in the figure indicate branches which start at c = 1 at σ (p,j) ±

Figure 4 :
Figure 4: Subharmonic instability at (from left to right) p = 0.175, p = 0.2, p = 0.392, and p = 0.5.Each upper panel shows the imaginary part of σ corresponding to the real part of σ shown in the lower panel.
shows the Krein signature for a range of mode frequencies over the range p ∈ [0, 1].All of the modes have positive Krein signature except for those corresponding to σ

Figure 5 :
Figure 5: The imaginary part of the growth rate and Krein signature.Eigenvalue collisions occur at the crossing of the lines: two of these are illustrated by solid disks and correspond to the second and third panels in figure 4. The solid lines indicate modes with positive Krein signature, s K = 1, and the broken lines indicate modes with negative Krein signature, s K = −1.The blue broken line indicates the mode σ (p,0) + and the red broken line indicates the mode σ (p,−1) − .

Figure 6 :
Figure 6: The eigenspectrum in the complex σ plane for a number of different values of the wave speed c.In each case the plot focuses on the spectrum near the origin.Purely imaginary eigenvalues also exist but may lie beyond the range shown.

Figure 7 :
Figure 7: The eigenspectrum in the complex σ plane for a number of different values of the wave speed c.In each case the plot focuses on the spectrum near the origin.Purely imaginary eigenvalues also exist but may lie beyond the range shown.