2.1. Test cases

The six test cases (table 1) discussed in the present study reproduce dam-break waves of non-Newtonian fluids with a high volume of release. The non-Newtonian fluids considered in the simulations are characterised by a power-law rheology with a consistency index, k, and power exponent, n. An additional test, Case 0, was conducted with a Newtonian fluid (water). In the following, k for Case 0 is equal to the water dynamic molecular viscosity.

Table 1. Main parameters of the test cases.

The computation domain in the simulations is L = $110h_{0}$ long, B = 1.5 $h_{0}$ wide and H = 1.5 $h_{0}$ high, where $h_{0}$ is the initial height of the lock region containing the wave fluid. The lock gate is situated at $x/h_{0}=0$ (figure 1). The length of the region containing the lock fluid at the start of the simulations is $60h_{0}$ , with the floodable area extending to the right of the lock gate over a length of 50 $h_{0}$ . This configuration allowed for the discharge to be maintained close to constant at the lock gate during the simulations. The spanwise width of the channel was sufficiently large for the width-averaged solution to be independent of the domain width and to capture the formation of a fairly large number of near-wall streaks.

Figure 1. Sketch of the computational domain and its dimensions. View of the channel showing the initial location and depth, $h_{0}$ , of the wave fluid and the boundary conditions.

As summarised in table 1, in the clear water case (Case 0), the density and dynamic molecular viscosity were $ \rho=1000 \,\textrm{kg m}^{-3}$ and $ k=0.0009\,\textrm{Pa s}$ , respectively. The non-Newtonian power-law fluids in the other six cases had a density $\rho = 1064 \,\textrm{kg m}^{-3}$ . Cases 1–4 considered power-law fluids characterised by four different rheological index values (n = 0.1, 0.3, 0.6, 0.9) and the same consistency index that was 100 times higher with respect to the water dynamic molecular viscosity (Case 0). Comparison of Cases 1, 2, 3 and 4 allows for studying the effect of varying the power-law index on the structure and propagation of the dam-break wave. Cases 5 and 6 have the same n value as for Case 2, but the k value is 10 and 100 times lower for Case 5 and Case 6, respectively, compared with Case 2. Comparison of Cases 2, 5 and 6 enables understanding the effect of varying the consistency index.

Results are presented in dimensionless form, considering the following length, time and velocity reference scales: $h_{0}$ , $t_{0}=\sqrt{h_{0}/g}$ and $U_{0}=\sqrt{gh_{0}}$ , respectively, where g is the gravitational acceleration.

The reference Reynolds number is defined as

(2.1) \begin{equation}{\textit{Re}}=\frac{\rho U_{0}h_{0}}{{\rm \mu} _{0}}=\frac{\rho U_{0}^{2-n}h_{0}^{n}}{k},\end{equation}

where ${\rm \mu} _{0}=k({U_{0}}/{h_{0}})^{n-1}$ denotes the reference dynamic molecular viscosity for power-law fluids. In the clear water case, ${\rm \mu} _{0}$ denotes the (water) dynamic viscosity. For uniform conditions in a wide channel, i.e. two-dimensional steady parallel flow conditions, a Reynolds number ${\textit{Re}}'= ({\rho V_{u}D})/({{\rm \mu} _{0}})$ can be defined with the hydraulics diameter D $=4h_{u}$ , where $h_{u}$ is the uniform flow depth and $V_{u}$ is the corresponding depth-averaged streamwise velocity. The critical Reynolds number value for which the flow remains laminar is (Mishra & Tripathi Reference Mishra and Tripathi1971)

(2.2) \begin{equation}{\textit{Re}}^{\prime}_{c}=2100\frac{\left(4n+2\right)(5n+3)}{3\left(3n+1\right)^{2}}.\end{equation}

Although (2.2) is strictly valid for uniform flow conditions, the ${\textit{Re}}^{\prime}_{c}$ value is assumed to be fairly representative of critical conditions even for unsteady (e.g. dam break) flows. The values of ${\textit{Re}}$ and ${\textit{Re}}^{\prime}_{c}$ for each case are also included in table 1.

2.2. LES solver, grid dependency study and model validation

The numerical model is the one described by Del Gaudio et al. (Reference Del Gaudio, Constantinescu, Di Cristo, De Paola and Vacca2024) who performed similar LES of dam-break waves in more complex geometries. The governing equations are the (filtered) Navier–Stokes equations and the advection equation for the volume fraction of wave fluid, which is used to track the interface between the wave fluid and the surrounding air in the VOF–LES simulations. To maintain a sharp interface, a high-resolution interface capturing scheme is employed (HRIC). The filtered continuity and momentum equations are

(2.3) \begin{align}&\qquad\qquad\qquad\qquad\qquad\qquad \frac{\partial u_{\!j}}{\partial x_{\!j}}=0, \end{align}

(2.4) \begin{align}&\rho \frac{\partial u_{i}}{\partial t}+\rho \frac{\partial u_{i}u_{\!j}}{\partial x_{\!j}}=-\frac{\partial p}{\partial x_{i}}+\frac{\partial }{\partial x_{\!j}}\left[{\rm \mu} \left(\frac{\partial u_{i}}{\partial x_{\!j}}+\frac{\partial u_{\!j}}{\partial x_{i}}\right)\right]-g\delta _{i,3}-\frac{\partial \sigma _{\textit{ij}}}{\partial x_{\!j}}, \end{align}

where $p$ and $u_{i}$ represent the dimensional filtered pressure and Cartesian velocity component in the $i$ direction, respectively. In (2.4), $\delta _{i,k}$ denotes the Kronecker delta function, and $x_{1}$ , $x_{2}$ , $x_{3}$ correspond to the $x$ , $y$ and $z$ variables, respectively. Assuming a power law as the constitutive relation, the dynamic molecular viscosity is

(2.5) \begin{equation}{\rm \mu} =k \left(S_{\textit{ij}}S_{\textit{ij}}\right)^{\frac{n-1}{2}},\end{equation}

where $S_{\textit{ij}}$ is the filtered velocity deformation (rate of strain) tensor:

(2.6) \begin{equation}S_{\textit{ij}}=\frac{1}{2}\left(\frac{\partial u_{i}}{\partial x_{\!j}}+\frac{\partial u_{\!j}}{\partial x_{i}}\right)\!.\end{equation}

The molecular viscosity tends to infinity near the free surface for a power-law fluid. Following Gnambonde et al. (Reference Gnambode, Orlandi, Ould-Rouiss and Nicolas2015), we limited the Frobenius norm of the velocity deformation tensor to bound the molecular viscosity. A cut off value of 10⁻⁴ was used. The subgrid stress tensor, $\sigma _{\textit{ij}}$ , is calculated using the dynamic Smagorinsky model (Rodi, Constantinescu & Stoesser Reference Rodi, Constantinescu and Stoesser2013; Gnambode et al. Reference Gnambode, Orlandi, Ould-Rouiss and Nicolas2015). Its expression is σ _ij $=$ −2 ${\nu}_t S_{\textit{ij}}=$ −2C _d $\varDelta$ ²|S|S _ij , where $\nu$ _t is the turbulent subgrid scale (eddy) viscosity, $\varDelta$ is the local grid spacing, |S| is the magnitude of $S_{\textit{ij}}$ and C _d is the model parameter that is dynamically calculated using information available from the smallest resolved scales (Rodi et al. Reference Rodi, Constantinescu and Stoesser2013).

