Numerical modelling of surface water wave interaction with a moving wall

TL;DR

This paper presents a numerical model for water wave interaction with a horizontally movable wall attached to springs, enabling wave energy damping through a semi-analytical finite difference method.

Contribution

It introduces a novel semi-analytical numerical approach for simulating wave interaction with a spring-mounted moving wall, allowing energy extraction from incident waves.

Findings

The method accurately models wave-wall interactions.

Spring parameters can be tuned for wave energy damping.

Validation confirms the model's reliability.

Abstract

In the present manuscript, we consider the practical problem of wave interaction with a vertical wall. However, the novelty here consists in the fact that the wall can move horizontally due to a system of springs. The water wave evolution is described with the free surface potential flow model. Then, a semi-analytical numerical method is presented. It is based on a mapping technique and a finite difference scheme in the transformed domain. The idea is to pose the equations on a fixed domain. This method is thoroughly tested and validated in our study. By choosing specific values of spring parameters, this system can be used to damp (or in other words to extract the energy of) incident water waves.

Tables

Table 1. Table 1. Matrix elements of the difference equation ( 5.7 ).

Node type	${𝒜}_{\mathrm{\hspace{0.17em}1},𝒋}$	${𝒜}_{\mathrm{\hspace{0.17em}2},𝒋}$	${𝒜}_{\mathrm{\hspace{0.17em}3},𝒋}$	${𝒜}_{\mathrm{\hspace{0.17em}4},𝒋}$	${𝒜}_{\mathrm{\hspace{0.17em}5},𝒋}$	${𝒜}_{\mathrm{\hspace{0.17em}6},𝒋}$	${𝒜}_{\mathrm{\hspace{0.17em}7},𝒋}$	${𝒜}_{\mathrm{\hspace{0.17em}8},𝒋}$
${𝒒}_{𝒋}\in {𝒬}_{\mathrm{\hspace{0.17em}0}}^h$	${\text{<span>textrho</span>}}_{\mathrm{\hspace{0.17em}1}}$	${\text{<span>textrho</span>}}_{\mathrm{\hspace{0.17em}2}}$	${\text{<span>textrho</span>}}_{\mathrm{\hspace{0.17em}3}}$	${\text{<span>textrho</span>}}_{\mathrm{\hspace{0.17em}4}}$	${\text{<span>textrho</span>}}_{\mathrm{\hspace{0.17em}5}}$	${\text{<span>textrho</span>}}_{\mathrm{\hspace{0.17em}6}}$	${\text{<span>textrho</span>}}_{\mathrm{\hspace{0.17em}7}}$	${\text{<span>textrho</span>}}_{\mathrm{\hspace{0.17em}8}}$
${𝒒}_{𝒋}\in {\gamma }_{\ell }^h$	$0$	$\frac{1}{2}{\rho }_{\mathrm{\hspace{0.17em}2}}-{\text{<span>textsigma</span>}}_{\mathrm{\hspace{0.17em}2}}$	${\rho }_{\mathrm{\hspace{0.17em}3}}$	$\frac{1}{2}{\rho }_{\mathrm{\hspace{0.17em}4}}-{\text{<span>textsigma</span>}}_{\mathrm{\hspace{0.17em}4}}$	$0$	${\rho }_{\mathrm{\hspace{0.17em}6}}$	${\rho }_{\mathrm{\hspace{0.17em}7}}$	$0$
${𝒒}_{𝒋}\in {\gamma }_b^h$	$\frac{1}{2}{\rho }_{\mathrm{\hspace{0.17em}1}}-{\text{<span>textsigma</span>}}_{\mathrm{\hspace{0.17em}1}}$	$0$	$\frac{1}{2}{\rho }_{\mathrm{\hspace{0.17em}3}}-{\text{<span>textsigma</span>}}_{\mathrm{\hspace{0.17em}3}}$	${\rho }_{\mathrm{\hspace{0.17em}4}}$	$0$	$0$	${\rho }_{\mathrm{\hspace{0.17em}7}}$	${\rho }_{\mathrm{\hspace{0.17em}8}}$
${𝒒}_{𝒋}\in {\gamma }_r^h$	${\rho }_{\mathrm{\hspace{0.17em}1}}$	$\frac{1}{2}{\rho }_{\mathrm{\hspace{0.17em}2}}+{\text{<span>textsigma</span>}}_{\mathrm{\hspace{0.17em}2}}$	$0$	$\frac{1}{2}{\rho }_{\mathrm{\hspace{0.17em}4}}+{\text{<span>textsigma</span>}}_{\mathrm{\hspace{0.17em}4}}$	${\rho }_{\mathrm{\hspace{0.17em}5}}$	$0$	$0$	${\rho }_{\mathrm{\hspace{0.17em}8}}$
${𝒒}_{\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}0}}$	$0$	$0$	$\frac{1}{2}{\rho }_{\mathrm{\hspace{0.17em}3}}-{\text{<span>textsigma</span>}}_{\mathrm{\hspace{0.17em}3}}$	$\frac{1}{2}{\rho }_{\mathrm{\hspace{0.17em}4}}-{\text{<span>textsigma</span>}}_{\mathrm{\hspace{0.17em}4}}$	$0$	$0$	${\rho }_{\mathrm{\hspace{0.17em}7}}$	$0$
${𝒒}_{N_1,\mathrm{\hspace{0.17em}0}}$	$\frac{1}{2}{\rho }_{\mathrm{\hspace{0.17em}1}}-{\text{<span>textsigma</span>}}_{\mathrm{\hspace{0.17em}1}}$	$0$	$0$	$\frac{1}{2}{\rho }_{\mathrm{\hspace{0.17em}4}}+{\text{<span>textsigma</span>}}_{\mathrm{\hspace{0.17em}4}}$	$0$	$0$	$0$	${\rho }_{\mathrm{\hspace{0.17em}8}}$

Equations

$$x=\xi (q,t),y=\eta (q,t),$$

$$\displaystyle\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}\$$

$$\displaystyle\boldsymbol{u}_{\,t}\ +\ (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})\,\boldsymbol{u}\ +\ \frac{1}{\rho}\;\boldsymbol{\nabla}p\$$

$$\mathbf{u}\cdot \mathbf{n}=0,\mathbf{x}\in {\Gamma }_b\cup {\Gamma }_r,t⩾0,$$

$$\mathbf{u}\cdot \mathbf{n}=\dot{s}(t),\mathbf{x}\in {\Gamma }_{\ell },t⩾0.$$

$$u=\dot{s}(t),\mathbf{x}\in {\Gamma }_{\ell },t⩾0.$$

$${\eta }_t+u\cdot {\eta }_x=v,\mathbf{x}\in {\Gamma }_s,t⩾0$$

$$p=p_a=\mathrm{const},\mathbf{x}\in {\Gamma }_s,t⩾0$$

$$\frac{\mathrm{d}\/\boldsymbol{x}}{\mathrm{d}\/t}\ =\ \boldsymbol{u}\,\bigl{(}\boldsymbol{x}(t),\,t\bigr{)}\,.$$

$$\displaystyle\frac{\mathrm{d}\/x}{\mathrm{d}\/t}\$$

$$\displaystyle\frac{\mathrm{d}\/y}{\mathrm{d}\/t}\$$

$$\displaystyle\frac{\partial\/\xi}{\partial\/t}\ +\ \frac{\partial\/\xi}{\partial\/q}\cdot\frac{\mathrm{d}\/q}{\mathrm{d}\/t}\ -\ u\$$

$$\displaystyle\frac{\partial\/\eta}{\partial\/t}\ +\ \frac{\partial\/\eta}{\partial\/q}\cdot\frac{\mathrm{d}\/q}{\mathrm{d}\/t}\ -\ v\$$

$$\mathrm{rot}\mathbf{u}=-\frac{\partial \text{/}u}{\partial \text{/}y}+\frac{\partial \text{/}v}{\partial \text{/}x}=0.$$

$$\mathbf{u}=\nabla \phi .$$

$$\Delta \phi =0,(\mathbf{x},t)\in \Omega (t)\times \text{mathds}{R}^{+}.$$

$${\eta }_t+{\phi }_x\cdot {\eta }_x={\phi }_y,y=\eta (x,t).$$

$${\phi }_t+\frac{1}{2}|\nabla \phi {|}^2+g\eta =0,y=\eta (x,t).$$

$$\nabla \phi \cdot \mathbf{n}=0,\mathbf{x}\in {\Gamma }_b\cup {\Gamma }_r,$$

$$\nabla \phi \cdot \mathbf{n}\equiv \frac{\partial \text{/}\phi }{\partial \text{/}x}=\dot{s}(t),\mathbf{x}\in {\Gamma }_{\ell }.$$

$$\xi (q,0)={\xi }_0(q),\eta (q,0)={\eta }_0(q),q\in [0,1].$$

$$\mathbf{u}(\mathbf{x},0)={\mathbf{u}}_0(\mathbf{x}),\mathbf{x}\in \Omega (0).$$

$$\phi (\mathbf{x},0)={\phi }_0(\mathbf{x})\equiv {\phi }_0({\mathbf{x}}_0)+{\int }_{\text{mbox}{\text{textgamma}}_{{\mathbf{x}}_0⇝\mathbf{x}}}{u}_0\mathrm{d}x+v_0\mathrm{d}y,$$

$$m\,\ddot{s}\ +\ k\,s\ =\ -\,\bigl{[}\,\mathscr{F}\,(t)\ -\ \mathscr{F}\,(0)\,\bigr{]}\,,$$

$$\mathscr{F}\,(t)\ =\ \int_{\,-h_{0}}^{\,\eta_{s}\,(t)}\,p\,\bigl{(}s(t),\,y,\,t\bigr{)}\;\mathrm{d}y\,,\qquad\eta_{\,s}\,(t)\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \eta\,\bigl{(}s\,(t),\,t\bigr{)}\,.$$

$$\mathcal{F}(0)=g\frac{h_0^2}{2}.$$

$$p(x,y,t)=-{\phi }_t-\frac{1}{2}|\nabla \phi {|}^2-gy.$$

$$s(0)=\dot{s}(0)=0.$$

$$x^{\ast }=\frac{x}{h_0},y^{\ast }=\frac{y}{h_0},t^{\ast }=t\sqrt{\frac{g}{h_0}}$$

$$s^{\ast }=\frac{s}{h_0},{\eta }^{\ast }=\frac{\eta }{h_0},\hat{A}{\phi }^{\ast }=\frac{\phi }{h_0\sqrt{g{h}_0}},p^{\ast }=\frac{p}{\rho g{h}_0},{\mathcal{F}}^{\ast }=\frac{\mathcal{F}}{\rho g{h}_0^3}$$

$${\ell }^{\ast }=\frac{\ell }{h_0},m^{\ast }=\frac{m}{\rho h_0^3},k^{\ast }=\frac{k}{\rho g{h}_0^2}.$$

$$\displaystyle\boldsymbol{\nabla}^{2}\,\varphi\$$

Full text

Key words and phrases: wave/body interaction; full Euler equations; movable wall; wave run-up; wave damping; wave energy extraction

MSC:

PACS:

∗ Corresponding author

Last modified:

1. 1 Introduction

2. 2 Problem formulation

3. 2.1 Potential flow model

4. 2.1.1 Initial conditions

5. 2.2 Piston motion

6. 2.3 Dimensionless problem

7. 3 Equations on a fixed domain

8. 3.1 On some properties of elliptic operators in curvilinear coordinates

9. 3.1.1 Example 1

10. 3.1.2 Example 2

11. 4 Finite difference scheme in curvilinear coordinates

12. 4.1 Approximation in interior nodes

13. 4.2 Approximation of boundary conditions

14. 4.2.1 Consistency

15. 4.3 Difference operator

16. 4.4 Mathematical study of the finite difference problem for the velocity potential

17. 4.4.1 Some discrete identities

18. 4.5 The main result

19. 4.5.1 Final remarks

20. 5 Numerical algorithm

21. 5.1 Approximation of the dynamic boundary condition (1)

22. 5.2 Computation of the velocity potential (2)

23. 5.3 Free surface motion (3)

24. 5.4 Moving wall displacement (4)

25. 5.5 Elliptic mapping construction (5)

26. 5.5.1 The re-computation

27. 5.6 Stability of the scheme

28. 5.6.1 Practical choice of the time step

29. 6 Numerical results

30. 6.1 Solitary wave run-up on a fixed wall

31. 6.2 Wave generation by a numerical wave maker

32. 6.3 Surface wave run-up on a movable wall

33. 6.3.1 Single pulse

34. 6.3.2 Solitary wave

35. 7 Discussion

36. 7.1 Conclusions

37. 7.2 Perspectives

38. 7.3 Side effects

39. Acknowledgments

40. A Transformation of coordinates

41. B Transformation of the Laplace operator

42. C Boundary conditions transformation

43. D Piston motion modelling

Introduction

The mathematical and numerical high fidelity modeling of water waves is a central topic in coastal and naval engineering. The incident waves come and interact with various coastal features. Nowadays the interaction of water waves with bathymetric features is relatively well understood (at least with some special features [49]). A more challenging situation is the interaction of waves with various coastal structures [53]. At this level the most studied situation is the wave/wall interaction and the wall is assumed to be fixed and impermeable. Violent wave impacts have to be modelled in general using more CFD-like methods in the Eulerian mesh-based [20, 9] or Lagrangian particle-based [21] approaches. In the present study we apply the free surface approximation by neglecting all the processes happening in the air above the free surface [10]. In this line of thinking, the interaction of periodic waves with a fixed vertical wall was recently studied in [5] in the framework of a fully nonlinear weakly dispersive wave model [24]. The conditions under which an extreme wave run-up on a vertical wall may happen were describe in [5] as well.

In this study, we focus on wave interactions with a movable vertical wall. The wall motion can be prescribed. In this case we model the wave generation process with a piston-type wave maker [38]. For instance, this problem was considered in the framework of Boussinesq-type equations in [63]. Otherwise, the wall can move under the action of incident waves. We can even assume that a system of horizontal tension/extension springs (with tunable rigidities) is attached to the wall. Thus, the wall may present a certain resistance to the action of waves. This problem can be regarded also as a piston’s free motion under the forcing of water waves. The piston mechanical energy conversion and recuperation is a different technological problem which is out of scope of the present study. The extraction of ocean wave energy on industrial scales is not yet a very common practice [33]. However, very active researches in this field are conducted since at least forty years [31, 32]. Consequently, the numerical methods developed below can be used to design and to optimize such Wave Energy Conversion (WEC) devices. Movable walls have been used, for example, in a triplate system proposed back in 1977 by Dr. Francis Farley. The rigidity of strings is chosen to minimize the reflected wave amplitude so that most of the wave energy be converted into the mechanical energy of the device. Mathematical and numerical modeling of this type of WEC devices is considered below. The main point is that (ordinary and partial) differential equations posed in time-varying domains are known to pose notorious theoretical and numerical difficulties [47].

