航空宇宙産業において、無人航空機(UAV: Unmanned Aerial Vehicle)や垂直離着陸機(VTOL: Vertical Take-Off and Landing Vehicle)は大きなトレンドとなっています。

これらの航空機は飛行のためにプロペラや動力源を動かす必要がありますが、特にプロペラの動きがノイズの大きな原因となっています。

ドローンなどの無人航空機のプロペラは、体積変位と定常的な翼負荷によって発生する多様なセルフノイズと、翼端渦が後続翼に衝突することで発生する相互作用がノイズの発生源と考えられます。

DJI Phantom 3

このケースでは、ドローンプロペラDJI-9450を対象とした空力音響ベンチマークをもとにTCAEによる解析ワークフローについて紹介します。

対象となるDJI-9450は、個人向けとして人気の高いドローン"DJI Phantom 3”に採用されてるものです。

このプロペラは、無響室で測定された実測データが公開されています。(参考文献[1][2][3])

ドローンモデル

DJI-9450のジオメトリ

実機には4つのプロペラが搭載されていますが、無響室での測定は単一のプロペラで行われたため、シミュレーションも単一プロペラでのアプローチとします。

プロペラのブレードは、EPPLER856(0-0.2r/R)とEPPLER63(0.2-1r/R)のプロファイルが組み合わされています。

プロペラの詳細は以下となります; 

DJI-9450の詳細データ

使用する形状データについて

TCAEでのシミュレーションに使用する典型的なフォーマットは、STLサーフェス形式となります。

CFDシミュレーションを行うためには、流体が流れるモデル内部が閉じられた水密性のあるモデルを使用する必要があります。

追加でFEAシミュレーションを行いたいという場合には、閉空間を持つソリッドデータが必要となります。

DJI-9450の厚みとねじれ分布

単一プロペラモデル

Acoustic Analogy(音響学的アナロジー)の理論

Acoustic Analogyは、高い解像度でナビエ・ストークス方程式を解くことに比べて、簡易的に圧力変動を計算する方法を提供します。

これにはいくつかの種類がありますが、TCAE内の音響解析モジュールでは、Ffowcs Williams-HawkingsのFarassat1A定式化を使用することができます。(参考文献[4][5])

Ffowcs Williams-Hawkings方程式の取得手順は以下のようになります。

まず重要なのは、ノイズが発生するサーフェスです。

これは、ノイズの発生源を囲むサーフェスや発生源そのものの表面(プロペラブレードなど)のことを指します。

今回のケースでは、回転するプロペラを囲む領域(円柱形状)を構成するサーフェスとなります。

ここでのボリュームメッシュやサーフェスメッシュは適切に細分化する必要があり、目安としては対応する周波数を捉えるためには、1つの波長で25[cell]です。(周波数1000[Hz]では、少なくとも0.01372[m]のセルサイズが必要)

流体力学方程式を組み替えて、左側の圧力に波動演算子を適用し、右側のソース項にグループ化されたその他の量を作成します。

右側には、モノポール(厚さ)/ダイポール(負荷)のサーフェス項と四重極ボリューム項が含まれています。

TCAEの実装では、乱気流からの非線形ソースである四重極ソースを外しています。

モノポールソース項は、物体が通過する際の流体の変位によって生成されるノイズを表します。

ダイポールソース項は、体表面でのForce分布による非定常運動から生じるノイズをモデル化します。

ここで残りの項を簡略化し、グリーン関数を用いて時間微分を空間微分に書き換えると、最終的に次の式が得られます。

P'_T はthickness noise(翼厚音)となります。

式中の記号は以下となります。

・ f : サーフェス記述関数

・Po: 基準密度

・r:  オブザーバとの距離

・Mr: 放射方向のマッハ数

・M: 局所的なマッハ数

・c: 音速

・dot(v)_n = dot(v)_i * ni , v dot(n) = v_i * dot(n)_i: v_iとn_iは、それぞれ速度と表面法線ベクトルのi番目の成分

・lr = li*ri: riは放射方向ベクトル成分、liは局所的な力の強さ成分、Pijは Lighthill圧縮テンソル成分、uf は表面速度

