A New Steam Medium Model for Fast and Accurate Simulation of District Heating Systems

In U.S. district heating (DH) systems, steam is the most common heat transport medium. Simulation of large steam DH systems requires a computationally e ﬃ cient and accurate steam model; however, the commonly adopted IF97 water model is not fast enough for district-scale simulations, and its discontinuous thermodynamic property functions have shown to cause simulation problems and sometimes failure. To address these issues, this work introduces a new steam medium model for heating applications. Further, we avoid the common numerical challenges at the phase change boundary by adopting a novel split-medium approach. We implemented the model in the equation-based Modelica language and evaluated the accuracy and numerical performance across multiple scales: from fundamental thermodynamic properties to complete heating districts of several sizes. The new implementation produces similar accuracy as the IF97 model, but at a speed that is 5.6 – 9.3 times faster for complete DH simulations. Moreover, for the DH system with the IF97 model, the dimension of the nonlinear system of equations increases linearly with the number of buildings, but it stays constant with the new model; this is critically important for large scale system simulations. The new steam medium model is available open-source in the Modelica IBPSA and Buildings Libraries.


Introduction
District heating (DH) can effectively reduce CO 2 emissions and enable communities to leverage economizes of scale benefits [1]. The global DH market is large, with an estimated 80,000 systems in operation that distribute hot water and steam through almost 600,000 km of distribution pipes [2]. In the United States, steam is the most common medium for DH, representing 97% of all installations [3]. Beyond heating buildings alone, steam DH provides beneficial waste-heat recovery opportunities when coupled with power plants (e.g., combined heat and power) or other industrial systems (e.g., wastewater treatment, metal refineries, etc.). While the conversion of steam DH systems (first generation) to hot water DH systems (second generation and later) is an important mechanism to realize deep carbon savings [2], high-quality energy such as steam will likely still be present in future energy grids due to its mutualistic benefits with power generation and high-heat industrial processes such as those listed in [4].
Modeling and simulation can improve the design and operation of DH systems while also aiding the transition of legacy steam-based systems to newer technologies. The literature on modeling and simulation of hot water DH systems is abundant [5,6,7,8] but out of the scope of this paper. Further, simulation studies of steam heating systems for combined heat and power plants [9,10] and other industrial applications [11] are prevalent within DH contexts; however, to our best knowledge, past studies have yet to model the thermo-fluid dynamics of steam heating systems at district scale. District-scale simulation capabilities are critical for future grid-interactive efficient buildings and for our ability to reach carbon targets through financially viable means. To accomplish this task, the first step is to efficiently evaluate fundamental steam properties for fast and accurate DH simulation.
Summarized in Table 1, several steam medium models are used in scientific and industrial practices. The International Association for the Properties of Water and Steam (IAPWS) develops formulations for thermodynamic properties of water and steam for various applications. Based on the Helmholtz free energy, the IAPWS-95 model covers the largest range of pressure p and temperature T and is intended for scientific applications where accuracy is of utmost importance [12]. The IAPWS also released their industrial formulation (IF97) for water and steam that is based on the Gibbs free energy [13]; this model reduces the computing time compared to the IAPWS-95 model with small accuracy sacrifices. The IF97 model is by far the most widely adopted medium model for simulation applications [17,18,19].
Largely based on the IF97 model, other researchers have developed approximation models for water and steam in order to further improve numerical performance. Åberg [14] fit fifth order polynomials to the IF97 models with spline functions to smooth the phase change transitions. Hofmann et  1 Region numbers correspond to definitions in [13]: 1 (liquid), 2 (vapor), 3 (supercritial), 4 (saturated), and 5 (high-temperature).
al. [15] implemented linear and quadratic approximation functions for pressure-based sub-regions with interpolation methods for in-between states. Affandi et al. [16] developed approximations for saturated steam functions based on fourth and fifth order power series with logarithmic transformation. In addition, there are other interpolation and table-based methods for reducing computing time and smoothing phase-change transitions. These include spline-based table look-up [20], look-up table interpolation [21], bi-cubic spline interpolation [22], biquadratic spline interpolation [23], and tabular Taylor series expansion [24,25]. However, these methods are data intensive and are often most suitable for computational fluid dynamics applications; thus they were excluded from summary Table 1.
While several steam models in Table 1 are freely available, they are not ideal for DH applications. First, the existing models cover a much larger range of p-T conditions than required for steam heating applications. Common absolute pressures for steam turbines are 0.45-1.8 MPa, and DH pressures are typically lower [26]. Steam pressures from 0.1 to 0.55 MPa are most common, but it can be as high as 2 MPa for some systems [27]. After preliminary evaluations, we selected the reduced p-T range that covers most applications while still enabling strong accuracy performance with low order polynomial fits. Second, most models span all five regions of water/steam as defined in [13], including two-phase state conditions. For steam heating, water is present as both liquid and vapor at different locations throughout the system, but equipment entraps the two-phase mixture by design. For example, steam is contained through steam traps, allowing only liquid condensate to pass through, while drip legs along distribution pipes collect and separate liquid condensate from the steam vapor. Thus, a reasonable modeling assumption for DH applications is to assume the steam is always at a vapor state (saturated or superheated), while condensate is always at a liquid state, and any two-phase conditions are contained within equipment (e.g., boilers, heat exchangers). Modeling two-phase flow within one medium model can be particularly burdensome for numerical solvers because the existing models introduce discontinuities at the phase change boundary, which can cause chattering [22], a well-known problem, and at times simulation fail-ure. Third, most of the above models contain separate forward and backward equations that are not analytical inverses, i.e., h(p, T (p, h * )) h * for some p, h * ∈ . While this can reduce computing time, it can cause serious error accumulation in lengthy calculations [24] and also be a source of chattering [28].
Others have attempted to address the numerical challenges associated with chattering to varying levels of success. Because chattering in two-phase flow systems is largely caused by large discontinuities in thermodynamic property derivatives at the phase change boundary, a popular solution is to smooth the problematic derivatives such that they are C 1 functions (aka, once continuously differentiable) [29]. This approach was taken in [18] using a bicubic spline interpolation method which eliminated the chattering problem, but resulted in poor accuracy with relative errors as high as 50%. Li et al. [28] similarly avoided the chattering problem when simulating refrigerant properties by implementing C 1 property functions by means of the splinebased table look-up method; this produced favorable results, but as previously mentioned, this is a data-intensive approach that we would like to avoid for large-scale DH system simulations. Another approach to chattering is a heuristic method that detects and minimizes the chattering problem [22]. This approach produced good accuracy. Moving boundary models also have been used to avoid chattering by identifying the phase boundary and adjusting the grid accordingly [30,31]. These have produced favorable accuracy results, but the resulting model is complex.
To address these gaps and limitations, we propose a new steam model with sufficient accuracy and numerical performance for large-scale thermo-fluid system modeling and simulation. While this model is most suitable for space and domestic water heating applications in the built environment, other applications in low and medium pressure ranges are also possible. To reiterate, the primary difference between this new model and previous steam models is the reduced p-T range, which enables lower order polynomial approximations for key thermodynamic functions. Our hypothesis with this implementation was that replacing frequently called functions with invertible polynomial approximations can produce sufficient accuracy while greatly reducing the computing time for large-scale thermo-fluid sys-tem simulations such as DH. In addition, we develop several new component, equipment, and system models that allow numerically efficient simulations of DH systems with our novel approach.
For simulations involving phase change, this improved steam medium can be coupled with a numerically-efficient liquid water model. As previously discussed, a reasonable assumption for DH systems is that steam is always at a vapor state, while condensate is always at a liquid state, and any two-phase conditions are contained within equipment. To our best knowledge, this split medium approach is atypical and has not yet been tested in the context of system simulation of steam heating and industrial processes. We hypothesize that splitting the liquid and vapor phases into two models (thus, eliminating the two-phase region) can help avoid the numerical challenges frequently experienced at the phase change boundary, such as chattering, and improve the overall numerical performance. Further, a split-medium approach allows the use of numerically efficient liquid water models, such as those in the Modelica Buildings Library [32] that improved computing performance by decoupling mass and energy balance equations. This paper demonstrates the performance of our proposed splitmedium approach.
The rest of this paper is organized as follows. In section 2, the mathematical formulations for the new polynomial approximations of thermodynamic property functions are presented as well as the Modelica implementation. In section 3, we present our approach for evaluating the water and steam implementation across several scales: at thermodynamic property scale, component scale, and district scale. The results and discussion for each of these three levels of evaluation are then presented in section 4. Lastly, final conclusions are provided in section 5.

