diff --git a/src/Allwmake b/src/Allwmake index 92a389d77e..300792371c 100755 --- a/src/Allwmake +++ b/src/Allwmake @@ -82,6 +82,7 @@ wmake $targetType sixDoFRigidBodyState wmake $targetType rigidBodyDynamics wmake $targetType rigidBodyMeshMotion wmake $targetType semiPermeableBaffle +wmake $targetType atmosphericModels # Needs access to Turbulence diff --git a/src/atmosphericModels/Make/files b/src/atmosphericModels/Make/files index a5a1c27f22..644149f555 100644 --- a/src/atmosphericModels/Make/files +++ b/src/atmosphericModels/Make/files @@ -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 diff --git a/src/atmosphericModels/atmosphericTurbulentTransportModels.C b/src/atmosphericModels/atmosphericTurbulentTransportModels.C new file mode 100644 index 0000000000..cf31e0c250 --- /dev/null +++ b/src/atmosphericModels/atmosphericTurbulentTransportModels.C @@ -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 . + +\*---------------------------------------------------------------------------*/ + +#include "turbulentTransportModels.H" + +// -------------------------------------------------------------------------- // +// RAS models +// -------------------------------------------------------------------------- // + +#include "kEpsilonLopesdaCosta.H" +makeRASModel(kEpsilonLopesdaCosta); + +// ************************************************************************* // diff --git a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.H b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.H index 9c37660b49..a7360f8d38 100644 --- a/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.H +++ b/src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.H @@ -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 | Copyright (C) 2015 OpenCFD Ltd. ------------------------------------------------------------------------------- License diff --git a/src/atmosphericModels/kEpsilonLopesdaCosta/kEpsilonLopesdaCosta.C b/src/atmosphericModels/kEpsilonLopesdaCosta/kEpsilonLopesdaCosta.C new file mode 100644 index 0000000000..805d0ca16b --- /dev/null +++ b/src/atmosphericModels/kEpsilonLopesdaCosta/kEpsilonLopesdaCosta.C @@ -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 . + +\*---------------------------------------------------------------------------*/ + +#include "kEpsilonLopesdaCosta.H" +#include "fvOptions.H" +#include "explicitPorositySource.H" +#include "bound.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace RASModels +{ + +// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // + +template +void kEpsilonLopesdaCosta::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 +void kEpsilonLopesdaCosta::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 +void kEpsilonLopesdaCosta::setPorosityCoefficients() +{ + fv::options::optionList& fvOptions(fv::options::New(this->mesh_)); + + forAll(fvOptions, i) + { + if (isA(fvOptions[i])) + { + const fv::explicitPorositySource& eps = + refCast(fvOptions[i]); + + if (isA(eps.model())) + { + const porosityModels::powerLawLopesdaCosta& pm = + refCast + ( + 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 +void kEpsilonLopesdaCosta::correctNut() +{ + this->nut_ = Cmu_*sqr(k_)/epsilon_; + this->nut_.correctBoundaryConditions(); + fv::options::New(this->mesh_).correct(this->nut_); + + BasicTurbulenceModel::correctNut(); +} + + +template +tmp kEpsilonLopesdaCosta::kSource +( + const volScalarField::Internal& magU, + const volScalarField::Internal& magU3 +) const +{ + return fvm::Su(CdSigma_*(betap_*magU3 - betad_*magU*k_()), k_); +} + + +template +tmp +kEpsilonLopesdaCosta::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 +kEpsilonLopesdaCosta::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> + ( + type, + alpha, + rho, + U, + alphaRhoPhi, + phi, + transport, + propertiesName + ), + + Cmu_ + ( + IOobject + ( + "Cmu", + this->runTime_.timeName(), + this->mesh_ + ), + this->mesh_, + dimensioned::lookupOrAddToDict + ( + "Cmu", + this->coeffDict_, + 0.09 + ) + ), + C1_ + ( + IOobject + ( + "C1", + this->runTime_.timeName(), + this->mesh_ + ), + this->mesh_, + dimensioned::lookupOrAddToDict + ( + "C1", + this->coeffDict_, + 1.44 + ) + ), + C2_ + ( + IOobject + ( + "C2", + this->runTime_.timeName(), + this->mesh_ + ), + this->mesh_, + dimensioned::lookupOrAddToDict + ( + "C2", + this->coeffDict_, + 1.92 + ) + ), + sigmak_ + ( + IOobject + ( + "sigmak", + this->runTime_.timeName(), + this->mesh_ + ), + this->mesh_, + dimensioned::lookupOrAddToDict + ( + "sigmak", + this->coeffDict_, + 1.0 + ) + ), + sigmaEps_ + ( + IOobject + ( + "sigmaEps", + this->runTime_.timeName(), + this->mesh_ + ), + this->mesh_, + dimensioned::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 +bool kEpsilonLopesdaCosta::read() +{ + if (eddyViscosity>::read()) + { + return true; + } + else + { + return false; + } +} + + +template +void kEpsilonLopesdaCosta::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>::correct(); + + volScalarField::Internal divU + ( + fvc::div(fvc::absolute(this->phi(), U))().v() + ); + + tmp 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 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 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 + +// ************************************************************************* // diff --git a/src/atmosphericModels/kEpsilonLopesdaCosta/kEpsilonLopesdaCosta.H b/src/atmosphericModels/kEpsilonLopesdaCosta/kEpsilonLopesdaCosta.H new file mode 100644 index 0000000000..b72992d9a4 --- /dev/null +++ b/src/atmosphericModels/kEpsilonLopesdaCosta/kEpsilonLopesdaCosta.H @@ -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 . + +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 kEpsilonLopesdaCosta +: + public eddyViscosity> +{ + // 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 kSource + ( + const volScalarField::Internal& magU, + const volScalarField::Internal& magU3 + ) const; + + virtual tmp 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 DkEff() const + { + return tmp + ( + new volScalarField + ( + "DkEff", + (this->nut_/sigmak_ + this->nu()) + ) + ); + } + + //- Return the effective diffusivity for epsilon + tmp DepsilonEff() const + { + return tmp + ( + new volScalarField + ( + "DepsilonEff", + (this->nut_/sigmaEps_ + this->nu()) + ) + ); + } + + //- Return the turbulence kinetic energy + virtual tmp k() const + { + return k_; + } + + //- Return the turbulence kinetic energy dissipation rate + virtual tmp 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 + +// ************************************************************************* // diff --git a/src/atmosphericModels/porosityModels/powerLawLopesdaCosta/powerLawLopesdaCosta.C b/src/atmosphericModels/porosityModels/powerLawLopesdaCosta/powerLawLopesdaCosta.C new file mode 100644 index 0000000000..e5d13f50e9 --- /dev/null +++ b/src/atmosphericModels/porosityModels/powerLawLopesdaCosta/powerLawLopesdaCosta.C @@ -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 . + +\*---------------------------------------------------------------------------*/ + +#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 groundPatches(coeffs.lookup("groundPatches")); + + // Functional form of the porosity surface area per unit volume as a + // function of the normalized vertical position + autoPtr> SigmaFunc + ( + Function1::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 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> groundSurfaceProcTris(Pstream::nProcs()); + List 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 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(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("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(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(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; +} + + +// ************************************************************************* // diff --git a/src/atmosphericModels/porosityModels/powerLawLopesdaCosta/powerLawLopesdaCosta.H b/src/atmosphericModels/porosityModels/powerLawLopesdaCosta/powerLawLopesdaCosta.H new file mode 100644 index 0000000000..f19677c808 --- /dev/null +++ b/src/atmosphericModels/porosityModels/powerLawLopesdaCosta/powerLawLopesdaCosta.H @@ -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 . + +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 + void apply + ( + scalarField& Udiag, + const scalarField& V, + const RhoFieldType& rho, + const vectorField& U + ) const; + + //- Apply resistance + template + 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 + +// ************************************************************************* // diff --git a/src/atmosphericModels/porosityModels/powerLawLopesdaCosta/powerLawLopesdaCostaTemplates.C b/src/atmosphericModels/porosityModels/powerLawLopesdaCosta/powerLawLopesdaCostaTemplates.C new file mode 100644 index 0000000000..be9142d191 --- /dev/null +++ b/src/atmosphericModels/porosityModels/powerLawLopesdaCosta/powerLawLopesdaCostaTemplates.C @@ -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 . + +\*---------------------------------------------------------------------------*/ + +#include "volFields.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +template +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 +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) + ); + } + } +} + + +// ************************************************************************* // diff --git a/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.H b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.H index bc7c496dd6..e3b00a535c 100644 --- a/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.H +++ b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.H @@ -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(); diff --git a/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModelI.H b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModelI.H index 63f28091b3..ecf384a274 100644 --- a/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModelI.H +++ b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModelI.H @@ -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_; diff --git a/src/fvOptions/sources/derived/explicitPorositySource/explicitPorositySource.C b/src/fvOptions/sources/derived/explicitPorositySource/explicitPorositySource.C index 8b8c37762c..a90b1903b3 100644 --- a/src/fvOptions/sources/derived/explicitPorositySource/explicitPorositySource.C +++ b/src/fvOptions/sources/derived/explicitPorositySource/explicitPorositySource.C @@ -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 diff --git a/src/fvOptions/sources/derived/explicitPorositySource/explicitPorositySource.H b/src/fvOptions/sources/derived/explicitPorositySource/explicitPorositySource.H index 576fac6843..c521774881 100644 --- a/src/fvOptions/sources/derived/explicitPorositySource/explicitPorositySource.H +++ b/src/fvOptions/sources/derived/explicitPorositySource/explicitPorositySource.H @@ -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 @@ -134,6 +134,11 @@ public: // Member Functions + const porosityModel& model() const + { + return porosityPtr_(); + } + // Add explicit and implicit contributions //- Add implicit contribution to momentum equation diff --git a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C index 1aaf5824f7..09bc91ed8d 100644 --- a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C +++ b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C @@ -2276,7 +2276,6 @@ Foam::triSurfaceTools::sideType Foam::triSurfaceTools::surfaceSide } -// triangulation of boundaryMesh Foam::triSurface Foam::triSurfaceTools::triangulate ( const polyBoundaryMesh& bMesh, @@ -2357,6 +2356,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 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 ( diff --git a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H index 84067be52a..a761186b90 100644 --- a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H +++ b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H @@ -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 | Copyright (C) 2016 OpenCFD Ltd. ------------------------------------------------------------------------------- License @@ -67,6 +67,7 @@ namespace Foam { // Forward declaration of classes +class boundBox; class edge; class labelledTri; class polyBoundaryMesh; @@ -76,7 +77,6 @@ class face; class Time; template class MeshedSurface; - /*---------------------------------------------------------------------------*\ Class triSurfaceTools Declaration \*---------------------------------------------------------------------------*/ @@ -484,6 +484,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