Analysis of Dynamic Stall on a Pitching Airfoil with Spectral Proper Orthogonal Decomposition

In the field of Helicopter Engineering the need for a simultaneous variation of the angle of attack of the blade during its rotation about the hub is the cause of a major issue known as dynamic stall. This event is responsible for significant structural problems due to the periodic presence of augmented loads conditions that may couple with the rotor dynamics and cause the earlier failure of the blade. Clearly, this represents a propulsive limit for the craft and a deeper knowledge of the phenomenon is necessary to devise methods for controlling or preventing it. Among the several techniques for analyzing it, the Proper Othogonal Decomposition has shown to provide an unbiased tool for the identification of coherent structures underlying the flow evolution. In this work the spectral variation of the decomposition, based on the formulation from Sieber (2016) et al. [37], is used to analyze the flow. This method takes advantage of the application of a filter that modifies the original formulation in order to extract new structures otherwise interpreted as noise of the considered data. The case study discussed here is that of a 2D pitching NACA 0012 airfoil, under deep stall conditions. The numerical data obtained by the CFD simulation with the DLR TAU code are decomposed with the stochastic tool. In particular two fields are decomposed: the vectorvalued velocity field and the scalar valued pressure field. New coherent structures are identified through the application of the filter and the comparison between the two analysis provides a criterion for choosing a proper filter in order to approach the dynamic stall phenomenon. The dynamics underlying the most correlated structures shows a higher sensitivity towards the filter adopted when one analyzes the pressure field, which depicts lower responses at higher frequencies if compared with the velocity field. Furthermore, the computation of the aerodynamic loads from a partial reconstruction of the pressure field is adopted to discuss the possible applicability of the spectral variation in the context of reduced order modeling.


