Networked oscillator based modeling and control of unsteady wakes

TL;DR

This paper introduces a networked oscillator model for analyzing and controlling unsteady wake flows, capturing energy exchanges among flow modes to improve flow control and reduce drag.

Contribution

It presents a novel networked oscillator framework that models modal interactions in unsteady flows, enabling effective control strategies for wake suppression.

Findings

Networked oscillator model outperforms empirical Galerkin models in capturing modal dynamics.

The approach successfully suppresses wake unsteadiness and reduces drag in flow over a circular cylinder.

Linear regression effectively identifies energy transfer networks among flow modes.

Abstract

A networked oscillator based analysis is performed for periodic bluff body flows to examine and control the transfer of kinetic energy. Spatial modes extracted from the flow field with corresponding amplitudes form a set of oscillators describing unsteady fluctuations. These oscillators are connected through a network that captures the energy exchanges amongst them. To extract the network of interactions among oscillators, amplitude and phase perturbations are impulsively introduced to the oscillators and the ensuing dynamics are analyzed. Using linear regression techniques, a networked oscillator model is constructed that reveals energy transfers and phase interactions among the modes. The model captures the nonlinear interactions amongst the modal oscillators through a linear approximation. A large collection of system responses are aggregated into a network model that captures interactions for general perturbations. The networked oscillator model describes the modal perturbation dynamics better than the empirical Galerkin reduced-order models. A model-based feedback controller is then designed to suppress modal amplitudes and the resulting wake unsteadiness leading to drag reduction. The strength of the proposed approach is demonstrated for a canonical example of two- dimensional unsteady flow over a circular cylinder. The present formulation enables the characterization of modal interactions to control fundamental energy transfers in unsteady vortical flows.

Equations

$$\mathbf{u}(\mathbf{x},t)\approx \bar{\mathbf{u}}(\mathbf{x})+\sum _{j=1}^N{a}_j(t){\mathbf{\varphi }}_j^{\mathbf{u}}(\mathbf{x}),$$

$$z_m=a_{2j-1}+i{a}_{2j}=r_m\exp (i{\theta }_m),$$

$$\mathbf{u}(\mathbf{x},t)=\bar{\mathbf{u}}(\mathbf{x})+\sum _{m=1}^M\left[\Re (z_m(t)){\mathbf{\varphi }}_{2m-1}^{\mathbf{u}}(\mathbf{x})+\Im (z_m(t)){\mathbf{\varphi }}_{2m}^{\mathbf{u}}(\mathbf{x})\right],$$

$${\dot{z}}_m^b=z_m^b({\lambda }_m-|z_m^b{|}^2+i{\Omega }_m^b)$$

$$z_m^{\prime }=a_{2j-1}^{\prime }+i{a}_{2j}^{\prime }={ϵ}_m^{\prime }{r}_m^b\exp (i[{\theta }_m^b+{\theta }_m^{\prime }]),$$

$$z_m=r_m\exp (i{\theta }_m)=z_m^b+z_m^{\prime }.$$

$${\zeta }_m=z_m^{\prime }/z_m^b={ϵ}_m^{\prime }\exp (i{\theta }_m^{\prime }).$$

$$\begin{array}{r}{E}_m^{\prime }=\frac{1}{2}\left(|z_m{|}^2-|z_m^b{|}^2\right)=\frac{1}{2}[(1+{ϵ}_m^{\prime }{)}^2-1](r_m^b{)}^2,\\ E_m^{\prime }(t_0)={\beta }_m{E}^b,\end{array}$$

$${\dot{\zeta }}_m=\sum _{n=\text{I}}^M[{\mathbf{A}}_{\mathcal{G}}{]}_{mn}\left({\zeta }_n-{\zeta }_m\right)=-\sum _{n=\text{I}}^M[{\mathbf{L}}_{\mathcal{G}}{]}_{mn}{\zeta }_n,$$

