Application of the Poor Man's Navier-Stokes Equations to Real-Time Control of Fluid Flow

Abstract

Control of fluid flow is an important, underutilized process possessing potential benefits ranging from avoidance of separation and stall on aircraft wings to reduction of friction in oil and gas pipelines to mitigation of noise from wind turbines. But the Navier-Stokes (N.-S.) equations, whose solutions describe such flows, consist of a system of time-dependent, multidimensional, nonlinear partial differential equations (PDEs) which cannot be solved in real time using current computing hardware. The poor man's Navier-Stokes (PMNS) equations comprise a discrete dynamical system that is algebraic—hence, easily (and rapidly) solved—and yet which retains many (possibly all) of the temporal behaviors of the PDE N.-S. system at specific spatial locations. Herein, we outline derivation of these equations and discuss their basic properties. We consider application of these equations to the control problem by adding a control force. We examine the range of behaviors that can be achieved by changing this control force and, in particular, consider controllability of this (nonlinear) system via numerical experiments. Moreover, we observe that the derivation leading to the PMNS equations is very general and may be applied to a wide variety of problems governed by PDEs and (possibly) time-delay ordinary differential equations such as, for example, models of machining processes.

1. Introduction

Given the ever-increasing capabilities of computers and electronic technology, the ability to control a fluid with a range of actuators is already being investigated [1–5] both in laboratory experiments and in numerical simulations. State-of-the-art industrial applications do not exhibit the full potential of control that researchers may envision. Fortunately, it is likely only a matter of time before phenomena such as separation on airfoils, effects of friction factors, and mitigation of turbulence effects are able to be attenuated or amplified to achieve desired outcomes. For example, ability to maintain laminar flow over an airfoil is clearly important, yet not only maintaining, but regaining laminar behavior from a separated flow via real-time control is most desirable. In the petroleum industry, the speed at which oil may be pumped is always subject to frictional losses in the pipelines and pumps. Reducing friction factors through real-time control adds to the general efficiency of such systems. On the other end of the power generation spectrum, control of blade interaction (e.g., pitch and rotation rate) with ambient wind can help reduce chatter and other efficiency losses as well as noise generation in wind farms.

Control schemes may be open loop or closed loop, having different methods of actuation or correction and different algorithms to provide instruction and commands—largely depending on whether open or closed. At present, some methods include microelectromechanical systems (MEMSs) sensors and actuators (as in [2]), local blowing and suction controllers (see, e.g., [4] and [6]), microscale flaps (e.g., [5]), controlled waving plates (as in Techet et al. [7]), and many others. A review of feedback-based methods is provided in Moin and Bewley [8]. It is not our goal to present a comprehensive fluid flow control system herein, or to restrict the study to a particular type of sensor or actuator. Rather, we focus on control of turbulence generically and consider application of forces to a low-order model of the N.-S. equations, the poor man’s Navier-Stokes (PMNS) equations, to accomplish this.

The notion of turbulence control has been researched for roughly the past 20 years, with a substantial mathematical formalism contributed by Abergel and Temam [9] in 1990. Though numerical simulations of turbulence were becoming more feasible at that time, the computing capability of even current modern supercomputers does not allow direct numerical simulations (DNSs) at the smallest of length scales for practical flow geometries and Reynolds numbers. Furthermore, it is well known that analytical solutions of the Navier-Stokes equations are still very much an open problem. Even with the increased processor speed and multicore parallel computing available to researchers today, wall functions and nonphysical subgrid-scale (SGS) quantities are still used to obtain solutions. In order to control turbulent behaviors in real time, speed and simplicity are essential. Furthermore, the amount of (physical) correction and actuation needed in order to produce the desired outcome should be easily determined and provided.

We propose herein the direct modeling of flow velocities through a discrete dynamical system (DDS) with the addition of an adjustable body—or control—force to achieve a desired system response. The behavior of this system has been shown to exhibit many—possibly all—of the temporal behaviors found in the full (PDE) N.-S. system at specific spatial locations. Through a straightforward process, the full N.-S. system is reduced to a coupled discrete dynamical system of equations resembling the much studied logistic map ([10]). It has been shown in 2D (by McDonough and Huang [11]) that the coupled system undergoes the usual bifurcation sequence exhibited by the one-dimensional map, and furthermore that the coupled 2D system contains additional regimes which are not possessed by one-dimensional maps and at the same time retains behaviors found in N.-S. type differential systems. In 3D, the DDS appears to be equally suitable for use as a velocity model (see [12]). The system has been shown to be a practical and realistic model of SGS velocities for large-eddy simulation as well as for the complete velocity. The system of PMNS equations can be fit to local velocimetry measurements of flow over a backward-facing step (see [13]); thus, the system response is sufficiently sophisticated to reproduce laboratory observations. Computationally, the equations are algebraic, and thus very efficiently evaluated, lending themselves naturally for application to control of small-scale turbulence.

Using a DDS to model naturally occurring phenomena is not a new idea—ecologist Robert May constructed the logistic map to model population dynamics [10]. Furthermore, the work of Lorenz [14] led to a “low-order” model consisting of ordinary differential equations (ODEs) via a procedure similar to what is employed herein. The seminal work of Ruelle and Takens [15] viewed the N.-S. equations as a dynamical system capable of generating chaotic, yet deterministic, solutions associated with a “strange attractor.” It has been generally accepted for at least the past decade that the description of turbulence as a deterministically chaotic system is justified (see, e.g., Frisch [16]). Thus, it is from this treatment of the N.-S. equations as a dynamical system that we consider the problem of turbulence control via a closely related DDS.

