On the accuracy of partially averaged Navier–Stokes resolution estimates

Partially Averaged Navier–Stokes computations can employ three different approaches for specifying the ratio of modelled-to-total turbulence kinetic energy fk. Use can be made of either a constant, a spatiallyor a spatiallyand temporally-varying value. This work compares different estimates for fk found in literature and evaluates them for two test-cases: a circular cylinder at = Re 3900 and a turbulent channel flow at = Re 395. Additionally, the estimates are compared to the a posteriori computed ratio of modelled-to-total turbulence kinetic energy ( f̃k) obtained from the PANS flow solution. The trends observed for the estimates are similar, although the magnitude varies significantly. All spatially varying fk approaches reduce the PANS model to a DES-like model, thereby entangling modelling and discretisation errors. At the same time, f̃k shows that the behaviour of these estimates is incorrect: fk becomes too large near the wall of the object and in the far field. It is observed that f̃k is always lower than the set value, when using fk fixed in space and time. Finally, it is clear that the estimates, applied to internal, boundary layer, flows yield too high values for fk. In order to minimise errors and increase the reliability of industrial CFD results, the approach with a constant fk is still preferable, assuming suitably fine grids are used.


Introduction
Scale-Resolving Simulations (SRS) allow the user to obtain a more accurate description of the flow compared to Reynolds Averaged Navier-Stokes (RANS) methods, by resolving part of the turbulence kinetic energy spectrum. At the same time, the computational cost is significantly lower than for Direct Numerical Simulation (DNS) which makes SRS attractive for industrial applications. SRS can be divided into three main categories: (1) Large Eddy Simulation (LES), where an (implicit or explicit) filter is applied throughout the entire domain, such that the large scale turbulent structures are resolved and the smaller scales are modelled using a 'sub-filter' model. The need to resolve a substantial part of the turbulence spectrum everywhere leads to excessive computational cost for industrial CFD. (2) 'Hybrid' methods, where RANS is used in near wall regions, and Large Eddy Simulation (LES) in the far field; and (3) 'Bridging' methods, such as Partially Averaged Navier-Stokes (PANS) (Girimaji and Abdol-Hamid, 2005). These Bridging methods consist of a blending of RANS and DNS, yet in contrast to Hybrid methods the blending is not location dependent. The advantage of a Bridging method is that the blending depends on user defined settings, and it allows a smooth transition between the turbulence modelling approaches. Bridging methods with a constant physical resolution, also do not suffer from commutation errors, which affect Hybrid methods due to the flow switching from RANS to LES regions. An advantage of PANS over LES is that the model offers a separation between modelling and discretisation error. This enables Verification and Validation processes, which are essential for the credibility of CFD solutions. In the particular case of the PANS model the blending between DNS and RANS depends on the modelled-to-total ratios of turbulence kinetic energy, f k , and dissipation, f ϵ .
For the usage of PANS two approaches can be distinguished: (1) the Constant f k approach, where a constant value of f k is used in the domain and throughout the simulation time. This approach was often used to verify the PANS model, but has mostly fallen out of favour recently since it is more computationally expensive in cases with a large range of different turbulent length scales. Theoretically, to use this approach in these cases a fine grid is required in the entire domain. In contrast, in approach (2), f k can vary in space allowing a coarser turbulent resolution in regions where large turbulent scales are dominant. This approach can be further subdivided into Static, where f k is fixed in time, or Dynamic for which f k can also vary in time. Between these approaches strong disagreements exist: advocates of the Constant approach claim that by using a varying f k one of the key advantages of the PANS model, the separation of modelling and discretisation error, is destroyed, and the model is reduced to a Hybrid model. On the other hand, advocates of the varying f k approach state that this way the grid, and therefore resources, can be used more optimally. It is argued that applying a constant f k is not reasonable due to the spatial and temporal variation in turbulence length scales and grid resolution. Instead, by varying it in the domain and simulation time, the length scales which can be resolved, are resolved. Note however that the spatial variation in f k reintroduces commutation error, since the PANS filtering operation does not commute with the spatial gradient (Girimaji and Wallin, 2013). Recent work such as Girimaji and Wallin (2013) and Davidson (2016) attempts to account for this error by adding a term in the k and momentum equations based on the gradient of f k .
In the case of Constant f k , it would be beneficial to have an a priori estimate of which physical parameter can be used for a particular flow on a given grid. In the case of varying f k , the need for a reliable estimate for f k is obvious. In the literature on the subject however, there is no consensus on which estimate to use. This paper aims to give an overview of several methods found in literature and their properties. Note that all these estimates only concern f k and all works assume = f 1.0. For more information on the effect of this parameter see Klapwijk et al. (2019b). The different estimates are compared for two test cases: a circular cylinder at = Re 3900, representative of a turbulent wake flow driven by spatially-developing coherent structures (Pereira et al., 2018a;2018b); and a turbulent channel flow at , = Re 395 representative of an internal wall-bounded flow. The results of both Static and Dynamic estimates are evaluated, and compared with the a posteriori computed modelled-to-total ratio of turbulence kinetic energy, f˜, k obtained from Constant f k computations.
The structure of the paper is as follows: Section 2 describes the PANS turbulence model in detail, the f k estimates are given in Section 3 and Section 4 describes the test cases and numerical setup. The applied estimates are compared with the a posteriori computed f˜k in Sections 5 and 6, followed by the conclusions.

