Active Flow Control

Recent Advances in Fundamentals and Applications

Edited by Hui Tang and Xin Wen

mdpi.com/journal/actuators

## **Active Flow Control: Recent Advances in Fundamentals and Applications**

## **Active Flow Control: Recent Advances in Fundamentals and Applications**

Editors

**Hui Tang Xin Wen**

Basel • Beijing • Wuhan • Barcelona • Belgrade • Novi Sad • Cluj • Manchester

*Editors* Hui Tang Department of Mechanical Engineering The Hong Kong Polytechnic University Hong Kong, China

Xin Wen School of Mechanical Engineering Shanghai Jiao Tong University Shanghai, China

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Actuators* (ISSN 2076-0825) (available at: https://www.mdpi.com/journal/actuators/special issues/active flow controls).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

Lastname, A.A.; Lastname, B.B. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-8900-8 (Hbk) ISBN 978-3-0365-8901-5 (PDF) doi.org/10.3390/books978-3-0365-8901-5**

© 2023 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license. The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons Attribution-NonCommercial-NoDerivs (CC BY-NC-ND) license.

## **Contents**


## *Article* **Deep Reinforcement Learning for Flow Control Exploits Different Physics for Increasing Reynolds Number Regimes**

**Pau Varela 1, Pol Suárez 2,3, Francisco Alcántara-Ávila 3, Arnau Miró 2, Jean Rabault 4, Bernat Font 2, Luis Miguel García-Cuevas 1, Oriol Lehmkuhl <sup>2</sup> and Ricardo Vinuesa 3,\***


**Abstract:** The increase in emissions associated with aviation requires deeper research into novel sensing and flow-control strategies to obtain improved aerodynamic performances. In this context, data-driven methods are suitable for exploring new approaches to control the flow and develop more efficient strategies. Deep artificial neural networks (ANNs) used together with reinforcement learning, i.e., deep reinforcement learning (DRL), are receiving more attention due to their capabilities of controlling complex problems in multiple areas. In particular, these techniques have been recently used to solve problems related to flow control. In this work, an ANN trained through a DRL agent, coupled with the numerical solver Alya, is used to perform active flow control. The Tensorforce library was used to apply DRL to the simulated flow. Two-dimensional simulations of the flow around a cylinder were conducted and an active control based on two jets located on the walls of the cylinder was considered. By gathering information from the flow surrounding the cylinder, the ANN agent is able to learn through proximal policy optimization (PPO) effective control strategies for the jets, leading to a significant drag reduction. Furthermore, the agent needs to account for the coupled effects of the friction- and pressure-drag components, as well as the interaction between the two boundary layers on both sides of the cylinder and the wake. In the present work, a Reynolds number range beyond those previously considered was studied and compared with results obtained using classical flow-control methods. Significantly different forms of nature in the control strategies were identified by the DRL as the Reynolds number *Re* increased. On the one hand, for *Re* ≤ 1000, the classical control strategy based on an opposition control relative to the wake oscillation was obtained. On the other hand, for *Re* = 2000, the new strategy consisted of energization of the boundary layers and the separation area, which modulated the flow separation and reduced the drag in a fashion similar to that of the drag crisis, through a high-frequency actuation. A crossapplication of agents was performed for a flow at *Re* = 2000, obtaining similar results in terms of the drag reduction with the agents trained at *Re* = 1000 and 2000. The fact that two different strategies yielded the same performance made us question whether this Reynolds number regime (*Re* = 2000) belongs to a transition towards a nature-different flow, which would only admits a high-frequency actuation strategy to obtain the drag reduction. At the same time, this finding allows for the application of ANNs trained at lower Reynolds numbers, but are comparable in nature, saving computational resources.

**Keywords:** numerical simulation; wake dynamics; flow control; machine learning; deep reinforcement learning

#### **1. Introduction**

In transport applications (especially, in aeronautics), drag reduction is directly related to a decrease in fuel use, which translates into reducing pollution and greenhouse-gas emis-

**Citation:** Varela, P.; Suárez, P.; Alcántara-Ávila, F.; Miró, A.; Rabault, J.; Font, B.; García-Cuevas, L.M.; Lehmkhul, O.; Vinuesa, R. Deep Reinforcement Learning for Flow Control Exploits Different Physics for Increasing Reynolds Number Regimes. *Actuators* **2022**, *11*, 359. https://doi.org/10.3390/ act11120359

Academic Editor: Luigi de Luca

Received: 31 October 2022 Accepted: 29 November 2022 Published: 2 December 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

sions [1]. In the past decades, several techniques have been developed and used to reduce drag, both passively (Bechert and Bartenwerfer [2]) and actively (Gad-el Hak [3]). Passive methods typically rely on fixed geometric changes without using actuators. An example of this technology would be the widespread use of winglets in aircraft. Inspired by bird wings, winglets consist of small wing extensions at the wing tip with an angle relative to the wing-span direction. Using winglets, lift-induced drag is decreased by reducing the size and formation of vortices at the wing tip [4], incurring (hopefully) small increases in the structural weight and parasitic drag. Regarding active methods, diverse techniques are used to reduce aircraft drag and associated emissions. In the work of Tiseira et al. [5] and Serrano et al. [6,7], the use of distributed electric propulsion is combined with boundarylayer ingestion, setting small propellers along the wing near the trailing edge. Thanks to both technologies, it is possible to increase the aerodynamic efficiency of small aircraft or unmanned air vehicles (UAVs). In the same way, the location and use of both jet pumps in wings and synthetic jets are studied to reduce the drag of different aerodynamic bodies. An example of blowing and suction used to control turbulent boundary layers is provided by Kametani and Fukagata [8], and a number of studies have shown the feasibility of this approach in turbulent wings [9–12].

Additional examples can be found in the work of Voevodin et al. [13] and Yousefi and Saleh [14], where the placement of the suction and ejection in wing airfoils is studied to reduce aerodynamic drag or achieve efficient control of the flow around the wing in specific flight-operating conditions. In the same way, the use of synthetic jets has been studied to achieve this same objective, as displayed in the work of Cui et al. [15] on an Ahmed body, or the investigation of Park et al. [16] using an array of synthetic jets applied to a car. A comprehensive review of active flow control in a turbulent flow is provided by Choi et al. [17].

Active-control methods rely on complicated control strategies due to their variable behaviors and dependencies. As shown by Muddada and Patnaik [18], even the control of simple actuators to reduce the drag of a cylinder immersed in a low Reynolds number flow can be complicated due to the interaction between boundary layers, separating shear layers and wake. In this case, by correctly developing a control algorithm, the drag is reduced by about 53 %, leaving ample room for improvement in the case of achieving better control. However, thanks to the application of artificial neural networks (ANNs) and deep reinforcement learning (DRL), it is possible to develop functional control strategies at an affordable computational cost in similar problems, as shown by the work by Rabault et al. [19]. An ANN is a non-parametric tool formed by layers of connected processing nodes or neurons and it can be trained to solve complex problems [20,21]. The ANN training can take place through different types of learning, where DRL is one of the fastest-growing in solving a wide range of cases [22,23], including flow control in complex geometries [24,25].

DRL is based on maximizing a reward by means of an agent interacting with the environment through actions, based on partial observations of the environment. Note that the combination of DRL and ANN was successfully carried out to solve active flow control in the work of Rabault et al. [19]. A proximal policy optimization (PPO) is used to obtain the control policy of two synthetic jets on a two-dimensional cylinder in a low Reynolds number flow. PPO parameterizes a policy function using an ANN with a set of neuron weights given an observation state. The ANN then produces a set of moments that describe a distribution function from which actions are sampled. It is possible to obtain an expression for the estimation of the gradient of the reward as a function of the neuron weights. Additionally, the PPO uses a critic network that estimates the new reward function, which is helpful with stochastic data. Moreover, there is a limit to the maximum update at every training step, preventing rare events from producing large rewards [26]. It was shown that DRL was feasible and enabled finding a control strategy, such that the cylinder drag was reduced by 8 %. Rabault and Kuhnle [27] developed a framework to parallelize environments, enabling speeding up computation and learning, thus reaching

a better solution faster. Additionally, Tang et al. [28] validated the DRL approach in a more complex problem using four jets on the cylinder and extending the Reynolds number range. The results of Rabault et al. in [19,27] have been utilized in a multitude of studies in recent years. For example, the works of Tokarev et al. [29] and Xu et al. [30], where oscillatory rotary controls were applied to the cylinder; the study of linear stability and low sensitivity that allows a better understanding of DRL by Li and Zhang [31], or the effort of Ren et al. [32] where the same DRL approach is used with a solver that allows calculations at higher Reynolds numbers, with the additional problem of controlling the increase in turbulence around the cylinder. Applications to square cylinders, which are relevant to civil engineering, have also been presented recently [**?** ], as well as using traditional modal analysis methods for defining effective reward functions [34]. Regarding these works, there is room for improvement through the application of high-performance computing (HPC), enabling conducting faster calculations at higher Reynolds numbers not previously calculated.

The main objective of this work is, indeed, to demonstrate the capabilities of the HPC resources available today to perform aerodynamic optimization with DRL-based control at much higher Reynolds numbers than the ones already studied in the literature. This project proposes the use of the numerical code Alya, developed in the Barcelona Supercomputer Center (BSC-CNS) [35], employing the HPC MareNostrum IV cluster. Alya is a tier-0 massively parallel code and is designed to solve discretized partial differential equations using finite elements. This code has been successfully used in different problems of fluid mechanics, including simulations of turbulent flows [36,37].

In the present work, the methodology is first detailed in Section 2, dedicated to the numerical setup and the application of the DRL and the ANN. The results are then summarized and discussed in Section 3, which starts with the validation of the Alya solver in the combination of an HPC cluster to solve the DRL problem in the cylinder. Then, the results derived from the increase in the Reynolds number, which exceed the conditions reported in the literature, are shown. Finally, the main conclusions are collected in Section 4.

#### **2. Methods**

This section is divided into two parts: (1) the problem description and the domain setup along the methodology (regarding the numerical simulations) and (2) the framework, including the DRL algorithm.

#### *2.1. Problem Configuration and Numerical Setup*

The domain simulation is two-dimensional (2D) and consists of a cylinder immersed in a flat rectangular channel without any curvature, in the same configuration as that described in the work by Rabault et al. [19]. The domain is normalized using a cylinder with diameter *D* as the reference length scale. The channel has a dimension of *L* = 22 in length, aligned with the main flow direction, and a height equal to *H* = 4.1. The origin of the coordinate system is set on the left-bottom corner of the channel, and the center of the cylinder displaced towards the bottom by 0.05*D* units, to develop the vortex-shedding behind it. A schematic representation of the domain is depicted in Figure 1. Aligned with the vertical axis and symmetrically on the cylinder wall, one at the top and one at the bottom, there are two synthetic jets with total openings of *ω* = 10° each. These positions of the jets were chosen so their actuations were normal to the flow and drag reduction came from effective actuations rather than the injection of momentum.

**Figure 1.** Main domain dimensions in terms of the cylinder diameter *D*, where *ω* represents the jet width and *θ*<sup>0</sup> is the angle of the center of the jet. In light blue, the parabolic velocity profile of the inlet and the jet are represented. The domain representation is not to scale.

Regarding the boundary conditions of the problem, an inlet parabolic velocity condition was imposed and expressed as in Equation (1):

$$\mathcal{U}\_{\rm in}(y) = \mathcal{U}\_{\rm max} \left[ 1 - \left( 2 \frac{y - H/2}{H} \right)^2 \right],\tag{1}$$

where [*U*in(*y*), *V*in(*y*) = 0] is the velocity vector and *U*max is the maximum velocity reached at the middle of the channel, which is equal to 1.5 times the mean velocity *U*, defined as shown in Equation (2):

$$
\Box = \frac{1}{H} \int\_0^H \mathcal{U}\_{\text{in}}(y) \text{d}y = \frac{2}{3} \mathcal{U}\_{\text{max}}.\tag{2}
$$

A value of *U*max = 1.5 is used so that the scaling velocity of the problem is *U* = 1. A no-slip condition is imposed on the cylinder solid wall, while a smooth-wall condition is imposed on the channel walls. The right boundary of the channel is set as a free outlet with zero velocity gradient and constant pressure. The jet velocity (*vi*) is a function of both the set jet angle (*θ*) and the mass flow rate (*Q*) determined by the ANN, as described in Equation (3):

$$v\_i = Q \frac{\pi}{2\omega\sigma R^2} \cos\left[\frac{\pi}{\omega}(\theta - \theta\_0)\right],\tag{3}$$

where *θ*<sup>0</sup> corresponds to the angle where the jet is centered and *R* is the cylinder radius. The scaling factor π/(2*ωR*2) is used so that the integration of the jet velocity over the jet width gives the desired mass flow rate *Q*. More details about the intensity parameters are provided in Section 2.2.

The Reynolds number *Re* = *UD*/*ν* of the simulation is varied between 100, 1000 and 2000, where *ν* is the fluid kinematic viscosity. For the mesh, an unstructured mesh of triangular elements with refinement near the cylinder wall and the jets was used, as shown in Figure 2. The number of elements changes with the Reynolds number as expressed in Table 1.


**Table 1.** Parameters of the simulations for each Reynolds number considered in this work.

(**b**) Detailed mesh around the cylinder.

**Figure 2.** Computational grid used for the *Re* = 100 calculations.

Alya is used to simulate the flow. This solver assumes that the flow is viscous and incompressible, where the governing Navier–Stokes equations can be written for a domain Ω as in Equation (5):

$$
\partial\_t \mathfrak{u} + (\mathfrak{u} \cdot \nabla)\mathfrak{u} - \nabla \cdot (2\nu \mathfrak{e}) + \nabla p = f \qquad \text{in } \Omega \in (t\_0, t\_f), \tag{4}
$$

$$
\nabla \cdot \mathbf{u} = 0 \qquad \text{in } \Omega \in (t\_0, t\_f), \tag{5}
$$

where  is a function of the velocity *u*, which defines the velocity strain-rate tensor (*-* = 1/2(∇*<sup>u</sup>* + ∇*Tu*) ) and *<sup>f</sup>* is the external body force. In Equation (4), the convective form of the nonlinear term, *C*nonc(*u*)=(*u* · ∇)*u*, is expressed as a term conserving energy, momentum and angular momentum [38,39]. This form is known as EMAC (energy-, momentum-, and angular-momentum-conserving equation), and its expression appears in Equation (6). The EMAC form (Equation (6)) adds to the equation the conservation of energy as well as linear and angular momentum at the discrete level:

$$\mathcal{L}\_{\text{emac}} = 2\boldsymbol{\mu} \cdot \boldsymbol{\varepsilon} + (\nabla \cdot \boldsymbol{\mu})\boldsymbol{\mu} - \frac{1}{2}\nabla |\boldsymbol{\mu}|^2. \tag{6}$$

The spatial discretization of the Navier-Stokes equations is performed by means of the finite-element method (FEM). Meanwhile, the time discretization uses a semi-implicit Runge–Kutta scheme of second order for the convective term, and a Crank–Nicolson scheme for the diffusive term [40]. In the time integration, Alya uses an eigenvalue timestep estimation as described by Trias and Lehmkuhl [41]. The complete formulation of the flow solver is described in the work by Lehmkuhl et al. [37].

At each time step, the numerical solution is obtained, and the drag *FD* and lift *FL* forces are integrated over the cylinder surface *S* as follows (Equation (7)):

*F* = (*ς* · *n*) · *e<sup>j</sup>* d*S* , (7)

where *ς* is the Cauchy stress tensor, *n* is the unit vector normal to the cylinder and *e<sup>j</sup>* is a unit vector in the direction of the main flow velocity when calculating the drag and a vector perpendicular to the velocity of flow for the calculation of the lift force. The drag *CD* and lift *CL* coefficients are computed as described in Equations (8) and (9):

$$\mathcal{C}\_{D} = \frac{2F\_{D}}{\rho \overline{\mathcal{U}}^{2} \mathcal{D}},\tag{8}$$

$$C\_L = \frac{2F\_L}{\rho \Pi^2 D}. \tag{9}$$

#### *2.2. DRL Setup*

As discussed in the introduction, the DRL interacts with the domain through three channels. The first channel is the observation state (*s*), based on the extraction of pressure values at a series of predefined points along the domain. These points, known as witness points or probes, are located in the same positions as in Rabault et al. [19]. Moreover, 151 witness points were distributed around the cylinder and along the wake, as shown in Figure 3. The values of the pressure obtained at the witness points are normalized by a factor *s*norm so that the state values given to the agent are between −1 and 1, approximately. The values of *s*norm for each Reynolds number case are given in Table 1.

**Figure 3.** Schematic representation of the computational domain, where the red dots correspond to the location of the probes. The position of two probes is remarked as black dots for further analysis in Section 3.3.

The second channel of interaction between the DRL and the numerical simulation is the action that is given by the control of the jets on the cylinder (*a*). The value of action *a* is directly related to the control value of the upper jet intensity *Q*1, while the bottom jet will do the opposite control to ensure a global zero mass flow rate between both jets, i.e., *Q*<sup>2</sup> = −*Q*1. This is a more realistic control and helps to make the numerical scheme more stable as reflected by Rabault et al. [19]. During the training, the maximum value of |*Q*1| is limited to |*Q*1| < 0.06 *Q*ref ≈ 0.88 for the *Re* = 100 case, as in the work by Rabault et al. [19], to avoid unrealistically large actuations. Note that *Q*ref is the mass flow rate intercepting the cylinder, and it is calculated as in Equation (10):

$$Q\_{\rm ref} = \int\_{-D/2}^{D/2} \rho l I\_{\rm in}(y) \, \mathrm{d}y \, . \tag{10}$$

For higher Reynolds numbers, this clipping value of |*Q*1| is reduced to 0.04. An observation should be made here about the energy consumption of the jets. If one reduces the drag, it would lead to a performance improvement. Nevertheless, the positive energy balance should be obtained when comparing the energy saved by reducing the drag compared with the energy used to actuate with the jets. However, one cannot determine the actual energy consumption of the jets in a numerical simulation, since this highly depends on the systems used on a real setup and, therefore, it is related to experiments and not numerical simulations. An estimation of the energy consumption would consist of a combination of the kinetic energy of the flow injected and the opposing to the local pressure when injecting. Nevertheless, by restricting the mass flow rate of the jets to only 6% of the mass flow intercepting the cylinder, the energy consumption of the jets is guaranteed to be negligible compared with the drag reduction obtained. In addition, even though the value of *Q* selected by the ANN is unique in each action, *Q* is not imposed as a constant value during the entire action duration. In particular, *Q* is assessed as a curve to prevent significant changes in the boundary condition between actions that could lead to numerical discrepancies, similar to how the smooth control was presented by Tang et al. [28]. This way, the imposed mass flow starts from the previous value *Q*<sup>0</sup> and increases or decreases linearly until the new value *Q*<sup>1</sup> is achieved during the entire action time *Ta* = *t*<sup>1</sup> − *t*0. Consequently, *Q* in Equation (3) is presented as a function of time in Equation (11) and is illustrated in Figure 4.

$$Q(t) = \frac{Q\_1 - Q\_0}{T\_a}(t - t\_0) + Q\_0. \tag{11}$$

It should be noted that the first action of each episode (an episode is understood as a sequence of interactions between the neural network and the simulation, which generates the input data for the agent algorithm) always starts from a baseline case, i.e., a periodically stable flow without jet actuation. This baseline case uses the same domain, and the flow is fully developed without applying jet control.

**Figure 4.** Flow rate of the jet *Q* smoothed and applied in the numerical simulation (green) versus the discrete *Q* decided by the DRL agent (red).

The third interaction channel between the DRL and the numerical simulation is the reward or goal (*r*), which in this case is aimed at minimizing the cylinder drag; it is defined in Equation (12):

$$r = r\_{\text{norm}} (-\langle \overline{\mathcal{C}\_D} \rangle - w \langle \overline{\mathcal{C}\_L} \rangle + \mathcal{C}\_{\text{offset}}), \tag{12}$$

where · indicates averaging over a baseline vortex-shedding period, *w* is a lift penalization factor, *C*offset is a coefficient to center the initial reward around 0, obtained from the value of *r* at the end of the baseline simulation, and *r*norm is used to normalize this reward between 0 and 1 approximately. In such a way, the agent receives the reward in an optimal range. Note that the value of *r*norm is set from an a priori guess of the expected maximum reward. The lift penalization factor is set in such a way that the drag is minimized, and at the same time, the effect of the possible growth of induced lift is mitigated. If this lift penalization is not introduced into the reward function, the agent can find a strategy where both jets blow in the same direction at their maximum strength, as discussed in Rabault and Kuhnle [27]. The values of *w*, *C*offset, and *r*norm are given in Table 1 for the different Reynolds numbers.

To summarize, the agent will choose an action (*a*) given a specific state (*s*) in order to maximize a reward (*r*). The function that determines what will be the reward is a normal distribution. During the training process, the agent will choose an action around the average of this normal distribution. This is known as the exploration noise and it helps the method to converge toward a better solution. Once the training has been performed, the agent can be tested in a deterministic mode. Here, the most probable action in the normal distribution is chosen in order to maximize the reward with the learning obtained during the training.

Following the work by Rabault et al. [19], the ANN is designed with two dense layers of 512 neurons and a PPO agent is selected to carry out the control. The PPO agent follows a policy-gradient method in order to obtain the weights of the ANN. The implementation of the DRL is done through the open-source Tensorforce library [42], which is built from the TensorFlow open-source library [43], and it includes defining and creating both the ANN and the control agent. The selection of the initial parameters of the DRL depends intrinsically on the problem itself. In this case, the total number of actions is related to the vortex-shedding period *Tk* = 1/ *fk* through the Strouhal number *St* = *fk* · *D*/*U*.

Figure 5 shows the lift coefficient of the baseline case. This figure shows that the vortex-shedding period is *Tk* = 3.37 time units for the *Re* = 100 case and *St* = 0.29; this is in agreement with both experiments [44] and simulations [19]. Taking the vortex-shedding period as reference, the action time *Ta* is defined as 7.5 % of *Tk*, as in Rabault et al. [19]. This period of actuation was found to be large enough so the consequences of the actuation can be perceived by the flow, and small enough so that the actuations can anticipate and adapt to the needs of the control. In the higher-*Re* cases, this actuation period is reduced to *Ta* = 0.2 due to their more chaotic flow behavior and, therefore, the necessity to adapt the actuations faster. Note that because this is a 2D flow, the word chaotic and not turbulent is

employed, which restricts more to 3D flows. Based on this simulation, it can be concluded that before starting the DRL, the baseline case needs to be run for over 50 time units in order to obtain a periodic stabilized flow. As suggested in the work of Rabault et al. [19], at *Re* = 100 the typical time of an episode should be between 6 and 8 vortex sheddings so that the agent has enough time to learn the suitable control policy. A total of 80 actuations are carried out on each episode, resulting in an episode duration of 20 time units. In the case of a higher Reynolds number, a total of 100 actions are conducted in each episode, since the duration of each actuation is reduced.

**Figure 5.** Temporal evolution of the lift coefficient in the baseline case (without jet actuation). (**a**) Baseline lift-coefficient results. The vortex-shedding period is 3.37 time units. (**b**) Zoom-in view of the vortex-shedding period in the baseline case.

For the high Reynolds numbers of 1000 and 2000, a parallel environment framework was adopted to speed up the learning process using a total of 20 environments. Moreover, in the case at *Re* = 1000, a comparison of the results was done when a single environment was employed. A total of 46 MareNostrum IV central-processing units (CPUs) were used in each environment. Therefore, a total of 46 or 920 CPUs were used depending on whether a single environment or 20 environments were considered in the case, respectively. All of these parameters are presented in Table 1. The general DRL framework is summarized in Figure 6.

**Figure 6.** General overview of the DRL-CFD (computational fluid dynamics) framework employed in this work. A multi-environment approach was used to parallelize the learning and speed up training.

#### **3. Results and Discussion**

This section is divided into three parts. First, the results obtained using the Alya solver and the DRL applied are validated using literature data. Once validated, the higher Reynolds number results are discussed, for which no precedent has been found in the literature. Finally, a cross-application of agents is analyzed to save computational resources in resolving cases with high Reynolds numbers.

#### *3.1. CFD and DRL Code Validation*

As can be seen in the work by Rabault et al. [19], for a *Re* = 100, through the correct actuation of the DRL, it was possible to achieve a decrease in the cylinder drag of 8 %. These results were validated in other simulations, such as Li and Zhang [31], and applied in similar research [29,45,46].

Figure 7 shows eleven different "trainings" launched and the average *CD* obtained for the last vortex-shedding period of the different episodes. It can be seen that all cases follow the same learning trend, where most of the said learning is observed in the first hundred episodes. The main trend is obtained by adjusting a fourth-degree polynomial fit using the average coefficient data, and its purpose is to help to visualize the data. From this point on, the slope of the *CD* decrease is not so pronounced, reaching a stable solution for 350 episodes.

**Figure 7.** Evolution of *CD* in the last vortex-shedding period as a function of the episode.

The obtained learning is then applied to the previously-shown baseline. To this end, one of the 11 trained agents is randomly selected and is run in a deterministic mode. This deterministic simulation is initiated from the conditions of the baseline case at 100 time units from the start. The *CD* and *CL* trends are shown in Figure 8, which indicates that a value of *CD* = 2.95 is obtained after applying the control, i.e., a decrease of 8.9 % compared with the baseline case. The *CD* improvement is calculated as a function of the baseline *CD*:

$$\text{importvement} = \left[1 - \frac{C\_D \text{ controlled}}{C\_D \text{ baseline}}\right] \times 100. \tag{13}$$

**Figure 8.** Temporal evolution of *CD* (**left**) and *CL* (**right**) obtained through the application of the DRL control (at *t* = 100), run in a deterministic mode, for *Re* = 100.

This result is rapidly achieved from the moment the control is applied. In the case of *CL*, the mean remains approximately at 0. However, the amplitude of the variation of this coefficient has been reduced, as shown by the lower standard deviation obtained. This may indeed have great beneficial consequences from the structural and stability points of view.

The control imposed on each jet is represented in Figure 9. As discussed above, the control starts at 100 time units of simulation, so the values of injected or suctioned mass prior to this time are null. Having imposed the synthetic-jet condition, everything injected by one jet (positive values) will be equivalent to what is suctioned by the other (negative values). In the first 10 time units approximately, the control of greater amplitude leads to a significant reduction in the drag, maximizing the reward. Next, the transitional control seeks to obtain a more stable actuation. After 25 time units, the solution is practically periodic.

**Figure 9.** Flow rate through each jet as a function of time after applying the DRL control.

More information about the control is obtained by observing the contours of instantaneous velocity and pressure for the controlled and uncontrolled cases, as shown in Figure 10. The DRL agent reduces *CD* by manipulating the wake vortex: it increases the size of the recirculation region and reduces both the frequency and the amplitude of the von Kármán street downstream of the cylinder. A similar conclusion is obtained from the instantaneous pressure field, which exhibits a decrease in the pressure maxima after applying the control. Comparing our results with those by Rabault et al. [19], it can be seen

that the response of the control is practically the same except for very small differences that can come from the different numerical schemes to resolve the flow.

**Figure 10.** Instantaneous flow fields at *Re* = 100, where for each pair of images, the baseline case without control is depicted on the top and the controlled case is depicted on the bottom. (**a**) Instantaneous magnitude of the velocity fields. (**b**) Instantaneous pressure fields.

Increasing the Reynolds number also requires reducing the period of actuation, as mentioned in Section 2.2. Additionally, increasing the Reynolds number leads to a more complex flow, thus a higher number of episodes are expected to be necessary to complete the training. Therefore, parallelization using 20 environments has been adopted when simulating

higher-*Re* cases. The parallelization using a multi-environment approach is fully detailed in the work of Rabault and Kuhnle [27].

In Figure 11, a comparison between employing a single environment and a 20 multienvironment approach is conducted for a case of *Re* = 1000. It can be observed that if the Reynolds increases an order of magnitude, in this case, the number of episodes for the DRL to act to the same degree is three times greater using only one environment, as compared with Figure 7. The same solution can be obtained by parallelizing through 20 environments, using 50 episodes in each environment, which helps to reduce the real calculation time significantly.

As mentioned before, the flow is more complex, and its structures are less periodic and more chaotic. To analyze the flow, the instantaneous and average velocity fields are represented in Figure 12.

**Figure 11.** Evolution of *CD* in the last vortex-shedding for each episode employing a single environment (**left**) or a 20-multi-environment approach (**right**) showing the learning curve indexed by the episode number for one of the 20 environments.

The obtained results are consistent with those of Ren et al. [32], which report an extensive flow analysis of the same Reynolds number. The control produces an elongation of the recirculation bubble behind the cylinder, reducing at the same time the hydrodynamic drag. For this Reynolds number, the control strategy found by the PPO agent is based on synchronous blowing. This strategy is similar to that found in lower Reynolds numbers: the jets produce an ejection or suction, generating a flow opposite to that caused by the wake.

#### *3.2. DRL Application at Reynolds Number 2000*

Through the use of parallelization with a multi-environment strategy, even higher Reynolds numbers can be achieved. The results obtained for training with *Re* = 2000 are discussed below and, to the authors' knowledge, there is no record of applying DRL control to such a high Reynolds number in the literature.

The *CD* and *CL* obtained through the application of learning in a deterministic mode are represented in Figure 13. This figure contrasts sharply with the one observed with *Re* = 100 (Figure 8). In this case, the starting baseline flow is much more chaotic, which can be seen in the magnitude of the peaks of both variables. The average drag coefficient of the non-controlled flow is *CD* = 3.39. After the control, the average in this period is *CD* = 2.79, which means more than 17 % of *CD* improvement. In this case, it has been chosen to represent the mean in the entire controlled period and not during the last vortex-shedding since, as the solution behaves more erratically, it could lead to an average that is far from reality. Regarding the lift coefficient, the applied control manages to set the average *CL* with a close value to 0, maximizing the reward function as desired and avoiding a systematic biased strategy (see Appendix B in Rabault and Kuhnle [27]). Moreover, the absolute value of the peaks during the oscillating periods is reduced but not as much as in the case of Reynolds number 100, due to the chaotic nature of the flow at this higher Reynolds number.

**Figure 12.** Instantaneous and average flow fields at *Re* = 1000, where for each pair of images, the baseline case without control is depicted on the top and the controlled case is depicted on the bottom. (**a**) Instantaneous magnitude of velocity fields. (**b**) Average magnitude of velocity fields, where the mean streamlines are indicated in black.

**Figure 13.** Temporal evolution of *CD* (**left**) and *CL* (**right**) obtained through the application of the DRL control (at *t* = 50), run in a deterministic mode, for *Re* = 2000.

In Figure 14 the average velocity and pressure fields are shown both for the baseline and the controlled deterministic case. As can be seen in the average velocity field comparison, the separation point in the controlled case is further downstream of the cylinder. This way, the average wake is narrowed faster in the controlled case, lowering the drag produced. This phenomenon is easier to observe if the difference between the baseline and the controlled case is considered as shown in Figure 15. This figure shows that the more significant difference between the controlled and uncontrolled cases occurs in the separation zone. At the same time, as observed for lower-Reynolds number control, the pressure drop behind the cylinder is reduced when controlled.

**Figure 14.** Average flow fields at *Re* = 2000, where each pair of images, the baseline case without control is depicted on the top and the controlled case is depicted on the bottom. (**a**) Average magnitude of velocity fields. (**b**) Average pressure fields.

**Figure 15.** Difference in the mean flow fields at *Re* = 2000, between the controlled and baseline case. (**a**) Difference in the mean magnitude of velocity fields. (**b**) Difference in mean pressure fields.

This control strategy is completely different from that used low at lower Reynolds numbers of values *Re* = 100 and 1000. In order to better understand this control strategy, a chronological velocity and pressure field snapshot is plotted in Figure 16. In this case, the control strategy does not involve the elongation of the recirculation bubble behind the cylinder. The flow separation is controlled and energized by a high-frequency actuation of the jets. The separation point is moved behind the cylinder lowering the drag of the cylinder, similar to the Eiffel paradox phenomenon, also known as the drag crisis, as defined by Stabnikov and Garbaruk [47]. The agent attempts to minimize the drag with this high frequency, breaking the vortices produced after the cylinder into smaller and less-energetic ones.

Additionally, one video (of the velocity and pressure fields) is shared to aid the visualization. The corresponding link can be found at the end of the document in the "Data Availability Statement".

**Figure 16.** Magnitude of the velocity (**left**) and pressure (**right**) field progression at *Re* = 2000. The control starts at *t* = 50, therefore the first two panels of each column correspond to the uncontrolled case and the rest (involving smaller vortices) correspond to the controlled flow.

#### *3.3. Cross-Application of Agents*

Once the *CD* improvement results have been obtained for all the calculated Reynolds numbers thanks to the use of DRL, the cross-application of agents for a flow at *Re* = 2000 is studied.

Cross-application of agents involves applying a previously trained ANN with a different Reynolds number flow to solve a similar problem, using the achieved learning in other conditions. The final objective is to reduce the computational cost of training the DRL agent, since training an ANN with a lower Reynolds number is less expensive. This approach can be very beneficial as long as the results produced by the agent are comparable, i.e., the physics are similar enough between the two different cases.

Here it is investigated the drag reduction obtained by the agents previously trained for Reynolds numbers 100 and 1000 when they are applied in a deterministic mode to the flow of *Re* = 2000, as can be seen in Figure 17.

**Figure 17.** Temporal evolution of *CD* obtained through cross-application of agents trained at *Re* = 100 (**left**) and *Re* = 1000 (**right**) in a deterministic mode to a case with a Reynolds number of 2000.

The cross-application of the agent trained at *Re* = 100 does not yield a result that reduces drag in a flow at *Re* = 2000. No drag improvement is obtained because the nature of the flow where the agent is applied is too different from the one where it was trained. However, when the trained agent at *Re* = 1000 is applied for the case of *Re* = 2000, a significant drag reduction is produced. An average *CD* value of 2.74 is obtained, translating into a 19 % drag reduction improvement compared with the baseline flow. This is a slightly larger drag reduction than that obtained by the agent trained at *Re* = 2000. This small difference can be attributed to the chaotic nature of the flow, from evaluation run to evaluation run, and sub-optimal trade-off with respect to the actual reward function due to less good control of the lift and drag fluctuations, as will be discussed later. The fact that two similar results are obtained with different natures of control strategies, as explained in Section 3.2, may indicate that the flow at *Re* = 2000 belongs to a transition regime toward a higher Reynolds number flow, which only admits a high-frequency control in order to obtain the drag reduction. Therefore, being in this transition regime would admit both controls yielding comparable results. In order to go deep into this topic, higher Reynolds numbers must be simulated, but that goes beyond the scope of this article. At the same time, this is a doubly positive result since it confirms the fact that it is possible to apply deep-learning models previously trained at lower Reynolds number but within the same *Re* regime (similar dynamics) while saving time and computational resources [48].

The average *CD* comparisons for each agent used in the cross-application at *Re* = 2000 are shown in Figure 18 (left). As detailed in the previous section, the control strategy varies significantly as the Reynolds number increases. At *Re* = 100 and 1000, the control attempts to keep the recirculation bubble behind the cylinder, at *Re* = 2000, high-frequency suction and blowing of the jets is applied to quickly force the flow transition, breaking the vortices into smaller ones. Despite these differences in the control strategies between the agents trained at *Re* = 1000 and 2000, both are capable to obtain a comparable drag reduction in a flow at *Re* = 2000. As shown in Figure 18 (right), the *CL* bias obtained using the policy from *Re* = 1000 is higher than that obtained using the policy of *Re* = 2000. Moreover, *CL*

fluctuations are larger when employing the agent trained at *Re* = 1000 (not shown here for the sake of brevity). Note that the *Re* = 2000 strategy reduces the drag almost as much as using cross-application of agents and, at the same time, exhibits a slightly better performance in the lift; since the lift is part of the reward function, the agent trained at *Re* = 2000 achieves a higher reward at *Re* = 2000 than the agent trained at *Re* = 1000. Nevertheless, the small differences between both "trainings" may be reversed for different baseline flows (changing the initial condition). This is similar to the results observed by Ren et al. [32].

**Figure 18.** Average *CD* (**left**) and average *CL* (**right**) comparison between three different crossapplication of agents.

In Figure 19, the average magnitude of velocity field at *Re* = 2000 is shown for the cross-application of the agent trained at *Re* = 1000. At first glance, the existent bias in the average velocity field, which is directly linked with the asymmetric lift behavior, is noticeable.

**Figure 19.** Average flow fields at *Re* = 2000, using the training obtained for *Re* = 1000.

In order to quantitatively observe the different control strategies, the power-spectral density (PSD) of the actions is plotted for the policy of *Re* = 2000, and the cross-application of agents, *Re* = 1000 and 100, in Figure 20 (left). The latter two agents share the same strategy, which consists of a low-frequency action with no significant content at medium and high frequencies. In contrast, with the policy of *Re* = 2000, the frequency spectrum is distributed, with a peak at a frequency of 1.1. In Figure 20 (right), the frequency spectrum of the pressure at the probe in the detachment region indicated in Figure 3 is shown. A low frequency of 0.24 governs the baseline flow, corresponding to the vortex-shedding frequency. This frequency is mitigated by the policy of *Re* = 2000 in higher frequencies. Specifically, the frequency peak at 1.1 is a consequence of the corresponding actuation at this frequency, which is responsible for breaking the flow into smaller vortices as seen in Figure 16. On the other hand, it can be seen how the agent trained at *Re* = 1000 tends towards a completely different strategy in terms of frequency, actuating with a lowfrequency strategy. The agent trained at *Re* = 100, also uses this low-frequency actuation, but it does not have the capabilities to influence the pressure field. As mentioned before, in those cases, the agent attempts an opposition control strategy increasing the recirculation bubble to minimize the drag.

**Figure 20.** Power-spectral density of the action output, *Q*, for the *Re* = 2000 case, with its own policy and two cross-application policies (**left**) and for pressure in the probe near the cylinder in the separation zone (**right**).

#### **4. Conclusions**

In this work, the high-performance CFD solver Alya was coupled with an ANN DRL agent to simulate and control the flow around a 2D cylinder with two jets attached to the cylinder surface. The main control objective was set on the reduction of the average drag. For the first time, the Reynolds number of the canonical 2D cylinder control problem from Rabault et al. [19] was extended to *Re* = 2000, providing new insights into the DRL control strategies for highly complex and dynamic flows. Additionally, the *Re* = 100 and *Re* = 1000 cases were studied for validation and comparison purposes.

The most striking outcome of the *Re* = 2000 case is that the DRL agent uses a radically different strategy from those obtained at lower Reynolds numbers while still being able to provide a 17 % drag reduction. In the new strategy, the agent attempts to delay the detachment point in the cylinder surface using a high-frequency signal in the actuation of the jets, similar to what can be observed in drag crisis phenomena. It is shown that the cylinder wake is narrowed by the breakdown of the detaching vortices into smaller structures. This strategy contrasts with the one obtained at a lower Reynolds number, where the agent acts at a lower frequency to perform opposition control and to elongate the recirculation bubble behind the cylinder. These results are further verified with the spectral analysis of the jet actuating signals and a pressure witness point in the wake. Table 2 summarizes the results from this paper and provides an overview of the results in previous works. Information about drag reduction, the control strategy, and the configuration of the jet location is also provided.

The applications of agents trained at *Re* = 100 and *Re* = 1000 on the *Re* = 2000 case were investigated (namely cross-applications). It was shown that the *Re* = 100 agent is not able to reduce the drag in the high Reynolds number regime due to the different dynamics of the system. On the other hand, the *Re* = 1000 provides satisfactory results, which even beat the *Re* = 2000 agent itself (19 % drag reduction), even though the mean wake displays an asymmetric pattern in this case, which corresponds to an overall reward (including lift bias penalization) that is lower than with the agent trained at *Re* = 2000. Still, this opens up the door to accelerating the training of an agent by first exposing it to a lower Reynolds number flow, which demands a lower computational effort, and then finalize the training at the target Reynolds number condition.

As a final observation, the actual time that the agent takes to calculate the output action once it receives the input state is negligible for an ANN consisting of 2 layers of 512 neurons each, and it could be considered instantaneous. Therefore, it will be completely possible to implement these schemes in real time. In the experiments, there can be some non-negligible delay times depending on the hardware systems used in the experiment. However, if the training is performed in the same conditions, i.e., the delay is also present during the training, the agent will learn a strategy that works under such circumstances. In addition, if the agent is trained in a numerical environment and applies it on the experimental side, then the state sent as an input must consider the known delay in the experiment.


**Table 2.** Summary of the results and strategies. In the strategy, column E states the elongation of the bubble of recirculation and D states the delay of the detachment point through a high-frequency actuation. Note that for the four cases of Tang et al. [28], a unique agent is trained through transfer learning for all four Reynolds numbers and jets have small *x*-velocity components.

**Author Contributions:** Conceptualisation: R.V., J.R. and O.L.; investigation: all authors; data curation: P.V., P.S., F.A.-Á., A.M., J.R. and B.F.; software: J.R., A.M., B.F. and O.L.; writing—original draft: P.V., P.S. and F.A.-Á.; editing: all authors; funding: L.M.G.-C., O.L. and R.V. All authors have read and agreed to the published version of the manuscript.

**Funding:** Pau Varela is partially supported through a grant for the mobility of doctoral students provided by Universitat Politècnica de València and the program Erasmus Prácticas E+ 2020-1. R.V. acknowledges funding from the ERC through grant no. "2021-CoG-101043998, DEEPCONTROL".

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available upon request from the corresponding author. The YouTube link to "Deep reinforcement learning for flow control on a cylinder at Re = 2000" is https://youtu.be/8R3adCQmeEA, accessed on 28 November 2022.

**Acknowledgments:** The authors acknowledge the contribution of Maxence Deferrez to this work. R.V. acknowledges funding from the ERC through grant no. "2021-CoG-101043998, DEEPCONTROL".

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **Abbreviations**

The following abbreviations are used in this manuscript:

#### **Abbreviations**


#### **Roman letters**



#### **Greek letters**


#### **References**


## *Article* **Experimental Study of Dynamical Airfoil and Aerodynamic Prediction**

**Zheyu Shi 1,2, Kaiwen Zhou 3,4,\*, Chen Qin <sup>2</sup> and Xin Wen 3,4,\***


**Abstract:** Dynamic stall is a critical limiting factor for airfoil aerodynamics and a challenging problem for active flow control. In this experimental study, dynamic stall was measured by high-frequency surface pressure tapes and pressure-sensitive paint (PSP). The influence of the oscillation frequency was examined. Dynamic mode decomposition (DMD) with time-delay embedding was proposed to predict the pressure field on the oscillating airfoil based on scattered pressure measurements. DMD with time-delay embedding was able to reconstruct and predict the dynamic stall based on scattered measurements with much higher accuracy than standard DMD. The reconstruction accuracy of this method increased with the number of delay steps, but this also prolonged the computation time. In summary, using the Koopman operator obtained by DMD with time-delay embedding, the future dynamic pressure on an oscillating airfoil can be accurately predicted. This method provides powerful support for active flow control of dynamic stall.

**Keywords:** dynamic mode decomposition; prediction; dynamical airfoil; time-delay embedding

#### **1. Introduction**

Stall is one of the most critical limiting factors for airfoil aerodynamics. Static stall occurs when the airfoil reaches a certain angle of attack. However, in real-world applications, dynamic stall caused by an unsteady motion is more common. An example of such motion is the rapid pitch-up of the airfoil in rotor blades and wind turbines. Compared with static stall, dynamic stall is a more complicated unsteady phenomenon due to the flow separation caused by the rapid change in the angle of attack and the variable stall frequency. Thus, it leads to an abrupt decay of aerodynamic forces and a fluctuation of the structural load on the airfoil. McCroskey et al. [1] comprehensively studied dynamic stall in an oscillating airfoil and dynamic vortex shedding from the airfoil. Further research [2–5] has since provided a detailed understanding of the characteristics of dynamic stall with respect to the lift and pitch moment coefficient. Carr et al. [6] reviewed the significant progress in dynamic stall and concluded that the key influencing factors of this process were oscillation amplitude and frequency, Mach number, and the three-dimensional effects of the airfoil.

Dynamic stall is common in practical engineering and detrimental to the performance of air vehicles. Many researchers employed vortex generators to realize effective control of the dynamic stall [7,8]. Vortex generator control is a passive control method, which is easy-implemented and does not require extra energy injection. However, passive control cannot adapt to the change in working conditions. Recently, many researchers focused on active flow control, which is adjustable based on a feedback control loop [9–12]. Therefore,

**Citation:** Shi, Z.; Zhou, K.; Qin, C.; Wen, X. Experimental Study of Dynamical Airfoil and Aerodynamic Prediction. *Actuators* **2022**, *11*, 46. https://doi.org/10.3390/ act11020046

Academic Editors: André Preumont and Luigi de Luca

Received: 6 December 2021 Accepted: 31 January 2022 Published: 2 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

real-time sensing of the complicated and unsteady flow, e.g., the dynamic stall, is critical in the control loop.

The use of pressure tapes for local measurements is considered a promising method to investigate the flow conditions on an airfoil due to their reliability and ease of implementation. To reduce the complexity of the algorithm or system, some researchers applied some conventional and easy-implemented methods to identify the aerodynamic parameters of the airfoils. Therefore, it is important to predict the flow dynamic based on scattered pressure data. Patel et al. [13] used a pressure sensor array on the upper surface of an airfoil to detect the onset of flow separation. Saini et al. [14] predicted the aerodynamic parameters (i.e., inflow velocity and angle of attack) using discrete surface pressures measured at ports in the vicinity of the leading edge of an airfoil. Provost et al. [15] used a time-shift linear model to predict the roll moment using a sparse set of surface pressure measurements. Juliano et al. [16] found that the pressure field on an oscillating airfoil can have distinct patterns at different locations that are difficult to capture using pressure sensors. Similarly, An et al. [17] develop a method for estimating the instantaneous lift coefficient on a rapidly pitching airfoil that uses a small number of pressure sensors and a measurement of the angle of attack. Gao et al. [18] detected unsteady boundary layer transition on a pitching airfoil using a statistical criterion calculated from thirty surface-pressure transducer measurements. With the development of artificial intelligence, the machine learning method has been successfully applied to sense or predict the aerodynamic parameters [19,20]. However, machine learning methods, such as a blacked box, cannot provide an interpretable physical explanation. Therefore, some researchers employed data-driven techniques to realize the system or parameters identification [21–23]. Zhou et al. [24] explained the physical relationship of pressure data and then developed an offline-online method to sense the real-time aerodynamic parameters. Saini et al. [25] developed the leading-edge flow sensing technique, which uses a few pressures in the airfoil leading-edge to identify and sense various flow events associated with vortex shedding in unsteady airfoils.

Most of the researches mentioned above focused on the flow field around static objects. The associated and available methods typically combine proper orthogonal decomposition (POD) with other algorithms, such as linear stochastic estimation [26] and particle filter [27]. However, it is difficult to predict the evolution of the aerodynamics of dynamic airfoil over time. In this regard, the dynamical system concept should be taken into consideration, and dynamic mode decomposition (DMD) can provide some inspiration, which is an equation-free and data-driven method to model dynamical systems. DMD was initially proposed by Schmid [28] to decompose the flow field into both spatially and temporally coherent structures and has the capability of predicting dynamic systems in the temporal domain. Berger et al. [29] successfully applied DMD to control robotic movement parameters by predicting human–robot interactions based on initial-state input data. However, it is difficult for standard DMD to make predictions based on scattered measurements because the method relies on an infinite-dimensional linear dynamic approximation for a finite-dimensional nonlinear dynamic system. It, therefore, may not provide an accurate result or may even fail when applied to highly nonlinear systems with low-dimensional measurement data [30]. To solve this problem, some data-driven methods have been proposed. The theory of DMD with delay embedding was proposed by Brunton et al. [31] to decompose a three-dimensional chaotic Lorenz system into a linear model in the delay coordinates using a single variable. Mohammad et al. [32] utilized DMD in time delay coordinates to capture the dynamics of a compressible signal containing oscillation perfectly. In addition, Susuki et al. [33] conducted this method to realize the prediction of the voltage dynamics of the rudimentary model directly from the single-bus measurement. The time delay embedding [34] has been investigated in detail on the structure and conditioning of linear models of nonlinear dynamics. Yuan et al. [35] developed DMD with a delay embedding method to predict high temporal and spatial resolved flow fields based on local particle image velocimetry (PIV) measurement. Arbabi et al. [36] applied the delay

embedding method to sparse measurements and realized predictive control of nonlinear flow based on a Koopman model framework.

In the current work, the scattered data obtained from pressure tapes on the surface of an airfoil were used to predict the surface pressure of the airfoil during the oscillation process. Pressure-sensitive paint (PSP) surface measurement provided a qualitative visualization to verify the dynamic stall on the oscillating airfoil. Based on the scattered pressure tape measurements, DMD with time-delay embedding was applied to reconstruct the spatiotemporal pressure distribution on the oscillating airfoil and predict the future data in real-time.

#### **2. Experimental Setup and Methodology**

#### *2.1. Measurement Setup*

The experiment was conducted in an open-circuit low-speed wind tunnel at the China Aerodynamics Research and Development Center. The test section of the wind tunnel has a cross-section of 1.8 m × 1.4 m. In the experiment, the freestream velocity was *U*<sup>∞</sup> = 34 m/s. The two-dimensional airfoil model OA309, made of high-strength carbon fiber composite materials, with a chord length of *c* = 0.4 m and span length of *s* = 1.8 m, was mounted on a pitching testing platform. The freestream Reynolds number based on the chord length *<sup>c</sup>* was Re = 0.9 × 106. The airfoil model was subjected to an unsteady periodic motion. As shown in Figure 1, the pitching motion mechanism was composed of a servo motor, planetary reducer, and encoder. The mechanism allowed control over motion parameters, such as mean angle of attack *α*0, pitching amplitude *α*1, and oscillation frequency *f*. In the experiment, *α*<sup>0</sup> was set as 10◦, *α*<sup>1</sup> was set as 10◦, and *f* was varied as 1, 2, 3, and 4 Hz. The reduced frequency was defined as *k* = *πfc*/*U*∞, with values of 0.037, 0.074, 0.111, and 0.148. The formula for the angle of attack *α* can be expressed as follows:

$$
\alpha = \alpha\_0 + \alpha\_1 \sin(2\pi f t) \tag{1}
$$

**Figure 1.** Schematic of a measurement system using pressure and PSP.

The PSP measurement system and pressure tapes are shown in Figure 1. Part of the upper surface of the airfoil was painted with PSP to evaluate the changes in pressure distribution when the airfoil was oscillating. The PSP was excited by a UV LED (UHP-T-385-EP from Prizmatix Ltd., Givat Shmuel, Israel), and images were captured by a camera (PCO 1600) with a filter. Using a synchronizer (BNC-575) and potentiometer, PSP

images were captured during certain phases of the oscillating airfoil. During the PSP measurement process, exposure signals from the camera were acquired simultaneously with measurements from the potentiometer, which contained the angle of attack of the airfoil. Using these two signals, every PSP frame was assigned to a certain angle of attack. Then, the phase-locked pressure distribution of the upper surface of the airfoil could be obtained by PSP measurement.

The pressure tape measurements and PSP measurements were taken simultaneously. Forty dynamic differential pressure sensors (8510B, ENDVECO) were arranged chordwise at the midspan of the airfoil model. The measurement range of these pressure sensors is 1 psi, and the measurement inaccuracy is ±0.05% of full scale. Locations of these 40 pressure sensors are detailed in Table 1. These pressure sensors adopted a dynamic sampling frequency to sample 256 data points per oscillation cycle of the dynamic airfoil. The acquisition time was approximately 30 s. The aerodynamic parameters, such as lift and pitching moments, were obtained by the integration of the pressure data. Aerodynamic loads can be calculated based on the integration of the surface pressure coefficient. The pressure coefficients can be calculated from the following expression:

$$C\_{p\_i} = \frac{P\_i - P\_{\infty}}{\frac{1}{2}\rho Ul\_{\infty}}\tag{2}$$

where *ρ* is the air density, *Pi* is the static pressure obtained from a single pressure tap, *P*<sup>∞</sup> is the static pressure of the incoming flow. The detailed process of the integration can refer to previous material [37].

**Table 1.** Locations of the pressure transducers on the OA209 airfoil.


#### *2.2. Dynamic Mode Decomposition with Time-Delay Embedding*

DMD is an efficient method of decomposing a dynamical system into specific spatial modes that evolve at a certain frequency and growth–decay rate in time. In DMD, the measurement data obtained by sensors at time *j* are arranged into a column vector *xj* = (*x*1,*j*, *x*2,*j*,... , *xn*,*j*), where *xj* is a data snapshot and *n* denotes the number of measurement points (in this case, the number of sensors). All *m* snapshots were sampled with a constant time interval (sampling time step) Δ*t*. DMD describes a discrete-time system as:

$$\mathbf{x}\_{\mathbf{j}\star\mathbf{1}} = A\mathbf{x}\_{\mathbf{j}} \tag{3}$$

where matrix *A* reflects the dynamical characteristics of the system. The first step of DMD is to arrange all *m* snapshots into two matrices, *X***<sup>1</sup>** = (*x***1**, *x***2**, ... , *xm−***1**) and *X***<sup>2</sup>** = (*x***2**, *x***3**,... , *xm*). Then, the linear approximation *A* may be written in terms of matrices based on Equation (4) as:

$$X\_2 \approx AX\_1\tag{4}$$

Using matrix *A*, known as the Koopman operator, future data can be predicted from current data. Thus, the core of DMD is to search for the optimal linear approximation *A* of the nonlinear dynamical system. The crucial step of the DMD algorithm is obtaining the rank-reduced representation of *A* in terms of a POD-projected matrix *A* , which can be written as follows:

$$\tilde{A} = \mathcal{U}^T A \mathcal{U} = \mathcal{U}^T \mathcal{X}\_2 V \Sigma^{-1} \tag{5}$$

where *<sup>U</sup>* ∈ **<sup>R</sup>***n*×*<sup>r</sup>* , <sup>Σ</sup> ∈ **<sup>R</sup>***r*×*<sup>r</sup>* , and *<sup>V</sup>* ∈ **<sup>R</sup>***m*×*<sup>r</sup>* are obtained by singular value decomposition (SVD). Here, *r* is defined as the truncated number and represents the rank of approximation to *X***<sup>1</sup>** using SVD reduction. The next step is to compute the eigen-decomposition of *A* , which is similar to *A*:

$$AW = W\Lambda\tag{6}$$

where the columns of *W* are eigenvectors of *A* and Λ is a diagonal matrix containing the eigenvalues *λk*, which capture the time dynamics of the corresponding DMD modes. These eigenvalues and eigenvectors of *A* can be related back to the similarity-transformed eigenvalues and eigenvectors of *A* to reconstruct the DMD modes given by:

$$
\Phi = X\_2 V \Sigma^{-1} W \tag{7}
$$

where Φ is the matrix whose columns each represent a DMD mode. Then, a snapshot at any time can be reconstructed from the DMD modes and eigenvalues:

$$\mathbf{x}\_{\mathbf{j}} = \Phi \Sigma^{j-1} \mathbf{b} \tag{8}$$

where *b* = *Φ−***<sup>1</sup>** *<sup>×</sup>* **<sup>1</sup>** contains the initial values of DMD modes. As shown above, the DMD method relies on the linear approximation of a nonlinear dynamic system. When the subdomains are of limited spatial dimensionality, it is difficult for DMD to provide accurate predictions. In this study, the direct use of discrete sensor measurements as the observables creates a system that is far removed from the infinite-dimensional linear approximation of highly nonlinear systems. To solve this problem, extended DMD [38] and kernel DMD [39] have previously been proposed. These methods use new observables to provide a set of coordinates in which the dynamical systems appear linear. However, a universal and systematic approach needs to be found. The details of time delay embedding methods can refer to the works of Brunton et al. [31]. Time-delay embedding augments the limited spatial observables by embedding future measurements into the current measurements. Therefore, more information can be captured to construct the linear Koopman operator.

A requirement of time-delay embedding is to have access to at least *nd* + 1 consecutive snapshots, where *nd* is the number of delay steps. The time-delay-embedded matrices *X***<sup>1</sup>** and *X***<sup>2</sup>** can be transformed into *H***<sup>1</sup>** and *H***<sup>2</sup>** as follows:

$$H\_1 = \begin{bmatrix} \mathbf{x\_1} & \mathbf{x\_2} & \cdots & \mathbf{x\_{m-n\_d}} \\ \mathbf{x\_2} & \mathbf{x\_3} & \cdots & \mathbf{x\_{m-n\_d+1}} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{x\_{Hd}} & \mathbf{x\_{Hd+1}} & \cdots & \mathbf{x\_{m-1}} \end{bmatrix} = \begin{bmatrix} y\_1, y\_2 \cdots, y\_{m-n\_d} \end{bmatrix} \tag{9}$$

$$H\_2 = \begin{bmatrix} \mathbf{x}\_2 & \mathbf{x}\_3 & \cdots & \mathbf{x}\_{m-n\_d+1} \\ \mathbf{x}\_3 & \mathbf{x}\_4 & \cdots & \mathbf{x}\_{m-n\_d+2} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{x}\_{n\_d+1} & \mathbf{x}\_{n\_d+2} & \cdots & \mathbf{x}\_m \end{bmatrix} = \begin{bmatrix} y\_2, y\_3 \cdots \ y\_{m-n\_d+1} \end{bmatrix},\tag{10}$$

where *yj* is a new snapshot after delay embedding, and *xm* occupies the last *n* rows *ym*−*nd*+1.

The dimensions of the matrices after delay embedding *H***<sup>1</sup>** and *H***2**, are (*n* × *nd*) × (*m* − *nd*). When *nd* = 1, the matrix after delay embedding is the same as the original matrix. Analogously to Equation (4), the linear mapping of *H***<sup>1</sup>** and *H***<sup>2</sup>** can be written as follows:

$$H\_2 \approx BH\_1\tag{11}$$

Then, the time-delay-embedded matrices *H***<sup>1</sup>** and *H***<sup>2</sup>** are used to build the DMD model. The matrix *B* is the Koopman operator and reflects the dynamical characteristics. Therefore, standard DMD can be applied to find the corresponding eigenvalues and eigenvectors of *B*. Using the matrix *B* as follows:

$$y\_{j+1} \approx By\_j \tag{12}$$

The current new snapshot *yj* at time *nd* + *j* − 1 can be used to predict the future new snapshot *yj*+1, which contains the data at time *nd* + *j*.

#### **3. Results and Discussion**

*3.1. Dynamic Stall*

#### 3.1.1. Flow Visualization

To examine the two-dimensional pressure field on the oscillating airfoil, the images from PSP measurement were phase-averaged over 15 periods. The evolution of the pressure field during the oscillation of the clean airfoil is shown in Figure 2 at two typical AoAs under an oscillation frequency of *f* = 1 Hz. It should be noted that these PSP measurements only provided a qualitative visualization because the PSP data was influenced by the temperature of the incoming flow.

The upper panels represent the downstroke, and the lower panels represent the upstroke. The low-pressure region is indicated in blue, and high pressures are shown in yellow. At *α* = 9◦, the pressure distribution in the downstroke is similar to that in the upstroke. However, when *α* increases to 13◦, an obvious low-pressure region forms at the leading edge of the airfoil in the upstroke, as can be seen in Figure 2d. The color intensity shows that the pressure at the leading edge in the downstroke is higher than that in the upstroke at the same angle of attack, which suggests the onset of flow separation and the decrease in lift. The phenomenon that flow separation occurs at the leading edge in the downstroke but does not occur in the upstroke is consistent with the dynamic stall was also investigated by Heine et al. using PIV measurement [8].

#### 3.1.2. Aerodynamic Loads

The angle of attack of the airfoil was varied to examine the influence of the VGs mounted at the leading edge of the airfoil on the static aerodynamic loads. In our study, the aerodynamic loads of static and dynamical airfoil were measured to verify the phenomenon of dynamic stall. The measurement results of static airfoil are shown in Figure 3. Note that there is no oscillation frequency for this static case. As shown in Figure 3, the curves of the lift and pitch moment show a relatively constant trend until *α* ≈ 13.5◦. With a further increase in the angle of attack, the lift and pitch moment drop abruptly. An increase of 0.5◦ in the angle of attack leads to a drop of approximately 0.39 in the lift coefficient and 0.08 in the pitch moment coefficient. Therefore, *α* = 13.5◦ can be considered as the stall angle with respect to flow separation. These results were in line with our common perception and are consistent with the findings of previous studies [8,40].

The influence of the oscillation frequency *f* of the dynamical airfoil on the aerodynamic loads was examined in this study. The lift and moment coefficients under different AoA for the oscillating airfoil at different *f* are shown in Figure 4. It should be reminded that the relevant aerodynamic load data were phase-averaged from 10 measurement periods, and the resolution was improved by interpolation. In the present work, instantaneous pressure data at different periods are similar, so phase-averaged data based on 10 periods is credible.

**Figure 2.** PSP visualization of the dynamic pressure on the oscillating airfoil with *f* = 1 Hz. (The horizontal direction is the full chord length, and the vertical direction is the spanwise direction).

**Figure 3.** Static stall characteristics are shown in coefficients of (**a**) lift (*Cl*) and (**b**) pitch moment (*Cm*) for the static airfoil.

**Figure 4.** Phase-averaged lift coefficient (*Cl*) at frequencies of (**a**) *f* = 1 Hz, (**b**) *f* = 2 Hz, (**c**) *f* = 3 Hz, and (**d**) *f* = 4 Hz for oscillating airfoils.

When *f* = 1 Hz, the lift coefficient of the upstroke is higher than that of the downstroke by up to 0.8 under the same angle of attack. There was a crucial characteristic feature in the curve of loop hysteresis, which was related to a certain *α* = 6◦ when *f* = 1 Hz. When α < 6◦, the lift coefficients of the upstroke is lower than that of the downstroke under the same AoA. Once α > 6◦, the lift coefficients of the upstroke is higher than that of the downstroke. Therefore, we can say that there are two hysteresis loops during the oscillating process. It can be assumed that the pressure of the lower surface of airfoil can be higher under small AoA due to the downstroke of the airfoil. When increasing the oscillation frequency to 2 Hz (Figure 4b), the certain AoA which decided the two-loop hysteresis became smaller. Continuing to improve the oscillation frequency *f* to 3 Hz (Figure 4c), only one loop hysteresis is left. However, when α= 0, 1, 2◦, the lift moments under upstroke and downstroke were close. When the oscillation frequency *f* was increased to 4 Hz (Figure 4d), the lift coefficients under even a small angle of attack during the downstroke were smaller than that during the upstroke. With increasing oscillation frequency *f*, the loop hysteresis becames more pronounced, and the shapes of the hysteretic loops changed. In detail, the increase of oscillation frequency *f* induced the slight increase of upstroke and downstroke in the lift coefficient under larger AoAs, but the lift coefficient under smaller AoAs in the downstroke decreased slightly.

In addition to the lift coefficients, the characteristics of pitch moment were also investigated, as shown in Figure 5. For *f* = 1 Hz, there is no abrupt drop in the moment coefficient such as the static case during the upstroke. The moment coefficient of both upstroke and downstroke was very close, i.e., the difference of moment coefficient between upstroke and downstroke was small. It is worth noting that compared to the static case, the moment coefficient under large AoAs decreased dramatically, which suggested that dynamic stall influenced the moment coefficient. For example, when α = 20◦, the moment coefficient in

the static case was about −0.055 but it decreased to about −0.085 in the dynamical case. With increasing oscillation frequency *f*, the curve of moment coefficient formed two-loop hysteresis. The moment coefficient of upstroke and downstroke showed a huge difference, suggesting that high oscillation frequency leads to a large fluctuation in pitch moment. Combined with the change in lift coefficient, it can be inferred that not only did the dynamic stall exert a negative influence on the aerodynamic load, the oscillation frequency *f* also did. More importantly, high oscillation frequency *f* makes the aerodynamic performance more intricate (i.e., two-loop hysteresis in a coefficient curve). This problem is detrimental to the performance of real-time active flow control. Therefore, real-time perception, even the prediction of aerodynamic parameters, is challenging and meaningful. In Section 3.2, we will use a novel data-driven method to reconstruct the time-resolved pressure data and then realize the prediction of future pressure.

**Figure 5.** Phase-averaged moment coefficient (*Cm*) for (**a**) *f* = 1 Hz, (**b**) *f* = 2 Hz, (**c**) *f* = 3 Hz, and (**d**) *f* = 4 Hz for the oscillating airfoil.

#### *3.2. Prediction of Future Data*

In this section, DMD with time delay embedding was used to reconstruct the history of pressure data of the oscillating airfoil. The high time-resolved pressure data collected by 20 sensors on the upper surface (see Table 1) were used in this method because the leading edge is the region with the most obvious pressure fluctuations. To examine the performance of this method, the error rate *e* between reconstructed data and true data is defined as:

$$\varepsilon = \frac{1}{m} \sum\_{t=1}^{m} \frac{||P\_{reconstructed,t} - P\_{true,t}||\_2}{||P\_{true,t}||\_2} \tag{13}$$

where *m* is the number of used snapshots, *Preconstructed*,*<sup>t</sup>* is the reconstructed pressure field at time t, and *Ptrue*,*t* is the true velocity field on the subdomain at time *t*.

POD analysis was applied to the time delay embedded dataset, and the mode energy spectrum is shown in Figure 6. The core of time delay embedding is the improvement of data dimensionality to capture more phase information. Therefore, the first three POD modes contained more than 98% energy. With the increase of the number of delay steps *nd*, more POD modes were needed to reconstruct the data accurately. It should be noted that when the number of delay steps was increased to 200, the POD energy spectrum tended to be similar. In this experiment, one period contains 256 snapshots. Therefore, when *nd* > 200, most of the phase information was included in the time delay embedded matrix.

**Figure 6.** POD mode energy spectrum of the time-delay embedded data with different *nd*.

The raw and reconstructed pressure data of sensor #1 were used to compare the reconstruction performance of the standard DMD and DMD with delay embedding. Sensor #1 is located at the leading edge, where the most obvious pressure fluctuations happen. The standard DMD (i.e., *nd* = 0) and DMD with time-delay embedding were then applied to reconstruct the whole-time sequence of the pressure data. As shown in Figure 7a, standard DMD cannot accurately predict the dynamic pressure on sensor #1 even at the initial stage of the oscillation cycle. With time-delay embedding, the prediction performance is clearly improved. With *nd* = 50 delay steps, DMD can capture the periodic trend, but with an underprediction of the oscillation amplitude and e ≈ 43%. With *nd* increased to 200, the predicted values agree closely with the ground truth both in periodic trend and amplitude for more than six cycles and e ≈ 7%. Thus, DMD with time-delay embedding demonstrates the ability to reconstruct the periodic single-point pressure data over a long period. With an increasing number of delay steps *nd*, the accuracy of the reconstruction is improved. However, an excessive number of delay steps will lead to too large a dimensionality of the matrix, which will not only consume computing and storage resources but negatively affect the real-time performance of the prediction algorithm. Therefore, the duration of the prediction cycle should be balanced with the requirement for computation time. In fact, the accurate prediction of just one cycle, that is, the next cycle, is usually helpful for effective decision-making in real-time control systems.

Based on the predicted pressures at single measurement points, DMD with time-delay embedding was then used to reconstruct the time-resolved pressure field on the upper surface of the airfoil. The ground truth of the pressure field obtained by 20 pressure sensors is shown in Figure 8a. The values reconstructed using DMD with time-delay embedding are shown in Figure 8b and match well with the measured data and accurately show the shape and peak of the three-dimensional surface.

**Figure 7.** Reconstruction of pressure data measured by sensor #1 (see Table 1) using (**a**) standard DMD and time-delay embedding DMD with number of delay steps (**b**) *nd* = 50, (**c**) *nd* = 100, and (**d**) *nd* = 200.

**Figure 8.** Pressure field on the upper surface of the oscillating airfoil: (**a**) ground-truth data from pressure sensors, (**b**) reconstruction using DMD with time-delay embedding.

The core of DMD is to find the Koopman operator, which allows the prediction of future data based on past data. The above section has demonstrated the accurate reconstruction of time-resolved pressure data using DMD with time-delay embedding, from which an optimal Koopman operator was obtained. Using the Koopman operator, real-time prediction can be realized. In Figure 9, the black line represents the ground-truth pressure data measured by sensor #1 over one oscillation period. The blue markers represent the past sampled pressure data that were used to build a new snapshot yj . The future new snapshot *yj*+<sup>1</sup> was obtained based on the Koopman operator *B* provided by DMD with time-delay embedding, which had been verified in the above cases. Then, the future pressure data at the position of sensor #1 could be extracted from *yj*+1. Furthermore, future snapshot *yj*+<sup>2</sup>

can also be calculated with the help of a predicted snapshot *yj*+<sup>1</sup> in the same way. When this operation is repeated, the data for a long period of time in the future can be predicted. In Figure 9, 256 snapshots (a whole oscillation period) were predicted. The predicted next-step value accurately matches the ground-truth value, indicating the feasibility and potential of DMD with time delay embedding. The computation was implemented in MATLAB, and all the cases are run on an Intel Xeon E5-1620 at 3.50 GHz. Every prediction took 0.09 s. Admittedly, it is difficult to meet real-time requirements, but high-performance industrial computers can realize real-time prediction in practical applications. Work is ongoing to further improve the rapidity and adaptability of this method.

**Figure 9.** Prediction of future pressure data based on Koopman operator obtained by DMD with time-delay embedding.

#### **4. Conclusions**

In this study, dynamic stall of oscillating airfoils was measured by dynamic pressure tapes and PSP visualization. The influence of the oscillation frequency was examined. Dynamic mode decomposition (DMD) with time-delay embedding was proposed to predict the pressure field on the oscillating airfoil based on scattered measurements. The major findings are as follows:

(1) Under the dynamic conditions of the oscillating airfoil, with increasing oscillation frequency, although lift coefficient can be improved at large AoAs, the aerodynamic characteristics showed a huge difference between upstroke and downstroke. The specific performance is that the loop hysteresis of both the lift and pitch moment coefficients were enlarged, which is detrimental to aerodynamic performance.

(2) DMD with time-delay embedding was able to reconstruct and predict dynamic stall based on scattered measurements with much higher accuracy than standard DMD. The reconstruction accuracy of this method increased with the number of delay steps, but this also prolonged the computation time. With the help of the Koopman operator obtained by DMD with time-delay embedding, the future dynamic pressure for a long time of an oscillating airfoil can be accurately predicted.

**Author Contributions:** Writing—original draft, Z.S.; Methodology and Investigation, K.Z.; Resources, C.Q.; Writing—review & editing, X.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China (grant number: 12072196), the Science and Technology Commission of Shanghai Municipality, China (grant number: 19JC1412900), the Oceanic Interdisciplinary Program of Shanghai Jiao Tong University (grant number: SL2020MS006), Advanced Jet Propulsion Innovation Center, China (grant number: AEAC HKCX2020-02-028), Aeronautical Science Foundation of China (grant number: 2020Z006057001) and State Key Laboratory of Aerodynamics, China (grant number: SKLA-20200303).

**Acknowledgments:** All individuals included in this section have consented to the acknowledgement. The authors gratefully thank for the support of Lingrui Jiao and Di Peng from Shanghai Jiao Tong University, who provided some help and suggestion for the pressure-sensitive paint measurement.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclature**


#### **References**


## *Article* **Experimental and Numerical Study on Incident Shock Wave/Boundary Layer Interaction Control**

**Chuanbiao Zhang, Yanhao Luo, Hua Liang \*, Shanguang Guo and Hesen Yang**

Science and Technology on Plasma Dynamics Laboratory, Air Force Engineering University, Xi'an 710038, China; zhangcb2020@126.com (C.Z.); lyh\_1415@163.com (Y.L.); shanguang\_guo@163.com (S.G.); yanghesen96@126.com (H.Y.)

**\*** Correspondence: lianghua82702@163.com

**Abstract:** This study was designed to explore the control effect of pulsed arc discharge plasma actuation on the incident shock wave/boundary layer interaction (ISWBLI). Research was conducted on an ISWBLI flow field with 10 kHz single-channel pulsed arc discharge plasma actuation and pulsed arc discharge plasma actuation array applied at Mach 2.0 experimentally and numerically. In the investigation, high-speed schlieren flow field visualization technology was adopted, focusing on the change in shock wave intensity caused by plasma actuation. Combined with the detached eddy simulation (DES) method, the numerical simulation focused on the regulating effect of plasma actuation on the separation zone. The key research results showed that, in terms of the spatial flow field, the simulation results were consistent with the experimental results. The single-channel actuation could only just achieve the control effect on the ISWBLI, while array actuation effectively weakened the shock wave intensity. Furthermore, the ISWBLI separation zone in the base flow field was crescent shaped. Its length at the middle of the flat plate was longer than that at the two sides. It was identified that, after applying single-channel actuation, the start of separation slightly moved forward. Similarly, after the application of array actuation, the start point of separation at the middle section in a spanwise direction moved forward by about 19 mm. The length of the separation zone increased by 30 mm but reduced at the two sides. Its influence, spanwise, was also significantly diminished.

**Keywords:** shock wave/boundary layer interaction; pulsed arc discharge plasma actuation; flow control; high-speed schlieren flow visualization; DES

### **1. Introduction**

The supersonic inlet is an important part of an aircraft's air intake system. When an aircraft is flying at supersonic speed, the inlet acts to slow down the airflow and provide subsonic airflow for the engine, ensuring the airflow is even and stable [1]. This process is of great importance as the quality of airflow directly affects the flight performance [2]. During the flight, the multiple compression configuration to slow down the airflow can significantly reduce the total pressure loss [3], but the ISWBLI triggered by compression is inevitable at the inlet lip. Moreover, the unsteadiness of the interaction region brings extra structural load to the aircraft [4,5], which threatens flight safety. In addition, with the adjustment of the flight attitude, the operating conditions of the inlet may change extremely. If that is the case, the ISWBLI can thicken the boundary layer and cause flow separation, resulting in total pressure loss and inlet unstart [6–8].

The ultimate approach to solve the ISWBLI issue is to address the separation of the boundary layer caused by the reverse pressure gradient. Therefore, flow control technologies based on boundary layer control are commonly used in the inlet [9]. In recent years, significant achievements have been made in solving the ISWBLI issue, regardless of their passive or active forms [10–12]. Various types of vortex generator (VG) were studied at

**Citation:** Zhang, C.; Luo, Y.; Liang, H.; Guo, S.; Yang, H. Experimental and Numerical Study on Incident Shock Wave/Boundary Layer Interaction Control. *Actuators* **2022**, *11*, 148. https://doi.org/10.3390/ act11060148

Academic Editor: Luigi de Luca

Received: 22 April 2022 Accepted: 26 May 2022 Published: 2 June 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

different incoming Mach numbers, configurations, and layouts, revealing groundbreaking outcomes [13–15]. Huang et al. [15] discovered that a VG can achieve better results in suppressing flow separation. Similarly, other studies [16–19] concluded that a micro-vortex generator (MVG), with a size of only a fraction of the thickness of the boundary layer, can achieve similar control effects to a VG. They found that an MVG can also reduce the drag of the aircraft during high-speed flights [20]. Furthermore, modern scholars focused on the control of the ISWBLI with MVGs through experiments and numerical simulation. Yan et al. [16] compared the ISWBLI with and without MVG control at Mach 2.5. In their studies, the MVG reduced the separation zone and improved the shape factor. Their experimental results also showed that vortex rings play an important role in the reduction of the separation zone. Mohd et al. [17] found that the presence of micro-ramps reduces the upstream interaction length by delaying the pressure rise, hence, suppressing the SBLI's effect. Their research proved that MVGs can also be efficiently applied in hypersonic ISWBLI control. However, with the continuous optimization of the aerodynamic shape of aircrafts, it is likely that the size of the inlet will be constantly reduced to reduce the drag and overall weight of the flight [21]. Besides, a flow control device with a convex shape could easily detach and damage the engine blades. For all of these reasons, higher requirements are imposed on flow control devices.

The plasma actuator is a kind of small-sized active flow control device without moving parts. It is easy to maintain and, therefore, meets the application needs in the inlet [22]. Plasma actuators have become a research hotspot in the field of flow control due to their advantageous features such as fast reactions [23,24], wide bandwidth, and tunable actuation intensity. The common forms of plasma actuator are dielectric barrier discharge (DBD), plasma synthetic jet (PSJ), and pulsed arc discharge (PAD). PSJs exhibit the unique ability to produce high-velocity (>300 m/s) pulsed jets at high frequency (>5 kHz); thus, they are widely used in high-speed flow control [25]. Venkat et al. [26] studied the control effect of a PSJ on the ISWBLI at Mach 3.0 and found that a PSJ can shift the dominant frequencies of the separated flow. Venkateswaran et al. [27] found that the position of the pulsed plasma jet can influence the control effect on the ISWBLI, suggesting that a more satisfactory effect is achieved when the actuator is located upstream of the separation zone rather than inside the separation zone. Kursat et al. [28] numerically tested the performance of a sweeping jet actuator on flow separation control. The simulation results suggested that the performance of the sweeping jet actuator, in terms of decreasing the size of the separation bubble, increases with the mass flow rate. Through numerical simulation, Kiyoshi et al. [29] studied two forms of DBD plasma actuators. They found that the electrode, in parallel with the incoming direction, increases the thermal effect of the plasma and exacerbates flow separation, but the form of the canted electrode generates a vortex and mitigates flow separation. Through experiments and simulations, Yan et al. [30] found that the thermal effect of arc plasma actuation can reduce the Mach number upstream of the shock wave. The shock wave intensity is also weakened. Upon the study of the streamwise array of pulsed spark discharge at different frequencies, Tang et al. [31] found that increased actuation frequency helps weaken shock intensity. Although studies on SWBLI control based on plasma actuators are widely available, many scholarly conclusions have focused on the configuration of the compression ramp. There is a research gap when it comes to the control of the ISWBLI based on plasma.

This paper focuses on the control effect of plasma actuation on shock intensity and the separation region. First, the control effect of the pulsed arc discharge plasma actuation on ISWBLI is verified by wind tunnel experiments. The difference in the control effect on the ISWBLI with single-channel and array pulsed arc discharge plasma actuation is studied by numerical simulation. Additionally, the spanwise control effect of plasma actuation on the separation region is also analyzed. The results will provide a reference for effective ISWBLI control in supersonic inlets.

#### **2. Experimental Setup**

#### *2.1. Wind Tunnel and Experimental Model*

The experiment was carried out in the supersonic wind tunnel in the Science and Technology on Plasma Dynamics Laboratory, China. The wind tunnel is a blowdown with a total temperature of 300 K, total pressure of 95.6 kPa, and unit Reynolds number of 1.19 × <sup>10</sup>7/m. It can provide Mach 2.0 supersonic flow for 2~3 s. The experimental section of the wind tunnel is equipped with a fixed bracket for model installation, allowing the position and angle of the model to be adjusted. There are two optical observation windows with a diameter of 250 mm on the side of the tunnel, providing light channels for observing the flow field. The installation of the experimental model in the wind tunnel is shown in Figure 1. The model consisted of a bottom flat plate and a top shock wave generator. The shock generator had a wedge angle of 15◦ at the front edge and was installed on the fixed bracket. It protruded downward through a strut and hung at the nozzle of the wind tunnel. Its lower surface formed an angle of 12.9◦ with the horizontal plane. The flat plate was made of acrylic plastic, and the upper surface was horizontal and remained at a spacing of 109.8 mm from the shock wave generator.

**Figure 1.** Schematic diagram of installation of experimental model in the wind tunnel.

As shown in Figure 2a, the plate model was 400 mm in length and 110 mm in width. The model was divided into three sections according to its functions, i.e., leading-edge section, actuation section, and interaction section. The first section was the zone from the leading edge of the plate to the most upstream plasma actuator. It was 70 mm in length and mainly used for adjusting the development of the boundary layer to control the downstream boundary layer state. The second section was 270 mm in length and was located in the downstream zone from the plasma actuator. The ISWBLI occurred in this zone by adjusting the angle of the shock generator. The coupling of plasma actuation with the interaction flow field also occurred in this area. The actuation section was the zone between the aforementioned two sections, and the pulsed arc discharge plasma actuator was equipped at this section. With reference to the relatively mature control mechanism developed previously through compression corner SWBLI control using pulsed arc discharge plasma, two kinds of actuator were set in the actuation section, i.e., a singlechannel pulsed arc discharge plasma actuator and a pulsed arc discharge plasma actuator array. The single-channel actuator was placed 70 mm away from the leading edge of the plate, and the distance between the positive and negative electrodes was 5 mm. Five

single-channel actuators formed an actuator array, and each kept a spacing of 15 mm in the flow direction. The connection of the pulsed arc discharge plasma actuator array is illustrated in Figure 2b. The positive electrode of the first actuator was connected to the output end of the plasma power supply while the negative electrode was connected to the positive electrode of the next actuator. The five actuators were connected in series to form a pulsed arc discharge plasma actuator array.

**Figure 2.** Schematic diagram of experimental model. (**a**) plate model (**b**) The connection of the pulsed arc discharge plasma actuator array.

#### *2.2. Plasma Actuation System and High-Speed Schlieren System*

As shown in Figure 3, the triggering of the pulsed arc discharge plasma actuation and the collection of electrical parameters were dependent on the plasma actuation system. Additionally, the high-speed schlieren system was used for the capture of the ISWBLI flow field. As mentioned above, the positive electrode of the array actuator was connected to the output end of the plasma power supply. The plasma power supply selected for this paper was a nanosecond pulse power supply that could provide a pulsed, high voltage ranging from 1 Hz to 100 kHz. The parameters of voltage and current were transmitted to the oscilloscope through the probes and then collected by the oscilloscope. The high-speed schlieren system was a typical, Z-type, optical path schlieren. Reflected by two concave mirrors, the beam was partially blocked at the knife edge and then received by a CCD camera. This process allowed the flow field information carried by the beam to be properly captured. The nanosecond pulse power supply was set with the output voltage as 20 kV, frequency as 10 kHz, and pulse width as 1 μs. The rising and falling edges were both 50 ns. The time interval between two consecutive arc plasma actuations was 100 μs. Additionally, the frame rate of the CCD camera was 25,000 fps; hence, the frame interval was 40 μs.

**Figure 3.** Schematic diagram of plasma actuation system and high-speed schlieren system.

#### **3. Experimental Results and Discussion**

As the schlieren image reflected the spanwise integral effect of the spatial flow field, the control effect of the single-channel pulsed arc discharge plasma actuation on the ISWBLI flow field could not be effectively refined in the schlieren image; it will not be listed herein. More flow field information is analyzed and discussed in the numerical simulation section. The experimental results of controlling the ISWBLI with array pulsed arc discharge plasma array actuation are shown in Figure 4. Figure 4a represents the instantaneous schlieren image of the base flow field. As illustrated, a strong incident shock wave was generated at the leading edge of the shock generator, and the incoming flow boundary layer was relatively stable. The interaction between the incident shock wave and the boundary layer induced a typical ISWBLI flow field structure; affected by the reverse pressure gradient after the shock wave, the boundary layer near the interaction point was thickened and lifted up, resulting in local flow field compression and a separation shock induced at the leading edge of this region. Then, the boundary layer was reattached to form local flow field compression at the wall surface, after which another shock wave appeared, namely, reattached shock.

**Figure 4.** Experimental schlieren results: (**a**,**c**,**e**) the base flow field; (**b**,**d**,**f**) the actuation flow field.

Figure 4b is the instantaneous schlieren image with 10 kHz array actuation. As shown in the figure, the high-frequency actuation induced a series of precursor blast waves in the flow field, and the entire region of the shock wave was covered with the precursor blast waves. Being constantly impacted by the continuous precursor blast waves, the structure of the separation shock wave and the reattached shock wave were furcated. Compared with the base flow field, the reattached shock wave almost disappeared, and the shock wave intensity was considerably weakened. On the other hand, the actuation created a thermal gas bulb. When passing through the interaction region, the thermal gas bulb was coupled with the shock wave and the boundary layer to make the structure divergent, showing that the thermal gas bulb plays a positive role in promoting flow mixing and increasing the kinetic energy of the boundary layer. Although the zone covered by the precursor blast waves included the foot of the incident shock wave, no significant change was observed in this part of the incident shock wave.

The structure of an ISWBLI flow field can be better displayed through time-averaged results; hence, 1000 schlieren images of the base flow field and the actuation flow field were averaged, respectively, as shown in Figure 4c,d. In the time-averaged results, after the actuation was applied, the reattached shock wave almost completely disappeared, and the gray value of the separation shock also decreased. Furthermore, the intensity of the two shock waves was weakened, and the size of the interaction region decreased, indicating that the thermal gas bulb successfully inhibited flow separation. As represented in Figure 4e,f, a flow field with or without actuation can be compared more intuitively according to the threshold value of the spatial gradient of schlieren intensity. In the base flow field, the separation shock wave and the reattached shock wave were shown as two discontinuous strips. The strips disappeared after actuation, indicating that the shock wave intensity was weakened.

#### **4. Numerical Results and Discussion**

In the previous section, the role of single-channel actuation and array actuation for ISWBLI was experimentally studied, and the results show that the high frequency pulsed arc plasma actuation array can effectively weaken the intensity of the separation shock wave. In this section, the action of pulsed arc plasma actuation in controlling the ISWBLI is further analyzed through numerical simulation. The analysis focuses on the influence of actuation on the separation zone in the spanwise direction. The ultimate goal is to obtain a comprehensive understanding of pulsed arc plasma actuation in controlling the ISWBLI.

#### *4.1. Calculation Method*

#### 4.1.1. Physical Model

The ISWBLI model studied in this section consisted of a shock wave generator and a flat plate, as shown in Figure 5. The length of the plate model was *L*s = 400 mm. The shock generator was located above the plate, and its height from the surface of the plate was *H*<sup>s</sup> = 110 mm. The width of the shock wave generator was *W*<sup>s</sup> = 110 mm, corresponding with that of the plate. The angle was also consistent with the shock wave generator model studied in the experiment. It is worth specifying that the origin of the coordinates was located at the center of the plate's leading edge. Furthermore, the positive direction of the *x*-axis was the incoming flow direction; the positive direction of the *y*-axis was the normal direction of the wall surface; and the *z*-axis was the spanwise direction of the plate.

**Figure 5.** Computational model of ISWBLI: (**a**) 3D view; (**b**) front view; (**c**) top view.

#### 4.1.2. Computational Grid and DES

ICEM was used to divide the computational grid into structured grids. At the same time, dense grids were located in the ISWBLI region, separation shock wave region, and plasma actuation region to obtain a more refined flow field structure. Three sets of grids were obtained. The number of the basic grids (Mesh\_1) was 4.21 million. On this basis, two sets of denser meshes (Mesh\_2, Mesh\_3) were obtained, and the number of grids was 6.27 million and 7.63 million, respectively. Comparing the skin friction coefficients calculated by three sets of grids, as shown in Figure 3, Mesh\_2 and Mesh\_3 were in good agreement; hence, Mesh\_2 (6.27 million) was chosen for further simulations, as shown in Figure 6. The height of the first layer of the skin was set to 0.007 mm to ensure y+ < 1.

**Figure 6.** Verification of mesh independence.

The spanwise width of the flow field entrance was 5.275 *L*s; the normal height was 5 *L*s; and the distance from the leading edge of the plate was 3.75 *L*s. The far-field boundary conditions of pressure were set based on the experimental conditions. As shown in Figure 7, the total pressure was 1.01 × <sup>10</sup><sup>5</sup> Pa; the total temperature was 300 K; and the incoming Mach number was Ma = 2.0. Adiabatic non-slip skin boundary conditions were adopted for the plate and shock wave generator.

**Figure 7.** Model grid and boundary conditions.

The DES method was used to calculate the three-dimensional ISWBLI flow field and the characteristics of the flow field under the control of plasma actuation. Figure 8 shows the comparison of the numerical schlieren and the experimental schlieren. Figure 8a shows the flow field structure of the experimental mean schlieren. As mentioned in the previous section, near the foot of the incident shock wave, the boundary layer thickened, and a separation shock wave appeared. As the boundary layer reattached to the wall, a weak reattachment shock wave appeared. It can be seen that the angle of the incident shock wave was about 41.2◦, that of the separation shock wave was about 38.4◦, and that of the reattachment shock wave was about 33.6◦. Upstream of the interaction region, the boundary layer was very thin and was a stable laminar boundary layer. The thickness of the boundary layer increased after passing through the interaction region. Figure 8b shows the simulation results calculated with the DES method. The figure clarifies that the angle of the incident shock wave was about 40.8◦. The difference between the simulation and the experiment was about 0.1% of the experimental results. The separation shock wave was slightly smaller, about 37◦, which is 3.6% less than the experimental results, and the reattachment shock wave was about 34.7◦, which is 3.3% higher than that in Figure 8a. Furthermore, the structure of the interaction region was well simulated in the numerical schlieren, and the boundary layer thickened, passing through the interaction region, which is consistent with the experimental results. Based on the analysis above, it can be concluded that the simulation result was credible.

#### 4.1.3. Phenomenological Model of Plasma Actuation

In-house and foreign studies showed that the control effect of arc plasma actuation over the flow field is mainly a thermal effect. The phenomenological model of plasma actuation is a standard structure used to study plasma simulation. Differing from other physical models, this one simplified plasma flow control into the injection of kinetic and thermal energy into the flow field due to plasma discharge, which considerably reduced the amount of calculation, and can be successfully coupled with the NS equation. Therefore, in this study, the arc plasma actuation was simplified as a heat source model, and it was added to the flow field in the form of a heat source term. The energy of a single actuation was set to 10 mJ, and the density of the discharge heating power was 1 × 1011 W/m3. This value was a rough estimate based on several wind tunnel experiments. The efficiency of the pulsed arc actuation in this paper was approximately 67%. It was the total thermal efficiency of actuators, including instantaneous rapid heating of gas, material heating of models, heat dissipation, and other thermal energies. The shape of the plasma heat source model was set as a square with a side length of 5 mm due to the fact that the spacing of the arc plasma actuation in the experiment was 5 mm.

In order to verify the accuracy of the plasma heat source model, the simulation results were compared with the experimental results obtained in the supersonic wind tunnel, as shown in Figure 9. To make the comparison more rigorous, the experimental five-channel array actuation was selected, as shown in Figure 9a. The precursor blast wave and thermal gas bulb produced by each actuation channel had similar structures, which showed that the experimental results were reliable. Figure 9b reveals that the calculated precursor blast wave and thermal gas bulb were consistent with the experiment. The simulation of the arc plasma actuation based on the proposed heat source model could achieve accurate results.

**Figure 9.** Plasma actuation evolution: (**a**) experimental schlieren; (**b**) numerical schlieren.

#### *4.2. Simulation Result and Disscussion*

#### 4.2.1. Base Flow Field

Figure 10 shows the density gradient contours of different spanwise cross sections. At the spanwise middle cross section (*z* = 0), the structure of the separation shock wave, reattached shock wave, and the separation zone could be clearly identified in the density gradient contour. Furthermore, the separation zone still existed at *z* = 0.025 m. The strength of the shock wave was lower than that at the spanwise middle cross section. At the cross section at the edge of the plate (*z* = 0.054 m), the reattached shock wave and the separation zone disappeared, and the ISWBLI was a regular reflection because of the strong turbulence in the boundary layer at the edge of the plate. Based on the flow field structures at the three cross sections, it can be concluded that the ISWBLI demonstrates significant, spanwise, three-dimensional characteristics.

**Figure 10.** Density gradient contour of different cross sections: (**a**) *z* = 0; (**b**) *z* = 0.025 m; (**c**) *z* = 0.054 m.

The velocity contour in the ISWBLI zone was further analyzed. Figure 11 shows the velocity contour at the y = 0 cross section. Velocity u represents the velocity in the *x*-axis direction of the flow field, and the region where Velocity u is less than 0 is defined as the separation zone, which is showed in blue in Figure 11. The calculation results show that the separation zone induced by the incident oblique shock wave was crescent shaped, and the length of the separation zone reached its peak at the central section of the plate and decreased spanwise. The separation zone had a width smaller than that of the shock wave generator and the plate and basically disappeared at the edge of the plate. These outcomes demonstrate that the ISWBLI flow field has significant, spanwise, three-dimensional characteristics, which is consistent with the conclusion extracted from Figure 10.

**Figure 11.** Velocity contour of ISWBLI at *y* = 0 cross section.

The specific location of the separation zone could be determined based on the skin friction coefficient Cf. Figure 12 represents the change curve of *C*<sup>f</sup> at the *z* = 0 cross section and, as can be seen, at the central cross section of the plate, Cf had a negative value from x = 0.126 m to x = 0.146 m in the flow direction. This indicates that the length of the separation zone at the central cross section of the plate was 20 mm.

**Figure 12.** Skin friction coefficient curve of base flow field at *z* = 0 cross section.

Figure 13a,b shows the skin friction coefficient curves at the *z* = 0.025 m cross section and *z* = 0.054 m cross section, respectively. It can be interpreted from the figure that, at *z* = 0.025 m, the separation zone was about 16 mm in length from x = 0.0136 m to x = 0.0152 m. At *z* = 0.054 m cross section, the separation zone was located at x = 0.146 m~x = 0.153 m with a length of 7 mm. The separation moved backward with a smaller size. This result shows that the size of the separation zone at the edge of the plate was smaller than that at the middle of the plate.

**Figure 13.** Skin friction coefficient curve of base flow field at different spanwise cross sections: (**a**) *z* = 0.025 m; (**b**) *z* = 0.054 m.

#### 4.2.2. Single-Channel Actuation Flow Field

In this section, the coupling of single-channel arc plasma actuation and the ISWBLI flow field was studied. The single-channel actuation was set at the central cross section of the plate, 70 mm from its leading edge. The actuation frequency was set to 10 kHz, and the actuation energy was 10 mJ. Figure 14 shows the pressure contour of the ISWBLI flow field controlled by plasma actuation at different moments. The start time of discharge was recorded as *t* = 0. When *t* = 0.5 μs, a region with increased pressure appeared at the actuation position, because a local, high-pressure region was induced by high temperature due to the discharge. The high-pressure region was a square, similar to the shape of the plasma source term. When *t* = 11 μs, the precursor blast wave structure induced by actuation appeared in the flow field, which was essentially generated due to the outward transfer of local high pressure. When *t* = 21 μs~*t* = 41 μs, the precursor blast wave continued to spread. It can be interpreted from the figure that the precursor blast wave weakened as it spread and almost entirely disappeared when *t* = 51 μs.

**Figure 14.** Pressure contour at *z* = 0 cross section at different moments.

Subsequently, the evolution of the thermal gas bulb generated by the arc plasma actuation was analyzed through the temperature contour of the flow field. As shown in Figure 15a, when *t* = 0.5 μs, a local, high-temperature zone appeared at the actuation position with a maximum temperature of 715 K. A thermal gas bulb was generated near the surface due to discharge heating. Then, the thermal gas bulb propagated downstream, accompanied by energy loss and structural diffusion. As shown in Figure 14c,d, since the structure of the thermal gas bulb was slightly higher than the boundary layer, the moving speed of the part far from the skin was higher while that of the part close to the skin was lower. This difference happened due to the mainstream shear force. When *t* = 100.5 μs, the second actuation was already triggered. At this time, the thermal gas bulb generated by the previous discharge propagated to the ISWBLI zone. As shown in Figure 15e, affected by the thermal gas bulb, the boundary layer was thickened. It can be concluded that the disturbance caused by the actuation was not only due to the shock, but also the thermal gas bulb.

**Figure 15.** Temperature contour at *z* = 0 cross section at different moments.

Figure 16 shows the numerical schlieren time-averaged flow field with and without the single-channel actuation flow field at different spanwise cross sections. At the three selected cross sections (*z* = 0, *z* = 0.025 m, *z* = 0.054 m), no obvious changes in the shock wave structure and the separation zone were identified after the actuation was applied.

**Figure 16.** Numerical schlieren of three cross sections: (**a**) *z* = 0 base flow field; (**b**) *z* = 0 actuation flow field; (**c**) *z* = 0.025 m base flow field; (**d**) *z* = 0.025 m actuation flow field; (**e**) *z* = 0.054 m base flow field; (**f**) *z* = 0.054 m actuation field.

Figure 17 shows the curves of the skin friction coefficient at different spanwise cross sections. At three selected cross sections (*z* = 0, *z* = 0.025 m, *z* = 0.054 m), the curves of the skin friction coefficient of the base flow field and the actuation flow field were basically the same. At the *z* = 0 cross section, the length of the separation zone slightly increased in the actuation flow field, and the start of separation moves forward moderately. At the cross section away from the center of plate (*z* = 0.025 m, *z* = 0.054 m), the skin friction coefficient curved of the base flow field, and the actuation flow field completely overlapped.

**Figure 17.** Skin friction coefficient curve at different spanwise cross sections: (**a**) *z* = 0; (**b**) *z* = 0.025 m; (**c**) *z* = 0.054 m.

The form of the separation zone with and without actuation was also taken into account. Figure 18 illustrates the velocity contour of the *y* = 0 cross section. Compared with the crescent-shaped separation zone in the base flow field, the length of the separation zone at the spanwise central cross section increased slightly with single-channel actuation, and the start of separation moved forward. However, no obvious changes were identified in the separation zone on both sides of the plate.

**Figure 18.** The form of separation zone at *y* = 0 cross section: (**a**) base flow field; (**b**) actuation flow field.

Based on the abovementioned results, when a 10 kHz single-channel actuation was applied, the intensity of the separation shock wave and the reattached shock wave in the ISWBLI flow field remained stable (barely weakened). Similarly, the size and shape of the separation zone did not change significantly either. Therefore, it was concluded that single-channel actuation can only just achieve the control effect over the ISWBLI flow field.

#### 4.2.3. Array Actuation Flow Field

The effect of array actuation on the ISWBLI was further simulated. As shown in Figure 19, array actuation was applied at the central cross section of the plate; the singlechannel actuation energy was set to 10 mJ; the actuation frequency was 10 kHz; the actuation point was 70 mm away from the plate's leading edge; and the actuation spacing was 15 mm. The array actuator was located upstream of the ISWBLI region, about 15~20 mm away from the interaction region.

**Figure 19.** The schematic diagram of array actuation model.

Figure 20 shows the numerical schlieren image of the time-averaged actuation flow field and the time-averaged base flow field at the *z* = 0 cross section. As shown in Figure 20a, at the *z* = 0 cross section, after the array actuation was applied, the separation shock wave basically disappeared, and the angle of the reattached shock wave was reduced from 34.6◦ in the base flow field to 30.5◦. These changes indicate that the array actuation can effectively weaken the intensity of the separation shock and reattached shock waves. Further, it was concluded that, after the application of actuation, the height of the separation zone at the *z* = 0 cross section was significantly increased. The height of the separation zone in the base flow field was about 5 mm, but, in the separation zone in the actuation flow field, it rose to 10 mm as the start of separation moved forward by about 19 mm. It can be concluded that the application of plasma actuation increases the size of the separation zone at the spanwise middle cross section and causes the start of separation to move forward. Upon further comparison of the time-averaged numerical schlieren, it was found that the intensity and angle of the incident shock wave remained unchanged, indicating that plasma actuation has no significant effect on the incident shock wave. This outcome is consistent with the experimental results presented in the previous section.

**Figure 20.** Numerical schlieren at *z* = 0 cross section: (**a**) 10 kHz actuation flow field; (**b**) base flow field.

The flow field at different cross sections was also analyzed. Figure 21a,b represents the numerical schlieren at the *z* = 0.025 m cross section and shows that the array actuation effectively weakened the intensity of the separation shock wave. Figure 21c,d show the numerical schlieren at the *z* = 0.054 m cross sections at the edge of the plate. The results show that the structure and intensity of the separation shock wave were not significantly affected after the actuation was applied. Hence, it can be concluded that actuation in the flow direction can hardly exert a full, spanwise control effect. In subsequent studies, the coupled layout of actuators in spanwise and flow directions should be considered to achieve a better control effect.

**Figure 21.** Numerical schlieren of two cross sections: (**a**) *z* = 0.025 m base flow field; (**b**) *z* = 0.025 m actuation flow field; (**c**) *z* = 0.054 m base flow field; (**d**) *z* = 0.054 m actuation field.

Complementarily, the interaction region of the ISWBLI at the *y* = 0 cross section was also examined. Figure 22a,b shows the velocity contours at the ISWBLI zone with and without actuation. The images show that, after the array actuation was applied, the start of separation at the spanwise middle position moved forward, and the length of the separation zone was increased. At the same time, the length of the separation zone near the edge of the plate was reduced, and the spanwise, affected area of the separation zone was also diminished.

**Figure 22.** Velocity contour of ISWBLI *y* = 0 cross section: (**a**) base flow field; (**b**) actuation flow field.

Figure 23 shows the curve of the skin friction coefficient at the *z* = 0 cross section in the *x*-direction. The red and blue lines represent the base flow field and the actuation flow field, respectively. After the array actuation was applied, the length of the separation zone at the *z* = 0 cross section increased from 20 mm to 50 mm, which was an increase of 125%. The start of separation in the base flow field was about 120 mm from the plate's leading edge, and the start of separation in the actuation flow field moved forward, about 100 mm from the plate's leading edge. These changes show that array actuation is likely to increase the length of the separation zone at the *z* = 0 cross section, making the start point of separation move forward.

**Figure 23.** Skin friction coefficient curve at *z* = 0 cross section.

The other spanwise cross sections of the flow field were the last to be examined. Figure 24a,b shows the curves of the skin friction coefficient of the flow field with and without actuation at the *z* = 0.025 m cross section and *z* = 0.054 m cross section, respectively. At the *z* = 0.025 m cross section, array actuation resulted in the decreased length of the separation zone, making the separation point move backward. At the *z* = 0.054 m cross section, the length of the separation zone was close to 0, and the Cf curves with and without actuation tended to be consistent. This result indicates that plasma actuation has a poor control effect over the separation zone at the edge of the plate.

**Figure 24.** Skin friction coefficient curve at different spanwise cross sections: (**a**) *z* = 0.025 m; (**b**) *z* = 0.054 m.

#### **5. Conclusions**

In this research, experiments and simulations were carried out to closely examine the control of the ISWBLI through pulsed arc discharge plasma actuation. The high-speed schlieren technology was used to study the control effect of actuation on the shock intensity. The spanwise, three-dimensional characteristics of the flow field with and without actuation were analyzed using the DES method. The major conclusions are presented below:


**Author Contributions:** C.Z. performed the experimental investigation. Y.L. performed the numerical simulations and analyzed the results; H.L. and S.G. supervised the work and reviewed and edited the manuscript; C.Z. and H.Y. wrote the main manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was founded by the National Science and Technology Major Project (Grant No. J2019-II-0014-0035).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data and materials are available upon request.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


### *Article* **Transonic Buffet Active Control with Local Smart Skin**

**Kai Ren 1, Chuanqiang Gao 1,\*, Fangqi Zhou <sup>2</sup> and Weiwei Zhang <sup>1</sup>**


**Abstract:** Transonic flight has high economic benefits, but the appearance of transonic buffet limits the flight envelope. The shock control bump currently used for transonic buffet suppression tends to degrade the aerodynamic performance of the non-buffeting state. In this study, a smart skin system is used to eliminate the fluctuating load of transonic buffet by measuring the airfoil lift coefficient as the feedback signal and adjusting the local skin height using data-driven, model-free adaptive control. Since the actuator height is dynamically adjusted only after the occurrence of transonic buffet, the smart skin can completely suppress the fluctuating load and does not affect the aerodynamic performance in the non-buffeting state. The suppression effect of the proposed smart skin on transonic buffet is verified by numerical simulation of the flow. The simulation results show that due to the introduction of closed-loop control, the fluctuating load of transonic buffet can be effectively suppressed for different positions and maximum heights of the actuator. Even when the flow state changes, the robust smart skin system can also achieve the control goal. Therefore, smart skins combining flexible materials and control technologies have the potential to effectively improve the aerodynamic performance of aircraft.

**Keywords:** smart skin; transonic buffet; closed-loop flow control; fluctuating load suppression

#### **1. Introduction**

Transonic flight has high economic benefits, but when the flow state enters above the transonic buffet boundary [1], which consists of Mach number (*M*) and angle of attack (*α*), transonic buffet flow [2] occurs on the wing surface. This is a kind of flow instability phenomenon caused by shock wave/boundary layer interaction, which is manifested as a large oscillation of shock wave on the wing surface, causing fluctuation of lift, affecting flight comfort and safety, easily causing fatigue damage to the wing structure, and limiting the flight envelope in the transonic section. In order to eliminate the fluctuating load of transonic buffet and improve the flow stability, various passive and active control schemes are used to modify the flow in the shock wave/boundary layer interaction region and trailing edge region [3]. The control devices mainly include vortex generator [4,5], shock control bump [6,7], streamwise slot [8], plasma actuator [9,10], and trailing edge flap [11,12]. Unfortunately, due to the complexity of shock buffet, a robust control scheme for eliminating shock buffet load has not been obtained in experiments. Researchers seek to analyze the advantages and disadvantages of different control schemes through numerical simulation, so as to provide some support for the realization of shock buffet control. The active control with the trailing edge flap as the actuator can effectively eliminate the fluctuating load, but it will occupy the control channel of the flap surface. Passive control represented by shock control bump (SCB) is widely studied because of its simple structure and effective improvement of transonic buffet boundary. However, the problem of SCB is that when the position and shape of the SCB are determined, it can only eliminate the buffeting load in the partial incoming flow state, and will reduce the aerodynamic performance in

**Citation:** Ren, K.; Gao, C.; Zhou, F.; Zhang, W. Transonic Buffet Active Control with Local Smart Skin. *Actuators* **2022**, *11*, 155. https:// doi.org/10.3390/act11060155

Academic Editors: Ioan Ursu, Hui Tang and Xin Wen

Received: 29 April 2022 Accepted: 8 June 2022 Published: 10 June 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the non-buffeting state. Therefore, the smart skin technology of a deformable wing [13] can combine the advantages of both control schemes mentioned above. The shape of the bump can be dynamically adjusted after transonic buffet occurs, which has the potential to effectively eliminate the transonic buffeting load without affecting the aerodynamic performance of the non-buffeting state.

With the development of new intelligent materials and distributed sensing and control technology, the deformable wing technology is evolving towards the direction of intelligence and flexibility. The use of intelligent materials such as piezoelectric materials [14] and shape memory alloys [15] to drive the deformation of aircraft and achieve the integration of flexible skin materials, structure, and actuator is a current research hotspot and a development direction for future applications. Table 1 compares the actuation performance of two intelligent materials: piezoelectric materials and shape memory alloys. Among the advantages of the piezoelectric materials is the fast response speed, while the shape memory alloy has the advantage of large output force and output displacement, compared to conventional wings that improve the aerodynamic performance of an aircraft under only a few specific flight conditions through the simple rigid deformation of flaps, leading edge slats, and other mechanisms. The intelligent deformable wings composed of intelligent deformable materials combined with distributed sensors and microcontrollers do not produce discontinuous changes and abrupt curvature changes, or reduce the separation of flow, and can optimize deformation parameters for each flight condition, which has better flight performance [16,17]. The deformable wing includes folding wing [18], flexible trailing edge [19], adaptive wingtip [20], active camber wing [21] and smart skin [22]. To enhance the performance of the aircraft, researchers use wings with variable camber [23,24], variable thickness [25], variable trailing edge deflection [19,26], and deformable winglets [27] to improve aerodynamic performance; use adaptive wings with adjustable maximum height and maximum height position to improve stall aerodynamic performance [28]; use deformable wings to adapt to different flight environments [29]; use flexible trailing edges to reduce structural weight of swept-wing aircraft [30]; use deformable trailing edges to control the boundary layer flow to reduce noise [31]; use adaptive wingtips to reduce fuel consumption and gust loads through cant, twist, and camber deformation [20]; use smart skin to adjust the inlet shape to make the engine have good performance over a wide range of rotation speed range [32]; and use smart skin to suppress panel flutter [33] and wing flutter [34].

**Table 1.** Comparison of intelligent material actuation performance.


In this study, smart skin technology is used to suppress the fluctuating load of transonic buffeting flow. Smart skin technology forms the neural network of aircraft by implanting intelligent structures within the aircraft skin and other components, including detection elements (sensors), microprocessor control systems (signal processors), drive elements (micro actuator), and interconnecting lines, etc. It can not only sense its own physical conditions (such as structural damage), but also be sensitive to the external environment (such as incoming flow condition). In this way, the components and even the whole aircraft body are endowed with self-monitoring, self-correction, self-adaptation, and memory, thinking, judgment, and response functions.

This paper is organized as follows: In Section 2, the composition of the smart skin system and the flow numerical simulation method to simulate the control performance are introduced. In Section 3, the adopted control method is introduced and simulations are performed to verify the effectiveness of the method in suppressing the transonic buffet fluctuating load. Section 4 analyzes the influence of smart skin system parameters on the control effect and tests the robustness of the system to different incoming flow states. Conclusion and discussion are presented in Section 5.

#### **2. Smart Skin System and Simulation Method**

#### *2.1. Smart Skin System*

The smart skin system includes sensors, microprocessor control systems, micro actuators, and interconnecting lines. The feedback signal used for the active control of transonic buffet is mainly the force coefficient and the wall pressure. Compared with the wall pressure, the lift coefficient (*Cl*) is obtained by the integration of the pressure distribution (*Cp*), which is smooth and beneficial to the stability of the system. Therefore, the lift coefficient is used as a feedback signal in the smart skin system. Since there is a direct proportional relationship between the lift force and the wing root bending moment, the lift coefficient is obtained by measuring the wing root bending moment using strain gauges. The design of the control law is described in Section 3.1.

The micro actuator uses the local deformation skin with adjustable height; the shape of deformation skin is a widely used Hicks–Henne shape function [7], as shown in Figure 1. The shape parameters include: (1) the length of deformation skin, *lb*, (2) the starting position of deformation skin, *x*0, (3) crest location, *xh*max*/lb*, and (4) the height of deformation skin, *hb* (defined as the ratio of the actual height to the chord length). The height of deformation skin is output by the control law and can be adjusted dynamically. The Hicks–Henne shape function is

$$\begin{array}{l} f\_b(\mathbf{x}\_b) = h\_b H(\mathbf{x}\_b) \\ H(\mathbf{x}\_b) = \sin^4(\pi \mathbf{x}\_b^m), \quad m = \ln(0.5) / \ln(\mathbf{x}\_{\text{hmax}} / l\_b) \end{array} \tag{1}$$

where, *xb* is the dimensionless chordwise position of deformation skin, which is defined as 0 ≤ *xb* = (*x* − *x0*)/*lb* ≤ 1. In this study, the length of deformation skin is taken to be 0.2 times the chord length, *lb* = 0.2 *c*, and the shape of deformation skin is symmetric, and the crest location is *xh*max*/lb* = 0.5.

**Figure 1.** Diagram of smart skin micro actuator.

#### *2.2. Flow Numerical Simulation Method*

In this study, the control performance of the smart skin system is analyzed by flow numerical simulation, and the numerical simulation uses the unsteady Reynolds-averaged Navier–Stokes (URANS) solver based on the finite volume method. The two-dimensional compressible URANS governing equations can be expressed as:

$$\frac{\partial}{\partial t} \iiint\limits\_{\Omega} \mathcal{W}dV + \iint\limits\_{\partial \Omega} F^{i}(\mathcal{W}, V\_{\text{grid}}) \cdot \mathfrak{n}dS = \iint\limits\_{\partial \Omega} F^{\sigma}(\mathcal{W}) \cdot \mathfrak{n}dS \tag{2}$$

where *W* is the flux variables, *F<sup>i</sup>* (*W*, *Vgrid*) is the inviscid flux, *Fv*(*W*) is the viscous flux, *Vgrid* is the grid velocity, Ω is the control volume, *∂*Ω is the control volume boundary, and *n* is the unit vector of normal direction outside the boundary.

The S-A turbulence model is introduced to close the Reynolds-averaged equation, and the spatial discretization and time integration of turbulent and mean flows are carried out in a loosely coupled way. The second-order advection upstream splitting method (AUSM) scheme is used to evaluate the inviscid flux with a reconstruction technique. The viscous flux term is discretized by the standard central scheme. In the turbulence model, the convective term is discretized by the second-order AUSM scheme and the destruction term by the second-order central scheme. For unsteady computations, the dual time stepping method is used to solve the governing equations. At the sub-iteration, the fourth-stage Runge–Kutta scheme is used with a local time stepping.

Due to the airfoil shape change during the numerical simulation, a grid deformation method is required to ensure the update and deformation of the airfoil surface and the solution domain grid. The grid deformation scheme uses the radial basis function (RBF) interpolation method [35].

The controlled object adopts the transonic buffet of RAE2822 airfoil, and the computational flow domain is discretized by a hybrid unstructured grid. The far field extends 20 chords away from the airfoil. There are 12,958 surface nodes and 40 layers of structured viscous grids around the airfoil. The distance between the first layer and the wall in the perpendicular direction is 5 × <sup>10</sup>−<sup>6</sup> chords (y<sup>+</sup> < 1). For a transonic flow: *<sup>M</sup>* = 0.73, *<sup>α</sup>* <sup>=</sup> 3.19◦, and *Re* = 6.5 × <sup>10</sup>6. Figure <sup>2</sup> shows the comparison between the calculated pressure distribution and the experimental result [36], which is basically consistent with each other. By fixing the Mach number of the incoming flow and simulating different angles of attack (the interval of angles of attack is 0.1◦) it can be determined whether the transonic buffet occurs. The buffet onset boundary obtained is shown in Figure 3, which is almost consistent with the onset boundary [7] calculated by Tian et al.

**Figure 2.** Comparison of pressure distribution on airfoil surface.

**Figure 3.** Comparison of transonic buffet onset boundary.

Choose the incoming flow state: *<sup>M</sup>* = 0.75, *<sup>α</sup>* = 4◦, and *Re* = 6.5 × 106 to simulate the time history of buffet. It can be seen from Figure 3 that this flow state is above the buffet onset boundary, and when the control is not turned on, the system will have a large lift fluctuation, as shown in Figure 4. Figure 5 shows the flow field pressure contours at different moments, corresponding to the four moments in Figure 4 (*a* and *c* are the middle positions when the shock wave moves upstream and downstream, respectively, and *b* and *d* are the most upstream and downstream positions of the shock wave, respectively). When the shock wave moves upstream, the length of the trailing edge separation bubble increases, and when the shock wave moves downstream, the separation bubble at the foot of the shock wave and the trailing edge separation bubble approach and merge gradually. From the pressure contours, it can be seen that the dimensionless chordwise range of the shock wave oscillation for this flow condition is 0.45~0.56, and the shock wave is at the most downstream position when the lift coefficient is maximum, and vice versa.

**Figure 4.** Lift coefficient time history.

**Figure 5.** Local pressure contours and streamlines of flow field at different moments: (**a**) moment *a*; (**b**) moment *b*; (**c**) moment *c*; (**d**) moment *d*.

#### **3. Transonic Buffet Control**

#### *3.1. Model-Free Adaptive Control*

Model-free adaptive control [37] is a data-driven control method. Its basic idea is to equivalently transform the discrete-time nonlinear system into a dynamic linearized data model based on the input and output data of the system at each moment using a pseudo partial derivative (PPD). Without establishing an accurate mathematical model of the flow system, the data-driven adaptive control of the flow system is achieved by estimating the system PPD online using the input and output data of the flow system, and then minimizing a given cost function.

For the transonic buffet flow system in this study, the feedback signal is the lift coefficient, *Cl*, the control signal is the actuator height, *hb*, and the controlled plant can be described by the following single-input, single-output, discrete-time nonlinear system:

$$\mathbb{C}\_{l}(k+1) = f(\mathbb{C}\_{l}(k), \cdot, \cdot, \mathbb{C}\_{l}(k - n\_{o}))h\_{b}(k), \cdot, \cdot, h\_{b}(k - n\_{i})) \tag{3}$$

where *ni* and *no* are unknown input and output delay orders, respectively, and *f*(·) is an unknown nonlinear function describing the relationship between lift coefficient and actuator height.

Based on PPD, the above flow plant is equivalently transformed into the following dynamic linearized data model:

$$\mathbf{C}\_{l}(k+1) = \mathbf{C}\_{l}(k) + \phi(k)\Delta h\_{b}(k)\tag{4}$$

where *φ*(*k*) is the PPD of system (3), it is estimated from the historical input and output data of the system by minimizing the following cost function:

$$J(\phi(k)) = \left| \mathbb{C}\_l(k) - \mathbb{C}\_l(k-1) - \phi(k)\Delta h\_b(k-1) \right|\_2 + \mu \left| \phi(k) - \hat{\phi}(k-1) \right|^2 \tag{5}$$

The estimation algorithm of the PPD (where *μ* > 0 is the weighting factor):

$$
\hat{\phi}(k) = \hat{\phi}(k-1) + \frac{\eta \Delta h\_b(k-1)}{\mu + \Delta h\_b(k-1)^2} (\Delta \mathcal{C}\_l(k) - \hat{\phi}(k-1) \Delta h\_b(k-1)) \tag{6}
$$

*η* ∈ (0, 1] is the update step of PPD, which can make the algorithm more flexible, and *η* = 1 is taken in all of this paper.

In order to eliminate fluctuating load of the transonic buffet flow and improve the flow stability, the following cost function are considered:

$$f(h\_b(k)) = \left| \mathbb{C}\_l^\*(k+1) - \mathbb{C}\_l(k+1) \right| \left| 2 + \lambda\_1 \right| h\_b(k) - h\_b(k-1) \left| 2 + \lambda\_2 \right| h\_b(k) - h\_{b0} \big|^2 \tag{7}$$

where *C*∗ *<sup>l</sup>* (*k +* 1) is the desired lift coefficient, *hb*<sup>0</sup> is the initial height of the actuator, and *λ*<sup>1</sup> and *λ*<sup>2</sup> are positive weighting factors. The first term of the cost function is used to limit the deviation between the lift coefficient and the desired lift coefficient, and the second term is used to limit the change of the actuator height. When the unstable steady flow of transonic buffet cannot be obtained in advance, the corresponding desired lift coefficient cannot be obtained, so the desired lift coefficient is taken as the time-averaged lift coefficient, and the deviation between the actuator height and the initial height of the actuator is limited by introducing the third term of the cost function.

Combined with the dynamic linearized data model (4) to minimize the cost function (7), the control law is:

$$h\_{\mathbb{B}}(k) = \frac{1}{\lambda\_1 + \lambda\_2 + \phi(k)^2} \left( \phi(k) \left( \mathbb{C}\_l^\*(k+1) - \mathbb{C}\_l(k) \right) + \left( \lambda\_1 + \phi(k)^2 \right) h\_{\mathbb{B}}(k-1) + \lambda\_2 h\_{\mathbb{B}0} \right) \tag{8}$$

The control system block diagram is shown in Figure 6. The loop in the upper part of the block diagram is an adaptive loop, indicating that the adaptive law uses the input and output data of the system to estimate the PPD online and adjust the controller parameters when the plant parameters change, so as to make the system work in an optimal or nearoptimal state. The loop of the lower part is a general feedback loop, which is obtained by minimizing the cost function.

**Figure 6.** Block diagram of model-free adaptive control system.

#### *3.2. Numerical Simulation and Mechanism Analysis of Control Law*

For the incoming flow state of *<sup>M</sup>* = 0.75, *<sup>α</sup>* = 4◦, and *Re* = 6.5 × <sup>10</sup>6, it can be seen from Section 2.2 that this state is above the buffet onset boundary. When the control is not turned on, the dimensionless chordwise range of the shock wave oscillation is 0.45~0.56. Smart skin system is used to eliminate the fluctuating load in this state. According to the analysis results by Tian et al., the fluctuating load can be effectively suppressed when the bump is located after the shock wave [7], so the range of actuator is chosen to be 0.55~0.75, that is, the starting position of deformation skin *x*<sup>0</sup> = 0.55, and the length of deformation skin *lb* is 0.2 *c*.

The transonic buffet fluctuating load in this state is suppressed by the above control system configuration. The time history of control signal (actuator height, *hb*), feedback signal (lift coefficient, *Cl*), and drag coefficient (*Cd*) are shown in Figure 7. It can be seen that when the control is turned on, the fluctuating load (*Cl* and *Cd*) has been completely eliminated after a few buffeting cycles, and the actuator height has returned to the initial 0, that is, the airfoil returns to the shape of RAE2822 airfoil. The maximum height of actuator during the control process is 0.00299; there is no additional increase in airfoil drag, and the final steady flow field pressure contour is shown in Figure 8.

**Figure 7.** Time history of the system: (**a**) lift coefficient; (**b**) drag coefficient; (**c**) actuator height.

**Figure 8.** Flow field pressure contour and streamlines after control convergence.

Figure 9 shows the time history of the shock wave dimensionless position *xsw/c* and the height of the actuator. It can be seen that the height of the actuator located at the downstream of the shock wave is close to the lowest when the shock wave is at the most downstream, and the height of the actuator is close to the highest when the shock wave is at the most upstream. By changing the height of the actuator, the oscillation of the shock wave is limited, thus eliminating the fluctuating load of the buffet. Figure 10 shows the pressure contour and streamlines for the four moments in Figure 9; *b'* and *d'* are the middle positions where the shock wave moves downstream and upstream, respectively, and *a'* and *c'* are the most upstream and downstream positions of the shock wave, respectively. When the shock wave moves from upstream to downstream, the height of the actuator decreases, which limits the movement of the separation bubble at the shock wave foot to downstream and prevents the merging of the separation bubble at the shock wave foot with the trailing edge separation bubble. When the shock wave moves upstream, the actuator height increases, which limits the growth of the trailing edge separation bubble upstream.

**Figure 9.** Time history of the shock wave position and the height of the actuator.

**Figure 10.** Local pressure contours and streamlines of flow field at different moments: (**a**) moment *a'*; (**b**) moment *b'*; (**c**) moment *c'*; (**d**) moment *d'*.

#### **4. Smart Skin Performance Analysis**

#### *4.1. The Influence of Actuator Position*

In the previous section, for the incoming flow state of *<sup>M</sup>* = 0.75, *<sup>α</sup>* = 4◦, and *Re* = 6.5 × 106, the dimensionless chordwise range of the shock wave oscillation is 0.45~0.56, and the actuator range is chosen to be 0.55~0.75. In this section, the influence of actuator position on the buffeting load suppression is analyzed. Two actuator positions are selected: (1) the actuator position is upstream of the shock wave oscillation range (actuator range: 0.25~0.45); (2) the actuator position coincides with the shock wave oscillation range (actuator range: 0.4~0.6). The input and output time history of control system are shown in Figures 11 and 12, respectively. It should be noted that in Section 3.2, we analyzed that for the actuator located downstream of the shock wave, when the shock wave was at the most downstream, the height of the actuator was close to the highest. For the actuator located upstream of the shock wave (actuator range: 0.25~0.45), the height of the actuator should be close to the lowest when the shock wave is at the most downstream. That is, the initial value of the PPD in the control law should be opposite to the sign in Section 3.2.

**Figure 11.** Time history of the system (actuator range: 0.25~0.45): (**a**) lift coefficient; (**b**) actuator height.

**Figure 12.** Time history of the system (actuator range: 0.4~0.6): (**a**) lift coefficient; (**b**) actuator height.

From Figures 11 and 12, we find that the required actuator height is smaller when the actuator position is close to the oscillation range of shock wave. As shown in Figure 13, we compare the maximum height of the actuator for different crest chordwise positions under the premise of basically the same settling time. Because the function of the actuator is to suppress the oscillation of the shock wave and to prevent the merging of the separation bubble at the shock wave foot and the trailing edge separation bubble, the position of the actuator far away from the shock wave needs a larger height to prevent the merging of the two separation bubbles, and to generate enough pressure disturbance to suppress the oscillation of the shock wave.

**Figure 13.** Relationship between maximum height of actuator and crest chordwise position.

#### *4.2. The Influence of Actuator Height*

Real wings design often has multiple constraints, such as the wing volume will limit the amplitude of the actuator subsidence, and actuator structure and drag constraints will limit the maximum height of the actuator. Therefore, this subsection will analyze the performance change of the smart skin when the actuator height has the maximum amplitude limit *hs* for the incoming flow state of *<sup>M</sup>* = 0.75, *<sup>α</sup>* = 4◦, and *Re* = 6.5 × 106, and the actuator range of 0.55~0.75. As can be seen in Figure 14, the smart skin system can eliminate the transonic buffet fluctuating load under different restrictions of maximum height of the actuator, but with the decrease of the maximum actuator height, the system settling time gradually increases. According to the previous smart skin control mechanism, the phase relationship between the actuator height and the shock wave position (or lift coefficient) plays a decisive role, while the height only affects the dynamic characteristics such as the settling time.

**Figure 14.** Influence of actuator height constraint *hs* on dynamic characteristics of the control system: (**a**) time history of lift coefficient; (**b**) time history of actuator height.

As shown in Table 2, the required actuator height for active control can be smaller than that for a passive bump device located immediately downstream of the shock wave, which reduces the requirement for the actuator and is conducive to engineering implementation. Although compared with passive control, active control requires energy input from the time history of input and output of the control system; the output displacement of actuator and settling time are small, so the energy input required is acceptable.


**Table 2.** Maximum height of actuator required for different control strategy.

#### *4.3. Robustness of The Control Law*

Due to the existence of atmospheric turbulence, the incoming flow state may change, and the dynamics characteristics of the transonic buffet flow system will change, which is manifested in the change of the oscillation range of shock wave, the mean value, and fluctuation amplitude of lift, such as for the other two states above the buffet onset boundary: (1) *<sup>M</sup>* = 0.73, *<sup>α</sup>* = 4.8◦, and *Re* = 6.5 × <sup>10</sup>6, and (2) *<sup>M</sup>* = 0.71, *<sup>α</sup>* = 5.5◦, and *Re* = 6.5 × <sup>10</sup>6. In state (1), oscillation range of shock wave is from 0.35~0.55, and fluctuation range of lift coefficient is 0.723~0.994. In state (2), oscillation range of shock wave is from 0.24~0.48, and fluctuation range of lift coefficient is 0.673~1.042. Figures 15 and 16 show the lift coefficient, drag coefficient, and actuator height time history of the system using smart skin to suppress the fluctuating load in the above two incoming flow states, respectively, and the actuator range from 0.55 to 0.75. It can be seen from the Figures 15 and 16 that the control strategy in this study has strong robustness, and after turning on the control, by adjusting the height

of smart skin, the large fluctuating load of transonic buffet can be completely eliminated, and there is no additional increase in airfoil drag. Compared with state (1), the angle of attack of state (2) increases, the shock wave position moves forward, and the fluctuating amplitude of lift coefficient increases, so the maximum height required to suppress the fluctuating load in state (2) increases.

**Figure 15.** Time history of the system (*M* = 0.73, *α* = 4.8◦): (**a**) lift coefficient; (**b**) drag coefficient; (**c**) actuator height.

**Figure 16.** Time history of the system (*M* = 0.71, *α* = 5.5◦): (**a**) lift coefficient; (**b**) drag coefficient; (**c**) actuator height.

#### **5. Conclusions**

In this study, the smart skin system is used to suppress the fluctuating load of transonic buffet flow for RAE2822 airfoil. The smart skin system uses the lift coefficient as feedback signal and adopts a model-free adaptive control to dynamically adjust the local skin height of airfoil. When the incoming flow state enters above the transonic buffet onset boundary, by preventing the merging of the separation bubble at the shock wave foot with the trailing edge separation bubble, the smart skin system can completely suppress the fluctuating load and does not affect the aerodynamic performance in the non-buffeting state. From the time history of the control system input and output, the system has a short settling time and does not need too large an output displacement, so the piezoelectric material with high precision of displacement control, fast response speed, low power consumption, and small size can be used as the actuator.

By analyzing the control performance of the smart skin system, the following conclusions were obtained:


(3) When the incoming flow state changes, although the oscillation range of shock wave and the fluctuation range of lift coefficient change greatly, the buffeting load in different incoming flow states can be suppressed by the model-free adaptive control with the mean value of the lift coefficient as the desired output of the control system.

The above research is carried out in a simulation environment in order to provide some theoretical support for the realization of shock buffeting control. To realize the control of shock buffet, more detailed parameters, such as the specification and distribution of sensors and the bandwidth and hysteresis of actuators, need to be further studied.

With the progress of material and control technology, new flexible controllers such as smart skin can have more flexible control methods without abrupt curvature changes on the surface of the aircraft, and have the ability to obtain the best control effect through adaptive and intelligent control for different flight conditions, which can greatly improve the aerodynamic performance of the aircraft.

**Author Contributions:** K.R. designed the controller and writing-original draft preparation.; C.G. conceived the original ideals, review; F.Z. proposed the comment for the paper; W.Z. provided research methods, supervision. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China (No. 11902269), the Aviation Science Foundation of China (No. 2019ZH053003), the Fundamental Research Funds for the Central Universities, and State Key Laboratory of Aerodynamics innovation Fund Project (JBKYC190102).

**Data Availability Statement:** Not applicable.

**Acknowledgments:** Fangqi Zhou is particularly grateful for the support of the State Key Laboratory of Aerodynamics innovation Fund Project.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Influence of Synthetic Jets on Multiscale Features in Wall-Bounded Turbulence**

**Biaohui Li 1, Jinhao Zhang <sup>1</sup> and Nan Jiang 1,2,\***


**Abstract:** This experimental research focuses on the impacts of submerged synthetic jets on a fullydeveloped turbulent boundary layer (TBL) under a drag reduction working case. Two-dimensional velocity vectors in the flow field are captured with the aid of a particle image velocimetry (PIV) system. Proper orthogonal decomposition (POD) analyses provide evidence that synthetic jets notably attenuate the induction effect of prograde vortex on the low-speed fluid in large-scale fluctuation velocity field, thereby weakening the bursting process of near-wall turbulent events. Furthermore, the introduced perturbance redistributes the turbulent kinetic energy (TKE) and concentrates the TKE onto small-scale coherent structures. Modal time coefficients in various orders of POD are divided into components of multiple frequency bands by virtue of complementary ensemble empirical mode decomposition (CEEMD). It is found that the turbulence signals are shifted from low-frequency to high-frequency bands thanks to synthetic jets, thus revealing the relationship between scales and frequency bands. One further method of scale decomposition is proposed, that is, the large-scale fluctuating flow field will be obtained after removing the high-frequency noise data with the help of continuous mean square error (CMSE) criterion.

**Keywords:** active control; multiscale; synthetic jet; turbulent boundary layer

#### **1. Introduction**

Wall turbulence drag reduction technologies have represented a dramatic topic in aerospace and other engineering applications over the past few decades as a result of the oil crisis in the early 1970s [1]. So far, various control schemes [2] have been put forward with remarkable success and are of great significance in improving the performance of vehicles and saving energy [3]. The consensus is that effective drag reduction is achieved by impairing coherent structures [4] in the turbulent boundary layer (TBL). Generally, passive approaches that only apply to a single design condition are reliable and well-studied. In contrast, active control techniques are too complicated to be attained, which has become a vital research direction to be conquered at present [5]. Among them, the accomplishment of skin friction drag reduction by employing blowing and suction at the wall is a widely studied strategy [6]. Fortunately, the synthetic jet actuator (SJA) has the ability to generate local blowing and suction without additional air supply. This facilitates interaction with vorticity structures in the TBL [7], which has tremendous potential.

Ricco et al. [8] reviewed sundry state-of-the-art techniques oriented toward turbulent drag reduction via spanwise actuation, including wall motion, plasma body forces and synthetic jets, pointing out the preponderance of synthetic jets for resistance reduction control. The interaction experiment of the inclined circular synthetic jets with turbulent channel flow was carried out by Iuso et al. [9]. Reductions in the average surface friction of up to 15% as well as in the turbulence pulsation of up to 12% were realized. Yao et al. [10] imposed a body force in the motion equation to simulate the spanwise motion induced

**Citation:** Li, B.; Zhang, J.; Jiang, N. Influence of Synthetic Jets on Multiscale Features in Wall-Bounded Turbulence. *Actuators* **2022**, *11*, 199. https://doi.org/10.3390/act11070199

Academic Editor: Luigi de Luca

Received: 4 July 2022 Accepted: 17 July 2022 Published: 18 July 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

by synthetic jets, successfully achieving an optimal net energy saving rate of 17%. In light of the research of Thomas et al. [11], spanwise synthetic jets generated by a pulsed DC plasma exciter enabled a resistance reduction margin of over 70%, with the drag reduction ratio increasing logarithmically as the number of streaks decrease. Cannata et al. [12] perturbed the TBL through a linear array of synthetic jets orthogonal to the streamwise direction and the coalescence of near-wall structures was observed. This meant a reduction in the number of coherent structures and a suppression of turbulent activity, enabling wall shear stress decreases. By means of direct numerical simulation (DNS), Xie et al. [13] verified that the wall shear stress aggrandized in the suction phase of synthetic jets and declined in the blowing phase, the overall drag reduction phenomenon was attributed to a finite counterflow resulting from the nonlinear interaction between the synthetic jets and mainstream.

Heretofore, the research on turbulent drag reduction control utilizing synthetic jets is primarily concentrated in the aviation industry. In the field of underwater drag reduction, traditional bionic technology is still adopted [14], which poses a huge challenge to researchers. In order to make up for the deficiency in this aspect, some scholars began to conduct studies about submerged synthetic jet drag reduction. The investigation of underwater synthetic jets distributed along the span and oriented in the wall-normal direction to simulate the spanwise force of an oscillating wall indirectly was completed by Segawa et al. [15]. They detected an appreciable drag reduction rate of up to 30% using the fiber Bragg grating (FBG) system. Spinosa et al. [16] adopted an array of round synthetic jets to modify the boundary layer flow field in a water channel. The artificial hairpin vortices of the same scale as that of the coherent structures in the near-wall region were generated by synthetic jets, accordingly, forming a drag reduction region with streamwise velocity deficit.

It is worth noting that a basic property of turbulence is its multiscale features [17], whereas current academia still lacks in-depth analyses on the regulation mechanism of synthetic jets on the scale characteristics of TBL. The present work revolves around how synthetic jets alter the multiscale properties of turbulence under drag reduction state. The snapshot proper orthogonal decomposition (POD) is employed to divide the scale of the flow field before and after the control is handled, then on this basis, a new approach of the scale definition is proposed on the grounds of complementary ensemble empirical mode decomposition (CEEMD). Eventually, the control effect of synthetic jets on turbulent burst events is further explored.

#### **2. Experimental Facilities and Procedure**

#### *2.1. Laboratory Equipment*

The active control investigations were accomplished at Tianjin University in a lowspeed recirculating water flume with a cross section of 6.2 m long, 0.3 m wide and 0.38 m high. The detailed schematic is exhibited in Figure 1. The dimensions of Plexiglas test plate placed horizontally was 2.9 m × 0.294 m × 0.02 m (length × width × thickness), and the front edge was trimmed with an ellipse of 4:1 along the long and short axes to prevent the flow separation. The tripwire, 2 mm in diameter, was mounted 100 mm downstream of the leading edge for accelerating the transition, and a trailing-edge flap at the tail of the platform was used to realize zero-pressure gradient flow. One 1 mm long, 10 mm wide spanwise slit was set at 2450 mm downstream of the central line of the front edge, and the SJA was located in the cavity below it. In the current work, the TBL flow developed at a mainstream velocity *u*∞ of 0.23 m/s with around 1.5% turbulence intensity, and the friction Reynolds number under the uncontrolled case was *Re<sup>τ</sup>* = 532. The coordinate axes *x*, *y* and *z* refer to the streamwise, wall-normal and spanwise directions, respectively, and *u*, *v* and *w* are the corresponding velocities. For convenience of contrast, the center of the synthetic jet outlet is designated as the origin of streamwise direction.

**Figure 1.** Experimental setup: (**a**) Front view; (**b**) Top view.

Figure 2 illustrates an enlarged view of the SJA. The rectangular slit is orientated vertical to the plate. A waterproof loudspeaker was selected as the exciter. It was attached to a cavity whose bottom was 38 mm in diameter and 18 mm deep. In a former study, we found that under the operating parameters of voltage 09 V, frequency 26 Hz and sinusoidal vibration waveform, SJA can produce a maximum drag reduction rate of about 17.27% [18]. Therefore, this working condition is elaborated in-depth below.

**Figure 2.** The synthetic jet actuator (SJA): (**a**) Front view; (**b**) Top view.

#### *2.2. Measurement Process*

The velocity vectors in the streamwise-wall-normal plane were acquired via the 2D-TR-PIV system. Note that the tracking of the fluid was realized by polyamide with a diameter of 20 μm and a density of 1.03 g/cm3. The *x*-*y* plane was illuminated by a laser with a maximum output power of 10 W and the light thickness of about 1 mm. The Speed Sense 9072 12-bit charge-coupled device (CCD) camera, faced the position of 50 mm downstream of the jet hole, was equipped with a 60 mm Nikon lens for recording, storing and processing particle images. Table 1 lists the parameters for measurement. The total sampling time was more than 68 s, ensuring statistical convergence. The transformation from particle images to velocity vectors was realized by dint of the Dynamic Studio software, and the pseudovector was expected to be less than 1%.


**Table 1.** Experimental parameters.

#### **3. Results and Analyses**

#### *3.1. Proper Orthogonal Decomposition*

For the two-dimensional two-component (2D-2C) velocity field, the snapshot POD [19] could be executed to availably separate the scale feature of flow. The specific idea is that the instantaneous field can be regarded as the superposition of the mean field and the fluctuation field:

$$\mathbf{x}'(t\_i) = \mathbf{x}(\boldsymbol{\xi}, t\_i) - \mathbf{\tilde{x}}(\boldsymbol{\xi}) \in \mathbb{R}^m, i = 1, 2, \dots, n, n \ll m \tag{1}$$

where *x* = (*u*, *v*), *m* and *n* are space and time series, respectively, and *m* = 92 × 156 × 2. Since POD is not sensitive to time resolution [20], and in order to reduce the computational burden, one snapshot is extracted for every five snapshots, so the time series is *n* = 8215. In view of this, the time scale equivalent is to be amplified by a factor of 5, and as a result, the value of frequency needs to be reduced by 5 times, accordingly, in the power spectral density (PSD) analysis of the time coefficients in the following. The specific process is as follows:

(1) The time series (fluctuation velocity) matrix is expressed as:

$$X = \left[ \mathbf{x}'(t\_1), \mathbf{x}'(t\_2), \mathbf{x}'(t\_3), \dots, \mathbf{x}'(t\_{n-1}), \mathbf{x}'(t\_n) \right] \in \mathbb{R}^{m \times n} \tag{2}$$

(2) Solving the eigenvectors *Ψ<sup>j</sup>* and eigenvalues *λ<sup>j</sup>* of the correlation matrix:

$$X^T X \Psi\_{\dot{\jmath}} = \lambda\_{\dot{\jmath}} \Psi\_{\dot{\jmath}\prime} \Psi\_{\dot{\jmath}} \in \mathbb{R}^n \tag{3}$$

(3) The basis function for POD mode is written as:

$$\Phi\_{\dot{j}} = \frac{X\Psi\_{\dot{j}}}{||X\Psi\_{\dot{j}}||} \in \mathbb{R}^{m}, j = 1, 2, \dots, n \tag{4}$$

(4) The time coefficient (expansion coefficient) of each mode will be realized as:

$$a\_{\dot{j}} = \Phi\_{\dot{j}}^T X \tag{5}$$

(5) Reconstructing the fluctuating velocity field:

$$\mathbf{x}'(\xi, t\_i) = \sum\_{j=1}^{n} a\_j(t\_i) \Phi\_j(\xi) \tag{6}$$

(6) Single-order mode turbulent kinetic energy ratio:

$$P(\lambda) = \frac{\lambda\_j}{\Sigma \lambda}, j = 1, 2, \dots, n \tag{7}$$

(7) The proportion of cumulative energy is defined as:

$$\mathcal{C}(\lambda) = \frac{\Sigma\_{j=1}^k \lambda\_j}{\Sigma \lambda}, 1 \le k \le n \tag{8}$$

The contribution of each POD mode to total turbulent kinetic energy (TKE) is given in Figure 3a. In the control state, the energy proportion is lower than that without control in the first four modes, and starting with the fifth mode, energy distribution of the two conditions is opposite. Thus, it can be inferred that the synthetic jets restrain the large-scale (low-order) TKE and strengthen the small-scale (high-order) turbulent events. Figure 3b plots the cumulative energy ratio. For the perturbed condition, more modes are required to achieve the same levels of TKE as the uncontrolled flow field, which is because the periodic disturbance generated by synthetic jets accelerates the migration of turbulent energy from large-scale to small-scale.

**Figure 3.** Flow field energy distributions: (**a**) Single energy; (**b**) Cumulative energy.

The levels of TKE of the aforementioned various modes are determined by normalized eigenvalues, so the amount of information of turbulent events contained will be attained by introducing information entropy [21]. The formula is as follows:

$$H(\lambda) = -\sum\_{j=1}^{k} P(\lambda\_j) \; \log P(\lambda\_j), 1 \le k \le n \tag{9}$$

where the base of the logarithm is two and the unit is bit. Information entropy is the expectation of information content, as a result, in turbulent events, low-order modes have higher probability, smaller information entropy, less uncertainty and events are ordered; high-order modes have lower probability, larger information entropy, more uncertainty and events are disordered.

The accumulated curve of information entropy along with modes is reported in Figure 4a. The cumulative *H*(*λ*) of the first seventeen modes in the forced case is smaller than that of the undisturbed state. Since the eighteenth mode, the *H*(*λ*) of the former exceeds the latter, and the difference between the two is further expanded as the number of modes increases. On the whole, the *H*(*λ*) of the flow field before and after the application of control is 6.13 and 6.53, respectively. In other words, the existence of the synthetic jet leads to disordered flow field. However, it should be pointed out that in the low-order modes range (1~17), the uncertainty of the turbulent events contained in the controlled condition is relatively low, and the flow field is more orderly, which is attributed to the perturbation exerting a positive modulation influence on the large-scale coherent structures.

**Figure 4.** The profiles of (**a**) information entropy and (**b**) relative entropy.

Based on the information entropy, Kullback et al. [22] proposed the concept of relative entropy, which can measure the disparity between two random distributions. The specific expression can be written as:

$$D\_{KL}(P\_0 \parallel P\_1) = \sum\_{j=1}^{k} P\_0(\lambda\_j) \log \frac{P\_0(\lambda\_j)}{P\_1(\lambda\_j)}, 1 \le k \le n \tag{10}$$

where *P*<sup>0</sup> and *P*<sup>1</sup> represent the energy content of single mode before and after exerting the synthetic jets, singly. Figure 4b displays the curve of relative entropy. The variation trend of the curve increases first and then decreases and reaches the peak at the fourth mode. This means that the synthetic jets cause the difference between the forced flow field and the initial case to aggrandize, and the increase is chiefly concentrated in the lower modes. Although the synthetic jets heighten the small-scale TKE, the *DKL* gradually diminishes in higher-order modes and infinitely approaches a constant value, hence certifying the dominant status of the large-scale coherent structures in the boundary layer flow field.

The scale features of each mode could be quantified by the autocorrelation analysis of the streamwise pulsation velocity along the flow direction, the spacing corresponding to the second peak of the correlation function is the streamwise scale. Figure 5 illustrates the spatial scale of the first 1000 modes in diverse wall-normal positions, and *δ* is the nominal thickness of the boundary layer in the uncontrolled case. In natural flow, the large-scale structures are mainly distributed in the first ten order modes, and some of the scale peaks appear near the wall. Once the synthetic jets were injected into the TBL, the spatial scale corresponding to the low-order modes is deadened to a certain extent at the same height while the status of small-scale structures is improved.

The PSD curves of the time coefficient are plotted in Figure 6. In the first mode, the energy is focused near the low-frequency region, and in the high-frequency region the energy accounts for a small proportion. In the low-frequency range, the PSD of unforced flow field is larger than that of the controlled state, reflecting the inhibition impact of synthetic jets on low-frequency signal. However, the consequence of Figure 6b is just inverse to that of Figure 6a, which is consistent with the conclusion of the energy contribution ratio for the aforementioned single-order mode. Starting from the thirteenth mode, the black curve has a peak at 26 Hz, the effect of external perturbance begins to come out. As for the eighteenth mode, the peak at the excitation frequency becomes larger, verifying that the perturbation introduced by synthetic jets belongs to high-frequency and small-scale.

**Figure 5.** The spatial scale in flow direction: (**a**) None; (**b**) Control.

**Figure 6.** Power spectral density: (**a**) Mode 1; (**b**) Mode 5; (**c**) Mode 13; (**d**) Mode 18.

In terms of the analysis process of Wu et al. [23], taking 50% of the TKE level as the threshold, the large- and small-scale fluctuation fields are reconstructed by Equation (5). For the variation of turbulent pulsation and spanwise vortex structure at different scales, the turbulent fluctuation velocity field can be studied by the conditional phase averaging. In this paper, the prograde vortex (*Λci* < 0) was adopted as the detection condition, the detection function expression is defined as:

$$
\begin{array}{l}
\langle u\_l'(\Delta \mathbf{x}, y) \rangle = \langle u\_l'(\mathbf{x} + \Delta \mathbf{x}, y) | \Lambda\_{\text{ci}}(\mathbf{x}, y\_r) < 0 \rangle \\
\langle u\_s'(\Delta \mathbf{x}, y) \rangle = \langle u\_s'(\mathbf{x} + \Delta \mathbf{x}, y) | \Lambda\_{\text{ci}}(\mathbf{x}, y\_r) < 0 \rangle
\end{array} \tag{11}
$$

where *u <sup>l</sup>* and *u <sup>s</sup>* stand for the large- and small-scale *u* , respectively, the wall-normal reference position *yr* = 0.15*δ*, where a complete spanwise vortex structure can be detected. The local swirling strength of the vorticity was identified by *λci* criterion [24], and the sign of local vorticity *ω* was assigned to *λci*, consequently, *Λci* = *λciω*/|*ω*|.

Figure 7 presents the condition average contours of large-scale streamwise fluctuating velocity field. The prograde vortex, whose vortex center is signified by the black dot, contributes great skin friction resistance due to the strong shear with the wall [25]. Compared with the undisturbed case, the large-scale positive fluctuation of the forced state is strengthened while the negative fluctuation is weakened. Moreover, the narrowing of the streamwise direction range in the low-speed region is conducive to cutting down the burst

intensity of turbulent events. The red point is the saddle point, where the large-scale highand low-speed fluids collide with each other, and the motion of the fluid is almost stagnant, contributing to the huge Reynolds shear stress [26]. In the disturbed case, the distance between left saddle point and vortex center becomes smaller and moves towards the far wall-normal direction, extremely lowering the mass, energy and momentum transport caused by the large-scale Reynolds shear stress near the wall.

**Figure 7.** Conditional averaging of large-scale flow field: (**a**) None; (**b**) Control.

According to the conditional mean contours of the small-scale flow fluctuation velocity field in Figure 8, both positive and negative pulsations in the controlled flow field are enhanced, and the small-scale vortex structure becomes more active. These results can be attributed to two aspects: (1) small-scale perturbance continuously arises from the synthetic jets; (2) the disturbance disrupts the large-scale coherent structures.

**Figure 8.** Conditional averaging of small-scale flow field: (**a**) None; (**b**) Control.

#### *3.2. Complementary Ensemble Empirical Mode Decomposition*

The time coefficients calculated by POD belong to one-dimensional time series, which can be analyzed and processed by referring to the empirical mode decomposition (EMD) brought forward by Huang et al. [27]. The core idea of EMD is to decompose one complex signal into several intrinsic mode functions (IMFs) and a residual. Each IMF contains the information of different time scales, and the residual denotes the slow change of the signal [28]. Considering the modal aliasing phenomenon in EMD, we adopt an improved approach, namely complementary ensemble empirical mode decomposition (CEEMD) [29]. The procedure is as hereunder mentioned:


(4) The mean value of the *h* group IMFs of the corresponding mode can be calculated to complete the decomposition.

The first mode time coefficient of POD is decomposed by CEEMD, and 12 IMFs are obtained, as manifested in Figure 9. These IMFs stand for pulsation information of different frequencies in the flow field and are arranged in order from high-frequency to low-frequency. Residual is not plotted in the graph. Regardless of the high-frequency or low-frequency IMFs, the fluctuation values are commonly concentrated near zero, making a small contribution to the original signal. The IMF amplitudes of the middle-frequency bands are large and show distinct periodicity. In the condition of synthetic jet flow, the amplitudes of high-frequency curves are markedly greater than that of the non-controlled case, and the amplitudes of the middle-frequency bands are attenuated, which implies that more high-frequency constituents are distributed in the perturbed flow field.

**Figure 9.** Intrinsic mode function: (**a**) None; (**b**) Control.

Figure 10 is a bar graph of the Pearson correlation coefficient (PCC) between IMFs and the residual with the original signal. In the natural flow, the contribution of IMFs to the time coefficient is relatively simple and focused in the low-frequency bands. By contrast, in the perturbed condition, the proportion of high-frequency IMFs cannot be ignored. This can be seen from the correlation between the first four IMFs with the initial signal is greater than that before applying control.

Regarding the reconstruction of the signal, the standard deviation of the correlation coefficient is usually taken as the threshold, and the part of the correlation coefficient greater than the standard deviation is retained and superimposed to get a denoised signal.

As can be seen from Figure 11, in the case of no control, the time series before and after reconstruction are almost the same, and the error between the two is small. Nevertheless, for the control condition, there is a certain deviation between the original and reconstructed time coefficient, the relative error is much larger than the former, meaning this reconstruction process is not suitable for a situation where high-frequency signals dominate.

**Figure 10.** Correlation coefficient of each IMF and original signal.

**Figure 11.** Correlation coefficient of each IMF and original signal: (**a**) None; (**b**) Control.

To prevent the useful information in the original sequence from being discarded, the noise data is discriminated by calculating the mean square error of two adjacent IMFs, that is, the continuous mean square error criterion (CMSE) [30]:

$$\text{CMSE}(c\_i, c\_{i+1}) = \frac{1}{n} \sum\_{t=1}^{n} [c\_i(t) - c\_{i+1}(t)]^2 = \frac{1}{n} \sum\_{t=1}^{n} [\text{AIMF}\_i(t)]^2, i = 1, 2, \dots, 11\tag{12}$$

where *c*(*t*) denotes the IMF. According to formula 12, a certain IMF with abrupt noise energy distribution can be calculated, the reconstructed signal will be acquired by summing up the remaining IMFs and residual quantity after abandoning the previous high-frequency components.

In Figure 12, for the first-order mode, before and after control, the minimum mean square error occurs at the 4th and 3rd IMFs, respectively, therefore, the signal will be reconstructed from the 5th and 4th IMFs. By the way, if the IMF of a certain order modal time coefficient does not have a sudden change in energy density, but shows a stepby-step decreasing trend, such as the 500th order mode shown in Figure 12b, which

means that turbulence signal in this mode is dominated by high-frequency components, in consequence, the time series is assigned a value of zero as the reconstructed signal.

**Figure 12.** The difference between adjacent intrinsic mode function: (**a**) Mode 1; (**b**) Mode 500.

Figure 13 indicates the comparison results of the reconstructed time sequence and the original time coefficient in the first-order mode. No matter what case, the reconstructed series shows a very high degree of coincidence with the original time coefficients. The correlation coefficient is close to one, and the relative error does not exceed 1%, which meets the expected requirements.

**Figure 13.** Contrast of reconstructed signal and original signal: (**a**) None; (**b**) Control.

Table 2 sums up the number of modes in which the reconstructed time coefficient is zero for different TKE levels. In the same case, the number of corresponding modes of the forced state exceeds that of the unperturbed flow, and the reconstructed signal of the latter does not appear to be zero until the energy ratio reaches 90%. The explanation for this phenomenon is that the periodic perturbation injected into the flow field causes the energy to be redistributed, and the low-frequency components in each order mode transition to the high-frequency components.


**Table 2.** The number of modes with zero signal reconstructed at various turbulent energy.

Then, the global frequency information of the time series within the sampling period is solved using the Hilbert transform [31], which is defined as the convolution of the signal itself and the reciprocal of the time domain:

$$y(t) = x(t) \* \frac{1}{\pi t} \tag{13}$$

where *x*(*t*) is the real part of the complex signal, and *y*(*t*) stands for the imaginary part.

Next, the instantaneous phase *θ*(*t*) and instantaneous frequency *f*(*t*) of the analytical signal will be represented by the following equation:

$$\theta(t) = \arctan \frac{y(t)}{x(t)},\\ f(t) = \frac{1}{2\pi} \frac{d\theta(t)}{dt} \tag{14}$$

The instantaneous frequencies of the original signal and the reconstructed sequence in the time domain are indicated in Figure 14. The time coefficients of the two working states have abrupt changes in high- and low-frequency components at certain moments, which is a typical feature of non-stationary signals. In the disturbed state, the time proportion of highfrequency bands is much higher than that in the unforced case. In comparison to the original sequence, the instantaneous frequency mutation of the reconstructed signal is immensely shortened because of filtering out high-frequency noise data, and it is maintained at a lower level. Some useful high-frequency constituents are retained in the reconstructed signal of disturbed condition, embodying the characteristics of external disturbances.

**Figure 14.** Instantaneous frequency distribution: (**a**) None; (**b**) Control.

Up to now, the time coefficient of each order is reconstructed, and then the pulsating velocity fields will be determined by formula 5. In this section, the original flow field and the reconstructed flow field at a certain instant are selected for comparative analysis (see Figure 15). In the original flow field, there are large-scale structures of red and blue alternating flow direction, and small-scale pulsations are distributed near the edge. As for the reconstructed flow field, on the basis of retaining the large-scale motion, the small-scale structures are filtered out, and the shape of the flow structures is relatively smooth.

**Figure 15.** Comparison of (**a**) original and (**b**) reconstructed fluid field.

Figure 16 manifests the relative probability statistics of the dimensionless large-scale streamwise and wall-normal pulsation velocities over the whole viewing window. The ejection event in the second quadrant and the sweep event in the fourth quadrant are the main constituents of the turbulent events, and the intensity of the burst events is directly related to the Reynolds shear stress [32]. In contrast, the relative probability of large-scale streamwise direction pulsation in the disturbed flow field decreases, which can be confirmed by the decline in the flow direction length of the contour. The enhancement in the area proportion of the red region indicates that the flow direction and wall-normal fluctuation of the smaller scale are heightened. Because the synthetic jet has a weak impact on the large-scale wall-normal pulsation, the variation of the outermost contour of the cloud image in the wall-normal direction is not prominent.

**Figure 16.** Relative probability of flow and wall-normal fluctuation velocity: (**a**) None; (**b**) Control.

The detection of burst events can also be implemented with the spatial local average velocity structure function, overcoming the limitation of traditional condition sampling method, which can characterize the relative motion and local deformation of turbulent vortex structure at specific spatial scales. The local mean structure function of the large-scale streamwise pulsation velocity along the flow direction is as follows:

$$
\delta\mu\_x = \overline{u\_r'(\mathbf{x}\_0 : \mathbf{x}\_0 + l\_\mathbf{x}, y\_0)} - \overline{u\_r'(\mathbf{x}\_0 - l\_\mathbf{x} : \mathbf{x}\_0, y\_0)}\tag{15}
$$

where *lx* represents the spatial scale of the streamwise direction. Yang et al. [33] selected four scales for analysis in the experiment, namely, the corresponding length of the four grid points of the flow direction. Since the fluctuation velocity has been filtered above, eight scales are adopted to explore here. The positive and negative values of *δux* signify that the fluid micro clusters are stretched and compressed along the streamwise in a region with a flow direction scale of 2*lx*, corresponding to ejection and sweeping events, respectively. In this study, the shear layer is obtained by detecting large-scale burst events, and the detection function is written as [34]:

$$D(\mathbf{x}\_{0\prime}y\_0, l\_{\mathbf{x}}) = \begin{cases} \mathbf{Q}2\_{\prime} & \delta u\_{\mathbf{x}} < 0\_{\prime} \ |\delta u\_{\mathbf{x}}| = \max |\delta u\_{\mathbf{x}}|\\ \mathbf{Q}4\_{\prime} & \delta u\_{\mathbf{x}} > 0\_{\prime} \ |\delta u\_{\mathbf{x}}| = \max |\delta u\_{\mathbf{x}}| \end{cases} \tag{16}$$

Conditional sampling is still performed at the wall-normal height of 0.15*δ*, and largescale ejection event (blue area), sweep event (red area) and inclined shear layer (white dashed lines) in the near-wall region are identified, as shown in Figure 17. Compared with the undisturbed case, the large-scale fluctuation amplitude of the controlled flow field is reduced, the inclination of the shear layer is strengthened, and the extension length of the streamwise direction is narrowed, demonstrating the radiation range of the collision between high- and low-speed fluids is diminished, and the Reynolds shear stress is extensively weakened.

**Figure 17.** Conditional average results of shear layer: (**a**) None; (**b**) Control.

#### **4. Conclusions**

The velocity vector field of the interaction between underwater synthetic jet and TBL was measured utilizing the TR-PIV technology, and the influences of periodic perturbation on the multi-scale features of wall turbulence was analyzed by means of various research methods, and the following main conclusions were drawn:


periodic disturbance generated by the synthetic jet accelerates the migration of lowfrequency signals to high-frequency signals.


**Author Contributions:** Conceptualization, B.L; methodology, B.L.; formal analysis, B.L.; investigation, B.L. and J.Z.; data curation, B.L.; writing—original draft preparation, B.L.; writing—review and editing, N.J.; visualization, J.Z.; supervision, N.J.; funding acquisition, N.J. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the National Natural Science Foundation of China (Grant Nos. 11732010, 11972251, 11872272, 11902218, 12172242), Chinesisch-Deutsche Zentrum für Wissenschaftsförderung (Grant No. GZ1575) and the Ministry of Industry and Information Technology (Grant No. (2019)360).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Dual Synthetic Jets Actuator and Its Applications—Part I: PIV Measurements and Comparison to Synthetic Jet Actuator**

**Zhenbing Luo \*, Zhijie Zhao, Xiong Deng, Lin Wang and Zhixun Xia**

College of Aerospace Science, National University of Defense Technology, Changsha 410073, China; 18041360895@163.com (Z.Z.); badi\_arg@126.com (X.D.); wanglin-2007@126.com (L.W.); zxxia@nudt.edu.cn (Z.X.) **\*** Correspondence: luozhenbing@163.com

**Abstract:** In order to understand the differences between dual synthetic jets (DSJs) and synthetic jets (SJs), particle image velocimetry (PIV) technology is used to capture the basic flow field characteristics of a dual synthetic jet actuator (DSJA) and a synthetic jet actuator (SJA), and then a careful comparison between them is implemented. The results indicate that a cycle of the DSJ is divided into two stages. In the near-field downstream, a pair of synthetic jets entrain fluid around them and interact with each other, making the flow field complex, and the time-periodic diaphragm dominates them. There is an unfavorable phenomenon of "self-support" between the two jets. In the far-field downstream, the two jets merge into a single, more stable SJ with a higher velocity and a double characteristic frequency. The DSJs have also shown good vectoring characteristics, with the vectoring deflection angle (VDA) changing from about −46◦ to 46◦. The above results demonstrate that the DSJA may replace the traditional SJA in all kinds of applications and extend the applying area of the SJ to more active flow control systems, which cannot be qualified by traditional SJA.

**Keywords:** dual synthetic jets; synthetic jet; vectoring characteristics; PIV

Wang, L.; Xia, Z. Dual Synthetic Jets Actuator and Its Applications—Part I: PIV Measurements and Comparison to Synthetic Jet Actuator. *Actuators* **2022**, *11*, 205. https://doi.org/ 10.3390/act11080205

**Citation:** Luo, Z.; Zhao, Z.; Deng, X.;

Academic Editors: Luigi de Luca, Hui Tang and Xin Wen

Received: 23 June 2022 Accepted: 17 July 2022 Published: 22 July 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

Synthetic jets have been widely researched in the past twenty years, and flow control means based on synthetic jet actuators have become a critical tool for the aerodynamicist [1–5]. The SJA is a zero-mass-flux device that converts electrical energy into momentum. A schematic of a traditional SJA and the schlieren flow visualization are shown in Figure 1 [6]. The traditional SJA is composed of a slot and an enclosed cavity, bounded by one or two piezoelectric (PZT) diaphragms. A piezoelectric element can vibrate by applying an alternating voltage to it, and the outside air is periodically pushed out or sucked into the chamber through the orifice with a rapid response. The SJA is deemed to be a "zero-mass-flux" device, because the net mass flowing from the jet to the outside air is zero. The mass of the drawn air is equal to the mass of the ejected air. As an example: imagine that you are sitting on a frictionless surface. Make an "O" shape with your mouth, and then breathe in. You do not move. Breathe out through your mouth, keeping an "O" shape, and you will move in the opposite direction to which you are breathing, based on Newton's Third Law of Motion. If you repeat these processes several hundred or thousand times per second, you will be an SJ. In fact, the SJA, with its advantages of simple structure, low cost, easy operation, no need of an air source and pipes, could establish an SJ more rapidly than other actuators that generate a steady jet or a pulsed jet. The distinctive characteristics of the synthetic jet flow field and the good working attributes of actuators make the SJA suitable for a variety of applications, such as active flow control, flight control, and thermal management in microelectronic devices [6–17].

The traditional SJA, however, with one side of the PZT diaphragm exposed to the air drawn from the base flow (BF) and the other surface exposed to the environment flow (EF), will have two problems of pressure failure and inefficient energy utilization, as shown in Figure 1. The pressure loading will cause a smaller deflection of the PZT diaphragm, supporting the fact that the SJA must require a higher power to drive the PZT diaphragm when the pressure differential between the BF and the EF is large. Moreover, half of the vibrating energy of the diaphragm is wasted in the environment, meaning the energy utilization is not efficient.

**Figure 1.** Schematic of a traditional SJA and schlieren flow visualization. The two surfaces of the diaphragm are exposed to the BF and EF, respectively [6].

A novel dual synthetic jet actuator was invented by us [18]. A schematic and a photograph of the DSJA are shown in Figure 2. The DSJA is made up of two enclosed cavities that are bounded by a single PZT diaphragm, with two emitting slots and a slide block (SB) placed in the middle of the two slots. In the evolutionary process of dual synthetic jets, two synthetic jets out of phase are forced out of the cavity in every oscillating cycle, instead of a fluid puff for every vibrating circle from a traditional SJA. Actuators with a similar physical structure driven by a loudspeaker, called twin synthetic jet actuators, have also been displayed and researched in references [19,20]. The unique merits of the DSJAs lie in the fact that the two chambers share one wall equipped with a single PZT diaphragm, and the middle slide block regulates the two synthetic jets. The two cavities sharing the same PZT diaphragm make the novel synthetic jet actuator not only double the function of the traditional SJA, but also resolves the problems of pressure failure and inefficient energy utilization of the SJA. As the diaphragm of the DSJA is completely surrounded by the fluid drawn from the BF, the DSJA avoids the high power required to drive the PZT diaphragm when the pressure differential is great. In addition, the slide block can adjust the two synthetic jets, enabling the DSJA to have a unique thrust-vectoring characteristic that the traditional SJA does not have [21,22].

To understand the detailed differences between the DSJ and the SJ, PIV technology is used to capture the basic flow field characteristics of the DSJA and the SJA, and then careful comparison between each other is implemented. In addition, the unique vectoring characteristic of the DSJA is also discussed to show the regulating ability of the middle slide block.

**Figure 2.** Schematic and photograph of a novel DSJA. The two cavities share the same PZT diaphragm, and the diaphragm is only surrounded by the fluid drawn from the BF [18].

#### **2. Experimental Technique**

The DSJA studied in this paper is "mini" in size, and is composed of two symmetrical, cylindrical chambers (inner diameter *D* = 46 mm, height *H* = 7 mm) bounded on one end by a PZT diaphragm (thickness = 0.2 mm) and driven by an electrical signal with voltage amplitude *UA* of 300 V and a frequency *f* of 500 Hz. A rigid wall (thickness *ht* = 4 mm) with two symmetrical slots (length *l* = 20 mm, width *h* = 2 mm) covers the upper end, and the distance between the two slots *d* is equal to 5 mm. The driving frequency is equal to the diaphragm resonance frequency, which can generate a higher velocity than that actuated at the Helmholtz frequency (1196.5 Hz). Figure 3 shows the overall structures and physical models of the SJA and the DSJA. It is worth noting that the left slot of the DSJA will be sealed to turn into a SJA under the test condition of the SJA. The key parameters of both modes are exhibited in Table 1.

**Figure 3.** Structure diagram of SJA and DSJA.


**Table 1.** Geometry and electrical parameters of SJA and DSJA.

PIV is an ideal technique that can characterize non-periodic and quasi-periodic phenomena temporally, but accurate phase referencing is not available. For fast-varying flow fields with a fine spatial structure, traditional PIV approaches cannot capture the short time and length scales of the flow field and its temporal–spatial evolution. In practical applications, it is necessary to address issues such as the transient response of both the fluid and the structure, the energy transfer process among the different scales in transition, turbulence, and so on. Considering that the SJ has a periodic unsteady flow and phase referencing rules are available, the transfer phase and sub-frequency (TPSF) technique is applied to capture the arbitrary phase of the DSJ, and the phase of the DSJ can be determined by the transfer-phase-to-equal (TPE) technique in previous work by us [14].

The PIV system comprises a twin-cavity laser (2 × 200 mJ), a light guide arm, a chargecoupled-device (CCD, 2456 × 2056 pixels, 12 bits) camera with a 24 mm F/2.8 lens, a narrow band-pass filter of 532 ± 5 nm, a sync controller and a PC software, as shown in Figure 4. The DSJA is fixed in a closed glass chamber with a size of 400 × 200 × 200 mm. The smoke, with an average diameter of less than 4μm, is adopted as tracer particles. A signal generator can provide two signals with the same phase and frequency. One is transmitted to the driving device of the DSJA, and the other one is provided to the sync controller as the external trigger signal, by which the locked phase measurements can be realized. The laser sheet, projected perpendicular to the slots at the midsection for x-y plane measurements, has a thickness of 0.5 mm and is employed to illuminate the two-dimensional flow fields. The sampling frequency is set at 10 Hz. The uncertainty of the tracer-particle displacement is less than 0.1 pixels. Therefore, the uncertainty of the measured velocity would be less than 0.7 m/s. A DSJ cycle is uniformly divided into 12 phase points. A PIV recording is implemented using TPSF and TPE technology, and then a phase-averaged velocity field is acquired based on 50 pairs of captured images at each phase point. Figure 3 also shows the DSJA coordinates. The flat plate containing the slots is parallel to the free surface of the water, and the jets are emitted into the fluid medium perpendicular to the plate. The camera is focused onto an area of 25 mm (horizontal) × 25 mm (vertical), with the lower border of the image coinciding with the flat plate.

**Figure 4.** Diagram of PIV system.

#### **3. Results and Discussion**

#### *3.1. Comparison between DSJA and SJA*

The tests of a novel DSJA and a traditional SJA at the same forcing conditions and some key parameters of the SJA are listed in Table 2.

**Table 2.** Key parameters: Comparison of a novel DSJA and a traditional SJA.


The Reynolds number (*Reu*) and the Strouhal number (*St*) of the SJ are defined as:

$$\text{Re}\_{\text{ul}} = \rho h u\_{amp} / \mu \tag{1}$$

$$St = hf/u\_{amp} \tag{2}$$

where *ρ* and *μ* is the air density and viscosity, respectively, and *uamp* is the velocity amplitude of the SJ.

The stroke length *L*<sup>0</sup> of the SJ and *Re* based on the blowing phase per unit width *ReI0* were provided by Smith and Glezer as follows [6]:

$$L\_0 = \int\_0^{T/2} u\_0(t)dt\tag{3}$$

$$\text{Re}\_{\text{I}\_{\text{0}}} = \text{I}\_{\text{0}} / \mu\text{h}, \text{ I}\_{\text{I}} = \rho\text{h} \int\_{0}^{T/2} u\_{\text{0}}^{2}(t)dt \tag{4}$$

where *u*0(*t*) is the instantaneous centerline velocity at the actuator slot exit and *T* is the time of a DSJ period. The dimensionless time *t \** is equal to *t/T-N* and *T = 1/f*, where *N* is the circles.

Figures 5 and 6 show the flow fields of the SJA and the DSJA at a different dimensionless time in an actuating cycle. Figure 7 shows time-locked velocity traces of one cycle at six measured points ((1.5, −3.5), (1.5, 0), (1.5, 3.5), (6, 0), (10, 0), and (20, 0)), and the points (1.5, −3.5) and (1.5, 3.5) are, respectively, located on the left and right slot exit centerline of the DSJA.

Figure 5 shows the PIV measurements of the SJA for comparison. Figures 6 and 7 demonstrate a cycle of the DSJ that can be divided into a "left" stage and a "right" stage. In the "left" stage, the flow field is dominated by the left jet (Jet1). Similarly, the flow field is dominated by the right jet (Jet2) in the "right" stage. This means that the time-periodic PZT diaphragm determines the flow field of the DSJ. It is worth noting that the two synthetic jets out of the phase entrain the air around them, and then interact with each other in the near-field downstream, as shown in Figure 6. When the Jet1 is in the suction stroke and the Jet2 is in the blowing stroke (Figure 6g–l), some fluid of the Jet2 is sucked into the left slot, while some fluid of the Jet1 is sucked into the right slot when Jet1 maintains the blowing stroke (Figure 6a–f), which supports that a phenomenon of "self-support" appears between the two synthetic jets. If all of the fluid of the jet is sucked into the adjacent slot, there will not be a stable jet in the far-field downstream of the DSJA. Therefore, the phenomenon of "self-support" between a pair of synthetic jets is unfavorable, which could weaken the strength of the dual synthetic jets downstream.

**Figure 5.** PIV measurements of the SJA.

**Figure 6.** PIV measurements of the DSJA.

**Figure 7.** Time-locked velocity traces on six measured points per unit circle.

Figure 6 also shows the vorticity magnitude and two counter-rotating vortex pairs produced by the DSJA. The vorticity of the vortex pairs contributes to the pressure difference. An "induced" velocity *V*, a translation velocity with which the vortex ring travels, is given by Hill and Saffman [15] as follows:

$$V = \frac{\Gamma}{2\pi a} \left[ \ln \frac{8a}{\left(4\nu t\right)^{1/2}} - 0.558 + O\left(\frac{\nu t}{a^2}\right)^{1/2} \right] \tag{5}$$

where Γ is the circulation and *a* is the radius of the vortex ring. The vortices decelerate not only by losing energy due to friction, but also through entrainment of fluid from the environment. The circulation Γ ∝ *ω* and the "induced" velocity *V* contributes to the pressure difference. The pressure in the area of the main vortex pairs is lower than in other flow domains. The stronger the vortex *ω*, the stronger the circulation Γ, the higher the velocity *V*, and the lower the pressure will get. Equation (5) also shows that the velocity *V* is strongly associated with the time *t*, indicating that the phase angle also contributes to the "induced" velocity *V* and the pressure difference. The results show that in the far field downstream, the fore vortex pairs entrain the aft vortex pairs, and the aft vortex pairs are deflected to the fore vortex pairs due to the lower pressure. Meanwhile, the fore vortex pairs will decelerate gradually by losing energy due to friction and entrainment of fluid from the environment, and the aft vortex pairs will also impact the fore vortex pairs. Lastly, a single, more stable synthetic jet could be realized by the merging of the two synthetic jets. Figure 7 also indicates that the two jets of the DSJA could merge into a single, more stable synthetic jet in the far field. Furthermore, Figure 7b,c shows that the DSJ downstream has a double characteristic frequency, which broadens the applying frequency domain.

The detailed comparisons between the DSJ and the SJ are shown in Figures 8–10. In the near-field downstream, the flow fields controlled by the DSJ and the SJ similarly emerge into four stages (accelerating blow, decelerating blow, accelerating suction and decelerating suction), while the flow characteristic of the DSJ is more complex than that of the SJ. Moreover, the blowing or sucking peak velocity of the DSJ is larger, showing a higher control ability. It is worth noting that in the mode of the DSJ, the left peak velocity is slightly lower than the peak velocity of the right exit, which may result from the preload differences of both sides. Considering that the PZT diaphragm is fixed by bolts, the clamping conditions on both sides of the diaphragm may be different, resulting in the inconsistent vibrating amplitude in both directions of the PZT diaphragm, thus showing the different velocity evolution. Furthermore, the two synthetic jets of the DSJ could merge into a single, more stable jet with a higher velocity in the far field downstream.

**Figure 8.** Streamline velocity of DSJ and SJ on x = 1 mm at different times.

**Figure 9.** Streamline velocity of DSJ and SJ on y = 0 at different times.

**Figure 10.** Mean velocity of DSJ (**left**) and SJ (**right**).

Based on the above analysis, the structure of the two cavities sharing a single diaphragm ensures that the DSJA not only has the unique property of zero mass flux, but also solves the problems of pressure failure and inefficient energy utilization existing in a traditional SJA, which extends the applications of the DSJA to more flow control systems, such as thrust vectoring control and flight control [23,24].

#### *3.2. Special Vectoring Characteristics of DSJA*

The slide block between the two slots can regulate the two synthetic jets by removing it left or right, making the DSJA have a unique vectoring characteristic. The DSJ are formed by the fusion of two SJs, and the momentum ratio and low-pressure zone strength of the two jets can be modulated by changing the area ratio of two slots, which indicates the vectoring deflection of the DSJ. The low-pressure zone and the inertia of the two SJs are the key parameters that demonstrate the vectoring deflection angle (VDA) of the DSJA [20]. Furthermore, it is the structure of the middle slide block that also determines the area and strength of the low-pressure zone, playing an important role in the vectoring characteristics. Based on the method proposed by Deng [19] to calculate the VDA, the PIV measurements have been carried out to show the vectoring characteristic of the DSJA, with different slot area ratios and changed structures of the slide block. The structure schematic of the slide block is shown in Figure 11, where the exit chamfer α and the length of the left exit *dL* are variable. In detail, α could be changed from 0◦ to 60◦ by transforming the slide block, and *dL* can be shifted from 0.4 mm to 3.6 mm with an interval of 0.4 mm. Additionally, the driving voltage (±250 V) and the frequency (550 Hz) are maintained constantly in this research.

**Figure 11.** Structure schematic of slide block.

The PIV test results have been shown in Figure 12, suggesting that the DSJ can deflect to the wider exit and the VDA increases gradually with the augmentation of *α* for the same *dL*. In addition, Figure 13 shows the VDA of the DSJ under different *α* and *dL*, obviously indicating that for a different α, the VDA basically presents an antisymmetric change relative to the point of *dL* = 2 mm. When α = 0◦ or 45◦, the VDA increases first and then decreases on the single side, and the extreme values are obtained when *dL* = 1.2 mm and 2.8 mm, while when α = 60◦, the VDA keeps decreasing. The bigger the *α*, the larger the VDA, which is consistent with the change trend of the PIV results. The larger *α* could generate a stronger Coanda effect, making the jet with a shorter exit width deflect to a greater angle, thus realizing a bigger VDA, although the jet with a larger exit could suppress the deflection. Due to the positioning error of the slide block, the two outlets are not completely symmetrical when *dL* = 2 mm, resulting in a certain deviation between the VDA and the theoretical value (0◦). It is worth noting that when α = 60◦, the VDA shows a linear relationship with the increase in *dL*, and the approximate variation range of the VDA (−46◦~46◦) is wider, which is more conducive to the design and optimization of the active flow control law.

**Figure 12.** PIV measurements of DSJA with different α and *dL*.

**Figure 13.** VDA measurements of DSJA with different α and *dL*.

#### **4. Conclusions**

To understand the detailed discrepancy between the DSJ and the SJ in quiescent surroundings, a novel DSJA was investigated using PIV, and its flow-field characteristics were compared with the traditional SJA. A transfer-phase and sub-frequency technique was introduced to capture the arbitrary phase of the DSJ, and a transfer-phase-to-equal technique was provided to determine the phase of the DSJ. In addition, the unique vectoring characteristics of the DSJA with a variable *α* and *dL* were also explored.

The results indicate that a cycle of the DSJ can be divided into a "left" stage and a "right" stage. In the "left" stage, the flow field is dominated by the left jet, while the flow field is dominated by the right jet in the "right" stage. In the near-field downstream of the DSJA, the two jets entrain air around them and interact with each other, and the flow field is more complex than that of the SJA. There is an unfavorable phenomenon of "self-support" between the two synthetic jets, weakening the strength of the DSJ downstream. In the-far field downstream of the DSJA, the two vortex pairs interact with each other and entrain fluid from the surroundings, and the DSJs merge into a single, more stable synthetic jet (similar to a steady jet) with a higher velocity and a double characteristic frequency. In addition, the DSJA has shown great vectoring characteristics, with the VDA changing from approximately −46◦ to 46◦, suggesting a broader applying prospect. These results indicate that the traditional SJA may be replaced by the DSJA in many applications, and the applying area of the synthetic jets could be extended to more active flow control systems that cannot be qualified by the traditional SJA.

In this paper, the comparison of the DSJ and the SJ in quiescent environments is explored, but its evolution process under the condition of incoming flow was not shown, which may have more engineering significance. Therefore, further work will focus on the comparison of the evolution process with the incoming flow, and wing tunnel tests will also be implemented.

**Author Contributions:** Conceptualization, Z.L. and Z.X.; methodology, Z.Z. and L.W.; formal analysis, Z.Z. and Z.L.; writing—original draft preparation, X.D. and Z.Z.; writing—review and editing, Z.L. and Z.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** The authors are grateful to the National Nature Science Foundation of China for its financial support for this project (Grants: U2141252, 11972369, 11872374).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclature**


#### **References**


## *Article* **Dual Synthetic Jets Actuator and Its Applications—Part II: Novel Fluidic Thrust-Vectoring Method Based on Dual Synthetic Jets Actuator**

**Jie-Fu Liu, Zhen-Bing Luo \*, Xiong Deng, Zhi-Jie Zhao, Shi-Qing Li, Qiang Liu and Yin-Xin Zhu**

College of Aerospace Science and Engineering, National University of Defense Technology, Changsha 410073, China; 18800100128@126.com (J.-F.L.); badi\_arg@126.com (X.D.); 18041360895@163.com (Z.-J.Z.); lishiqing18@nudt.edu.cn (S.-Q.L.); liuqiang12@nudt.edu.cn (Q.L.); yinxin.zhu@nudt.edu.cn (Y.-X.Z.) **\*** Correspondence: luozhenbing@163.com

**Abstract:** A novel fluidic thrust-vectoring (FTV) control method based on dual synthetic jets actuator (DSJA) is proposed and evaluated. Numerical simulations are governed by the compressible Unsteady Reynolds-Averaged Navier–Stokes (URANS) equations. According to the results, DSJA is capable of deflecting a primary jet with a velocity of 100 m/s and a height of 50 mm by approximately 18 degrees with a momentum coefficient of 1.96%. It produces comparatively linear control characteristics in almost all deflection angles evaluated (0~23 degrees). The low pressure generated by DSJA, the ejecting enhanced by DSJA, and the co-flow effect produced by the accelerated secondary jet all play roles in the deflection of the primary jet. Since the primary jet is strong enough, the potential of DSJA to provide thrust vector control is revealed.

**Keywords:** fluidic thrust vectoring; dual synthetic jets actuator; active flow control

### X.; Zhao, Z.-J.; Li, S.-Q.; Liu, Q.; Zhu, Y.-X. Dual Synthetic Jets Actuator and Its Applications—Part II: Novel Fluidic Thrust-Vectoring Method Based on Dual Synthetic Jets Actuator. *Actuators* **2022**, *11*, 209. https://doi.org/10.3390/

**Citation:** Liu, J.-F.; Luo, Z.-B.; Deng,

Academic Editor: Luigi de Luca

Received: 16 June 2022 Accepted: 18 July 2022 Published: 29 July 2022

act11080209

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

Fluidic thrust-vectoring (FTV) technology is a kind of thrust-vectoring technology developed to control the direction of thrust without using the conventional moving surfaces, so as to reduce the structural complexity, weight, and volume, and maintaining the cost of the thrust vector control nozzle. Numerous researchers have investigated the use of FTV systems for aircraft control, and many FTV methods have been developed.

Until now, the most researched idea is to change the flow state in nozzles through the proper set of high-speed secondary jets (or suction) with designed nozzle shape, such as the shock-vector control method [1], throat-shifting method [2], dual-throat method [3] and counterflow method [4]. The co-flow method has been investigated and verified by experiments to control the subsonic and supersonic primary jets [5,6]. The co-flow method has been proven capable of realizing flight control by the flight test of MAGMA aircraft [7]. However, this method requires a high-pressure air source. Carrying an air supply would result in the massive occupation of fuselage space and weight. A high-pressure air source could also be supplied from power plant bleeding, leading to an impact on the condition of the engine. The mass of MAGMA aircraft's flow control module occupies 40% of the total system [7] (providing circle control and FTV control). The pressure supply and the supporting pipeline system cause a large space and weight burden for aircraft.

Wen [8] strengthened the co-flow FTV control method by upgrading the secondary co-flow from a conventional steady jet to a sweeping jet. The sweeping jet was found to be more capable of enhancing the flow mixing than the steady jet, which leads to a higher efficiency in co-flow vectoring control. Using a sweeping jet greatly reduces the mass flow rate of the secondary jet to provide thrust vector control, and hopefully will alleviate the burden brought by a high-pressure gas source to a certain extent.

A bypass dual-throat nozzle [9] is developed from conventional dual-throat methods. This introduces secondary flows from a bypass set between the upstream minimum area

and the upstream convergent section instead of air supply. This improvement reduces the layout difficulty of using a dual-throat vectoring nozzle and achieves an equivalent or even better performance than the conventional dual-throat vectoring nozzle. In this study, the realizable *k* − *ε* turbulent model with the standard wall function is selected.

Meanwhile, predecessors have developed some FTV methods with no dependence on air source, called passive FTV methods. Using an ejecting-mixing fluidic vectoring nozzle [10,11] takes advantage of the entrainment effect of the primary jet on the surrounding air. In detail, the deflection angle of the primary jet could be controlled by adjusting the mass rate of the entrainment secondary jet through the control valves. The results of a supersonic primary jet numerical simulation [10], a flight test powered by electric duct fan [12], and a 6 kg turbojet engine ground experiment [13] show that this method is suitable to control the direction of subsonic and supersonic primary jets.

As one of the lightest active-flow control means, synthetic jet (SJ) actuators have been widely applied in many fields, such as separation control [14], vortex-induced vibration control [15], and flight control [16]. Researchers have found that with the help of zeromass-flux synthetic jets, the direction of the primary jet can be deflected even without any moving components. SJ actuators can be driven by piezoelectric (PZT) diaphragms or electrodes fixed in a cavity, capable of switching work state rapidly in response to different electrical signals.

Smith [17,18] found that a synthetic jet actuator installed close to the primary jet outlet could deflect the primary jet with a low momentum coefficient. However, the deflection becomes undetectable when the velocity of the primary jet increases to 30 m/s. Pack [19] carried out an experimental study with a diffuser attached to the jet exit, and the excitation was introduced at the junction between the jet exit and the diffuser inlet, demonstrating that the presence of the wide-angle diffuser increases the effectiveness of the added periodic momentum caused by a favorable interaction among the excitation, the jet shear layer, and the diffuser wall. However, Pack [19] and Chiekh [20] found that the direct deployment of the diffuser wall will lead to no-control primary flow into a bistable situation, which is intolerable in an FTV system as it will result in an unstable control moment.

The PZT diaphragm on the conventional SJ actuators is set between different pressure conditions of basic flow (the flow field to be controlled) and environment flow (usually the equipment cabin), which indicates that half of the radiation energy of the vibrating diaphragm is injected into somewhere unwanted. It also means that the PZT must overcome the pressure difference on two sides of the diaphragm to vibrate. When the pressure differential is large, the driving vibration will be difficult or even unachievable. In short, the structural design of conventional PZT-driven SJ actuators limits their energy efficiency and environmental adaptability. Subsequently, Luo [21] invented a new-generation synthetic jet actuator, known as a dual synthetic jet actuator (DSJA). Different from conventional SJ actuators, DSJ actuators have two cavities on both sides of the PZT diaphragm. Both cavities are connected to the controlled flow field, which guarantees that the pressure difference between the two sides of the PZT diaphragm will not be too large. Therefore, the energy efficiency and environmental adaptability of DSJ actuators are improved. Moreover, DSJ actuators have greater flexibility and stronger control ability. For example, DSJ could perform both functions of "push" and "pull" to vector a primary jet [22]. However, according to the latest explorations, DSJ actuators or SJ actuators are only capable of deflecting a primary jet with low velocity and small height.

To sum up, plenty of significant studies have investigated the feasibility of synthetic jet application to FTV. As a control method without introducing an air source, FTV by DSJ is a technology worthy of further in-depth research and development. However, existing studies have shown that the simple deployment of DSJ is insufficient to provide effective FTV control. This work explores the potential of dual synthetic jets for thrust vector control by proposing a novel method and its configuration.

#### **2. Materials and Methods**

#### *2.1. Method and Configuration*

The configuration of this novel FTV method is shown in Figure 1. The new method was created from the reference and improvement of existing DSJ actuators jet-deflecting method and co-flow FTV nozzle. The height of the primary jet inlet is 50 mm. Avoiding the problems of pressure loading and energy inefficiency, DSJ actuators were selected (comparing to SJ actuators) as the provider of zero-mass jet to adapt to the changing external environment in flight. In this configuration, the DSJ actuators are installed on the edge of the primary nozzle outlet, with a total thickness of 10 mm. To prevent the unexpected deflection from damaging the practicability of the nozzle, passive secondary jet inlets are arranged as gaps between the primary jet inlet wall and the nozzle wall. The gaps are connected to the outside atmosphere and are approximately 8 mm high at the downstream edge. The exits of DSJ actuators are arranged at the two downstream corners of the edge and are set 45 degrees up and down from the mainstream direction, respectively.

**Figure 1.** Configuration design (**a**) and numerical calculation mesh configuration (**b**).

#### *2.2. Simulation Setup*

To explore flow details in nozzles and the control mechanism, 2D-simulations were applied by ANSYS Fluent. The flow is governed by the compressible Unsteady Reynolds-Averaged Navier–Stokes equations (URANS):

$$\frac{\partial \rho}{\partial t} + \frac{\partial}{\partial \varkappa\_i} (\rho u\_i) = 0 \tag{1}$$

$$\frac{\partial}{\partial t}(\rho u\_i) + \frac{\partial}{\partial x\_j}(\rho u\_i u\_j) = -\frac{\partial p}{\partial x\_i} + \frac{\partial}{\partial x\_j} \left[ \mu \left( \frac{\partial u\_i}{\partial x\_j} + \frac{\partial u\_j}{\partial x\_i} - \frac{2}{3} \delta\_{ij} \frac{\partial u\_l}{\partial x\_l} \right) \right] + \frac{\partial}{\partial x\_j} \left( -\rho \overline{u\_i' u\_j'} \right) \tag{2}$$

where the density *ρ* is defined by the ideal gas equation:

$$P = \rho RT\tag{3}$$

and the molecular viscosity *μ* is assumed to be constant. In addition, the gravity effect is neglected.

Considering the complex flow conditions such as separation, attachment, and vortex, the realizable *k* − *ε* turbulent model (RKE) is selected to model the Reynolds stresses −*ρu i u j* . This model is also applied in the studies conducted by Gu [9], in which the turbulent viscosity is calculated by using an improved method. The exact transport equation of the fluctuating component vorticity is used to derive the dissipation rate equation.

Non-Equilibrium Wall Function is used to compute complex flows such as separation, reattachment, and impact near the wall. It further extends the applicability of the wall function approach by including the effects of pressure gradient and is found to accurately predict the reattachment length, skin friction, and static pressure coefficient with the RKE model [23].

Semi-Implicit Method for Pressure-Linked Equation (SIMPLE) is used in the coupling of velocity and pressure. The spatial discretization is second order for pressure and second order upwind for density, momentum, and energy.

The configuration of mesh and boundary conditions are shown in Figure 1. Wall y+ is approximately 1. No slip wall is applied for the nozzle and actuators' wall, pressure inlet for the primary jet, or passive secondary jet. As to environment, pressure inlet is applied for upstream boundaries, and pressure outlet for the side and downstream boundaries.

The initial condition is the transient simulation result of the baseline condition (where the primary jet is present alone and dual synthetic jets are off). All the flow field results shown were captured at the time when dual synthetic jets present for 0.1 s on the baseline of the initial condition. Applied first-order implicit time-integration schemes were used for temporal discretization. A timestep of 5 × <sup>10</sup>−<sup>5</sup> s (1/40 of the DSJ driving cycle) was applied for all transient simulations. The maximum iteration of each timestep is 50 and the convergence criterion is 1 × <sup>10</sup>−<sup>6</sup> for continuity and energy. In addition, DSJ exits were set as velocity inlet boundaries to simulate the periodic jets. Referring to the driven model of synthetic jets [24], the velocity is defined by the UDF formula shown as follows:

$$V(t) = V\_{\dot{\jmath}} \times \sin(2\pi t/T) \tag{4}$$

where *Vj* is the peak velocity of synethtic jets and *T* denotes the cycle of DSJ.

Time averaging is used to make the concerned controlled-flow field more intuitive. The specific process is to average the data of the flow field of 40 time-steps (corresponding to flow time of 0.002 s, single actuators cycle).

Three sets of meshes with grid quantities from 0.13 million to 0.21 million were evaluated for grid independence verification, with the deflection angle selected as the judgment index. As shown in Figure 2, the mesh with grid quantities of 0.18 million is competent for this simulation with sufficient accuracy.

**Figure 2.** Computed results of meshes with different grid amounts.

To validate the accuracy of the numerical method, simulations and experiments were conducted and compared. Limited to the experimental conditions, the primary jet velocity *Vp* is 34 m/s. The nozzle time-average pressure coefficient *Cp* distributed along the upper nozzle wall is selected as the judgment index when the deflection angle is 13 degrees. As shown in Figure 3, the maximum difference between the experiment and simulation result is less than 0.02. It indicates that the numerical method applied can accurately predict the actual flow situation.

**Figure 3.** Comparison of time-average pressure coefficient *Cp* along the upper wall.

#### *2.3. Relevant Definitions*

The method applied to determine the deflection angle *δ* is the same as all evaluated simulation cases in this work. In detail, we selected the point with the maximum velocity on the section of *x* = 0.4 m based on time-average flow field, and marked it as point A. The included angle between velocity direction at point A and *x* direction is defined as the deflection angle. Thus, the deflection angle *δ* can be defined as:

$$\delta = \tan^{-1}(\frac{V\_{Ax}}{V\_{Ay}}) \tag{5}$$

where *VAx* denotes *x* velocity component at point A and *VAy* denotes the *y* velocity component.

Momentum coefficient *C<sup>μ</sup>* is an important dimensionless parameter to evaluate the efficiency of active control means. For zero-mass periodic jets with sinusoidal variation (including DSJ), considering air as an incompressible gas, the momentum coefficient can be defined as:

$$\mathcal{C}\_{\mu} = \frac{n \times V\_{j}^{2} h\_{j}}{V\_{p}^{2} h\_{p}} \tag{6}$$

where *Vj* is the peak velocity of synthetic jets, *hj* is the height of dual synthetic jet exits, *Vp* is the velocity of the primary jet and *hp* denotes the height of the primary jet outlet. *n* denotes the number of zero-mass synthetic jets activated for control and *n* = 2 for a single DSJ actuator.

#### **3. Results**

#### *3.1. Control Characteristics of Thrust Vectoring*

To illustrate the effectiveness of DSJ in this configuration, Figure 4 shows the velocity contours of the flow field under different conditions. The primary jet velocity *Vp* is 100 m/s (based on the height of primary jet inlet, corresponding to Re = 3.9 × 104). The DSJ peak velocity *Vj* in Figure 5b is set as 70 m/s and its frequency *f* is 500 Hz (based on the height of primary jet inlet, corresponding to St = 0.25).

**Figure 4.** Velocity contours of flow field with no control (**a**), co-flow effect (**b**) and DSJ control (**c**) (*Vp* = 100 m/s).

The activation of DSJ makes the primary jet deflect upward; meanwhile, the secondary jet is accelerated. The primary jet finally gets a stable vector downstream. The primary jet is vectored upside with an eventual mean angle of approximately 18 degrees. The secondary jet will attain a mass flow of 0.541 kg/s by the ejection of the vectored primary jet and DSJ. As a contrast, injecting a secondary jet with the same mass flow can lead to a deflection of 4 degrees alone by the co-flow effect, as shown in Figure 4b. In this configuration, the properties of DSJ are critical in the deflection of the primary jet. Meanwhile, through the ejection effect enhanced by DSJ, the deflected primary jet activates the secondary jet to be a co-flow, offering an improvement in deflection performance.

To investigate the control characteristics of this configuration, simulations are implemented with different DSJ peak velocities (10 m/s, 20 m/s, ... , 100 m/s). The velocity distribution, deflection angles, and velocity distribution, are shown and analyzed as follows.

Figure 5 shows the time-average velocity distribution of controlled flow with different DSJ peak velocities. The DSJ frequency *f* is maintained at 500 Hz in this work. As shown in the velocity contours, the vectoring becomes observable when peak velocity increases to 30 m/s and increases as DSJ peak velocity increases.

**Figure 5.** Time-average velocity contours of controlled flow with DSJ of different intensities. (**a**) *Vj* = 10 m/s ( *C<sup>μ</sup>* =0.04%), (**b**) *Vj* = 30 m/s ( *C<sup>μ</sup>* =0.36%), (**c**) *Vj* = 50 m/s ( *C<sup>μ</sup>* =1.00%), (**d**) *Vj* = 70 m/s ( *C<sup>μ</sup>* =1.96%), (**e**) *Vj* = 90 m/s ( *C<sup>μ</sup>* =3.24%), (**f**) *Vj* = 100 m/s ( *C<sup>μ</sup>* =4.00%).

Figure 6 shows the variation of the deflection angle with DSJ peak velocity and the momentum coefficient. To most FTV means, the changing law between the deflection angle and DSJ peak velocity was supposed to be characterized by a dead zone, a linear control region, and saturation (according to Mason [5]). In this configuration, the establishment of deflection occurs when DSJ peak velocity is less than 20 m/s, and the control characteristic remains linear with peak velocity until the deflection angle reaches 23 degrees. Considering that the geometric exit angle of the nozzle is 28 degrees (which limits the maximum of the deflection angle), it shows a comparatively linear control characteristics in almost all control ranges. The highest efficiency occurs at *Cμ*, ranging from 0 to 0.5%, and the deflection angle reaches 18 degrees at *Cμ* of 1.96%. This indicates that this configuration has good efficiency performance in FTV control.

**Figure 6.** Variation of deflection angle with DSJ peak velocity (**a**) and DSJ momentum coefficient (**b**).

#### *3.2. Pressure Characteristics of Flow Field*

This section investigates the controlled flow in pressure characteristics. The timeaverage pressure distributions of the controlled flow with different DSJ peak velocities are shown in Figure 7.

**Figure 7.** Time-average pressure contours of the controlled flow with DSJ of different intensities. (**a**) *Vj* = 10 m/s ( *C<sup>μ</sup>* =0.04%), (**b**) *Vj* = 30 m/s ( *C<sup>μ</sup>* =0.36%), (**c**) *Vj* = 50 m/s ( *C<sup>μ</sup>* =1.00%), (**d**) *Vj* = 70 m/s ( *C<sup>μ</sup>* =1.96%), (**e**) *Vj* = 90 m/s (*C<sup>μ</sup>* =3.24%), (**f**) *Vj* = 100 m/s ( *C<sup>μ</sup>* =4.00%).

Based on time-average pressure, the iconic flow field structure is the low-pressure zone formed by DSJ. There are two cores of the low-pressure zone, one near the DSJ actuator exits and the other located at the wall of the throat. The greater the DSJ peak velocity, the larger the low-pressure zone. In addition, DSJ also reduces the pressure in the primary jet outlet. In Figure 8, the upper surface and outlet section are selected to visually compare the pressure change under different peak velocities.

Corresponding to two low-pressure core areas, two peaks of low-pressure zones along the nozzle wall are shown in Figure 7. Peak 1 is close to the exits of DSJ, directly caused by DSJ. With the augmentation of DSJ peak velocity, the growth rate of the low pressure in peak 1 decreases gradually, and the position of peak 1 moves downstream slightly. Peak 2 of the low-pressure zone appears on the wall of the throat, further downstream of peak 1.

With small DSJ velocity (below 50 m/s), the augmentation of peak 1 is significant, while the primary jet is far from the wall and peak 2 is not yet established. When the DSJ peak velocity is above 60 m/s, peak 2 begins to be established with an obvious augmentation while the augmentation of peak 1 is descending. This shows the significant effect of the wall on deflection when the primary jet is close enough to the wall. Eventually, the establishment of a low-pressure zone maintains a relatively stable level in all deflection angles. This indicates that the direct effects of DSJ (corresponding to peak 1) and the evolution of flow

of the controlled primary jet (corresponding to peak 2) provide comparable contributions in deflection.

**Figure 8.** Time-average pressure coefficient *Cp* distributed along the upper wall.

#### *3.3. Periodic Evolution of Flow Field*

To investigate the evolution of controlled flow, the characteristics of flow field in a DSJ cycle are explored. Figure 9 shows the velocity and pressure distribution of the flow field near the outlet at four typical moments (0T, T/4, T/2, 3T/4) of the DSJ cycle, starting from 0.1 s. To control flow, at the T/4 moment, the upper exit of DSJ is blowing at peak velocity, while the lower exit is sucking at peak velocity. Parts of the primary jet are absorbed in the DSJ cavity, and a strong low-pressure zone is formed downstream, which causes a substantial deflection of the primary jet. The blowing at the upper exit does not have a significant impact at this stage. At the T/2 moment, the upper exit finishes blowing and is going to absorb, while the lower exit is in the opposite situation, and the velocity of both exit inlets is almost zero. The primary jet is deflected at the previous moment develop to attach to the wall and the low-pressure zone generated at previous moment extends in the primary jet tunnel. At the 3T/4 moment, the upper exit of DSJ is absorbing at peak velocity, while the lower exit is blowing at peak velocity. The primary jet deflected at the moment of T/4 begins to attach to the wall. The pressure rises at the attachment area as the wall blocks the y-direction velocity of the primary jet, indicating that the primary jet is continuously deflecting when passing through the low-pressure zone. The low-pressure zone extends and moves downstream. The secondary jet is accelerated by the ejecting of the attached primary jet and the absorbing of the upper exit. An adverse pressure gradient occurs between the low-pressure zone and high-pressure attachment, and a separation vortex occurs at the throat. At the 0T moment, the upper exit finishes absorbing and is going to blow, while the lower exit is in the opposite situation. The separation vortex moves downstream. Throughout the cycle of flow-field evolution, the low-pressure generated by DSJ, the ejecting enhanced by DSJ, and the co-flow effect produced by the accelerated secondary jet all seem to play roles in the deflection of the primary jet.

**Figure 9.** *Cont*.

**Figure 9.** Periodic variation of velocity and pressure near the DSJ (*Vj* = 70 m/s). (**a**) T/4 velocity, (**b**) T/2 velocity, (**c**) 3T/4 velocity, (**d**) 0T velocity, (**e**) T/4 pressure, (**f**) T/2 pressure, (**g**) 3T/4 pressure, (**h**) 0T pressure.

#### **4. Conclusions**

A novel method of thrust vectoring based on DSJ and its configuration are proposed and evaluated by simulations. Different from previous schemes, it has two passive secondary flow outlets to stabilize the primary jet, and a finely designed nozzle wall geometric shape to optimize deflection capability. Numerical results reveal that this configuration has a high efficiency in deflecting a primary jet with a velocity of 100 m/s and a height of 50 mm and shows comparatively linear control characteristics in almost all control ranges.

For the momentum coefficient *Cμ*, the highest efficiency occurs at *Cμ*, ranging from 0 to 0.5%, and the deflection angle reaches approximate 18 degrees at a *Cμ* of 1.96%. The low pressure generated by DSJ, the ejecting enhanced by DSJ, and the co-flow effect produced by accelerated secondary jet all seem to play roles in the deflection of the primary jet.

Since the primary jet is strong enough to drive a small aircraft, the potential of dual synthetic jets to provide FTV control is initially revealed. The demonstrated control capability of DSJ makes it possible to provide thrust vector control without any moving parts or air supply. The detailed structure of flow field and the mechanism of this method need to be further studied.

In the future, relative experimental observations will be conducted to further investigate this mechanism, and flight tests will be implemented to further verify its effectiveness and control characteristics.

**Author Contributions:** Writing: original draft, J.-F.L.; writing: review and editing, Z.-B.L., X.D., Z.-J.Z., S.-Q.L., Q.L. and Y.-X.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China, grant numbers U2141252, 11972369, 11872374, the Excellent Youth Foundation of Hu'nan Scientific Committee, grant number 2020JJ2031 and the Cultivation Program for Young Talents in College of Aerospace Science and Engineering, grant numbers ZD-22-02.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Dual Synthetic Jets Actuator and Its Applications—Part III: Impingement Flow Field and Cooling Characteristics of Vectoring Dual Synthetic Jets**

**Xiong Deng 1,2,\*, Zhaofeng Dong 1, Qiang Liu 1, Can Peng 1, Wei He <sup>1</sup> and Zhenbing Luo 1,\***


**Abstract:** In order to understand the impingement flow field and cooling characteristics of vectoring dual synthetic jets (DSJ), an experimental investigation was performed to analyze the parameter effects. With the variation of the slot location, the vectoring angle of DSJ can be adjusted from 34.5◦ toward the left to 29.5◦ toward the right. The vectoring function can greatly extend the length of impingement region. There are three local peaks both for the local cooling performance (*Nu*) and the whole cooling performance (*Nuavg*). Although the peak *Nu* at a certain location of the slider is higher than that at the center, the corresponding *Nuavg* is lower. As for different driving frequencies, the vectoring angle reaches its minimum of 9.7◦ at 350 Hz, but the *Nu* is obviously improved. There is one local peak of *Nuavg* values at 350 Hz rather than three local peaks at 250 Hz and 450 Hz. The slot locations where the *Nuavg* of 250 Hz and 450 Hz reach maximum are different. With the increase in driving voltage from ±100 V to ±200 V, the vectoring angle drops from 46.9◦ to 22.2◦, but both *Nu* and *Nuavg* are improved. The maximum *Nuavg* of each driving voltage occurs at the center location of the slider. The choking effect and the cross flow have dominated the vectoring angle and the cooling performance of impingement DSJ. Vectoring DSJ will give impetus to the thermal management of large-area electric devices in spaced-constrained cooling and removing dynamic hotspots.

**Keywords:** dual synthetic jets; impingement flow field; vectoring angle; impingement cooling; thermal management

#### **1. Introduction**

Due to the unique capability of additional momentum injecting into ambient fluid without complex plumbing, synthetic jet technology shows a significant potential for active flow control. Synthetic jets have been engaged in applications of aerodynamic force control [1–5], flow separation control [6–9], mixing enhancement [10–12], jet vectoring control [13–15], heat transfer enhancement [16–19], etc. With the development of integrated electrical systems in miniaturization and high power, the cooling problem becomes a greater challenge. Traditional cooling techniques such as fans and heat pipes cannot satisfy the demands of compact cooling systems. Thus, there is an urgent need to seek a simple and high-efficiency cooling method.

Impinging jets have been widely used in electronic cooling due to the high heat transfer coefficient [20]. However, there are some obvious disadvantages for impinging jets, such as high energy consumptions and needs of external fluid plumbing. Thus, the development of an impinging synthetic jet which needs no fluid supply system has drawn researchers' attention. Because of the interactions of coherent structures with thermal boundaries [21,22], the impingement cooling performance of a synthetic jet is much better than a continuous jet [23]. Luis and Alfonso [24] revealed the mechanism of the heat transfer enhancement of an impinging synthetic jet by analyzing the vortex dynamics. They

**Citation:** Deng, X.; Dong, Z.; Liu, Q.; Peng, C.; He, W.; Luo, Z. Dual Synthetic Jets Actuator and Its Applications—Part III: Impingement Flow Field and Cooling Characteristics of Vectoring Dual Synthetic Jets. *Actuators* **2022**, *11*, 376. https://doi.org/10.3390/ act11120376

Academic Editor: Luigi de Luca

Received: 7 November 2022 Accepted: 10 December 2022 Published: 15 December 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

discovered that a secondary vortex with an opposite circulation was generated after the vortex pair arrived at the heated wall. Greco et al. [25] investigated the combined effect of impingement distance and stroke length on the cooling performance of an impinging synthetic jet. Both inner and outer ring-shaped regions occurred at a short impingement distance under the highest dimensionless stroke length. As the impingement spacing increased, the outer ring-shaped regions disappeared owing to the weakening of secondary vortex rings. Qiu et al. [26] analyzed the effect of the interactions between a circular synthetic jet and a cross flow in a microchannel on the heat transfer process. They pointed out that the heat transfer enhancement could be divided into two parts of an impinging region and an entraining region. Lyu et al. [27] investigated the heat transfer characteristics of a single synthetic jet with two planar-lobed orifices impinging on flat and concave surfaces. It was demonstrated that theN=6 petal-shaped orifice was a most promising orifice configuration. For the baseline round orifice, the spatially-averaged *Nu* was reduced up to 25% on the concave surface with respect to the flat surface. Lau et al. [28] analyzed the effects of nanoparticle types and volume concentrations on the heat transfer in a threedimensional microchannel with a synthetic jet based on a Eulerian approach. With the increase in the volume concentration, the viscosity became higher, resulting in a poorer cooling performance. The overall cooling performance of Al2O3-water nanofluid was the best at the volume concentration of 5%. Wang et al. [29] investigated the flow characteristics and unsteady heat transfer of noncircular synthetic jets impinging on a heated plate. Compared with the circular orifice, the jet penetration from the square orifice into the wall shear layer was deeper, and the maximum stagnation cooling coefficient increased 42%. The axis-switching phenomenon of elliptic and rectangular orifices enhanced the near-wall mixing and the turbulent kinetic energy. The heat transfer was improved compared with the circular orifice. The flow and thermal behavior of impingement synthetic jet driven by sinusoidal and square signals were analyzed by Singh et al. [30]. It was found that the cooling performance in the stagnation region with square signal is 18.94% higher than that with sinusoidal signal, but the averaged Nusselt number with the sinusoidal signal is 20.86% higher.

As above mentioned, the impingement cooling based on a traditional synthetic jet cannot vary the jet direction, and thus the cooling area is certain and limited. It is not applicable to the large-area cooling and the removal of dynamic hotspots. Considering this, vectoring dual synthetic jets (DSJ) was developed realized by controlling the movement of a slider based on a synthetic jet [31–34]. It has been found that DSJ has exhibited an excellent performance in impingement cooling and spray cooling [35–38]. Moreover, vectoring DSJ has been proven to be a useful cooling technology in space-constrained and large-area electronic cooling [39].

The vectoring function and the doubled energy utilization efficiency of vectoring DSJ make it potentially useful for cooling applications. The aim of the present work is to investigate the effects of slider location, driving voltage and driving frequency on the impingement flow field and cooling characteristics of vectoring DSJ. A compact thermal management strategy is provided for space-restricted electronic cooling based on a vectoring DSJ.

#### **2. Experimental Device and Methods**

#### *2.1. Experimental Device*

The experiment set up, as sketched in Figure 1, included a vectoring DSJ actuator, an impinged plate, a particle image velocimetry (PIV) system, and an infrared camera. The vectoring DSJ actuator driven by a sine signal was installed under the plate with a constant distance of 40 mm. The configurations of two slots are shown in Figure 1a, and each width of two slots could be adjusted by changing the slider location through a step motor (seeing Figure 1b). A stainless-steel foil (196 mm × 100 mm × 0.08 mm) was selected as the impinged plate.

**Figure 1.** Experimental device.

A 2D PIV system of MicroVec Corp. was used to measure the impingement flow field of vectoring DSJ. A sheet light with a thickness of 1 mm was generated from a twin-cavity Nd: YAG laser with a separation of 12 μs. The output energy of each laser pulse was about 200 mJ with a duration of 7 ns. The impingement model of vectoring DSJ was covered in a Plexiglas® container. The *xoz* plane was brightened by a laser sheet of 1 mm thickness at a frequency of 10 Hz. The particle images with a resolution of 0.035 mm/pixel were recorded by a charge-coupled device (CCD) through an optical filter. The cross-frame time was set according to the DSJ velocity. Phase-locked and phase-shifting technologies were used to capture the images at 16 equal-divided phases of the DSJ period. A velocity field has been processed using the MicroVec 3.0 software package based on a two-frame crosscorrelation algorithm with an interrogation window of 32 × 16 pixels and 50% overlap. A phase-averaged flow field was calculated according to 24 pairs of velocity fields at each phase point. Then the period-averaged flow field of impingement DSJ was achieved by averaging the 16 phase-averaged flow fields.

A FLIR infrared camera was installed above to measure the temperature distribution on the top surface of the impinged plate. The measurement range was from 5 ◦C to 2500 ◦C with an accuracy of ±1 ◦C. The impinged plate was steadily and uniformly heated by a direct-current source. Since the maximum Biot number was 0.0013, which was much less than the critical value of 0.1 along the thickness direction, the temperature of the top surface and impinged surface at the same location could be deemed to be consistent. To improve the radiation feature and the measurement accuracy, a thin and uniform layer of a black matted paint with an emissivity of 0.95 was painted on both sides of the plate. When the plate reached its thermal equilibrium, fifty thermal images were recorded with a sampling frequency of 20 Hz. The steady-status temperature distribution was obtained by averaging the fifty thermal images.

#### *2.2. Vectoring Evaluation Method of Impinging DSJ*

In order to explicitly quantify the vectoring angle of DSJ, an evaluation method was proposed based on the potential core. Firstly, a potential core was selected to cover the vast bulk of jet momentum in an evaluated flow field. Secondly, we considered a line and projected the velocity vectors in the potential core onto the line. An objective function was defined as the sum of the projection values of the velocity vectors in the potential core. Thus, the objective function could be deemed as the impinging strength of the jet momentum in the potential core along the line. Finally, the maximum of the objective function could be achieved for a specific projected line. The inclination angle of the specific line was prescribed as the vectoring angle of flow field. The detailed algorithm is explained as follows.

In this paper, the critical velocity *Vc* is defined as half of the maximum velocity magnitude in the period-averaged flow field of the vectoring DSJ. The region with velocity magnitude beyond the critical velocity is regarded as the potential core.

$$V\_{\mathfrak{c}} = \frac{V\_{\text{max}}}{2} = \frac{\max(|V\_{\mathfrak{i}}|)}{2} = \frac{\max(\sqrt{u\_{\mathfrak{i}}^2 + w\_{\mathfrak{i}}^2})}{2} \tag{1}$$

where *V<sup>i</sup>* denotes the velocity vectors with components of *ui* and *wi* along *x*-axis and *z*-axis, respectively.

According to the definition above, the objective function *So* is expressed as

$$S\_{\mathcal{O}} = \sum\_{i=1}^{n} (\mathbf{V}\_{i} \cdot \mathbf{e}\_{\mathcal{O}}) \Bigg|\_{||\mathbf{V}\_{i}|| \ge V\_{\mathcal{C}}} \tag{2}$$

where *e<sup>θ</sup>* is the unit vector along an assuming line with an inclination angle *θ* with respect to *z*-axis, depicted as

$$\varepsilon\_{\theta} = (\frac{\tan \theta}{\sqrt{1 + \tan^2 \theta}}, \frac{1}{\sqrt{1 + \tan^2 \theta}}) \tag{3}$$

Then, the objective function can be written as

$$\begin{aligned} S\_o &= \sum\_{i=1}^n \left( V\_i \cdot \mathbf{e}\_\theta \right) \Big|\_{\substack{|V\_i| \ge V\_c \\ |V\_i| \ge V\_c}} \\ &= \left( \sum\_{i=1}^n \frac{u\_i \tan \theta}{\sqrt{1 + \tan^2 \theta}} + \sum\_{i=1}^n \frac{w\_i}{\sqrt{1 + \tan^2 \theta}} \right) \Big|\_{\substack{|V\_i| \ge V\_c \\ |V\_i| \ge V\_c}} \\ &= \frac{\tan \theta \sum\_{i=1}^n u\_i + \sum\_{i=1}^n w\_i}{\sqrt{1 + \tan^2 \theta}} \Bigg|\_{\substack{|V\_i| \ge V\_c \\ |V\_i| \ge V\_c}} \end{aligned} \tag{4}$$

It can be concluded that there must be a specific angle *θ*0, at which the objective function reaches its maximum. The mathematical relation is given by

$$(S\_{\mathcal{O}})\_{\text{max}} = \lim\_{\vartheta \to \vartheta\_0} S\_{\mathcal{O}} \tag{5}$$

As mentioned above, *θ*<sup>0</sup> represents the direction of the jet momentum in the potential core, and also the vectoring angle of flow field.

#### *2.3. Cooling Performance Evaluation Method of Impinging DSJ*

The local cooling performance of impingement DSJ, namely the Nusselt number (*Nu*), is evaluated by

$$Nu = \frac{hd\_c}{\lambda} \tag{6}$$

where *λ* is the thermal conductivity of air at the average temperature of the impinged surface. The slider width of 5 mm is chosen as the characteristic length *de*. *h* is the local heat transfer coefficient of impingement DSJ expressed as

$$h = \frac{q\_{nct}}{T\_s - T\_{\dot{j}}} \tag{7}$$

where *qnet* is the net heat flux removed by impingement jet, *Ts* and *Tj* are the local temperature of impingement surface and jet temperature, respectively.

The net heat flux is estimated as follows

$$q\_{net} = \frac{Q\_{elv} - Q\_{loss}}{A} \tag{8}$$

where *Qele* is the heating power, *Qloss* is the total heat loss, *A* is the effective heating area of impinged surface.

The area-averaged Nusselt number (*Nuavg*) of the whole impingement surface is calculated by:

$$Nu\_{avg} = \frac{\int\_{\Omega} Nudxdy}{A\_{\Omega}} \tag{9}$$

where *A*Ω is the image area recorded by the infrared camera.

In the experiments, the heating power is kept a constant of 40 W with a measurement error less than 1.4%. The total loss heat resulted from the radiation and the lateral is estimated to be within 6.4% of the heating power. The uncertainty of *Nu* is less than 8% through the evaluation method [40].

#### **3. Results and Discussion**

#### *3.1. Effect of Slider Location*

The slider location determines the area ratio of two slots which influences the momentum ratio of two synthetic jets and the strength of the low-pressure region. Thus, different momentums of two synthetic jets result in the deflection of DSJ. Comparing to a traditional synthetic jet, vectoring DSJ can greatly extend the impingement area. The PIV results with different dimensionless slider location *d*\*, defined as *d*/*de*, are exhibited in Figure 2. The driving voltage and driving frequency are set to ±150 V and 450 Hz, respectively.

**Figure 2.** Velocity vector maps, contours, and vectoring angles at different slider locations.

As demonstrated in Figure 2, the DSJ deflects toward the side of the slot with a greater cross section. The reason is that the jet momentum from a smaller width slot is higher, which can be deemed as a primary jet. Conversely, the other jet with a lower momentum is

regarded as a disturbing jet. As the primary jet issues from the slot, the other jet is in the suction phase. A low-pressure region is generated near the slot, which causes the primary jet deflects. As the disturbing jet is formed in the blowing phase, the primary jet has moved downstream. The entrainment of vortex structures in the disturbing jet accelerates the deflection of the downstream primary jet. So, the DSJ deflects toward the slot with a larger cross section. With *d*\* varying from −0.2 to 0.2, the vectoring DSJ can deflect from 34.5◦ toward the left to 29.5◦ toward the right. Since the momentums of two synthetic jets are equal at *d*\* = 0, there is no deflection. The vectoring angle of 1.3◦ at *d*\* = 0 (in Figure 2c) may be attributed to the asymmetry of two slots resulting from machining errors. It is observed that the stagnation region length is about 8*de* (from −4*de* to 4*de*) for the normal DSJ (seeing Figure 2c), while for the vectoring DSJ it increases to 15*de* (from more than −8*de* in Figure 2a to 7*de* in Figure 2e). Additionally, the jet velocity becomes a little higher at an off-center slider location. It has been proved that the heat transfer coefficient decreases with the increase in the inclination angle of synthetic jet, and reaches the maximum for the normal impingement [41]. Considering the vectoring characteristic, the vectoring DSJ can be obliquely installed to make its impingement normal. This feature will be greatly useful in spaced-constrained cooling applications.

Figure 3 shows the *Nu* distributions on the impinged surface at different slider locations. The driving voltage and the driving frequency are kept at ±150 V and 500 Hz, respectively. It is indicated that the *Nu* distribution at *d*\* = 0 is self-symmetrical and two *Nu* distributions at a same |*d*\*| are symmetrical about the line *x*/*de* = 0. The high *Nu* area gradually moves away from the center with the increase in |*d*\*|. The reason is that with the increase in |*d*\*|, the vectoring angle of DSJ becomes larger which causes the stagnation region moving away from the center.

**Figure 3.** *Nu* contours on the impingement surface at different slider locations.

Figure 4a draws the *Nu* on the centerline *y*/*de* = 0 along the *x*-axis. The solid dot represents the peak *Nu* on the curves. It can be seen that the curves at the same |*d*\*| are symmetrical with the line *x*/*de* = 0. With the increase in |*d*\*|, the curves and peak point of *Nu* gradually move away from the center, which is consistent with the variation trend of the vectoring angle (shown in Figure 2). The peak *Nu* reaches its maximum at |*d*\*| = 0.08 rather than at the center location. The reason is that the DSJ velocity increases with the increase in |*d*\*|, which improves the local cooling performance. As |*d*\*| continues to increase, the thickness of the viscous layer in the slot neck gradually grows and eventually fills the entire neck. As a result, the choking effect occurs in the slot neck and the DSJ velocity reduces. Moreover, the vectoring angle increases with the increase in |*d*\*| (shown in Figure 2) and the impingement DSJ will gradually transform to a cross flow, resulting in the weakness of the cooling performance. There are three local peaks of *Nuavg* at *d*\* = 0 and ±0.28 (shown in Figure 4b). The reason is that at the beginning of |*d*\*| increase, the DSJ velocity increases, which improves the whole cooling performance. As |*d*\*| continues to increase, the choking effect impairs the DSJ velocity, and meanwhile the increase in vectoring angle enhances the cross-flow effect of the DSJ. Thus, the whole cooling performance reduces rapidly. Although the maximum values of *Nu* at *d*\* = ±0.08 and ±0.16 are larger than that at *d*\*=0 (shown in Figure 4a), the values of *Nuavg* are smaller (shown in Figure 4b). Since the DSJ velocity *d*\* = ±0.08 and ±0.16 are higher than that at *d*\* = 0 owing to the smaller cross sections, and the local cooling performance. However, the vectoring angles at *d*\* = ±0.08 and ±0.16 enhance the cross-flow effect and reduce the impingement stagnation region of the DSJ, which causes the decrease in overall cooling performance.

#### *3.2. Effect of Driving Frequency*

Previous work [42] has shown that the frequency response of a piezoelectric-driven synthetic jet actuator is consistent with that of a damped fourth order system governed by two key resonance frequencies—the natural frequency of the vibrating diaphragm and the Helmholtz resonance frequency. There are two local peaks of jet velocity and energy efficiency at the two resonance frequencies. In this section, experiments on the effect of driving frequency are carried out at a constant driving voltage of ±150 V and *d*\* = 0.2. The PIV results are illustrated in Figure 5.

It is observed from Figure 5 that the vectoring angle shows a non-monotonic variation trend with the increase in driving frequency. The minimum vectoring angle is 9.7◦ at the driving frequency of 350 Hz. The natural frequency of vibrating diaphragm is about 360 Hz. Since the amplitude of the vibrating diaphragm will increases significantly near the natural frequency, the DSJ velocity at 350 Hz is much higher than that at the other two frequencies. A higher DSJ velocity requires a greater force to deflect. Thus, the vectoring angle and the impingement stagnation region become small. To prevent electrical and mechanical failure of the vibrating diaphragm, the driving frequency should not be set at the resonant

frequency. The driving frequency should be chosen by balancing the cooling performance and the service life of the DSJ actuator.

**Figure 5.** Velocity vector maps, contours and vectoring angles at different driving frequencies.

*Nu* contours on the impinged surface at different driving frequencies and slider locations are displayed in Figure 6. It can be seen that the peak *Nu* and the high *Nu* area of 350 Hz are obviously larger than those of 250 Hz and 450 Hz for each slider location. It is attributed to the frequency of 350 Hz close to the resonant frequency where the jet velocity is very higher (shown in Figure 5). As a result, the local heat transfer coefficients are higher than those of other two frequencies which are far away from the resonant frequency. For the driving frequency of 350 Hz, the peak *Nu* at *d*\* = 0 is greater than those at *d*\* = ±0.24. The tendencies are contrary at driving frequency of 250 Hz and 450 Hz. The reason is that the DSJ velocity of 350 Hz is high enough at *d*\* = 0, and the choking effect occurs in the narrow slot neck during the slider moving away from the center. It means although the cross section of the slot becomes smaller, the jet velocity may decrease, such as the slider locations of *d*\* = ±0.24. For the driving frequency of 250 Hz and 450 Hz, the DSJ velocities are low at *d*\* = 0 and increase with the decrease in the slot cross section at *d*\* = ±0.24. A higher DSJ velocity produces a better local cooling performance and a higher local *Nu*.

**Figure 6.** *Nu* contours at different driving frequencies and slider locations.

*Nuavg* with different drive frequencies and slider locations are shown in Figure 7. Results in Figure 7 indicate that there are three local peaks at *d*\* = 0 and *d*\* = ±0.24 both for 250 Hz and 450 Hz, which is similar to Figure 4b. The distinction is that the maximum *Nuavg* of 450 Hz occurs at *d*\* = 0 rather than that of 250 Hz at *d*\* = ±0.24. The reason is that at

*d*\* = 0 the DSJ velocity of 250 Hz is smaller than that of 450 Hz and increases more obviously at *d*\* = ±0.24. In addition, the impingement area is larger at *d*\* = ±0.24 than that at *d*\*=0. As a result, the cooling performance at *d*\* = ±0.24 is better than that at *d*\* = 0. There is only one local peak with the drive frequency of 350 Hz, which is different from the drive frequency of 250 Hz and 450 Hz. It is attributed to that the DSJ velocity at 350 Hz becomes small owing to the choking effect in the slot neck when the slider moves away from the center. The decrease in DSJ velocity is more remarkable than the increase in impingement stagnation region, and thus the whole cooling performance gradually becomes weak. For the driving frequency of 250 Hz, although the DSJ velocity at *d*\* = ±0.08 is a little higher than that at *d*\* = 0, the cross-flow effect is more obvious due to the vectoring deflection. Overall, the whole cooling performance decreases a bit. As the slider location moves to *d*\* = ±0.24, the increase in the DSJ velocity is more remarkable than the enhancement of the cross-flow effect. Thus, the whole cooling performance gradually improves. As |*d*\*| continues to increase, the choking effect occurs and impairs the DSJ velocity. The whole cooling performance reduces. The explanation for the variation trend of *Nuavg* at 350 Hz is similar to the above mentioned about 250 Hz. The differences are that the maximum *Nuavg* of 450 Hz occurs at *d*\* = 0 rather than that of 250 Hz at *d*\* = ±0.24. The reason is that the DSJ velocity of 250 Hz can be obviously enhanced with the increase in |*d*\*|, which improves the whole cooling performance.

**Figure 7.** *Nuavg* variations with the driving frequency and the slider location.

#### *3.3. Effect of Driving Voltage*

The driving voltage influences the input power imposed on the vibrating diaphragm, and then regulates the velocities of the primary jet and the disturbing jet. To investigate the effect of driving voltage experiments under three different driving voltages are executed. The corresponding velocity vector maps and contours at the driving frequency of 450 Hz and *d*\* = 0.2 are shown in Figure 8.

It can be seen that the velocity becomes high with the increase in the driving voltage, but the vectoring angle of DSJ decreases from 46.9◦ to 22.2◦. It is attributed to that a high driving voltage can improve the amplitude of vibrating diaphragm, as well as the jet velocity. As the driving voltage increases, both of the primary jet and the disturbing jet are enhanced. Thus, the primary jet requires more force to deflect. However, the deflection force generated by the suction phase of the disturbing jet increases a little. The reason is that the fluid away from the primary jet is easily inhaled due to the three-dimensional effect of the suction process. Thus, the vectoring angle of DSJ becomes small. It is also observed that in view of the deflection symmetry, the impingement region of vectoring DSJ

can be obviously extended. This feature is useful for improving the impingement cooling uniformity and removing dynamic hotspots.

**Figure 8.** Velocity vector maps, contours and vectoring angles at different driving voltages.

Figure 9 shows *Nu* distributions on impinged surface at different driving voltages and slider locations. It is indicated that for each slider location, the peak *Nu* and the high *Nu* region gradually increase with the increase in the driving voltage. The reason is that the primary jet velocity becomes high with the increase in the driving voltage, and thus the local heat transfer coefficient increases. The core of the high *Nu* region can move toward the center location, which is agreement with the variation trend of the vectoring angle with the increase in the driving voltage.

**Figure 9.** *Nu* contours at different driving voltages and slider locations.

*Nuavg* at different driving voltages and slider locations are plotted in Figure 10. The curves indicate that *Nuavg* improves with the increase in the driving voltage at each slider location. It is attributed to the obvious increases in the peak *Nu* and the high *Nu* region (shown in Figure 9). For each driving voltage, there is one peak of *Nuavg* at the center location of the slider and decreases when the slider moves away from the center. The reason is that the primary jet velocity is so high at the center location of the slider that the choking effect occurs when the slider moves away from the center. As a result, the DSJ velocity reduces, as well as the heat transfer coefficient. Moreover, the impingement DSJ can transform to a cross flow owing to the vectoring deflection. Above two reasons make the whole cooling performance weak. Although a high driving can improve the cooling performance, the energy consumption also increases. The driving voltage should be selected by considering the energy utilization efficiency in the cooling applications.

**Figure 10.** *Nuavg* variations with driving voltage and slider location.

#### **4. Conclusions**

The effects of slot location, driving voltage, and driving frequency on impingement flow field and cooling characteristics of vectoring DSJ are experimentally investigated. Conclusions are summarized as follows:

(1) The DSJ deflects toward the side of the slot with a greater cross section. The vectoring angle can be adjusted from 34.5◦ toward left to 29.5◦ toward right. The length of the impingement region is over two times as large as that of the normal DSJ. With the increase in |*d*\*|, the maximum *Nu* increases first and then decreases. There are three local peaks of *Nuavg* at *d*\* = 0 and ±0.28. Although maximum *Nu* at *d*\* = ±0.08 and ±0.16 are larger than that at *d*\* = 0, the *Nuavg* are smaller.

(2) The vectoring angle is lowest of 9.7◦ at the driving frequency of 350 Hz, but the peak *Nu* and high *Nu* area are obviously augmented. There is one local peak of *Nuavg* at 350Hz rather than three local peaks of *Nuavg* at 250 Hz and 450 Hz. The difference is that the maximum *Nuavg* of 250 Hz and 450 Hz, respectively, occurs at *d*\* = 0 and at *d*\* = ±0.24.

(3) With the increase in the driving voltage, the vectoring angle decreases from 46.9◦ to 22.2◦, but the length of the impingement area is extended. There is one local peak of *Nuavg* at *d*\* = 0 for each driving voltage.

**Author Contributions:** Conceptualization, X.D. and Z.L.; methodology, X.D. and W.H.; formal analysis, X.D. and Z.D.; investigation, X.D. and Q.L.; writing—original draft, X.D.; writing—review and editing, C.P. and Z.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China (grant numbers 11972369, 11872374) and the Science and Technology Innovation Program of Hunan Province (2021RC3075).

**Data Availability Statement:** Data available on request from the authors.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclature**


#### **References**


## *Article* **Dual Synthetic Jet Actuator and Its Applications—Part IV: Analysis of Heat Dissipation and Entropy Generation of Liquid Cooling with Dual Synthetic Jet Actuator**

**Ying Kang, Zhenbing Luo \*, Xiong Deng, Yinxin Zhu and Zhixun Xia**

College of Aerospace Science and Engineering, National University of Defense Technology, Changsha 410073, China **\*** Correspondence: luozhenbing@163.com

**Abstract:** Increasing heat flux restricts the development of the miniaturization of electronic devices. There is an urgent need for a heat dissipation method that will efficiently cool the chip. This paper presents a novel liquid cooling device based on dual synthetic jets actuator (DSJA) technology. The characteristics of the temperature and velocity field of the device are numerically studied by a three-dimensional coupled heat transfer model. The entropy generation rate caused by heat transfer and fluid friction was studied to analyze the effective work loss and irreversibility of the heat transfer process. When the DSJA is turned on, the temperature of the heat source with a heat flux of 200 W/cm2 is 73.07 ◦C, and the maximum velocity is 24.32 m/s. Compared with the condition when the the DSJA is closed, the temperature decreases by 25.15 ◦C, and the velocity increases by nearly 20 m/s. At this time, the total inlet flow is 1.26 L/min. The larger frictional entropy generation is mainly distributed near the inlet and outlet of the channel and the jet orifice. The higher the velocity is, the more obvious the frictional entropy generation is. Due to the large temperature gradient, there is a large thermal entropy generation rate at the fluid–solid interface.

**Keywords:** liquid cooling; dual synthetic jets actuator; entropy generation

### **1. Introduction**

The continuous development of micro-electronic devices will depend on efficient heat dissipation technology. At present, it is a great challenge to design a heat sink with a heat flux exceeding 1 MW/m2 [1]. Many electronic devices are damaged by overheating. When the temperature of electronic equipment increases by 10 ◦C, the reliability decreases by 50% [2]. Therefore, it is urgent that we seek an efficient heat dissipation technology.

Heat dissipation technology is divided into natural convection technology and forced convection technology. Forced convection technology can rapidly reduce the surface temperature of electronic devices, including spray cooling technology, microchannel cooling technology, synthetic jet technology, etc. The spray cooling device has high heat dissipation efficiency but the system is complex. The microchannel technology has the advantage of small volume and high heat transfer. Many researchers have studied the structural parameters of microchannels [3–13]. Rhombus fractal-like units [5] and spider-netted microchannels [7] have excellent heat dissipation performance. Rectangular microchannels are widely used in our lives due to their good heat dissipation and hydraulic performance. Many studies have shown that microchannels have a strong heat dissipation capability. However, to meet the heat dissipation demand, they often need to provide more pump power [14].

Synthetic jet technology [15] is an active flow control technology, which can be used in the field of heat transfer enhancement. The synthetic jet actuator (SJA) comprises a diaphragm, a cavity, and a jet orifice. The diaphragm vibrates periodically, changing the volume of the cavity. When the cavity volume is compressed, the fluid in the cavity is ejected through the jet orifice, forming a jet at the outlet of the actuator, and developing

**Citation:** Kang, Y.; Luo, Z.; Deng, X.; Zhu, Y.; Xia, Z. Dual Synthetic Jet Actuator and Its Applications—Part IV: Analysis of Heat Dissipation and Entropy Generation of Liquid Cooling with Dual Synthetic Jet Actuator. *Actuators* **2022**, *11*, 382. https://doi.org/10.3390/ act11120382

Academic Editor: Luigi de Luca

Received: 8 November 2022 Accepted: 16 December 2022 Published: 19 December 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

continuously downstream. When the cavity volume is expanded, the external fluid is sucked into the cavity through both sides of the jet orifice. So far, the SJ has completed a periodic movement [15,16]. Different parameters have different effects on the heat transfer of the SJ. There is an optimal frequency to realize the maximum Nusselt number of SJs impinging on a constant heat flux disk. In addition, when the jet-to-surface distance is small, the Nusselt number increases with the jet-to-surface distance [17]. SJs can significantly destroy the thermal boundary layer to enhance heat transfer [18,19]. However, one side of the diaphragm of the SJA is located at the air side, which can easily cause ballast failure and low energy utilization.

Based on SJ technology, Luo et al. [20] invented the dual synthetic jet (DSJ) technology. The DSJ has superior performance and overcomes the problems of SJs [21]. The DSJ technology consists of two cavities, one diaphragm, and two jet outlets. The diaphragm separates the cavity. When the diaphragm vibrates, the cavity on one side expands and the cavity on the other side compresses, the jet outlet on the expanding side inhales the surrounding fluid, and the compressed side ejects the fluid in the cavity. The alternately ejected/inhaled jets interact near the outlet and merge into the DSJ, which continues to develop downstream. Research on SJ and SJA technology mostly focuses on the field of air cooling. The SJA technology has been proven to be well applied to heat dissipation in confined spaces [22].

For the heat dissipation of higher heat flux, the ability of air cooling is not enough. Therefore, some scholars have preliminarily studied the heat transfer mechanism of SJs in the liquid cooling field. Combining SJs with microchannels, the influence of the interaction between the jet and incoming flow on heat transfer is studied [23–28]. The PIV technique was used to observe the flow field of SJs underwater. Compared with the condition without SJs, the effect of an SJ makes the heat transfer increase by about 4.3 times [23]. The SJ promotes fluid mixing in the channel and enhances heat transfer [29]. It is worth noting that the SJ contributes little to the pressure drop in the channel. With a low flow rate, heat transfer can be enhanced by 130% [30]. Multiple synthetic jets can be combined with microchannels at the same time for better heat dissipation [26]. The convective heat transfer capacity of two SJs is more effective when 180◦ is out of phase [31,32]. The higher the frequency and amplitude of the diaphragm, the higher the convective heat transfer capacity. There is an optimal fluid inlet temperature of 297.15 K for achieving maximum heat transfer. The overall performance of the heat sink increases significantly with the Reynolds number at the inlet [32]. The underwater PIV experiment of DSJ technology shows that the jet strength is the highest when the driving frequency is 30 Hz [33]. The heat transfer characteristics and physical properties of the SJ flow field in the channel are very complex, so a comprehensive parameter study of SJAs is required [34].

Some other scholars have carried out entropy generation analysis on the studied device. Entropy generation analysis is helpful to understand the source of irreversible loss. Luis et al. [35] analyzed local entropy generation and global entropy generation of different parts of a hybrid microjet heat sink. It is found that the irreversibility of the system mainly comes from the heat transfer of the copper plate. The total entropy generation consists of thermal entropy generation and friction entropy generation. Mehdi et al. [36] found that the latter has little influence on the irreversible loss of the heat sink. Omid et al. [37] reviewed the entropy generation theory of nanofluids and introduced two calculation methods for entropy generation.

In our previous research, a side-mounted liquid cooling device was proposed [38]. Good heat transfer and flow performance are obtained. To further simplify the fluid circuit and reduce the flow resistance, this paper proposes a front-placement liquid cooling device. In this paper, the DSJ is combined with a channel for liquid cooling. The entropy generation of the cooling device is analyzed from the point of view of the second law of thermodynamics. The novelty of this study is that a new type of forwarding compound DSJ liquid cooling device is proposed. It is of great significance to solve the heat dissipation of electronic devices in a high heat flux confined space. At the same time, this study fills in

the blanks that remain from the study of DSJs by the entropy generation analysis method. The research contents of this paper are as follows: Section 2 mainly introduces the relevant models and numerical methods; Section 3 is the result discussion. The temperature field, velocity field, entropy generation, and the effect of diaphragm frequency have been studied. Section 4 is the conclusion.

#### **2. Numerical Methods**

#### *2.1. Physical Model*

In this paper, the DSJA cooling device consists of a DSJA, a channel, and an aluminum plate, as shown in Figures 1 and 2. The jet orifices are directly arranged on the upper side of the cavity to connect the channel and the actuator cavity. The channel consists of an inlet and an outlet. Jet orifices set to 5 × 5 arrays. The fluid enters the cooling device in three ways. One fluid enters the channel, and the other two fluids enter the channel through the jet orifices, through the actuator. The three converge in the channel and move downstream. The yellow arrow in Figure 1 indicates the path of fluid flow. The jet orifices are row 1, row 2, row 3, row 4, and row 5 along the positive direction of the *y*-axis.

**Figure 2.** Schematic diagram of channel: (**a**) Vertical view; (**b**) Side view.

When the DSJA diaphragm is not vibrating, the actuator only acts as a fluid circuit. At this time, jets from the jet orifices can be regarded as steady jets. When the DSJA diaphragm vibrates upward, the cavity on the upper side of the actuator shrinks, and fluid is ejected from the third, fourth, and fifth-row jet orifices. At the same time, the cavity at the lower side of the actuator expands, and the fluid in the channel flows into the cavity at the lower side through the 1st and 2nd discharge orifices. When the DSJA diaphragm vibrates downward, the flow situation is opposite to the above process. Different jet outlets alternately eject/inhale fluid, and the formed jets fuse and develop with each other, and interacts with the incoming flow of the channel to form a violent disturbance.

#### *2.2. Governing Equations*

In this paper, a three-dimensional fluid–solid coupling model is established. The thermophysical properties of fluids and solids are considered constant. The fluid is assumed to be incompressible. The working fluid is water, and the solid is aluminum. The mass, momentum, and energy equations of the problem [25] can be written as

$$\frac{\partial \rho}{\partial t} + \nabla \cdot \rho \mathbf{U} = 0 \tag{1}$$

$$\frac{\partial}{\partial t}(\rho lI) + \nabla \cdot (\rho llII) = -\nabla p + \nabla \cdot (\overline{\mathbf{r}}) \tag{2}$$

$$\frac{\partial}{\partial t}(\rho E) + \nabla \cdot (\mathcal{U}(\rho E + p)) = \nabla \cdot \left(\lambda\_{\rm eff} \nabla T + \left(\stackrel{\rm \rm \\_}{r\_{\rm eff}} \mathcal{U}\right)\right) \tag{3}$$

#### *2.3. Boundary Conditions and Numerical Solver Setting*

The simulation software used in this research is Fluent. Pressure-based solver is selected for calculation. The energy equation is opened. Considering the need to solve the jet impact problem and large-scale flow changes, the SST k-ω model is selected. The calculation includes a fluid domain and a solid domain. The material of the fluid domain is water, and the material of the solid domain is aluminum.

The boundary conditions of the cooling device are set as shown in Figure 3. Considering that the cooling device is completely symmetrical, only half of the model is calculated. Symmetrical boundary conditions are adopted for the model section. The entrance of the channel and the entrance of the actuator adopt the velocity inlet boundary. The inlet velocity is 1.5 m/s and 0.3 m/s, respectively. The channel outlet is set as the pressure outlet. The chip uses a constant heat flux boundary. The inlet water temperature is 25 ◦C. The fluid–solid interface is set as the coupling surface. The rest are set as non-sliding walls.

**Figure 3.** Boundary conditions of the calculation model.

The actuator diaphragm uses velocity inlet boundary conditions of user-defined functions (UDF) as

$$\begin{cases} \boldsymbol{u}(t)\_I = 0.02\sin(2\pi ft + 0^\circ) \\ \boldsymbol{u}(t)\_I = 0.02\sin(2\pi ft + 180^\circ) \end{cases} \tag{4}$$

The pressure and velocity are coupled by a SIMPLE algorithm. The relaxation factors for pressure, density, and momentum are set by default to 0.3, 1, and 0.7, respectively. The criterion of the energy equation is 10−6, and the convergence criterion of mass and momentum is 10<sup>−</sup>3. In this paper, the frequency is selected as *f* = 30 Hz [33]. The time step is set to 2.778 × <sup>10</sup>−<sup>4</sup> s, and the total time is 6.67 s.

#### *2.4. Numerical Method Validation*

Grid independence verification is performed on the model, as shown in Table 1. Because the model is completely symmetrical, only half of the model is used for structural

mesh generation. Three sets of grids are divided, and the *y+* value is less than 1. The error between the number of grids (8,555,843) and the number of grids (11,713,823) is less than 8%. To ensure the calculation accuracy and velocity, 8,555,843 grids are selected for calculation, as shown in Figure 4.

**Table 1.** Grid independence verification.


**Figure 4.** Structural grid of the model: (**a**) Vertical view; (**b**) Front view; (**c**) Side view; (**d**) Local details of (**c**).

The experimental model in reference [33] is selected to verify the numerical method in this paper. As shown in Figure 5, the RNG k-ε model and SST k-ω model are selected. The velocity distribution at 2 mm away from the DSJ orifice was studied. The results show that the SST k-ω model is in better agreement with the experiment.

#### *2.5. Relevant Definitions*

The hydraulic diameter of the channel is

$$D\_h = \frac{2w}{w+h\_0} \tag{5}$$

where *w*, and *h*<sup>0</sup> are the channel width and height, respectively. In this paper, *Dh* = 1.8182 mm.

The channel Reynolds number *Rec* [28] and the DSJ Reynolds number *Redsj* are defined as follows:

$$Re\_{\mathfrak{c}} = \frac{\rho u\_{\rm in} D\_{\rm h}}{\mu} \tag{6}$$

$$Re\_{dsj} = \frac{\rho u\_{dsj} d}{\mu} \tag{7}$$

where *uin* = 1.5 m/s, *udsj* = 24.32 m/s, and *d* = 0.36 mm are the flow velocity of the channel inlet, the maximum flow velocity at the exit of the actuator, and the diameter of the jet orifice, respectively. The value of the channel Reynolds number is 2714. The range of jet Reynolds numbers *Redsj* is 43,998.192.

The Nusselt number of a channel is defined as

$$Nu = \frac{D\_h h}{\lambda} \tag{8}$$

where *h* = *<sup>q</sup>* (*Tw*−*Tf*) is the convective heat transfer coefficient of DSJA, where *Tw* and *Tf* are the solid surface temperature and the fluid temperature, respectively.

The volume entropy production rate includes the entropy production rate caused by irreversible heat transfer and the irreversible flow loss caused by fluid friction. All formulas [36] are as follows:

$$
\dot{\boldsymbol{S}}\_{\text{gen,total}}^{\prime\prime\prime} = \dot{\boldsymbol{S}}\_{\text{gen,thermal}}^{\prime\prime\prime} + \dot{\boldsymbol{S}}\_{\text{gen,friction}}^{\prime\prime} \tag{9}
$$

$$\dot{S}\_{gen, thermal}^{\prime\prime} = \frac{\lambda\_f}{T\_f^2} \left[ \left( \frac{\partial T}{\partial \mathbf{x}} \right)^2 + \left( \frac{\partial T}{\partial y} \right)^2 + \left( \frac{\partial T}{\partial z} \right)^2 \right] \tag{10}$$

$$\begin{split} \dot{S}\_{gen,frictional}^{\prime\prime} &= \frac{\mu}{T\_f} \left\{ 2 \left[ \left( \frac{\partial \mu}{\partial x} \right)^2 + \left( \frac{\partial \underline{v}}{\partial y} \right)^2 + \left( \frac{\partial \underline{w}}{\partial z} \right)^2 \right] + \left( \frac{\partial \underline{u}}{\partial y} + \frac{\partial \underline{v}}{\partial x} \right)^2 + \left( \frac{\partial \underline{u}}{\partial z} + \frac{\partial \underline{w}}{\partial x} \right)^2 \\ &\quad + \left( \frac{\partial \underline{v}}{\partial z} + \frac{\partial \underline{w}}{\partial y} \right)^2 \right\} \end{split} \tag{11}$$

*u*, *v*, *w* are velocity components in *x*, *y*, *z* directions respectively.

The total entropy generation is the integral sum of thermal entropy generation and frictional entropy generation in the entire cooling device.

$$\int\int\int\_{\Omega} \dot{\mathbf{S}}\_{\text{gen}}dV = \int\int\int\_{\Omega} \dot{\mathbf{S}}\_{\text{gen},\text{thermal}}dV + \int\int\int\_{\Omega} \dot{\mathbf{S}}\_{\text{gen},friction}dV\tag{12}$$

**Figure 5.** Numerical simulation and experimental [33] verification of DSJ underwater (*z* = 2 mm).

#### **3. Results**

#### *3.1. Temperature Field*

To study the enhanced heat transfer mechanism of the DSJ, the central section of the channel (section *x* = 28.5 mm) is taken as an example for analysis. The inlet velocity of the channel is 1.5 m/s, and the inlet velocity of DSJA is 0.3 m/s. The heat flux of the heat source is 200 W/cm2. The heat source is located at −5.25 mm < *<sup>y</sup>* < 5.25 mm.

The temperature field of the channel is shown in Figure 6. The average temperature of the chip is 120.35 ◦C. It is observed that the substrate temperature near the outlet is higher than that of the inlet. This is because the fluid takes away heat along the flow direction, increasing the fluid temperature. A thermal boundary layer attached to the substrate is gradually thickening along the flow direction, weakening the heat transfer capability. There are partitions between adjacent channels to enhance the heat conduction effect.

**Figure 6.** Temperature distribution of channel's center section (*x* = 28.5 mm).

The temperature fields of the DSJ when off and the DSJ when on are shown in Figures 7 and 8, respectively. In Figure 7, the maximum temperature of the heat source is 98.22 ◦C. The thermal boundary layer gradually develops and thickens along the channel's top wall from the channel inlet. It is destroyed in the jet impact area. Compared with other areas, a part of the thermal boundary layer is observed above the fourth and fifth row of jet orifices, considering the jet deflection caused by the weak jet of the upper cavity itself and the incoming flow. The impact damage ability to the thermal boundary layer is weakened. The existence of the thermal boundary layer will deteriorate the heat transfer, which will help to reduce the temperature of the actuator.

**Figure 7.** Temperature distribution (DSJ off) of channel's center section (*x* = 28.5 mm).

**Figure 8.** Temperature distribution of channel's center section (*x* = 28.5 mm) (DSJ on): (**a**) 0 T (1 T); (**b**) 0.25 T; (**c**) 0. 5 T; (**d**) 0.75 T.

Compared with the DSJ when off, the maximum temperature drop of the heat source is 25.15 ◦C when the DSJ is on, as shown in Figure 8. The alternating vibration of the diaphragm of the DSJA has a periodic effect on the temperature distribution of the aluminum plate. At 0 T (1 T), the thermal boundary layer was observed above the third, fourth, and fifth-row jet orifice. Consider that the upper cavity is in the suction stage at this time. The attachment of the thermal boundary layer results in an obvious high-temperature region. The temperature of the aluminum plate above the first and second rows of jet orifices is obviously lower, which is caused by jet impingement cooling. At 0.25 T, the temperature contour shows an obvious low-temperature area, and the thermal boundary layer has been destroyed. At this time, the upper cavity is in the ejection stage, and the jet ability is strong. At 0.5 T, the temperature distribution near the impact surface is relatively uniform. Consider the weakening stage of the jet from the upper cavity. At 0.75 T, the temperature of the corresponding impact area of the first and second row of jet orifices decreased significantly. This is due to the impact of the jet at the outlet of the lower cavity. Influenced by the flow of the jets and incoming flow in the suction channel of the cavity, a partially disordered thermal boundary layer was observed. Then, the cooling device enters the stage of 0 T (1 T) again, and the above process is repeated. The periodic vibration of the diaphragm brings about the periodic change of the temperature of the aluminum plate, and the reciprocating disturbance of the jet and its interaction with the incoming flow strengthen the convection heat transfer.

The distribution of Nu on the impact surface (*z* = 44 mm) along the y direction is shown in Figure 9. When the DSJ is off, the distribution of Nu is low in the middle and high on both sides. This is because the middle area is close to the heat source and has a higher temperature. Both sides have low-temperature regions because they are far away from heat sources and close to the inlet and outlet of the channel. When the DSJ is on, the distribution of Nu varies at different times. The distribution of Nu is opposite to the temperature distribution of the DSJA in Figure 8. The area where the thermal boundary layer is damaged has a large Nu number. The importance of breaking the thermal boundary layer for enhancing heat transfer is further illustrated. In addition, at the moment of 0.75 T~0 T (1 T), a larger Nu number is observed in the area of *y* < 0. It indicates that the lower cavity is in the injection stage at this time. Since there are only two rows of jet outlets in the lower cavity, the jet velocity is higher. At 0.75 T, the impact area is wider, corresponding to the wider Nu distribution in Figure 9. In a word, the effect of the DSJ increases the Nu number and presents a "wavy" Nu distribution curve. It reflects the effective disturbance of the DSJ to the flow field.

**Figure 9.** The distribution of *Nu* on the impact surface (*z* = 44 mm) along the *y* direction.

#### *3.2. Flow Field Characteristics*

This section studies the flow field characteristics of cooling devices, as shown in Figures 10–13. It can be observed from Figure 10 that when there is only a channel, the streamlines in the channel are parallel to each other along the streamlines. The fluid in the channel does not affect each other. Thick velocity boundary layers were observed near the upper and lower walls of the channel.

When the DSJ is off, the maximum velocity in the flow field reaches 4.5 m/s. The fluid entering the channel flows parallel to the channel at the beginning, and its flow path is disturbed when it passes near the jet orifice under the influence of the jet. The fluid becomes disordered. The channel flow is "squeezed" upward by the jet, and a low-velocity area (the blue part in Figure 11) is observed between adjacent jets. At the same time, the jet deflects along the flow direction under the action of channel flow. The fourth and fifth rows of jets deflect most obviously. It is consistent with the phenomenon described in Section 3.1.

When the DSJ is on, the flow field distribution of one cycle of the DSJA is shown in Figures 12 and 13. At 0 T, the actuator diaphragm vibrates downward. The upper cavity volume expands while the lower cavity volume compresses. An obvious vortex structure was observed in the upper cavity. Since the vibration of the diaphragm is small at this time, the flow from the actuator will inhibit the suction of fluid in the channel. In addition, when the jet interacts with the incoming flow, the flow in the channel is disordered, while when the jet is far away, the flow in the channel remains stable. The region of a steady fluid flow corresponds to the gradually thickened thermal boundary layer in Figure 8a. It further shows that the enhancement of fluid kinetic energy destroys the thermal boundary layer and enhances heat transfer. At 0.25 T, the upper cavity ejects fluid, and the jet impinges on the aluminum plate surface to form an obvious low-temperature area. At 0.5 T, the diaphragm is in a balanced position. The ability of the upper cavity to eject the jet is weakened. The jet deflects under the influence of incoming flow. At 0.75 T, the lower cavity ejects a strong jet. An obvious vortex structure is formed in the channel, and the temperature of the jet impingement zone is obviously reduced. At the same time, it was observed that the fluid in some channels was sucked into the upper cavity. The large

vortex structure was observed above the third, fourth, and fifth jet orifices, because it was affected by the three actions of ejected fluid, inhaled fluid, and incoming flow. The resulting temperature distribution is consistent with that observed in Figure 8.

**Figure 10.** Velocity distribution of channel's center section (*x* = 28.5 mm).

**Figure 11.** Velocity distribution (DSJ off) of channel's center section (*x* = 28.5 mm): (**a**) overall view; (**b**) partial enlarged view of jet orifice area.

**Figure 12.** Velocity distribution of channel's center section (*x* = 28.5 mm, overall view) (DSJ on): (**a**) 0 T (1 T); (**b**) 0.25 T; (**c**) 0.5 T; (**d**) 0.75 T.

To further analyze the flow field characteristics of the DSJA, the vorticity distribution of the DSJA is shown in Figure 14. When the DSJ is off, only large vortices are observed in the channel, especially near the jet orifices, as shown in Figure 14a. The vorticity distribution of one cycle after the DSJ is turned on is shown in Figure 14b–e. It can be found that the vorticity increases significantly, and the vorticity distribution also changes significantly. Under the action of the DSJ, the fluid in the channel is affected by alternate blowing and suction, forming a periodic disturbance. A larger range of large vorticity distribution was observed in the channel. In addition, the reciprocating up and down the vibration of the diaphragm increases the vorticity in the cavity. It further shows that the fluid mixing under the action of the DSJ is enhanced, which enhances the convective heat transfer capability of the fluid.

#### *3.3. Entropy Generation Analysis*

The entropy generation rate is used to reflect the irreversibility of the heat transfer and flow process. This section studies the entropy generation characteristics of the DSJA from the perspective of the second law of thermodynamics. The total entropy production rate includes the entropy production rate caused by irreversible heat transfer and the irreversible flow loss caused by fluid friction.

**Figure 13.** Velocity distribution of channel's center section (*x* = 28.5 mm, partial enlarged view of jet orifice area) (DSJ on): (**a**) 0 T (1 T); (**b**) 0.25 T; (**c**) 0.5 T; (**d**) 0.75 T.

**Figure 14.** X vorticity distribution of channel's center section (*x* = 28.5 mm, overall view): (**a**) DSJ off; (**b**) 0 T (1 T); (**c**) 0.25 T; (**d**) 0.5 T; (**e**) 0.75 T (**d**,**e**, DSJ on).

The thermal entropy generation rate of the DSJ when off and the DSJ when on are shown in Figures 15 and 16. When the DSJ is off, the higher thermal entropy generation rate is mainly located at the interface between the fluid domain and the solid domain. This is mainly due to the large temperature gradient in these places. The water temperature at the inlet and outlet of the channel is low, the temperature gradient is large, and the thermal entropy generation rate is higher. When the DSJ is on, a high thermal entropy generation rate is also observed at the fluid–solid interface. The heat transfer in the region where the thermal boundary layer exists is poor and the thermal entropy generation rate is high. Under the disturbance of the DSJ, the distribution of the thermal entropy generation rate near the jet orifice is obviously disordered. At the same time, the impingement of the jet makes the thermal boundary layer thinner and the thermal entropy generation rate lower.

**Figure 16.** Thermal entropy generation rate in the fluid domain (*x* = 28.5 mm, 0.75 T of DSJ on).

The frictional entropy generation rate of the DSJ when off and the DSJ when on are shown in Figures 17 and 18. The higher frictional entropy generation rate is mainly located near the inlet and outlet of the channel and the jet orifice, especially in the jet impingement area. This is mainly due to the large velocity gradient at the inlet and outlet of the channel. The jet orifice diameter is small, and the velocity will change suddenly. When the DSJ is off, the jet deflects under the influence of incoming flow. The deflection area has a low-velocity area. Therefore, the frictional entropy generation rate is large. This is consistent with the phenomenon observed in Figure 11.

When the DSJ is on, Figure 18 observed the frictional entropy output value in a larger area than when the DSJ is off. This is because the diaphragm is at a limited position of downward vibration at this moment. The lower cavity ejects the jet and the upper cavity inhales the jet. The high strength of the jet results in a large velocity gradient. Observing the DSJ-off and DSJ-on cases, it can be found that the frictional entropy generation increases with the increase of the velocity of the DSJA. This is due to the increase in velocity gradient and the decrease in temperature. In the whole heat dissipation device, the greater the

favor greater velocity gradients.

**Figure 18.** Frictional entropy generation rate in the fluid domain (*x* = 28.5 mm, 0.75 T of DSJ on).

velocity, the more obvious the frictional entropy generation. Small channels and jet orifices

The total entropy generation rate of the DSJ when off and the DSJ when on are shown in Figures 19 and 20. It can be seen that the high total entropy generation rate is mainly located near the jet orifice and the inlet and outlet of the channel. We use Formula (12) to calculate the overall value of entropy, as shown in Table 2. It can be seen that the total entropy generation is dominated by the value of thermal entropy generation because the contribution of frictional entropy generation is very small. Since there is no velocity in the solid domain, its frictional entropy generation rate is 0. The total entropy generation in the fluid domain is much greater than that in the solid domain. Moreover, the total entropy generation of the DSJ when off is 0.01309, while the value of the DSJ when on is only 0.00414. The entropy production is reduced, and the latter is more than three times lower than the former. It indicates that the irreversible loss of the DSJA is greatly reduced.

**Figure 19.** Total entropy generation rate in the fluid domain (*x* = 28.5 mm, DSJ off).

**Figure 20.** Total entropy generation rate in the fluid domain (*x* = 28.5 mm, 0.75 T of DSJ on).

**Table 2.** Thermal entropy generation, frictional entropy generation, and total entropy generation in different areas of the cooling device.


#### *3.4. Effect of Diaphragm Frequency*

The influence of the diaphragm frequency on the performance of the cooling device was studied, keeping other settings unchanged.

The change of the periodic average temperature of the wall along the y direction along the channel with the diaphragm frequency is shown in Figure 21. The temperature distribution of the impact surface is not monotonic with the frequency. The impact surface temperature is the highest when the frequency is 10 Hz. When the frequency is 30 Hz, the impact surface has the lowest temperature and the temperature change is more obvious. The temperature near the jet impingement point is lower than that in other areas. This is because the effect of the dual synthetic jets strengthens the fluid-mixing and improves the convection heat transfer capacity of the fluid.

**Figure 21.** The change of the periodic average temperature of the wall along the y direction along the channel with the diaphragm frequency (*x* = 28.5 mm).

The channel inlet and outlet pressure drop changes with the diaphragm frequency as shown in Figure 22. The channel inlet and outlet pressure drop is the lowest when the diaphragm frequency is 10 Hz and the highest when the diaphragm frequency is 40 Hz. It shows that the pressure drop at the inlet and outlet of the channel is not monotonous with the change in the diaphragm frequency. It can be seen from Figures 21 and 22 that the cooling device has the best performance due to the optimal frequency. In the frequency range studied in this paper, the best comprehensive performance of the heat sink is obtained when the frequency is 30 Hz.

**Figure 22.** The channel inlet and outlet pressure drop changes with the diaphragm frequency (*x* = 28.5 mm).

#### **4. Conclusions**

In this paper, combined with the DSJ technology, a forward liquid cooling device is designed. The heat transfer and flow characteristics of the device are analyzed. The entropy generation of the DSJA is analyzed from the point of view of the second law of thermodynamics. The main conclusions are as follows


**Author Contributions:** Writing: original draft, Y.K.; Writing: review and editing, Z.L., X.D., Y.Z. and Z.X. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China, grant numbers 11872374, 11972369; the Science and Technology Innovation Program of Hunan Province, grant number 2021RC3075; and the Natural Science Foundation of Hunan Province of China, grant number 2021JJ40672.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Hongxin Wang 1,2, Degang Xu 1,\*, Linwen Li 2, Kaiwen Zhou 3, Xin Wen 3,\* and Hui Tang <sup>4</sup>**

<sup>1</sup> School of Automation, Central South University, Changsha 410083, China


**Abstract:** This paper aims to further the understanding of the mixing process of in-line twin synthetic jets (SJs) and their impact in the near-wall region in a flat-plate laminar boundary layer. A numerical study has been carried out, in which colored fluid particles and the *Q* criterion are used to track the SJ-induced vortex structures at the early stage of the evolution. Interacting vortex structures at four selected phase differences are presented and analyzed. It is found that the fluid injected at the early stage of the blowing stroke mainly contributes to the formation of the hairpin legs, the fluid injected near the maximum blowing mainly contributes to the formation of the hairpin head, and the fluid injected at the late stage of the blowing stroke contributes very little to the formation of the hairpin vortex. It is also confirmed that, irrespective of the phase difference, the hairpin vortex issued from the upstream actuator is more capable of maintaining its coherence than its counterpart issued from the downstream actuator. The influence of the interacting vortex structures on the boundary layer is also studied through investigating excess wall shear stress. In all cases, a pair of streaks of high wall shear stress can be observed with similar size. Among them, the streaks have the strongest wall shear stress, with the largest gap at phase difference 0 when partially interacting vortex structures are produced. The findings can provide valuable guiding information for the applications of synthetic jets in heat transfer, mixing control, and flow control in a crossflow.

**Keywords:** twin synthetic jets; laminar boundary layer; vortex interaction; simulation

#### **1. Introduction**

A synthetic jet (SJ) actuator can produce a train of vortical structures by the periodical motion of an oscillatory diaphragm on a cavity [1–3]. These vortical structures can effectively enhance fluid mixing in the crossflow, thereby providing an effective mechanism of boundary layer flow control [4–9].

Previous works were carried out to investigate the impact of the SJ-induced vortical structures in the near-wall region of a boundary layer, where effective flow separation control and wall heat transfer are ultimately desired. Using dye visualization, Zhong et al. [10] and Jabbal and Zhong [11] looked at more detailed interactions and identified three key SJ-induced vortex structures, i.e., hairpin vortices, stretched vortex rings, and tilted vortex rings, under different jet operation and boundary layer flow conditions. Different thermal footprints on the wall are correlated to the vortical structures induced by SJs. Using high-frame-rate particle image velocimetry (PIV) experiments, Jabbal and Zhong [12] confirmed that the wall shear stress shares similar patterns with the thermal footprints, and found that hairpin vortices and stretched vortex rings are able to offer enhanced heat transfer in the near-wall region. Using numerical methods, Lee et al. [13] confirmed the effective wall cooling of single SJs in micro-channels. They also found that

**Citation:** Wang, H.; Xu, D.; Li, L.; Zhou, K.; Wen, X.; Tang, H. Evolution and Near-Wall Effect of the Vortex Structures Induced by In-Line Twin Synthetic Jets in a Crossflow. *Actuators* **2022**, *11*, 234. https:// doi.org/10.3390/act11080234

Academic Editor: André Preumont

Received: 22 July 2022 Accepted: 12 August 2022 Published: 16 August 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the flow structures induced by stronger SJs in the micro-channel can exert larger mixing regions and stronger convection at the heated wall surface.

To increase the impact region, SJ arrays are usually applied. For example, in the application of fluid cooling, Lee et al. [14] found that twin SJs in out-of-phase operation can significantly enhance wall surface cooling due to the enhanced fluid mixing. However, unlike single SJs, the interactions between multiple circular SJs and cross-flows, as well as among multiple SJs, have not been fully studied and understood. Very limited investigations have been reported on this topic in the literature, especially on the near-wall effect. Based on oil visualization, Watson et al. [15] examined the in-teraction of two neighboring SJs in a boundary layer flow, and pointed out that the twin SJs can perform in either a constructive (enhancing the strength of the SJs) or a destructive (weakening the SJs) manner under different configurations. Using hot-wire measurements, Liddle and Wood [16] and Liddle et al. [17] confirmed that the phase difference plays a very important role in the interaction of twin SJs in a crossflow. At a 2700 phase difference, the SJ actuators in the upstream and downstream locations yielded a more concentrated movement of near-wall fluid towards the vortex core than at a 900 phase difference. With PIV measurements of in-line twin SJs in a channel flow, Honami and Motosuke [18] found that the twin SJs can have different vortex strengths with different configurations. Although these studies revealed useful physical insights, they are not enough to give a holistic picture of the flow structures induced by multiple circular SJs in a crossflow and their near-wall effect. Therefore, more detailed studies are highly desired.

Although these works have revealed some useful physical insights into the interaction of in-line twin SJs in boundary layers, the detailed vortex interaction and evolution of twin SJs are not well understood. Limited by the spatial and temporal resolution, it is difficult for current experimental methods to capture the three-dimensional flow structures of twin SJs. In addition, it is also difficult to examine the evolution of the interacting flow structure issued from different SJ actuators in Cartesian coordinates. Therefore, two research questions still remain: First, how does the mixing process of the twin-SJ-induced flow structures occur in detail? Second, what is the influence of these flow structures on the boundary layer flow, especially in the near-wall region? By examining the flow structures and their impact in the boundary layer, we can obtain results which are valuable for flow control and heat transfer applications. For example, strong vortices and high wall shear stress can produce high flow mixing in the boundary layer, delaying flow separation and enhancing heat transfer on the wall. Therefore, to address these questions in the present research, a numerical study is conducted. A cluster of colored fluid particles are used to track the evolution and mixing process of the single and twin SJs at the early stage. In the Lagrangian perspective, fluid particles are applied to examine the flow dynamical system rather than a continuum [19]. It is widely used to present high-quality visible details of time-dependent vortex structures, which may not be apparent in the classical Eulerian fields [20,21]. In addition, to capture the associated vortex structures, *Q* iso-surfaces are defined [22]. Finally, the excess wall shear stress is eval-uated and compared.

#### **2. Methodology**

#### *2.1. Problem Description*

The twin SJ actuators used in the present numerical study are identical to those used in our previous experiment [23,24]. The SJ cavity has an orifice with a diameter of *Do* = 5 mm. The centers of the two SJ actuator orifices are at a distance of *d* = 10 mm. Similar to previous experiments [23], the oscillation diaphragm of each actuator has a diameter of *Dd* = 45 mm, which oscillates like a piston, with velocity varying in a sinusoidal manner.

$$
\mu\_d(t) = \pi \Delta f \sin(2\pi ft) \tag{1}
$$

where Δ is the peak-to-peak displacement of the diaphragm and f the oscillation frequency. During an actuation cycle, the beginning of blowing is defined as *t*/*T* = 0 (and 1). In the

present study, the working flow is water, which is incompressible. According to Tang and Zhong [25], the jet velocity *Uo* averaged over an entire actuation cycle can be calculated as

$$\overline{\mathcal{U}}\_o = \frac{1}{T} \int\_0^{T/2} u\_o(t)dt = \Delta f \left(\frac{D\_d}{D\_o}\right)^2\tag{2}$$

where the instantaneous jet velocity is *uo*(*t*).

In current study, the bulk flow has a velocity of *U*∞ = 0.055 m/s. The bulk flow has a boundary layer thickness of 15 mm, i.e., δ = 3*Do*, at the actuator orifices. Based on Reynolds number, the momentum thickness is *Re*<sup>θ</sup> ≈ 106. The actuators have a diaphragm displacement of Δ = 0.158 mm with a frequency of *f* = 1 Hz. To characterize the performance of SJs in a laminar boundary layer, the dimensionless stroke length *L* and velocity ratio *VR* are also defined according to Wen and Tang [24]. The stroke length *L* is defined as

$$L = \frac{L\_o}{D\_o} = \frac{\overline{\mathcal{U}}\_o}{fD\_o} \tag{3}$$

where *Lo* is the jet column length during the SJ blowing stroke. The jet strength is characterized by the velocity ratio *VR* as

$$VR = \frac{\mathcal{U}\_o}{\mathcal{U}\_\infty} \tag{4}$$

From Equations (2)–(4), the settings in the present study yield *L* = 2.6 and *VR* = 0.23, to produce hairpin vortices according to the parametric map revealed by Jabbal & Zhong [11]. These hairpin vortices can stay in the boundary layer and hence are very desired in SJ-based flow separation delay.

In the present study, the two SJ actuators are arranged in line with the crossflow and are separated by a fixed distance *d* = 2*Do*. Phase difference between the two SJs are investigated by Δ*φ* = 0, *π*/2, *π* and 3*π*/2 (or −*π*/2). Here Δ*φ* = 0 means that the twin SJs are in phase, and the Δ*φ* = *π*/2 means the downstream actuator is *π*/2 behind, and so on.

#### *2.2. Numerical Approaches*

Figure 1 shows the computational geometry and boundary conditions, which are similar to those used in Wen and Tang [26] and Zhou and Zhong [27]. In the *x* (streamwise) direction, the simulation domain has a length of 55Do. In the *y* (wall-normal) and *z* (spanwise) directions, the domain is 30*Do* and 20*Do*, respectively. The origin of the coordinate system is set at the midpoint of two actuator orifice centers, which are 15*Do* downstream from the crossflow inlet. Although not presented in Figure 1, the cavities of the two SJ actuators are included in the simulation. The crossflow in-let on the left of the computational domain is set as the velocity inlet with a laminar Blasius velocity profile. The wall of flat plate in the boundary layer and the SJ actuators are set as no-slip conditions. Both lateral sides of the crossflow domain are set as periodic conditions. The upper face and the downstream boundary on the right of the crossflow domain are set as the pressure outlet. On the diaphragms of the actuators, a time-dependent velocity condition is set according to Equation (1).

To activate the SJ actuator, the motion of the diaphragm in the cavity is controlled by setting a time-dependent velocity boundary condition according to Equation (1). In this study, water is used, which is incompressible. The fluid is governed by mass con-servation and momentum conservation. The flow is laminar; therefore, no turbulent model is used. The second-order implicit scheme is used in time, and the second-order upwind scheme is used in space. For the pressure–velocity coupling, the pressure-implicit with splitting of operators method is used. The unsteady, three-dimensional incompressible Navier–Stokes equations are solved by a commercial CFD code, ANSYS FLUENT 2020 R2.

$$\frac{\partial \mu\_i}{\partial \mathbf{x}\_i} = 0 \tag{5}$$

$$\frac{\partial \mu\_i}{\partial t} + \mu\_j \frac{\partial \mu\_i}{\partial \mathbf{x}\_j} = -\frac{1}{\rho} \frac{\partial p}{\partial \mathbf{x}\_i} + v \frac{\partial^2 \mu\_i}{\partial \mathbf{x}\_j \partial \mathbf{x}\_j} \tag{6}$$

where ρ is the fluid density and ν is viscosity, and the summation notion is applied. Sensitivity studies have been conducted to ensure the mesh number and time-step size are optimum. Three different mesh numbers, i.e. 2.1 × 106, 2.8 × 106, and 3.7 × 106, and three different time-step sizes, i.e. 1/240 s, 1/360 s, and 1/480 s, are examined. Mesh and time sensitivity studies have been conducted for the case of twin SJs operating in phase (i.e., Δφ = 0). With different mesh sizes and time steps, the time-averaged velocity profiles along a spanwise, wall-parallel line (i.e., *x*/*Do* = 5 and *y*/*Do* = 2) within the boundary layer are compared. As shown in Figure 2, the simulation result seems independent of the selected mesh size and time step, confirming the mesh and time convergence in the current study. The total number of meshes used in the simulation is about 2.8 million. In one diaphragm oscillation cycle, there are 240 time steps as a compromise between computational accuracy and time.

**Figure 1.** Computational geometry and boundary conditions for current CFD simulations. The actuator cavities are included in the simulation (although the cavities are not presented). The orifices of the two actuators are presented as dash circles located close to the origin of the coordination.

**Figure 2.** Comparison of the time-averaged velocity profiles along a spanwise, wall-parallel line (i.e., *<sup>x</sup>*/*Do* = 5 and *<sup>y</sup>*/*Do* = 2) within the boundary layer for the (**a**) mesh size (chosen value 2.8 <sup>×</sup> 106) and (**b**) time step (chosen value 1/240) sensitivity studies for the twin SJs operating in phase.

In the present study, the SJ-induced vortex structures are identified using the *Q* criterion [22]. A vortex structure can be captured by an iso-surface of a positive *Q* value. To differentiate individual SJs in the present simulation, particles of two different colors are continuously and uniformly released from the two actuator orifices during the blowing cycle once the flow field reaches its steady state. The density of these particles is set to that of water, and the discrete particle model (DPM) is used to track them. Similar to our dye visualization experiment [23,24], the particles issued from the upstream actuator are colored in red, whereas those from the downstream actuator are in green.

To reflect the capability of the SJ-induced vortex structures in flow separation control, wall shear stress is calculated based on the least squares cell-based gradient discretization proposed by Anderson and Bonhaus [28] to ensure a second-order interpolation, with a spatial resolution of 0.5% of the boundary layer thickness or *y*+ = 0.2.

#### *2.3. Validation*

The present CFD framework has been well validated with a single SJ in our previous work [26]. To further validate this framework for the current twin-SJ configuration, experiments are conducted using color dye visualization and PIV measurement. In dye visualization, food dye with two different colors fills the two actuator cavities, i.e. red color for the upstream actuator and green color for the down-stream. Methanol is premixed with the dye to achieve a similar density to water. Two cameras are used to capture images from the side and bottom views. In the PIV measurement, the Dantec system is used. Polyamide seeding particles are used to track the flow, and a FlowSense 2M charge-coupled-device camera is used to capture the particle images. A double pulsed Nd:YAG laser is used to illuminate the seeding particles. The sampling frequency of PIV measurement is 10 Hz. For each case, 300 image pairs are taken to achieve convergence. The PIV measurement results are verified with Blasius solution. More details of the dye visualization and PIV measurement can be found in our previous work [23]. To validate the CFD framework, flow patterns represented by colored fluid particles in the simulation are first compared to those in the color dye visualization experiment. As shown in Figure 3, the simulated particle patterns match with the dye visualization in a good manner. As shown in Figure 4, instantaneous and time-averaged velocity profiles are also compared at *x* = 4*Do*, with a good agreement between the simulation results and the PIV data.

**Figure 3.** CFD Particle patterns (**left**) and experimental dye visualization (**right**) for the twin SJs at Δ*φ* = *π*/2 and *t*/*T* = 0.25 (The lines indicate the edge of the boundary layer).

**Figure 4.** Comparison of instantaneous and time-averaged streamwise velocity profiles taken at *x* = 4*Do* from the CFD and PIV data for twin SJs operating at Δ*φ* = 0. The *t*/*T* = 0.45, 0.7, 0.95 profiles and the mean profile are shifted by *u*/*U*∞ = 0.5, 1.0, 1.7 and 2, respectively.

#### **3. Evolution and Interaction of In-Line Twin SJ-Induced Flow Structures**

#### *3.1. Flow Structures of a Single SJ in a Crossflow*

The single-SJ-induced flow structures visualized by the fluid particles, and the *Q* iso-surfaces are presented first to serve as the baseline case. As shown in Figure 5a, iso-surfaces of *Q* = 5 are superimposed over the particles, showing that the flow structures represented by the fluid particles coincide with the primary hairpin vortices captured by the *Q* iso-surfaces in terms of the shape, size, and location. This confirms that, in addition to the *Q* iso-surfaces, the SJ-induced hairpin vortices can also be examined through the pattern and motion of fluid particles. As shown in Figure 5b, where only the fluid particles are presented, each hairpin vortex is well represented by a strong head and a pair of slim, stretched legs. The hairpin head also carries a spiral of fluid particles when seen from the side cross-section view, indicating very strong rotation inside. However, no particles are observed in the induced secondary vortices that are captured by the *Q* iso-surfaces as shown in Figure 5a. This is not surprising because the particles are all released from the actuator orifice, while the secondary vortices form with ambient fluid that does not carry any particles. These secondary vortices are very close to the wall and hence play an important role in SJ-related flow separation control [26].

**Figure 5.** (**a**) Fluid particles plus *Q* = 5 iso-surfaces and (**b**) fluid particles for a single SJ at *t*/*T* = 0.

#### *3.2. Evolution of Single and In-Phase Twin SJs*

By using a sequence of snapshots of fluid particles, the evolution of the single SJ is presented in the upper two rows of Figure 6. It can be seen that, after expulsion from the actuator, the SJ-induced vortex ring gradually tilts downstream because of the velocity gradient in the boundary layer. Its upstream part is weakened by the boundary layer. As time advances, a hairpin vortex is formed by the shear stress in the boundary layer. It should be noted that the hairpin vortex is formed in a boundary layer with a flat plate. The wall-surface shape and obstacles in the boundary layer can significantly affect the evolution of hairpin vortices and horseshoe vortices, which is comprehensively discussed in previous works [29–31]. To further reveal the contribution of various SJ portions to the resulting hairpin vortex, three groups of fluid particles are marked in the snapshots, which are released from the orifice exit at three consecutive instants, i.e., *t* = 0.0625*T* (group one), 0.25*T* (group two), and 0.4375*T* (group three). Each group consists of four selected particles, representing the upstream, central, downstream, and side portions of the SJ cross-section, respectively. By tracking these three groups of particles, it is found that the hairpin head contains fluid particles mainly from group two, while the hairpin legs contain fluid particles mainly from group one. As for the fluid particles in group three, most of them are inhaled back into the actuator during the SJ ingestion process (*t* = 0.5*T*~*T*). This result indicates that the fluid injected at the early stage of the blowing stroke (t ≈ 0) mainly contributes to the formation of hairpin legs, the fluid injected near the maximum blowing (*t* ≈ 0.25*T*) mainly contributes to the formation of the hairpin head, and the fluid injected at the late stage of the blowing stroke (t ≈ 0.5*T*) contributes very little to the formation of the hairpin vortex.

**Figure 6.** Evolution of fluid particles of the single SJ (**upper** two rows) and in-phase twin SJs (**lower** two rows). Group 1 particles (in pink) are released at *t*/*T* = 0.0625, Group 2 particles (in yellow) are released at *t*/*T* = 0.25, and Group 3 particles (in blue) are released at *t*/*T* = 0.4375. All particles are released at the actuator orifice exit.

If another SJ is introduced near the single SJ, the interaction between these two SJs will result in the formation of new flow structures. Take the in-phase, in-line twin SJs as an example, whose evolution is presented in the lower two rows of Figure 6. The SJ issued from the upstream actuator (in red, also called the trailing vortex in this case) is similar to the single SJ before it entrains its downstream counterpart, while the SJ issued from the downstream actuator (in green, also called the leading vortex) is significantly affected by its upstream counterpart. The interaction between these two SJs is evidenced by the lift-up of the leading vortex legs and their entrainment into the trailing vortex. Furthermore, since it is the leading vortex that is significantly affected by the interaction, the same three groups of fluid particles as in the single-SJ case are released from the downstream orifice exit. By tracking these particles and comparing their traces with those in the single SJ, it can be seen that during the blowing stroke (*t* = 0~0.5*T*), the evolution of these particles is similar to that in the single-SJ case, except that the particles released from the downstream portion of the SJ cross-section at the early stage of the blowing stroke (marked as "1D") are stripped out from the main vortex structure under the influence of the preceding vortices. When the twin SJs start interacting with each other in the ingestion period (*t* = 0.5*T*~*T*), particles originally in the lower portion of the vortex leg, e.g. the particles "1U" and "3C", are significantly lifted up by the trailing vortex, while particles originally in the vortex head, e.g. the particles in group two, still remain close in the highly rotating vortex head, indicating a well-maintained coherence in the head even under the strong interaction.

#### *3.3. Impact of Phase Differences*

With the simulated colored fluid particles, three different vortical structures are observed at the four selected phase differences. At Δφ = π/2 (Figure 7b), the two SJs merge together to form one combined vortex, which is much stronger than single hairpin vortices and can pass through the boundary layer quickly. At Δφ = 3π/2 (Figure 7d), the two hairpin vortices are separate from each other without obvious interaction, forming two completely separated vortices. At Δφ = 0 (Figure 7a) and π (Figure 7c), the legs of one hairpin vortex are entrained into the head of the other, forming partially interacting vortex structures. As such, the twin SJs have distinct vortex strengths at various phase differences, as revealed in the authors' previous PIV experiments [23,24]. Furthermore, it is seen that, no matter which type of interaction occurs, the green particles issued from the downstream actuator are mainly scattered in the outer layer of the interacting structures, whereas the red particles is-sued from the upstream actuator are mainly enclosed inside. The different particle distributions can be caused by the different strengths of the vortical structures issued from the twin SJ actuators. The flow structures issued from the upstream actuator can hold vortical structures better than those from the downstream. Therefore, the parti-cles issued from the upstream actuator are more concentrated inside of the interacting flow structures. This is especially so when the interacting structures develop into the second actuation period, as shown in the plane view of Figure 7a–c. It is also seen that at Δφ = 3π/2 (Figure 7d), unlike its counterpart, the head of the SJ issued from the upstream actuator is not affected at all by the near-wall particles. All these observations support our previous statement, i.e., in any case, the SJ issued from the up-stream actuator can better maintain its vortical structure than the SJ issued from the downstream actuator [23,24], and hence it dominates the interaction. The reason is obvious: the vortex produced by the upstream actuator has sufficient time to develop itself before the interaction, while the SJ issued from the downstream actuator is disturbed from the very beginning of its emergence.

The different capability of the two hairpin vortices in maintaining their own coherence can be further demonstrated in snapshots from the streamwise cross-section view. Figure 8 shows the snapshots from the Δφ = 0, π/2 and π cases, at the same instant *t*/*T* = 0 as in Figure 7. In all three snapshots, the boundaries between the particles from different actuators are clear. At Δφ = π/2 (Figure 8b), the head of the combined vortex includes two cores of fluid particles issued from the upstream actuator (in red) which are surrounded by scattered fluid particles issued from the downstream actua-tor (in green). At Δφ = 0, the twin SJs produce partially interacting vortex structures (Figure 8a). The fluid particles issued from the upstream actuator still appear in two cores, forming the head of the trailing vortex, which is surrounded by the fluid particles issued from the downstream actuator. Moreover, two legs of the leading vortex with green particles appear near the mid-span

plane, whose induced flow pushes the upper portion of the trailing vortex heads (two red cores) away from each other, which is prominent if compared with the combined vortex head shown in Figure 8b. In the partially interacting vortex structures at Δφ = π (Figure 8c), the head of the trailing vortex (in green) no longer appears as two solid cores. Instead, it appears as two hollow tubes enclosing the two counter-rotating legs of the leading vortex with red particles, and its shape differs greatly from those in Figure 8a,b.

**Figure 7.** Colored fluid particles of the in-line twin SJs (*t*/*T* = 0), (**a**) Δ*φ* = 0, (**b**) Δ*φ* = *π*/2, (**c**) Δ*φ* = *π*, (**d**) Δ*φ* = 3*π*/2.

**Figure 8.** Snapshots of colored fluid particles from the streamwise cross-section view at (**a**) *x* = 11*Do* for the Δ*φ* = 0 case, (**b**) *x* = 11*Do* for the Δ*φ* = *π*/2 case and (**c**) *x* = 9.5*Do* for Δ*φ* = *π* case (*t*/*T* = 0). These streamwise cross-section planes are indicated with dashed lines in Figure 7.

Particles issued from selected positions of the SJ actuator orifices are also tracked to study the interactions of SJ vortex structures, as shown in Figure 9. The single SJ and the twin SJs at Δφ = π/2 and 3π/2 are selected for presentation. Compared to those in the single SJ, the particles issued at the center of the orifices at Δφ = π/2 are lifted to a much higher position at x = 3Do outside of the boundary layer. This is because a stronger vortex is merged from the twin SJs. The two particles issued from different actuators have strong interactions and rotate around the merged vortex head. In addition, the particles issued at *z*/*Do* = 0.4 have a slightly faster rotation than that in a single SJ with a shorter traveling distance. On the other hand, the trajectories of particles at Δφ = 3π/2 are similar to those in the single SJ because the two separated hairpin vortices are produced, with the exception that the particle issued at z/Do = 0.4 from the downstream orifice is trapped in the tip of the hairpin leg with a straighter trajectory.

**Figure 9.** Tracking of particles that are released from different positions of the actuator orifices, top row: at center, middle & buttom rows: at *z*/*Do* = 0.4 (side and plan views)). Vortex structures are also added at *t*/*T* = 0.25 indicated by the particles with larger diameters.

Figure 10 shows the correlations of the streamwise and wall-normal velocity components at three selected stations in the mid-span plane. In the single-SJ case, the upwash flow of hairpin vortices can increase the velocity in the wall-normal direction and decrease the velocity in the streamwise direction in the boundary layer [23,24]. For the twin SJs at Δφ = π/2, the disturbance of the boundary layer is much stronger, indicated by the larger correlation loop of the streamwise and wall-normal velocity components. The shape of correlation loop is similar to that in the single SJ, because only the stronger hairpin vortices are merged with similar vortex structures. On the other hand, the shape of the correlation loop becomes very different at Δφ = 3π/2. At a phase difference of Δφ = 3π/2, the SJ flow structures are produced alternately from the twin actuators with an almost equal time interval. Therefore, the twin SJs have a weak interaction but with a doubled frequency compared to single SJs. On the other hand, the twin SJs merge into a combined vortex at Δφ = π/2 with a high strength but the same frequency to the single SJs. Differing from the single loop in the single-SJ case, two loops are observed with substantial differences. This indicates that the two separated hairpin vortices at Δφ = 3π/2 also have notable differences in vortex strength. In addition, the smaller loops also imply that the hairpin vortices at Δφ = 3π/2 are slightly weaker than that in the single-SJ case.

**Figure 10.** Correlations of streamwise and wall-normal velocity components in the mid-span plane at three selected positions, top: (*x* = 5*Do*, *y* = 3.4*Do*), middle: (*x* = 5*Do*, *y* = 2.4*Do*), bottum: (*x* = 5*Do*, *y* = 1.4*Do*).

#### **4. Impact of SJ Vortices on Wall Shear Stress**

Previous studies have revealed that the SJ-induced near-wall streamwise vortices, i.e., the hairpin legs and induced secondary vortices, play an important role in flow separation control or enhancement of heat and mass transfer [12]. In this section, in order to reveal details of the evolution of these streamwise vortices, iso-surface of a smaller value of *Q* = 1 are colored by streamwise vorticity. To examine the impact of the SJ vortices on the boundary layer, the excess wall shear stress is presented and calculated as (*τ*<sup>w</sup> − *τ*w,nojet)/*τ*w,nojet, which is a useful indicator of the SJ effectiveness in heat transfer enhancement [12].

#### *4.1. Instantaneous Influence*

Figure 11 shows the single-SJ-induced vortices and the corresponding contours of excess wall shear stress at the instant *t*/*T* = 0. Evolution of the primary hairpin vortices and secondary vortices has been discussed in our previous study [26]. These vortex structures influence the boundary layer as well as the wall shear stress. A pair of streaks of positive excess wall shear stress is observed outside of the SJs. This streak pair becomes wider near the regions between the legs of hairpin vortices and the secondary vortices (at about 8 < *x*/*Do* < 10), where the downwash flow is strengthened by the interaction of the vortices.

**Figure 11.** Vortex structures in the single-SJ case, captured using *Q* = 1 iso-surfaces (**a**) and excess wall shear stress (**b**) at *t*/*T* = 0.

The two SJs merge into a stronger vortex when their phase difference is Δ*φ* = *π*/2. As shown in Figure 12a, both the hairpin vortices and the induced vortices, including the secondary vortices and horseshoe vortex, become stronger compared to those induced by the single SJ. The evolution of these vortices presented in the upper two rows of Figure 12b reveals that, at *t*/*T* = 0.25 when the upstream SJ is at its maximum blowing, a vortex ring (labeled as A) is issued the upstream actuator. It lifts up the horseshoe vortex produced in the previous cycle. In the downstream part of this horseshoe vortex a pair of secondary vortices begin to form. At *t*/*T* = 0.5, the vortex ring A begins to evolve into a hairpin vortex due to its interaction with the crossflow. Meanwhile a new vortex ring (labeled as B) emerges at the downstream orifice, breaking the horseshoe vortex. In addition, the downstream secondary vortices become more obvious. At *t*/*T* = 0.75 & 1, the twin SJs combine into a stronger hairpin vortex, which is accompanied by stronger secondary vortices and a larger horseshoe vortex.

The enhancement of wall shear stress is more significant, especially near to the orifice (*x/Do* < 12) in Figure 12a. However, the combined vortices are so strong that they quickly penetrate the boundary layer, resulting in significant decay of the excess wall shear stress in the far downstream. In addition, some deficit of the wall shear stress is found outboard of the secondary vortices, which is associated with the induced upwash flow [26]. Detailed examination in the near-orifice region shows that changing of the excess wall shear stress is obviously periodic during one actuation cycle, as shown in the bottom row of Figure 12b. At *t*/*T* = 0.25, besides the area around the upstream orifice, high excess wall shear stress is also found near the downstream orifice, which is induced by the legs of the preceding combined vortex and the associated secondary vortices. At *t*/*T* = 0.5, high excess wall shear stress is induced around the downstream orifice due to the emergence of the vortex ring B. The excess wall shear stress downstream of this orifice becomes weaker. At *t*/*T* = 0.75, as the two SJs begin to combine, the regions of high excess wall shear stress around the two orifices connect. At downstream of the orifices, a large area of negative excess wall shear stress forms between the new forming vortex and the preceding one. At *t*/*T* = 1, as the newly formed combined vortex convects downstream, positive excess wall shear stress recovers near the orifices.

**Figure 12.** Vortex structures and excess wall shear stress for the twin SJs operating with Δ*φ* = *π*/2. (**a**) Vortex structures (**left**) and excess wall shear stress (**right**) at *t*/T = 0. (**b**) Evolution of vortices (**upper** two rows) and corresponding excess wall shear stress (**bottom** row). The color map can be referred to Figure 11.

Two completely separated vortices are produced at Δ*φ* = 3*π*/2. Primary hairpin vortex and outboard secondary vortices are induced as shown in left subplot of Figure 13a. However, their evolution processes start at different locations and instants. As shown in the upper two rows of Figure 13b, a pair of secondary vortices (labeled as #1) begins to form between the two orifices at *t*/*T* = 0.25. At *t*/*T* = 0.5, these secondary vortices become more obvious, as the horseshoe vortex is eliminated by the suction of the downstream actuator. They further grow at *t*/*T* = 0.75 and 1, and are sandwiched by the head of the SJ issued from the upstream actuator and the legs of the SJ issued from the downstream actuator. At *t*/*T* = 1, a vortex ring emerges from the downstream orifice and begins to evolve into a hairpin vortex at *t*/*T* = 1.25 (equivalently *t*/*T* = 0.25). By interacting with the preceding hairpin legs, it induces a pair of secondary vortices (labeled as #2) downstream. This pair of secondary vortices grows at *t*/*T* = 1.5 (equivalently *t*/*T* = 0.5), and finally becomes comparable to its counterpart (i.e., #1) at *t*/*T* = 1.75 (equivalently *t*/*T* = 0.75). Therefore, both the primary hairpin vortices and the induced secondary vortices double their appearance frequency (in other words, double their number) at far downstream. On

the other hand, they are all shorter than those in single SJ case. The strip pair of enhanced wall shear stress also reflects this doubled frequency, as shown in the right subplot of Figure 13a. The strips are narrower compared to single-SJ case. Also, due to the doubled frequency effect, the strips in the near-orifice region can maintain the high excess wall shear stress, as revealed in the bottom row of Figure 13b.

**Figure 13.** Vortex structures and excess wall shear stress for the twin SJs operating with Δ*φ* = 3*π*/2. (**a**) Vortex structures (**left**) and excess wall shear stress (**right**) at *t*/T = 0. (**b**) Evolution of vortices (**upper** two rows) and corresponding excess wall shear stress (**bottom** row). The color map can be referred to Figure 11.

Partially interacting vortex structures are produced at Δ*φ* = 0 and *π*. The *Q* isosurfaces reveal that, no matter whether it is the trailing vortex (left subplot of Figure 14a) or leading vortex (left subplot of Figure 15a), the vortex produced from the upstream actuator dominates the interacting vortex structure in terms of vortex coherence. This is consistent to what has been revealed with the colored fluid particles in Section 3.3. For this reason, the vortex produced from the upstream actuator is also called the *dominant hairpin vortex* hereinafter, despite the phase difference. In both cases, a pair of much longer secondary vortices is induced by the interacting vortex structure, compared to the single-SJ case. Although appearing similar, these two partially interacting vortex structures show an interesting difference. That is, the lateral distance between the two legs of the partially

interacting SJs at Δ*φ* = 0 is much larger than that at Δ*φ* = *π*. This is because in both cases the legs of the interacting vortex structure are from the dominant hairpin vortex. At Δ*φ* = 0, the leading vortex is entrained through the gap between the legs of trailing dominant vortex. Therefore, the gap is enlarged during the process. On the other hand, at Δ*φ* = π, the legs of dominant leading vortex is pushed together by the trailing vortex, resulting in a smaller gap. This difference also affects the pattern of excess wall shear stress. Although the patterns in both cases are similar to that in the Δ*φ* = *π*/2 case in far downstream, the deficit in the centerline low stress region in the Δ*φ* = 0 case is significantly larger (right subplot of Figure 14a). Furthermore, at Δ*φ* = 0, the side view reveals that the legs of the dominant hairpin vortex are gradually "pulled" down toward the wall as the reaction to their entrainment of the leading hairpin legs, which are then able to exert more influence on the wall shear stress, especially in the far-downstream region. As such, the two high stress streaks in the mid- and far-downstream regions are stronger than those in the near-downstream region, which is quite different from what is observed in the other cases.

**Figure 14.** Vortex structures and excess wall shear stress for the twin SJs operating with Δ*φ* = 0. (**a**) Vortex structures (**left**) and excess wall shear stress (**right**) at *t*/T = 0. (**b**) Evolution of vortices (**upper** two rows) and corresponding excess wall shear stress (**bottom** row). The color map can be referred to Figure 11.

**Figure 15.** Vortex structures and excess wall shear stress for the twin SJs operating with Δ*φ* = *π*. (**a**) Vortex structures (**left**) and excess wall shear stress (**right**) at *t*/T = 0. (**b**) Evolution of vortices (**upper** two rows) and corresponding excess wall shear stress (**bottom** row). The color map can be referred to Figure 11.

#### *4.2. Time-Averaged Influence*

The time-averaged influence of the single and twin SJs on the wall shear stress is compared in Figure 16. It is seen that about 1*Do* downstream of each operating actuator there exists an obvious region of negative excess stress, which is caused by the blockage of the SJs and the roll-up of the downstream branch of the SJ-induced vortex rings. Further downstream, the averaging removes the periodic patterns of the two streaks in all cases. Regardless the phase difference, the twin SJs can generate significantly more excess stress than the single SJ. But the persistency of the influence varies with phase difference. At Δ*φ* = *π*/2 the combined vortex structures induce a pair of wall shear stress streaks that is strong only within a short distance downstream from the orifices and starts decaying further downstream (see Figure 16c). This is because the combined vortices can penetrate

the boundary layer deeply. The vortices move away from the wall as they travel in the crossflow, resulting in a fast decaying influence on the wall in the far-downstream region. This is evidenced by the particle tracking (Figure 7b) and instant flow structures (Figure 12) On the contrary, at Δ*φ* = 3*π*/2 the two completely separated vortices can stay close to the wall (Figures 7d and 13), and produce a pair of high stress streaks that sustain their strength for a much longer distance. The partially interacting vortices at Δ*φ* = 0 induce two streaks of high wall shear stress with gradually increased strength as shown in Figure 11b, between which a much wider strip of low stress value is observed. However, the time averaged high stress streaks produced by the other partially interacting vortices at Δ*φ* = *π* (Figure 16d) appears similar to that in the Δ*φ* = *π*/2 case.

**Figure 16.** Time-averaged excess wall shear stress, (**a**) single SJ, (**b**) Δ*φ* = 0, (**c**) Δ*φ* = *π*/2, (**d**) Δ*φ* = *π*, (**e**) Δ*φ* = 3*π*/2.The color map can be referred to Figure 11.

In Figure 17, the excess wall shear stress is time-averaged and compared at two positions, i.e., *x* = 4*Do* and 18*Do*. In all cases, the two high stress streaks are represented by two peaks. At *x* = 4*Do*, the two peaks induced by the single SJ are symmetrically located about 0.8*Do* away from the mid-span plane, whereas the two peaks induced by all the twin SJs are slightly closer. It is also seen that the Δ*φ* = 0 case has a peak value close to that in the single SJ case. On the other hand, the other three twin SJ cases all have a much larger peak value (at least 20% higher). This is because in the Δ*φ* = 0 case the legs of the leading vortex (issued from the downstream actuator) are lifted up immediately after their emergence (as shown in Figures 6 and 7a), and hence the near-field wall shear stress is only affected by the legs of the trailing vortex (issued from the upstream actuator), while in the other three twin-SJ cases, i.e., Δ*φ* = *π*/2, *π* and 3*π*/2, the near-field wall shear stress is affected by the legs of both vortices.

**Figure 17.** Time-averaged excess wall shear stress along spanwise direction at (**a**) *x* = 4*Do* and (**b**) *x* = 18*Do*.

At *x* = 18*Do* that is far from the two orifices, due to the spreading of the vortex structures and their penetrating in the boundary layer, the two stress peaks in all cases are separated from each other with a larger distance, and, except for the Δ*φ* = 0 case, their peak values reduce quite significantly if compared with those at *x* = 4*Do*. In the Δ*φ* = 0 case, wall shear stress gets larger as the flow structures traveling to downstream. This may be caused by the closer distance of the hairpin vortex to the wall and the increasing strength of the hairpin vortex during the merging process of the twin SJs. In addition, since the two legs of the dominant vortex is pushed away from each other by its counterpart vortex, the distance of the corresponding peaks becomes much larger than that in the other cases including the single SJ case. This is also consistent with what is observed in Figure 16. In addition, negative values of excess wall shear stress are also found at *x* = 18*Do* in |*z*/*Do*|>2regions, except at Δ*φ* = 3*π*/2. From the instantaneous patterns of excess wall shear stress (Figures 12a, 14a and 15a), these negative values are ascribed to the deficit of wall shear stress outboard of the secondary vortices.

The excess wall shear stress is time-averaged, and then spanwise averaged within a range of *z* = −1.6*Do*~1.6*Do* where most of the non-zero excess stress is covered as revealed in Figure 17, and the streamwise distribution of the newly calculated quantity is then plotted and compared in Figure 18. It is seen that this time- and spatial-averaged excess stress in the single SJ case reaches a gentle peak at about *x* = 3*Do*, and then decreases along the streamwise direction. In all the twin SJ cases, the curve begins with a significant rise, which corresponds to the region of negative excess stress right next to the downstream orifice as shown in Figure 16. Compared to the single SJ, in the Δ*φ* = *π*/2 case the combined vortex structures produce a peak of more than 40% larger at about *x* = 3.5*Do*, after which the mean excess stress decreases with a higher rate, reaching a value of 26% larger at *x* = 22*Do*. The trend is similar in the Δ*φ* = 3*π*/2 case where the two completely separated hairpin vortices are generated. Because of the doubling frequency effect in this case, the mean excess stress becomes almost unchanged after *x* = 12*Do*, which turns to be larger than that in the Δ*φ* = *π*/2 case.

Although some difference is expected between the two cases that generate partially interacting vortex structures, it is still surprising to see a huge difference. The curve in the Δ*φ* = *π* case follows a similar trend to that in the Δ*φ* = *π*/2 case, but with a slightly higher peak. But the curve in the Δ*φ* = 0 case varies quite differently: it reaches a plateau first at about *x* = 4*Do* where the excess stress is only 15% larger than that in the single SJ, and then increases further and reaches a gentle peak at about *x* = 19*Do* where the excess stress is about 50% greater. This huge discrepancy stems from the difference in the evolution of the dominant vortex structure in these two cases. At Δ*φ* = 0, the legs of the dominant hairpin vortex are gradually "pulled" down toward the wall, which can exert stronger influence on the wall shear stress, especially in the far-downstream region.

**Figure 18.** Mean excess wall shear stress along streamwise direction between *z* = −1.6*Do*~1.6*Do*.

#### **5. Conclusions**

The interaction of in-line twin SJs is examined in a laminar boundary layer numerically, that operate at various phase differences along a flat plate. Fluid particles of different colors released from the actuator orifices are tracked so as to show the mixing process of SJs in the boundary layer. In addition, the *Q* criterion is used to identify the twin SJ-induced vortex structures including both the primary and secondary vortices. The impact of these vortex structures on the boundary layer is also investigated through studying the excess wall shear stress. The major findings are as follows:


both the Δ*φ* = 0 and Δ*φ* = *π* cases produce partially interacting vortex structures, the induced two excess stress streaks are significantly different, which stems from the difference in the evolution of the vortex structures, especially of the dominant hairpin vortex, in these two cases.

**Author Contributions:** Writing—original draft preparation, H.W.; writing—review and editing, X.W.; supervision and project administration, D.X. and H.T.; methodology, L.L. and K.Z.; funding acquisition, X.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** The authors gratefully acknowledge the financial support provided by the National Natural Science Foundation of China, China (grant Nos. 12072196) and the Science and Technology Commission of Shanghai Municipality, China (grant No. 19JC1412900), Advanced Jet Propulsion Innovation Center, China (AEAC HKCX2020-02-028), Aeronautical Science Foundation of China (2020Z006057001) and State Key Laboratory of Aerodynamics, China (SKLA-20200303).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclature**


#### **References**


**Zhuoqi Liu 1, Tianyu Pan 1,\*, Shiqi Wang 2,\* and Zhaoqi Yan <sup>3</sup>**


**Abstract:** Active flow control methods are commonly used in expanding the operating range of compressors. Indeed, unsteady active control methods are the main focus of researchers due to their effectiveness. For constructing an unsteady active control system, reliable actuators are significant. To compare with conventional actuators such as synthetic jet actuators and rotating valves, fluidic oscillators have structurally robust characteristics and can generate self-excited and self-sustained oscillating jets, which leads to its higher applicability in compressors under severe working conditions. Thus, to explore the feasibility of unsteady active control systems by the usage of fluidic oscillators, a low-frequency and low-speed oscillator is first designed and experimentally studied for improving the stability of a low-speed axial flow compressor. During the experiments, a special casing is designed to install 15 uniformly distributed oscillators in the tip region of compressor. Based on the unsteady micro injections of the rotor tip with rotor rotation frequency, the results indicate that the frequency/period of oscillators are flexible, in which the values are decoupled with the variation of inlet pressure. When the inlet-to-outlet pressure ratio of the oscillator is in the range of 1.1~2.0, the maximum velocity ranges from 30 m/s to 80 m/s. Moreover, the mass flow rate of the single oscillator only varies from 0.017‰ to 0.059‰ from the designed compressor mass flow rate. For the improvement of the compressor stall margin, the value is 3.45% when the total mass flow of oscillators is 0.08% of the designed compressor mass flow.

**Keywords:** fluidic oscillator; pulsating jets; axial compressor; unsteady stall margin improvement

#### **1. Introduction**

The compressor is one of the key components of aero engines; it is expected to achieve not only high efficiency and a high-pressure ratio but also high stability. There are two ways to improve the stable operating margin, which are passive flow control and active flow control methods. [1] Passive flow control methods can remarkably extend the operating range of compressors. By changing the structure of the casing or blades, they reallocate the flow rate and blade loading in the tip region and extend the operating range of compressors. Typical passive flow control methods include casing treatment [2–4], blade tip winglets [5,6], and vortex generators [7,8]. However, passive control methods will inevitably cause a significant penalty in the efficiency and the pressure ratio [9]. Thus, the concept of active flow control was put forward by Epstein et al. [10] and Ffowcs Williams et al. [11]; by monitoring the flow field parameters and using the actuating mechanisms, signals can be added to the flow field as needed to inhibit the development of pre-stall disturbances. In this way, active flow control methods can not only delay the instability phenomenon but also reduce the efficiency loss.

Active flow control methods, such as tip injection [12,13], variable guide vanes [14,15], and synthetic jets [16,17], can be divided into steady and unsteady active control methods according to the temporal characteristics of the perturbation signals. Steady tip injection

**Citation:** Liu, Z.; Pan, T.; Wang, S.; Yan, Z. Characteristics of a Fluidic Oscillator with Low Frequency and Low Speed and Its Application to Stall Margin Improvement. *Actuators* **2022**, *11*, 341. https://doi.org/ 10.3390/act11120341

Academic Editor: Luigi de Luca

Received: 26 October 2022 Accepted: 21 November 2022 Published: 22 November 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

<sup>1</sup> Research Institute of Aero-Engine, Beihang University, Beijing 100191, China

is a widely studied steady active control method. It has many technical advantages. By changing the jet parameters [18] or using the over-tip recirculation system [19], steady tip injection can act on the tip region and expand the operating range. Based on the total inject momentum, Lin et al. [20] and Li et al. [21] categorized two kinds of stall margin improvement mechanisms. According to the range of spanwise blade load changes caused by air injection from the steady aspect, these two kinds of mechanisms are micro and macro tip injection. Unsteady active flow control methods have higher control efficiency than steady control methods because they can obtain control effects that are several orders of magnitude higher than the percentage of excitation energy of the mainstream [22]. Li et al. [23] used acoustic excitation, Kefalakis et al. [24] used unsteady jets, and Saddoughi et al. [25] used plasma excitation to conduct unsteady active control of compressors, all of which achieved certain stall margin expansion with small energy inputs. However, there has not been a common understanding of the unsteady stall margin expansion mechanism caused by the unsteady nature of excitation.

Due to the harsh working conditions inside compressors, the difficulty of unsteady active flow control lies in the lack of reliable actuators. The wall-attachment oscillator [26] can generate self-excited and self-sustained unsteady jets based on the Coanda effect. It has a simple internal structure and contains no moving parts. Compared with traditional unsteady actuators, such as synthetic jet actuators, rotating valves, or plasma actuators, fluidic oscillators are structurally robust and resistant to electromagnetic interference. In addition, the wall-attachment oscillator has a large jet velocity range and less added mass flow [27]. Thus, it is a good choice to promote unsteady excitation engineering applications, especially when high frequencies are required. In recent years, wall-attachment oscillators have been widely studied in airfoil lift augmentation [28], cavity noise suppression [29], and heat conduction enhancement [30,31]. In the field of aero-engines, wall-attachment oscillators have been applied to adjust the thrust vector [32]; control flow separation at the s-shaped inlet [33], corner of compressor stator [34,35], and outlet guide vane of high-load turbine [36]; and reduce loss caused by turbine tip clearance leakage vortex [37].

Through careful design, wall-attachment oscillators can generate unsteady micro jets with very little additional flow rate; thus, they have the potential to realize effective unsteady active control while further reducing the efficiency penalty. Moreover, wallattachment oscillators are more robust and reliable than traditional unsteady actuators. Therefore, when facing severe working environments inside aeroengines, it is essential to explore the feasibility of unsteady stall margin expansion based on fluidic oscillators. Thus, a low-frequency and low-speed oscillator configuration for a low-speed compressor was redesigned based on the models of Wang et al. [38] and Zhou et al. [39,40]. A series of experiments were conducted to investigate the characteristics of this configuration in detail by using hot-wire anemometer and other measurement methods. A specially designed casing was used to anchor 15 uniformly distributed oscillators in the circumferential direction in the tip region of a low-speed axial compressor rotor. Then, the micro tip injection experiments with the unsteady excitation frequency of rotor rotation were conducted.

This paper is organized as follows. First, a brief introduction to the research subjects and methodology is presented in Section 2. The oscillator characteristics and experimental results in active stall control are presented and analyzed in Section 3, which is followed by a discussion in Section 4. Finally, the main findings are concluded in Section 5.

#### **2. Research Subjects and Methodology**

#### *2.1. Research Subjects*

2.1.1. Low-Speed Axial Compressor Test Rig TA 36

This paper expects unsteady micro tip injection experiments based on fluidic oscillators on the low-speed compressor TA 36, so the test rig will be introduced first.

The schematic of TA 36 and the measurement sections are shown in Figure 1. At Plane A (compressor inlet) and Plane C (50 mm downstream of stator outlet), 4 pressure sensors are embedded in the casing wall at the position shown in Figure 1 to measure the static pressure of the inlet and outlet. In the meantime, 2 pressure rakes are respectively placed in the inlet and outlet ducts to measure total pressure. At Plane B (0.5 chord length upstream of the rotor inlet), 8 high-frequency response pressure sensors are used on TA 36 in order to capture the pressure perturbations near the rotor tip. All of the pressure sensors are equally spaced around the circumference.

**Figure 1.** Schematic of TA 36 and measurement sections.

To design a fluidic oscillator suitable for TA 36, three compressor parameters are highlighted in the present paper, which are the design flow rate, the tip rotational velocity, and the characteristic frequencies associated with unstable phenomena.

Tables 1 and 2 shows the main design parameters of the TA 36 compressor [41]. The design mass flow is 6.5 kg/s. According to the design rotational speed (2930 rpm) and the tip diameter (600 mm), the tip rotational velocity is approximately 92 m/s. Figure 2 shows the static pressure signals of Plane B at 100% design rotational speed. From the change in the propagation velocity of pressure disturbances, it can be judged that TA36 is a spike-type stall. The calculated stall inception frequency and stall cell propagation frequency are about 38 Hz and 27 Hz, respectively. The rotor rotation frequency is 49 Hz since the rotating speed is 2930 rpm (the design speed).

**Table 1.** The aerodynamic design parameters of TA36 [41].



**Table 2.** The structural design parameters of TA36 [41].

**Figure 2.** Time-resolved static pressure disturbance of Plane B at 100% design working speed.

#### 2.1.2. Oscillator Configuration

Combined with the three key parameters of TA36 described above, a low-speed and low-frequency fluidic oscillator was redesigned based on models of Wang et al. [38] and Zhou et al. [39,40] in order to achieve unsteady micro tip injection on TA36.

This oscillator is a wall-attachment pulsed oscillator, and its interior and exterior structure are shown in Figure 3. The corresponding key parameters are listed in Table 3; all of the parameters are nondimensionalized by the depth of the interior structure, which is 2 mm. In Figure 3a, the blue arrows mark the planar core channel, and the orange tubes in Figure 3c are feedback channels. Together these two parts are named as the feedback loop. To meet the requirements of subsequent experiments, the oscillator's response frequency needs to be adjustable. Studies have shown that the response frequency of the wall-attachment pulsed oscillator is mainly controlled by the total length of the feedback loop [38,42]. Therefore, this configuration uses plastic tubes as feedback channels, which makes it easy to adjust its length by changing tubes. The interior structures are processed by the high-precision computerized numerical control method in order to ensure the consistency of the oscillators. In addition, feedback channels connect to the oscillator's main body by threaded pneumatic connectors.

**Table 3.** The key dimensions of the fluidic oscillator.


**Figure 3.** The interior and exterior structure of the fluidic oscillator: (**a**) interior structure; (**b**) critical region; and (**c**) exterior structure.

#### *2.2. Experiment Methods*

#### 2.2.1. Oscillator Test Bench

As shown in Figure 4, the oscillator test bench includes an air pump (JUBA, TJ1600X6- 3201, Wenling, Zhejiang, China) providing high-pressure air, a pressure regulator valve (SMC, AW20-02BG, Tokyo, Japan), a flow meter (Sevenstar, CS100, Tokyo, China), and a pressure gauge (Asmik, MIK-Y290, Hangzhou, Zhejiang, China). The test system is the StreamLine Pro CTA system. The hotwire (Dantec, 55P11, Skovlunde, Danmark) is placed at the center of the left outlet, less than 1 mm away from the outlet section along the depth direction. The hotwire calibration velocity range is 5~100 m/s and the sampling rate is 10 KHz.

**Figure 4.** The oscillator test bench.

The oscillator outlet pressure is equal to the atmospheric pressure *Patm*. The dimensionless parameter inlet-to-outlet pressure ratio is *Pi*/*Patm*, where *Pi* is the total inlet pressure of the oscillator. The feedback channel length is *Lf* and the feedback channel inner diameter is *df* . During the test, *df* is 4 mm and *Lf* ranges from 2.5 m to 7 m. At each length, *Pi*/*Patm* gradually increases from 1.1 to 2.0 in order to obtain thorough frequency and velocity characteristics of this oscillator.

#### 2.2.2. Test Casing

Restricted by the outlet size, pulsating jets from a single oscillator can be quite limited. It is essential to arrange a certain number of oscillators on the casing to achieve the desired control effect. Figure 5 presents the test casing specially designed to fix oscillators. Thirty oscillators are anchored on the test casing, and the number of oscillators can be varied by disconnecting the air supply tubes of some of them. The axial position of oscillators is adjusted by switching the relative position between the inner ring and axial positioning rings. As shown in Figure 5, unsteady jets are 25◦ away from the mainstream and the axial position of actuation is 0.2 cord length downstream of the rotor tip leading edge.

**Figure 5.** The test casing.

Meanwhile, the height of the oscillator is relatively short, so after the oscillator array is installed on the test casing, slot-like chambers still remain in the casing, as marked by the red boxes in Figure 5, in order to protect the rotor.

#### 2.2.3. Oscillator–Compressor Experimental System

The composition of the oscillator–compressor experimental system is shown in Figure 6. It includes an air pump (JUBA, TJ1600X6-3201, China) providing high-pressure air, an air tank (JUBA, JB1.0-8, Wenling, Zhejiang, China) storing the compressed air, a strainer (JUBA, JB-30, Wenling, Zhejiang, China) filtering imported air, a manifold dividing the air 30 ways, a pressure regulator valve (SMC, AW20-02BG, Tokyo, Japan), a flow meter (SMC, PFMB7501, Tokyo, Japan), and a pressure gauge (Asmik, MIK-Y290, Hangzhou, Zhejiang, China). The data acquisition system is a NI-PXI system, which contains 48 channels, and its maximum sampling rate is 100 KHz. TA 36 has a controllable bleed valve to reach the pre-stall region as slowly as possible while recovering from the stall as quickly as possible. This setup ensures that all the steady operating characteristics of TA 36 can be recorded by the standard time–mean instruments.

During the experiment, the throttling process was stopped when the working condition of TA 36 reached the near-stall condition of baseline, the pressure regular valve was opened, and after the oscillators worked stably, we continued the throttling process. Pressure data was continuously collected during the whole process.

**Figure 6.** The oscillator–compressor experimental system.

#### **3. Results**

*3.1. Oscillator Characteristics*

#### 3.1.1. Frequency Response

Hotwire tests were conducted under the condition that *Pi*/*Patm* ranged from 1.1 to 2.0. The feedback channel's inner diameter *df* was 4 mm, and the length *Lf* range was 2.5–7 m.

The oscillator characteristic test results of some feedback channel lengths had very little difference, which is difficult to present. Therefore, in order to show the frequency characteristics that changed with *Lf* as clearly as possible, Figure 7 shows the experimental results with the feedback channel lengths *Lf* of 2.8 m, 3.0 m, 3.2 m, 3.8 m, 4.6 m, 5.2 m, and 6.2 m. With the decrease of *Lf* , the response frequency of the oscillator increased. All frequency–pressure response curves of different *Lf* show similar trends with the increase of *Pi*/*Patm*. When *Pi*/*Patm* increased from 1.1 to 1.25, the response frequency increased quickly, and after *Pi*/*Patm* reached 1.25, the frequency basically did not change but floated around a certain value. According to the research of Cerretelli et al. [43], wall-attachment pulsed fluidic oscillators have the potential to realize frequency–pressure decoupling by designing interior configurations rationally. When *Pi*/*Patm* was 1.25 to 2.0, the oscillator in this paper had response frequencies insensitive to the pressure change, that is, it had the frequency–pressure decoupling characteristics, which can effectively reduce the difficulty of frequency regulation in future experimental work.

**Figure 7.** The oscillator frequency response spectrum.

In addition, three characteristic frequencies related to unstable phenomena are also marked in Figure 7. This configuration can achieve these frequencies by selecting the feedback channel length *Lf* . *Lf* corresponding to the rotating stall frequency (RSF, 27 Hz), the stall inception frequency (SIF, 38 Hz), and the rotor rotating frequency (RRF, 49 Hz) were respectively 5.2 m, 3.8 m, and 3.0 m. For each *Lf* , frequency deviations to the target frequency were less than 1 Hz, which means the error was no more than 3.70% and is in line with the experimental expectations.

When *Pi*/*Patm* is from 1.25 to 2.0, the average frequency on each curve is defined as the average decoupled frequency. Further revealing the frequency characteristics of this oscillator, the average decoupled frequency and the average decoupled period *T* related to *Lf* of all test results are plotted in Figure 8. An empirical formula was obtained by linear fitting, where *c* represents the local speed of sound. Compared with the formulas proposed by Wang et al. [38] and Zhou et al. [39], the *T* of three curves are linear with *Lf* , but the slopes and the intercepts are different. The slope of the formula of this configuration is between Wang's and Zhou's formulas, and the intercept is negative while in the other two formulas, it is zero. These differences indicate that the slope and intercept of an oscillator are sensitive to its internal structure. In future work, if target frequencies are known, the feedback channel length *Lf* can be approximately predicted by this empirical formula.

$$T = 1/f \approx \frac{2.6L\_f}{c} - \left(3.14 \times 10^{-3}\right) \tag{1}$$

**Figure 8.** Changes of decoupled frequency and period *T* with *Lf* [38,39].

#### 3.1.2. Velocity Response

Figure 9 shows the typical velocity response when *Lf* is 3 m and *Pi*/*Patm* is 1.3. In each period *T*, the time of high pulsed speed is *Ton* and the time of low pulsed speed is *Toff* . The duty cycle, which is the ratio of *Ton* and *T*, was about 56%. In addition, the low pulsed speed was close to 0 and barely fluctuated, while the high pulsed frequency had obvious fluctuation. For the convenience of the following description, the average value of the high pulse velocity of each waveform is defined as *Vmax*.

**Figure 9.** *Pi*/*Patm* = 1.3, the typical velocity response of the oscillator.

Figure 10 shows the comparison of velocity responses with *Pi*/*Patm* of 1.3, 1.5, 1.7, and 1.9, respectively. As the inlet pressure rose, the waveform of the velocity response was almost unchanged, and the *Vmax* achieved increased as well as the fluctuating amplitude of the high pulsed velocity. Even if *Pi*/*Patm* was the same, there were still some differences in neighboring waveforms due to velocity fluctuations. Therefore, the maximum velocity response of the oscillator was estimated by the arithmetic average *Vmax* of 10 continuous waveforms, which was 30 m/s to 80 m/s, as shown by the blue curve in Figure 11. The unsteady jet speed generated by this oscillator configuration was lower than the tangential velocity of the rotor tip of TA36, which was larger than 92 m/s, so this configuration meets the requirements of the unsteady micro tip injection used in this paper.

**Figure 10.** The oscillator velocity response at: (**a**) *Pi*/*Patm* = 1.3; (**b**) *Pi*/*Patm* = 1.5; (**c**) *Pi*/*Patm* = 1.7; and (**d**) *Pi*/*Patm* = 1.9.

**Figure 11.** The oscillator's average *Vmax* and mass flow rate response.

#### 3.1.3. Mass Flow Response

The mass flow response of this oscillator is shown by the red curve in Figure 11. Restricted by the throat area, the flow rate increases monotonically with the increase of the inlet pressure. When the range of *Pi*/*Patm* was 1.1 to 2.0, a single oscillator mass flow rate was from 0.115 g/s to 0.388 g/s. The dimensionless parameter *Mn* is the per millage of the oscillator flow rate to the design flow rate of the compressor TA36. In Figure 11, the *Mn* of a single oscillator ranges from 0.017 to 0.059, which is far less than the design flow rate of the compressor and meets the micro-injection requirement of the experiments.

In summary, this new oscillator configuration can realize the three frequencies related to the unstable phenomena of TA 36 by adjusting the feedback channel length. It can also achieve response frequency/period-pressure decoupling as an empirical formula between frequency/period and the feedback channel length was established, which is helpful to predict and select the corresponding target frequency. The *Pi*/*Patm* range was between 1.1 and 2.0, and the maximum jet velocity range was 30 m/s to 80 m/s, which was smaller than the flow velocity at the rotor tip region. The mass flow range of the single oscillator was only 0.017~0.059‰ of the designed compressor flow rate, which ensures that the excitations were micro-injections, and the jet energy basically does not cause efficiency penalty.

#### *3.2. Stall Control Application*

In this study, 15 oscillators distributed evenly above the compressor rotor tip along the circumference were used, with an axial position of 0.2 times the tip cord length downstream of the leading edge. The total inlet pressure of all oscillators was 1.7 times atmospheric pressure, the jet frequency is rotor rotating frequency (RRF), and the maximum actuate velocity is 70 m/s. With a jet flow of only 0.08% of the compressor design mass flow and the maximum jet momentum of 0.27% of the compressor momentum corresponding to design mass flow, the active stall control experiments were conducted on TA 36 to explore the influence on the stability of the low-speed compressor.

#### 3.2.1. Parameter Definition and Experimental Repeatability Verification

In order to quantitatively describe the effect of the unsteady excitation generated by oscillators on compressor stability, the following several parameters are defined.

The total-to-static pressure rise coefficient *ψ* is usually used to describe the low-speed compressor's pressure rise characteristics, which is defined as Formula (2). The flow coefficient *φ* is defined as Formula (3):

$$
\psi = \frac{P\_2 - P\_{t1}}{0.5 \rho l I\_m^2} \tag{2}
$$

$$\phi = \frac{V\_x}{lI\_m} = \frac{q\_v}{lI\_m A\_0} \tag{3}$$

where *P*<sup>2</sup> is the static pressure of the stator exit, *Pt*<sup>1</sup> is the total inlet pressure, *Um* is the tangential velocity of the mid-span, *qv* is volume flow, and *A*<sup>0</sup> is the inlet area.

Efficiency is expressed as Formula (4), *Pt*<sup>2</sup> is the total pressure at stator exit, and *τ* is the torque.

$$\eta = \frac{q\_v (P\_{t2} - P\_{t1})}{\tau} \tag{4}$$

The comprehensive stall margin is adopted in this study to evaluate the compressor stability. The stall margin (SM) and stall margin improvement (SMI) are described as follows:

$$SM = \left(\frac{\psi\_s / \phi\_s}{\psi\_d / \phi\_d} - 1\right) \times 100\% \tag{5}$$

$$SMI = SM\_{exp} - SM\_{ref} \tag{6}$$

where subscripts *s* and *d* refer to stall and design operation points, respectively, and the working conditions corresponding to the subscripts *exp* and *ref* are selected as needed.

In Figure 12, a repetitive result of the baseline characteristics tested at the design speed is shown. The baseline cases are tested with the unmodified TA 36 casing. The results show good consistency, and the subsequent experimental results can be considered to be true and credible.

**Figure 12.** The repeated experiment results of TA36 under baseline condition: (**a**) pressure rise characteristics; and (**b**) efficiency characteristics.

#### 3.2.2. Effect of Oscillators on Compressor Stability

Results of the unsteady actuation using oscillators are shown in Figure 13. In legends, EXP refers to the usage of test casing only, due to the assembly with oscillators; it has 30 slot-like chambers over the rotor tip region compared with the baseline case. EXP\_RRF refers to the usage of test casing and RRF unsteady excitation, and SP refers to stall point.

**Figure 13.** The results of the oscillator–compressor experiments: (**a**) pressure rise characteristics; and (**b**) efficiency characteristics.

In Figure 13a, compared with the baseline (SP 1), test casing (SP 2) itself can obviously improve the stall margin; meanwhile, it also has a significant efficiency loss of 2.77%. This is because of chambers in the rotor tip region after the casing and oscillators are assembled, as shown in Figure 5. The special structure of the test casing is designed by considering the safety and the related parameters of tip injection and casing treatment, which mainly include circumferential coverage and porosity. Studies [13,44] have shown that with the increase of the circumferential coverage of the tip injection, the SMI also increases, and the trend is basically first fast and then slow. The oscillators used in this experiment have two outlets, and the width of each is only 2 mm. When 15 oscillators are applied, the circumferential coverage of air injection is only 1.592%, and considering that the flow rate of air injection is only 0.08% of the compressor's designed flow rate, chambers are left on the casing. In these chambers, the momentum of the injected air is averaged, thereby expanding the circumferential coverage of the excitation to 20 times the original. At the same time, chambers on the casing are similar to the axially slotted casing treatment. According to previous research experience [45], with the increase of the porosity, the SMI and efficiency loss of the slotted casing treatment increase significantly. As mentioned above, 30 oscillators are always installed on the test casing, and the number of oscillators is adjusted by adjusting the air supply. Therefore, from the perspective of casing treatment, the porosity of the casing is independent of the number of oscillators used, which is always 63.67%, and the efficiency loss may be large. Therefore, it is necessary to reasonably design the depth of the oscillator insertion into the casing to reduce the efficiency loss as much as possible. At present, this paper only explores the feasibility of the application of fluidic oscillators on the stall margin control. This form of test casing was selected after considering both safety and its possible stall margin improvement ability. In future work, the influence of unsteady excitation of the oscillator directly acting on the rotor tip will be explored.

Moreover, in Figure 13a, after performing an unsteady excitation with the frequency of RRF, the flow coefficient of SP 3 was slightly smaller than that of SP 2, which was significantly less than the flow coefficient difference between SP 1 and SP 2. This illustrates that unsteady excitation generated by oscillators can have a certain effect. However, it is weaker than the test casing itself. Although the efficiency of case EXP\_RRF is obviously lower than the baseline case, there is barely any efficiency loss compared with the case of EXP, and even some efficiency gain in near-choke and near-stall areas.

In this study, the combination of casing treatment and unsteady tip injection is technically realizing the decoupling of steady and unsteady stall margin improvement effects. When only using the test casing, the stall margin expansion (SP 1 to SP 2) completely

reflected the steady mechanism, and adding the unsteady excitation on this basis, the extra stall margin expansion gain (SP 2 to SP 3) fully reflected the unsteady mechanism. To quantify the steady and unsteady effects in this work, the three stall points in Figure 13a are individually plotted in Figure 14. According to Formula (5), the stall margin of case Baseline was 32.68%, the stall margin of case EXP was 51.65%, and the stall margin of case EXP\_RRF was 55.10%. According to Formula (6), if you select *SMref* as the stall margin of case Baseline and *SMexp* as the stall margin of case EXP, the steady stall margin improvement of the case EXP is 18.97%. Select *SMref* as the stall margin of case EXP and *SMexp* as the stall margin of case EXP\_RRF, and the unsteady stall margin improvement of case EXP\_RRF is 3.45%. Select *SMref* as the stall margin of case Baseline and *SMexp* as the stall margin of case EXP\_RRF, the steady and unsteady coupling improvement is 22.42%.

**Figure 14.** The summary of stall points of TA36.

#### **4. Discussion**

This oscillator configuration has the following advantages when applied to the lowspeed compressor TA 36:


At the same time, because of the special structure of the test casing, the unsteady stall margin expanded by oscillators can be compared with that of the chambers. In Figure 13, the chambers are similar to the conventional passive flow control method and can produce up to 18.97% stall margin expansion, but at a cost of 2.77% efficiency loss. Unsteady micro jets have a weaker effect on stall margin than that of steady methods, but they produce basically no efficiency loss, perhaps because the influence mechanism is changed. Based on previous research experience, [20,21,23], the unsteady micro jets in this study input very small total energy, which can only act on the tip leakage vortex. Therefore, on the one hand, the jets themself suppress the advance and development of the leakage vortex. On the other hand, the unsteady frequency in the jet may interact with some frequencies in the

leakage flow, stimulating the low-energy fluid in the tip region and delaying the leakage vortex generation, thus improving the stable working margin of the compressor.

This paper only verifies the feasibility of the innovative application of the oscillator in the compressor's active stall control field, and there are still many shortcomings. First, the unsteady excitation was not directly acting on the rotor tip flow field, and the stall margin expansion limit was not reached. Second, only the RRF unsteady excitation experiment was carried out, and the influence of other frequencies has not been clarified. Third, the data tested in the flow field were very limited, so it is difficult to fully explore and reveal the unsteady stall margin expansion mechanism. Lastly, compared with synthetic jet actuators or rotating valve actuators, oscillators have more application advantages in generating high excitation frequencies in the harsh operating environment of aeroengines, and their effects on the stability of high-speed compressors have not been explored. On the one hand, future work will optimize the appearance design of the oscillator and explore the influence of other excitation parameters on stability. Numerical simulation will be used to explore the unsteady stall margin mechanism combined with experimental results. On the other hand, oscillators will be applied on a high-speed compressor rig to explore the application value.

#### **5. Conclusions**

In this paper, the characteristics of a new pulsed fluidic oscillator with low frequency and low speed are experimentally studied. Fifteen unsteady oscillators are applied at the rotor tip of a low-speed axial flow compressor (TA36). Based on the usage of low-response steady transducers and high-response dynamic transducers, the influence of the unsteady excitation of oscillators on compressor stability is explored. The concluding remarks can be summarized as follows:


**Author Contributions:** Writing—original draft preparation, Z.L.; writing—review and editing, T.P., S.W., and Z.Y.; supervision and project administration, T.P. and S.W.; methodology, Z.L. and Z.Y.; funding acquisition, T.P. and S.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** The authors acknowledge the support of the National Natural Science Foundation of China (No. 51976005 and No. 52006002), Science Center for Gas Turbine Project (P2022-B-II-004- 001), National Science and Technology Major Project (2017-II-0005-0018), Advanced Jet Propulsion Innovation Center (AEAC HKCX2019-01-016), the Fundamental Research Funds for the Central Universities, and Beijing Nova Program.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland www.mdpi.com

*Actuators* Editorial Office E-mail: actuators@mdpi.com www.mdpi.com/journal/actuators

Disclaimer/Publisher's Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Academic Open Access Publishing

mdpi.com ISBN 978-3-0365-8901-5