CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

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

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

Like Tree1Likes
  • 1 Post By najimaddin96

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   December 26, 2019, 09:49
Default 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: 7
najimaddin96 is on a distinguished road
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  |
-------------------------------------------------------------------------------
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/>.

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

#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)),
    kappa_(readScalar(dict.lookup("kappa"))),
    Cmu_(readScalar(dict.lookup("Cmu"))),
    Uref_(readScalar(dict.lookup("Uref"))),
    z0_(readScalar(dict.lookup("z0")))
    //Zref_(readScalar(dict.lookup("Zref"))),
    //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,
    //const labelList& addr
//)
//{
    //z0_.rmap(blptf.z0_, addr);
    //zGround_.rmap(blptf.zGround_, addr);
    //Ustar_.rmap(blptf.Ustar_, addr);
//}


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

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


and .H header file
Code:
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2014-2016 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/>.

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


Thank you for reading and hoping for your quick responses.
Tolga KURUMUS likes this.

Last edited by najimaddin96; December 26, 2019 at 09:52. Reason: error in formula
najimaddin96 is offline   Reply With Quote

Reply

Tags
atmboundarylayerinlet

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Long output in terminal. ssa_cfd OpenFOAM Running, Solving & CFD 1 March 18, 2019 06:25
Difficulty in calculating angular velocity of Savonius turbine simulation alfaruk CFX 14 March 17, 2017 07:08
Inlet & Outlet Velocity BC issue naiter OpenFOAM 3 December 19, 2012 08:14
SimpleFoam k and epsilon bounded nedved OpenFOAM Running, Solving & CFD 1 November 25, 2008 21:21
Velocity inlet K and Epsilon Daniel Bruno FLUENT 6 June 17, 2001 03:51


All times are GMT -4. The time now is 17:25.