SENAC Vol. 53, No. （1 2020. 1）― 15 ―［共同研究成果］

DBD プラズマアクチュエータを用いた自動車後流制御による抵抗低減― HPC 利用による高解像度シミュレーションに基づいて ―藤井 孝藏：東京理科大学工学部情報工学科
浅田 健吾：東京理科大学工学部情報工学科

1.

はじめに

現代社会において自動車は生活に欠かすことのできない移動手段の一つとなっている．19 世紀後半から開発されてきた自動車は各種要素技術の発展とともに性能は向上し，近年発売される自動車の燃費は 20km/L 以上にも達している．しかしながら，昨今環境問題への意識が高まる中，さらなる燃費向上は依然として工学分野として取り組むべき課題である．自動車の燃費向上には空気抵抗の低減が不可欠であり，これまで形状を工夫することによって抵抗低減を実現してきた．
近年航空機のように滑らかな形状をした自動車が増加しているのはこのためである．一方で，自動車開発において，燃費以上にデザイン性，すなわち“見た目”は販売台数を左右する重要な要素である．ここで問題となるのは空力的に優れた形状が必ずしもデザインとして優れているとは限らないことである．このため，デザインの自由度を保ちながら抵抗低減を実現するには形状工夫によらない手段が必要である．
そんな中，近年航空工学分野において，新たに DBD(Dielectric Barrier Discharge：誘電バリア放電)プラズマアクチュエータ（以下 PA）と呼ばれる流体制御デバイス[1]を用いた抵抗低減手法が注目されている．PA は，対象とする流れ場に局所的な変動を与えることで大きな流れを制御するマイクロ流体制御デバイスで，2 枚の電極とそれに挟まれた誘電体からなる単純な構造を持つ（図 1）．全体の薄さは数百マイクロメートルと非常に薄く軽量であるため，取り付けることによる形状変化が極めて小さいのが特徴である．PA を駆動するには電極間に高周波・高電圧の交流電圧を印加する．電圧印加によって露出電極表面にプラズマが発生し，露出電極周辺に非定常な流体変動が誘起され，その変動が全体の流れ場に作用することで流体機器の性能を向上させる[2]．
PA は，単純な構造，少ない消費電力，既存のシステムに対して導入が簡単であるという特徴を持つため，航空機の翼[3-10]，風力発電タービンにおける剥離流れの制御[11,12]や，ガスタービンエンジンの漏れ流れ抑制[13]など，様々な分野において実用化にむけた基盤研究が進められている．図 2 に，PA を用いた翼周り流れの剥離制御効果を示す．翼前縁に設置された PA を駆動させることで，翼面から大きく剥離した流れが翼面に付着し，翼性能が向上する．

図 1: DBD プラズマアクチュエータの構造と翼型への設置例.

― 16 ―

SENAC Vol. 53, No. （1 2020. 1）

図 2: DBD プラズマアクチュエータを用いた NACA0015 翼周り流れの剥離制御（制御／非制御）[4].
以上のように PA は対象物の形状変更を必要としないにもかかわらず多くの利点を有するために自動車メーカを含めた研究者が試行的な実験を開始している[14-16]．しかしながら，基本的な流れ制御の方針自体が不明確なままであり，翼型以外への適用は未だ手探りの状況にある．そこ
で，本課題では，高解像度シミュレーションを行うことで自動車周り流れの特性を理解し，PA による流れ制御の指針を示すことを目指す．特に本研究では Ahmed model[17,18]と呼ばれる簡易自動車モデル周り流れを対象とし，Large-eddy simulation (LES)を行うことで，剥離流れ制御において効果的な PA 位置を明らかにする．Ahmed model 周り流れの PA 取り付け位置に関する議論は文献[17]において実験を用いてなされているが，本研究では高解像度シミュレーションを行うことでより詳細に PA による空気抵抗低減と流れ場との関係を明らかにする．

2.

問題設定

