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