Variational finite element methods for waves in a Hele-Shaw tank☆

The damped motion of driven water waves in a Hele-Shaw tank is investigated variationally and numerically. The equations governing the hydrodynamics of the problem are derived from a variational principle for shallow water. The variational principle includes the effects of surface tension, linear momentum damping due to the proximity of the tank walls and incoming volume flux through one of the boundaries representing the generation of waves by a wave pump. The model equations are solved numerically using (dis)continuous Galerkin finite element methods and are compared to exact linear wave sloshing and driven wave sloshing results. Numerical solutions of the nonlinear shallow water-wave equations are also validated against laboratory experiments of artificially driven waves in the Hele-Shaw tank.


Introduction
The study of water waves has been an important area of research for years; their significance becomes obvious when looking at ocean and offshore engineering or naval architecture. A well known approach when studying the generation and evolution of nonlinear water waves is to use potential flow theory; that is, to assume that the fluid velocity is incompressible and irrotational, which allows to define it as the gradient of a scalar called the velocity potential. The potential flow model results from the incompressible Euler equations (in which the effect of water viscosity is neglected) and consists of a Laplace equation for the velocity potential coupled with nonlinear free surface boundary conditions. The dynamics of water waves is Hamiltonian and the governing equations can be also obtained from Luke's variational principle [1] or from its dynamical equivalent presented by Miles [2] .
A simplification of the classical potential flow theory can be obtained for shallow water waves. The shallow water theory excludes dispersion and reduces the dynamics to the free surface only, by providing solutions for the free surface elevation and the velocity potential at this free surface. There are, however, asymptotic theories which extend the shallow water theory to include dispersive effects, for example the Boussinesq and Benney-Luke water-wave equations [3][4][5][6][7] .
The variational principle for water waves can be extended to include forcing and/or dissipation, essentially by adding time-dependent internal or boundary conditions. This is necessary for most practical problems involving waves generated by a driven wavemaker or the removal of a sluice gate; for example, see Bokhove and Kalogirou [8] for a problem concerning water waves generated by removing a sluice gate in a wave tank. Likewise, dissipative wave problems were studied by Thornton et al. [9] and Bokhove et al. [10] , who investigated the dynamics of gravity water waves in a vertical Hele-Shaw cell with a moving interface between water and air. The Hele-Shaw cell [11] is a narrow wave tank consisting of two closely spaced glass plates. The motion of waves in this cell is susceptible to momentum damping due to friction, which becomes important because of the proximity of the two plates. Thornton et al. [9] presented results regarding damped wave sloshing, while Bokhove et al. [10] also included a wave pump and studied the damped motion of driven water waves.
Hele-Shaw flows are nearly two-dimensional and this is their great advantage considering that they have been mainly used for visualisation of fluid flows. Such flows have become of interest in connection with applications in coastal engineering [12] , biomedical engineering [13] and the oil industry [14] . Several researchers have investigated the stability of the interface between two viscous fluids in a horizontal Hele-Shaw cell (e.g. Saffman and Taylor [15] , Park and Homsy [16] , Park et al. [17] ), but others also studied free surface waves (e.g. Rajchenbach et al. [18] ) or the Kelvin-Helmholtz instability (e.g. Plouraboué and Hinch [19] ) in vertical Hele-Shaw cells. Classical Hele-Shaw flows are based on the assumption of low Reynolds number (Taylor and Saffman [20] ) and the dominant flow velocity profile between the plates is parabolic, i.e. of Poiseuille-type [21,22] . Based on the parabolic flow profile approximation in the lateral direction, the governing equations in a vertical Hele-Shaw cell can be derived by width-averaging the Navier-Stokes equations.
In this study, in contrast to the classical case, we retain inertial effects due to the proximity of the glass plates and we extend the variational principle to include an exponential time-dependent term representing linear damping effects. The use of linear damping is validated in Thornton et al. [9] , where the authors compared numerical solutions to damped wave sloshing experiments. The dynamics thus become non-autonomous due to the explicit time-dependence and therefore this work is based on the theory developed in Bokhove and Kalogirou [8] for non-autonomous variational principles in time. Starting from the one-dimensional shallow water equations, we add the effects of surface tension, linear damping and forcing through the boundary (representing the motion of the wave pump). The continuum model is discretised in space and time using a (dis)continuous Galerkin finite element method and the numerical model is then implemented using the automated system Firedrake [23] . Preliminary results are first compared to exact solutions for validation purposes and then to wave structures observed during laboratory experiments in a Hele-Shaw cell ( Fig. 1 ).
The structure of the rest of this paper is as follows: Section 2 presents mathematical models for waves with surface tension and forced-dissipative shallow water waves. In Section 3 , we provide the numerical methods and results, including numerical validation for linear undriven and driven wave sloshing. Section 4 contains a description of the set-up and design of the experiments and details about the analysis of the experimental data. We also compare results obtained numerically to experimental observations. We conclude in Section 5 by summarising the outcomes of this study.