PANS model
The Partially Averaged Navier-Stokes equations are obtained by filtering the continuity and momentum equations, thereby decomposing all instantaneous quantities, Φ, into a resolved, ⟨Φ⟩, and a modelled (unresolved) component, ϕ, according to = + (Girimaji and Abdol-Hamid, 2005). The PANS equations for incompressible, single-phase, Newtonian flow are In these equations U i denotes the velocity components, P the pressure, ν the kinematic viscosity, ρ the density and τ(U i , U j ) the sub-filter stress tensor which is modelled using Boussinesq's hypothesis, with ν t the turbulence viscosity, k the modelled turbulence kinetic energy, δ ij the Kronecker delta and ⟨S ij ⟩ the resolved strain-rate tensor, defined as To close the set of equations a RANS model is used. The PANS model in this work is based on the k SST model (Pereira et al., 2018b;Menter et al., 2003). The transport equations of the SST model are reformulated to include the modelled-to-total ratio of turbulence kinetic energy and dissipation rate, This leads to the modified transport equations For the model constants and auxiliary functions, F 1 and F 2 , see Menter et al. (2003), while for more details on the implementation of the PANS model used here, the reader is referred to Pereira et al. (2015).

f k estimates in literature
The f k estimates found in literature are divided according to category (Static, based on an a priori RANS computation, and Dynamic, computed during a PANS computation). Within this manuscript the original notation is modified to maintain consistency between the different estimates and to properly compare them. Some general definitions are the grid sizes For clarity, a distinction is made between L t , based on total (modelled plus resolved) quantities and l t , based on modelled quantities. In the case of Static approaches = L l , t t while in the case of Dynamic approaches, estimates based on both length scales can be found in literature. Note that while l t can be obtained directly from the PANS transport equations, L t must be obtained by computing the instantaneous velocity fluctuations, making the numerical implementation more difficult. Since = f 1.0, dissipation occurs entirely at the smallest scales ( = E and = ). Fig. 1 shows a summary of the f k estimates, sorted per approach, indicating that there is no clear relationship between estimation method and Reynolds number or number of cells. In the literature, the estimates are applied to a range of test cases, including a turbulent jet, swirl in expansion, channel flow, open cavity, backward facing step, bluff bodies, square and circular cylinders, hill and hump flows. There is often little reasoning as to why a particular estimation is applied to a certain test case. This is surprising since the performance of turbulence models is in general highly case dependent. An exception is the work by Luo (2019), in which results using Dynamic f k are compared to those using a Constant f k , as well as from Detached Eddy Simulation (DES), for a backward facing step. The author claims that the Dynamic results are 'almost comparable to the DES computation', with the Constant f k underperforming in predicting skin friction and Reynolds stress profiles. However the applied grid is rather coarse in the wall-normal direction ( + y 1), meaning that the M. Klapwijk, et al. International Journal of Heat and Fluid Flow 80 (2019) 108484 Constant f k computation is not able to properly resolve the boundary layer, leading to poor results. The DES and Dynamic PANS both apply RANS in the boundary layer, leading to superior results. It is also shown in the paper that the difference in results is smaller on a finer grid, indicating that numerical errors may play a role and therefore makingit difficult to generalise this conclusion. Note that in the literature some estimates are explicitly bounded to the interval [0,1], whereas other papers do not mention this. In this manuscript, such explicit bounds are not included to highlight the differences between estimates; of course in the implementation of these methods such bounds should be included.

Static PANS estimates
These estimates are based on an a priori RANS computation, so all turbulence kinetic energy is modelled, i.e. = K k, Girimaji (2004): Abdol-Hamid and Girimaji estimated f k based on the ratio between the unresolved turbulent lengthscale and characteristic grid size using

Abdol-Hamid and
C h is a model coefficient which must be calibrated; in the original paper a value of 1.0 is used. Girimaji and Abdol-Hamid (2005): Girimaji and Abdol-Hamid use an estimate very similar to that of Abdol-Hamid and Girimaji (2004), but with a different constant, and replacing Δ max with Δ min . The estimate is given as In the limit = f 1.0, Eq. (13) reduces to Eq. (12), therefore this estimate is not addressed any further in the current work. Jeong and Girimaji (2010): Jeong and Girimaji define the estimate as with λ T the Taylor scale of turbulence and Δ the grid size. The precise definition of the grid size is not given; in the current work Δ is taken as Δ avg . Surprisingly this estimate uses the grid size in the denominator, while all other methods use a ratio with the grid size in the numerator. This choice is questionable, since this implies that grid refinement leads to a grid which is less capable of resolving structures, which is counterintuitive. The authors do not actually use this estimate in their work, they use a constant f k in the domain.
For a coarse grid l t ≪ Δ, so f k goes to 1. Eq. (16) does satisfy the definition that f k should be bounded between 0.0 and 1.0, which is not a common property of the estimates addressed.

Dynamic PANS estimates
These estimates are evaluated during a PANS computation. At very time step the used f k is computed based on the instantaneous flow field, i.e. f k is updated per time step and is therefore spatially and temporally varying. Since the estimates are evaluated during a PANS computation part of the turbulence spectrum is being resolved, so k ≤ K, ⟨K⟩ > 0 and L t > l t . Elmiligui et al.(2004): Elmiligui et al. use a variable f k defined as following the damping function from a Hybrid model. The turbulent length scale is defined here as = f k is defined in such a way that it equals 1.0 in the viscous sub layer, where the unresolved characteristic length scale tends to be small. f k is also bounded to 1.0. Basu et al.(2007): Basu et al. use Interestingly the definition of the grid size is slightly different by taking the time step and velocity into account: Fig. 1. Literature overview with a selection of the available f k estimates (Abdol-Hamid and Girimaji, 2004;Girimaji and Abdol-Hamid, 2005;Frendi et al., 2007;Jeong and Girimaji, 2010;Foroutan and Yavuzkurt, 2014;Han et al., 2013;Elmiligui et al., 2004;Basu et al., 2007;Song and Park, 2009;Basara et al., 2011;Luo et al., 2014;Davidson, 2014;Basara et al., 2018;Davidson and Friess, 2019;Luo, 2019). The Reynolds numbers are based on the free stream velocity and characteristic length scale; the year indicates the year of publication. The results are shown as a function of the number of grid cells and the approach.
M. Klapwijk, et al. International Journal of Heat and Fluid Flow 80 (2019) Theoretically, not only grid resolution but also temporal resolution determines how much of the turbulence spectrum can be resolved. Most estimates implicitly rely on the user to ensure a sufficient temporal resolution based on Courant number and/or + t . In contrast Eq. (19) incorporates this explicitly, which increases the robustness. Nevertheless, for the grids and time steps considered in this manuscript = in which η indicates the Kolmogorov length scale, defined as = .
3 1/4 An advantage of the approximation, II, is the absence of the singularity, which is present in formulation I. In the current work, the difference between the formulation (I) and the approximation (II) was investigated. For an external flow the difference was found to be negligible, but the singularity in formulation I leads to additional peaks inside a boundary layer. Consequently in the current work the approximation is applied. It is mentioned that for Δ either the maximum or volumetric average can be used, although the authors do not specify which one is used in their work. The differences were found to be marginal, and therefore in the current work Δ avg is employed for this estimate.
This estimate is also used by Davidson (2014), who recognised that on a coarse grid the f k obtained is too high, leading to dissipation of the turbulent fluctuations. Luo et al.(2014): Luo et al. use the same estimate as Abdol-Hamid and Girimaji (2004), but with the inclusion of a different constant, C PANS , which is taken as 0.1. The authors state that the constant should be further calibrated. The estimate is given by Note that this implies that an estimate formulated for Static PANS is employed in a Dynamic approach, resulting in a different estimate for the unresolved length scale, l t , due to the reduction in k. This explains the need for the inclusion of the constant C PANS .
Luo (2019): Luo uses the same estimate as Luo et al. (2014). It is remarked that the estimate is not rigorous in theory and needs additional validation. To this end, for the constant C PANS three values are employed (0.1, 0.3 and 0.5). The value determines the extent of the near wall RANS region, and it is concluded that 0.3 yields the best results. Since this estimate with Luo et al. (2014), and otherwise just 3 or 5 times higher, it is not addressed further in the current work. Davidson and Friess (2019): Davidson and Friess derive an estimate based on the equivalence criterion between DES and PANS (Friess et al., 2015). The estimate is formulated as The estimate is designed to make the PANS model behave as a DES model. According to the authors, the estimate self-adapts, by forcing f˜k towards f k , without the need for computing f˜k. This feature is designated 'passive control' by the authors. Basara et al.(2018): Finally Basara et al. (2018) employ the estimate of Basara et al. (2011), but with L t defined differently, here designated L t . L t is defined based on k t instead of K, where k t is defined as with k ssv the 'scale-supplying' resolved kinetic energy, which requires an additional transport equation. k t is the total kinetic energy (modelled plus resolved), but obtained solely from the additional transport equation and the k equation in the PANS model. This implies that k t ≈ K. An advantage of this is that no (expensive) averaging operations are needed to obtain the resolved, and thereby total, kinetic energy. However, it must be noted that this only works if this extra equation is solved. In Basara et al. (2018) this equation is formulated in the context of the four equation PANS k f model.

