A Grey-Box Dynamic Model of Plate Heat Exchangers Used in an Urban Heating System

This paper deals with the establishment of a grey-box model of plate heat exchangers. First, a second-order state space model of plate heat exchangers is developed. Then, the model is solved by the white-box method and grey-box method, respectively. In the white-box method, the linearization of the state space model is carried out by implicit difference scheme and the overall heat transfer coefficient was calculated with empirical correlation of the Nusselt number. In the grey-box method, a newly developed parameter identification method was established. The simulation results of two outlet temperatures by the grey-box and white-box method, respectively, are compared with the test data. It is shown that the grey-box method is more suitable and accurate for the simulations of outlet temperatures than the white-box method. Thus, it can provide practitioners and heating engineers with tools to design controllers for district heating systems.


Introduction
Heat exchangers are commonly used in many fields such as power, district heating, heat recovery, industrial process cooling and heating, hydraulic oil cooling or cooling of machine tools.Plate heat exchangers (PHEs) are probably the most common type of heat exchangers due to their high heat transfer coefficients, compactness, convenience for cleaning, and flexibility.
Due to the importance of PHEs, many researchers have investigated modeling and the optimum design parameters for the system.In Gut and Pinto's work [1], a mathematical model was developed in an algorithmic form for the steady-state simulation of gasketed PHEs with generalized configurations.The configuration influence on the exchanger performance was studied and a method for configuration optimization was developed.Gut et al. [2] presented a parameter estimation procedure for PHEs that handles experimental data from multiple configurations.The heat transfer correlations obtained for PHEs are intimately associated with the configurations experimentally tested and the corresponding flow distribution patterns.Gut and Pinto [3] proposed a screening method to find the optimal set of configurations of a given PHE.Mota et al. [4] proposed an algorithm based on the screening method, for the optimization of heat exchange area of PHEs.In the work of Zhang et al. [5], a general three-dimensional distributed parameter model (DPM) was developed for optimal design of the plate-fin heat exchanger.A three-dimensional (3D) numerical model of the plate-type heat exchanger was established, and the flow pattern and heat transfer effect of the cold and hot fluids were analyzed by using ANSYS (ANSYS, Inc., Energies 2017, 10, 1398 2 of 16 Canonsburg, PA, USA) [6].It can effectively solve the problems of optimization design of a plate-type heat exchanger.Yang et al. [7] proposed a mathematical model to quantitatively evaluate the effect of flow maldistribution in a multi-channel heat exchanger by computational fluid dynamics (CFD) technology.
In district heating system, PHEs and their networks frequently undergo transients resulting from heat load variations and regulations.Consequently, it is very important to know the behavior of PHEs when they are subjected to transient flow.
Luo et al. [8] have developed a new linearized one-dimensional model to solve non-linear problems with small non-linear disturbances.Dwivedi and Das [9] presented a predictive model to suggest the transient response of PHEs, subjected to a step flow variation.They indicated that the theoretical model developed for flow transient is capable of predicting the transient performance of the PHEs satisfactorily.Laszczyk [10] proposed a simplified cross-convection model for a wide class of heat exchangers.The simplified tuning procedure of the model was suggested and verified based on experimental data.The validation results show that the suggested modeling and tuning method is useful for practical applications.In Turhan et al.'s work [11], a data-driven grey-box identification approach provided the physically meaningful nonlinear continuous-time model.The validated results by the Monte Carlo tool serve to demonstrate the effectiveness and the viable application of the suggested control solutions to air conditioning systems.
In the general case, the convective heat transfer coefficient was calculated with empirical correlations of the Nusselt number in the above literature.The Nusselt number is a function of the Reynolds number and the Prandtl number.The influence of the Reynolds and the Prandtl numbers on heat transfer in a wide range of their values cannot be accounted for by simple power equations with a constant using the Prandtl number exponent, which is usually used in empirical correlations.For this reason, a semi-empirical equation was proposed and the Prandtl number's influence on heat transfer was discussed in [12].
In addition to the convective heat transfer coefficient, fouling resistance also impacts the overall heat transfer coefficient of PHE.Fouling of PHE leads to a rise in pressure drop across the exchanger and the decrease of the overall heat transfer coefficient.Jonsson et al. [13] proposed a method which used nonlinear physical state space model for detection of the fouling of heat exchangers in transient states.The model parameters were estimated by using an extended Kalman filter and measurements of inlet and outlet temperatures and mass flow rates.Mahdi et al. [14] proposed a two-dimensional dynamic model for milk fouling in a PHE and stated that the aggregation rate of unfolded protein is found to increase exponentially with wall temperatures rising and can be accompanied by a substantial reduction in the heat-transfer coefficient.Grijspeerdt et al. [15] carried out two-dimensional (2D) and 3D CFD simulations of flow between two corrugated plates to investigate the flow distribution for milk processing with respect to fouling.
Optimal operation and real-time control require a more accurate description of the time-domain behavior of heat exchangers.Therefore, it is necessary to develop a mathematical model of a PHE, which is easy to parameterize and accurate enough to capture the essential dynamics and nonlinearities.A simplified and efficient grey-box method is established in this paper, which can avoid the calculation of the overall heat transfer coefficient.