To implement a control force in the system of PMNS equations, we add a body force to the right-hand side of the N.-S. equations. Since the PMNS equations model local velocities directly, a positive constant control force generally increases the value of the velocity in the respective direction. From studies such as [6], it is apparent that wall-bounded turbulence can be controlled via suction and blowing in the direction normal to the wall. The study herein proposes a numerically analogous approach by varying the value of this additional control force, thus adding to, or taking from, the velocity component magnitude generated by the PMNS equations corresponding to wall-normal directions.

Turbulent flows are known to have relatively high dimension for even low Reynolds number due to the spatial and temporal scales over which they occur (see [8]). Furthermore, without considering variable control forces, the PMNS dynamical system contains nine bifurcation parameters. The codimension of the system is therefore sufficiently high that analytical approaches to studying its behavior quickly become intractable. Though it has already been shown that the PMNS system is capable of chaotic behavior (see [11, 12]), it is worthwhile to mention that chaotic attractors are known to exist within systems with codimension as low as one (e.g., Hopf bifurcations of the N.-S. equations). The system is therefore very complex, and though previous studies have shed light on some of its features, this investigation continues an ongoing effort to explore the gamut of potentially useful characteristics it possesses.

In the following, sections we will outline the derivation of the PMNS equations beginning with the full PDE system of equations. The behavior of the DDS in general is the subject of continuous research, a precis of which will be included herein. We then present regime maps illustrating the types of behaviors that the PMNS equations are capable of reproducing and the corresponding bifurcation parameter and body force values for each regime. The aforementioned system complexity due to high codimension motivates the creation of the regime maps. In doing so, the system behavior is understood throughout the domain of bifurcation parameter space. From these regime maps, we deduce appropriate control forces and then implement these while iterating the PMNS equations to demonstrate controllability of the system toward a desired behavior. Time series of these results are shown.

2. Analysis

Here, we describe the treatment applied to the N.-S. equations by which we derive the PMNS equations. After the derivation, we present a brief discussion of the features that are common to PMNS and full N.-S. equations, specifically the symmetry between equations and highly coupled nature of the full PDE system. Then we discuss the extension of these equations to a control context.

2.1. Derivation of PMNS Equations

We begin with the incompressible N.-S. equations, 𝐔_𝑡+𝐔⋅∇𝐔=−∇𝑃+𝜈Δ𝐔,∀(𝐱,𝑡)∈Ω×𝑡₀,𝑡_𝑓,Ω⊆ℝ³,(2.1) and the divergence-free constraint, ∇⋅𝐔=0,(2.2) so that there are as many equations as necessary to obtain the velocity components and pressure 𝑃. In (2.1), the velocity vector in the 3D domain Ω is given as 𝐔=(𝑢,𝑣,𝑤)^𝑇; 𝜈 is the kinematic viscosity; the Laplace and gradient operators are given as Δ and ∇, respectively, within the appropriate coordinate system. The 𝑡 subscript denotes partial differentiation with respect to time; 𝑡₀ and 𝑡_𝑓 are initial and final times, respectively.

As is frequently used in theoretical analysis of N.-S. equations, we employ a Leray projection to a divergence-free subspace of solutions, thus eliminating the pressure gradient in (2.1). We refer the reader to Constantin and Foias [17] or Foias et al. [18] for details. After the typical scaling of dependent and independent variables and Leray projection, the 3D N.-S. equations take the dimensionless form 𝑢_𝑡+𝑢𝑢_𝑥+𝑣𝑢_𝑦+𝑤𝑢_𝑧=1𝑣ReΔ𝑢,(2.3)_𝑡+𝑢𝑣_𝑥+𝑣𝑣_𝑦+𝑤𝑣_𝑧=1𝑤ReΔ𝑣,(2.4)_𝑡+𝑢𝑤_𝑥+𝑣𝑤_𝑦+𝑤𝑤_𝑧=1ReΔ𝑤,(2.5) which is a 3D system of Burgers’ equations. The 𝑥, 𝑦, and 𝑧 subscripts denote partial differentiation with respect to the spatial variables, and Re is the Reynolds number: Re=𝑈𝐿/𝜈, with 𝑈 and 𝐿 denoting the appropriate velocity and length scales, respectively.

Consistent with the mathematical understanding of the N.-S. equations (see, e.g., [17]), we begin with representing velocity components in Fourier space: 𝑢(𝐱,𝑡)=^∞_𝐤_=−∞𝑣(𝐱,𝑡)=_∞𝑎_𝐤(𝑡)𝜑_𝐤(𝐱),(2.6)_𝐤_=−∞𝑤(𝐱,𝑡)=_∞𝑏_𝐤(𝑡)𝜑_𝐤(𝐱),(2.7)_𝐤_=−∞𝑐_𝐤(𝑡)𝜑_𝐤(𝐱),(2.8) with the 3D wavevector 𝐤≡(𝑘₁,𝑘₂,𝑘₃)^𝑇. We assume that the tensor product basis set {𝜑_𝐤} is complete in 𝐿²(Ω) and orthonormal. Furthermore, we note that construction of such a basis could be difficult for some computational purposes—for example, in constructing proper orthogonal decomposition (POD) algorithms. It is, however, not necessary to do so for the study herein (which distinguishes our approach from POD). For brevity, the following details of the derivation will be provided for only the 𝑥-momentum equation, with treatment of the 𝑦- and 𝑧-momentum equations being identical.

The Galerkin procedure is applied to (2.3), (2.4), (2.5) resulting in an infinite system of ODEs for the Fourier coefficients contained within. Substitution of (2.6), (2.7), and (2.8) into (2.3) leads to _𝐥𝑎̇_𝐥𝜑_𝐥+𝑖_𝐥,𝐦𝑙₁+𝑚₁𝑎_𝐥𝑎_𝐦𝜑_𝐥𝜑_𝐦+𝑖_𝐥,𝐦𝑙₂+𝑚₂𝑎_𝐥𝑏_𝐦𝜑_𝐥𝜑_𝐦+𝑖_𝐥,𝐦𝑙₃+𝑚₃𝑎_𝐥𝑐_𝐦𝜑_𝐥𝜑_𝐦=−1Re_𝐥𝑙²¹+𝑙²²+𝑙²³𝑎_𝐥𝜑_𝐥.(2.9)

