Commit e5a710d3 authored by mattijs's avatar mattijs
Browse files

Merge branch 'master' of /home/noisy3/OpenFOAM/OpenFOAM-dev

parents eba6c985 2c62c83b
#!/bin/sh
cd ${0%/*} || exit 1 # run from this directory
if [ "$PWD" != "$WM_PROJECT_DIR" ]
then
echo "Error: Current directory is not \$WM_PROJECT_DIR"
echo " The environment variables are inconsistent with the installation."
echo " Check the OpenFOAM entries in your dot-files and source them."
exit 1
fi
# wmake is required for subsequent targets
( cd wmake/src && make )
......
#!/bin/sh
cd ${0%/*} || exit 1 # run from this directory
if [ "$PWD" != "$WM_PROJECT_DIR/applications" ]
then
echo "Error: Current directory is not \$WM_PROJECT_DIR/applications"
echo " The environment variables are inconsistent with the installation."
echo " Check the OpenFOAM entries in your dot-files and source them."
exit 1
fi
set -x
wmake all solvers
......
......@@ -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"
......
......@@ -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
);
{
// 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
// 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;
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);
......@@ -82,7 +82,7 @@
(
"kappa",
wallFunctionDict,
0.4187
0.41
)
);
......@@ -92,7 +92,7 @@
(
"E",
wallFunctionDict,
9.0
9.8
)
);
......
......@@ -304,7 +304,7 @@
(
"kappa",
wallFunctionDict,
0.4187
0.41
)
);
......@@ -314,7 +314,7 @@
(
"E",
wallFunctionDict,
9.0
9.8
)
);
......
#include "Xstream.H"
#include "GLstream.H"
#include "PSstream.H"
#include "shapes2D.H"
#include "IStringStream.H"
using namespace Foam;
int main()
{
colour mauve("mauve", 1, 0, 1);
lineStyle solid("Solid", 2.0, 2.0, IStringStream("1(1.0)")());
lineStyle broken("Broken", 2.0, 10.0, IStringStream("4(1 1 4 1)")());
//Xstream wind
GLstream wind
//PSstream wind
(
"GsTest",
primary("Black"),
primary("White"),
0.5, 0.5, 0.5, 0.5, 500, 500
);
do
{
wind << rectangle2D(point2D(0.0, 0.0), point2D(100.0, 100.0));
wind.setColour(mauve);
wind.setLineStyle(solid);
wind << line2D(point2D(0.0, 0.0), point2D(0.0, 200.0));
wind << line2D(point2D(0.0, 200.0), point2D(200.0, 200.0));
//wind.setLineStyle(broken);
wind << line2D(point2D(200.0, 200.0), point2D(200.0, 0.0));
wind << line2D(point2D(200.0, 0.0), point2D(0.0, 0.0));
wind << string2D(point2D(200.0, 0.0), "Hi there");
} while (wind.waitForEvent());
return 0;
}
GsTest.C
EXE = $(FOAM_USER_APPBIN)/GsTest
EXE_INC = -I$(LIB_SRC)/Gstream/lnInclude
EXE_LIBS = -lGstream -lGL $(XLIBS)
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ 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 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
Description
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "IStringStream.H"
#include "myBoundBox.H"
#include "myBoundBoxList.H"
#include "octree.H"
#include "octreeData.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Main program:
int main(int argc, char *argv[])
{
argList::validOptions.insert("x1", "X1");
argList::validOptions.insert("y1", "Y1");
argList::validOptions.insert("z1", "Z1");
argList::validOptions.insert("x2", "X2");
argList::validOptions.insert("y2", "Y2");
argList::validOptions.insert("z2", "Z2");
# include "setRootCase.H"
# include "createTime.H"
# include "createMesh.H"
// Calculate BB of all cells
myBoundBoxList allBb(mesh.nCells());
const pointField& allPoints = mesh.points();
vectorField bbMin(mesh.nCells());
bbMin = vector(GREAT, GREAT, GREAT);
vectorField bbMax(mesh.nCells());
bbMax = vector(-GREAT, -GREAT, -GREAT);
const labelListList& pCells = mesh.pointCells();
forAll(pCells, pointi)
{
const point& vertCoord = allPoints[pointi];
const labelList& cells = pCells[pointi];
forAll(cells, celli)
{
label cellNum = cells[celli];
bbMin[cellNum].x() = min(bbMin[cellNum].x(), vertCoord.x());
bbMin[cellNum].y() = min(bbMin[cellNum].y(), vertCoord.y());
bbMin[cellNum].z() = min(bbMin[cellNum].z(), vertCoord.z());
bbMax[cellNum].x() = max(bbMax[cellNum].x(), vertCoord.x());
bbMax[cellNum].y() = max(bbMax[cellNum].y(), vertCoord.y());
bbMax[cellNum].z() = max(bbMax[cellNum].z(), vertCoord.z());
}
}
forAll(allBb, celli)
{
allBb[celli] = myBoundBox(bbMin[celli], bbMax[celli]);
}
myBoundBox meshBb(allPoints);
scalar typDim = meshBb.minDim()/111;
myBoundBox shiftedBb
(
meshBb.min(),
point
(
meshBb.max().x() + typDim,
meshBb.max().y() + typDim,
meshBb.max().z() + typDim
)
);
Info<< "Mesh" << endl;
Info<< " bounding box :" << shiftedBb << endl;
Info<< " typical dimension:" << shiftedBb.typDim() << endl;
/*
* Now we have allBb and shiftedBb
*/
// Construct table of subset of cells
labelList cellIndices(10);
cellIndices[0] = 1433;
cellIndices[1] = 1434;
cellIndices[2] = 1435;
cellIndices[3] = 1436;
cellIndices[4] = 1437;
cellIndices[5] = 1438;
cellIndices[6] = 1439;
cellIndices[7] = 1440;
cellIndices[8] = 1441;
cellIndices[9] = 1442;
// Get the corresponding bounding boxes
forAll(cellIndices, i)
{
allBb[i] = allBb[cellIndices[i]];
}
allBb.setSize(cellIndices.size());
// Wrap indices and mesh information into helper object