## Abstract

Plasma-based flow control was explored as a means of enhancing the performance of a maneuvering flat-plate wing. For this purpose, a numerical investigation was conducted via large-eddy simulation (LES). The wing has a rectangular planform, a thickness to chord ratio of 0.016, and an aspect ratio of 2.0. Computations were carried out at a chord-based Reynolds number of 20,000, such that the configuration and flow conditions are typical of those commonly utilized in a small unmanned air system (UAS). Solutions were obtained to the Navier–Stokes equations, that were augmented by source terms used to represent body forces imparted by plasma actuators on the fluid. A simple phenomenological model provided these body forces resulting from the electric field generated by the plasma. The numerical method is based upon a high-fidelity time-implicit scheme and an implicit LES approach, which were applied to obtain solutions on an overset mesh system. Specific maneuvers considered in the investigation all began at 0 deg angle of attack, and consisted of a pitch-up and return, a pitch-up and hold, and a pitch-up to 60 deg. The maximum angle of attack for the first two maneuvers was 35 deg, which is well above that for static stall. Two different pitch rates were imposed for each of the specified motions. In control situations, a plasma actuator was distributed in the spanwise direction along the wing leading edge, or extended in the chordwise direction along the wing tip. Control solutions were compared with baseline results without actuation in order to assess the benefits of flow control and to determine its effectiveness. In all cases, it was found that plasma control can appreciably improve the time integrated lift over the duration of the maneuvers. The wing-tip actuator could achieve up to a 40% increase in the integrated lift, above that of the baseline value.

## Introduction

A wide range of nontraditional intelligence, surveillance, and reconnaissance missions are currently being proposed for small unmanned air systems. Because of their size and high degree of maneuverability, small UASs are able to perform operations that cannot be carried out by more conventional configurations. At the low speeds associated with these applications, the flight Reynolds number is typically less than 60,000. In order to maximize the lifting surface area for the small geometric size dictated by these UAS planforms, the aspect ratio is generally kept below 2.0. Among a number of different planforms being considered for small UAS operations, are simple rectangular wings with flat cross sections. An overview for the development and design of these vehicles is given by Mueller et al. [1].

Due to their inherent simplicity, various aerodynamic aspects of flat-plate low-Reynolds number rectangular wings have been studied for many years [2–18]. These studies include experimental investigations [2–4,6–9,11–14,16], analytic analysis [5], and numerical simulations [9,10,15,17,18]. While most of the efforts have focused on steady applications, some have also considered plunging motions [13–18]. Because of a loss in aerodynamic efficient at low-Reynolds number, particularly at higher angles of attack, some recent explorations have been devoted to flow control [10–12,19]. It is indicated that active control offers the potential for greatly improving the aerodynamic performance for the common operation of small unmanned air systems.

The use of plasma-based actuators have been found to be an effective means of actively controlling a variety of diverse flow situations. As an alternative to more complex mechanical or pneumatic systems, such devices have no moving parts, low power requirements, a fast time response, are surface adapting, and can operate over a broad range of actuation frequencies. These considerations have spawned experimental explorations for a number of different applications [20–34]. Several related numerical computations have also been performed, where modeling was typically used to represent plasma-induced body forces imparted by actuators on the fluid field [19,35–44].

Dielectric-barrier-discharge (DBD) plasma actuators typically operate in the low radio frequency range (1–10 kHz) with voltage amplitudes of 5–10 kV. An overview of the design, optimization, and application of these actuators has been given by Corke and Post [31]. Experimental measurements indicate that time-averaged body forces generated by such DBD devices are the dominant mechanism for exerting control. Because of this, simple models which account for these forces may be effectively utilized in numerical calculations. This approach is particularly useful in the large-eddy simulation of complex transitional flows, which place severe demands upon computational resources.

The present investigation employs large-eddy simulation to explore the use of plasma-based actuation for aerodynamic efficiency improvement of a rectangular flat-plate wing with an aspect ratio of 2.0 at a chord-based Reynolds number of 20,000, which undergoes unsteady maneuvers that typify small UAS operations. This work compliments the plasma control experiments of Vey et al. [11,12], that were conducted for fixed wings at the same Reynolds number and aspect ratio. For the current simulations, the plate thickness is taken as 1.6% of the chord. The configuration corresponds to wing pitching motion that was studied experimentally by Yilmaz and Rockwell [13,16], and to the computations of Visbal and Garmann [15,17,18] who considered both pitching and perching maneuvers. The wing geometry and flow conditions are also identical to the plasma flow control simulations of Rizzetta and Visbal [19] for the stationary case, and extend those results to the maneuvering situation. In the sections to follow, the governing equations, empirical plasma-force model, numerical method, and LES approach are described. Details of the computations are summarized, and physical characteristics of the resulting simulations are elucidated. All maneuvers began with the wing at 0° angle of attack. The specific cases that were considered include a pitch-up to 35 deg and return (PUAR), a pitch-up to 35 deg and hold (PUAH), and a pitch-up to 60 deg (PU). These motions were selected to be representative of those that might occur during perching or when encountering gusts. Two different pitch rates were specified for each of the three maneuvers. In addition to baseline cases without control, plasma actuators were simulated both along the leading edge and along the wing tip of the configuration. These orientations allowed the plasma force to be directed either in the streamwise direction across the span, or in the outboard spanwise direction across the chord. Results of simulations with actuators are compared with baseline cases in order to assess control effectiveness.