Modeling of Plate Heat Exchanger
In order to derive the mathematical model of the PHE, the following assumptions are considered.
(1) The fluid is assumed to be incompressible.
(2) Due to the short width of the heat transfer channels, the flow rate component perpendicular to the flow direction of the plate approximates to 0, and the temperature distribution in this direction is neglected.(3) Heat loss to the surroundings is ignored.(4) The fluid velocity in the main flow direction is high enough to justify that heat conduction along the flow direction can be neglected as compared to heat convection.
Energies 2017, 10, 1398 3 of 16 (5) No phase change of the fluid occurs inside the heat exchanger.(6) The end effects of the first and last channels are neglected.(7) The thermal capacity of the plates is neglected.
The fluid inside a channel exchanges heat with the neighbor channels through the plates, as shown in Figure 1.The arrangement of flow channels in each side is parallel.According to Assumption (2), the temperatures of cold and hot fluid in their respective channels are constant over the cross-section of channels.The control volumes dV h (in the hot water channel) and dV c (in the cold water channel) are shown in Figure 1, where, Energies 2017, 10, 1398 3 of 17 (6) The end effects of the first and last channels are neglected.(7) The thermal capacity of the plates is neglected.
The fluid inside a channel exchanges heat with the neighbor channels through the plates, as shown in Figure 1.The arrangement of flow channels in each side is parallel.According to Assumption (2), the temperatures of cold and hot fluid in their respective channels are constant over the cross-section of channels.The control volumes Applying the thermal balance to the control volumes shown in Figure 1, it is possible to derive the partial differential equations for the channel temperatures h T and c T .
( ) ( )  1) and ( 2) can be simplified as: ( ) ( ) Applying the thermal balance to the control volumes shown in Figure 1, it is possible to derive the partial differential equations for the channel temperatures T h and T c .
where ρ w is the fluid density (kg/m 3 ); c p is specific heat (J/(kg•K)); δ p is the distance between two neighboring plates (m); b p is the effective channel width which is obtained with the effective plate exchange area dividing by the length of the flow channel (m); k denotes the overall heat transfer coefficient between the hot and cold fluid (W/(m 2 •K)); T c and T h are respectively the channel temperatures of cold and hot fluid (K); and u c and u h are the flow velocity of cold and hot fluid (m/s), respectively.Equations ( 1) and ( 2) can be simplified as: The overall heat transfer coefficient k is obtained from Equation ( 5), which represents a series association of thermal resistances.
where h h and h c denote the hot fluid side and cold fluid side convective heat transfer coefficient, respectively (W/(m 2 •K)); H p is the plate thickness (m); λ p is the plate thermal conductivity (W/(m•K)); R f ,h and R f ,c are the hot and cold fluid fouling factor, respectively (m 2 •K/W); d e = 2δ p is the hydraulic diameter (m); and λ w is the fluid thermal conductivity (W/(m•K)).
As a general rule, the calculation of the convective heat transfer coefficient is required for an accurate knowledge of the fluid velocity field which in turn leads to high-order models with high computational effort.Therefore, the approach shown below relies on a suitable approximation of the functional dependence of the Nusselt number Nu, on the Prandtl number Pr and the Reynolds number Re.For Newtonian fluids under turbulent flow, the correlation whose general form is presented in Equation ( 7) is commonly used, where c Nu , c Re , and c Pr are empirical parameters [16].
The fluid parameters λ w , µ, ρ, c p are evaluated at the average of inlet and outlet temperatures.
Integrating Equations ( 3) and ( 4) over the length of the channel yields the partial differential equations where l p is the effective channel length (m).According to Assumption (3), the inlet temperatures of channels at the same side are identical.In addition, fluid flows and geometrical parameters of each channel of cold fluid and hot fluid are the same.Therefore, the overall heat transfer coefficient k is constant along the flow direction under the given geometrical parameters.
In the next step, Equations ( 9) and ( 10) are discretized spatially to obtain the ordinary differential equations (odes).To acquire a low-order model of the PHE, the fluid channel is considered as a whole control volume.Time derivatives of the temperature are assumed to be uniform inside the control volumes.Therefore the integral for time derivative terms in Equations ( 9) and ( 10) can be written as: The spatial derivative terms in Equations ( 9) and ( 10) can be derived as: Energies 2017, 10, 1398 5 of 16 where T c,1 and T c,2 are respectively the inlet and outlet channel temperatures of cold fluid (K); and T h,1 and T h,2 are respectively the inlet and outlet channel temperatures of hot fluid (K).The last terms in Equations ( 9) and ( 10) denote the transient heat transfer between the two fluids, where the integral can be simplified via the trapezoidal rule.This simplification relies on a small temperature deviation between the inlet and outlet temperatures, together with a linear temperature distribution along the channel.Thus, the last integral terms can be simplified as follows: Substituting Equations ( 11)-( 16) into Equations ( 9) and ( 10), the simplified ordinary differential equations for PHE are as follows: Equations ( 17) and ( 18) can be written as the state space form: . T h,2 .