In the present work we consider a two-dimensional non-stationary problem of surface wave motion in a domain with one moving wall. We assume that the wall remains vertical during its motion. Moreover, at least in a vicinity of the moving wall the bathymetry has to be flat to allow its free motion under the action of waves or to follow a prescribed trajectory. The fluid flow is assumed to be potential and we address the complete (i.e. fully nonlinear and fully dispersive) water wave problem [73]. We propose first a reformulation of this problem on a fixed and domain (a square) by introducing a new curvilinear coordinate system. Then, we propose a robust finite difference discretization of this problem with mathematically proven good qualities. The performance of this algorithm is illustrated on several examples of practical interest. In particular, we study the influence of springs rigidity on the oscillatory motion of the wave/wall system.

The results presented in this study can be viewed under a different angle and, thus, they can be applied to another problem of coastal structures protection from wave loads and impacts. In particular, the catastrophic consequences of 2011 Tohoku earthquake and tsunami are widely known nowadays [61]. The protecting structures are made very solid by implementing various security coefficients. However, more economic protections can be designed if we allow them to move under the action of waves. The moving parts can be related to fixed ones by a system of flexible strings and by tuning their rigidity one can reduce significantly the wave run-up on moving parts as well as wave loads on fixed solid structural elements.

We would like to mention also related works. The generation of water waves by an accelerating moving vertical plate in a channel of constant depth was studied in [77]. The Authors measured the free surface elevation and pressure distributions on the wall for three different values of plate’s accelerations. In another recent work [76] a nonlinear model for incompressible free surface flows was proposed. Then, this model was used to study wave run-up on a sloping beach and on a moving vertical wall. The two-dimensional unsteady problem of wave interaction with a moving vertical wall was studied analytically in [48]. Moreover, the comparisons with numerical solutions obtained with the complex boundary element method was presented as well. Another numerical model was presented in [40] as well. In last two works the waves were generated by an initial disturbance located close to the moving wall.

The present manuscript is organized as follows. The physical and mathematical problems are formulated in Section 2 and then it is reformulated on a fixed domain in Section 3. The proposed finite difference approximation is described in Section 4 and studied mathematically in Section 4.4. The construction of the mapping in the fixed domain is discussed in Section 5.5 and the numerical algorithm in general is described in Section 5. Some numerical results are presented in Section 6. Finally, the main conclusions and perspectives are outlined in Section 7.

Problem formulation

Consider a two-dimensional motion of an incompressible, homogeneous and ideal fluid with free surface. The sketch of the fluid domain is depicted in Figure 1. The Cartesian coordinate system $O\,x\,y\,$, $\boldsymbol{x}\ =\ (x,\,y)$ is chosen such that the horizontal axis $O\,x$ is directed along the undisturbed free surface and the axis $O\,y$ points vertically upwards. The fluid domain $\Omega\,(t)$ is assumed to be simply connected and bounded from below by a horizontal bottom $\Gamma_{b}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \{\boldsymbol{x}\ \in\ \Omega\,(t)\ |\ y\ =\ -h_{0}\}$ and from above by the free surface — $\Gamma_{s}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \{\boldsymbol{x}\ \in\ \Omega\,(t)\ |\ y\ =\ \eta(x,\,t)\}\,$. On the sides the domain $\Omega\,(t)$ is bounded by two vertical walls. The right wall $\Gamma_{r}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \{\boldsymbol{x}\ \in\ \Omega\,(t)\ |\ x\ =\ \ell\}$ is assumed to be fixed111Thus, the length of the undisturbed domain is equal to $\ell\ >\ 0\,$., while the left one $\Gamma_{\ell}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \{\boldsymbol{x}\ \in\ \Omega\,(t)\ |\ x\ =\ s\,(t)\}$ is connected to a system of springs and can move horizontally under the wave action. By $s\,(t)$ we denote the displacement of the moving wall222We assume additionally that the wall is non-deformable and remains vertical during the interaction with waves. with respect to its undisturbed position $x\ =\ 0\,$. The closed domain $\Omega\,(t)$ together with its boundaries is denoted by $\bar{\Omega}\,(t)\,$.

Remark 1**.**

In order to be able to simulate water waves with a shape, which is not necessarily a graph, we can adopt a parametric representation of the free surface $\Gamma_{s}\,$:

[TABLE]

where $q\ \in\ \mathds{R}$ is a real parameter and $t$ is the time evolution variable. Functions $\xi\,(\cdot,\,\cdot)\,:\ \mathds{R}\times\mathds{R}^{\,+}\ \mapsto\ \mathds{R}$ and $\eta\,(\cdot,\,\cdot)\,:\ \mathds{R}\times\mathds{R}^{\,+}\ \mapsto\ \mathds{R}$ are assumed to be sufficiently smooth. Without any restriction we can assume that $q\ \in\ [\,0,\,1\,]\,$. Below in the text (see Section 3) this parametrization will appear more naturally when we map the unknown fluid domain $\Omega\,(t)$ to a fixed reference domain.

The fluid flow is described by Euler equations [52]:

[TABLE]

where $\boldsymbol{u}\ =\ (u,\,v)$ is the velocity vector, $\rho\ =\ \mathrm{const}\ >\ 0$ is the fluid density, $p$ is the total pressure and $\mathbf{g}\ =\ (0,\,-\,g)$ is the gravity acceleration. The subscripts such as $\boldsymbol{u}_{\,t}$ denote partial derivatives, i.e. $\boldsymbol{u}_{\,t}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \frac{\partial\/\boldsymbol{u}}{\partial\/t}\,$. These equations have to be supplemented by appropriate boundary conditions in order to have a well-posed problem. The impermeability condition on fixed solid surfaces (the right vertical wall $\Gamma_{r}$ and the bottom $\Gamma_{b}$) reads

[TABLE]

where $\boldsymbol{n}$ is the vector of unitary exterior normal to the corresponding boundary of the domain $\Omega\,(t)\,$. The left wall is impermeable as well and the boundary condition reads

[TABLE]

Taking into account the horizontal character of the wall motion, this condition can be further simplified to give

[TABLE]

Traditionally on the free surface we prescribe the kinematic

[TABLE]

and dynamic

[TABLE]

boundary conditions [73]. Here $p_{a}$ is the constant atmospheric pressure. The boundary conditions mentioned above are enough to obtain a closed system of equations which describe the motion of an ideal fluid. One has to prescribe also the initial conditions, but their form is highly problem-dependent. Consequently, it will not be discussed at the current level of generality.

Remark 2**.**

Let us obtain also kinematic free surface boundary conditions when the free surface $\Gamma_{s}$ is given in a parametric form (2.1). The free surface $\Gamma_{s}$ is characterized by the fact that its velocity is determined by the velocity of fluid particles constituting this boundary. In other words, the point $\boldsymbol{x}\ \in\ \Gamma_{s}$ moves with the velocity of the fluid particle located in this point, i.e.

[TABLE]

By taking the material derivative of both sides of the parametric form (2.1), we have

[TABLE]

Thus, by taking into account the physical information (2.6) we obtain the following kinematic boundary conditions for $\forall\,\boldsymbol{x}\ \in\ \Gamma_{s}$ and for $\forall\,t\ \geqslant\ 0\,$:

[TABLE]

In the conditions above the quantity $\frac{\mathrm{d}\/q}{\mathrm{d}\/t}$ is the rate of change of the parameter $q$ for fluid particles located at the free surface. If one uses the Lagrangian description of a flow, fluid particles keep their labels and, thus, $\frac{\mathrm{d}\/q}{\mathrm{d}\/t}\ \equiv\ 0\,$. However, we use an arbitrary parametrization and in our case $\frac{\mathrm{d}\/q}{\mathrm{d}\/t}\ \neq\ 0$ in general.

Potential flow model

The formulation presented above is still too general. We shall simplify it further by assuming the flow to be irrotational, i.e.

[TABLE]

It implies the existence of a function $\varphi\,:\ \mathds{R}^{2}\times\mathds{R}^{+}\ \mapsto\ \mathds{R}$ which is called the velocity potential such that the velocity field is given by

[TABLE]

Substituting this form into the incompressibility condition (2.2) yields that the velocity potential has to be a harmonic function

[TABLE]

The kinematic boundary condition is obtained straightforwardly by substituting into (2.4) the velocity components by their representations in terms of the velocity potential:

[TABLE]

We note also that velocity components $u$ and $v$ in the kinematic boundary conditions (2.7), (2.8) are expressed in terms of the velocity potential $\varphi$ according to (2.9).

The momentum equation (2.3) can be integrated and combined with the dynamic boundary condition (2.5) to give the so-called Euler–Lagrange integral333The Euler–Lagrange integral becomes the Bernoulli integral for steady flows. We implicitly chose the gauge where the Bernoulli ‘constant’ is identically zero. evaluated at the free surface:

[TABLE]

On solid stationary walls the velocity potential satisfies the homogeneous Neumann condition

[TABLE]

(with $\boldsymbol{n}$ defined above) while on the moving wall $\Gamma_{\ell}$ it satisfies the following non-homogeneous Neumann condition

[TABLE]

So, instead of seeking for two components of the velocity field \boldsymbol{u}(\boldsymbol{x},\,t)\ =\ \bigl{(}u(\boldsymbol{x},t),\,v(\boldsymbol{x},t)\bigr{)} and the pressure function $p(\boldsymbol{x},\,t)\,$, the formulation was simplified to one unknown function, i.e. the velocity potential $\varphi(\boldsymbol{x},\,t)\,$, thanks to the irrotationality assumption. The potential flow formulation with free surface is known as the full water wave problem and it was shown in numerous studies to be an excellent model for water waves [73, 41] (see [12] for water wave problem history review).

2.1.1  Initial conditions

In order to obtain a well-posed problem, we have to specify also the appropriate initial conditions. The free surface elevation in the parametric form is initially given by two functions:

[TABLE]

Let us assume that the initial distribution of fluid particles is known as well:

[TABLE]

Then, we can construct an initial condition for the velocity potential $\varphi$ in the following way:

[TABLE]

where $\varphi_{\,0}\,(\boldsymbol{x}_{\,0})$ is the value of the velocity potential in an arbitrary point $\boldsymbol{x}_{\,0}\ \in\ \bar{\Omega}\,(0)$ and $\mbox{\textgamma}_{\boldsymbol{x}_{\,0}\,\leadsto\boldsymbol{x}}\ \subseteq\ \bar{\Omega}\,(0)$ is an arbitrary piecewise smooth curve connecting points $\boldsymbol{x}_{\,0}$ with $\boldsymbol{x}\,$. We assume that the initial distribution of the velocities $\boldsymbol{u}_{\,0}\,(\boldsymbol{x})$ is potential. Thus, the curvilinear integral above does not depend on the path $\mbox{\textgamma}_{\boldsymbol{x}_{\,0}\,\leadsto\boldsymbol{x}}\,$.

Piston motion

In order to describe the horizontal motion of the piston we adopt the following very simple model based on the second law of Newton [62, 50]:

[TABLE]

where $m$ is the wall mass and $k$ is the stiffness coefficient of springs. Finally, $\mathscr{F}\,(t)$ is the force acting on the left moving wall. It can be computed by integrating the contributions of the whole water column

[TABLE]

If initially the fluid was at rest, the force $\mathscr{F}(0)$ consists only of the hydrostatic loading on the wall, i.e.

[TABLE]

The pressure $p\,(\boldsymbol{x},\,t)\ \equiv\ p\,(x,\,y,\,t)$ can be computed at any point inside the fluid thanks to the Euler–Lagrange integral:

[TABLE]

A more detailed derivation of equation (2.10) can be found in Appendix D.

The second order nonlinear and non-autonomous Ordinary Differential Equation (ODE) (2.10) has to be completed with two initial conditions. If nothing is indicated, we start the integration from the rest state, i.e.

[TABLE]

Dimensionless problem

Equations given in previous Section can be further simplified by choosing appropriate scaled variables. Namely, we introduce the following scaled independent

[TABLE]

and dependent

[TABLE]

variables. The dimensional coefficients and parameters appearing in governing equations scale as follows:

[TABLE]

Finally, dimensionless governing equations read (where for simplicity we drop out all asterisk $*$ symbol from superscripts):

[TABLE]

In numerical simulations below we shall solve this dimensionless system of equations. It is equivalent to setting simply $h_{0}\ =\ 1$ and $g\ =\ 1$ in the computer code.

Equations on a fixed domain

The main difficulty of the problem described above is that the computational domain $\Omega\,(t)$ is time dependent. First of all, it is due to the motion of the free surface, but also due to the motion of the left wall. Consequently, the strategy adopted in this study consists in transforming the problem to a fixed domain $\mathcal{Q}^{\,0}\,$. Obviously, this transformation will be time dependent due to the unsteady character of the problem. In the past; the idea of using conformal mappings from Complex Analysis has become popular in 2D Hydrodynamics [52]. For the 2D water wave problem it was proposed by L. Ovsyannikov (1974) [64] and implemented much later numerically by A. Dyachenko et al. (1996) [29]. Strictly speaking, in our developments we do not need the conformal property of the underlying mapping. Consequently we shall relax this assumption. Moreover, it will allow us to have a fixed and finite transformed domain $\mathcal{Q}^{\,0}\,$. For instance, we can choose $\mathcal{Q}^{\,0}$ as a unit square

[TABLE]

and we consider a smooth bijective mapping (i.e. a diffeomorphism) $\boldsymbol{x}\,:\ \mathcal{Q}^{\,0}\ \mapsto\ \Omega\,(t)\,$:

[TABLE]

[TABLE]

Additionally we assume that the Jacobian matrix $\mathcal{J}\,(\boldsymbol{q},\,t)$ of transformation (3.1) is non-singular, i.e.

[TABLE]

This situation is schematically depicted in diagram (3.2) and in Figure 2. For more details on the computation of partial derivatives we refer to Appendix A.

It is assumed that left $\Gamma_{\ell}$ and right $\Gamma_{r}$ walls are images of left and right sides of the square $\mathcal{Q}^{\,0}\,$, i.e. $q^{\,1}\ =\ 0$ and $q^{\,1}\ =\ 1$ correspondingly ($0\ \leqslant\ q^{\,2}\ \leqslant\ 1$). Similarly, the upper ($q^{\,2}\ =\ 1$) and lower ($q^{\,2}\ =\ 0$) sides of the square are respectively mapped on the free surface $\Gamma_{s}$ and bottom $\Gamma_{b}\,$. The boundary of the square $\mathcal{Q}^{\,0}$ will be denoted by $\gamma\,$, i.e.

[TABLE]

The practical construction of the mapping $\boldsymbol{x}\,(\boldsymbol{q},\,t)$ is described in Section 5.5. At the present stage of the exposition we shall assume that mapping $\boldsymbol{x}\,(\boldsymbol{q},\,t)$ is simply known.

In this Section we shall skip the intermediate computations, which are rather standard in the Differential Geometry, for example. We detail them in Appendices A and B. The system of governing equations (2.12) – (2.20) posed on the fixed domain $\mathcal{Q}^{\,0}$ reads:

[TABLE]

where the implicit convention about the summation over repeated indices is adopted. The transformation of the Laplace equation is detailed in Appendix B.

Remark 3**.**

If the free surface is given in the parametric form (2.1), then we have to add one more kinematic free surface boundary condition:

[TABLE]

See also Appendix C for some hints on the derivation of boundary conditions in the transformed domain.

Notice that equation (2.18) is not affected by transformation (3.1). That is why we do not repeat it here. Above we introduced the following notations. The coefficients in transformed Laplacian are defined as

[TABLE]

where $\mathcal{G}=(g_{\alpha\beta})_{\alpha,\beta\,=\,1,\,2}$ is the familiar metric tensor

[TABLE]

In fact, tensor $\mathcal{K}\ =\ (\mathcal{K}_{\,\alpha\beta})_{\alpha,\,\beta\,=\,1,2}$ can be expressed in a compact matrix form through the metric tensor $\mathcal{G}$ as

[TABLE]

The Cartesian $(u,\,v)$ and the contravariant $v^{\,1}$ components of the velocity field are defined as

[TABLE]

In order to compute the derivatives of the inverse mapping, we use the following expressions:

[TABLE]

[TABLE]

The mathematical problem of a potential flow formulation with free surface in curvilinear coordinates consists in determining the velocity potential $\varphi\,(q^{\,1},\,q^{\,2},\,t)$ and the free surface profile $\eta\,(q^{\,1},\,t)\,$. The velocity potential satisfies in the domain $\mathcal{Q}^{\,0}$ an elliptic equation (3.4) with non-constant coefficients along with boundary conditions (3.5) – (3.9) on $\partial\mathcal{Q}^{\,0}\,$. In order to obtain a well-posed problem we have to complete the formulation with two initial conditions at $t\ =\ 0$ directly in the fixed domain:

[TABLE]

Earlier the initial conditions were discussed in Section 2.1.1. They can be transposed on the fixed domain $\mathcal{Q}^{\,0}$ using the inverse mapping $\boldsymbol{x}^{\,-1}\,(\boldsymbol{q},\,t)\,$.

On some properties of elliptic operators in curvilinear coordinates

It can be noticed that the Laplace equation in curvilinear coordinates (3.4) has a more complex form comparing to the initial Cartesian coordinates $O\,x\,y\,$. In particular, constant coefficients become variable in space and time. Moreover, the mixed derivatives appear as well. However, some properties are important to construct qualitatively correct numerical discretizations. For example, below we shall construct a finite difference scheme for equation (3.4) with a positive definite difference operator. The proof of this fact relies heavily on the uniform ellipticity property of operator (3.4).

Lemma 1**.**

Partial differential equation (3.4) is uniformly elliptic.

Proof.

By our assumptions the mapping (3.1) is differentiable with bounded derivatives in $\mathcal{Q}^{\,0}$ and non-degenerate with a positive Jacobian $\mathcal{J}\ >\ 0\,$, $\forall\boldsymbol{q}\ \in\ \mathcal{Q}^{\,0}\,$, $\forall t\ >\ 0\,$. Thanks to the identity

[TABLE]

the metric components $g_{\,1\,1}$ and $g_{\,2\,2}$ take strictly positive values. Thanks to definitions (3.13), it is straightforward to conclude that $\mathcal{K}_{\,1\,1}\,$, $\mathcal{K}_{\,2\,2}$ are also positive. The matrix

[TABLE]

is symmetric and its determinant is equal to $\det\mathcal{K}\ =\ 1$ due to (3.13). Consequently, it possesses two orthonormal eigenvectors corresponding to eigenvalues

[TABLE]

where the discriminant is defined as

[TABLE]

Then, for any real numbers $\forall\xi\,$, $\zeta\ \in\ \mathds{R}$ and for any point $\forall\boldsymbol{q}\ \in\ \mathcal{Q}^{\,0}$ we have the following inequalities:

[TABLE]

where $\mathscr{Q}$ is the quadratic form defined as

[TABLE]

The constants $c_{\,1,\,2}$ are defined as

[TABLE]

The positivity of $c_{\,1}\ >\ 0$ implies the uniform ellipticity property. ∎

The constants $c_{\,1,\,2}$ defined in the proof above can be used to estimate the convergence speed of iterative methods, which depends directly on the conditioning number [19] of the matrix corresponding to the difference operator [69]. The conditioning number depends on the ratio $\dfrac{c_{2}}{c_{1}}\,$. High ratio implies a poorly conditioned difference operator and, thus, more iterations are needed to converge to the solution within prescribed accuracy. We underline also that constants $c_{\,1,\,2}$ depend only on the mapping (3.1). Below we provide two examples which illustrate situations where the ratio $\dfrac{c_{2}}{c_{1}}\ \gg\ 1\,$.

3.1.1  Example 1

Let us take a physical domain in the form of a rectangle:

[TABLE]

with substantially different sizes of the walls, e.g. we can assume $0\ <\ \ell_{\,1}\ \ll\ \ell_{\,2}\,$. The mapping (3.1) from $\mathcal{Q}^{\,0}\ \mapsto\ \Omega_{\,\blacksquare}$ is given explicitly by formulas

[TABLE]

Then, we can compute explicitly the metric coefficients

[TABLE]

the Jacobian

[TABLE]

and eigenvalues of the matrix $\mathcal{K}\,$:

[TABLE]

Thus, the conditioning of the finite difference operator will scale with

[TABLE]

Thus, for highly distorted domains the ratio $\dfrac{c_{\,2}}{c_{\,1}}$ becomes large. If on the reference domain $\mathcal{Q}^{\,0}$ we use a uniform square grid, than mapping (3.24) generates a uniform grid with distorted cells. The convergence of iterative methods on such grids will be slowed down. Thus, in real computations we should avoid highly distorted cells since they will penalize the convergence of linear solvers.

3.1.2  Example 2

Let us take another physical domain $\Omega_{\,\square}$ having the shape of a parallelogram. It can be obtained as an image of the reference domain $\mathcal{Q}^{\,0}$ under the following mapping

[TABLE]

where $0\ <\ \mbox{\textpsi}\ \ll\ \dfrac{\pi}{2}$ is a small angle. Then, the metric coefficients of this mapping are

[TABLE]

The Jacobian $\mathcal{J}\ =\ \sin\mbox{\textpsi}$ and eigenvalues of the matrix $\mathcal{K}$ are

[TABLE]

and henceforth

[TABLE]

So, we showed that the ratio $\dfrac{c_{\,2}}{c_{\,1}}$ becomes large for domains $\Omega_{\,\square}$ featuring a small angle. The problem is that the grid generated by mapping (3.25) is substantially non-orthogonal, since it consists of parallelograms with an acute angle \textpsi . On such grids we should expect some reduction of iterative methods convergence speed. Thus, in our computations we should avoid grids featuring small angles.

Finite difference scheme in curvilinear coordinates

In the numerical simulation of free surface potential flows of an ideal fluid, an elliptic equation to determine the velocity potential $\varphi$ has to be solved at every time step. According to our numerical algorithm (see Section 5 for more details), we determine first the velocity potential value on the free surface $\Gamma_{s}$, then using this value we reconstruct the velocity potential $\varphi$ in the whole fluid domain (by taking into account other boundary conditions on $\Gamma_{b}\,$, $\Gamma_{r}$ and $\Gamma_{\ell}$). So, in this Section we assume that we know the values of the velocity potential in nodes which constitute the pre-image of the free surface $\Gamma_{s}\,$. The goal is to determine the values of the velocity potential in all remaining grid nodes. These values will be determined by solving a system of difference equations constructed below.

Moreover, we assume that the curvilinear grid $\Omega_{\,h}^{\,n}$ on the n$${}^{\mathrm{\small\textsf{th}}} time layer is already constructed. The nodes $\boldsymbol{x}_{\,\boldsymbol{j}}^{\,n}$ are images under the mapping (3.1) of fixed nodes $\boldsymbol{q}_{\,\boldsymbol{j}}\,$, which constitute the uniformly distributed grid $\mathcal{Q}^{\,0}_{\,h}\ \subseteq\ \mathcal{Q}^{\,0}\,$. Here $\boldsymbol{j}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ (j_{\,1},\,j_{\,2})$ is a multi-index and \boldsymbol{q}_{\,\boldsymbol{j}}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \bigl{(}q_{\,j_{\,1}}^{\,1},\,q_{\,j_{\,2}}^{\,2}\bigr{)}\,. The uniform grid is traditionally defined as

[TABLE]

We make an additional geometrical (and not very restrictive) assumption on the mapping (3.1): boundary components $\gamma_{\,\ell}\,$, $\gamma_{\,b}\,$, $\gamma_{\,r}$ and $\gamma_{\,s}$ are mapped on corresponding boundary components $\Gamma_{\ell}\,$, $\Gamma_{b}\,$, $\Gamma_{r}$ and $\Gamma_{s}$ of the fluid domain $\Omega(t)$ (see Figure 2 for an illustration). The boundary of the reference domain $\mathcal{Q}^{\,0}$ after the discretization is denoted as

[TABLE]

The sketch of the discretized reference domain $\mathcal{Q}^{\,0}_{\,h}$ is depicted in Figure 3.

In order to construct a finite difference approximation we employ the integro-interpolation method using the terminology of Tikhonov & Samarskii [75]. In the western literature this method is closer to the finite volume/conservative finite difference methods [36, 35]. The choice for finite differences seems to be natural since we would like to solve an elliptic equation on a simple Cartesian domain [68]. For this purpose we replace equation (3.4) by the following integral relation

[TABLE]

where $\mathcal{C}\ \subseteq\ \mathcal{Q}^{\,0}$ is an arbitrary sub-domain with a piece-wise smooth boundary. We introduced above the notation

[TABLE]

Below we construct discrete approximations based on this integral formulation.

Remark 4**.**

Please, notice that the transformed Laplace equation (3.4) can be compactly rewritten using the new notation:

[TABLE]

Approximation in interior nodes

Let us consider an interior grid $\mathcal{Q}^{\,0}_{\,h}$ node marked with symbol [math] in Figure 4(a) along with the integration contour $A\,B\,C\,D\,$. For this contour the integral equation (4.1) takes the form

[TABLE]

Now, by applying a quadrature rule to all integrals above, one can obtain a discrete analogue of this integral relation. In the present study we employ the quadrature formula of rectangles. For the integral over the segment $B\,C$ we have

[TABLE]

where we used the grid numbering shown in Figure 4(a) and we introduced the following notation

[TABLE]

for left- and right-sided finite difference operators. By using similar approximations to discretize other integrals over sides $A\,D\,$, $D\,C$ and $A\,B$ we obtain a fully discrete analogue of the integral equation (4.1):

[TABLE]

By dividing the both parts of the last equation by the area $h_{\,1}\cdot h_{\,2}$ of the rectangle $A\,B\,C\,D\,$, we obtain more compact difference equations in any interior node $\boldsymbol{q}_{\,\boldsymbol{j}}\ \in\ \mathcal{Q}_{\,h}^{\,0}\,$:

[TABLE]

These equations approximate the original differential equation (3.4) to the order \mathcal{O}\bigl{(}h_{\,1}^{\,2}\ +\ h_{\,2}^{\,2}\bigr{)} provided that solutions are sufficiently smooth. This is achieved by using the $9-$point stencil shown in Figure 4(a).

Approximation of boundary conditions

For the sake of manuscript conciseness we explain the boundary conditions treatment on the example of condition (3.8) imposed on the fixed444The left wall might be moving in the physical space. However, it becomes fixed in the transformed domain. boundary $\gamma_{\,\ell}\,$. Let \gamma_{\,\ell}^{\,h}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \bigl{\{}\boldsymbol{q}_{\,0,\,j_{\,2}}\ \in\ \gamma_{\,h}\ |\ j_{\,2}\ =\ 1,\,2,\,\ldots,\,N_{2}-1\bigr{\}} be the set of ‘interior’ grid nodes belonging to the left boundary $\gamma_{\,\ell}$ except two angular nodes $\boldsymbol{q}_{\,0,\,N_{2}}$ and $\boldsymbol{q}_{\,0,\,0}\,$, which deserve a special consideration. For a boundary node $\boldsymbol{q}_{\,\boldsymbol{j}}\ \in\ \gamma_{\,\ell}^{\,h}$ we choose the integration contour as it is shown in Figure 4(b). The quadrature rule for the integral over side $B\,C$ is given in equation (4.3) (one has just to substitute $j_{\,1}\ \hookleftarrow\ 0$). On the segment $A\,D$ we have to use the boundary condition (3.8) that we can rewrite in a compact form:

[TABLE]

Using this information, we have directly the following approximation:

[TABLE]

For the side $D\,C$ the formula of rectangles yield the following approximation of the integral:

[TABLE]

For the segment $A\,B$ we give only the final result without intermediate calculations:

[TABLE]

Now we substitute all these approximations into the integral equation (4.2) and after dividing its both parts by $h_{\,2}$, we obtain the following difference relation in any node $\gamma_{\,\ell}^{\,h}\,$:

[TABLE]

Here we introduced the following notation:

[TABLE]

This difference equation is inhomogeneous because of the wall motion. If the (left) wall is fixed, then $\mbox{\textmugreek}^{\,0}\ \equiv\ 0$ and this relation becomes homogeneous.

Remark 5**.**

Please, notice that in all formulas in this Section $j_{\,1}\ \equiv\ 0\,$.

4.2.1  Consistency

In this Section we show that equation (4.5) approximates the boundary condition (3.8) to the order $\mathcal{O}(h_{\,1}^{\,2}\ +\ h_{\,2}^{\,2})\,$. Indeed, let us substitute a smooth solution $\varphi\,(q^{\,1},\,q^{\,2},\,t)$ into equation (4.5) and perform local Taylor expansions. Then, the consistency error \textpsi can be computed:

[TABLE]

In the very last expression all quantities are evaluated in the same node $\boldsymbol{q}_{\,0,\,j_{\,2}}\,$, $j_{\,2}\ =\ 1,\,2,\,\ldots,\,N_{2}-1\,$. Moreover, taking into account the boundary condition (4.4) and assuming that Laplace equation (3.4) (see also Remark 4) is fulfilled up to the boundary $\gamma_{\,\ell}\,$, we obtain that

[TABLE]

Thus, the boundary condition is approximated to the second order in space. In a similar way we can construct finite difference approximations in boundary nodes \gamma_{\,b}^{\,h}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \bigl{\{}\boldsymbol{q}_{\,j_{\,1},\,0}\ \in\ \gamma^{\,h}\ |\ j_{\,1}\ =\ 1,\,2,\,\ldots,\,N_{1}-1\bigr{\}} and \gamma_{\,r}^{\,h}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \bigl{\{}\boldsymbol{q}_{\,N_{1},\,j_{\,2}}\ \in\ \gamma^{\,h}\ |\ j_{\,2}\ =\ 1,\,2,\,\ldots,\,N_{2}-1\bigr{\}}\,. In the derivation of these boundary conditions one has to use the following impermeability conditions (compare with conditions (3.7) and (3.9)):