## The Governing Equations

were also employed, and Stokes’ hypothesis for the bulk viscosity coefficient has been invoked.

## The Empirical Plasma Model

Many quantitative aspects of the fundamental processes governing plasma/fluid interactions remain unknown or computationally prohibitive, particularly for transitional and turbulent flows. These circumstances have given rise to the development of a wide spectrum of models with varying degrees of sophistication, that may be employed for more practical simulations. Among the simplified methods focused specifically on discharge/fluid coupling is that of Roth et al. [45,46], who associated transfer of momentum from ions to neutral particles based upon the gradient of electric pressure. A more refined approach, suitable for coupling with fluid response was an empirical model proposed by Shyy et al. [47], using separate estimates for the charge distribution and electric field. Known plasma physics parameters were linked to experimental data. This representation has been successfully employed for several previous simulations of plasma-controlled flows [19,35–40,43,44], and its basic formulation was also adopted in the present investigation.

A schematic representation of a typical single asymmetric DBD plasma actuator is depicted in Fig. 1. The actuator consists of two electrodes that are separated by a dielectric insulator, and mounted on a body surface. An oscillating voltage, in the 1–10 kHz frequency range, is applied across the electrodes, developing an electric field about the actuator. When the imposed voltage is sufficiently high, the dielectric produces a barrier discharge, that weakly ionizes the surrounding gas. Momentum acquired by the resulting charged particles from the electric field, is transferred to the primary neutral molecules by a combination of electrodynamic body forces and poorly understood complex collisional interactions. Because the bulk fluid cannot respond rapidly to the high frequency alternating voltage, the dominant effect of actuation is to impose a time-mean electric field on the external flow. In the numerical simulation of control applications, the entire process may be modeled as a body force vector acting on the net fluid adjacent to the actuator, which produces a flow velocity.

The model for the geometric extent of the plasma field generated by such an actuator is indicated in Fig. 2. The triangular region defined by the line segments OA, OB, and AB constitutes the plasma boundary. Outside of this region the electric field is not considered strong enough to ionize the air [47]. The electric field has its maximum value at point O, and varies linear within OAB. The peak value of the electric field can be estimated from the applied voltage and the spacing between the electrodes. Along the segment AB, the electric field diminishes to its threshold value, which was taken as 30 kV/cm [47]. The electric body force is equal to $qcE$ and provides coupling from the plasma to the fluid, resulting in the source vector $S$ appearing in Eq. (1). Some uncertainty exists regarding the direction of the force vector, which was related to the ratio OA/OB in the original work of Shyy et al. [47]. Within the region OAB, the charge density $qc$ is taken to be constant. The plasma scale parameter $Dc$ arises from nondimensionalization of the governing equations, and represents the ratio of the electrical force of the plasma to the inertial force of the fluid.

Some specific details of the plasma model incorporated in the present simulations, were similar to those of the original experiment of Shyy et al. [47], and to the previous control computations of Rizzetta and Visbal for airfoil sections and wings in both steady flow [19, 44] and for flapping motion [43]. The ratio of the threshold electric field magnitude to its peak value was set to 0.133, with OA (nondimensionalized by the chord) taken as 0.03 and OA/OB = 0.5. For the purposes of the present computations, it is assumed that the plasma force is tangential to the actuator surface. This force is oriented in the streamwise direction across the span for the leading-edge actuator, and in the outboard direction across the chord for the wing-tip actuator. Due to empiricism of the formulation, there is ambiguity regarding the value of the scale parameter $Dc$, which can be increased or decreased to produce more or less force.

DBD actuators are inherently unsteady devices. As mentioned previously, within the context of the empirical model, the body force imposed on the fluid is assumed to be steady owing to the high frequency of the applied voltage. It should be noted that the body force ($qcE$) indicated in Fig. 2 may be directed in a specific direction by proper orientation of the triangle OAB. In the present applications, this force is acting in the streamwise direction for an actuator along the wing leading edge, or in the outboard direction for a wing-tip actuator.

## The Numerical Method

In this expression, which is employed to advance the solution in time, $Qp+1$ is the $p+1$ approximation to $Q$ at the $n+1$ time level $Qn+1$, and $\Delta Q=Qp+1-Qp$. For $p=1$, $Qp=Qn$. Second-order-accurate backward-implicit time differencing was used to obtain temporal derivatives.

