On the use of synthetic inflow turbulence for scale-resolving simulations of wetted and cavitating flows

Abstract The Delft Twist 11 Hydrofoil is a common test case for investigating the interaction between turbulence and cavitation modelling in computational fluid dynamics. Despite repeated investigations, results reported for the lift and drag coefficient are accompanied by significant uncertainties, both in experimental and numerical studies. When using scale-resolving approaches, it is known that turbulent fluctuations must be inserted into the domain in order to prevent the flow from remaining laminar around the body of interest, although this has been overlooked until now for the present test case. This work investigates the errors occurring when a laminar inflow is applied for mildly separated or attached flows, by employing the partially averaged Navier–Stokes equations with varying values for the ratio of modelled-to-total turbulence kinetic energy, and with varying grid densities. It is shown that depending on the grid resolution laminar leading edge separation can occur. When turbulent fluctuations are added to the inflow, the leading edge separation is suppressed completely, and the turbulent separation zone near the trailing edge reduces in size. The inflow turbulence has a large effect on the skin friction, which increases with increasing turbulence intensity to a limit determined by the grid resolution. In cavitating conditions the integral quantities are dominated by the shedding sheet cavity. The turbulence intensity has little effect on the pressure distribution, leading to a largely unaffected sheet cavitation, although the shedding behaviour is affected. It is shown that, especially in wetted flow conditions, with scale-resolving methods inflow turbulence is necessary to match the experimental flow field.


Introduction
In both research and industrial contexts, computational fluid dynamics (CFD) is being increasingly used to resolve (a part of) the turbulence spectrum, through the application of Scale-Resolving Simulations (SRS), instead of more traditional Reynolds averaged Navier-Stokes (RANS) modelling methods. SRS consist of several approaches, such as large eddy simulation (LES), partially averaged Navier-Stokes (PANS) or RANS-LES Hybrid methods, such as detached eddy simulation (DES). The increased physical resolution in the obtained solutions is necessary for flow cases exhibiting important time-varying stochastic features, such as strongly separated flows, cavitation, broadband noise, etc. It is well known from literature that for attached flows, SRS require synthetic inflow turbulence or leading edge roughness to induce transition to turbulence, in order to prevent laminar solutions (Tabor and Baba-Ahmadi, 2010;Klapwijk et al., 2020). Nevertheless, it is clear from literature that for maritime applications including cavitation, the use * Corresponding author at: Delft University of Technology, Mekelweg 2, 2628 CD, Delft, The Netherlands. E-mail address: m.d.klapwijk@tudelft.nl (M. Klapwijk). of a resolved turbulent inflow is often neglected, potentially leading to large modelling errors.
The current study aims to explore the reduction in modelling errors when using resolved inflow turbulence, and the application of a synthetic turbulence generator for simulating cavitation dynamics. Such a methodology can result in an increased physical resolution while using smaller domains, i.e. less computational cells. Potential use cases include: simulating a propeller in a cavitation tunnel, simulating a propeller in behind condition without resolving the entire flow around the hull upstream, or predicting the interaction between two wings, while only simulating the downstream wing. The test case of choice in the current work, is the 3D twisted hydrofoil studied by Foeth et al. (2006). This is a well documented test case, exhibiting a shedding sheet cavity representative of a ship propeller, while avoiding the additional complications due to a rotating body. The test case was experimentally studied in both wetted and cavitating conditions, and is a common https://doi.org/10.1016/j.oceaneng.2021.108860 Received 13 October 2020; Received in revised form 29 January 2021; Accepted 7 March 2021 numerical test case (Hoekstra et al., 2011). Due to the cavitation occurring inside the boundary layer, the cavitation behaviour is strongly dependent on the interaction between turbulence and cavitation modelling, emphasizing the need for proper inflow boundary conditions. An overview of previous investigations into this test case is given in Section 3.2.
The test case was designed to study cavitation behaviour, which naturally has led to a focus on attempting to capture the cavitating behaviour. Unfortunately, an effect of this is that some of the difficulties in simulating this test case have been overlooked. Although it is known from literature that inflow turbulence is necessary for SRS methods, all SRS results in the open literature for this test case do not employ such methods. Indeed, they obtain reasonable results for the cavity length and shedding behaviour. It seems that the presence of a sheet cavity on the wing surface obscures some of the difficulties in simulating this flow. The laminar inflow leads to an incorrect boundary layer prediction, but still the presence of cavitation leads to a production of turbulence and vorticity due to the occurrence of shear layers in the flow. Due to these turbulent perturbations, turbulence-like structures arise when solving the Navier-Stokes equations, masking the modelling errors due to the laminar inflow. Due to the complex interaction between modelling and numerical errors such results strongly depend on the used grid, time step and turbulence model. Commonly, coarser grids results in lower eddy-viscosity levels (see e.g. Diskin et al., 2015), leading to the occurrence of increased dynamics. In such cases grid refinement can suppress dynamics, thereby leading to an increasing comparison error with grid refinement. In contrast, under wetted flow conditions, the lack of inflow turbulence structures can lead to unrealistic flow results when using SRS approaches.
In the current work we attempt to show the mismatch between computations and experiments which can occur for this case when using a steady inflow for SRS. The focus will therefore be on wetted flow conditions, to prevent the presence of turbulence due to cavitation obscuring the errors. Nevertheless, results for cavitating conditions are also presented. All results are obtained using PANS, in order to be able to systematically vary the turbulence resolution between RANS and (underresolved) direct numerical simulation (DNS), while simultaneously utilising several grid densities. A theoretical advantage of PANS is that due to the explicit setting of the filter between RANS and DNS, the discretisation and modelling errors are decoupled (Pereira et al., 2018;Klapwijk et al., 2020). The use of a single formulation ranging from RANS to DNS prevents ad hoc behaviour when switching between resolving and modelling turbulence, as can occur for Hybrid methods Klapwijk et al., 2020). Iterative, statistical and discretisation errors are assessed. The turbulent unsteady inflow is provided by a modified version of the body-forcing approach developed in Klapwijk et al. (2020), which is based on the digital filtering method by Xie and Castro (2008).
The paper is organised as follows. Section 2 describes the different mathematical approaches for simulating turbulence, phase change due to cavitation and synthetic inflow turbulence; Section 3 describes the test case and numerical setup, while Section 4 details the estimation of the numerical errors. After this, the results are addressed in two sections: wetted flow in Section 5, and cavitating flow in Section 6. The Discussion and Conclusions in Sections 7 and 8 complete the study.

Mathematical modelling approaches
The approach to resolving and modelling turbulence employed in this work is the partially averaged Navier-Stokes (PANS) method (Girimaji and Abdol-Hamid, 2005). In this methodology, the instantaneous quantities, , are decomposed into a resolved, ⟨ ⟩, and a modelled (unresolved) component, . Applying this decomposition to the conservation of mass and Navier-Stokes equations for a Newtonian fluid, including phase change by cavitation, yields, In these equations denotes the velocity components, the static pressure, the kinematic viscosity and the density. The subscripts and indicate liquid and vapour phases respectively, symbols without subscript refer to the mixture quantities, defined according to anḋindicate the vapour volume fraction and source term provided by the cavitation model (see Section 2.2), and ( , ) denotes the modelled Reynolds stress tensor, which is computed using Boussinesq's hypothesis, with the eddy-viscosity, the modelled turbulence kinetic energy, the Kronecker delta and 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 − SST model (Pereira et al., 2018;Menter et al., 2003). The transport equations of the SST model are reformulated to include the ratio of modelled-to-total turbulence kinetic energy and dissipation rate, This leads to the modified transport equations and For the model constants and auxiliary functions, 1 and 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. (2005). The filtering of the Navier-Stokes equations depends on the values chosen for and . determines the physical resolution of the flow, i.e. to what extent the turbulence spectrum is resolved. determines the overlap between the energy-containing and the dissipation ranges. Following Klapwijk et al. (2019b), = 1.0 to avoid adding excessive dissipation.
can either vary in space and time, or be kept constant in the domain. The drawback of the first approach is that this entangles the modelling error and numerical error, thereby destroying one of the key advantages of the PANS model. As shown in Klapwijk et al. (2019a), there is no consensus on what estimate for to use, and the estimates are unreliable in boundary layers and problematic to apply in combination with inflow turbulence. Consequently, in this work constant, a priori chosen, values of (1.00, 0.75, 0.50, 0.25 and 0.00) are employed. By definition, PANS computations using = 0.00 are equal to implicit LES computations, also sometimes referred to as underresolved DNS, where it is assumed that the added numerical diffusion due to the use of coarse(r) grids and low order (upwinding) schemes acts as a sub-filter model. Previous studies in Klapwijk et al. (2019a) showed that the real ratio of modelled-to-total levels of turbulence solved is lower than the a priori chosen values, which is a desirable conservative characteristic of the approach.

Synthetic turbulence generator
Synthetic inflow turbulence is generated using the method developed in Klapwijk et al. (2020), which is a modified version of the digital filter method by Xie and Castro (2008) and Kim et al. (2013). In this method, random numbers, , , , with a zero mean and unit variance, are generated on a Cartesian grid at each time step. Here , indicate the position indices and the velocity component. These numbers are spatially correlated using after which they are temporally correlated with the numbers generated during the previous time step using Here  = ∕ is the Lagrangian time scale based on the desired integral length scale . and are filter coefficients used to generate spatial correlations. The spatially and temporally correlated numbers are transformed to velocity fluctuations using in which indicates the Lund transformation matrix, which is based on a Cholesky decomposition of the desired Reynolds stress tensor . Following Klapwijk et al. (2020), the velocity fluctuations are transformed to body-forces term in the momentum equations, explicitly added on the right hand side of the equations. Note that in the current work a string of pseudo-random numbers, i.e. the same range of random numbers for every computation, is used, to compare different computations.
To improve the iterative convergence, a divergence-free velocity field is enforced by adding a source term, , to the mass equation every non-linear loop. This source term is computed by computing the flux, , through each face for every cell of the turbulence generator, The source term equals a summation of the fluxes over the faces, and is added to maintain mass conservation in the cells of the turbulence generator.

Cavitation model
The multiphase flow is modelled using the homogeneous mixture Eulerian Volume of Fluid approach. An additional transport equation is solved for the vapour volume fraction, =  ∕( +  ), with  indicating the phase volume. The transport equation is formulated as From the mixture properties can be calculated using Eq. (3), under the condition that The source terṁis modelled by the Schnerr-Sauer cavitation model (Schnerr and Sauer, 2001;Hoekstra and Vaz, 2009), which is based on the linearised Rayleigh-Plesset equation for bubble dynamics and reads, The dependence on the local pressure leads to an vastly simplified inception criterion compared to reality, but the method is widely applied in literature. indicates the maximum bubble radius, the bubble concentration, and the vapour pressure. In the present study, the bubble radius and concentration were set to = 1 × 10 8 and = 1 × 10 −5 , respectively, based on Vaz et al. (2017).

Test case description
The Delft Twist 11 Hydrofoil is a NACA0009 profile with a spanwise-varying angle of attack from 0 • at the sides to 11 • at midspan, mounted with an angle of attack at the wall = −2 • . The chord length = 0.15 m and the span = 2 . The spanwise twist is given by (Foeth et al., 2006;Foeth, 2008 The Reynolds number is = ∞ ∕ = 1 × 10 6 , leading to an uniform inflow velocity of ∞ = 6.97 m∕s. In cavitating conditions the cavitation number is = In the experiments sand roughness with a grain size of 10 −4 m was applied at the leading edge at ∕ = 0.04 to force transition to turbulence, but also leading to an increase in drag. Foeth (2008) measured the lift force, the pressure at different locations on the suction side and applied particle image velocimetry (PIV) to analyse the developing cavity shape and shedding behaviour. The drag was not measured. For the lift and pressures, the calibration errors of the sensors are reported, but no total uncertainties are given. Regarding the setup, an uncertainty of 2.7% is reported for the inflow velocity, and 5% for the cavitation number. Inflow turbulence levels of the cavitation tunnel were 2%-3% at the location of the hydrofoil. This is higher than the more recent values reported for the same cavitation tunnel by Varadharajan (2019) (≈ 1.5%), but those were obtained with a 50% lower mean tunnel velocity. From both sources no information is available on the integral length scale, which makes matching the experimental setup and quantitative validation difficult. For more details the reader is referred to Foeth (2008).

Literature overview
The Delft Twist 11 Hydrofoil is a common test case for investigating the interaction between turbulence and cavitation modelling in CFD. In literature, a number of different turbulence modelling approaches for this case have been applied, such as Euler equations (Koop, 2008;Schnerr et al., 2008), unsteady RANS (Oprea, 2013;Bensow, 2011;Ji et al., 2014a;Whitworth, 2011;Vaz et al., 2017), unsteady RANS with an ad-hoc eddy-viscosity correction (such as the so-called Reboud correction Reboud et al., 1998) (Bensow, 2011;Li et al., 2010;Vaz et al., 2017;Hong and Zhang, 2020), Reynolds Stress Models (RSM) (Whitworth, 2011), LES (Bensow, 2011;Asnaghi et al., 2017;Chen et al., 2017;Ji et al., 2013;Long et al., 2018;Lu et al., 2010), Hybrid models (Bensow, 2011;Whitworth, 2011;Vaz et al., 2017) and PANS (Ji et al., 2014b). The use of unsteady RANS typically suppresses the cavity dynamics, while LES results show a shape and shedding behaviour which is more in agreement with the experiments, although the cavity length is underpredicted. Hybrid models, such as DES, shed less cavity structures due to the sheet cavity being close to the wing surface, i.e. in a RANS region. The use of the Reboud correction leads to increased cavity dynamics and cavity length, however its ad hoc empirical nature introduces additional modelling error in the results. The PANS results reported show a shedding behaviour comparable to experiments, however only one ratio of modelled-to-total turbulent kinetic energy was investigated. From literature it is know that the results can vary significantly based on this ratio (Pereira et al., 2018;Klapwijk et al., 2019aKlapwijk et al., , 2020. Most investigations (Oprea, 2013;Li et al., 2010;Bensow, 2011;Whitworth, 2011;Ji et al., 2013Ji et al., , 2014aLidtke et al., 2016;Vaz et al., 2017) make use of the Schnerr-Sauer cavitation model, which is also employed in the current work.
An overview of the reported mean lift coefficient = ∕( ∞ ) and Strouhal number = ∕ ∞ as function of the number of grid cells in cavitating conditions can be found in Fig. 1. The coefficients are defined using the lift force and the shedding frequency , and divided by the turbulence model. The mean lift coefficient is underpredicted by every simulation found in the open literature, which indicates that the cavitation extent is underestimated (Lidtke et al., 2016). A wide variety of grid types and densities is employed in the literature. It must be noted that, in general, the grids are relatively coarse considering the Reynolds number, with high wall-normal resolutions. The stream-and spanwise resolutions are regularly not reported, and the sensitivity of the results towards the grid is often not investigated. This leads to few reported uncertainties for integral values.
As mentioned in the introduction, the literature focusses on cavitation conditions. Nevertheless, some results for wetted flow are reported (Hoekstra et al., 2011;Lidtke et al., 2016;Vaz et al., 2017;Asnaghi et al., 2018). These results are mostly limited to the forces and the pressure distribution on the centreline, no flow visualizations are given. Hoekstra et al. (2011) do report that within their workshop 'all participants report an attached boundary layer on the entire foil', i.e. no boundary layer separation is observed on the wing surface. See Fig. 2 for an overview of the reported mean lift and drag coefficient = ∕( ∞ ) as function of the number of grid cells and turbulence model. Some results match the experimental lift coefficient well, however again the spread in results is significant.

Numerical setup
The computational domain extends two chord lengths upstream of the leading edge and five chord lengths downstream. Half of the wing is modelled with a symmetry plane at the centre of the wing. The use of a symmetry plane in the setup can lead to modelling errors when attempting to resolve turbulence using SRS and when adding synthetic turbulence. Nevertheless, it was employed to reduce numerical cost. It must be remarked that most results in literature, including the LES results, employ the same simplification. The top and bottom of the domain are located one chord length from the wing, matching the dimensions of the cavitation tunnel used by Foeth et al. (2006). The boundaries of the domain are modelled as slip-walls, whereas on the wing a no-slip boundary condition is applied. At the inlet an inflow boundary condition is used, fixed pressure is defined at the outlet. Fig. 3 shows the wing geometry, computational domain and grid topology. The boundary conditions at the inlet are a Dirichlet condition for all velocity components and turbulence quantities, and at the outlet a Dirichlet condition for the pressure. The wing's surface is modelled as a non-slip wall, with a symmetry boundary condition at the centre of the wing. All other domain boundaries are modelled as slip walls. Four geometrically similar, multi-block hexahedral structured grids ( 1, 2, 3 and 4) are used, but for comparison purposes also an additional grid ( 0 * ) is employed. This grid is identical to grid 1, but it includes a local refinement box with dimensions 4∕3 × × 2∕3 surrounding the wing. In this refinement box the grid is refined by a factor 2 with respect to grid 1. The average values of the non-dimensional wall cellsizes + , + and + (normal, chordwise and spanwise) on the suction side centreline of the different grids are given in Table 1. Due to the 3D geometry and grid topology, the averaged values do not give a good indication of the cell distribution over the wing. Fig. 4 shows the + , + and + distribution over the chord at the wing centreline on the suction side. Over the entire wing + is well below 1 for all grids. On the suction side, + reaches high values at the leading edge, drops immediately after the leading edge and then increases again. It remains relatively constant but decreases again towards the trailing edge. Finally, + is high at the leading edge, and decreases along the wing. The values for + at the wing midspan are lower than the values towards the wall of the cavitation tunnel (where = −2 • ).
When comparing these resolutions to guidelines from literature for a well-resolved LES or PANS, + < 1, 50 < + < 150 and 15 < + < 40, Piomelli and Balaras (2002), it is clear that in wall-normal direction, all grids are sufficiently fine. Grids 0 * and 1 do comply with the required + , however none of the applied grids reaches the required resolution in spanwise direction. Due to the computational cost, the grids are not refined below these levels.
For computations with a resolved turbulent inflow, the turbulence generator is located at = −2 with respect to the leading edge of the wing, which is close to the inflow boundary condition.

Solver
The numerical solver used for all simulations in this work is Re-FRESCO ), a SIMPLE based, multiphase unsteady incompressible viscous flow solver using RANS and SRS methods, complemented with cavitation models and volume-fraction transport equations for different phases. The code traditionally focusses strongly on code and solution Verification and Validation, as one can see by several studies Pereira et al., 2018;Liebrand et al., 2020;Vaz et al., 2015;Lloyd et al., 2017). For the current work, time integration is performed using a second-order implicit three time level scheme, the convection terms in the momentum equations are discretised using a second-order accurate flux limited Quadratic Upstream Interpolation for Convective Kinematics (QUICK) scheme, the turbulence and cavitation equations are discretised using a firstorder upwind scheme. Diffusion is central second-order accurate, and non-orthogonality and eccentricity of the grid are considered by extra deferred corrections.