The VOF free surface tracking technique was successfully applied to simulate dam-break and flood-wave flows (e.g. see Horna-Munoz & Constantinescu Reference Horna-Munoz and Constantinescu2018, Reference Horna-Munoz and Constantinescu2020). A pure-advection equation is solved for the volume fraction of water ( $\eta$ ): $({\partial \eta }/{\partial t})+u_{\!j}({\partial \eta }/{\partial x_{\!j}})=0$ . In cells containing only water, $\eta = 1$ . In cells containing only air, $\eta = 0$ . A value of 0.5 was chosen for $\eta$ to visualise the interface and determine the front position in LES.

The SIMPLE (semi-implicit method for pressure linked equations) algorithm is used to integrate the discretised Navier–Stokes equations that are advanced in time using a semi-implicit, iterative method. The third-order MUSCL scheme is used to discretise the advective terms, the second-order central-differences scheme is used to discretise the diffusive and pressure gradients terms, and a second-order accurate scheme is used to discretise the time derivatives in the governing equations. A pressure outlet boundary condition with a zero volume fraction of the liquid phase and a pressure of $0\,\textrm{Pa}$ is specified at the top boundary. No-slip conditions (i.e. zero velocity) are imposed at the bottom boundary (horizontal channel bed) and at the left boundary. The flow is assumed to be periodic in the spanwise (y) direction. The initial velocity field is set to zero and the lock gate is removed instantaneously at $t=0$ . The time step is the same for all simulations and it corresponds to a Courant number close to 0.2 away from the channel bed.

A finer mesh was used in regions situated close to the channel bottom, such that the attached boundary layers and the velocity gradients were well resolved. The Cartesian-like, unstructured mesh was refined near this boundary using 15 prism layers. This allowed for placing the first point off the surface of the channel bottom at approximately three wall units for the highest Reynolds number simulation. No wall functions were used. The spanwise and streamwise sizes of a typical computational cell were of the order of 20–30 wall units. Two hundred grid points were used to discretise the domain along the spanwise direction. The total number of grid cells was close to 18 million for all the simulations.

A grid sensitivity analysis was conducted. As an example, results for Case 5 are discussed. Figure 2 shows the temporal evolution of the dimensionless front position, $ x_{f}/h_{0} $ , and the profile of the spanwise averaged velocity at $ x/h_{0} $ $=$ 12.5 and $ x_{f}/h_{0} $ $=$ 25 in wall coordinates (U ⁺ $=$ $U/u_{\tau}$ and z ⁺ $=$ $z u_{\tau}/\nu_{w}$ , where U(x, z) is the dimensional spanwise averaged velocity profile, $u_{\tau}(x)$ is the local bed friction velocity and ${\nu}_w(x)$ is the kinematic molecular viscosity at the channel bed). The temporal evolution of the front position (figure 2 a) is close to identical in the simulations conducted with grids containing 18 million cells and 25 million cells. Significant differences are observed between the velocity profiles in figure 2(b) corresponding to the simulation performed using a coarse mesh containing close to 9 million cells and the two other simulations conducted on meshes containing more than 18 million cells. The velocity profiles inside the dam-break wave and the front velocity in the two finer-mesh simulations are very close. So, one can conclude that simulations performed using meshes with at least 18 million cells are grid independent. The slope of the logarithmic region in the two finer-mesh simulations is larger than the one expected for a Newtonian fluid (e.g. water). This behaviour is consistent with previous findings by Gnambonde et al. (Reference Gnambode, Orlandi, Ould-Rouiss and Nicolas2015) for non-Newtonian shear-thinning fluids.

Figure 2. Grid dependency study for Case 5. (a) Non-dimensional temporal variation of the front position; (b) non-dimensional streamwise velocity profile at section $x/h_{0}$ $=$ 12.5 when the front is situated at $ x_{f}/h_{0} $ $=$ 25. The very fine mesh contains 25 million cells, the fine mesh contains 18 million cells and the coarse mesh contains 9 million cells.

The solver (STARCCM+) was extensively validated for turbulent flows involving Newtonian fluids (water). In particular, the solver was shown to accurately predict unsteady free surface flows in RANS-VOF simulations with a deformable free surface (e.g. Horna-Munoz & Constantinescu Reference Horna-Munoz and Constantinescu2018, Reference Horna-Munoz and Constantinescu2020). Del Gaudio et al. (Reference Del Gaudio, Constantinescu, Di Cristo, De Paola and Vacca2024) compared free surface deformation and wall pressures from 3-D LES of a water dam-break wave impacting a vertical wall with laboratory experimental data. Very good agreement between experimental data and LES predictions was observed for these variables. In the same study, the non-Newtonian fluid rheology solver was validated using data from the dam-break experiments conducted by Minussi & Maciel (Reference Minussi and Maciel2012). In terms of non-dimensional wall units, the level of mesh refinement in the present simulations is similar to that used by Del Gaudio et al. (Reference Del Gaudio, Constantinescu, Di Cristo, De Paola and Vacca2024) in their dam-break simulations conducted with Newtonian and non-Newtonian fluids.

2.3. Shallow water models

2.3.1. Governing equations

Considering a homogeneous layer of fluid flowing over a horizontal bed, without lateral inflow or outflow, assuming that spatial variations occur over scales larger than the flow depth, and that flow resistance due to the sidewalls and surface tension effects are negligible, the dimensional depth-averaged momentum and mass conservation equations are (Hogg & Pritchard Reference Hogg and Pritchard2004)

(2.9) \begin{align}&\frac{\partial V}{\partial t}+\beta \frac{\partial hV^{2}}{\partial x}+gh\frac{\partial h}{\partial x}+\frac{\tau _{b}}{\rho }=0, \end{align}

(2.10) \begin{align}&\qquad\quad \frac{\partial h}{\partial t}+\frac{\partial\!\left(Vh\right)}{\partial x}=0, \end{align}

where h is the flow depth, V is the depth-averaged velocity, $\beta$ is the shape factor and $\tau_{b}$ is the bottom shear stress. Independently of the fluid rheology and flow conditions, i.e. laminar or turbulent, $\beta$ and $\tau_{b}$ are evaluated assuming that the shape of vertical velocity profile corresponds to uniform conditions. For a power-law fluid and laminar flow conditions, the expressions for the shape factor and the bottom shear stress are (Ng & Mei Reference Ng and Mei1994)

(2.11) \begin{equation}\beta =2\frac{2n+1}{3n+2}\quad\tau _{b}=k\left(\frac{2n+1}{n}\frac{V}{h}\right)^{n}.\end{equation}

For clear water and turbulent flow conditions, a unitary value of the shape factor is assumed, and the bottom shear stress is expressed as

(2.12) \begin{equation}\tau _{b}=\rho C_{\!f}V^{2},\end{equation}

where C _f is the friction coefficient. For turbulent flow of a Newtonian fluid over a rough boundary and at high Reynolds number, a constant value is generally assumed for the friction coefficient (Hogg & Pritchard Reference Hogg and Pritchard2004). By contrast, for turbulent flow over a smooth bottom, C _f may be evaluated using the empirical Blasius formula (Blasius, Reference Blasius1913). In the case of a wide-open domain, i.e. assuming the hydraulic diameter is equal to four times the flow depth, it reads

