2. Formulation of the problem

The flow domain consists of two semi-infinite regions, $\varOmega _+$ and $\varOmega _-$ , occupied by fluids of densities $\rho _+$ and $\rho _-$ ( $\rho _+ \gt \rho _-$ ), respectively. These regions are connected through the interface $\varGamma : y = \eta (x,t)$ , whose equilibrium position is at $y = 0$ , as shown in figure 1. A circular cylinder of radius $r$ starts moving in one of the layers without intersecting the interface, following the prescribed trajectory $(x_c(t),y_c(t))$ for $t \gt 0$ .

In subsequent derivations and discussion, we will use non-dimensional variables. The density of the lower layer $\rho _+$ , the initial submergence depth $h$ and the characteristic speed of the cylinder $u_0$ form a dimensionally independent set. Consequently, the quantities $h/u_0$ and $\rho u_0^2$ determine the units of time and pressure, respectively. The non-dimensional density of the upper layer, which characterises the density ratio of the two fluids, is denoted by $\rho = \rho _-/\rho _+$ . By definition, $0 \leqslant \rho \lt 1$ . For subsequent discussion, it will be useful to introduce the Atwood number, which characterises the density difference of the two layers:

(2.1) \begin{equation} a = \frac {\rho _+ - \rho _-}{\rho _+ + \rho _-} = \frac {1 - \rho }{1+\rho }. \end{equation}

In the special case, when $a = 1$ , we deal with the free surface flow of a one-layer fluid. The limiting case of $a = 0$ , when a two-layer fluid becomes homogeneous and the interface disappears, has no meaning in the current formulation.

The fluid velocity vector $\boldsymbol{U}_{\pm } = (U_{\pm },V_{\pm })$ and pressure $p_{\pm }$ in both layers should satisfy the irrotational Euler equations

(2.2) \begin{align} \cfrac {\partial U_{\pm }}{\partial x} + \cfrac {\partial V_{\pm }}{\partial y} = 0, \quad \cfrac {\partial U_{\pm }}{\partial y} - \cfrac {\partial V_{\pm }}{\partial x} &= 0, \quad (x,y) \in \varOmega _{\pm }, \end{align}

(2.3) \begin{align} \cfrac {\partial U_\pm }{\partial t} + U_\pm \cfrac {\partial U_\pm }{\partial x} + V_\pm \cfrac {\partial U_\pm }{\partial y} + \cfrac 1{\rho _{\ast }}\cfrac {\partial p_\pm }{\partial x} &= 0,\quad (x,y) \in \varOmega _\pm , \end{align}

(2.4) \begin{align} \cfrac {\partial V_\pm }{\partial t} + U_\pm \cfrac {\partial V_\pm }{\partial x} + V_\pm \cfrac {\partial V_\pm }{\partial y} + \cfrac 1{\rho _{\ast }}\cfrac {\partial p_\pm }{\partial y} &= -\lambda ,\!\!\!\!\quad (x,y) \in \varOmega _\pm . \end{align}

Here, the non-dimensional density parameter $\rho _{\ast }$ is equal to $\rho$ in the upper layer and to unity in the lower layer. The important physical parameter $\lambda = gh/u_0^2$ , defined with the gravity acceleration $g$ , is the square of the inverse Froude number. The fully nonlinear kinematic and dynamic boundary conditions at the interface $\varGamma (t): y = \eta (x,t)$ are

(2.5) \begin{align} \cfrac {\partial \eta }{\partial t} + U_{\pm } \cfrac {\partial \eta }{\partial x} &= V_{\pm }, \quad (x,y) \in \varGamma , \end{align}

(2.6) \begin{align} p_+ &= p_-, \quad (x,y) \in \varGamma . \end{align}

The boundary condition on the moving cylinder situated in the lower layer takes the form

(2.7) \begin{equation} (\boldsymbol{U}_+ - \boldsymbol{U}_c)\boldsymbol{\cdot }\boldsymbol n = 0, \quad (x,y) \in S, \end{equation}

where the vector $\boldsymbol n$ is the unit normal to the surface of the cylinder and the vector $\boldsymbol{U}_c(t) = (\dot x_c(t),\dot y_c(t))$ is the prescribed non-dimensional velocity of the cylinder. In the case of a cylinder in the upper layer, the subscript $_+$ in (2.7) must be changed to $_-$ . We assume that both fluids are at rest at infinity, as described by the condition

(2.8) \begin{equation} \eta \to 0, \quad \boldsymbol{U}_{\pm } \to 0, \quad |\boldsymbol x| \to \infty . \end{equation}

