SENAC Vol. 55, No. 2（2022. 4）― 13 ―［共同研究成果］

Computational analysis of the body-force models forplasma actuator-assisted flow control simulations
Di Chen*: Dept. Information and Comp. Tech., Tokyo University of Science
Kengo Asada: Dept. Information and Comp. Tech., Tokyo University of Science
Satoshi Sekimoto**: Dept. Information and Comp. Tech., Tokyo University of Science
Kozo Fujii: Dept. Information and Comp. Tech., Tokyo University of Science
Hiroyuki Nishida: Tokyo University of Agriculture and Technology*

Currently, Institute of Industrial Science, University of Tokyo
Currently, Tokyo University of Agriculture and Technology**

1．
． Introduction
The well-known “ion wind” induced by a dielectric barrier discharge plasma actuator (DBDPA) has been extensively used as an active flow control device to avoid flow separation thanks to thehigh availability and feasibility. DBD-PA only consists of two electrodes with a dielectric layer
between them, see Fig. 1(a). Due to the thin and light structure, DBD-PA can be attached on any flator curved surfaces, corner, or edge, where flow separation control1, 2is considered without changing
the original shape of the body surface. A review conducted by Wang et al.3 shows a number of flow
control applications of DBD-PA till the year 2013. The applications include not only on separation
control shown in Fig. 1(b,c)4,5 but also on noise reduction,6 skin-friction reduction,7and else.
DBD-PA for flow separation control is typically installed on the suction surface near the leadingedge1, 2 It can generate a wall-jet flows with the maximum velocity up to around 10m/s, when thepeak-to-peak voltage (𝑉𝑉𝑝𝑝𝑝𝑝 ) and the base frequency of operating alternating current (AC) are 5-20kV
and 1-10kHz, respectively. The ionization effect is largely determined by the applied voltage as wellas electric permittivity and thickness of the dielectric layer. The mechanism of plasma-assisting flow
control is described as the interaction of the ionized gas and the neutral air, which results in anelectrohydrodynamic body-force field coupling with the momentum transfer in the external flow
downstream of the exposed electrode.8-11
Massines et al.12 and Roth et al.,13 proposed the early model of body force field generated by asingle DBD-PA, which was one-dimensional (1-D) based on static formulation and neglected thepresence of the charged particles, therefore it barely fitted two-dimensional (2-D) or threedimensional (3-D) applications. Semi-empirical models of 2-D plasma flow modelling, using linear,14
exponential functions,10, 11 and Gaussian distribution2 of the spatial decay for the 2-D body force
component were reviewed by Corke et al.15 Numerous studies have made great efforts on the bodyforce modelling.2, 12-14

― 14 ―

SENAC Vol. 55, No. （2 2022. 4）

Figure 1. (a) Sketch of DBD-PA. (b) Wind tunnel test on airfoil4. (c) DBD-PA attached on theleading edge of a small model plane (wingspan of 3m) in a flight test5.
The phenomenological such as Suzen-Huang (S-H) model2 and Shyy model14 were largely
employed because of the simplicity in the computational flow simulations of Asada et al.16,17 andVisbal et al.,1 respectively. The well-known analytical S-H model is an electrostatic (semi-empirical)
model derived from Enloe et al.’s work,10 the body force vector can be expressed as,(1)
𝒇𝒇𝑖𝑖 = 𝑄𝑄𝑐𝑐 𝑬𝑬𝑖𝑖 = 𝑄𝑄𝑐𝑐 (−∇𝜙𝜙),