2 . 1 簡 略 化自 動 車 モ デ ル
図 3 に対象とする Ahmed model の概略図を示す．Ahmed model は簡略化した自動車モデルで，過去に多くの実験，シミュレーションが行われている．前方は R をもった鈍頭形状（フロント部）で，後方は自動車のリアガラスを模擬した角柱を切り落としたような形状（スラント部）になっている．図 4 に示す計算領域を含む各寸法はモデル高さ H で無次元化されており，それぞれ，
L1/H=2.928, L2/H=0.697, B/H=1.35, R/H=0.347, c/H=0.174, X1=7.3, X2=21, Y1=6.493, Z1=4.861
である．形状パラメータであるスラント部の角度αは 25 度とした．また，位置座標（x, y, z）
はスラント部後縁を x=0,モデル中央を y=0, 計算領域底面を z=0 とする．
図 5 に Ahmed model 周り流れの概略を示す．前方から衝突した流れが逆圧力勾配でわずかに剥離後，乱流遷移し，後方で大きく流れが剥離する．後方の上面からはスパン方向に軸をもつ大規模な渦が，物体側面と上面の角からは流れ方向に軸をもった大規模な渦が生成される．これらの大規模渦が低圧領域を形成し，本物体流れの空力抵抗を増加させている．本課題ではこれらの大規模渦をプラズマアクチュエータで制御することで空力抵抗低減を目指している．

図 3：簡略化自動車モデル（Ahmed model）.

DBD プラズマアクチュエータを用いた自動車後流制御による抵抗低減― 17 ―

図 4：Ahmed model と計算領域.

図 5：Ahmed model 周り流れ[18].

3.

計算条件

3.1

計 算 条件 お よ び 計 算 手 法

一様流とモデル高さ H 基準のレイノルズ数は文献[19]に合わせ，ReH = 2.0 x 10-5 とし，マッハ数は計算効率を考慮し，実際よりも高い 0.2 とした．比熱比は空気のγ＝1.4，プラントル数はPr＝0.72 とする．
表 1 に本研究で用いる計算手法を示す．３次元圧縮性 Navier-Stokes 方程式を支配方程式として空間差分に高解像度スキームである 6 次精度 compact 差分法[20]を 10 次精度 3 重対角フィルタとともに用いた．時間積分には ADI-SGS[21]を使用し，音速基準の無次元時間刻みを 5 x 10-4（クーラン数は約 11）とした．
PA による流れ場への運動量供給は，支配方程式に体積力項を導入することで模擬する．体積力のモデル化によるプラズマ効果と流体支配方程式のカップリングには，幾つかの種類が存在するが，本研究では比較的単純なモデルである Suzen と Huzng ら[22]の定常体積モデル（以下 Suzenモデル）を非定常に拡張した体積モデルを用いる．この修正モデルの数値計算への導入方法やその信頼性に関する詳細の記述，および数値計算の妥当性については既出論文を参照されたい[23,24]．アクチュエータの出力を決めるパラメータ Dc を 0.4 とした．この値はこれまで行ってきた翼流れ制御[8]において一様流 U∞とアクチュエータの誘起流速 UDBD の比が UDBD/U∞~1 となる現実的な値である．また，PA は連続駆動を用い，周波数は一様流と代表長さを基準とした無次元周波数で 60 である．

― 18 ―

SENAC Vol. 53, No. （1 2020. 1）

表 1: 計算手法.
支配方程式

3 次元圧縮性ナビエ・ストークス方程式

空間差分

6 次精度コンパクト中心差分 ＋3 重対角フィルター（αf＝0.40）

時間積分

2 次精度 ADI-SGS 陰解法 ＋内部反復 5 回乱流モデル

Implicit LES

3.2 計算格子
図 6 に示すように格子はモデル近傍の境界層を解像する格子（Zone 1：赤），外部領域（Zone 2：青）の計 2 ゾーンからなっている．計算領域は参照実験[25]の風洞と同等の領域とし，流入領域から総格子点数は約 4400 万点（表 2）である．物体表面および計算領域の低面のみ境界層の発達を解像するために滑りなし条件とし，その他，計算領域側面および上面は格子点数を削減するために滑り条件としている．

(a)

(b)
図 6：計算格子.

表 2：計算格子点数.
jmax x kmax x lmax

Number of grid points

Zone 1511 x 510 x 101

26,321,610

Zone 2619 x 201 x 141

17,543,079

Total

N/A

43,864,689

