PR

【OpenFOAM v2006】LES deltaの実装

2020.11.09

OpenFOAMにおけるLES deltaの実装について説明する

ソースコードはGithubを参照
≫mtkbirdman.com/OpenFOAM/

キーワード：LES delta, implementation

スポンサーリンク

はじめに

OpenFOAMで，LESのdeltaを実装する

実装するdeltaは以下の論文に出てくるもの

≫Abe, K. I. (2013). An improved anisotropy-resolving subgrid-scale model with the aid of a scale-similarity modeling concept. International Journal of Heat and Fluid Flow, 39, 42-52. https://doi.org/10.1016/j.ijheatfluidflow.2012.12.001

C言語については以下のサイトがおすすめ
≫苦しんで覚えるＣ言語 – 苦しんで覚えるC言語

C++は次のサイトを読んで完全に理解した
≫C++入門

このLES deltaを使った乱流モデルの実装はこちら
≫【OpenFOAM v2006】非線形SGSモデルの実装

概要

上の論文で使用されている\(\Delta\)は次のように定式化されている

\begin{eqnarray}
\Delta=\sqrt{the~maximum~area~among~the~faces~of~a~cell} \tag{14}
\end{eqnarray}

これをmaxSqrtFaceDeltaと名付けて実装していく

実装

ベースモデルのコピー

次のコマンドをたたいてmaxDeltaxyzを$WM_PROJECT_USER_DIR/src/myTurbulenceModelsにコピーする

>> cp -r $WM_PROJECT_DIR/src/TurbulenceModels/turbulenceModels/LES/LESdeltas/maxDeltaxyz $WM_PROJECT_USER_DIR/src/myTurbulenceModels

次のコマンドで，乱流モデルの名前をmaxDeltaxyzからmaxSqrtFaceDeltaに変更する

>> cd $WM_PROJECT_USER_DIR/src/myTurbulenceModels/
>> mv maxDeltaxyzmaxSqrtFaceDelta
>> cd maxSqrtFaceDelta
>> mv maxDeltaxyz.H maxSqrtFaceDelta.H
>> mv maxDeltaxyz.C maxSqrtFaceDelta.C
>> sed -i -e "s/maxDeltaxyz/maxSqrtFaceDelta/g" maxSqrtFaceDelta.*

プログラムの変更

maxSqrtFaceDelta.H

maxSqrtFaceDelta.HのコードはGithubのコードを参照
≫mtkbirdman.com/OpenFOAM/src/myTurbulenceModels/maxSqrtFaceDelta/maxSqrtFaceDelta.H

注意点は以下の通り

  • TypeName("maxDeltaxyz");は上のコマンドでは変更されないので手動でTypeName("maxSqrtFace");に変更する
maxSqrtFaceDelta.C

maxSqrtFaceDelta.CのコードはGithubのコードを参照
≫mtkbirdman.com/OpenFOAM/src/myTurbulenceModels/maxSqrtFaceDelta/maxSqrtFaceDelta.C

重要な点のみ説明する

void Foam::LESModels::maxSqrtFaceDelta::calcDelta()
{
    const fvMesh& mesh = turbulenceModel_.mesh();

    label nD = mesh.nGeometricD();

    const cellList& cells = mesh.cells();
    scalarField sqrtSmax(cells.size());

    forAll(cells, celli)
    {
        scalar deltaMaxTmp = 0.0;
        const labelList& cFaces = cells[celli];

        forAll(cFaces, cFacei)
        {
            label facei = cFaces[cFacei];

            vector S = mesh.faceAreas()[facei];
            scalar tmp = sqrt(mag(S));
            if (tmp > deltaMaxTmp)
            {
                deltaMaxTmp = tmp;
            }
        }

        sqrtSmax[celli] = deltaCoeff_*deltaMaxTmp;
    }

    if (nD == 3)
    {
        delta_.primitiveFieldRef() = sqrtSmax;
    }
    else if (nD == 2)
    {
        WarningInFunction
            << "Case is 2D, LES is not strictly applicable" << nl
            << endl;

        delta_.primitiveFieldRef() = sqrtSmax;
    }
    else
    {
        FatalErrorInFunction
            << "Case is not 3D or 2D, LES is not applicable"
            << exit(FatalError);
    }

    // Handle coupled boundaries
    delta_.correctBoundaryConditions();
}

基本的にはmaxDeltaxyzのループ構造を利用している

セルのラベルのリストを返すmesh.cells()や面の面積に対応する外積ベクトルを返すmesh.faceAreas()はprimitiveMeshというクラスで実装されている

コンパイル・計算・結果

親記事(≫【OpenFOAM v2006】非線形SGSモデルの実装)を参考

まとめ

最初mesh.faceAreas()がセルの面の面積のスカラーを集めたベクトルを返していると勘違いしていたためてこずったが，元のソースコード(≫primitiveMeshFaceCentresAndAreas.C)を確認したことでミスに気づくことができた

公式のドキュメントやコードを読む癖がついてきて何よりである

スポンサーリンク

コメント