Data-Driven Discovery of a New Ginzburg-Landau
Reduced-Order Model for Vortex Shedding

The simplest case of vortex shedding over a cylinder can be modeled by the incompressible, two-dimensional Navier-Stokes equations, together with mass continuity:

	$\displaystyle u_{x}+v_{y}$	$\displaystyle=0,$		(1a)
	$\displaystyle u_{t}+uu_{x}+vu_{y}$	$\displaystyle=-\frac{1}{\rho}p_{x}+\nu(u_{xx}+u_{yy}),$		(1b)
	$\displaystyle v_{t}+uv_{x}+vv_{y}$	$\displaystyle=-\frac{1}{\rho}p_{y}+\nu(v_{xx}+v_{yy}),$		(1c)

where $\rho$ is the fluid’s density, $\nu$ is the fluid’s kinematic viscosity, and the state variables are $(\textbf{u},p)=(u,v,p)$, the fluid’s horizontal and vertical velocities and pressure at every point in spacetime. At the fluid-body interface on the cylinder’s surface, the boundary conditions are the no-slip condition on the surface of the body, $\textbf{u}=0$ on $\Omega$. The nondimensional parameter Reynolds number, $\text{Re}={U_{\infty}L}/{\nu}$, characterizes the dynamics, where $U_{\infty}$ is the freestream velocity of the flow at the inlet before interacting with the body and $L$ is a characteristic length scale, in this case the diameter of the cylinder. The behavior of the system is then parameterized by Re and displays a wide variety of behaviors from creeping flow at low Re to fully turbulent flow at high Re [1]. Relevant for this work, for $\text{Re}\lesssim 47$, the solution $(u(x,y,t),v(x,y,t),p(x,y,t))$ is stable and steady, and presents as a long, narrow tail behind the cylinder that grows with increasing Re [2]. Once Re increases past a critical value $\text{Re}_{c}\approx 47$, the wake begins oscillating as discrete vortices begin to break free and are advected downstream [2]. Eventually the behavior saturates onto a limit cycle – the periodic vortex shedding. Figure 1 (left) shows the saturated flowfield at values of Re below, at, and above the critical threshold.

This behavior of transitioning from steady flow to periodic flow as Re increases past $\text{Re}_{c}$ is a hallmark of a supercritical Hopf bifurcation [63]. A Hopf bifurcation occurs in a system of differential equations when the real part of a pair of complex eigenvalues (the critical eigenvalues) switches from negative to positive as a critical parameter increases past its critical value [63]. The real part of the critical eigenvalues determines the linear growth rate, so a positive real part indicates that even infinitesimal perturbations will grow in amplitude and the steady state of the system is unstable [63]. In fluid flows, such perturbations may arise as small fluctuations in the velocity or pressure fields, for example due to physically imperfect boundaries or numerical round-off errors [2, 11, 12]. In a supercritical Hopf bifurcation, the growth of the perturbation is countered by nonlinear effects, which ultimately causes the amplitude of the perturbation to saturate onto a stable limit cycle.

The competing influence of the linear and nonlinear effects is captured by the Stuart-Landau equation (SLE) [64, 65],

$$\frac{d}{dt}A=\lambda A-\mu|A|^{2}A,$$

(2)

where $\lambda=\lambda_{r}+\mathbf{i}\lambda_{i}$ and $\mu=\mu_{r}+\mathbf{i}\mu_{i}$ are complex parameters. $A$ is known as the Landau variable, or order parameter [20], and in the vicinity of the bifurcation it represents the complex amplitude of the perturbation of the system away from the unstable equilibrium. Given that the perturbation starts with small amplitude, $A>|A|^{2}A$ for small enough $t$ and the linear term will dominate the early-time behavior of the solution. Positive $\lambda_{r}$ leads to exponential growth (i.e. linear instability) and $\lambda_{i}\neq 0$ causes oscillations. Once $A$ grows large enough, then the nonlinear terms begin to influence the solution in a significant way. When $\mu_{r}<0$, the nonlinear terms will counteract the exponential growth of the linear terms, resulting in a stable oscillatory solution. When $\mu_{r}$ is also positive, however, then the nonlinear terms further the growth of $A$. This is the subcritical case, where perturbations will grow to large amplitude and no nearby stable limit cycle exists [63]. Figure 1 (middle) shows the phase space of the SLE $(\mathbf{Re}(A),\mathbf{Im}(A))$ parameterized by the critical parameter $\lambda_{r}$, with $\mu_{r}<0$. The different flowfields (left) are linked to different behaviors in phase space depending on the value of the control parameter.

For a general nonlinear system, the SLE can be derived by expanding the dynamics for the amplitude of the critical eigenvectors (associated with the critical eigenvalues) in a Taylor series. Then to leading order, the supercritical Hopf bifurcation is modeled in the normal form by the SLE. The components of the complex amplitude are then the most critical spatial modes of the system, dominating the saturated oscillatory dynamics. These modes can be found, for example, via POD [27]. In the specific case of vortex shedding past a cylinder, much experimental, numerical, and analytic work has been done to show and justify the SLE as a ROM of vortex shedding [2, 11, 12, 27, 61]. Early investigations focused on the wake of a three-dimensional cylinder, using experimental fluid flow measurements [2, 11, 16, 66] and numerical simulations [12] to derive effective coefficients for the SLE governing the local oscillations at some point in the wake.

Noack et al. [27] developed a hierarchy of ROMs for the transient and steady cylinder wake using POD and Galerkin projection, first finding the dominant modes of the data, and then projecting the original data onto this lower-order subspace. The POD modes of flowfield data are the most critical spatial modes and represent the dominant oscillatory dynamics. Subsequent modes contribute progressively less energy to the overall dynamics, and higher modes capture nonlinear interactions and less dominant features [27, 30]. The models presented by Noack et al. capture both the saturated and transient behavior of vortex shedding. Their low-order model achieves this by incorporating a “shift” mode to account for the underlying unstable steady flow in addition to incorporating the critical modes capturing the dominant oscillatory behavior. Higher-order models include higher-order effects from the less-dominant modes, increasing the accuracy at the cost of additional computations and model complexity [27]. The critical POD modes for the flowfield from $\text{Re}=50$ can be found in Figure 1 (right). The POD modes, being coupled, display similar but offset periodic patterns. The work of Noack et al., along with the rest discussed in this section, establishes that the amplitudes of these modes can be effectively modeled by the SLE [11, 27, 61].

In [43], Barkley studies the time-averaged mean flow of the cylinder wake and undertakes a two-dimensional linear stability analysis of the base and mean flows for increasing Re. Contrasting the base and mean flows, he confirms via the leading eigenvalues of the base flow that the onset of oscillation corresponds to a Hopf bifurcation in which the growth rate crosses zero at $\text{Re}_{c}$, whereas the growth rate of the mean flow remains zero. Sipp and Lebedev [61] expand upon Barkley’s mean flow analysis by undertaking a global weakly nonlinear multiple scales analysis of the two-dimensional Navier-Stokes equations, developing a fast time coordinate $t_{1}=\epsilon t$. From this, they derive the SLE with the Landau variable as the slow time behavior of the critical global mode. Although the purely temporal SLE does not incorporate a spatial coordinate directly like the spatiotemporal CGLE does, their weakly nonlinear analysis is global in nature and provides the coefficients of the SLE governing the local oscillator at any point in space [61].

Although the SLE is a purely temporal ordinary differential equation, it was found that the measured coefficients varied depending on spatial location within the wake [2, 11], hinting at deeper spatial dependence [61]. It was later found that the values of the coefficients of the SLE vary depending both on measurements’ spatial locations within the wake and very strongly on end effects from finite cylinder length, and hence a more complete model was sought in the spatiotemporal one-dimensional complex Ginzburg-Landau equation (CGLE) [13, 14, 55, 67],

$$\frac{\partial}{\partial t}A=\lambda A-\mu|A|^{2}A+\kappa\frac{\partial}{\partial x}A+\gamma\frac{\partial^{2}}{\partial x^{2}}A,$$

(3)

where in general $A$ and all parameters are again complex. The CGLE is a partial differential equation (PDE), whose solution varies with both time and space. The CGLE may be thought of as an advective and diffusive coupling of an infinite field of Stuart-Landau oscillators. It is derived under the assumption of a long wavelength, supercritical Hopf instability in a homogeneous medium. However, the translational invariance and symmetry assumptions of the CGLE are not satisfied by the vortex shedding system [23], hence significant effort has been made in establishing a connection between vortex shedding and the CGLE.

In first proposing the CGLE as a ROM of vortex shedding past a cylinder, the one-dimensional CGLE was proposed with its spatial coordinate corresponding to the $z$-coordinate of the flowfield i.e. the spanwise coordinate of the cylinder; various boundary conditions were imposed to model different experimental setups [13, 14]. With coefficients derived from measurements of three-dimensional cylinder wake flows, simulations of the CGLE accurately recreated many of the patterns and regimes of behavior observed in three-dimensional cylinder flow, such as chevrons [62]. Some works also attempted to incorporate the streamwise- as well as the spanwise-coordinate into a two-dimensional CGLE for the three-dimensional cylinder wake, mainly in the form of spatially varying coefficients and under certain assumptions such as slow variations in the streamwise direction [56, 57]. These works expand on the applicability of the CGLE as a ROM for three-dimensional vortex shedding past a cylinder and further motivate the incorporation of streamwise dependence.

Analytically, via weakly nonlinear asymptotic expansion, Le Dizes et al. [68] deduced the one-dimensional CGLE from the Navier-Stokes equations as the equation governing the amplitude of the critical global mode, again in the spanwise direction. Noack [15] used a low-dimensional Galerkin method to analyze the stability of both steady and periodic cylinder wakes, finding that both the onset of periodicity in two dimensions at $\text{Re}_{c}$ and the onset of periodicity in three dimensions at $\text{Re}_{c,2}\approx 170$ are supercritical Hopf bifurcations modeled by the Landau equation. Zuccoli et al. [21] also used a multiple scales analysis on the Navier-Stokes equations for Taylor-Couette flow, using both slow time and space scales to derive a form of the one-dimensional CGLE, with the spatial coordinate as the spanwise coordinate and the Landau variable as the amplitude of a perturbation.

The discussion here and in the previous subsection firmly establishes that the Bénard-von Karman instability is characterized by the Hopf bifurcation, justifying the SLE and CGLE as ROMs of vortex shedding. From this, a hierarchy of ROMs for vortex shedding past a cylinder emerges, from an SLE description of the two-dimensional wake to local SLE and global CGLE descriptions of the three-dimensional wake. Thus the SLE has been rigorously linked to two-dimensional and three-dimensional vortex shedding, but the CGLE has only been linked to three-dimensional vortex shedding, mainly through the spanwise $z$-coordinate. Our contribution to this hierarchy of ROMs will be the data-driven discovery and development of a one-dimensional CGLE for the two-dimensional cylinder wake; by incorporating the streamwise $x$-coordinate into a novel Landau variable, we will answer whether streamwise dependence could be included in a CGLE ROM for vortex shedding.

More generally, ROMS capture the dominant features of the flow in a simpler mathematical framework. Modeling and solving the full Navier-Stokes equations, even with simplifications such as two-dimensionality and incompressibility as in Eqns. (1a)–(1c), is significantly more involved and expensive than modeling the SLE or CGLE. With reduced-order modeling, we capture the essential features of vortex shedding with a single complex variable rather than with a solution to the Navier-Stokes equations, and a better understanding of ROMs may lead to a better understanding of the original phenomenon

The CGLE is a complex partial differential equation. From Eqn. (3), the CGLE can be expanded into a system of two coupled, real differential equations; the two equations in the system describe the dynamics of the real and imaginary parts of the complex amplitude. We perform this derivation, first for clarity expanding the complex terms:

$$\frac{\partial(A_{r}+\mathbf{i}A_{i})}{\partial t}=(\lambda_{r}+\mathbf{i}\lambda_{i})(A_{r}+\mathbf{i}A_{i})-(\mu_{r}+\mathbf{i}\mu_{i})|A_{r}+\mathbf{i}A_{i}|^{2}(A_{r}+\mathbf{i}A_{i})+(\kappa+\mathbf{i}\kappa)\frac{\partial(A_{r}+\mathbf{i}A_{i})}{\partial x}+(\gamma_{r}+\mathbf{i}\gamma_{i})\frac{\partial^{2}(A_{r}+\mathbf{i}A_{i})}{\partial x^{2}}.$$

(4)

Then, expanding and separating the real and imaginary parts of the equation and relabeling the complex amplitude as $A=\omega+\mathbf{i}\eta$, we obtain:

	$\displaystyle A_{r}=\omega\;\;\Rightarrow\;\;\frac{\partial\omega}{\partial t}=\lambda_{r}\omega-\lambda_{i}\eta-(\omega^{2}+\eta^{2})(\mu_{r}\omega-\mu_{i}\eta)+\kappa_{r}\omega_{x}-\kappa_{i}\eta_{x}+\gamma_{r}\omega_{xx}-\gamma_{i}\eta_{xx},$		(5a)
	$\displaystyle A_{i}=\eta\;\;\Rightarrow\;\;\frac{\partial\eta}{\partial t}=\lambda_{i}\omega+\lambda_{r}\eta-(\omega^{2}+\eta^{2})(\mu_{i}\omega+\mu_{r}\eta)+\kappa_{i}\omega_{x}+\kappa_{r}\eta_{x}+\gamma_{i}\omega_{xx}+\gamma_{r}\eta_{xx}.$		(5b)

Finally, we expand the nonlinear terms:

	$\displaystyle\Rightarrow\;\;\frac{\partial\omega}{\partial t}=\lambda_{r}\omega-\lambda_{i}\eta-\mu_{r}\omega^{3}+\mu_{i}\omega^{2}\eta-\mu_{r}\omega\eta^{2}+\mu_{i}\eta^{3}+\kappa_{r}\omega_{x}-\kappa_{i}\eta_{x}+\gamma_{r}\omega_{xx}-\gamma_{i}\eta_{xx},$		(6a)
	$\displaystyle\Rightarrow\;\;\frac{\partial\eta}{\partial t}=\lambda_{i}\omega+\lambda_{r}\eta-\mu_{i}\omega^{3}-\mu_{r}\omega^{2}\eta-\mu_{i}\omega\eta^{2}-\mu_{r}\eta^{3}+\kappa_{i}\omega_{x}+\kappa_{r}\eta_{x}+\gamma_{i}\omega_{xx}+\gamma_{r}\eta_{xx}.$		(6b)

Notice the symmetry in the coefficients, with the real and imaginary parts of $\lambda$, $\mu$, $\kappa$ and $\gamma$ appearing on alternate terms between Eqns. (6a) and (6b) (e.g. $\gamma_{r}$ is on $\omega_{xx}$ for (6a), but on $\eta_{xx}$ for (6b)). This system of two real PDEs is the form of the CGLE we will derive from our data and numerically solve.

More generally, any generic two-dimensional system undergoing a Hopf bifurcation can be converted to its normal form, from which essential features of the dynamics can be determined and compared across systems. The conversion process involves applying a series of near-identity transformations, rescalings, and coefficient normalizations to the system in order to distill the normal form parameters [63]. The normal form of the Hopf bifurcation is:

$$\begin{bmatrix}\dot{y}_{1}\\ \dot{y}_{2}\end{bmatrix}=\begin{bmatrix}\beta&-1\\ 1&\beta\end{bmatrix}\begin{bmatrix}y_{1}\\ y_{2}\end{bmatrix}+\sigma^{\prime}(y_{1}^{2}+y_{2}^{2})\begin{bmatrix}y_{1}\\ y_{2}\end{bmatrix}+D_{x}(y_{1},y_{2}),$$

(7)

where $\beta\equiv{\zeta}/{\phi}$ and $\zeta$ and $\phi$ are respectively the real and imaginary parts of the complex eigenvalues. $D_{x}(y_{1},y_{2})$ represents the spatial derivative terms. $\sigma$ is the first Lyapunov coefficient and is computed as a function of second- and third-order derivatives of the nonlinear terms of the original system; reference Kuznetsov [63] for further details. The parameter $\beta$ determines whether a perturbation experiences growth ($\beta>0$) or decay ($\beta<0$), and the parameter $\sigma^{\prime}=\operatorname{sgn}(\sigma)=\pm 1$ determines whether the nonlinearity will saturate a perturbation’s growth ($\sigma^{\prime}=-1$) or amplify it ($\sigma^{\prime}=+1$). The case where a perturbation tends to grow but its growth is saturated by the nonlinear terms has parameters ($\beta>0$ , $\sigma<0$) and is the supercritical Hopf bifurcation. This bifurcation results in a stable limit cycle, while in the subcritical case, an unstable limit cycle exists for $\beta<0$. By computing the parameters ($\beta$, $\sigma$), we can distill and compare basic features across different systems; in § 4 of this work, we will compute ($\beta$, $\sigma$) for a variety of systems generated with SINDy in order to compare and discuss them.

The data is obtained through a computational fluid dynamics (CFD) simulation using the immersed boundary projection method (IBPM) [69, 70]. IBPM uses 3rd order Runge-Kutta integration to solve the vorticity-stream function formulation of the two-dimensional Navier-Stokes equations, simulated here in a domain of $x\in[-2,52]$ and $y\in[-8,8]$ with a cylinder of diameter $D=1$ situated at the origin. We used a grid of $dx=dy=0.02$ resulting in 2.16 million mesh points, and the time step was $dt=0.02$ with an initial condition of zero velocity and pressure everywhere. We impose no-slip boundary conditions on the surface of the body, far-field freestream conditions on the vertical domain boundaries, and inflow and outflow conditions on the left and right horizontal domain boundaries, respectively.

A nondimensional measure of the frequency of the system is the Strouhal number, $\text{St}=fL/U_{\infty}$, where $f$ is the frequency of vortex shedding. St can be measured from the coefficient of lift $C_{l}(t)$ (discussed more later), and for $\text{Re}=50$, we obtain $\text{St}=0.126$. This is in-line with Williamson [72] and Park et al. [73] who give 0.12 and 0.124, respectively. The specific flowfield data used throughout this work in order to learn a sparse, interpretable model is vorticity $\omega(x,y,t)$ at $\text{Re}=50$, specifically a single period of flowfield data after saturation onto the limit cycle.

We seek to transform our flowfield data into a form that can be modeled as a one-dimensional Ginzburg-Landau system depending on $(x,t)$, hence we must find a way to transform our two-dimensional data $\omega(x,y,t)$ into one-dimensional data $\overline{\omega}(x,t)$. There are many reasonable ways to do this. For example, taking the vorticity along a single horizontal line at a specified $y-$coordinate is simple, but neglects to incorporate the vertical distribution of vorticity in any way. Instead, we choose to integrate along the $y$-coordinate for each $(x,t)$,

