1 Introduction

The investigation and understanding of wall bounded turbulent flows is of great importance to both academic and industrial institutions, with many practical applications ranging from aerospace to nuclear power. The understanding of flow phenomena such as thermal fatigue in T-juctions [27], traverse vibration due to generation of Dean vortex [21] switching in pipe elbows [127] and the application of ribbed ducts to enhance heat transfer [64] have benefitted from CFD of one form or another for several decades now. A commonly used approach in industrial CFD involves solving the discretised Reynolds averaged Navier–Stokes (RANS) equations, with closure provided by a turbulence model. This approach allows for reduced grid resolution; particularly within the near wall regions. While this approach allows relatively cheap simulations, the accuracy limitations of RANS methods can’t be ignored.

Direct numerical simulation (DNS) entails resolution of the full bandwidth of turbulence, and thus eliminates the accuracy issues of RANS. However, this imposes substantial demands upon grid resolution and overall costs. Large eddy simulation (LES) alleviates some of the expense associated with DNS by removing the computation of the smallest, dissipative scales that account for 99% of the overall computation cost [74]. Sub Grid Scale (SGS) turbulence models, first introduced by Smagorinsky [89], account for the effect of the small scale structures whilst the large scales of motion are resolved. Despite this, sufficient computational power is still unavailable to match the high resolution demands required for wall resolved LES for complex cases at high Reynolds Numbers.

Spalart [91] suggests the readiness of fully resolved LES is still a considerable time away for simulations of complex geometry. Recent studies have suggested a deceleration in Moore’s law. Moreover Waldrop [119] states that the growth factor of Moore’s law has slowed greatly with current computing techniques reaching their maximal operating limits. This was more recently emphasised by Larsson and Wang [52] who, based on recent technological advancements, projected that LES will not, and simply cannot, replace RANS solely in the design process within the next 30 years. Instead Larsson suggests it is more appropriate to consider questioning how and where LES could be utilised in conjunction with RANS methods within the design process.

It is these considerations that have led to growing recognition for the coupling of RANS and LES, in attempt to unite the advantages of each method [14]. Models of this nature have been generalised as hybrid RANS LES Methods (HRLM), with one prominent subcategory being Embedded LES—where regions of LES are embedded within a broader RANS domain (Fig. 1)

Fig. 1
figure 1

Example of embedded LES zonal arrangement with RANS/LES coupling

Despite LES being deemed more accurate than RANS, it does not mean that this accuracy is guaranteed. The sensitivity of LES imposes a well defined QA criteria that stipulates a strong dependence on numerics and grid generation, and it is possible that a weakly prescribed LES can produce worse results than a standard RANS. This conditioning is arguably controlled by the expertise and experience of the user. The risks associated with low experience users could potentially be alleviated by adopting an autonomous/semi-autonomous CFD approach. and it is argued that such software could become as significant as the improvement of computational speeds within industrial design [52].

The purpose of the following paper is to highlight the challenges associated with modelling internal turbulent flows, notably turbo-machinery and nuclear components found within industry; and to explore the potential benefits that hybrid RANS–LES methods can provide. Through initially presenting the currently available approaches of hybrid methods whilst also addressing the notable limitations on interface coupling, the advantages of the hybrid approaches becomes clear, with a comprehensive review of their application to industry style cases being presented. Furthermore comments on the cost perspective of both LES and hybrid RANS–LES approaches is given in Sect. 4, which provides a distinct outlook on the benefits to be gained. The paper concludes by introducing two new approaches to CFD within industry, aimed at progressing the continuing development and readiness of high accuracy computational methods outside of academia.

2 Hybrid RANS–LES

For even relatively simple geometries, the use of LES for internal flow is restricted to moderate Reynolds numbers, due to the mesh requirements in the near wall region. HRLMs can generally be separated into two distinct classes: (1) Global and (2) Zonal and shall be briefly introduced in the following section. The reader is also briefly referred to Table 1 for a simple comparison between the two classes.

2.1 Global RANS–LES

Global approaches are based on the continuous treatment of flow variables at model defined interfaces [25] which are generally prescribed from the solution. These methods follow a unified approach and generally solve one set of transport equations throughout the entire domain. An alternate formulation for the turbulent length scale forces a switch between the RANS and LES, which typically occurs within the near wall regions.

Table 1 Advantages and disadvantages of global and zonal approaches to hybrid RANS–LES

A widely used global method is detached eddy simulation (DES), originally proposed two decades ago by Spalart [94]. DES is a single grid approach where an unsteady RANS simulation, originally the Spalart-Allmaras model [92] but later adapted for others such as \(k-\omega \,SST\) [59, 95] and the elliptic relaxation \(\varphi -f\), is solved globally with a switching term enabling the model to change to LES when the local turbulent length scale exceeds the local grid dimensions [94]. DES has been implemented with great success for a number of different applications, whilst generally performing better for flows dominated by large regions of separation. A wide range of successful studies employing DES have been reported in the literature, beyond the scope of the present work; instead the reader is encouraged to refer to the following examples [15, 56, 101, 113]. This approach does however suffer from a couple notable issues as described by Deck [23], Tucker and Davidson [105] and Spalart et al. [93] (amongst others). Most notably amongst these is Modelled Stress Depletion (MSD) and Grid Induced Separation (GIS), where inconsistent handling between the filtered LES and time-averaged stresses, which are assumed equal at the RANS – LES interface, can lead to a reduction in the skin friction and ultimately initiates premature separation. It is documented that this issue is typically a result of an unwanted switch to LES within an attached boundary layer due to local grid refinement.

The more recently developed Delayed DES (DDES) by Spalart et al. [93] and later Improved DDES (IDDES) of Shur et al. [86], use blending functions to prevent the early switching to LES in attached boundaries. Meanwhile IDDES further extends the model’s capabilities by combining the approach with WMLES to resolve the mismatch between the modelled and resolved log layers.

In addition to DES, there has been a number of global methods developed over the past decade which have no explicit definition of the RANS–LES interface, instead are based upon a modified RANS model that adapts to scale resolving simulation under specific conditions. Such examples of this class include the Partially Averaged Navier–Stokes (PANS) [30] and the Partially Integrated Transport Method (PITM) [10]. The former uses a sensor based upon the ratio of modelled to total turbulent kinetic energy to determine which turbulence solving mode the model acts. The PANS ratio was later improved by Basara et al. [8] to allow for run time adaptation.

An alternate approach to the general consensus of global methods has been proposed by Xiao and Jenny [125], in the form of a consistent dual-mesh hybrid RANS - LES, with the model more recently being generalised and applied to wall-bounded turbulent flow by Tunstall et al. [108]. Although not strictly of the natural form of the global methods the duel-mesh approach solves both RANS and LES models simultaneously on separate grids, with added drift terms in the momentum equation used to relax the flow statistics towards a consistent solution, or that which is the most relevant at a particular space and time. This framework alleviates the fundamental problems encountered along the interfaces in general global methods such as DES, allowing a clean solution between near-wall RANS and LES free-stream conditions. In a similar approach Uribe et al. [110], a two-velocity field splits the residual stress tensor into two parts, locally isotropic and inhomogeneous, with a blending function between the two.

Fig. 2
figure 2

Example of a two region embedded LES pipe simulation showing the transition between RANS and LES (left), using the field of 3-axis ellipsoid synthetic eddies randomly generated at the RANS–LES interface represented as a schematic (right)

2.2 Zonal RANS–LES

Restricting LES to a specific portion of the domain, whilst maintaining a RANS solution elsewhere offers potential for gains in efficiency. Generally, zonal approaches adopt segregated LES regions embedded within a RANS solution. Each subdomain is then solved individually with a separate set of transport equations. For example, to a two-region Embedded LES simulation is shown in Fig. 2; this figure illustrates the transition region between an upstream RANS and downstream LES through the injection of synthetic fluctuations. Unlike global methods, there is a clear transitional boundary between RANS and LES regions, although this imposes further complications which shall be considered later.

As an extension to DDES, Deck proposed a zonal formulation of the standard model [22, 23]. Zonal DES (ZDES) acted to utilise a-priori knowledge of the flow field in order to strategically prescribe specific RANS and DES regions in the domain. The method proved advantages for flows characterised by strongly developing turbulence inherited from upstream boundary layers. ZDES differs from Spalart’s DES through the definition of a ZDES length scale, sub-grid scale length and treatment of the near wall function in LES mode [50]. Improved formulations of ZDES were later released again by Deck [24, 26].

Poletto et al. [72] investigated the use of Embedded DDES regions whilst computing flow around a 2D hump, with an emphasis on a synthetic RANS to LES interface condition. They noted a considerable reduction in CPU time between an embedded DDES approach and a full domain DDES method could be achieved, with an overall 30% speed-up in CPU requirement per time step.

In addition Davidson and Peng have worked towards producing a PANS based ELES model [17, 18] again with an emphasis on the interface conditions, where Partially Averaged Navier–Stokes (PANS) is a modified \(k-\epsilon\) model that can operate in both RANS and LES modes. However the studies indicate that the results within the LES region are highly sensitive to the interface conditions and upstream turbulent quantities, with the recovery of the recirculation regions downstream of the hump flow slightly over-predicting the measured value.

An early implementation of ELES is presented by Cokljat et al. [14] who considered a fully developed pipe, channel flow and flow over a backwards-facing step. This work was one of the first to demonstrate the potential of ELES whilst offering notable advantages to industrial CFD. A number of more recent implementations have followed, most notably Li et al. al. [54] validated a two-region ELES solver that consisted of single RANS and implicit LES zones with the transfer of variables handled at the interface. They applied this implementation to both channel flow and flow over periodic hills noting good agreement to full LES data. Anupindi and Sandberg [2] implemented an ELES model within OpenFOAM based again on separate regions for the RANS and LES zones and the transfer of flow variables via the interfaces. In a slightly different approach, Vonlanthen et al. [118] produced a one-way nesting procedure which embeds a highly resolved LES within a low resolution LES.

