1. Introduction
Increasingly, research efforts in the aviation community are focusing on reducing emissions and improving the efficiency of aircraft, to meet more stringent regulations associated with current noise and climate targets (European Commission 2011). Despite significant progress in addressing aircraft noise over the last half-century, turbomachinery remains one of the leading contributors to the acoustic signature of aviation to this day. While some of the acoustic sources, such as jet noise, have been reduced considerably, other sources, such as fan and combustion noise, are starting to dominate in terms of their relative importance. Addressing these is more challenging, however, because in some cases it requires a detailed understanding of flow instabilities and how they interact with the surface geometries and the flow field.
More powerful supercomputer resources are now allowing some of these issues to be addressed in a systematic way by facilitating direct calculations of flows and their acoustic footprint within the framework of one computational code. Computational aeroacoustics has thus re-entered a ‘second golden age’ (Lele & Nichols Reference Lele and Nichols2014) where the interaction of hydrodynamic structures and sound fields can be simulated with a high degree of precision and fidelity.
Despite remarkable progress, however, a significant portion of such investigations has relied on performing direct numerical simulations (DNS) and processing of the results using standard data analysis techniques. Even though it is acknowledged that many of the known compressible instabilities observed in these set-ups are due to an interplay between hydrodynamic instabilities and acoustic feedback, the majority of analyses have concentrated solely on the solution of the governing equations (see, e.g. Colonius & Lele Reference Colonius and Lele2004; Sandberg Reference Sandberg2013, Reference Sandberg2015; Sandberg et al. Reference Sandberg, Michelassi, Pichler, Chen and Johnstone2015) or the formulation of a global stability problem (see, e.g. Fosas de Pando et al. Reference Fosas de Pando, Schmid and Sipp2014, Reference Fosas de Pando, Schmid and Sipp2017; Sidharth et al. Reference Sidharth, Dwivedi, Nichols, Jovanović and Candler2021; Glazkov et al. Reference Glazkov, Fosas De Pando, Schmid and He2023a , Reference Glazkov, Fosas De Pando, Schmid and Heb ), as local stability efforts, that focus on the hydrodynamic and acoustic parts separately, often fail to describe complex behaviour or identify frequential selection principles. Only a global point of view, where instabilities, feedback and receptivities work together as a single, closed-loop unit, can provide insight into the aeroacoustic flow generated by a rotor–stator configuration. Within this viewpoint, even a weak feedback by acoustic waves can constitute a substantial component in the overall instability, when coupled to a strong receptivity at an upstream location.
While a global analysis of closed-loop behaviour is paramount for any proper analysis, an even more detailed picture arises from a dual approach. It effectively assesses the sensitivity to parameter, flow field or geometric modifications of each constitutive part of the global instability. As such, it presents a quantification of robustness and resilience of the underlying global instability and – for unsteady flows – allows for the dissection of observed flow behaviour into a chain of individual subprocesses.
This dual (or adjoint) approach to describing instabilities of unsteady flows will be taken and advocated for in this article and demonstrated on the aeroacoustic flow through a rotor–stator model. To this end, the forward DNS code will be coupled with an adjoint component, which arises when formulating the flow model through an optimisation problem (see e.g. Gunzburger Reference Gunzburger2002). In so doing, solving the dual problem for the adjoint variables furnishes gradient information for the flow with respect to the chosen control parameters. This information is critical for systematising the identification and quantification of sensitive structures, and the decomposition of an overall intricate instability mechanism into its simpler, but globally linked, components.
The traditional approach to this kind of analysis is typically to select a linearisation point around which the dynamics is studied. For unsteady, high-Reynolds-number flows, past studies (see, e.g. Sipp & Lebedev Reference Sipp and Lebedev2007; Oberleithner, Rukes & Soria Reference Oberleithner, Rukes and Soria2014; Beneddine et al. Reference Beneddine, Sipp, Arnault, Dandois and Lesshafft2016) have often chosen a time-averaged solution (i.e. a mean flow) as a linearisation point. In addition, in the context of rotor–stator interactions, they mostly concentrated on purely aerodynamic optimisations (see, e.g. Thomas, Kenneth & Earl Reference Thomas, Kenneth and Earl2005; Papadimitriou & Giannakoglou Reference Papadimitriou and Giannakoglou2006; Wang & He Reference Wang and He2010a , Reference Wang and Heb ) which are less afflicted by transient fluid motion and feedback effects. In contrast, in our case, we consider the unsteady flow through a moving geometry and hence must abandon the concept of a steady linearisation point. Instead, a linearisation about nonlinear flow snapshots is performed, which effectively focuses on the sensitivities of a linear-tangent model, i.e. a linear approximation of the flow in tangent paths along a nonlinear phase-space trajectory. The linkage of the time-varying (primal, direct) linearisation point to the (dual, adjoint) sensitivities requires a special computational method referred to as checkpointing (Wang, Moin & Iaccarino Reference Wang, Moin and Iaccarino2009), which addresses the issue of storing the full nonlinear primal solution by instead storing an optimal subset of snapshots from which missing portions of the trajectory can be reconstructed as needed.
Besides the optimisation formalism, the realisation of this dual analysis puts overly stringent requirements on the numerical schemes. The stator and rotor sections of the compressor stage model introduce a sliding mesh across which the flow and acoustic waves must propagate (see figure 1 for a sketch). This propagation has to be highly accurate, so that the dual flow fields (sensitivities) are not contaminated by numerical artefacts, which may degrade the physical sensitivity signal. A special sliding-interface scheme has thus been developed that enables our dual analysis across a moving mesh and avoids unphysical distortion and dissipation effects in the dual fields (see Glazkov et al. Reference Glazkov, Fosas de Pando, Schmid and He2025 for details).
Our focus is on the type of analysis for complex unsteady, moving-geometry flows with inherent acoustic feedback loops. For this reason, the stage geometry is simplified by considering linear cascades consisting of slender ellipsoids for the stator and controlled-diffusion airfoils for the rotor. Despite this geometric curtailment, this configuration exhibits complex flow interactions that are present in more complex configurations. This configuration has previously been investigated in low-Mach-number aeroelasticity experiments and in numerical studies at low Reynolds numbers (see Glazkov et al. Reference Glazkov, Fosas De Pando, Schmid and He2023a , Reference Glazkov, Fosas De Pando, Schmid and Heb ). In these investigations, it has been shown that each single passage is characterised by an oscillating laminar separation bubble shedding eddies downstream, which later interact with vortices generated at the trailing edge and create upstream-propagating acoustic waves. This latter component interacts with the region of maximum receptivity, upstream of the separation bubble and near the leading edge on the pressure surface, thus closing the hydrodynamicacoustic feedback loop. In a multi-passage layout these feedback loops manifest themselves as resonances and large-wavelength disturbances throughout the domain.
The dual analysis that follows is devoted to (i) identifying the flow and instability features arising as a result of the hydrodynamic–acoustic interaction processes and (ii) correlating these with adjoint sensitivities in a qualitative way.
2. Numerical method
2.1. Aeroacoustics flow solver
The computational aeroacoustics solver used in this work is based on the pseudo-wave method of Sesterhenn (Reference Sesterhenn2000) and implements the non-dimensional compressible gas flow equations for a calorically perfect gas. These equations are then discretised on a smooth structured grid, with a differentiable grid metric. A fifth-order compact, low-dispersion upwind scheme (Adams & Shariff Reference Adams and Shariff1996) is used for the differentiation of the inviscid terms. The dissipation terms are differentiated using a third-order compact scheme (Lele Reference Lele1992). These schemes are implemented with their locally compact variant, with sixteen ghost cells around each domain region, to facilitate the required parallel scaling performance on parallel clusters. The nonlinear solver and the checkpointing algorithm (Griewank & Walther Reference Griewank and Walther2000) both use a low-storage variant of the Runge–Kutta 4 algorithm RK4(3)5[2R+]C (Kennedy, Carpenter & Lewis Reference Kennedy, Carpenter and Lewis2000) to integrate the system of equations in time. The resulting system of equations can be represented in concise form as the initial value problem
\begin{gather} \frac {{\rm d}\boldsymbol{v}}{{\rm d}t} = \boldsymbol{F}(\boldsymbol{v},t), \qquad \text{with initial condition} \qquad \boldsymbol{v}(t=0) = \boldsymbol{\eta }. \end{gather}
In the above expression, the state-vector components
$p, \boldsymbol{u}, s$
have been cast into a composite vector
$\boldsymbol{v}.$
Further details of the underlying code and implementation details may be found in the published literature (Fosas de Pando Reference Fosas de Pando2012; Fosas de Pando et al. Reference Fosas de Pando, Schmid and Sipp2014, Reference Fosas de Pando, Schmid and Sipp2017; Glazkov et al. Reference Glazkov, Fosas De Pando, Schmid and He2023a
,
Reference Glazkov, Fosas De Pando, Schmid and Heb
, Reference Glazkov, Fosas de Pando, Schmid and He2025).
2.2. Linearisation
Before describing the details of the linearisation algorithm, it is important to stress that we employ a linear-tangent approximation about the current flow state
$\boldsymbol{v},$
rather than a steady or averaged base flow. The conjugate transform of the local (linear-tangent) Jacobian, associated with (2.1), is used for the adjoint analysis.
At the core of the method described in this paper is the extraction of the linearised dynamics from (2.1) through the system
\begin{gather} \frac {{\rm d}\tilde {\boldsymbol{v}}}{{\rm d}t} = \left . \frac {\partial \boldsymbol{F}}{\partial \boldsymbol{v}} \right \vert _{\boldsymbol{v}} \tilde {\boldsymbol{v}} \quad \equiv \quad \boldsymbol{A}(\boldsymbol{v},t) \ \tilde {\boldsymbol{v}}, \end{gather}
and its adjoint equivalent, obtained through the operator
$\boldsymbol{A}^H$
. While Jacobian-vector products can be obtained by simple quasi-linearisations of the form
\begin{gather} \boldsymbol{A}(\boldsymbol{v},t) \tilde {\boldsymbol{v}} \approx \frac {\boldsymbol{F}(\boldsymbol{v}+\epsilon \tilde {\boldsymbol{v}}, t) - \boldsymbol{F}(\boldsymbol{v}, t)}{\epsilon }, \end{gather}
direct access to the linear operator
$\boldsymbol{A}$
and, crucially,
$\boldsymbol{A}^H$
is not possible with this approach. Furthermore, the numerical scheme’s compact differentiation stencils introduce spatial coupling that requires sequential perturbations to every degree of freedom for element-wise constructions of
$\boldsymbol{A}$
.
The directed-graph-based linearisation method used in this work instead breaks down
$\boldsymbol{F}$
into a logical sequence of constituent nonlinear sub-blocks that depend only on a quantity vector
$\boldsymbol{q}$
, and its spatial derivatives
${\partial \boldsymbol{q}}/{\partial \boldsymbol{x}}$
, where the quantity vector is either the state-vector
$\boldsymbol{v}$
, or the first spatial derivative of the state-vector
${\partial \boldsymbol{v}}/{\partial \boldsymbol{x}}$
. With the differentiation operators moved outside of the sub-blocks, each component is now spatially decoupled and can be linearised at all points simultaneously, with the solution assembled by propagating
$\tilde {\boldsymbol{v}}$
through the sequence of linearised sub-blocks and differentiation operators. This method of extracting the linearised operator
$\boldsymbol{A}$
on-the-fly also allows for immediate access to the corresponding adjoint operator
$\boldsymbol{A}^H$
, by simply reversing the direction of the logical sequence and taking the transpose of all linearised sub-blocks. For efficiency reasons, only matrix-vector products are constructed to avoid dense linear operators. A detailed exposition of this methodology can be found in Fosas de Pando, Sipp & Schmid (Reference Fosas de Pando, Sipp and Schmid2012).
2.3. Sliding interface and unsteady adjoint
The critical enabling method used in this work is the treatment of the sliding interface between the rotor and stator mesh blocks. Through the application of the high-order interpolation method of Delfs (Reference Delfs2001), and exploiting the block structure of the underlying interior numerical scheme (Fosas de Pando et al. Reference Fosas de Pando, Sipp and Schmid2012), numerical dissipation is minimised in the forward problem. The interior and sliding-plane schemes are constructed so as to share a common numerical structure from which the adjoint implementation follows directly through the transposition of the underlying linear algebra and processor communications logic, thus ensuring that no additional sensitivities are introduced into the adjoint (dual) fields.
Functions responsible for the formation of the matrix-vector products
$\boldsymbol{Av}$
and
$\boldsymbol{A}^H\boldsymbol{v}$
, and the evaluation of the nonlinear right-hand side
$\boldsymbol{F}(\boldsymbol{v},t)$
, are integrated into wrappers in PETSc (Abhyankar et al. Reference Abhyankar, Brown, Constantinescu, Ghosh, Smith and Zhang2018; Balay et al. Reference Balay, Gropp, McInnes, Smith, Arge, Bruaset and Langtangen1997, Reference Balay2019, Reference Balay2020) to make use of the high-performance linear algebra methods and checkpointing schemes implemented therein (Zhang, Constantinescu & Smith Reference Zhang, Constantinescu and Smith2022). Nonlinear-adjoint looping with checkpointing is used to compute the adjoint sensitivity fields for the moving-interface problem. Complete implementation details and validations for the sliding plane and the nonlinear-adjoint looping are provided in Glazkov et al. (Reference Glazkov, Fosas de Pando, Schmid and He2025).
2.4. Flow domain and conditions
Figure 1. Geometry of the rotor–stator demonstration case. The stator stage consists of two slender elliptical blades, while the rotor stage is made up of three controlled-diffusion airfoils. This configuration requires a sliding interface due to the relative motion of the rotor grid with respect to the static stator mesh. Periodic boundary conditions are imposed in the cross-stream coordinate direction. Sponge layers are placed at the inlet and outlet to suppress reflections at the boundaries.
The numerical domain for our set-up consists of a periodic mesh with two stators placed upstream of three rotor blades (see figure 1 for a sketch). The stator blades are
$1$
-to-
$10$
aspect-ratio ellipses placed in a uniform axial flow and staggered so that they are at a
$3^\circ$
angle of attack relative to the incoming flow, resulting in unsteady von Kármán wakes shedding from the stator trailing edges. We assume that the flow turning introduced by these stator blades is insignificant, and that the resulting flow remains predominantly axial before interacting with the rotor blades. The rotor blades are controlled-diffusion airfoils developed by Sanger & Shreeve (Reference Sanger and Shreeve1986) and Sasaki & Breugelmans (Reference Sasaki and Breugelmans1998). All simulations presented here are performed at a chord-based Reynolds number of
${Re} = 100\, 000$
, an exit Mach number at the rotor blade row of
$M = 0.3$
and the quantities are non-dimensionalised by the outlet flow conditions in the rotor domain. The inlet flow angle in the stationary reference frame of the blades is set to
$\theta _{in} = 37.5^\circ$
to match the inlet angle of the rotor blades, and the non-dimensional pitch length for the periodic boundary of the linear cascade is set to
$h_p = 0.6$
. These flow conditions are adapted from an experiment and previous numerical studies (He & Yi Reference He and Yi2017; He Reference He2021). By decreasing the Reynolds number from
${Re} = 195\,000$
to
${Re} = 100\,000$
, and increasing the Mach number from
$M = 0.05$
to
$M = 0.3$
, a coarser mesh, and consequently faster time stepping, can be used for our demonstration case – a choice that reduces the computational load without compromising the macroscopic flow behaviour. At these conditions, the flow exhibits a laminar separation bubble, located at approximately
${C_{{ax}}}/{2}$
on the suction side, where
$C_{{ax}}$
is the axial chord of the rotor blade, that acts as an amplifier for boundary-layer disturbances. The gap between the rotor and stator blade rows is set to
$C_{{ax}}/2$
.
The flow in the upstream region of the rotor blade row is used as the basis for setting the initial conditions for the stators. The axial projection of the velocity in the rotor reference frame is used as the axial velocity in the stator mesh, while the vertical component is extracted for use as the non-dimensional mesh velocity,
$v_{{sp}}=-0.791$
. The case is then run for
$20$
time units so that the initial transients are removed from the flow field, to produce the undisturbed initial condition
$\boldsymbol{v}_0$
, before initialising the dual (adjoint) sensitivity analysis. The grid parameters for the rotor and stator meshes are listed in table 1.
Table 1. Main grid parameters used for the rotor–stator demonstration case. The parameters
$n_\star$
denote the number of mesh points in the labelled direction
$\star$
. Here,
$L_\star$
denotes axial length of region
$\star$
,
$h_{{passage}}$
stands for the height of the passage between two rotor blades in the cascade and
$\min \backslash \max {\Delta s_i}$
define the minimum and maximum arclength between grid points along the curvilinear grid in the two directions
$i=0,1$
.
3. Results
3.1. Nonlinear forward problem
The cost functional considered in the context of the Lagrange optimisation problem, by means of which this adjoint sensitivity analysis is formulated, is taken to be
\begin{gather} \mathcal{J} = \frac {\langle \boldsymbol{v}(\tau ) - \boldsymbol{v}_{\textit{targ}},\boldsymbol{v}(\tau ) - \boldsymbol{v}_{\textit{targ}} \rangle _{\boldsymbol{M}}}{2}, \end{gather}
where
$\boldsymbol{M}$
is a matrix representing the discretised Chu-norm inner product (see Chu Reference Chu1965, Hanifi, Schmid & Henningson Reference Hanifi, Schmid and Henningson1996 and Fosas de Pando Reference Fosas de Pando2012 for further details). In this expression
$\boldsymbol{v}_{{targ}}$
and
$\boldsymbol{v}$
represent the end states, at
$\tau =2$
, of two nonlinear simulations run within the time interval
$t\in [0, \tau ]$
. The interval length of
$\tau$
has been chosen to ensure a meaningful extraction of sensitivity information from the back-propagated adjoint fields and to avoid the necessity for corrective measures, such as Least-Squares Shadowing (see also § 3.2). At
$t=0,$
the initial conditions,
$\boldsymbol{v}(t=0) = \boldsymbol{\eta }$
, for these simulations are
\begin{align} \boldsymbol{\eta }_{{targ}} = \boldsymbol{v}_0, \\[-29pt] \nonumber \end{align}
\begin{align} \boldsymbol{\eta } = \boldsymbol{v}_0 + \mathcal{G}\left (\alpha \boldsymbol{\psi } \boldsymbol{\cdot }\boldsymbol{m}\right ), \end{align}
where
$\boldsymbol{\psi }\in [-1, 1]$
is a uniformly distributed random variable that is adjusted so that the pointwise sum vanishes, i.e.
$\sum _{i} \boldsymbol{\psi }_i=0$
, and is then multiplied by the amplitude
$\alpha =0.01$
. The mask function
$\boldsymbol{m}\in \{0, 1\}$
restricts the random perturbation to
$-1.5\leqslant x \leqslant -1.0$
, and sets to zero a halo of thirty grid points in the vicinity of the stator blades, to ensure that the perturbation does not affect any impermeability and no-slip conditions at the blade surfaces. A two-dimensional Gaussian filter
$\mathcal{G}$
is then applied to this field twice to smooth out the perturbation field, thus rendering all initial condition fields differentiable. The initial condition and the final state are illustrated for the entropy field in figure 2.
Figure 2. The initial and final states of the entropy component of the state vector showing the propagation of a random perturbation through the blade passage. (a) s, t = 0.00, and (b) s, t = 2.00.
In effect, we follow a block of random perturbations, superimposed on an unsteady but equilibrated flow field, through the stator and rotor section of the blade passages to probe and assess any hydrodynamic and acoustic flow features they may trigger over a full flow-through horizon. The dual analysis of this forward experiment will then give insight into essential and dominant components of the instability sequence.
3.2. Sensitivity analysis
We proceed by considering the sensitivity of the flow to initial conditions, as measured by the cost functional
$\mathcal{J}$
.
The adjoint equations can be derived by adding the governing equations (2.1) via Lagrange multipliers (or dual/adjoint variables
$\boldsymbol{w}$
) to the cost functional (3.1) to form an augmented Lagrangian
${\mathcal {L}}.$
This composite functional
$\mathcal {L}$
is then rendered stationary with respect to all its dependent variables, namely, the state-vector
$\boldsymbol{v},$
the dual state-vector
$\boldsymbol{w}$
and the control variable
$\boldsymbol{v}_{{targ}}.$
Setting the first variation to zero for all three dependencies yields three governing equations: one nonlinear evolution equation for the state-vector
$\boldsymbol{v},$
one linear but variable-coefficient evolution equation for the dual state-vector
$\boldsymbol{w}$
and one algebraic optimality condition from where we deduce the cost functional gradient. The dual state-vector
$\boldsymbol{w}$
consists of adjoint variables corresponding to primal state vector, i.e. the adjoint pressure
$p^\dagger ,$
the adjoint velocity field
$\boldsymbol{u}^\dagger $
and the adjoint entropy
$s^\dagger .$
In what follows, we commence with random-noise input to trigger instabilities in our rotor–stator configuration, which are subsequently analysed by an adjoint sensitivity investigation. For this reason, only a single primal-dual pair is required, and the final optimality condition can be safely discarded for our purpose. The coupled set of primal-dual evolution equations is then given as
\begin{align} \frac {{\rm d}\boldsymbol{v}}{{\rm d}t} = \boldsymbol{F}(\boldsymbol{v}, t), \qquad \text{with initial condition} \qquad \boldsymbol{v}(t=0) &= \boldsymbol{\eta }, \end{align}
\begin{align} -\frac {{\rm d}\boldsymbol{w}}{{\rm d}t} = \boldsymbol{A}^H(\boldsymbol{v}) \boldsymbol{w}, \qquad \text{with terminal condition} \qquad \boldsymbol{w}(t=\tau ) &= {{\rm d}\mathcal{J}}/{{\rm d}\boldsymbol{\eta }}. \\[10pt] \nonumber \end{align}
The primal equation (3.4) is solved forward in time over the time interval
$t \in [0, \ 2],$
starting from the initial condition
$\boldsymbol{\eta },$
while the dual equation (3.5) is evolved backward in time from
$t=\tau =2$
to
$t=0,$
starting with the terminal condition
$\boldsymbol{w}(\tau ) = {d\mathcal{J}}/{d\boldsymbol{\eta }}.$
The dual variables carry the sensitivities of all corresponding primal variables with respect to variations in
${\mathcal {J}}.$
We again note that the coefficients in the dual equation (3.5) depend on the solution
$\boldsymbol{v}$
of the primal equation (3.4). The implemented checkpointing algorithm (Griewank & Walther Reference Griewank and Walther2000) relies on an efficient flow-field management system that gathers solutions of the primal problem for injection into the coefficients of the dual problem.
Chaotic flow dynamics result in adjoint fields that grow exponentially in time, and therefore the range of the adjoint fields needs to be corrected to remove this exponential growth. The infinity norm of the adjoint pressure
$p^\dagger$
is shown in figure 3(a), effectively displaying the divergence due to a positive Lyapunov coefficient (Wang Reference Wang2013). The apparent growth of this norm can be approximated by the exponential
$10^{-5.208t}$
(note that time decreases from
$t=2$
to
$t=0$
in the adjoint). This does mean, however, that certain features are non-resolvable with this normalisation procedure. Where it is appropriate, this is adjusted and clearly stated so that observations may be made for all features of interest. A corrected or normalised adjoint pressure signal is shown in figure 3(b).
Figure 3. Evolution of the infinity norm of the adjoint pressure field
$p^\dagger$
. (a) Infinity norm (maximum value) of the raw adjoint pressure signal, showing a divergence due to a positive Lyapunov exponent stemming from the chaotic nature of the underlying flow field. (b) Corrected adjoint pressure signal, where the diverging part of the signal has been removed.
The peaks in figure 3 correlate with rotor–wake interaction events, when the stator wakes communicate with the leading edges of the rotor. These features are associated with short-time-scale pressure sensitivities that can only act within a compact time interval of an interaction event, due to faster propagation speed of the pressure component.
Trajectory divergence issues, and, consequently, positive Lyapunov exponents in the adjoint, are common features of dynamical systems exhibiting chaotic structures. Unbounded exponential growth resulting from this will degrade solution quality, eventually distorting the relevant physical structures and obscuring them below the resolving power of a finite precision machine. By considering laminar cases, or cases with low levels of turbulence, and integrating them over short time scales, this exponential growth is not given the time to degrade the solution and the growth can be factored out to give a bounded representation of the solution dynamics that remains quantitatively valid. As this analysis is pushed towards higher turbulence levels, or configurations that require longer integration times to capture larger feedback loops, techniques such as least-squares shadowing (Ni & Wang Reference Ni and Wang2017; Blonigan Reference Blonigan2017) are potentially viable alternatives to the approach presented here.
3.2.1. Rotor–wake interaction
In the forward problem, starting with random perturbations about a reference flow field, interaction events between the rotor and the stator wake are observed to be the dominant disturbance amplification mechanism. These events consist of the rotor’s leading edge slicing the vortices in the stator wake, after which a pressure pulse emanates from the leading edge. Instabilities are subsequently triggered within the boundary layers of the rotor blade as the sliced vortices lift up low-velocity fluid from the blade surface. This interaction process is shown in the panels of figure 4, visualised by the dilatation field
$\vert \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} \vert$
with
$\boldsymbol{u}$
as the velocity vector. It shows the collision of shed vortices with the leading edge of the rotor (panel (a)), followed by the generation of a pressure pulse (panel (b)) which, in turn, feeds back disturbances toward the stator (panel (c)). The impingement of wake vortices on the suction-side boundary layer is also clearly visible (panel (c)).
Figure 4. Impact of a vortex, shed from the stator blade, with the rotor’s leading edge (a), triggering a pressure wave that propagates omnidirectionally (b, c). The flow field is visualised by the dilatation field
$\vert \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} \vert$
with
$\boldsymbol{u}$
as the velocity field. (a)
$\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}$
, t = 1.70, (b)
$\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}$
, t = 1.75 and (c)
$\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}$
, t = 1.90.
Combined with the observed interactions in the forward vortex-impingement problem, the dual/adjoint analysis identifies structures of sensitivity that would, if applied, at later times maximise the result of the interaction in the final state
$\boldsymbol{v}(\tau )$
. It is worth stressing also that the adjoints propagate and grow backwards in time because smaller perturbations to the forward simulation at earlier times will result in larger changes in the final state, and the adjoint sensitivity decreases in time as the perturbation is less able to affect the flow at the end of the simulation
$(t=\tau ).$
Moreover, we choose to display the primal and dual solutions synchronised in time, despite the fact that the dual solutions propagate retrograde in time. As a result, the disappearance of sensitivities (adjoint fields) immediately pinpoints significant events in the primal fields.
To begin with, tilted Orr-like structures are first observed in the vicinity of the trailing-edge boundary layers on the stator blades (see figure 5
a). These structures represent optimal forcing modes that would maximise the disturbance to the stator wakes for the portion of the wake that later impacts the leading edge. As the structures tilt into the wake’s shear direction, they gain in kinetic energy, before impinging on the rotor. A sensitivity link to the rotor’s leading edge arises (see the thin adjoint filament in all panels of figure 5), anticipating the eventual collision of a particular vortex within the vortex street with the leading edge of the rotor at time
$t\approx 1.7.$
In other words, the adjoint flow field targets the vortical structure and concentrates on the shear layer of the specific vortex that later will impact the leading edge and will contribute maximally to the acoustic feedback via a radiating pressure wave.
Figure 5. Propagation of adjoint sensitivities (illustrated by the adjoint entropy field) from Orr-like structures at the stator trailing edge and to the resulting vortices in the free stream. Vortex cores have been highlighted with black contours. (a) s †, t = 1.00, (b) s †, t = 1.30, and (c) s †, t = 1.50.
This concentration in the adjoint is argued to have two reasons. First, optimising the vortex position and its vorticity will result in more energetic interaction events with the leading edges, the result of which is an increase of the cost functional
$\mathcal{J}$
. Secondly, it is known, from previous numerical investigations (Glazkov et al. Reference Glazkov, Fosas De Pando, Schmid and He2023a
), that at these flow conditions the sensitivities on the suction surface of the rotor blades outweigh the sensitivities on the pressure surface by orders of magnitude. This can be attributed to the amplifying effects of the laminar separation bubble located approximately at the axial mid-chord on the suction side. Furthermore, it was also established that this separation bubble, and a portion of the flow on the pressure surface near the leading edge, act as wavemaker regions that are pivotal to sustaining an internal feedback loop that links convective boundary-layer structures, their scattering into acoustic waves by the trailing edge and their retriggering of upstream boundary-layer perturbations. Modifying any component of this feedback loop would have significant global effects on the entire flow field.
Figure 6. Adjoint pressure pulse targeting the rotor’s leading edge at the moment of vortex impact. The pressure wave consists of a concentric and a non-concentric component (see (a)): while the former gradually disappears (b), the latter persists until impact with the leading edge (b, c). (a) p †, t = 1.65, (b) p †, t = 1.70, and (c) p †, t = 1.76.
Immediately prior to a collision, an adjoint pressure wave converges on the leading edge, as shown in figure 6. This phenomenon is observed for every wake–rotor interaction event over the course of the simulation and gives rise to the distinct peaks in figure 3. By comparing figure 6 with the forward snapshots in figure 4, the moment of convergence of the adjoint pressure wave at
$t\approx 1.75$
is not commensurate with the earlier appearance of the pressure wave at
$t\approx 1.70$
, suggesting that this sensitivity is not just attempting to energise the initial interaction event, but also a subsequent interaction of the wake vortices with the leading edge, that occurs later at
$t\approx 1.75$
. Indeed, this is further corroborated by observing that the adjoint pressure wave consists of two components: a concentric component, as seen at
$t=1.65$
in figure 6, that disappears at
$t=1.70$
, and a non-concentric component that persists until
$t\approx 1.75$
. In the forward problem of figure 4, the pressure waves from the leading edge are largely concentric, even in the presence of the velocity gradients around the two surfaces of the rotor airfoil.
Based on our earlier observations, we continue by evaluating adjoint data in the vicinity of the leading edge along a near-surface contour defined in the inset of figure 7. It shows that the interaction sensitivities are dominated by a stationary pressure component located at the point where the vortex impacts the pressure surface. In our case, this location, measured in chordwise arclength coordinates, falls at
$s=-0.01$
(with
$s=0$
signifying the stagnation point), which is also near the point of greatest curvature of the surface. In this plot, the acoustic propagation timescale is essentially instantaneous – and so is not observed directly. The convection of the vortex downstream along the boundary layer is captured, however, by a propagating wave seen for
$t \geqslant 1.70.$
Incidentally, the suction surface is also seen to exhibit sensitivities at the mirrored image point
$s=0.01$
, albeit at far lower levels. Consequently, this indicates a more robust and resilient flow structure there. The above observation, confirming a localisation of the dual fields on the pressure surface, also explains the propagation delay seen in the adjoint pressure, shown in figure 6, which cannot be attributed to a different velocity gradient on the two sides of the airfoil.
Figure 7. Waterfall plot showing near-surface adjoint pressure along the red–blue contour wrapped around the leading edge, sketched in the inset. The check mark on the pressure side identifies the chord arclength coordinate
$s=-0.01,$
where the adjoint pressure attains a local maximum.
3.2.2. Interactions from neighbouring blades
Finally, we again concentrate on the acoustic pressure component of the flow and assume that it propagates essentially under a wave equation. A wave equation is self-adjoint, owing to the second derivative in time and space: if we impose a pressure disturbance at some time to optimise an acoustic response at some later time and place, it follows that the form and location of this pressure disturbance must propagate at the acoustic wave speed, and moreover must interact with obstacles in a similar way to the forward problem. Indeed, this can be seen in figure 8, where the adjoint pressure sensitivities discussed earlier are observed to reflect off of neighbouring blade surfaces, in addition to being diffracted and scattered by the leading edges of the rotor stage. The sequence of panels demonstrates the interaction of pressure waves emanating from a neighbouring blade with the suction-side boundary layer of a specific foil, stimulating boundary-layer perturbations that subsequently travel downstream toward the blade’s trailing edge. Starting from figure 8(a), the fully developed boundary-layer instabilities, together with the pressure wave that generated them, are observed. Following this development backward in time via the adjoint equations (see figure 8 b,c), the boundary-layer structures subside and the pressure wave retracts to a concentrated pressure pulse, thus identifying the location on the neighbouring blade and time instant that ultimately cause the triggering of the suction-side boundary layer. We like to stress that a simple backward-in-time approach of the governing equations will not properly identify the same scenario (even though the propagation part of the pressure wave is nearly self-adjoint); rather, the dual equations have to invoked to connect an upstream cause to a downstream effect.
Figure 8. Reflections of dual/adjoint pressure sensitivities generated from neighbouring blades in the rotor section. (a) p †, t = 1.50, (b) p †, t = 1.55, and (c) p †, t = 1.60.
Figure 9. Distortion of adjoint pressure sensitivities by vortices in the wake. (a) p †, t = 1.05, (b) p †, t = 1.07, and (c) p †, t = 1.11.
3.2.3. Acoustic interactions with wake vortices
The overall adjoint pressure dynamics are not just restricted to interactions with solid surfaces, however. For a wavelength comparable to that of the vortex core size, any acoustic pressure disturbance can be partially scattered by the stator wakes, the result of which is a scattering pattern seen at
$t=1.05$
in figure 9(a). In this panel, the dual field identifies localised wavepackets of adjoint pressure that eventually impact the boundary layer of the rotor blade at mid-chord. The acoustic wave that underlies this scattering process emanates from the rotor’s leading edge (see figure 9
c) and interacts with the larger pressure wave from the neighbouring passage. In this scenario, the dual field, while propagating gradient or sensitivity information from the stator, is blocked by the vortices shedding from the stator’s trailing edge. This blockage significantly influences subsequent interactions with the primal flow fields, as only perturbations that avoid a collision with the wake vortices will impact the overall dynamics of the flow.
It is remarkable that the dual field accounts for this interaction/feedback mechanism, even though it does not result in substantially increased sensitivities for the wake vortices themselves, despite the fact that the positioning of the vortices is critical for the scattering properties of the vortex sheet.
3.2.4. Wake–boundary-layer interactions
A final analysis of the complex hydrodynamic–acoustic interaction in our rotor–stator stage concerns the receptivity of the separation bubble on the suction side of the rotor blade to impinging wake vortices from the stator wake.
Solving the dual equation backwards in time and concentrating on the instabilities supported by the separated suction-side region, we observe that select vortices, passing in the vicinity of the suction-side boundary layers of the rotor, carry shear-promoting sensitivities, as shown for a particular interaction event in figure 10; the leading vortex is rotating clockwise. In most cases, these vortices in the forward problem do not merge with the boundary-layer flow around the blade, unless the impact with the leading edge deposits the vortex directly into the boundary layer. Instead, the interaction is more subtle. The global shearing and swirling flow around the blade deflects the wake–vortex pair, behind which the breakup point of the suction-side shear layer of the separation bubble temporarily moves upstream.
Figure 10. Sensitivity of wake vortices when approaching the rotor suction side. Note that these structures untwist as the vortex passes the most sensitive regions located within the separated flow regions. This vortex has a clockwise rotation. (a) s †, t = 1.20, (b) s †, t = 1.25, and (c) s †, t = 1.30.
Links between the vortex and the stability properties of the shear layer are further strengthened by the adjoint, which demonstrates a rapid decay in sensitivity once the targeted wake vortex passes past the separation bubble. Moreover, the shear sensitivity structures surrounding the vortex, in figure 10, untwist as they pass this region. That is to say, the radial wavenumber decreases.
4. Conclusions
The hydro-acoustic interactions within a rotor–stator configuration are an example of a highly complex fluid problem that requires sophisticated tools and formalisms to analyse in detail. In this paper, we have performed an unsteady dual analysis of a model problem that displayed most features of a full realistic configuration. This analysis has been based on a time-domain aeroacoustics code, which is capable of resolving instabilities in the boundary layer and acoustic wave propagation and is additionally able to pass dual/adjoint information cleanly through the sliding mesh interface.
Our simplified model consists of two slender ovate stator blades, placed upstream of three controlled-diffusion rotor blades. The latter blade row moves relative to the stators to ensure a zero angle of attack on the rotor blade from the resulting velocity fields in the rotor reference frame. The flow has been left to evolve until a limit-cycle behaviour was reached, after which this established base flow was seeded with a randomised perturbation field to trigger instabilities.
By casting this layout into a Lagrangian optimisation framework, a set of dual equations and initial conditions are generated that were then integrated backwards in time, using a checkpointing algorithm to account for changes to the tangent-linear operator as the base flow evolves. Through studying the evolution of these dual fields, it was possible to show how flow interactions between the rotor and stator in the forward problem could be traced back to sensitive flow structures which originated at earlier times within the stator boundary layers.
In this way, the acoustic pulses, which occur when the leading edge interacts with the stator wake, can be traced back to Orr-like structures within the stator boundary layer. Further optimisation pathways are also seen by the introduction of short-time-scale adjoint pressure waves which converge on the leading edge during the interaction event, though these appear to be transient rather than the result of an optimisation of a particular physical flow feature.
Observations were also made of stator wake vortices being targeted with shear-promoting perturbations in the region of the laminar separation bubble. Through earlier work, the dynamics of the separation bubble are known to be coupled to the trailing-edge acoustic source, and this link between the stator vortices and the bubble provides evidence that these local instability dynamics and feedback loops on the stator are being affected by the flows from the rotor.
While much has been learned about the specific aeroacoustics within a rotor–stator set-up, the dual analysis presented here also serves as a demonstration for the proper analysis of general fluid problems that are characterised by a complex network of instabilities, omni-directional wave propagation, receptivities and feedback mechanisms. Flow features that emerge under these conditions cannot be analysed by local or classical techniques, but rather require a global viewpoint, to which dual fields contribute valuable sensitivity information. The mutual interplay of individual instability components can be qualitatively and quantitatively assessed and described by tracing dual information of a given feature backward in time to its origin. In this manner, cause-and-effect details about the flow can be gained that support our understanding of the intrinsic flow mechanisms and lay the foundation for performance improvements via passive design modifications or active flow control measures. However, the interpretation of sensitive zones, responsible for the emergence of observable flow features, must be carried out with great circumspection due to the unsteady nature of the underlying process. Concepts from structural sensitivity analysis, such as the wavemaker (also using primal-dual concepts), do not readily generalise to the unsteady case.
While the dual perspective is computationally more involved and, particularly for unsteady flows, like our case, can be challenging to interpret, the added value of possessing a tool that objectively dissects a complex flow into its essential mechanisms can hardly be overestimated. Many complex fluid problems that are driven by a set of complicated, multi-scale and multi-physics interactions are ideal candidates for an analysis based on the techniques applied to our demonstration case.
It is worth mentioning that both the forward (primal) and adjoint (dual) simulations of our rotor–stator flow are rich and highly dynamic in their coherent flow features, and panels of static snapshots cannot do justice to these unsteady motions. To fully appreciate the complexity of the flow and the efficacy of the dual fields, the reader is urged to consult the animations that are available as supplementary material.
Supplementary movie
Supplementary movie is available at https://doi.org/10.1017/jfm.2025.10338.
Declaration of interests
The authors report no conflict of interest.
Open access