where 𝑬𝑬𝑖𝑖 is the electric field vector. The force is contributed by two different parts: the external
electric potential 𝜙𝜙 and the electric field created by the net charged density 𝑄𝑄𝑐𝑐 . The net body force
obtained by the closed-form solution of the S-H model was well validated in the experimental study,21
by calibrating the nondimensional plasma parameter 𝐷𝐷𝑐𝑐 representing the scaling of the electrical toinertial forces, which also describes the strength and the scale of body force from anotherperspective.16,17, 21,22 These studies indicated that the flow control authority of the plasma actuators
are well realized by the simulations using these models. Airfoil stall characteristics are well simulated,
and the computational and experimental pressure distributions over an airfoil agree pretty well.
Despite the low cost of using the analytical model, it is true that the S-H model neglects thecomplex plasma chemistry which usually causes the highly unsteady forcing on the plasma flow. On
the other hand, the charged-particle models associated with the fully-coupled approaches consider theordinary force diffusion, drift motion, Coulomb acceleration of electrons, and positive and negativeions, respectively. Drift-Diffusion (D-D) model was first developed from the physical modelling,23-25
which focused on the electric-field effects on the charged particles. More recently the D-D model wasused to simulate the discharge plasma evolution and the DBD-PA-induced body-force field.9,26
However, due to the time-consuming computation, few studies applied the D-D model in flow control
simulations. Gaitonde et al.27 conducted plasma-based stall control simulations with coupledapproaches that largely reduced the complexity of the broad-spectrum problem, nonetheless, theinduced flow field of high temporal resolution during a single discharge cycle still remains unclear.
In the present study, the body-force field was obtained by the S-H model, as well as by the D-Dcomputation of high temporal resolution. A high-fidelity Navier-Stokes solver is employed to simulate
the DBD-PA induced flow with the incorporating body force field. Based on the published results by
Chen et al.,28 this report summarizes the nature of the S-H model and the D-D model, clarifies thedifference between the two models from the viewpoint of the DBD-PA flow control authority, andfinally provides several suggestions in using the body-force models for plasma actuator-assisted flow
control simulations.

Computational analysis of the body-force models forplasma actuator-assisted flow control simulations― 15 ―

2．
． Methodology
2.1 Body force modelling
2.1.1 Suzen-Huang model
As we introduced in Eqn. (1), the body force vector is computed by multiplying the charge
density 𝑄𝑄𝑐𝑐 and electric field vector 𝑬𝑬𝑖𝑖 , which are solved in the Maxwell’s equations of the external
electric potential 𝜙𝜙 and the charged particle potential, respectively, as follows∇ ∙ (𝜀𝜀𝑟𝑟 ∇𝜙𝜙) = 0,(2)
𝑄𝑄𝑐𝑐
∇ ∙ (𝜀𝜀𝑟𝑟 ∇𝑄𝑄𝑐𝑐 ) = 2 ,(3)𝜆𝜆𝑑𝑑
where 𝜀𝜀𝑟𝑟 denotes the relative permittivity of the dielectric layer, and 𝜆𝜆𝑑𝑑 denotes the Debye length.
In Fig. 2(a), the boundary conditions for solving Eqn. (2) and (3) on the exposed electrode and thewall above the covered electrode can be written as𝜙𝜙(𝑡𝑡) = 𝜙𝜙 𝑚𝑚𝑚𝑚𝑚𝑚 𝑓𝑓(𝑡𝑡),(4)(5)
𝑄𝑄𝑐𝑐,𝑤𝑤 (𝑥𝑥, 𝑡𝑡) = 𝑄𝑄𝑐𝑐𝑚𝑚𝑚𝑚𝑚𝑚 𝐺𝐺(𝑥𝑥)𝑓𝑓(𝑡𝑡),respectively, where 𝜙𝜙 𝑚𝑚𝑚𝑚𝑚𝑚 and 𝑄𝑄𝑐𝑐𝑚𝑚𝑚𝑚𝑚𝑚 are the maximum values of the external electric potential andthe charge density, respectively. Consequently, the time variation of the non-dimensional S-H bodyforce is given as follows(6)
𝑺𝑺𝑖𝑖 (𝑥𝑥, 𝑦𝑦, 𝑧𝑧, 𝑡𝑡) = 𝐷𝐷𝑐𝑐 sin2 (2𝜋𝜋𝐹𝐹𝑏𝑏𝑏𝑏𝑏𝑏𝑏𝑏 𝑡𝑡) 𝑺𝑺𝑆𝑆𝑆𝑆 (𝑥𝑥, 𝑦𝑦, 𝑧𝑧),
where 𝐹𝐹𝑏𝑏𝑏𝑏𝑏𝑏𝑏𝑏 is the nondimensional base frequency of the AC power, 𝑺𝑺𝑆𝑆𝑆𝑆 is the non-dimensional
body-force vector of S-H model. 𝑺𝑺𝑆𝑆𝑆𝑆 follows a half Gaussian function 𝐺𝐺(𝑥𝑥) in two-dimensional
spatial distribution described in Fig. 2(b).2 The magnitude of the S-H body force is manipulated by
the nondimensional plasma parameter 𝐷𝐷𝑐𝑐 , which is given empirically as
𝑄𝑄𝑐𝑐,𝑟𝑟𝑟𝑟𝑟𝑟 𝐸𝐸𝑟𝑟𝑟𝑟𝑟𝑟 𝐿𝐿𝑟𝑟𝑟𝑟𝑟𝑟
𝐷𝐷𝑐𝑐 =.
(7)𝜌𝜌𝑟𝑟𝑟𝑟𝑟𝑟 𝑈𝑈𝑟𝑟𝑟𝑟𝑟𝑟 2
However, it needs further calibration by experiments and numerical simulations according to the
induced flow field.21
2.1.2 Drift-Diffusion modelNotwithstanding the capability of the D-D model in the 3-D body force simulation by Nishida
et al.,26 the model comparison is conducted in 2-D simulations for convenience and simplicity in thisstudy. The electron, the positive ion, and the negative ion are considered with basic plasma chemistry
including electron impact ionization, attachment, and recombination. The ionization and attachment
coefficients and electron mobility are calculated by BOLSIG+29 a free electron Boltzmann equation
solver, assuming the ambient gas is air (N2:O2 = 0.8:0.2).
Here the governing equations to be solved are explained in detail in the previous studies,9, 24-26,
28, 30the computational electrohydrodynamic (EHD) force regarded as the two-dimensional bodyforce is expressed as the rate of momentum transfer per unit volume due to collisions,
𝐷𝐷𝑝𝑝
𝐷𝐷𝑒𝑒
𝐷𝐷𝑛𝑛
𝒇𝒇 = 𝑒𝑒𝑐𝑐 (𝑛𝑛𝑝𝑝 − 𝑛𝑛𝑒𝑒 − 𝑛𝑛𝑛𝑛 )𝑬𝑬 − [ ∇𝑛𝑛𝑝𝑝 + ∇𝑛𝑛𝑒𝑒 +∇𝑛𝑛 ],
(8)
𝜇𝜇𝑝𝑝
𝜇𝜇𝑒𝑒
𝜇𝜇𝑛𝑛 𝑛𝑛where 𝑛𝑛 is the plasma density, 𝜇𝜇 is the charged particle mobility, and their subscript e, p, n denotethe electron, the positive ion, and the negative ion, D is the coefficients of diffusion. The first term on
the right side of Eqn. (8) plays a dominant role and corresponds to the expression used in S-H modelin Eqn. (1), while the second term represents the diffusion effect. The obtained 2-D body force in Eqn.
(8) is first averaged in phases, then input into the N-S solver for the 2-D simulation of the induced
flow, however the spanwise uniformity of body-force distribution is assumed in the 3-D simulation
for simplicity.