The original concept for embedded simulation arose broadly from the need of LES users to define unsteady inlet boundary conditions from mean flow data, but the challenges facing a more general use of ELES are significantly greater. In theory the LES region itself could be entirely engulfed within a global RANS region, as depicted in Fig. 1. For instance a cuboidal LES region with 6 sides will have a different combination of inlet and outlets depending on the flow direction, and in more complex flows, these definitions may change during the simulation; it is indeed quite likely that instantaneous flow will travel in both directions across any boundary [4]. It is these interfacing features which represent the majority of new developments in the community, along with the development of best practice guidelines to facilitate it’s usage [60].

2.2.1 RANS to LES Interface

Tabor and Baba-Ahmadi [96] presented a comprehensive review on the subject of inlet conditions for LES, separating the methods into two distinct categories; (1) Precursor/recycling methods and (2) Synthesised Methods. The former involves computing a precursor simulation either prior to or concurrently with the target LES to generate a library of turbulent content that can then be mapped, and scaled if required, onto the inlet of the LES domain. In contrast, synthesised methods involve the injection of synthetic fluctuations upon a prescribed mean flow. More recently Wu [124] provided an updated review on LES inlet techniques, whilst also beginning to describe their application within an hybrid RANS–LES framework. Following the classification proposed by Tabor, Wu highlights that research has led to the emergence of several synthetic techniques. Wu argues that recycling methods are less efficient than the synthetic counterparts and despite being highly accurate their additional costs and lack of generality restricts use with hybrid models [39].

The easiest approach would be to inject white noise upon the mean flow of the upstream RANS. Although simple this method fails to provide the spatial or temporal correlations required to sustain turbulent generation downstream, as shown by Schluter et al. [78] and illustrated in Fig. 3 (top). Here the generated velocity field disappears very quickly with the flow immediately becoming laminar.

In attempt to alleviate the lack of spatial and temporal coherence Jarrin et al. [38] presented the Synthetic Eddy Method (SEM). SEM is derived from the vortex method originally outlined in a PhD thesis by Sergent [81]. Unlike white noise, SEM generates a fluctuation field through the superposition of coherent structures across an allotted plane, each with a prescribed weighted velocity distribution. The generated velocity fluctuations are normalised to reproduce the required second order turbulent statistics, based upon an upstream RANS field, via a Cholesky decomposition of the Reynolds stress tensor.

Improvements to the SEM formulation have since materialised, aiming to reduce the recovery/development length of the resolved turbulence whilst potentially alleviating some of the computational costs. Firstly an improved version by Adamian et al. [1] documented a decrease in error of the averaged flow parameters whilst also shortening the transition region, through a modification of the determining of the linear scale of generated eddy structures. Meanwhile Skillen et al. [88] proposed an alternate fluctuation normalisation factor, based upon a running average of the eddy concentration that guarantees the desired statistical properties, regardless of the spatial distribution or length-scale of the eddies, thus correctly allowing for an inhomogeneous distribution of eddy size. Figure 3 (bottom) illustrates the generated velocity field within OpenFOAM acting on a straight pipe flow. The computation here is one-way coupled such that the instantaneous turbulent fluctuations are superimposed onto a mean velocity profile obtained from an upstream Reynolds Stress Model (RSM).

Fig. 3
figure 3

Comparison between the flow development within an LES pipe using different synthetic inlet conditions, images taken from a slice of both the inlet plane and stream-wise centreline. (top) White noise inlet. (bottom) Improved synthetic Eddy method

Pamies et al. [66] also provided subsequent improvements to the original SEM formulation by modifying the eddy shape function to take into account the variance in structures across the turbulent boundary layers observed in experiments. Each randomly generated eddies’ shape, length and time scales are dependable on their ’altitude’ or distance from the wall, and are adjusted depending on the specification of its given mode value.

Meanwhile the most recent improvements to SEM incorporates the use of anisotropic forcing of the Reynolds Stresses within an overlapping region of RANS and LES. This method was originally outlined in a PhD thesis by De Meux et al. [19] although the reader is also directed to a later publication by the same author [20], both documenting a decrease in development length and a faster recovery of the turbulent statistics.

Despite positive results perhaps the most documented inherent weakness of the various forms of SEM is that the generated velocity field does not satisfy the divergence-free condition, leading to an introduction of significant pressure fluctuations at the inlet plane [71]. This can have a detrimental effect upon the development length, and expectedly pressure fluctuations are unwanted when one is concerned with noise prediction. Divergence Free SEM (DFSEM) is an alternative version to the standard approach that proposed a change in definition of the velocity fluctuations associated with the eddies, to provide a divergence free velocity field. Poletto et al. [73] reports the development of DFSEM, and is largely based on the original methodology by Jarrin et al. [38], however instead of applying the fluctuations to the velocity field directly they are first applied to the vorticity field and then transformed back to velocity by taking the curl of the vorticity.

In a separate category to the synthetic eddy variants there exists a range of synthetic turbulence models based upon the superposition of weighted spatiotemporal Fourier model, reviewed in detail by Tabor and Baba-Ahmadi [96]. A notable example developed for specific use in hybrid RANS–LES is the NTS Synthetic Turbulence Generator (STG), documented by Shur et al. [87].

Similarly Auerswald and Bange [5] developed an anisotropic synthetic generator based on probability of Fourier modes. This approach was later extended by Auerswald et al. [6] to account for the integral length scales in y and z directions.

Additionally a divergence free approach to Fourier mode analysis has been published by Patruno and Ricci [69].

Despite the somewhat documented difficulty and importance of a well developed synthetic fluctuation method for the inlet to LES simulation (both coupled and uncoupled RANS–LES), the recent publications reviewed above are proving to be reasonbaly successful. Both the NTS STG and SEM variants have been proven to produce ’life-like’ turbulence consistent with the desired statistics, whilst also maintaining a short development/lag phase with minimal computational expense.

As the case with ELES it may be required that the computation switch back to RANS mode from an upstream LES. Thus introducing another interface with differing demands to those when passing from RANS to LES. A review and explanation of the methods involved is provided in the next section.

2.2.2 LES to RANS Interface

On returning to the RANS domain it is required that one recovers the mean flow data from the resolved turbulent descriptions produced from LES. Thus the interface must aim to imitate a statistical averaging operation, where the LES solution is averaged both spatially and temporally to remove the unnecessary turbulent data. Additionally the interface must ensure that the fluctuations from the LES solution leaving the LES domain do so without incorrectly producing turbulent anomalies in the RANS domain [46]. This averaging process takes place at the interface which itself houses an overlapping region between the two meshes [79], this enables for flow variables to spatially adapt to the model change, whilst trilinear interpolation is used to exchange flow quantities between meshes.

However the process is not as trivial when reestablishing the modelled turbulent quantities, since these don’t appear within the LES formulation. Instead values such as eddy viscosity, turbulent kinetic energy and turbulent frequency have to be reconstructed. Konig et al. [46] introduced three different methods to accomplish this, testing each case upon a turbulent channel flow. Each method involves calculating the turbulent kinetic energy directly from the LES solution as, \(k=\frac{1}{2}\bar{u_i''u_i''}\) for the turbulent viscosity, method A and B firstly determines the turbulent dissipation, \(\epsilon\) from the transport equation of the turbulent kinetic energy with method A taking into account both the resolved and unresolved scales whilst method B neglects the unresolved motions. Method C takes a slightly different approach, and instead calculated the turbulent viscosity via computing the turbulent frequency, \(\omega\). This approach (method C in particular) has recently been utilised by Anupindi and Sandberg [2] for use in an ELES solver implemented in OpenFOAM with reasonable success documented. Meanwhile Gritskevich et al. [34] document the results of two different approaches to the recovery of the mean flow data when passing from LES to RANS. The first, labelled a 1 stage approach continuously solves the chosen RANS model throughout the whole computational domain. However, within the LES region the RANS model runs passively in the background using the LES flow field but having no effect upon the solution; upon entering the RANS domain the solution to the RANS model is reinstated. The second or 2 stage approach initially computes the whole domain with the RANS model, the solution for the turbulent variables i.e. turbulent kinetic energy (k), turbulent dissipation (\(\epsilon\)) are then frozen within the LES domain during the second computation. This frozen solution is then re-activated at the inlet to the RANS domain.

3 HRLM Applied to Internal Flows

In the following section we review examples of (1) the benefits of using LES over RANS, and (2) the application of hybrid RANS–LES methods, both in the context of internal flows. All the papers discussed in this section are summarised in Table  2 for clarity and to facilitate comparison.

Table 2 Summary of industrial relevant flow applications reviewed in this paper and the modelling method applied, also shown are the grid densities and Reynolds number used for each case respectively

3.1 Turbo-Machinery

3.1.1 Internal Cooling Ducts

Turbine blades are generally operating at temperatures above the operational limit of the blade material, and small fluctuations above the design temperature can reduce a systems operational life by 50% [111]. These facts would suggest that the cooling mechanism employed is of the upmost importance. Cooling is achieved through internal ducts housed within the blade and employ forced convection through induced separation, promoting highly turbulent/unsteady flow [114]. The flow characteristics and prediction of heat transfer have been extensively studied on these geometries using RANS, LES and increasingly hybrid RANS–LES. In addition the angle alignment and size of the ribs have also been altered in attempt to assess the variation in heat transfer capabilities.

Table 2 documents a number of the reviewed application computations, evidently many are available for internal cooling systems/ribbed ducts. The RANS studies listed report incorrect predictions of the Nusselt number [90] demonstrating a variation of around 250% [102]. Tucker [104] provided a compilation of some of the more notable RANS models’ prediction of Nusselt number distribution across a ribbed channel, see Fig. 4. In this case it is likely that diffusion terms in the RANS models would act to suppress the high levels of near wall turbulence, having a direct impact on the heat transfer in this region.

