A hybrid transient model for simulation of air-cooled refrigeration systems: Description and experimental validation

Accepted Manuscript A hybrid transient model for simulation of air-cooled refrigeration systems: Description and experimental validation Jing Wu, Emilie Gagnière, Françoise Couenne, Boussad Hamroun, Thierry Latour, Christian Jallut PII:

S0140-7007(14)00269-2

DOI:

10.1016/j.ijrefrig.2014.10.008

Reference:

JIJR 2896

To appear in:

International Journal of Refrigeration

Received Date: 11 July 2014 Revised Date:

13 September 2014

Accepted Date: 1 October 2014

Please cite this article as: Wu, J., Gagnière, E., Couenne, F., Hamroun, B., Latour, T., Jallut, C., A hybrid transient model for simulation of air-cooled refrigeration systems: Description and experimental validation, International Journal of Refrigeration (2014), doi: 10.1016/j.ijrefrig.2014.10.008. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights A hybrid dynamic model for refrigeration systems is proposed.



Flows through heat exchangers are approximated by a cascade of CSTRs.



All the refrigerant thermodynamic properties are analytically determined.



Phase transition within heat exchangers is performed by matrix operations.



Pressure losses in evaporator are taken into account by lumped models.

AC C

EP

TE D

M AN U

SC

RI PT



ACCEPTED MANUSCRIPT

A hybrid transient model for simulation of air-cooled refrigeration systems: description and experimental validation

RI PT

Jing Wu1,*, Emilie Gagnière1, Françoise Couenne1, Boussad Hamroun1, Thierry Latour2, Christian Jallut1 1 Université de Lyon, Université Lyon 1, Laboratoire d’Automatique et de Génie des Procédés, UMR 5007, CNRS/UCBL, 43 Bd du 11 Novembre 1918, 69622 Villeurbanne Cedex, France 2 CIAT, Avenue Jean Falconnier, 01350 Culoz, France * (corresponding author. Tel.: +33 472 431 856; E-mail: [email protected])

M AN U

SC

Abstract In this paper are described a hybrid dynamic model for transient simulation of refrigeration systems as well as dynamic experiments that have been performed on an air/water heap pump. The machine under consideration is made of an evaporator, a condenser, an expansion valve, a variable speed scroll compressor and a receiver. The refrigerant and second fluid flows in heat exchangers are approximated by a cascade of Continuous Stirred Tank Reactors (CSTRs). This model is quite flexible since a unique structure is used for the evaporator and the condenser models according to different boundary conditions. This is due to the use of a switching procedure between different configurations based on a phase stability test that is designed to ensure the continuity of the system simulation. An analytical thermodynamic model of the refrigerant based on an equation of state is used. Good agreement between simulation results and experimental data is achieved.

1. Introduction

TE D

Keywords: Refrigeration system; Heat exchanger; Hybrid dynamic model; Thermodynamics.

EP

In the field of buildings refrigeration, the coupling of a Vapour Compression System (VCS) with a cold thermal storage can lead to significant cost savings (Rismanchi et al., 2012). Daytime electricity consumption can be shifted to night time when electricity is cheaper. Cold energy is stored when the cooling demand of the building is low, whereas the storage works in parallel with the VSC when the cooling demand is high. Consequently, the VCS can be subject to a wide range of operating conditions and transient periods. In order to manage automatically and safely such a situation, it is recommended to use first principles dynamic models for the design of automatic controllers (Rasmussen and Alleyne, 2004; Schurt et al., 2009; Wallace et al., 2012; Catano et al., 2013a,b) and to optimize the dynamic operating conditions of such coupled systems (Wang et al., 2007b).

AC C

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

There are other circumstances where VCS are operated under transient conditions so that dynamic modelling can be necessary. Cycling conditions (Ndiaye and Bernier, 2012) as well as shut-down or start-up severe transients (Li and Alleyne, 2010; Uhlmann and Bertsch, 2012; Schalbart and Haberschill, 2013) are examples of such situations. The objective of this study is to propose a first principles dynamic model of a VCS provided by the CIAT company. In future works, this VCS will be coupled to a PCM (Phase Change Materials) cold thermal storage and then be used for experimental control studies of building refrigeration systems. This model will be used to design a multivariable controller allowing regulating the VCS cooling capacity and the evaporator superheating by using the compressor rotation speed and the valve position as control variables (Romero et al., 2011). 1

ACCEPTED MANUSCRIPT Such a control is all the more necessary that the VCS is coupled to a thermal storage in order to satisfy a varying thermal demand.

M AN U

SC

RI PT

A detailed review of recent literature on VCS dynamic modelling is presented by Rasmussen (2012a). The primary dynamics of these systems are driven by fluctuating states of refrigerant in evaporators and condensers as well as by the thermal inertia of the exchangers’ walls. VCS evaporators and condensers are generally represented by an equivalent tubular heat-exchanger, the equivalent parameters of which allowing to represent the real exchanger properties (Rasmussen and Alleyne, 2004; Eldredge et al., 2008). Two approaches are then generally used to derive evaporators and condensers dynamic models (Bendapudi et al., 2008). The first one consists in deriving 1D models from refrigerant mass, energy balances and sometimes momentum balances (Zhang et al., 2009). The resulting PDEs are discretized in order to be numerically solved (Bendapudi et al., 2008; Hermes and Melo, 2008; Ndiaye and Bernier, 2012). A classical approach used in chemical engineering to perform such a discretization consists in approximating a flow by a cascade of Continuous Stirred Tank Reactors (CSTRs). Such a method has already been applied to VCS heat exchangers modelling (Shalbart and Haberschill, 2013). The second approach, the movingboundary method, has been proposed to get a low order lumped parameter model easier to solve (Willatzen et al., 1998a; Wang et al., 2007a; McKinley and Alleyne, 2008). It consists in considering separately single and two-phase zones of the refrigerant and to calculate the time dependent position of the boundaries between these zones. In some circumstances, one boundary can be removed from the heat exchanger refrigerant volume and switching between different sub-models is necessary (Willatzen et al., 1998b; McKinley and Alleyne, 2008). This situation, that is very well documented in automatic control literature (Song et al., 2012), is said to be hybrid in the sense that the dynamic evolution of a system is described by different sets of differential equations according to discrete events.

