Run-Up Simulation of Impulse-Generated Solitary Waves

Impulse waves in natural lakes or reservoirs are generated by landslides or avalanches and may lead to great damage in Alpine valleys (Miller 1960; Schnitter and Weber 1964; Fuchs and Boes 2010; Brideau et al. 2012; Carey et al. 2012). Plane impulse waves have been investigated at the Laboratory of Hydraulics, Hydrology and Glaciology (VAW), ETH Zurich, leading to a general description of the wave generation and propagation processes used to predict the main wave features as a function of the governing slide impact parameters. Fuchs and Hager (2012), Fuchs (2013), Fuchs and Hager (2015), and Hafsteinsson (2014), generated the waves using a pneumatic piston-type wave generator. The wave heights and wave run-up heights on a linearly inclined slope were measured optically and pointwise using ultrasonic sensors and capacitance wave gauges.

This paper numerically simulates these paddle-generated impulse waves. Because the breaking of these water waves can be observed, Eulerian methods have been widely used, either volume of fluid (VOF) methods (Hirt and Nichols 1981; Gopala and van Wachem 2008; Pilliod and Puckett 2004; Pringle et al. 2016; Wroniszewski et al. 2014), or level-set methods (Sethian 1999; Osher and Fedkiw 2006). Lagrangian meshless methods have also been proven to be a robust alternative (Khayyer et al. 2008; Xu 2016). For Eulerian methods, the use of a fine grid is required to compute the liquid–air interface accurately, which may lead to prohibitive memory requirements and computing time. This limitation is overcome by using dynamic adaptive meshing. A natural and efficient approach to achieve adaptive meshing in the framework of structured grids is to use an octree structure (Samet 2005). Free-surface flow octree solvers have been combined either with the level-set method (Strain 1999; Losasso et al. 2004; Nikitin et al. 2011) or with a VOF approach (Popinet 2003). In Laurmaa et al. (2016), for a given velocity field, a semi-Lagrangian octree scheme was described for the displacement of free surfaces, avoiding the Courant–Friedrichs–Lewy (CFL) condition. this paper couples their approach to a finite-element Navier–Stokes solver in order to compute the velocity and to perform these impulse wave simulations.

The outline of the paper is as follows. First, the model is presented. The Navier–Stokes equations (NSEs) are coupled with the transport equation for the volume fraction of liquid $\varphi$. A splitting scheme (e.g., Glowinski 2003) enables transport and diffusion to be decoupled. Transport of the volume fraction of liquid $\varphi$ and of the velocity $v$ is performed on an adaptive octree grid using a semi-Lagrangian scheme, thus enabling CFL numbers larger than 1. Diffusion of the velocity $v$ is performed on a fixed, unstructured finite-element tetrahedral mesh. This method has already been successfully applied, without the adaptive octree feature, to Newtonian, viscoelastic, multiphase flows by Maronnier et al. (1999, 2003), Bonito et al. (2006, 2008); Caboussat (2006), Caboussat et al. (2012, 2013), and James et al. (2014). Finally, paddle-generated water waves in a tilted cavity are simulated and compared with experimental data available in Hafsteinsson (2014) and Hafsteinsson et al. (2017). Wave generation, propagation, and breaking are all successfully simulated.

Let $\mathrm{\Lambda }\subset {\mathbb{R}}^3$ be the cavity containing the liquid at all times $0\le t\le T$, where $T>0$ is the final time of simulation. Let $\mathrm{\Omega }(t)\subset \mathrm{\Lambda }$ be the region occupied by the liquid at time $t$, $\mathrm{\Omega }(t)=\{x\in \mathrm{\Lambda };\varphi (x,t)=1\}$, where $\varphi$: $\mathrm{\Lambda }\times [0,T]\to \mathbb{R}$ is the volume fraction of liquid. Let $Q_T$ be the liquid space-time domain $Q_T=\{(x,t)\in \mathrm{\Lambda }\times (0,T);\varphi (x,t)=1\}$. Let $v$: $Q_T\to {\mathbb{R}}^3$ be the velocity field and let $p$: $Q_T\to \mathbb{R}$ be the pressure field. Note that $v$ and $p$ are defined only in the liquid domain. The time-dependent incompressible Navier–Stokes equations are solvedwhere $\rho$ = density; $\mu$ = dynamic viscosity of the fluid; $g$ = gravitational field; and $\epsilon (v)=(\nabla v+\nabla v^T)/2$ is the strain rate tensor. The NSEs are coupled with the advection equationwhich translates the fact that the fluid particles move at velocity $v$. Given the initial volume fraction of liquid $\varphi (x,0)$, the initial liquid domain is $\mathrm{\Omega }(0)=\{x\in \mathrm{\Lambda }:\varphi (x,0)=1\}$, the initial velocity then has to be prescribed in $\mathrm{\Omega }(0)$. The boundary conditions for Eq. (1) and (2) are then as follows. Neglecting forces due to pressure from the gas outside the liquid region and any capillary forces, it is assumed that the interface is stress free. Therefore the following boundary condition holds at the boundary of the liquid region $\mathrm{\Omega }(t)$ that is not in contact with the boundary of the cavity $\partial \mathrm{\Lambda }$where $n$ = outward unit normal to the free surface. If the liquid is in contact with $\partial \mathrm{\Lambda }$, either Dirichlet conditions (no-slip or inflow, for instance) or slip conditions can be imposed.