[TABLE]

All these difference equations have six point stencils to achieve the second order accuracy: three nodes lie in the interior of the domain $\mathcal{Q}^{\,0}_{\,h}$ and three on the boundary $\gamma^{\,h}\,$. Finally, there are also two angular nodes $\boldsymbol{q}_{\,0,\,0}$ and $\boldsymbol{q}_{\,N_{1},\,0}\,$. The stencils contain four nodes: one interior and three on the boundaries.

Difference operator

In order to determine the discrete values of the velocity potential \bigl{\{}\varphi_{\,\boldsymbol{j}}\bigr{\}}\,, we constructed above a finite difference problem with the uniform second order accuracy. In this difference problem one has to find the values of the velocity potential $\varphi_{\,\boldsymbol{j}}$ in all nodes except those on the free surface, where the Dirichlet-type condition is imposed:

[TABLE]

In this way we have to determine only \bigl{(}N_{1}\ +\ 1\bigr{)}\cdot N_{2} the values \bigl{\{}\varphi_{\,\boldsymbol{j}}\bigr{\}}_{\,\boldsymbol{j}} in the nodes \bigl{\{}\boldsymbol{x}_{\,j_{\,1},\,j_{\,2}}\bigr{\}} with $j_{\,1}\ =\ 0,\,1,\,\ldots,\,N_{1}$ and $j_{\,2}\ =\ 0,\,1,\,\ldots,\,N_{2}-1\,$. In the resulting problem we have precisely \bigl{(}N_{1}\ +\ 1\bigr{)}\cdot N_{2} equations.

In order to study theoretically the finite difference problem, it will be convenient to make a change of variables in order to get homogeneous Dirichlet’s boundary conditions in nodes $\boldsymbol{q}_{\,j_{\,1},\,N_{2}}\,$, $j_{\,1}\ =\ 0,\,1,\,\ldots,\,N_{1}\,$. For this purpose we introduce a new grid function

[TABLE]

and a new dependent variable

[TABLE]

Then, in the nodes $\boldsymbol{q}_{\,\boldsymbol{j}}\ \in\ \gamma^{\,h}_{\,s}$ on the free surface the function $\mbox{\textphi}_{\,\boldsymbol{j}}$ will take zero values. The difference problem for \textphi will differ from the one for $\varphi$ by inhomogeneous right hand sides in the nodes close to the free surface pre-image $\gamma^{\,h}_{\,s}\,$. For the new problem we introduce the operator notation

[TABLE]

where \textnu is the vector of the right hand sides. The difference operator $\Lambda$ can be decomposed in sub-operators for our convenience:

[TABLE]

and on one more level:

[TABLE]

Operators \bigl{\{}\Lambda_{\,\alpha\,\beta}\bigr{\}}_{\,1\,\leqslant\,\alpha,\,\beta\,\leqslant\,2} are defined below as follows.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Above we used the notation

[TABLE]

with $\alpha,\,\beta\ =\ 1,\,2\,$. Moreover, all the operators vanish in the nodes of the upper boundary $\gamma_{\,s}^{\,h}\,$, i.e.

[TABLE]

Below we study the properties of the operator $\Lambda$ to show that the problem (4.6) is well-posed.

Mathematical study of the finite difference problem for the velocity potential

Let us define by $\mathcal{H}$ the set of functions defined on the grid $\mathcal{Q}^{\,0}_{\,h}$ and which take zero values on $\gamma_{\,s}^{\,h}\,$, i.e.

[TABLE]

We introduce also on space $\mathcal{H}$ the scalar product

[TABLE]

This scalar product induces a natural norm

[TABLE]

Let us define a new operator $\mathfrak{A}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ -\Lambda$ and study its properties. The operator $\Lambda$ was introduced in (4.6) and rigorously defined earlier.

For any fixed value of index $j_{\,2}\,$, a function $u\ \in\ \mathcal{H}$ can be considered as a function being defined on a 1D grid:

[TABLE]

The class of all grid functions defined on $\mathcal{Q}_{\,h}^{\,1}$ will be denoted by $\mathcal{H}_{\,1}\,$. For two arbitrary functions $u,\,v\ \in\ \mathcal{H}$ we define the scalar product on the layer $j_{\,2}$ as

[TABLE]

Then, it is not difficult to see that

[TABLE]

However, we can proceed in a different way by considering instead ‘vertical’ slices. Let $\mathcal{H}_{\,2}$ be the set of 1D functions defined on the grid

[TABLE]

and taking zero value for $j_{\,2}\ =\ N_{\,2}\,$, i.e.

[TABLE]

Similarly, in $\mathcal{H}$ we can introduce a scalar product on the layer $j_{\,1}$ as follows:

[TABLE]

and the scalar product in $\mathcal{H}$ can be similarly expressed as

[TABLE]

It is not difficult to show that formulas (4.16) and (4.18) provide the same value for $\langle\,u,\,v\,\rangle_{\,\mathcal{H}}$ and, thus, they are nothing else than two different ways to express the scalar product in $\mathcal{H}\,$. Intuitively, we can understand it in the following way. Imagine that one wants to compute the sum of all elements of a matrix. To do it, one can move along the lines and then along the columns, which corresponds to (4.18), or vice versa (4.16). Since the addition is commutative, the resulting sum will be the same.

4.4.1  Some discrete identities

Let us remind that for any discrete functions \bigl{\{}\mbox{\textalpha}_{\,j}\bigr{\}}_{\,j\,=\,0}^{\,N}\,, \bigl{\{}u_{\,j}\bigr{\}}_{\,j\,=\,0}^{\,N} and \bigl{\{}v_{\,j}\bigr{\}}_{\,j\,=\,0}^{\,N} defined on a 1D grid \mathcal{Q}_{\,h}^{\,1}\ =\ \bigl{\{}\,q_{\,j}\ \in\ [\,0,\,1\,]\,\bigr{\}}_{\,j\,=\,0}^{\,N}\,, the first discrete Green identity can be shown:

[TABLE]

and also

[TABLE]

where for the sake of simplicity we introduced the following notations

[TABLE]

[TABLE]

Moreover, the following summation by parts formulas hold as well

[TABLE]

All these properties mentioned above show that the discretization proposed in the present study belongs to the class of operational finite difference schemes [70], which became later known as mimetic methods [54].

The main result

Theorem 1**.**

The operator $\mathfrak{A}\ \equiv\ -\Lambda$ is self-adjoint in the space $\mathcal{H}\,$.

Proof.

For two arbitrary grid functions $u,\,v\ \in\ \mathcal{H}$ we consider their scalar product on a layer $0\ \leqslant\ j_{\,2}\ <\ N_{2}\,$:

[TABLE]

Taking into account definition (4.7) of the operator $\Lambda_{\,1\,1}$ and Green’s identity (4.19) we obtain

[TABLE]

or simply

[TABLE]

By substituting the last expression into formula (4.16), we obtain the equality

[TABLE]

The last identity is absolutely symmetric with respect to the functions $u$ and $v\,$, i.e.

[TABLE]

Consequently, the operator $\Lambda_{\,1\,1}$ is self-adjoint in space $\mathcal{H}$, i.e.

[TABLE]

Using the summation by parts formulas (4.21), (4.22) along with the definition (4.8) of the operator $\Lambda_{\,1\,2}$ in interior and boundary nodes, we obtain the following expression for the scalar product

[TABLE]

It is obvious that

[TABLE]

Thus, the operator $\Lambda_{\,1\,2}$ is not self-adjoint. Nevertheless, we continue the proof without any deception. Consider right now the operator $\Lambda_{\,2\,1}$ defined in (4.9). Using similar techniques we obtain the following expression for the scalar product of $\Lambda_{\,2\,1}\,u$ with $v\,$:

[TABLE]

From the last formula it can be clearly seen that

[TABLE]

Thus, the operator $\Lambda_{\,2\,1}$ as well as $\Lambda_{\,1\,2}$ is not self-adjoint. However, we can study the sum $\Lambda_{\,2\,1}\ +\ \Lambda_{\,1\,2}\,$. By regrouping the terms in (4.24) and (4.25) we obtain

[TABLE]

Above we used also the fact that $\mathcal{K}_{\,1\,2}\ \equiv\ \mathcal{K}_{\,2\,1}\,$. From the last formula (4.26) we readily obtain that

[TABLE]

Thus, the operator $\Lambda_{\,2\,1}\ +\ \Lambda_{\,1\,2}$ is self-adjoint in the space $\mathcal{H}\,$, i.e.

[TABLE]

even if constitutive operators $\Lambda_{\,1\,2}$ and $\Lambda_{\,2\,1}$ are not self-adjoint if taken separately.

Finally, consider the operator $\Lambda_{\,2\,2}\,$, which was defined in (4.10). We employ first the discrete Green identities (4.19), (4.20) to evaluate scalar products on layers j_{\,1}\ \in\ \bigl{\{}0,\,1,\,\ldots,\,N_{1}\bigr{\}}\,:

[TABLE]

Using this partial result we can now compute the scalar product in the $\mathcal{H}$ space as well:

[TABLE]

By changing symmetrically grid functions $u$ and $v\,$, we obtain

[TABLE]

Consequently, the operator $\Lambda_{\,2\,2}$ is self-adjoint in the space $\mathcal{H}$, i.e.

[TABLE]

Since \Lambda\ \equiv\ \Lambda_{\,1\,1}\ +\ \bigl{(}\Lambda_{\,1\,2}\ +\ \Lambda_{\,2\,1}\bigr{)}\ +\ \Lambda_{\,2\,2} and each term is an self-adjoint operator, the operator $\Lambda$ is self-adjoint as well. Hence, the same holds for $\mathfrak{A}\ \equiv\ -\Lambda\,$. ∎

Below we are going to show also that the operator $\mathfrak{A}$ is positive definite. The proof will be based on the uniform ellipticity property of equation (3.4), which was established in Lemma 1. Another key ingredient consists in the properties of the finite difference approximation to the following mixed Boundary Value Problem (BVP) for the Laplace equation posed on $\mathcal{Q}_{\,0}\,$:

[TABLE]

with boundary conditions

[TABLE]

This problem is a particular case of the BVP considered earlier for elliptic equation (3.4) with the following coefficients:

[TABLE]

In this case the finite difference operator $\Lambda$ admits a much simpler form since finite difference operators corresponding to mixed derivatives vanish, i.e.

[TABLE]

The finite difference scheme for the BVP (4.28) – (4.30) can be written as

[TABLE]

with operators $\Delta_{\,h}^{\,(1,\,2)}$ defined as

[TABLE]

[TABLE]

Below we introduce also the operators $\mathfrak{B}_{\,1}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ -\Delta_{\,h}^{\,(1)}\,$, $\mathfrak{B}_{\,2}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ -\Delta_{\,h}^{\,(2)}$ and $\mathfrak{B}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \mathfrak{B}_{\,1}\ +\ \mathfrak{B}_{\,2}\,$, i.e. $\mathfrak{B}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ -\,\Delta_{\,h}\,$. By applying Theorem 1 we conclude that the operator $\mathfrak{B}\,:\ \mathcal{H}\ \mapsto\ \mathcal{H}$ is self-adjoint in the space $\mathcal{H}\,$.

It is not difficult to see that the family of grid functions

[TABLE]

satisfies the following identities:

[TABLE]

Henceforth, the functions $\{\mbox{\textpsi}^{\,(k,\,l)}\}_{\,k,\,l}$ are eigenfunctions555Please, notice that each eigenfunction $\mbox{\textpsi}^{\,(k,\,l)}$ is a grid function taking in each node \boldsymbol{j}\ =\ \bigl{(}j_{\,1},\,j_{\,2}\bigr{)} the value $\mbox{\textpsi}^{\,(k,\,l)}_{\,j_{\,1},\,j_{\,2}}$ given in equation (4.31). of the operator $\mathfrak{B}$ and numbers $\{\lambda^{\,(k,\,l)}\}_{\,k,\,l}\ \in\ \mathds{R}$ are its eigenvalues:

[TABLE]

If we choose the normalization factor $\alpha_{\,k}$ as

[TABLE]

then the eigenfunctions $\{\mbox{\textpsi}^{\,(k,\,l)}\}_{\,k,\,l}$ of the operator $\mathfrak{B}$ form an orthonormal basis in $\mathcal{H}\,$, i.e.

[TABLE]

In other words, an arbitrary function $u\ \in\ \mathcal{H}$ can be represented as a sum of a finite Fourier series:

[TABLE]

where the numbers \bigl{\{}c_{\,k\,l}\ \equiv\ \langle\,u,\,\mbox{\textpsi}^{\,(k,\,l)}\,\rangle_{\,\mathcal{H}}\bigr{\}}_{\,k,\,l} are called the Fourier coefficients. The Parseval theorem holds in the finite dimensional setting as well:

[TABLE]

Since all eigenvalues of the operator $\mathfrak{B}$ are positive, the following Lemma holds:

Lemma 2**.**

The finite difference operator $\mathfrak{B}$ is positive definite in space $\mathcal{H}$ and for any function $u\ \in\ \mathcal{H}$ we have the following estimations

[TABLE]

with

[TABLE]

Proof.

Let us take an arbitrary function $u\ \in\ \mathcal{H}$ and we expand it (4.33) on the basis of orthonormal eigenfunctions of the operator $\mathfrak{B}\,$. Then we have

[TABLE]

where we used the fact that $\mbox{\textpsi}^{\,(k,\,l)}$ is an eigenfunction of $\mathfrak{B}$ (4.32). Then, by taking into account the orthonormality property we can greatly simplify the scalar product:

[TABLE]

Obviously, among a finite number of eigenvalues $\{\lambda^{\,(k,\,l)}\}_{\,k,\,l}$ we can always choose the minimal and maximal ones. Thus,

[TABLE]

Taking into account the last estimation along with the Parseval identity (4.34), we obtain immediately the estimations (4.35). ∎

Thanks to the property of self-adjointness and positive definiteness of the operator $\mathfrak{B}\,$, we can define another scalar product $\langle\,u,\,v\,\rangle_{\,\mathfrak{B}}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \langle\,\mathfrak{B}\,u,\,v\,\rangle_{\,\mathcal{H}}\,$, which implies the corresponding energy norm:

[TABLE]

generated by the operator $\mathfrak{B}\,$. Below we use an equivalent representation of the energy norm, which is more suitable for our purposes:

[TABLE]

Lemma 3**.**

For any function $u\ \in\ \mathcal{H}$ the following identities hold:

[TABLE]

[TABLE]

Proof.

Taking into account equation (4.15) and discrete identities shown in Section 4.4.1, we obtain the equality

[TABLE]

Then, from the last result and relation (4.16) follows immediately the requested identity (4.37). Similarly, the discrete identities along with definition (4.17) yield

[TABLE]