Conclusions based on literature
An overview of the required input quantities and properties of the estimates is given in Table 1. A comparison of the formulation of the Static estimates shows that the formulations by Abdol-Hamid and Girimaji (2004), Girimaji and Abdol-Hamid (2005), Frendi et al. (2007) and Han et al. (2013) are essentially the same estimate. The magnitude can differ due to the application of different constants, but the trend is the same. This also shows that there is no consensus on how to define the grid size, which is also true for LES (Pope, 2000). The grid definition could have a large effect on strongly anisotropic grids. The estimate by Jeong and Girimaji (2010) appears to be incorrect due to the use of the grid size in the denominator, while that by Foroutan and Yavuzkurt (2014) is the only one which by definition keeps f k bounded between 0.0 and 1.0, which is a theoretical advantage. The other estimates are most likely explicitly bounded to a maximum value of 1.0, although this is not always clear in literature.
More variation can be found between the Dynamic estimates. Firstly, it is observed that the estimates by Elmiligui et al. (2004), Basu et al. (2007), Luo et al. (2014), Luo (2019) and Davidson and Friess (2019) are all based on l t , so only on modelled quantities. This is questionable since for low f k , the RANS model has little effect on the solution. The reasoning behind this dependence on l t instead of L t is related to the difficulties in obtaining K for statistically unsteady flows, as recognised by Basara et al. (2018) and Pereira (2018). In the case of a statistically steady flow = + K K k, whereby ⟨K⟩ can be obtained from the difference between the instantaneous and mean velocity; for a statistically unsteady flow the difference between instantaneous and mean velocity leads to an overprediction of K due to the energy contained in the large scale motions. An overprediction in K results in reduced values for f k (Pereira, 2018).
Nevertheless, the estimates of Song and Park (2009), Basara et al. (2011) and Basara et al. (2018) are based on L t , and therefore require computing K. An interesting exception to this is the method of Basara et al. (2018), where the total k t is obtained using an additional transport equation. However this estimate therefore only works in the context of a specific PANS formulation, and is thus not applicable in a general PANS formulation. Therefore it is not applied in the current work.
It is observed that due to their formulation, for all estimates (with the exception of Jeong and Girimaji (2010)) = + f lim 1 y k 0 . Some arguments for this behaviour can be found in the occurrence of the smallest length scales at the wall. This implies however that the PANS model is reduced to a Hybrid model with a behaviour similar to DESlike models. For some methods this is mentioned as a goal, while formulating the estimate (Elmiligui et al., 2004;Davidson and Friess, 2019). This behaviour does not happen however, with a Constant f k approach.
Furthermore, it is obvious that due to the application of the different empirical constants any result can be obtained using the different estimates. The authors are therefore of the opinion that the magnitude of the estimation is less relevant than the trend of the estimation methods. All estimates are proportional to Δ n , with often = n 2/3, so grid refinement only affects the magnitude. Consequently, only results for a single grid are shown in this work.

