Many commercial solar thermal power plants rely on indirect thermal storage systems in order to provide a stable and reliable power supply, where the working uid is commonly thermal oil and the storage uid is molten salt. The thermal oil - molten salt heat exchanger control strategies, to charge and discharge the thermal storage system, strongly a ect the performance of the whole plant. Shell-and-tube heat exchangers are the most common type of heat exchangers used in these facilities. With the aim of developing advanced control strategies accurate and fast dynamic models of shell-and-tube heat exchangers are essential. For this reason, several shell-and-tube heat exchanger models with di erent degrees of complexity have been studied, analyzed and validated against experimental data from the CIEMAT-PSA molten salt test loop for thermal energy systems facility. Simulation results are compared in steady-state as well as transient predictions in order to determine the required complexity of the model to yield accurate results. Keywords: Multi-pass shell-and-tube heat exchanger, dynamic simulation, transient response, thermal energy storage, molten salt, Modelica. 1. Introduction Heat exchanger modeling and simulation has been extensively addressed in the literature due to its im- portance in industrial applications. Many kinds of heat exchanger models have been developed targeting di erent aspects of interest such as, steady state and transient predictions, one-dimensional process and three-dimensional detailed models, simpli ed analyt- ical models, linearized models, etc. In [1], it is pre- sented an exhaustive review of heat exchanger models dealing with di erent elds of interest. Special men- tion deserve shell-and-tube heat exchangers for their Corresponding authorrole in process engineering, thus being these the sub- ject of numerous publications dealing with their e- cient modeling. In [2], a detailed modeling review of this kind of heat exchangers is presented. This work focuses on the comparison and evalu- ation of several shell-and-tube heat exchanger mod- els with di erent level of complexity and with the purpose of process simulation and control design. In both areas, not only reliable and accurate but also fast heat exchanger models are required. The simulation results from these models have been compared against experimental data from a multi-pass shell-and-tube heat exchanger. This heat exchanger is part of the CIEMAT-PSA molten salt test loop for thermal energy systems (MOSA) facil- Preprint submitted to Applied Thermal Engineering September 19, 2018 Figure 1: CIEMAT-PSA MOSA facility ity, a facility with the aim of studying Thermal En- ergy Storage (TES) systems in solar thermal power plants. Many factors such as, environmental issues and concerns about sustainability are presently encourag- ing investment and research into alternative sources that might provide renewable, clean, sustainable and ecient energy. Solar thermal power plants are ap- propriate for large-scale energy production due to their dispatchability on demand, since they eciently store heat in TES systems. Therefore, many commer- cial solar thermal power plants rely on TES systems to provide a stable and reliable power supply [3, 4, 5]. The performance of solar thermal power plants with TES systems is highly in uenced by the heat exchanger control strategies applied in the charging and discharging processes [6]. Therefore, advanced control strategies may improve the performance of the whole plant. Accurate and fast dynamic plant models are required in order to design and test con- trol strategies. For this reason, the heat exchanger model is an important part in the design, testing and validating of advanced control strategies for TES sys- tems in solar thermal power plants. This paper is organized as follows, section 2 intro- Figure 2: Two-unit multi-pass shell-and-tube heat exchanger duces the experimental plant and the heat exchanger under consideration, section 3 describes the di erent heat exchanger models considered in this work, sec- tion 4 presents simulation results from several heat exchanger models against experimental data at de- sign conditions, steady-state and transient predic- tions. Finally, section 5 draws the main conclusions and summarizes future work. 2. Experimental Plant: CIEMAT-PSA MOSA Facility A molten salt testing facility was set up at Plataforma Solar de Almer a (PSA), division of CIEMAT, the public research center for Energy, En- vironmental and Technological Research, which is owned by the Spanish government. The aim of this facility is to study TES systems, evaluate and con- trol the heat exchange between the molten salt and di erent kinds of potential heat transfer uid in solar thermal power plants (i.e. thermal oil and pressur- ized gases). In order to study pressurized gases, the MOSA facility is coupled to the innovative uids test loop facility [7] by means of a CO 2- molten salt heat exchanger. This last facility comprises two parabolic- trough collectors. The main components of the MOSA facility ( g- ure 1) are brie y described in the following list. Con- sult [8] and [9] for further details. •Molten salt tanks. Hot and cold molten salt tanks reproduce the sensible-heat TES systems 2 Shell-side inlet nozzle Shell-side outlet nozzleTube-side outlet nozzle Tube-side inlet nozzleVertical segmental ba es Longitudinal ba eFigure 3: Heat exchanger schematic representation [10] of commercial solar power plants. The hot tank is at the ground level whereas the cold tank is under ground level. •Molten salt air cooler. This air cooler replicates the molten salt discharging process by cooling down the molten salt. •CO2- molten salt heat exchanger. This heat ex- changer allows exchanging heat from pressurized gases from the innovative uids test loop facility. •Two anged pipe sections . Components can be installed in the two anged pipe sections in order to be tested in the molten salt circuit under real working conditions. •Electrical heat-tracing system. Its purpose is to prevent salt freezing. •Thermal oil loop. It allows storing and releasing thermal energy to/from the molten salt. This loop includes the following components: a ther- mal oil expansion tank, a centrifugal pump, anoil heater, a thermal oil - molten salt heat ex- changer, thermal oil air cooler, an expansion tank and nitrogen bottles to render the molten salt and thermal oil inert. The purpose of the oil heater is to provide the same amount of heat than parabolic-trough collectors. Therefore, the oil heater can be used to emulate them and repli- cate transients such as, start-ups, shutdowns and cloud disturbances. The multipurpose MOSA facility is exible and can work in four di erent operating modes, they are sum- marized as follows. Further details about operating modes can be consulted in [9]. •Mode 1. In this mode, energy coming from the innovative uids test loop is used to charge the molten salt TES system. •Mode 2. The molten salt is cooled down by means of the air cooler system in this mode. •Mode 3. In mode 3, the TES system is charged from thermal energy of the thermal oil loop by means of the thermal oil - molten salt heat ex- changer. •Mode 4. This mode discharges the TES system and thus heating up thermal oil by means of the same heat exchanger than in mode 3. 2.1. Two-unit Multi-Pass Shell-and-Tube Heat Ex- changer This work focuses on the modeling of the thermal oil loop heat exchanger. This heat exchanger is com- posed of two counter- ow multi-pass shell-and-tube Table 1: Heat exchanger nominal conditions in mode 3 Feature Shell side Tube side Fluid Solar salt VP-1 Inlet mass ow (kg =s) 2.08 1.57 Inlet pressure (bar) 2 14 Outlet pressure (bar) 1.6 13.97 Inlet temperature (C) 290 380 Outlet temperature (C) 373 313 3 Table 2: Heat exchanger design parameters Design parameter Value Area of heat transfer 83:02 m2 Mean heat transfer coef. 238:5 W=(m2K) Tube-side volume 0:111 m3 Shell-side volume 0:255 m3 Tubes per shell pass 153 Tube outer diameter 12:7 mm Tube inner diameter 12:5 mm Tube length per shell pass 3:4 m Shell-side passes per unit 2 Ttube-side passes per unit 2 Baes per shell pass 39 Bae spacing 80 cm Last bae spacing 90 cm Bae cut 17.35 % units, see gure 2 for a real picture of the system and gure 3 for a schematic representation. The shell- side uid is molten salt, solar salt in particular (60 % NaNO 3and 40 % KNO 3) whereas the tube-side uid, due to its high pressure, is the commercial Therminol VP-1 thermal oil, see table 1 for nominal operating conditions. Each unit of the heat exchanger was designed fol- lowing a Tubular Exchanger Manufacturers Associa- tion (TEMA) design, in paticular a N-type front end stationary head, F-type shell and U-type rear end sta- tionary head (NFU) design. Both units are tilted 2 ° in order to facilitate their drainage. The F-type shell has two shell passes de ned by a longitudinal bae as well as two tube passes in U shape (see gure 3). The F-type shell is the most common and economical heat exchanger design used at commercial parabolic- trough solar thermal power plants [11]. Thirty-nine vertical segmental baes per shell pass, with 17.35 % vertical bae cuts (see gure 4), force the shell-side uid to follow a S-shaped path (see gure 5) in or- der to increase the convective heat transfer coecient which has its highest value in cross ow. In counter ow (see gure 3), the tube-side uid enters the inlet nozzle, ows along the tube bundle turning around due to the longitudinal bae and the U-tube design, nally leaving the heat exchanger through the outlet nozzle. Consult table 2 to obtain information about the main design parameters of the heat exchanger. Vertical ba e cut Longitudinal ba eVertical segmental ba e Tube bundleSealing stripDrainage channel Figure 4: Heat exchanger unit cross-sectional area [10] Longitudinalba e Vertical segmental ba esU shape S shapeVertical ba e cut Figure 5: Heat exchanger unit cross-sectional area [10] 4 3. Shell-and-Tube Heat Exchanger Models with Di erent Degree of Complexity Several heat exchanger models have been studied and analyzed. General assumptions considered in all the models are summarized in the following list. •Heat conduction and radiation are negligible in the uids. Axial heat ow is also negligible in both uids. •Thermal conductivity in the tube walls is in - nite in normal direction of ow and zero in ow direction. •The thermal capacitance of the heat exchanger shell is neglected. •One-dimensional treatment in the direction of ow. The Modelica language [12] has been used for the modeling of the heat exchanger. Modelica was devel- oped and is maintained by the Modelica Association, a non-pro t and non-governmental international as- sociation. This modeling language has been designed to model conveniently complex physical systems be- cause the language supports the object-oriented and equation-based paradigms. Any medium model from the Modelica Media li- brary [13], or new implemented models which takes into account this library interface, can be used to calculate the uid thermodynamic properties in the shell side as well as in the tube side. Both thermody- namic property mediums have been implemented in Modelica because they are not available in the Mod- elica Media library. Solar salt has been implemented according to the thermodynamic properties available in [14, 15] and Therminol VP-1 thermodynamic prop- erties according to [16]. Both uids depend only on temperature. The studied models are listed as follows, they are explained in the following subsections. Nomenclature is shown in table 3. •Quasi-Steady-state model (QSS model). •Dynamic model (DY model).•Dynamic with Thermal Capacitance Model (DYTC model). •Dynamic with Thermal Capacitance and Cell Method Model (DYCM model). 3.1. Quasi-Steady-state Model (QSS model) The Quasi-Steady-state (QSS) model is the easi- est heat exchanger model. It is an algebraic model centered on a speci c operating point for general counter- ow heat exchangers. It is based on the mathematical development presented in [17]. This model simpli es the heat exchanger considering it as a concentric tube system where the uids ow in op- posite directions while transferring sensible heat. It is a lumped parameter model which assumes that lo- cally every state at the heat exchanger is a steady- state. The main particular assumptions, besides those stated in section 3, are the following. •Thermal capacitance of the tube bundle metal structure has been neglected. •Fluid friction against shell and tube walls is ne- glected. There is neither pressure drop at the shell side nor the tube side of the heat exchanger. •Fluid incomprehensibility is assumed. Fluid mass storage at the heat exchanger is neglected, thus inlet and outlet mass ow rate are equal. •Uniform thermodynamic properties are as- sumed. Speci c heat capacities are evaluated at the mean temperature between inlet and outlet temperatures. In order to break the algebraic loop between the in- let and outlet temperatures and to prevent zero mass ow rate, eqs. 1 and 2 are considered. Ttb;out =Ttb;intb(Ttb;inTsh;in); (1) Tsh;out =Tsh;in +sh(Ttb;inTtb;out ); (2) wheretbandshare two dimensionless factors whose values are non-zero when the mass ow rates ( _mtband _msh) have positive values. They are calcu- lated by eqs. 3 and 4. 5 Table 3: Nomenclature Latin letters Var Description Units A Surface area [m2] C Heat capacity [J =K] cp Speci c heat capacity [J =(kg K)] d Diameter [m] D Characteristic dimension [m] f Friction factor [-] G Mass velocity [kg =(m2s)] _H Enthalpy ow rate [W] j Colburn j factor [-] K Thermal conductivity [W =(m K)] l Length [m] m Mass [kg] _m Mass ow rate [kg =s] n Number [-] Nu Nusselt number [-] p Pressure [Pa] Pr Prandtl number [-] _Q Heat ow rate [W] Re Reynolds number [-] T Temperature [K] t Time [s] U Internal Energy [J] Greek letters Var. Description Units Heat transfer coecient [W =(m2K)]  Dimensionless factor [-]  Dimensionless factor [-] Subs. Description Subs. Description amb Ambient av Average b Bae he Heat exchanger in Inlet loss Thermal losses pa Pass out Outlet sh Shell side tb Tube side w Wall tb=1ehe 1_mtbcp;tb _mshcp;shehe; (3) sh=_mtbcp;tb _mshcp;sh: (4)The dimensionless factor ( he) is calculated as fol- lows, he= heAw1 _mtbcp;tb1 _mshcp;sh ; (5) where the overall heat transfer coecient ( he) can be calculated by eq. 6. 1 he=1 sh+1 tb: (6) If the thermal capacitance of the tube walls is high, the model cannot capture the transient behav- ior of the system. Furthermore, this model neglects the uids mass storage in the heat exchanger, and therefore their thermal capacitance. Eqs. 1 and 4 do not consider heat losses between the shell-side uid and the environment. However, thermal losses can be included in the model by means of the follow- ing mathematical reformulation. Heat losses are as- sumed as a ctitious mass ow rate ( _ msh;loss ) with the same thermodynamic properties as the shell-side uid. This mass ow rate is added to the total mass ow rate to properly calculate the temperature. _msh= _msh;in + _msh;loss: (7) In order to assure a correct energy balance, this mass ow rate is calculated by eq. 8. _msh;loss = lossAamb cp;shTsh;avTamb Tsh;outTsh;in ;(8) where heat losses are calculated by Newton's law of cooling between the ambient temperature ( Tamb), and the average shell temperature ( Tsh;av = (Tsh;in+ Tsh;out )=2). Heat transfer coecients ( sh, tband loss) are values that need to be calibrated with ex- perimental data. 3.2. Dynamic Model (DY model) The Dynamic (DY) model described in this sec- tion was developed considering that in [18]. Same assumptions than those presented in the QSS model hold for this model, besides uid mass storage at the heat exchanger which is taken into account in the 6 DY model. This is a dynamic distributed parameter model, where each cell or Control Volume (CV) is a small lumped parameter counter- ow heat exchanger model. Therefore, eqs. 9 and 10 represent the energy balance for the tube side and the shell side respec- tively in each cell of the DY model. CtbdTtb;out dt= _mtbcp;tb(Ttb;inTtb;out ) +_Qtb;(9) CshdTsh;out dt= _mshcp;sh(Tsh;inTsh;out ) +_Qsh; (10) where heat capacities are de ned by eqs. 11 and 12. Masses depend on volumes, which as constant values, and densities which depend on the mean temperature of the CVs. Speci c heat capacities also depend on the mean temperature of CVs. Ctb=mtbcp;tb; (11) Csh=mshcp;sh: (12) Heat ow rates are de ned by eqs. 13 and 14, _Qtb= heAw(Tsh;outTtb;out ); (13) _Qsh= heAw(Ttb;outTsh;out )_Qloss; (14) where the overall heat transfer coecient ( he) can be calculated by eq. 6 and thermal losses by eq. 15. _Qloss= lossAamb(Tsh;outTamb): (15) Heat transfer coecients ( sh, tband loss) are values that need to be calibrate with experimental data. 3.3. Dynamic with Thermal Capacitance Model (DYTC model) The Dynamic with thermal capacitance (DYTC) model is an extension of the DY model consider- ing the heat capacity of the tube bundle metallic parts, therefore assuming that the wall temperature is the average value between the tube-side and shell- side uid temperatures in each cell ( Tw= (Tsh;out +Ttb;out )=2), the tube-side and shell-side heat capac- ities can be then reformulated according to eqs. 16 and 17. Ctb=mtbcp;tb+1 2mwcp;w; (16) Csh=mshcp;sh+1 2mwcp;w; (17) where the masses and speci c heat capacities are parameters of the model. 3.4. Dynamic with Thermal Capacitance and Cell Method Model (DYCM model) The Dynamic with thermal capacitance and cell method (DYCM) model has been developed following an object-oriented methodology based on rst prin- ciples. A preliminary version of this model was pre- sented in [10], some extensions have been introduced with respect to this preliminary version. One of them is heat losses to the environment and the others are more detailed shell and tube models considering the cell method introduced in [18]. The Modelica icon of this model is shown in g- ure 3, whereas the component diagram can be seen in gure 6. This diagram shows that each component has been modeled independently (shell, tube bundle, tube bundle walls and convection processes). The model considers the uids ow in the tube side, as well as, in the shell side and convective heat trans- fer processes between the tube-side uid, the tube bundle walls, the shell-side uid and the ambient. Exchange of information between components is per- formed by means of connectors. Round connectors exchange information about uids (mass ow rate, pressure, speci c enthalpy and mass fraction, the last one is optional), whereas square connectors exchange information about heat (temperature and heat ow rate). The main components of the model are de- scribed in the following paragraphs. 3.4.1. Convection model. Newton's law of cooling has been considered to model the convective heat transfer. The heat transfer coecient can be dynamically calculated from any 7 system gdefaultsShellShell-side inlet Shell-side outletTube-side outlet Tube-side inletTube bundle walls Tube bundleShell-side outletShell-side inlet Tube bank temperatureHeatportconvection_shellSideFluid_tubeWalls Tube-side inletTube-side outlet convection_tubeWalls_tubeSideFluidconvection_shellSideFluid_shellWalls b_shella_shell b_tube a_tubeAmbientFigure 6: Heat exchanger Modelica component diagram correlation implemented in the model. Convection model components can be seen in gure 6: heat transfer from the tube-side uid to inner tube walls (convection tubeWalls tubeSideFluid component), heat transfer from outer tube walls to the shell-side uid ( convection shellSideFluid tubeWalls compo- nent) and heat transfer from shell-side uid to the environment ( convection shellSideFluid shellWalls component). Heat transfer coecients ( sh, tband loss) are values that must be calibrated for the particular heat exchanger under consideration. 3.4.2. Tube bundle wall model. This model is shown in gure 6 as the tube bundle walls component. It is discretized in CVs, the energy balance is modeled in each one of them according to eq. 18. The number of parallel tubes is a parameter of the model. Density and speci c heat capacity val- ues can be con gured as average constant values or interpolated values as function of wall temperature. mwcp;wdTw dt=_Qw: (18) Figure 7: TEMA E-type shell heat exchanger (lateral view) 3.4.3. Shell model. The shell component in gure 6 belongs to this class. Shell volume and dimensions have been mod- eled according to manufacturer data (see table 2). The ow of the shell-side uid has been discretized in CVs where the one-dimensional dynamic mass (eq. 19), dynamic energy (eq. 20) and static mo- mentum (eq. 21) balance equations have been con- sidered. These equations apply for each CV in the model. Mass and energy balance equations have been discretized according to the Finite Volume Method (FVM) [19], whereas the static momentum balance equation is lumped in a global equation for the whole shell. Therefore, pressure drop is equally distributed among the CVs in the shell model. Pressure drop is calculated according to [20, 21], where pressure drop at cross- ow and window sections, as well as pressure drop at the inlet and outlet nozzles, for shell-and-tube heat exchangers with segmental baes, are taken into account. dmsh dt= _msh;in_msh;out; (19) dUsh dt=_Hsh;in_Hsh;out +_Qsh; (20) psh=psh;inpsh;out: (21) This model has been extended with respect to 8 Figure 8: Cell method in the shell side of a E-type shell heat exchanger (lateral view) Figure 9: Cell method in the shell side of a F-type shell heat exchanger (upper view) the one presented in [10] by considering the cell method introduced in [18]. Such work focused in heat exchangers with TEMA E-type shell con gurations (one shell pass). Figure 7 shows a lateral view of a E-type shell heat exchanger with 6 vertical segmental baes and 2 tube passes. The shell is divided in cells to account for the di erent ow patterns by means of linking the cells to represent the ow of the shell- side uid. The number of shell cells is the number of baes plus one multiplied by the number of tube passes, see gure 8. The same idea can be applied to F-type shells, as presented in [2]. Therefore, this approach was consid- ered in the present model. In [2], the cell method was applied to each shell pass. Note that there are 2 shellpasses due to the longitudinal bae (see gure 3). Additionally, two di erent kinds of cells with di erent ow conditions where de ned in that work: window and cross- ow cells. In window cells, the shell uid ows in counter ow with respect to the tube-side uid, whereas it ows in cross ow in the other kind of cells. Figure 9 shows the cell distribution, where window and cross- ows cells can be seen. Window cells have a diagonal line pattern background. Note that this is an upper view of the heat exchanger and only one shell pass is represented in this gure. 3.4.4. Tube bundle model. The tube bundle component in gure 6 belongs to this class. The ow of the uid in a single tube of the tube bundle has been modeled as in the shell model, i.e. considering the one-dimensional dynamic mass (eq. 19), dynamic energy (eq. 20) and static momen- tum (eq. 21) balance equations. Only one tube is modeled per CV or cell, and then mass and heat ow rates are scaled considering the number of parallel tubes in each particular cell. Pressure drop in the tube bundle has been mod- eled considering all the possible cases in the Moody diagram [22] using a model available in the Modelica Fluid library [13]. Pressure drop at nozzles has been calculated according to [21], whereas pressure drop in the U-shaped elbows turned out to be small and has been neglected. As in the shell model, the tube bundle model has been extended with respect to the one presented in [10] by considering the cell method. In [18], the cell method was applied to a E-type shell heat exchanger (many tube passes). Figure 10 shows a lateral view of a E-type shell heat exchanger with two tube passes, where the tube bundle is divided in cells in order to be linked to shell cells (see gure 8) and to establish heat exchange between them. In [2], the same method was applied to a F-type shell heat exchanger and this variant was applied in this model. Figure 11 shows the upper view of a F- type shell heat exchanger where only one tube pass is represented. Segmental vertical baes are included in Figure 11 although they do not in uence the tube- side uid ows. Three di erent ow cells are de ned 9 Figure 10: Cell method in the tube side of a E-type shell heat exchanger (lateral view) Figure 11: Cell method in the tube side of a F-type shell heat exchanger (upper view) (A,B and C cells) in parallel. The cell numbers de- note the ow order in this representation. In each one of the ow cells, one tube is considered and therefore heat ow rates are scaled as a function of the num- ber of parallel tubes in the tube bundle. Additionally, the tube-side inlet and outlet channels in each unit are taken into account in the model (D cells). They dispatch or gather the tube-side uid to or from the tubes in the tube bundle. In [18], the number of cells was the number of baf- es plus one multiply by the number of tube passes for a TEMA E-type shell heat exchanger. Consid- ering an object-oriented modeling approach it is not required that the number of shell cells would be the same as the number of tube cells, since using heatports multiple cells can be linked. This is the case for shell cells 1 and 13 (see gure 9), shell cell 1 is connected with tube cells A4 and B7, and shell cell 13 is connected with tube cells B1 and C1 (see g- ure 11). For F-type shell heat exchangers, the number of shell and tube cells are given by eqs. 22 and 23, re- spectively, where nbis the number of baes and npa the number of shell or tube passes which in a F-type shell are the same. nsh=npa(2nb+ 1); (22) ntb=npa(2nb+ 3): (23) Therefore, in our case that could make 158 shell CVs and 162 tube CVs, since the studied heat ex- changer has 39 baes per unit with two passes per unit. If we consider the remaining components in g- ure 6 (tube bundle walls and convection components) and their 158 CVs per component, that could make a total of 952 CVs. In order to reduce the time required for the simulation, the number of CVs has been re- duced. Section 4.1 shows an study on the impact of the reduction in the number of CVs with respect to the simulation results. 4. Simulation, Calibration and Validation Against Experimental Data This section concerns with the simulation of the di erent heat exchanger models. Section 4.1 stud- ies the in uence of the number of CVs with respect to simulation results in the DYCM model. In or- der to verify the models, design performance cal- culations provided by the manufacturer were com- pared to simulation results at design conditions in section 4.2. Thermocouple measurement uncertain- ties are presented in section 4.3. After that, some values and parameters of the models were calibrated in section 4.4 according to data from experimental campaigns. Once the models were calibrated, simula- tions results were validated against experimental data in section 4.5 in steady state and at several transient responses to check the model performances under a wide range of operating conditions. 10 Dymola [23] is the tool used for the Modelica implementations and simulations. The numerical solver used for the dynamic simulations has been DASSL [24], a variable-step variable-order multistep solver in the numerical integration of di erential alge- braic equation systems, where the Newton iteration method is used to solve the resulting nonlinear al- gebraic system of equations at each time step. The absolute and relative tolerances were set to 104. In order to perform the heat exchanger simulations, the following inputs must be provided to the models, as shown in gure 12. •Tube-side and shell-side inlet mass ow rates ( _mtb;in, _msh;in). •Tube-side and shell-side inlet uid temperatures (Ttb;in,Tsh;in). •Tube-side and shell-side outlet pressures ( _ ptb;out , _psh;out ). •Ambient temperature ( Tamb). 4.1. Level of Discretization In this section, the DYCM model is simulated con- sidering di erent numbers of CVs in order to study its convergence. Nominal operating conditions have been taken into account in the simulation according to table 1. The mean overall heat transfer coecient has been set according to the heat exchanger design parameters (see table 2). The tube-side and shell-side heat transfer coecients have been calculated consid- ering eq. 6 and assuming that both coecients have the same value. With the aim of studying the convergence not only in steady-state but also in transients, the initial tem- perature of the whole heat exchanger was set to 290C. Table 4 shows molten salt and thermal oil outlet temperatures in steady state, di erences with respect to the most detailed model (158 CVs) and execution times (compilation plus simulation times) for di erent cases. The number of CVs refers to the shell component. From these results a good compromise between ac- curacy and execution time would be 50 CVs since the TInMS mﬂowInMSm source_molten_salts m_ﬂow sink_molten_saltspOutMS m source_thermal_oilm_ﬂow mﬂowInOil TInOilsink_thermal_oilp pOutOil Tube-side inlet nozzleTube-side outlet nozzle Shell-side inlet nozzle Shell-side outlet nozzleKAmbientTemperature tAmbT p TTsh,in ptb,out Tambmsh,in psh,out Ttb,inmtb,in..Figure 12: Heat exchanger model inputs error committed is negligible (0 :12C) with respect to the most detailed model while being close to four times faster. Execution times are low in this simula- tion since it is a simple step response. However, when considering more complex simulations, the simulation time is a factor than must be taken into account. 4.2. Design Conditions All the models at nominal operating conditions (see table 1) have been compared in this section in the same terms as in the simulation performed in sec- tion 4.1. The DYCM model has been discretized in 50 CVs, whereas the DY and DYTC models consider 160 CVs, according to the formula given in [18] (see section 3.4.4) for the given number of baes, tubes passes and units. Figures 13 and 14 show the molten salt and thermal oil outlet temperatures for the di er- ent models and also includes the nominal operating conditions provided by the manufacturer. Table 5 11 Table 4: Level of discretization CVsTemp. (C) Di . (C) Execution MS Oil MS Oil time (s) 158 370:20 315:18 - - 40:10 110 370:18 315:19 0:02 0:01 22:12 90 370:15 315:19 0:05 0:01 18:46 50 370:08 315:19 0:12 0:01 10:73 30 369:97 315:22 0:23 0:04 7:34 10 369:20 315:65 1:00 0:47 4:70 shows molten salt and thermal oil outlet tempera- tures in steady state and execution times for all the models plus manufacturer design conditions. The simple algebraic QSS model provides a good approximation. The DY and DYTC models pre- dict higher molten salt outlet temperatures than the DYCM model since they assume a pure counter- ow heat exchanger whereas the DYCM model considers two di erent ow sections: cross- ow and window sections, and a more precise approximation in theory with respect to the remaining models and the manu- facturer calculations. The e ect of the tube bundle thermal capacitance in the DYTC model with respect to the DY model can be clearly seen in Figure 14 in the initial tem- perature change from 290C to nominal operating temperature. Table 5: Design conditions ModelTemperature (C) Execution MS Oil time (s) Manufacturer 373:00 313:00 - QSS 374:18 315:06 2:60 DY 373:20 312:58 16:67 DYTC 373:20 312:57 16:70 DYCM 370:08 315:19 10:73 4.3. Thermocouple meausurement uncertainties According to the International Electrotechnical Commission (IEC) 584.3 norm, the allowable man- ufacturing tolerance of the K-type class 2 thermo- couples is up to 3C. However, the thermocouples 0 500 1000 1500 2000 2500 3000 3500 4000 Simulation time (s)290300310320330340350360370380Molten salt outlet temperature ( °C) Manufacturer QSS DY DYTC DYCMFigure 13: Molten salt outlet temperature from di erent mod- els at design conditions 0 500 1000 1500 2000 2500 3000 3500 4000 Simulation time (s)290295300305310315320Thermal oil outlet temperature ( °C) Manufacturer QSS DY DYTC DYCM Figure 14: Thermal oil outlet temperature from di erent mod- els at design conditions in the facility were calibrated with respect to a digi- tal thermometer in a temperature dry-well calibrator. With the aim of tting thermocouple temperatures to the digital thermometer temperatures, polynomial functions are used. Therefore, the accuracy of the thermocouples is improved, the combined standard uncertainty at 400C involve the following uncertain- ties, at lower temperatures they can be slightly lower. 4.4. Calibration Several heat transfer correlations, besides the mean overall heat transfer coecient provided by the man- ufacturer, have been implemented in the models and compared against data from experimental campaigns. In the shell side: Gaddis and Gnielinski [21], the Bell- Delaware method [25] and a correlation proposed in [26] which is a curve t from data provided in [27], whereas in the tube side: [28], [29] and [30] correla- tions have been also tested. 12 However, simulation results did not agree with ex- perimental data. This is because there is a perfor- mance detriment in the heat exchanger [10]. In order to identify the origin of the performance detriment and to develop an accurate heat exchanger model, thermal losses as well as heat transfer correlations on both uid sides were calibrated against experimental data. The applied methodology is explained in detail in [31]. Heat transfer coecients can be calculated from the Nusselts number by means of eq. 24, =NuK D: (24) A commonly used expression to calculate the heat transfer coecient in the shell side of heat exchang- ers, and analogous to eq. 24, is given by eq. 25, =jcpG Pr2=3: (25) 4.4.1. Thermal losses. As a result of the calibration process using exper- imental data as described in [31], the following Col- burn j factor for thermal losses, eq. 26, was obtained. Therefore, the heat transfer coecient for thermal losses is calculated by means of eq. 25 and eq. 26. This heat transfer correlation can be used in all the developed models. jloss= 1:1858Re0:9545 loss: (26) 4.4.2. Tube-side heat transfer correlation. The QSS, DY and DYTC models do not model the tube-side inlet and outlet channels in each heat ex- changer unit, for this reason correlations in the liter- ature cannot be directly used. For instance, Gnielin- ski's correlation [28] (see eq. 28) was derived con- sidering uid ow in straight ducts. Although this correlation is a good approximation for the tube side of heat exchangers, the coecients appearing on it (1000, 12.7) can be adjusted experimentally, since uid ow path in heat exchangers is commonly com- plex. Therefore, those coecients where calibrated in eq. 27 obtaining (1792, 29.93). The calibrationwas performed using experimental data and the cal- ibration process is described in [31], as a result the tube-side heat transfer coecient in the QSS, DY and DYTC models can be calculated by means of eq. 24 and eq. 27. Nutb=fw 8(Retb1792)Prtb 1 + 29:93p fw=8(Pr2=3 tb1)" 1 +dw lw2=3# : (27) Since the DYCM model has into account the tube- side inlet and outlet channels in each heat exchanger unit, the calibration of the previously mentioned co- ecients yields similar values to the coecients in the original Gnielinski's correlation. For this reason, the original correlation has been used, eq. 28. Nutb=fw 8(Retb1000)Prtb 1 + 12:7p fw=8(Pr2=3 tb1)" 1 +dw lw2=3# : (28) The friction factor ( fw) in both equations has been calculated considering the Filonenko's correla- tion [32], eq. 29, fw= (1:82 lnRetb1:64)2: (29) 4.4.3. Shell-side heat transfer correlation. Shell-side heat transfer correlation was also cali- brated for the QSS, DY and DYTC models as de- scribed in [31], therefore the shell-side heat transfer coecient can be calculated by means of eq. 25 and eq. 30. jsh= 3:2470Re1:1077 sh; (30) Even though the DYCM has a di erent and more complex shell-side model structure, the calibrated co- ecients, eq. 31, are similar to those obtained for the other models, eq. 30. jsh= 3:2609Re1:1074 sh: (31) 13 Nevertheless, the calibrated shell-side correlations (eq. 30 and eq. 31) do not agree with correlations in the literature: Gaddis and Gnielinski [21], the Bell- Delaware method [25] and [26], therefore this suggests that the performance detriment must be related to the heat transfer in the shell side. The most common causes for deterioration in per- formance of F-shell heat exchangers are thermal leak- age or physical leakage due to the longitudinal baf- e [33] together with fouling, corrosion, design errors and fabrication issues. Additionally, two potential is- sues were identi ed with this heat exchanger, as pre- sented in [9]. One of them is the bypass of molten salt through the drainage channels and the other one is the nitrogen accumulation inside the shell due to the heat exchanger tilt angle. Further investigation and inspection in the physical device is necessary to address the performance issue. 4.5. Experimental validation In this section data from experimental campaigns is compared against simulation results from the mod- els previously introduced in section 3, three experi- ments are considered. One of them in steady-state conditions and the two remaining considering tran- sients: replication of cloud disturbances in the solar eld and molten salts, as well as thermal oil, mass ow rate steps. Again, the DYCM model has been discretized in 50 CVs, whereas the DY and DYTC models consider 160 CVs. 4.5.1. Steady-state experiment. A steady-state experiment was performed in the MOSA facility close to nominal operating conditions (see table 1). Figure 15 shows the inlet molten salt and thermal oil mass ow rates in the experiment. All the models were simulated considering the cal- ibrated heat transfer correlations (see section 4.4), for comparison purposes the DYCM model was sim- ulated also considering the overall heat transfer coef- cient provided by the manufacturer (DYCM manu- facturer), as well as a shell-side heat transfer correla- tion available in the literature, Gaddis and Gnielin- ski's correlation [21] (DYCM correlation). Experimental inlet and simulated outlet molten salt and thermal oil temperatures can be seen in 11:30 12:00 12:30 13:00 13:30 14:001.31.41.51.61.71.81.922.1 Local timeMass flow rate (kg/s) Molten salt Thermal oilFigure 15: Steady-state experiment: mass ow rates 11:30 12:00 12:30 13:00 13:30 14:00 Local time290300310320330340350360370Molten salt temperature ( °C) Exp. inlet Exp. outlet QSS calibrated DY calibrated DYTC calibrated DYCM calibrated DYCM manufacturer DYCM correlation Figure 16: Steady-state experiment: molten salt temperatures 11:30 12:00 12:30 13:00 13:30 14:00 Local time300310320330340350360370380Thermal oil temperature ( °C)Exp. inlet Exp. outlet QSS calibrated DY calibrated DYTC calibrated DYCM calibrated DYCM manufacturer DYCM correlation Figure 17: Steady-state experiment: thermal oil temperatures gures 16 and 17, respectively. Table 6 shows the maximum di erence between experimental and sim- ulated thermal oil and molten salt temperatures to- gether with the Root Mean Square Error (RMSE) and the execution time for each model. Simulation re- 14 Table 6: Simulation results in the steady-sate experiment ModelMolten salt Thermal oil Execution Maximum di erence RMSE Maximum di erence RMSE time QSS calibrated 5:51C 3:87C 6:17C 4:22C 3:06 s DY calibrated 3:21C 1:77C 7:25C 2:09C 12:09 s DYTC calibrated 3:37C 1:80C 7:25C 2:14C 11:60 s DYCM calibrated 1:22C 0:54C 3:04C 1:78C 54:45 s DYCM manufacturer 20:82C 19:56C 15:75C 12:81C 20:03 s DYCM correlation 22:00C 20:73C 17:02C 13:91C 59:57 s Table 7: Simulation results in the cloud disturbance experiment ModelMolten salt Thermal oil Execution Maximum di erence RMSE Maximum di erence RMSE time QSS 13:58C 4:30C 14:14C 6:26C 2:94 s DY 6:70C 2:68C 12:01C 3:76C 9:90 s DYTC 4:28C 2:61C 10:59C 3:35C 9:61 s DYCM 4:29C 1:51C 9:85C 2:68C 66:07 s sults from calibrated models agree with experimental data, being the DYCM results the most precise but at expense of a higher computation time since the DYCM model is more complex than the other mod- els. It can be seen in gures 16 and 17 that for the DYCM manufacturer and DYCM correlation models, the simulated performance is much higher than the real one, this again suggests issues in the shell-side heat transfer of the heat exchanger, as commented in section 4.4. 4.5.2. Cloud disturbance experiment. Clouds in the solar elds were replicated in an op- erating mode 3 experiment (see section 2) by reducing the oil heater temperature and then setting it back to its original value, as it can be seen in the ther- mal oil inlet temperature in gure 19. This gure also shows the experimental and simulated thermal oil outlet temperatures given by each model. Ther- mal oil and molten salt mass ow rates were kept constant to their nominal values during the whole experiment. It can be seen in gure 19 that the sim- ulated thermal oil outlet temperature given by the models react faster to the disturbance that the realsystem, this must be further studied and it might be related to unmodeled dynamics or to issues in the thermal oil outlet temperature measurement. 11:30 12:00 12:30 13:00 13:30 14:00 Local time290300310320330340350360Molten salt temperature ( °C) Exp. inlet Exp. outlet QSS DY DYTC DYCM Figure 18: Cloud experiment: molten salt temperatures On the other hand, gure 18 shows the experimen- tal inlet and outlet as well as simulated outlet molten salt temperatures given by each model. In this case, all the models provides a good agreement with re- spect to experimental data. This fact reinforces the idea that the disagreement between the reaction time of the real system and the models with respect to the thermal oil outlet temperature should be related to 15 11:30 12:00 12:30 13:00 13:30 14:00 Local time300310320330340350360370380Thermal oil temperature ( °C) Exp. inlet Exp. outlet QSS DY DYTC DYCMFigure 19: Cloud experiment: thermal oil temperatures issues in the thermocouple, thermowell, or in the pip- ing between the outlet of the heat exchanger and the thermocouple, because otherwise the outlet molten salt dynamics would not show such a good agreement. Table 7 shows the same information than Table 6 but for the cloud disturbance experiment. It can be seen that the DYCM model provides slightly better results than the DYTC model but at the expense of a much higher computation time. 4.5.3. Mass ow rate experiment. Some thermal oil and molten salt mass ow rate steps were applied in this experiment as shown in gure 20. Same information as for the previous ex- periment is also shown in this case, where gure 21 shows the experimental inlet and outlet as well as simulated outlet molten salt temperatures for each model, gure 22 shows same variables for thermal oil and table 8 shows the maximum di erence between experimental and simulated thermal oil and molten salt temperature, the RMSE and execution times. There is a good agreement between experimental and simulated outlet thermal oil and outlet molten salt temperatures, as shown in table 8, being the RMSE lower than 3 :50C for all the evaluated mod- els. The maximum di erence between experimental and simulated temperatures is mainly due to the ini- tialization process, where the DYCM model provides the best initialization as also shown in table 8. In this experiment, it was found out that the mod- els slightly underpredict heat transfer when the ther- mal oil mass ow rate is higher than its nominal value while keeping the molten salt mass ow rate at nom- 11:30 12:00 12:30 13:00 13:30 14:00 14:30 Local time1.61.822.22.42.62.8Mass flow rate (kg/s)Molten salt Thermal oilFigure 20: Mass ow rate steps experiment: mass ow rates 11:30 12:00 12:30 13:00 13:30 14:00 14:30 Local time300310320330340350360370Molten salt temperature ( °C) Exp. inlet Exp. outlet QSS DY DYTC DYCM Figure 21: Mass ow rate steps experiment: molten salt tem- peratures inal conditions. This is depicted in gure 20, at ap- proximately 13:15, where the thermal oil mass ow rate is 2:02 kg=s, being its nominal value 1 :57 kg=s (see table 1), and keeping the molten salt mass ow rate at its nominal value, 2 :08 kg=s. Figure 22 shows that heat transfer is slightly underpredicted, since the simulated thermal oil outlet temperatures are higher that the experimental one. In gure 21, this can be also appreciated since simulated molten salt temper- atures are lower that the experimental one. Never- theless, the discrepancy between experimental and simulated temperatures is still low in this part of the experiment, the highest one is found in the ther- mal oil temperature, being this lower than 2 :90C for the DY, DYTC and DYCM models, and lower than 3:60C for the QSS model. Again the DYCM model provides slightly better results at the expense of a much higher computation time (see table 8). 16 Table 8: Simulation results in the mass ow rate steps experiment ModelMolten salt Thermal oil Execution Maximum di erence RMSE Maximum di erence RMSE time QSS 12:33C 3:42C 12:35C 3:44C 2:94 s DY 9:35C 1:74C 17:02C 2:39C 12:96 s DYTC 9:35C 1:72C 17:02C 2:36C 13:22 s DYCM 3:15C 1:59C 7:68C 2:27C 72:30 s 11:30 12:00 12:30 13:00 13:30 14:00 14:30 Local time320330340350360370380Thermal oil temperature ( °C) Exp. inlet Exp. outlet QSS DY DYTC DYCM Figure 22: Mass ow rate steps experiment: thermal oil tem- peratures 5. Conclusions and Ongoing Work This paper has presented a comparison of di er- ent shell-and-tube heat exchanger models with dif- ferent degree of complexity. Simulations results have been compared against experimental data in terms of accuracy and simulation time. The following para- graphs summarize the main conclusions drawn from this study. With respect to the tube-side modeling, it is im- portant that the heat exchanger model takes into account the inlet and outlet channels connected to the tube bundle in order to obtain accurate results if heat transfer correlations in the literature for the tube bundle are planning to be used, as it is the case for the DYCM model where the Gnielinski's correla- tion was considered. Another option is to use general heat exchanger models in the literature and calibrate the tube-side heat transfer correlation, this is the case for the QSS, DY and DYTC models considered in this work, where the tube-side inlet and outlet channels are not taken into account.All the models helped to identify that there is a performance detriment in the heat exchanger studied, but the most complex model, the DYCM model, sug- gests that the performance detriment is related to the heat transfer in the shell side, i.e between the shell- side uid and the tube bundle, since heat transfer correlations in the literature can be used in the tube- side and simulation results agree with experimental data, whereas the shell-side heat transfer correlation had to be calibrated. According to the results, the QSS is a fast and simple algebraic model which may be used in ad- vanced control techniques which require fast and mul- tiple model evaluations, for instance model predictive control. The DYCM model provides slightly better results with respect to experimental data, however the computation time to perform the simulations is much higher than with the other models; if this sim- ulation time can be assumed, this is the most pre- cise model for process simulations as well as to test control strategies in simulation. The DYTC model provides better results than the QSS and DY mod- els and results are close to those obtained with the DYCM model but in lower computation times. Ongoing work includes inspecting the real system in order to identify the physical causes for the heat exchanger performance detriment. Additionally, the QSS and DYCM models are being used for the de- sign, testing, validation and optimization of control strategies in heat exchangers for TES systems in so- lar thermal power plants [34]. In particular, the QSS model is being used in a feed-forward model predic- tive control and the DYCM model as the dynamic model of the real system. 17 Acknowledgements This research has been funded by the EU 7th Framework Programme (Theme Energy 2012.2.5.2) under grant agreement 308912 - HYSOL project - Innovative Con guration of a Fully Renewable Hy- brid CSP Plant, the Spanish Ministry of Economy and Competitiveness through ERDF, PLAN E funds (C.N. SolarNOVA ICT-CEPU 2009-02) and the Na- tional R+D+i Plan Project DPI2014-56364-C2-1/2- R (ENERPRO-EFFERDESAL).