第 26 回数値流体力学シンポジウム
講演番号 D10-3

DBD プラズマアクチュエータを用いた翼剥離流れ制御の非定常流れ場解析
Unsteady flow analysis of DBD plasma actuator-based separation control around an airfoil○ 佐藤允，ISAS/JAXA, 神奈川県相模原市中央区由野台 3-1-1, sato@flab.isas.jaxa.jp
岡田浩一，菱友システムズ, 愛知県名古屋市港区大江町 6-19, okada@flab.isas.jaxa.jp
阿部圭晃，東大・院, 神奈川県相模原市中央区由野台 3-1-1, abe@flab.isas.jaxa.jp
青野光，ISAS/JAXA, 神奈川県相模原市中央区由野台 3-1-1, aono@flab.isas.jaxa.jp
野々村拓, ISAS/JAXA, 神奈川県相模原市中央区由野台 3-1-1, nonomura@flab.isas.jaxa.jp
藤井孝藏，ISAS/JAXA, 神奈川県相模原市中央区由野台 3-1-1, fujii@flab.isas.jaxa.jp
Makoto Sato, ISAS/JAXA, 3-1-1, Yoshinodai, Chuo-ku, Sagamihara, Kanagawa, 252-5210
Koichi Okada, Ryoyu systems Co.ltd, Ryoyu Systems, 6-19 Oye-cho, Minato-Ku, Nagoya , Aichi, 455-0024
Yoshiaki Abe, The University of Tokyo, 3-1-1, Yoshinodai, Chuo-ku, Sagamihara, Kanagawa, 252-5210
Hikaru Aono, ISAS/JAXA, 3-1-1, Yoshinodai, Chuo-ku, Sagamihara, Kanagawa, 252-5210
Taku Nonomura, ISAS/JAXA, 3-1-1, Yoshinodai, Chuo-ku, Sagamihara, Kanagawa, 252-5210
Kozo Fujii, ISAS/JAXA, 3-1-1, Yoshinodai, Chuo-ku, Sagamihara, Kanagawa, 252-5210
The flow-fields controlled by DBD plasma actuator around the NACA0015 airfoil are simulated with implicit large-eddy
simulation. The DBD plasma actuator is installed at the 5 % chord length from the leading edge and operated in the
burst modulation with the non-dimensional burst frequency (F+) of 1 and 6, burst ratio of 10 % and Dc=1.The
suppression of massive leading edge separation is observed in both non-dimensional burst frequencies. Common
feature in transient flow near the leading edge controlled by the DBD plasma actuator is two-dimensional vortices
generation from the separated shear layer. This leads the increase of the periodic component of Reynolds shear stressnear the shear layer region. For F+ =1, these vortices not only become three-dimensional structures by the vortex
paring but also generate the longitudinal vortices, which produce the non-periodic component of Reynolds shear stress.
While for F+ =6, these vortices are broken down by the transition and then forms three-dimensional structures.
Separation region near the leading edge is smaller in F+ =6. Hence better aerodynamic performance is attained in F+=6.
１．序論
流れの剥離は航空機の翼など流体機械の性能を低下させるだけでなく失速など航空機の安全性にも関わるため，その制御方法が古くから数多く研究されてきた．例えばボルテックスジェネレー
ターのような受動的な剥離制御は，装置の構造が簡易であるので，実機でも多く用いられている．しかし受動的な剥離制御は設計点
での流れは制御できるが，設計点外の条件下では抵抗を増加させてしまう等の問題がある．そのため，能動的な剥離制御の研究も多くなされてきた．特に定常な吸い込みや吹き出しを用いた定常ジェットの様な能動的な剥離制御はその有効性が示されていたが，機構が複雑で重量増加に見合う効果が得られないため，実用化された例は限られている．
これらの問題点を解決する方法として，近年スケールの小さいデバイスを用いて，局所的に小さな変動を流れに加える事で流れを制御する，誘電体バリア放電（Dielectric Barrier Discharge）プラ
ズマアクチュエータ(1),(2)（以下プラズマアクチュエータまたはア
クチュエータ）やシンセティックジェット(3),(4),(5)が世界中で盛んに研究されている．これらのデバイスは少ない投入エネルギで流れ場を大規模に変えることができるので,これまでのように翼型な
どの形状工夫に頼らずに,流体機械を設計できる可能性がある．例
えば平板にアクチュエータを付けて流れ場を制御することで，翼型相当の流れ場を実現することが可能になり，従来の様な翼型の設計が不要となる．
本研究では，電気的な駆動が可能で，構成装置が単純，素早い応答性という特徴があるプラズマアクチュエータに注目した．ジェットフラップやシンセティックジェットのように翼表面に穴を設ける必要がなく構造面でも有利である．プラズマアクチュエー
タは Fig. 1 に示すように誘電体を挟む 2 つの電極で構成される．

両電極に交流高電圧をかけることにより非定常な速度変動を流体に誘起し，時間平均を取ると露出電極から被覆電極へと向かう空
気の流れ（約数 m/sec）が誘起される．プラズマアクチュエータについての研究は実験，数値シミュレーションともに多く行われている．基礎研究として，プラズマアクチュエータが雰囲気流体に速度を誘起するメカニズムに関する研究(6)，誘起速度に関する研究(7)等がある．応用研究としては，翼面上でのプラズマアクチ
ュエータの設置位置の剥離制御に及ぼす影響に関する研究(8)，Fig.
2（横軸に時間，縦軸に電圧を示す）に示すように入力電圧波形として、バースト波を用いて剥離制御をより効率的に行おうとする
研究(11),(12),(13),(14)などが挙げられる．
バースト波形は入力電圧の基本
周波数より低周波数のステップ入力波形を掛け合わせ，オンとオ
フを周期的に切り替えるもので，バースト周波数を適切に選ぶ事
で通常の交流電圧入力では剥離制御できなかった迎角でも，投入エネルギを抑えつつ剥離制御が行えることが報告されている(11),(12),(13),(14)．

Fig. 1 Configuration of DBD plasma actuator.
1

Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 D10-3   uk(1)
0
 t  xk ui  ( ui uk  pik ) 1   ik(2)
 Dcqc E it xkRe  x k  e  ((e  p)uk ) 1  ul  kl


(3)
t
x
Re  xTbase=1/fbase
Ton (DBDon)

Toff (DBDoff)

Voltagetimek
T=1/f+k1
 qk

 Dcqc uk E k( 1)PrReM 2  x k

ui, qi, ρ, p, e, τij, δij 及び t はそれぞれ無次元化された速度ベクトル，
熱流束ベクトル，密度，圧力，単位体積あたりの全エネルギ，応 クロネッカーのデルタ及び時間である．式(2), (3)で，力テンソル，
右辺最後の項はプラズマアクチュエータが単位体積あたりに流体へ与える体積力 (DcqcEi) とエネルギ (DcqcukEk) をそれぞれ表している．プラズマアクチュエータによる項については次節で詳しく述べる．
式中に現れている 3 つの基本的な無次元パラメータ Re, M∞, Prはレイノルズ数，一様流マッハ数，プラントル数であり，以下のように定義される．
c
 uc
u
(4)
Re    , M   , Pr   pak
ここで μ, c, a, c p そして k は，それぞれ粘性係数，翼弦長，音速，
定圧比熱であり，添字“∞”は，一様流状態を表している．実際
の計算では，これらを一般座標へと変換したものを用いた．

Fig. 2 Burst wave image.
我々のグループでも現象理解を目的として，計算と実験両方のアプローチで研究を進めている(12),(13),(14),(15),(16), (17),(19),(20)．限られた条件下ではあるものの現象理解が進み，バースト波を用いた剥離制
御には少なくとも 2 つの剥離制御メカニズムがあること(14),(16)，必要に応じた剥離制御方法を選択できる様な評価指標を選ぶ必要があること(16)が分かっている．最近では，次世代スーパーコンピュ
ータ戦略分野のプロジェクト(16)において，「形状工夫による流体機器設計」を越える「マイクロデバイス流体制御を活用した流体機器設計」という新たな概念を提案し，実用化に向けた道を示す
ために積極的に研究を進めている．具体的には，デバイス作動パ
ラメータの最適化や実機を想定した高レイノルズ数，高動圧，高マッハ数での制御効果を理解する上で，低レイノルズ数領域での解析は重要になると考え，著者ら(19),(20)がレイノルズ数 Re=6.3×
104 の NACA0015 翼型に対して，プラズマアクチュエータの位置や作動条件を変えた剥離制御の LES 解析を用いた大規模パラメ
トリックスタディを行い，プラズマアクチュエータの剥離抑止メ
カニズムだけでなく制御デバイスの設計という観点からも様々な知見を得ている．例えば，プラズマアクチュエータの取り付け位
置や作動パラメータである投入運動量の大きさと剥離抑止効果の感度を明らかした．またこのレイノルズ数領域ではプラズマアク
チュエータにより乱流遷移が促進されることが剥離制御に重要であることを明らかにしている．
特にバースト波の重要なパラメータであるバースト周波数 f+を
ストロハル数と同様に主流速度u∞と翼弦長c で無次元化した無次元バースト周波数 F+ の最適な値に関する研究が多くなされている． 例えば，Corke らは F+=1 付近の値が剥離制御に最も効果的であると主張(9),(10)し，一方で F+=10 付近の値がより効果的であるという報告(11),(12), (13)もなされている．
このようにアクチュエータによる流れ制御のメカニズムについての現象理解は進んできているが，依然不明な点が多く，特にバースト波を用いた剥離制御は有用性が知られているにも関わらず，
未だに詳細な剥離制御メカニズムは分かっていない．そこで，本研究では著者らの解析(19) ,(20)の中から特徴的な2ケース（F+=1と6）について注目し，プラズマアクチュエータのバースト駆動を用いた翼剥離制御が，どのような現象によってバースト周波数の最適値が決定されているのか及びその剥離制御メカニズムを明らかにする事を目的とし，非定常流れ場について周波数解析，位相平均による流れ場可視化，及び線形安定性解析などの解析を行い流れ場の理解を目指す．

（２）プラズマアクチュエータのモデル化
前節において，プラズマアクチュエータの効果が，運動量保存
則に流体に与えられる力として DcqcEi, エネルギ保存則に流体に与えられるエネルギとして DcqcukEk という形で支配方程式に含まれることを述べた．qc 及び Ek はアクチュエータが空間に作り出す電荷分布と電場ベクトルをそれぞれ表している．また Dc は電磁力に関する無次元数で次のように表される．
qc,ref E ref c qc,ref ref
(5)
Dc 

 u2u2
ここで，添字“ref”は代表値を表し，φは静電ポテンシャル（電圧）である．上式より，Dc の物理的な意味はアクチュエータから
流体に加えられる電気的なエネルギと流体の一様流状態における動圧との比である．つまり，入力電圧の大きさが大きい程 Dc の値
は大きくなるため，Dc は入力電圧の大きさを表すパラメータとして用いる事ができる．本研究においても入力電圧の大きさを表す
のにこの Dc を用いることにする．具体的に Dc の値を決めるには
qc,ref とφref の値を決める必要があるが，φref は入力電圧の最大値等で簡単に決める事ができるのに対して qc,ref を与えるためには空間の電荷分布を計測等によりあらかじめ知っておく必要がある．しかし，qc,ref の値を計測するのは困難であるため，本研究では一様流に対して誘起流速がどれほど生じるかという事からこの Dc の値を現実的な値として見積もることにする．
次に実際に計算を行うにはqc及びEiの空間分布を求める必要があるが，本研究では Suzen らによって提案された体積力モデル(21)を用いた．この Suzen らによるモデルは Dc の値を適切に決める事で，ある程度信頼性のある誘起速度分布が得られる(8)．Fig. 3 にSuzen らによるモデルの体積力分布を示す．面塗りは力の絶対値
を，黄色い矢印は力の方向を示している．また，白い部分は露出
電極及び，被服電極である．このモデルは図の赤い矢印で示すように物体表面に向かう力と物体面に沿う 2 つの大きな力が働くのが特徴である．

２．解析モデル（１）支配方程式
支配方程式として 3 次元圧縮性 Navier-Stokes 方程式にプラズマアクチュエータによる効果をソース項として加えたものを用いた．
質量保存則，運動量保存則，エネルギ保存則からなり，一様流条
件及び翼弦長 c で無次元化を行った形で以下のように表される．
2

Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 D10-30（４）計算格子
計算格子を Fig. 5, Fig. 6 に示す．Zone1, 2 は各方向 2 点おき,
Zone3 は各方向 50 点おきに表示している．流体計算には Fig. 5 に示すような C 型格子を用い，外部境界は翼型前縁を中心に翼弦長
c に対して 25c とし，スパン長は 0.2c とした．計算は，翼周りの流れを解くための格子 Zone1 (青) とアクチュエータからの擾乱を解くための解像度の高い格子であるZone2 (赤) の2 つのゾーンからなる重合格子法を用いている．各格子間の物理量のやり取りには解強制置換法(31)を用いた．Zone2 には前節で述べたプラズマ
アクチュエータをモデル化した体積力をモデル用格子 (緑) から前処理として内挿している．その際，プラズマアクチュエータの
体積力分布は 2 次元仮定基づき得られたものであり，本研究ではスパン方向に同一分布を仮定している．格子点数は Zone1:
759×134×179 ， Zone2: 129×134×61 及 び ， モ デ ル 用 格 子 :1201(×134)×801 である．LES 計算は Zone1, Zone2 で行い，総格子
点数は約 2,000 万点程度，最小格子幅は Zone1, Zone2 ともに，境界層の厚み方向に 1.2×10-4 である．コンパクト差分法の空間解像度は高く，我々の経験では渦流れを例にとると，計算効率は 50
倍以上の格子点数を用いた汎用的な 2 次精度数値計算に匹敵する(32),(33) ,(34)．一番条件が厳しいと考えられる流れが付着している迎角
9 度の場合において，一般座標として翼面に沿った方向へ ξ, スパ
ン方向へ η, 翼面から垂直な方向へ ζ とすると，(Δξ+, Δη+, Δζ+min) =(8, 9, 1)である．これは Teramoto(29), Kawai ら(34)の乱流境界層を解像するのに必要な基準を十分満たす．内挿による空間精度の低下を防ぐため，各格子は接合部で格子を 10 点ずつ一致させている．
また，座標軸は翼の翼弦長方向に x， スパン方向に y， 垂直方向に z 軸とする．

21424

force magnitudeFig. 3 Force image of Suzen model.
また，本研究で使用するモデルの具体的な体積力分布を Fig. 4に示す．翼面に沿った方向に 1.5×10-2, 翼面から垂直な方向に
1×10-3 の領域をとっている．また，物体のちょうど中心に物体面があり，その上下に白で示されている露出電極と被服電極が存在する．体積力は物体表面から 5×10-3 の領域にのみ働いている．この領域を赤い矢印で示した点を中心に翼面の接線方向に物体面が重なる様に回転させ翼面上に配置する．また，配置する位置であるが，
前述の文献(15)において採用されている翼前縁から 5 %の位置に配置する．

