SENAC Vol. 56, No. （2 2023. 4）― 13 ―[ 共同研究成果 ]低レイノルズ数領域における円弧翼周りの非定常流体解析佐々木大輔
大阪公立大学髙瀬拓海，岡本正人
金沢工業大学

本研究では，直交格子法の一つである Building-Cube Method を用いて，レイノルズ数 1,000と 10,000 における円弧翼周りの非定常流体解析を行った．解析の結果，それぞれのレイノルズ数において，風洞実験で得られた空力係数と同様の傾向を得ることができ，レイノルズ数効果を確認することができた．更に，迎角を上げた際の流れ場の変化について，可視化図より明らかにした．

1.

緒言

近年，様々な宇宙探査計画が各国で計画されている．その一つに火星探査が挙げられる．火星は，過去にも多くの探査計画が実施されてきた．例えば，無人探査ローバー「キュリオシティ」やランダー型探査機「インサイト」等，様々な手法で数多くの探査が行われてきた．その結果，火星に対する理解が進む一方，火星で起きる事象の普遍的な理解には程遠い状態にある．したがって，現在も火星は研究の対象とされている．例えば，宇宙航空研究開発機構（JAXA）は，地表面を詳細にかつ広範囲に探査可能な「火星探査飛行機」を提唱している．
火星探査飛行機を実現するためには，火星の大気環境が地球と異なることから，火星大気環境に適した飛行機を設計する必要がある．例えば，火星大気の空気密度は地球の 1/100 程度しかないため，揚力を得るのが困難となる．その大気条件で飛行する場合，レイノルズ数は 104～105 のオーダーとなることから，低レイノルズ数領域において空力性能に優れた翼型の設計が必要となる[1]．レイノルズ数とは，慣性力と粘性力の比で表される無次元量であり，火星飛行機のレイノルズ数域は粘性力の強い流れと言える．地球上において，このレイノルズ数域を飛行するものは鳥や昆虫であり，これらの翼型を模倣することで火星大気環境下でも飛行可能な飛行機が開発できると考えられる．
岡本らは火星大気環境を模した低圧風洞を用いることで，種々の翼型に対して風洞試験を実施しており，その知見から低レイノルズ数領域においては一般的に薄翼が適していることがわかっている．また，岡本らはコガネムシに着目し，様々なキャンバーを持つ円弧翼に対して風洞試験を行うことで低レイノルズ数領域における空力特性を明らかにしている[2]．円弧翼の空力性能にはレイノルズ数依存性があることが風洞実験で確認されているため，異なるレイノルズ数領域での解析を実施することで，Computational Fluid Dynamics (CFD)においても同様の傾向が得られるかを本研究では確認する．そして空力性能のレイノルズ数依存性と流れ場の関係性について考察を行う．今回，風洞試験で空力性能の差が顕著であった Re=1,000 と Re=10,000 について解析を実施する．一般的に Re=1,000 では流れが 2 次元的であることが知られていることから 2 次元ソルバーを用いて解析を実施する．また Re=10,000 においては流れに 3 次元性が強く表れることが知られているため，3 次元ソルバーを用いて大規模非定常流体解析を実施する．

― 14 ―

2.

SENAC Vol. 56, No. （2 2023. 4）

解析手法

解析には非圧縮性流体解析ソルバーの Building-Cube Method（BCM）を用いる[3]．BCM は，ブロック型直交格子に分類される CFD 手法である[4,5]．BCM の特徴として，複雑形状に対しても格子生成が容易であること，計算領域を Cube と呼ばれる領域に分割し，Cube 内を等間隔・同数の Cell で分割することにより，均等な負荷分散を達成できること，また空間高次精度解法を採用できること等が挙げられる．図 1 に物体近傍の Cube および Cell を示す．
本研究で計算に用いた支配方程式は非圧縮性 Navier-Stokes 方程式と連続の式である．
ここで，時間積分法は 2 次精度 Adams-Bashforth 法，空間離散化手法は 3 次精度風上差分法（対流項）および 2 次精度中心差分法（拡散項）を使用する．今回，境界埋め込み法等は使用せず，物体形状は階段状表現として解析を行う．

図 1 Cube と Cell3.

解析対象及び解析条件

本研究で使用する解析モデルはコガネムシの翼型を模した薄い円弧翼である．解析モデルを図2 に示す．解析モデルの翼弦長 c を基準として，キャンバーを 3％c，翼の板厚を 0.01c としている．解析領域は翼弦長 c を基準に x，y 方向に 48 倍とする．3 次元解析の際は，z 方向に翼を延伸し， 1.5 倍の解析領域とする．流入条件は層流とし，レイノルズ数 1,000 では 2 次元解析，レイノルズ数 10,000 では 3 次元解析を行う．なお，3 次元解析の際，計算領域の両端には周期境界条件を適用し，2 次元性を確保している．

図 2 円弧翼

低レイノルズ数領域における円弧翼周りの非定常流体解析― 15 ―

本研究では異なるレイノルズ数での解析を実施するに当たって，2 次元ソルバーと 3 次元ソルバーを用いた．どちらも Cube と呼ばれる立方体(2 次元では正方形)の中に等間隔の Cell を配置している．今回は翼型に迎角を付けた解析を複数回行うため，それぞれの迎角に適した格子を作成する必要がある．本解析では表 1 に示すような 4 種類の格子を使用し，低迎角と高迎角で異なる格子を使用する．2 次元解析では約 100 万格子，3 次元解析では最大で約 1 億格子となる．図 3に 2 次元における円弧翼周りの低迎角用 Cube 配置を示す．本研究では，AOBA-B を使用した並列計算を実施することで，約 1 億格子規模の非定常流体解析を行い，平均場を取得する．
表 1 格子情報Type of Airfoil
Type of mesh
Number of cubes
Number of cells in a cube
Number of cells
Min. cell size
Outer size

2D modelLow AOAHigh AOA
3,547
5,149
162
162
908,032
1,318,144
0.0007324
0.0007324
48×4848×48

3D modelLow AOAHigh AOA
10,53015,808163
183
43,130,880
92,192,256
0.002930
0.002604
48×48×1.5 48×48×1.5

図 3 円弧翼周りの Cube 配置

4.

解析結果

揚力係数，抗力係数，モーメント係数の解析値と実験値を比較したグラフを図 4，図 5，図 6にそれぞれ示す．図 4 より，実験では Re=1,000 に比べて Re=10,000 の方が揚力傾斜が大きいことが分かる．解析でも同様の傾向が得られている．図 5 より，抗力係数の実験値は低迎角において，レイノルズ数依存性が顕著に確認できる．本解析で得た解析値においても同様のレイノルズ数依存性が確認でき，低迎角では Re=1,000 に比べて Re=10,000 の方が小さい．
また，
図 6 より，Re=1,000においてモーメント係数が広い範囲で一定となっているのに対して，Re=10,000 では高迎角でモーメント係数が減少する傾向にあることが確認できる．
以下では，揚力係数の迎角変化に着目し，流れ場と対比させることで，流れ場の理解を試みる．
2 次元解析を実施した Re=1,000 において，迎角 6 度までは揚力係数がよく一致しているが，それ以降では実験値との差が広がる傾向にある．図 7 に迎角 3 度と迎角 9 度における可視化図を示す．迎角 3 度では，定常的な流れ場となっており，前縁から後縁まで流れが剥がれることなく翼

― 16 ―

SENAC Vol. 56, No. （2 2023. 4）型に沿って流れている．一方，迎角 9 度では剥離した流れが循環渦を形成して後縁付近に停滞している．更に，迎角 12 度の場合，翼型下面から上面に回り込もうとする流れが後縁で停滞していた渦と干渉することで後流域に非対称な渦を発生させ，非定常流れとなる．これらの高迎角域では流れが 2 次元的ではないと考えられ，揚力係数が大きく実験値と乖離している原因と考えられる．
また，3 次元解析を実施した Re=10,000 において，揚力係数の解析値は迎角 3 度以下で実験値と近い値を取っており，迎角 6 度以上では実験値と離れていることがわかる．実験値では迎角 6度付近で揚力傾斜が減少し始めており，揚力係数の分布が非線形となっている．実験値同様，解析値でも同様に揚力係数の減少傾向が確認できる．図 8(a)に 0.75 度，図 8(b)に 3 度の迎角における瞬間的な 3 次元渦構造を示す．これらの低迎角では 2 次元的な流れ場となっており，翼型後縁から 2 次元的な非対称渦が放出されていることが確認できる．一方，高迎角における瞬間的な3 次元渦構造を図 8(c)(d)に示す．図 8(c)は迎角 6 度，図 8(d)は迎角 9 度の可視化図である．どちらの迎角においても非周期的な 3 次元流れ場となっていることが確認できる．揚力係数が非線形になっている迎角と一致していることから，非線形性は 3 次元流れによる影響が大きいと考えられる．
以上のことから，本研究で得られた解析値においても実験値と同様のレイノルズ数依存性が確認できた．そのため，3%キャンバーを持つ円弧翼の性能は，Re=1,000 に比べて Re=10,000 の方が引き出せることがわかった．

低レイノルズ数領域における円弧翼周りの非定常流体解析

図 4 揚力係数

図 5 抗力係数

図 6 モーメント係数― 17 ―

― 18 ―

SENAC Vol. 56, No. （2 2023. 4）(a) 迎角 3 度(b) 迎角 6 度(c) 迎角 9 度(d) 迎角 12 度
図 7 Re=1,000 における瞬間場の速度コンター図

低レイノルズ数領域における円弧翼周りの非定常流体解析(a) 迎角 0.75 度(b) 迎角 3 度(c) 迎角 6 度
図 8 Re=10,000 における瞬間場の速度コンター図（続く）― 19 ―

― 20 ―

SENAC Vol. 56, No. （2 2023. 4）(d) 迎角 9 度
図 8 Re=10,000 における瞬間場の速度コンター図

5.

結言

本研究では 2 次元／3 次元非圧縮性 BCM ソルバーを用いた数値計算により，低レイノルズ数領域における 3%キャンバーを持つ円弧翼周りの流れ場の解析を実施した．その結果，異なるレイノルズ数領域(Re=1,000 と Re=10,000)においても，実験値と同様の空力係数を得られることを示し，円弧翼の空力特性がレイノルズ数に依存することを明らかにした．また，流れ場の可視化を通して，円弧翼周りの空力現象を明らかにした．

謝辞
本研究は，東北大学サイバーサイエンスセンターの共同研究制度により，スーパーコンピ
ュータ AOBA-B を利用することで実現することができた．また，研究にあたっては同センター関係各位に有益なご指導とご協力をいただいた．
参考文献[1] H. Hidaka and M. Okamoto, “An Experimental Study of Triangular Airfoils for Mars
Airplane”, Aerospace Technology Japan, Vol. 12, No. ists29, Pk_21-Pk_27, 2014.
https://doi.org/10.2322/tastj.12.Pk_21[2] M. Okamoto and K. EBINA, “Effectiveness of Large-Camber Circular Arc Airfoil at Very
Low Reynolds Numbers”, Transactions of the Japan Society for Aeronautical and Space
Sciences, Vol.59, No.5, pp295-304, 2016. https://doi.org/10.2322/tjsass.59.295[3] 坂井玲太郎，大林茂，松尾裕一，中橋和博，“Building-Cube Method を用いた実用的な大規模乱流解析手法の構築”，JSASS-2013-2116-A，第 45 回流体力学講演会/航空宇宙数値シミュレーション技術シンポジウム 2013 講演論文集，船堀，2013.
[4] K. Nakahashi and L. S. Kim, “Building-Cube Method for Large-Scale, High Resolution
Flow Computations”, AIAA Paper 2004-0423, 2004.
[5] D. Sasaki and K. Nakahashi, “Rapid Large-scale Cartesian Meshing for AerodynamicComputations”, ICAS Paper 2014-0842, Proceedings of 29th Congress of the International
Council of the Aeronautical Sciences (ICAS2014), St. Petersburg, Russia, 2014.