Numerical errors
As is generally accepted, numerical errors can be divided into input, round-off, iterative, discretisation, and, in the case of unsteady computations, statistical errors (Roache, 2009). The round-off error is assumed to be negligible due to the use of double precision arithmetic (Eça and Hoekstra, 2006). One of the sources of input error are boundary conditions. The effect of changing the inflow boundary condition is assessed in Section 5. In this section, the iterative, discretisation and statistical error are assessed.
The iterative convergence is assessed based on the residuals, normalised by the diagonal element of the left-hand-side matrix of the  (Oprea, 2013;Li et al., 2010;Bensow, 2011;Whitworth, 2011;Ji et al., 2013Ji et al., , 2014aLidtke et al., 2016;Long et al., 2018;Vaz et al., 2017;Hong and Zhang, 2020;Asnaghi et al., 2018).   linear system of equations. All wetted flow calculations with steady inflow condition are converged until ∞ is below 10 −6 . Computations with a resolved turbulent inflow are converged until 2 is below 10 −4 and ∞ below 10 −2 . The largest residuals occur near the turbulence generator and the wing leading edge. Residuals can be reduced further by increasing the number of outerloops, with the exception of the equation which stagnates. For cavitating computations, the 2 -norm is generally in the order of 10 −3 , while the ∞ -norm is in the order of 10 −1 for all equations except for the equation. For this equation the 2 -norm is in the order of 10 −3 and the ∞ -norm in the order of 10 0 -10 −1 .
The finite length of a CFD simulation introduces a random uncertainty in the mean of the signal. To estimate the statistical uncertainty, and to remove the start-up effect, we use the Transient Scanning Technique (TST), developed by Brouwer et al. (2015), and applied before in Klapwijk et al. (2020). For wetted flow cases with a steady inflow or low levels of inflow turbulence intensity ( < 10%), the statistical uncertainty for integral quantities is below 1%, for higher inflow turbulence intensities ( > 10%) it is below 5%. For cavitating computations, the statistical uncertainty for the mean lift and drag forces is below 3%.
The discretisation error, ( ), is estimated using a power series expansion (Roache, 2009;Pereira et al., 2018) where 0 indicates the estimated solution for zero discretisation error and is the refinement ratio, based on with ℎ the cell length. Based on the data it is not possible to obtain an accurate estimation of , so both a value of 1.0 and 2.0 are used. = 1.0 leads to a conservative estimate, , whereas = 2.0 yields a lower value, , .

Wetted flow results
This section describes the results for wetted flow simulations, i.e. without cavitation. Among the results reported are the lift coefficient = ∕( ∞ ), the drag coefficient = ∕( ∞ ), the pressure coefficient = ( − ∞ )∕ ∞ and the streamwise skin friction coefficient , = , ∕ ∞ at the wing midspan. The results are compared against experimental results by Aït Bouziad (2005) (designated EPFL), and against results obtained with the panel code XFoil (Drela, 1989). With XFoil the 2D NACA0009 cross-section was computed at an angle of attack at the centreline ( = 9 • ). Since the 3D effects due to the twist are not included in XFoil, differences in the magnitude of are to be expected. Nevertheless, it gives a good indication of the pressure distribution along the chord.

