【OpenFOAM v2006】低Re数型非線形RANSモデルの実装

OpenFOAM v2006で低Re数型非線形RANSモデルと壁関数を実装する

計算に使ったファイルはGithubにアップしてある
mtkbirdman.com/OpenFOAM/

はじめに

参考にしたサイトは以下の通り

OpenFOAMカスタマイズ入門:レイノルズ平均型乱流モデルのカスタマイズに挑戦 – Qiita
Turbulence Models in OpenFOAM

今回実装する低Re数型非線形RANSモデルはこれ

Ken-ichi ABE, A Hybrid LES/RANS Approach Using an Anisotropy-Resolving Algebraic Turbulence Model, International Journal of Heat and Fluid Flow, Vol. 26, pp. 204-222, 2005.
Ken-ichi ABE, Yong-Jun JANG and Michael A. LESCHZINER, An Investigation of Wall-Anisotropy Expressions and Length-Scale Equations forNon-Linear Eddy-Viscosity Models, International Journal of Heat and Fluid Flow, Vol. 24, pp. 181-198, 2003.

2003年の論文のモデルを改良した2005年の論文のモデル(AJL2005)を実装する.記述は2003年の論文の方が詳しい

乱流モデルAJL2005についての説明は面倒なので紙面の都合上割愛させていただく

C++については次のサイトを見て完全に理解した
C++入門

OpenFOAMのテンソル計算については次のサイトが詳しい
Tensor Operations in OpenFOAM _ CFD WITH A MISSION

OpenFOAMにおけるRANSについてのドキュメントはこちら
OpenFOAM User Guide Reynolds Averaged Simulation (RAS)

乱流モデル実装までの手順

乱流モデルの実装は以下の手順で行う

  1. 実装したい乱流モデルを見つける
  2. OpenFOAMに実装されている乱流モデルの中から実装したいモデルに近いものを見つける (このモデルをベースに実装する)
  3. ベースとなる乱流モデルのファイルをコピーして乱流モデルの名前を書き換える
  4. 名前を書き換えた乱流モデルをコンパイルして計算してみる
  5. 乱流モデルのファイルの中身を書き換える
  6. 書き換えた乱流モデルをコンパイルして計算する

手順1はすでに完了しているので,手順2から説明していく

ベースとなる乱流モデルの決定

AJL2005は低Re数型非線形RANSモデルである

OpenFOAMで実装されているNon-linear eddy viscosity modelsには

の2つがあるが,低Re数型であるのはLienCubicKEだけなのでこれをベースとする

ちなみになんとLienCubicKEの参考論文はインターネット上では手に入らない (筆者調べ2020年9月時点)

どうしても読みたい人は近くの図書館で探してみよう
CiNii 図書 – Engineering turbulence modelling and experiments 3 _ proceedings of the Third International Symposium on Engineering Turbulence Modelling and Measurements, Heraklion-Crete, Greece, 27-29.html

ベースとなる乱流モデルのコピー・名前の書き換え

ベースとなる乱流モデルをコピーしてくる

LienCubicKEに関するファイルは$WM_PROJECT_DIR/src/TurbulenceModels/incompressible/turbulentTransportModels/RASというディレクトリにあるので,これを$WM_PROJECT_USER_DIR/src/myTurbulenceModelsにコピーする

myTurbulenceModelsにあたるディレクトリ名はわかりやすければ何でもいい

>> mkdir -p $WM_PROJECT_USER_DIR/src/myTurbulenceModels
>> cp -r $WM_PROJECT_DIR/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/LienCubicKE $WM_PROJECT_USER_DIR/src/myTurbulenceModels

ちなみに$WM_PROJECT_DIR$WM_PROJECT_USER_DIR$WM_PROJECT_DIR/etc/bashrcで定義されている変数で,echoを使って中身を確認できる

>> echo $WM_PROJECT_DIR
/opt/OpenFOAM/OpenFOAM-v2006
>> echo $WM_PROJECT_USER_DIR
/home/mtkbirdman/OpenFOAM/mtkbirdman-v2006

コピーしてきたすべてのファイルでLienCubicKEAJL2005に書き換える

>> cd $WM_PROJECT_USER_DIR/src/myTurbulenceModels
>> mv LienCubicKE AJL2005
>> cd AJL2005
>> mv LienCubicKE.C AJL2005.C
>> mv LienCubicKE.H AJL2005.H
>> sed -i -e "s/LienCubicKE/AJL2005/g" AJL2005.*

乱流モデルのコンパイル(1回目)

