看到 cfd online 的后处理代码：compressibleyPlusLES，感觉修改下作为后处理工具，不失为一个更简单直接的方法。
volScalarField yPlusTemp
(
IOobject
(
"yPlusTemp",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("yPlusTemp", dimless, 0.0)
);
volVectorField uPlus
(
IOobject
(
"uPlus",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedVector("uPlus", dimless, vector::zero)
);
volScalarField uTau
(
IOobject
(
"uTau",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("uTau", dimVelocity, 0.0)
);
volScalarField yPlus
(
IOobject
(
"yPlus",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("yPlus", dimless, 0.0)
);
Info<< "Reading field UMean\n" << endl;
volVectorField UMean
(
IOobject
(
"UMean",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
),
mesh
);
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
),
mesh
);
Info<< "Reading field rho\n" << endl;
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
),
mesh
);
# include "compressibleCreatePhi.H"
Info<< "successs1\n" << endl;
autoPtr<basicThermo> pThermo
(
basicThermo::New(mesh)
);
basicThermo& thermo = pThermo();
Info<< "successs2\n" << endl;
autoPtr<compressible::LESModel> sgsModel
(
compressible::LESModel::New(rho, UMean, phi, thermo)
);
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
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 3 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, see <http://www.gnu.org/licenses/>.
Application
yPlusLES
Description
Calculates and reports yPlus and uPlus for all wall patches, for the specified times, for compressible LES.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
//#include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
//aali
#include "compressible/LES/LESModel/LESModel.H"
#include "basicThermo.H"
//ali
#include "wallFvPatch.H"
#include "LESModel.H"
#include "nearWallDist.H"
#include "wallDist.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);
Info<< "Time = " << runTime.timeName() << endl;
fvMesh::readUpdateState state = mesh.readUpdate();
// Wall distance
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"
volScalarField::GeometricBoundaryField d = nearWallDist(mesh).y();
volScalarField muEff(sgsModel->muEff());
const fvPatchList& patches = mesh.boundary();
dimensionedScalar uTauAvg("uTauAvg", dimVelocity, 0);
const volScalarField muLam(sgsModel->mu());
scalar nPatch = 0;
Info<< "Summary: " << nl << endl;
forAll(patches, patchi)
{
const fvPatch& currPatch = patches[patchi];
if (typeid(currPatch) == typeid(wallFvPatch))//isA<wallFvPatch>(currPatch))
{
yPlusTemp.boundaryField()[patchi] =
d[patchi]
*sqrt
(
muEff.boundaryField()[patchi]/rho.boundaryField()[patchi]
*mag(UMean.boundaryField()[patchi].snGrad())
)
/muLam.boundaryField()[patchi]/rho.boundaryField()[patchi];
const scalarField& YpTemp = yPlusTemp.boundaryField()[patchi];
//effective viscosity used in the calculation of u*
uTau.boundaryField()[patchi] =
sqrt
(
sgsModel->muEff()
*mag(UMean.boundaryField()[patchi].snGrad())
/sgsModel->rho()
);
const fvPatchScalarField& uTauWall = uTau.boundaryField()[patchi];
dimensionedScalar uTauTmp("uTauTmp", dimVelocity, average(uTauWall));
uTauAvg += uTauTmp;
nPatch ++;
Info<< "Patch " << patchi
<< " named " << currPatch.name()
<< " y+ : 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 is = u* x y/nu
yPlus = y.y() * uTauAvg / ((sgsModel->mu()) / (sgsModel->rho()));
uPlus = UMean / uTauAvg;
Info << "Writing yPlus and uPlus to corresponding fields." << nl <<endl;
yPlus.write();
uPlus.write();
// Info<< "Writing yPlus to field "
// << yPlus.name() << nl << endl;
//
// yPlus.write();
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //
