Optimal waveform for fast synchronization of airfoil wakes

Abstract We obtain an optimal actuation waveform for fast synchronization of periodic airfoil wakes through the phase reduction approach. Using the phase reduction approach for periodic wake flows, the spatial sensitivity fields with respect to the phase of the vortex shedding are obtained. The phase sensitivity fields can uncover the synchronization properties in the presence of periodic actuation. This study seeks a periodic actuation waveform using phase-based analysis to minimize the time for synchronization to modify the wake shedding frequency of NACA0012 airfoil wakes. This fast synchronization waveform is obtained theoretically from the phase sensitivity function by casting an optimization problem. The obtained optimal actuation waveform becomes increasingly non-sinusoidal for higher angles of attack. Actuation based on the obtained waveform achieves rapid synchronization within as low as two vortex shedding cycles irrespective of the forcing frequency, whereas traditional sinusoidal actuation requires ${O}(10)$ shedding cycles. Further, we analyse the influence of actuation frequency on the vortex shedding and the aerodynamic coefficients using force-element analysis. The present analysis provides an efficient way to modify the vortex lock-on properties in a transient manner with applications to fluid–structure interactions and unsteady flow control.


Introduction
Unsteady periodic fluid flows are common in nature and engineering settings including vortex shedding over flapping wings, bluff bodies and airfoils.Modifying the vortex shedding behavior of wake flows is of high relevance to developing efficient engineering systems.However, controlling such flows is challenging owing to their periodically varying base states [Colonius and Williams, 2011].For the time-periodic base state, the timing of actuation becomes important.For this purpose, it is necessary to characterize the perturbation dynamics with respect to the time-periodic base state, which can be achieved using a phase reduction technique [Winfree, 1967, Kuramoto, 1984].The phase reduction approach expresses the perturbation dynamics using a single scalar phase variable.Recently, it has been used for studying periodic flows to reveal the phase sensitivity fields [Kawamura and Nakao, 2013, 2015, Khodkar and Taira, 2020, Kawamura et al., 2022, Loe et al., 2021, Iima, 2021], synchronization characteristics to external forcing [Taira and Nakao, 2018, Khodkar and Taira, 2020, Khodkar et al., 2021, Skene and Taira, 2022] and flow control [Nair et al., 2021, Loe et al., 2023] in a computationally inexpensive way.
Examining synchronization properties for periodic wakes can offer insights to modify the vortex shedding behavior and has several applications in unsteady flow control and fluid-structure interactions.Control of vortex shedding of wake flows has direct implications towards modifying the aerodynamic characteristics, reduction of structural vibration and noise emissions.Synchronization control has been studied in the context of vortex-induced vibrations for bluff-body wakes [Feng andWang, 2010, Konstantinidis andBouris, 2016].In addition, actuating a flow by taking advantage of synchronization can be efficient in enhancing the aerodynamic performance [Pastoor et al., 2008, Joe et al., 2011, Wang and Tang, 2018, Asztalos et al., 2021].Further, hydrodynamic synchronization is shown to result in efficient swimming in microscale swimmers at lower Reynolds number [Golestanian et al., 2011, Kawamura andTsubaki, 2018].Hence, it is beneficial to analyse the parameters that result in optimal synchronization in fluid flows.While most synchronization studies for fluid flows have characterized this asymptotic synchronization process to external sinusoidal actuation [Taira and Nakao, 2018, Khodkar and Taira, 2020, Herrmann et al., 2020, Giannenas et al., 2022], it is often desirable to modify the vortex shedding as quickly as possible for flow control to take effect.This study considers the fast synchronization of wakes to an external forcing signal for wake flows.
The concept of fast synchronization has been studied in biology to promote the rapid adjustment of the biological clock to jet lag and facilitate treatments for cardiac arrhythmias [Guevara andGlass, 1982, Granada andHerzel, 2009].In the context of biological and simpler oscillatory systems, Zlotnik et al. [2013] and Takata et al. [2021] applied the phase reduction approach to analytically obtain the optimal waveform for fast synchronization, maximizing the synchronization speed to external periodic forcing.In dynamical systems, synchronization is also referred to as entrainment [Strogatz, 1994].
In this study, we apply such a phase reduction approach to perform the phase-based fast synchronization for NACA0012 airfoil wakes at post-stall angles of attack with leading-and trailing-edge actuation.The overview of the fast synchronization analysis is shown in figure 1.We analytically find the optimal actuation waveforms, and the airfoil wake flows are actuated numerically using the optimal and sinusoidal waveforms at various forcing frequencies.The respective synchronization speeds are compared to validate the theoretical results.Further, we investigate the influence of actuation on the flow fields and the lift coefficients.The paper is organized as follows.The phase-based description and the framework to obtain the waveform for fast synchronization are presented in Section 2. The current approach is demonstrated with an example of NACA0012 airfoil wakes in Section 3. Conclusions are offered in Section 4.

Fast synchronization analysis through the phase-reduction approach
To obtain the optimal actuation waveform for fast synchronization for periodic fluid flows, we use phase reduction analysis [Taira andNakao, 2018, Kawamura et al., 2022].We identify the phase sensitivity fields that encode the effect of timing of actuation and then analytically solve an optimization problem to obtain the synchronization waveform in terms of the phase sensitivity function.

Phase reduction approach
We consider incompressible time-periodic fluid flows governed by the Navier-Stokes equations  = N ((, )), where  is the flow state.These equations are given by where Re is the Reynolds number and  is the pressure.For a time-periodic flow  0 (, ), it satisfies  0 (, +) =  0 (, ), where  is the time period of the limit cycle and Ω  = 2/ is the natural frequency of the system.Here, we define a phase  such that With the definition of , we can identify the full state vector of the limit cycle solution  0 (, ) at every .Given a stable limit cycle solution, with the frequency of the limit cycle being Ω  , the phase in the vicinity of the limit cycle can be described using the generalized phase variable Θ((, )).Thus, the generalized phase dynamics is described as Leveraging the phase dynamics, we can derive the phase response to sufficiently small perturbations, which provides the corresponding change to phase dynamics as Here, (, ) = ∇  Θ()| = 0 ( ,  ) is the spatial phase sensitivity field as it quantifies the phase response of the system to any given small perturbation and D is the considered spatial domain.This spatial phase sensitivity field can be obtained using either a direct impulse-based method [Taira and Nakao, 2018, Khodkar and Taira, 2020, Nair et al., 2021, Loe et al., 2021], or an adjoint-based approach [Kawamura and Nakao, 2013, 2015, Kawamura et al., 2022], or a Jacobian-free approach [Iima, 2021].We use the adjoint-based phase reduction framework to obtain the phase sensitivity fields in the present study as we can obtain the high-fidelity (, ) in all the flow variables by solving a single pair of forward and adjoint simulations.

Adjoint-based approach for phase sensitivity fields
We utilize the adjoint-based formulation to find the phase sensitivity fields  for time-periodic wakes.Let us consider the dynamics of a small perturbation  ′ (, , ) by linearizing the Navier-Stokes equations about the periodic base state  0 (, ).The perturbation dynamics is given by  ′ = L (, ) ′ , where L (, ) is the linearized Navier-Stokes operator.
To obtain the phase dynamics of the dominant limit cycle oscillation, we consider Floquet eigenfunction  and adjoint eigenfunction  * corresponding to the zero eigenvalue.This Floquet-zero eigenvalue corresponds to the phase degree of freedom for the stable limit cycle dynamics.Thus phase dynamics is obtained by projecting the perturbed dynamics in equation 4 on the adjoint eigenfunction as We note that the norm of  * is arbitrary and the normalization condition of  * is appropriately chosen to satisfy More details on the properties of  * are given in our earlier work [Kawamura et al., 2022].Here, by comparing equations 6 & 5, for perturbations in the form of velocity, we obtain (, ) =  * (, ).Hence, phase sensitivity is the adjoint-zero eigenfunction of the linearized Navier-Stokes operator.The spatial phase sensitivity fields can be obtained by solving the dynamics, which in two dimensions is governed by the linearized adjoint equations of where  * = ( * ,  * ) and  = (, ).Thus, the phase sensitivity fields with respect to perturbations in the velocity field are obtained by seeking a periodic solution for equation 1 and solving the system of adjoint equations (equation 7).Since, the adjoint equations are analogous to the Navier-Stokes equations, the same numerical scheme can be used to solve them.An overview of the phase description for airfoil wakes is shown in figure 2. The phase is defined based on the lift coefficient   −   plane, where  = 0,  correspond to mean   ,  = /2 corresponds to maximum   and  = 3/2 corresponds to minimum   [Taira and Nakao, 2018].

Synchronization analysis to external periodic forcing
Establishing this oscillator dynamics enables us to study the synchronization characteristics of the system to an external periodic forcing signal with a frequency Ω  different from the wake shedding frequency Ω  .We introduce a localized periodic forcing at location  0 at a forcing frequency Ω  .The spatial profile of forcing is given by a Dirac delta function ℎ() = ( −  0 ).Using equation 5, the governing phase dynamics becomes where (, ) =  (Ω  )ℎ() and the local phase sensitivity function is given by  () = ∫ D (, )ℎ() = ( 0 , ).To characterize the synchronization of the system to external forcing, we consider the relative phase () between the phase of the system  () and that of the forcing signal Ω   as () =  () − Ω  .The dynamics of the relative phase is provided as where ΔΩ = Ω  − Ω  .The asymptotic behaviour of relative phase dynamics can be obtained by averaging over a period of forcing [Kuramoto, 1984, Ermentrout andKopell, 1991], where is the phase coupling function and ΔΩ = Ω  − Ω  .Synchronization occurs if the relative phase becomes a constant, i.e.,  → 0. Hence, the synchronization condition is given as The synchronization condition determines the forcing frequency required to synchronize the dynamics to the external actuation based on the phase coupling function.
We aim to identify the optimal periodic actuation to synchronize the system to a forcing frequency as quickly as possible.
Hence, the rate of convergence of  to a fixed point  * should be maximized to satisfy Therefore, we can formulate an optimization problem to maximize | |, which occurs when − Γ ′ ( * ) is large.Here −Γ ′ ( * ) is the synchronization speed .The cost function J is therefore formulated as where  and  are Lagrangian multipliers and The first term corresponds to maximizing the synchronization speed, the second term constrains the energy of actuation, and the third term directly follows from equation 13.
Since the synchronization is independent of the initial phase, without loss of generality, we consider the fixed point,  * = 0.This optimization can be solved analytically using the calculus of variations [Zlotnik et al., 2013].The optimal waveform for fast synchronization can then be derived as Hence, once we compute the local phase sensitivity function  (), the optimal waveform for fast synchronization can be analytically found using equation 15 for various Ω  and .The optimal speed of synchronization is characterized by then computing −Γ ′ (0) using the optimal waveform given by equation 15.Even though the current formulation to obtain the optimal waveform is defined for a single pointwise actuation, this directly extends to the case with multiple pointwise actuators or spatially distributed actuators following from equation 5.This is reflected as a change in the inner product ⟨•⟩ in equations 14 and 15, where the inner product would be computed as an integration for the multiple local phase sensitivities and actuation waveforms.Next, we uncover these optimal waveforms for the airfoil wakes using the local phase sensitivity functions and assess their performance for fast synchronization.
3 Phase synchronization analysis of airfoil wakes

Computational set-up
This study considers the two-dimensional incompressible laminar flow over NACA0012 airfoils at angles of attack,  = 35 • , 45 • and 55 • and chord-based Reynolds number of Re =  ∞ / = 100, where  ∞ ,  and  are the free-stream velocity, airfoil chord length and kinematic viscosity, respectively.The flow dynamics is governed by incompressible Navier-Stokes equations 1 and the obtained flow fields present with periodic vortex shedding [Kawamura et al., 2022].The actuation in equation 4 is introduced as a localized force with the form (, Ω  ) =  (Ω  )( −  0 ), where  0 is a forcing location.The Dirac delta function is approximated with a three-cell discrete delta function [Roma et al., 1999].
The periodic flows over the airfoil are computed numerically through the immersed boundary projection method [Taira andColonius, 2007, Kajishima andTaira, 2016].For the numerical simulation, we consider a computational domain 30,30].The quarter-chord of the airfoil is placed at the origin.The smallest grid size is set to Δ min / = 0.02, and the time step is chosen to be Δ = 0.005.The present computational setup has been validated and is the same as that used in Kawamura et al. [2022].The same computational setup is used for adjoint simulations of the phase sensitivity fields.

Synchronization analysis for airfoil wakes
The spatial phase sensitivity fields with respect to the streamwise and transverse velocity components   and   for NACA0012 airfoils at  = 35 • , 45 • and 55 • are obtained through the adjoint-based approach described in Section 2.2.Using the obtained spatial phase sensitivity fields, we can compute the optimal waveform for fast synchronization at each grid point as per equation 15, which is then used to obtain the optimal synchronization speed at each grid point.We investigate the effect of the angle of attack on the synchronization speed and waveforms of NACA0012 airfoil wakes, as shown in figure 3. We consider the case when the forcing frequency Ω  = Ω  and ΔΩ = 0.It follows from equation 15 that the optimal actuation waveform at each point is proportional to the corresponding derivative of the local phase sensitivity function  ′ ().The spatial distributions of synchronization speed  around the airfoil found using the optimal waveform for  = 35 • , 45 • and 55 • are depicted in figures 3()-().As the angle of attack increases, the overall magnitude of the synchronization speed decreases, indicating an increased difficulty in synchronization for higher post-stall angles of attack.With an increase in , we observe stronger and larger leading-and trailing-edge vortex structures.To achieve synchronization with external forcing, the vortex formation time and the length scale have to be modified.This therefore becomes challenging with higher , which is reflected in reduced synchronization speed.Further, we also note that the white region around the airfoil corresponds to a small optimal synchronization speed, indicating that, irrespective of the actuation energy, these spatial locations are not conducive for flow modification.This means that the actuation effort must penetrate the outside of the boundary layer.
We also observed that, for all , the local maxima in the synchronization speed are attained near the leading and trailing edges, suggesting them as optimal actuation locations for synchronization (indicated as black dots).The leading and trailing edges are the most sensitive regions since they are specific regions in the flow field with high curvature.Further, the regions with high synchronization speed become more compact and are concentrated at the leading and trailing edges with an increase in  due to the earlier flow separation and the concentration of gradients at the leading and trailing edges at higher .Even though we considered three angles of attack, we expect this trend in synchronization speed and optimal waveform to hold true for much higher angles of attack.For the NACA0012 airfoil at  = 100, we observe a steady wake until  ≈ 20 • .This results in a constant lift coefficient and, hence, an ill-defined phase.Further, as  → 90 • , we approach a zero mean lift coefficient.However, for 30 • ≤  ≤ 90 • , we observe periodic vortex shedding and we can leverage the optimal waveform analysis.We expect a similar trend in the optimal waveform and in the synchronization speed with an increase in the angle of attack.An increase in the angle of attack results in the formation of stronger leading-and trailing-edge vortices, thereby increasing the difficulty in synchronization and the reduction in synchronization speed.The asymmetry in the vortex formation and roll-up between the leading-and trailing-edge vortices also increases with most , resulting in a non-sinusoidal optimal waveform.However, for  → 90 • , we approach a symmetric bluff-body vortex shedding.Hence, overall, the optimal waveform outperforms the synchronization speed of a sinusoidal waveform at most higher angles of attack.It is noteworthy that, due to the difficulty in synchronization at higher angles of attack, we will require a larger actuation effort to synchronize the wake to a different frequency.
The optimal actuation waveforms in the  and  velocity directions, at the leading and trailing edges for various , are shown in figures 3()-().As  increases, the optimal waveform becomes increasingly non-sinusoidal, due to the asymmetry in the vortex formation and shedding process near the leading and trailing edges at higher angles of attack.Further, the optimal waveform at the trailing edge at higher angles of attack suggests a smaller time duration where actuation is significant (for 0 <  < 3/4 in figure 3(e)), in comparison with the leading-edge optimal waveform.This is in line with flow physics, as we observe a more compact and stronger vortex roll-up at the trailing edge when compared to the vortex formation at the leading edge.We would like to point out that these optimal waveforms are obtained by independently maximizing the synchronization speed using the respective local phase sensitivity functions.We can also obtain the optimal waveforms at the leading and trailing edges by optimizing the synchronization speed using the local phase sensitivity functions simultaneously.For this present study, both these cases lead to similar results with minimal modification, where simultaneous optimization results in the same waveforms but with more actuation energy at the trailing edge than at the leading edge.This difference should be carefully considered for more complex flow fields when using multiple actuation locations.Next, we numerically validate the synchronization analysis by introducing actuation at the optimal actuation locations near the leading and trailing edges (as shown in figure 3)(d)-(e)).Here, we consider the optimal waveform and a sinusoidal waveform with the same averaged actuation direction at different forcing frequencies.We present the numerical results at  = 45 • as a representative case.The numerical results of synchronization for a forcing frequency within 5% of the natural frequency are shown in figure 4. Here, we choose an actuation amplitude of  = 0.1 to achieve synchronization for laminar flows at higher angles of attack.
To assess the synchronization speed, we consider cycle-to-cycle variations of the   coefficient and measure the inter-peak phase difference Δ  for each cycle as shown in figure 4().The optimal waveform actuation achieves synchronization in two shedding cycles, in comparison to O (10) shedding cycles for the sinusoidal waveform (see figure 4()) for different forcing frequencies.Since the optimal waveform is based on the phase sensitivity function, it can efficiently identify the "when" and "how" to efficiently synchronize the system to an external forcing signal, thus achieving fast synchronization.The effect of actuation frequency on flow physics is examined using the instantaneous vorticity fields of synchronized and unperturbed in figure 4().We observe streamwise elongation of the leading and trailing edge vortices for lower frequency actuation, Ω  = 0.95Ω  when compared with the unperturbed case.
On the other hand, we observe more compact leading and trailing edge vortices for higher actuation frequencies, Ω  = 1.05Ω  .Hence, the modification of vortex shedding frequency through optimal waveform actuation is achieved by modifying the vortex formation length scale near the leading and trailing edges.Thus, the phase-sensitivity-based optimal waveform deviates from the sinusoidal waveform to target more actuation energy at the right time to achieve rapid flow modification.
To further examine the effect of the present actuation over the lift coefficients, let us monitor the force elements [Chang, 1992].Force element theory enables us to identify the flow structures responsible for lift generation.We compute an auxiliary potential function   that satisfies the Laplace equation ∇ 2   = 0, with the boundary condition − • ∇  =  •   on the airfoil surface, where   is the unit vector in the lift direction.The lift force is obtained by taking the inner product of ∇  with the momentum equation and integrating with D in two dimensions as where the first term denotes the surface integral and the second term denotes the line integral on the airfoil surface.The first integrand herein referred to the lift element   , and is used to monitor the effect of vortical structures on the lift force.
The lift coefficient   for a vortex shedding period for the unperturbed and the actuation frequencies Ω  = 0.95Ω  and 1.05Ω  are shown in figure 5.The snapshots are shown corresponding to the unperturbed flow fields (black box), synchronized flow fields at Ω  = 0.95Ω  (blue box) and Ω  = 1.05Ω  (red box) at max   (Δ) and min   (∇).Owing to the actuation, we notice a significant change in   compared to the unperturbed case for both frequencies, especially for Ω  = 1.05Ω  .For a 5% increase in frequency (Ω  = 1.05Ω  ), we observe a 17% increase in max   and a 8% increase in mean   compared with the unperturbed case.However, we do note that a similar amount of actuation is introduced to the flowfield.It is noteworthy that the swift modification of the shedding timing is achieved by the lift increases for high-frequency actuation.We further analyse the wake with the lift elements   (Δ) for unperturbed (black box) and high-frequency actuation (red box) as shown in figure 5. We observe a strong compact positive   near the leading and trailing edges.This suggests that the increased strength and compactness of the vortex increases the local circulation, and thereby the lift force [Eldredge and Jones, 2019].
We now consider the low-frequency actuation (Ω  = 0.95Ω  ), where we do not observe a significant change in mean   in comparison with the unperturbed case.In contrast to the high-frequency actuation, the optimal waveform actuation achieves a reduction in the wake shedding frequency through a reduction in min   .The lift force elements   corresponding to this case (blue box,∇) show a streamwise-elongated positive force element effectively pushing away the shear layer from the airfoil surface, thereby reducing the overall lift force.Overall, through the lift element theory, we identified that high-frequency actuation using the optimal waveform results in compact vortices at the leading and trailing edges, and the lower wake shedding frequency is achieved by streamwise elongation of the vortices at the leading and trailing edges.Through a high-frequency actuation using the optimal waveform, a transient increase in lift is observed, albeit with a considerable actuation effort.By demonstrating the effectiveness of optimal waveform analysis for  = 45 • , we show the potential of this method for the analysis of a wide range of periodic fluid flows and their control in a transient manner.

Conclusions
We presented a theoretical framework to find an optimal actuation waveform for maximizing the synchronization speed for periodic fluid flows.This was demonstrated for periodic post-stall airfoil wakes using localized forcing.We leveraged the phase reduction approach to identify the sensitivity with respect to the vortex shedding phases, thereby identifying the right time and direction of actuation for efficient synchronization.The optimal actuation waveform for fast synchronization departs from a sinusoidal waveform for higher angles of attack.We showed that the optimal waveform significantly outperforms the sinusoidal waveform in terms of synchronization speed.We further identified that the modification of wake shedding frequency is achieved by the elongation of vortical structures, whereas synchronization to a higher frequency is achieved by compacting vortical structures near the leading and trailing edges.The present study based on phase reduction with an optimal waveform approach holds potential to develop transient flow control strategies that produce a quick response.

Figure 1 :
Figure 1: Fast synchronization analysis of flow over NACA0012 airfoil at  = 55 • and Re = 100.() Periodic actuation using the optimal waveform.() Comparison of synchronized flowfield for a forcing frequency Ω  = 0.9Ω  with the baseline vorticity field .() Lift coefficient   () when actuated with the fast synchronization and sinusoidal waveforms.

Figure 2 :
Figure 2: An overview of the phase reduction approach for flow over a NACA0012 airfoil at  = 45 • and Re = 100.() Definition of phase based on the lift coefficient,   −   plane.() Instantaneous spatial phase sensitivity fields   and   with respect to the perturbations in velocity fields.

Figure 5 :
Figure 5: Characterization of   for the NACA0012 airfoil at  = 45 • and Re = 100.Variation of   when actuated with forcing frequencies Ω  = 1.05Ω  and 0.95Ω  .Instantaneous lift force elements   are shown for the actuated and unperturbed cases at max   (Δ) and min   (∇).