Large-eddy simulation of free-surface turbulent flow in a non-prismatic channel

Traditional studies of channel transitions take the empirical approach, which consider only 1D time-averaged free-surface flow (Chow 1959). They ignore 3D instantaneous turbulent features and fail to achieve optimal hydraulic efficiency and minimise erosion risks. The purpose of this study is to predict the turbulent features by means of two-phase large-eddy simulation (LES). The model domain is illustrated in Figure 1. The scope of work includes rigorous validations of LES predictions using newly available experimental data and benchmark solutions. These represent unique and new contributions. In the broad sense, this study demonstrates the two-phase LES approach as a modern approach to future hydraulic design, offering a good complement to the experimental approach.

Unlike the Reynolds-averaged equation models that simulate time-averaged flow field, LES provides instantaneous flow field. Knowledge about the instantaneous flow velocity and the associated peak strength is most relevant to addressing channel erosion risks. Previously, some researchers (e.g., Nashta & Garde 1988) obtained semi-analytical functions for describing transition profiles. Others (Murty Bhallamudi & Hanif Chaudhry 1992; Younus & Hanif Chaudhry 1994; Rahman & Hanif Chaudhry 1997) numerically simulated two-dimensional (2D) depth-averaged flow in expansion. These studies do not produce detailed structures of the flow field or separation features. This is one of the significant gaps addressed by this study.

Compared to other types of transitions (Chow 1959), the warped transition is a superior choice (USCAE 1991). Its sidewalls are flush with both the upstream and downstream channel sections (Li 2022). The seamless joints help reduce turbulent disturbances. A few recent studies (Choi et al. 2022; Zeng & Li 2022) compute instantaneous flows in a straight-wall transition rather than a warped transition.

We address the question of a realistic flow condition at the inlet (Figure 1), with required fluctuations in LES. Periodic conditions often used in LES studies (Xie et al. 2013; Nikora et al. 2019; Ramos et al. 2019) are invalid for a non-prismatic channel. However, inlet conditions may be generated using another method such as precursor method, recycling and rescaling method, or a synthetic method. The 1st method is reportedly the most accurate (Tabor & Baba-Ahmadi 2010), and the 2nd method is effective in producing a fully developed turbulent boundary layer, but both incur high computing costs and can possibly produce a periodicity (Dhamankar et al. 2018). The synthetic methods are a trade-off between computational accuracy and cost. The synthetic methods including a spectral synthesiser (SS) and a vortex method (VM) have been applied in recent LES studies of one-phase air flow (Antoniou et al. 2017).

We also address the question of the possible reflection of surface waves from the outlet (Figure 1) and the influence of initial conditions in LES on the flow development. The results from this study include 3D velocity, vortical motions and first- and second-order statistics of turbulence.

Few applications of SS and VM have been reported in the literature, and little is known about suitable values for parameters as input to the SS and VM. Previously, Liu et al. (2021) applied an SS in their study of flow in a rectangular vegetated channel, and Zeng & Li (2022) applied VM in their investigation of flow choking in a circular partially full pipe. The rectangular channel and circular pipe are prismatic channels. Non-prismatic channels are a new case, where more complicated mechanisms of flow separation, eddy formation and turbulence interaction are involved. This is one important difference of this study from the previous studies. Another important difference of this study is a systematic assessment of proper parameters in the SS and VM and their values. The contributions from this study include improved strategies for computing open channel turbulence and for visualising and analysing the computational results, with relevance to practical applications. The contributions will be discussed later in this paper.

We considered an incompressible flow of an air–water mixture, with air above the free surface and water below. The position of the free-surface was tracked using the volume of fluid method. The density, ⁠, of the mixture was calculated based on the volume fraction, ⁠, as ⁠, where the subscripts 1 and 2 refer to air and water, respectively. For a given cell, The viscosity, ⁠, of the mixture was calculated as ⁠.

We assess three ways to specify the inlet condition for the distribution of mean-inflow velocity ⁠:

• (1)

  1/7th power-law (Table 1, labelled as ‘1/7’);

• (2)

  1/14th power-law (labelled as ‘1/14’);

• (3)

  mapped 3D velocity from the laboratory measurements (labelled as ‘EXP’).

Table 1

Boundary condition, initial condition (IC) and hydraulic condition for LES runs

Run .	.	.	.	.	ITF .	IC .	SGS .	.	.	Nodes .	NBT .	Re .	Fr .	.	.
ID .	(m/s) .	(L/s) .	.	.	.	.	.	.	(m) .	.	.	.	.	.	.
T1	0.487	22.6	2.67	1/7th	VM	‘L’	DSM	5.13	0.2813	3763707	–	111230	0.33	–	–
T2					VM			7.7			–	111230		–	–
T3					VM			11.6			–	111230		–	–
T4					SS			5.13			–	111230		–	–
T5					SS			7.7			–	111230		–	–
T6					SS			11.6			–	111230		–	–
T7	0.487	22.6	2.67	1/7th	SS	‘L’	DSM	7.7	10a	3763707	–	111230	0.33	–	–
D1	0.511	23.7b	2.67	1/7th	VM	‘L’	DSM	7.7	0.2813	3763707	–	116712	0.34	0.74	0.76
D2											Yes	116712	0.34	0.72	0.70
D3	0.602	28.0	8c	1/7th	VM	‘L’	DSM	7.7	0.2813	4021768	–	137497	0.40	0.73	0.77
D4	0.487	22.6		1/14th							–	111231	0.33	0.54	0.57
D5	0.487	22.6	2.67	1/7th	SS	‘N’	DSM	7.7	0.2813	3763707	–	111231	0.33	–	–
D6						‘L’					–	111231	0.33	0.69	0.68
E1	0.481	22.3	2.67	EXP	VM	‘S’	DSM	7.7	0.2813	3763707	–	109860	0.32	0.78	0.77
E2					SS		WALE	2.5	0.2284		–	109860		0.85	0.82
E3					VM		WALE	2.5	0.2284		–	109860		0.82	0.81
E4					VM		WALE	2.5	0.2284		Yes	109860		0.80	0.78
E5					VM		WALE	7.7	0.2813		–	109860		0.77	0.77
E6					SS		WALE	7.7	0.2813		–	109860		0.84	0.81
E7					VM		DKE	7.7	0.2813		–	109860		0.78	0.78
E8					SS		DKE	7.7	0.2813		–	109860		0.87	0.81
E9	0.481	22.3	8c	EXP	VM	‘S’	DSM	7.7	0.2813	4021768	–	109860	0.32	0.85	0.84
P1	0.481	22.3	2.67	INSd	–	‘S’	DSM	–	–	3763707	–	109860	0.32	0.85	0.84
P2				INSe	–			–	–		–	109860	0.32	0.88	0.85