Then, using relation (4.18) we obtain the second requested identity (4.38). ∎

Now we can come back to the original finite difference operator $\mathfrak{A}$ and prove the following

Theorem 2**.**

The operator $\mathfrak{A}\ \equiv\ -\,\Lambda$ is positive definite in the space $\mathcal{H}\,$, and we have the following estimations:

[TABLE]

where the constants $c_{\,1,\,2}$ were defined in (3.22), (3.23) and $\lambda_{\,\min,\,\max}$ — in (4.36).

Proof.

As in the first step of the proof we take $v\ \equiv\ u\ \in\ \mathcal{H}$ and substitute it into equations (4.23). After regrouping interior and boundary nodes, we obtain:

[TABLE]

We perform the same operation with equation (4.27) as well:

[TABLE]

Similarly, we set $v\ \equiv\ u$ in equation (4.26) and using the last two formulas we have

[TABLE]

where we used the quadratic form $\mathscr{Q}(\cdot,\,\cdot)$ defined earlier in (3.21). By taking into account the uniform ellipticity property (3.20), we can estimate the quantity $\langle\,\mathfrak{A}\,u,\,u\,\rangle_{\,\mathcal{H}}$ from below:

[TABLE]

Then, in the expression above for $\mathscr{B}_{\,\mathrm{low}}$ we replace666We can do it since the right finite difference can be considered as the left one for the subsequent value of the index. Thus, this operation can be seen as a change of index by one in all summations. all right derivatives $u_{\,q^{\,\alpha}}^{\,\natural}$ by their left counterparts $u_{\,q^{\,\alpha}}^{\,\flat}\,$:

[TABLE]

Now we can make use of Lemma 3, where we showed identities (4.37) and (4.38), which yield:

[TABLE]

In other words, we just showed the following estimation from below:

[TABLE]

Finally, by applying Lemma 2 we obtain that

[TABLE]

The two last inequalities correspond precisely to the required lower bound in (4.39). By departing again from equation (4.40) and using the same techniques we can show the upper bound in (4.39). This completes the proof of this Theorem. ∎

Remark 6**.**

During the proof of Theorem 2 we established also the following estimates:

[TABLE]

or equivalently

[TABLE]

Such operators $\mathfrak{A}$ and $\mathfrak{B}$ satisfying the inequalities above are usually referred to as energetically equivalent [69].

4.5.1  Final remarks

From the positive definiteness property of the operator $\mathfrak{A}$ follows the unique solvability property of the finite difference problem (4.6). To solve this problem numerically one can use almost any iterative technique. It is natural to use the conjugate gradient method whose convergence is guaranteed if the operator is positive definite and self-adjoint [67]. The speed of convergence depends essentially on the conditioning number $\kappa$ of the difference operator $\mathfrak{A}\,$, which can be estimated as

[TABLE]

For larger values of $\kappa$ the convergence will be slower. As we showed earlier in Section 3.1, the first factor $\dfrac{c_{\,2}}{c_{\,1}}$ is completely determined by the properties of the geometric mapping (3.1). This ratio becomes large for the meshes with highly elongated or highly distorted elements. The second factor $(\star)$ depends on the number of nodes $N_{1,\,2}$ chosen in the directions $O\,q^{1,\,2}$ correspondingly and $(\star)$ becomes larger when we refine the grid (i.e. when $N_{1,\,2}$ increase).

Numerical algorithm

In this Section we describe briefly how the numerical code is organized and we provide the details on the treatment of kinematic and dynamic boundary conditions on the free surface and on the lateral fixed and moving walls.

Assume that the position of the wall $x\ =\ s^{\,n}$ on the time layer $t\ =\ t^{\,n}$ is known along with all other fields. The grid $\Omega_{\,h}^{\,n}$ with nodes $\boldsymbol{x}_{\,\boldsymbol{j}}^{\,n}\,$, ordered using a multi-index \boldsymbol{j}\ =\ \bigl{(}j_{\,1},\,j_{\,2}\bigr{)} in the physical space is constructed. On this grid we know the values of grid functions \bigl{\{}\xi_{\,j_{\,1}}^{\,n}\bigr{\}}_{\,j_{\,1}}\,, \bigl{\{}\eta_{\,j_{\,1}}^{\,n}\bigr{\}}_{\,j_{\,1}} and \bigl{\{}\varphi_{\,\boldsymbol{j}}^{\,n}\bigr{\}}_{\,\boldsymbol{j}}\,. The time marching requires to compute these quantities on the following time layer $t\ =\ t^{\,n+1}\,$. It is done in several stages:

1. (1)

First, we compute the values of the velocity potential \bigl{\{}\varphi_{\,j_{\,1},\,N_{2}}^{\,n+1}\bigr{\}}_{\,j_{\,1}} at the free surface node $\boldsymbol{q}_{\,j_{\,1},\,N_{2}}\,$. 2. (2)

These values are used as the Dirichlet data in the linear problem (4.6) to determine the velocity potential \bigl{\{}\varphi_{\,j_{\,1},\,j_{\,2}}^{\,n+1}\bigr{\}}_{\,j_{\,1},\,j_{\,2}} in all other nodes $\boldsymbol{q}_{\,\boldsymbol{j}}$ with $0\ \leqslant\ j_{\,1}\ \leqslant\ N_{1}$ and $0\ \leqslant\ j_{\,2}\ <\ N_{2}\,$. 3. (3)

Then, we find new positions \bigl{\{}\xi_{\,j_{\,1}}^{\,n+1}\bigr{\}}_{\,j_{\,1}}\,, \bigl{\{}\eta_{\,j_{\,1}}^{\,n+1}\bigr{\}}_{\,j_{\,1}} of the free surface by using the finite difference approximation of the kinematic boundary conditions (3.5) and (3.12). 4. (4)

The new position $x\ =\ s^{\,n+1}$ of the moving wall is found from the discretized version of the Equation (LABEL:eq:sping). 5. (5)

Finally, the new grid $\bar{\Omega}_{\,h}^{\,n+1}$ is constructed on the following time layer $t\ =\ t^{\,n+1}\,$.

Every step above involves grid functions on various layers and to achieve a better understanding of all the stages of our algorithm, below we provide a detailed description of every item.

Approximation of the dynamic boundary condition (1)

The finite difference approximation to the dynamic boundary condition (3.6) takes the following form:

[TABLE]

where \boldsymbol{j}\ =\ \bigl{(}j_{\,1},\,j_{\,2}\bigr{)} is a multi-index with $j_{\,1}\ =\ 1,\,2,\,\ldots,\,N_{1}\,-\,1$ and $j_{\,2}\ =\ N_{2}\,$. In order to compute the speeds of grid nodes we employ finite differences applied to two last grids $\bar{\Omega}_{\,h}^{\,n}$ and $\bar{\Omega}_{\,h}^{\,n-1}\,$:

[TABLE]

where $\tau_{\,n-1}\ \equiv\ t^{\,n}\ -\ t^{\,n-1}$ is the local time step, \lvert\,\boldsymbol{u}_{\,\boldsymbol{j}}^{\,n}\,\rvert^{\,2}\ \equiv\ \bigl{(}u_{\,\boldsymbol{j}}^{\,n}\bigr{)}^{\,2}\ +\ \bigl{(}v_{\,\boldsymbol{j}}^{\,n}\bigr{)}^{\,2}\,. In formula (5.2) we use only the first order approximation in time. This is done for the sake of algorithm memory efficiency. Otherwise, to achieve at least the second order accuracy in $\tau\,$, we would have to keep in memory the grids for three time layers. In the current implementation we keep track of the preceding grid only.

The Cartesian components $u_{\,\boldsymbol{j}}^{\,n}$ and $v_{\,\boldsymbol{j}}^{\,n}$ of the velocity vector $\boldsymbol{u}_{\,\boldsymbol{j}}^{\,n}$ are computed after approximating formulas (3.15) and (3.16):

[TABLE]

Above one has to approximate also the partial derivatives $\varphi_{\,q^{\,1}}\,$, $\boldsymbol{x}_{\,q^{\,1}}$ over the independent variable $q^{\,1}$ along the upper side of the square $\mathcal{Q}^{\,0}$. These derivatives are computed with standard central finite differences up to the second order accuracy. The derivatives $\varphi_{\,q^{\,2}}\,$, $\boldsymbol{x}_{\,q^{\,2}}$ are computed with one-sided finite differences of the second order as well (in order to have the uniform second order accuracy in space).

We use equation (5.1) in order to compute the updated values of the velocity potential only in interior nodes of the grid $\gamma_{\,s}^{\,h}\,$. For two nodes located in the upper corners $\boldsymbol{q}_{\,0,\,N_{2}}\ =\ (0,\,1)$ and $\boldsymbol{q}_{\,N_{1},\,N_{2}}\ =\ (1,\,1)$ we use modified approximation formulas which take into account the lateral walls $\gamma_{\,\ell,\,r}$ impermeability. In the computational domain the left wall has the fixed coordinate $q^{\,1}\ =\ 0\,$. The intersection point of the free surface with the left boundary is a triple point and it has the coordinate $\boldsymbol{x}_{\,0,\,N_{2}}\,$. This point is permanently moving up and down the left wall.

In the coordinate system $O\,q^{\,1}\,q^{\,2}$ the Cartesian components $\boldsymbol{u}\ =\ \dfrac{\mathrm{d}\/\boldsymbol{x}}{\mathrm{d}\/t}$ become the so-called contravariant components of the velocity $\boldsymbol{v}\ =\ \dfrac{\mathrm{d}\/\boldsymbol{q}}{\mathrm{d}\/t}$ with

[TABLE]

Consequently, if a fluid particle moves such that its coordinate $q^{\,\alpha}$ does not change, then the corresponding contravariant velocity component $v^{\,\alpha}$ has to vanish. For instance, the fluid particles, which constitute the free surface, have the coordinate $q^{\,2}\ =\ 1\,$, consequently

[TABLE]

Similarly, for fluid particles sliding along the left vertical wall we always have $q^{\,1}\ =\ 0\,$, which yields

[TABLE]

The point $\boldsymbol{x}_{\,0,\,N_{2}}\ \in\ \gamma_{\,s}^{\,h}\ \cap\ \gamma_{\,\ell}^{\,h}\,$, thus

[TABLE]

From formula (5.4) and using also equations (3.18), (3.19) we obtain the following expressions for contravariant components of the velocity:

[TABLE]

From these expressions and using boundary condition777We implicitly use also the fact that the transformation (3.1) is non-degenerate. (5.5) for the triple point $\boldsymbol{x}_{\,0,\,N_{2}}$ we obtain that

[TABLE]

Henceforth, $\lvert\,\boldsymbol{u}\,\rvert^{\,2}\ \equiv\ u^{\,2}\ +\ v^{\,2}\ =\ u\,x_{\,t}\ +\ v\,y_{\,t}$ and the free surface dynamic boundary condition (3.6) in the corner $\boldsymbol{q}_{\,0,\,N_{2}}$ becomes:

[TABLE]

The difference equation (5.1) correspondingly takes the form:

[TABLE]

In a similar way one can derive the finite difference discretization of the boundary condition in the right triple point. We give only the final result. It coincides with equation (5.6) provided that we make a substitution $j_{\,1}\ =\ N_{1}\,$.

Computation of the velocity potential (2)

On the second stage of the numerical algorithm we solve the finite difference problem (4.6) in order to find the values of the velocity potential \bigl{\{}\varphi_{\,\boldsymbol{j}}^{\,n+1}\bigr{\}}_{\,\boldsymbol{j}} in the nodes $\boldsymbol{q}_{\,\boldsymbol{j}}\ \in\ \bar{\mathcal{Q}}_{\,h}^{\,0}\ \setminus\ \gamma_{\,s}^{\,h}$. In other words, in the nodes \boldsymbol{j}\ =\ \bigl{(}j_{\,1},\,j_{\,2}) with $0\ \leqslant\ j_{\,1}\ \leqslant\ N_{1}$ and $0\ \leqslant\ j_{\,2}\ <\ N_{2}\,$.

In order to invert the linear system (4.6) we tested several iterative methods and no method clearly outperformed the others. In the final implementation of the code we used the Successive Over Relaxation (SOR) method [78] that we remind briefly here. The choice of the SOR method can be explained by two main reasons:

1. (1)

We would like to have a sufficiently simple to implement and sufficiently efficient method to find approximate solutions to relatively large linear systems 2. (2)

The method does not have to rely heavily on the matrix structure (such as Thomas’s algorithm, for example). For instance, if we include obstacles in the fluid domain, its topology might change, which implies some drastic changes in the pattern of non-zero elements of the matrix.

Thus, to our opinion SOR method represents the best trade-off among the efficiency, generality and ease of implementation.

We use a nine-point stencil in the finite difference scheme. Thus, equation $\boldsymbol{j}$ has the following general form:

[TABLE]

The last condensed form is obtained from (4.6) by regrouping all the terms in front of unknowns \bigl{\{}\varphi_{\,\boldsymbol{j}}^{\,n+1}\bigr{\}}_{\,\boldsymbol{j}}\,. The corresponding coefficient is denoted by $\mathscr{A}_{\,k,\,\boldsymbol{j}}\,$. In equation (5.7) we used a local numeration of indices depicted in Figure 4(a). The right hand side is given by

[TABLE]

Then, one step of the SOR method reads:

[TABLE]

[TABLE]

where $(m)$ is the iteration number and \texttheta is the relaxation parameter. The value of this parameter weakly influences the iterative process speed of convergence. The question of finding the optimal value of the relaxation parameter \texttheta is the key to the efficiency of the SOR method. This question can be studied theoretically for the Dirichlet problem of the Poisson equation in a square domain $[\,0,\,\ell\,]^{\,2}$ with uniform grids. We know also that the optimal value $\mbox{\texttheta}^{\,\star}$ always belongs to the interval:

[TABLE]

For the Poisson equation on a square domain with the uniform discretization $h_{\,1}\ =\ h_{\,2}\ \equiv\ h\,$, one can show

[TABLE]

In the last formula for $\ell\ =\ 1$ and $h\ =\ \frac{1}{30}$ we obtain the value $\mbox{\texttheta}^{\,\star}\ \approx\ 1.81\,$. In our discrete problem for the velocity potential $\varphi_{\,\boldsymbol{j}}$ we cannot determine theoretically the optimal value $\mbox{\texttheta}^{\,\star}\,$. This value was determined experimentally for each problem under consideration. In numerical computations below we always took $\mbox{\texttheta}^{\,\star}\ \in\ [\,1.85,\,1.95\,]\,$, depending on the discretization parameters $h_{\,1,\,2}$ as well.

In boundary nodes the stencil contains six points and in corner points only four. However, it does not change anything for the SOR scheme since extra coefficients $\mathscr{A}_{\,k,\,\boldsymbol{j}}$ can be set to zero for our convenience. The expressions of coefficients $\mathscr{A}_{\,k,\,\boldsymbol{j}}$ are given in Table 1. In this Table we use the following notation:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

