Line data Source code
1 : // EnergyPlus, Copyright (c) 1996-2025, The Board of Trustees of the University of Illinois,
2 : // The Regents of the University of California, through Lawrence Berkeley National Laboratory
3 : // (subject to receipt of any required approvals from the U.S. Dept. of Energy), Oak Ridge
4 : // National Laboratory, managed by UT-Battelle, Alliance for Sustainable Energy, LLC, and other
5 : // contributors. All rights reserved.
6 : //
7 : // NOTICE: This Software was developed under funding from the U.S. Department of Energy and the
8 : // U.S. Government consequently retains certain rights. As such, the U.S. Government has been
9 : // granted for itself and others acting on its behalf a paid-up, nonexclusive, irrevocable,
10 : // worldwide license in the Software to reproduce, distribute copies to the public, prepare
11 : // derivative works, and perform publicly and display publicly, and to permit others to do so.
12 : //
13 : // Redistribution and use in source and binary forms, with or without modification, are permitted
14 : // provided that the following conditions are met:
15 : //
16 : // (1) Redistributions of source code must retain the above copyright notice, this list of
17 : // conditions and the following disclaimer.
18 : //
19 : // (2) Redistributions in binary form must reproduce the above copyright notice, this list of
20 : // conditions and the following disclaimer in the documentation and/or other materials
21 : // provided with the distribution.
22 : //
23 : // (3) Neither the name of the University of California, Lawrence Berkeley National Laboratory,
24 : // the University of Illinois, U.S. Dept. of Energy nor the names of its contributors may be
25 : // used to endorse or promote products derived from this software without specific prior
26 : // written permission.
27 : //
28 : // (4) Use of EnergyPlus(TM) Name. If Licensee (i) distributes the software in stand-alone form
29 : // without changes from the version obtained under this License, or (ii) Licensee makes a
30 : // reference solely to the software portion of its product, Licensee must refer to the
31 : // software as "EnergyPlus version X" software, where "X" is the version number Licensee
32 : // obtained under this License and may not use a different name for the software. Except as
33 : // specifically required in this Section (4), Licensee shall not use in a company name, a
34 : // product name, in advertising, publicity, or other promotional activities any name, trade
35 : // name, trademark, logo, or other designation of "EnergyPlus", "E+", "e+" or confusingly
36 : // similar designation, without the U.S. Department of Energy's prior written consent.
37 : //
38 : // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
39 : // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY
40 : // AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
41 : // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
42 : // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
43 : // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
44 : // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
45 : // OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
46 : // POSSIBILITY OF SUCH DAMAGE.
47 :
48 : // C++ Headers
49 : #include <algorithm>
50 :
51 : // ObjexxFCL Headers
52 : #include <ObjexxFCL/Fmath.hh>
53 :
54 : // EnergyPlus Headers
55 : #include <EnergyPlus/Data/EnergyPlusData.hh>
56 : #include <EnergyPlus/DataEnvironment.hh>
57 : #include <EnergyPlus/DataGlobals.hh>
58 : #include <EnergyPlus/DataIPShortCuts.hh>
59 : #include <EnergyPlus/DataReportingFlags.hh>
60 : #include <EnergyPlus/General.hh>
61 : #include <EnergyPlus/GroundTemperatureModeling/FiniteDifferenceGroundTemperatureModel.hh>
62 : #include <EnergyPlus/GroundTemperatureModeling/KusudaAchenbachGroundTemperatureModel.hh>
63 : #include <EnergyPlus/InputProcessing/InputProcessor.hh>
64 : #include <EnergyPlus/UtilityRoutines.hh>
65 : #include <EnergyPlus/WeatherManager.hh>
66 :
67 : namespace EnergyPlus {
68 :
69 : namespace GroundTemp {
70 : //******************************************************************************
71 :
72 : // Finite difference model factory
73 1 : FiniteDiffGroundTempsModel *FiniteDiffGroundTempsModel::FiniteDiffGTMFactory(EnergyPlusData &state, const std::string &objectName)
74 : {
75 : // SUBROUTINE INFORMATION:
76 : // AUTHOR Matt Mitchell
77 : // DATE WRITTEN Summer 2015
78 :
79 : // PURPOSE OF THIS SUBROUTINE:
80 : // Read input and creates instance of finite difference ground temp model
81 :
82 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
83 1 : bool found = false;
84 : int NumNums;
85 : int NumAlphas;
86 : int IOStat;
87 :
88 : // New shared pointer for this model object
89 1 : auto *thisModel = new FiniteDiffGroundTempsModel();
90 :
91 1 : GroundTemp::ModelType modelType = GroundTemp::ModelType::FiniteDiff;
92 :
93 : // Search through finite diff models here
94 1 : std::string_view const cCurrentModuleObject = GroundTemp::modelTypeNamesUC[(int)modelType];
95 1 : const int numCurrModels = state.dataInputProcessing->inputProcessor->getNumObjectsFound(state, cCurrentModuleObject);
96 :
97 1 : for (int modelNum = 1; modelNum <= numCurrModels; ++modelNum) {
98 :
99 2 : state.dataInputProcessing->inputProcessor->getObjectItem(state,
100 : cCurrentModuleObject,
101 : modelNum,
102 1 : state.dataIPShortCut->cAlphaArgs,
103 : NumAlphas,
104 1 : state.dataIPShortCut->rNumericArgs,
105 : NumNums,
106 : IOStat);
107 :
108 1 : if (objectName == state.dataIPShortCut->cAlphaArgs(1)) {
109 : // Read input into object here
110 :
111 1 : thisModel->modelType = modelType;
112 1 : thisModel->Name = state.dataIPShortCut->cAlphaArgs(1);
113 1 : thisModel->baseConductivity = state.dataIPShortCut->rNumericArgs(1);
114 1 : thisModel->baseDensity = state.dataIPShortCut->rNumericArgs(2);
115 1 : thisModel->baseSpecificHeat = state.dataIPShortCut->rNumericArgs(3);
116 1 : thisModel->waterContent = state.dataIPShortCut->rNumericArgs(4) / 100.0;
117 1 : thisModel->saturatedWaterContent = state.dataIPShortCut->rNumericArgs(5) / 100.0;
118 1 : thisModel->evapotransCoeff = state.dataIPShortCut->rNumericArgs(6);
119 :
120 1 : found = true;
121 1 : break;
122 : }
123 : }
124 :
125 1 : if (found) {
126 1 : state.dataGrndTempModelMgr->groundTempModels.push_back(thisModel);
127 :
128 : // Simulate
129 1 : thisModel->initAndSim(state);
130 :
131 : // Return the pointer
132 1 : return thisModel;
133 : }
134 :
135 0 : ShowFatalError(state, fmt::format("{}--Errors getting input for ground temperature model", GroundTemp::modelTypeNames[(int)modelType]));
136 0 : return nullptr;
137 : }
138 :
139 : //******************************************************************************
140 :
141 1 : void FiniteDiffGroundTempsModel::initAndSim(EnergyPlusData &state)
142 : {
143 : // SUBROUTINE INFORMATION:
144 : // AUTHOR Matt Mitchell
145 : // DATE WRITTEN Summer 2015
146 :
147 : // PURPOSE OF THIS SUBROUTINE:
148 : // Initializes and simulates finite difference ground temps model
149 :
150 1 : getWeatherData(state);
151 :
152 1 : developMesh();
153 :
154 1 : performSimulation(state);
155 1 : }
156 :
157 : //******************************************************************************
158 :
159 3 : void FiniteDiffGroundTempsModel::getWeatherData(EnergyPlusData &state)
160 : {
161 : // SUBROUTINE INFORMATION:
162 : // AUTHOR Matt Mitchell
163 : // DATE WRITTEN Summer 2015
164 :
165 : // PURPOSE OF THIS SUBROUTINE:
166 : // Finds correct environment for reading all weather data. Loops over all weather data in weather file
167 : // and data structure containing daily average of required weather data.
168 :
169 : // SUBROUTINE ARGUMENT DEFINITIONS:
170 : // an environment is available to process
171 :
172 : // Save current environment so we can revert back when done
173 3 : const int Envrn_reset = state.dataWeather->Envrn;
174 3 : const Constant::KindOfSim KindOfSim_reset = state.dataGlobal->KindOfSim;
175 3 : const int TimeStep_reset = state.dataGlobal->TimeStep;
176 3 : const int HourOfDay_reset = state.dataGlobal->HourOfDay;
177 3 : const bool BeginEnvrnFlag_reset = state.dataGlobal->BeginEnvrnFlag;
178 3 : const bool EndEnvrnFlag_reset = state.dataGlobal->EndEnvrnFlag;
179 3 : const bool EndMonthFlag_reset = state.dataEnvrn->EndMonthFlag;
180 3 : const bool WarmupFlag_reset = state.dataGlobal->WarmupFlag;
181 3 : const int DayOfSim_reset = state.dataGlobal->DayOfSim;
182 3 : const std::string DayOfSimChr_reset = state.dataGlobal->DayOfSimChr;
183 3 : const int NumOfWarmupDays_reset = state.dataReportFlag->NumOfWarmupDays;
184 3 : const bool BeginDayFlag_reset = state.dataGlobal->BeginDayFlag;
185 3 : const bool EndDayFlag_reset = state.dataGlobal->EndDayFlag;
186 3 : const bool BeginHourFlag_reset = state.dataGlobal->BeginHourFlag;
187 3 : const bool EndHourFlag_reset = state.dataGlobal->EndHourFlag;
188 :
189 3 : if (!state.dataWeather->WeatherFileExists) {
190 2 : ShowSevereError(state,
191 : "Site:GroundTemperature:Undisturbed:FiniteDifference -- using this model requires specification of a weather file.");
192 2 : ShowContinueError(
193 : state, "Either place in.epw in the working directory or specify a weather file on the command line using -w /path/to/weather.epw");
194 3 : ShowFatalError(state, "Simulation halted due to input error in ground temperature model.");
195 : }
196 :
197 : // We add a new period to force running all weather data
198 2 : int originalNumOfEnvrn = state.dataWeather->NumOfEnvrn;
199 2 : ++state.dataWeather->NumOfEnvrn;
200 2 : ++state.dataWeather->TotRunPers;
201 2 : state.dataWeather->Environment.redimension(state.dataWeather->NumOfEnvrn);
202 2 : state.dataWeather->RunPeriodInput.redimension(state.dataWeather->TotRunPers);
203 2 : state.dataWeather->Environment(state.dataWeather->NumOfEnvrn).KindOfEnvrn = Constant::KindOfSim::ReadAllWeatherData;
204 2 : state.dataWeather->RPReadAllWeatherData = true;
205 2 : state.dataGlobal->WeathSimReq = true;
206 : // RunPeriod is initialized to be one year of simulation
207 : // RunPeriodInput(TotRunPers).monWeekDay = 0; // Why do this?
208 :
209 2 : Weather::SetupEnvironmentTypes(state);
210 :
211 : // We reset the counter to the original number of run periods, so that GetNextEnvironment will fetch the one we added
212 2 : state.dataWeather->Envrn = originalNumOfEnvrn;
213 2 : bool Available = true;
214 2 : bool ErrorsFound = false;
215 2 : Weather::GetNextEnvironment(state, Available, ErrorsFound);
216 2 : if (ErrorsFound) {
217 0 : ShowFatalError(state, "Site:GroundTemperature:Undisturbed:FiniteDifference: error in reading weather file data");
218 : }
219 :
220 2 : if (state.dataGlobal->KindOfSim != Constant::KindOfSim::ReadAllWeatherData) {
221 : // This shouldn't happen
222 0 : ShowFatalError(state, "Site:GroundTemperature:Undisturbed:FiniteDifference: error in reading weather file data, bad KindOfSim.");
223 : }
224 :
225 2 : weatherDataArray.dimension(state.dataWeather->NumDaysInYear);
226 :
227 2 : state.dataGlobal->BeginEnvrnFlag = true;
228 2 : state.dataGlobal->EndEnvrnFlag = false;
229 2 : state.dataEnvrn->EndMonthFlag = false;
230 2 : state.dataGlobal->WarmupFlag = false;
231 2 : state.dataGlobal->DayOfSim = 0;
232 2 : state.dataGlobal->DayOfSimChr = "0";
233 2 : state.dataReportFlag->NumOfWarmupDays = 0;
234 :
235 2 : Real64 annualAveAirTemp_num = 0.0;
236 :
237 732 : while ((state.dataGlobal->DayOfSim < state.dataWeather->NumDaysInYear) || (state.dataGlobal->WarmupFlag)) { // Begin day loop ...
238 :
239 730 : ++state.dataGlobal->DayOfSim;
240 :
241 : // Reset daily values
242 730 : Real64 outDryBulbTemp_num = 0.0;
243 730 : Real64 relHum_num = 0.0;
244 730 : Real64 windSpeed_num = 0.0;
245 730 : Real64 horizSolarRad_num = 0.0;
246 730 : Real64 airDensity_num = 0.0;
247 730 : int denominator = 0;
248 :
249 730 : auto &tdwd = weatherDataArray(state.dataGlobal->DayOfSim); // "This day weather data"
250 :
251 730 : state.dataGlobal->BeginDayFlag = true;
252 730 : state.dataGlobal->EndDayFlag = false;
253 :
254 18250 : for (state.dataGlobal->HourOfDay = 1; state.dataGlobal->HourOfDay <= 24; ++state.dataGlobal->HourOfDay) { // Begin hour loop ...
255 :
256 17520 : state.dataGlobal->BeginHourFlag = true;
257 17520 : state.dataGlobal->EndHourFlag = false;
258 :
259 87600 : for (state.dataGlobal->TimeStep = 1; state.dataGlobal->TimeStep <= state.dataGlobal->TimeStepsInHour; ++state.dataGlobal->TimeStep) {
260 :
261 70080 : state.dataGlobal->BeginTimeStepFlag = true;
262 :
263 : // Set the End__Flag variables to true if necessary. Note that
264 : // each flag builds on the previous level. EndDayFlag cannot be
265 : // .TRUE. unless EndHourFlag is also .TRUE., etc. Note that the
266 : // EndEnvrnFlag and the EndSimFlag cannot be set during warmup.
267 : // Note also that BeginTimeStepFlag, EndTimeStepFlag, and the
268 : // SubTimeStepFlags can/will be set/reset in the HVAC Manager.
269 :
270 70080 : if (state.dataGlobal->TimeStep == state.dataGlobal->TimeStepsInHour) {
271 17520 : state.dataGlobal->EndHourFlag = true;
272 17520 : if (state.dataGlobal->HourOfDay == 24) {
273 730 : state.dataGlobal->EndDayFlag = true;
274 730 : if (!state.dataGlobal->WarmupFlag && (state.dataGlobal->DayOfSim == state.dataGlobal->NumOfDayInEnvrn)) {
275 2 : state.dataGlobal->EndEnvrnFlag = true;
276 : }
277 : }
278 : }
279 :
280 70080 : Weather::ManageWeather(state);
281 :
282 70080 : outDryBulbTemp_num += state.dataEnvrn->OutDryBulbTemp;
283 70080 : airDensity_num += state.dataEnvrn->OutAirDensity;
284 70080 : relHum_num += state.dataEnvrn->OutRelHumValue;
285 70080 : windSpeed_num += state.dataEnvrn->WindSpeed;
286 70080 : horizSolarRad_num += max(state.dataEnvrn->SOLCOS(3), 0.0) * state.dataEnvrn->BeamSolarRad + state.dataEnvrn->DifSolarRad;
287 :
288 70080 : state.dataGlobal->BeginHourFlag = false;
289 70080 : state.dataGlobal->BeginDayFlag = false;
290 70080 : state.dataGlobal->BeginEnvrnFlag = false;
291 70080 : state.dataGlobal->BeginSimFlag = false;
292 :
293 70080 : ++denominator;
294 :
295 : } // TimeStep loop
296 :
297 17520 : state.dataGlobal->PreviousHour = state.dataGlobal->HourOfDay;
298 :
299 : } // ... End hour loop.
300 :
301 730 : tdwd.dryBulbTemp = outDryBulbTemp_num / denominator;
302 730 : tdwd.relativeHumidity = relHum_num / denominator;
303 730 : tdwd.windSpeed = windSpeed_num / denominator;
304 730 : tdwd.horizontalRadiation = horizSolarRad_num / denominator;
305 730 : tdwd.airDensity = airDensity_num / denominator;
306 :
307 : // Log data for domain initialization using KA model
308 730 : annualAveAirTemp_num += tdwd.dryBulbTemp;
309 :
310 730 : if (tdwd.dryBulbTemp < minDailyAirTemp) {
311 10 : minDailyAirTemp = tdwd.dryBulbTemp;
312 10 : dayOfMinDailyAirTemp = state.dataGlobal->DayOfSim;
313 : }
314 :
315 730 : if (tdwd.dryBulbTemp > maxDailyAirTemp) {
316 32 : maxDailyAirTemp = tdwd.dryBulbTemp;
317 : }
318 :
319 : } // ... End day loop.
320 :
321 2 : annualAveAirTemp = annualAveAirTemp_num / state.dataWeather->NumDaysInYear; // Used for initializing domain
322 :
323 : // Reset Environment when done reading data
324 2 : --state.dataWeather->NumOfEnvrn; // May need better way of eliminating the extra environment that was added to read the data
325 2 : --state.dataWeather->TotRunPers;
326 2 : state.dataGlobal->KindOfSim = KindOfSim_reset;
327 2 : state.dataWeather->RPReadAllWeatherData = false;
328 2 : state.dataWeather->Environment.redimension(state.dataWeather->NumOfEnvrn);
329 2 : state.dataWeather->RunPeriodInput.redimension(state.dataWeather->TotRunPers);
330 2 : state.dataWeather->Envrn = Envrn_reset;
331 2 : state.dataGlobal->TimeStep = TimeStep_reset;
332 2 : state.dataGlobal->HourOfDay = HourOfDay_reset;
333 2 : state.dataGlobal->BeginEnvrnFlag = BeginEnvrnFlag_reset;
334 2 : state.dataGlobal->EndEnvrnFlag = EndEnvrnFlag_reset;
335 2 : state.dataEnvrn->EndMonthFlag = EndMonthFlag_reset;
336 2 : state.dataGlobal->WarmupFlag = WarmupFlag_reset;
337 2 : state.dataGlobal->DayOfSim = DayOfSim_reset;
338 2 : state.dataGlobal->DayOfSimChr = DayOfSimChr_reset;
339 2 : state.dataReportFlag->NumOfWarmupDays = NumOfWarmupDays_reset;
340 2 : state.dataGlobal->BeginDayFlag = BeginDayFlag_reset;
341 2 : state.dataGlobal->EndDayFlag = EndDayFlag_reset;
342 2 : state.dataGlobal->BeginHourFlag = BeginHourFlag_reset;
343 2 : state.dataGlobal->EndHourFlag = EndHourFlag_reset;
344 3 : }
345 :
346 : //******************************************************************************
347 :
348 2 : void FiniteDiffGroundTempsModel::developMesh()
349 : {
350 : // SUBROUTINE INFORMATION:
351 : // AUTHOR Matt Mitchell
352 : // DATE WRITTEN Summer 2015
353 :
354 : // PURPOSE OF THIS SUBROUTINE:
355 : // Creates static mesh used for model
356 :
357 : // Surface layer parameters
358 2 : constexpr Real64 surfaceLayerThickness = 2.0;
359 2 : constexpr Real64 surfaceLayerCellThickness = 0.015;
360 2 : constexpr int surfaceLayerNumCells = static_cast<int>(surfaceLayerThickness / surfaceLayerCellThickness);
361 :
362 : // Center layer parameters
363 2 : constexpr int centerLayerNumCells = 80;
364 :
365 : // Deep layer parameters
366 2 : constexpr Real64 deepLayerThickness = 0.2;
367 2 : constexpr Real64 deepLayerCellThickness = surfaceLayerCellThickness;
368 2 : constexpr int deepLayerNumCells = static_cast<int>(deepLayerThickness / deepLayerCellThickness);
369 :
370 : // Other
371 2 : Real64 currentCellDepth = 0.0;
372 :
373 2 : totalNumCells = surfaceLayerNumCells + centerLayerNumCells + deepLayerNumCells;
374 :
375 : // Allocate arrays
376 2 : cellArray.allocate(totalNumCells);
377 2 : cellDepths.allocate(totalNumCells);
378 :
379 454 : for (int i = 1; i <= totalNumCells; ++i) {
380 :
381 : // Reference to thisCell
382 452 : auto &thisCell = cellArray(i);
383 :
384 : // Set the index
385 452 : thisCell.index = i;
386 :
387 : // Give thickness to the cells
388 452 : if (i <= surfaceLayerNumCells) {
389 : // Constant thickness mesh here
390 266 : thisCell.thickness = surfaceLayerCellThickness;
391 :
392 186 : } else if (i <= (centerLayerNumCells + surfaceLayerNumCells)) {
393 : // Geometric expansion/contraction here
394 160 : int numCenterCell = i - surfaceLayerNumCells;
395 :
396 160 : if (numCenterCell <= (centerLayerNumCells / 2)) {
397 80 : Real64 centerLayerExpansionCoeff = 1.10879;
398 80 : thisCell.thickness = surfaceLayerCellThickness * std::pow(centerLayerExpansionCoeff, numCenterCell);
399 : } else {
400 80 : thisCell.thickness =
401 80 : cellArray((surfaceLayerNumCells + (centerLayerNumCells / 2)) - (numCenterCell - (centerLayerNumCells / 2))).thickness;
402 : }
403 : } else {
404 : // Constant thickness mesh here
405 26 : thisCell.thickness = deepLayerCellThickness;
406 : }
407 :
408 : // Set minimum z value
409 452 : thisCell.minZValue = currentCellDepth;
410 :
411 : // Populate depth array for use later when looking up temperatures
412 452 : cellDepths(i) = currentCellDepth + thisCell.thickness / 2.0;
413 :
414 : // Update local counter
415 452 : currentCellDepth += thisCell.thickness;
416 :
417 : // Set maximum z value
418 452 : thisCell.maxZValue = currentCellDepth;
419 :
420 : // Set base properties
421 452 : thisCell.props.conductivity = baseConductivity;
422 452 : thisCell.props.density = baseDensity;
423 452 : thisCell.props.specificHeat = baseSpecificHeat;
424 452 : thisCell.props.diffusivity = baseConductivity / (baseDensity * baseSpecificHeat);
425 : }
426 2 : }
427 :
428 : //******************************************************************************
429 :
430 2 : void FiniteDiffGroundTempsModel::performSimulation(EnergyPlusData &state)
431 : {
432 : // SUBROUTINE INFORMATION:
433 : // AUTHOR Matt Mitchell
434 : // DATE WRITTEN Summer 2015
435 :
436 : // PURPOSE OF THIS SUBROUTINE:
437 : // Simulates model, repeating years, until steady-periodic temperatures are determined.
438 :
439 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
440 2 : timeStepInSeconds = Constant::rSecsInDay;
441 2 : bool convergedFinal = false;
442 :
443 2 : initDomain(state);
444 :
445 : // Loop until converged
446 : do {
447 :
448 : // loop over all days
449 5124 : for (state.dataGlobal->FDsimDay = 1; state.dataGlobal->FDsimDay <= state.dataWeather->NumDaysInYear; ++state.dataGlobal->FDsimDay) {
450 :
451 5110 : bool iterationConverged = false;
452 :
453 5110 : doStartOfTimeStepInits();
454 :
455 : // Loop until iteration temperature converges
456 : do {
457 :
458 : // For all cells
459 826225293 : for (int cell = 1; cell <= totalNumCells; ++cell) {
460 :
461 822585534 : if (cell == 1) {
462 3639759 : updateSurfaceCellTemperature(state);
463 818945775 : } else if (cell > 1 && cell < totalNumCells) {
464 815306016 : updateGeneralDomainCellTemperature(cell);
465 3639759 : } else if (cell == totalNumCells) {
466 3639759 : updateBottomCellTemperature();
467 : }
468 : }
469 :
470 : // Check iteration temperature convergence
471 3639759 : iterationConverged = checkIterationTemperatureConvergence();
472 :
473 3639759 : if (!iterationConverged) {
474 : // Shift temperatures for next iteration
475 3634649 : updateIterationTemperatures();
476 : }
477 :
478 3639759 : } while (!iterationConverged);
479 :
480 : // Shift temperatures for next timestep
481 5110 : updateTimeStepTemperatures(state);
482 : }
483 :
484 : // Check final temperature convergence
485 14 : convergedFinal = checkFinalTemperatureConvergence(state);
486 :
487 14 : } while (!convergedFinal);
488 2 : }
489 :
490 : //******************************************************************************
491 :
492 3639759 : void FiniteDiffGroundTempsModel::updateSurfaceCellTemperature(const EnergyPlusData &state)
493 : {
494 : // SUBROUTINE INFORMATION:
495 : // AUTHOR Matt Mitchell
496 : // DATE WRITTEN Summer 2015
497 :
498 : // PURPOSE OF THIS SUBROUTINE:
499 : // Determines heat transfer to surface. Updates surface cell temperature.
500 :
501 : // FUNCTION LOCAL VARIABLE DECLARATIONS:
502 3639759 : Real64 numerator = 0.0;
503 3639759 : Real64 denominator = 0.0;
504 3639759 : Real64 resistance = 0.0;
505 : Real64 G_hr;
506 : Real64 Cd;
507 :
508 3639759 : Real64 constexpr rho_water = 998.0; // [kg/m3]
509 : // evapotranspiration parameters
510 3639759 : Real64 constexpr absor_Corrected = 0.77;
511 3639759 : constexpr Real64 convert_Wm2_To_MJhrmin = 3600.0 / 1000000.0;
512 3639759 : constexpr Real64 convert_MJhrmin_To_Wm2 = 1.0 / convert_Wm2_To_MJhrmin;
513 :
514 3639759 : const auto &thisCell = cellArray(1);
515 3639759 : const auto &cellBelow_thisCell = cellArray(2);
516 3639759 : const auto &cwd = weatherDataArray(state.dataGlobal->FDsimDay); // "Current Weather Day"
517 :
518 : // Add effect from previous time step
519 3639759 : numerator += thisCell.temperature_prevTimeStep;
520 3639759 : ++denominator;
521 :
522 : // Conduction to lower cell
523 3639759 : resistance = thisCell.thickness / 2.0 / (thisCell.props.conductivity * thisCell.conductionArea) +
524 3639759 : cellBelow_thisCell.thickness / 2.0 / (cellBelow_thisCell.props.conductivity * cellBelow_thisCell.conductionArea);
525 3639759 : numerator += thisCell.beta / resistance * cellBelow_thisCell.temperature;
526 3639759 : denominator += thisCell.beta / resistance;
527 :
528 : // Convection to atmosphere
529 3639759 : if (cwd.windSpeed > 0.1) {
530 3639759 : Real64 constexpr airSpecificHeat = 1003; // '[J/kg-K]
531 3639759 : resistance = 208.0 / (cwd.airDensity * airSpecificHeat * cwd.windSpeed * thisCell.conductionArea);
532 : } else {
533 : // Future development should include additional natural convection effects here
534 : }
535 3639759 : numerator += thisCell.beta / resistance * cwd.dryBulbTemp;
536 3639759 : denominator += thisCell.beta / resistance;
537 :
538 : // For convenience convert to Kelvin once
539 3639759 : const Real64 currAirTempK = cwd.dryBulbTemp + 273.15;
540 :
541 : // Convert input solar radiation [w/m2] into units for ET model, [MJ/hr-min]
542 : // Diffuse + Direct Beam Radiation
543 3639759 : const Real64 incidentSolar_MJhrmin = cwd.horizontalRadiation * convert_Wm2_To_MJhrmin;
544 :
545 : // Absorbed solar radiation, [MJ/hr-min]
546 3639759 : const Real64 absorbedIncidentSolar_MJhrmin = absor_Corrected * incidentSolar_MJhrmin;
547 :
548 : // Calculate saturated vapor pressure, [kPa]
549 3639759 : const Real64 vaporPressureSaturated_kPa = 0.6108 * std::exp(17.27 * cwd.dryBulbTemp / (cwd.dryBulbTemp + 237.3));
550 :
551 : // Calculate actual vapor pressure, [kPa]
552 3639759 : const Real64 vaporPressureActual_kPa = vaporPressureSaturated_kPa * cwd.relativeHumidity;
553 :
554 : // Calculate another Q term, [MJ/m2-hr]
555 3639759 : const Real64 QRAD_NL = 2.042E-10 * pow_4(currAirTempK) * (0.34 - 0.14 * std::sqrt(vaporPressureActual_kPa));
556 :
557 : // Calculate another Q term, [MJ/hr]
558 3639759 : const Real64 netIncidentRadiation_MJhr = absorbedIncidentSolar_MJhrmin - QRAD_NL;
559 :
560 : // constant
561 3639759 : constexpr Real64 CN = 37.0;
562 :
563 : // Check whether there was sun
564 3639759 : if (netIncidentRadiation_MJhr < 0.0) {
565 713385 : G_hr = 0.5 * netIncidentRadiation_MJhr;
566 713385 : Cd = 0.96;
567 : } else {
568 2926374 : G_hr = 0.1 * netIncidentRadiation_MJhr;
569 2926374 : Cd = 0.24;
570 : }
571 :
572 3639759 : const Real64 slope_S = 2503.0 * std::exp(17.27 * cwd.dryBulbTemp / (cwd.dryBulbTemp + 237.3)) / pow_2(cwd.dryBulbTemp + 237.3);
573 3639759 : constexpr Real64 pressure = 98.0;
574 3639759 : constexpr Real64 psychrometricConstant = 0.665e-3 * pressure;
575 :
576 : // Evapotranspiration constant, [mm/hr]
577 3639759 : const Real64 evapotransFluidLoss_mmhr =
578 3639759 : (evapotransCoeff * slope_S * (netIncidentRadiation_MJhr - G_hr) +
579 3639759 : psychrometricConstant * (CN / currAirTempK) * cwd.windSpeed * (vaporPressureSaturated_kPa - vaporPressureActual_kPa)) /
580 3639759 : (slope_S + psychrometricConstant * (1 + Cd * cwd.windSpeed));
581 :
582 : // Convert units, [m/hr]
583 3639759 : const Real64 evapotransFluidLoss_mhr = evapotransFluidLoss_mmhr / 1000.0;
584 :
585 : // Calculate latent heat, [MJ/kg]
586 : // Full formulation is cubic: L(T) = -0.0000614342 * T**3 + 0.00158927 * T**2 - 2.36418 * T + 2500.79[5]
587 : // In: Cubic fit to Table 2.1,p.16, Textbook: R.R.Rogers & M.K. Yau, A Short Course in Cloud Physics, 3e,(1989), Pergamon press
588 : // But a linear relation should suffice;
589 : // note-for now using the previous time step temperature as an approximation to help ensure stability
590 3639759 : const Real64 latentHeatVaporization = 2.501 - 2.361e-3 * thisCell.temperature_prevTimeStep;
591 :
592 : // Calculate evapotranspiration heat loss, [MJ/m2-hr]
593 3639759 : const Real64 evapotransHeatLoss_MJhrmin = rho_water * evapotransFluidLoss_mhr * latentHeatVaporization;
594 :
595 : // Convert net incident solar units, [W/m2]
596 3639759 : const Real64 netIncidentRadiation_Wm2 = netIncidentRadiation_MJhr * convert_MJhrmin_To_Wm2;
597 :
598 : // Convert evapotranspiration units, [W/m2]
599 3639759 : const Real64 evapotransHeatLoss_Wm2 = evapotransHeatLoss_MJhrmin * convert_MJhrmin_To_Wm2;
600 :
601 : // Calculate overall net heat ?gain? into the cell, [W]
602 3639759 : const Real64 incidentHeatGain = (netIncidentRadiation_Wm2 - evapotransHeatLoss_Wm2) * thisCell.conductionArea;
603 :
604 : // Add any solar/evapotranspiration heat gain here
605 3639759 : numerator += thisCell.beta * incidentHeatGain;
606 :
607 : // Calculate the return temperature and leave
608 3639759 : cellArray(1).temperature = numerator / denominator;
609 3639759 : }
610 :
611 : //******************************************************************************
612 :
613 815306016 : void FiniteDiffGroundTempsModel::updateGeneralDomainCellTemperature(int const cell)
614 : {
615 : // SUBROUTINE INFORMATION:
616 : // AUTHOR Matt Mitchell
617 : // DATE WRITTEN Summer 2015
618 :
619 : // PURPOSE OF THIS SUBROUTINE:
620 : // Update cell temperature based on HT from cells above and below
621 :
622 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
623 815306016 : Real64 numerator = 0.0;
624 815306016 : Real64 denominator = 0.0;
625 815306016 : Real64 resistance = 0.0;
626 :
627 815306016 : auto &thisCell = cellArray(cell);
628 815306016 : const auto &cellAbove_thisCell = cellArray(cell - 1);
629 815306016 : const auto &cellBelow_thisCell = cellArray(cell + 1);
630 :
631 : // add effect from cell history
632 815306016 : numerator += thisCell.temperature_prevTimeStep;
633 815306016 : ++denominator;
634 :
635 : // Conduction resistance between this cell and above cell
636 815306016 : resistance = thisCell.thickness / 2.0 / (thisCell.conductionArea * thisCell.props.conductivity) +
637 815306016 : cellAbove_thisCell.thickness / 2.0 / (cellAbove_thisCell.conductionArea * cellAbove_thisCell.props.conductivity);
638 :
639 815306016 : numerator += thisCell.beta / resistance * cellAbove_thisCell.temperature;
640 815306016 : denominator += thisCell.beta / resistance;
641 :
642 : // Conduction resistance between this cell and below cell
643 815306016 : resistance = thisCell.thickness / 2.0 / (thisCell.conductionArea * thisCell.props.conductivity) +
644 815306016 : cellBelow_thisCell.thickness / 2.0 / (cellBelow_thisCell.conductionArea * cellBelow_thisCell.props.conductivity);
645 :
646 815306016 : numerator += thisCell.beta / resistance * cellBelow_thisCell.temperature;
647 815306016 : denominator += thisCell.beta / resistance;
648 :
649 : // now that we have passed all directions, update the temperature
650 815306016 : thisCell.temperature = numerator / denominator;
651 815306016 : }
652 :
653 : //******************************************************************************
654 :
655 3639759 : void FiniteDiffGroundTempsModel::updateBottomCellTemperature()
656 : {
657 : // SUBROUTINE INFORMATION:
658 : // AUTHOR Matt Mitchell
659 : // DATE WRITTEN Summer 2015
660 :
661 : // PURPOSE OF THIS SUBROUTINE:
662 : // Updates bottom cell temperature based on earth heat flux HT from cell above
663 :
664 : // REFERENCES:
665 : // Fridleifsson, I.B., R. Bertani, E.Huenges, J.W. Lund, A. Ragnarsson, L. Rybach. 2008
666 : // 'The possible role and contribution of geothermal energy to the mitigation of climate change.'
667 : // IPCC scoping meeting on renewable energy sources: 59-80.
668 :
669 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
670 3639759 : Real64 numerator = 0.0;
671 3639759 : Real64 denominator = 0.0;
672 3639759 : Real64 resistance = 0.0;
673 3639759 : Real64 constexpr geothermalGradient = 0.025; // C/m
674 :
675 3639759 : auto &thisCell = cellArray(totalNumCells);
676 3639759 : auto &cellAbove_thisCell = cellArray(totalNumCells - 1);
677 :
678 : // Initialize
679 3639759 : numerator += thisCell.temperature_prevTimeStep;
680 3639759 : ++denominator;
681 :
682 : // Conduction resistance between this cell and above cell
683 3639759 : resistance = ((thisCell.thickness / 2.0) / (thisCell.conductionArea * thisCell.props.conductivity)) +
684 3639759 : ((cellAbove_thisCell.thickness / 2.0) / (cellAbove_thisCell.conductionArea * cellAbove_thisCell.props.conductivity));
685 :
686 3639759 : numerator += (thisCell.beta / resistance) * cellAbove_thisCell.temperature;
687 3639759 : denominator += thisCell.beta / resistance;
688 :
689 : // Geothermal gradient heat transfer
690 3639759 : Real64 HTBottom = geothermalGradient * thisCell.props.conductivity * thisCell.conductionArea;
691 :
692 3639759 : numerator += thisCell.beta * HTBottom;
693 :
694 3639759 : cellArray(totalNumCells).temperature = numerator / denominator;
695 3639759 : }
696 :
697 : //******************************************************************************
698 :
699 14 : bool FiniteDiffGroundTempsModel::checkFinalTemperatureConvergence(const EnergyPlusData &state)
700 : {
701 : // SUBROUTINE INFORMATION:
702 : // AUTHOR Matt Mitchell
703 : // DATE WRITTEN Summer 2015
704 :
705 : // PURPOSE OF THIS SUBROUTINE:
706 : // Checks final temperature convergence
707 :
708 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
709 14 : bool converged = true;
710 14 : Real64 constexpr finalTempConvergenceCriteria = 0.05;
711 :
712 14 : if (state.dataGlobal->FDnumIterYears == maxYearsToIterate) return converged;
713 :
714 2951 : for (int cell = 1; cell <= totalNumCells; ++cell) {
715 :
716 2938 : auto &thisCell = cellArray(cell);
717 :
718 2938 : if (std::abs(thisCell.temperature - thisCell.temperature_finalConvergence) >= finalTempConvergenceCriteria) {
719 1021 : converged = false;
720 : }
721 :
722 2938 : thisCell.temperature_finalConvergence = thisCell.temperature;
723 : }
724 :
725 13 : ++state.dataGlobal->FDnumIterYears;
726 :
727 13 : return converged;
728 : }
729 :
730 : //******************************************************************************
731 :
732 3639759 : bool FiniteDiffGroundTempsModel::checkIterationTemperatureConvergence()
733 : {
734 : // SUBROUTINE INFORMATION:
735 : // AUTHOR Matt Mitchell
736 : // DATE WRITTEN Summer 2015
737 :
738 : // PURPOSE OF THIS SUBROUTINE:
739 : // Checks iteration temperature convergence
740 :
741 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
742 3639759 : bool converged = true;
743 3639759 : Real64 constexpr iterationTempConvergenceCriteria = 0.00001;
744 :
745 29476507 : for (int cell = 1; cell <= totalNumCells; ++cell) {
746 :
747 29471397 : if (std::abs(cellArray(cell).temperature - cellArray(cell).temperature_prevIteration) >= iterationTempConvergenceCriteria) {
748 3634649 : converged = false;
749 3634649 : break;
750 : }
751 : }
752 :
753 3639759 : return converged;
754 : }
755 :
756 : //******************************************************************************
757 :
758 2 : void FiniteDiffGroundTempsModel::initDomain(EnergyPlusData &state)
759 : {
760 : // SUBROUTINE INFORMATION:
761 : // AUTHOR Matt Mitchell
762 : // DATE WRITTEN Summer 2015
763 :
764 : // PURPOSE OF THIS SUBROUTINE:
765 : // Initializes model using Kusuda-Achenbach model.
766 : // Average ground temp initialized to average annual air temperature
767 :
768 : // Temporary KA model for initialization
769 2 : auto tempModel = std::make_unique<KusudaGroundTempsModel>(); // (AUTO_OK) Why does this have to be a unique_ptr?
770 :
771 2 : tempModel->Name = "KAModelForFDModel";
772 2 : tempModel->modelType = GroundTemp::ModelType::Kusuda;
773 2 : tempModel->aveGroundTemp = annualAveAirTemp;
774 4 : tempModel->aveGroundTempAmplitude =
775 2 : (maxDailyAirTemp - minDailyAirTemp) / 4.0; // Rough estimate here. Ground temps will not swing as far as the air temp.
776 2 : tempModel->phaseShiftInSecs = dayOfMinDailyAirTemp * Constant::rSecsInDay;
777 2 : tempModel->groundThermalDiffusivity = baseConductivity / (baseDensity * baseSpecificHeat);
778 :
779 : // Initialize temperatures and volume
780 454 : for (int cell = 1; cell <= totalNumCells; ++cell) {
781 452 : auto &thisCell = cellArray(cell);
782 :
783 452 : Real64 depth = (thisCell.maxZValue + thisCell.minZValue) / 2.0;
784 :
785 : // Initialize temperatures
786 452 : if (tempModel) {
787 452 : thisCell.temperature = tempModel->getGroundTempAtTimeInSeconds(state, depth, 0.0); // Initialized at first day of year
788 : }
789 452 : thisCell.temperature_finalConvergence = thisCell.temperature;
790 452 : thisCell.temperature_prevIteration = thisCell.temperature;
791 452 : thisCell.temperature_prevTimeStep = thisCell.temperature;
792 :
793 : // Set cell volume
794 452 : thisCell.volume = thisCell.thickness * thisCell.conductionArea;
795 : }
796 :
797 : // Initialize freezing calculation variables
798 2 : evaluateSoilRhoCpInit();
799 :
800 : // Initialize the groundTemps array
801 2 : groundTemps.dimension({1, state.dataWeather->NumDaysInYear}, {1, totalNumCells}, 0.0);
802 :
803 2 : tempModel.reset();
804 2 : }
805 :
806 : //******************************************************************************
807 :
808 3634649 : void FiniteDiffGroundTempsModel::updateIterationTemperatures()
809 : {
810 : // SUBROUTINE INFORMATION:
811 : // AUTHOR Matt Mitchell
812 : // DATE WRITTEN Summer 2015
813 :
814 : // PURPOSE OF THIS SUBROUTINE:
815 : // Updates iteration temperatures for convergence checks
816 :
817 825065323 : for (int cell = 1; cell <= totalNumCells; ++cell) {
818 821430674 : cellArray(cell).temperature_prevIteration = cellArray(cell).temperature;
819 : }
820 3634649 : }
821 :
822 : //******************************************************************************
823 :
824 5110 : void FiniteDiffGroundTempsModel::updateTimeStepTemperatures(const EnergyPlusData &state)
825 : {
826 : // SUBROUTINE INFORMATION:
827 : // AUTHOR Matt Mitchell
828 : // DATE WRITTEN Summer 2015
829 :
830 : // PURPOSE OF THIS SUBROUTINE:
831 : // Updates timestep temperatures for convergence checks.
832 :
833 1159970 : for (int cell = 1; cell <= totalNumCells; ++cell) {
834 :
835 1154860 : auto &thisCell = cellArray(cell);
836 :
837 1154860 : thisCell.temperature_prevTimeStep = thisCell.temperature;
838 :
839 : // Log temps for later use
840 1154860 : groundTemps(state.dataGlobal->FDsimDay, cell) = thisCell.temperature;
841 : }
842 5110 : }
843 :
844 : //******************************************************************************
845 :
846 5110 : void FiniteDiffGroundTempsModel::doStartOfTimeStepInits()
847 : {
848 : // SUBROUTINE INFORMATION:
849 : // AUTHOR Matt Mitchell
850 : // DATE WRITTEN Summer 2015
851 :
852 : // PURPOSE OF THIS SUBROUTINE:
853 : // Updates cell properties for each timestep
854 :
855 1159970 : for (int cell = 1; cell <= totalNumCells; ++cell) {
856 :
857 1154860 : auto &thisCell = cellArray(cell);
858 :
859 1154860 : evaluateSoilRhoCpCell(cell);
860 :
861 1154860 : thisCell.beta = (timeStepInSeconds / (thisCell.props.rhoCp * thisCell.volume));
862 : }
863 5110 : }
864 :
865 : //******************************************************************************
866 :
867 49 : Real64 FiniteDiffGroundTempsModel::interpolate(Real64 const x, Real64 const x_hi, Real64 const x_low, Real64 const y_hi, Real64 const y_low)
868 : {
869 49 : return (x - x_low) / (x_hi - x_low) * (y_hi - y_low) + y_low;
870 : }
871 :
872 : //******************************************************************************
873 :
874 20 : Real64 FiniteDiffGroundTempsModel::getGroundTemp(EnergyPlusData &state)
875 : {
876 :
877 : // SUBROUTINE INFORMATION:
878 : // AUTHOR Matt Mitchell
879 : // DATE WRITTEN Summer 2015
880 :
881 : // PURPOSE OF THIS SUBROUTINE:
882 : // Interpolates between days and depths to find correct ground temperature
883 :
884 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
885 : // Interpolation variables
886 : int i0; // First day
887 : int i1; // Next day
888 : int j1; // Next cell index (with depth greater than y-depth
889 : // Real64 T_i0_j0; // Temp at int( x-day ); cell lower_bound( y-depth )
890 : // Real64 T_i1_j0; // Temp at int( x-day ) + 1; cell lower_bound( y-depth )
891 : // Real64 T_i0_j1; // Temp at int( x-day ); cell lower_bound( y-depth ) + 1
892 : // Real64 T_i1_j1; // Temp at int( x-day ) + 1; cell lower_bound( y-depth ) + 1
893 : // Real64 T_ix_j0; // Temp at x-day; cell lower_bound( y-depth )
894 : // Real64 T_ix_j1; // Temp at x-day; cell lower_bound( y-depth ) + 1
895 : // RETURNS: Final Temperature--Temp at x-day; y-depth
896 :
897 20 : if (depth < 0.0) {
898 0 : depth = 0.0;
899 : }
900 :
901 : // Get index of nearest cell with depth less than depth
902 20 : auto it = std::lower_bound(cellDepths.begin(), cellDepths.end(), depth);
903 20 : int j0 = static_cast<int>(std::distance(cellDepths.begin(), it)); // Cell index with depth less than y-depth
904 :
905 : // Compensate for 1-based array
906 20 : ++j0;
907 :
908 : // Fraction of day
909 20 : const Real64 dayFrac = simTimeInDays - static_cast<int>(simTimeInDays); // Fraction of day
910 :
911 20 : if (j0 < totalNumCells - 1) {
912 : // All depths within domain
913 14 : j1 = j0 + 1;
914 :
915 14 : if (simTimeInDays <= 1 || simTimeInDays >= state.dataWeather->NumDaysInYear) {
916 : // First day of year, last day of year, and leap day
917 : // Interpolate between first and last day
918 1 : i0 = state.dataWeather->NumDaysInYear;
919 1 : i1 = 1;
920 :
921 : // Lookup ground temps
922 1 : const Real64 T_i0_j0 = groundTemps(i0, j0);
923 1 : const Real64 T_i0_j1 = groundTemps(i0, j1);
924 1 : const Real64 T_i1_j0 = groundTemps(i1, j0);
925 1 : const Real64 T_i1_j1 = groundTemps(i1, j1);
926 :
927 : // Interpolate between days holding depth constant
928 1 : const Real64 T_ix_j0 = interpolate(dayFrac, 1, 0, T_i1_j0, T_i0_j0);
929 1 : const Real64 T_ix_j1 = interpolate(dayFrac, 1, 0, T_i1_j1, T_i0_j1);
930 :
931 : // Interpolate to correct depth now that we're at the right time
932 1 : return interpolate(depth, cellDepths(j1), cellDepths(j0), T_ix_j1, T_ix_j0);
933 : }
934 :
935 : // All other days
936 13 : i0 = static_cast<int>(simTimeInDays);
937 13 : i1 = i0 + 1;
938 :
939 : // Lookup ground temps
940 13 : const Real64 T_i0_j0 = groundTemps(i0, j0);
941 13 : const Real64 T_i0_j1 = groundTemps(i0, j1);
942 13 : const Real64 T_i1_j0 = groundTemps(i1, j0);
943 13 : const Real64 T_i1_j1 = groundTemps(i1, j1);
944 :
945 : // Interpolate between days holding depth constant
946 13 : const Real64 T_ix_j0 = interpolate(dayFrac, 1, 0, T_i1_j0, T_i0_j0);
947 13 : const Real64 T_ix_j1 = interpolate(dayFrac, 1, 0, T_i1_j1, T_i0_j1);
948 :
949 : // Interpolate to correct depth now that we're at the right time
950 13 : return interpolate(depth, cellDepths(j1), cellDepths(j0), T_ix_j1, T_ix_j0);
951 : }
952 :
953 : // Requesting a temperature deeper than domain. Pass deepest point in domain.
954 6 : j0 = totalNumCells;
955 6 : j1 = j0;
956 :
957 6 : if (simTimeInDays <= 1 || simTimeInDays >= state.dataWeather->NumDaysInYear) {
958 : // First day of year, last day of year, and leap day
959 : // Interpolate between first and last day
960 1 : i0 = state.dataWeather->NumDaysInYear;
961 1 : i1 = 1;
962 :
963 : // Lookup ground temps
964 1 : const Real64 T_i0_j1 = groundTemps(i0, j1);
965 1 : const Real64 T_i1_j1 = groundTemps(i1, j1);
966 :
967 : // Interpolate between days holding depth constant
968 1 : return interpolate(dayFrac, 1, 0, T_i1_j1, T_i0_j1);
969 : }
970 :
971 : // All other days
972 5 : i0 = static_cast<int>(simTimeInDays);
973 5 : i1 = i0 + 1;
974 :
975 : // Lookup ground temps
976 5 : const Real64 T_i0_j1 = groundTemps(i0, j1);
977 5 : const Real64 T_i1_j1 = groundTemps(i1, j1);
978 :
979 : // Interpolate between days holding depth constant
980 5 : return interpolate(dayFrac, 1, 0, T_i1_j1, T_i0_j1);
981 : }
982 :
983 : //******************************************************************************
984 :
985 10 : Real64 FiniteDiffGroundTempsModel::getGroundTempAtTimeInSeconds(EnergyPlusData &state, Real64 const _depth, Real64 const seconds)
986 : {
987 : // SUBROUTINE INFORMATION:
988 : // AUTHOR Matt Mitchell
989 : // DATE WRITTEN Summer 2015
990 :
991 : // PURPOSE OF THIS SUBROUTINE:
992 : // Retrieves ground temperature when input time is in seconds
993 :
994 10 : depth = _depth;
995 :
996 10 : simTimeInDays = seconds / Constant::rSecsInDay;
997 :
998 10 : if (simTimeInDays > state.dataWeather->NumDaysInYear) {
999 1 : simTimeInDays = remainder(simTimeInDays, state.dataWeather->NumDaysInYear);
1000 : }
1001 :
1002 10 : return getGroundTemp(state);
1003 : }
1004 :
1005 : //******************************************************************************
1006 :
1007 10 : Real64 FiniteDiffGroundTempsModel::getGroundTempAtTimeInMonths(EnergyPlusData &state, Real64 const _depth, int const month)
1008 : {
1009 : // SUBROUTINE INFORMATION:
1010 : // AUTHOR Matt Mitchell
1011 : // DATE WRITTEN Summer 2015
1012 :
1013 : // PURPOSE OF THIS SUBROUTINE:
1014 : // Returns ground temperature when input time is in months
1015 :
1016 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
1017 : // TODO: Fixing this to be floating point 12.0 causes diffs and failed tests
1018 10 : Real64 const aveDaysInMonth = state.dataWeather->NumDaysInYear / 12;
1019 :
1020 10 : depth = _depth;
1021 :
1022 : // Convert months to days. Puts time in middle of specified month
1023 10 : simTimeInDays = aveDaysInMonth * ((month - 1) + 0.5);
1024 :
1025 10 : if (simTimeInDays > state.dataWeather->NumDaysInYear) {
1026 1 : simTimeInDays = remainder(simTimeInDays, state.dataWeather->NumDaysInYear);
1027 : }
1028 :
1029 : // Get and return ground temperature
1030 10 : return getGroundTemp(state);
1031 : }
1032 :
1033 : //******************************************************************************
1034 :
1035 1154860 : void FiniteDiffGroundTempsModel::evaluateSoilRhoCpCell(int const cell)
1036 : {
1037 : // SUBROUTINE INFORMATION:
1038 : // AUTHOR Edwin Lee
1039 : // DATE WRITTEN Summer 2011
1040 :
1041 : // PURPOSE OF THIS SUBROUTINE:
1042 : // Evaluates the soil properties on a single cell
1043 :
1044 : // Real64 rhoCP_soil;
1045 :
1046 : // These vary by domain now, so we must be careful to retrieve them every time
1047 :
1048 1154860 : auto &thisCell = cellArray(cell);
1049 :
1050 : // set some temperatures here for generalization -- these could be set in the input file
1051 : // constexpr Real64 frzAllIce = -0.5;
1052 : // constexpr Real64 frzIceTrans = -0.4;
1053 : // constexpr Real64 frzLiqTrans = -0.1;
1054 : // constexpr Real64 frzAllLiq = 0.0;
1055 :
1056 : // calculate this cell's new Cp value based on the cell temperature
1057 : // if (thisCell.temperature >= frzAllLiq) {
1058 : // rhoCP_soil = rhoCp_soil_liq_1;
1059 : // } else if (thisCell.temperature <= frzAllIce) {
1060 : // rhoCP_soil = rhoCP_soil_ice;
1061 : // } else if (thisCell.temperature > frzLiqTrans) {
1062 : // rhoCP_soil = rhoCp_soil_liq_1 + (rhoCP_soil_transient - rhoCP_soil_liq) / (frzAllLiq - frzLiqTrans) * (frzAllLiq -
1063 : // thisCell.temperature);
1064 : // } else if (thisCell.temperature >= frzIceTrans) {
1065 : // rhoCP_soil = rhoCP_soil_transient;
1066 : // } else {
1067 : // rhoCP_soil = rhoCP_soil_ice + (rhoCP_soil_transient - rhoCP_soil_ice) / (frzIceTrans - frzAllIce) * (thisCell.temperature - frzAllIce);
1068 : // }
1069 :
1070 : // TODO: The calculated rhoCP_soil is commented on this line and never used. Curious.
1071 1154860 : thisCell.props.rhoCp = baseDensity * baseSpecificHeat; // rhoCP_soil;
1072 :
1073 1154860 : thisCell.props.specificHeat = thisCell.props.rhoCp / thisCell.props.density;
1074 1154860 : }
1075 :
1076 2 : void FiniteDiffGroundTempsModel::evaluateSoilRhoCpInit()
1077 : {
1078 : // SUBROUTINE INFORMATION:
1079 : // AUTHOR Edwin Lee
1080 : // DATE WRITTEN Summer 2011
1081 :
1082 : // PURPOSE OF THIS SUBROUTINE:
1083 : // Evaluates the soil properties
1084 :
1085 : // These vary by domain now, so we must be careful to retrieve them every time
1086 2 : const Real64 Theta_liq = waterContent;
1087 2 : const Real64 Theta_sat = saturatedWaterContent;
1088 :
1089 : // Assumption
1090 2 : const Real64 Theta_ice = Theta_liq;
1091 :
1092 : //'Cp (freezing) calculations
1093 2 : constexpr Real64 rho_ice = 917.0; //'Kg / m3
1094 2 : constexpr Real64 rho_liq = 1000.0; //'kg / m3
1095 2 : rhoCp_soil_liq_1 = 1225000.0 / (1.0 - Theta_sat); // J/m3K
1096 : // from( " An improved model for predicting soil thermal conductivity from water content at room temperature, Fig 4" )
1097 2 : constexpr Real64 CP_liq = 4180.0; //'J / KgK
1098 2 : constexpr Real64 CP_ice = 2066.0; //'J / KgK
1099 2 : constexpr Real64 Lat_fus = 334000.0; //'J / Kg
1100 2 : constexpr Real64 Cp_transient = Lat_fus / 0.4 + (0.5 * CP_ice - (CP_liq + CP_ice) / 2.0 * 0.1) / 0.4;
1101 : // from( " Numerical and experimental investigation of melting and freezing processes in phase change material storage" )
1102 2 : rhoCP_soil_liq = rhoCp_soil_liq_1 * (1.0 - Theta_sat) + rho_liq * CP_liq * Theta_liq;
1103 2 : rhoCP_soil_transient = rhoCp_soil_liq_1 * (1.0 - Theta_sat) + ((rho_liq + rho_ice) / 2.0) * Cp_transient * Theta_ice;
1104 2 : rhoCP_soil_ice = rhoCp_soil_liq_1 * (1.0 - Theta_sat) + rho_ice * CP_ice * Theta_ice; //'!J / m3K
1105 2 : }
1106 :
1107 : //******************************************************************************
1108 :
1109 : } // namespace GroundTemp
1110 : } // namespace EnergyPlus
|