Wave models with surface tension
We consider an inviscid fluid in a one-dimensional domain x ∈ [0, L ] with hard walls at x = 0 and x = L . A variational principle for the Lagrangian density L describing the dynamics of shallow water waves with surface tension can be obtained by extending the well-known Luke's variational principle [1] and is given by: with water depth h = h (x, t ) and velocity potential φ = φ(x, t ) as functions of the spatial coordinate x and time t . Partial derivatives are denoted by ∂ x and ∂ t . Here, g is the acceleration of gravity, ρ is the density, σ is the surface tension, h 0 is the rest water depth over a flat bottom and T > 0 is the final time. Dependence on the contact line angle between water, air and the solid glass plates is imposed by prescribing constants E 0 and E L at x = 0 and x = L, respectively. The line element x lies along the free surface and L 0 σ d S is the energy stored by surface tension.
The variational principle in (1) is defined as: and requires the action δ of the integral S [ h, φ] to be stationary for arbitrary small changes in h or φ. Applying the variational operator in (1) thus yields: where we have used the boundary conditions ∂ x φ = 0 at x = 0 and x = L, and end-point conditions δh = 0 at t = 0 and t = T . Assuming that the variations δh, δφ are arbitrary, the dynamics emerge as: Eq. (4) are the standard shallow water equations but with added surface tension [24,25] .

Forced-dissipative shallow water flow
We modify the variational principle (1) to add in-and out-flow at x = 0 and linear momentum damping. The former is a simple model of a wave pump and the latter models the bulk friction in a Hele-Shaw cell, following the derivation in Thornton et al. [9] and Bokhove et al. [10] . This Hele-Shaw cell consists of two narrowly-placed glass plates of 6 mm thickness and it has an open top with a water-air interface. A periodic wave motion is created via pumping of water into the tank by a stilling-well pump. The centrepiece of the wave tank is 57 cm in length, 20 cm in height and 2.2 ± 0.1 mm in width (the glass plates are clamped into a perspex frame using spacers of 2.2 mm, to secure a constant width at every height). A summary of the dimensions of the Hele-Shaw cell, as well as the fluid's physical properties as used in the experiments, can be found in Table 1 .
We non-dimensionalise the problem by rescaling horizontal and vertical lengths with typical length scales L x and H 0 , respectively. The following transformations are thus introduced: where superscript stars are used to denote dimensionless coordinates and variables. A modified variational principle for the Hele-Shaw dynamics then arises by integrating the term ( φ∂ t h )e γ t by parts (the first term in (1) with added damping γ ) and is given below (note that for simplicity the stars are dropped from all terms): Here, γ = 3 ν/ 2 is the damping factor, in which ν is the viscosity and 2 = W is the width of the Hele-Shaw cell. In addition, Q ( t ) is the volume flux in or out of the Hele-Shaw cell at x = 0 , such that the volume flow rate is 2 Q . A time dependent contact angle E ( t ) is also used at x = 0 , while at the other end x = L we set E L = 0 . We further replace the term Using suitable boundary and end-point conditions, variation of (6) gives: The equations of motion follow from the arbitrariness of the variations in (7) and are given by: The second term in the momentum Eq. (8a) represents the effect of linear damping, causing the wave energy to damp out in time.

