Explicit equations for two-dimensional water waves with constant vorticity

TL;DR

This paper derives exact equations for 2D water waves with constant vorticity over uneven bottoms and introduces a highly accurate numerical method for solving these equations.

Contribution

It provides a novel, compact formulation of water wave equations with constant vorticity and develops an efficient numerical approach.

Findings

Exact conformal equations for water waves with vorticity

High-accuracy numerical method demonstrated

Applicable to arbitrary bottom profiles

Abstract

Governing equations for two-dimensional inviscid free-surface flows with constant vorticity over arbitrary non-uniform bottom profile are presented in exact and compact form using conformal variables. An efficient and very accurate numerical method for this problem is developed.

Equations

$${\phi }_t-\Omega \theta +\Omega y{\phi }_x+({\phi }_x^2+{\phi }_y^2)/2+gy+p/{\rho }_{\ast }=const,$$

$$X^{(s)}(u,t)+i{Y}^{(s)}(u,t)\equiv Z^{(s)}(u,t)=z(u+i,t),$$

$$X^{(b)}(u,t)+i{Y}^{(b)}(u,t)\equiv Z^{(b)}(u,t)=z(u,t).$$

$$\varphi (u+i,t)\equiv {\Phi }^{(s)}(u,t),\varphi (u,t)\equiv {\Phi }^{(b)}(u,t).$$

$${\Phi }^{(s)}(u,t)=e^{-\hat{k}}{\Phi }^{(b)}(u,t),$$

$$V_{(x)}-i{V}_{(y)}$$

$$-\text{mbox}Im\left(Z_t^{(s)}{\bar{Z}}_u^{(s)}\right)$$

$$0$$

$$\text{mbox}Re\left({\Phi }_t^{(s)}-{\Phi }_u^{(s)}{Z}_t^{(s)}/Z_u^{(s)}\right)+|{\Phi }_u^{(s)}/Z_u^{(s)}{|}^2/2+g\text{mbox}Im{Z}^{(s)}$$

$$-\Omega \text{mbox}Im{\Phi }^{(s)}+\Omega \text{mbox}Im{Z}^{(s)}\text{mbox}Re\left({\Phi }_u^{(s)}/Z_u^{(s)}\right)=0.$$

$${\Phi }^{(b)}=\hat{S}\psi -i(\Omega /2)(1-i\hat{R})[Y^{(b)}{]}^2,$$

$${\Phi }^{(s)}=(1+i\hat{R})\psi -i(\Omega /2)\hat{S}[Y^{(b)}{]}^2.$$

$$\zeta (w,t)=\int \frac{a_k(t)}{\cosh (k)}{e}^{ikw}\frac{dk}{2\pi },a_{-k}={\bar{a}}_k,$$

$$Z^{(b)}(u,t)=Z(\hat{S}a(u,t)),$$

$$\zeta (u+i,t)\equiv \xi (u,t)=(1+i\hat{R})a(u,t),$$

$$Q\equiv \frac{\{\hat{R}\psi +(\Omega /2)([\text{mbox}ImZ(\xi ){]}^2-\hat{S}[\text{mbox}ImZ(\hat{S}a){]}^2){\}}_u}{|Z_{\xi }(\xi ){\xi }_u{|}^2}.$$

$$a_t=-\text{mbox}Re[{\xi }_u(\hat{T}+i)Q].$$

$${\psi }_t$$

$$a(\vartheta ,t)$$

$$\xi (\vartheta ,t)$$

$$\psi (\vartheta ,t)$$

$${\rho }_t$$

$${\psi }_t$$

$$\Phi$$

$$\mathsf{Q}$$

$$\dot{\alpha }(t)=-\frac{1}{2\pi }{\int }_0^{2\pi }\mathsf{Q}(\vartheta )d\vartheta .$$

$$E$$

$$\mathcal{B}(q)=q-i\Delta \text{mbox}Ln\left[\left(i\sin q+\sqrt{ϵ+{\cos }^{2}q}\right)(1+ϵ{)}^{-1/2}\right],$$

