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