3.3 計算ケース
表 3 に本研究の計算ケースを示す．本研究ではルーフ部側に PA を設置した場合とスラント部側に PA を設置した場合についてどちらに配置する方がスラント部の剥離領域抑制に効果的であるか議論を行う．翼流れの剥離制御において，これまで PA 位置に関しての議論はおこなわれており，剥離点付近に PA を設置することが効果的であることがわかっている[8]．しかしながら翼型では

― 19 ―DBD プラズマアクチュエータを用いた自動車後流制御による抵抗低減

剥離点付近で形状は滑らかに変化するのに対して，今回対象としている Ahmed モデルでは剥離点付近（スラント部前縁付近）で形状が急激に変化する．そのため，剥離点の上流に PA を設置する場合と下流に設置する場合とで制御効果に違いがみられると予測される．そこで本研究では図 7に示すようにスラント部前縁の上流と下流に PA を設置したケースの LES を行い，どちらが剥離制御に効果的であるか議論を行う．表 3 にそれぞれの計算ケースにおける PA の露出電極端の位置（下流側）を示す．

(a)PA location 1(b) PA location 20100
Sx
図 7: PA の体積力分布.
表 3: 計算ケース.

計算ケース
PA 設置位置（x/H）

3.4

Baseline
N/A

PA location 1-0.705

PA location 2-0.692

計 算 機に 関 し て計算機には，東北大学サイバーサイエンスセンターのスーパーコンピュータ SX-ACE を利用した．
実際の計算では，計算領域を分割して各ノードに割り当てて計算を行う．利用ノード数は 8 ノードである．ノード内はスレッド並列（自動並列，一部 Open MP を利用），ノード間は MPI によるハイブリッド並列を用いた．1 ケースの計算時間はおよそ 70 時間である．

4.

計算結果

4 . 1 基 本 流れ
今回行った計算の検証として，Ahmed model スラント部分の圧力係数（Cp）分布を Lienhart とBecker の風洞試験データ[25]と比較する．図 8 にスパン中央部（y/H=0），中間部（y/H=0.35），
スパン角部付近（y/H=0.63）の Cp 分布を示す．いずれの場所においても x/H=-0.7 付近に大きな負圧部分が見られるが，これはスラント部が始まる位置であり，流れが曲げられ加速されるためである．各位置において実験値とのわずかな差異はみられるが圧力分布のおおよその傾向は予測できており，得られた結果は妥当であると考えられる．
図 9 に Baseline ケースの瞬間流れ場を速度勾配テンソルの第二不変量による渦構造を示す．流― 20 ―

SENAC Vol. 53, No. （1 2020. 1）

れは図中左下から流入し，フロント R 部が終わる付近で側面，上面ともに剥離する．下流で再付着後乱流へと遷移し，下流へ行くにしたがい乱流境界層は発達していく．スラント部にさしかかると上面で流れは再度剥離し，角部からは長手方向に軸をもつ大規模な渦が生成される．これらの特徴は Ahmed ら[18] によって報告されている内容と一致する．
次節からはこの流れを基準とし，PA による流れ制御を行い，PA 設置位置が流れ制御効果に与える影響を議論する．

<Cp>ty/H=0
Baseline Exp.
Baseline Comp.

-2-10-0.8-0.7-0.6-0.5-0.4-0.3-0.2-0.10<Cp>ty/H=0.35-2-10

Baseline Exp.
Baseline Comp.
-0.8-0.7-0.6-0.5-0.4-0.3-0.2-0.10<Cp>ty/H=0.63
Baseline Exp.
Baseline Comp.

-2-10-0.8-0.7-0.6-0.5-0.4-0.3-0.2-0.1x/H

図 8: スラント部圧力係数分布（Baseline）.

0u/U∞

1.4

図 9: 瞬間流れ場（baseline）.

0

― 21 ―DBD プラズマアクチュエータを用いた自動車後流制御による抵抗低減