Maronnier et al. (2003) proposed an implicit time splitting scheme to decouple advection and diffusion phenomena in Eqs. (1)–(3). Advection was handled by a semi-Lagrangian scheme on a structured grid, whereas diffusion (a time-dependent Stokes problem) was solved on a tetrahedral grid with standard continuous, piecewise linear finite elements. Bonito et al. (2006) and James et al. (2014) successfully extended this splitting scheme to more complex viscoelastic and multiphase flows.

Let $N>0$ be the number of time steps and $\mathrm{\Delta }t=T/N$ be the time step, such that $t_n=n\mathrm{\Delta }t$, $n=0,1,\dots ,N$. Assume that the approximation ${\varphi }^n$: $\mathrm{\Lambda }\to \mathbb{R}$ of the volume fraction of liquid is known at time step $t_n$, which defines ${\mathrm{\Omega }}^n=\{x\in \mathrm{\Lambda };{\phi }^n(x)=1\}$, the approximation of $\mathrm{\Omega }(t_n)$, the liquid domain at time $t_n$. Assume that the approximation $v^n$: ${\mathrm{\Omega }}^n\to {\mathbb{R}}^3$ of the velocity field $v$ at time $t_n$ is also available. Then ${\phi }^{n+1}$ and $v^{n+1}$ are computed as shown in the following subsections.

The splitting scheme [Eqs. (6)–(9)] allows for the use of different grids/methods to solve the prediction step [Eqs. (6) and (7)] and the correction step [Eqs. (8) and (9)]. The prediction step (advection) is implemented using the adaptive octree, whereas the correction step (diffusion, Stokes) is solved on a fixed, unstructured tetrahedral mesh (Fig. 1). As stated by Laurmaa et al. (2016), the choice of the adaptive octree is motivated by the capture of fine details at the interface while keeping coarser cells in the bulk of the fluid for efficiency.

A paddle-generated wave in a tilted 3D cavity was used as a benchmark to determine the accuracy of the numerical scheme for solving free surface flows governed by the NSE. Although the wave propagated in a 3D cavity, the cavity was narrow enough that the wave only exhibited 2D features. The experimental wave profile measurements (Hafsteinsson 2014; Hafsteinsson et al. 2017) were provided by the Laboratory of Hydraulics, Hydrology and Glaciology (VAW), ETH Zürich.

The experimental cavity was 14 m long, 0.5 m wide, and 0.7 m high. One of the side walls was lined with observational glass windows and the other wall was lined with smooth polyvinyl chloride along with the bottom. A pneumatic wave-generating piston was mounted at one end of the cavity and was controlled with electrical impulses sent by a computer. The setup was successfully used in several experiments by VAW (Fuchs and Hager 2012; Fuchs 2013; Fuchs and Hager 2015).

Fig. 6 illustrates the experimental setup. The cavity is tilted at angle $\beta$. A plate pushes the water parallel to the bottom of the cavity with a velocity profile designed to generate solitary waves with a given height relative to the still water depth. The relative wave height $R_h=a/h_0$ is then the ratio of the initial wave height $a$ to the still water depth $h_0=0.2\mathrm{m}$. The run-up height is $r$. In order to define more accurately what is meant by still water depth, define the $x$-axis as parallel to the bottom of the cavity, directed from the paddle to the other end of the cavity. The paddle makes a sweeping motion from $x=-x_0$ to $x=x_0$. The still water depth is the maximum water depth when the plate is at position $x=0$ in the middle of the plate sweep. The start and end coordinates of the plate depend on the desired relative wave height $R_h$. In the following analysis, all combinations of parameters $\beta =1°,6°$, and $R_h=0.3,0.5$, and 0.7 were simulated.

Two types of sensors were used to measure the wave profiles, ultrasonic distance sensors (USDS) and capacitance wave gauges (CWGs). The USDS were placed at the top of the cavity and sent an acoustic signal downward toward the water surface. The signal reflects off the water free surface and the free surface height was deduced from the time taken for the signal to return to the sensor. Hafsteinsson (2014) gave an estimate of 3 mm for expected measurement errors. However, this sensor can give spurious values if the signal is reflected away from the sensor. This can occur in particular if the free surface under the sensor is at a steep angle. Because this happens when $\beta =6°$, the less accurate CWG sensors were also used.

Capacitance wave gauges are vertical enamel-coated wires attached to the top of the cavity whose capacitance varies linearly with respect to water depth. After calibration, they can be used to determine water depth. The obtained wave profiles are noisy and show jumps of up to 8 mm between sampling points, and Hafsteinsson (2014) estimates the measurement error to be also up to 8 mm. For this reason, they were smoothed with a Savitzky–Golay filter (Savitzky and Golay 1964). The Savitzky–Golay filter smooths the signal by fitting low-level polynomials with the linear least-squares method with successive subsets of adjacent data points. The polynomials were also used to derive a smooth velocity which was used as a Dirichlet boundary condition at the paddle location for the numerical simulation. As parameters for the filter, a data point window size of 71 data points and polynomials of degree 3 were used. Numerical wave profiles were also smoothed for visualization clarity with the Savitzky-Golay filter but to a much lesser extent, with a data point window size of 7 data points and polynomials of degree 3.

A high-speed camera captured photographs which were used to compare static-in-time wave profiles with the simulated wave profiles. The camera was located between Sensors 3 and 4; the CWG wires can be seen hanging from the top of the cavity. The photographs were fitted together by matching initial still water heights and CWG wires with their known positions in the numerical cavity. Wave slices of the simulated wave matched the camera-side window of the experiment cavity. Although the water level at rest is not shown, the initial water level of the simulated wave and the water-air interface of the real wave at the camera-side window did indeed match. Distortion and refraction effects due to the viewing angle of the camera were less than 1 mm between the middle of the photo and its peripheral areas.

The physical properties of the water for the Navier–Stokes simulations were $\rho ={10}^3\mathrm{kg}/{\mathrm{m}}^3$ and $\mu {=10}^{-3}\mathrm{kg}/(\mathrm{ms})$. Because the wave only exhibited two-dimensional features, the cavity was truncated to a width of 0.125 m (instead of 0.5 m). The base oct of the cavity corresponded to parameters $B_x=235$, $B_y=1$, and $B_z=12$, where $B_x$, $B_y$, and $B_z$ are as defined in Fig. 3. Three meshes/octrees—coarse, medium, and fine—were used to check convergence. The coarse setting corresponded to a finite-element tetrahedral mesh of size $H=0.0102$; to an octree with parameters $l_{\mathrm{liquid}}=1$ and $l_{\max }=3$, which resulted in an octree of minimum size $h=3.79{10}^{-3}$ (the octree size is $4h$ in the bulk of the liquid); and to a time step $\mathrm{\Delta }t=0.0125\mathrm{s}$. This setting corresponded to a CFL number close to 4 on the octree, which was a good trade-off between accuracy and computing time. To check convergence, the mesh size, octree size, and time step were divided by two (medium setting) and four (fine setting). Thus the medium setting corresponded to $H=0.0051$, $l_{\mathrm{liquid}}=2$, $l_{\max }=4$, $h=1.89{10}^{-3}$, and $\mathrm{\Delta }t=0.00625$, and the fine setting corresponded to $H=0.00255$, $l_{\mathrm{liquid}}=3$, $l_{\max }=5$, $h=0.949{10}^{-3}$, and $\mathrm{\Delta }t=0.003125$. Because $l_{\mathrm{liquid}}=l_{\max }-2$, the cells in the bulk of the liquid were approximately of the same size as tetrahedra, whereas the cells on the liquid–air interface were four times smaller. Boundary conditions for the velocity were set to the imposed Dirichlet conditions at the paddle location and to slip conditions at the walls of the cavity. The paddle movement was simulated in practice by computing at each step which nodes were intersecting the region that was swept by the paddle and imposing the paddle velocity at those points. The octree was aligned with the water level and therefore not with the cavity. It was refined to level $l_{\max }$ at the bottom of the cavity to capture the cavity slope at the wall in contact with the paddle and at the water free surface, but not on the lateral sides which were aligned with the octree, and hence captured exactly. Fig. 7 justifies the choice of the octree parameters $l_{\mathrm{liquid}}$ and $l_{\max }$. The computed wave corresponding to $R_h=0.7$ and $\beta =6°$ was reported at time $t=1\mathrm{s}$ for several values these two parameters. Results indicate that the choice advocated in Fig. 1 yields accurate results.