The implicit segment of the algorithm (left-hand side of Eq. (12)) incorporates second-order-accurate centered differencing for all spatial derivatives, and utilizes nonlinear artificial dissipation [50] to augment stability. For simplicity, the dissipation terms are not shown in Eq. (12). Efficiency is enhanced by solving this implicit portion of the factorized equations in diagonalized form [51]. Temporal accuracy, which can be degraded by use of the diagonal form, is maintained by utilizing subiterations within a time step. This technique has been commonly invoked in order to reduce errors due to factorization, linearization, diagonalization, and explicit application of boundary conditions. It is useful for achieving temporal accuracy on overset zonal mesh systems, and for a domain decomposition implementation on parallel computing platforms. Any deterioration of the solution caused by use of artificial dissipation and by lower-order spatial resolution of implicit operators is also reduced by the procedure. Three subiterations per time step have been applied in the current simulations to preserve second-order temporal accuracy.

where $Q\u2227$ designates the filtered value of $Q$. It is noted that the filtering operation is a postprocessing technique, applied to the evolving solution in order to regularize features that are captured but poorly resolved. Equation (14) represents a one-parameter family of eighth-order filters, where numerical values for the $an$’s may be found in Ref. [55]. The filter coefficient $\alpha f$ is a free adjustable parameter which may be selected for specific applications, where $|\alpha f|<0.5$. The value of $\alpha f$ determines sharpness of the filter cutoff and has been set to 0.30 for the present simulations.

The aforementioned features of the numerical algorithm are embodied in a parallel version of the time-accurate three-dimensional computer code FDL3DI [55], which has proven to be reliable for steady and unsteady fluid flow problems, including vortex breakdown [56,57], transitional wall jets [58], synthetic jet actuators [59], roughness elements [60], plasma flows [35–40,43,44,61], and direct numerical and large-eddy simulations of subsonic [62,63], and supersonic flow fields [64,65].

## The LES Approach

In the LES approach, physical dissipation at length scales smaller than those in the inertial range is not resolved, thereby allowing for less spatial resolution and a savings in computational resources. For nondissipative numerical schemes, without the use of subgrid-scale (SGS) models, this leads to an accumulation of energy at high mesh wave numbers, and ultimately to numerical instability. Traditionally, explicitly added SGS models are then employed as a means to dissipate this energy. In the present methodology, the effect of the smallest fluid structures is accounted for by a high-fidelity implicit large-eddy simulation (HFILES) technique, which has been successfully utilized for a number of turbulent and transitional computations. The present HFILES approach was first introduced by Visbal et al. [66,67] as a formal alternative to conventional methodologies, and is predicated upon the high-order compact differencing and low-pass spatial filtering schemes, without the inclusion of additional SGS modeling. This technique is similar to monotonically integrated large-eddy simulation (MILES) [68] in that it relies upon the numerical solving procedure to provide the dissipation that is typically supplied by conventional SGS models. Unlike MILES however, dissipation is contributed by the aforementioned high-order Pade-type low-pass filter only at high spatial wavenumbers where the solution is poorly resolved. This provides a mechanism for the turbulence energy to be dissipated at scales that cannot be accurately represented on a given mesh system, in a fashion similar to subgrid modeling. For purely laminar flows, filtering may be required to maintain numerical stability and preclude a transfer of energy to high-frequency spatial modes due to spurious numerical events. The HFILES methodology thereby permits a seamless transition from large-eddy simulation to direct numerical simulation as the resolution is increased. In the HFILES approach, the unfiltered governing equations may be employed, and the computational expense of evaluating subgrid models, which can be substantial, is avoided. This procedure also enables the unified simulation of flow fields where laminar, transitional, and turbulent regions simultaneously coexist.

It should also be noted that the HFILES technique may be interpreted as an approximate deconvolution SGS model [69], which is based upon a truncated series expansion of the inverse filter operator for the unfiltered flow field equations. Mathew et al. [70] have shown that filtering provides a mathematically consistent approximation of unresolved terms arising from any type of nonlinearity. Filtering regularizes the solution, and generates virtual subgrid model terms that are equivalent to those of approximate deconvolution.

## The Unsteady Maneuvers

where $a=11.0$ is a smoothing parameter, $t1=0.5$ is the start-up time, $t2=t1+$ the pitch time, $t3=t2+$ the pause time, $t4=t3+$ the finish-up time.

In the above expressions, $\omega $ is the nondimensional pitch rate, the pitch time = $\alpha max/\omega $, the pause time = 2.0, and the finish-up time = 0.5.

With a pitch rate $\omega =0.05$, maneuvers for a pitch-up and return (PUAR), a pitch-up and hold (PUAH), and a constant pitch-up (PU) are illustrated in Fig. 3. For the PUAR and PUAH maneuvers, $\alpha max$ was 35.0 deg. With the constant pitch-up, $\alpha max$ was set to 70 deg, but the solution was not evolved beyond $\alpha =60.0deg$. This maintained the same start-up description as the other cases. In the PUAH maneuver, $\alpha $ was set to a constant value after the maximum angle of attack was achieved ($\alpha max=35.0deg$), and the function defined in Eq. (15) was abandoned. Corresponding maneuvers for the more rapid pitch rate $\omega =0.20$ are given in Fig. 4. The hold time for the PUAH motions was sufficiently long, so that the duration of the entire maneuver was identical to the PUAR total time (see Figs. 3 and 4).

