Fix Bugs in Emmel and Blocken Convection Algorithms #9129

mitchute · 2021-10-13T02:32:04Z

Pull request overview

Fixes Theta in Emmel Vertical Convection Coefficients Model is incorrectly determined #9106
Emmel and Blocken convection correlations for external wind-driven convection have some bugs in the correlation values and the calculated angle between the surface normal and the wind direction.
This adds new unit tests and corrects the errors.

Emmel et al. Table 3 - Convection correlations for walls and roofs

Blocken et al. Table 6 - Convection correlations for windward surfaces

mitchute

Ready for review.

mitchute · 2021-10-13T19:22:07Z

src/EnergyPlus/ConvectionCoefficients.cc

-        Theta = state.dataEnvrn->WindDir - Surface(SurfNum).Azimuth - 90.0; // TODO double check theta
+        Theta = CalcWindSurfaceTheta(state.dataEnvrn->WindDir, Surface(SurfNum).Azimuth);


This is the same theta angle calculation that was suspect before. Replacing it with the new unit-tested function.

mitchute · 2021-10-13T19:22:14Z

src/EnergyPlus/ConvectionCoefficients.cc

-        Theta = state.dataSurface->SurfOutWindDir(SurfNum) - Surface(SurfNum).Azimuth - 90.0; // TODO double check theta
+        Theta = CalcWindSurfaceTheta(state.dataEnvrn->WindDir, Surface(SurfNum).Azimuth);


This argument should be SurfOutWindDir?

mitchute · 2021-10-13T19:22:43Z

src/EnergyPlus/ConvectionCoefficients.cc

+Real64 CalcWindSurfaceTheta(Real64 const WindDir, Real64 const SurfAzimuth)
+{
+    // Computes the angle theta between the wind direction and the surface azimuth
+    // Should always be a value between 0-180 deg
+
+    Real64 theta = std::abs(WindDir - SurfAzimuth);
+    if (theta > 180) {
+        return abs(theta - 360);
+    } else {
+        return theta;
+    }
+}
+


New function for calculating the angle theta between the wind direction and the surface azimuth.

mitchute · 2021-10-13T19:23:01Z

src/EnergyPlus/ConvectionCoefficients.cc

