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