Fig. 4
figure 4

Nusselt number, Nu distribution and subsequent heat transfer capability across a ribbed channel for various RANS models reproduced from [56, 104] (left), and instantaneous streamlines in side and top down views for Zonal LES computation from [56] (right)

The relatively poor performance of RANS models is not entirely unexpected, particularly their application to complex three-dimensional structures or regions of flow separation. Heat transfer is directly linked to the structure of the flow, and so if the turbulence is not correctly captured then one cannot expect to calculate the correct levels of heat transfer.

Increased accuracy of these flows are achieved when adopting LES. Tafti [97] conducted a computation of a periodic ribbed duct at a Reynolds number 20,000, using a relatively large grid resolution both the heat transfer and skin friction coefficient was predicted within 5% of experimental values. However the computation is then much more costly with almost 20 hours required to compute a single flow through.

In addition, Labbe [49] (conducting a study based on the experimental set-up by Casarsa et al. [9]) noted that the LES flow filed is strongly linked to the inlet conditions applied. Here slight discrepancies were found due to a poor handling of the inlet.

As yet there are, to our knowledge, no published documents using ELES for an internal cooling duct. Instead most have chosen to adopt global hybrid methods, most notably DES. However in a similar context Jorgensen et al. [42] applied ELES to flow around a wall-bounded cube. A schematic of the domain is shown in Fig. 5. Synthetic fluctuations in the form of the 2D vortex method of Mathey al el. [58] were applied at the inlet to the LES domain. The mean and standard deviations of the stream-wise flow component showing good agreement with experimental data and also the body induced turbulence being well captured, despite this the simulations experienced some difficulties due to poor control of the upstream turbulent structures.

In order to mitigate the high costs reported above, Viswanathan and Tafti [113] adopted DES on the same domain used by Tafti [97] achieving only a 9% deviation with respect to fine grid LES at 1 / 10th of the cost. Additionally several studies by Vishwanathan and Tafti [112, 115, 117], Liu et al. [56] and Kubacki et al. [47] successfully implemented hybrid based methods to predict the turbulent and mean flow features in stationary, rotating and roughened ribbed ducts – all documenting reasonable success.

Again Tucker [104] has compiled a comparative graph of the Nusselt number development across a single rib for a number of turbulence models, including hybrid RANS–LES show in Fig. 4 (right). It is shown that the hybrid models provide a much more favourable result than any RANS alternative, despite slight over-predictions estimated by ZLES after re-attachment of the flow.

One could argue that global hybrid methods may be more suitable for flows of this nature given the proximity of subsequent ribs, where a clear interface cannot be defined for ELES as the high accuracy detail is required throughout the whole domain, rather than within a particular confinement of the geometry.

Fig. 5
figure 5

(top) Schematic of the embedded LES domain for the floor-mounted cube used by Jorgensen et al. [42]. (bottom) Coherent small scale vortex structures generated within the LES domain [42]

3.1.2 Two-Pass Cooling Duct

This case is an extension of the the single periodic ribbed duct flow including two channels connected via a \(180^0\) bend. Sewall and Tafti [82] conducted an LES study of a stationary two pass duct, with and without a rib in the bend, whilst Viswanathan et al. [116] conducted the same simulation using a \(k-\omega\) based DES model. The former LES study required a grid resolution of 0.7 million cells greater than the DES of Viswanathan et al. nearly an order of magnitude less than Sewall’s simulation. Despite the DES computation over-predicting the flows developments length at the inlet there is a good quantitate comparison with both LES and experimental results.

3.2 Nuclear Applications

The development and operation of nuclear reactors demand accurate understanding of the complex thermal hydraulics involved in order to ensure rigorous safety regulations are satisfied. For example, excessive thermal fatigue resulting from the turbulent mixing of fluids at differing temperatures can lead to a complete destruction of the system [12]. It is argued that thermal fatigue is the most common cause of unexpected failures within Nuclear Power Plants (NPP’s) particularly within piping components such as mixing tee’s [29]. Flows throughout complex piping geometries are highly unsteady and involve broad bandwidth, low frequency fluctuations that are poorly represented through standard URANS [31]. Westin et al. [122] and Walker et al. [120] illustrated that both steady and unsteady RANS calculations are unable predict realistic mixing between two flows [122].

Additional structural vibration and fluid-born noise can be generated in pipe regions involving sharp bends due to the presence of large pressure gradients [129]. The capturing of these secondary flow features has been studied extensively by Takamura et al. [98], Yuki et al. [127], Tanaka and Ohshima [99] amongst others for the Japan Sodium cooled Fast Reactor (JSFR).

The complex and transient characteristics of the flow regimes found in the nuclear industry exemplifies the requirement for higher order turbulence models within constraints the constraints of resources.

3.2.1 Thermal Mixing in a T-Junction