0（５）計算条件
一様流マッハ数は計算効率の点から，圧縮性の無視できる範囲
でできるだけ大きい値を取るために M∞=0.2 とし， 一様流レイノ
ルズ数は翼弦長を代表長として Rec=6.3×104 とした．この値は参考実験(15)と同じである．また，迎角 α は前縁剥離を生じる α=12°とした．

5
force magnitudeFig. 4 Body force distribution of Suzen model (Dc=8).
（３）数値計算手法
ISAS/JAXAで開発された圧縮性流体解析ソルバーLANS3D(23),(24)を用いた．
以下に今回用いた手法を示す．
3次元圧縮性Navier-Stokes
方程式を支配方程式とし，一様流の音速a∞と代表長である翼弦長cで無次元化した．誘起速度領域が小さいプラズマアクチュエータの変動および剥離せん断層を効率的に解くために，メトリックス，
ジャコビアン，対流項および粘性項の離散化には6次精度のコンパ
クト差分法(25)を用いた．壁面付近で計算が不安定になり発散することを防ぐため，境界から1点目と2点目は通常用いられている片側4次精度の陽的差分ではなく，片側2次精度の陽的差分を用いた．
高周波の数値振動を抑えるために10次精度3重対角フィルター(26)を使用した（フィルター係数はαf=0.495）．時間積分法には2次精度3点後退差分に，内部反復(27)を5回入れたADI-SGS陰解法(28)を用いた．無次元時間刻み幅は，Teramotoの乱流境界層LES解析の研究(29)を参考に，最大クーラン数が1.7程度になる2.0×10-4とした．通常のLESであれば陽的なサブグリッドスケールが必要であるが，今回は乱流のサブグリッドスケールの効果は前述したフィルターが担うと考え，陽的なモデルを用いないImplicit LES(30)を用いた．
流出境界は静圧が静止気体に回復，その他は外挿とした．壁面は滑り無しの条件を用いた．また解析対象が2次元形状なので側面の両端は周期境界条件(10点重ね合わせ)を課した．

ξ
η
ζ

z
y

x

Fig. 5 Computational grids whole image.

Fig. 6 Computational grids near the leading edge: Zone1 (blue), Zone2(red) and model grid (green).
（６）プラズマアクチュエータの作動条件
プラズマアクチュエータの作動条件は前述の文献(9), (15)を参3

Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 D10-3
考にし，バースト波を用いた剥離制御において剥離制御効果が現れたと報告のあった 2 つの周波数について計算を行った．バース
ト波の周波数を f +，バースト波の基となる sine 波の周波数を fbaseとし，バースト比率 BR を
T
(6)
BR  onTとする．ここで，BR=100%は通常の交流入力を意味する．バース
ト波一周期中の波数を n とすると
f baseBR 1

n
T

f 

（６）解析データ取得区間
まず剥離制御をおこなっていない条件で 40 万 step（16 flow
through time）計算を行い,それを初期解としてアクチュエータを駆動させた.50 万 step（20 flow through time 後）までは準定常状態に達していないとみなし，そこから110 万step までの計60 万step (24flow through time)分を解析データ取得区間として用いる．DBD-offのケースも同様に処理した．

(7)
３．計算結果及び考察
３－１．時間平均流れ場
バースト周波数を変えた場合の空力性能及び流れ場の違いを示すために，時間平均流れ場の比較を行う．流れ場は時間方向だけなくスパン空間方向にも平均してある．
（１）空力性能の比較
まず Table 2 に DBD-off，F+=1，6 の 3 ケースの時間平均空力係数を示す．揚力係数，抗力係数，L/D ともに剥離制御を行ったことで大幅に改善されている．F+=6 のケースの揚力係数は F+=1 のケースと比べて同程度であるが，抗力係数は大きく改善している．
特に L/D で評価すると顕著である．過去の多くの研究は剥離制御効果の指標として揚力係数のみに注目している場合が多い．しかし本研究結果からも推測できるように，空力特性の評価として重要な抵抗係数や揚抗比も剥離制御効果の指標にいれた，プラズマ
アクチュエータの剥離制御の最適な制御パラメータの議論を行う必要がある．
Table 2 Aerodynamic coefficients.

となる．また f +を次の様に無次元化する．


F 

f c
u

(8)

次にこれらのパラメータを用いてプラズマアクチュエータを作動させた時の体積力分布の時間変動について定式化する．通常，
プズマアクチュエータを交流電圧により作動させた場合，交流一周期中に 2 回の放電が起こっていることが知られている．
そこで，本研究では，その 2 回の放電の両方で体積力が起こっていると仮定し，さらに両者とも流れ方向（露出電極から被服電極の方向）に同様の大きさの体積力が発生すると仮定する．具体的には時間変動する体積力分布を S(x,y,z,t)，時間平均した Suzen らの体積力空間分布モデルを SSuzen(x,y,z)で表すと次の様に表される．

S(x, y,z,t)  SSuzen(x, y,z)sin2 (2f baset)(9)

この式において sin2(2πfbaset)の部分は上記の仮定によるものである．
また，この力を時間平均した時の力を Save(x,y,z)とし，この値に対して前述したDc の値を定義するとSave x, y,z  

Ton
0Sx, y,z,t dt(10)1
 Dc SSuzenx, y,z
2

fbase

f+

BR

location

[Hz]

[Hz]

[%]

Fbase

F+

Cμ

CD

L/D

DBD-off

0.427

0.151

2.8

F+=1

0.842

0.069

12.2

F+=6

0.895

0.058

15.5

（２）圧力係数分布と摩擦係数分布
Fig. 7 に各ケースの翼面上圧力係数分布を示す．黒の線は制御なしのケース（DBD-off）を，青は F+=1 のケース，赤は F+=6 のケースを示している．DBD-off のケースは前縁から完全に剥離しており，翼弦から x/c=5 %の周辺からフラットな分布をとなる．
また，翼後縁の Cp の値も 0 から負圧側にずれており，前縁剥離を生じた流れ場特有のCp分布になる．
一方剥離制御した場合はF+=1，
6 共に前縁付近で効果的に剥離を抑えられおり，翼前縁で加速された流れにより作られる負圧のピークは大きな値となる．その後
翼上面の Cp 分布はフラットな領域を形成することから翼面上に層流剥離泡を形成していることが分かる．更に Cp 分布の傾きが急激に変わることから，F+=1 のケースは x/c=12％で，F+=6 のケー
スは x/c=15-20％付近で遷移していると予想される(35)．また F+=6のケースは前縁の曲率の大きい部分（x/c=0-5％）で鋭い負圧のピークを持ち，Table 2 に示したように揚力係数増加と抗力係数減少に寄与している．