Full text

Most theories for surface water waves have been developed under assumption of irrotational flows (see, for example, Refs.Zakharov67 ; Z1999 ; LZ2005 , and references therein). The reason is that irrotational incompressible flows are completely determined by the boundary distribution of the velocity potential through solution of the Laplace equation, and this fact effectively reduces spatial dimensionality of the problem. However, in many real situations waves propagate on shear currents, and over non-uniform depth. Interaction between waves and currents is important in many aspects, for example, as a mechanism of giant wave formation Peregrine ; Kharif-Pelinovsky ; LP2006 . The fully nonlinear problem of vortical flows under free surface is very complicated. Various simplifications are assumed in theoretical studies, as weak nonlinearity, mild slope, long wave length, and so on (see, for examples, Refs.B1962 ; FJ1970 ; Y1972 ; J1990 ; Kirby ; B2007 ). The only possibility for non-potential flows, when dimensionality is reduced, is the two-dimensional (2D) problem with constant vorticity, since in this case all perturbations are irrotational. In particular, relatively long solitary waves in water of finite uniform depth with constant vorticity were calculated analytically Miroshnikov2002 ; Choi2003 . There are also fully nonlinear numerical results for steady flows obtained by a boundary integral method (see Refs.SilvaPeregrine1988 ; PullinGrimshaw1988 ; VdB1994-1996 , and references therein). At the same time, unsteady rotational waves have not been extensively simulated. The present work is intended to fill this gap. Here for the first time an exact and compact formulation of the fully nonlinear problem of waves with constant vorticity is presented in terms of so called conformal variables. The conformal variables were introduced earlier and successfully used to describe purely potential flows on a constant depth or on infinitely deep water DKSZ96 ; DZK96 ; ZDV2002 ; DZ2005Pisma ; ZDP2006 . Later, the description has been generalized by the present author to the case of potential flows over arbitrary non-uniform and time-dependent bottom profile R2004PRE ; R2005PLA . One of the main advantages of equations in conformal variables is an easy numerical implementation with fast Fourier transform (FFT) subroutines. In this work, after derivation of exact evolutionary equations, some illustrative numerical results will be presented. Here, for simplicity, a non-uniform bottom profile does not depend on time, though a generalization to non-static bottom is straightforward.

Derivation of exact equations. In what follows, style and notations will be the same as in Ref.R2005PLA . We consider here a 2D incompressible inviscid non-stationary flow in $(x,y)$-plane, bounded by a given bottom profile from below and by an unknown free surface from above. The flow is rotational, with a constant vorticity field, $\partial_{x}V_{(y)}-\partial_{y}V_{(x)}=-\Omega=const.$ The velocity field thus can be represented as follows, ${\bf V}=(\Omega y+\varphi_{x},\quad\varphi_{y}),$ where the potential $\varphi(x,y,t)$ satisfies the Laplace equation $\varphi_{xx}+\varphi_{yy}=0$. Let us also introduce a harmonically conjugate function $\theta(x,y,t)$: $\varphi_{x}=\theta_{y},\quad\varphi_{y}=-\theta_{x}.$ Then it is easy to check that two components of the vector Euler equation are equivalent to a single scalar equation, analogous to the Bernoulli equation,

[TABLE]

where $g$ is the gravity acceleration, $\rho_{*}$ is a constant density of the fluid, and $p(x,y,t)$ is the pressure. It is Eq.(1) that makes possible reduction of dimensionality of the problem, in the same manner as for purely potential flows.

Important point is that complex combination $\tilde{\phi}(z,t)=\varphi(x,y,t)+i\theta(x,y,t)$ is an analytic function of the complex argument $z=x+iy$. Analyticity is preserved under conformal coordinate transforms $z=z(w,t)$, where $w=u+iv$ is a new complex variable. We choose an analytic function $z(w,t)$ in such a manner that $w=u$ at the bottom and $w=u+i$ at the free surface. Shape of the surface will be given in a parametric form,