(2.13) \begin{equation}C_{\!f}=\frac{0.316}{8}\left[\frac{{\rm \mu} _{0}}{4\rho Vh}\right]^{1/4}.\end{equation}

Similarly to the clear water case, to deal with power-law fluids in turbulent flow conditions, we assume a unitary value of the shape factor. Moreover, the bottom shear stress is evaluated using (2.12) with the friction coefficient given by the correlation of Dodge & Metzner (Reference Dodge and Metzner1959) modified for wide open channels:

(2.14) \begin{equation}C_{\!f}=\frac{a_{n}}{2}\left[\frac{k}{\rho }\frac{\left(6+\frac{2}{n}\right)^{n}}{{2^{2n+3}}V^{2-n}h^{n}}\right]^{{b_{n}}}\end{equation}

with

(2.15) \begin{equation}a_n=0.0665+0.01175n,\quad b_n=0.365 - 0.177n+0.062n^2.\end{equation}

2.3.2. Numerical solutions of the SWE

The SWE solution of the dam-break problem has been found by solving the nonlinear hyperbolic system of (2.9)–(2.10), assuming either laminar or turbulent flow conditions inside the dam-break wave. A second-order (in space and in time) finite-volume scheme was employed. The numerical fluxes were evaluated using the Harten–Lax–Van Leer (HLL) scheme (Harten, Lax & van Leer Reference Harten, Lax and van Leer1983). To guarantee the second-order spatial accuracy of the scheme, the values of the conserved variables at the two sides of the interface were estimated using a piecewise linear reconstruction scheme with a nonlinear min-mod limiter (LeVeque, Reference LeVeque, George and Berger2011). Time integration was performed using the second-order, two-step Runge–Kutta scheme (Gottlieb & Shu Reference Gottlieb and Shu1998). At the outlet, a nonlinear non-reflecting characteristic boundary condition was prescribed (Nycander, Hogg & Frankcombe Reference Nycander, Hogg and Frankcombe2008). Additional details about the numerical scheme can be found from Campomaggiore et al. (Reference Campomaggiore, Di Cristo, Iervolino and Vacca2016). In what follows, the numerical results of the SWE model corresponding to laminar and turbulent flow conditions inside the wave are denoted SM-LN and SM-TN, respectively (table 2).

Table 2. Acronyms of the different numerical and theoretical models.

2.3.3. Novel analytical solution for a turbulent dam-break wave of a power-law fluid

The theoretical modelling of dam-break flows in canonical configurations (e.g. wide horizontal channel, 2-D flow) based on analytically integrating the SWE equations (2.9)–(2.10) with some additional assumptions has been the object of several studies. The first approximate analytical solution, derived by Ritter (1982), neglects the bottom resistance and assumes $\beta$ = 1. Approximate analytical solutions for turbulent clear-water dam-break waves have been proposed (e.g. Hogg & Pritchard Reference Hogg and Pritchard2004; Chanson Reference Chanson2009). Moreover, an analytical solution for a laminar dam break wave of a power-law fluid was derived by Hogg & Pritchard (Reference Hogg and Pritchard2004). In the next sections, the analytical solution by Hogg & Pritchard (Reference Hogg and Pritchard2004) for a dam-break wave of a laminar power-law fluid will be identified with the label HM-LA.

In the present section, an approximate analytical solution for an instantaneous-release dam-break wave of turbulent power-law fluid propagating in a horizontal, wide channel is deduced. The new solution is based on the conceptual model of Whitham (Reference Whitham1955) and can be thought as the generalisation for turbulent power-law fluids of Chanson’s (Reference Chanson2009) solution valid for a turbulent dam-break wave of clear water. To this aim, considering the reference length, time and velocity scales previously indicated (i.e. $h_{0}$ , $t_{0}=\sqrt{h_{0}/g}$ and $U_{0}=\sqrt{gh_{0}}$ ), the following dimensionless quantities are defined:

(2.16) \begin{equation}\tilde{x}=\frac{x}{h_{0}} ;\quad \tilde{t}=\frac{t}{t_{0}} ;\quad \tilde{h}=\frac{h}{h_{0}} ;\quad \tilde{V}=\frac{V}{U_{0}}.\end{equation}

The SWE model for a turbulent power-law fluid, consisting of (2.9) (assuming $\beta$ = 1), (2.10), (2.12), (2.14) and (2.15) is rewritten in dimensionless form as follows:

(2.17) \begin{align}&\frac{\partial \tilde{V}}{\partial \tilde{t}}+\tilde{V}\frac{\partial \tilde{V}}{\partial \tilde{x}}+\frac{\partial \tilde{h}}{\partial \tilde{x}}+\tilde{C}_{\!f}\frac{\tilde{V}^{\lambda }}{\tilde{h}^{\phi }}= 0, \end{align}

(2.18) \begin{align}&\qquad \frac{\partial \tilde{h}}{\partial \tilde{t}}+\tilde{V}\frac{\partial \tilde{h}}{\partial \tilde{x}}+\tilde{h}\frac{\partial \tilde{V}}{\partial \tilde{x}}=0, \\[9pt] \nonumber \end{align}

in which

(2.19) \begin{equation}\tilde{C}_{\!f}=\frac{a_{n}}{2}\left[\frac{k}{\rho }\frac{\left(6+\frac{2}{n}\right)^{n}}{2^{2n+3} U_{0}^{2-n}h_{0}^{n}}\right]^{{b_{n}}} ;\quad \lambda =b_{n}\left(n-2\right)+2 ;\quad \phi =nb_{n}+1,\end{equation}

and the expressions of the coefficients a _n and b _n are given by (2.15).

Following Whitham (Reference Whitham1955), the turbulent dam-break wave of a power-law fluid is analysed by splitting the solution domain into an outer region, $-\infty \lt \tilde{x}\leq \tilde{x}_{1}$ , in which the friction is neglected (inviscid fluid), and an inner region, close to the wave front, ${\tilde{x}_{1}}\lt \tilde{x}\leq \tilde{x}_{s},$ where the effects of bottom friction are important. The streamwise lengths $\tilde{x}_{1}$ and $\tilde{x}_{s}$ denote the dimensionless position of the inviscid–turbulent interface and the wave front position, respectively. Both quantities are unknown and vary with time.

In the outer region, where bottom friction is neglected (inviscid fluid), the solution of (2.17) and (2.18) is the well-known Ritter’s solution (Ritter Reference Ritter1892):

(2.20a) \begin{align}\tilde{h}_{I}&=\frac{1}{9}\left(2-\frac{\tilde{x}}{\tilde{t}}\right)^{2}\qquad -\!\infty \lt \tilde{x}\leq \tilde{x}_{1}, \end{align}

(2.20b) \begin{align} \tilde{V}_{I}&=\frac{2}{3}\left(1+\frac{\tilde{x}}{\tilde{t}}\right)\!, \end{align}

with the wave front travelling with a dimensionless celerity $\tilde{U}_{f,I}$ equal to two. The subscript ‘I’ denotes the inviscid region. Considering that the velocity does not spatially vary inside the inner region, its value can be assumed to be equal to the (unknown) wave front velocity $\tilde{U}_{\!f}$ (Chanson Reference Chanson2009). Therefore, assuming that the acceleration and the inertial terms are negligible with respect to the flow resistance inside the inner region, (2.17) reduces to the diffusive equation:

(2.21) \begin{equation}\frac{\partial \tilde{h}_{D}}{\partial \tilde{x}}+\tilde{C}_{\!f}\frac{\tilde{U}_{\!f}^{\lambda }}{\tilde{h}_{D}^{\phi }}=0, \quad {\tilde{x}_{1}}\leq \tilde{x}\leq \tilde{x}_{s}.\end{equation}

Integrating (2.21) between $\tilde{x}_{1}$ and $\tilde{x}_{s}$ , and imposing the boundary condition $\tilde{h}_{D}(\tilde{x}_{s})=0,$ leads to the following inner region approximate solution:

(2.22a) \begin{equation}\tilde{h}_{D}=\left[\tilde{C}_{\!f}\left(\phi +1\right)\tilde{U}_{\!f}^{\lambda } \right]^{\frac{1}{1+\phi }}\left(\tilde{x}_{s}-\tilde{x}\right)^{\frac{1}{1+\phi }} {\tilde{x}_{1}}\leq \tilde{x}\leq \tilde{x}_{s}\end{equation}

assuming

(2.22b) \begin{equation}\tilde{V}_{D}=\tilde{U}_{\!f}.\end{equation}

The subscript ‘D’ is used for variables inside the turbulent region that are calculated by solving the diffusive (2.21). Imposing that the inner and outer domain solutions coincide at the interface between the two regions ( ${\tilde{h}_{I}}(\tilde{x}_{1})=\tilde{h}_{D}(\tilde{x}_{1})$ and ${\tilde{V}_{I}}(\tilde{x}_{1})=\tilde{V}_{D}(\tilde{x}_{1})=\tilde{U}_{\!f}$ ) leads to the following expressions for $\tilde{x}_{1}$ and $\tilde{x}_{s}$ :

(2.23a) \begin{align}&\qquad\qquad\qquad\qquad\qquad\quad \tilde{x}_{1}=\frac{1}{2}\big(3\tilde{U}_{\!f}-2\big)\tilde{t},\end{align}

(2.23b) \begin{align}&\tilde{x}_{s}=\frac{1}{2\tilde{C}_{\!f}\left(\phi +1\right)}\big[\tilde{U}_{\!f}^{-\lambda }2^{-2\phi -1}\big(2-\tilde{U}_{\!f}\big)^{2\phi +2}+\tilde{C}_{\!f}\big(\phi +1\big)\big(3\tilde{U}_{\!f}-2\big)\tilde{t}\big].\end{align}

Equations (2.23a ) and (2.23b ) allow for calculating $\tilde{x}_{1}(t)$ and $\tilde{x}_{s}$ (t), provided that the wave front velocity $\tilde{U}_{\!f}$ is known $.$ The latter is evaluated by imposing mass conservation, i.e. the mass of wave fluid in the inner region ( ${\tilde{x}_{1}}\lt \tilde{x}\leq \tilde{x}_{s}$ ) must be equal to the mass of ideal (inviscid) fluid in the corresponding front region $\tilde{x}_{1}\leq \tilde{x}\leq \tilde{U}_{f,I}\tilde{t}$ :

(2.24) \begin{equation}\int _{\tilde{x}_{1}}^{\tilde{x}_{s}}\tilde{h}_{D}\,{\rm d}x=\int _{\tilde{x}_{1}}^{\tilde{U}_{f,I}\tilde{t}}\tilde{h}_{I}\,{\rm d}x,\end{equation}

with $\tilde{U}_{f,I}=2$ .

Using (2.20a ), (2.22a ), (2.23a ) and (2.23b ), the following relationship between the dimensionless time and the wave front velocity $\tilde{U}_{\!f}$ is obtained:

(2.25) \begin{equation}\tilde{C}_{\!f}\left(2+\phi \right)\tilde{t}=\tilde{U}_{\!f}^{-\lambda }\left(1-\frac{\tilde{U}_{\!f}}{2}\right)^{2\phi +1}.\end{equation}

For each time $\tilde{t}$ , the iterative solution of (2.25) provides the value of $\tilde{U}_{\!f}$ . Then, $\tilde{x}_{s}$ and $\tilde{x}_{1}$ are calculated using (2.23b ) and (2.23a ), respectively. Once $\tilde{U}_{\!f},\tilde{x}_{s}$ and $\tilde{x}_{1}$ are known, the free-surface profiles can be calculated using (2.20a ) for $-\infty \lt \tilde{x}\leq \tilde{x}_{1}$ and (2.22a ) for ${\tilde{x}_{1}}\leq \tilde{x}\leq \tilde{x}_{s}$ .

The present analytical solution for turbulent clear water can be applied for analysing the wave propagating over a smooth wall by accounting for the friction coefficient variability along the current (Chanson Reference Chanson2009). Indeed, by adopting the Blasius formula (2.13) to estimate the friction coefficient, the analytical solution of a dam-break wave propagating over a smooth bed is given by (2.25), (2.20a ) and (2.22a ) assuming

(2.26) \begin{equation}\tilde{C}_{\!f}=\frac{0.316}{8}\left[\frac{{\rm \mu} _{0}}{4\rho U_{0}h_{0}}\right]^{1/4} \lambda =7/4 \quad \phi =5/4.\end{equation}

In the next sections, the SW analytical solution of a turbulent dam-break wave of clear water and power-law fluid propagating over a smooth bed will be identified with the label SM-TA.

As a final remark, we point out that because $0\leq \tilde{U}_{\!f}\leq 2$ , one can deduce two asymptotic expressions of the front velocity valid for high and respectively for low values of the wave front velocity. It is worth pointing out that high/low values of the wave front celerity correspond to small/high values of the time. Indeed, for high $\tilde{U}_{\!f}$ values (i.e. close to 2), the second term on the right-hand side of (2.25) tends to vanish, while the first term remains bounded. For such conditions, (2.25) may be approximated as follows:

(2.27) \begin{equation}\tilde{C}_{\!f}\left(2+\phi \right)\tilde{t} =\left(1-\frac{\tilde{U}_{\!f}}{2}\right)^{2\phi +1}.\end{equation}

The front velocity can then be calculated as

(2.28) \begin{equation}\tilde{U}_{\!f}=2-\left[\tilde{C}_{\!f}\left(2+\phi \right)\tilde{t}\right]^{\tfrac{1}{2\phi +1}}.\end{equation}

Meanwhile, for small $\tilde{U}_{\!f}$ . values (i.e. close to 0), the second term on the right-hand side of (2.25) remains finite, while the first term tends to diverge. For such conditions, the approximation of (2.25) becomes

(2.29) \begin{equation}\tilde{C}_{\!f}\left(2+\phi \right)\tilde{t} =\tilde{U}_{\!f}^{-\lambda },\end{equation}

which leads to the following expression of the wave front velocity:

(2.30) \begin{equation}\tilde{U}_{\!f}=\left(\frac{1}{\tilde{C}_{\!f}\left(2+\phi \right)\tilde{t}}\right)^{\tfrac{1}{\lambda }}.\end{equation}

To better illustrate this behaviour, one can consider the ratio between the two terms on the right-hand side of (2.25):