で表される．Dc の値は文献(8)を参考にし，現実的な値として 1 とした．これは過去に我々のグループで行ってきた値である 8 よりも小さい．上記の体積力はfbase, BR, f+で表現可能である．今回計
算するケースは効果的に剥離制御が行えるケースとして報告されている F+=1 と F+=6 の 2 つの場合を計算し，両者の流れ場を比較する比較のため，アクチュエータがない条件（DBD-off）の結果
も合わせて示す．
BR の値は文献(15)で最も制御効果があったとされる 10 %とした．また，アクチュエータの取り付け位置は前述の様に文献(15)と同じ翼前縁から 5 %の位置に設置した場合で計算を行った．シンセティックジェットにはアクチュエータの生み出
す運動量の強さを決める無次元値である，無次元運動量係数比 Cμがある．無次元運動量係数比は，アクチュエータの生み出す運動量と一様流の運動量の比である．プラズマアクチュエータの場合
は，予め静止気体中にアクチュエータを駆動させた解析を行い，アクチュエータが生みだす体積力から値を概算した．F+=1 とF+=6の場合共にDc=1の場合でCμ=0.0086%とシンセティックジェ
ットと比べると小さい値である(4),(5)．詳細は文献(22)に示す．Table
1 に各作動条件をまとめる． 参考実験(13),(14),(15)の f +と fbase も示している
Table 1. Plasma actuator parameters.
DBD

CL

Dc

[%]

DBD-off

N/A

N/A

N/A

N/A

N/A

N/A

N/A

N/A

F+=1

5%

6,000

100

10

60

1

0.0086

1

F+=6

5%

6,000

600

10

60

6

0.00861

4

Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 D10-3-4.0（３）翼弦方向速度及び流線
Fig. 9に各計算ケースの翼弦方向（x方向）速度（x方向の流速uを一様流音速 a∞で割ったもの）と流線を拡大図と共に示す．灰色の流線が付着流れを，黒色の流線が循環流れをそれぞれ示す．赤色，オレンジ色，黄色，黄緑色，水色，紫色の点はそれぞれx/c=5，
10，15，20，25，30％位置を示す．DBD-OFFのケースは明らかに前縁から剥離しており，翼のほぼ全領域で大きな循環領域を生成している．またF+=1のケースは流れ場では前縁付近で剥離するものの，x/c=5-20％の領域に循環領域を形成しその後再付着していることから前縁失速を抑制できている．一方F+=6のケースは流れ場
では前縁付近で剥離するものの5-12 ％の領域にF+=1のケースよりも小さい循環領域を形成しその後再付着することから，Table 2に示したようにF+=1のケースよりも効果的に前縁失速を抑制できている．なおFig. 7に示した剥離泡と循環領域は一致する．

DBD-off
-3.0

F+=1
F+=6

Cp

-2.0

-1.0

0.0

1.0
0.0

0.2

0.4

0.6

0.8

x/c

Fig. 7 Surface pressure coefficient.
Fig. 8 に各ケースの翼上面における摩擦係数分布を示す．黒の線は制御なしのケース（DBD-off）を，青は F+=1 のケース，赤はF+=6 のケースを示している．off のケースは前縁から完全に剥離し，その後再付着しない．一方 F+=1，6 のケースは再付着している．
5％で小さいピークが立つのはアクチュエータが誘起する体積力の影響である．F+=1 のケースは 10-25％で，F+=6 のケースは10-15％で負をとる領域があるが，後に示すように剥離流れが形成する循環領域の影響である．また Cf のピーク値で傾きが正に変わり，その後高い正値を取ることから，F+=1 のケースは 12％で，F+=6 のケースは 20％付近で乱流遷移している．これは前述した
Cp 分布からの予測とおおよそ一致する．特に F+=6 で値が高いのは次節に示すように壁面付近に強い循環領域が存在するためである．その後ともに値が正値となり付着するが，F+=6 のケースの方が狭い領域で急激に値が変化する．また各ケースの剥離点，付着点及び循環領域を Table 3 に示す．剥離点及び付着点位置の定義
は摩擦係数がそれぞれ正から負へ，負から正へ 0 を通過する位置である．F+=1 のケースでは DBD-off のケースよりも遅く剥離し，24％位置で再付着するが，F+=6 のケースでは F+=6 と同程度の位置で剥離し 15％位置で再付着することから循環領域が小さい．

DBDoff

F+ = 1Act.

Act.

F+ = 6

Act.

0.02

Act.

DBD-off

-0.05

F+=1

0.01

Cf

F+=6

Fig. 9 Time-averaged and spanwise-averaged chord direction velocity
distributions and stream lines.

0.00（４）レイノルズ応力分布
各計算ケースにおけるレイノルズ応力の分布について議論する．
Fig. 10は各計算ケースのレイノルズ応力のうち-u’w’の分布を拡大図と共に示す．ここで，それぞれの速度変動は翼弦方向と翼面垂直方向成分であり，一様流音速で無次元化されている．レイノル
ズ応力の-u’w’成分は流れの混合の強い位置で大きな値をとり，一様流からせん断層への主流方向運動量供給の度合いを表している．
DBD-offのケースはx/c= 20 %から50 %付近の剥離せん断層でレイ
ノルズ応力が大きいことが分かる．この領域で運動量の交換はしているものの剥離点(x/c=3 %)より遠方でかつ翼面から離れているので，再付着には至っていない．またF+=1のケースではx/c=15 %から30 %の翼面に近い領域に，レイノルズ応力が高い領域が広く存在している． つまりDBD-offのケースと比較してアクチュエー
タ後方(x/c=15 -30 %)で運動量の交換が強く，より翼面に近いので剥離が抑えられている．一方F+=6のケースでは翼前縁から
x/c=10 %から15 %付近にF+=1のケースよりも前縁と壁面に近い領

-0.01

-0.02
0.0

0.2

0.4

0.6

0.8

x/c

Fig. 8 Surface skin friction coefficient.
Table 3 Separation and reattached locations.
Separation

Reattached

Recirculation

location

location

RegionDBD-off

2.8%

N/A

N/A

F+=1

3.5 %

24%

10-25%

F+=6

3.3 %

15%

10-15%

0.25
u/a ∞

5

Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 D10-3
域でレイノルズ応力が高い領域が狭く存在する．つまりDBD-offのケースやF+=1のケースと比較して，より前方で運動量の交換が起きるので剥離が抑えられている．

1.0 10-2
DBD-off
8.0 10-3

F+=1

maxDBDoffF+=6

6.0 10-3

4.0 10-3

2.0 10-3

F+ = 1

0.0
0.0

0.20

0.40

0.60

0.80

1.0

x/c

Fig. 11 TKE maximum evolution.

Act.

３－２．非定常流れ特性（１）瞬間流れ場
前節まで議論してきた時間平均流れ場や空力特性に対して，実際の瞬間流れ場がどの様な流れ場になっているか議論する．Fig.
12 は各ケースにおける瞬間場とその拡大図を示している．左半分における背景の面塗りは翼弦方向（x 軸方向）の速度を表し，翼の奥方向半面は翼面近傍での翼弦方向速度を表わしている．等値
面は速度勾配テンソルの第 2 不変量 Q を表わし，等値面の面塗り
は x 軸周りの渦度を表し時計周りが赤，反時計回りが青である．
また，スパン方向に引かれた赤色の線はアクチュエータの設置位置を示してしる．
アクチュエータによる制御を行わないDBD-OFFの流れは翼前縁から剥離し，剥離せん断層を形成しているのがわかる．またせん断層の不安定性からスパン方向に軸を持つ渦（V1，
V2）が放出され，下流に流れるにつれ渦同士の合体を起こし，その後崩壊して 3 次元的な乱流構造になり，下流では複雑な渦構造を形成する．また F+=1 のケースでは，アクチュエータを駆動したことによりせん断層から 2 次元的な渦（V1，V2，V3）が多く放出され，スパン方向に不安定になり（V2），渦同士の間でリブ構造の様な縦渦（ribs）が発生している．その後 2 次元的な渦構造
が縦渦の影響で崩壊して，3 次元的な渦構造になり最終的には乱流に遷移する．また翼の奥方向半面に示すように縦渦が発生している領域で，壁面付近で速度変動が大きい．結果として剥離制御を行わないケースよりも早い段階で流れが付着している様子がわ
かる．F+=6 のケースでは，F+=1 と異なりより乱流境界層に近い
非常に細かな 3 次元的な渦を伴う流れとなる．また F+=1 のケー
スよりも，せん断層から発生する 2 次元的な渦が早くスパン方向に不安定になる．更に翼の奥方向半面に示すように全体的に壁面付近で速度変動が大きいことも特徴である．
DBD-off や F+=1 のケースよりも早い段階乱流へ遷移することで流れが付着している．

