1. Introduction
The problem of a vortex interacting with a wing or an appendage is common in many applications from aerodynamics and naval hydrodynamics. Different forms of such an interaction have been the topic of numerous experimental and computational studies in the literature (see, for instance, Lee, Burggraf & Conlisk Reference Lee, Burggraf and Conlisk1998; Rockwell Reference Rockwell1998; McKenna, Bross & Rockwell Reference McKenna, Bross and Rockwell2017; Yang, Veldhuis & Eitelberg Reference Yang, Veldhuis and Eitelberg2017a; Yang et al. Reference Yang, Zhou, Sciacchitano, Veldhuis and Eitelberg2017b). In general, a vortex impacting a wing or hydrofoil experiences a deformation and a cross-stream shift, due to both potential and viscous effects. In the case of an orthogonal collision these phenomena are coupled to a modification of the diameter of the vortex core, which is different on the pressure and suction side branches. In addition, its interaction with the boundary layer on the body can trigger transition or separation of the latter and the generation of additional structures, influencing this way the evolution and the stability of the primary colliding vortex itself. If the vortex is strong enough to induce separation of the boundary layer and ejection of its vorticity, the latter will significantly disrupt the primary vortex as well, even before its actual collision with the foil. It is worth noting that the above phenomena are of great practical importance, since significant forces and pressure fluctuations are generated by the impact of the vortex on the body and its ability to generate a disturbance on the boundary layer. As a consequence, they impact significantly the structural integrity and the acoustic signature of the overall system (Roger, Schram & Moreau Reference Roger, Schram and Moreau2014; Quaglia et al. Reference Quaglia, Léonard, Moreau and Roger2017; Yang et al. Reference Yang, Veldhuis and Eitelberg2017a; Wang et al. Reference Wang, Guo, Xu and Su2019).
Investigating these interactions via computations and/or experiments is challenging, and a number of early studies exploited analytical or simplified numerical methodologies in idealized configurations: vortex filaments techniques and inviscid vortex assumptions were utilized in a number of works (see, for example, Marshall Reference Marshall1994; Marshall & Yalamanchili Reference Marshall and Yalamanchili1994; Marshall & Grant Reference Marshall and Grant1996; Marshall & Krishnamoorthy Reference Marshall and Krishnamoorthy1997; Marshall Reference Marshall2002). For a columnar vortex impacting orthogonally on a wing or a cylinder, for example, a decrease of the core radius on the suction side and an increase on the pressure side were observed. In particular, Marshall & Krishnamoorthy (Reference Marshall and Krishnamoorthy1997) found that the impulsive cutting of a vortex with a non-zero axial velocity produces phenomena resembling those occurring in the motion of a piston within a tube. On the suction side of the vortex an ‘expansion wave’ is produced and the radius of the vortex core experiences a decrease in the vicinity of the cutting surface, with rings of tangential vorticity moving away from the surface along the direction of the axial velocity. The behaviour of the compression side of the vortex was instead found to be a function of its axial velocity and how it compares to the speed of the waves of tangential vorticity within the vortex core. Based on that, they identified subcritical and supercritical cases: in the former, rings of tangential vorticity propagate upstream in the direction opposite to the axial velocity within the vortex core, with an increase of the vortex radius near the cutting surface; in the latter, the upstream propagation of tangential vorticity cannot occur and the pressure side branch of the vortex is similar to an impinging non-swirling jet on a flat plate. Experiments carried out in the same configuration demonstrated that the assumption of impulsive cutting is acceptable only when the circulation of the vortex is small enough to avoid separation of the boundary layer. Marshall (Reference Marshall2002) focused on the scenario of breakdown of the inviscid assumption, featuring boundary layer separation and ejection of secondary vorticity, influencing pressure fluctuations and the evolution and stability of the primary vortex. In this scenario the secondary vorticity wraps around the primary vortex, promoting additional vorticity ejection from it and substantially increasing the complexity of the flow. He developed analytical models aimed at representing the behaviour of the secondary vorticity field in presence of orthogonal vortex–blade interaction. The above results were also supported by the measurements by Kim & Komerath (Reference Kim and Komerath1995) and Green, Coton & Early (Reference Green, Coton and Early2006) in similar configurations. In particular, the particle imaging velocimetry (PIV) experiments by Green et al. (Reference Green, Coton and Early2006) demonstrated a decrease of the axial velocity within the vortex core on the pressure side of the blade and an increase on the suction side, consistent with the expansion and contraction of the size of the vortex core, respectively. On the pressure side a realignment of the velocity and vorticity vectors occurred, missing instead on the suction side.
Higher fidelity computations, able to account for the boundary layer–vortex interactions identified above, have also been reported recently, although still for simplified configurations. For instance, Gordnier & Visbal (Reference Gordnier and Visbal1999) and Garmann & Visbal (Reference Garmann and Visbal2015) utilized implicit large eddy simulation (LES) at values of the Reynolds number of order 
$ {O}(10^{4})$ in canonical flow problems involving a delta wing vortex and a streamwise-oriented vortex, respectively, impinging on a flat plate. In the former case the interaction was the source of separation phenomena, formation of hairpin vortices and spiral vortex breakdown, with pressure fluctuations on the plate corresponding to the modes of the vortex breakdown. Liu & Marshall (Reference Liu and Marshall2004) solved the full Navier–Stokes equations for Reynolds numbers of order 
$ {O}(10^{3})$, considering the orthogonal collision of a columnar vortex on an airfoil. They found that both the cutting process and the lift force on the airfoil, generated by asymmetries of the impinging vortex arising on the suction and pressure sides, were almost independent of the Reynolds number. They also showed that the process of vortex cutting is always incomplete and the colliding vortex keeps wrapping around the leading edge of the airfoil. The presence of an axial velocity within the vortex core is the source of the asymmetries observed on the suction and pressure sides, producing a lift force, in agreement with results from earlier inviscid computations. Recent computations by Saunders & Marshall (Reference Saunders and Marshall2015, Reference Saunders and Marshall2017), in a similar configuration, identified three phases of the vortex cutting process: (i) vorticity is generated at the leading edge of the blade, due to the velocity induced by the approaching vortex; such vorticity is oriented in the opposite direction, compared to that within the vortex core; (ii) vorticity diffusion and destruction occur between the vortex core and the boundary layer on the leading edge of the blade, due to the opposite orientation of vorticities; (iii) the process of vorticity destruction within the vortex core stops, due to reorientation of the leading edge vorticity within the boundary layer; as a consequence the primary vortex keeps wrapping around the leading edge, stretching while moving downstream across the blade.
All the above studies provide a glimpse into the complex physics of vortex/body interactions in canonical configurations. More complex configurations, which are typical in many practical applications, have not been considered at least in this context. A characteristic example is that of helical vortices, which are typically shed from a rotating blade, and are robust coherent structures defining the wake of the rotating foil. For the case of a propeller or a wind turbine these vortices may co-exist and interact with each other and with additional vortical structures (for instance, the hub vortex), as well as the wake of the blades. In many situations downstream bodies are also present (such as a rudder in the near wake of a propeller), which is similar to the canonical configurations reviewed above. In this case, however, their interaction with the boundary layer on the foil is complex, leading to a faster transition and triggering higher levels of turbulent fluctuations near its surface, compared to canonical scenarios where an isolated vortex collides with a body within a laminar flow.
We should note that this interaction is of major practical significance, due to the fatigue loads induced on the foil and the vibrations and erosion associated with cavitation. This is documented in a number of studies correlating the stress and damage of marine rudders to the trajectory and intensity of the main coherent structures in the wake of propellers (Van Terwisga et al. Reference Van Terwisga, Fitzsimmons, Ziru and Foeth2009; Choi, Kim & Lee Reference Choi, Kim and Lee2010; Paik et al. Reference Paik, Kim, Kim, Park, Heo and Yu2010, Reference Paik, Kim, Kim, Kim and Suh2012; Rhee et al. Reference Rhee, Lee, Lee and Oh2010; Felli, Falchi & Pereira Reference Felli, Falchi and Pereira2011).
Felli, Camussi & Guj (Reference Felli, Camussi and Guj2009) and Felli & Falchi (Reference Felli and Falchi2011) conducted detailed experiments for the case of a naval propeller (the INSEAN E779A), operating upstream of a hydrofoil (a NACA0020) with and without offset relative to the propeller axis, respectively. Using laser Doppler velocimetry and PIV techniques they mainly focused on the visualization and dynamics of the tip/hub vortices and reported statistics in the near wake. In particular, Felli & Falchi (Reference Felli and Falchi2011) found that the interaction of the hub vortex with the boundary layer on the hydrofoil had a profound impact on the turbulent kinetic energy across the rudder and the near wake. A significant spanwise shift of the hub vortex, opposite to that experienced by the tip vortices, was also reported. These results were confirmed by Felli et al. (Reference Felli, Falchi and Pereira2011), who measured the pressure fluctuations on the surface of the hydrofoil. They were significantly higher in the area of the rudder populated by the hub vortex than across the trajectory of the tip vortices. In contrast, Felli et al. (Reference Felli, Camussi and Guj2009), studying a configuration where the axis of the propeller was not aligned with the downstream hydrofoil, reported no spanwise shift of the hub vortex and a smaller, albeit still significant, area of high turbulence populating the axis of the propeller wake. This can be attributed to the offset between the hub vortex and the hydrofoil, resulting in reduced shear between the hub vortex and the boundary layer on the rudder. These experiments point out to the complexity of the flow in such systems and motivated a series of computations, where eddy resolving techniques were adopted to capture the dynamics of the dominant coherent structures.
Muscari, Dubbioso & Di Mascio (Reference Muscari, Dubbioso and Di Mascio2017), for example, utilized detached eddy simulation (DES), coupled with a Chimera technique, to compute the wake of a simplified, single-bladed, INSEAN E779A propeller operating upstream of a NACA0020 rudder. The adopted numerical resolution was relatively coarse, utilizing 
$\sim$15 million finite volumes across 160 overlapping blocks. They observed the formation of hairpin vortices on the surface of the rudder, due to the interaction with the hub vortex, contributing to reconnection between hub and tip vortices downstream of the hydrofoil. In agreement to the above experiments, a significant spanwise shift of the tip vortices was produced across the rudder, with the pressure and suction side branches moving outward and inward, respectively. Reconnection was also observed between the pressure and suction side branches of the tip vortices. DES computations were also reported by Wang et al. (Reference Wang, Guo, Xu and Su2019), who considered the same propeller–rudder configuration as in Felli & Falchi (Reference Felli and Falchi2011). An unstructured mesh composed of a total of 26 million finite volumes was utilized in the framework of an overset grid technique. They demonstrated an increasing complexity of the wake topology at higher propeller loads and with larger incidence angles of the rudder. Higher loads promoted stronger vortices shed by the propeller, a more intense vorticity exchange between tip vortices and rudder boundary layer and the formation of larger structures connecting the tip vortices to each other. Those structures originated from the interaction of the tip vortices with the rudder boundary layer. More recently Hu et al. (Reference Hu, Zhang, Sun and Guo2019) reported LES results for the same configuration. They were in broad qualitative agreement with the reference experiment by Felli & Falchi (Reference Felli and Falchi2011) and the DES computations by Muscari et al. (Reference Muscari, Dubbioso and Di Mascio2017). All the above computations mainly focused on the topology and dynamics of the dominant coherent structures. The effect of this interaction on turbulence within the wake, as well as the rudder boundary layers, documented by detailed statistics, was not analysed, maybe due to the lack of proper spatial and/or temporal resolution. This information is especially important in many fronts: (i) validation of lower fidelity modelling strategies; (ii) prediction of cavitation; (iii) definition of the acoustic signature of the overall propeller–rudder system. We should note that in the naval hydrodynamics context the latter depends not only on the pressure fluctuations on the surface of bodies, as typical in aerodynamics applications, but especially on wake turbulence. Ianniello (Reference Ianniello2016) indicated that this is due to the quick decay of the surface acoustic terms for low values of Mach number, making the nonlinear volume terms in the near wake the main source of noise.
High resolution LES can certainly close this gap and provide a detailed characterization of the complex turbulent wakes generated in configurations such as the above. Several examples can be found for the case of propellers in open-water conditions (i.e. Di Felice et al. Reference Di Felice, Felli, Liefvendahl and Svennberg2009; Liefvendahl Reference Liefvendahl2010; Liefvendahl, Felli & Tröeng Reference Liefvendahl, Felli and Tröeng2010; Lu, Bensow & Bark Reference Lu, Bensow and Bark2014; Balaras, Schroeder & Posa Reference Balaras, Schroeder and Posa2015; Asnaghi, Svennberg & Bensow Reference Asnaghi, Svennberg and Bensow2018, Reference Asnaghi, Svennberg and Bensow2020; Posa et al. Reference Posa, Broglia, Felli, Falchi and Balaras2018, Reference Posa, Broglia, Felli, Falchi and Balaras2019b), but to the best of our knowledge, none for propellers interacting with downstream hydrofoils, probably due to cost considerations. In this work we report LES of the INSEAN E1658 propeller operating upstream of a NACA0020 hydrofoil, mimicking a rudder. The level of resolution of the computational grid goes beyond what is currently available in the literature. This was achieved via a highly efficient structured solver in cylindrical coordinates, where the boundary conditions on the rotating propeller and the downstream rudder were imposed using an immersed-boundary (IB) formulation. The overall method is a good fit for the particular configuration, which can be tackled on a cylindrical coordinates grid at a fraction of the cost of classical overset boundary conforming approaches. The present computations capture in great detail the physics of the interactions involving tip and hub vortices and the boundary layer on the surface of the rudder and – most important – how this is reflected in the turbulence statistics across the wake. In particular, we will demonstrate that the hub vortex is actually the main source of turbulence across the hydrofoil and experiences a significant spanwise displacement, which is opposite to that of the tip vortices, resulting in high levels of shear on the suction sides of the hydrofoil. This interaction was found to produce a substantial increase of turbulence levels, both across the hydrofoil itself and the near wake, leading to a non-monotonic evolution of turbulent stresses as one moves away from the propeller–rudder system. This is the first study to clearly identify the importance of the hub vortex in the overall wake dynamics for this class of flows via detailed turbulence statistics.
The paper is structured as follows: methodology (§ 2), computational set-up (§ 3), validation and analysis of the flow physics, pointing to the interaction of the wake structures from the propeller with the downstream hydrofoil (§ 4) and conclusions (§ 5).
2. Methodology
This study was carried out via LES, solving the filtered Navier–Stokes equations for incompressible flows in non-dimensional form:
\begin{gather} \mathcal{D}( \tilde{\boldsymbol{u}}) = 0, \end{gather}
\begin{gather}\frac{\partial \tilde{\boldsymbol{u}} }{\partial t} + \mathcal{D}( \tilde{\boldsymbol{u}} \tilde{\boldsymbol{u}}) = - \mathcal{G}( \tilde{p}) - \mathcal{D}(\boldsymbol{\tau}) + \frac{1 }{Re} \mathcal{L}( \tilde{\boldsymbol{u}}) + \boldsymbol{f}, \end{gather}
where 
$\tilde {\boldsymbol {u}}$ and 
$\tilde {p}$ are the resolved (filtered) velocity vector and pressure, respectively, with the filter operator indicated by 
$\sim$ and implicitly defined by the adopted computational grid, whose local spacing defines the filter size; 
$\mathcal {D}(\cdot )$, 
$\mathcal {G}(\cdot )$ and 
$\mathcal {L}(\cdot )$ stand for the divergence, gradient and Laplace operators, respectively. The first term of (2.2) is the time derivative of 
$\tilde {\boldsymbol {u}}$. The Reynolds number is defined as 
$Re=ul/\nu$, where 
$u$ is a reference velocity, 
$l$ a reference length and 
$\nu$ the kinematic viscosity of the fluid. The quantity 
$\boldsymbol {\tau }=\widetilde {\boldsymbol {uu}}- \tilde {\boldsymbol {u}} \tilde {\boldsymbol {u}}$ is the subgrid-scale (SGS) stress tensor, originating from filtering the nonlinear terms of the Navier–Stokes equations and representing the effect of the small, unresolved scales on the large, resolved scales of the flow. In this work the SGS stresses are parameterized using the wall-adapting local eddy-viscosity (WALE) model by Nicoud & Ducros (Reference Nicoud and Ducros1999). The WALE model computes the eddy viscosity using the square of the velocity gradient tensor, which takes into account both the deformation and rotation tensors of the resolved velocity field. The accuracy and robustness of the WALE model as implemented in the solver utilized in the present study has been established in a variety of complex flows, from naval hydrodynamics (Balaras et al. Reference Balaras, Schroeder and Posa2015; Posa & Balaras Reference Posa and Balaras2016, Reference Posa and Balaras2018, Reference Posa and Balaras2020; Posa, Broglia & Balaras Reference Posa, Broglia and Balaras2019a) to turbomachinery (Posa, Lippolis & Balaras Reference Posa, Lippolis and Balaras2015, Reference Posa, Lippolis and Balaras2016; Posa & Lippolis Reference Posa and Lippolis2018, Reference Posa and Lippolis2019) applications. Posa & Balaras (Reference Posa and Balaras2018), for example, reported LES of a self-propelled submarine and verified that the model accurately reproduces the proper near wall limiting behaviour in complex geometries and non-equilibrium boundary layers, without requiring special corrections. This feature is inherent to the WALE model, as discussed in detail by Nicoud & Ducros (Reference Nicoud and Ducros1999).
The quantity 
$\boldsymbol {f}$ in (2.2) is a volume force, utilized to enforce the boundary conditions on solid bodies via the IB methodology proposed by Balaras (Reference Balaras2004) and Yang & Balaras (Reference Yang and Balaras2006). This forcing function is computed at the interface nodes of the computational grid, which are those outside of the body, but having at least one interior neighbour, and is zero everywhere else. Linear velocity reconstructions are utilized to compute 
$\boldsymbol {f}$, while the ‘field extension’ strategy proposed by Yang & Balaras (Reference Yang and Balaras2006) is adopted to address the problem of grid cells moving from the solid to the fluid regions, due to the motion of the propeller. It is worth noting that in the framework of the IB methodology, also the eddy-viscosity needs to be properly reconstructed in the vicinity of bodies, as discussed by Balaras (Reference Balaras2004), since its computation by SGS models relies on the derivatives of the velocity field, which are in general not reconstructed at the interface nodes.
All spatial derivatives are approximated using second-order central differences on a staggered grid. The advancement in time is carried out using an exact projection method (Van Kan Reference Van Kan1986), where the convective, viscous and SGS terms along the radial and streamwise directions are advanced with an explicit three-step Runge–Kutta scheme. All azimuthal derivatives are treated implicitly with the second-order Crank–Nicolson scheme to relax the stability constraints at the axis of the cylindrical grid. The resulting Poisson equation, utilized to enforce the divergence-free condition, is solved by a direct approach. In particular, the eptadiagonal system is converted into a series of pentadiagonal problems, one for each meridian slice of the cylindrical grid, using fast Fourier transforms along the azimuthal direction. Each pentadiagonal system is then solved using the divide-and-concur algorithm developed by Rossi & Toivanen (Reference Rossi and Toivanen1999).
3. Geometry and computational set-up
The propeller geometry in this work is the same as that considered in earlier studies by Posa et al. (Reference Posa, Broglia, Felli, Falchi and Balaras2018, Reference Posa, Broglia, Felli, Falchi and Balaras2019b) for the INSEAN E1658 propeller in isolated conditions. The reported global performance parameters and flow fields in the wake were in very good agreement with the experimental results by Felli & Falchi (Reference Felli and Falchi2018). Posa et al. (Reference Posa, Broglia, Felli, Falchi and Balaras2018, Reference Posa, Broglia, Felli, Falchi and Balaras2019b) analysed values of the advance coefficient equal to 
$J=0.56$, 
$J=0.65$ and 
$J=0.74$ (
$J=U_\infty /nD$, where 
$U_\infty$ is the free-stream velocity, 
$n$ the propeller rotational frequency and 
$D$ its diameter). Here the same values of advance coefficient are considered, with the purpose of enabling comparisons with the wake flow of the same propeller in isolated conditions. The INSEAN E1658 was placed upstream of a NACA0020 hydrofoil of infinite spanwise extent, in order to avoid end effects (see figure 1). The Reynolds number relative to the propeller is
\begin{equation} Re_p = \frac{c^{b}_{0.70}\sqrt{(\varOmega 0.70R)^{2}+U_\infty^{2}} }{\nu},\end{equation}
where 
$c^{b}_{0.70}$ is the chord length of the propeller blades at 
$70\,\%$ of their radius and 
$\varOmega =2 {\rm \pi}n$ is the angular speed of the propeller. 
$Re_p$ was set to 
$Re_p (J=0.56) \approx 360\ 000$, 
$Re_p (J=0.65) \approx 310\ 000$ and 
$Re_p (J=0.74) \approx 280\ 000$. The Reynolds number relative to the hydrofoil, 
$Re_h$, is
\begin{equation} Re_h = \frac{c^{h} U_\infty }{\nu},\end{equation}
where 
$c^{h}$ is the chord of the hydrofoil; 
$Re_h$ was 
$Re_h \approx 500\ 000$. It is useful to point out that 
$Re_h$ was constant across cases, since 
$J$ was modified by changing the propeller rotational speed. The choice of the set-up, in terms of location of the hydrofoil downstream of the propeller and its chord length, was based on the earlier experimental study by Felli & Falchi (Reference Felli and Falchi2011). However, it should be noted that they considered a different propeller geometry, the four-bladed INSEAN E779A, whose wake was studied via experiments conducted at values of advance coefficient 
$J=0.88$ and Reynolds number, based on the chord of the rudder and the free-stream velocity, 
$Re_h=1.36 \times 10^{6}$. The hydrofoil was aligned with the free stream and the axis of the upstream propeller. Its leading edge was placed at 
$0.5D$ from the propeller plane and its chord length was equal to 
$c=0.833D$. This configuration was selected, since typical of ships in straight-ahead conditions, in terms of both relative size and position between propeller and rudder.
Set-up of the present flow problem: the INSEAN E1658 propeller placed upstream of a NACA0020 hydrofoil of infinite spanwise extent; 
$D$ is the diameter of the propeller.
The same cylindrical coordinates grid, composed of 900, 1026 and 2050 (1.9 billion) nodes along the radial, azimuthal and streamwise directions, was utilized across the three load conditions considered in this study. A meridian slice of the grid is shown in figure 2. Compared to the grid adopted for the open-water computations involving the same propeller (see Posa et al. Reference Posa, Broglia, Felli, Falchi and Balaras2018, Reference Posa, Broglia, Felli, Falchi and Balaras2019b), the present one was kept unchanged within the propeller region as well as upstream. On the wake side of the propeller the resolution was increased to resolve the boundary layer on the surface of the downstream hydrofoil and its interaction with the propeller wake. Specifically, 560 points were distributed across the chord of the hydrofoil, while the grid was also refined at outer radii, beyond the tip of the propeller blades, adding 100 nodes along the radial direction between 
$r/D=0.6$ and 
$r/D=1.0$. It will be shown in the results section that the tip vortices on the pressure sides of the hydrofoil experience a chordwise shift towards outer radii, leading to a radial expansion of the region of the computational domain requiring finer levels of resolution. The Kolmogorov length scale, 
$\eta$, was estimated and the grid was verified fine enough for a well-resolved LES. For instance, in the region of the hydrofoil the grid step across the radial direction was 
${\rm \Delta} r/\eta \approx 18$ and the one across the axial direction was 
${\rm \Delta} z/\eta \approx 34$. Of course, the resolution across the azimuthal direction was a function of the radial coordinate, with the azimuthal spacing ranging from 0 at 
$r/D=0.0$ to 
$r{\rm \Delta} \vartheta /\eta \approx 70$ at 
$r/D=0.5$. Another useful quantity to verify the suitability of the adopted resolution is the ratio between the eddy-viscosity and the molecular viscosity, 
$\nu _t/\nu$. In the region between the propeller plane (
$z/D=0.0$) and 
$z/D=2.3$, encompassing also the downstream hydrofoil, spanning from 
$z/D=0.5$ to 
$z/D=1.33$, the ratio 
$\nu _t/\nu$ kept within 5, reaching values still lower than 10 at downstream locations, due to grid coarsening. Additional evidence of the quality of the computational grid is provided in figure 3, where the power spectral densities for the time history of the streamwise velocity component are shown, referring to probes located on the hydrofoil plane (
$x/D=0$) at spanwise coordinates 
$y/D=0$ and 
$y/D=0.5$ on the left and right sides, respectively, and streamwise coordinates 
$z/D=0.45$, 
$z/D=1.40$ and 
$z/D=2.00$ from top to bottom. Frequencies on the horizontal axes were scaled using the frequency of the passage of the propeller blades, 
$f_{bp}$. In front of the hydrofoil (first row) the spectra at the wake axis (left) follow closely the 
$-5/3$ slope of the Kolmogorov cascade and no dominant frequency can be identified. This is not the case in the vicinity of the tip vortices (right), where the main frequency is that of the blade passage and the energy is concentrated at low frequencies, since the tip vortices are still stable. Due to their coherence there, turbulence is still very far from homogeneity and isotropy. Just downstream of the hydrofoil (central row) the slope of the power spectral density is very close to that of the Kolmogorov cascade both at the wake axis and at the outer radii, where no dominant frequency is visible anymore. Further downstream, the panels in the bottom row of figure 3 still demonstrate the presence of the inertial range of frequencies, although in this region the computational grid is already coarser than across the hydrofoil. The extent of the computational domain was enlarged, compared to the earlier open-water computations, for which we verified the first- and second-order statistics insensitive to the location of the boundaries of the domain. The inflow and lateral boundaries were kept at 
$2.5D$ and 
$5.0D$ from the propeller plane and the propeller axis, respectively, while the outflow boundary was shifted downstream, from 
$5.0D$ to 
$6.5D$, relative to the streamwise position of the propeller plane, taking into account that the trailing edge of the hydrofoil locates at 
$z/D\approx 1.33$. On the inflow, lateral and outflow boundaries uniform velocity, homogeneous Neumann and convective conditions were enforced, respectively. The Lagrangian grid defining the propeller geometry was identical to that adopted by Posa et al. (Reference Posa, Broglia, Felli, Falchi and Balaras2018, Reference Posa, Broglia, Felli, Falchi and Balaras2019b) and was composed of 
$\sim$166 000 triangles, whereas the Lagrangian grid for the rudder was built using 
$\sim$160 000 triangles. Note that the hydrofoil does not extend all the way to the boundary (see figure 2a), which is tied to the current implementation of the IB method. The velocity adjusts to its free-stream value over a few grid cells, while the results in the region of interest are not affected.
Meridian slice of the computational grid. Only 1 of every 64 nodes shown in the 
$yz$ plane. Dotted-dashed line representing the location of the axis of the grid.
Power spectral densities of the time history of streamwise velocity at probes placed on the plane of the hydrofoil (
$x/D=0$) at the spanwise coordinates 
$y/D=0$ (a,c,e) and 
$y/D=0.5$ (b,d,f) and the streamwise coordinates 
$z/D=0.45$ (a,b), 
$z/D=1.40$ (c,d) and 
$z/D=2.00$ (e,f). The dashed line in each panel standing for the 
$-5/3$ slope of the Kolmogorov cascade. Horizontal axes scaled by the frequency of the passage of the propeller blades, 
$f_{bp}$.
Computations were carried out with a constant Courant–Friedrichs–Lewy number, CFL=1, to meet the stability requirements of the three-step Runge–Kutta scheme adopted for the approximation of the explicit terms. This resulted in a very high resolution in time, equivalent to a rotation of approximately 
$0.057^{\circ }$ per time step. For each case of advance coefficient the simulations were run for two flow-through times to remove the memory of the initial transient. Some probes were placed within the wake, where first- and second-order statistics were computed during the simulations, to verify the establishment of statistically steady conditions. Then, statistics of the flow were sampled during 28 additional propeller rotations across the whole computational domain. First- and second-order statistics were compared with those computed during the first half of the same period and differences were found usually within a few per cent.
All computations were carried out on 
$1024$ Intel Xeon Phi7250 (KnightLandings) cores at 1.4 GHz of Marconi KNL at CINECA, Italy, requiring 
$\sim$10 million CPU hours per each case of advance coefficient.
4. Results
Results will be presented in this section as ensemble-averaged and phase-averaged statistics. Ensemble averages were computed on the fly as time averages from all available instantaneous realizations of the flow, after initial transients were removed
\begin{equation} \bar{g} ( \boldsymbol{x} ) = \frac{1 }{N} \sum^{N}_{i=1} g ( \boldsymbol{x}, t_i ),\end{equation}
where 
$g (\boldsymbol {x}, t )$ stands for any physical quantity, 
$\bar {g} ( \boldsymbol {x} )$ for its ensemble average in time, 
$N$ for the number of instantaneous realizations of the flow after reaching statistically steady conditions, 
$\boldsymbol {x}$ for the position vector of any point in space and 
$t$ for time. The ensemble-averaged root mean squares were therefore computed as
\begin{equation} \bar{g}' ( \boldsymbol{x} ) = \sqrt{\frac{1 }{N} \sum^{N}_{i=1} [ g ( \boldsymbol{x}, t_i ) ]^{2} - [ \bar{g} ( \boldsymbol{x} ) ]^{2} }.\end{equation} Phase averages were defined considering instantaneous realizations associated with the same propeller–rudder configuration, identified here by the angle 
$\varphi$
\begin{equation} \hat{g} ( \boldsymbol{x}, \varphi ) = \frac{1 }{\hat{N}} \sum_{i=1}^{\hat{N}} g [ \boldsymbol{x}, t_0+(i-1) \mathcal{T}_{bp} ],\end{equation}
where 
$\hat {N}$ is the number of the available instantaneous realizations of the flow corresponding to the particular propeller–rudder configuration defined by 
$\varphi$, 
$t_0$ is the time instant when computation of the statistics is started and 
$\mathcal {T}_{bp}$ is the period of the blade passage, 
$\mathcal {T}_{bp}=2{\rm \pi} /\varOmega$, with 
$\varOmega$ standing for the propeller angular speed. Therefore, phase-averaged root mean squares were computed as
\begin{equation} \hat{g}' ( \boldsymbol{x}, \varphi ) = \sqrt{ \frac{1 }{\hat{N}} \sum_{i=1}^{\hat{N}} \{ g [ \boldsymbol{x}, t_0+(i-1) \mathcal{T}_{bp} ] \}^{2} - \{ \hat{g} ( \boldsymbol{x}, \varphi )\}^{2} }. \end{equation}In all results presented below the coordinates are defined in a reference frame placed at the intersection between the propeller plane and its axis.
4.1. Open-water configuration and baseline differences
The results of LES computations of the INSEAN E1658 propeller in the open-water configuration at the same loading conditions will be frequently utilized in the discussion below to better illuminate the complex flow patterns induced by the presence of the rudder. Those computations were reported in the works by Posa et al. (Reference Posa, Broglia, Felli, Falchi and Balaras2018, Reference Posa, Broglia and Balaras2019a,Reference Posa, Broglia, Felli, Falchi and Balarasb). It is worth noting that the numerical resolution in the area around the propeller in the present LES is identical to that adopted for the open-water computations. The grid was significantly refined around the rudder, however, to capture the thin boundary layers forming on its surface. The open-water computations were designed to match the conditions in the experiments by Felli & Falchi (Reference Felli and Falchi2018). Direct comparison of integral quantities such as thrust, torque and efficiency were within 
$3\,\%$ for all loading conditions. Velocity statistics at various locations across the wake were also in very good agreement with the PIV measurements by Felli & Falchi (Reference Felli and Falchi2018). Details can be found in Posa et al. (Reference Posa, Broglia, Felli, Falchi and Balaras2018, Reference Posa, Broglia and Balaras2019a,Reference Posa, Broglia, Felli, Falchi and Balarasb).
The comparison with the open-water case is presented in terms of global parameters of performance in table 1, where the coefficients of thrust and torque as well as the propeller efficiency are reported. They are defined as
\begin{equation} K_T = \frac{T }{\rho n^{2} D^{4}},\quad K_Q = \frac{Q }{\rho n^{2} D^{5}},\quad \eta = \frac{K_T J }{2 {\rm \pi}K_Q},\end{equation}
where 
$T$ and 
$Q$ are the dimensional values of thrust and torque, respectively, and 
$\rho$ is the density of the fluid. All other quantities were defined above. In table 1 also the variation of each quantity from the computations in the open-water configuration to those with downstream hydrofoil is provided. It is shown that, although no dramatic change of thrust and torque occurs, due to the presence of the hydrofoil, the blockage it generates downstream of the propeller alters the pressure field on the pressure side of its blades, producing an increase of both 
$K_T$ and 
$K_Q$. However, the increase of thrust is more significant than that of torque, leading to an overall rise also for the efficiency of propulsion. This is rather uniform across values of advance coefficient, albeit slightly higher at lower propeller loads.
Integral parameters of propeller performance: thrust coefficient, 
$K_T$, torque coefficient, 
$K_Q$, and efficiency, 
$\eta$.
Also the drag and lift coefficients on the hydrofoil were computed, across a spanwise extent defined by 
$-1.0 < y/D < 1.0$, as
\begin{equation} \mathcal{C_D} = \frac{\mathcal{D} }{0.5 \rho U_\infty^{2} c l},\quad \mathcal{C_L} = \frac{\mathcal{L} }{0.5 \rho U_\infty^{2} c l},\end{equation}
where 
$\mathcal {D}$ and 
$\mathcal {L}$ stand for the drag and lift forces, while 
$l=2D$ is the considered spanwise length. Results are reported in table 2. Of course, the average of 
$\mathcal {C_L}$ is not provided, since it is equal to 0. Quite surprisingly the time average of the drag coefficient, 
$\overline {\mathcal {C_D}}$, was found to be a decreasing function of the load on the propeller. This can be explained considering the suction generated by the propeller itself, due to the acceleration of the flow, decreasing the pressure component of the drag acting on the hydrofoil. As a consequence, higher propeller loads correlate with lower values of pressure in front of the hydrofoil, thus with lower values of 
$\mathcal {D}$. Also the standard deviation of the drag coefficient, 
$\widehat {\mathcal {C_D}}$, is provided in table 2. In contrast with the mean value, 
$\widehat {\mathcal {C_D}}$ is an increasing function of the propeller load, likely as a result of the higher levels of turbulence of the propeller wake for lower advance coefficients. For the standard deviation of 
$\mathcal {C_L}$, which is an order of magnitude higher than that of 
$\mathcal {C_D}$, the dependence on the propeller load conditions is not clear. This could be the consequence of both a weaker dependence of 
$\widehat {\mathcal {C_L}}$ on the propeller load and the lift coefficient being statistically equal to 0, resulting in a stronger dependence of 
$\widehat {\mathcal {C_L}}$ on the size of the statistical sample.
Time average of the drag coefficient, 
$\overline {\mathcal {C_D}}$, and standard deviations of the drag and lift coefficients, 
$\widehat {\mathcal {C_D}}$ and 
$\widehat {\mathcal {C_L}}$, acting on the hydrofoil placed downstream of the propeller.
A general overview of the phase-averaged flow is shown in figure 4, where the coherent structures populating the wake are visualized using isosurfaces of the second invariant of the velocity gradient tensor, 
$Q$ (Hunt, Wray & Moin Reference Hunt, Wray and Moin1988). On the left and right sides of figure 4 results for the open-water configuration and those with the downstream hydrofoil are shown, respectively. Visualizations refer to the case at 
$J=0.65$. Views in the top and bottom panels of figure 4 deal with the two opposite sides of the hydrofoil. The computations capture both the tip and hub vortices, as well as smaller streamlined vortices populating the wake of each propeller blade. Although we plan to conduct in the future more detailed computations on the source of these structures, in our earlier work (Posa et al. Reference Posa, Broglia, Felli, Falchi and Balaras2019b) it is pointed out that they do not originate from roll-up within the vorticity sheet released by each blade, nor from separation phenomena experienced by the boundary layer. They appear instead associated with instability phenomena occurring across the blades, since their presence was observed already at their trailing edge, within the boundary layer over their suction side, which keeps attached across the whole extent of the blades. They are instead missing within the boundary layer on the pressure side. The dynamics of the vortices shed by the propeller is substantially affected by the presence of the rudder. After the initial contraction of the radius of the helical trajectory of the tip vortices, due to the flow acceleration produced by the propeller, in the vicinity of the leading edge of the rudder they shift towards the outer radii. Across the rudder they undergo a spanwise displacement, depending on their position across the hydrofoil. It is worth noting that the swirl of the propeller wake, having an azimuthal velocity component, produces a significant spanwise gradient of pressure on both faces of the rudder. The pressure sides (
${RPS}$) are those where the azimuthal velocity is directed towards the surface of the hydrofoil, whereas conditions of lower pressure (
${RSS}$) are generated where the azimuthal velocity within the propeller wake is directed away from the surface of the hydrofoil. Therefore, anti-symmetric conditions are generated on the two opposite surfaces of the rudder, as we are going to demonstrate in detail in the next sections. As shown in figure 4, the behaviour of the pressure and suction side branches of the tip vortices across the rudder is different. They move towards outer and inner spanwise coordinates, respectively. These phenomena are consistent with the experiments by Felli & Falchi (Reference Felli and Falchi2011) and the DES computations by Muscari et al. (Reference Muscari, Dubbioso and Di Mascio2017). Downstream of the rudder the right panels of figure 4 display the signature of the pressure side branches of the tip vortices coming from the two faces of the hydrofoil, whereas the suction side branches undergo an earlier instability, due to interaction with the several structures populating the core of the propeller wake, as discussed more in detail later. At the leading edge of the rudder the hub vortex splits between its two sides, experiencing a spanwise displacement as well, although oriented in the direction opposite to that of the tip vortices. Due to the level of 
$Q$ selected for the visualizations in figure 4, the two branches of the hub vortex are not visible downstream of the rudder. It should be pointed out that in figure 4 the isosurfaces on the left side are less noisy, since in open-water conditions it was possible to exploit the symmetry of the geometry for increasing the size of the statistical sample.
Visualization of the coherent structures in the wake of the propeller for 
$J=0.65$. Phase-averaged isosurfaces based on the Q-criterion, with 
$QD^{2}/U_\infty ^{2}=100$, coloured using the azimuthal vorticity. Left and right panels refer to the conditions without and with downstream hydrofoil, respectively. Top and bottom panels for the two sides of the hydrofoil. The white solid arrows indicating the direction of the propeller rotation. Vorticity scaled by 
$U_\infty /D$.
More details are provided in figure 5, where the top and bottom panels show contours of phase-averaged streamwise velocity and azimuthal vorticity, respectively. On the left, central and right panels the PIV experiments by Felli & Falchi (Reference Felli and Falchi2018), the LES computations in open water-conditions by Posa et al. (Reference Posa, Broglia, Felli, Falchi and Balaras2019b) and the present LES computations with downstream rudder are reported, respectively. The orientation of the meridian semi-plane considered in figure 5 is shown in figure 6, where its projection is visualized via a dashed arrow. The solid arrow 
${A}$ in figure 6 is utilized to indicate the orientation of the viewer in figure 5. In comparison with the open-water conditions the presence of the downstream hydrofoil causes an acceleration of the flow and a slight increase of the diameter of the wake, which are due to the blockage generated at the wake axis. The contours of azimuthal vorticity in the bottom part of figure 5 isolate the coherent structures populating the propeller wake and the signature of the wake of each blade. It is possible to see that the acceleration of the flow at the streamwise coordinates of the rudder produces a stronger bending of the vorticity sheets shed from the trailing edge of each propeller blade. The origin of all the above differences and how they are reflected in the statistics of the flow will be discussed in detail in the following sections. In the bottom part of figure 5 the signature of the tip vortices in the central and right panels is very similar, besides differences tied to the size of the statistical sample. We exploited the symmetry of the open-water computations to increase the size of the sample for the results shown in the central panels. However, downstream of the propeller a better resolution was utilized for the present computations with downstream hydrofoil. The computational grids start differing at 
$z/D=0.15$, with that of the open-water case stretching at a faster rate away from the propeller. For instance, at the location of the leading edge of the hydrofoil, which is 
$z/D=0.5$, the grid spacing across the streamwise direction is 
$40\,\%$ of that utilized for the simulations in open-water conditions. Further downstream the grid of the open-water case keeps coarsening, while uniform grid spacing is adopted across the hydrofoil. Therefore, the comparison of the central and right panels of figure 5 suggests a good grid independence of the signature of the tip vortices. When compared to the results on the left panel (experiments) the footprint of the tip vortices from both computations appears more elongated. Having observed grid independence of this behaviour, this result could be the consequence of the computational grid being more able than the experimental one to resolve the vorticity layer connecting the tip vortices with the wake of each blade, since the spatial resolution of the experiments was lower than that of the simulations. For instance, at the propeller plane the resolution of the computations across the streamwise direction is 16 times higher than in the experiments and across the hydrofoil is still 6 times better. It is also worth noting that the comparison between experiments and computations in open-water conditions shows differences at the wake axis, with values of azimuthal vorticity higher in the former case. As discussed in Posa et al. (Reference Posa, Broglia, Felli, Falchi and Balaras2018, Reference Posa, Broglia, Felli, Falchi and Balaras2019b), this is likely due to inaccuracy of the measurements in that region. It was indeed verified that the ensemble average of the radial velocity from the computations was equal to zero at the wake axis, as expected, due to the symmetry of the average wake system in the configuration with the isolated propeller. This was not the case from the experiments. This discrepancy may be tied to the selection of the parameters of the PIV set-up by Felli & Falchi (Reference Felli and Falchi2018). In that study PIV was mainly tailored to capture the tip vortices at the outer radii, rather than the hub vortex at the inner radii, where flow conditions, for instance in terms of velocity magnitude and direction, are very different.
Phase-averaged contours of streamwise velocity (a–c) and azimuthal vorticity (d–f) from the PIV experiments by Felli & Falchi (Reference Felli and Falchi2018) (a,d), the LES computations by Posa et al. (Reference Posa, Broglia, Felli, Falchi and Balaras2019b) in open-water conditions (b,e) and the present LES computations with downstream hydrofoil (c,f). Contours on the semi-plane shown by the dashed arrow in figure 6. Velocity and vorticity scaled by 
$U_\infty$ and 
$U_\infty /D$, respectively.
Visualization from downstream of the orientation of the meridian semi-plane relative to the contours in figure 5. The arrow 
${A}$ showing the orientation of the viewer in figure 5. Also the positions of the pressure and suction sides of the rudder are indicated as 
${RPS}$ and 
${RSS}$, respectively. For visibility of the whole propeller geometry the rudder is visualized with a reduced opacity.
4.2. Propeller wake upstream of the rudder
In this section results for the propeller wake will be presented upstream of the rudder, with a special focus on the influence of the leading edge on the behaviour of the main coherent structures, especially the hub vortex and the tip vortices.
Figure 7 reports results for the phase-averaged azimuthal component of vorticity for all three simulated cases, in comparison to the corresponding open-water configuration. The contours are shown over the 
$x=0$ plane, which is the one aligned with the downstream rudder. The displacement of the tip vortices towards outer radii is evident and is due to both the blockage effect generated by the hydrofoil and their interaction with the surface of the hydrofoil. This observation is consistent with the findings by Felli & Falchi (Reference Felli and Falchi2011) in their PIV visualizations of the flow and with the computations by Wang et al. (Reference Wang, Guo, Xu and Su2019) and Hu et al. (Reference Hu, Zhang, Sun and Guo2019). Actually, the influence of the downstream rudder extends to all radii, affecting also the bending of the trailing wake from the propeller blades and the hub vortex, which experiences a deflection and an expansion towards outer radii. This phenomenon enhances the interaction of the hub vortex with the vortices populating the inner radii of the propeller wake, coming from the root of its blades (for more details, see the earlier study in open-water conditions by Posa et al. Reference Posa, Broglia, Felli, Falchi and Balaras2019b). We will show that the additional shear generated in this region substantially affects the turbulence levels. In addition, it is interesting to note that the wake of the propeller blades experiences an increase of the levels of vorticity in the vicinity of the leading edge, when compared to the open-water case. The results are similar in a qualitative sense across the three simulated values of the advance coefficient, with higher loads producing larger values of vorticity and more profound local maxima/minima through the span of the wake of each blade.
Phase-averaged contours of azimuthal vorticity on the 
$x=0$ plane for 
$J=0.74$ (a), 
$J=0.65$ (b) and 
$J=0.56$ (c). Top: rudder placed downstream. Bottom: open-water condition. Vorticity scaled by 
$U_\infty /D$. Arrows indicating the location of the tip vortices. Dotted-dashed line for the propeller axis.
The influence of the downstream rudder is also visible in figure 8, where the azimuthal vorticity is shown in the 
$y=0$ plane, orthogonal to the axis of the rudder. The contraction of the propeller wake in open-water conditions, due to the acceleration of the flow through the propeller plane, is overtaken by the blockage effect produced by the rudder, causing wake expansion. Also in this plane the bending of the vorticity sheet shed by each propeller blade is modified, compared to open-water conditions, producing stronger gradients of its pitch across the span: the acceleration of the flow is felt especially at intermediate radii, while the inner radii undergo a deceleration, due to both the potential and viscous interaction with the leading edge of the rudder.
Phase-averaged contours of azimuthal vorticity on the 
$y=0$ plane for 
$J=0.74$ (a), 
$J=0.65$ (b) and 
$J=0.56$ (c). Top: rudder placed downstream. Bottom: open-water condition. Vorticity scaled by 
$U_\infty /D$. Arrows indicating the location of the tip vortices. Dotted-dashed line for the propeller axis.
Figure 9 reports a comparison across advance coefficients on the levels of phase-averaged turbulent kinetic energy on a cross-section just upstream of the leading edge of the hydrofoil (
$z/D=0.49$). At the azimuthal locations aligned with the rudder the signature of the tip vortices wrapping around the leading edge is visible, confirming their displacement towards outer radii, especially at higher loads. At the inner radii, values of the turbulent kinetic energy are significantly higher, compared to both tip vortices and shear layers coming from the trailing edge of the propeller blades. This behaviour affects both the hub vortex at the wake axis and the inner region of the trailing wake from the propeller blades. The blockage produced by the rudder is forcing the large hub vortex towards both right and left sides in each panel of figure 9, likely enhancing also its destabilization. It should be pointed out that the statistics in figure 9 were computed synchronized with the frequency of the blade passage. Therefore, the large values of turbulent kinetic energy observed at the wake axis are associated with fluctuations from the phase-averaged fields. This way it is expected that any periodicity of the motion of the hub vortex, tied to the frequency of the blade passage, is removed. In addition, it was verified that the ensemble-averaged statistics are actually not substantially different from the phase-averaged ones, demonstrating that most turbulence comes from small scales fluctuations. This point is very clear, for instance, in the following figure 15, where both phase-averaged turbulence and ensemble-averaged turbulence are visualized just downstream of the leading edge of the hydrofoil. Significant differences are visible at the outer radii, due to the azimuthal displacement of the tip vortices as the propeller blades change their position in time, while this is not the case at the inner radii, dominated by the two branches of the hub vortex. A similar result was also shown above via the power spectral densities in front of the leading edge of the hydrofoil (figure 3a,b), where the energy at large frequencies, corresponding to ‘actual turbulence’, was verified much higher at the axis than at the outer radii of the propeller wake. The shear produced with the leading edge and with the structures coming from the root of the propeller blades (see Posa et al. Reference Posa, Broglia, Felli, Falchi and Balaras2018, Reference Posa, Broglia, Felli, Falchi and Balaras2019b for more details) is the source of large values of turbulent kinetic energy. It is also worth noting that the hub vortex is splitting into two structures. We will show below that they move over the surface of the hydrofoil on its two sides, without reconnecting downstream. Although values of turbulent kinetic energy at inner radii are dependent on load conditions, the flow physics of the hub vortex was found qualitatively similar across simulated advance coefficients.
Phase-averaged contours of turbulent kinetic energy on the cross-section at 
$z/D=0.49$ for 
$J=0.74$ (a), 
$J=0.65$ (b) and 
$J=0.56$ (c) with rudder placed downstream of the propeller. Turbulent kinetic energy scaled by 
$U_\infty ^{2}$. Solid arrows indicating the signature of the tip vortices wrapping around the leading edge of the hydrofoil. View from downstream. Rudder removed from the visualizations for visibility purposes. The dotted-dashed line standing for the projection of the trajectory swept by the tip of the propeller blades.
Figure 10 shows radial profiles of ensemble-averaged quantities at the same streamwise location as in figure 9 along the 
$y$ direction aligned with the leading edge of the hydrofoil (indicated as 
${\alpha }$ in figure 9). Comparisons against results from the open-water computations are included. For all advance coefficients the streamwise velocity component in figure 10(a) displays a sharp decrease, compared to the open-water configuration, with also small negative values at the axis, due to the collision of the large hub vortex with the leading edge of the rudder. The presence of local maxima/minima at intermediate radii, due to the several vortices populating the wake of each blade, is still distinguishable upstream of the hydrofoil. It is interesting to note the behaviour of the radial velocity in figure 10(b). It demonstrates that for the propeller in isolated conditions the wake contraction is still ongoing at this streamwise location, as the radial velocity is negative across all radii. The values are substantially higher (in magnitude) in the vicinity of the hydrofoil and are strongly negative at inner radii. We will demonstrate later that this result is tied to a spanwise deflection of the hub vortex, which becomes very evident across the rudder. At the outer radii the radial velocity is positive, with a sharp maximum at the edge of the propeller wake, which is consistent with the radial displacement of the tip vortices seen in figure 7. The azimuthal velocity component (figure 10c) experiences a substantial increase, due to the blockage produced by the leading edge of the rudder and the consequent lateral acceleration it generates. The only exception is found at the wake axis, where the flow physics of the hub vortex is substantially modified by the presence of the hydrofoil. Again, at the intermediate radii the presence of local maxima/minima is still visible. In figure 10(d) the comparison of the profiles of turbulent kinetic energy between configurations shows: (i) an increase of the area of large values of turbulent kinetic energy at inner radii, due to the expansion of the hub vortex upstream of the hydrofoil; (ii) local maxima at intermediate radii, tied to the small vortices shed across the span of the propeller blades; (iii) an outward displacement of the peak at the edge of the propeller wake, associated with the tip vortices, displaying also a wider radial extent, due to the overlapping effects of multiple tip vortices. In addition, a slight increase of turbulence within the core of the tip vortices occurs, due to the shear with the boundary layer on the leading edge of the rudder. The dependence on load conditions in all panels of figure 10 is evident, with higher values of all variables for lower advance coefficients, although trends are similar across values of 
$J$.
Ensemble-averaged profiles of streamwise velocity (a), radial velocity (b), azimuthal velocity (c) and turbulent kinetic energy (d) at 
$x/D=0$ and 
$z/D=0.49$, along the direction aligned with the leading edge of the hydrofoil (
${\alpha }$ in figure 9). For comparison, also results in open-water conditions are reported. Legends referring to values of advance coefficient (
$J=0.74, 0.65, 0.56$) and propeller configurations (
$i$: isolated propeller; 
$r$: rudder placed downstream).
Similar radial profiles are provided in figure 11 along a direction orthogonal to the hydrofoil, with 
$y/D=0$ and 
$z/D=0.49$ (marked as line 
${\beta }$ in figure 9). Figure 11(a) shows a slight streamwise acceleration of the flow within the wake core, but the blockage generated by the leading edge of the rudder is mainly balanced by a radial expansion of the propeller wake, compared to the case with the isolated propeller. This expansion is confirmed by the behaviour of the radial velocity in figure 11(b), where the peak at inner radii is due to the lateral deflection produced in the vicinity of the leading edge of the rudder. In contrast with the result seen in figure 10(c), the azimuthal velocity away from the rudder is decreased, compared to all cases in open-water conditions (see figure 11c), especially at the inner radii, which are those closer to the leading edge of the rudder. The profiles for the ensemble-averaged turbulent kinetic energy in figure 11(d) are similar to those in figure 10(d), with the lateral deflection of the hub vortex producing even wider areas of large turbulent kinetic energy at inner radii. However, at the edge of the wake the peaks tied to the tip vortices, although displaced outward also in this case, are slightly lower, compared to those in the open-water configuration, maybe because of the acceleration experienced by the flow. This behaviour differs from that observed at outer radii directly in front of the leading edge of the rudder in figure 10(d).
Ensemble-averaged profiles of streamwise velocity (a), radial velocity (b), azimuthal velocity (c) and turbulent kinetic energy (d) at 
$y/D=0$ and 
$z/D=0.49$, along a direction orthogonal to the downstream hydrofoil (
${\beta }$ in figure 9). For comparison, also results in open-water conditions are reported. Legends referring to values of advance coefficient (
$J=0.74, 0.65, 0.56$) and propeller configurations (
$i$: isolated propeller; 
$r$: rudder placed downstream).
4.3. Propeller wake across the rudder
Figure 12 shows the coherent structures populating the wake of the propeller, using the Q-criterion computed from the phase-averaged flow fields. Only the case at 
$J=0.65$ is considered, since results are qualitatively similar across propeller load conditions. Surfaces are coloured using the streamwise vorticity component. Several structures populate the propeller wake: the tip vortices at the outer radii, the hub vortex at the wake axis and additional helical vortices shed across the span of each blade. They were clearly identified for the same propeller and load conditions in our earlier works dealing with the open-water configuration (Posa et al. Reference Posa, Broglia, Felli, Falchi and Balaras2018, Reference Posa, Broglia, Felli, Falchi and Balaras2019b). In figure 12(a) the tip vortices moving along the suction (lower) side of the rudder experience a gradual displacement towards the wake axis. The tip vortices located on the pressure (upper) side of the rudder surface display the opposite behaviour. Actually this spanwise displacement is faster on the pressure side than on the suction side. In addition, on both upper and lower areas of the hydrofoil a significant deformation of the tip vortices occurs at the leading edge, where the two branches on the two opposite sides of the rudder keep connected. All these results are consistent with earlier experimental and numerical studies in the literature for similar configurations (see, for instance, Felli & Falchi Reference Felli and Falchi2011; Muscari et al. Reference Muscari, Dubbioso and Di Mascio2017; Hu et al. Reference Hu, Zhang, Sun and Guo2019; Wang et al. Reference Wang, Guo, Xu and Su2019), where the spanwise displacement of the tip vortices was mainly explained via the image vortex model. Moving across the chord of the hydrofoil the tip vortices tend to lose their coherence, compared to the upstream locations. Nonetheless, even downstream of the trailing edge of the hydrofoil the weakened signature of the vortices shed from the tip of the propeller blades is still distinguishable. In figure 12(b) the opposite side of the rudder surface is visualized. It is evident that the topology of the tip vortices is anti-symmetric, compared to that in figure 12(a), since the suction and pressure sides are respectively located in the upper and lower areas of the surface of the hydrofoil, relative to the axis of the propeller wake.
Visualization of the coherent structures in the wake of the propeller for 
$J=0.65$. Phase-averaged isosurfaces based on the Q-criterion, with 
$QD^{2}/U_\infty ^{2}=50$, coloured using the streamwise vorticity component. In panels (a) and (b) the two sides of the rudder are shown. Vorticity scaled by 
$U_\infty /D$.
To better illuminate the dynamics of the hub vortex the visualizations in figure 13 were extracted from those in figure 12, removing the outer region of the propeller wake. As reported above, at the leading edge of the rudder it splits into two structures, whose spanwise displacement across the hydrofoil is opposite to that observed for the tip vortices. This is actually consistent with the image vortex model. The hub vortex is counter-rotating, compared to the tip vortices, as demonstrated by the sign of the streamwise vorticity. However, at downstream locations the spanwise velocity of the two branches of the hub vortex switches its orientation, which is very evident downstream of the hydrofoil. The source of this counter-deflection of the two branches of the hub vortex on both sides of the rudder will become clearer in the following discussion. It is worth noting that downstream of the rudder the two branches of the hub vortex keep separated, moving towards opposite spanwise directions. Figure 13 also displays several streamlined vortices populating the wake of each blade. As in the open-water case, their orientation is roughly consistent with the geometrical pitch of the propeller blades, which means that their angle relative to the azimuthal direction is increasing from outer to inner radii. These vortices feature negative values of streamwise vorticity, thus they are counter-rotating, relative to the tip vortices, shown in figure 12.
Visualization of the coherent structures in the wake of the propeller for 
$J=0.65$. Phase-averaged isosurfaces based on the Q-criterion, with 
$QD^{2}/U_\infty ^{2}=50$, coloured using the streamwise vorticity component. In panels (a) and (b) the two sides of the rudder are shown. The outer region of the propeller wake was removed for better clarity of its inner structures. Vorticity scaled by 
$U_\infty /D$.
The coherent structures shown in figure 13 were also visualized in figure 14 for a lower value of 
$Q$, which makes the core of the inner vortices less distinguishable, but better highlights the interaction between the outer and inner structures at downstream locations. Figure 14(a) demonstrates that, while the vortices in the upper region of the inner wake experience almost no interaction with other structures, besides that with the rudder boundary layer, in the lower region they are pushed upward by the tip vortices. The positive values of streamwise vorticity in the lower part of figure 14(a) come indeed from the tip vortices moving upward and affecting all other structures of the propeller wake. This is also the likely source of the counter-deflection of the hub vortex in the downstream region of the rudder, observed in figure 13. It is clear that the above interaction triggers the destabilization of the coherent structures of the propeller wake and an increase of the levels of turbulence. Figure 14(b) shows similar phenomena, but in the upper region of the propeller wake, due to the anti-symmetry of the overall wake system. In both panels of figure 14 it is also interesting to see that the vortices populating the intermediate radii do not experience the same spanwise displacement as the tip vortices and the two branches of the hub vortex, except for that due to the interaction with the tip vortices moving towards the wake axis. This behaviour may be due to the lower intensity of the former, compared to the latter, making their displacement much smaller.
Visualization of the coherent structures in the wake of the propeller for 
$J=0.65$. Phase-averaged isosurfaces based on the Q-criterion, with 
$QD^{2}/U_\infty ^{2}=25$, coloured using the streamwise vorticity component. In panels (a) and (b) the two sides of the rudder are shown. The outer region of the propeller wake was removed for better clarity of its inner structures. Vorticity scaled by 
$U_\infty /D$.
The signature of the deformation experienced by the tip vortices, interacting with the leading edge of the rudder, is evident in figure 15, via both phase-averaged (top) and ensemble-averaged (bottom) contours of turbulent kinetic energy at the cross-section of streamwise coordinate 
$z/D=0.51$, just downstream of the leading edge. In the top part of figure 15 the arrows point to the legs of the tip vortices, wrapping around the leading edge of the rudder. They were already distinguishable in figure 9 just upstream of the hydrofoil, but become even more obvious here, since the stretching of these vortices along the streamwise direction and their interaction with the boundary layer produce higher values of turbulent kinetic energy within their cores. It is also evident in figure 15 that the stretching of the tip vortices around the leading edge of the hydrofoil is coupled with their outward spanwise displacement along the leading edge of the hydrofoil and increases as they move downstream (see figure 12). For clarity, in each panel of figure 15 the wake of each blade, indicated via a lower case letter, was correlated with the signature of the corresponding tip vortex, indicated via an upper case letter. We should recall that the pitch of the tip vortices depends on the advance coefficient. Their delay relative to the wake of each blade increases at higher loads. This is shown via the comparison across the top panels of figure 15. It is also interesting to note that the wake 
${g}$ is experiencing in all cases a deformation due to the presence of the hydrofoil, increasing its azimuthal distance from the wake 
${a}$ and reducing that from the wake 
${f}$. Such influence on the topology of the propeller wake will become more dramatic as it develops downstream across the hydrofoil. Ensemble averages in the bottom part of figure 15 confirm the local increase of turbulent kinetic energy in the areas of interaction between the leading edge of the rudder and the tip vortices. This is especially clear for 
$J=0.65$, because of the choice of the colour scale, selected for comparison purposes across the three simulated load conditions.
Contours of turbulent kinetic energy on the cross-section at 
$z/D=0.51$ for 
$J=0.74$ (a), 
$J=0.65$ (b) and 
$J=0.56$ (c). Top: phase-averaged statistics. Bottom: ensemble-averaged statistics. Turbulent kinetic energy scaled by 
$U_\infty ^{2}$. Grey area showing the cross-section of the hydrofoil. The dotted-dashed line standing for the projection of the trajectory swept by the tip of the propeller blades.
The modified topology of the propeller wake across the rudder is evident in figure 16, where both phase-averaged and ensemble-averaged contours of turbulent kinetic energy are reported over the cross-section located at 
$z/D=0.8$. As for the open-water condition, moving away from the propeller plane the footprint of both tip and hub vortices is evident, whereas that of the smaller vortices shed across the span of the propeller wake declines at a much faster rate (see also the discussion in Posa et al. Reference Posa, Broglia, Felli, Falchi and Balaras2018, Reference Posa, Broglia, Felli, Falchi and Balaras2019b). Both suction sides of the rudder experience a contraction of the propeller wake, with the tip vortices and the hub vortex moving closer to each other at a faster rate for higher propeller loads. The opposite occurs on both pressure sides of the rudder. Note that in the upper and lower parts of each panel of figure 16 the suction side is on the right and the left, respectively.
Contours of turbulent kinetic energy on the cross-section at 
$z/D=0.80$ for 
$J=0.74$ (a), 
$J=0.65$ (b) and 
$J=0.56$ (c). Top: phase-averaged statistics. Bottom: ensemble-averaged statistics. Turbulent kinetic energy scaled by 
$U_\infty ^{2}$. Grey area showing the cross-section of the hydrofoil. The dotted-dashed line standing for the projection of the trajectory swept by the tip of the propeller blades.
Radial profiles were extracted from the ensemble-averaged statistics in the bottom part of figure 16 at 
$45^{\circ }$ (figure 17a) and 
$-45^{\circ }$ (figure 17b), relative to the 
$xy$ reference frame, for a more quantitative comparison across cases. Those directions are indicated by 
${\gamma }$ and 
${\delta }$ in the bottom panels of figure 16. The strong dependence on load conditions is obvious. In figure 17(a,b) the outer maxima, tied to the tip vortices, are very similar, although their radial location is different and in the former case the peak value is slightly higher. This trend is confirmed for the local maxima/minima at intermediate radii, that are the ensemble-averaged signature of the trailing wake from the propeller blades. At inner radii differences between the two panels of figure 17 are more interesting. In figure 17(a) the inner peak is higher and broader than in figure 17(b), due to the branch of the hub vortex, which is displaced upwards on this side of the hydrofoil. In figure 17(b) the inner peak is produced within the rudder boundary layer and is much sharper than that seen in figure 17(a).
Radial profiles of ensemble-averaged turbulent kinetic energy extracted from the bottom panels of figure 16 across the directions 
${\gamma }$ (a) and 
${\delta }$ (b) at 
$45^{\circ }$ and 
$-45^{\circ }$, relative to the 
$xy$ reference frame.
Further downstream the deformation of the propeller wake leads to higher levels of turbulent kinetic energy, because of the increased shear between the structures populating the flow. In figure 18 the location at 
$z/D=1.2$ is shown, in the vicinity of the trailing edge, whose streamwise coordinate is 
$z/D=1.33$. At this location the tip vortices, the trailing wake from the propeller blades and each branch of the hub vortex come very close to each other on both suction sides of the rudder. On the pressure sides of the rudder turbulence is significantly lower than on the suction sides, due to the lower levels of shear stress. Nonetheless additional areas of high turbulent kinetic energy are visible in the outer region of both rudder pressure sides in the vicinity of the surface of the hydrofoil. There the tip vortices are stretched in such a way to be almost aligned with the surface of the body, as reported in the visualizations of figure 12, increasing this way the region of shear with its boundary layer. In contrast, on the suction sides of the rudder the most significant increase of turbulence levels occurs within the shrinking area between them and the hub vortex.
Contours of turbulent kinetic energy on the cross-section at 
$z/D=1.20$ for 
$J=0.74$ (a), 
$J=0.65$ (b) and 
$J=0.56$ (c). Top: phase-averaged statistics. Bottom: ensemble-averaged statistics. Turbulent kinetic energy scaled by 
$U_\infty ^{2}$. Grey area showing the cross-section of the hydrofoil. The dotted-dashed line standing for the projection of the trajectory swept by the tip of the propeller blades.
The above discussion is reinforced by figure 19, where the phase-averaged and ensemble-averaged turbulent shear stresses, 
$\widehat {u_r'u_z'}$ and 
$\overline {u_r'u_z'}$, are shown at 
$z/D=1.2$ for the three values of advance coefficient. Their correlation with the contours of turbulent kinetic energy in figure 18 is clear. While the shear on both pressure sides of the rudder, where the tip vortices move outwards, is relatively small, it is obviously increased on the suction sides, where the propeller wake experiences a contraction. However, as reported above, on both pressure sides the stretching of the tip vortices along the surface of the rudder enhances the interaction with the boundary layer, producing large values of shear stress and turbulent fluctuations at the most outward locations. In addition, in all panels of figure 19 we can see that the turbulent shear stress is significant at the outer radii across the whole azimuthal extent of the propeller wake, due to the mutual interaction of the tip vortices, which is consistent with the results in open-water conditions reported in our earlier studies.
Contours of turbulent shear stress, 
$\widehat {u_r'u_z'}$ and 
$\overline {u_r'u_z'}$, on the cross-section at 
$z/D=1.20$ for 
$J=0.74$ (a), 
$J=0.65$ (b) and 
$J=0.56$ (c). Top: phase-averaged statistics. Bottom: ensemble-averaged statistics. Shear stress scaled by 
$U_\infty ^{2}$. Grey area showing the cross-section of the hydrofoil. The dotted-dashed line standing for the projection of the trajectory swept by the tip of the propeller blades.
Radial profiles were also extracted from the ensemble-averaged statistics reported in figures 18 and 19, across the directions indicated by 
${\gamma }$ and 
${\delta }$ in the bottom panels of figures 18 and 19. In the top part of figure 20 it is clear that the turbulent kinetic energy on the suction side (top panel of figure 20a), when compared to that on the pressure side (top panel of figure 20b), displays: (i) a higher outer peak; (ii) a broad additional maximum at intermediate radii; (iii) a sharper maximum within the rudder boundary layer, where fluctuations are enhanced by the presence of the branch of the hub vortex, displaced towards the suction side of the rudder; (iv) a smaller radial extent of the wake, due to the displacement of the tip vortices across the span of the rudder. This occurs at all simulated values of the advance coefficient, whose effect consists in increasing the levels of turbulent fluctuations. The profiles of ensemble-averaged 
$\overline {u_r'u_z'}$ turbulent shear stress are consistent with those for the turbulent kinetic energy. In the bottom panel of figure 20(a) four maxima/minima can be identified. The one indicated as 
${A}$ is located within the rudder boundary layer, whereas the one in 
${B}$ is attributable to the interaction of the hub vortex with the turbulent boundary layer. In 
${C}$ the increase of the shear is likely due to the interaction of the hub vortex with the shrinking propeller wake on the suction side, producing the peak of turbulent kinetic energy seen in the top panel of figure 20(a) at 
$r/D \approx 0.2$. The maximum in 
${D}$ is located within the region of the tip vortices and is mainly due to their mutual interaction, being rather uniform across the whole azimuthal extent of the propeller wake, as seen in figure 19. The comparison with the lower panel of figure 20(b) demonstrates that this maximum is very similar on the pressure side of the rudder, although located further away from it. However, overall the results shown in figure 20(b) (bottom) are very different from those in figure 20(a) (bottom). Only the peaks in 
${A}$ (rudder boundary layer) and in 
${D}$ (tip vortices) can be identified, due to the lack of the contribution by the hub vortex. It is worth noting that in such comparison across cases and locations within the propeller wake the turbulent shear stress 
$\overline {u_r'u_z'}$ was reported, since it was found to be the most significant one across the streamwise extent of the hydrofoil, especially for characterizing the interactions involving the coherent structures populating the wake flow.
Radial profiles of ensemble-averaged turbulent kinetic energy (top) and turbulent shear stress 
$\overline {u_r'u_z'}$ (bottom) extracted from the bottom panels of both figures 18 and 19 across the directions 
${\gamma }$ (a) and 
${\delta }$ (b) at 
$45^{\circ }$ and 
$-45^{\circ }$, relative to the 
$xy$ reference frame.
4.4. Flow downstream of the rudder
The rudder has a significant impact on the flow downstream. The left side of figure 21 displays contours of phase-averaged turbulent kinetic energy on the 
$x=0$ plane. The expansion of the propeller wake is evident and increases with the propeller load. The maxima seen in 
${A}$ are due to the shear between the branches of the tip vortices coming from the two pressure sides of the rudder surface and the shear layer from the trailing edge of the rudder. Although the tip vortices are not anymore as coherent as at upstream locations, their signature at outer radii is still clear, with maxima increasing at downstream locations in 
${B}$. Away from the rudder an increase of turbulent kinetic energy is observed also at inner radii in 
${C}$. This phenomenon can be better understood looking at the contours in the orthogonal plane 
$y=0$, as reported on the right side of figure 21. There a contraction of the propeller wake occurs downstream of the rudder. The branches of the tip vortices on the suction sides of the rudder move closer to the wake axis and to each branch of the hub vortex, leading downstream to an increase of the values of turbulent kinetic energy (
${D}$), due to the shear generated between these counter-rotating structures. Further downstream the stress generated between the coherent structures coming from the two sides of the hydrofoil, converging towards the wake axis, leads to high levels of turbulent kinetic energy at inner radii (
${E}$). The dependence of the values of turbulent kinetic energy at downstream locations on the advance coefficient is very strong.
Contours of phase-averaged turbulent kinetic energy on the 
$x=0$ (a,c,e) and 
$y=0$ (b,d,f) meridian planes for 
$J=0.74$, 
$J=0.65$ and 
$J=0.56$ from top to bottom. Turbulent kinetic energy scaled by 
$U_\infty ^{2}$.
The above discussion is reinforced by the ensemble-averaged fields of 
$\overline {u_r'u_z'}$ turbulent shear stress, which are reported in figure 22, demonstrating a strong correlation with the fields of turbulent kinetic energy shown in figure 21. On the 
$x=0$ plane (figure 22a,c,e) the wide areas of positive shear stress in 
${A}$ are the most obvious feature at the trailing edge of the rudder and correlate with the maxima of turbulent kinetic energy seen in figure 21. The behaviours in 
${B}$ and 
${C}$ are also consistent with those observed for the normal turbulent stresses, although these phenomena are more obvious at higher load conditions. In 
${B}$ the destabilization of the pressure side branches of the tip vortices is the likely source of the increasing values of shear stress, while in 
${C}$ the major role of the two branches of the hub vortex will be demonstrated later. On the 
$y=0$ plane (figure 22b,d,f) the region of negative values of shear stress downstream of the propeller is associated with the hub vortex. Then it splits and moves towards the outer radii downstream of the rudder, coming closer to the suction side branches of the tip vortices, which experience a simultaneous displacement towards the inner radii. This causes the increase of turbulent kinetic energy and turbulent shear stress observed in 
${D}$. In 
${E}$ the shear between the two branches of the hub vortex and the suction side branches of the tip vortices produces the sharp increase of turbulent kinetic energy, discussed above (this will become clear in the discussion below). It is also interesting to point out in figure 22 that the streamwise locations of the regions indicated by 
${D}$ and 
${E}$ move upstream at higher loads, suggesting an earlier interaction between the coherent structures populating the wake of the overall system at lower values of the advance coefficient.
Contours of ensemble-averaged turbulent shear stress 
$\overline {u_r'u_z'}$ on the 
$x=0$ (a,c,e) and 
$y=0$ (b,d,f) meridian planes for 
$J=0.74$, 
$J=0.65$ and 
$J=0.56$ from top to bottom. Shear stress scaled by 
$U_\infty ^{2}$.
To better illuminate the complex flow phenomena downstream of the overall system, figure 23 shows the ensemble-averaged turbulent kinetic energy at the cross-sections of streamwise coordinates 
$z/D=1.4$, 
$z/D=2.0$ and 
$z/D=3.0$ (locations are marked in figure 21 by the dashed lines 
${\varepsilon }$, 
${\zeta }$ and 
${\eta }$, respectively). Note that the trailing edge of the rudder is located at 
$z/D=1.33$. For brevity only one load condition (
$J=0.65$) is considered, since the physics is similar across advance coefficients. The wake of the hydrofoil, initially aligned with the 
$y$ direction, experiences an azimuthal deflection, because of the swirl coming from the propeller wake. Figure 23 also demonstrates the gradual displacement towards the wake axis of the branches of the tip vortices coming from the suction sides of the rudder. Their interaction with the two branches of the hub vortex leads to an increase of turbulence from 
$z/D=1.4$ to 
$z/D=2.0$. The simultaneous displacement of the branches of the hub vortex upwards (on the left side) and downwards (on the right side) enhances the anti-symmetry of the overall wake and the shear between its several structures. For instance, the two sides of the hub vortex come closer to both the wake of the hydrofoil and the pressure side branches of the tip vortices. As a consequence, a further increase of turbulent kinetic energy occurs from 
$z/D=2.0$ to 
$z/D=3.0$.
Contours of ensemble-averaged turbulent kinetic energy on the cross sections at 
$z/D=1.4$ (a), 
$z/D=2.0$ (b) and 
$z/D=3.0$ (c) for 
$J=0.65$. The three streamwise locations are indicated by the dashed lines 
${\varepsilon }$, 
${\zeta }$ and 
${\eta }$ in the middle row of figure 21. Turbulent kinetic energy scaled by 
$U_\infty ^{2}$. The dotted-dashed line standing for the projection of the trajectory swept by the tip of the propeller blades.
The cross-stream displacement along the streamwise evolution of the structures populating the wake of the overall system is better characterized in figure 24. There the ensemble-averaged streamwise component of vorticity is shown for the case at 
$J=0.65$. In figure 24(a) the cores of the two branches of the hub vortex are distinguishable and their topology relative to the tip vortices is in good qualitative agreement with the experimental observation by Felli & Falchi (Reference Felli and Falchi2011). It is also interesting to note the presence of a thin layer of opposite vorticity between the two branches of the hub vortex, which is mainly due to the shear layer shed from the trailing edge of the rudder. In the case of an infinite foil in uniform flow the latter will dissipate quickly. In contrast, in the present case it experiences the large-scale disturbance from the hub vortices, leading to its roll-up, as visible in figure 24(b), which is similar to the spiral roll-up in the case of a finite wing. In figure 24(c) these regions of positive vorticity eventually merge with the tip vortices further downstream, which dominate the outer region of the overall wake.
Contours of ensemble-averaged streamwise vorticity on the cross sections at 
$z/D=1.4$ (a), 
$z/D=2.0$ (b) and 
$z/D=3.0$ (c) for 
$J=0.65$. The three streamwise locations are indicated by the dashed lines 
${\varepsilon }$, 
${\zeta }$ and 
${\eta }$ in the middle row of figure 21. Vorticity scaled by 
$U_\infty /D$. Note that the colour scale of panel (c) differs from those of panels (a) and (b). The dotted-dashed line standing for the projection of the trajectory swept by the tip of the propeller blades.
5. Conclusions
Large eddy simulations were carried out on a notional propeller (the INSEAN E1658) upstream of a hydrofoil (an infinite NACA0020), mimicking a rudder at zero incidence, relative to the free-stream. Three values of the propeller load were considered, corresponding to advance coefficients equal to 
$J=0.74, 0.65, 0.56$. An immersed-boundary methodology was adopted to enforce no-slip boundary conditions on the surface of the two bodies, immersed within a cylindrical grid composed of 
$\sim$1.9 billion nodes. The overall approach was validated via comparisons against measurements by Felli & Falchi (Reference Felli and Falchi2018) on the same propeller in open-water conditions, using a computational grid with the same resolution as in the present work (Posa et al. Reference Posa, Broglia, Felli, Falchi and Balaras2018, Reference Posa, Broglia, Felli, Falchi and Balaras2019b). The reported results were also in good qualitative agreement with earlier experimental and computational studies in similar configurations (see Felli & Falchi Reference Felli and Falchi2011; Muscari et al. Reference Muscari, Dubbioso and Di Mascio2017; Hu et al. Reference Hu, Zhang, Sun and Guo2019; Wang et al. Reference Wang, Guo, Xu and Su2019).
Compared to earlier studies in the literature, the present computations point to the importance of the hub vortex and its interaction with the boundary layer on the hydrofoil on the dynamics of the overall wake. The hub vortex is split into two branches, with a lateral expansion, enhancing its shear with the vortices coming from the inner radii of the propeller blades and increasing the turbulence levels within the core of the propeller wake. Across the hydrofoil the two branches of the hub vortex experience a spanwise shift opposite to that of the tip vortices, coming closer to their suction side branches. This behaviour promotes shear and triggers higher levels of turbulent fluctuations. In addition, as a consequence of this interaction, the two branches of the hub vortex eventually switch their spanwise velocity. These phenomena persist downstream of the hydrofoil, reinforcing the anti-symmetry of the overall wake and producing a non-monotonic evolution of turbulent kinetic energy as the two branches of the hub vortex approach the pressure side branches of the tip vortices. Although turbulence levels were found to be an increasing function of the load on the propeller, the qualitative dependence of the flow physics on the advance coefficient was rather limited.
Despite the complexity, the baseline dynamics above is fairly robust as our observations are in good qualitative agreement with experiments and computations in the literature involving different geometries and working conditions. In all cases, at the leading edge of the rudder the tip vortices experience an outward displacement, stretching over the two sides of the hydrofoil, as they remain connected downstream. While the suction side branches move towards inner radial coordinates across the hydrofoil, the pressure side branches experience a deflection towards outer radii. As a consequence, the overall wake undergoes an expansion within the plane of the hydrofoil, due to the opposite spanwise displacement of the propeller wake moving across the two opposite sides of the rudder. Ongoing work focuses on investigating the impact of such phenomena on the acoustic signature of the overall propeller/hydrofoil system.
Acknowledgements
A.P. was supported by the European Union Horizon 2020 research and innovation programme HOLISHIP ‘Holistic Optimisation of Ship Design and Operation for Life Cycle’, under grant agreement n. 689074. E.B. was partially supported by the Office of Naval Research Grant N00014-18-1-2671, monitored by Dr K.-H. Kim. The authors are grateful to M. Brown (Naval Surface Warfare Center, Carderock Division, Maryland), for generating the Lagrangian grid of the INSEAN E1658 propeller, and to M. Felli and M. Falchi (CNR-INM), for providing their experimental results. We acknowledge PRACE for awarding us access to Marconi-KNL at CINECA, Italy, with an allocation granted to the project ‘LESRudProp – LES of the Interaction Between a Rudder and the Wake of a Propeller’ (project no. 2018184409) in the framework of the 17th PRACE Call for Proposals for Project Access.
Declaration of interests
The authors report no conflict of interest.
