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. Keywords: Steam, District Heating, Numerical E ciency, Modeling, Simulation, Modelica 1. Introduction District heating (DH) can e ectively reduce CO 2emissions and enable communities to leverage economizes of scale ben- eﬁts [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, rep- resenting 97% of all installations [3]. Beyond heating build- ings alone, steam DH provides beneﬁcial waste-heat recovery opportunities when coupled with power plants (e.g., combined heat and power) or other industrial systems (e.g., wastewater treatment, metal reﬁneries, etc.). While the conversion of steam DH systems (ﬁrst generation) to hot water DH systems (sec- ond 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 mutu- alistic beneﬁts with power generation and high-heat industrial processes such as those listed in [4]. Modeling and simulation can improve the design and oper- ation 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 abun- dant [5, 6, 7, 8] but out of the scope of this paper. Further, simulation studies of steam heating systems for combined heat Corresponding author. Email address: wangda.zuo@colorado.edu (Wangda Zuo)and power plants [9, 10] and other industrial applications [11] are prevalent within DH contexts; however, to our best knowl- edge, past studies have yet to model the thermo-ﬂuid dynamics of steam heating systems at district scale. District-scale simula- tion capabilities are critical for future grid-interactive e cient buildings and for our ability to reach carbon targets through ﬁ- nancially viable means. To accomplish this task, the ﬁrst step is to e ciently evaluate fundamental steam properties for fast and accurate DH simulation. Summarized in Table 1, several steam medium models are used in scientiﬁc 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 pand temperature Tand is intended for scientiﬁc ap- plications where accuracy is of utmost importance [12]. The IAPWS also released their industrial formulation (IF97) for wa- ter 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 sacriﬁces. The IF97 model is by far the most widely adopted medium model for simulation applica- tions [17, 18, 19]. Largely based on the IF97 model, other researchers have developed approximation models for water and steam in or- der to further improve numerical performance. Åberg [14] ﬁt ﬁfth order polynomials to the IF97 models with spline func- tions to smooth the phase change transitions. Hofmann et Preprint submitted to Energy October 23, 2021 Table 1: Review of existing steam medium models. NameApplicable RangeRegions1 Intended FreelyRef. p (MPa) T (C) Application Available IAPWS-95 0–1000 22–1000 1–5 Scientiﬁc Yes [12] IF97 0–100 0–800 1–5 All Industrial Yes [13] Polynomial Fits with Splines0.08–100 0–800 1–5 All Industrial No [14] Polynomial Fits 0.01–100 0–800 1–5 All Industrial No [15] Power Series with Log Transformations0–21 0–370 4 All Industrial Yes [16] Desired Steam Model 0.1–0.55 0–160 2 Heating Only Yes – 1Region numbers correspond to deﬁnitions in [13]: 1 (liquid), 2 (vapor), 3 (supercritial), 4 (saturated), and 5 (high-temperature). al. [15] implemented linear and quadratic approximation func- tions for pressure-based sub-regions with interpolation methods for in-between states. A andi et al. [16] developed approxima- tions for saturated steam functions based on fourth and ﬁfth order power series with logarithmic transformation. In addi- tion, there are other interpolation and table-based methods for reducing computing time and smoothing phase-change transi- tions. These include spline-based table look-up [20], look-up table interpolation [21], bi-cubic spline interpolation [22], bi- quadratic spline interpolation [23], and tabular Taylor series ex- pansion [24, 25]. However, these methods are data intensive and are often most suitable for computational ﬂuid 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 mod- els cover a much larger range of p-Tconditions than required for steam heating applications. Common absolute pressures for steam turbines are 0.45–1.8 MPa, and DH pressures are typi- cally 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 sys- tems [27]. After preliminary evaluations, we selected the re- duced p-Trange that covers most applications while still en- abling strong accuracy performance with low order polynomial ﬁts. Second, most models span all ﬁve regions of water /steam as deﬁned in [13], including two-phase state conditions. For steam heating, water is present as both liquid and vapor at dif- ferent locations throughout the system, but equipment entraps the two-phase mixture by design. For example, steam is con- tained 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 su- perheated), 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 ﬂow within one medium model can be particularly burdensome for numer- ical solvers because the existing models introduce discontinu- ities at the phase change boundary, which can cause chatter- ing [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2 < . While this can re- duce 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. Be- cause chattering in two-phase ﬂow 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 C1functions (aka, once continuously di erentiable) [29]. This approach was taken in [18] using a bicubic spline interpolation method which elimi- nated 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 C1property functions by means of the spline- based 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 de- tects 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 bound- ary 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 su cient accuracy and numerical perfor- mance for large-scale thermo-ﬂuid system modeling and simu- lation. While this model is most suitable for space and domestic water heating applications in the built environment, other ap- plications in low and medium pressure ranges are also possible. To reiterate, the primary di erence between this new model and previous steam models is the reduced p-Trange, which enables lower order polynomial approximations for key thermodynamic functions. Our hypothesis with this implementation was that re- placing frequently called functions with invertible polynomial approximations can produce su cient accuracy while greatly reducing the computing time for large-scale thermo-ﬂuid sys- 2 tem simulations such as DH. In addition, we develop several new component, equipment, and system models that allow nu- merically e cient simulations of DH systems with our novel approach. For simulations involving phase change, this improved steam medium can be coupled with a numerically-e cient liquid wa- ter 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 con- ditions are contained within equipment. To our best knowl- edge, this split medium approach is atypical and has not yet been tested in the context of system simulation of steam heat- ing and industrial processes. We hypothesize that splitting the liquid and vapor phases into two models (thus, eliminat- ing the two-phase region) can help avoid the numerical chal- lenges frequently experienced at the phase change boundary, such as chattering, and improve the overall numerical perfor- mance. Further, a split-medium approach allows the use of numerically e cient liquid water models, such as those in the Modelica Buildings Library [32] that improved computing per- formance by decoupling mass and energy balance equations. This paper demonstrates the performance of our proposed split- medium approach. The rest of this paper is organized as follows. In section 2, the mathematical formulations for the new polynomial approx- imations 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 implementa- tion 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, ﬁnal conclusions are provided in section 5. 2. Mathematical Formulation and Implementation 2.1. Mathematical Formulation To ﬁnd polynomial approximations for inherently complex, natural thermo-physical phenomena, a least squares approach is one of the best and most commonly employed methods. With A2Rmnandm>n, least squares is a special minimization problem that solves the overdetermined matrix problem Ax=b by minimizing the residual r(x)=bAx, where Ais the in- dependent data points, xis the vector of the independent vari- ables, and bis the solution data points. In principle, any norm could be used to ﬁnd an x2<nthat minimizeskr(x)k; however with least squares, the two-norm is used because (1) kxk2 2is a smooth function for any x2<n, and (2) the problem leads to a linear system of equations. Conversely, other norms like the 1 or one-norm lead to nonlinear, non-di erentiable optimization problems. Thus, the least squares problem is formulated as min x2<nkr(x)k2 2=min x2<nmX i=1r2 i (1) where mis the number of data points in the ﬁt. There are many di erent ways to solve this least squares optimization problem.In this work, we employ a linear least squares method for two- dimensional polynomial surface ﬁts 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 thermody- namic properties to be ill-conditioned because of the structure of the relationship between di erent 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., speciﬁc enthalpy). Thus, before solving the least squares problem if Ais ill-conditioned, one can improve the conditioning of Ato minimize the sensitivity to round o errors and alleviate numerical problems. In this work, we im- prove the conditioning of input variables pressure and tempera- ture by centering and scaling the inputs to standard-normal dis- tributions using p=ppmean psdand (2) T=TTmean Tsd; (3) where pandTare the normalized pressure and temperature, respectively; and subscripts mean andsdsignify 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 speciﬁc en- thalpy and entropy, presented in the following sections, were ﬁt using the least squares approach described above. Independent variables x=h pTiTwere selected. Data tables for saturated steam along the dew curve and superheated steam were gener- ated over the p-Trange 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 ﬁrst through ﬁfth were evaluated be- fore selecting the lowest order ﬁts with acceptable accuracies for DH applications. Table 2 presents the normalization pa- rameters used that pertain to the selected p-Tinput grid for the polynomial approximations. Table 2: Normalization parameters for polynomial approximations. Parameter Value Units pmean 1.235e +06 Pa psd 8.326e +05 Pa Tmean 6.776e +02 K Tsd 1.611e +02 K 2.1.1. Property Selections Because of the highly nonlinear structure of thermodynamic properties of steam, ﬁtting all properties with low order poly- nomials was not possible. Furthermore, replacing all functions is not necessary because (1) some functions are infrequently 3 called and thus do not have signiﬁcant impact on computing time, and (2) some functions do not exhibit problematic nonlin- earities 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: speciﬁc enthalpy and speciﬁc 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 beneﬁts of higher order polynomial approximations for density functions in the new steam model. 2.1.2. Speciﬁc Enthalpy The current IF97 formulation for speciﬁc enthalpy is rather complex with 5 equations and 147 parameters (coe cients 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 in- verted analytically, i.e., h(p;T(p;h)=hfor any p;h2<. For reference, the IF97 Region 2 (vapor) implementations for speciﬁc enthalpy are given in Appendix A. With our improved implementation for DH applications, the linear approximation for speciﬁc enthalpy is given as h(p;T)=a1+a2p+a3T; (4) where his speciﬁc enthalpy and a1through a3are regression coecients. Setting the matrix problem as h Abi =266666666666666664p1T1h1 p2T2h2 ::::::::: pmTmhm377777777777777775; (5) we then use Equation 1 and formulate the least squares problem as min amX i=1(hi(a1+a2pi+a3Ti))2(6) and solve for coe cients a=(a1;a2;a3). Once the coe cients are determined, the backward equation can be formulated di- rectly from Equation 4 by solving for Tas T(p;h)=c1+c2p+c3h; (7) where c=[c1;c2;c3]="a1Tsd a3Tmean;a2Tsd a3;Tsd a3# : (8) Table 3 lists the regression coe cients for Equation 4 found using this least squares approach.Table 3: Regression coe cients for speciﬁc enthalpy h(p;T) and T(p;h). Parameter Value Units a1 3.272e +06 kJ /kg a2 -1.838e +04 kJ /kg-Pa a3 3.504e +05 kJ /kg-K 2.1.3. Speciﬁc Entropy Similar to speciﬁc enthalpy, the IF97 formulation for speciﬁc entropy is also rather complex and contains inconsistent for- ward and backward equations. The IF97 formulation for spe- ciﬁc enthalpy (provided in Appendix A) contains 7 equations and 294 parameters (coe cients and exponents) for the forward function, while the backward function contains 138 unique pa- rameters. While a linear surface ﬁt did not produce su cient accuracy for speciﬁc entropy functions, an invertible polyno- mial approximation was still possible. Using a function that is quadratic in pressure and linear in temperature, speciﬁc entropy can be calculated as s(p;T)=d1+d2p+d3T+p(d4p+d5T); (9) where sis speciﬁc entropy and d=(d1;d2;d3;d4;d5) are regres- sion coe cients. Note that Equation 9 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 ﬁts, Equation 9 was selected because its accuracy was su cient while still meeting our goal of having invertible functions. Similar to speciﬁc enthalpy, we formulate the matrix problem h Abi =266666666666666664p1T1s1 p2T2s2 ::::::::: pmTmsm377777777777777775; (10) and solve the least squares problem min dmX i=1(si(d1+d2pi+d3Ti+pi(d4pi+d5Ti)))2:(11) Solving for Tdirectly from Equation 9, we write the back- ward function as T(p;s)=(sd1p(d2+d4p))Tsd d3+d5p+Tmean: (12) Using this approach, we obtained regression coe cients for Equation 9, which are listed in Table 4. 2.2. 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, 4 Table 4: Regression coe cients for speciﬁc entropy s(p;T) and T(p;s). Parameter Value Units d1 7530 kJ /kg-K d2 -636.9 kJ /kg-K-Pa d3 532.8 kJ /kg-K2 d4 159.8 kJ /kg-K-Pa2 d5 4.130 kJ /kg-K2-Pa a static version of the Media.Steam package that was used in this paper can be found on GitHub at https: //github.com /lbl- srg/modelica-buildings, commit 8489367. Because the commonly-adopted IF97 water medium was de- veloped 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 formula- tions, as implemented in the Modelica Standard Library ver- sion 2.3 Media.Water.StandardWaterOnePhase package. Be- cause phase change was not to be included in the steam model, we used the typical single-phase implementation with pandT as independent variables. However, in order to improve the numerical e ciency of this medium for large-scale thermo-ﬂuid system simulations, we made some important modiﬁcations. First, we reduced the ap- plicable p-Trange 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 hand sfunctions 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 veriﬁed by users. Second, we replaced commonly called functions for energy and exergy analysis with ﬁrst and second order polynomial approximations. These two functions are speciﬁc enthalpy h(p;T) (Equation 4) and speciﬁc entropy s(p;T) (Equation 9), and their corresponding backward functions Equation 7 and Equation 12, which are analytical in- verses of the forward functions. Third, we implemented ex- pressions for the derivatives to facilitate symbolic processing of function derivatives in the numerical solver. Further, we couple the new steam medium with a numerically-e cient 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 chal- lenges 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. 3. 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 conﬁgura- tions are included in the component and district scale evalua- tions, which are summarized in Figure 1. Two cases involve the IF97 model with two-phase ﬂow – 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 mod- els from both of these libraries can be applied for steam DH applications, their original design intentions di er, and cor- respondingly, they produce di erent 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 evalua- tion. The remaining three cases are common to both. Further details regarding the model setup for each simulation case are provided in the following sections. Liquid Phase Vapor Phase Buildings.Media.Water ( Constant Water Density) Buildings.Media. Steam S+CW S+VW Buildings.Media.Specialized.Water. TemperatureDependentDensity ( Variable Water Density) Buildings.Media. Steam Case Key Details not modeled I(MBL) and I(MSL) Modelica.Media.Water.StandardWater (IF97 Model) Two-Phase Mixture I+VW Buildings.Media.Specialized.Water. TemperatureDependentDensity ( Variable Water Density) Modelica.Media.Water. StandardWaterOnePhase (IF97 Model) Figure 1: Medium models for component and district scale evaluations. Several new component and system models for DH model- ing 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 reﬁned and open-source released in the Modelica Buildings Li- brary. Beyond newly developed models, all components used are existing in the Modelica Standard Library v3.2 and Model- ica Buildings Library v7.0.0. 3.1. Thermodynamic Property Scale The accuracy and numerical performance of the proposed speciﬁc enthalpy and speciﬁc entropy approximations in Equa- tion 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 ef- fects from components, equipment, and systems. For accuracy, we calculated the residuals between the new approximation and the IF97 model across the entire p-Trange. In addition, we evaluated the Root Mean Square Error (RMSE) and Coe cient of Variation of the Root Mean Square Error (CVRMSE) using RMS E =rP(yiˆyi)2 Nand (13) 5 CVRMS E =RMS E ¯y; (14) where yiis the individual reference data generated by the IF97 model, ˆ yiis the corresponding evaluation data predicted by the new model, ¯ yis the mean of the reference dataset, and Nis the total number of data points. For numerical performance, handswere calculated over a dense grid of p-Tstates. The same simulation settings and computer were used across all evaluations in this paper (de- tails in subsection 3.4). The range of p-Tstates corresponded to the design limits of the new steam model. Modelica code for the speciﬁc enthalpy evaluation is given in Listing 1 as an example. With a ﬁxed pstate, 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 108for the new steam model). This model is simulated across the range of p, sweeping from 100 kPa to 550 kPa at in- crements of 250 Pa, for a total of 1,801 simulations. For com- parison with a reference model, the same script is repeated with the IF97 model by setting the Medium declaration to Model- ica.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 108and needed to be increased to 101for completion. Note, these error tolerance values are not simulation tolerances, but absolute error thresholds (parameter errin Listing 1). 3.2. Component Scale Because of the discontinuities at the phase-change barrier, numerical solvers have shown to experience di culties solv- ing thermo-ﬂuid problems with two-phase ﬂow. 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 us- ing IF97 and the standard mixing volume model in the Mod- elica Standard Library and Modelica Buildings Library. This component-scale experiment was selected to isolate some of the common numerical challenges of thermo-ﬂuid system mod- eling 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 m3that 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 pres- sure is set to 200 kPa. The mixing volume is replaceable in order to allow four di erent simulation runs to be completed with the same boundary conditions. Four cases from Figure 1 are included in the component- scale evaluation: I(MBL), I(MSL), S +CW, S +VW. First, the Modelica Buildings Library case with the IF97 model – I(MBL) – includes the mixing volume designed for single-phase ﬂuid:1model EvaluateEnthalpy 2 package Medium = Buildings.Media.Steam 3 " Steam medium model "; 4 parameter Modelica.SIunits.AbsolutePressure p =100000 5 " Pressure "; 6 parameter Modelica.SIunits.Temperature TMin = Medium.saturationTemperature (p) 7 " Minimum temperature for the simulation "; 8 parameter Modelica.SIunits.Temperature TMax =273 .15 +160 9 " Maximum temperature for the simulation "; 10 parameter Modelica.SIunits.MassFraction X[ Medium.nX ]= Medium.X_default 11 " Mass fraction "; 12 parameter Real err =1E -8 " Error "; 13 Medium.Temperature T " Temperature "; 14 Medium.ThermodynamicState state_pTX 15 " Medium state "; 16 Medium.ThermodynamicState state_phX 17 " Medium state "; 18 Modelica.SIunits.SpecificEnthalpy h 19 " Specific enthalpy "; 20 Modelica.SIunits.AbsolutePressure pMed 21 " Pressure "; 22 Medium.Temperature TMed " Temperature "; 23 Medium.BaseProperties basPro " Base properties "; 24protected 25 constant Real conv ( unit ="1/s")=1 26 " Conversion factor to satisfy unit check "; 27equation 28 // Compute temperatures used as input 29 T= TMin + conv * time *( TMax - TMin ); 30 // Check setting the state 31 state_pTX = Medium.setState_pTX (p=p,T=T,X=X); 32 state_phX = Medium.setState_phX (p=p,h=h,X=X); 33 // Compute properties 34 pMed = Medium.pressure ( state_pTX ); 35 TMed = Medium.temperature ( state_pTX ); 36 h= Medium.specificEnthalpy ( state_pTX ); 37 // Check the implementation of the functions 38 assert ( abs (p- pMed )<err ," Error in p."); 39 assert ( abs (T- TMed )<err ," Error in T."); 40 assert ( abs ( Medium.pressure ( state_pTX )- Medium.pressure ( state_phX ))<err ," Error in p state. "); 41 assert ( abs ( Medium.temperature ( state_pTX )- Medium.temperature ( state_phX ))<err ," Error in T state. "); 42 // Check base properties implementation 43 assert ( abs (h- basPro.h )<err ," Error in h."); 44 basPro.state.p =p; 45 basPro.state.T =T; 46annotation ( experiment ( Tolerance =1e -6, 47 StopTime =1 .0)); 48end EvaluateEnthalpy ; Listing 1: Modelica code for speciﬁc enthalpy evaluation. Buildings.Fluid.MixingVolumes.MixingVolume . Second, the Modelica Standard Library case with the IF97 model – I(MSL) – implements a two-phase equilibrium boiler model [36]: Modelica.Fluid.Examples.DrumBoiler.BaseClasses.Equilibri- umDrumBoiler . This model was designed for water-steam phase change and provides the physically-closest basis to the new implementation. The third and fourth cases involve the new split-medium implementation, but with di erent liquid water models. The new implementation includes a newly developed MixingVolumeEvaporation model as well as the new steam model presented in section 2. Both a simple liquid water model with constant density (S +CW) and a liquid water model with a variable temperature-dependent-density (S +VW) are implemented for comparison purposes. While the Modelica Buildings and Modelica Standard Library model details can be found in their respective open-source libraries, the MixingVol- umeEvaporation model is also a new contribution from this 6 Q_flow G=10 W/K T=300 °CK sV=VFixed TemperatureThermal ConductorHeat Flow Rate SensorReplaceable Mixing V olume Specific Entropy SensorBoundary with Reference Pressure(1) I(MBL) (2) I(MSL) (3) S+CW (4) S+VW}Figure 2: Diagram of Modelica model for mixing volume evaluations. work. Details regarding this model’s design principles and its mathematical formulation are included in Appendix B. 3.3. District Scale The objective of the district scale evaluation was to assess the accuracy and numerical performance of the new implemen- tation in a typical DH system design across districts of several sizes. An overarching description of the selected system is pre- sented next, followed by the Modelica implementation. 3.3.1. System Description Figure 3 depicts the schematic diagram for the DH evalua- tion. 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 de- signs seen in central plants [37], this conﬁguration was selected because it represents real-world control dynamics while repre- senting 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 ﬁxed pressure drops without any heat transfer. Lastly, the buildings contain a steam heat exchanger, a steam trap, and a conden- sate return pump. The building-side components are assumed to have steady energy and mass balances, with the mass ﬂow and heat transfer determined based on tabulated load data. For this DH case study, we adopted tabulated heating load proﬁles from the common exercise of IBPSA Project 1 [35] WP3.1, one of the work packages of the international collabora- tion 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]. 3.3.2. Modelica Implementation To evaluate the steam medium implementation across a va- riety of district sizes, a vector-style DH system model was de- veloped in Modelica (Figure 4). This top-level Modelica di- agram has a clear one-to-one relationship with the subsystemsdepicted in Figure 3. A parameter Nrepresenting the total num- ber of buildings can be adjusted to represent districts of multi- ple sizes. The tabulated heating load proﬁle for each building is scaled by 1=Nin order to keep the central plant sizing constant across all simulation cases. This DH system is simulated for the ﬁrst 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 de- signed 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 or- der to prevent unintended reverse ﬂow. The boiler model con- tains the MixingVolumeEvaporation to represent to ﬂuid vol- ume. The boiler also includes pressure drop in the working ﬂuid and represents the heat transfer with the ambient environment and with the fuel through heat ports. Both the thermal conduc- tance 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 theTabulated Heating Load data reader, the variable heat ﬂow rate ˙Qfrom IBPSA Project 1 WP3.1 is input directly. A con- densate return pump is modeled with a prescribed mass ﬂow rate of ˙ m=˙Q=hf g;std, where hf g;stdis the standard enthalpy of vaporization at the design pressure. The Heat Exchanger Volume model is conﬁgured similarly to the MixingVolumeE- vaporation model, except it represents a condensation process instead, and the steady state balance is adopted for this instance. Because the heat ﬂow rate is directly proportional to the mass ﬂow rate (given that the change in enthalpy is essentially ﬁxed), a heat port is not permissible with the steady balance. Lastly, theSteam Trap represents an isenthalpic process where liquid condensate is discharged at atmospheric pressure, while any ﬂashed steam vapor due to the drop in pressure is returned to a liquid state before discharge. Similarly to subsection 3.2, four medium model conﬁgura- tions 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 TemperatureDependentDen- sityfrom the Modelica Buildings Library. 3.4. Simulation Settings All simulations ran in Dymola 2021 on a Windows 10 work- station 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 com- pared to other numerical solvers. The simulation tolerance was set to 106for all cases as well. 4. Evaluation Results and Discussion This section presents the results from simulations across the three scales: thermodynamic properties, component, and dis- 7 Pump speed Measured pressure Part load ratio Steam Boiler Condensate Return Pump Building N ... ...Steam Supply Pipe (lossless) Condensate Return Pipe (fixed resistance) Condensate Return Pump (Ideal mass flow control) Trap TrapHeat Exchanger Heat Exchanger Feedwater Tank VentFeedwater Pump s wMeasured liquid volume Building 1 Central Plant Distribution Network Liquid volume setpoint Pressure setpoint Figure 3: Schematic diagram of DH steam loop with a central plant, distribution network, and Nnumber of interconnected buildings. Central Plant s wBuildings [:] Distribution Network Figure 4: Top level diagram of Modelica model for the DH system. trict scales. First the accuracy of the model will be presented, followed by the numerical performance. 4.1. 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. 4.1.1. Accuracy Property evaluations for speciﬁc enthalpy and entropy pro- duced acceptable accuracy across the entire p-Trange. As shown in Figure 7, the largest residuals for h(p;T) occurred near the saturation curve and at low p-Tconditions. The reason for this is that the nonlinearities in h(p;T) are higher in these re- gions 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-Tconditions and the middle of the saturation curve, respec- tively. Across the entire range, the RMSE was 0.76 kJ /kg and MyPIk=1/VBoiWatSet PI k=1/pSteSetpk=1 s sPump Controller Feedwater TankInflowing Port (Liquid Water) Feedwater PumpCheck Valve Boiler ControllerNormalized Liquid Volume Normalized PressureBoilerOutflowing Port (Steam Vapor)Pump Power Fuel Heat Flow Rate s wFigure 5: Diagram of Modelica model for the central plant. I k=1 Mm_flowk=1/dhVapStd s sTabulated Heating Load Heat Exchanger V olume (Steady State)Steam TrapCondensate Return PumpMass Flow RateIntegralHeat Flow Rate Heating Energy Inflowing Port (Steam Vapor)Outflowing Port (Liquid Water) s w Figure 6: Diagram of Modelica model for the interconnected building. 8 CVRMSE was less than 0.0001%. Fitting the linear approxi- mation using a least square approach, the R2value was 0.999. Overall, a linear surface ﬁt is able to capture h(p;T) well in this p-Trange. 100 110 120 130 140 150 160 Temperature (C) 100150200250300350400450500550Pressure (KPa)Specific Enthalpy Residuals (kJ/kg) Saturation Curve 2.4 1.8 1.2 0.6 0.00.61.21.82.4 Figure 7: Residuals of the polynomial approximation for h(p;T). As shown in Figure 8, the largest residuals for s(p;T) oc- curred at high pressure-temperature conditions. The shape of the residuals indicate that the errors are relatively consistent with respect to temperature, but the quadratic ﬁt in pressure underestimates sslightly in the middle of the range while it overestimates sat the high and low limits. This shape is due to the fact that a low order polynomial is being ﬁt to a nonlin- ear 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 lin- ear approximation using a least square approach, the R2value was 0.996. It is important to note that the distribution of p-T grid points selected for the least squares optimization inﬂuences the shape of these ﬁts. Based on the relatively high accuracy achieved across the entire range of conditions, we can conclude that the selected grid was appropriate. 100 110 120 130 140 150 160 Temperature (C) 100150200250300350400450500550Pressure (KPa)Specific Entropy Residuals (kJ/kg-K) Saturation Curve 0.1500 0.1125 0.0750 0.0375 0.00000.03750.07500.11250.1500 Figure 8: Residuals of the polynomial approximation for s(p;T).4.1.2. Numerical Performance The results in Table 5 show that the new polynomial imple- mentations 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 equa- tions 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 simpliﬁed implementa- tions for hand sproduced numerical savings at the thermo- dynamic property scale, the savings are often compounded for larger thermo-ﬂuid system models that involve nonlinear sys- tems of equations. This will be evaluated with the district scale models. However, the computing time savings and calcula- tion accuracies demonstrate how fundamental improvements in steam property modeling can be achieved through simple func- tion replacement. Table 5: Computing time savings at thermodynamic property scale. Results given are average computing times over 10 simulation runs. PropertyMedium Computing Time Savings Model Time (s) (%) hIF97 16.5 – Linear Eq. (4) 10.0 39.4 sIF97 15.1 – Quadratic Eq. (9) 13.5 10.3 4.2. 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. 4.2.1. Accuracy Figure 9 depicts the temperature and mass results of the mix- ing volume evaluation cases. In this simulation, the mixing volume ﬂuid is heated via a constant, high temperature bound- ary. With the I(MBL) case, undesirable erratic behavior is clear. At initialization, the liquid in the volume is at 40C. For the ﬁrst 28 minutes, the temperature gradually increases. However, when water starts to boil, the ﬂuid temperature oscillates be- tween the saturation temperate and the boundary temperature of 300C. Mass and density are also oscillating in the I(MBL) case, as the ﬂuid 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 ﬂuid (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, 9 this numerical problem can be avoided if an explicit, ﬁxed- time step method is employed (i.e., Euler) rather than an im- plicit, variable-time step method (i.e., DASSL). With Euler, the I(MBL) case performs physically correct and avoids the numer- ical chattering issue. Further, the I(MBL) mixing volume can simulate not only a saturated ﬂuid, but supercooled and super- heated ﬂuids as well, which may be of interest for some use cases. However, implicit, variable-time step methods are often used for large thermo-ﬂuid system simulations because of their ability to deal with sti systems. Thus, the Euler results were excluded from this paper. Figure 9: Evolution of (a) ﬂuid temperature and (b) ﬂuid density, including both the liquid and vapor components, and (c) ﬂuid mass in the mixing volume through the boiling process. 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 ﬂuid at a satu- rated 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 liq- uid into vapor. In the I(MSL) case, the liquid water continues into a ”negative” mass state after all of the liquid has boiled o . This is physically not correct. Conversely, the new mix- ing volume implementation is designed to stop the simulation when no more liquid water is present because it is conﬁgured 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 thisassumption is valid for steam boilers during normal operation in most heating applications [36]. With the new numerically e cient implementation, some mi- nor accuracy errors are introduced. Since the steam implemen- tation 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 com- mon 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 tem- perature range, causing higher than normal errors to be intro- duced. However, with that said the S +CW case did produce less than 1% errors in terms of CVRMSE for most property evalu- ations. The largest errors were introduced in the S +CW case forwandm, with CVRMSE values of 5.9% and 7.6% respec- tively. 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. Table 6: RMSE and CVRMSE of mixing volume evaluations relative to the I(MSL) drum boiler with standard water. PropertyRMSE CVRMSE (%) S+CW S +VW Units S +CW S +VW w 58.6 1.9 kg /cm35.9 0.2 s<0.001<0.001 kg /cm3<0.001<0.001 h 1.5 1.5 kJ /kg 0.06 0.06 u 2.5 2.3 kJ /kg 0.4 0.4 s 0.005 0.005 kJ /kg-K 0.07 0.07 T 0.005 0.005 K 0.001 0.001 m 0.59 0.02 kg 7.6 0.3 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 mod- eling sensitivity at low mass quantities and mass ﬂow rates, par- ticularly with large, complex, DH system models where prob- lems can be di cult to diagnose. One solution to this problem can be to include smoothing functions that account for regular- ization at low mass conditions to avoid diverging total property evaluations such as U. This will be implemented in future ver- sions of the steam DH components in the Modelica Buildings Library. 4.2.2. Numerical Performance With this small simulation case, there were negligible dif- ferences in computing time between the four cases. However, the structure of the di erential algebraic system of equations indicates the likelihood of numerical performance di erences 10 for larger thermal-ﬂuid system models. While each of the four cases contain one continuous time state variable, the nonlin- ear systems of equations varies. The I(MBL) case contains one nonlinear system involving two iteration variables: du=dtand d=dt. The I(MSL) case also contains a nonlinear system but with one iteration variable: dVw=dt. Conversely, the new imple- mentations 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 experi- ences, the time required to iteratively solve the nonlinear sys- tems largely contribute to the total computing time. Thus, this likely will have an impact on computing time for DH applica- tions, which will be quantitatively tested in the next section. 4.3. 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. 4.3.1. Accuracy Figure 10 shows that the measured heat ﬂow rate at the build- ing for each of the four cases followed the input data ﬁle. Given that the input to the building model was the data table and the condensate return pumps controlled the mass ﬂow rate ide- ally, this was generally expected. The error di erences in the S+CW and S +VW cases are due to the building model imple- mentation. Because the mass ﬂow rate is prescribed based on ˙m=˙Q=hf g;std, errors in the hf g;stdcalculation a ect the mass ﬂow rate. This error could be resolved by changing the conden- sate return pump control to use a PI controller that adjusts the measured heat ﬂow rate to track the tabulated input value. 170 180 190 200 210 220 230 240 Time (h)02468Heat Flow Rate (kW) Input Data I(MSL) I+VW S+CW S+VW Figure 10: Heat ﬂow rate compared to the input data ﬁle for the four cases. As seen in Table 7 and Table 8, the accuracy of thermody- namic property evaluations had insigniﬁcant di erences across all cases. Here, the building 1 results are used as proxies for all simulation cases, since the accuracy did not signiﬁcantlychange building to building, with respect to N, and at di er- ent points throughout the system. Subscripts inandoutrefer to the inlet and outlet ports of the building model, where the inlet state is steam vapor and the outlet state is liquid water con- densate. 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 pro- duced some of the higher errors, all errors were within accept- able ranges. For example, the CVRMSE values are within 1% for all property evaluations. With this split medium approach and simple polynomial approximations for hands, these DH simulations were able to produce acceptable accuracy for all thermodynamic property states. Table 7: RMSE of building 1 thermodynamic properties from DH evaluations. Property S +CW S +VW I +VW Units Tin 0.007 0.007 0.001 K pin 7.49 8.58 5.79 Pa hin 1.74 1.74 0.001 kJ /kg sin 0.007 0.007 0.0 kJ /kg-K Tout 0.017 0.017 0.017 K pout 7.49 8.58 5.79 Pa hout 2.79 2.79 2.79 kJ /kg sout 0.007 0.007 0.007 kJ /kg-K Table 8: CVRMSE (%) of building 1 thermodynamic properties from DH eval- uations. Property S +CW S +VW I +VW Tin 0.002 0.002 0.0 pin 0.002 0.003 0.002 hin 0.064 0.064 0.0 sin 0.11 0.11 0.0 Tout 0.004 0.004 0.004 pout 0.002 0.003 0.002 hout 0.50 0.50 0.50 sout 0.40 0.40 0.40 Throughout the 15 day simulation, the di erence in total fuel consumption by the boiler was less than 2% for all cases. How- ever, interestingly, slight di erences in thermodynamic states caused short term variance (e.g., on the order of 1 to 5 hours) in boiler fuel ﬂow rate (see Figure 11). With this system con- ﬁguration, the boiler fuel ﬂow 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 rela- tive pressure at the boiler outlet does not follow the same exact trend line for all four simulation cases. In this ﬁgure, relative pressure prelis calculated as prel=pmea pset1; (15) where pmeais the measured pressure at the boiler outlet and pset 11 -1.0e-04-5.0e-050.0e+005.0e-051.0e-041.5e-04Relative Pressure(a) Setpoint 170 180 190 200 210 220 230 240 Time (h)0.02.55.07.510.012.5Fuel Flow Rate (kW)(b)I(MSL) I+VW S+CW S+VWFigure 11: Boiler control (a) input relative pressures and (b) output fuel ﬂow rates. The relative pressure is calculated using Equation 15. is the pressure setpoint of 300 kPa. Note that prelin 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, particu- larly with realistic variable-load control. This DH evaluation revealed the importance of well-tuned feedwater pump and boiler controls on power and fuel ﬂow rates. One major beneﬁt of Modelica is the ability to simulate realistic control of complex thermo-ﬂuid system models. While these controlability aspects did not greatly inﬂuence the over- all 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. 4.3.2. Numerical Performance Table 9 presents the structure of the translated model for each of the four cases. The number of equations varied slightly be- tween 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 nonlin- ear systems present in the district models. For the I(MSL) case, both the number of linear systems and the dimension of the ﬁrst nonlinear increase with the number of buildings. Conversely, the linear and nonlinear systems do not change with respect to Nfor the three split medium cases. This is advantageous par- ticularly for simulations of large DH systems. As seen in Figure 12, computing times were signiﬁcantly re- 0.0 0.2 0.4 0.6 0.8 1.0 log10(Number of Buildings)0.500.751.001.251.501.752.00log10(Computing Time) (a)I(MSL) I+VW S+VW S+CW 1.2 1.3 1.4 1.5 1.6 1.7 log10(Number of Continuous Time States)0.500.751.001.251.501.752.00log10(Computing Time) (b)Figure 12: Computing time in seconds, presented on a log-log scale (base 10), with respect to (a) the total number of interconnected buildings and (b) the number of continuous time states. 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 stud- ies involving small districts, where linear and quadratic scaling was present [6, 39]. However most importantly, the rate of com- puting 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 quantiﬁed computing time savings in this study. We are able to achieve this result with the new steam model, the liquid water models in the Model- ica 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 equa- tions 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. 5. Conclusion This work introduces a new steam medium model that en- ables numerically e cient DH simulations. Accuracy and nu- merical performance were aptly balanced by reducing the ap- 12 Table 9: Translation results with respect to the number of buildings in the DH system model. Results given are after Dymola’s built-in model manipulation. NI(MSL) /I+VW /S+VW /S+CW Equations Linear Systems Nonlinear Systems 1 767 /769/769/773 f2,2g/f2g/f2g/f2g f 2,4g/f3g/f3g/f1g 2 1,132 /1,134 /1,134 /1,140f2,2,2g/f2g/f2g/f2g f 3,4g/f3g/f3g/f1g 3 1,497 /1,499 /1,499 /1,507f2,2,2,2g/f2g/f2g/f2g f 4,4g/f3g/f3g/f1g 4 1,862 /1,864 /1,864 /1,874f2,2,2,2,2g/f2g/f2g/f2g f 5,4g/f3g/f3g/f1g 5 2,227 /2,229 /2,229 /2,241f2,2,2,2,2,2g/f2g/f2g/f2g f 6,4g/f3g/f3g/f1g 6 2,592 /2,594 /2,594 /2,608f2,2,2,2,2,2,2g/f2g/f2g/f2g f 7,4g/f3g/f3g/f1g 7 2,957 /2,959 /2,959 /2,975f2,2,2,2,2,2,2,2g/f2g/f2g/f2g f 8,4g/f3g/f3g/f1g 8 3,322 /3,324 /3,324 /3,342f2,2,2,2,2,2,2,2,2g/f2g/f2g/f2g f 9,4g/f3g/f3g/f1g 9 3,687 /3,689 /3,689 /3,709f2,2,2,2,2,2,2,2,2,2 g/f2g/f2g/f2g f 10,4g/f3g/f3g/f1g 10 4,052 /4,054 /4,054 /4,076f2,2,2,2,2,2,2,2,2,2,2 g/f2g/f2g/f2g f 11,4g/f3g/f3g/f1g plicable p-Trange to suit the majority of DH application; re- placing handsfunctions with invertible, linear and quadratic polynomial ﬁts; and developing component, equipment, and system models that implement a split-medium approach with separate models for liquid and vapor phases. Modelica models were developed and simulated to evaluate the steam medium at thermodynamic property, component, and district scales. At the thermodynamic property scale, handscalculations fell within 2.1 kJ /kg (0.01%) and 0.13 kJ /kg-K (2.0%), respectively, while reducing computing time by 34%-40%. At the component scale, the new implementation avoided a chattering behavior while most property evaluations were within 0.4% of the more computationally-expensive IF97 case. At the district scale, energy consumption for the most numerically-e cient S +CW was within 2% of the I(MSL) case. Signiﬁcant computing time savings on the order of 5.6 to 9.3 times were achieved through the new implementation, while all thermodynamic properties were calculated within 0.5% of the IF97 implementation of I(MSL). Furthermore, the nonlinear system of equations scaled linearly with problem size for the I(MSL) case while it is held constant with our new implementation, which is critically im- portant for large district simulations. The new steam model is open-source released in the Mod- elica 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 be- tween 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 vari- able temperature-dependent-density water model is more suit- able for hourly optimization of steam plants and districts. Re- gardless of the speciﬁc application, this new steam model can help aid the transition of legacy DH systems to newer sustain- able designs, while providing a pathway for district-scale anal-ysis and optimization in the many steam heating applications that are likely to remain. 6. 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), O ce of Energy E - ciency and Renewable Energy, and Building Technologies Of- ﬁce. This program is managed by Oak Ridge National Labora- tory (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 reﬂect the policies and views of DOE, ORNL, ORAU, or ORISE. In addition, this material is based upon work supported by the DOE’s O ce of Energy Ef- ﬁciency and Renewable Energy under the Advanced Manufac- turing O ce, award number DE-EE0009139, and the Building Technologies O ce, contract number DE-AC02-05CH11231. Further, this work emerged from the IBPSA Project 1, an in- ternational project conducted under the umbrella of the Interna- tional Building Performance Simulation Association (IBPSA). Project 1 will develop and demonstrate a BIM /GIS and Mod- elica Framework for building and community energy system design and operation. 7. 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 as- sumes any legal liability or responsibility for the accuracy, com- pleteness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not in- fringe privately owned rights. Reference herein to any speciﬁc commercial product, process, or service by trade name, trade- mark, manufacturer, or otherwise does not necessarily consti- tute or imply its endorsement, recommendation, or favoring 13 by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not neces- sarily state or reﬂect 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, rel- evant formulations for speciﬁc enthalpy and entropy are repro- duced here. In the IF97 formulation for Region 2 (vapor), basic equations for Gibbs free energy gare expressed in a dimension- less form (; )=g(p;T) RT= 0(; )+ r(; ); (A.1) where 0represents the ideal-gas part and rrepresents the residual part. This equation is formulated in terms of the spe- ciﬁc gas constant of ordinary water R=0:461526 kJ=(kgK), the reduced pressure =p=p, and the inverse reduced tem- perature =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 0is 0(; )=ln+9X i=1n0 i J0 i; (A.2) where each of the 9 coe cients for n0 iandJ0 iare tabulated in Table 10 of [40]. Similarly, the form of the residual part ris given as r(; )=43X i=1niIi( 0:5)Ji; (A.3) where each of the 43 coe cients for niand exponents IiandJi are tabulated in Table 11 of [40]. In terms of Gibbs free energy, speciﬁc 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 0and rwith respect to at constantare "@ 0 @ # =9X i=1n0 iJ0 i J0 i1and (A.5) "@ r @ # =43X i=1niIiJi( 0:5)Ji1: (A.6) In terms of Gibbs free energy, speciﬁc 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 coe cients, and 95 exponents, while s(p;T) includes 7 equations, 104 coe cients, and 190 exponents. In addition to the forward equations, separate backward equations are similarly implemented. However, Region 2 is further di- vided into three subregions 2a, 2b, and 2c (see [40] for sub- region divisions). Thus, three T(p;h) and three T(p;s) cover Region 2, but only subregion 2a is applicable for the reduced p-Trange speciﬁed for DH applications in Table 1. The dimen- sionless backward equation T(p;h) for subregion 2a is T(p;h) T=34X i=1niIi(2:1)Ji; (A.8) where=h=h,h=2000 kJ /kg, and each of the 34 coe - cients for niand exponents IiandJiare tabulated in Table 20 of [40]. Lastly, the backward function T(p;s) for subregion 2a is T(p;s) T=46X i=1niIi(2)Ji; (A.9) where=s=s,s=2 kJ=(kgk), and each of the 46 coe - cients for niand exponents IiandJiare tabulated in Table 25 of [40]. In total, T(p;h) includes 34 coe cients and 68 expo- nents, while T(p;s) includes 46 coe cients and 92 exponents. Modelica implementations of the above forward and backward equations for speciﬁc 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, rep- resents the evaporation process of water with liquid and vapor phases in equilibrium. It is designed for the ﬂuid to be at a sat- urated state only, which is a reasonable assumption for steam boilers operating under normal conditions. The fundamental equations are as follows. The ﬂuid mass min the volume is calculated as m=sVs+wVw; (B.1) whereis density, Vis volume, and subscripts sandwrepre- sent the steam and liquid water components. The total internal energy Uis U=sVshs+wVwhwpV; (B.2) where his speciﬁc enthalpy, pis pressure, and V=Vs+Vwis the total volume of the ﬂuid. The MixingVolumeEvaporation is conﬁgured to allow both steady state and dynamic mass and energy balances. The dy- namic mass balance is given as dm dt=˙ms+˙mw; (B.3) 14 where ˙ msand ˙mware the mass ﬂow 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+˙mshs+˙mwhw; (B.4) where ˙Qis the net heat ﬂow 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 qual- ity of 1 is discharged from the outlet port with hs= hv. Second, we assume that the liquid and vapor compo- nents in the volume are at equilibrium. Further, it is im- portant to note that the new MixingVolumeEvaporation is very similar to Modelica.Fluid.Examples.DrumBoiler. Base- Classes.EquilibriumDrumBoiler [36] with the following excep- tions: 1. Rather than a two-phase medium, ﬂuid mediums are mod- eled as two single-state ﬂuids, with liquid water at the up- stream 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.