EP

TE D

In the model that we propose here, the compressor and expansion valve of the VCS are represented statically as it is classically done (Rasmussen, 2012a). The main effort has been devoted to the derivation of a dynamic model for the heat exchangers. We propose a 1D model that is approximated by a cascade of CSTRs (Shalbart and Haberschill, 2013). The main difficulty for evaporators and condensers dynamic modelling is due to the fact that the equations to be solved, which are derived from mass and energy balances, are different according to the refrigerant physical state (vapour, liquid or two-phase mixture). Furthermore, this physical state varies with position and time. This is a hybrid system: a switching procedure between the different sets of equations requires the consideration of discrete events associated with transitions between single- and two-phase refrigerant configuration within each CSTR and at each time. We propose here such a switching procedure that allows deriving a very compact and flexible model of heat exchangers in the sense that no assumption has to be made about the physical state of the refrigerant according to time and position. The model allows to determine this state and to switch automatically between the different configurations. It is then applicable for both the evaporator and the condenser. The compactness and flexibility characters of the model that we propose for the heat exchangers is very important for control purpose of a VCS subject to a wide range of operating conditions and transients. Such a hybrid formulation of the heat exchangers model is achievable since we use an analytical approach for the thermodynamic properties calculations based on an equation of state (Lemmon, 2003). As far as we know from open literature, this approach is original within the framework of VCS and heat exchangers dynamic modelling.

AC C

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

2

ACCEPTED MANUSCRIPT As mentioned by Rasmussen and Shenoy (2012b), few publications provide detailed experimental validation due to the challenge of proper initializing and tuning of the model. Most of them give either only simulation results or qualitative comparisons. In this paper, the model is validated by comparison with experimental data provided by the CIAT Company.

2. Experimental set-up

RI PT

The system under consideration is operated by the CIAT Company: it is shown on figure 1. The heat pump is mainly made of a scroll compressor, an air-cooled finned-tube condenser, a plate evaporator, an Electronic Expansion Valve (EEV) and a refrigerant storage receiver. The working refrigerant is R-410A.

EP

TE D

M AN U

SC

In order to perform experimental tests, the refrigeration system is placed in a 300 m 3 climate chamber to impose external conditions to the condenser. The ambient temperature can be varied between -15 °C and 46 °C and the relative humidity up to 95 %. Energy is provided to the evaporator through an external heat exchanger. In order to perform dynamic experiments, the machine is operated manually (the regulators of the machine are disconnected) and open loop dynamic experiments can be performed by manipulating the compressor rotation speed comp, the EEV position , the fan rotation speed f and the external liquid flow rate qa.

AC C

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Fig. 1: flow sheet of the CIAT system

3. Main features of the model 3.1. Global structure

In order to get a more efficient simulator for the model, a proper interconnection between the sub-models of each element of the system has to be chosen (Rasmussen and Shenoy, 2012b). Since the refrigerant total mass remains constant during the simulation, the compressor and expansion valve refrigerant mass flow rates qcompt  and qEEV t  are used as boundary conditions for the refrigerant mass balances of the two exchangers (see figure 1). As a matter of fact, the overall refrigerant mass balances for the two exchangers are as follows:

 3



ACCEPTED MANUSCRIPT Lcond

qcomp t   qEEV t  

  z ,t S

d

cond

dz

dt

Leva

qEEV t   qcomp t  

cond

0

d

  z ,t S eva

(1) eva

dz 

0

dt

dmrec ( t ) dt

Lcond

  z ,t S cond

0

cond

dz 

Leva

  z ,t S eva

eva

RI PT

Since refrigerant mass accumulation is assumed to occur only within the two exchangers and the receiver, it is evident from equations (1) to conclude that the total refrigerant mass

dz  mrec ( t ) is constant if qcompt  and qEEV t  are

0

M AN U

SC

imposed to the two boundaries of each heat exchanger. These mass flow rates are computed from EEV and compressor static models that also calculate the refrigerant thermodynamic  are used as inlet these conditions conditions at outlet of the compressor and of the EEV: boundary conditions for the evaporator and condenser model. The refrigerant thermodynamic conditions at outlet of the evaporator and condenser are then computed by using the heat exchangers models to serve as inlet conditions for the compressor and for the EEV. Inlet secondary fluids temperature and flow rates are given as boundary conditions for the secondary fluids models of the evaporator and condenser.

TE D

To sum up, the inputs of the model are the inlet secondary fluids temperatures and mass flow rates, the compressor rotation speed and the EEV position. The total refrigerant mass is given through the initial conditions holding in the heat exchangers and in the receiver and remains constant. As outputs, time variations of the refrigerant thermodynamic state are calculated at the inlet and outlet of each sub-system as well as according to the position inside the heat exchangers. Time variations of walls and secondary fluids temperatures profiles in the heat exchangers are also computed. 3.2. Thermodynamic model

EP

3.2.1. Single phase domain

In the single-phase domain, we use the Lemmon equation of state that has been designed by considering the R-410A refrigerant as a pseudo-pure fluid (Lemmon, 2003). The independent variables of this thermodynamic model are the mass density m and the temperature T. These variables are also the independent variables of our dynamic model.

AC C

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58  59 60 61 62 63 64 65

The Lemmon model provides analytical expressions for pressure, specific internal energy and enthalpy: Pm m ,T  , um m ,T  and hm m ,T  . From these expressions, we have Pm obtained derivatives expressions that are needed for the simulations – i.e. – m ,T  ,