Run .	.	.	.	.	ITF .	IC .	SGS .	.	.	Nodes .	NBT .	Re .	Fr .	.	.
ID .	(m/s) .	(L/s) .	.	.	.	.	.	.	(m) .	.	.	.	.	.	.
T1	0.487	22.6	2.67	1/7th	VM	‘L’	DSM	5.13	0.2813	3763707	–	111230	0.33	–	–
T2					VM			7.7			–	111230		–	–
T3					VM			11.6			–	111230		–	–
T4					SS			5.13			–	111230		–	–
T5					SS			7.7			–	111230		–	–
T6					SS			11.6			–	111230		–	–
T7	0.487	22.6	2.67	1/7th	SS	‘L’	DSM	7.7	10a	3763707	–	111230	0.33	–	–
D1	0.511	23.7b	2.67	1/7th	VM	‘L’	DSM	7.7	0.2813	3763707	–	116712	0.34	0.74	0.76
D2											Yes	116712	0.34	0.72	0.70
D3	0.602	28.0	8c	1/7th	VM	‘L’	DSM	7.7	0.2813	4021768	–	137497	0.40	0.73	0.77
D4	0.487	22.6		1/14th							–	111231	0.33	0.54	0.57
D5	0.487	22.6	2.67	1/7th	SS	‘N’	DSM	7.7	0.2813	3763707	–	111231	0.33	–	–
D6						‘L’					–	111231	0.33	0.69	0.68
E1	0.481	22.3	2.67	EXP	VM	‘S’	DSM	7.7	0.2813	3763707	–	109860	0.32	0.78	0.77
E2					SS		WALE	2.5	0.2284		–	109860		0.85	0.82
E3					VM		WALE	2.5	0.2284		–	109860		0.82	0.81
E4					VM		WALE	2.5	0.2284		Yes	109860		0.80	0.78
E5					VM		WALE	7.7	0.2813		–	109860		0.77	0.77
E6					SS		WALE	7.7	0.2813		–	109860		0.84	0.81
E7					VM		DKE	7.7	0.2813		–	109860		0.78	0.78
E8					SS		DKE	7.7	0.2813		–	109860		0.87	0.81
E9	0.481	22.3	8c	EXP	VM	‘S’	DSM	7.7	0.2813	4021768	–	109860	0.32	0.85	0.84
P1	0.481	22.3	2.67	INSd	–	‘S’	DSM	–	–	3763707	–	109860	0.32	0.85	0.84
P2				INSe	–			–	–		–	109860	0.32	0.88	0.85

All runs have a large Reynolds number (Re = ⁠) and Froude number (Fr⁠) below unity.

^aViscosity ratio (⁠⁠) is used instead of l or D.

^bA larger value is used to compensate lower velocities in sidewall boundary layers.

^cThe inlet in Figure 1 is moved to ⁠.

^dInlet instantaneous velocities are constructed from the LES results for E1.

^eInlet instantaneous velocities are constructed from experimental data.

View Large

The measurements were made from 70 points in the plane (Figure 1). Inlet turbulent fluctuations (ITFs) are usually needed for the superimposition to ⁠. We assess two ways to meet this need: using SS to generate ITF from a total of 100 Fourier harmonics; using 2D VM to superimpose a fluctuating vorticity field from a total of 190 vortices. In both, either the turbulence intensity and the eddy-to-molecular viscosity ratio ⁠, and the length scale l, and the hydraulic diameter D, or TKE k and turbulent dissipation rate were supplied. Estimation from the measurements (Ramamurthy et al. 2017) from two planes (⁠ and 0.885) gave ⁠; ⁠; ⁠. A Reynolds-averaged simulation suggested ⁠.

We assess two ways to specify outlet condition: (1) pressure outlet ⁠, with the flow depth equal to 0.2371 m, and (2) pressure outlet, in combination with numerical beach technique (NBT) for damping waves (Park et al. 1999). Other boundary conditions were: on the no-slip bottom and at the no-slip sidewalls. The pressure is zero at the free surface.

The model domain (Figure 1(a)) was discretised into a mesh of hexahedrons. The regions near the bottom, sidewalls were refined using 30 inflation layers. The first layer off a wall had a wall distance ⁠, where y is the distance to the wall, and is the friction velocity. The mesh resolves the viscous sublayer. The between-layer spacing varied smoothly. The mesh sizes satisfy Zang's (1991) guidelines and/or Piomelli's (1993) criteria. The growth rate of sizes was kept within 1.2. The mesh had a low skewness; the average value is below 0.17 and the maximum is below 0.65. With low skewness (Chung 2002), the mesh quality is good for maintaining numerical accuracy and stability.

We assess three different ways to specify the initial conditions:

• (1)

  The left flank had ⁠; the middle region and right flank had water velocity and ⁠, where was uniform at = constant and satisfied (Figure 1(a)). The free surface rose linearly along the transition length, from at to at ⁠, where is the flow depth at the outlet; for ⁠; for ⁠. was uniform in the -direction. The pressure p was hydrostatic at ⁠; and at ⁠. This set of initial conditions is labelled as ‘L’ in Table 1.

• (2)

  For the entire channel, was uniform at a given cross-section (⁠ = constant) and satisfied ⁠. Other remarks are similar to those given in initial condition ‘L’. This set of initial conditions is labelled as ‘N’.

• (3)

  The ensemble averages of from run D1 (Table 1) were used as initial conditions for runs E1 − E8. This set of initial conditions is labelled as ‘S’.

The time step was = 0.01 s for all the LES runs (Table 1), which satisfies the Courant–Friedrichs–Lewy criterion for numerical stability. The model equations were solved in a pressure–velocity coupling manner, using the SIMPLE algorithm. Iterations continued until the convergence criterion of 10⁻⁶ or a maximum of 120 iterations per time step was reached. The LES runs were carried out using Ansys FLUENT software (Ansys 2018) and the outputs were visualised and processed using MATLAB (MathWorks 2022) scripts.

Runs T1 − T3 differed only in turbulence intensity as an input to VM (Table 1). The middle value of was an estimate from experimental data (Ramamurthy et al. 2017). The other input to VM was the hydraulic diameter D. The mean-inflow velocity distribution matched the experimental discharge. Runs T4 − T6 matched the conditions of runs T1 − T3, respectively, but used SS rather than VM. Run T7 matched run T5 but used the viscosity ratio (⁠⁠) rather than D.

Using SS, run T6 produced profiles of similar accuracy (Figure 2(f)–2(j)) as run T2. However, the use of SS surprisingly underpredicted (Figure 2(p)–2(t)). The turbulence intensity was actually large for run T6 than run T2 (Table 1). Among runs T1–T7, run T6 gave the largest underpredictions of ⁠. Run T6 has an average for the five selected (⁠⁠) profiles. Among the same seven runs, run T3 gave the lowest value (0.67) for the coefficient of correlation, among the five selected profiles. The results for run T7 had similar quality as runs T4 − T6. Ai & Mak (2015) reported severe underpredictions by SS compared to VM in their LES study of air flow over a building, without offering validation. We validate the advantage of VM with turbulence intensity (⁠⁠) and hydraulic diameter (⁠⁠), and use them for most of the subsequent runs D1–E9 (Table 1).

In Figure 3(a)–3(c), it is shown that the primary flow was stronger than the secondary flow, i.e., and ⁠, as expected. The data series of both and had a nearly zero mean. The primary flow had a reverse velocity as large as 0.2⁠. The data series of ⁠, and all exhibited significant fluctuations; had larger RMSDs, with a maximum of 0.25⁠.

We examined a near-bed plane and noticed stronger fluctuations near the left sidewall, compared to those at the middle and those near the right sidewall (Figure 1(c), green squares). Near the water surface, near the left sidewall had a mean of ⁠, much smaller than at the middle and near the right sidewall, although the left sidewall location is in a narrower section closer to the transition entrance.

The Kolmogorov −5/3 law is a well-established spectral distribution of turbulent energy in the inertial subrange (Kolmogorov 1941). This study used the law as a benchmark solution for validating LES predictions. We extracted data series of instantaneous velocities over a period of 10 s from the predictions for a series of cells and carried out a Fourier analysis of the data series. Examples of output velocity spectral density, E, from the analysis are plotted in Figure 4. The distributions show a typical energy spectrum for a turbulent flow, the existence of energy containing eddies and a well-fitted −5/3 slope (dashed lines) in the inertial subrange. We conclude that the LES runs properly resolved a major part of turbulence energy and captured the cascading of energy to smaller eddies.

Averaging instantaneous LES output is needed for practical use of LES results. What is the required minimum length of data series or minimum sample size? We answer this question by evaluating the accumulative average and the autocorrelation function over a time period T elapsed following a spin-up period ⁠. Examples of are plotted in Figure 3(d)–3(f). Take as an example = 80 in Figure 3(d). was the average of 1,000 instantaneous values from = 70 to 80. At ≥ 80, changed very little. This is true for other locations in the model channel. Thus, 10 or 10 s is suitable for ensemble averaging, corresponding to sample size = 1,000.

