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