m

Pm  u u  ,T . m ,T  and m  m ,T , m  T T m m

 3.2.2. Two-phase domain   The Lemmon model provides expressions for the dew and bubble pressures. On the figure 2 are shown the evolutions of the two pressures that are almost identical. The pseudo-pure 4

ACCEPTED MANUSCRIPT fluid assumption for the R-410A is very good since a unique two-phase equilibrium pressure dPd can be considered Pd T  as well as its derivative T  that we also need for the dT simulations. In order to simplify the two-phase domain calculations, we have determined empirical polynomials from available data (ASHRAE, 2001) for uv T  , ul T  , v T  , l  T  , hv T  , hl T  and, if x is the vapour mass fraction, the two-phase system  properties are given as functions of the two independent variables T and x: 

RI PT





1 x 1 x     d T , x   v T   l T  ud T , x   xu v T   1  x ul T 

TE D

M AN U

SC

hd T , x   xhv T   1  x hl T 

(2)

EP

Fig. 2: bubble and dew point curves of R-410A calculated with the Lemmon model

3.3. Network representation of the equivalent tubular heat exchanger Only one model is derived that is applied to both condenser and evaporator. This model is based on the equivalent tubular heat exchanger 1D model that is widely used for VCS modelling (Rasmussen, 2012a). The parameters describing the equivalent tubular heat exchanger are determined from the real geometry of the heat exchangers such as the heat transfer surface area, the mass of the heat exchanger wall, the refrigerant volume, etc.

AC C

1 2 3 4 5 6 7 8 9 10  11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

The main assumptions of the model are: (a) the refrigerant pressure is assumed to be uniform so that momentum balance equation is not included; (b) in two-phase zone of the refrigerant, the liquid and vapour are in thermodynamic equilibrium; (c) the axial conduction effects are neglected; (d) the secondary fluid is assumed to remain as a single phase fluid. Its density and heat capacity are assumed to be constant. 5

ACCEPTED MANUSCRIPT The use of networks for dynamic modelling is very well documented in chemical engineering (Gilles, 1998; Mangold et al., 2002; Couenne et al., 2008). In our case, the network is made of lumped elements: such a network has been used for heat exchangers and VCS modelling (Schalbart and Haberschill, 2013) as well as chemical processes (see for example Choulak et al., 2004; Goma Bilongo et al., 2012).

M AN U

SC

RI PT

The network that we use for the modelling of the equivalent tubular heat exchanger is represented on the figure 3 in the case of a co-current configuration. The refrigerant and secondary fluid flows are assumed to be equivalent to a serial arrangement of N identical CSTRs. The radial transfer through the heat exchanger wall is represented by using the resistance-capacitance electrical analogy (see for example Bond et al., 2013 or Maestre et al., 2013 for recent applications of this analogy).

Fig. 3: equivalent tubular heat exchanger network

EP

4. Model equations

TE D

Once the network is defined, the balance equations are written for each CSTR and capacitance element of the network. The spatial precision of the network is related to the number of CSTRs and capacitances that are used.

4.1. Heat exchanger model

AC C

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

4.1.1. Balance equations

The mass and energy balance equations for the refrigerant in the kth CSTR are (see figure

3):

dk  qk1  qk dt  d (Vk  k uk )  qk 1hk 1   iw Aiw (Twk  Tk )  qk hk  dt  Vk  d k uk   k duk   qk 1hk 1   iw Aiw (Tw  Tk )  qk hk k dt    dt

Vk



where Vk is the volume of the kth CSTR.



6

(3)

(4)

ACCEPTED MANUSCRIPT The energy balance equations for the tube wall and the secondary fluid kth volumes are: m w c pw ma c pa

dTwk dt dTak dt

  iw Aiw ( Tk  Twk )   ow Aow( Tak  Twk )

(5)

 c pa qa ( Tak 1  Tak )   ow Aow( Twk  Tak )

(6)

RI PT

The heat transfer coefficients  iw and  ow are calculated by using available correlations in open literature (see Appendix 1). 4.1.2. Structure of the heat exchanger model



1

1

M AN U

SC

The solution of equations (3) and (4) depends on the fact that, within the kth CSTR, the refrigerant is a single- or a two-phase system. At each time, a test has to be performed to answer this question and then the simulation can continue according to the resulting configuration. We describe here a unique mathematical structure that can handle globally this problem. For such hybrid systems, the time evolution of the system state variables depends on continuous equations as well as discrete events. In order to derive a unique global numerical solution representing all the configurations, we define a unique state vector X T  1 T1 x1 Tw Ta N TN x N Tw Ta where  k , Tk and x k are respectively the refrigerant mass density, the temperature and the vapour mass fraction within the kth CSTR and N is the total number of CSTRs. The last state variables are the wall and second fluid temperatures, the evolution of which being represented continuously by the balance equations (5) and (6). N

N



TE D

Let us first consider separately the single- and the two-phase configurations for the refrigerant. 4.1.2.1. Single-phase configuration

EP

Let us consider that, at time t, the refrigerant content of the kth CSTR is a single phase represented by the Lemmon equation of state (Lemmon, 2003). Due to the Lemmon thermodynamic model structure, in single-phase configuration, mass density and temperature are selected to describe the refrigerant thermodynamic state, which are also employed to express the derivatives of the specific internal energy:

AC C

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29  30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

dum ,k  um ,k   dt   m ,k

 d m ,k  um ,k   dt   T  k T

 dTk    dt

(7)

Another advantage to consider temperature as one of the state variables is that it has to be known to calculate heat fluxes. By combining equations (4) and (7), one can easily derive the following equation involving the temperature and mass density time derivatives:  d m ,k   u   u  u k   m ,k  m ,k T    m ,k  m ,k Vk        dt   Tk  m ,k      q k 1 hm ,k 1   iw Aiw ( Twk  Tk )  q k hm ,k 7

 dTk     dt 