― 16 ―

SENAC Vol. 55, No. （2 2022. 4）

Figure 2. (a) Boundary condition of charged particle in S-H model, following a half
Gaussian distribution.2 (b) Spatial and temporal distribution of body force in S-Hmodel.
2.2 Computational fluid dynamics (CFD)
Two typical flows are considered. First, the flow field is assumed to be globally quiescent andlaminar over a flat plate, then the strength of the induced flow is investigated. Second, the body-force
models are applied and examined in practical flow control simulations, where the separated flow overa NACA0015 airfoil with a light-stall angle of attack of 12° in turbulence transition state (Reynolds
number is 63,000) is investigated. The Mach number is 0.2, the specific heat ratio (𝛾𝛾) is 1.4, and the
Prandtl number (Pr) is 0.72, those are keeping the same as the previous simulation16 and theexperiment setup.21 On the surface of the flat plate and the airfoil, no-slip and adiabatic conditions
are imposed.
The flow field is described by the 3-D compressible Navier-Stokes equations as below,
augmented by the term 𝑺𝑺𝑖𝑖 representing the local forcing on the ionized region by DBD-PA. The nondimensional forms of the continuity, momentum, energy equations, as well as ideal gas equation are
written as follows:𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕𝒖𝒖𝑘𝑘
+
= 0,
(9)
𝜕𝜕𝜕𝜕
𝜕𝜕𝑥𝑥𝑘𝑘𝜕𝜕𝜕𝜕𝒖𝒖𝑖𝑖 𝜕𝜕(𝜌𝜌𝒖𝒖𝑖𝑖 𝒖𝒖𝑘𝑘 + 𝑝𝑝𝜹𝜹𝑖𝑖𝑖𝑖 )1 𝜕𝜕𝝉𝝉𝑖𝑖𝑖𝑖
+=+ 𝑺𝑺𝑖𝑖 ,(10)
𝜕𝜕𝒙𝒙𝑘𝑘𝑅𝑅𝑅𝑅 𝜕𝜕𝒙𝒙𝑘𝑘
𝜕𝜕𝜕𝜕
1 𝜕𝜕𝒖𝒖𝑙𝑙 𝝉𝝉𝑘𝑘𝑘𝑘
1
𝜕𝜕𝒒𝒒𝑘𝑘𝜕𝜕𝜕𝜕 𝜕𝜕((𝑒𝑒 + 𝑝𝑝)𝒖𝒖𝑘𝑘 )+=(−) + 𝑺𝑺𝑘𝑘 𝒖𝒖𝑘𝑘
(11)2(𝛾𝛾 − 1)𝑃𝑃𝑃𝑃𝑀𝑀∞ 𝜕𝜕𝒙𝒙𝑘𝑘
𝑅𝑅𝑅𝑅 𝜕𝜕𝒙𝒙𝑘𝑘
𝜕𝜕𝒙𝒙𝑘𝑘𝜕𝜕𝜕𝜕
𝑝𝑝 = (𝛾𝛾 − 1)(𝑒𝑒 − 0.5𝜌𝜌𝒖𝒖𝑘𝑘 𝒖𝒖𝑘𝑘 )(12)
where 𝒙𝒙𝑖𝑖 is the position vector, 𝒖𝒖𝑖𝑖 is the velocity vector, 𝒒𝒒𝑘𝑘 is the heat flux vector, 𝜌𝜌 is thedensity, 𝑝𝑝 is the static pressure, 𝑒𝑒 is the total energy per unit volume, 𝝉𝝉𝑖𝑖𝑖𝑖 is the viscous stresstensor, 𝜹𝜹𝑖𝑖𝑖𝑖 is the Kronecker delta, 𝑺𝑺𝑖𝑖 is the body force vector, 𝛾𝛾 is the ratio of specific heats, 𝑡𝑡 is
the time. In S-H model, 𝑺𝑺𝑖𝑖 is directly given by Eqn. (6) with a proper 𝐷𝐷𝑐𝑐 value, while 𝒇𝒇𝑖𝑖 in D-D
model (Eqn. (8)) is normalized as
1
𝒇𝒇 .
(13)
𝑺𝑺𝑖𝑖 =2 𝑖𝑖
𝜌𝜌∞ 𝑈𝑈∞
2.3 Numerical schemes
To solve the governing equations in Eqn. (9) – (12), this research employs a compressible NavierStokes flow solver LANS3D,31,32 which is extensively used and well validated in recent 30 years. The
implicit large-eddy simulations (LES) are conducted for solving the more complex flow field above
the airfoil in turbulent transition state.16 All the spatial derivatives are obtained with a sixth-order

