Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
No results found
Show changes
Showing
with 617 additions and 370 deletions
......@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2017 OpenFOAM Foundation
Copyright (C) 2019 OpenCFD Ltd.
Copyright (C) 2019-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -92,7 +92,6 @@ Usage
spectralDistribution is the energy spectral distribution of the collimated
external beam.
SourceFiles
fvDOM.C
......@@ -121,8 +120,7 @@ class fvDOM
:
public radiationModel
{
// Private data
// Private Data
//- Incident radiation [W/m2]
volScalarField G_;
......@@ -184,6 +182,9 @@ class fvDOM
//- Spectral energy distribution for the external beam
scalarList spectralDistribution_;
//- Time-dependent spectral distributions
autoPtr<TimeFunction1<scalarField>> spectralDistributions_;
//- Solar calculator
autoPtr<solarCalculator> solarCalculator_;
......@@ -222,10 +223,10 @@ public:
//- Destructor
virtual ~fvDOM();
virtual ~fvDOM() = default;
// Member functions
// Member Functions
// Edit
......
......@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2015 OpenFOAM Foundation
Copyright (C) 2016-2020 OpenCFD Ltd.
Copyright (C) 2016-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -37,7 +37,6 @@ License
#include "wallPolyPatch.H"
#include "constants.H"
using namespace Foam::constant;
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
......@@ -83,7 +82,7 @@ void Foam::radiation::solarLoad::updateReflectedRays
{
if (includePatches[patchID])
{
for (label bandI = 0; bandI < nBands_; bandI++)
for (label bandI = 0; bandI < nBands_; ++bandI)
{
qrBf[patchID] +=
reflectedFaces_->qreflective(bandI).boundaryField()[patchID];
......@@ -94,7 +93,7 @@ void Foam::radiation::solarLoad::updateReflectedRays
const scalarField& sf = mesh_.magSf().boundaryField()[patchID];
const labelUList& cellIs = patches[patchID].faceCells();
for (label bandI = 0; bandI < nBands_; bandI++)
for (label bandI = 0; bandI < nBands_; ++bandI)
{
forAll(cellIs, i)
{
......@@ -109,9 +108,6 @@ void Foam::radiation::solarLoad::updateReflectedRays
}
}
}
}
......@@ -133,7 +129,7 @@ bool Foam::radiation::solarLoad::updateHitFaces()
}
case solarCalculator::mSunDirTracking:
{
label updateIndex = label
const label updateIndex = label
(
mesh_.time().value()/solarCalc_.sunTrackingUpdateInterval()
);
......@@ -167,7 +163,7 @@ void Foam::radiation::solarLoad::updateAbsorptivity
for (const label patchID : includePatches)
{
absorptivity_[patchID].setSize(nBands_);
for (label bandI = 0; bandI < nBands_; bandI++)
for (label bandI = 0; bandI < nBands_; ++bandI)
{
absorptivity_[patchID][bandI] =
boundaryRadiation.absorptivity(patchID, bandI);
......@@ -189,7 +185,7 @@ void Foam::radiation::solarLoad::updateDirectHitRadiation
// Reset qr and qrPrimary
qrBf = 0.0;
for (label bandI = 0; bandI < nBands_; bandI++)
for (label bandI = 0; bandI < nBands_; ++bandI)
{
volScalarField::Boundary& qprimaryBf =
qprimaryRad_[bandI].boundaryFieldRef();
......@@ -199,7 +195,7 @@ void Foam::radiation::solarLoad::updateDirectHitRadiation
forAll(hitFacesId, i)
{
const label faceI = hitFacesId[i];
label patchID = patches.whichPatch(faceI);
const label patchID = patches.whichPatch(faceI);
const polyPatch& pp = patches[patchID];
const label localFaceI = faceI - pp.start();
const vector qPrim =
......@@ -303,7 +299,7 @@ void Foam::radiation::solarLoad::updateSkyDiffusiveRadiation
const label cellI = cellIds[faceI];
if (includeMappedPatchBasePatches[patchID])
{
for (label bandI = 0; bandI < nBands_; bandI++)
for (label bandI = 0; bandI < nBands_; ++bandI)
{
qrBf[patchID][faceI] +=
(Ed + Er)
......@@ -313,7 +309,7 @@ void Foam::radiation::solarLoad::updateSkyDiffusiveRadiation
}
else
{
for (label bandI = 0; bandI < nBands_; bandI++)
for (label bandI = 0; bandI < nBands_; ++bandI)
{
Ru_[cellI] +=
(Ed + Er)
......@@ -328,6 +324,7 @@ void Foam::radiation::solarLoad::updateSkyDiffusiveRadiation
break;
case solarCalculator::mSunLoadConstant:
case solarCalculator::mSunLoadTimeDependent:
{
for (const label patchID : includePatches)
{
......@@ -340,7 +337,7 @@ void Foam::radiation::solarLoad::updateSkyDiffusiveRadiation
const label cellI = cellIds[faceI];
if (includeMappedPatchBasePatches[patchID])
{
for (label bandI = 0; bandI < nBands_; bandI++)
for (label bandI = 0; bandI < nBands_; ++bandI)
{
qrBf[patchID][faceI] +=
solarCalc_.diffuseSolarRad()
......@@ -350,7 +347,7 @@ void Foam::radiation::solarLoad::updateSkyDiffusiveRadiation
}
else
{
for (label bandI = 0; bandI < nBands_; bandI++)
for (label bandI = 0; bandI < nBands_; ++bandI)
{
Ru_[cellI] +=
(
......@@ -370,28 +367,24 @@ void Foam::radiation::solarLoad::updateSkyDiffusiveRadiation
void Foam::radiation::solarLoad::initialise(const dictionary& coeffs)
{
coeffs.readEntry("spectralDistribution", spectralDistribution_);
nBands_ = spectralDistribution_.size();
qprimaryRad_.setSize(nBands_);
spectralDistributions_.reset
(
new TimeFunction1<scalarField>
(
mesh_.time(), "spectralDistribution", coeffs
)
);
if (coeffs.readIfPresent("gridUp", verticalDir_))
{
verticalDir_.normalise();
}
else
{
const uniformDimensionedVectorField& g =
meshObjects::gravity::New(mesh_.time());
verticalDir_ = (-g/mag(g)).value();
}
spectralDistribution_ =
spectralDistributions_->value(mesh_.time().value());
coeffs.readEntry("useReflectedRays", useReflectedRays_);
nBands_ = spectralDistribution_.size();
spectralDistribution_ =
spectralDistribution_/sum(spectralDistribution_);
qprimaryRad_.setSize(nBands_);
forAll(qprimaryRad_, bandI)
{
qprimaryRad_.set
......@@ -413,9 +406,21 @@ void Foam::radiation::solarLoad::initialise(const dictionary& coeffs)
);
}
if (coeffs.readIfPresent("gridUp", verticalDir_))
{
verticalDir_.normalise();
}
else
{
const uniformDimensionedVectorField& g =
meshObjects::gravity::New(mesh_.time());
verticalDir_ = (-g/mag(g)).value();
}
coeffs.readIfPresent("solidCoupled", solidCoupled_);
coeffs.readIfPresent("wallCoupled", wallCoupled_);
coeffs.readIfPresent("updateAbsorptivity", updateAbsorptivity_);
coeffs.readEntry("useReflectedRays", useReflectedRays_);
}
/*
......@@ -721,6 +726,7 @@ void Foam::radiation::solarLoad::calculateQdiff
Foam::radiation::solarLoad::solarLoad(const volScalarField& T)
:
radiationModel(typeName, T),
solarCalc_(coeffs_, mesh_),
dict_(coeffs_),
qr_
(
......@@ -750,17 +756,18 @@ Foam::radiation::solarLoad::solarLoad(const volScalarField& T)
mesh_,
dimensionedScalar(dimMass/dimLength/pow3(dimTime), Zero)
),
solarCalc_(coeffs_, mesh_),
verticalDir_(Zero),
useReflectedRays_(false),
absorptivity_(mesh_.boundaryMesh().size()),
spectralDistribution_(),
nBands_(0),
spectralDistributions_(),
qprimaryRad_(0),
verticalDir_(Zero),
nBands_(0),
updateTimeIndex_(0),
solidCoupled_(true),
absorptivity_(mesh_.boundaryMesh().size()),
wallCoupled_(false),
updateAbsorptivity_(false),
firstIter_(true),
updateTimeIndex_(0)
useReflectedRays_(false),
firstIter_(true)
{
initialise(coeffs_);
}
......@@ -773,6 +780,7 @@ Foam::radiation::solarLoad::solarLoad
)
:
radiationModel(typeName, dict, T),
solarCalc_(dict, mesh_),
dict_(dict),
qr_
(
......@@ -802,18 +810,18 @@ Foam::radiation::solarLoad::solarLoad
mesh_,
dimensionedScalar(dimMass/dimLength/pow3(dimTime), Zero)
),
solarCalc_(dict, mesh_),
verticalDir_(Zero),
useReflectedRays_(false),
absorptivity_(mesh_.boundaryMesh().size()),
spectralDistribution_(),
nBands_(0),
spectralDistributions_(),
qprimaryRad_(0),
verticalDir_(Zero),
nBands_(0),
updateTimeIndex_(0),
solidCoupled_(true),
wallCoupled_(false),
absorptivity_(mesh_.boundaryMesh().size()),
updateAbsorptivity_(false),
firstIter_(true),
updateTimeIndex_(0)
useReflectedRays_(false),
firstIter_(true)
{
initialise(dict);
}
......@@ -821,12 +829,6 @@ Foam::radiation::solarLoad::solarLoad
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::label Foam::radiation::solarLoad::nBands() const
{
return nBands_;
}
bool Foam::radiation::solarLoad::read()
{
if (radiationModel::read())
......@@ -870,13 +872,25 @@ void Foam::radiation::solarLoad::calculate()
updateAbsorptivity(includePatches);
}
bool facesChanged = updateHitFaces();
const bool facesChanged = updateHitFaces();
if (facesChanged)
const bool timeDependentLoad =
solarCalc_.sunLoadModel() == solarCalculator::mSunLoadTimeDependent;
if (facesChanged || timeDependentLoad)
{
// Reset Ru
Ru_ = dimensionedScalar("Ru", dimMass/dimLength/pow3(dimTime), Zero);
solarCalc_.correctDirectSolarRad();
solarCalc_.correctDiffuseSolarRad();
spectralDistribution_ =
spectralDistributions_->value(mesh_.time().value());
spectralDistribution_ =
spectralDistribution_/sum(spectralDistribution_);
// Add direct hit radiation
const labelList& hitFacesId = hitFaces_->rayStartFaces();
updateDirectHitRadiation(hitFacesId, includeMappedPatchBasePatches);
......
......@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2015 OpenFOAM Foundation
Copyright (C) 2018 OpenCFD Ltd.
Copyright (C) 2018-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -31,7 +31,7 @@ Group
grpRadiationModels
Description
The solar load radiation model includes Sun primary hits, their
The \c solarLoad radiation model includes Sun primary hits, their
reflective fluxes and diffusive sky radiative fluxes.
The primary hit rays are calculated using a face shading algorithm.
......@@ -43,14 +43,52 @@ Description
Fair Weather Conditions Method from the ASHRAE Handbook.
By default the energy is included in cells adjacent to the patches into
the energy Equation (wallCoupled = false). On coupled patches the flux is
the energy equation (\c wallCoupled=false). On coupled patches the flux is
by default added to the wall and considered into the solid
(solidCoupled = true).
The solarLoad model can be used in conjuntion with fvDOM and viewFactor
radiation models. The flag useSolarLoad must be true on the rediation
dictionary.
(\c solidCoupled=true).
The \c solarLoad model can be used in conjuntion with \c fvDOM and
\c viewFactor radiation models. The flag \c useSolarLoad must be
\c true on the \c radiationProperties dictionary.
Usage
Minimal examples by using \c constant/radiationProperties:
\verbatim
solarLoadCoeffs
{
// Mandatory entries
useReflectedRays true;
spectralDistribution (1 5 1 2);
// Optional entries
solidCoupled true;
wallCoupled false;
updateAbsorptivity true;
// Mandatory/Optional (inherited) entries
...
}
\endverbatim
where the entries mean:
\table
Property | Description | Type | Reqd | Deflt
useReflectedRays | Flag to use reflected rays | bool | yes | -
spectralDistribution | Spectral distribution for the integrated <!--
--> solar heat flux | TimeFunction1\<scalarField\> | yes | -
solidCoupled | Flag to couple solids through mapped <!--
--> boundary patch using qr | bool | no | true
wallCoupled | Flag to couple wall patches using qr <!--
--> | bool | no | false
updateAbsorptivity | Flag to enable absorptivity updates <!--
--> | bool | no | false
\endtable
The inherited entries are elaborated in:
- \link radiationModel.H \endlink
- \link solarCalculator.H \endlink
- \link TimeFunction1.H \endlink
SourceFiles
solarLoad.C
......@@ -65,6 +103,7 @@ SourceFiles
#include "faceShading.H"
#include "faceReflecting.H"
#include "solarCalculator.H"
#include "TimeFunction1.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -81,10 +120,10 @@ class solarLoad
:
public radiationModel
{
// Private Data
private:
// Private data
//- Solar calculator
solarCalculator solarCalc_;
//- Dictionary
dictionary dict_;
......@@ -102,46 +141,46 @@ private:
//- and wallCoupled false
DimensionedField<scalar, volMesh> Ru_;
//- Solar calculator
solarCalculator solarCalc_;
//- Vertical direction (Default is g vector)
vector verticalDir_;
//- Include reflected rays from specular surfaces
bool useReflectedRays_;
//- Absorptivity list
List<List<tmp<scalarField>>> absorptivity_;
//- Spectral distribution for the integrated solar heat flux
scalarList spectralDistribution_;
//-Number of bands
label nBands_;
//- Time-dependent spectral distributions
autoPtr<TimeFunction1<scalarField>> spectralDistributions_;
//- Primary solar radiative heat flux per band [W/m2]
PtrList<volScalarField> qprimaryRad_;
//- Couple solids through mapped boundary patch using qr (default:true)
//- Vertical direction (default is g vector)
vector verticalDir_;
//- Number of bands
label nBands_;
//- Update Sun position index
label updateTimeIndex_;
//- Couple solids through mapped boundary patch using qr
bool solidCoupled_;
//- Couple wall patches using qr (default:false)
//- Couple wall patches using qr
bool wallCoupled_;
//- Absorptivity list
List<List<tmp<scalarField>>> absorptivity_;
//- Update absorptivity
bool updateAbsorptivity_;
//- Include reflected rays from specular surfaces
bool useReflectedRays_;
//- First iteration
bool firstIter_;
//- Update Sun position index
label updateTimeIndex_;
// Private Member Functions
//- Initialise
//- Initialise model parameters
void initialise(const dictionary&);
//- Update direct hit faces radiation
......@@ -166,6 +205,7 @@ private:
//- Update absorptivity
void updateAbsorptivity(const labelHashSet& includePatches);
//- No copy construct
solarLoad(const solarLoad&) = delete;
......@@ -192,15 +232,18 @@ public:
virtual ~solarLoad() = default;
// Member functions
// Member Functions
// Edit
// Evaluation
//- Solve
//- Read radiationProperties dictionary
bool read();
//- Solve radiation equations
void calculate();
//- Read radiation properties dictionary
bool read();
// Access
//- Source term component (for power of T^4)
virtual tmp<volScalarField> Rp() const;
......@@ -208,17 +251,17 @@ public:
//- Source term component (constant)
virtual tmp<DimensionedField<scalar, volMesh>> Ru() const;
// Access
//- Number of bands
label nBands() const;
//- Primary solar heat flux
const volScalarField& qprimaryRad(const label bandI) const
{
return qprimaryRad_[bandI];
}
//- Return const access to the number of bands
label nBands() const noexcept
{
return nBands_;
}
//- Return const access to the primary solar heat flux
const volScalarField& qprimaryRad(const label bandI) const
{
return qprimaryRad_[bandI];
}
};
......
......@@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2015-2020 OpenCFD Ltd.
Copyright (C) 2015-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -46,8 +46,12 @@ const Foam::Enum
>
Foam::solarCalculator::sunDirectionModelTypeNames_
({
{ sunDirModel::mSunDirConstant, "constant" },
{ sunDirModel::mSunDirTracking, "tracking" },
// old long names (v2012 and earlier)
{ sunDirModel::mSunDirConstant, "sunDirConstant" },
{ sunDirModel::mSunDirTracking, "sunDirTracking" },
{ sunDirModel::mSunDirTracking, "sunDirTracking" }
});
......@@ -55,14 +59,20 @@ const Foam::Enum
<
Foam::solarCalculator::sunLModel
>
Foam::solarCalculator::sunLoadModelTypeNames_
Foam::solarCalculator::sunLModelTypeNames_
({
{ sunLModel::mSunLoadConstant, "constant" },
{ sunLModel::mSunLoadTimeDependent, "timeDependent" },
{ sunLModel::mSunLoadFairWeatherConditions, "fairWeather" },
{ sunLModel::mSunLoadTheoreticalMaximum, "theoreticalMaximum" },
// old long names (v2012 and earlier)
{ sunLModel::mSunLoadConstant, "sunLoadConstant" },
{
sunLModel::mSunLoadFairWeatherConditions,
"sunLoadFairWeatherConditions"
},
{ sunLModel::mSunLoadTheoreticalMaximum, "sunLoadTheoreticalMaximum" },
{ sunLModel::mSunLoadTheoreticalMaximum, "sunLoadTheoreticalMaximum" }
});
......@@ -70,41 +80,34 @@ Foam::solarCalculator::sunLoadModelTypeNames_
void Foam::solarCalculator::calculateBetaTheta()
{
scalar runTime = 0.0;
switch (sunDirectionModel_)
scalar runTime = 0;
if (sunDirectionModel_ == mSunDirTracking)
{
case mSunDirTracking:
{
runTime = mesh_.time().value();
break;
}
case mSunDirConstant:
{
break;
}
runTime = mesh_.time().value();
}
scalar LSM = 15.0*(dict_.get<scalar>("localStandardMeridian"));
const scalar LSM = 15.0*(dict_.get<scalar>("localStandardMeridian"));
scalar D = dict_.get<scalar>("startDay") + runTime/86400.0;
scalar M = 6.24004 + 0.0172*D;
scalar EOT = -7.659*sin(M) + 9.863*sin(2*M + 3.5932);
const scalar D = dict_.get<scalar>("startDay") + runTime/86400.0;
const scalar M = 6.24004 + 0.0172*D;
const scalar EOT = -7.659*sin(M) + 9.863*sin(2*M + 3.5932);
dict_.readEntry("startTime", startTime_);
scalar LST = startTime_ + runTime/3600.0;
const scalar LST = startTime_ + runTime/3600.0;
scalar LON = dict_.get<scalar>("longitude");
const scalar LON = dict_.get<scalar>("longitude");
scalar AST = LST + EOT/60.0 + (LON - LSM)/15;
const scalar AST = LST + EOT/60.0 + (LON - LSM)/15;
scalar delta = 23.45*sin(degToRad((360*(284 + D))/365));
const scalar delta = 23.45*sin(degToRad((360*(284 + D))/365));
scalar H = degToRad(15*(AST - 12));
const scalar H = degToRad(15*(AST - 12));
scalar L = degToRad(dict_.get<scalar>("latitude"));
const scalar L = degToRad(dict_.get<scalar>("latitude"));
scalar deltaRad = degToRad(delta);
const scalar deltaRad = degToRad(delta);
beta_ = max(asin(cos(L)*cos(deltaRad)*cos(H) + sin(L)*sin(deltaRad)), 1e-3);
theta_ = acos((sin(beta_)*sin(L) - sin(deltaRad))/(cos(beta_)*cos(L)));
......@@ -150,7 +153,7 @@ void Foam::solarCalculator::calculateSunDirection()
}
void Foam::solarCalculator::init()
void Foam::solarCalculator::initialise()
{
switch (sunDirectionModel_)
{
......@@ -165,7 +168,6 @@ void Foam::solarCalculator::init()
calculateBetaTheta();
calculateSunDirection();
}
break;
}
case mSunDirTracking:
......@@ -197,6 +199,32 @@ void Foam::solarCalculator::init()
dict_.readEntry("diffuseSolarRad", diffuseSolarRad_);
break;
}
case mSunLoadTimeDependent:
{
directSolarRads_.reset
(
new TimeFunction1<scalar>
(
mesh_.time(),
"directSolarRad",
dict_
)
);
diffuseSolarRads_.reset
(
new TimeFunction1<scalar>
(
mesh_.time(),
"diffuseSolarRad",
dict_
)
);
directSolarRad_ = directSolarRads_->value(mesh_.time().value());
diffuseSolarRad_ = diffuseSolarRads_->value(mesh_.time().value());
break;
}
case mSunLoadFairWeatherConditions:
{
dict_.readIfPresent
......@@ -207,7 +235,8 @@ void Foam::solarCalculator::init()
dict_.readEntry("A", A_);
dict_.readEntry("B", B_);
dict_.readEntry("C", C_);
dict_.readEntry("groundReflectivity", groundReflectivity_);
if (!dict_.readIfPresent("beta", beta_))
{
calculateBetaTheta();
......@@ -216,17 +245,16 @@ void Foam::solarCalculator::init()
directSolarRad_ =
(1.0 - 0.75*pow(skyCloudCoverFraction_, 3.0))
* A_/exp(B_/sin(beta_));
dict_.readEntry("groundReflectivity", groundReflectivity_);
break;
}
case mSunLoadTheoreticalMaximum:
{
dict_.readEntry("Setrn", Setrn_);
dict_.readEntry("SunPrime", SunPrime_);
directSolarRad_ = Setrn_*SunPrime_;
dict_.readEntry("groundReflectivity", groundReflectivity_);
dict_.readEntry("C", C_);
directSolarRad_ = Setrn_*SunPrime_;
break;
}
}
......@@ -243,26 +271,32 @@ Foam::solarCalculator::solarCalculator
:
mesh_(mesh),
dict_(dict),
direction_(Zero),
directSolarRad_(0.0),
diffuseSolarRad_(0.0),
groundReflectivity_(0.0),
A_(0.0),
B_(0.0),
beta_(0.0),
theta_(0.0),
skyCloudCoverFraction_(0.0),
Setrn_(0.0),
SunPrime_(0.0),
C_(dict.get<scalar>("C")),
sunDirectionModel_
(
sunDirectionModelTypeNames_.get("sunDirectionModel", dict)
),
sunLoadModel_(sunLoadModelTypeNames_.get("sunLoadModel", dict)),
coord_()
sunLoadModel_(sunLModelTypeNames_.get("sunLoadModel", dict)),
direction_(Zero),
sunTrackingUpdateInterval_(0),
startTime_(0),
gridUp_(Zero),
eastDir_(Zero),
coord_(),
directSolarRad_(0),
diffuseSolarRad_(0),
directSolarRads_(),
diffuseSolarRads_(),
skyCloudCoverFraction_(0),
groundReflectivity_(0),
A_(0),
B_(0),
beta_(0),
theta_(0),
C_(0.058),
Setrn_(0),
SunPrime_(0)
{
init();
initialise();
}
......@@ -270,19 +304,29 @@ Foam::solarCalculator::solarCalculator
void Foam::solarCalculator::correctSunDirection()
{
switch (sunDirectionModel_)
if (sunDirectionModel_ == mSunDirTracking)
{
case mSunDirConstant:
{
break;
}
case mSunDirTracking:
{
calculateBetaTheta();
calculateSunDirection();
directSolarRad_ = A_/exp(B_/sin(max(beta_, ROOTVSMALL)));
break;
}
calculateBetaTheta();
calculateSunDirection();
directSolarRad_ = A_/exp(B_/sin(max(beta_, ROOTVSMALL)));
}
}
void Foam::solarCalculator::correctDirectSolarRad()
{
if (sunLoadModel_ == mSunLoadTimeDependent)
{
directSolarRad_ = directSolarRads_->value(mesh_.time().value());
}
}
void Foam::solarCalculator::correctDiffuseSolarRad()
{
if (sunLoadModel_ == mSunLoadTimeDependent)
{
diffuseSolarRad_ = diffuseSolarRads_->value(mesh_.time().value());
}
}
......
......@@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2015-2017 OpenCFD Ltd.
Copyright (C) 2015-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -27,58 +27,188 @@ Class
Foam::solarCalculator
Description
The solar calculator model provides information about the Sun direction
and Sun load model. The available models are:
For the Sun direction:
1) SunDirConstant : the direction is given in 'sunDirection'
2) SunDirTracking : the direction is calculated from the following
parameters:
localStandardMeridian : GMT (Local Zone Meridian) in hours
startDay : day from 1 to 365)
startTime: in hours
longitude: in degrees
latitude: in degrees
gridUp: grid orientation upwards
gridEast grid orientation eastwards
This model should be use in transient calculations.
The keyword 'sunTrackingUpdateInterval' (in hours) specifies on which
interval is the Sun direction updated.
Solar Load models available:
1) SunLoadConstant: direct and diffusive heat fluxes are provided by the
entries 'directSolarRad' and 'diffuseSolarRad'
2) SunLoadFairWeatherConditions: The solar fluxes are calculated following
the Fair Weather Conditions Method from the ASHRAE Handbook. The entries
are:
skyCloudCoverFraction: Fraction of sky covered by clouds (0-1)
A : Apparent solar irradiation at air mass m = 0
B : Atmospheric extinction coefficient
beta: Solar altitude (in degrees) above the horizontal. This
can be read or calculated providing the respective parameters
for Sun position explained above.
groundReflectivity : ground reflectivity
In this model the flux is calculated as:
directSolarRad =
(1 - 0.75*skyCloudCoverFraction^3)*A/exp(B/sin(beta));
3) SunLoadTheoreticalMaximum: The entries are:
Setrn
SunPrime:
groundReflectivity : ground reflectivity
In this model the flux is calculated as:
directSolarRad = Setrn*SunPrime;
The diffuse on vertical/horizontal walls and ground-reflected radiation are
A solar calculator model providing models
for the solar direction and solar loads.
Available models for the solar direction:
- \c constant: Constant sunbeam direction.
- \c tracking: Transient model calculating sunbeam direction
based on a given set of parameters.
Available models for the solar load:
- \c constant: Constant solar load.
- \c timeDependent: Time-dependent solar load.
- \c fairWeather: Solar fluxes are calculated following
the "Fair Weather Conditions Method from the ASHRAE Handbook".
- \c theoreticalMaximum: Theoretically maximum solar load.
Usage
Minimal examples by using \c constant/radiationProperties:
\c sunDirectionModel - Option-1:
\verbatim
solarLoadCoeffs
{
sunDirectionModel constant;
sunDirection (1 0 -1);
}
\endverbatim
where the entries mean:
\table
Property | Description | Type | Reqd | Deflt
sunDirection | Sunbeam direction | vector | no | calculated
\endtable
\c sunDirectionModel - Option-2:
\verbatim
solarLoadCoeffs
{
sunDirectionModel tracking;
sunTrackingUpdateInterval 800;
localStandardMeridian 9;
startDay 204;
startTime 15;
longitude 139.74;
latitude 35.658;
gridUp (0 0 1);
gridEast (1 0 0);
}
\endverbatim
where the entries mean:
\table
Property | Description | Type | Reqd | Deflt
sunTrackingUpdateInterval | Interval to update the Sun direction <!--
--> [decimal hours] | scalar | yes | -
localStandardMeridian | GMT (Local Zone Meridian) [decimal hours]<!--
--> | scalar | yes | -
startDay | Day from 1 to 365 | scalar | yes | -
startTime | Start time for the Sun position [decimal hours] <!--
--> | scalar | yes | -
longitude | Geographic coordinate specifying the east–west <!--
--> position of a point on the surface of a planetary <!--
--> body [degree] | scalar | yes | -
latitude | Geographic coordinate specifying the north–south <!--
--> position of a point on the surface of a planetary <!--
--> body [degree] | scalar | yes | -
gridUp | Grid orientation upwards | vector | yes | -
gridEast | Grid orientation eastwards | vector | yes | -
\endtable
\c sunLoadModel - Option-1:
\verbatim
solarLoadCoeffs
{
sunLoadModel constant;
directSolarRad 100;
diffuseSolarRad 0;
}
\endverbatim
where the entries mean:
\table
Property | Description | Type | Reqd | Deflt
directSolarRad | Direct solar irradiation [W/m2] | scalar | yes | -
diffuseSolarRad | Diffuse solar irradiation on vertical surfaces <!--
--> [W/m2] | scalar | yes | -
\endtable
\c sunLoadModel - Option-2:
\verbatim
solarLoadCoeffs
{
sunLoadModel timeDependent;
directSolarRad <TimeFunction1<scalar>>;
diffuseSolarRad <TimeFunction1<scalar>>;
}
\endverbatim
where the entries mean:
\table
Property | Description | Type | Reqd | Deflt
directSolarRad | Time-series of direct solar irradiation <!--
--> [W/m2] | TimeFunction1\<scalar\> | yes | -
diffuseSolarRad | Time-series of diffuse solar irradiation on <!--
--> vertical surfaces [W/m2] <!--
--> | TimeFunction1\<scalar\> | yes | -
\endtable
\c sunLoadModel - Option-3:
\verbatim
solarLoadCoeffs
{
sunLoadModel fairWeather;
skyCloudCoverFraction 0.25;
groundReflectivity 1.0;
A 0.1;
B 0.2;
C 0.058;
beta 0.15;
}
\endverbatim
where the entries mean:
\table
Property | Description | Type | Reqd | Deflt
A | Apparent solar irradiation at air mass m = 0 <!--
--> | scalar | yes | -
B | Atmospheric extinction coefficient <!--
--> | scalar | yes | -
C | Solar diffusivity constant | scalar | yes | -
groundReflectivity | Ground reflectivity | scalar | yes | -
skyCloudCoverFraction | Fraction of sky covered by clouds [0,1] <!--
--> | scalar | no | 0
beta | Solar altitude (in degrees) above the horizontal <!--
--> | scalar | no | calculated
\endtable
In this model the flux is calculated as:
\verbatim
directSolarRad = (1 - 0.75*skyCloudCoverFraction^3)*A/exp(B/sin(beta));
\endverbatim
\c sunLoadModel - Option-4:
\verbatim
solarLoadCoeffs
{
sunLoadModel theoreticalMaximum;
Setrn 1.0;
SunPrime 4.0;
groundReflectivity 1.0;
C 0.058;
}
\endverbatim
where the entries mean:
\table
Property | Description | Type | Reqd | Deflt
Setrn | Parameter in maximum theoretical direct solar <!--
--> model | scalar | yes | -
SunPrime | Parameter in maximum theoretical direct solar <!--
--> model | scalar | yes | -
groundReflectivity | Ground reflectivity | scalar | yes | -
C | Solar diffusivity constant | scalar | yes | -
\endtable
In this model the flux is calculated as:
\verbatim
directSolarRad = Setrn*SunPrime;
\endverbatim
Note
- The \c sunDirectionModel:tracking can only be used
in transient calculations.
- The keyword \c sunTrackingUpdateInterval (in hours) specifies on which
interval is the Sun direction updated.
- The diffuse on vertical/horizontal walls and ground-reflected radiation are
calculated following the ASHRAE Handbook.
- The range of \c skyCloudCoverFraction is [0,1].
SourceFiles
solarCalculator.C
......@@ -93,6 +223,7 @@ SourceFiles
#include "DynamicField.H"
#include "HashSet.H"
#include "coordinateSystem.H"
#include "TimeFunction1.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -107,37 +238,34 @@ class solarCalculator
{
public:
// Public enumeration
// Public Enumeration
//- Sun direction models
//- Options for the Sun direction models
enum sunDirModel
{
mSunDirConstant,
mSunDirTracking
};
//- Direct sun load models
//- Names for sunDirModel
static const Enum<sunDirModel> sunDirectionModelTypeNames_;
//- Options for the Sun load models
enum sunLModel
{
mSunLoadConstant,
mSunLoadTimeDependent,
mSunLoadFairWeatherConditions,
mSunLoadTheoreticalMaximum
};
protected:
//- Sun direction models
static const Enum<sunDirModel> sunDirectionModelTypeNames_;
//- Sun load models
static const Enum<sunLModel> sunLoadModelTypeNames_;
//- Names for sunLModel
static const Enum<sunLModel> sunLModelTypeNames_;
private:
// Private data
// Private Data
//- Reference to mesh
const fvMesh& mesh_;
......@@ -145,56 +273,78 @@ private:
//- Dictionary
dictionary dict_;
//- Direction
vector direction_;
//- Sun direction model
sunDirModel sunDirectionModel_;
//- Direct solar irradiation
scalar directSolarRad_;
//- Sun load model
sunLModel sunLoadModel_;
//- Diffuse solar irradiation on vertical surfaces
scalar diffuseSolarRad_;
//- Ground reflectivity
scalar groundReflectivity_;
// sunDirectionModel = constant
//- Fair weather direct solar load model parameters
scalar A_;
scalar B_;
scalar beta_;
scalar theta_;
//- Sunbeam direction
vector direction_;
//- Sky cloud cover fraction [0-1]
scalar skyCloudCoverFraction_;
// sunDirectionModel = tracking
//- Maximum theoretical direct solar load model parameters
scalar Setrn_;
scalar SunPrime_;
//- Interval to update the Sun direction [decimal hours]
scalar sunTrackingUpdateInterval_;
//- Start time for the Sun position [decimal hours]
scalar startTime_;
//- Diffusive solar load model parameters
scalar C_;
//- Up grid orientation
vector gridUp_;
//- Sun direction model
sunDirModel sunDirectionModel_;
//- East grid orientation
vector eastDir_;
//- Sun load model
sunLModel sunLoadModel_;
//- Grid coordinate system
autoPtr<coordinateSystem> coord_;
// sunLoadModel = constant
//- Direct solar irradiation
scalar directSolarRad_;
//- Diffuse solar irradiation on vertical surfaces
scalar diffuseSolarRad_;
// sunLoadModel = timeDependent
//- Time-series of direct solar irradiation
autoPtr<TimeFunction1<scalar>> directSolarRads_;
//- Time-series of diffuse solar irradiation on vertical surfaces
autoPtr<TimeFunction1<scalar>> diffuseSolarRads_;
// sunLoadModel = fairWeather
//- Sky cloud cover fraction [0-1]
scalar skyCloudCoverFraction_;
//- Grid coordinate system
autoPtr<coordinateSystem> coord_;
//- Ground reflectivity
scalar groundReflectivity_;
//- East grid orientation
vector eastDir_;
//- Fair weather direct solar load model parameters
scalar A_;
scalar B_;
scalar beta_;
scalar theta_;
//- Up grid orientation
vector gridUp_;
//- Diffusive solar load model parameter
scalar C_;
//- Interval in decimal hours to update Sun direction for SunDirTracking
scalar sunTrackingUpdateInterval_;
//- Start time for the Sun position (decimal hours)
scalar startTime_;
// sunLoadModel = theoreticalMaximum
//- Maximum theoretical direct solar load model parameters
scalar Setrn_;
scalar SunPrime_;
//- No copy construct
......@@ -204,15 +354,15 @@ private:
void operator=(const solarCalculator&) = delete;
// Private members
// Private Member Functions
//- Init
void init();
//- Initialise model parameters
void initialise();
//- Calculate beta and theta angles
void calculateBetaTheta();
//- Calculate Sun direction
//- Calculate the Sun direction
void calculateSunDirection();
......@@ -234,101 +384,109 @@ public:
// Member Functions
// Evaluation
//- Correct the Sun direction
void correctSunDirection();
//- Correct direct solar irradiation
void correctDirectSolarRad();
//- Correct diffuse solar irradiation
void correctDiffuseSolarRad();
// Access
//- const access to direction
const vector& direction() const
//- Return const access to the Sun direction model
const sunDirModel& sunDirectionModel() const noexcept
{
return direction_;
return sunDirectionModel_;
}
//- Return const access to the Sun load model
const sunLModel& sunLoadModel() const noexcept
{
return sunLoadModel_;
}
//- Non-const access to direction
//- Return non-const access to the Sun direction
vector& direction()
{
return direction_;
}
//- Return direct solar irradiation
//- Return const access to the Sun direction
const vector& direction() const noexcept
{
return direction_;
}
//- Return non-const access to the direct solar irradiation
scalar& directSolarRad()
{
return directSolarRad_;
}
//- Return const access to direct solar irradiation
const scalar& directSolarRad() const
//- Return const access to the direct solar irradiation
const scalar& directSolarRad() const noexcept
{
return directSolarRad_;
}
//- Return diffuse solar irradiation
//- Return non-const access to the diffuse solar irradiation
scalar& diffuseSolarRad()
{
return diffuseSolarRad_;
}
//- Return diffuse solar irradiation
const scalar& diffuseSolarRad() const
//- Return const access to the diffuse solar irradiation
const scalar& diffuseSolarRad() const noexcept
{
return diffuseSolarRad_;
}
//- Return C constant
scalar C()
//- Return const access to the C constant
scalar C() const noexcept
{
return C_;
}
//- Return beta
scalar beta()
//- Return const access to beta
scalar beta() const noexcept
{
return beta_;
}
//- Return theta
scalar theta()
//- Return const access to theta
scalar theta() const noexcept
{
return theta_;
}
//- Return Sun direction model
sunDirModel sunDirectionModel() const
{
return sunDirectionModel_;
}
//- Return Sun load model
sunLModel sunLoadModel() const
{
return sunLoadModel_;
}
//- Return groundReflectivity
scalar groundReflectivity()
//- Return const access to the ground reflectivity
scalar groundReflectivity() const noexcept
{
return groundReflectivity_;
}
//- Return coordinateSystem
const coordinateSystem& coord()
//- Return const access to the coordinate system
const coordinateSystem& coord() const noexcept
{
return *coord_;
}
//- Return sunTrackingUpdateInterval
scalar sunTrackingUpdateInterval()
//- Return const access to sunTrackingUpdateInterval
scalar sunTrackingUpdateInterval() const noexcept
{
return sunTrackingUpdateInterval_;
}
//- Return startTime
scalar startTime()
//- Return const access to startTime
scalar startTime() const noexcept
{
return startTime_;
}
//- Recalculate
void correctSunDirection();
};
......
......@@ -249,8 +249,13 @@ void Foam::isoAdvection::boundFlux
// First try to pass surplus fluid on to neighbour cells that are
// not filled and to which dVf < phi*dt
while (mag(alphaOvershoot) > aTol && nFacesToPassFluidThrough > 0)
for (label iter=0; iter<10; iter++)
{
if(mag(alphaOvershoot) < aTol || nFacesToPassFluidThrough == 0)
{
break;
}
DebugInfo
<< "\n\nBounding cell " << celli
<< " with alpha overshooting " << alphaOvershoot
......
......@@ -233,8 +233,7 @@ void Foam::reconstruction::plicRDF::setInitNormals(bool interpolate)
void Foam::reconstruction::plicRDF::calcResidual
(
Map<scalar>& normalResidual,
Map<scalar>& avgAngle
List<normalRes>& normalResidual
)
{
exchangeFields_.setUpCommforZone(interfaceCell_,false);
......@@ -246,7 +245,6 @@ void Foam::reconstruction::plicRDF::calcResidual
const labelListList& stencil = exchangeFields_.getStencil();
normalResidual.clear();
forAll(interfaceLabels_, i)
{
......@@ -260,13 +258,13 @@ void Foam::reconstruction::plicRDF::calcResidual
scalar maxDiffNormal = GREAT;
scalar weight= 0;
const vector cellNormal = normal_[celli]/mag(normal_[celli]);
for (const label gblIdx : stencil[celli])
forAll(stencil[celli],j)
{
const label gblIdx = stencil[celli][j];
vector normal =
exchangeFields_.getValue(normal_, mapNormal, gblIdx);
if (mag(normal) != 0 && i != 0)
if (mag(normal) != 0 && j != 0)
{
vector n = normal/mag(normal);
scalar cosAngle = max(min((cellNormal & n), 1), -1);
......@@ -291,8 +289,9 @@ void Foam::reconstruction::plicRDF::calcResidual
vector newCellNormal = normalised(interfaceNormal_[i]);
scalar normalRes = (1 - (cellNormal & newCellNormal));
avgAngle.insert(celli, avgDiffNormal);
normalResidual.insert(celli, normalRes);
normalResidual[i].celli = celli;
normalResidual[i].normalResidual = normalRes;
normalResidual[i].avgAngle = avgDiffNormal;
}
}
......@@ -399,14 +398,14 @@ void Foam::reconstruction::plicRDF::reconstruct(bool forceUpdate)
// nextToInterface is update on setInitNormals
const boolList& nextToInterface_ = RDF_.nextToInterface();
labelHashSet tooCoarse;
bitSet tooCoarse(mesh_.nCells(),false);
for (int iter=0; iter<iteration_; ++iter)
{
forAll(interfaceLabels_, i)
{
const label celli = interfaceLabels_[i];
if (mag(interfaceNormal_[i]) == 0 || tooCoarse.found(celli))
if (mag(interfaceNormal_[i]) == 0 || tooCoarse.test(celli))
{
continue;
}
......@@ -439,8 +438,7 @@ void Foam::reconstruction::plicRDF::reconstruct(bool forceUpdate)
normal_.correctBoundaryConditions();
centre_.correctBoundaryConditions();
Map<scalar> residual;
Map<scalar> avgAngle;
List<normalRes> normalResidual(interfaceLabels_.size());
surfaceVectorField::Boundary nHatb(mesh_.Sf().boundaryField());
nHatb *= 1/(mesh_.magSf().boundaryField());
......@@ -456,43 +454,42 @@ void Foam::reconstruction::plicRDF::reconstruct(bool forceUpdate)
);
RDF_.updateContactAngle(alpha1_, U_, nHatb);
gradSurf(RDF_);
calcResidual(residual, avgAngle);
calcResidual(normalResidual);
}
label resCounter = 0;
scalar avgRes = 0;
scalar avgNormRes = 0;
Map<scalar>::iterator resIter = residual.begin();
Map<scalar>::iterator avgAngleIter = avgAngle.begin();
while (resIter.found())
forAll(normalResidual,i)
{
if (avgAngleIter() > 0.26 && iter > 0) // 15 deg
const label celli = normalResidual[i].celli;
const scalar normalRes= normalResidual[i].normalResidual;
const scalar avgA = normalResidual[i].avgAngle;
if (avgA > 0.26 && iter > 0) // 15 deg
{
tooCoarse.set(resIter.key());
tooCoarse.set(celli);
}
else
{
avgRes += resIter();
avgRes += normalRes;
scalar normRes = 0;
scalar discreteError = 0.01*sqr(avgAngleIter());
scalar discreteError = 0.01*sqr(avgA);
if (discreteError != 0)
{
normRes= resIter()/max(discreteError, tol_);
normRes= normalRes/max(discreteError, tol_);
}
else
{
normRes= resIter()/tol_;
normRes= normalRes/tol_;
}
avgNormRes += normRes;
resCounter++;
}
++resIter;
++avgAngleIter;
}
reduce(avgRes,sumOp<scalar>());
......
......@@ -78,6 +78,14 @@ class plicRDF
{
// Private Data
//- residuals storage
struct normalRes
{
label celli;
scalar normalResidual;
scalar avgAngle;
};
//- Reference to mesh
const fvMesh& mesh_;
......@@ -129,8 +137,7 @@ class plicRDF
//- compute the normal residuals
void calcResidual
(
Map<scalar>& normalResidual,
Map<scalar>& avgAngle
List<normalRes>& normalResidual
);
//- interpolation of the normals from the previous time step
......
......@@ -10,7 +10,6 @@ FoamFile
version 2.0;
format ascii;
class volVectorField;
location "0";
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......
......@@ -10,7 +10,6 @@ FoamFile
version 2.0;
format ascii;
class volScalarField;
location "1";
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -21,27 +20,7 @@ internalField uniform 0;
boundaryField
{
patch0_half0
{
type cyclic;
}
patch1_half0
{
type cyclic;
}
patch2_half0
{
type cyclic;
}
patch2_half1
{
type cyclic;
}
patch1_half1
{
type cyclic;
}
patch0_half1
".*"
{
type cyclic;
}
......
......@@ -3,7 +3,6 @@ cd "${0%/*}" || exit # Run from this directory
. ${WM_PROJECT_DIR:?}/bin/tools/CleanFunctions # Tutorial clean functions
#------------------------------------------------------------------------------
cleanCase
rm -f 0/enstrophy
cleanCase0
#------------------------------------------------------------------------------
......@@ -3,9 +3,14 @@ cd "${0%/*}" || exit # Run from this directory
. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions # Tutorial run functions
#------------------------------------------------------------------------------
restore0Dir
runApplication blockMesh
runApplication boxTurb
runApplication $(getApplication)
runApplication -s enstrophy postProcess -func enstrophy
runApplication postProcess -func enstrophy
#------------------------------------------------------------------------------
......@@ -10,7 +10,6 @@ FoamFile
version 2.0;
format ascii;
class dictionary;
location "constant";
object boxTurbDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......
......@@ -10,7 +10,6 @@ FoamFile
version 2.0;
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......
......@@ -10,7 +10,6 @@ FoamFile
version 2.0;
format ascii;
class dictionary;
location "constant";
object turbulenceProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......
......@@ -10,7 +10,6 @@ FoamFile
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......
......@@ -10,7 +10,6 @@ FoamFile
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......
......@@ -10,7 +10,6 @@ FoamFile
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......
......@@ -53,4 +53,5 @@ boundaryField
}
}
// ************************************************************************* //
......@@ -52,4 +52,5 @@ boundaryField
}
}
// ************************************************************************* //