Steady inflow condition
Most of the computations with a steady inflow converge to a steady solution, with the exception of flow solutions which show leading edge separation extending along the chord. Consequently, the results in this section are obtained with a larger time step which is kept fixed to * = ∞ × ∕ = 6.97 × 10 −3 , leading to a maximum and average Courant number of 48 and 2.3 on the finest grid, and 22 and 1.2 on the coarsest grid, respectively. These large Courant numbers were deemed acceptable based on preliminary computations with a smaller timestep, which showed the same solutions with a similar iterative convergence. Fig. 5 shows the lift and drag coefficient as a function of the grid refinement ℎ ∕ℎ 1 and the physical resolution , together with the discretisation error as error bars. The plotted values are also given in the Appendix. Both the lift and drag have a low grid dependency, leading to small discretisation errors.
has little influence on the lift coefficient in the range 1.00−0.50, however the lift coefficient suddenly drops by almost 10% for = 0.25. For = 0.00, lift increases again. We will comment on this later. The drag coefficient shows a linear decrease with decreasing , with the exception for = 0.25 on the two finest grids, and = 0.00, where the drag suddenly increases. A division of the drag into pressure and friction drag shows that the friction drags linearly decreases with . However for = 0.00 and 0.25 on the fine grids the pressure drag doubles, due to flow separation occurring at the leading edge, as will be explained later based on flow visualizations.
The effect of reducing is visualized in Fig. 6 by showing the timeaveraged streamwise skin friction coefficient, limiting streamlines and time-averaged normalized streamwise velocity ∕ ∞ . With a reduction in , the streamwise skin friction decreases. For = 0.25, the flow starts to change. Depending on the grid resolution, two solution regimes are predicted. For the coarser grids ( 3 and 4), turbulent separation is observed in the streamlines near the trailing edge along the entire span of the wing. Along the spanwise position the angle of attack varies, the separation region moves towards the leading edge with increasing angle of attack. This is accompanied by a negative friction coefficient on the wing surface, indicating recirculating flow over the wing. On the finer grids ( 1 and 2), the separation location moves further upstream, and the flow exhibits laminar separation close to the leading edge. The flow solution is now also unsteady, which together with the separation at the leading edge, is a clear indication of a stall condition. After the laminar separation bubble the flow periodically detaches and reattaches, leading to a time-averaged positive skin friction, which increases the pressure drag as observed in Fig. 5. The grid sensitivity of the flow separation is likely related to the + resolution behind the leading edge (see Fig. 4). While it is low at the leading edge, on grids 3 and 4 it more than doubles over a chordwise distance of ≈ ∕30, thereby suppressing the flow separation. On the finer grids, + still increases, but the increase is less drastic, allowing the grid to resolve the flow separation. The results for = 0.00 show this stall condition on all grids; due to the absence of a turbulence model, the flow instability is not suppressed by any eddy-viscosity. Again the timeaveraged skin friction is positive behind the separation bubble at the leading edge, leading to an increased pressure drag. For this value, the largest differences occur more towards the tunnel side wall where the angle of attack is lower. On the coarse grid there is a large region with negative friction, which decreases in size on the finer grids.
To quantify the differences in time-averaged skin friction coefficient, Fig. 7 shows , at the wing centreline. , decreases with decreasing . For = 1.00, , is close to the XFoil prediction, the deviation in slope is likely related to 3D effects which are not included in the XFoil prediction. The trend remains the same in the range 0.5 ≤ ≤ 1.0, only the magnitude varies. For = 0.25 on the coarse grids, , shows a trough for 0.05 < ∕ < 0.1, after which it increases until ∕ = 0.2 and then decreases again until it becomes negative for ∕ > 0.4. On these coarser grids, the low values are likely related to an absence of transition, i.e. the flow remains laminar along the chord. On the finer grids, the flow separates at the leading edge, so , shows a peak at ∕ ≈ 0.05, after which it becomes negative until ∕ ≈ 0.4, and then remains positive until the trailing edge. The leading edge separation induces transition at ∕ ≈ 0.4, leading to a higher than observed for = 0.25 on the coarser grids. Fig. 6 already showed that this pattern varies significantly along the span of the wing. Finally for = 0.00, the pattern is similar on all grids. There is a large separation region at the leading edge, inducing transition. , on the fine grids ( 2 and 1) equals the , predicted by = 0.25, which is in line with the earlier visual observations.
The pressure distribution at the wing surface is less affected by the flow pattern, as shown in Fig. 8. The different computations show little difference, with the exception of the = 0.00 and 0.25 on grids 2 and Values are also given in Table 5 in the Appendix.
1, which is of course directly linked to the unsteady flow separation. The other computations predict similar distributions, but for = 0.25, is slightly lower in the range 0.2 ≲ ∕ ≲ 0.7, and higher for ∕ ≳ 0.7. Two observations can now be made. Firstly, both , and are strongly affected by , and vary for each value; while and show two sets of solutions based on being above or below a certain threshold. Secondly, in the range ∕ ≤ 0.4, does not depend on , which is favourable for cavitation inception predictions. The cavitation behaviour and attached sheet cavity extent will be insensitive to , and consequently to the lack of inflow turbulence.
What happens when we reduce ? The change in the flow is related to changes in the turbulence intensity, both modelled and resolved. Due to the reduction of the effect of modelled turbulence is reduced, while due to the steady inflow boundary condition the resolved flow is laminar (i.e. the resolved turbulence intensity is zero). This is clearly visible in Fig. 9, which shows , and ∕ between the location where we will later insert resolved inflow turbulence and the leading edge. Both and decrease downstream, but since ∕ depends on the ratio of and (Eq. (9)), ∕ is relatively constant along the streamwise direction. However with decreasing , ∕ decreases, leading to the flow becoming laminar with decreasing . Note that for = 0.00, ∕ = 0 by definition. Due to low inflow turbulence level, bypass transition on the wing surface is delayed, and the flow separates at the trailing edge. For even lower , the flow becomes completely laminar and separates at the leading edge. This was already observed by Foeth (2008), who recognised that for the limited Reynolds numbers typical for model scale hydrofoils transition to turbulence does not occur at the leading edge unless the boundary layer is locally disturbed, hence the application of surface roughness in the experiments. It is known from literature that laminar boundary layers separate earlier than turbulent boundary layers (Rist and Augustin, 2006). The occurrence of turbulent trailing edge separation leads to the decrease in lift coefficient, and the decrease in drag coefficient is proportional to the reduction in skin friction due to the laminar boundary layer along the chord. With leading edge separation the lift and drag coefficient increase again, due to a different pressure distribution. Unsteady structures can be observed (e.g. based on the −criterion), potentially leading to a wrong conclusion that the approach is resolving a turbulent boundary layer.
Up to now, with the exception of the work by Foeth (2008) who employed surface roughness, separation has not been reported in literature for this test case. A number of reasons can be identified to explain this. Simulations using RANS or Hybrid models do not show flow separation due to the inherent assumptions employed in RANS in the boundary layer, leading to transition to turbulent flow at the leading edge. For SRS results, such as LES and PANS, it is shown that for intermediate values, such as 0.25, the behaviour is grid sensitive and only occurs on fine grids. The literature overview in Figs. 1 and 2 showed that most LES results reported in literature are obtained using coarser grids than the ones here, thereby potentially hiding this behaviour. It must also be remarked that most investigations in literature focus on cavitation dynamics. However, due to the nature of the test case, when cavitation is included, this separation is again masked since due to the cavitation inception criterion used in CFD, the sheet cavity occurs at the same location as the flow separation at the leading edge. The use of XFoil further confirms the hypothesis that the leading edge separation with a laminar inflow is physical. When using a turbulent inflow, similar , and distributions are found, however when reducing the inflow turbulence (when moving from what is known in XFoil as a 'dirty' towards a 'clean' wind tunnel) the solution becomes more and more difficult to converge, again pointing to unsteady flow separation at the leading edge.
To ensure that further grid refinement does not yield different flow behaviour, additional = 0.25 and 0.00 computations were performed M. Klapwijk et al. Fig. 6. Limiting streamlines and time-averaged skin friction coefficient , on the wing surface, and streamwise velocity ∕ ∞ , for different values, using different grids with steady inflow. Results for 4 are similar to 3 and are therefore not shown. Flow from left to right. on a grid with a local refinement box surrounding the wing (grid 0 * ). The integral values are = 0.35 and 0.42, × 10 2 = 1.00 and 0.85 for = 0.25 and 0.00 respectively. The drag appears to continue to increase with increasing simulation time. Fig. 10 shows , and at the midspan, along with a flow visualization for = 0.00. For this additionally refined grid, the wing does not show such a large region of recirculating flow at the leading edge, as for grid 1, leading to a pressure distribution closer matching the experiments. At a first glance it appears that refining the grid to this level solves the problems observed previously in this section. However, investigation of shows that leading edge separation still occurs. The difference is that here becomes positive again at ∕ ≈ 0.1, after which it remains positive until close towards the trailing edge. On grid 1, also seemed to recover at ∕ ≈ 0.07, but then decreases again and does not become positive until ∕ ≈ 0.4. The increased grid resolution, decreases the leading edge separation, and therefore leads to a reasonable prediction of the drag. The lift however, is even lower than for the PANS results on coarser grids. Potentially, the leading edge separation is reduced even more on a further refined grid, however it is important to realise that this grid is already refined to a level which is currently unaffordable for industrially oriented cases. This grid is significantly finer than LES grids employed in the literature for this test case.
To summarise, it is clear that flow predictions with SRS with a steady inflow do not match experimental conditions, but it could be argued that this would have little effect on cavitation predictions due to the inception criterion used in CFD.
does not vary, so the inception behaviour will be the same. However, we emphasize that the laminarization of the flow is a problem, as also known from literature (Tabor and Baba-Ahmadi, 2010;Klapwijk et al., 2020). Firstly, does change towards the trailing edge, potentially affecting the development and dynamics of shed cavities. Secondly the integral quantities (lift and drag) differ significantly from experiments, making validation impossible. It is interesting to note that the lift decreases with , away from the experimental value, with the exception of = 0.00, which gives M. Klapwijk et al. higher values than = 1.00. From this, the incorrect conclusion could be made that = 0.00, also known as Implicit LES, does not suffer from these problems and yields the best match for experimental values. Nevertheless, as shown, , , and are different, and in fact the entire flow changes, from a steady attached flow to an unsteady flow separating at the leading edge. Indeed, the lift is better predicted, but for the wrong reasons. Instead it is necessary to introduce synthetic inflow turbulence to obtain a turbulent boundary layer from the leading edge onwards, to better match the experimental results.