Considering the equations without the imaginary factor (which can be removed via methods analogous to those of elementary ODE analysis), inner products of the above equation and each basis function are formed. Employing orthonormality then reduces the 𝑥-momentum equation to 𝑎̇_𝐤+_𝐥,𝐦𝐴^{(1)𝐥𝐦,𝐤}𝑎_𝐥𝑎_𝐦+_𝐥,𝐦𝐴^{(2)𝐥𝐦,𝐤}𝑎_𝐥𝑏_𝐦+_𝐥,𝐦𝐴^{(3)𝐥𝐦,𝐤}𝑎_𝐥𝑐_𝐦=−𝜂^(𝑎)Re||𝐤||²𝑎_𝐤,(2.10) where 𝐤=−∞,…,∞. The factor 𝜂^(𝑎) is a normalization constant arising from the fact that derivatives of basis functions may not possess the same normalization as the functions themselves. Similarly, there are 𝜂^(𝑏) and 𝜂^(𝑐) terms in the omitted 𝑦- and 𝑧-equations, respectively. The 𝐴^{(𝑖)𝐥𝐦,𝐤} factors, with 𝑖=1,2,3, are the Galerkin triple products, for example, 𝐴^{(1)𝐥𝐦,𝐤}=𝑙₁+𝑚₁𝜑_𝐥𝜑_𝐦𝜑_𝐤𝑑𝐱,(2.11) with analogous terms 𝐵^{(𝑖)𝐥𝐦,𝐤} and 𝐶^{(𝑖)𝐥𝐦,𝐤} used in the 𝑦- and 𝑧-equations, respectively.

Removing all but a single arbitrary wavevector in (2.10) results in a single equation, with analogous results for the 𝑦- and 𝑧-momentum equations also shown: 𝑎̇=−𝐴⁽¹⁾𝑎²−𝐴⁽²⁾𝑎𝑏−𝐴⁽³⁾𝑎𝑐−𝜂^(𝑎)Re||𝐤||²𝑏̇=−𝐵₍₁₎𝑎,(2.12)𝑏𝑎−𝐵⁽²⁾𝑏²−𝐵⁽³⁾𝑏𝑐−𝜂^(𝑏)Re||𝐤||²𝑐̇=−𝐶₍₁₎𝑏,(2.13)𝑐𝑎−𝐶⁽²⁾𝑐𝑏−𝐶⁽³⁾𝑐²−𝜂^(𝑐)Re||𝐤||²𝑐,(2.14) with wavevector subscript notation suppressed, and dot (⋅) notation indicating ordinary differentiation with respect to time. The above equations are a system of three nonlinear ODEs. We discretize these using an explicit single-step forward Euler integration method in time and denote the time-step parameter as 𝜏. This leads to 𝑎^𝑛+1=𝑎^𝑛−𝜏𝐴⁽¹⁾(𝑎^𝑛)²+𝐴⁽²⁾𝑎^𝑛𝑏^𝑛+𝐴⁽³⁾𝑎^𝑛𝑐^𝑛+𝜂^(𝑎)Re||𝐤||²𝑎^𝑛𝑏,(2.15)_𝑛+1=𝑏^𝑛−𝜏𝐵⁽¹⁾𝑏^𝑛𝑎^𝑛+𝐵⁽²⁾(𝑏^𝑛)²+𝐵⁽³⁾𝑏^𝑛𝑐^𝑛+𝜂^(𝑏)Re||𝐤||²𝑏^𝑛𝑐,(2.16)_𝑛+1=𝑐^𝑛−𝜏𝐶⁽¹⁾𝑐^𝑛𝑎^𝑛+𝐶⁽²⁾𝑐^𝑛𝑏^𝑛+𝐶⁽³⁾(𝑐^𝑛)²+𝜂^(𝑐)Re||𝐤||²𝑐^𝑛.(2.17) Again for brevity, we omit analysis of the 𝑦- and 𝑧-direction equations and rearrange the 𝑥-momentum equation as 𝑎^𝑛+1=𝜏𝐴⁽¹⁾𝑎^𝑛1−𝜂^(𝑎)/Re||𝐤||²𝜏𝜏𝐴⁽¹⁾−𝑎^𝑛−𝜏𝐴⁽²⁾𝑎^𝑛𝑏^𝑛−𝜏𝐴⁽³⁾𝑎^𝑛𝑐^𝑛.(2.18) Recalling the logistic map 𝑎^𝑛+1=𝛽𝑎^𝑛(1−𝑎^𝑛),(2.19) we note that to recover this form in the above equation requires 1−𝜂^(𝑎)/Re||𝐤||²𝜏𝜏𝐴⁽¹⁾=1,(2.20) which implies 1−𝜂^(𝑎)Re||𝐤||²𝜏=𝜏𝐴⁽¹⁾.(2.21) Having recovered the form of the logistic map via the requirement in (2.20), it is observed that as Re→∞, the left-hand side of (2.21) approaches unity from below. As this term approaches one, the first term on the right-hand side of (2.18) becomes the logistic map, as shown in the following expressions.