Mathematical Formulation
To find polynomial approximations for inherently complex, natural thermo-physical phenomena, a least squares approach is one of the best and most commonly employed methods. With A ∈ R m×n and m > n, least squares is a special minimization problem that solves the overdetermined matrix problem Ax = b by minimizing the residual r(x) = b − Ax, where A is the independent data points, x is the vector of the independent variables, and b is the solution data points. In principle, any norm could be used to find an x ∈ n that minimizes r(x) ; however with least squares, the two-norm is used because (1) x 2 2 is a smooth function for any x ∈ n , and (2) the problem leads to a linear system of equations. Conversely, other norms like the ∞ or one-norm lead to nonlinear, non-differentiable optimization problems. Thus, the least squares problem is formulated as where m is the number of data points in the fit. There are many different ways to solve this least squares optimization problem.
In this work, we employ a linear least squares method for twodimensional polynomial surface fits of several orders using the CurveFit toolbox with MATLAB R2019b. More information on the particular algorithms is available in [33]. It is common for matrix problems involving steam thermodynamic properties to be ill-conditioned because of the structure of the relationship between different parameters. Conceptually, ill-conditioned means that small percentage errors in the input variables (e.g., pressure) can lead to large percentage errors in the solution (e.g., specific enthalpy). Thus, before solving the least squares problem if A is ill-conditioned, one can improve the conditioning of A to minimize the sensitivity to round off errors and alleviate numerical problems. In this work, we improve the conditioning of input variables pressure and temperature by centering and scaling the inputs to standard-normal distributions using p = p − p mean p sd and (2) where p and T are the normalized pressure and temperature, respectively; and subscripts mean and sd signify the mean and standard deviation, respectively, of the selected input data set. For interested readers, much more information on conditioning and matrix problems in the context of numerical methods can be found in [34]. The invertible polynomial approximations for specific enthalpy and entropy, presented in the following sections, were fit using the least squares approach described above. Independent variables x = p T T were selected. Data tables for saturated steam along the dew curve and superheated steam were generated over the p-T range from the IF97 formulations (Appendix A). A dense input grid was developed at standard increments of 0.1 K and 1 kPa, resulting in m = 118, 982 points. Several polynomial orders from first through fifth were evaluated before selecting the lowest order fits with acceptable accuracies for DH applications. Table 2 presents the normalization parameters used that pertain to the selected p-T input grid for the polynomial approximations.

