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 (Reference Jeffreys1925) and Dressler (Reference Dressler1949). Jeffreys (Reference Jeffreys1925) studied the linear stability of the steady-uniform flow down an inclined plane. Dressler (Reference Dressler1949) 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 (Reference Zanuttigh and Lamberti2002), Que & Xu (Reference Que and Xu2006) and Ivanova et al. (Reference Ivanova, Gavrilyuk, Nkonga and Richard2017) conducted numerical simulations for the nonlinear development of roll waves using the shallow-water equations. Richard & Gavrilyuk (Reference Richard and Gavrilyuk2012) introduced a modified shallow-water model to produce wave profiles in better agreement with the laboratory observation by Brock (Reference Brock1967).
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 (Reference Ng and Mei1994) determined the shortest wavelength corresponding to zero energy dissipation across the shock. Balmforth & Mandre (Reference Balmforth and Mandre2004) 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 (Reference Yu and Chu2022), 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.
The most noticeable feature of the wave packet is its front runner, which asymptotically approaches the long-wavelength limit of the Dressler (Reference Dressler1949) analytical solution according to Yu & Chu (Reference Yu and Chu2022). The front runner in a wave packet also was observed in a mudflow simulation by Liu & Mei (Reference Liu and Mei1994) and by Meza & Balakotaiah (Reference Meza and Balakotaiah2008), who studied the nonlinear waves on vertically falling films dominated by viscosity and surface tension. Meza & Balakotaiah (Reference Meza and Balakotaiah2008) 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 (Reference Yu and Chu2022) 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. (Reference Aureli, Dazzi, Maranzoni, Mignosa and Vacondio2015) 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 (Reference Xie and Chu2019) 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:
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
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:
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 (Reference Yu and Chu2022) 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 (Reference Yu and Chu2022). 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 Reference Criminale, Jackson and Joslin2003). 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 (Reference Dressler1949) solution. The green dashed lines in the figure indicate the asymptotic solutions.
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 (Reference Brock1967) on a smooth channel bed. For a given $c_f$ and $Fr$, Brock (Reference Brock1967) 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.
We used the Gerris Flow Solver by Popinet (Reference Popinet2003, Reference Popinet2011) 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. (Reference Liang, Zang, Borthwick and Taylor2007) and An & Yu (Reference An and Yu2012) simulated shock reflection by a circular cylinder. They found approximately 80 % saving of the CPU time compared with using uniform mesh. Popinet (Reference Popinet2015) 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.
The shock-capturing scheme for spatial discretization is a generalized MINMOD limiter (GML) which uses the limiter function
to reconstruct the cell-interface flux with a free parameter varied from $\beta =1$ to 2 (Sweby Reference Sweby1984; Nessyahu & Tadmor Reference Nessyahu and Tadmor1990; Kurganov & Tadmor Reference Kurganov and Tadmor2002). 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 (Reference Popinet2011).
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 Reference Gottlieb and Shu1998). The second-order accurate Cartesian cut-cell method of An & Yu (Reference An and Yu2012) and Causon, Ingram & Mingham (Reference Causon, Ingram and Mingham2001) 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:
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
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
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).
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. (Reference Stern, Wilson, Coleman and Paterson2001) to determine the accuracy of the simulation. The fractional error of the peak force coefficient is
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.
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 (Reference Forbes and Schwartz1981) 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.
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 (Reference Brock1967). The channel slope $S_o$ is selected by Brock (Reference Brock1967) 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 (Reference Brock1967) 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 (Reference Xie and Chu2019, Reference Xie and Chu2020).
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 (Reference Potter, Wiggert and Ramadan2017), 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.
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.
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$.
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.
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:
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.
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 Reference Wright1992). 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
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.
Acknowledgements
We thank the reviewers for their most helpful comments on an earlier version of our paper.
Funding
This research is supported by the Natural Sciences and Engineering Research Council of Canada under Grant RGPIN-2019-05776 and a MEDA Scholarship to Boyuan Yu.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Direct numerical simulation for Dressler's periodic solution
We reproduced the Dressler (Reference Dressler1949) analytical solution in this appendix by direct numerical simulation. Dressler (Reference Dressler1949) fitted the steep wavefronts between periodic smooth profiles. He obtained the wave profiles from the numerical integration of an ordinary differential equation (ODE) in a coordinate system moving with the celerity of the waves. We captured the shock waves between the smooth profiles from numerical simulation of the shallow-water equations, which are partial differential equations (PDE).
Figure 12 shows (a) the depth of the wave crest $\hat {h}/H$ and the depth of the wave trough $\check {h}/H$ at the steep wavefront, (b) the velocity $\hat {u}/U$ at the wavefront, and (c) the celerity $c/U$ of the periodic roll waves and their dependence on the wavelength $S_o \lambda /H$ of the perturbation. Note that the wave trains initiated by the periodic perturbation in space are not determined until the perturbation wavelength is specified.
The lines in the figure are Dressler's analytical solution obtained from integrating the ODE. On the other hand, the symbols are numerical solutions of the partial differential equations (2.1), (2.2) and (2.3). The solver of the PDE captured the steep wavefronts of Dressler's periodic roll waves shown in figure 12. It also accurately captures the front runners of all roll-wave packets considered in this paper.
Appendix B. Mesh refinement study
We conducted a mesh refinement study for the peak force coefficient $\hat {C}_{HU^2}$ obtained for the impact of roll waves against a CP with an undisturbed-flow Froude number of ${Fr}=3.71$. The CP of the diameter $D/H=75$ is at a position $S_ox_b/H=1005$. Figure 13(a) shows the transient variation of the force coefficient $\hat {C}_{HU^2}$ obtained for four levels of mesh refinements. The logarithmic plot in figure 13(b) shows the order of convergence obtained from mesh refinement calculations.
The dynamic quadtree adaptive mesh has a maximum mesh size $\varDelta _{max}$ and a minimum size $\varDelta _{min}$ over the computational domain of $L_x \times L_y= 162\ {\rm m} \times 3\ {\rm m}$. The numbers of computational cells over the channel length for the four refinement levels are $L_x/\varDelta _{min} = 6912$, 13 824, 27 648 and 55 296. Therefore, $r = {\Delta x_{k}}/{\Delta x_{k+1}} = 2$. The estimated values of the peak force coefficients are $\hat {C}_{HU^2}= 24.29$, 21.03, 20.56 and 20.36 for the four levels of refinements. Let three of the four estimates be $(\hat {C}_{k-1}, \hat {C}_{k}, \hat {C}_{k+1})$ for three sizes of the cell. The method of Stern et al. (Reference Stern, Wilson, Coleman and Paterson2001) determines the true value by extrapolation to zero cell size using the formula
where $r = {\Delta x_{k}}/{\Delta x_{k+1}}$ and $\hat {P}$ is the order of convergence determined by the formula
The true value $\hat {C}_{\Delta x \rightarrow 0}$ was determined by extrapolation using (B1) and (B2) to zero grid size. The error relative to this true value is $| { \hat {C}_{k} - \hat {C}_{\Delta x \rightarrow 0} } |$. In percentage, the fractional error is
For this example, the calculations in table 2 give an estimated true value of $\hat {C}_{HU^2}= {20.56}$, and the order of convergence $\hat {P}= 1.558$. The numbers of computational cells used in the present simulation are $L_x/\varDelta _{min}=27648$, $L_y/\varDelta _{min} =512$, $L_x/\varDelta _{max}=1728$ and $L_y/\varDelta _{max}= 32$. Therefore, the fractional error for the peak force coefficient for $Fr = 3.71$ and $S_ox_b/H=1005$ is $FE(\%)=1.175$%.
Appendix C. Momentum flux and alternative force coefficient
Roll waves are fast-moving instabilities. The momentum flux of the front runner is significant both for its amplitude and its speed. Figure 14 shows the momentum flux of the front runner relative to the momentum flux of the undisturbed flow ${h}{u}^2/(HU^2)$. For the undisturbed flow with $Fr = 5.60$, the relative momentum fluxes are $hu^2/(HU^2) = 14.85$, 25.41, 32.87 and 39.73 at position $S_o {x} /H = 400$, 800, 1200 and 1600, respectively.
The appropriate force coefficient is the one normalized by the front-runner kinetic energy $\frac {1}{2} \rho \hat {h} \hat {u}^2$, as follows:
Figure 15 shows that this alternative coefficient $\hat {C}_{hu^2}$ has a more moderate range of value, demonstrating that the peak momentum flux of the front runner $\rho \hat {h} \hat {u}^2$ is the correct scaling for the impact force.
The depth and velocity of the front runner $\hat {h}$ and $\hat {u}$ can be determined from the simulation without the obstacle. Then, the impact force $F_w$ is determined using the alternative force coefficient in figure 15. This way, we could estimate the peak impact force of the roll waves on various shapes and orientations of the obstacles without using a refined mesh to resolve the flow near the solid surface.
Appendix D. The minor effects of channel friction and obstacle size
We examine in this appendix the effects of bed friction and the obstacle size. But these are minor effects – relative to the dominant influences of the undisturbed Froude number $Fr$ and the obstacle position $S_o x_b/H$ studied in §§ 5 and 6.
First, let us examine further the bed friction effect for the undisturbed-flow Froude number of $Fr = 3.71$. Figure 16(a,b) shows the dependence of the peak force coefficient $\hat {C}_{HU^2}$ on friction coefficient. The solid symbols denote the reference case 1 for $c_f= 0.00728$ and $S_o=0.0501$.
When the friction coefficient is doubled changing from the reference value of $c_f = 0.00728$ for smooth channel bed to $c_f= 0.01456$, the percentage change to the wave force coefficient is a negligible 2 % and 1 % at positions $S_ox_b/H= 502.4$ and 1005, respectively. However, in the extreme when the friction coefficient increases ten times from the reference to $c_f=0.0728$, the increase of wave force coefficient is 16 % and 10 %, respectively. For the same Froude number $Fr = 3.71$, the channel slope also rose ten times from $S_o= 0.0501$ to 0.501 as required by (2.4).
The goal of this calculation is to show that friction is a minor effect. Note that Brock's laboratory channel with a slope $S_o=0.0501$ in reference case 1 was already on a steep slope. In reality, the extreme friction coefficient $c_f = 0.0728$ – ten times greater than the value of a smooth channel bed – is not likely a possibility for a rough natural channel.
The other minor effect is the obstacle's size. Figure 16(c,d) shows the peak force coefficient $\hat {C}_{HU^2}$ for the diameter-to-depth ratios $D/H= 1$, 5, 25 and 75, at locations $S_ox_b/H=502.4$ and 1005, respectively. The value of the coefficient almost reached its asymptotic values for reference case 1 when $D/H= 75$.