At the initial time $t = 0$ , the cylinder centre is placed at the point $(x_c(0),y_c(0) = (0,\pm 1)$ , where plus and minus signs correspond to the cylinder location in the lower and upper layer, respectively. The initial velocity field $\boldsymbol{U}_{\pm }(0) = (U_\pm (0),V_\pm (0))$ must satisfy the compatibility conditions

(2.9) \begin{align} \cfrac {\partial U_{\pm }}{\partial x} + \cfrac {\partial V_{\pm }}{\partial y} = 0, \quad \cfrac {\partial U_{\pm }}{\partial y} - \cfrac {\partial V_{\pm }}{\partial x} &= 0, \quad (x,y) \in \varOmega _{\pm }, \end{align}

(2.10) \begin{align} (\boldsymbol{U}_\pm (0) - \boldsymbol{U}_c(0))\boldsymbol{\cdot }\boldsymbol n &= 0, \quad (x,y) \in S. \end{align}

The initial data determine velocity circulation around the contour of the cylinder as

(2.11) \begin{equation} \gamma = \oint _S U_\pm (0)\, {\rm d}x + V_\pm (0)\, {\rm d}y, \end{equation}

which remains constant at $t \gt 0$ due to the Kelvin theorem. Here, the index $+$ or $-$ is chosen, depending on the location of the cylinder below or above the interface, respectively. The non-zero velocity $U_{\pm }(0)$ at the initial time corresponds to the problem of an impulsive start of the body starting to move with non-zero speed. In this study, however, we consider the smooth start of the cylinder with constant acceleration from rest, corresponding to the initial data $\boldsymbol{U}_{\pm }(0) = (0,0)$ and $\boldsymbol{U}_c(0) = (0,0)$ , for which the compatibility conditions are automatically satisfied.

3. Reduction to the boundary

A challenging aspect of water wave problems in a nonlinear formulation is that the boundary conditions are enforced on an unknown boundary, the position of which is part of the solution. The mathematical problem considered here allows for a closed formulation in the form of a system of integro-differential equations valid on the moving boundary $\varGamma (t): y = \eta (x,t)$ . The tangential and normal fluid velocities on the interface in both layers are determined as projections of the fluid-particle velocity vector $\boldsymbol{U}_{\pm } = (U_\pm , V_\pm )$ onto the tangent $(1,\eta _x)$ and normal $(-\eta _x,1)$ vectors at the corresponding point on the interface $y = \eta (x,t)$ :

(3.1) \begin{equation} u_{\pm } = U_{\pm } + \eta _x V_{\pm }, \quad v_{\pm } = V_{\pm } - \eta _x U_{\pm }, \quad y = \eta (x,t). \end{equation}

We assume that two fluids do not mix and therefore the components of fluid-particle velocities in two layers normal to the interface $y = \eta (x,t)$ must be equal: $v_+ = v_-$ . To minimise the use of subscripts, we adopt the notation $u = u_+$ , $w = u_-$ , $v = v_+ = v_-$ .

In the new variables, the kinematic boundary condition (2.5) reduces to a single equation

(3.2) \begin{equation} \eta _t = v. \end{equation}

Projecting the momentum equations (2.3) and (2.4) in both fluids onto the tangent line to the interface $\boldsymbol \tau = (1,\eta _x)$ , and using the dynamic boundary condition (2.6) allow us to exclude fluid pressure $p_{\pm }$ from the problem formulation:

(3.3) \begin{equation} u_t-\rho w_t + \frac 12\cfrac {\partial }{\partial x} \left (\cfrac {u^2 - \rho w^2 - 2\eta _x v(u - \rho w) - (1-\rho ) v^2}{1+\eta _x^2}\right ) + (1-\rho )\lambda \eta _x = 0. \end{equation}

According to the method developed by Ovsyannikov et al. (Reference Ovsyannikov, Makarenko, Nalimov, Liapidevskii, Plotnikov, Sturova, Bukreev and Vladimirov1985), differential equations (3.2) and (3.3) must be complemented by two boundary integral equations following from the Cauchy integral formula for the complex velocity function $F_{\pm }(z) = U_\pm - iV_\pm$ , which is analytic in $\varOmega _{\pm }$ . The analyticity follows from (2.2), equivalent to the Cauchy–Riemann conditions for the velocity components $U_\pm$ and $V_\pm$ . When the cylinder is located in the lower layer, the boundary integral equation for the lower layer $\varOmega _+$ has the form

(3.4) \begin{equation} 2\pi i F_+(z,t) = \int \limits _{\varGamma (t)} \!\!\frac {F_+(\zeta ,t){\rm d}\zeta }{\zeta - z} + \frac {r^2}{(z - z_c(t))^2}\overline {\int \limits _{\varGamma (t)}\!\!\frac {F_+(\zeta ,t){\rm d}\zeta }{\zeta - z_{\ast }(t)}} + \frac {\gamma }{z - z_c(t)} + \frac {2\pi ir^2\dot z_c(t)}{(z-z_c(t))^2}, \end{equation}

where $z_{\ast }(t) = z_c(t) + r^2/\overline {z-z_c}$ is the inverse image of the complex point $z$ with respect to the circle centred at the point $z_c(t) = x_c(t) + iy_c(t)$ . The constant $\gamma$ is the velocity circulation around the cylinder, as defined by (2.11). The boundary integral equation (3.4), which contains integrals over the free boundary $\varGamma (t)$ only, was first derived by Ovsyannikov et al. (Reference Ovsyannikov, Makarenko, Nalimov, Liapidevskii, Plotnikov, Sturova, Bukreev and Vladimirov1985) using the Cauchy integral formula for the analytic function $F_+(z,t)$ and the Milne–Thomson circle theorem. The derivation of this integral equation has been reproduced multiple times by Makarenko (Reference Makarenko2003), Kostikov & Makarenko (Reference Kostikov and Makarenko2018) and Pardo & Nedic (Reference Pardo and Nedic2021). A combination of real and imaginary parts of (3.4), evaluated at the interface, gives the real-valued version of the boundary integral equation. This can be written in operator form as follows:

(3.5) \begin{equation} [I + A(\eta )+ r^2A^r(\eta )] v - [B(\eta ) + r^2B^r(\eta )] u = v^d. \end{equation}

Here, the integral operators $A$ and $B$ correspond to the problem of waves on the interface without a cylinder:

(3.6) \begin{align} A(\eta )v(x) &= \frac 1\pi \mathop {\int \!\!\!\!\!\!\boldsymbol{-}}\limits _{-\infty }^\infty \frac {\eta (x) - \eta (s) - (x-s)\eta '(x)}{(x-s)^2 + [\eta (x) - \eta (s)]^2} v(s){\rm d}s, \end{align}

(3.7) \begin{align} B(\eta )u(x) &= \frac 1\pi \mathop {\int \!\!\!\!\!\!\boldsymbol{-}}\limits _{-\infty }^\infty \frac {x - s + [\eta (x)-\eta (s)]\eta '(x)}{(x-s)^2 + [\eta (x) - \eta (s)]^2} u(s){\rm d}s, \end{align}

and the integral operators with superscripts $A^r$ and $B^r$ are responsible for interaction of the cylinder with the interface:

(3.8) \begin{align} A^r(\eta ) v (x) &= \frac 1\pi \mathop {\int \!\!\!\!\!\!\boldsymbol{-}}\limits _{-\infty }^\infty \frac {\left[\eta (s) - y_c - r^2P(x)\right]Q'(x) - \left[s - x_c - r^2Q(x)\right]P'(x)}{\left[\eta (s) - y_c - r^2P(x)\right]^2 + \left[s - x_c - r^2Q(x)\right]^2} v(s){\rm d}s, \end{align}

(3.9) \begin{align} B^r(\eta )u(x) &= \frac 1\pi \mathop {\int \!\!\!\!\!\!\boldsymbol{-}}\limits _{-\infty }^\infty \frac {\left[\eta (s) - y_c - r^2P(x)\right]P'(x) + \left[s - x_c - r^2Q(x)\right]Q'(x)}{\left[\eta (s) - y_c - r^2P(x)\right]^2 + \left[s - x_c - r^2Q(x)\right]^2} u(s)\,{\rm d}s. \end{align}

It should be noted that the above-mentioned dashed integrals are to be considered in the sense of the Cauchy principal value to account for the discontinuities of the kernel functions. Here, the functions $P(x)$ and $Q(x)$ are the generalised versions of the Poisson kernels:

(3.10) \begin{equation} P(x) = \frac {\eta (x) - y_c}{(x-x_c)^2 + (\eta (x) - y_c)^2}, \quad Q(x) = \frac {x-x_c}{(x-x_c)^2 + (\eta (x)-y_c)^2}, \end{equation}

which depend on the interface shape $\eta (x,t)$ and the coordinates of the cylinder $(x_c(t),y_c(t))$ . The term $v^d$ is the normal component of the fluid velocity induced at the interface by a dipole of unit strength:

(3.11) \begin{equation} v^d(x) = 2r^2\left (\dot y_c Q'(x) - \dot x_c P'(x)\right )\!. \end{equation}

Here and in what follows, a prime denotes the derivative of an ordinary function with respect to the spatial coordinate $x$ and an upper dot signifies the derivative of a one-variable function with respect to the time variable $t$ .

The analogue of the integral equation for the upper layer without a body can be obtained from (3.5) by omitting the terms $A^r v$ , $B^ru$ and $v^d$ , which are responsible for the interaction of the cylinder with the interface, and by inversing the signs preceding the integrals due to the reversed direction of integration:

(3.12) \begin{equation} [I - A(\eta )]v + B(\eta )w = 0. \end{equation}

The system of equations (3.2), (3.3), (3.5) and (3.12) forms a complete set of equations for the unknown interface $\eta$ , the normal velocity $v$ , and the tangential velocities $u$ and $w$ on the interface. This system fully defines the problem of interface waves generated by the arbitrary motion of a cylinder in the lower layer.

When the cylinder moves in the upper layer, the differential equations (3.2) and (3.3) on the interface remain the same, as they carry no information about the obstacle. However, the integral equations (3.5) and (3.12) exchange their roles due to the changed position of the cylinder relative to the interface. Consequently, the following integral equations must be satisfied on the interface:

(3.13) \begin{align} \left[I - A(\eta )- r^2A^r(\eta )\right] v + \left[B(\eta ) + r^2B^r(\eta )\right] w &= v^d, \\[-28pt] \nonumber \end{align}

(3.14) \begin{align} [I + A(\eta )]v - B(\eta )u &= 0. \end{align}

Equations (3.13) and (3.14) can be obtained from (3.5) and (3.12) by interchanging the functions $u$ and $w$ , and inverting the signs of the integral operators. The dependence on time is omitted in the above-mentioned formulae, as it appears only as a parameter. In the next section, we will consider a cylinder moving in the lower layer below the interface for clarity of mathematical expressions, keeping in mind that the case of a cylinder in the upper layer can be treated similarly.

4. Small-time asymptotic solution

The boundary integro-differential equations (3.2), (3.3), (3.5) and (3.12) form a nonlinear evolution system, which is difficult to solve analytically. The integral operators in (3.5) and (3.12), while acting linearly on the unknown functions $u$ , $w$ and $v$ , have kernels that depend nonlinearly on the interface displacement $\eta$ . This system can be solved for the initial stage of motion using a recursive scheme derived by expanding the involved variables into small-time series.

We consider the unsteady flow at the interface, induced by the smooth start of the cylinder, which is characterised by zero initial conditions:

(4.1) \begin{equation} u(x,0) = w(x,0) = v(x,0) = \eta (x,0) = 0. \end{equation}

The trajectory of the cylinder, moving in the lower layer of a two-layered fluid with constant acceleration from rest at an angle $\theta$ with respect to the horizontal axis, is defined by

(4.2) \begin{equation} (x_c(t), y_c(t)) = \left(t^2\cos \theta , -1+t^2\sin \theta \right)\!. \end{equation}

Kinematic boundary condition (3.2) and zero initial conditions (4.1) dictate the following structure of the power series of unknown functions with respect to the time variable $t$ :

(4.3) \begin{align} \eta (x,t) &= t^2 \eta _2(x) + t^3 \eta _3(x) + \cdots , \qquad u(x,t) = tu_1(x) + t^2 u_2(x) + \cdots , \end{align}

(4.4) \begin{align} w(x,t) &= tw_1(x) + t^2w_2(x) + \cdots , \qquad v(x,t) = tv_1(x) + t^2 v_2(x) + \cdots. \end{align}

Substituting series (4.3)–(4.4) into (3.2)–(3.3) and grouping the terms of equal order, we obtain the following set of recursive relations between terms:

(4.5) \begin{align} \eta _{n+1}(x) &= \frac 1{n+1}v_n(x), \\[-28pt] \nonumber \end{align}

(4.6) \begin{align} u_{n}(x) - \rho w_{n}(x) &= \phi _{n}(x), \end{align}

where function

(4.7) \begin{align} \phi _n(x) &= - \frac 1{2n}\frac {\partial }{\partial x}\left [ \sum \limits _{j=1}^{n-1} (u_j u_{n-j-1} - \rho w_j w_{n-j-1}) \right. \nonumber \\&\quad \left. -\, 2\sum \limits _{j+k+l=0}^{n-2} \eta _j(u_k - \rho w_k)v_l - (1-\rho )\sum \limits _{j=1}^{n-1} v_j v_{n-j-1} + 2(1-\rho )\lambda \eta _{n-1} \right ]\!, \end{align}

inherits the structure of linear and nonlinear terms from boundary differential equation (3.3).

The leading-order terms in the time expansions of the generalised Poisson kernels $P(x)$ and $Q(x)$ , defined by (3.10), can be found explicitly as

(4.8) \begin{align} P(\eta ) &= p(x) + \Bigl [ -p'(x)\cos \theta + q'(x)\sin \theta - q'(x)\eta _2(x)\Bigr ]t^2 + \cdots , \end{align}

(4.9) \begin{align} Q(\eta ) &= q(x) + \Bigl [ -p'(x)\sin \theta - q'(x)\cos \theta + p'(x)\eta _2(x)\Bigr ]t^2 + \cdots . \end{align}

The pair of functions

(4.10) \begin{equation} p(x) = \frac 1{1+x^2}, \quad q(x) = \frac {x}{1+x^2} \end{equation}

forms a conjugate pair with respect to the Hilbert transform, which is defined by the convolution with Cauchy’s kernel as

(4.11) \begin{equation} Hu(x) = \frac 1\pi \mathop {\int \!\!\!\!\!\!\boldsymbol{-}}\limits _{-\infty }^\infty \frac {u(s)}{x-s}{\rm d}s, \end{equation}

meaning that $Hp(x) = q(x)$ and $Hq(x) = -p(x)$ . The functions $p(x)$ and $q(x)$ , referred to here as Poisson kernels, will play a crucial role in constructing the solution to the problem. These functions decay at infinity and, thus, the Hilbert transform (4.11) is well-defined for them. The property of the Hilbert transform, that it coincides with its inverse up to a sign $(-H = H^{-1})$ , will be used extensively in our study.

The time series for the integral operators, defined by (3.6)–(3.9), have the form

(4.12) \begin{align} A(\eta ) &= t^2 A_2 + t^3 A_3 + \cdots , \quad B(\eta ) = B_0 + t^4 B_4 + \cdots , \end{align}

(4.13) \begin{align} A^r(\eta ) &= A^r_0+ t^2 A^r_2 + \cdots , \quad B^r(\eta ) = B^r_0 + t^2 B^r_2 + \cdots , \end{align}

where the missing terms are accounted for by the form of the expansion (4.3) for $\eta$ and the smooth acceleration of the cylinder with zero initial speed. The coefficients of the leading-order in expansions (4.12)–(4.13), expressed via the integrals

(4.14) \begin{align} A_2v(x) &= \frac 1\pi \mathop {\int \!\!\!\!\!\!\boldsymbol{-}}\limits _{-\infty }^\infty \frac {\eta _2(x) - \eta _2(s)}{(x-s)^2}v(s)\,{\rm d}s - \frac 1\pi \eta ^{\prime}_2(x)\mathop {\int \!\!\!\!\!\!\boldsymbol{-}}\limits _{-\infty }^\infty \frac {v(s)}{x-s}{\rm d}s, \end{align}

(4.15) \begin{align} B_0u(x) &= Hu(x), \end{align}

(4.16) \begin{align} A^r_0 v(x) &= \frac 1\pi \mathop {\int \!\!\!\!\!\!\boldsymbol{-}}\limits _{-\infty }^\infty \frac {\left(1-r^2p(x)\right)q'(x) - \left(s - r^2q(x)\right)p'(x)}{\left(1-r^2p(x)\right)^2 + \left(s-r^2q(x)\right)^2}v(s){\rm d}s, \end{align}

(4.17) \begin{align} B^r_0 u(x) &= \frac 1\pi \mathop {\int \!\!\!\!\!\!\boldsymbol{-}}\limits _{-\infty }^\infty \frac {\left(1-r^2p(x)\right)p'(x) + \left(s - r^2q(x)\right)q'(x)}{\left(1-r^2p(x)\right)^2 + \left(s-r^2q(x)\right)^2}u(s)\,{\rm d}s \end{align}

are the crucial integral components, contributing to the small-time asymptotic solution.

The normal fluid velocity on the interface $v^d$ induced by a moving dipole has the power series expansion

(4.18) \begin{equation} v^d(\eta ) = t v^d_1 + t^3 v^d_3 + \cdots \end{equation}

with coefficients

(4.19) \begin{align} v^d_1(x)& = 4r^2(q'(x)\sin \theta - p'(x)\cos \theta ), \end{align}

(4.20) \begin{align} v^d_3(x)& = 4r^2(p''(x)\cos 2\theta - q''(x)\sin 2\theta ) + 4r^2 \bigl [\eta _2(x)(p'(x)\sin \theta + q'(x)\cos \theta )\bigr ]'. \end{align}

Replacing the unknown functions with their time expansions in the integral equations (3.5) and (3.12) for the case of a cylinder moving in the lower layer leads to the following system of recursive equations:

(4.21) \begin{align} \left[I + r^2 A^r_0 \right]v_n(x) - \left[H + r^2 B^r_0\right] u_n(x) &= \varphi _n(x), \\[-28pt] \nonumber \end{align}

(4.22) \begin{align} v_n(x) + H w_n(x) &= \psi _n(x) , \end{align}

where $I$ is the identity operator, and the non-homogeneous terms $\varphi _n$ and $\psi _n$ can be evaluated explicitly via the coefficients up to the $(n-1)$ th order by

(4.23) \begin{align} \varphi _n(x) &= v^d_n(x) + \sum \limits _{j=1}^n \Bigl (B_ju_{n-j}(x)-A_jv_{n-j}(x)\Bigr ) + r^2\sum \limits _{j=1}^n \Bigl (B^r_ju_{n-j}(x) - A^r_jv_{n-j}(x)\Bigr ), \\[-28pt] \nonumber \end{align}

(4.24) \begin{align} \psi _n(x) &= \sum \limits _{j=1}^n \Bigl (A_jv_{n-j}(x)-B_jw_{n-j}(x)\Bigr ). \\[0pt] \nonumber \end{align}

Equations (4.5), (4.6), (4.21) and (4.22) form the full set of recursive equations for the unknown coefficients $(\eta _n,u_n,w_n,v_n)$ . Using (4.6), the pair of integral equations (4.21) and (4.22) can be rewritten in operator form for a vector-valued unknown function $\boldsymbol v_n = (v_n,w_n)$ :

(4.25) \begin{equation} \Bigl (I + \frac {r^2}{1+\rho }\boldsymbol A_0^r\Bigr )\boldsymbol v_n(x) = \boldsymbol f_n(x), \end{equation}

where

(4.26) \begin{equation} \boldsymbol A_0^r = \begin{pmatrix} A_0^r & -\rho B_0^r \\ {\textit{HA}}_0^r & -\rho {\textit{HB}}_0^r \end{pmatrix}, \quad \boldsymbol f_n(x) = \frac 1{1+\rho }\begin{pmatrix} H\phi _n+ \varphi _n + \rho \psi _n + r^2 B_0^r\phi _n \\ -\phi _n + H\varphi _n - H\psi _n +r^2{\textit{HB}}_0^r\phi _n \end{pmatrix}. \end{equation}

Equation (4.25) is the Fredholm integral equation of the second kind, and is uniquely solvable for small values of the parameter $r$ and a given density parameter $0\leqslant \rho \lt 1$ (Trikomi Reference Trikomi1957), that is, in the case of a sufficiently large submergence depth of the cylinder. The solution is found iteratively using the convergent Neumann series

(4.27) \begin{equation} \boldsymbol v_n(x) = \sum \limits _{k = 0}^\infty \frac {(-1)^k r^{2k}}{(1+\rho )^k} \Bigl \{ \boldsymbol A_0^r\Bigr \}^k \boldsymbol f_n(x). \end{equation}

We now have a recursive procedure to construct the small-time asymptotic solution. By solving the system of equations (4.5), (4.6), (4.21) and (4.22) sequentially for different orders $n$ , we find the first four coefficients for the interface displacement:

(4.28) \begin{align} \eta _1(x) &= 0, \\[-12pt] \nonumber \end{align}

(4.29) \begin{align} \eta _2(x) &= (1+a)\left (r^2 - \frac {ar^4}4\right )(q'(x)\sin \theta - p'(x)\cos \theta ) + O(r^6), \end{align}

(4.30) \begin{align} \eta _3(x) &= 0, \end{align}

(4.31) \begin{align} \quad \eta _4(x) &=\frac {1+a}2 \left (r^2 - \frac {ar^4}4\right )\!\left [\left (\!\frac {\lambda a}6\cos \theta - \sin 2\theta \!\right )q''(x) + \left ( \cos 2\theta + \frac {\lambda a}6 \sin \theta \!\right )\!p''(x)\!\right ] \nonumber \\ &\quad + \frac {a(1+a)}8r^4\left [\left (\cos 2\theta + \frac {\lambda (1+a)}6\sin \theta \right )q'(x) \right . \nonumber\\ &\quad+ \left . \left ( \sin 2\theta - \frac {\lambda (1+a)}6 \cos \theta \right )p'(x)\right ] \nonumber \\ &\quad + \frac {1+a}8 r^4\left [ \left (\frac 13 + \frac 13a - \frac 29a^2\right )\Bigl ( p^{(4)}(x)\cos 2\theta - q^{(4)}(x)\sin 2\theta \Bigr )\right . \nonumber\\ &\quad+ \left . a\left (a+\frac 13\right )p''(x) - a\left (a+\frac 43\right ) q'(x)\right ] + O(r^6). \end{align}

Detailed calculation for the above-mentioned solution can be found in Appendices A and B. It is worth noting that for $a = 1$ , (4.28)–(4.31) describe the shape of the free surface, as predicted by Kostikov & Makarenko (Reference Kostikov and Makarenko2018) for a cylinder moving in a one-layer fluid. In the limit of a weakly stratified fluid, when the densities of the two fluids are almost equal ( $a \to 0$ ), the contribution of the nonlinear terms of order $r^4$ is reduced to zero, while the leading linear terms of order $r^2$ are reduced by half compared with the limiting case of free surface flow ( $a = 1$ ).

In the case of a cylinder moving in the upper layer of a two-layer fluid, the coefficients can be derived following the same principle:

(4.32) \begin{align} \eta _1(x) &= 0, \end{align}

(4.33) \begin{align} \eta _2(x) &= (1-a)\left (r^2 + \frac {ar^4}4\right )(q'(x)\sin \theta + p'(x)\cos \theta ), \end{align}

(4.34) \begin{align} \eta _3(x) &= 0, \end{align}

(4.35) \begin{align} \eta _4(x) &=-\frac {1-a}2 \left (r^2 + \frac {ar^4}4\right )\!\left [\!\left (\sin 2\theta + \frac {\lambda a}6\cos \theta \!\right )q''(x) + \left ( \cos 2\theta - \frac {\lambda a}6 \sin \theta \!\right )p''(x)\!\right ] \nonumber \\ &\quad + \frac {a(1-a)}8r^4\left [\left (\cos 2\theta + \frac {\lambda (1-a)}6\sin \theta \right )q'(x) \right . \nonumber\\ &\quad - \left . \left ( \sin 2\theta - \frac {\lambda (1-a)}6 \cos \theta \right )p'(x)\right ] \nonumber\\ &\quad + \frac {1-a}8 r^4\left [ \left (-\frac 13 + \frac 13a + \frac 2{9}a^2\right )\Bigl ( p^{(4)}(x)\cos 2\theta + q^{(4)}(x)\sin 2\theta \Bigr )\right . \nonumber\\ &\quad - \left . a\left (a-\frac 13\right )p''(x) + a\left (a-\frac 43\right ) q'(x)\right ]\!. \end{align}

It is evident that in the case of a one-layer fluid with a free surface ( $a = 1$ ), the coefficients (4.33) and (4.35) turn to zero and the interface disappears. In a weakly stratified fluid ( $a \to 0$ ), contribution of nonlinear terms with $r^4$ is reduced to zero, and the leading linear terms with $r^2$ start to dominate the flow with increasing magnitude.

Ideally, the equivalence of solutions (4.28)–(4.31) and (4.32)–(4.35) must obey a certain conversion rule, in accordance with the inversion of the direction of the gravity force relative to the interface. Comparing the above-mentioned coefficients, we can conclude that to the first order of accuracy, the solution for the cylinder in the upper layer (4.33) can be obtained from the solution for the cylinder in the lower layer (4.29) by inverting the signs of parameter $a$ and the differentiated kernel function $p'(x)$ . A different principle applies to the next order of accuracy. Similarly, coefficient (4.35) differs from coefficient (4.31) in the signs of parameter $a$ , and the differentiated kernel functions $p'(x)$ and $q''(x)$ , but not in a general way. This is due to a cancellation effect occurring as a result of multiple conversions at this order.

Equations (4.28)–(4.35) provide explicit analytical representations of the instantaneous shape of the interface of a two-layer fluid, when a circular cylinder starts to move with constant acceleration from rest in one of the layers. The flow regimes are determined by four parameters: the radius of the cylinder ( $r$ ), the direction of motion ( $\theta$ ), the rate of acceleration ( $\lambda$ ) and the relative density ratio (the Atwood number $a$ ). Certain limitations on the range of values for parameter $\lambda$ should be imposed to avoid a growing cavity behind the body. A situation when the fluid detaches from the body surface, forming the cavity, can occur when a cylinder is abruptly displaced from its initial position. The safe range for the Froude number, for which the fluid pressure does not decrease below a certain level around a circular cylinder in a one-layer fluid, was established by Norkin (Reference Norkin2012). In our notation, this bound corresponds to the condition $\lambda \gt 0.78$ . To the best of our knowledge, there is no known study of this kind of cavitation phenomenon for a two-layer fluid. In the subsequent analysis, we assume that the bound $\lambda \gt 0.78$ holds for flow without cavity formation.

The asymptotic solution (4.28)–(4.35) has been obtained under the assumption of small time and small radius, using infinite series with respect to the corresponding parameters. Hence, the validity of the results is limited to the analysis of small obstacles within a short time interval. Earlier, Pardo & Nedic (Reference Pardo and Nedic2021) extended our asymptotic solution for a one-layer fluid (Kostikov & Makarenko Reference Kostikov and Makarenko2018) to study cylinders of arbitrary radius. They used an iterative numerical algorithm to solve the boundary integral equation on the free surface to calculate coefficients up to the eighth order, valid without any assumptions on the size of the cylinder. By comparing with their numerical calculations, they established that the fourth-order analytical solution (4.28)–(4.31) with $a = 1$ works well for cylinders of radius $r \leqslant 0.5$ and for the time interval $t \lt 0.8$ . They investigated the temporal validity of the time expansions depending on the number of calculated terms and concluded that the the fourth-order solution deviates from the eighth-order solution by less than 5 % for $t \lt 0.6$ when $\lambda = 5$ and for $t \lt 0.85$ when $\lambda = 0.1$ . Later, Pardo et al. (Reference Pardo, Barua, Lisak and Nedic2022) validated the analytical model of Kostikov & Makarenko (Reference Kostikov and Makarenko2018) and semi-analytical model of Pardo & Nedic (Reference Pardo and Nedic2021) experimentally. Solution (4.28)–(4.35) generalises the results of Kostikov & Makarenko (Reference Kostikov and Makarenko2018) to a two-layer fluid. Based on analogy with the one-layer free surface problem, a similar range of applicability with respect to the radius $r$ and time $t$ may be anticipated. However, this validity range remains an inferred estimate, since it has not yet been independently demonstrated or validated for the present two-layer problem.

5. Hydrodynamic forces

The vector of hydrodynamic forces $(X,Y)$ acting on the cylinder, accelerating in the lower layer of a two-layer fluid, can be found from Sedov’s formula (Sedov Reference Sedov1965):

(5.1) \begin{equation} X(t) + i Y(t) = \pi r^2 \ddot z_c(t) + i \frac {\rm d}{{\rm d} t}\underset {S}{\int } z F_+(z,t) \,{\rm d} z + \overline {\frac {i}{2} \underset {S}{\int } (F_+(z,t))^2\, {\rm d}z}, \end{equation}

where $S$ is the contour of the cylinder and $F(z,t) = U_+ - i V_+$ is the complex velocity of the fluid in the lower layer $\varOmega _+$ . Since Sedov’s formula (5.1) and its subsequent evaluation are invariant with respect to the location of the cylinder relative to the interface, the function $F(z,t)$ will herafter be written without a subscript.

To calculate the contour integrals in (5.1), we use representation (3.4) for function $F(z,t)$ , introduced in § 3. For clarity of the subsequent discussion, we reproduce this formula here:

(5.2) \begin{equation} 2\pi i F(z,t) = \int \limits _{\varGamma (t)} \frac {F(\zeta ,t)\,{\rm d}\zeta }{\zeta - z} + \frac {r^2}{(z - z_c(t))^2}\overline {\int \limits _{\varGamma (t)}\frac {F(\zeta ,t)\,{\rm d}\zeta }{\zeta - z_{\ast }(t)}} + \frac {\gamma }{z - z_c(t)} + \frac {2\pi ir^2\dot z_c(t)}{(z-z_c(t))^2}. \end{equation}

By expanding the complex quantity $z_\ast (t) = z_c(t) + r^2/\overline {z - z_c(t)}$ in powers of $r$ in the last integral term of formula (5.2), it can be rewritten as a multipole expansion:

(5.3) \begin{equation} F(z,t) = F^{(0)}(z,t) + \frac 1{2\pi i}\frac {\gamma }{z - z_c(t)}+ \frac {r^2\dot z_c(t)}{(z-z_c(t))^2}- \frac {r^2}{(z-z_c(t))^2} \sum \limits _{n=0}^{\infty }\frac {r^{2n}\overline {F^{(n)}(z_c,t)}}{(z-z_c(t))^n}, \end{equation}

where

(5.4) \begin{equation} F^{(n)}(z,t) = \frac {1}{2\pi i} \underset {\varGamma (t)}{\int }\frac {F(\zeta ,t)\,{\rm d}\zeta }{(\zeta - z)^{n+1}} = \frac 1{n!}\frac {{\rm d}^n}{{\rm d}z^n} F^{(0)}(z,t), \quad n = 0,1, \ldots \end{equation}

For each $n = 0,1,\ldots$ , the function $F^{(n)}(z,t)$ is analytic in the fluid layer, occupied by the cylinder, that is, for $ z \in \{z \,|\, \textrm {Im}z\lt \eta (x,t)\}$ . From definition (5.4), it follows by the residue theorem:

(5.5) \begin{equation} \frac {1}{2 \pi i}\underset {S}{\int } \frac {F^{(0)}(\zeta ,t)\,{\rm d}\zeta }{(z-z_c(t))^{n+1}} = \textrm {Res}_{z = z_c(t)} \left (\frac {F^{(0)}(z,t)}{(z - z_c(t))^{n+1}}\right ) =F^{(n)}(z_c,t). \end{equation}

Substituting the multipole expansion (5.3) into Sedov’s formula (5.1) and using the residue theorem, we obtain the following formula for zero circulation ( $\gamma = 0$ ):

(5.6) \begin{align} \frac {X(t) + i Y(t)}{2\pi r^2} &= -\frac {\ddot z_c(t)}2 + \frac {{\rm d}\overline {F^{(0)}(z_c(t),t)}}{{\rm d}t} - \overline {\dot z_c(t)F^{(1)}(z_c(t),t)} \nonumber \\&\quad +\sum _{n=0}^\infty r^{2n}F^{(n)}(z_c(t),t),\overline {F^{(n+1)}(z_c(t),t)}. \end{align}

Golikov & Makarenko (Reference Golikov and Makarenko2022, Reference Golikov and Makarenko2024) introduced the representation (5.6) for the unsteady forces in powers of $r$ for both forced and free motions of the cylinder in a one-layer fluid. They observed that the first term in the total force in formula (5.6) is responsible for the added mass of the cylinder during its unsteady motion in an unbounded fluid. The second term corresponds to the integral of $zF(z,t)$ in (5.1) and reflects the unsteady nature of the process. The last term, which depends bilinearly on the coefficients $F^{(n)}$ , represents the integral of $(F(z,t))^2$ in (5.1) and is also present in the Blasius–Chaplygin formula for nonlinear loads in the case of stationary motion. A key property of the representation (5.6) is that it depends only on the instantaneous shape of the interface $\eta (x,t)$ , the velocity components $u(x,t)$ , $w(x,t)$ , $v(x,t)$ on it, and the trajectory of the cylinder $z_c(t)$ .

The function $F^{(n)}(z,t)$ , defined by the boundary integral (5.4), can be expanded in powers of $t$ , using the series for the velocity componenets $u_n(x)$ and $v_n(x)$ and the interface displacement $\eta _n(x)$ :

(5.7) \begin{align} 2 \pi i F^{(n)}(z,t) &= t \int \limits _{-\infty }^{+\infty } \frac {(i v_1(x) - u_1(x))\,{\rm d}x}{(x - z)^{n+1}} + t^3 \int \limits _{-\infty }^{+\infty } \frac {( iv_3(x) - u_3(x)) \,{\rm d} x}{(x - z)^{n+1}} \nonumber\\&\quad + (n+1) t^3 \int \limits _{-\infty }^{+\infty } \frac {(v_1(x) + iu_1(x)) \eta _2(x) \,{\rm d} x}{(x - z)^{n+2}} + O\big(t^5\big). \end{align}

By calculating the integrals in (5.7) up to third order in time $t$ and fourth order in radius $r$ , and substituting the result into (5.6), we obtain the following asymptotic formula for the hydrodynamic force on a cylinder moving in the lower layer of a two-layer fluid at an arbitrary angle to the interface:

(5.8) \begin{align} \frac {X(t) + i Y(t)}{2 \pi r^2} &= -\,e^{i \theta } - 2 r^2 a\left (1-\frac {a r^2}{4} \right )\left [ \frac {e^{ i \theta }}{(t^2 e^{-i \theta } +2i)^2} + t^2 \frac {6e^{ 2 i \theta } -i\lambda (1+a) e^{ i \theta }}{(t^2 e^{- i \theta } +2i)^3}\right ] \nonumber\\&\quad -\frac { ar^4}{2}\Bigg [ \frac {6ia e^{ 2 i \theta } + \lambda a(1+a) e^{ i \theta }-2i(3a^2+7a+1)}{2(t^2 e^{- i \theta } +\!2 i)^2} + \frac {2(3a^2+4a+1)}{(t^2 e^{- i \theta }+\!2 i)^3} \nonumber\\&\quad + \frac {16a}{(t^2 e^{- i \theta }+ \!2 i)^3 (t^2 e^{i \theta }-\! 2 i)^2}+\frac {16 (-a^2+2a+3) e^{2 i \theta }}{(t^2 e^{- i \theta } +\!2i)^5}\Bigg ]t^2 +\! O\big(t^4 + r^6\big). \end{align}

The detailed calculation of the integrals in (5.7) can be found in Appendices A and B.

Using the same methodology, the asymptotic formula for the hydrodynamic forces on a cylinder in the upper layer can be derived as

(5.9) \begin{align} \frac {X(t) + i Y(t)}{2 \pi r^2} &= -\,e^{i \theta } + 2 ar^2\left (1+\frac {ar^2}{4} \right )\left [ \frac {e^{ i \theta }}{(t^2 e^{-i \theta } -2i)^2} + t^2 \frac {6e^{ 2 i \theta } -i\lambda (1-a) e^{ i \theta }}{(t^2 e^{- i \theta } -2i)^3}\right ] \nonumber\\&\quad + \frac {ar^4}{2}\Bigg [ \frac {6 i ae^{ 2 i \theta } + \lambda a(1-a)e^{ i \theta }+2i(3a^2-7a+1)}{2(t^2 e^{- i \theta } -2 i)^2} + \frac {2(3a^2-4a+1)}{(t^2 e^{- i \theta }-2 i)^3} \nonumber\\&\quad - \frac {16a}{(t^2 e^{- i \theta }-2 i)^3 (t^2 e^{i \theta }+ 2 i)^2}+\frac {16 (a^2-2a+3) e^{-2 i \theta }}{(t^2 e^{- i \theta } -\!2i)^5}\Bigg ]t^2 +\! O\big(t^4 + r^6\big). \end{align}

Comparing the asymptotics (5.8) and (5.9) for the hydrodynamic forces, we conclude that (5.8) can be transformed into (5.9) almost completely by changing the sign of the parameter $a$ . The transformation of the higher-order terms in the expansion formulae is difficult to trace because of the multiple overlaps of various effects, such as gravity, the position of the cylinder relative to the interface and the direction of motion.

Neglecting the terms of the order $r^4$ in (5.8) and expanding the remaining terms in powers of $t$ , we obtain the leading-order asymptotic representation for the hydrodynamic force:

(5.10) \begin{equation} \frac {X(t) + iY(t)}{2\pi r^2} = -e^{i\theta } + \frac {ar^2}2e^{i\theta } - \frac {3ar^2}2\left (ie^{2i\theta } + \frac {(1+a)\lambda }6 e^{i\theta }\right )t^2 + O(t^4+r^4). \end{equation}

It is interesting to observe that the structure of the asymptotics (5.10) inherits the structure of the coefficients (4.29) and (4.31) in the time expansion of the interface deformations. The second and third terms in (5.10) represent the wave-induced forces on the cylinder and, in the limit of a homogeneous fluid ( $a \to 0$ ), vanish, making the force almost uniform. Analysis of (5.10) shows that the hydrodynamic force acting on the accelerated cylinder is generally reduced in a two-layer fluid ( $0 \lt a \lt 1$ ) compared with a homogeneous one-layer fluid ( $a = 1$ ). The specific influence of the density parameter $a$ on the forces for various directions of cylinder motion will be examined in detail in § 7.

6. Comparisons with alternative theories

The amount of research on unsteady interface flow in a two-layered fluid in the presence of an obstacle is limited and, to the best of our knowledge, there is no rigorous experimental or numerical validation of the asymptotic solution obtained here. Theoretical models for a one-layered fluid, which are more extensively represented in the literature, are mostly concerned with steady-state solutions to diffraction and radiation problems. Most existing models for the initial stage of wave motion investigate free surface deformations that occur when a cylinder starts to move impulsively with constant velocity (Haussling & Coleman Reference Haussling and Coleman1979) or are produced by constant body motion after an initial acceleration period (Telste Reference Telste1987; Gorlov Reference Gorlov2000). Laboratory experiments focus mainly on water-entry and water-exit problems (Greenhow Reference Greenhow1988; Moyo & Greenhow Reference Moyo and Greenhow2000; Vincent et al. Reference Vincent, Rivero-Rodriguez, Falla, Ashraf, Terrapon, Dorbolo and Scheid2026). It is therefore problematic to bring our theoretical model, specifically designed for initial perturbations of the interface, in strict correspondence with the limited available data, even for a one-layer fluid. Only recently, Pardo et al. (Reference Pardo, Barua, Lisak and Nedic2022) conducted the first experiments to estimate the prediction time of their small-time asymptotic solution obtained by a recursive semi-analytical method (Pardo & Nedic Reference Pardo and Nedic2021). Their method extended the asymptotic solution of Kostikov & Makarenko (Reference Kostikov and Makarenko2018) to cylinders of arbitrary radius under general submersion motions and thus provided solid experimental validation of the theoretical approach used in this paper. In this section, we provide the primary consistency check of our results for a two-layer fluid by comparing them with the exact solution of linear wave theory.

6.1. Linear water wave theory

In the theoretical solution of the linear problem of internal waves generated by a moving circular cylinder, the body can be modelled as a point dipole with a time-dependent moment, turning on at time $t = 0$ . The position of the dipole is determined by the trajectory of the cylinder $(x_c(t), y_c(t))$ and its intensity is given by $\mu (t) = - 2\pi r^2 U_c(t)$ , where $U_c(t)$ is the speed of the cylinder. Following the approach proposed by Bukreev et al. (Reference Bukreev, Gusev and Sturova1983, Reference Bukreev, Gusev and Sturova1986), the equations of motion for a cylinder below the interface in a two-layer fluid can be formulated in terms of the vertical fluid velocities in two layers:

(6.1) \begin{align}&\qquad\qquad\qquad \Delta V_- =0, \quad y \gt 0, \end{align}

(6.2) \begin{align}& \Delta V_+ = \mu (t)\delta (x - x_c(t))\frac {\partial }{\partial y}\delta (y - y_c(t)), \quad y \lt 0 \end{align}

(6.3) \begin{align}&\qquad\qquad\qquad V_- = V_+, \quad y = 0, \end{align}

(6.4) \begin{align}& \cfrac {\partial ^3 V_-}{\partial t^2 \partial y} - \lambda \cfrac {\partial ^2 V_-}{\partial x^2} = \frac 1\rho \left (\cfrac {\partial ^3 V_+}{\partial t^2 \partial y} - \lambda \cfrac {\partial ^2 V_+}{\partial x^2}\right ), \quad y = 0, \end{align}

(6.5) \begin{align}&\qquad\qquad\qquad V_{\pm } \to 0. \quad y \to \pm \infty . \end{align}

Here, $\delta$ is the Dirac delta function and the parameter $\lambda$ appears as a non-dimensional replacement for the gravitational acceleration $g$ . Once the system of equations (6.1)–(6.5) has been solved, the interface deformation can be found from the formula:

(6.6) \begin{equation} \frac {\partial \eta }{\partial t} = V_-, \quad y = 0. \end{equation}

To solve the initial-boundary value problem (6.1)–(6.5), we apply the Fourier transform with respect to the spatial variable $x$ and the Laplace transform with respect to the time variable $t$ :

(6.7) \begin{equation} \hat V(k,y,s) = \int \limits _0^\infty e^{-st}\, {\rm d}t \int \limits _{-\infty }^{\infty } V(x,y,t) e^{-ikx}\,{\rm d}x, \end{equation}

where the variable $k$ is real and the real part of the complex variable $s$ is positive. The transformed problem takes the form of a second-order boundary value problem for $\hat V(k,y,s)$ :

(6.8) \begin{align}&\qquad\qquad\qquad\qquad\qquad \hat V^{\prime\prime}_- - k^2 \hat V_- = 0, \quad y \gt 0 \end{align}

(6.9) \begin{align}& \hat V^{\prime\prime}_+ - k^2 \hat V_+ = \int _0^\infty \mu (t) \exp (-st - ikx_c(t))\frac {\partial }{\partial y}\delta (y - y_c(t))\,{\rm d}t, \quad y \lt 0 \end{align}

(6.10) \begin{align}&\qquad\qquad\qquad\qquad\qquad\quad \hat V_- = \hat V_+, \quad y = 0 \end{align}

(6.11) \begin{align}&\qquad\qquad\quad s^2 \hat V^{\prime}_- + \lambda k^2 \hat V_- = \frac 1\rho \left(s^2 \hat V^{\prime}_+ + \lambda k^2 \hat V_+\right), \quad y = 0 \end{align}

(6.12) \begin{align}&\qquad\qquad\qquad\qquad\qquad\,\, \hat V_{\pm } \to 0, \quad y \to \pm \infty . \end{align}

Here, the prime denotes differentiation with respect to $y$ . Once the problem (6.8)–(6.12) is solved, the transformed displacement of the interface from the unperturbed flat position can be found from the formula

(6.13) \begin{equation} \hat \eta (k,s) = \frac {\hat V_-(k,0,s)}s = \frac {(1+a)s}{s^2 + a\lambda k}\int \limits _0^\infty \mu (t) \exp {(-st - ikx_c(t) + k y_c(t))}\, {\rm d}t. \end{equation}

The solution of the original problem is found by taking the inverse Fourier and Laplace transforms of (6.13)

(6.14) \begin{equation} \eta (x,t) = \frac 1{2\pi }\int \limits _{-\infty }^\infty {\rm d}k \int \limits _0^\infty \mu (\tau ) \exp (ik(x-x_c(\tau ))-ky_c(\tau ))H(k,t,\tau )\,{\rm d}\tau , \end{equation}

where

(6.15) \begin{equation} H(k,t,\tau ) = \frac 1{2\pi i}\int \limits _{c-i\infty }^{c+i\infty } \frac {(1+a)se^{s(t-\tau )}}{s^2 + a\lambda k}\,{\rm d}s = a\cos \Bigl ( \sqrt {a\lambda k}(t - \tau )\Bigr ) .\end{equation}

The integral in (6.15) converges if the real part of $s$ exceeds some minimum value, set as $c_0$ , such that all singularities of the integrand $s = \pm i a\lambda k$ lie to the left of the line $\textrm {Re}(s) = c_0$ . The inverse Laplace transform given in (6.15) then exists for all $c \gt c_0$ (see Andrews & Shivamoggi Reference Andrews and Shivamoggi1998). Substituting the expression (6.15) into the integral formula (6.14) and taking into account that function $H(k,t,\tau )$ is non-zero only for $\tau \lt t$ , we obtain the final expression for the interface displacement of the linear theory:

(6.16) \begin{align} \eta (x,t) &= 2r^2(1+a)\int _0^\infty k\,{\rm d}k \boldsymbol{\cdot }\int _0^t e^{-ky_c(\tau )}\nonumber\\&\quad\times \Bigl [\dot x_c(\tau ) \sin k(x - x_c(\tau )) - \dot y_c(\tau )\cos k(x - x_c(\tau ))\Bigr ] \cos \sqrt {a\lambda k }(t - \tau )\,{\rm d}\tau . \end{align}

Following the same principle, we can derive the interface displacement for a cylinder moving in the upper layer, see Bukreev et al. (Reference Bukreev, Gusev and Sturova1983). In this case, the resulting interface displacement differs from (6.16) only by the multiplication constant $1-a$ .

Stepanyants & Sturova (Reference Stepanyants and Sturova2021) used the same principle of replacing a submerged body with a dipole to describe the elastic waves on the floating elastic ice sheet caused by the motion of a circular cylinder. They showed that small-amplitude waves in the linear approximation are not too sensitive to the shape of a moving body, and the correction to the circular shape of the cylinder by an ice sheet is negligible provided that the cylinder moves relatively slow or sufficiently far from the interface. Similarly, the dipole approximation can be justified for the motion of a cylinder below the free surface or in a two-layered fluid. In fact, this approach has been used extensively (Bukreev et al. Reference Bukreev, Gusev and Sturova1983, Reference Bukreev, Gusev and Sturova1986) to describe wave generation by the unsteady motion of a circular cylinder in a two-layer fluid. Earlier, Havelock (Reference Havelock1927) and Sretensky (Reference Sretensky1937) applied the method of images in the linear problem of water waves to describe the free surface flow caused by the horizontal motion of a cylindrical body, and arrived at the same integral formula (6.16) for the case $a = 1$ .

A comparison of the derived asymptotic formula for the interface displacement (4.28)–(4.31) with its analogue from the linear wave theory (6.16) reveals that, in both cases, the parameter $(1+a)$ appears as a multiplication factor for the entire solution. Furthermore, the parameter $a$ specifically affects the term with $\lambda$ , which is related to the Froude number.

Figure 2.

Snapshots of the interface for a cylinder of radius $r = 0.5$ , moving in different directions in the lower layer of a two-layered fluid with $\lambda = 3.5$ , predicted by higher-order asymptotic solution (4.29)–(4.31) (solid line) and by solution of the linear theory (6.16) (dashed line).

Figure 2 compares the higher-order asymptotic solution given by (4.28)–(4.31) with the solution from linear wave theory (6.16) for a cylinder moving in the lower layer of a two-layer fluid. The figure shows the interface displacement at regular time intervals for three directions of cylinder motion: horizontal ( $\theta = 0$ ), vertical submersion ( $\theta = -\pi /2$ ), vertical emersion ( $\theta = \pi /2$ ). Two distinct cases are displayed: a small density difference ( $a = 0.21$ ) and free surface flow ( $a = 1$ ), to highlight features specific to a two-layered fluid. Figure 2 demonstrates that at early times ( $t \lt 0.5$ ), the agreement between the two solutions is reasonably good. A further increase in time leads to growing discrepancies, which may be attributed to the increasing influence of nonlinear effects or the limited validity range of the asymptotic solution. As noted previously, agreement persists for a longer duration in a two-layered fluid with a small density difference, which is consistent with the structure of the asymptotic solution. It is also noteworthy that good agreement is maintained until later times for the case of horizontal motion near the wave crest pushed ahead by the cylinder.

6.2. The Boussinesq approximation

In the Boussinesq approximation, commonly used to describe the internal waves in a weakly stratified fluid, the densimetric parameter $\lambda ^\ast = (1-\rho )\lambda$ is introduced, where the non-dimensional parameter $\rho = \rho _-/\rho _+$ is the density ratio. Taking the limit $\rho \to 1$ , we can reproduce the boundary differential equations (3.2) and (3.3) for a weakly stratified two-layer fluid in the form:

(6.17) \begin{equation} \eta _t = v, \quad u_t-w_t + \frac 12\cfrac {\partial }{\partial x} \left (\cfrac {u^2 - 2\eta _x v(u - w) - w^2}{1+\eta _x^2}\right ) + \lambda ^\ast \eta _x = 0. \end{equation}

The integral equations (3.5) and (3.12), which complement the system (6.17), remain unchanged, as they contain no information about the density difference. Following the theoretical approach, discussed earlier, we obtain the small-time asymptotic solution for the interface displacement in the Bouassinesq approximation. The coefficients in the time expansion of this solution are given by

(6.18) \begin{align} \eta _1(x) &= 0, \quad \eta _2(x) = r^2(q'(x)\sin \theta - p'(x)\cos \theta ), \quad \eta _3(x) = 0, \nonumber \\ \eta _4(x) &=\frac {r^2}2\left [\left (\frac {\lambda ^\ast }{12}\cos \theta - \sin 2\theta \right )q''(x) + \left ( \cos 2\theta + \frac {\lambda ^\ast }{12} \sin \theta \right )p''(x)\right ] \nonumber \\ &\quad + \frac {r^4}{24}\left [\frac {\lambda ^\ast }4\left (q'(x)\sin \theta - p'(x)\cos \theta \right )+ \Bigl ( p^{(4)}(x)\cos 2\theta - q^{(4)}(x)\sin 2\theta \Bigr )\right ]. \end{align}

The vector of hydrodynamic forces (5.8) in the Boussinesq approximation reduces to

(6.19) \begin{equation} \frac {X + iY}{2\pi r^2} = - e^{i\theta } + it^2 r^2 \frac {\lambda _\ast e^{i\theta }}{(t^2e^{-i\theta }+2i)^3}. \end{equation}

Figure 3.

Interface $y = \eta (x,t)$ of a two-layer fluid with different density difference for a circular cylinder of radius $r = 0.4$ and parameter $\lambda = 10$ , moving vertically upwards in the lower layer ( $x_c = 0$ , $y_c = -1 + t^2$ ).

Figure 3 illustrates the good agreement between the two high-order asymptotic solutions, obtained with and without Boussinesq approximation, when the density difference is sufficiently small. As the density stratification intensifies, as shown in figure 3(b), the asymptotic solution based on the Boussinesq approximation becomes ineffective, as it substantially overestimates the amplitude of the interface flow.

7. Analysis of the asymptotic solution

In this section, we study the factor of a two-layer fluid in the initial stage of the wave–structure interaction. We investigate the effect of the density stratification parameter, namely Atwood number $a$ , on the dynamics of the interface flow of a two-layer fluid and the hydrodynamic forces exerted on the moving body. It is interesting to compare the initial perturbations formed on the interface of a two-layer fluid as a result of the unsteady motion of a cylinder in the lower and upper layers.

7.1. Motion away from the interface

Figures 4(a, b) and 4(c, d) show the evolution of the interface with time when the cylinder starts to move vertically away from the interface in the upper and lower layer, respectively. Two distinct configurations of a two-layer fluid with intermediate ( $a = 0.5$ ) and strong ( $a = 0.95$ ) stratification are displayed. Figure 4(a, b) shows that the initially planar interface is entrained by the elevating cylinder in the upper layer and forms a central elevation zone, the magnitude of which is governed by the density difference. When the density difference is large, the interface is almost insensitive to the action of a moving body in the upper layer. Similarly, a central depression zone forms above a submerging cylinder in the lower layer, but the morphology of the wave is quite different. According to figure 4(c, d), at later times, the interface flattens and evolves into a vertical jet, which occurs as the result of the collision of small-scale waves propagating towards each other along the moving curvilinear interface. When the density difference is small, the jet forms at later times and cannot be fully observed within the validity period of the asymptotic solution. The variation in the wave response observed in figure 4 for the two different stratification cases is due to the fact that the interface of a two-layer fluid experiences different excessive pressures from the top layer, which alleviates the wave motion generated by a moving obstacle.

Figure 4.

Snapshots of the interface $y = \eta (x,t)$ of a two-layer fluid at three successive moments of time for the cylinder of radius $r = 0.4$ and parameter $\lambda = 4$ , moving (a, b) vertically upwards in the upper layer ( $x_c = 0$ , $y_c = 1 + t^2$ ) and (c, d) vertically downwards in the lower layer ( $x_c = 0$ , $y_c = -1 - t^2$ ).

Figures 5(a) and 5(b) compare the hydrodynamic forces predicted by (5.8) and (5.9), respectively, acting on the cylinder in the situations depicted in figure 4. The vertical foce acting on the lower cylinder grows with time, reaching an extremum at the point where the second-order and fourth-order terms come into balance and the jet forms at the interface. The rate of growth is larger for stronger stratification. The force on the upper cylinder experiences a similar growing trend, but without reaching an extremum.

Figure 5.

Vertical force on a cylinder of radius $r = 0.4$ and $\lambda = 4$ , moving vertically (a) downwards in the lower layer ( $x_c = 0, y_c = -1 - t^2$ ) and (b) upwards in the upper layer ( $x_c = 0, y_c = 1 + t^2$ ) of a two-layered fluid with different density difference.

Figure 6.

Snapshots of the interface $y = \eta (x,t)$ of a two-layer fluid ( $a = 0.4$ ) generated by vertical submersion of the cylinder in the lower layer ( $x_c = 0, y_c = -1-t^2$ ) of different radii $r$ and parameter $\lambda$ at time $t = 0.7$ .

Figure 6 shows representative instantaneous profiles of the interface as affected by variation in input parameters, such as the square of the inverse Froude number $\lambda$ and the non-dimensional cylinder radius $r$ . According to figure 6, the parameter $\lambda$ , which is related to the acceleration rate of the body and its submergence depth, determines the characteristic shape of the interface. Over time, the restoring force of gravity competes with the intertia of the flow, resulting in flattening of the initial elevation of the interface. An increase in cylinder size leads to a larger affected area of the interface and to the splitting of the splash-jet.

Jetting at the interface may be considered an undesirable event, as it can lead to instability of the interface flow. For example, in applications such as magnetised target fusion technology (Pardo et al. Reference Pardo, Barua, Lisak and Nedic2022), the interface between the magnetised plazma and the cavity must be maintained as smooth as possible. Therefore, it is beneficial to know under which conditions the jetting regime occurs. Figure 7 shows the perturbation height $\eta _0 = \eta (0,t)$ at the point on the interface above the centre of the cylinder $x = 0$ for various values of the density difference parameter ( $a$ ), acceleration rate ( $\lambda$ ) and cylinder size ( $r$ ). The plots are characterised by a smooth descent followed by a reversal, where the interface stops sinking after the submerging body and changes direction. As shown in figure 7, the presence of a second layer ( $a \lt 1$ ), as well as an incresed cylinder velocity, results in a delayed and suppressed splash response of the interface. An increase in cylinder radius leads, as expected, to deeper deflection but has no effect on the turnaround point.

Investigating the extrema of the interface displacement as a function of time, $y_0 = \eta (0,t)$ , yields an explicit representation for the critical time, at which the depression of the interface evolves into a reversed motion, forming a vertical jet:

(7.1) \begin{equation} t_c = \sqrt {\frac {1-ar^2/4}{2+\lambda a/3-r^2(2-13a/4+\lambda a/24+\lambda a^2/8-7a^2/12)}}. \end{equation}

Neglecting terms with $r^2$ in (7.1), we obtain an approximation for the critical time valid for small-sized cylinders:

(7.2) \begin{equation} t_c = \sqrt {\frac 1{2+\lambda a/3}}. \end{equation}

Equation (7.2) demonstrates that, to leading order, the critical time is inversely proportional to the product of parameters $\lambda$ and $a$ , meaning that the change of the characteristic speed of the cylinder can be compensated by proportional change in the density difference.

Figure 7.

Vertical submersion of a circular cylinder of radius $r = 0.4$ in the lower layer of a two-layer fluid ( $x_c = 0, y_c = -1 - t^2$ ): (a) displacement of the interface at zero point $\eta _0 = \eta (0,t)$ for $\lambda = 4$ ; (b) contour plots of critical time $t_c$ as a function of parameters $a$ and $\lambda$ .

7.2. Motion to the interface

Figures 8(a, b) and 8(c, d) show the initial dynamics of the interface when the cylinder starts to move vertically towards the interface in the upper and lower layers, respectively. Two distinct configurations of a two-layer fluid with intermediate ( $a = 0.5$ ) and strong ( $a = 0.95$ ) stratification are displayed. As seen in figure 8(a, b), the downward motion of the cylinder in the upper layer induces a vertical deflection of the interface without subsequent relaxation. In figure 8(c, d), the vertical uplift of an added fluid mass by the upward motion of the cylinder in the lower layer can be observed. The disturbances induced on the interface by the upper and lower cylinders are not vertical reflections of each other. The inertial fluid layer formed on top of the cylinder can split into small-scale waves flowing around the sides of the lifting body, which does not occur when the cylinder approaches the interface from above. The overall magnitude of the disturbances produced on the interface is larger for a smaller density difference.

Figure 8.

Snapshots of the interface $y = \eta (x,t)$ of a two-layer fluid at three successive moments of time for the cylinder of radius $r = 0.4$ and parameter $\lambda = 10$ , moving (a, b) vertically downwards in the upper layer ( $x_c = 0$ , $y_c = 1 - t^2$ ); (c, d) vertically upwards in the lower layer ( $x_c = 0$ , $y_c = -1 + t^2$ ).

Figure 9.

Vertical force on a cylinder of radius $r = 0.4$ and $\lambda = 10$ , moving vertically (a) upwards in the lower layer ( $x_c = 0, y_c = -1 + t^2$ ) and (b) downwards in the upper layer ( $x_c = 0, y_c = -1 - t^2$ ) of a two-layered fluid with different density difference.

Figures 9(a) and 9(b) show the hydrodynamic forces predicted by (5.8) and (5.9), respectively, on the cylinder moving in the upper and lower layer of a two-layer fluid. According to figure 9, the force acting on the cylinder approaching the interface from both the upper and lower layers exhibits a similar growing trend, the rate of which is highly dependent on the density stratification. When the density difference is sufficiently small, the hydrodynamic force on the cylinder moving in the lower fluid remains almost constant during the initial phase of the wave flow.

Figure 10 shows instantaneous profiles of the interface for variable input parameters $\lambda$ and $r$ when the cylinder is lifted up towards the interface in the lower layer. The general form for all conditions is a swelling of the interface as the cylinder approaches it. The volume of the swell and its surge speed depend on the radius of the cylinder $r$ . According to figure 10, for a comparatively fast motion of the cylinder ( $\lambda \lt \lambda _c$ ), inertial forces dominate and the interface always rises. For slow motion of the cylinder ( $\lambda \gt \lambda _c$ ), the contribution of the gravity-induced outflow becomes critical and the elevation reaches a local maximum. The critical value of $\lambda _c$ is determined by the density stratification through the parameter $a$ .

Figure 10.

Snapshots of the interface $y = \eta (x,t)$ of a two-layer fluid ( $a = 0.4$ ) generated by vertical elevation of the cylinder in the lower layer ( $x_c = 0, y_c = -1+t^2$ ) of different radius $r$ and parameter $\lambda$ at time $t = 0.7$ .