mirror of
https://github.com/OpenFOAM/OpenFOAM-6.git
synced 2025-12-08 06:57:46 +00:00
atmosphericModels: Added Lopes da Costa porosity and turbulence models
Specialized variants of the power law porosity and k epsilon turbulence models
developed to simulate atmospheric flow over forested and non-forested complex
terrain.
Class
Foam::powerLawLopesdaCosta
Description
Variant of the power law porosity model with spatially varying
drag coefficient
given by:
\f[
S = -\rho C_d \Sigma |U|^{(C_1 - 1)} U
\f]
where
\vartable
\Sigma | Porosity surface area per unit volume
C_d | Model linear coefficient
C_1 | Model exponent coefficient
\endvartable
Reference:
\verbatim
Costa, J. C. P. L. D. (2007).
Atmospheric flow over forested and non-forested complex terrain.
\endverbatim
Class
Foam::RASModels::kEpsilonLopesdaCosta
Description
Variant of the standard k-epsilon turbulence model with additional source
terms to handle the changes in turbulence in porous regions represented by
the powerLawLopesdaCosta porosity model.
Reference:
\verbatim
Costa, J. C. P. L. D. (2007).
Atmospheric flow over forested and non-forested complex terrain.
\endverbatim
The default model coefficients are
\verbatim
kEpsilonLopesdaCostaCoeffs
{
Cmu 0.09;
C1 1.44;
C2 1.92;
sigmak 1.0;
sigmaEps 1.3;
}
\endverbatim
Tutorial case to follow.
This commit is contained in:
Henry Weller
@ -81,6 +81,7 @@ wmake $targetType rigidBodyDynamics
|
||||
wmake $targetType rigidBodyMeshMotion
|
||||
wmake $targetType waves
|
||||
wmake $targetType semiPermeableBaffle
|
||||
wmake $targetType atmosphericModels
|
||||
|
||||
|
||||
#------------------------------------------------------------------------------
|
||||
|
||||
@ -4,4 +4,7 @@ derivedFvPatchFields/atmBoundaryLayerInletK/atmBoundaryLayerInletKFvPatchScalarF
|
||||
derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.C
|
||||
derivedFvPatchFields/nutkAtmRoughWallFunction/nutkAtmRoughWallFunctionFvPatchScalarField.C
|
||||
|
||||
atmosphericTurbulentTransportModels.C
|
||||
porosityModels/powerLawLopesdaCosta/powerLawLopesdaCosta.C
|
||||
|
||||
LIB = $(FOAM_LIBBIN)/libatmosphericModels
|
||||
|
||||
35
src/atmosphericModels/atmosphericTurbulentTransportModels.C
Normal file
35
src/atmosphericModels/atmosphericTurbulentTransportModels.C
Normal file
@ -0,0 +1,35 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2018 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 "turbulentTransportModels.H"
|
||||
|
||||
// -------------------------------------------------------------------------- //
|
||||
// RAS models
|
||||
// -------------------------------------------------------------------------- //
|
||||
|
||||
#include "kEpsilonLopesdaCosta.H"
|
||||
makeRASModel(kEpsilonLopesdaCosta);
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -2,7 +2,7 @@
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2014-2016 OpenFOAM Foundation
|
||||
\\ / A nd | Copyright (C) 2014-2018 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
|
||||
@ -0,0 +1,477 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2018 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 "kEpsilonLopesdaCosta.H"
|
||||
#include "fvOptions.H"
|
||||
#include "explicitPorositySource.H"
|
||||
#include "bound.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
namespace RASModels
|
||||
{
|
||||
|
||||
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
|
||||
|
||||
template<class BasicTurbulenceModel>
|
||||
void kEpsilonLopesdaCosta<BasicTurbulenceModel>::setPorosityCoefficient
|
||||
(
|
||||
volScalarField::Internal& C,
|
||||
const porosityModels::powerLawLopesdaCosta& pm
|
||||
)
|
||||
{
|
||||
if (pm.dict().found(C.name()))
|
||||
{
|
||||
const labelList& cellZoneIDs = pm.cellZoneIDs();
|
||||
|
||||
const scalar Cpm = readScalar(pm.dict().lookup(C.name()));
|
||||
|
||||
forAll(cellZoneIDs, zonei)
|
||||
{
|
||||
const labelList& cells =
|
||||
this->mesh_.cellZones()[cellZoneIDs[zonei]];
|
||||
|
||||
forAll(cells, i)
|
||||
{
|
||||
const label celli = cells[i];
|
||||
C[celli] = Cpm;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template<class BasicTurbulenceModel>
|
||||
void kEpsilonLopesdaCosta<BasicTurbulenceModel>::setCdSigma
|
||||
(
|
||||
volScalarField::Internal& C,
|
||||
const porosityModels::powerLawLopesdaCosta& pm
|
||||
)
|
||||
{
|
||||
if (pm.dict().found(C.name()))
|
||||
{
|
||||
const labelList& cellZoneIDs = pm.cellZoneIDs();
|
||||
const scalarField& Sigma = pm.Sigma();
|
||||
|
||||
const scalar Cpm = readScalar(pm.dict().lookup(C.name()));
|
||||
|
||||
forAll(cellZoneIDs, zonei)
|
||||
{
|
||||
const labelList& cells =
|
||||
this->mesh_.cellZones()[cellZoneIDs[zonei]];
|
||||
|
||||
forAll(cells, i)
|
||||
{
|
||||
const label celli = cells[i];
|
||||
C[celli] = Cpm*Sigma[celli];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template<class BasicTurbulenceModel>
|
||||
void kEpsilonLopesdaCosta<BasicTurbulenceModel>::setPorosityCoefficients()
|
||||
{
|
||||
fv::options::optionList& fvOptions(fv::options::New(this->mesh_));
|
||||
|
||||
forAll(fvOptions, i)
|
||||
{
|
||||
if (isA<fv::explicitPorositySource>(fvOptions[i]))
|
||||
{
|
||||
const fv::explicitPorositySource& eps =
|
||||
refCast<const fv::explicitPorositySource>(fvOptions[i]);
|
||||
|
||||
if (isA<porosityModels::powerLawLopesdaCosta>(eps.model()))
|
||||
{
|
||||
const porosityModels::powerLawLopesdaCosta& pm =
|
||||
refCast<const porosityModels::powerLawLopesdaCosta>
|
||||
(
|
||||
eps.model()
|
||||
);
|
||||
|
||||
setPorosityCoefficient(Cmu_, pm);
|
||||
setPorosityCoefficient(C1_, pm);
|
||||
setPorosityCoefficient(C2_, pm);
|
||||
setPorosityCoefficient(sigmak_, pm);
|
||||
setPorosityCoefficient(sigmaEps_, pm);
|
||||
|
||||
setCdSigma(CdSigma_, pm);
|
||||
setPorosityCoefficient(betap_, pm);
|
||||
setPorosityCoefficient(betad_, pm);
|
||||
setPorosityCoefficient(C4_, pm);
|
||||
setPorosityCoefficient(C5_, pm);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template<class BasicTurbulenceModel>
|
||||
void kEpsilonLopesdaCosta<BasicTurbulenceModel>::correctNut()
|
||||
{
|
||||
this->nut_ = Cmu_*sqr(k_)/epsilon_;
|
||||
this->nut_.correctBoundaryConditions();
|
||||
fv::options::New(this->mesh_).correct(this->nut_);
|
||||
|
||||
BasicTurbulenceModel::correctNut();
|
||||
}
|
||||
|
||||
|
||||
template<class BasicTurbulenceModel>
|
||||
tmp<fvScalarMatrix> kEpsilonLopesdaCosta<BasicTurbulenceModel>::kSource
|
||||
(
|
||||
const volScalarField::Internal& magU,
|
||||
const volScalarField::Internal& magU3
|
||||
) const
|
||||
{
|
||||
return fvm::Su(CdSigma_*(betap_*magU3 - betad_*magU*k_()), k_);
|
||||
}
|
||||
|
||||
|
||||
template<class BasicTurbulenceModel>
|
||||
tmp<fvScalarMatrix>
|
||||
kEpsilonLopesdaCosta<BasicTurbulenceModel>::epsilonSource
|
||||
(
|
||||
const volScalarField::Internal& magU,
|
||||
const volScalarField::Internal& magU3
|
||||
) const
|
||||
{
|
||||
return fvm::Su
|
||||
(
|
||||
CdSigma_
|
||||
*(C4_*betap_*epsilon_()/k_()*magU3 - C5_*betad_*magU*epsilon_()),
|
||||
epsilon_
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
|
||||
|
||||
template<class BasicTurbulenceModel>
|
||||
kEpsilonLopesdaCosta<BasicTurbulenceModel>::kEpsilonLopesdaCosta
|
||||
(
|
||||
const alphaField& alpha,
|
||||
const rhoField& rho,
|
||||
const volVectorField& U,
|
||||
const surfaceScalarField& alphaRhoPhi,
|
||||
const surfaceScalarField& phi,
|
||||
const transportModel& transport,
|
||||
const word& propertiesName,
|
||||
const word& type
|
||||
)
|
||||
:
|
||||
eddyViscosity<RASModel<BasicTurbulenceModel>>
|
||||
(
|
||||
type,
|
||||
alpha,
|
||||
rho,
|
||||
U,
|
||||
alphaRhoPhi,
|
||||
phi,
|
||||
transport,
|
||||
propertiesName
|
||||
),
|
||||
|
||||
Cmu_
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"Cmu",
|
||||
this->runTime_.timeName(),
|
||||
this->mesh_
|
||||
),
|
||||
this->mesh_,
|
||||
dimensioned<scalar>::lookupOrAddToDict
|
||||
(
|
||||
"Cmu",
|
||||
this->coeffDict_,
|
||||
0.09
|
||||
)
|
||||
),
|
||||
C1_
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"C1",
|
||||
this->runTime_.timeName(),
|
||||
this->mesh_
|
||||
),
|
||||
this->mesh_,
|
||||
dimensioned<scalar>::lookupOrAddToDict
|
||||
(
|
||||
"C1",
|
||||
this->coeffDict_,
|
||||
1.44
|
||||
)
|
||||
),
|
||||
C2_
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"C2",
|
||||
this->runTime_.timeName(),
|
||||
this->mesh_
|
||||
),
|
||||
this->mesh_,
|
||||
dimensioned<scalar>::lookupOrAddToDict
|
||||
(
|
||||
"C2",
|
||||
this->coeffDict_,
|
||||
1.92
|
||||
)
|
||||
),
|
||||
sigmak_
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"sigmak",
|
||||
this->runTime_.timeName(),
|
||||
this->mesh_
|
||||
),
|
||||
this->mesh_,
|
||||
dimensioned<scalar>::lookupOrAddToDict
|
||||
(
|
||||
"sigmak",
|
||||
this->coeffDict_,
|
||||
1.0
|
||||
)
|
||||
),
|
||||
sigmaEps_
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"sigmaEps",
|
||||
this->runTime_.timeName(),
|
||||
this->mesh_
|
||||
),
|
||||
this->mesh_,
|
||||
dimensioned<scalar>::lookupOrAddToDict
|
||||
(
|
||||
"sigmaEps",
|
||||
this->coeffDict_,
|
||||
1.3
|
||||
)
|
||||
),
|
||||
|
||||
CdSigma_
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"CdSigma",
|
||||
this->runTime_.timeName(),
|
||||
this->mesh_
|
||||
),
|
||||
this->mesh_,
|
||||
dimensionedScalar("CdSigma", dimless/dimLength, 0)
|
||||
),
|
||||
betap_
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"betap",
|
||||
this->runTime_.timeName(),
|
||||
this->mesh_
|
||||
),
|
||||
this->mesh_,
|
||||
dimensionedScalar("betap", dimless, 0)
|
||||
),
|
||||
betad_
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"betad",
|
||||
this->runTime_.timeName(),
|
||||
this->mesh_
|
||||
),
|
||||
this->mesh_,
|
||||
dimensionedScalar("betad", dimless, 0)
|
||||
),
|
||||
C4_
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"C4",
|
||||
this->runTime_.timeName(),
|
||||
this->mesh_
|
||||
),
|
||||
this->mesh_,
|
||||
dimensionedScalar("C4", dimless, 0)
|
||||
),
|
||||
C5_
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"C5",
|
||||
this->runTime_.timeName(),
|
||||
this->mesh_
|
||||
),
|
||||
this->mesh_,
|
||||
dimensionedScalar("C5", dimless, 0)
|
||||
),
|
||||
|
||||
k_
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
IOobject::groupName("k", alphaRhoPhi.group()),
|
||||
this->runTime_.timeName(),
|
||||
this->mesh_,
|
||||
IOobject::MUST_READ,
|
||||
IOobject::AUTO_WRITE
|
||||
),
|
||||
this->mesh_
|
||||
),
|
||||
epsilon_
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
IOobject::groupName("epsilon", alphaRhoPhi.group()),
|
||||
this->runTime_.timeName(),
|
||||
this->mesh_,
|
||||
IOobject::MUST_READ,
|
||||
IOobject::AUTO_WRITE
|
||||
),
|
||||
this->mesh_
|
||||
)
|
||||
{
|
||||
bound(k_, this->kMin_);
|
||||
bound(epsilon_, this->epsilonMin_);
|
||||
|
||||
if (type == typeName)
|
||||
{
|
||||
this->printCoeffs(type);
|
||||
}
|
||||
|
||||
setPorosityCoefficients();
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
template<class BasicTurbulenceModel>
|
||||
bool kEpsilonLopesdaCosta<BasicTurbulenceModel>::read()
|
||||
{
|
||||
if (eddyViscosity<RASModel<BasicTurbulenceModel>>::read())
|
||||
{
|
||||
return true;
|
||||
}
|
||||
else
|
||||
{
|
||||
return false;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template<class BasicTurbulenceModel>
|
||||
void kEpsilonLopesdaCosta<BasicTurbulenceModel>::correct()
|
||||
{
|
||||
if (!this->turbulence_)
|
||||
{
|
||||
return;
|
||||
}
|
||||
|
||||
// Local references
|
||||
const alphaField& alpha = this->alpha_;
|
||||
const rhoField& rho = this->rho_;
|
||||
const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
|
||||
const volVectorField& U = this->U_;
|
||||
volScalarField& nut = this->nut_;
|
||||
fv::options& fvOptions(fv::options::New(this->mesh_));
|
||||
|
||||
eddyViscosity<RASModel<BasicTurbulenceModel>>::correct();
|
||||
|
||||
volScalarField::Internal divU
|
||||
(
|
||||
fvc::div(fvc::absolute(this->phi(), U))().v()
|
||||
);
|
||||
|
||||
tmp<volTensorField> tgradU = fvc::grad(U);
|
||||
volScalarField::Internal G
|
||||
(
|
||||
this->GName(),
|
||||
nut.v()*(dev(twoSymm(tgradU().v())) && tgradU().v())
|
||||
);
|
||||
tgradU.clear();
|
||||
|
||||
// Update epsilon and G at the wall
|
||||
epsilon_.boundaryFieldRef().updateCoeffs();
|
||||
|
||||
volScalarField::Internal magU(mag(U));
|
||||
volScalarField::Internal magU3(pow3(magU));
|
||||
|
||||
// Dissipation equation
|
||||
tmp<fvScalarMatrix> epsEqn
|
||||
(
|
||||
fvm::ddt(alpha, rho, epsilon_)
|
||||
+ fvm::div(alphaRhoPhi, epsilon_)
|
||||
- fvm::laplacian(alpha*rho*DepsilonEff(), epsilon_)
|
||||
==
|
||||
C1_*alpha()*rho()*G*epsilon_()/k_()
|
||||
- fvm::SuSp(((2.0/3.0)*C1_)*alpha()*rho()*divU, epsilon_)
|
||||
- fvm::Sp(C2_*alpha()*rho()*epsilon_()/k_(), epsilon_)
|
||||
+ epsilonSource(magU, magU3)
|
||||
+ fvOptions(alpha, rho, epsilon_)
|
||||
);
|
||||
|
||||
epsEqn.ref().relax();
|
||||
fvOptions.constrain(epsEqn.ref());
|
||||
epsEqn.ref().boundaryManipulate(epsilon_.boundaryFieldRef());
|
||||
solve(epsEqn);
|
||||
fvOptions.correct(epsilon_);
|
||||
bound(epsilon_, this->epsilonMin_);
|
||||
|
||||
// Turbulent kinetic energy equation
|
||||
tmp<fvScalarMatrix> kEqn
|
||||
(
|
||||
fvm::ddt(alpha, rho, k_)
|
||||
+ fvm::div(alphaRhoPhi, k_)
|
||||
- fvm::laplacian(alpha*rho*DkEff(), k_)
|
||||
==
|
||||
alpha()*rho()*G
|
||||
- fvm::SuSp((2.0/3.0)*alpha()*rho()*divU, k_)
|
||||
- fvm::Sp(alpha()*rho()*epsilon_()/k_(), k_)
|
||||
+ kSource(magU, magU3)
|
||||
+ fvOptions(alpha, rho, k_)
|
||||
);
|
||||
|
||||
kEqn.ref().relax();
|
||||
fvOptions.constrain(kEqn.ref());
|
||||
solve(kEqn);
|
||||
fvOptions.correct(k_);
|
||||
bound(k_, this->kMin_);
|
||||
|
||||
correctNut();
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
} // End namespace RASModels
|
||||
} // End namespace Foam
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -0,0 +1,245 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2018 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::RASModels::kEpsilonLopesdaCosta
|
||||
|
||||
Group
|
||||
grpRASTurbulence
|
||||
|
||||
Description
|
||||
Variant of the standard k-epsilon turbulence model with additional source
|
||||
terms to handle the changes in turbulence in porous regions represented by
|
||||
the powerLawLopesdaCosta porosity model.
|
||||
|
||||
Reference:
|
||||
\verbatim
|
||||
Costa, J. C. P. L. D. (2007).
|
||||
Atmospheric flow over forested and non-forested complex terrain.
|
||||
\endverbatim
|
||||
|
||||
The default model coefficients are
|
||||
\verbatim
|
||||
kEpsilonLopesdaCostaCoeffs
|
||||
{
|
||||
Cmu 0.09;
|
||||
C1 1.44;
|
||||
C2 1.92;
|
||||
sigmak 1.0;
|
||||
sigmaEps 1.3;
|
||||
}
|
||||
\endverbatim
|
||||
|
||||
See also
|
||||
Foam::RASModels::kEpsilon
|
||||
Foam::porosityModels::powerLawLopesdaCosta
|
||||
|
||||
SourceFiles
|
||||
kEpsilonLopesdaCosta.C
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#ifndef kEpsilonLopesdaCosta_H
|
||||
#define kEpsilonLopesdaCosta_H
|
||||
|
||||
#include "RASModel.H"
|
||||
#include "eddyViscosity.H"
|
||||
#include "powerLawLopesdaCosta.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
namespace RASModels
|
||||
{
|
||||
|
||||
/*---------------------------------------------------------------------------*\
|
||||
Class kEpsilonLopesdaCosta Declaration
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
template<class BasicTurbulenceModel>
|
||||
class kEpsilonLopesdaCosta
|
||||
:
|
||||
public eddyViscosity<RASModel<BasicTurbulenceModel>>
|
||||
{
|
||||
// Private Member Functions
|
||||
|
||||
// Disallow default bitwise copy construct and assignment
|
||||
kEpsilonLopesdaCosta(const kEpsilonLopesdaCosta&);
|
||||
void operator=(const kEpsilonLopesdaCosta&);
|
||||
|
||||
|
||||
protected:
|
||||
|
||||
// Protected data
|
||||
|
||||
// Standard model coefficients
|
||||
|
||||
volScalarField Cmu_;
|
||||
volScalarField::Internal C1_;
|
||||
volScalarField::Internal C2_;
|
||||
volScalarField sigmak_;
|
||||
volScalarField sigmaEps_;
|
||||
|
||||
// Lopes da Costa porosity coefficients
|
||||
|
||||
volScalarField::Internal CdSigma_;
|
||||
volScalarField::Internal betap_;
|
||||
volScalarField::Internal betad_;
|
||||
volScalarField::Internal C4_;
|
||||
volScalarField::Internal C5_;
|
||||
|
||||
|
||||
// Fields
|
||||
|
||||
volScalarField k_;
|
||||
volScalarField epsilon_;
|
||||
|
||||
|
||||
// Protected Member Functions
|
||||
|
||||
void setPorosityCoefficient
|
||||
(
|
||||
volScalarField::Internal& C,
|
||||
const porosityModels::powerLawLopesdaCosta& pm
|
||||
);
|
||||
|
||||
void setCdSigma
|
||||
(
|
||||
volScalarField::Internal& C,
|
||||
const porosityModels::powerLawLopesdaCosta& pm
|
||||
);
|
||||
|
||||
void setPorosityCoefficients();
|
||||
|
||||
virtual void correctNut();
|
||||
|
||||
virtual tmp<fvScalarMatrix> kSource
|
||||
(
|
||||
const volScalarField::Internal& magU,
|
||||
const volScalarField::Internal& magU3
|
||||
) const;
|
||||
|
||||
virtual tmp<fvScalarMatrix> epsilonSource
|
||||
(
|
||||
const volScalarField::Internal& magU,
|
||||
const volScalarField::Internal& magU3
|
||||
) const;
|
||||
|
||||
|
||||
public:
|
||||
|
||||
typedef typename BasicTurbulenceModel::alphaField alphaField;
|
||||
typedef typename BasicTurbulenceModel::rhoField rhoField;
|
||||
typedef typename BasicTurbulenceModel::transportModel transportModel;
|
||||
|
||||
|
||||
//- Runtime type information
|
||||
TypeName("kEpsilonLopesdaCosta");
|
||||
|
||||
|
||||
// Constructors
|
||||
|
||||
//- Construct from components
|
||||
kEpsilonLopesdaCosta
|
||||
(
|
||||
const alphaField& alpha,
|
||||
const rhoField& rho,
|
||||
const volVectorField& U,
|
||||
const surfaceScalarField& alphaRhoPhi,
|
||||
const surfaceScalarField& phi,
|
||||
const transportModel& transport,
|
||||
const word& propertiesName = turbulenceModel::propertiesName,
|
||||
const word& type = typeName
|
||||
);
|
||||
|
||||
|
||||
//- Destructor
|
||||
virtual ~kEpsilonLopesdaCosta()
|
||||
{}
|
||||
|
||||
|
||||
// Member Functions
|
||||
|
||||
//- Re-read model coefficients if they have changed
|
||||
virtual bool read();
|
||||
|
||||
//- Return the effective diffusivity for k
|
||||
tmp<volScalarField> DkEff() const
|
||||
{
|
||||
return tmp<volScalarField>
|
||||
(
|
||||
new volScalarField
|
||||
(
|
||||
"DkEff",
|
||||
(this->nut_/sigmak_ + this->nu())
|
||||
)
|
||||
);
|
||||
}
|
||||
|
||||
//- Return the effective diffusivity for epsilon
|
||||
tmp<volScalarField> DepsilonEff() const
|
||||
{
|
||||
return tmp<volScalarField>
|
||||
(
|
||||
new volScalarField
|
||||
(
|
||||
"DepsilonEff",
|
||||
(this->nut_/sigmaEps_ + this->nu())
|
||||
)
|
||||
);
|
||||
}
|
||||
|
||||
//- Return the turbulence kinetic energy
|
||||
virtual tmp<volScalarField> k() const
|
||||
{
|
||||
return k_;
|
||||
}
|
||||
|
||||
//- Return the turbulence kinetic energy dissipation rate
|
||||
virtual tmp<volScalarField> epsilon() const
|
||||
{
|
||||
return epsilon_;
|
||||
}
|
||||
|
||||
//- Solve the turbulence equations and correct the turbulence viscosity
|
||||
virtual void correct();
|
||||
};
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
} // End namespace RASModels
|
||||
} // End namespace Foam
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#ifdef NoRepository
|
||||
#include "kEpsilonLopesdaCosta.C"
|
||||
#endif
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#endif
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -0,0 +1,419 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2018 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 "addToRunTimeSelectionTable.H"
|
||||
#include "powerLawLopesdaCosta.H"
|
||||
#include "geometricOneField.H"
|
||||
#include "fvMatrices.H"
|
||||
#include "triSurfaceMesh.H"
|
||||
#include "triSurfaceTools.H"
|
||||
|
||||
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
namespace porosityModels
|
||||
{
|
||||
defineTypeNameAndDebug(powerLawLopesdaCosta, 0);
|
||||
addToRunTimeSelectionTable(porosityModel, powerLawLopesdaCosta, mesh);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
|
||||
|
||||
Foam::porosityModels::powerLawLopesdaCostaZone::powerLawLopesdaCostaZone
|
||||
(
|
||||
const word& name,
|
||||
const word& modelType,
|
||||
const fvMesh& mesh,
|
||||
const dictionary& dict
|
||||
)
|
||||
:
|
||||
zoneName_(name + ":porousZone")
|
||||
{
|
||||
dictionary coeffs(dict.optionalSubDict(modelType + "Coeffs"));
|
||||
|
||||
// Vertical direction within porous region
|
||||
vector zDir(coeffs.lookup("zDir"));
|
||||
|
||||
// Span of the search for the cells in the porous region
|
||||
scalar searchSpan(readScalar(coeffs.lookup("searchSpan")));
|
||||
|
||||
// Top surface file name defining the extent of the porous region
|
||||
word topSurfaceFileName(coeffs.lookup("topSurface"));
|
||||
|
||||
// List of the ground patches defining the lower surface
|
||||
// of the porous region
|
||||
List<wordRe> groundPatches(coeffs.lookup("groundPatches"));
|
||||
|
||||
// Functional form of the porosity surface area per unit volume as a
|
||||
// function of the normalized vertical position
|
||||
autoPtr<Function1<scalar >> SigmaFunc
|
||||
(
|
||||
Function1<scalar>::New("Sigma", dict)
|
||||
);
|
||||
|
||||
// Searchable triSurface for the top of the porous region
|
||||
triSurfaceMesh searchSurf
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
topSurfaceFileName,
|
||||
mesh.time().constant(),
|
||||
"triSurface",
|
||||
mesh.time()
|
||||
)
|
||||
);
|
||||
|
||||
// Limit the porous cell search to those within the lateral and vertical
|
||||
// extent of the top surface
|
||||
|
||||
boundBox surfBounds(searchSurf.points());
|
||||
boundBox searchBounds
|
||||
(
|
||||
surfBounds.min() - searchSpan*zDir, surfBounds.max()
|
||||
);
|
||||
|
||||
const pointField& C = mesh.cellCentres();
|
||||
|
||||
// List of cells within the porous region
|
||||
labelList porousCells(C.size());
|
||||
label porousCelli = 0;
|
||||
|
||||
forAll(C, celli)
|
||||
{
|
||||
if (searchBounds.contains(C[celli]))
|
||||
{
|
||||
porousCells[porousCelli++] = celli;
|
||||
}
|
||||
}
|
||||
|
||||
porousCells.setSize(porousCelli);
|
||||
|
||||
pointField start(porousCelli);
|
||||
pointField end(porousCelli);
|
||||
|
||||
forAll(porousCells, porousCelli)
|
||||
{
|
||||
start[porousCelli] = C[porousCells[porousCelli]];
|
||||
end[porousCelli] = start[porousCelli] + searchSpan*zDir;
|
||||
}
|
||||
|
||||
// Field of distance between the cell centre and the corresponding point
|
||||
// on the porous region top surface
|
||||
scalarField zTop(porousCelli);
|
||||
|
||||
// Search the porous cells for those with a corresponding point on the
|
||||
// porous region top surface
|
||||
List<pointIndexHit> hitInfo;
|
||||
searchSurf.findLine(start, end, hitInfo);
|
||||
|
||||
porousCelli = 0;
|
||||
|
||||
forAll(porousCells, celli)
|
||||
{
|
||||
const pointIndexHit& hit = hitInfo[celli];
|
||||
|
||||
if (hit.hit())
|
||||
{
|
||||
porousCells[porousCelli] = porousCells[celli];
|
||||
|
||||
zTop[porousCelli] =
|
||||
(hit.hitPoint() - C[porousCells[porousCelli]]) & zDir;
|
||||
|
||||
porousCelli++;
|
||||
}
|
||||
}
|
||||
|
||||
// Resize the porous cells list and height field
|
||||
porousCells.setSize(porousCelli);
|
||||
zTop.setSize(porousCelli);
|
||||
|
||||
// Create a triSurface for the ground patch(s)
|
||||
triSurface groundSurface
|
||||
(
|
||||
triSurfaceTools::triangulate
|
||||
(
|
||||
mesh.boundaryMesh(),
|
||||
mesh.boundaryMesh().patchSet(groundPatches),
|
||||
searchBounds
|
||||
)
|
||||
);
|
||||
|
||||
// Combine the ground triSurfaces across all processors
|
||||
if (Pstream::parRun())
|
||||
{
|
||||
List<List<labelledTri>> groundSurfaceProcTris(Pstream::nProcs());
|
||||
List<pointField> groundSurfaceProcPoints(Pstream::nProcs());
|
||||
|
||||
groundSurfaceProcTris[Pstream::myProcNo()] = groundSurface;
|
||||
groundSurfaceProcPoints[Pstream::myProcNo()] = groundSurface.points();
|
||||
|
||||
Pstream::gatherList(groundSurfaceProcTris);
|
||||
Pstream::scatterList(groundSurfaceProcTris);
|
||||
Pstream::gatherList(groundSurfaceProcPoints);
|
||||
Pstream::scatterList(groundSurfaceProcPoints);
|
||||
|
||||
label nTris = 0;
|
||||
forAll(groundSurfaceProcTris, i)
|
||||
{
|
||||
nTris += groundSurfaceProcTris[i].size();
|
||||
}
|
||||
|
||||
List<labelledTri> groundSurfaceTris(nTris);
|
||||
label trii = 0;
|
||||
label offset = 0;
|
||||
forAll(groundSurfaceProcTris, i)
|
||||
{
|
||||
forAll(groundSurfaceProcTris[i], j)
|
||||
{
|
||||
groundSurfaceTris[trii] = groundSurfaceProcTris[i][j];
|
||||
groundSurfaceTris[trii][0] += offset;
|
||||
groundSurfaceTris[trii][1] += offset;
|
||||
groundSurfaceTris[trii][2] += offset;
|
||||
trii++;
|
||||
}
|
||||
offset += groundSurfaceProcPoints[i].size();
|
||||
}
|
||||
|
||||
label nPoints = 0;
|
||||
forAll(groundSurfaceProcPoints, i)
|
||||
{
|
||||
nPoints += groundSurfaceProcPoints[i].size();
|
||||
}
|
||||
|
||||
pointField groundSurfacePoints(nPoints);
|
||||
|
||||
label pointi = 0;
|
||||
forAll(groundSurfaceProcPoints, i)
|
||||
{
|
||||
forAll(groundSurfaceProcPoints[i], j)
|
||||
{
|
||||
groundSurfacePoints[pointi++] = groundSurfaceProcPoints[i][j];
|
||||
}
|
||||
}
|
||||
|
||||
groundSurface = triSurface(groundSurfaceTris, groundSurfacePoints);
|
||||
}
|
||||
|
||||
// Create a searchable triSurface for the ground
|
||||
triSurfaceSearch groundSearch(groundSurface);
|
||||
|
||||
// Search the porous cells for the corresponding point on the ground
|
||||
|
||||
start.setSize(porousCelli);
|
||||
end.setSize(porousCelli);
|
||||
|
||||
forAll(porousCells, porousCelli)
|
||||
{
|
||||
start[porousCelli] = C[porousCells[porousCelli]];
|
||||
end[porousCelli] = start[porousCelli] - searchSpan*zDir;
|
||||
}
|
||||
|
||||
groundSearch.findLine(start, end, hitInfo);
|
||||
|
||||
scalarField zBottom(porousCelli);
|
||||
|
||||
forAll(porousCells, porousCelli)
|
||||
{
|
||||
const pointIndexHit& hit = hitInfo[porousCelli];
|
||||
|
||||
if (hit.hit())
|
||||
{
|
||||
zBottom[porousCelli] =
|
||||
(C[porousCells[porousCelli]] - hit.hitPoint()) & zDir;
|
||||
}
|
||||
}
|
||||
|
||||
// Create the normalized height field
|
||||
scalarField zNorm(zBottom/(zBottom + zTop));
|
||||
|
||||
// Create the porosity surface area per unit volume zone field
|
||||
Sigma_ = SigmaFunc->value(zNorm);
|
||||
|
||||
// Create the porous region cellZone and add to the mesh cellZones
|
||||
|
||||
cellZoneMesh& cellZones = const_cast<cellZoneMesh&>(mesh.cellZones());
|
||||
|
||||
label zoneID = cellZones.findZoneID(zoneName_);
|
||||
|
||||
if (zoneID == -1)
|
||||
{
|
||||
zoneID = cellZones.size();
|
||||
cellZones.setSize(zoneID + 1);
|
||||
|
||||
cellZones.set
|
||||
(
|
||||
zoneID,
|
||||
new cellZone
|
||||
(
|
||||
zoneName_,
|
||||
porousCells,
|
||||
zoneID,
|
||||
cellZones
|
||||
)
|
||||
);
|
||||
}
|
||||
else
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "Unable to create porous cellZone " << zoneName_
|
||||
<< ": zone already exists"
|
||||
<< abort(FatalError);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
Foam::porosityModels::powerLawLopesdaCosta::powerLawLopesdaCosta
|
||||
(
|
||||
const word& name,
|
||||
const word& modelType,
|
||||
const fvMesh& mesh,
|
||||
const dictionary& dict,
|
||||
const word& dummyCellZoneName
|
||||
)
|
||||
:
|
||||
powerLawLopesdaCostaZone(name, modelType, mesh, dict),
|
||||
porosityModel
|
||||
(
|
||||
name,
|
||||
modelType,
|
||||
mesh,
|
||||
dict,
|
||||
powerLawLopesdaCostaZone::zoneName_
|
||||
),
|
||||
Cd_(readScalar(coeffs_.lookup("Cd"))),
|
||||
C1_(readScalar(coeffs_.lookup("C1"))),
|
||||
rhoName_(coeffs_.lookupOrDefault<word>("rho", "rho"))
|
||||
{}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
|
||||
|
||||
Foam::porosityModels::powerLawLopesdaCosta::~powerLawLopesdaCosta()
|
||||
{}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
const Foam::scalarField&
|
||||
Foam::porosityModels::powerLawLopesdaCostaZone::Sigma() const
|
||||
{
|
||||
return Sigma_;
|
||||
}
|
||||
|
||||
|
||||
void Foam::porosityModels::powerLawLopesdaCosta::calcTransformModelData()
|
||||
{}
|
||||
|
||||
|
||||
void Foam::porosityModels::powerLawLopesdaCosta::calcForce
|
||||
(
|
||||
const volVectorField& U,
|
||||
const volScalarField& rho,
|
||||
const volScalarField& mu,
|
||||
vectorField& force
|
||||
) const
|
||||
{
|
||||
scalarField Udiag(U.size(), 0.0);
|
||||
const scalarField& V = mesh_.V();
|
||||
|
||||
apply(Udiag, V, rho, U);
|
||||
|
||||
force = Udiag*U;
|
||||
}
|
||||
|
||||
|
||||
void Foam::porosityModels::powerLawLopesdaCosta::correct
|
||||
(
|
||||
fvVectorMatrix& UEqn
|
||||
) const
|
||||
{
|
||||
const vectorField& U = UEqn.psi();
|
||||
const scalarField& V = mesh_.V();
|
||||
scalarField& Udiag = UEqn.diag();
|
||||
|
||||
if (UEqn.dimensions() == dimForce)
|
||||
{
|
||||
const volScalarField& rho =
|
||||
mesh_.lookupObject<volScalarField>(rhoName_);
|
||||
|
||||
apply(Udiag, V, rho, U);
|
||||
}
|
||||
else
|
||||
{
|
||||
apply(Udiag, V, geometricOneField(), U);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void Foam::porosityModels::powerLawLopesdaCosta::correct
|
||||
(
|
||||
fvVectorMatrix& UEqn,
|
||||
const volScalarField& rho,
|
||||
const volScalarField& mu
|
||||
) const
|
||||
{
|
||||
const vectorField& U = UEqn.psi();
|
||||
const scalarField& V = mesh_.V();
|
||||
scalarField& Udiag = UEqn.diag();
|
||||
|
||||
apply(Udiag, V, rho, U);
|
||||
}
|
||||
|
||||
|
||||
void Foam::porosityModels::powerLawLopesdaCosta::correct
|
||||
(
|
||||
const fvVectorMatrix& UEqn,
|
||||
volTensorField& AU
|
||||
) const
|
||||
{
|
||||
const vectorField& U = UEqn.psi();
|
||||
|
||||
if (UEqn.dimensions() == dimForce)
|
||||
{
|
||||
const volScalarField& rho =
|
||||
mesh_.lookupObject<volScalarField>(rhoName_);
|
||||
|
||||
apply(AU, rho, U);
|
||||
}
|
||||
else
|
||||
{
|
||||
apply(AU, geometricOneField(), U);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
bool Foam::porosityModels::powerLawLopesdaCosta::writeData(Ostream& os) const
|
||||
{
|
||||
os << indent << name_ << endl;
|
||||
dict_.write(os);
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -0,0 +1,229 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2018 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::powerLawLopesdaCosta
|
||||
|
||||
Description
|
||||
Variant of the power law porosity model with spatially varying
|
||||
drag coefficient
|
||||
|
||||
given by:
|
||||
|
||||
\f[
|
||||
S = -\rho C_d \Sigma |U|^{(C_1 - 1)} U
|
||||
\f]
|
||||
|
||||
where
|
||||
\vartable
|
||||
\Sigma | Porosity surface area per unit volume
|
||||
C_d | Model linear coefficient
|
||||
C_1 | Model exponent coefficient
|
||||
\endvartable
|
||||
|
||||
Reference:
|
||||
\verbatim
|
||||
Costa, J. C. P. L. D. (2007).
|
||||
Atmospheric flow over forested and non-forested complex terrain.
|
||||
\endverbatim
|
||||
|
||||
See also
|
||||
Foam::RASModels::kEpsilonLopesdaCosta
|
||||
|
||||
SourceFiles
|
||||
powerLawLopesdaCosta.C
|
||||
powerLawLopesdaCostaTemplates.C
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#ifndef powerLawLopesdaCosta_H
|
||||
#define powerLawLopesdaCosta_H
|
||||
|
||||
#include "porosityModel.H"
|
||||
#include "Function1.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
namespace porosityModels
|
||||
{
|
||||
|
||||
/*---------------------------------------------------------------------------*\
|
||||
Class powerLawLopesdaCostaZone Declaration
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
class powerLawLopesdaCostaZone
|
||||
{
|
||||
protected:
|
||||
|
||||
// Protected data
|
||||
|
||||
//- Automatically generated zone name for this porous zone
|
||||
const word zoneName_;
|
||||
|
||||
//- Porosity surface area per unit volume zone field
|
||||
scalarField Sigma_;
|
||||
|
||||
|
||||
public:
|
||||
|
||||
//- Constructor
|
||||
powerLawLopesdaCostaZone
|
||||
(
|
||||
const word& name,
|
||||
const word& modelType,
|
||||
const fvMesh& mesh,
|
||||
const dictionary& dict
|
||||
);
|
||||
|
||||
// Member Functions
|
||||
|
||||
//- Return the porosity surface area per unit volume zone field
|
||||
const scalarField& Sigma() const;
|
||||
};
|
||||
|
||||
|
||||
/*---------------------------------------------------------------------------*\
|
||||
Class powerLawLopesdaCosta Declaration
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
class powerLawLopesdaCosta
|
||||
:
|
||||
public powerLawLopesdaCostaZone,
|
||||
public porosityModel
|
||||
{
|
||||
// Private data
|
||||
|
||||
//- Cd coefficient
|
||||
scalar Cd_;
|
||||
|
||||
//- C1 coefficient
|
||||
scalar C1_;
|
||||
|
||||
//- Name of density field
|
||||
word rhoName_;
|
||||
|
||||
|
||||
// Private Member Functions
|
||||
|
||||
//- Apply resistance
|
||||
template<class RhoFieldType>
|
||||
void apply
|
||||
(
|
||||
scalarField& Udiag,
|
||||
const scalarField& V,
|
||||
const RhoFieldType& rho,
|
||||
const vectorField& U
|
||||
) const;
|
||||
|
||||
//- Apply resistance
|
||||
template<class RhoFieldType>
|
||||
void apply
|
||||
(
|
||||
tensorField& AU,
|
||||
const RhoFieldType& rho,
|
||||
const vectorField& U
|
||||
) const;
|
||||
|
||||
//- Disallow default bitwise copy construct
|
||||
powerLawLopesdaCosta(const powerLawLopesdaCosta&);
|
||||
|
||||
//- Disallow default bitwise assignment
|
||||
void operator=(const powerLawLopesdaCosta&);
|
||||
|
||||
|
||||
public:
|
||||
|
||||
//- Runtime type information
|
||||
TypeName("powerLawLopesdaCosta");
|
||||
|
||||
//- Constructor
|
||||
powerLawLopesdaCosta
|
||||
(
|
||||
const word& name,
|
||||
const word& modelType,
|
||||
const fvMesh& mesh,
|
||||
const dictionary& dict,
|
||||
const word& cellZoneName
|
||||
);
|
||||
|
||||
//- Destructor
|
||||
virtual ~powerLawLopesdaCosta();
|
||||
|
||||
|
||||
// Member Functions
|
||||
|
||||
//- Transform the model data wrt mesh changes
|
||||
virtual void calcTransformModelData();
|
||||
|
||||
//- Calculate the porosity force
|
||||
virtual void calcForce
|
||||
(
|
||||
const volVectorField& U,
|
||||
const volScalarField& rho,
|
||||
const volScalarField& mu,
|
||||
vectorField& force
|
||||
) const;
|
||||
|
||||
//- Add resistance
|
||||
virtual void correct(fvVectorMatrix& UEqn) const;
|
||||
|
||||
//- Add resistance
|
||||
virtual void correct
|
||||
(
|
||||
fvVectorMatrix& UEqn,
|
||||
const volScalarField& rho,
|
||||
const volScalarField& mu
|
||||
) const;
|
||||
|
||||
//- Add resistance
|
||||
virtual void correct
|
||||
(
|
||||
const fvVectorMatrix& UEqn,
|
||||
volTensorField& AU
|
||||
) const;
|
||||
|
||||
|
||||
// I-O
|
||||
|
||||
//- Write
|
||||
bool writeData(Ostream& os) const;
|
||||
};
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
} // End namespace porosityModels
|
||||
} // End namespace Foam
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#ifdef NoRepository
|
||||
#include "powerLawLopesdaCostaTemplates.C"
|
||||
#endif
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#endif
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -0,0 +1,87 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2018 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 "volFields.H"
|
||||
|
||||
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
||||
|
||||
template<class RhoFieldType>
|
||||
void Foam::porosityModels::powerLawLopesdaCosta::apply
|
||||
(
|
||||
scalarField& Udiag,
|
||||
const scalarField& V,
|
||||
const RhoFieldType& rho,
|
||||
const vectorField& U
|
||||
) const
|
||||
{
|
||||
const scalar C1m1b2 = (C1_ - 1.0)/2.0;
|
||||
|
||||
forAll(cellZoneIDs_, zonei)
|
||||
{
|
||||
const labelList& cells = mesh_.cellZones()[cellZoneIDs_[zonei]];
|
||||
|
||||
forAll(cells, i)
|
||||
{
|
||||
const label celli = cells[i];
|
||||
|
||||
Udiag[celli] +=
|
||||
V[celli]*rho[celli]
|
||||
*Cd_*Sigma_[i]*pow(magSqr(U[celli]), C1m1b2);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template<class RhoFieldType>
|
||||
void Foam::porosityModels::powerLawLopesdaCosta::apply
|
||||
(
|
||||
tensorField& AU,
|
||||
const RhoFieldType& rho,
|
||||
const vectorField& U
|
||||
) const
|
||||
{
|
||||
const scalar C1m1b2 = (C1_ - 1.0)/2.0;
|
||||
|
||||
forAll(cellZoneIDs_, zonei)
|
||||
{
|
||||
const labelList& cells = mesh_.cellZones()[cellZoneIDs_[zonei]];
|
||||
|
||||
forAll(cells, i)
|
||||
{
|
||||
const label celli = cells[i];
|
||||
|
||||
AU[celli] =
|
||||
AU[celli]
|
||||
+ I
|
||||
*(
|
||||
0.5*rho[celli]*Cd_*Sigma_[i]
|
||||
*pow(magSqr(U[celli]), C1m1b2)
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -2,7 +2,7 @@
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2012-2015 OpenFOAM Foundation
|
||||
\\ / A nd | Copyright (C) 2012-2018 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
@ -223,6 +223,9 @@ public:
|
||||
//- Return const access to the cell zone IDs
|
||||
inline const labelList& cellZoneIDs() const;
|
||||
|
||||
//- Return dictionary used for model construction
|
||||
const dictionary& dict() const;
|
||||
|
||||
//- Transform the model data wrt mesh changes
|
||||
virtual void transformModelData();
|
||||
|
||||
|
||||
@ -2,7 +2,7 @@
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation
|
||||
\\ / A nd | Copyright (C) 2012-2018 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
@ -35,6 +35,12 @@ inline bool Foam::porosityModel::active() const
|
||||
}
|
||||
|
||||
|
||||
inline const Foam::dictionary& Foam::porosityModel::dict() const
|
||||
{
|
||||
return dict_;
|
||||
}
|
||||
|
||||
|
||||
inline const Foam::labelList& Foam::porosityModel::cellZoneIDs() const
|
||||
{
|
||||
return cellZoneIDs_;
|
||||
|
||||
@ -2,7 +2,7 @@
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2012-2016 OpenFOAM Foundation
|
||||
\\ / A nd | Copyright (C) 2012-2018 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
@ -61,13 +61,6 @@ Foam::fv::explicitPorositySource::explicitPorositySource
|
||||
{
|
||||
read(dict);
|
||||
|
||||
if (selectionMode_ != smCellZone)
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "selection mode is " << selectionModeTypeNames_[selectionMode_]
|
||||
<< exit(FatalError);
|
||||
}
|
||||
|
||||
porosityPtr_.reset
|
||||
(
|
||||
porosityModel::New
|
||||
|
||||
@ -2,7 +2,7 @@
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2012-2016 OpenFOAM Foundation
|
||||
\\ / A nd | Copyright (C) 2012-2018 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
@ -131,6 +131,11 @@ public:
|
||||
|
||||
// Member Functions
|
||||
|
||||
const porosityModel& model() const
|
||||
{
|
||||
return porosityPtr_();
|
||||
}
|
||||
|
||||
// Add explicit and implicit contributions
|
||||
|
||||
//- Add implicit contribution to momentum equation
|
||||
|
||||
@ -32,7 +32,6 @@ License
|
||||
#include "plane.H"
|
||||
#include "geompack.H"
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
|
||||
|
||||
const Foam::label Foam::triSurfaceTools::ANYEDGE = -1;
|
||||
@ -2328,7 +2327,6 @@ Foam::triSurfaceTools::sideType Foam::triSurfaceTools::surfaceSide
|
||||
}
|
||||
|
||||
|
||||
// triangulation of boundaryMesh
|
||||
Foam::triSurface Foam::triSurfaceTools::triangulate
|
||||
(
|
||||
const polyBoundaryMesh& bMesh,
|
||||
@ -2415,6 +2413,96 @@ Foam::triSurface Foam::triSurfaceTools::triangulate
|
||||
}
|
||||
|
||||
|
||||
Foam::triSurface Foam::triSurfaceTools::triangulate
|
||||
(
|
||||
const polyBoundaryMesh& bMesh,
|
||||
const labelHashSet& includePatches,
|
||||
const boundBox& bBox,
|
||||
const bool verbose
|
||||
)
|
||||
{
|
||||
const polyMesh& mesh = bMesh.mesh();
|
||||
|
||||
// Storage for surfaceMesh. Size estimate.
|
||||
DynamicList<labelledTri> triangles
|
||||
(
|
||||
mesh.nFaces() - mesh.nInternalFaces()
|
||||
);
|
||||
|
||||
label newPatchi = 0;
|
||||
|
||||
forAllConstIter(labelHashSet, includePatches, iter)
|
||||
{
|
||||
const label patchi = iter.key();
|
||||
const polyPatch& patch = bMesh[patchi];
|
||||
const pointField& points = patch.points();
|
||||
|
||||
label nTriTotal = 0;
|
||||
|
||||
forAll(patch, patchFacei)
|
||||
{
|
||||
const face& f = patch[patchFacei];
|
||||
|
||||
if (bBox.containsAny(points, f))
|
||||
{
|
||||
faceList triFaces(f.nTriangles(points));
|
||||
|
||||
label nTri = 0;
|
||||
|
||||
f.triangles(points, nTri, triFaces);
|
||||
|
||||
forAll(triFaces, triFacei)
|
||||
{
|
||||
const face& f = triFaces[triFacei];
|
||||
|
||||
triangles.append(labelledTri(f[0], f[1], f[2], newPatchi));
|
||||
|
||||
nTriTotal++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (verbose)
|
||||
{
|
||||
Pout<< patch.name() << " : generated " << nTriTotal
|
||||
<< " triangles from " << patch.size() << " faces with"
|
||||
<< " new patchid " << newPatchi << endl;
|
||||
}
|
||||
|
||||
newPatchi++;
|
||||
}
|
||||
triangles.shrink();
|
||||
|
||||
// Create globally numbered tri surface
|
||||
triSurface rawSurface(triangles, mesh.points());
|
||||
|
||||
// Create locally numbered tri surface
|
||||
triSurface surface
|
||||
(
|
||||
rawSurface.localFaces(),
|
||||
rawSurface.localPoints()
|
||||
);
|
||||
|
||||
// Add patch names to surface
|
||||
surface.patches().setSize(newPatchi);
|
||||
|
||||
newPatchi = 0;
|
||||
|
||||
forAllConstIter(labelHashSet, includePatches, iter)
|
||||
{
|
||||
const label patchi = iter.key();
|
||||
const polyPatch& patch = bMesh[patchi];
|
||||
|
||||
surface.patches()[newPatchi].name() = patch.name();
|
||||
surface.patches()[newPatchi].geometricType() = patch.type();
|
||||
|
||||
newPatchi++;
|
||||
}
|
||||
|
||||
return surface;
|
||||
}
|
||||
|
||||
|
||||
// triangulation of boundaryMesh
|
||||
Foam::triSurface Foam::triSurfaceTools::triangulateFaceCentre
|
||||
(
|
||||
|
||||
@ -2,7 +2,7 @@
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
|
||||
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
@ -55,6 +55,7 @@ class edge;
|
||||
class labelledTri;
|
||||
class polyBoundaryMesh;
|
||||
class plane;
|
||||
class boundBox;
|
||||
|
||||
/*---------------------------------------------------------------------------*\
|
||||
Class triSurfaceTools Declaration
|
||||
@ -473,6 +474,16 @@ public:
|
||||
const bool verbose = false
|
||||
);
|
||||
|
||||
|
||||
static triSurface triangulate
|
||||
(
|
||||
const polyBoundaryMesh& bMesh,
|
||||
const labelHashSet& includePatches,
|
||||
const boundBox& bBox,
|
||||
const bool verbose = false
|
||||
);
|
||||
|
||||
|
||||
//- Face-centre triangulation of (selected patches of) boundaryMesh.
|
||||
// Needs
|
||||
// polyMesh (or polyBoundaryMesh) since only at this level are the
|
||||
|
||||
Reference in New Issue
Block a user