## Details of the Computations

Details of the computations are provided in following subsections.

### Flow Conditions.

The flow conditions for the computations corresponded to a Mach number $M\u221e=0.1$ and Reynolds number $Re=20,000$, based upon the freestream conditions and the wing chord. These duplicate the flow investigated by Rizzetta and Visbal [19], and are representative of those for small UAS applications.

### Computational Meshes.

Geometry of the flat-plate wing configuration appears in Fig. 5. Because the wing is symmetric about the midspan, only one-half of the configuration is simulated numerically. Moreover, experiments for pitching and heaving motions have indicated flow field symmetry about the midspan, [13,16] so that this representation should provide physically meaningful results. It is perhaps the forced motions of the experiments and the present investigation that preclude the appearance of asymmetries that may evolve in unforced situations due to spurious events. As noted previously, the nondimensional plate thickness $h=0.016$. The computational flow field surrounding the wing is described by an overset mesh system consisting of eight distinct grids, which are shown in Fig. 6. Only a fraction of the mesh points are displayed in the figure. The minimal grid spacing normal to all planes of the plate surface is equal to 0.0001, and the outer domain boundary extended 100 chord lengths upstream, downstream, above, and below the wing configuration. The outboard boundary was located 5.0 chords from the wing tip. On the wing upper surface, the plate was represented by $(331\xd7171)$ mesh points in the streamwise and spanwise directions, respectively. A uniform surface grid spacing of 0.001 was employed in the leading-edge and wing-tip regions, extending from the plate boundaries to beyond the actuator widths. Fewer grid points were distributed over the lower surface of the wing in order to conserve computational resources. In the vertical direction across the thickness of the plate, 61 mesh point were employed. The entire system consisted of approximately $54.9\xd7106$ points, and each of the eight distinct grids included a five-plane overlap in the overset regions. This maintained the high-order differencing and filtering schemes. Because all mesh points from different grids in the overlapping regions were coincident, no interpolation was required.

Domain decomposition was utilized for parallel processing, and was implemented by partitioning each of the eight fundamental grids into smaller subregions. Once again, a five-plane overlap was applied in overlapping subregions, and the entire domain was distributed across 616 computing cores. In conjunction with the aforementioned computational mesh system, an auxiliary system was also developed for the purpose of assessing the effect of spatial resolution on one of the resulting solutions. This auxiliary system maintained the identical regions of clustering and grid stretching ratios as those of its more refined counterpart. The sizes of all meshes on each system are listed in Table 1.

Grid designation | ||
---|---|---|

Mesh number | Fine grid size | Coarse grid size |

1 | ($551\xd7251\xd7231$) | ($276\xd7126\xd7116$) |

2 | ($111\xd767\xd7231$) | ($56\xd738\xd7116$) |

3 | ($111\xd767\xd7231$) | ($56\xd738\xd7116$) |

4 | ($339\xd767\xd761$) | ($174\xd738\xd731$) |

5 | ($119\xd7201\xd7231$) | ($64\xd7101\xd7116$) |

6 | ($120\xd7201\xd7231$) | ($65\xd7101\xd7116$) |

7 | ($322\xd7201\xd771$) | ($157\xd7101\xd741$) |

8 | ($162\xd7201\xd775$) | ($81\xd7101\xd737$) |

Total points | 55,046,049 | 7,196,882 |

Grid designation | ||
---|---|---|

Mesh number | Fine grid size | Coarse grid size |

1 | ($551\xd7251\xd7231$) | ($276\xd7126\xd7116$) |

2 | ($111\xd767\xd7231$) | ($56\xd738\xd7116$) |

3 | ($111\xd767\xd7231$) | ($56\xd738\xd7116$) |

4 | ($339\xd767\xd761$) | ($174\xd738\xd731$) |

5 | ($119\xd7201\xd7231$) | ($64\xd7101\xd7116$) |

6 | ($120\xd7201\xd7231$) | ($65\xd7101\xd7116$) |

7 | ($322\xd7201\xd771$) | ($157\xd7101\xd741$) |

8 | ($162\xd7201\xd775$) | ($81\xd7101\xd737$) |

Total points | 55,046,049 | 7,196,882 |

### Boundary Conditions.

On all planar boundaries of the wing surface, the no slip condition was enforced, along with an isothermal wall, and a third-order accurate implementation of a zero normal pressure gradient. At far field boundaries (upstream, downstream, upper, lower, and outboard), freestream conditions were specified for all dependent variables. A simple extrapolation of variables was applied at the downstream outflow. Grid stretching in far field regions (seen in Fig. 6) transfers information to high spatial wave numbers, and it is then dissipated by the low-pass numerical filter [72]. This technique prevents any spurious reflections, particularly in the outflow area of the computational domain.