Figs. 8–19 show results for the three mesh settings: coarse (fineness 0), medium (fineness 1), and fine (fineness 2). For each combination of experimental parameters $\beta =1°,6°$ and $R_h=0.3,0.5$, and 0.7 the evolution of the wave height over time was first plotted at four fixed points in the cavity and snapshots of the numerical wave were overlaid on photographs. Four pairs (USDS and CWG) of sensors were placed at coordinates $x_1=0.677$, $x_2=1.177$, $x_3=1.677$, and $x_4=2.177\mathrm{m}$, the last two locations being visible on the photographs. The sensors provided water height data from times $t=0–4\mathrm{s}$ at a sampling rate of 100 Hz. The pneumatic-driven paddle moved with a slight delay compared with the electrical input it received. Wave profiles were therefore shifted to correct for this fact. Despite this, measurements were available at four different points in the cavity, which allows comparing the wave profiles from numerical simulations with those from the experimental profile.

All computations were performed sequentially on a 3.3 Ghz Intel (Santa Clara, California) Xeon E5-2643 with 32 Gb RAM. The computing time when $\beta =6°$ and $R_h=0.7$ was approximately 1 h for the coarse setting, 3 h for the medium setting, and 47 h for the fine setting. The $47/3$ ratio is close to $16=2^4$, which is the theoretical ratio corresponding to refining the mesh size by a factor of 2 (thus multiplying the number of vertices by a factor $2^3$) and the time step by a factor of 2.

Figs. 8–13 show results for $\beta =1°$. For these experimental parameters, the wave does not break in the data shown. However, for $R_h=0.7$, the wave does break soon after Sensor 4. Very good agreement was achieved between numerical and experimental profiles. For $R_h=0.7$, where the wave presented a sharper and higher crest, the mesh needs to be refined more than for lower, smoother waves in order for the wave profile to be captured. Ultrasonic distance sensors are more accurate than CWGs and should be considered the baseline except when the wave is too steep, making the ultrasonic signal bounce away from the sensor, resulting in spurious artifacts as in the two middle frames in Fig. 12. The simulated wave crest was slightly lower than the real one (Fig. 13), but seemed to converge to the measured wave profile (Fig. 12). Fig. 20 presents the vector plot of the velocity corresponding to Fig. 11(c). Fig. 21 presents the plot of the pressure field corresponding to Fig. 17(b); the pressure differed very little from the hydrostatic pressure.

For $\beta =6°$, the waves broke and caused a splash. These types of simulations are typically more difficult and cannot be simulated by shallow water models. Turbulence, fine-scale drops of water, and air bubbles occur during wave breaking, and these aspects are not captured by the proposed model. Figs. 14–19 show results for these parameters. Because the USDS 4 sensor was located above (right) the still water surface and was dry at the start of the test, its curve should be disregarded. Because the CWG 4 sensor was not calibrated for measurements in this range and is also limited in measuring air–water mixtures after wave breaking (its curves significantly undercut the USDS curves), its curve should also be disregarded.

In some of the results reported when $\beta =6°$ [e.g., Fig. 16(b)], the free surface appears to contain unphysical perturbation (presence of a ragged interface). Fig. 22 shows a close-up of the free surface, cells, and velocity at the vertices of the finite-element mesh. These oscillations remained localized within one tetrahedron and did not spread during the simulation.

Numerical wave profiles showed excellent agreement with the experimental wave profiles despite the fact that the model does not capture turbulence and the air cushion formed beneath the breaking wave (Fig. 19). For $R_h=0.7$, once again the higher, sharper wave crest was more difficult to capture and the slight error in the solitary wave heights caused a more significant tongue shape during breaking. Despite this, water wave profiles quickly matched again after wave breaking.

Table 1 reports the run-up height—denoted $r$ in Fig. 6—for $\beta =6°$ and $R_h=0.3,0.5$, and 0.7. The computed values correspond to the finest setting and are 10% smaller than the experimental values. The coarse and middle meshes yielded 40 and 20% errors, respectively. Thus a finer mesh would be needed to obtain more accurate results.

This paper presented a numerical method to accurately simulate paddle-generated impulse waves. Two space discretizations were used. A dynamic octree was used to solve advection, and a fixed, unstructured finite-element tetrahedral mesh was used to solve diffusion. The method allows the use of CFL numbers larger than 1, although the octree was refined in the neighborhood of the liquid–air interface. Interpolation issues between octree and finite elements were discussed. Experiments using paddle-generated waves in a tilted cavity conducted at VAW were reproduced with accuracy.

The authors thank the Swiss National Science Foundation and the Commission for Technology and Innovation for their financial support. Viljami Laurmaa was supported by SNF grant 143470. Gilles Steiner was supported by CTI grant 14359.1 PFES-ES. All the computations were performed using the industrial cfsFlow software developed by EPFL and Ycoor Systems. Helgi J. Hafsteinsson is acknowledged for conducting the experiments at VAW. The (unknown) reviewers are thanked for constructive remarks.