Computational analysis of the body-force models forplasma actuator-assisted flow control simulations― 17 ―

compact difference scheme.33 Lower-upper symmetric alternating direction implicit and symmetric
Gauss-Seidel (ADI-SGS) method is utilized for time integration. 10th order filtering33 is applied witha filtering coefficient of 0.42.
The current numerical schemes in D-D force computation largely follow Nishida et al.’s previous
studies.9, 26, 30 given by the formula in Eqn. (6). The input time interval of the S-H and the D-D force
match the time step of CFD simulation in real scale.
2.4 Computational grids
The information of the computational grid systems for the body force computation and the CFD
are listed in Table 1. In the 2-D simulation of the quiescent air, the 2-D computational domains andgrids are shown in Fig. 3(a), the D-D force (red zone) is computed and interpolated into the CFD grids(white background). In the 3-D simulations of the separated flow above the airfoil, the C-type grid isdesigned, the near and far field are shown in Fig. 3(b) and (c), respectively. The spanwise length ofthe computational region is set to 0.2 times the chord length. The D-D force is interpolated into theactuator zone (highlighted in green). A message passing interface (MPI) between the fluid zones is
employed for the parallel computing in CFD simulations.

Figure 3. Domain incorporation of the D-D body force (in red) into CFD (in white) for the flat
plate (a), and the NACA0015 airfoil (b), the green grids denote the overset actuator zone, theexposed electrode is highlighted in the zoom-in snip; (c) Computational domain of the airfoil andthe schematic of inflow, 𝐿𝐿𝑟𝑟𝑟𝑟𝑟𝑟 = 𝑐𝑐 = 0.1m.

Table 1. Grid density of the body-force computation and CFD
Computation domain𝜉𝜉𝜂𝜂𝜁𝜁D-D body force (red in Fig. 3a, b)
600
3
250
Flat plate (white in Fig. 3a)
873
3
416
Airfoil surface (white in Fig. 3b)
759
134
179

Total points
450,000
1,089,504
18,205,374

― 18 ―