4.2 PA の設置位置と剥離制御効果
表 4 に各ケースにおける抵抗係数を示す．基準となる制御なし流れ（Baseline）の抵抗係数に対して，PA をスラント部前縁上流に配置したケース(PA location 1)は 10.8%抵抗係数が減少している．一方で，スラント部前縁より下流に PA を配置したケースでは効力係数は 13.1%減少しており，スラント角部より下流に PA を配置したほうが抵抗低減には効果的であることがわかる．
文献[15]において本研究と同様にスラント部前縁前後に PA を配置した場合の抵抗低減効果に関する議論がなされているが，スラント部前縁の上流に PA を配置したほうが効果的であると今回とは異なる結果が報告されている．考えられる原因のひとつは PA 取り付け位置の微妙な違いである．今回行ったシミュレーションではスラント部前縁の上流に PA を配置したケース（PA location1）の露出電極の下流端はスラント部前縁とは一致しておらず，体積力の最も強く発生する被覆電極の一部はルーフ部に存在する．そのため，PA によって誘起される体積力の方向は主にルーフ部の接線方向である．一方で，文献[15]では，露出電極の下流端はほぼスラント部前縁に一致しており，被覆電極は完全にスラント部に沿って配置されている．そのため，体積力は主にスラント部の接線方向に誘起されており，今回の PA location 1 のケースよりもむしろ PA location 2 のケースに近い体積力が誘起されているのではないかと推測される．このような誘起体積力の微妙な方向の違いが文献[15]と異なる結果を示した原因なのではないかと考えられる．

表 4: 抵抗係数.

抵抗係数（CD）
ΔCD [%]

Baseline

PA location 1

PA location 20.3824

0.34120.3323

N/A-10.8-13.1

図 10 に各ケースのスラント部圧力係数分布を示す．PA による制御ありのケースはいずれもx/H=-0.7 から x/H=0 にかけて Baseline ケースよりも負圧の値が抑えられているのがわかる．このように PA によって剥離領域が抑制されることでスラント部の負圧が抑えられ，その結果表 4で抗力が減少したと考えられる．また，PA ありの 2 ケースを比較した場合，PA をスラント角部より下流に設置した方（PA location 2）が負圧の値を抑えられているのがわかる．
図 11 にスラント部の摩擦係数分布を示す．摩擦係数は壁面垂直方向の速度によって正負が変わるために正であれば流れは付着しており，ゼロまたは負である場合，流れは剥離していると判断できる．Baseline ケースでは x/H = -0.65 から x/H = -0.2 付近まで負の値をとり，流れが剥離しているのがわかる．また，PA をスラント角部上流に配置したケース（PA location 1）の剥離位置は x/H = -0.65 付近で Baseline ケースと同等であるが，再付着点は x/H = -0.35 付近と剥離領域が縮小しているのがわかる．一方で PA を下流に設置した PA location 2 のケースでは剥離位置は x/H=-0.6 付近に後退し，再付着点も x/H = -0.45 付近へと上流側に移動し，剥離領域が最も小さくなっているのがわかる．

― 22 ―

SENAC Vol. 53, No. （1 2020. 1）

-4

0.06

Baseline Comp.
PA location 1
PA location 2-3.5-30.04

-2

<Cf>t

<Cp>t-2.5

-1.5
-1

0.03
0.02
0.01

-0.5

000.5

Baseline Comp.
PA location 1
PA location 2

0.05

-0.8

-0.7

-0.6-0.5-0.4-0.3-0.2-0.1-0.010-0.8-0.7x/H-0.6-0.5-0.4-0.3-0.2-0.10x/H

図 10: スラント部圧力係数分布（y/H=0）.

図 11: スラント部摩擦係数分布（y/H=0）.

図 12 にスラント部を側面から見た時間平均流れ場を示す．長手方向（x）速度 u で面塗りし，流線とともに示している．図中，青で示される部分は剥離領域である．図 11 で議論したようにBaseline ケースではスラント角部から剥離した流れがスラント後縁付近まで再付着せず，大きな剥離領域が形成されているのがわかる．これに対して PA を用いた場合，剥離領域は縮小しており，特にスラント角部より下流に PA をつけた場合に青色で示される剥離領域はほとんどなくなっているのがわかる．

(a)Baseline(b)PA location 10u/U∞(c) PA location 21.4

図 12 時間平均流れ場（y/H=0）.

4.3

今 後 の研 究 計 画

本研究では PA の駆動方法として連続駆動のみに着目し，効果的な PA 設置位置に関して議論を行なった．翼流れをはじめいくつかの流れ場の制御に関しては，デューティサイクルにより電圧印加のオンオフを繰り返すバースト駆動がより効果的であることが示されている．今回の応用に関して，バースト駆動を行なった場合に今回の知見が有効であるかどうかは今後の検討による．
また，Ahmed モデルの後流は 3 次元的に複雑な流れであるのに対し，本研究では PA の設置方法はスパン方向に一様な単純な設置方法であり，最適な設置位置に関しても検討の余地がある．今後は後流制御に対してより有効な PA の駆動方法と設置位置との組み合わせを検討していく予定である．

