標準的van Driest減衰関数型LESフィルタ幅の実装
本記事は，オープンCAE学会V&V委員会の第2回 OpenFOAMコード検証勉強会での「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フィルタ幅のコードのコピー
まず， カスタムライブラリのソースコードを置くディレクリを作成する．なお， 場所やディレクトリ名は任意で良い．
run
cd ..
mkdir src
cd src
元のvanDriestDeltaのコードをコピーし， ファイル名を変更する． なお，#で始まる行は前のコマンドの出力である．
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ディレクトリがあるので， これをコピーする．
cp -a $FOAM_SRC/TurbulenceModels/turbulenceModels/Make .
gitを用いたコード履歴管理
これからソースコードやMake内の設定ファイルを修正していくので，変更点を記録するために， gitを用いて管理する．
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が無視するように設定しておくと便利である．
lnInclude/
Make/*/
.gitignore についてもgitで管理するように設定する．
git add .gitignore
git commit -m "Add .gitignore."
コンパイルするファイルのリストと生成ライブラリ名の定義
Make/files にはコンパイルするファイル名リストと， 生成ライブラリ名を定義するので，例えば以下のように変更する． ここで， $(FOAM_USER_LIBBIN)は，ユーザーのカスタムライブラリを置くディレクトリを示す環境変数である．なお， ライブラリ名は任意である．
vanDriestStandardDelta.C
//LIB = $(FOAM_LIBBIN)/libturbulenceModels
LIB = $(FOAM_USER_LIBBIN)/libvanDriestStandardDelta
クラス名の変更
ソースコード中のクラス名をvanDriestDeltaからvanDriestStandardDeltaに変更する．また，ランタイム型名をvanDriestStandardDeltaからvanDriestStandardに変更する．
sed -i s/vanDriest/vanDriestStandard/g vanDriestStandardDelta.H
sed -i s/vanDriest/vanDriestStandard/g vanDriestStandardDelta.C
gitで変更箇所を確認しておく．
git diff
インクルードパスの設定
本格的に実装を変更する前に， コンパイル可能か調べる．
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が無いというエラーになるので， このファイルを捜す．
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を定義にしておく．
EXE_INC = \
-DDEBUG \
-I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude
LIB_LIBS = \
-lfiniteVolume \
-lmeshTools
コード変更前のコンパイル・実行テスト
コンパイル可能かテストする．
wmake
実装を修正する前に， vanDriestStandardDeltaがソルバーで使用可能かテストするため，
channel395のチュートリアル・ケースでテストする． まず，channel395チュートリアルをコピーする．
cp -a $FOAM_TUTORIALS/incompressible/pimpleFoam/LES/channel395 .
pushd channel395
次に， system/controlDictのlibでカスタムライブリlibvanDriestStandardDelta.soを指定する． なお， .soは共有ライブラリ(shared object)の拡張子である． また，1時間ステップで終了して結果を出力するように， endTimeとwriteInterval を変更する．
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を指定する．
LES
{
// LESModel WALE;
LESModel Smagorinsky;
turbulence on;
printCoeffs on;
// delta cubeRootVol;
delta vanDriestStandard;
また， vanDriestStandardの係数設定用のvanDriestStandardCoeffsを用意する．
// vanDriestCoeffs
vanDriestStandardCoeffs
{
ソルバーpimpleFoamが実行できれば， vanDriestStandardのカスタムライブラリが有効になっている事が確認できる．
blockMesh
pimpleFoam
LESフィルタ幅の実装変更
vanDriestStandard のソースコードのディレクトリに復帰する．
popd
先の資料で示した通り，一般的にはvan Driestの減衰関数型フィルター幅Δvは以下で定義される．
Δv = [ 1 − exp (−y+/A+)) ] Δgeo
ここで， y+は無次元壁座標，A+は係数でOpenFOAMのデフォルト値は26，Δgeoは幾何的グリッドフィルタ幅であり，通常は格子体積の1/3乗である． 上式を実装するには，vanDriestDeltaのソースを例えば以下のように修正すれば良い． なお，#ifdef DEBUGと#endif間はデバッグ用コードであるため， 必須でないが，このコードはグリッドフィルタ幅deltaの最小値と最大値を表示し，場のファイルを出力すると共に， 格子毎の各項を表示する．
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_は不要となったので， この変数に関連する部分を修正する．なお， 以下では修正がわかりやすいようにコメントアウトしているが，最終的にはコメントアウトしている行は消去して構わない．
// scalar Cdelta_;
/*
Cdelta_
(
dict.subDict(type() + "Coeffs").lookupOrDefault<scalar>("Cdelta", 0.158)
),
*/
(snip)
// coeffsDict.readIfPresent<scalar>("Cdelta", Cdelta_);
実行テスト
channel395チュートリアルでソルバーを実行する．
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の減衰関数fv = 1 − exp (−y/(y*A+))，およびLESグリッドフィルタ幅Δ = fvΔgeoが正しく計算されているか確認する．なお，channel395 チュートリアルでのA+は26である．
LESグリッドフィルタ幅の分布確認
paraviewを用いて， LESグリッドフィルタ幅delta の分布を確認したのが， 以下の図である．
図：グリッドフィルタ幅の分布
デバッグコードを除いたコンパイル
LESグリッドフィルタ幅deltaの算出結果が正しい事を確認したら，デバッグコードを除いて再コンパイルするため，ソースコードのディレクトリに復帰する．
popd
Make/optionsのEXE_INCを以下のように変更する．
EXE_INC = \
/* -DDEBUG */ \
-I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude
再コンパイルする．
wmake
最後に，gitを用いてコミットする．
git add Make/ vanDriestStandardDelta.*
git commit -m "Implement standard van Driest damping function type."
まとめ
van Driest 減衰関数型 LES フィルタ幅として一般的な実装を行なった．OpenFOAMにおいて，乱流モデル関係ライブラリのカスタマイズは容易な部類に入るので，カスタマイズ初心者の方は是非試してみてください．