Act.

F+ = 6Act.

Act.

-0.001

0.001-u’w’ /a∞2

Fig. 10 Time-and spanwise-averaged Reynolds stress distributions.
（５）乱流運動エネルギ分布
翼面の圧力係数と摩擦係数のプロットから乱流遷移を示したが，乱流遷移の位置を示すために，ζ 方向(翼から離れる方向)に対して乱流運動エネルギ(TKE)が最大値を示す場所の chord 方向位置を
Fig. 11 に示す．よってケースごとに同じ x 座標位置でも z 座標は異なる．ここで乱流運動エネルギ算出に用いた速度変動は一様流
音速で無次元化されている．黒の線は制御なしのケース（DBD-OFF）を，青は F+=1 のケース，赤は F+=6 のケースを示している．物体後流に生じる自由せん断乱流に関する過去の研究
から，流れ方向速度の変動値は徐々に増加していき，一度ピークを示した後減少していくことが示されている．この変化は流れが層流領域から線形領域，非線形領域，偶然化領域を経て乱流状態へと遷移することを示しており，本研究でもこの知見に則って乱流遷移していると考える．つまり DBD-off のケースは 30％，F+=1のケースは 15-20％，F+=6 のケースは 12％の各ピーク位置で乱流へ遷移していると考えられる．また F+=6 のケースでは，ピークが鋭く F+=1 よりも小さい TKEmax の最大値で乱流へと遷移していることから，より効率的に遷移していると考えられる．また後の位相平均解析で示すように，TKEmax の最大値をとる領域が渦の合体する領域にほぼ一致する．

6

Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 D10-3
している．F+=1，6 共に Fbase=60 で BR=10%なので，1 バーストにつき，F+=1 は 12 回，F+=6 は 2 回 Fbase の sine 波が入る．
F+=1 のケースはプラズマアクチュエータを駆動している1/60T
から 6/60T の 6 区間とアクチュエータを off にした 7/60T から
12/60T の 6 区間計 12 区間を示す．プラズマアクチュエータを駆動し始めた 1/60T から 2/60T までの区間はせん断層の x/c=15％か
ら25％の位置でアクチュエータ駆動よりも前に発生した渦が3個（SL1，2，3)が放出されるだけであるが，3/60T から 12/60T までの区間はせん断層ではアクチュエータを駆動したことにより
x/c=5％から 10％位置で渦が 6 個（V1，2，3，4，5，6）放出される．その後，それらが 10-20％の区間で 3 回合体（pairing1，3，4）する．(9/60T で 1 回，12/60T で 2 回)その後合体を繰り返した後に3 次元的な渦構造となり，崩壊して乱流へと遷移する．また 1/60Tのときに剥離せん断層から生じた 3 つの渦（SL1，2，3)はその後
2 つ目（SL2）と 3 つ目（SL1）が 9/60T で合体（pairing2）し，1
つ目（SL3）はそのまま合体することなく，下流にいくにつれて減衰する．
一方 F+=6 のケースはプラズマアクチュエータを駆動している
1/10T の 1 区間とアクチュエータを off にした 2/10T から 10/10T の
10 区間を示している．最後の 2 区間は最初の 1/10T と 2/10T を繰り返して表示している． DBD プラズマアクチュエータを駆動したことにより 3/10T でせん断層から渦が 2 個（V1，2）放出されて，渦が崩壊する．また 2 つ目（V1）はその後合体することなく，渦度が小さくなりながら移流して消えていくことから，3 次元的な渦へと崩壊して，乱流へ遷移していくと考えられる．
F+=1 は F+=6 よりも 1 バースト周期で渦が多く放出されるが，これは 1 バースト周期にそれぞれ入る sine 波の個数（1 バーストにつき，F+=1 は 12 回，F+=6 は 2 回 fbase の sine 波が入る．
）と関連がある．つまり F+と fbase の比がアクチュエータを駆動した時にできる渦の個数に，影響していると予想される．また F+の値を変えると，せん断層から発生する渦のスケールが異なるが，これは別途行っている実験結果や我々の経験と同様の傾向となる．

DBD-offribs
pairing

V1 V2

F+= 1

V3 ribs
V1V2

F+ = 60.3

0pairing

V1

u/a∞

-10

10

ωx

Fig. 12 Iso-surfaces of 2nd invariant of the velocity gradient tensors
and chord direction velocity distributions(Q/ a∞=10).
（２）位相平均流れ場
瞬間流れ場の理解を深めるために位相平均を用いた流れ場の可視化を行った．各物理量を--~f = f + f’ = f + f + f”

(12)-

として 3 成分分解した．ここで，瞬間値はf，時間平均はf，変動~
値は f’，周期変動成分fは，非周期変動成分は f”である．位相平均には後の周波数解析に示すように，流れ場の最も支配的な周波数である，プラズマアクチュエータの無次元バースト周波数 F+を用いた．Fig. 13，Fig. 14 示すのは，それぞれ F+=1，6 の場合の，位相平均した y 軸周りの渦度の面塗りと速度勾配テンソルの第 2 不変量の等高線を，各位相で時系列で示したものである．瞬間流れ場から，アクチュエータを駆動し始めてから終わるまでの区間で多くの渦構造が観察されており，この渦構造が流れ制御に対して重要と考えられる．位相平均を行う場合には，F+=1 の場合は 1
周期を 60 分割，F+=6 の場合は 1 周期を 10 分割してそれぞれの区間を位相平均した．これによりそれぞれの位相平均区間の時間ス
ケール（両ケース共に 1 位相平均区間の間隔は一様流速と翼弦長で無次元化した無次元時間で 0．083）が同じとすることで，現象を理解しやすくした．流れ場はそれぞれスパン方向にも空間平均してある．なお，速度勾配テンソルの第 2 不変量の正値は渦を，
負値はせん断を示している．赤色，オレンジ色，黄色，黄緑色，
水色，紫色の点はそれぞれ x/c=5，10，15，20，25，30％位置を示
7Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 D10-3
SL1V1

SL2 SL3

1/60T(DBDon)

1/10T(DBDon)

2/60T(DBDon)

2/10T(DBDoff)

V2
V1

V2

3/60T(DBDon)

3/10T(DBDoff)

4/60T(DBDon)

4/10T(DBDoff)

5/60T(DBDon)

5/10T(DBDoff)

6/60T(DBDon)

6/10T(DBDoff)

7/60T(DBDoff)

7/10T(DBDoff)

8/60T(DBDoff)

8/10T(DBDoff)

9/60T(DBDoff)

9/10T(DBDoff)

10/60T(DBDoff)

10/10T(DBDoff)

11/60T(DBDoff)

1/10T(DBDon)

12/60T(DBDoff)

2/10T(DBDoff)

V3

V4

V5

pairing1

pairing2

V6

pairing3-50

ωypairing4-50

50

Fig. 13 Phase and span averaged y-vorticity and second invariant of thevelocity gradient tensor for each phases at F+=1 (black lines, contour
range:0.1-50 with 20 lines)