In the definition of quantities $\mbox{\textrho}_{\,k}$ we used the notation introduced in equations (4.11) – (4.14). Finally, the coefficient $\mathscr{A}_{\,0,\,\boldsymbol{j}}$ is determined as

[TABLE]

This concludes the description of our velocity potential solver.

Free surface motion (3)

The free surface position on the following time layer $t\ =\ t^{\,n+1}$ is found by integrating in time kinematic conditions (2.7), (2.8). Here we consider the simplest approximation of these equations using an explicit upwind scheme:

[TABLE]

where we introduced the vector \mbox{{\texteta}}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \bigl{(}\eta,\,\xi\bigr{)} and index $j_{\,1}\ =\ 1,\,2,\,\ldots,\,N_{1}\,-\,1\,$. The contravariant velocity component $v_{\,\boldsymbol{j}}^{\,1}$ is approximated as

[TABLE]

The components of the (Cartesian) velocity vector $\boldsymbol{u}_{\,\boldsymbol{j}}^{\,n}$ are computed using formulas (5.3) except the fact that we use the values of the velocity potential \bigl{\{}\varphi_{\,\boldsymbol{j}}^{\,n+1}\bigr{\}}_{\,\boldsymbol{j}} on the subsequent time layer $t\ =\ t^{\,n+1}\,$.

We employ formula (5.9) in order to compute the free surface position in ‘interior’ nodes $\boldsymbol{q}_{\,j_{\,1},\,N_{1}}$ with $j_{\,1}\ =\ 1,\,2,\,\ldots,\,N_{1}\,-\,1\,$, which lie on the upper side of the square $\mathcal{Q}_{\,h}^{\,0}\,$. Thanks to impermeability conditions (5.5), kinematic conditions (2.7), (2.8) take a very simple form:

[TABLE]

Taking into account this information, the finite difference formula (5.9) takes a much simpler form as well:

[TABLE]

The fully discrete scheme to compute the free surface transport is described.

Moving wall displacement (4)

The movable wall position is given by the function $x\ =\ s\,(t)\,$, which is a solution to equation (2.10). The fluid pressure $p$ is determined from the Cauchy–Lagrange integral given in dimensionless variables in equation (2.20). The force acting on the vertical wall is given in dimensionless variables in equation (2.19). In curvilinear coordinates the pressure is given in equation (3.11) and it can be used to compute the force in curvilinear coordinates as well (3.10). The advantage of this formulation is that the wall has a fixed position $q^{\,1}\ =\ 0$ in the transformed domain. The integral in equation (3.10) is discretized using the trapezoidal quadrature formula:

[TABLE]

In order to compute the time derivative of the velocity potential $\varphi_{\,t}$ appearing in the pressure term \bigl{\{}p_{\,0,\,j_{\,2}}\bigr{\}}_{j_{\,2}\,=\,0}^{\,N_{\,2}}\,, we use the values of the velocity potential on three time layers \bigl{\{}\varphi_{\,\boldsymbol{j}}^{\,n-1}\bigr{\}}_{\,\boldsymbol{j}}\,, \bigl{\{}\varphi_{\,\boldsymbol{j}}^{\,n}\bigr{\}}_{\,\boldsymbol{j}} and \bigl{\{}\varphi_{\,\boldsymbol{j}}^{\,n+1}\bigr{\}}_{\,\boldsymbol{j}}\,. So, the time derivative of the velocity potential is approximated as:

[TABLE]

However, due to the change of variables in (3.10) one has to compute also the velocities of grid nodes. This is done as specified in (5.2). Cartesian velocities computed as specified in (5.3) with the only difference — we use the values of velocity potential at the new time level \bigl{\{}\varphi_{\,\boldsymbol{j}}^{\,n+1}\bigr{\}}_{\,\boldsymbol{j}}\,.

In order to integrate numerically the nonlinear differential equation (2.18), we rewrite this equation as a system of two first order differential equations:

[TABLE]

together with initial conditions:

[TABLE]

The last system of equations can be rewritten in the vectorial form for our convenience:

[TABLE]

where we introduced the following notations:

[TABLE]

The Cauchy problem (5.11) is solved using the so-called modified Euler method888This scheme is also called the second order explicit Runge–Kutta scheme [4].:

[TABLE]

If we exclude the intermediate variable $\upupsilon^{\,\star}$ from last equations and rewrite the obtained equations in the component-wise form, we obtain the following scheme:

[TABLE]

The last scheme is implemented in our numerical code.

Elliptic mapping construction (5)

In previous sections we assumed that mapping (3.1) was known. In order to complete our numerical method description we have to describe how to construct this mapping in practice. Namely, we shall use the so-called equidistribution method [45]. A typical grid generated using this method is shown in Figure 5. Recently this method was successfully applied in 1D to the simulation of conservation laws [44]. Below we provide a detailed description of this method to two spatial dimensions.

According to the equidistribution method, the coordinates of grid nodes \bigl{\{}\boldsymbol{x}_{\,\boldsymbol{j}}^{\,n+1}\bigr{\}}_{\,\boldsymbol{j}}\ \subseteq\ \Omega\,(t) in the physical space on the subsequent time layer $t\ =\ t^{\,n+1}$ are determined by solving the following vectorial parabolic equation:

[TABLE]

where $\mbox{\textvarsigma}\ >\ 0$ is a positive smoothing parameter whose value is chosen to minimize the oscillations in grid nodes trajectories. Above $\varpi\ =\ \varpi\,(\boldsymbol{x},\,t)$ is the so-called monitoring function, which determines the local density of the grid nodes and \bigl{\{}g_{\,\alpha\,\alpha}\bigr{\}}_{\,\alpha\,=\,1}^{\,2} are metric tensor components defined in (3.14).

The equidistribution method based on the solution of parabolic equations for grid nodes coordinates was proposed in [71]. A few years later this method was applied to the simulation of 1D (compressible) gas dynamics problems in [1, 72]. This method has been especially designed for non-stationary problems and its application allows to avoid (or simply reduce) abrupt changes in the nodes positions when we move from one time layer to the following. In this method one constructs the new mesh by taking into account the position of nodes at the last time layer. The new coordinates of grid nodes are determined by solving a discretized version of the parabolic problem (5.12). This modification of the classical equidistribution method (proposed in [17]) of the construction of dynamically adapted moving grids has been used to solve numerous 1D problems [14, 13, 15, 18, 22, 80, 58], this list is not being exhaustive.

The system of nonlinear equations (5.12) is solved using an iterative method of alternating directions. The monitoring function $\varpi$ was computed at the known solution on the time layer $t\ =\ t^{\,n}\,$. As the initial guess for this iterative process we take the known mesh \bigl{\{}\boldsymbol{x}_{\,\boldsymbol{j}}^{\,n}\bigr{\}}_{\,\boldsymbol{j}} from the previous time layer as well. We mention that at the initial moment of time $t\ =\ 0$ we employ the traditional (steady) equidistribution method as it was described in [46]. It is equivalent to solve the system (5.12) with $\mbox{\textvarsigma}\ \equiv\ 0\,$.

Since the position of boundary nodes varies999Indeed, we are dealing with free surface flows in the presence of a moving wall. from one time step to another, before solving equations (5.12) we have to determine the grid \bigl{\{}\boldsymbol{x}_{\,\boldsymbol{j}}^{\,n+1}\bigr{\}}_{\,\boldsymbol{j}} along the boundary $\partial\,\mathcal{Q}_{\,h}^{\,0}\,$. Let us assume that a portion \mathring{\Gamma}\ \in\ \bigl{\{}\Gamma_{\,\ell},\,\Gamma_{\,r},\,\Gamma_{\,b},\,\Gamma_{\,s}\bigr{\}}\ \equiv\ \partial\,\mathcal{Q}^{\,0} of the boundary $\Gamma^{\,n+1}\ \equiv\ \Gamma(t^{\,n+1})$ is given in the following parametric form:

[TABLE]

where \textzeta is a real parameter. Then, the coordinates of grid nodes \bigl{\{}\boldsymbol{x}_{\,\boldsymbol{j}}^{\,n+1}\bigr{\}}_{\,\boldsymbol{j}}\ \subseteq\ \mathring{\Gamma}^{\,h} on the boundary can be determined by formulas:

[TABLE]

We assume that the nodes on the boundary under consideration are re-ordered with a single scalar index $j$ for the sake of notation compactness. So, the problem is to determine the discrete scalar function \bigl{\{}\mbox{\textzeta}_{\,j}^{\,n+1}\bigr{\}}_{\,j} in the parameter space. To do it, we solve numerically the following scalar nonlinear difference equation:

[TABLE]

where h\ \in\ \bigl{\{}h_{\,1},\,h_{\,2}\bigr{\}} is the discrete step in space, which depends on the boundary under consideration. The function \tilde{\mathcal{J}}_{\,j}^{\,n+1}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \Bigl{\{}\sqrt{x_{\,\zeta}^{\,2}\ +\ y_{\,\zeta}^{\,2}}\Bigr{\}}\Bigr{|}_{\,j}^{\,n+1} is the trace of the transformation Jacobian at this boundary. Please, notice that we compute the monitoring function $\varpi^{\,n}$ on the known position of the nodes. The difference equation (5.13) is solved on four sides (left and right walls, bottom and free surface) of the discretized rectangle $\mathcal{Q}_{\,h}^{\,0}\,$. The distribution of boundary nodes obtained in this way is used as the Dirichlet-type boundary condition for the 2D grid generation, which is achieved by solving numerically equation (5.12). By the end of this step we know the grid $\bar{\Omega}_{\,h}^{\,n+1}$ nodes coordinates \bigl{\{}\boldsymbol{x}_{\,\boldsymbol{j}}^{\,n+1}\bigr{\}}_{\,\boldsymbol{j}}\,.

5.5.1  The re-computation

At the next step we recompute quantities $\mbox{{\texteta}}^{\,n+1}$ and \bigl{\{}\varphi^{\,n+1}_{\,\boldsymbol{j}}\bigr{\}}_{\,\boldsymbol{j}} to take into account the new mesh. For this, we solve again equation (5.1) (or an analogue of (5.6) on the boundaries) in which the Cartesian components of the velocity (5.3) are still computed on the previous grid $\bar{\Omega}_{\,h}^{\,n}$ and using the velocity potential values from the previous time layer $t\ =\ t^{\,n}\,$. However, in formulas (5.3) we take the new velocity of grid nodes. In other words, formula (5.2) is replaced by

[TABLE]

During the second solution of system (5.7) the right hand side member (5.8) is modified as well to take into account the new velocity of the left wall and, obviously, the new location of mesh nodes on this boundary. Another modification concerns the velocity components $u\,$, $v$ when they are used in the computation of the contravariant component $v_{\,\boldsymbol{j}}^{\,1,\,n}$ whose discretized formula is given in (5.10). The quantity $v_{\,\boldsymbol{j}}^{\,1,\,n}$ is needed to determine the new position of the free surface elevation $\mbox{{\texteta}}^{\,n+1}\,$. During the present step, the components $u\,$, $v$ are computed as

[TABLE]

The last formulas replace equations (5.3) during this stage.

Stability of the scheme

In order to study the stability of the proposed scheme, we consider the linearized governing equations. Moreover, we consider only the Initial Value Problem101010In fact, we have the so-called Cauchy–Poisson problem, which is of mixed type between IVP and BVP. We say that we have an IVP since there are no boundary conditions on lateral boundaries, i.e. the domain is unbounded in the horizontal directions. However, there are still boundary conditions to be satisfied on the impermeable bottom and on the (linearized) free surface. (IVP) in order to remove the complexity of (lateral) boundary conditions treatment. The linear water wave problem is known as the Cauchy–Poisson problem [73, 59, 27, 42] since the pioneering works of Augustin Louis Cauchy [6] and Siméon Denis Poisson [66]. Let us formulate this problem in precise mathematical terms. Consider a fluid layer of infinite horizontal extent ($-\infty\ <\ x\ <\ +\infty$) over a solid bottom of uniform depth $y\ =\ -h_{\,0}\ =\ \mathrm{const}\,$. The fluid domain is bounded from above by the free surface, whose location is assumed to be at $y\ =\ 0$ after the linearization. The Cauchy–Poisson problem consists in finding the free surface elevation $\eta$ and the velocity potential $\varphi$ by solving the Laplace equation111111This equation is to be compared with equation (2.12). in a strip:

[TABLE]

The last equation has to be completed by the following free surface boundary conditions:

[TABLE]

and by one bottom impermeability condition:

[TABLE]

It is not difficult to see that free surface boundary conditions given above are nothing else but linearized versions of equations (2.13) and (2.14). A particular solution to the Cauchy–Poisson problem can be easily obtained using some elementary Fourier analysis [73, 59, 27, 42]:

[TABLE]

where $k\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \dfrac{2\,\pi}{\lambda}$ is the wave number, $\lambda$ is the wave length, $\omega\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \dfrac{2\,\pi}{T}$ is the wave frequency and $T$ is its period. Wave amplitudes $\eta_{\,0}\ \in\ \mathds{R}$ and $\varphi_{\,0}\ \in\ \mathds{R}$ are some real numbers. The wave frequency $\omega$ is related to the wave number $k$ through the so-called dispersion relation of gravity waves:

[TABLE]

For the sake of simplicity we consider only waves moving rightwards. It fixes the branch $+$ in the relation above. We reiterate the fact that the dispersion relation (5.17) is a necessary condition for the existence of solutions (5.16).

The numerical scheme considered above can be applied to the Cauchy–Poisson problem as well. The semi-discretization in time of free surface boundary conditions (5.14), (5.15) reads:

[TABLE]

From now on, for the sake of simplicity we take the time step $\tau_{\,n}\ \equiv\ \tau\ >\ 0$ to be constant. The elementary solution (5.16) can be semi-discretized as well:

[TABLE]

where we introduced the notation $\uprho\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \mathrm{e}^{-\,\mathrm{i}\,\omega\,\tau}\,$. By substituting this semi-discrete solution ansatz into relations (5.18), (5.19) we obtain the following linear system of equations with respect to wave amplitudes $\eta_{\,0}$ and $\varphi_{\,0}\,$:

[TABLE]

In order to have non-trivial solutions, the determinant of this system has to vanish. It gives us the following quadratic equation with respect to the transfer coefficient $\uprho\,$:

[TABLE]

To have the linear stability property of our scheme, it is necessary that both roots of this equation verify the inequality $\lvert\,\rho^{\,\pm}\,\rvert\ \leqslant\ 1\,$. To meet this requirement, it is sufficient to ask that the discriminant of the quadratic equation (5.20) is not positive, i.e.

[TABLE]

or equivalently

[TABLE]

where the wave frequency was defined in equation (5.17) (with the sign $+$ by our convention).