We confirmed the suitability by further calculating ⁠. Examples of values are plotted in Figure 3(g)–3(i). If two quantities separated by k are essentially uncorrelated, the selected number of ensembles N is sufficient and suitable for averaging (Box et al. 2015). measures the correlation between the univariate data series at time t and at time ⁠, where k increases from 0, … , 10 s. The wavy curves (Figure 3(g)–3(i)) fluctuated in wavelength and amplitude, depending on the direction and location. However, they had some common features at due to a perfect self-correlation; r_k decayed with increasing from zero and became bounded by for s. The curves oscillated within a range of small values for ⁠. The oscillations might be associated with the Hurst phenomenon (Nordin et al. 1972; Meneveau & Sreenivasan 1991; Dimitriadis et al. 2021). These low values of towards mean that the instantaneous velocities separated by a 10 s time span were uncorrelated. The preceding discussion served the purpose of confirming the suitability of using 10 s for averaging.

Hereinafter, data sampling from the LES output started after the spin-up period (or ⁠) at a sampling frequency and continued until ⁠. In summary, the sampling time duration is ⁠. This choice is justified on the following grounds: by an order of magnitude; is 8 times the advection time scale ⁠; the cumulative average of velocity supports the choice s; the sufficiency of s is manifested in the autocorrelation function. The ensemble-averaged quantities presented hereinafter use ⁠.

Note that in a large water body, turbulent motions can contain slow oscillations at large scales. To capture such oscillations, the sampling time duration needs to be longer. For example, Soulsby (1980) and Walter et al. (2011) used as long as 8–12 min for oceanic motions in order to capture the desired turbulent length scales and obtain quasi-stationary statistics.

The distributions displayed complex 3D flow structures. The near-bottom flow had small pockets of strong velocities in the transition and nearby regions. There was a hysteresis of the occurrence of separation; the further away from the water surface, the further to the transition entrance. The near-surface flow separated in the transition along both sidewalls and triggered eddies in both flanks. The reverse flow had a maximum magnitude of −0.46⁠.

Extending the upstream section length (from for run D1 to 8 for run D3, Table 1) improves the accuracy of first- and second-order statistics of turbulence, as will be demonstrated later. Here, we compared the distributions of ensemble-averaged TKE in the aforementioned three selected horizontal planes. An example of the distributions is shown in Figure 5(b). In the transition, compared to near-surface and near-bottom, the mid-depth had large values, with a maximum of 0.11 (Figure 5(b)). Large pockets of strong occupied the core of recirculation zones in the left flank. Near the surface, many small pockets of strong flanked a strong velocity core, within which the magnitude of was relatively weak. Near the bottom, there were pockets of weak magnitude in the two flanks due to the proximity of channel boundaries.

We quantitatively located turbulent eddies using the Okubo–Weiss parameter ⁠, where is the normal strain component, is the shear strain component and is the relative vorticity. Horizontal areas of ⁠, where is a certain threshold, are dominated by eddies, with a core in the local minima of ⁠. The areas exhibit large vorticity and strong circulation density of the velocity field. Areas of are dominated by strain. Areas of have background vorticity and strain. We chose ⁠, where is the standard deviation of values in the flow area. For more details, refer to Zeng & Li (2022).

As an example, contours of normalised for D3 are shown in Figure 5(c). The normalisation used an eddy length scale m. Eddies larger than can be discerned from the computed flow field. For D3, based on and ⁠; for the aforementioned three horizontal planes.

In the transition, a string of strong alternating eddy and strain pockets originated near the transition entrance and developed along the left- and right side of the core of strong velocity (Figure 5(a)). These eddies and strains were coherent turbulence structures and local-scale flow features. The string on the left side traced a meandering path, whose cross-stream amplitude increased as eddy patterns developed towards downstream. At the transition exit, the cross-stream amplitude of the eddy patterns had length scales comparable to ⁠. The patterns' development was accompanied by the weakening of local rotations.