A number of reports have successfully adopted LES to the problem of thermal mixing in a T-junction. Lee et al. [53] studied the temperature fluctuations and structural response of a standard mixing tee configuration by coupling LES with a newly developed thermal stress fatigue analyses. Meanwhile Ndombo and Howard [62] coupled LES with the inclusion of synthetic turbulence at the inlet. In addition to standard T-junction flows model multiple studies have also looked at the effects of an upstream pipe bend prior to the mixing zone. Pasutto et al. [68], Aulery et al. [7] and Tunstll et al. [107] successfully applied LES to predict the effects of upstream bends on the mixing of the two flow streams. All showed good agreement to experimental data whilst the later also adopted the use of DFSEM at the inlet to the two pipe streams.

Despite the obvious accuracy advantages provided from LES models their applicability within an industrial framework is limited. The computations disclosed in the previous paragraph were expensive and required access to significant High Performance Computing (HPC) resource. It can be estimated that an LES simulation for a straight inlet T-Junction at a moderate Reynolds number of \(10^4\) with a maximum grid resolution of 20 million required, at minimum, \(>300\) cpu-hours [63] to compute; with this increasing further if higher order statistics are desired. In contrast a RANS computation could be completed in less than one day, whilst still providing accurate predictions of the time-averaged thermal-hydraulic characteristics [55]. It is not uncommon for flows in nuclear applications to reach Reynolds numbers of \(10^6\) or \(10^7\), and so measures must be taken to reduce the cost overheads whilst maintaining a noteworthy degree of accuracy.

Gritskevich et al. [31] investigated the application of hybrid RANS–LES approaches for predicting heat transfer at the walls in a T-junction’s mixing zone. Of particular interest here is the success of ELES, with Fig. 6 illustrating a typical domain decomposition for ELES in a T-Junction. In Gritskevich’s computation there are two RANS to WMLES interfaces, both have been treated with the vortex method, whilst a single WMLES to RANS interface downstream of the mixing zone has no additional treatments.

Fig. 6
figure 6

(top) Decomposition of computational domain for ELES T-junction calculations. (bottom) Iso-surfaces coloured with velocity, with BCD boundary conditions [31]

It is suggested by Gritskevich et al. [31] that the use of global schemes, such as Scale Adaptive Simulation (SAS) and DDES, could be problematic for the computation of flows where there is no natural geometrical or flow induced conversion to Scale Resolving Simulation (SRS) – i.e. there is no seperated flow. This issue is highlighted by Deck et al. [26] who describes the different modes of separated flows for DES and uses this as a justification for the development of ZDES. More concisely, global models will suffer greatly if the initial instabilities of the flow aren’t sufficiently strong or, if there is little variation in turbulent length and time scales across the domain. Menter et al. [60] has previously emphasised this, stating that flows deemed to be stable are characterised by a continuous development of turbulence and are largely dependent upon the upstream flow conditions; this is typical of majority of channel and pipe networks, where the flow is likely to remain attached over much of the domain. The use of Embedded LES is therefore more promising than global type models for these types of flow. The rigid interface to the LES allows for the injection of turbulent fluctuations to encourage a change to SRS, whereas there exists an ambiguity in the transfer to SRS in DES variants.

Turbulent instabilities can be subdued by numerical dissipation, and so a poor choice of convection scheme can have a particularly damning effect on the level of unsteadiness in the flow. This is apparent in Fig. 6, where the dissipative nature of the second order Boundary Central Difference (BCD) scheme [40] impedes upon the formation of turbulence on both global models. However, the most striking effect occurs with SAS. A high sensitivity to the choice of velocity interpolation scheme causes a conversion back to a URANS formulation shortly downstream of the junction. This subsequently leads to a drastic under-representation of the RMS velocity magnitudes. On the contrary, the ELES result appears to be little affected by the convection scheme, and instead the introduction of resolved turbulence at the RANS–LES boundaries is accounted for downstream at the formation of new turbulence in the mixing zone.

The work by Gritskevich et al. [31] is later extended to include the combination of ELES with IDDES and SAS models so that the previous WMLES zones are occupied with either IDDES or SAS [33]; additionally the RANS to SRS interface is computed with both vortex method and a Generator of Synthetic Turbulence (GST) [1]. Each method was able to predict the flow and thermal mixing with reasonable accuracy with respect to experimental data, however they appear to slightly under perform when compared to the ELES case in Gritskevich et al. [31], which employed a WMLES model, albeit at additional cost. Overall the findings of the two papers for ELES is encouraging.

3.3 Capturing Secondary Flow Structures in Pipe Bends

T-junctions represent only one component of an entire complex piping system of a nuclear power plan or gas networks. Pipe systems commonly include a number of 90 degree bends to redirect the flow or to artificially enhance the erosion process in helium stratification [45, 77]. Although common, the flow characteristics of a 90 degree bend typically includes the pressure-induced separation at the apex, which is shown to be sensitive to the location of the upstream conditions. The separation can lead to a strong unsteady shear layer and a significant momentum deficit [77]. The flow exhibits several other distinguishing features, including: longitudinal streamline curvature, velocity profile inhomogeneity and development of secondary, swirling motions named Dean vortices [21]. These features present a significant challenge to RANS models.