Introduction
The relative motion of the flow past a helicopter rotor during forward flight causes a differential load distribution over the blades, due to a dynamic pressure that increases for the blades undergoing the advancing phase and, simultaneously, decreases for the blades undergoing the retreating phase [47]. The roll moment that would develop accordingly is therefore prevented by means of a change of the angle of attack of the blades, the pitching motion being one of the solutions applied. The angle of attack is increased and decreased, respectively, for the retreating blade and the advancing blade, describing a periodic motion that develops with the same period as the rotor revolution, producing rapid variations of the incidence of the flow past the airfoils [10]. These variations, when leading the airfoil to exceed the static stall angle of attack, cause the occurrence of dynamic stall: this event is characterized by sudden excursions in the aerodynamic loads, which reach considerably higher peaks if compared with the behavior of the static stall and, more, can couple with structure dynamics (e.g. emphasizing the dynamic flutter), thus anticipating its fatigue failure [13]. Hence, this phenomenon represents the main issue in the design of a helicopter rotor and also defines its propulsive limit: this gives an idea on why the airspeed performances in helicopters engineering is so dependent on a good comprehension of the dynamic stall phenomenon and justifies such a great effort of the research world toward this topic [22]. Many experiments were performed so far in order to achieve a better description of the events occurring during dynamic stall: among them, the results that deserve a particular mention are the ones from Carr et al. (1977) [7] and McCroskey et al. (1978) [28], as well as the ones from Gerontakos in more recent years (2004) [18]. Based on the former, a detailed description of the physics evolving in the context of dynamic stall over pitching airfoils is provided by McCroskey (1981) [29]. With regard to numerical investigation, the unsteadiness typical of this phenomenon, as well as the turbulence and the laminar-to-turbulent transition involved, represent a significant issue in terms of computational costs. This fact explains why satisfying results were made possible only in recent years, thanks to the improvements made either in the computers performances and in the numerical techniques applied to the Computational Fluid Dynamics (CFD) [10]. The need to achieve a better ability to simulate and, so, predict the phenomenon led, over the last decades, to a rich production of publications, whose main target is to provide the tools for distinguishing the best approach when one has to deal with dynamic stall. The various turbulence models are intensively tested and discussed in literature together with the space-time discretization of the numerical model they are applied to. A detailed study can be found in Singh et al. (2019) [38], where the authors compare the behavior of seven different turbulence models for a 2D case; another accurate analysis is the one performed by Wang et al. (2010) [45] and Wang et al. (2012) [46] in which attention is paid to dynamic stall occurring at low Reynolds numbers and the results are meticulously described for either a 2D case and a 3D case. In Kai et al. (2017) [21] the model is first improved to obtain more accurate results than the ones from Wang et al. (2010) [45] and then a plunging motion is coupled with the pitching motion in order to analyze the possible effects.
Marchetto et al. (2019) [27] study the impact of a Gurney flap on dynamic stall, after validating their turbulence model against experiments, for a baseline airfoil, by evidencing strengths and weaknesses of the CFD approach. Understanding the peculiarities of dynamic stall is just the drawback of a broader problem that has committed the fluid dynamic community for many decades: understanding the essence of turbulence. Since the early 1920s a great effort has been put by researchers in the study of this topic, but the fact that a general model for turbulence is not yet available suggests that the phenomenon itself is not yet sufficiently understood. The original concept of turbulence as a stochastic phenomenon characterized by the composition of a mean replicable flow field and a randomly fluctuating one has been gradually refined during the years with further inspections that have determined the presence of a certain organization even within the chaotic component [6]. In order to better identify these features, since its introduction to the turbulence context thanks to Lumley (1967) [25], the Proper Orthogonal Decomposition (POD) has represented a powerful stochastic tool, able to recover the so called coherent structures within time dependent turbulent flows, first of all because of its unbiased approach which made no more necessary to know these structures a priori [1,16,2]. Since then the method has been widely adopted, either in experimental and numerical analysis, for the study of turbulent flows and many results can be found in literature. During the years, POD has also undergone different interpretations, among which it is worth to mention the one from    [23] for the case of a jet flow in the context of swirl-stabilized combustors. Ricciardi et al. (2018) [33] applied SPOD to the numerical study of the fluid stream past a landing gear, in order to identify the turbulent structures causing noise in the internal cavities of the wheels. In the context of flows past wings, analysis were made for different cases. Ribeiro et al. (2017Ribeiro et al. ( , 2019 [31,32] studied the flow past a blade with a NACA 0012 airfoil section: the authors applied the tool to their numerical data, analyzing coherent structures either for the 3D data and the the 2D data obtained by spanwise averaging the former. In Steinfurth et al. (2018) [43] SPOD is exploited to compare the different behaviors of the flow, in the experimental case of a blade subject to different active wake control configurations. Dynamic stall was investigated by means of SPOD by Wen et al. (2018) [47]: the numerical data for the analysis were taken from simulations of a blade and then spanwise averaged, for three different configurations of motion, namely: pitching-surging, pitching-surging-rotation, pitching-surging-yawing. The purpose of this work is to analyze the dynamic stall phenomenon by applying the SPOD to the numerical data obtained by CFD simulations of a 2D pitching NACA 0012 airfoil. The analysis takes into consideration two different types of flow field: the vector-valued velocity field and the scalar-valued pressure field. A comparison between the different results allows for a better understanding of either the peculiarities of each field and how a proper SPOD filter should be chosen in order to extract new structures from fields of the same kind in the context of deep dynamic stall. In particular, the discussion is thus organized: in Chapter 2 the phenomenon of dynamic stall is identified as a consequence of the helicopter rotor dynamics and described in terms of physical events. Chapter 3 provides an introduction to the POD theory and then SPOD main concepts are drawn in order to highlight similarities and differences with the former and to show its potential.
Chapter 4 presents the numerical model adopted. In particular, it is described in terms of geometry and spatial discretization; its suitability to reproduce the dynamic stall phenomenon is addressed by comparing its prediction of the loads coefficients with the experimental measurements from Gerontakos (2004) [18], for two different turbulence models. In the last section another flow configuration is described according to the flow field behavior and corresponding effects on the loads curves. The numerical data obtained from this latter simulation represent the input for the modal analysis performed with the SPOD. In Chapter 5 the stochastic tool is applied to the data related to a selection of the whole domain, corresponding to 2 cycles out of the 4 performed. The time and spatial correlations of the two fields are identified and discussed. The loads coefficients are computed from a partial reconstruction of the pressure field, with two different application of the SPOD, and compared with the original data from the CFD solver. In the concluding part the results are further discussed and compared.

The Rotor-Craft Dynamics
Differently from a fixed wing operating conditions, which are usually characterized by a steady freestream, a rotating wing mostly operates within unsteady inlet flows [47]. This is even more clear if one considers the motion of a rotor blade in forward flight ( Fig. 2.1): the local velocity of the relative wind at a given section is obtained (if neglecting other effects) by the composition of the circular velocity and the induced velocity of the advancing motion of the rotor. While the former is fixed for a chosen radial position, the latter changes according to the phase angle of the rotating motion of the blade, namely the blade azimuth ϕ [11]. As regards the blade undergoing the advancing phase (0°< ϕ < 180°), therefore called advancing blade, the relative wind results in a velocity higher than the flight speed U f f due to the addiction of the circular velocity, with the peak velocities experienced by the tip region of the blade at ϕ = 90°. Accordingly, the opposite behavior occurs for the retreating phase (180°< ϕ < 360°): here, the freestream velocity is lower than the flight speed until the point that a reverse flow region arises close to the hub, at approximately ϕ = 270°, with the fluid flowing from the trailing edge to the leading edge [13]. These conclusions find validation if one considers the definition of the lift produced by a blade during forward flight, L b [5]: where: ρ is the fluid density, c is the chord length, U T is the hover tip speed, e is the root cutout radius, R is the blade radius, c l is the section lift coefficient, r is the blade radial location and µ ∞ is the advance ratio, defined as [47]: Ω being the rotor speed. In fact, while the quantities before the integral are constant during the revolution, the integral itself, which can be defined as the mean blade lift coefficient C l [5], is dependent on the blade azimuth. In particular it is easy to state that this dependency presents a maximum and a minimum coincident with the angles previously mentioned, causing different lift forces on two opposite blades. This unbalance of the air loads occurring at diametrically opposed positions, together with the lever arm made by the blades, generates a roll moment on the whole system. A trim is then necessary to prevent this couple, which is done by means of a change of the angle of attack α of the blades, the section lift coefficient c l being strictly dependent on it [11]. Hence, in order to recover the same lift on the two blades, the angle of attack must be increased for the retreating blade and, accordingly reduced for the advancing blade: this is usually done by means of a pitching motion applied to the blades, whose incidence is changed periodically during the rotation and can be described with the following [13]: where: α 0 is the mean angle of attack and α 1 is the pitch amplitude.
When an unsteady motion applied to a lifting surface, such as pitching motion, causes it to exceed its static stall angle of attack, the dynamic stall phenomenon occurs, producing strong excursions in the aerodynamic loads which can couple with the structure dynamics and lead it to earlier failure [22,13]. Furthermore, due to the presence of the reverse flow region, a negative lift region arises close to the hub: although the values of the force here produced are lower than the ones detectable at the tip, it is important to account for this loss of lift. In fact, according to this behavior of the flow, the blades are usually twisted such that the angle of attack at the tip is higher than one at the hub, even though the pitching motion amplitude decreases from the hub to the tip [10]. Hence, it is reasonable to expect dynamic stall to occur close to the tip, that is where the angle of attack reaches the highest values within the evolution of the motion. In particular, Bousman (1998) [4], as part of an experimental investigation on the behavior of the air loads for three different flight conditions of the UH-60A helicopter, showed that the region of the blade mostly affected by this phenomenon was that one in the range between r = 0.77R and r = 0.92R ( Fig. 2.2) during the retreating phase. In fact, the author discovered that the airloads on the outboard section of the blade, although characterized by higher incidence angles that lead the flow to be supercritical, is mostly affected by the separation caused by the blade tip's vortex than by dynamic stall.

The Flow Morphology
As aforementioned, the dynamic stall is a phenomenon that occurs when a lifting surface undergoes a rapid change of its incidence up to angles of attack that are beyond its static stall limit. When static stall occurs, the aerodynamic coefficients show sudden excursions due to a separation of the flow over the suction surface, with a decay of the lifting load and, accordingly, a rise of the drag force. While this phenomenon is more typical of low airspeeds, the dynamic stall affects rotor blades undergoing fast streams since they operate close to the limits for the detachment of the flow [22]. One of the main consequences of the phenomenon, i.e. the increase of the maximum lift achievable with an airfoil in quasi-static conditions, was already observed in early 1930s but the first experimental data on the effects of dynamic stall were obtained in 1960 [22,7]. Since then, the dynamic stall has been involved in an intense process of investigation and many results produced during the 1970s and 1980s are still today the reference points for the research environment. The majority of the knowledge on this phenomenon is to be attributed to 2D pitching airfoils experiments, even though dynamic stall has a 3D character [22]; however, since the present work focuses onto the bi-dimensional effects, it is not its purpose to describe the physics undergoing the 3D behavior of the event.
The structure at the basis of dynamic stall is a vortex-like disturbance, which develops and sheds on the suction side, causing significant alterations to the forces and moments acting of the airfoil [29]. As the vortex leaves the airfoil, the flow becomes fully separated: this condition persists for a part of the cycle, until it reattaches again. This behavior, when viewed against the angle of attack, defines clear hysteresis curves in the airloads plots [7]. The effects thus produced allows to define two different regimes of dynamic stall (Fig. 2.3): • Light stall: it can be detected when the peak angle of the motion is only slightly greater than the static stall angle [49]. The flow operates only partially in conditions of full separation, which allows to reach airloads peaks higher than those of the steady case, but still similar. Clear hysteresis loops are detectable, especially for the moment coefficient which shows considerably negative values and negative damping phases (described further on) [10]. Moreover, the vertical dimension of the boundary layer is on the order of the airfoil thickness and smaller than the that one associated with static stall, and it is particularly influenced by the airfoil geometry, the reduced frequency k and the Mach number M [29,49]; • Deep stall: it occurs when the motion develops with a mean angle of attack close to the static stall angle and the amplitude of the motion is such that the latter is exceeded by at least 5° [49]. The flow is mainly under conditions of separation. This fact produces severe peaks in the airloads coefficients, as well as large hysteresis loops which result in a high aerodynamic damping impact [10]. As regards this regime, the dimension of the boundary layer increases up to the order of the chord length of the airfoil, showing a minor dependence on the geometric and motion parameters [29,49].  [29]. The plots of forces and moment coefficients are shown for both the light stall regime (left) and deep stall regime (right). The mean angle is written above each set of plots, while the pitch amplitude is always 10°.
In general, although more clearly noticeable from the deep stall graphs, the dynamic stall process develops through six main stages which can be derived from the trend of the airloads during a pitching cycle, as proposed by McCroskey (1981) [29]. Carr et al. (1977) [7], based on 2D experiments on an oscillating NACA 0012 airfoil, provide a thorough description of these events as well as a careful analysis of flow behavior from one step to the following. In Figure 2.4 every particular flow condition is referred with a letter and depicted in terms of either the effect on normal force and moment coefficients (respectively, c N and c M ) and the flow distribution on the airfoil. The first part of the cycle, called upstroke, is dominated by the linear trend of the thin airfoil theory in static conditions, the flow being unseparated even when the angle of attack exceeds the static stall condition in (a) [29]. This delay in the separation is reported by Corke et al. (2015) [13] as a consequence of two effects: the increase in the effective camber, occurring when the rate of change of the incidence angle is positive, and the Magnus effect caused by the pitching motion of the leading edge. The boundary layer begins to thicken in (b), showing the first appearance of flow reversal regions arising close to the trailing edge until in (c) the first eddies develop according to the wavy pattern in the furthest zone. The spread of the flow reversal over the suction side (d) until the initiation of the dynamic stall process is determined by the sudden separation of the boundary-layer near the leading edge, which states the formation of the so called Leading Edge Vortex (LEV). This vortex, as reported by McCroskey (1981) [29], starts to shed rearward at a speed lower than half the freestream velocity of the 2D configuration, U ∞ , causing the typical excursions in the airloads coefficients plots. In fact the core of the vortex, while shifting toward the trailing edge, acts as a low-pressure source which is the cause of either an abrupt increase of the normal force (f) and a sudden drop of the moment curve, referred to as moment stall (g). This latter is due to the nose-down effect produced on the airfoil as the shedding vortex changes the pressure distribution on the suction side and represents a clear difference with respect to the static stall, in which moment stall and lift stall occur simultaneously. In fact, lift stall (h) begins only later, as the vortex core nears the midchord: starting from this event the behavior of both the curves shows a similar dropping trend, until the LEV core is almost at the trailing edge. Here, the moment stall reaches its maximum negative value (i) and starts to rise rapidly to values similar to the ones of the static condition, while the normal force is still decaying steeply. The vortex sheds into the wake right after, in the first moments of the downstroke phase. The flow enters the fully separated condition in (j) showing a behavior comparable with that of an airfoil under static stall condition, unless secondary or even tertiary vortices arise, causing additional fluctuations on the airloads trends. The boundary layer then begins to reattach (k) gradually from the leading edge, until the reattachment is complete: this occurs during the downstroke, but the linear behavior is recovered only some moments after the initiation of the upstroke (l). From this event, the cycle starts again and repeats with the same pattern [13,7,29]. One of the main drawbacks of this highly nonlinear behavior of the airloads is related to the stall flutter which can be defined as a single-degree-of-freedom oscillatory motion, differently from the classical flutter occurring in unseparated flows which usually involves more degrees of freedom [29]. This phenomenon excites the natural frequency of the torsional mode and can produce substantial increases of the limit-cycle vibrations amplitudes up to values that can lead to divergent, and so fatal, oscillations [22]. The stall flutter directly relates to the work exerted by an airfoil on the fluid due to its motion, which, for a pitching profile, depends on the pitching moment about the axis of rotation [29]. This value is usually nondimensionalized, which leads to the definition of the damping factor symbolized by D.F. or C W and given by [22]: Considering, as above, the plot of c M vs α, it is straightforward to conclude that Equation 2.4 has the meaning of the area enclosed by the pitching moment curve throughout a cycle and, in fact, the damping coefficient is experimentally calculated as [13]: D and U denoting respectively the downstroke and upstroke phases.
In accordance with the definition of C W , it is clear that a positive damping means an introduction of energy made by the body onto the fluid which, in fact, acts as an actual damping for the airfoil motion; furthermore, from Equation 2.4 one can conclude that a positive damping is represented by a counterclockwise closed loop in the pitching moment trace or, put another way, a region enclosed by an upstroke phase segment, as upper bound, and a downstroke phase segment, as lower bound [10]. A negative damping coefficient, instead, has exact opposite features so the fluid acts as an energy source for the body motion, promoting the divergence of aeroelastic oscillations [22]. Throughout dynamic stall, both the conditions occur. If the contribute of the positive damping phases is such that the negative damping condition is not balanced, new energy is added to the airfoil motion after every cycle, resulting in a limit-cycle growth up to exaggerated amplitudes [13]. It has been noticed that an increase in the mean angle of attack causes an enlargement of the area related to negative damping; however, this behavior holds until the values of α 0 are consistent with the onset of deep stall: in fact, deep stall penetration introduces a new counterclockwise loop which increases the amount of positive damping, thus stabilizing the motion. Then, the reduced frequency is increased in order to delay the occurrence of stall at high incidence angles, which reduces the portion of flow separation and justifies a higher aeroelastic stability of the deep stall regime if compared with the light stall regime ( Fig. 2.5) [22,13].  [13]. Negative damping (orange) and positive damping (gray) regions are identified within a cycle for both the light stall (left) and the deep stall (right).

Numerical Modeling
Due to the difficulty in reproducing experimentally every nonlinear aspect of the flight conditions of helicopters and the regimes of large turbulent separated regions appearing throughout dynamic stall, reasonable predictions on the phenomenon can be achieved only by means of a numerical solution of the equations underlying it physics. As a fluid dynamic problem, dynamic stall is governed by the Navier-Stokes equations, which describe the behavior of a fluid, encompassing the conservation of mass and momentum of the domain studied and the way energy interchanges within it, as follows [22,30]: where u = u i = (u, v, w), with i = 1, 2, 3, is the velocity vector defined for the three directions of the space x = x i = (x, y, z), p is the pressure, T = τ i , j is the stress tensor, h is the specific enthalpy and q = q i is the specific heat flux vector; Re is the dimensionless Reynolds number, defined as: with V being a characteristic velocity of the flow, being a characteristic length of the flow, µ being the dynamic viscosity of the fluid and ν being the kinematic viscosity of the fluid; while Pr is the dimensionless Prandtl number, defined as: with c p being the specific heat of the fluid and K being the thermal conductivity of the fluid. Although analytical solutions are available for System 2.6, they are only applicable to simple cases, which makes it impossible to exploit them for solving practical problems. These latter, instead, are solved with a numerical approach that requires the space and time discretization of the investigated physical domain [3]. The solution thus obtained is then discrete: differently from an analytical solution, which provides a function describing the behavior of the dependent variables continuously throughout the space-time domain, a numerical solution can be defined as a set of "numbers" relating to specific points (called grid or mesh points) at given time steps, according to the discretization adopted [48]. Despite Direct Numerical Simulation (DNS) being a possible way to solve numerically the Navier-Stokes equations, it is known that the number of grid points needed to achieve enough space resolution scales as Re 9 4 and the computational time as Re 3 ; therefore, according to the typical operating ranges of the Reynolds number in the helicopters environment (10 5 <Re< 10 6 [22]), it appears clear that such an approach is not affordable even by the modern supercomputers [3]. For such a reason, a common way to approach the solution of System 2.6 is by modeling certain features of the flow through new sets of equations that require the definition of specific constants for the closure of the problem, which have to be determined: an example is given by the Boussinesq assumption for modeling the Reynolds stress tensor [10]. A first class of approximation strategies includes the so-called Large Eddy Simulations (LES) approach, which takes advantage of the Kolmogorov theory (1941) for modeling the small turbulence scales, characterized by a more universal nature: the large scales are instead accurately resolved, which requires a discrete resolution of the grid but still less demanding than a DNS computational domain. However, its inherent three-dimensional and unsteady character makes it still too computationally onerous if one considers the accuracy required in many engineering applications, but it represents a promising tool if specific conditions were to need a more careful investigation [3].
Instead, the most common method in engineering environment, according to the averaging process is based on, goes by the name of Reynolds-Averaged Navier-Stokes (RANS). In fact, Navier-Stokes equations are reassessed after applying the Reynolds decomposition to the flow variables that are separated into a mean and fluctuating components and then time averaged. Then, the closure approaches are distinguished depending on the number of transport equations adopted. Amongst other, one-equation and two-equations models are the most diffused techniques. As regards the former, the Spalart-Allmaras (SA) model utilizes only one transport equation for an eddy viscosity variable and was developed for the aerospace industry in order to reach satisfactory capabilities to accurately predict turbulent phenomena involving adverse pressure gradients and to recover smooth laminar-to-turbulent transition at selected positions [3]. Different versions have been created to improve specific aspects of the original model dating back to 1992: as an example, the Spalart-Allmaras with Edwards modification (SAE) has shown to provide a versatile tool for a satisfactory analysis of different applications [10]. Two-equations models are usually named according to the variables solved by the additional transport equations: among them, the two most adopted are the so-called K − ε and K − ω, solving for the kinetic turbulent energy K and, respectively, the turbulent dissipation rate ε or the specific turbulent dissipation rate ω. Despite the first being still the most widely used, it shows a degraded accuracy when flows involve adverse pressure gradients: to this end, the K − ω model, especially in its modified version from Menter (1994Menter ( , 2003 known as K − ω Menter Shear-Stress Transport (SST), provides an improved tool, combining the advantages of both the models [3]. Although superior than one-equation closures, two-equations models are still insufficient to accurately predict flow physics, most of all when strong separations occurs [10]. The phenomenon of dynamic stall has been, and still is, widely investigated by the research community that has already produced many publications comparing the turbulence models against the several experimental data available, particularly favoring the use of the K − ω SST model [19]. Although promising results have been obtained either for the 2D and the 3D cases numerically analyzed, the determination of the most appropriate turbulence model is still an open endeavor due to the aforementioned complexity in obtaining validating data from experiments under the same conditions [10,22,19].

Stochastic Tools for Turbulent Flows
The hope for the formulation of an universal theory of turbulence indissolubly lies in a deeper understanding of the physics undergoing turbulent flows [6]. What makes turbulence such a challenging and, hence, fascinating problem is its inherent complexity, which forces researchers to continuously devise new approaches for facing it. In fact, its innermost nature relies on nonlinearity, unsteadiness, rotationality and dimensionality. This explains why, although being normally approached as a stochastic issue, turbulence absolutely can not be treated with the usual simplifications adopted in statistical mechanics, for the linearization usually applied to fluid mechanics problems would lead to destroy the real essence of turbulence [20]. After being considered for many years mainly as a random phenomenon, turbulence was probably firstly attributed an organized feature in 1956 by Townsend, while at that time some experiments had already agreed upon its intermittent nature [6,24]. However, the current concept of coherent structures in turbulent phenomena became widespread only in 1971 with Brown and Roshko investigation [15]. According to Berkooz et al. (1993) [2], coherent structures in turbulent flows can be physically described as "organized spatial features which repeatedly appear (often in flows dominated by local shear) ad undergo a characteristic temporal life cycle". Clearly, a hint of order in the apparent randomness of the turbulent motion aroused suddenly a great interest in the fluid dynamics community, for it represented a first step towards a better understanding, and so handling of the phenomenon [15]. In this respect, an important contribution was given by Lumley (1967) [25] by introducing to the turbulence community the so-called Proper Orthogonal Decomposition (POD), which, henceforth, provided a mathematical shape to the concept of coherent structures as well as an unbiased method for identifying them [1,40,8].

Proper Orthogonal Decomposition
The purpose of POD is to perform a modal decomposition of a set of numerical or experimental data, by extracting a basis that is indeed optimal. Let the data set consist of N s vector-valued realizations and be defined as is usually referred to as snapshot and X = (x, t ) ∈ D = R 3 × R + ; the optimality condition reduces to find a basis {φ j (X)} +∞ j =1 in the infinite set of vector-valued functions spanned by {θ j (X)} +∞ j =1 , such that the error between the original signal and its projection onto the basis is minimum in some average sense, that is: where (·, ·), · and 〈·〉 are respectively the inner product, the induced norm and the ensemble average in the L 2 (D) space [20,12,1,40]. More, one is only concerned on the parallelism between the basis and the snapshots, so the amplitude of φ is normalized; this, bearing in mind that the statement in Equation 3.1 is equivalent to maximizing the projection, allows to reassess the issue as a classical constrained optimization problem [12,26]: Then, one can prove that the calculus of variations reduces this problem to an eigenvalue problem in the shape of a Fredholm integral equation, as follows [20,1,12]: where the kernel of the operator F , R, is the autocorrelation matrix, defined as [12,1]: and * denotes the conjugate complex.
If the integration domain D is bounded, one can evoke the Hilbert-Schmidt theory, which yields important properties of the POD [16, 12, 1, 2]: (i) Equation 3.3 has a denumerable infinity of solutions, which are represented by a discrete set of POD eigenvalues, {λ (n) } +∞ n=1 , and POD eigenfunctions or modes, {φ (n) } +∞ n=1 , and are solutions to 3.2; (ii) since R is non-negative, F is self-adjoint and non-negative, so the eigenvalues are real and positive, and then can be ordered: and their series converges: where δ mn is the Kronecker delta; (iv) the collection of eigenfunctions forms a complete orthogonal set, so the original data set can be reconstructed as follows: where the POD modes coefficients a (k) are obtained by projection of the original signal onto the basis:

Proper Orthogonal Decomposition
(v) according to Mercer's theorem, R can be decomposed as: (vi) putting together statements (iii), (iv) and (v) it is possible to prove that even the modes coefficients are uncorrelated (i.e. orthogonal) and represent, if averaged, the eigenvalues: (vii) from statements (iii) and (v) it is possible to prove that each eigenvalue is the measure of the relative contribution of the corresponding structure of order (n), and the following holds: where tr(·) denotes the trace of the matrix. If u(X) is a velocity field, then E is equal to twice the turbulent kinetic energy K of the flow, otherwise it represents just a measure of the information of the initial signal carried by the modes, and for such a reason it will be referred to as POD modes energy (though it may have no connections with energy in the mechanical sense).
When the input data arise from experiments or numerical simulations the N s snapshots are vectors which is the number of points (or nodes) where the variable is defined [20]. In particular, in the case of numerical data, or even certain kinds of experimental measurements (e.g. Particle Image Velocimetry (PIV)), the resolution of the computational domain is much greater than the number of time realizations, i.e. N p N s ; the resolution time for such eigenvalue problems of size N s can be reduced by means of the snapshot POD proposed by   [40,41,42], for which the variable X is assimilated to time and so the correlation is actually a temporal correlation [12,37]. In the present work, as often found in literature, the decomposition is applied to the fluctuating part of the data u (x, t ) that is extracted from the original signal by subtraction of the mean component u(x), which, in the cases with finite dimensions, is calculated with the arithmetic mean [9,40,20]: If the snapshots matrix U of the fluctuating components is defined as a column-wise collection of the snapshots, as follows [20]: the solution to problem 3.2 can be computed by directly applying the Singular Value Decomposition (SVD) to U ; however, as described in the next section, SPOD is based on a manipulation of the correlation matrix that arises from the original formulation of the snapshot method. In fact, this approach stems from a particular definition of the modes in terms of the original data set, that is [12]: Then, by rearranging Equations 3.3 and 3.14 it is possible to prove that the sufficient condition needed for the problem to have a solution is [12]: which is a N s -dimensional eigenvalue problem for the correlation matrix C of size N s × N s , defined as [37,44]: W being a N p × N p diagonal matrix accounting for the cell volumes of the grid (i.e. a spatial weighting located at the nodes) and a (n) , n = 1, · · · , N s are vectors of the same size N s as the time discretization, and for this reason are interpreted as coefficients purely dependent on time [12]: Then, the collection of modes φ (n) is interpreted as a purely space-dependent basis and each vector of size N p is computed according to the following [12]: where the factor before the sum is needed for the vectors to be normalized, while already being orthogonal by definition; also, statement (vi), concerning the orthogonality of the modal coefficients, is still valid and is reported below for the finite-dimensional case: As a result, the original data set can be recovered from [12]: At this point, for the sake of clarity, two aspects should be specified. The first regards the inherent linearity of the method. In fact, the decomposition is assessed as a linear procedure and according to Equation 3.20 the original nonlinear problem is reformatted by means of a basis that arises from a linear process. This, clearly, does not mean in any way that the problem has been linearized but, rather, that it has been projected onto a finite-dimensional hyperplane through the definition of the correlation matrix that is, indeed, the result of a nonlinear phenomenon, and, furthermore, this hyperplane is, on average, the closest (i.e. optimal) to the original data set [39]. Although the linear character allows to take advantage of the knowledge on the theory of linear operators for defining the properties of the decomposition, the drawback that one must bear in mind is that the optimality statement is implied only with regard to other linear representations [20]. The second remark concerns the separation of variables in Equation 3.20, which is a difference with respect to the formulation of statement (iv) for the same concept. In fact, the general treatment reported at the beginning of this section assesses the problem for X, regardless of its dependence on one or more variables, and, in the same way, even the reconstruction of the signal is made through eigenfunctions depending on that variable, while the modal coefficients are pure constants. However, in common applications of modal decompositions (e.g. Fourier series) there is a tendency to split variables and, so, to identify space-dependent modes and time-dependent coefficients, which leads to expansions in the form of 3.20. It should be stressed that there is no a priori justification for doing that, when data depend on more variables, and such an approach may not be suitable for all the cases but, rather, its adoption should be examined according to the features of the input signal and the information one expects to extract from it. In what follows, the attention is addressed to the identification of space modes and therefore the separated-variables approach is adopted [44,20].
It is important to notice that, in the discrete problem, Equation 3.3 reduces to a finite-dimensional problem on the autocorrelation matrix R: as a result, it is possible to interpret the eigenvectors of the matrix, in a geometrical sense, as the "principal axes" of the data set, thus representing the most correlated structures of the signal. This suggests the interpretation of the modes as a mathematical shape of the coherent structures [20]. Although the correlation matrix C is not indeed equal to R, the latest consideration regarding the modes is still valid for the snapshot method, since, as suggested by its derivation, the problem that is being solved is again 3.3 [39]. As a concluding remark, another useful aspect of the application of the POD is related to the possibility to build a reduced order model, by taking advantage of the basis provided by the modes. Whether its purpose is to reconstruct the original field or to provide a projection space for the governing equations in a numerical model, the important benefit given by the optimality is that the first modes are able to recover the greatest part of the initial signal. In fact, after deciding the percentage of original information to be retained for the reconstruction (δ), the number of retained modes N r can be derived from the relative information content, according to [12]: In particular, the experience suggests that very few modes are needed to recover at least the 90% of the original signal, which results in extremely reduced dimensions of the problem investigated [20].

Spectral Proper Orthogonal Decomposition
The key concept at the basis of the Spectral variation of the POD proposed by Sieber et al. (2016) [37] is the application of a low-pass filter along the diagonals of the correlation matrix C, which leads to the definition of the filtered correlation matrix S, whose elements are given by: where the filter coefficient g k is the component of a vector g of length 2N f + 1 and N f is the filter width. Then the problem is reassessed for the filtered matrix, in the same way as described above for the snapshot approach, the new eigenvalue problem being: where µ (n) and b (n) , n = 1, · · · , N s are respectively the eigenvalues and time coefficients of the SPOD [37]. The orthogonality of the mode coefficients is still valid, but now they are scaled with the new eigenvalues, such that Equation 3.19 is restated as follows: More, the energy carried by each mode is still represented by the corresponding eigenvalue, such that: Again, the spatial modes are obtained in a similar way as in the snapshot POD, but in accordance with the new variables, as follows: It is worth noting that, because of the filter action on the matrix, the new modes computed are not orthogonal, with the inevitable drawback that they are no more uncorrelated [37]. For a further description of that, the reader is referred to the paper. Then, the original data set is reconstructed accordingly as: Moreover,   [37] propose a method for coupling the modes. In fact, when periodic coherent structures are present, their behavior is described by couples of modes that show similar spectral content, in the same way as sine and cosine components compose modes of a signal decomposition through the Fourier series. Usually, the coupling procedure is made through the visualization of the phase portraits and the spatial modes plots, even though this technique has not any character of objectivity. Although the approach is not reported here, the intention is to take advantage of the Dynamic Mode Decomposition (DMD) applied to the time coefficients in order to evaluate the harmonic correlation, which is based on the cross-spectra of the coefficients, evaluated at a quarter period phase lag. The coupled modes, then, appear as maxima in the integrated cross-spectral density matrix obtained through the DMD modes and the indexes of the corresponding i -th row and j -th column define the modes constituting the coupleb, as follows: The same indexes also allow to compute the relative energy content of the pair as: The method even provides an average frequency of each modes couple, weighted with the corresponding correlation values extracted from the the cross-spectral matrix. Although the approach gives an objective procedure for the identification of the couples, a validation by means of visualization of the results is all the same necessary since, as reported by the authors, it may even lead to unphysycal pairs [37].
Since the essential aspect of the SPOD is the application of the filter, its features deserve a particular discussion. As described by Sieber et al. (2018) [35] in the generalized derivation of the method, the filter is actually a window (Gaussian) function w(τ) applied to the time series which result to be represented by a set of windows centered at the given time step t k corresponding to the snapshots, with τ being a short time scale within which the data are correlated. When the correlation matrix is computed in the discrete case, the result is Equation 3.22, where, in fact, the filter coefficients are defined as: Then, the effect of the application of the filter should be considered with respect to the quantities of the correlation matrix it modifies. According to Figure 3.1, it is clear how the diagonals contain the information of the dynamics of the signal, that is the dependence on time, while the anti-diagonals represent, for a given time t k , the effect of the time delay between the time series, which can be related to τ as follows [37]: As a result, changes of the diagonals elements refer to a variation in the amplitudes of the modes coefficients and, differently, changes of the anti-diagonals are related to their phase. This justifies that, by applying the filter as described, for it acts as a constraint for the temporal variation, the elements in the diagonals get smoothened in such a way that the diagonal similarity augments as the filter width increases; this leads to an equalization of the anti-diagonals as well, which results in a fewer distance between the two curves in the figure. Therefore, the effect of the filter on the modes coefficients is to constrain their frequency and rate of change of the amplitude, with the same consequences of a low-pass filter applied to a signal [37]. The decomposition, instead, shows a band-pass effect on the coefficients, for which the bandwidth is controlled by the cutoff frequency of the low-pass filter, given by (for a Gaussian filter): More, the center frequency of the bandwidth is not dependent on the filter characteristics, but rather is data driven, that is the filter results to "stick" to a particular frequency related to coherent structures of the data. These reasons provide a criterion for selecting the correct width of the filter, which, in fact, the authors suggest to choose with values in the same time scale as the coherent structures of interest [37]. Another important aspect of the method is that the filter operation is applied before the decomposition, therefore the whole dynamic information is still contained in the decomposition, which differs from what would happen if the filter was applied to the coefficients after the decomposition. As a result, the filter application results in a shifting of the dynamic content from one mode to other modes: this can also be verified by analyzing the eigenvalues µ, which define the energy of the modes. In fact, if compared with the eigenvalues of the non-filtered matrix, they show a lower energy content for the first modes, and contextually a higher content for the successive ones: this means that the dynamic content is extracted from these modes and transferred to less energetic ones. The advantage of this process is that new modes can be detected, which was not possible with the classical approach since they were hidden as noise in the most energetic modes [37]. The results obtained with the SPOD also depend on the way the diagonals are smoothened, that is the shape of filter used, which, in fact, can be described by any kind of finite impulse response function. In particular, Sieber et al. (2016) [37] consider two possibilities: a Gaussian filter, for a smoother response, and a box filter for a more intense effect. The former can be described by the classical Gaussian function, whose maximum should be shifted according to the chosen value for the filter width; then, the following can be considered for the discrete case: where At this point it is worth noting that the application of these filters makes it necessary to define boundary conditions for the filtered matrix. In fact, when looking at Equation 3.22, it is possible to conclude that, as the series index k varies, the indexes of the non-filtered matrix, for certain diagonals, may take values greater than the size of the matrix itself. In this context, the boundary conditions define how these values are treated; in particular   [37] refer to two possible choices: either padding them to zero or mapping them back into the domain, by adding or subtracting the number of snapshots, which is called periodic condition since in this case the time series are considered as cyclic.
In particular, when the latter is chosen together with a box filter with the same width as the number of snapshots, the filtered matrix reduces to a Toeplitz matrix, whose eigenvalues and eigenvectors are obtained by a Fourier transform of the first row. Therefore, this limiting case actually coincides with the Discrete Fourier Transform (DFT) of the data. For this reason, as stated by Sieber et al.
(2016) [37], "the SPOD allows for continuous fading from the energetically optimal POD ..." (N f = 0) "... to a purely spectral DFT" (N f = N s ), which decomposes the snapshots purely as sine and cosine combinations and this is why all the modes show the same level of correlation ( Fig. 3.2).

Numerical Model
The modal decomposition is applied to a flow field generated by a NACA 0012 airfoil undergoing a pitching motion, under deep dynamic stall conditions. The data analyzed are provided by numerical simulations performed through the Deutsches Zentrum für Luft-und Raumfahrt (DLR) TAU-code. Although the purpose of this work is not to provide a complete validation of the investigated model, but rather to inspect the modal features of the flow, nevertheless it is worth being aware of its ability to simulate the phenomenon and for this reason two simulations are run in order to compare the results with the ones from an available experimental data set.

TAU-Code
The DLR TAU-code is a 3D hybrid RANS solver composed by a collection of algorithms coded in C language [17]. The solver has been under development at the DLR for more than two decades and it represents a fundamental CFD tool for the analysis performed at the DLR itself and for the academic environment of the German universities, as well as for the aerospace industry, which has given important contributions to the validation and improvement of the code [34]. Although the code does not include a grid generation environment, it counts modules for the grid manipulation, adaption and deformation, one of them implementing the so called Chimera technique. This useful method allows to deal with moving surfaces, which are discretized through multi-block overlapping grids where the information is transferred from one block to another by interpolating between the corresponding elements on the common boundaries. The pre-processing module is responsible for partitioning the computational domain into a number of blocks equal to the number of physical processors used during the simulation and for creating a dual grid on which the Navier-Stokes equations are solved; this new grid stems from the mesh provided by the user, whose cells mid points and faces determine the new triangular facets. The high efficiency on parallel computations, then, represents an important feature of the solver, which takes advantage of the Message Passing Interface (MPI) protocol in order to solve the equations on the domains defined at the previous step. The several computational schemes and turbulence models available are given as input for the code, together with the other parameters of the simulation, through an ASCII-file compiled by the user [34,14].

Computational Domain
The chosen airfoil is a NACA 0012 characterized by a chord length c of 0.15 m and subject to a pitching motion about its quarter chord location ( Fig. 4.1). The profile is described by a symmetric closed curve with the maximum thickness located at 30%c and equal to 12%c; among many other, the configuration adopted here presents a blunt shape at the trailing edge.  The computational grid for the fluid domain surrounding the present airfoil is generated with Pointwise. In order to perform the unsteadiness of the body, it has been chosen to adopt the Chimera technique, according to which a multi-block approach is necessary. In particular, two different blocks lying on the x-z plane are created: a inner block for the region in the neighborhood of the airfoil, and a outer block for the farfield, which acts as a fixed reference frame for the rotation of the other block (Fig 4.2). The former is built as a O-grid block, that consists on a circular region spanning a radius equal to 1.74c, centered at the origin of the frame that coincides with the quarter chord point (Fig. 4.2b). Then, the whole block is subject to a pure rotation about the y axis passing through that point, which results to be the center of rotation of the pitching motion of the airfoil, described by Equation 2.3. Whereas the second block is built as O-grid of radius 100c centered at the origin of the frame, with a circular hole corresponding to the location of the inner block ( Fig. 4.2a).
As concerns the kind of elements used, it has been chosen to adopt quadrilateral elements for the inner block, where the mesh is extruded in the direction normal to the airfoil surface for the first 40 layers and then extended to the circular boundary with a structured grid. The airfoil outline is discretized with 690 steps, with an increased density at the leading and trailing edge, symmetrically distributed over the two sides, while the spacing of the first layer from the surface is calculated according to the estimate of the dimensionless wall distance y + = 0.65 together with the flow conditions of the case study. In particular the value thus obtained is equal to 1.5 · 10 −6 mm or, alternatively, to 1 · 10 −5 in dimensionless units and would correspond to an estimate of y + < 0.1 if the flow conditions of the considered experiments are taken into account. Then, the grid is further extruded radially for 5 more layers in order to generate the overset region, emphasized in Figure 4.3b as a red annulus. In fact, this latter marks the starting grid of the outer fixed block: the remaining part of its mesh is then generated mainly as an unstructured grid, except for a more refined region behind the airfoil which extends for more than 6c in order to better "capture" the wake behavior ( Fig.4.3b).
The resulting mesh counts, respectively: 102512 points and 101824 cells for the inner block and 103282 points and 126594 cells for the outer block. Although the grid is not assembled with Pointwise, but rather with a dedicated function of the TAU-code, Figure 4.3 allows to understand the actual domain over which the flow equations are solved. Furthermore, the domain has been so far described with a 2D character. In fact, the solver requires a 3D domain, so the mesh described above is extruded in the y axis such that the depth in this direction is of a single cell. Then, if the two parallel planes separated by this distance present a constant offset for the corresponding nodes and are given a symmetry boundary condition, the grid is reduced to a real 2D domain and the flow equations are solved accordingly [14]. This consideration justifies the bi-dimensional treatment of the results that follow.

Comparison with Experiments
Two different numerical simulations are run with TAU and the corresponding results are compared with the experimental data produced under the same conditions by Gerontakos (2004)   where the reduced frequency k is given by: According to the concepts outlined above, such conditions are consistent with the onset of deep dynamic stall. Furthermore, the referred Mach number suggests that no compressible effects should take place throughout the simulation. The starting condition of the simulations is provided by a converged steady case previously solved under a convergence criterion on the relative error for the density ρ, which has to be lower than 1 · 10 −5 in order for the iterations to be stopped. Then, the unsteady simulations are run for 4 complete pitching cycles, each of which is discretized into 3600 time steps, i.e. a dimensionless time step of 2.78 · 10 −4 ; each time step, in turn, is performed for 300 inner iterations, which provides a reasonable trade-off between the convergence reached and the calculation time required. The forces and moment coefficients are compared below. Whereas the experimental curves are

Comparison with Experiments
directly taken from the reference publication [18], the numerical data reported result from the averaging of the last 3 cycles since they show barely distinguishable differences among each other, while being significantly different from the first. The results in Figure 4.4 show a very similar behavior between the two models considered, and even a good agreement with experiments, when the upstroke part of the motion is considered.    However, the two curves separate at about 12°and from this point the SAE trend shows higher values than the other which, instead, is closer to the experiments. Both the models anticipate the lift stall condition, which occurs at α = 24.3°according to the experimental data but is predicted at α = 22.34°and α = 22.55°respectively by the SAE and the SST. This behavior can be also noticed in the other two plots. In fact, the beginning of the moment stall ( Fig. 4.6) is clearly anticipated by the CFD (17.5°from the SAE and 18.6°from the SST, compared to 20.26°from measurements), while the drag coefficient curves (Fig. 4.5) separate at around 13°, showing, hereafter, considerably amplified values with respect to the experiments. Furthermore, the two models show substantial discrepancies in the prediction of the loads peaks. While this difference is not so marked in the lift curve, for which the SAE provides a closer value to experiments than the SST, it is more clear when drag and moment coefficients plots are considered. In particular, for these results the SST predicts values that are not as overestimated as the ones given by the SAE and for the moment coefficient both the models show very marked differences. After the lift stall occurs, the formation of a second vortex is suggested by the presence of a second peak shown by both the models. In particular the flow results show that it forms at the trailing edge and it is therefore called Trailing Edge Vortex (TEV). Although this aspect is by no means detectable in the experimental curves, the two model show a similar behavior in this respect. In fact, the three plots show this sudden change of trend occurring at very similar angles, whereas the load values are comparable only for the lift plot. The drag and especially the moment coefficients curves present a second peak considerably lower for the SAE than for the SST and, in particular, this is the condition of the maximum (negative) moment acting on the airfoil for both the cases. The beginning of the downstroke phase of the SAE is characterized by another strong peak, such that it defines the maximum values reached by either the lift and the drag plots. Although present also in the SST case, this event produces barely considerable values and for this reason this condition marks a first distinguishable element of difference between the two models. In the experimental data, instead, the occurrence of this structure only determines a smooth change of the curves trend around 21.8°. Then, the other part of the downstroke phase is characterized by the continuous succession of LEV and TEV formation and shedding, evidenced by the several peaks in the plots of both the models. In this context, the SST model shows a smoother behavior, while the SAE model predicts another great peak in the lift coefficient curve at an angle of attack of 14.44°w hich is totally missing in either the experimental results and the the other model. This behavior is also the cause of an increased value of the area enclosed by the moment coefficient curve such that the associated negative damping results to be much higher than the positive one. Another great difference from the experimental data can be seen in the phase facing negative angles of attack: this part of the motion clearly shows another hysteresis when the angle reaches its minimum value and then starts to increase again. This behavior is totally lacking in the CFD results which show, for this phase, a very similar plot. This evident difference, though, is not noticeable from the other coefficients plots, which instead show a fair agreement with the experiments. The results presented show that, overall, the SST model recover a solution that better approximates the experimental results, since the discrepancies are clearly smaller than for the other model. More, the results recovered here by this model seem to agree with the ones presented for the same case by other authors in literature [45,21,27], in terms of either the behavior and the values of the plots, with a clear tendency to anticipate the onset of the lift stall and general inability to predict the exact values of the loads, especially for the phase of the flow under fully separated conditions. For the SAE model, instead, no results have been found yet in literature for the same case. For such a reason, and considering the results presented above, the SST model, in this context, is considered to be the one that better predicts the behavior of the dynamic stall, if compared with the other model. This consideration has led to the choice of this model for the simulation of the case study, whose results are reported in the next section.

Case Study Results
The grid described above is used together with the SST model in order to produce the results that are analyzed with the SPOD. The data for the simulation are provided by the Institute of Helicopter Technology of the Technical University of Munich, and since no experimental data are yet available the following results are not validated but only discussed in terms of numerical data. The configuration parameters are summarized in Table 4 Even for this case, 4 complete pitching cycles are performed and the time discretization, as well as the convergence criteria, are the same reported in the previous section. The results clearly show the typical behavior detectable during dynamic stall phenomena, with distinct hysteresis in the loads coefficients plots (Fig. 4.7 The coefficients plots present circled numbers that refer to the corresponding flow conditions in Figure 4.10, which depicts, throughout one pitching cycle, some instantaneous contours of the curl field (y component) normalized as c U ∞ and related flow streamlines.
The upstroke phase starts with the linear behavior 1 of the lift curve ( Fig.4.7) to which corresponds the condition of fully attached flow (Fig. 4.10a). As a result, even the trend of the other loads plots (Fig. 4.8, 4.9) follows the trace of the typical results obtainable from the steady simulations. The flow evolves with the same features, even though a small laminar separation bubble arises close to the leading edge at an angle of attack of about 12°and starts to grow hereafter. This process clearly produces the thickening of the boundary layer in the suction side but, however, it does not even detach from the surface when the static stall angle of attack is approached at 14.94° (  Fig. 4.10b). Near this condition, a slight change of the trend can be seen in the lift plot, while the drag curve does not show to be affected appreciably. As concerns the moment plot, a more pronounced negative trace can be detected exactly at this angle of attack, since the core of the high vorticity region has just passed the center of pressure, thus inducing a negative moment on the airfoil. The beginning of the moment stall 2 , approximately at 16.1°, marks a sudden change in the trend of the moment coefficient curve which starts to drop rapidly. Contextually, a change of the slope can be noticed even in the lift plot, whereas the drag coefficient still shows a smooth behavior. As Figure 4.10c suggests, this event is again connected with the motion and the suction effect of the separation bubble. In particular, at this point it can be clearly recognized as the Leading Edge Vortex (LEV) and as its core passes well beyond the quarter chord location a stronger nose-down pitching moment is exerted onto the airfoil, which justifies the steep decay of c M . The LEV grows as it sheds rearward, causing the continuous drop of the moment curve and, concurrently, the increase of the lift force due to the augmented suction effect performed by the low pressure region within the vortex. Exactly before the LEV sheds into the wake, the lift curve reaches its peak value 3 , which defines the occurrence of the lift stall at α = 19.7° (Fig. 4.10d). This event is also the cause of an inversion of the slope in the moment curve, and even the gradient of the drag plot is  clearly reduced although the maximum drag is reached only slightly later at α = 20.3°. During this phase the LEV reaches its maximum size of about 90% c and, as a result, the suction side appears to be completely covered by a circulation region. However, the streamlines still evidently reattach to the surface, until the full detachment condition is reached and the typical trend of the stalled flow can be seen in the lift curve. As the vortex detaches, suddenly a counter-clockwise circulation region arises at the trailing edge at about 20.1°and grows hereafter. The induced suction effect due to low pressure region within this bubble, produces another negative moment on the airfoil and, consequently, the slope of c M becomes negative again at α = 20.7°and, contextually, the drag coefficient shows a smooth rising while the lift curve does not show to be remarkably affected. The growth process continues and leads the moment to reach its maximum (negative) value 4 at α = 21.23° (Fig. 4.10e). At this point, the circulation region can be clearly distinguished as the Trailing Edge Vortex (TEV) which has also the effect of amplifying the negative trend of the drag curve. At the same time, the rise of another smaller counter-clockwise circulation region can be detected just before the quarter chord: its suction effect, then, induces a positive moment on the airfoil, which seems to be the reason why the moment curve starts to rise hereafter. Either this negative curl bubble and the TEV continue to grow attached to suction side while a clockwise circulation region clearly starts to be enclosed between them (Fig. 4.10f). However, during this phase no particular changes seem to be exerted by the flow onto the airfoil in terms of aerodynamic loads. As soon as the TEV leaves the airfoil at about 22.55°, the positive curl region, which is indeed a second vortex, suddenly attaches to the airfoil in such a way that its core is located beyond the quarter chord ( Fig. 4.10g). Meanwhile, the circulation bubble near the front is still attached and at this point has reached its maximum size, after moving rearward past the quarter chord. As a result, all the three loads clearly experience an abrupt change of slope in correspondence of this event. Then, their trend remains unchanged even after the downstroke phase has begun. In fact, the second vortex sheds into the wake 5 only after the maximum angle of the motion (22.73°) is reached: in particular, this happens at the angle of attack of 22.71° (Fig. 4.10h). Consequently, as for the primary vortex, the three loads plots show again a sudden change of trend. However, it should be emphasized that this second vortex produces effects on the loads that are significantly lower than the ones related to the primary. This event marks the beginning of a fully detached flow condition that causes the mainly monotonic trend of the curves of the aerodynamic coefficients, except for some smoother variations. This wavy pattern of the downstroke phase is to be attributed to the continuous succession of clockwise and counter-clockwise circulation regions that form respectively at the front and at the back and attach to the surface for a short period before shedding into the wake, as suggested by Figure 4.10i and Figure 4.10j. The reattachment process starts from the leading edge at about 15°and moves rearward. In particular, from the angle of attack of 10.83°, whose flow behavior is depicted in Figure 4.10k, the coefficients plots show an almost constant trend 6 until the end of the downstroke. Then, although the flow is mainly attached (Fig. 4.10l), it recovers the same conditions as the steady case only after the following cycle has begun and it repeats with the discussed features.

Modal Analysis of Simulation Data
In this chapter, the Spectral Proper Orthogonal Decomposition (SPOD) is applied to the numerical results described in the last section of the previous chapter. The decomposition is performed for two different cases in order to show the advantages and disadvantages of each one; in particular the comparison concerns the modal analysis of a vector-valued field and a scalar-valued field, respectively the velocity field and the pressure field. The data used for the study refer to the inner pitching block. The reason of this choice relies on the time saving due to the handling a lower amount of data related to a smaller domain. The analysis was also tested on the whole domain and the results showed that the flow dynamics is distributed in a higher number of modes, even though the main features do not seem to be highly affected. In order to avoid the influence of the wake dynamics on the structures connected to the wall effects, Wen et al (2018) [47] selected a confined region close to the airfoil, from which they extracted the data for the decomposition. This consideration supports the choice made in the present work. The snapshots are taken every 36 time steps within a cycle which is discretized into 3600 time intervals. As a result, 100 outputs are produced for every cycle and the analysis is performed over the last 2 complete cycles, which means 200 snapshots in all. The values of the studied variables are located at the nodes of the referred region. In particular, since the block rotates at every time step, its coordinates in the geodesic reference frame are different from one time interval to another while the coordinates of the body-fixed frame are inertial with respect to the pitching motion and, therefore, unchanged throughout the simulation. In order for the decomposition to use a pure time correlation, all the snapshots values must be described in the same reference frame and for this reason the variables are extracted with respect to the body-fixed frame coordinates.

Code Implementation
A Python code, based on Sieber et al (2016) [37], was independently written expressly for the analysis of the investigated case and was validated against the Matlab code and corresponding data made available online by the paper authors. The routine applies the SPOD to data sets arising from 2D results of a pitching NACA 0012, which have to be stored as arrays for the coordinates of the domain discretization and the variables one wants to investigate, be they scalar-valued or vector-valued. Then, the specific parameters of the motion have to be set in the first lines of the code, in order for it to properly scale the quantities used, according to the case. During the run, the user is asked to input the kind of field that is to be analyzed and, subsequently, the value of the relative energy content δ for the reconstruction of the signal and the filter width N f to be used for the analysis. The code, then, is responsible for printing all the plots that will be discussed in the next sections and for storing them in a folder specified by the user; in the event that the pressure field is chosen for the study, the routine also produces the plots of the estimate of the aerodynamic coefficients (c l , c d , c M ) against the angle of attack during a cycle, which are only based on the pressure contribution since no information is given on the skin friction coefficient c f from the single pressure field. The kind of boundary conditions is chosen by the user between periodic and zero-padded. The code adopts a Gaussian filter with a constant K equal to 3 (see Equation 3.33), thus given as: However, if the filter width is set equal to the number of snapshots, the code automatically applies a box filter with periodic boundary conditions in order to recover a DFT of the data set.
In what follows, the SPOD results are obtained with periodic boundary conditions. As additional note, the ordinal indexes of the modes appear as subscripts instead of as superscripts, as reported in Chapter 3, for the sake of a more compact representation.

Vector-Valued Variable Decomposition: Velocity Field
In this work, as commonly found in literature, the decomposition is applied to the data set deprived of the time-averaged component, that is to the fluctuations of the variable throughout the time evolution of the phenomenon. However, by analyzing the mean flow it is possible to investigate the dominant behavior of the variable and, in fact, it can be thought as the 0-th mode of the POD when the original field (without subtraction of the mean component) is decomposed. In Figure 5.1a it is possible to notice how the velocity time average shows a clear deficit in a wide region close to the suction side, while the remaining part of the domain is characterized by magnitudes more similar to the free stream velocity, except for the zone of the stagnation point at the leading edge. In addition, the streamlines trend over the upper region depicts a curvature more marked than that one would expect by looking at the shape and orientation of the airfoil. This leads to conclude that this region is dominated by thick boundary layer flow conditions and it is confirmed by the description of the flow evolution in the previous chapter. The spectral content of the flow can be analyzed through the spatiallyaveraged Power Spectral Density (PSD) of the velocity field. Figure 5.1b depicts the frequency response obtained with the Welch's method, expressed in terms of the Strouhal number S t based on the boundary layer maximum thickness which is estimated as the projection of the chord in the vertical axis when the angle of attack is maximum within the pitching cycle. Then, the following holds: where f is the frequency. The PSD plot, although characterized by many strong fluctuations, probably due to the graphical effect of a reduced size of the time domain, suggests that a strong and narrow peak is present at S t = 0.017, which corresponds to the frequency of the pitching motion. Two weaker peaks, featuring a much broader band, can be identified at S t = 0.083 and S t = 0.182 that are respectively equal to 5 and 11 times the main frequency. Furthermore, if the DFT is performed on the data set, the energy distribution of the modes pairs in the frequency domain reflects the PSD curve, as already evidenced, for different case studies, by Sieber et al (2016) [37]. In order to investigate the structures of the energetically optimal modal basis, the fluctuating component of the flow is then decomposed by means of POD, which is recovered by setting the SPOD filter width as N f = 0. The correlation matrix is computed according to  In Figure 5.2 the first 50 modes are represented in terms of the energetic content of their eigenvalues, either for single modes (blue) and cumulative contribution (red). The number of retained modes for the coupling is chosen such that the cumulative energy content is at least 99%, which represents a reasonable trade-off between the intention of describing the most significant part of the flow and the aim of rejecting the smallest contribution, which may affect the modes as a noise source. In the case of POD, the condition is guaranteed by the first 16 modes. Then, 8 couples are identified and plotted with respect to their modal energy against the Strouhal number in Figure  5 As the plot suggests, the POD is not able to extract the modes pair referring to the 11-th harmonic, as predicted in Figure 5.1b: this action is instead performed by the SPOD, as will be seen. The couple connected to the first harmonic (1) describes, by itself, 65% of the initial data set but it is not the one with the best correlation that is instead shown by the second pair (2), which contains 21% of the total modes energy. The third couple (3), featuring a relative energy of 6%, is chosen together with the other two already mentioned for in the SPOD they will constitute the most correlated pairs of the three main frequencies derived from the spatially-averaged PSD. Their coherent structures are described below. In the same way the time coefficients are linked, the corresponding spatial modes are constituted by a real and an imaginary component and they are both reported. The correlation of the modes can be seen through the most defined structures and their periodicity in the space domain. Clear correlations in time usually go in hand with as much distinguishable correlations in space, which depict structures of similar dimensions, shifted in the traveling direction of a fraction of the wavelength. The PSD of the coupled time coefficients for each pair is plotted in order to show the response in the whole frequency domain and to determine the influence of other harmonics than the dominant one, which can be done by comparing the average frequency in Figure 5.3 with the peak frequency of the corresponding PSD in Figure 5 The first pair (modes 4 and 1) shows a distinct peak at S t = 0.017 (Fig 5.4f), i.e. the flow dominant harmonic, even though the comparison with the average frequency S t = 0.026 suggests that the influence of the higher harmonics cannot be neglected. The strong peak in the PSD is due to the predominant time response of the time coefficients of the first mode (Fig. 5.5a) which makes the effect of the fourth mode visible only from the second weaker peak. In fact, by analyzing the time evolution of the two coefficients, it may be inferred that they do not represent the two components of the first harmonic; however, they both show their main peak at the same frequency and this may justify why the coupling method correlates them. In this regard, one may say that the pair created is not physical, which is an eventuality already considered by the authors of the reference paper. This is also justified by the inspection of the spatial modes. The most energetic single mode, ψ 1 (Fig. 5.4b), clearly shows a broad spatial distribution of the structures, which makes it difficult to identify characteristic dimensions and wavelength and, then, the connection with the real component of the pair. The most defined regions are located near the leading and the trailing edge, though the first ones evidently show different dimensions. The fourth mode (Fig. 5.4a), instead, presents more defined coherent structures with a clearer periodicity in the spatial domain behind the airfoil. However, the structure at the leading edge and the ones with blue contour still show a wide distribution which reduces the distinction. As regards the second pair (modes 2 and 3), its frequency response (Fig. 5.4f) depicts a broader band around the peak at S t = 0.083, i.e. the fifth harmonic. The influence of other harmonics can be even inferred by considering the average frequency of the couple, S t = 0.136; still, both the time responses of the coupled coefficients show a clear similarity, as the Lissajous figure (Fig. 5.5b) suggests: in fact, the spiral trend reflects the similitude even at different frequencies which justifies the high correlation obtained. In terms of spatial eigenfunctions, defined coherent structures can be identified in the wake, starting from a limited region over the trailing edge. The dimensions and wavelength are more evident than in the previous pair, but the broadness of the spatial distribution is still clearly present, most of all for the real part of the couple over the suction side (Fig 5.4d). The third single mode (Fig. 5.4e) shows a similar sequence of the structures as ψ 2 , but shifted in the streamwise direction of a size comparable with half the wavelength of the pair, which means that the modes are fairly correlated together even in the space domain. The presence of background noise can be seen even in ψ 3 contour plot and, in particular, attention is drawn towards the two weak (blue) structures over the mid chord location, since their behavior with respect to the filter size allows to better understand the action of the filter itself. The third pair (modes 5 and 6) shows a PSD peak (Fig. 5.4i) corresponding to the second harmonic, i.e. S t = 0.033, featuring narrower band than the second pair; however, a second substantial peak can be detected around S t = 0.18 which suggests the influence of the eleventh harmonic. As a result, the average frequency is greater than that related to the first peak: in particular, it is equal to S t = 0.075. The time response of the coefficients is dominated by the values of the real component of the pair, which causes such an irregular shape of the phase portrait in Figure 5.5c, while the frequency distribution of both shows a similar trend. The spatial correlation can be seen in the shifted configuration of the coherent structures in the imaginary component with respect to the real one. As regards the definition of the structures, ψ 5 (Fig. 5.4g) shows fairly clear coherence with a barely viewable background noise, while the peak values show a broadness that makes them less distinct and covers all the suction side. The imaginary component shows instead a greater influence of the background noise that occupies the main part of the domain. Although the structures are still distinguishable, the shape is more irregular than in the other mode: this is clearer when one considers the valley which, in addition, seems to be modified by another peak value near the trailing edge, resulting in a weaker spatial periodicity.
A filter size N f = 30 is applied in order to perform the SPOD. The chosen value reflects the time needed for the main circulation flow structure (interpreted as constituted by the LEV and its counterpart, the TEV, together as one) to arise and travel out of the considered domain. It should be stressed that such a value is close to 20, i.e. the time scale connected with the frequency of the fifth harmonic; the results were obtained and analyzed even for that filter width, but they are not reported here. In fact, although not showing marked differences with the ones discussed below, they showed a better description of the coherent structures related to the modes 3 and 4, which represent indeed that harmonic; for the other two pairs, instead, the higher broadness and background noise influence made it harder to identify clear spatial coherence, which led to the choice of discussing the results of the mentioned filter. As suggested by Figure 5.6, the mode energy is redistributed among the less energetic modes. As a result, the distance between the single mode energies (blue stems) is reduced and the cumulative energy curve shows a smoother slope than in the POD case. The main consequence can be seen in modes from 3 to 10 that now appear as couples according to the mutual similar energetic contents and the first mode experiences a reduction up to 10% than the value obtained with the POD. Accordingly, the number of modes required for a reconstruction of the 99% of the initial data set is now increased up to 26 and they are coupled as previously described. The same six modes analyzed for the POD are considered. The three pairs show in this case an increased similarity of the distribution on the frequency domain with that shown by the spatiallyaveraged PSD. The filter operation clearly performs a reduction of the influence of different harmonics on the couples that show an augmented correlation as well as an average frequency closer to the peak one. In particular, these considerations are evidently observable for the third pair which is now located at S t = 0.023 and contains 8% of the original data. The modes energy of the other two couples is reduced, more clearly for the second one which now features the 12%, and only slightly for the first that carries the 63%, which is only smoothly affected since the loss of energy of the first single mode is partially contained by the increase of the eigenvalue connected with the second single mode. In fact, as can be seen from Figure 5.8, the first couple is now constituted by the two first single modes which causes a clear change in the trend of the frequency response of the time coefficients (Fig. 5.8c). The PSD shows augmented response values and the band is much narrower, which suggests a much lower influence of the other harmonics; in addition, the trend continues in a smoother way and the previous second peak, as well as other peaks, is completely missing. In fact, the analysis of the phase portrait ( Fig. 5.9a) indicates that the time coefficients can be nearly identified as sine and cosine components, the distortion being caused by the difference of magnitude of the two series (b 1 shows greater values in the whole set) and the evident, yet weak, oscillation at a higher frequency shown by the second mode. The coherent structures of ψ 1 (Fig.  5.8a) show a high similarity with the ones obtained with the POD, which suggests that, although a considerable amount of energy has been extracted, still the distribution of the spatial correlation is not appreciably affected by the application of the filter. The correlation with the imaginary component (Fig. 5.8b) is still unidentifiable. In fact, although the time evolution represents structures that develop throughout the time interval with similar dynamics, however the flow features described clearly refer to different phenomena, which suggests that in this case they are not representative of the real and imaginary components of traveling structures but rather different phenomena occur- ring with a similar time periodicity. This consideration is also supported by the evident difference of energetic levels of the corresponding eigenvalues. A further analysis on the different filter widths ranging from 0 to 30 showed that the red structure at the trailing edge of ψ 2 is gradually recovered from the the small fraction on the top of the red region of the corresponding POD mode in Figure  5.4d: the remaining part of that originally single structure is the the other red zone beyond the blue valley in the wake. In this case the SPOD shows an important result. In fact, the two modes evolving similarly in time, though phase shifted, but showing clear differences in the spatial correlations, can be directly connected to the primary and secondary vortexes, respectively represented by ψ 1 and ψ 2 , which are not indeed periodic coherent structures. Even for the second pair (modes 4 and 3) the PSD (Fig. 5.8f)) shows a reduced bandwidth in the surrounding of a more distinct peak, at the dominant frequency of the fifth harmonic. In this context, the action of the filter produces a drastic reduction of the components at higher frequencies that showed a weaker but not negligible effect with the POD, but still it does not affect the bandwidth components to the point of modifying the dynamics underlying its frequencies. Again, the phase portrait ( Fig. 5.9b) reflects the presence of different frequencies in the time evolution of the time coefficients, while being better correlated especially in terms of the main frequency. The coherent structures depict a high influence towards the filter application. In fact, if compared with the POD, the distribution of the spatial coherence of mode 4 ( Fig. 5.8d) displays now reduced dimensions of the red regions in favor of augmented sizes of the blue valleys: in particular, that one present in the studied region shows the rise of a distortion which produces an enlargement of the structure. The background noise can be clearly identified and this holds even for the third mode (Fig. 5.8e).
In particular, this latter evidently depicts the rise of a new coherent structure in the suction side, approximately at the mid chord: this result clearly represents the action of filter, for this region could only be interpreted as a noise of the first structure of the wake. This can be seen either in a reduced noise around this structure and weaker connection with the "parent" correlated region; in addition, the structure at the trailing edge shows an increased coherence according to the reduced number of contour levels inside it. The spatial periodicity can be inferred from the shifting of the structures in the wake direction from ψ 4 to ψ 3 , which supports the similarity already discussed in terms of time coefficients. The effects on the third pair are even more marked. In fact, in addition to the reduction of the bandwidth of the frequency response, the filter action completely changes the peak frequency that is now centered at the eleventh harmonic value, while the peak frequency recovered from the POD (at the second harmonic) is still present but featuring a lower response (Fig. 5.8i). A better correlation of the time evolution is also depicted by the phase portrait in Figure 5.9c, which now describes a more regular trend, with characters similar to the second pair. The spatial coherence is clearly outlined by both the single modes, which also feature an almost negligible impact of the background noise, thus letting the highest correlated regions emerge and be identified clearly. The coherent structures show a defined size with a reduced presence of wide boundaries, especially for the real component (Fig. 5.8g). The length scales, smaller than in the other modes, can be precisely identified, which leads to easily detect the relation of the distance with the wavelength in the shifted configuration. The action of the SPOD allows for a better matching of the coherent structures with the corresponding physical phenomena, which is hardly applicable to the POD results. By virtue of the pitching motion, the wake described in the body-fixed reference frame experiences a change of orientation, according to which the features of the flow connected to the major angles of the motion follow a path more inclined than the other. As a result, the first pair shows a distribution of the spatial correlation (more evident from the imaginary part) directed towards a steeper direction if compared with the structures in the other two couples that are, instead, slightly parted from the orientation of the chord. This fact, supported by the PSDs that do not show clear narrow peaks, suggests that pairs 2 and 3 are not representative of higher harmonics of the main vortex structure, but rather distinct flow features. In additional support of this consideration, the flow field inspection showed that the influence of the main vortex affects only the second vortex. In conclusion, the other two pairs represent the fluctuations of the flow connected with the evolution under full separation conditions; in particular, the second couple carries the information of the several eddies developing and shedding into the wake during the concluding part of this phase and it appears partially influenced by the second vortex according to the new structure extracted by the SPOD in ψ 3 . The third pair, instead, on the basis of the similar dimensions of the structures, but reduced wavelength, with respect to the previous couple, can be interpreted as the alterations of a higher harmonic of the latter. A further analysis is conducted by applying a filter characterized by a width N f = 63 which is equal to the time spent by the flow under fully separated conditions that produce the highest fluctuations due to the continuous formation and shedding of circulation zones. The modes energy redistribution (Fig. 5.10) performed by the filter produces again a strong effect on the first eigenvalue; however, the content added to the second eigenvalue almost equals the quantity extracted from the first. As a result, the introduction of energy into the following modes is done at the expense of the energy of the modes from 3 to 6 which feature reduced energies while showing, in twos, similar contents. For this case, the condition on the relative energy content equal to 99% is achieved by retaining 34 modes. The spectrum of the pairs (Fig. 5.11) shows the introduction of three new highly correlated couples, featuring better harmonic coherence than the first one. The coupled modes previously identified are still detectable located at the main frequencies depicted by the spatiallyaveraged PSD, experiencing a less marked reduction of total energy, respectively decreased to 61%, 8% and 6%. The first pair does not show substantial changes towards the increased filter width. The spectral content is further strengthened around the dominant peak and contextually reduced for the higher harmonics, resulting in a narrower band (Fig. 5.12c). In support of this behavior, the phase portrait of the coefficients (Fig. 5.13a) depicts a nearly circular shape which suggests the high level of correlation between the time series. In addition, the modes (Fig. 5.12a and 5.12b) reflect a barely distinguishable dependence on the filter, with slight influences on the outer boundaries of the structures. More distinct effects can be instead recognized in the second pair. The first consequence can be seen in the frequency response (Fig. 5.12f), where the not negligible contribution of the second peak viewable in the previous analysis (around S t = 0.165) is here clearly suppressed, which depicts a trend more typical of mono-frequent modes. In this regard, the action of the filter suggests a fading of the decomposition towards a behavior more representative of the purely spectral study, with the resultant loss of energetic optimality character. As a result, the real component of the pair (Fig. 5.12d) depicts an evident reduction of regularity of the structures shapes, while the background noise shows a reduced influence. The spatial correlation appears more fragmented with an associated decreased definition of the wavelength. The imaginary component instead shows to be only slightly affected, with a barely perceivable reduction of the background noise and clearer enhancement of the correlation near the trailing edge.  As regards the third pair, the adoption of a higher filter width provides a stronger frequency response at the second peak shown by the previous analysis, now located at S t = 0.165, while the response at the first peak is clearly reduced (Fig. 5.12i). Furthermore, the increased regularity of the phase portrait suggests an enhanced similarity of the time evolution of the coefficients series ( Fig. 5.13c). The broad band of the frequency response, still showing a peak, supports the conclusions drawn on the previous analysis regarding the physical meaning of this pair, which are further enforced by the peak frequency value, now equal to the second harmonic of the dominant frequency of the second pair. The effect on the eigenfunctions can be seen in a drastic reduction of the fragmentation of the modes, especially for the structure over the suction side either for ψ 6 ( Fig. 5.12g) and ψ 5 (Fig. 5.12h). In particular, the former depicts an evident suppression of the background noise, while some distortion is introduced into the border of the red structure. The latter, instead, shows a higher influence of the background noise than in the previous study. However, the definition of the structures is still preserved and the spatial periodicity can be still clearly recognized through the wavelength.

Scalar-Valued Variable Decomposition: Pressure Field
In this case the analysis is conducted in the same pattern as the previous. The first results discussed are the mean field and the spatially-averaged PSD (Fig. 5.14) from which it is possible to infer the general features of the flow. The spatial distribution of the time-averaged pressure field depicts a clear widely spread pressure deficit in the region over the suction side, where the lowest values are close to the airfoil surface and affect the most of it in the chord direction up to the trailing edge. The blue region shows an irregular shape with the maximum thickness located beyond the midchord point, and spreading up to the wake, which suggests the high influence of the continuous presence of circulation zones encompassing low-pressure structures. The pressure distribution on the pressure side is, instead, more typical of stationary flow conditions, with high values near the stagnation point and contour levels distributed in a more regular shape. However, the trailing edge region shows the influence of low-pressure structures, which can be related to the presence of the trailing edge vortexes affecting this part of the domain due to the rotation of the grid. As regards the spectral response, this field shows a smoother trend. In fact, despite the good agreement with the previous field upon the dominant peak at S t = 0.017, no other distinguishable responses can be detected except for slight changes of trend at the 5 th , 10 th , 16 th and 22 nd harmonics. As a result the trend of the whole response can be fitted with a linear approximation in the semi-logarithmic scale. This suggests that there are not significant responses of the pressure field at the highest frequencies; however, the decomposition is able to find pairs characterized by a clear correlation, even though their energetic content is much lower than that of the dominant harmonic. The energy plot in this case (Fig. 5.15) depicts a lower prevalence of the first eigenvalue on the other, if compared with the velocity decomposition. In addition, the second and the third modes show a higher energy content, though at much different levels than in the previous section. In general, the first modes show an enhanced energetic content; hence, 14 modes are required to reconstruct the 99% of the original data set, i.e. 2 modes less than in the previous case. Then, seven pairs are identified by the coupling method and they are reported in the spectrum plot ( Fig.  5.16). As can be seen, the reduced number of dots makes it difficult to recognize any connection with the spatially-averaged PSD, although the linear trend is almost recovered. The three modes chosen for the coherent structures analysis are the ones showing the highest energetic content and, in this case, they also feature distinct correlations of the time coefficients. In particular, the first pair contains, by itself, more than 75% of the total energy, while the second and the third are respectively representative of the 14% and 5%. In terms of spectral response, although showing almost regular spacing that may be interpreted as related to the number of the harmonic, however they cannot be related to any characteristic frequency of the flow. In fact, as will be seen below, their average frequencies show high influences of different responses. The dynamics of the first pair (Fig. 5.17c) shows a narrow peak at S t = 0.017, which suggests its connection with the main vortex shedding behavior. However, the high influence of the subsequent peaks cannot be neglected for it leads to an average value of S t = 0.030. This is also supported by the inspection of the time evolution of the two time series (Fig. 5.18a): in fact, although the first mode shows clearly stronger responses, mainly centered at the dominant frequency, the second mode has a dynamics split into several harmonics, which causes the higher value of the average frequency of the pair. The spatial distribution of the eigenfunctions shows, in this case, a much enhanced definition if compared with the previous analysis. In fact, although clearly affected by broadness, the coherent structures are better identifiable and their location can be inferred more distinctly, as well as the wavelength of the spatial periodicity. The behavior of the first two eigenmodes is similar to that shown by the previous analysis. In fact, even here, the pair predicts a time correlation that does not go in hand with the spatial correlation. The first single mode (Fig. 5.17a) reports only one coherent structure, as well as a high impact of the background noise; in this case that region is connected with the main vortex shedding, but also features the presence of the second vortex that is not recognized as a distinct structure. The imaginary component (Fig. 5.17b) is instead connected with the main vortex formation (blue region) and subsequent wake convection (red region) and no clear effects regarding the second vortex can be detected. This analysis, if compared with the velocity field decomposition, allows for a better identification of the connection between the main vortex shedding and its corresponding mode. In fact, the spatial distribution and the spectral content can be easily connected to the frequency and shedding characters of the vortex that originates upstream and travels downstream, leaving the airfoil while, contextually, the trailing edge vortex strengthens and sheds into the wake later.
The second pair shows no dominant frequencies and an evident broad bandwidth (Fig. 5.17f), with the highest response barely distinguishable at the sixth harmonic (S t = 0.099) and another significant value of the PSD located at the dominant frequency of the flow. Furthermore, the phase portrait of the pair (Fig. 5.18b) suggests a low level of correlation. This aspect can be observed in spatial eigenfunctions as well. In fact, it is hard to define a clear scale and shape of the structures, which depict a mixed character of single and double nature. This is clearer in ψ 3 (Fig. 5.17d), where the distribution of the blue valleys is more typical of single structures while the red regions depict a behavior more representative of double structures. A similar consideration can be applied to the imaginary component of the pair (Fig. 5.17e), in particular if one considers that the blue region does not define clearly the shape of the structure and it is not evident whether it represents a single structure or not. More, an evident wavelength cannot be recovered since the spatial periodicity is corrupted by the number and distribution of the correlated zones, most of all because the structures in the suction side of ψ 4 are hardly referable to a clear counterpart in ψ 3 . For this reason it is difficult to draw a conclusion on the connection between the pair content and a physical aspect of the flow. As regards the third pair, the spectral response is not, again, centered at a specific frequency, even though the highest value of the PSD (Fig. 5.17i) is located at S t = 0.165, i.e. the tenth harmonic of the shedding frequency. However, many other peaks can be identified: the two strongest are respectively located at the eleventh and fifth harmonic; as a result, the response spreads over a fairly broad bandwidth. The correlation of the time coefficients (Fig. 5.18c) depicts, as well, the influence of several frequencies in the time series and the prevalence of the values related to the first mode which show a stronger response. Although the correlation of the time domain may suggest an as much uncorrelated spatial eigenfunctions, the contour plots reveal an almost clear distribution of the structures, as well as a fairly defined wavelength beyond the trailing edge. Both the singles modes are evidently affected by the presence of background noise which, however, does not prevent the recognition of the shapes of the structures. The configuration present in the real component (Fig. 5.17g) can be clearly observed in the imaginary component (Fig. 5.17h), shifted in the wake direction and following a curved path connected with the pitching motion. The similar shape of the leading edge structures supports this consideration; however, the trailing edge structure in ψ 6 has clearly reduced sizes if compared with its counterpart in ψ 5 and, in addition, it introduces a fragmented region on the pressure side that is hardly referred to a specific structure in the real component. Although this mode represents for sure an alteration of the pressure field, however, the observation of the spectral response suggests that no conclusions can be drawn on the connection with physical phenomena, even though the higher wavelength of the structures on the suction side may be representative of the influence of the first mode, that is the one connected with the shedding vortex. The SPOD is performed by applying a filter with size equal to the time that the low pressure core of the main vortex needs to form at the leading edge and shed out of the domain of interest, i.e. N f = 35. As already mentioned, this action evidently reduces the energy featured by the first single mode, while the second one results to be enhanced (Fig. 5.19). This causes an increased distance between the second eigenvalue and the third, which depicts a slightly reduced energy content that makes it comparable with the fourth. In addition, the fifth and the sixth modes shows a stronger similarity at a barely augmented energetic level, while the following two eigenvalues experience a more distinct introduction of modes energy. As a results, the retained modes in this case, according to the chosen percentage, are 24. The 12 couples are plotted in the SPOD spectrum in Figure 5.20. As can be seen, the dots still follow a linear trend and no clear dominant pairs can be identified. The three chosen modes represent even in this analysis, the most energetic and most correlated ones, even though the energy contents of each are changed according to the previous considerations: in particular, they feature, respectively, 63%, 15% and 7%. The marked effect on the average frequencies can be seen in the general translation of the pairs to the left, which suggests a reduced influence on the higher harmonics; this is much more evident on the third pair that experiences a reduction of nearly half of the frequency it shows in the POD. A similar effect of the filter was observed on the same pair in the velocity field decomposition where, instead, the adoption of SPOD led to an increase of the average frequency. The first pair shows a clear narrower response at the dominant peak of the flow, while the influences of the other frequencies are evidently suppressed (Fig. 5.21c). In terms of time coefficients, both the modes show much similar values throughout the series; however, the influence of a higher frequency in the second mode causes the irregular shape of the phase portrait which, although smoother, still cannot be considered as a circle (Fig. 5.22a). As regards the spatial modes, the characters shown in the previous decomposition are still recognizable, although the fragmentation spreading into the wake is now stronger. The red structure in the first single mode (Fig. 5.21a) appears with reduced sizes, which lead the strongest region to cover half the chord length; the trailing edge structure is, instead, closer to the airfoil surface. The introduction of a higher correlated part directed towards the wake, as well as the tendency of the structure to increase as they move downstream, may be well representative of the behavior of the LEV that forms at the leading edge and thickens as it sheds rearward. More evident effects of the filter can be observed in the second pair. The spectral response is clearly smoothed, with a reduced bandwidth and a clearer peak centered at the fifth harmonic, which determines an important change if compared with the POD results that estimated the highest value occurring at the sixth harmonic ( Fig. 5.21f). The influence of the response at the first harmonic is still viewable, though with a lower impact on the average frequency that is now closer to the peak one (S t = 0.08). The phase portrait ( may be taken as a further evidence of the presence of this flow features oscillating at the given frequency, while a similar conclusion, based on the POD results, would not have been possible because of either the different modes coupled and the already discussed difficulty in stating the spatial correlation. The dynamics of the third pair depicts a drastic reduction of components with frequencies greater than S t = 0.20 and the peak response is now located at the sixth harmonic, even though significant contributions can be identified, most of all at lower frequencies ( Fig. 5.21i). The phase portrait in Figure 5.22c depicts an enhanced regularity, but still the limited circular portion suggests a good correlation only for some components of the main frequency. In general, the coherent structures of this pair define a clear repeating scheme that allows to infer the character of double structures of the spatial correlation. The real component (Fig. 5.21g), although characterized by a reduced background noise, still shows an incomplete separation of the red structures in the rear part of the domain, while the structures in the first part of the suction side are neither well defined nor referable to the other correlated regions. The imaginary component (Fig. 5.21h), instead, shows clearer definition and separation of the coherent structures, while the background noise is more present in the remaining part of the domain. As in ψ 6 , the correlated regions move towards the upper part of the wake, which helps to identify the spatial periodicity, while the defined structures over the suction side represent a clear element of difference with the other component. In particular, the directionality of the structures suggests that this pair contains the information of flow features typical of the fully separated conditions, occurring in the first part of the downstroke. The flow visualization supports this conclusion and, in particular, reveals that the second and the third pairs may represent, respectively, the upstream and downstream behavior of the phenomena evolving during this phase. In fact, all the low pressure regions originate on the suction side as single structures and travel downward until they shed into the wake where they induce higher pressure limited zones that form a double structure traveling away at the same velocity. A second SPOD analysis is performed even for this field and the adopted filter width is the same as in the previous section, N f = 63, i.e. the time spent by the flow under conditions of full separation. The extraction of energy from the first eigenvalue and the concomitant increased content of the second one, provide in this case a clearer similarity of the two modes ( Fig. 5.23). The modes from 3 to 6 show, in pairs, comparable energetic levels and a clearly enhanced prevalence on the subsequent eigenvalues which are still recognizable as coupled features. The redistribution of the total energy leads to the necessity of retaining the first 34 modes in order to reconstruct the 99% of the initial data. The spectrum of the coupled modes ( Fig. 5.24) shows distinct consequences of the action of the filter. As regards the pairs so far discussed, especially the second and the third, it can be observed that the average frequency is further, markedly, reduced which suggests a high impact on the dynamics of the structures; in addition, the energetic content of the two mentioned couples is more similar and the linear trend is more comparable with the behavior of the spatially-averaged PSD curve. However, a significantly correlated pair, featuring a distinct energy content, can be detected at S t = 0.164 (4), i.e. the tenth harmonic. A similar effect of the filter was not observed in the discussion developed so far, and then the mentioned pair will be discussed with the same attention as the other three. The four identified pairs are now representative, respectively, of the 59%, 12%, 8% and 4% of the original data set. The dynamics of the first pair is further centered at the first harmonic of the flow, with a response featuring a almost mono-frequency character and, then, negligible influences of the responses at higher frequencies (Fig. 5.25c). As a result, the phase portrait (Fig. 5.26a) approximates fairly well a circular shape expectable from a purely spectral analysis. In addition, no clear differences with the previous analysis can be detected in the spatial distribution of the coherent structures. The effects on the second pair are more marked and such that the dynamics shows a behavior closer to the DFT. In fact, when a purely spectral analysis is performed, this pair shows a single frequency response at the second harmonic; in this analysis, the peak is further reduced to coincide with the third harmonic, which denotes the gradual suppression of the influence of different frequencies responses made by the action of the filer (Fig. 5.25f). The time evolution of the coefficients shows an enhanced regularity, but still the influence of higher frequencies is recognizable, although with a lower impact (Fig. 5.26b). Even for this pair the effects on the coherent structures are not significant. In fact, a general reduction of the background noise and broadness of the single structures can be observed in both the real and imaginary components (Figs. 5.25d and 5.25e). However, in the former, increased regularity and similarity of the shapes allow for a better identification of the characteristic length scales of the pair. The spectral response of the third pair is still maximum at the sixth harmonic and even though some different peaks are suppressed, however two of them show a significant contribution (Fig. 5.25i). This multi-frequency character is also evidenced in the phase portrait of the real and imaginary components of the time coefficients (Fig. 5.26c). This pair shows a greater sensitivity towards the filter width, in terms of spatial eigenfunctions. In fact, the real component depicts a better distinction of the blue valleys, and the structures near the leading edge are gathered into a single structure (Fig. 5.25g). This aspect emphasizes the similarity with the imaginary component where, in fact, only one structure is present at the leading edge ( Fig. 5.25h). However, both the single modes show an increased impact of the background noise. As anticipated, a fourth pair is chosen for this case, i.e. the one referring to the fourth distinct dot in the SPOD spectrum. In particular, this couple shows a strong correlation even for lower and higher filter widths than that one used here, even though the modes constituting it change from time to time. This fact suggests the presence of fairly correlated flow structures occurring with this frequency, which has an average value of S t = 0.164. The peak frequency at S t = 0.165 ( Fig. 5.27c) supports the consideration that the impact of different frequencies is low, but still not negligible as reported in the phase diagram (Fig. 5.28). The band of the PSD also depicts a partition into two main ranges, where the second one is centered at the peak frequency and shows a higher response. In terms of spatial correlation, the eigenfunctions (Figs. 5.27a and 5.27b) report clearly (

Data Reconstruction
The reconstruction of the original data is an important analysis that allows for a better understanding of how new models can be built, based on the basis computed, with a reduced number of modes. This aspect is fundamental in the reduced order modeling, where predictions on given phenomena can be drawn by solving problems that use the POD space as input. As a result, there is no need to perform complete simulations for every condition of interest, since when one is performed directly the other can be modeled by taking advantage of its decomposition. A further benefit of this approach consists on the optimality of the decomposition, which allows for a high level of confidence even with a limited number of modes that results in a reduced computational cost.
Instead of comparing the reconstruction of the flow field of some snapshots, which may lead to cumbersome conclusions, it has been chosen to compare the reconstruction of the loads acting on the airfoil during a pitching cycle. In particular, by decomposing and reconstructing the pressure field, it is possible to integrate the values over the airfoil surface in order to obtain the loads coefficients and compare their values, quantitatively, with the corresponding original data from the TAU solution. The values reported and discussed below are obtained from the reconstruction of the pressure field with a number of modes such that the relative energy content is over the 95%; two cases are compared: the POD and the SPOD with N f = 35, which account for, respectively, 7 and 12 modes. It should be stressed that, since the contribution of the skin friction coefficient is not added to the pressure loads, the plots are only meant to represent an estimate of the real loads and part of the discrepancies can be attributed to this approximation. However, since the pressure loads represent the dominant effect, sufficient insights can be extracted from their behavior. Another aspect that should be clarified regards the sampling. In fact, the original solution stems from the monitored data obtained from the solver after every time step (3600 per cycle), while the reconstructed values arise from the integration of the reconstructed field, based on the snapshots that are taken every 36 time steps (100 per cycle). When the lift coefficient is considered, one can observe that the chosen reconstruction percentage provides, in general, a complete and precise estimate. The POD (Fig. 5.29a) shows a reduced influence of the oscillations during the upstroke phase, if compared with the SPOD (Fig. 5.29b) which, in turn, accentuates the effect of the highest harmonics. In fact, it can be noticed that the linear trend is reconstructed better by the POD, even though the wavy outline shown by the SPOD presents values that are closer to the original data. The change of slope occurring at 15°in the red curve appears anticipated in the POD reconstruction, which causes a further separation between the two lines. However, the SPOD depicts a better agreement during a major part of the extra lift phase. As regards the peak value, it is better recovered by the SPOD, even though this method anticipates its occurrence, while the POD evaluates its estimate at the same angle as the original plot. The downstroke phase until the second vortex shedding shows a similar behavior among all the curves, while the remaining part is better recovered by the SPOD. In fact, although showing values more smoothed than the original solution, it predicts fairly satisfactorily the occurrence of the several oscillations, while in the POD they appear closer in magnitude, but clearly delayed. The plateau preceding the reattachment of the flow shows a similar behavior between the two methods.
The drag coefficient shows enhanced properties of the SPOD (Fig. 5.30b) in the reconstruction of the upstroke phase, since it better recovers the curvature of the curve, which denotes a better

Further Considerations
description of the fluctuations occurring at lower frequencies. The POD (Fig. 5.30a), instead, depicts a closer result in the last part of the upstroke, and predicts correctly the position of the peak. However, it tends to overestimate the peak value, as well as the SPOD which, in addition, shows a delayed prediction of the angle. In both the reconstructions, the bifurcation of the peak, which has its counterpart in the moment curve, is completely missed. An increased number of the retained modes has shown that this behavior is connected with much weaker fluctuations that are recovered only with an almost full reconstruction, rather than with the filter adopted. This plot shows again very similar results for the first part of the downstroke. Clearer differences appear around α = 17.5°, starting from which the SPOD depicts a smoother trend, unable to recover the several fluctuations and generally overestimating the reference values. The POD, instead, reports a more wavy pattern but, still, with higher values. The end of the downstroke and beginning of the upstroke is underestimated by both the methods, which show similar behaviors. The highest discrepancies can be observed in the moment coefficient. The beginning part of the upstroke shows clearly higher values for both the methods but, while the POD (Fig. 5.31a) recovers a good agreement with the reference curve from about 12.5°, the SPOD (Fig. 5.31b) presents a wavy trend, connected with that shown in the lift curve, which produces clearly amplified values. The behavior is then similar from 17°. In fact, they both predict the minimum at the same angle (around 20.9°) that coincides, in the reference data, with a small peak between the two valleys.
In this regard, they show a similar behavior to that described in the drag curve, i.e. the smaller oscillations require more modes to be recovered. The downstroke phase depicts, as in the previous cases, a good proximity of the data, until the occurrence of the second vortex shedding for which, however, the POD slightly overestimates the moment value. The remaining part of the downstroke shows great discrepancies, especially as regards the ending phase. The SPOD tends to overestimate the values until 12.5°, though with a good prediction of the angles related to the fluctuations; from that angle to the end the values are clearly underestimated. Similar considerations can be drawn for the POD which, instead, shows generally undestimated values even for the first phase after the shedding of the secondary vortex and a overall tendency to delay the estimate of the angles.

Further Considerations
This section is intended to provide additional details for the analysis reported in the previous sections and a comparison between the two different field decompositions. As already reported, the first pair showed, in both the cases, clear uncorrelated spatial structures in the face of a fairly strong time correlation. The reason of this behavior was found in the content of the single modes constituting the couples, which describe both the vortexes forming throughout the pitching cycle and both the fields seemed to converge on this consideration. Furthermore, the results showed that the velocity field decomposition provides a better extraction of the structures connected with the second vortex, if compared with the other analysis. In fact, it can be noticed that it is able to describe also the small counterclockwise region arising at the quarter chord, as consequence induced by the low pressure core of the secondary vortex. The first mode shows, instead, a fairly confusing distribution of the structures connected with the primary vortex since they are depicted by the first single mode as strong fluctuations of the time averaged flow field of the w component (not reported). In contrast, the second single mode of the pressure decomposition provides a distinct description of the structures related to the main vortex, while the secondary vortex is not recognized as a different structure from the primary and are both represented as a single structure in the first single mode. The results so far discussed showed a greater sensitivity of the velocity field towards the fluctuations at the higher harmonics, while the pressure field depicted a clear response only at the dominant frequency of the flow. This aspect, originally suggested by the spatially averaged PSDs of the fields, was evidenced even by the analysis of the modes and related coefficients. In fact, the harmonics that depicted stronger responses showed even good time and space correlation, made further enhanced by the action of the filter. This was clear in the velocity field analysis, where different filtering of the correlation matrix did not cause drastic changes in those pairs that showed already high levels of correlation, but, rather, performed augmented definition of the structures and reduced presence of the background noise. In contrast, the pressure field analysis showed that the action of the filter is harbinger of strong changes in the dynamics of the identified modes, which depict noticeable variations of either the peak and average frequencies. In particular it was observed that the fourth pair analyzed at the end of the third section is already detected by filters with a width in the range between 35 and 63. However, in those cases almost the same features can be seen as referred to a couple constituted by two completely different modes (7 and 8, in the investigated situations). This interesting finding may be interpreted in two different ways: as a further evidence of the existence and importance of those flow structures and as additional prove of the strong impact of the filter on the dynamics of the pressure field. In particular, this latter suggests that a greater attention should be paid when the SPOD is applied to such fields, as the pressure field, in the context of dynamic stall. In the cases presented in this work, the filter sizes are chosen such that they correspond to the time needed by the main vortex structure to form and leave the measurement domain, interpreted according to the flow features that describe it in each field. A second filter is chosen, independently of the field considered, as the time spent by the flow under conditions of full separation. In particular, the first class of filters provided, in general, better descriptions of the coherent structures connected with the second pair, while the second class, in this sense, produced better results for the third couple. However, this consideration is more suitable if the velocity field is considered, since the dynamics of a chosen couple is preserved between different applications, and the differences can be properly discussed. In case of the pressure field, instead, the strong variations experienced by the the dynamics of the pairs make it difficult to judge whether the filter is improving or not the extraction of the new structures, since there is no way one could be aware whether the differences relate to the change of the spectral content or, rather, to an enhanced correlation. As a concluding remark, the data reconstruction showed that, in general, the SPOD tends to produce results that better predict the occurrence of the flow fluctuations. However, for the same reconstruction percentage, the POD evidenced a greater ability in recovering proper slopes and magnitudes of the loads. In view of building an efficient reduced order model, this consideration gives a way to decide how to approach the choice of the basis. Clearly, a good trade-off allows for the identification of the most suitable function space. In fact, reduced dimensions allow for a saving of computational resources but do not take into account the weaker variations, while more sensitive models may be too onerous in terms of computational costs. For this reason, a proper analysis of the filter width is required. In fact, one should select the filter in such a way that only the structures of interest are amplified so that the extraction of energy from the dominant modes does not reach a condition where more modes are needed to compensate for this deficit.

Conclusions and Outlooks
Spectral Proper Orthogonal Decomposition (SPOD) is applied to investigate the coherent structures connected with dynamic stall. The data for the analysis were obtained from the CFD solver DLR TAU, which solved the RANS equations with the SST turbulence model for a pitching NACA 0012 airfoil at Re= 1.569 · 10 6 and M = 0.445, oscillating about the quarter chord with k= 0.135, α 0 = 13.18°and α 1 = 9.548°. For each cycle 100 outputs (or snapshots) are evenly taken. The last two cycles, i.e. 200 snapshots, are decomposed and analyzed with the stochastic tool. The decomposition is performed on two different field types: the vector-valued velocity field and the scalar-valued pressure field. Pure POD is initially adopted in order to consider the energetically optimal basis and the corresponding spectral response. Then its results are compared with the ones from SPOD that is applied by filtering the POD correlation matrix along the diagonals with a Gaussian distribution: two different filter widths are chosen for each field. A reconstruction of the field is performed in order to consider strengths and weaknesses of the possible application of SPOD in the context of reduced order modeling. This analysis is done by comparison of the airloads originally obtained by TAU with the ones computed from a partial reconstruction of the pressure field, for POD and a case of SPOD.

Spectral POD and Dynamic Stall
The results showed that SPOD is a suitable tool for the analysis of coherent structures underlying flows under conditions of deep dynamic stall. The analysis proved that a good choice for the filter width is of a time scale of characteristic flow events. In particular, in this work the first filter width was chosen equal to 30 and 35, respectively for the velocity and pressure field, which correspond to the time needed by the main vortex structure to form and shed out of the measurement domain, basing on the flow visualization of the corresponding field. The second filter size was chosen, in both the cases, equal to 63, which represents the time spent by the flow under full separation conditions. The scalar-valued field showed a nearly linear spectral response, with a single clear peak at the natural frequency of the flow, S t = 0.017, corresponding to the pitching frequency. This aspect seemed to be connected with the high sensitivity of the pressure field towards the filter adopted, since significant changes in the dynamics of the coherent structures were depicted by subsequent changes of the filter size. In particular, it was noticed that the same structures were observable in different pairs of modes, as the filter width was increased: this evidence raised the question whether the variations were to be attributed to an enhanced correlation rather than to a progressive "migration" of the dynamic content towards other modes. The vector-valued field, instead, showed a higher stability towards the filter size, which in fact allowed for a strengthened correlation of the pairs around their dominant frequency, without clear shifting of the dynamic content. This resulted in a general reduction of the background noise and an associated increased definition of the shapes and periodicity of the coherent structures. Two modes were identified at the fifth and eleventh harmonic through the action of the filter that increased the phase correlation of the corresponding pairs. They were interpreted, respectively, as the upstream and downstream evolution of the several eddies developing during the phase of the flow under fully separated conditions. For both the analysis the first mode was found to be connected with the primary and secondary vortexes, where the velocity decomposition provided a clearer description of the secondary vortex while the evolution of the primary was better represented by the pressure decomposition. The analysis of the reconstructed field showed that, for the same fraction of relative energy content, the POD requires a fewer number of modes than the SPOD, due to the redistribution of the energy among the less energetic modes performed by the latter. Furthermore, it was generally observed that the SPOD better predicted the angles corresponding to the several fluctuations of the coefficients curves, while closer values to the true solution were recovered by the POD.

Future Works
The adoption of the RANS technique employed in the present work clearly acted as an additional filter since it restricted the modal analysis to the biggest scales of the considered flow. The inability of this method to recover the last detail of the phenomenon is also clearly evidenced by the aforementioned literature. In order to account for the influence of the smallest scales and the importance they may have on the biggest, a deeper CFD analysis should be considered, e.g. by means of finer techniques such as LES. Such an analysis would also provide further proves on the dynamics described here for the main structures or, in contrast, evidence the decisive impact of the dissipation occurring at the smaller scales. In addition, the complete 3D character of the dynamic stall phenomenon over the blade may be investigated through the decomposition of full rotor simulation data such as to obtain additional insights on the structures evolving in the physical domain.
In conclusion, the building of a reduced order model based on the SPOD basis would represent an interesting benchmark for the determination of the method suitability to such applications. The subsequent validation of the results with the data obtained, for the same case, by numerical simulations or even experimental measurements would give a proper yardstick for the strengths and weaknesses of such an approach an the benefits it may produce in terms of computational resources.