Time discrete weak formulations
To discretise the variational principle (6) , a Galerkin-Ritz finite element expansion is introduced by writing: with basis functions ϕ k ( x ) and indices k , covering the degrees of freedom. The Einstein summation convention over repeated indices is also adopted. The resulting discretisation is implemented in the finite element package Firedrake [23] . The weak formulations arise from the discrete variational principle (6) , but this time we do not integrate by parts like before; we have the following instead: Note that here we choose to discretise the variational principle directly and not the resulting system of nonlinear evolution equations. The two approaches can be equivalent; in our chosen way, the weak formulations arise by taking the arbitrary variations δh h , δφ h to be equal to the test functions, i.e. δh h = ϕ , δφ h = 0 or δh h = 0 , δφ h = ϕ .
As in Bokhove and Kalogirou [8] , the time stepping scheme is put in the standard "Kamiltonian" form for variational principles with explicit time dependence by a coordinate transformation P h = φ h e γ t and Q h = h h . On implementing the scheme, let t be the time step, T the final time and N = T / t; then t n = n t, for n = 0 , 1 , 2 , . . . , N, are the time levels at which the solutions φ h and h h will be calculated. A Störmer-Verlet time discretisation (Hairer et al. [26] ) for this "Kamiltonian" formulation then consists of a half step to update φ h in time, and a final half step for φ h , Note that the test functions δh h and δP h are space dependent functions only. The conserved Hamiltonian, in the absence of the pump when Q (t) = 0 , at each time step is defined by the sum of kinetic and potential energies and is given by:

Numerical validation
To verify the numerical implementation, we compare numerical solutions with a standing wave solution of the linearised equations. The linear code is further verified in comparison with an exact steady state solution emerging from the forced shallow water equations with a (nearly) harmonic forcing for Q ( t ). The effect of surface tension on such driven waves is also investigated. Note that in the rest of this paper we use the condition ∂ x h = 0 at x = 0 , such that the contact angle E(t ) = 0 . All simulations were done on a uniform mesh, to third-order in space using quadratic Lagrange polynomials and second-order in time using Störmer-Verlet or modified-midpoint symplectic schemes. Unless otherwise stated, 50 elements are used with time step t = 0 . 0028 (chosen such that to satisfy the stability criterion | ω t | < 2, where ω is the maximum frequency for a wave to travel from one grid point to another).

Damped wave sloshing
We linearise around the rest solution h = h 0 , φ = 0 and use a hard-wall boundary condition ∂ x φ = 0 at x = 0 , so that Q (t) = 0 . In the absence of damping, i.e. for γ = 0 , the linear equations can be solved exactly and a standing wave solution is found to be: with k = 2 π m L (for integer m > 0), wave amplitude A and frequency ω 2 = h 0 k 2 (1 + σ k 2 ) .
The above solution is used in a comparison with the numerical solution of the linearised undamped shallow water equations (with γ = 0 ). The visual comparison after two wave periods is illustrated in Fig. 2 , where several snapshots of the numerical solution (plotted with solid line) and the exact solution (plotted with circles) can be seen. Moreover, the deviation of the total energy H ( t ) (as defined in Eq. (12) ) from the initial energy H(0) = H(t = 0) can be seen in the bottom panel of the figure, which remains bounded and shows no drift, as expected.
In the presence of damping, dissipation of energy is observed due to the side-wall friction. This can be explained by looking at the dispersion relation, which is now given by ω 2 + iγ ω − h 0 k 2 (1 + σ k 2 ) = 0 and admits the complex solutions: (assuming that the determinant is positive). The negative imaginary part of the dispersion coefficient causes a decay in the wave amplitude over time, and therefore it is responsible for the loss of energy. Mathematically, the energy can be cast into a conserved form by adding an exponential factor e γ t and calculating the Hamiltonian as in (12) . Fig. 3 demonstrates a comparison between the energies for linear undamped waves with γ = 0 and linear damped waves with γ = 0 . 058 ; the former is shown in the top panel and the latter in the middle and bottom panels. The dissipation of energy due to damping is obvious when looking at the middle panel, where the energy is computed without the exponential factor e γ t . Nevertheless, the energy including the exponential factor (as in (12) ) is conserved, up to bounded oscillations vanishing when t → 0, as shown in the bottom panel.   (12) ).

