Line data Source code
1 : // EnergyPlus, Copyright (c) 1996-2025, The Board of Trustees of the University of Illinois,
2 : // The Regents of the University of California, through Lawrence Berkeley National Laboratory
3 : // (subject to receipt of any required approvals from the U.S. Dept. of Energy), Oak Ridge
4 : // National Laboratory, managed by UT-Battelle, Alliance for Sustainable Energy, LLC, and other
5 : // contributors. All rights reserved.
6 : //
7 : // NOTICE: This Software was developed under funding from the U.S. Department of Energy and the
8 : // U.S. Government consequently retains certain rights. As such, the U.S. Government has been
9 : // granted for itself and others acting on its behalf a paid-up, nonexclusive, irrevocable,
10 : // worldwide license in the Software to reproduce, distribute copies to the public, prepare
11 : // derivative works, and perform publicly and display publicly, and to permit others to do so.
12 : //
13 : // Redistribution and use in source and binary forms, with or without modification, are permitted
14 : // provided that the following conditions are met:
15 : //
16 : // (1) Redistributions of source code must retain the above copyright notice, this list of
17 : // conditions and the following disclaimer.
18 : //
19 : // (2) Redistributions in binary form must reproduce the above copyright notice, this list of
20 : // conditions and the following disclaimer in the documentation and/or other materials
21 : // provided with the distribution.
22 : //
23 : // (3) Neither the name of the University of California, Lawrence Berkeley National Laboratory,
24 : // the University of Illinois, U.S. Dept. of Energy nor the names of its contributors may be
25 : // used to endorse or promote products derived from this software without specific prior
26 : // written permission.
27 : //
28 : // (4) Use of EnergyPlus(TM) Name. If Licensee (i) distributes the software in stand-alone form
29 : // without changes from the version obtained under this License, or (ii) Licensee makes a
30 : // reference solely to the software portion of its product, Licensee must refer to the
31 : // software as "EnergyPlus version X" software, where "X" is the version number Licensee
32 : // obtained under this License and may not use a different name for the software. Except as
33 : // specifically required in this Section (4), Licensee shall not use in a company name, a
34 : // product name, in advertising, publicity, or other promotional activities any name, trade
35 : // name, trademark, logo, or other designation of "EnergyPlus", "E+", "e+" or confusingly
36 : // similar designation, without the U.S. Department of Energy's prior written consent.
37 : //
38 : // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
39 : // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY
40 : // AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
41 : // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
42 : // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
43 : // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
44 : // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
45 : // OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
46 : // POSSIBILITY OF SUCH DAMAGE.
47 :
48 : // C++ Headers
49 : #include <cassert>
50 : #include <cmath>
51 : #include <string>
52 :
53 : // ObjexxFCL Headers
54 : #include <ObjexxFCL/Array.functions.hh>
55 : #include <ObjexxFCL/Fmath.hh>
56 :
57 : // EnergyPlus Headers
58 : #include <AirflowNetwork/Solver.hpp>
59 : #include <EnergyPlus/Construction.hh>
60 : #include <EnergyPlus/Data/EnergyPlusData.hh>
61 : #include <EnergyPlus/DataEnvironment.hh>
62 : #include <EnergyPlus/DataHeatBalFanSys.hh>
63 : #include <EnergyPlus/DataHeatBalSurface.hh>
64 : #include <EnergyPlus/DataHeatBalance.hh>
65 : #include <EnergyPlus/DataIPShortCuts.hh>
66 : #include <EnergyPlus/DataMoistureBalance.hh>
67 : #include <EnergyPlus/DataSurfaces.hh>
68 : #include <EnergyPlus/EMSManager.hh>
69 : #include <EnergyPlus/General.hh>
70 : #include <EnergyPlus/HeatBalFiniteDiffManager.hh>
71 : #include <EnergyPlus/InputProcessing/InputProcessor.hh>
72 : #include <EnergyPlus/Material.hh>
73 : #include <EnergyPlus/MoistureBalanceEMPDManager.hh>
74 : #include <EnergyPlus/OutputProcessor.hh>
75 : #include <EnergyPlus/PhaseChangeModeling/HysteresisModel.hh>
76 : #include <EnergyPlus/PluginManager.hh>
77 : #include <EnergyPlus/UtilityRoutines.hh>
78 : #include <EnergyPlus/ZoneTempPredictorCorrector.hh>
79 :
80 : namespace EnergyPlus {
81 :
82 : namespace HeatBalFiniteDiffManager {
83 :
84 : // Module containing the heat balance simulation routines
85 :
86 : // MODULE INFORMATION:
87 : // AUTHOR Richard J. Liesen
88 : // DATE WRITTEN October 2003
89 : // RE-ENGINEERED Curtis Pedersen, 2006, Changed to Implicit FD calc for conduction.
90 : // and included enthalpy formulations for phase change materials
91 : // PURPOSE OF THIS MODULE:
92 : // To encapsulate the data and algorithms required to
93 : // manage the finite difference heat balance simulation on the building.
94 :
95 : // REFERENCES:
96 : // The MFD moisture balance method
97 : // C. O. Pedersen, Enthalpy Formulation of conduction heat transfer problems
98 : // involving latent heat, Simulation, Vol 18, No. 2, February 1972
99 : // Fan system Source/Sink heat value, and source/sink location temp from CondFD
100 :
101 10240031 : void ManageHeatBalFiniteDiff(EnergyPlusData &state,
102 : int const SurfNum,
103 : Real64 &SurfTempInTmp, // INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
104 : Real64 &TempSurfOutTmp // Outside Surface Temperature of each Heat Transfer Surface
105 : )
106 : {
107 : // SUBROUTINE INFORMATION:
108 : // AUTHOR Richard Liesen
109 : // DATE WRITTEN May 2000
110 :
111 : // PURPOSE OF THIS SUBROUTINE:
112 : // This subroutine manages the moisture balance method. It is called
113 : // from the HeatBalanceManager at the time step level.
114 : // This driver manages the calls to all of
115 : // the other drivers and simulation algorithms.
116 :
117 : // Get the moisture balance input at the beginning of the simulation only
118 10240031 : if (state.dataHeatBalFiniteDiffMgr->GetHBFiniteDiffInputFlag) {
119 : // Obtains conduction FD related parameters from input file
120 0 : GetCondFDInput(state);
121 0 : state.dataHeatBalFiniteDiffMgr->GetHBFiniteDiffInputFlag = false;
122 : }
123 : // Solve the zone heat & moisture balance using a finite difference solution
124 10240031 : CalcHeatBalFiniteDiff(state, SurfNum, SurfTempInTmp, TempSurfOutTmp);
125 10240031 : }
126 :
127 12 : void GetCondFDInput(EnergyPlusData &state)
128 : {
129 : // SUBROUTINE INFORMATION:
130 : // AUTHOR Curtis Pedersen
131 : // DATE WRITTEN July 2006
132 : // MODIFIED Brent Griffith Mar 2011, user settings
133 :
134 : // PURPOSE OF THIS SUBROUTINE:
135 : // This subroutine is the main driver for initializations for the variable property CondFD part of the
136 : // MFD algorithm
137 :
138 : static constexpr std::string_view routineName = "GetCondFDInput";
139 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
140 : int IOStat; // IO Status when calling get input subroutine
141 12 : Array1D_string MaterialNames(3); // Number of Material Alpha names defined
142 12 : Array1D_string ConstructionName(3); // Name of Construction with CondFDsimplified
143 : int MaterialNumAlpha; // Number of material alpha names being passed
144 : int MaterialNumProp; // Number of material properties being passed
145 12 : Array1D<Real64> MaterialProps; // Temporary array to transfer material properties (allocated based on user input)
146 12 : bool ErrorsFound(false); // If errors detected in input
147 : int propNum;
148 : int pcMat;
149 : int vcMat;
150 : int inegptr;
151 : bool nonInc;
152 :
153 12 : auto &s_ip = state.dataInputProcessing->inputProcessor;
154 12 : auto &s_ipsc = state.dataIPShortCut;
155 12 : auto &s_hbfd = state.dataHeatBalFiniteDiffMgr;
156 12 : auto &s_mat = state.dataMaterial;
157 :
158 : // user settings for numerical parameters
159 12 : s_ipsc->cCurrentModuleObject = "HeatBalanceSettings:ConductionFiniteDifference";
160 :
161 12 : if (s_ip->getNumObjectsFound(state, s_ipsc->cCurrentModuleObject) > 0) {
162 : int NumAlphas;
163 : int NumNumbers;
164 2 : s_ip->getObjectItem(state,
165 1 : s_ipsc->cCurrentModuleObject,
166 : 1,
167 1 : s_ipsc->cAlphaArgs,
168 : NumAlphas,
169 1 : s_ipsc->rNumericArgs,
170 : NumNumbers,
171 : IOStat,
172 1 : s_ipsc->lNumericFieldBlanks,
173 1 : s_ipsc->lAlphaFieldBlanks,
174 1 : s_ipsc->cAlphaFieldNames,
175 1 : s_ipsc->cNumericFieldNames);
176 :
177 1 : if (!s_ipsc->lAlphaFieldBlanks(1)) {
178 : {
179 2 : s_hbfd->CondFDSchemeType =
180 1 : static_cast<CondFDScheme>(getEnumValue(CondFDSchemeTypeNamesUC, Util::makeUPPER(s_ipsc->cAlphaArgs(1))));
181 1 : if (s_hbfd->CondFDSchemeType == CondFDScheme::Invalid) {
182 0 : ShowSevereError(state,
183 0 : format("{}: invalid {} entered={}, must match CrankNicholsonSecondOrder or FullyImplicitFirstOrder.",
184 0 : s_ipsc->cCurrentModuleObject,
185 0 : s_ipsc->cAlphaFieldNames(1),
186 0 : s_ipsc->cAlphaArgs(1)));
187 0 : ErrorsFound = true;
188 : }
189 : }
190 : }
191 :
192 1 : if (!s_ipsc->lNumericFieldBlanks(1)) {
193 1 : s_hbfd->SpaceDescritConstant = s_ipsc->rNumericArgs(1);
194 : }
195 1 : if (!s_ipsc->lNumericFieldBlanks(2)) {
196 1 : state.dataHeatBal->CondFDRelaxFactorInput = s_ipsc->rNumericArgs(2);
197 1 : state.dataHeatBal->CondFDRelaxFactor = state.dataHeatBal->CondFDRelaxFactorInput;
198 : }
199 1 : if (!s_ipsc->lNumericFieldBlanks(3)) {
200 1 : state.dataHeatBal->MaxAllowedDelTempCondFD = s_ipsc->rNumericArgs(3);
201 : }
202 :
203 : } // settings object
204 :
205 12 : pcMat = s_ip->getNumObjectsFound(state, "MaterialProperty:PhaseChange");
206 12 : vcMat = s_ip->getNumObjectsFound(state, "MaterialProperty:VariableThermalConductivity");
207 :
208 12 : int numProps = setSizeMaxProperties(state);
209 12 : MaterialProps.allocate(numProps);
210 :
211 12 : s_hbfd->MaterialFD.allocate(s_mat->materials.size());
212 :
213 : // Load the additional CondFD Material properties
214 12 : s_ipsc->cCurrentModuleObject = "MaterialProperty:PhaseChange"; // Phase Change Information First
215 :
216 12 : if (pcMat != 0) { // Get Phase Change info
217 : // CondFDVariableProperties = .TRUE.
218 2 : for (int Loop = 1; Loop <= pcMat; ++Loop) {
219 :
220 : // Call Input Get routine to retrieve material data
221 3 : s_ip->getObjectItem(state,
222 1 : s_ipsc->cCurrentModuleObject,
223 : Loop,
224 : MaterialNames,
225 : MaterialNumAlpha,
226 : MaterialProps,
227 : MaterialNumProp,
228 : IOStat,
229 1 : s_ipsc->lNumericFieldBlanks,
230 1 : s_ipsc->lAlphaFieldBlanks,
231 1 : s_ipsc->cAlphaFieldNames,
232 1 : s_ipsc->cNumericFieldNames);
233 :
234 : // Load the material derived type from the input data.
235 1 : ErrorObjectHeader eoh{routineName, s_ipsc->cCurrentModuleObject, MaterialNames(1)};
236 1 : int matNum = Material::GetMaterialNum(state, MaterialNames(1));
237 1 : if (matNum == 0) {
238 0 : ShowSevereItemNotFound(state, eoh, s_ipsc->cAlphaFieldNames(1), MaterialNames(1));
239 0 : ErrorsFound = true;
240 0 : continue;
241 : }
242 1 : auto const *mat = state.dataMaterial->materials(matNum);
243 :
244 1 : if (mat->group != Material::Group::Regular) {
245 0 : ShowSevereError(state,
246 0 : format("{}: Reference Material is not appropriate type for CondFD properties, material={}, must have regular "
247 : "properties (L,Cp,K,D)",
248 0 : s_ipsc->cCurrentModuleObject,
249 0 : mat->Name));
250 0 : ErrorsFound = true;
251 : }
252 :
253 : // Once the material derived type number is found then load the additional CondFD variable material properties
254 : // Some or all may be zero (default). They will be checked when calculating node temperatures
255 1 : auto &matFD = s_hbfd->MaterialFD(matNum);
256 1 : matFD.tk1 = MaterialProps(1);
257 1 : matFD.numTempEnth = (MaterialNumProp - 1) / 2;
258 1 : if (matFD.numTempEnth * 2 != (MaterialNumProp - 1)) {
259 0 : ShowSevereError(state, format("GetCondFDInput: {}=\"{}\", mismatched pairs", s_ipsc->cCurrentModuleObject, MaterialNames(1)));
260 0 : ShowContinueError(state, format("...expected {} pairs, but only entered {} numbers.", matFD.numTempEnth, MaterialNumProp - 1));
261 0 : ErrorsFound = true;
262 : }
263 1 : matFD.TempEnth.dimension(2, matFD.numTempEnth, 0.0);
264 1 : propNum = 2;
265 : // Temperature first
266 5 : for (int pcount = 1, pcount_end = matFD.numTempEnth; pcount <= pcount_end; ++pcount) {
267 4 : matFD.TempEnth(1, pcount) = MaterialProps(propNum);
268 4 : propNum += 2;
269 : }
270 1 : propNum = 3;
271 : // Then Enthalpy
272 5 : for (int pcount = 1, pcount_end = matFD.numTempEnth; pcount <= pcount_end; ++pcount) {
273 4 : matFD.TempEnth(2, pcount) = MaterialProps(propNum);
274 4 : propNum += 2;
275 : }
276 1 : nonInc = false;
277 1 : inegptr = 0;
278 4 : for (int pcount = 1, pcount_end = matFD.numTempEnth - 1; pcount <= pcount_end; ++pcount) {
279 3 : if (matFD.TempEnth(1, pcount) < matFD.TempEnth(1, pcount + 1)) {
280 3 : continue;
281 : }
282 0 : nonInc = true;
283 0 : inegptr = pcount + 1;
284 0 : break;
285 : }
286 1 : if (nonInc) {
287 0 : ShowSevereError(state,
288 0 : format("GetCondFDInput: {}=\"{}\", non increasing Temperatures. Temperatures must be strictly increasing.",
289 0 : s_ipsc->cCurrentModuleObject,
290 : MaterialNames(1)));
291 0 : ShowContinueError(state,
292 0 : format("...occurs first at item=[{}], value=[{:.2R}].", fmt::to_string(inegptr), matFD.TempEnth(1, inegptr)));
293 0 : ErrorsFound = true;
294 : }
295 1 : nonInc = false;
296 1 : inegptr = 0;
297 4 : for (int pcount = 1, pcount_end = matFD.numTempEnth - 1; pcount <= pcount_end; ++pcount) {
298 3 : if (matFD.TempEnth(2, pcount) <= matFD.TempEnth(2, pcount + 1)) {
299 3 : continue;
300 : }
301 0 : nonInc = true;
302 0 : inegptr = pcount + 1;
303 0 : break;
304 : }
305 1 : if (nonInc) {
306 0 : ShowSevereError(state,
307 0 : format("GetCondFDInput: {}=\"{}\", non increasing Enthalpy.", s_ipsc->cCurrentModuleObject, MaterialNames(1)));
308 0 : ShowContinueError(state, format("...occurs first at item=[{}], value=[{:.2R}].", inegptr, matFD.TempEnth(2, inegptr)));
309 0 : ShowContinueError(state, "...These values may be Cp (Specific Heat) rather than Enthalpy. Please correct.");
310 0 : ErrorsFound = true;
311 : }
312 : }
313 : }
314 : // Get CondFD Variable Thermal Conductivity Input
315 :
316 12 : s_ipsc->cCurrentModuleObject = "MaterialProperty:VariableThermalConductivity"; // Variable Thermal Conductivity Info next
317 12 : if (vcMat != 0) { // variable k info
318 : // CondFDVariableProperties = .TRUE.
319 2 : for (int Loop = 1; Loop <= vcMat; ++Loop) {
320 :
321 : // Call Input Get routine to retrieve material data
322 3 : s_ip->getObjectItem(state,
323 1 : s_ipsc->cCurrentModuleObject,
324 : Loop,
325 : MaterialNames,
326 : MaterialNumAlpha,
327 : MaterialProps,
328 : MaterialNumProp,
329 : IOStat,
330 1 : s_ipsc->lNumericFieldBlanks,
331 1 : s_ipsc->lAlphaFieldBlanks,
332 1 : s_ipsc->cAlphaFieldNames,
333 1 : s_ipsc->cNumericFieldNames);
334 :
335 1 : ErrorObjectHeader eoh{routineName, s_ipsc->cCurrentModuleObject, MaterialNames(1)};
336 : // Load the material derived type from the input data.
337 1 : int matNum = Material::GetMaterialNum(state, MaterialNames(1));
338 1 : if (matNum == 0) {
339 0 : ShowSevereItemNotFound(state, eoh, s_ipsc->cAlphaFieldNames(1), MaterialNames(1));
340 0 : ErrorsFound = true;
341 0 : continue;
342 : }
343 :
344 1 : auto *mat = s_mat->materials(matNum);
345 :
346 1 : if (mat->group != Material::Group::Regular) {
347 0 : ShowSevereError(state,
348 0 : format("{}: Reference Material is not appropriate type for CondFD properties, material={}, must have regular "
349 : "properties (L,Cp,K,D)",
350 0 : s_ipsc->cCurrentModuleObject,
351 0 : mat->Name));
352 0 : ErrorsFound = true;
353 : }
354 :
355 : // Once the material derived type number is found then load the additional CondFD variable material properties
356 : // Some or all may be zero (default). They will be checked when calculating node temperatures
357 :
358 1 : auto &matFD = s_hbfd->MaterialFD(matNum);
359 1 : matFD.numTempCond = MaterialNumProp / 2;
360 1 : if (matFD.numTempCond * 2 != MaterialNumProp) {
361 0 : ShowSevereError(state, format("GetCondFDInput: {}=\"{}\", mismatched pairs", s_ipsc->cCurrentModuleObject, MaterialNames(1)));
362 0 : ShowContinueError(state, format("...expected {} pairs, but only entered {} numbers.", matFD.numTempCond, MaterialNumProp));
363 0 : ErrorsFound = true;
364 : }
365 1 : matFD.TempCond.dimension(2, matFD.numTempCond, 0.0);
366 1 : propNum = 1;
367 : // Temperature first
368 5 : for (int pcount = 1, pcount_end = matFD.numTempCond; pcount <= pcount_end; ++pcount) {
369 4 : matFD.TempCond(1, pcount) = MaterialProps(propNum);
370 4 : propNum += 2;
371 : }
372 1 : propNum = 2;
373 : // Then Conductivity
374 5 : for (int pcount = 1, pcount_end = matFD.numTempCond; pcount <= pcount_end; ++pcount) {
375 4 : matFD.TempCond(2, pcount) = MaterialProps(propNum);
376 4 : propNum += 2;
377 : }
378 1 : nonInc = false;
379 1 : inegptr = 0;
380 4 : for (int pcount = 1, pcount_end = matFD.numTempCond - 1; pcount <= pcount_end; ++pcount) {
381 3 : if (matFD.TempCond(1, pcount) < matFD.TempCond(1, pcount + 1)) {
382 3 : continue;
383 : }
384 0 : nonInc = true;
385 0 : inegptr = pcount + 1;
386 0 : break;
387 : }
388 1 : if (nonInc) {
389 0 : ShowSevereError(state,
390 0 : format("GetCondFDInput: {}=\"{}\", non increasing Temperatures. Temperatures must be strictly increasing.",
391 0 : s_ipsc->cCurrentModuleObject,
392 : MaterialNames(1)));
393 0 : ShowContinueError(state, format("...occurs first at item=[{}], value=[{:.2R}].", inegptr, matFD.TempCond(1, inegptr)));
394 0 : ErrorsFound = true;
395 : }
396 : }
397 : }
398 :
399 403 : for (auto &matFD : s_hbfd->MaterialFD) {
400 391 : if (matFD.numTempEnth == 0) {
401 390 : matFD.numTempEnth = 3;
402 390 : matFD.TempEnth.dimension(2, 3, -100.0);
403 : }
404 391 : if (matFD.numTempCond == 0) {
405 390 : matFD.numTempCond = 3;
406 390 : matFD.TempCond.dimension(2, 3, -100.0);
407 : }
408 : }
409 :
410 12 : if (ErrorsFound) {
411 0 : ShowFatalError(state, "GetCondFDInput: Errors found getting ConductionFiniteDifference properties. Program terminates.");
412 : }
413 :
414 12 : InitialInitHeatBalFiniteDiff(state);
415 12 : }
416 :
417 12 : int setSizeMaxProperties(EnergyPlusData &state)
418 : {
419 : int numArgs;
420 : int numAlphas;
421 : int numNumerics;
422 12 : int maxTotalProps = 0;
423 :
424 12 : auto &s_ip = state.dataInputProcessing->inputProcessor;
425 :
426 12 : s_ip->getObjectDefMaxArgs(state, "MaterialProperty:PhaseChange", numArgs, numAlphas, numNumerics);
427 12 : maxTotalProps = max(maxTotalProps, numNumerics);
428 :
429 12 : s_ip->getObjectDefMaxArgs(state, "MaterialProperty:VariableThermalConductivity", numArgs, numAlphas, numNumerics);
430 12 : maxTotalProps = max(maxTotalProps, numNumerics);
431 :
432 12 : return maxTotalProps;
433 : }
434 :
435 158532 : void InitHeatBalFiniteDiff(EnergyPlusData &state)
436 : {
437 :
438 : // SUBROUTINE INFORMATION:
439 : // AUTHOR Richard J. Liesen
440 : // DATE WRITTEN Oct 2003
441 : // RE-ENGINEERED C O Pedersen 2006
442 : // B. Griffith May 2011 move begin-environment and every-timestep inits, cleanup formatting
443 :
444 : // PURPOSE OF THIS SUBROUTINE:
445 : // This subroutine sets the initial values for the FD moisture calculation
446 :
447 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
448 : bool ErrorsFound;
449 :
450 158532 : auto &s_hbfd = state.dataHeatBalFiniteDiffMgr;
451 :
452 158532 : if (s_hbfd->GetHBFiniteDiffInputFlag) {
453 : // Obtains conduction FD related parameters from input file
454 12 : GetCondFDInput(state);
455 12 : s_hbfd->GetHBFiniteDiffInputFlag = false;
456 : }
457 :
458 158532 : auto &SurfaceFD = s_hbfd->SurfaceFD;
459 158532 : ErrorsFound = false;
460 :
461 : // now do begin environment inits.
462 158532 : if (state.dataGlobal->BeginEnvrnFlag && s_hbfd->MyEnvrnFlag) {
463 1112 : for (int SurfNum = 1; SurfNum <= state.dataSurface->TotSurfaces; ++SurfNum) {
464 1036 : if (state.dataSurface->Surface(SurfNum).HeatTransferAlgorithm != DataSurfaces::HeatTransferModel::CondFD) {
465 68 : continue;
466 : }
467 968 : if (state.dataSurface->Surface(SurfNum).Construction <= 0) {
468 0 : continue; // Shading surface, not really a heat transfer surface
469 : }
470 968 : int ConstrNum = state.dataSurface->Surface(SurfNum).Construction;
471 968 : if (state.dataConstruction->Construct(ConstrNum).TypeIsWindow) {
472 0 : continue; // Windows simulated in Window module
473 : }
474 968 : auto &thisSurface = SurfaceFD(SurfNum);
475 968 : thisSurface.T = TempInitValue;
476 968 : thisSurface.TOld = TempInitValue;
477 968 : thisSurface.TT = TempInitValue;
478 968 : thisSurface.Rhov = RhovInitValue;
479 968 : thisSurface.RhovOld = RhovInitValue;
480 968 : thisSurface.RhoT = RhovInitValue;
481 968 : thisSurface.TD = TempInitValue;
482 968 : thisSurface.TDT = TempInitValue;
483 968 : thisSurface.TDTLast = TempInitValue;
484 968 : thisSurface.TDOld = TempInitValue;
485 968 : thisSurface.TDreport = TempInitValue;
486 968 : thisSurface.RH = 0.0;
487 968 : thisSurface.RHreport = 0.0;
488 968 : thisSurface.EnthOld = EnthInitValue;
489 968 : thisSurface.EnthNew = EnthInitValue;
490 968 : thisSurface.EnthLast = EnthInitValue;
491 968 : thisSurface.QDreport = 0.0;
492 968 : thisSurface.CpDelXRhoS1 = 0.0;
493 968 : thisSurface.CpDelXRhoS2 = 0.0;
494 968 : thisSurface.TDpriortimestep = 0.0;
495 968 : thisSurface.PhaseChangeState = Material::Phase::Transition;
496 968 : thisSurface.PhaseChangeStateOld = Material::Phase::Transition;
497 968 : thisSurface.PhaseChangeStateOldOld = Material::Phase::Transition;
498 968 : thisSurface.PhaseChangeStateRep = Material::phaseInts[(int)Material::Phase::Transition];
499 968 : thisSurface.PhaseChangeStateOldRep = Material::phaseInts[(int)Material::Phase::Transition];
500 968 : thisSurface.PhaseChangeStateOldOldRep = Material::phaseInts[(int)Material::Phase::Transition];
501 968 : thisSurface.PhaseChangeTemperatureReverse = 50;
502 :
503 968 : state.dataMstBal->TempOutsideAirFD(SurfNum) = 0.0;
504 968 : state.dataMstBal->RhoVaporAirOut(SurfNum) = 0.0;
505 968 : state.dataMstBal->RhoVaporSurfIn(SurfNum) = 0.0;
506 968 : state.dataMstBal->RhoVaporAirIn(SurfNum) = 0.0;
507 968 : state.dataMstBal->HConvExtFD(SurfNum) = 0.0;
508 968 : state.dataMstBal->HMassConvExtFD(SurfNum) = 0.0;
509 968 : state.dataMstBal->HConvInFD(SurfNum) = 0.0;
510 968 : state.dataMstBal->HMassConvInFD(SurfNum) = 0.0;
511 968 : state.dataMstBal->HSkyFD(SurfNum) = 0.0;
512 968 : state.dataMstBal->HGrndFD(SurfNum) = 0.0;
513 968 : state.dataMstBal->HAirFD(SurfNum) = 0.0;
514 : }
515 76 : s_hbfd->WarmupSurfTemp = 0;
516 76 : s_hbfd->MyEnvrnFlag = false;
517 : }
518 158532 : if (!state.dataGlobal->BeginEnvrnFlag) {
519 158456 : s_hbfd->MyEnvrnFlag = true;
520 : }
521 :
522 : // now do every timestep inits
523 :
524 1930572 : for (int SurfNum = 1; SurfNum <= state.dataSurface->TotSurfaces; ++SurfNum) {
525 1772040 : if (state.dataSurface->Surface(SurfNum).HeatTransferAlgorithm != DataSurfaces::HeatTransferModel::CondFD) {
526 117714 : continue;
527 : }
528 1654326 : if (state.dataSurface->Surface(SurfNum).Construction <= 0) {
529 0 : continue; // Shading surface, not really a heat transfer surface
530 : }
531 1654326 : int ConstrNum = state.dataSurface->Surface(SurfNum).Construction;
532 1654326 : if (state.dataConstruction->Construct(ConstrNum).TypeIsWindow) {
533 0 : continue; // Windows simulated in Window module
534 : }
535 1654326 : auto &thisSurface = SurfaceFD(SurfNum);
536 1654326 : thisSurface.T = thisSurface.TOld;
537 1654326 : thisSurface.Rhov = thisSurface.RhovOld;
538 1654326 : thisSurface.TD = thisSurface.TDOld;
539 1654326 : thisSurface.TDT = thisSurface.TDreport; // PT changes from TDold to TDreport
540 1654326 : thisSurface.TDTLast = thisSurface.TDOld;
541 1654326 : thisSurface.EnthOld = thisSurface.EnthOld;
542 1654326 : thisSurface.EnthNew = thisSurface.EnthOld;
543 1654326 : thisSurface.EnthLast = thisSurface.EnthOld;
544 1654326 : thisSurface.TDpriortimestep = thisSurface.TDreport; // Save TD for heat flux calc
545 : }
546 158532 : }
547 :
548 12 : void InitialInitHeatBalFiniteDiff(EnergyPlusData &state)
549 : {
550 :
551 : // SUBROUTINE INFORMATION:
552 : // AUTHOR Linda Lawrie
553 : // DATE WRITTEN March 2012
554 :
555 : // PURPOSE OF THIS SUBROUTINE:
556 : // This routine performs the original allocate, inits and setup output variables for the
557 : // module.
558 :
559 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
560 : Real64 dxn; // Intermediate calculation of nodal spacing. This is the full dx. There is
561 : // a half dxn thick node at each surface. dxn is the "capacitor" spacing.
562 : int Ipts1; // Intermediate calculation for number of full thickness nodes per layer. There
563 : // are always two half nodes at the layer faces.
564 : int Layer; // Loop counter
565 : int Delt;
566 :
567 : Real64 Alpha;
568 : Real64 mAlpha;
569 : Real64 StabilityTemp;
570 : Real64 StabilityMoist;
571 : Real64 kt;
572 : Real64 RhoS;
573 : Real64 Por;
574 : Real64 Cp;
575 : Real64 Dv;
576 : Real64 DeltaTimestep; // zone timestep in seconds, for local check of properties
577 : Real64 ThicknessThreshold; // min thickness consistent with other thermal properties, for local check
578 :
579 12 : auto &s_hbfd = state.dataHeatBalFiniteDiffMgr;
580 12 : auto &ConstructFD = s_hbfd->ConstructFD;
581 12 : auto &SigmaR = s_hbfd->SigmaR;
582 12 : auto &SigmaC = s_hbfd->SigmaC;
583 12 : auto &SurfaceFD = s_hbfd->SurfaceFD;
584 12 : auto &QHeatInFlux = s_hbfd->QHeatInFlux;
585 12 : auto &QHeatOutFlux = s_hbfd->QHeatOutFlux;
586 :
587 12 : ConstructFD.allocate(state.dataHeatBal->TotConstructs);
588 12 : SigmaR.allocate(state.dataHeatBal->TotConstructs);
589 12 : SigmaC.allocate(state.dataHeatBal->TotConstructs);
590 :
591 12 : SurfaceFD.allocate(state.dataSurface->TotSurfaces);
592 12 : QHeatInFlux.allocate(state.dataSurface->TotSurfaces);
593 12 : QHeatOutFlux.allocate(state.dataSurface->TotSurfaces);
594 :
595 : // And then initialize
596 160 : for (int SurfNum = 1; SurfNum <= state.dataSurface->TotSurfaces; ++SurfNum) {
597 148 : QHeatInFlux(SurfNum) = 0.0;
598 148 : QHeatOutFlux(SurfNum) = 0.0;
599 148 : state.dataHeatBalSurf->SurfOpaqInsFaceCondFlux(SurfNum) = 0.0;
600 148 : state.dataHeatBalSurf->SurfOpaqOutFaceCondFlux(SurfNum) = 0.0;
601 : }
602 :
603 : // Setup Output Variables
604 :
605 : // set a Delt that fits the zone time step and keeps it below 200s.
606 :
607 12 : s_hbfd->fracTimeStepZone_Hour = 1.0 / double(state.dataGlobal->TimeStepsInHour);
608 :
609 12 : for (int index = 1; index <= 20; ++index) {
610 12 : Delt = (s_hbfd->fracTimeStepZone_Hour * Constant::rSecsInHour) / index; // TimeStepZone = Zone time step in fractional hours
611 12 : if (Delt <= 200) {
612 12 : break;
613 : }
614 : }
615 :
616 63 : for (int ConstrNum = 1; ConstrNum <= state.dataHeatBal->TotConstructs; ++ConstrNum) {
617 51 : auto const &thisConstruct = state.dataConstruction->Construct(ConstrNum);
618 51 : auto &thisConstructFD = ConstructFD(ConstrNum);
619 : // Need to skip window constructions, IRT, air wall and construction not in use.
620 : // Need to also skip constructions for surfaces that do not use CondFD.
621 51 : if (thisConstruct.TypeIsWindow) {
622 6 : continue;
623 : }
624 45 : if (thisConstruct.TypeIsIRT) {
625 0 : continue;
626 : }
627 45 : if (thisConstruct.TypeIsAirBoundary) {
628 0 : continue;
629 : }
630 45 : if (!thisConstruct.IsUsed) {
631 0 : continue;
632 : }
633 45 : if (!findAnySurfacesUsingConstructionAndCondFD(state, ConstrNum)) {
634 3 : continue;
635 : }
636 :
637 42 : thisConstructFD.Name.allocate(thisConstruct.TotLayers);
638 42 : thisConstructFD.Thickness.allocate(thisConstruct.TotLayers);
639 42 : thisConstructFD.NodeNumPoint.allocate(thisConstruct.TotLayers);
640 42 : thisConstructFD.DelX.allocate(thisConstruct.TotLayers);
641 42 : thisConstructFD.TempStability.allocate(thisConstruct.TotLayers);
642 42 : thisConstructFD.MoistStability.allocate(thisConstruct.TotLayers);
643 :
644 42 : int TotNodes = 0;
645 42 : SigmaR(ConstrNum) = 0.0;
646 42 : SigmaC(ConstrNum) = 0.0;
647 :
648 165 : for (Layer = 1; Layer <= thisConstruct.TotLayers; ++Layer) { // Begin layer loop ...
649 :
650 : // Loop through all of the layers in the current construct. The purpose
651 : // of this loop is to define the thermal properties and to.
652 : // determine the total number of full size nodes in each layer.
653 : // The number of temperature points is one more than this
654 : // because of the two half nodes at the layer faces.
655 : // The calculation of dxn used here is based on a standard stability
656 : // criteria for explicit finite difference solutions. This criteria
657 : // was chosen not because it is viewed to be correct, but rather for
658 : // lack of any better criteria at this time. The use of a Fourier
659 : // number based criteria such as this is probably physically correct.
660 : // Change to implicit formulation still uses explicit stability, but
661 : // now there are special equations for R-only layers.
662 :
663 123 : int CurrentLayer = thisConstruct.LayerPoint(Layer);
664 123 : auto *mat = state.dataMaterial->materials(CurrentLayer);
665 :
666 123 : thisConstructFD.Name(Layer) = mat->Name;
667 123 : thisConstructFD.Thickness(Layer) = mat->Thickness;
668 :
669 : // Do some quick error checks for this section.
670 :
671 123 : if (mat->ROnly) { // Rlayer
672 :
673 : // These values are only needed temporarily and to calculate flux,
674 : // Layer will be handled
675 : // as a pure R in the temperature calc.
676 : // assign other properties based on resistance
677 :
678 8 : mat->SpecHeat = 0.0001;
679 8 : mat->Density = 1.0;
680 8 : mat->Thickness = 0.1; // arbitrary thickness for R layer
681 8 : mat->Conductivity = mat->Thickness / mat->Resistance;
682 8 : kt = mat->Conductivity;
683 8 : thisConstructFD.Thickness(Layer) = mat->Thickness;
684 :
685 8 : SigmaR(ConstrNum) += mat->Resistance; // add resistance of R layer
686 8 : SigmaC(ConstrNum) += 0.0; // no capacitance for R layer
687 :
688 8 : Alpha = kt / (mat->Density * mat->SpecHeat);
689 :
690 8 : mAlpha = 0.0;
691 :
692 115 : } else if (mat->group == Material::Group::AirGap) { // Group 1 = Air
693 :
694 : // Again, these values are only needed temporarily and to calculate flux,
695 : // Air layer will be handled
696 : // as a pure R in the temperature calc.
697 : // assign
698 : // other properties based on resistance
699 :
700 0 : mat->SpecHeat = 0.0001;
701 0 : mat->Density = 1.0;
702 0 : mat->Thickness = 0.1; // arbitrary thickness for R layer
703 0 : mat->Conductivity = mat->Thickness / mat->Resistance;
704 0 : kt = mat->Conductivity;
705 0 : thisConstructFD.Thickness(Layer) = mat->Thickness;
706 :
707 0 : SigmaR(ConstrNum) += mat->Resistance; // add resistance of R layer
708 0 : SigmaC(ConstrNum) += 0.0; // no capacitance for R layer
709 :
710 0 : Alpha = kt / (mat->Density * mat->SpecHeat);
711 0 : mAlpha = 0.0;
712 115 : } else if (thisConstruct.TypeIsIRT) { // make similar to air? (that didn't seem to work well)
713 0 : ShowSevereError(state,
714 0 : format("InitHeatBalFiniteDiff: Construction =\"{}\" uses Material:InfraredTransparent. Cannot be used currently "
715 : "with finite difference calculations.",
716 0 : thisConstruct.Name));
717 0 : if (thisConstruct.IsUsed) {
718 0 : ShowContinueError(state, "...since this construction is used in a surface, the simulation is not allowed.");
719 : } else {
720 0 : ShowContinueError(state, "...if this construction were used in a surface, the simulation would be terminated.");
721 : }
722 0 : continue;
723 : } else {
724 115 : auto *matReg = mat;
725 115 : assert(matReg != nullptr);
726 : // Regular material Properties
727 115 : kt = matReg->Conductivity;
728 115 : RhoS = matReg->Density;
729 115 : Por = matReg->Porosity;
730 115 : Cp = matReg->SpecHeat;
731 : // Need Resistance for reg layer
732 115 : matReg->Resistance = matReg->Thickness / matReg->Conductivity;
733 115 : Dv = matReg->VaporDiffus;
734 115 : SigmaR(ConstrNum) += matReg->Resistance; // add resistance
735 115 : SigmaC(ConstrNum) += matReg->Density * matReg->SpecHeat * matReg->Thickness;
736 115 : Alpha = kt / (RhoS * Cp);
737 115 : mAlpha = 0.0;
738 :
739 : // check for Material layers that are too thin and highly conductivity (not appropriate for surface models)
740 115 : if (Alpha > DataHeatBalance::HighDiffusivityThreshold) {
741 4 : DeltaTimestep = state.dataGlobal->TimeStepZoneSec;
742 4 : ThicknessThreshold = std::sqrt(Alpha * DeltaTimestep * 3.0);
743 4 : if (mat->Thickness < ThicknessThreshold) {
744 0 : ShowSevereError(
745 : state,
746 0 : format(
747 : "InitialInitHeatBalFiniteDiff: Found Material that is too thin and/or too highly conductive, material name = {}",
748 0 : mat->Name));
749 0 : ShowContinueError(state,
750 0 : format("High conductivity Material layers are not well supported by Conduction Finite Difference, "
751 : "material conductivity = {:.3R} [W/m-K]",
752 0 : mat->Conductivity));
753 0 : ShowContinueError(state, format("Material thermal diffusivity = {:.3R} [m2/s]", Alpha));
754 0 : ShowContinueError(
755 0 : state, format("Material with this thermal diffusivity should have thickness > {:.5R} [m]", ThicknessThreshold));
756 0 : if (mat->Thickness < DataHeatBalance::ThinMaterialLayerThreshold) {
757 0 : ShowContinueError(state,
758 0 : format("Material may be too thin to be modeled well, thickness = {:.5R} [m]", mat->Thickness));
759 0 : ShowContinueError(state,
760 0 : format("Material with this thermal diffusivity should have thickness > {:.5R} [m]",
761 : DataHeatBalance::ThinMaterialLayerThreshold));
762 : }
763 0 : ShowFatalError(state, "Preceding conditions cause termination.");
764 : }
765 : }
766 :
767 : } // R, Air or regular material properties and parameters
768 :
769 : // Proceed with setting node sizes in layers
770 :
771 123 : dxn = std::sqrt(Alpha * Delt * s_hbfd->SpaceDescritConstant); // The Fourier number is set using user constant
772 :
773 : // number of nodes=thickness/spacing. This is number of full size node spaces across layer.
774 123 : Ipts1 = int(mat->Thickness / dxn);
775 : // set high conductivity layers to a single full size node thickness. (two half nodes)
776 123 : if (Ipts1 <= 1) {
777 80 : Ipts1 = 1;
778 : }
779 123 : if (mat->ROnly || mat->group == Material::Group::AirGap) {
780 :
781 8 : Ipts1 = 1; // single full node in R layers- surfaces of adjacent material or inside/outside layer
782 : }
783 :
784 123 : dxn = mat->Thickness / double(Ipts1); // full node thickness
785 :
786 123 : StabilityTemp = Alpha * Delt / pow_2(dxn);
787 123 : StabilityMoist = mAlpha * Delt / pow_2(dxn);
788 123 : thisConstructFD.TempStability(Layer) = StabilityTemp;
789 123 : thisConstructFD.MoistStability(Layer) = StabilityMoist;
790 123 : thisConstructFD.DelX(Layer) = dxn;
791 :
792 123 : TotNodes += Ipts1; // number of full size nodes
793 123 : thisConstructFD.NodeNumPoint(Layer) = Ipts1; // number of full size nodes
794 : } // end of layer loop.
795 :
796 42 : thisConstructFD.TotNodes = TotNodes;
797 42 : thisConstructFD.DeltaTime = Delt;
798 :
799 : } // End of Construction Loop. TotNodes in each construction now set
800 :
801 : // now determine x location, or distance that nodes are from the outside face in meters
802 63 : for (int ConstrNum = 1; ConstrNum <= state.dataHeatBal->TotConstructs; ++ConstrNum) {
803 51 : auto &thisConstructFD = ConstructFD(ConstrNum);
804 51 : auto const &thisConstruct = state.dataConstruction->Construct(ConstrNum);
805 51 : if (thisConstructFD.TotNodes > 0) {
806 42 : thisConstructFD.NodeXlocation.allocate(thisConstructFD.TotNodes + 1);
807 42 : thisConstructFD.NodeXlocation = 0.0; // init them all
808 42 : Ipts1 = 0; // init counter
809 165 : for (Layer = 1; Layer <= thisConstruct.TotLayers; ++Layer) {
810 123 : int OutwardMatLayerNum = Layer - 1;
811 411 : for (int LayerNode = 1; LayerNode <= thisConstructFD.NodeNumPoint(Layer); ++LayerNode) {
812 288 : ++Ipts1;
813 288 : if (Ipts1 == 1) {
814 42 : thisConstructFD.NodeXlocation(Ipts1) = 0.0; // first node is on outside face
815 :
816 246 : } else if (LayerNode == 1) {
817 81 : if (OutwardMatLayerNum > 0 && OutwardMatLayerNum <= thisConstruct.TotLayers) {
818 : // later nodes are Delx away from previous, but use Delx from previous layer
819 81 : thisConstructFD.NodeXlocation(Ipts1) =
820 81 : thisConstructFD.NodeXlocation(Ipts1 - 1) + thisConstructFD.DelX(OutwardMatLayerNum);
821 : }
822 : } else {
823 : // later nodes are Delx away from previous
824 165 : thisConstructFD.NodeXlocation(Ipts1) = thisConstructFD.NodeXlocation(Ipts1 - 1) + thisConstructFD.DelX(Layer);
825 : }
826 : }
827 : }
828 42 : Layer = thisConstruct.TotLayers;
829 42 : ++Ipts1;
830 42 : thisConstructFD.NodeXlocation(Ipts1) = thisConstructFD.NodeXlocation(Ipts1 - 1) + thisConstructFD.DelX(Layer);
831 : }
832 : }
833 :
834 160 : for (int Surf = 1; Surf <= state.dataSurface->TotSurfaces; ++Surf) {
835 148 : if (!state.dataSurface->Surface(Surf).HeatTransSurf) {
836 0 : continue;
837 : }
838 148 : if (state.dataSurface->Surface(Surf).Class == DataSurfaces::SurfaceClass::Window) {
839 6 : continue;
840 : }
841 142 : if (state.dataSurface->Surface(Surf).HeatTransferAlgorithm != DataSurfaces::HeatTransferModel::CondFD) {
842 5 : continue;
843 : }
844 137 : int ConstrNum = state.dataSurface->Surface(Surf).Construction;
845 137 : int TotNodes = ConstructFD(ConstrNum).TotNodes;
846 137 : int TotLayers = state.dataConstruction->Construct(ConstrNum).TotLayers;
847 :
848 : // Allocate the Surface Arrays
849 137 : SurfaceFD(Surf).T.allocate(TotNodes + 1);
850 137 : SurfaceFD(Surf).TOld.allocate(TotNodes + 1);
851 137 : SurfaceFD(Surf).TT.allocate(TotNodes + 1);
852 137 : SurfaceFD(Surf).Rhov.allocate(TotNodes + 1);
853 137 : SurfaceFD(Surf).RhovOld.allocate(TotNodes + 1);
854 137 : SurfaceFD(Surf).RhoT.allocate(TotNodes + 1);
855 137 : SurfaceFD(Surf).TD.allocate(TotNodes + 1);
856 137 : SurfaceFD(Surf).TDT.allocate(TotNodes + 1);
857 137 : SurfaceFD(Surf).TDTLast.allocate(TotNodes + 1);
858 137 : SurfaceFD(Surf).TDOld.allocate(TotNodes + 1);
859 137 : SurfaceFD(Surf).TDreport.allocate(TotNodes + 1);
860 137 : SurfaceFD(Surf).RH.allocate(TotNodes + 1);
861 137 : SurfaceFD(Surf).RHreport.allocate(TotNodes + 1);
862 137 : SurfaceFD(Surf).EnthOld.allocate(TotNodes + 1);
863 137 : SurfaceFD(Surf).EnthNew.allocate(TotNodes + 1);
864 137 : SurfaceFD(Surf).EnthLast.allocate(TotNodes + 1);
865 137 : SurfaceFD(Surf).QDreport.allocate(TotNodes + 1);
866 137 : SurfaceFD(Surf).CpDelXRhoS1.allocate(TotNodes + 1);
867 137 : SurfaceFD(Surf).CpDelXRhoS2.allocate(TotNodes + 1);
868 137 : SurfaceFD(Surf).TDpriortimestep.allocate(TotNodes + 1);
869 137 : SurfaceFD(Surf).PhaseChangeState.allocate(TotNodes + 1);
870 137 : SurfaceFD(Surf).PhaseChangeStateOld.allocate(TotNodes + 1);
871 137 : SurfaceFD(Surf).PhaseChangeStateOldOld.allocate(TotNodes + 1);
872 137 : SurfaceFD(Surf).PhaseChangeStateRep.allocate(TotNodes + 1);
873 137 : SurfaceFD(Surf).PhaseChangeStateOldRep.allocate(TotNodes + 1);
874 137 : SurfaceFD(Surf).PhaseChangeStateOldOldRep.allocate(TotNodes + 1);
875 137 : SurfaceFD(Surf).PhaseChangeTemperatureReverse.allocate(TotNodes + 1);
876 137 : SurfaceFD(Surf).condMaterialActuators.allocate(TotLayers);
877 137 : SurfaceFD(Surf).specHeatMaterialActuators.allocate(TotLayers);
878 137 : SurfaceFD(Surf).condNodeReport.allocate(TotNodes + 1);
879 137 : SurfaceFD(Surf).specHeatNodeReport.allocate(TotNodes + 1);
880 137 : SurfaceFD(Surf).heatSourceFluxMaterialActuators.allocate(TotLayers - 1);
881 137 : SurfaceFD(Surf).heatSourceInternalFluxLayerReport.allocate(TotLayers - 1);
882 137 : SurfaceFD(Surf).heatSourceInternalFluxEnergyLayerReport.allocate(TotLayers - 1);
883 137 : SurfaceFD(Surf).heatSourceEMSFluxLayerReport.allocate(TotLayers - 1);
884 137 : SurfaceFD(Surf).heatSourceEMSFluxEnergyLayerReport.allocate(TotLayers - 1);
885 :
886 : // Initialize the allocated arrays.
887 137 : SurfaceFD(Surf).T = TempInitValue;
888 137 : SurfaceFD(Surf).TOld = TempInitValue;
889 137 : SurfaceFD(Surf).TT = TempInitValue;
890 137 : SurfaceFD(Surf).Rhov = RhovInitValue;
891 137 : SurfaceFD(Surf).RhovOld = RhovInitValue;
892 137 : SurfaceFD(Surf).RhoT = RhovInitValue;
893 137 : SurfaceFD(Surf).TD = TempInitValue;
894 137 : SurfaceFD(Surf).TDT = TempInitValue;
895 137 : SurfaceFD(Surf).TDTLast = TempInitValue;
896 137 : SurfaceFD(Surf).TDOld = TempInitValue;
897 137 : SurfaceFD(Surf).TDreport = TempInitValue;
898 137 : SurfaceFD(Surf).RH = 0.0;
899 137 : SurfaceFD(Surf).RHreport = 0.0;
900 137 : SurfaceFD(Surf).EnthOld = EnthInitValue;
901 137 : SurfaceFD(Surf).EnthNew = EnthInitValue;
902 137 : SurfaceFD(Surf).EnthLast = EnthInitValue;
903 137 : SurfaceFD(Surf).QDreport = 0.0;
904 137 : SurfaceFD(Surf).CpDelXRhoS1 = 0.0;
905 137 : SurfaceFD(Surf).CpDelXRhoS2 = 0.0;
906 137 : SurfaceFD(Surf).TDpriortimestep = 0.0;
907 137 : SurfaceFD(Surf).PhaseChangeState = Material::Phase::Transition;
908 137 : SurfaceFD(Surf).PhaseChangeStateOld = Material::Phase::Transition;
909 137 : SurfaceFD(Surf).PhaseChangeStateOldOld = Material::Phase::Transition;
910 137 : SurfaceFD(Surf).PhaseChangeStateRep = Material::phaseInts[(int)Material::Phase::Transition];
911 137 : SurfaceFD(Surf).PhaseChangeStateOldRep = Material::phaseInts[(int)Material::Phase::Transition];
912 137 : SurfaceFD(Surf).PhaseChangeStateOldOldRep = Material::phaseInts[(int)Material::Phase::Transition];
913 137 : SurfaceFD(Surf).PhaseChangeTemperatureReverse = 50;
914 137 : SurfaceFD(Surf).condNodeReport = 0.0;
915 137 : SurfaceFD(Surf).specHeatNodeReport = 0.0;
916 137 : SurfaceFD(Surf).heatSourceInternalFluxLayerReport = 0.0;
917 137 : SurfaceFD(Surf).heatSourceInternalFluxEnergyLayerReport = 0.0;
918 137 : SurfaceFD(Surf).heatSourceEMSFluxLayerReport = 0.0;
919 137 : SurfaceFD(Surf).heatSourceEMSFluxEnergyLayerReport = 0.0;
920 :
921 : // Setup EMS data
922 567 : for (int lay = 1; lay <= TotLayers; ++lay) {
923 : // Setup material layer names actuators
924 430 : int matLay = state.dataConstruction->Construct(ConstrNum).LayerPoint(lay);
925 : // Actuator name format: "{SurfName}:{MaterialLayerName}"
926 430 : std::string actName = fmt::format("{}:{}", state.dataSurface->Surface(Surf).Name, state.dataMaterial->materials(matLay)->Name);
927 430 : SurfaceFD(Surf).condMaterialActuators(lay).actuatorName = actName;
928 430 : SurfaceFD(Surf).specHeatMaterialActuators(lay).actuatorName = actName;
929 :
930 : // only setup for heat source actuator for layers 1 to N-1
931 430 : if (lay != TotLayers) {
932 293 : SurfaceFD(Surf).heatSourceFluxMaterialActuators(lay).actuatorName = actName;
933 : }
934 430 : }
935 : }
936 :
937 160 : for (int SurfNum = 1; SurfNum <= state.dataSurface->TotSurfaces; ++SurfNum) {
938 148 : if (!state.dataSurface->Surface(SurfNum).HeatTransSurf) {
939 0 : continue;
940 : }
941 148 : if (state.dataSurface->Surface(SurfNum).Class == DataSurfaces::SurfaceClass::Window) {
942 6 : continue;
943 : }
944 142 : if (state.dataSurface->Surface(SurfNum).HeatTransferAlgorithm != DataSurfaces::HeatTransferModel::CondFD) {
945 5 : continue;
946 : }
947 :
948 137 : SetupOutputVariable(state,
949 : "CondFD Inner Solver Loop Iteration Count",
950 : Constant::Units::None,
951 137 : SurfaceFD(SurfNum).GSloopCounter,
952 : OutputProcessor::TimeStepType::Zone,
953 : OutputProcessor::StoreType::Sum,
954 137 : state.dataSurface->Surface(SurfNum).Name);
955 :
956 : // Setup EMS Material Actuators for Conductivity and Specific Heat
957 137 : int ConstrNum = state.dataSurface->Surface(SurfNum).Construction;
958 :
959 137 : auto const &thisConstruct = state.dataConstruction->Construct(ConstrNum);
960 : // Setup internal heat source output variables
961 : // Only setup for layers 1 to N-1
962 430 : for (int lay = 1; lay < thisConstruct.TotLayers; ++lay) {
963 879 : SetupOutputVariable(state,
964 586 : format("CondFD Internal Heat Source Power After Layer {}", lay),
965 : Constant::Units::W,
966 293 : SurfaceFD(SurfNum).heatSourceInternalFluxLayerReport(lay),
967 : OutputProcessor::TimeStepType::Zone,
968 : OutputProcessor::StoreType::Average,
969 293 : state.dataSurface->Surface(SurfNum).Name);
970 879 : SetupOutputVariable(state,
971 586 : format("CondFD Internal Heat Source Energy After Layer {}", lay),
972 : Constant::Units::J,
973 293 : SurfaceFD(SurfNum).heatSourceInternalFluxEnergyLayerReport(lay),
974 : OutputProcessor::TimeStepType::Zone,
975 : OutputProcessor::StoreType::Sum,
976 293 : state.dataSurface->Surface(SurfNum).Name);
977 : }
978 :
979 137 : if (state.dataGlobal->AnyEnergyManagementSystemInModel) {
980 50 : for (int lay = 1; lay <= thisConstruct.TotLayers; ++lay) {
981 76 : EnergyPlus::SetupEMSActuator(state,
982 : "CondFD Surface Material Layer",
983 38 : SurfaceFD(SurfNum).condMaterialActuators(lay).actuatorName,
984 : "Thermal Conductivity",
985 : "[W/m-K]",
986 38 : SurfaceFD(SurfNum).condMaterialActuators(lay).isActuated,
987 38 : SurfaceFD(SurfNum).condMaterialActuators(lay).actuatedValue);
988 76 : EnergyPlus::SetupEMSActuator(state,
989 : "CondFD Surface Material Layer",
990 38 : SurfaceFD(SurfNum).specHeatMaterialActuators(lay).actuatorName,
991 : "Specific Heat",
992 : "[J/kg-C]",
993 38 : SurfaceFD(SurfNum).specHeatMaterialActuators(lay).isActuated,
994 38 : SurfaceFD(SurfNum).specHeatMaterialActuators(lay).actuatedValue);
995 : }
996 :
997 : // Setup EMS Actuator and Output Variables for Heat Flux
998 : // Only setup for layers 1 to N-1
999 38 : for (int lay = 1; lay < thisConstruct.TotLayers; ++lay) {
1000 52 : EnergyPlus::SetupEMSActuator(state,
1001 : "CondFD Surface Material Layer",
1002 26 : SurfaceFD(SurfNum).heatSourceFluxMaterialActuators(lay).actuatorName,
1003 : "Heat Flux",
1004 : "[W/m2]",
1005 26 : SurfaceFD(SurfNum).heatSourceFluxMaterialActuators(lay).isActuated,
1006 26 : SurfaceFD(SurfNum).heatSourceFluxMaterialActuators(lay).actuatedValue);
1007 78 : SetupOutputVariable(state,
1008 52 : format("CondFD EMS Heat Source Power After Layer {}", lay),
1009 : Constant::Units::W,
1010 26 : SurfaceFD(SurfNum).heatSourceEMSFluxLayerReport(lay),
1011 : OutputProcessor::TimeStepType::Zone,
1012 : OutputProcessor::StoreType::Average,
1013 26 : state.dataSurface->Surface(SurfNum).Name);
1014 78 : SetupOutputVariable(state,
1015 52 : format("CondFD EMS Heat Source Energy After Layer {}", lay),
1016 : Constant::Units::J,
1017 26 : SurfaceFD(SurfNum).heatSourceEMSFluxEnergyLayerReport(lay),
1018 : OutputProcessor::TimeStepType::Zone,
1019 : OutputProcessor::StoreType::Sum,
1020 26 : state.dataSurface->Surface(SurfNum).Name,
1021 : Constant::eResource::Electricity,
1022 : OutputProcessor::Group::Building,
1023 : OutputProcessor::EndUseCat::Heating);
1024 : }
1025 : }
1026 :
1027 137 : int TotNodes = ConstructFD(state.dataSurface->Surface(SurfNum).Construction).TotNodes; // Full size nodes, start with outside face.
1028 1378 : for (int node = 1; node <= TotNodes + 1; ++node) { // include inside face node
1029 3723 : SetupOutputVariable(state,
1030 2482 : format("CondFD Surface Temperature Node {}", node),
1031 : Constant::Units::C,
1032 1241 : SurfaceFD(SurfNum).TDreport(node),
1033 : OutputProcessor::TimeStepType::Zone,
1034 : OutputProcessor::StoreType::Average,
1035 1241 : state.dataSurface->Surface(SurfNum).Name);
1036 3723 : SetupOutputVariable(state,
1037 2482 : format("CondFD Surface Heat Flux Node {}", node),
1038 : Constant::Units::W_m2,
1039 1241 : SurfaceFD(SurfNum).QDreport(node),
1040 : OutputProcessor::TimeStepType::Zone,
1041 : OutputProcessor::StoreType::Average,
1042 1241 : state.dataSurface->Surface(SurfNum).Name);
1043 1241 : SetupOutputVariable(state,
1044 2482 : format("CondFD Phase Change State {}", node),
1045 : Constant::Units::None,
1046 1241 : (int &)SurfaceFD(SurfNum).PhaseChangeStateRep(node),
1047 : OutputProcessor::TimeStepType::Zone,
1048 : OutputProcessor::StoreType::Average,
1049 1241 : state.dataSurface->Surface(SurfNum).Name);
1050 1241 : SetupOutputVariable(state,
1051 2482 : format("CondFD Phase Change Previous State {}", node),
1052 : Constant::Units::None,
1053 1241 : (int &)SurfaceFD(SurfNum).PhaseChangeStateOldRep(node),
1054 : OutputProcessor::TimeStepType::Zone,
1055 : OutputProcessor::StoreType::Average,
1056 1241 : state.dataSurface->Surface(SurfNum).Name);
1057 3723 : SetupOutputVariable(state,
1058 2482 : format("CondFD Phase Change Node Temperature {}", node),
1059 : Constant::Units::C,
1060 1241 : SurfaceFD(SurfNum).TDT(node),
1061 : OutputProcessor::TimeStepType::Zone,
1062 : OutputProcessor::StoreType::Average,
1063 1241 : state.dataSurface->Surface(SurfNum).Name);
1064 3723 : SetupOutputVariable(state,
1065 2482 : format("CondFD Phase Change Node Conductivity {}", node),
1066 : Constant::Units::W_mK,
1067 1241 : SurfaceFD(SurfNum).condNodeReport(node),
1068 : OutputProcessor::TimeStepType::Zone,
1069 : OutputProcessor::StoreType::Average,
1070 1241 : state.dataSurface->Surface(SurfNum).Name);
1071 3723 : SetupOutputVariable(state,
1072 2482 : format("CondFD Phase Change Node Specific Heat {}", node),
1073 : Constant::Units::J_kgK,
1074 1241 : SurfaceFD(SurfNum).specHeatNodeReport(node),
1075 : OutputProcessor::TimeStepType::Zone,
1076 : OutputProcessor::StoreType::Average,
1077 1241 : state.dataSurface->Surface(SurfNum).Name);
1078 1241 : if (state.dataGlobal->DisplayAdvancedReportVariables) {
1079 303 : SetupOutputVariable(state,
1080 202 : format("CondFD Surface Heat Capacitance Outer Half Node {}", node),
1081 : Constant::Units::W_m2K,
1082 101 : SurfaceFD(SurfNum).CpDelXRhoS1(node),
1083 : OutputProcessor::TimeStepType::Zone,
1084 : OutputProcessor::StoreType::Average,
1085 101 : state.dataSurface->Surface(SurfNum).Name);
1086 303 : SetupOutputVariable(state,
1087 202 : format("CondFD Surface Heat Capacitance Inner Half Node {}", node),
1088 : Constant::Units::W_m2K,
1089 101 : SurfaceFD(SurfNum).CpDelXRhoS2(node),
1090 : OutputProcessor::TimeStepType::Zone,
1091 : OutputProcessor::StoreType::Average,
1092 101 : state.dataSurface->Surface(SurfNum).Name);
1093 : }
1094 : }
1095 :
1096 : } // End of the Surface Loop for Report Variable Setup
1097 :
1098 12 : ReportFiniteDiffInits(state); // Report the results from the Finite Diff Inits
1099 12 : }
1100 :
1101 13 : int numNodesInMaterialLayer(EnergyPlusData &state, std::string const &surfName, std::string const &matName)
1102 : {
1103 13 : auto &s_hbfd = state.dataHeatBalFiniteDiffMgr;
1104 91 : for (auto const &surface : state.dataSurface->Surface) {
1105 78 : if (surface.Name == surfName) {
1106 13 : int constrNum = surface.Construction;
1107 66 : for (int lay = 1; lay <= state.dataConstruction->Construct(constrNum).TotLayers; ++lay) {
1108 53 : int matLay = state.dataConstruction->Construct(constrNum).LayerPoint(lay);
1109 53 : if (state.dataMaterial->materials(matLay)->Name == matName) {
1110 0 : return s_hbfd->ConstructFD(constrNum).NodeNumPoint(lay);
1111 : }
1112 : }
1113 : }
1114 13 : }
1115 :
1116 13 : return 0;
1117 : }
1118 :
1119 1374679 : void relax_array(Array1D<Real64> &a, // Array to relax
1120 : Array1D<Real64> const &b, // Array to relax towards
1121 : Real64 const r // Relaxation factor [0-1]
1122 : )
1123 : {
1124 1374679 : assert(equal_dimensions(a, b));
1125 1374679 : assert((0.0 <= r) && (r <= 1.0));
1126 1374679 : Real64 const q(1.0 - r);
1127 17266133 : for (int i = a.l(), e = a.u(); i <= e; ++i) {
1128 15891454 : a(i) = r * b(i) + q * a(i);
1129 : }
1130 1374679 : }
1131 :
1132 14564705 : Real64 sum_array_diff(Array1D<Real64> const &a, Array1D<Real64> const &b)
1133 : {
1134 14564705 : assert(equal_dimensions(a, b));
1135 14564705 : Real64 s(0.0);
1136 165114072 : for (int i = a.l(), e = a.u(); i <= e; ++i) {
1137 150549367 : s += a(i) - b(i); //? Should this be in abs?
1138 : }
1139 14564705 : return s;
1140 : }
1141 :
1142 10240031 : void CalcHeatBalFiniteDiff(EnergyPlusData &state,
1143 : int const Surf, // Surface number
1144 : Real64 &SurfTempInTmp, // INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
1145 : Real64 &TempSurfOutTmp // Outside Surface Temperature of each Heat Transfer Surface
1146 : )
1147 : {
1148 :
1149 : // SUBROUTINE INFORMATION:
1150 : // AUTHOR Richard J. Liesen
1151 : // DATE WRITTEN Oct 2003
1152 : // MODIFIED Aug 2006 by C O Pedersen to include implicit solution and variable properties with
1153 : // material enthalpy added for Phase Change Materials.
1154 : // Sept 2010 B. Griffith, remove allocate/deallocate, use structure variables
1155 : // March 2011 P. Tabares, add relaxation factor and add surfIteration to
1156 : // update TD and TDT, correct interzone partition
1157 : // May 2011 B. Griffith add logging and errors when inner GS loop does not converge
1158 : // November 2011 P. Tabares fixed problems with adiabatic walls/massless walls and PCM stability problems
1159 :
1160 : // PURPOSE OF THIS SUBROUTINE:
1161 : // this routine controls the calculation of the fluxes and temperatures using
1162 : // finite difference procedures for
1163 : // all building surface constructs.
1164 :
1165 10240031 : int const ConstrNum = state.dataSurface->Surface(Surf).Construction;
1166 :
1167 10240031 : auto &s_hbfd = state.dataHeatBalFiniteDiffMgr;
1168 10240031 : auto &constructFD = s_hbfd->ConstructFD(ConstrNum);
1169 :
1170 10240031 : int const TotNodes = constructFD.TotNodes;
1171 10240031 : int const TotLayers = state.dataConstruction->Construct(ConstrNum).TotLayers;
1172 :
1173 10240031 : SurfTempInTmp = 0.0;
1174 10240031 : TempSurfOutTmp = 0.0;
1175 :
1176 10240031 : int const Delt = constructFD.DeltaTime; // (seconds)
1177 :
1178 : // Aliases
1179 10240031 : auto &surfaceFD = s_hbfd->SurfaceFD(Surf);
1180 :
1181 10240031 : Real64 HMovInsul = 0;
1182 10240031 : if (state.dataSurface->AnyMovableInsulation) {
1183 0 : HMovInsul = state.dataSurface->extMovInsuls(Surf).H; // Even if this is not a movable insulation surface?
1184 : }
1185 : // Start stepping through the slab with time.
1186 20480062 : for (int J = 1, J_end = nint(state.dataGlobal->TimeStepZoneSec / Delt); J <= J_end; ++J) { // PT testing higher time steps
1187 :
1188 : int GSiter; // iteration counter for implicit repeat calculation
1189 35044771 : for (GSiter = 1; GSiter <= s_hbfd->MaxGSiter; ++GSiter) { // Iterate implicit equations
1190 35044767 : surfaceFD.TDTLast = surfaceFD.TDT; // Save last iteration's TDT (New temperature) values
1191 35044767 : surfaceFD.EnthLast = surfaceFD.EnthNew; // Last iterations new enthalpy value
1192 :
1193 : // Original loop version
1194 35044767 : int i(1); // Node counter
1195 153647293 : for (int Lay = 1; Lay <= TotLayers; ++Lay) { // Begin layer loop ...
1196 :
1197 : // For the exterior surface node with a convective boundary condition
1198 118602526 : if ((i == 1) && (Lay == 1)) {
1199 35044767 : ExteriorBCEqns(state,
1200 : Delt,
1201 : i,
1202 : Lay,
1203 : Surf,
1204 35044767 : surfaceFD.T,
1205 35044767 : surfaceFD.TT,
1206 35044767 : surfaceFD.Rhov,
1207 35044767 : surfaceFD.RhoT,
1208 35044767 : surfaceFD.RH,
1209 35044767 : surfaceFD.TD,
1210 35044767 : surfaceFD.TDT,
1211 35044767 : surfaceFD.EnthOld,
1212 35044767 : surfaceFD.EnthNew,
1213 : TotNodes,
1214 : HMovInsul);
1215 : }
1216 :
1217 : // For the Layer Interior nodes. Arrive here after exterior surface node or interface node
1218 118602526 : if (TotNodes != 1) {
1219 324615375 : for (int ctr = 2, ctr_end = constructFD.NodeNumPoint(Lay); ctr <= ctr_end; ++ctr) {
1220 207129160 : ++i;
1221 207129160 : InteriorNodeEqns(state,
1222 : Delt,
1223 : i,
1224 : Lay,
1225 : Surf,
1226 207129160 : surfaceFD.T,
1227 207129160 : surfaceFD.TT,
1228 207129160 : surfaceFD.Rhov,
1229 207129160 : surfaceFD.RhoT,
1230 207129160 : surfaceFD.RH,
1231 207129160 : surfaceFD.TD,
1232 207129160 : surfaceFD.TDT,
1233 207129160 : surfaceFD.EnthOld,
1234 207129160 : surfaceFD.EnthNew);
1235 : }
1236 : }
1237 :
1238 118602526 : if ((Lay < TotLayers) && (TotNodes != 1)) { // Interface equations for 2 capacitive materials
1239 83557759 : ++i;
1240 83557759 : IntInterfaceNodeEqns(state,
1241 : Delt,
1242 : i,
1243 : Lay,
1244 : Surf,
1245 83557759 : surfaceFD.T,
1246 83557759 : surfaceFD.TT,
1247 83557759 : surfaceFD.Rhov,
1248 83557759 : surfaceFD.RhoT,
1249 83557759 : surfaceFD.RH,
1250 83557759 : surfaceFD.TD,
1251 83557759 : surfaceFD.TDT,
1252 83557759 : surfaceFD.EnthOld,
1253 83557759 : surfaceFD.EnthNew,
1254 : GSiter);
1255 35044767 : } else if (Lay == TotLayers) { // For the Interior surface node with a convective boundary condition
1256 35044767 : ++i;
1257 35044767 : InteriorBCEqns(state,
1258 : Delt,
1259 : i,
1260 : Lay,
1261 : Surf,
1262 35044767 : surfaceFD.T,
1263 35044767 : surfaceFD.TT,
1264 35044767 : surfaceFD.Rhov,
1265 35044767 : surfaceFD.RhoT,
1266 35044767 : surfaceFD.RH,
1267 35044767 : surfaceFD.TD,
1268 35044767 : surfaceFD.TDT,
1269 35044767 : surfaceFD.EnthOld,
1270 35044767 : surfaceFD.EnthNew,
1271 35044767 : surfaceFD.TDreport);
1272 : }
1273 :
1274 : } // layer loop
1275 :
1276 : // Apply Relaxation factor for stability, use current (TDT) and previous (TDTLast) iteration temperature values
1277 : // to obtain the actual temperature that is going to be used for next iteration. This would mostly happen with PCM
1278 : // Tuned Function call to eliminate array temporaries and multiple relaxation passes
1279 35044767 : if (GSiter > 15) {
1280 24696 : relax_array(surfaceFD.TDT, surfaceFD.TDTLast, 0.9875);
1281 35020071 : } else if (GSiter > 10) {
1282 190771 : relax_array(surfaceFD.TDT, surfaceFD.TDTLast, 0.875);
1283 34829300 : } else if (GSiter > 5) {
1284 1149270 : relax_array(surfaceFD.TDT, surfaceFD.TDTLast, 0.5);
1285 : }
1286 :
1287 : // the following could blow up when all the node temps sum to less than 1.0. seems poorly formulated for temperature in C.
1288 : // PT delete one zero and decrease number of minimum iterations, from 3 (which actually requires 4 iterations) to 2.
1289 :
1290 35044767 : if ((GSiter > 2) && (std::abs(sum_array_diff(surfaceFD.TDT, surfaceFD.TDTLast) / sum(surfaceFD.TDT)) < 0.00001)) {
1291 10240027 : break;
1292 : }
1293 :
1294 : } // End of Gauss Seidell iteration loop
1295 :
1296 10240031 : surfaceFD.GSloopCounter = GSiter; // outputs GSloop iterations, useful for pinpointing stability issues with condFD
1297 10240031 : if (state.dataHeatBal->CondFDRelaxFactor != 1.0) {
1298 : // Apply Relaxation factor for stability, use current (TDT) and previous (TDreport) temperature values
1299 : // to obtain the actual temperature that is going to be exported/use
1300 9942 : relax_array(surfaceFD.TDT, surfaceFD.TDreport, 1.0 - state.dataHeatBal->CondFDRelaxFactor);
1301 9942 : surfaceFD.EnthOld = surfaceFD.EnthNew;
1302 : }
1303 :
1304 115353574 : for (int I = 1; I <= (TotNodes + 1); I++) {
1305 : // When the phase change process reverses its direction while melting or freezing (without completing its phase
1306 : // to either liquid or solid), the temperature at which it changes its direction is saved
1307 : // in the variable PhaseChangeTemperatureReverse, and this variable will hold the value of the temperature until
1308 : // the next reverse in the process takes place.
1309 105113543 : if (((surfaceFD.PhaseChangeStateOld(I) == Material::Phase::Freezing &&
1310 4214 : surfaceFD.PhaseChangeState(I) == Material::Phase::Transition) ||
1311 105113543 : (surfaceFD.PhaseChangeStateOld(I) == Material::Phase::Transition &&
1312 315196288 : surfaceFD.PhaseChangeState(I) == Material::Phase::Freezing)) ||
1313 105113543 : ((surfaceFD.PhaseChangeStateOld(I) == Material::Phase::Melting && surfaceFD.PhaseChangeState(I) == Material::Phase::Transition) ||
1314 105113543 : (surfaceFD.PhaseChangeStateOld(I) == Material::Phase::Transition &&
1315 104969202 : surfaceFD.PhaseChangeState(I) == Material::Phase::Melting))) {
1316 0 : surfaceFD.PhaseChangeTemperatureReverse(I) = surfaceFD.TDT(I);
1317 : }
1318 : }
1319 :
1320 10240031 : surfaceFD.PhaseChangeStateOldOld = surfaceFD.PhaseChangeStateOld;
1321 10240031 : surfaceFD.PhaseChangeStateOld = surfaceFD.PhaseChangeState;
1322 :
1323 10240031 : surfaceFD.PhaseChangeStateOldOldRep = surfaceFD.PhaseChangeStateOldRep;
1324 10240031 : surfaceFD.PhaseChangeStateOldRep = surfaceFD.PhaseChangeStateRep;
1325 : } // Time Loop //PT solving time steps
1326 :
1327 10240031 : TempSurfOutTmp = surfaceFD.TDT(1);
1328 10240031 : SurfTempInTmp = surfaceFD.TDT(TotNodes + 1);
1329 10240031 : state.dataMstBal->RhoVaporSurfIn(Surf) = 0.0;
1330 :
1331 : // For ground surfaces or when raining, outside face inner half-node heat capacity was unknown and set to -1 in ExteriorBCEqns
1332 : // Now check for the flag and set equal to the second node's outer half-node heat capacity if needed
1333 10240031 : if (surfaceFD.CpDelXRhoS2(1) == -1.0) {
1334 1206364 : surfaceFD.CpDelXRhoS2(1) = surfaceFD.CpDelXRhoS1(2); // Set to node 2's outer half node heat capacity
1335 : }
1336 10240031 : CalcNodeHeatFlux(state, Surf, TotNodes);
1337 :
1338 : // Determine largest change in node temps
1339 10240031 : Real64 MaxDelTemp = 0.0;
1340 115353574 : for (int NodeNum = 1; NodeNum <= TotNodes + 1; ++NodeNum) { // need to consider all nodes
1341 105113543 : MaxDelTemp = max(std::abs(surfaceFD.TDT(NodeNum) - surfaceFD.TDreport(NodeNum)), MaxDelTemp);
1342 : }
1343 10240031 : surfaceFD.MaxNodeDelTemp = MaxDelTemp;
1344 10240031 : surfaceFD.TDreport = surfaceFD.TDT;
1345 10240031 : surfaceFD.EnthOld = surfaceFD.EnthNew;
1346 10240031 : }
1347 :
1348 12 : void ReportFiniteDiffInits(EnergyPlusData &state)
1349 : {
1350 :
1351 : // SUBROUTINE INFORMATION:
1352 : // AUTHOR Richard Liesen
1353 : // DATE WRITTEN November 2003
1354 : // MODIFIED B. Griffith, May 2011 add reporting of node x locations
1355 :
1356 : // PURPOSE OF THIS SUBROUTINE:
1357 : // This routine gives a detailed report to the user about
1358 : // the initializations for the Finite Difference calculations
1359 : // of each construction.
1360 :
1361 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
1362 : bool DoReport;
1363 :
1364 12 : auto &s_hbfd = state.dataHeatBalFiniteDiffMgr;
1365 :
1366 : // Formats
1367 : static constexpr std::string_view Format_702(" ConductionFiniteDifference Node,{},{:.8R},{},{},{}\n");
1368 :
1369 12 : print(state.files.eio,
1370 : "! <ConductionFiniteDifference HeatBalanceSettings>,Scheme Type,Space Discretization Constant,Relaxation Factor,Inside Face Surface "
1371 : "Temperature Convergence Criteria\n");
1372 12 : print(state.files.eio,
1373 : " ConductionFiniteDifference HeatBalanceSettings,{},{:.2R},{:.2R},{:.4R}\n",
1374 12 : CondFDSchemeTypeNamesCC[static_cast<int>(s_hbfd->CondFDSchemeType)],
1375 12 : s_hbfd->SpaceDescritConstant,
1376 12 : state.dataHeatBal->CondFDRelaxFactorInput,
1377 12 : state.dataHeatBal->MaxAllowedDelTempCondFD);
1378 :
1379 36 : General::ScanForReports(state, "Constructions", DoReport, "Constructions");
1380 :
1381 12 : if (DoReport) {
1382 :
1383 : // Write Descriptions
1384 5 : print(state.files.eio, "{}\n", "! <Construction CondFD>,Construction Name,Index,#Layers,#Nodes,Time Step {hours}");
1385 5 : print(state.files.eio,
1386 : "{}\n",
1387 : "! <Material CondFD Summary>,Material Name,Thickness {m},#Layer Elements,Layer Delta X,Layer Alpha*Delt/Delx**2,Layer Moisture "
1388 : "Stability");
1389 :
1390 : // HT Algo issue
1391 5 : if (state.dataHeatBal->AnyCondFD) {
1392 5 : print(state.files.eio,
1393 : "{}\n",
1394 : "! <ConductionFiniteDifference Node>,Node Identifier, Node Distance From Outside Face {m}, Construction Name, Outward Material "
1395 : "Name (or Face), Inward Material Name (or Face)");
1396 : }
1397 :
1398 25 : for (int ThisNum = 1; ThisNum <= state.dataHeatBal->TotConstructs; ++ThisNum) {
1399 20 : auto &construct = state.dataConstruction->Construct(ThisNum);
1400 :
1401 20 : if (construct.TypeIsWindow) {
1402 4 : continue;
1403 : }
1404 19 : if (construct.TypeIsIRT) {
1405 0 : continue;
1406 : }
1407 19 : if (construct.TypeIsAirBoundary) {
1408 0 : continue;
1409 : }
1410 19 : if (!construct.IsUsed) {
1411 0 : continue;
1412 : }
1413 19 : if (!findAnySurfacesUsingConstructionAndCondFD(state, ThisNum)) {
1414 3 : continue;
1415 : }
1416 :
1417 16 : auto &constructFD = s_hbfd->ConstructFD(ThisNum);
1418 : static constexpr std::string_view Format_700(" Construction CondFD,{},{},{},{},{:.6R}\n");
1419 16 : print(state.files.eio,
1420 : Format_700,
1421 16 : construct.Name,
1422 : ThisNum,
1423 16 : construct.TotLayers,
1424 0 : int(constructFD.TotNodes + 1),
1425 16 : constructFD.DeltaTime / Constant::rSecsInHour);
1426 :
1427 52 : for (int Layer = 1; Layer <= construct.TotLayers; ++Layer) {
1428 : static constexpr std::string_view Format_701(" Material CondFD Summary,{},{:.4R},{},{:.8R},{:.8R},{:.8R}\n");
1429 36 : print(state.files.eio,
1430 : Format_701,
1431 : constructFD.Name(Layer),
1432 : constructFD.Thickness(Layer),
1433 : constructFD.NodeNumPoint(Layer),
1434 : constructFD.DelX(Layer),
1435 : constructFD.TempStability(Layer),
1436 : constructFD.MoistStability(Layer));
1437 : }
1438 :
1439 : // now list each CondFD Node with its X distance from outside face in m along with other identifiers
1440 16 : int Inodes = 0;
1441 :
1442 52 : for (int Layer = 1; Layer <= construct.TotLayers; ++Layer) {
1443 36 : int OutwardMatLayerNum = Layer - 1;
1444 110 : for (int LayerNode = 1; LayerNode <= constructFD.NodeNumPoint(Layer); ++LayerNode) {
1445 74 : ++Inodes;
1446 74 : if (Inodes == 1) {
1447 32 : print(state.files.eio,
1448 : Format_702,
1449 32 : format("Node #{}", Inodes),
1450 : constructFD.NodeXlocation(Inodes),
1451 16 : construct.Name,
1452 : "Surface Outside Face",
1453 : constructFD.Name(Layer));
1454 :
1455 58 : } else if (LayerNode == 1) {
1456 :
1457 20 : if (OutwardMatLayerNum > 0 && OutwardMatLayerNum <= construct.TotLayers) {
1458 40 : print(state.files.eio,
1459 : Format_702,
1460 40 : format("Node #{}", Inodes),
1461 : constructFD.NodeXlocation(Inodes),
1462 20 : construct.Name,
1463 : constructFD.Name(OutwardMatLayerNum),
1464 : constructFD.Name(Layer));
1465 : }
1466 38 : } else if (LayerNode > 1) {
1467 38 : OutwardMatLayerNum = Layer;
1468 76 : print(state.files.eio,
1469 : Format_702,
1470 76 : format("Node #{}", Inodes),
1471 : constructFD.NodeXlocation(Inodes),
1472 38 : construct.Name,
1473 : constructFD.Name(OutwardMatLayerNum),
1474 : constructFD.Name(Layer));
1475 : }
1476 : }
1477 : }
1478 :
1479 16 : int Layer = construct.TotLayers;
1480 16 : ++Inodes;
1481 32 : print(state.files.eio,
1482 : Format_702,
1483 32 : format("Node #{}", Inodes),
1484 : constructFD.NodeXlocation(Inodes),
1485 16 : construct.Name,
1486 : constructFD.Name(Layer),
1487 : "Surface Inside Face");
1488 : }
1489 : }
1490 12 : }
1491 :
1492 2350788 : Real64 terpld(Array2<Real64> const &a, Real64 const x1, int const nind, int const ndep)
1493 : {
1494 : // author:c. o. pedersen
1495 : // purpose:
1496 : // this function performs a linear interpolation
1497 : // on a two dimensional array containing both
1498 : // dependent and independent variables.
1499 :
1500 : // inputs:
1501 : // a = two dimensional array
1502 : // nind=row containing independent variable
1503 : // ndep=row containing the dependent variable
1504 : // x1 = specific independent variable value for which
1505 : // interpolated output is wanted
1506 : // outputs:
1507 : // the value of dependent variable corresponding
1508 : // to x1
1509 : // routine returns first or last dependent variable
1510 : // for out of range x1.
1511 :
1512 2350788 : int const first(a.l2());
1513 :
1514 2350788 : assert(a.size() > 0u);
1515 2350788 : Array2<Real64>::size_type l(1);
1516 2350788 : Real64 r(a[0]);
1517 2350788 : int last(first);
1518 9403152 : for (int i1 = first + 1, e1 = a.u2(); i1 <= e1; ++i1, ++l) {
1519 7052364 : if (a[l] > r) {
1520 7052364 : r = a[l];
1521 7052364 : last = i1;
1522 : }
1523 : }
1524 :
1525 2350788 : Array2<Real64>::size_type lind(a.index(nind, 0));
1526 2350788 : Array2<Real64>::size_type ldep(a.index(ndep, 0));
1527 2350788 : if ((a.size2() == 1u) || (x1 <= a[lind + first])) { // [ lind + first ] == ( nind, first )
1528 0 : return a[ldep + first]; // [ ldep + first ] == ( ndep, first )
1529 2350788 : } else if (x1 >= a[lind + last]) { // [ lind + last ] == ( nind, last )
1530 0 : return a[ldep + last]; // [ ldep + last ] == ( ndep, last )
1531 : } else {
1532 : int i;
1533 2350788 : int i1(first);
1534 2350788 : int i2(last);
1535 6179781 : while ((i2 - i1) > 1) {
1536 3828993 : i = i1 + ((i2 - i1) >> 1); // Tuned bit shift replaces / 2
1537 3828993 : if (x1 < a[lind + i]) { // [ lind + i ] == ( nind, i )
1538 1010936 : i2 = i;
1539 : } else {
1540 2818057 : i1 = i;
1541 : }
1542 : }
1543 2350788 : i = i2;
1544 2350788 : lind += i;
1545 2350788 : ldep += i;
1546 2350788 : Real64 const fract((x1 - a[lind - 1]) / (a[lind] - a[lind - 1])); // [ lind ] == ( nind, i ), [ lind - 1 ] == ( nind, i - 1 )
1547 2350788 : return a[ldep - 1] + fract * (a[ldep] - a[ldep - 1]); // [ ldep ] == ( ndep, i ), [ ldep - 1 ] == ( ndep, i - 1 )
1548 : }
1549 : }
1550 :
1551 35044767 : void ExteriorBCEqns(EnergyPlusData &state,
1552 : int const Delt, // Time Increment
1553 : int const i, // Node Index
1554 : int const Lay, // Layer Number for Construction
1555 : int const Surf, // Surface number
1556 : [[maybe_unused]] Array1D<Real64> const &T, // Old node Temperature in MFD finite difference solution
1557 : Array1D<Real64> &TT, // New node Temperature in MFD finite difference solution.
1558 : [[maybe_unused]] Array1D<Real64> const &Rhov, // MFD Nodal Vapor Density[kg/m3] and is the old or last time step result.
1559 : Array1D<Real64> &RhoT, // MFD vapor density for the new time step.
1560 : [[maybe_unused]] Array1D<Real64> &RH, // Nodal relative humidity
1561 : Array1D<Real64> const &TD, // The old dry Temperature at each node for the CondFD algorithm..
1562 : Array1D<Real64> &TDT, // The current or new Temperature at each node location for the CondFD solution..
1563 : Array1D<Real64> &EnthOld, // Old Nodal enthalpy
1564 : Array1D<Real64> &EnthNew, // New Nodal enthalpy
1565 : int const TotNodes, // Total nodes in layer
1566 : Real64 const HMovInsul // Conductance of movable(transparent) insulation.
1567 : )
1568 : {
1569 :
1570 : // SUBROUTINE INFORMATION:
1571 : // AUTHOR Richard Liesen
1572 : // DATE WRITTEN November, 2003
1573 : // MODIFIED B. Griffith 2010, fix adiabatic and other side surfaces
1574 : // May 2011, B. Griffith, P. Tabares
1575 : // November 2011 P. Tabares fixed problems with adiabatic walls/massless walls
1576 : // November 2011 P. Tabares fixed problems PCM stability problems
1577 : // RE-ENGINEERED Curtis Pedersen 2006
1578 :
1579 35044767 : auto &s_hbfd = state.dataHeatBalFiniteDiffMgr;
1580 35044767 : auto const &surface(state.dataSurface->Surface(Surf));
1581 35044767 : int const surface_ExtBoundCond(surface.ExtBoundCond);
1582 :
1583 : Real64 Tsky;
1584 : Real64 QRadSWOutFD; // Short wave radiation absorbed on outside of opaque surface
1585 35044767 : Real64 QRadSWOutMvInsulFD(0.0); // SW radiation at outside of Movable Insulation
1586 35044767 : if (surface_ExtBoundCond == DataSurfaces::OtherSideCondModeledExt) {
1587 : // CR8046 switch modeled rad temp for sky temp.
1588 0 : Tsky = state.dataSurface->OSCM(surface.OSCMPtr).TRad;
1589 0 : QRadSWOutFD = 0.0; // eliminate incident shortwave on underlying surface
1590 : } else { // Set the external conditions to local variables
1591 35044767 : QRadSWOutFD = state.dataHeatBalSurf->SurfOpaqQRadSWOutAbs(Surf);
1592 35044767 : QRadSWOutMvInsulFD = state.dataHeatBalSurf->SurfQRadSWOutMvIns(Surf);
1593 35044767 : Tsky = state.dataEnvrn->SkyTemp;
1594 : }
1595 :
1596 35044767 : if (surface_ExtBoundCond == DataSurfaces::Ground || state.dataEnvrn->IsRain) {
1597 5858343 : TDT(i) = TT(i) = state.dataMstBal->TempOutsideAirFD(Surf);
1598 5858343 : RhoT(i) = state.dataMstBal->RhoVaporAirOut(Surf);
1599 5858343 : s_hbfd->SurfaceFD(Surf).CpDelXRhoS1(i) = 0.0; // Outside face does not have an outer half node
1600 5858343 : s_hbfd->SurfaceFD(Surf).CpDelXRhoS2(i) = -1.0; // Set this to -1 as a flag, then set to node 2's outer half node heat capacity
1601 29186424 : } else if (surface_ExtBoundCond > 0) {
1602 : // this is actually the inside face of another surface, or maybe this same surface if adiabatic
1603 : // switch around arguments for the other surf and call routines as for interior side BC from opposite face
1604 :
1605 11812751 : int const ext_bound_construction(state.dataSurface->Surface(surface_ExtBoundCond).Construction);
1606 11812751 : int const LayIn(state.dataConstruction->Construct(ext_bound_construction).TotLayers); // layer number for call to interior eqs
1607 11812751 : int const NodeIn(s_hbfd->ConstructFD(ext_bound_construction).TotNodes + 1); // node number "I" for call to interior eqs
1608 11812751 : int const TotNodesPlusOne(TotNodes + 1);
1609 11812751 : if (surface_ExtBoundCond == Surf) { // adiabatic surface, PT added since it is not the same as interzone wall
1610 : // as Outside Boundary Condition Object can be left blank.
1611 :
1612 916581 : auto &surfaceFD = s_hbfd->SurfaceFD(Surf);
1613 916581 : InteriorBCEqns(state,
1614 : Delt,
1615 : NodeIn,
1616 : LayIn,
1617 : Surf,
1618 916581 : surfaceFD.T,
1619 916581 : surfaceFD.TT,
1620 916581 : surfaceFD.Rhov,
1621 916581 : surfaceFD.RhoT,
1622 916581 : surfaceFD.RH,
1623 916581 : surfaceFD.TD,
1624 916581 : surfaceFD.TDT,
1625 916581 : surfaceFD.EnthOld,
1626 916581 : surfaceFD.EnthNew,
1627 916581 : surfaceFD.TDreport);
1628 916581 : TDT(i) = surfaceFD.TDT(TotNodesPlusOne);
1629 916581 : TT(i) = surfaceFD.TT(TotNodesPlusOne);
1630 916581 : RhoT(i) = surfaceFD.RhoT(TotNodesPlusOne);
1631 :
1632 916581 : surfaceFD.CpDelXRhoS1(i) = 0.0; // Outside face does not have an outer half node
1633 916581 : surfaceFD.CpDelXRhoS2(i) = surfaceFD.CpDelXRhoS1(TotNodesPlusOne); // Save this for computing node flux values
1634 :
1635 : } else {
1636 :
1637 : // potential-lkl-from old CALL InteriorBCEqns(Delt,nodeIn,LayIn,Surf,SurfaceFD(Surface(Surf)%ExtBoundCond)%T, &
1638 10896170 : auto &surfaceFDEBC = s_hbfd->SurfaceFD(surface_ExtBoundCond);
1639 10896170 : InteriorBCEqns(state,
1640 : Delt,
1641 : NodeIn,
1642 : LayIn,
1643 : surface_ExtBoundCond,
1644 10896170 : surfaceFDEBC.T,
1645 10896170 : surfaceFDEBC.TT,
1646 10896170 : surfaceFDEBC.Rhov,
1647 10896170 : surfaceFDEBC.RhoT,
1648 10896170 : surfaceFDEBC.RH,
1649 10896170 : surfaceFDEBC.TD,
1650 10896170 : surfaceFDEBC.TDT,
1651 10896170 : surfaceFDEBC.EnthOld,
1652 10896170 : surfaceFDEBC.EnthNew,
1653 10896170 : surfaceFDEBC.TDreport);
1654 :
1655 10896170 : TDT(i) = surfaceFDEBC.TDT(TotNodesPlusOne);
1656 10896170 : TT(i) = surfaceFDEBC.TT(TotNodesPlusOne);
1657 10896170 : RhoT(i) = surfaceFDEBC.RhoT(TotNodesPlusOne);
1658 :
1659 10896170 : s_hbfd->SurfaceFD(Surf).CpDelXRhoS1(i) = 0.0; // Outside face does not have an outer half node
1660 10896170 : s_hbfd->SurfaceFD(Surf).CpDelXRhoS2(i) = surfaceFDEBC.CpDelXRhoS1(TotNodesPlusOne); // Save this for computing node flux values
1661 : }
1662 :
1663 11812751 : Real64 const QNetSurfFromOutside(state.dataHeatBalSurf->SurfOpaqInsFaceCondFlux(surface_ExtBoundCond)); // filled in InteriorBCEqns
1664 : // QFluxOutsideToOutSurf(Surf) = QnetSurfFromOutside
1665 11812751 : state.dataHeatBalSurf->SurfOpaqOutFaceCondFlux(Surf) = -QNetSurfFromOutside;
1666 11812751 : s_hbfd->QHeatOutFlux(Surf) = QNetSurfFromOutside;
1667 :
1668 : } else { // regular outside conditions
1669 17373673 : Real64 TDT_i(TDT(i));
1670 17373673 : Real64 const TDT_p(TDT(i + 1));
1671 :
1672 17373673 : Real64 Tgndsurface = 0.0;
1673 17373673 : if (state.dataSurface->Surface(Surf).UseSurfPropertyGndSurfTemp) {
1674 0 : Tgndsurface = state.dataSurface->GroundSurfsProperty(Surf).SurfsTempAvg;
1675 : } else {
1676 17373673 : Tgndsurface = state.dataMstBal->TempOutsideAirFD(Surf);
1677 : }
1678 :
1679 : // Boundary Conditions from Simulation for Exterior
1680 17373673 : Real64 const hconvo(state.dataMstBal->HConvExtFD(Surf));
1681 :
1682 17373673 : Real64 const hrad(state.dataMstBal->HAirFD(Surf));
1683 17373673 : Real64 const hsky(state.dataMstBal->HSkyFD(Surf));
1684 17373673 : Real64 const hgnd(state.dataMstBal->HGrndFD(Surf));
1685 17373673 : Real64 const Toa(state.dataMstBal->TempOutsideAirFD(Surf));
1686 17373673 : Real64 const Tgnd(Tgndsurface);
1687 :
1688 17373673 : if (surface.HeatTransferAlgorithm == DataSurfaces::HeatTransferModel::CondFD) {
1689 :
1690 17373673 : int const ConstrNum(surface.Construction);
1691 17373673 : int const MatLay(state.dataConstruction->Construct(ConstrNum).LayerPoint(Lay));
1692 17373673 : auto *mat = state.dataMaterial->materials(MatLay);
1693 17373673 : auto const &matFD(s_hbfd->MaterialFD(MatLay));
1694 17373673 : auto const &condActuator = s_hbfd->SurfaceFD(Surf).condMaterialActuators(Lay);
1695 17373673 : auto const &specHeatActuator = s_hbfd->SurfaceFD(Surf).specHeatMaterialActuators(Lay);
1696 :
1697 : // regular outside conditions
1698 :
1699 : // Calculate the Dry Heat Conduction Equation
1700 :
1701 17373673 : if (mat->ROnly || mat->group == Material::Group::AirGap) { // R Layer or Air Layer **********
1702 : // Use algebraic equation for TDT based on R
1703 1116311 : Real64 const Rlayer(mat->Resistance);
1704 1116311 : TDT_i = (TDT_p + (QRadSWOutFD + hgnd * Tgnd + (hconvo + hrad) * Toa + hsky * Tsky) * Rlayer) /
1705 1116311 : (1.0 + (hconvo + hgnd + hrad + hsky) * Rlayer);
1706 :
1707 1116311 : } else { // Regular or phase change material layer
1708 :
1709 : // Set Thermal Conductivity. Can be constant, simple linear temp dep or multiple linear segment temp function dep.
1710 16257362 : auto const &matFD_TempCond(matFD.TempCond);
1711 16257362 : assert(matFD_TempCond.u2() >= 3);
1712 16257362 : Real64 const lTC(matFD_TempCond.index(2, 1));
1713 : Real64 kt;
1714 16257362 : if (matFD_TempCond[lTC] + matFD_TempCond[lTC + 1] + matFD_TempCond[lTC + 2] >= 0.0) { // Multiple Linear Segment Function
1715 : // Use average temp of surface and first node for k
1716 75304 : kt = terpld(matFD_TempCond, (TDT_i + TDT_p) / 2.0, 1, 2); // 1: Temperature, 2: Thermal conductivity
1717 : } else {
1718 16182058 : kt = mat->Conductivity; // 20C base conductivity
1719 16182058 : Real64 const kt1(matFD.tk1); // linear coefficient (normally zero)
1720 16182058 : if (kt1 != 0.0) {
1721 0 : kt = +kt1 * ((TDT_i + TDT_p) / 2.0 - 20.0);
1722 : }
1723 : }
1724 :
1725 : // Check for phase change material
1726 16257362 : Real64 const TD_i(TD(i));
1727 16257362 : Real64 const Cpo(mat->SpecHeat); // Specific heat from idf
1728 16257362 : Real64 Cp(Cpo); // Specific heat modified if PCM, otherwise equal to Cpo // Will be changed if PCM
1729 16257362 : auto const &matFD_TempEnth(matFD.TempEnth);
1730 16257362 : assert(matFD_TempEnth.u2() >= 3);
1731 16257362 : Real64 const lTE(matFD_TempEnth.index(2, 1));
1732 16257362 : Real64 RhoS(mat->Density);
1733 16257362 : if (mat->hasPCM) {
1734 0 : auto *matPC = dynamic_cast<Material::MaterialPhaseChange *>(mat);
1735 0 : adjustPropertiesForPhaseChange(state, i, Surf, matPC, TD_i, TDT_i, Cp, RhoS, kt);
1736 0 : s_hbfd->SurfaceFD(Surf).EnthalpyF = matPC->enthalpyF;
1737 0 : s_hbfd->SurfaceFD(Surf).EnthalpyM = matPC->enthalpyM;
1738 16257362 : } else if (matFD_TempEnth[lTE] + matFD_TempEnth[lTE + 1] + matFD_TempEnth[lTE + 2] >=
1739 : 0.0) { // Phase change material: Use TempEnth data to generate Cp
1740 : // Enthalpy function used to get average specific heat. Updated by GS so enthalpy function is followed.
1741 0 : EnthOld(i) = terpld(matFD_TempEnth, TD_i, 1, 2); // 1: Temperature, 2: Enthalpy
1742 0 : EnthNew(i) = terpld(matFD_TempEnth, TDT_i, 1, 2); // 1: Temperature, 2: Enthalpy
1743 0 : if (EnthNew(i) != EnthOld(i)) {
1744 0 : Cp = max(Cpo, (EnthNew(i) - EnthOld(i)) / (TDT_i - TD_i));
1745 : }
1746 : } // Phase Change Material option
1747 :
1748 : // EMS Conductivity Override
1749 16257362 : if (condActuator.isActuated) {
1750 0 : kt = condActuator.actuatedValue;
1751 : }
1752 :
1753 : // EMS Specific Heat Override
1754 16257362 : if (specHeatActuator.isActuated) {
1755 0 : Cp = specHeatActuator.actuatedValue;
1756 : }
1757 :
1758 : // Update EMS internal variables
1759 16257362 : s_hbfd->SurfaceFD(Surf).condNodeReport(i) = kt;
1760 16257362 : s_hbfd->SurfaceFD(Surf).specHeatNodeReport(i) = Cp;
1761 :
1762 : // Choose Regular or Transparent Insulation Case
1763 16257362 : Real64 const DelX = s_hbfd->ConstructFD(ConstrNum).DelX(Lay);
1764 16257362 : Real64 const Delt_DelX(Delt * DelX);
1765 16257362 : s_hbfd->SurfaceFD(Surf).CpDelXRhoS1(i) = 0.0; // Outside face does not have an outer half node
1766 16257362 : s_hbfd->SurfaceFD(Surf).CpDelXRhoS2(i) = (Cp * DelX * RhoS) / 2.0; // Save this for computing node flux values
1767 :
1768 16257362 : if (HMovInsul <= 0.0) { // Regular case
1769 :
1770 16257362 : if (s_hbfd->CondFDSchemeType == CondFDScheme::CrankNicholsonSecondOrder) { // Second Order equation
1771 0 : Real64 const Cp_DelX_RhoS_2Delt(Cp * DelX * RhoS / (2.0 * Delt));
1772 0 : Real64 const kt_2DelX(kt / (2.0 * DelX));
1773 0 : Real64 const hsum(0.5 * (hconvo + hgnd + hrad + hsky));
1774 0 : TDT_i = (QRadSWOutFD + Cp_DelX_RhoS_2Delt * TD_i + kt_2DelX * (TDT_p - TD_i + TD(i + 1)) + hgnd * Tgnd +
1775 0 : (hconvo + hrad) * Toa + hsky * Tsky - hsum * TD_i) /
1776 0 : (hsum + kt_2DelX + Cp_DelX_RhoS_2Delt);
1777 16257362 : } else if (s_hbfd->CondFDSchemeType == CondFDScheme::FullyImplicitFirstOrder) { // First Order
1778 16257362 : Real64 const Two_Delt_DelX(2.0 * Delt_DelX);
1779 16257362 : Real64 const Cp_DelX2_RhoS(Cp * pow_2(DelX) * RhoS);
1780 16257362 : Real64 const Two_Delt_kt(2.0 * Delt * kt);
1781 16257362 : TDT_i = (Two_Delt_DelX * (QRadSWOutFD + hgnd * Tgnd + (hconvo + hrad) * Toa + hsky * Tsky) + Cp_DelX2_RhoS * TD_i +
1782 16257362 : Two_Delt_kt * TDT_p) /
1783 16257362 : (Two_Delt_DelX * (hconvo + hgnd + hrad + hsky) + Two_Delt_kt + Cp_DelX2_RhoS);
1784 : }
1785 :
1786 : } else { // HMovInsul > 0.0: Transparent insulation on outside
1787 : // Transparent insulation additions
1788 :
1789 : // Movable Insulation Layer Outside surface temp
1790 :
1791 0 : Real64 const TInsulOut((QRadSWOutMvInsulFD + hgnd * Tgnd + HMovInsul * TDT_i + (hconvo + hrad) * Toa + hsky * Tsky) /
1792 0 : (hconvo + hgnd + HMovInsul + hrad + hsky)); // Temperature of outside face of Outside Insulation
1793 0 : Real64 const Two_Delt_DelX(2.0 * Delt_DelX);
1794 0 : Real64 const Cp_DelX2_RhoS(Cp * pow_2(DelX) * RhoS);
1795 0 : Real64 const Two_Delt_kt(2.0 * Delt * kt);
1796 :
1797 : // Wall first node temperature behind Movable insulation
1798 0 : if (s_hbfd->CondFDSchemeType == CondFDScheme::CrankNicholsonSecondOrder) {
1799 0 : TDT_i = (Two_Delt_DelX * (QRadSWOutFD + HMovInsul * TInsulOut) + Cp_DelX2_RhoS * TD_i + Two_Delt_kt * TDT_p) /
1800 0 : (Two_Delt_DelX * HMovInsul + Two_Delt_kt + Cp_DelX2_RhoS);
1801 0 : } else if (s_hbfd->CondFDSchemeType == CondFDScheme::FullyImplicitFirstOrder) {
1802 : // Currently same as Crank Nicholson, need fully implicit formulation
1803 0 : TDT_i = (Two_Delt_DelX * (QRadSWOutFD + HMovInsul * TInsulOut) + Cp_DelX2_RhoS * TD_i + Two_Delt_kt * TDT_p) /
1804 0 : (Two_Delt_DelX * HMovInsul + Two_Delt_kt + Cp_DelX2_RhoS);
1805 : } else {
1806 0 : assert(false); // Illegal CondFDSchemeType
1807 : }
1808 :
1809 : } // Regular layer or Movable insulation cases
1810 :
1811 : } // R layer or Regular layer
1812 :
1813 17373673 : CheckFDNodeTempLimits(state, Surf, i, TDT_i);
1814 :
1815 17373673 : TDT(i) = TDT_i;
1816 :
1817 : } // regular detailed FD part or SigmaR SigmaC part
1818 :
1819 : // Determine net heat flux to outside face
1820 : // One formulation that works for Fully Implicit and CrankNicholson and massless wall
1821 :
1822 17373673 : Real64 const Toa_TDT_i(Toa - TDT_i);
1823 17373673 : Real64 const QNetSurfFromOutside(QRadSWOutFD + (hgnd * (-TDT_i + Tgnd) + (hconvo + hrad) * Toa_TDT_i + hsky * (-TDT_i + Tsky)));
1824 :
1825 : // Same sign convention as CTFs
1826 17373673 : state.dataHeatBalSurf->SurfOpaqOutFaceCondFlux(Surf) = -QNetSurfFromOutside;
1827 :
1828 : // Report all outside BC heat fluxes
1829 17373673 : state.dataHeatBalSurf->SurfQdotRadOutRepPerArea(Surf) = -(hgnd * (TDT_i - Tgnd) + hrad * (-Toa_TDT_i) + hsky * (TDT_i - Tsky));
1830 17373673 : state.dataHeatBalSurf->SurfQdotRadOutRep(Surf) = surface.Area * state.dataHeatBalSurf->SurfQdotRadOutRepPerArea(Surf);
1831 17373673 : state.dataHeatBalSurf->SurfQRadOutReport(Surf) = state.dataHeatBalSurf->SurfQdotRadOutRep(Surf) * state.dataGlobal->TimeStepZoneSec;
1832 :
1833 : } // regular BC part of the ground and Rain check
1834 35044767 : }
1835 :
1836 207129160 : void InteriorNodeEqns(EnergyPlusData &state,
1837 : int const Delt, // Time Increment
1838 : int const i, // Node Index
1839 : int const Lay, // Layer Number for Construction
1840 : int const Surf, // Surface number
1841 : [[maybe_unused]] Array1D<Real64> const &T, // INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
1842 : [[maybe_unused]] Array1D<Real64> &TT, // INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
1843 : [[maybe_unused]] Array1D<Real64> const &Rhov, // INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
1844 : [[maybe_unused]] Array1D<Real64> &RhoT, // INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
1845 : [[maybe_unused]] Array1D<Real64> &RH, // INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
1846 : Array1D<Real64> const &TD, // INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
1847 : Array1D<Real64> &TDT, // INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
1848 : Array1D<Real64> &EnthOld, // Old Nodal enthalpy
1849 : Array1D<Real64> &EnthNew // New Nodal enthalpy
1850 : )
1851 : {
1852 :
1853 : // SUBROUTINE INFORMATION:
1854 : // AUTHOR Richard Liesen
1855 : // DATE WRITTEN November, 2003
1856 : // MODIFIED May 2011, B. Griffith and P. Tabares
1857 : // RE-ENGINEERED C. O. Pedersen, 2006
1858 :
1859 207129160 : int const ConstrNum(state.dataSurface->Surface(Surf).Construction);
1860 :
1861 207129160 : auto &s_hbfd = state.dataHeatBalFiniteDiffMgr;
1862 :
1863 207129160 : int const MatLay(state.dataConstruction->Construct(ConstrNum).LayerPoint(Lay));
1864 207129160 : auto *mat = state.dataMaterial->materials(MatLay);
1865 207129160 : auto const &matFD = s_hbfd->MaterialFD(MatLay);
1866 207129160 : auto const &condActuator = s_hbfd->SurfaceFD(Surf).condMaterialActuators(Lay);
1867 207129160 : auto const &specHeatActuator = s_hbfd->SurfaceFD(Surf).specHeatMaterialActuators(Lay);
1868 :
1869 207129160 : Real64 const TD_i(TD(i));
1870 :
1871 207129160 : Real64 const TDT_m(TDT(i - 1));
1872 207129160 : Real64 TDT_i(TDT(i));
1873 207129160 : Real64 const TDT_p(TDT(i + 1));
1874 207129160 : Real64 const TDT_mi((TDT_m + TDT_i) / 2.0);
1875 207129160 : Real64 const TDT_ip((TDT_i + TDT_p) / 2.0);
1876 :
1877 : // Set Thermal Conductivity. Can be constant, simple linear temp dep or multiple linear segment temp function dep.
1878 207129160 : auto const &matFD_TempCond(matFD.TempCond);
1879 207129160 : assert(matFD_TempCond.u2() >= 3);
1880 207129160 : Real64 const lTC(matFD_TempCond.index(2, 1));
1881 : Real64 ktA1; // Variable Outer Thermal conductivity in temperature equation
1882 : Real64 ktA2; // Thermal Inner conductivity in temperature equation
1883 207129160 : if (matFD_TempCond[lTC] + matFD_TempCond[lTC + 1] + matFD_TempCond[lTC + 2] >= 0.0) { // Multiple Linear Segment Function
1884 225912 : ktA1 = terpld(matFD.TempCond, TDT_ip, 1, 2); // 1: Temperature, 2: Thermal conductivity
1885 225912 : ktA2 = terpld(matFD.TempCond, TDT_mi, 1, 2); // 1: Temperature, 2: Thermal conductivity
1886 : } else {
1887 206903248 : ktA1 = ktA2 = mat->Conductivity; // 20C base conductivity
1888 206903248 : Real64 const kt1(matFD.tk1); // temperature coefficient for simple temp dep k. // linear coefficient (normally zero)
1889 206903248 : if (kt1 != 0.0) {
1890 0 : ktA1 += kt1 * (TDT_ip - 20.0);
1891 0 : ktA2 += kt1 * (TDT_mi - 20.0);
1892 : }
1893 : }
1894 :
1895 207129160 : Real64 const Cpo(mat->SpecHeat); // Const Cp from input
1896 207129160 : Real64 Cp(Cpo); // Cp used // Will be changed if PCM
1897 207129160 : Real64 kt(0.0);
1898 207129160 : auto const &matFD_TempEnth(matFD.TempEnth);
1899 207129160 : assert(matFD_TempEnth.u2() >= 3);
1900 207129160 : Real64 const lTE(matFD_TempEnth.index(2, 1));
1901 207129160 : Real64 RhoS(mat->Density);
1902 207129160 : if (mat->hasPCM) {
1903 329910 : auto *matPC = dynamic_cast<Material::MaterialPhaseChange *>(mat);
1904 329910 : adjustPropertiesForPhaseChange(state, i, Surf, matPC, TD_i, TDT_i, Cp, RhoS, kt);
1905 329910 : ktA1 = matPC->getConductivity(TDT_ip);
1906 329910 : ktA2 = matPC->getConductivity(TDT_mi);
1907 206799250 : } else if (matFD_TempEnth[lTE] + matFD_TempEnth[lTE + 1] + matFD_TempEnth[lTE + 2] >= 0.0) { // Phase change material: Use TempEnth data
1908 0 : EnthOld(i) = terpld(matFD_TempEnth, TD_i, 1, 2); // 1: Temperature, 2: Enthalpy
1909 0 : EnthNew(i) = terpld(matFD_TempEnth, TDT_i, 1, 2); // 1: Temperature, 2: Enthalpy
1910 0 : if (EnthNew(i) != EnthOld(i)) {
1911 0 : Cp = max(Cpo, (EnthNew(i) - EnthOld(i)) / (TDT_i - TD_i));
1912 : }
1913 : } // Phase Change case
1914 :
1915 : // EMS Conductivity Override
1916 207129160 : if (condActuator.isActuated) {
1917 0 : kt = condActuator.actuatedValue;
1918 0 : ktA1 = kt;
1919 0 : ktA2 = kt;
1920 : }
1921 :
1922 : // EMS Specific Heat Override
1923 207129160 : if (specHeatActuator.isActuated) {
1924 0 : Cp = specHeatActuator.actuatedValue;
1925 : }
1926 :
1927 : // Update EMS internal variables
1928 207129160 : s_hbfd->SurfaceFD(Surf).condNodeReport(i) = kt;
1929 207129160 : s_hbfd->SurfaceFD(Surf).specHeatNodeReport(i) = Cp;
1930 :
1931 207129160 : Real64 const DelX(s_hbfd->ConstructFD(ConstrNum).DelX(Lay));
1932 207129160 : Real64 const Cp_DelX_RhoS_Delt(Cp * DelX * RhoS / Delt);
1933 :
1934 207129160 : switch (s_hbfd->CondFDSchemeType) {
1935 0 : case CondFDScheme::CrankNicholsonSecondOrder: { // Adams-Moulton second order
1936 0 : Real64 const inv2DelX(1.0 / (2.0 * DelX));
1937 0 : TDT_i = ((Cp_DelX_RhoS_Delt * TD_i) + ((ktA1 * (TD(i + 1) - TD_i + TDT_p) + ktA2 * (TD(i - 1) - TD_i + TDT_m)) * inv2DelX)) /
1938 0 : (((ktA1 + ktA2) * inv2DelX) + Cp_DelX_RhoS_Delt);
1939 0 : } break;
1940 207129160 : case CondFDScheme::FullyImplicitFirstOrder: { // Adams-Moulton First order
1941 207129160 : Real64 const invDelX(1.0 / DelX);
1942 207129160 : TDT_i = ((Cp_DelX_RhoS_Delt * TD_i) + ((ktA2 * TDT_m) + (ktA1 * TDT_p)) * invDelX) / (((ktA1 + ktA2) * invDelX) + Cp_DelX_RhoS_Delt);
1943 207129160 : } break;
1944 0 : default:
1945 0 : assert(false); // Illegal CondFDSchemeType
1946 : }
1947 :
1948 207129160 : CheckFDNodeTempLimits(state, Surf, i, TDT_i);
1949 :
1950 207129160 : TDT(i) = TDT_i;
1951 207129160 : s_hbfd->SurfaceFD(Surf).CpDelXRhoS1(i) = s_hbfd->SurfaceFD(Surf).CpDelXRhoS2(i) =
1952 207129160 : (Cp * DelX * RhoS) / 2.0; // Save this for computing node flux values, half nodes are the same here
1953 207129160 : }
1954 :
1955 83557759 : void IntInterfaceNodeEqns(EnergyPlusData &state,
1956 : int const Delt, // Time Increment
1957 : int const i, // Node Index
1958 : int const Lay, // Layer Number for Construction
1959 : int const SurfNum, // Surface number
1960 : [[maybe_unused]] Array1D<Real64> const &T, // INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
1961 : [[maybe_unused]] Array1D<Real64> &TT, // INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
1962 : [[maybe_unused]] Array1D<Real64> const &Rhov, // INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
1963 : [[maybe_unused]] Array1D<Real64> &RhoT, // INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
1964 : [[maybe_unused]] Array1D<Real64> &RH, // RELATIVE HUMIDITY.
1965 : Array1D<Real64> const &TD, // OLD NODE TEMPERATURES OF EACH HEAT TRANSFER SURF IN CONDFD.
1966 : Array1D<Real64> &TDT, // NEW NODE TEMPERATURES OF EACH HEAT TRANSFER SURF IN CONDFD.
1967 : [[maybe_unused]] Array1D<Real64> const &EnthOld, // Old Nodal enthalpy
1968 : Array1D<Real64> &EnthNew, // New Nodal enthalpy
1969 : [[maybe_unused]] int const GSiter // Iteration number of Gauss Seidel iteration
1970 : )
1971 : {
1972 :
1973 : // SUBROUTINE INFORMATION:
1974 : // AUTHOR Richard Liesen
1975 : // DATE WRITTEN November, 2003
1976 : // MODIFIED May 2011, B. Griffith, P. Tabares, add first order fully implicit, bug fixes, cleanup
1977 : // RE-ENGINEERED Curtis Pedersen, Changed to Implicit mode and included enthalpy. FY2006
1978 :
1979 : // PURPOSE OF THIS SUBROUTINE:
1980 : // calculate finite difference heat transfer for nodes that interface two different material layers inside construction
1981 :
1982 83557759 : auto &s_hbfd = state.dataHeatBalFiniteDiffMgr;
1983 83557759 : auto const &surface = state.dataSurface->Surface(SurfNum);
1984 83557759 : auto &surfFD = s_hbfd->SurfaceFD(SurfNum);
1985 :
1986 83557759 : if (surface.HeatTransferAlgorithm == DataSurfaces::HeatTransferModel::CondFD) { // HT Algo issue
1987 :
1988 83557759 : int const ConstrNum(surface.Construction);
1989 83557759 : auto const &construct(state.dataConstruction->Construct(ConstrNum));
1990 :
1991 83557759 : int const MatLay(construct.LayerPoint(Lay));
1992 83557759 : auto *mat = state.dataMaterial->materials(MatLay);
1993 :
1994 83557759 : int const MatLay2(construct.LayerPoint(Lay + 1));
1995 83557759 : auto *mat2 = state.dataMaterial->materials(MatLay2);
1996 :
1997 83557759 : auto const &condActuator1 = surfFD.condMaterialActuators(Lay);
1998 83557759 : auto const &condActuator2 = surfFD.condMaterialActuators(Lay + 1);
1999 :
2000 83557759 : auto const &specHeatActuator1 = surfFD.specHeatMaterialActuators(Lay);
2001 83557759 : auto const &specHeatActuator2 = surfFD.specHeatMaterialActuators(Lay + 1);
2002 :
2003 83557759 : auto const &heatFluxActuator = surfFD.heatSourceFluxMaterialActuators(Lay);
2004 :
2005 83557759 : Real64 const TDT_m(TDT(i - 1));
2006 83557759 : Real64 const TDT_p(TDT(i + 1));
2007 :
2008 83557759 : bool const RLayerPresent(mat->ROnly || mat->group == Material::Group::AirGap);
2009 83557759 : bool const RLayer2Present(mat2->ROnly || mat2->group == Material::Group::AirGap);
2010 :
2011 83557759 : Real64 const Rlayer(mat->Resistance); // Resistance value of R Layer
2012 83557759 : Real64 const Rlayer2(mat2->Resistance); // Resistance value of next layer to inside
2013 :
2014 83557759 : if (RLayerPresent && RLayer2Present) {
2015 :
2016 0 : TDT(i) = (Rlayer2 * TDT_m + Rlayer * TDT_p) / (Rlayer + Rlayer2); // Two adjacent R layers
2017 :
2018 : } else {
2019 :
2020 83557759 : auto const &matFD(s_hbfd->MaterialFD(MatLay));
2021 83557759 : auto const &matFD2(s_hbfd->MaterialFD(MatLay2));
2022 83557759 : Real64 TDT_i(TDT(i));
2023 :
2024 : // Set Thermal Conductivity. Can be constant, simple linear temp dep or multiple linear segment temp function dep.
2025 :
2026 83557759 : Real64 kt1(0.0);
2027 83557759 : if (!RLayerPresent) {
2028 83364011 : auto const &matFD_TempCond(matFD.TempCond);
2029 83364011 : assert(matFD_TempCond.u2() >= 3);
2030 83364011 : Real64 const lTC(matFD_TempCond.index(2, 1));
2031 83364011 : if (matFD_TempCond[lTC] + matFD_TempCond[lTC + 1] + matFD_TempCond[lTC + 2] >= 0.0) { // Multiple Linear Segment Function
2032 0 : kt1 = terpld(matFD.TempCond, (TDT_i + TDT_m) / 2.0, 1, 2); // 1: Temperature, 2: Thermal conductivity
2033 : } else {
2034 83364011 : kt1 = mat->Conductivity; // 20C base conductivity
2035 83364011 : Real64 const kt11(matFD.tk1); // temperature coefficient for simple temp dep k. // linear coefficient (normally zero)
2036 83364011 : if (kt11 != 0.0) {
2037 0 : kt1 += kt11 * ((TDT_i + TDT_m) / 2.0 - 20.0);
2038 : }
2039 : }
2040 : }
2041 :
2042 83557759 : Real64 kt2(0.0);
2043 83557759 : if (!RLayer2Present) {
2044 83364011 : auto const &matFD2_TempCond(matFD2.TempCond);
2045 83364011 : assert(matFD2_TempCond.u2() >= 3);
2046 83364011 : Real64 const lTC2(matFD2_TempCond.index(2, 1));
2047 83364011 : if (matFD2_TempCond[lTC2] + matFD2_TempCond[lTC2 + 1] + matFD2_TempCond[lTC2 + 2] >= 0.0) { // Multiple Linear Segment Function
2048 0 : kt2 = terpld(matFD2_TempCond, (TDT_i + TDT_p) / 2.0, 1, 2); // 1: Temperature, 2: Thermal conductivity
2049 : } else {
2050 83364011 : kt2 = mat2->Conductivity; // 20C base conductivity
2051 83364011 : Real64 const kt21(matFD2.tk1); // temperature coefficient for simple temp dep k. // linear coefficient (normally zero)
2052 83364011 : if (kt21 != 0.0) {
2053 0 : kt2 += kt21 * ((TDT_i + TDT_p) / 2.0 - 20.0);
2054 : }
2055 : }
2056 : }
2057 :
2058 83557759 : Real64 RhoS1(mat->Density);
2059 83557759 : Real64 const Cpo1(mat->SpecHeat); // constant Cp from input file
2060 83557759 : Real64 Cp1(Cpo1); // Will be reset if PCM
2061 83557759 : Real64 const Delx1(s_hbfd->ConstructFD(ConstrNum).DelX(Lay));
2062 :
2063 83557759 : Real64 RhoS2(mat2->Density);
2064 83557759 : Real64 const Cpo2(mat2->SpecHeat);
2065 83557759 : Real64 Cp2(Cpo2); // will be reset if PCM
2066 83557759 : Real64 const Delx2(s_hbfd->ConstructFD(ConstrNum).DelX(Lay + 1));
2067 :
2068 : // Calculate the Dry Heat Conduction Equation
2069 :
2070 : // Source/Sink Flux Capability ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2071 :
2072 83557759 : Real64 QSSFlux = 0.0;
2073 83557759 : if ((surface.Area > 0.0) && (construct.SourceSinkPresent && Lay == construct.SourceAfterLayer)) {
2074 : // Source/Sink flux value at a layer interface // Includes QPV Source
2075 5787375 : QSSFlux = (state.dataHeatBalFanSys->QRadSysSource(SurfNum) + state.dataHeatBalFanSys->QPVSysSource(SurfNum)) / surface.Area;
2076 : }
2077 :
2078 : // update report variables
2079 83557759 : auto &surfFD = s_hbfd->SurfaceFD(SurfNum);
2080 :
2081 : // only includes internal heat source
2082 83557759 : surfFD.heatSourceInternalFluxLayerReport(Lay) = QSSFlux * surface.Area;
2083 83557759 : surfFD.heatSourceInternalFluxEnergyLayerReport(Lay) = QSSFlux * surface.Area * state.dataGlobal->TimeStepZoneSec;
2084 :
2085 : // Add EMS actuated value
2086 83557759 : if (heatFluxActuator.isActuated) {
2087 160612 : Real64 actuatedVal = heatFluxActuator.actuatedValue;
2088 160612 : if (actuatedVal >= 0) {
2089 160612 : QSSFlux += heatFluxActuator.actuatedValue;
2090 : } else {
2091 0 : ShowSevereError(state, fmt::format("Surface: {}, Material: {}", surface.Name, mat->Name));
2092 0 : ShowContinueError(state, "EMS Actuator does not support negative values");
2093 0 : ShowFatalError(state, "Program terminates due to preceding conditions.");
2094 : }
2095 :
2096 : // Update report variables
2097 : // Only includes the EMS values
2098 160612 : surfFD.heatSourceEMSFluxLayerReport(Lay) = heatFluxActuator.actuatedValue * surface.Area;
2099 160612 : surfFD.heatSourceEMSFluxEnergyLayerReport(Lay) =
2100 160612 : heatFluxActuator.actuatedValue * surface.Area * state.dataGlobal->TimeStepZoneSec;
2101 : }
2102 :
2103 : //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2104 :
2105 83557759 : Real64 const TD_i(TD(i));
2106 :
2107 83557759 : auto const &matFD_TempEnth(matFD.TempEnth);
2108 83557759 : assert(matFD_TempEnth.u2() >= 3);
2109 83557759 : Real64 const lTE(matFD_TempEnth.index(2, 1));
2110 83557759 : Real64 const matFD_sum(matFD_TempEnth[lTE] + matFD_TempEnth[lTE + 1] + matFD_TempEnth[lTE + 2]);
2111 :
2112 83557759 : auto const &matFD2_TempEnth(matFD2.TempEnth);
2113 83557759 : assert(matFD2_TempEnth.u2() >= 3);
2114 83557759 : Real64 const lTE2(matFD2_TempEnth.index(2, 1));
2115 83557759 : Real64 const matFD2_sum(matFD2_TempEnth[lTE2] + matFD2_TempEnth[lTE2 + 1] + matFD2_TempEnth[lTE2 + 2]);
2116 :
2117 83557759 : if (RLayerPresent && !RLayer2Present) { // R-layer first
2118 :
2119 : // Check for PCM second layer
2120 193748 : if (mat2->hasPCM) {
2121 0 : auto *matPC2 = dynamic_cast<Material::MaterialPhaseChange *>(mat);
2122 0 : assert(matPC2 != nullptr);
2123 0 : adjustPropertiesForPhaseChange(state, i, SurfNum, matPC2, TD_i, TDT_i, Cp2, RhoS2, kt2);
2124 193748 : } else if ((matFD_sum < 0.0) && (matFD2_sum > 0.0)) { // Phase change material Layer2, Use TempEnth Data
2125 0 : Real64 const Enth2Old(terpld(matFD2_TempEnth, TD_i, 1, 2)); // 1: Temperature, 2: Thermal conductivity
2126 0 : Real64 const Enth2New(terpld(matFD2_TempEnth, TDT_i, 1, 2)); // 1: Temperature, 2: Thermal conductivity
2127 0 : EnthNew(i) = Enth2New; // This node really doesn't have an enthalpy, this gives it a value
2128 0 : if ((std::abs(Enth2New - Enth2Old) > smalldiff) && (std::abs(TDT_i - TD_i) > smalldiff)) {
2129 0 : Cp2 = max(Cpo2, (Enth2New - Enth2Old) / (TDT_i - TD_i));
2130 : }
2131 : }
2132 :
2133 : // EMS Conductivity 2 Override
2134 193748 : if (condActuator2.isActuated) {
2135 0 : kt2 = condActuator1.actuatedValue;
2136 : }
2137 :
2138 : // EMS Specific Heat 2 Override
2139 193748 : if (specHeatActuator2.isActuated) {
2140 0 : Cp2 = specHeatActuator1.actuatedValue;
2141 : }
2142 :
2143 : // Update EMS internal variables
2144 193748 : surfFD.condNodeReport(i) = kt1;
2145 193748 : surfFD.specHeatNodeReport(i) = Cp1;
2146 193748 : surfFD.condNodeReport(i + 1) = kt2;
2147 193748 : surfFD.specHeatNodeReport(i + 1) = Cp2;
2148 :
2149 : // R layer first, then PCM or regular layer
2150 193748 : Real64 const Delt_Delx2(Delt * Delx2);
2151 193748 : Real64 const Cp2_fac(Cp2 * pow_2(Delx2) * RhoS2 * Rlayer);
2152 193748 : Real64 const Delt_kt2_Rlayer(Delt * kt2 * Rlayer);
2153 193748 : if (s_hbfd->CondFDSchemeType == CondFDScheme::CrankNicholsonSecondOrder) {
2154 0 : TDT_i = (2.0 * Delt_Delx2 * QSSFlux * Rlayer + (Cp2_fac - Delt_Delx2 - Delt_kt2_Rlayer) * TD_i +
2155 0 : Delt_Delx2 * (TD(i - 1) + TDT_m) + Delt_kt2_Rlayer * (TD(i + 1) + TDT_p)) /
2156 0 : (Delt_Delx2 + Delt_kt2_Rlayer + Cp2_fac);
2157 193748 : } else if (s_hbfd->CondFDSchemeType == CondFDScheme::FullyImplicitFirstOrder) {
2158 193748 : Real64 const Two_Delt_Delx2(2.0 * Delt_Delx2);
2159 193748 : Real64 const Two_Delt_kt2_Rlayer(2.0 * Delt_kt2_Rlayer);
2160 193748 : TDT_i = (Two_Delt_Delx2 * (QSSFlux * Rlayer + TDT_m) + Cp2_fac * TD_i + Two_Delt_kt2_Rlayer * TDT_p) /
2161 193748 : (Two_Delt_Delx2 + Two_Delt_kt2_Rlayer + Cp2_fac);
2162 : }
2163 :
2164 193748 : CheckFDNodeTempLimits(state, SurfNum, i, TDT_i);
2165 :
2166 193748 : surfFD.CpDelXRhoS1(i) = 0.0; // - rlayer has no capacitance, so this is zero
2167 193748 : surfFD.CpDelXRhoS2(i) = (Cp2 * Delx2 * RhoS2) / 2.0; // Save this for computing node flux values
2168 :
2169 83557759 : } else if (!RLayerPresent && RLayer2Present) { // R-layer second
2170 :
2171 : // Check for PCM layer before R layer
2172 193748 : if (mat->hasPCM) {
2173 0 : auto *matPC = dynamic_cast<Material::MaterialPhaseChange *>(mat);
2174 0 : adjustPropertiesForPhaseChange(state, i, SurfNum, matPC, TD_i, TDT_i, Cp1, RhoS1, kt1);
2175 193748 : } else if ((matFD_sum > 0.0) && (matFD2_sum < 0.0)) { // Phase change material Layer1, Use TempEnth Data
2176 0 : Real64 const Enth1Old(terpld(matFD_TempEnth, TD_i, 1, 2)); // 1: Temperature, 2: Thermal conductivity
2177 0 : Real64 const Enth1New(terpld(matFD_TempEnth, TDT_i, 1, 2)); // 1: Temperature, 2: Thermal conductivity
2178 0 : EnthNew(i) = Enth1New; // This node really doesn't have an enthalpy, this gives it a value
2179 0 : if ((std::abs(Enth1New - Enth1Old) > smalldiff) && (std::abs(TDT_i - TD_i) > smalldiff)) {
2180 0 : Cp1 = max(Cpo1, (Enth1New - Enth1Old) / (TDT_i - TD_i));
2181 : }
2182 : }
2183 :
2184 : // EMS Conductivity 1 Override
2185 193748 : if (condActuator1.isActuated) {
2186 0 : kt1 = condActuator1.actuatedValue;
2187 : }
2188 :
2189 : // EMS Specific Heat 1 Override
2190 193748 : if (specHeatActuator1.isActuated) {
2191 0 : Cp1 = specHeatActuator1.actuatedValue;
2192 : }
2193 :
2194 : // Update EMS internal variables
2195 193748 : surfFD.condNodeReport(i) = kt1;
2196 193748 : surfFD.specHeatNodeReport(i) = Cp1;
2197 193748 : surfFD.condNodeReport(i + 1) = kt2;
2198 193748 : surfFD.specHeatNodeReport(i + 1) = Cp2;
2199 :
2200 193748 : Real64 const Delt_Delx1(Delt * Delx1);
2201 193748 : Real64 const Cp1_fac(Cp1 * pow_2(Delx1) * RhoS1 * Rlayer2);
2202 193748 : Real64 const Delt_kt1_Rlayer2(Delt * kt1 * Rlayer2);
2203 193748 : if (s_hbfd->CondFDSchemeType == CondFDScheme::CrankNicholsonSecondOrder) {
2204 0 : TDT_i = (2.0 * Delt_Delx1 * QSSFlux * Rlayer2 + (Cp1_fac - Delt_Delx1 - Delt_kt1_Rlayer2) * TD_i +
2205 0 : Delt_Delx1 * (TD(i + 1) + TDT_p) + Delt_kt1_Rlayer2 * (TD(i - 1) + TDT_m)) /
2206 0 : (Delt_Delx1 + Delt_kt1_Rlayer2 + Cp1_fac);
2207 193748 : } else if (s_hbfd->CondFDSchemeType == CondFDScheme::FullyImplicitFirstOrder) {
2208 193748 : Real64 const Two_Delt_Delx1(2.0 * Delt_Delx1);
2209 193748 : Real64 const Two_Delt_kt1_Rlayer2(2.0 * Delt_kt1_Rlayer2);
2210 193748 : TDT_i = (Two_Delt_Delx1 * (QSSFlux * Rlayer2 + TDT_p) + Cp1_fac * TD_i + Two_Delt_kt1_Rlayer2 * TDT_m) /
2211 193748 : (Two_Delt_Delx1 + Two_Delt_kt1_Rlayer2 + Cp1_fac);
2212 : }
2213 :
2214 193748 : CheckFDNodeTempLimits(state, SurfNum, i, TDT_i);
2215 :
2216 193748 : surfFD.CpDelXRhoS1(i) = (Cp1 * Delx1 * RhoS1) / 2.0; // Save this for computing node flux values
2217 193748 : surfFD.CpDelXRhoS2(i) = 0.0; // - rlayer has no capacitance, so this is zero
2218 :
2219 193748 : } else { // Regular or Phase Change on both sides of interface
2220 :
2221 : // Consider the various PCM material location cases
2222 83170263 : if ((matFD_sum > 0.0) && (matFD2_sum > 0.0)) { // Phase change material both layers, Use TempEnth Data
2223 :
2224 0 : Real64 const Enth1Old(terpld(matFD_TempEnth, TD_i, 1, 2)); // 1: Temperature, 2: Thermal conductivity
2225 0 : Real64 const Enth2Old(terpld(matFD2_TempEnth, TD_i, 1, 2)); // 1: Temperature, 2: Thermal conductivity
2226 0 : Real64 const Enth1New(terpld(matFD_TempEnth, TDT_i, 1, 2)); // 1: Temperature, 2: Thermal conductivity
2227 0 : Real64 const Enth2New(terpld(matFD2_TempEnth, TDT_i, 1, 2)); // 1: Temperature, 2: Thermal conductivity
2228 :
2229 0 : EnthNew(i) = Enth1New; // This node really doesn't have an enthalpy, this gives it a value
2230 :
2231 0 : if ((std::abs(Enth1New - Enth1Old) > smalldiff) && (std::abs(TDT_i - TD_i) > smalldiff)) {
2232 0 : Cp1 = max(Cpo1, (Enth1New - Enth1Old) / (TDT_i - TD_i));
2233 : }
2234 :
2235 0 : if ((std::abs(Enth2New - Enth2Old) > smalldiff) && (std::abs(TDT_i - TD_i) > smalldiff)) {
2236 0 : Cp2 = max(Cpo2, (Enth2New - Enth2Old) / (TDT_i - TD_i));
2237 : }
2238 :
2239 : // if
2240 :
2241 83170263 : } else if ((matFD_sum > 0.0) && (matFD2_sum < 0.0)) { // Phase change material Layer1, Use TempEnth Data
2242 :
2243 139981 : Real64 const Enth1Old(terpld(matFD_TempEnth, TD_i, 1, 2)); // 1: Temperature, 2: Thermal conductivity
2244 139981 : Real64 const Enth1New(terpld(matFD_TempEnth, TDT_i, 1, 2)); // 1: Temperature, 2: Thermal conductivity
2245 139981 : EnthNew(i) = Enth1New; // This node really doesn't have an enthalpy, this gives it a value
2246 :
2247 139981 : if ((std::abs(Enth1New - Enth1Old) > smalldiff) && (std::abs(TDT_i - TD_i) > smalldiff)) {
2248 98121 : Cp1 = max(Cpo1, (Enth1New - Enth1Old) / (TDT_i - TD_i));
2249 : }
2250 :
2251 83170263 : } else if ((matFD_sum < 0.0) && (matFD2_sum > 0.0)) { // Phase change material Layer2, Use TempEnth Data
2252 :
2253 297108 : Real64 const Enth2Old(terpld(matFD2_TempEnth, TD_i, 1, 2)); // 1: Temperature, 2: Thermal conductivity
2254 297108 : Real64 const Enth2New(terpld(matFD2_TempEnth, TDT_i, 1, 2)); // 1: Temperature, 2: Thermal conductivity
2255 297108 : EnthNew(i) = Enth2New; // This node really doesn't have an enthalpy, this gives it a value
2256 :
2257 297108 : if ((std::abs(Enth2New - Enth2Old) > smalldiff) && (std::abs(TDT_i - TD_i) > smalldiff)) {
2258 214156 : Cp2 = max(Cpo2, (Enth2New - Enth2Old) / (TDT_i - TD_i));
2259 : }
2260 :
2261 : } // Phase change material check
2262 :
2263 83170263 : if (mat->hasPCM) {
2264 0 : auto *matPC = dynamic_cast<Material::MaterialPhaseChange *>(mat);
2265 0 : adjustPropertiesForPhaseChange(state, i, SurfNum, matPC, TD_i, TDT_i, Cp1, RhoS1, kt1);
2266 : }
2267 83170263 : if (mat2->hasPCM) {
2268 0 : auto *matPC2 = dynamic_cast<Material::MaterialPhaseChange *>(mat2);
2269 0 : adjustPropertiesForPhaseChange(state, i, SurfNum, matPC2, TD_i, TDT_i, Cp2, RhoS2, kt2);
2270 : }
2271 :
2272 : // EMS Conductivity 1 Override
2273 83170263 : if (condActuator1.isActuated) {
2274 477915 : kt1 = condActuator1.actuatedValue;
2275 : }
2276 :
2277 : // EMS Conductivity 2 Override
2278 83170263 : if (condActuator2.isActuated) {
2279 477915 : kt2 = condActuator2.actuatedValue;
2280 : }
2281 :
2282 : // EMS Specific Heat 1 Override
2283 83170263 : if (specHeatActuator1.isActuated) {
2284 477915 : Cp1 = specHeatActuator1.actuatedValue;
2285 : }
2286 :
2287 : // EMS Specific Heat 2 Override
2288 83170263 : if (specHeatActuator2.isActuated) {
2289 477915 : Cp2 = specHeatActuator2.actuatedValue;
2290 : }
2291 :
2292 : // Update EMS internal variables
2293 83170263 : surfFD.condNodeReport(i) = kt1;
2294 83170263 : surfFD.specHeatNodeReport(i) = Cp1;
2295 83170263 : surfFD.condNodeReport(i + 1) = kt2;
2296 83170263 : surfFD.specHeatNodeReport(i + 1) = Cp2;
2297 :
2298 83170263 : Real64 const Delt_Delx1(Delt * Delx1);
2299 83170263 : Real64 const Delt_Delx2(Delt * Delx2);
2300 83170263 : Real64 const Delt_Delx1_kt2(Delt_Delx1 * kt2);
2301 83170263 : Real64 const Delt_Delx2_kt1(Delt_Delx2 * kt1);
2302 83170263 : Real64 const Delt_sum(Delt_Delx1_kt2 + Delt_Delx2_kt1);
2303 83170263 : Real64 const Cp1_fac(Cp1 * pow_2(Delx1) * Delx2 * RhoS1);
2304 83170263 : Real64 const Cp2_fac(Cp2 * Delx1 * pow_2(Delx2) * RhoS2);
2305 83170263 : Real64 const Cp_fac(Cp1_fac + Cp2_fac);
2306 83170263 : if (s_hbfd->CondFDSchemeType ==
2307 : CondFDScheme::CrankNicholsonSecondOrder) { // Regular Internal Interface Node with Source/sink using Adams Moulton second
2308 : // order
2309 0 : TDT_i = (2.0 * Delt_Delx1 * Delx2 * QSSFlux + (Cp_fac - Delt_sum) * TD_i + Delt_Delx1_kt2 * (TD(i + 1) + TDT_p) +
2310 0 : Delt_Delx2_kt1 * (TD(i - 1) + TDT_m)) /
2311 0 : (Delt_sum + Cp_fac);
2312 83170263 : } else if (s_hbfd->CondFDSchemeType == CondFDScheme::FullyImplicitFirstOrder) { // First order adams moulton
2313 83170263 : TDT_i = (2.0 * (Delt_Delx1 * Delx2 * QSSFlux + Delt_Delx2_kt1 * TDT_m + Delt_Delx1_kt2 * TDT_p) + Cp_fac * TD_i) /
2314 83170263 : (2.0 * (Delt_Delx2_kt1 + Delt_Delx1_kt2) + Cp_fac);
2315 : }
2316 :
2317 83170263 : CheckFDNodeTempLimits(state, SurfNum, i, TDT_i);
2318 :
2319 83170263 : surfFD.CpDelXRhoS1(i) = (Cp1 * Delx1 * RhoS1) / 2.0; // Save this for computing node flux values
2320 83170263 : surfFD.CpDelXRhoS2(i) = (Cp2 * Delx2 * RhoS2) / 2.0; // Save this for computing node flux values
2321 :
2322 83170263 : if (construct.SourceSinkPresent && (Lay == construct.SourceAfterLayer)) {
2323 5787375 : state.dataHeatBalFanSys->TCondFDSourceNode(SurfNum) = TDT_i; // Transfer node temp to Radiant System
2324 5787375 : state.dataHeatBalSurf->SurfTempSource(SurfNum) = TDT_i; // Transfer node temp to DataHeatBalSurface module
2325 5787375 : surfFD.QSource = QSSFlux;
2326 5787375 : surfFD.SourceNodeNum = i;
2327 : }
2328 :
2329 83170263 : if (construct.SourceSinkPresent && (Lay == construct.TempAfterLayer)) {
2330 5787375 : state.dataHeatBalSurf->SurfTempUserLoc(SurfNum) = TDT_i; // Transfer node temp to DataHeatBalSurface module
2331 : }
2332 :
2333 : } // End of R-layer and Regular check
2334 :
2335 83557759 : TDT(i) = TDT_i;
2336 : }
2337 :
2338 : } // End of the CondFD if block
2339 83557759 : }
2340 :
2341 46857518 : void InteriorBCEqns(EnergyPlusData &state,
2342 : int const Delt, // Time Increment
2343 : int const i, // Node Index
2344 : int const Lay, // Layer Number for Construction
2345 : int const SurfNum, // Surface number
2346 : [[maybe_unused]] Array1D<Real64> const &T, // INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF (Old).
2347 : [[maybe_unused]] Array1D<Real64> &TT, // INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF (New).
2348 : [[maybe_unused]] Array1D<Real64> const &Rhov, // INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
2349 : [[maybe_unused]] Array1D<Real64> &RhoT, // INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
2350 : [[maybe_unused]] Array1D<Real64> &RH, // INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
2351 : Array1D<Real64> const &TD, // INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
2352 : Array1D<Real64> &TDT, // INSIDE SURFACE TEMPERATURE OF EACH HEAT TRANSFER SURF.
2353 : Array1D<Real64> &EnthOld, // Old Nodal enthalpy
2354 : Array1D<Real64> &EnthNew, // New Nodal enthalpy
2355 : Array1D<Real64> &TDreport // Temperature value from previous HeatSurfaceHeatManager iteration's value
2356 : )
2357 : {
2358 : // SUBROUTINE INFORMATION:
2359 : // AUTHOR Richard Liesen
2360 : // DATE WRITTEN November, 2003
2361 : // MODIFIED B. Griffith, P. Tabares, May 2011, add first order fully implicit, bug fixes, cleanup
2362 : // November 2011 P. Tabares fixed problems with adiabatic walls/massless walls
2363 : // November 2011 P. Tabares fixed problems PCM stability problems
2364 : // RE-ENGINEERED C. O. Pedersen 2006
2365 :
2366 : // PURPOSE OF THIS SUBROUTINE:
2367 : // Calculate the heat transfer at the node on the surfaces inside face (facing zone)
2368 :
2369 46857518 : auto &s_hbfd = state.dataHeatBalFiniteDiffMgr;
2370 :
2371 46857518 : auto const &surface(state.dataSurface->Surface(SurfNum));
2372 46857518 : auto &surfFD = s_hbfd->SurfaceFD(SurfNum);
2373 46857518 : int const ConstrNum(surface.Construction);
2374 :
2375 : // Set the internal conditions to local variables
2376 : Real64 const NetLWRadToSurfFD(
2377 46857518 : state.dataHeatBalSurf->SurfQdotRadNetLWInPerArea(SurfNum)); // Net interior long wavelength radiation to surface from other surfaces
2378 46857518 : Real64 const QRadSWInFD(state.dataHeatBalSurf->SurfOpaqQRadSWInAbs(SurfNum)); // Short wave radiation absorbed on inside of opaque surface
2379 : Real64 const SurfQdotRadHVACInPerAreaFD(
2380 46857518 : state.dataHeatBalSurf->SurfQdotRadHVACInPerArea(SurfNum)); // Total current radiant heat flux at a surface
2381 46857518 : Real64 const QRadThermInFD(state.dataHeatBal->SurfQdotRadIntGainsInPerArea(SurfNum)); // Thermal radiation absorbed on inside surfaces
2382 :
2383 : // Boundary Conditions from Simulation for Interior
2384 46857518 : Real64 hconvi(state.dataMstBal->HConvInFD(SurfNum));
2385 :
2386 46857518 : Real64 const Tia(state.dataZoneTempPredictorCorrector->zoneHeatBalance(surface.Zone).MAT);
2387 :
2388 : //++++++++++++++++++++++++++++++++++++++++++++++++++++++
2389 : // Do all the nodes in the surface Else will switch to SigmaR,SigmaC
2390 46857518 : Real64 TDT_i(TDT(i));
2391 46857518 : Real64 const QFac(NetLWRadToSurfFD + QRadSWInFD + QRadThermInFD + SurfQdotRadHVACInPerAreaFD);
2392 46857518 : if (surface.HeatTransferAlgorithm == DataSurfaces::HeatTransferModel::CondFD) {
2393 46857518 : int const MatLay(state.dataConstruction->Construct(ConstrNum).LayerPoint(Lay));
2394 46857518 : auto *mat = state.dataMaterial->materials(MatLay);
2395 46857518 : auto const &matFD(s_hbfd->MaterialFD(MatLay));
2396 46857518 : auto const &condActuator = surfFD.condMaterialActuators(Lay);
2397 46857518 : auto const &specHeatActuator = surfFD.specHeatMaterialActuators(Lay);
2398 :
2399 : // Calculate the Dry Heat Conduction Equation
2400 :
2401 46857518 : if (mat->ROnly || mat->group == Material::Group::AirGap) { // R Layer or Air Layer
2402 : // Use algebraic equation for TDT based on R
2403 1116311 : Real64 constexpr IterDampConst(
2404 : 5.0); // Damping constant for inside surface temperature iterations. Only used for massless (R-value only) Walls
2405 1116311 : Real64 const Rlayer(mat->Resistance);
2406 1116311 : if ((i == 1) && (surface.ExtBoundCond > 0)) { // this is for an adiabatic partition
2407 0 : TDT_i = (TDT(i + 1) + (QFac + hconvi * Tia + TDreport(i) * IterDampConst) * Rlayer) / (1.0 + (hconvi + IterDampConst) * Rlayer);
2408 : } else { // regular wall
2409 1116311 : TDT_i = (TDT(i - 1) + (QFac + hconvi * Tia + TDreport(i) * IterDampConst) * Rlayer) / (1.0 + (hconvi + IterDampConst) * Rlayer);
2410 : }
2411 1116311 : surfFD.CpDelXRhoS1(i) = 0.0; // Save this for computing node flux values - rlayer has no capacitance
2412 1116311 : surfFD.CpDelXRhoS2(i) = 0.0; // Inside face does not have an inner half node
2413 :
2414 1116311 : } else { // Regular or PCM
2415 45741207 : Real64 const TDT_m(TDT(i - 1));
2416 :
2417 : // Set Thermal Conductivity. Can be constant, simple linear temp dep or multiple linear segment temp function dep.
2418 45741207 : auto const &matFD_TempCond(matFD.TempCond);
2419 45741207 : assert(matFD_TempCond.u2() >= 3);
2420 45741207 : Real64 const lTC(matFD_TempCond.index(2, 1));
2421 : Real64 kt;
2422 45741207 : if (matFD_TempCond[lTC] + matFD_TempCond[lTC + 1] + matFD_TempCond[lTC + 2] >= 0.0) { // Multiple Linear Segment Function
2423 : // Use average of surface and first node temp for determining k
2424 75304 : kt = terpld(matFD_TempCond, (TDT_i + TDT_m) / 2.0, 1, 2); // 1: Temperature, 2: Thermal conductivity
2425 : } else {
2426 45665903 : kt = mat->Conductivity; // 20C base conductivity
2427 45665903 : Real64 const kt1(matFD.tk1); // linear coefficient (normally zero)
2428 45665903 : if (kt1 != 0.0) {
2429 0 : kt = +kt1 * ((TDT_i + TDT_m) / 2.0 - 20.0);
2430 : }
2431 : }
2432 :
2433 45741207 : Real64 RhoS(mat->Density);
2434 45741207 : Real64 const TD_i(TD(i));
2435 45741207 : Real64 const Cpo(mat->SpecHeat);
2436 45741207 : Real64 Cp(Cpo); // Will be changed if PCM
2437 45741207 : auto const &matFD_TempEnth(matFD.TempEnth);
2438 45741207 : assert(matFD_TempEnth.u2() >= 3);
2439 45741207 : Real64 const lTE(matFD_TempEnth.index(2, 1));
2440 45741207 : if (mat->hasPCM) {
2441 219940 : auto *matPC = dynamic_cast<Material::MaterialPhaseChange *>(mat);
2442 219940 : adjustPropertiesForPhaseChange(state, i, SurfNum, matPC, TD_i, TDT_i, Cp, RhoS, kt);
2443 45521267 : } else if (matFD_TempEnth[lTE] + matFD_TempEnth[lTE + 1] + matFD_TempEnth[lTE + 2] >=
2444 : 0.0) { // Phase change material: Use TempEnth data
2445 437089 : EnthOld(i) = terpld(matFD_TempEnth, TD_i, 1, 2); // 1: Temperature, 2: Enthalpy
2446 437089 : EnthNew(i) = terpld(matFD_TempEnth, TDT_i, 1, 2); // 1: Temperature, 2: Enthalpy
2447 437089 : if ((std::abs(EnthNew(i) - EnthOld(i)) > smalldiff) && (std::abs(TDT_i - TD_i) > smalldiff)) {
2448 331648 : Cp = max(Cpo, (EnthNew(i) - EnthOld(i)) / (TDT_i - TD_i));
2449 : }
2450 : } // Phase change material check
2451 :
2452 : // EMS Conductivity Override
2453 45741207 : if (condActuator.isActuated) {
2454 0 : kt = condActuator.actuatedValue;
2455 : }
2456 :
2457 : // EMS Specific Heat Override
2458 45741207 : if (specHeatActuator.isActuated) {
2459 0 : Cp = specHeatActuator.actuatedValue;
2460 : }
2461 :
2462 : // Update EMS internal variables
2463 45741207 : surfFD.condNodeReport(i) = kt;
2464 45741207 : surfFD.specHeatNodeReport(i) = Cp;
2465 :
2466 45741207 : Real64 const DelX(s_hbfd->ConstructFD(ConstrNum).DelX(Lay));
2467 45741207 : Real64 const Delt_DelX(Delt * DelX);
2468 45741207 : Real64 const Two_Delt_DelX(2.0 * Delt_DelX);
2469 45741207 : Real64 const Delt_kt(Delt * kt);
2470 45741207 : Real64 const Cp_DelX2_RhoS(Cp * pow_2(DelX) * RhoS);
2471 45741207 : if ((surface.ExtBoundCond > 0) && (i == 1)) { // this is for an adiabatic or interzone partition
2472 0 : if (s_hbfd->CondFDSchemeType == CondFDScheme::CrankNicholsonSecondOrder) { // Adams-Moulton second order
2473 0 : TDT_i = (Two_Delt_DelX * (QFac + hconvi * Tia) + (Cp_DelX2_RhoS - Delt_DelX * hconvi - Delt_kt) * TD_i +
2474 0 : Delt_kt * (TD(i + 1) + TDT(i + 1))) /
2475 0 : (Delt_DelX * hconvi + Delt_kt + Cp_DelX2_RhoS);
2476 0 : } else if (s_hbfd->CondFDSchemeType == CondFDScheme::FullyImplicitFirstOrder) { // Adams-Moulton First order
2477 0 : Real64 const Two_Delt_kt(2.0 * Delt_kt);
2478 0 : TDT_i = (Two_Delt_DelX * (QFac + hconvi * Tia) + Cp_DelX2_RhoS * TD_i + Two_Delt_kt * TDT(i + 1)) /
2479 0 : (Two_Delt_DelX * hconvi + Two_Delt_kt + Cp_DelX2_RhoS);
2480 : }
2481 0 : } else { // for regular or interzone walls
2482 45741207 : if (s_hbfd->CondFDSchemeType == CondFDScheme::CrankNicholsonSecondOrder) {
2483 0 : TDT_i = (Two_Delt_DelX * (QFac + hconvi * Tia) + (Cp_DelX2_RhoS - Delt_DelX * hconvi - Delt_kt) * TD_i +
2484 0 : Delt_kt * (TD(i - 1) + TDT_m)) /
2485 0 : (Delt_DelX * hconvi + Delt_kt + Cp_DelX2_RhoS);
2486 45741207 : } else if (s_hbfd->CondFDSchemeType == CondFDScheme::FullyImplicitFirstOrder) {
2487 45741207 : Real64 const Two_Delt_kt(2.0 * Delt_kt);
2488 45741207 : TDT_i = (Two_Delt_DelX * (QFac + hconvi * Tia) + Cp_DelX2_RhoS * TD_i + Two_Delt_kt * TDT_m) /
2489 45741207 : (Two_Delt_DelX * hconvi + Two_Delt_kt + Cp_DelX2_RhoS);
2490 : }
2491 : }
2492 45741207 : surfFD.CpDelXRhoS1(i) = (Cp * DelX * RhoS) / 2.0; // Save this for computing node flux values
2493 45741207 : surfFD.CpDelXRhoS2(i) = 0.0; // Inside face does not have an inner half node
2494 :
2495 : } // Regular or R layer
2496 :
2497 46857518 : CheckFDNodeTempLimits(state, SurfNum, i, TDT_i);
2498 :
2499 46857518 : TDT(i) = TDT_i;
2500 :
2501 : } // End of Regular node or SigmaR SigmaC option
2502 :
2503 46857518 : Real64 const QNetSurfInside(-(QFac + hconvi * (-TDT_i + Tia)));
2504 : // Pass inside conduction Flux [W/m2] to DataHeatBalanceSurface array
2505 46857518 : state.dataHeatBalSurf->SurfOpaqInsFaceCondFlux(SurfNum) = QNetSurfInside;
2506 46857518 : }
2507 :
2508 : // todo - function not used
2509 0 : void CheckFDSurfaceTempLimits(EnergyPlusData &state,
2510 : int const SurfNum, // surface number
2511 : Real64 const CheckTemperature // calculated temperature, not reset
2512 : )
2513 : {
2514 :
2515 : // SUBROUTINE INFORMATION:
2516 : // AUTHOR Linda Lawrie
2517 : // DATE WRITTEN August 2012
2518 :
2519 : // PURPOSE OF THIS SUBROUTINE:
2520 : // Provides a single entry point for checking surface temperature limits as well as
2521 : // setting up for recurring errors if too low or too high.
2522 :
2523 : // METHODOLOGY EMPLOYED:
2524 : // Use methodology similar to HBSurfaceManager
2525 0 : auto &s_hbfd = state.dataHeatBalFiniteDiffMgr;
2526 :
2527 0 : int ZoneNum = state.dataSurface->Surface(SurfNum).Zone;
2528 :
2529 0 : if (state.dataGlobal->WarmupFlag) {
2530 0 : ++s_hbfd->WarmupSurfTemp;
2531 : }
2532 0 : if (!state.dataGlobal->WarmupFlag || s_hbfd->WarmupSurfTemp > 10 || state.dataGlobal->DisplayExtraWarnings) {
2533 0 : if (CheckTemperature < DataHeatBalSurface::MinSurfaceTempLimit) {
2534 0 : if (state.dataSurface->SurfLowTempErrCount(SurfNum) == 0) {
2535 0 : ShowSevereMessage(state,
2536 0 : format("Temperature (low) out of bounds [{:.2R}] for zone=\"{}\", for surface=\"{}\"",
2537 : CheckTemperature,
2538 0 : state.dataHeatBal->Zone(ZoneNum).Name,
2539 0 : state.dataSurface->Surface(SurfNum).Name));
2540 0 : ShowContinueErrorTimeStamp(state, "");
2541 0 : if (!state.dataHeatBal->Zone(ZoneNum).TempOutOfBoundsReported) {
2542 0 : ShowContinueError(state, format("Zone=\"{}\", Diagnostic Details:", state.dataHeatBal->Zone(ZoneNum).Name));
2543 0 : if (state.dataHeatBal->Zone(ZoneNum).FloorArea > 0.0) {
2544 0 : ShowContinueError(
2545 : state,
2546 0 : format("...Internal Heat Gain [{:.3R}] W/m2",
2547 0 : state.dataHeatBal->Zone(ZoneNum).InternalHeatGains / state.dataHeatBal->Zone(ZoneNum).FloorArea));
2548 : } else {
2549 0 : ShowContinueError(
2550 0 : state, format("...Internal Heat Gain (no floor) [{:.3R}] W", state.dataHeatBal->Zone(ZoneNum).InternalHeatGains));
2551 : }
2552 0 : if (state.afn->simulation_control.type == AirflowNetwork::ControlType::NoMultizoneOrDistribution) {
2553 0 : ShowContinueError(state,
2554 0 : format("...Infiltration/Ventilation [{:.3R}] m3/s", state.dataHeatBal->Zone(ZoneNum).NominalInfilVent));
2555 0 : ShowContinueError(state, format("...Mixing/Cross Mixing [{:.3R}] m3/s", state.dataHeatBal->Zone(ZoneNum).NominalMixing));
2556 : } else {
2557 0 : ShowContinueError(state, "...Airflow Network Simulation: Nominal Infiltration/Ventilation/Mixing not available.");
2558 : }
2559 0 : if (state.dataHeatBal->Zone(ZoneNum).IsControlled) {
2560 0 : ShowContinueError(state, "...Zone is part of HVAC controlled system.");
2561 : } else {
2562 0 : ShowContinueError(state, "...Zone is not part of HVAC controlled system.");
2563 : }
2564 0 : state.dataHeatBal->Zone(ZoneNum).TempOutOfBoundsReported = true;
2565 : }
2566 0 : ShowRecurringSevereErrorAtEnd(state,
2567 0 : "Temperature (low) out of bounds for zone=" + state.dataHeatBal->Zone(ZoneNum).Name +
2568 0 : " for surface=" + state.dataSurface->Surface(SurfNum).Name,
2569 0 : state.dataSurface->SurfLowTempErrCount(SurfNum),
2570 : CheckTemperature,
2571 : CheckTemperature,
2572 : _,
2573 : "C",
2574 : "C");
2575 : } else {
2576 0 : ShowRecurringSevereErrorAtEnd(state,
2577 0 : "Temperature (low) out of bounds for zone=" + state.dataHeatBal->Zone(ZoneNum).Name +
2578 0 : " for surface=" + state.dataSurface->Surface(SurfNum).Name,
2579 0 : state.dataSurface->SurfLowTempErrCount(SurfNum),
2580 : CheckTemperature,
2581 : CheckTemperature,
2582 : _,
2583 : "C",
2584 : "C");
2585 : }
2586 : } else {
2587 0 : if (state.dataSurface->SurfHighTempErrCount(SurfNum) == 0) {
2588 0 : ShowSevereMessage(state,
2589 0 : format("Temperature (high) out of bounds ({:.2R}] for zone=\"{}\", for surface=\"{}\"",
2590 : CheckTemperature,
2591 0 : state.dataHeatBal->Zone(ZoneNum).Name,
2592 0 : state.dataSurface->Surface(SurfNum).Name));
2593 0 : ShowContinueErrorTimeStamp(state, "");
2594 0 : if (!state.dataHeatBal->Zone(ZoneNum).TempOutOfBoundsReported) {
2595 0 : ShowContinueError(state, format("Zone=\"{}\", Diagnostic Details:", state.dataHeatBal->Zone(ZoneNum).Name));
2596 0 : if (state.dataHeatBal->Zone(ZoneNum).FloorArea > 0.0) {
2597 0 : ShowContinueError(
2598 : state,
2599 0 : format("...Internal Heat Gain [{:.3R}] W/m2",
2600 0 : state.dataHeatBal->Zone(ZoneNum).InternalHeatGains / state.dataHeatBal->Zone(ZoneNum).FloorArea));
2601 : } else {
2602 0 : ShowContinueError(
2603 0 : state, format("...Internal Heat Gain (no floor) [{:.3R}] W", state.dataHeatBal->Zone(ZoneNum).InternalHeatGains));
2604 : }
2605 0 : if (state.afn->simulation_control.type == AirflowNetwork::ControlType::NoMultizoneOrDistribution) {
2606 0 : ShowContinueError(state,
2607 0 : format("...Infiltration/Ventilation [{:.3R}] m3/s", state.dataHeatBal->Zone(ZoneNum).NominalInfilVent));
2608 0 : ShowContinueError(state, format("...Mixing/Cross Mixing [{:.3R}] m3/s", state.dataHeatBal->Zone(ZoneNum).NominalMixing));
2609 : } else {
2610 0 : ShowContinueError(state, "...Airflow Network Simulation: Nominal Infiltration/Ventilation/Mixing not available.");
2611 : }
2612 0 : if (state.dataHeatBal->Zone(ZoneNum).IsControlled) {
2613 0 : ShowContinueError(state, "...Zone is part of HVAC controlled system.");
2614 : } else {
2615 0 : ShowContinueError(state, "...Zone is not part of HVAC controlled system.");
2616 : }
2617 0 : state.dataHeatBal->Zone(ZoneNum).TempOutOfBoundsReported = true;
2618 : }
2619 0 : ShowRecurringSevereErrorAtEnd(state,
2620 0 : "Temperature (high) out of bounds for zone=" + state.dataHeatBal->Zone(ZoneNum).Name +
2621 0 : " for surface=" + state.dataSurface->Surface(SurfNum).Name,
2622 0 : state.dataSurface->SurfHighTempErrCount(SurfNum),
2623 : CheckTemperature,
2624 : CheckTemperature,
2625 : _,
2626 : "C",
2627 : "C");
2628 : } else {
2629 0 : ShowRecurringSevereErrorAtEnd(state,
2630 0 : "Temperature (high) out of bounds for zone=" + state.dataHeatBal->Zone(ZoneNum).Name +
2631 0 : " for surface=" + state.dataSurface->Surface(SurfNum).Name,
2632 0 : state.dataSurface->SurfHighTempErrCount(SurfNum),
2633 : CheckTemperature,
2634 : CheckTemperature,
2635 : _,
2636 : "C",
2637 : "C");
2638 : }
2639 : }
2640 : }
2641 0 : }
2642 :
2643 354918110 : void CheckFDNodeTempLimits(EnergyPlusData &state,
2644 : int surfNum, // surface number
2645 : int nodeNum, // node number
2646 : Real64 &nodeTemp // calculated temperature, not reset
2647 : )
2648 : {
2649 354918110 : auto &s_hbfd = state.dataHeatBalFiniteDiffMgr;
2650 :
2651 354918110 : auto &surfFD = s_hbfd->SurfaceFD(surfNum);
2652 354918110 : auto &surfName = state.dataSurface->Surface(surfNum).Name;
2653 354918110 : auto &minTempLimit = DataHeatBalSurface::MinSurfaceTempLimit;
2654 354918110 : auto &maxTempLimit = state.dataHeatBalSurf->MaxSurfaceTempLimit;
2655 354918110 : if (nodeTemp < minTempLimit) {
2656 0 : if (surfFD.indexNodeMinTempLimit == 0) {
2657 0 : ShowSevereMessage(state,
2658 0 : format("Node temperature (low) out of bounds [{:.2R}] for surface={}, node={}", nodeTemp, surfName, nodeNum));
2659 0 : ShowContinueErrorTimeStamp(state, "");
2660 0 : ShowContinueError(state, format("Value has been reset to the lower limit value of {:.2R}.", minTempLimit));
2661 : }
2662 0 : ShowRecurringSevereErrorAtEnd(
2663 0 : state, "Node temperature (low) out of bounds for surface=" + surfName, surfFD.indexNodeMinTempLimit, nodeTemp, nodeTemp, _, "C", "C");
2664 0 : nodeTemp = minTempLimit;
2665 354918110 : } else if (nodeTemp > maxTempLimit) {
2666 0 : if (surfFD.indexNodeMaxTempLimit == 0) {
2667 0 : ShowSevereMessage(state,
2668 0 : format("Node temperature (high) out of bounds [{:.2R}] for surface={}, node={}", nodeTemp, surfName, nodeNum));
2669 0 : ShowContinueErrorTimeStamp(state, "");
2670 0 : ShowContinueError(state, format("Value has been reset to the upper limit value of {:.2R}.", maxTempLimit));
2671 : }
2672 0 : ShowRecurringSevereErrorAtEnd(state,
2673 0 : "Node temperature (high) out of bounds for surface=" + surfName,
2674 0 : surfFD.indexNodeMaxTempLimit,
2675 : nodeTemp,
2676 : nodeTemp,
2677 : _,
2678 : "C",
2679 : "C");
2680 0 : nodeTemp = maxTempLimit;
2681 : }
2682 354918110 : }
2683 :
2684 10240031 : void CalcNodeHeatFlux(EnergyPlusData &state,
2685 : int const SurfNum, // surface number
2686 : int const TotNodes // number of nodes in surface
2687 : )
2688 : {
2689 :
2690 : // SUBROUTINE INFORMATION:
2691 : // AUTHOR M.J. Witte
2692 : // DATE WRITTEN Sept-Nov 2015
2693 : // PURPOSE OF THIS SUBROUTINE:
2694 : // Calculate flux at each condFD node
2695 :
2696 10240031 : auto &s_hbfd = state.dataHeatBalFiniteDiffMgr;
2697 10240031 : auto &surfFD = s_hbfd->SurfaceFD(SurfNum);
2698 :
2699 : // SurfaceFD.QDreport( n ) is the flux at node n
2700 : // When this is called TDT( NodeNum ) is the new node temp and TDpriortimestep( NodeNum ) holds the previous node temp
2701 : // For the TDT and TDpriortimestep arrays, Node 1 is the outside face, and Node TotNodes+1 is the inside face
2702 :
2703 : // Last node is always the surface inside face. Start calculations here because the outside face is not defined for all surfaces.
2704 : // Note that TotNodes is the number of nodes in the surface including the outside face node, but not the inside face node
2705 : // so the arrays are all allocated to Totodes+1
2706 :
2707 : // Heat flux at the inside face node (TotNodes+1)
2708 10240031 : surfFD.QDreport(TotNodes + 1) = state.dataHeatBalSurf->SurfOpaqInsFaceCondFlux(SurfNum);
2709 :
2710 : // Heat flux for remaining nodes.
2711 105113543 : for (int node = TotNodes; node >= 1; --node) {
2712 : // Start with inside face (above) and work outward, positive value is flowing towards the inside face
2713 : // CpDelXRhoS1 is outer half-node heat capacity, CpDelXRhoS2 is inner half node heat capacity
2714 : Real64 interNodeFlux; // heat flux at the plane between node and node+1 [W/m2]
2715 : Real64 sourceFlux; // Internal source flux [W/m2]
2716 94873512 : if (surfFD.SourceNodeNum == node) {
2717 1184661 : sourceFlux = surfFD.QSource;
2718 : } else {
2719 93688851 : sourceFlux = 0.0;
2720 : }
2721 94873512 : interNodeFlux = surfFD.QDreport(node + 1) + surfFD.CpDelXRhoS1(node + 1) * (surfFD.TDT(node + 1) - surfFD.TDpriortimestep(node + 1)) /
2722 94873512 : state.dataGlobal->TimeStepZoneSec;
2723 94873512 : surfFD.QDreport(node) = interNodeFlux - sourceFlux +
2724 94873512 : surfFD.CpDelXRhoS2(node) * (surfFD.TDT(node) - surfFD.TDpriortimestep(node)) / state.dataGlobal->TimeStepZoneSec;
2725 : }
2726 10240031 : if (state.dataEnvrn->IsRain) {
2727 0 : state.dataHeatBalSurf->SurfOpaqOutFaceCondFlux(SurfNum) = -surfFD.QDreport(1); // Update the outside flux if it is raining
2728 : }
2729 10240031 : }
2730 :
2731 549850 : void adjustPropertiesForPhaseChange(EnergyPlusData &state,
2732 : int finiteDifferenceLayerIndex,
2733 : int surfNum,
2734 : Material::MaterialPhaseChange *mat,
2735 : Real64 temperaturePrevious,
2736 : Real64 temperatureUpdated,
2737 : Real64 &updatedSpecificHeat,
2738 : Real64 &updatedDensity,
2739 : Real64 &updatedThermalConductivity)
2740 : {
2741 549850 : auto &s_hbfd = state.dataHeatBalFiniteDiffMgr;
2742 549850 : auto &surfFD = s_hbfd->SurfaceFD(surfNum);
2743 :
2744 549850 : updatedSpecificHeat = mat->getCurrentSpecificHeat(temperaturePrevious,
2745 : temperatureUpdated,
2746 549850 : surfFD.PhaseChangeTemperatureReverse(finiteDifferenceLayerIndex),
2747 549850 : surfFD.PhaseChangeStateOld(finiteDifferenceLayerIndex),
2748 : surfFD.PhaseChangeState(finiteDifferenceLayerIndex));
2749 :
2750 549850 : surfFD.PhaseChangeStateRep(finiteDifferenceLayerIndex) = Material::phaseInts[(int)surfFD.PhaseChangeState(finiteDifferenceLayerIndex)];
2751 549850 : updatedDensity = mat->getDensity(temperaturePrevious);
2752 549850 : updatedThermalConductivity = mat->getConductivity(temperatureUpdated);
2753 549850 : }
2754 :
2755 64 : bool findAnySurfacesUsingConstructionAndCondFD(EnergyPlusData const &state, int const constructionNum)
2756 : {
2757 269 : for (auto const &thisSurface : state.dataSurface->Surface) {
2758 263 : if (thisSurface.Construction == constructionNum) {
2759 68 : if (thisSurface.HeatTransferAlgorithm == DataSurfaces::HeatTransferModel::CondFD) {
2760 58 : return true;
2761 : }
2762 : }
2763 122 : }
2764 6 : return false;
2765 : }
2766 :
2767 : } // namespace HeatBalFiniteDiffManager
2768 :
2769 : } // namespace EnergyPlus
|
