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