$$\overline{\omega}(x,t)=\int_{-H}^{+H}\omega(x,y,t)dy,$$

(8)

which incorporates the vertical distribution of vorticity while still producing one-dimensional data in a simple manner. The vortices spread out vertically as they dissipate downstream, hence the ideal bounds of integration would be $-\infty\leq y\leq\infty$; in practice the vertical movement of a vortex is slow relative to its horizontal movement, and a vertical limit of $H=8D$ at $x=50D$ downstream is more than sufficient. We perform this integration at every point $(x,t)$, effectively removing the $y$-coordinate and thus reducing the dimension of our data from $(x,y,t)$ to $(x,t)$. Figure 2a depicts a snapshot of the vorticity flowfield data from the CFD simulation overlaid with the bounds of integration for a particular $x_{i}$; Figure 2b shows the instantaneous coarse-grained vorticity $\overline{\omega}(x;t=0)$. Figures 2c and 2d repeat this process through time to produce the full spatiotemporal coarse-grained vorticity $\overline{\omega}(x,t)$.

The coarse-grained system captures the spatiotemporal regularity and pattern-forming behavior of the original vortex shedding flow. Each peak or trough in the periodic signal respectively corresponds to a vortex of positive or negative vorticity in the original flowfield. The signal is zero upstream of the cylinder where the flow is unperturbed, then begins to oscillate sharply in the region of the cylinder $x\in[-0.5,0.5]$. As the vortices detach and shed downstream, the coarse-grained signal smoothly propagates, its amplitude growing and peaking in the mid-wake region, after which it slowly decays. Since vorticity is diffusive and the signal is clearly decaying, in a semi-infinite domain $x\in[0,+\infty)$ we expect that $\overline{\omega}(x,t)\to 0$ as $x\to\infty$.

We propose this method of coarse-graining for two reasons: first, it is a simple method to reduce the spatiotemporal dimension of the data while incorporating information from the original flowfield; and second, it is in analogy to an analysis of the time-averaged mean flow, as in Barkley [43]. In this work, we propose a spatial analog of mean flow analysis that leads to novel, interpretable results.

We seek to discover the real and imaginary components of the CGLE in the 2-equation form of Eqns. (6a)–(6b), hence we require a second signal in addition to our coarse-grained vorticity $\overline{\omega}(x,t)$. As with coarse-graining, there are many reasonable ways to lift the dimension of a system, in this case from 1 signal to 2 signals. For example, Loiseau et al. based their system off of the coefficient of lift $C_{l}(t)=2L/(\rho U_{\infty}^{2}L)$, a nondimensional measure of the lift force on the body [38]. Reference Figure 3 (left) for a plot of $C_{l}(t)$ from our data at $\text{Re}=50$. Loiseau et al. then lifted the dimension of their system by differentiating $C_{l}(t)$ with respect to time; differentiating the signal $C_{l}(t)$ produces a second periodic signal of the same frequency and similar (but not the same) amplitude but with a phase shift. Loiseau went on to use SINDy on these signals to discover a Stuart-Landau-like model for the coefficient of lift. Our discovery of a Ginzburg-Landau-like model from spatiotemporal data is a natural extension of that work; here however, we produce this second signal via time-delay embedding, which also produces a second periodic signal of the same frequency and same amplitude, and which is less prone to amplifying numerical errors than derivative embedding.

Specifically, if our signal is $C_{l}(t)$, then the time-delay embedded system is $\mathcal{U}=[C_{l}(t),C_{l}(t+\tau)]$, where $\tau$ is the degree of time-delay embedding [45, 46]. Choosing a good $\tau$ is crucial to recovering a good model, but it’s not obvious a priori what that choice should be. Figure 3 shows the phase space of the system under different time-delay embeddings and under time-differentiation.

Given a period of oscillation $\mathcal{T}$ of $C_{l}(t)$, the orbit in phase space of the time-delay embedded system is most circular for $\tau=\frac{1}{4}\mathcal{T}$, hence maximizing the independence of the signals. As the time-delay embedding increases or decreases, the signals become increasingly in-phase until a trivial time-delay embedding of $\tau=\frac{1}{2}n\mathcal{T},n=0,1,2...$ and we recover perfectly colinear features. If instead of time-delay embedding, we pair $C_{l}(t)$ with its time derivative, as can be seen in Figure 3(d), the effect is to produce an unrotated ellipsoidal orbit in phase space, resulting in a similar effect of obtaining mostly independent features. However, finding the numerical derivative of $C_{l}(t)$ will introduce and amplify numerical errors. Recalling Figure 1, the real and imaginary parts of the Stuart-Landau equation also produce a circular orbit in phase space.

Returning to our spatiotemporal coarse-grained vorticity data $\overline{\omega}(x,t)$, we make an initial choice of $\tau=\frac{1}{4}\mathcal{T}$ as this produced maximally independent features in the related system $C_{l}(t)$. Doing so produces the system of time-delay embedded coarse-grained data $\mathcal{U}=[\overline{\omega}(x,t)\,,\,\eta(x,t)]$ where $\eta(x,t)\equiv\overline{\omega}(x,t+\tau)$ is the time-delay embedded signal. Figure 4a shows line plots of the instantaneous system of coarse-grained data $\mathcal{U}$. We seek a model for the dynamics $\mathcal{U}_{t}$ of this system, and the instantaneous dynamics are shown in Figure 4b. Time-differentiating the spatiotemporal coarse-grained vorticity produces higher-order behavior evident at the peaks of the oscillations in $\mathcal{U}_{t}$ at $10\lesssim x\lesssim 30$, in which the curve doubles back on itself. This is also the area of the global maximum of $\overline{\omega}$ and $\eta$, occurring at $x\approx 20$. Higher-order behavior is less easily explained by a sparse system of PDEs, hence the inclusion of a time-differentiated signal in the base system of signals $\mathcal{U}$ would produce more complex dynamics. By time-delay embedding our data, we avoid introducing additional complexity into the signals and dynamics. Figures 4c and 4d show the full spatiotemporal plots of the coarse-grained vorticity $\overline{\omega}(x,t)$ and its time derivative $\frac{\partial}{\partial t}\overline{\omega}(x,t)$. Note that the color plots of the time-delay embedded data (not shown) would look identical except for a phase shift of $\frac{1}{4}\mathcal{T}$.

