solarCalculator.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2015-2021 OpenCFD Ltd.
9-------------------------------------------------------------------------------
10License
11 This file is part of OpenFOAM.
12
13 OpenFOAM is free software: you can redistribute it and/or modify it
14 under the terms of the GNU General Public License as published by
15 the Free Software Foundation, either version 3 of the License, or
16 (at your option) any later version.
17
18 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 for more details.
22
23 You should have received a copy of the GNU General Public License
24 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25
26\*---------------------------------------------------------------------------*/
27
28#include "solarCalculator.H"
29#include "Time.H"
30#include "unitConversion.H"
31#include "constants.H"
32
33using namespace Foam::constant;
34
35// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
36
37namespace Foam
38{
40}
41
42
43const Foam::Enum
44<
46>
48({
49 { sunDirModel::mSunDirConstant, "constant" },
50 { sunDirModel::mSunDirTracking, "tracking" },
51
52 // old long names (v2012 and earlier)
53 { sunDirModel::mSunDirConstant, "sunDirConstant" },
54 { sunDirModel::mSunDirTracking, "sunDirTracking" }
55});
56
57
58const Foam::Enum
59<
61>
63({
64 { sunLModel::mSunLoadConstant, "constant" },
65 { sunLModel::mSunLoadTimeDependent, "timeDependent" },
66 { sunLModel::mSunLoadFairWeatherConditions, "fairWeather" },
67 { sunLModel::mSunLoadTheoreticalMaximum, "theoreticalMaximum" },
68
69 // old long names (v2012 and earlier)
70 { sunLModel::mSunLoadConstant, "sunLoadConstant" },
71 {
72 sunLModel::mSunLoadFairWeatherConditions,
73 "sunLoadFairWeatherConditions"
74 },
75 { sunLModel::mSunLoadTheoreticalMaximum, "sunLoadTheoreticalMaximum" }
76});
77
78
79// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
80
81void Foam::solarCalculator::calculateBetaTheta()
82{
83 scalar runTime = 0;
84
85 if (sunDirectionModel_ == mSunDirTracking)
86 {
87 runTime = mesh_.time().value();
88 }
89
90 const scalar LSM = 15.0*(dict_.get<scalar>("localStandardMeridian"));
91
92 const scalar D = dict_.get<scalar>("startDay") + runTime/86400.0;
93 const scalar M = 6.24004 + 0.0172*D;
94 const scalar EOT = -7.659*sin(M) + 9.863*sin(2*M + 3.5932);
95
96 dict_.readEntry("startTime", startTime_);
97
98 const scalar LST = startTime_ + runTime/3600.0;
99
100 const scalar LON = dict_.get<scalar>("longitude");
101
102 const scalar AST = LST + EOT/60.0 + (LON - LSM)/15;
103
104 const scalar delta = 23.45*sin(degToRad((360*(284 + D))/365));
105
106 const scalar H = degToRad(15*(AST - 12));
107
108 const scalar L = degToRad(dict_.get<scalar>("latitude"));
109
110 const scalar deltaRad = degToRad(delta);
111 beta_ = max(asin(cos(L)*cos(deltaRad)*cos(H) + sin(L)*sin(deltaRad)), 1e-3);
112 theta_ = acos((sin(beta_)*sin(L) - sin(deltaRad))/(cos(beta_)*cos(L)));
113
114 // theta is the angle between the SOUTH axis and the Sun
115 // If the hour angle is lower than zero (morning) the Sun is positioned
116 // on the East side.
117 if (H < 0)
118 {
119 theta_ += 2*(constant::mathematical::pi - theta_);
120 }
121
123 << tab << "altitude : " << radToDeg(beta_) << nl
124 << tab << "azimuth : " << radToDeg(theta_) << endl;
125}
126
127
128void Foam::solarCalculator::calculateSunDirection()
129{
130 gridUp_ = normalised(dict_.get<vector>("gridUp"));
131 eastDir_ = normalised(dict_.get<vector>("gridEast"));
132
133 coord_.reset
134 (
135 new coordinateSystem("grid", Zero, gridUp_, eastDir_)
136 );
137
138 // Assuming 'z' vertical, 'y' North and 'x' East
139 direction_.z() = -sin(beta_);
140 direction_.y() = cos(beta_)*cos(theta_); // South axis
141 direction_.x() = cos(beta_)*sin(theta_); // West axis
142
143 direction_.normalise();
144
146 << "Sun direction in absolute coordinates : " << direction_ <<endl;
147
148 // Transform to actual coordinate system
149 direction_ = coord_->transform(direction_);
150
152 << "Sun direction in the Grid coordinates : " << direction_ <<endl;
153}
154
155
156void Foam::solarCalculator::initialise()
157{
158 switch (sunDirectionModel_)
159 {
160 case mSunDirConstant:
161 {
162 if (dict_.readIfPresent("sunDirection", direction_))
163 {
164 direction_.normalise();
165 }
166 else
167 {
168 calculateBetaTheta();
169 calculateSunDirection();
170 }
171 break;
172 }
173 case mSunDirTracking:
174 {
175 if (word(mesh_.ddtScheme("default")) == "steadyState")
176 {
178 << " Sun direction model can not be sunDirtracking if the "
179 << " case is steady " << nl << exit(FatalError);
180 }
181
182 dict_.readEntry
183 (
184 "sunTrackingUpdateInterval",
185 sunTrackingUpdateInterval_
186 );
187
188 calculateBetaTheta();
189 calculateSunDirection();
190 break;
191 }
192 }
193
194 switch (sunLoadModel_)
195 {
196 case mSunLoadConstant:
197 {
198 dict_.readEntry("directSolarRad", directSolarRad_);
199 dict_.readEntry("diffuseSolarRad", diffuseSolarRad_);
200 break;
201 }
202 case mSunLoadTimeDependent:
203 {
204 directSolarRads_.reset
205 (
207 (
208 "directSolarRad",
209 dict_,
210 &mesh_
211 )
212 );
213
214 diffuseSolarRads_.reset
215 (
217 (
218 "diffuseSolarRad",
219 dict_,
220 &mesh_
221 )
222 );
223
224 directSolarRad_ =
225 directSolarRads_->value(mesh_.time().timeOutputValue());
226 diffuseSolarRad_ =
227 diffuseSolarRads_->value(mesh_.time().timeOutputValue());
228 break;
229 }
230 case mSunLoadFairWeatherConditions:
231 {
232 dict_.readIfPresent
233 (
234 "skyCloudCoverFraction",
235 skyCloudCoverFraction_
236 );
237
238 dict_.readEntry("A", A_);
239 dict_.readEntry("B", B_);
240 dict_.readEntry("C", C_);
241 dict_.readEntry("groundReflectivity", groundReflectivity_);
242 if (!dict_.readIfPresent("beta", beta_))
243 {
244 calculateBetaTheta();
245 }
246
247 directSolarRad_ =
248 (1.0 - 0.75*pow(skyCloudCoverFraction_, 3.0))
249 * A_/exp(B_/sin(beta_));
250 break;
251 }
252 case mSunLoadTheoreticalMaximum:
253 {
254 dict_.readEntry("Setrn", Setrn_);
255 dict_.readEntry("SunPrime", SunPrime_);
256 dict_.readEntry("groundReflectivity", groundReflectivity_);
257 dict_.readEntry("C", C_);
258
259 directSolarRad_ = Setrn_*SunPrime_;
260 break;
261 }
262 }
263}
264
265
266// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
267
269(
270 const dictionary& dict,
271 const fvMesh& mesh
272)
273:
274 mesh_(mesh),
275 dict_(dict),
276 sunDirectionModel_
277 (
278 sunDirectionModelTypeNames_.get("sunDirectionModel", dict)
279 ),
280 sunLoadModel_(sunLModelTypeNames_.get("sunLoadModel", dict)),
281 direction_(Zero),
282 sunTrackingUpdateInterval_(0),
283 startTime_(0),
284 gridUp_(Zero),
285 eastDir_(Zero),
286 coord_(),
287 directSolarRad_(0),
288 diffuseSolarRad_(0),
289 directSolarRads_(),
290 diffuseSolarRads_(),
291 skyCloudCoverFraction_(0),
292 groundReflectivity_(0),
293 A_(0),
294 B_(0),
295 beta_(0),
296 theta_(0),
297 C_(0.058),
298 Setrn_(0),
299 SunPrime_(0)
300{
301 initialise();
302}
303
304
305// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
306
308{
309 if (sunDirectionModel_ == mSunDirTracking)
310 {
311 calculateBetaTheta();
312 calculateSunDirection();
313 directSolarRad_ = A_/exp(B_/sin(max(beta_, ROOTVSMALL)));
314 }
315}
316
317
319{
320 if (sunLoadModel_ == mSunLoadTimeDependent)
321 {
322 directSolarRad_ = directSolarRads_->value(mesh_.time().value());
323 }
324}
325
326
328{
329 if (sunLoadModel_ == mSunLoadTimeDependent)
330 {
331 diffuseSolarRad_ = diffuseSolarRads_->value(mesh_.time().value());
332 }
333}
334
335
336// ************************************************************************* //
scalar delta
#define M(I)
Enum is a wrapper around a list of names/values that represent particular enumeration (or int) values...
Definition: Enum.H:61
Top level data entry class for use in dictionaries. Provides a mechanism to specify a variable as a c...
Definition: Function1.H:96
Base class for coordinate system specification, the default coordinate system type is cartesian .
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
bool readEntry(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX, bool mandatory=true) const
const Type & value() const
Return const reference to value.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:91
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:290
A solar calculator model providing models for the solar direction and solar loads.
void correctSunDirection()
Correct the Sun direction.
sunDirModel
Options for the Sun direction models.
static const Enum< sunLModel > sunLModelTypeNames_
Names for sunLModel.
void correctDiffuseSolarRad()
Correct diffuse solar irradiation.
void correctDirectSolarRad()
Correct direct solar irradiation.
static const Enum< sunDirModel > sunDirectionModelTypeNames_
Names for sunDirModel.
sunLModel
Options for the Sun load models.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
A class for handling words, derived from Foam::string.
Definition: word.H:68
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
dynamicFvMesh & mesh
engineTime & runTime
volScalarField H(IOobject("H", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE), mesh, dimensionedScalar(dimLength, Zero))
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
#define DebugInfo
Report an information message using Foam::Info.
constexpr scalar pi(M_PI)
Different types of constants.
Namespace for OpenFOAM.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
dimensionedScalar asin(const dimensionedScalar &ds)
dimensionedScalar exp(const dimensionedScalar &ds)
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
Definition: quaternionI.H:680
dimensionedScalar sin(const dimensionedScalar &ds)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
constexpr scalar degToRad() noexcept
Multiplication factor for degrees to radians conversion.
constexpr scalar radToDeg() noexcept
Multiplication factor for radians to degrees conversion.
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
error FatalError
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
dimensionedScalar cos(const dimensionedScalar &ds)
dimensionedScalar acos(const dimensionedScalar &ds)
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
constexpr char tab
The tab '\t' character(0x09)
Definition: Ostream.H:52
dictionary dict
const dimensionedScalar & D
volScalarField & e
Definition: createFields.H:11
Unit conversion functions.
const vector L(dict.get< vector >("L"))