$$[{\mathbf{A}}_{\mathcal{G}}{]}_{mn}=\{\begin{array}{ll}{w}_{mn}=|w_{mn}|\exp (i\angle w_{mn})& \text{if there is an edge from }n\text{ to }m\\ 0& \text{otherwise.}\end{array}$$

$$\begin{array}{r}\dot{a_j}=\sum _{k=0}^N{\psi }_{jk}{a}_k+\sum _{k,l=0}^N{\chi }_{jkl}{a}_k{a}_l,j=1,2,\dots N,\end{array}$$

$$C_D=\frac{F_D}{\frac{1}{2}\rho U^{2}d}\text{and}{C}_L=\frac{F_L}{\frac{1}{2}\rho U^{2}d},$$

$$\begin{array}{r}{△}_m=\frac{1}{\text{max}(|{\dot{\zeta }}_m|)}\sqrt{\frac{{\sum }_{k=1}^{n_t}|{\hat{\dot{\zeta }}}_m(k)-{\dot{\zeta }}_m(k){|}^2}{n_t}},\end{array}$$

$$\dot{\mathbf{\zeta }}=-{\mathbf{L}}_{\mathcal{G}}\mathbf{\zeta }+\mathbf{B}\mathbf{v},$$

$$\dot{\mathbf{\zeta }}=(-{\mathbf{L}}_{\mathcal{G}}-\mathbf{B}\mathbf{K})\mathbf{\zeta },$$

$$J={\int }_0^{\infty }\left[\mathbf{\zeta }(t{)}^T\mathbf{Q}\mathbf{\zeta }(t)+\mathbf{v}(t{)}^T\mathbf{S}\mathbf{v}(t)\right]\mathrm{d}t,$$

$${\dot{z}}_m=-\sum _{n=\text{I}}^M[{\mathbf{L}}_{\mathcal{G}}-\mathbf{B}\mathbf{K}{]}_{mn}{z}_n.$$

$$P_D=\frac{F_{D}U}{\frac{1}{2}\rho U^{3}d}\text{and}{P}_a=\frac{⟨\mathbf{f}\cdot \mathbf{u}⟩}{\frac{1}{2}\rho U^{3}d},$$

$$\mathbf{u}(\mathbf{x},t)={\mathbf{u}}^{\ast }(\mathbf{x})+\sum _{j=1}^N{a}_j(t){\mathbf{\varphi }}_j^{\mathbf{u}}(\mathbf{x}).$$

Full text

††volume: ???

1 Introduction

Oscillations play an important role in many physical and biological systems. These oscillations often result from a set of self-sustaining (autonomous) entities called oscillators. Biological oscillators, including neurons and heart cells are integral to the various rhythms and regulatory systems of the human body. Such collective rhythms arise from the coupling of multiple oscillators with the physics encapsulated by the transfer of energy between them (Strogatz, 2004). There has been a rich history of studies on the collective dynamics of oscillators, in particular by Kuramoto (1984) and Strogatz (2014). The foundational work laid out by Kuramoto (1975) elegantly describes the interactive phase dynamics between oscillators. Mutual synchronization of a system occurs when interacting oscillators affect their phases so as to spontaneously lock on to a particular frequency or phase (Pikovsky et al., 2003). In the works of Aizawa (1976) and Mirollo & Strogatz (1990), the oscillator phase interactions are generalized to incorporate amplitude variation effects. Amplitude and phase coupling between oscillators result in a variety of interesting physical phenomena.

Unsteady fluid flows are often characterized by temporal oscillations. In flows past bluff bodies, oscillatory behavior of the flows is exhibited through shedding of coherent vortices observed in the wake. Such periodic shedding generates unsteady forces on the body which can lead to detrimental increase in drag associated structural fatigue due to the emergence of flow-induced vibrations (Williamson & Govardhan, 2004; Sarpkaya, 2004). In the work of Roshko (1954), the relationship between form drag and vortex shedding was explored in detail. It was demonstrated that unsteady force oscillations and drag can be reduced by mitigating the wake unsteadiness. Since then, a myriad of studies using active and passive flow control strategies have focussed on controlling bluff-body wake vortex shedding and the resulting unsteady forces, summarized in a review by Choi et al. (2008). Although there have been tremendous breakthroughs in applying flow control for drag reduction, only a few studies make use of the fundamental energy transfer mechanisms and interactions in unsteady fluid flows and controlling the flow unsteadiness therein.

The wake unsteadiness stems from the inherent vortical oscillations. Using modal decomposition techniques, oscillations embedded in fluid flows can be extracted naturally as spatial structures (modes) and their associated temporal weights (Taira et al., 2017). For time-periodic flows, individual coherent structures are described by conjugate mode pairs extracted from these modal decomposition techniques. These mode pairs can be viewed as a set of modal oscillators exhibiting periodic fluctuations. In particular, proper orthogonal decomposition (POD) (Sirovich, 1987; Aubry et al., 1988; Berkooz et al., 1993) and dynamic mode decomposition (DMD) (Schmid, 2010; Rowley et al., 2009; Kutz et al., 2016) techniques can decompose snapshot flow field data into modal oscillators based on energy and dynamics of the flow, respectively. The interactions between these modal oscillators govern the behavior of unsteady fluid flows. The general behavior of nonlinear flows can also be described by spectral analysis of the linear, infinite-dimensional Koopman operator which yields Koopman modes (Mezić, 2005; Rowley et al., 2009; Mezić, 2013).

The interactions between the modes in unsteady fluid flows can be encapsulated by reduced-order models which may dramatically reduce the computational complexity to analyze fluid flows of interest (Aubry et al., 1988; Noack et al., 2003). The projection of the Navier–Stokes equations onto the modal basis results in an empirical Galerkin formulation (Holmes et al., 2012). The empirical Galerkin models contain linear and quadratic nonlinear terms which govern the dynamics of the modes. However, strong nonlinear interactions between modes in unsteady fluid flow pose a significant challenge in the accuracy of reduced-order models (Lassila et al., 2014; Brunton & Noack, 2015). The study by Sapsis (2013) uses dynamically orthogonal equations which relates the low-dimensionality of fluid flow attractors with energy exchanges between the mean and dynamical components of the flow. In the current study, we examine how modal oscillators in unsteady fluid flows interact to distribute energy in non-equilibrium conditions. This involves a perturbation-response based approach to characterize the interactions and extract both phase and amplitude coupling between modal oscillators in a fluid flow. In this approach, perturbations are introduced to the modal oscillators and the ensuing energy transfer dynamics are analyzed using a graph-theoretic framework.

A graph $\mathcal{\boldsymbol{G}}=\{\mathcal{V},\mathcal{E},\mathcal{W}\}$ consists of a set of nodes $\mathcal{V}$ connected by edges $\mathcal{E}$ with associated edge weights $\mathcal{W}$ (Newman, 2010; Nair & Taira, 2015). The network nodes form the quantities of interest with the interactions between them as edges. Such a simplistic viewpoint of analyzing collections of interactions in network science has had far-reaching socio-economic as well as scientific impacts (Barabási, 2016). Network analysis is primarily concerned with extracting the interactions between quantities of interest (Bollobás, 1998; Newman, 2010) and has found widespread applications in social sciences (Otte & Rousseau, 2002), biological sciences (Zhu et al., 2007) and many other fields (Dorogovtsev & Mendes, 2013). The studies of various technological networks like the Internet and World Wide Web have improved the efficiency of global communications tremendously (Newman, 2010). The Human Connectome project (Wedeen et al., 2005; Hagmann et al., 2007; Sporns, 2011) provides access to big neural data to explore connectivity pathways in the brain and has led to a wealth of new discoveries. In epidemiology, network analysis has aided in forecasting epidemics and designing appropriate containment and control measures (Salathé & Jones, 2010; Robinson et al., 2012; Dudas et al., 2016).

The application of network analysis has recently been extended to represent vortical interactions in fluid flows (Nair & Taira, 2015). The network-theoretic framework is comprised of discrete point vortices as nodes and interactions between them as edges. The network representation allows for the utilization of techniques such as spectral sparsification (Spielman & Srivastava, 2011) to identify key vortical interactions and maintain the spectral properties of such interactions. The application of such spectral sparsifiers has led to the development of sparsified-dynamics models, which preserve the invariants of discrete vortex dynamics. These models are capable of tracking bulk motion of vortical clusters in a computationally effective manner. Moreover, the extraction of the vortical network structure of turbulent flows has revealed the scale-free network property of decaying two-dimensional isotropic turbulence (Taira et al., 2016). The resulting framework enables the assessment of the resilience of turbulent flow structures. As such, network analysis has demonstrated the ability to capture fundamental vortical interactions in fluid flows, leveraging graph-theoretic principles. In the present work, we extend network analysis to describe and control modal interactions in fluid flows, by casting the fluid flow in terms of a networked oscillator system. In this approach, we view the modal oscillators as nodes and interactions between them as edges, highlighting complex energy transfer dynamics.

We utilize modal decomposition techniques in conjuction with coupled oscillator models to capture the interactive physics involved in unsteady fluid flows. Our objective in this work is three-fold; (1) characterize the nonlinear energy transfer between modes and construct a networked dynamics model for amplitude and phase perturbations in unsteady fluid flows, (2) describe interactive dynamics between modes from a network-theoretic perspective and (3) control the perturbations with respect to the limit cycle state of periodic flows as well as the full state itself. To accomplish these goals, mode pairs describing individual coherent structures in baseline (unforced and unperturbed) time-periodic flows are considered as a set of oscillators. A graph representation with the oscillators as nodes and coupling interactions between them as edges is constructed to form a set of networked oscillators.

Perturbations are introduced in the unsteady fluid flow to examine the associated energy transfer among the modes which are tracked over the networked oscillator based framework. A networked oscillator model herein describes the temporal dynamics of modes in unsteady fluid flow with a network structure embedded in it. The associated network structure is extracted using a linear regression procedure (Kutner et al., 2004). With the network dynamics model established, we are able to study the amplitude and phase dynamics of these perturbations on the modal interaction network. Thus, we arrive at a reduced-order network model using modal oscillators that highlight interaction dynamics among the modes which can subsequently be used for the control of the flow unsteadiness. We design flow control strategies to suppress modal fluctuations using the network based insights that consequently leads to drag reduction. In what follows, we first lay the theoretical foundation of this work in §2. We then demonstrate the strength of our approach using a canonical example of two-dimensional cylinder flow in §3. We end the paper with concluding remarks in §4.

2 Formulation

In this study, we are interested in modeling and controlling the unsteady fluctuations of flows with well-identified coherent structures captured as spatial modes. Various modal decompositions techniques can be used to extract coherent structures in fluid flows as mentioned in §1. In the present work, we are interested in capturing the energy transfer dynamics between those modes. In this regard, POD provides an optimal set of minimum number of modes capturing maximum energy content of the dynamical system. We first describe the procedure to extract these modes in unsteady fluid flows using the POD technique.

2.1 Oscillator representation

We consider a baseline case corresponding to time-periodic flow without any forcing or perturbations introduced in the Navier–Stokes equations. To perform the POD analysis, we use the method of snapshots (Sirovich, 1987). Using this approach, the unsteady velocity field $\boldsymbol{u}$ can be approximated by a finite expansion in terms of a mean (time-averaged) velocity field $\bar{\boldsymbol{u}}$ and $N$ orthonormal spatial POD modes $\boldsymbol{\phi}_{j}^{\boldsymbol{u}}$ as

[TABLE]

where $a_{j}(t)=\left<\boldsymbol{u}(\boldsymbol{x},t)-\bar{\boldsymbol{u}}(\boldsymbol{x}),\boldsymbol{\phi}_{j}^{\boldsymbol{u}}(\boldsymbol{x})\right>$ are the temporal coefficients and $\left<\cdot,\cdot\right>$ denotes the inner product over the computational domain. With $||\boldsymbol{\phi}_{j}^{\boldsymbol{u}}||=1$, the kinetic energy of the modes are given by $a^{2}_{j}/2$, amounting to a total modal energy of $E=\sum_{j=1}^{N}a^{2}_{j}/2$. As mentioned earlier, we obtain modes in conjugate pairs from POD for time-periodic flows. Each conjugate mode pair, describing periodic coherent structures in the baseline case define an oscillator in our analysis. From a set of $N$ POD modes, we obtain $M=N/2$ such oscillators. These oscillators can also be identified using the procedure described in the work of Sieber et al. (2016).

The temporal oscillations associated with each oscillator can be described by their coefficients in the complex plane. Conjugate mode pairs $(\boldsymbol{\phi}_{2j-1}^{\boldsymbol{u}},\boldsymbol{\phi}_{2j}^{\boldsymbol{u}})$ with temporal coefficients $(a_{2j-1},a_{2j})$ can be represented in the complex plane as,

[TABLE]

where $m=\text{I},\text{II}\dots M$, $j=1,2\dots N/2$, $r_{m}=|z_{m}|$, and $\theta_{m}=\angle z_{m}$. Throughout this work, the oscillators will be numbered in roman numerals, $\{$I, II$~{}...~{}M\}$ denoted by index $m$ to distinguish from mode numbering, $\{1,2,~{}...~{}N\}$ denoted by index $j$. Each oscillator $m$ is associated with a temporal coefficient defined on the complex plane ($z_{m}$) which consists of an odd-number mode with coefficient $a_{2j-1}$ and an even-number mode with coefficient $a_{2j}$. That is, modes $1$ and $2$ constitute oscillator I, modes $3$ and $4$ constitute oscillator II, and so on. The temporal coefficient corresponding to the mean flow ($\bar{\boldsymbol{u}}$) is fixed to unity, $z_{0}=1$. Equivalent to Eq. (1), we can recover the velocity field in terms of oscillators and their associated temporal coefficients as

[TABLE]

where $\Re(z_{m})$ and $\Im(z_{m})$ represent real and imaginary components of the temporal coefficient $z_{m}$, respectively.

In the baseline case, oscillators follow a natural limit-cycle, the generalized dynamics of which is described by the Stuart–Landau equation as

[TABLE]

with $z_{m}^{b}=r_{m}^{b}\exp(i\theta_{m}^{b})$, $\lambda_{m}=(\overline{r^{b}_{m}})^{2}$ and $\Omega_{m}^{b}$ is the oscillator frequency, where $\overline{\cdot}$ denotes the time average. Here, the superscript $()^{b}$ denotes the baseline case. Based on Eq. (4), the conjugate mode pairs of each oscillator exchanges energy among themselves to describe the self-sustaining equilibrium behavior of each coherent structure in unsteady fluid flows. Stable limit cycles are a characteristic feature of systems exhibiting such self-sustaining oscillations (Strogatz, 2014). This self-sustaining behavior may not be true for systems with non-periodic or chaotic attractors, such as turbulent flows.

The baseline temporal dynamics of the individual oscillators for periodic flow are independent of each other as given by Eq. (4) and does not involve any coupling between oscillators. To highlight interactions between oscillators, we introduce perturbations impulsively in DNS to the baseline temporal coefficients of oscillators at $t=t_{0}$. These perturbations introduced to the baseline shedding state leads to the emergence of nonlinear interactions and energy exchange among the oscillators resulting in fluctuations described by

[TABLE]

where $()^{\prime}$ denotes the perturbation quantity and $\epsilon_{m}^{\prime}(t_{0})$ and $\theta_{m}^{\prime}(t_{0})$ are the amplitude and phase perturbations for the $m$th oscillator, respectively. It should be noted that $\epsilon_{m}^{\prime}$ is normalized by the baseline amplitude of each oscillator.

The total temporal coefficient for each oscillator in the perturbed case can then be described by combining the baseline temporal coefficients and fluctuation as

[TABLE]

The overall oscillator model highlighting perturbations introduced is shown in figure 1 (a). The blue circle in the figure indicates the natural limit-cycle of oscillator $m$ with baseline temporal coefficient ($z_{m}^{b}$) lying on it. The perturbations in amplitude ($\epsilon_{m}^{\prime}$) and phase ($\theta_{m}^{\prime}$) result in a total temporal coefficient ($z_{m}$) off the limit cycle. We need to consider two main factors to characterize the perturbations introduced for each oscillator $m$; (1) the initial amplitude of perturbation $\epsilon_{m}^{\prime}(t_{0})$ and (2) the initial phase perturbation size $\theta_{m}^{\prime}(t_{0})$. Once these factors are determined, using Eq. (5) and (6) in Eq. (3), the initial velocity field for DNS of the perturbed case is prescribed.

Once the perturbations are introduced, the perturbation energy is distributed among the oscillators according to the flow physics. The temporal dynamics of the modes are attributed to the interactions between them resulting primarily from the convective term of the Navier–Stokes equation. These interactions cause energy transfers among the modes in unsteady fluid flow. As the POD modes are orthogonal to each other, the temporal coefficients in the perturbed case can be extracted by projection, as $a_{j}=(\boldsymbol{u}-\bar{\boldsymbol{u}},\boldsymbol{\phi}_{j}^{\boldsymbol{u}})$. We then construct the temporal coefficients in the complex plane $z_{m}$ for each oscillator $m$ using Eq. (2). To capture the fluctuating amplitude and phase of the temporal coefficients, we track the normalized fluctuation

[TABLE]

Thus, the normalized fluctuation tracks the amplitude $\epsilon_{m}^{\prime}(t)$ and phase $\theta_{m}^{\prime}(t)$ fluctuations of the modal oscillators due to interactions.

The amplitude fluctuation of the oscillators is related to the variation in oscillator energy $E_{m}^{\prime}(t)$ compared to the baseline, given by

[TABLE]

where $E_{m}^{\prime}(t_{0})$ represents the oscillator energy input, $E^{b}$ is the total baseline modal energy given by $E^{b}=\sum_{j=1}^{N}(a_{j}^{b})^{2}/2=\sum_{m=\text{I}}^{M}(r^{b}_{m})^{2}/2$ and $\beta_{m}$ represents the factor of input energy introduced to oscillator $m$ at initial time $t_{0}$ compared to the total baseline energy. If no amplitude perturbation is introduced in the simulation, $E_{m}^{\prime}(t_{0})=0$. If only a phase perturbation is introduced in the simulation, variation in oscillator energy may be observed due to oscillator interactions, although the oscillator energy input is zero. To be consistent with the tracking of normalized fluctuation $\zeta_{m}$ in Eq. (7), the normalized quantities $E_{m}^{\prime}/(r^{b}_{m})^{2}$ and $a_{j}^{\prime}/r^{b}_{m}$ are used to track the energy perturbations of the oscillators and amplitude perturbations of the modes, respectively while $\theta_{m}^{\prime}$ tracks the phase of the oscillators. This normalization is adopted to facilitate the tracking of the energy transfer between modes. Multiple combinations of oscillator energy input and phase perturbations can be used to quantify energy exchange for each perturbed case.

2.2 Networked oscillator representation

Using the oscillators extracted from POD and their normalized fluctuations extracted from the impulse response about the baseline flow, we can now create a network-theoretic representation of the unsteady fluid flow system. In the present work, we consider the oscillators (modes pairs) characterizing unsteady fluctuations to be the nodes ($\mathcal{V}$) of the network and the interaction between them as edges ($\mathcal{E}$). To characterize the interactions between oscillators, we need a model that entails coupling between oscillators. This motivates a more general networked oscillator model of $M$ linearly coupled oscillators given by

[TABLE]

where $m=\text{I},\text{II}\dots M$, $[\boldsymbol{A}_{\mathcal{G}}]_{mn}$ and $[\boldsymbol{L}_{\mathcal{G}}]_{mn}$ are the adjacency and the in-degree graph Laplacian matrices, respectively. The networked oscillator model considers propagation processes on networks and is used as a simple model for examining the spread of kinetic energy among oscillators. It should be emphasized that this is a linear model.

The dynamics of normalized fluctuations of temporal coefficients are nonlinear in general due to convective physics in unsteady fluid flows. The nonlinear change of variables into oscillator magnitude and phase coordinates enables this linear representation of the nonlinear flow dynamics. In general, coordinate transformations that embed nonlinear dynamics in a linear framework are related closely to the Koopman operator (Koopman, 1931; Mezić, 2005, 2013). Obtaining useful coordinate transformations that approximate intrinsic Koopman coordinates remains an open challenge in data-driven dynamical systems (Kutz et al., 2016; Brunton et al., 2016a). Current data-driven methods to identify Koopman eigenfunctions include the variational approach of conformation dynamics (VAC) (Noé & Nuske, 2013; Nüske et al., 2016), the equivalent extended dynamic mode decomposition (EDMD) (Williams et al., 2015a, b; Klus et al., 2016), and the use of delay coordinates (Brunton et al., 2017; Arbabi & Mezić, 2016). The proposed transformation into oscillator coordinates in the present study provides a promising approximate linear representation of nonlinear dynamics.

For a set of $M$ oscillators, the above adjacency matrix $\boldsymbol{A}_{\mathcal{G}}\in\mathbb{C}^{M\times M}$ concisely describes the network connectivity given as

[TABLE]

In our formulation, no self-loops are present in the network connectivity of oscillators. This is primarily attributed to the nature of the networked oscillator model relying on interactions between oscillators to describe the dynamics of oscillator $m$. This is reflected from Eq. (9), for $n=m$, the contribution to the dynamics of oscillator $m$ is zero. The rows of the adjacency matrix indicates the dependence of the oscillators $n$ on the dynamics of oscillators in column $m$, i.e., $w_{mn}$. The in-degree ($k_{m}$) of a node $m$ represents the summation of the incoming weights of the edges connected to it given by $k_{m}=\sum_{n=\text{I}}^{M}[{\boldsymbol{A}}_{\mathcal{G}}]_{mn}$. The (in-degree) graph Laplacian is closely related to the adjacency matrix as $\boldsymbol{L}_{\mathcal{G}}=\boldsymbol{D}_{\mathcal{G}}-\boldsymbol{A}_{\mathcal{G}}$, where $\boldsymbol{D}_{\mathcal{G}}$ is a diagonal matrix with elements equal to the in-degree of the nodes, $\boldsymbol{D}_{{\mathcal{G}}}=\text{diag}([k_{m}]_{m=\text{I}}^{M})$. The graph Laplacian is a discrete analog of the continuous Laplacian operator. In the current study, it encodes structural properties of the networked oscillator model. Propagation (diffusion) processes on networks are neatly described by this graph Laplacian.

Using the knowledge of time evolution of normalized fluctuations from DNS for each oscillator, (i.e., $\boldsymbol{\zeta}$ and $\dot{\boldsymbol{\zeta}}$), we use a linear regression procedure to determine the adjacency matrix weights. These weights could also be obtained through a Galerkin regression approach (Loiseau & Brunton, 2016). Linear regression allows us to examine the relationships between predictive (independent) and response (dependent) variables. It enables us to find the best possible linear model fit for estimating the response of a system by regressing the dynamics onto the predictor variables. Here, the time derivative of the normalized fluctuations of the individual oscillators ($\dot{\zeta}_{m}$) are the response variables and $\left(\zeta_{n}-\zeta_{m}\right)$ are the predictor variables. The result of such a multiple regression over $m$ oscillators yields the adjacency matrix weights for the networked oscillator model. It must be noted that as the temporal coefficient associated with the mean flow is constrained to unity, $z_{0}^{\prime}=0$ and does not contribute to normalized fluctuations of the modes. Thus, mean flow is considered as an isolated node in this formulation.

The oscillator network representation is illustrated in figure 1 (b), where the oscillators are taken to be nodes of the network. The oscillator model in figure 1 (a) highlights the temporal coefficient dynamics associated with the nodes of the network, e.g., oscillator II (shown in blue). The complex adjacency weights obtained from linear regression $w_{mn}$ can be decomposed into a magnitude ($|w_{mn}|$) and phase ($\angle{w_{mn}}$). The magnitude of the edge weights signify the influence of oscillator $n$ on oscillator $m$ as illustrated in figure 1 (b). The phase of the edge weights represent the individual modal contributions in oscillator phase interactions. In particular, it highlights the phase advances or delays imposed between modes of interacting oscillators. The odd-odd mode interactions and even-even mode interactions are given by $|w_{n,m}|\cos{\angle{w_{n,m}}}$ while the odd-even and even-odd mode interactions are given by $|w_{n,m}|\sin{\angle{w_{n,m}}}$ and $-|w_{n,m}|\sin{\angle{w_{n,m}}}$, respectively.

For comparison with the networked oscillator model, we also construct the Galerkin projection model for the incompressible Navier-Stokes equations, by projecting the equations onto the POD modes to construct the POD-Galerkin reduced-order model (Noack et al., 2011; Holmes et al., 2012). The resulting Galerkin system is given by

[TABLE]

where $\psi_{jk}$ and $\chi_{jkl}$ are the linear and quadratic coefficient terms, respectively. Here, as mentioned before, the temporal coefficient associated with the mean is fixed to unity, $a_{0}=1$. The linear term contains diffusive physics of the modes while the quadratic term contains the convective physics. The networked oscillator model captures the effect of both diffusion and convection in the network structure that is extracted. Both the network and Galerkin models require only the specification of initial condition of the temporal coefficients of the mode pairs only.

The networked oscillator model, as discussed above, describes both amplitude and phase dynamics of collective oscillation and captures the energy exchange and phase dynamics among the modes. As discussed earlier, the present model is derived from the analysis of impulse response in DNS. In §3, we describe the application of the oscillator network-based formalism to a canonical problem of two-dimensional incompressible flow over a cylinder. The model can also serve as a foundation to perform network-informed feedback control of unsteady flow as will be discussed in §3.4.

3 Application to cylinder flow

Two-dimensional flow over a circular cylinder is a canonical problem in fluid mechanics and serves as a model to capture fundamental wake dynamics for many flows. The coherent structures identified are well-reported in literature, and thus cylinder flow forms an excellent testbed for our analysis. Moreover, for the cylinder flow problem, POD modes are well characterized and are the same as DMD modes as demonstrated in the work by Chen et al. (2012). Thus, as an alternative to POD for the current approach, DMD can be also considered. To ensure orthogonality of the modes required for the current approach, recursive DMD can be performed to obtain orthogonal pure frequency modes (Noack et al., 2016).

3.1 Computational approach

For this study, we perform DNS of incompressible flow past a circular cylinder using the immersed boundary projection method (Taira & Colonius, 2007; Colonius & Taira, 2008; Kajishima & Taira, 2017) at a diameter-based Reynolds number of $Re=100$. This method employs a Cartesian grid with the immersed body formulation to generate the cylinder. We take advantage of the multi-domain technique to simulate cylinder flow in free space. The innermost domain is chosen as $-1\leq x/d\leq 29,-15\leq y/d\leq 15$ with a resolution of $600\times 600$ grid points where $d$ is the cylinder diameter. Here, $x$ is the streamwise direction and $y$ is the cross-stream direction. The outermost domain is chosen as $-16\leq x/d\leq 44,-30\leq y/d\leq 30$, far enough so as to not affect the results in the near-field. Uniform flow is prescribed at the domain boundaries. For time integration, the method uses an implicit Crank-Nicholson scheme for the viscous term and an Adam-Bashforth discretization for the convective term.

From the simulation, the drag coefficient $(C_{D})$ and lift coefficient $(C_{L})$ are computed as

[TABLE]

where, $F_{D}$ and $F_{L}$ are the drag and lift forces on the cylinder, $\rho$ is the freestream density, and $U$ is the freestream velocity field. These are obtained by integrating the boundary force representing the cylinder in the immersed boundary formulation. The Strouhal number for the flow is defined as $St=f_{n}d/U$, where $f_{n}$ is the natural shedding frequency. We obtain a Strouhal number of $St=0.164$, drag coefficient $C_{D}=1.35\pm 0.009$ and lift coefficient $C_{L}=\pm 0.325$ which agree well with the literature (Liu et al., 1998; Taira & Colonius, 2007; Canuto & Taira, 2015). The flow exhibits vortex shedding behavior in the cylinder wake, as shown by the instantaneous vorticity field in figure 2 (a). Such vortex shedding characterizes a von Kármán vortex street due to the repetitive pattern of vortices in the unsteady wake. The time-averaged (mean) vorticity field is shown in figure 2 (b). We then perform POD with the method of snapshots (Sirovich, 1987) using the velocity field data ($\boldsymbol{u}$) gathered from DNS. The modes and temporal coefficients obtained are in agreement with those from the work by Noack et al. (2003).

3.2 Baseline (unperturbed flow)

The extracted spatial modes and associated temporal coefficients ($a_{j}^{b}$) are shown in figure 2 (c) and (d), respectively, in terms of the vorticity field, $\boldsymbol{\phi}_{j}^{\boldsymbol{\omega}}=\nabla\times\boldsymbol{\phi}_{j}^{\boldsymbol{u}}$. We define the conjugate mode pairs as independent oscillators in our formulation. Oscillator I constitutes mode pair ($1,2$), oscillator II constitutes mode pair ($3,4$) and so on. The mode pairs (oscillators) are ordered in terms of decreasing energy content, $E_{j}=(a_{j}^{b})^{2}/2$, shown in figure 2 (top to bottom). As the first eight POD modes capture $99.98\%$ of the modal kinetic energy, we choose $N=8$ for our analysis. The energy content of oscillators I, II, III and IV are $96.87\%,2.18\%,0.88\%$ and $0.05\%$, respectively. The oscillator temporal coefficients in the complex plane ($z_{m}^{b}$) corresponding to temporal coefficients of the individual mode pairs from Eq. (2) are shown in figure 2 (e). For this canonical problem, the frequency associated with the higher-order oscillators are harmonics of the lowest-order oscillator, $\Omega_{m}^{b}=m\Omega_{\text{I}}^{b}$. The frequency of oscillator I corresponds to the Strouhal number associated with the natural shedding cycle of the flow, $\Omega_{\text{I}}=0.164$. As the frequency ($\Omega_{m}^{b}$) associated with the temporal coefficients of the mode pairs increases, the size of the spatial modal structures decreases. Again, the oscillators in the baseline flow are associated with limit cycle temporal dynamics, which can be described by Eq. (4), independent of each other with no coupling between them. This lack of coupling in the generalized limit-cycle dynamics does not capture interactions between oscillators in unsteady fluid flow.

3.3 Perturbed flow

To capture interactions between oscillators, additional energy and phase perturbations are impulsively introduced to the simulation through the initial condition. These perturbations cause added fluctuations in the temporal coefficients of the modes. The projected coefficients $a_{n}=\left<\boldsymbol{u}-\bar{\boldsymbol{u}},\boldsymbol{\phi}_{n}^{\boldsymbol{u}}\right>$ from the perturbed case are extracted from DNS and the normalized fluctuation $\zeta_{m}(t)$ is tracked using the networked oscillator model discussed in §2. We expect that as the perturbation convects downstream, the normalized fluctuation of the oscillators will decay to zero and the perturbed field will return to the baseline limit-cycle behavior. However, the addition of a perturbation to the flow results in a constant phase shift of each oscillator in the final limit cycle, compared to the unperturbed limit cycle.

In two-dimensional unsteady cylinder flow, the leading POD mode pair corresponding to oscillator I contains maximum energy and any deviations in the final limit cycle is immediately detected by the oscillator I phase ($\theta_{\text{I}}$). Thus, to construct the normalized fluctuation time history for each oscillator, we align the phase of oscillator I for the perturbed case and the baseline case as $\zeta_{m}=(z_{m}-z_{m}^{b}|_{\theta=\theta_{\text{I}}})/z_{m}^{b}|_{\theta=\theta_{\text{I}}}$. Once the normalized fluctuation $\zeta_{m}$ is constructed, we also compute its time derivative $\dot{\zeta}_{m}$. We then construct the library of functions ($\zeta_{n}-\zeta_{m}$) for each oscillator $m$. We perform a simple linear regression procedure to obtain the network structure $[\boldsymbol{A}_{\mathcal{G}}]_{mn}$. Once the coupling function (network structure) is obtained, we solve the linear networked oscillator model in Eq. (9) for prediction with a prescribed initial condition and compare the fluctuations with those obtained in DNS. While it is not necessary, one could also consider interaction terms that are quadratic or higher (Brunton et al., 2016b; Loiseau & Brunton, 2016).

3.3.1 Amplitude perturbations

We first introduce an amplitude perturbation to oscillator II, which is associated with the first harmonic of the natural shedding frequency of the flow. Introducing perturbations at this frequency perturbs not only the natural shedding frequency but also the oscillators associated with higher harmonics of the flow due to interactions present in the flow field. Let us consider the addition of $20\%$ of baseline modal energy to oscillator II ($\beta_{\text{II}}=0.2$). No perturbation in phase is added, i.e., $\theta_{\text{II}}^{\prime}(t_{0})=0$. With these initial perturbations determined, we can construct the initial flow field required for DNS of the perturbed case. The initial vorticity field, as shown in figure 3 (a) for the perturbed case, is different from figure 2 (a) for the baseline. Here, as oscillator II is energized, we notice the prevalence of spatial vortical structures associated with mode pair ($3,4$). The vorticity magnitude associated with these structures increases due to the added energy and correspondingly appears in the initial vorticity field for the perturbed case. As the introduced perturbation in oscillator II convects downstream in the numerical simulation, it interacts with other oscillators. Using the procedure highlighted above, we construct the normalized fluctuations $\boldsymbol{\zeta}$ and extract the adjacency matrix $\boldsymbol{A}_{\mathcal{G}}$ that highlights network interactions.

The magnitude and phase of the adjacency matrix are shown in figure 3 (b). In this example, the dynamics of oscillator I (leading oscillator) is not affected noticeably by the other oscillators. The dynamics of oscillators II and III show high levels dependence on oscillator I. This is consistent with our expectation as most of the energy in the flow is held in oscillator I and passes down to higher-order oscillators. Also, if we evaluate the strength of the incoming edge weights of each oscillator ($\sum_{n}w_{mn}$) for this perturbed case, oscillator III shows maximum influence from the other oscillators. The complex matrix structure of the extracted network summarizes the interactions in the flow upon introduction of impulse perturbations. To reveal oscillator phase interactions, we evaluate the phase of the edge weights in figure 3 (b). For the dynamics of oscillator II, i.e., mode pair $(3,4)$, the phase of the edge weights indicates that interactions between modes $1\rightarrow 3$, $2\rightarrow 4$, $8\rightarrow 3$ and $7\rightarrow 4$ are larger than other modal interactions. For this perturbed case, in general, phase advancing effects are predominant for the oscillators affecting each other.

Once the network structure is obtained, we can solve the linear networked oscillator model with the prescribed initial condition corresponding to the perturbation. For comparison, we also construct the empirical Galerkin reduced-order model. The oscillator dynamics from DNS (reference) and those predicted from the networked oscillator and POD-Galerkin models are shown in figure 3 (c). The initial amplitude of the perturbed temporal coefficients for oscillator II corresponds to the red dot, $r_{\text{II}}(t_{0})$. However, no phase perturbation is introduced in this example. The initial amplitude and phase of the other oscillators are unchanged, as no perturbations are introduced in these oscillators. It can be seen in figure 3 (c) that the oscillator dynamics in the networked oscillator model agree well with DNS trajectories, particularly the dynamics of oscillators I, II and III. The networked oscillator model and DNS show negligible differences in tracking the trajectory of oscillator II from the perturbed initial condition to its final limit cycle. The model trajectory of oscillator IV shows some similarities with DNS. Due to its low energy content, any small deviations in the dynamics of lower-order oscillators causes comparable changes in the trajectory of oscillator IV. On the other hand, the POD-Galerkin model overpredicts the fluctuations in oscillators II, III and IV. It can be seen that the networked oscillator model tracks the fluctuations better than the POD-Galerkin model.

To further compare the details of the predicted trajectories, we track the fluctuations in modal amplitude and oscillator energy in figure 4 (a) and (b), respectively. The networked oscillator model shows excellent agreement with DNS, tracking the amplitude perturbation $(a_{2n-1}^{\prime})$ and energy transfers $(E_{m}^{\prime})$. Moreover, we see agreement in the long-time behavior of the fluctuations as the flow returns to the baseline state for the networked oscillator model. In contrast, the POD-Galerkin model is not well-designed for modeling the long-time behavior of modal fluctuations and hence is not expected to work well as time progresses. This is especially true as POD modes are used for the empirical Galerkin formulation. As indicated by the green dashed line, the POD-Galerkin model overpredicts these fluctuations and their associated time scales.

3.3.2 Phase perturbations

Similar to tracking amplitude perturbations, we can also analyze the flow response to phase perturbations. These phase perturbations are introduced by shifting the initial phase of the oscillators at $t_{0}$. As energy input is related to the amplitude of the perturbations, phase perturbations do not introduce any additional energy to the flow. Here, let us introduce a phase perturbation of size $\theta_{\text{II}}^{\prime}(t_{0})=-\pi/2$ by shifting the phase associated with oscillator II at $t_{0}$ in the baseline case in the clockwise direction by $\pi/2$. The initial vorticity field resulting from the phase perturbation is shown in figure 5 (a). Over time, both the amplitude and phase of the oscillators are affected with the introduction of phase perturbations in DNS. These fluctuations are tracked and the network structure is extracted as shown in figure 5 (b). As seen from the magnitude of the edge weights, similar to the amplitude perturbation case, oscillator I dynamics are negligibly affected by the other oscillators. Furthermore, oscillator IV does not contribute much to the dynamics of oscillators II and III. Deductions of modal interactions can be made from the phase corresponding to the edge weights.

The oscillator dynamics for the phase perturbation case are shown in figure 5 (c). We notice that the phase perturbation introduced in oscillator II at $t_{0}$ results in a phase shift to an initial condition indicated by the red dot. Similar to the amplitude perturbation case, the dynamics in the networked oscillator model show better agreement with the DNS than the POD-Galerkin model. This is confirmed on tracking the amplitude fluctuations $(a_{2n-1}^{\prime})$ of the individual modes and corresponding phase dynamics $(\theta_{m}^{\prime})$ as shown in figure 6 (a) and (b), respectively. The POD-Galerkin model is unable to correctly track the phase of the oscillators and hence overpredicts the amplitude of the fluctuations.

3.3.3 Aggregate network model

The perturbed flow cases examined above are limited to amplitude and phase perturbations introduced to oscillator II only. Similar analyses can be performed for oscillators I, III and IV with various combinations of amplitude and phase perturbations introduced to each oscillator. For each case, the network structure of interactions can be extracted individually by tracking the trajectory of the oscillator fluctuations from DNS as demonstrated in the approach above. While the corresponding linear networked oscillator models can be built for each specific case, the resulting network structure obtained depends on the initial perturbations introduced. Thus, one may argue that such individually tuned model does not necessarily capture the overall interactions in a generally effective manner for any perturbation. Such a general network interaction model is desired not only to capture interactions for any combination of perturbations but also to design effective strategies for flow control. In this section, we describe our approach to build an aggregate network model that captures the fluid flow response to a variety of perturbation inputs.

We consider a range of perturbed flow cases by varying the amplitude and phase perturbation size and to which combination of oscillators to perturb. For a particular oscillator perturbed, we vary the energy input ($\beta_{m}$) from $0.05$ to $1$ ($5\%$ to $100\%$ additional baseline modal energy in the flow). We also vary the phase perturbation size ($\theta_{m}^{\prime}$) from $-2\pi/3$ to $2\pi/3$ to cover a broad coverage of plausible initial perturbations introduced. The range of amplitude and phase perturbations considered here are quite large. These large perturbations lead to the emergence of strong nonlinear interactions between the oscillators. However, we are interested in extracting a linear networked oscillator model to approximate these nonlinear interactions for the purpose of control. For the multitude of perturbation cases, we collect the oscillator fluctuation data from DNS. For the aggregate network model, we concatenate the trajectories of all the perturbed cases obtained from DNS. We then compute the normalized fluctuation $\zeta_{m}$ and its time derivative $\dot{\zeta}_{m}$, instead of individually tracking perturbations and monitoring the oscillator fluctuations in each perturbed case.

We then segregate the combined data into training and test sets for performing cross-validation and evaluating the predictive capabilities of the network model. Varying fractions of the combined input-output data are randomly chosen as the training sets. For each training set, regression analysis is performed to extract a corresponding network model. We then examine the in-degree ($k_{m}=\sum_{n=\text{I}}^{M}[{\boldsymbol{A}}_{\mathcal{G}}]_{mn}$) and out-degree ($k_{n}=\sum_{m=\text{I}}^{M}[{\boldsymbol{A}}_{\mathcal{G}}]_{mn}$) of the network nodes for each model extrated. The variation of the network degrees with respect to the fraction of the chosen training set is shown in figure 7 (a). We observe that as the fraction of the training data used increases, convergence of the network degree is obtained. We also learn that the average in-degree increases from low-order oscillators to higher-order oscillators and the average out-degree decreases. Oscillators I and IV have maximum out-degree and in-degree respectively so that oscillator I influences the other oscillators most while oscillator IV is the most influenced. Oscillators II and III have a more balanced in and out-degree indicating more balanced energy transfer for each mode.

Based on the network degree convergence, a fraction of training data $f_{\text{train}}=0.8$ corresponding to $80\%$ of the combined input-output data is randomly chosen as the training set and the remaining $20\%$ of the data is used to assess the performance of the network model extracted. Using the training set, the aggregate network model extracted is shown in figure 7 (b). The magnitude of the network reveals that the lower-order oscillators have more influence on the dynamics of the higher-order oscillators. This follows from our earlier discussion regarding the network node degrees and also agrees with our intuition that energy is passed from lower-order oscillators to higher-order oscillators. Alternatively, we can construct the aggregate network model by considering complete trajectories of $80\%$ of randomly chosen perturbed cases yielding a very similar aggregate network model.

We also compare the aggregate network model with network models built with non-randomly chosen training sets. In particular, we consider trajectories of perturbed cases of individual oscillators. Combining the perturbed trajectories for various amplitude and phase perturbation sizes for each oscillator, we build network models corresponding to oscillators I, II, III and IV. We use the same test data set as before to assess the performance of each of these models along with the aggregate model. We predict the time derivative of the normalized fluctuation $\dot{\zeta}_{m}$ with each of the models on the test data set and compare with DNS reference values. In figure 7 (c), we assess the prediction error through a normalized root mean square deviation ($\triangle_{m}$), given by

[TABLE]

where $n_{t}$ is the length of test data set. Here, $\dot{\zeta}_{m}(k)$ is the time derivative of oscillator $m$ obtained from DNS and $\hat{\dot{\zeta}}_{m}(k)$ is the predicted time derivative value based on oscillator network model for the $k$th test data set. The root mean square deviation is normalized by the range of the measured data. We can observe from figure 7 (c) that the aggregate model achieves the least error in predicting the normalized fluctuation time derivative. As the aggregate model contains more information of the general interactions based on various oscillator perturbations, it yields enhanced predictive capabilities for energy transfer amongst all oscillators.

Using the aggregate network structure shown in figure 7 (b), we evaluate the performance of the networked oscillator model in tracking the fluctuations of the perturbation cases examined in this work. We compare the DNS trajectories of the oscillator dynamics with the aggregate network oscillator model with introduction of amplitude perturbations corresponding to $\beta_{m}=0.2$ in figure A1 in the Appendix. The aggregate networked oscillator model is able to track the amplitude and phase dynamics for each case with reasonable accuracy. We can see that the introduction of an amplitude perturbation in oscillator I does not cause large deviations in the higher-order oscillators. Moreover, the prescribed energy input requires only slight deviation from the limit cycle orbit for perturbations introduced on oscillator I as it possesses high level of energy to begin with. A perturbation introduced to oscillator III results large fluctuations in oscillators II and IV, while a perturbation introduced to oscillator IV results in comparable fluctuations in oscillator III but only small deviations in oscillator I and II dynamics. These features are well captured by the aggregate networked oscillator model.

We also show the comparison of predicted trajectories of the aggregate networked oscillator model for phase perturbation cases shown in figure A2 of the Appendix. As oscillator I contains large amounts of energy, a phase perturbation of $\theta_{\text{I}}^{\prime}$ causes an effective phase perturbation of $m\theta_{n}^{\prime}$ for oscillators II, III and IV. Thus, the phase perturbation on oscillator I induces a combined phase perturbation on the higher-order oscillators as they synchronize their phases with respect to oscillator I. Phase perturbations introduced in oscillator III cause fluctuations in oscillators II and IV, while phase perturbations introduced in oscillator IV do not produce any fluctuations in lower-order oscillators. As with the amplitude perturbation cases, the aggregate network model also captures the features associated with the phase perturbation cases. Thus, one aggregate network structure is able to capture a multitude of perturbation cases, yielding a global oscillator interaction model.

The aggregate network model is not only useful to gain insights into energy transfer, but is also used to design feedback control laws to manipulate the flow unsteadiness. In this work, we only characterize the interactions between the baseline POD modes in the flow induced by the perturbations introduced. However, these perturbations result in the emergence of residual flow structures not captured by the baseline POD modes. Appendix B discusses how to capture the emerging modes. Based on the modal interactions characterized, we design flow control strategies to suppress modal fluctuations. In the next section, we discuss such flow control design strategies using the networked oscillator model in conjunction with a linear quadratic regulator (LQR), and discuss how the residual behavior not captured by the modal dynamics influence the controlled flow.

3.4 Feedback control

With the modal network interactions characterized and the aggregate network model developed in the previous section, we now design a feedback control strategy to suppress modal oscillations. As mentioned previously, modal oscillations in the flow are reflected in the temporal coefficients associated with the modes. Suppression of these modal oscillations is critical in reducing the wake unsteadiness. Roshko (1954) and Mao et al. (2015) reported the strong relationship between unsteadiness in the wake and the drag force acting on the bluff body. If the application of control forces the modal temporal coefficients to zero, the flow will approach the mean flow. However, the mean flow is not a steady solution to the Navier–Stokes equation in general, and hence the flow diverges from the mean towards the unstable steady state in this cylinder flow case.

In the seminal work of Noack et al. (2003), it was shown that the mean flow and the unstable steady state are connected by a shift mode. The shift mode for the cylinder flow problem is shown in figure 8 (top left). This shift mode captures the transient dynamics between the onset of vortex shedding near the unstable steady state shown in figure 8 (bottom left) and the baseline mean flow of the periodic von Kármán vortex street in the globally stable limit cycle. This evolution takes place along the parabolic inertial manifold shown in grey in figure 8. Here, the limit cycle of oscillator I (modes $1$ and $2$) is illustrated in the $a_{1}$-$a_{2}$ plane. The vertical axis represents the change in the temporal dynamics associated with the shift mode (as shown in red) with application of control. The drag coefficient also varies between the mean flow and unstable fixed point. A minimum drag coefficient is attained at the unstable steady state $C_{D}^{*}$, which also gives zero lift force on the cylinder. The kinetic energy associated with the shift mode varies as $a_{\triangle}^{2}$ and the drag force on the cylinder in the mean shift regime scales as $\sqrt{C_{D}^{\prime}}\propto a_{\triangle}$, as shown in figure 8 (blue) where, $C_{D}^{\prime}=C_{D}-C_{D}^{*}$. Thus, a realization of mean shift towards the unstable steady state yields a reduction in drag force.

We first develop a low-dimensional control framework based on the networked oscillator model to attenuate perturbations in the flow. We then extend the formulation to suppress the overall flow unsteadiness. In the present work, the aggregate networked oscillator model is a linear ordinary differential equation given by Eq. (9) governing the modal perturbation dynamics in the flow. As the networked oscillator model is linear, we can exploit the use of the LQR to control the modal fluctuations. Adding a forcing input to the networked oscillator model (in vector form), we arrive at

[TABLE]

where $\boldsymbol{\zeta}=[\zeta_{\text{I}},\zeta_{\text{II}}\dots\zeta_{M}]^{T}$, $\boldsymbol{v}\in\mathbb{C}^{M\times 1}$ is the forcing input and $\boldsymbol{B}\in\mathbb{R}^{M\times 1}$ is the actuation input matrix. The $m^{th}$ entry of $\boldsymbol{B}$ corresponds to forcing added to oscillator $m$.

We implement optimal full-state control with $\boldsymbol{v}=-\boldsymbol{K}\boldsymbol{\zeta}$ such that

[TABLE]

where the gain matrix $\boldsymbol{K}$ is determined from the Riccati equation for LQR. An optimal control strategy using LQR minimizes the quadratic cost function of the form,

[TABLE]

where $\boldsymbol{Q}$ and $\boldsymbol{S}$ are the state deviation and input usage weights respectively. In the control studies considered henceforth, we set $\boldsymbol{Q}=\boldsymbol{I}$ and $\boldsymbol{S}=\sigma\boldsymbol{I}$ and consider a range of values for $\sigma$. The case corresponding to $\sigma=1$ weighs the state deviation cost and input usage cost equally while $\sigma<1$ weighs the state deviation cost more.

We can force the individual oscillators or a combination of oscillators (reflected in $\boldsymbol{B}$). To aid the selection of which oscillators to force, we examine the movement of the eigenvalues of the Laplacian matrix with and without control for different combinations of input matrix $\boldsymbol{B}$ as shown in figure 9 (a). Eigenvalues of Laplacian matrix ($-\boldsymbol{L}_{\mathcal{G}}$) reveal the dynamical characteristics of the system. As four oscillators are considered in this work, we obtain four eigenvalues ($\lambda_{\text{I}}\dots\lambda_{\text{IV}}$) denoted by the use of star symbols in figure 9 (a) for varying $\sigma$. Note that $\lambda_{m}\leq 0$ for all eigenvalues which is a characteristic feature of Laplacian-based systems. It should be noted that the presence of $\lambda_{m}=0$ poses difficulties in computing the empirical Grammians to quantify the degrees of controllability. This issue requires careful assessment for network-based control design and will be investigated in upcoming studies.

We compare the response for input matrix $\boldsymbol{B}$ for both single oscillator input and some multiple oscillator input cases. For each input case, we determine the LQR gain matrix $\boldsymbol{K}$ for values of $\sigma$ ranging from $0.1$ to $1000$ and examine the movement of eigenvalues of ($-\boldsymbol{L}_{\mathcal{G}}-\boldsymbol{B}\boldsymbol{K}$) which govern the behavior of the controlled system. For the single-oscillator input cases, we observe large movement in the eigenvalue corresponding to the forced oscillator as $\sigma$ is decreased, i.e., the real part of $\lambda_{m}$ decreases when more forcing input is provided to oscillator $m$. We also notice that an input in oscillator I affects all eigenvalues while inputs to the higher-order oscillators do not move the $\lambda_{\text{I}}$ eigenvalue. This is expected as oscillator I has maximum out-degree and correspondingly has the highest influence in the network.

As input on oscillator I is required to affect the $\lambda_{\text{I}}$ eigenvalue, we consider multiple oscillator input cases including oscillator I. A noticeable movement in the eigenvalues is observed with inputs on oscillators I and IV. For forcing inputs added to all oscillators, only $\lambda_{\text{I}}$ eigenvalue is affected with no influence on the other system eigenvalues. To summarize the effectiveness of forcing input on controlling the system behavior, we track the $\text{max}(\Re(\lambda))$ for multiple oscillator input cases in figure 9 (b). We observe that for small $\sigma$, oscillator input combinations of I and IV outperform the other input combinations while, for mid-range values of $\sigma$, oscillator combinations of I and III exhibit better control performance.

We design the LQR controller such that we move system eigenvalues towards the left side of the complex plane. Thus, we select the input matrix $\boldsymbol{B}=[1~{}0~{}0~{}1]^{T}$ which adds forcing input to oscillators I and IV. Considering full-state feedback control with this choice of $\boldsymbol{B}$, we compute the control gain matrix $\boldsymbol{K}$ using LQR. In the following, we examine two control scenarios. First, we attenuate modal perturbations introduced in the flow. Then, we apply the LQR-based feedback control to suppress all modal amplitudes altogether. This entails reducing the unsteadiness in the flow due to modal oscillations and achieve drag reduction.

3.4.1 Disturbance rejection

We first illustrate the control of modal disturbances introduced to the cylinder flow. This entails the control of fluctuations in the modal temporal coefficients ($z_{m}^{\prime}$). To demonstrate the control, we consider random amplitude perturbations added to all oscillators with $\text{max}(\beta_{m})=0.02$. The temporal coefficients of the oscillators for this perturbed case obtained from DNS are shown in figure 10 (top). The networked oscillator control system given by Eq. (15) describes control of temporal coefficients associated with the modes. To implement control design in DNS, a body force corresponding to $-\boldsymbol{B}\boldsymbol{K}\boldsymbol{\zeta}$ is added to the momentum equation with the spatial modal information is incorporated along with the temporal coefficients. The addition of body force is The results with the application of control in DNS for $\sigma=0.1$ are shown in figure 10 (bottom). We can drive the oscillators to the natural limit cycle much faster with control. We notice that as control is introduced in oscillators I and IV, the effectiveness of control is more pronounced with these oscillators. This also follows from figure 9 (a) where the correspondence between the Laplacian eigenvalues and oscillator inputs was discussed. Similar analysis can be performed for controlling a multitude of modal perturbations introduced.

3.4.2 Total oscillation control

In this section, we consider the control of total modal oscillations in the flow associated with wake unsteadiness. This requires the control of the temporal coefficients associated with the oscillators ($z_{m}$) instead of just the modal fluctuations ($z_{m}^{\prime}$). Though the oscillator interactions characterized in this work are based on fluctuations with respect to the natural limit cycle (baseline) state, we consider these interactions to be characteristic of the modal oscillations in the flow. In fact, the energy transfer mechanism should be similar for $z_{m}^{\prime}$ and $z_{m}^{\prime}+z_{m}^{b}$. We realize that by suppressing $z_{m}$, the flow is attracted towards the unstable steady state. Since we can assume that the above networked oscillator based LQR control is applicable near the baseline limit cycle, we expect that there is some region of validity of control to achieve drag reduction. Inhibiting nonlinear energy transfer should remove the energy input to sustain wake oscillations. Hence, we expect to use the same aggregate model extracted to suppress the limit cycle along with fluctuations of the modes to some degree.

Our objective is to control the temporal coefficients associated with the modes, yielding a control system of the form,

[TABLE]

In DNS, this amounts to adding a body force corresponding to $-\boldsymbol{B}\boldsymbol{K}\boldsymbol{z}$. Note that this control is similar to perturbation control, except here the target state for the oscillator temporal dynamics is zero as opposed to the baseline limit cycle. Hence, this control input steers the modal amplitudes corresponding to the baseline limit cycle to zero as shown in figure 11 (a). Here, we show the control performance for $\sigma=0.1$ and $1$. With control, the temporal coefficient associated with oscillator I first decays to zero followed by the higher-order oscillators. The energy lost by oscillator I is compensated with an initial increase in energy associated with oscillator IV. This is expected as oscillator IV is the most influenced node in the network. This highlights the energy transfers from lower-order oscillators to higher-order oscillators. As control input is added to oscillator IV, we suppress the corresponding modal oscillations. We show the input cost of control in figure 11 (b). We see that the input cost for control, $J_{\text{in}}=\int_{0}^{\infty}\boldsymbol{v}(t)^{T}\boldsymbol{S}\boldsymbol{v}(t)\rm{d}{\it t}$, is lower for the $\sigma=0.1$ case compared to the $\sigma=1$ case.

We show the suppression of total modal energy with control in figure 11 (c). Initially, we see that modal energy decreases faster with $\sigma=0.1$. However, higher level of modal energy is suppressed with $\sigma=1$ using higher input cost. This is also supported by our observation of the maximum pole movement in figure 9 (b). Once the modal amplitudes are forced to zero, the total modal energy reduces to zero. That is, the flow returns to the mean flow along with remnants of small structures corresponding to oscillator IV associated with increase in energy as shown in figure 11 (c). This is expected as the unsteady flow field is decomposed into mean flow and modal components, so that forcing the modal components to zero reduces the flow to the mean. However, the mean flow is not an equilibrium and a shift in the mean flow is observed as time progresses. This mean flow deformation is attributed to the Reynolds stress generated by the modal fluctuations which modifies the base flow (Brunton & Noack, 2015). The change in the base flow leads to a corresponding decrease in modal energy until equilibrium is achieved. In the case of the cylinder flow problem, the flow tends towards the unstable steady state.

As discussed earlier, the mean shift between the mean flow and the unstable steady state can be described by the shift mode. The temporal coefficient ($a_{\triangle}$) corresponding to the shift mode ($\boldsymbol{u}_{\triangle}$) can be obtained by projection as $a_{\triangle}=\left<\boldsymbol{u}-\bar{\boldsymbol{u}},\boldsymbol{u}_{\triangle}\right>$. We also perform control considering $\sigma=10$ and $100$. The variation of shift mode temporal coefficients with respect to temporal coefficients of the dominant POD modes for the range of $\sigma$ considered is shown in figure 12 (a). For reference, the evolution of the flow from unstable steady state to the mean flow that follows a parabolic inertial manifold is shown in grey (Noack et al., 2003). The darker grey region of the manifold indicates the region of effective control using the networked oscillator model. In this region, the nonlinear energy transfers are inhibited successfully to reduce wake oscillations. Thus, in suppressing the modal oscillations, we achieve a mean shift in the flow exhibited by the shift mode temporal coefficient.

We also present the temporal dynamics associated with the shift mode in figure 12 (b). It can be observed that initially $a_{\triangle}=0$ in the limit cycle and as control is applied, the shift mode temporal coefficient decreases and approaches the unstable steady state. As $\sigma$ decreases, we achieve a steeper decrease in $a_{\triangle}$. Examining the long time history of the temporal variation of the shift mode, the choice of lower $\sigma$ results in low frequency oscillations in $a_{\triangle}$. These low frequency oscillations are attributed to the exchange of energy between the low and moderate drag states (Munday & Taira, 2013). For control effort corresponding to $\sigma=100$, a steady state with application of control is achieved more rapidly.

We also evaluate the drag power $P_{D}$ and the invested actuation power for control $P_{a}$ given by

[TABLE]

where $\boldsymbol{f}$ is the actuation force added in DNS. The drag and actuation power required for $\sigma=0.1$ and $1$ are shown in figure 13 (a) and (b), respectively. As expected, the drag power for the controlled cases are lower than the baseline. We observe that the actuation power needed at steady state with control is significantly lower that the drag power. This results from suppression of the modal oscillations associated with the flow.

The networked oscillator control framework considered here is based on interactive dynamics of the baseline POD modes and does not include any other modes. We note that the shift mode itself can be incorporated into the controller (Noack et al., 2011; Luchtenburg et al., 2009). The shift mode is not added to the current formulation but such extension can be realized. To force the flow completely to the unstable steady state, additional feedback controller about the unstable steady state can be introduced. Details on such a patch controller design is given in appendix B.

We evaluate the unsteady forces on the cylinder with application of control compared to the baseline drag and lift coefficient variation in figures 12 (c) and (d), respectively. The variation in the drag coefficient is similar to the shift mode variation as discussed previously. Almost a $12\%$ reduction in drag and $70\%$ reduction in lift fluctuation are achieved with the application of control. We can thus reduce unsteady forces on the cylinder by suppressing modal oscillations in the flow and inhibiting energy transfers therein. Using a networked oscillator model in conjunction with optimal control, we can control energy transfer dynamics effectively for unsteady fluid flows.

4 Concluding Remarks

We construct a networked oscillator model for describing modal interactions in unsteady fluid flows. The modal oscillators comprised of POD conjugate mode pairs constitute the nodes on the modal network. The interactions between the oscillators form edges of the network, which are characterized by analyzing impulse-response to the fluid flow. Perturbations are introduced in the modal oscillators and their transfer of perturbation energy over the network is studied. The amplitude and phase perturbations introduced in the Navier–Stokes equations are tracked using regression techniques to develop the networked oscillator framework. The networked oscillator model is a linear approximation to the nonlinear modal interactions in unsteady fluid flows.

Using a canonical example of unsteady flow over a cylinder, the energy transfer dynamics are highlighted. Agreement of the model with DNS is observed for both amplitude and phase perturbation cases. The networked oscillator model captures the modal perturbation dynamics better compared to the empirical Galerkin formulation. An aggregate network model capturing general oscillator interactions in the flow is obtained using combined regression techniques over a collection of perturbed cases. The aggregate network model is attributed with the least error in terms of prediction of the temporal dynamics of the modal oscillators. The degree of the nodes corresponding to the aggregate network structure provides insights on the importance of individual oscillators to energy transfers and the overall system dynamics.

With the knowledge of network interactions between oscillators, an optimal feedback control strategy is designed to suppress oscillator fluctuations with respect to the natural limit cycle faster. A judicious choice of input oscillators forced is made by examining the movement of the graph Laplacian poles with LQR design. On control of modal perturbations in the flow, faster suppression of the oscillator temporal dynamics compared to the baseline limit cycle is observed. Controlling the overall fluctuations of the oscillators resulted not only in the suppression of the modal oscillations but also in a mean shift towards the unstable steady state exhibited by the shift mode. As energy transfers that sustain wake oscillations were inhibited, a mean shift of the flow to a lower energy state inside the parabolic inertial manifold is observed. The mean shift correspondingly led to a reduction in the unsteady forces on the cylinder. A $12\%$ reduction in drag force and a $70\%$ reduction in lift force are achieved with feedback flow control. In the long-time behavior of the controlled flow, a competition between low and moderate drag states are observed.

The networked oscillator modeling and control approach shown here has leveraged the knowledge of modal interactions providing insights beyond traditional approaches. Controlling the modal interactions at a fundamental level motivates analogous studies using localized actuators and limited sensors for modeling and controlling unsteady fluid flows.

Acknowledgements

We are grateful for the support by the U.S. Air Force Office of Scientific Research (Award number: FA9550-16-1-0650; Program Manager: Dr. Douglas R. Smith).

Appendix A

The performance of the aggregate model is shown in figures A1 and A2. DNS comparison of the oscillator dynamics using the aggregate network oscillator model with introduction of amplitude perturbations corresponding to $\beta_{m}=0.2$ on each oscillation is shown in figure A1. The comparison of predicted trajectories for phase perturbation of $\theta_{m}^{\prime}(t_{0})=-\pi/2$ on each oscillator is shown in figure A2.

Appendix B

To force the flow to the unstable steady state, we first perform a POD decomposition of the collected velocity snapshots $\boldsymbol{u}$ after the application of control near the unstable steady state $\boldsymbol{u}^{*}$ as

[TABLE]

In particular, here we consider the $\sigma=1$ case and collect velocity snapshots in the time interval between $t=19.3$ and $80$. During this time interval, we start to observe oscillations in the shift mode temporal coefficient. The extracted spatial modes and temporal coefficients are shown in figure B1 (a) and (b) respectively. Here, we show the dominant POD modes. We can see that mode $1$ shares similarities with the shift mode. Although modes $2$ and $3$ correspond to a frequency similar to oscillator I, the spatial modal structure is different. Here, we consider the first $10$ modes for control. We can construct a linearized model by curve-fitting the temporal coefficients. Designing an additional controller (patch controller) based on these coefficients, we can force the flow to the unstable steady state and reach minimum drag performance as shown in figure B1 (c) and (d).