As a final note, our system of data $\mathcal{U}=[\overline{\omega},\eta]$ was a result of specific choices: what CFD flowfield data to use, how to coarse-grain it, and how to produce a second commensurate signal. There are other choices that may also produce a system of signals that could be used as a basis for developing a ROM for cylinder vortex shedding. For example, one might directly incorporate the time derivative of the data by using a system like $\mathcal{U}=[\overline{\omega},\overline{\omega}_{t}]$, though with the drawback that numerical differentiation will introduce and amplify numerical errors. Alternately, one might change the underlying flowfield data by coarse-graining the horizontal $u(x,y,t)$ and vertical $v(x,y,t)$ components of the velocity flowfield, $\mathcal{U}=(\overline{u},\overline{v})$.

Having obtained our flowfield data $\omega(x,y,t)$ and post-processed it into our system of coarse-grained, time-delay embedded signals $\mathcal{U}$, we now develop a sparse, interpretable system of PDEs to describe the dynamics of our spatiotemporal data using SINDy [36, 37, 71]. Given a library of candidate functions, SINDy leverages sparsity-promoting and regularizing techniques to determine which subset of candidate terms explains most of the dynamics: the reconstruction error is minimized under some norm while a sparsity-promoting regularization term ensures a final model which is sparse and hence interpretable. In this work, we specifically use an extension of SINDy known as PDE-FIND [37], which extends the SINDy framework from ODEs to PDEs.

Given data $\mathcal{U}$, we propose to model it as a dynamical system

$$\mathcal{U}_{t}=\mathcal{N}(f_{1}(\mathcal{U}),...,f_{i}(\mathcal{U}),...),$$

(9)

where the library terms $f_{i}$ are a set of functions which we apply to the original data, for example powers, spatial derivatives, and products thereof. For a homogeneous system, the coefficients are assumed constant. We can frame this as a matrix-vector product:

$$\mathcal{U}_{t}=\Theta\>\cdot\>\xi,$$

(10)

where $\Theta$ is the library of candidate functions $f_{i}(\mathcal{U}$) and $\xi$ is the vector of coefficients for each function. Given some data $\mathcal{U}$ we wish to model as a dynamic system as in (9), we find the numerical time-derivative of the data $\mathcal{U}_{t}$ and evaluate the candidate functions $f_{i}$ within $\Theta$. The choice of candidate function library $\Theta$ is significant. In theory, one may fill the library with as many functions as desired, but in practice, prior knowledge of the system and dynamics helps to select the most appropriate functions. Previous work has established a clear relationship between vortex shedding and the Stuart- and Ginzburg-Landau equations [13, 68, 11, 66]. Given this, we choose a library of functions inspired by the form of the CGLE, including polynomial terms up to third order, spatial derivatives up to second order, and variations thereof; the exact libraries of functions used throughout this study can be found in Appendix A. With many possible candidate functions, the regression problem for $\xi$ is underconstrained. We desire that $\xi$ be sparse (i.e. it has many zero entries), since having fewer terms with nonzero values results in a more interpretable dynamical system for $\mathcal{U}_{t}$. This is achieved by solving an optimization problem,

$$\hat{\xi}=\text{argmin}||\Theta\xi-\mathcal{U}_{t}||_{2}^{2}+\lambda||\xi||_{0},$$

(11)

where $\hat{\xi}$ is the sparse version of $\xi$ returned from SINDy. In practice, the use of the $l_{0}$ norm to regularize $\xi$ results in an $np$-hard optimization problem, so we solve the optimization problem with a convex relaxation, such as the sequentially thresholded least-squares (STLS) regression algorithm [36].

We can then compute the matrix-vector product $\hat{\mathcal{U}}_{t}=\Theta\cdot\hat{\xi}$, obtaining SINDy’s reconstruction of the time-derivative of the data. We then calculate the error between the data and the SINDy reconstruction:

$$\epsilon=\mathcal{U}_{t}-\hat{\mathcal{U}}_{t}=\mathcal{U}_{t}-\Theta\cdot\hat{\xi}.$$

(12)

We then find the normalized root-mean-square error (NRMSE), defined as:

$$\text{NRMSE}=\sqrt{\frac{\displaystyle\sum_{t}\sum_{x}\epsilon^{2}}{\displaystyle\sum_{t}\sum_{x}\mathcal{U}_{t}^{2}}}.$$

(13)

The NRMSE is computed for each signal within $\mathcal{U}$, and the total NRMSE is the $l^{2}$ norm of the individual NRMSEs. This single error measurement for the entire system is a metric of the model’s quality. After obtaining a PDE with low NRMSE, the final step in this analysis is to numerically solve the system of PDEs and compare the result to the original data, observing the quality and similarities of the models.

To begin our analysis, we use SINDy to fit a Ginzburg-Landau model to the coarse-grained vorticity. We first find a model for the local dynamics of the wake near the cylinder, $x\in[2,9]$. Then we find a global model for a much larger wake of $x\in[2,50]$. For both spatial domains, Figure 5 shows the input dynamics $\overline{\omega}_{t}$, the SINDy fit $\hat{\overline{\omega}}_{t}$, and their error $\epsilon$; as before, the accompanying signal $\eta_{t}$ (not shown) is similar to $\overline{\omega}_{t}$ but with a phase shift of $\frac{1}{4}\mathcal{T}$. The identified coefficients of the CGLE as well as the NRMSE and normal form parameters $(\sigma,\beta)$ are shown in Table 1. $\overline{\omega}_{t}$ and $\hat{\overline{\omega}}_{t}$ overlap on the line plots, with the error at any given point being at most one order of magnitude smaller than $\overline{\omega}_{t}$ and $\hat{\overline{\omega}}_{t}$ at that point. The NRMSE is NRMSE = $6.7\times 10^{-3}$ for the near-wake and NRMSE = $3.2\times 10^{-2}$ for the entire wake. SINDy is clearly able to reconstruct the dynamics in the near-wake and the entire wake with little error.