After the transition, local rotations persisted longer near the surface than near the bottom and diminished at about ⁠. The string on the right side virtually hugged the sidewall, with little cross-stream shifting of eddy positions at the mid-depth (Figure 5(c)). This was not the case near the bottom and near the surface, where the cross-stream amplitude increased slightly.

The hit-rate and a factor of two (Okaze et al. 2021) are useful for assessing the quality of alternative LES set-ups, parameters and assigned parametric values. The better the predicted values of a flow variable compared with measured values, the higher the and scores. The quality is acceptable (i.e., the prediction is in reasonable agreement with the measurement, which is the goal) if and ⁠. A short description of the metrics is given below. Let ⁠, ⁠, …, denote n measured values of from experiments, and ⁠, ⁠, …, denote the corresponding computed values. The validation metrics are expressed as: ⁠, where if or and otherwise, is an allowed relative deviation, and is an absolute deviation; ⁠, where = 1 if and otherwise. According to Schatzmann et al. (2010), = 0.25. The threshold is 0.20 for and 0.06 for ⁠, as estimated following previous studies. Table 1 lists the scores for the runs, which are based on comparisons with 1,184 measured values of streamwise velocities, lumped together from the 11 dotted cross-sections (CSs) in Figure 1(c). The scores are interpreted as follows:

About initial condition, the uniform distribution (for run D5) is unsatisfactory because both scores are out of acceptable range. A plausible reason is that the turbulence development is history-dependent, as demonstrated in Zeng & Li (2022). Once established, the flow tends to hug the left or right wall, depending on the flow history. In the broad sense, the initial condition used in a LES should reflect the flow history of the experiment in question. Either the initial condition ‘L’ or ‘S’ does in this paper.

For inlet mean-flow, a 1/14th power-law (for run D4) is unsatisfactory because the score is too low. A 1/7th law (for D3) achieves acceptable scores. In fact, this law is a good approximation of turbulent boundary layer profiles at high Re, as is the case in this study. Also, the 1/7th law is easier/more efficient to implement than the well-established logarithmic law.

Extending the upstream section length to eight times the depth (for run E9) from 2.67 times (run E1) improves scores. In particular, the extension improves the prediction accuracy of second-order turbulence statistics. It is worth mentioning that the extension increases computational nodes by 7% and computing time by 37% in this study. For similar LES problems of free-surface turbulence, we recommend an adaption length of about eight times the depth, when the inlet conditions reflect the real flow reasonably well. Some previous researchers suggest an adaptation length of three to four times the building height in air flow problems.

The VM (for run D1) slightly improves scores over the SS (for run D6), when the generated ITF is superimposed on a mean-flow profile (given by e.g., a 1/7th law) that is laterally uniform in the absence of secondary flow and eddy motion. A VM generates stronger ITF than a SS, for the same values of input parameters (⁠ and ⁠). The VM is more reliable in building secondary eddies and associated rotations, representing secondary flows of Prandtl's second kind and causing velocity shear. A similar discussion is given in previous LES studies of air flow turbulence (Thordal et al. 2019).

For the first time, this study explores mapping 3D inlet mean-flow from experimental data. In this new case, the SS seems slightly better than the VM, as indicated by comparing run E2 to run E3, run E6 to run E5 and run E8 to run E7. Plausible explanations are that the mapped velocity already contains secondary flow eddies, and it becomes less critical or even redundant to add turbulent eddies from the VM. The same redundancy does not occur in the SS because it yields little influence on secondary velocity.

All the three SGS stresses models: DSM, WALE and DKE (for runs E1, E5 and E7) enjoy a similar degree of success. This is regardless of whether combined with l or with D is used as parameters, as can be seen from comparing run E3 to run E5. The most likely reason is that the mesh for the runs resolves all the important large eddies. The SGS model for smaller eddies has no critical influence on the LES output.

It is constructive to discuss the and scores for CSs in different zones (Figure 1(c)). For the three CSs in zone 1 before the transition (Figure 1(c)), the scores (>0.95) are much better for runs E1 − E8 than runs D1 − D6, due to the fact that runs E1 − E8 map inlet mean-flow from 3D experimental data, preserving the observed flow behaviour before the transition. The scores are not as high just after the transition to some distance downstream (⁠2.49). The influence of the inlet condition on the local flow behaviour weakens and that of separation strengthens, and it is more challenging to accurately predict the behaviour of stronger separated flow.

After the transition, the scores for zones 2, 4 and 6 are good when the upstream section length is eight times the depth (for run D3). The scores change little by using NBT, which is clear from comparing run E4 to run E3 and run D2 to run D1. It is difficult to achieve high scores for zones 4 and 5, where separation is vigorous. The scores barely pass the acceptable threshold for almost all the runs in Table 1. The score is below the acceptable threshold for zone 6 when using 3D experimental data to give inlet mean-flow and the VM to generate fluctuations (for runs E1, E3 − E5 and E7).

The recommended length of an upstream channel section of eight times the flow depth has practical implications. In the studies of Nashta & Garde (1988) and Shettar & Keshava Murthy (1996), the depths are, respectively, 0.059 − 0.12 m and 0.05–0.06 m. The corresponding recommended lengths are 0.47 − 0.96 m and 0.4 − 0.48 m. In Chow (1959), the flume-to-canal expansion has a depth of 1.15 m, the recommended length being 9.2 m. For a concrete siphon expansion with a design depth of 1.25 m (Hinds 1928), the recommended length is 10 m. This study provides some guidelines for the practical choice of channel lengths.

For those runs where an SS was used to generate inlet velocity fluctuations, the LES methods give underestimations of RMSD in the upstream channel expansion before the expansion, when compared to measured values on a point-to-point basis. Improvement can be achieved by mapping inlet mean-flow from experimental data of 3D velocity, combined with a VM (for runs E1, E5 and E7). The issue is that such data is often not available in many cases. An alternative is to use a 1/7th power-law together with the extended approach section length (for run D3) or a precursor simulation. One limitation is excessively high computing costs. Future work should investigate the influence of varying downstream channel length on LES results.

This study has made unique contributions in two aspects. First, the two-phase flow with a free surface advances from the commonly used rigid-lid approximation. Free-surface variations are important to non-prismatic channel design. Second, the combination of cumulative averaging, autocorrelation and spectrum analysis is an improved method for determining the sample size for ensemble averaging. It avoids trial-and-error and is more efficient and accurate, compared to the traditional method of trying different averaging times.

In this paper, two-phase LES of turbulent flow is simulated in a non-prismatic open channel. The significance of this work lies in progressively exploring strategies for LES model set-up, selection of proper parameters and parametric value assignments. We validate the LES predictions using experimental data and benchmark solution. The validated strategies can be used to study turbulence under other similar conditions with good accuracy. We show the practical results exerted from this study, which are related to the velocity distributions, spectral distribution of flow energy and hit-rate for validations in other similar studies.

This study reveals 3D complex turbulence structures in the non-prismatic channel, as a result of processes mingled together, involving eddy motions translated from upstream and locally generated by asymmetrical separation. It also reveals the influence of flow history on the turbulence development. One implication is that there may be a need to test a subtle approach to specify the initial conditions in a two-phase LES. This newly recognised challenge has not yet received much attention in the literature.

We show that mapping mean-velocity distribution from measurements, combined with the spectral synthesiser approach to velocity fluctuations, gives satisfactory inlet conditions. An alternative is a 1/7th power-law, combined with the VM. To improve the prediction and control of turbulent behaviour, we recommend adopting a channel length of eight times the depth of the channel.

For the first time, we report instantaneous TKE and coherent turbulence structure in a non-prismatic channel. The practical significance is that the former is of relevance to turbulence control and hydraulic efficiency improvements, whereas the latter helps address channel erosion and instability risks.

This study received financial support from the Natural Sciences and Engineering Research Council of Canada through Discovery Grants held by S.S.L.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.