### Temporal Considerations.

All computations were performed with a nondimensional time step of $\Delta t=0.0001$. This value was constrained by accuracy considerations, and provided 10,000 times steps for each chord length traveled by any flow structures which convect at the freestream velocity. A similar time step was employed in previous computations [19,43,44].

### Actuator Configurations.

As seen in Fig. 5 are the actuator configurations and plasma force directions considered in the investigation. Along the wing leading edge, the actuator occupies the region defined by $0.025\u2264x\u22640.055$ and $0.0\u2264z\u22640.940$. For the tip region, the actuator is specified within $0.06\u2264x\u22641.0$ and $0.945\u2264z\u22640.975$. It is evident that both actuators have a width of 0.03, and that they do not overlap. These actuator configurations are identical to the ones previously considered by Rizzetta and Visbal [19] for the same wing with steady flow conditions at 8.0 and 25.0 degree angle of attack. With a value of $Dc=1000.0$, effective flow control was achieved in those simulations. The results also compared favorably with the flow-control experiments of Vey et al. [11,12]. The value $Dc=1000.0$ is also adopted in the present simulations. During the investigation of Ref. [19], both continuous and pulsed operation of actuators was explored. Because it was found that pulsing provided less control, only continuously operated actuators are considered here. In the results to follow, solutions for the leading-edge actuator will be referred to with the designation LE, while those for the wing-tip actuator will be termed WT.

### Initial Flow Fields.

The flow fields for all unsteady maneuvers were initialized from solutions with $\alpha =0deg$. Because of the blunt leading and trailing edges of the wing configuration, these results were not steady flows. In each case, the solutions were allowed to evolve sufficiently in time such that an equilibrium state had been attained. These initial flow fields are pictured in Fig. 7 in terms of the $Q$-criterion vortex identification function [73], and contour lines of pressure coefficient at the midspan. The $Q$-criterion has commonly been utilized to represent vortical fluid structures, and has been colored by the streamwise velocity in the figure. Note that the leading-edge actuator (LE) has precluded the formation of small but coherent vortices at the leading edge. With wing-tip actuation (WT), weak tip vortices are generated; in addition to a region of disturbed flow outboard of the wing tip.

## Aerodynamic Force Coefficients for the Maneuvers

Time histories of the aerodynamic force coefficients ($Cl$, $Cd$, $Cm$) are used as metrics to evaluate the performance of the plasma actuators. For each maneuver, results of the LE and WT actuators are compared to the baseline computation without control. The force coefficients also may be integrated in time to provide an additional comparison. Final values of the integrated coefficients represent the resultant forces incurred over the duration of the maneuver. The final values are tabulated in Table 2.

Maneuver | $\omega $ | Configuration | Grid | $Cli(tf)$ | $Cdi(tf)$ | $Cmi(tf)$ |
---|---|---|---|---|---|---|

PUAR | 0.05 | baseline | fine | 25.490 | 13.058 | −3.455 |

PUAR | 0.05 | LE | fine | 26.737 | 14.943 | −7.598 |

PUAR | 0.05 | WT | fine | 32.643 | 16.244 | −4.450 |

PUAH | 0.05 | baseline | fine | 34.241 | 21.453 | −6.290 |

PUAH | 0.05 | LE | fine | 39.232 | 25.702 | −3.209 |

PUAH | 0.05 | WT | fine | 48.186 | 30.348 | −8.756 |

PU | 0.05 | baseline | fine | 26.182 | 21.080 | −5.814 |

PU | 0.05 | LE | fine | 27.517 | 22.831 | −3.369 |

PU | 0.05 | WT | fine | 31.201 | 24.702 | −5.789 |

PUAR | 0.20 | baseline | fine | 10.034 | 5.646 | −1.155 |

PUAR | 0.20 | LE | fine | 10.034 | 6.177 | −3.361 |

PUAR | 0.20 | WT | fine | 12.240 | 6.493 | −1.349 |

PUAH | 0.20 | baseline | fine | 16.176 | 10.308 | −2.582 |

PUAH | 0.20 | LE | fine | 17.253 | 11.158 | −1.373 |

PUAH | 0.20 | WT | fine | 19.355 | 12.369 | −3.029 |

PU | 0.20 | baseline | fine | 8.621 | 6.701 | −1.401 |

PU | 0.20 | baseline | coarse | 8.563 | 6.595 | −1.385 |

PU | 0.20 | LE | fine | 9.377 | 7.564 | −1.959 |

PU | 0.20 | WT | fine | 9.808 | 7.516 | −1.400 |

Maneuver | $\omega $ | Configuration | Grid | $Cli(tf)$ | $Cdi(tf)$ | $Cmi(tf)$ |
---|---|---|---|---|---|---|

PUAR | 0.05 | baseline | fine | 25.490 | 13.058 | −3.455 |

PUAR | 0.05 | LE | fine | 26.737 | 14.943 | −7.598 |

PUAR | 0.05 | WT | fine | 32.643 | 16.244 | −4.450 |

PUAH | 0.05 | baseline | fine | 34.241 | 21.453 | −6.290 |

PUAH | 0.05 | LE | fine | 39.232 | 25.702 | −3.209 |

PUAH | 0.05 | WT | fine | 48.186 | 30.348 | −8.756 |

PU | 0.05 | baseline | fine | 26.182 | 21.080 | −5.814 |

PU | 0.05 | LE | fine | 27.517 | 22.831 | −3.369 |

PU | 0.05 | WT | fine | 31.201 | 24.702 | −5.789 |

PUAR | 0.20 | baseline | fine | 10.034 | 5.646 | −1.155 |

PUAR | 0.20 | LE | fine | 10.034 | 6.177 | −3.361 |

PUAR | 0.20 | WT | fine | 12.240 | 6.493 | −1.349 |

PUAH | 0.20 | baseline | fine | 16.176 | 10.308 | −2.582 |

PUAH | 0.20 | LE | fine | 17.253 | 11.158 | −1.373 |

PUAH | 0.20 | WT | fine | 19.355 | 12.369 | −3.029 |

PU | 0.20 | baseline | fine | 8.621 | 6.701 | −1.401 |

PU | 0.20 | baseline | coarse | 8.563 | 6.595 | −1.385 |

PU | 0.20 | LE | fine | 9.377 | 7.564 | −1.959 |

PU | 0.20 | WT | fine | 9.808 | 7.516 | −1.400 |

### Results for $\omega =0.05$.

Results for all maneuvers with $\omega =0.05$ are found in Figs. 8–13. Displayed in Fig. 8 are time histories of the force coefficients for the PUAR maneuver. It is evident in the figure that the WT actuator appreciably increases the lift. Most of this benefit occurs during the early return portion of the maneuver. As observed in Fig. 9, the integrated lift improves throughout the motion, attaining a 28% improvement over the baseline case by the end of the maneuver with the WT actuator (see Table 2). For the LE actuator, lift enhancement is minimal. It should be noted that use of either actuator results in an increase in drag. This is not unexpected, and although increased drag is typically considered undesirable at cruise conditions, it may be of benefit during perching maneuvers to slow a vehicle while landing. Figures 8 and 9 also indicate that the LE actuator generates a significant amount of nose-up moment. This may adversely affect stability during the course of the maneuver.

Presented in Figs. 10 and 11 are results for the PUAH maneuver. Here, increased lift is seen for both actuators, with the WT actuator being more effective (40.7% for the WT actuator, and 14.6% for the LE actuator). This improvement takes place almost entirely during the hold phase of the maneuver, and is similar to that found previously for the stationary case at $\alpha $ = 25 deg (see Ref. [19]). Behavior of the results with regard to drag and moment is much like that of the PUAR case.

Solutions for the PU maneuver appear in Figs. 12 and 13. Here, the WT actuator has increased the integrated lift by 19.1%, while a more modest improvement occurs for the LE actuator (5.1%). Both actuators produce little additional drag. And while the WT actuator has a negligible effect on the moment compared to the baseline, the LE actuator maintains a constant value of $Cm$ ($Cm\u22480$) for an extended portion of the maneuver.

### Results for $\omega =0.20$.

Figures 14–19 show corresponding results of all maneuvers for the more rapid pitch rate $\omega =0.20$. Time histories of the aerodynamic coefficients for the PUAR case are given in Fig. 14. These results are similar to those for $\omega =0.05$, with no increase of the integrated lift for the LE actuator, but a 22.0% improvement for the WT actuator. The LE actuator also produces a greater variation in $Cm$.

In the PUAH case (Figs. 15 and 16), control is less effective at the higher pitch rate ($Cl$ increase of 19.7% for the WT actuator) as the flow is inertially dominated and viscous effects cannot be easily influenced by plasma forces. Once again, the LE actuator produces a considerable variation in the pitching moment. Despite the inertial dominance at this high pitch rate, for the PU maneuver illustrated in Figs. 19 and 20, there were integrated lift improvements of 13.8% and 8.8%, respectively, for the WT and LE actuators.

### Effect of Grid Resolution.

Because of the large number of cases that that were considered in this investigation, it was not possible to perform a grid resolution study for all maneuvers and pitch rates. However, solutions for the baseline PU maneuver with $\omega =0.20$ were obtained on both the fine and coarse computational mesh systems given in Table 1. Results of the simulations are provided in Figs. 20 and 21. It is evident from the figures that there are only minor differences between the two computations. Although, as mentioned earlier, this case tends to be inertially dominated, the results provide a level of confidence for the adequacy of grid resolution for other cases. It should also be noted that previous solutions for steady cases were obtained on these grids (see Ref. [19]), which also exhibited sufficient spatial resolution.

## Instantaneous Flow Fields for the Maneuvers

