Line data Source code
1 : // EnergyPlus, Copyright (c) 1996-2023, The Board of Trustees of the University of Illinois,
2 : // The Regents of the University of California, through Lawrence Berkeley National Laboratory
3 : // (subject to receipt of any required approvals from the U.S. Dept. of Energy), Oak Ridge
4 : // National Laboratory, managed by UT-Battelle, Alliance for Sustainable Energy, LLC, and other
5 : // contributors. All rights reserved.
6 : //
7 : // NOTICE: This Software was developed under funding from the U.S. Department of Energy and the
8 : // U.S. Government consequently retains certain rights. As such, the U.S. Government has been
9 : // granted for itself and others acting on its behalf a paid-up, nonexclusive, irrevocable,
10 : // worldwide license in the Software to reproduce, distribute copies to the public, prepare
11 : // derivative works, and perform publicly and display publicly, and to permit others to do so.
12 : //
13 : // Redistribution and use in source and binary forms, with or without modification, are permitted
14 : // provided that the following conditions are met:
15 : //
16 : // (1) Redistributions of source code must retain the above copyright notice, this list of
17 : // conditions and the following disclaimer.
18 : //
19 : // (2) Redistributions in binary form must reproduce the above copyright notice, this list of
20 : // conditions and the following disclaimer in the documentation and/or other materials
21 : // provided with the distribution.
22 : //
23 : // (3) Neither the name of the University of California, Lawrence Berkeley National Laboratory,
24 : // the University of Illinois, U.S. Dept. of Energy nor the names of its contributors may be
25 : // used to endorse or promote products derived from this software without specific prior
26 : // written permission.
27 : //
28 : // (4) Use of EnergyPlus(TM) Name. If Licensee (i) distributes the software in stand-alone form
29 : // without changes from the version obtained under this License, or (ii) Licensee makes a
30 : // reference solely to the software portion of its product, Licensee must refer to the
31 : // software as "EnergyPlus version X" software, where "X" is the version number Licensee
32 : // obtained under this License and may not use a different name for the software. Except as
33 : // specifically required in this Section (4), Licensee shall not use in a company name, a
34 : // product name, in advertising, publicity, or other promotional activities any name, trade
35 : // name, trademark, logo, or other designation of "EnergyPlus", "E+", "e+" or confusingly
36 : // similar designation, without the U.S. Department of Energy's prior written consent.
37 : //
38 : // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
39 : // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY
40 : // AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
41 : // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
42 : // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
43 : // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
44 : // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
45 : // OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
46 : // POSSIBILITY OF SUCH DAMAGE.
47 :
48 : // C++ Headers
49 : #include <cmath>
50 :
51 : // ObjexxFCL Headers
52 : #include <ObjexxFCL/Fmath.hh>
53 :
54 : // EnergyPlus Headers
55 : #include <EnergyPlus/Construction.hh>
56 : #include <EnergyPlus/ConvectionCoefficients.hh>
57 : #include <EnergyPlus/Data/EnergyPlusData.hh>
58 : #include <EnergyPlus/DataEnvironment.hh>
59 : #include <EnergyPlus/DataHVACGlobals.hh>
60 : #include <EnergyPlus/DataHeatBalSurface.hh>
61 : #include <EnergyPlus/DataHeatBalance.hh>
62 : #include <EnergyPlus/DataSurfaces.hh>
63 : #include <EnergyPlus/DataWater.hh>
64 : #include <EnergyPlus/EcoRoofManager.hh>
65 : #include <EnergyPlus/Material.hh>
66 : #include <EnergyPlus/OutputProcessor.hh>
67 : #include <EnergyPlus/SolarShading.hh>
68 : #include <EnergyPlus/UtilityRoutines.hh>
69 : #include <EnergyPlus/WeatherManager.hh>
70 : #include <EnergyPlus/ZoneTempPredictorCorrector.hh>
71 :
72 : namespace EnergyPlus {
73 :
74 : namespace EcoRoofManager {
75 : // Module containing the heat balance simulation routines
76 : // calculation (initialization) routines
77 :
78 : // MODULE INFORMATION:
79 : // AUTHOR David Sailor and Toan Pham, Portland State University
80 : // DATE WRITTEN Jan 2007
81 : // MODIFIED Oct 2010
82 : // RE-ENGINEERED na
83 :
84 : // PURPOSE OF THIS MODULE:
85 : // Module for implementing an ecoroof (aka Green Roof)
86 :
87 : // METHODOLOGY EMPLOYED:
88 : // Vikram Madhusudan's Portland State Univ. MS Thesis (Dec 2005) based on FASST model
89 : // of Frankenstein and Koenig (2004) - DRDC/CRREL Technical Report TR-04-25.
90 : // Precipitation schedules and irrigation schedules can be used to define hourly moisture
91 : // inputs (m). Moisture transport updated Oct 2010.
92 : // REFERENCES:
93 : // OTHER NOTES:
94 :
95 : // USE STATEMENTS:
96 : // Use statements for data only modules
97 : // Using/Aliasing
98 : using namespace DataSurfaces;
99 : using namespace DataHeatBalance;
100 :
101 : // Functions
102 :
103 443820 : void CalcEcoRoof(EnergyPlusData &state,
104 : int const SurfNum, // Indicator of Surface Number for the current surface
105 : int const ZoneNum, // Indicator for zone number where the current surface
106 : int const ConstrNum, // Indicator for construction index for the current surface
107 : Real64 &TempExt // Exterior temperature boundary condition
108 : )
109 : {
110 : // SUBROUTINE INFORMATION
111 : // AUTHOR David Sailor and Toan Pham
112 : // DATE WRITTEN January 2007
113 : // MODIFIED David Sailor - to fix initialization between DD runs and during warm-up
114 : // RE-ENGINEERED na
115 :
116 : // PURPOSE OF THIS MODULE:
117 :
118 : // To calculate the heat balance for surfaces with eco roof specified as outside surface
119 : // Note that only ONE ecoroof construction can be employed at present time. If multiple
120 : // surfaces have ecoroof as the outside layer the energy balance is only calculated for
121 : // the first such surface.
122 :
123 : // METHODOLOGY EMPLOYED:
124 : // Vikram Madhusudan's Portland State Univ. MS Thesis (Dec 2005) based on FASST model
125 : // of Frankenstein and Koenig (2004) - DRDC/CRREL Technical Report TR-04-25.
126 : // Some data used herein are from: European Centre for Medium-Range Weather Forecasts (ECMWF)
127 : // IFS Documentation, CY25R1 (April 2002), www.ecmwf.int/research/ifsdocs/CY24r1/Physics/
128 : // Physics-08-03.html.
129 : // The Atmospheric Boundary Layer - by J.R. Garratt (Cambridge Atmos. & Space Science Series), 316pp.
130 : // Using/Aliasing
131 : using namespace DataEnvironment;
132 : using namespace DataHeatBalance;
133 : using namespace DataHeatBalSurface;
134 : using namespace DataSurfaces;
135 : using ConvectionCoefficients::InitExteriorConvectionCoeff;
136 : using ConvectionCoefficients::SetExtConvectionCoeff;
137 : using ConvectionCoefficients::SetIntConvectionCoeff;
138 :
139 : // SUBROUTINE PARAMETER DEFINITIONS:
140 443820 : Real64 constexpr Kv(0.4); // Von Karmen's constant (source FASST)
141 443820 : Real64 constexpr rch(0.63); // Turbulent Schimdt Number
142 443820 : Real64 constexpr rche(0.71); // Turbulent Prandtl Number
143 443820 : Real64 constexpr Rair(0.286e3); // Gas Constant of air J/Kg K
144 443820 : Real64 constexpr g1(9.81); // Gravity. In m/sec^2.
145 443820 : Real64 constexpr Sigma(5.6697e-08); // Stefan-Boltzmann constant W/m^2K^4
146 443820 : Real64 constexpr Cpa(1005.6); // Specific heat of Water Vapor. (J/Kg.K)
147 :
148 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
149 : int EcoLoop; // an integer loop variable for the simultaneous solution iteration
150 :
151 : Real64 AbsThermSurf; // Thermal absoptance of the exterior surface
152 : DataSurfaces::SurfaceRoughness RoughSurf; // Roughness index of the exterior (ecoroof) surface.
153 : Real64 HMovInsul; // "Convection" coefficient of movable insulation
154 : Real64 Tgk; // Ground temperature in Kelvin
155 : Real64 Ta; // current air temperature
156 : Real64 Ws; // Wind Speed (m/s)
157 : Real64 Waf; // Windspeed within canopy (m/s)
158 : Real64 Latm; // Long Wave Radiation (W/m^2)
159 : Real64 qaf; // mixing ratio of air near canopy
160 : Real64 qg; // mixing ratio of air at surface.
161 : Real64 RS; // shortwave radiation
162 : Real64 EpsilonOne;
163 : Real64 eair;
164 : Real64 Rhoa;
165 : Real64 Tak;
166 : Real64 qa; // mixing ratio of air
167 : Real64 Tafk;
168 : Real64 Taf;
169 : Real64 Rhof;
170 : Real64 Rhoaf; // Average air density
171 : Real64 sigmaf;
172 : Real64 Zd; // zero displacement height (m)
173 : Real64 Zo; // foliage roughness length (m)
174 : Real64 Cfhn; // transfer coefficient at near-neutral conditions
175 : Real64 Cf; // bulk Transfer coefficient, equation 10 page 6 (FASST).
176 : Real64 ra; // Aerodynamic Resistance
177 : Real64 f2inv; // intermediate calculation variable
178 : Real64 f1; // intermediate calculation variable
179 : Real64 f2; // intermediate calculation variable
180 : Real64 r_s; // Minimum Stomatal Resistance, specific to each plant
181 : Real64 Mg; // Surface soil moisture content m^3/m^3 (Moisture / MoistureMax)
182 : Real64 dOne; // intermediate calculation variable
183 : Real64 esf; // the saturation vapor pressure (Pa)
184 : Real64 qsf; // Saturation specific humidity at leaf temperature (qfsat)
185 : Real64 Lef; // Latent heat of vaporation at leaf surface temperature (J/kg)
186 : Real64 Desf; // Derivative of Saturation vapor pressure
187 : Real64 dqf; // Derivative of saturation specific humidity
188 : Real64 dqg; // this is given by Clausius-Clapeyron equation
189 : Real64 esg; // Saturation vapor pressure (Pa)
190 : Real64 qsg; // Saturation specific humidity(mixing ratio?) at ground surface temperature
191 : Real64 Leg; // Latent heat vaporization at the ground temperature (J/kg)
192 : Real64 Desg; // derivative of esg Saturation vapor pressure(?)
193 : Real64 F1temp; // intermediate variable in computing flux from the soil layer
194 : Real64 P1; // intermediate variable in the equation for Tf
195 : Real64 P2; // intermediate variable in the equation for Tf and Tg
196 : Real64 P3; // intermediate variable in the equation for Tf and Tg
197 : Real64 Rhog; // Density of air at the soil surface temperature
198 : Real64 Rhoag; // Average density of air with respect to ground surface and air temperature
199 : Real64 Rib; // Richardson Number
200 : Real64 Chng; // bulk transfer coefficient near ground
201 : Real64 Ce; // bulk transfer coefficient (this is in fact Ceg in equation 28 main report)
202 : Real64 Gammah; // latent heat exchange stability correction factor
203 : Real64 Chg; // in fact it is the same as Ce (=Ceg) is transfer coefficient (but wot?)
204 : Real64 T3G; // intermediate variable in the equation for Tg
205 : Real64 T2G; // intermediate variable in the equation for Tg
206 : Real64 LeafTK; // the current leaf's temperature (Kelvin)
207 : Real64 SoilTK; // the current soil's temperature (Kelvin)
208 : Real64 Chne; // is similar to near ground bulk transfer coefficient for latent heat flux (at neutral condition)
209 : Real64 Tif; // previous leaf temperature
210 : Real64 rn; // rn is the combined effect of both stomatal and aerodynamic resistances
211 : // in fact this is called r'' in the main report
212 : Real64 T1G; // intermediate variable in the equation for Tg
213 : Real64 Qsoilpart1; // intermediate variable for evaluating Qsoil (part without the unknown)
214 : Real64 Qsoilpart2; // intermediate variable for evaluating Qsoil (part coeff of the ground temperature)
215 :
216 : // INTEGER,EXTERNAL :: GetNewUnitNumber ! external function to return a new (unique) unit for ecoroof writing
217 443820 : int unit(0); // not actually used in the function it is passed into
218 :
219 443820 : Ws = DataEnvironment::WindSpeedAt(state, state.dataSurface->Surface(SurfNum).Centroid.z); // use windspeed at Z of roof
220 443820 : if (Ws < 2.0) { // Later we need to adjust for building roof height...
221 443820 : Ws = 2.0; // Set minimum possible wind speed outside vegetation to 2.0 m/s
222 : // consistent with FASST TR-04-25 p. x (W' = 2.0)
223 : }
224 :
225 443820 : RoughSurf = state.dataMaterial->Material(state.dataConstruction->Construct(ConstrNum).LayerPoint(1)).Roughness;
226 443820 : AbsThermSurf = state.dataMaterial->Material(state.dataConstruction->Construct(ConstrNum).LayerPoint(1)).AbsorpThermal;
227 443820 : HMovInsul = 0.0;
228 :
229 443820 : if (state.dataSurface->Surface(SurfNum).ExtWind) {
230 443820 : InitExteriorConvectionCoeff(state,
231 : SurfNum,
232 : HMovInsul,
233 : RoughSurf,
234 : AbsThermSurf,
235 443820 : state.dataHeatBalSurf->SurfOutsideTempHist(1)(SurfNum),
236 443820 : state.dataHeatBalSurf->SurfHcExt(SurfNum),
237 443820 : state.dataHeatBalSurf->SurfHSkyExt(SurfNum),
238 443820 : state.dataHeatBalSurf->SurfHGrdExt(SurfNum),
239 443820 : state.dataHeatBalSurf->SurfHAirExt(SurfNum));
240 : }
241 :
242 887640 : Latm = 1.0 * Sigma * 1.0 * state.dataSurface->Surface(SurfNum).ViewFactorGround * pow_4(state.dataEnvrn->GroundTempKelvin) +
243 443820 : 1.0 * Sigma * 1.0 * state.dataSurface->Surface(SurfNum).ViewFactorSky * pow_4(state.dataEnvrn->SkyTempKelvin);
244 :
245 443820 : if (state.dataEcoRoofMgr->EcoRoofbeginFlag)
246 4 : initEcoRoofFirstTime(state, SurfNum, ConstrNum); // Initialization statements for first entry into ecoroof routines
247 :
248 443820 : initEcoRoof(state, SurfNum, ConstrNum);
249 :
250 : // If current surface is = FirstEcoSurf then for this time step we need to update the soil moisture
251 443820 : if (SurfNum == state.dataEcoRoofMgr->FirstEcoSurf) {
252 976404 : UpdateSoilProps(state,
253 88764 : state.dataEcoRoofMgr->Moisture,
254 88764 : state.dataEcoRoofMgr->MeanRootMoisture,
255 88764 : state.dataEcoRoofMgr->MoistureMax,
256 88764 : state.dataEcoRoofMgr->MoistureResidual,
257 88764 : state.dataEcoRoofMgr->SoilThickness,
258 88764 : state.dataEcoRoofMgr->Vfluxf,
259 88764 : state.dataEcoRoofMgr->Vfluxg,
260 : ConstrNum,
261 88764 : state.dataEcoRoofMgr->Alphag,
262 : unit,
263 88764 : state.dataEcoRoofMgr->Tg,
264 88764 : state.dataEcoRoofMgr->Tf,
265 88764 : state.dataEcoRoofMgr->Qsoil);
266 :
267 88764 : Ta = OutDryBulbTempAt(state,
268 88764 : state.dataSurface->Surface(SurfNum).Centroid.z); // temperature outdoor - Surface is dry, use normal correlation
269 88764 : state.dataEcoRoofMgr->Tg = state.dataEcoRoofMgr->Tgold;
270 88764 : state.dataEcoRoofMgr->Tf = state.dataEcoRoofMgr->Tfold;
271 :
272 88764 : if (state.dataConstruction->Construct(ConstrNum).CTFCross(0) > 0.01) {
273 0 : state.dataEcoRoofMgr->QuickConductionSurf = true;
274 0 : F1temp = state.dataConstruction->Construct(ConstrNum).CTFCross(0) /
275 0 : (state.dataConstruction->Construct(ConstrNum).CTFInside(0) + state.dataHeatBalSurf->SurfHConvInt(SurfNum));
276 0 : Qsoilpart1 =
277 0 : -state.dataHeatBalSurf->SurfCTFConstOutPart(SurfNum) +
278 0 : F1temp * (state.dataHeatBalSurf->SurfCTFConstInPart(SurfNum) + state.dataHeatBalSurf->SurfOpaqQRadSWInAbs(SurfNum) +
279 0 : state.dataHeatBal->SurfQdotRadIntGainsInPerArea(SurfNum) +
280 0 : state.dataConstruction->Construct(ConstrNum).CTFSourceIn(0) * state.dataHeatBalSurf->SurfQsrcHist(SurfNum, 1) +
281 0 : state.dataHeatBalSurf->SurfHConvInt(SurfNum) * state.dataZoneTempPredictorCorrector->zoneHeatBalance(ZoneNum).MAT +
282 0 : state.dataHeatBalSurf->SurfQdotRadNetLWInPerArea(SurfNum));
283 : } else {
284 88764 : Qsoilpart1 = -state.dataHeatBalSurf->SurfCTFConstOutPart(SurfNum) +
285 88764 : state.dataConstruction->Construct(ConstrNum).CTFCross(0) * state.dataHeatBalSurf->SurfTempIn(SurfNum);
286 88764 : F1temp = 0.0;
287 : }
288 :
289 88764 : Qsoilpart2 =
290 88764 : state.dataConstruction->Construct(ConstrNum).CTFOutside(0) - F1temp * state.dataConstruction->Construct(ConstrNum).CTFCross(0);
291 :
292 88764 : state.dataEcoRoofMgr->Pa = state.dataEnvrn->StdBaroPress; // standard atmospheric pressure (apparently in Pascals)
293 88764 : Tgk = state.dataEcoRoofMgr->Tg + DataGlobalConstants::KelvinConv;
294 88764 : Tak = Ta + DataGlobalConstants::KelvinConv;
295 :
296 88764 : sigmaf = 0.9 - 0.7 * std::exp(-0.75 * state.dataEcoRoofMgr->LAI); // Fractional veg cover based on (2) from FASST TR-04-25
297 : // Formula for grasses modified to incorporate limits from
298 : // Table 1 for sigmaf_max and min (0.20 to 0.9)
299 :
300 177528 : EpsilonOne = state.dataEcoRoofMgr->epsilonf + state.dataEcoRoofMgr->epsilong -
301 88764 : state.dataEcoRoofMgr->epsilong * state.dataEcoRoofMgr->epsilonf; // Checked (eqn. 6 in FASST Veg Models)
302 88764 : state.dataEcoRoofMgr->RH = state.dataEnvrn->OutRelHum; // Get humidity in % from the DataEnvironment.cc
303 88764 : eair = (state.dataEcoRoofMgr->RH / 100.0) * 611.2 * std::exp(17.67 * Ta / (Tak - 29.65));
304 88764 : qa = (0.622 * eair) / (state.dataEcoRoofMgr->Pa - 1.000 * eair); // Mixing Ratio of air
305 88764 : Rhoa = state.dataEcoRoofMgr->Pa / (Rair * Tak); // Density of air. kg/m^3
306 88764 : Tif = state.dataEcoRoofMgr->Tf;
307 :
308 : // Air Temperature within the canopy is given as
309 : // (Deardorff (1987)). Kelvin. based of the previous temperatures
310 88764 : Tafk = (1.0 - sigmaf) * Tak + sigmaf * (0.3 * Tak + 0.6 * (Tif + DataGlobalConstants::KelvinConv) + 0.1 * Tgk);
311 :
312 88764 : Taf = Tafk - DataGlobalConstants::KelvinConv; // Air Temperature within canopy in Celcius (C).
313 88764 : Rhof = state.dataEcoRoofMgr->Pa / (Rair * Tafk); // Density of air at the leaf temperature
314 88764 : Rhoaf = (Rhoa + Rhof) / 2.0; // Average of air density
315 88764 : Zd = 0.701 * std::pow(state.dataEcoRoofMgr->Zf, 0.979); // Zero displacement height
316 88764 : Zo = 0.131 * std::pow(state.dataEcoRoofMgr->Zf, 0.997); // Foliage roughness length. (m) Source Ballick (1981)
317 88764 : if (Zo < 0.02) Zo = 0.02; // limit based on p.7 TR-04-25 and Table 2
318 :
319 : // transfer coefficient at near-neutral condition Cfhn
320 88764 : Cfhn = pow_2(Kv / std::log((state.dataEcoRoofMgr->Za - Zd) / Zo)); // Equation 12, page 7, FASST Model
321 88764 : Waf = 0.83 * std::sqrt(Cfhn) * sigmaf * Ws + (1.0 - sigmaf) * Ws; // Wind Speed inside foliage. Equation #6, FASST model
322 88764 : Cf = 0.01 * (1.0 + 0.3 / Waf); // The bulk Transfer coefficient, equation 10 page 6.
323 88764 : state.dataEcoRoofMgr->sheatf = state.dataEcoRoofMgr->e0 + 1.1 * state.dataEcoRoofMgr->LAI * Rhoaf * Cpa * Cf *
324 : Waf; // Intermediate calculation for Sensible Heat Transfer
325 88764 : state.dataEcoRoofMgr->sensiblef =
326 177528 : state.dataEcoRoofMgr->sheatf *
327 88764 : (Taf - state.dataEcoRoofMgr->Tf); // DJS Jan 2011 sensible flux TO foliage into air (Frankenstein 2004, eqn7)
328 : // sourced from Frankenstein et al (2004a). Added e0 windless correction factor.
329 : // an early version had (1.0-0.7)*sigmaf in calc of sensiblef... how did that get there!?! Fixed.
330 :
331 : // These parameters were taken from "The Atm Boundary Layer", By J.R. Garratt
332 : // NOTE the Garratt eqn. (A21) gives esf in units of hPA so we have multiplied
333 : // the constant 6.112 by a factor of 100.
334 88764 : esf = 611.2 * std::exp(17.67 * Tif / (Tif + DataGlobalConstants::KelvinConv - 29.65));
335 :
336 : // From Garratt - eqn. A21, p284. Note that Tif and Tif+DataGlobalConstants::KelvinConv() usage is correct.
337 : // Saturation specific humidity at leaf temperature again based on previous temperatures
338 :
339 88764 : qsf = 0.622 * esf / (state.dataEcoRoofMgr->Pa - 1.000 * esf); // "The Atm Boundary Layer", J.R Garrat for Saturation mixing ratio
340 : // Calculate stomatal resistance and atmospheric resistance Part
341 88764 : ra = 1.0 / (Cf * Waf); // Aerodynamic Resistance. Resistance that is caused
342 : // by the boundary layer on a leaf surface to transfer water vapor. It is measured in
343 : // s/m and depends on wind speed, leaf's surface roughness,
344 : // and stability of atsmophere.
345 :
346 88764 : CalculateEcoRoofSolar(state, RS, f1, SurfNum);
347 88764 : if (state.dataEcoRoofMgr->MoistureMax == state.dataEcoRoofMgr->MoistureResidual) {
348 0 : f2inv = 1.0e10;
349 : } else {
350 177528 : f2inv = (state.dataEcoRoofMgr->MeanRootMoisture - state.dataEcoRoofMgr->MoistureResidual) /
351 88764 : (state.dataEcoRoofMgr->MoistureMax - state.dataEcoRoofMgr->MoistureResidual); // Equation 19 p. 9 FASST model
352 : }
353 :
354 : // In FASST, Eq 20 is used to compute Moisture.
355 :
356 88764 : f2 = 1.0 / f2inv; // In dry areas f2 --> LARGE so that r_s --> LARGE
357 : // and both rn and Latent flux --> 0.0
358 88764 : state.dataEcoRoofMgr->f3 = 1.0 / (std::exp(-0.0 * (esf - eair))); // Note the 0.0 here is gd which has value of 0.0
359 : // for most plants and 0.03 for trees (see ECMWF
360 : // note given above.
361 177528 : r_s = state.dataEcoRoofMgr->StomatalResistanceMin * f1 * f2 * state.dataEcoRoofMgr->f3 /
362 88764 : state.dataEcoRoofMgr->LAI; // Stomatal Resistance (r_s)
363 88764 : rn = ra / (ra + r_s); // rn is foliage surface wetness ... NOT a resistance
364 :
365 : // This routine is to calculate ground moisture factor. This factor is from *****
366 88764 : Mg = state.dataEcoRoofMgr->Moisture / state.dataEcoRoofMgr->MoistureMax; // m^3/m^3.
367 88764 : dOne = 1.0 - sigmaf * (0.6 * (1.0 - rn) + 0.1 * (1.0 - Mg));
368 :
369 : // Latent heat of vaporation at leaf surface temperature. The source of this
370 : // equation is Henderson-Sellers (1984)
371 88764 : Lef = 1.91846e6 * pow_2((Tif + DataGlobalConstants::KelvinConv) / (Tif + DataGlobalConstants::KelvinConv - 33.91));
372 : // Check to see if ice is sublimating or frost is forming.
373 88764 : if (state.dataEcoRoofMgr->Tfold < 0.0) Lef = 2.838e6; // per FASST documentation p.15 after eqn. 37.
374 :
375 : // Derivative of Saturation vapor pressure, which is used in the calculation of
376 : // derivative of saturation specific humidity.
377 :
378 177528 : Desf = 611.2 * std::exp(17.67 * (state.dataEcoRoofMgr->Tf / (state.dataEcoRoofMgr->Tf + DataGlobalConstants::KelvinConv - 29.65))) *
379 177528 : (17.67 * state.dataEcoRoofMgr->Tf * (-1.0) * std::pow(state.dataEcoRoofMgr->Tf + DataGlobalConstants::KelvinConv - 29.65, -2) +
380 88764 : 17.67 / (DataGlobalConstants::KelvinConv - 29.65 + state.dataEcoRoofMgr->Tf));
381 88764 : dqf = ((0.622 * state.dataEcoRoofMgr->Pa) / pow_2(state.dataEcoRoofMgr->Pa - esf)) * Desf; // Derivative of saturation specific humidity
382 177528 : esg = 611.2 * std::exp(17.67 * (state.dataEcoRoofMgr->Tg /
383 88764 : ((state.dataEcoRoofMgr->Tg + DataGlobalConstants::KelvinConv) - 29.65))); // Pa saturation vapor pressure
384 : // From Garratt - eqn. A21, p284.
385 : // Note that Tg and Tg+DataGlobalConstants::KelvinConv() usage is correct.
386 88764 : qsg = 0.622 * esg / (state.dataEcoRoofMgr->Pa - esg); // Saturation mixing ratio at ground surface temperature.
387 :
388 : // Latent heat vaporization at the ground temperature
389 88764 : Leg = 1.91846e6 * pow_2(Tgk / (Tgk - 33.91));
390 : // Check to see if ice is sublimating or frost is forming.
391 88764 : if (state.dataEcoRoofMgr->Tgold < 0.0) Leg = 2.838e6; // per FASST documentation p.15 after eqn. 37.
392 :
393 177528 : Desg = 611.2 * std::exp(17.67 * (state.dataEcoRoofMgr->Tg / (state.dataEcoRoofMgr->Tg + DataGlobalConstants::KelvinConv - 29.65))) *
394 177528 : (17.67 * state.dataEcoRoofMgr->Tg * (-1.0) * std::pow(state.dataEcoRoofMgr->Tg + DataGlobalConstants::KelvinConv - 29.65, -2) +
395 88764 : 17.67 / (DataGlobalConstants::KelvinConv - 29.65 + state.dataEcoRoofMgr->Tg));
396 88764 : dqg = (0.622 * state.dataEcoRoofMgr->Pa / pow_2(state.dataEcoRoofMgr->Pa - esg)) * Desg;
397 :
398 : // Final Ground Atmosphere Energy Balance
399 : // Density of air at the soil surface temperature
400 88764 : Rhog = state.dataEcoRoofMgr->Pa / (Rair * Tgk);
401 :
402 : // Average density of air with respect to ground surface and air temperature
403 88764 : Rhoag = (Rhoa + Rhog) / 2.0;
404 88764 : Rib = 2.0 * g1 * state.dataEcoRoofMgr->Za * (Taf - state.dataEcoRoofMgr->Tg) / ((Tafk + Tgk) * pow_2(Waf)); // Richardson Number
405 :
406 : // Compute the stability factor Gammah
407 88764 : if (Rib < 0.0) {
408 88748 : Gammah = std::pow(1.0 - 16.0 * Rib, -0.5);
409 : } else {
410 16 : if (Rib >= 0.19) {
411 0 : Rib = 0.19;
412 : }
413 16 : Gammah = std::pow(1.0 - 5.0 * Rib, -0.5);
414 : }
415 :
416 88764 : if (RoughSurf == DataSurfaces::SurfaceRoughness::VerySmooth) {
417 0 : state.dataEcoRoofMgr->Zog = 0.0008;
418 88764 : } else if (RoughSurf == DataSurfaces::SurfaceRoughness::Smooth) {
419 0 : state.dataEcoRoofMgr->Zog = 0.0010;
420 88764 : } else if (RoughSurf == DataSurfaces::SurfaceRoughness::MediumSmooth) {
421 88764 : state.dataEcoRoofMgr->Zog = 0.0015;
422 0 : } else if (RoughSurf == DataSurfaces::SurfaceRoughness::MediumRough) {
423 0 : state.dataEcoRoofMgr->Zog = 0.0020;
424 0 : } else if (RoughSurf == DataSurfaces::SurfaceRoughness::Rough) {
425 0 : state.dataEcoRoofMgr->Zog = 0.0030;
426 : } else { // VeryRough
427 0 : state.dataEcoRoofMgr->Zog = 0.005;
428 : } // TODO: fix this after creating FindEnumeratedValueIndex()
429 :
430 88764 : Chng = pow_2(Kv / std::log(state.dataEcoRoofMgr->Za / state.dataEcoRoofMgr->Zog)) / rch; // bulk transfer coefficient near ground
431 88764 : Chg = Gammah * ((1.0 - sigmaf) * Chng + sigmaf * Cfhn);
432 88764 : state.dataEcoRoofMgr->sheatg = state.dataEcoRoofMgr->e0 + Rhoag * Cpa * Chg * Waf; // added the e0 windless correction
433 88764 : state.dataEcoRoofMgr->sensibleg =
434 177528 : state.dataEcoRoofMgr->sheatg *
435 88764 : (Taf - state.dataEcoRoofMgr->Tg); // sensible flux TO soil (W/m^2) DJS Jan 2011 (eqn. 32 in Frankenstein 2004)
436 :
437 88764 : Chne = pow_2(Kv / std::log(state.dataEcoRoofMgr->Za / state.dataEcoRoofMgr->Zog)) / rche;
438 88764 : Ce = Gammah * ((1.0 - sigmaf) * Chne + sigmaf * Cfhn); // this is in fact Ceg in eq (28)
439 :
440 : // we can approximate Gammae by Gammah (Supported by FASST Veg Models p. 15)
441 177528 : qaf = ((1.0 - sigmaf) * qa + sigmaf * (0.3 * qa + 0.6 * qsf * rn + 0.1 * qsg * Mg)) /
442 88764 : (1.0 - sigmaf * (0.6 * (1.0 - rn) + 0.1 * (1.0 - Mg)));
443 88764 : qg = Mg * qsg + (1.0 - Mg) * qaf; // eq main report (13)
444 : // According to FASST documentation this is correct.
445 : // The following also used Rhoaf and Rhoag respectively... shouldn't these be water densities??!!
446 88764 : state.dataEcoRoofMgr->Lf = Lef * state.dataEcoRoofMgr->LAI * Rhoaf * Cf * Waf * rn * (qaf - qsf); // This had Leg but should be Lef...
447 88764 : state.dataEcoRoofMgr->Lg = Ce * Leg * Waf * Rhoag * (qaf - qg) * Mg; // In the FASST documentation there is NO Mg. However, in looking
448 : // back at the Deardorff 1978 paper it appears that an alpha = Mg term is
449 : // used to distinguish from POTENTIAL and ACTUAL ground surface evaporation...
450 : // the Lf and Lg calculations are NOT used in this formulation
451 : // rather, the effects are included in terms of dqg and dqf !
452 : // These equations for Lf and Lg are based on Deardorff's paper, but there is
453 : // a sign difference !!! (qsf -qaf) and (qg - qaf) ?!
454 : // These may be useful, however for telling the UpdateSoilProps routine
455 : // how much evaporation has taken place...
456 88764 : state.dataEcoRoofMgr->Vfluxf = -1.0 * state.dataEcoRoofMgr->Lf / Lef / 990.0; // water evapotranspire rate [m/s]
457 88764 : state.dataEcoRoofMgr->Vfluxg = -1.0 * state.dataEcoRoofMgr->Lg / Leg / 990.0; // water evapotranspire rate [m/s]
458 88764 : if (state.dataEcoRoofMgr->Vfluxf < 0.0)
459 64662 : state.dataEcoRoofMgr->Vfluxf = 0.0; // According to FASST Veg. Models p. 11, eqn 26-27, if Qfsat > qaf the actual
460 88764 : if (state.dataEcoRoofMgr->Vfluxg < 0.0)
461 44354 : state.dataEcoRoofMgr->Vfluxg = 0.0; // evaporative fluxes should be set to zero (delta_c = 1 or 0).
462 :
463 : // P1, P2, P3 correspond to first, second and third terms of equation 37 in the main report.
464 :
465 : // Note: the FASST model has a term -gamma_p*(1.0-exp...) in first line for P1 (c1_f) where gamma_p is
466 : // a precipitation variable. So, if we assume no precip this term vanishes. We should
467 : // revisit this issue later.
468 : // Implement an iterative solution scheme to solve the simultaneous equations for Leaf and Soil temperature.
469 : // Prior experience suggests that no more than 3 iterations are likely needed
470 88764 : LeafTK = state.dataEcoRoofMgr->Tf + DataGlobalConstants::KelvinConv;
471 88764 : SoilTK = state.dataEcoRoofMgr->Tg + DataGlobalConstants::KelvinConv;
472 :
473 355056 : for (EcoLoop = 1; EcoLoop <= 3; ++EcoLoop) {
474 798876 : P1 = sigmaf * (RS * (1.0 - state.dataEcoRoofMgr->Alphaf) + state.dataEcoRoofMgr->epsilonf * Latm) -
475 532584 : 3.0 * sigmaf * state.dataEcoRoofMgr->epsilonf * state.dataEcoRoofMgr->epsilong * Sigma * pow_4(SoilTK) / EpsilonOne -
476 266292 : 3.0 *
477 532584 : (-sigmaf * state.dataEcoRoofMgr->epsilonf * Sigma -
478 532584 : sigmaf * state.dataEcoRoofMgr->epsilonf * state.dataEcoRoofMgr->epsilong * Sigma / EpsilonOne) *
479 532584 : pow_4(LeafTK) +
480 532584 : state.dataEcoRoofMgr->sheatf * (1.0 - 0.7 * sigmaf) * (Ta + DataGlobalConstants::KelvinConv) +
481 532584 : state.dataEcoRoofMgr->LAI * Rhoaf * Cf * Lef * Waf * rn * ((1.0 - 0.7 * sigmaf) / dOne) * qa +
482 266292 : state.dataEcoRoofMgr->LAI * Rhoaf * Cf * Lef * Waf * rn * (((0.6 * sigmaf * rn) / dOne) - 1.0) * (qsf - LeafTK * dqf) +
483 266292 : state.dataEcoRoofMgr->LAI * Rhoaf * Cf * Lef * Waf * rn * ((0.1 * sigmaf * Mg) / dOne) * (qsg - SoilTK * dqg);
484 798876 : P2 = 4.0 * (sigmaf * state.dataEcoRoofMgr->epsilonf * state.dataEcoRoofMgr->epsilong * Sigma) * pow_3(SoilTK) / EpsilonOne +
485 266292 : 0.1 * sigmaf * state.dataEcoRoofMgr->sheatf +
486 266292 : state.dataEcoRoofMgr->LAI * Rhoaf * Cf * Lef * Waf * rn * (0.1 * sigmaf * Mg) / dOne * dqg;
487 532584 : P3 = 4.0 *
488 532584 : (-sigmaf * state.dataEcoRoofMgr->epsilonf * Sigma -
489 532584 : (sigmaf * state.dataEcoRoofMgr->epsilonf * Sigma * state.dataEcoRoofMgr->epsilong) / EpsilonOne) *
490 532584 : pow_3(LeafTK) +
491 266292 : (0.6 * sigmaf - 1.0) * state.dataEcoRoofMgr->sheatf +
492 266292 : state.dataEcoRoofMgr->LAI * Rhoaf * Cf * Lef * Waf * rn * (((0.6 * sigmaf * rn) / dOne) - 1.0) * dqf;
493 :
494 : // T1G, T2G, & T3G corresponds to first, second & third terms of equation 38
495 : // in the main report.
496 : // as with the equations for vegetation the first term in the ground eqn in FASST has a
497 : // term starting with gamma_p --- if no precip this vanishes. Again, revisit this issue later.
498 :
499 798876 : T1G = (1.0 - sigmaf) * (RS * (1.0 - state.dataEcoRoofMgr->Alphag) + state.dataEcoRoofMgr->epsilong * Latm) -
500 532584 : (3.0 * (sigmaf * state.dataEcoRoofMgr->epsilonf * state.dataEcoRoofMgr->epsilong * Sigma) / EpsilonOne) * pow_4(LeafTK) -
501 266292 : 3.0 *
502 532584 : (-(1.0 - sigmaf) * state.dataEcoRoofMgr->epsilong * Sigma -
503 532584 : sigmaf * state.dataEcoRoofMgr->epsilonf * state.dataEcoRoofMgr->epsilong * Sigma / EpsilonOne) *
504 532584 : pow_4(SoilTK) +
505 532584 : state.dataEcoRoofMgr->sheatg * (1.0 - 0.7 * sigmaf) * (Ta + DataGlobalConstants::KelvinConv) +
506 532584 : Rhoag * Ce * Leg * Waf * Mg * ((1.0 - 0.7 * sigmaf) / dOne) * qa +
507 532584 : Rhoag * Ce * Leg * Waf * Mg * (0.1 * sigmaf * Mg / dOne - Mg) * (qsg - SoilTK * dqg) +
508 532584 : Rhoag * Ce * Leg * Waf * Mg * (0.6 * sigmaf * rn / dOne) * (qsf - LeafTK * dqf) + Qsoilpart1 +
509 266292 : Qsoilpart2 * (DataGlobalConstants::KelvinConv); // finished by T1G
510 :
511 532584 : T2G = 4.0 *
512 532584 : (-(1.0 - sigmaf) * state.dataEcoRoofMgr->epsilong * Sigma -
513 532584 : sigmaf * state.dataEcoRoofMgr->epsilonf * state.dataEcoRoofMgr->epsilong * Sigma / EpsilonOne) *
514 532584 : pow_3(SoilTK) +
515 532584 : (0.1 * sigmaf - 1.0) * state.dataEcoRoofMgr->sheatg + Rhoag * Ce * Leg * Waf * Mg * (0.1 * sigmaf * Mg / dOne - Mg) * dqg -
516 : Qsoilpart2;
517 :
518 798876 : T3G = (4.0 * (sigmaf * state.dataEcoRoofMgr->epsilong * state.dataEcoRoofMgr->epsilonf * Sigma) / EpsilonOne) * pow_3(LeafTK) +
519 532584 : 0.6 * sigmaf * state.dataEcoRoofMgr->sheatg + Rhoag * Ce * Leg * Waf * Mg * (0.6 * sigmaf * rn / dOne) * dqf;
520 :
521 266292 : LeafTK = 0.5 * (LeafTK + (P1 * T2G - P2 * T1G) / (-P3 * T2G + T3G * P2)); // take avg of old and new each iteration
522 266292 : SoilTK = 0.5 * (SoilTK + (P1 * T3G - P3 * T1G) / (-P2 * T3G + P3 * T2G)); // take avg of old and new each iteration
523 : // in earlier implementations we simply assigned new Leaf and Soil temps once (no iteration -- ecoloop =1,1, but
524 : // with LeafTK = ( P1*T2G - P2*T1G )/( -P3*T2G + T3G*P2 ) and
525 : // SoilTK =(SoilTK+( P1*T3G - P3*T1G )/( -P2*T3G + P3*T2G ). This iterative solution was initially attempted to
526 : // deal with stability issues associated with the CTF. It has virtually no impact on temperatures and it is not
527 : // clear how much it helped with stability problems. Eventually when CTF solution is replaced with a finite
528 : // difference scheme this loop structure should be removed.
529 :
530 : } // This loop does an iterative solution of the simultaneous equations
531 88764 : state.dataEcoRoofMgr->Qsoil =
532 88764 : -1.0 * (Qsoilpart1 - Qsoilpart2 * (SoilTK - DataGlobalConstants::KelvinConv)); // This is heat flux INTO top of the soil
533 88764 : state.dataEcoRoofMgr->Tfold = LeafTK - DataGlobalConstants::KelvinConv;
534 88764 : state.dataEcoRoofMgr->Tgold = SoilTK - DataGlobalConstants::KelvinConv;
535 :
536 : } // if firstecosurface (if not we do NOT need to recalculate ecoroof energybalance as all ecoroof surfaces MUST be the same
537 : // this endif was moved here from the if statement regarding whether we are looking at the first ecoroof surface or not.
538 :
539 443820 : state.dataHeatBalSurf->SurfOutsideTempHist(1)(SurfNum) = state.dataEcoRoofMgr->Tgold; // SoilTemperature
540 443820 : TempExt = state.dataEcoRoofMgr->Tgold;
541 443820 : }
542 :
543 4 : void initEcoRoofFirstTime(EnergyPlusData &state, int const SurfNum, int const ConstrNum)
544 : {
545 4 : auto &thisMat = state.dataMaterial->Material(state.dataConstruction->Construct(ConstrNum).LayerPoint(1));
546 4 : auto &thisEcoRoof = state.dataEcoRoofMgr;
547 :
548 4 : thisEcoRoof->EcoRoofbeginFlag = false;
549 :
550 4 : if (state.dataSurface->Surface(SurfNum).HeatTransferAlgorithm != DataSurfaces::HeatTransferModel::CTF)
551 0 : ShowSevereError(state,
552 : "initEcoRoofFirstTime: EcoRoof simulation but HeatBalanceAlgorithm is not ConductionTransferFunction(CTF). EcoRoof model "
553 : "currently works only with CTF heat balance solution algorithm.");
554 :
555 : // ONLY READ ECOROOF PROPERTIES IN THE FIRST TIME
556 4 : thisEcoRoof->Zf = thisMat.HeightOfPlants; // Plant height (m)
557 4 : thisEcoRoof->LAI = thisMat.LAI; // Leaf Area Index
558 4 : thisEcoRoof->Alphag = 1.0 - thisMat.AbsorpSolar; // albedo rather than absorptivity
559 4 : thisEcoRoof->Alphaf = thisMat.Lreflectivity; // Leaf Reflectivity
560 4 : thisEcoRoof->epsilonf = thisMat.LEmissitivity; // Leaf Emisivity
561 4 : thisEcoRoof->StomatalResistanceMin = thisMat.RStomata; // Leaf min stomatal resistance
562 4 : thisEcoRoof->epsilong = thisMat.AbsorpThermal; // Soil Emisivity
563 4 : thisEcoRoof->MoistureMax = thisMat.Porosity; // Max moisture content in soil
564 4 : thisEcoRoof->MoistureResidual = thisMat.MinMoisture; // Min moisture content in soil
565 4 : thisEcoRoof->Moisture = thisMat.InitMoisture; // Initial moisture content in soil
566 4 : thisEcoRoof->MeanRootMoisture = thisEcoRoof->Moisture; // DJS Oct 2007 Release --> all soil at same initial moisture for Reverse DD fix
567 :
568 4 : thisEcoRoof->SoilThickness = thisMat.Thickness; // Total thickness of soil layer (m)
569 :
570 4 : thisEcoRoof->FirstEcoSurf = SurfNum; // this determines WHEN to updatesoilProps
571 :
572 : // DJS NOVEMBER 2010 - Make calls to SetupOutput Variable to allow for reporting of ecoroof variables
573 12 : SetupOutputVariable(state,
574 : "Green Roof Soil Temperature",
575 : OutputProcessor::Unit::C,
576 4 : thisEcoRoof->Tg,
577 : OutputProcessor::SOVTimeStepType::Zone,
578 : OutputProcessor::SOVStoreType::State,
579 8 : "Environment");
580 12 : SetupOutputVariable(state,
581 : "Green Roof Vegetation Temperature",
582 : OutputProcessor::Unit::C,
583 4 : thisEcoRoof->Tf,
584 : OutputProcessor::SOVTimeStepType::Zone,
585 : OutputProcessor::SOVStoreType::State,
586 8 : "Environment");
587 12 : SetupOutputVariable(state,
588 : "Green Roof Soil Root Moisture Ratio",
589 : OutputProcessor::Unit::None,
590 4 : thisEcoRoof->MeanRootMoisture,
591 : OutputProcessor::SOVTimeStepType::Zone,
592 : OutputProcessor::SOVStoreType::State,
593 8 : "Environment");
594 12 : SetupOutputVariable(state,
595 : "Green Roof Soil Near Surface Moisture Ratio",
596 : OutputProcessor::Unit::None,
597 4 : thisEcoRoof->Moisture,
598 : OutputProcessor::SOVTimeStepType::Zone,
599 : OutputProcessor::SOVStoreType::State,
600 8 : "Environment");
601 12 : SetupOutputVariable(state,
602 : "Green Roof Soil Sensible Heat Transfer Rate per Area",
603 : OutputProcessor::Unit::W_m2,
604 4 : thisEcoRoof->sensibleg,
605 : OutputProcessor::SOVTimeStepType::Zone,
606 : OutputProcessor::SOVStoreType::State,
607 8 : "Environment");
608 12 : SetupOutputVariable(state,
609 : "Green Roof Vegetation Sensible Heat Transfer Rate per Area",
610 : OutputProcessor::Unit::W_m2,
611 4 : thisEcoRoof->sensiblef,
612 : OutputProcessor::SOVTimeStepType::Zone,
613 : OutputProcessor::SOVStoreType::State,
614 8 : "Environment");
615 12 : SetupOutputVariable(state,
616 : "Green Roof Vegetation Moisture Transfer Rate",
617 : OutputProcessor::Unit::m_s,
618 4 : thisEcoRoof->Vfluxf,
619 : OutputProcessor::SOVTimeStepType::Zone,
620 : OutputProcessor::SOVStoreType::State,
621 8 : "Environment");
622 12 : SetupOutputVariable(state,
623 : "Green Roof Soil Moisture Transfer Rate",
624 : OutputProcessor::Unit::m_s,
625 4 : thisEcoRoof->Vfluxg,
626 : OutputProcessor::SOVTimeStepType::Zone,
627 : OutputProcessor::SOVStoreType::State,
628 8 : "Environment");
629 12 : SetupOutputVariable(state,
630 : "Green Roof Vegetation Latent Heat Transfer Rate per Area",
631 : OutputProcessor::Unit::W_m2,
632 4 : thisEcoRoof->Lf,
633 : OutputProcessor::SOVTimeStepType::Zone,
634 : OutputProcessor::SOVStoreType::State,
635 8 : "Environment");
636 12 : SetupOutputVariable(state,
637 : "Green Roof Soil Latent Heat Transfer Rate per Area",
638 : OutputProcessor::Unit::W_m2,
639 4 : thisEcoRoof->Lg,
640 : OutputProcessor::SOVTimeStepType::Zone,
641 : OutputProcessor::SOVStoreType::State,
642 8 : "Environment");
643 :
644 12 : SetupOutputVariable(state,
645 : "Green Roof Cumulative Precipitation Depth",
646 : OutputProcessor::Unit::m,
647 4 : thisEcoRoof->CumPrecip,
648 : OutputProcessor::SOVTimeStepType::Zone,
649 : OutputProcessor::SOVStoreType::State,
650 8 : "Environment");
651 12 : SetupOutputVariable(state,
652 : "Green Roof Cumulative Irrigation Depth",
653 : OutputProcessor::Unit::m,
654 4 : thisEcoRoof->CumIrrigation,
655 : OutputProcessor::SOVTimeStepType::Zone,
656 : OutputProcessor::SOVStoreType::State,
657 8 : "Environment");
658 12 : SetupOutputVariable(state,
659 : "Green Roof Cumulative Runoff Depth",
660 : OutputProcessor::Unit::m,
661 4 : thisEcoRoof->CumRunoff,
662 : OutputProcessor::SOVTimeStepType::Zone,
663 : OutputProcessor::SOVStoreType::State,
664 8 : "Environment");
665 12 : SetupOutputVariable(state,
666 : "Green Roof Cumulative Evapotranspiration Depth",
667 : OutputProcessor::Unit::m,
668 4 : thisEcoRoof->CumET,
669 : OutputProcessor::SOVTimeStepType::Zone,
670 : OutputProcessor::SOVStoreType::State,
671 8 : "Environment");
672 12 : SetupOutputVariable(state,
673 : "Green Roof Current Precipitation Depth",
674 : OutputProcessor::Unit::m,
675 4 : thisEcoRoof->CurrentPrecipitation,
676 : OutputProcessor::SOVTimeStepType::Zone,
677 : OutputProcessor::SOVStoreType::Summed,
678 8 : "Environment");
679 12 : SetupOutputVariable(state,
680 : "Green Roof Current Irrigation Depth",
681 : OutputProcessor::Unit::m,
682 4 : thisEcoRoof->CurrentIrrigation,
683 : OutputProcessor::SOVTimeStepType::Zone,
684 : OutputProcessor::SOVStoreType::Summed,
685 8 : "Environment");
686 12 : SetupOutputVariable(state,
687 : "Green Roof Current Runoff Depth",
688 : OutputProcessor::Unit::m,
689 4 : thisEcoRoof->CurrentRunoff,
690 : OutputProcessor::SOVTimeStepType::Zone,
691 : OutputProcessor::SOVStoreType::Summed,
692 8 : "Environment");
693 12 : SetupOutputVariable(state,
694 : "Green Roof Current Evapotranspiration Depth",
695 : OutputProcessor::Unit::m,
696 4 : thisEcoRoof->CurrentET,
697 : OutputProcessor::SOVTimeStepType::Zone,
698 : OutputProcessor::SOVStoreType::Summed,
699 8 : "Environment");
700 4 : }
701 :
702 443820 : void initEcoRoof(EnergyPlusData &state, int const SurfNum, int const ConstrNum)
703 : {
704 : // Using/Aliasing
705 : using namespace DataEnvironment;
706 :
707 443820 : auto &thisMat = state.dataMaterial->Material(state.dataConstruction->Construct(ConstrNum).LayerPoint(1));
708 443820 : auto &thisEcoRoof = state.dataEcoRoofMgr;
709 443820 : auto &thisSurf = state.dataSurface->Surface(SurfNum);
710 :
711 : // DJS July 2007
712 : // Make sure the ecoroof module resets its conditions at start of EVERY warmup day and every new design day
713 : // for Reverse DD testing
714 443820 : if (state.dataGlobal->BeginEnvrnFlag || state.dataGlobal->WarmupFlag) {
715 380460 : thisEcoRoof->Moisture = thisMat.InitMoisture; // Initial moisture content in soil
716 380460 : thisEcoRoof->MeanRootMoisture = thisEcoRoof->Moisture; // Start the root zone moisture at the same value as the surface.
717 380460 : thisEcoRoof->Alphag = 1.0 - thisMat.AbsorpSolar; // albedo rather than absorptivity
718 : }
719 :
720 443820 : if (state.dataGlobal->BeginEnvrnFlag && thisEcoRoof->CalcEcoRoofMyEnvrnFlag) {
721 36 : thisEcoRoof->Tgold = OutDryBulbTempAt(state, thisSurf.Centroid.z); // OutDryBulbTemp initial guess
722 36 : thisEcoRoof->Tfold = OutDryBulbTempAt(state, thisSurf.Centroid.z); // OutDryBulbTemp initial guess
723 36 : thisEcoRoof->Tg = 10.0;
724 36 : thisEcoRoof->Tf = 10.0;
725 36 : thisEcoRoof->Vfluxf = 0.0;
726 36 : thisEcoRoof->Vfluxg = 0.0;
727 36 : thisEcoRoof->CumRunoff = 0.0;
728 36 : thisEcoRoof->CumET = 0.0;
729 36 : thisEcoRoof->CumPrecip = 0.0;
730 36 : thisEcoRoof->CumIrrigation = 0.0;
731 36 : thisEcoRoof->CurrentRunoff = 0.0;
732 36 : thisEcoRoof->CurrentET = 0.0;
733 36 : thisEcoRoof->CurrentPrecipitation = 0.0;
734 36 : thisEcoRoof->CurrentIrrigation = 0.0;
735 36 : thisEcoRoof->CalcEcoRoofMyEnvrnFlag = false;
736 : }
737 :
738 443820 : if (!state.dataGlobal->BeginEnvrnFlag) thisEcoRoof->CalcEcoRoofMyEnvrnFlag = true;
739 443820 : }
740 :
741 88764 : void UpdateSoilProps(EnergyPlusData &state,
742 : Real64 &Moisture,
743 : Real64 &MeanRootMoisture,
744 : Real64 const MoistureMax,
745 : Real64 const MoistureResidual,
746 : Real64 const SoilThickness,
747 : Real64 const Vfluxf, // Water mass flux from vegetation [m/s]
748 : Real64 const Vfluxg, // Water mass flux from soil surface [m/s]
749 : int const ConstrNum, // Indicator for construction index for the current surface
750 : Real64 &Alphag,
751 : [[maybe_unused]] int const unit, // unused1208
752 : [[maybe_unused]] Real64 const Tg, // unused1208
753 : [[maybe_unused]] Real64 const Tf, // unused1208
754 : [[maybe_unused]] Real64 const Qsoil // unused1208
755 : )
756 : {
757 : // SUBROUTINE INFORMATION
758 : // AUTHOR David Sailor
759 : // DATE WRITTEN Jan 2007
760 : // MODIFIED Stephen Forner, Portland State University (SF); 7/15/2010
761 : // RE-ENGINEERED na
762 :
763 : // PURPOSE OF THIS MODULE:
764 :
765 : // Track moisture input/output to ecoroof soil media (precipitation, irrigation, evapotranspiration, runoff)
766 : // Update soil thermal properties associated with variations in soil moisture and update CTF calculations
767 : // for the ecoroof construction layer.
768 :
769 : // METHODOLOGY EMPLOYED:
770 : // Define 2 soil layers (top and root). Moisture redistribution between top and root layers is proportional
771 : // to moisture differences. Input and Output of moisture are partitioned between layers.
772 : // Soil thermal properties vary based on non-dimensionalization of experimental data for 8 typical soils.
773 : // Specifically, THERMAL PROPERTY = Dry Value + (fraction of moisture content)*Wet Value
774 :
775 : // Using/Aliasing
776 : using namespace DataEnvironment;
777 : using namespace DataSurfaces;
778 :
779 : // SUBROUTINE PARAMETER DEFINITIONS:
780 : static Real64 const depth_fac((161240.0 * std::pow(2.0, -2.3)) / 60.0);
781 :
782 : // Soil Parameters from Reference listed in the code:
783 88764 : Real64 constexpr alpha(23.0); // These parameters are empirical constants
784 88764 : Real64 constexpr n(1.27); // These parameters are empirical constants
785 88764 : Real64 constexpr lambda(0.5); // These parameters are empirical constants
786 : // This is another parameter of the soil which describes the soil conductivity at the saturation point (m/s)
787 88764 : Real64 constexpr SoilConductivitySaturation(5.157e-7);
788 :
789 : // INTERFACE BLOCK SPECIFICATIONS:
790 : // na
791 :
792 : // DERIVED TYPE DEFINITIONS:
793 : // na
794 :
795 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
796 : Real64 RatioMax;
797 : Real64 RatioMin;
798 : Real64 MoistureDiffusion; // moisture transport down from near-surface to root zone
799 : Real64 SoilConductivity; // Moisture dependent conductivity to be fed back into CTF Calculator
800 : Real64 SoilSpecHeat; // Moisture dependent Spec. Heat to be fed back into CTF Calculator
801 : Real64 SoilAbsorpSolar; // Moisture dependent Solar absorptance (1-albedo)
802 : Real64 SoilDensity; // Moisture dependent density to be fed back into CTF Calculator
803 : Real64 SatRatio;
804 : Real64 TestRatio; // Ratio to determine if timestep change in properties is too abrupt for CTF
805 : Real64 AvgMoisture; // Average soil moisture over depth of ecoroof media
806 : int index1;
807 :
808 : // NOTE: As Energyplus calls the energy balance manager (and hence CalcEcoroof)
809 : // once for each surface within each zone that has an ecoroof
810 : // --- the CALCECOROOF routine is called multiple times within each time step
811 : // So, it is important that the UpdateSoilProps subroutine is called ONLY ONCE for each time step!!!
812 :
813 : // Recall Moisture = near-surface moisture value (m^3/m^3)
814 : // Recall MeanRootMoisture = root zone moisture value (m^3/m^3)
815 :
816 : // DJS 2009 set the ratiomax and ratiomin values in the code now (rather than as parameters) so that
817 : // DJS 2009 we can link them to timesteps and make these limits apply to actual RATES...
818 : // DJS 2009 reasonable rates are +/- 10% change in properties per 15 minute period... Otherwise we have
819 : // DJS 2009 stability issues.
820 : // DJS 2011 FEB - Since we no longer use CTF with soil-dependent properties (Do not RECALL INITCONDUCTION...
821 : // DJS 2011 FEB - we may be able to get away with NO limits on rates of change when using CFD routine.
822 : // DJS 2011 FEB - for now we stick with 20% per quarter hour.
823 88764 : RatioMax = 1.0 + 0.20 * state.dataGlobal->MinutesPerTimeStep / 15.0;
824 88764 : RatioMin = 1.0 - 0.20 * state.dataGlobal->MinutesPerTimeStep / 15.0;
825 :
826 88764 : if (state.dataEcoRoofMgr->UpdatebeginFlag) {
827 :
828 : // SET dry values that NEVER CHANGE
829 4 : state.dataEcoRoofMgr->DryCond = state.dataMaterial->Material(state.dataConstruction->Construct(ConstrNum).LayerPoint(1)).Conductivity;
830 4 : state.dataEcoRoofMgr->DryDens = state.dataMaterial->Material(state.dataConstruction->Construct(ConstrNum).LayerPoint(1)).Density;
831 4 : state.dataEcoRoofMgr->DryAbsorp = state.dataMaterial->Material(state.dataConstruction->Construct(ConstrNum).LayerPoint(1)).AbsorpSolar;
832 4 : state.dataEcoRoofMgr->DrySpecHeat = state.dataMaterial->Material(state.dataConstruction->Construct(ConstrNum).LayerPoint(1)).SpecHeat;
833 :
834 : // DETERMINE RELATIVE THICKNESS OF TWO LAYERS OF SOIL (also unchanging)
835 4 : if (SoilThickness > 0.12) {
836 4 : state.dataEcoRoofMgr->TopDepth = 0.06; // For now use 6cm as top depth - everything else is root layer
837 : } else {
838 0 : state.dataEcoRoofMgr->TopDepth = 0.5 * SoilThickness; // In unusual case of very thin soil make topdepth half of total
839 : }
840 : // This loop outputs the minimum number of time steps needed to keep the solution stable
841 : // The equation is minimum timestep in seconds=161240*((number of layers)**(-2.3))*(Total thickness of the soil)**2.07
842 4 : if (state.dataMaterial->Material(state.dataConstruction->Construct(ConstrNum).LayerPoint(1)).EcoRoofCalculationMethod == 2) {
843 2 : Real64 const depth_limit(depth_fac * std::pow(state.dataEcoRoofMgr->TopDepth + state.dataEcoRoofMgr->RootDepth, 2.07));
844 2 : for (index1 = 1; index1 <= 20; ++index1) {
845 2 : if (double(state.dataGlobal->MinutesPerTimeStep / index1) <= depth_limit) break;
846 : }
847 2 : if (index1 > 1) {
848 0 : ShowWarningError(state, "CalcEcoRoof: Too few time steps per hour for stability.");
849 0 : if (ceil(60 * index1 / state.dataGlobal->MinutesPerTimeStep) <= 60) {
850 0 : ShowContinueError(
851 : state,
852 0 : format("...Entered Timesteps per hour=[{}], Change to some value greater than or equal to [{}] for assured stability.",
853 0 : state.dataGlobal->NumOfTimeStepInHour,
854 0 : 60 * index1 / state.dataGlobal->MinutesPerTimeStep));
855 0 : ShowContinueError(state, "...Note that EnergyPlus has a maximum of 60 timesteps per hour");
856 0 : ShowContinueError(state,
857 : "...The program will continue, but if the simulation fails due to too low/high temperatures, instability "
858 : "here could be the reason.");
859 : } else {
860 0 : ShowContinueError(state,
861 0 : format("...Entered Timesteps per hour=[{}], however the required frequency for stability [{}] is over the "
862 : "EnergyPlus maximum of 60.",
863 0 : state.dataGlobal->NumOfTimeStepInHour,
864 0 : 60 * index1 / state.dataGlobal->MinutesPerTimeStep));
865 0 : ShowContinueError(state, "...Consider using the simple moisture diffusion calculation method for this application");
866 0 : ShowContinueError(state,
867 : "...The program will continue, but if the simulation fails due to too low/high temperatures, instability "
868 : "here could be the reason.");
869 : }
870 : }
871 : }
872 :
873 4 : state.dataEcoRoofMgr->RootDepth = SoilThickness - state.dataEcoRoofMgr->TopDepth;
874 : // Next create a timestep in seconds
875 4 : state.dataEcoRoofMgr->TimeStepZoneSec = state.dataGlobal->MinutesPerTimeStep * 60.0;
876 :
877 4 : state.dataEcoRoofMgr->UpdatebeginFlag = false;
878 : }
879 :
880 88764 : state.dataEcoRoofMgr->CurrentRunoff = 0.0; // Initialize current time step runoff as it is used in several spots below...
881 :
882 : // FIRST Subtract water evaporated by plants and at soil surface
883 88764 : Moisture -= (Vfluxg)*state.dataGlobal->MinutesPerTimeStep * 60.0 / state.dataEcoRoofMgr->TopDepth; // soil surface evaporation
884 88764 : MeanRootMoisture -= (Vfluxf)*state.dataGlobal->MinutesPerTimeStep * 60.0 / state.dataEcoRoofMgr->RootDepth; // plant extraction from root zone
885 :
886 : // NEXT Update evapotranspiration summary variable for print out
887 88764 : state.dataEcoRoofMgr->CurrentET = (Vfluxg + Vfluxf) * state.dataGlobal->MinutesPerTimeStep * 60.0; // units are meters
888 88764 : if (!state.dataGlobal->WarmupFlag) {
889 12672 : state.dataEcoRoofMgr->CumET += state.dataEcoRoofMgr->CurrentET;
890 : }
891 :
892 : // NEXT Add Precipitation to surface soil moisture variable (if a schedule exists)
893 88764 : Moisture += state.dataEcoRoofMgr->CurrentPrecipitation / state.dataEcoRoofMgr->TopDepth; // x (m) evenly put into top layer
894 :
895 88764 : if (!state.dataGlobal->WarmupFlag) {
896 12672 : state.dataEcoRoofMgr->CumPrecip += state.dataEcoRoofMgr->CurrentPrecipitation;
897 : }
898 :
899 : // NEXT Add Irrigation to surface soil moisture variable (if a schedule exists)
900 88764 : state.dataEcoRoofMgr->CurrentIrrigation = 0.0; // first initialize to zero
901 88764 : state.dataWaterData->Irrigation.ActualAmount = 0.0;
902 88764 : if (state.dataWaterData->Irrigation.ModeID == DataWater::IrrigationMode::IrrSchedDesign) {
903 0 : state.dataEcoRoofMgr->CurrentIrrigation = state.dataWaterData->Irrigation.ScheduledAmount; // units of m
904 0 : state.dataWaterData->Irrigation.ActualAmount = state.dataEcoRoofMgr->CurrentIrrigation;
905 : // elseif (Irrigation%ModeID ==IrrSmartSched .and. moisture .lt. 0.4d0*MoistureMax) then
906 177524 : } else if (state.dataWaterData->Irrigation.ModeID == DataWater::IrrigationMode::IrrSmartSched &&
907 88760 : Moisture < state.dataWaterData->Irrigation.IrrigationThreshold * MoistureMax) {
908 : // Smart schedule only irrigates when scheduled AND the soil is less than 40% saturated
909 0 : state.dataEcoRoofMgr->CurrentIrrigation = state.dataWaterData->Irrigation.ScheduledAmount; // units of m
910 0 : state.dataWaterData->Irrigation.ActualAmount = state.dataEcoRoofMgr->CurrentIrrigation;
911 : } else { // no schedule
912 88764 : if (state.dataWaterData->RainFall.ModeID == DataWater::RainfallMode::EPWPrecipitation) {
913 40334 : state.dataEcoRoofMgr->CurrentIrrigation = 0; // units of m
914 40334 : state.dataWaterData->Irrigation.ActualAmount = state.dataEcoRoofMgr->CurrentIrrigation;
915 : }
916 : }
917 :
918 88764 : Moisture += state.dataEcoRoofMgr->CurrentIrrigation / state.dataEcoRoofMgr->TopDepth; // irrigation in (m)/timestep put into top layer
919 :
920 88764 : if ((state.dataEnvrn->RunPeriodEnvironment) && (!state.dataGlobal->WarmupFlag)) {
921 0 : state.dataEcoRoofMgr->CumIrrigation += state.dataEcoRoofMgr->CurrentIrrigation;
922 : // aggregate to monthly for reporting
923 0 : int month = state.dataEnvrn->Month;
924 0 : state.dataEcoRoofMgr->MonthlyIrrigation.at(month - 1) += state.dataWaterData->Irrigation.ActualAmount * 1000.0;
925 : }
926 :
927 : // Note: If soil top layer gets a massive influx of rain &/or irrigation some of
928 : // the water will simply run right off the top and not penetrate at all!
929 : // At the present time this limit is fairly small due to some minor stability issues
930 : // in EnergyPlus. If the moisture changes too rapidly the code cannot handle the rapid changes in
931 : // surface characteristics and heat fluxes. The result that I've noticed is a non-physical fluctation
932 : // in ground surface temperature that oscillates up to 10 deg C from one hour to the next until the
933 : // code catches up. The temporary solution is to simply limit how much moisture can enter the soil
934 : // in any time step to 0.5"/hour. In the future this might be fixed by running with finer time steps
935 : // or by using a finite difference temperature solution rather than the CTF.
936 : // I suspect that 15 minute intervals may be needed. Another option is to have an internal moisture
937 : // overflow bin that will hold extra moisture and then distribute it in subsequent hours. This way the
938 : // soil still gets the same total moisture... it is just distributed over a longer period.
939 177528 : if (state.dataEcoRoofMgr->CurrentIrrigation + state.dataEcoRoofMgr->CurrentPrecipitation >
940 88764 : 0.5 * 0.0254 * state.dataGlobal->MinutesPerTimeStep / 60.0) {
941 0 : state.dataEcoRoofMgr->CurrentRunoff = state.dataEcoRoofMgr->CurrentIrrigation + state.dataEcoRoofMgr->CurrentPrecipitation -
942 0 : (0.5 * 0.0254 * state.dataGlobal->MinutesPerTimeStep / 60.0);
943 : // If we get here then TOO much moisture has already been added to soil (must now subtract excess)
944 0 : Moisture -= state.dataEcoRoofMgr->CurrentRunoff / state.dataEcoRoofMgr->TopDepth; // currently any incident moisture in excess of 1/4 "
945 : // per hour
946 : // simply runs off the top of the soil.
947 : }
948 : // Now, if top layer is beyond saturation... the excess simply runs off without penetrating into the lower
949 : // layers.
950 88764 : if (Moisture > MoistureMax) {
951 0 : state.dataEcoRoofMgr->CurrentRunoff += (Moisture - MoistureMax) * state.dataEcoRoofMgr->TopDepth;
952 0 : Moisture = MoistureMax;
953 : }
954 :
955 88764 : if (state.dataMaterial->Material(state.dataConstruction->Construct(ConstrNum).LayerPoint(1)).EcoRoofCalculationMethod == 1) {
956 :
957 : // THE SECTION BELOW WAS THE INITIAL MOISTURE DISTRIBUTION MODEL.
958 : // Any line with "!-" was code. A line with "!" was just a comment. This is done in case this code needs to be resurected in the future.
959 : // See below this commented out code for the new moisture distribution model.
960 : //*********************************************************************************************************
961 : //*********************************************************************************************************
962 : // NEXT Redistribute moisture based on moisture diffusion.
963 : // The effective diffusivities should be revisted when better moisture transport data in ecoroof soils are
964 : // available.
965 : // Here the diffusion rate is in units of [1/s]
966 : // A value of 0.0001 would be ~ 36% / hour
967 : // A value of 0.00005 would be ~ 18%/hour change in moisture
968 : // A value of 0.000025 gives about a 9%/hour change in moisture
969 : // a value of 0.0000125 gives 5%/hour...
970 : // Note: This formulation allows moisture to have a directional bias - ie, it can sink into soil easier than
971 : // it can be brought up.
972 8094 : if (Moisture > MeanRootMoisture) {
973 : // move moisture from top layer into root zone
974 2994 : MoistureDiffusion = min((MoistureMax - MeanRootMoisture) * state.dataEcoRoofMgr->RootDepth,
975 2994 : (Moisture - MeanRootMoisture) * state.dataEcoRoofMgr->TopDepth);
976 2994 : MoistureDiffusion = max(0.0, MoistureDiffusion); // Safety net to keep positive (not needed?)
977 : // at this point moistureDiffusion is in units of (m)/timestep
978 2994 : MoistureDiffusion *= 0.00005 * state.dataGlobal->MinutesPerTimeStep * 60.0;
979 2994 : Moisture -= MoistureDiffusion / state.dataEcoRoofMgr->TopDepth;
980 2994 : MeanRootMoisture += MoistureDiffusion / state.dataEcoRoofMgr->RootDepth;
981 5100 : } else if (MeanRootMoisture > Moisture) {
982 : // move moisture to top layer from root zone
983 1862 : MoistureDiffusion =
984 1862 : min((MoistureMax - Moisture) * state.dataEcoRoofMgr->TopDepth, (MeanRootMoisture - Moisture) * state.dataEcoRoofMgr->RootDepth);
985 1862 : MoistureDiffusion = max(0.0, MoistureDiffusion); // Safety net (not needed?)
986 : // at this point moistureDiffusion is in units of (m)/timestep
987 1862 : MoistureDiffusion *= 0.00001 * state.dataGlobal->MinutesPerTimeStep * 60.0;
988 1862 : Moisture += MoistureDiffusion / state.dataEcoRoofMgr->TopDepth;
989 1862 : MeanRootMoisture -= MoistureDiffusion / state.dataEcoRoofMgr->RootDepth;
990 : }
991 : } else {
992 : //********************************************************************************************************
993 : //********************************************************************************************************
994 : // THE SECTION ABOVE WAS THE MOISTURE DISTRIBUTION MODEL. REPLACED SF-7/21/2010
995 :
996 : // NEXT redistribute the moisture in the soil based on:
997 : // Marcel G Schaap and Martinus Th van Genuchten, 2006, 'A modified Maulem-van
998 : // Genuchten Formulation for Improved Description of the Hydraulic Conductivity Near Saturation'.
999 : // Written in MATLAB by Vishal Sharma (of Portland State) and modified for FORTRAN by Stephen Forner Summer 2010
1000 : // This model is based on curve fit data that describes the capillary motion of the water in combination with the gravitational
1001 : // forces on the water.
1002 : // This set of equations is unstable if the time step is larger than given by a curve fit equation. This first DO loop is to
1003 : // see how many time steps are needed to be stable.
1004 : // This method of moisture distribution relies on variables which are experimentally determined: alpha, lambda, n and the
1005 : // hydraulic conductivity at saturation.
1006 :
1007 : // Now, solve for the soil parameters for of the top soil layer
1008 :
1009 80670 : state.dataEcoRoofMgr->RelativeSoilSaturationTop = (Moisture - MoistureResidual) / (MoistureMax - MoistureResidual);
1010 80670 : if (state.dataEcoRoofMgr->RelativeSoilSaturationTop < 0.0001) {
1011 0 : if (state.dataEcoRoofMgr->ErrIndex == 0) {
1012 0 : ShowWarningMessage(state,
1013 0 : format("EcoRoof: UpdateSoilProps: Relative Soil Saturation Top Moisture <= 0.0001, Value=[{:.5R}].",
1014 0 : state.dataEcoRoofMgr->RelativeSoilSaturationTop));
1015 0 : ShowContinueError(state, "Value is set to 0.0001 and simulation continues.");
1016 0 : ShowContinueError(state, "You may wish to increase the number of timesteps to attempt to alleviate the problem.");
1017 : }
1018 0 : ShowRecurringWarningErrorAtEnd(state,
1019 : "EcoRoof: UpdateSoilProps: Relative Soil Saturation Top Moisture < 0. continues",
1020 0 : state.dataEcoRoofMgr->ErrIndex,
1021 0 : state.dataEcoRoofMgr->RelativeSoilSaturationTop,
1022 0 : state.dataEcoRoofMgr->RelativeSoilSaturationTop);
1023 0 : state.dataEcoRoofMgr->RelativeSoilSaturationTop = 0.0001;
1024 : }
1025 80670 : state.dataEcoRoofMgr->SoilHydroConductivityTop =
1026 161340 : SoilConductivitySaturation * std::pow(state.dataEcoRoofMgr->RelativeSoilSaturationTop, lambda) *
1027 80670 : pow_2(1.0 - std::pow(1.0 - std::pow(state.dataEcoRoofMgr->RelativeSoilSaturationTop, n / (n - 1.0)), (n - 1.0) / n));
1028 80670 : state.dataEcoRoofMgr->CapillaryPotentialTop =
1029 80670 : (-1.0 / alpha) * std::pow(std::pow(1.0 / state.dataEcoRoofMgr->RelativeSoilSaturationTop, n / (n - 1.0)) - 1.0, 1.0 / n);
1030 :
1031 : // Then the soil parameters for the root soil layer
1032 80670 : state.dataEcoRoofMgr->RelativeSoilSaturationRoot = (MeanRootMoisture - MoistureResidual) / (MoistureMax - MoistureResidual);
1033 80670 : state.dataEcoRoofMgr->SoilHydroConductivityRoot =
1034 161340 : SoilConductivitySaturation * std::pow(state.dataEcoRoofMgr->RelativeSoilSaturationRoot, lambda) *
1035 80670 : pow_2(1.0 - std::pow(1.0 - std::pow(state.dataEcoRoofMgr->RelativeSoilSaturationRoot, n / (n - 1.0)), (n - 1.0) / n));
1036 80670 : state.dataEcoRoofMgr->CapillaryPotentialRoot =
1037 80670 : (-1.0 / alpha) * std::pow(std::pow(1.0 / state.dataEcoRoofMgr->RelativeSoilSaturationRoot, n / (n - 1.0)) - 1.0, 1.0 / n);
1038 :
1039 : // Next, using the soil parameters, solve for the soil moisture
1040 80670 : state.dataEcoRoofMgr->SoilConductivityAveTop =
1041 80670 : (state.dataEcoRoofMgr->SoilHydroConductivityTop + state.dataEcoRoofMgr->SoilHydroConductivityRoot) * 0.5;
1042 80670 : Moisture +=
1043 161340 : (state.dataEcoRoofMgr->TimeStepZoneSec / state.dataEcoRoofMgr->TopDepth) *
1044 161340 : ((state.dataEcoRoofMgr->SoilConductivityAveTop *
1045 242010 : (state.dataEcoRoofMgr->CapillaryPotentialTop - state.dataEcoRoofMgr->CapillaryPotentialRoot) / state.dataEcoRoofMgr->TopDepth) -
1046 80670 : state.dataEcoRoofMgr->SoilConductivityAveTop);
1047 :
1048 : // Now limit the soil from going over the moisture maximum and takes excess to create runoff
1049 80670 : if (Moisture >= MoistureMax) { // This statement makes sure that the top layer is not over the moisture maximum for the soil.
1050 0 : Moisture = 0.9999 * MoistureMax; // then it takes any moisture over the maximum amount and makes it runoff
1051 0 : state.dataEcoRoofMgr->CurrentRunoff += (Moisture - MoistureMax * 0.9999) * state.dataEcoRoofMgr->TopDepth;
1052 : }
1053 :
1054 : // Now make sure that the soil does not go below the moisture minimum
1055 80670 : if (Moisture <= (1.01 * MoistureResidual)) {
1056 0 : Moisture = 1.01 * MoistureResidual;
1057 : }
1058 :
1059 : // Next, solve the parameters for the bottom layer
1060 80670 : state.dataEcoRoofMgr->SoilConductivityAveRoot = state.dataEcoRoofMgr->SoilHydroConductivityRoot;
1061 :
1062 : // Now make sure the rate of liquid leaving the soil is more than one drop per hour
1063 80670 : if ((state.dataEcoRoofMgr->SoilConductivityAveRoot * 3600.0) <= (2.33e-7)) {
1064 80670 : state.dataEcoRoofMgr->SoilConductivityAveRoot = 0.0;
1065 : }
1066 :
1067 : // Using the parameters above, distribute the Root Layer moisture
1068 80670 : state.dataEcoRoofMgr->TestMoisture = MeanRootMoisture;
1069 80670 : MeanRootMoisture +=
1070 161340 : (state.dataEcoRoofMgr->TimeStepZoneSec / state.dataEcoRoofMgr->RootDepth) *
1071 161340 : ((state.dataEcoRoofMgr->SoilConductivityAveTop *
1072 242010 : (state.dataEcoRoofMgr->CapillaryPotentialTop - state.dataEcoRoofMgr->CapillaryPotentialRoot) / state.dataEcoRoofMgr->RootDepth) +
1073 161340 : state.dataEcoRoofMgr->SoilConductivityAveTop - state.dataEcoRoofMgr->SoilConductivityAveRoot);
1074 :
1075 : // Limit the moisture from going over the saturation limit and create runoff:
1076 80670 : if (MeanRootMoisture >= MoistureMax) {
1077 0 : MeanRootMoisture = 0.9999 * MoistureMax;
1078 0 : state.dataEcoRoofMgr->CurrentRunoff += (Moisture - MoistureMax * 0.9999) * state.dataEcoRoofMgr->RootDepth;
1079 : }
1080 :
1081 : // Limit the soil from going below the soil saturation limit:
1082 80670 : if (MeanRootMoisture <= (1.01 * MoistureResidual)) {
1083 0 : MeanRootMoisture = 1.01 * MoistureResidual;
1084 : }
1085 :
1086 : // Next, track runoff from the bottom of the soil:
1087 80670 : state.dataEcoRoofMgr->CurrentRunoff += state.dataEcoRoofMgr->SoilConductivityAveRoot * state.dataEcoRoofMgr->TimeStepZoneSec;
1088 :
1089 : //~~~END SF EDITS
1090 : }
1091 :
1092 : // NEXT Limit moisture values to saturation (create RUNOFF that we can track)
1093 : // CurrentRunoff is sum of "overwatering" in a timestep and excess moisture content
1094 88764 : if (!state.dataGlobal->WarmupFlag) {
1095 12672 : state.dataEcoRoofMgr->CumRunoff += state.dataEcoRoofMgr->CurrentRunoff;
1096 : }
1097 :
1098 88764 : if (MeanRootMoisture <= MoistureResidual * 1.00001) {
1099 0 : Moisture -= (MoistureResidual * 1.00001 - MeanRootMoisture) * state.dataEcoRoofMgr->RootDepth / state.dataEcoRoofMgr->TopDepth;
1100 : // If the plant has extracted more moisture than is in the root zone soil, then make it come from
1101 : // the top layer rather than the root zone... unless top is also dry...
1102 0 : if (Moisture < MoistureResidual * 1.00001) Moisture = MoistureResidual * 1.00001;
1103 0 : MeanRootMoisture = MoistureResidual * 1.00001; // Need to make sure we do not divide by zero later.
1104 : }
1105 :
1106 : // ------------------------------------------------------------------------------------------
1107 : // Having updated moisture values now we modify soil thermal properties based on moisture content
1108 :
1109 : // NOTE: Variables SoilAbsorpSolar, SoilDensity, SoilSpecHeat, and SoilConductivity are the values
1110 : // that the soil properties OUGHT to attain for the current moisture level. These values are then
1111 : // moderated using the TestRatio variable so that from one time step to the next no variable
1112 : // changes by more than a certain percentage (typically 5-20%).
1113 :
1114 : // Note wet soil absorptance is generally 25-50% higher than dry soil absorptance (assume linear)
1115 177528 : SoilAbsorpSolar = state.dataEcoRoofMgr->DryAbsorp +
1116 88764 : (0.92 - state.dataEcoRoofMgr->DryAbsorp) * (Moisture - MoistureResidual) / (MoistureMax - MoistureResidual);
1117 : // Limit solar absorptivity to 95% so soil abledo is always above 5%
1118 88764 : if (SoilAbsorpSolar > 0.95) SoilAbsorpSolar = 0.95;
1119 : // Limit solar absorptivity to greater than 20% so that albedo is always less than 80%
1120 88764 : if (SoilAbsorpSolar < 0.20) SoilAbsorpSolar = 0.20;
1121 :
1122 : // Need to use for albedo in CalcEcoroof
1123 88764 : TestRatio = (1.0 - SoilAbsorpSolar) / Alphag;
1124 88764 : if (TestRatio > RatioMax) TestRatio = RatioMax;
1125 88764 : if (TestRatio < RatioMin) TestRatio = RatioMin;
1126 88764 : Alphag *= TestRatio; // included 1.0 - to make it albedo rather than SW absorptivity
1127 : // Note wet soil density is calculated by simply adding the mass of water...
1128 88764 : AvgMoisture = (state.dataEcoRoofMgr->RootDepth * MeanRootMoisture + state.dataEcoRoofMgr->TopDepth * Moisture) / SoilThickness;
1129 88764 : SoilDensity = state.dataEcoRoofMgr->DryDens + (AvgMoisture - MoistureResidual) * 990.0;
1130 : // Note 990 kg/m^3 is water density and the moisture is depth-averaged
1131 :
1132 : // Note wet soil has specific heat that is 40% higher than dry soil (assume linear)
1133 : // OLD :: SoilSpecHeat = DrySpecHeat*(1.0d0+ 0.4d0*(AvgMoisture-MoistureResidual)/(MoistureMax-MoistureResidual))
1134 : // This is now based on Melos Hagos's results for C (March 2009)
1135 : // SoilSpecHeat = DrySpecHeat + 3.09*(AvgMoisture) CLEARLY in ERROR BY FACTOR of 1000
1136 : // DJS - Melos report has Spec = Cdry + 1.9 theta (where C is in kJ/kg/K), so...
1137 88764 : SoilSpecHeat = state.dataEcoRoofMgr->DrySpecHeat + 1900.0 * AvgMoisture;
1138 :
1139 : // Note wet soil has thermal conductivity that is up to 3 times that of dry soil ...
1140 : // For now simply let it DOUBLE over the range of moisture
1141 :
1142 : // Note wet soil has thermal conductivity that is up to 3 times that of dry soil ...
1143 : // OLD :: SoilConductivity = DryCond* (1.0d0 + 1.0d0 * (AvgMoisture-MoistureResidual)/(MoistureMax-MoistureResidual))
1144 : // This is now based on Melos Hagos's results for k/kdry (March 2009)
1145 88764 : SatRatio = (AvgMoisture - MoistureResidual) / (MoistureMax - MoistureResidual);
1146 88764 : SoilConductivity = (state.dataEcoRoofMgr->DryCond / 1.15) * (1.45 * std::exp(4.411 * SatRatio)) / (1.0 + 0.45 * std::exp(4.411 * SatRatio));
1147 : // DJS 2009 - note, this allows the actual conductivity to dip a little below the dry value... corresponding to
1148 : // DJS 2009 - "bone dry" if you will, when moisture --> residual value.
1149 :
1150 : // HERE WE RE-RUN THE CONDUCTION TRANSFER FUNCTION (CTF) CALCULATIONS
1151 :
1152 : // NOTE: CTF uses the original dataMaterial.Material( )%xxx variables, so must update these for current construction and
1153 : // time step...
1154 : // TestRatio variable is available just in case there are stability issues. If so, we can limit the amount
1155 : // by which soil properties are allowed to vary in one time step (10% in example below).
1156 :
1157 88764 : TestRatio = SoilConductivity / state.dataMaterial->Material(state.dataConstruction->Construct(ConstrNum).LayerPoint(1)).Conductivity;
1158 88764 : if (TestRatio > RatioMax) TestRatio = RatioMax;
1159 88764 : if (TestRatio < RatioMin) TestRatio = RatioMin;
1160 88764 : state.dataMaterial->Material(state.dataConstruction->Construct(ConstrNum).LayerPoint(1)).Conductivity *= TestRatio;
1161 88764 : SoilConductivity = state.dataMaterial->Material(state.dataConstruction->Construct(ConstrNum).LayerPoint(1)).Conductivity;
1162 :
1163 88764 : TestRatio = SoilDensity / state.dataMaterial->Material(state.dataConstruction->Construct(ConstrNum).LayerPoint(1)).Density;
1164 88764 : if (TestRatio > RatioMax) TestRatio = RatioMax;
1165 88764 : if (TestRatio < RatioMin) TestRatio = RatioMin;
1166 88764 : state.dataMaterial->Material(state.dataConstruction->Construct(ConstrNum).LayerPoint(1)).Density *= TestRatio;
1167 88764 : SoilDensity = state.dataMaterial->Material(state.dataConstruction->Construct(ConstrNum).LayerPoint(1)).Density;
1168 :
1169 88764 : TestRatio = SoilSpecHeat / state.dataMaterial->Material(state.dataConstruction->Construct(ConstrNum).LayerPoint(1)).SpecHeat;
1170 88764 : if (TestRatio > RatioMax) TestRatio = RatioMax;
1171 88764 : if (TestRatio < RatioMin) TestRatio = RatioMin;
1172 88764 : state.dataMaterial->Material(state.dataConstruction->Construct(ConstrNum).LayerPoint(1)).SpecHeat *= TestRatio;
1173 88764 : SoilSpecHeat = state.dataMaterial->Material(state.dataConstruction->Construct(ConstrNum).LayerPoint(1)).SpecHeat;
1174 :
1175 : // Now call InitConductionTransferFunction with the ConstrNum as the argument. As long as the argument is
1176 : // non-zero InitConductionTransferFunction will ONLY update this construction. If the argument is 0 it will
1177 : // rerun the ENTIRE InitConductionTransferFunction on all constructions (as in initial code start-up mode).
1178 :
1179 : // DJS The following works for most simulations, but has stability issues in some cases.
1180 : // DJS - in the present version it seems best to NOT update soil thermal properties.
1181 : // Call InitConductionTransferFunctions(ConstrNum)
1182 : // DJS In future revision we will implement these modified soil thermal properties in the finite difference
1183 : // solution scheme.
1184 :
1185 : // DJS - Note the following write/format statements should be commented out prior to releasing within EnergyPlus
1186 : // DJS - they are handy in doing hourly validation/evaluations, though, so leave them in for future development.
1187 :
1188 : // write(unit,799) DayofYear, HourOfDay, Qsoil,Tg, Tf, Moisture, MeanRootMoisture,CumPrecip &
1189 : // ,CumET,CumRunoff, CumIrrigation, SoilDensity, SoilSpecHeat,SoilConductivity,Alphag
1190 : // 799 format(' ',I3,' ',I3,' ',' ',f9.3,' ',f6.2,' ',f6.2,' ',f5.3,' ',f5.3,' ',f6.4, ' ' &
1191 : // f7.3, ' ', f7.3, ' ',f7.3, ' ',f6.1,' ',f7.1,' ',f6.3,' ',f6.2)
1192 88764 : }
1193 :
1194 88764 : void CalculateEcoRoofSolar(EnergyPlusData &state,
1195 : Real64 &RS, // Solar on roof (assumed horizontal)
1196 : Real64 &f1, // Solar term for Stomatal Resistance
1197 : int const SurfNum)
1198 : {
1199 :
1200 : // Use SOLCOS(3) here to calculate RS since as stated in comments above this is only done for one surface currently.
1201 : // So, it is better to assume that the roof is flat until multiple surfaces can be handled. Then, the next line can
1202 : // use state.dataHeatBal->SurfCosIncAng(state.dataGlobal->HourOfDay, state.dataGlobal->TimeStep, SurfNum) instead of
1203 : // SOLCOS(3).
1204 177528 : RS = max(state.dataEnvrn->SOLCOS(3), 0.0) * state.dataEnvrn->BeamSolarRad +
1205 88764 : state.dataSolarShading->SurfAnisoSkyMult(SurfNum) * state.dataEnvrn->DifSolarRad;
1206 88764 : Real64 f1inv = min(1.0, (0.004 * RS + 0.005) / (0.81 * (0.004 * RS + 1.0)));
1207 88764 : f1 = 1.0 / f1inv;
1208 88764 : }
1209 :
1210 : } // namespace EcoRoofManager
1211 :
1212 2313 : } // namespace EnergyPlus
|