(8)

ACCEPTED MANUSCRIPT As far as we have defined a unique state vector including the vapour mass fraction, we have to include the following dummy equation (9): dxm ,k 0 dt

(9)

dX  Bm dt

(10)

SC

Am ( X ) 

RI PT

Finally, equations (3), (5), (6), (8) and (9) written for k = 1 to N can be represented according to the following implicit system of 5N ordinary differential equations (see Appendix 2 for the content of the matrix Am and the vector Bm):

4.1.2.2. Two-phase configuration

M AN U

In two-phase configuration, the thermodynamic properties of each phase only depend on one intensive variable since thermodynamic equilibrium is assumed. Therefore, the two-phase refrigerant properties can be described using the temperature and the vapour mass fraction as well as their time derivatives according to equations (2):

TE D



(11) (12)

By combining equations (4), (11) and (12), one can obtain the following equation involving the vapour mass fraction and temperature time derivatives:     d ,k V k   u k     x k 

  u    k  d ,k T  x k

  T

EP



dd,k d,k  dx k d,k  dTk       dt  x k T dt  Tk x dt dud,k ud,k  dx k ud,k  dTk       dt  x k T dt  Tk x dt

 dx k    d ,k   u   dt  k  T k   

 q k 1 hd ,k 1   iw Aiw ( Twk  Tk )  q k hd ,k

  u    k  d ,k x  Tk

  dTk     x  dt

   

(13)

AC C

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Finally, by combining equations (3), (5), (6) (11) and (13), the following implicit system of 5N ordinary differential equations is obtained (see Appendix 2 for the content of the matrix Ad and the vector Bd): Ad ( X ) 

dX  Bd dt

(14)

Let us notice that equation systems (9) and (14) depend on the refrigerant mass flow rates that have to be calculated. To this end, we consider the time evolution of the pressure within the kth CSTR according to the single or two-phase configurations. 4.1.2.3. Mass flow rates calculations 8

ACCEPTED MANUSCRIPT If the kth CSTR configuration is single-phase, the pressure is function of the temperature and the mass density so that:

dPm,k Pm,k  dk Pm,k  dTk       dt  k T dt  Tk  dt

By combining equations (3), (8) and (15), one can express the pressure time derivative within the kth CSTR as function the two adjacent flow rates:

RI PT



(15)

dPm ,k  c1k qk 1  c2k qk  c3k dt

(16)

SC

If the kth CSTR configuration is two-phase, the pressure is the liquid-vapour equilibrium pressure that is only function of the temperature:



(17)

M AN U

dPd,k dPd dTk    dt dTk  dt

By combining equations (3), (11), (13) and (17), the pressure time derivative is expressed as a function of the two adjacent flow rates: dPd ,k  d1k qk 1  d 2k qk  d 3k dt

(18)

TE D

c1k , c2k , c3k , d 1k , d 2k and d 3k are functions of the refrigerant properties within the kth CSTR

as well as the temperature difference (Tw  Tk ) (see Appendix 2). k

The equations (16) and (18) can be rewritten according to a unique form:

 dPk  v kc1  (1 v k )d1 qk1  v kc 2  (1 v k )d2 qk  v kc 3  (1 v k )d3 dt k



EP



k



k

k



k

k

(19)

where vk is 1 or 0 according to the fact that the kth CSTR configuration is single- or two-phase. 

AC C

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Since we have assumed uniform pressure in the heat exchangers, the following constraint can be assumed whatever the configuration:

dPk dPk1  dt dt



(20)

The combination of equation (16) or (18) and (19) gives N-1 relations to calculate the N-1 unknown mass flow rates qk (k=1 ··· N-1). 4.1.2.4. Global and unique structure of the heat exchanger model

9

ACCEPTED MANUSCRIPT Let us first consider the dynamical part of the model. The equations systems (10) and (14) can be expressed according to a unique global generic model:

A(X) 

dX B dt

(21)

RI PT

with A  E  Am  ( I  E )  Ad and B  E  Bm  (I  E)  Bd where E is a (5N)×(5N) diagonal matrix with vk as the kth  diagonal element and I is a (5N)×(5N) identity matrix. During the computation, if at least one CSTR configuration is modified, A and B are updated.

 As far as the mass flow rates calculations are concerned, a unique structure can also be proposed. Under the assumption that the pressure is uniform throughout the heat exchanger, the constraint equation (20) can be expressed as follows: k

k

k

k1

k1

k1

k1

k1

k

(22)

k

When applied to the CSTRs 1 to N-1, the resulting N-1 equations (22) can be expressed as follows:

M AN U



k1

k

SC

vkc1  (1 vk )d1 qk1  vkc2  (1 vk )d2  vk1c1  (1 vk1)d1 qk v k1c 2  (1 v k1 )d2 qk1  v k1c 3  (1 v k1 )d3  v k c 3  (1 v k )d3

CQ D

(23)

The matrix Am, Ad, C, the vectors Bm, Bd, D and the expressions of c1k , c2k , c3k , d 1k , d 2k

d 3k are given in the Appendix 2.

4.1.3. Detection of the configuration change

EP



TE D

qN 1. Let us recall that q0 and qN are known. They are given by the where QT  q1 q2  compressor and EEV static models as boundary conditions for the heat exchangers models within the framework of the heat pump model (see section 3.1).