Test cases and numerical setup
The estimates are applied to two canonical test cases: one representative of a turbulent wake flow with coherent structures, and one of an internal boundary layer flow. The numerical solver used for all simulations in this work is ReFRESCO (Vaz et al., 2009), a multiphase unsteady incompressible viscous flow solver using RANS and Scale-Resolving Simulation models, complemented with cavitation models and volume-fraction transport equations for different phases.
The selected test-case for the turbulent wake flow is the flow around a circular cylinder at = Re 3900. This flow was thoroughly investigated using PANS by Pereira et al. (2018a,b). In the current work the finest grid, and set-up, as employed by Pereira et al. (2018a,b) are used. All terms in the equations are discretised with second-order accurate schemes. The rectangular computational domain measures 22D in transverse and 3D in span-wise direction, with an inflow located 10D upstream of the cylinder and the outflow 40D downstream, as shown in In order to investigate the effect of f k estimates inside a boundary layer, a second test case is used: a turbulent channel flow at = = Re u / 395. The setup as employed by Klapwijk et al. (2019a) is used. Computations are made using a rectangular domain, with two noslip walls oriented normal to the y-axis, as shown in Fig. 3. The remaining boundaries are connected using periodic boundary conditions in order to approximate an infinite channel. The Cartesian grid density M. Klapwijk, et al. International Journal of Heat and Fluid Flow 80 (2019)  Time integration is performed using a second-order implicit scheme, the convection terms are discretised using a second-order accurate central differencing scheme. The turbulence equations are discretised using a first-order upwind scheme. Iterative and discretization errors were shown to be negligible in Klapwijk et al. (2019a). The f k values employed in the current work are 0.15, 0.10 and 0.05, following the results obtained in Klapwijk et al. (2019a).
Based on the results previously obtained for these test cases by Pereira et al. (2018b) and Klapwijk et al. (2019a) the employed grids are judged to have sufficient resolution to support the applied f k values.
In order to validate the different f k estimates, the outcomes are compared with the a posteriori computed ratio of modelled-to-total turbulence kinetic energy from a Constant f k computation, designated For a channel flow the computation of u i is straightforward. For a cylinder however, due to the statistically unsteady flow, it is difficult to distinguish between the time-varying mean velocity and the ensemble averaged turbulent velocity. In the results presented here this difference is neglected, leading to an overpredicted value for u i and consequently a reduced f˜k.