Fig. 7
figure 7

Schematic of the 90 degree pipe domain of Rohrig et al. [77] (top left) and snapshot of the computational domain both across the bend and cross-section views (top right). Results for reduced inlet size of LES region (btm left) and instantaneous velocity for case with 1D upstream inlet section (btm right)

Despite this, the study of flow around a 90 degree pipe elbows has predominantly adopted stead-state RANS models. A few recent publications such as Rohrig et al. [77] have begun to embrace the application of scale resolving methods such as LES. Rohrig’s study involved the benchmarking of a number of RANS models against LES for the domain seen in Fig. 7. Ultimately the study confirmed the overwhelming superiority of LES. It is found that the selected RANS models fail to correctly capture the mean velocity profile within the momentum deficit region of the bend, resulting in a underestimation in wall shear stress. It should however be noted that more favourable performance is possible when using more advances RANS models, such as the more recently developed Elliptic Blending Reynolds Stress Model (EBRSM) [51], as evidenced in the work by Tunstall et al. [106]. That said, significant deviations in the mean flow profiles remain even with EBRSM.

If employing LES for this flow, for example, at a reasonable Reynolds number of \(10^4\), one requires approx. 20 million cells [77] to compute the flow around the bend, whereas only 1.5 million might be required to obtain a grid-converged RANS solution [106]. Another key point to note is the inlet condition. For the LES by Rohrig a fully developed profile at the inlet was applied through an a-priori LES of a straight pipe section, which itself required a computational domain of 5 million cells. The LES requires HPC to compute. Thus the costs associated with computing even a simple pipe elbow with LES is formidable for industry.

An possible alternative here is Hybrid RANS–LES or WMLES, as Tunstall et al. [106] also presented. An example of ELES was presented by Holgate et al. [35], where a reduced upstream domain size was presented as shown in Fig. 7. It can be seen that results even with an inlet length equal to 1 diameter provide good results; significantly reducing the cost of LES simulation for this flow and removing the requirement of a precursor computation entirely.

4 Cost Perspective and Emerging Use Modes

It is possible to identify an order-of-magnitude cost associated with internal flow simulations in terms of the spatial and temporal resolution requirements, and extrapolate this as a function of the flow Reynolds number. Figure 8 is an adaptation of cost estimates for large eddy simulation based on approximations proposed by Piomelli and Balaras [70] and Tucker [104]; largely based on the findings of Chapman [11]. The dashed line represents the LES grid estimation and is suggested to scale with \(Re^{2/5}\). Meanwhile the dash-dot line, scaled with \(Re^{9/5}\), represents the wall-modelled LES with resolution down to \(y+<100\). A more recent paper by Choi and Moin [13], who revisited the work by Chapman, suggests that theses relationships are much more severe i.e. scaling with \(Re^{13/7}\) and Re for wall-resolved and wall-modelled LES respectively. Approximations therein provide a broad understanding of the grid size N required for conventional finite-difference/finite-volume approaches for a range of Reynolds numbers. While loosely factored into the formulation, elements such as geometry complexity and physical domain size are not explicitly considered, and thus the estimate is taken only as a rough rule of thumb. A comparison of computational requirement for a pipe flow at a Reynolds number of around \(10^6\) is provided in Table 3.

Table 3 Estimation of CPU resources required for RANS, hybrid and wall resolved LES for internal flow applications i.e. pipe flow at \(Re\approx 10^6\)
Fig. 8
figure 8

Grid requirements for LES, DES and ELES with differing applications listed. Costs estimated based upon usage of UKs largest HPC facility as a non-contributing member. T = T-junction; RD = ribbed duct; Ch = channel flow

As discussed in the previous section it is not uncommon for the Reynolds numbers to reach \(10^6\)\(10^7\) for industrial cases. This would suggest that grid sizes for LES computations would exceed 1 billion cells in the higher Reynolds range, which is beyond that generally available for research studies, let alone industrial work. Clearly then, this cost could be reduced substantially if a hybrid scheme was adopted, for example with a Re of \(10^7\) the estimated number of nodes for WMLES is between \(10^7\) and \(10^8\), placing the notional computational expenses within a more realistic region. Also listed in Fig. 8 are the estimated notional cost of calculations on ARCHER, the UK’s national HPC facility,Footnote 1 for completing a computation with respect to the corresponding grid resolution.

A chapter by Fureby written in ‘Whither Turbulence and Big Data in the 21st Century?’ [28] discusses the increasing need for LES in engineering flows, but highlights the challenge posed by the generation of extremely large data-sets such as those from LES/DNS of high Reynolds number flows. As an example, consider the pipe bend simulation of Rohrig et al. [77], reviewed in Sect. 3.3. Using 23 million cells the resulting memory required to store both the flow (up and k) and statistical variables is around 2.3 GB per time step. A typical time step is 1.5 ms and 6000 time-steps required per flow through time unit, with over 30 such units needed to obtain time-averaged data; which translates to around \(1.8 \times 10^{5}\) time-steps. For memory and running averaging purposes each case is run with a purged time-step that provides 20 time-steps to be stored. In this case, a total memory of around 46 GB purely for the data would be required, (plus additional storage for the mesh and case files). This storage requirement would be expected to increase dramatically upon post-processing, for instance if high quality images and animations are required. Delving deeper into the running of this simulation, it was found that on 18 nodes (432 processors, \(\approx 55000\) cells per processor), \(\approx 5s\) real-time was required per time-step. Therefore to complete the case and allow for higher order statistics to be correctly averaged then a total run-time would exceed 10 days, costing in excess of 1.5 MAUs. Such resources will remain beyond the reach of a large section of industrial CFD users for some time to come.