名前を書き換えたAJL2005(実際はLienCubicKE)をコンパイルするために$WM_PROJECT_USER_DIR/src/myTurbulenceModelsに2つのファイルMake/filesMake/optionsを作成する

>> cd $WM_PROJECT_USER_DIR/src/myTurbulenceModels
>> mkdir Make
>> cd Make
>> touch files
>> touch options

それぞれのファイルに以下を書きこむ

Make/files
AJL2005/AJL2005.C
LIB = $(FOAM_USER_LIBBIN)/libmyTurbulenceModels
Make/options
EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude \
    -I$(LIB_SRC)/fvOptions/lnInclude \
    -I$(LIB_SRC)/transportModels \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \

LIB_LIBS = \
    -lincompressibleTurbulenceModels \
    -lturbulenceModels \
    -lfvOptions \
    -lfiniteVolume \
    -lmeshTools

ここまででディレクトリ構成は以下のようになる

myTurbulenceModels
    ├── AJL2005
    │   ├── AJL2005.C
    │   └── AJL2005.H
    └── Make
        ├── files
        └── options

AJL2005をコンパイルする

>> cd $WM_PROJECT_USER_DIR/src/myTurbulenceModels
>> wmake libso
wmake libso (myTurbulenceModels)
Make/options:9:12: warning: backslash-newline at end of file
 LIB_LIBS = \

    ln: ./lnInclude
Making dependency list for source file AJL2005.C
g++ -std=c++11 -m64 -pthread -DOPENFOAM=2006 -DWM_DP -DWM_LABEL_SIZE=32 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -Wno-attributes -Wno-unknown-pragmas  -O3  -DNoRepository -ftemplate-depth-100 -I/opt/OpenFOAM/OpenFOAM-v2006/src/finiteVolume/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/meshTools/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/fvOptions/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/transportModels -I/opt/OpenFOAM/OpenFOAM-v2006/src/TurbulenceModels/turbulenceModels/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/TurbulenceModels/incompressible/lnInclude  -iquote. -IlnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/OpenFOAM/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/OSspecific/POSIX/lnInclude   -fPIC -c AJL2005/AJL2005.C -o Make/linux64Gcc63DPInt32Opt/AJL2005/AJL2005.o
g++ -std=c++11 -m64 -pthread -DOPENFOAM=2006 -DWM_DP -DWM_LABEL_SIZE=32 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -Wno-attributes -Wno-unknown-pragmas  -O3  -DNoRepository -ftemplate-depth-100 -I/opt/OpenFOAM/OpenFOAM-v2006/src/finiteVolume/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/meshTools/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/fvOptions/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/transportModels -I/opt/OpenFOAM/OpenFOAM-v2006/src/TurbulenceModels/turbulenceModels/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/TurbulenceModels/incompressible/lnInclude  -iquote. -IlnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/OpenFOAM/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/OSspecific/POSIX/lnInclude   -fPIC -shared -Xlinker --add-needed -Xlinker --no-as-needed  Make/linux64Gcc63DPInt32Opt/AJL2005/AJL2005.o -L/opt/OpenFOAM/OpenFOAM-v2006/platforms/linux64Gcc63DPInt32Opt/lib \
    -lincompressibleTurbulenceModels -lturbulenceModels -lfvOptions -lfiniteVolume -lmeshTools  -o /home/mtkbirdman/OpenFOAM/mtkbirdman-v2006/platforms/linux64Gcc63DPInt32Opt/lib/libmyTurbulenceModels.so

プログラムに問題があればここで大量のエラーを吐き出してくれる

※ちなみに,Make/filesMake/optionsの編集をVS codeでやりたいときはWSLのhomeディレクトリをエクスプローラーで開く必要がある.そのときのhomeディレクトリのパスを以下に示す

C:\Users\<ユーザー名>\AppData\Local\Packages\CanonicalGroupLimited.Ubuntu20.04onWindows_79rhkp1fndgsc\LocalState\rootfs\home\<WSLのユーザー名>

テストケースの準備

テストケースとして/tutorials/incompressible/pimpleFoam/LES/channel395$WM_PROJECT_USER_DIRにコピーしてくる

>> cp -r $WM_PROJECT_DIR/tutorials/incompressible/pimpleFoam/LES/channel395 $WM_PROJECT_USER_DIR

もともとあった0を削除して0/nuTildaを0/epsilonに変更し,0.origをコピーする

>> cd $WM_PROJECT_USER_DIR/channel395
>> rm -r 0
>> mv 0.orig/nuTilda 0.orig/epsilon
>> cp -r 0.orig 0

いくつかのファイルに修正を加える

0/epsilon
/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2006                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      epsilon;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField   uniform 1;