@@ -8787,16 +8800,15 @@ Real64 CalcBlockenWindward(Real64 const WindAt10m,

    Real64 Theta; // angle between wind and surface azimuth

-    Theta = WindDir - SurfAzimuth - 90.0; // TODO double check theta
-    if (Theta > 180.0) Theta -= 360.0;
+    Theta = CalcWindSurfaceTheta(WindDir, SurfAzimuth);


Update suspect theta calculation.

Could have combined 8001-8003 as Real64 Theta = ...

Thanks. Fixed.

mitchute · 2021-10-13T19:27:26Z

src/EnergyPlus/ConvectionCoefficients.cc

    if (Theta <= 11.25) {
        Hf = 4.6 * std::pow(WindAt10m, 0.89);
-    } else if ((11.25 < Theta) && (Theta <= 33.75)) {
+    } else if (Theta <= 33.75) {
        Hf = 5.0 * std::pow(WindAt10m, 0.8);
-    } else if ((33.75 < Theta) && (Theta <= 56.25)) {
+    } else if (Theta <= 56.25) {
        Hf = 4.6 * std::pow(WindAt10m, 0.84);
-    } else if ((56.25 < Theta) && (Theta <= 100.0)) {
+    } else if (Theta <= 90.0) {


A few things here:

Don't need multiple bounding checks on theta at each condition. It should successively fall through until it hits the right one.

This is a "windward" function, so an upper bound of 100 deg seems inappropriate. 90 deg should be the max a "windward" surface could see.

upper bound, rounding, what if Theta comes out as 90.0000000001 ?

Agreed, that makes sense. So... 91? Back to 100?

And YOU get a unit test! That's up to you, I guess 91 should be adequate but 100 should never hurt anything, unless you want to use 91 and assert on an else.

Thanks. I backed it down to 100 deg. Who knows how often these things are actually used, but 100 isn't unreasonable.

mitchute · 2021-10-13T19:28:53Z

src/EnergyPlus/ConvectionCoefficients.cc

    if (Theta <= 22.5) {
        Hf = 5.15 * std::pow(WindAt10m, 0.81);
-    } else if ((22.5 < Theta) && (Theta <= 67.5)) {
+    } else if (Theta <= 67.5) {
        Hf = 3.34 * std::pow(WindAt10m, 0.84);
-    } else if ((67.5 < Theta) && (Theta <= 112.5)) {
+    } else if (Theta <= 112.5) {
        Hf = 4.78 * std::pow(WindAt10m, 0.71);
-    } else if ((112.5 < Theta) && (Theta <= 157.5)) {
+    } else if (Theta <= 157.5) {
        Hf = 4.05 * std::pow(WindAt10m, 0.77);
-    } else if ((157.5 < Theta) && (Theta <= 180.0)) {
+    } else if (Theta <= 180.0) {


Remove duplicative checks here.

mitchute · 2021-10-13T19:29:39Z

src/EnergyPlus/ConvectionCoefficients.cc

@@ -8894,19 +8905,18 @@ Real64 CalcEmmelRoof(EnergyPlusData &state,

    Real64 Theta; // angle between wind and surface azimuth

-    Theta = WindDir - LongAxisOutwardAzimuth - 90.0; // TODO double check theta
-    if (Theta > 180.0) Theta -= 360.0;
+    Theta = CalcWindSurfaceTheta(WindDir, LongAxisOutwardAzimuth);


Update suspect theta calculation here.

mitchute · 2021-10-13T19:31:08Z

src/EnergyPlus/ConvectionCoefficients.cc

-        Hf = 5.15 * std::pow(WindAt10m, 0.81);
-    } else if ((22.5 < Theta) && (Theta <= 67.5)) {
-        Hf = 3.34 * std::pow(WindAt10m, 0.84);
-    } else if ((67.5 < Theta) && (Theta <= 112.5)) {
-        Hf = 4.78 * std::pow(WindAt10m, 0.71);
-    } else if ((112.5 < Theta) && (Theta <= 157.5)) {
-        Hf = 4.05 * std::pow(WindAt10m, 0.77);
-    } else if ((157.5 < Theta) && (Theta <= 180.0)) {
-        Hf = 3.54 * std::pow(WindAt10m, 0.76);
+        Hf = 5.11 * std::pow(WindAt10m, 0.78);
+    } else if (Theta <= 67.5) {
+        Hf = 4.60 * std::pow(WindAt10m, 0.79);
+    } else if (Theta <= 112.5) {
+        Hf = 3.67 * std::pow(WindAt10m, 0.85);
+    } else if (Theta <= 157.5) {
+        Hf = 4.60 * std::pow(WindAt10m, 0.79);
+    } else if (Theta <= 180.0) {
+        Hf = 5.11 * std::pow(WindAt10m, 0.78);


Looks like a pasteo error. Based on the paper, these are the numbers for a vertical surface, not a roof surface.

mitchute · 2021-10-13T19:32:30Z

src/EnergyPlus/ConvectionCoefficients.cc

@@ -8919,7 +8929,7 @@ Real64 CalcEmmelRoof(EnergyPlusData &state,
                                      "CalcEmmelRoof: Convection model wind angle calculation suspect and high theta correlation",
                                      state.dataConvectionCoefficient->CalcEmmelRoofErrorIDX);

-        Hf = 3.54 * std::pow(WindAt10m, 0.76);
+        Hf = 3.67 * std::pow(WindAt10m, 0.85);


I'm not sure it can even make it here, TBH. I updated this to pull from the roof correlation, but this probably needs to be double-checked (by me) and removed if it's inaccessible.

I updated the error handling in these functions and fixed the theta angle calculation so that it's always bounded between 0-180.

mitchute · 2021-10-13T19:33:09Z

tst/EnergyPlus/unit/ConvectionCoefficients.unit.cc

@@ -2630,3 +2631,178 @@ TEST_F(ConvectionCoefficientsFixture, TestSetAdaptiveConvectionAlgoCoefficient)
    expected_curve = UtilityRoutines::FindItemInList("NUSSELTJURGESDUPCURVE", state->dataConvectionCoefficient->HcOutsideUserCurve);
    ASSERT_EQ(algorithm_identifier, expected_curve);
 }
+


You get a unit test, and you get a unit test. Everybody gets a unit test.

…alculations

mitchute · 2021-10-14T16:37:05Z

@rraustad Thanks for the reviews. I made a few updates based on your comments. Let me know if you have anything else.

rraustad · 2021-10-15T17:03:13Z

Yeah, looks like 0 diffs here and a final review did not find anything that stands out. Good to go.

Myoldmopar · 2021-11-01T16:04:12Z

Another one that looks fully ready. Thanks for all the review here @rraustad. I'll pull develop into this one as well for due diligence. Assuming it is clean, I'll get it in the queue.

Myoldmopar · 2021-11-02T16:29:33Z

I pulled develop in locally, ran all the tests and everything passed. I ran regressions against develop so that I didn't have to wait on CI and it came out perfect, no diffs. Merging this, thanks @mitchute

Fix Bugs in Emmel and Blocken Convection Algorithms

mitchute added 3 commits October 12, 2021 19:10

add new unit test for EmmelVertical, currently fails

7924f09

add CalcWindSurfaceTheta and unit test

714607d

fix EmmelVertical convection correlation

8537c55

mitchute added the Defect Includes code to repair a defect in EnergyPlus label Oct 13, 2021

mitchute added this to the EnergyPlus 2022.1 milestone Oct 13, 2021

mitchute requested a review from Myoldmopar October 13, 2021 02:32

mitchute self-assigned this Oct 13, 2021

mitchute marked this pull request as draft October 13, 2021 02:32

mitchute removed the request for review from Myoldmopar October 13, 2021 02:32

mitchute added 4 commits October 13, 2021 09:04

add unit test for emmel roof

f3f12b5

add unit test for BlockenWindward

0ab8110

correct surface theta for blocken windward

935ea17

fix a couple of other surface to wind angle angle calcs

4385d3e

mitchute marked this pull request as ready for review October 13, 2021 19:21

mitchute commented Oct 13, 2021

View reviewed changes

address reviewer comments, update error handling, bound theta angle c…

ae5c103

…alculations

cleanup warnings

7043253

Merge remote-tracking branch 'origin/develop' into emmel_conv

fb0a920

Myoldmopar merged commit 739c93d into develop Nov 2, 2021

Myoldmopar deleted the emmel_conv branch November 2, 2021 16:29

yujiex pushed a commit that referenced this pull request Jan 27, 2022

Merge pull request #9129 from NREL/emmel_conv

d8602da

Fix Bugs in Emmel and Blocken Convection Algorithms

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Fix Bugs in Emmel and Blocken Convection Algorithms #9129

Fix Bugs in Emmel and Blocken Convection Algorithms #9129

mitchute commented Oct 13, 2021 •

edited

Loading

mitchute left a comment

mitchute Oct 13, 2021

mitchute Oct 13, 2021

rraustad Oct 13, 2021 •

edited

Loading

mitchute Oct 13, 2021

mitchute Oct 13, 2021

rraustad Oct 13, 2021

mitchute Oct 14, 2021

mitchute Oct 13, 2021

rraustad Oct 13, 2021

mitchute Oct 13, 2021

rraustad Oct 13, 2021

mitchute Oct 14, 2021

mitchute Oct 13, 2021

mitchute Oct 13, 2021

rraustad Oct 13, 2021

mitchute Oct 14, 2021

mitchute Oct 13, 2021

mitchute Oct 13, 2021

mitchute Oct 14, 2021

mitchute Oct 13, 2021

mitchute commented Oct 14, 2021

rraustad commented Oct 15, 2021

Myoldmopar commented Nov 1, 2021

Myoldmopar commented Nov 2, 2021

		Theta = state.dataEnvrn->WindDir - Surface(SurfNum).Azimuth - 90.0; // TODO double check theta
		Theta = CalcWindSurfaceTheta(state.dataEnvrn->WindDir, Surface(SurfNum).Azimuth);

		Theta = state.dataSurface->SurfOutWindDir(SurfNum) - Surface(SurfNum).Azimuth - 90.0; // TODO double check theta
		Theta = CalcWindSurfaceTheta(state.dataEnvrn->WindDir, Surface(SurfNum).Azimuth);

Fix Bugs in Emmel and Blocken Convection Algorithms #9129

Fix Bugs in Emmel and Blocken Convection Algorithms #9129

Conversation

mitchute commented Oct 13, 2021 • edited Loading

Pull request overview

mitchute left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

rraustad Oct 13, 2021 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mitchute commented Oct 14, 2021

rraustad commented Oct 15, 2021

Myoldmopar commented Nov 1, 2021

Myoldmopar commented Nov 2, 2021

mitchute commented Oct 13, 2021 •

edited

Loading

rraustad Oct 13, 2021 •

edited

Loading