(2.31) \begin{equation}P=\frac{\left(1-\frac{U_{\!f}}{2}\right)}{U_{\!f}^{\lambda }}^{2\phi +1}.\end{equation}

The value of this ratio indicates which of the two terms is predominant and, based on this, which of the two approximate solutions given by (2.27) and (2.29) should be applied. Figure 3 reports the P ratio as a function of the wave front velocity, $\tilde{U}_{\!f}$ , for several values of the rheological index, n. Results indicate only a weak variability of P with n. Moreover, figure 3 can be used to estimate the applicability range of the two limiting expressions deduced for $\tilde{U}_{\!f}$ , i.e. (2.28) and (2.30). Assuming that a difference of two orders of magnitude allows neglecting one term with respect to the other one on the right-hand side of (2.25), the lower and upper bounds for the applicability of (2.28) and (2.30) are $ \tilde{U}_{\!f}^{\textit{min}}=1.4\div 1.5$ and $\tilde{U}_{\!f}^{\textit{max}}=0.03\div 0.06$ , respectively. In terms of the propagation time, the maximum time when the asymptotic solution given by (2.28) is valid, $\tilde{t}^{\textit{max}}$ , and the minimum time after which the asymptotic solution given by (2.30) is valid, $\tilde{t}^{\textit{min}}$ , can be found using (2.25) and setting $\tilde{U}_{\!f}=\tilde{U}_{\!f}^{\textit{min}}$ and $\tilde{U}_{\!f}=\tilde{U}_{\!f}^{\textit{max}}$ , respectively.

Figure 3. Variation of the P ratio (2.31) with the wave front velocity, $\tilde{\textrm{U}}_{{f}},$ for different values of the rheological index, n.

Compared with the SM-TN numerical solutions that are based on solving the full shallow water equations, the newly proposed analytical model shows significantly better agreement in terms of the predicted free-surface profiles compared with those given by the classical Ritter solution for non-Newtonian cases. The new analytical model slightly underpredicts the front speed given by the SM-TN solution. The effect is stronger for low values of n. The new analytical model correctly predicts the shape of the dam-break wave near its front for all n values. Moreover, for values of n between 0.3 and 0.9, the free surface profiles predicted by the new analytical model and by the corresponding SM-TN solution are in good quantitative agreement over the whole length of the wave.

3.1. Newtonian dam-break wave

Case 0 is used to illustrate the internal flow structure for a dam-break wave of Newtonian fluid (water) where the flow is strongly turbulent inside the wave. The streamwise velocity increases fairly monotonically inside the wave as the front is approached (figure 5 a), which means the turbulence will be stronger in regions situated closer to the wave front. Three-dimensional turbulent flow structures are present along the entire wave at the air–fluid interface. These structures are driven by the mean shear between the wave fluid and the surrounding air. As for the well-documented case of gravity currents (e.g. see Constantinescu Reference Constantinescu2014; Tokyay, Constantinescu & Meiburg Reference Tokyay, Constantinescu and Meiburg2012; Dai & Huang Reference Dai and Huang2022), the instabilities present near the front of the wave generate lobe and cleft structures (figures 5 b and 5 c). Such structures were already observed for high-Reynolds-number dam-break waves by Del Gaudio et al. (Reference Del Gaudio, Constantinescu, Di Cristo, De Paola and Vacca2024) for the case of waves propagating into a channel of limited width. The absence of sidewalls in the present simulations conducted with periodic boundary conditions in the spanwise direction does not have a significant effect on the development of these structures at the front of the wave.

Figure 5. Visualisation of the instantaneous flow structure of the dam-break wave for Case 0. (a) Non-dimensional streamwise velocity in the vertical plane $ y/B $ = 0.5; (b) streamwise non-dimensional velocity component in a horizontal plane situated at approximately 20 wall units from the channel bed; (c) near-wall streaks visualised using a Q iso-surface (Q = 500). The black arrow shows the position of the front of the wave ( $x_{\!f}/h_0 = 25$ ). The red arrows point towards the lobes and cleft structures at the front. The vertical red dashed lines show the locations of the velocity and turbulent shear stresses profiles in Figure 6.

Figure 6. Vertical profiles of the non-dimensional mean streamwise velocity and primary turbulent shear stress, $\overline{u^{\prime}w^{\prime}}$ , for Case 0. Results are shown when the front is situated at $ x_{f}/h_{0} $ = 28. (a) Streamwise (spanwise-averaged) velocity profile at $x/h_0 = 0$ (lock gate), $\kappa$ = 0.28; (b) streamwise (spanwise-averaged) velocity profile at $x/h_{0}$ = 8, $\kappa$ = 0.28; (c) streamwise (spanwise-averaged) velocity profile at $x/h_0 = 12.5$ , $\kappa$ = 0.40; (d) comparison of Case 0 velocity profile at $x/h_0 = 12.5$ with velocity profile measured by Wuthrich et al. (Reference Wüthrich, Pfister, Nistor and Schleiss2018). The blue dashed line corresponds to a power-law profile, $ U/U_{B} $ = ( $ h/H $ )^1/8; (e) $\overline{u^{\prime}w^{\prime}}$ profiles; ( f) $\overline{u^{\prime}w^{\prime}}$ profiles in wall coordinates. The horizontal dotted lines in panels (e) and ( f) correspond to the boundary between the logarithmic layer and the constant velocity layer of the dam-break wave.

Figures 5(b) and 5(c) also reveal the presence of near-bed coherent structures in the form of parallel streaks of higher and lower streamwise velocity near the channel bed, similar to those observed in turbulent boundary layers and fully developed turbulent flows in open channels. In the case of a turbulent dam-break flow, the streaks extend for some finite distance behind the front and gradually vanish at large distances from the front (e.g. see figures 5(b) and 5(c) where streaks are not present for $x/h_{0}$ $ \lt $ 5). This happens because, far from the wave front, the mean streamwise velocity inside the wave is not large enough to generate a sufficiently strong turbulent flow that can result in the formation of near-bed streaks. The average width of the streaks is close to 500 wall units (∼0.05B). So, these streaks are much larger in terms of their wall unit dimensions compared with those observed in fully turbulent boundary layers and channel flows. This is expected, as the flow conditions behind the front of the current are not exactly equivalent to those in the aforementioned statistically steady turbulent flows. Such large streaks were also observed behind the front of Newtonian high-Reynolds-number gravity currents in the large eddy simulations of Ooi, Constantinescu & Weber (Reference Ooi, Constantinescu and Weber2009) and Tokyay & Constantinescu (Reference Tokyay and Constantinescu2015). Given the reduced height of the wave close to the front, the streaks of high and low velocity affect the flow structure until close to the wave interface inside this region.

The flow structure inside turbulent dam-break waves is investigated based on the spanwise-averaged streamwise velocity profiles and primary turbulent shear stress profiles, $\overline{u'w'}$ , where the overbar denotes spanwise averaging and the velocity fluctuations in the streamwise and vertical directions are denoted $u^{\prime}$ and $w^{\prime}$ , respectively. A log-law region is present inside the turbulent dam-break wave in Case 0 (figure 6 a–6 c). The non-dimensional spanwise-averaged streamwise velocity profiles are plotted in inner coordinates at three streamwise locations extending from the lock gate to $x/h_0 = 12.5$ at the time when the front is situated at $x_{\!f}/h_0 = 25$ . The vertical coordinate in wall units is $z^+ = z^.u_{\tau}(x)/{\nu}_w(x)$ . In Case 0, $\nu_{w}(x)$ is a constant equal to the kinematic molecular viscosity of the water. Each velocity profile displays two main distinct regions: a boundary layer containing the wall viscous, buffer and log-law (sub) layers followed by a constant velocity layer that extends up to the wave interface with the air (see also discussion of figure 4). At the lock gate (figure 6 a), where turbulence is weak, the thickness of the viscous sublayer is large as it extends up to approximately 20 wall units away the bed surface. The viscous sublayer extends up to approximately 5 wall units at the other two streamwise locations where the turbulence is much stronger (figure 6 e) and velocity streaks are present. The main observation related to the presence of a log-law layer in Newtonian dam-break flows is that the slope of the logarithmic velocity profile (i.e. $1/\kappa$ ) decreases with the distance from the lock gate. The velocity profiles at $x/h_0 = 0$ and 8 suggest a value of $\kappa$ close to 0.28. However, in the near-front region of high streamwise velocity and strong turbulence, $\kappa$ approaches to the von Kármán constant value associated with steady turbulent flows (e.g. $\kappa$ ≈ 0.4 at $x/h_0 = 12.5$ in figure 6 c). In this region of strong turbulence, the log-layer extends over close to 40 % of the wave height, z_max (see also Case 0 velocity profile in figure 4).

Figure 6(d) shows that the predicted non-dimensional velocity profile at $x/h_0 = 12.5$ in Case 0 is in good qualitative agreement with that measured in a lab experiment conducted by Wuthrich et al. (Reference Wüthrich, Pfister, Nistor and Schleiss2018) for a turbulent dam-break wave propagating over a smooth, horizontal channel. Wuthrich et al. (Reference Wüthrich, Pfister, Nistor and Schleiss2018) showed that their streamwise velocity profiles were well approximated by a power-law function with an exponent of $ 1/8 $ , which also corresponds to a log-law region with a value of $\kappa$ close to 0.4. At larger distances from the bed, the velocity inside the wave was close to constant. The $ 1/8 $ power-law profile also provides a good approximation of the numerically predicted velocity profiles inside the log-law layer in the region where the flow inside the wave is strongly turbulent (e.g. see figure 6 c).

Near the lock gate, where turbulence is weak, turbulent shear stresses are found to be negligible (see profile of the primary turbulent shear stress $\overline{u'w'}$ at $x/h_0 = 0$ in figure 6 e). In the near-front region, where the turbulence inside the wave is strong (e.g. at $x/h_0 = 12.5$ in figure 6 f), the peak of $\overline{u'w'}$ occurs at approximately 100 wall units, which is comparable to the distance where the log-law layer starts of approximately 80 wall units (figure 6 c). This behaviour is consistent with what is observed in turbulent boundary layers and fully developed open channel flows. At higher elevations, the turbulent shear stresses decay monotonically with increasing distance from the smooth bed and reach a minimum value slightly above the interface between the log-law layer and the constant velocity layer. Interestingly, despite the fact that the mean streamwise velocity is close to constant, the shear stress inside the constant velocity layer increases as the wave–air interface is approached. This is consistent with the presence of smaller-scale deformations along this interface in figure 5(a).

3.2. Non-Newtonian dam-break wave

Case 5, which is characterised by a relatively high Reynolds number (Re = 7.2 $\times$ 10⁴, Re_MR = Re = 2.73 $\times$ 10⁵), is used to discuss the spanwise-averaged flow and turbulence structure inside non-Newtonian dam-break waves. Case 5 shows many similarities with Case 0 representing a Newtonian dam-break wave. For example, the friction-velocity Reynolds numbers at $x/h_0 = 12.5$ are very close in the two simulations (Re $\tau$ ≈ 500, where Re $_{\tau}$ = $u_{\tau} \delta/{\nu}_w$ and $\delta$ is the local thickness of the boundary layer).

Figures 7(a) and 7(b) show the presence of a log-law region inside the non-Newtonian dam-break wave. As for the case of a wave of Newtonian fluid, turbulence is very weak near the lock gate in Case 5, as illustrated by the very low values of the primary turbulent shear stress at $x/h_0 = 0$ in figure 7(c). The velocity profiles are qualitatively similar at $x/h_0 = 0$ in the Case 5 and Case 0 simulations with a thick viscous sublayer extending until $\textit{z}$ ⁺ ≈ 20, a relatively thin log-law sublayer where $\kappa$ ≈ 0.28, and a thick constant velocity layer extending until the air–water interface. It is relevant to point out that the velocity profiles near the lock-gate have a similar behaviour in the non-Newtonian and Newtonian cases even though the dynamic molecular viscosity, $\mu$ , varies by more than one order of magnitude in between the channel bed and the air–water interface (e.g. see profile of ${\mu}/{\mu}_w$ at $ x/h_{0} $ = 0 in figure 7(e), where $\mu_{w}$ is the dynamic molecular viscosity at the channel bed) in Case 5.

Figure 7. Vertical structure of the mean flow and primary turbulent shear stress, $\overline{u^{\prime}w^{\prime}}$ , for Case 5. Results are shown when the front is situated at $x_{\!f}/h_0 = 25$ . (a) Streamwise (spanwise-averaged) velocity profile at $x/h_0 = 0$ (lock gate), $\kappa$ = 0.28; (b) streamwise (spanwise-averaged) velocity profile at $x/h_0 = 12.5$ , $\kappa$ = 0.26; (c) $\overline{u^{\prime}w^{\prime}}$ profile at $x/h_0 = 12.5$ , (d) $\overline{u^{\prime}w^{\prime}}$ profile at $x/h_0 = 12.5$ in wall coordinates; (e) non-dimensional (spanwise-averaged) dynamic molecular viscosity profiles. The horizontal dotted lines in panels (c), (d) and (e) correspond to the boundary between the logarithmic layer and the constant velocity layer of the dam-break wave.

The predicted profiles at $ x/h_{0} $ = 12.5 when the front is situated at $ x/x_{f} $ = 25 in figures 6 and 7 are used to illustrate differences between the structure of non-Newtonian and Newtonian dam-break waves in the regions where the wave flow is strongly turbulent. In both Case 0 and Case 5 simulations, the viscous sublayer extends until $\textit{z}$ ⁺ ≈ 5. Though a well-developed log-law sublayer is present in both simulations, this layer starts at z ⁺ ≈ 80 in Case 0 and at z ⁺ ≈ 100 in Case 5. Also, the thickness of the log-law sublayer decays in the non-Newtonian case. For example, the interface with the constant velocity layer is situated around $ z/z_{\textit{max}} $ = 0.41 (z ⁺ ≈ 850) in Case 0 and $ z/z_{\textit{max}} $ = 0.32 (z ⁺ ≈ 1000) in Case 5. This effect is consistent with the observed behaviour of non-Newtonian shear-thinning fluids that are characterised by a flattened velocity profile due to the increase of the dynamic molecular viscosity away from the wall (e.g. see profile of $\mu/{\mu}_w$ at $x/h_0 = 12.5$ in figure 7 e).