Let $\Delta\,x$ be the discretization step size in the horizontal direction. The minimal wave length $\lambda_{\,\inf}\,$, which can be represented on the grid with spacing $\Delta\,x$ is $\lambda_{\,\inf}\ \equiv\ 2\,\Delta\,x\,$. All other waves satisfy the inequality $\lambda\ \geqslant\ \lambda_{\,\inf}\,$. These considerations on the wave length can be translated into the language of wave numbers, i.e.

[TABLE]

Henceforth, we can derive the following estimation for the wave frequency using the dispersion relation (5.17):

[TABLE]

Consequently, in order to satisfy the stability condition (5.21), it is sufficient to impose the following restriction on the time step $\tau$ magnitude:

[TABLE]

An analogue of this stability condition will be used below during the simulation of fully nonlinear problems. The most important conclusion of this Section is that we have a hyperbolic-type stability Courant–Friedrichs–Lewy condition [11] — the time step is a linear function of the mesh spacing.

5.6.1  Practical choice of the time step

In the previous Section we studied the scheme stability in the linear case. However, in numerical simulations presented below, we solve the nonlinear problem. Henceforth, one may ask the question how to choose the time step in practical nonlinear simulations. Below we explain our approach to this problem.

In order to cover the whole family of problems, we choose to work in scaled variables. The CFL condition (5.22) can be rewritten in dimensionless variables as

[TABLE]

where $h_{\,1}^{\,\ast}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \dfrac{h_{\,1}}{h_{\,0}}\,$. As grid spacing $h_{\,1}\,$, we take the smallest horizontal spacing along the free surface (since most important variations take place there):

[TABLE]

Consequently, the dimensionless CFL condition can be rewritten as

[TABLE]

However, our numerical simulations show that this estimation is too pessimistic. Consequently, in all simulations presented below we took the time step according to the following less restrictive formula:

[TABLE]

with $\kappa\ =\ 0.95\,$. Most probably, the last condition is pessimistic as well. However, it guaranteed the stability of our nonlinear computations.

Numerical results

Above we described the proposed finite difference scheme and our resolution algorithm on a fixed reference domain. Below we present several validation tests and numerical experiments which show the performance and abilities of our numerical approach.

Solitary wave run-up on a fixed wall

In order to illustrate the applicability of our numerical algorithm, we consider the classical problem of the solitary wave/fixed wall interaction. Due to symmetry considerations, this set-up is equivalent to the head-on collision of two equal solitary waves. This problem is well-studied in the literature [57, 60, 74, 10, 7] and it can serve as the first validation test.

Consider a solitary wave of amplitude $\upalpha$ moving in the leftward direction. The channel has a constant depth and the wall is assumed to be fixed in order to be able to perform comparisons with previous investigations. On the right the channel is also bounded by a fixed vertical wall. The total length of the channel is equal to $\ell\ =\ 20\,$. In this Section we provide all values in dimensionless variables as it was explained earlier in Section 2.3. Moreover, we assume in these computations that waves do not overturn. In other words, the free surface is traditionally given as the graph of a function $y\ =\ \eta\,(x,\,t)\,$. The initial condition for the free surface elevation $\eta_{\,0}\,(x)\ \equiv\ \eta\,(x,\,0)$ and velocity field $\boldsymbol{u}_{\,0}\,(\boldsymbol{x})\ \equiv\ \boldsymbol{u}\,(\boldsymbol{x},\,0)$ are given by the following approximate formulas [65]:

[TABLE]

[TABLE]

[TABLE]

where $x_{\,0}$ is the wave crest initial position.

The interaction of a solitary wave with amplitude $\upalpha\ =\ 0.4$ with a vertical (fixed) wall is shown in Figure 6 at different moments of time. Under a wave we show also the adapted grid. On this Figure we can see that after the reflection the solitary wave does not recover completely its initial shape. In particular, behind the main wave we observe a slight dispersive tail. The numerical simulations on refined grids show that it is not a numerical effect. This dispersive tail appears on all grids and it reflects an intrinsic property of the full Euler equations — their non-integrability [30, 28]. In this example we use the following monitoring function in order to adapt the grid to the solution:

[TABLE]

where $\upvartheta\ >\ 0$ is a positive ad-hoc parameter. In computations presented in the Section we use $\upvartheta\ =\ 10\,$. The main rationale behind this choice of the monitoring function $\varpi\,(\boldsymbol{x},\,t)$ is to put more nodes in areas where the waves are large.

In Figure 6 we can see that the refined area follows somehow the solitary wave crest during its motion towards and fromwards the vertical wall. From equation (6.4) it can be seen that the monitoring function does not depend on the vertical variable $y\,$. As a result, we obtain the grids with almost vertical lines. Consequently, in accordance with stability condition (5.22) to determine an admissible local time step $\tau_{\,n}$ we can use the following formula:

[TABLE]

where $0\ <\ \upnu\ \leqslant\ 1$ is the security factor and $\Delta\,x_{\,\min}^{\,n}$ is the minimal mesh spacing on the free surface at time $t\ =\ t^{\,n}\,$, i.e.

[TABLE]

One of the main characteristics of the wave/wall interaction is the maximal amplitude of the wave on the wall. This quantity is called the maximal run-up and will be denoted in our study as $\mathcal{R}_{\,\max}\,$. Obviously, in our experimental conditions this quantity depends on the amplitude of the incident solitary wave, i.e. $\mathcal{R}_{\,\max}\ =\ \mathcal{R}_{\,\max}\,(\upalpha)\,$. In Figure 7 we represent this dependence according to our numerical simulations (solid black line), experimental data (filled markers), other computations (empty markers) and the following asymptotic analytical prediction [74] (dashed line):

[TABLE]

In Figure 7 we can see that there is an overall good agreement among all presented data up to the amplitude $\upalpha\ \lesssim\ 0.4\,$. For higher waves some divergences start to appear. However, the agreement between our numerical model with other potential flow solvers [34, 10] continues up to $\upalpha\ \lesssim\ 0.6\,$. The experimental points go rather below our predictions. It can be easily explained by the neglection of viscous and friction effects in our numerical (and mathematical) model. Nevertheless, we would like to mention that our numerical results agree particularly well with experimental data reported in [16].

Wave generation by a numerical wave maker

The study of water wave interaction with movable partially submerged bodies and objects is a very important problem of the modern computational hydrodynamics. In the previous Section 6.1 the wall was considered to be fixed in order to compare our numerical predictions with available data. Starting from this Section we allow the left wall as a movable object. As the first step towards the freely moving solid boundary, we consider first the situation where left wall motion is prescribed by a given law. Physically, it corresponds to the vertical piston motion used in many experimental facilities. Traditionally, the generated waves are understood using linear or weakly nonlinear theories [55]. Recently this situation was modelled numerically with Boussinesq-type equations in [63]. The results presented below are fully nonlinear and fully dispersive.

Consider a numerical wave tank with horizontal bottom and of uniform depth $d\ =\ 1\ \mathsf{m}$ (when the water is unperturbed). The left wall is initially located in the point $x\ =\ 0\,$. For times $t\ >\ 0$ the wall moves according to the following law:

[TABLE]

Thus, during the initial times the wall moves to the right. In this way, the wall motion is completely determined by three parameters $\alpha\ >\ 0\,$, $\beta\ >\ 0$ and $\omega\,$. The first one ($\alpha$) specifies the maximal wall oscillation amplitude (in the horizontal extent), the second parameter $\beta$ controls the speed of the relaxation121212It is easy to see that the wall oscillation amplitude tends to its maximum value $\alpha$ as $t\ \to\ +\infty\,$. This justifies the term ‘relaxation’ associated to this parameter $\beta\,$. The speed of this convergence depends on the value of parameter $\beta\,$, i.e. bigger is faster. towards the stationary periodic regime and $\omega$ is the frequency of wall oscillations (or equivalently it controls the oscillation period).

In Figure 8 we depict a few simulations results, which were obtained by varying the parameter $\omega$ for fixed values of two other parameters:

[TABLE]

Figure 8 shows clearly that the generated wave amplitude (measured at the moving wall) depends essentially on the wall oscillation period, which is determined by parameter $\omega\,$. This effect is illustrated in Figure 8 for three different values of the parameter \omega\ \in\ \bigl{\{}0.5\;\mathsf{Hz},\,1.0\;\mathsf{Hz},\,3.0\;\mathsf{Hz}\bigr{\}}\,. For each wall trajectory we show also the corresponding free surface excursion on the moving vertical wall. The increase in $\omega$ results in the reduction of the oscillation period. We can also witness the relaxation effect on the free surface oscillation. Indeed, as the time $t$ goes on, the wave amplitude increases and becomes practically stationary. The speed of relaxation depends on the parameter $\beta\,$. We underline the fact that for all three (considered) values of the wall oscillation frequency $\omega\,$, the amplitude of oscillations $\alpha$ was the same ($\alpha\ \equiv\ 0.4\;\mathsf{m}$). Nevertheless, the wave amplitudes registered on the wall were significantly different. This is the first practical conclusion that we can draw from our numerical simulations: the piston oscillation period has a much bigger influence on the generated wave amplitude, than the amplitude of these oscillations. We can explain this observation with a simple physical argument: during ‘fast’ piston oscillations, the energy is pumped into the fluid close to the piston faster than waves can evacuate this energy by propagating into the wave tank.

The second important observation that we can make based on our simulations is that the maximal wave run-up does not take place at the moments where the wall displacement is maximal. A more detailed investigation showed that the wave amplitude on the wall is maximal when the wall acceleration is maximal. It happens when the wall passes by its initial (or mean, or unperturbed) position $x\ =\ 0\,$. After passing this point, the wall decelerates and the wave has time to flow away from the piston, which reduces the value of run-up.

Finally, we can also observe that the wave run-up value on the wall is asymmetric with respect to the still (unperturbed) water level. This effect is much better visible for fast wall motions. Indeed, for $\omega\ =\ 3.0\;\mathsf{Hz}$ the maximal run-up is $\ \approx\ 43\ \mathsf{cm}\,$, while the maximal run-down is $\ \approx\ 34\ \mathsf{cm}\,$.

Surface wave run-up on a movable wall

In this Section we consider the wave/wall interaction problem, where the wall motion is not prescribed, but it is determined as a part of the problem solution. The left wall is attached to a spring system, which can deform under the wave action.

6.3.1  Single pulse

To begin, we consider first the case of a single localized wave impulse interacting with a moving wall. The channel depth is taken to be $d\ =\ 1\ \mathsf{m}$ and the channel length is $\ell\ =\ 20\ \mathsf{m}\,$. At the initial moment of time $t\ =\ 0\,$, the velocity field is taken to be quiescent and the free surface shape is given by the following formula:

[TABLE]

For subsequent times $t\ >\ 0$ this initial condition is separated in two waves of the amplitude $\ \approx\ \frac{\upalpha}{2}\,$, which move in opposite directions. The wave moving leftwards interacts with the moving wall. We would like to underline the fact that the initial left wall position ($x\ =\ 0$) does not coincide with the wall equilibrium position in the absence of water. Even if the free surface coincides with the still water level, the springs are deformed under the action of the hydrostatic pressure. Thus, the equilibrium wall position $x\ =\ x_{\,e}$ with undeformed springs is located somewhere to the right from its initial position, i.e. $x_{\,e}\ >\ 0\,$. We performed a series of numerical simulations by varying the parameters $\upalpha\,$, $x_{\,0}\,$ and $\Upsilon$ of the initial condition (6.6). Additionally we varied also the wall mass $m$ and the springs rigidity coefficient $k\,$.

First, we conducted a series of numerical experiments for small initial wave amplitudes $\upalpha$ in view of comparisons against the numerical simulations reported in [40]. The Authors of the preceding study employed the Boundary Integral Equations Method (BIEM) [2, 3]. Unfortunately, the Authors of [40] did not publish any tabulated data to perform quantitative comparisons. However, we performed qualitative comparisons of numerical results. We choose the same values as in [40], i.e. $x_{\,0}\ =\ 0.7\ \mathsf{m}\,$, $\Upsilon\ =\ 0.5\ \mathsf{m}$ and dimensionless spring rigidity coefficient $k\ =\ 10\,$. The wall dimensionless mass is $m\ =\ 1.5\,$.

In Figure 9(a) we show the free surface oscillations recorded at the initial pulse location131313Here we mean the center or pulse crest. $x\ =\ x_{\,0}\ =\ 0.7\ \mathsf{m}\,$. We show the results for two initial wave amplitudes. The line (1) corresponds to $\upalpha\ =\ 1\ \mathsf{cm}\,$, while $\upalpha\ =\ 10\ \mathsf{cm}$ is depicted by line (2). Please, notice that curves presented in Figure 9 are shown relatively to the initial wave amplitude $\upalpha\,$. These curves compare quantitatively and qualitatively well to those presented in [40]. To the graphical resolution the oscillation decay rates and oscillation periods are the same. In Figure 9(b) we show moving wall trajectories in time for two initial wave amplitudes \upalpha\ \in\ \bigl{\{}1\;\mathsf{cm},\,10\;\mathsf{cm}\bigr{\}}\,. These curves also compare well with previous computations reported in [40]. We can also see that for small times the solid and dashed curves shown in Figure 9(a,b) almost coincide. It means that for small initial wave amplitudes $\upalpha$ the results of computations depend almost linearly on $\upalpha\,$. This observation supports the applicability of the linear theory [48].

Let us study now the influence of the springs rigidity parameter $k$ on the resulting wave field. We took the following parameters in this computation:

[TABLE]

The parameter $k$ is variable. In Figure 10(a) we can see that the first run-up on the wall is large. Then, the following run-ups are generally decaying in amplitude since the wave energy spreads over the whole wave tank during the dispersive (and possibly nonlinear) effects. Moreover, for $k\ =\ 5$ (curve 2) the decay rate is the strongest since a larger part of the wave energy is converted into the wall motion. Thus, we have the wave energy transformation into the elastic energy of the spring system. When the springs rigidity is high, the wall performs small later motions, as it can be seen in Figure 10(b). The increase in $k$ results also in the increase of the wall oscillation frequency around the initial equilibrium position $x\ =\ 0\,$. In the same time, the wave run-up records on the wall are practically identical with the fixed wall case (cf. curves 3 and 4 in Figure 10(a)). In other words, for high rigidities $k\ \gg\ 1\,$, the ‘moving’ wall interacts in a way very similar to the fixed wall. On the other hand, when the rigidity decreases, the wall oscillation period and amplitude both increase. Wave run-up on the wall decay rate also increases with the parameter $k\,$. However, in contrast to large values of $k\,$, a flexible spring system can allow the maximal wave run-up value during the second (or even further subsequent) wave interactions with the wall. In Figure 10(a) we can see that the value of the maximal wave run-up on the wall is generally lower for weak springs comparing to quasi-fixed wall results ($k\ \gg\ 1$). From these simulations we can conjecture that flexible spring systems should be used for more efficient wave run-up reduction.

