V23.2.0-IOFreeze: IDD and transition fixes

design/FY2023/NFP-EarthTube-1DEnhancement.md

@@ -1,11 +1,13 @@
-Earth Tube 1-D Conduction Enhancement
+Earth Tube 1-D Conduction Enhancement (Design Document)
 ================
 
 **Rick Strand, University of Illinois at Urbana-Champaign**
 
  - April 2023
  - Revision Date (Version 2): May 11, 2023
  - Revision Date (Version 3): May 12, 2023
+ - Revision Date (Version 4): June 7, 2023 (converted to a design document, minor grammatical edits, addition of the Code Development setion)
+ - Revision Date (Version 5): June 27, 2023 (switch to LU decomposition)
 
 
 ## Justification for New Feature ##
@@ -17,7 +19,9 @@ The improvement of the earth tube model will be done in stages.  The first step
 
 Technicalities Meeting (5/3/23): team review the document briefly and stated the desire for this to get additional feedback in the next technicalities meeting.  Two initial requests from the team included modifying the proposal to factor in the impact of the soil at the depth of the earth tube since this would be more realistic of the situation (as in, the earth tube in not infinite in width) and to contact John Nelson who has expertise in earth tube modeling and has requested improved earth tube simulation capabilities in EnergyPlus in the past.  This update includes a modification of the finite difference scheme to address the concern regarding what is happening at the depth of the earth tube itself.  In addition, GitHub was used to contact John Nelson and that conversation is contained in the comments for the GitHub issue regarding earth tubes (https://github.com/NREL/EnergyPlus/issues/6627).  John suggested that the model at some point include the ability to model temperature variation axially along the earth tube (something already logged as a potential future enhancement) and also to provide some improved controls where air might bypass the earth tube and use heat recovery in certain situations to save the earth tube’s ability to cool.  This is an interesting idea but may also be fairly complicated because it may deal with specific control algorithms and bridging between the zone and air loop simulations.  Such a concept is different from what is currently in EnergyPlus and there doesn’t appear to be another model in E+ that bridges with air flow at the zone and air loop level.  Such a concept is definitely beyond the scope of this initial enhancement.
 
-Slack Discussion (5/11/23): Neal Kris brought up two issues in this NFP.  First, there is the issue of the infinitely wide domain which is unrealistic.  Second, there is the question of what air temperature is being used for heat transfer between the earth tube and the air passing through it.  The second issue will be resolved by adding the discussion of the effectiveness-NTU heat exchanger model that is being used.  This will be similar in the assumptions as the heat exchanger algorithm used in the low temperature radiant system model (solid side with a “constant” temperature and a fluid side which has a temperature which varies throughout).  The first issue will be addressed by modifying the solution domain to have a user-defined width with adiabatic boundary conditions on either side.  In reality, at some distance from the earth tube there will be negligible heat transfer toward/away from the earth tube.  Making this a user-defined input parameter will allow the width to be studied.
+Slack Discussion (5/11/23): Neal Kruis brought up two issues in this NFP.  First, there is the issue of the infinitely wide domain which is unrealistic.  Second, there is the question of what air temperature is being used for heat transfer between the earth tube and the air passing through it.  The second issue will be resolved by adding the discussion of the effectiveness-NTU heat exchanger model that is being used.  This will be similar in the assumptions as the heat exchanger algorithm used in the low temperature radiant system model (solid side with a “constant” temperature and a fluid side which has a temperature which varies throughout).  The first issue will be addressed by modifying the solution domain to have a user-defined width with adiabatic boundary conditions on either side.  In reality, at some distance from the earth tube there will be negligible heat transfer toward/away from the earth tube.  Making this a user-defined input parameter will allow the width to be studied.
+
+Technicalities Meeting (6/14/23): Jason DeGraw expressed specific concern about using matrix inversion as part of the solution technique due to concerns about slow execution speeds.  He recommended that an LU decomposition technique be used as a much faster approach and noted that there were several already in EnergyPlus.  As a result, this document was modified to remove inversion as a solution approach and replace it with LU decomposition, using the existing routines in the Window Manager.
 
 ## Overview ##
 
@@ -37,7 +41,7 @@ Overall, the approach will be to add a new 1-D heat transfer model that will run
 
 The 1-D heat conduction model itself will be implemented using a 1-D implicit finite difference scheme.  The solution dimension will be bounded by a temperature just below grade at the top of the solution space and another temperature boundary condition at a depth below the earth tube.  The temperatures at these upper and lower boundaries will be set using the existing equation from the literature for undisturbed soil temperature.  At the sides of the solution space, adiabatic conditions will be assumed.  The user will be given the option to specify the overall thickness of the solution region, but there will only be heat transfer vertically (none horizontally since it will be adiabatic on the “sides”).
 
-The finite difference grid will include an upper portion from the upper temperature boundary to the node just above the earth tube and also a lower portion from the node just below the earth tube to the lower boundary condition.  The node at the earth tube itself will also be modeled as an “all soil” node but also include a connection to the earth tube itself.  In a one dimension model like the one proposed for this enhancement, it is more realistic to model the level at the earth tube to be soil because in all directions there will be soil and not an entire layer of air/earth tube.  Thus, this earth tube node will capture the conditions of the soil at that level, and a heat exchanger model that connects to the  temperature of the earth tube node will be used, and checks will be made to insure that the heat loss of the air is equal to the total heat gain of the soil.
+The finite difference grid will include an upper portion from the upper temperature boundary to the node just above the earth tube and also a lower portion from the node just below the earth tube to the lower boundary condition.  The node at the earth tube itself will also be modeled as an “all soil” node but also include a connection to the earth tube itself.  In a one dimensional model like the one proposed for this enhancement, it is more realistic to model the level at the earth tube to be soil because in all directions there will be soil and not an entire layer of air/earth tube.  Thus, this earth tube node will capture the conditions of the soil at that level, and a heat exchanger model that connects to the  temperature of the earth tube node will be used, and checks will be made to insure that the heat loss of the air is equal to the total heat gain of the soil.
 
 The solution for this finite difference model will use an implicit scheme.  Implicit schemes are known for being inherently stable.  Since it is possible that the number of nodes may be small and thus the nodes large, it will be important to avoid potential stability issues in the solution.  The user will be given some flexibility in the input to specify the number of nodes above and below the earth tube.  While some reasonable limits will be applied to those inputs to avoid too few or too many nodes, the control of the number of nodes will also help during the testing and verification process to make sure that the model is producing results that make sense and are believable.  An example of the node layout is shown in the figure below.
 
@@ -202,7 +206,7 @@ This parameter sets the number of nodes below the earth tube, between the earth
 This parameter sets the dimensionless distance above the earth tube for the solution space.  The maximum value is 1.0, and the minimum value is 0.25.  When this parameter is set to 1.0, the upper boundary is set to be half of the diameter below the ground surface and the solution space thickness above the earth tube is the depth of the earth tube minus the earth tube diameter.  This maximum distance (earth tube depth minus diameter) is multiplied by this parameter to constrain the solution space to less than the maximum (when the parameter is less than 1.0).    The default value for this parameter is 1.0 (the maximum value).
 
 \paragraph{Field: Earth Tube Dimensionless Boundary Below}\label{field-earth-tube-dimensionless-boundary-below}
-This parameter sets the dimensionless distance below the earth tube for the solution space.  The maximum value is 1.0, and the minimum value is 0.25.  This parameter is interpretted in a similar fashion as the previous parameter where the depth of the solution space below the earth tube is determined by the maximum distance above the earth tube (earth tube depth minus diameter).  This allows the user to have different thickness for the modeled portion of the ground above and below the earth tube.  The default value for this parameter is 0.25 (the minimum value).
+This parameter sets the dimensionless distance below the earth tube for the solution space.  The maximum value is 1.0, and the minimum value is 0.25.  This parameter is interpreted in a similar fashion as the previous parameter where the depth of the solution space below the earth tube is determined by the maximum distance above the earth tube (earth tube depth minus diameter).  This allows the user to have different thickness for the modeled portion of the ground above and below the earth tube.  The default value for this parameter is 0.25 (the minimum value).
 
 \paragraph{Field: Earth Tube Dimensionless Solution Space Width}\label{field-earth-tube-dimensionless-solution-space-width}
 This parameter sets the dimensionless width of the solution space horizontally as a function of the earth tube radius as defined in the main earth tube input syntax.  The maximum value is 20.0, and the minimum value is 3.0.  The default value for this parameter is 4.0 which means that the width of the solution space is four times the radius.  In other words, this would include soil one radius length beyond the edges of the tube on either side of the earth tube.
@@ -374,6 +378,28 @@ There is already an existing EarthTubeSimpleTest.idf file that is used to test t
 
 Due to the fact that the two new fields in the EarthTube input are at the end of the current description and are both optional, no transition changes will be needed.  Existing files will continue to run without problem and will default to the simple solution method.  In addition, the new input syntax EarthTube:Parameters is also optional and only needed for the 1-D (Vertical) model.  Thus, current user files will not need this new input to continue to run their earth tube models.
 
+
+## Code Development ##
+
+The current code has one main subroutine that is tasked with the calculation of the earth tube performance.  First, after some preliminary initializations, the air flow related terms are calculated for the air side of the earth tube.  Then, the average earth temperature between the depth of the top of the earth tube and the bottom of the earth tube.  Next, given the air flow and thermal properties of the soil, the heat transfer coefficient including the various thermal resistances are calculated.  Finally, given the various conditions already calculated, an outlet air temperature is calculated for the earth tube.
+
+All of this methodology must be preserved since the goal of this implementation is to allow the modeing of both the existing model and the new enhancement so that the results can be compared side-by-side.  Most of the auxiliary functions that will be shared between the existing and the new model (air flow terms, heat transfer coefficients, etc.) will be converted to subroutines, reducing the clutter in the main calculation routine.  The existing model will also be converted into a function that will calculate the average temperature for various depths.  The main calculation routine will then chose one of two options: the existing model or the "vertical" model.  The existing model will continue to calculate the ground tempeature as is currently done, using the undisturbed temperature correlation using the new subroutine for average temperature between two depths.  The new "vertical" model will call out to a new subroutine that will implement the implicit finite difference scheme presented above.
+
+The new subroutine that will actually contain new code and not simply code moved from the main routine will deal with the calculation of the new vertical model being installed for this work.  The main purpose of this new subroutine will be the calculation of the node temperatures using the implicit finite difference scheme.  There will likely also need to be a one-time initialization routine for the various terms in the node equations since they will be based on thermo-physical properties of the unique situation (based on user input) and thus constant (assumption).
+
+Another initialization that only needs to be done once per day is the calculation of the soil temperatures.  In the current implementation, the soil temperature at the depth of the earth tube is calculated every zone time step.  This is unnecessary because the equation for the soil temperature at a given depth is set based on the day of the year, not on the hour of the day.  So, this initialization could be calculated once each day, increasing the efficiency of the routine.
+
+The new routine for the vertical model will have the following sequence/components:
+
+1. Initialize the node equation terms (once per simulation).  Also calculate the depth of the various nodes for the vertical solution.  Also make initial guesses for the various node temperatures based on the undisturbed ground temperature model.
+2. Calculate the air flow related terms for the air side of the earth tube (impacting the air heat balance of the zone).
+3. Calculate the soil temperature at the various node depths used for the model (once per day).  This includes the temperature at the depth of the tube as well as at the depth of the top and the bottom of the solution space.
+4. Calculate the heat transfer coefficients and thermal resistances that will be used for the heat exchanger model.  This will have an impact on the earth tube node equation since there will be heat exchange between the earth tube node, surrounding nodes, and the air passing through the earth tube.
+5. Calculate the node temperatures and outlet air temperature using an implicit finite difference scheme.  The node equations will be set up in a matrix format that follows the general layout of Ax = b.  This will be solved via the LU Decomposition strategy using existing EnergyPlus routines (see subroutine LUdecomposition and LUsolution in WindowManager.cc--will need to update so that it is no longer limited to an index of 10).  This will avoid the need for iteration to obtain a solution, but it will require the decomposition and solution to be calculated each time step because one of the terms in the A matrix will include the effectiveness (based on NTU as shown above) of the heat transfer between the air passing through the earth tube and the ground.  The remaining terms in the A matrix will include thermo-physical properties of the ground and the earth tube material itself which will all be assumed to be constant.  Experience with the existing earth tube model and the low temperature hydronic radiant system model has shown that in many situations, the effectiveness is very close or equal to 1.  So, during initializations, the decomposed matris for A when the effectiveness is 1 will be calculated and stored to hopefully reduce the time required to obtain a solution.  In addition, the decomposition of the A matrix for when the earth tube is not operating (zero flow conditions).  This is because the earth tube could spend considerable time being turned off (not needed during cooling season or not being run during heating season, for example).  However, due to the transient nature of the model, it will need to model the ground temperature nodes even when the earth tube is not running.  Storing A for zero flow and an effectiveness of one will increase data requirements but this will provide a faster solution since during those times, the A matrix will already have been LU decomposed and thus only the solution routine will need to be run.
+6. Record and report output from the vertical model.  Output for the vertical model will include all of the existing output for the current model.  In addition, output variables will be established for the node temperatures both above the earth tube, below the earth tube, and at the location of the earth tube.
+
+One potential improvememt that could be considered in the future for the model is the ability to specify starting ground temperatures rather than calculate them based on the undisturbed ground temperature as the initial starting point.  The point of this would be to allow the user to run multiple simulations sequentially.  The last node temperatures of one annual simulation could then be plugged into the next annual simulation as a way to get a better starting point or to see what happens when the earth tube is present for multiple years.  Eventually the starting conditions would be no longer be a factor, but at this point it is uncertain whether such a feature is necessary.  It is possible that a single annual simulation might be long enough to establish a realistic trend, but this is not known at this time and would be studied at a later date in a future FY.
+
 ## References ##
 
 Existing EnergyPlus Earth Tube Model, EnergyPlus Engineering Reference, EnergyPlus Input-Output Reference.