Lakshmi N. Sridhar
Chemical Engineering Department, University of Puerto Rico, Mayaguez, PR, USA.
DOI: 10.4236/ajcm.2026.163011   PDF    HTML   XML   2 Downloads   20 Views   Citations

Abstract

This study presents a reduced-order modeling and optimal control framework for a shear-flow instability system derived from the incompressible Navier-Stokes equations via linear stability analysis of the Orr-Sommerfeld equation and subsequent center-manifold reduction. The resulting four-dimensional nonlinear dynamical system captures the essential interactions among the dominant Tollmien-Schlichting disturbance mode, the mean-flow correction induced by nonlinear Reynolds-stress feedback, and the actuator dynamics driven by an external control input. The model exhibits Hopf bifurcation, marking the transition from a stable equilibrium to sustained oscillations in the form of a stable limit cycle. Bifurcation analysis using numerical continuation reveals the existence of a Hopf point at (−0.000147, 0.052789, −0.001858, 0.105579, 0.105579), with a negative first Lyapunov coefficient indicating a supercritical bifurcation. Eigenvalue analysis of the Jacobian matrix at the bifurcation point confirms the presence of a conjugate imaginary pair crossing the stability boundary, validating the onset of oscillatory instability. An optimal control problem is formulated to minimize the energy of the dominant instability mode while penalizing control effort, with the control input simultaneously acting as a bifurcation parameter. The problem is solved using PYOMO.DAE with IPOPT under both unconstrained and Hopf-bifurcation-aware formulations. Results show that incorporating a Hopf constraint significantly reduces the objective function value and suppresses oscillatory control behavior. The study demonstrates that integrating bifurcation information into optimal control design enhances stability, reduces disturbance energy, and improves overall control efficiency in nonlinear fluid systems.

Keywords

Navier-Stokes Equations, Orr-Sommerfeld Stability, Hopf Bifurcation, Reduced-Order Modeling, Optimal Control

Share and Cite:

1. Introduction

The transition from laminar to turbulent flow remains one of the most fundamental and challenging problems in fluid mechanics because of its strong influence on drag generation, heat transfer, mixing efficiency, and aerodynamic performance. The nonlinear dynamics of fluid systems governed by the Navier-Stokes equations produce a wide range of instability mechanisms, including modal instabilities, transient growth, nonlinear interactions, and transition pathways. Understanding these mechanisms has therefore remained a central objective of theoretical, computational, and experimental fluid mechanics.

Early developments in hydrodynamic stability theory established the importance of nonmodal mechanisms in shear-flow transition by demonstrating that disturbances can experience significant transient amplification even when classical eigenvalue analysis predicts asymptotic stability [1]. The limitations of traditional eigenvalue approaches were further highlighted through investigations showing that hydrodynamic instability can arise from mechanisms beyond exponentially growing eigenmodes, emphasizing the importance of transient energy growth and non-normality of the governing operators [2]. Subsequent studies demonstrated that secondary instability mechanisms play a critical role in the breakdown of boundary layers and the progression toward turbulence [3].

The transition to turbulence in wall-bounded flows has been extensively investigated through nonlinear instability analysis, with particular emphasis on the interaction between disturbances, coherent structures, and nonlinear flow dynamics [4]. The Orr-Sommerfeld equation has remained one of the primary mathematical frameworks for analyzing viscous shear-flow instability, providing a foundation for predicting disturbance amplification and critical transition conditions [5]. Critical-layer dynamics were also identified as an important mechanism controlling instability development in shear flows, providing further insight into the interaction between disturbance waves and the underlying velocity field [6].

Advances in global instability analysis expanded the understanding of spatially developing flows by identifying instability mechanisms that cannot be captured using only local stability approaches [7]. Investigations of instability in pressure-driven channel flows and non-Newtonian fluids further demonstrated the applicability of Orr-Sommerfeld-based stability analysis to complex fluid systems with modified rheological behavior [8]. Comprehensive studies of stability and transition in shear flows subsequently provided a unified framework connecting linear stability theory, transient growth mechanisms, and nonlinear transition processes [9].

Numerical investigations continued to extend classical stability approaches to more complex configurations. Modified Orr-Sommerfeld formulations were employed to analyze stability characteristics of boundary-layer and channel-flow systems, demonstrating the adaptability of classical stability theory to generalized flow configurations [10]. Mathematical analyses of Prandtl boundary layers further advanced the understanding of spectral stability and the relationship between boundary-layer structures and instability development [11]. Historical studies of Sommerfeld’s contributions emphasized the long-lasting influence of classical mathematical approaches on modern fluid mechanics and stability analysis [12].

Global linear stability theory has continued to provide important insights into instability mechanisms associated with finite domains and complex flow geometries [13]. Further investigations of hydrodynamic instability mechanisms revealed the importance of nonlinear interactions, energy transfer, and multiple instability pathways in shear flows [14]. Energy amplification analysis demonstrated that different flow components contribute differently to transient growth mechanisms, providing deeper understanding of disturbance amplification in channel flows [15].

The role of Tollmien-Schlichting waves in boundary-layer transition has remained a major research topic, with reduced descriptions providing insight into the evolution of finite-amplitude disturbances and transition mechanisms [16]. Global stability analyses of boundary layers have further demonstrated the importance of spatial effects and complex eigenvalue structures in predicting instability onset [17]. Studies of open shear flows emphasized the distinction between convective and absolute instabilities and their influence on the downstream development of disturbances [18].

With increasing computational capabilities, high-fidelity numerical formulations and reduced-order descriptions have become important tools for analyzing fluid instabilities. Advanced formulations for viscous shear-flow stability problems have improved the computational treatment of complex eigenvalue problems and stability boundaries [19]. Modern linear stability theory has further integrated classical modal approaches with transient growth, nonlinear mechanisms, and transition prediction methodologies [20]. Mathematical studies of Prandtl-type boundary layers have continued to improve the theoretical understanding of instability mechanisms associated with high-Reynolds-number flows [21].

Recent investigations have extended stability analysis to more challenging nonlinear and transitional regimes. Stability studies of Prandtl-type shear flows for the Navier-Stokes equations have provided additional theoretical understanding of boundary-layer instability and spectral properties [22]. Energy growth mechanisms in shear flows have continued to be examined as a key factor influencing transition in nominally stable configurations [23]. Research on transition mechanisms in wall-bounded flows has highlighted the importance of nonlinear interactions and coherent structures in the evolution from organized disturbances to turbulent states [24]. Nonlinear boundary-layer instability studies have further demonstrated that finite-amplitude effects play a crucial role in determining transition behavior beyond linear predictions [25].

More recent developments have focused on integrating computational methods, nonlinear modeling, and advanced stability analysis techniques. Investigations of shear-flow instability and transition have emphasized the need for combined approaches that incorporate linear theory, nonlinear dynamics, and computational modeling [26]. Advances in computational hydrodynamic stability methods have enabled more efficient analysis of complex instability problems and improved prediction of transition characteristics [27]. Studies of absolute and convective instabilities have further clarified the mechanisms governing disturbance propagation and instability development in open shear flows [28].

Despite significant progress in hydrodynamic stability theory, several challenges remain in connecting instability analysis with active control strategies. Traditional stability analyses generally identify instability thresholds and dominant mechanisms but do not directly provide strategies for modifying system behavior after the onset of instability. Reduced-order models provide a promising approach by preserving dominant physical mechanisms while enabling efficient nonlinear analysis, bifurcation characterization, and control optimization.

The present work develops a reduced-order dynamical framework to analyze and control instability-driven oscillations in shear-flow systems. The model is derived from the dominant instability dynamics of the underlying fluid system and is formulated to capture nonlinear saturation, oscillatory behavior, and stability transitions. By integrating reduced-order modeling, bifurcation analysis, and optimal control, this framework provides a computationally efficient approach for understanding and manipulating transition phenomena in engineering flow systems.

The rest of this paper is organized as follows. First, the model equations are presented, followed by a description of the numerical procedures used. The results discussion and conclusions are then presented.

2. Model Equations