A more recent alternative to in-house commercial CFD involves the outsourcing of larger computations to online or cloud based services. Wu et al. [123] illustrates the advantages gained from using cloud based environments. A good overview on the definition, characteristic and overall potential impact of cloud computing can be found in a publication by Shawish and Salama [85].

Fig. 9
figure 9

Flowchart for automating the embedded simulation process

4.1 Automated Embedded Simulation

In light of what has been presented in this paper, and more importantly the information learned from previous studies surrounding the application of hybrid RANS–LES models, a workflow can be proposed for Embedded Simulation. A semi-automating capability is highly attractive for industry, based on the use of RANS as a precursor to ELES, where the users input and required knowledge is substantially reduced. Not only does this enable reduced workflow times, but it also increases simulation repeatability, which in many ways is of greater importance than overall simulation accuracy; since in many industrial CFD studies it is the qualitative information which is important.

Figure 9 illustrates the process for which the automated system could follow. This process automatically implements the general capabilities outlined by Cokljat et al. [14] for which an ELES model is required to provide in an industrial environment. A pre-cursor RANS calculation allows for sufficient knowledge of the flow-field to be obtained, so that a suitable LES mesh can be generated and the extent of regions requiring further detailed simulation can be identified. The overall aim of automation is to allow for a streamlined process chain for computing internal fluid flow. Such a process remains to be clearly defined although preliminary work by Prost et al. [76] has demonstrated the concept of automated RANS/LES interface control by using algebraic sensors to and distinguish between attached and detached flow regions.

4.2 Multi-Dimensional Coupling

Finally we note that in many industrial sectors the need to retain low spatial dimension simulations will likely remain for some time. This is particularly relevant for systems where the flow length scale, e.g. pipe diameter, is much smaller than the overall system dimension e.g. pipe length. For example, gas network companies who consider flow analysis over vast lengths of pipelines simply cannot afford full 3D or even 2D simulation of the whole system; and in most instances 1D modelling is sufficient. However there are invariably instances where more detailed simulation is required and hence the industry must step from 1D to 3D analysis. The idea of multi-dimensional coupling involves a continuous workflow or seamless transition between differing physical platforms could be achieved. This is portrayed in Fig. 10 which illustrates a transition from fast 1D interactive tools to high-detail 3D complex physics tools. For the analysis of pipe networks the 1D studies can be used to model flow through an expansive network, and to quickly cycle through a range of scenarios in order to identify the critical cases. A subset of these cases could then be studied in 3D, combining RANS with LES where needed; for example to provide unsteady 3D insight to a range of piecewise components or larger combined scenarios (such as bends, valves and compressors) to examine the detailed impact of surges through precise representation of local pipework and equipment. For more details see [35].

This concept is already receiving interest. A recently coupled 1D - 3D solver implemented in OpenFOAM uses Riemann invariants to ensure the 1D boundaries are updated each time step from the computed 3D region, whilst the 3D inlet boundaries apply the upstream 1D pressure head for the pressure whilst the velocity is set from a homogenous Neumann condition; this method was successfully applied to transient pipe flows with only minor losses noted [121]. In addition, an earlier publication by Prince [75] focussed on the coupling of 1D and 3D systems for computing flows through subway transit networks.

Fig. 10
figure 10

Concept for multi-dimensional coupling process for pipe networks (top), mesh used for 1D–3D of water hammer in pipes from  [121]

5 Conclusion

It has been demonstrated that hybrid RANS–LES approaches have a role to play within industry to facilitate scale resolution of turbulent flows. Furthermore the benefits of embedded LES, where reduced computational costs are clear compared to full LES, whilst providing an equivalent level of fidelity within chosen regions of the domain. This paper has established the feasibility of applying hybrid RANS–LES methods to industrial CFD applications, whilst documenting the separate ingredients involved within the approaches. It is clear that the use of a hybrid solver would allow for a reduction in total grid resolution over traditional LES, thus reducing both the computational expense and simulation time, whilst also potentially increasing the overall simulation accuracy when compared to RANS.

Future trends have been discussed and greater efficiency and repeatability could be attained by automating the computational process when introducing the two concepts outlined above. For example by reducing user input required to generate suitable meshes for turbulence simulation and making use of synthetic turbulence generators to pass from RANS to LES regions. The development of multi-dimensional coupling methods is particularly useful for industries which consider pipelines or systems which operate over a significant separation of scales.