・下付き文字は遅延時間で評価される量

(上記の説明は簡易的な内容となる、参考文献[6]を参照)

各オブザーバの音圧値の時系列が計算されると(Acoustic Analogyのシミュレーションが完了)、次に周波数領域への変換を行います。

そして、結果データとして各オブザーバの音圧レベル、パワースペクトル密度を取得します。

これらの量は、各ノイズの大きさ(音圧変動周波数)とオブザーバ位置での全体の大きさの指標として活用することができます。

信号処理と呼ばれるこのプロセスは、以下のステップで実行されています。

選択したフィルタ(なし、Butterworth、Chebyshev1、Chebyshev2)により不要な周波数から信号を除去します。

これらは標準的なフィルタで、TCAEではPython SciPyライブラリを実装しています。(参考文献[7])

クリーンアップされた信号は、フーリエ変換により周波数領域に変換されます。

適切な結果を得るためには、データ入力の役割を果たすCFDの結果が十分に正確であることが重要となります。

また、結果が安定している部分のCFDデータのみを使用するように設定することが望ましいです(プロペラの回転初期のデータは含めない)。

また、非定常計算のタイムステップサイズΔtも重要な項目です。

音響解析を行う場合には、CFDのタイムステップサイズを固定して実行することが仕様となります。

これは、フーリエ変換への入力が正しく、サンプリング周波数f s = 1 tが定数 fixed tに対してのみ明確に定義されているためです。

サンプリング周波数f sは、信号処理結果の精度に制限を与えます。

ナイキスト周波数f N = f s 2を超える周波数は正確ではありません。(エイリアシング効果によって歪みが発生)

信号が十分長い場合には、より高度にデータを取り扱うことが可能となります。

ここで用いられるWelch法(参考文献[8])は、選択した数のセグメントで信号を分割することができます。

これは、各セグメントにWindow Functionを適用し、それぞれの結果をフーリエ変換することで、各セグメントを平均化し、最終的な結果を得るという流れとなります。

Window Functionの役割は、フーリエ変換で区間の端の値よりも中心部のデータに注目させることです。

最も重要な理由としては、フーリエ変換では最後の値が最初の値と一致することを暗黙の仮定としているため、これは整数周期の長さを持つ完全な周期信号では正しいですが、現実のアプリケーションでは当てはまりません。

つまりフィルタを使用しない場合、フーリエ変換の入力がおかしくなる可能性があることを意味します。

フーリエ変換が終了すると、考慮した各周波数の値が得られます。

これによって、音圧レベルとパワースペクトル密度を計算されることになります。

どちらも「音の大きさ」を周波数の関数として表現していますが、その定義や意味は同じではありません。

音圧レベル(SPL)の定義: 

Po: 基準音圧、選択されたWindow Functionに依存する正規化係数

パワースペクトル密度は、音圧レベルとは異なり力であるため、ここでは2乗の表記となり、サンプリング周波数で割った密度となります。

CAAシミュレーションへの入力

音響解析用TCAAモジュールへの入力は、ノイズ発生源を取り囲む円柱などのサーフェス上で評価された流体量となります。

これらの量は密度、圧力、速度であり、CFDシミュレーションの実行中に準備されており、今回のケースではプロペラを囲む円柱をソースサーフェスとして使用するため、別途作成して追加をしています。

ここでは、すべてのノイズ発生源がソースサーフェスの内側にあることが重要となります。

このようなケースでの適切な半径は、プロペラ半径の1.2倍程度が推奨されています。

結果精度に影響を与えるもう一つの要因としては、表面解像度があります。

解像度が高ければ高いほど、より正確な結果が得られますが、理論上はローカルのCFDメッシュの粗さよりも高い解像度を持つソースサーフェスは有効ではありません。

解析領域の全体図

回転領域とプロペラ

TCAAソースサーフェスとなる円柱サーフェス

メッシュモデルの作成

CFD解析のメッシュモデルは、メッシュ作成モジュールTMESH内のsnappyHexMeshを用いて作成されています。

メッシュ作成からシミュレーション、後処理までのすべての作業を同一GUI上で実施可能となっています。

ソフトウェアGUI

プロペラハブ付近のメッシュ分布

メッシュはモデル壁面に沿って徐々に細分化されるように、レベルを設定しています。

高精度なCAAシミュレーションを行うためには、メッシュの細分化や境界層レイヤーを適切に追加することが重要となります。

領域全体のメッシュ分布

このケースでは、200万から1200万までのメッシュセル数によるスタディを実行しましたが、結果に大きな影響を与えることはありませんでした。

今回の設定では、最小セルは0.013m以下であるべきという基本理念に則していますが、正確なメッシュ設定は常に検討課題となります。

CFDシミュレーションの設定

CFDシミュレーションは、流体解析モジュールTCFDで設定を行います。

シミュレーションのセットアップは、GUI上のトップメニューから順に行うのみの単純な作業となります。

ソフトウェアGUI(TCFD)

・シミュレーションタイプ: Propeller

・非定常解析

・圧縮性を考慮

・壁面粗さ: なし

・回転数: 6000 [RPM]

・Outlet: 静圧

・乱流モデル: k-omega SST

・壁関数を適用

・乱流強度: 5 [%]

・作業流体: 空気

・基準圧力: 1 [atm]

・動粘性係数: 2.144×10E-5 [PaS]

・基準密度: 1.2 [kg/m3]

・CPU時間(CFD): 25 [core.hours/point/round]

CAAシミュレーションの設定

CAAシミュレーションは、音響解析モジュールTCAAで設定を行います。

TCAAはlibAcoustics OpenFOAMライブラリ(参考文献[10])をベースに構築されており、このライブラリはTCAEに最適な適合となるように大幅な調整が行われています。

オリジナルとの最大の違いは、TCAAはディスクからのデータを使用するため、1つのデータセットをより多くのCAAシミュレーションに繰り返し使用することができることです。(オリジナルのlibAcousticsライブラリは、OpenFOAMのランタイムでテンポラリーデータを使用して動作します。)

ソフトウェアGUI(TCAA)

・音速: 343 [m/s] 

・ソースデータ: TCFD結果

・基準圧力: 1 [atm]

・物質密度: 7800 [kg/m3]

・作動流体: 空気

・Acoustic Analogy: Ffowcs Williams-Hawkings

・基準音圧: 2e-5 [Pa]

・回転数: 6000 [RPM]

シミュレーション結果

シミュレーションはTMESH、TCFD、TCAAの順で実行され、レポート作成までを自動実行します。

すべての解析データは、対応するCSVファイルに保存され、自動レポートに記載するグラフやポスト処理に使用されます。

タイムステップごとの結果グラフ

以下の結果がプロジェクトディレクトリ内に保存されるものとなります。

結果の可視化を行う際はParaViewに搭載された数多くのフィルタ機能を活用して、データの処理を行うことができます。

プロペラ近傍の速度

プロペラ近傍の圧力

乱流運動エネルギー

乱流散逸率

Blade to Bladeの結果

音響解析の結果

音響解析の結果で最も有用なデータは、オブザーバ位置におけるタイムステップごとの音圧、音圧レベル、パワースペクトル密度です。

以下のグラフは、オブザーバ位置を実測時に合わせた、原点をプロペラの中心、回転軸をZ軸として下を向いた位置[0.128, 1.423, 0.779]での結果となります。

実測データとの比較

グラフ"Unfiltered signal"は、Acoustic Analogyによるシミュレーション結果を示しており、オブザーバー位置での音圧推移を表しています。

このケースでは、シミュレーション時間の0.3秒~0.5秒のデータを使っており、プロペラの回転数に換算すると、30〜50回転になります。

グラフ”Power spectral density”と”Sound pressure level”は、フーリエ変換された時間信号から計算された周波数領域の結果です。

パワースペクトル密度については、実際のデータと比較を行っています。

この結果から、今回使用した解析モデルは、ブレード通過周波数(bpf = rotations per second)を最も大きな周波数として正しく認識していることが分かります。

bpfの倍数もピークとして捕らえることができていますが、参考データとの差が存在してい部分も見られます。

これは以下のような理由が考えられます。

・参考データの研究で使用したプロペラと異なる形状を使用した可能性があるため(視覚的には近似しているものの、微妙な差異が異なる結果を引き起こす可能性があります。)

・参考データでは、オブザーバの正確な位置がプロペラ中心ではなく、ドローン中心基準で記載されているため、ドローン寸法を推定した設定としているため

・参考シミュレーション(参考文献[2])では、メッシュ数もタイムステップデータも豊富であり、TCAEケースでは実際のエンジニアリングを考慮した軽量な解析モデルとしているため

これらの点を考慮すると、軽量化されたTCAEシミュレーションは、DJIプロペラ9450の音響挙動を適切に予測していると考えられます。

まとめ

このケースでは、ドローンプロペラを対象としたCFD&CAAシミュレーションについて紹介しました。

CFDの解析手法についてはチューニングの余地があると考えられますが、特別な調整をすることなく、実測データの傾向を掴むことができました。

ターボ機械全般のエンジニアリングシミュレーションツールとして、効果的に活用することが可能であることを、この記事では示しています。

参考文献

[1] Intaratep, N., Alexander, W.N., Devenport, W.J., Grace, S.M., Dropkin, A. (2016) Experimental Study of Quadcopter Acoustics and Performance at Static Thrust Conditions, 22nd AIAA/CEAS Aeroacoustics Conference.

[2] Mankbadi, R.R., Afari, S.O., and Golubev, V.V (2020) Simulations of Broadband Noise of a Small UAS’ Propeller, 10.2514/6.2020-1493.

[3] Zawodny, N.S., Boyd, D.D., and Burley, C.L. (2013) Acoustic Characterization and Prediction of Representative, Small-Scale Rotary-Wing Unmanned Aircraft System Components, Proceedings of the American Helicopter Society 72nd Annual Forum.

[4]  Glegg, S. and Devenport, W. (2017) Aeroacoustics of Low Mach Number Flows, Academic Press

[5] Farassat, F. (2007) Derivation of Formulations 1 and 1A of Farassat

[6] Brentner, K.S., Farassat, F. (2003) Modeling aerodynamically generated sound of helicopter rotors, Progress in Aerospace Sciences, Volume 39, Issues 2–3

[7] Virtanen, P., Gommers, R., Oliphant, T.E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., van der Walt, S. J., Brett, M., Wilson, J., Millman, K.J., Mayorov, Nelson, A.R.J.,  Jones, E., Kern, R., Larson, E., Carey, C.J., Polat, İ., Feng, Y., Moore, E.W., VanderPlas, J., Laxalde, D., Perktold, J., Cimrman, R., Henriksen, I., Quintero, E.A., Harris, Ch.R., Archibald, A.B., Ribeiro, A.H., Pedregosa, F., van Mulbregt, P., (2020) SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods, 17(3), 261-272.

[8] Welch, P.D. (1967) The Use of Fast Fourier Transform for the Estimation of Power Spectra: A Method Based on Time Averaging Over Short, Modified Periodograms, IEEE Transactions on Audio and Electroacoustics, vol. 15, no. 2, pp. 70-73, doi: 10.1109/TAU.1967.1161901.

[9] Heinzel, G., Rüdiger, A. and Schilling R. (2002) Spectrum and spectral density estimation by the Discrete Fourier transform (DFT), including a comprehensive list of window functions and some new flat-top windows.

[10] Epikhin, A., Evdokimov, I., Kraposhin, M., Kalugin, M., Strijhak, S. (2015) Development of a Dynamic Library for Computational Aeroacoustics Applications Using the OpenFOAM Open Source Package, Procedia Computer Science, DOI: 10.1016/j.procs.2015.11.018

[11] Moslem, F., Masdari, M., Fedir, K., Moslem, B., 2022/08/02, SP  – 095441002210913, Experimental investigation into the aerodynamic and aeroacoustic performance of bioinspired small-scale propeller planforms, DO  – 10.1177/09544100221091322, JO  – Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering