Merge pull request #10040 from NREL/10022ShadingSched

Fix overlapping shadow calculations for shading surface transmittance >0 and <1

doc/engineering-reference/src/climate-sky-and-solar-shading-calculations/shading-module.tex

@@ -406,39 +406,93 @@ \subsubsection{Overlapping Shadows}\label{overlapping-shadows}
 \caption{Complex Overlapping Condition \protect \label{fig:complex-overlapping-condition}}
 \end{figure}
 
-If two shadows overlap the receiving surface, they may also overlap each other as in Figure~\ref{fig:multiple-shadow-overlaps}. The vertices of this overlap can be computed.~ The areas of all overlaps can be computed.~ The total sunlit area can be expressed as the sum of all polygon areas given a proper sign on each of the areas.
+If two shadows overlap the receiving surface, they may also overlap each other as in Figure~\ref{fig:multiple-shadow-overlaps}. The vertices of this overlap can be computed, and the areas of all overlaps can be computed. For opaque shadows, the total sunlit area can be expressed as the sum of all polygon areas using the sign convention shown in Table \ref{table:surface-area-characteristic-convention}:
 
-The following convention was adopted:
+\begin{equation}
+SunlitArea = {\sum\limits_{i = 1}^n {A_i}}
+\end{equation}
 
-% table 24
-\begin{longtable}[c]{@{}ll@{}}
-\caption{Surface / Area Characteristic / Convention \label{table:surface-area-characteristic-convention}} \tabularnewline
+Each shadow's area is subtracted from the receiving surface and so on through multiple overlaps where the sign of the overlap area is the product of the signs of the overlapping areas. For the shadows in Figure~\ref{fig:multiple-shadow-overlaps}, start with the receiving surface area A, subtract shadow area B, then subtract shadow area C. Because shadows B and C overlap, the overlap area D has been subtracted twice, so add back area D to get the final sunlit area on the receiving surface.
+
+
+\begin{longtable}[c]{>{\raggedright}p{2.9in}p{1.5in}p{1.59in}}
+\caption{Overlapping Shadow Surface Area Convention \label{table:surface-area-characteristic-convention}} \tabularnewline
 \toprule
-Surface Characteristic & Area Convention \tabularnewline
+Surface Characteristic & Area Convention & Sunlit Area \tabularnewline
 \midrule
 \endfirsthead
 
-\caption[]{Surface / Area Characteristic / Convention} \tabularnewline
+\caption[]{Overlapping Shadow Surface Area Convention} \tabularnewline
 \toprule
-Surface Characteristic & Area Convention \tabularnewline
+Surface Characteristic & Area Convention & Sunlit Area \tabularnewline
 \midrule
 \endhead
 
-receiving surface & positive (A) \tabularnewline
-overlap between shadow and receiving & negative (B  \&  C) \tabularnewline
-overlap between two shadows & positive (D) \tabularnewline
+receiving surface & positive A & A \tabularnewline
+overlap between shadow and receiving & negative B & A-B \tabularnewline
+overlap between shadow and receiving & negative C & A-B-C \tabularnewline
+overlap between two shadows & positive D & A-B-C+D \tabularnewline
 \bottomrule
 \end{longtable}
 
-and so on through multiple overlaps where the sign of the overlap area is the product of the signs of the overlapping areas.
-
 \begin{figure}[hbtp] % fig 47
 \centering
 \includegraphics[width=0.9\textwidth, height=0.9\textheight, keepaspectratio=true]{media/image642.png}
 \caption{Multiple Shadow Overlaps \protect \label{fig:multiple-shadow-overlaps}}
 \end{figure}
 
-Partially transparent shadowing surfaces can also be modeled by giving a transparency ($\tau$) to every shadowing polygon. Let $\tau$ of the receiving polygon be one. Then the $\tau$ of every overlap of polygons i and j is the product of $\tau$\(_{i}\)and $\tau$\(_{j}\). The shaded area is then computed by summing A\(_{i}\)*(1 - $\tau$\(_{i}\)) for all overlap polygons.
+Partially transparent shadowing surfaces can also be modeled by giving a transparency ($\tau$) to every shadowing polygon. The sunlit area is computed by:
+
+\begin{equation}
+SunlitArea = A_1 + {\sum\limits_{i = 2}^n {A_i}*(1-\tau_i)}
+\end{equation}
+
+The actual $\tau$ of overlapping polygons i and j is the product of $\tau_i$ and $\tau_j$. For example, if a wall is fully shaded by two partially transparent shades, one with a transmittance of 0.8 and one with a transmittance of 0.5, the resulting sunlit fraction would be $0.8 * 0.5 = 0.4$. Because each shadow is first applied with its own transmittance, the ``transmittance'' for the overlap correction ($\tau_k$) in the sunlit area summation is derived from the individual surface transmittance values:
+
+\begin{equation}
+\tau_k = (\tau_{i} + \tau_{j}) - (\tau_{i} * \tau_{j})
+\end{equation}
+
+where $\tau_k$ is the correction ``transmittance''  for the shadow overlap area (D in Figure~\ref{fig:multiple-shadow-overlaps}).
+
+For example, using the overlapping shadows in Figure~\ref{fig:multiple-shadow-overlaps}, the sunlit area can be calculated two ways as shown in Table~\ref{table:sunlit-area-calculation-overlapping-shadows}. Method 1 avoids any overlaps by calculating adjusted areas. Method 2 uses the EnergyPlus summation approach and shows the transmittance correction calculation for the overlap area to arrive at the same result. Stepping through Method 2: start with the full surface area A, subtract the shaded area $B*(1-\tau_b)$, subtract the shaded area $C*(1-\tau_c)$. At this point, the overlapping area D has been subtracted twice, once using $\tau_b$ and once using $\tau_c$. To correct for this, area D is added back in with the adjusted ``transmittance'' using the equation shown above for $\tau_k$. Once the sunlit area is known for each surface, the sunlit fraction is calculated, which is $13.56/20=0.678$ for this example.
+
+\begin{longtable}[c]{>{\raggedright}p{1.2in}p{1.5in}p{2.0in}p{2.0in}}
+\caption{Sunlit Area Calculations with Overlapping Partially Transmitting Shadows\label{table:sunlit-area-calculation-overlapping-shadows}} \tabularnewline
+\toprule
+Region & Area & Transmittance & Sunlit Area \tabularnewline
+\midrule
+\endfirsthead
+
+\caption[]{Sunlit Area Calculations with Overlapping Partially Transmitting Shadows} \tabularnewline
+\toprule
+Region & Area & Transmittance & Sunlit Area \tabularnewline
+\midrule
+\endhead
+
+Assumptions: \tabularnewline
+A & 20 & 1.0 \tabularnewline
+B & 10 & 0.8 \tabularnewline
+C & 8 & 0.4 \tabularnewline
+D & 3 &  \tabularnewline
+\midrule
+Method 1 & Adjusted Areas:\tabularnewline
+$A-B-C+D$ & $20-10-8+3 = 5$ & $1.0$               & $5*1.0 = 5.00$ \tabularnewline
+$B-D$     & $10-3=7$        & $0.8$               & $7*0.8 = 5.60$ \tabularnewline
+$C-D$     & $8-3=5$         & $0.4$               & $5*0.4 = 2.00$ \tabularnewline
+$D$       & $3$             & $0.8 * 0.4 = 0.32$  & $3*0.32 = 0.96$ \tabularnewline
+Total     & $20$            &                     & $13.56$ \tabularnewline
+\midrule
+Method 2 & Full Areas:\tabularnewline
+$A$   & $20$         & n/a                             & $20.00$ \tabularnewline
+$B$   & $-10$        & $0.8$                           & $-10*(1-0.8) = -2.00$ \tabularnewline
+$C$   & $-8$         & $0.4$                           & $-8*(1-0.4) = -4.80$ \tabularnewline
+$D$   & $3$          & $(0.8+0.4)-(0.8 * 0.4) = 0.88$  & $+3*(1-0.88) = 0.36$ \tabularnewline
+Total &              &                                 & $13.56$ \tabularnewline
+\bottomrule
+\end{longtable}
+
+
 
 It is easy to determine the sunlit area of a window once all the shadow and overlap vertices on the wall have been computed. Consider wall 2 of Figure~\ref{fig:overall-shadowing-scheme-depiction}. First, the wall is considered a simple rectangle and the window on it is ignored. The shadow overlapping is performed and the sunlit portion of the gross wall area is computed. Then the window rectangle is overlapped with the shadow to determine its sunlit area. The sunlit area of the window is subtracted from the gross wall sunlit area to determine the net wall sunlit area. During this calculation it is not necessary to recompute the shadows, because they were precisely determined on the wall.
 

src/EnergyPlus/SolarShading.cc

@@ -4637,8 +4637,6 @@ void DeterminePolygonOverlap(EnergyPlusData &state,
 
     // SUBROUTINE INFORMATION:
     //       AUTHOR         Legacy Code
-    //       DATE WRITTEN
-    //       MODIFIED       na
     //       RE-ENGINEERED  Lawrie, Oct 2000
 
     // PURPOSE OF THIS SUBROUTINE:
@@ -4667,20 +4665,11 @@ void DeterminePolygonOverlap(EnergyPlusData &state,
     // REFERENCES:
     // BLAST/IBLAST code, original author George Walton
 
-    // Using/Aliasing
-
-    int N;    // Loop index
-    int NV1;  // Number of vertices of figure NS1
-    int NV2;  // Number of vertices of figure NS2
-    int NV3;  // Number of vertices of figure NS3 (the overlap of NS1 and NS2)
-    int NIN1; // Number of vertices of NS1 within NS2
-    int NIN2; // Number of vertices of NS2 within NS1
-
-    // Check for exceeding array limits.
 #ifdef EP_Count_Calls
     ++state.dataTimingsData->NumDetPolyOverlap_Calls;
 #endif
 
+    // Check for exceeding array limits.
     if (NS3 > state.dataSolarShading->MaxHCS) {
 
         state.dataSolarShading->OverlapStatus = TooManyFigures;
@@ -4705,9 +4694,11 @@ void DeterminePolygonOverlap(EnergyPlusData &state,
     }
 
     state.dataSolarShading->OverlapStatus = PartialOverlap;
-    NV1 = state.dataSolarShading->HCNV(NS1);
-    NV2 = state.dataSolarShading->HCNV(NS2);
-    NV3 = 0;
+    int NV1 = state.dataSolarShading->HCNV(NS1); // Number of vertices of figure NS1
+    int NV2 = state.dataSolarShading->HCNV(NS2); // Number of vertices of figure NS2
+    int NV3 = 0;                                 // Number of vertices of figure NS3 (the overlap of NS1 and NS2)
+    int NIN1 = 0;                                // Number of vertices of NS1 within NS2
+    int NIN2 = 0;                                // Number of vertices of NS2 within NS1
 
     if (!state.dataSysVars->SutherlandHodgman) {
         INCLOS(state, NS1, NV1, NS2, NV2, NV3, NIN1); // Find vertices of NS1 within NS2.
@@ -4740,14 +4731,14 @@ void DeterminePolygonOverlap(EnergyPlusData &state,
         CLIPPOLY(state, NS1, NS2, NV1, NV2, NV3);
     }
 
-    if (NV3 < state.dataSolarShading->MaxHCV && NS3 <= state.dataSolarShading->MaxHCS) {
+    if (NV3 < state.dataSolarShading->MaxHCV) {
 
         if (!state.dataSysVars->SutherlandHodgman) {
             ORDER(state, NV3, NS3); // Put vertices in clockwise order.
         } else {
             assert(equal_dimensions(state.dataSolarShading->HCX, state.dataSolarShading->HCY));
             int l = state.dataSolarShading->HCX.index(NS3, 1);
-            for (N = 1; N <= NV3; ++N, ++l) {
+            for (int N = 1; N <= NV3; ++N, ++l) {
                 state.dataSolarShading->HCX[l] = nint64(state.dataSolarShading->XTEMP(N)); // [ l ] == ( N, NS3 )
                 state.dataSolarShading->HCY[l] = nint64(state.dataSolarShading->YTEMP(N));
             }
@@ -4759,20 +4750,23 @@ void DeterminePolygonOverlap(EnergyPlusData &state,
         if (std::abs(state.dataSolarShading->HCAREA(NS3)) * HCMULT < std::abs(state.dataSolarShading->HCAREA(NS1))) {
             state.dataSolarShading->OverlapStatus = NoOverlap;
         } else {
-            if (state.dataSolarShading->HCAREA(NS1) * state.dataSolarShading->HCAREA(NS2) > 0.0)
+            if (state.dataSolarShading->HCAREA(NS1) * state.dataSolarShading->HCAREA(NS2) > 0.0) {
                 state.dataSolarShading->HCAREA(NS3) = -state.dataSolarShading->HCAREA(NS3); // Determine sign of area of overlap
-            Real64 const HCT_1(state.dataSolarShading->HCT(NS1));
-            Real64 const HCT_2(state.dataSolarShading->HCT(NS2));
-            Real64 HCT_3(HCT_2 * HCT_1); // Determine transmission of overlap
-            if (HCT_2 >= 0.5 && HCT_1 >= 0.5) {
-                if (HCT_2 != 1.0 && HCT_1 != 1.0) {
-                    HCT_3 = 1.0 - HCT_3;
-                }
             }
-            state.dataSolarShading->HCT(NS3) = HCT_3;
+            Real64 const HCT_1 = state.dataSolarShading->HCT(NS1);
+            Real64 const HCT_2 = state.dataSolarShading->HCT(NS2);
+            if (HCT_2 == 1.0 || HCT_1 == 1.0) {
+                state.dataSolarShading->HCT(NS3) = HCT_1 * HCT_2;
+            } else {
+                // Determine transmission of overlap which corrects for prior shadows
+                // The resulting transmission of overlapping shadows is HCT_1 * HCT_2
+                // Shadows with HCT_1 and HCT_2 have already been applied in the overlapping area
+                // so the correction is the difference between (HCT_1+HCT_2) and (HCT_1*HCT_2)
+                state.dataSolarShading->HCT(NS3) = (HCT_1 + HCT_2) - HCT_1 * HCT_2;
+            }
         }
 
-    } else if (NV3 > state.dataSolarShading->MaxHCV) {
+    } else {
 
         state.dataSolarShading->OverlapStatus = TooManyVertices;
 
@@ -4791,26 +4785,6 @@ void DeterminePolygonOverlap(EnergyPlusData &state,
             state.dataSolarShading->TrackTooManyVertices(state.dataSolarShading->NumTooManyVertices).SurfIndex2 =
                 state.dataSolarShading->CurrentSurfaceBeingShadowed;
         }
-
-    } else if (NS3 > state.dataSolarShading->MaxHCS) {
-
-        state.dataSolarShading->OverlapStatus = TooManyFigures;
-
-        if (!state.dataSolarShading->TooManyFiguresMessage && !state.dataGlobal->DisplayExtraWarnings) {
-            ShowWarningError(state,
-                             format("DeterminePolygonOverlap: Too many figures [>{}]  detected in an overlap calculation. Use "
-                                    "Output:Diagnostics,DisplayExtraWarnings; for more details.",
-                                    state.dataSolarShading->MaxHCS));
-            state.dataSolarShading->TooManyFiguresMessage = true;
-        }
-
-        if (state.dataGlobal->DisplayExtraWarnings) {
-            state.dataSolarShading->TrackTooManyFigures.redimension(++state.dataSolarShading->NumTooManyFigures);
-            state.dataSolarShading->TrackTooManyFigures(state.dataSolarShading->NumTooManyFigures).SurfIndex1 =
-                state.dataSolarShading->CurrentShadowingSurface;
-            state.dataSolarShading->TrackTooManyFigures(state.dataSolarShading->NumTooManyFigures).SurfIndex2 =
-                state.dataSolarShading->CurrentSurfaceBeingShadowed;
-        }
     }
 }
 

src/EnergyPlus/SolarShading.hh

@@ -371,18 +371,18 @@ struct SolarShadingData : BaseGlobalStruct
     Array1D<Real64> SurfMultCircumSolar;        // Contribution to eff sky view factor from circumsolar brightening
     Array1D<Real64> SurfMultHorizonZenith;      // Contribution to eff sky view factor from horizon or zenith brightening
 
-    int FBKSHC;                     // HC location of first back surface
-    int FGSSHC;                     // HC location of first general shadowing surface
-    int FINSHC;                     // HC location of first back surface overlap
-    int FRVLHC;                     // HC location of first reveal surface
-    int FSBSHC;                     // HC location of first subsurface
+    int FBKSHC = 0;                 // HC location of first back surface
+    int FGSSHC = 0;                 // HC location of first general shadowing surface
+    int FINSHC = 0;                 // HC location of first back surface overlap
+    int FRVLHC = 0;                 // HC location of first reveal surface
+    int FSBSHC = 0;                 // HC location of first subsurface
     int LOCHCA = 0;                 // Location of highest data in the HC arrays
-    int NBKSHC;                     // Number of back surfaces in the HC arrays
-    int NGSSHC;                     // Number of general shadowing surfaces in the HC arrays
-    int NINSHC;                     // Number of back surface overlaps in the HC arrays
-    int NRVLHC;                     // Number of reveal surfaces in HC array
-    int NSBSHC;                     // Number of subsurfaces in the HC arrays
-    bool CalcSkyDifShading;         // True when sky diffuse solar shading is
+    int NBKSHC = 0;                 // Number of back surfaces in the HC arrays
+    int NGSSHC = 0;                 // Number of general shadowing surfaces in the HC arrays
+    int NINSHC = 0;                 // Number of back surface overlaps in the HC arrays
+    int NRVLHC = 0;                 // Number of reveal surfaces in HC array
+    int NSBSHC = 0;                 // Number of subsurfaces in the HC arrays
+    bool CalcSkyDifShading = false; // True when sky diffuse solar shading is
     int ShadowingCalcFrequency = 0; // Frequency for Shadowing Calculations
     int ShadowingDaysLeft = 0;      // Days left in current shadowing period