Let us consider now the case of the weak springs system ($k\ =\ 1$) and a moving wall with mass $m\ =\ 1\,$. The question we would like to investigate is how the wave run-up and wall trajectory depend on the incident wave amplitude $\upalpha\,$? In Figure 11(a) we show the wave run-up value on the moving wall for a moderate amplitude $\upalpha\ =\ 20\ \mathsf{cm}$ (curve 1) and a fairly strong amplitude $\upalpha\ =\ 60\ \mathsf{cm}$ (curve 2). The initial condition is given by the single pulse formula (6.6). The initial impulse center is $x_{\,0}\ =\ 12\ \mathsf{m}$ and its half-length $\Upsilon\ =\ 5\ \mathsf{m}\,$. In Figure 11(b) we show the wall displacement $s\,(t)\,$. As expected, the wall displacement becomes larger as the incident wave amplitude $\upalpha$ increases. It can be seen that for weak springs the wall first retracts under the incident wave action. Then, it comes back towards its initial position and slightly oscillates around this point. It is very interesting to observe the high amplitude pulse (curve 2) interaction with the moving wall depicted in Figure 11(a). Here the wave run-up value grows until the wall attains its left-most position. Then, on the way back the wall motion provokes the second maximal run-up value. It explains the non-monotonic behaviour of the curve 2 around its maximum. This behaviour is possible only for high amplitude waves and weak (i.e. easily deformable) springs.

We naturally come to the question of the wall mass $m$ and springs rigidity $k$ on the maximal wave run-up value of the impulse (6.6) on the movable wall. We use the following parameters of the incident pulse:

[TABLE]

In Figure 12(a) we depict the maximal run-up value as the function of the wall mass $m\,$. It can be easily seen that for a fixed rigidity parameter $k\,$, the maximal run-up is a non-monotonic function of $m\,$. Even more, for every value of $k\,$, one can find the value $m^{\,\star}$ for which the maximal wave run-up will be minimal. It implies that the wave run-up can always be reduced by choosing the appropriate wall parameters $k^{\,\star}$ and $m^{\,\star}\,$. In Figure 12(b) we show the maximal wave run-up dependence on $k$ for fixed dimensionless wall masses m\ \in\ \bigl{\{}1,\,3,\,5\bigr{\}}\,. So, it can be seen that the run-up depends monotonically on $k\,$. Namely, it increases with wall rigidity $k\,$. As we increase $k\ \to\ +\infty\,$, the maximal run-up tends to its limiting value, which depends on the incident wave parameters and which is independent on the dimensionless wall mass $m\,$.

6.3.2  Solitary wave

Finally, we show some numerical results on the solitary wave interaction with a moving wall. This test case is a natural extension of the case considered earlier in Section 6.1. The initial condition is provided by the same formulas (6.1) – (6.3). In order to avoid any interactions of reflected waves from the right wall (during the simulation time), we increase the computational domain length to $\ell\ =\ 40\ \mathsf{m}$ (from $\ell\ =\ 20\ \mathsf{m}$). The computations were performed for a solitary wave with moderate amplitude $\upalpha\ =\ 20\ \mathsf{cm}\,$. We would like to study the influence of moving wall parameters on the maximal solitary wave run-up. In Figure 13 we show the maximal run-up as a function of the wall mass $m$ (a) and as a function of springs rigidity $k$ (b) (the other parameter being fixed). With a dashed line we show the maximal run-up value on a fixed wall. The comparison of this result with Figure 12 indicates that there are qualitative similarities in the behaviour of a single pulse with a solitary wave (cf. Figure 13). So, for solitary waves the run-up can be significantly reduced as well. All differences observed in Figures 12 and 13 can be solely explained by the incident wave amplitude and its shape.

In Figure 14(a) we show the wave run-up record under the action of a solitary wave on the moving wall as a function of time. This computation was performed for the wall mass $m\ =\ 5$ and springs rigidity is $k\ =\ 1\,$. In Figure 14(b) we report the wall displacement during the wave/wall interaction process. In Figure 14(a) we depict also the wave run-up record for a fixed wall as well. It can be clearly seen that the presence of a moving wall reduces significantly wave run-up height on it. We can also notice that for a movable wall the wave run-up happens twice in contrast to the fixed wall case. Indeed, the first (and usually maximal) run-up takes place when the wall is retracting under the incident wave action. Then, at some moment the springs accumulate enough elastic energy in order to push back the wall towards its equilibrium position. Right after this turning point the second run-up takes place. However, its value is usually lower than the first one.

In Figure 15 we show also the free surface evolution in space and in time. The left panel 15(a) shows the classical fixed wall case for the sake of comparison. In this case the wave field essentially consists of one incident and one reflected waves. Eventual inelastic collision effects are negligible for a solitary wave of amplitude $\upalpha\ =\ 20\ \mathsf{cm}\,$. A much more interesting wave field is generated by the interaction with a moving wall. It is depicted in Figure 15(b). The wall motion generates multiple reflected wave, which travel with lower speed (in agreement with their lower amplitude).

Discussion

The present study was devoted to the simulation of free surface waves in a two-dimensional tank of variable size. Namely, a vertical wall is attached to a system of springs and it can move under the action of incident water waves. Modelling and numerical simulation of this problem was discussed in this study. Below we outline the main conclusions and perspectives of the present study.

Conclusions

In the present study we considered the problem of surface water wave modelling in domains with variable geometry. In general, problems on time-varying are known to be notoriously difficult [47] (both theoretically and numerically). Here we considered a special case of a numerical wave tank with a moving wall. This wall can move according to the prescribed law. In this case we model the wave generation process by a piston-type wave maker [38]. In our formulation the wall can also freely move under the wave action. It is realized by connecting the wall to a system of elastic springs. In this way, the wall mass $m$ and springs rigidity $k$ are taken into account in our model by writing an additional nonlinear second order differential equation for the wall position. Surface waves are described mathematically using the full water wave problem, i.e. the irrotational incompressible perfect fluid flow described by Euler equations [73]. No approximation was made and, thus, the presented results are fully nonlinear and fully dispersive. Moreover, we proposed a numerical algorithm on adaptive moving grids. The adaptation was achieved using the so-called equidistribution principle [43, 44]. The discrete problem was rigorously studied and the proposed discrete operator was shown to be self-adjoint on adaptive grids as well as the continuous problem. Thus, the numerical method preserves qualitatively some properties of the continuous operator, which is referred to as the ‘structure preserving’ property [39] in modern numerical analysis.

The proposed algorithm was first validated on the well-known case of a fixed rigid vertical wall. There is enough analytical, numerical [34, 10] and experimental [16] data available to validate the solver. Then, this numerical tool was applied to generate periodic waves in silico by a moving piston-type wave maker. Finally, the wall motion and incident wave run-up on this moving wall was studied using our numerical means. Our results compared quite well with some published data [40]. The influence and effect of various parameters on the wave run-up was investigated.

Perspectives

In the present article the geometry, mathematical formulation and the numerical method have been presented for the two-dimensional configuration. In future works we plan to extend this methodology for 3D flows. A priori, this generalization is going to be tedious, but straightforward.

As another possible direction for future research, we would like to underline that the moving (solid non-deformable) part in our study was a wall. This object is very important in coastal engineering. However, in future works we would like to incorporate moving objects of more general geometrical shape. As a minor point of improvement, one can notice that in all simulations presented earlier the bottom was taken to be flat. It was done on purpose to isolate the effect of the moving wall on incident waves. However, it is not a limitation of the method. The formulation presented earlier works for general bottoms and the simultaneous effect of the moving boundary and local bathymetric features on the wave run-up on this wall has to be investigated.

Finally, the formulation we used was fully conservative. In other words, no dissipation mechanism was included in our study except for the motion of the wall under the action of water waves (however, the whole wave/wall system was conservative). Of course, this model is an idealization of the reality. In real flows some dissipative effects are also present due to the molecular viscosity, friction, turbulence or other mechanisms. For this purpose a suitable visco-potential formulation has been proposed [23, 26, 25], which occupies the intermediate level between the potential flows and Navier–Stokes equations.

Side effects

Our work has also at least one nice side effect. Namely, if we simplify the problem considered in the present study by fixing the wall position, e.g. $s(t)\ \equiv\ 0\,$, we obtain a robust numerical wave tank for 2D surface waves similar to [37]. However, the technique presented here is fully adaptive. By the way, we are not aware of any recent study using similar adaptive redistribution methods for the full Euler equations with free surface. So, this aspect might be new.

Acknowledgments

This research was supported by RSCF project No 14-17-00219. D. Dutykh acknowledges the support of the CNRS under the PEPS InPhyNiTi project FARA and project $\mathrm{N}^{\circ}$EDC26179 — “Wave interaction with an obstacle” as well as the hospitality of the Institute of Computational Technologies SB RAS during his visit in October 2015. D. Dutykh would like to acknowledge the help of his colleague Tom Hirschowitz who timely provided tea leaves in cases of emergency.

Appendix A Transformation of coordinates

In this Appendix we explain some calculation details of coordinate transformations which allow to pass from the system (2.12) – (2.20) posed on a variable domain $\Omega\,(t)$ to the system (3.4) – (3.11) on the fixed domain $\mathcal{Q}^{\,0}\,$. Consider a smooth bijective map (3.1) (see also Figure 2 for an illustration):

[TABLE]

We consider that the time variable is the same in both coordinate systems. Two maps written above are mutually inverse. So, let us write this condition

[TABLE]

Then, we differentiate (A.2) and (A.3) with respect to $x\,$. We thus obtain two relations

[TABLE]

The last two relations can be considered as a system of two linear equations with respect to $\dfrac{\partial\/q^{\,1}}{\partial\/x}$ and $\dfrac{\partial\/q^{\,2}}{\partial\/x}$ as unknowns. By trivially applying the Cramer rule we obtain

[TABLE]

where $\mathcal{J}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ x_{q^{\,1}}\,y_{q^{\,2}}\ -\ x_{q^{\,2}}\,y_{q^{\,1}}$ is the Jacobian defined in (3.3). The Jacobian can be also seen as the determinant of a matrix:

[TABLE]

Similarly, by differentiating (A.2) and (A.3) with respect to $y$ and using Cramer’s rule, one can show that

[TABLE]

Now, by differentiating (A.2) and (A.3) with respect to $t$ we obtain the following two equations

[TABLE]

Solving the last linear system with respect to $\dfrac{\partial\/q^{\,1}}{\partial\/t}$ and $\dfrac{\partial\/q^{\,2}}{\partial\/t}$ yields

[TABLE]

This concludes the computation of partial derivatives of the inverse mapping $\boldsymbol{q}\ =\ \boldsymbol{q}\,(\boldsymbol{x},\,t)\,$.

Appendix B Transformation of the Laplace operator

Now we can rewrite the Laplace operator (2.12) in curvilinear coordinates \bigl{(}\boldsymbol{q},\,t\bigr{)}\,. The idea consists in viewing the velocity potential as

[TABLE]

The first derivatives can be easily computed:

[TABLE]

and using the just derived relation (A.4) we can express everything in terms of derivatives solely with respect to $q_{1}$, $q_{2}$:

[TABLE]

Similarly, we can compute the other partial derivative of the velocity potential $\varphi$ with respect to $y$ in new coordinates:

[TABLE]

Now we can apply recursively the just obtained results for $\frac{\partial\/\varphi}{\partial\/x}\,(\boldsymbol{q},\,t)$ and $\frac{\partial\/\varphi}{\partial\/y}\,(\boldsymbol{q},\,t)$ to compute the second derivatives as well:

[TABLE]

By taking into account the relations derived earlier, we obtain that

[TABLE]

Using similar methods we can obtain the following expression for $\frac{\partial\/^{\,2}\,\varphi}{\partial\/y^{\,2}}\,(\boldsymbol{q},\,t)\,$:

[TABLE]

The Laplace operator can be particularly compactly written in new variables if we introduce the metric tensor components:

[TABLE]

along with the implied anisotropic diffusion coefficients:

[TABLE]

Using this notation we can finally write the Laplace operator in transformed coordinates:

[TABLE]

As a result, the Laplace equation reads:

[TABLE]

where we employed an implicit summation over repeated indices141414In the literature this rule is referred to as the Einstein summation convention.. This result coincides exactly with equation (3.4) presented earlier in the manuscript.

Appendix C Boundary conditions transformation

Finally, let us derive some relations which turn out to be useful in the transformation of boundary conditions. In particular, in this Appendix we focus on the kinematic condition (2.4). First of all, during the derivation of boundary conditions one needs time derivatives of the inverse mapping in (A.1):

[TABLE]

These formulas were given also earlier in (A.6) with more details on their derivation. Let us compute Cartesian components of the velocity vector:

[TABLE]

The free surface elevation $\eta$ on the transformed domain becomes a function of the variable $q^{\,1}$ and time $t\,$:

[TABLE]

Now, we can compute partial derivatives of the function $\eta\,$, which appear in the boundary conditions:

[TABLE]

By substituting all these elements into equation (2.4), we obtain the following condition:

[TABLE]

where $v^{\,1}$ is a contravariant component of the velocity, which can be computed according to this formula:

[TABLE]

The last equation (C.1) completes the derivation of the free surface kinematic boundary conditions. Other boundary conditions can be derived in a similar manner.

Remark 7**.**

Above, in the main text of our manuscript we do not use the tilde notation $\tilde{\eta}$ in order to avoid overloading of the text with multiple notations. It is assumed that the reader can make the difference between $\eta\,(x,\,t)$ and $\tilde{\eta}\,(q^{\,1},\,t)$ depending on the context.

Appendix D Piston motion modelling

Let us introduce a Cartesian coordinate system with the horizontal axis $O\,\chi$ looking rightwards such that the left moving wall is located at the point $\chi\ =\ 0$ in the absence of any kind of loading (only the hydrostatic atmospheric pressure that we neglect here). In this case the wave tank is dry. Now we fill it with still water of the uniform depth $h_{\,0}\ >\ 0\,$. The hydrostatic force $\mathscr{F}\,(0)\ >\ 0$ (defined in equations (2.11)) appears on the left wall and the springs accumulate some elastic deformation. It results in the wall displacement of magnitude $\Uptheta\,(0)\ >\ 0$ in the negative direction of the axis $O\,\chi\,$. The displacement $\Uptheta\,(0)$ can be found by applying the Hook law [51]:

[TABLE]

where $k$ is the springs rigidity. For $t\ >\ 0$ we have another force $\mathscr{F}\,(t)$ acting on the left wall due to the wave motion. This force acts against the wall. Thus, in the second law of Newton we take it with the opposite sign [62, 50]:

[TABLE]

where $m$ is the wall mass. The initial conditions for this ODE are:

[TABLE]

Now, we make a change of coordinates $0\,\chi\ \rightsquigarrow\ O\,x$ so that the wall is located initially at the point $x\ =\ 0\,$, i.e.

[TABLE]

Then, the initial conditions become by construction:

[TABLE]

Equation of motion (D.1) transforms to:

[TABLE]

The last equation can be rewritten as

[TABLE]

Finally, by denoting the wall displacement by $s\,(t)\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ x\,(t)$ we recover equation (2.10).