The spatial behavior of the coarse-grained dynamics and their time derivative have been discussed in § 3; like $\overline{\omega}_{t}$, $\epsilon$ also varies with $x$, but whereas $\overline{\omega}_{t}$ attains its global maximum at $x\approx 20$ and decays from there, $\epsilon$ has a global maximum in the mid-wake ($x\approx 10$) and a secondary local maximum in the far-wake ($x\approx 30$); these are areas in which a greater amount of the dynamics are not explained by the terms in $\Theta(\overline{\omega},\eta)$. These areas respectively coincide with the beginning and end of the higher-order behavior in $\overline{\omega}_{t}$ and $\hat{\overline{\omega}}_{t}$, and $\epsilon$ reaches a local minimum at $x\approx 20$, in between its primary global maximum and the secondary local maximum. This coincides with the global maximum of $\overline{\omega}_{t}$ and the peak of its higher-order behavior. This indicates that the two areas of transition between lower- and higher-order behaviors and back are least explained by the modeled CGLE. This supports the notion that the wake may be separated into three regimes of behavior: the near-wake, the mid-wake, and the far-wake. These subdomains will be explored more in the following subsections, where we both vary the subdomain of the data and tune the library terms and sparsity used to generate the model.

For both spatial domains, SINDy naturally identified coefficients whose values and signs show the same symmetry pattern expected from Eqns. (6a)–(6b), for example identical coefficients ($\lambda_{r}$) between the $\overline{\omega}$ term in the equation for $\overline{\omega}_{t}$ and the $\eta$ term in the equation for $\eta_{t}$. There is one key difference however: in a given equation, all four cubic terms have different coefficients, instead of two terms with coefficient $\mu_{r}$ and two with $\mu_{i}$. However, the four terms in one equation are still matched with four opposing terms in the other equation, hence instead of a single $\mu=\mu_{r}+\mathbf{i}\mu_{i}$ suggested by Eqns. (6a)–(6b), we instead have $\mu_{1}=\mu_{1_{r}}+\mathbf{i}\mu_{1_{i}}$ and $\mu_{2}=\mu_{2_{r}}+\mathbf{i}\mu_{2_{i}}$. While this means the generated equation cannot be neatly factored back into the form of Eqn. (3), we may still compute the normal form parameters $(\sigma,\beta)$. The normal form characterizes the differences in behaviors of different systems arising from differences in signs and magnitudes of some terms between the different spatial domains (e.g. $\lambda_{r}$ and $\lambda_{i}$); variations across the terms give rise to different system behaviors. For example, $(\sigma<0,\beta>0)$ for $x\in[2,9]$ corresponds to the classic supercritical Hopf bifurcation with an unstable fixed point and a stable limit cycle, consistent with previous conclusions that vortex shedding is characterized by a supercritical Hopf bifurcation [11, 12, 66]. Conversely, $(\sigma>0,\beta<0)$ for $x\in[2,50]$ corresponds to a subcritical Hopf bifurcation with a stable fixed point and unstable limit cycle. Differences in the local and global ROMs will be explore further in § 4.2.

Having found a sparse system of PDEs that closely fits our data with minimal reconstruction error, we use a pseudospectral integration method to compute a numerical solution for the SINDy model generated from the near-wake. We focus on this model for now because of its lower error; local models for different subregions of the wake and global models for the entire wake will be explored more later. Figure 6 compares the original data $\mathcal{U}$ in the top row to the numerical solution of the generated system of PDEs in the bottom row. The numerical solution shows a close agreement with the data in the near-wake, reinforcing the validity of the CGLE as a model for our data and as a ROM for vortex shedding.

The numerical solution of our generated model closely matches the original data and lends great insight into its behavior, but some of its features are an approximation. Although convenient and simple to implement, a pseudospectral method supposes periodic boundaries which our original data clearly does not have. Our data originally comes from the vortex shedding wake behind a cylinder, and once coarse-grained, there is clearly strongly time-dependent behavior at the boundaries. In particular, the behavior of $\overline{\omega}$ near the surface of the cylinder is nearly discontinuous. This poses a challenge, because the model found by SINDy is generated without consideration of well-posedness and boundary conditions.

Further, periodic boundary conditions allow our initial pulse to eventually self-interact once it has propagated the length of the domain. To overcome this, we embed the solution in a large spatial domain and investigate a spatiotemporal subdomain comparable to the original data, early enough in time that the initial pulse has not yet significantly interacted with itself. Although this strategy works well, it works better when we also include regularization in the form of hyperdiffusion terms (fourth-order spatial derivatives e.g. $\eta_{xxxx}$) in each equation with coefficient $R$. Hyperdiffusion serves to damp high frequency oscillations that result from self-interaction in a periodic domain. With distinct $\mu_{1}$ and $\mu_{2}$ and the inclusion of hyperdiffusion, our final equations then become the following:

	$\displaystyle\frac{\partial\overline{\omega}}{\partial t}=\lambda_{r}\overline{\omega}-\lambda_{i}\eta-\mu_{1,r}\overline{\omega}^{3}+\mu_{2,i}\overline{\omega}^{2}\eta-\mu_{2,r}\overline{\omega}\eta^{2}+\mu_{1,i}\eta^{3}+\kappa_{r}\overline{\omega}_{x}-\kappa_{i}\eta_{x}+\gamma_{r}\overline{\omega}_{xx}-\gamma_{i}\eta_{xx}+R\overline{\omega}_{xxxx},$		(14a)
	$\displaystyle\frac{\partial\eta}{\partial t}=\lambda_{i}\overline{\omega}+\lambda_{r}\eta-\mu_{1,i}\overline{\omega}^{3}-\mu_{2,r}\overline{\omega}^{2}\eta-\mu_{2,i}\overline{\omega}\eta^{2}-\mu_{1,r}\eta^{3}+\kappa_{i}\overline{\omega}_{x}+\kappa_{r}\eta_{x}+\gamma_{i}\overline{\omega}_{xx}+\gamma_{r}\eta_{xx}+R\eta_{xxxx}.$		(14b)

While we use a large spatial domain in order to prevent self-interaction for as long as possible, too large of a domain requires infeasibly long calculation time. Hyperdiffusion with sufficiently large $R$ serves to damp the higher-order oscillations that begin once the signal begins self-interacting; however, we generally want a value of $R$ as low as possible as it will damp the main oscillations as well. For the numerical solution in the bottom row of Figure 6, we use $R=10^{-3}$. This value of $R$ was found to be optimal, with larger values of $R$ resulting in the same signal but with lower amplitude, and smaller values of $R$ resulting in quickly unstable solutions. Regularization in the form of hyperdiffusion is easy to include in the framework of a pseudospectral numerical solution and is highly effective.