5.

まとめ

高解像度スキームを用いた自動車モデル流れの LES を実施し，DBD プラズマアクチュエータ（PA）

DBD プラズマアクチュエータを用いた
自動車後流制御による抵抗低減― 23 ―による後流制御において効果的な設置位置の検討を行なった．
まず，基準となる制御なし流れのシミュレーションで得られたスラント部の圧力係数を実験値と比較し，妥当な結果が得られていることを確認した．また，乱流遷移や長手方向に軸をもつ大規模な渦等，特徴的な流れ場を予測できていることを示した．設置位置のことなる PA を用いた流れ制御を行なった 2 ケースは両者とも制御なしのケースよりも抵抗係数を 10%以上低減できたが，今回用いた PA の配置に関しては PA をスラント角部より上流に配置したものよりも下流に配置したものの方が剥離を抑制し，抵抗を低減できることがわかった．
今後は PA の設置位置や駆動条件をより広く調査し，同デバイスのより効果的な利用方法を提案していく．

謝辞
ここに記載させて頂いた成果は,東北大学サイバーサイエンスセンターのスーパーコンピュータを利用することで実現することができたものである．また，研究にあたっては同センター関係各位に有益なご指導とご協力を頂いた．センターの皆様にこの場を借りて謝意を示したい．

参考文献[1] T. C. Corke, C. L. Enloe, and S. P. Wilkinson, “Dielectric Barrier Discharge Plasma
Actuators for Flow Control,” Annual Review of Fluid Mechanics, Vol. 42, pp. 505–529, 2010.
[2] T. C. Corke, M. L. Post, and D. M. Orlov, “Single Dielectric Barrier Discharge Plasma
Enhanced Aerodynamics: Physics, Modeling and Applications,” Experiments in Fluids, Vol.
46, No. 1, pp. 1–26, 2009.
[3] M. L. Post and T. C. Corke, “Separation Control on High Angle of Attack Airfoil Using
Plasma Actuators,” AIAA Journal, Vol. 42 No.11,2004.
[4] N. Benard, J. Jolibois, and E. Moreau, “Lift and Drag Performances of an AxisymmetricAirfoil Controlled by Plasma Actuator,” Journal of Electrostatics, Vol. 67, No. 2-3, pp. 133–139, 2009.
[5] 二宮由光, 藤井孝藏, “プラズマアクチュエータを用いた失速制御における周波数の影響,” 第45 回飛行機シンポジウム講演集, 1F1, October 2007.
[6] K. Asada, Y. Ninomiya, A. Oyama, and K. Fujii, “Airfoil Flow Experiment on the Duty
Cycle of DBD Plasma Actuator,” 47th AIAA Aerospace Sciences Meeting, The New Horizons
Forum and Aerospace Exposition, Orlando, Florida, January 2009.
[7] K. Fujii “High-Performance Computing Based Exploration of Flow Control with MicroDevices” Philosophical Transaction A, The Royal Society, Vol. 372, Article ID 20130326, 2014.
[8] M. Sato, H. Aono, A. Yakeno, T. Nonomura, K. Fujii, K. Okada, and K. Asada,“Multifactorial Effects of Operating Conditions of Dielectric Barrier Discharge Plasma
Actuator on Laminar Separated Flow Control,” AIAA Journal, Vol. 53, No. 9, 2015.
[9] H. Aono, S. Kawai, T. Nonomura, M. Sato, K. Fujii and K. Okada, “Plasma-Actuator
Burst-Mode Frequency Effects on Leading-Edge Flow-Separation Control at Reynolds
Number 2.6·105,” AIAA Journal Vol. 55, pp. 3789-3806, 2017.
[10] K. Fujii, “Three Flow Features behind the Flow Control Authority of DBD Plasma
Actuator: Result of High-Fidelity Simulations and the Related Experiments,” Applied
Science 2018, Vol. 8, Issue 4, 2018.
[11] H. Matsuda, M. Tanaka, S. Goshima, K. Amemori, M. Nomura and T. Osako,“Experimental Study on Plasma Aerodynamic Control for Improving Wind Turbine
Performance,” Asian Congress on Gas Turbins 2012, Shanghai, P. R. China, August 2012.