Evaluation of Static estimates
The Static estimates are applied to results obtained with the k SST RANS model (Menter et al., 2003). All plots of the estimates are limited between 0.0 and 1.0, even if the estimate itself is not necessarily bounded between these limits. For the cylinder case the results are shown in a contour plot at the centre of the domain in Fig. 4 and the estimates, averaged in spanwise direction, are quantitatively compared on an axial line located on the domain centreline in the vertical direction in Fig. 5. Fig. 4 also shows the time-averaged axial velocity. For the channel flow case, due to the statistical stationarity, only a quantitative comparison is given. Technically, the estimates should be applied to a steady-state computation, however for the cylinder case the flow is inherently unsteady. Therefore the time-averaged quantities are used, immediatelt highlighting a limitation of using Static estimates.
The estimates of Abdol-Hamid and Girimaji (2004), Girimaji and Abdol-Hamid (2005) and Han et al. (2013) vary in magnitude due to the different constants and/or grid sizes but overall show a similar behaviour (see Figs. 4a, b and e). f k is 1.0 (or larger) upstream and near the wall and decreases towards 0.0 in the wake. The lowest values can be found for Foroutan and Yavuzkurt (2014), the largest for Han et al. (2013). It is clear that the estimate of Jeong and Girimaji (2010) is incorrect, f k is 1.0 in the entire domain, except in the first layer of cells near the wake (not visible in the figure). Finally the estimate of Abdol-Hamid and Girimaji (2004) is similar to the estimate of Foroutan and Yavuzkurt (2014), although the values in the wake are somewhat higher and the RANS region near the wall is thicker. These estimates show both a wider wake region where f k < 1.0, and maintain these low values further downstream, compared to the other estimates. The plots in Fig. 5 show that, with the exception of the estimate of Foroutan and Yavuzkurt (2014)

Evaluation of Dynamic estimates
The Dynamic estimates are applied to instantaneous flow fields from PANS computations performed with f k fixed in time and space, denoted as f k,c . This is not how a true Dynamic approach should work, since this way the flow field does not depend on the estimate. The advantage of this approach is that oscillations in the estimates are suppressed. Consequently, the different estimates can be compared more objectively. In the contour plots, the results are again bounded between 0.0 and 1.0, even if the estimate itself is not. The a posteriori computed value f˜k is also shown for comparison. Fig. 6 shows the values of the estimates applied to the cylinder in a contour plot at the centre of the domain, and Fig. 7 shows f k on axial lines located on the domain centreline in the vertical direction, averaged in spanwise direction.
All estimates show an increase in estimated f k with decreasing f k,c , indicating that in a Dynamic approach f k should converge to a target value. If f k,c is larger than this target value, f k is smaller than the target value, and vice versa. As observed by Davidson and Friess (2019), this implies that the estimated f k is implicitly linked to f˜k. Note that due to the spatial and temporal variation of the flow field, the target f k will also vary, leading to potentially oscillatory behaviour for f k .
There is little difference between all the estimates whether  Basara et al. (2011) shows an interesting trend; because of the dependence on L t the wake shows high f k values in the wake centre, but lower values surrounding the wake centre. The difference is clear when comparing the estimate to the one of Luo et al. (2014), which has an almost identical formulation but depends on l t . The estimate of Luo et al. (2014) is unaffected by f k,c in the wake, but increases in the far-field and upstream. Note that the low values in the entire domain for this estimate are mostly related to the small constant (C PANS ) used in the formulation. Finally the estimate of Davidson and Friess (2019) shows a comparable trend, but there is less ambiguity in f k . It is either 0.0 in the wake, or 1.0 elsewhere. The formulation therefore ensures a DES-like behaviour, as was desired in formulating the estimate.
The behaviour observed for the a posteriori computed ratio, f˜, k differs from the estimates. Firstly, the effect of f k,c is clearly visible; as expected with decreasing f k,c , f˜k decreases. Secondly, it can be seen that in general f˜k is significantly lower than f k,c . It appears that modifying f k,c has little effect on f˜k in the entire domain. Instead, it mainly affects the peak values of f˜k occurring in the domain. Thirdly, due to the laminar flow upstream and in the far-field, both k and ⟨K⟩ ≈ 0, leading to f˜0, k which is in strong contrast to the results of the estimated f k values. Finally, f˜k is also low in the near-wall regions, as opposed to the estimates which all give f k ≥ 1 due to The peaks in the wake seem to be best predicted by the estimate of Davidson and Friess (2019), most likely due to the dependence on L t . However outside of the wake the estimate deviates from f˜k. Fig. 7 shows that only the estimates of Elmiligui et al. (2004), Basu et al. (2007) and Davidson and Friess (2019)  the estimates are all larger. It is important to note that only the estimates of Song andPark (2009) andLuo et al. (2014) yield values significantly smaller than 1.0 upstream. This is relevant for cases when synthetic turbulence is added at the inflow, since the introduced fluctuations should not be dissipated before they reach the object of interest. However close to the wall the estimate is still significantly larger than 1.0.

