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