While inside the regions of the dam-break wave generated in Case 0, where turbulence is strong the slope of the log law corresponds to $\kappa$ ≈ 0.4, the predicted value of $\kappa$ is close to 0.26 in the non-Newtonian dam-break wave generated in Case 5 (figure 7 b). This value is slightly lower than the value observed at $x/h_0 = 0$ in Case 5 (figure 7 a). The fact that the slope of the log law (i.e. $1/\kappa$ ) in regions where the flow is strongly turbulent is significantly higher for non-Newtonian dam-break waves compared with the value observed for Newtonian waves is not surprising and is consistent with the behaviour observed for simpler steady turbulent flows. For example, values of $\kappa$ close to 0.26 were reported for non-Newtonian fluids with a power law index n = 0.5 in steady turbulent pipe-flow experiments by Rudman et al. (Reference Rudman, Blackburn, Graham and Pullum2004) and in LES of the same type of flows by Gnambode et al. (Reference Gnambode, Orlandi, Ould-Rouiss and Nicolas2015).

In the region of strong turbulence inside the wave, the peak of $\overline{u'w'}$ is situated at approximately 100–200 wall units from the channel bed in Case 5 (figure 7 d), which is comparable to the value (∼100 wall units) predicted for the Newtonian dam-break wave in Case 0. Moreover, in both simulations, the primary shear stress decays monotonically inside the log-law sublayer and reaches its minimum slightly above the boundary between the log-law sublayer and the constant velocity layer. The values of $\overline{u'w'}$ are negligible inside the lower part of the constant velocity layer (figure 7 d), which was not the case for the Newtonian dam-break wave (figure 6 e). Though in both simulations, $\overline{u'w'}$ increases rapidly inside of the top part of the constant velocity layer ( $ z/z_{\textit{max}} $ > 0.8), the rate of increase is lower for the non-Newtonian dam-break wave due to the weaker instabilities developing along the air–fluid interface in this case.

4.1. Newtonian dam-break waves

Case 0 is used to illustrate the temporal evolution of the front position, the wave shape and the bed shear stress characteristics in turbulent dam-break waves of Newtonian fluids. It also serves as a limiting case to discuss dam-break waves of non-Newtonian fluids. The performance of turbulent SWE analytical (SM-TA) and numerical (SM-TN) models is discussed. For the sake of completeness, the comparison with the Ritter solution is also included.

Figure 8(a) compares the temporal evolution of the front position predicted by LES with those given by the SM-TA and SM-TN models and Ritter’s solution for Case 0. LES results show two distinct phases. During the initial acceleration phase, when the fluid inside the wave increases its velocity by releasing potential energy, the turbulence and viscous dissipation are negligible. The accelerating phase is followed by a deceleration phase, starting at approximately $ t/t_{0} $ = 8 (figure 8 b), during which bed friction effects are important. In contrast, the SW-TA and SW-TN models predict that the deceleration phase starts very soon after the release of the lock. This happens because the bottom turbulent shear stress is applied starting at $ t/t_{0} $ = 0.

Figure 8. Temporal evolution of the (a) non-dimensional front position and (b) front velocity for Case 0. Results of the shallow-water turbulent numerical model (SM-TN) and analytical solutions (SM-TA, Ritter) are compared with LES predictions.

During the deceleration phase, the temporal evolution of the front position predicted by LES can be described by a power-law relationship $ x_{f}/h_{0} $ = $ \gamma (t/t_{0})^{\alpha}$ . Since in practical applications one is generally interested in characterising the dam-break evolution at large times after the gate opening, the predictive capabilities of the shallow models are mainly assessed based on how well they can predict the values of $ \alpha $ and $ \gamma $ inferred from LES (i.e. $ \alpha $ = 0.885 and $ \gamma $ = 2.38). The SWE analytical and numerical models predict fairly close values of the power-law exponent (0.8 $ \lt $ $ \alpha $ $ \lt $ 0.87, see table 3).

Table 3. Power-law coefficients predicted by the different simulations and analytical solutions during the deceleration phase. Values predicted by the Hogg & Pritchard (Reference Hogg and Pritchard2004) solution (HM-LA) are also reported. Re _MR is the spanwise- and streamwise-averaged Metzner Reynolds number (Dodge & Metzner Reference Dodge and Metzner1959) calculated when $x_{\!f}/h_0 = 25$ .

Another important requirement for SWE numerical simulations is to accurately predict the sediment entrainment capacity in applications where the dam-break wave propagates over an erodible bed. Figure 9 compares the distributions of the non-dimensional bed shear stress predicted by LES and the SM-TN model in between the lock gate and the front when x_f = 25h ₀. LES predicts a close to linear decay of the bed shear stress from the lock gate until close to the front, with local variations due to the resolved turbulent eddies inside the wave. Except for the larger amplification of the bed shear stress very close to the front, the use of a turbulent-flow friction coefficient leads to very close agreement of the SM-TN model with LES.

Figure 9. Non-dimensional (spanwise-averaged) bed shear stress distribution for Case 0. The model predictions of the shallow-water turbulent numerical model (SM-TN) are compared with LES predictions when $x_{\!f}/h_0 = 25$ .

Given that the passage of a dam-break wave can induce severe flooding, it is important to also discuss how well SWE models predict the free surface profile. In addition to the SM-TN and SM-TA profiles, figure 10 also includes Ritter’s solution. Consistent with figure 8, the front in Ritter’s solution has travelled a considerably larger distance compared with the wave predicted by LES. However, the water depth in Ritter’s solution is very small in between $x/h_{0}$ ≈ 19 and the front so, overall, Ritter’s solution provides a good approximation to the free surface profile predicted by LES. Although the propagation length of the wave is underpredicted in the SM-TN simulation, the flow depth is reproduced quite accurately until the front region. Qualitatively, the analytical model SM-TA shows a similar behaviour, though it slightly overestimates the flow depth over part of the wave body compared with SM-TN. Overall, the use of a friction formula based on turbulent flows results in fairly accurate predictions of the water–air interface and bed shear stress distributions beneath the deeper parts of the wave.

Figure 10. Top (free surface) boundary of the dam break wave for Case 0. LES predictions are compared with the shallow-water turbulent numerical model (SM-TN) and analytical solutions (SM-TA, Ritter) at $ t/t_{0} $ = 11.2. The vertical arrows show the position of the front predicted by the different models.

4.2. Non-Newtonian dam-break waves: power-law index effect

As shown in figure 11, the temporal evolution of the non-dimensional front position in the LES (Case 1 to Case 4) conducted with a constant value of the consistency index (k = 0.09) is strongly influenced by n. Lower values of n, corresponding to a stronger shear-thinning behaviour of the fluid, result in a faster wave propagation after the end of the acceleration phase. These findings are consistent with the theoretical understanding of non-Newtonian rheology, where shear-thinning fluids experience reduced resistance under high-shear conditions. In particular, near the front, where strong shear is generated inside the wave, a faster propagation of the front with decreasing n is observed in the simulations. LES results show that the temporal evolution of the front position during the deceleration phase is well described by the power-law relationship discused for waves of Newtonian fluids and that $ \alpha $ increases monotonically with decreasing n (table 3). Moreover, for very small values of n, the power law exponent approaches a value of one ( $ \alpha $ = 0.968 for n = 0.1). Varying n while keeping the consistency index constant has only a small effect on the value of the power-law parameter, $ \gamma $ (table 3).

