UplusとYplusを算出するユーティリティ

[mei]

皆様，こんばんは

CFDonlineよりUplusYplusを算出するユーティリティを見つけ使ってみたのですが，Uplusが対数速度分布にのらず困っています．
別途壁面せん断応力を計算しており，それと比較してみたのですが，摩擦速度が計算できてないように思えました．

プログラムにおいて確かに
                uTau.boundaryField()[patchi] =
                    sqrt
                    (
            LESModel->nu()
                   *mag(UMean.boundaryField()[patchi].snGrad())
                    );
としてあるのですが，逆算すると摩擦速度が0.4ほど．別途計算した壁面せん断応力から算出したものが0.6ほどでして・・・
どうすればいいか皆目検討もつきません．
どなたかご教授お願いいたします．

以下プログラムです．あと，もとはRANSであるのをLES用に変更しています．多分・・・
License
    This file is part of OpenFOAM.

    OpenFOAM is free software; you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by the
    Free Software Foundation; either version 2 of the License, or (at your
    option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM; if not, write to the Free Software Foundation,
    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA

Application
    plusPostLES

Description

    calculates y+ and u+ fields for wall-bounded flows computed with
    one of the available low-Re LES (no wall function!) turbulence
    models. More specifically it

    ::     calculates and outputs y+ (avg., min., max.) based on the
    velocity gradient at the wall 

    ::    calculates and outputs the wall averaged friction velocity

    ::  writes fields of y+ and U+ to the corresponding time directory

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"
#include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
#include "incompressible/LES/LESModel/LESModel.H"
#include "RASModel.H"
#include "LESModel.H"
#include "nearWallDist.H"
#include "wallDist.H"
#include "wallFvPatch.H"

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

int main(int argc, char *argv[])
{
    timeSelector::addOptions();

    #include "setRootCase.H"
    #include "createTime.H"
    instantList timeDirs = timeSelector::select0(runTime, args);
    #include "createMesh.H"

    forAll(timeDirs, timeI)
    {
        runTime.setTime(timeDirs[timeI], timeI);

        fvMesh::readUpdateState state = mesh.readUpdate();

        wallDist y(mesh, true);

        if (timeI == 0 || state != fvMesh::UNCHANGED)
        {
            Info<< "Calculating wall distance\n" <<endl;
            Info<< "Writing wall distance to field " << y.name() << nl << endl;
            y.write();
        }

    #include "createFields.H"

        const fvPatchList& patches = mesh.boundary();
       
        dimensionedScalar uTauAvg("uTauAvg", dimVelocity, 0);

        scalar nPatch = 0;
               
    Info<< "Summary: " << nl << endl;

        forAll(patches, patchi)
        {
            const fvPatch& currPatch = patches[patchi];
            const tmp<volScalarField> nu = LESModel->nu();

            if (typeid(currPatch) == typeid(wallFvPatch))
            {
                yPlusTemp.boundaryField()[patchi] =
                    d[patchi]
                   *sqrt
                    (
                       nu
                        *mag(U.boundaryField()[patchi].snGrad())
                    )
                   /nu;
               
        const scalarField& YpTemp = yPlusTemp.boundaryField()[patchi];

                uTau.boundaryField()[patchi] =
                    sqrt
                    (
            LESModel->nu()
                   *mag(UMean.boundaryField()[patchi].snGrad())
                    );

                const fvPatchScalarField& uTauWall = uTau.boundaryField()[patchi];
               
                dimensionedScalar uTauTmp("uTauTmp", dimVelocity, average(uTauWall));
               
                uTauAvg += uTauTmp;

                nPatch ++;

                Info<< "  y+ for Patch " << patchi
                    << " named " << currPatch.name() << ":"
                    << " min: " << min(YpTemp) << " max: " << max(YpTemp)
                    << " average: " << average(YpTemp)
                    << " avgUGradWall: " <<  average(mag(U.boundaryField()[patchi].snGrad())) << nl << endl;
            }
        }

        uTauAvg /= nPatch;
       
        Info << "  avg. friction velocity uTau is: "
             << uTauAvg.value() << " (averaged over " << nPatch << " wall(s))" << nl <<endl;

        yPlus = y.y() * uTauAvg / LESModel->nu();
        uPlus = UMean / uTauAvg;
       
        Info << "Writing yPlus and uPlus to corresponding fields." << nl <<endl;
        yPlus.write();
        uPlus.write();

    }

    Info<< "End\n" << endl;

    return 0;
}

[ONO Hiroki]

plusPostLESでは，乱流モデルに依存する変数(nuSgs,k(sgs)など)を使用せずに，直接，壁面と壁面第1セルの速度勾配をとっていることから，
「壁面第1セルが粘性底層内に位置し，かつnuSgsが壁面で(ほぼ)0へ減衰していること」を前提として作られていると推定できます。
(Descriptionの"for wall-bounded flows computed with one of the available low-Re LES (no wall function!) turbulence models."からもその趣旨が読み取れます。 )


mei様が別途計算に使用したものは，おそらくOF標準のwallShearStressをLES対応したものと思われますが，
こちらはReffをもとに計算するので，もしnuSgsがnuと比べて十分小さくなければ，その影響が入るのでplusPostLESとは違った値となるでしょう。

壁面近傍でnuSgsが十分に小さければ，上記に起因する差異はほとんどないはずですが，もう一点，
plusPostLESではuTauを計算する際に速度の計算結果として，UではなくUMeanを使用しています。
createFields.Hを見てみないとこのUMeanの正体はわかりませんが，おそらくfieldAverageで生成された時間平均値でしょうから，
もう一方の計算でも同様にUMeanを使わないと，結果は完全には一致しないでしょう。




2014年1月15日水曜日 3時44分00秒 UTC+9 mei:

[mei]

ONO Hiroki様

ご返信ありがとうございます．
壁面第1セルは粘性底層内に位置していますが，nuSgsが壁面で(ほぼ)0へ減衰しているか確認してみます．
また，ご指摘の通りにwallShearStressをLESに修正したものを使っています．確認したらUを使っていましたので，こちらを変更して再度計算してみます．

ありがとうございます．道が出来ました．