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