The configuration change detection is based on thermodynamic stability tests. Such stability tests are commonly used for mixtures at given pressure and temperature (Michelsen, 1982). In our case, a phase stability test at given temperature and volume is more suitable (Souza et al., 2006) since it is based on the value of the refrigerant density as it is calculated from the balance equations. For a single phase at a given temperature and a given density, if the tangent plane to the Helmoltz free energy is somewhere greater than the Helmoltz free energy, the single phase under consideration is unstable. This formulation is equivalent to the fact that a single-phase system having a mass density located between the saturated liquid and vapour densities is less stable than a system made of the two phases. According to the initial situation, the vapour mass fraction or the mass density is used to apply this stability test to each CSTR and at each time.

AC C

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

From an initial two-phase configuration, if the vapour mass fraction reaches a value close to 0 or 1, the configuration becomes single phase. From an initial single-phase configuration, if the mass density becomes greater than the saturated vapour mass density or lesser than the saturated liquid mass density at the current pressure, the configuration becomes two-phase.

10

ACCEPTED MANUSCRIPT 4.1.4. Heat exchanger model simulation examples

RI PT

We provide here some simulation examples of the heat exchanger model in order to illustrate its flexibility. A fourth-fifth order Runge-Kutta method with adaptive step size (ode45) is used within the framework of the Matlab software. The detection of configuration change is achieved by using the option “event” of the function ode45. In fact, the integration process is stopped when a configuration change event occurs: the state variables as well as the model matrices are updated before the integration to continue.

EP

TE D

M AN U

SC

On figure 4, the dynamical responses of the heat exchanger to step changes of the secondary fluid inlet temperature are shown. We can observe that the heat exchanger runs initially as an evaporator and finally becomes a condenser when the secondary fluid inlet temperature is less than the refrigerant temperature over 300 s.

AC C

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Fig. 4: dynamical responses to step changes of the secondary fluid inlet temperature

The responses to step changes of the outlet mass flow rate are represented on figure 5. The heat exchanger runs initially as an evaporator and then becomes a condenser when the outlet mass flow rate is less than the inlet mass flow rate for 20 s. This is due to the fact that the refrigerant temperature has increased with the pressure increase. In order for the system to go back to its initial state, the accumulated mass of refrigerant has to be removed by increasing the outlet mass flow rate for 20 s.

11

SC M AN U

TE D

Fig. 5: dynamical responses to step changes of the outlet mass flow rate

EP

These simulation examples show the high flexibility of the heat exchanger model that we have proposed. It can be used without any prior assumption concerning the heat transfer direction. By applying the same switching procedure to the secondary fluid flow direction, one can also plan to simulate the switching between co- and counter-current configurations. Such a switching procedure can for example be used to simulate a change of the heat pump configuration between refrigeration and heating. 4.2. Scroll compressor static model

AC C

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

RI PT

ACCEPTED MANUSCRIPT

From data provided by DANFOSS Company, we have determined empirical correlations for the volumetric v and isentropic is efficiencies of the compressor. For confidentiality reasons, we cannot provide these correlations but simply explain the way there are included in the model.

 scroll compressor  model allows calculating the compressor mass flow rate The static according to the following equation: qcomp  aVacompv





(24)

where Va is the vapour pocket volume at the compressor inlet and  a the compressor inlet refrigerant mass density.

12



ACCEPTED MANUSCRIPT The compressor outlet enthalpy is calculated by using the isentropic efficiency is :

is 

hbis  ha hb  ha



(25)

The outlet enthalpy for the isentropic case is obtained by solving the following equations:

 

  

 sbis bis ,Tbis  sa a ,Ta  is is is  Pb b ,Tb  Pb



RI PT



so that



hbis  hbis bis ,Tbis



(27)

SC



(26)



M AN U

The Lemmon equation (Lemmon, 2003) is used for these calculations. 4.3. Receiver

EP

TE D

If the pressure in the receiver can be assumed to be close to that of the evaporator, the receiver can be integrated in the evaporator model (Eldredge et al., 2008). Within the framework of our heat exchanger model, it is simply a supplementary CSTR that is assumed to run adiabatically (see figure 6).

Fig. 6: evaporator and receiver network

The equations (21) and (23) are extended in order to include the receiver mass and energy balances. The receiver inlet mass flow rate becomes the new first boundary condition for the mass balances, the evaporator outlet mass flow rate remaining the second one.

AC C

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

4.4. Electronic expansion valve The refrigerant mass flow rate through the expansion valve is calculated by:

qEEV  Cd AEEV 2c Pc  Pd 



(28)

where Cd is the mass flow coefficient. Cd is represented by using an empirical correlation proposed by Chen et al. (2009): Cd  a  b qEEV

(29) 13



ACCEPTED MANUSCRIPT The expansion is assumed to be an isenthalpic process, so that:

hc  hd

(30)

The parameters a and b are estimated from the experimental results.



RI PT

5. Experimental results and model validation 5.1. Experimental results analysis

TE D

M AN U

SC

In order to validate our model, we have used experimental data obtained by the CIAT Company by manipulating the scroll compressor rotation speed according to the time variations represented on figure 7. The experimental data and simulation results are represented in dimensionless forms for confidentiality reasons.

Fig. 7: compressor rotation speed input

EP

One can see on figure 8 that the inlet and outlet condenser pressures are very close so that the pressure uniformity assumption is valid for the condenser. On the figure 9 are represented the receiver inlet and outlet pressures as well the evaporator outlet pressure. It can be seen that the pressure loss due to the evaporator is significant while the one due to the receiver is comparatively low. Including the receiver into the evaporator model as described in section 4.3 is then justified.

AC C

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

14

Fig. 8: dimensionless condenser pressure Fig. 9: dimensionless evaporator and receiver pressure

M AN U

SC

In order to keep the benefit of the heat exchanger model flexibility as it has been illustrated in section 4.1.4, we propose to represent the pressure losses due to the evaporator by using two lumped pressures losses models according to the figure 10, the pressure within the receiver and the evaporator remaining uniform. These pressure losses are then included in the static part of the model. Empirical correlations for these pressure losses are obtained from experimental results.

TE D

