2.1 Mathematical model

The numerical simulation of the flow field is performed by the integration of the Reynolds-averaged Navier–Stokes equations written in non-dimensional form in the frame of reference fixed to the propeller in terms of absolute velocity:

(2.1) $$\begin{eqnarray}\displaystyle \left\{\begin{array}{@{}l@{}}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0\\ \displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D74E}\times \boldsymbol{u}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }[(\boldsymbol{u}-\unicode[STIX]{x1D74E}\times \boldsymbol{r})\otimes \boldsymbol{u}]+\unicode[STIX]{x1D735}p-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D70F}=0,\end{array}\right. & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{u}$ is the absolute velocity, $\unicode[STIX]{x1D74E}$ is the (constant) angular velocity, $\boldsymbol{r}=\boldsymbol{x}-\boldsymbol{x}_{0}$ is the position of point $\boldsymbol{x}$ with respect to $\boldsymbol{x}_{0}$ (origin of the non-inertial frame of reference), $p$ is the pressure, $\unicode[STIX]{x1D749}=\unicode[STIX]{x1D708}_{T}(\unicode[STIX]{x1D735}\boldsymbol{u}+\unicode[STIX]{x1D735}\boldsymbol{u}^{\text{T}})$ the stress tensor, $\unicode[STIX]{x1D708}_{T}=1/Re+\unicode[STIX]{x1D708}_{t}$ , $\unicode[STIX]{x1D708}_{t}$ being the turbulent viscosity calculated by the detached eddy simulation (DES) approach described in Spalart (Reference Spalart2009) and $Re$ the Reynolds number (later defined).

As boundary condition, at solid walls the velocity of the fluid is set equal to the local wall velocity, i.e. $\boldsymbol{u}=\unicode[STIX]{x1D74E}\times \boldsymbol{r}$ , at the inflow the velocity is set to the value of the undisturbed flow, $\boldsymbol{u}=\boldsymbol{U}_{\infty }$ , whereas at the outflow the velocity is extrapolated from inner points.

2.2 Numerical model

The numerical integration of the governing equations (2.1) is performed by means of an in-house solver. It is thoroughly described in Muscari et al. (Reference Muscari, Di Mascio and Verzicco2013) and the references cited therein, and we will briefly report here its main features for the sake of completeness.

The numerical algorithm is based on a finite volume technique with pressure and velocity co-located at cell centres. The computational domain is partitioned into structured blocks, each one subdivided into hexahedrons. Conservation laws are then applied to each hexahedron:

(2.2) $$\begin{eqnarray}\displaystyle \left\{\begin{array}{@{}l@{}}\displaystyle \mathop{\sum }_{i=1}^{6}\int _{{\mathcal{S}}_{i}}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}S=0\\ \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}\int _{{\mathcal{V}}}\boldsymbol{u}\,\text{d}V+\mathop{\sum }_{i=1}^{6}\int _{{\mathcal{S}}_{i}}(\boldsymbol{{\mathcal{F}}}_{c}-\boldsymbol{{\mathcal{F}}}_{d})\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}S+\int _{{\mathcal{V}}}\unicode[STIX]{x1D74E}\times \boldsymbol{u}\,\text{d}V=0,\end{array}\right. & & \displaystyle\end{eqnarray}$$

where ${\mathcal{S}}_{i}$ is the $i$ th face of the generic finite volume ${\mathcal{V}}$ and the fluxes are decomposed into a convective part $\boldsymbol{{\mathcal{F}}}_{c}$ and a diffusive term $\boldsymbol{{\mathcal{F}}}_{d}$ .

A standard second-order centred scheme is used for the computation of viscous terms, whereas a centred fourth-order scheme has been applied for the computation of Eulerian terms. Because of the treatment of the viscous terms, the scheme is formally second-order accurate in space, although previous numerical experience proved that the use of a high-order approximation of the Eulerian terms can remarkably reduce the actual error.

The physical time derivatives are approximated by a second-order accurate, three-point backward finite difference formula. The resulting scheme is fully implicit in time and a dual time stepping approach is adopted for its integration.

The equations are discretized on a block structured grid with partial overlapping. This approach renders domain discretization and grid quality control much easier than with structured meshes with abutting blocks when dealing with complex geometries like the one under investigation. High performance computing is achieved by a hybrid shared (OMP) and distributed (MPI) memory parallelization.

5.1 Proper orthogonal decomposition

POD is a mathematical procedure that, given an ensemble of data preliminary removed of its mean, constructs a basis for the ensemble that is optimal, in the sense that no other decomposition of the same order captures an equivalent amount of kinetic energy. This technique uses a second-order statistics, namely the two-point correlation tensor, in order to extract organized large-scale structures from turbulent flows. The most energetic modes can be viewed as ‘dominant’ flow structures.

For our case, the data set is formed by the values of velocity and pressure in the $M=171\times 32\times 64$ nodes of the cylindrical block described above, for $N=1200$ instants, arranged to form a snapshot matrix:

(5.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D64E}=[\boldsymbol{u}(\boldsymbol{x},t_{1}),\boldsymbol{u}(\boldsymbol{x},t_{2}),\ldots ,\boldsymbol{u}(\boldsymbol{x},t_{N})]. & & \displaystyle\end{eqnarray}$$

The POD modes are obtained by searching the $N$ functions $\unicode[STIX]{x1D753}_{i}$ that yield the best representation of the data set, that is that maximize the inner products

(5.2) $$\begin{eqnarray}\displaystyle I=\frac{\langle \boldsymbol{u}(\boldsymbol{x},t_{i}),\unicode[STIX]{x1D753}_{i}(\boldsymbol{x})\rangle ^{2}}{\Vert \unicode[STIX]{x1D753}_{i}(\boldsymbol{x})\Vert ^{2}}. & & \displaystyle\end{eqnarray}$$

The resolution of the maximization problem is performed using the calculus of variations and the result is an integral equation:

(5.3) $$\begin{eqnarray}\displaystyle \int _{\unicode[STIX]{x1D6FA}}\unicode[STIX]{x1D64D}(\boldsymbol{x},\boldsymbol{x}^{\prime })\unicode[STIX]{x1D753}(\boldsymbol{x}^{\prime })\,\text{d}\boldsymbol{x}^{\prime }=\unicode[STIX]{x1D706}\unicode[STIX]{x1D753}(\boldsymbol{x}). & & \displaystyle\end{eqnarray}$$

The kernel of this equation is the two-point correlation tensor $\unicode[STIX]{x1D64D}(\boldsymbol{x},\boldsymbol{x}^{\prime })$ that is

(5.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D64D}(\boldsymbol{x},\boldsymbol{x}^{\prime })=\int _{0}^{\text{T}}\boldsymbol{u}(\boldsymbol{x},t)\boldsymbol{u}(\boldsymbol{x}^{\prime },t)\,\text{d}t. & & \displaystyle\end{eqnarray}$$

Since the autocorrelation tensor is self-adjoint and positive semidefinite (Holmes et al. Reference Holmes, Lumley, Berkooz and Rowley2012), the POD modes are orthogonal and satisfy the following orthonormality condition

(5.5) $$\begin{eqnarray}\displaystyle \int _{\unicode[STIX]{x1D6FA}}\unicode[STIX]{x1D753}_{n}(\boldsymbol{x})\unicode[STIX]{x1D753}_{m}(\boldsymbol{x})\,\text{d}\boldsymbol{x}=\unicode[STIX]{x1D6FF}_{nm}. & & \displaystyle\end{eqnarray}$$

The eigenvalues $\unicode[STIX]{x1D706}_{i}$ provide the energy contribution associated with the corresponding $i$ th POD mode $\unicode[STIX]{x1D753}_{i}$ ; the energy percentage content of the $i$ th mode and the cumulative energy up to the $i$ th mode are defined as

(5.6a,b ) $$\begin{eqnarray}\displaystyle E_{i}=\frac{\unicode[STIX]{x1D706}_{i}}{\mathop{\sum }_{j=1}^{N}\unicode[STIX]{x1D706}_{j}},\quad E_{i}^{cum}=\mathop{\sum }_{k=1}^{i}E_{k}. & & \displaystyle\end{eqnarray}$$

One of the major benefits of POD comes from its ability to sort the obtained modal basis in relation to their relative importance; in this way the dominant features may be selected to obtain a reduced-order representation of the full problem.

The temporal coefficients $a_{i}(t)$ are obtained via a projection of the instantaneous signal onto the POD eigenfunctions:

(5.7) $$\begin{eqnarray}\displaystyle a_{i}(t)=\langle \boldsymbol{u}(\boldsymbol{x},t),\unicode[STIX]{x1D753}_{i}(\boldsymbol{x})\rangle . & & \displaystyle\end{eqnarray}$$

The size of the correlation tensor is based on the spatial degrees of freedom $N_{DOF}$ of the data. For numerical data, $N_{DOF}$ is generally large and is equal to the number of grid points $M$ times the number of variables $n_{c}$ in the array $\boldsymbol{u}$ used for the decomposition, i.e. $N_{DOF}=M\times n_{c}$ . When $N_{DOF}$ is significantly higher than the number of ensemble members $N$ , in order to reduce the computational cost, Sirovich (Reference Sirovich1987) proposes the method of snapshots which requires the solution of an $N\times N$ eigenvalue problem instead of the original eigenvalue problem of size $N_{DOF}\times N_{DOF}$ . This can be accomplished by analysing the correlation of instantaneous snapshots of the flow, defined as

(5.8) $$\begin{eqnarray}\displaystyle C(t,t^{\prime })=\frac{1}{N}\int _{\unicode[STIX]{x1D6FA}}\boldsymbol{u}(\boldsymbol{x},t)\boldsymbol{u}(\boldsymbol{x},t^{\prime })\,\text{d}\boldsymbol{x}. & & \displaystyle\end{eqnarray}$$

$C(t,t^{\prime })$ is then used as a kernel to formulate a eigenvalue problem similar to (5.3). The solution of the eigenvalue problem will result in the temporal coefficients $a_{i}(t)$ , which are the eigenvectors, and the eigenvalues $\unicode[STIX]{x1D706}_{i}$ . The POD basis is then calculated by post-multiplying the snapshot matrix with the obtained eigenvectors. The reconstructed field can be obtained as:

(5.9) $$\begin{eqnarray}\displaystyle \boldsymbol{u}(\boldsymbol{x},t)=\mathop{\sum }_{n=1}^{N}a_{i}(t)\unicode[STIX]{x1D753}_{i}(\boldsymbol{x}). & & \displaystyle\end{eqnarray}$$

5.2 Dynamic mode decomposition

Dynamic mode decomposition (DMD) is a relatively recent development in the field of modal decomposition (Schmid Reference Schmid2010; Rowley & Dawson Reference Rowley and Dawson2017).

Given a sequence of snapshots $\boldsymbol{v}_{i}$ , DMD assumes the existence of a constant linear mapping $\unicode[STIX]{x1D63C}$ that connects the flow field $\boldsymbol{v}_{i}$ to the next $\boldsymbol{v}_{i+1}$ ,

(5.10) $$\begin{eqnarray}\displaystyle \boldsymbol{v}_{i+1}=\unicode[STIX]{x1D63C}\boldsymbol{v}_{i}, & & \displaystyle\end{eqnarray}$$

such that one may write the sequence of the flow fields, the snapshots matrix $\unicode[STIX]{x1D651}_{1}^{N}$ , as a Krylov sequence,

(5.11) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D651}_{1}^{N}=\{\boldsymbol{v}_{1},\unicode[STIX]{x1D63C}\boldsymbol{v}_{1},\ldots ,\unicode[STIX]{x1D63C}^{N-1}\boldsymbol{v}_{1}\}. & & \displaystyle\end{eqnarray}$$

The goal of DMD is to determine the characteristics of the dynamical process described by $\unicode[STIX]{x1D63C}$ without the explicit knowledge of $\unicode[STIX]{x1D63C}$ . DMD can be interpreted as an extension of the classical Arnoldi algorithm, where the eigenvalues and eigenvectors of $\unicode[STIX]{x1D63C}$ can be evaluated through the identification of the eigenvalues and the eigenvectors of a companion matrix $\unicode[STIX]{x1D64E}$

(5.12) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D64E}=\left[\begin{array}{@{}ccccc@{}}0 & & & & c_{1}\\ 1 & 0 & & & c_{2}\\ & \ddots & \ddots & & \vdots \\ & & 1 & 0 & c_{N-2}\\ & & & 1 & c_{N-1}\end{array}\right], & & \displaystyle\end{eqnarray}$$

where the only unknown are $(c_{1},c_{2},\ldots ,c_{N-1})$ , which are the coefficients of the polynomial

(5.13) $$\begin{eqnarray}\displaystyle \boldsymbol{v}_{N}=c_{1}\boldsymbol{v}_{1}+c_{2}\boldsymbol{v}_{2}+\cdots +c_{N-1}\boldsymbol{v}_{N-1}+\boldsymbol{r}, & & \displaystyle\end{eqnarray}$$

that express the vector $\boldsymbol{v}_{N}$ as a linear combination of the previous $N-1$ snapshots plus the residual vector $\boldsymbol{r}$ . By combining (5.10) and (5.13), it follows:

(5.14) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D63C}\unicode[STIX]{x1D651}_{1}^{N-1}=\unicode[STIX]{x1D651}_{2}^{N}=\unicode[STIX]{x1D651}_{1}^{N-1}\unicode[STIX]{x1D64E}+\boldsymbol{r}\boldsymbol{e}_{N-1}^{\text{T}}, & & \displaystyle\end{eqnarray}$$

with $\boldsymbol{e}_{N-1}$ the $(N-1)$ unit vector and $\unicode[STIX]{x1D651}_{2}^{N}=\{\boldsymbol{v}_{2},\boldsymbol{v}_{3},\ldots ,\boldsymbol{v}_{N}\}$ . The coefficient $c_{i}$ are to be evaluated using a QR-decomposition of the snapshots matrix $\unicode[STIX]{x1D651}_{1}^{N-1}$ . If $\boldsymbol{y}_{i}$ are the eigenvectors of $\unicode[STIX]{x1D64E}$ , it can be easily demonstrated that $\unicode[STIX]{x1D719}_{i}=\unicode[STIX]{x1D651}_{1}^{N-1}\boldsymbol{y}_{i}$ are the approximated eigenvectors of $\unicode[STIX]{x1D63C}$ .

Unfortunately, the implementation of the classical procedure yields an ill-conditioned algorithm. Therefore, the procedure used for the evaluation of the DMD mode is that based on the solution of the eigenvalue problem of a full matrix $\tilde{\unicode[STIX]{x1D64E}}$ which is more numerically stable (Schmid Reference Schmid2010). The algorithm consists of an economy-size singular value decomposition (SVD) of the snapshots matrix of $\unicode[STIX]{x1D651}_{1}^{N-1}=\unicode[STIX]{x1D650}\unicode[STIX]{x03A3}\unicode[STIX]{x1D652}^{\ast }$ ; substituting the economy-size SVD of $\unicode[STIX]{x1D651}_{1}^{N-1}$ in (5.14) leads to a full matrix $\tilde{\unicode[STIX]{x1D64E}}=\unicode[STIX]{x1D650}^{H}\unicode[STIX]{x1D651}_{2}^{N}\unicode[STIX]{x1D652}\unicode[STIX]{x03A3}^{-1}$ , where $\unicode[STIX]{x1D650}^{H}$ is a complex conjugate transpose of the basis $\unicode[STIX]{x1D650}$ . The DMD modes are then obtained from the eigenvectors $\boldsymbol{y}_{i}$ of $\tilde{\unicode[STIX]{x1D64E}}$ using the following expression:

(5.15) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D753}_{i}=\unicode[STIX]{x1D650}\boldsymbol{y}_{i}. & & \displaystyle\end{eqnarray}$$

The snapshot matrix $\unicode[STIX]{x1D651}_{1}^{N}$ can be decomposed in a general matrix form:

(5.16) $$\begin{eqnarray}\displaystyle [\boldsymbol{v}_{1},\ldots ,\boldsymbol{v}_{N}]=[\unicode[STIX]{x1D753}_{1},\ldots ,\unicode[STIX]{x1D753}_{N}]\left[\begin{array}{@{}ccccc@{}}\unicode[STIX]{x1D6FC}_{1} & & & \\ & \unicode[STIX]{x1D6FC}_{2} & & \\ & & & \ddots & \\ & & & & \unicode[STIX]{x1D6FC}_{N}\end{array}\right]\left[\begin{array}{@{}cccc@{}}1 & \unicode[STIX]{x1D707}_{1} & \ldots \, & \unicode[STIX]{x1D707}_{1}^{N-1}\\ 1 & \unicode[STIX]{x1D707}_{2} & \ldots \, & \unicode[STIX]{x1D707}_{2}^{N-1}\\ \vdots & \vdots & \vdots & \vdots \\ 1 & \unicode[STIX]{x1D707}_{N} & \ldots \, & \unicode[STIX]{x1D707}_{N}^{N-1}\end{array}\right], & & \displaystyle\end{eqnarray}$$

where the DMD mode $\unicode[STIX]{x1D753}_{i}$ identifies a spatial structure of the flow and the temporal evolution (either grow/decay rate and frequency) are governed by the Vandermonde matrix $\unicode[STIX]{x1D651}_{and}$ of the eigenvalues $\unicode[STIX]{x1D707}_{i}$ . The amplitude $\unicode[STIX]{x1D6FC}_{i}$ is obtained from the following optimization problem:

(5.17) $$\begin{eqnarray}\displaystyle \min _{\unicode[STIX]{x1D6FC}}\Vert {\mathcal{J}}(\unicode[STIX]{x1D6FC})\Vert =\min _{\unicode[STIX]{x1D6FC}}\Vert \unicode[STIX]{x1D651}_{1}^{N-1}-\unicode[STIX]{x03A6}\unicode[STIX]{x1D63F}_{\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D651}_{and}\Vert _{F}^{2}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x03A6}=[\unicode[STIX]{x1D753}_{1}(x),\unicode[STIX]{x1D753}_{2}(x),\ldots ,\unicode[STIX]{x1D753}_{N}(x)]$ and $\unicode[STIX]{x1D63F}_{\unicode[STIX]{x1D6FC}}=\text{diag}(\unicode[STIX]{x1D6FC})$ . The reader is referred to Rowley et al. (Reference Rowley, Mezic, Bagheri, Schlatter and Henningson2009) and Schmid (Reference Schmid2010) for further theoretical and implementation details.

It is worth noticing that, with respect to POD, in the DMD algorithm the mean is not first subtracted to flow data. In fact, it can be shown that the removal of the mean flow makes DMD identical to a temporal discrete Fourier transform (Chen, Tu & Rowley Reference Chen, Tu and Rowley2012), which is generally undesirable for non-periodic flows.

To select a subset of DMD modes that gives the most relevant contribution to the flow dynamics, Jovanović et al. (Reference Jovanović, Schmid and Nichols2014) proposed a sparsity-constrained variant of DMD, i.e. the sparsity promoting (SP) DMD. The implementation of the sparsity constraints involves the identification of the optimal amplitudes of the DMD modes by solving a convex optimization problem where the error norm (5.17) is minimized with a penalty parameter that induces sparsity in structure of the amplitude vector. The optimal amplitude of the extracted non-zero modes is evaluated by solving the following minimization problem:

(5.18) $$\begin{eqnarray}\displaystyle \min _{\unicode[STIX]{x1D6FC}}\left({\mathcal{J}}(\unicode[STIX]{x1D6FC})+\unicode[STIX]{x1D6FE}\mathop{\sum }_{i=1}^{N}|\unicode[STIX]{x1D6FC}_{i}|\right), & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FE}$ is the sparsity regularization parameter, such that higher values of $\unicode[STIX]{x1D6FE}$ give rise to a significant reduction of non-zero amplitude DMD modes, which has as counterpart a higher approximation error than that observed with lower $\unicode[STIX]{x1D6FE}$ values. The key element in the method is the definition of a trade-off between the number of extracted modes and the accuracy of the approximation of the reconstructed solution. More details about the implementation of the SP algorithm can be found in Jovanović et al. (Reference Jovanović, Schmid and Nichols2014).

As an example, for the more loaded propeller case ( $J=0.38$ ), in figure 9 the absolute value of the amplitude of the modes obtained without applying the SP algorithm is compared with the amplitude for different values of the SP parameter. In the discussion, the symbol $\unicode[STIX]{x1D705}$ indicates the frequency of the modes and is expressed in terms of blade passing frequency (BPF). As said above, the higher the value of $\unicode[STIX]{x1D6FE}$ , the smaller the number of relevant modes selected by the algorithm. However, it is possible to highlight a common trend for the different cases. In particular, starting from the almost uniform distribution of modes for $\unicode[STIX]{x1D6FE}=0$ , for a more selective choice of $\unicode[STIX]{x1D6FE}$ the spectra highlight the persistence of some dominant modes in the frequency range $\unicode[STIX]{x1D705}<4$ , whereas for $\unicode[STIX]{x1D705}>4$ only the multiples of the blade passing frequency outlive, with amplitudes that do not depend on the SP constant.

Figure 9.

DMD spectrum for $J=0.38$ .

Figure 10.

Performance loss (5.19) versus modes.

Obviously, the smaller the number of modes (i.e. the higher the value of $\unicode[STIX]{x1D6FE}$ ) the more the quality of the least square approximation (5.17) is degraded. In figure 10 the performance loss due to the SP algorithm, defined as

(5.19) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F1}_{loss}=100\sqrt{\frac{{\mathcal{J}}(\unicode[STIX]{x1D6FC})}{{\mathcal{J}}(\unicode[STIX]{x1D6FC})|_{\unicode[STIX]{x1D6FE}=0}}}, & & \displaystyle\end{eqnarray}$$

is depicted versus the per cent of non-zero amplitude modes. The symbols on the $J=0.38$ curve represent the three values of $\unicode[STIX]{x1D6FE}$ considered in figure 9. With the sole inclusion of the mean flow the performance loss is approximately 13 %, then it decreases and for $\unicode[STIX]{x1D6FE}=1400$ (that, in this case, corresponds to 93 modes) is 10 %. In figure 10, curve $J=0.71$ is also reported for comparison; it can be seen that $\unicode[STIX]{x1D6F1}_{loss}$ is around 3.6 % also with the mean flow alone, and weakly depends on the number of non-zero modes.

In order to remove transient modes which have high amplitudes but high decay rates, we further multiply each amplitude resulting from the SP algorithm by its damping factor $|\unicode[STIX]{x1D707}_{i}|^{N-1}$ (Statnikov, Meinke & Schröder Reference Statnikov, Meinke and Schröder2017). In fact, strongly damped modes influence only the early stage of time evolution. Since the aim of ROM is the identification of the most important modes for the entire available time history, the effect of these damped modes is limited and the use of their undamped amplitude can be misleading.

6.1 Case $J=0.71$

At $J=0.71$ , the propeller wake keeps stable for a distance which is longer than the computational domain considered in the analysis; this is consistent with the experimental visualizations presented in Felli et al. (Reference Felli, Camussi and Di Felice2011), where the onset of instability is observed at a distance of 6–8 diameters downstream the hub. In fact, as a consequence of the moderate loading, the distance between two consecutive tip vortices is sufficiently large to delay the disruption of the blade vortex sheet with the consequent inception of the destabilization.

At this loading condition, the flow past the propeller consists of organized and ordered motion of coherent structures that preserve their morphology during the downstream convection. DMD and POD analyses correctly identify this aspect of the flow. In figure 11, the amplitude of the computed DMD modes are depicted. The spectrum of the DMD shows the dominance of frequencies that corresponds to the BPF and its multiples; a set of modes in the range between the first and the second BPF with reduced amplitudes can also be observed.

Similarly, the POD analysis confers to only a few modes the most relevant contribution to the energy content of the flow. In particular, from figure 12(b), where the cumulative energy of the POD modes is reported, it can be noted that more than 85 % of the energy of the flow has to be ascribed to the first six modes.

Figure 11.

DMD spectrum for $J=0.71$ .

Figure 12.

POD relative (a) and cumulative (b) energy of the modes.

These modes of the POD decomposition identify the evolution of the coherent structures in the outer portion of the propeller slipstream and are associated with the dynamics of the tip vortices. Figures 13 and 14 show the three-dimensional representation of modes 1 and 3 in terms of constant pressure surfaces, and the correspondent auto spectral density (ASD), that is the Fourier transform of the autocorrelation function (Newland Reference Newland1993), of the temporal coefficients $a_{1}(t)$ and $a_{3}(t)$ . Clearly, the ASD associates the aforementioned modes to the first and second BPF and, consistently, the POD mode 3 is characterized by half the spatial wavelength of mode 1. Moreover, the shape of the mode is conserved for the entire length of the domain as a consequence of the stable nature of the downstream motion of the tip vortices. POD modes 5 and 7 (not shown for the sake of conciseness) are associated with the third and fourth BPF and show similar spatial features but length scales equal to $1/3$ and $1/4$ of the first mode.

DMD modes correspondent to tonal frequencies show the same spatial characteristics of the POD ones; in particular the richest dynamic content of the wake is limited to its outer portion affected by the passage of the tip vortices. For example, in figures 15 and 16 the DMD modes associated with the first and second BPF are shown; their similarity to POD mode 1 and 3, respectively, is evident.

Figure 13.

Case $J=0.71$ ; POD mode 1.

Figure 14.

Case $J=0.71$ ; POD mode 3.

Figure 15.

Case $J=0.71$ ; DMD mode $\unicode[STIX]{x1D705}=4$ .

Figure 16.

Case $J=0.71$ ; DMD mode $\unicode[STIX]{x1D705}=8$ .

As shown by Rempfer & Fasel (Reference Rempfer and Fasel1994), the POD eigenvalues occur in pairs of almost equal value. Each pair of eigenfunctions $\unicode[STIX]{x1D753}_{2n}$ and $\unicode[STIX]{x1D753}_{2n+1}$ provides the shape of a related coherent structure and its evolution while moving downstream. When the coherent structure is travelling at constant speed with constant shape, the eigenvalue problem degenerates, giving rise to pairs of identical eigenvalues. This is the case of the first six POD modes, which are characterized by equal energy content, spatial distribution and spectral content of temporal coefficient in pairs. The only significant difference is a temporal shift of their related coefficients, as shown, for example, in the case of $a_{1}(t)$ and $a_{2}(t)$ in figure 17 and by the phase portrait in figure 18. In this second figure, the data describe a circle, because they represent an oscillation made of two sine waves with equal frequency but a phase shift of $90^{\circ }$  (Oberleithner et al. Reference Oberleithner, Sieber, Nayeri, Paschereit, Petz, Hege, Noack and Wygnanski2011). This behaviour was already found in Sarmast et al. (Reference Sarmast, Dadfar, Mikkelsen, Schlatter, Ivanelli, Sørensen and Henningson2014) in the analysis of the flow behind a wind turbine. When the pair of eigenvalues are not exactly the same, because of the variation of the structures as they propagate, the related temporal coefficients show not only a time shift but also a small variation, as shown for example in figures 19 and 20. Moreover, in Schmid, Violato & Scarano (Reference Schmid, Violato and Scarano2012) it is shown that, for periodic and convective phenomena, a composition of a pair of (real) POD modes $\unicode[STIX]{x1D753}_{2n}+\text{i}\unicode[STIX]{x1D753}_{2n+1}$ coincides with a complex DMD mode, while for non-periodic or transient flows, this coincidence does not hold.

Figure 17.

Case $J=0.71$ POD coefficients $a_{1}$ and $a_{2}$ .

Figure 18.

Case $J=0.71$ phase portrait of coefficients $a_{1}$ and $a_{2}$ .

Figure 19.

Case $J=0.71$ POD coefficients $a_{9}$ and $a_{10}$ .

Figure 20.

Case $J=0.71$ phase portrait of coefficients $a_{9}$ and $a_{10}$ .

The remaining contribution to the modal description of the wake is related to flow structures that are poorly correlated to the tonal frequencies associated with the blade passage. This is evidenced by modes correspondent to non-integer frequencies, almost contiguously distributed in the frequency range between the first and second BPF in the DMD spectrum (figure 11). The amplitude of these modes is noticeably smaller than the tonal ones and, as alternatively described by the POD decomposition, provides a minimal contribution to the energetic content of the flow. In fact, modes higher than 8 contributes collectively to 15 % of the total energetic content of the flow (figure 12 b), clearly indicating that the flow is dominated by large-scale coherent motion.

The dynamic evolution of the hub region is described by POD modes starting from the ninth, with spatial and temporal characteristics comparable to those yielded by DMD. These modes consists of two systems of helical structures (figures 21–24 for the DMD and POD decomposition, respectively). Close to the hub, the modes look like a single helix with increasing pitch that, starting from the second half of the domain, evolve into a double helix with reduced pitch (i.e. reduced spatial wavelength). From the physical point of view, these modes show that in the inner region of the slipstream the perturbations are propagated by different wavelengths whose growth and decay rates overlap along the $x$ axis, confirming the results presented in Viola et al. (Reference Viola, Iungo, Camarri, Porté-Agel and Gallaire2014), Ashton et al. (Reference Ashton, Viola, Camarri, Gallaire and Iungo2016) for a wind turbine operating at equivalent loading condition. This aspect is supported by the presence of different frequencies in the ASD of the POD temporal coefficients $a_{9}(t)$ and $a_{11}(t)$ . Moreover, from these spectra it can be deduced that the secondary helix develops with greater strength, in particular for mode 9.

In addition to what already evidenced by the POD outcome, the DMD approach identifies a mode that describes the mutual interaction between the hub vortex and the tip vortex (figure 22). This mode is associated with $\unicode[STIX]{x1D705}=6.23$ (i.e. a frequency close to the average of the first and second BPF) and highlights the onset of a resonant phenomenon between the inner and the outer portions of the slipstream and the oscillatory nature of the hub vortex at frequency higher than the tonal one.

Figure 21.

Case $J=0.71$ ; DMD mode $\unicode[STIX]{x1D705}=6.23$ .

Figure 22.

Case $J=0.71$ ; DMD mode $\unicode[STIX]{x1D705}=6.00$ .

Figure 23.

Case $J=0.71$ ; POD mode 9.

Figure 24.

Case $J=0.71$ ; POD mode 11.

Figure 25.

Case $J=0.71$ ; correspondence of coherent structures and frequency content.

In order to inspect the multi-frequency nature of the hub vortex region, the POD analysis was also carried out for a shorter domain that excluded the region close to the propeller hub. This analysis (not entirely shown, for conciseness) permitted to associate the correct frequencies to particular structures of the modes, as shown in figure 25(a,b), that display, on a two-dimensional $x$ – $z$ section, the vertical velocity $w$ corresponding to mode 10 and 12. The coherent structure described by POD mode 10 is strong immediately past the propeller, because of the abrupt reduction of the hub sectional area. Slightly downstream, $w$ is markedly weakened for almost half of the domain and, then, reinforces with a different morphology and spatial wavelength that looks compatible with the double helix shape obtained by the three-dimensional analysis. Specifically, the peak close to $\unicode[STIX]{x1D705}=5.8$ can be ascribed to the perturbation close to the hub and the one close to $\unicode[STIX]{x1D705}=6.5$ to the fluctuations in the far field. Conversely, mode 12 in figure 25(b) evidences a different strength of the spiral: in this case, the fluctuations are stronger about the middle of the domain where the pattern resembles a single helix that gradually changes into a weaker double helix shape downstream.

The modal description, integrated by simple momentum based considerations, can give some details on the hub vortex dynamics and its relation with the flow inside the propeller slipstream. In particular, immediately past the propeller, the flow around the hub vortex is strongly perturbed due to the radial displacement of the root vortices that follow the hub surface and wrap around a region of decelerated flow around the propeller axis. These perturbations are immediately weakened by the momentum imparted by the propeller to the incoming flow and by the concurrent contraction of the propeller slipstream that further limits the radial deflections of the vortical structures. The acceleration of the flow inside the propeller slipstream is consistent with the fact that the single helix develops with an increasing pitch angle. However, as a consequence of the mild propeller loading, the length scale of the momentum gain is relatively short and, hence, the hub vortex fluctuations restart around the middle of the shown domain in a different form compatible to the double helix oscillations.

6.2 Case $J=0.38$

At this condition, the flow complexity enriches the modal representation obtained by the POD and DMD decompositions. In particular, POD analysis highlights a significant redistribution of energy content among the modes with respect to $J=0.71$ , as depicted in figure 12(a). It is also interesting to notice that the cumulative energy distribution (figure 12 b) at $J=0.38$ is characterized by a wider spread of energy among the modes: more than 120 modes are necessary to attain 90 % of the total energy (only 16 modes are required by the more stable condition to cumulate the same amount of energy). These trends are a clear indication that, in spite of the persistence of large structures, a significant part of the kinetic energy stems from high frequency and small-scale turbulence.

The complexity of the flow condition is even better captured by the DMD, which highlights additional characteristics of the flow. The analysis of the spectrum of the SP DMD, in particular the one in figure 9(c) obtained with a high but not too selective value of $\unicode[STIX]{x1D6FE}$ , allows us to identify five different set of modes, of which one is related to the blade passing frequency and the others are associated with specific aspects of the destabilization mechanism, as shown by their spatial morphology. The reader is referred to the temporal sequence in figure 7 and the slipstream flow depicted in figure 8 to better follow the analysis. The five different sets of modes identified in the DMD analysis are presented and compared/integrated with the results obtain by POD.

6.2.1 Case $\unicode[STIX]{x1D705}=4$ and multiples

This group of modes is related with the tonal frequencies, associated with the passage of tip vortices detached from the propeller blades. For these modes, the results obtained using both DMD and POD analyses are in perfect agreement. In particular, both approaches identify as the most energetic mode the first blade passing frequency $(\unicode[STIX]{x1D705}=4)$ . In figure 26 the pressure isosurface of the DMD mode at this frequency is depicted. Both DMD and POD analyses point out a contraction in length of the mode, when compared to the one with the same frequency in the more stable condition $J=0.71$ . In fact, this mode is stronger in the near wake of the propeller (up to 2.5 propeller radii) and weakens in the far field, where the wake is altered due to the destabilization (downstream the red line in figure 26).

When compared to the higher advance ratio condition, for $J=0.38$ the POD analysis reveals a significant loss of the energy associated with the two POD modes corresponding to $\unicode[STIX]{x1D705}=4$ . Indeed, their contribution to the whole energy content is around 32 % (see figure 12 a) against 70 % of the more stable condition. Also the energy associated at higher harmonic is significantly lower: 8 % and 4.6 % for $\unicode[STIX]{x1D705}=8$ and $12$ , respectively.

Figure 26.

Case $J=0.38$ ; DMD mode $\unicode[STIX]{x1D705}=4$ .

Figure 27.

Case $J=0.38$ ; DMD set of modes around $\unicode[STIX]{x1D705}=2$ .

6.2.2 Case $\unicode[STIX]{x1D705}=2$

The modes at this frequency are associated with the pairing process of two consecutive tip vortices that occurs in the transition region of the domain, as it is evident from the inspection of figure 27(a). The visualizations by $\unicode[STIX]{x1D706}_{2}$ criterion (figure 6) show a very complex evolution of the tip vortices (i.e. propagation of small waves along the vortex filament, with abrupt change of pitch); all these events are partially accounted for by the modes close to the reference one at $\unicode[STIX]{x1D705}=2$ . In fact, their shape highlights only the perturbation of the spirals at the boundary of the slipstream, the perturbation of the hub vortex being not visible. The modes at frequency higher than 2 (figure 27 c,d) show a reduced pitch and a more evident increase of the external diameter. These modes can be associated with the global burst of the propeller slipstream occurring after the first coupling process develops (see figure 9 in Felli et al. (Reference Felli, Camussi and Di Felice2011)). The gradual increase of the sectional area stretches the vortex filaments outwards and consequently the pitch reduces; this effect is also concomitant to the averaged reduction of the speed within the propeller slipstream to satisfy mass conservation. The reference mode at $\unicode[STIX]{x1D705}=2$ clearly shows a spatial frequency of the tip vortices doubled with respect to the first BPF and, in particular, the pitch of the spirals is regular and azimuthally constant. On the other hand, the modes associated with neighbouring frequencies describe the displacements about this mean condition, and can be related to pairing mechanism described in the sequence of figure 7. In particular, the sequence shows that the pitch varies azimuthally for both the merged and individual structures. This variable pitch distribution is accounted for by modes at neighbouring frequencies (figure 27 b–d). All these modes recall the non-constant distribution of the pitch during a complete period of the helix but centred with respect the one at $\unicode[STIX]{x1D705}=2$ . The marked alteration of the pitch is obviously related to the fact that the deformation, straining and merging of the interacting vortices are localized processes.

The POD analysis adds some details to the above description. In particular, modes 3 and 4 exhibit two similar peaks at 1.95 and 2.05, correlated with the pairing of adjacent tip vortices (figure 28). The energy content associated with these modes is 8 % of the total energy (figure 12 a), that is a fourth of that associated with the first two modes.

Figure 28.

Case $J=0.38$ ; POD mode 3.

Figure 29.

Case $J=0.38$ ; DMD set of modes around $\unicode[STIX]{x1D705}=1$ .

6.2.3 Case $\unicode[STIX]{x1D705}=1$

This set of modes is associated with the destabilization of the hub vortex and its effects on the outer structures of the wake. These modes highlight an increasing oscillation of the hub vortex (figure 29), in accordance with the recent stability analysis for the wake of wind turbines (Ashton et al. Reference Ashton, Viola, Camarri, Gallaire and Iungo2016). This mode is characterized by a considerably larger displacement of the outer structures of the wake (resembling the merged tip-vortex structures) with respect to the modes at $\unicode[STIX]{x1D705}=2$ . These observations relate this mode to the precession motion of the propeller slipstream observed in the experiments by Felli et al. (Reference Felli, Camussi and Di Felice2011) and confirmed by experiments on wind turbines. The POD mode that is closer to this frequency is the seventh mode (see figure 30), which is characterized by a double peak in the coefficient spectra at $\unicode[STIX]{x1D705}=0.9$ and 1.05, but the associated coherent structure has a different morphology with respect to the corresponding DMD mode. This is due to the multi-component frequency spectra of the POD coefficient, as it can be evidenced also by the presence of a components at $\unicode[STIX]{x1D705}=1.3$ and 2.

Figure 30.