Methods
In order to solve Equation (19), the overall heat transfer coefficient k must be obtained.Generally, it is calculated according to Equation (5); this method called the white-box method.

White-Box Method
Applying the implicit difference scheme, the state-space equation Equation ( 19) can be expressed as: where the superscript "j" denotes the jth instant, and ∆τ is the time interval.Equation ( 20) takes the form: with The current outlet temperature vectors can be calculated by the next expression: Energies 2017, 10, 1398 6 of 16 The accuracy of the solution of Equation ( 19) directly depends on the parameters in the Nusselt correlation.Precisions of the empirical parameters c Nu , c Re and c Pr will influence on the overall heat transfer coefficient and consequently lead to simulation errors of the outlet temperatures.Besides, fouling resistance also impacted the overall heat transfer coefficient of PHE.Therefore, an alternative approach to solve Equation (19), which does not depend on the parameters in the Nusselt correlation, will be applied in the following.

Model Linearization
In Equations ( 17) and ( 18), there are product terms of the temperature and velocity, which need to be linearized to implement the grey-box approach.For those simple variables, e.g., temperature and velocity, they can be linearized using the following equation: where θ o and δθ denote the value of variables on the balanced point and the increment of temperature and velocity variables, respectively.Complex variables, e.g., the heat transfer coefficient, are normally considered as the function of flow rates of hot side and cold side fluids.They ought to be linearized with the Taylor formula.In this paper, a one-order Taylor series is applied to obtain the linearization: The value of variables on the balanced point satisfy the following formulas: Equations ( 24) and ( 25) are substituted into Equations ( 17) and (18).Neglecting the high-order terms, and dropping the increment symbol "δ" for simplication, Equations ( 17) and ( 18) can be linearized as: .
where T, u, w are the following vectors: The system matrices of Equations (28) are as follows: Energies 2017, 10, 1398 7 of 16

Parameter Matrix Identification
In this section, a parameter identification method is developed to calculate the unknown parameter matrix.The method makes good use of the model structure derived in previous sections.The model has been developed is a second-order linear state space model and all the state variables are measureable.On account of the test data of system variables being merely available at discrete-time instants, the model is transformed into discrete-time state-space equations [17].Because .
The parameter matrix identification is to find the parameter matrices A id , B id , F id ∈ R 2×2 , which minimize the F-norm of the error matrix ε(M).
The solution of Equation ( 41) is as follows: where the symbol " †" represents the Moore-Penrose inverse of matrix, and for full rank the matrix can be calculated as follows:

Test Procedure
In order to validate the presented methods, a district heating system in Tianjin was tested.The system schematic is shown in Figure 2. The heat source and heating equipment in this district heating system are a ground-source heat pump (GSHP) and a fan coil (FC), respectively.In heating substations, the hot water from the heat source transfers heat to the cold water from the FC through the PHE (HISAKA PHE, RX-396B-TNUP-269, HISAKA Works (China) Co. Ltd., Shanghai, China).A picture of the PHE in the district heating system is shown in Figure 3.
For the inlet and outlet temperatures of PHE, the primary and secondary flow rates were measured every minute during the heating period of 2015-2016.
The inlet and outlet temperatures of PHE were measured by platinum resistance thermometers (PT100, ELECALL Electric Co. Ltd., Yueqing, China).The accuracy of the thermometers is ±(0.15 + 0.002 |t|) • C. The two flow rates were measured by means of Ultrasonic flow meters (Fluke Corporation, Washington, DC, USA).The accuracy of the flow meters is ±2%.The test data is recorded by a data logger.
According to [18], uncertainty analysis is conducted to validate the test results.The uncertainty of temperature is established to be ±2.06%.The uncertainty of the convective heat transfer coefficients is calculated to be ±3.12%.The overall heat transfer coefficient calculated by the temperature and flow rate is approximately ±4.57%, an acceptable uncertain value range in engineering applications [19].
The parameters of the PHE are presented in Table 1.The constant empirical parameters c Nu , c Re and c Pr are estimated by the manufacturer's data.

Test Procedure
In order to validate the presented methods, a district heating system in Tianjin was tested.The system schematic is shown in Figure 2. The heat source and heating equipment in this district heating system are a ground-source heat pump (GSHP) and a fan coil (FC), respectively.In heating substations, the hot water from the heat source transfers heat to the cold water from the FC through the PHE (HISAKA PHE, RX-396B-TNUP-269, HISAKA Works (China) Co. Ltd., Shanghai, China).A picture of the PHE in the district heating system is shown in Figure 3.
For the inlet and outlet temperatures of PHE, the primary and secondary flow rates were measured every minute during the heating period of 2015-2016.
The inlet and outlet temperatures of PHE were measured by platinum resistance thermometers (PT100, ELECALL Electric Co. Ltd., Yueqing, China).The accuracy of the thermometers is ±(0.15 + 0.002 |t|) °C.The two flow rates were measured by means of Ultrasonic flow meters (Fluke Corporation, Washington, DC, USA).The accuracy of the flow meters is ±2%.The test data is recorded by a data logger.
According to [18], uncertainty analysis is conducted to validate the test results.The uncertainty of temperature is established to be ±2.06%.The uncertainty of the convective heat transfer coefficients is calculated to be ±3.12%.The overall heat transfer coefficient calculated by the temperature and flow rate is approximately ±4.57%, an acceptable uncertain value range in engineering applications [19].
The parameters of the PHE are presented in Table 1.The constant empirical parameters Nu c , Re c and Pr c are estimated by the manufacturer's data.

Parameter Matrix Identification
During the cold period, test data over 6 days (10-15 January 2016) were chosen to identify the parameter matrices.The data were rearranged as Equations ( 35)-(38), and parameter matrix identification was performed by applying Equations ( 42) and (43).
After obtaining the parameter matrices, the two outlet temperatures T h,2 and T c,2 were calculated from Equation (34), where the two inlet temperatures T h,1 and T c,1 and the two flow rates u h and u c were substituted with the test data as model inputs.The calculated two outlet temperatures T h,2 and T c,2 were compared with the test data.Figures 4 and 5 show the comparison results in 20 h.The fit errors of outlet temperatures are calculated, and they are less than 0.2 • C. calculated from Equation (34), where the two inlet temperatures

Validation of the Proposed Methods
In order to validate the proposed methods, test data over 4 days (6-9 January 2016) were selected.The hot side flow rates ranged from about 0.15 to 0.4 m/s in this period.
For white-box method, the two outlet temperatures For grey-box method, according to Equation (34) the two outlet temperatures were calculated.The model inputs in 20 h, i.e., the two inlet temperatures

Validation of the Proposed Methods
In order to validate the proposed methods, test data over 4 days (6-9 January 2016) were selected.The hot side flow rates ranged from about 0.15 to 0.4 m/s in this period.
For white-box method, the two outlet temperatures T h,2 and T c,2 were calculated based on Equation (23).
For grey-box method, according to Equation (34) the two outlet temperatures were calculated.The model inputs in 20 h, i.e., the two inlet temperatures T h,1 and T c,1 and the two flow rates u h and u c , are depicted in Figures 6a,b and 7a,b.
For white-box method, the two outlet temperatures For grey-box method, according to Equation (34) the two outlet temperatures were calculated.The model inputs in 20 h, i.e., the two inlet temperatures  The plots of simulated vs. experimental outlet temperatures are presented in Figures 8a and 9a.It can be seen that the simulated results by the grey-box method are closer to the test data than those obtained by the white-box method.
The errors and of two outlet temperatures by the two methods are shown in Figures 8b  and 9b.It can be observed in the Figure 8b that the simulation error of hot outlet temperature T h,2 by grey-box method is restricted within 0.5 • C, and the error by white-box method can be to a maximum of 2.6 • C. Figure 9b shows that the simulation error of cold outlet temperature T c,2 by grey-box method is restricted within 1 • C, and the error by white-box method can be to a maximum of 2.2 • C.
The relative errors of the two methods with respect to the experimental data are shown in Figures 8c and 9c.It can be seen that the relative errors of T h,2 by grey-box method are less than 1%, and the error by white-box method is less than 5% in the majority from Figure 8c. Figure 9c shows that the relative error of T c,2 by the grey-box method is also less than 2%, and the error by white-box method is less than 4% in the majority.