Fig. 10: evaporator pressure losses representation

5.2. Comparison between experiments and simulations

EP

The comparisons between experimental data and model outputs are presented on figures 11, 12, 13 and 14. It can be seen from these figurs that good agreement between simulation results and experimental data is achieved. Globally, the simulation gives reasonable results with errors being less than 10% at most of time.

AC C

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

RI PT

ACCEPTED MANUSCRIPT

Fig. 11: dimensionless condenser pressure

Fig. 12: dimensionless condenser temperature

15

Fig. 13: dimensionless evaporator pressure

Fig. 14: dimensionless evaporator temperature

TE D

M AN U

SC

We can see from the figure 15 that the total mass of system remains almost the same throughout the whole simulation period. Actually, the tiny inaccuracies are due to the state variables reinitialization in the configuration change operations. The calculated mass flow rates through the compressor and the EEV are depicted on the figure 16. The two flow rates are identical in stationary conditions, whereas they separate for transient changes.

6. Conclusion

EP

Fig. 15: dimensionless simulation mass distribution Fig. 16: dimensionless simulation mass flow rates

AC C

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

RI PT

ACCEPTED MANUSCRIPT

In this paper, a model for simulating the dynamic behaviour of refrigeration systems is described. The assumptions that have been made for compressor and EEV are classical: they have been considered statically. The main originality of the proposed model is related to the hybrid 1D heat exchanger model. This model has been discretized by using a network of CSTRs for the fluid phases and of resistances and capacitances for the heat exchanger wall. As far as the refrigerant is concerned, a hybrid model is proposed. Within each CSTR, the refrigerant fluid state is either single- or two-phase. According to the configuration, the equations derived from mass and energy balances are different and a switching procedure is organised by matrix operations since all the thermodynamic properties are analytically calculated. This procedure allows solving the balance equations in all the configurations. Separated simulations of this heat exchanger model have been provided in order to illustrate its high flexibility. It has been shown for example that the model is able to simulate

16

ACCEPTED MANUSCRIPT automatically a change in the way the heat exchanger runs, from an evaporator use to a condenser use. As far as the heat pump dynamic model is concerned, comparisons of simulated data with experimental results are given. Provided that lumped pressure losses are included in the evaporator model, the model simulations and experimental results are close one to each other.

RI PT

7. Acknowledgements

M AN U

SC

This work has been funded by the French National Research Agency (ANR) within the framework of the project ANR-11-SEED-0004-02 ACLIRSYS (Advanced Control for Low Inertia Refrigeration Systems). This project includes 4 industrial companies (CIAT, Cristopia, DANFOSS, and CMDL MANASLU Ing.) and 3 public research laboratories (LAGEPUniversité Lyon 1, LCIS-Grenoble-INP, LGP2S-CNAM Paris). The authors gratefully acknowledge the technical support provided by DANFOSS Company who provided the compressor static model.

Appendix 1: heat transfer correlations

The correlation proposed by Wang et al. (2000) is employed for calculating air-side heat transfer coefficient for the outside of the finned tubes: cPair Prair

TE D

 air   f  airVair

2/3

(0.086 Re Dc

P3

N P4 (

Fp Dc

) P5 (

Fp Dh

) P6 (

Fp Pt

) 0.93 )

(A1-1)

The fin efficiency ηf is determined by the correlation developed by Chen and Hsu (2007):

  f  (0.208  1)(0.40783  4.17947105  Re d

EP

S

9

2.598 10 Re d 2 )(1



(A1-2)

25.61 ) (S /do ) 0.55 Re d

AC C

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

The Shah correlation for condensation (1979) is used to calculate the refrigerant heat transfer coefficients in the two-phase zones of the condenser:



 cd   l (1  x ) 0.8  

3.8 x 0.76 (1  x ) 0.04   Pr 0.38 

(A1-3)

where αl is the heat transfer coefficient assuming all the mass flowing as liquid and is calculated by the Dittus-Boelter equation as:

 l  0.023 Rel 0.8 Prl 0.4 l / D

(A1-4)

17

ACCEPTED MANUSCRIPT The equation (A1-4) is also used to calculate the refrigerant side heat transfer coefficients in single-phase zones and water-side heat transfer coefficient in evaporator. The refrigerant heat transfer coefficients of single-phase zones in condenser are computed by the manufacturer. In two-phase zones of evaporator, the refrigerant heat transfer coefficient is calculated by the correlation proposed by Hsieh and Lin (2002):

 ed   l (88Bo 0.5 )

RI PT

(A1-5)

Appendix 2: heat exchanger model matrices

M AN U

SC

0  0   Am1  Ad1  Bm1   B d1      B  B  Am2 Ad 2 Am    , Ad    , Bm   m2  , Bd   d 2              B AmN  Ad N  m  0  0  Bd N   N 0 0 0 1 0 1 0 0 0 0  0 W W 0 0  Tk xk uk  a k  k bk  k 0 0 0  with Amk   0 0 1 0 0 , Ad k  1  RTk  R xk 0 0 ,    0 0 0 1 0 0 0 0 1 0  0 0 0 0 1  0 0 1 0 0

( qk 1  qk ) / V     hm ,k 1qk 1  hm ,k qk   iw Aiw ( Twk  Tk ) / V   Bmk  0   A ( T  T )   A ( T  T ) /( m c ) wk ow ow ak wk w Pw  iw iw k   c Pa qa ( Tak 1  Tak )   ow Aow ( Twk  Tak ) /( ma c Pa ) ( qk 1  qk ) / V     hd ,k 1qk 1  hd ,k qk   iw Aiw ( Twk  Tk ) / V   Bd k  0   A ( T  T )   A ( T  T ) /( m c ) wk ow ow ak wk w Pw  iw iw k  c q ( T  T )   A ( T  T ) /( m  ak ow ow wk ak a c Pa )  Pa a ak 1





 





 

