zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

corrected y+ calc, re-structured, better handling for walls

Browse Source
This commit is contained in:
andy
2009-08-03 09:53:13 +01:00
parent 35500b5dc3
commit 827f834a93
5 changed files with 154 additions and 114 deletions
Download Patch File Download Diff File Expand all files Collapse all files

62
applications/solvers/incompressible/boundaryFoam/boundaryFoam.C
View File

@ -26,8 +26,8 @@ Application
boundaryFoam
Description
Steady-state solver for 1D turbulent flow, typically to generate boundary
layer conditions at an inlet, for use in a simulation.
Steady-state solver for incompressible, 1D turbulent flow, typically to
generate boundary layer conditions at an inlet, for use in a simulation.
Boundary layer code to calculate the U, k and epsilon distributions.
Used to create inlet boundary conditions for experimental comparisons
@ -42,7 +42,6 @@ Description
#include "wallFvPatch.H"
#include "makeGraph.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
@ -52,6 +51,7 @@ int main(int argc, char *argv[])
#include "createTime.H"
#include "createMesh.H"
#include "createFields.H"
#include "interrogateWallPatches.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -74,66 +74,24 @@ int main(int argc, char *argv[])
UEqn.solve();
// Correct driving force for a constant mass flow rate
// Correct driving force for a constant volume flow rate
dimensionedVector UbarStar = flowMask & U.weightedAverage(mesh.V());
U += (Ubar - UbarStar);
gradP += (Ubar - UbarStar)/(1.0/UEqn.A())().weightedAverage(mesh.V());
label id = y.size() - 1;
scalar wallShearStress =
flowDirection & turbulence->R()()[id] & wallNormal;
scalar yplusWall
// = ::sqrt(mag(wallShearStress))*y[id]/laminarTransport.nu()()[id];
= ::sqrt(mag(wallShearStress))*y[id]/turbulence->nuEff()()[id];
Info<< "Uncorrected Ubar = " << (flowDirection & UbarStar.value())<< tab
<< "pressure gradient = " << (flowDirection & gradP.value()) << tab
<< "min y+ = " << yplusWall << endl;
turbulence->correct();
Info<< "Uncorrected Ubar = " << (flowDirection & UbarStar.value())
<< ", pressure gradient = " << (flowDirection & gradP.value())
<< endl;
#include "evaluateNearWall.H"
if (runTime.outputTime())
{
volSymmTensorField R
(
IOobject
(
"R",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
turbulence->R()
);
runTime.write();
const word& gFormat = runTime.graphFormat();
makeGraph(y, flowDirection & U, "Uf", gFormat);
makeGraph(y, laminarTransport.nu(), gFormat);
makeGraph(y, turbulence->k(), gFormat);
makeGraph(y, turbulence->epsilon(), gFormat);
//makeGraph(y, flowDirection & R & flowDirection, "Rff", gFormat);
//makeGraph(y, wallNormal & R & wallNormal, "Rww", gFormat);
//makeGraph(y, flowDirection & R & wallNormal, "Rfw", gFormat);
//makeGraph(y, sqrt(R.component(tensor::XX)), "u", gFormat);
//makeGraph(y, sqrt(R.component(tensor::YY)), "v", gFormat);
//makeGraph(y, sqrt(R.component(tensor::ZZ)), "w", gFormat);
makeGraph(y, R.component(tensor::XY), "uv", gFormat);
makeGraph(y, mag(fvc::grad(U)), "gammaDot", gFormat);
#include "makeGraphs.H"
}
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

64
applications/solvers/incompressible/boundaryFoam/createFields.H
View File

@ -49,74 +49,14 @@
)
);
dimensionedVector Ubar
(
transportProperties.lookup("Ubar")
);
dimensionedVector Ubar(transportProperties.lookup("Ubar"));
vector flowDirection = (Ubar/mag(Ubar)).value();
tensor flowMask = sqr(flowDirection);
// Search for wall patches faces and store normals
scalar nWallFaces(0);
vector wallNormal(vector::zero);
const fvPatchList& patches = mesh.boundary();
forAll(patches, patchi)
{
const fvPatch& currPatch = patches[patchi];
if (isType<wallFvPatch>(currPatch))
{
forAll(currPatch, facei)
{
nWallFaces++;
if (nWallFaces == 1)
{
wallNormal =
- mesh.Sf().boundaryField()[patchi][facei]
/mesh.magSf().boundaryField()[patchi][facei];
}
else if (nWallFaces == 2)
{
vector wallNormal2 =
mesh.Sf().boundaryField()[patchi][facei]
/mesh.magSf().boundaryField()[patchi][facei];
//- Check that wall faces are parallel
if
(
mag(wallNormal & wallNormal2) > 1.01
||mag(wallNormal & wallNormal2) < 0.99
)
{
Info<< "boundaryFoam: wall faces are not parallel"
<< endl
<< abort(FatalError);
}
}
else
{
Info<< "boundaryFoam: number of wall faces > 2"
<< endl
<< abort(FatalError);
}
}
}
}
//- create position array for graph generation
scalarField y = wallNormal & mesh.C().internalField();
dimensionedVector gradP
(
"gradP",
dimensionSet(0, 1, -2, 0, 0),
vector(0, 0, 0)
vector::zero
);

35
applications/solvers/incompressible/boundaryFoam/evaluateNearWall.H Normal file
View File

@ -0,0 +1,35 @@
{
// Evaluate near-wall behaviour
scalar nu = turbulence->nu().boundaryField()[patchId][faceId];
scalar nut = turbulence->nut()().boundaryField()[patchId][faceId];
symmTensor R = turbulence->devReff()().boundaryField()[patchId][faceId];
scalar epsilon = turbulence->epsilon()()[cellId];
// scalar omega = turbulence->omega()()[cellId];
scalar k = turbulence->k()()[cellId];
scalar Up =
flowDirection & (U[cellId] - U.boundaryField()[patchId][faceId]);
scalar tauw = flowDirection & R & wallNormal;
scalar uTau = ::sqrt(mag(tauw));
scalar yPlus = uTau*y[cellId]/(nu + ROOTVSMALL);
scalar uPlus = Up/(uTau + ROOTVSMALL);
scalar nutPlus = nut/nu;
scalar kPlus = k/(sqr(uTau) + ROOTVSMALL);
scalar epsilonPlus = epsilon*nu/(pow4(uTau) + ROOTVSMALL);
// scalar omegaPlus = omega*nu/(sqr(uTau) + ROOTVSMALL);
scalar Rey = Up*y[cellId]/nu;
Info<< "Rey = " << Rey << ", uTau = " << uTau << ", nut+ = " << nutPlus
<< ", y+ = " << yPlus << ", u+ = " << uPlus
<< ", k+ = " << kPlus << ", epsilon+ = " << epsilonPlus
<< endl;
}

74
applications/solvers/incompressible/boundaryFoam/interrogateWallPatches.H Normal file
View File

@ -0,0 +1,74 @@
// Search for wall patches faces and store normals
label faceId(-1);
label patchId(-1);
label nWallFaces(0);
vector wallNormal(vector::zero);
const fvPatchList& patches = mesh.boundary();
forAll(patches, patchi)
{
const fvPatch& currPatch = patches[patchi];
if (isType<wallFvPatch>(currPatch))
{
const vectorField nf = currPatch.nf();
forAll(nf, facei)
{
nWallFaces++;
if (nWallFaces == 1)
{
wallNormal = -nf[facei];
faceId = facei;
patchId = patchi;
}
else if (nWallFaces == 2)
{
const vector wallNormal2 = -nf[facei];
//- Check that wall faces are parallel
if
(
mag(wallNormal & wallNormal2) > 1.01
|| mag(wallNormal & wallNormal2) < 0.99
)
{
FatalErrorIn(args.executable())
<< "wall faces are not parallel for patches "
<< patches[patchId].name() << " and "
<< currPatch.name() << nl
<< exit(FatalError);
}
}
else
{
FatalErrorIn(args.executable()) << "number of wall faces > 2"
<< nl << exit(FatalError);
}
}
}
}
if (nWallFaces == 0)
{
FatalErrorIn(args.executable()) << "No wall patches identified"
<< exit(FatalError);
}
else
{
Info<< "Generating wall data for patch: " << patches[patchId].name() << endl;
}
// store local id of near-walll cell to process
label cellId = patches[patchId].faceCells()[faceId];
// create position array for graph generation
scalarField y =
wallNormal
& (mesh.C().internalField() - mesh.C().boundaryField()[patchId][faceId]);
Info<< " Height to first cell centre y0 = " << y[cellId] << endl;

33
applications/solvers/incompressible/boundaryFoam/makeGraphs.H Normal file
View File

@ -0,0 +1,33 @@
volSymmTensorField R
(
IOobject
(
"R",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
turbulence->R()
);
runTime.write();
const word& gFormat = runTime.graphFormat();
makeGraph(y, flowDirection & U, "Uf", gFormat);
makeGraph(y, turbulence->nu(), gFormat);
makeGraph(y, turbulence->k(), gFormat);
makeGraph(y, turbulence->epsilon(), gFormat);
makeGraph(y, flowDirection & R & flowDirection, "Rff", gFormat);
makeGraph(y, wallNormal & R & wallNormal, "Rww", gFormat);
makeGraph(y, flowDirection & R & wallNormal, "Rfw", gFormat);
makeGraph(y, sqrt(mag(R.component(tensor::XX))), "u", gFormat);
makeGraph(y, sqrt(mag(R.component(tensor::YY))), "v", gFormat);
makeGraph(y, sqrt(mag(R.component(tensor::ZZ))), "w", gFormat);
makeGraph(y, R.component(tensor::XY), "uv", gFormat);
makeGraph(y, mag(fvc::grad(U)), "gammaDot", gFormat);