boundaryField
{
    bottomWall
    {
        type            epsilonWallFunction;
        value           $internalField;
    }
    topWall
    {
        type            epsilonWallFunction;
        value           $internalField;
    }
    sides1_half0
    {
        type            cyclic;
    }
    sides2_half0
    {
        type            cyclic;
    }
    inout1_half0
    {
        type            cyclic;
    }
    inout2_half0
    {
        type            cyclic;
    }
    sides2_half1
    {
        type            cyclic;
    }
    sides1_half1
    {
        type            cyclic;
    }
    inout1_half1
    {
        type            cyclic;
    }
    inout2_half1
    {
        type            cyclic;
    }
}

// ************************************************************************* //

内部領域の値は1を指定し,壁関数はepsilonWallFunctionを用いる

0/k
/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2006                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "1";
    object      k;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField   uniform 1;

boundaryField
{
    bottomWall
    {
        type            fixedValue;
        value           uniform 1e-32;
    }
    topWall
    {
        type            fixedValue;
        value           uniform 1e-32;
    }
    sides1_half0
    {
        type            cyclic;
    }
    sides2_half0
    {
        type            cyclic;
    }
    inout1_half0
    {
        type            cyclic;
    }
    inout2_half0
    {
        type            cyclic;
    }
    sides2_half1
    {
        type            cyclic;
    }
    sides1_half1
    {
        type            cyclic;
    }
    inout1_half1
    {
        type            cyclic;
    }
    inout2_half1
    {
        type            cyclic;
    }
}


// ************************************************************************* //

内部領域の値は1を指定し,壁面上ではできるだけ小さい値(1e-32)を指定する

本来は滑りなし条件のため壁面上でk=0であるが,0を指定するとプログラム上のどこかで0割りが発生してしまう

0/nut
/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2006                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "1";
    object      nut;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField   uniform 0;

boundaryField
{
    bottomWall
    {
        type            nutLowReWallFunction;
        value           uniform 0;
    }
    topWall
    {
        type            nutLowReWallFunction;
        value           uniform 0;
    }
    sides1_half0
    {
        type            cyclic;
    }
    sides2_half0
    {
        type            cyclic;
    }
    inout1_half0
    {
        type            cyclic;
    }
    inout2_half0
    {
        type            cyclic;
    }
    sides2_half1
    {
        type            cyclic;
    }
    sides1_half1
    {
        type            cyclic;
    }
    inout1_half1
    {
        type            cyclic;
    }
    inout2_half1
    {
        type            cyclic;
    }
}

// ************************************************************************* //

壁関数はnutLowReWallFunctionを使う

constant/turbulenceProperties
/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2006                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      turbulenceProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

simulationType  RAS;

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

// ************************************************************************* //

RASモデルのAJL2005を指定する

system/controlDict
/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2006                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application     simpleFoam;
startFrom       startTime;
startTime       0;
stopAt          endTime;
endTime         5000;
deltaT          1;
writeControl    timeStep;
writeInterval   5000;
purgeWrite      0;
writeFormat     ascii;
writePrecision  6;
writeCompression off;
timeFormat      general;
timePrecision   6;
runTimeModifiable true;

functions
{
    #includeFunc    residuals
    #includeFunc    "sample"
    #includeFunc    "turbulenceFields"
    #includeFunc    "wallShearStress"
}

libs
(
    "libmyTurbulenceModels.so"
);

// ************************************************************************* //

ALJ2005のライブラリを読み込むために,libs("libmyTurbulenceModels.so");を加えている

functionsについてはGithubのコードを参照
mtkbirdman.com/OpenFOAM/channel395/

system/fvSchemes
/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2006                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default         steadyState;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default         none;
    div(phi,U)      bounded Gauss linearUpwind grad(U);
    div(phi,k)      bounded Gauss limitedLinear 1;
    div(phi,epsilon) bounded Gauss limitedLinear 1;
    div((nuEff*dev2(T(grad(U))))) Gauss linear;
    div(nonlinearStress) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear uncorrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         uncorrected;
}

wallDist
{
    method          meshWave;
}

// ************************************************************************* //

div(nonlinearStress) Gauss linear;を加えるなど

system/fvSolution
/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2006                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    p
    {
        solver          GAMG;
        tolerance       1e-10;
        relTol          0.1;
        smoother        DICGaussSeidel;
        nPreSweeps      0;
        nPostSweeps     2;
        cacheAgglomeration on;
        agglomerator    faceAreaPair;
        nCellsInCoarsestLevel 10;
        mergeLevels     1;
    }

    "(U|k|epsilon|R|nuTilda)"
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-10;
        relTol          0.1;
    }
    omega
    {
        solver          PBiCGStab;
        preconditioner  DILU;
        tolerance       1e-10;
        relTol          0.1;
    }
}

SIMPLE
{
    nNonOrthogonalCorrectors 0;
    consistent          yes;
    pRefCell        1001;
    pRefValue       0;
}

relaxationFactors
{
    equations
    {
        U               0.9;
        k               0.8;
        epsilon         0.7;
        nuTilda         0.8;
        omega           0.8;
    }
}

cache
{
    grad(U);
}

// ************************************************************************* //

PIMPLE→SIMPLE(あんまりわかっていない)

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

runApplication blockMesh
runApplication decomposePar -cellDist
runParallel $(getApplication)
runApplication reconstructPar
runApplication postChannel

#------------------------------------------------------------------------------

計算の実行(1回目)

計算を実行して,コンパイルがうまくできていることを確認する

ちなみに計算はこんな感じ

結果が正確かどうかは置いといて,とりあえず正常に計算は回っているので良しとする

プログラムの修正

プログラムの全体はGithubにあるコードを確認してもらうとして,ここでは重要な部分のみ説明する
mtkbirdman.com/OpenFOAM/src/myTurbulenceModels/AJL2005/

AJL2005.C

メインの部分がここ

// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //

tmp<volScalarField> AJL2005::Rt() const
{
    return sqr(k_)/(nu()*epsilon_);
}

tmp<volScalarField> AJL2005::fw(const dimensioned<scalar> xi) const
{
    return exp(-sqr((pow(nu()*epsilon_,0.25)*y_/nu())/xi));
}

tmp<volScalarField> AJL2005::fMu(const volScalarField& Rt,const volScalarField& fw) const
{
    return (scalar(1) + (35.0/pow(Rt,0.75))*exp(-pow(Rt/30.0,0.75)))
            *(scalar(1) - fw);
}

tmp<volScalarField> AJL2005::f2(const volScalarField& Rt,const volScalarField& fw) const //feps
{
    return (scalar(1) - 0.3*exp(-sqr(Rt/6.5)))
            *(scalar(1) - fw);
}

void AJL2005::correctNut()
{
    correctNonlinearStress(fvc::grad(U_));
}


void AJL2005::correctNonlinearStress(const volTensorField& gradU)
{
    nut_ = Cmu_*fMu(Rt(),fw(26.0))*sqr(k_)/epsilon_;
    nut_.correctBoundaryConditions();

    volSymmTensorField S(symm(gradU));
    volTensorField W(-skew(gradU));

    volScalarField tau(nut_/k_);
    
    volScalarField fB(scalar(1) + Ceta_*(CD_*tau)*(mag(W) - mag(S)));
    volScalarField CB(1.0/(scalar(1) + (22.0/3.0)*sqr(CD_*tau)*magSqr(W) + (2.0/3.0)*sqr(CD_*tau)*(magSqr(W) - magSqr(S))*fB));
    
    volScalarField fr1((magSqr(W) - magSqr(S))/max(magSqr(W) + magSqr(S),dimensionedScalar("minWS",dimensionSet(0,0,-2,0,0,0,0),SMALL)));
    volScalarField fr2(magSqr(S)/max(magSqr(W) + magSqr(S),dimensionedScalar("minWS",dimensionSet(0,0,-2,0,0,0,0),SMALL)));
    
    volScalarField fs1(fr1*fr2*(0.15*Ceta_)*sqr(CD_*tau)*(magSqr(W) - magSqr(S)));
    volScalarField fs2(-fr1*fr2*(scalar(1)+(0.07*Ceta_)*(CD_*tau)*(mag(W) - mag(S))));
    
    volScalarField nStar((nu()*epsilon_,0.25)*y_/nu());
    volVectorField N(fvc::grad(nStar));
    volVectorField d(N/mag(N));
    volSymmTensorField dd(symm(d*d));

    volScalarField taud1((scalar(1) - fw(15.0))*(k_/epsilon_) + fw(15.0)*(1.0)*sqrt(nu()/epsilon_));
    volSymmTensorField bw1(
        fw(26.0)*(
            - (1.0)*(1.0/2.0)*dev(dd)
            + (scalar(1) - sqr(fr1))*sqr(taud1)*(
                ((0.25*0.5)/(scalar(1) + 0.5*sqr(taud1)*mag(S)*mag(W)))*twoSymm(S&W)
                + ((1.5*0.5)/(scalar(1) + 0.5*sqr(taud1)*magSqr(S)))*dev(innerSqr(S))
            )
        )
    );
    
    volScalarField taud2((scalar(1)-fw(15.0))*(k_/epsilon_) + fw(15.0)*(3.0)*sqrt(nu()/epsilon_));
    volSymmTensorField bw2(
        fw(26.0)*(
            - (0.0)*(1.0/2.0)*dev(dd)
            + (scalar(1) - sqr(fr1))*sqr(taud2)*(
                (((13.0/30.0)*1.0)/(scalar(1) + 1.0*sqr(taud2)*mag(S)*mag(W)))*twoSymm(S&W)
                + ((0.6*1.0)/(scalar(1) + 1.0*sqr(taud2)*magSqr(S)))*dev(innerSqr(S))
            )
        )
    );

    volScalarField ftau(exp(-3.0*(scalar(1) - fw(26.0))*(k_/epsilon_)*mag(S)));
    volSymmTensorField bw(ftau*bw1 + (scalar(1) - ftau)*bw2);
    
    nonlinearStress_ = 
        2.0*nut_*S // Cancel linear term
        +2.0*(k_/CD_)*(
            - CB*(tau*CD_)*S //b1
            + (scalar(1) - fw(26.0))*(
                2.0*CB*sqr(tau*CD_)*(twoSymm(S&W) + dev(innerSqr(S))) //b2
                - CB*fs1*(tau*CD_)*S + 2.0*fs2*sqr(tau*CD_)*dev(innerSqr(S)) //bs
            )
            + CD_*bw //bw
        );
}

注意すべき点は以下の通り

  1. \({\bf \Omega}={\bf -}{\sf skew(gradU)}\)
  2. nonlinearStress_volSymmTensorField
  3. nonlinearStress_はレイノルズ応力の非線形項のみ

\({\bf \Omega}={\bf -}{\sf skew(gradU)}\)についてはTensor Operations in OpenFOAMに書いてある

レイノルズ応力テンソルReddyViscosity.Cで次のように定義されている

template<class BasicTurbulenceModel>
Foam::tmp<Foam::volSymmTensorField>
Foam::eddyViscosity<BasicTurbulenceModel>::R() const
{
    tmp<volScalarField> tk(k());

    // Get list of patchField type names from k
    wordList patchFieldTypes(tk().boundaryField().types());

    (中略)

    return tmp<volSymmTensorField>
    (
        new volSymmTensorField
        (
            IOobject
            (
                IOobject::groupName("R", this->alphaRhoPhi_.group()),
                this->runTime_.timeName(),
                this->mesh_,
                IOobject::NO_READ,
                IOobject::NO_WRITE,
                false
            ),
            ((2.0/3.0)*I)*tk() - (nut_)*dev(twoSymm(fvc::grad(this->U_))),
            patchFieldTypes
        )
    );
}

この時点でRは次の形になっている

\begin{eqnarray}
{\bf R_{ij}}={\overline {u’_{i} u’_{j}}}=\frac{2}{3}k\delta_{ij}-\nu_{t}(\frac{\partial U_{i}}{\partial x_{j}}+\frac{\partial U_{j}}{\partial x_{i}})
\end{eqnarray}

そして,nonlinearStress_nonlinearEddyViscosity.Cで次のように定義され,レイノルズ応力テンソルRに加えられている

    nonlinearStress_
    (
        IOobject
        (
            IOobject::groupName("nonlinearStress", alphaRhoPhi.group()),
            this->runTime_.timeName(),
            this->mesh_
        ),
        this->mesh_,
        dimensionedSymmTensor(sqr(dimVelocity), Zero)
    )
template<class BasicTurbulenceModel>
Foam::tmp<Foam::volSymmTensorField>
Foam::nonlinearEddyViscosity<BasicTurbulenceModel>::R() const
{
    tmp<volSymmTensorField> tR
    (
        eddyViscosity<BasicTurbulenceModel>::R()
    );
    tR.ref() += nonlinearStress_;
    return tR;
}

このようにnonlinearStress_volSymmTensorFieldでレイノルズ応力の非線形項のみなので,nonlinearStress_ = ...の右辺のテンソルはすべてvolSymmTensorFieldである必要があるし,AJL2005では線形項をキャンセルする必要がある

輸送方程式の部分はほとんど変更していない

void AJL2005::correct()
{
    if (!turbulence_)
    {
        return;
    }

    nonlinearEddyViscosity<incompressible::RASModel>::correct();

    tmp<volTensorField> tgradU = fvc::grad(U_);
    const volTensorField& gradU = tgradU();

    volScalarField G
    (
        GName(),
        (nut_*twoSymm(gradU) - nonlinearStress_) && gradU
    );


    // Update epsilon and G at the wall
    epsilon_.boundaryFieldRef().updateCoeffs();

    // Dissipation equation
    tmp<fvScalarMatrix> epsEqn
    (
        fvm::ddt(epsilon_)
      + fvm::div(phi_, epsilon_)
      - fvm::laplacian(DepsilonEff(), epsilon_)
      ==
        Ceps1_*G*epsilon_/k_
      - fvm::Sp(Ceps2_*f2(Rt(),fw(3.3))*epsilon_/k_, epsilon_)
    );

    epsEqn.ref().relax();
    epsEqn.ref().boundaryManipulate(epsilon_.boundaryFieldRef());
    solve(epsEqn);
    bound(epsilon_, epsilonMin_);


    // Turbulent kinetic energy equation
    tmp<fvScalarMatrix> kEqn
    (
        fvm::ddt(k_)
      + fvm::div(phi_, k_)
      - fvm::laplacian(DkEff(), k_)
      ==
        G
      - fvm::Sp(epsilon_/k_, k_)
    );

    kEqn.ref().relax();
    solve(kEqn);
    bound(k_, kMin_);


    // Re-calculate viscosity and non-linear stress
    correctNonlinearStress(gradU);
}
AJL2005.H

大事なのはここらへん

//- Return the effective diffusivity for k
tmp<volScalarField> DkEff() const
{
    return tmp<volScalarField>
    (
        new volScalarField(
            "DkEff",
            nut_/(sigmak_/(scalar(1)+5.0*fw(5.0))) + nu()
        )
    );
}

//- Return the effective diffusivity for epsilon
tmp<volScalarField> DepsilonEff() const
{
    return tmp<volScalarField>
    (
        new volScalarField(
            "DepsilonEff",
            nut_/(sigmaEps_/(scalar(1)+5.0*fw(5.0))) + nu()
        )
    );
}

ヘッダーファイルの中でもfwは使えた

低Re数型RANSモデル用の壁関数

OpenFOAMに実装されているepsilonの壁関数はepsilonWallFunctionしかないが,なんか無駄な機能がごちゃごちゃしているのでシンプルな低Re数型RANSモデル用の壁関数も実装しておく

>> mkdir -p $WM_PROJECT_USER_DIR/src/myWallFunctions
>> cp -r $WM_PROJECT_DIR/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction $WM_PROJECT_USER_DIR/src/myWallFunctions
>> cd $WM_PROJECT_USER_DIR/src/myWallFunctions
>> mv epsilonWallFunction epsilonLowReWallFunction
>> cd $WM_PROJECT_USER_DIR/src/myWallFunctions/epsilonLowReWallFunction
>> mv epsilonWallFunctionFvPatchScalarField.C epsilonLowReWallFunctionFvPatchScalarField.C
>> mv epsilonWallFunctionFvPatchScalarField.H epsilonLowReWallFunctionFvPatchScalarField.H
>> sed -i s/epsilonWallFunctionFvPatchScalarField/epsilonLowReWallFunctionFvPatchScalarField/g epsilonLowReWallFunctionFvPatchScalarField.*
>> mkdir $WM_PROJECT_USER_DIR/src/myWallFunctions/Make
>> cd $WM_PROJECT_USER_DIR/src/myWallFunctions/Make
>> touch files
>> touch options

epsilonLowReWallFunctionFvPatchScalarField.C

大事な部分だけ

void Foam::epsilonLowReWallFunctionFvPatchScalarField::calculate
(
    const turbulenceModel& turbModel,
    const List<scalar>& cornerWeights,
    const fvPatch& patch,
    scalarField& epsilon0
)
{
    const label patchi = patch.index();

    const scalarField& y = turbModel.y()[patchi];

    const tmp<volScalarField> tk = turbModel.k();
    const volScalarField& k = tk();

    const tmp<scalarField> tnuw = turbModel.nu(patchi);
    const scalarField& nuw = tnuw();

    const tmp<scalarField> tnutw = turbModel.nut(patchi);
    const scalarField& nutw = tnutw();

    const fvPatchVectorField& Uw = turbModel.U().boundaryField()[patchi];

    const scalarField magGradUw(mag(Uw.snGrad()));
    
    forAll(nutw, facei)
    {
        const label celli = patch.faceCells()[facei];
        const scalar w = cornerWeights[facei];
        epsilon0[celli] += w*2.0*k[celli]*nuw[facei]/sqr(y[facei]);
    }
}
epsilonLowReWallFunctionFvPatchScalarField.H

Githubのコードを参照
mtkbirdman.com/OpenFOAM/src/myWallFunctions/epsilonLowReWallFunction/

Make/files
epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.C
LIB = $(FOAM_USER_LIBBIN)/libmyWallFunctions
Make/options

もしかしたら必要ないものが混じっているかも

EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude \
    -I$(LIB_SRC)/fvOptions/lnInclude \
    -I$(LIB_SRC)/transportModels \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
    -I$(LIB_SRC)/finiteVolume/lnInclude

LIB_LIBS = \
    -lincompressibleTurbulenceModels \
    -lturbulenceModels \
    -lfvOptions \
    -lfiniteVolume \
    -lmeshTools

コンパイル・計算の実行(2回目)

コンパイル

AJL2005とepsilonLowReWallFunctionをそれぞれコンパイルする

>> cd $WM_PROJECT_USER_DIR/src/myTurbulenceModels
>> wclean libso
>> wmake libso
wmake libso (myTurbulenceModels)
g++ -std=c++11 -m64 -pthread -DOPENFOAM=2006 -DWM_DP -DWM_LABEL_SIZE=32 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -Wno-attributes -Wno-unknown-pragmas  -O3  -DNoRepository -ftemplate-depth-100 -I/opt/OpenFOAM/OpenFOAM-v2006/src/finiteVolume/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/meshTools/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/fvOptions/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/transportModels -I/opt/OpenFOAM/OpenFOAM-v2006/src/TurbulenceModels/turbulenceModels/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/TurbulenceModels/incompressible/lnInclude  -iquote. -IlnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/OpenFOAM/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/OSspecific/POSIX/lnInclude   -fPIC -c AJL2005/AJL2005.C -o Make/linux64Gcc63DPInt32Opt/AJL2005/AJL2005.o
g++ -std=c++11 -m64 -pthread -DOPENFOAM=2006 -DWM_DP -DWM_LABEL_SIZE=32 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -Wno-attributes -Wno-unknown-pragmas  -O3  -DNoRepository -ftemplate-depth-100 -I/opt/OpenFOAM/OpenFOAM-v2006/src/finiteVolume/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/meshTools/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/fvOptions/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/transportModels -I/opt/OpenFOAM/OpenFOAM-v2006/src/TurbulenceModels/turbulenceModels/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/TurbulenceModels/incompressible/lnInclude  -iquote. -IlnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/OpenFOAM/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/OSspecific/POSIX/lnInclude   -fPIC -shared -Xlinker --add-needed -Xlinker --no-as-needed  Make/linux64Gcc63DPInt32Opt/AJL2005/AJL2005.o -L/opt/OpenFOAM/OpenFOAM-v2006/platforms/linux64Gcc63DPInt32Opt/lib \
    -lincompressibleTurbulenceModels -lturbulenceModels -lfvOptions -lfiniteVolume -lmeshTools  -o /home/mtkbirdman/OpenFOAM/mtkbirdman-v2006/platforms/linux64Gcc63DPInt32Opt/lib/libmyTurbulenceModels.so
>> cd $WM_PROJECT_USER_DIR/src/myWallFunctions
>> wmake libso
wmake libso (myWallFunctions)
g++ -std=c++11 -m64 -pthread -DOPENFOAM=2006 -DWM_DP -DWM_LABEL_SIZE=32 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -Wno-attributes -Wno-unknown-pragmas  -O3  -DNoRepository -ftemplate-depth-100 -I/opt/OpenFOAM/OpenFOAM-v2006/src/finiteVolume/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/meshTools/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/fvOptions/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/transportModels -I/opt/OpenFOAM/OpenFOAM-v2006/src/TurbulenceModels/turbulenceModels/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/TurbulenceModels/incompressible/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/finiteVolume/lnInclude -iquote. -IlnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/OpenFOAM/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/OSspecific/POSIX/lnInclude   -fPIC -c epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.C -o Make/linux64Gcc63DPInt32Opt/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.o
g++ -std=c++11 -m64 -pthread -DOPENFOAM=2006 -DWM_DP -DWM_LABEL_SIZE=32 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -Wno-attributes -Wno-unknown-pragmas  -O3  -DNoRepository -ftemplate-depth-100 -I/opt/OpenFOAM/OpenFOAM-v2006/src/finiteVolume/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/meshTools/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/fvOptions/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/transportModels -I/opt/OpenFOAM/OpenFOAM-v2006/src/TurbulenceModels/turbulenceModels/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/TurbulenceModels/incompressible/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/finiteVolume/lnInclude -iquote. -IlnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/OpenFOAM/lnInclude -I/opt/OpenFOAM/OpenFOAM-v2006/src/OSspecific/POSIX/lnInclude   -fPIC -shared -Xlinker --add-needed -Xlinker --no-as-needed  Make/linux64Gcc63DPInt32Opt/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.o -L/opt/OpenFOAM/OpenFOAM-v2006/platforms/linux64Gcc63DPInt32Opt/lib \
    -lincompressibleTurbulenceModels -lturbulenceModels -lfvOptions -lfiniteVolume -lmeshTools  -o /home/mtkbirdman/OpenFOAM/mtkbirdman-v2006/platforms/linux64Gcc63DPInt32Opt/lib/libmyWallFunctions.so

