diff --git a/applications/solvers/incompressible/boundaryFoam/boundaryFoam.C b/applications/solvers/incompressible/boundaryFoam/boundaryFoam.C index 2ee69c120e..62a8db992b 100644 --- a/applications/solvers/incompressible/boundaryFoam/boundaryFoam.C +++ b/applications/solvers/incompressible/boundaryFoam/boundaryFoam.C @@ -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" diff --git a/applications/solvers/incompressible/boundaryFoam/createFields.H b/applications/solvers/incompressible/boundaryFoam/createFields.H index b7b8b32229..e99a990013 100644 --- a/applications/solvers/incompressible/boundaryFoam/createFields.H +++ b/applications/solvers/incompressible/boundaryFoam/createFields.H @@ -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(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 ); diff --git a/applications/solvers/incompressible/boundaryFoam/evaluateNearWall.H b/applications/solvers/incompressible/boundaryFoam/evaluateNearWall.H new file mode 100644 index 0000000000..b9683b9e68 --- /dev/null +++ b/applications/solvers/incompressible/boundaryFoam/evaluateNearWall.H @@ -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; +} \ No newline at end of file diff --git a/applications/solvers/incompressible/boundaryFoam/interrogateWallPatches.H b/applications/solvers/incompressible/boundaryFoam/interrogateWallPatches.H new file mode 100644 index 0000000000..27ddd24544 --- /dev/null +++ b/applications/solvers/incompressible/boundaryFoam/interrogateWallPatches.H @@ -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(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; + diff --git a/applications/solvers/incompressible/boundaryFoam/makeGraphs.H b/applications/solvers/incompressible/boundaryFoam/makeGraphs.H new file mode 100644 index 0000000000..7367dc7c7d --- /dev/null +++ b/applications/solvers/incompressible/boundaryFoam/makeGraphs.H @@ -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);