. Property Selections
Because of the highly nonlinear structure of thermodynamic properties of steam, fitting all properties with low order polynomials was not possible. Furthermore, replacing all functions is not necessary because (1) some functions are infrequently called and thus do not have significant impact on computing time, and (2) some functions do not exhibit problematic nonlinearities at the phase change boundary and thus do not typically contribute to the numerical problems. Therefore, we selected two thermodynamic property functions to replace due to their structure and impact potential: specific enthalpy and specific entropy. We also considered replacing the density functions ρ(p, T ) and T (p, ρ); however, because invertible polynomials could not accurately represent these functions over the entire range, we decided to keep the original IF97 implementation. In the future, we plan to also evaluate the potential benefits of higher order polynomial approximations for density functions in the new steam model.

Specific Enthalpy
The current IF97 formulation for specific enthalpy is rather complex with 5 equations and 147 parameters (coefficients and exponents). In addition, the implementation has inconsistent forward and backward equations since the forward function h(p, T ) is non-invertible, and an additional backward equation with 102 unique parameters is also required. Thus, the goal here is to replace the complex IF97 equations for h(p, T ) and T (p, h) with simple polynomial approximations that can be inverted analytically, i.e., h(p, T (p, h * ) = h * for any p, h * ∈ . For reference, the IF97 Region 2 (vapor) implementations for specific enthalpy are given in Appendix A.
With our improved implementation for DH applications, the linear approximation for specific enthalpy is given as where h is specific enthalpy and a 1 through a 3 are regression coefficients. Setting the matrix problem as we then use Equation 1 and formulate the least squares problem as and solve for coefficients a = (a 1 , a 2 , a 3 ). Once the coefficients are determined, the backward equation can be formulated directly from Equation 4 by solving for T as where Table 3 lists the regression coefficients for Equation 4 found using this least squares approach.

Specific Entropy
Similar to specific enthalpy, the IF97 formulation for specific entropy is also rather complex and contains inconsistent forward and backward equations. The IF97 formulation for specific enthalpy (provided in Appendix A) contains 7 equations and 294 parameters (coefficients and exponents) for the forward function, while the backward function contains 138 unique parameters. While a linear surface fit did not produce sufficient accuracy for specific entropy functions, an invertible polynomial approximation was still possible. Using a function that is quadratic in pressure and linear in temperature, specific entropy can be calculated as where s is specific entropy and d = ( is written as products of sums rather than sums of products where possible to improve its numerical computing performance. While higher order models could produce better fits, Equation 9 was selected because its accuracy was sufficient while still meeting our goal of having invertible functions. Similar to specific enthalpy, we formulate the matrix problem and solve the least squares problem Solving for T directly from Equation 9, we write the backward function as Using this approach, we obtained regression coefficients for Equation 9, which are listed in Table 4.

Modelica Implementation
The new steam medium model is implemented in the equation-based, object-oriented modeling language Modelica. A complete steam vapor model is open-source released in the Modelica IBPSA [35] and Buildings Libraries [32]. While these open-source libraries are continually under development, Because the commonly-adopted IF97 water medium was developed for industrial applications with simulation speed in mind, it was selected as the basis for our improved steam model. Most property functions are calculated using the IF97 formulations, as implemented in the Modelica Standard Library version 2.3 Media.Water.StandardWaterOnePhase package. Because phase change was not to be included in the steam model, we used the typical single-phase implementation with p and T as independent variables.
However, in order to improve the numerical efficiency of this medium for large-scale thermo-fluid system simulations, we made some important modifications. First, we reduced the applicable p-T range appropriately for DH applications. For this implementation, the valid temperature range is 100 • C ≤ T ≤ 160 • C, and the valid pressure range is 100 kPa ≤ p ≤ 550 kPa. It is important to note that these ranges do not encompass all DH applications, but represent the majority of installations. The reduced range was selected so that invertible polynomial functions can accurately represent the commonly-called h and s functions and still meet most DH modelers' needs. While the steam model can be used outside of this p-T , the accuracy is not guaranteed and must be verified by users. Second, we replaced commonly called functions for energy and exergy analysis with first and second order polynomial approximations. These two functions are specific enthalpy h(p, T ) (Equation 4) and specific entropy s(p, T ) (Equation 9), and their corresponding backward functions Equation 7 and Equation 12, which are analytical inverses of the forward functions. Third, we implemented expressions for the derivatives to facilitate symbolic processing of function derivatives in the numerical solver.
Further, we couple the new steam medium with a numerically-efficient liquid water model, which we refer as a "split-medium" approach. Because the steam model is vapor only and the liquid water model is liquid only, the two-phase region is avoided in our implementation. Thus, numerical challenges caused by the two-phase region (i.e., chattering) are also eliminated. With this Modelica implementation, we evaluate the performance of the new steam medium model as well as our proposed split-medium approach across several scales.

Evaluation Approach
To evaluate the accuracy and numerical performance across multiple scales, Modelica models are developed at thermody-namic property scale (simple function evaluations), component scale (a mixing volume), and district scale (complete heating districts of several sizes). Several medium model configurations are included in the component and district scale evaluations, which are summarized in Figure 1. Two cases involve the IF97 model with two-phase flow -I(MBL) and I(MSL) -while the remaining three use our split-medium approach -S+CW, S+VW, and I+VW. With the two IF97 cases, models from two separate libraries are evaluated: the Modelica Buildings Library (MBL) and Modelica Standard Library (MSL). While the models from both of these libraries can be applied for steam DH applications, their original design intentions differ, and correspondingly, they produce different simulation results. The I(MBL) case is only used in the component-scale evaluation, while the I+VW case is only used in the district-scale evaluation. The remaining three cases are common to both. Further details regarding the model setup for each simulation case are provided in the following sections.  Several new component and system models for DH modeling were developed for these case studies. While presenting the complete details of these models is out of the scope of this paper, a general description of each model's physics principles is provided herein. These models are in the process of being refined and open-source released in the Modelica Buildings Library. Beyond newly developed models, all components used are existing in the Modelica Standard Library v3.2 and Modelica Buildings Library v7.0.0.

Thermodynamic Property Scale
The accuracy and numerical performance of the proposed specific enthalpy and specific entropy approximations in Equation 4 and Equation 9 are evaluated with respect to the IF97 model. The objective of this evaluation was to evaluate solely the thermodynamic properties separate from higher-level effects from components, equipment, and systems. For accuracy, we calculated the residuals between the new approximation and the IF97 model across the entire p-T range. In addition, we evaluated the Root Mean Square Error (RMSE) and Coefficient of Variation of the Root Mean Square Error (CVRMSE) using where y i is the individual reference data generated by the IF97 model,ŷ i is the corresponding evaluation data predicted by the new model,ȳ is the mean of the reference dataset, and N is the total number of data points. For numerical performance, h and s were calculated over a dense grid of p-T states. The same simulation settings and computer were used across all evaluations in this paper (details in subsection 3.4). The range of p-T states corresponded to the design limits of the new steam model. Modelica code for the specific enthalpy evaluation is given in Listing 1 as an example. With a fixed p state, properties are calculated across the range of T , evaluating state properties in the forward (state_pTX) and backward (state_phX) directions to check for consistency. Assert statements notify if the error in any of the property evaluation exceeds the predetermined tolerance (set to 10 −8 for the new steam model). This model is simulated across the range of p, sweeping from 100 kPa to 550 kPa at increments of 250 Pa, for a total of 1,801 simulations. For comparison with a reference model, the same script is repeated with the IF97 model by setting the Medium declaration to Modelica.Media.Water.StandardWaterOnePhase. Because in the IF97 model, the backward functions are not analytical inverses of the original functions, the simulation fails with an error tolerance of 10 −8 and needed to be increased to 10 −1 for completion. Note, these error tolerance values are not simulation tolerances, but absolute error thresholds (parameter err in Listing 1).

Component Scale
Because of the discontinuities at the phase-change barrier, numerical solvers have shown to experience difficulties solving thermo-fluid problems with two-phase flow. Chattering is one example of a well-known issue that has been demonstrated previously with a Modelica-based simulation of a boiler pipe model featuring the IF97 medium [22]. For this case, we use water and steam mixing volumes to evaluate the performance of the new implementation compared to existing models using IF97 and the standard mixing volume model in the Modelica Standard Library and Modelica Buildings Library. This component-scale experiment was selected to isolate some of the common numerical challenges of thermo-fluid system modeling involving phase change that may not appear at smaller scales while being harder to diagnose at larger scales.
Shown in Figure 2, the Modelica-based evaluation features a water mixing volume of 0.1 m 3 that is exposed to a constant temperature boundary of 300 • C via a thermal conductor with a constant thermal conductance of 10 W/K. The reference pressure is set to 200 kPa. The mixing volume is replaceable in order to allow four different simulation runs to be completed with the same boundary conditions. Four cases from Figure 1    work. Details regarding this model's design principles and its mathematical formulation are included in Appendix B.

District Scale
The objective of the district scale evaluation was to assess the accuracy and numerical performance of the new implementation in a typical DH system design across districts of several sizes. An overarching description of the selected system is presented next, followed by the Modelica implementation. Figure 3 depicts the schematic diagram for the DH evaluation. This DH system is broken into three subsystems: a central plant, the distribution network, and building end users. The central plant features a feedwater tank, feedwater pump, and a single steam boiler. Saturated steam is discharged from the boiler at 300 kPa. The fuel load ratio for the boiler is controlled to maintain the boiler discharge pressure, while the feedwater pump speed is controlled to maintain the liquid water level in the boiler. While there are several mechanical and control designs seen in central plants [37], this configuration was selected because it represents real-world control dynamics while representing one of the more simple designs. For the distribution network, we assume there are no mass nor energy losses in the steam supply pipe, while the condensate return pipes have fixed pressure drops without any heat transfer. Lastly, the buildings contain a steam heat exchanger, a steam trap, and a condensate return pump. The building-side components are assumed to have steady energy and mass balances, with the mass flow and heat transfer determined based on tabulated load data.

System Description
For this DH case study, we adopted tabulated heating load profiles from the common exercise of IBPSA Project 1 [35] WP3.1, one of the work packages of the international collaboration project that focuses on district energy system simulations. This building heating load features a single-family residential building in the heating dominated climate of Belgium. More information regarding development and testing of the building model can be found in Saelens et al. [38].

Modelica Implementation
To evaluate the steam medium implementation across a variety of district sizes, a vector-style DH system model was developed in Modelica (Figure 4). This top-level Modelica diagram has a clear one-to-one relationship with the subsystems depicted in Figure 3. A parameter N representing the total number of buildings can be adjusted to represent districts of multiple sizes. The tabulated heating load profile for each building is scaled by 1/N in order to keep the central plant sizing constant across all simulation cases. This DH system is simulated for the first 15 days of the year.
The plant model diagram is shown in Figure 5. In the central plant, the feedwater pump and boiler both have dynamic energy and mass balances, and PI controllers are used to maintain the water level and pressure setpoints. This plant model was designed to have a consistent interface for several district energy plant varieties, including chilled and hot water production. A check valve was also included in the central plant model in order to prevent unintended reverse flow. The boiler model contains the MixingVolumeEvaporation to represent to fluid volume. The boiler also includes pressure drop in the working fluid and represents the heat transfer with the ambient environment and with the fuel through heat ports. Both the thermal conductance and heat capacity of the boiler drum metal and insulation are included in the model. The building model diagram is shown in Figure 6. Through the Tabulated Heating Load data reader, the variable heat flow rateQ from IBPSA Project 1 WP3.1 is input directly. A condensate return pump is modeled with a prescribed mass flow rate ofṁ =Q/∆h f g,std , where h f g,std is the standard enthalpy of vaporization at the design pressure. The Heat Exchanger Volume model is configured similarly to the MixingVolumeEvaporation model, except it represents a condensation process instead, and the steady state balance is adopted for this instance. Because the heat flow rate is directly proportional to the mass flow rate (given that the change in enthalpy is essentially fixed), a heat port is not permissible with the steady balance. Lastly, the Steam Trap represents an isenthalpic process where liquid condensate is discharged at atmospheric pressure, while any flashed steam vapor due to the drop in pressure is returned to a liquid state before discharge.
Similarly to subsection 3.2, four medium model configurations were included in this evaluation. From Figure 1, the four cases are: I(MSL), S+CW, S+VW, and I+VW. Relative to the component-scale evaluation, an additional split-medium case (I+VW) was included that features the IF97 steam model Stan-dardWaterOnePhase from the Modelica Standard Library and the variable density water model TemperatureDependentDensity from the Modelica Buildings Library.

Simulation Settings
All simulations ran in Dymola 2021 on a Windows 10 workstation with a Intel ® Xeon ® 3.60GHz CPU and 32.0GB of RAM. The DASSL solver was selected for all case studies, after preliminary tests demonstrated its lower computing time compared to other numerical solvers. The simulation tolerance was set to 10 −6 for all cases as well.

Evaluation Results and Discussion
This section presents the results from simulations across the three scales: thermodynamic properties, component, and dis-   trict scales. First the accuracy of the model will be presented, followed by the numerical performance.

Thermodynamic Property Scale
Following the methodology prescribed in subsection 3.1, the results for model accuracy and numerical performance at the thermodynamic property scale evaluation are as follows.

Accuracy
Property evaluations for specific enthalpy and entropy produced acceptable accuracy across the entire p-T range. As shown in Figure 7, the largest residuals for h(p, T ) occurred near the saturation curve and at low p-T conditions. The reason for this is that the nonlinearities in h(p, T ) are higher in these regions than other regions. The residuals ranged from −2.1 kJ/kg (0.075%) to 2.1 kJ/kg (0.079%), which correspond to the low p-T conditions and the middle of the saturation curve, respectively. Across the entire range, the RMSE was 0.76 kJ/kg and     CVRMSE was less than 0.0001%. Fitting the linear approximation using a least square approach, the R 2 value was 0.999. Overall, a linear surface fit is able to capture h(p, T ) well in this p-T range.  As shown in Figure 8, the largest residuals for s(p, T ) occurred at high pressure-temperature conditions. The shape of the residuals indicate that the errors are relatively consistent with respect to temperature, but the quadratic fit in pressure underestimates s slightly in the middle of the range while it overestimates s at the high and low limits. This shape is due to the fact that a low order polynomial is being fit to a nonlinear function, but it still is able to capture the general functional form. The residuals ranged from −0.058 kJ/kg-K (0.82%) to 0.13 kJ/kg-K (1.97%). Across the entire range, the RMSE was 0.013 kJ/kg-K and CVRMSE was < 0.0001%. Fitting the linear approximation using a least square approach, the R 2 value was 0.996. It is important to note that the distribution of p-T grid points selected for the least squares optimization influences the shape of these fits. Based on the relatively high accuracy achieved across the entire range of conditions, we can conclude that the selected grid was appropriate.

Numerical Performance
The results in Table 5 show that the new polynomial implementations reduce the computing time by 10%-39%. These time savings are average values from 10 simulation runs. In the problem formulation, both the number of variables and equations were lower in the new models than the IF97 models. The number of variables decreased from 91 to 29 (a 68% reduction) for both properties, and the number of equations decreased from 28 to 23 (18%) for enthalpy and 27 to 22 (19%) for entropy. In addition to many other factors, these can partially explain the computing time reductions. While the simplified implementations for h and s produced numerical savings at the thermodynamic property scale, the savings are often compounded for larger thermo-fluid system models that involve nonlinear systems of equations. This will be evaluated with the district scale models. However, the computing time savings and calculation accuracies demonstrate how fundamental improvements in steam property modeling can be achieved through simple function replacement.

Component Scale
Following the methodology prescribed in subsection 3.2, the results for model accuracy and numerical performance at the component scale evaluation are as follows. Figure 9 depicts the temperature and mass results of the mixing volume evaluation cases. In this simulation, the mixing volume fluid is heated via a constant, high temperature boundary. With the I(MBL) case, undesirable erratic behavior is clear. At initialization, the liquid in the volume is at 40 • C. For the first 28 minutes, the temperature gradually increases. However, when water starts to boil, the fluid temperature oscillates between the saturation temperate and the boundary temperature of 300 • C. Mass and density are also oscillating in the I(MBL) case, as the fluid switches back and forth between one-phase (liquid or vapor) and two-phase states (liquid-vapor mixture). The mixing volume in the MBL case was designed for single phase fluid (air, water, etc.), and functions correctly in those instances. However, its performance with a two-phase IF97 medium is undesirable with the selected numerical solver of DASSL. It is important to note that the MBL behavior seen here is a numerical problem, not a physics problem. Indeed, this numerical problem can be avoided if an explicit, fixedtime step method is employed (i.e., Euler) rather than an implicit, variable-time step method (i.e., DASSL). With Euler, the I(MBL) case performs physically correct and avoids the numerical chattering issue. Further, the I(MBL) mixing volume can simulate not only a saturated fluid, but supercooled and superheated fluids as well, which may be of interest for some use cases. However, implicit, variable-time step methods are often used for large thermo-fluid system simulations because of their ability to deal with stiff systems. Thus, the Euler results were excluded from this paper. The later three cases -I(MSL), S+CW, and S+VW -perform with the correct physics for most of the simulation. Because the MSL mixing volume and new MixingVolumeEvaporation were designed for boiling processes with the fluid at a saturated state, each of these three cases initialize at the saturation temperature. Throughout the boiling process, the liquid volume gradually decreases as more of the water is converted from liquid into vapor. In the I(MSL) case, the liquid water continues into a "negative" mass state after all of the liquid has boiled off. This is physically not correct. Conversely, the new mixing volume implementation is designed to stop the simulation when no more liquid water is present because it is configured to always be at a saturated state. Therefore, supercooled liquid nor superheated vapor are not permissible to be within the Mix-ingVolumeEvaporation as is currently designed. Note that this assumption is valid for steam boilers during normal operation in most heating applications [36].

Accuracy
With the new numerically efficient implementation, some minor accuracy errors are introduced. Since the steam implementation is intended for normal operating conditions, the RMSE and CVRMSE are presented in Table 6 for mixing volumes with liquid-to-vapor ratios within the ratio of 1 to 9 which are common in real operations. Compared to the I(MSL) case -an IF97 steam medium with the MSL equilibrium drum boiler model -the S+VW case generally produces lower errors than then S+CW case. This was to be expected, because the liquid water with constant density is being applied outside of its design temperature range, causing higher than normal errors to be introduced. However, with that said the S+CW case did produce less than 1% errors in terms of CVRMSE for most property evaluations. The largest errors were introduced in the S+CW case for ρ w and m, with CVRMSE values of 5.9% and 7.6% respectively. However, these errors were greatly reduced by changing the liquid water model from a constant density (S+CW case) to a variable temperature-dependent-density model (S+VW case). However, because boiling a volume of water completely is not a typical normal-operation scenario for DH applications, both the S+CW and S+VW models can be deemed acceptable. When evaluating accuracy with the mixing volumes, we also discovered that errors in internal energy and mass calculations increase as mass goes to zero. This is an atypical scenario in normal operation, because boilers are controlled to maintain a water level setpoint. However, it is important to recognize modeling sensitivity at low mass quantities and mass flow rates, particularly with large, complex, DH system models where problems can be difficult to diagnose. One solution to this problem can be to include smoothing functions that account for regularization at low mass conditions to avoid diverging total property evaluations such as U. This will be implemented in future versions of the steam DH components in the Modelica Buildings Library.

Numerical Performance
With this small simulation case, there were negligible differences in computing time between the four cases. However, the structure of the differential algebraic system of equations indicates the likelihood of numerical performance differences for larger thermal-fluid system models. While each of the four cases contain one continuous time state variable, the nonlinear systems of equations varies. The I(MBL) case contains one nonlinear system involving two iteration variables: du/dt and dρ/dt. The I(MSL) case also contains a nonlinear system but with one iteration variable: dV w /dt. Conversely, the new implementations S+CW and S+VW contain no nonlinear systems. While the nonlinear systems are small for this component-scale evaluation, their sizes and quantities will grow as DH system models increase in the number of buildings. From experiences, the time required to iteratively solve the nonlinear systems largely contribute to the total computing time. Thus, this likely will have an impact on computing time for DH applications, which will be quantitatively tested in the next section.

District Scale
Following the methodology prescribed in subsection 3.3, the results for model accuracy and numerical performance at the district scale evaluation are as follows. To evaluate impacts with respect to the number of buildings N, simulations are repeated with 1 to 10 buildings. Figure 10 shows that the measured heat flow rate at the building for each of the four cases followed the input data file. Given that the input to the building model was the data table and the condensate return pumps controlled the mass flow rate ideally, this was generally expected. The error differences in the S+CW and S+VW cases are due to the building model implementation. Because the mass flow rate is prescribed based oṅ m =Q/∆h f g,std , errors in the h f g,std calculation affect the mass flow rate. This error could be resolved by changing the condensate return pump control to use a PI controller that adjusts the measured heat flow rate to track the tabulated input value. As seen in Table 7 and Table 8, the accuracy of thermodynamic property evaluations had insignificant differences across all cases. Here, the building 1 results are used as proxies for all simulation cases, since the accuracy did not significantly change building to building, with respect to N, and at different points throughout the system. Subscripts in and out refer to the inlet and outlet ports of the building model, where the inlet state is steam vapor and the outlet state is liquid water condensate. As expected, the errors generally increased inversely to the steam/water model complexity. In this experiment, the I+VW case produced the lowest errors compared to the I(MSL) case, followed by the S+VW case. While the S+CW case produced some of the higher errors, all errors were within acceptable ranges. For example, the CVRMSE values are within 1% for all property evaluations. With this split medium approach and simple polynomial approximations for h and s, these DH simulations were able to produce acceptable accuracy for all thermodynamic property states.  Throughout the 15 day simulation, the difference in total fuel consumption by the boiler was less than 2% for all cases. However, interestingly, slight differences in thermodynamic states caused short term variance (e.g., on the order of 1 to 5 hours) in boiler fuel flow rate (see Figure 11). With this system configuration, the boiler fuel flow rate is controlled to maintain the pressure at the boiler's outlet. This is a common design for steam DH system. As seen in Figure 11(a), the measured relative pressure at the boiler outlet does not follow the same exact trend line for all four simulation cases. In this figure, relative pressure p rel is calculated as

Accuracy
where p mea is the measured pressure at the boiler outlet and p set is the pressure setpoint of 300 kPa. Note that p rel in this case is controlled within ±0.005%. This minuscule variation in the measured input signal for the PI controller caused the boilers' transient fuel burn rates to vary across the four cases. Based on these results, we can conclude that the fuel consumption rate is sensitive to the thermodynamic property evaluation, particularly with realistic variable-load control. This DH evaluation revealed the importance of well-tuned feedwater pump and boiler controls on power and fuel flow rates. One major benefit of Modelica is the ability to simulate realistic control of complex thermo-fluid system models. While these controlability aspects did not greatly influence the overall energy consumption in this case study, it did heavily impact the transient fuel consumption rate. For simulation scenarios that deal with minute to daily time scales in particular, a major takeaway is that the fuel rate greatly depends on the model's ability to calculate thermodynamic properties accurately. If this is the case, more accurate steam and water models such as in the I+VW case would be favorable. Table 9 presents the structure of the translated model for each of the four cases. The number of equations varied slightly between the I(MSL) case, the variable water I+VW and S+VW cases, and the constant water S+CW case. Notably, the split medium approaches reduced the quantity of linear and nonlinear systems present in the district models. For the I(MSL) case, both the number of linear systems and the dimension of the first nonlinear increase with the number of buildings. Conversely, the linear and nonlinear systems do not change with respect to N for the three split medium cases. This is advantageous particularly for simulations of large DH systems.

Numerical Performance
As seen in Figure 12, computing times were significantly re- duced with the split medium approach and with the improved steam medium model. With the new steam model, the S+CW case was on average 8.5 times faster than the I(MSL) case, with the largest speed up of 9.3 times occurred with N = 10. The S+VW case and I+VW case were on average 6.5 and 4.5 times faster than the I(MSL) case as well. As the number of buildings and number of continuous time states increased, the computing time increased sub-linearly from 1 to 10 buildings. This is equal to or better than some previous district energy simulation studies involving small districts, where linear and quadratic scaling was present [6,39]. However most importantly, the rate of computing time change begins to increase in the I(MSL) case with 9 and 10 buildings. This along with the results found in Table 9 indicate that the computing savings for larger district may be orders-of-magnitude greater than the quantified computing time savings in this study. We are able to achieve this result with the new steam model, the liquid water models in the Modelica Buildings Library, and our new component models because in them, density is decoupled from pressure. This allows the Modelica tools to decouple the energy and mass balance equations and eliminate nonlinear systems that are computationally costly. In the future, we will expand upon this work to evaluate large DH systems featuring hundreds of buildings.

Conclusion
This work introduces a new steam medium model that enables numerically efficient DH simulations. Accuracy and numerical performance were aptly balanced by reducing the ap- The new steam model is open-source released in the Modelica IBPSA and Buildings Libraries, while new component, equipment, and system models for steam-based DH systems are in the process of being publicly released. With district-scale simulations growing in importance as societies strive to meet energy and climate targets, this work enables the modeling and simulation of complete steam-based DH systems that, for all intents and purposes, was not previously possible. Depending on the target application, several options are available for steam and water medium models that can strike favorable balances between numerical performance and accuracy. For example, we suggest coupling the new steam model with a constant density water model for large-scale energy simulations, while a variable temperature-dependent-density water model is more suitable for hourly optimization of steam plants and districts. Regardless of the specific application, this new steam model can help aid the transition of legacy DH systems to newer sustainable designs, while providing a pathway for district-scale anal-ysis and optimization in the many steam heating applications that are likely to remain.

Acknowledgements
This research was supported in part by an appointment with IBUILD Graduate Student Research Program sponsored by the U.S. Department of Energy (DOE), Office of Energy Efficiency and Renewable Energy, and Building Technologies Office. This program is managed by Oak Ridge National Laboratory (ORNL). This program is administered by the Oak Ridge Institute for Science and Education (ORISE) for the DOE. ORISE is managed by ORAU under DOE contract number DESC0014664. All opinions expressed in this paper are the author's and do not necessarily reflect the policies and views of DOE, ORNL, ORAU, or ORISE. In addition, this material is based upon work supported by the DOE's Office of Energy Efficiency and Renewable Energy under the Advanced Manufacturing Office, award number DE-EE0009139, and the Building Technologies Office, contract number DE-AC02-05CH11231. Further, this work emerged from the IBPSA Project 1, an international project conducted under the umbrella of the International Building Performance Simulation Association (IBPSA). Project 1 will develop and demonstrate a BIM/GIS and Modelica Framework for building and community energy system design and operation.

Disclaimer
This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.

Appendix A. IF97 Region 2 Equations
Complete details for the IF97 formulation of water and steam are available in [13,40]. However, to provide a reference point for the complex functions we replaced for DH applications, relevant formulations for specific enthalpy and entropy are reproduced here. In the IF97 formulation for Region 2 (vapor), basic equations for Gibbs free energy g are expressed in a dimensionless form γ(π, τ) = g(p, T ) RT = γ 0 (π, τ) + γ r (π, τ), (A.1) where γ 0 represents the ideal-gas part and γ r represents the residual part. This equation is formulated in terms of the specific gas constant of ordinary water R = 0.461526 kJ/(kg · K), the reduced pressure π = p/p * , and the inverse reduced temperature τ = T * /T , where superscript * indicates the reducing quantity. For these equations, p * = 1 MPa and T * = 540 K. The form of the ideal-gas part γ 0 is where each of the 9 coefficients for n 0 i and J 0 i are tabulated in Table 10 of [40]. Similarly, the form of the residual part γ r is given as where each of the 43 coefficients for n i and exponents I i and J i are tabulated in Table 11 of [40]. In terms of Gibbs free energy, specific enthalpy h = g − T (∂g/∂T ) p . Thus, the dimensionless relation in terms of the ideal-gas and residual parts can be expressed as h(π, τ) RT = τ ∂γ 0 ∂τ π + ∂γ r ∂τ π , (A. 4) where the partial derivatives of γ 0 and γ r with respect to τ at constant π are ∂γ 0 ∂τ π = In terms of Gibbs free energy, specific entropy s = −(∂g/∂T ) p , and the IF97 relation is s(π, τ) R = τ ∂γ 0 ∂τ π + ∂γ r ∂τ π − γ 0 + γ r . (A.7) This completes the forward equations for h(p, T ) and s(p, T ) in the IF97 formulation, Region 2. In total, h(p, T ) includes 5 equations, 52 coefficients, and 95 exponents, while s(p, T ) includes 7 equations, 104 coefficients, and 190 exponents. In addition to the forward equations, separate backward equations are similarly implemented. However, Region 2 is further divided into three subregions 2a, 2b, and 2c (see [40] for subregion divisions). Thus, three T (p, h) and three T (p, s) cover Region 2, but only subregion 2a is applicable for the reduced p-T range specified for DH applications in Table 1. The dimensionless backward equation T (p, h) for subregion 2a is where η = h/h * , h * = 2000 kJ/kg, and each of the 34 coefficients for n i and exponents I i and J i are tabulated in Table 20 of [40]. Lastly, the backward function T (p, s) for subregion 2a is where σ = s/s * , s * = 2 kJ/(kg · k), and each of the 46 coefficients for n i and exponents I i and J i are tabulated in Table 25 of [40]. In total, T (p, h) includes 34 coefficients and 68 exponents, while T (p, s) includes 46 coefficients and 92 exponents. Modelica implementations of the above forward and backward equations for specific enthalpy and entropy along with the other thermodynamic properties are used for comparative evaluation in this study.

Appendix B. Mixing Volume Evaporation Model
The new MixingVolumeEvaporation model, which will be open-source released in the Modelica Buildings Library, represents the evaporation process of water with liquid and vapor phases in equilibrium. It is designed for the fluid to be at a saturated state only, which is a reasonable assumption for steam boilers operating under normal conditions. The fundamental equations are as follows. The fluid mass m in the volume is calculated as m = ρ s V s + ρ w V w , (B.1) where ρ is density, V is volume, and subscripts s and w represent the steam and liquid water components. The total internal energy U is where h is specific enthalpy, p is pressure, and V = V s + V w is the total volume of the fluid. The MixingVolumeEvaporation is configured to allow both steady state and dynamic mass and energy balances. The dynamic mass balance is given as whereṁ s andṁ w are the mass flow rates of steam and liquid water, respectively. Note that with the split medium approach, the pure liquid water is always assigned at the upstream port while the pure steam vapor is always at the downstream port. The dynamic energy balance is given as dU dt =Q +ṁ s h s +ṁ w h w , (B.4) whereQ is the net heat flow rate through the boiler's enclosure and from the fuel. Two principal assumptions are made with this model. First, we assume that the water is always at a saturated state within the boiler, and saturated steam vapor with a quality of 1 is discharged from the outlet port with h s = h v . Second, we assume that the liquid and vapor components in the volume are at equilibrium. Further, it is important to note that the new MixingVolumeEvaporation is very similar to Modelica.Fluid.Examples.DrumBoiler. Base-Classes.EquilibriumDrumBoiler [36] with the following exceptions: 1. Rather than a two-phase medium, fluid mediums are modeled as two single-state fluids, with liquid water at the upstream port, and steam vapor at the downstream port. 2. The metal drum and insulation are excluded from the mass and energy balances of the mixing volume. Instead, these are included in the boiler model, which instantiates the mixing volume. Note in the component-scale experiment, the mass of the drum metal was set to zero in the MSL case for consistency purposes. 3. The new model protects against incorrect physics where the steam and liquid water volumes must always be equal to or greater than zero. 4. Lastly, the steady state balances accurately hold mass and internal energy constant.