Thus, the advanced time step equations can be expressed as 𝑎^𝑛+1=1−𝜂^(𝑎)Re||𝐤||²𝑏𝜏𝑎^𝑛(1−𝑎^𝑛)−𝜏𝐴⁽²⁾𝑎^𝑛𝑏^𝑛−𝜏𝐴⁽³⁾𝑎^𝑛𝑐^𝑛,(2.22)_𝑛+1=1−𝜂^(𝑏)Re||𝐤||²𝑐𝜏𝑏^𝑛(1−𝑏^𝑛)−𝜏𝐵⁽¹⁾𝑏^𝑛𝑎^𝑛−𝜏𝐵⁽³⁾𝑏^𝑛𝑐^𝑛,(2.23)_𝑛+1=1−𝜂^(𝑐)Re||𝐤||²𝜏𝑐^𝑛(1−𝑐^𝑛)−𝜏𝐶⁽¹⁾𝑐^𝑛𝑎^𝑛−𝜏𝐶⁽²⁾𝑐^𝑛𝑏^𝑛.(2.24)

Now define the following bifurcation parameters: 𝛽₁=1−𝜂^(𝑎)Re||𝐤||²𝛽𝜏,(2.25)₂=1−𝜂^(𝑏)Re||𝐤||²𝛽𝜏,(2.26)₃=1−𝜂^(𝑐)Re||𝐤||²𝛾𝜏,(2.27)₁₂𝛾=𝜏𝐴⁽²⁾𝛾₁₃=𝜏𝐴⁽³⁾,(2.28)₂₁𝛾=𝜏𝐵⁽¹⁾𝛾₂₃=𝜏𝐵⁽³⁾,(2.29)₃₁=𝜏𝐶⁽¹⁾𝛾₃₂=𝜏𝐶⁽²⁾.(2.30) Then the equations of motion collapse to the algebraic system 𝑎𝑏^(𝑛+1)=𝛽₁𝑎^(𝑛)1−𝑎^(𝑛)−𝛾₁₂𝑎^(𝑛)𝑏^(𝑛)−𝛾₁₃𝑎^(𝑛)𝑐^(𝑛),(2.31)_(𝑛+1)𝑐=𝛽₂𝑏^(𝑛)1−𝑏^(𝑛)−𝛾₂₁𝑏^(𝑛)𝑎^(𝑛)−𝛾₂₃𝑏^(𝑛)𝑐^(𝑛),(2.32)_(𝑛+1)=𝛽₃𝑐^(𝑛)1−𝑐^(𝑛)−𝛾₃₁𝑐^(𝑛)𝑎^(𝑛)−𝛾₃₂𝑐^(𝑛)𝑏^(𝑛).(2.33) Here, we employ parentheses with the time-level indices to emphasize that these equations comprise an iterated map. Following [11], and in deference to Frisch [16], we term (2.31), (2.32), and (2.33) the “poor man’s Navier-Stokes (PMNS) equations.”

For the study performed herein, we consider the isotropic case, which means that the 𝛾_𝑖𝑗 values are all equal to one another, and the 𝛽_𝑖 values are also equal to one another. Note that if 𝛾_𝑖𝑗=0, for all 𝑖,𝑗, then we recover the logistic map (2.19). Also notice that as Re approaches infinity, the value of 𝛽_𝑖 approaches unity from below. In the DDS of PMNS equations, we limit the starting values between zero and one. Noting that the function 𝑓(𝑎)≡𝑎(1−𝑎) has a maximum range of 1/4 for 𝑎∈[0,1], we rescale the 𝛽_𝑖 values by a factor of 4 (as done by May [10]) to obtain the typical scaling for DDSs. Thus, as in the logistic map, the range of 𝛽_𝑖 values is between zero and four. It is well known that as the bifurcation parameter of the logistic map increases, the system behavior becomes increasingly chaotic. This is in accord with our formulation of the bifurcation parameter 𝛽, which is a function of the Reynolds number. As Re increases, so does the value of 𝛽, and the system behaves correspondingly chaotically.

Note the similar structure of all three components of (2.31), (2.32), and (2.33). This feature is shared with the full N.-S. system of PDEs. Moreover, each equation is nonlinear, and this property is not shared with other often-studied DDSs shown to exhibit chaotic behaviors, for example, discretization of the Lorenz equations (see [14]) or the related 2D DDS due to Hénon [19]. An in-depth investigation of the coupled nature of the PMNS equations is not performed herein, yet may be found in [11, 12]. At present, it is simply important to note that “symmetrical” structure and nonlinearity are desirable qualities preserved in the derivation of the PMNS system from the N.-S. equations.

2.2. Control of PMNS Equations

In order to carry out the intended computational experiments, we must begin with valid ranges for the bifurcation parameters 𝛽_𝑖 and 𝛾_𝑖𝑗. Consideration of isotropy, mentioned above, greatly reduces the space of bifurcation parameter values to consider. Instead of a nine-dimensional space of parameters, we set all 𝛽_𝑖, 𝑖=1,2,3, values equal to one another. The bifurcation parameters 𝛾_𝑖𝑗, 𝑖,𝑗=1,2,3 and 𝑖≠𝑗, are also set equal. From the studies contained in [11, 12], we know which values of bifurcation parameters produce behavior that is of interest in the context of this study of turbulence control and collect data beginning in this subset of bifurcation parameter space.

The following section contains results collected via numerical experiments allowing bifurcation parameters to vary while simultaneously applying a “control force.” The control force is simply a constant added to the right-hand side of the PMNS equations. This is analogous to adding a force at specified wavevectors as is commonly done in DNS and may be viewed as suction or blowing, akin to the study performed by Moin and Bewley [6]. With 𝑎_CF, 𝑏_CF, and 𝑐_CF, the control—or body—forces for each equation, and including the simplification of isotropy, we will solve the following DDS: 𝑎𝑏^(𝑛+1)=𝛽𝑎^(𝑛)1−𝑎^(𝑛)−𝛾𝑎^(𝑛)𝑏^(𝑛)−𝛾𝑎^(𝑛)𝑐^(𝑛)+𝑎_CF,(2.34)_(𝑛+1)𝑐=𝛽𝑏^(𝑛)1−𝑏^(𝑛)−𝛾𝑏^(𝑛)𝑎^(𝑛)−𝛾𝑏^(𝑛)𝑐^(𝑛)+𝑏_CF,(2.35)_(𝑛+1)=𝛽𝑐^(𝑛)1−𝑐^(𝑛)−𝛾𝑐^(𝑛)𝑎^(𝑛)−𝛾𝑐^(𝑛)𝑏^(𝑛)+𝑐_CF.(2.36)