Case $J=0.38$ ; POD mode 7.

6.2.4 Case $1<\unicode[STIX]{x1D705}<2$

Modes at intermediate frequencies, between the previous two last sets, can be ascribed to a sort of resonance caused by the interaction of the merged tip vortex system with the hub vortex. The information synthesized by these modes allows us to gain a deeper insight on the propagation of the destabilization from the outer to the inner structures of the propeller slipstream.

The morphology of these modes changes remarkably with the frequency (figure 31). Modes close to $\unicode[STIX]{x1D705}=1$ (e.g. $\unicode[STIX]{x1D705}=1.22$ ) are characterized by external structures with small pitch, resembling vortex rings, and a strong fluctuating structure in the hub region. With the increase of $\unicode[STIX]{x1D705}$ , the central structure progressively weakens and the external ones exhibit strong variations of the pitch along the azimuthal position. This latter characteristic underlines the deformation and the straining of the tip structures already observed for the modes at $\unicode[STIX]{x1D705}=2$ and associated with their interaction and pairing.

These modes confirm what suggested by both numerical and experimental visualizations, namely that the core of the slipstream starts to oscillate after the destabilization of the outer region. In fact, as the pairing process starts, the circumferential variation of the pitch of outer filament causes the unbalance of the induced velocities to the hub vortex and hence promotes its destabilization. At first, this mechanism slowly builds up and the fluctuations of the hub become non-negligible downstream of tip vortices pairing, as shown by the systematic analysis provided in Felli et al. (Reference Felli, Camussi and Di Felice2011), figure 10. This delay can be ascribed to the stabilizing effect exerted by the acceleration within the propeller slipstream and its consequent contraction. Due to the higher loading, this effect is stronger than for $J=0.71$ . This inertial effect persists for a longer period, because, when the tip vortices merge, the hub vortex evolve as an independent structure, the blade vortex sheet being cut (Felli & Di Felice Reference Felli and Di Felice2005; Felli et al. Reference Felli, Camussi and Di Felice2011); on the contrary, the outer region results more prone to destabilization because of the shear effect caused by the contact with the undisturbed, slower flow.

On the basis of the above description, the mode at $\unicode[STIX]{x1D705}=1.46$ can be associated with the mutual and self-sustained excitation between the merged tip vortices and the hub vortex. These interactions lead to the physical behaviour described at lower frequencies: the oscillation induced on the outer structure amplifies the shear at the boundaries of the propeller slipstream that, in its turn, dampens out the hub vortex oscillation and frequency.

Interestingly, the modes closer to the resonant frequency are the most important ones, as it can be deduced by the disappearing of modes at frequencies at the extremities of the interval (see spectra for $\unicode[STIX]{x1D6FE}=2400$ ). Moreover, the modes that excite the hub vortex are those that show a non-constant pitch in the azimuthal direction, that cause an unbalanced perturbation to the inner vortex structure. Approaching the frequency of the second BPF, the shape of the correspondent modes reveals the progressive weakening of the hub vortex.

Figure 31.

Case $J=0.38$ ; DMD modes with frequencies between $\unicode[STIX]{x1D705}=1$ and $\unicode[STIX]{x1D705}=2$ .

6.2.5 Case $\unicode[STIX]{x1D705}<1$

Modes at frequencies lower than 1 appear in the DMD analysis only for moderate values of $\unicode[STIX]{x1D6FE}$ , and also in the POD but with a very small energy content. Representative modes of this group are shown in figure 32(a,b); their features suggest their connection with the precession motion of the slipstream considered as a vortex tube. In particular, in the outer region of the slipstream, the tip vortices are not distinguishable; instead, a homogeneous sheet generated by the shear forces (arising from the contact with the external flow) can be observed. The shear, on its turn, dampens the precession motion excited by the hub vortex. The POD analysis (figure 33) allows us to capture this particular feature of the flow, and provides also an information concerning the associated energy, which is quite small (0.7 %). Moreover, it can be observed that the perturbation of the hub vortex is visible only at the end of the domain, further downstream with respect to modes related to hub vortex excitation by the tip vortex system (i.e. $1<\unicode[STIX]{x1D705}<2$ ).

Figure 32.

Case $J=0.38$ ; DMD modes with frequencies $\unicode[STIX]{x1D705}<1$ .

Figure 33.

Case $J=0.38$ ; POD mode 21.

6.3 Further comments on the hub vortex dynamics

It is interesting to note the different shape of the modes associated with the dynamics of the hub vortex between the two loading conditions. At $J=0.71$ the modes of the hub vortex show a transition from a single helix close to the propeller to a double helix close to the end of the domain. On the other hand, at $J=0.38$ the modes are shaped like a single helix in the inner region (and someway coupled with the outer region of the slipstream). This difference can be explained by considering the relation between the advance coefficient (propeller loading) and the geometry (pitch angle) of the vortex structures shed by the blades, in addition to the comments on the momentum/energy transfer imparted by the propeller discussed above. Specifically, the pitch angle is strictly related to the three-dimensional features of the induced velocity field within the slipstream. For example, during very high loading conditions the wake structures resemble vortex rings and the induced velocity has a dominant axial component; conversely, as the propeller loading is gradually reduced, the helix elongates and the tangential component increases. At $J=0.71$ , the oscillations of the hub vortex are not damped, because of the weak acceleration of the flow; moreover, they can be amplified by a non-negligible tangential velocity induced by the vortex structures in the outer region of the slipstream. On the contrary, at $J=0.38$ , the momentum transfer is remarkably stronger in the wake and no evidence of such a transition is visible in the inner region, the prevailing disturbance being caused by the coupling process of the tip vortices.