ωy50Fig. 14 Phase and span averaged y-vorticity and second invariant of thevelocity gradient tensor for each phases at F+=6 (black lines, contour
range:0.1-50 with 20 lines)
8

Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 D10-3（３）レイノルズ応力の周期，非周期成分分解
Fig. 15，Fig. 16はそれぞれF+=1と6の場合における流れ場の，
時間方向及びスパン空間方向に平均したレイノルズせん断応力を周期成分と非周期成分に分けて拡大図と共に示す．等高線は速度
勾配テンソルの第 2 不変量である．周期成分は DBD プラズマア
クチュエータの駆動周波数 F+の周期である．Fig. 10 に示した時間方向およびスパン空間方向に平均したレイノルズ応力は，式(12)に示す．3 成分分解を用いることで，周期成分と非周期成分に分ける．赤色，オレンジ色，黄色，黄緑色，水色，紫色の点はそれ
ぞれ x/c=5，10，15，20，25，30％位置を示している．共に周期成
分よりも非周期成分の方が支配的であることから，剥離せん断層から発生するアクチュエータの周波数に同期した渦（F+=1 のpairing5 で生成された大規模な渦または F+=6 の V1）による運動量の交換作用はあまり強くないと言える．F+=1 のケースでは周期成分，非周期成分共に x/c=15 %から 30 %付近に値の高い領域があるが，非周期成分の方が高い領域が壁面に近い．一方 F+=6 のケースでは周期成分は x/c=10 %から 20 %付近に値の高い領域が，非周期成分では 15 %付近に値の高い領域がある．
次に Fig. 17 と Fig. 18 にそれぞれ F+=1 と 6 の場合のレイノルズ
せん断応力の周期成分と非周期成分を位相ごとに示す．等高線は速度勾配テンソルの第 2 不変量である．F+=1 は 6/60T から 21/60Tに示すように，周期成分が 2 次元的な大規模渦構造（黒線で示す
渦），特に F+=1 の周波数に同期した渦構造によって生成されている．一方非周期成分は 1/60T では流れ場が乱流へ遷移しているので全体的に生成されているが，6/60T から 21/60T では 2 次元的な大規模渦構造（特に 16/60T）と渦同士の間の領域によって強く生成されている．非周期成分が F+=1 の周波数に同期してない大規模渦構造で強くなるのはFig. 12 に示したように渦が2 次元的な渦が，スパン方向の不安定性から変形して 3 次元的な構造になるからであると考えられる．また，非周期成分が渦同士の間の領域で生成されるのは，これも Fig. 12 に示したように，渦同士の間の領域でリブ構造の様な縦渦が発生するためである．一方 F+=6 もF+=1 同様に周期成分は F+=6 に同期した 2 次元的な大規模渦構造（黒線で示す渦）によって生成されているが，流れが x/c=12%付
近で遷移するので渦構造が崩壊し，x/c=20%以降では F+=1 ほどは周期成分は強く生成されない．非周期成分も同様に流れ場が乱流へ遷移しているので全体的に生成されているが，F+=1 と異なり渦同士の間の領域ではあまり強く生成されず，渦構造自身の方が強く生成されている．これは x/c=12%付近で既に遷移しているために，Fig. 12 に示すように大規模な渦構造があるものの瞬間場では 3 次元的な渦構造をしている為である．

Periodic Reynolds shear stress
Act.

Act.

Nonperiodic Reynolds shear stress
Act.

Act.

~ a2
 u~w


-0.0010.001
 u "w" a2

Fig. 15 Phase and span averaged total component, periodical component
and nonperiodical component of Reynolds shear stress distributions at
F+=1.

Periodic Reynolds shear stress
Act.

Act.

Nonperiodic Reynolds shear stress
Act.

Act.

~ a2
 u~w


-0.0010.001
 u "w" a2

Fig. 16 Phase and span averaged total component, periodical component
and nonperiodical component of Reynolds shear stress distributions at
F+=6.

9

Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 D10-3
1/60T(DBDon)

1/60T(DBDon)PSD*  PSD 
Act.

Act.

6/60T(DBDon)

6/60T(DBDon)

11/60T(DBDoff)

11/60T(DBDoff)

16/60T(DBDoff)

16/60T(DBDoff)

21/60T(DBDoff)

21/60T(DBDoff)

Periodic Reynolds shear stress

Non-periodic Reynolds shear stress0.004

 u "w" a2

Fig. 17 Phase and span averaged total component, periodical component
and nonperiodical component of Reynolds shear stress and second
invariant of the velocity gradient tensor distributions at F+=1.
1/10T(DBDon)

1/10T(DBDon)

Act.

Act.

3/10T(DBDoff)

3/10T(DBDoff)

5/10T(DBDoff)

5/10T(DBDoff)

7/10T(DBDoff)

7/10T(DBDoff)

9/10T(DBDoff)

9/10T(DBDoff)

Periodic Reynolds shear stress

Non-periodic Reynolds shear stress

10-2

~ a2
 u~w


-0.004(13)

まず DBD-off のケースは x/c=15％位置で St=10 付近にピークが見えることから，剥離せん断層からこの周波数で渦が放出されていることが分かる．その後 x/c=20％位置で急激に高周波成分が卓越
し，x/c=25％位置で低調波成分である St=5 にピークが観察されるが，これは渦の合体によるものと考えられる．
また F+=1 のケースでは 5 箇所の計測位置すべてで駆動周波数であるSt=1 の成分が卓越している．
x/c=5％位置ではSt=120 でFbase（Fbase=60 であるが式（9）に示すように交流一周期中に 2 回の放電が起こっていると仮定している）のピークが明確に見える．
F+=1 の場合は 1 バースト周期で基本周波数の sine 波が 12 周期入るので，流体が基本周波数にある程度追従していると考えられる．
x/c=10％位置で，x/c=5％に比べて高周波が成長するのはアクチュ
エータを駆動することで x/c=5-10%の領域で剥離せん断層から細かい渦（Fig. 13 の V1，2，3，4，5，6）が発生するためである．
その後 x/c=15％位置で St=1 のピークが強くなるが，これは剥離せん断層から出てくる渦が合体をくり返して，最終的に 1 つの渦構造（pairing4 によって生成される渦）になるためである．またこの位置でスペクトルの傾きが-5/3 乗に近くなり，低周波から高周波まですべての成分が大きくなることから，渦が流れ方向に進む
につれ，2 次元的な渦が崩壊することで x/c=15％付近で乱流へ遷移する．
一方 F+=6 のケースでは駆動周波数である St=6 の成分以外にも高調波や低調波成分のピークが F+=1 よりも強く見える．x/c=5％
位置ではF+=1のようにSt=120でFbaseの明確なピークは見えない．
これは F+=6 の場合は 1 バースト周期で基本周波数の sine 波が 2
回しか入らないので，1 バースト周期で 12 sine 波が入る F+=1 の場合と異なり，流体が追従しきれていないためと考えられる．また x/c=5%では駆動周波数とその高調波のピークを持つこと以外
は層流流れのスペクトルの傾向に近いが，x/c=10%でスペクトルの傾きが急激に-5/3 乗に近くなり，低周波から高周波成分まですべての成分が大きくなることから，この付近で乱流へ遷移していると考えられる．これは渦崩壊によって生じる縦渦構造やスケー
ルの小さい乱流渦生成による偶然化過程を経て，さまざまな周波数のエネルギ成分が増加するためで，
Fig. 11 に示したように，
TKE
最大値の分布が増加後に急激に減少したことに対応する．また低調波のピークが見えることから渦構造が 3 次元的になって崩壊していると考えられる．