3. Results and Discussion

Here, we present results from iteration of the DDS of (2.34), (2.35), and (2.36) for a range of bifurcation parameter values. We will focus on the isotropic case as already noted above. The bifurcation parameters will vary as 𝛽∈[0,4], with assigned values of 𝛾 (which can typically range from −0.8 to 0.8, with the range depending on values of 𝛽). With varying 𝛽 (related to Re) and fixed 𝛾, we will allow the control force 𝑎_CF to vary and observe the range of behavior types which the system of PMNS equations will produce. Following the analysis techniques in [11, 12], we first describe the types of discernible behavior the system exhibits. Our method of categorizing behaviors—or regime types—is via processing and analysis of the temporal power spectra. Observing the power spectra of the time series output from the DDS and classifying the behavior types into regimes will allow a comparative study of the system behavior, as well as provide insight into the transitions between the regimes as the bifurcation parameter and control force are varied. Having discriminated the regime types, we further explore the ability to control the time series of the DDS by “turning on” the control force to obtain a desired system behavior. Thus, given a specific behavior type corresponding to a pair of bifurcation parameters, we attempt to control the system by means of the control force 𝑎_CF.

The majority of the data presented in this work was calculated on a 376-node Dell high-performance computing cluster at the University of Kentucky Computing Center. Regime map calculations were performed in parallel using OpenMP with three cores from a 2.66 GHz Intel Xeon X5650 six-core processor. Less intensive calculations (i.e., calculating a single time series) were performed on a desktop computer with a 3.60 GHz Intel Pentium 4 and are computed nearly instantaneously. Computations were performed using double precision (64-bit) Fortran. In regime map cases, for each point corresponding to a pair of 𝛽 and 𝑎_CF values, the DDS was iterated 1×10⁴ times. A radix-2 fast-Fourier transform was applied to the last 2¹² (or 4096) iterations to produce power spectra for regime identification. In cases where a single time series was computed, the system was iterated 2×10⁴ times with the control force activated at iteration number 10⁴. Though the computing resources used to conduct these experiments are somewhat powerful, it is important to observe that the PMNS equations may be iterated 𝒪(10⁵) times nearly instantly using a modern laptop computer. Advances in mobile computing hardware and chip technology will soon render the computation of numerous iterations of an algebraic map (such as the PMNS equations) a trivial task for even the smallest and most ubiquitous of devices.

3.1. Regime Maps

As in the studies [11, 12], we present regime maps containing the behavior types—or regimes—which the system of PMNS equations is capable of producing. In the aforementioned investigations, various regime types were distinguished from one another in an automatic and objective fashion. Thus, given a pair of 𝛽 and 𝛾 bifurcation parameters and a value of the control force 𝑎_CF, the system response can be characterized. Furthermore, from those studies, the range of values for which to allow the control force to vary, given a pair of bifurcation parameters, is known via numerical experimentation. The behaviors exhibited by the 3D PMNS equations include all of those behaviors produced by the 1D family of logistic maps and the types of N.-S. regimes found in laboratory experiments. Regime types are distinguished from one another via the power spectral density (PSD), and thus their classification requires human definition of power and frequency criteria for which one regime is distinguished from another. Some discussion of this procedure is appropriate here; however, more detail is included in the two aforementioned studies.

Our goal is to distinguish the various behavior types which the DDS will produce, and do so automatically and objectively. Though some insight may be gained by calculating the fractal dimension, Lyapunov exponents, or other statistical quantities (or by constructing bifurcation diagrams), these will not allow automatic regime classification of the chaotic behaviors of which the DDS is capable. (Some comparison of these methods and their effectiveness is given in [12], and we refer the reader to this citation for more detail, though it suffices to say that the method we present is the most effective for our intents and purposes.) Thus, in order to understand the system behavior as bifurcation parameters—including control forces—are varied, we perform computations across the full range of values with 𝒪(10¹⁰) realizations of the DDS contained within some regime maps. Furthermore, we also compute regime maps over restricted subdomains of bifurcation parameter and control force such that we can clearly study the transitions and localized—in bifurcation parameter space—behaviors. With so many instances of the equations, our means of classifying regime types must be efficient, automatic, and objective.

Identification of the regime type needs to be only qualitative. Thus, the PSD appears to contain sufficient information to permit distinguishing one regime from an inherently different one. The following is a list of the regime types which we have identified in the time series produced by the DDS via this technique: (i) steady, (ii) periodic, (iii) periodic with different fundamental frequency, (iv) subharmonic, (v) phase locked, (vi) quasiperiodic, (vii) noisy subharmonic, (viii) noisy phase locked, (ix) noisy quasiperiodic with fundamental frequency, (x) noisy quasiperiodic without fundamental frequency, (xi) broadband with fundamental frequency, (xii) broadband with different fundamental frequency, (xiii) broadband without fundamental frequency, and (xiv) divergent, as displayed in Figure 1.

Figure 1 

(a) 2501 × 1251 regime map for 𝛾=−0.05 with Δ𝛽=0.0016=Δ𝑎CF, (b) legend of regime types.