[TABLE]

The bottom profile will be determined by

[TABLE]

Thus, we have an analytic function $\phi(w,t)=\tilde{\phi}(z(w,t),t)$ defined in the stripe $0\leq v\leq 1$ in the $(u,v)$-plane. Let us designate boundary values of this function as written below,

[TABLE]

Since $\Phi^{(s)}(u,t)$ and $\Phi^{(b)}(u,t)$ are values of the same analytic function at points $u$ and $u+i$, they are related to each other by a linear transform (see R2005PLA ),

[TABLE]

with $e^{-\hat{k}}\equiv\exp(i\hat{\partial}_{u})$. That means $\Phi^{(s)}_{k}(t)=e^{-k}\Phi^{(b)}_{k}(t)$ for the corresponding Fourier images.

The velocity components are determined by the following relations

[TABLE]

Now we are going to write equations of motion in the conformal variables. First, we have two kinematic conditions which in our case take form

[TABLE]

where $\bar{Z}$ denotes complex conjugate value, and the subscripts denote the corresponding partial derivatives. At free surface the pressure is constant (we neglect here surface tension $\sigma$, otherwise $p^{(s)}=\sigma\kappa+const$, where $\kappa$ is the surface curvature). Therefore from Eq.(1) we have the dynamic boundary condition in conformal variables,

[TABLE]

Taking into account Eq.(8), it is convenient to represent $\Phi^{(b)}(u,t)$ in the form

[TABLE]

where $\psi(u,t)$ is some unknown real function, and the linear operators $\hat{S}$ and $\hat{R}$ are diagonal in Fourier representation: $S_{k}={1}/{\cosh(k)},\quad R_{k}=i\tanh(k).$ It is essential that when $\hat{S}$ or $\hat{R}$ acts on a purely real function, the result is also real. We will also need operator $T_{k}=-i\coth(k)=R_{k}^{-1}$. Since $e^{-\hat{k}}\hat{S}=(1+i\hat{R})$ and $e^{-\hat{k}}(1-i\hat{R})=\hat{S}$, we have from Eqs.(5) and (10) the following formula for $\Phi^{(s)}(u,t)$,

[TABLE]

Now we should take into account that the function $z(w,t)$ can be represented as a composition of two functions (see R2004PRE ; R2005PLA ), that is $z(w,t)=Z(\zeta(w,t)),$ where a known analytic function $Z(\zeta)=X(\zeta)+iY(\zeta)$ determines bottom shape. The conformal mapping $Z(\zeta)$ does not have any singularities within a sufficiently wide horizontal stripe above the real axis in $\zeta$-plane. An intermediate analytic function $\zeta(w,t)$ takes real values at the real axis, and therefore

[TABLE]

where $a_{k}(t)$ is Fourier transform of a real function $a(u,t)$. On the bottom $\zeta(u,t)=\hat{S}a(u,t)$, therefore

[TABLE]

At the free surface we have relations

[TABLE]

and $Z^{(s)}=Z(\xi),\quad Z^{(s)}_{u}=Z_{\xi}(\xi)\xi_{u},\quad Z^{(s)}_{t}=Z_{\xi}(\xi)\xi_{t}.$

Thus we have in our system two unknown real functions, $\psi(u,t)$ and $a(u,t)$. All the other quantities are expressed through these two. Our purpose now is to derive equations determining time derivatives $\psi_{t}$ and $a_{t}$. To do this, we divide Eq.(7) by $|Z^{(s)}_{u}|^{2}$ and obtain that $\mbox{Im}({\xi_{t}}/{\xi_{u}})=-Q$, where

[TABLE]

Since $\xi_{t}/\xi_{u}=\zeta_{t}(w,t)/\zeta_{w}(w,t)|_{w=u+i}$, there exists a relation between the real and imaginary parts: $\mbox{Im}(\xi_{t}/\xi_{u})=\hat{R}\,\mbox{Re}(\xi_{t}/\xi_{u})$, so $\mbox{Im}({\xi_{t}}/{\xi_{u}})=-Q$ means $\xi_{t}=-\xi_{u}(\hat{T}+i)Q$ and it gives us equation determining $a_{t}$,

[TABLE]

After that, Eqs.(9) and (11) allow us to express $\psi_{t}$:

[TABLE]

Now, exact and explicit evolutionary equations have been derived. In the next section it is explained how one can numerically simulate them with a high accuracy.

Numerical method and example. Let us consider the case when a bottom profile is periodic on $x$. Obviously, there exist solutions with velocity field and surface elevation having the same spatial period $L$. Without loss of generality, the potential $\varphi(x,y,t)$ can be assumed periodic (the part of $\varphi$ which is proportional to $x$ and corresponds to a constant horizontal velocity, can be excluded by a redefinition of $y$ zero level). Making a proper choice for the length and time scales, we may write $g=1$, and $Z(\zeta+2\pi)=2\pi+Z(\zeta)$. Direct substitution into Eqs.(16) and (17) shows the solutions are $2\pi$-periodic on the variable $\vartheta=u\alpha(t)$, where a real function $\alpha(t)$ depends on time in a non-trivial manner in order to cancel non-periodic terms in Eqs.(16) and (17). The non-periodic terms take place because the operator $\hat{T}$ is singular at small $k$, and its action on a constant function $Q_{0}$ is non-periodic: $\hat{T}Q_{0}=Q_{0}u$ (see Refs.R2004PRE ; R2005PLA ). So, we can write for rescaled dimensionless quantities

[TABLE]

where $\rho_{m}(t)$ and $\psi_{m}(t)$ are Fourier coefficients of $2\pi$-periodic real functions $\rho(\vartheta,t)$ and $\psi(\vartheta,t)$. As a result, in variables $(\vartheta,t)$ equations for $\rho_{t}$ and $\psi_{t}$ look similar to Eqs.(16) and (17), but all the $u$-derivatives should be replaced by $\vartheta$-derivatives, and operators $\hat{R}$, $\hat{S}$, and $\hat{T}$ should be everywhere replaced by new operators $\hat{\mathsf{R}}_{\alpha}$, $\hat{\mathsf{S}}_{\alpha}$, and $\hat{\mathsf{T}}_{\alpha}$ respectively:

[TABLE]

where $\xi=\vartheta+i\alpha+(1+i\hat{\mathsf{R}}_{\alpha})\rho$, and

[TABLE]

These new operators are diagonal in the discrete Fourier representation: ${\mathsf{R}}_{\alpha}(m)=i\tanh(\alpha m)$, ${\mathsf{S}}_{\alpha}(m)=1/\cosh(\alpha m)$, and ${\mathsf{T}}_{\alpha}(m)=-i\coth(\alpha m)$ for $m\not=0$, ${\mathsf{T}}_{\alpha}(0)=0$. The system of equations is closed by the following condition for $\dot{\alpha}(t)$, which ensures cancellation of the non-periodic terms in Eqs.(16) and (17),

[TABLE]

The above system of equations has two apparent integrals of motion, namely the area $A$ occupied by fluid, $A=A_{0}+\int Y(\xi)\,\mbox{Re}[Z_{\xi}(\xi)\xi_{\vartheta}]\,d\vartheta,$ and the total energy $E$ (kinetic energy plus potential energy in the gravitational field),

[TABLE]

where $A_{0}$ ans $E_{0}$ are constant, and all the integrals are in the limits from [math] to $2\pi$.

