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