― 24 ―SENAC Vol. 53, No. （1 2020. 1）[12] D. Greenblatt, A. B. Harav, and H. M. Vahl, “Dynamic Stall Control on a Vertical Axis
Wind Turbine Using Plasma Actuators,” AIAA Journal, Vol. 52, No. 2, pp. 456–461, 2014.
[13] D. P. Rizzetta and M. R. Visbal, “Simulation of Plasma-based Flow Control Strategies forTransitional Highly Loaded Low- Pressure Turbines,” 37th AIAA Fluid Dynamics Conference
and Exhibit, Fluid Dynamics and Co-located Conferences, Miami, Florida, June 2007.
[14] K. Shimizu, T. Nakajima, S. Sekimoto, K. Fujii, T. Hiraoka, Y. Nakamura, T. Nouzawa,J. Ikeda and M. Tsubokura, “Aerodynamic drag reduction of a simplified vehicle model bypromoting flsw separation using plasma actuator, JSME Mechanical Engineering Letters,
Bulletin of the JSME, Vol.5, No. 19-00354,2019.
[15] Z. Hui, X. Hu, P. Guo, Z. Wang and J. Wang, “Separation Flow Control of a Generic
Ground Vehicle Using an SDBD Plasma Actuator,” MDPI, Open Access Journal, vol. 12, issue
20, pp. 1-14, 2019.
[16] S. Shadmani, S. M. Mousavi Nainiyan, M. Mirzaei, R. Ghasemiasl and S. G. Pouryoussefi,“Experimental Investigation of Flow Control over an Ahmed Body using DBD PlasmaActuator,” Journal of Applied Fluid Mechanics, Vol. 11, No. 5, pp. 1267-1276, 2018.
[17] C. Hinterberger, M. García-Villalba, W. Rodi, Large Eddy Simulation of flow around theAhmed body. In “Lecture Notes in Applied and Computational Mechanics / The
Aerodynamics of Heavy Vehicles: Trucks, Buses, and Trains”, R. McCallen, F. Browand, J.
Ross (Eds.), Springer Verlag, ISBN: 3-540-22088-7, 2004.
[18] S. Ahmed, G. Ramm, and G. Faltin, "Some Salient Features of The Time-Averaged
Ground Vehicle Wake," SAE Technical Paper 840300, 1984.
[19] S. Krajnović, L. Davidson, “Flow Around a Simplified Car, Part 1: Large Eddy
Simulation.” ASME. J. Fluids Eng., Vol. 127, Issue 5, pp. 907-918, 2005.
[20] S. K. Lele, “Compact finite difference schemes with spectral-like resolution,” J. Comput.
Phys., Vol. 103, pp. 16–42, 1992.
[21] K. Fujii, “Simple Ideas for the Accuracy and Efficiency Improvement of the Compressible
FlowSimulationMethods.” Paper presented at International CFD Workshop for Super- Sonic
Transport Design, Tokyo, March, 1998.
[22] Y. B. Suzen and P. G. Huang, “Simulations of Flow Separation Control using PlasmaActuators,” 44th AIAA Aerospace Sciences Meeting and Exhibit, Aerospace Sciences Meetings,Reno, Nevada, January 2006.
[23] K. Asada, T. Nonomura, H. Aono, M. Sato, K. Okada, K. Fujii, “LES of Transient Flows
Controlled by DBD Plasma Actuator over a Stalled Airfoil,” International Journal ofComputational Fluid Dynamics, Vol. 29, 2015.
[24] H. Aono, S. Sekimoto, M. Sato, A. Yakeno, T. Nonomura, and K. Fujii, “Computationaland Experimental Analysis of Flow Structures Induced by a Plasma Actuator with Burst
Modulations in Quiescent Air” Bulletin of the JSME Mechanical Engineering Journal Vol. 2,No. 4, 2015.
[25] H. Lienhart, and S. Becker, “Flow and Turbulent Structure in the Wake of a Simplified
Car Model,” SAE paper no. 2003-01-0656, 2003.