更新履歴

2025/01/22: 初版
2025/02/20: 計算環境とOpenFOAM設定を修正
2025/04/19: boundaryの設定を修正

はじめに

本稿ではこちらのチュートリアルを参考に、"NACA 0012"翼型モデルの非圧縮性乱流モデルを非定常で計算する備忘を綴る。

環境

AWS Parallel Cluster: 3.9.20
ヘッドノードインスタンスタイプ：r6g.medium
計算ノードインスタンスタイプ: c7gn.xlarge
OpenFOAM: OpenFOAM-v2306
スケジューラ: Slurm

大体これで12hくらい計算時間がかかります。こちらのEC2インスタンスタイプ料金表を基に一回あたりの計算コストを算出：

(0.0608_{[\$]}+0.3148_{[\$]})\times12_{[h]}=4.5072_{[\$]}

よって、1計算あたり約700円弱となる。

チュートリアルケースの取得

以下のチュートリアルケースを作業ディレクトリにコピーしてくる。

cp -r $FOAM_TUTORIALS/incompressible/simpleFoam/airFoil2D .

treeコマンドでファイルとフォルダを確認すると以下のようになる。

% tree .
.
├── 0.orig
│   ├── U
│   ├── nuTilda
│   ├── nut
│   └── p
├── Allclean
├── Allrun
├── constant
│   ├── polyMesh.orig
│   │   ├── boundary.gz
│   │   ├── cells.gz
│   │   ├── faces.gz
│   │   ├── neighbour.gz
│   │   ├── owner.gz
│   │   └── points.gz
│   ├── transportProperties
│   └── turbulenceProperties
└── system
    ├── controlDict
    ├── decomposeParDict
    ├── fvSchemes
    └── fvSolution

4 directories, 18 files

主な各ファイルの役割は次の通り。

0.orig/
場の変数の設定ファイル（速度U、圧力pなど）が含まれるディレクトリ
Allrun
1. 境界条件や初期値の設定を行う0を0.origから、そしてconstant/polymeshをconstant/polymesh.origから参照
2. 並列計算を実行
constant/
メッシュや物性値、乱流モデルなどの設定が含まれるディレクトリ
1. polyMesh/
  メッシュ情報を含むディレクトリ。チュートリアル内でplot3dToFoamを用いたメッシュ分割結果が格納
2. transportProperties/, turbulenceProperties/
  物性値の設定ファイル
system/
計算の制御の設定ファイルなどが含まれるディレクトリ

Paraviewによるメッシュの確認

作業ディレクトリ直下にcase.foamなどを作成して、ParaViewでメッシュを確認する。

メッシュの置き換え

メッシュは上のOpenFOAMチュートリアルで用意されているものではなく、NASAが提供するNACA 0012翼型に特化した計算格子があるので、それを用いる。メッシュの詳細についてはこちらを参考

本稿でもそれに倣って、(2-D) 225 x 65 (129 points on airfoil surface)を用いる。作業ディレクトリ以下で以下のコマンドを実行

cd constant/
wget https://turbmodels.larc.nasa.gov/NACA0012_grids/n0012_225-65.p2dfmt
cd ..
plot3dToFoam constant/n0012_225-65.p2dfmt -2D 1 -noBlank

constant/polyMesh/が生成されているので、改めてメッシュを可視化して確認してみる。

【遠距離ビュー】

【近距離ビュー】

パッチの分割

上記の操作により作成されたconstant/polyMesh/boundaryを見るとdefaultFacesのみになっているので、以下のコマンドでパッチを自動分割する。

autoPatch -overwrite 30

こうするとauto_*というパッチ名が自動で割り当てられているので、分かり易く次のように書き換える。面倒であるが、ParaViewを使って1パッチずつ対応を確認した。

constant/polyMesh/boundary

5
(
    back
    {
        type            empty;
        physicalType    empty;
        nFaces          14336;
        startFace       28432;
    }

    front
    {
        type            empty;
        physicalType    empty;
        nFaces          14336;
        startFace       42768;
    }

    outlet
    {
        type            patch;
        physicalType    outlet;
        nFaces          128;
        startFace       57104;
    }

    aerofoil
    {
        type            wall;
        physicalType    wall;
        nFaces          128;
        startFace       57232;
    }
    
    inlet
    {
        type            patch;
        physicalType    inlet;
        nFaces          224;
        startFace       57360;
    }
)

パッチが定義できたところで以下のコマンドでpolyMesh.origを置き換える。

rm -rf constant/polyMesh.orig/*
cp -r constant/polyMesh/* constant/polyMesh.orig/

空力係数の計算

forceCoeffs関数を用いてCl値/Cd値など空力係数を計算する。

system/forceCoeffs

forces
{
    type            forceCoeffs;
    libs            ( forces );
    writeControl    timeStep;
    writeInterval   1;
    patches         ( aerofoil );
    rho             rhoInf;
    log             true;
    rhoInf          1;
    liftDir         ( 0 1 0 );
    dragDir         ( 1 0 0 );
    CofR            ( 0.25 0 0 ); // 回転中心座標。翼先頭から1/4地点
    pitchAxis       ( 0 0 1 );
    magUInf         51.4815;
    lRef            1; // 翼の長さ (NASA 2DグリッドでX方向=1)
    Aref            1; // =翼の長さ×スパン長（NASA 2DグリッドでZ方向=1）
}

PISO法による非定常計算

OpenFOAMには非定常解析ソルバーとして、PISO法やPIMPLE法などが用意されている。本章ではPISO法アルゴリズム概要について簡単に記載して置く。

ナビエ・ストークス（Navier-Stokes/NS）方程式

一般に、流速がマッハ数0.3以下（密度変化が5%以下）であれば非圧縮性流れとして考える。後節で記述するが、今回の翼形状周りの流れシミュレーションには流速$=51.4815_{[m/sec]}\simeq0.15_{Ma}$を与えるので、非圧縮性方程式を考える。

流速を$\boldsymbol{u}$、圧力を$\boldsymbol{p}$、流体密度を$\rho$、粘性応力テンソルを$\boldsymbol{\tau}$として、ナビエ-ストークス方程式は次式で表される：

\cfrac{\partial\boldsymbol{u}}{\partial t}
+
\left(
    \boldsymbol{u}\cdot\nabla
\right)
\boldsymbol{u}
=
-\cfrac{1}{\rho}\nabla\boldsymbol{p}
+
\nabla\cdot\boldsymbol{\tau}

ここで、粘性応力テンソル$\boldsymbol{\tau}$はせん断粘性係数を$\mu$、粘性係数を$\lambda$として以下のように分解される：

\boldsymbol{\tau}
=\mu
\lbrack
    \nabla\boldsymbol{u}
    +(\nabla\boldsymbol{u})^\mathsf{T}
\rbrack
+\lambda
(
    \nabla\cdot\boldsymbol{u}
)
\boldsymbol{1}

さらに多くの場合、体積粘性率$\chi=\lambda+(2/3)\mu$はストークス近似のもとで0とされる。
故に動粘性係数を$\nu:=\mu/\rho$として元の式に代入すると：

\cfrac{\partial\boldsymbol{u}}{\partial t}
+
\left(
    \boldsymbol{u}\cdot\nabla
\right)
\boldsymbol{u}
=
-\cfrac{1}{\rho}\nabla\boldsymbol{p}
+
\nabla\cdot
\lbrack
    \nu\lbrace
        \nabla\boldsymbol{u}
        +(\nabla\boldsymbol{u})^\mathsf{T}
    \rbrace
\rbrack
-
\nabla
\left(
    \frac{2}{3}\nu\nabla\boldsymbol{u}
\right)

PISO法による圧力-流速連成

OpenFOAMのPISO法では以下の計算手順により圧力と流速を連成させる：

予測速度場の計算
ナビエ-ストークス方程式を空間離散化して次のような式を構築する：
```
A_P\boldsymbol{u}_P + \sum A_N\boldsymbol{u}_N = -\nabla p
```
ここで、$A_P$および$A_N$はそれぞれ注目セルと隣接セルの係数である。
これを次の形に分離する：
```
A\boldsymbol{u} = \boldsymbol{H} - \nabla p
```
圧力項を一旦無視して予測速度場$\hat{\boldsymbol{u}}$を計算する：
```
A\hat{\boldsymbol{u}} = \boldsymbol{H}
```
OpenFOAMにおけるmomentumPredictorオプションは、この仮速度場の計算を実行するスイッチである。
圧力修正方程式の構築
非圧縮性条件$\nabla\cdot\boldsymbol{u}=0$を満たすため、圧力補正方程式を構築する。
仮速度場$\hat{\boldsymbol{u}}$を基に、次のような式が導かれる：
```
\nabla\cdot\left(\frac{1}{\rho}\nabla\boldsymbol{p}'\right) = \nabla\cdot\hat{\boldsymbol{u}}
```
ここで、$\boldsymbol{p}'$は圧力補正項を表す。
圧力場の更新
圧力補正項$\boldsymbol{p}'$を用いて圧力場を更新する：
```
\boldsymbol{p}^{(n+1)} = \boldsymbol{p}^{(n)} + \boldsymbol{p}'
```
速度場の修正
修正された圧力場$\boldsymbol{p}'$を用いて速度場を次のように修正する：
```
A\boldsymbol{u}^{(n+1)} = \boldsymbol{H} - \nabla\boldsymbol{p}^{(n+1)}
```
この修正はPISO法の反復処理中にnCorrectorsオプションによって指定された回数だけ行われる。
残差の監視
各計算ステップの残差（例えば速度、圧力の収束状況）はsolverInfo.datファイルに出力される。このファイルには、各反復ステップの残差値が記録され、収束基準を確認するために使用する。
反復処理
圧力補正と速度修正を複数回繰り返すことで、連続の式$\nabla\cdot\boldsymbol{u}=0$をより高い精度で満たすように調整する。必要に応じて、圧力補正方程式を再構築し、精度を向上させる。

非定常計算への変更

コピーした定常計算チュートリアルとNASA提供のグリッドを元に非定常計算を行う変更を行う。

`controlDict`

前述の通り、今回は代数方程式ソルバーにpisoFoamを用いる。

system/controlDict

application     pisoFoam;
startFrom       startTime;
startTime       0;
stopAt          endTime;
endTime         5.0;
deltaT          1e-4;
writeControl    runTime;
writeInterval   1e-4;
purgeWrite      2001;
writeFormat     ascii;
writePrecision  10;
writeCompression off;
timeFormat      general;
timePrecision   6;
runTimeModifiable true;
functions
{
    #include "forceCoeffs"
    #includeFunc solverInfo(U,p,k,omega)   
}

`deltaT`と`writePrecision`について

非定常シミュレーション結果として、システムの振動が観測できる値にする。例えば、以下のようにグラフ可視化しても、データを直接表示しても同じ値を取るケースなどが生じうる。

【postProcessing/forces/0/coefficient.datの描画結果】

【postProcessing/forces/0/coefficient.datのデータ中身】

このように、系の振動を観測できない場合はdeltaTを小さく、writePrecisionを大きくするのを試してみると良さそうである。

`endTime`について

後述する仰角alphaDegにもよるが、coefficient.datを観察してみると大5秒もあれば収束することが分かる。あるいは、残差solverInfo.datの収束を見て設定しても良い。

`fvScheme`

以下の点だけ修正を加える。backwardで後方二次精度となる。

ddtSchemes
{
    default         backward;
}

`fvSolution`

Field	Linear Solver	Smoother	Tolerance (rel)
U	Smooth solvers	Gauss Seidel Smoother	0.01
p	GAMG Solver	Gauss Seidel Smoother	0.01
nuTilda	Smooth solvers	Gauss Seidel Smoother	0.01

数値スキーム設定`system/fvSolution`

system/fvSolution

solvers
{
    p
    {
        solver          GAMG;
        smoother        GaussSeidel;
        tolerance       1e-09;
        relTol          1e-02;
    }
    pFinal
    {
        $p;
        relTol          0;
    }
    U
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        nSweeps         2;
        relTol          1e-02;
        tolerance       1e-09;
    }
    nuTilda
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        nSweeps         2;
        relTol          1e-02;
        tolerance       1e-09;
    }
}

PISO
{
    momentumPredictor        yes;
    nCorrectors              3;
    nNonOrthogonalCorrectors 0;
}

relaxationFactors
{
    fields
    {
        p               0.3;
    }
    equations
    {
        U               0.7;
        nuTilda         0.7;
    }
}

`tolerance`について

計算結果として出力されるpostProcessing/solverInfo(U,p,k,omega)/0/solverInfo.datを見てみると、以下のようなグラフを可視化することができる。

【postProcessing/solverInfo(U,p,k,omega)/0/solverInfo.datの描画結果】

まだ十分に収束しているとは言えないが、設定した$1e-8$付近のオーダーで残差が推移していることが分かる。
加えて、この時のcoefficient.datを併せて見てみると

【postProcessing/forces/0/coefficient.datの描画結果】

このように観測データである空力係数Cl, Cdが一定時間経過後に収束していることが読み取れる。（前述した振動はグラフでは見えないレベルで生じている）一方で、torelanceの値が大き過ぎる、或いは設定しないと以下のように空力係数の値が収束（厳密には微振動）せず、大きな振れ幅のまま発散してしまう。

【空力係数が発散するケース】

実際、このようなケースは後述するCd, Clプロットで比較データと大きく差異が生じてしまったため、良いとは推測し難い。このように、空力係数と残差の時間推移を観察しながら適切なtorelanceの決定することで、PISO法による代数方程式計算の打ち切り（精度）を調整する。

乱流モデル`constant/turbulenceProperties`の設定

2次元シミュレーションでのLESモデルの使用に関して、LESモデルは通常、3次元非定常計算で用いられるが、2次元問題や定常計算でも選択可能となっている。ところが、2次元でのLESモデルの適用は物理的に適切でないため、今回はSpalart-Allmarasモデルを採用する。参考はこちら。

constant/turbulenceProperties

simulationType          RAS;

RAS
{
    RASModel            SpalartAllmaras;
    turbulence          on;
    printCoeffs         on;
}

場の設定

チュートリアルを再現するために、場の設定を変更する。

`Allrun`

今後の拡張として仰角alphaDegを変えてCFDを流したいので、foamDictionaryコマンドを利用して場や制御の変数を書き換えるようにする。

Allrun

#!/bin/sh
cd "${0%/*}" || exit                                # Run from this directory
. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions

#------------------------------------------------------------------------------
# Variables

parallel=true # flag to enable computations in parallel mode
alphaDeg=0 # Angle in degrees (adjust as needed)
speed=51.4815 # Speed (adjust as needed)

# Calculate cos(α) and sin(α)
alphaRad=$(echo "$alphaDeg * 3.141592653589793238462643 / 180" | bc -l)
cosAlpha=$(echo "c($alphaRad)" | bc -l)
sinAlpha=$(echo "s($alphaRad)" | bc -l)

#------------------------------------------------------------------------------
# Set initial speed field

# Modify `0.orig/U` internalField
foamDictionary 0.orig/U \
    -entry internalField \
    -set "uniform ($(
        echo "$speed * $cosAlpha" | bc -l
    ) $(
        echo "$speed * $sinAlpha" | bc -l
    ) 0)"

# Modify `0.orig/U` boundaryField inlet
foamDictionary 0.orig/U \
    -entry boundaryField.inlet.freestreamValue \
    -set "uniform ($(
        echo "$speed * $cosAlpha" | bc -l
    ) $(
        echo "$speed * $sinAlpha" | bc -l
    ) 0)"

# Modify `0.orig/U` boundaryField outlet
foamDictionary 0.orig/U \
    -entry boundaryField.outlet.freestreamValue \
    -set "uniform ($(
        echo "$speed * $cosAlpha" | bc -l
    ) $(
        echo "$speed * $sinAlpha" | bc -l
    ) 0)"

#------------------------------------------------------------------------------
# Direction setting of force coefficient (lift and drag)

# Modify `system/forceCoeffs` liftDir
foamDictionary system/forceCoeffs \
    -entry forces.liftDir \
    -set "($(
        echo "-1 * $sinAlpha" | bc -l
    ) $(
        echo "$cosAlpha" | bc -l
    ) 0)"

# Modify `system/forceCoeffs` dragDir
foamDictionary system/forceCoeffs \
    -entry forces.dragDir \
    -set "($(
        echo "$cosAlpha" | bc -l
    ) $(
        echo "$sinAlpha" | bc -l
    ) 0)"

#------------------------------------------------------------------------------
# Initialization

restore0Dir
[ -d constant/polyMesh ] || cp -rf constant/polyMesh.orig constant/polyMesh

#------------------------------------------------------------------------------
#  Running Parallel Computations

if [ "$parallel" = true ]
then
    runApplication decomposePar
    runParallel $(getApplication)
else
    runApplication $(getApplication)
fi

`0.orig`

速度場`U`

Patch	Condition	Value [$m\cdot s^{-1}$]
Inlet	freestreamVelocity	51.4815
Outlet	freestreamVelocity	51.4815
Sides (y-dir)	empty	-
Aerofoil	fixedValue	(0.0, 0.0, 0.0)

上記のAllrunによりinternalFieldフィールドは書き換えられるが、安全のために51.4815で初期化しておく。

0.orig/U

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (51.4815 0 0);

boundaryField
{
    inlet
    {
        type            freestreamVelocity;
        freestreamValue $internalField;
    }
    outlet
    {
        type            freestreamVelocity;
        freestreamValue $internalField;
    }
    front
    {
        type            empty;
    }
    back
    {
        type            empty;
    }
    aerofoil
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
}

圧力場`p`

Patch	Condition	Value [$m^2\cdot s^{-2}$]
Inlet	freestreamPressure	0
Outlet	freestreamPressure	0
Sides (y-dir)	empty	-
Aerofoil	zeroGradient	-

0.orig/p

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    inlet
    {
        type            freestreamPressure;
        freestreamValue $internalField;
    }
    outlet
    {
        type            freestreamPressure;
        freestreamValue $internalField;
    }
    front
    {
        type            empty;
    }
    back
    {
        type            empty;
    }
    aerofoil
    {
        type            zeroGradient;
    }
}

乱流動粘性係数`nut`

Patch	Condition	Value [$m^2\cdot s^{-1}$]
Inlet	freestream	$8.58\times10^{-6}$
Outlet	freestream	$8.58\times10^{-6}$
Sides (y-dir)	empty	-
Aerofoil	fixedValue	0

0.orig/nut

dimensions      [0 2 -1 0 0 0 0];

internalField   uniform 8.56e-6;

boundaryField
{
    inlet
    {
        type            freestream;
        freestreamValue $internalField;
    }
    outlet
    {
        type            freestream;
        freestreamValue $internalField;
    }
    front
    {
        type            empty;
    }
    back
    {
        type            empty;
    }
    aerofoil
    {
        type            fixedValue;
        value           uniform 0.0;
    }
}

乱流動粘性`nuTilda`

Patch	Condition	Value [$m^2\cdot s^{-1}$]
Inlet	freestream	$3.432\times10^{-5}$
Outlet	freestream	$3.432\times10^{-5}$
Sides (y-dir)	empty	-
Aerofoil	fixedValue	0

0.orig/nuTilda

dimensions      [0 2 -1 0 0 0 0];

internalField   uniform 3.432e-5;

boundaryField
{
    inlet
    {
        type            freestream;
        freestreamValue $internalField;
    }
    outlet
    {
        type            freestream;
        freestreamValue $internalField;
    }
    front
    {
        type            empty;
    }
    back
    {
        type            empty;
    }
    aerofoil
    {
        type            fixedValue;
        value           uniform 0.0;
    }
}

シミュレーション結果

本シミュレーション結果の検証には、NACA 0012翼型の揚力係数Clおよび抗力係数Cdと迎角（Angle of Attack）alphaDegとの関係を示したデータセットを使用した。このデータセットは、以下のような実験的/数値的な結果に基づいく。

実験データ: Ladsonらによる実験データ

Ladsonらの実験データは、NASA Langley 研究センターで取得されたもので、風洞試験による高精度なClとCdの測定値を含む。
数値データ: 数値シミュレーション (CFL3D ソフトウェア)
NASAが開発したCFL3Dソフトウェアを使用して計算された数値データ。この数値解は、実験データとの整合性を評価するための基準として使用される。

本シミュレーションでは、実験データとの比較を行い、揚力係数および抗力係数の計算結果がこれらの既知データとどの程度一致しているかを検証した。特に、迎角の増加に伴う非線形挙動や失速時の特性を検証の焦点に置く。

1. 揚力係数：`Cl`値

2. 抗力係数：`Cd`値

このように揚力係数Cl値、抗力係数Cd値いずれのケースにおいても実験データの特徴を良く再現出来ているように見える。仰角alphaDeg=16付近を頭打ちにそれ以上の仰角だと、翼状後流の剥離が顕著となり、揚力/抗力共にドラスティックに変化する。

流れ場の可視化

ParaViewを用いて翼状まわりの流れ場を可視化する。手順は以下の通り：

processor*/ディレクトリがある同じ階層にcase.foamを作成し、ParaViewでそれを読み込む。
画像内①にあるCase Typeを"Decomposed Case"に設定
画像内②にあるScalar Sizeを64bitに設定
Applyを押下
画像内④で描画する場を選択

仰角`alphaDeg`ごとの流速場分布

1. `alphaDeg`=0

2. `alphaDeg`=10

3. `alphaDeg`=15

剥離現象の観察

ParaViewの視覚化手法を"Surface LIC (Line Integral Convolution)"を選択して、場の流れを直感的に視覚化してみる。alphaDeg=19と仰角が大き過ぎて揚力が落ち込んでいるケースで試してみる。

おわりに

いい感じに再現できた。OpenFOAMは情報共有がネット上にもあまりなく、苦戦したが本稿を作成するに際して一番役に立ったのはこの30秒のYouTube動画だった。意外とYouTubeで調べて見るのもアリだなと思った。