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