SENAC Vol. 55, No. （2 2022. 4）

Airfoil actuator (green in Fig. 3b)
149
134
111
2,216,226* In curvilinear coordinates (𝜉𝜉, 𝜂𝜂, 𝜁𝜁), 𝜉𝜉 and 𝜁𝜁 denote wall-parallel and wall-normal directions,
respectively, 𝜂𝜂 denotes the spanwise direction in the airfoil case.

2.5 DBD-PA setup
Fig. 4(a) shows the configuration of DBD-PA in the D-D computational domain, which isconsistent with Sekimoto et al.’s experiment.34 The physical parameters as shown in Fig. 4(a) areprecalculated in the previous numerical studies using D-D model.26, 29 The plasma actuator consistsof a 160𝜇𝜇𝜇𝜇 thick Kapton dielectric layer and two 50𝜇𝜇𝜇𝜇 thick copper electrodes. The thickness ofthe exposed electrode can be neglected in CFD simulations. A sinusoidal form 0.5 𝑉𝑉𝑝𝑝𝑝𝑝 sin (2𝜋𝜋𝑓𝑓𝑏𝑏𝑏𝑏𝑏𝑏𝑏𝑏 𝑡𝑡)
of AC power is applied on the electrodes, the peak-to-peak voltage (𝑉𝑉𝑝𝑝𝑝𝑝 ) amplitude is from 7kV to
20kV and the base frequency (𝑓𝑓𝑏𝑏𝑏𝑏𝑏𝑏𝑏𝑏 ) is set to 10kHz. In the quiescent flow simulation, the continuousmode is applied, therefore the flow is permanently driven by the DBD-PA.
In the separated flow over the airfoil, the exposed electrode of DBD-PA is placed at 5% chord
length from the leading edge, the DBD-PA is periodically activated in the duty cycle as shown in Fig.
4(b). It is called burst mode actuation compared to the continuous mode actuation generally used. The
burst frequency 𝑓𝑓 + = 500Hz and the burst ratio BR = 0.1. The burst mode with a proper 𝐹𝐹 + is wellknown to have better performance in flow separation control than the continuous mode.15,16, 19, 22, 34
𝐹𝐹 + and 𝐹𝐹𝑏𝑏𝑏𝑏𝑏𝑏𝑏𝑏 are the normalized values of burst and base frequency, respectively.

Figure 4. (a) PA configuration and physical parameters in the D-D computation, including
electrodes and dielectric layer highlighted in yellow and grey, respectively (𝐿𝐿𝑟𝑟𝑟𝑟𝑟𝑟 = 0.1𝑚𝑚);(b) Schematic diagram of a bursting wave in burst actuation mode.
2.6 Computational resources
To obtain the time-varying D-D force fields which are necessary as an input data for the flow
field simulations, an in-house computation server is used. It takes around 90 hours using one node
per case, which computes 4.5 discharge cycles. For the CFD simulations, the supercomputer
subsystem named AOBA-A of the Cyberscience Center in Tohoku University is used. 8 nodes areused for computing the flow fields in the divided subdomains with MPI techniques, it takes around
26 hours per case (non-dimensional time 𝑡𝑡𝑈𝑈𝑟𝑟𝑟𝑟𝑟𝑟 /𝐿𝐿𝑟𝑟𝑟𝑟𝑟𝑟 = 80) of the flat-plate flow simulation, andaround 70 hours per case (non-dimensional time 𝑡𝑡𝑈𝑈∞ /𝑐𝑐 = 9.6) of the airfoil flow simulation.
． Results
3．

3.1 Body force fields of DBD-PA
The body force terms on the right side of Eqn. (10) and (11) are computed by the D-D modeland the S-H model. Fig. 5 show the time history of the wall-parallel EHD net force within 3 cycles,Computational analysis of the body-force models forplasma actuator-assisted flow control simulations― 19 ―

the D-D force computations are conducted at the 𝑉𝑉𝑝𝑝𝑝𝑝 from 7kV to 20kV. In the case of S-H model,
the nondimensional plasma parameter 𝐷𝐷𝑐𝑐 is set to 0.0117 referred to the previous parametricstudies,21 in which the near-field induced flow agrees well with the experiment of 𝑉𝑉𝑝𝑝𝑝𝑝 = 7𝑘𝑘𝑘𝑘. Dueto the lack of strong evidence relating the controlling parameter 𝐷𝐷𝑐𝑐 to the applied voltage, we only
plot the results of the S-H model of 𝐷𝐷𝑐𝑐 = 0.0117 (SH-Dc0.0117) at the corresponding voltage of
7kV for all the following comparisons.
As shown in Fig. 5 the time variation of body force in S-H model which was developed by Asadaet al.16, 17 assumes to follow the function of sin(2𝜋𝜋𝐹𝐹𝑏𝑏𝑏𝑏𝑏𝑏𝑏𝑏 𝑡𝑡) in Eqn. (14). The net force is positive all
the time and symmetrical between negative- and positive-going phases, that probably overestimatesthe input power when the body force peak is approximately same with that in the D-D case of 7kV.
The time-averaged value of the net body force in Fig. 5(c) are 10.21mN for D-D model at 7kV,
18.45mN for 10kV, 21.07mN for S-H model with 𝐷𝐷𝑐𝑐 = 0.0117. The power law of 𝑉𝑉𝑝𝑝𝑝𝑝1.49 agrees
with Nishida et al.’s 2-D simulation results26 in a range of 𝑉𝑉𝑝𝑝𝑝𝑝 from 10kV to 24kV when the relative
permittivity of the dielectric layer (glass epoxy) is set to 5.0.

Figure 5. Time history (three periods) of the dimensional EHD net force computed in the entire
body-force domain using the D-D model at 7kV, 8kV, 9kV, 10kV (a), and 12kV, 15kV, 20kV (b),as well as the S-H model with 𝐷𝐷𝑐𝑐 = 0.0117; (c) Time-averaged D-D body force .vs. peak-topeak voltage and S-H body force.

3.2 Induced flow in quiescent field

When the body force fields are obtained by D-D or S-H model, it is quite straightforward to
relate the force production to the induced flow structure in quiescent air over a flat plate. The timevarying body force terms are incorporated into the right-hand side of the N-S equations, the induced
flow fields are then computed by the CFD solver. The time-averaged wall-parallel velocity contours
using the D-D model from 7kV to 9kV are compared with the corresponding experimental results,34

― 20 ―

SENAC Vol. 55, No. （2 2022. 4）

as well as the S-H model with 𝐷𝐷𝑐𝑐 = 0.0117 in Fig. 6(a), which is determined based on the previous
study, where Aono et al. conducted showed the guideline of Dc values by the comparison ofexperimental and computational result with S-H model.21 Fig. 6(b) shows the maximum velocity
measured at 20mm downstream (position A), the far-field numerical and experimental results arecompared. The current results are obtained in the flow field averaged from 0.4s to 0.8s.

Figure 6. (a) Wall-parallel velocity fields in the present CFD simulations (left) and our previous
experiments34 (right); (b) Far-field velocity (𝑢𝑢𝐴𝐴 ) at 20mm downstream.

The induced flow structures and the downstream velocity (𝑢𝑢𝐴𝐴 ) computed by the D-D model show
very good agreement with the experiment. Compared to 8kV, the velocity decline at 9kV is probably
caused by the less body force production in Fig. 5(c) due to the unstable EHD force pushing fluid
towards downstream in positive-going phase. However, the induced flow of S-H model is muchstronger than those of D-D models from 7kV to 9kV, the local velocity at 20mm downstream is evenhigher than the D-D case of 12kV in Fig. 6(b).
Fig. 7 further shows the near-field velocity profiles in a range of 𝑥𝑥/𝐿𝐿𝑟𝑟𝑟𝑟𝑟𝑟 = 0.02 to 0.08 (2mmto 8mm) in the downstream of the exposed electrode, as shown in the top of Fig. 7(b). The experimentof 9kV shows the remarkable unsteadiness, which is the potential cause of the reversed induced flowat 𝑥𝑥/𝐿𝐿𝑟𝑟𝑟𝑟𝑟𝑟 = 0.02, and the exceeding velocity compared to the D-D result at the other downstream
positions. On the other hand, the near-field velocity profiles of the S-H model show the much smaller
gradient of the velocity field, compared to the D-D results.

Computational analysis of the body-force models forplasma actuator-assisted flow control simulations― 21 ―

Figure 7. (a) Wall-parallel velocity profiles in the near field; (b) Thrust at 6mm downstream.
3.3 Flow separation control over NACA0015 airfoil
The PA-assisted flow control have gained great interest in a variety of separated flow fields overthe Reynolds number from 103 to 106 ,36and the angle of attack of pre-to-post stall from 10 degrees
to 16 degrees. In this report, a comparative study among the D-D model, S-H model and experimentis conducted in the separated flow over a NACA0015 airfoil at a light-stall angle of attack of 12°
and a transitional Reynolds number of 63,000. Burst actuation (𝐹𝐹 + = 5, BR = 0.1) is set to maximize
the separation control effect.35, 37 In the following results, “EXP-7kV” denotes the experimental resultwith the same setups of flow filed and DBD-PA in Asada et al.’s work,37 “NOACT” is the case without
control.
Fig. 8 shows that the flows are almost completely attached in all the three controlled cases, italso shows little difference between the S-H and the D-D cases either in velocity field or in TKE(turbulent kinetic energy) field. The TKE has overshoot in the transition position of the shear layer,
and the overshoot becomes larger when the transitional shear layer is thick.36 The chordwise locationsof 𝑇𝑇𝑇𝑇𝑇𝑇𝑚𝑚𝑚𝑚𝑚𝑚 peaks in Fig. 8(c) show very good agreement with the reattachment positions in Fig.
9(b). However, the S-H case shows the relatively larger TKE overshoot compared to the D-D cases,
due to the delay of turbulent transition. On the contrary, the early and smooth transition process in theD-D cases results in the smaller TKE overshoots.

― 22 ―

SENAC Vol. 55, No. （2 2022. 4）

Figure 8. Distribution of chordwise velocity (a), turbulence kinetic energy (b); Spatial
development of maximum turbulent kinetic energy (c) obtained from (b).
Distributions of the aerodynamic coefficients in the chordwise direction are further shown in Fig.
9. In the cases with PA control, the time- and span-averaged pressure coefficient (𝐶𝐶𝑝𝑝 ) shows the
suction peaks as well as the plateau distributions which corresponds to the so-called laminar
separation bubble. The 𝐶𝐶𝑝𝑝 curves indicate the lift increase in the cases with PA control, while theskin-friction coefficient (𝐶𝐶𝑓𝑓 ) shows the reduction of pressure drag, as well as the flow-separated
regions shown in the bottom of Fig. 9(b). Better agreement of 𝐶𝐶𝑝𝑝 distribution can be observed
between the D-D cases and the experiment at 7kV, compared with the S-H case. As to the 𝐶𝐶𝑓𝑓 curves
in Fig. 9(b), it is found that the reattachment points in the D-D cases move closer towards the leadingedge compared to the S-H case. The discrepancy of the reattachment position is partly caused by thedifferent location of the maximum induced velocity in the quiescent flow simulation, which is
discussed in detail by Chen et al.28

Computational analysis of the body-force models forplasma actuator-assisted flow control simulations― 23 ―

Figure 9. Pressure coefficient (a), skin-friction coefficient (b) along the airfoil chord.
4．
． Summary and suggestions
In computational studies of PA-assisted flow control, S-H model can be quickly implemented
because of the low cost, it has been widely used and also well validated by experiments. On the other
hand, D-D model as a high-fidelity body force model has never been employed before in practicalflow control simulations. The study shown in Ref. 28 and summarized in this report by the present
authors first shows the usage of the two models, as well as the result comparison in the quiescent flow
and airfoil flow simulation. Using the different body force model, the effect of separation control isalmost the same in the separated airfoil flow, despite the large discrepancy of the induced flow in thequiescent air.
When we consider the cost reduction of the body force computation, for the current flow
condition in light separated and turbulent transition state, S-H model is acceptable within thereasonable range of computational accuracy so far as discussions are qualitative. However, it isprobably necessary to use the D-D model of high temporal fidelity in fully turbulent flows, as the
body-force unsteadiness within the duty cycle may not be neglected.
Acknowledge
The flow field computations presented in this report is performed on the “SX-Aurora TSUBASA”
in the Cyberscience Center, Tohoku University. Supports by the center staff are very muchacknowledged. This study is also partly supported by the Grants-in-Aid for Scientific Research A(JSPS 18H03816).

References
1
2
3
4
5
6
7

M. Visbal and D. Gaitonde, AIAA paper 2006-3230 (2006).
Y. Suzen and G. Huang, AIAA paper 2005-4633 (2005).
J.-J. Wang, K.-S. Choi, L.-H. Feng, T.N. Jukes, and R.D. Whalley, Prog. Aerospace Sci. 62, 52(2013).
S. Sekimoto, T. Nonomura, and K. Fujii, AIAA J. 55, 1385 (2017).
S. Sekimoto, K. Fujii, H. Yoneda, AIAA paper 2020-0825 (2019).
E. Peers, X. Huang, and X. Luo, IEEE Trans. Plasma Sci. 37, 2250 (2009).
S. Roy, P. Zhao, A. DasGupta, and J. Soni, AIP Adv. 6, (2016).

