calculating velocity and epsilon with the variation of height at the inlet bc

 Register Blogs Members List Search Today's Posts Mark Forums Read

December 26, 2019, 08:49
calculating velocity and epsilon with the variation of height at the inlet bc
#1
New Member

Najmiddin
Join Date: Dec 2018
Posts: 16
Rep Power: 4
Hello community,

The next question you are about to read is most probably asked somewhere else before, but I could not find the necessary answer and it might be because of my lack of knowledge in C++

So, I am trying to analyze a 2D rectangular geometry (adapted from a paper). Inlet boundary condition for velocity, k and epsilon is atmospheric boundary layer, as said in the paper. The already existing bc atmBoundaryLayer has different formulae for velocity, k and epsilon. Below is the original formulae used in the above mentioned boundary condition:

Quote:
 U =U = \frac{U^*}{\kappa} ln\left(\frac{z - z_g + z_0}{z_0}\right) U^* = \kappa\frac{U_{ref}}{ln\left(\frac{Z_{ref} + z_0}{z_0}\right)}
and so on, for k and epsilon you can look up, the main point is they are different.

Equations I want to implement are as follows:

Quote:
 U = Uref_*pow(patch().Cf().component(1)/Zref_,0.14) k = 1.5*pow((U_*0.05),2) epsilon = pow(Cmu_,0.75)*pow(k_,1.5)/(kappa_*patch().Cf().component(1))

(patch().Cf().component(1) is the height. I saw the usage of it in some other boundary condition. Basically, my velocity changes with the height. I don't know how to declare what (patch().Cf().component(1) means in the code. I tried to include libraries I thought were necessary for(patch().Cf().component(1) to work. But after each wmake I consistently get the error of "patch is not declared in this scope".

My second question is, as seen in the k and epsilon, k is velocity dependent, which is calculated with the U formula and epsilon is k dependent, and get an error of U and k not being declared.

Below attached are codes and errors:
.C file

Code:
/*---------------------------------------------------------------------------*\
=========                 |
\\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
\\    /   O peration     |
\\  /    A nd           | Copyright (C) 2014-2016 OpenFOAM Foundation
\\/     M anipulation  |
-------------------------------------------------------------------------------
This file is part of OpenFOAM.

OpenFOAM is free software: you can redistribute it and/or modify it
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/>.

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

#include "ABLformulized.H"
#include "volFields.H"

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

namespace Foam
{

// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

atmBoundaryLayer::atmBoundaryLayer()
:
flowDir_(Zero),
//zDir_(Zero),
kappa_(0),
Cmu_(0),
Uref_(0),
//Zref_(0),
z0_(0)
//zGround_(0),
//Ustar_(0)
{}

atmBoundaryLayer::atmBoundaryLayer(const vectorField& p, const dictionary& dict)
:
flowDir_(dict.lookup("flowDir")),
//zDir_(dict.lookup("zDir")),
//kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
//Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
//z0_("z0", dict, p.size()),
//zGround_("zGround", dict, p.size()),
//Ustar_(p.size())
//{
//if (mag(flowDir_) < SMALL || mag(zDir_) < SMALL)
//{
//FatalErrorInFunction
//<< "magnitude of n or z must be greater than zero"
//<< abort(FatalError);
//}

// Ensure direction vectors are normalized
//flowDir_ /= mag(flowDir_);
//zDir_ /= mag(zDir_);

//Ustar_ = kappa_*Uref_/(log((Zref_ + z0_)/z0_));
//}
{}

atmBoundaryLayer::atmBoundaryLayer
(
const atmBoundaryLayer& ptf,
const fvPatchFieldMapper& mapper
)
:
flowDir_(ptf.flowDir_),
//zDir_(ptf.zDir_),
kappa_(ptf.kappa_),
Cmu_(ptf.Cmu_),
Uref_(ptf.Uref_),
z0_(ptf.z0_)
//Zref_(ptf.Zref_),
//z0_(ptf.z0_, mapper),
//zGround_(ptf.zGround_, mapper),
//Ustar_(ptf.Ustar_, mapper)
{}

atmBoundaryLayer::atmBoundaryLayer(const atmBoundaryLayer& blpvf)
:
flowDir_(blpvf.flowDir_),
//zDir_(blpvf.zDir_),
kappa_(blpvf.kappa_),
Cmu_(blpvf.Cmu_),
Uref_(blpvf.Uref_),
//Zref_(blpvf.Zref_),
z0_(blpvf.z0_)
//zGround_(blpvf.zGround_),
//Ustar_(blpvf.Ustar_)
{}

// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

//void atmBoundaryLayer::autoMap(const fvPatchFieldMapper& m)
//{
//z0_.autoMap(m);
//zGround_.autoMap(m);
//Ustar_.autoMap(m);
//}

//void atmBoundaryLayer::rmap
//(
//const atmBoundaryLayer& blptf,
//)
//{
//}

tmp<vectorField> atmBoundaryLayer::U(const vectorField& p) const
{

//vectorField n(patch().nf());

scalarField Un
(
Uref_*pow(patch().Cf().component(1)/z0_,0.14) //(Ustar_/kappa_)*log(((zDir_ & p) - zGround_ + z0_)/z0_)
);

return flowDir_*Un;
}

tmp<scalarField> atmBoundaryLayer::k(const vectorField& p) const
{
return 1.5*pow((U_*0.05),2); //sqr(Ustar_)/sqrt(Cmu_);
}

tmp<scalarField> atmBoundaryLayer::epsilon(const vectorField& p) const

{
//vectorField n(patch().nf());

return pow(Cmu_,0.75)*pow(k_,1.5)/(kappa_*patch().Cf().component(1)); //pow3(Ustar_)/(kappa_*((zDir_ & p) - zGround_ + z0_));
}

void atmBoundaryLayer::write(Ostream& os) const
{
//z0_.writeEntry("z0", os) ;
os.writeKeyword("flowDir")
<< flowDir_ << token::END_STATEMENT << nl;
//os.writeKeyword("zDir")
//<< zDir_ << token::END_STATEMENT << nl;
os.writeKeyword("kappa")
<< kappa_ << token::END_STATEMENT << nl;
os.writeKeyword("Cmu")
<< Cmu_ << token::END_STATEMENT << nl;
os.writeKeyword("Uref")
<< Uref_ << token::END_STATEMENT << nl;
os.writeKeyword("z0")
<< z0_ << token::END_STATEMENT << nl;
//os.writeKeyword("Zref")
//<< Zref_ << token::END_STATEMENT << nl;
//zGround_.writeEntry("zGround", os) ;
}

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

} // End namespace Foam

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


Code:
/*---------------------------------------------------------------------------*\
=========                 |
\\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
\\    /   O peration     |
\\  /    A nd           | Copyright (C) 2014-2016 OpenFOAM Foundation
\\/     M anipulation  |
-------------------------------------------------------------------------------
This file is part of OpenFOAM.

OpenFOAM is free software: you can redistribute it and/or modify it
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/>.

Class
Foam::atmBoundaryLayer

Group
grpRASBoundaryConditions grpInletBoundaryConditions

Description
This class provides functions to evaluate the velocity and turbulence
distributions appropriate for atmospheric boundary layers (ABL).

The profile is derived from the friction velocity, flow direction and
"vertical" direction:

\f[
U = \frac{U^*}{\kappa} ln\left(\frac{z - z_g + z_0}{z_0}\right)
\f]

\f[
k = \frac{(U^*)^2}{\sqrt{C_mu}}
\f]

\f[
\epsilon = \frac{(U^*)^3}{\kappa(z - z_g + z_0)}
\f]

where
\vartable
U^*     | Friction velocity
\kappa  | von Karman's constant
C_mu    | Turbulence viscosity coefficient
z       | Vertical coordinate
z_0     | Surface roughness height [m]
z_g     | Minimum z-coordinate [m]
\endvartable
and
\f[
U^* = \kappa\frac{U_{ref}}{ln\left(\frac{Z_{ref} + z_0}{z_0}\right)}
\f]
where
\vartable
U_{ref} | Reference velocity at \f$Z_{ref}\f$ [m/s]
Z_{ref} | Reference height [m]
\endvartable

Use in the atmBoundaryLayerInletVelocity, atmBoundaryLayerInletK and
atmBoundaryLayerInletEpsilon boundary conditions.

Reference:
D.M. Hargreaves and N.G. Wright,  "On the use of the k-epsilon model
in commercial CFD software to model the neutral atmospheric boundary
layer", Journal of Wind Engineering and Industrial Aerodynamics
95(2007), pp 355-369.

Usage
\table
Property     | Description                      | Required  | Default
flowDir      | Flow direction                   | yes       |
zDir         | Vertical direction               | yes       |
kappa        | von Karman's constant            | no        | 0.41
Cmu          | Turbulence viscosity coefficient | no        | 0.09
Uref         | Reference velocity [m/s]         | yes       |
Zref         | Reference height [m]             | yes       |
z0           | Surface roughness height [m]     | yes       |
zGround      | Minimum z-coordinate [m]         | yes       |
\endtable

Example of the boundary condition specification:
\verbatim
ground
{
type            atmBoundaryLayerInletVelocity;
flowDir         (1 0 0);
zDir            (0 0 1);
Uref            10.0;
Zref            20.0;
z0              uniform 0.1;
zGround         uniform 0.0;
}
\endverbatim

Note
D.M. Hargreaves and N.G. Wright recommend Gamma epsilon in the
k-epsilon model should be changed from 1.3 to 1.11 for consistency.
The roughness height (Er) is given by Er = 20 z0 following the same
reference.

SourceFiles
atmBoundaryLayer.C

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

#ifndef atmBoundaryLayer_H
#define atmBoundaryLayer_H

#include "fvPatchFields.H"

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

namespace Foam
{

/*---------------------------------------------------------------------------*\
Class atmBoundaryLayer Declaration
\*---------------------------------------------------------------------------*/

class atmBoundaryLayer
{
// Private data

//- Flow direction
vector flowDir_;

//- Direction of the z-coordinate
//vector zDir_;

//- Von Karman constant
scalarField kappa_;

//- Turbulent viscosity coefficient
scalarField Cmu_;

//- Reference velocity
scalarField Uref_;

//- Surface roughness height
//const scalar Zref_;

//- Reference height
scalarField z0_;

//- Minimum coordinate value in z direction
//scalarField zGround_;

//- Friction velocity
//scalarField Ustar_;

public:

// Constructors

//- Construct null
atmBoundaryLayer();

//- Construct from the coordinates field and dictionary
atmBoundaryLayer(const vectorField& p, const dictionary&);

//- Construct by mapping given
// atmBoundaryLayer onto a new patch
atmBoundaryLayer
(
const atmBoundaryLayer&,
const fvPatchFieldMapper&
);

//- Construct as copy
atmBoundaryLayer(const atmBoundaryLayer&);

// Member functions

// Access

//- Return flow direction
const vector& flowDir() const
{
return flowDir_;
}

//- Return z-direction
//const vector& zDir() const
//{
//return zDir_;
//}

//- Return friction velocity
//const scalarField& Ustar() const
//{
//return Ustar_;
//}

// Mapping functions

//- Map (and resize as needed) from self given a mapping object
void autoMap(const fvPatchFieldMapper&);

//- Reverse map the given fvPatchField onto this fvPatchField
void rmap(const atmBoundaryLayer&, const labelList&);

// Evaluate functions

//- Return the velocity distribution for the ATM
tmp<vectorField> U(const vectorField& p) const;

//- Return the turbulent kinetic energy distribution for the ATM
tmp<scalarField> k(const vectorField& p) const;

//- Return the turbulent dissipation rate distribution for the ATM
tmp<scalarField> epsilon(const vectorField& p) const;

//- Write
void write(Ostream&) const;
};

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

} // End namespace Foam

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

#endif

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

Errors:
Code:
~/ABLformulized\$ wmake
wmakeLnInclude: linking include files to ./lnInclude
Making dependency list for source file ABLformulizedEpsilon/ABLformulizedEpsilonFvPatchScalarField.C
Making dependency list for source file ABLformulizedK/ABLformulizedKFvPatchScalarField.C
Making dependency list for source file ABLformulizedInletVelocity/ABLformulizedInletVelocityFvPatchVectorField.C
Making dependency list for source file ABLformulized/ABLformulized.C
g++ -std=c++0x -m64 -Dlinux64 -DWM_ARCH_OPTION=64 -DWM_DP -DWM_LABEL_SIZE=32 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -O3  -DNoRepository -ftemplate-depth-100 -I/opt/openfoam4/src/finiteVolume/lnInclude -I/opt/openfoam4/src/meshTools/lnInclude -IlnInclude -I. -I/opt/openfoam4/src/OpenFOAM/lnInclude -I/opt/openfoam4/src/OSspecific/POSIX/lnInclude   -fPIC -c ABLformulized/ABLformulized.C -o Make/linux64GccDPInt32Opt/ABLformulized/ABLformulized.o
ABLformulized/ABLformulized.C: In member function ‘Foam::tmp<Foam::Field<Foam::Vector<double> > > Foam::atmBoundaryLayer::U(const vectorField&) const’:
ABLformulized/ABLformulized.C:142:25: error: ‘patch’ was not declared in this scope
Uref_*pow(patch().Cf().component(1)/z0_,0.14)
^
ABLformulized/ABLformulized.C: In member function ‘Foam::tmp<Foam::Field<double> > Foam::atmBoundaryLayer::k(const vectorField&) const’:
ABLformulized/ABLformulized.C:151:21: error: ‘U_’ was not declared in this scope
return 1.5*pow((U_*0.05),2); //sqr(Ustar_)/sqrt(Cmu_);
^
ABLformulized/ABLformulized.C: In member function ‘Foam::tmp<Foam::Field<double> > Foam::atmBoundaryLayer::epsilon(const vectorField&) const’:
ABLformulized/ABLformulized.C:160:31: error: ‘k_’ was not declared in this scope
return pow(Cmu_,0.75)*pow(k_,1.5)/(kappa_*patch().Cf().component(1)); //pow
^
ABLformulized/ABLformulized.C:160:53: error: ‘patch’ was not declared in this scope
return pow(Cmu_,0.75)*pow(k_,1.5)/(kappa_*patch().Cf().component(1)); //pow
^
/opt/openfoam4/wmake/rules/General/transform:8: recipe for target 'Make/linux64GccDPInt32Opt/ABLformulized/ABLformulized.o' failed
make: *** [Make/linux64GccDPInt32Opt/ABLformulized/ABLformulized.o] Error 1


Last edited by najimaddin96; December 26, 2019 at 08:52. Reason: error in formula

 Tags atmboundarylayerinlet