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