― 24 ―

8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

32
33
34
35
36
37

SENAC Vol. 55, No. （2 2022. 4）

J.P. Boeuf and L.C. Pitchford, J. Appl. Phys. 97, 103307 (2005).
H. Nishida and T. Abe, J. Appl. Phys. 110, 013302 (2011).
C.L. Enloe, T.E. McLaughlin, R.D. VanDyken, K.D. Kachner, E.J. Jumper, T.C. Corke, M. Post,
and O. Haddad, AIAA J. 42, 595 (2004).
D.M. Orlov, Modelling and Simulation of Single Dielectric Barrier Discharge Plasma
Actuators, Thesis (Ph.D.), University of Notre Dame (2006).
F. Massines, A. Rabehi, P. Decomps, R. Ben Gadri, P. Ségur, and C. Mayoux, J. Appl. Phys. 83,
2950 (1998).
J.R. Roth, D.M. Sherman, and S.P. Wilkinson, AIAA J. 38, 1166 (2000).
W. Shyy, B. Jayaraman, and A. Andersson, J. Appl. Phys. 92, 6434 (2002).
T.C. Corke, C.L. Enloe, and S.P. Wilkinson, Annal Rev. Fluid Mech. 42, 505 (2009).
K. Asada and K. Fujii, in AIAA 5th Flow Control Conf. Reston, Virginia (2010), pp. 1–12.
K. Asada, T. Nonomura, H. Aono, M. Sato, K. Okada, and K. Fujii, Int. J. Comp. Fluid Dyn. 29,
215 (2015).
S. Kawai, T. Bouwhuis, Y. Abe, A. Yakeno, T. Nonomura, H. Aono, A. Oyama, H.W.M.M.
Hoeijmakers, and K. Fujii, Theor. Comput. Fluid Dyn. 32, 805 (2018).
S. Sekimoto, T. Nonomura, and K. Fujii, AIAA J. 55, 1385 (2017).
M. Kotsonis, S. Ghaemi, L. Veldhuis, and F. Scarano, J. Phys. D. Appl. Phys. 44, 045204 (2011).
H. Aono, S. Sekimoto, M. Sato, A. Yakeno, T. Nonomura, and K. Fujii, Mech. Eng. J. 2, 15 (2015).
K. Fujii, Appl. Sci. 8, 546 (2018).
V.I. Gibalov and G.J. Pietsch, J. Phys. D. Appl. Phys. 37, 2082 (2004).
J.P. Boeuf and L.C. Pitchford, J. Appl. Phys. 97, 103307 (2005).
A. V Likhanskii, M.N. Shneider, S.O. Macheret, R.B. Miles, A. V Likhanskii, M.N. Shneider, S.O.
Macheret, and R.B. Miles, J. Appl. Phys. 103, 13 (2008).
H. Nishida, T. Nonomura, and T. Abe, AIAA J. 54, 659 (2016).
D.V. Gaitonde, M.R. Visbal, and S. Rouy, in Proc. ASME/FEDSM2006 ( 2006), pp. 1339–1356.
D. Chen, K. Asada, S. Sekimoto, K. Fujii, and H. Nishida, Phys. Fluids 33, 037115 (2021).
URL: http://www.bolsig.laplace.univ-tlse.fr/ for more information about the BOLSIG+ database.
K. Nakai, A. Nakano, and H. Nishida, AIAA J. 1 (2021).
K. Fujii, H. Endo and M. Yasuhara, Activities of computational fluid dynamics in Japan:
compressible flow simulations, In High performance computing, 139, John Wiley and Sons Ltd.
(1992).
K. Fujii and S. Obayashi, J. Aircraft. 26, 1123 (1989).
S.K. Lele, J. Comp. Phys. 103, 16 (1992).
S. Sekimoto, K. Fujii, S. Hosokawa, and H. Akamatsu, Sensors Actuators, A Phys. 306, 111951(2020).
H. Aono, S. Kawai, T. Nonomura, M. Sato, K. Fujii, and K. Okada, AIAA J. 55, 3789 (2017).
M. Sato, K. Okada, K. Asada, H. Aono, T. Nonomura, and K. Fujii, Phys. Fluids 32, 025102 (2020).
K. Asada, Y. Ninomiya, K. Fujii, and A. Oyama, in 47th AIAA Aerospace. Sci. Meet. Reston,
Virginia, (2009), pp. 1–14.