【OpenFOAM v2006】LES deltaの実装

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 Flow39, 42-52. https://doi.org/10.1016/j.ijheatfluidflow.2012.12.001

C言語については以下のサイトがおすすめ
苦しんで覚える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)を確認したことでミスに気づくことができた

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

Sponsored Link

コメント

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