The analysis of this section has found close agreement between the original data, the model generated by SINDy, and the numerical solution of that model. As such, our method of coarse-graining, time-delay embedding, and model discovery with SINDy has resulted in a novel one-dimensional CGLE for two-dimensional cylinder vortex shedding, incorporating for the first time the streamwise coordinate into the Landau variable. In the following subsection, we will explore the spatial variation of local models found at different points in the wake.

In the previous subsection, we found a sparse model in the form of the CGLE to explain the dynamics of the coarse-grained vorticity in both the near-wake and the entire wake. Now we refine our method by varying and tuning certain parameters of our data-driven method. First we vary the spatial subdomain of interest and show how the error of the fitted model varies as the subdomain moves downstream. Then we expand the library of terms used to generate the model in order to better capture higher-order behavior in our data. We discuss each analysis in turn.

Regarding the spatial domain of the data, the initial work of § 4.1 investigated both a specific flowfield subdomain of $x\in[2,9]$ inspired by Rudy et al. [37], and a much larger flowfield of $x\in[2,50]$ to fully investigate the global wake dynamics. From this, it was clear that the data, the model coefficients, and the error are spatially dependent, so we next investigate how generated CGLE models vary depending on the location within the wake they are modeling. Specifically, for a subdomain of length $W_{L}$, SINDy generates a model from a specified library. The spatial subdomain has a midpoint $x_{mid}$ that then increments downstream, and the SINDy fit process is repeated. At each station within the wake, the NRMSE of the fitted model and its normal form parameters are computed. This entire process is then repeated for different values of $W_{L}$. In the limit $W_{L}\rightarrow 0$, the spatial coordinate of the coarse-grained data vanishes and our spatiotemporal data $\mathcal{U}(x,t)$ collapses to purely temporal data $\mathcal{U}(t;x_{mid})$ at that point $x_{mid}$. For this limiting case, we use SINDy to fit an ordinary differential equation in the form of the SLE (Eq. (2)), which mirrors the CGLE (Eqs. (6a)–(6b)) but without the spatial derivative terms. This is “Library SL” and the full list of terms is shown in Appendix A.

Regarding the library of model terms, the initial work applied SINDy to the data with Library GL, a judiciously chosen library of candidate functions $\Theta$ with terms exactly matching those in the CGLE. Now we will allow for additional nonlinear terms in three forms: higher-order polynomial terms up to fifth order, products of the polynomial and derivative terms (e.g. $\eta^{2}\eta_{xx}$), and combinations of those. These libraries are named for the total number of terms they feature, and respectively are referred to as “GL₂₄”, “GL₄₉”, and “GL₁₀₄”. For reference, the SL and GL Libraries have 9 and 13 terms. The full list of terms for all libraries are shown in Appendix A. These libraries contain many terms, and some of the larger libraries will have terms that contribute very little to the observed dynamics. Later we will tune the sparsity parameter $\lambda$ from Eq. (11) in order to eliminate some terms from these larger libraries, leaving behind only those additional nonlinear terms that better explain the dynamics of the coarse-grained vorticity. The best sparse models will explain as much of the dynamics with as few terms as possible. For now, the NRMSE of the nonsparse systems will provide a baseline to which we may later compare the sparse models.

As before, we quantify the goodness of fit of the model returned by SINDy with the NRMSE. Figure 7 shows how the local models’ NRMSE varies spatially in the wake, with parameterization by the spatial subdomain size and by the choice of library. Figure 7(a) shows how the NRMSE from models generated with our initial GL Library varies with spatial location in the wake, plotted specifically against $x_{mid}$, the midpoint of the spatial domain. In dark blue is the base case using (GL,7). In a lighter shade of blue and with different linestyles are the NRMSE for cases with larger and smaller $W_{L}$, (GL,15) and (GL,2). In red is the case (SL,0), which does not incorporate any spatial derivative terms. Figure 7(b) shows the spatial variation of NRMSE from models generated with equal subdomain sizes $W_{L}=7$ and increasingly larger libraries “GL₂₄”, “GL₄₉”, and “GL₁₀₄”. The dark blue case of (GL,7) is identical on both plots of 7.

In both Figures 7(a)-(b), the NRMSE shows clear spatial dependence in for all $W_{L}$ and model library. The error is moderate in the near-wake, then grows to a global maximum at $x\approx 20$ in the mid-wake before a sharp decline in the far-wake. Indeed, the NRMSE is 1-2 orders of magnitude lower in the far-wake than at its maximum in the mid-wake. In the mid-wake, the dynamics are characterized by the nonlinear growth and saturation of the SLE and CGLE and potentially by higher-order effects, but in the far-wake, the dynamics are primarily dominated by diffusion. Since the NRMSE is lower where diffusive dynamics dominate, this behavior is better explained by the terms in GL Library. Furthermore, in Figure 7(a), the NRMSE generally decreases as $W_{L}$ decreases but variation with $x_{mid}$ remains similar, suggesting more accurate models can be fit to a more narrow spatial slice of the data. However, the limiting case of $W_{L}=0$, which necessarily has no spatial component and hence no spatial derivatives in its model term library, has larger NRMSE in the far-wake region than does the base case of (GL,$W_{L}=7$), reinforcing that derivative and specifically diffusive terms dominate the far-wake coarse-grained dynamics.

In Figure 7(b), the NRMSE of all local models improves as the library of candidate functions $\Theta$ becomes larger. Interestingly, the peak at $x\approx 20$ is much less pronounced for Library GL₁₀₄ (the largest library), though the error in the far-wake is still an order of magnitude lower. This suggests that the mid-wake region of the flowfield is more influenced by higher-order effects, hence when higher-order terms are included in the model, the error at the mid-wake greatly decreases. In the far-wake, where the dynamics are dominated by diffusion, the inclusion of additional higher-order terms helps the least. The minimum NRMSE achieved anywhere across this parameter and hyperparameter tuning is $\approx 2\times 10^{-5}$, occurring far downstream for a model generated from Library GL₁₀₄. This is over 2 magnitudes of an improvement in NRMSE over the initial case (Table 1). However, this has come at the cost of interpretability, as with zero sparsity ($\lambda=0$) this model has over 100 terms and cannot feasibly be numerically simulated; in Appendix B, we will tune our sparsity parameter $\lambda$ to remove those terms which have the least effect on the error.

The local models clearly vary spatially, hence by investigating the coefficients of the models in each subdomain and computing their normal form parameters $(\sigma,\beta)$, we can draw conclusions about the behavior of that local model and the spatial characteristics of the coarse-grained data. Figure 8 shows $\sigma$ and $\beta$ vs. $x_{mid}$ for both the spatiotemporal ($W_{L}\neq 0$) CGLE models from Library GL and for purely temporal ($W_{L}=0$) SLE models from Library SL. The normal form parameters were also computed for models (GL,15) and (GL,2) from Figure 7(a), but were found to deviate only marginally from the base case (GL,7), and hence are not shown. This suggests that $W_{L}$ does not have a strong effect on the normal form parameters of the local model, at least not until it vanishes for (SL,0).

The behavior of a Hopf bifurcation depends on the signs of $(\sigma,\beta)$; specifically, the sign of $\beta$ determines whether the fixed point is stable ($\beta<0$) or unstable ($\beta<0$), and the sign of $\sigma$ determines whether the nonlinear terms cause a stable limit cycle ($\sigma<0$) or unstable limit cycle ($\sigma>0$). As the spatial domain moves downstream, the local model transitions between several distinct regions of behavior, delineated in Figures 8(a)-(b). For each of these regions, Table 2 shows the normal form parameters $(\sigma,\beta)$, the stability of the equilibrium at the origin, the existence and stability of the limit cycle, and figures of the normal form phase portraits.

Near the cylinder for $0\lesssim x\lesssim 5$, Region I has normal form parameters ($\sigma<0$, $\beta>0$). Moving downstream, $\sigma$ and $\beta$ undergo a variety of sign and magnitude changes, transitioning through three unique sets of ($\sigma$, $\beta$). Throughout the wake, the normal form parameters from the GL Library and SL Library fits and generally agree with each other. While derivative terms are excluded from the SL Library, the nature of the data clearly suggests advective and diffusive behavior. In the absence of derivative terms, the values of the remaining coefficients compensate, and their changing values produce different $(\sigma,\beta)$. The main distinction is the differing behavior in the far wake, namely the $(SL,0)$ fits never experience Region III and instead experience a brief return of Region I at $x\approx 30$. Both regions feature $\beta>0$ and are distinguished by whether $\sigma<0$ (Region I) or $\sigma>0$ (Region III). For the (GL,7) model, $\beta>0$ for $18\lesssim x\lesssim 25$ while $\sigma$ remains positive but decreasing to an asymptotic value. For the $(SL,0)$ model however, $\sigma$ becomes negative and $\beta$ simultaneously becomes positive, then both quickly switch back. Although difficult to determine visually, close inspection reveals that the bounds of these regions for both parameters coincide exactly, and both parameters are $\mathcal{O}(10^{-4})$. Recalling Figure 7, the NRMSE produced by the model $(SL,0)$ sees a pause in its decline at $x\approx 30$ while other models’ NRMSEs continue to decrease. This indicates that the fit from $(SL,0)$ in this region is relatively worse and the dynamics here may be due to higher-order effects not captured by the CGLE and not represented by its normal form parameters.

Further regarding the far-wake behavior, there is suggestive behavior indicative of long-range spatial variation in both parameters. While $\sigma$ reaches an asymptotic value in the far-wake for Library GL fits, $\sigma$ from Library SL briefly becomes negative in the far-wake before increasing and possibly attaining a secondary local maxima at $x\approx 50$. Regarding $\beta$, for all three fits the parameter decreases moving downstream after its maximum, but appears to begin increasing again far downstream at $x\approx 45$. Although the present spatial domain of analysis is insufficient to clearly capture the behavior, it appears that the local models’ evolve on a large length scale and the normal forms have long-range spatial variation.

The spatial variation of normal forms may also suggest a connection with the wavemaker region, a region of local absolute instability known to exist in the wake close to the cylinder, $0\lesssim x\lesssim 4$ [19, 58]. In this region, perturbations grow and ultimately affect the entire domain. This near-wake region corresponds to Region I in Figure 8, where $\beta>0$ indicating that the fixed point is unstable and perturbations to this state grow. Sufficiently far downstream, the local flowfield is convectively but not absolutely unstable, and correspondingly, $\beta<0$, signaling that perturbations from the stable fixed point decay, consistent with a flowfield in which the vorticity slowly decays away far downstream. The recurrence of regions of $\beta>0$ amid the downstream diffusion may again indicate longer-range behavior or a secondary wavemaker region.

The analysis of the previous section explored generated local models’ spatial behaviors and variations. In this section, we further iterate on the form of the model and minimize the NRMSE by fitting a global heterogeneous CGLE to the data of the entire wake $x\in[2,50]$. This is done by introducing $x-$dependent coefficients to replace constant coefficients on the model terms. By explicitly allowing for spatial variation in the coefficients of the terms, we find a global model for the entire wake that captures the spatially changing dynamics. The general form of our dynamical system becomes:

$$\mathcal{U}_{t}=\mathcal{N}(g_{1}(x)\cdot f_{1}(\mathcal{U}),...,g_{i}(x)\cdot f_{i}(\mathcal{U}),...),$$

(15)

where our data $\mathcal{U}=[\overline{\omega},\eta]$ is again our coarse-grained vorticity and its time-delay embedding. First, we will try using $g_{i}(x)$ as polynomials in $x$ up to a specified order as an approximation to the Taylor series of the true, unknown coefficient functions for each term in the Library GL. In Figure 9(a), the NRMSE of a generated model is plotted against the order $n$ of the Taylor polynomials. The different libraries are:

• •

  H-GL₁ – Taylor series in $x$ on the linear terms of the GL library

• •

  H-GL₃ – Taylor series in $x$ on the cubic terms of the GL library

• •

  H-GL_1,3 – Taylor series in $x$ on both the linear and cubic terms of the GL library

At $n=0$, the coefficients are constant and we reproduce the initial case with Library GL on the global wake. The NRMSE then improves with Taylor polynomials on the linear terms, more with Taylor polynomials on the cubic terms, and even more with Taylor polynomials on both the linear and cubic terms. The NRMSE generally improves with higher order $n$ of the Taylor polynomials, although improvements generally diminish around $n=7$, after which the NRMSE of the generated model experiences only marginal improvements. At $n=9$, Library H-GL_1,3 achieves a minimum NRMSE of $1.2\times 10^{-3}$, an order of magnitude improvement over the initial case (NRMSE = $3.2\times 10^{-2}$). The heterogeneous system’s improvement of the global model’s NRMSE further reinforces that the vortex shedding wake and its ROMs have spatially-dependent behavior that a homogeneous global model fails to capture.