Conclusions
In this paper, a second-order state space model of PHEs was developed.Then the model was solved by the white-box method and grey-box method, respectively.In the white-box method, the linearization of the state space model was carried out by implicit difference scheme and the overall heat transfer coefficient was calculated with empirical correlation of Nusselt number.In the grey-box method, one-order Taylor series were used to make the linearization of the nonlinear state space model and a newly developed parameter matrix identification method was established.
The simulation results of outlet temperatures by the grey-box and white-box method respectively were compared with the test data.The plots of predicted outlet temperatures by the grey-box method were closer to the test data than those obtained by the white-box method.The comparison of the FIT by the two methods indicates that the grey-box method performs better than the white-box method.
When the PHE is in operation for a period in actual project, fouling resistance also impacts the overall heat transfer coefficient of PHEs.The grey-box method does not introduce any approximation on the heat transfer model, thus it can refrain from errors caused by differences of parameters in the Nusselt correlation.In addition, Equation ( 7) consists of a steady-state correlation, but the white-box method uses the Nusselt in a dynamic simulation of the operative variables of the PHE.Therefore the grey-box model can be used in transient simulation in a more appropriate way and it is more suitable and accurate for the predication of outlet temperatures than the white-box method for that.The presented results can produce useful guidelines for optimal control of heating system.
In our future works, based on the proposed grey-box method, the controller design and its application for the regulation of the flow rates of district heating systems will be investigated.

Figure 1 .
Figure 1.Heat transfer process of the plate heat exchanger (PHE).

where wρ
is the fluid density (kg/m 3 ); p c is specific heat (J/(kg•K)); p δ is the distance between two neighboring plates (m); p b is the effective channel width which is obtained with the effective plate exchange area dividing by the length of the flow channel (m); k denotes the overall heat transfer coefficient between the hot and cold fluid (W/(m 2 •K)); c T and h T are respectively the channel temperatures of cold and hot fluid (K); and c u and h u are the flow velocity of cold and hot fluid (m/s), respectively.Equations (

Figure 1 .
Figure 1.Heat transfer process of the plate heat exchanger (PHE).

Figure 3 .
Figure 3. HISAKA PHE in the district heating system.

Figure 3 .
Figure 3. HISAKA PHE in the district heating system.
u were substituted with the test data as model inputs.The calculated two outlet temperatures the test data.Figures4 and 5show the comparison results in 20 h.The fit errors of outlet temperatures are calculated, and they are less than 0.2 °C.

Figure 4 .
Figure 4. Fit results of hot side outlet temperature by grey-box method.Figure 4. Fit results of hot side outlet temperature by grey-box method.

Figure 4 . 17 Figure 5 .
Figure 4. Fit results of hot side outlet temperature by grey-box method.Figure 4. Fit results of hot side outlet temperature by grey-box method.Energies 2017, 10, 1398 11 of 17

TFigure 5 .
Figure 5. Fit results of cold side outlet temperature by grey-box method.

Figure 8 .
Figure 8.(a) Comparison of simulated ,2 h T and grey-box method; (c) The relative error of

T
by white-box and grey-box method.

Figure 8 .
Figure 8.(a) Comparison of simulated T h,2 with test data; (b) Simulation error of T h,2 by white-box and grey-box method; (c) The relative error of T h,2 by white-box and grey-box method.

Figure 9 .
Figure 9. (a) Comparison of simulated,2 c T with measurement data; (b) Simulation error of ,2 c T

T
by white-box and grey-box method.

Figure 9 .
Figure 9. (a) Comparison of simulated T c,2 with measurement data; (b) Simulation error of T c,2 by white-box and grey-box method; (c) The relative error of T c,2 by white-box and grey-box method.

Table 1 .
Parameter of the plate heat exchanger (PHE).