EP

 

TE D

 

 0   R1  L1q0   M 1  O1 0    L2 M 2  O2    R2     C 0   , D  0       0 LN  2 M N  2  O N  2 RN 2    0 0  0 LN 1 M N 1    RN 1  O N 1q N  with

AC C

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Lk  vk c1k  ( 1 vk )d1k , M k  vk c2k  ( 1  vk )d 2k  vk 1c1k 1  ( 1  vk 1 )d1k 1 ,

Ok  vk 1c2k 1  (1  vk 1 )d 2k 1 , Rk  vk 1c3k 1  ( 1  vk 1 )d 3k 1  vk c3k  ( 1  vk )d 3k ,

d k ( uk  ak  k ) hk 1d k d ( u  ak  k ) hk d k   , c2 k   ck  k k , bk  k bk  k bk  k bk  k d k  iw Aiw ( Twk  Tk ) dP hk 1  Wxk / Rxk dP Wxk / Rxk  hk c3k  , d1k  , d 2k  , bk  k dTk WTk  Wxk RTk / Rxk dTk WTk  Wxk RTk / Rxk c1k  ck 

d 3k 

dP  iw Aiw ( Twk  Tk ) , dTk WTk  Wxk RTk / Rxk

18

ACCEPTED MANUSCRIPT uk u P P )T , bk  ( k ) , ck  ( )T , d k  ( ) , Tk  k  k Tk ( u )k   ( u )k Rxk  ( k )T , RTk  ( k )x , Wxk  ( )T , WTk  ( )x xk Tk xk Tk