Figure 11. Effect of the power law index, n, on the temporal evolution of the non-dimensional front position. Vertical arrows show the onset of the deceleration phase.

In the following discussion, results of the SM-LN and SM-LA are included to highlight the inconsistencies that may arise from neglecting the presence of turbulence in SWE models. The use of a laminar friction coefficient in the SM-LN simulation leads to an increased velocity of the front compared with LES, while the opposite is observed for the SM-TN simulation using a turbulent friction coefficient (see results in figure 12 for Case 1 and Case 3). The predictions of the SM-TA solution using the turbulent friction coefficient are close to those given by the corresponding SM-TN numerical model. For all non-Newtonian dam-break cases, SM-TA is slightly more dissipative than SM-TN, which explains the faster propagation of the front predicted by SM-TN in figure 12(b).

Figure 12. Temporal evolution of the non-dimensional front position predicted by LES, the shallow-water laminar (SM-LN) and turbulent (SM-TN) numerical models, and by the analytical solution (SM-TA). (a) Case 1; (b) Case 3.

The implementation of the friction coefficient formula (2.14) for turbulent condition significantly improves the accuracy of the SWE model in predicting the propagation of the front especially for lower values of n. For example, results for Case 3 (n = 0.3) in figure 12(b) show that compared with LES, the overprediction of the front position by SM-LN is approximately 50 % larger than the underprediction of the front position by SM-TN. For higher values of n like Case 1 (n = 0.9) in figure 12(a), SWE model predictions are in good agrement with LES for both types of friction coefficients, with SM-TN and SM-TA showing slightly better agreement with LES than SM-LN.

The better performance of the SM-TN model to predict the front propagation is primarily due to ability of the Dodge and Metzner formula to more accurately capture the distribution of the friction coefficient beneath the turbulent dam break wave, as illustrated in figure 13 for Case 1 (n = 0.9). Conversely, the SM-LN model underestimates the resistance.

Figure 13. Friction coefficient distribution for Case 1. The LES predictions are compared with those of the SM-LN and SM-TN numerical models when x_f /h ₀ = 22.5.

For non-Newtonian fluid cases, the relative performance of the SWE models with respect to LES in predicting the free surface profile and the depth levels inside the wave is similar to that observed for Newtonian fluid cases (figure 10). For example, the SM-TN model gives the closest predictions of the depth levels inside the wave generated in Case 2 compared with LES (figure 14). The SM-LN model underestimates the wave depth over an important fraction of the wave’s length. Meanwhile, the SM-TA model overestimates the wave depth. Ritter’s solution significantly overestimates the propagation distance of the front. This is fully expected, given that there is no deceleration phase in Ritter’s solution ( $ \alpha $ = 1). However, it is interesting that the wave-depth levels predicted by Ritter’s solution away from the front region provide a good approximation to LES predictions, better than that shown by the SM-LN model.

Figure 14. Top (free surface) boundary of the dam break-wave for Case 2 shown at $ t/t_{0} $ = 14. The vertical arrows show the front position predicted by LES and the different SWE models.

The power-law relationship $x_{\!f}/h_0 = \gamma(t/t_0)^{\alpha}$ that provides a good approximation of the front evolution during the decelleration phase in the LES performed with non-Newtonian fluids can also be used to approximate the front evolution predicted by the SWE models. Figure 15 summarises the variation of $ \alpha $ with n obtained using LES, SWE numerical models, and the analytical outputs of the model proposed in this work (SM-TA) and that of Hogg & Pritchard (Reference Hogg and Pritchard2004) for laminar dam-break flows (HM-LA).

Figure 15. Power-law coefficient during the deceleration phase as a function of the power law index, n, for the different models.

SWE models that assume a turbulent friction coefficient (SM-TA and SM-TN) show a monotonically decreasing $ \alpha $ with increasing n, as also predicted by LES. The main difference is the rate of decrease, with the SM-TA model showing the smallest rate of change of $ \alpha $ with n. The largest differences in the predicted $ \alpha $ values by the three models are observed for low values of n. This is explained by the inability of SWE models to predict the correct value of the friction coefficient in very shear-thinning fluids. Both models that assume laminar condition, SM-LN and HM-LA, predict that $ \alpha $ decreases with n for low values of the power law index, while for $n \gt 0.3$ , an increase is observed.

4.3. Non-Newtonian dam-break waves: consistency index effect

The consistency index, k, characterises the fluid’s resistance to deformation. Lower values of k indicate reduced molecular viscosity and an increased contribution of inertia in non-Newtonian fluid flows. Figure 16 compares the front propagation in LES conducted with n = 0.6. Compared with Case 2, the value of k was 10 times smaller in Case 5 and 100 times smaller in Case 6, resulting in higher Reynolds numbers (table 1). As shown in figure 16, the main effect of decreasing k is an acceleration of the front after the end of the acceleration phase once inertia effects dominate inside the body of the turbulent wave. The highest value of $ \alpha $ ( $=$ 0.957) is observed for the highest Reynolds number simulation (Case 6, see table 3). This suggests that in the limit of a very small k, or infinite Reynolds number, $ \alpha $ will approach one, which is consistent with Ritter’s inviscid solution.

Figure 16. Effect of the consistency index, k, on the evolution of the non-dimensional front position. Vertical arrows show the onset of the deceleration phase.

Figure 17. Temporal evolution of the non-dimensional front position predicted by LES and some of the SWE numerical models and analytical solutions. (a) Case 2; (b) Case 6.

As observed in § 4.2, the use of a laminar friction coefficient results is a faster propagation of the front compared with LES, while the opposite is true for SWE models using a turbulent friction coefficient. For both Case 2 and Case 6, the SM-LN model consistently shows larger differences with LES compared with those shown by the SM-TN model (figure 17). Moreover, the predictions of the analytical model using a friction coefficient based on (2.14) are very close to those obtained using a full SWE numerical model, as can be seen from comparing SM-TA and SM-TN predictions in figure 17. An important finding is that for a constant n, the absolute diference between the predictions of the numerical and analytical SWE models and those of LES is increasing with increasing Reynolds number (i.e. with decreasing k), as can be seen by comparing results for Case 2 and Case 6 in figures 17(a) and 17(b), respectively.

Similar to LES, all SWE models predict a monotonic decrease of $ \alpha $ with increasing k, as shown in figure 18. The largest differences with LES are shown by the SM-LN model while, as indicated in table 3, the effect of k on $ \alpha $ is negligible for the HM-LA solution ( $ \alpha $ ≈ 0.82 for the cases with n = 0.6). Except for very high Reynolds numbers, the $ \alpha $ ₂ values predicted by the numerical and analytical SWE models using a turbulent friction coefficient are very close (i.e. less than 3 % difference) to the value obtained from LES. For both SM-TN and SM-TA, the difference increases to 6 %–7 % for k = 0.0009 (Pa s)^0.6. This value is significantly lower than the difference of approximately 14 % observed between the $ \alpha $ values predicted by LES and the SM-LN model for k = 0.0009 (Pa s)^0.6.

Figure 18. Power-law coefficient during the deceleration phase as a function of the consistency index, k, for the different models.