The model represents a low-dimensional dynamical reduction of boundary-layer transition dynamics derived from the incompressible Navier-Stokes equations through a modal decomposition based on the leading instability mechanism captured by the Orr-Sommerfeld equation. The governing idea is that near the critical Reynolds number, the flow is dominated by a small number of coherent structures associated with Tollmien-Schlichting wave instabilities, allowing the full infinite-dimensional fluid system to be approximated by a finite set of interacting modes. In this reduced framework, the variables $x_v\left(t\right)$ and $y_v\left(t\right)$ represent the real and imaginary components of the dominant unstable Orr-Sommerfeld eigenmode, corresponding physically to the amplitude and phase dynamics of the primary disturbance wave. Their evolution captures the onset of linear instability through exponential growth or decay determined by the effective growth rate, as well as nonlinear saturation arising from convective self-interaction, which limits unbounded amplification and leads to finite-amplitude oscillatory states (limit cycles) characteristic of Hopf bifurcation behavior. The variable $m_v\left(t\right)$ represents a slowly evolving mean-flow correction induced by Reynolds-stress feedback from the fluctuating disturbance field. This component models the redistribution of momentum in the boundary layer due to finite-amplitude wave activity and provides a coupling mechanism through which disturbances modify the base flow that supports them. The actuator state $u_c\left(t\right)$ represents a first-order dynamical model of external control action, introducing a finite response time between the applied input $v\left(t\right)$ and its effect on the flow field. Importantly, the control input $v\left(t\right)$ plays a dual role: it acts both as a direct forcing term on the actuator dynamics and as a parameter that modifies the effective linear stability of the flow by shifting the growth rate of the dominant Orr-Sommerfeld mode. As a result, the model captures not only natural transition dynamics but also controlled modification of the underlying instability mechanism. Overall, the coupled system describes the interaction between linear instability growth, nonlinear saturation, mean-flow distortion, and external control within a unified framework, providing a tractable representation of transition dynamics near a Hopf bifurcation in shear flows. Starting from the incompressible Navier-Stokes equations, consider a parallel shear base flow $U\left(y\right)$ representing a boundary-layer profile.

The present model considers an incompressible two-dimensional viscous boundary-layer flow developing over a flat plate with a parallel base velocity profile (U(y)). The analysis is restricted to a neighborhood of the critical Reynolds number, where the leading Orr-Sommerfeld eigenvalue approaches the imaginary axis and the first Tollmien-Schlichting instability becomes dominant. Disturbances are assumed to possess a single streamwise wavenumber, corresponding to the most unstable linear mode, while spanwise variations are neglected for simplicity. The classical no-slip and impermeability conditions are imposed at the wall, and all disturbance quantities decay to zero in the free stream. The control input represents an externally applied actuation, such as wall blowing/suction or localized body forcing, whose primary effect is to modify the effective linear growth rate of the dominant Orr-Sommerfeld mode. This same variable serves both as the bifurcation parameter governing the onset of instability and as the control variable used to suppress oscillatory flow behavior.

The total velocity field is decomposed as

$u=U\left(y\right)+u^{\prime }\left(x,y,t\right)$ (1)

where $U\left(y\right)$ is the base flow and $u^{\prime }$ denotes small perturbations. Substituting this decomposition into the Navier-Stokes equations and linearizing about the base flow yields the governing equations for infinitesimal disturbances. Assuming normal-mode perturbations of the form

$\begin{array}{l}{u}^{\prime }\left(x,y,t\right)=\widehat{u}\left(y\right)\text{exp}\left(i\alpha x+\lambda t\right)+\text{c}\text{.c}.\\ \lambda =\mu +i\omega \end{array}$(2)

and the system reduces to the Orr-Sommerfeld equation, which takes the form ${ℒ}_{OS}\varphi =\lambda ℳ\varphi$ , where the eigenvalue $\lambda =\mu +i\omega$ governs temporal growth and oscillation of disturbances. The condition $\mu =0$ corresponds to the onset of instability. Hence we have

$\begin{array}{l}{u}^{\prime }\left(x,y,t\right)\approx A\left(t\right){\varphi }_{OS}\left(y\right){\text{e}}^{i\alpha x}+\text{c}\text{.c}\\ A\left(t\right)=x_v\left(t\right)+i{y}_v\left(t\right)\end{array}$ (3)

Near the critical Reynolds number where $\mu \approx 0$ , the flow is dominated by a single unstable Orr-Sommerfeld eigenmode, allowing the disturbance field to be approximated as

$u^{\prime }\left(x,y,t\right)\approx A\left(t\right){\varphi }_{OS}\left(y\right){\text{e}}^{i\alpha x}+\text{c}\text{.c}$ (4)

Projecting the Navier-Stokes equations onto this dominant eigenfunction and retaining leading nonlinear interactions yields a center-manifold reduction governing the complex amplitude $A\left(t\right)$ . The resulting evolution equation is written using standard time-derivative notation as

$\frac{\text{d}A}{\text{d}t}=\left(\mu +i\omega \right)A-\beta {\left|A\right|}^{2}A+\gamma u_c,$ (5)

where the cubic term arises from nonlinear convective effects in Navier-Stokes, and $u_c$ represents external actuation that modifies the disturbance dynamics. Writing $A\left(t\right)=x_v\left(t\right)+i{y}_v\left(t\right)$ and separating real and imaginary parts yields two coupled real equations describing the evolution of the Tollmien-Schlichting wave components.

To capture feedback of disturbances on the mean flow, a slowly varying correction $m_v\left(t\right)$ is introduced via projection of the Reynolds-averaged Navier-Stokes equations, leading to

$\frac{\text{d}{m}_v}{\text{d}t}=-\delta m_v-\kappa \left(x_v^2+y_v^2\right),$ (6)

where the quadratic term represents energy transfer from fluctuations to the mean flow. Finally, actuator dynamics are modeled using a first-order lag equation,

$\frac{\text{d}{u}_c}{\text{d}t}=-\sigma u_c+v\left(t\right),$ (7)

where $v\left(t\right)$ is the control input. In addition, control modifies the effective linear growth rate through $\mu ={\mu }_0+v\left(t\right)$ , representing base-flow modification that shifts the Orr-Sommerfeld eigenvalue.

The next step that involves collecting all projected contributions. In this step, the full Navier-Stokes perturbation dynamics are consistently reduced into a finite-dimensional ODE system by projecting onto selected modes and retaining only leading-order physically relevant terms.

Starting from the incompressible Navier-Stokes equations, we write the decomposition

$u\left(x,y,t\right)=U\left(y\right)+u^{\prime }\left(x,y,t\right),$ (8)

substitute into Navier-Stokes, and separate linear and nonlinear contributions. After linearization, the dynamics of infinitesimal perturbations are governed by the Orr-Sommerfeld equation, which defines a spectrum of eigenmodes ${\varphi }_j\left(y\right)$ with eigenvalues ${\lambda }_j$ .

The key reduction step is to assume that near the instability threshold only one critical eigenmode dominates, so the perturbation field is expanded as

$u^{\prime }\left(x,y,t\right)\approx A\left(t\right){\varphi }_{OS}\left(y\right){\text{e}}^{i\alpha x}+\text{c}\text{.c}$ (9)

and possibly additional slow modes representing mean-flow correction. The phrase “collecting projected contributions” refers to applying a Galerkin projection (or center-manifold projection) onto these chosen basis functions. Concretely, this means taking the inner product of the Navier-Stokes residual with the adjoint eigenfunctions and integrating over the spatial domain, thereby converting the PDE into evolution equations for the modal amplitudes.

