Line data Source code
1 : // EnergyPlus, Copyright (c) 1996-2025, The Board of Trustees of the University of Illinois,
2 : // The Regents of the University of California, through Lawrence Berkeley National Laboratory
3 : // (subject to receipt of any required approvals from the U.S. Dept. of Energy), Oak Ridge
4 : // National Laboratory, managed by UT-Battelle, Alliance for Sustainable Energy, LLC, and other
5 : // contributors. All rights reserved.
6 : //
7 : // NOTICE: This Software was developed under funding from the U.S. Department of Energy and the
8 : // U.S. Government consequently retains certain rights. As such, the U.S. Government has been
9 : // granted for itself and others acting on its behalf a paid-up, nonexclusive, irrevocable,
10 : // worldwide license in the Software to reproduce, distribute copies to the public, prepare
11 : // derivative works, and perform publicly and display publicly, and to permit others to do so.
12 : //
13 : // Redistribution and use in source and binary forms, with or without modification, are permitted
14 : // provided that the following conditions are met:
15 : //
16 : // (1) Redistributions of source code must retain the above copyright notice, this list of
17 : // conditions and the following disclaimer.
18 : //
19 : // (2) Redistributions in binary form must reproduce the above copyright notice, this list of
20 : // conditions and the following disclaimer in the documentation and/or other materials
21 : // provided with the distribution.
22 : //
23 : // (3) Neither the name of the University of California, Lawrence Berkeley National Laboratory,
24 : // the University of Illinois, U.S. Dept. of Energy nor the names of its contributors may be
25 : // used to endorse or promote products derived from this software without specific prior
26 : // written permission.
27 : //
28 : // (4) Use of EnergyPlus(TM) Name. If Licensee (i) distributes the software in stand-alone form
29 : // without changes from the version obtained under this License, or (ii) Licensee makes a
30 : // reference solely to the software portion of its product, Licensee must refer to the
31 : // software as "EnergyPlus version X" software, where "X" is the version number Licensee
32 : // obtained under this License and may not use a different name for the software. Except as
33 : // specifically required in this Section (4), Licensee shall not use in a company name, a
34 : // product name, in advertising, publicity, or other promotional activities any name, trade
35 : // name, trademark, logo, or other designation of "EnergyPlus", "E+", "e+" or confusingly
36 : // similar designation, without the U.S. Department of Energy's prior written consent.
37 : //
38 : // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
39 : // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY
40 : // AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
41 : // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
42 : // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
43 : // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
44 : // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
45 : // OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
46 : // POSSIBILITY OF SUCH DAMAGE.
47 :
48 : // EnergyPlus Headers
49 : #include <EnergyPlus/Construction.hh>
50 : #include <EnergyPlus/Data/EnergyPlusData.hh>
51 : #include <EnergyPlus/DataConversions.hh>
52 : #include <EnergyPlus/DisplayRoutines.hh>
53 : #include <EnergyPlus/Material.hh>
54 : #include <EnergyPlus/OutputReportPredefined.hh>
55 : #include <EnergyPlus/UtilityRoutines.hh>
56 :
57 : namespace EnergyPlus::Construction {
58 :
59 300 : void ConstructionProps::calculateTransferFunction(EnergyPlusData &state, bool &ErrorsFound, bool &DoCTFErrorReport)
60 : {
61 :
62 : // SUBROUTINE INFORMATION:
63 : // AUTHOR Russ Taylor
64 : // DATE WRITTEN June 1990
65 : // MODIFIED July 1994, LKL, cosmetic and improve execution time
66 : // Dec 1995, Apr 1996, RKS, cosmetic and clean-up changes, changes to allow proper
67 : // handling of resistive layers
68 : // June 2000, RKS, addition of QTFs (both 1- and 2-D solutions for constructions
69 : // with embedded/internal heat sources/sinks)
70 : // July 2010-August 2011, RKS, R-value only layer enhancement
71 : // RE-ENGINEERED June 1996, February 1997, August-October 1997, RKS; Nov 1999, LKL
72 :
73 : // PURPOSE OF THIS SUBROUTINE:
74 : // This subroutine serves as the main drive for the
75 : // calculation of Conduction Transfer Functions (CTFs)
76 : // using the state space method.
77 :
78 : // METHODOLOGY EMPLOYED:
79 : // The basic steps of this routine (which may be a little difficult
80 : // to decipher until another major revision is done) are:
81 : // 1. Determine if enough material info has been entered
82 : // 2. Determine whether construct is (a) all resistive,
83 : // (b) the reverse of a previously calculated construct, or
84 : // (c) neither (a) nor (b), i.e. a layer for which CTFs must
85 : // be calculated.
86 : // 3. If the answer to 2 is (a), calculate the overall resistance
87 : // and use this as the CTF (steady state conduction).
88 : // 4. If the answer to 2 is (b), transfer the CTFs for the reverse
89 : // construction to the CTF arrays for this construct (reversing
90 : // the inside and outside terms).
91 : // 5. If the answer to 2 is (c), calculate the CTFs using the state
92 : // space method described below.
93 : // The state space method of calculating CTFs involves
94 : // applying a finite difference grid to a multilayered
95 : // building element and performing linear algebra on the
96 : // resulting system of equations (in matrix form).
97 : // CTFs must be calculated for non-reversed layers which
98 : // have an appreciable thermal mass. A conversion from
99 : // SI units to English units is made due to concerns
100 : // about round off problems noted in earlier version of
101 : // this subroutine.
102 :
103 : // REFERENCES:
104 : // Seem, J.E. "Modeling of Heat Transfer in Buildings",
105 : // Department of Mechanical Engineering, University of
106 : // Wisconsin-Madison, 1987.
107 : // Strand, R.K. "Testing Design Description for the CTF
108 : // Calculation Code in BEST", BSO internal document,
109 : // May/June 1996.
110 : // Strand, R.K. "Heat Source Transfer Functions and Their
111 : // Applicatoin to Low Temperature Radiant Heating System",
112 : // Ph.D. Dissertation, Department of Mechanical and
113 : // Industrial Engineering, University of Illinois at
114 : // Urbana-Champaign, 1995.
115 :
116 : // SUBROUTINE PARAMETER DEFINITIONS:
117 300 : constexpr Real64 PhysPropLimit(1.0e-6); // Physical properties limit.
118 : // This is more or less the traditional value from BLAST.
119 :
120 300 : constexpr Real64 RValueLowLimit(1.0e-3); // Physical properties limit for R-value only layers
121 : // This value was based on trial and error related to CR 7791 where a
122 : // user had entered a "no insulation" layer with an R-value of 1.0E-05.
123 : // Some trial and error established this as a potential value though
124 : // there is no guarantee that this is a good value.
125 :
126 300 : constexpr int MinNodes(6); // Minimum number of state space nodes
127 : // per layer. This value was chosen based on experience with IBLAST.
128 :
129 300 : constexpr Real64 MaxAllowedCTFSumError(0.01); // Allow a 1 percent
130 : // difference between the CTF series summations. If the difference is
131 : // greater than this, then the coefficients will not yield a valid steady
132 : // state solution.
133 :
134 300 : constexpr Real64 MaxAllowedTimeStep = 7.0; // Sets the maximum allowed time step
135 : // for CTF calculations to be 7 hours. This is done in response to some
136 : // rare situations where odd or faulty input will cause the routine to
137 : // go off and get some huge time step (in excess of 20 hours). This value
138 : // is a compromise that does not really solve any input problems. One run
139 : // indicated that 2 meters of concrete will result in a time step of slightly
140 : // more than 3 hours. So, 7 hours was arbitrarily picked as a ceiling for
141 : // time steps so that an error message can be produced to warn the user
142 : // that something isn't right. Note that the 7 hour limit does not guarantee
143 : // that problems won't exist and it does not necessarily avoid any problems
144 : // that interpolated temperature histories might cause.
145 :
146 300 : this->CTFCross.fill(0.0);
147 300 : this->CTFFlux.fill(0.0);
148 300 : this->CTFInside.fill(0.0);
149 300 : this->CTFOutside.fill(0.0);
150 300 : this->CTFSourceIn.fill(0.0);
151 300 : this->CTFSourceOut.fill(0.0);
152 300 : this->CTFTimeStep = 0.0;
153 300 : this->CTFTSourceOut.fill(0.0);
154 300 : this->CTFTSourceIn.fill(0.0);
155 300 : this->CTFTSourceQ.fill(0.0);
156 300 : this->CTFTUserOut.fill(0.0);
157 300 : this->CTFTUserIn.fill(0.0);
158 300 : this->CTFTUserSource.fill(0.0);
159 300 : this->NumHistories = 0;
160 300 : this->NumCTFTerms = 0;
161 300 : this->UValue = 0.0;
162 :
163 300 : if (!this->IsUsedCTF) {
164 54 : return;
165 : }
166 :
167 246 : Array1D<Real64> cp(Construction::MaxLayersInConstruct); // Specific heat of a material layer
168 246 : Array1D<Real64> dl(Construction::MaxLayersInConstruct); // Thickness of a material layer
169 246 : Array1D<Real64> dx(Construction::MaxLayersInConstruct); // Distance between nodes in a particular material layer
170 246 : Array1D<Real64> lr(Construction::MaxLayersInConstruct); // R value of a material layer
171 246 : Array1D_int Nodes(Construction::MaxLayersInConstruct); // Array containing the number of nodes per layer
172 246 : Array1D_bool ResLayer(Construction::MaxLayersInConstruct); // Set true if the layer must be handled as a resistive
173 246 : Array1D<Real64> rho(Construction::MaxLayersInConstruct); // Density of a material layer
174 246 : Array1D<Real64> rk(Construction::MaxLayersInConstruct); // Thermal conductivity of a material layer
175 246 : Array1D_int AdjacentResLayerNum(Construction::MaxLayersInConstruct); // Layers that are adjacent to each other which are resistive
176 :
177 : Real64 amatx; // Intermediate calculation variable
178 : Real64 amatxx; // Intermediate calculation variable
179 : Real64 amaty; // Intermediate calculation variable
180 : Real64 BiggestSum; // Largest CTF series summation (maximum of SumXi, SumYi, and SumZi)
181 : Real64 cap; // Thermal capacitance of a node (intermediate calculation)
182 : Real64 capavg; // Thermal capacitance of a node (average value for a node at an interface)
183 : Real64 cnd; // Total thermal conductance (1/Rtot) of the bldg element
184 : Real64 dtn; // Intermediate calculation of the time step
185 : Real64 dxn; // Intermediate calculation of nodal spacing
186 : Real64 dxtmp; // Intermediate calculation variable ( = 1/dx/cap)
187 : Real64 dyn; // Nodal spacing in the direction perpendicular to the main direction
188 : Real64 SumXi; // Summation of all of the Xi terms (inside CTFs) for a construction
189 : Real64 SumYi; // Summation of all of the Xi terms (cross CTFs) for a construction
190 : Real64 SumZi; // Summation of all of the Xi terms (outside CTFs) for a construction
191 :
192 : Real64 DeltaTimestep; // zone timestep in seconds, for local check of properties
193 :
194 246 : this->CTFTimeStep = state.dataGlobal->TimeStepZone;
195 246 : Real64 rs = 0.0;
196 246 : int LayersInConstruct = 0;
197 246 : int NumResLayers = 0;
198 246 : ResLayer = false;
199 246 : AdjacentResLayerNum = 0; // Zero this out for each construct
200 :
201 647 : for (int Layer = 1; Layer <= this->TotLayers; ++Layer) { // Begin layer loop ...
202 :
203 : // Loop through all of the layers in the current construct. The purpose
204 : // of this loop is to define the thermal properties necessary to
205 : // calculate the CTFs.
206 :
207 401 : int CurrentLayer = this->LayerPoint(Layer);
208 :
209 401 : ++LayersInConstruct;
210 :
211 : // Obtain thermal properties from the Material derived type
212 :
213 401 : auto *thisMaterial = state.dataMaterial->materials(CurrentLayer);
214 401 : assert(thisMaterial != nullptr);
215 :
216 401 : dl(Layer) = thisMaterial->Thickness;
217 401 : rk(Layer) = thisMaterial->Conductivity;
218 401 : rho(Layer) = thisMaterial->Density;
219 401 : cp(Layer) = thisMaterial->SpecHeat; // Must convert
220 : // from kJ/kg-K to J/kg-k due to rk units
221 :
222 401 : if (this->SourceSinkPresent && !thisMaterial->WarnedForHighDiffusivity) {
223 : // check for materials that are too conductive or thin
224 0 : if ((rho(Layer) * cp(Layer)) > 0.0) {
225 0 : Real64 Alpha = rk(Layer) / (rho(Layer) * cp(Layer));
226 0 : if (Alpha > DataHeatBalance::HighDiffusivityThreshold) {
227 0 : DeltaTimestep = state.dataGlobal->TimeStepZoneSec;
228 0 : Real64 const ThicknessThreshold = std::sqrt(Alpha * DeltaTimestep * 3.0);
229 0 : if (thisMaterial->Thickness < ThicknessThreshold) {
230 0 : ShowSevereError(
231 : state,
232 0 : format(
233 : "InitConductionTransferFunctions: Found Material that is too thin and/or too highly conductive, material name = {}",
234 0 : thisMaterial->Name));
235 0 : ShowContinueError(state,
236 0 : format("High conductivity Material layers are not well supported for internal source constructions, "
237 : "material conductivity = {:.3R} [W/m-K]",
238 0 : thisMaterial->Conductivity));
239 0 : ShowContinueError(state, format("Material thermal diffusivity = {:.3R} [m2/s]", Alpha));
240 0 : ShowContinueError(state,
241 0 : format("Material with this thermal diffusivity should have thickness > {:.5R} [m]", ThicknessThreshold));
242 0 : if (thisMaterial->Thickness < DataHeatBalance::ThinMaterialLayerThreshold) {
243 0 : ShowContinueError(state,
244 0 : format("Material may be too thin to be modeled well, thickness = {:.5R} [m]", thisMaterial->Thickness));
245 0 : ShowContinueError(state,
246 0 : format("Material with this thermal diffusivity should have thickness > {:.5R} [m]",
247 : DataHeatBalance::ThinMaterialLayerThreshold));
248 : }
249 0 : thisMaterial->WarnedForHighDiffusivity = true;
250 : }
251 : }
252 : }
253 : }
254 401 : if (thisMaterial->Thickness > 3.0) {
255 0 : ShowSevereError(state, "InitConductionTransferFunctions: Material too thick for CTF calculation");
256 0 : ShowContinueError(state, format("material name = {}", thisMaterial->Name));
257 0 : ErrorsFound = true;
258 : }
259 :
260 401 : if (rk(Layer) <= PhysPropLimit) { // Thermal conductivity too small,
261 : // thus this must be handled as a resistive layer
262 135 : ResLayer(Layer) = true;
263 : } else {
264 266 : lr(Layer) = dl(Layer) / rk(Layer);
265 266 : ResLayer(Layer) = (dl(Layer) * std::sqrt(rho(Layer) * cp(Layer) / rk(Layer))) < PhysPropLimit;
266 : }
267 :
268 : // If not a resistive layer, nothing further is required
269 : // for this layer.
270 :
271 401 : if (ResLayer(Layer)) { // Resistive layer-check for R-value, etc.
272 135 : ++NumResLayers; // Increment number of resistive layers
273 135 : lr(Layer) = thisMaterial->Resistance; // User defined thermal resistivity
274 135 : if (lr(Layer) < RValueLowLimit) { // User didn't define enough
275 : // parameters to calculate CTFs for a building element
276 : // containing this layer.
277 :
278 0 : ShowSevereError(state, format("InitConductionTransferFunctions: Material={}R Value below lowest allowed value", thisMaterial->Name));
279 0 : ShowContinueError(state, format("Lowest allowed value=[{:.3R}], Material R Value=[{:.3R}].", RValueLowLimit, lr(Layer)));
280 0 : ErrorsFound = true;
281 :
282 : } else { // A valid user defined R-value is available.
283 : // If this is either the first or last layer in the construction,
284 : // then assign other properties based on air at 1 atm, 300K.
285 : // Reference for air properties: Incropera and DeWitt,
286 : // Introduction to Heat Transfer, Appendix A, Table A.4,
287 : // John Wiley & Sons, New York, 1985.
288 : // If this is not the first or last layer in the construction,
289 : // then use the "exact" approach to model a massless layer
290 : // based on the node equations for the state space method.
291 :
292 135 : if ((Layer == 1) || (Layer == this->TotLayers) || (!state.dataMaterial->materials(this->LayerPoint(Layer))->ROnly)) {
293 117 : cp(Layer) = 1.007;
294 117 : rho(Layer) = 1.1614;
295 117 : rk(Layer) = 0.0263;
296 117 : dl(Layer) = rk(Layer) * lr(Layer);
297 : } else {
298 18 : cp(Layer) = 0.0;
299 18 : rho(Layer) = 0.0;
300 18 : rk(Layer) = 1.0;
301 18 : dl(Layer) = lr(Layer);
302 : }
303 : }
304 : } // ... end of resistive layer determination IF-THEN block.
305 : } // ... end of layer loop.
306 :
307 : // If errors have been found, just return
308 :
309 246 : if (ErrorsFound) return;
310 :
311 : // Combine any adjacent resistive-only (no mass) layers together
312 : // to avoid a divide by zero error in the CTF calculations below.
313 : // Since the inner and outer layers cannot be resistive layers
314 : // (inner and outer layer still converted to equivalent air layer)
315 : // there can only be resistive layers adjacent to one another if
316 : // there are more than three total layers and more than one
317 : // resistive layer.
318 246 : if ((LayersInConstruct > 3) && (NumResLayers > 1)) {
319 1 : int NumAdjResLayers = 0;
320 3 : for (int Layer = 2; Layer <= LayersInConstruct - 2; ++Layer) {
321 2 : if ((ResLayer(Layer)) && (ResLayer(Layer + 1))) {
322 0 : ++NumAdjResLayers;
323 : // There is method to the next assignment statement. As the layers get shifted, the layer
324 : // numbers will also shift. Thus, we have to also shift which layer we are dealing with.
325 0 : AdjacentResLayerNum(NumAdjResLayers) = Layer + 1 - NumAdjResLayers;
326 : }
327 : }
328 1 : for (int AdjLayer = 1; AdjLayer <= NumAdjResLayers; ++AdjLayer) {
329 0 : int Layer = AdjacentResLayerNum(AdjLayer);
330 : // Double check to make sure we are in the right place...
331 0 : if ((ResLayer(Layer)) && (ResLayer(Layer + 1))) {
332 : // Shift layers forward after combining two adjacent layers. Then
333 : // restart the do loop.
334 0 : cp(Layer) = 0.0;
335 0 : rho(Layer) = 0.0;
336 0 : rk(Layer) = 1.0;
337 0 : lr(Layer) += lr(Layer + 1);
338 0 : dl(Layer) = lr(Layer);
339 0 : --NumResLayers; // Combining layers so decrease number of resistive layers
340 0 : for (int Layer1 = Layer + 1; Layer1 <= LayersInConstruct - 1; ++Layer1) {
341 0 : lr(Layer1) = lr(Layer1 + 1);
342 0 : dl(Layer1) = dl(Layer1 + 1);
343 0 : rk(Layer1) = rk(Layer1 + 1);
344 0 : rho(Layer1) = rho(Layer1 + 1);
345 0 : cp(Layer1) = cp(Layer1 + 1);
346 0 : ResLayer(Layer1) = ResLayer(Layer1 + 1);
347 : }
348 : // Then zero out the layer that got shifted forward
349 0 : cp(LayersInConstruct) = 0.0;
350 0 : rho(LayersInConstruct) = 0.0;
351 0 : rk(LayersInConstruct) = 0.0;
352 0 : lr(LayersInConstruct) = 0.0;
353 0 : dl(LayersInConstruct) = 0.0;
354 : // Now reduce the number of layers in construct since merger is complete
355 0 : --LayersInConstruct;
356 : // Also adjust layers with source/sinks if two layers are merged
357 0 : if (this->SourceSinkPresent) {
358 0 : --this->SourceAfterLayer;
359 0 : --this->TempAfterLayer;
360 : }
361 : } else { // These are not adjacent layers and there is a logic flaw here (should not happen)
362 0 : ShowFatalError(state, format("Combining resistance layers failed for {}", this->Name));
363 0 : ShowContinueError(state, "This should never happen. Contact EnergyPlus Support for further assistance.");
364 : }
365 : }
366 : }
367 :
368 : // Convert SI units to English. In theory, conversion to English
369 : // units is not necessary; however, Russ Taylor noted that some
370 : // numerical problems when SI units were used and decided to continue
371 : // calculating CTFs in English units.
372 :
373 647 : for (int Layer = 1; Layer <= LayersInConstruct; ++Layer) { // Begin units conversion loop ...
374 :
375 401 : lr(Layer) *= DataConversions::CFU;
376 401 : dl(Layer) /= DataConversions::CFL;
377 401 : rk(Layer) /= DataConversions::CFK;
378 401 : rho(Layer) /= DataConversions::CFD;
379 401 : cp(Layer) /= (DataConversions::CFC * 1000.0);
380 :
381 : } // ... end of layer loop for units conversion.
382 :
383 246 : if (this->SolutionDimensions == 1) {
384 246 : dyn = 0.0;
385 : } else {
386 0 : dyn = (this->ThicknessPerpend / DataConversions::CFL) / double(NumOfPerpendNodes - 1);
387 : }
388 :
389 : // Compute total construct conductivity and resistivity.
390 :
391 647 : for (int Layer = 1; Layer <= LayersInConstruct; ++Layer) {
392 401 : rs += lr(Layer); // Resistances in series sum algebraically
393 : }
394 :
395 246 : cnd = 1.0 / rs; // Conductivity is the inverse of resistivity
396 :
397 246 : bool RevConst = false;
398 :
399 246 : if (LayersInConstruct > NumResLayers) {
400 :
401 : // One or more are not simple resistive layers so CTFs will have to be
402 : // calculated unless this is a reverse of a previously defined
403 : // this->
404 :
405 : // Check for reversed construction of interzone surfaces by checking
406 : // previous constructions for same number of layers as first indicator.
407 :
408 : // previously this loop would go from 1..currentConstructionIndex-1
409 : // instead of that, we'll loop through the list and stop when we get to the current construction
410 : // should be the same behavior, we're just checking it by address
411 352 : for (auto &otherConstruction : state.dataConstruction->Construct) {
412 352 : if (&otherConstruction == this) break;
413 :
414 : // If a source or sink is present in this construction, do not allow any
415 : // checks for reversed constructions, i.e., always force EnergyPlus to
416 : // calculate CTF/QTFs. So, don't even check for reversed constructions.
417 206 : if (this->SourceSinkPresent) break; // Constr DO loop
418 :
419 206 : if (this->TotLayers == otherConstruction.TotLayers) { // Same number of layers--now | check for reversed construct.
420 :
421 80 : RevConst = true;
422 :
423 90 : for (int Layer = 1; Layer <= this->TotLayers; ++Layer) { // Begin layers loop ...
424 :
425 : // RevConst is set to FALSE anytime a mismatch in materials is found.
426 : // This will exit this DO immediately and go on to the next construct
427 : // (if any remain).
428 :
429 87 : int OppositeLayer = this->TotLayers - Layer + 1;
430 :
431 87 : if (this->LayerPoint(Layer) != otherConstruction.LayerPoint(OppositeLayer)) {
432 :
433 77 : RevConst = false;
434 77 : break; // Layer DO loop
435 : }
436 :
437 : } // ... end of layers loop.
438 :
439 : // If the reverse construction isn't used by any surfaces then the CTFs
440 : // still need to be defined.
441 80 : if (RevConst && !otherConstruction.IsUsedCTF) {
442 0 : RevConst = false;
443 : }
444 :
445 80 : if (RevConst) { // Curent construction is a reverse of
446 : // construction Constr. Thus, CTFs do not need to be re-
447 : // calculated. Copy CTF info for construction Constr to
448 : // construction ConstrNum.
449 :
450 3 : this->CTFTimeStep = otherConstruction.CTFTimeStep;
451 3 : this->NumHistories = otherConstruction.NumHistories;
452 3 : this->NumCTFTerms = otherConstruction.NumCTFTerms;
453 :
454 : // Transfer the temperature and flux history terms to CTF arrays.
455 : // Loop through the number of CTF history terms ...
456 29 : for (int HistTerm = 0; HistTerm <= this->NumCTFTerms; ++HistTerm) {
457 :
458 26 : this->CTFInside[HistTerm] = otherConstruction.CTFOutside[HistTerm];
459 26 : this->CTFCross[HistTerm] = otherConstruction.CTFCross[HistTerm];
460 26 : this->CTFOutside[HistTerm] = otherConstruction.CTFInside[HistTerm];
461 26 : if (HistTerm != 0) this->CTFFlux[HistTerm] = otherConstruction.CTFFlux[HistTerm];
462 :
463 : } // ... end of CTF history terms loop.
464 :
465 3 : break; // Constr DO loop
466 :
467 : } // ... end of reversed construction found block
468 :
469 : } // ... end of reversed construct (same number of layers) block.
470 :
471 : } // ... end of construct loop (check reversed--Constr)
472 :
473 149 : if (!RevConst) { // Calculate CTFs (non-reversed constr)
474 :
475 : // Estimate number of nodes each layer of the construct will require
476 : // and calculate the nodal spacing from that
477 :
478 444 : for (int Layer = 1; Layer <= LayersInConstruct; ++Layer) { // Begin loop thru layers ...
479 :
480 : // The calculation of dxn used here is based on a standard stability
481 : // criteria for explicit finite difference solutions. This criteria
482 : // was chosen not because it is viewed to be correct, but rather for
483 : // lack of any better criteria at this time. The use of a Fourier
484 : // number based criteria such as this is probably physically correct,
485 : // though the coefficient (2.0) may not be.
486 :
487 : // If this is a "resistive" layer, only need a single node
488 298 : if ((ResLayer(Layer)) && (Layer > 1) && (Layer < LayersInConstruct)) {
489 18 : Nodes(Layer) = 1;
490 18 : dx(Layer) = dl(Layer);
491 : } else {
492 280 : dxn = std::sqrt(2.0 * (rk(Layer) / rho(Layer) / cp(Layer)) * this->CTFTimeStep);
493 :
494 280 : int ipts1 = int(dl(Layer) / dxn); // number of nodes=thickness/spacing
495 :
496 : // Limit the upper and lower bounds of the number of
497 : // nodes to MaxCTFTerms and MinNodes respectively.
498 :
499 280 : if (ipts1 > Construction::MaxCTFTerms) { // Too many nodes
500 0 : Nodes(Layer) = Construction::MaxCTFTerms;
501 280 : } else if (ipts1 < MinNodes) { // Too few nodes
502 264 : Nodes(Layer) = MinNodes;
503 : } else { // Calculated number of nodes ok
504 16 : Nodes(Layer) = ipts1;
505 : }
506 :
507 280 : dx(Layer) = dl(Layer) / double(Nodes(Layer)); // calc node spacing
508 : }
509 :
510 : } // . .. end of layers in construction loop (calculating #nodes per layer)
511 :
512 : // Determine the total number of nodes (rcmax)
513 :
514 146 : this->rcmax = 0;
515 444 : for (int Layer = 1; Layer <= LayersInConstruct; ++Layer) {
516 298 : this->rcmax += Nodes(Layer);
517 : }
518 :
519 : // Nodes are placed throughout layers and at the interface between
520 : // layers. As a result, the end layers share a node with the adjacent
521 : // layer-leaving one less node total for all layers.
522 :
523 146 : --this->rcmax;
524 146 : if (this->SolutionDimensions > 1) this->rcmax *= NumOfPerpendNodes;
525 :
526 : // This section no longer needed as rcmax/number of total nodes is allowed to float.
527 : // If reinstated, this node reduction section would have to be modified to account for
528 : // the possibility that a 2-D solution is potentially being performed.
529 : // Check to see if the maximum number of nodes for the construct has
530 : // been exceeded. Reduce the nodes per layer if necessary, but only
531 : // if the number of nodes in a particular layer is greater than the
532 : // minimum node limit.
533 :
534 : // DO WHILE (rcmax > MaxTotNodes) ! Begin total node reduction loop ...
535 :
536 : // rcmax = 0
537 :
538 : // DO Layer = 1, LayersInConstruct ! Begin layer node reduction ...
539 :
540 : // ! If more nodes than the minimum limit for a layer, reduce the
541 : // ! number of nodes.
542 :
543 : // IF (Nodes(Layer) > MinNodes) THEN
544 : // Nodes(Layer) = Nodes(Layer)-1
545 : // dx(Layer) = dl(Layer)/DBLE(Nodes(Layer)) ! Recalc node spacing
546 : // END IF
547 :
548 : // rcmax = rcmax + Nodes(Layer) ! Recalculate total number of nodes
549 :
550 : // END DO ! ... end of layer loop for node reduction.
551 :
552 : // rcmax = rcmax-1 ! See note above on counting rcmax
553 :
554 : // END DO ! ... end of total node reduction loop.
555 :
556 : // For constructions that have sources or sinks present, determine which
557 : // node the source/sink is applied at and also where the temperature
558 : // calculation has been requested.
559 146 : this->setNodeSourceAndUserTemp(Nodes);
560 :
561 : // "Adjust time step to ensure stability." If the time step is too
562 : // small, it will result in too many history terms which can lead to
563 : // solution instability. The method used here to determine whether or
564 : // not the time step will produce a stable solution is based on a pure
565 : // Fourier number calculation (Fo = 1) and has not proven to be
566 : // completely effective. If too many history terms are calculated,
567 : // the time step is adjusted and the CTFs end up being recalculated
568 : // (see later code in this routine).
569 :
570 146 : dtn = 0.0;
571 146 : this->CTFTimeStep = 0.0;
572 444 : for (int Layer = 1; Layer <= LayersInConstruct; ++Layer) {
573 298 : if (Nodes(Layer) >= Construction::MaxCTFTerms) {
574 0 : if (this->SolutionDimensions == 1) {
575 0 : dtn = rho(Layer) * cp(Layer) * pow_2(dx(Layer)) / rk(Layer);
576 : } else { // 2-D solution requested-->this changes length parameter in Fourier number calculation
577 0 : dtn = rho(Layer) * cp(Layer) * (pow_2(dx(Layer)) + pow_2(dyn)) / rk(Layer);
578 : }
579 0 : if (dtn > this->CTFTimeStep) this->CTFTimeStep = dtn;
580 : }
581 : }
582 :
583 : // If the user defined time step is significantly different than the
584 : // calculated time step for this construct, then CTFTimeStep must be
585 : // revised.
586 :
587 146 : if (std::abs((state.dataGlobal->TimeStepZone - this->CTFTimeStep) / state.dataGlobal->TimeStepZone) > 0.1) {
588 :
589 146 : if (this->CTFTimeStep > state.dataGlobal->TimeStepZone) {
590 :
591 : // CTFTimeStep larger than TimeStepZone: Make sure TimeStepZone
592 : // divides evenly into CTFTimeStep
593 0 : this->NumHistories = int((this->CTFTimeStep / state.dataGlobal->TimeStepZone) + 0.5);
594 0 : this->CTFTimeStep = state.dataGlobal->TimeStepZone * double(this->NumHistories);
595 :
596 : } else {
597 :
598 : // CTFTimeStep smaller than TimeStepZone: Set to TimeStepZone
599 146 : this->CTFTimeStep = state.dataGlobal->TimeStepZone;
600 146 : this->NumHistories = 1;
601 : }
602 : }
603 :
604 : // Calculate the CTFs using the state space method
605 : // outlined in Seem's dissertation. The main matrices
606 : // AMat, BMat, CMat, and DMat must be derived from
607 : // applying a finite difference network to the layers of
608 : // each bldg element.
609 :
610 : // This section must continue looping until the CTFs
611 : // calculated here will produce a stable solution (less
612 : // history terms than MaxCTFTerms).
613 :
614 : // This first subsection calculates the elements of AMat
615 : // which characterizes the heat transfer inside the
616 : // building element.
617 :
618 : // Initialize loop control logical
619 : // Set after CTFs are calculated, based on whether there are too many CTF terms
620 146 : bool CTFConvrg = false;
621 :
622 146 : this->AExp.allocate(this->rcmax, this->rcmax);
623 146 : this->AExp = 0.0;
624 146 : this->AMat.allocate(this->rcmax, this->rcmax);
625 146 : this->AMat = 0.0;
626 146 : this->AInv.allocate(this->rcmax, this->rcmax);
627 146 : this->AInv = 0.0;
628 146 : this->IdenMatrix.allocate(this->rcmax, this->rcmax);
629 146 : this->IdenMatrix = 0.0;
630 1731 : for (int ir = 1; ir <= this->rcmax; ++ir) {
631 1585 : this->IdenMatrix(ir, ir) = 1.0;
632 : }
633 146 : this->e.dimension(this->rcmax, 0.0);
634 146 : this->Gamma1.allocate(3, this->rcmax);
635 146 : this->Gamma1 = 0.0;
636 146 : this->Gamma2.allocate(3, this->rcmax);
637 146 : this->Gamma2 = 0.0;
638 146 : this->s.allocate(3, 4, this->rcmax);
639 146 : this->s = 0.0;
640 :
641 292 : while (!CTFConvrg) { // Begin CTF calculation loop ...
642 :
643 146 : this->BMat(3) = 0.0;
644 :
645 146 : if (this->SolutionDimensions == 1) {
646 :
647 : // Set up intermediate calculations for the first layer.
648 146 : cap = rho(1) * cp(1) * dx(1);
649 146 : cap *= 1.5; // For the first node, account for the fact that the
650 : // half-node at the surface results in a "loss" of some
651 : // thermal mass. Therefore, for simplicity, include it
652 : // at this node. Same thing done at the last node...
653 146 : dxtmp = 1.0 / dx(1) / cap;
654 :
655 146 : this->AMat(1, 1) = -2.0 * rk(1) * dxtmp; // Assign the matrix values for the
656 146 : this->AMat(2, 1) = rk(1) * dxtmp; // first node.
657 146 : this->BMat(1) = rk(1) * dxtmp; // Assign non-zero value of BMat.
658 :
659 146 : int Layer = 1; // Initialize the "layer" counter
660 :
661 146 : int NodeInLayer = 2; // Initialize the node (in a layer) counter (already
662 : // on the second node for the first layer
663 :
664 1439 : for (int Node = 2; Node <= this->rcmax - 1; ++Node) { // Begin nodes loop (includes all nodes except the
665 : // first/last which have special equations) ...
666 :
667 1293 : if ((NodeInLayer == Nodes(Layer)) && (LayersInConstruct != 1)) { // For a node at
668 : // the interface between two adjacent layers, the
669 : // capacitance of the node must be calculated from the 2
670 : // halves which may be made up of 2 different materials.
671 :
672 152 : cap = (rho(Layer) * cp(Layer) * dx(Layer) + rho(Layer + 1) * cp(Layer + 1) * dx(Layer + 1)) * 0.5;
673 :
674 152 : this->AMat(Node - 1, Node) = rk(Layer) / dx(Layer) / cap; // Assign matrix
675 152 : this->AMat(Node, Node) = -1.0 * (rk(Layer) / dx(Layer) + rk(Layer + 1) / dx(Layer + 1)) / cap; // values for | the current
676 152 : this->AMat(Node + 1, Node) = rk(Layer + 1) / dx(Layer + 1) / cap; // node.
677 :
678 152 : NodeInLayer = 0; // At an interface, reset nodes in layer counter
679 152 : ++Layer; // Also increment the layer counter
680 :
681 : } else { // Standard node within any layer
682 :
683 1141 : cap = rho(Layer) * cp(Layer) * dx(Layer); // Intermediate
684 1141 : dxtmp = 1.0 / dx(Layer) / cap; // calculations.
685 1141 : this->AMat(Node - 1, Node) = rk(Layer) * dxtmp; // Assign matrix
686 1141 : this->AMat(Node, Node) = -2.0 * rk(Layer) * dxtmp; // values for the
687 1141 : this->AMat(Node + 1, Node) = rk(Layer) * dxtmp; // current node.
688 : }
689 :
690 1293 : ++NodeInLayer; // Increment nodes in layer counter
691 1293 : if (Node == this->NodeSource) this->BMat(3) = 1.0 / cap;
692 :
693 : } // ... end of nodes loop.
694 :
695 : // Intermediate calculations for the last node.
696 146 : cap = rho(LayersInConstruct) * cp(LayersInConstruct) * dx(LayersInConstruct);
697 146 : cap *= 1.5; // For the last node, account for the fact that the
698 : // half-node at the surface results in a "loss" of some
699 : // thermal mass. Therefore, for simplicity, include it
700 : // at this node. Same thing done at the first node...
701 146 : dxtmp = 1.0 / dx(LayersInConstruct) / cap;
702 :
703 146 : this->AMat(this->rcmax, this->rcmax) = -2.0 * rk(LayersInConstruct) * dxtmp; // Assign matrix
704 146 : this->AMat(this->rcmax - 1, this->rcmax) = rk(LayersInConstruct) * dxtmp; // values for the
705 146 : this->BMat(2) = rk(LayersInConstruct) * dxtmp; // last node.
706 :
707 146 : this->CMat(1) = -rk(1) / dx(1); // Compute the necessary elements
708 146 : this->CMat(2) = rk(LayersInConstruct) / dx(LayersInConstruct); // of all other
709 146 : this->DMat(1) = rk(1) / dx(1); // matrices for the state
710 146 : this->DMat(2) = -rk(LayersInConstruct) / dx(LayersInConstruct); // space method
711 :
712 : } else { // 2-D solution requested (assign matrices appropriately)
713 :
714 : // As with the 1-D solution, we are accounting for the thermal mass
715 : // of the half-node at the surface by adding it to the first row
716 : // of interior nodes at both sides of the construction. This is not
717 : // exact, but it does take all of the thermal mass into account.
718 0 : amatx = rk(1) / (1.5 * rho(1) * cp(1) * dx(1) * dx(1));
719 0 : amaty = rk(1) / (1.5 * rho(1) * cp(1) * dyn * dyn);
720 :
721 : // FIRST ROW OF NODES: This first row within the first material layer
722 : // is special in that it is exposed to a boundary condition. Thus,
723 : // the equations are slightly different.
724 : // Note also that the first and last nodes in a row are slightly
725 : // different from the rest since they are on an adiabatic plane in
726 : // the direction perpendicular to the main direction of heat transfer.
727 0 : this->AMat(1, 1) = -2.0 * (amatx + amaty);
728 0 : this->AMat(2, 1) = 2.0 * amaty;
729 0 : this->AMat(this->NumOfPerpendNodes + 1, 1) = amatx;
730 :
731 0 : for (int Node = 2; Node <= this->NumOfPerpendNodes - 1; ++Node) {
732 0 : this->AMat(Node - 1, Node) = amaty;
733 0 : this->AMat(Node, Node) = -2.0 * (amatx + amaty);
734 0 : this->AMat(Node + 1, Node) = amaty;
735 0 : this->AMat(Node + this->NumOfPerpendNodes, Node) = amatx;
736 : }
737 :
738 0 : this->AMat(this->NumOfPerpendNodes, this->NumOfPerpendNodes) = -2.0 * (amatx + amaty);
739 0 : this->AMat(this->NumOfPerpendNodes - 1, this->NumOfPerpendNodes) = 2.0 * amaty;
740 0 : this->AMat(this->NumOfPerpendNodes + this->NumOfPerpendNodes, this->NumOfPerpendNodes) = amatx;
741 :
742 0 : BMat(1) = amatx;
743 :
744 0 : int Layer = 1;
745 0 : int NodeInLayer = 2;
746 0 : amatx = rk(1) / (rho(1) * cp(1) * dx(1) * dx(1)); // Reset these to the normal capacitance
747 0 : amaty = rk(1) / (rho(1) * cp(1) * dyn * dyn); // Reset these to the normal capacitance
748 0 : assert(this->NumOfPerpendNodes > 0); // Autodesk:F2C++ Loop setup assumption
749 0 : int const Node_stop(this->rcmax + 1 - 2 * this->NumOfPerpendNodes);
750 0 : for (int Node = this->NumOfPerpendNodes + 1; Node <= Node_stop; Node += this->NumOfPerpendNodes) {
751 : // INTERNAL ROWS OF NODES: This is the majority of nodes which are all within
752 : // a solid layer and not exposed to a boundary condition.
753 0 : if ((LayersInConstruct == 1) || (NodeInLayer != Nodes(Layer))) {
754 : // Single material row: This row of nodes are all contained within a material
755 : // and thus there is no special considerations necessary.
756 0 : if (NodeInLayer == 1) {
757 : // These intermediate variables only need to be reassigned when a new layer is started.
758 : // When this is simply another row of the same material, these have already been assigned correctly.
759 0 : amatx = rk(Layer) / (rho(Layer) * cp(Layer) * dx(Layer) * dx(Layer));
760 0 : amaty = rk(Layer) / (rho(Layer) * cp(Layer) * dyn * dyn);
761 : }
762 :
763 : // Note that the first and last layers in a row are slightly different
764 : // from the rest since they are on an adiabatic plane in the direction
765 : // perpendicular to the main direction of heat transfer.
766 0 : this->AMat(Node, Node) = -2.0 * (amatx + amaty);
767 0 : this->AMat(Node + 1, Node) = 2.0 * amaty;
768 0 : this->AMat(Node - this->NumOfPerpendNodes, Node) = amatx;
769 0 : this->AMat(Node + this->NumOfPerpendNodes, Node) = amatx;
770 :
771 0 : for (int NodeInRow = 2; NodeInRow <= this->NumOfPerpendNodes - 1; ++NodeInRow) {
772 0 : int Node2 = Node + NodeInRow - 1;
773 0 : this->AMat(Node2 - 1, Node2) = amaty;
774 0 : this->AMat(Node2, Node2) = -2.0 * (amatx + amaty);
775 0 : this->AMat(Node2 + 1, Node2) = amaty;
776 0 : this->AMat(Node2 - this->NumOfPerpendNodes, Node2) = amatx;
777 0 : this->AMat(Node2 + this->NumOfPerpendNodes, Node2) = amatx;
778 : }
779 :
780 0 : int Node2 = Node - 1 + this->NumOfPerpendNodes;
781 0 : this->AMat(Node2, Node2) = -2.0 * (amatx + amaty);
782 0 : this->AMat(Node2 - 1, Node2) = 2.0 * amaty;
783 0 : this->AMat(Node2 - this->NumOfPerpendNodes, Node2) = amatx;
784 0 : this->AMat(Node2 + this->NumOfPerpendNodes, Node2) = amatx;
785 :
786 : } else { // Row at a two-layer interface (half of node consists of one layer's materials
787 : // and the other half consist of the next layer's materials)
788 0 : capavg = 0.5 * (rho(Layer) * cp(Layer) * dx(Layer) + rho(Layer + 1) * cp(Layer + 1) * dx(Layer + 1));
789 0 : amatx = rk(Layer) / (capavg * dx(Layer));
790 0 : amatxx = rk(Layer + 1) / (capavg * dx(Layer + 1));
791 0 : amaty = (rk(Layer) * dx(Layer) + rk(Layer + 1) * dx(Layer + 1)) / (capavg * dyn * dyn);
792 :
793 0 : this->AMat(Node, Node) = -amatx - amatxx - 2.0 * amaty;
794 0 : this->AMat(Node + 1, Node) = 2.0 * amaty;
795 0 : this->AMat(Node - this->NumOfPerpendNodes, Node) = amatx;
796 0 : this->AMat(Node + this->NumOfPerpendNodes, Node) = amatxx;
797 :
798 0 : for (int NodeInRow = 2; NodeInRow <= this->NumOfPerpendNodes - 1; ++NodeInRow) {
799 0 : int Node2 = Node + NodeInRow - 1;
800 0 : this->AMat(Node2 - 1, Node2) = amaty;
801 0 : this->AMat(Node2, Node2) = -amatx - amatxx - 2.0 * amaty;
802 0 : this->AMat(Node2 + 1, Node2) = amaty;
803 0 : this->AMat(Node2 - this->NumOfPerpendNodes, Node2) = amatx;
804 0 : this->AMat(Node2 + this->NumOfPerpendNodes, Node2) = amatxx;
805 : }
806 :
807 0 : int Node2 = Node - 1 + this->NumOfPerpendNodes;
808 0 : this->AMat(Node2, Node2) = -amatx - amatxx - 2.0 * amaty;
809 0 : this->AMat(Node2 - 1, Node2) = 2.0 * amaty;
810 0 : this->AMat(Node2 - this->NumOfPerpendNodes, Node2) = amatx;
811 0 : this->AMat(Node2 + this->NumOfPerpendNodes, Node2) = amatxx;
812 :
813 0 : if (Node == this->NodeSource) BMat(3) = 2.0 * double(this->NumOfPerpendNodes - 1) / capavg;
814 0 : NodeInLayer = 0;
815 0 : ++Layer;
816 : }
817 0 : ++NodeInLayer;
818 : }
819 :
820 : // LAST ROW OF NODES: Like the first row of nodes, this row is exposed to a boundary
821 : // condition and thus has slightly modified nodal equations.
822 :
823 : // As with the 1-D solution, we are accounting for the thermal mass
824 : // of the half-node at the surface by adding it to the first row
825 : // of interior nodes at both sides of the construction. This is not
826 : // exact, but it does take all of the thermal mass into account.
827 0 : amatx /= 1.5;
828 0 : amaty /= 1.5;
829 :
830 0 : int Node = this->rcmax + 1 - this->NumOfPerpendNodes;
831 0 : this->AMat(Node, Node) = -2.0 * (amatx + amaty);
832 0 : this->AMat(Node + 1, Node) = 2.0 * amaty;
833 0 : this->AMat(Node - this->NumOfPerpendNodes, Node) = amatx;
834 :
835 0 : for (int thisNode = this->rcmax + 2 - this->NumOfPerpendNodes; thisNode <= this->rcmax - 1; ++thisNode) {
836 0 : this->AMat(thisNode - 1, thisNode) = amaty;
837 0 : this->AMat(thisNode, thisNode) = -2.0 * (amatx + amaty);
838 0 : this->AMat(thisNode + 1, thisNode) = amaty;
839 0 : this->AMat(thisNode - this->NumOfPerpendNodes, thisNode) = amatx;
840 : }
841 :
842 0 : this->AMat(this->rcmax, this->rcmax) = -2.0 * (amatx + amaty);
843 0 : this->AMat(this->rcmax - 1, this->rcmax) = 2.0 * amaty;
844 0 : this->AMat(this->rcmax - this->NumOfPerpendNodes, this->rcmax) = amatx;
845 :
846 0 : this->BMat(2) = amatx;
847 :
848 0 : this->CMat(1) = -rk(1) / dx(1) / double(this->NumOfPerpendNodes - 1);
849 0 : this->CMat(2) = rk(LayersInConstruct) / dx(LayersInConstruct) / double(this->NumOfPerpendNodes - 1);
850 :
851 0 : this->DMat(1) = rk(1) / dx(1) / double(this->NumOfPerpendNodes - 1);
852 0 : this->DMat(2) = -rk(LayersInConstruct) / dx(LayersInConstruct) / double(this->NumOfPerpendNodes - 1);
853 : }
854 :
855 : // Calculation of the CTFs based on the state space
856 : // method. This process involves finding the exponential
857 : // and inverse of AMat and using these results to
858 : // determine the CTFs. The Gammas are an intermediate
859 : // calculations which are necessary before the CTFs can
860 : // be computed in TransFuncCoeffs.
861 146 : DisplayString(state, "Calculating CTFs for \"" + this->Name + "\"");
862 :
863 : // CALL DisplayNumberAndString(ConstrNum,'Matrix exponential for Construction #')
864 146 : this->calculateExponentialMatrix(); // Compute exponential of AMat
865 :
866 : // CALL DisplayNumberAndString(ConstrNum,'Invert Matrix for Construction #')
867 146 : this->calculateInverseMatrix(); // Compute inverse of AMat
868 :
869 : // CALL DisplayNumberAndString(ConstrNum,'Gamma calculation for Construction #')
870 146 : this->calculateGammas();
871 : // Compute "gamma"s from AMat, AExp, and AInv
872 :
873 : // CALL DisplayNumberAndString(ConstrNum,'Compute CTFs for Construction #')
874 146 : this->calculateFinalCoefficients(); // Compute CTFs
875 :
876 : // Now check to see if the number of transfer functions
877 : // is greater than MaxCTFTerms. If it is, then increase the
878 : // time step and the number of history terms and
879 : // recalculate. Whether or not it will be necessary to
880 : // recalculate the CTFs is controlled by this DO WHILE
881 : // loop and the logical CTFConvrg.
882 :
883 146 : CTFConvrg = true; // Assume solution convergence
884 :
885 : // If too many terms, then solution did not converge. Increase the
886 : // number of histories and the time step. Reset CTFConvrg to continue
887 : // the DO loop.
888 146 : if (this->NumCTFTerms > (Construction::MaxCTFTerms - 1)) {
889 0 : ++this->NumHistories;
890 0 : this->CTFTimeStep += state.dataGlobal->TimeStepZone;
891 0 : CTFConvrg = false;
892 : }
893 :
894 : // If the number of terms is okay, then do a further check on the summation of
895 : // the various series summations. In theory, Sum(Xi) = Sum(Yi) = Sum(Zi). If
896 : // this is not the case, then the terms have not reached a valid solution, and
897 : // we need to increase the number of histories and the time step as above.
898 146 : if (CTFConvrg) {
899 146 : SumXi = this->s0(2, 2);
900 146 : SumYi = this->s0(1, 2);
901 146 : SumZi = this->s0(1, 1);
902 1057 : for (int HistTerm = 1; HistTerm <= this->NumCTFTerms; ++HistTerm) {
903 911 : SumXi += this->s(2, 2, HistTerm);
904 911 : SumYi += this->s(1, 2, HistTerm);
905 911 : SumZi += this->s(1, 1, HistTerm);
906 : }
907 146 : SumXi = std::abs(SumXi);
908 146 : SumYi = std::abs(SumYi);
909 146 : SumZi = std::abs(SumZi);
910 146 : BiggestSum = max(SumXi, SumYi, SumZi);
911 146 : if (BiggestSum > 0.0) {
912 292 : if (((std::abs(SumXi - SumYi) / BiggestSum) > MaxAllowedCTFSumError) ||
913 146 : ((std::abs(SumZi - SumYi) / BiggestSum) > MaxAllowedCTFSumError)) {
914 0 : ++this->NumHistories;
915 0 : this->CTFTimeStep += state.dataGlobal->TimeStepZone;
916 0 : CTFConvrg = false;
917 : }
918 : } else { // Something terribly wrong--the surface has no CTFs, not even an R-value
919 0 : ShowFatalError(state, format("Illegal construction definition, no CTFs calculated for {}", this->Name));
920 : }
921 : }
922 :
923 : // Once the time step has reached a certain point, it is highly likely that
924 : // there is either a problem with the input or the solution. This should
925 : // be extremely rare since other checks should flag most bad user input.
926 : // Thus, if the time step reaches a certain point, error out and let the
927 : // user know that something needs to be checked in the input file.
928 146 : if (this->CTFTimeStep >= MaxAllowedTimeStep) {
929 0 : ShowSevereError(state, format("CTF calculation convergence problem for Construction=\"{}\".", this->Name));
930 0 : ShowContinueError(state, "...with Materials (outside layer to inside)");
931 0 : ShowContinueError(state, format("(outside)=\"{}\"", state.dataMaterial->materials(this->LayerPoint(1))->Name));
932 0 : for (int Layer = 2; Layer <= this->TotLayers; ++Layer) {
933 0 : if (Layer != this->TotLayers) {
934 0 : ShowContinueError(state, format("(next)=\"{}\"", state.dataMaterial->materials(this->LayerPoint(Layer))->Name));
935 : } else {
936 0 : ShowContinueError(state, format("(inside)=\"{}\"", state.dataMaterial->materials(this->LayerPoint(Layer))->Name));
937 : }
938 : }
939 0 : ShowContinueError(state,
940 : "The Construction report will be produced. This will show more details on Constructions and their materials.");
941 0 : ShowContinueError(state, "Attempts will be made to complete the CTF process but the report may be incomplete.");
942 0 : ShowContinueError(state, "Constructs reported after this construction may appear to have all 0 CTFs.");
943 0 : ShowContinueError(state, "The potential causes of this problem are related to the input for the construction");
944 0 : ShowContinueError(state, "listed in the severe error above. The CTF calculate routine is unable to come up");
945 0 : ShowContinueError(state, "with a series of CTF terms that have a reasonable time step and this indicates an");
946 0 : ShowContinueError(state, "error. Check the definition of this construction and the materials that make up");
947 0 : ShowContinueError(state, "the construction. Very thin, highly conductive materials may cause problems.");
948 0 : ShowContinueError(state, "This may be avoided by ignoring the presence of those materials since they probably");
949 0 : ShowContinueError(state, "do not effect the heat transfer characteristics of the construction. Highly");
950 0 : ShowContinueError(state, "conductive or highly resistive layers that are alternated with high mass layers");
951 0 : ShowContinueError(state, "may also result in problems. After confirming that the input is correct and");
952 0 : ShowContinueError(state, "realistic, the user should contact the EnergyPlus support team.");
953 0 : DoCTFErrorReport = true;
954 0 : ErrorsFound = true;
955 0 : break;
956 : // CALL ShowFatalError(state, 'Program terminated for reasons listed (InitConductionTransferFunctions) ')
957 : }
958 :
959 : } // ... end of CTF calculation loop.
960 :
961 : } // ... end of IF block for non-reversed constructs.
962 :
963 : } else { // Construct has only resistive layers (no thermal mass).
964 : // CTF calculation not necessary, overall resistance
965 : // (R-value) is all that is needed.
966 :
967 : // Set time step for construct to user time step and the number of
968 : // inter-time step interpolations to 1
969 97 : this->CTFTimeStep = state.dataGlobal->TimeStepZone;
970 97 : this->NumHistories = 1;
971 97 : this->NumCTFTerms = 1;
972 :
973 97 : this->s0(1, 1) = cnd; // CTFs for current time
974 97 : this->s0(2, 1) = -cnd; // step are set to the
975 97 : this->s0(1, 2) = cnd; // overall conductance
976 97 : this->s0(2, 2) = -cnd; // of the construction.
977 :
978 97 : this->e.allocate(1);
979 97 : this->e = 0.0;
980 97 : this->s.allocate(2, 2, 1);
981 97 : this->s = 0.0;
982 97 : this->s(1, 1, 1) = 0.0; // CTF temperature
983 97 : this->s(2, 1, 1) = 0.0; // and flux
984 97 : this->s(1, 2, 1) = 0.0; // history terms
985 97 : this->s(2, 2, 1) = 0.0; // are all
986 97 : this->e(1) = 0.0; // zero.
987 :
988 97 : if (this->SourceSinkPresent) {
989 0 : ShowSevereError(state, format("Sources/sinks not allowed in purely resistive constructions --> {}", this->Name));
990 0 : ErrorsFound = true;
991 : }
992 :
993 97 : RevConst = false; // In the code that follows, handle a resistive
994 : // layer as a non-reversed this->
995 :
996 : } // ... end of resistive construction IF block.
997 :
998 : // Transfer the CTFs to the storage arrays for all non-reversed
999 : // constructions. This transfer was done earlier in the routine for
1000 : // reversed constructions.
1001 :
1002 246 : if (!RevConst) { // If this is either a new construction or a non-
1003 : // reversed construction, the CTFs must be stored
1004 : // in the proper arrays. If this is a reversed
1005 : // construction, nothing further needs to be done.
1006 :
1007 : // Copy the CTFs into the storage arrays, converting them back to SI
1008 : // units in the process. First the "zero" terms and then the history terms...
1009 243 : this->CTFOutside[0] = this->s0(1, 1) * DataConversions::CFU;
1010 243 : this->CTFCross[0] = this->s0(1, 2) * DataConversions::CFU;
1011 243 : this->CTFInside[0] = -this->s0(2, 2) * DataConversions::CFU;
1012 243 : if (this->SourceSinkPresent) {
1013 : // QTFs...
1014 0 : this->CTFSourceOut[0] = this->s0(3, 1);
1015 0 : this->CTFSourceIn[0] = this->s0(3, 2);
1016 : // QTFs for temperature calculation at source/sink location
1017 0 : this->CTFTSourceOut[0] = this->s0(1, 3);
1018 0 : this->CTFTSourceIn[0] = this->s0(2, 3);
1019 0 : this->CTFTSourceQ[0] = this->s0(3, 3) / DataConversions::CFU;
1020 0 : if (this->TempAfterLayer != 0) {
1021 : // QTFs for user specified interior temperature calculations...
1022 0 : this->CTFTUserOut[0] = this->s0(1, 4);
1023 0 : this->CTFTUserIn[0] = this->s0(2, 4);
1024 0 : this->CTFTUserSource[0] = this->s0(3, 4) / DataConversions::CFU;
1025 : }
1026 : }
1027 :
1028 1251 : for (int HistTerm = 1; HistTerm <= this->NumCTFTerms; ++HistTerm) {
1029 : // "REGULAR" CTFs...
1030 1008 : this->CTFOutside[HistTerm] = this->s(1, 1, HistTerm) * DataConversions::CFU;
1031 1008 : this->CTFCross[HistTerm] = this->s(1, 2, HistTerm) * DataConversions::CFU;
1032 1008 : this->CTFInside[HistTerm] = -this->s(2, 2, HistTerm) * DataConversions::CFU;
1033 1008 : if (HistTerm != 0) this->CTFFlux[HistTerm] = -e(HistTerm);
1034 1008 : if (this->SourceSinkPresent) {
1035 : // QTFs...
1036 0 : this->CTFSourceOut[HistTerm] = this->s(3, 1, HistTerm);
1037 0 : this->CTFSourceIn[HistTerm] = this->s(3, 2, HistTerm);
1038 : // QTFs for temperature calculation at source/sink location
1039 0 : this->CTFTSourceOut[HistTerm] = this->s(1, 3, HistTerm);
1040 0 : this->CTFTSourceIn[HistTerm] = this->s(2, 3, HistTerm);
1041 0 : this->CTFTSourceQ[HistTerm] = this->s(3, 3, HistTerm) / DataConversions::CFU;
1042 0 : if (this->TempAfterLayer != 0) {
1043 : // QTFs for user specified interior temperature calculations...
1044 0 : this->CTFTUserOut[HistTerm] = this->s(1, 4, HistTerm);
1045 0 : this->CTFTUserIn[HistTerm] = this->s(2, 4, HistTerm);
1046 0 : this->CTFTUserSource[HistTerm] = this->s(3, 4, HistTerm) / DataConversions::CFU;
1047 : }
1048 : }
1049 : }
1050 :
1051 : } // ... end of the reversed construction IF block.
1052 :
1053 246 : this->UValue = cnd * DataConversions::CFU;
1054 :
1055 246 : if (allocated(this->AExp)) this->AExp.deallocate();
1056 246 : if (allocated(this->AMat)) AMat.deallocate();
1057 246 : if (allocated(this->AInv)) this->AInv.deallocate();
1058 246 : if (allocated(this->IdenMatrix)) this->IdenMatrix.deallocate();
1059 246 : if (allocated(this->e)) this->e.deallocate();
1060 246 : if (allocated(this->Gamma1)) this->Gamma1.deallocate();
1061 246 : if (allocated(this->Gamma2)) this->Gamma2.deallocate();
1062 246 : if (allocated(this->s)) this->s.deallocate();
1063 246 : }
1064 :
1065 146 : void ConstructionProps::calculateExponentialMatrix()
1066 : {
1067 :
1068 : // SUBROUTINE INFORMATION:
1069 : // AUTHOR Russ Taylor
1070 : // DATE WRITTEN June 1990
1071 : // MODIFIED Dec 1995, Apr 1996, RKS; June 2000 RKS
1072 : // RE-ENGINEERED June 1996, RKS; Nov 1999, LKL;
1073 :
1074 : // PURPOSE OF THIS SUBROUTINE:
1075 : // This subroutine computes the exponential matrix exp(AMat*delt) for
1076 : // use in the state space method for the calculation of CTFs.
1077 :
1078 : // METHODOLOGY EMPLOYED:
1079 : // Uses the method of Taylor expansion combined with scaling and
1080 : // squaring to most efficiently compute the exponential matrix. The
1081 : // steps in the procedure are outlined in Seem's dissertation in
1082 : // Appendix A, page 128. Exponential matrix multiplication modified
1083 : // to take advantage of the characteristic form of AMat. AMat starts
1084 : // out as a tri-diagonal matrix. Each time AMat is raised to a higher
1085 : // power two extra non-zero diagonals are added. ExponMatrix now
1086 : // recognizes this. This should speed up the calcs somewhat. Also, a
1087 : // new cut-off criteria based on the significant figures of double-
1088 : // precision variables has been added. The main loop for higher powers
1089 : // of AMat is now stopped whenever these powers of AMat will no longer
1090 : // add to the summation (AExp) instead ofstopping potentially at the
1091 : // artifical limit of AMat**100.
1092 :
1093 : // REFERENCES:
1094 : // Seem, J.E. "Modeling of Heat Transfer in Buildings",
1095 : // Department of Mechanical Engineering, University of
1096 : // Wisconsin-Madison, 1987.
1097 : // Strand, R.K. "Testing Design Description for the CTF
1098 : // Calculation Code in BEST", BSO internal document,
1099 : // May/June 1996.
1100 :
1101 146 : constexpr Real64 DPLimit(1.0e-20);
1102 : // This argument is nice, but not sure it's accurate -- LKL Nov 1999.
1103 : // Parameter set to the significant figures limit of double
1104 : // precision variables plus a safety factor.- The argument for setting this parameter to 1E-20 involves the
1105 : // number of significant figures for REAL(r64) variables which is 16 and the largest power to which
1106 : // AMat will be raised which is 100. This would be a factor of 1E-18. A factor of "safety" of another 100
1107 : // arrives at the value chosen. It is argued that if one number is 1E-16 larger than a second number, then
1108 : // adding the second to the first will not effect the first. However, on the conservative side, there could
1109 : // be up to 100 numbers which might, added together, still could effect the original number. Each
1110 : // successive power of AMat will have terms smaller than the previous power. Thus, when the ratio between
1111 : // the terms of the latest power of AMat and the total (AExp) is less than DPLim, all further powers of
1112 : // AMat will have absolutely no effect on the REAL(r64) value of AExp. Thus, there is no need to
1113 : // continue the calculation. In effect, AExp has "converged". In REAL(r64)ity, 1E-16 would probably guarantee
1114 : // convergence since AMat terms drop off quickly, but the extra powers allows for differences between
1115 : // computer platforms.
1116 :
1117 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
1118 : Real64 AMatRowNorm; // Row norm for AMat
1119 : Real64 AMatRowNormMax; // Largest row norm for AMat
1120 146 : Array2D<Real64> AMat1; // AMat factored by (delt/2^k)
1121 146 : Array2D<Real64> AMato; // AMat raised to the previous power (power of AMat1-1)
1122 146 : Array2D<Real64> AMatN; // Current value of AMat raised to power n (n = 1,2...)
1123 : Real64 CheckVal; // Used to avoid possible overflow from Double->REAL(r64)->Integer
1124 : Real64 fact; // Intermediate calculation variable (delt/2^k)
1125 : int i; // Loop counter
1126 : int ic; // Loop counter
1127 : int ict; // Loop counter
1128 : int idm; // Loop counter
1129 : int ir; // Loop counter
1130 : int isq; // Loop counter
1131 : int j; // Loop counter
1132 : int k; // Power of 2 which is used to factor AMat
1133 : int l; // Theoretical power to which the A matrix must be
1134 : // raised to accurately calculate the exponential matrix
1135 : // when exponential calculation loop can be exited (i.e.
1136 : // the significant figure limit for REAL(r64)
1137 : // variables reached)
1138 :
1139 146 : AMat1.allocate(this->rcmax, this->rcmax);
1140 146 : AMato.allocate(this->rcmax, this->rcmax);
1141 146 : AMatN.allocate(this->rcmax, this->rcmax);
1142 :
1143 : // Subroutine initializations. AMat is assigned to local variable AMat1 to
1144 : // avoid the corruption of the original AMat 2-d array.
1145 146 : AMat1 = AMat;
1146 :
1147 : // Other arrays are initialized to zero.
1148 146 : this->AExp = 0.0;
1149 146 : AMato = 0.0;
1150 146 : AMatN = 0.0;
1151 :
1152 : // Step 1, page 128 (Seem's thesis): Compute the matrix row norm.
1153 : // See equation (A.3) which states that the matrix row norm is the
1154 : // maximum summation of the elements in a row of AMat multiplied by
1155 : // the time step.
1156 :
1157 : // Note With change to row-major arrays "row" here now means "column"
1158 :
1159 146 : AMatRowNormMax = 0.0; // Start of Step 1 ...
1160 :
1161 1731 : for (i = 1; i <= this->rcmax; ++i) {
1162 :
1163 1585 : AMatRowNorm = 0.0;
1164 24582 : for (j = 1; j <= this->rcmax; ++j) {
1165 22997 : AMatRowNorm += std::abs(AMat1(j, i));
1166 : }
1167 :
1168 1585 : AMatRowNorm *= this->CTFTimeStep;
1169 :
1170 1585 : AMatRowNormMax = max(AMatRowNormMax, AMatRowNorm);
1171 :
1172 : } // ... end of Step 1.
1173 :
1174 : // Step 2, page 128: Find smallest integer k such that
1175 : // AMatRowNormMax< = 2^k
1176 :
1177 146 : k = int(std::log(AMatRowNormMax) / std::log(2.0)) + 1; // Autodesk:Num Handle AMatRowNormMax=0
1178 :
1179 : // Step 3, page 128: Divide (AMat*delt) by 2^k. This section of code
1180 : // takes advantage of the fact that AMat is tridiagonal. Thus, it
1181 : // only factors the elements of the AMat that are known to be non-zero.
1182 :
1183 146 : fact = this->CTFTimeStep / std::pow(2.0, k); // Start of Step 3 ...
1184 146 : AMat1 *= fact; // ... end of Step 3.
1185 :
1186 : // Step 4, page 128: Calculate l, the highest power to which AMat
1187 : // must be taken theoretically to accurately calculate its exponential.
1188 : // This is based on a paper by Cadzow and Martens ("Discrete-Time and
1189 : // Computer Control Systems",Prentice-Hall, pp. 389-390, 1970). This
1190 : // number is now used as the maximum power to which AMat must be
1191 : // raised in order to calculate the exponential matrix. A new cut-off
1192 : // criteria based on the number of significant figures in a double-
1193 : // precision variable is used as a more practical limit on the
1194 : // exponentiation algorithm. One thing to note is that the EnergyPlus
1195 : // code is slightly different here than what is presented in Seem's
1196 : // dissertation. In Seem's dissertation, AMatRowNormMax (delta, the
1197 : // timestep, is already factored into this term above) is divided by
1198 : // 2^k with k defined in the code above. Dividing AMatRowNormMax by
1199 : // 2^k would have the net effect of decreasing the value of
1200 : // AMatRowNormMax and thus potentially CheckVal as well. l, the integer
1201 : // value of CheckVal, is used to define the number of terms in the
1202 : // exponential matrix of AMat that are required for an accurate estimation
1203 : // of that matrix. In practice, additional terms probably won't have
1204 : // a large effect but in general should make the calculation MORE accurate.
1205 : // Also, if the new terms are too small, then as noted above the cut-off
1206 : // criteria will not use them.
1207 :
1208 146 : CheckVal = min(3.0 * AMatRowNormMax + 6.0, 100.0);
1209 146 : l = int(CheckVal);
1210 :
1211 : // Step 5, page 128: Calculate the exponential. First, add the
1212 : // linear term to the identity matrix.
1213 146 : this->AExp = AMat1 + this->IdenMatrix; // Start of Step 5 ...
1214 :
1215 : // Now, add successive terms to the expansion as per the standard
1216 : // exponential formula. AMato contains the last "power" of AMat
1217 : // which saves the program from having to remultiply the entire power
1218 : // of AMat each time. Since this is still the linear power of AMat,
1219 : // AMat1 is still tridiagonal in nature.
1220 146 : AMato = AMat1;
1221 :
1222 146 : i = 1; // Initialize the counter for the following DO loop
1223 :
1224 : // The following DO WHILE loop continues to raise AMat to successive
1225 : // powers and add it to the exponential matrix (AExp).
1226 8935 : while (i < l) { // Begin power raising loop ...
1227 :
1228 8789 : ++i; // Increment the loop counter
1229 8789 : bool SigFigLimit = true; // Set the significant factor limit flag
1230 :
1231 135813 : for (ir = 1; ir <= this->rcmax; ++ir) { // Begin matrix multiplication loop ...
1232 : // The following matrix multiplication could be "optimized" since
1233 : // for one-dimensional heat transfer AMat is 3-diagonal, AMat squared
1234 : // is 5-diagonal, etc. However, the code can be much simpler if we
1235 : // ignore this fact and just do a generic matrix multiplication.
1236 : // For 2-D heat transfer, the number of off-diagonal non-zero terms
1237 : // is slightly more complicated as well.
1238 2239182 : for (ic = 1; ic <= this->rcmax; ++ic) {
1239 2112158 : AMatN(ic, ir) = 0.0;
1240 41050656 : for (ict = 1; ict <= this->rcmax; ++ict) {
1241 : // Make sure the next term won't cause an underflow. If it will end up being
1242 : // so small as to go below TinyLimit, then ignore it since it won't add anything
1243 : // to AMatN anyway.
1244 38938498 : if (std::abs(AMat1(ic, ict)) > Constant::rTinyValue) {
1245 6082426 : if (std::abs(AMato(ict, ir)) > std::abs(double(i) * Constant::rTinyValue / AMat1(ic, ict)))
1246 282571 : AMatN(ic, ir) += AMato(ict, ir) * AMat1(ic, ict) / double(i);
1247 : }
1248 : }
1249 : }
1250 : } // ... end of matrix multiplication loop.
1251 :
1252 : // Update AMato and AExp matrices
1253 8789 : AMato = AMatN;
1254 8789 : this->AExp += AMato;
1255 :
1256 : // The next DO loop tests the significant figures limit criteria to
1257 : // see if any values in AExp are still changing appreciably.
1258 9152 : for (ir = 1; ir <= this->rcmax; ++ir) {
1259 60671 : for (ic = 1; ic <= this->rcmax; ++ic) {
1260 : // Test of limit criteria:
1261 60308 : if (std::abs(this->AExp(ic, ir)) > Constant::rTinyValue) { // Next line divides by AExp entry so it
1262 : // must be checked to avoid dividing by zero.
1263 : // If the ratio between any current element in the power
1264 : // of AMat and its corresponding element in AExp is
1265 : // greater than the number which might effect the overall
1266 : // exponential matrix based on stability criteria, then
1267 : // continue raising AMat to another power (SigFigLimit = false).
1268 :
1269 53066 : if (std::abs(AMato(ic, ir) / this->AExp(ic, ir)) > DPLimit) {
1270 1489 : SigFigLimit = false;
1271 1489 : break; // DO loop (anytime SigFigLimit is false, AMat must continue to be raised another power)
1272 : }
1273 :
1274 : } else { // There are still elements of AExp which are zero, so
1275 : // the raising of AMat to higher powers should continue.
1276 :
1277 7242 : SigFigLimit = false;
1278 7242 : break; // DO loop (anytime SigFigLimit is false, AMat must continue to be raised another power)
1279 : }
1280 : }
1281 9094 : if (!SigFigLimit) break; // DO loop (anytime SigFigLimit is false, AMat must continue to be raised another power)
1282 : }
1283 :
1284 : // Compute next term, only if necessary. If SigFigLimit is still true,
1285 : // then all of the new terms being added to AExp are too small to
1286 : // affect it. Thus, there is no need to continue this do loop further.
1287 :
1288 8789 : if (SigFigLimit) i = 100; // SigFigLimit is still true, set i to maximum possible
1289 : // value of l (100).
1290 :
1291 : } // ... end of power raising loop and Step 5.
1292 :
1293 : // Step 6, page 128:
1294 : // Square AExp "k times" to obtain the actual exponential matrix
1295 : // (remember that AExp was scaled earlier in this routine).
1296 :
1297 1459 : for (isq = 1; isq <= k; ++isq) { // Begin squaring DO loop and Step 6 ...
1298 :
1299 : // Use AMato to store the old values of AExp
1300 1314 : AMato = this->AExp;
1301 1314 : bool Backup = true; // Used when numerics get to small in Exponentiation
1302 1314 : this->AExp = 0.0;
1303 :
1304 : // Multiply the old value of AExp (AMato) by itself and store in AExp.
1305 17927 : for (ir = 1; ir <= this->rcmax; ++ir) {
1306 268680 : for (ic = 1; ic <= this->rcmax; ++ic) {
1307 4644252 : for (idm = 1; idm <= this->rcmax; ++idm) {
1308 4392185 : if (std::abs(AMato(idm, ir) * AMato(ic, idm)) > Constant::rTinyValue) {
1309 2168531 : this->AExp(ic, ir) += AMato(idm, ir) * AMato(ic, idm);
1310 2168531 : Backup = false;
1311 : }
1312 : }
1313 : }
1314 : }
1315 : // Backup is true when every item of AExp didnt pass the TinyLimit test
1316 1314 : if (Backup) {
1317 1 : this->AExp = AMato;
1318 1 : break;
1319 : }
1320 :
1321 : } // ... end of squaring loop and Step 6.
1322 :
1323 146 : AMat1.deallocate();
1324 146 : AMato.deallocate();
1325 146 : AMatN.deallocate();
1326 146 : }
1327 :
1328 146 : void ConstructionProps::calculateInverseMatrix()
1329 : {
1330 :
1331 : // SUBROUTINE INFORMATION:
1332 : // AUTHOR Rick Strand
1333 : // DATE WRITTEN Dec 1995
1334 : // MODIFIED June 2000 RKS (made routine generic to allow for 2-D solutions)
1335 : // RE-ENGINEERED June 1996, February 1997 RKS
1336 :
1337 : // PURPOSE OF THIS SUBROUTINE:
1338 : // This subroutine computes the inverse of AMat for use
1339 : // in the calculation of the CTFs.
1340 :
1341 : // METHODOLOGY EMPLOYED:
1342 : // Uses row elimination to zero the off-diagonal terms of
1343 : // AMat while performing the same operations on another
1344 : // matrix which starts as the identity matrix. Once AMat
1345 : // has been converted to an identity matrix(I), the other
1346 : // matrix which started as the I will then be the inverse
1347 : // of A. This algorithm has been customized for a
1348 : // tri-diagonal matrix.
1349 :
1350 : // REFERENCES:
1351 : // Any linear algebra test (this is a generic routine).
1352 :
1353 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
1354 146 : Array2D<Real64> AMat1; // Intermediate calculation matrix equivalent at first to AMat
1355 :
1356 : // Subroutine initializations ...
1357 146 : AMat1.allocate(this->rcmax, this->rcmax);
1358 :
1359 146 : AMat1 = AMat; // Set AMat1 = AMat to avoid AMat changes
1360 146 : this->AInv = this->IdenMatrix; // Set AInv to Identity Matrix
1361 :
1362 : // Use Gaussian elimination to zero all of the elements of AMat left
1363 : // of the diagonal.
1364 : // This DO loop will cycle through each of the rows in AMat except the
1365 : // last row which is handled later because it does not have to be used
1366 : // to eliminate any other rows. The index ir is the current row
1367 : // number and also the column of the current diagonal element.
1368 :
1369 1585 : for (int ir = 1; ir <= this->rcmax - 1; ++ir) { // Begin forward elimination loop ...
1370 :
1371 : // Factor all of the elements of the row being used to zero the next
1372 : // row in both AMat and AInv by the diagonal element of this row.
1373 : // We should only need to factor the elements to the right of the
1374 : // diagonal since those to the right of it should be zero.
1375 12145 : for (int ic = ir + 1; ic <= this->rcmax; ++ic) {
1376 10706 : AMat1(ic, ir) /= AMat1(ir, ir);
1377 : }
1378 :
1379 : // In the forward elimination process, all the elements in AInv to the
1380 : // right of the diagonal are zero so they do not need to be factored.
1381 12145 : for (int ic = 1; ic <= ir; ++ic) {
1382 10706 : this->AInv(ic, ir) /= AMat1(ir, ir);
1383 : }
1384 :
1385 1439 : AMat1(ir, ir) = 1.0; // By definition, the diagonal of AMat is now 1.
1386 :
1387 : // Use this factored row to eliminate the off-diagonal element of the
1388 : // rows below the current one (ir)...
1389 :
1390 12145 : for (int irr = ir + 1; irr <= this->rcmax; ++irr) { // Start of row reduction loop...
1391 :
1392 133842 : for (int ic = ir + 1; ic <= this->rcmax; ++ic) {
1393 123136 : AMat1(ic, irr) -= AMat1(ir, irr) * AMat1(ic, ir);
1394 : }
1395 :
1396 : // Now, determine the effect on the next row of AInv. Again, all of
1397 : // the elements in AInv to the right of the diagonal are zero, so they
1398 : // can be ignored.
1399 :
1400 77627 : for (int ic = 1; ic <= ir; ++ic) {
1401 66921 : this->AInv(ic, irr) -= AMat1(ir, irr) * this->AInv(ic, ir);
1402 : }
1403 :
1404 10706 : AMat1(ir, irr) = 0.0; // By definition, the element to the left of the
1405 : // diagonal in the next row of AMat is now zero.
1406 :
1407 : } // ...end of row reduction loop
1408 :
1409 : } // ... end of the forward elimination loop.
1410 :
1411 : // Factor the last row of AInv by the current value of the last
1412 : // diagonal element of AMat. After this is done, all of the diagonal
1413 : // elements of AMat are unity and all of the elements in AMat left of
1414 : // the diagonal are zero.
1415 :
1416 1731 : for (int ic = 1; ic <= this->rcmax; ++ic) {
1417 1585 : this->AInv(ic, this->rcmax) /= AMat1(this->rcmax, this->rcmax);
1418 : }
1419 146 : AMat1(this->rcmax, this->rcmax) = 1.0;
1420 :
1421 : // Now, use back substitution to eliminate the elements to the right
1422 : // of the diagonal in AMat. The procedure is similar to the forward
1423 : // elimination process except that we only have to operate on AInv,
1424 : // though now all of the columns of AInv may be non-zero.
1425 :
1426 : // This DO loop will cycle through the remaining rows which are not
1427 : // yet diagonalized in reverse order. Note that the only effect on
1428 : // AMat is that the off-diagonal element is zeroed. The diagonal
1429 : // (which has already been set to unity) is not effected by this row
1430 : // elimination process.
1431 : // In the following code ir is the column being zeroed and irr is the
1432 : // row being worked on
1433 :
1434 1585 : for (int ir = this->rcmax; ir >= 2; --ir) { // Begin reverse elimination loop ...
1435 12145 : for (int irr = 1; irr <= ir - 1; ++irr) {
1436 200763 : for (int ic = 1; ic <= this->rcmax; ++ic) {
1437 190057 : this->AInv(ic, irr) -= AMat1(ir, irr) * this->AInv(ic, ir);
1438 : }
1439 10706 : AMat1(ir, irr) = 0.0;
1440 : }
1441 : } // ... end of reverse elimination loop.
1442 :
1443 : // At this point, AMat1 is equal to the identity matrix (I)
1444 : // and AInv is equal to the inverse of AMat.
1445 :
1446 146 : AMat1.deallocate();
1447 146 : }
1448 :
1449 146 : void ConstructionProps::calculateGammas()
1450 : {
1451 :
1452 : // SUBROUTINE INFORMATION:
1453 : // AUTHOR Russ Taylor
1454 : // DATE WRITTEN June 1990
1455 : // RE-ENGINEERED July 1996, RKS
1456 :
1457 : // PURPOSE OF THIS SUBROUTINE:
1458 : // Compute gammas as defined in Seem's dissertation.
1459 : // Runs as a subroutine of the conduction transfer
1460 : // function solver (InitializeCTFs).
1461 :
1462 : // METHODOLOGY EMPLOYED:
1463 : // Determine the Gamma1 and Gamma2 based on the results
1464 : // from the ExponMatrix and InvertMatrix subroutines.
1465 : // This routine is specialized to take advantage of the
1466 : // fact that most of BMat consists of zeroes.
1467 :
1468 : // REFERENCES:
1469 : // The state space method of calculating CTFs is
1470 : // outlined in the doctoral dissertation of John Seem,
1471 : // "Modeling of Heat Transfer in Buildings", Department
1472 : // of Mechanical Engineering, University of Wisconsin-
1473 : // Madison, 1987.
1474 :
1475 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
1476 146 : Array2D<Real64> ATemp; // Intermediate variable equal to AExp - I
1477 :
1478 : // Compute Gamma1 from equation (2.1.12) in Seem's dissertation which
1479 : // states that: Gamma1 = [AInv] * ([AExp]-[I]) * [BMat]
1480 : // noting that BMat contains only the non-zero values of the B Matrix.
1481 :
1482 146 : ATemp.allocate(this->rcmax, this->rcmax);
1483 146 : ATemp = this->AExp - this->IdenMatrix;
1484 146 : Gamma1 = 0.0;
1485 :
1486 1731 : for (int i = 1; i <= this->rcmax; ++i) {
1487 :
1488 24582 : for (int is1 = 1; is1 <= this->rcmax; ++is1) {
1489 :
1490 22997 : if (this->SolutionDimensions == 1) {
1491 22997 : this->Gamma1(1, i) += this->AInv(is1, i) * ATemp(1, is1) * this->BMat(1);
1492 22997 : this->Gamma1(2, i) += this->AInv(is1, i) * ATemp(this->rcmax, is1) * this->BMat(2);
1493 : } else { // SolutionDimensions = 2
1494 0 : for (int SurfNode = 1; SurfNode <= this->NumOfPerpendNodes; ++SurfNode) {
1495 0 : this->Gamma1(1, i) += this->AInv(is1, i) * ATemp(SurfNode, is1) * this->BMat(1);
1496 0 : this->Gamma1(2, i) += this->AInv(is1, i) * ATemp(this->rcmax + 1 - SurfNode, is1) * this->BMat(2);
1497 : }
1498 : }
1499 :
1500 22997 : if (this->NodeSource > 0) {
1501 0 : this->Gamma1(3, i) += this->AInv(is1, i) * ATemp(this->NodeSource, is1) * this->BMat(3);
1502 : }
1503 : }
1504 : }
1505 :
1506 146 : ATemp.deallocate();
1507 : // Compute Gamma2 from equation (2.1.13) in Seem's dissertation which
1508 : // states that: Gamma2 = [AInv] * ([Gamma1]/delt - [BMat])
1509 : // again noting that BMat contains only the non-zero values of B.
1510 146 : Gamma2 = 0.0;
1511 :
1512 1731 : for (int i = 1; i <= this->rcmax; ++i) {
1513 :
1514 6340 : for (int j = 1; j <= 3; ++j) {
1515 :
1516 73746 : for (int is1 = 1; is1 <= this->rcmax; ++is1) {
1517 :
1518 68991 : if (this->SolutionDimensions == 1) {
1519 68991 : if ((j == 1) && (is1 == 1)) {
1520 1585 : this->Gamma2(j, i) += this->AInv(is1, i) * (this->Gamma1(j, is1) / this->CTFTimeStep - this->BMat(1));
1521 67406 : } else if ((j == 2) && (is1 == this->rcmax)) {
1522 1585 : this->Gamma2(j, i) += this->AInv(is1, i) * (this->Gamma1(j, is1) / this->CTFTimeStep - this->BMat(2));
1523 65821 : } else if ((j == 3) && (is1 == this->NodeSource)) {
1524 0 : this->Gamma2(j, i) += this->AInv(is1, i) * (this->Gamma1(j, is1) / this->CTFTimeStep - this->BMat(3));
1525 : } else { // the element of the actual BMat is zero
1526 65821 : this->Gamma2(j, i) += this->AInv(is1, i) * (this->Gamma1(j, is1) / this->CTFTimeStep);
1527 : }
1528 : } else { // SolutionDimensions = 2
1529 0 : if ((j == 1) && ((is1 >= 1) && (is1 <= this->NumOfPerpendNodes))) {
1530 0 : this->Gamma2(j, i) += this->AInv(is1, i) * (this->Gamma1(j, is1) / this->CTFTimeStep - this->BMat(1));
1531 0 : } else if ((j == 2) && ((is1 <= this->rcmax) && (is1 >= this->rcmax + 1 - this->NumOfPerpendNodes))) {
1532 0 : this->Gamma2(j, i) += this->AInv(is1, i) * (this->Gamma1(j, is1) / this->CTFTimeStep - this->BMat(2));
1533 0 : } else if ((j == 3) && (is1 == this->NodeSource)) {
1534 0 : this->Gamma2(j, i) += this->AInv(is1, i) * (this->Gamma1(j, is1) / this->CTFTimeStep - this->BMat(3));
1535 : } else { // the element of the actual BMat is zero
1536 0 : this->Gamma2(j, i) += this->AInv(is1, i) * (this->Gamma1(j, is1) / this->CTFTimeStep);
1537 : }
1538 : }
1539 : }
1540 : }
1541 : }
1542 146 : }
1543 :
1544 146 : void ConstructionProps::calculateFinalCoefficients()
1545 : {
1546 :
1547 : // SUBROUTINE INFORMATION:
1548 : // AUTHOR Russ Taylor
1549 : // DATE WRITTEN June 1990
1550 : // MODIFIED Apr 96, RKS, cosmetic, algorithm neutral changes
1551 : // RE-ENGINEERED July 1996, RKS; Nov 1999, LKL (allocatable arrays)
1552 :
1553 : // PURPOSE OF THIS SUBROUTINE:
1554 : // Subprogram to calculate the Sj and ej coefficients of Seem's
1555 : // dissertation. Follows Seem's technique to compute the coefficients
1556 : // in order with minimum storage requirements.
1557 :
1558 : // METHODOLOGY EMPLOYED:
1559 : // Combine the results of the ExponMatrix, InvertMatrix, and
1560 : // CalculateGammas routines together to arrive at the temperature
1561 : // coefficients (s, s0) and the heat flux history coefficients (e) of
1562 : // the CTFs. The outline of this subroutine is based on step 5 of
1563 : // Seem's suggested implementation of the state space method found on
1564 : // pages 26+27 of his dissertation.
1565 :
1566 : // REFERENCES:
1567 : // The state space method of calculating CTFs is outlined in the
1568 : // doctoral dissertation of John Seem, "Modeling of Heat Transfer in
1569 : // Buildings", Department of Mechanical Engineering, University of
1570 : // Wisconsin-Madison, 1987. In particular, the equations used for
1571 : // these calculations are equations (2.1.24) through (2.1.26) in Seem's
1572 : // dissertation.
1573 :
1574 : // SUBROUTINE PARAMETER DEFINITIONS:
1575 146 : constexpr Real64 ConvrgLim(1.0e-13); // Convergence limit (ratio) for cutting off the calculation of further
1576 : // CTFs. This value was found to give suitable accuracy in IBLAST.
1577 :
1578 : Real64 avg; // Intermediate calculation variable (average)
1579 : bool CTFConvrg; // Set after CTFs are calculated, based on whether there are too many CTFs terms
1580 146 : Array2D<Real64> PhiR0; // Product of Phi( = AExp) and R0 matrices from the state space method
1581 : Real64 rat; // Intermediate calculation variable (ratio of flux history terms)
1582 146 : Array2D<Real64> Rnew; // Current R matrix
1583 146 : Array2D<Real64> Rold; // R matrix from the last iteration
1584 : Real64 SurfNodeFac; // Multiplying factor applied to various surface nodes
1585 : Real64 trace; // Trace of the product of Phi( = AExp) and R0
1586 :
1587 : // Subroutine initializations
1588 146 : PhiR0.allocate(this->rcmax, this->rcmax);
1589 146 : Rnew.allocate(this->rcmax, this->rcmax);
1590 146 : Rold.allocate(this->rcmax, this->rcmax);
1591 146 : PhiR0 = 0.0;
1592 146 : Rold = 0.0;
1593 :
1594 146 : this->s0 = 0.0;
1595 146 : this->s = 0.0;
1596 146 : this->e = 0.0;
1597 146 : Rnew = this->IdenMatrix; // Rnew initialized to the identity matrix
1598 :
1599 : // Calculate Gamma1-Gamma2. Gamma1 is not used by itself in the
1600 : // equations, only Gamma1-Gamma2. Thus, reset Gamma1 to:
1601 : // Gamma1-Gamma2
1602 1731 : for (int i = 1; i <= this->rcmax; ++i) {
1603 6340 : for (int j = 1; j <= 3; ++j) {
1604 4755 : this->Gamma1(j, i) -= this->Gamma2(j, i);
1605 : }
1606 : }
1607 :
1608 : // Compute s0. See Seem's thesis equation (2.1.24) which states that:
1609 : // s0 = (CMat*R0*Gamma2) + (DMat)
1610 : // Note that for a two-dimensional solution, there is more than one
1611 : // node at the surface and the effect of each of these must be added
1612 : // together.
1613 146 : if (this->SolutionDimensions == 1) {
1614 146 : this->s0(1, 1) = this->CMat(1) * this->Gamma2(1, 1) + this->DMat(1);
1615 146 : this->s0(2, 1) = this->CMat(1) * this->Gamma2(2, 1);
1616 146 : this->s0(3, 1) = this->CMat(1) * this->Gamma2(3, 1);
1617 146 : this->s0(1, 2) = this->CMat(2) * this->Gamma2(1, this->rcmax);
1618 146 : this->s0(2, 2) = this->CMat(2) * this->Gamma2(2, this->rcmax) + this->DMat(2);
1619 146 : this->s0(3, 2) = this->CMat(2) * this->Gamma2(3, this->rcmax);
1620 : } else { // SolutionDimensions = 2
1621 0 : for (int SurfNode = 1; SurfNode <= this->NumOfPerpendNodes; ++SurfNode) {
1622 0 : if ((SurfNode == 1) || (SurfNode == this->NumOfPerpendNodes)) {
1623 0 : SurfNodeFac = 0.5;
1624 : } else {
1625 0 : SurfNodeFac = 1.0;
1626 : }
1627 0 : this->s0(1, 1) += SurfNodeFac * this->CMat(1) * this->Gamma2(1, SurfNode);
1628 0 : this->s0(2, 1) += SurfNodeFac * this->CMat(1) * this->Gamma2(2, SurfNode);
1629 0 : this->s0(3, 1) += SurfNodeFac * this->CMat(1) * this->Gamma2(3, SurfNode);
1630 0 : this->s0(1, 2) += SurfNodeFac * this->CMat(2) * this->Gamma2(1, this->rcmax + 1 - SurfNode);
1631 0 : this->s0(2, 2) += SurfNodeFac * this->CMat(2) * this->Gamma2(2, this->rcmax + 1 - SurfNode);
1632 0 : this->s0(3, 2) += SurfNodeFac * this->CMat(2) * this->Gamma2(3, this->rcmax + 1 - SurfNode);
1633 : }
1634 0 : this->s0(1, 1) += double(this->NumOfPerpendNodes - 1) * this->DMat(1);
1635 0 : this->s0(2, 2) += double(this->NumOfPerpendNodes - 1) * this->DMat(2);
1636 : }
1637 :
1638 146 : if (this->NodeSource > 0) {
1639 0 : this->s0(1, 3) = this->Gamma2(1, this->NodeSource);
1640 0 : this->s0(2, 3) = this->Gamma2(2, this->NodeSource);
1641 0 : this->s0(3, 3) = this->Gamma2(3, this->NodeSource);
1642 : }
1643 146 : if (this->NodeUserTemp > 0) {
1644 0 : this->s0(1, 4) = this->Gamma2(1, this->NodeUserTemp);
1645 0 : this->s0(2, 4) = this->Gamma2(2, this->NodeUserTemp);
1646 0 : this->s0(3, 4) = this->Gamma2(3, this->NodeUserTemp);
1647 : }
1648 :
1649 : // Check for and enforce symmetry in the cross term (Y)
1650 146 : if (std::abs(this->s0(2, 1)) != std::abs(this->s0(1, 2))) {
1651 138 : avg = (std::abs(this->s0(2, 1)) + std::abs(this->s0(1, 2))) * 0.5;
1652 138 : this->s0(2, 1) *= avg / std::abs(this->s0(2, 1));
1653 138 : this->s0(1, 2) *= avg / std::abs(this->s0(1, 2));
1654 : }
1655 :
1656 : // Compute S's and e's from 1 to n-1. See equations (2.1.25) and
1657 : // (2.1.26) and Appendix C.
1658 146 : int inum = 1; // Set history term counter
1659 146 : CTFConvrg = false; // Set the convergence logical to false
1660 :
1661 : // The following DO WHILE loop calculates each successive set of time
1662 : // history terms until there are rcmax number of history terms or the
1663 : // latest flux history term is negligibly small compared to the first
1664 : // flux history term.
1665 1003 : while ((!CTFConvrg) && (inum < this->rcmax)) { // Begin CTF calculation loop ...
1666 :
1667 : // Compute e(inum) based on Appendix C (Seem's dissertation). First,
1668 : // compute the new PhiR0 and its trace.
1669 :
1670 857 : trace = 0.0;
1671 :
1672 11804 : for (int ir = 1; ir <= this->rcmax; ++ir) {
1673 :
1674 187140 : for (int ic = 1; ic <= this->rcmax; ++ic) {
1675 176193 : PhiR0(ic, ir) = 0.0;
1676 3433992 : for (int is = 1; is <= this->rcmax; ++is) {
1677 : // Make sure the next term won't cause an underflow. If it will end up being so small
1678 : // as to go below TinyLimit, then ignore it since it won't add anything to PhiR0 anyway.
1679 3257799 : if (std::abs(Rnew(ic, is)) > Constant::rTinyValue) {
1680 2788317 : if (std::abs(this->AExp(is, ir)) > std::abs(Constant::rTinyValue / Rnew(ic, is)))
1681 2430221 : PhiR0(ic, ir) += this->AExp(is, ir) * Rnew(ic, is);
1682 : }
1683 : }
1684 : }
1685 :
1686 10947 : trace += PhiR0(ir, ir);
1687 : }
1688 :
1689 : // Now calculate ej from the trace. According to Appendix C:
1690 : // e(j) = -Trace[AExp*R(j-1)]/j
1691 :
1692 857 : this->e(inum) = -trace / double(inum);
1693 :
1694 : // Update Rold and compute Rnew. Note: PhiR0 = AExp*R(j-1) here.
1695 : // According to Appendix C: R(j) = AExp*R(j-1) + e(j-1)
1696 :
1697 11804 : for (int ir = 1; ir <= this->rcmax; ++ir) {
1698 187140 : for (int ic = 1; ic <= this->rcmax; ++ic) {
1699 176193 : Rold(ic, ir) = Rnew(ic, ir);
1700 176193 : Rnew(ic, ir) = PhiR0(ic, ir);
1701 : }
1702 10947 : Rnew(ir, ir) += this->e(inum);
1703 : }
1704 :
1705 : // Compute S(inum) based on eq.(2.1.25) which states:
1706 : // S(j) = CMat*[R(j-1)*(Gamma1-Gamma2)+R(j)*Gamma2]
1707 : // + e(j)*DMat
1708 857 : if (this->SolutionDimensions == 1) {
1709 3428 : for (int j = 1; j <= 3; ++j) {
1710 35412 : for (int is2 = 1; is2 <= this->rcmax; ++is2) {
1711 32841 : this->s(j, 1, inum) += this->CMat(1) * (Rold(is2, 1) * this->Gamma1(j, is2) + Rnew(is2, 1) * this->Gamma2(j, is2));
1712 65682 : this->s(j, 2, inum) +=
1713 32841 : this->CMat(2) * (Rold(is2, this->rcmax) * this->Gamma1(j, is2) + Rnew(is2, this->rcmax) * this->Gamma2(j, is2));
1714 32841 : if (this->NodeSource > 0) {
1715 0 : this->s(j, 3, inum) +=
1716 0 : (Rold(is2, this->NodeSource) * this->Gamma1(j, is2) + Rnew(is2, this->NodeSource) * this->Gamma2(j, is2));
1717 : }
1718 32841 : if (this->NodeUserTemp > 0) {
1719 0 : this->s(j, 4, inum) +=
1720 0 : (Rold(is2, this->NodeUserTemp) * this->Gamma1(j, is2) + Rnew(is2, this->NodeUserTemp) * this->Gamma2(j, is2));
1721 : }
1722 : }
1723 2571 : if (j != 3) this->s(j, j, inum) += this->e(inum) * this->DMat(j);
1724 : }
1725 : } else { // SolutionDimensions = 2
1726 0 : for (int j = 1; j <= 3; ++j) {
1727 0 : for (int is2 = 1; is2 <= this->rcmax; ++is2) {
1728 0 : for (int SurfNode = 1; SurfNode <= this->NumOfPerpendNodes; ++SurfNode) {
1729 0 : if ((SurfNode == 1) || (SurfNode == this->NumOfPerpendNodes)) {
1730 0 : SurfNodeFac = 0.5;
1731 : } else {
1732 0 : SurfNodeFac = 1.0;
1733 : }
1734 0 : this->s(j, 1, inum) +=
1735 0 : SurfNodeFac * this->CMat(1) * (Rold(is2, SurfNode) * this->Gamma1(j, is2) + Rnew(is2, SurfNode) * this->Gamma2(j, is2));
1736 0 : this->s(j, 2, inum) += SurfNodeFac * this->CMat(2) *
1737 0 : (Rold(is2, this->rcmax + 1 - SurfNode) * this->Gamma1(j, is2) +
1738 0 : Rnew(is2, this->rcmax + 1 - SurfNode) * this->Gamma2(j, is2));
1739 : }
1740 0 : if (this->NodeSource > 0) {
1741 0 : this->s(j, 3, inum) +=
1742 0 : (Rold(is2, this->NodeSource) * this->Gamma1(j, is2) + Rnew(is2, this->NodeSource) * this->Gamma2(j, is2));
1743 : }
1744 0 : if (this->NodeUserTemp > 0) {
1745 0 : this->s(j, 4, inum) +=
1746 0 : (Rold(is2, this->NodeUserTemp) * this->Gamma1(j, is2) + Rnew(is2, this->NodeUserTemp) * this->Gamma2(j, is2));
1747 : }
1748 : }
1749 : }
1750 0 : this->s(1, 1, inum) += this->e(inum) * this->DMat(1) * double(this->NumOfPerpendNodes - 1);
1751 0 : this->s(2, 2, inum) += this->e(inum) * this->DMat(2) * double(this->NumOfPerpendNodes - 1);
1752 : }
1753 :
1754 : // Check for and enforce symmetry in the cross term (Y)
1755 857 : if (std::abs(s(2, 1, inum)) != std::abs(s(1, 2, inum))) {
1756 855 : avg = (std::abs(s(2, 1, inum)) + std::abs(s(1, 2, inum))) * 0.5;
1757 855 : this->s(2, 1, inum) *= avg / std::abs(s(2, 1, inum));
1758 855 : this->s(1, 2, inum) *= avg / std::abs(s(1, 2, inum));
1759 : }
1760 :
1761 : // Check for convergence of the CTFs.
1762 857 : if (e(1) == 0.0) {
1763 :
1764 0 : this->NumCTFTerms = 1; // e(1) is zero, so there are no history terms.
1765 0 : CTFConvrg = true; // CTF calculations have converged--set logical.
1766 :
1767 : } else {
1768 : // e(1) is non-zero -- Calculate and compare the ratio of the flux
1769 : // terms to the convergence limit.
1770 857 : rat = std::abs(e(inum) / this->e(1));
1771 :
1772 857 : if (rat < ConvrgLim) {
1773 :
1774 : // If the ratio is less than the convergence limit, then any other
1775 : // terms would have a neglible impact on the CTF-based energy balances.
1776 92 : this->NumCTFTerms = inum;
1777 92 : CTFConvrg = true; // CTF calculations have converged--set logical.
1778 : }
1779 : } // ... end of convergence check block.
1780 :
1781 857 : ++inum;
1782 :
1783 : } // ... end of CTF calculation loop.
1784 : // Continue to the next coefficient if the solution has not converged
1785 146 : if (!CTFConvrg) { // Compute last e and S, if still unconverged.
1786 :
1787 : // Compute e(inum) based on Appendix C (Seem's dissertation) or see
1788 : // equation above. First compute the new PhiR0 and its trace.
1789 :
1790 54 : trace = 0.0;
1791 :
1792 336 : for (int ir = 1; ir <= this->rcmax; ++ir) {
1793 1800 : for (int is = 1; is <= this->rcmax; ++is) {
1794 1518 : trace += this->AExp(is, ir) * Rnew(ir, is);
1795 : }
1796 : }
1797 :
1798 54 : this->e(this->rcmax) = -trace / double(this->rcmax); // Now calculate ej from the trace.
1799 :
1800 : // Compute S(inum) based on eq.(2.1.25) which states:
1801 : // S(last) = CMat*R(last-1)*(Gamma1-Gamma2)+e(last)*DMat
1802 :
1803 54 : if (this->SolutionDimensions == 1) {
1804 216 : for (int j = 1; j <= 3; ++j) {
1805 1008 : for (int is2 = 1; is2 <= this->rcmax; ++is2) {
1806 846 : this->s(j, 1, this->rcmax) += this->CMat(1) * Rnew(is2, 1) * this->Gamma1(j, is2);
1807 846 : this->s(j, 2, this->rcmax) += this->CMat(2) * Rnew(is2, this->rcmax) * this->Gamma1(j, is2);
1808 846 : if (this->NodeSource > 0) {
1809 0 : this->s(j, 3, this->rcmax) += Rnew(is2, this->NodeSource) * this->Gamma1(j, is2);
1810 : }
1811 846 : if (this->NodeUserTemp > 0) {
1812 0 : this->s(j, 4, this->rcmax) += Rnew(is2, this->NodeUserTemp) * this->Gamma1(j, is2);
1813 : }
1814 : }
1815 : }
1816 54 : this->s(1, 1, this->rcmax) += this->e(this->rcmax) * this->DMat(1);
1817 54 : this->s(2, 2, this->rcmax) += this->e(this->rcmax) * this->DMat(2);
1818 54 : this->NumCTFTerms = this->rcmax;
1819 : } else { // SolutionDimensions = 2
1820 0 : for (int j = 1; j <= 3; ++j) {
1821 0 : for (int is2 = 1; is2 <= this->rcmax; ++is2) {
1822 0 : for (int SurfNode = 1; SurfNode <= this->NumOfPerpendNodes; ++SurfNode) {
1823 0 : if ((SurfNode == 1) || (SurfNode == this->NumOfPerpendNodes)) {
1824 0 : SurfNodeFac = 0.5;
1825 : } else {
1826 0 : SurfNodeFac = 1.0;
1827 : }
1828 0 : this->s(j, 1, this->rcmax) += SurfNodeFac * this->CMat(1) * Rnew(is2, SurfNode) * this->Gamma1(j, is2);
1829 0 : this->s(j, 2, this->rcmax) += SurfNodeFac * this->CMat(2) * Rnew(is2, this->rcmax + 1 - SurfNode) * this->Gamma1(j, is2);
1830 : }
1831 0 : if (this->NodeSource > 0) {
1832 0 : this->s(j, 3, this->rcmax) += Rnew(is2, this->NodeSource) * this->Gamma1(j, is2);
1833 : }
1834 0 : if (this->NodeUserTemp > 0) {
1835 0 : this->s(j, 4, this->rcmax) += Rnew(is2, this->NodeUserTemp) * this->Gamma1(j, is2);
1836 : }
1837 : }
1838 : }
1839 0 : this->s(1, 1, this->rcmax) += this->e(this->rcmax) * this->DMat(1) * double(this->NumOfPerpendNodes - 1);
1840 0 : this->s(2, 2, this->rcmax) += this->e(this->rcmax) * this->DMat(2) * double(this->NumOfPerpendNodes - 1);
1841 : }
1842 :
1843 : // Check for and enforce symmetry in the cross term (Y)
1844 :
1845 54 : if (std::abs(s(2, 1, this->rcmax)) != std::abs(s(1, 2, this->rcmax))) {
1846 54 : avg = (std::abs(s(2, 1, this->rcmax)) + std::abs(s(1, 2, this->rcmax))) * 0.5;
1847 54 : this->s(2, 1, this->rcmax) *= avg / std::abs(s(2, 1, this->rcmax));
1848 54 : this->s(1, 2, this->rcmax) *= avg / std::abs(s(1, 2, this->rcmax));
1849 : }
1850 :
1851 : } // ... end of IF block for calculation of last e and S.
1852 :
1853 146 : PhiR0.deallocate();
1854 146 : Rnew.deallocate();
1855 146 : Rold.deallocate();
1856 146 : }
1857 :
1858 74 : void ConstructionProps::reportTransferFunction(EnergyPlusData &state, int const cCounter)
1859 : {
1860 :
1861 : static constexpr std::string_view Format_700{" Construction CTF,{},{:4},{:4},{:4},{:8.3F},{:15.4N},{:8.3F},{:8.3F},{:8.3F},{:8.3F},{}\n"};
1862 74 : print(state.files.eio,
1863 : Format_700,
1864 74 : this->Name,
1865 : cCounter,
1866 74 : this->TotLayers,
1867 74 : this->NumCTFTerms,
1868 74 : this->CTFTimeStep,
1869 74 : this->UValue,
1870 74 : this->OutsideAbsorpThermal,
1871 74 : this->InsideAbsorpThermal,
1872 74 : this->OutsideAbsorpSolar,
1873 74 : this->InsideAbsorpSolar,
1874 74 : Material::surfaceRoughnessNames[(int)this->OutsideRoughness]);
1875 :
1876 164 : for (int I = 1; I <= this->TotLayers; ++I) {
1877 90 : int Layer = this->LayerPoint(I);
1878 90 : auto const *thisMaterial = state.dataMaterial->materials(Layer);
1879 90 : switch (thisMaterial->group) {
1880 1 : case Material::Group::AirGap: {
1881 : static constexpr std::string_view Format_702(" Material:Air,{},{:12.4N}\n");
1882 1 : print(state.files.eio, Format_702, thisMaterial->Name, thisMaterial->Resistance);
1883 1 : } break;
1884 89 : default: {
1885 : static constexpr std::string_view Format_701(" Material CTF Summary,{},{:8.4F},{:14.3F},{:11.3F},{:13.3F},{:12.4N}\n");
1886 89 : Material::MaterialBase const *mp = thisMaterial;
1887 89 : print(state.files.eio, Format_701, mp->Name, mp->Thickness, mp->Conductivity, mp->Density, mp->SpecHeat, mp->Resistance);
1888 89 : } break;
1889 : }
1890 : }
1891 :
1892 380 : for (int I = this->NumCTFTerms; I >= 0; --I) {
1893 306 : if (I != 0) {
1894 : static constexpr std::string_view Format_703(" CTF,{:4},{:20.8N},{:20.8N},{:20.8N},{:20.8N}\n");
1895 232 : print(state.files.eio, Format_703, I, this->CTFOutside[I], this->CTFCross[I], this->CTFInside[I], this->CTFFlux[I]);
1896 : } else {
1897 : static constexpr std::string_view Format_704(" CTF,{:4},{:20.8N},{:20.8N},{:20.8N}\n");
1898 74 : print(state.files.eio, Format_704, I, this->CTFOutside[I], this->CTFCross[I], this->CTFInside[I]);
1899 : }
1900 : }
1901 :
1902 74 : if (this->SourceSinkPresent) {
1903 : // QTFs...
1904 0 : for (int I = this->NumCTFTerms; I >= 0; --I) {
1905 : static constexpr std::string_view Format_705(" QTF,{:4},{:20.8N},{:20.8N}\n");
1906 0 : print(state.files.eio, Format_705, I, this->CTFSourceOut[I], this->CTFSourceIn[I]);
1907 : }
1908 : // QTFs for source/sink location temperature calculation...
1909 0 : for (int I = this->NumCTFTerms; I >= 0; --I) {
1910 : static constexpr std::string_view Format_706(" Source/Sink Loc Internal Temp QTF,{:4},{:20.8N},{:20.8N},{:20.8N}\n");
1911 0 : print(state.files.eio, Format_706, I, this->CTFTSourceOut[I], this->CTFTSourceIn[I], this->CTFTSourceQ[I]);
1912 : }
1913 0 : if (this->TempAfterLayer != 0) {
1914 : // QTFs for user specified interior temperature calculation...
1915 0 : for (int I = this->NumCTFTerms; I >= 0; --I) {
1916 : static constexpr std::string_view Format_707(" User Loc Internal Temp QTF,{:4},{:20.8N},{:20.8N},{:20.8N}\n");
1917 0 : print(state.files.eio, Format_707, I, this->CTFTUserOut[I], this->CTFTUserIn[I], this->CTFTUserSource[I]);
1918 : }
1919 : }
1920 : }
1921 74 : }
1922 :
1923 249 : void ConstructionProps::reportLayers(EnergyPlusData &state)
1924 : {
1925 : // Report the layers for each opaque construction in predefined tabular report
1926 : // J. Glazer March 2024
1927 249 : if (state.dataOutRptPredefined->pdchOpqConsLayCol.size() > 0) {
1928 612 : for (int i = 1; i <= this->TotLayers; ++i) {
1929 382 : int layerIndex = this->LayerPoint(i);
1930 382 : auto const *mat = state.dataMaterial->materials(layerIndex);
1931 382 : OutputReportPredefined::PreDefTableEntry(state, state.dataOutRptPredefined->pdchOpqConsLayCol[i - 1], this->Name, mat->Name);
1932 : }
1933 : }
1934 249 : }
1935 :
1936 4 : bool ConstructionProps::isGlazingConstruction(EnergyPlusData &state) const
1937 : {
1938 : // SUBROUTINE INFORMATION:
1939 : // AUTHOR Simon Vidanovic
1940 : // DATE WRITTEN September 2016
1941 :
1942 : // PURPOSE OF THIS SUBROUTINE:
1943 : // Commonly used routine in several places in EnergyPlus which examines if current
1944 : // construction is glazing construction
1945 4 : const Material::Group group = state.dataMaterial->materials(LayerPoint(1))->group;
1946 0 : return (group == Material::Group::Glass) || (group == Material::Group::Shade) || (group == Material::Group::Screen) ||
1947 4 : (group == Material::Group::Blind) || (group == Material::Group::GlassSimple);
1948 : }
1949 :
1950 14 : Real64 ConstructionProps::setThicknessPerpendicular(EnergyPlusData &state, Real64 userValue)
1951 : {
1952 : // Limits set here are similar to limits used in HydronicSystemBaseData::sizeRadiantSystemTubeLength
1953 : // which is for autosizing tube length. The limits are not as strictly enforced here because they are
1954 : // not being used for autosizing so more latitude with the values is desired. Warnings will still alert
1955 : // the user if values seem outside of the ordinary range anticipated.
1956 14 : Real64 returnValue = userValue / 2.0; // Divide by two because internally the half distance will be used to calculate CTFs
1957 14 : if (returnValue <= 0.001) { // lowest reasonable value for "half" the tube spacing in meters
1958 4 : ShowWarningError(state, "ConstructionProperty:InternalHeatSource has a tube spacing that is less than 2 mm. This is not allowed.");
1959 4 : ShowContinueError(
1960 : state,
1961 4 : format("Construction={} has this problem. The tube spacing has been reset to 0.15m (~6 inches) for this construction.", this->Name));
1962 4 : ShowContinueError(state, "As per the Input Output Reference, tube spacing is only used for 2-D solutions and autosizing.");
1963 2 : returnValue = 0.075; // default "half" tube spacing in meters (roughly equivalent to 15cm or 6 inches of tube spacing)
1964 12 : } else if (returnValue < 0.005) { // below this value for "half" the tube spacing in meters throw a warning
1965 2 : ShowWarningError(state, "ConstructionProperty:InternalHeatSource has a tube spacing that is less than 1 cm (0.4 inch).");
1966 1 : ShowContinueError(state, format("Construction={} has this concern. Please check this construction to make sure it is correct.", this->Name));
1967 3 : ShowContinueError(state, "As per the Input Output Reference, tube spacing is only used for 2-D solutions and autosizing.");
1968 11 : } else if (returnValue > 0.5) { // above this value for "half" the tube spacing in meters throw a warning
1969 2 : ShowWarningError(state, "ConstructionProperty:InternalHeatSource has a tube spacing that is greater than 1 meter (39.4 inches).");
1970 1 : ShowContinueError(state, format("Construction={} has this concern. Please check this construction to make sure it is correct.", this->Name));
1971 3 : ShowContinueError(state, "As per the Input Output Reference, tube spacing is only used for 2-D solutions and autosizing.");
1972 : }
1973 14 : return returnValue;
1974 : }
1975 :
1976 12 : Real64 ConstructionProps::setUserTemperatureLocationPerpendicular(EnergyPlusData &state, Real64 userValue)
1977 : {
1978 12 : if (userValue < 0.0) {
1979 2 : ShowWarningError(state, "ConstructionProperty:InternalHeatSource has a perpendicular temperature location parameter that is less than zero.");
1980 1 : ShowContinueError(state, format("Construction={} has this error. The parameter has been reset to 0.", this->Name));
1981 1 : return 0.0;
1982 11 : } else if (userValue > 1.0) {
1983 2 : ShowWarningError(state,
1984 : "ConstructionProperty:InternalHeatSource has a perpendicular temperature location parameter that is greater than one.");
1985 1 : ShowContinueError(state, format("Construction={} has this error. The parameter has been reset to 1.", this->Name));
1986 1 : return 1.0;
1987 : } else { // Valid value between 0 and 1
1988 10 : return userValue;
1989 : }
1990 : }
1991 :
1992 150 : void ConstructionProps::setNodeSourceAndUserTemp(Array1D_int &Nodes)
1993 : {
1994 150 : this->NodeSource = 0;
1995 150 : this->NodeUserTemp = 0;
1996 150 : if (!this->SourceSinkPresent) return;
1997 :
1998 10 : for (int Layer = 1; Layer <= this->SourceAfterLayer; ++Layer) {
1999 7 : this->NodeSource += Nodes(Layer);
2000 : }
2001 :
2002 3 : if ((this->NodeSource > 0) && (this->SolutionDimensions > 1)) this->NodeSource = (this->NodeSource - 1) * this->NumOfPerpendNodes + 1;
2003 :
2004 13 : for (int Layer = 1; Layer <= this->TempAfterLayer; ++Layer) {
2005 10 : this->NodeUserTemp += Nodes(Layer);
2006 : }
2007 :
2008 3 : if ((this->NodeUserTemp > 0) && (this->SolutionDimensions > 1))
2009 2 : this->NodeUserTemp = (this->NodeUserTemp - 1) * this->NumOfPerpendNodes +
2010 2 : round(this->userTemperatureLocationPerpendicular * (this->NumOfPerpendNodes - 1)) + 1;
2011 : }
2012 :
2013 798 : void ConstructionProps::setArraysBasedOnMaxSolidWinLayers(EnergyPlusData &state)
2014 : {
2015 798 : this->AbsDiff.allocate(state.dataHeatBal->MaxSolidWinLayers);
2016 798 : this->AbsDiffBack.allocate(state.dataHeatBal->MaxSolidWinLayers);
2017 798 : this->layerSlatBlindDfAbs.allocate(state.dataHeatBal->MaxSolidWinLayers);
2018 798 : this->AbsBeamCoef.allocate(state.dataHeatBal->MaxSolidWinLayers);
2019 798 : this->AbsBeamBackCoef.allocate(state.dataHeatBal->MaxSolidWinLayers);
2020 798 : this->tBareSolCoef.allocate(state.dataHeatBal->MaxSolidWinLayers);
2021 798 : this->tBareVisCoef.allocate(state.dataHeatBal->MaxSolidWinLayers);
2022 798 : this->rfBareSolCoef.allocate(state.dataHeatBal->MaxSolidWinLayers);
2023 798 : this->rfBareVisCoef.allocate(state.dataHeatBal->MaxSolidWinLayers);
2024 798 : this->rbBareSolCoef.allocate(state.dataHeatBal->MaxSolidWinLayers);
2025 798 : this->rbBareVisCoef.allocate(state.dataHeatBal->MaxSolidWinLayers);
2026 798 : this->afBareSolCoef.allocate(state.dataHeatBal->MaxSolidWinLayers);
2027 798 : this->abBareSolCoef.allocate(state.dataHeatBal->MaxSolidWinLayers);
2028 :
2029 3603 : for (int Layer = 1; Layer <= state.dataHeatBal->MaxSolidWinLayers; ++Layer) {
2030 2805 : this->AbsDiff(Layer) = 0.0;
2031 2805 : this->AbsDiffBack(Layer) = 0.0;
2032 : }
2033 3603 : for (int Layer = 1; Layer <= state.dataHeatBal->MaxSolidWinLayers; ++Layer) {
2034 2805 : auto &slatBlindDfAbs = this->layerSlatBlindDfAbs(Layer);
2035 510510 : for (int iSlatAng = 0; iSlatAng < Material::MaxSlatAngs; ++iSlatAng) {
2036 507705 : auto &blindDfAbs = slatBlindDfAbs[iSlatAng];
2037 507705 : blindDfAbs.Sol.Ft.Df.Abs = 0.0;
2038 507705 : blindDfAbs.Sol.Ft.Df.AbsGnd = 0.0;
2039 507705 : blindDfAbs.Sol.Ft.Df.AbsGnd = 0.0;
2040 507705 : blindDfAbs.Sol.Bk.Df.Abs = 0.0;
2041 : }
2042 : }
2043 3603 : for (int Layer = 1; Layer <= state.dataHeatBal->MaxSolidWinLayers; ++Layer) {
2044 2805 : std::fill(this->AbsBeamCoef(Layer).begin(), this->AbsBeamCoef(Layer).end(), 0.0);
2045 2805 : std::fill(this->AbsBeamBackCoef(Layer).begin(), this->AbsBeamBackCoef(Layer).end(), 0.0);
2046 2805 : std::fill(this->tBareSolCoef(Layer).begin(), this->tBareSolCoef(Layer).end(), 0.0);
2047 2805 : std::fill(this->tBareVisCoef(Layer).begin(), this->tBareVisCoef(Layer).end(), 0.0);
2048 2805 : std::fill(this->rfBareSolCoef(Layer).begin(), this->rfBareSolCoef(Layer).end(), 0.0);
2049 2805 : std::fill(this->rfBareVisCoef(Layer).begin(), this->rfBareVisCoef(Layer).end(), 0.0);
2050 2805 : std::fill(this->rbBareSolCoef(Layer).begin(), this->rbBareSolCoef(Layer).end(), 0.0);
2051 2805 : std::fill(this->rbBareVisCoef(Layer).begin(), this->rbBareVisCoef(Layer).end(), 0.0);
2052 2805 : std::fill(this->afBareSolCoef(Layer).begin(), this->afBareSolCoef(Layer).end(), 0.0);
2053 2805 : std::fill(this->abBareSolCoef(Layer).begin(), this->abBareSolCoef(Layer).end(), 0.0);
2054 : }
2055 798 : }
2056 :
2057 : } // namespace EnergyPlus::Construction
|
