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 1 : 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 1 : const int Envrn_reset = state.dataWeather->Envrn;
174 1 : const Constant::KindOfSim KindOfSim_reset = state.dataGlobal->KindOfSim;
175 1 : const int TimeStep_reset = state.dataGlobal->TimeStep;
176 1 : const int HourOfDay_reset = state.dataGlobal->HourOfDay;
177 1 : const bool BeginEnvrnFlag_reset = state.dataGlobal->BeginEnvrnFlag;
178 1 : const bool EndEnvrnFlag_reset = state.dataGlobal->EndEnvrnFlag;
179 1 : const bool EndMonthFlag_reset = state.dataEnvrn->EndMonthFlag;
180 1 : const bool WarmupFlag_reset = state.dataGlobal->WarmupFlag;
181 1 : const int DayOfSim_reset = state.dataGlobal->DayOfSim;
182 1 : const std::string DayOfSimChr_reset = state.dataGlobal->DayOfSimChr;
183 1 : const int NumOfWarmupDays_reset = state.dataReportFlag->NumOfWarmupDays;
184 1 : const bool BeginDayFlag_reset = state.dataGlobal->BeginDayFlag;
185 1 : const bool EndDayFlag_reset = state.dataGlobal->EndDayFlag;
186 1 : const bool BeginHourFlag_reset = state.dataGlobal->BeginHourFlag;
187 1 : const bool EndHourFlag_reset = state.dataGlobal->EndHourFlag;
188 :
189 1 : if (!state.dataWeather->WeatherFileExists) {
190 0 : ShowSevereError(state,
191 : "Site:GroundTemperature:Undisturbed:FiniteDifference -- using this model requires specification of a weather file.");
192 0 : 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 0 : 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 1 : int originalNumOfEnvrn = state.dataWeather->NumOfEnvrn;
199 1 : ++state.dataWeather->NumOfEnvrn;
200 1 : ++state.dataWeather->TotRunPers;
201 1 : state.dataWeather->Environment.redimension(state.dataWeather->NumOfEnvrn);
202 1 : state.dataWeather->RunPeriodInput.redimension(state.dataWeather->TotRunPers);
203 1 : state.dataWeather->Environment(state.dataWeather->NumOfEnvrn).KindOfEnvrn = Constant::KindOfSim::ReadAllWeatherData;
204 1 : state.dataWeather->RPReadAllWeatherData = true;
205 1 : state.dataGlobal->WeathSimReq = true;
206 : // RunPeriod is initialized to be one year of simulation
207 : // RunPeriodInput(TotRunPers).monWeekDay = 0; // Why do this?
208 :
209 1 : 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 1 : state.dataWeather->Envrn = originalNumOfEnvrn;
213 1 : bool Available = true;
214 1 : bool ErrorsFound = false;
215 1 : Weather::GetNextEnvironment(state, Available, ErrorsFound);
216 1 : if (ErrorsFound) {
217 0 : ShowFatalError(state, "Site:GroundTemperature:Undisturbed:FiniteDifference: error in reading weather file data");
218 : }
219 :
220 1 : 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 1 : weatherDataArray.dimension(state.dataWeather->NumDaysInYear);
226 :
227 1 : state.dataGlobal->BeginEnvrnFlag = true;
228 1 : state.dataGlobal->EndEnvrnFlag = false;
229 1 : state.dataEnvrn->EndMonthFlag = false;
230 1 : state.dataGlobal->WarmupFlag = false;
231 1 : state.dataGlobal->DayOfSim = 0;
232 1 : state.dataGlobal->DayOfSimChr = "0";
233 1 : state.dataReportFlag->NumOfWarmupDays = 0;
234 :
235 1 : Real64 annualAveAirTemp_num = 0.0;
236 :
237 366 : while ((state.dataGlobal->DayOfSim < state.dataWeather->NumDaysInYear) || (state.dataGlobal->WarmupFlag)) { // Begin day loop ...
238 :
239 365 : ++state.dataGlobal->DayOfSim;
240 :
241 : // Reset daily values
242 365 : Real64 outDryBulbTemp_num = 0.0;
243 365 : Real64 relHum_num = 0.0;
244 365 : Real64 windSpeed_num = 0.0;
245 365 : Real64 horizSolarRad_num = 0.0;
246 365 : Real64 airDensity_num = 0.0;
247 365 : int denominator = 0;
248 :
249 365 : auto &tdwd = weatherDataArray(state.dataGlobal->DayOfSim); // "This day weather data"
250 :
251 365 : state.dataGlobal->BeginDayFlag = true;
252 365 : state.dataGlobal->EndDayFlag = false;
253 :
254 9125 : for (state.dataGlobal->HourOfDay = 1; state.dataGlobal->HourOfDay <= 24; ++state.dataGlobal->HourOfDay) { // Begin hour loop ...
255 :
256 8760 : state.dataGlobal->BeginHourFlag = true;
257 8760 : state.dataGlobal->EndHourFlag = false;
258 :
259 17520 : for (state.dataGlobal->TimeStep = 1; state.dataGlobal->TimeStep <= state.dataGlobal->TimeStepsInHour; ++state.dataGlobal->TimeStep) {
260 :
261 8760 : 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 8760 : if (state.dataGlobal->TimeStep == state.dataGlobal->TimeStepsInHour) {
271 8760 : state.dataGlobal->EndHourFlag = true;
272 8760 : if (state.dataGlobal->HourOfDay == 24) {
273 365 : state.dataGlobal->EndDayFlag = true;
274 365 : if (!state.dataGlobal->WarmupFlag && (state.dataGlobal->DayOfSim == state.dataGlobal->NumOfDayInEnvrn)) {
275 1 : state.dataGlobal->EndEnvrnFlag = true;
276 : }
277 : }
278 : }
279 :
280 8760 : Weather::ManageWeather(state);
281 :
282 8760 : outDryBulbTemp_num += state.dataEnvrn->OutDryBulbTemp;
283 8760 : airDensity_num += state.dataEnvrn->OutAirDensity;
284 8760 : relHum_num += state.dataEnvrn->OutRelHumValue;
285 8760 : windSpeed_num += state.dataEnvrn->WindSpeed;
286 8760 : horizSolarRad_num += max(state.dataEnvrn->SOLCOS(3), 0.0) * state.dataEnvrn->BeamSolarRad + state.dataEnvrn->DifSolarRad;
287 :
288 8760 : state.dataGlobal->BeginHourFlag = false;
289 8760 : state.dataGlobal->BeginDayFlag = false;
290 8760 : state.dataGlobal->BeginEnvrnFlag = false;
291 8760 : state.dataGlobal->BeginSimFlag = false;
292 :
293 8760 : ++denominator;
294 :
295 : } // TimeStep loop
296 :
297 8760 : state.dataGlobal->PreviousHour = state.dataGlobal->HourOfDay;
298 :
299 : } // ... End hour loop.
300 :
301 365 : tdwd.dryBulbTemp = outDryBulbTemp_num / denominator;
302 365 : tdwd.relativeHumidity = relHum_num / denominator;
303 365 : tdwd.windSpeed = windSpeed_num / denominator;
304 365 : tdwd.horizontalRadiation = horizSolarRad_num / denominator;
305 365 : tdwd.airDensity = airDensity_num / denominator;
306 :
307 : // Log data for domain initialization using KA model
308 365 : annualAveAirTemp_num += tdwd.dryBulbTemp;
309 :
310 365 : if (tdwd.dryBulbTemp < minDailyAirTemp) {
311 7 : minDailyAirTemp = tdwd.dryBulbTemp;
312 7 : dayOfMinDailyAirTemp = state.dataGlobal->DayOfSim;
313 : }
314 :
315 365 : if (tdwd.dryBulbTemp > maxDailyAirTemp) {
316 16 : maxDailyAirTemp = tdwd.dryBulbTemp;
317 : }
318 :
319 : } // ... End day loop.
320 :
321 1 : annualAveAirTemp = annualAveAirTemp_num / state.dataWeather->NumDaysInYear; // Used for initializing domain
322 :
323 : // Reset Environment when done reading data
324 1 : --state.dataWeather->NumOfEnvrn; // May need better way of eliminating the extra environment that was added to read the data
325 1 : --state.dataWeather->TotRunPers;
326 1 : state.dataGlobal->KindOfSim = KindOfSim_reset;
327 1 : state.dataWeather->RPReadAllWeatherData = false;
328 1 : state.dataWeather->Environment.redimension(state.dataWeather->NumOfEnvrn);
329 1 : state.dataWeather->RunPeriodInput.redimension(state.dataWeather->TotRunPers);
330 1 : state.dataWeather->Envrn = Envrn_reset;
331 1 : state.dataGlobal->TimeStep = TimeStep_reset;
332 1 : state.dataGlobal->HourOfDay = HourOfDay_reset;
333 1 : state.dataGlobal->BeginEnvrnFlag = BeginEnvrnFlag_reset;
334 1 : state.dataGlobal->EndEnvrnFlag = EndEnvrnFlag_reset;
335 1 : state.dataEnvrn->EndMonthFlag = EndMonthFlag_reset;
336 1 : state.dataGlobal->WarmupFlag = WarmupFlag_reset;
337 1 : state.dataGlobal->DayOfSim = DayOfSim_reset;
338 1 : state.dataGlobal->DayOfSimChr = DayOfSimChr_reset;
339 1 : state.dataReportFlag->NumOfWarmupDays = NumOfWarmupDays_reset;
340 1 : state.dataGlobal->BeginDayFlag = BeginDayFlag_reset;
341 1 : state.dataGlobal->EndDayFlag = EndDayFlag_reset;
342 1 : state.dataGlobal->BeginHourFlag = BeginHourFlag_reset;
343 1 : state.dataGlobal->EndHourFlag = EndHourFlag_reset;
344 1 : }
345 :
346 : //******************************************************************************
347 :
348 1 : 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 1 : constexpr Real64 surfaceLayerThickness = 2.0;
359 1 : constexpr Real64 surfaceLayerCellThickness = 0.015;
360 1 : constexpr int surfaceLayerNumCells = static_cast<int>(surfaceLayerThickness / surfaceLayerCellThickness);
361 :
362 : // Center layer parameters
363 1 : constexpr int centerLayerNumCells = 80;
364 :
365 : // Deep layer parameters
366 1 : constexpr Real64 deepLayerThickness = 0.2;
367 1 : constexpr Real64 deepLayerCellThickness = surfaceLayerCellThickness;
368 1 : constexpr int deepLayerNumCells = static_cast<int>(deepLayerThickness / deepLayerCellThickness);
369 :
370 : // Other
371 1 : Real64 currentCellDepth = 0.0;
372 :
373 1 : totalNumCells = surfaceLayerNumCells + centerLayerNumCells + deepLayerNumCells;
374 :
375 : // Allocate arrays
376 1 : cellArray.allocate(totalNumCells);
377 1 : cellDepths.allocate(totalNumCells);
378 :
379 227 : for (int i = 1; i <= totalNumCells; ++i) {
380 :
381 : // Reference to thisCell
382 226 : auto &thisCell = cellArray(i);
383 :
384 : // Set the index
385 226 : thisCell.index = i;
386 :
387 : // Give thickness to the cells
388 226 : if (i <= surfaceLayerNumCells) {
389 : // Constant thickness mesh here
390 133 : thisCell.thickness = surfaceLayerCellThickness;
391 :
392 93 : } else if (i <= (centerLayerNumCells + surfaceLayerNumCells)) {
393 : // Geometric expansion/contraction here
394 80 : int numCenterCell = i - surfaceLayerNumCells;
395 :
396 80 : if (numCenterCell <= (centerLayerNumCells / 2)) {
397 40 : Real64 centerLayerExpansionCoeff = 1.10879;
398 40 : thisCell.thickness = surfaceLayerCellThickness * std::pow(centerLayerExpansionCoeff, numCenterCell);
399 : } else {
400 40 : thisCell.thickness =
401 40 : cellArray((surfaceLayerNumCells + (centerLayerNumCells / 2)) - (numCenterCell - (centerLayerNumCells / 2))).thickness;
402 : }
403 : } else {
404 : // Constant thickness mesh here
405 13 : thisCell.thickness = deepLayerCellThickness;
406 : }
407 :
408 : // Set minimum z value
409 226 : thisCell.minZValue = currentCellDepth;
410 :
411 : // Populate depth array for use later when looking up temperatures
412 226 : cellDepths(i) = currentCellDepth + thisCell.thickness / 2.0;
413 :
414 : // Update local counter
415 226 : currentCellDepth += thisCell.thickness;
416 :
417 : // Set maximum z value
418 226 : thisCell.maxZValue = currentCellDepth;
419 :
420 : // Set base properties
421 226 : thisCell.props.conductivity = baseConductivity;
422 226 : thisCell.props.density = baseDensity;
423 226 : thisCell.props.specificHeat = baseSpecificHeat;
424 226 : thisCell.props.diffusivity = baseConductivity / (baseDensity * baseSpecificHeat);
425 : }
426 1 : }
427 :
428 : //******************************************************************************
429 :
430 1 : 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 1 : timeStepInSeconds = Constant::rSecsInDay;
441 1 : bool convergedFinal = false;
442 :
443 1 : initDomain(state);
444 :
445 : // Loop until converged
446 : do {
447 :
448 : // loop over all days
449 4026 : for (state.dataGlobal->FDsimDay = 1; state.dataGlobal->FDsimDay <= state.dataWeather->NumDaysInYear; ++state.dataGlobal->FDsimDay) {
450 :
451 4015 : bool iterationConverged = false;
452 :
453 4015 : doStartOfTimeStepInits();
454 :
455 : // Loop until iteration temperature converges
456 : do {
457 :
458 : // For all cells
459 666326947 : for (int cell = 1; cell <= totalNumCells; ++cell) {
460 :
461 663391586 : if (cell == 1) {
462 2935361 : updateSurfaceCellTemperature(state);
463 660456225 : } else if (cell > 1 && cell < totalNumCells) {
464 657520864 : updateGeneralDomainCellTemperature(cell);
465 2935361 : } else if (cell == totalNumCells) {
466 2935361 : updateBottomCellTemperature();
467 : }
468 : }
469 :
470 : // Check iteration temperature convergence
471 2935361 : iterationConverged = checkIterationTemperatureConvergence();
472 :
473 2935361 : if (!iterationConverged) {
474 : // Shift temperatures for next iteration
475 2931346 : updateIterationTemperatures();
476 : }
477 :
478 2935361 : } while (!iterationConverged);
479 :
480 : // Shift temperatures for next timestep
481 4015 : updateTimeStepTemperatures(state);
482 : }
483 :
484 : // Check final temperature convergence
485 11 : convergedFinal = checkFinalTemperatureConvergence(state);
486 :
487 11 : } while (!convergedFinal);
488 1 : }
489 :
490 : //******************************************************************************
491 :
492 2935361 : 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 2935361 : Real64 numerator = 0.0;
503 2935361 : Real64 denominator = 0.0;
504 2935361 : Real64 resistance = 0.0;
505 : Real64 G_hr;
506 : Real64 Cd;
507 :
508 2935361 : Real64 constexpr rho_water = 998.0; // [kg/m3]
509 : // evapotranspiration parameters
510 2935361 : Real64 constexpr absor_Corrected = 0.77;
511 2935361 : constexpr Real64 convert_Wm2_To_MJhrmin = 3600.0 / 1000000.0;
512 2935361 : constexpr Real64 convert_MJhrmin_To_Wm2 = 1.0 / convert_Wm2_To_MJhrmin;
513 :
514 2935361 : const auto &thisCell = cellArray(1);
515 2935361 : const auto &cellBelow_thisCell = cellArray(2);
516 2935361 : const auto &cwd = weatherDataArray(state.dataGlobal->FDsimDay); // "Current Weather Day"
517 :
518 : // Add effect from previous time step
519 2935361 : numerator += thisCell.temperature_prevTimeStep;
520 2935361 : ++denominator;
521 :
522 : // Conduction to lower cell
523 2935361 : resistance = thisCell.thickness / 2.0 / (thisCell.props.conductivity * thisCell.conductionArea) +
524 2935361 : cellBelow_thisCell.thickness / 2.0 / (cellBelow_thisCell.props.conductivity * cellBelow_thisCell.conductionArea);
525 2935361 : numerator += thisCell.beta / resistance * cellBelow_thisCell.temperature;
526 2935361 : denominator += thisCell.beta / resistance;
527 :
528 : // Convection to atmosphere
529 2935361 : if (cwd.windSpeed > 0.1) {
530 2935361 : Real64 constexpr airSpecificHeat = 1003; // '[J/kg-K]
531 2935361 : resistance = 208.0 / (cwd.airDensity * airSpecificHeat * cwd.windSpeed * thisCell.conductionArea);
532 : } else {
533 : // Future development should include additional natural convection effects here
534 : }
535 2935361 : numerator += thisCell.beta / resistance * cwd.dryBulbTemp;
536 2935361 : denominator += thisCell.beta / resistance;
537 :
538 : // For convenience convert to Kelvin once
539 2935361 : 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 2935361 : const Real64 incidentSolar_MJhrmin = cwd.horizontalRadiation * convert_Wm2_To_MJhrmin;
544 :
545 : // Absorbed solar radiation, [MJ/hr-min]
546 2935361 : const Real64 absorbedIncidentSolar_MJhrmin = absor_Corrected * incidentSolar_MJhrmin;
547 :
548 : // Calculate saturated vapor pressure, [kPa]
549 2935361 : const Real64 vaporPressureSaturated_kPa = 0.6108 * std::exp(17.27 * cwd.dryBulbTemp / (cwd.dryBulbTemp + 237.3));
550 :
551 : // Calculate actual vapor pressure, [kPa]
552 2935361 : const Real64 vaporPressureActual_kPa = vaporPressureSaturated_kPa * cwd.relativeHumidity;
553 :
554 : // Calculate another Q term, [MJ/m2-hr]
555 2935361 : 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 2935361 : const Real64 netIncidentRadiation_MJhr = absorbedIncidentSolar_MJhrmin - QRAD_NL;
559 :
560 : // constant
561 2935361 : constexpr Real64 CN = 37.0;
562 :
563 : // Check whether there was sun
564 2935361 : if (netIncidentRadiation_MJhr < 0.0) {
565 701358 : G_hr = 0.5 * netIncidentRadiation_MJhr;
566 701358 : Cd = 0.96;
567 : } else {
568 2234003 : G_hr = 0.1 * netIncidentRadiation_MJhr;
569 2234003 : Cd = 0.24;
570 : }
571 :
572 2935361 : const Real64 slope_S = 2503.0 * std::exp(17.27 * cwd.dryBulbTemp / (cwd.dryBulbTemp + 237.3)) / pow_2(cwd.dryBulbTemp + 237.3);
573 2935361 : constexpr Real64 pressure = 98.0;
574 2935361 : constexpr Real64 psychrometricConstant = 0.665e-3 * pressure;
575 :
576 : // Evapotranspiration constant, [mm/hr]
577 2935361 : const Real64 evapotransFluidLoss_mmhr =
578 2935361 : (evapotransCoeff * slope_S * (netIncidentRadiation_MJhr - G_hr) +
579 2935361 : psychrometricConstant * (CN / currAirTempK) * cwd.windSpeed * (vaporPressureSaturated_kPa - vaporPressureActual_kPa)) /
580 2935361 : (slope_S + psychrometricConstant * (1 + Cd * cwd.windSpeed));
581 :
582 : // Convert units, [m/hr]
583 2935361 : 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 2935361 : const Real64 latentHeatVaporization = 2.501 - 2.361e-3 * thisCell.temperature_prevTimeStep;
591 :
592 : // Calculate evapotranspiration heat loss, [MJ/m2-hr]
593 2935361 : const Real64 evapotransHeatLoss_MJhrmin = rho_water * evapotransFluidLoss_mhr * latentHeatVaporization;
594 :
595 : // Convert net incident solar units, [W/m2]
596 2935361 : const Real64 netIncidentRadiation_Wm2 = netIncidentRadiation_MJhr * convert_MJhrmin_To_Wm2;
597 :
598 : // Convert evapotranspiration units, [W/m2]
599 2935361 : const Real64 evapotransHeatLoss_Wm2 = evapotransHeatLoss_MJhrmin * convert_MJhrmin_To_Wm2;
600 :
601 : // Calculate overall net heat ?gain? into the cell, [W]
602 2935361 : const Real64 incidentHeatGain = (netIncidentRadiation_Wm2 - evapotransHeatLoss_Wm2) * thisCell.conductionArea;
603 :
604 : // Add any solar/evapotranspiration heat gain here
605 2935361 : numerator += thisCell.beta * incidentHeatGain;
606 :
607 : // Calculate the return temperature and leave
608 2935361 : cellArray(1).temperature = numerator / denominator;
609 2935361 : }
610 :
611 : //******************************************************************************
612 :
613 657520864 : 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 657520864 : Real64 numerator = 0.0;
624 657520864 : Real64 denominator = 0.0;
625 657520864 : Real64 resistance = 0.0;
626 :
627 657520864 : auto &thisCell = cellArray(cell);
628 657520864 : const auto &cellAbove_thisCell = cellArray(cell - 1);
629 657520864 : const auto &cellBelow_thisCell = cellArray(cell + 1);
630 :
631 : // add effect from cell history
632 657520864 : numerator += thisCell.temperature_prevTimeStep;
633 657520864 : ++denominator;
634 :
635 : // Conduction resistance between this cell and above cell
636 657520864 : resistance = thisCell.thickness / 2.0 / (thisCell.conductionArea * thisCell.props.conductivity) +
637 657520864 : cellAbove_thisCell.thickness / 2.0 / (cellAbove_thisCell.conductionArea * cellAbove_thisCell.props.conductivity);
638 :
639 657520864 : numerator += thisCell.beta / resistance * cellAbove_thisCell.temperature;
640 657520864 : denominator += thisCell.beta / resistance;
641 :
642 : // Conduction resistance between this cell and below cell
643 657520864 : resistance = thisCell.thickness / 2.0 / (thisCell.conductionArea * thisCell.props.conductivity) +
644 657520864 : cellBelow_thisCell.thickness / 2.0 / (cellBelow_thisCell.conductionArea * cellBelow_thisCell.props.conductivity);
645 :
646 657520864 : numerator += thisCell.beta / resistance * cellBelow_thisCell.temperature;
647 657520864 : denominator += thisCell.beta / resistance;
648 :
649 : // now that we have passed all directions, update the temperature
650 657520864 : thisCell.temperature = numerator / denominator;
651 657520864 : }
652 :
653 : //******************************************************************************
654 :
655 2935361 : 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 2935361 : Real64 numerator = 0.0;
671 2935361 : Real64 denominator = 0.0;
672 2935361 : Real64 resistance = 0.0;
673 2935361 : Real64 constexpr geothermalGradient = 0.025; // C/m
674 :
675 2935361 : auto &thisCell = cellArray(totalNumCells);
676 2935361 : auto &cellAbove_thisCell = cellArray(totalNumCells - 1);
677 :
678 : // Initialize
679 2935361 : numerator += thisCell.temperature_prevTimeStep;
680 2935361 : ++denominator;
681 :
682 : // Conduction resistance between this cell and above cell
683 2935361 : resistance = ((thisCell.thickness / 2.0) / (thisCell.conductionArea * thisCell.props.conductivity)) +
684 2935361 : ((cellAbove_thisCell.thickness / 2.0) / (cellAbove_thisCell.conductionArea * cellAbove_thisCell.props.conductivity));
685 :
686 2935361 : numerator += (thisCell.beta / resistance) * cellAbove_thisCell.temperature;
687 2935361 : denominator += thisCell.beta / resistance;
688 :
689 : // Geothermal gradient heat transfer
690 2935361 : Real64 HTBottom = geothermalGradient * thisCell.props.conductivity * thisCell.conductionArea;
691 :
692 2935361 : numerator += thisCell.beta * HTBottom;
693 :
694 2935361 : cellArray(totalNumCells).temperature = numerator / denominator;
695 2935361 : }
696 :
697 : //******************************************************************************
698 :
699 11 : 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 11 : bool converged = true;
710 11 : Real64 constexpr finalTempConvergenceCriteria = 0.05;
711 :
712 11 : if (state.dataGlobal->FDnumIterYears == maxYearsToIterate) {
713 1 : return converged;
714 : }
715 :
716 2270 : for (int cell = 1; cell <= totalNumCells; ++cell) {
717 :
718 2260 : auto &thisCell = cellArray(cell);
719 :
720 2260 : if (std::abs(thisCell.temperature - thisCell.temperature_finalConvergence) >= finalTempConvergenceCriteria) {
721 687 : converged = false;
722 : }
723 :
724 2260 : thisCell.temperature_finalConvergence = thisCell.temperature;
725 : }
726 :
727 10 : ++state.dataGlobal->FDnumIterYears;
728 :
729 10 : return converged;
730 : }
731 :
732 : //******************************************************************************
733 :
734 2935361 : bool FiniteDiffGroundTempsModel::checkIterationTemperatureConvergence()
735 : {
736 : // SUBROUTINE INFORMATION:
737 : // AUTHOR Matt Mitchell
738 : // DATE WRITTEN Summer 2015
739 :
740 : // PURPOSE OF THIS SUBROUTINE:
741 : // Checks iteration temperature convergence
742 :
743 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
744 2935361 : bool converged = true;
745 2935361 : Real64 constexpr iterationTempConvergenceCriteria = 0.00001;
746 :
747 16530796 : for (int cell = 1; cell <= totalNumCells; ++cell) {
748 :
749 16526781 : if (std::abs(cellArray(cell).temperature - cellArray(cell).temperature_prevIteration) >= iterationTempConvergenceCriteria) {
750 2931346 : converged = false;
751 2931346 : break;
752 : }
753 : }
754 :
755 2935361 : return converged;
756 : }
757 :
758 : //******************************************************************************
759 :
760 1 : void FiniteDiffGroundTempsModel::initDomain(EnergyPlusData &state)
761 : {
762 : // SUBROUTINE INFORMATION:
763 : // AUTHOR Matt Mitchell
764 : // DATE WRITTEN Summer 2015
765 :
766 : // PURPOSE OF THIS SUBROUTINE:
767 : // Initializes model using Kusuda-Achenbach model.
768 : // Average ground temp initialized to average annual air temperature
769 :
770 : // Temporary KA model for initialization
771 1 : auto tempModel = std::make_unique<KusudaGroundTempsModel>(); // (AUTO_OK) Why does this have to be a unique_ptr?
772 :
773 1 : tempModel->Name = "KAModelForFDModel";
774 1 : tempModel->modelType = GroundTemp::ModelType::Kusuda;
775 1 : tempModel->aveGroundTemp = annualAveAirTemp;
776 2 : tempModel->aveGroundTempAmplitude =
777 1 : (maxDailyAirTemp - minDailyAirTemp) / 4.0; // Rough estimate here. Ground temps will not swing as far as the air temp.
778 1 : tempModel->phaseShiftInSecs = dayOfMinDailyAirTemp * Constant::rSecsInDay;
779 1 : tempModel->groundThermalDiffusivity = baseConductivity / (baseDensity * baseSpecificHeat);
780 :
781 : // Initialize temperatures and volume
782 227 : for (int cell = 1; cell <= totalNumCells; ++cell) {
783 226 : auto &thisCell = cellArray(cell);
784 :
785 226 : Real64 depth = (thisCell.maxZValue + thisCell.minZValue) / 2.0;
786 :
787 : // Initialize temperatures
788 226 : if (tempModel) {
789 226 : thisCell.temperature = tempModel->getGroundTempAtTimeInSeconds(state, depth, 0.0); // Initialized at first day of year
790 : }
791 226 : thisCell.temperature_finalConvergence = thisCell.temperature;
792 226 : thisCell.temperature_prevIteration = thisCell.temperature;
793 226 : thisCell.temperature_prevTimeStep = thisCell.temperature;
794 :
795 : // Set cell volume
796 226 : thisCell.volume = thisCell.thickness * thisCell.conductionArea;
797 : }
798 :
799 : // Initialize freezing calculation variables
800 1 : evaluateSoilRhoCpInit();
801 :
802 : // Initialize the groundTemps array
803 1 : groundTemps.dimension({1, state.dataWeather->NumDaysInYear}, {1, totalNumCells}, 0.0);
804 :
805 1 : tempModel.reset();
806 1 : }
807 :
808 : //******************************************************************************
809 :
810 2931346 : void FiniteDiffGroundTempsModel::updateIterationTemperatures()
811 : {
812 : // SUBROUTINE INFORMATION:
813 : // AUTHOR Matt Mitchell
814 : // DATE WRITTEN Summer 2015
815 :
816 : // PURPOSE OF THIS SUBROUTINE:
817 : // Updates iteration temperatures for convergence checks
818 :
819 665415542 : for (int cell = 1; cell <= totalNumCells; ++cell) {
820 662484196 : cellArray(cell).temperature_prevIteration = cellArray(cell).temperature;
821 : }
822 2931346 : }
823 :
824 : //******************************************************************************
825 :
826 4015 : void FiniteDiffGroundTempsModel::updateTimeStepTemperatures(const EnergyPlusData &state)
827 : {
828 : // SUBROUTINE INFORMATION:
829 : // AUTHOR Matt Mitchell
830 : // DATE WRITTEN Summer 2015
831 :
832 : // PURPOSE OF THIS SUBROUTINE:
833 : // Updates timestep temperatures for convergence checks.
834 :
835 911405 : for (int cell = 1; cell <= totalNumCells; ++cell) {
836 :
837 907390 : auto &thisCell = cellArray(cell);
838 :
839 907390 : thisCell.temperature_prevTimeStep = thisCell.temperature;
840 :
841 : // Log temps for later use
842 907390 : groundTemps(state.dataGlobal->FDsimDay, cell) = thisCell.temperature;
843 : }
844 4015 : }
845 :
846 : //******************************************************************************
847 :
848 4015 : void FiniteDiffGroundTempsModel::doStartOfTimeStepInits()
849 : {
850 : // SUBROUTINE INFORMATION:
851 : // AUTHOR Matt Mitchell
852 : // DATE WRITTEN Summer 2015
853 :
854 : // PURPOSE OF THIS SUBROUTINE:
855 : // Updates cell properties for each timestep
856 :
857 911405 : for (int cell = 1; cell <= totalNumCells; ++cell) {
858 :
859 907390 : auto &thisCell = cellArray(cell);
860 :
861 907390 : evaluateSoilRhoCpCell(cell);
862 :
863 907390 : thisCell.beta = (timeStepInSeconds / (thisCell.props.rhoCp * thisCell.volume));
864 : }
865 4015 : }
866 :
867 : //******************************************************************************
868 :
869 2790720 : Real64 FiniteDiffGroundTempsModel::interpolate(Real64 const x, Real64 const x_hi, Real64 const x_low, Real64 const y_hi, Real64 const y_low)
870 : {
871 2790720 : return (x - x_low) / (x_hi - x_low) * (y_hi - y_low) + y_low;
872 : }
873 :
874 : //******************************************************************************
875 :
876 930240 : Real64 FiniteDiffGroundTempsModel::getGroundTemp(EnergyPlusData &state)
877 : {
878 :
879 : // SUBROUTINE INFORMATION:
880 : // AUTHOR Matt Mitchell
881 : // DATE WRITTEN Summer 2015
882 :
883 : // PURPOSE OF THIS SUBROUTINE:
884 : // Interpolates between days and depths to find correct ground temperature
885 :
886 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
887 : // Interpolation variables
888 : int i0; // First day
889 : int i1; // Next day
890 : int j1; // Next cell index (with depth greater than y-depth
891 : // Real64 T_i0_j0; // Temp at int( x-day ); cell lower_bound( y-depth )
892 : // Real64 T_i1_j0; // Temp at int( x-day ) + 1; cell lower_bound( y-depth )
893 : // Real64 T_i0_j1; // Temp at int( x-day ); cell lower_bound( y-depth ) + 1
894 : // Real64 T_i1_j1; // Temp at int( x-day ) + 1; cell lower_bound( y-depth ) + 1
895 : // Real64 T_ix_j0; // Temp at x-day; cell lower_bound( y-depth )
896 : // Real64 T_ix_j1; // Temp at x-day; cell lower_bound( y-depth ) + 1
897 : // RETURNS: Final Temperature--Temp at x-day; y-depth
898 :
899 930240 : if (depth < 0.0) {
900 0 : depth = 0.0;
901 : }
902 :
903 : // Get index of nearest cell with depth less than depth
904 930240 : auto it = std::lower_bound(cellDepths.begin(), cellDepths.end(), depth);
905 930240 : int j0 = static_cast<int>(std::distance(cellDepths.begin(), it)); // Cell index with depth less than y-depth
906 :
907 : // Compensate for 1-based array
908 930240 : ++j0;
909 :
910 : // Fraction of day
911 930240 : const Real64 dayFrac = simTimeInDays - static_cast<int>(simTimeInDays); // Fraction of day
912 :
913 930240 : if (j0 < totalNumCells - 1) {
914 : // All depths within domain
915 930240 : j1 = j0 + 1;
916 :
917 930240 : if (simTimeInDays <= 1 || simTimeInDays >= state.dataWeather->NumDaysInYear) {
918 : // First day of year, last day of year, and leap day
919 : // Interpolate between first and last day
920 930240 : i0 = state.dataWeather->NumDaysInYear;
921 930240 : i1 = 1;
922 :
923 : // Lookup ground temps
924 930240 : const Real64 T_i0_j0 = groundTemps(i0, j0);
925 930240 : const Real64 T_i0_j1 = groundTemps(i0, j1);
926 930240 : const Real64 T_i1_j0 = groundTemps(i1, j0);
927 930240 : const Real64 T_i1_j1 = groundTemps(i1, j1);
928 :
929 : // Interpolate between days holding depth constant
930 930240 : const Real64 T_ix_j0 = interpolate(dayFrac, 1, 0, T_i1_j0, T_i0_j0);
931 930240 : const Real64 T_ix_j1 = interpolate(dayFrac, 1, 0, T_i1_j1, T_i0_j1);
932 :
933 : // Interpolate to correct depth now that we're at the right time
934 930240 : return interpolate(depth, cellDepths(j1), cellDepths(j0), T_ix_j1, T_ix_j0);
935 : }
936 :
937 : // All other days
938 0 : i0 = static_cast<int>(simTimeInDays);
939 0 : i1 = i0 + 1;
940 :
941 : // Lookup ground temps
942 0 : const Real64 T_i0_j0 = groundTemps(i0, j0);
943 0 : const Real64 T_i0_j1 = groundTemps(i0, j1);
944 0 : const Real64 T_i1_j0 = groundTemps(i1, j0);
945 0 : const Real64 T_i1_j1 = groundTemps(i1, j1);
946 :
947 : // Interpolate between days holding depth constant
948 0 : const Real64 T_ix_j0 = interpolate(dayFrac, 1, 0, T_i1_j0, T_i0_j0);
949 0 : const Real64 T_ix_j1 = interpolate(dayFrac, 1, 0, T_i1_j1, T_i0_j1);
950 :
951 : // Interpolate to correct depth now that we're at the right time
952 0 : return interpolate(depth, cellDepths(j1), cellDepths(j0), T_ix_j1, T_ix_j0);
953 : }
954 :
955 : // Requesting a temperature deeper than domain. Pass deepest point in domain.
956 0 : j0 = totalNumCells;
957 0 : j1 = j0;
958 :
959 0 : if (simTimeInDays <= 1 || simTimeInDays >= state.dataWeather->NumDaysInYear) {
960 : // First day of year, last day of year, and leap day
961 : // Interpolate between first and last day
962 0 : i0 = state.dataWeather->NumDaysInYear;
963 0 : i1 = 1;
964 :
965 : // Lookup ground temps
966 0 : const Real64 T_i0_j1 = groundTemps(i0, j1);
967 0 : const Real64 T_i1_j1 = groundTemps(i1, j1);
968 :
969 : // Interpolate between days holding depth constant
970 0 : return interpolate(dayFrac, 1, 0, T_i1_j1, T_i0_j1);
971 : }
972 :
973 : // All other days
974 0 : i0 = static_cast<int>(simTimeInDays);
975 0 : i1 = i0 + 1;
976 :
977 : // Lookup ground temps
978 0 : const Real64 T_i0_j1 = groundTemps(i0, j1);
979 0 : const Real64 T_i1_j1 = groundTemps(i1, j1);
980 :
981 : // Interpolate between days holding depth constant
982 0 : return interpolate(dayFrac, 1, 0, T_i1_j1, T_i0_j1);
983 : }
984 :
985 : //******************************************************************************
986 :
987 930240 : Real64 FiniteDiffGroundTempsModel::getGroundTempAtTimeInSeconds(EnergyPlusData &state, Real64 const _depth, Real64 const seconds)
988 : {
989 : // SUBROUTINE INFORMATION:
990 : // AUTHOR Matt Mitchell
991 : // DATE WRITTEN Summer 2015
992 :
993 : // PURPOSE OF THIS SUBROUTINE:
994 : // Retrieves ground temperature when input time is in seconds
995 :
996 930240 : depth = _depth;
997 :
998 930240 : simTimeInDays = seconds / Constant::rSecsInDay;
999 :
1000 930240 : if (simTimeInDays > state.dataWeather->NumDaysInYear) {
1001 0 : simTimeInDays = remainder(simTimeInDays, state.dataWeather->NumDaysInYear);
1002 : }
1003 :
1004 930240 : return getGroundTemp(state);
1005 : }
1006 :
1007 : //******************************************************************************
1008 :
1009 0 : Real64 FiniteDiffGroundTempsModel::getGroundTempAtTimeInMonths(EnergyPlusData &state, Real64 const _depth, int const month)
1010 : {
1011 : // SUBROUTINE INFORMATION:
1012 : // AUTHOR Matt Mitchell
1013 : // DATE WRITTEN Summer 2015
1014 :
1015 : // PURPOSE OF THIS SUBROUTINE:
1016 : // Returns ground temperature when input time is in months
1017 :
1018 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
1019 : // TODO: Fixing this to be floating point 12.0 causes diffs and failed tests
1020 0 : Real64 const aveDaysInMonth = state.dataWeather->NumDaysInYear / 12;
1021 :
1022 0 : depth = _depth;
1023 :
1024 : // Convert months to days. Puts time in middle of specified month
1025 0 : simTimeInDays = aveDaysInMonth * ((month - 1) + 0.5);
1026 :
1027 0 : if (simTimeInDays > state.dataWeather->NumDaysInYear) {
1028 0 : simTimeInDays = remainder(simTimeInDays, state.dataWeather->NumDaysInYear);
1029 : }
1030 :
1031 : // Get and return ground temperature
1032 0 : return getGroundTemp(state);
1033 : }
1034 :
1035 : //******************************************************************************
1036 :
1037 907390 : void FiniteDiffGroundTempsModel::evaluateSoilRhoCpCell(int const cell)
1038 : {
1039 : // SUBROUTINE INFORMATION:
1040 : // AUTHOR Edwin Lee
1041 : // DATE WRITTEN Summer 2011
1042 :
1043 : // PURPOSE OF THIS SUBROUTINE:
1044 : // Evaluates the soil properties on a single cell
1045 :
1046 : // Real64 rhoCP_soil;
1047 :
1048 : // These vary by domain now, so we must be careful to retrieve them every time
1049 :
1050 907390 : auto &thisCell = cellArray(cell);
1051 :
1052 : // set some temperatures here for generalization -- these could be set in the input file
1053 : // constexpr Real64 frzAllIce = -0.5;
1054 : // constexpr Real64 frzIceTrans = -0.4;
1055 : // constexpr Real64 frzLiqTrans = -0.1;
1056 : // constexpr Real64 frzAllLiq = 0.0;
1057 :
1058 : // calculate this cell's new Cp value based on the cell temperature
1059 : // if (thisCell.temperature >= frzAllLiq) {
1060 : // rhoCP_soil = rhoCp_soil_liq_1;
1061 : // } else if (thisCell.temperature <= frzAllIce) {
1062 : // rhoCP_soil = rhoCP_soil_ice;
1063 : // } else if (thisCell.temperature > frzLiqTrans) {
1064 : // rhoCP_soil = rhoCp_soil_liq_1 + (rhoCP_soil_transient - rhoCP_soil_liq) / (frzAllLiq - frzLiqTrans) * (frzAllLiq -
1065 : // thisCell.temperature);
1066 : // } else if (thisCell.temperature >= frzIceTrans) {
1067 : // rhoCP_soil = rhoCP_soil_transient;
1068 : // } else {
1069 : // rhoCP_soil = rhoCP_soil_ice + (rhoCP_soil_transient - rhoCP_soil_ice) / (frzIceTrans - frzAllIce) * (thisCell.temperature - frzAllIce);
1070 : // }
1071 :
1072 : // TODO: The calculated rhoCP_soil is commented on this line and never used. Curious.
1073 907390 : thisCell.props.rhoCp = baseDensity * baseSpecificHeat; // rhoCP_soil;
1074 :
1075 907390 : thisCell.props.specificHeat = thisCell.props.rhoCp / thisCell.props.density;
1076 907390 : }
1077 :
1078 1 : void FiniteDiffGroundTempsModel::evaluateSoilRhoCpInit()
1079 : {
1080 : // SUBROUTINE INFORMATION:
1081 : // AUTHOR Edwin Lee
1082 : // DATE WRITTEN Summer 2011
1083 :
1084 : // PURPOSE OF THIS SUBROUTINE:
1085 : // Evaluates the soil properties
1086 :
1087 : // These vary by domain now, so we must be careful to retrieve them every time
1088 1 : const Real64 Theta_liq = waterContent;
1089 1 : const Real64 Theta_sat = saturatedWaterContent;
1090 :
1091 : // Assumption
1092 1 : const Real64 Theta_ice = Theta_liq;
1093 :
1094 : //'Cp (freezing) calculations
1095 1 : constexpr Real64 rho_ice = 917.0; //'Kg / m3
1096 1 : constexpr Real64 rho_liq = 1000.0; //'kg / m3
1097 1 : rhoCp_soil_liq_1 = 1225000.0 / (1.0 - Theta_sat); // J/m3K
1098 : // from( " An improved model for predicting soil thermal conductivity from water content at room temperature, Fig 4" )
1099 1 : constexpr Real64 CP_liq = 4180.0; //'J / KgK
1100 1 : constexpr Real64 CP_ice = 2066.0; //'J / KgK
1101 1 : constexpr Real64 Lat_fus = 334000.0; //'J / Kg
1102 1 : constexpr Real64 Cp_transient = Lat_fus / 0.4 + (0.5 * CP_ice - (CP_liq + CP_ice) / 2.0 * 0.1) / 0.4;
1103 : // from( " Numerical and experimental investigation of melting and freezing processes in phase change material storage" )
1104 1 : rhoCP_soil_liq = rhoCp_soil_liq_1 * (1.0 - Theta_sat) + rho_liq * CP_liq * Theta_liq;
1105 1 : rhoCP_soil_transient = rhoCp_soil_liq_1 * (1.0 - Theta_sat) + ((rho_liq + rho_ice) / 2.0) * Cp_transient * Theta_ice;
1106 1 : rhoCP_soil_ice = rhoCp_soil_liq_1 * (1.0 - Theta_sat) + rho_ice * CP_ice * Theta_ice; //'!J / m3K
1107 1 : }
1108 :
1109 : //******************************************************************************
1110 :
1111 : } // namespace GroundTemp
1112 : } // namespace EnergyPlus
|