~ a2
 u~w


-0.0041 cu2 u

-5/3

0.004
 u "w" a2

10-4

Fig. 18 Phase and span averaged total component, periodical component
and nonperiodical component of Reynolds shear stress and second
invariant of the velocity gradient tensor distributions at F+=6.

10-6
10-8

5%
10%

（３）周波数解析
流れ場の周波数特性と渦構造を定量的に考察するために，周波数解析を行った．データはスパン空間方向に 5 点平均してある．
計測位置は 5，10，15，20，25％位置における乱流運動エネルギが最大となる位置で，Fig. 11 で用いた位置と同じである．Fig.
19,Fig. 20,Fig 21 に各ケースの各計測位置の翼弦方向（x 方向）速度変動のパワースペクトル密度 (PSD) を示す．
横軸はストロハル
数 St で，定義は F+と同じく一様流速と翼弦長で無次元化してある．縦軸は流れ方向速度変動のパワースペクトル密度(PSD)で，は以下の式に示すように，一様流速の 3 乗と翼弦長で無次元化した．

10-10
10-12

15%
20%
25%

10-14
0.1

1

10

100

1000

Fig. 19 Power Spectrum Density of streamwise velocity fluctuation at
DBD-off.
10

Copyright © 2012 by JSFM

10-2

第 26 回数値流体力学シンポジウム
D10-3
講演番号
講演番号
えられる．さらに，低調波成分を励起することで，不安定周波数で生成される渦のペアリングが促進され，運動量の混合・取り入
れを増幅させることも剥離抑制の一因として考えられる．ここでは示さないが，ほぼ同条件において F+=10,20 といった最も不安定な成分に近いバースト周波数を用いた場合よりも F+≒5 の方がより乱流遷移は促進され，空力性能が向上したことから，翼面剥離
では上記のようなメカニズムに基づいたアクチュエータの駆動条
件を設定することが鍵となるものと考えられる．

-5/3

10-4
10-6
5%

10-8

10%

10-10

15%

100.0

20%

10-12

25%

10-14
0.1

1

10

100

1000

FFT
60.0

Power Spectrum Density of streamwise velocity fluctuation at-iFig. 20
F+=1.

LSA

80.0

40.0

10-2

-5/3

20.0

10-4
0.0
0.0

10-6

20.0
St

30.0

40.0

5%

10-8

10%

10-10

Fig. 22

Space amplitude developing of DBD-off case.

15%

４．まとめ
レイノルズ数63,000，一様流マッハ 数0.2，迎角12 °におけるNACA0015 翼型周りに対して， DBDプラズマアクチュエータを適用した剥離制御を行い，位相平均解析や線形安定性解析，周波数解析などの非定常解析を用いて流れ場の理解を行った．DBDア
クチュエータは翼前縁から翼弦長のx/c=5%位置に配置し，従来よりも小さい投入エネルギであるDc=1の条件下において， F+=1,6のバースト周波数を用いた．共に剥離を抑えることができ，F+=6のケースはF+=1よりも空力特性が良いことが分かった．
F+=1のケースは流れが剥離した後に，
x/c=15%で乱流遷移することで生じる壁面付近のレイノルズ応力による運動量の交換作用で再付着することで，剥離を抑えることができた．次に瞬間流れ場
からF+=1のケースはアクチュエータの駆動により多くの渦を誘起し，それらが合体を繰り返すことで大規模な渦構造が生じ，さらに2次元的な渦の間の領域に縦渦（リブ構造）が発生することが分かった．また位相平均解析より，強い2次元的な渦がレイノルズせん断応力の周期的な成分を，2次元的な渦の間の領域に発生する縦渦と遷移によって生じるスケールの小さい渦が非周期的なレイノ
ルズせん断応力を生じさせ，特にこの非周期的な乱流成分が、周期的な成分と比べて支配的であることが分かった．
またF+=6のケースは流れが剥離した後にx/c=10%で乱流遷移することで生じる壁面付近のレイノルズ応力による運動量の交換作用で再付着するが，よりスムーズに乱流へ遷移することでF+=1のケースよりも早く再付着することが分かった．次に瞬間流れ場からF+=1のケースに比べて誘起された2次元的な渦構造は早く崩壊
し，スムーズに3次元的な大規模構造へと遷移することが分かった．
更に， 2次元的な渦がレイノルズせん応力の周期的な成分を，3次元的な大規模渦構造と遷移によって生じるスケールの小さい渦がレイノルズせん応力の非周期的な成分を生成するが，F+=1同様に非周期的成分が支配的であることが分かった．またDBD-offのケースの流れ場を用いて線形安定性解析を行い，剥離せん断層の不

20%

10-12

25%

10-14
0.1
Fig 21
F+=6.

10.0

1

10

100

1000

Power Spectrum Density of streamwise velocity fluctuation at（４）線形安定性解析
本節では，線形安定性解析を用いてバースト周波数と剥離せん断
層の不安定周波数の関係を調べる．剥離せん断層における乱流遷
移は変曲点不安定性から生じると考え，非粘性，非圧縮及び平行流を仮定した，一般的なレイリー方程式を用いて線形安定性解析を行った．解析で用いた基本流として，DBD-off の場合における流れ方向 10%位置の速度分布を tanh で近似した式を用いた．本解析では擾乱の時間発展項は固定し，周波数をパラメータとして擾
乱の空間増幅率を調べた．また，シミュレーション結果の周波数解析より，各周波数帯における流れ方向速度の空間増幅率を調べ，それぞれを比較した．Fig. 14 の横軸はストロハル数，縦軸は空間増幅率であり線形安定性解析の結果（LSA）とシミュレーション
結果（FFT）を示している．両分布の傾向は一致しており，最も空間増幅率の大きいストロハル数はLSAでSt=16,シミュレーショ
ンで St=14 となっている．ただし，シミュレーションでは St=6 付
近でもう一つの山が観察される．このピークの要因として，流れの非線形作用や翼面曲率の影響が考えられる．ただし，St=16 の高調波成分にはピークが見られないため，非線形性の影響は疑問の余地があり，今後の検討課題である．これらの結果と，本研究で示したバースト周波数の効果について考える．バースト周波数として F+=6 を用いた場合，St=6 付近の不安定となる成分を励起することに加え，その高調波である St≒15 付近も励起されるため，
不安定擾乱の成長がより促進され，乱流遷移が速やかに進むと考11

Copyright © 2012 by JSFM

安定周波数との関係からF+=6のケースが効果的に遷移を促進する
可能性を示唆した．
今後解析を行う予定であるより実用的な高レイノルズ数流れでは，DBD-off の場合は層流のまま剥離する低レイノルズ数流れと異なり，乱流へ遷移してから前縁ないし後縁から剥離すると予想される．よって以上のことから，乱流遷移によって流れを再付着させる F+=6 のケースよりも，渦の合体によって大規模構造を作ることで剥離を抑える F+=1 のケースの方が，剥離を抑えることができる可能性があると推測できる．

(14)

(15)

(16)
謝辞
宇宙航空研究開発機構宇宙科学研究所の焼野氏，東京大学大学
院の浅田氏及び関本氏から有意義なコメントを頂いた．本研究の解析は東京大学情報基盤センターのFX10 スーパーコンピュータ
システム（Oakleaf-FX）を使用した．また文部科学省「高性能汎用計算機高度利用事業「HPCI 戦略プログラム」分野４ 次世代
ものづくり」課題１「輸送機器・流体機器の流体制御による革新
的効率化・低騒音化に関する研究開発の支援」のもとに理化学研
究所計算科学研究機構が保有する京速コンピュータ「京」を使用した．ここに記して謝意を表する．

(17)

(18)
(19)

参考文献
(1) Post, M. L. and Corke, T. C., “Separation Control on High Angle of
Attack Airfoil Using Plasma Actuators,” AIAA Journal, Vol.42
No.11 pp.2177-2184, 2004.
(2) Corke, T. C., Lon Enloe, C. L. and Wilkinson, S. P., “Dielectric
Barrier Discharge Plasma Actuators for Flow Control,” Annual
Review of Fluid Mechanics, Vol. 42, pp. 505-529, 2010.
(3) Glezer, A. and Amitay, M.: Synthetic Jets, Annual Review of Fluid
Mechanics, 34 (2002), pp. 503-529.
(4) Okada, K., Nonomura, T., Fujii, K. and Miyaji, K: Computational
Analysis of Vortex Structures Induced by a Synthetic Jet to Control
Separated Flows, International Journal of Flow Control,（投稿中）(5) 阿部ら, “シンセティックジェットによる翼周り剥離流れの制御に関するLES：スパン方向渦の解析，” 第26 回数値流
体力学シンポジウム論文集，(2012).
(6) Font, G. and Morgan, W., “Plasma Discharges in AtmosphericPressure Oxygen for Boundary Layer Separation Control” AIAAPaper 2005-4632.
(7) Roth, J. R. and Dai, X., “Optimization of the Aerodynamic Plasma
Actuator as an Electrohydrodynamic (EHD) Electrical Device,”
AIAA Paper 2006-1203.
(8) Tsubakino, D. and Fujii, K., “Effective Lay-out of plasma
Actuators for a Flow Separation Control on a Wing,” AIAAPaper 2007-474.
(9) Corke, T. C., Chuan, H. and Patel, M. P., “Plasma flaps and slats :An application of weakly-ionized plasma actuators,” AIAA Paper,
2004-2127 (2004).
(10) Patel, M. P., Terry Ng, T., Vasudevan, S., Corke, T. C., Post, M. L.,
McLaughlin, T. E., Suchomel, C. F. “Scaling Effects of an
Aerodynamic Plasma Actuator,” AIAA Paper 2007-635.
(11) Sidorenko, A. A. et al., “Pulsed Discharge Actuators for
Rectangular Wings Separation Control,” AIAA 2007-941.
(12) Asada, K. and Fujii, K., “Computational Analysis of UnsteadyFlow-field Induced by Plasma Actuator in Burst Mode,” AIAAPaper 2010-5090.
(13) Sekimoto, S., Asada, K., Usami, T., Ito, S., Nonomura, T., Ooyama,
A. and Fujii, K., “Experimental study of effects of frequency for

(20)

(21)
(22)

(23)

(24)
(25)

(26)

(27)
(28)

(29)

(30)

(31)

(32)

(33)
12

第 26 回数値流体力学シンポジウム
講演番号 D10-3
burst wave on a DBD plasma actuator for separation control”,
AIAA Paper, 2011-3989 (2011).
Sekimoto, S.; Asada, K.; Anyoji, M.; Nonomura, T. & Fujii, K.
(2012), Comparative Study of Co-flow and Counter Blowing
DBD Plasma Actuators for Separated Flow over an Airfoil, in'AIAA 2012-1137'.
Asada, K., Ninomiya, Y., Oyama, A., Fujii, K., “Airfoil FlowExperiment on the Duty Cycle of DBD Plasma Actuator,” AIAAPaper 2009-531.
Asada, K., and Fujii, K., Burst Frequency Effect of DBD Plasma
Actuator on the Control of Separated Flow over an Airfoil, AIAA
paper 2012-3054.
Kaneda, I.; Sekimoto, S.; Nonomura, T.; Asada, K.; Oyama, A. &Fujii, K. (2012), 'An Effective Three-Dimensional Layout ofActuation Body Force for Separation Control', International
Journal of Aerospace Engineering 2012.
http://www.mext.go.jp/b_menu/houdou/21/07/attach/1282504.htm
Sato, M. et al., “Parametric study on separation control by DBD
plasma actuator over NACA0012 and NACA0015 airfoil at
Reynolds number 63,000”, 65th Annual Meeting of the APS
Division of Fluid Dynamics, Vol. 57, No. 17, (2012)
佐藤ら, “DBD プラズマアクチュエータを用いた翼剥離制御に関する大規模パラメトリックスタディ，” 第 26 回数値流
体力学シンポジウム論文集，(2012).
Suzen, Y. B. and Huang, P. G., “Simulations of Flow Separation
Control using Plasma Actuator,” AIAA Paper 2006-877.
青野,野々村,藤井,DBD プラズマアクチュエータモデルを用いた静止気体中における誘起流れの数値解析,日本機械学会
流体工学部門論文集,252（2012）
Fujii, K., “Progress and Future Prospects of CFD in
Aerospace-Wind Tunnel and Beyond”, Progress in Aerospace
Sciences, Vol. 41, No. 6, (2005), pp. 455-470.
Fujii, K., “CFD contributions to high-speed shock-related
problems”, Shock Waves, Vol. 18, No.2, (2008), pp145-154.
Lele, S.K., “Compact Finite Difference Scheme with Spectral-LikeResolution,” Journal of Computational Physics, Vol.103, pp. 16-22,
1992.
Gaitonde, D. V. and Visbal, R. M., “Pade Type Higher-Order
Boundary Filters for the Navier-Stokes Equations,” AIAA Journal,
Vol. 38, No. 11, pp. 2103-2112, 2000.
Chakravarthy, S. R., “Relaxation Methods for Unfactored Implicit
Upwind Schemes,” AIAA Paper 84 - 0165.
Nishida, H. and Nonomura, T., “ADI-SGS Scheme on Ideal
Magnetohydrodynamics”, Journal of Computational Physics Vo.
228, Issue 9, (2006), pp. 3182-3188.
Teramoto, S., “Large-Eddy Simulation of Transitional Boundary
Layer with Impinging Shock wave,” AIAA Journal, Vol. 38, No.
11, pp. 2103-2112, 2000.
Visbal, R. M. and Rizzetta, P. D., “Large-eddy simulation on
general geometries using compact differencing and filtering
schemes,” AIAA Paper 2002-288.
Fujii, K., “Unified Zonal Method Based on the Fortifies Solution
Algorithm,” Journal of Computational Physics, Vol.118, pp.92-108,
1995.
Arasawa, T., Fujii, K. and Miyaji, K., “High-Order Compact
Difference Scheme Applied to Double-Delta Wing Vortical Flows”,
Journal of Aircraft, Vol.41, No.4, (2004), pp. 953-957.
Kawai, S. and Fujii, K., “Analysis and Prediction of Thin-Airfoil
Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 D10-3
Stall Phenomena Using Hybrid Turbulent Methodology”, AIAA
Journal, Vol.43, No.5, (2005), pp. 953-961.
(34) Kawai, S. and Fujii, K. “Compact Scheme with Filtering forLarge-Eddy Simulation of Transitional Boundary Layer,” AIAA
Journal, Vol. 46, No. 2, Mar. (2008), pp. 690-700.
(35) Gerakupulos, R., Boutilier, M. S. H., and Yarusevych, S.,“Aerodynamic Characterization of a NACA 0018 Airfoil at Low
Reynolds Numbers,” AIAA 2010-4629, 2010.

13

Copyright © 2012 by JSFM