Driven wave sloshing
We solve the linearised shallow water equations while imposing a flux boundary condition at x = 0 given by Q (t) = q sin (ωt) , such that Q (t) = h 0 ∂ x φ at x = 0 . The exact standing wave solution in this case is found to be: with k = (2 m +1) 2 π L (for integer m > 0) and again ω 2 = h 0 k 2 (1 + σ k 2 ) . This exact solution is compared against the numerical solution of the linearised and driven shallow water equations and the result of the comparison is shown in Fig. 4 . The fluctuations in the plot of the total energy are again small and bounded.

Driven waves with surface tension
We also solve the forced nonlinear shallow water equations (undamped) for zero and positive values of the surface tension in order to investigate its effect on the flow. The fluid initially is at rest and then the wavemaker is switched on at  t > 0. The generated wave comes in from the left and propagates to the right. The result can be seen in Fig. 5 . The exact standing wave solutions of the linearised equations are also plotted at a specific time (different in each case because the wave frequency changes for σ = 0, see Eq. (15) ) to show that even though the numerical solution is nonlinear, it still tends to behave like the analytical solution. We have further verified that the energy variations remain small and bounded. Note that simulating water waves with surface tension comes with a time step restriction of the order of ( x ) 2 .
The observed effect of surface tension follows the Young-Laplace equation p = σ κ, which states that the variations in hydrostatic pressure are proportional to the curvature of a surface. In the absence of surface tension (top panel), the free surface at the right end of the tank remains almost flat before the wave reaches that area, because in this case there are no variations in the pressure and p = 0 . The free surface in the same region has a bigger curvature if surface tension is present (bottom panel), and in this case the hydrostatic pressure is increased. Note that in the latter case the value of surface tension used is relatively large compared to the surface tension of water (about 10 −5 in dimensionless units), chosen such that to study its effect on the free surface.

Experimental set-up
To validate our variational finite element model, we will compare simulations thereof against laboratory measurements of a new Hele-Shaw cell for wave-beach interactions, see Fig. 1 (the new tank is an alteration of the Hele-Shaw tank in [9,10] ). In this study, we are only interested in the hydrodynamics and not the beach dynamics, so we operate it without the beach particles and only with water. A similar comparison has been presented by Sun and Hsu [27] who also carried out experiments in a vertical Hele-Shaw device and computed the wave number and rate of decay along propagation distance of free surface waves. However, their setting was different from ours in the sense that the waves were generated far away from the Hele-Shaw cell (which was installed at the centre of a wind-wave channel), and this resulted in the creation of incident, reflected and diffracted waves. Here we conduct the experiments in a way such that the waves come from one direction only and reflect on the hard walls of the wave tank.
The generation of waves in our Hele-Shaw tank is controlled by a stilling-well pump described next. On the left and right sides of the tank there are two stilling-well chambers of width 13.5 ± 0.1 mm, connected to two hoses. The left stilling well is connected to two aquarium pumps via hoses, placed in parallel in a bucket of water hanging below the table on which the wave tank has been placed, such that the height between bucket and wave tank can be adjusted. The two aquarium pumps are connected to a motor controller, fed by a power supply of 12 V and 2 A maximum, which can be loaded onto an Arduino board and programmed by a laptop using the Arduino program. These pumps operate one way and the pump action is repeated in loops. When the pumps do not operate, the wave tank empties out till the pumps get turned on again. Their strength is determined by analog dials for the current and voltage on the power supply. The mean water level can be set by adjusting the frequency of the pump, the current and voltage applied, and the distance of the water level in the bucket below the wave tank. Every time we set up the wave tank, this requires some heuristic adjustment: when the pump action is too strong the wave tank floods, meaning that either the current, the voltage and/or the bucket must be lowered. An equilibrium state with a mean water level can be found within a certain parameter range.
We performed several runs lasting between 20-60 s each. The frequency of the generated waves was in the range of 0.07-0.83 Hz, with the low frequency waves being smooth and of small amplitude, while the ones with higher frequency were much steeper and even breaking waves in some cases. We therefore choose to concentrate on a case with observed waves being smooth (non-breaking) and of relatively high amplitude.

Data analysis
The goal is to validate the mathematical model by comparing our numerical results with the experimental ones. To this end, we need to estimate the volume flow rate by analysing the snapshots from the experiment and accurately compute the volume flux function Q ( t ) used in the model (e.g. see Eq. (6) ). This was achieved by using image processing techniques in Matlab. Knowing the camera's resolution (1080 × 1920 pixels) allowed us to estimate the number of pixels that define the Hele-Shaw cell in the horizontal and vertical directions. We then used an edge detection code to detect discontinuities in brightness, essentially specifying the coordinates of the free surface for every snapshot. The integral of the line defining the free surface (i.e. the area under the free surface) times the width of the tank, defines the volume of water contained in the tank at each time. We finally calculated the displaced area at sufficient discrete times t i and over a few wave periods.
The plot of the displaced area of water A against time can be seen in Fig. 6 . From the figure we can extract useful information about the experiment, like the wave frequency F and the mean water depth h 0 . For this particular case, we find F = 0 . 7042 Hz, while the mean water level is determined by the time average of A over a suitable number of periods.

Comparison with numerical simulations
To determine the boundary conditions for the numerical model, we calculate the area A = A (t) and the water height h (0, t ) on the left side of the wave tank with the operable stilling-well pump. We can then determine the volume flux Q ( t ) ∝ d A /d t at x = 0 into the tank. We finally solve the mathematical model numerically with a volume flux boundary where Q a is the amplitude and F is the frequency of the observed wave during the experiment. To find the amplitude Q a of the volume flux, we take the derivative of the artificial signal ˜ A (t) and thus find Q a = A a × 2 π F = 0 . 0204 m 2 /s. We perform numerical simulations using the finite element methods presented in Section 3 with 100 elements and time step t = 0 . 001 ; the methods utilise quadratic Lagrange polynomials, yielding third-order in space, and the second-order Störmer-Verlet time integrator. We also employ a first-order Godunov finite volume method for the shallow water equations, applying the HLL-flux and using 800 cells (the finite volume solver uses h (0, t ) and Q (t) = hu | x =0 as boundary conditions at x = 0 ). We concentrate on the experiment associated with Fig. 6 , where the wave frequency is F = 0 . 7042 Hz. The result of the comparison between experiment and numerical simulations is demonstrated in Fig. 7 , where on the left we show three snapshots from the experiment and on the right we plot the corresponding simulated wave profiles. The snapshots on the left are chosen from Fig. 6 and are associated with states of roughly minimum, average and maximum volume of water in the cell, respectively. In particular, the presented snapshots correspond to approximate times t = 1 . 5 , 1 . 75 , 2 . 0 s and illustrate three stages of the wave propagation: first a smooth wave comes in from the left and the water level slowly increases; the wave then travels along the tank and when it reaches the right end, it reflects on the wall and travels to the opposite direction.
The numerically computed wave profiles were chosen after the required tuning to match the position of the wave front. We plot the profiles obtained via the finite element method, as well as the finite volume solver to confirm the steepness of the wave, which is close to becoming a shock. The observed ripples ahead of the wave front are resolved, and we believe these to be due to the effect of surface tension (absent in the finite volume computations). The simulation is left to run for approximately 4 wave periods and we plot the numerically computed total energy in Fig. 8 , where the dashed and solid curves correspond to the energy obtained using the finite element method and the finite volume method, respectively. The energy computed with the finite element method is slightly higher than the one we get with the finite volume method, due to the presence of surface tension which adds the last term in Eq. (12) . Note that the total energy is computed without the exponential term e γ t and clearly remains bounded due to the balance between forcing and damping.