Features of the unsteady flow fields are represented in Figs. 22–29. For each value of $\omega $, the pitch-up portions of all the maneuvers are identical up to $\alpha =35.0deg$. Therefore, solutions during the pitch-up for the baseline, and LE and WT actuators can be compared at any instantaneous value of $\alpha $. Visual comparisons are made in terms of planar contours of the streamwise velocity and spanwise vorticity at the midspan symmetry plane, and isosurfaces of the $Q$-criterion.

### Results for $\omega =0.05$.

Planar contours of the streamwise velocity at the symmetry plane are displayed in Fig. 22, at four different values of $\alpha $. At this plane, it is evident that the LE actuator has decreased the amount of low-speed flow and instantaneous separation at all angles of attack. It should be noted however, that the WT actuator is not observed to be effective at the symmetry plane, due to its location along the wing tip. Overall, the WT actuator was found to be quite effective based upon the integrated lift, and this is even somewhat noticeable for $\alpha =35.0deg$.

Corresponding values of instantaneous planar contours of the spanwise vorticity are pictured in Fig. 23. Here the complexity of the flow field, characterized by small-scale structures, is evident. A decrease in the extent of the highly disrupted flow is indicated for the LE actuator. Isosurfaces of the $Q$-criterion, colored by the streamwise velocity, are depicted in Fig. 24. It is indicated that the LE actuator has helped minimize the vertical extent of disrupted flow, as previously noted. Large regions of vortical structures, emanating from the wing tips, are apparent for the WT actuator. The increased lift in the WT case is attributed not only these tip vortices, which are forced outboard of the wing, but also to the general spanwise-directed flow created by the actuator.

The PU maneuver differs from the PUAR and PUAH maneuvers in that it attains a much higher angle of attack. Seen in Fig. 25 are results of the PU maneuver at $\alpha =60.0deg$. It is observed that the flow field disruption is even more extensive at this angle of attack, and the solutions for each case are similar.

### Results for $\omega =0.20$.

Solutions for $\omega =0.20$ are illustrated in Figs. 26–29. All of the instantaneous results shown in Figs. 22–29 at both values of $\omega $, were obtained by periodically recording information during the course of the simulation. Because of the difference in pitch rates, the instantaneous values of $\alpha $ are slightly different in Figs. 26–29, from those in Figs. 22–25. The values are nonetheless approximately equal, so that a comparison between the different pitch rates can be made.

Shown in Fig. 26 are instantaneous contours of the streamwise velocity at the symmetry plane for the pitch-up segment of the maneuvers with $\omega =0.20$. As noted earlier, because viscous forces cannot respond readily to the higher pitch rate, inertia terms of the governing equations tend to dominate the composition of the flow field. Upon comparing Fig. 26 with Fig. 22 at the lower pitch rate, it is quite evident that the flow is less disturbed at the same angles of attack for $\omega =0.20$ than it is for $\omega =0.05$. This angular delay or dynamic-stall lag, is commonly observed for higher pitch rates (e.g., see Ref. [17]), and is brought about by unsteady effects. Apart from the delay, trends with regard to flow control are quite similar for both values of $\omega $. The LE actuator appears to be quite effective, but again it must be remembered that the view in the figure is at the plane of symmetry. Corresponding contours of the spanwise vorticity appear in Fig. 27. The effect of control, the dynamic lagging (compared to Fig. 23), and the presence of fine-scale structures are evident in the figure. Provided in Fig. 28 are the isosurface of the $Q$-criterion for the pitch-up with $\omega =0.20$. Once again, the inertially dominated characteristics of the flow are exhibited upon comparison with Fig. 24.

Displayed in Fig. 29 are the instantaneous results of the PU maneuver at $\alpha =60.0deg$ for $\omega =0.20$. In this situation, due to the extremely high angle of attack, the dynamic-stall lag and angular delay of the flow with $\omega =0.20$ has caught up with that for $\omega =0.05$. Therefore the flow field states for both cases (Fig. 29 and Fig. 25) appear quite similar.

## Features of the Flow Fields and Mechanisms for Control

Due to the sharp plate boundaries, flow separation occurs at the leading edge. As the angle of attack increases during unsteady maneuvers, this results in a massive disruption of the flow just downstream. The situation is then characterized by a chaotic complex flow field, with a large unsteady reverse-flow region. Plasma actuation is able to mitigate effects of separation, or to delay them to higher angles of attack, thereby increasing lift. For the leading-edge actuator, any mitigation is achieved by providing a strong jetlike flow near the wing surface. But because the actuator is located downstream of the leading edge, where separation originates, the flow is difficult to modify by this form of control. In addition, most of the massive disruption takes place away from the surface.

With the wing-tip actuator, the plasma-induced force is able to move some of the disrupted flow outboard of the tip. It also generates large tip vortices. Although these vortices lie somewhat outboard of the wing tip, the circulation about them also facilitates mitigation of the reversed surface flow. It is probably for these reasons that the tip actuator is somewhat more effective than the leading-edge actuator for lift enhancement.