Resolved turbulent inflow condition
To trigger the flow to become turbulent and suppress unphysical separation, synthetic turbulence is added at the inflow. The need for such boundary conditions is well known for LES and PANS (Tabor and Baba-Ahmadi, 2010;Klapwijk et al., 2020), however LES results in the open literature for this specific test case do not employ such methods.
The turbulence generator is located at = −2 , with turbulent fluctuations inserted in a plane perpendicular to the flow with dimensions × and a thickness in primary flow direction of ∕15, which corresponds to 2 − 3 cells in streamwise direction. A number of conditions with varying levels of turbulence intensity and integral length scale are investigated, the settings at the location of the turbulence generator are summarized in Table 2. The decay in turbulence kinetic energy is compared to the slope of theoretical decay of isotropic homogeneous turbulence (Comte-Bellot and Corrsin, 1966) In this equation indicates the measurement location and is approximately equal to 1.27, based on experiments. Table 2 also shows the turbulence intensity,  = √ 2∕3⟨ ⟩∕ ∞ = √ ∕ ∞ , expected at the leading edge based on the theoretical decay. For Case A and C the turbulence intensity is comparable to the experimental setup; no integral length scale was reported for the experiments.
From a theoretical perspective, it is incorrect to insert velocity fluctuations in a computation with = 1.00, since for such a value all turbulence should be modelled. In this work it is still done, to investigate the trends. Due to the use of pseudo-random numbers the curves for the different computations overlap.
The time step for these computations is kept fixed to * = ∞ × ∕ = 6.97 ⋅ 10 −4 , with maximum Courant numbers of approximately 5 occurring close to the leading edge. In the remainder of the domain the Courant number is well below 1. The 2 norms of the residuals are (10 −4 ), (10 −6 ), (10 −7 ) and (10 −7 ), for momentum, pressure, turbulence kinetic energy and dissipation, respectively. Of these equations only the equation stagnates, the other equations keep on converging and decrease at least one order of magnitude during a time step. The ∞ norms for the same equations are (10 −2 ), (10 −3 ), (10 −4 ) and (10 −3 ). The largest residuals occur in the cells where the turbulence generator is located. During the computation the flow passes the chord seven times, with averaging applied during the final four time units ( = ∕ ∞ ), leading to a maximum statistical error of 5%, which occurs for the case with highest inflow turbulence intensity.
First, cases A and B are compared on all grids with 0.0 ≤ ≤ 1.0. The streamwise development of the synthetic isotropic turbulence between the generator and the leading edge, is shown for = 0.25 in Fig. 11. The Reynolds stresses decrease around two orders of magnitude on the coarse grid, while on the fine grid the decrease is about one order of magnitude, which is comparable to Eq. (21). An initially surprising observation is that the decay of resolved turbulence is almost    insensitive to the value (not shown in the graph), suggesting no effect of the chosen on the computations. However, the statistical representation of turbulence (the 'RANS' contribution) does vary with , a higher leads to a higher eddy-viscosity. As observed earlier, on finer grids the eddy-viscosity levels also increase.
These results, in combination with results from literature (Pereira et al., 2005;Klapwijk et al., 2020) suggest that while the selected value of strongly affects the production of modelled turbulence, it does not affect the decay of resolved turbulence in this streamwise range. While this might seem counterintuitive, it can be explained from the equations which are being solved. The filtered Navier-Stokes equations are independent of , only affects the turbulence transport equations for and . Due to the formulation of the turbulence generator fluctuations are introduced regardless of , which can lead to the unphysical combination of = 1.00, i.e. RANS, with resolved turbulent fluctuations. After their addition to the flow, the development of the fluctuations is described by the Navier-Stokes equations. The difference between computations with different values, is the level of ∕ . However in this short range the effect of the increased eddy-viscosity is too small to significantly dampen the velocity fluctuations, and therefore the decay of resolved turbulence is comparable. A similar observation was made when using synthetic turbulence for a channel flow in Klapwijk et al. (2020), close to the turbulence generator the solutions for different values are comparable, and velocity fluctuations are being dampened only further downstream. The comparable decay is also an indication that the grids are fine enough to not add excessive diffusion.
In Figs. 12 and 13, the flow is visualized using the time-averaged streamwise velocity and friction coefficient at the wing surface, for different values on grids 4 and 1. The visualization shows that while the flow at the surface is affected by the inflow turbulence, for = 0.25 a region of separated flow still occurs near the trailing edge. The size of the separated flow region is reduced compared to the results obtained with a steady inflow, especially on the finer grid. The separation is also smaller for Case B than for Case A. When = 0.00, the region of separation is larger, especially on the coarser grid. This is an indication of the increased grid sensitivity of ILES, due to the absence of a sub-filter model. Note that for both these levels of inflow turbulence, for all values, the large flow separation at the leading edge as observed for low values without inflow turbulence (Fig. 6) disappears.
The effect on the skin friction coefficient is further analysed in Fig. 14, by plotting , versus ∕ at midspan for varying for the two cases. A comparison with the steady inflow results shows that , increases with the addition of inflow turbulence. An increase in also leads to an increase in , , independent of used in the computation. As expected, = 1.00 gives a RANS distribution, where , has a peak at the leading edge and then decreases along the chord, while remaining positive. The magnitude is comparable to the XFoil prediction. With decreasing , the magnitude of , decreases along the chord. For = 0.50, there is a slight kink at ∕ = 0.1, but the profile remains similar to the RANS results. However the result for = 0.25 clearly shows not only a much lower , , but also drops after the leading edge to a negative value, after which it increases again to a constant value along the wing. This shows that the flow is still laminar at the leading edge, but transitions to turbulent flow at ∕ = 0.1, indicating that the inflow turbulence is not sufficient to trigger a turbulent flow at the leading edge. For Case A, it can also be observed that for ∕ > 0.9, , still becomes negative, indicating turbulent flow separation. It is noteworthy that on finer grids , always increases, independent of . For = 0.25, this has the effect that, while still present, the separation region at the trailing edge is reduced (as was observed previously in Fig. 12). On the coarse grid, , for = 0.00 remains close to zero for both cases, indicating that transition does not occur and that the flow remains laminar. On the finer grid it increases, but remains close to zero along the chord.
When comparing Case B to Case A, we can see that the increase in inflow turbulence leads to an increase in , along the chord for all values, except for = 0.00 on the coarse grid. For Case B, on the coarse grid there is still a small region where , < 0 close to the leading edge, while on the fine grid , remains positive along the wing.
The effect on the pressure coefficient is limited, and is therefore not shown. There is no difference between different values, except at the trailing edge ( ∕ > 0.9) in the location of turbulent flow separation. As expected, this difference is larger on the coarse grid, 4, and larger for Case A than for Case B. For Case A, there is a still a difference between = 0.25 and 0.00 and the other values on the finest grid; for Case B there is no longer an observable difference between the different values on the finest grid.
The effect on the integral quantities is presented in Fig. 15, where the time-averaged lift and drag coefficient are given as a function of the grid refinement and . The predicted trend for the lift for Case A is similar to what was observed without inflow turbulence (Fig. 5). There is a sudden decrease in for = 0.25, but the predicted lift by = 0.00 is again approximately 5% larger than for = 1.00, which is caused by a higher pressure on the pressure side of the wing between    the leading edge and midchord, close to the tunnel wall. For Case B the decrease at = 0.25 is smaller, and is more constant across the range, especially on the finest grid. The addition of inflow turbulence leads to a temporally varying lift coefficient, since now the effective angle of attack varies in time. As expected the standard deviation for Case B is larger than for Case A due to the higher turbulence intensity. For Case A, × 10 3 ≈ 6 while for Case B, × 10 3 ≈ 23. For the drag coefficient, in both cases A and B, the trend is the same as for the case with a steady inflow: decreases with decreasing . For both cases, increases across the range compared to the steady inflow, and , for Case B is larger than for Case A. Again the standard deviation in the drag coefficient is larger for Case B than for Case A, i.e.
× 10 3 ≈ 1.3 versus × 10 3 ≈ 0.4. It is again emphasized, that, while ILES ( = 0.00) theoretically should involve the least modelling of turbulence, it is not the best approach. It does yield the highest lift coefficient, but at the same time the predicted drag force is lower than all other results. The skin friction shows that the flow is still laminar, but due to the now absent leading edge separation the lift coefficient is reduced compared to the results with steady inflow. The use of a different convection scheme in the momentum equation might improve these results. It is known from literature that ILES requires a convection scheme which adds enough dissipation to act as a sub-filter model. In the current work QUICK is used, which might not fulfil that requirement by either adding too much, or not sufficient dissipation. An investigation into different convection schemes is outside of the scope of the current work. The difficulties in predicting integral quantities do indicate that the use of ILES is rather sensitive to the setup, next to the entanglement of modelling and discretisation errors, which is inherent to the method . Consequently, in the remainder of this work = 0.25 is employed.

Increasing  versus 
When varying both  and  systematically (cases A, B, C and D), a distinction can be made between the effect of  versus  . These variations are only investigated for = 0.25, on the finest grid, since here the largest difference between results from cases with and without inflow turbulence occur. Varying  has almost no effect on the timeaveraged lift and drag coefficient. However, with a smaller integral length scale, the standard deviation of the signal increases, see Table 3. Increasing  leads to an increase of 5 − 10% in time-averaged quantities; as expected the increase in standard deviation is again significantly larger (between 2 and 3 times larger). The streamwise decay of  (not shown) is as expected comparable for all computations, but varies in magnitude depending on the inflow  .
The flow visualization in Fig. 16 is in line with previous results, modifying  shows only small differences, while increasing  decreases the size of separation region at the trailing edge. This is also clearly visible in Fig. 17, where , increases with  . It is important to realise that the location of transition is not affected by the settings of the turbulence generator, only the magnitude of the skin friction. For  = 2.0%, there is only a difference between ∕ = 0.013 and 0.033 near the trailing edge, whereas for  = 6.3%, ∕ = 0.013 yields a higher skin friction along the chord. The increase in , manifests itself already at the leading edge. All computations show a trough at ∕ ≈ 0.05, but the minimum , varies. For  = 2.0%, for both integral length scales min ( , ) × 10 3 ≈ −1, indicating locally recirculating flow. The lack of difference between the two integral values indicates that while the inflow turbulence intensity is high enough to prevent leading edge separation on the wing, due to the turbulence decay it is too low at the leading edge to further affect the boundary layer flow. When the inflow turbulence intensity is increased to  = 6.3%, min ( , ) increases. When ∕ = 0.033 the minimum is close to 0, while for the smaller integral length scale the minimum is approximately +1. Here the inflow turbulence intensity does affect the boundary layer flow, the smaller integral length is closer in magnitude to the turbulent length scales occurring closer to the wing surface, leading to an increasingly turbulent boundary layer and therefore increased skin friction. For the pressure distribution, again no differences are observed, except near the trailing edge. These observations imply that a sufficiently small integral length scale can trigger transition directly at the leading edge, similarly to roughness applied to the wing. However, such a small length scale must be supported by the grid to be convected from the turbulence generator until the leading edge and enter the boundary layer. Due to the employed grid resolutions in this work, no integral length scales below ∕ = 0.013 are investigated. For the remainder of this work, this smaller length scale is employed.

The effect of increasing inflow turbulence intensity
To investigate the effect of increasing the inflow turbulence intensity, cases B, C, E, F, G and H are compared. The integral length scale is kept fixed to ∕ = 0.013, and only grid 1 with = 0.25 is investigated. The streamwise development of follows the theoretical decay for all cases, and an increase in also leads to an increase in ∕ .
The effect on the integral quantities is shown in Fig. 18. Both the time-averaged lift and drag increase with increasing inflow turbulence, and both quantities increase by approximately 10% between the lowest and the highest inflow turbulence levels. However the difference between a RANS solution and a SRS solution without inflow turbulence is much larger for the drag than for the lift, as seen in previous sections.
So while the increase in inflow turbulence results in a lift force which is comparable to the RANS solution, the drag force is still significantly underpredicted. Note also that for both lift and drag, the increase with increasing  seems to converge, meaning even higher inflow turbulence intensities will likely not result in a significantly higher mean force value. Note that for  = 15.5%, ∕ 2 ∞ = 0.60, which is very high, but it is employed for the sake of completeness. The standard deviation of the signal has not converged. Furthermore, the large value for at  = 0.0% is caused by the separation at the leading edge, so this should not be interpreted as the 'correct' result.
Finally, the effect of increasing inflow turbulence intensity on , and is investigated in Fig. 19. The , profiles show a small separation at the leading edge for non-zero ∕ 2 ∞ , leading to transition. This transition location does not vary with the inflow turbulence intensity, although an increase in  does lead to an increase in skin friction along the wing, which is also observed in RANS computations in literature (Lopes et al., 2020). For , , the results with the highest  seem to converge, indicating a limit to the skin friction which can be obtained using the turbulence generator on this grid, with this value for the integral length scale. The minimum value of , at ∕ ≈ 0.05 also increases, until a limit is reached at approximately , ×10 3 = 1.1. A consequence of the increasing inflow turbulence intensity is a change of the slope of , , in the range 0.1 ≤ ∕ ≤ 0.7. The peak at ∕ ≈ 0.1 reaches a maximum for  = 6.3%, after which it decreases again, while at ∕ ≈ 0.7 the , keeps increasing with increasing  . Finally, again for , there is little difference between the computations with inflow turbulence.

Conclusions based on wetted flow computations with inflow turbulence
Based on these comparisons, the use of synthetic turbulence at the inflow can suppress the large leading edge separation as observed for PANS with a low value and usual steady inflow conditions. Provided that the grid can support the selected integral length scale, a decrease in integral length scales increases the skin friction. The effect of varying  is larger. No inflow turbulence leads to either a large leading edge separation zone, or a laminar flow along the wing. A small amount of inflow turbulence already suppresses the leading edge separation, and leads to a small separation at the leading edge, inducing transition. To increase the lift and drag coefficient, higher inflow turbulence intensities are required.
, shows that even with the highest levels of turbulence attempted in this work, still laminar to turbulent transition is visible downstream of the leading edge. To trigger transition further towards the leading edge, the integral length scale should be reduced. This does require the grids to be refined further, due to the computational costs this was not attempted in the current work. The application of inflow turbulence does also decrease the turbulent separation region at the trailing edge, however it was not completely removed. Separation at the trailing edge is never mentioned in literature, however there is an indication that this phenomenon was also present in the experiments. When we compare one of the few reported wetted flow velocity fields obtained with PIV by Foeth (2008) with the current computations using the same colour scale, it becomes clear that the turbulent separation at the leading edge also occurred in the experiments. Fig. 20 shows the same wavy region of low velocity at the trailing edge, although the velocity scale prohibits the observation of negative velocities at the trailing edge. This implies that turbulent separation at the trailing edge is a feature of this test case, and we should not aim to attempt to remove this by increasing the inflow turbulence intensity.
is only affected near the trailing edge; closer to the leading edge it matches the limited experimental data well. This indicates that the   presence of synthetic turbulence should have little effect on cavitation inception (this will be investigated in the next section). This does not mean that cavitation is not affected at all. The presence of resolved turbulence in the flow can potentially disturb the formed cavity, leading to additional dynamics, varying wing loading or noise. Together with the ability to tune the turbulence inflow statistics, the procedure here used is a promising method to compute noise due to cavitation dynamics, or dynamically varying blade loading for cavitating or non-cavitating propellers in a non-uniform wake field.

Cavitating results
The final investigation in this work concerns the application of inflow turbulence to a case with cavitation. Thus far in literature PANS is rarely combined with multi-phase flows, a notable exception being the work by Ji et al. (2014b). However this work focusses on analysing the cavity dynamics, only a single is attempted with a steady inflow, and the associated modelling errors are not investigated. To the knowledge of the authors the combination of cavitation with inflow turbulence has not been attempted before. The computations are compared against numerical results obtained with DDES and IDDES; and against experimental results by Aït Bouziad (2005) and Foeth (2008), designated EPFL and Delft, respectively.
In the computations we match the cavitation number, = 1.07, while the settings for the inflow turbulence correspond to Case B of the wetted flow computations, meaning ∕ = 0.013 and  = 6.3%; the turbulence intensity is comparable to the experiments. For this level of inflow turbulence, the lift coefficient is relatively insensitive to the selected value. The relatively low turbulence intensities limit the risk of numerical instabilities when cavitation modelling is included. Higher inflow turbulence intensities lead to improved drag predictions, but in cavitating computations can also lead to increased computational instability due to the combination of shedding sheet cavitation and inflow turbulence. The integral length scale in the experiments is unknown, the currently employed value is selected based on numerical reasons. The Courant number is well below 1, with the exception of some cells at the trailing edge where a maximum value of 4 is reached. The computations are started from a wetted flow computation with cavitation introduced over a period corresponding to around 0.5 shedding cycles. Based on the TST it is found that the first five shedding cycles must be removed to eliminate the start-up effects. Computations are then continued for an additional six cycles in the stationary range, reducing the statistical uncertainty for the mean lift and drag forces below 3%. A shedding cycle is 1∕ larger than the wetted flow time unit, so = ∕ . The normalised residuals reach at least 10 −3 for 2 for all equations. The ∞ norms are higher, and occur near the leading edge. Fig. 21 shows the developed time signal and power spectral density (PSD) of the lift coefficient and vapour volume, compared with DDES and IDDES results, following removal of the start-up effects, as discussed in Section 4. The PSD is computed using the pwelch algorithm and applying a Hann window with 50% overlap, resulting in averaging over 6 segments. Note the significantly larger variations in for IDDES. The PSDs for the PANS computations show no clear dominant shedding frequency due to the windowing, but show that the resolved turbulent inflow leads to a higher PSD at higher frequencies.
The integral values can be found in Table 4. As for wetted flow, the application of a resolved turbulent inflow condition increases the lift for the PANS computations; is below the lift predicted by the Hybrid models (DDES and IDDES). When comparing to experiments, all models underpredict the lift. Fig. 1 already showed that this is common in literature, indicating that the discrepancy might also be attributed to unknown experimental uncertainty. It is remarkable that the fluctuations in lift are also significantly smaller for PANS than for the Hybrid models, we will comment on this later. Note that the standard deviation of the lift coefficient is an order of magnitude larger than for the wetted flow computations (Section 5), indicating that the fluctuations in the lift due to the inflow turbulence are negligible compared to the cavitation induced fluctuations. Unsurprisingly, the predicted drag is higher for all models which employ RANS close the wall, the fluctuations in drag however are comparable for all models. For PANS, the standard deviation of the drag has increased by a factor 2 compared to wetted flow, again indicating that the fluctuations due to Table 4 Integral values of cavitating computations with PANS with a steady and a resolved turbulent inflow and = 0.25 on 1. DDES results were previously reported in Vaz et al. (2017).
computed from PSDs without windowing (not shown in the paper). cavitation are dominant. This is confirmed by the observation that the standard deviation for a steady inflow condition is slightly higher than for a resolved turbulent inflow. The average vapour volume increases for PANS when the resolved turbulent inflow condition is employed, but the mean and standard deviation of the vapour volume remain below the predictions done by DDES and IDDES. Finally, the Strouhal number is computed from a single segment PSD of the lift coefficient and vapour volume to emphasize the lower frequencies of the shedding sheet cavity (not shown in the paper). For PANS without inflow turbulence, is larger than , implying different mechanisms are responsible for the vapour and lift fluctuations respectively, and indicating additional dynamics in the sheet cavity. For PANS with a resolved inflow, the shedding frequency halves although it must be remarked that this shedding frequency is difficult to discern. Fig. 22 shows both the instantaneous and time-averaged streamwise skin friction and velocity at the centreline for both cases with PANS. The trailing edge separation can still be observed towards the tunnel wall, where the angle of attack is smaller, but in the regions downstream of the sheet cavity it is absent. The size decreases again when inflow turbulence is added. An interesting observation is that the low momentum inside the cavity leads to a time-averaged velocity field, which looks similar to the time-averaged velocity field for wetted flow with low without inflow turbulence, i.e. when flow separation occurs at the leading edge (see Fig. 6). In this case the presence of cavitation masks the leading edge separation which was observed in the wetted flow computation. Similar observations can be made from the time-averaged skin friction and pressure coefficients at midspan, see Fig. 23. The skin friction coefficient is based on the density of the liquid to enable comparing to the wetted flow results. Again the similarity with wetted flow cases without inflow turbulence is clear in the shape of , (compare for example Figs. 23 and 7). The application of inflow turbulence leads to a higher , in the range ∕ > 0.5, but towards the leading edge (in the cavitation region), little differences occur. This is not surprising, since in this region the local flow is dominated by the presence of the low momentum fluid inside the cavity. In the range 0.25 ≤ ∕ ≤ 0.35, the steady inflow results show more re-circulation than the resolved turbulent inflow results, indicating that the application of inflow turbulence affects the boundary layer in the region where the cavity is periodically being shed.
The time-averaged pressure coefficient clearly shows the presence of cavitation compared to the results obtained for wetted flow. At the leading edge, the pressure coefficient shows the suction peak due to the presence of a stagnation point. The magnitude of this peak depends on the value of 0 chosen in the Schnerr-Sauer cavitation model, but does not affect the dynamic cavitation behaviour .  cavity. When comparing to the experimental results, the cavity length is underpredicted by all numerical approaches. The application of a resolved turbulent inflow has little effect on the mean cavity length, however, the mean cavity length predicted by PANS is lower than the cavity length predicted by DDES or IDDES, related to the higher shedding frequency. This also explains the lower predicted vapour volume. Again note the similarity between for cavitating conditions, with the laminar wetted flow condition on grid 1 which also exhibits a horizontal plateau at − ≈ −1.1 in the range 0 < ∕ ≤ 0.2 (see Fig. 8).
The effect of the resolved turbulent inflow is most visible in the standard deviation of the streamwise skin friction and pressure coefficients near the leading edge. For the skin friction in the region ∕ ≤ 0.2, the standard deviation fluctuates but is at some points almost 30% higher, while the standard deviation of the pressure coefficient increases almost 50% in the same region. Along the remainder of the chord the standard deviation for both PANS computations is similar, and in magnitude, is comparable to the results obtained with DDES and IDDES. However, the chordwise variation is sensitive to the model choice. For IDDES two clear peaks can be observed, the largest is related to the varying length of the sheet, while the peak at the leading edge is due to the cavity detaching from the leading edge and growing again. A somewhat similar behaviour is observed for PANS with a steady inflow. However, for the resolved turbulent inflow reduces at the leading edge, indicating a different shedding behaviour. Instead of the periodically growing and detaching cavity as observed in experiments and obtained with IDDES, for PANS only parts of the sheet are shed, similar to what was observed for DDES . For PANS, this happens at a higher shedding frequency, leading to a lower standard deviation of the lift coefficient and a lower time-averaged vapour volume.

Discussion
In this work, a synthetic turbulence generator is used to reduce the modelling errors in SRS incurred by using a steady (RANS-like) inflow, both in wetted and cavitating conditions. For the test case considered, when applying ≤ 0.5, it is necessary to insert turbulent fluctuations upstream of the wing in order to obtain a flow field predicted by SRS which is close to the experimentally observed flow field. Since lowering reduces the amount of modelled inflow turbulence, this must be replaced by resolved turbulence in order to maintain a physically correct flow field approaching the object of interest. When a steady inflow is applied, the flow around the object remains laminar along most or all of the chord, and -for certain combinations of sufficiently fine grid resolution and low -separates already at the leading edge. Flow separation was not previously reported for numerical results found in literature, although an indication of this phenomenon can be found in the experimental work by Foeth (2008). The reasons for this oversight vary: the results obtained with RANS or Hybrid models do not exhibit flow separation, due to the RANS assumption of a fully-developed turbulent boundary layer. The LES results reported in literature also did not exhibit this, although this could be related to the employed grid resolution at the wall; not just in wall-normal direction, but also in the directions parallel to the wall (which is in literature much coarser than the ones here used). We observed in the current work that the inclusion of cavitation in the computation tends to mask the leading M. Klapwijk et al. Fig. 22. Instantaneous (top row) and time-averaged (bottom row) limiting streamlines and friction coefficient , , on the wing surface, and streamwise velocity, ∕ ∞ , for = 0.25 on grid 1, for steady (left) and a resolved turbulent inflow (right) (∕ = 0.013,  = 6.3%). The instantaneous iso-contour of the vapour volume is indicated in green ( = 0.1). edge separation, which is likely also the case for the LES results from literature.
Using synthetic inflow turbulence, transition is triggered closer to the leading edge, resulting in turbulent flow along the wing, thereby avoiding laminar flow separation. This behaviour is reminiscent of the application of vortex generators on aerofoils, where the vortex generators energize the boundary layer, thereby improving the stall behaviour by preventing leading edge separation (Clancy, 2006;Rist and Augustin, 2006;Sreejith and Sathyabhama, 2018). The turbulent flow separation at the trailing edge seen in certain simulations is, however, thought to be physically correct, since this was also observable in the experiments. It is noted that the sensitivity to inflow turbulence quantities for SRS is analogous to what is commonly observed when using RANS with transition models, where modest differences in inflow turbulence intensity can dramatically affect boundary layer development (see e.g. Liebrand et al., 2019Liebrand et al., , 2020Lopes et al., 2020).
An alternative approach for obtaining a turbulent boundary layer in the simulations is to trip the boundary layer close to the leading edge. This is commonly applied for measurements of turbulent boundary layers performed at moderate Reynolds numbers (10 4 -10 5 ) or for flow control purposes, in order to reduce the size of laminar separation bubbles, and reduce the associated drag. Applications include aircraft wings, wind turbines or blades of turbomachinery (Rist and Augustin, 2006;Sreejith and Sathyabhama, 2018). The sand grain roughness applied in the experiments of the present test case could be reproduced numerically in a number of ways: by geometrically resolving the roughness (Asnaghi and Bensow, 2020), applying a simplified trip in the geometry (e.g. Sreejith and Sathyabhama, 2018;Winkler et al., 2020), or by adding wall-normal velocity fluctuations (e.g. Rist and Augustin, 2006;Schlatter and Örlü, 2012;Wolf et al., 2012). An advantage of a trip is that there is no need to tune the turbulence generator to obtain a turbulent boundary layer. Nevertheless, there are several challenges involved in applying this type of approach. Geometrically resolved sand grain roughness is far from standard practice in CFD, and will in any case only quantitatively agree with what is used in experiments, while wall-normal velocity fluctuations and a trip -despite the much simpler geometry -must also be tuned to obtain a flow disturbance equivalent to the roughness applied in the experiments. Wall-normal velocity fluctuations could also have a detrimental effect on iterative convergence. Although these approaches do reduce the need for a refined grid upstream of the object, leading to slightly reduced numerical cost per simulation, the required modification to the grid makes it a less general methodology, which is more difficult to apply across a range of test cases or operating conditions (Sreejith and Sathyabhama, 2018). In real-life engineering applications it can be expected that inflow turbulence is present. Consequently, the use of a steady inflow in combination with SRS always introduces a certain level of modelling error. When trips are used, there is also the risk that a too thick turbulent boundary layer develops, which may affect the flow prediction downstream. It can also introduce a pressure jump (Sreejith and Sathyabhama, 2018), which can have implications for numerical cavitation inception (which is typically based on the simplified criterion < − ). Therefore, when the development of a turbulent boundary layer is the primary goal, applying leading edge roughness is probably the better approach. However, should the focus of simulations be more on large-scale dynamics, such as noise generation of propellers due to inflow turbulence (Lloyd et al., 2014) or cavitation dynamics, the development of the boundary layer is of less importance compared to the interaction of turbulence with the leading edge. In such cases the inflow turbulence generator is the more appropriate choice.
It was shown in this work that, despite specifying unphysically high inflow turbulence intensities in the simulations (significantly higher than the reported 2-3% in the experiments), it was not possible to obtain transition directly at the leading edge, since the inflow turbulence fluctuations do not enter the boundary layer due to insufficient (spanwise) grid resolution. This agrees with the findings of Tangermann and Klein (2020), who studied the effect of inflow turbulence for a NACA profile at moderate angle of attack using DDES. While using inflow turbulence with similar ∕ to that used in the present work, they observed that smaller integral length scales are needed to enter the boundary layer and cause transition, while larger length scales mainly affect loading, as they induce angle of attack fluctuations. However, to be able to properly resolve a fully turbulent boundary layer from the leading edge onwards, the grid resolution has to be increased by several orders of magnitude to support these smaller length scales, as well as the transition process they initiate. The lack of resolution leads to a limit in skin friction: for computations with low values, the skin friction remains significantly lower than the results with higher values. It is noted, however, that when applying a trip it might also be difficult to numerically reproduce the transition behaviour behind the trip due to the fine grid resolution required.
As already mentioned in the Introduction, the Delft Twist 11 Hydrofoil was used as a test case representative of a simplified propeller. The results in this work indicated that the fundamental shedding frequency of the sheet cavity was not affected by the synthetic inflow turbulence. For a propeller operating in a wake field, this finding likely holds, since the dynamics of a sheet cavity at the blade passage frequency are governed by the rotation of the blade in a non-uniform mean velocity field, although cavity shedding and higher-frequency dynamics are expected to be somewhat sensitive to inflow turbulence. Simulating a non-uniform wake field in SRS leads to the additional requirement of specifying inhomogeneous anisotropic inflow turbulence. This is already possible using the current methodology , making the extension to engineering applications feasible. The current setup is also well suited to investigate tip vortex cavitation dynamics for wings or propellers, where the inflow turbulence can perturb the cavitating vortex (Bosschers, 2018). Using this approach, underwater radiated noise could also be assessed.

Conclusion
The Delft Twist 11 Hydrofoil was evaluated in wetted and cavitating flow conditions, using the PANS methodology with a varying ratio of modelled-to-total kinetic energy, on different grids and using both a steady and a resolved turbulent inflow. It was shown that modelling errors can occur when applying a steady inflow condition, as used for RANS, due to the flow around the object remaining laminar. This leads to underpredicted lift and drag forces, making validation impossible. Depending on the grid resolution, the wing can even exhibit laminar leading edge separation, i.e. stall.
When employing a resolved turbulent inflow, where homogeneous isotropic turbulence is inserted into the flow, the leading edge separation disappears regardless of the chosen input values for the turbulence generator. The region of trailing edge separation is reduced by increasing the inflow turbulence intensity, which also increases the skin friction along the wing. The separation at the trailing edge is also observed in the experimental results, meaning that, although this was not reported previously, we do obtain a proper match with the experimental flow field. Increasing the inflow turbulence intensity also increases the mean lift and drag force, though only till a certain limit. The fluctuations in lift and drag keep increasing. An important observation is that the inflow turbulence has almost no effect on the pressure distribution. This implies that inception remains unaffected when using the in CFD commonly applied cavitation inception criterion. Varying the integral length scale has less effect on the flow separation. When using a resolved turbulent inflow, it is necessary that SRS grids are refined upstream to support the convection and development of inflow turbulence. For cavitating conditions, the variations in lift and drag due to inflow turbulence are significantly smaller than the fluctuations due to the shedding sheet cavity. The addition of inflow turbulence does affect the shedding behaviour: smaller parts of the sheet cavity are being shed at a higher frequency. The predicted mean lift and drag match the numerical results reported in the literature, while the predicted Strouhal number is higher than the experimental value. To the knowledge of the authors, this is the first application of synthetic inflow turbulence to a test case including cavitation.
As part of future work it is desirable to improve the iterative convergence of the computations, which was seen to reduce due to the interaction of cavitation and inflow turbulence. Secondly, the influence of applying a symmetry plane at the foil midspan in order to reduce the domain size (and thereby also the required number of grid cells) should be further investigated. When resolving turbulence, the symmetry plane introduces a modelling error, and with the inclusion of inflow turbulence this error becomes even larger, being already present upstream of the wing. It is also clear that further refined grids are a necessity to obtain a proper resolved boundary layer and properly capture transition behaviour with SRS methods. With respect to the validation, the inclusion of more physics in the CFD computations leads to more stringent demands for experimental data. It is recommended to not only measure integral quantities, but also characteristics of the setup such as turbulence intensity and integral length scale, and the flow in the boundary layer. The absence of this information for the Delft Twist foil highlights the need for new experimental test cases aimed at validation of SRS investigations of multiphase flows.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.