where ak  (

Nomenclature heat exchange surface area (m2)

Bo:

boiling number, Bo 

CSTR:

Continuous Stirred Tank Reactor

Cd : cp:

mass flow coefficient  specific heat capacity (J K-1 kg-1)

Dc:

Fin collar outside diameter (Do+2δ) (m)

Do:

tube outside diameter (m)

Dh:

hydraulic diameter (m)

Fp:

fin pitch (m)

G:

mass flux (kg m-2 s-1)

h:

specific enthalpy (J kg-1)

L:

channel length (m)

m:

mass (kg)

q:

mass flow rate (kg s-1)

N:

number of tube row, number of CSTRs

P:

pressure (Pa)

Pr:

Prandtl number

Pt:

transverse tube pitch (m)

Re:

Reynolds number

s: S: T: t:

RI PT

A:



EP

TE D

M AN U

SC

hlv

AC C

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

fin spacing (m)

transversal section (m2) temperature (K) time (s)

u:

specific internal energy (J kg-1)

vm:

specific volume of the vapour-liquid mixture (m3 kg-1)

V:

volume of CSTR (m3)

Vair:

air velocity (m s-1)

x:

vapour mass fraction

19

ACCEPTED MANUSCRIPT Greek symbols heat flux (W m-2)

α:

heat transfer coefficient (W m-2 K-1)

δ:

fin thickness (m)

∆hlv:

vaporization enthalpy (J kg-1)

:

rotation speed (tr min-1)

:

EEV position

RI PT

:

Subscripts a, b, c, d: characteristic points secondary fluid

comp:

compressor

cond:

condenser

d:

two-phase system

eva:

evaporator

EEV:

electronic expansion valve

f:

fan

iw:

inside of tube

l:

liquid

m:

one-phase system

ow:

outside of tube

rec:

receiver

v:

vapour

w:

tube wall

EP

TE D

M AN U

SC

a:

AC C

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

References

ASHRAE, 2001. Fundamentals Handbook. Bendapudi, S., Braun, J.E., Groll, E.A., 2008. A comparison of moving-boundary and finite-volume formulations for transients in centrifugal chillers. Int. J. Refrigeration 31(8), 1437-1452. Bond D.M.E., Clark W.M., Mark Kimber M., 2013. Configuring wall layers for improved insulation performance. Appl. Energ. 112, 235-245. Catano J., Zhang T., Wen J.T., Jensen M.K., Peles Y., 2013a. Vapor compression refrigeration cycle for electronics cooling – Part I: Dynamic modeling and experimental validation. Int. J. Heat Mass Transfer 66, 911–921. Catano J., Lizarralde F., Zhang T., Wen J.T., Jensen M.K., Peles Y., 2013b. Vapor compression refrigeration cycle for electronics cooling – Part II: gain-scheduling control for critical heat flux avoidance. Int. J. Heat Mass Transfer 66, 922–929.

20

ACCEPTED MANUSCRIPT

EP

TE D

M AN U

SC

RI PT

Chen, H.T. and Hsu, W.L., 2007. Estimation of heat-transfer characteristics on a vertical annular circular fin of finned-tube heat exchangers in forced convection. Int. J. Heat Mass Transfer 51(7-8), 1920-1932. Chen, L., Liu, J.H., Chen, J.P., Chen Z.J., 2009. A new model of mass flow characteristics in electronic expansion valves considering metastability. Int. J. Therm. Sci. 48(6), 1235-1242. Choulak S., Couenne F., Le Gorrec Y., Jallut C., Cassagnau P., Michel A., 2004, A generic dynamic model for simulation and control of reactive extrusion, Ind. Eng. Chem. Res. 43(23), 7373-7382. Couenne F., Jallut C., MAschke B., Tayakout M., Breddveld P., 2008. Structured modeling for processes: A thermodynamical network theory. Comput. Chem. Eng. 32(6), 1120-1134. Eldredge, B.D., Rasmussen, B.P., Alleyne, A.G., 2008. Moving-boundary heat exchanger models with variable outlet phase. J. Dyn. Syst.-T. ASME 130(6), 1003-1012. Gilles E.D., 1998. Network theory for chemical processes. Chem. Eng. Technol. 21(2), 121-132. Goma Bilongo, T., Couenne, F., Jallut, C., Le Gorrec, Y., Di Martino, A., 2012. Dynamic modeling of the reactive twin-screw co-rotating extrusion process: experimental validation by using inlet glass fibers injection response and application to polymers degassing. Ind. Eng. Chem. Res. 51(35), 11381-11388. Hermes, C. J. L., Melo, C., 2008. A first-principles simulation model for the start-up and cycling transients of household refrigerators. Int. J. Refrigeration 31(8), 1341-1357. Hsieh, Y.Y. and Lin, T.F., 2002. Saturated flow boiling heat transfer and pressure drop of refrigerant R-410A in a vertical plate heat exchanger. Int. J. Heat Mass Transfer 45(5), 1033-1044. Lemmon, E. W., 2003. Pseudo-pure fluid equations of state for the refrigerant blends R-410A, R- 404A, R-507A and R-407A. Int. J. Thermophys. 24(4), 991-1006. Li B., Alleyne A.G., 2010, A dynamic model of a vapour compression cycle with shut-down and start-up operations. Int. J. Refrigeration 33, 538-552. Maestre R.I., Gonzalez Gallero F.J., Alvarez Gomez P., Mena Baladés J.D., 2013. Performance assessment of a simplified hybrid model for a verticalground heat exchanger. Energ. Buildings 66, 437 – 444. Mangold M., Motz S., Gilles E.D., 2002. A network theory for the structured modelling of chemical processes. Chem. Eng. Sci. 57, 4099-4116. McKinley, T. L. and Alleyne, A.G., 2008. An advanced nonlinear switched heat exchanger model for vapor compression cycles using the moving-boundary method. Int. J. Refrigeration 31(7), 1253-1264. Michelsen, M. L., 1982. The isothermal flash problem. Part I. Stability. Fluid Phase Equilibria, 9, 1-19. Ndiaye D., Bernier M., 2012, Transient model of a geothermal heat pump in cycling conditions -Part A: The model. Int. J. Refrigeration 35, 2110-2123. Rasmussen B.P., Alleyne A.G., 2004, Control-Oriented Modeling of Transcritical Vapor Compression Systems. J. Dyn. Syst.-T. ASME 126, 54-64 Rasmussen B.P., 2012a. Dynamic modeling for vapor compression systems-Part 1: Literature review. HVAC&R Res. 18(5), 934-955. Rasmussen B.P. and Shenoy, B., 2012b. Dynamic modeling for vapor compression systems-Part II: Simulation turorial. HVAC&R Res. 18(5), 956-973 Rismanchi B., Saidur R., Masjuki H.H., T.M.I. Mahlia T.M.I., 2012. Energetic, economic and environmental benefits of utilizing the ice thermal storage systems for office building applications. Energ. Buildings 50, 347354. Romero J.A., Navarro-Esbri J., Belman-Flores J.M., 2011, A simplified black-box model oriented to chilled water temperature control in a variable speed vapour compression system. Appl. Therm. Eng. 31, 329-335. Shah, M.M., 1979. A general correlation for heat transfer during film condensation inside pipes. Int. J. Heat Mass Transfer 22(4), 547-556. Shalbart P., Haberschill P., 2013. Simulation of the behaviour of a centrifugal chiller during quick start-up. Int. J. Refrigeration 36, 222-236. Schurt L. C., Hermes C. J. L., Neto A. T., 2009, A model-driven multivariable controller for vapor compression refrigeration systems. Int. J. Refrigeration 32, 1672-1682. Song, C., Wu, B., Li, P., 2012. A hybrid model-based optimal control method for nonlinear systems using simultaneous dynamic optimization strategies. J. Process Contr. 22(5), 852-860. Souza, A. T., Cardozo-Filho, L., Wolff, F., Guirardello, R., 2006. Application of interval analysis for Gibbs and Helmoltz free energy global minimization in phase stability analysis. Braz. J. Chem. Eng. 23(1), 117-124. Uhlmann M., Bertsch S. S., 2012, Theoretical and experimental investigation of startup and shutdown behavior of residential heat pumps. Int. J. Refrigeration 35, 2138-2149. Wallace M., McBride R., Aumi S., Mhaskar P., House J., Salsbury T., 2012, Energy efficient model predictive building temperature control. Chem. Eng. Sci. 69, 45–58. Wang, C.C., Chi, K.Y., Chang, C.J., 2000. Heat transfer and friction characteristics of plain fin-and-tube heat exchangers, par II: Correlation. Int. J. Heat Mass Transfer 43(15), 2693-2700.

AC C

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

21

ACCEPTED MANUSCRIPT

RI PT

SC M AN U TE D EP AC C

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Wang, F. Q., Maidment, G. G., Missenden, J. F., Tozer, R. M., 2007a. A novel special distributed method for dynamic refrigeration system simulation. Int. J. Refrigeration 30(5), 887-903. Wang, F. Q., Maidment, G. G., Missenden, J. F., Tozer, R. M., 2007b. The novel use of phase change materials in refrigeration plant. Part 2: Dynamic simulation model for the combined system. Appl. Therm. Eng. 27, 2902-2910. Willatzen, M., Pettit, N.B.O.L., Ploug-Sørensen, L., 1998a. A general dynamic simulation model for evaporators and condensers in refrigeration. Part I: moving-boundary formulation of two-phase flows with heat exchange. Int. J. Refrigeration 21(5), 398–403. Willatzen, M., Pettit, N.B.O.L., Ploug-Sørensen, L., 1998b. A general dynamic simulation model for evaporators and condensers in refrigeration. Part II: simulation and control of an evaporator, Int. J. Refrigeration 21(5), 404–414. Zhang W.-J., Zhang C.-L., Ding G.-L., 2009, On three forms of momentum equation in transient modeling of residential refrigeration systems. Int. J. Refrigeration 32, 938-944.

22