いくつかのファイルを変更する

0/epsilon
/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2006                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      epsilon;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField   uniform 1;

boundaryField
{
    bottomWall
    {
        type            epsilonLowReWallFunction;
        value           $internalField;
    }
    topWall
    {
        type            epsilonLowReWallFunction;
        value           $internalField;
    }
    sides1_half0
    {
        type            cyclic;
    }
    sides2_half0
    {
        type            cyclic;
    }
    inout1_half0
    {
        type            cyclic;
    }
    inout2_half0
    {
        type            cyclic;
    }
    sides2_half1
    {
        type            cyclic;
    }
    sides1_half1
    {
        type            cyclic;
    }
    inout1_half1
    {
        type            cyclic;
    }
    inout2_half1
    {
        type            cyclic;
    }
}

// ************************************************************************* //

壁関数をepsilonLowReWallFunctionに変更する

system/controlDict
/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2006                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application     simpleFoam;
startFrom       startTime;
startTime       0;
stopAt          endTime;
endTime         5000;
deltaT          1;
writeControl    timeStep;
writeInterval   5000;
purgeWrite      0;
writeFormat     ascii;
writePrecision  6;
writeCompression off;
timeFormat      general;
timePrecision   6;
runTimeModifiable true;

functions
{
    #includeFunc    residuals
    #includeFunc    "sample"
    #includeFunc    "turbulenceFields"
    #includeFunc    "wallShearStress"
}

libs
(
    "libmyTurbulenceModels.so"
    "libmyWallFunctions.so"
);

// ************************************************************************* //

libs("libmyWallFunctions.so");を追加する

system/blockMeshDict
/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2006                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

scale   1;

vertices
(
    (0 0 0)
    (4 0 0)
    (0 1 0)
    (4 1 0)
    (0 2 0)
    (4 2 0)
    (0 0 2)
    (4 0 2)
    (0 1 2)
    (4 1 2)
    (0 2 2)
    (4 2 2)
);

blocks
(
    hex (0 1 3 2 6 7 9 8)   (40 25 30) simpleGrading (1  25 1)
    hex (2 3 5 4 8 9 11 10) (40 25 30) simpleGrading (1 -25 1)
);

edges
(
);

boundary
(
    bottomWall
    {
        type            wall;
        faces           ((0 1 7 6));
    }
    topWall
    {
        type            wall;
        faces           ((4 10 11 5));
    }

    sides1_half0
    {
        type            cyclic;
        neighbourPatch  sides1_half1;
        faces           ((0 2 3 1));
    }
    sides1_half1
    {
        type            cyclic;
        neighbourPatch  sides1_half0;
        faces           ((6 7 9 8));
    }

    sides2_half0
    {
        type            cyclic;
        neighbourPatch  sides2_half1;
        faces           ((2 4 5 3));
    }
    sides2_half1
    {
        type            cyclic;
        neighbourPatch  sides2_half0;
        faces           ((8 9 11 10));
    }

    inout1_half0
    {
        type            cyclic;
        neighbourPatch  inout1_half1;
        faces           ((1 3 9 7));
    }
    inout1_half1
    {
        type            cyclic;
        neighbourPatch  inout1_half0;
        faces           ((0 6 8 2));
    }

    inout2_half0
    {
        type            cyclic;
        neighbourPatch  inout2_half1;
        faces           ((3 5 11 9));
    }
    inout2_half1
    {
        type            cyclic;
        neighbourPatch  inout2_half0;
        faces           ((2 8 10 4));
    }
);

mergePatchPairs
(
);

// ************************************************************************* //

メッシュのGradingを変えてy+の値を小さくする

計算の実行

計算を実行する

結果はこんな感じ

きちんと論文に書いてあるのと同じような分布を得ることができた

まとめ

今回は乱流モデルの実装に挑戦してみた

低Re数型RANSモデルの実装についての資料はあまりなかったのでどこかの誰かの参考になればうれしい

この記事でグラフを書くのに使ったExcelファイルはこちら

Sponsored Link

コメント

タイトルとURLをコピーしました