Some discussion of the nomenclature and criteria for categorizing the regimes is due. Several of the regimes are referred to as “noisy” along with another description which identifies a distinguishing or primary characteristic of the behavior. This is to imply some broadband features of the PSD along with a salient distinguishing characteristic. Differing from laboratory experiments, the numerical DDS evaluation is of course done entirely within the computer. Thus, there is no “noise” in the sense of the type of signal fluctuation coming from instrumentation or sensors as would be a concern in laboratory experiments. The broadband noise which is observed in the PSDs we consider herein is an actual feature of the DDS. It has been shown in [11, 12] that this noisy behavior exhibits at least a mild sensitivity to initial conditions (SIC) which is known as the “hallmark of chaos” (see, e.g., [20]), thus implying existence of a strange attractor in dissipative systems as considered herein. For further information on the nature of the PMNS equations themselves, the regime types, and the process by which these regimes are distinguished, we defer to the studies in [11, 12].

Figure 1(a) shows the regime map corresponding to a fixed value of 𝛾=−0.05 for 𝛽 ranging from zero to four with the forcing term 𝑎_CF between −0.56 and 1.44. Other values of 𝛾 will also be used, but the results presented herein focus on bifurcation parameter 𝛾 near zero. This is not done for simplification; rather it is a choice which permits a more complete investigation of the behaviors achievable via control (and is typical in LES runs where the 𝛽_𝑖 and 𝛾_𝑖𝑗 are computed rather than assigned). In the studies of [11, 12], no control force terms were used, and the regime maps were computed allowing both 𝛽 and 𝛾 to vary. The types of behavior which the PMNS equations exhibit over the range of 𝛽 values with 𝛾 fixed near zero include nearly all possible regimes, so this choice of 𝛾 results in a thorough investigation.

A key of the colors corresponding to each behavior type is presented in Figure 1(b). There are several features of the regime map worth indicating specifically. Attention should be given to the bifurcation sequence which occurs as the 𝛽 value and the 𝑎_CF value are increased. For most of the map (𝛽<3.3), the sequence corresponds to that of the usual 1D logistic map, (2.19): steady→periodic→subharmonic→chaotic.(3.1) For 𝛽 values larger than three, the system begins to show more erratic behavior with a less clearly defined bifurcation sequence. This region of the map for values of 𝛽>3.3 will be of special interest later as we consider regime maps corresponding to different values of 𝛾. Notice that the regions of the map corresponding to steady and periodic, even including the regions of subharmonic and quasiperiodic, are clearly demarcated. The boundaries between the more noisy and chaotic regime types are less clearly defined and thus complicate the matter of choosing a control parameter value.

After the system bifurcates from subharmonic behavior, as the value of control force and bifurcation parameter are increased, there are several locations where the boundaries between the regimes are not well defined. This result is shared with the 2D and 3D studies of the PMNS equations in [11, 12] and is not unique to the implementation of a control force. The unclear boundary between regimes is what will be referred to as a “fractal region” and is an interesting and important result which is addressed in the two aforementioned investigations. This unclear boundary persists regardless of how finely we are able to resolve bifurcation parameter space, that is, very small steps in the 𝛽 and 𝑎_CF directions. We will provide a high-resolution “zoom-in” of fractal and other interesting regions below.

The same type of behavior is observed in Figures 2(a), 2(b), and 2(c), where the same values of 𝛽 and 𝑎_CF are allowed to vary with 𝛾 fixed at 𝛾=0.0, 𝛾=0.05, and 𝛾=0.1, respectively. We present these regime maps to show the numerous possibilities of controlling the PMNS equations—within a small range of control force values—given a particular pair of bifurcation parameters. To better illustrate this, it is clear from Figure 1, corresponding to 𝛾=−0.05, that for 𝛽>3.3 there are at least eight unique regime types available within a small range of 𝛽 and 𝑎_CF values. The sequence of bifurcations, as 𝛽 and/or 𝑎_CF is increased, is somewhat unstructured. This region for 𝛽>3.3 becomes much more ordered when 𝛾 is held fixed at zero, as shown in Figure 2(a). In this case—as shown in the studies of [11, 12]—the aforementioned bifurcation sequence corresponding to the logistic map (2.19) is sustained across the entire range of 𝛽 values.

(a)

(b)

(c)

Figure 2 

2501×1251 regime maps (a), (b), and (c) for which 𝛾=0.00, 𝛾=0.05, and 𝛾=0.10, respectively, with Δ𝛽=0.0016=Δ𝑎CF.

The ordered sequence of bifurcations presented in Figure 2(a) shows that there are well-defined regions for which the values of control force and bifurcation parameter can be determined and thus used in a control scheme. Given a pair of bifurcation parameters, the control force required to achieve a particular behavior type is immediately known from the regime maps presented herein. As the value of 𝛾 is increased, as in Figures 2(b) and 2(c), a “window” of periodic behavior appears near 𝛽=3.0, and for 𝛽=3.5, there is an interesting bifurcation sequence as 𝛽 is increased: steady→quasiperiodic→subharmonic→phaselock→chaos.(3.2) Figure 2(c), however, shows that the aforementioned region of subharmonic behavior is consumed by the adjacent quasiperiodicity and even some “islands” of phase-lock behavior. The change from 𝛾=0.05 to 𝛾=0.1 shown in Figures 2(b) and 2(c) illustrates the wide variety of achievable behaviors which will be available via changing the control force.

To better exhibit the range of control possibilities within a small range of bifurcation parameter and control force, a localized regime map is presented in Figure 3. It should be noted that this regime map is not simply a visual magnification of the corresponding subregion in Figure 1 for which 𝛽>3.3. Rather, it is produced from a higher-resolution calculation consisting of 1121² points with Δ𝛽=0.000625=Δ𝑎_CF, that is, the bifurcation parameter grid is much finer. Herein, many fractal boundaries exist between the regimes; yet there are still very clear regions of behavior types, and the SIC characteristic of the PMNS equations does not present a problem with respect to determining a suitable control force. Furthermore, there exist ample regions of noisy behavior which may be useful for triggering turbulent activity when desired (e.g., increased thermal advection, avoidance of boundary-layer separation, etc.).

Figure 3 

