Impact force of roll waves against obstacles

Abstract

The roll waves in open-channel flow on steep slopes can strike an obstacle with great force. We conducted two-dimensional shallow-water simulations to study the impact force of the waves against structures of various shapes and orientations. The focus is on the front runner of a wave packet developed from spatial instability. The numerical results include the stand-off distance of the bow shock wave, the front face's run-up height and the wave force on the obstacle. The strength of the impact depends on the Froude number of the undisturbed flow and the obstacle's distance from the local disturbance but not much on the form of the perturbation that initiates the instability. The wave force could reach a peak of more than an order of magnitude greater than the force on the structure without the roll waves. However, an obstacle with a sharp and pointy front can deflect the incident waves, significantly reducing the impact force.

1. Introduction

Roll waves are instabilities in open-channel flow on steep slopes. The abrupt wavefront of the roll waves moving at a higher speed than the undisturbed flow can strike an obstacle in its path with great force.

The analytical study of roll waves began with the classical works of Jeffreys (1925) and Dressler (1949). Jeffreys (1925) studied the linear stability of the steady-uniform flow down an inclined plane. Dressler (1949) fitted shocks between periodic smooth profiles for finite-amplitude development of the waves. According to Dressler's approach, one must specify the wavelength to determine the waves’ amplitude and celerity.

Zanuttigh & Lamberti (2002), Que & Xu (2006) and Ivanova et al. (2017) conducted numerical simulations for the nonlinear development of roll waves using the shallow-water equations. Richard & Gavrilyuk (2012) introduced a modified shallow-water model to produce wave profiles in better agreement with the laboratory observation by Brock (1967).

Roll waves’ coarsening toward longer wavelengths was a common observation in laboratory, analytical and numerical investigations. However, searching for a preferred wavelength in the coarsening process remains challenging. Ng & Mei (1994) determined the shortest wavelength corresponding to zero energy dissipation across the shock. Balmforth & Mandre (2004) included the bottom topography in their nonlinear model of roll waves; they found an inverse cascade that does not continue to the longest spatial scale but becomes interrupted over intermediate wavelengths.

Many existing works have focused on the temporal development of roll waves. In this and a previous paper by Yu & Chu (2022), we study the spatial evolution of an instability initiated by a local disturbance. The example in figure 1 shows the contrast between (a) the ‘periodic wave trains’ developed from a temporal instability and (b) the ‘wave packet’ evolving from a spatial instability. The wave trains in (a) initiated by a periodic perturbation in space are not determined until the perturbation wavelength is specified. See a discussion of this issue in Appendix A. On the other hand, the wave packet shown in (b), produced by a local disturbance, is determined by the Froude number of the undisturbed flow and the distance of the wave packet's advance.

Figure 1. The contrast between (a) the ‘periodic wave trains’ developed from temporal instability and (b) the ‘wave packet’ characterized by the front runner ($FR$) developing from spatial instability. The amplitudes of the roll waves on the top and bottom of panel (a) are $(\hat {h}-\check {h})/H=1.20$ and 2.51 depending on the perturbation wavelengths, which are $S_o \lambda /H=12$ and 120, respectively. The undisturbed flow is on a slope $S_o=0.0501$ with a velocity $U=1.03$ m s$^{-1}$, a depth $H=0.00798$ m and a Froude number $Fr = U/\sqrt {g'H}=3.71$ in these contrasting examples of instabilities. The symbols in the figure are $H=$ undisturbed-flow depth, $\hat {h}=$ depth of the roll-wave crest, $\check {h}=$ depth of the trough, $\hat {h}_{FR}=$ the peak of the front-runner depth, $S_o=$ channel slope, $x=$ longitudinal distance from the channel inlet, $t=$ time, $\lambda =$ wavelength and $\lambda _p=$ length of the wave packet.

The most noticeable feature of the wave packet is its front runner, which asymptotically approaches the long-wavelength limit of the Dressler (1949) analytical solution according to Yu & Chu (2022). The front runner in a wave packet also was observed in a mudflow simulation by Liu & Mei (1994) and by Meza & Balakotaiah (2008), who studied the nonlinear waves on vertically falling films dominated by viscosity and surface tension. Meza & Balakotaiah (2008) used the word ‘tsunami’ in exaggeration to describe the front runner in their study of the falling film.

The wave impact against structures is strictly a three-dimensional (3-D) problem. But our simulation for the impact force was two-dimensional (2-D) using the depth-averaged shallow-water equations. We extended the one-dimensional (1-D) shock-capturing scheme used by Yu & Chu (2022) to capture the 2-D steep wavefront and the force acting on the structures of various shapes and orientations. The numerical simulations predict the features of the unsteady flow, including the stand-off distance of the bow shock wave, the front face's run-up height and the impact force on the structure. The simulations also determine the effect of the structure's shape and orientation on the impact force. Finally, we show how a sharp and pointy obstacle front can deflect the incident wave, reducing the impact significantly.

Our 2-D shallow-water modelling of the wave force neglects the viscous and turbulence stresses. The viscous drag due to the boundary layer development around the body can be significant when the wave drag becomes negligible, for example, on a submerged body in deep water. On the other hand, a full 3-D Navier–Stokes simulation is computationally expensive, although it can account for the viscous drag and non-hydrostatic effects. The calculation of the roll-wave impact force against structures – conducted in this study for the first time – is computationally demanding because the computational domain must be large to accommodate the nonlinear evolution of spatial instability. Limited by the finite computational resources, only 2-D simulation is currently feasible. However, the efficacy of the 2-D model is noteworthy – noting that the forces due to waves’ impact often dominate over viscous drag. Aureli et al. (2015) found the 3-D Navier–Stokes simulation determined the dam-break wave force slightly better than the 2-D model compared with their laboratory measurement. Xie & Chu (2019) simulated the tsunami impact force on coastal structures by both 2-D and 3-D models. Their 3-D simulation yielded only a slightly higher force coefficient than their 2-D model.

2. Formulation

The 2-D problem is to find the depth $h$ and velocity components $u$ and $v$ in the $x$ and $y$ directions using the following depth-averaged equations:

(2.1)$$\begin{gather} \frac{\partial h}{\partial t}+\frac{\partial({uh})}{\partial x}+\frac{\partial({vh})}{\partial y}=0, \end{gather}$$

(2.2)$$\begin{gather}\frac{\partial{(uh)}}{\partial t}+\frac{\partial( {u^2 h}+\frac{1}{2}g'h^2)}{\partial x}+\frac{\partial({uvh})}{\partial y}= g'hS_o-\frac{c_f}{2}u\sqrt{u^2+v^2}, \end{gather}$$

(2.3)$$\begin{gather}\frac{\partial {(vh)}}{\partial t}+\frac{\partial({vuh})}{\partial x}+\frac{\partial( {v^2 h}+\frac{1}{2}g'h^2)}{\partial y}={-}\frac{c_f}{2}v\sqrt{u^2+v^2}, \end{gather}$$

where $g$ is the gravitational acceleration, $c_f$ is the quadratic friction coefficient and $S_o=\tan \theta$ is the channel slope. The components of the gravity acceleration are $g'=g \cos \theta$ and $g \sin \theta =g'S_o$ in the $z$ and $x$ directions, respectively. The $x$-axis is downward at an angle $\theta$ relative to the horizontal and the $z$-axis is upward perpendicular to the inclined $x$–$y$ plane. The water depth $h$ is measured along the $z$-axis.

The initial condition is an undisturbed flow of uniform depth $H$ with a longitudinal velocity $U$ downward parallel to the channel bed in the $x$ direction. The balance between gravity and bed friction in the undisturbed flow yields the dimensionless relation

(2.4)\begin{equation} S_o = \tfrac{1}{2} {c_f} {Fr}^2, \end{equation}

where $Fr =U/\sqrt {g'H}$ is the undisturbed-flow Froude number.

3. Roll-wave packet and its front runner

A roll-wave packet, such as the one shown in figure 1(b), is a spatial instability developed from a local disturbance. To reproduce this instability, we introduced a slight modification to the uniform flow, with the depth $h$ and unit discharge $q_x =uh$ specified at the inlet as follows:

(3.1a)$$\begin{gather} h(x=0,y,t)=\begin{cases} H\left[1+\epsilon\sin\left(\dfrac{2{\rm \pi} t}{T}\right)\right] & {\rm for} \ 0< t\leq\dfrac{T}{2}\\ H & {\rm for} \ t>\dfrac{T}{2} \end{cases}, \end{gather}$$

(3.1b)$$\begin{gather}q_x (x=0, y, t) = {Fr} \sqrt{g' [h(x=0,y,t)]^3},\quad q_y (x=0, y, t ) = 0. \end{gather}$$

The depth of the inlet disturbance $h(x=0,y, t)$ is the upper half of a sinusoid of period $T$; the amplitude of the sinusoidal disturbance is $\epsilon$. The inlet $q_x(x=0,y, t)$ is selected to keep the Froude number of the disturbance equal to the Froude number of the undisturbed flow.

This inlet perturbation defined by (3.1) is a ‘type-c disturbance of constant $Fr$’. The six inlet disturbances examined by Yu & Chu (2022) were type-a, type-b and type-c of constant $Fr$ or constant $q$. Type-a is a periodic disturbance of constant amplitude. Type-b is one full sinusoid from $t = 0$ to $T$, while type-c is one half of a sinusoid from $t = 0$ to $T/2$. All six inlet disturbances produced the same front runner according to Yu & Chu (2022). This study considers only the roll waves produced by ‘type-c disturbance of constant $Fr$’.

Figure 1(b) shows an example of a roll-wave packet developed from the type-c disturbance of constant $Fr$ with $\epsilon =0.2$ and $T=0.94$ s on an undisturbed flow with a Froude number $Fr = 3.71$. The flow is convectively unstable as both the front runner and the trailing edges of the roll-wave packet move downstream from the source (Criminale, Jackson & Joslin 2003). The front-runner amplitude $\hat {h}_{FR}/H$ increases with time as the packet advances and disperses in space. The spatial extent of the wave packet $S_o\lambda _p/H$ and the number of steep wavefronts within the wave packet $N_p$ increase with the dimensionless time $S_otU/H$. In this example, $N_p= 3$, 5, 7 and 9, $S_o\lambda _p/H= 52.75$, 137.9, 232.9 and 333.2 and $\hat {h}_{FR}/H=2.79$, 3.62, 3.95 and 4.07, at times $S_otU/H=162.95$, 293.31, 436.71 and 586.63, respectively.

We carried out a series of simulations of the roll waves for undisturbed-flow Froude number $Fr = 3.71$, 4.63 and 5.60. Figure 2 shows all simulation results for (a) the amplitude $\hat {h}_{FR}/H$, (b) the celerity of the front runner ${c}_{FR}/U$ and (c) the spatial extent $S_o \lambda _p /H$ of the wave packet, and their dependence on the front-runner position $S_o x_{FR}/H$ and the Froude number $Fr$. The amplitude and celerity of the front runner approach the long-wavelength limit of the Dressler (1949) solution. The green dashed lines in the figure indicate the asymptotic solutions.

Figure 2. The dependence of (a) the amplitude $\hat {h}_{FR}/H$ and (b) the celerity of the front runner ${c}_{FR}/U$ and (c) the spatial extent of roll-wave packet $S_o \lambda _p/H$ on Froude number $Fr = U/\sqrt {g'H}$ and front-runner distance of advance $S_o x_{FR}/H$. (d) The three tables on the right-hand side list the conditions of four simulation cases for verifying the effects of the perturbation amplitude $\epsilon$ and the friction coefficient $c_f$. The lines and the $\circ$, $\square$ and $\triangle$ symbols denote case 1, case 2, case 3 and case 4, respectively. The black, red and blue colours of the lines and symbols refer to $Fr = 3.71$, 4.63 and 5.60, respectively. The green dashed lines in (a) and (b) represent the long-wavelength limit of Dressler's solution.

The tables embedded in figure 2(d) show the simulation conditions for different undisturbed-flow Froude numbers. Simulation cases 1, 2, 3 and 4 were conducted to examine the effect of bed friction $c_f$ and the inlet perturbation amplitude $\epsilon$.

For the reference case 1, the friction coefficient and channel slope are $(c_f, S_o)=(0.00728, 0.0501)$, (0.00786, 0.0843) and (0.00760, 0.119) for $Fr = 3.71$, 4.63 and 5.60, respectively. These are the same conditions as in the laboratory experiment by Brock (1967) on a smooth channel bed. For a given $c_f$ and $Fr$, Brock (1967) selected the channel slope $S_o$ to satisfy (2.4).

In simulation case 2, we doubled the friction and channel slope $(c_f, S_o)$. Remarkably, the front runner amplitude $\hat {h}_{FR}/H$ and the extent of the wave packet $S_o \lambda _p /H$, marked by the open circles in figure 2, are nearly identical to the results for reference case 1 denoted by the solid lines.

The amplitude of the inlet disturbance for both cases 1 and 2 is $\epsilon = 0.2$. For cases 3 and 4, the inlet perturbations are $\epsilon =0.1$ and 0.05. The change in perturbation amplitude by a factor of 4 has only a negligibly small effect on the roll-wave packet and its front runner, as shown by the triangular and square symbols in the figure.

Therefore, we conclude that the nonlinear development of the roll-wave packet depends mainly on two dominant parameters: the undisturbed-flow Froude number $Fr$ and the front-runner distance of advance $S_o x_{FR}/H$.

4. Numerical method and adaptive mesh

Figure 3 sketches the 2-D simulation problem to determine the impact force of the front runner on a circular prism of diameter $D$. The top panels of the figure show the depth profiles on the centre plane, while the bottom panels show the depth-contour maps. Panels (a) and (b) refer to two different times, one before and the other after the impact of the front runner against the structure.

Figure 3. The depth profiles on the centre $x$–$z$ plane (top) and the depth-contour maps on the inclined $x$–$y$ plane (bottom) (a) before the arrival of the front runner at time $S_o (t-t_{arr}) U/H = -4.02$ and (b) after the impact when the front runner just past the circular prism positioned at a distance $S_ox_b/H=251.2$ ($x_b/D=66.7$) from the inlet at time $S_o (t-t_{arr}) U/H= 1.59$ where the arrival time $S_o t_{arr} U/H=166.32$. The simulation results in this example are for an undisturbed-flow Froude number of ${Fr}=3.71$. The green stars mark the prism's position $(x_b, y_b)$. The yellow lines are iso-depths of $h/H=2$, 4, 6, 8.

We used the Gerris Flow Solver by Popinet (2003, 2011) on an adaptive quadtree mesh that dynamically refines and coarsens depending on the flow feature. The computational efficiency of the quadtree adaptive mesh for 2-D shallow-water flow simulation has been reported in several publications. Liang et al. (2007) and An & Yu (2012) simulated shock reflection by a circular cylinder. They found approximately 80 % saving of the CPU time compared with using uniform mesh. Popinet (2015) simulated Tohoku tsunami and he found 97 % saving of CPU time.

Figure 4(a) shows the mesh on an $x$–$y$ plane for the entire length of the channel at time $S_o(t-t_{arr})U/H=-2.68$ for the undisturbed-flow Froude number $Fr = 3.71$. Panels (b) and (d) of the figure show the close-up views of the mesh at time $S_o(t-t_{arr})U/H=-2.68$ and $+$2.68, before and after the arrival of the front runner, respectively. Panel (c) shows the adaptive mesh at time $S_o(t-t_{arr})U/H=0$ when the impact force reaches its peak.

Figure 4. The adaptive quadtree mesh used in the simulation for ${{Fr}=3.71}$ and $x_b/D=66.7$ ($S_o x_b/H= 251.2$). (a) Global view of the mesh. Panels (b) and (d) are the close-up views of the mesh before and after the impact by the front runner and panel (c) is the close-up view at the instant when the impact force reached its peak.

The shock-capturing scheme for spatial discretization is a generalized MINMOD limiter (GML) which uses the limiter function

(4.1)\begin{equation} \psi = {\rm minmod}\left(\beta\frac{q_i-q_{i-1}}{\varDelta},\frac{q_{i+1}-q_{i-1}}{2\varDelta}, \beta\frac{q_{i+1}-q_{i}}{\varDelta}\right) \end{equation}

to reconstruct the cell-interface flux with a free parameter varied from $\beta =1$ to 2 (Sweby 1984; Nessyahu & Tadmor 1990; Kurganov & Tadmor 2002). In the limiter function, $q_i$ is the cell-centred value of the unknown variable and $\varDelta$ is the size of the $i$th cell. The GML reduces to the dissipative classical MINMOD limiter when $\beta =1$. It becomes the least dissipative SUPERBEE limiter if $\beta =2$. We used $\beta =1.50$, balancing dissipation and computational stability. The implementation of the GML for an adaptive mesh is due to Popinet (2011).

The temporal discretization was a second-order predictor–corrector scheme for the hyperbolic part of the equations and a third-order total variation diminishing Runge–Kutta scheme for the source term (Gottlieb & Shu 1998). The second-order accurate Cartesian cut-cell method of An & Yu (2012) and Causon, Ingram & Mingham (2001) defines the boundary of the circular prism.

The impact on the obstacle is a nonlinear problem of wave scattering. The integration of the hydrostatic pressure force per unit length $f_{hs}=\frac {1}{2} \rho g' h^2$ along the closed path $s$ following the boundary of the obstacle gives the wave force $F_w$ as follows:

(4.2)\begin{equation} F_w ={-} \oint_s f_{hs} \,\boldsymbol{i} \boldsymbol{\cdot} \boldsymbol{n} \, {\rm d}s, \end{equation}

where $\boldsymbol {{i}}$ is the unit vector in the direction of the undisturbed flow and $\boldsymbol {{n}}$ is the outward normal unit vector perpendicular to the obstacle's solid boundary. We ignored the viscous drag due to the turbulent boundary layer around the obstacle in calculating the ‘wave’ force by (4.2).

The wave force $F_w$, adimensionalized by the product of the frontal projection area $HD$ and the stagnation pressure $\frac {1}{2} \rho U^2$, is the wave force coefficient

(4.3)\begin{equation} C_{HU^2} = \frac{F_w}{ \frac{1}{2} \rho U^2 H D }. \end{equation}

Before the arrival of the front runner, the wave force coefficient has an undisturbed value $\bar {C}_{HU^2}$. The force coefficient reaches its peak value of $\hat {C}_{HU^2}$ at the arrival of the front runner. We define the impact-duration coefficient

(4.4)\begin{equation} C_{\mathcal{T}}=\frac{ S_o{\mathcal{T}} U}{H} = \frac{\displaystyle\int_{t_{b}^*}^{t_{i}^*} [C_{HU^2}-\bar{C}_{HU^2}]\, {\rm d}t^*}{\hat{C}_{HU^2}-\bar{C}_{HU^2}}, \end{equation}

where ${\mathcal {T}}$ is the duration of the front-runner impact. The integration in (4.4) is over a period of time from $t_{b}^* = S_o t_{b} U/H$ to $t_{i}^* = S_o t_{i} U/H$ when $C_{HU^2}$ is greater than the undisturbed value $\bar {C}_{HU^2}$. The dimensionless parameter $S_o{\mathcal {T}}U/H$ reflects the impact duration of the front runner against the obstacle. The pink-coloured rectangle in figure 5(a) defines the peak force coefficient $\hat {C}_{HU^2}$ and impact-duration coefficient $S_o{\mathcal {T}}U/H$. It is equal to the area below the solid line (representing the time series) and above the dotted line (indicating the undisturbed-flow value).

Figure 5. (a) The wave force coefficient $C_{HU^2}$, (b) shock stand-off distance $x_s/D$ and (c) centre-plane run-up height $R/D$ as a function of the dimensionless time $S_o(t-t_{arr})U/H$. The red dashed line denotes the analytical solution obtained by Forbes & Schwartz (1981).

The total number of computation cells in the adaptive mesh varies with time. But there are a maximum mesh size $\varDelta _{max}$ and a minimum size $\varDelta _{min}$. In our numerical simulation for the impact force over a computational domain of $L_x \times L_y=162\ {\rm m} \times 3\ {\rm m}$, the numbers of computational cells are $L_x/\varDelta _{min}=27\,648$, $L_y/\varDelta _{min} =512$, $L_x/\varDelta _{max}=1728$ and $L_y/\varDelta _{max}= 32$. We conducted a mesh refinement study using the method of Stern et al. (2001) to determine the accuracy of the simulation. The fractional error of the peak force coefficient is

(4.5)\begin{equation} FE(\%) = \frac{|\hat{C}_{HU^2}-[\hat{C}_{HU^2}]_{\varDelta \rightarrow 0}|}{[\hat{C}_{HU^2}]_{\varDelta \rightarrow 0}} \times 100. \end{equation}

The exact value $[\hat {C}_{HU^2}]_{\varDelta \rightarrow 0}$ was estimated by extrapolation. The detailed calculation for a simulation with $Fr = 3.71$ and $S_ox_b/H=1005$ in Appendix B gives the fractional error $FE(\%)=1.175$%.

5. Impact on a circular prism for $Fr = 3.71$

We analyse the impact on a circular prism at a distance $S_ox_b/H = 251.2$ for an undisturbed-flow Froude number $Fr = 3.71$. In this example, the inlet perturbation has a period $S_oTU/H= 6.08$ and an amplitude of $\epsilon = 0.2$. The smooth channel bed has a friction coefficient $c_f= 0.00728$. The dimensionless diameter of the circular prism is $D/H= 75$. Figure 5 shows the temporal variation of (a) the wave force coefficient $C_{HU^2}$, (b) the shock stand-off distance $x_s/D$ and (c) the centre-plane run-up height $R/D$. The impact force peaks when the front runner arrives at time $S_ot_{arr}U/H= 166.32$. The front runner has an amplitude $\hat {h}_{FR}/H= 2.94$, a velocity of $\hat {u}_{FR}/U = 1.44$ and a celerity of ${c}_{FR}/U = 1.64$ before its impact on the circular prism.

The nonlinear dynamics is the impingement of the front runner on the cylinder, the rapid rise of the impact force to a first peak, then a second peak and eventually back to the undisturbed flow. These stages for $Fr = 3.71$ are identified in figure 5 by the labels a, b, c, d, e, f, g, h, i and j for their occurrences at dimensionless times $t^*=S_o (t-t_{arr}) U/H = -2.70$, $-$0.75, $-$0.41, 0.00, 0.44, 0.65, 1.66, 4.37, 8.43 and 32.7, respectively. Figure 6(a–i) shows the corresponding depth-contour maps of the bow shock waves and the front runners on the $x$–$y$ plane.

Figure 6. Panels (a–i) show the depth $h/H$ contours of the bow shock waves on the $x$–$y$ plane at times $t^*=S_o (t-t_{arr}) U/H = -2.70$, $-$0.75, $-$0.41, 0.00, 0.44, 0.65, 1.66, 4.37 and 8.43, respectively. The corresponding wave force coefficient $C_{HU^2}$, stand-off distance $x_s/D$ and run-up height $R/D$ are shown in figure 5.

The stages represented in figure 6 are: (a) before the arrival of the front runner, (b) the beginning of the front-runner interaction with the bow shock wave, (c) the rise of the run-up to its maximum height, (d) the maximum compression of the stand-off distance and the rise of the impact force to its peak, (e–i) the relaxation of the bow shock wave with the reduction of the impact force and (j) the arrival of the second waves in the packet. We provide further description of these stages in the following subsections.

5.1. Compression and rebound of the bow shock wave

The arrival of the front runner compresses the bow shock wave, leading to a sharp reduction of the stand-off distance from the undisturbed-flow value of $\bar {x}_s/D= 0.203$ to a minimum value ${x}_s/D= 0.137$ that occurs at time $t^*_d= 0.0$. The compression is associated with a rapid rise of the run-up height from the undisturbed-flow value of $\bar {R}/D= 0.0756$ to a peak value of $\hat {R}/D= 0.215$. The wave force coefficient reaches its peak value $\hat {C}_{HU^2}= 10.70$ at the arrival time $t^*_d= 0.0$, but the run-up height reaches its peak $\hat {R}/D=0.215$ earlier at time $t^*_c = -0.41$.

The stand-off distance of the bow shock wave rebounds from the compressed minimum to approach ${x}_s/D= 0.270$, denoted by the red dashed line. The stand-off distance ${x}_s/D= 0.270$ was the analytical solution obtained by Forbes & Schwartz (1981) ignoring the friction effect.

The compression of the bow shock wave occurs from $t^*_b$ to $t^*_d$ and the expansion from $t^*_d$ to $t^*_i$, as shown in figure 5(b). The stand-off distance $x_s$ is defined at the cell centre. Therefore, the curve for $x_s/D$ vs $S_o(t-t_{arr})U/H$ in figure 5(b) manifests itself in a jagged manner as the shock position jumps from one cell to the other. Remarkably, a similar compression and expansion of the stand-off distance occurs at a later time $t^*_j$ on the arrival of the second wavefront in the wave packet.

5.2. Collision of the front runner with the bow shock wave

The front runner and the bow shock wave are both steep wavefronts. Figure 6(a) shows the depth contours of the front runner and the bow shock wave at the time $t_a^* = -2.70$. The depths of the wavefronts at this time are $h_{FR}/H = 2.96$ for the front runner and $\hat {h}_{BS}/H = 4.76$ for the leading edge of the bow shock wave. Then, the depth rises sharply at the times $t^* = -0.75$, $-$0.41, 0.00, as shown in figures 6(b), 6(c) and 6(d), respectively. As shown in figure 5, the depths in front of the circular prism are $R/D={0.129,\, 0.215,\, 0.201}$, which are ${R}/H = 9.68$, 16.1, 15.1, at these times varying from $t^* = t^*_b$ to $t^*_d$ when the front runner crashes into the bow shock wave. The dynamic merger of the two steep wavefronts is not additive. It produces a significantly higher bow shock wave, leading to the run-up ${R}/H=16.1$ at time $t_c^* = -0.41$ and the peak of the impact force $\hat {C}_{HU^2} = 10.7$ at time $t_d^* = 0.00$. The production of the significant bow shock wave may explain the exceedingly large impact force acting on the obstacle.

5.3. Wave impact coefficient

Before the arrival of the front runner, the wave force coefficient has a value $\bar {C}_{HU^2}= 1.52$. As the result of the front-runner impact, the wave force coefficient rises to a peak value of $\hat {C}_{HU^2} = 10.7$, more than six times greater than the value before the impact.

The integration of the force coefficient $C_{HU^2}$ using (4.4) over the time from $t^*_b = -0.75$ to ${t_{i}^*}= 4.35$ gives an impact-duration coefficient of $C_{\mathcal {T}}={ S_o{\mathcal {T}} U}/{H}= 2.142$. The pink-coloured rectangle in figure 5(a) defines the impact equivalent with a peak force coefficient of $\hat {C}_{HU^2} = 10.7$ and an impact-duration coefficient of $C_{\mathcal {T}} = 2.142$.

The second front in the wave packet arrives much later. It is associated with a much smaller peak coefficient of $\hat {C}^{2nd}_{HU^2}= 5.17$ that occurred at time $t^*_j = {32.7}$.

6. Dependence on Froude number and obstacle's location

We conducted simulations of the impact on a circular prism to find the impact coefficients on two dominant influencing parameters. Figure 7 shows the dependence of (a) the peak wave force coefficient $\hat {C}_{HU^2}$, and (b) the impact-duration coefficient ${ S_o{\mathcal {T}} U}/{H}$ on the undisturbed-flow Froude number $Fr$ and the location of the obstacle $S_ox_b/H$. The colours of the lines and symbols black, red and blue in the figure define the Froude numbers $Fr = 3.71$, 4.63 and 5.60, respectively. Four series of simulations designated as cases 1, 2, 5 and 6 are conducted. The tables embedded in the figure list the condition of the simulations.

Figure 7. Dependence of (a) the peak wave force coefficients $\hat {C}_{HU^2}$ and (b) the impact-duration coefficients $S_o{\mathcal {T}}U/H$ for the circular prism on two parameters: $Fr$ and $S_ox_b/H$. The black, red and blue colours refer to $Fr = 3.71$, 4.63 and 5.60, respectively. The green dashed lines in (a) and (b) represent the asymptotic limits. (c) The tables list the conditions of the simulations for the effects of the perturbation period $S_oTU/H$ and the friction coefficient $c_f$.

6.1. Reference case 1

The solid lines in figure 7 represent the simulation results of reference case 1 for $Fr = 3.71$, 4.63 and 5.60. In this reference case 1, the inlet perturbation has an amplitude of $\epsilon = 0.2$ and the period of $T_1^*=S_oT_1U/H=6.08$, 13.4 and 19.8 for $Fr = 3.71$, 4.63 and 5.60, respectively. The friction coefficients and channel slopes for the reference case 1 are $(c_f, S_o)=(0.00728, 0.0501)$, (0.00786, 0.0843) and (0.00760, 0.119) for $Fr = 3.71$, 4.63 and 5.60, respectively. The friction coefficients in this reference case are the same as the smooth channel bed in the laboratory experiments by Brock (1967). The channel slope $S_o$ is selected by Brock (1967) to satisfy (2.4) for balancing channel slope and friction in the undisturbed uniform flow. The reference periods $T_1$ are the same as the periods in the experiment by Brock (1967) who produced the wave train in the laboratory for comparison with Dressler's analytical solution.

6.2. Cases 5 and 6 for the effect of inlet perturbation period

The simulations of cases 5 and 6 are conducted to examine the effect of the inlet perturbation period. For simulation case 5, the perturbation period is $S_oTU/H=2 T_1^*$, two times longer than the reference period. For simulation case 6, the period is $S_oTU/H=\frac {1}{2} T_1^*$, which is one half of the reference period. We find little change in the impact force when the perturbation period, $T$ in (3.1), changes by a factor of four. The triangle symbols, representing the simulation results for cases 5 and 6, closely follow the solid lines for reference case 1 in figure 7.

6.3. Effect of bed friction

The channel slope and rough channel bed friction coefficient for simulation case 2 are $S_o= 0.100$ and $c_f= 0.0146$; they are greater than the reference slope and friction of case 1 by a factor of two. Figure 7 shows that the results of case 2 (denoted by the circle symbols) are essentially the same as reference case 1 (represented by the solid lines). Therefore, friction relatively is a minor effect that does not change much the dependence of the impact on the two dominant parameters.

The friction effect may be negligible on a channel of a moderate slope but we cannot completely rule out a friction dependence. The additional simulations presented in Appendix D show the force coefficient can be as much as 16 % higher if the friction and the channel slope are increased from the reference value by a factor of ten. The gravity components would be $g' = g \cos \theta =8.74$ m s$^{-2}$ and $g'S_o = g \sin \theta =4.91$ m s$^{-2}$ on such a steep slope of angle $\theta = 27^\circ$. The higher impact force coefficient for this case of extreme friction is due to the significant change in the gravity components. However, in reality, such extreme friction coefficient $c_f = 0.0728$ – ten times greater than the value of a smooth channel bed – is not likely for a rough bed on a natural channel.

6.4. Effect of obstacle size

Appendix D also examines the dependence of the wave force on obstacle size. Changing the relative size of the obstacle from $D/H= 75$ to 1 would lead to an 8 % reduction in the force coefficient, according to the calculation in the appendix. A similar minor dependence of the force coefficient on the obstacle size was observed in the study of the tsunami impact force on the coastal structure by Xie & Chu (2019, 2020).

6.5. Asymptotic limits

The front-runner amplitude increases toward the long-wavelength limit of Dressler's solution, which, as indicated by the green dashed lines in figure 2, are $\hat {h}_{FR}/H= 4.20$, 6.48 and 9.25 and $c_{FR}/U= 1.89$, 2.08 and 2.26 for $Fr = 3.71$, 4.63 and 5.60, respectively. The corresponding asymptotic limits of the wave force coefficient and the impact-duration coefficient, represented by the green dashed lines in figure 7, are ${\hat {C}}_{HU^2} \simeq 22.3$, 46.9 and 88.4 and $S_o{\mathcal {T}}U/H \simeq 3.0$, 3.8 and 4.1, for $Fr = 3.71$, 4.63 and 5.60, respectively.

To determine the asymptotic limits, we have conducted 2-D numerical simulation following the wave packet over a long distance as far as $S_ox_b/H= 4000$. Such simulation is computationally demanding because a vast computational domain is needed to fully accommodate the wave packet's spatial development to the nonlinear stage while maintaining the numerical stability of the computation.

6.6. Fast-moving instabilities

The wave force coefficients in figure 7 are significant because the roll waves are fast-moving instabilities. The front runner not only has a considerable wave height, it also has a large velocity occurring simultaneously. The combined effect of the wave height and speed of the front runner produces the large momentum flux ${h}{u}^2$ as demonstrated in Appendix C. For ${Fr}=5.60$, the relative momentum flux to the undisturbed flow $hu^2/(HU^2)$ of the front runner could reach 39.73 at the asymptotic limit. The momentum flux is so significant that we must consider the roll waves when evaluating the loading force on a structure. Determining the structural loading based on the undisturbed flow would significantly underestimate the impact.

6.7. Wave drag vs viscous drag

We ignored viscous drag in calculating the wave force using (4.2). In the absence of waves, the drag is due to boundary layer development and the formation of a turbulence wake. On the other hand, the wave force – due to the bow shock wave and the water run-up on the front face of the structure – can be more than an order of magnitude greater than the viscous drag on a submerged body.

In a textbook such as Potter, Wiggert & Ramadan (2017), the drag coefficient for turbulent flow around a fully submerged circular cylinder is less than one. This is to compare with the asymptotic limits of the wave force coefficients ${\hat {C}}_{HU^2} \simeq 22.3$, 46.9 and 88.4 for $Fr = 3.71$, 4.63 and 5.60, respectively, represented by the green dashed lines in figure 7.

7. Shapes and orientations of the structures

We used the same numerical scheme to calculate the wave impact force on structures of different shapes and orientations. The obstacles examined are (a) square prism (SP), (b) circular prism (CP), (c) rotated-square prism (RP) and (d) triangular prism (TP). The base of the TP is an equilateral triangle. The depth-contour maps of $h/H$ in figure 8 show the patterns of the scattering waves around these obstacles on the $x$-$y$ plane at different times.

Figure 8. Details of the contour maps of $h/H$ on the $x$–$y$ plane for the wave impact on (a) square prism, (b) circular prism, (c) rotated-square prism and (d) triangular prism located at $x_b/D=266.67$ ($S_ox_b/H=1004.7$), at times $S_o(t-t_{arr})tU/H = -2.68$, 0, $+$2.68. The simulations were conducted for the dimensionless width of the obstacles $D/H =75$ and the undisturbed Froude number ${Fr} = 3.71$.

The stand-off distance $x_s$ of the bow shock wave and the run-up height $R$ depend on the obstacles’ shape. We found a stand-off distance particularly small for structures with a sharp and pointy front, such as the RP and TP. The small stand-off distance leads to a small run-up height and a small pile-up of water in front of the obstacle. Therefore, the wave impact force on the TP can be significantly smaller compared with the impact on a SP with the same dimensionless width $D/H$.

Figure 9(a–c) shows the dependence of the peak wave force coefficient $\hat {C}_{HU^2}$ on the two parameters $Fr$ and $S_o x_b/H$ for the obstacles of the four different cross-sections. All things being equal, the highest impact force is on the SP. The CP is the second, and the RP is the third in the wave impact force ranking. Finally, the least impact force is on the TP.

Figure 9. The dependence of the peak wave force coefficients $\hat {C}_{HU^2}$ (a–c) and the impact-duration coefficient $S_o{\mathcal {T}} U/H$ (d–f), on $Fr$ and $S_o x_b/H$ for obstacles of different shapes and orientations. The symbols $\square$, $\circ$, $\diamond$ and $\triangleleft$ refer to the results obtained for the SP, CP, RP and TP, respectively.

Depending on the obstacle's shapes and orientations, the wave force coefficient can be as high as $\hat {C}_{HU^2} \simeq 46.2$ for the front runner with ${Fr} = 3.71$ impacting on SP and as low as $\hat {C}_{HU^2} \simeq 8.12$ for the impact on TP located at $S_ox_b/H = 1000$. At the same location $S_ox_b/H = 1000$, the force coefficients on SP and TP are $\hat {C}_{HU^2} \simeq 84.0$ and 12.1 for ${Fr} = 4.63$, and $\hat {C}_{HU^2} \simeq 145$ and 15.6 for ${Fr} = 5.60$, respectively. The peak wave force coefficients $\hat {C}_{HU^2}$ of all obstacles at $S_ox_b/H = 1000$ are more than an order of magnitude greater than the textbook values of the drag coefficients of turbulent flow on obstacles.

Figure 9(d–f) shows impact-duration coefficients $C_{\mathcal {T}}={S_o{\mathcal {T}} U}/{H}$. The longest duration coefficient is associated with the impact on TP, while the SP has the shortest impact duration. Table 1 summarizes the numerical values for the peak wave force coefficient $\hat {C}_{HU^2}$ and the impact-duration coefficient $C_{\mathcal {T}}$ and the product of the two coefficients $\hat {C}_{HU^2} C_{\mathcal {T}}$ for the obstacles SP, CP, RP and TP located at $S_ox_b/H=1000$.

Table 1. The peak wave force coefficient $\hat {C}_{HU^2}$, impact-duration coefficient $C_{\mathcal {T}}$ and the product of $\hat {C}_{HU^2}$ and $C_{\mathcal {T}}$ for SP, CP, RP and TP located at $S_ox_b/H=1000$.

7.1. Direction cosine of the oblique surface

The orientation of the obstacle surface is provided by the direction cosine, $\cos \phi = - \boldsymbol {{i}\,{\cdot }\,{n}}$, defined by the unit vector $\boldsymbol {{i}}$ in the direction of the incoming waves and the outward unit vector $\boldsymbol {{n}}$ perpendicular to the oblique surface. As shown in figure 10, the direction cosines are (a) $-\boldsymbol {i} \boldsymbol {\cdot } \boldsymbol {n}=\cos 45^\circ$ and (b) $-\boldsymbol {i} \boldsymbol {\cdot } \boldsymbol {n}=\cos 60^\circ$ for the oblique surfaces of the RP and the TP, respectively.

Figure 10. The incident waves and reflected waves on the oblique surfaces of (a) the RP and (b) TP. The direction cosines are (a) $-\boldsymbol {i} \boldsymbol {\cdot } \boldsymbol {n}=\cos 45^\circ$ and (b) $-\boldsymbol {i} \boldsymbol {\cdot } \boldsymbol {n}=\cos 60^\circ$ for the oblique surfaces of the RP and the TP, respectively.

We propose the following heuristic interpretation for impact reduction. The pressure force per unit length $f_{hs}$ in (4.2) is assumed proportional to the normal component of the incident and reflected waves. If the incident wave momentum flux is $\dot {m} \, \boldsymbol {i}$, the contribution of the normal component to the pressure on the frontal surface of the obstacle would be $- \dot {m} \, \boldsymbol {i} \boldsymbol {\cdot } \boldsymbol {n}=\dot {m} \cos \phi$. That is, $f_{hs} \propto \dot {m} [\cos \phi ]$, because the pressure is proportional to the component of the incident and reflected waves perpendicular to the oblique surface. Then the force $F_w$ is the component of the pressure force in the direction of the wave. The result is the following approximation for the wave force:

(7.1)\begin{equation} F_w ={-} \oint_{s} f_{hs} \, \boldsymbol{i} \boldsymbol{\cdot} \boldsymbol{n} \, {\rm d}s \propto \oint_{s} \dot{m} \, [\cos \phi]^2 \,{\rm d}s. \end{equation}

As shown in figure 8, the back end of the obstacle remains dry. Therefore, the contribution of the rear surface to the wave force is negligible.

On the oblique surface of the RP, the direction-cosine square is $[\cos 45^\circ ]^2= 0.5$. Therefore, the roll-wave impact force on the RP prism would be 50 % of the impact force on the SP if the cosine square were the only factor. On the other hand for the oblique surface of the TP, $[\cos 60^\circ ]^2= 0.25$. The wave impact force on the TP would be 25 % of the force on the SP.

This heuristic interpretation does not account for the different lengths of the line integration between RP and TP. It also does not consider the possible difference in the pressure distribution over the oblique surfaces. Nevertheless, the cosine-square correction has brought the data of the impact-force coefficients for the RS prism denoted by the diamond symbols to closely follow the data for the TP represented by the triangular symbols as shown in figure 11.

Figure 11. The relative wave force coefficient ${\hat {C}_{HU^2}}/ \{ [\hat {C}_{HU^2}]_{SP} [\cos \phi ]^2 \}$ for the RP and TP, corrected by the direction-cosine square. The dependence of the relative coefficient on (a) the position $S_o x_b/H$ and undisturbed-flow Froude number $Fr$, and (b) the amplitude of the front runner $\hat {h}_{FR}/H$. The symbols $\diamond$ and $\triangleleft$ denote the data for RP and TP, respectively.

7.2. Impact by front runner of finite amplitude

A similar radiation pressure dependence on the cosine square is known for electromagnetic waves impinging on oblique plane surfaces (Wright 1992). Unlike electromagnetic waves, the roll-wave dynamics is nonlinear. Hence, the reduction is further dependent on the amplitude of the roll waves. Figure 11(a) shows the reduction in terms of the relative coefficient that depends on the obstacle's position and undisturbed-flow Froude number, and figure 11(b) on the amplitude of the front runner, $\hat {h}_{FR}/H$. Note that, as shown in figure 11(b), the ratio ${\hat {C}_{HU^2}}/{ [\hat {C}_{HU^2}]_{SP} [\cos \phi ]^2} \simeq 1$ as front runner's amplitude $\hat {h}_{FR}/H$ approaches unity for small-amplitude roll waves. Therefore, the cosine square will be the only reduction to the force on the oblique surface if the waves’ amplitude is sufficiently small.

We may use the relationship in figure 11(b) to find the wave impact on the structure with a sharp front pointing to the incoming waves of finite amplitude. As an example, for a front runner with an amplitude of $\hat {h}_{FR}/H= 2$, the relative coefficient varies in the range

(7.2)\begin{equation} \frac{\hat{C}_{HU^2}}{ [\hat{C}_{HU^2}]_{SP} [\cos \phi]^2} \simeq 0.72 \ {\rm to} \ 0.83 \end{equation}

in the range of Froude number varying from $Fr = 5.60$ to 3.71. The selection of a TP would have $[\cos \phi ]^2 = 0.25$ and hence a relative coefficient $[\hat {C}_{HU^2}]_{TP} / [\hat {C}_{HU^2}]_{SP} \simeq 0.18$ to 0.21. This calculation shows the effect of the front runner's amplitude and the structure's shape and orientation. The wave force on the TP pointing at the incoming wave is $[\hat {C}_{HU^2}]_{SP} / [\hat {C}_{HU^2}]_{TP} \simeq 4.8$ to 5.6 times smaller compared with the impact on a blunt object of the SP having the same width in this example with a front-runner amplitude of $\hat {h}_{FR}/H= 2$.

8. Summary and conclusion

In this paper, we studied the spatial development of roll waves produced by a local disturbance focusing on the impact force of the waves against obstacles. Numerical simulations were conducted using the shallow-water model and an adaptive quadtree mesh to determine the evolution of the roll-wave packet and the impact force against obstacles of various shapes and orientations. We found the wave impact force raised sharply to a peak on the arrival of the front runner of the wave packet.

The peak force coefficient depends on two dominant influencing parameters: the Froude number of the undisturbed flow $Fr$ and the dimensionless obstacle distance from the local disturbance $S_o x_b/H$. The front-runner amplitude increases toward the long-wavelength limit of Dressler's solution, which are $\hat {h}_{FR}/H= 4.20$, 6.48 and 9.25 and $c_{FR}/U= 1.89$, 2.08 and 2.26 for $Fr = 3.71$, 4.63 and 5.60, respectively. The corresponding asymptotic limits of the wave force coefficient for the impact on the CP are ${\hat {C}}_{HU^2} \simeq 22.3$, 46.9 and 88.4 for $Fr = 3.71$, 4.63 and 5.60, respectively. These wave force coefficients are much greater than the steady force coefficient without the roll waves. Therefore, determining the structural loading from the steady uniform flow and ignoring the roll-wave instability would significantly underestimate the impact.

The roll waves are destructive for the waves’ large momentum flux compared with the uniform flow. The wave impact against blunt objects like the SP and CP is most significant. On the other hand, pointy objects such as the RP and TP can re-direct the incident waves to reduce the impact force. For example, the rotation of a SP by 45$^\circ$ reduces the impact force coefficient by more than a factor of two. We suggest using a pointy front with a sharp angle to minimize the impact and propose a heuristic interpretation to evaluate the impact reduction. The impact force reduces proportionately to the direction-cosine square $[\cos \phi ]^2$ of the oblique surface for small-amplitude waves. A further reduction is related to the finite amplitude of the roll waves.

The simulation for the impact force is computationally demanding because a large computation domain is required to fully accommodate the roll-wave packet's development to the nonlinear stage. Limited by the finite computational resources, we have ignored the viscous and turbulence stresses in our shallow-water 2-D modelling of the wave force. The present simulation is the first-ever numerical study of the roll-wave impact force on structures. However, our 2-D simulation results are comprehensive, covering a good range of conditions for various influencing factors, setting the stage for future studies of the 3-D aspects of the problem.