1. Introduction
Dynamic stall is a critical flow phenomenon that significantly limits the performance of unsteady wings due to large-amplitude fluctuations in lift force beyond the static limit (McCroskey Reference McCroskey1981; Carr Reference Carr1988). This phenomenon is intrinsically linked to the unsteady separation of flow over the wing surface, culminating in the formation of a dynamic stall vortex (DSV) (McCroskey, Carr & McAlister Reference McCroskey, Carr and McAlister1976). The DSV is a dominant flow feature of dynamic stall flows, generated by the accumulation of surface vorticity flux concentrated near the leading edge during the wing’s motion (Eldredge & Jones Reference Eldredge and Jones2019). The severe unsteady aerodynamic loads associated with dynamic stall can impose significant performance and structural constraints, underscoring the importance of effective mitigation strategies with active flow control (Karim & Acharya Reference Karim and Acharya1994; Greenblatt & Wygnanski Reference Greenblatt and Wygnanski2001; Weaver, McAlister & Tso Reference Weaver, McAlister and Tso2004; Kissing et al. Reference Kissing, Stumpf, Kriegseis, Hussong and Tropea2021).
While active flow control has demonstrated promising capabilities in mitigating dynamic stall, its design can be challenging due to the highly unsteady and transient nature of the flow (Gardner et al. Reference Gardner, Opitz, Wolf and Merz2017; Zhu et al. Reference Zhu, Qiu, Feng, Wang and Li2022; Gardner et al. Reference Gardner, Jones, Mulleners, Naughton and Smith2023). Particularly, the transient nature of the flow makes it difficult, yet even crucial, to determine an effective timing and location for control inputs to be introduced. The timing of actuation for dynamic stall flows is unclear due to the intricate interplay between the time scale of the flow dynamics and that of the wing motion (Feszty, Gillies & Vezza Reference Feszty, Gillies and Vezza2004; Lombardi, Bowles & Corke Reference Lombardi, Bowles and Corke2013; Tadjfar & Asgari Reference Tadjfar and Asgari2018; de Souza et al. Reference de Souza, Wolf, Safari and Yeh2025). Such interplay also makes actuator placement a daunting question since it causes the separation point to travel over a great extent of the suction surface (Weaver et al. Reference Weaver, McAlister and Tso2004; Post & Corke Reference Post and Corke2006; Singh et al. Reference Singh, Peake, Kokkalis, Khodagolian, Coton and Galbraith2006; Kissing et al. Reference Kissing, Kriegseis, Li, Feng, Hussong and Tropea2020). Therefore, there is a critical need for formal theoretic approaches to determine effective timing and location of actuation (Natarajan, Freund & Bodony Reference Natarajan, Freund and Bodony2016; Nair et al. Reference Nair, Taira, Brunton and Brunton2021; Godavarthi, Kawamura & Taira Reference Godavarthi, Kawamura and Taira2023) for mitigating dynamic stall.
Complex network analysis emerges as an appealing tool to address the need, since analogous questions arise in the field of network science (Albert & Barabási Reference Albert and Barabási2002; Newman Reference Newman2010; Iacobello, Ridolfi & Scarsoglio Reference Iacobello, Ridolfi and Scarsoglio2021). Determining a location for actuator placement is similar to identifying a hub for disease transmission over a human or animal network (Nande et al. Reference Nande, Adlam, Sheen, Levy and Hill2021; Safari et al. Reference Safari, Fleming, Galvis, Deka, Machado and Yeh2025). Modulating the timing of introducing control input is also analogous to determining when to implement preventive measures for disease outbreaks, such as social distancing and isolation (Schlosser et al. Reference Schlosser, Maier, Jack, Hinrichs, Zachariae and Brockmann2020; Vinceti et al. Reference Vinceti, Filippini, Rothman, Ferrari, Goffi, Maffeis and Orsini2020). Therefore, not only does network science offer a broad and mature set of tools to address the aforementioned challenges for flow control (Sujith & Unni Reference Sujith and Unni2020; Yeh, Meena & Taira Reference Yeh, Meena and Taira2021), but its use by the fluid dynamics community has also been witnessed in a growing number of studies, including model reduction (Kaiser et al. Reference Kaiser, Noack, Cordier, Spohn, Segond, Abel, Daviller, Östh, Krajnović and Niven2014), force prediction (Iacobello, Kaiser & Rival Reference Iacobello, Kaiser and Rival2022) and analysing turbulent interactions (Iacobello et al. Reference Iacobello, Scarsoglio, Kuerten and Ridolfi2019; Li et al. Reference Li, Fernex, Semaan, Tan, Morzyński and Noack2021; Tandon & Sujith Reference Tandon and Sujith2023).
Constructing a network model requires the definitions of nodes and edge weights that quantify the interactions between the nodes. Following the classic perspective in aerodynamics where fluid flows can be modelled as an ensemble of discrete vortices (Katz Reference Katz1981; Saffman Reference Saffman1995; Anderson Reference Anderson2016), vortical elements have been adopted by many fluid dynamics studies as the nodes of fluid networks (Lander & Holland Reference Lander and Holland1993; Bai et al. Reference Bai, Erichson, Meena, Taira and Brunton2019; Meena & Taira Reference Meena and Taira2021). To quantify the interactions between the vortical nodes, discrete vortex methods (DVMs) (Cottet & Koumoutsakos Reference Cottet and Koumoutsakos2000; Eldredge Reference Eldredge2007; Eldredge, Wang & Michael Reference Eldredge, Wang and Michael2009) become a natural choice due to their extensive applications to aerodynamic studies (Ramesh et al. Reference Ramesh, Gopalarathnam, Granlund, Ol and Edwards2014; SureshBabu, Ramesh & Gopalarathnam Reference SureshBabu, Ramesh and Gopalarathnam2019; Alvarez & Ning Reference Alvarez and Ning2020; Herndon & Jaworski Reference Herndon and Jaworski2023). In general, DVMs model a flow by time evolving the locations of discrete vortices in a Lagrangian manner. Building on this idea, Ansari, Żbikowski & Knowles (Reference Ansari, Żbikowski and Knowles2006) demonstrated that DVMs can effectively capture the wake dynamics and unsteady aerodynamic forces on an impulsively started flat plate. The work of Eldredge (Reference Eldredge2007) developed a DVM to model the flow around an unsteady body, formulating the approach based on the principles of vorticity generation and diffusion described by Koumoutsakos, Leonard & Pepin (Reference Koumoutsakos, Leonard and Pepin1994). The DVM framework was later improved by Xia & Mohseni (Reference Xia and Mohseni2013), who enforced a Kutta condition at the leading edge for high angles of attack and relaxed this condition at lower angles to accurately model the leading-edge flow over a pitching flat plate. Subsequent improvements have further expanded the applicability of DVMs to more complex flows (Xia & Mohseni Reference Xia and Mohseni2017; Darakananda & Eldredge Reference Darakananda and Eldredge2019; Dumoulin, Eldredge & Chatelain Reference Dumoulin, Eldredge and Chatelain2023).
In general, DVMs model unsteady aerodynamic flows by introducing new vortical elements into the flow while determining their circulation strength upon formation (Clements Reference Clements1973; Sarpkaya Reference Sarpkaya1975; Ramesh et al. Reference Ramesh, Gopalarathnam, Granlund, Ol and Edwards2014). Such introduction of finite-circulation vortices into the flow can be translated to vorticity injection that accounts for the formation of a DSV. Thus, in the context of DVMs, the objective of dynamic stall mitigation can be realised by modulating the strength of discrete vortices introduced to the flow, and such modulation can be achieved by perturbing the vortical elements present in the flow. Also, through an objective-informed combination of the DVM and vortical network analysis, the effective timing and location of actuation can be identified by finding a critical set of vortical nodes at which an introduced perturbation at the identified timing can result in a high level of modification in vorticity injection.
In this study we present a DVM-based network analysis to guide the timing, location and direction of actuation for mitigating dynamic stall. To this end, we model the dynamic stall flow as a time-varying network of discrete vortices and use a DVM (Ramesh et al. Reference Ramesh, Gopalarathnam, Granlund, Ol and Edwards2014; Lee et al. Reference Lee, Suresh Babu, Bryant and Gopalarathnam2022) to extract the interactions between the vortical nodes, as shown in figure 1. Motivated by the objective of flow control, we quantify these vortical interactions (edge weights) by tracking the spread of displacement perturbations between the vortical elements with the use of the DVM. Moreover, this integration of the DVM and network analysis allows us to recast the questions of actuation timing, direction and actuator location into an optimisation problem that seeks an optimal set of vortical nodes to perturb and the timing to perturb them in a direction such that vorticity injection is minimised during dynamic stall and, hence, lift fluctuations can be reduced.
Figure 1. An overview of the present study: (a) a flow undergoing dynamic stall is modelled by a DVM (Ramesh et al. Reference Ramesh, Gopalarathnam, Granlund, Ol and Edwards2014); (b) the vortical elements are treated as network nodes to form a vortical network; (c) their interactions are extracted from the DVM as the edge weights to enable network analysis for mitigating dynamic stall. The vortices shown in (b) comprise of leading-edge vortices (LEVs) and trailing-edge vortices (TEVs).
With such an integration of the DVM and network analysis, this study extends the broadcast mode analysis of Yeh et al. (Reference Yeh, Meena and Taira2021) in several different aspects. While Yeh et al. (Reference Yeh, Meena and Taira2021) considered a set of Eulerian grid points as the vortical nodes, this work directly utilises the Lagrangian vortical elements in the DVM as the nodes of the network. Moreover, motivated by the objective of flow control, the present study quantifies the interactions between the vortical nodes by tracking the spread of displacement perturbation among them. This edge-weight definition distinguishes the present study from Yeh et al. (Reference Yeh, Meena and Taira2021) and many previous ones (Nair & Taira Reference Nair and Taira2015; De et al. Reference De, Gupta, Unni, Ravindran, Kasthuri, Marwan, Kurths and Sujith2023; Wang et al. Reference Wang, Zhang, Sun and Zhou2025) that adopted the induced velocity between vortical elements as the edge weights of the vortical networks, according to the Biot–Savart law. The broadcast mode analysis by Yeh et al. (Reference Yeh, Meena and Taira2021) used a singular value decomposition to seek the maximum amplification of vortical perturbations for modifying the evolution of two-dimensional (2-D) turbulence by targeting dynamically influential regions. Compared with Yeh et al. (Reference Yeh, Meena and Taira2021), this work generalises the analysis by formulating an optimisation problem that seeks a vortical perturbation not only for high amplification, but also in a manner that the identified perturbation maximises the reduction of vorticity injection, given the objective of reducing lift fluctuation by weakening the DSV. Combining with the control-motivated edge weights, such a generalisation of broadcast mode analysis provides insights into the optimal direction of actuation, in addition to identifying dynamically influential regions in the flow for actuator placement. Moreover, it also allows for the identification of effective timing of actuation by considering different time horizons over which the maximum reduction of vorticity injection is sought.
In what follows, we discuss the problem set-up for the dynamic stall flows in § 2. The construction of the DVM-based network model with the goal of guiding circulation control is elaborated in § 3. In § 3.3 we present the broadcast mode analysis (Grindrod & Higham Reference Grindrod and Higham2014; Yeh et al. Reference Yeh, Meena and Taira2021) to identify the optimal location and timing for actuation with the objective of dynamic stall mitigation by modulating vorticity injection. Section 4 leverages the insights from the analysis to inform flow control design and demonstrates the effectiveness of the network-based approach for mitigating dynamic stall.
2. Problem set-up
We consider two incompressible flows undergoing dynamic stall, as shown in figure 2. The first case (case A) is a flow over a 2-D flat plate aerofoil in pitch-up motion at a chord-based Reynolds number,
$ \textit{Re} \equiv u_{\infty }L_{c}/\nu _{\infty } = 1\,000$
. Here,
$u_{\infty }$
is the free-stream velocity,
$L_{c}$
is the chord length and
$\nu _{\infty }$
is the kinematic viscosity. For this case of pitching aerofoil, the instantaneous angle of attack is prescribed via
\begin{align} {\alpha }(t) = \alpha _0 + \frac {\alpha _f}{b} \ln \left [ \frac {\cosh (a (t - t_1))} {\cosh (a (t - t_2))}\right ]\!, \end{align}
following Eldredge et al. (Reference Eldredge, Wang and Michael2009) and Ramesh et al. (Reference Ramesh, Gopalarathnam, Granlund, Ol and Edwards2014). We consider this pitch-up and hold motion, since the onset of dynamic stall occurs during the pitch-up process, and the flow exhibits rich physics involving the formation and convection of the DSV. The second case (case B) is a spanwise-periodic three-dimensional (3-D) turbulent flow over an SD
$7003$
aerofoil undergoing periodic heaving-plunging motion at
$ \textit{Re} = 60\,000$
, where the instantaneous vertical velocity of the aerofoil is given by
\begin{align} \dot {h}(t) = -v_0 \sin \big (2\pi t / T_{p}\big ), \end{align}
following Visbal (Reference Visbal2011), Ramos et al. (Reference Ramos, Wolf, Yeh and Taira2019). The motion parameters in (2.1) and (2.2) are provided in the caption of figure 2, along with the motion profiles, instantaneous lift coefficients and flow visualisations at instants of the DSV departure for both cases.
Figure 2. The two dynamic stall flows considered in this study: (a) a flow over a pitching flat plate aerofoil (case A) – the red circles (
) are the results by Ramesh et al. (Reference Ramesh, Gopalarathnam, Granlund, Ol and Edwards2014); (b) flow over a periodic heaving-plunging SD
$7003$
aerofoil (case B) – the red circles () are the results by Visbal (Reference Visbal2011). Contour lines of spanwise vorticity
$\omega _{z}L_{c}/u_{\infty } \in [-20, 20]$
is shown for case (a). Isosurface of
$QL^{2}_{c}/u^{2}_{\infty } = 80$
coloured by spanwise vorticity is shown for case (b). Motion parameters in (2.1) for case A:
$\alpha _0 = \alpha _f = 22.5^\circ$
,
$aL_c/u_\infty = 11$
,
$b = 10.78$
,
$[t_1,\,t_2]/(L_c/u_\infty ) = [1,\,1.98]$
; motion parameters in (2.2) for case B:
$\nu _{0}/u_{\infty } = 0.25$
,
$T_{p}u_{\infty }/L_{c} = 4\pi$
. The DVM states at key instants are overlaid on the contour lines of spanwise vorticity for comparison. In the DVM, blue circles (
) indicate LEVs and red circles (
) indicate TEVs. The close-up view for these key instances are presented to highlight the comparison of leading-edge flow between DVM and CFD. For the DVM states, the circulation strength of discrete vortices is indicated by transparency – the lower the circulation strength, the greater the transparency, and vice versa.
Baseline (uncontrolled) flows of both cases are obtained using large-eddy simulations (LES) and the DVM. The LES are performed using the incompressible flow solver package, CharLES (Khalighi et al. Reference Khalighi, Ham, Nichols, Lele and Moin2011; Brès et al. Reference Brès, Ham, Nichols and Lele2017). This solver uses a node-based, second-order finite-volume method for spatial discretisation and a fractional-step scheme for time stepping. To account for the unsteady motions of the aerofoil, a time-dependent fictitious forcing is added to the right-hand side of the momentum equation (Kim & Choi Reference Kim and Choi2006). Further details about the LES solver, governing equations and computational meshes are provided in Appendix A.
The DVM developed by Ramesh et al. (Reference Ramesh, Gopalarathnam, Granlund, Ol and Edwards2014) and Lee et al. (Reference Lee, Suresh Babu, Bryant and Gopalarathnam2022) is employed in the present study. It models unsteady flows over aerofoils by introducing discrete point vortices into the flow, whose circulations are modulated by an empirical leading-edge suction parameter. The time integration of the DVM adopts a fourth-order Runge–Kutta scheme, with a time step
$\Delta t u_{\infty }/L_{c} = 0.01$
. The comparisons in lift obtained from the DVM and LES are also shown in figure 2. Although over-predictions in the DVM lift are seen in both cases, we observe qualitative agreements in their temporal profiles. Quantitatively, the levels of agreement in lift between the DVM and computational fluid dynamics (CFD) are also similar to those reported in previous studies (Ramesh et al. Reference Ramesh, Gopalarathnam, Granlund, Ol and Edwards2014; Lee et al. Reference Lee, Suresh Babu, Bryant and Gopalarathnam2022). The DVM states at representative instants are also shown on top of the spanwise vorticity fields obtained from the flow simulations. We observe that the onset of the DSV core formation, which occurs at
$tu_{\infty }/L_{c} = 1.35$
for the flat plate aerofoil and
$t/T_{p} = 0.21$
for the SD
$7003$
aerofoil, is captured in agreement with the flow simulations. For the plunging SD
$7003$
case, the over-prediction in lift by the DVM for
$t/T_{p} \gt 0.25$
can be attributed to the DSV that is located closer to the aerofoil, compared with that in the flow simulation. The results of the DVM base flow will be used to construct the network model of discrete vortices, which is to be detailed in the following section.
3. Construction of the discrete vortex network
Our goal of mitigating dynamic stall is approached by modelling the stall flows as networks of vortical elements, which enables a control-focused network analysis to provide guidance to effective actuation timing, location and direction. The first step of the network analysis is to construct an adjacency matrix,
$\unicode{x1D63C}$
, with each of its entries,
$A_{\textit{ij}}$
, representing a quantifiable interaction (edge weight) between a pair of vortical nodes
$i$
and
$j$
. Below, we discuss how these vortical interactions are quantified with the use of a DVM.
3.1. The discrete vortex method
Here, we provide a concise overview of the DVM developed by Ramesh et al. (Reference Ramesh, Gopalarathnam, Granlund, Ol and Edwards2014) and Lee et al. (Reference Lee, Suresh Babu, Bryant and Gopalarathnam2022), which is employed in this study to quantify the interactions between the vortical nodes. Although some details below may be specific to the adopted DVM, our discussion will emphasise how it facilitates the construction of the vortical network in a manner that the general concept is applicable to any DVMs (Cottet & Koumoutsakos Reference Cottet and Koumoutsakos2000; Eldredge Reference Eldredge2007). For more details about the DVM, the readers are referred to the original works by Ramesh et al. (Reference Ramesh, Gopalarathnam, Granlund, Ol and Edwards2014) and Lee et al. (Reference Lee, Suresh Babu, Bryant and Gopalarathnam2022).
Figure 3. The flow state in DVM is illustrated in (a) with arrangement of bound, trailing-edge and leading-edge vortical elements and its comparison with an actual flow simulation is presented in (b).
The DVM simulates 2-D flows over an aerofoil using three types of discrete vortices: leading-edge vortices (LEVs), trailing-edge vortices (TEVs) and bound vortices (BVs), as illustrated in figure 3. The BVs are distributed along the mean camber line of the aerofoil and represent its location. The LEVs represent the flow separated from the leading edge of the aerofoil and the TEVs capture the flow past the trailing edge. Together, LEVs and TEVs are referred to as free vortices (FVs), since their locations vary in time. Here, we denote the locations and circulations of all FVs at time
$t = t_{k}$
by
\begin{align} \boldsymbol{\varXi }_{k}^{\textit{f}v} = \left[ \boldsymbol{\xi }_{1, k}^{\textit{f}v},\, \boldsymbol{\xi }_{2, k}^{\textit{f}v},\,\ldots,\, \boldsymbol{\xi }_{n_{\textit{f}v}, k}^{\textit{f}v} \right] \quad \text{and} \quad \boldsymbol{\varGamma }^{\textit{f}v}_{k} = \left[ \gamma ^{\textit{f}v}_{1},\, \gamma ^{\textit{f}v}_{2},\,\ldots,\, \gamma ^{\textit{f}v}_{n_{\textit{f}v}} \right]\!, \end{align}
where
$n_{\textit{f}v}$
is the number of FVs present in the flow at
$t_k$
. Note that
$n_{\textit{f}v}$
increases in time, i.e.
$n_{\textit{f}v} = n_{\textit{f}v}(t)$
, since the DVM introduces new FVs to the flow at every time step. Also, once a FV is introduced to the flow, its circulation is held constant in time; therefore, the time index
$k$
is not needed for
$\gamma ^{\textit{f}v}$
values. Similarly, the circulations of the BVs at
$t = t_{k}$
are denoted by
\begin{align} \boldsymbol{\varGamma }^{{bv}}_{k} = \left[ \gamma ^{{bv}}_{1, k},\, \gamma ^{{bv}}_{2, k},\,\ldots,\, \gamma ^{{bv}}_{n_{{bv}}, k} \right]\!, \end{align}
where
$n_{{bv}}$
is the number of BVs distributed along the aerofoil. We also note that the locations of the BVs,
$[\boldsymbol{\xi }^{{bv}}_{1, k},\,\boldsymbol{\xi }^{{bv}}_{2, k},\,\ldots,\,\boldsymbol{\xi }^{{bv}}_{n_{{bv}}, k}]$
, are given by a prescribed motion of the aerofoil and, hence, are not the state variables to be determined by the DVM.
For each time step from
$t_{k}$
, the DVM seeds a TEV and a LEV, respectively, at the trailing edge (
$\boldsymbol{\xi }^{\textit{te}v}_k$
) and leading edge (
$\boldsymbol{\xi }^{\textit{le}v}_k$
). (In practice, the introduction of a LEV at a given instant is conditioned by the instantaneous leading-edge suction parameter (Ramesh et al. Reference Ramesh, Gopalarathnam, Granlund, Ol and Edwards2014). However, for an instant where a LEV is not introduced, this can be alternatively viewed as introducing a null LEV with zero circulation without loss of generality.) The circulations of the newly seeded TEV and LEV, denoted respectively as
$\gamma ^{\textit{te}v}_{k+1}$
and
$\gamma ^{\textit{le}v}_{k+1}$
, are determined along with the circulations of all BVs,
$\boldsymbol{\varGamma }^{{bv}}_{k+1}$
, by imposing the no-penetration condition over the aerofoil and Kelvin’s circulation theorem. The new TEV and LEV are appended to the state vectors of FVs such that
$\boldsymbol{\varXi }^{\textit{f}v}_{k'} = [\boldsymbol{\varXi }_{k}^{\textit{f}v},\, \boldsymbol{\xi }^{\textit{te}v}_k,\, \boldsymbol{\xi }^{\textit{le}v}_k]$
and
$\boldsymbol{\varGamma }^{\textit{f}v}_{k'} = [\boldsymbol{\varGamma }^{\textit{f}v}_{k},\, \gamma ^{\textit{te}v}_{k+1},\, \gamma ^{\textit{le}v}_{k+1}]$
, thereby increasing the number of FVs by 2. With the circulations for all BVs and new FVs determined, the DVM then updates the locations of all FVs in
$\boldsymbol{\varXi }^{\textit{f}v}_{k'}$
by
\begin{align} \boldsymbol{\xi }^{\textit{f}v}_{i, k+1} = \boldsymbol{\xi }^{\textit{f}v}_{i, k} + \Delta t \boldsymbol{\dot {\xi }}^{\textit{f}v}_{i, k} \quad \text{for all }i = 1,\,2,\,\ldots ,\,n_{\textit{f}v}(t_{k})+2, \end{align}
where
$\Delta t$
is the time increment between
$t_{k}$
and
$t_{k+1}$
, and
$\boldsymbol{\dot {\xi }}^{\textit{f}v}_{i, k}$
is the velocity of the
$i$
th FV induced by all other discrete vortices. This induced velocity is given by
\begin{align} \begin{split} \boldsymbol{\dot {\xi }}^{\textit{f}v}_{i, k} = \sum ^{n_{\textit{f}v}+2}_{j=1} \mathcal{V}_{{bs}}\left (\boldsymbol{\xi }^{\textit{f}v}_{i, k},\,\boldsymbol{\xi }^{\textit{f}v}_{j, k},\,\gamma ^{\textit{f}v}_{j}\right ) &+ \sum ^{n_{{bv}}}_{j=1} \mathcal{V}_{{bs}}\left (\boldsymbol{\xi }^{\textit{f}v}_{i, k},\, \boldsymbol{\xi }^{{bv}}_{j, k+1},\,\gamma ^{{bv}}_{j, k+1}\right )\!, \end{split} \end{align}
where
$\mathcal{V}_{{bs}}$
is the Biot–Savart operator written as
\begin{align} \mathcal{V}_{{bs}}\left (\boldsymbol{\xi }_{i},\,\boldsymbol{\xi }_{j},\,\gamma _{j}\right ) = \frac {\gamma _{j}}{2\pi } \frac {\boldsymbol{\hat {e}}_{z} \times (\boldsymbol{\xi }_{i} - \boldsymbol{\xi }_{j})}{|\boldsymbol{\xi }_{i} - \boldsymbol{\xi }_{j}|^{2} + r^{2}_{c}}, \end{align}
with the core size
$r_{c} = 1.3 \Delta t u_{\infty }$
(Lee et al. Reference Lee, Suresh Babu, Bryant and Gopalarathnam2022). Once the locations of all FVs are updated, we obtain
$ \boldsymbol{\varXi }^{\textit{f}v}_{k+1} = [ \boldsymbol{\xi }^{\textit{f}v}_{1, k+1},\, \boldsymbol{\xi }^{\textit{f}v}_{2, k+1},\, \ldots ,\, \boldsymbol{\xi }^{\textit{f}v}_{n_{{fv}(t_k)}, k+1},\, \boldsymbol{\xi }^{\textit{te}v}_{k+1},\, \boldsymbol{\xi }^{\textit{le}v}_{k+1}]$
, and the time iteration of the DVM is concluded with
$\boldsymbol{\varXi }^{\textit{f}v}_{k+1}$
,
$ \boldsymbol{\varGamma }^{\textit{f}v}_{k+1} = \boldsymbol{\varGamma }^{\textit{f}v}_{k'} = [\boldsymbol{\varGamma }^{\textit{f}v}_{k},\,\gamma ^{\textit{te}v}_{k+1},\,\gamma ^{\textit{le}v}_{k+1}]$
and
$\boldsymbol{\varGamma }^{{bv}}_{k+1}$
. Summarising the discussion above, we can express the time-stepping operation of the DVM in brief as
\begin{align} \left [ \boldsymbol{\varXi }^{\textit{f}v}_{k+1},\,\gamma ^{\textit{te}v}_{k+1},\,\gamma ^{\textit{le}v}_{k+1},\, \boldsymbol{\varGamma }^{{bv}}_{k+1} \right ] = \mathcal{N} \left (\boldsymbol{\varXi }^{\textit{f}v}_{k},\,\boldsymbol{\varGamma }^{\textit{f}v}_{k}\right )\!, \end{align}
where
$\mathcal{N}$
can be viewed as the discrete-time nonlinear state propagator of the DVM, whose nonlinearity stems from that in the Biot–Savart law.
Observing (3.6), we make a few important notes on how the DVM can be integrated with network analysis to strategise dynamic stall mitigation: first, (3.6) suggests that the circulation introduced from the leading edge (
$\gamma ^{\textit{le}v}_{k+1}$
) can be modified by perturbing the locations of FVs (
$\boldsymbol{\varXi }^{\textit{f}v}_{k}$
) present in the flow; second, if the location of a FV is perturbed at a time instant, the locations of all FVs will be perturbed at all future instants successively, resulting in continuous modification in the circulation introduced from the leading edge. Since this circulation introduction can be translated to a concentrated vorticity injection from the leading edge, which results in the formation of a DSV, (3.6) also implies that dynamic stall mitigation can be potentially achieved by perturbing the FVs in the flow. Thus, the objective of the network analysis is to identify a set of FVs and an effective timing at which introduced perturbations to their locations can result in favourable modifications to the circulations of LEVs. Below, we discuss how this objective can be achieved by appropriately quantifying the interactions between the vortical elements.
3.2. The time-evolving vortical network: adjacency and communicability matrices
Motivated by the observations on (3.6), we treat FVs in the DVM as the nodes in our network analysis. To quantify the interactions between them, we focus on those that take place due to displacement perturbations introduced to the vortical nodes, since such perturbations result in a continuous change in the circulation introduced from the leading edge.
In this study we construct the adjacency matrix by measuring the displacement perturbation passed between a pair of vortical nodes. Given a priori an instantaneous FV state at
$t_{k}$
where
$\boldsymbol{\varXi }^{\textit{f}v}_{k}$
and
$\boldsymbol{\varGamma }^{\textit{f}v}_{k}$
are known, we define the edge weight,
$A_{\textit{ij}, k}$
, as the resulting displacement of vortical node
$i$
over a time increment due to a small displacement perturbation introduced to vortical node
$j$
. That is,
\begin{align} A_{\textit{ij}, k} = \frac {\partial \boldsymbol{\xi }_{i, k+1}^{\textit{f}v}}{\partial \boldsymbol{\xi }_{j, k}^{\textit{f}v}}. \end{align}
Using the state propagator of the DVM (3.6), this edge weight can be computed by
\begin{align} A_{\textit{ij}, k} = \frac {1}{\epsilon } \left .\left [\mathcal{N} \left (\boldsymbol{\varXi }^{\textit{f}v}_{k} + \epsilon \boldsymbol{\varXi }^{\prime }(j),\,\boldsymbol{\varGamma }^{\textit{f}v}_{k}\right ) -\mathcal{N} \left (\boldsymbol{\varXi }^{\textit{f}v}_{k},\,\boldsymbol{\varGamma }^{\textit{f}v}_{k}\right )\right ] \right \rvert _{\boldsymbol{\xi }_{i, k+1}^{\textit{f}v}}\!, \end{align}
where
$\epsilon \boldsymbol{\varXi }^{\prime }(j)$
introduces an
$\boldsymbol{\epsilon }$
perturbation to the location of the
$j$
th FV and zeros for all other FVs. The calculation for edge weights is demonstrated in figure 4. The translation of the
$i$
th FV over a time step is changed from
$\boldsymbol{\xi }^{\textit{f}v}_{i, k+1}$
to
$\boldsymbol{\xi '}^{\textit{f}v}_{i, k+1}$
due to an
$\boldsymbol{\epsilon }$
perturbation that displaces
$\boldsymbol{\xi }^{\textit{f}v}_{j, k}$
to
$\boldsymbol{\xi '}^{\textit{f}v}_{j, k}$
, as shown in figure 4(b). Thus, following (3.8), the edge weight,
$A_{\textit{ij}, k}$
, is computed by
\begin{align} A_{\textit{ij}, k} = \frac {1}{|\boldsymbol{\epsilon }|}\left (\boldsymbol{\xi '}^{\textit{f}v}_{i, k+1} - \boldsymbol{\xi }^{\textit{f}v}_{i, k+1}\right )\!. \end{align}
The value of
$|\boldsymbol{\epsilon }|$
is set to
$10^{-7} r_{c}$
, where linear convergence of each
$A_{\textit{ij}}$
is observed. We note that, since both
$\boldsymbol{\epsilon }$
and
$\boldsymbol{\xi }$
are vectors, the resulting
$A_{\textit{ij}, k}$
is in fact a 2-by-2 matrix for the 2-D setting of the present study. More specifically, the
$\boldsymbol{\epsilon }$
perturbation introduced to the
$j$
th FV is applied independently in the
$x$
and
$y$
directions. The
$x$
-direction perturbation introduced to the
$j$
th FV will induce a displacement of the
$i$
th FV in both the
$x$
and
$y$
directions. Likewise, perturbing the same
$j$
th FV in the
$y$
direction yields the corresponding displacement of the
$i$
th FV in both directions. For conceptual simplicity, we refer to the direction-embedded pairwise interaction as the edge weight. Once the edge weights for all pairs of vortical nodes,
$A_{\textit{ij}, k}$
, are computed, the adjacency matrix,
$\unicode{x1D63C}_{k} \in \mathbb{R}^{2n_{\textit{f}v}\times 2n_{\textit{f}v}}$
, can be fully constructed with a given
$\boldsymbol{\varXi }^{\textit{f}v}_{k}$
and
$\boldsymbol{\varGamma }^{\textit{f}v}_{k}$
at
$t = t_{k}$
.
Figure 4. Demonstration of the edge-weight calculation: (a) the translations of the
$i$
th and
$j$
th FVs over a time step without perturbations; (b) an
$\boldsymbol{\epsilon }$
perturbation is introduced to displace the
$j$
th FV and changes the translations of the
$i$
th and
$j$
th FVs over a time step; (c) the edge weight,
$A_{\textit{ij}, k}$
, is obtained by calculating the differences between the perturbed and unperturbed translations for all FVs.
Note that the locations of FVs,
$\boldsymbol{\varXi }^{\textit{f}v}_{k}$
, vary in time. This results in a time-varying vortical network, where the edge weights and the constructed
$\unicode{x1D63C}_{k}$
become time dependent. In order to track the perturbation spread over such a network, a communicability matrix is required to account for the time-varying interactions between the vortical nodes (Grindrod & Higham Reference Grindrod and Higham2013, Reference Grindrod and Higham2014). The communicability matrix,
$\unicode{x1D64E}(t_{n}, t_{m})$
, describes the interactions between the network nodes that occur over a finite-time horizon from
$t_{n}$
to
$t_{m}$
. To construct
$\unicode{x1D64E}(t_{n}, t_{m})$
, we first collect a series of snapshots for
$\boldsymbol{\varXi }^{\textit{f}v}_{k}$
and
$\boldsymbol{\varGamma }^{\textit{f}v}_{k}$
at each time step of the DVM over the time interval
$t \in [t_n, t_m]$
. These snapshots result in a series of adjacency matrices,
$\unicode{x1D63C}_{k}$
, where
$k$
is the time index of the FV state about which the adjacency matrix is constructed. Then, we follow Grindrod et al. (Reference Grindrod, Parsons, Higham and Estrada2011) and Grindrod & Higham (Reference Grindrod and Higham2013) to construct the communicability matrix
$\unicode{x1D64E}(t_{n}, t_{m})$
by recurrently using
\begin{align} \unicode{x1D64E}(t_{n}, t_{k+1}) = \big [\boldsymbol{I} + \unicode{x1D63C}_{k}\big ]\big [\boldsymbol{I} + e^{-1} \big(\unicode{x1D64E}(t_{n}, t_{k}) - \boldsymbol{I} \big)\big ], \end{align}
for
$k = n,\,n+1,\,\ldots ,\,m-1$
, with
$\unicode{x1D64E}(t_{n}, t_{n}) = \boldsymbol{I}$
, the identity matrix. In (3.10) the time-downweighting factor,
$e^{-1}$
, enforces a finite-time decay of vortical perturbations and suppresses unbounded growth of their magnitudes. This is conceptually similar to the exponential discounting in resolvent analysis (Jovanović Reference Jovanović2004; Yeh et al. Reference Yeh, Benton, Taira and Garmann2020) and can be motivated by the flow hysteresis during dynamic stall, where the memory effect of wing motion eventually diminishes (Corke & Thomas Reference Corke and Thomas2015). For the present vortical perturbation network, each entry of the communicability matrix, or
$S_{\textit{ij}}$
of
$\unicode{x1D64E}(t_{n}, t_{m})$
, represents the resulting displacement of the
$i$
th vortical node at
$t_{m}$
due to the displacement introduced at the
$j$
th vortical node at an earlier time
$t_{n}$
. That is, given a FV displacement perturbation at
$t_n$
,
$\boldsymbol{\varXi }'_n$
, the resulting displacement at
$t_m$
,
$\boldsymbol{\varXi }'_m$
is related to
$\boldsymbol{\varXi }'_n$
by
$\unicode{x1D64E}(t_{n}, t_{m})$
through
\begin{align} \boldsymbol{\varXi }'_m = \unicode{x1D64E}(t_n, t_m)\boldsymbol{\varXi }'_n. \end{align}
Thus, the communicability matrix,
$\unicode{x1D64E}(t_n, t_m)$
, can be viewed as a state-transition operator of the displacement perturbation of FVs,
$\boldsymbol{\varXi }'$
, over a finite-time horizon between
$t_n$
and
$t_m$
.
Figure 5. The adjacency and communicability matrices constructed about and between three representative time instances. (a–c) Locations of FVs at each instance. (d–f) Adjacency matrices constructed about each instance. (g–i) Communicability matrices constructed between the combinations of the instances. For visual clarity, the matrices in (d–i) are subsampled such that only the interactions between the vortical elements highlighted by solid colours in (a–c) are shown.
Let us examine in figure 5 the adjacency and communicability matrices constructed at and between a few representative time instants, respectively. These are accompanied by the instantaneous
$\boldsymbol{\varXi }^{\textit{f}v}$
for the same set of instants. For the adjacency matrices at earlier times, we zero pad the rows and columns corresponding to the vortical nodes that have yet to be seeded into the flow. Also, each matrix is dissected into four blocks according to either TEVs or LEVs the rows and columns are associated with. At the earliest instant
$t_1$
, where the LEV introduction is just initiated by the DVM, we observe high levels of intra-community interactions within the TEV and LEV groups from the diagonal blocks of
$\unicode{x1D63C}(t_{1})$
. Meanwhile, the off-diagonal blocks reveal weaker inter-community interactions between TEVs and LEVs, which is attributed to the longer distances between the TEVs and LEVs. At
$t = t_2$
, more LEVs and TEVs are introduced to the flow. Their interactions populate the inactive entries of
$\unicode{x1D63C}(t_1)$
and emerge in the off-diagonal blocks of
$\unicode{x1D63C}(t_2)$
. Notably, while the interactions in the TEVs-to-TEVs block remain limited to the entries near the diagonal (close neighbours of vortical elements), we observe significantly richer interactions within the LEV community, which can potentially be leveraged to mitigate dynamic stall. Similar observations can be made for
$\unicode{x1D63C}(t_3)$
, at which instant the DVM state is characterised by the formation of a DSV. At this instant, we also observe that a LEV has advected to the trailing edge, which results in high levels of interactions between this LEV and the TEVs just seeded into the flow.
Next, let us move our attention to the communicability matrices constructed between these three representative time instants. Since each
$S_{\textit{ij}}$
element of
$\unicode{x1D64E}(t_{n}, t_{m})$
represents the displacement of the
$i$
th vortex at
$t_{m}$
due to a perturbation introduced to the
$j$
th vortex at
$t_{n}$
, the
$j$
th column of
$\unicode{x1D64E}(t_{n}, t_{m})$
represents the resulting displacements of all the vortices present at
$t_{m}$
due to the perturbation of the
$j$
th vortex at
$t_{n}$
. Therefore, the columns of non-zero elements in
$\unicode{x1D64E}(t_{n}, t_{m})$
are the same as those in
$\unicode{x1D63C}(t_{n})$
, and the rows of non-zero elements in
$\unicode{x1D64E}(t_{n}, t_{m})$
are the same as those in
$\unicode{x1D63C}(t_{m})$
. Such a structure of
$\unicode{x1D64E}(t_{n}, t_{m})$
suggests that only the vortical elements present at
$t_n$
can be perturbed, and the perturbation can spread to all vortical elements present at
$t_m$
. Comparing
$\unicode{x1D64E}(t_{1}, t_{2})$
with
$\unicode{x1D63C}(t_{1})$
, we observe higher levels of interactions in the off-diagonal entries for both intra- and inter-community interactions within and between the TEV and LEV groups, respectively. The same observations can be made by comparing
$\unicode{x1D64E}(t_{1}, t_{3})$
with
$\unicode{x1D63C}(t_{1})$
and
$\unicode{x1D64E}(t_{2}, t_{3})$
with
$\unicode{x1D63C}(t_{2})$
, which is attributed to the spread of perturbations over time. Similar to
$\unicode{x1D63C}$
, stronger intra-community interactions within the LEV group are observed for
$\unicode{x1D64E}$
. Most importantly, all communicability matrices show stronger LEVs-to-TEVs interactions than TEVs-to-LEVs interactions. This implies that perturbing LEVs is more effective in resulting in a global modification of discrete vortex dynamics compared with perturbing TEVs.
3.3. The network broadcast analysis
Recall the remarks we made on the state propagator of the DVM (3.6): it suggests that dynamic stall mitigation can potentially be achieved by identifying a set of FVs and an effective timing at which introduced perturbations result in desirable modifications to the LEV circulations. This motivates the use of network broadcast analysis (Grindrod & Higham Reference Grindrod and Higham2014; Yeh et al. Reference Yeh, Meena and Taira2021) to identify the optimal timing and the corresponding set of FVs at which introduced perturbations effectively reduce the circulations of the LEVs. In this analysis, the capability of each node in spreading information over a time-varying network is quantified by a broadcast mode, and that of each node in collecting information at a later time is distilled into a receiving mode. Traditionally, the broadcast/receiving modes are obtained from either a column/row sum or a singular value decomposition of the communicability matrix (Grindrod & Higham Reference Grindrod and Higham2014; Yeh et al. Reference Yeh, Meena and Taira2021). Here, we customise its formulation by focusing on how introduced vortical perturbations result in desirable changes of the flow.
In this study we generalise the network broadcast analysis by formulating an optimisation problem that seeks a set of vortical elements where introduced perturbations minimise the circulation generation from the leading edge over a specified time horizon. Recall that the communicability matrix,
$\unicode{x1D64E} (t_n, t_k)$
, serves as a finite-time state propagator for a displacement perturbation,
$\boldsymbol{\varXi }'_n$
, over
$t \in [t_n, t_k]$
, as discussed in (3.11). The resulting displacement perturbation at
$t_k$
, or
$\boldsymbol{\varXi }'_k$
, will result in a change in circulation of the new LEV introduced at
$t_{k+1}$
. Therefore, we can relate the change of LEV circulation at
$t_{k+1}$
, or
$\Delta \gamma ^{\textit{le}v}_{k+1}$
, to the displacement perturbation introduced at
$t_n$
via a state-space representation, written as
\begin{align} \boldsymbol{\varXi }'_k &= \unicode{x1D64E} (t_n, t_k) \boldsymbol{\varXi }'_n, \\[-12pt] \nonumber \end{align}
\begin{align} \Delta \gamma ^{\textit{le}v}_{k+1} &= \sum _{i=1}^{n_{\textit{f}v}(t_k)}\unicode{x1D63E}_k \boldsymbol{\varXi }'_k, \\[9pt] \nonumber \end{align}
where
$\unicode{x1D63E}_{k}$
is a diagonal matrix in which each of its diagonal elements,
\begin{align} C_{\textit{ii}, k} = -\text{sgn}({\gamma ^{\textit{le}v}_{k+1}})\frac {\partial \gamma ^{\textit{le}v}_{k+1}}{\partial \boldsymbol{\xi }^{\textit{f}v}_{i, k}}, \end{align}
quantifies the reduction of LEV circulation due to a displacement perturbation introduced at the
$i$
th vortical node at
$t_k$
. Using the DVM state propagator in (3.6), we can compute these diagonal elements via
\begin{align} C_{\textit{ii}, k} = -\frac {\text{sgn}({\gamma ^{\textit{le}v}_{k+1}})}{\epsilon } \left .\left [\mathcal{N} \left (\boldsymbol{\varXi }^{\textit{f}v}_{k} + \epsilon \boldsymbol{\varXi }^{\prime }(i),\,\boldsymbol{\varGamma }^{\textit{f}v}_{k}\right ) - \mathcal{N} \left (\boldsymbol{\varXi }^{\textit{f}v}_{k},\,\boldsymbol{\varGamma }^{\textit{f}v}_{k}\right )\right ] \right \rvert _{\gamma _{k+1}^{\textit{le}v}}\!. \end{align}
Revisiting (3.12), now we can generalise the network broadcast analysis by solving an optimisation problem
to seek a broadcast mode,
$\boldsymbol{\varXi }^{{b}}(t_n, t_m)$
, as an optimal displacement perturbation introduced to vortical elements at
$t_n$
that maximises the sum of the reduction in LEV circulations over all later times until
$t_m$
. Meanwhile, the receiving mode,
$\boldsymbol{\varXi }^{\text{r}}(t_n, t_m)$
, is calculated via
\begin{align} \boldsymbol{\varXi }^{\text{r}}(t_n, t_m) = \unicode{x1D63E}_k \unicode{x1D64E}(t_n, t_m) \boldsymbol{\varXi }^{{b}}, \end{align}
which quantifies the receptivity of each vortical element present at
$t_m$
to the perturbation introduced at
$t_n$
in the shape of the broadcast mode,
$\boldsymbol{\varXi }^{{b}}$
, with the objective of resulting in the maximum reduction in LEV circulations. Note that the obtained broadcast mode,
$\boldsymbol{\varXi }^{{b}}(t_n, t_m)$
, is optimal for the chosen time horizon
$t \in [t_n, t_m]$
. Furthermore, the resulting scalar objective function,
$\mathcal{J}(\boldsymbol{\varXi }^{{b}}(t_n, t_m))$
, corresponds to the total reduction in circulations of LEVs introduced to the flow over
$t \in [t_n, t_m]$
. Thus, the optimal timing of introducing perturbations can be identified by the perturbation time
$t_n$
that results in the maximum value of
$\mathcal{J}(\boldsymbol{\varXi }^{{b}}(t_n, t_m))$
for a given
$t_m$
, which can be set to a critical instant associated with the DSV formation. This is demonstrated in figure 6, where the broadcast analysis is conducted between different combinations of
$t_n$
and
$t_m$
, with
$t_m \gt t_n$
. In particular, we sweep over different values of the perturbation time
$t_n$
while keeping
$t_m$
at two representative values, as shown in figure 6(a). For three representative
$t_n$
-
$t_m$
combinations, marked as (i), (ii) and (iii), the broadcast and receiving modes are visualised by colouring the FVs present at
$t_n$
and
$t_m$
, respectively, with the corresponding broadcast/receiving strengths, as quantified by the broadcast/receiving modes. Note that the locations of the FVs of high broadcast strength will provide insights into the location of actuation. Since broadcast modes also contain information about the directions along which the vortices are perturbed to result in the highest reduction of LEV circulation, they can also inform the design of momentum-based control inputs in nonlinear flows. Moreover, the perturbation time
$t_n$
corresponding to high levels of
$\Delta \varGamma (t_n, t_m) = \mathcal{J}(\boldsymbol{\varXi }^{{b}}(t_n, t_m))$
will be used to guide the timing of introducing control inputs for mitigating dynamic stall. The use of the network broadcast analysis to guide active flow control will be discussed in detail in the next section.
Figure 6. A demonstration of the results of the broadcast mode analysis. The analysis yields three components: (a) the total reduction in LEV circulation between
$t_n$
(perturbation time) and
$t_m$
(a later time); (b) the broadcast modes
$\boldsymbol{\varXi }^{{b}}$
, which identify the optimal set of vortical nodes where introduced perturbations result in the highest reduction of LEV circulation; and (c) receiving modes
$\boldsymbol{\varXi }^{{r}}$
, which highlight the vortical nodes most influenced by the perturbations originating from the broadcast nodes. Three pairs of
$[t_n, t_m]$
are chosen for the visualisation of broadcast and response modes.
4. Results
4.1. Baseline flow physics
To understand the dynamic stall flow physics of the two cases considered, the spanwise vorticity fields at a few representative instants are co-plotted along with the lift coefficients in figure 7. The instantaneous angle of attack and the plunge velocity corresponding to the selected instants are marked along the top horizontal axis for case A (pitching flat plate) and case B (plunging SD
$7003$
), respectively. The instantaneous separation point at each instant is also marked by a green dot on the suction side to indicate the flow evolution to leading-edge separation as a signature of dynamic stall (Carr, McAlister & McCroskey Reference Carr, McAlister and McCroskey1977; Leishman Reference Leishman1990).
Figure 7. The lift fluctuations for (a) case A and (b) case B are shown alongside the spanwise vorticity fields at representative time instants. For case A, the close-up view of the flow near the leading edge is presented in the dashed square for
$tu_{\infty }/L_{c} = 1.35$
. For case B, the vorticity fields are averaged in the spanwise direction. The separation point is marked with a green dot (
) in each vorticity field to identify leading-edge separation.
In the beginning of pitching and downstroke motions, respectively, for cases A and B, we observe that the separation point on the suction surface migrates upstream from the trailing edge, as indicated by instants (o) in figures 7(a) and 7(b). At a later instant (i), leading-edge separation occurs for both cases. At the same time, the stagnation point moves downstream along the pressure side of the aerofoil. Consequently, the flow directed toward the suction surface from the stagnation point negotiates the leading-edge curvature. This accelerates the flow near the leading edge, getting a strong suction force that contributes to a rise in lift force. Moreover, at instant (i) negative vorticity accumulates over the suction surface near the leading edge, which marks the onset of DSV formation (Karim & Acharya Reference Karim and Acharya1994). Starting from instant (ii) for case A, the aerofoil is held at
$\alpha = 45^{\circ }$
. The separation point remains near the leading edge, while the DSV convects away from the aerofoil surface and results in decreasing lift, as indicated by instants (ii) and (iii) in figure 7(a). For case B, the shear layer rolls up into smaller vortices due to the Kelvin–Helmholtz instability (Visbal Reference Visbal2011; Ramos et al. Reference Ramos, Wolf, Yeh and Taira2019). These vortices subsequently coalesce into the DSV, which is visible at instant (ii) in figure 7(b). As the DSV departs from the aerofoil surface, a reduction in the lift force from its peak value is observed. During the upstroke phase (
$t/T_{p} \in [0.5, 1]$
) for case B, the DSV passes beyond the trailing edge, accompanied by the shift of the separation point downstream along the suction surface. The flow gradually becomes fully attached, as seen at time instants (iii) and (iv) in figure 7(b).
From the preceding discussion, it is evident that the formation of the DSV core is preceded by instants (i) in figures 7(a) and 7(b). This instant captures critical flow features that contribute to the formation of the DSV core and eventually DSV in both cases. For case A, we characterise this instant by shear layer ’pinch-up’, where the separated shear layer from the leading edge immediately reattaches to the suction surface at
$x/L_{c} = 0.22$
. At this instant, the accumulation of negative vorticity near the leading edge has initiated, which precedes the onset of shear layer roll-up that forms the DSV core. For the turbulent flow case B, the mechanism of DSV core formation differs from that of case A. Here, the shear layer roll-up into spanwise vortices has initiated in the high-Reynolds-number flow, as observed from the instant (i) in figure 7(b). Moreover, we observe that the formation of the DSV core is driven by the merging of these spanwise vortices. This instant (i) for both cases marks the onset of DSV core formation, and we will revisit them when discussing the results of the network analysis in the next section.
Since the dynamics of the DSV play a crucial role in the pressure distribution on the aerofoil surface and lift fluctuation, interfering with its formation process by modulating the vorticity injection from the leading edge can hold a key to mitigating dynamic stall. Next, we perform the broadcast mode analysis on the vortical network to identify the effective actuation timing and location, with the goal of achieving effective modulation of vorticity injection to mitigate dynamic stall.
4.2. Broadcast mode analysis
As discussed in § 3.3, the network broadcast mode analysis seeks a broadcast mode,
$\boldsymbol{\varXi }^{{b}}(t_n, t_m)$
, as an optimal displacement perturbation introduced to vortical elements at
$t_n$
that maximises the reduction in LEV circulations over all later times until
$t_m$
. The objective function of the optimisation problem is defined as the sum of the reduction in LEV circulations introduced to the flow over
$t \in [t_n, t_m]$
due to a perturbation seeded at
$t_n$
. Therefore, the dependency of
$\mathcal{J}(\boldsymbol{\varXi }^{{b}}(t_n, t_m)) = \Delta \varGamma (t_n, t_m)$
on the time of perturbation,
$t_n$
, and the resulting broadcast mode,
$\boldsymbol{\varXi }^{{b}}(t_n, t_m)$
, provide valuable insights into the effective timing and location for flow actuation.
Figure 8. The broadcast mode analyses for case A (a) and case B (b), sweeping over different combinations of perturbation time
$t_n$
and a later time
$t_m$
. For both cases, the reduction in LEV circulation,
$\Delta \varGamma (t_{n}, t_{m})$
, is shown on the right-hand side, and the broadcast modes at three representative instances are visualised on the left-hand side. The broadcast modes are visualised by colouring the corresponding FV nodes by their broadcast strength.
Figure 9. The profiles of
$\Delta \varGamma (t_n, t_m)$
for a fixed
$t_m$
for (a) case A at
$t_{m}u_{\infty }/L_{c} = 3$
and (b) case B at
$t_{m}/T_{p} = 0.38$
. The peak value of
$\Delta \varGamma$
occurs at
$t_{n}u_{\infty }/L_{c} = 1.35$
for case A and at
$t_{n}/T_{p} = 0.21$
for case B, both marked by a red circle (). The corresponding broadcast and receiving modes are inserted for case A with
$[t_{n}, t_{m}]/(L_{c}/u_{\infty })= [1.35, 3.00]$
and, for case B, with
$[t_{n}, t_{m}]/T_{p} = [0.21, 0.38]$
.
This is shown in figure 8, where the broadcast mode analysis is conducted by sweeping over different combinations of perturbation time
$t_n$
and a later time
$t_m$
. We obtain the contour plots for
$\Delta \varGamma (t_n, t_m)$
below the main diagonal of figures 8(a) and 8(b), since
$t_m \gt t_n$
. Note that the introduction of LEVs in the DVM does not occur until
$tu_\infty /L_c = 1.05$
for case A and
$t/T_p = 0.11$
for case B. Consequently, the FV states at early times, as represented by instants (o) for both cases, are characterised by the presence of TEVs only. Hence, the broadcast modes at instants (o) are obtained only from these TEV nodes present in the flow. Moreover, the absence of LEVs at early times of perturbation (
$t_n u_\infty /L_c \lt 1.05$
for case A and
$t/T_p \lt 0.11$
for case B) leads to low levels of
$\Delta \varGamma$
for both cases. This indicates that introducing perturbations at time instances characterised solely by the presence of TEVs is ineffective in modulating the circulation injected from the leading edge. This can be attributed to the convection of TEVs away from the aerofoil surface that results in low levels of interactions with the aerofoil surface.
Starting from the instants where LEV shedding is initiated, a rapid increase in
$\Delta \varGamma$
is observed for both cases. The peaks of
$\Delta \varGamma$
occur at perturbation time
$t_{n}u_{\infty }/L_{c} = 1.35$
for case A and
$t_{n}/T_{p} = 0.21$
for case B, both after the initiation of LEV shedding. These peak instants are marked by instants (i) in figures 8(a) and 8(b). We make an important observation that these specific instants are corresponding to those leading to the DSV core formation, as previously discussed in figure 7. This suggests that the most effective timing of introducing perturbations for reducing LEV circulation corresponds to the early stages of DSV formation. Past instant (i), a decrease in
$\Delta \varGamma$
is observed along with the concluded formation of the DSV and its departure from the suction surface. This can be seen at the representative instants (ii) in figures 8(a) and 8(b) for both cases. At these later time instants, while the DSV continues to accumulate circulations of weaker LEVs, the DSV has become a robust structure and begins to convect away from the aerofoil surface. These observations underscore that once the DSV formation is concluded, the effectiveness of flow perturbation for modulating LEV circulation is significantly reduced.
Let us take a closer look at the
$\Delta \varGamma (t_n,\,t_m)$
profiles with respect to different perturbation times (
$t_n$
) and a fixed
$t_m$
in figure 9. Here, we choose
$t_{m}u_{\infty }/L_{c} = 3$
for case A and
$t_{m}/T_{p} = 0.38$
for case B and extract the
$\Delta \varGamma$
values along the constant-
$t_m$
lines from figure 8. Note that both of the chosen instants are associated with the concluded formation of the DSV over the suction surface. The resulting trends of
$\Delta \varGamma (t_n)$
are shown in figures 9(a) and 9(b) for case A and case B, respectively. In particular, the broadcast and receiving modes obtained from the peak
$t_n$
, highlighted by a red dot on the
$\Delta \varGamma (t_n)$
profile, are also inserted for both cases. Once again, we observe that the broadcast modes,
$\boldsymbol{\varXi }^{{b}}$
, highlight the vortical elements residing along the shear layer emanating from the leading edge. This suggests that perturbations introduced to the shear layer can potentially induce high levels of LEV circulation reduction, as quantified by
$\Delta \varGamma$
. Moreover, the receiving modes,
$\boldsymbol{\varXi }^{{r}}$
, obtained at the chosen later times
$t_{m}$
, reveal that the LEVs residing within the DSV core are those perturbed the most by the perturbations seeded into the flow at the earlier perturbation time
$t_{n}$
that results in the highest
$\Delta \varGamma$
. These observations indicate that the optimal timing for introducing perturbations for active flow control occurs slightly after the onset of leading-edge separation and prior to the formation of a DSV core, as identified by the LEV shedding in the DVM. With the identified timing for effective actuation, next we leverage the insights from the broadcast modes to identify an effective location on the aerofoil surface for actuator placement.
Note that the broadcast mode represents the optimal displacement perturbation applied to the Lagrangian vortical nodes. We can translate this Lagrangian perturbation to an Eulerian one by computing the induced velocity perturbation at an Eulerian grid point due to a displacement perturbation introduced at a Lagrangian vortical node by the Biot–Savart operator,
$\mathcal{V}_{{bs}}$
, used in the DVM. In this study we aim to place an actuator on the surface of the aerofoil. To this end, we compute the velocity perturbation on the aerofoil surface due to a broadcast mode perturbation via
\begin{align} \Delta \boldsymbol{u}(\boldsymbol{x}_{s}, k) = \sum ^{n_{\textit{f}v}}_{j = 1} \left [\mathcal{V}_{{bs}} \left (\,\boldsymbol{x}_{s},\,\boldsymbol{\xi }^{\textit{f}v}_{j, k} + \epsilon \boldsymbol{\varXi }_{j, k}^{{b}},\,\gamma ^{\textit{f}v}_{j}\right ) - \mathcal{V}_{{bs}}\left (\boldsymbol{x}_{s},\,\boldsymbol{\xi }^{\textit{f}v}_{j, k},\,\gamma ^{\textit{f}v}_{j}\right )\right ]\!, \end{align}
where
$\boldsymbol{x}_{s}$
parametrises a set of Eulerian points on the aerofoil surface at which the velocity perturbation is sought, and the time index
$k$
is associated with the time instants where
$\Delta \varGamma$
reaches its maximum. This will allow us to choose the actuator locations according to the magnitude of the velocity perturbation,
$\lvert \Delta \boldsymbol{u} \rvert$
. This is shown in figure 10, where we observe a region near the leading edge that exhibits high levels of
$\lvert \Delta \boldsymbol{u} \rvert$
for both cases. This occurs around
$x/L_{c} = 0.02$
from the pitching flat plate (case A) and
$x/L_{c} = 0.04$
for the plunging SD7003 (case B), indicating potentially effective locations to introduce control inputs. Moreover, when momentum-based actuators (Cattafesta III & Sheplak Reference Cattafesta III and Sheplak2011) are deployed for flow control, we can also inform the direction of actuation according to the orientation of the velocity perturbation at the peak location. These directions for both cases are also inserted in figures 10(a) and 10(b) for case A and case B, respectively. Note that the optimal direction shows a combination of both wall-normal and tangential components for both cases.
Figure 10. The induced velocity on the aerofoil surface due to the broadcast mode perturbation. The close-up view of the aerofoil section with high magnitude of perturbation velocity,
$\lvert \Delta \boldsymbol{u}\rvert$
, is presented along with its direction. (a) Case A: pitching flat plate and (b) Case B: plunging SD7003.
As a short summary, the broadcast mode analysis provided insights into the effective timing, location and direction for flow actuation. It suggests that applying actuation near the leading edge shortly before the formation of the DSV core can favourably reduce vorticity injection, potentially mitigating dynamic stall. We also note that, due to the use of small perturbation in the construction of the vortical network, the broadcast mode analysis conducted in this section is linear by nature. Although it offers guidance to the effective perturbations that result in a desirable initial departure of the flow states, the long-time nonlinear evolution of the flow under such control inputs needs to be further examined with great care. Therefore, to further interrogate these insights in their effectiveness of reducing the lift fluctuation and vorticity injection, we perform nonlinear flow simulations of the controlled flows, as discussed in the following section.
4.3. Broadcast mode analysis guided active flow control
In this section we design active flow control based on the results of broadcast mode analysis and examine its effectiveness in mitigating dynamic stall using LES of controlled flows. Here, we employ a flow actuator (Cattafesta III & Sheplak Reference Cattafesta III and Sheplak2011) on the aerofoil surface, according to the insights obtained from figure 10. This actuator is implemented as a time-varying Dirichlet boundary condition on the aerofoil surface, where the velocity profile is given by
\begin{align} \boldsymbol{u}_{\textit{act}} (\boldsymbol{x}_s, t) = u_{\textit{amp}} \times g(\boldsymbol{x}_s, \boldsymbol{x}_{\textit{act}}) \times \phi (t) \, \hat {\boldsymbol{e}}_{\textit{act}}. \end{align}
Here,
$u_{\textit{amp}}$
controls the actuation strength,
$g(\boldsymbol{x}_s, \boldsymbol{x}_{\textit{act}})$
provides a compact spatial support for the actuator placed at
$\boldsymbol{x}_{\textit{act}}$
on the aerofoil surface,
$\phi (t)$
controls the temporal duty cycle of the actuation input and
$\hat {\boldsymbol{e}}_{\textit{act}}$
denotes the direction of the momentum-based actuation. The specific forms and parameters for
$g(\boldsymbol{x}_s, \boldsymbol{x}_{\textit{act}})$
and
$\phi (t)$
are provided in table 1 for the actuator model:
\begin{align}\text{Actuator model:}\quad \boldsymbol{u}_\textit{act} (\boldsymbol{x}_s, t) {} {}& = u_\textit{amp} {} {}\times {} {}\underbrace{ {} {}\exp\left(-\frac{\left|\boldsymbol{x}_s - \boldsymbol{x}_\textit{act}\right|^2}{\sigma^2} \right)}_{g(\boldsymbol{x}_s)} {} {}\nonumber \\ & \quad \times {} {}\underbrace{\left[ {} {}\frac{e^{\tilde{\nu}(t - t_{\textit{on}})}}{1 + e^{\tilde{\nu}(t - t_{\textit{on}})}} - \frac{e^{\tilde{\nu}(t - t_{\textit{off}})}}{1 + e^{\tilde{\nu}(t - t_{\textit{off}})}} {} {}\right]}_{\phi(t)} \, \hat{\boldsymbol{e}}_{\textit{act}}\end{align}
Table 1. The velocity boundary condition for the actuator in the simulation set-up. The analytical forms, parameters and shape for the spatial and temporal profiles of actuation are provided in figures 11(a) and 11(b), respectively.
Figure 11. The spatial and temporal profiles of the actuator velocity boundary condition used in the simulation set-up are shown in (a) and (b), respectively.
We leverage the insights from figures 9 and 10 to inform the designs of actuator location,
$\boldsymbol{x}_{\textit{act}}$
, actuation direction,
$\hat {\boldsymbol{e}}_{\textit{act}}$
, and the temporal profile of actuation,
$\phi (t)$
. Here,
$\boldsymbol{x}_{\textit{act}}$
is determined according to the location where the induced velocity,
$\Delta \boldsymbol{u}$
, exhibits the maximum magnitude on the aerofoil surface, as indicated by figure 10. This results in the choices of
$\boldsymbol{x}_{\textit{act}}$
to be
$2\,\%$
- and
$4\,\%$
-chord downstream of the leading edge for case A and case B, respectively. The actuation direction,
$\hat {\boldsymbol{e}}_{\textit{act}}$
, is also determined in correspondence to the direction of
$\Delta \boldsymbol{u}$
at the same hot spot, as indicated by the arrows in figure 10. Lastly, the temporal profile is chosen such that the actuation strength reaches
$99\,\%$
of its maximum at the optimal timing according to the broadcast mode analysis in figure 9 (
$tu_\infty /L_c = 1.35$
for case A and
$t/T_p = 0.21$
for case B). We also quantify the actuation input using the momentum coefficient
$C_{\mu }$
, defined as
\begin{align} C_{\mu } = \frac {\iint \rho \lvert \boldsymbol{u}_{\textit{act}}(\boldsymbol{x}_s,t) \lvert \, \left [\boldsymbol{u}_{\textit{act}}(\boldsymbol{x}_s,t) \boldsymbol{\cdot }\hat {\boldsymbol{e}}_{n}(\boldsymbol{x}_{s})\right ] \textrm {d}\boldsymbol{x}_s \textrm {d}t}{\rho _{\infty } u^{2}_{\infty } L_{c} T^*}, \end{align}
where
$\hat {\boldsymbol{e}}_{n}(\boldsymbol{x}_{s})$
is the local surface normal vector of the aerofoil and
$T^*$
is the time scale of the wing motion. For case A (pitching flat plate),
$T^*$
is chosen to be the ramp-up time interval such that
$T^*u_\infty /L_c = 1$
. For case B (plunging SD7003),
$T^*$
is chosen to be the plunging period of the motion, i.e.
$T^* = 0.5T_p$
. The momentum coefficients for both cases are also reported in table 1, whose values are chosen in accordance with those reported in previous studies (Chen, Li & Hu Reference Chen, Li and Hu2014; Zigunov, Sellappan & Alvi Reference Zigunov, Sellappan and Alvi2022; Lin, Tsai & Tsai Reference Lin, Tsai and Tsai2023; de Souza et al. Reference de Souza, Wolf, Safari and Yeh2025).
The results of broadcast mode analysis guided flow control for both cases are shown in figures 12 and 13. For both cases, the lift fluctuations in time for baseline and controlled flows are provided with instantaneous flow visualisations at three representative instants. Note that the temporal profiles of control inputs are shown by the shaded regions for both cases. The instantaneous angle of attack and plunge velocity are also shown along the horizontal axis on the top. For both cases, we observe a significant decrease in lift fluctuation. For case A (pitching flat plate), a
$21\,\%$
reduction in peak lift is observed, as shown in figure 12(a). More excitingly, we observe a
$14\,\%$
decrease in peak lift for the 3-D turbulent flow over the plunging SD
$7003$
aerofoil, as shown in figure 12(b). Additionally, a reduction in moment fluctuations is observed in both cases, as shown in figure 13. For case A, a
$10.2\,\%$
reduction in the moment coefficient is observed around
$tu_{\infty }/L_{c} = 1.75$
. The local negative moment peak at
$tu_{\infty }/L_{c} = 1.75$
arises from the influence of the DSV, whereas the earlier peak around
$tu_{\infty }/L_{c} = 1$
is attributed to the added-mass effects. Furthermore, the
$13.2\,\%$
reduction in peak moment coefficient observed for case B (figure 13) is achieved. These reductions in lift and moment fluctuations for both cases suggest mitigation of dynamic stall using the broadcast mode analysis guided flow control.
Figure 12. Comparison of lift fluctuations between the baseline and controlled flows for (a) case A and (b) case B. The spanwise vorticity fields for case A and the
$Q$
isosurface coloured by spanwise vorticity for case B are provided at representative time instants for both baseline and controlled flows. For case A, a vortical structure in the controlled flow is highlighted by a green circle (
) to illustrate its evolution. For case B, the lift coefficient,
$C_{L}(t)$
, for both the baseline and control cases is cycle averaged over four periods of heaving-plunging motion.
In addition to the successful reduction of lift fluctuations achieved by the analysis guided flow control, the flow visualisations inserted in figure 12 also show that the physical mechanism that results in these reductions agrees with the ideation of the broadcast mode analysis and the insights of receiving modes. Flow modifications can be observed starting from instants (i) for both cases, which mark the optimal timings indicated by the broadcast mode analysis. At instants (i) for both cases, we observe that the flow physics near the leading edge are significantly disrupted by the control inputs. This is particularly clear in the visualisations for case A, where we observe that the leading-edge separation at instant (i) (
$tu_\infty u/L_c = 1.35$
) is effectively suppressed in the controlled flow. This critical modification interferes with the formation process of the DSV core. This is evident at instant (ii) (
$tu_\infty u/L_c = 1.75$
) in figure 12(a). At this instant, a DSV has formed near the leading edge in the baseline flow. However, in the controlled flow we observe that the vorticity accumulated near the leading edge at instant (i) continues to convect along the suction surface and forms a small vortex near the mid-chord at instant (ii). This is evident from figure 12(a), where the small vortex formed by accumulated vorticity in the controlled flow is marked by a circle. Moreover, this small vortex does not develop into the DSV observed at instant (iii). Instead, the DSV at instant (iii) is developed by accumulating the vorticity in between the leading edge and the small vortex at instant (ii). These observations indicate that the control input interferes with the DSV formation by dividing the accumulated vorticity into two vortical structures. Notably, instant (ii) is also the instant where the lift reaches its maximum value. For instant (iv) at
$tu_{\infty }/L_{c} = 2.2$
, the formations of the DSV are concluded for both baseline and controlled flows. At this instant, we observe a significantly smaller DSV in the controlled flow, suggesting that the reduction in lift fluctuation is indeed achieved by the modification of the formation process of the DSV.
Figure 13. Comparison of the
$z$
direction moment fluctuations, computed about the quarter chord, between baseline and controlled flows for (a) case A and (b) case B. For case B, the moment coefficient
$C_{M}(t)$
for both the baseline and controlled flow is cycle averaged over four periods of heaving-plunging motion.
Similar observations can be made for case B (plunging SD7003) presented in figure 12(b). For this high-Reynolds-number 3-D turbulent flow, the control input disrupts both the roll-up and vortex merging processes, as indicated at instant (i) (
$t/T_p = 0.21$
). In the baseline flow, the spanwise vortical structures developed from the shear layer roll-up coalesce and form the DSV core. In the controlled flow the shear layer rolls into smaller vortical structures over the suction surface. Moreover, we observe that the control inputs discourage the spanwise rollers from immediately merging and forming the DSV core. Instead, these spanwise vortices continue to convect downstream past the trailing edge. As a consequence, the formation of the DSV core is postponed due to the control input, similar to what we observed in case A. In particular, the DSV in the control flow observed at instant (ii) (
$t/T_p = 0.39$
) forms from the breakdown of the leading-edge flow into turbulent structures, bypassing the merging process between laminar spanwise rollers in the baseline flow. The DSV at instant (ii) also shows a smaller size in the controlled flow compared with that in the baseline. The strong DSV in the baseline flow triggers the formation of a coherent TEV of positive vorticity, which is not observed in the controlled flow until instant (iii). At instant (iii), we observe the departure of the DSV from the aerofoil surface in the baseline flow. Following the DSV departure, gradual flow reattachment is observed during the upstroke motion. For the controlled flow, the higher lift after
$t/T_p = 0.39$
than the baseline flow can be attributed to the delayed DSV departure due to its delayed formation.
Figure 14. The evolution of suction-surface pressure distribution,
$C_{p}(x) = (p(x) - p_{\infty })/0.5 \rho _{\infty } u^{2}_{\infty }$
, for the baseline flows and the difference in the distribution between the baseline and controlled flows for case A (a,b) and case B (c,d) is presented. Representative time instants indicated in figure 12 are marked along the horizontal axis in each panel. Note that the actuator over the suction surface is highlighted by a shaded magenta region for both cases in (b) and (d).
The modified DSV dynamics leaves a clear signature in the suction-surface pressure distribution, which leads to the change in lift fluctuation. This is shown in figure 14, where the temporal evolution of the suction-surface pressure distribution (
$C_{p}$
) for the baseline flows and the change in the suction-surface pressure distribution (
$\Delta C_{p}$
) for the controlled flows are shown for both cases. For case A, the baseline flow at
$tu_{\infty }/L_{c} = 1.35$
, marked as (i) in figure 14(a), exhibits a suction peak near the leading edge. This is the same time instant corresponding to the shear layer pinch-up instant, also marked as (i) in figure 12(a). This suction results from the accelerated flow around the leading edge directed from the stagnation point on the bottom surface of the aerofoil. The suction region progressively expands over the suction surface, particularly in the region
$x/L_{c} \lt 0.3$
during the interval
$tu_{\infty }/L_{c} \lt 2$
, coinciding with the DSV formation above the suction surface. This suction induced by the formation of the DSV is modified in the controlled flow, as indicated in figure 14(b). Downstream of the actuator region, a positive
$\Delta C_{p}$
indicates the pressure recovery, which leads to the reduction in lift. This recovery is related to both the delayed formation of the DSV and also its weaker strength.
We observe similar patterns for the suction-surface pressure distribution for case B, as shown in figure 14(c). Between
$tu_{\infty }/L_{c} = 0$
and
$0.21$
(instant i), the spanwise vortices from shear layer roll-up, as observed in figure 7(b), give rise to the low amplitude pressure fluctuations travelling downstream, as shown in figure 14(c). These rollers merge at instant (i) and form the core of the DSV, which builds up strong suction over the region of
$x/L_c \lt 0.3$
. After this instant, the DSV convects over the suction surface at a lower speed than the spanwise rollers observed in
$tu_{\infty }/L_{c} = [0, 0.21]$
, as indicated by the dark blue region in figure 14(c). Under the influence of control, we again observe pressure recovery starting from instant (i), due to the delayed formation of the DSV. Moreover, the reduced suction over the region
$x/L_{c} \gt 0.5$
during
$t/T_{p} \in [0.21, 0.39]$
reflects the diminished impact of the DSV in the controlled flow compared with baseline. Note here that for
$t/T_{p} \in [0.21, 0.39]$
, there is an enhanced suction for
$x/L_{c} \lt 0.5$
, which is due to the delayed presence of the DSV. However, the overall change in pressure distribution is favourable, resulting in the reduced lift fluctuation observed in figure 12.
Figure 15. The temporal evolution of negative vorticity within the bounding box for both the baseline and control flows in cases A and B. The bounding boxes are chosen such that the entire suction surfaces of the aerofoils are covered.
Recall that the objective function of the broadcast mode analysis is to reduce circulation injection during the shedding of LEVs. Here, we analyse the physical flows obtained from the LES to show that this objective is indeed achieved through the model-informed flow control. Following Visbal & Garmann (Reference Visbal and Garmann2018), we quantify the vorticity injection that forms the DSV by integrating the negative vorticity within a chosen bounding box. This is computed via
\begin{align} \varGamma (t) = \int _{\boldsymbol{x}\in {\it bounding\,box}} \omega _{z}^*(\boldsymbol{x}, t) \textrm{d}v(\boldsymbol{x}), \end{align}
where
$\omega _{z}^*$
represents negative spanwise vorticity that characterises the DSV. A higher negative value of
$\varGamma (t)$
indicates higher levels of negative vorticity introduced into the bounding box, which can be attributed to a stronger DSV.
The modifications in
$\varGamma (t)$
and the choices of the bounding boxes for both cases are shown in figure 15. We observe significantly decreased vorticity injection in the controlled flows for both cases. For the baseline flow in case A, the non-zero vorticity within the bounding box is observed prior to the onset of pitching at
$tu_{\infty }/L_{c} \lt 1$
. This is primarily due to the presence of the boundary layer along the aerofoil surface. Following the initiation of pitching at
$tu_{\infty }/L_{c} = 1$
, the amount of negative vorticity within the bounding box increases due to accumulation of vorticity within the bounding box to form the DSV core. This is similar to the observations in figure 12(a), where the shear layer continuously feeds negative vorticity into the DSV. With the model-informed control, a reduction in
$\varGamma (t)$
is observed at
$tu_{\infty }/L_{c} = 1.35$
for case A. This is marked as instant (i) in figure 12(a), which is the same instant in figure 7(a) that marks the initiation of the DSV core formation in the baseline flow. The reduction in
$\varGamma (t)$
becomes more pronounced at
$tu_{\infty }/L_{c} = 1.75$
, marked as instant (ii), where the lift reaches its peak in the baseline flow. At this instant, a
$26 \,\%$
reduction in vorticity within the bounding box is observed, suggesting that the reduction in lift fluctuation (shown in figure 14
b) is attributed to a weakened DSV in the controlled flow, agreeing with the objective of the broadcast mode analysis.
For case B presented by figure 15(b), we observe a similar reduction in the vorticity content within the bounding box but of a different control mechanism. For instant (i) at
$t/T_{p} = 0.21$
, the control input does not result in a notable decrease in
$\varGamma (t)$
. However, the small spanwise rollers in the controlled flow continue to travel along the suction surface without immediately merging near the leading edge to form the DSV core. This leads to a rapid decrease in
$\varGamma (t)$
past
$t/T_p = 0.21$
. The formation of the DSV is delayed in the controlled flow for instant (ii), where we observe a significant decrease in
$\varGamma (t)$
within the bounding box over
$t/T_{p} \in [0.21, 0.35]$
. The weaker DSV formed in the controlled flow suggests a lower injection of negative vorticity, which is also consistent with the objective of the broadcast mode analysis of the vortical network.
Figure 16. The time histories of the lift coefficient for the case of the pitching flat plate aerofoil (case A) are presented for five actuation cases: (a) actuation at the optimal time,
$tu_{\infty }/L_{c} = 1.35$
; (b) early actuation at
$tu_{\infty }/L_{c} = 1$
; (c) delayed actuation at
$tu_{\infty }/L_{c} = 1.5$
; (d) wall-normal suction; (e) wall-normal blowing. The corresponding time histories of negative circulation in the flow are shown in (f), (g), (h), (i) and (j), respectively, along with the spanwise vorticity field at
$tu_{\infty }/L_{c} = 1.75$
for each case.
Furthermore, we consider in figure 16 three different timings to turn on the control inputs to highlight the effectiveness of the model-guided flow control. Here, we focus on the case of 2-D flow over a pitching flat plate (case A). In addition to the optimal timing guided by the broadcast mode analysis (
$tu_{\infty }/L_{c} = 1.35$
, figures 16
a and 16
f), LES of controlled flows are also performed with an early actuation at the beginning of pitching (
$tu_{\infty }/L_{c} = 1$
, figures 16
b and 16
g) and a delayed actuation when the aerofoil is midway through its pitching motion (
$tu_{\infty }/L_{c} = 1.5$
, figures 16
c and 16
h). For these three cases, we examine the control effects by comparing their lift fluctuations (figure 16
a–c) and the circulations in the bounding box (figure 16
f–h) to those in the baseline flow. Flow visualisations using the vorticity field at the instant
$tu_{\infty }/L_{c} = 1.75$
, where the baseline flow exhibits peak lift, are also inserted in figure 16(f–h) for comparison.
Compared with the case of optimal timing for actuation, the cases with both early and delayed actuation exhibit suboptimal control effects with respect to the reduction of lift fluctuation. For the case with early actuation, the lift recovers to the same level as that in the baseline flow after the control input is turned off. This case exhibits the smallest reduction in circulation within the bounding box among the three cases at
$tu_{\infty }/L_{c} = 1.75$
and, as a result, fails to disrupt the formation of the DSV and consequently to reduce the peak lift. These observations are consistent with the insights from the broadcast mode analysis, which indicates that perturbing the vortices prior to the optimal time of
$tu_{\infty }/L_{c} = 1.35$
has minimal effect on reducing the LEV strengths and, consequently, the development of the DSV. On the other hand, for the case with delayed actuation, the initial formation of the DSV was not interfered. Instead, the delayed control input prevents the shear layer from continuously feeding vorticity into the DSV, leading to the earlier departure of the DSV compared with the other two cases. Comparing this case with the optimal one, although the same levels of lift reduction are observed at
$tu_{\infty }/L_{c} = 1.75$
, the lift increase in the early stage of DSV formation is not circumvented and results in suboptimal control effectiveness in terms of peak lift. Overall, we observe that the case with optimal timing not only achieves high reduction in peak lift, but the circulation in the bounding box at
$tu_{\infty }/L_{c} = 3$
as the DSV convects away is also the lowest among the three cases.
Let us consider two more actuation strategies in figure 16 to highlight the effectiveness of the model-guided flow control. In figures 16(d) and 16(e), respectively, wall-normal suction and wall-normal blowing are considered, since their uses have been widely examined for dynamic stall control in prior studies (McCloud, Hall & Brady Reference McCloud, Hall and Brady1960; Karim & Acharya Reference Karim and Acharya1994; Alrefai & Acharya Reference Alrefai and Acharya1996; Atik et al. Reference Atik, Kim, Van Dommelen and Walker2005; Gardner et al. Reference Gardner, Richter, Mai and Neuhaus2014; Müller-Vahl et al. Reference Müller-Vahl, Nayeri, Paschereit and Greenblatt2016). For these two cases, the flow actuation begins at the onset of pitching (
$tu_{\infty }/L_{c} = 1$
) and continues until
$tu_{\infty }/L_{c} = 3$
. Here, the momentum coefficient is fixed at
$C_{\mu } = 5.5\,\%$
across all cases. The cases of wall-normal suction and blowing exhibit suboptimal control effects relative to the network-guided strategy in reducing lift fluctuations (figure 16). For wall-normal suction (figures 16
d and 16
i), the reduction in peak lift coefficient is less pronounced than in the network-guided case. Although suction delays the formation of the DSV at
$tu_{\infty }/L_{c} = 1.75$
(figure 16
i), the injection of negative vorticity into the flow is slightly higher compared with the network-guided actuation. The case of wall-normal blowing produces the weakest modification to the peak lift coefficient (figure 16
e) and allows a DSV to form near the leading edge at
$tu_{\infty }/L_{c} = 1.75$
, as in figure 16(j). Collectively, these results show that the network-guided actuation achieves a more effective reduction in peak lift by suppressing circulation injection. This outcome further supports the use of broadcast mode analysis for designing flow control strategies for dynamic stall.
The broadcast mode analysis was designed to reduce the strength of LEVs. Its results motivate us to target the flow near the leading edge at the critical instants associated with the shear layer pinch-up. The receiving mode revealed that the circulation of LEVs residing in the DSV core at later times can be reduced. For both cases considered, a reduction in peak lift is observed in the LES of controlled flows, as shown in figure 12. This results from the reduction in the suction induced by the DSV, as in figure 14, attributed to the formation of a weaker vortex core from reduced injection of negative vorticity, as confirmed by figure 15. These observations are in strong agreement with predictions from the network-based analysis of discrete vortices.
5. Conclusions
We performed broadcast mode analysis of the vortical networks arising from dynamic stall flows to identify the optimal actuation location, timing and direction for suppressing the dynamic stall. This approach is demonstrated on two flows: a 2-D flow over a flat plate aerofoil at
$ \textit{Re} = 1000$
and a heaving-plunging spanwise-periodic 3-D turbulent flow over an SD
$7003$
aerofoil at
$ \textit{Re} = 60\,000$
. The framework leverages the time varying free-vortex states obtained from the DVM to understand the evolution of displacement perturbations to influence the circulation strength of LEVs within the vortical network. The broadcast mode analysis is reformulated as an optimisation problem that seeks the effective vortical nodes and the optimal timing in reducing the circulation strength of LEVs, with the objective of mitigating dynamic stall. The results of the analysis showed that the vortical nodes along the shear layer are the most effective nodes that appear at the instant of shear layer pinch-up. The analysis also showed that the broadcast mode based perturbation obtained from the optimal timing propagates to the vortical nodes residing in the DSV core at a later time. These insights from the broadcast mode analysis are used to inform the choice of actuator location, timing and direction in the simulation.
To validate the effectiveness of the network-guided actuation strategy, we implemented an actuator positioned near the leading edge in the LES and activated it at the optimal timing in the direction suggested by the network analysis. Simulation results show a significant reduction in peak lift by
$21 \,\%$
for the 2-D flows over the pitching flat plate aerofoil and
$14 \,\%$
for the 3-D turbulent flow over the plunging SD7003 aerofoil. The effect of flow control is found to modulate the negative vorticity injection into the DSV, consistent with the prediction of the broadcast mode analysis that the vortices within the DSV core are affected by vortical perturbations.
Lastly, we also note that the accuracy of the present broadcast mode analysis can be limited by that of the DVM employed in the present study (Ramesh et al. Reference Ramesh, Gopalarathnam, Granlund, Ol and Edwards2014; Lee et al. Reference Lee, Suresh Babu, Bryant and Gopalarathnam2022). We observed discrepancies in the lift coefficient and the departure of the DSV. These discrepancies could be attributed to the thin-aerofoil theory adopted by the DVM, wherein the aerofoil is represented by a distribution of BVs along the camber line. Such an assumption can be inaccurate for aerofoils with finite thickness, and the results can be potentially improved by representing the aerofoil surface with discrete vortex panels (Morino & Kuo Reference Morino and Kuo1974; Katz Reference Katz1981; Zhu et al. Reference Zhu, Wolfgang, Yue and Triantafyllou2002; Xia & Mohseni Reference Xia and Mohseni2017). The discrepancies could also be attributed to the setting that the critical leading-edge suction parameter was held constant in the present DVM. In such a case, potential improvements can be made by considering a time-varying critical leading-edge suction parameter, which has been shown to be more capable of capturing unsteady separation dynamics over aerofoils (Darakananda et al. Reference Darakananda, Eldredge, da Silva, Colonius and Williams2018; Narsipur et al. Reference Narsipur, Hosangadi, Gopalarathnam and Edwards2020, Reference Narsipur, Ramesh, Gopalarathnam and Edwards2023). Although improvement of the present DVM is possible, it reasonably captures the onset of DSV core formation and provides a simple yet practical framework for constructing a network-based model for the design of flow control. Its insights were also supported by the results of nonlinear simulations of controlled flows, demonstrating its capability to inform actuation strategies for dynamic stall mitigation.
In summary, the broadcast mode analysis guided flow control successfully reduced the lift fluctuation for both cases. This underscores the effectiveness of the network-based approach for determining the actuator location, timing and direction for dynamic stall mitigation by targeting the vorticity injection into the DSV. This study establishes a robust framework for active flow control design in the context of unsteady and transient flow separation, highlighting the potential of network-based analysis for flow control design.
Acknowledgements
We are grateful for the support given by the Army Research Office (W911NF-23-1-0109), monitored by program managers Drs Joseph Myers, Jack Edwards, and Kenneth Granlund. We acknowledge the computing resources provided by North Carolina State University High Performance Computing Services Core Facility. The 3-D turbulent flow simulations were performed using computational resources sponsored by the Department of Energy’s Office of Energy Efficiency and Renewable Energy and located at the National Renewable Energy Laboratory.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Solver and mesh set-up
We used the incompressible flow solver in the CharLES software package to simulate the flow in a body-fixed non-inertial reference frame. The solver solves the incompressible Navier–Stokes equations that write
\begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot } \boldsymbol{u} &= 0, \\[-12pt] \nonumber \end{align}
\begin{align} \frac {\partial \boldsymbol{u}}{\partial t} + (\boldsymbol{u} \boldsymbol{\cdot } \boldsymbol{\nabla }) \boldsymbol{u} &= -\boldsymbol{\nabla }p + 2 \boldsymbol{\nabla }\boldsymbol{\cdot } \left [ \left (\nu + \nu _{t}\right ) \unicode{x1D64E} \right ]+ \boldsymbol{f}, \\[9pt] \nonumber \end{align}
where the fictitious forcing term
$\boldsymbol{f}$
introduced for the non-inertial body-fixed reference (Agarwal & Deese Reference Agarwal and Deese1987; Beddhu, Taylor & Whitfield Reference Beddhu, Taylor and Whitfield1996; Kim & Choi Reference Kim and Choi2006) frame writes
\begin{align} \boldsymbol{f} = - 2 \dot {\alpha }(t) \times \boldsymbol{u} - \ddot {\alpha }(t) \times \boldsymbol{x} - \dot {\alpha }(t) \times \left ( \dot {\alpha }(t) \times \boldsymbol{x} \right ) - \ddot {{h}}(t). \end{align}
Here,
$\boldsymbol{u}$
denotes the fluid velocity relative to the non-inertial body-fixed reference frame,
$p$
is the pressure,
$\nu$
and
$\nu _t$
are respectively the reference viscosity and the subgrid-scale (SGS) eddy viscosity,
$\unicode{x1D64E} \equiv ( {1}/{2}) [\boldsymbol{\nabla }\boldsymbol{u} + (\boldsymbol{\nabla }\boldsymbol{u})^T ]$
is the strain rate tensor and
$\boldsymbol{x}$
is the position vector of the non-inertial body-fixed reference frame. The boundary condition imposed on the aerofoil is a no-slip condition. At the far field, the free-stream velocity is prescribed as
\begin{align} \boldsymbol{u} = -\dot {h}(t) - \dot {\alpha }(t) \times \boldsymbol{x}. \end{align}
Figure 17. Computational meshes for (a) the pitching flat plate and (b) the periodically heaving-plunging SD
$7003$
aerofoil is presented. For the flat plate aerofoil, the near-field mesh is shown together with the instantaneous spanwise vorticity field. For the SD
$7003$
aerofoil, the near-field mesh is shown with an isosurface of the instantaneous
$Q$
criterion, coloured by pressure.
For case A, 2-D simulations of flow around the flat plate aerofoil were performed by setting the SGS eddy viscosity to zero. An O-shaped computational grid, shown in figure 17(a), was employed with the plate pivoted about its leading edge, located at
$(x/L_{c}, y/L_{c}) = (0, 0)$
. The computational domain extends over
$x/L_{c}, y/L_{c} \in [-30, 30]$
in both the streamwise and transverse directions. The mesh was designed such that the first-cell thickness is
$\Delta y/L_{c} = 1 \times 10^{-3}$
, containing approximately
$5.5 \times 10^{5}$
grid nodes. Validation of the solver and mesh configuration was performed by comparing the lift coefficients obtained from CharLES with the numerical results of Ramesh et al. (Reference Ramesh, Gopalarathnam, Granlund, Ol and Edwards2014), as shown in figure 2(a).
For case B, 3-D LES were performed on the computational grid shown in figure 17(b). In this case, we used Vreman’s SGS model to calculate the SGS eddy viscosity
$\nu _t$
(Vreman Reference Vreman2004). The SD
$7003$
aerofoil is positioned such that its leading edge is located at
$(x/L_{c}, y/L_{c}) = (0, 0)$
. The computational domain extends over
$x/L_{c}, y/L_{c} \in [-20, 20]$
in the streamwise and transverse directions, and
$z/L_{c} \in [0, 0.4]$
in the spanwise direction. The mesh was designed with a first-cell thickness of
$\Delta x/L_{c} = 2.5 \times 10^{-3}$
,
$\Delta y/L_{c} = 2.5 \times 10^{-4}$
and
$\Delta z/L_{c} = 1 \times 10^{-2}$
, resulting in approximately
$3$
million mesh nodes. Validation of the solver and mesh setting for case B is performed by comparing the lift coefficient obtained from CharLES with the numerical results of Visbal (Reference Visbal2011), as shown in figure 2(b).
Open access