Discussion and conclusions
The review of modelled-to-total kinetic energy f k estimates presented in this paper makes clear that there is no consensus on how to estimate f k from a given flow field on a given grid, both for Static and Dynamic PANS. These approaches are both strongly dependent on this estimate, potentially leading to significant modelling errors. Due to differences in the definition of the characteristic grid dimension and the   M. Klapwijk, et al. International Journal of Heat and Fluid Flow 80 (2019) 108484 application of empirical constants, it is clear that the absolute values of the estimates should be treated with care. Instead more emphasis should be placed on the predicted trends. An issue unaddressed in literature is that Static estimates should be applied to a steady computation; however for statistically unsteady flows, such solutions are unobtainable. In this work, the mean flow field was used. Both the Static and Dynamic estimates do not yield reasonable results for the channel flow case and significantly overpredict f k . For the  M. Klapwijk, et al. International Journal of Heat and Fluid Flow 80 (2019) 108484 cylinder case with a Static computation, the estimate of Foroutan and Yavuzkurt (2014) seems most appropriate, since it is the only one which is properly bounded between 0.0 and 1.0. In case of Dynamic PANS, only the estimates of Elmiligui et al. (2004), Basu et al. (2007) and Davidson and Friess (2019) are bounded between these limits, although that of Basu et al. (2007) generally predicts too high values. It is observed that estimates based on K instead of k generally lead to better predictions, however K is difficult to obtain in statistically unsteady flows. It is shown that the f k value employed in a Constant f k computation, f k,c , mostly affects the peak values of f˜k in the field, and generally < f f k k c , . This difference is sufficiently large that the authors are of the opinion that even if f˜k would be corrected for the energy contained in the large scale motions, still f f , k k c , which is a favourable property of the PANS model. Generally the estimates tend to give values of f k which are significantly larger than f˜k. Aside from the difference in magnitude, the trends observed for the estimates differ in key aspects from the computed f˜, k indicating more fundamental issues. Most estimates are constructed such that . A comparison with f˜k shows that although this principle is correct, the region in which it is applied is not. In the Constant f k computations = f 1.0 k only in the first layer of cells near the wall, whereas in the estimates this occurs in the entire boundary layer. This behaviour also explains the failure of the estimates for the channel flow case, and it gives rise to the belief that the estimates should not be applied inside boundary layers. A consequence of this behaviour is that the PANS model behaves more like a DES model. This is sometimes described in literature as an advantage or a goal in the derivation of the estimate, although this does imply that the unfavourable properties of DES, such as error entanglement, are then also incorporated. A second issue with the estimates is that they all yield = f 1.0 k if the flow is laminar (upstream and in the far-field). This implies that in case of laminar flow, the PANS model resorts to the RANS parent model. For Static computations this becomes problematic if during the subsequent PANS computation synthetic turbulence is added at the inflow, since the introduced fluctuations might be dissipated before they reach the object of interest. In the opinion of the authors, it is not possible to design a general estimate (applicable in the entire domain) which does not suffer from this problem. Upstream of the object no information is available concerning the resolution which can be supported, except for the grid size. The estimates found in literature which depend on k are strongly dependent on values set at the inflow boundary condition, and the decay of modelled turbulence. The estimates depending on K suffer from the fact that it is not possible to estimate f k , unless synthetic turbulence is added. However to resolve the synthetic turbulence, f k should be below 1.0. It seems that the estimates are only valid for cases which show strongly separated vortical structures; and even then only in the wake of the object. To enable the usage for other cases, it is beneficial to limit f k in laminar regions to a certain threshold and only apply the estimate in the wake of the object. For this threshold, currently no definition is available.
Finally it must be remarked that Dynamic PANS computations run the risk that f k will show an oscillatory behaviour due to the strong spatial variation of the estimates. Not only is f k temporally and spatially varying, but also the flow field upon which it is based. This combination might negatively influence the results. Although this hypothesis is not investigated in this work, it contributes to the opinion of the authors that despite potential theoretical advantages to the usage of Dynamic PANS, the Constant PANS approach, with a f k fixed in time and space, is still preferable in order to minimise errors in CFD results and increase the reliability of industrial CFD.

Declaration of Competing Interest
The authors declare that they do not have any financial or nonfinancial conflict of interests.