Conclusions
The motion of the free water surface in a Hele-Shaw cell is studied by variational modelling and experiments. Despite the popularity of horizontal Hele-Shaw cells, the use of a vertical Hele-Shaw cell to study the hydrodynamics of water waves is fairly novel [18,27] , and takes advantage of the fact that the flow becomes nearly two-dimensional and offers better visualisation of the wave phenomena. In this paper, we derive a mathematical model to describe the hydrodynamics of shallow water waves in such a Hele-Shaw tank. Novel features in this model (compared to the usual shallow water equations) include the effect of surface tension, linear momentum damping and in-or out-flow conditions at one boundary representing the motion of a driven wavemaker.
The model is discretised in space using finite element methods and in time using symplectic integrators. As a validation, we compare numerical solutions to wave sloshing and driven wave sloshing exact solutions. We finally simulate the motion of driven and damped water waves in a Hele-Shaw tank and compare against wave shapes observed in experiments using our own device. A qualitative comparison between the numerical and experimental waves shows a reasonable agreement between the two, even though the numerically computed profiles have a steeper wave front (confirmed by comparisons with a finite volume solver) and capillary waves -the characteristic surface tension ripples -appear ahead of the wave. Nonetheless, this preliminary validation of simulations against the observed waves in the Hele-Shaw cell is satisfactory, and could be improved by using the actual boundary data for Q ( t ) and h (0, t ) in the simulations instead of the artificial data used here. Furthermore, the difference between the observed profile, the computed hydraulic jump in the shallow water model and the dispersed capillary waves in the finite element model suggests that long-wave dispersion may be needed, and possibly also damping due to the contact line.
During the experiments, we noticed that for frequencies higher than about 1 Hz, the free surface was very steep and for even higher values of the frequency breaking waves emerged. The mathematical model breaks down in these cases because it is derived based on the assumption of irrotational waves with a smooth free surface. To overcome this problem, the twophase Navier-Stokes equations or other approximations thereof should be considered. This will be essential also in terms of adding wave dispersion, which is absent in shallow water modelling. Furthermore, in order to incorporate the assumption of a Poiseuille-type lateral flow profile in the width-averaged Navier-Stokes equations [9,18] , the factor of β = 6 5 should appear in front of the nonlinear terms in the Navier-Stokes and continuity equations. We leave this as future work.
A further extension of this subject is the study of beach formation by breaking waves in a Hele-Shaw cell [12] . This can be achieved by adding small particles in the tank, essentially representing a "slice" of beach. The Hele-Shaw beach configuration can be useful in benchmark modelling of one-dimensional breaking waves and also in the study of dynamic shingle beaches. studentship. Special regards go to Wout Zweers for constructing the Hele-Shaw tank used in the experiments, and Gareth Keevil and Rob Thomas for their help with setting up and conducting preliminary experiments. We also thank the Firedrake project team, in particular Colin Cotter, David Ham and Lawrence Mitchell, for assistance in using Firedrake.

Supplementary material
Supplementary material associated with this article can be found, in the online version, at 10.1016/j.apm.2016.02.036