In this projection, each physical mechanism contributes a distinct term in the reduced ODE. The linear growth and oscillation terms come directly from the Orr-Sommerfeld eigenvalue $\lambda =\mu +i\omega$ , producing the term $\left(\mu +i\omega \right)A$ . The nonlinear cubic term arises from quadratic convective interactions $\left(u^{\prime }\cdot \nabla \right)u^{\prime }$ , which, after projection, generate a self-interaction proportional to ${\left|A\right|}^{2}A$ with coefficient $\beta$ . This term represents energy transfer among modes and ultimately saturates growth. The mean-flow correction term originates from Reynolds stress feedback in the averaged Navier-Stokes equations. When projecting onto the mean mode, quadratic products of fluctuations such as $x_v^2+y_v^2$ drive slow evolution of the base flow, producing a damping-type equation for $m_v$ . Finally, the control contribution enters as an external forcing term added before projection; depending on how it modifies the momentum equations, it either projects onto the same unstable eigenmode (entering as $\gamma u_c$ ) or modifies the base flow, effectively shifting the eigenvalue as $\mu ={\mu }_0+v\left(t\right)$ . Summarizing this amounts to summing all these physically distinct projected effects—linear stability from Orr-Sommerfeld, nonlinear self-interaction from Navier-Stokes convection, mean-flow feedback from Reynolds stresses, and external forcing/control—into a closed low-dimensional dynamical system governing $\left(x_v,y_v{m}_v,u_c\right)$ .

Additionally we define (in order to be able to control $\mu$ )

$\mu ={\mu }_0+v\left(t\right)$ (10)

Collecting all projected contributions and making the substitution for $\mu$ yields the reduced four-dimensional nonlinear system

$\begin{array}{l}\frac{\text{d}{x}_v}{\text{d}t}=\left({\mu }_0+v\right)x_v-\omega y_v-\beta \left(x_v^2+y_v^2\right)x_v+\gamma u_c,\\ \frac{\text{d}{y}_v}{\text{d}t}=\omega x_v+\left({\mu }_0+v\right)y_v-\beta \left(x_v^2+y_v^2\right)y_v,\\ \frac{\text{d}{m}_v}{\text{d}t}=-\delta m_v-\kappa \left(x_v^2+y_v^2\right),\\ \frac{\text{d}{u}_c}{\text{d}t}=-\sigma u_c+v\left(t\right),\end{array}$ (11)

which constitutes a center-manifold reduction of the Navier-Stokes system near an Orr-Sommerfeld instability, retaining nonlinear saturation, mean-flow feedback, and control-induced bifurcation shifting. Table 1 and Table 2 give the variable and parameter details.

The mean-flow correction state $m_v$ is included to quantify the slow distortion of the base flow produced by Reynolds-stress feedback from the finite-amplitude disturbance field. In the present reduced-order model, this correction is treated as a passive diagnostic variable whose evolution is driven by the disturbance energy but whose feedback on the dominant Orr-Sommerfeld mode is assumed to be of higher order near the Hopf bifurcation. Consequently, a one-way coupling is adopted in which the instability dynamics determine the mean-flow evolution, while the influence of the mean-flow correction on the leading instability mode is neglected. This assumption is consistent with a first-order center-manifold reduction in the vicinity of the critical Reynolds number, where the dominant instability dynamics are governed primarily by the critical Orr-Sommerfeld mode and higher-order mean-flow feedback effects are comparatively small. The variable $m_v$ nevertheless provides useful information regarding the evolution of the mean-flow distortion and can be incorporated as a two-way coupling term in higher-order reduced models if required.

More details of the derivation can be found in Appendix.

Table 1. Model variables.

Symbol	Variable Name	Physical Interpretation	Units
$x_v$	Real modal amplitude	Real component of the dominant Orr-Sommerfeld (Tollmien-Schlichting) instability mode	Dimensionless
$y_v$	Imaginary modal amplitude	Imaginary component of the dominant Orr-Sommerfeld (Tollmien-Schlichting) instability mode	Dimensionless
$m_v$	Mean-flow correction amplitude	Slow mean-flow distortion induced by Reynolds-stress feedback	Dimensionless
$u_c$	Actuator state	Dynamic response of the flow-control actuator	Dimensionless
$v$	Control input/bifurcation parameter	External control action that shifts the effective growth rate of the instability	Dimensionless
$t$	Time	Independent time variable	h

Table 2. Model parameters.

Symbol	Parameter Name	Value	Units
${\mu }_0$	Base linear growth-rate parameter	−0.1	h−1
$\omega$	Oscillation frequency of the dominant Orr-Sommerfeld mode	1.0	h−1
$\beta$	Nonlinear saturation coefficient	1.0	h−1
$\gamma$	Control coupling coefficient	0.5	h−1
$\delta$	Mean-flow relaxation coefficient	0.3	h−1
$\kappa$	Mean-flow feedback coefficient	0.2	h−1
$\sigma$	Actuator decay coefficient	1.0	h−1

The reduced-order model is obtained through modal projection and center-manifold reduction of the Navier-Stokes and Orr-Sommerfeld equations. Consequently, the modal amplitudes $x_v$ , $y_v$ , $m_v$ , and $u_c$ are normalized quantities and are therefore dimensionless. The coefficients ${\mu }_0$ , $\omega$ , $\beta$ , $\gamma$ , $\delta$ , $\kappa$ , and $\sigma$ have units of inverse time to maintain dimensional consistency in the governing ordinary differential equations.

The reduced-order model (ROM) developed in this study is intended as a nonlinear dynamical framework for investigating stability transitions, bifurcation behavior, and control strategies rather than as a direct replacement for the full Orr-Sommerfeld stability formulation. The reduced model captures the essential dynamical features associated with instability onset through the evolution of dominant modes and the resulting nonlinear amplitude interactions. In particular, the linearization of the reduced-order system near the equilibrium state provides the local growth rate and oscillation frequency through the eigenvalues of the reduced Jacobian,

${\lambda }_{ROM}={\sigma }_{ROM}\pm i{\omega }_{ROM}$ (12)

where ${\sigma }_{ROM}$ represents the modal growth or decay rate and ${\omega }_{ROM}$ represents the oscillation frequency. These quantities define the local stability characteristics and identify the transition from stable to oscillatory dynamics through the Hopf bifurcation condition,

${\sigma }_{ROM}=0,{\omega }_{ROM}\ne 0$ (13)

The reduced model should therefore be interpreted as a bifurcation-consistent normal-form representation of the dominant instability mechanism. Quantitative prediction of the complete spectrum of the underlying Navier-Stokes operator would require direct comparison with the leading Orr-Sommerfeld eigenvalues, which is beyond the scope of the present reduced-order framework. The primary objective of the ROM is to provide an efficient computational platform for nonlinear stability analysis and optimal control near critical operating conditions.

Additional details of the center-manifold reduction and Galerkin projection leading from the Navier-Stokes and Orr-Sommerfeld equations to the reduced-order model are provided in Appendix A.

3. Bifurcation Analysis and Optimal Control

Bifurcation Analysis

Continuation and bifurcation computations were carried out using the MATLAB-based software MATCONT. Bifurcation analysis provides important insights into the mechanisms underlying the existence of multiple steady states and self-sustained oscillations in nonlinear dynamical systems. In particular, branch points and limit points are associated with the emergence of multiple steady-state solutions, whereas oscillatory dynamics and periodic solutions originate from Hopf bifurcations. The package, originally developed and subsequently enhanced by several researchers, enables the systematic identification of limit points (LP), branch points (BP), and Hopf bifurcation points (H) in nonlinear dynamical models [29]-[31]. This program identifies Limit points (LP), branch points (BP), and Hopf bifurcation points(H) for a system of ordinary differential equations

$\frac{\text{d}x}{\text{d}t}=f\left(x,\alpha \right)$ (14)

$x\in R^n$ where the bifurcation parameter is $\alpha$ .

Optimal Control

Pyomo. dae [32] is used for the Optimal Control calculations. Pyomo with its Pyomo. DAE module provides an efficient framework for the formulation and solution of dynamic optimization problems involving systems of differential and algebraic equations. In Pyomo. DAE, the continuous-time dynamic model is transformed into a nonlinear programming (NLP) problem through discretization techniques such as orthogonal collocation, allowing the resulting optimization problem to be solved using standard NLP solvers. The NLP is solved using IPOPT [33].

The optimal control problem was formulated using Pyomo. DAE framework, where the nonlinear dynamic model was coupled with the manipulated control variable to determine the optimal operating trajectory. The general dynamic optimization problem was formulated as

