3.1 LES formulation

Mass-weighted low-pass spatial filtering is applied to obtain the Cartesian tensor version of the compressible LES equations. The present flow is subsonic and adiabatic with jet total and ambient static temperatures equal. Thus, the same treatment of the energy equation has been adopted as in the LES simulations of Karabasov et al.^{(Reference Karabasov, Afsar, Hynes, Dowling, McMullan, Pokora, Page and McGuirk4)} (Ma_J = 0.9) and Wang and McGuirk^{(Reference Wang and McGuirk5)} (supersonic jet), namely that the steady state solution of constant total enthalpy (or temperature) is a good approximation even in an unsteady flow. This produced good agreement with measurements of both mean and turbulence properties in these earlier works and is assumed here for the chosen NPR = 1.5 test case (Ma_J = 0.71). The static temperature ( ${\tilde T_s}$ ) is calculated from the specified inlet total temperature ( ${\tilde T_t}$ ) and predicted velocity field; the density is evaluated from the equation of state for an ideal gas. The governing equations solved were:

(1) \begin{equation}\frac{{\partial \bar \rho }}{{\partial t}} + \frac{\partial }{{\partial {x_j}}}(\bar \rho {\tilde u_i}) = 0\end{equation}

(2) \begin{equation}\frac{\partial }{{\partial t}}(\bar \rho {\tilde u_i}) + \frac{\partial }{{\partial {x_j}}}(\bar \rho {\tilde u_j}{\tilde u_i}) = - \frac{{\partial \bar p}}{{\partial {x_i}}} + \frac{\partial }{{\partial {x_j}}}\left[ {\mu \left( {\frac{{\partial {{\tilde u}_i}}}{{\partial {x_j}}} + \frac{{\partial {{\tilde u}_j}}}{{\partial {x_i}}}} \right) - \frac{1}{3}\frac{{\partial {{\tilde u}_k}}}{{\partial {x_k}}}{\delta _{ij}}} \right] - \frac{\partial }{{\partial {x_j}}}(\overline {\rho {u_i}{u_j}} - \bar \rho {\tilde u_i}{\tilde u_j})\end{equation}

(3) \begin{equation}\bar p = \bar \rho R{\tilde T_s}\end{equation}

where $(\overline {\rho {u_i}{u_j}} - \bar \rho {\tilde u_i}{\tilde u_j})$ is the residual (SGS) stress tensor $\tau_{ij}^{sgs}$ . To close the momentum equations two SGS models were used; the first is the standard Smagorinsky^{(Reference Smagorinsky32)} model:

(4) \begin{equation}\tau _{_{ij}}^{sgs} = ({\overline {\rho {u_i}u} _j} - \bar \rho {\tilde u_i}{\tilde u_j}) = - 2{\mu _{sgs}}{\tilde S_{ij}} + \frac{2}{3}\bar \rho {{{\tilde {\rm k}}}_{sgs}}{\delta_{ij}}\end{equation}

(5) \begin{equation}{{{\tilde {\rm S}}}_{ij}} = \frac{1}{2}\left(\frac{{\partial {{\tilde u}_i}}}{{\partial {x_j}}} + \frac{{\partial {{\tilde u}_j}}}{{\partial {x_i}}}\right)\quad {\mu _{sgs}} = \bar \rho l_{sgs}^2S\ S = {(2{\tilde S_{ij}}{\tilde S_{ij}})^{\frac{1}{2}}}\ \end{equation}

(6) \begin{equation} l_{\textit{sgs}}=C_{s}\Delta\quad \Delta=(\delta x\delta y\delta z)^{1/3} \qquad \text{and near walls:}\quad l_{\textit{sgs}}=C_{s}\Delta(1-e^{n^{+}/A^{+}}) \end{equation}

where l _sgs is the SGS length scale and Δ is a mesh-based length scale (cube root of local cell volume). In the near wall region van Driest^{(Reference Van Driest33)} damping is used: n⁺ is wall normal distance and A⁺ = 26. Literature reveals a wide range (0.09–0.18) can be found for the Smagorinsky constant C _S with a tendency for the upper end; in the present work C_S was set to 0.15 as in earlier LES studies Refs. (Reference Karabasov, Afsar, Hynes, Dowling, McMullan, Pokora, Page and McGuirk4) and (Reference Wang and McGuirk5).

The second SGS model follows Piomelli and Geurts^{(Reference Piomelli, Geurts, Leschziner, Bontoux, Geurts, Launder and Tropea34)}; this retains the SGS viscosity definition but calculates l _sgs differently. The motivation is to avoid the link between l _sgs and mesh size, which has several undesirable features (Piomelli and Geurts^{(Reference Piomelli, Rouhi and Geurts35)}). In the Piomelli and Geurts^{(Reference Piomelli, Geurts, Leschziner, Bontoux, Geurts, Launder and Tropea34)} model l _sgs is set proportional to a fraction of the (estimated) integral length scale L; dimensional analysis relates this to the resolved scale turbulence energy and vorticity magnitude:

(7) \begin{equation} L\propto\frac{k^{1/2}}{\omega}\quad \textit{where}: \quad k=\frac{1}{2}\left\langle(\tilde{\text{u}}_{i}-<\tilde{\text{u}}_{i}>)(\tilde{\text{u}}_{i}-<\tilde{\text{u}}_{i}>)\right\rangle\quad \omega=\langle{\tilde\omega}_{i}{\tilde\omega}_{i}\rangle^{1/2}\end{equation}

Both k and ω are statistically averaged since L is a statistically averaged quantity. Resolved rather than fluctuating vorticity is used Ref. (Reference Piomelli, Geurts, Leschziner, Bontoux, Geurts, Launder and Tropea34) to ensure that in a laminar boundary layer the turbulence length scale is zero. Piomelli et al.^{(Reference Piomelli, Rouhi and Geurts35)} have noted that the calculation of L can be improved by including SGS quantities and proposed a modification along these lines. This version has, however, so far only been tested in homogeneous isotropic turbulence and fully developed channel flow, and the earlier Ref. (Reference Piomelli, Geurts, Leschziner, Bontoux, Geurts, Launder and Tropea34) model was adopted in the current work. SGS length scale and viscosity are calculated from:

(8) \begin{equation}\;{l_{sgs}} \propto L \propto \frac{{{k^{\frac{1}{2}}}}}{\omega } = {C_{PG}}\frac{{{k^{\frac{1}{2}}}}}{\omega }\quad \quad {\mu _{sgs}} = \bar \rho C_{PG}^2\frac{k}{{{\omega ^2}}}S\end{equation}

C_PG is a model coefficient, determined by Piomelli and Geurts^{(Reference Piomelli, Geurts, Leschziner, Bontoux, Geurts, Launder and Tropea34)} in a parametric study of a fully developed channel flow to minimise error in the skin friction coefficient. The minimum was found for C_PG = 0.12, although, in the flow considered, the error proved insensitive to variations over the range of 0.05–0.15. The model has not previously been applied to a convergent duct flow and, given the wide range of optimum value for C_PG, further consideration in choosing a value for C_PG was deemed appropriate. Piomelli et al^{(Reference Piomelli, Rouhi and Geurts35)} suggested that: “the optimum value of the coefficient may be chosen to minimise the total simulation error in prediction of a set of selected quantities known from experiments.” Since the integral length scale L plays a dominant role in the model, the value of C_PG is here determined by minimising the error in prediction of L in a zero pressure gradient boundary layer (see Section 4.2).

3.2 CFD code details

The CFD code used adopts a hexahedral, cell-centred, multi-block, structured mesh and a pressure-based approach. Originally written for RANS calculations of compressible flows it was converted to LES for prediction of multiple impinging jets in Page et al.^{(Reference Page, Li and McGuirk36)}. The same spatial and temporal discretisation schemes were adopted as in the LES studies of Karabasov et al.^{(Reference Karabasov, Afsar, Hynes, Dowling, McMullan, Pokora, Page and McGuirk4)} and Wang and McGuirk^{(Reference Wang and McGuirk5)}. Second order spatial discretisation was employed, and temporal advancement made using a low-storage Runge-Kutta 4th-order (5-stage) scheme; in all simulations the maximum CFL was always less than 0.25. References (Reference Karabasov, Afsar, Hynes, Dowling, McMullan, Pokora, Page and McGuirk4) and (Reference Wang and McGuirk5) provide evidence of acceptable numerical robustness and dispersion/dissipation characteristics.

3.3 R²M approach for LES inlet condition generation

