標準的van Driest減衰関数型LESフィルタ幅の実装

# 標準的van Driest減衰関数型LESフィルタ幅の実装

本記事は，[オープンCAE学会V&V委員会](http://www.opencae.or.jp/activity/committee/vandv/)の[第2回 OpenFOAMコード検証勉強会](https://opencae-kanto.connpass.com/event/68544/)での「OpenFOAMにおけるvan Driest減衰関数の標準的実装とSmagorinsky LESモデルによるチャネル流れの解析講習」(今野雅)の発表資料(LuaTex形式)の一部をpandocを用いてMarkdown形式に変換し，修正したものです．

## はじめに

OpenFOAM における LES フィルタ幅 vanDriestDeltaは，幾何的な LES フィルタ幅に van Driest 減衰関数を乗ずる一般的な van Driest 減衰関数型 LES フィルタ幅の定義と異なる．ここでは，OpenFOAM-v1612+ に対して，一般的な van Driest 減衰関 数型 LES フィルタ幅の実装を行い，チャネル流のチュートリアルで動作テストする．

## van Driest減衰関数型LESフィルタ幅のコードのコピー

まず， カスタムライブラリのソースコードを置くディレクリを作成する．なお， 場所やディレクトリ名は任意で良い．

```bash
run
cd ..
mkdir src
cd src
```

元のvanDriestDeltaのコードをコピーし， ファイル名を変更する． なお，<span>\#</span>で始まる行は前のコマンドの出力である．

```bash
find $FOAM_SRC -name vanDriestDelta
# /FullPathToSrc/TurbulenceModels/turbulenceModels/LES/LESdeltas/vanDriestDelta

cp -a `find $FOAM_SRC -name vanDriestDelta` vanDriestStandardDelta
cd vanDriestStandardDelta
ls
# vanDriestDelta.C  vanDriestDelta.H

mv vanDriestDelta.C vanDriestStandardDelta.C
mv vanDriestDelta.H vanDriestStandardDelta.H
```

wmakeを用いてソースコードをコンパイルするにはMakeディレクトリが必要となるが，元のvanDriestDeltaのディレクトリの上流のディレクトリのどこかにvanDriestDelta用のMakeディレクトリがあるので， これをコピーする．

```bash
cp -a $FOAM_SRC/TurbulenceModels/turbulenceModels/Make .
```

## gitを用いたコード履歴管理

これからソースコードやMake内の設定ファイルを修正していくので，変更点を記録するために， gitを用いて管理する．

```bash
git init
# Initialized empty Git repository in /FullPath/src/vanDriestStandardDelta/.git/

git add .
git commit -m "Initial commit."
# master (root-commit) f6455fc] Initial commit.
#  4 files changed, 363 insertions(+), 0 deletions(-)
#  create mode 100644 Make/files
#  create mode 100644 Make/options
#  create mode 100644 vanDriestStandardDelta.C
#  create mode 100644 vanDriestStandardDelta.H
```

ここで， mオプションで指定するコミットメッセージは任意である． また，OpenFOAMではコンパイルする際に， lnIncludeディレクトリを作成し，その下にソースコードがインクルードするファイルへのシンボリックを作成する．さらに， Makeディレクトリの下に，linux64GccDPInt32Optのようなオブジェクトファイルなどを置くディレクトリを作成する．これらのディレクトリ下のファイルは自動で生成されるものなので，以下のような内容で.gitignoreを作成して，gitが無視するように設定しておくと便利である．

```git:.gitignore
lnInclude/
Make/*/
```

.gitignore についてもgitで管理するように設定する．

```bash
git add .gitignore
git commit -m "Add .gitignore."
```

## コンパイルするファイルのリストと生成ライブラリ名の定義

Make/files にはコンパイルするファイル名リストと， 生成ライブラリ名を定義するので，例えば以下のように変更する． ここで， $(FOAM_USER_LIBBIN)は，ユーザーのカスタムライブラリを置くディレクトリを示す環境変数である．なお， ライブラリ名は任意である．

```C++:Make/files
vanDriestStandardDelta.C

//LIB = $(FOAM_LIBBIN)/libturbulenceModels
LIB = $(FOAM_USER_LIBBIN)/libvanDriestStandardDelta
```

## クラス名の変更

ソースコード中のクラス名をvanDriestDeltaからvanDriestStandardDeltaに変更する．また，ランタイム型名をvanDriestStandardDeltaからvanDriestStandardに変更する．

```bash
sed -i s/vanDriest/vanDriestStandard/g vanDriestStandardDelta.H
sed -i s/vanDriest/vanDriestStandard/g vanDriestStandardDelta.C
```

gitで変更箇所を確認しておく．

```bash
git diff
```

## インクルードパスの設定

本格的に実装を変更する前に， コンパイル可能か調べる．

```bash
wmake
# wmakeLnInclude: linking include files to ./lnInclude
# Making dependency list for source file vanDriestStandardDelta.C
# could not open file LESdelta.H for source file vanDriestStandardDelta.C due to No such file or directory
```

LESdelta.Hが無いというエラーになるので， このファイルを捜す．

```bash
find $FOAM_SRC -name LESdelta.H
# /FullPath/src/TurbulenceModels/turbulenceModels/LES/LESdeltas/LESdelta/LESdelta.H
# /FullPath/src/TurbulenceModels/turbulenceModels/lnInclude/LESdelta.H
```

2つのファイルが見付かっているが， 1つめがファイルの実体であり，2つめのlnInclude以下のファイルは，実体ファイルにシンボリックリンクしたものである．コンパイル時にLESdelta.Hが見付かるよう， Make/filesのEXE_INCに，LESdelta.Hのインクルードパスを追加すれば良いのであるが， 以下に示すように，通常はlnIncludeのディレクトリを指定する． これは， lnIncludeディレクトリには，他のディレクトリに散在するインクルードファイルを，シンボリックリンクという形で集約しているので，実体ファイルのディレクトリを指定するより，インクルードパスの指定が少なくて済むためである． また，デバッグコードを有効にするために， -DDEBUG を追加して， DEBUGを定義にしておく．

```C++:Make/files
EXE_INC = \
    -DDEBUG \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude

LIB_LIBS = \
    -lfiniteVolume \
    -lmeshTools
```

## コード変更前のコンパイル・実行テスト

コンパイル可能かテストする．

```bash
wmake
```

実装を修正する前に， vanDriestStandardDeltaがソルバーで使用可能かテストするため，
channel395のチュートリアル・ケースでテストする． まず，channel395チュートリアルをコピーする．

```bash
cp -a $FOAM_TUTORIALS/incompressible/pimpleFoam/LES/channel395 .
pushd channel395
```

次に， system/controlDictのlibでカスタムライブリlibvanDriestStandardDelta.soを指定する． なお， .soは共有ライブラリ(shared object)の拡張子である． また，1時間ステップで終了して結果を出力するように， endTimeとwriteInterval を変更する．

```C++:system/controlDict
libs ("libvanDriestStandardDelta.so");

application     pimpleFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

//endTime         1000;                                                                       
endTime         0.2;

deltaT          0.2;

writeControl    timeStep;

//writeInterval   200;                                                                        
writeInterval   1;
```

さらに， constant/turbulencePropertiesにおいて， LESModelにSmagorinskyを指定し， deltaにvanDriestStandardを指定する．

```C++:constant/turbulenceProperties
LES
{
    //    LESModel        WALE;                                                               
    LESModel        Smagorinsky;

    turbulence      on;

    printCoeffs     on;

    //    delta           cubeRootVol;                                                        
    delta           vanDriestStandard;
```

また， vanDriestStandardの係数設定用のvanDriestStandardCoeffsを用意する．

```C++:constant/turbulenceProperties
    //    vanDriestCoeffs                                                                     
    vanDriestStandardCoeffs
    {
```

ソルバーpimpleFoamが実行できれば， vanDriestStandardのカスタムライブラリが有効になっている事が確認できる．

```bash
blockMesh
pimpleFoam
```

## LESフィルタ幅の実装変更

vanDriestStandard のソースコードのディレクトリに復帰する．

```bash
popd
```

先の資料で示した通り，一般的にはvan Driestの減衰関数型フィルター幅*Δ*<sub>*v*</sub>は以下で定義される．

*Δ*<sub>*v*</sub> = [ 1 − exp (−*y*<sup>+</sup>/*A*<sup>+</sup>)) ] *Δ<sub>geo</sub>*

 ここで， *y*<sup>+</sup>は無次元壁座標，*A*<sup>+</sup>は係数でOpenFOAMのデフォルト値は26，*Δ<sub>geo</sub>*は幾何的グリッドフィルタ幅であり，通常は格子体積の1/3乗である． 上式を実装するには，vanDriestDeltaのソースを例えば以下のように修正すれば良い． なお，#ifdef DEBUGと#endif間はデバッグ用コードであるため， 必須でないが，このコードはグリッドフィルタ幅deltaの最小値と最大値を表示し，場のファイルを出力すると共に， 格子毎の各項を表示する．

```C++:vanDriestStandardDelta.C
    scalar cutOff = wallPointYPlus::yPlusCutOff;
//    wallPointYPlus::yPlusCutOff = 500;
    wallPointYPlus::yPlusCutOff = HUGE;
    wallDistData<wallPointYPlus> y(mesh, ystar);
    wallPointYPlus::yPlusCutOff = cutOff;

/*
    delta_ = min
    (
        static_cast<const volScalarField&>(geometricDelta_()),
        (kappa_/Cdelta_)*((scalar(1) + SMALL) - exp(-y/ystar/Aplus_))*y
    );
*/

    const volScalarField& geometricDelta = geometricDelta_();
    delta_ = (scalar(1) - exp(-y/ystar/Aplus_)) * geometricDelta;

#ifdef DEBUG
    Info << "vanDriestDelta min/max = " << min(delta_).value() << ", "
         << max(delta_).value() << endl;
    delta_.write();
    forAll(delta_, deltaI)
    {
        Info << "C: " << mesh.C()[deltaI]
            << " y: " << y[deltaI]
            << " ystar: " << ystar[deltaI]
            << " fv: " << scalar(1)-exp(-y[deltaI]/ystar[deltaI]/Aplus_)
            << " deltaGeo: " << geometricDelta[deltaI]
            << " delta: " << delta_[deltaI]
            << endl;
    }
#endif
```

また， 変数Cdelta_は不要となったので， この変数に関連する部分を修正する．なお， 以下では修正がわかりやすいようにコメントアウトしているが，最終的にはコメントアウトしている行は消去して構わない．

```C++:vanDriestStandardDelta.H
    //        scalar Cdelta_;
```

```C++:vanDriestStandardDelta.H
/*                                                                                                                      
    Cdelta_                                                                                                             
    (                                                                                                                   
        dict.subDict(type() + "Coeffs").lookupOrDefault<scalar>("Cdelta", 0.158)                                        
    ),                                                                                                                  
*/
(snip)
   //    coeffsDict.readIfPresent<scalar>("Cdelta", Cdelta_);  
```

## 実行テスト

channel395チュートリアルでソルバーを実行する．

```bash
pushd channel395
pimpleFoam > log.pimpleFoam
more log.pimpleFoam
# vanDriestDelta min/max = 1.04361e-29, 0.088154
# C: (0.05 0.00480001 0.0333333) y: 0.00480001 ystar: 0.00120188 fv: 0.14239 deltaGeo: 0.04 delta: 0.0056956
# C: (0.15 0.00480001 0.0333333) y: 0.00480001 ystar: 0.00120847 f: 0.141671 deltaGeo: 0.04 delta: 0.00566684
```

格子中心の座標ベクトルCの*y*成分と， 壁座標のyが一致することを確認する．
また， van Driestの減衰関数*f*<sub>*v*</sub> = 1 − exp (−*y*/(*y*<sup>\*</sup>*A*<sup>+</sup>))，およびLESグリッドフィルタ幅*Δ* = *f*<sub>*v*</sub>*Δ*<sub>*g**e**o*</sub>が正しく計算されているか確認する．なお，channel395 チュートリアルでの*A*+は26である．

##LESグリッドフィルタ幅の分布確認

paraviewを用いて， LESグリッドフィルタ幅delta の分布を確認したのが， 以下の図である．

![グリッドフィルタ幅の分布](https://qiita-image-store.s3.amazonaws.com/0/131980/06bc4ada-50ff-163d-f5fe-a9ac30b43f2b.png)
図：グリッドフィルタ幅の分布

デバッグコードを除いたコンパイル
--------------------------------

LESグリッドフィルタ幅deltaの算出結果が正しい事を確認したら，デバッグコードを除いて再コンパイルするため，ソースコードのディレクトリに復帰する．

```bash
popd
```

Make/optionsのEXE_INCを以下のように変更する．

```C++:Make/options
EXE_INC = \
    /* -DDEBUG */ \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude
```

再コンパイルする．

```bash
wmake
```

最後に，gitを用いてコミットする．

```bash
git add Make/ vanDriestStandardDelta.*
git commit -m "Implement standard van Driest damping function type."
```

## まとめ
van Driest 減衰関数型 LES フィルタ幅として一般的な実装を行なった．OpenFOAMにおいて，乱流モデル関係ライブラリのカスタマイズは容易な部類に入るので，カスタマイズ初心者の方は是非試してみてください．