$\underset{u\left(t\right)}{\text{min}}J=\Phi \left(x\left(t_f\right)\right)+{\displaystyle \underset{0}{\overset{t_f}{\int }}L\left(x\left(t\right),u\left(t\right)\right)\text{d}t}$ (15)

subject to the nonlinear dynamic system

$\frac{\text{d}x\left(t\right)}{\text{d}t}=f\left(x\left(t\right),u\left(t\right),p\right),0\le t\le t_f$ (16)

where $x\left(t\right)$ represents the vector of dynamic states, $u\left(t\right)$ represents the control variable, and $p$ represents the model parameters including the bifurcation parameter.

The optimization was performed over the finite time horizon $0\le t\le t_f$ where $t_f$ is the final optimization time. The initial conditions were imposed as equality constraints:

$x_i\left(0\right)=x_{i,0},i=1,2,\cdots ,n$ where $x_{i,0}$ represents the initial value of the $i^{\text{th}}$ state variable.

The manipulated control variable was restricted within physically meaningful operating limits $u_{\text{min}}\le u\left(t\right)\le u_{\text{max}}$ where $u_{\text{min}}$ and $u_{\text{max}}$ represent the lower and upper allowable control limits, respectively.

State variables were constrained through path constraints when required:

$x_{i,\text{min}}\le x_i\left(t\right)\le x_{i,\text{max}},i=1,2,\cdots ,n$

No terminal constraints were imposed in the present optimization formulation. Therefore, the final state was determined by the optimal solution:

$x\left(t_f\right)=x_f$

where $x_f$ represents the terminal state obtained from the optimal trajectory.

The continuous-time optimization problem was transformed into a nonlinear programming (NLP) problem using Pyomo. DAE orthogonal collocation method on finite elements. The time domain was discretized as

$0=t_0<t_1<\cdots <t_{n_{fe}}=t_f$

where $n_{fe}$ represents the number of finite elements. Within each finite element, the differential variables were approximated using collocation polynomials with $n_{cp}$ collocation points per finite element.

The resulting NLP problem was expressed as

$\underset{x,u}{\text{min}}J\left(x,u\right)$ (17)

subject to

$F\left(x,u\right)=0;G\left(x,u\right)\le 0$ (18)

where $F$ represents the discretized dynamic equations generated by Pyomo. DAE and $G$ represent the inequality constraints associated with state and control limitations. The resulting NLP problem was solved using IPOPT. The stability constraint was incorporated through the neural-network approximation of the maximum real eigenvalue:

${\lambda }_{\text{max}}=\text{max}\left(\text{Re}\left(\lambda \left(J\left(x\right)\right)\right)\right)$ (19)

which was replaced by the differentiable neural-network approximation ${\widehat{\lambda }}_{\text{max}}=NN\left(x,p\right)$

The stability-aware objective function was then formulated as

$J_{\text{total}}=J+\eta \text{smooth\_max}\left(0,{\widehat{\lambda }}_{\text{max}}-ϵ\right)$ (20)

where $\eta$ is the stability penalty coefficient and $ϵ$ is a small positive margin introduced to improve numerical robustness near the Hopf boundary. The smooth_max function is a differentiable approximation of the maximum operator that provides a continuous transition between stable and unstable regions, allowing the instability penalty to be smoothly activated when the predicted maximum eigenvalue exceeds the stability threshold while maintaining the differentiability required by IPOPT.

This formulation avoids repeated eigenvalue calculations during optimization and provides a smooth, differentiable stability measure compatible with IPOPT-based gradient optimization.

4. Formation of Stability Dataset from MATCONT Results

A stability dataset was developed from numerical continuation calculations performed in MATCONT. The stability dataset consists of rows, each representing a continuation point from an equilibrium branch. Each row contains the state variables, the bifurcation parameter, and a stability measure. The stability measure is a numerical value derived from the Jacobian matrix. The Jacobian matrix is computed numerically at each equilibrium point. The eigenvalues are then computed automatically using MATLAB. The maximum value of the real part of these eigenvalues is then computed as a scalar stability measure.

The stability measure is computed using “eig_real_max = max(real(eigvals));” in MATLAB. The stability measure is a quantitative metric in which negative values indicate locally asymptotically stable equilibria, positive values indicate instability, and a zero crossing indicates a Hopf bifurcation. The stability dataset is then saved as a CSV file. The dataset can then be used in subsequent computational calculations to perform classification or regression to identify stability boundaries or approximate bifurcations.

Neural Network Surrogate for Stability Prediction

Direct embedding of eigenvalue calculations into IPOPT-based optimal control is impractical for several reasons: 1) computing eigenvalues at each time step is computationally expensive, 2) the mapping from states to the maximum eigenvalue is non-smooth near eigenvalue crossings, and 3) symbolic differentiation of eigenvalues is challenging.

Prior to neural-network training, all input variables were standardized to improve numerical conditioning and training stability. Let $x_{\text{raw}}$ denote the vector of state variables and bifurcation parameters = obtained from the stability dataset. For each input variable $j$ , the training mean ${\mu }_j$ and training standard deviation ${\sigma }_j$ were computed over all training samples. The training mean for the input variable $j$ is defined as the arithmetic average over all training samples as

${\mu }_j=\frac{1}{N}{\displaystyle \sum _{k=1}^N{x}_j^{\left(k\right)}}$ and the training standard deviation for the input variable $j$ is defined as: ${\sigma }_j=\frac{1}{N}{\left({\displaystyle \sum _{k=1}^N{x}_j^{\left(k\right)}-{\mu }_j}\right)}^2$

The standardized inputs were defined as

$x_j=\frac{x_{\text{raw},j}-{\mu }_j}{\max \left({\sigma }_j,{\epsilon }_{\text{scale}}\right)}$ (21)

where ${\epsilon }_{\text{scale}}$ is a small positive regularization parameter introduced to prevent numerical singularities associated with extremely small variances. This transformation ensures that all inputs remain properly scaled while avoiding excessively large neural-network activation arguments during optimization. The normalization procedure improves neural-network conditioning and enhances the robustness of gradient-based optimization. The vectors $\mu$ and $\sigma$ computed during training were stored and embedded identically within the Pyomo optimal-control formulation to ensure consistency between neural-network training and deployment.

To overcome these limitations, a feedforward neural network is trained to approximate the maximum eigenvalue as a smooth function of the system state and bifurcation parameter. A typical architecture employs the hyperbolic tangent (tanh) as a smooth activation function. If the input vector is denoted by $x$ , which represents the scaled variables, then the network is defined as:

The vectors $\mu ,\sigma$ computed during training were stored and embedded identically within the Pyomo optimal control formulation to ensure consistency between neural network training and deployment.

To avoid repeated eigenvalue computations during optimization, a feedforward neural network is trained to approximate the maximum real eigenvalue as a smooth function of the system states and bifurcation parameter. Using hyperbolic tangent activation functions ensures smooth differentiability required by IPOPT. If the input vector is denoted by $x$ , which are the scaled variables, the network architecture is defined as

$\begin{array}{l}{z}_1=\tanh \left(W_{1}x+b_1\right)\\ z_2=\tanh \left(W_2{z}_1+b_2\right)\\ \lambda \text{\_max\_NN}=W_3{z}_2+b_3\end{array}$ (22)

where $W_i$ and $b_i$ denote the weight matrices and bias vectors, respectively. The hidden-layer variables $z_1$ and $z_2$ represent nonlinear transformations of the input variables and intermediate features. The final output ${\widehat{\lambda }}_{\text{max}}$ provides a smooth approximation of the spectral abscissa. Because tanh is infinitely differentiable, the network is fully smooth, guaranteeing the availability of first and second derivatives required by IPOPT. Without biases, the network output would be constrained to pass through the origin, limiting flexibility.

The hidden-layer outputs $z_1$ , $z_2$ , represent nonlinear combinations of the inputs and previous-layer features, respectively. Each element of $z_1$ is a smoothed combination of the original inputs, while each element of $z_2$ encodes more abstract patterns extracted from $z_1$ . The final output $\lambda \text{\_max\_NN}$ provides a smooth approximation of the maximum real eigenvalue, enabling efficient and differentiable stability evaluation within the optimal control problem. The integration into optimal control is done using a soft penalty formulation where we use a smooth_max function that converts $\lambda$ into a smooth, nonnegative penalty that only “activates” when the system is unstable:

$s_{\max }\left(\lambda +{\epsilon }_1\right)=\text{smooth\_max}\left(\lambda +{\epsilon }_1\right)=\frac{\left(\lambda +{\epsilon }_1\right)+\sqrt{{\left(\lambda +{\epsilon }_1\right)}^2+{\epsilon }_1}}{2}$ (23)

$\lambda$ is the neural network’s predicted maximum eigenvalue at the current state and parameter, while $\epsilon$ is a small positive safety margin to ensure differentiability. The soft penalty formulation involves the new objective function, where the original objective function ${\displaystyle \sum optv\left(t\right)}$ is modified to ${\displaystyle \sum {\left(optv\left(t\right)+{\alpha }_{\text{hopf}}\cdot s_{\max }\left(\lambda +{\epsilon }_1\right)\right)}^2}$ ${\alpha }_{\text{Hopf}}$ controls how aggressively instability is penalized, and prevents numerical issues at exactly $\lambda =0$ and slightly shifts the stability boundary. When ${\alpha }_{\text{Hopf}}$ is 0, no Hopf constraint is implemented. The goal is to minimize the objective function value while ensuring differentiability for IPOPT and avoiding non-smoothness in the optimization. This approach avoids repeated eigenvalue computations and provides a smooth, differentiable surrogate suitable for gradient-based optimization. Furthermore, this enables the incorporation of stability constraints into optimal control without explicitly computing eigenvalues during optimization.

To facilitate integration into optimal control, a stability dataset is constructed from MATCONT continuation results. At each equilibrium point, the Jacobian is evaluated and its eigenvalues computed. The stability metric is defined as the maximum real part of the eigenvalues (${\lambda }_{\text{max}}=\max \Re \left(\lambda \left(J\right)\right)$ ). Negative values indicate stability, positive values indicate instability, and zero crossings correspond to Hopf bifurcations.

5. Results

The bifurcation performed with MATCONT revealed the existence of a Hopf bifurcation point at $\left(x_v,y_v,m_v,u_c,v\right)$ values of (−0.000147, 0.052789, −0.001858, 0.105579, 0.105579) with a first Lyapunov coefficient = −1.999622e+00. Figure 1 shows the limit cycle.

Figure 1. Limit cycle.

The analytical Jacobian matrix is

$Ja=\left[\begin{array}{cccc}\mu -\beta \left(3x_v^2+y_v^2\right)& -\omega -2\beta x_v{y}_v& 0& \gamma \\ \omega -2\beta x_v{y}_v& \mu -\beta \left(x_v^2+3y_v^2\right)& 0& 0\\ -2\kappa x_v& -2\kappa y_v& -\delta & 0\\ 0& 0& 0& -\sigma \end{array}\right].$ (24)

Using the parameter values ${\mu }_0=-0.1,\omega =1,\beta =1,\gamma =0.5,\delta =0.3,\kappa =0.2,\sigma =1$ , and the variables at the Hopf point $\left(x_v,y_v,m_v,u_c,v\right)=\left(-0.000147,0.052789,-0.001858,0.105579,0.105579\right)$ , we obtain$\mu ={\mu }_0+v=-0.1+0.105579=0.005579$ . Substituting these values into the Jacobian gives the numerical Jacobian matrix

$Jn=\left[\begin{array}{cccc}0.00279226& -0.99998448& 0& 0.5\\ 1.00001552& -0.00278106& 0& 0\\ 0.00005880& -0.02111560& -0.3& 0\\ 0& 0& 0& -1\end{array}\right].$ (25)

The corresponding eigenvalues are −1, −0.3, 0.999i and −0.999i. The presence of the conjugate imaginary eigenvalues confirm the existence of the Hopf bifurcation point. Table 3 summarizes the CSV file that was developed after the bifurcation analysis.

Table 3. Contents of CSV file developed after bifurcation analysis.

$x_v$	$y_v$	$m_v$	$u_c$	$\nu$	Real max (eig)
−0.0027345	0.071565	−0.0034193	0.14334	0.14334	0.033081
−0.0021854	0.068269	−0.0031103	0.13668	0.13668	0.027347
−0.0016757	0.064965	−0.0028155	0.13002	0.13002	0.02157
−0.0012055	0.061654	−0.0025351	0.12336	0.12336	0.01575
−0.00077535	0.058337	−0.0022692	0.11669	0.11669	0.009887
−0.00038547	0.055014	−0.0020178	0.11003	0.11003	0.00398
−0.00014739	0.052789	−0.0018578	0.10558	0.10558	5.2868e−06
−3.6213e−05	0.051686	−0.001781	0.10337	0.10337	−0.0019708
0.00027212	0.048354	−0.0015588	0.096711	0.096711	−0.0079658
0.00053923	0.045018	−0.0013513	0.090049	0.090049	−0.014005
0.00076486	0.041679	−0.0011585	0.083387	0.083387	−0.020089
0.00094875	0.038339	−0.0009805	0.076724	0.076724	−0.026217
0.0010907	0.034997	−0.00081731	0.070062	0.070062	−0.03239
0.0011904	0.031654	−0.00066894	0.063398	0.063398	−0.038609
0.0012477	0.028312	−0.00053543	0.056735	0.056735	−0.044872
0.0012624	0.024971	−0.00041677	0.05007	0.05007	−0.05118
0.0012344	0.021632	−0.00031299	0.043406	0.043406	−0.057533
0.0011636	0.018296	−0.00022407	0.03674	0.03674	−0.063932
0.0010497	0.014963	−0.00015	0.030074	0.030074	−0.070376

The optimal control minimizes the function ${\displaystyle \sum optv\left(t\right)}={\displaystyle \sum {\left(x_v\left(t\right)\right)}^2}+{\left(y_v\left(t\right)\right)}^2+0.5\ast {\left(v\left(t\right)\right)}^2$ with and without the Hopf bifurcation constraint. The objective function minimizes the energy of the dominant Tollmien-Schlichting (Orr-Sommerfeld) instability mode, represented by $x_v^2+y_v^2$ , while simultaneously minimizing the control effort required to modify the flow stability through the bifurcation parameter $v\left(t\right)$ . me. To investigate the effect of bifurcation-aware operation, an optimal control problem was solved both with and without a Hopf-bifurcation-avoidance constraint using PYOMO.DAE coupled with IPOPT. PYOMO.DAE with IPOPT was used. When no Hopf constraint was implemented, ${\alpha }_{\text{Hopf}}=0$ , in ${\displaystyle \sum {\left(optv\left(t\right)+{\alpha }_{\text{hopf}}\cdot s_{\max }\left(\lambda +{\epsilon }_1\right)\right)}^2}$ the obtained value of ${\displaystyle \sum optv\left(t\right)}$ was 49.427. For a Hopf constraint, ${\alpha }_{\text{Hopf}}=0.5$ , the obtained value of ${\displaystyle \sum optv\left(t\right)}$ was 5.0. This shows a considerable decrease in the minimized function, indicating the effect of the Hopf bifurcation constraint. This corresponds to an 89.9% reduction in the minimized objective function value. (Figures 2(a)-(d)) show the optimal control profiles without and with the Hopf constraint. A comparison of (Figures 2(b)-(d)) shows the considerable reduction in the oscillatory spikes in the control profiles.

Figure 2. (a): State variables for optimal control without Hopf constraint; (b): Control variables optimal control without Hopf constraint; (c): State variables for optimal control with Hopf constraint; (d): Control variables optimal control with Hopf constraint.

6. Discussion

The bifurcation analysis performed using MATCONT confirmed that the reduced-order model derived from the Navier-Stokes and Orr-Sommerfeld equations exhibits a Hopf bifurcation. The computed Hopf point was located at $\left(x_v,y_v,m_v,u_c,v\right)=\left(-0.000147,0.052789,-0.001858,0.105579,0.105579\right)$ , indicating that the instability occurs when the control/bifurcation parameter reaches a value of approximately $v=0.105579$ . The existence of this bifurcation demonstrates that the reduced model successfully captures the transition from a steady equilibrium state to sustained oscillatory behavior, which is characteristic of hydrodynamic instabilities in shear flows. Since the model was obtained through a center-manifold reduction of the Navier-Stokes equations near an Orr-Sommerfeld instability, the detected Hopf point represents the onset of oscillatory Tollmien-Schlichting wave dynamics within the reduced-order framework.

The first Lyapunov coefficient at the Hopf point was found to be $l_1=-1.999622$ . The negative value of the first Lyapunov coefficient is particularly significant because it indicates that the Hopf bifurcation is supercritical. In a supercritical Hopf bifurcation, the equilibrium loses stability smoothly and gives rise to a stable limit cycle of small amplitude. This behavior is consistent with the numerical continuation results shown in Figure 1, where a stable periodic orbit emerges from the bifurcation point. Physically, this implies that once the critical instability threshold is crossed, the disturbance amplitudes do not grow without bound. Instead, nonlinear effects represented by the cubic saturation terms stabilize the oscillations and lead to finite-amplitude periodic behavior. Such saturation mechanisms are well known in fluid mechanics and represent the balance between linear instability growth and nonlinear energy redistribution.

Additional confirmation of the Hopf bifurcation was obtained through examination of the Jacobian matrix evaluated at the bifurcation point. Substitution of the parameter values and state variables into the analytical Jacobian produced the numerical Jacobian matrix from which the eigenvalues were calculated. The resulting eigenvalue spectrum contained two negative real eigenvalues, −1 and −0.3, together with a purely imaginary conjugate pair, $0.999i$ and $-0.999i$ . The appearance of a conjugate pair on the imaginary axis while the remaining eigenvalues remain in the left-half complex plane is the defining local condition for a Hopf bifurcation. The eigenvalue analysis therefore independently verifies the bifurcation detected by MATCONT. Furthermore, the imaginary part of the critical eigenvalues provides the oscillation frequency associated with the emerging limit cycle, linking the bifurcation directly to periodic Tollmien-Schlichting-type wave behavior represented by the variables $x_v$ and $y_v$ .

From a physical perspective, the Hopf bifurcation marks the transition from a stable flow regime to an oscillatory regime dominated by the leading Orr-Sommerfeld instability mode. The variables $x_v$ and $y_v$ represent the real and imaginary components of the dominant disturbance mode, while $m_v$ captures mean-flow modification induced by Reynolds-stress feedback and $u_c$ represents actuator dynamics. The emergence of a stable limit cycle therefore reflects the nonlinear interaction between disturbance growth, mean-flow distortion, and control dynamics. The reduced-order model successfully reproduces these interactions while remaining computationally tractable compared with the full Navier-Stokes equations.

The optimal control study further demonstrated the importance of accounting for bifurcation structure during control design. The objective function was formulated to minimize the energy of the dominant Orr-Sommerfeld instability mode together with the control effort required to modify the bifurcation parameter. Consequently, the optimization problem seeks a balance between suppressing instability-induced oscillations and avoiding excessive control action. This objective is particularly relevant in flow-control applications because aggressive control strategies may suppress disturbances but often require impractically large inputs. The inclusion of a control penalty therefore promotes more realistic operating conditions.

The comparison between the unconstrained and bifurcation-constrained optimal control solutions provides important insight into the role of Hopf bifurcation information in dynamic optimization. When no Hopf constraint was imposed, corresponding to $\alpha =0$ , the optimal objective function value was found to be 49.427. In this case, the optimizer seeks to minimize the objective without explicitly considering the proximity of the system to the instability boundary. Although a feasible solution is obtained, the resulting trajectories may operate near the bifurcation region where oscillatory dynamics are present. Such operation can lead to increased disturbance amplitudes, larger control fluctuations, and reduced robustness to parameter variations or external disturbances.

When the Hopf constraint was introduced, corresponding to $\alpha =5$ , the optimal objective function value decreased dramatically to 5.0. This reduction corresponds to approximately 89.9% relative improvement compared with the unconstrained case. The substantial decrease in the objective function indicates that incorporating bifurcation information into the optimization process enables the controller to steer the system away from regions of oscillatory instability. By maintaining operation within a dynamically favorable region of the parameter space, the optimization simultaneously reduces disturbance energy and control effort. The resulting operating conditions are therefore both more stable and more efficient.

The control profiles shown in Figures 2(a)-(d) provide further evidence of the benefits of bifurcation-aware optimization. In the absence of the Hopf constraint, the control signal exhibits stronger oscillatory variations and larger spikes. These fluctuations are symptomatic of operation near an instability boundary, where the controller must react more aggressively to suppress disturbance growth. In contrast, the control profiles obtained with the Hopf constraint are significantly smoother and exhibit reduced oscillatory behavior. The reduction in control spikes suggests that the system remains farther from the unstable region and therefore requires less corrective action. Such behavior is highly desirable in practical control applications because smoother control signals generally improve actuator longevity, reduce energy consumption, and enhance overall system reliability.

The results collectively demonstrate the value of integrating bifurcation analysis and optimal control within a unified computational framework. MATCONT provides detailed information regarding the location and nature of critical bifurcations, while PYOMO. DAE and IPOPT enable the incorporation of this information into a dynamic optimization problem. The resulting bifurcation-aware control strategy not only suppresses instability-related oscillations but also achieves substantial improvements in control performance. More broadly, the study illustrates how reduced-order models derived from the Navier-Stokes and Orr-Sommerfeld equations can bridge the gap between nonlinear stability analysis and optimal control design. The significant reduction in the objective function value achieved through Hopf-constrained optimization highlights the potential of this approach for controlling oscillatory fluid-dynamical systems and for developing operating strategies that explicitly account for underlying nonlinear stability boundaries.

The present study demonstrates a systematic framework that combines fluid-dynamical stability analysis, nonlinear bifurcation theory, and optimal control within a single computational methodology. By linking the Navier-Stokes equations to a reduced-order model through the Orr-Sommerfeld equation and center-manifold reduction, the work establishes a direct connection between fundamental flow-instability mechanisms and practical control design. This integration enables the identification of critical operating conditions associated with Hopf bifurcations and provides a means of actively steering the system away from undesirable oscillatory regimes.

A major contribution of this research is the incorporation of bifurcation information directly into the optimal control formulation. Traditional control strategies often focus solely on minimizing a performance objective without explicitly accounting for the proximity of the system to instability boundaries. In contrast, the proposed approach utilizes knowledge of the Hopf bifurcation structure to guide the optimization process toward dynamically stable operating regions. The resulting reduction in oscillatory behavior and control effort demonstrates that bifurcation-aware optimization can significantly improve system performance while enhancing operational robustness.

From a fluid-mechanics perspective, the study provides insight into the control of Tollmien-Schlichting-wave-driven instabilities in shear flows. The results show that enforcing a Hopf-bifurcation constraint can substantially reduce disturbance energy and suppress the emergence of large-amplitude oscillations. Such capabilities are important in applications involving laminar-flow preservation, transition delay, flow stabilization, and drag reduction. Although the present work employs a reduced-order model, the methodology can serve as a foundation for more detailed investigations involving higher-dimensional fluid models and computational fluid dynamics simulations.

The research also contributes to the broader field of nonlinear dynamical systems by illustrating how bifurcation analysis can be integrated with dynamic optimization. Many engineering systems, including fluid flows, chemical reactors, power systems, thermoacoustic combustors, biological systems, and environmental processes, exhibit Hopf bifurcations that lead to self-sustained oscillations. The framework developed in this study is therefore not limited to boundary-layer flows but may be adapted to a wide range of nonlinear systems in which oscillatory instabilities affect performance, safety, or efficiency.