1121×1121 regime map for 𝛾=−0.05 with Δ𝛽=0.000625=Δ𝑎CF.

There is a tendency for regime types to orient themselves into regions aligned along the direction of varying control force (vertically in the regime maps). When attempting to vary the control force to achieve a different regime type, this presents a problem. In these regions, shown clearly in Figure 3, adjusting the control force 𝑎_CF will lead to a bifurcation in behavior due to the “shape” of the regime in bifurcation parameter space. In addition to this problem, a change of control force will often times lead from one chaotic regime to another, or, worse yet, to a nonphysical result such as divergence. Note that this pattern of regime types oriented along the direction of the control force occurs for ranges of 𝛽 corresponding to chaotic behaviors. That said, the degree to which the regime types are structured in this manner is not consistent for all combinations of bifurcation parameters. The case of Figure 3 considers 𝛾=−0.05. The tendency for regimes to be structured in this manner depends largely on the value of 𝛾 and needs to be studied in further detail. It may in fact be the case that certain values of 𝛾 exhibit different bifurcation sequences, and thus control is limited to certain regions of bifurcation parameter space in which more desirable regime combinations and bifurcation sequences may be found.

Despite this obstacle, there are still regions where the control force may be varied to achieve a particular regime with a priori knowledge of the appropriate control force for a given set of bifurcation parameters. This is shown in the preceding discussion. From (2.25), (2.26), (2.27), (2.28), (2.29), and (2.30), it is clear that the bifurcation parameters contain physical quantities which govern the behavior of the N.-S. equations, that is, the bifurcation parameters correspond to physical flow scenarios. The regime maps presented herein show that given a set of bifurcation parameters—or a physical scenario—employing a control force will alter the system, allowing for a robust control to nearly any desired outcome.

3.2. Control of Time Series

In this subsection, we use contents of the regime map in Figure 1 to help select bifurcation parameters leading to a chaotic behavior. Given the bifurcation parameters corresponding to a particular chaotic regime type, we successfully employ a control force to induce bifurcation to other less chaotic regimes. In the context of applying the PMNS equations to a real-time control system, the bifurcation parameters corresponding to a given scenario can, in principle, be computed via physical data collected from sensors. Thus, the range of control force values appropriate for achieving the desired behavior type could be easily computed via on-board hardware in many specific systems.

Figure 4 shows the time series of the PMNS equations while allowing the value of the control force to change during system evolution. The DDS was iterated 1×10⁴ times prior to applying the control force required to obtain the desired system response. Prior to application of the external force, the methods we employ for detecting the regime type distinguish the behavior as broadband without a fundamental frequency. In starting with this regime type, the PMNS equations produce time series which are erratic and chaotic—essentially “white noise.” Controlling the time series of this most chaotic regime is considered herein, though it should be noted that the PMNS equations can be controlled between essentially any regime types. Thus, the results presented in Figure 4 are only a small fraction of the possible control scenarios that can be achieved via this technique. We focus on controlling the broadband without fundamental frequency regime to show that our procedure is capable of handling even the most chaotic behaviors.

Figure 4 

Time series response of the PMNS equations. Control force 𝑎CF is changed half way through 2×104 iterations with 𝛽=3.3, 𝛾=−0.05, and 𝑎CF=0.17. Control of broadband w/o fund. freq. to (a) quasiperiodicity, (b) phase lock, (c) subharmonic, (d) periodic, and (e) steady, by setting 𝑎CF equal to 0.035, 0.098, 0.06, 0.00, and −0.2, respectively.

The regime map in Figure 1 was constructed with 𝛾=−0.05 and varying values of 𝛽 and 𝑎_CF. Though there are many combinations possible, we choose 𝛽=3.3 and 𝑎_CF=0.17 so that the behavior of the 𝑥-component of velocity is broadband without a fundamental frequency. Many other combinations of 𝛽 and 𝑎_CF would have produced similar results; thus, there is nothing special about these values. Given the known behavior of the PMNS equations using these values of 𝛽 and 𝑎_CF, we then allow the DDS to evolve additional 1×10⁴ iterations after changing the value of 𝑎_CF. The second value of 𝑎_CF also comes from the regime map in Figure 1. For example, to produce a quasiperiodic behavior from the chaotic broadband, the value of 𝑎_CF is changed to 0.035. Again, there are many other values to which the control force can be changed so that the DDS will produce quasiperiodicity from the broadband signal, but a value must be assigned a priori.

The ability to control the 𝑥-direction component of velocity is clearly demonstrated in Figure 4, and the same would also be true of the 𝑦- and 𝑧-direction components. Ability to control a velocity component with wall-normal direction forcing as has been experimentally investigated via blowing and suction (see [4, 6]) and is easily simulated with the PMNS equations, the results of which are presented in Figure 5. Regime maps as constructed in Figures 2 and 3 can be constructed for the case of varying wall-normal control force 𝑏_CF, as performed in previous studies (see [12] for details). From those regime maps, appropriate values of the control force 𝑏_CF are obtained and used to calculate the time series shown in Figure 5.

Figure 5 

Time series response (𝑥-direction) of the PMNS equations with 𝑦-direction control force 𝑏CF is changed half way through 2×104 iterations with 𝛽=2.552, 𝛾=−0.05, and 𝑏CF=0.5872. Control of broadband w/o fund. freq. to (a) quasiperiodicity, (b) phase lock, (c) subharmonic, (d) periodic, and (e) steady, by setting 𝑏CF equal to 0.03888, 0.548, 0.38, 0.25, and 0.12, respectively.

It is clear from both Figures 4 and 5 that the system responds very quickly to a change in control force and does not require substantial time to stabilize. This is in contrast with the initial behavior of the DDS, where many iterations may be required before initial transient behavior decays and the salient behavior of the system can be determined (as explained in [11, 12]). Not only can velocity be controlled within very few iterations, but the efficiency of the PMNS equations allows these sorts of computations to be performed quickly in wall-clock time.