Eqs.(21)-(25) are easy for numerical simulation if the function $Z(\zeta)$ is given by a simple formula as it takes place for many interesting bottom profiles. The numerical method employed here is naturally based on the discrete Fourier representation, since all the linear operators in the equations are efficiently computed with modern FFT subroutines in $m$-representation, while all the nonlinear operations are simple in $\vartheta$-representation. As primary dynamical variables, the quantities $\alpha(t)$, $\rho_{m}(t)$, and $\psi_{m}(t)$ are taken, with $0\leq m<M$ (for negative $m$ the relations $\rho_{-m}=\bar{\rho}_{m}$ and $\psi_{-m}=\bar{\psi}_{m}$ are used). After each step of a Runge-Kutta 4-th order procedure, only spectral components with $|m|<M_{eff}$ are kept, where $M_{eff}\approx(1/4)N$, $M\approx(3/8)N$, and $N=2^{12...19}$ is the size of arrays for the fast Fourier transform (during computations, $N$ is doubled several times as small-scale structures develop). As a result of the adaptive increasing of $N$, the right hand sides of Eqs.(21)-(22) can be computed with nearly the same numerical error $\delta_{0}<N10^{-18}$ as it is for the FFT subroutine using C-type double. Since the time step is decreased as $\tau\sim 1/N$ for the stability reasons, an error for the free surface position at $t\sim 1$ can be estimated as $\delta\lesssim N^{2}10^{-18}$. Practically, $A$ and $E$ are conserved up to $10$ decimal digits for most part of the evolution. In a final stage, the larger $N_{final}$ is used, the later time moment is when the high accuracy is lost.

Here an example is given which demonstrates potentialities of the method. Let the bottom profile be determined by formula $Z(\zeta)=iY_{0}+{\cal B}(\zeta-i\alpha_{0})$, where $Y_{0}=0.02\pi$, $\alpha_{0}=0.02\pi$, and

[TABLE]

with $\Delta=0.6$ and $\epsilon=0.02$. At $t=0$ we put $\alpha=\alpha_{0}$ and $\rho=0.05\tanh[15\sin(\vartheta-\gamma)]\exp[-2(1-\cos(\vartheta-\gamma))],$ where $\gamma=0.5$. In Fig.1-a, the bottom profile and the initial surface shape are presented normalized to the spatial period $L=100$ m and shifted appropriately in the vertical direction. Dimensionless vorticity $\Omega=2.0$ together with the parameter $Y_{0}$ give the rotational part of the horizontal velocity field $V^{(\Omega)}_{(x)}\approx 1.6(y+1)$ m/s. This corresponds to a backward flow along the bottom in the deeper regions. We choose $\psi(\vartheta,0)$ in such a way that at $t=0$ the normal component of the total velocity field at the free surface is zero: $\psi(\vartheta,0)=-({\Omega}/{2})\hat{\mathsf{T}}_{\alpha}[Y^{2}(\xi)-\hat{\mathsf{S}}_{\alpha}Y^{2}(\vartheta+\hat{\mathsf{S}}_{\alpha}\rho)].$ Some results of the computation are presented in Fig.1, where also a comparison is made to the case $\Omega=0$. The initial hump at the surface decays into two oppositely propagating solitary waves having different speeds and different profiles (crest of the right-propagating wave is more sharp; however, a maximum curvature is finite). In this example the right-propagating wave first meets the region of relatively shallow depth, where its crest becomes more and more steep, and finally the wave profile overturns, as it is seen in Fig.2. Velocity distribution along the overturning wave is shown in Fig.3. Another interesting phenomenon observed here is the wave blocking near $x\approx 60$ m, where an average horizontal flow velocity ($\sim 1.6$ m/s for $\Omega=2.0$) approaches speed of typical waves $v_{ph}\approx\sqrt{gh}$ ($h$ is the local depth), and therefore waves cannot enter the shallow region from the right. In different simulations, for $\Omega\gtrapprox 2.4$, the blocking was so strong that wave height near the point $x=60$ was comparable to the local depth (not shown). However, an extensive discussion of this phenomenon is not possible in this Brief Report.

Acknowledgments.These investigations were supported by RFBR Grant 06-01-00665, by the Program “Fundamental Problems of Nonlinear Dynamics” from the RAS Presidium, and by Grant “Leading Scientific Schools of Russia”.