## Summary and Conclusions

Large-eddy simulations were carried out in order to explore the use of plasma-based actuators for lift enhancement during unsteady motions. For these simulations, the subsonic flow past a rectangular wing with an aspect ratio of 2.0 were considered at a flight Reynolds number of 20,000. Three different maneuvers were investigated at two different unsteady rates, such that the configuration, flow conditions, and prescribed motions were representative of those employed during small UAS operations. Solutions were obtained on an overset mesh system utilizing a high-fidelity numerical scheme and an implicit LES approach. A phenomenological model was used to represent plasma forces imparted by actuators upon the fluid field, thus making the simulations tractable.

Two different control strategies were investigated, consisting of actuators that were distributed either spanwise along the wing leading edge or chordwise along the wing tip. Temporal histories of the aerodynamic force coefficients, and integrated values of these coefficients over the duration of the maneuvers, were utilized as metrics by which to judge performance improvements. Both actuators were found to enhance lift during all the maneuvers, but generally, use of the wing-tip actuator led to greater improvement. At the higher of the two pitch rates considered in the study, flow fields were dominated by inertial forces, and plasma control was less effective.

Motions specified in this investigation were defined to approximate those encountered during gust interactions or perching situations when landing. It was found that plasma actuation could increase the integrated lift by over 40%. Use of control also appreciably increased drag, which might be beneficial for perching maneuvers. In addition, the controlling mechanisms considerably influenced the pitching moment. Although the WT actuator generated a greater increase in lift, the LE actuator provided less variation in moment. This might be useful for improving stability characteristics during small UAS applications. Although it is beyond the scope of the present work, simultaneous operation of both a leading-edge and a wing-tip actuator might provide an optimal control strategy.

## Acknowledgment

The work presented here was sponsored by the U. S. Air Force Office of Scientific Research, under a task monitored by D. Smith. Computational resources were supported in part by a grant of supercomputer time from the U. S. Department of Defense Supercomputing Resource Centers at the Stennis Space Center, MS and Wright-Patterson AFB, OH.

### Nomenclature

- $a$ =
smoothing parameter, 11.0

- $c$ =
wing chord

- $Cd$ =
wing drag coefficient

- $Cl$ =
wing lift coefficient

- $Cm$ =
wing pitching moment coefficient about $x=0.25$

- $Dc$ =
plasma scale parameter

- $ec$ =
electron charge, $1.6\xd710-19$coulomb

- $E$ =
nondimensional electric field vector

- $E$ =
total specific energy

- $Er$ =
reference electric field magnitude

- $Ex,Ey,Ez$ =
nondimensional components of the electric field vector

- $F,G,H$ =
inviscid vector fluxes

- $Fv,Gv,Hv$ =
viscous vector fluxes

- $h$ =
nondimensional plate thickness

- $J$ =
Jacobian of the coordinate transformation

- $M$ =
Mach number

- $p$ =
nondimensional static pressure

- $Pr$ =
Prandtl number, 0.73 for air

- $qc$ =
nondimensional charge density

- $Q$ =
vector of dependent variables

- $Qi$ =
components of the heat flux vector

- $Re$ =
reference Reynolds number, $\rho \u221eu\u221ec/\mu \u221e$

- $s$ =
nondimensional wing semispan

- $S$ =
source vector

- $t$ =
nondimensional time

- $tf$ =
value of $t$ at end of maneuver

- $T$ =
nondimensional static temperature

- $u,v,w$ =
nondimensional Cartesian velocity components in the $x,y,z$ directions

- $u1,u2,u3$ =
$u,v,w$

- $U,V,W$ =
contravariant velocity components

- $x,y,z$ =
nondimensional Cartesian coordinates in the streamwise, vertical, and spanwise directions

- $x1,x2,x3$ =
$x,y,z$

- $\alpha $ =
instantaneous of attack

- $\gamma $ =
specific heat ratio, 1.4 for air

- $\delta ij$ =
Kronecker δ function

- $\delta \xi 2,\delta \eta 2,\delta \zeta 2,\u2003$$\delta \xi 6,\delta \eta 6,\delta \zeta 6$ =
second-order and sixth-order finite-difference operators in $\xi ,\eta ,\zeta $

- $\Delta Q$ =
$Qp+1-Qp$

- $\Delta t$ =
time step size

- $\mu $ =
nondimensional molecular viscosity coefficient

- $\xi ,\eta ,\zeta $ =
nondimensional body-fitted computational coordinates

- $\xi t,\xi x,\xi y,\xi z,\u2003$$\eta t,\eta x,\eta y,\eta z,\u2003$$\zeta t,\zeta x,\zeta y,\zeta z$ =
metric coefficients of the coordinate transformation

- $\rho $ =
nondimensional fluid density

- $\rho c$ =
electron charge number density, $1\xd71011$/$cm3$

- $\tau ij$ =
components of the viscous stress tensor

- $\omega $ =
nondimensional pitch rate