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