An additional impact of the work lies in its computational methodology. The combination of MATCONT for bifurcation detection and PYOMO.DAE with IPOPT for optimal control provides a practical workflow for identifying critical stability boundaries and incorporating them into optimization-based decision-making. This workflow can be readily extended to more complex systems involving multiple bifurcations, parameter uncertainties, and advanced control architectures. Consequently, the study offers a useful template for future research aimed at coupling nonlinear dynamics, bifurcation theory, and optimal control.

Overall, the results demonstrate that explicit consideration of Hopf bifurcation behavior during optimization can lead to substantial improvements in system performance, stability, and control efficiency. By bridging the disciplines of fluid mechanics, nonlinear dynamics, and optimization, this research advances the development of bifurcation-aware control strategies and provides a foundation for the design of more reliable and efficient engineering systems operating near critical stability boundaries.

7. Conclusions

This study developed a bifurcation-aware optimal control framework for a reduced-order model of shear-flow instability derived from the incompressible Navier-Stokes equations using Orr-Sommerfeld stability theory and center-manifold reduction. The resulting four-dimensional system captures the essential physics of boundary-layer transition, including Tollmien-Schlichting wave growth, nonlinear saturation, mean-flow modification, and actuator-induced forcing.

From a fluid engineering perspective, the primary contribution of this work is the explicit connection between nonlinear stability theory and practical flow control design in systems relevant to boundary-layer transition, aerodynamic surfaces, and internal shear flows. In these engineered configurations, the onset of oscillatory instabilities is directly associated with increased drag, loss of efficiency, and degraded flow performance. By identifying and controlling the Hopf bifurcation structure of the reduced system, the present framework provides a systematic method for delaying or suppressing transition-related oscillations.

Bifurcation analysis using MATCONT revealed a Hopf bifurcation point separating steady and oscillatory regimes, with the computed first Lyapunov coefficient indicating a supercritical transition to stable limit-cycle oscillations. This behavior is representative of oscillatory boundary-layer instabilities observed in practical fluid systems, where small changes in operating conditions can trigger sustained unsteady flow.

The optimal control problem, implemented in PYOMO.DAE and solved using IPOPT, demonstrated that incorporating bifurcation information significantly improves control performance. When a Hopf-bifurcation-aware constraint was included, the system exhibited an approximately 89.9% reduction in the objective function compared to the unconstrained case. More importantly, the controlled dynamics showed suppressed oscillatory energy, smoother actuator effort, and operation maintained away from the instability boundary.

These results highlight that embedding nonlinear stability structure directly into the control formulation is not only mathematically consistent but also practically valuable for engineered fluid systems. In particular, the approach provides a pathway for designing control strategies that are informed by the underlying physics of transition rather than relying solely on linearized or heuristic stabilization methods.

Overall, the study demonstrates that reduced-order bifurcation modeling combined with optimal control offers a powerful framework for improving flow stability and performance in engineering applications where shear-flow instabilities govern transition and efficiency losses.

Data Availability Statement

All data used are presented in the paper.

Acknowledgements

Dr. Sridhar thanks Dr. Carlos Ramirez for encouraging him to write single-author papers.

Appendix

Near the critical Reynolds number $R{e}_c$ , the Orr-Sommerfeld operator possesses a critical eigenvalue pair

${\lambda }_c=\mu +i\omega$ (26)

with corresponding direct and adjoint eigenfunctions $\varphi \left(y\right)$ and ${\varphi }^{†}\left(y\right)$ . The perturbation velocity field is expanded as

$u^{\prime }\left(x,y,t\right)=A\left(t\right)\varphi \left(y\right){\text{e}}^{i\alpha x}+\overline{A}\left(t\right)\overline{\varphi }\left(y\right){\text{e}}^{-i\alpha x}+u_s$ (27)

where $A\left(t\right)$ is the complex amplitude of the critical Tollmien-Schlichting mode and $u_s$ contains higher-order stable modes. The bar on top indicates the complex conjugate.

The direct eigenfunction $\varphi \left(y\right)$ satisfies the Orr-Sommerfeld eigenvalue problem${ℒ}_{OS}\varphi =\lambda \varphi$ , while the corresponding adjoint eigenfunction satisfies ${ℒ}_{OS}^{†}{\varphi }^{†}=\overline{\lambda }{\varphi }^{†}$ . The direct and adjoint eigenfunctions are normalized according to $\langle {\varphi }^{†},\varphi \rangle =1$ where the inner product is defined as $\langle f,g\rangle ={\displaystyle \underset{0}{\overset{\infty }{\int }}{f}^{*}\left(y\right)g\left(y\right)\text{d}y}$ . This normalization uniquely fixes the amplitude of the reduced-order model and removes the arbitrary scaling of the Orr-Sommerfeld eigenfunctions.

Since all remaining modes are linearly stable near the bifurcation point, the center-manifold theorem implies that their dynamics are slaved to the critical amplitude $A$ . Consequently,

$u_s=H\left(A,\overline{A}\right)$ (28)

where $H$ is a smooth nonlinear function satisfying $H=O\left({\left|A\right|}^2\right)$ . Substituting the modal expansion into the Navier-Stokes equations and projecting onto the adjoint eigenfunction using the inner product $\langle f,g\rangle ={\displaystyle \underset{\Omega }{\int }{f}^{*}g\text{d}\Omega }$ , yields the solvability condition

$\langle {\varphi }^{†},\frac{\partial u^{\prime }}{\partial t}-L{u}^{\prime }-N\left(u^{\prime },u^{\prime }\right)\rangle =0$ (29)

where $L$ denotes the linearized Navier-Stokes operator and $N$ denotes the quadratic convective nonlinearity. Evaluation of the projection gives

$\frac{\text{d}A}{\text{d}t}=\left(\mu +i\omega \right)A-c{\left|A\right|}^{2}A+\gamma u_c+O\left({\left|A\right|}^5\right)$ (30)

where $c=\langle {\varphi }^{†},N_3\rangle$ is the Landau coefficient arising from the self-interaction of the critical mode.

Retaining the leading nonlinear contribution yields the Stuart-Landau amplitude equation

$\frac{\text{d}A}{\text{d}t}=\left(\mu +i\omega \right)A-c{\left|A\right|}^{2}A+\gamma u_c$ (31)

Writing $A=x_v+i{y}_v$ gives ${\left|A\right|}^2=x_v^2+y_v^2$ . Substituting into the Stuart-Landau equation and separating real and imaginary parts yields

$\begin{array}{l}\frac{\text{d}{x}_v}{\text{d}t}=\mu x_v-\omega y_v-\beta \left(x_v^2+y_v^2\right)x_v+\gamma u_c\\ \frac{\text{d}{y}_v}{\text{d}t}=\omega x_v+\mu y_v-\beta \left(x_v^2+y_v^2\right)y_v\end{array}$ (32)

where $\beta =\text{Re}\left(c\right)$ . The finite-amplitude disturbance modifies the base flow through Reynolds stresses. Averaging the Navier-Stokes equations over the fast oscillation period gives

$\frac{\partial U_m}{\partial t}=R\left(A,\overline{A}\right)-\delta U_m$ (33)

where $U_m$ denotes the mean-flow correction and R represents the Reynolds-stress forcing. Instead of solving for the entire function $U_m\left(y,t\right)$ , we assume that its shape is approximately fixed and can be represented by a dominant mode:

$U_m\left(y,t\right)=m_v\left(t\right)\psi \left(y\right)$ (34)

where $\psi \left(y\right)$ = chosen mean-flow basis function (mode shape), $m_v\left(t\right)$ = amplitude of that mode.

Now substitute this expansion into the Reynolds-averaged equation and take the inner product with an adjoint mode ${\psi }^{†}\left(y\right)$ ; $\langle {\psi }^{†},\text{mean-flow equation}\rangle$ . This operation removes the spatial dependence and leaves an ODE for the amplitude $m_v\left(t\right)$ . Projecting the equation for $\frac{\partial U_m}{\partial t}$ onto the dominant mean-flow mode $\psi \left(y\right)$ gives

$\frac{\text{d}{m}_v}{\text{d}t}=-\delta m_v-\kappa {\left|A\right|}^2$ (35)

where $\kappa =\langle {\psi }^{+},R\rangle$ . The actuator dynamics are modeled as

$\frac{\text{d}{u}_c}{\text{d}t}=-\sigma u_c+v$ (36)

Finally, assuming that the control modifies the effective growth rate of the Orr-Sommerfeld mode, $\mu ={\mu }_0+v$ , the reduced system becomes

$\begin{array}{l}\frac{\text{d}{x}_v}{\text{d}t}=\left({\mu }_0+v\right)x_v-\omega y_v-\beta \left(x_v^2+y_v^2\right)x_v+\gamma u_c\\ \frac{\text{d}{y}_v}{\text{d}t}=\omega x_v+\left({\mu }_0+v\right)y_v-\beta \left(x_v^2+y_v^2\right)y_v\\ \frac{\text{d}{m}_v}{\text{d}t}=-\delta m_v-\kappa \left(x_v^2+y_v^2\right)\\ \frac{\text{d}{u}_c}{\text{d}t}=-\sigma u_c+v\end{array}$ (37)

All parameter and variable details are in Table 1 and Table 2.

Conflicts of Interest

The author, Dr. Lakshmi N Sridhar, has no conflicts of interest.

References

[1]	Schmid, P.J. (2007) Nonmodal Stability Theory. Annual Review of Fluid Mechanics, 39, 129-162.[CrossRef]
[2]	Trefethen, L.N., Trefethen, A.E., Reddy, S.C. and Driscoll, T.A. (2001) Hydrodynamic Stability without Eigenvalues. Science, 261, 578-584.[CrossRef] [PubMed]
[3]	Herbert, T. (2002) Secondary Instability of Boundary Layers. Annual Review of Fluid Mechanics, 20, 487-526.[CrossRef]
[4]	Kerswell, R.R. (2005) Recent Progress in Understanding the Transition to Turbulence in a Pipe. Nonlinearity, 18, R17-R44.[CrossRef]
[5]	Orszag, S.A. (2006) Accurate Solution of the Orr-Sommerfeld Stability Equation. Journal of Fluid Mechanics, 50, 689-703.[CrossRef]
[6]	Maslowe, S.A. (2006) Critical Layers in Shear Flows. Annual Review of Fluid Mechanics, 18, 187-214.
[7]	Alizard, F. and Robinet, J. (2007) Spatially Convective Global Modes in a Boundary Layer. Physics of Fluids, 19, Article ID: 114105.[CrossRef]
[8]	Sahu, K.C., Valluri, P., Spelt, P.D.M. and Matar, O.K. (2007) Linear Instability of Pressure-Driven Channel Flow of a Newtonian and a Herschel-Bulkley Fluid. Physics of Fluids, 19, Article 122101.[CrossRef]
[9]	Schmid, P.J. and Henningson, D.S. (2008) Stability and Transition in Shear Flows. Springer.
[10]	Monwanou, A.V., Miwadinou, C.H. and Chabi Orou, J.B. (2013) Stability Analysis of Boundary Layer in Poiseuille Flow through a Modified Orr-Sommerfeld Equation. Applied Physics Research, 4, 138-148.[CrossRef]
[11]	Grenier, E., Guo, Y. and Nguyen, T.T. (2014) Spectral Stability of Prandtl Boundary Layers: An Overview. Analysis & PDE, 7, 1045-1075.
[12]	Eckert, M. (2015) Fluid Mechanics in Sommerfeld’s School. Annual Review of Fluid Mechanics, 47, 1-20.[CrossRef]
[13]	Theofilis, V. (2015) Global Linear Instability. Annual Review of Fluid Mechanics, 43, 319-352.[CrossRef]
[14]	Haller, G. (2015) Lagrangian Coherent Structures. Annual Review of Fluid Mechanics, 47, 137-162.[CrossRef]
[15]	Jovanović, M.R. and Bamieh, B. (2005) Componentwise Energy Amplification in Channel Flows. Journal of Fluid Mechanics, 534, 145-183.[CrossRef]
[16]	Cossu, C. and Brandt, L. (2016) On Tollmien-Schlichting Wave Transition in Boundary Layers. European Journal of Mechanics—B/Fluids, 23, 815-833. https://www.sciencedirect.com/science/article/abs/pii/S0997754604000469
[17]	Huang, Y., Ou, W., Chen, M., Lu, Z., Jiang, N., Liu, Y., et al. (2017) Taylor Dispersion in Two-Dimensional Bacterial Turbulence. Physics of Fluids, 29, Article 051901.[CrossRef]
[18]	Huerre, P. (2000) Open Shear Flow Instabilities. In: G.K., Moffatt, H.K. and Worster, M.G., Eds., Perspectives in Fluid Dynamics: A Collective Introduction to Current Research, Cambridge University Press, 159-229.
[19]	Sherwin, S.J. and Blackburn, H.M. (2005) Three-Dimensional Instabilities and Transition of Steady and Pulsatile Axisymmetric Stenotic Flows. Journal of Fluid Mechanics, 533, 297-327.[CrossRef]
[20]	Jovanović, M.R., Schmid, P.J. and Nichols, J.W. (2014) Sparsity-Promoting Dynamic Mode Decomposition. Physics of Fluids, 26, 024103.[CrossRef]
[21]	Grenier, E. and Nguyen, T.T. (2024) On Nonlinear Instability of Prandtl’s Boundary Layers: The Case of Rayleigh’s Stable Shear Flows. Journal de Mathématiques Pures et Appliquées, 184, Article No. 104523.
[22]	Chen, Q., Wu, D. and Zhang, Z. (2023) On the Stability of Shear Flows of Prandtl Type for the Steady Navier-Stokes Equations. Science China Mathematics, 66, 679-722.[CrossRef]
[23]	Reddy, S.C. and Henningson, D.S. (1993) Energy Growth in Viscous Channel Flows. Journal of Fluid Mechanics, 252, 209-238.[CrossRef]
[24]	Pringle, C.C.T. and Kerswell, R.R. (2010) Using Nonlinear Transient Growth to Construct the Minimal Seed for Shear Flow Turbulence. Physical Review Letters, 105, Article 154502.[CrossRef] [PubMed]
[25]	Karp, M. and Hack, M.J.P. (2020) Optimal Suppression of a Separation Bubble in a Laminar Boundary Layer. Journal of Fluid Mechanics, 892, A23.[CrossRef]
[26]	Schmid, P.J. (2011) Application of the Dynamic Mode Decomposition to Experimental Data. Experiments in Fluids, 50, 1123-1130.[CrossRef]
[27]	Blackburn, H.M., Sherwin, S.J. and Barkley, D. (2008) Convective Instability and Transient Growth in Steady and Pulsatile Stenotic Flows. Journal of Fluid Mechanics, 607, 267-277.[CrossRef]
[28]	Huerre, P. and Monkewitz, P.A. (1985) Absolute and Convective Instabilities in Free Shear Layers. Journal of Fluid Mechanics, 159, 151-168.[CrossRef]
[29]	Dhooge, A., Govaerts, W. and Kuznetsov, Y.A. (2003) MATCONT: A MATLAB Package for Numerical Bifurcation Analysis of ODEs. ACM Transactions on Mathematical Software, 29, 141-164.[CrossRef]
[30]	Kuznetsov, Y.A. (1998) Elements of Applied Bifurcation Theory. Springer.
[31]	Govaerts, W.J.F. (2000) Numerical Methods for Bifurcations of Dynamical Equilibria. Society for Industrial and Applied Mathematics.[CrossRef]
[32]	Hart, W.E., Laird, C.D., Watson, J.-P., Woodruff, D.L., Hackebeil, G.A., Nicholson, B.L. and Siirola, J.D. (2017) Pyomo—Optimization Modeling in Python, 2nd Edition, Vol. 67, Springer.
[33]	Wächter, A. and Biegler, L.T. (2006) On the Implementation of an Interior-Point Filter Line-Search Algorithm for Large-Scale Nonlinear Programming. Mathematical Programming, 106, 25-57.[CrossRef]