The Xiao et al.^{(Reference Xiao, Dianat and McGuirk22)} R²M technique was applied in the present work. This attaches an inlet condition (IC) domain to a main simulation (MS) domain (i.e. the nozzle of Fig. 1, see Fig. 4). Rescaling/recycling are used to match specified target statistics throughout the IC domain, within which statistical properties are assumed streamwise homogeneous, effectively implying that the MS domain inflow conditions are changing only slowly axially. The IC domain length is chosen purely to ensure streamwise turbulent structures (axial 2-point correlations) are not constrained by IC domain end boundary conditions (after statistical convergence of the simulation this has been retrospectively checked (see Section 4.2). Figure 4 shows reference plane labelling: IC domain inlet/outlet: A (x/D = 0), C (x/D = 2.5) (also MS domain inlet); re-cycling plane B (x/D = 2.08), and D, E, and F at convergence start/end (x/D = 3.0, x/D = 3.64), and nozzle exit (x/D = 4.17).

Fig. 4. Inlet Condition and Main Simulation (nozzle) domains.

The method generates spatially and temporally correlated time-series, which are consistent with specified (target) statistical properties (mean velocity $U_i^T$ and three turbulence normal stresses (rms values $u_{rms,i}^T$ ) in Cartesian components) chosen to match measurements. The nozzle inlet duct is axisymmetric, and the mean velocity in the measurements of Ref. (Reference Trumper, Behrouzi and McGuirk29) is thus available on a single radial line in cylindrical polar co-ordinates (U_x = f(r), U_r = U_θ = 0); radial profiles for the associated turbulence rms values (which were not measured in Ref. (Reference Trumper, Behrouzi and McGuirk29), were created by rescaling boundary layer DNS data from Spalart^{(Reference Spalart37)} as described next. (NB, this therefore constitutes an approximation in the inlet data.)

A low-Re k-ε RANS prediction of flow in a long pipe of the same diameter as the inlet duct was first carried out starting from a flat velocity profile corresponding to the experimentally measured mass flow for NPR = 1.5. A downstream location in this pipe was then selected where the predicted mean axial velocity best matched the measured target data profile. The Spalart data were then rescaled to create axisymmetric target turbulence rms data using Equation (9) below (only radial component shown for illustration).

(9) \begin{equation}\left[\sqrt{\overline{u^{2}_{r}}}\right]_{target} (r)=\left[\sqrt{\overline{u^{2}_{r}}}\right]_{SPALART}(r)\left[\frac{k^{1/2}_{RANS}(r)}{k^{1/2}_{SPALART}(r)}\right] \end{equation}

Finally, since the LES code solves for Cartesian components, the statistical values for the 1st and 2nd moments in cylindrical polar form were transformed to Cartesian components using standard relations.

The R²M technique was then implemented as follows:

1. (i) Every nth time step new estimates of statistical moments of velocity were evaluated:

   (10) \begin{equation}{\left\langle {{{\tilde u}_i}} \right\rangle ^{new}}\left( {y,z} \right) = \alpha {\left\langle {\tilde u_i^n(x,y,z,{t_n})} \right\rangle _x} + \left( {1 - \alpha } \right){\left\langle {{{\tilde u}_i}} \right\rangle ^{old}}(y,z)\end{equation}
   (11) \begin{equation}u_{rms,i}^{new}\left( {y,z} \right) = {\left( {\alpha {{\left\langle {{{[\tilde u_i^n\left( {x,y,z,{t_n}} \right) - {{\left\langle {{{\tilde u}_i}} \right\rangle }^n}(y,z)]}^2}} \right\rangle }_x} + (1 - \alpha ){{\left[ {u_{rms,i}^{old}(y,z)} \right]}^2}} \right)^{1/2}}\end{equation}

2. (ii) The instantaneous velocity at all mesh points within the IC domain was then rescaled using (10) and (11) and the target statistical data:

   (12) \begin{equation}\tilde u_i^n\left( {x,y,z,{t_n}} \right) = \frac{{u_{rms,i}^T(y,z)}}{{u_{rms,i}^{new}\left( {y,z} \right)}}\left[ {\tilde u_i^n\left( {x,y,z,{t_n}} \right) - {{\left\langle {{{\tilde u}_i}} \right\rangle }^{new}}(y,z)} \right] + U_i^T(y,z)\end{equation}

3. (iii) Finally, the velocity field obtained from (12) was recycled from plane B to A and the next time step begun.

When the IC simulation has reached a statistically stationary state, the transfer of IC domain data to the MS domain inlet (plane C) ensures a time-varying inlet condition for the main simulation with the desired statistics and self-consistent spatial/temporal correlations.

3.4 Mesh design

Figure 5 presents the multi-block mesh structure adopted (alternate lines only); the block topology of central H-grid and outer O-grids (four blocks) is optimal for axisymmetric geometry and avoidance of centreline stiffness. LES mesh design followed guidelines using a-priori RANS solutions (low Re Launder-Sharma model used here) as in the approach proposed by Gant^{(Reference Gant38)} and Dianat et al.^{(Reference Dianat, McGuirk, Fokeer and Spencer39)}. The first guideline was used to select a mesh density appropriate for resolution of turbulence structures in regions not too close to a wall. Pope^{(Reference Pope40)} has proposed that for a “well-resolved simulation” in such regions the mesh should capture greater than 80% of fluctuating energy, a criterion also used by Celik et al.^{(Reference Celik, Cehreli and Yavuz41)} and shown in Dianat et al.^{(Reference Dianat, McGuirk, Fokeer and Spencer39)} to imply that L/Δ should be > 12 (L = k^3/2/ε and Δ = mesh size).

Fig. 5. (a) Nozzle geometry/mesh. (b) IC domain mesh.

The second guideline was used to design the near-wall mesh design. Choi and Moin^{(Reference Choi and Moin42)} have reported that typical near-wall cell sizes for wall-resolved simulations have varied widely in published LES studies: Δx⁺ = 50−130, Δy⁺ <1 with 10–20 nodes in the boundary layer, and Δz⁺ = 15–30 (for streamwise, wall normal, wall parallel directions).

Preliminary RANS calculations on various meshes were carried out to produce a final mesh that satisfied the L/Δ criterion at mesh nodes greater than 10%D from the nozzle wall, had a near-wall mesh with Δx⁺ ~200, Δz⁺ ~30 and a first node next to the wall at Δy+~0.3. With 65 nodes across the boundary layer at nozzle inlet, 24 at end of convergence, and 55 at nozzle exit the resolution of the boundary layer is significantly finer than that used by Bres et al.^{(Reference Bres, Jordan, Jaunet, Le Rallic, Cavalier, Towne, Lele, Colonius and Schmidt28)} (only 10–20 cells inside the boundary layer). The streamwise mesh is larger than recommended, but examination of results at nozzle inlet (see Section 4.2) showed good agreement with mean velocity, turbulence normal stress and integral length scale measurements from equilibrium boundary layer data. Further, Warnack and Fernholz^{(Reference Warnack and Fernholz20)} observed streamwise length scales grow significantly in accelerated boundary layers (by a factor of 4 in the experiments of Ref. (Reference Warnack and Fernholz20) which had a slightly lower K_max than in the Trumper et al.^{(Reference Trumper, Behrouzi and McGuirk29)} data), and a moderate increase in Δx⁺ was thus considered acceptable. The final nozzle mesh had 100 × 90 × 90 nodes (x, y, z directions) in the H-grid and 100 × 60 × 90 nodes in each of the four outer blocks, a total of 3 million cells. The IC domain length was chosen to be five inlet boundary layer thicknesses, and his mesh increased total grid size to 4.6 million cells. Overall this is a similar mesh density as used by Bres et al.^{(Reference Bres, Jaunet, Le Rallic, Jordan, Colonius and Lele6)}, whose nozzle geometry was four times as long as that considered here.

At nozzle inlet a spatially and temporally varying total pressure was applied; this was evaluated at each time step by extrapolating the static pressure upstream from cells just inside the inlet plane and converting this to P_T using an instantaneous Mach No. calculated from a velocity extracted from the R²M procedure. Simulations were run for 10T_FT (calculated using U₀) before sampling for statistical properties (every 50 time steps) was carried out for a further 10T_FT (a 20% increase in averaging time showed no influence on statistics).

4.1 Optimisation of Piomelli and Geuerts SGS model constant

The Piomelli and Geurts SGS coefficient C_PG was optimised for boundary layer flow using simulations carried out in the IC domain alone. Target data comprised the measured inlet mean velocity and turbulence normal stress profiles as detailed in Section 3.3. Simulations were conducted with several C_PG values. When a statistically stationary state had been reached the spatial 2-point axial correlation with varying axial separation R ₁₁ was evaluated on a plane within the IC domain (at x/D = 0.415). The axial integral scale L ₁₁ was then calculated from R ₁₁:

(13) \begin{equation}{R_{11}}(x,y,z,dx) = \frac{{\left\langle {u(x,y,z,t)u(x + dx,y,z,t)} \right\rangle }}{{\sqrt {\left\langle {{u^2}(x,y,z,t} \right\rangle } \sqrt {\left\langle {{u^2}(x + dx,y,z,t} \right\rangle } }}\end{equation}

(14) \begin{equation}{L_{11}}(x,y,z) = \int\limits_0^\infty {{R_{11}}(x,y,z,dx')dx'} \end{equation}

Figure 6 shows the correlation function for two points in the boundary layer one near the wall (y⁺ ~ 20) and one near the boundary layer outer edge (y’/δ = 0.95).

Fig. 6. Two-point spatial correlation near wall (black), near outer edge (red) (x measured in m from IC domain inlet).

The correlation functions have their first zero crossing well within the IC domain length (this was used as the integration upper limit in evaluating length scales). Thus, the IC domain size seems acceptable for calculating the integral scale, although to capture the correlation shape at large separations as the correlation asymptotes to zero a longer domain would be preferable. Figure 7 shows the predicted variation of L ₁₁ across the boundary layer for two values of C_PG (five values in all were explored between 0.005 and 0.12 spanning the low error range observed in Piomelli and Guerts^{(Reference Van Driest33)}, although only two cases are shown here for clarity). Figure 7 indicates that C_PG = 0.01 gives best agreement for the integral length scale when judged against the boundary layer data of Antonia and Luxton^{(Reference Antonia and Luxton43)}. In terms of near wall gradient and the value at the boundary layer outer edge, C_PG = 0.01 clearly shows best performance. Note the value of L ₁₁ at the layer edge is quite sensitive to C_PG, unlike the skin friction parameter in the channel flow case. Based on these observations, all simulations below were conducted with C_PG = 0.01.

Fig. 7. Axial integral length scale symbols – Antonia and Luxton^{(Reference Antonia and Luxton43)} data lines – LES (C_PG: blue = 0.12, red = 0.01).

Fig. 8. LES R²M (P+G SGS) predictions, mean axial velocity, (left) plane z/D = 0, (right) plane C.

4.2 Assessment of R²M method performance against measured data

Figure 8 shows predicted mean axial velocity contours within the IC domain and part of the MS domain on plane z/D = 0 (left) and plane C (right). Contours on the first plane illustrate that streamwise homogeneous flow is produced within the whole IC domain, with constant boundary layer thickness, showing no sign of distortion as the solution passes through the IC/MS domain interface (x/D = 2.5). The cross-stream contours on the radial/azimuthal plane C (MS inlet conditions) display an axisymmetric flow with the correct boundary layer thickness (see Table 2 below). Results for the Piomelli and Guerts SGS model simulation are shown in this Section, but there is little difference between the SGS models, both perform well for a zero pressure gradient boundary layer.

Fig. 9. LES R²M (P+G SGS) inlet mean axial velocity profile.

Quantitative assessment of the ability of the R²M technique to generate velocity time series with statistics that match the target data is provided in Figs. 9 and 10. In Fig. 9 the mean axial velocity profile produced for the nozzle inlet boundary condition is plotted in log-law format together with the experimental data. The LES prediction agrees well with the measured data, in both log-law region portion (~250 < y⁺ < ~5000) and low Reynolds number buffer region (y⁺ < 100).

In Fig. 10 the predicted normal stress profiles are plotted along four radial lines (0⁰, 90⁰, 180⁰, and 270⁰). The excellent axisymmetry of these 2nd order statistics demonstrates a sufficiently long averaging time. Anisotropy between streamwise, wall normal, and wall parallel components has been captured well (significantly better than in Bres et al.^{(Reference Bres, Jaunet, Le Rallic, Jordan, Colonius and Lele6)}). Peak stress values in the near wall region (y/D < 0.05) are slightly overpredicted for axial (14%) and transverse (10%) components, but underpredicted for the wall normal stress (25%). The likely cause for these errors is the use of Equation (9), which effectively assumes the entire profile obeys outer layer similarity. Fernholz and Finlay^{(Reference Fernholz and Finley30)} have analysed experimental data to show that whilst this is true for Re _θ > 5000 and for y^’/δ > 0.1, Reynolds number effects are still important for y^’/δ < 0.1 and inner scaling using wall units should have been used for this region. The near wall agreement could have been improved by extracting rms data using separate inner/outer similarity relations as given in Fernholz and Finlay^{(Reference Fernholz and Finley30)}. An alternative would be to use experimental boundary layer data taken at higher Reynolds number as provided by Carlier and Stanislas^{(Reference Carlier and Stanislas44)}.

Table 2 Global boundary layer parameters at nozzle inlet for NPR = 1.5

Fig. 10. LES R²M (P+G SGS) inlet Reynolds normal stress profiles: blue: axial (n = 1), red: radial (wall normal, n = 2), green: azimuthal (wall parallel, n = 3) symbols are target data.

Figure 11 presents a power spectral density plot for axial turbulence fluctuations at a point in the log-law region (y⁺ ~1000). This shows that the R²M technique has enabled realistic energy spectra to be created, with the expected −5/3 behaviour displayed for around two decades of PSD (10⁻²–10⁻⁴) before enhanced fall-off at higher frequencies where SGS viscosity is effective. There is no evidence of the recycling frequency errors reported by Morgan et al.^{(Reference Morgan, Larsson, Kawai and Lele10)}; these must be present as recycling establishes a perfect correlation between inlet and recycling planes, but this artificial frequency component is clearly in the present case at significantly lower energy than the resolved scale turbulence.

Finally, Table 2 provides predicted overall boundary layer parameters at nozzle inlet. Comparison between measurements and the LES solution indicates good agreement for all parameters.

Fig. 11. LES R²M (P+G SGS) PSD at nozzle inlet.

Fig. 12. Instantaneous u velocity contours on three cross-sections between planes C and D LES-R²M (Smagorinski SGS).

In summary, examination of statistical properties of the nozzle inlet time-series generated by the R²M approach – mean velocity (Fig. 9 and Table 2), 2nd moment turbulence quantities (Fig. 10), and integral length scale (Fig. 7) – shows that it has performed well, producing boundary conditions for the nozzle flow simulation that are in close accord with measured data with time-series data over the inlet plane which possess physically plausible spatial correlations (Fig. 6) and energy spectrum (Fig. 11).

4.3 Flow development within the nozzle

Figure 12 provides an illustration of the importance of an LES inlet condition method that creates physically realistic and self-consistent turbulence structures. Instantaneous axial velocity contours are displayed on three planes between nozzle inlet and start of convergence. Large-scale energy-containing flow structures are clearly visible within the boundary layer and sustained throughout the inlet section up to the start of nozzle convergence. Without the benefit of spatial and temporal internal consistency, the simulation would lead to modification of the inlet structures – in all probability their decay – with the flow having insufficient time for fluctuations to be re-established before the acceleration region is entered. LES results are shown this time from the Smagorinski SGS model, but again the influence of the SGS model at this stage of the flow is relatively small.

Figure 13 indicates how the turbulent structures develop under the influence of the nozzle acceleration. Instantaneous U and W velocity contours are shown on plane z/D = 0.0 (this time from the LES-R²M (P+G SGS) simulation). Axial velocity contours (Fig. 13 (top)) are shown for a region from the centre of the IC domain through to nozzle exit. The imposition of a favourable pressure gradient is seen to lead to a significant reduction in the size and amplitude of the turbulent eddies, and rapid reduction in boundary layer thickness. W contours (Fig. 13 bottom) (plotted over the entire IC and MS domains) characterise resolved-scale large-vortical structures revealed by regions of alternating W-velocity sign. These structures are angled away from the wall as observed in boundary layer measurements and are stretched longitudinally as they pass through the convergence, reflecting the Warnack and Fernholz^{(Reference Warnack and Fernholz19)} observation of increased streamwise integral scales. Unlike the U-structures, the W-structures persist to nozzle exit.

Fig. 13. LES R²M (P+G SGS) instantaneous U (top) and W (bottom) on plane z/D = 0.

Fig. 14. LES R²M (P+G SGS) predicted time-mean axial velocity on z/D = 0 plane.

The predicted time-averaged axial velocity field within the nozzle on a vertical diametral plane (z/D = 0) is displayed in Fig. 14, again for the P+G SGS model. Flow acceleration already starts to appear ~0.4D upstream of the convergence, with the thick inlet boundary layer drastically reducing in thickness by half-way through the convergence section. Evidence of a much thinner wall boundary layer re-appears in the nozzle exit parallel extension growing from the convergence/parallel extension corner. The contour plot also shows a small local region of acceleration/deceleration associated with flow curvature caused by the same corner; the effect of this is very significant and will be discussed further below.

It is in the convergent nozzle section where the favourable pressure gradient will exert a strong influence on the turbulence structure; the normal stress anisotropy may depart significantly from its inlet equilibrium structure (Fig. 10). Evidence of this is revealed by in Fig. 15, which displays turbulence statistics at start (plane D) and end (plane E) of the convergence section; results with both SGS models are presented. Note that profiles are non-dimensionalised using local boundary layer edge freestream velocity and boundary layer thickness to allow easier comparison, since acceleration/recovery processes change boundary layer size significantly – δ = 21mm (at C), = 0.3mm (close to plane E), and = 12mm (at F).

Fig. 15. LES R²M with P+G (left) and Smagorinski (right) SGS models Reynolds normal stress profiles, $\text{z/D}=\text{0} $ plane, (a) plane D, (b) plane E, blue $-\overline{\textit{u}^{\text{2}}}$ axial ( $\textit{n}=\text{1}$ ), red $-\overline{\textit{v}^{\text{2}}}$ radial (wall normal, $\textit{n}=\text{2}$ ), green $-\overline{\textit{w}^{\text{2}}}$ azimuthal (wall parallel, $\textit{n}=\text{3}$ ).

The most noticeable deviation from equilibrium turbulence characteristics at the start of convergence (Fig. 15(a)) is an increase in the near-wall $\overline {{u^2}} $ peak value compared to the equilibrium level (by ~85%, Figs 10 and 15(a)). This response is consistent with the measurements of Bourassa and Thomas^{(Reference Bourassa and Thomas45)} for a highly accelerated turbulent boundary layer in a plane convergent duct. Their data showed that as the flow approached the convergence the near wall axial stress increased by ~70%, agreeing well with Fig. 15(a). The explanation is that a favourable pressure gradient is also induced in the constant area duct upstream of convergence (see Fig. 14). Although in this region K remains comparatively small (here as in the Bourassa and Thomas^{(Reference Bourassa and Thomas45)} case only ~1.0 × 10⁻⁶), thinning of the boundary layer is still observed. The extra strain rate associated with this then leads to enhanced turbulence production of $\overline {{u^2}} $ . Figure 15(a) shows this effect is restricted to the near wall region (as also in Bourassa and Thomas^{(Reference Bourassa and Thomas45)}). Both SGS models produce similar $\overline {{u^2}} $ profiles; the models differ only in the response of the wall parallel stress $\overline {{w^2}} $ , where the Smagorinski prediction shows an essentially unchanged $\overline {{w^2}} $ peak, whereas the P + G model produces an extra peak very close to the wall.

The flow experiences much stronger acceleration further downstream in the convergence section. The expected re-laminarisation process will lead to diminished turbulence levels. The strength of reduction due to flow acceleration may be estimated by examining the Piomelli and Yuan^{(Reference Piomelli and Yuan15)} data, where the peak K value (4.0 × 10⁻⁶) was similar to that in the current nozzle flow. Comparison of Figs. 10 and 15(b) shows that non-dimensional turbulence levels have certainly decreased in the nozzle flow. More detailed comparison, however, reveals significant differences to the simulation results of Ref. (Reference Piomelli and Yuan15). Using peak turbulence levels as a measure, the present flow shows maximum values of $\overline {{u^2}} /\overline {{v^2}} /\overline {{w^2}} $ have reduced to around 20%/75%/80% of their inlet values. This is in stark contrast to behaviour reported in Ref. (Reference Piomelli and Yuan15), where these numbers were around 95%/15%/15%. Individual normal stresses in the current flow are clearly behaving quite differently than expected purely due to favourable pressure gradient effects. The nozzle flow predictions are essentially the same with both SGS models, although a slightly lower $\overline {{u^2}} $ peak and higher $\overline {{w^2}} $ peak are observed in the Smagorinski prediction.

It is clear that some other factor in addition to bulk flow acceleration is active in the current flow. The study reported in Ref. (Reference Piomelli and Yuan15) focused on the isolated imposition of a favourable pressure gradient on a boundary layer developing on a flat surface, whereas the present flow has a different geometrical arrangement and in particular includes a corner between convergence and parallel extension sections. It is possible – as noted in Ref. (Reference Trumper, Behrouzi and McGuirk29) – that corner flow effects have exerted a strong influence on the near-wall flow and turbulence development. The observations made on Fig. 15(b) are consistent with earlier comments where attention was drawn to a small region of high velocity in the corner location (the near wall region of high velocity at x/D = 3.65 in Fig. 14). This suggests that near-wall flow development has been modified both upstream and downstream of the corner. The upstream effect is shown in Fig. 16, which presents a predicted time-averaged axial velocity profile just upstream of the corner location (from both SGS models and normalised by U_CL). A local maximum very close to the wall and a wall-jet like profile in the region within 10% of the wall have appeared. The near wall flow structure no longer resembles a boundary layer-like shape. Both SGS models display the same flow feature; the P+G SGS model indicates a substantially stronger acceleration, whilst the Smagorinski SGS result has a similar profile shape but a smaller peak velocity.

Fig. 16. LES R²M predicted mean axial velocity profiles near convergence end.

Fig. 17. Near-wall mean axial velocity contours, start of nozzle exit parallel section, (top) LES R²M (P+G SGS), (bottom) LES R²M (Smagorinski SGS).

To examine effects downstream of the corner, Fig. 17 presents a zoomed-in view of a region very close to the wall and immediately downstream (radially only 0.004D (~0.25mm) and axially only 0.2D (~12mm). The P+G SGS model reveals the appearance of a small recirculation region, whereas in the Smagorinski SGS simulation no negative velocity is apparent, rather just enhanced growth of the near wall low velocity region. In fact, closer examination revealed a single row of cells immediately next to the wall (the region x/D = 3.65–3.7) does contain very small velocities (U/U₀ < 1%). Further inspection of the instantaneous Smagorinsky SGS time-series in this region showed this was characterised by generally positive axial velocity with occasional bursts of negative velocity.

Fig. 18. Near-wall resolved axial normal stress, start of nozzle exit parallel section P+G SGS model (top); Smagorinski SGS model (bottom).

Clearly, the Smagorinsky solution avoids time-mean separation because its SGS eddy viscosity is sufficiently larger than in the P+G SGS model to promote more momentum diffusion into the near-wall region just after the sharp corner. Both flow patterns in Fig. 17 are consistent with the presence of a blockage, which accelerates (see Fig. 16) and deflects the upstream bulk flow radially inwards. The backflow region in the P+G SGS simulation causes stronger flow deviation than the weaker blockage in the Smagorinski SGS simulation and is consistent with the higher peak velocity in Fig. 16. The differing profile shapes in Fig. 16 imply stronger strain rates with the P+G SGS model and hence the cause of slightly different turbulence profiles between SGS model predictions at plane E, for example, higher $\overline {{u^2}} $ peak with the P+G SGS model.

These velocity differences inevitably cause different turbulence development in the parallel extension section, as illustrated in Fig. 18 via contours of resolved axial normal stress (for a similar near wall zone as Fig. 17). The higher strain rates associated with the presence of a recirculation zone in the P+G SGS model generate much higher levels of turbulence (and much closer to the wall) than in the Smagorinsky SGS simulation. The P + G peak value is in the shear layer on the edge of the recirculation zone, whereas in the Smagorinski calculation the turbulence level is highest in the central region of the thickening boundary layer and increases gradually downstream, reaching a peak close to nozzle exit. A high spot of turbulence exists close to the wall just downstream of the corner in the Smagorinski SGS results; this is evidence of the transitory separation described above. It is important to note the contour scale difference in Fig. 18, the peak value with the P+G SGS model is an order of magnitude greater than predicted by the Smagorinski SGS.

Given the large differences indicated in Fig. 18, it is inevitable that predicted turbulence conditions at nozzle exit differ markedly between the two SGS models (Fig. 19).

Fig. 19. LES R²M predictions of Reynolds normal stress profiles at nozzle exit Left – P+G GS model, Right – Smagorinski SGS model blue $-\overline{\textit{u}^{\text{2}}} $ axial ( $\textit{n}=\text{1}$ ), red $-\overline{\textit{v}^{\text{2}}} $ radial (wall normal $\textit{n}=\text{2} $ ), green $-\overline{\textit{w}^{\text{2}}}$ azimuthal (wall parallel, $\textit{n}=\text{3}$ ).

For both SGS models the axial $\overline {{u^2}} $ stress has returned to being the largest, although peak values are quite different between the two models (NB ${\overline {{u^2}} } {/{2U_e^2}}$ is plotted for this model). In the Smagorinsky SGS simulation the peak axial stress is almost the same as in the inlet boundary layer; in the P+G SGS simulation, however, it is twice that in the inlet profile. For both SGS models $\overline {{v^2}} $ and $\overline {{w^2}} $ stresses have essentially recovered to close to their inlet profile levels, although the extra near wall $\overline {{w^2}} $ peak noted for the P+G model in Fig. 15(a) has now returned. The development of turbulence in the parallel extension is clearly very different between the two SGS model, only comparison with nozzle exit measurements can determine which captures the flow physics better, and this is addressed in Section 4.3.

4.4 Flow conditions at nozzle exit and comparison with measurements

Table 3 first presents measured and predicted global boundary layer parameters at nozzle exit; at this level of integrated mean velocity properties, both SGS models provide solutions close to the measured data.

Table 3 Global boundary layer parameters at nozzle exit

Exit profiles for the mean velocity (in Clauser-chart format) are examined in more detail in Fig. 20, which compares experimental data, LES results from both SGS models and results from a low Re Launder-Sharma RANS prediction.

Fig. 20. Predictions of nozzle exit mean axial velocity profile.

Both LES predictions show close agreement with experimental data, with the P+G model slightly better close to the wall; a definite improvement is seen when LES is compared to low Re RANS predictions. Based on mean velocity accuracy alone, it is difficult to distinguish between the two SGS models.

Fig. 21. LES R²M predicted nozzle exit axial turbulence intensity compared to measurements. Top: Smagorinski SGS model. Bottom: Piomelli and Geurts SGS model.

A quite different picture emerges when exit turbulence profiles are examined. Figure 21 compares predicted axial rms using both SGS models against nozzle exit LDA measurements.

For the Smagorinsky model the measurements are significantly under-predicted – a peak value less than half that measured. The location of the peak stress is also substantially further away from the wall than in the measured data, with the profile from the simulation showing a much more rounded shape. All these aspects are considerably improved when the Piomelli and Guerts model results are examined. Data for two values of the constant C_PG are shown, covering the entire range explored during the optimisation study reported in Section 4.1. Both values show superior performance compared to the Smagorinsky results. The value of C_PG identified by optimising the prediction of boundary layer inlet integral length scale returns very close agreement with the LDA measurements, predicting the measured peak value, its location and the general profile shape well. The sensitivity to C_PG is in fact not very great and the different profile shapes between the SGS models is entirely consistent with the observations made when discussing Fig. 18. It is clear that the corner effects have played a dominant role in determining the exit turbulence characteristics. It is also notable that the level of turbulence in the core flow outside the boundary layer is not well predicted by any of the simulations. The probable cause for this is that when the target data for the nozzle inlet conditions were specified, emphasis was placed primarily on the high turbulence region inside the boundary layer. If the actual turbulence in the outer core flow of the Trumper et al.^{(Reference Trumper, Behrouzi and McGuirk29)} experiment had been known, a better fit to this part of the profile might have been achieved.