A noticeable difference between control of the PMNS system with streamwise control force (Figure 4) and wall-normal control force (Figure 5) is the increased time required to stabilize in the cases of phase-locked and steady behaviors. Notice that among each behavior type, the response to change in streamwise control force (Figure 4) requires similar times to stabilize. The wall-normal control force case differs from this in that the time required for the phase-locked and steady behaviors to stabilize is slightly longer than in other cases. This difference is quite negligible from the perspective of system control as only few additional iterations are required. That having been said, this observation may allude to more interesting characteristics of the phase lock regime which has been previously studied in [11, 12] and maybe important from the perspective of dynamical systems analysis.

Significant research has been conducted in recent years to control dynamical systems (e.g., [21–23]). The majority of this work, however, does not pertain to systems which model physics, especially fluid motion. The results presented herein are unique in two ways. First, the dynamical system modeled is directly derived from the equations of fluid motion. Second, the control parameter analogously corresponds to a method which is frequently used in laboratory experiments of turbulence control, namely, suction and blowing.

Though the PMNS dynamical system models the motions of fluids, the attempts to control this system may be compared with other studies in which a more general dynamical system is controlled. In [22], a 4D chaotic system is controlled by means of recursive backstepping, and it is shown that chaotic behavior may be eradicated in a similar fashion as in the present study. Control of the system from chaotic to stable states of steady or periodic behavior was exhibited, and the system response times were relatively short. A control problem of a different nature is studied in [23], where a finite-time control scheme is implemented to synchronize two different dynamical systems with chaotic behavior. Effective control is demonstrated as the trajectories of one system converge to those of the other. Furthermore, convergence occurs rapidly once control is implemented.

Results from those general dynamical systems previously discussed have much in common with the time series results presented in Figures 4 and 5. The primary shared features are the ability to control the system to periodic or other stable trajectories. Moreover, rapid system response is observed in each. With the exception of the transition to phase-locked and steady behaviors in the case using a wall-normal control force, the response times herein only span a few iterations. These are expected to be quicker than the results contained in [22, 23], where an algorithm depending on the system output is implemented. No such algorithm is yet developed for the PMNS equations, and the control is achieved with a priori knowledge of the appropriate control force required for a given set of bifurcation parameters from the regime maps shown in Figures 2 and 3. In the following section, we discuss a simple outline for a potential control scheme.

3.3. A Control Algorithm

Here, we present one (of probably many) possible algorithms by means of which the PMNS equations might be used to control turbulent fluid flow. It is important to observe that these equations can be rapidly evaluated, as already emphasized, and their simple algebraic structure is such as to permit easy implementations on microprocessors.

For simplicity, we assume that the desired flow regime is known, although one can easily envision situations where this might not be true. Then to achieve this desired behavior, carry out the following steps. (1)Collect flow data, possibly (but not necessarily) at many locations. (2)Use these data to construct PMNS equation bifurcation parameters. (3)Run the PMNS equations and use the regime map algorithm to identify the current physical flow regime. (4)If regime is the desired one, return to 1. If it is not, begin systematic perturbation of control force parameters with PMNS equation and regime map calculations until desired flow regime is obtained. (5)Convert PMNS equation control force to physical one and send corresponding signal to controller.

We remark that there are a number of details still to be treated in this control procedure, not the least of which is the need to perform implementations on microprocessors. But in this regard, the only part of the algorithm involving more than a few lines of code is the regime map program, and even this is relatively small and should easily fit on modern microprocessors. The more difficult parts of the algorithm are physically collecting sufficient data to build PMNS equation bifurcation parameters with adequate accuracy to be of use, and similarly inverting the process to convert PMNS equation control parameters to physical ones. Clearly, these issues are system dependent; some work is currently in progress to permit investigation of these ideas for a particularly simple system.

4. Summary and Conclusions

In this study, we have derived a discrete dynamical system directly from the 3D, incompressible Navier-Stokes equations and investigated the use of these equations within the context of turbulence control. We have studied the possible behaviors which the PMNS equations may exhibit and included computational results for 𝒪(10¹⁰) instances. Having a priori knowledge of the system behavior, we are able to select bifurcation parameters and control forces such that time series of the PMNS equations are controlled during their evolution. In doing this, we demonstrated the ability of the equations to control a velocity field, and we propose that due to their efficiency, the PMNS equations are well suited for turbulence control.

Derivation of the PMNS equations is, in principle, quite general and can be applied to a wide variety of problems governed by PDEs and (possibly) time-delay ODEs (e.g., models of machining processes). The derivation does not introduce any nonphysical quantities or attempt to model any physical ones. The PMNS equations have been shown to have significant potential as a “synthetic velocity” model, and herein, the ability to manipulate velocity fields across a wide variety flow behaviors was shown. The subject of ongoing investigations and future work will be implementation of the PMNS equations into hardware for use within a real-time control system similar to that described above.

Nomenclature

𝐔=(𝑢,𝑣,𝑤)𝑇:	3D velocity vector
𝜈:	Kinematic viscosity
𝑃:	Pressure
Re:	Reynolds number
𝑡:	Time
Ω:	Spatial domain in ℝ3
𝐤:	Wavevector
𝜂𝑖:	Normalization constant, 𝑖=𝑎,𝑏,𝑐
𝜏:	Time step parameter
𝛽𝑖:	Bifurcation parameter, 𝑖=1,2,3
𝛾𝑖𝑗:	Bifurcation parameter, 𝑖,𝑗=1,2,3 with 𝑖≠𝑗
𝑛:	Iteration counter
𝑎CF:	𝑢-direction forcing term
𝑏CF:	𝑣-direction forcing term
𝑐CF:	𝑤-direction forcing term.