- Email: [email protected]
Dynamic simulation of electronic expansion valve controlled refrigeration system under different heat transfer conditions
Accepted Manuscript Title: Dynamic simulation of electronic expansion valve controlled refrigeration system under different heat transfer conditions Author: Yujia Shang, Aiguo Wu, Xing Fang, Yuwen You PII: DOI: Reference:
S0140-7007(16)30231-6 http://dx.doi.org/doi: 10.1016/j.ijrefrig.2016.07.020 JIJR 3393
To appear in:
International Journal of Refrigeration
Received date: Revised date: Accepted date:
18-2-2016 31-5-2016 25-7-2016
Please cite this article as: Yujia Shang, Aiguo Wu, Xing Fang, Yuwen You, Dynamic simulation of electronic expansion valve controlled refrigeration system under different heat transfer conditions, International Journal of Refrigeration (2016), http://dx.doi.org/doi: 10.1016/j.ijrefrig.2016.07.020. 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.
Dynamic simulation of electronic expansion valve controlled refrigeration system under different heat transfer conditions Yujia Shanga,
Aiguo Wua,
Xing Fanga,
Yuwen Youb
College of Electric Engineering and Automation, Tianjin University, No. 92, Weijin Road, Nankai District, Tianjin 300072, China b College of Automation, Tianjin urban construction University, Tianjin 300384, China Correspondence should be addressed to Xing Fang: [email protected] a
Highlights
Based on bubble theory, evaporator is divided into two heat transfer conditions. A switch criteria between two heat transfer conditions is presented. We develop a dynamic model that can predict the instability of system control loop. The oscillation of the refrigeration system is studied using the proposed model.
Abstract: This work presents a new dynamic lumped parameter model that is able to simulate the oscillation of electronic expansion valve (EEV) controlled refrigeration system due to unsuitable control or abrupt change of external parameters. Based on bubble dynamics theory, the heat transfer characteristic is used to divide the evaporator into the stable and unstable conditions, and a matching switch criteria is presented. Together with the model of main components, the dynamic behavior of oscillation of system parameters is simulated successfully. Model validation against experimental data demonstrates the capabilities of the modeling approach in predicting the instability of system control loop. This study provides the basis for further developing advanced control algorithm which will ensure that an EEV-controlled system refrigeration system has better stability. Keywords: Refrigeration system, Simulation, Superheat, Control, Stability
1.Introduction In vapor compression refrigeration systems, the instability of system control loop is the phenomena of oscillation of all parameters such as the refrigerant flow rate, refrigerant pressures and temperatures at different points. The instability will lead lower safety, lower life-span and higher energy consumption (Liang et al., 2010; Mithraratne and Wijeysundera, 2001). Dynamic simulation can be used to study the instability of control loop and provide guidance for design and control strategy of refrigeration system. The simulation of the instability on refrigeration system control characteristics have been noted in some studies for
thermostatic expansion valve (TEV) controlled evaporator refrigeration systems (Yasuda and Toufer, 1983; Ibrahim, 2001; Tahat et al., 2001; Lenger et al., 1998). With the progressively application of EEV, there have been a number of reported studies on the instability in EEV-controlled refrigeration systems (Li et al., 2004; Aprea and Mastrubllo, 2002; Chen et al., 2008; Fallahsohi et al., 2010; Leonardo et al., 2010). Among these researches, a model-driven multivariable controller for vapor compression-EEV refrigeration systems was presented (Leonardo et al., 2010). Fallahsohi et al. (2010) developed a predictive functional control method to control the degree of evaporator superheat (DS) with an
Page 1 of 17
electronic expansion valve. However, the system oscillation can be still observed in these advanced control refrigeration systems when the DS setting is low enough. Therefore, there is a need for a model, with low computational cost, to simulate the dynamic behavior of oscillation for an EEV-controlled refrigeration system and this model should be easily used to develop advanced control algorithms. The key point of modeling is to define the unstable condition, using a macroscopical Nomenclature
parameter of the expansion valve or evaporator, which is the cause of instability. With the frequency-domain analysis method, several characteristic parameters of TEV, such as the thermal capacitance of the bulb and thermal conductance between the bulb and wall were found to be associated with instability (Broersen and Van der Jagt, 1980; Yasuda and Toufer, 1983; Ibrahim, 2001). For the evaporator, there is an evaporator critical time constant, below which the system operates in a stable manner. And
A
state-space matrices
A
exchange surface area[m2]
rcrit
coefficient in Traviss, Baron and Rohsenow equation critical bubble radius
Av
valve opening [%]
void fraction
Cp
specific heat [J(kg K)-1]
com
global efficiency of the compressor
Cv
coefficient of valve discharge
v
compressor volumetric efficiency
C v, p
specific heat at constant pressure [J(kg K)-1]
surface tension [N m-1]
density [kg m-3]
Subscripts pts
two-phase region of evaporator and
K)-1]
C v ,v
specific heat at constant volume [J(kg
d
Diameter [m]
DS
degree of refrigerant superheat Reviewer #3: Response 2
f
compressor frequency [Hz]
H
Enthalpy [kJ kg-1)
H os
outlet enthalpy in isentropic compression
1
superheat region of condenser 2
process [J kg-1]
superheat region of evaporator and two-phase region of condenser
L
length [m]
3
subcooled region of condenser
m
mass flow rate [kg s-1]
avg
average
P
pressure [k Pa)
c
condenser
Pr
Prandtl number
cu
copper
Q
heat transfer rate [W)
chw
chilled water
q
output heat flux [W m-2]
clw
cooled water
q ONB
minimum heat flux for the onset of nucleate
com
compressor
conv
convective
e
evaporator
s
surface exchange by length unit
m-1)
f
liquid refrigerant
Se
heat transfer area of the evaporator [m2]
g
gas refrigerant
i
inner
o
outer
ONB
onset of nucleate boiling
boiling [W Re
m-2)
Reynold number
T
temperature
u
controllable inputs vectors
[m2
[oC]
U
overall heat-transfer coefficient [W
V
volume [m3]
(m2 K) -1]
r sat
refrigerant
Greek symbols
wt
water
w all
Wall metal
heat transfer coefficient [W (m2 K) -1]
this critical time constant is considered to change with the variation of the partial load and two-phase region heat transfer coefficient (Tahat et al., 2001). The dynamic behavior of oscillations was also explained by the strong non-linear course of the
saturation
heat transfer coefficient (Gruhle and Isermann, 1985; Lars, 1999). Chen et al. (2002) demonstrated the influence of heat transfer coefficient of evaporator on the DS. Considering that the variation of heat transfer coefficient may actually
Page 2 of 17
be an integrated result of density-wave oscillation and other oscillations like pressure drop oscillation, Liang et al. (2011) investigated the two-phase flow oscillations of R-22 in a refrigeration system with a horizontal straight tube evaporator. Among these parameters, the heat transfer coefficient impacts the control behavior more Reviewer #1: Response 2 directly and makes the extrapolation of the proposed model to another EEV-controlled refrigeration system more easily. Therefore, this paper developed a dynamic simulation model of the EEV-controlled refrigeration system where heat exchangers adopted the lumped parameter model with the moving-boundary hypothesis. The heat transfer characteristic is used to divide the evaporator into the stable and unstable conditions, and a matching switch criteria is presented. The dynamic model could be used to simulate the oscillation of system parameters by the switch between the stable and unstable condition. This study provides the basis for further developing advanced control algorithm which will ensure that an EEV-controlled system refrigeration system has better stability. The rest of this paper is organized as follows. In Section 2, the component model is presented. Section 3 shows the unstable condition and switching criteria. In Section 4, the experimental test facility is described. Section 5 shows the results and discussion. Finally, in Section 6, the main conclusions of the work are summarized.
2. Component model In this section the modeling of the refrigeration system consisting of the variable speed compressor, EEV, evaporator and condenser is presented. Firstly, the sub-models of four components are developed and then integrating them together. It can be observed that the dynamics of the compressor and EEV are much faster than that of the two heat exchangers. Thus a common assumption is that the time dependency of the compression and expansion processes is not important and the compressor and EEV will be treated as static components. 2.1 Compressor The compressor sub-model provides the refrigerant mass flow rate suctioned from the evaporator and discharged to the condenser, mass flux in the evaporator (Gosney, 1982):
m com f com v vV com
C v ,v
v ( S t - aS t ((
pc
)
Cv , p
- 1))
pe
H com , o H com , i ( H os H com ,i ) / com where
C v , v and C v , p
(2) (3)
are calculated for the saturated
vapor at the evaporating pressure, a and S t are fitted to experimental data as a part of the compressor model validation exercise. 2.2 Expansion valve Liquid refrigerant flowing through an expansion valve can simply be modeled as the following orifice equation: m v Ave C v
c , o ( Pc - Pe )
H v ,i H v ,o
(4) (5)
where Ave Ao Av / 100 is the effective passage flow area, Av is the valve opening[%], Ao is the nominal orifice cross-sectional area, and assuming the enthalpy value is constant in expansion process. 2.3 Evaporator The heat exchanger is divided into control volumes or zones (shown in Fig. 1). According to moving-boundary method, interface boundaries between zones vary over time and are tracked by the model. The governing equations are derived by considering conservation of mass and energy in each control volume. Several assumptions are made in order to reduce these equations: (1) One-dimensional flow for refrigerant in tubes. (2) Pressure drops are assumed to be negligible. (3) Thermodynamic equilibrium is assumed for two-phase fluid. (4) The refrigerants superheated vapour region are assumed as homogeneous flow. (5) The dynamic storage of mass and energy in superheated vapour region is neglected. (6) Exchangers are fully insulated. (7) The water of the heat exchangers is uniform. On the basis of the above, a set of partial derivative equations based on mass energy balance conservation principle can be expressed as:
(1) A
( ) t
(m ) z
0
(6)
Page 3 of 17
A
( h)
t
(m h) z
H e , i , T ch w ] , Eq. (9), (10), (11) can be expressed as s i i (T w T r )
(7)
the matrix equation: 1
(C p A ) w
(T w ) t
x e A e fe (x e , u e ) s o o (T w t T w ) s i i (T w T r )
It can be observed that in evaporator, two-phase refrigerant dominate the major portion of the heat exchangers. Thus a simplified model of the evaporator was proposed by LEDUCQ et al. (2002). These generalized governing equations are applied only on two-phase region, expressed as following Eqs. (9)–(11). Refrigerant mass balance: Ae L e ,1
d e ,1 dt
Ae ( e ,1 e , g )
d L e ,1 dt
m e ,i m e , g
(9)
d e , f he , f dPe
(1 e )
( e , f h e , f e , g h e , g )
d L e ,1 dt
d e , g he , g dPe
e 1)
d Pe dt
Ae (1 e )
m e , i he , i m e , g he , g
dt
e , i1 d e , i (T er 1 T e , w ) e , o 1 d e , o (T ch w T e , w )
a e ,11 Ae a e ,21 0
a e ,12 a e ,22 0
0 0 a e ,33
(17)
(11)
The condenser behaviour is modelled by dividing the heat exchanger into three zones: superheat region, two-phase region and subooled region. Eq. (6), (7) and (8) are applied on all these regions to obtain the condenser model (He, 1996). 1
Static balance equations for the superheated region are written as Eqs. (12)–(14): Refrigerant mass balance: m e ,o m e , g
(12)
Refrigerant energy balance equation in steady state: m co m C p , g d T U e , 2 D e , i d l (T ch w - T e )
(13)
Integrating energy balance equation along the superheated zone, the degree of refrigerant superheat at the evaporator outlet (DS) can be expressed as Eq. (14):
D S (Tchw Te )(1 e
( Le Le ,1 )
U e ,2 d e ,i m com C p , g
)
(16)
2.4 Condenser
Tube wall energy balance: d Te , w
m e ,i m e , o f e ( x e , u e ) m e , i H e , i m e , o H e , g e , i1 d e , i L e ,1 (Te , w Te ) e , i1 d e , i (Ter 1 Te , w ) e , o 1 d e , o (Tchw Te , w )
(10)
e , i 1 d e , i L e ,1 (T e , w T e )
cu C p , cu Ae , cu
Detailed expressions of the evaporator are embodied as Eqs. (16)–(17).
The indication expressions of A e are shown in Table 1. The correlations used in the model for estimating the heat transfer coefficients in the heat exchangers are detailed in Table 2.
Refrigerant energy balance: Ae L e ,1 (
(15)
(8)
(14)
Defining the state variables as x e [ Pe , L e ,1 , T w , e ] and the vector of inputs as u e [ m e , i , m e , o ,
x c A c fc (x c , u c )
(18)
x c [ L c ,1 , L c , 2 , Pc , H c , o , T c , w 1 , T c , w 2 , T c , w 3 ] a c ,11 a c , 21 a c ,31 A c a c , 41 a c ,51 0 a c ,71
0
a c ,13
0
0
0
a c , 22
a c , 23
0
0
0
a c ,32
a c ,33
a c ,34
0
0
a c , 42
a c , 43
0
0
0
0
0
0
a c ,55
0
0
0
0
0
a c ,66
a c ,72
0
0
0
0
(19) 0 0 0 0 0 0 a c ,77
u c [ m c , i , m c , o , H c , i , Tclw ] m c , i H c , i m c , i H c , g c , i 1 d c , i L c ,1 (Tc , w 1 Tc ,1 ) m H m c , o H c , f c , i 2 d c , i L c , 2 (T c , w 2 T c ) c ,o c , g m c , o H c , f m c , o H c , o c , i 3 d c , i L c ,3 (Tc , w 3 Tc ,3 ) m c ,i m c ,o fe (x e , u e ) d (T Tc ,1 ) c , o d c , o (Tclw T c , w 1 ) c ,i1 c ,i c , w1 e , i1 d e , i (Ter 1 Te , w ) c , o d c , o (Tclw Tc , w 2 ) d (T T ) d (T T ) e,w c ,o c ,o clw c ,w 3 e , i1 e , i er 1
(20)
(21)
(22)
Page 4 of 17
rcrit 0 .3 1 0 (1 .5 - Pe / 2 .2 1 0 ) -6
Expressions of all elements in matrix A c are given in Table 3. The condenser and evaporator mean void fraction are fitted to experimental data, being both considered to be equal to 0.8.
Fig.2 portrays the solution procedure of the heat exchanger. The two heat exchangers are combined with the variable speed compressor and the EEV which can provide the time varying inputs. In Fig.2, the time derivative of the state vector is a function ( A -1 f ) of the state vector ( x ), the model inputs, and
6
.
During practical operation, the changes of the EEV opening and compressor speed always increase the refrigerant mass in the evaporator, which causes that the given heat flux is relatively low. When q ONB is higher than the output heat flux of the evaporator ( q ), nucleus boiling will disappears and the vapor is generated by evaporation only on the interface by convective mechanism. Then, it is difficult to keep the dry-out point inside the evaporator because of the disappearance of nucleus boiling and the deterioration of the heat transfer.
numerous physical parameters. The function ( A -1 f ) can be acquired by Eqs. (15) for the evaporator and Eqs. (18) for the condenser. The solution is advanced in time by integrating the derivative from known initial conditions ( x(0) ). The heat exchanger model also provides the refrigerant pressure and outlet refrigerant enthalpy to the compressor and EEV.
3. Unstable condition and switch criteria In saturated flow boiling, the heat transfer coefficient at the nucleate boiling are very high because of the nature of nucleate boiling, and the nucleate and convective components are superimposed by a very complex mechanism. Therefore, the nucleate boiling in evaporator must be activated in order to maintain high heat change efficiency. Due to the surface tension forces, a given amount of heat flux ( q ONB ) is needed for the onset of nucleate boiling. Predictive equations for q ONB , which depends upon a number of operating parameters, are developed from bubble dynamics theory. The following equation was used by Steiner and Taborek, (1992) and is recommended for use with this method: qONB
2 T sat conv rcrit v ( H g - H f )
In order to simulate the dynamic behavior of oscillation, the idea is to divide the evaporator into nucleate and pure convective boiling condition, shown as Fig. 3. Through the switch between two heat transfer conditions, the step change of the heat transfer coefficient can be seen as a strong disturbance to the controller and the cause of control instability. Therefore, the heat transfer coefficient in the two-phase region is expressed as the sum of the convective and the nucleate boiling terms (Kandlikar, 1990): e , i1 C 1C o conv C 3 Bo F fl conv C2
C4
(24)
being C 1 C 4 constant parameters, Co convection number, Bo boiling number, F fl fluid dependent parameter, and conv the liquid convective heat transfer coefficient, computed with the Dittus-Boelter correlation. When the nucleate boiling condition is activated, both of the nucleate and convective boiling term in Eq. (24) can be used, which stated the heat transfer of the evaporator is at stable condition. However, if the pure convective boiling condition is activated, the nucleate boiling term in Eq. (24) is zero and the equation remains valid for convective boiling only: e , i1 C 1C o conv C2
(25)
(23)
where is the surface tension force of R22. It seems that the constant critical radius in Steiner-Taborek model is not applicable for the present condition. Analysis shows that a better fitting to the experimental data can be obtained by multiplying the constant value by a correction factor, and the modified critical radius
The activation of the pure convective boiling condition deteriorates the heat transfer and leads to a sudden descent of the DS, which reduces the safety and efficiency of the refrigeration system. Therefore, the pure convective boiling condition can be seen as unstable condition, and the DS at critical nucleate boiling condition is known as minimal stable value (MSV). Fig. 4 shows the events that trigger the change
Page 5 of 17
of the heat transfer condition. The heat flux discriminant is used when evaporator works at the nucleate boiling condition, and the pure convective boiling condition is activated when Eq. (26) becomes true: qONB q
(26)
where q denotes a threshold that specifies the maximum value of q O N B for the onset of nucleate boiling: q
Q
(27)
Se
being Q the output cooling capacity of the evaporator, S e the heat transfer area of the evaporator ( S e 3 .0 2 m 2 , manufacturer data for the evaporator). It is found that when Eq. (26) is satisfied, the disturbance on q ONB induced by the discontinuous decrease of the heat transfer coefficient leads to the switch data oscillation. In the open-loop simulation (abrupt changes of the compressor frequency or EEV opening), the switching oscillation will eventually vanish. However, the repeating switch between two conditions will not disappear spontaneously during the simulation embedded control algorithm, which causes the pure convective boiling condition cannot be activated normally. In order to expand the application of the proposed model, the DS at the outlet of evaporator is used to judge whether the nucleate boiling condition is activated again when Eq. (28) becomes true: DS MSV
(28)
where MSV is the DS at critical nucleate boiling condition at which q O N B is equal to q and the nucleation is just activated. Therefore, it is reasonable to use the MSV as the dividing point between two different heat transfer conditions. In addition, a dead zone technology has been adopted in order to provide a smooth switching progress during the simulation. For instance, if q O N B q , the pure convective boiling condition is active. Then, the DS discriminant will not be judged in the next few seconds. Fig. 6 shows the flow diagram for dynamic simulation of the evaporator model. Considering that the refrigeration system is slow in response
speed, the simulation step size should be as large as possible without affecting the solution precision. On the other hand, the heat transfer coefficient is discontinuous with the change of heat transfer condition. Discontinuity can be troublesome in the case of large simulation step size, and it may lead to chattering or convergence problems. A dead zone method has been adopted to provide a smooth switching progress during the simulation. Therefore, the variable time-step Runge-Kutta solver is used in order to quicken simulation when ensuring the stability of the proposed model. At the beginning of simulation, the evaporator works at the nucleate boiling condition. For a given time step,#3: the coefficients are first calculated with Reviewer Response 9 the input parameters. According to the heat transfer condition at the preceding time step, the Eq. (26) or Eq. (28) would be used to judge which heat transfer condition in the evaporator is activated at the current time step. When the heat transfer condition is chosen, the thermo-dynamic properties of R22 are looked up from REFPROP7. At the end of the time step, the output parameters of differential form are obtained with matrix calculation and are calculated after integral over time. The output data of parameters are then deposited in storage database. Proceed down the next time step with the same way till the appointed time has been attained.
4. Experimental work 4.1. Experimental facility The experiment to validate the model presented above is carried out in a liquid-to-liquid chiller whose schematic diagram is shown in Fig. 5. The experimental refrigeration system is able to work with the refrigerants R134a and R22, but the refrigerant used for the tests presented here is R22. The compressor is of hermetic variable speed reciprocating type and the expansion device is a pulse-width-modulation (PWM) electric expansion device with variable opening from 0 to 100%. The evaporator and condenser are both double-pipe exchanger heat exchangers, connected to two distinct temperature and flow controlled secondary heat transfer loops. Pure water is used for the evaporator and condenser. The cooling load system regulates the chilled water temperature using immersed electrical resistances, and the nomial output cooling capacity is 35kW. The thermodynamic states of the refrigerant
Page 6 of 17
and secondary fluids at each point are calculated using measurements of pressure and temperature at the inlet and outlet of each element of the facility. These measurements are obtained using Immersion type-T thermocouples (measurement uncertainty of ±0.2 o C ) and pressure transducers (measurement uncertainty of ±0.1%). The refrigerant mass flow rate is measured by a Coriolis-effect mass flow meter with a certified accuracy within ±0.25% of the reading, and the secondary fluids mass flow rates are measured with volumetric flow meters, introducing a maximum error of ± 0.33%. Furthermore, a capacitive sensor is installed to obtain the compressor rotation speed, with a maximum error after calibration of ±1 %. A control and data acquisition system is used for monitoring the experimental variables and for setting the compressor speed and valve opening. 4.2. Identifying the output heat flux model In the refrigeration system, both of the variations of compressor speed ( f ) and EEV opening ( VO ) may effect on the cooling capacity and output heat flux. The disturbance of cooling capacity ( Δ Q ) caused by the variation of the compressor speed ( Δf ) and EEV opening ( ΔVO ) can be represented using a first-order transfer function model (Li et al., 2009). Therefore, the disturbance of output heat flux ( Δq ) is expressed as follows: Δq=G 1 (s) Δf G 2 (s) ΔV O
where
G1
and
G2
(29)
are both first-order transfer
function model: G (s)
K τs+1
G1
and
G2
can be
estimated. Considering that the change of EEV opening has little effect on the output heat flux, G 2 is approximated at zero for simplicity.
G 1 (s)=
Δq
=
0.12
Δ f 32s+1 Δq G 2 (s)= 0 ΔVO
(31) (32)
5. Results and discussion In this section the validation of the model using experimental measurements and correlation analysis are presented. To this end, two sets of experiments have been undertaken: (a) sudden change of EEV opening and compressor frequency, and (b) oscillation of the system when the DS is regulated by P–I controller. Based on the compressor frequency, the output heat flux can be obtained with Eq. (31). The evaporating pressure, condensing pressure and the DS at evaporator outlet are chosen for validation. Finally, the instability of system when the system is affected by external disturbance, is simulated. 5.1. Sudden change of EEV opening and compressor speed The comparison between simulation and experimental results of dynamic process are depicted in Fig. 8 and 9. The test is initiated when the steady-state condition is achieved for the following operating conditions in both tests: EEV opening of 40%, compressor speed of 50 Hz.
(30)
where K represents the static, τ the time constant and s the Laplace operator. In order to identify the parameters of G 1 and G2
Matlab, the parameters of
, two sets of tests have been performed. At the
beginning, fixing the compressor speed at 45Hz and EEV opening at 40% in both tests. Once in steady state, the disturbances of EEV opening (40–44) and compressor speed (45–50) are imposed respectively. The heat flux response due to the sudden reduction of compressor speed is presented in Fig. 7. Using the System Identification Toolbox in
Independent disturbances of EEV opening and compressor frequency are imposed following the sequence shown in Table 4. It is noteworthy that the disturbances of EEV opening in Fig. 8 cause a decrease of the DS at t=200s, and the DS quickly falls to zero at t=600s. The similar results can be seen in Fig. 9 at t=100s and 300s, respectively. During the tests, the change in DS is influenced by two factors. One is the total mass of refrigerant in evaporator, and the other is the heat transfer between the refrigerant and chilled water. The first step increase of the EEV opening and decrease of the compressor speed imply more refrigerant mass
Page 7 of 17
in the evaporator, which results in a longer two-phase region for evaporator, and thus the DS drops. Since the second disturbances are just a little bit stronger than the first, the differences of the the heat transfer between the refrigerant and chilled water. The comparison results show that a good matching between the simulated and experimental data is achieved not only in the nucleate boiling condition but also the pure convective boiling condition. Fig. 8.b and 9.b show the relationship between the actual heat flux and q O N B as the changes of EEV opening and compressor frequency, respectively. The output heat flux is a constant value corresponding to the compressor frequency of 50Hz in Fig. 8.b, and the variation of the output heat flux caused by the change of the compressor frequency is obtained with Eq. (31) in Fig. 9.b. Based on Eq. (23), the given amount of heat flux needed for the onset of nucleate boiling ( q O N B ) can be calculated with the changing DS. At the beginning, the DS is relatively high, and the output heat flux is clearly higher than the value needed. Therefore, the nucleate boiling condition is activated. When the DS reduces at its MSV, the output heat flux is lower than q O N B , which cause that the pure convective boiling condition is activated according to Eq. (26). As the DS is higher than MSV, the nucleate boiling condition is activated again based on Eq. (28) and deterioration of the heat transfer finally disappear. This result provides a further understanding of the stability of control loop. The difference between the output heat flux and q O N B can be seen as the stability margin, which should be positive in order to activate the nucleate boiling condition. It is seen from Fig. 8.b and 9.b that q O N B increases with the increase of EEV opening and decrease of compressor frequency. And the decrease of compressor frequency also causes the decrease of output heat flux. In practical operation, in response to system load changes, compressor speed may be continuously altered to adjust output cooling capacity to match the changes in cooling load. Due to the interaction, the EEV opening should be regulated systematically in order to drive the DS towards its reference. The analysis indicates that the decrease of compressor speed leads to the decrease of the stability margin, which makes the control of DS more difficult. begins. It is worth noting that the steady-state MSV (after the third switch) is almost the reference of
refrigerant mass in the evaporator are little. Therefore, the sudden reductions in the DS in Fig. 8.a and 9.a are due to the sudden deterioration of
5.2. Oscillation of the control loop The oscillations of system parameters when a decreasing step change on the DS is imposed, are depicted in Fig. 10. The compressor speed is 50 Hz, and the EEV opening is derived by P–I algorithm to drive the DS towards the desired references of 8 o C . For the EEV controller, a set of constant P–I values, 1.0 and 0.5, respectively, is adopted. The temperatures of the chilled and cooled water are 12oC and 30oC, respectively. The reference of the DS is step decreased to 5.6 o C at t=50s, and kept on this value. It can be seen that the oscillations of the DS near the reference can be observed after about t=100s, and the system parameters also fluctuate due to the oscillation of the EEV opening. In Fig. 10, a switch from the nucleate boiling condition to pure convective boiling condition can be seen at t=70s. Such an effect is induced by the large magnitude of the imposed change. Due to the over modulation, the DS is lower than MSV during the transient. Consequently, the pure convective boiling condition is activated, which implies the refrigerant in the evaporator cannot absorb enough heat for phase change, and then the DS drops suddenly. However, a further decline of the DS is rejected by the controller which closes the EEV to a smaller position in order to maintain the DS to its reference. As the DS is higher than MSV, the nucleate boiling condition is activated again and deterioration of the heat transfer finally disappears. However, when the DS is close to the reference of 5.6 o C at about t=150s, the pure convective boiling condition is triggered due to the insufficient of the output heat flux to support nucleate boiling, and an abruptly decrease of the DS occurs one more time. The root reason can be concluded that the new balance stability margin is too small, which causes the switch from nucleate to pure convective boiling condition is inevitable with the effect of the controller. Consequently, the periodical change of the heat transfer condition leads to the disturbance on the DS, the controller has to pursue the DS ceaselessly, and the oscillation of the control loop
the DS of 5.6 o C in the example. However, the oscillation of the system can be still emulated
Page 8 of 17
successfully because of the over modulation, as long as the new balance q O N B is close enough to the output heat flux. It is clear that actual fluctuations of the DS control loop and pressures of heat exchangers are consistent with the simulated ones. However, special attention should be paid to the large difference between the tested and simulated condensing pressure at about t=100s. The difference can be concluded that the change of EEV opening in simulation is obviously faster and larger than the actual experimental one within the time interval 50s
The system is first working under the same steady-state condition with the example in section 5.2, and the water temperature goes up 2.5 o C at t=50s. Such variation of the surrounding condition leads to the increase of q O N B and sudden drop of the DS. Albeit the disturbances on the DS have been recognized by the controller, the final stability margin is so small that the periodical switch between two heat transfer conditions cannot be eliminated.
6. Conclusion In this paper, A dynamic lumped parameter model for EEV-controlled refrigeration system has been developed. Two types of heat transfer condition in the evaporator are defined according to different heat transfer mechanisms of refrigerant. Based on the theory of bubble dynamics, a novel technique for the switch between different heat transfer conditions is presented. The difference between the output heat flux and q O N B is seem as Reviewer #1: Response 4 the stability margin. The analysis provides detailed insights to the operational stability due to the change of heat transfer mechanism of refrigerant. In addition, the dynamic behavior of oscillation of system parameters is simulated successfully by the switch between stable and unstable condition. The predicted results in comparison to the corresponding experimental data fall into a reasonable range. The discontinuous changes of the heat transfer coefficient can be used to quantify the system’s inherent characteristic which causes the instability, which is benefit to the design of controller. This study provides the basis for further developing advanced control algorithm which will ensure that an EEV-controlled system refrigeration system has better stability.
Acknowledgements
5.3. Sudden change of the temperature of cooling water
This work is supported by the Natural Science Foundation of China under Grants 61403274.
Finally, the effect of sudden change in the temperature of cooling water on the behavior of the system is simulated by the proposed model, and the transient responses of system parameters when the step changes are imposed on the cooling water temperature, are shown in Fig. 11.
References: Aprea, C., Mastrubllo, R., 2002. Experimental evaluation of electronic and thermostatic expansion valves performance using R22 and R407C. Appl Therm Eng 22, 205–218. Broersen, P.M.T., Van der Jagt M.F.G., 1980. Hunting
Page 9 of 17
of evaporator controlled by a thermostatic expansion valve. Trans ASME, J Dynam Syst Meas Control 102, 130–5. Chen, W., Chen, Z.J., Zhu, R.Q., Wu, Y.Z., 2002. Experimental investigation of a minimum stable superheat control system of an evaporator. Int. J. Refrigeration 25, 1137–1142. Chen, Y., Deng, S., Xu, X., Chan, M., 2008. A study on the operational stability of a refrigeration system having a variable speed compressor. Int. J. Refrigeration 31, 1368–1374. Gosney, W.C., 1982. Principles of Refrigeration, Cambridge. Gruhle, W.D., Isermann, H., 1985. Modeling and control of a refrigerant evaporator. J. Dyn. Syst. Meas. Control 107, 235-240. He, X.D., 1996. Dynamic Modeling and Multivariable Control of Vapor Compression Cycles in Air Conditioning Systems, PhD Thesis, Massachusetts Institute of Technology, America. Ibrahim., G.A., 2001. Effect of sudden changes in evaporator external parameters on a refrigeration system with an evaporator controlled by a thermostatic expansion valve. Int. J. Refrigeration 24, 566–576. Koury, R.N.N., Machado, L., Ismail, K.A.R., 2001. Numerical simulation of a variable speed refrigeration system, Int. J. Refrigeration 24 (2), 192–200. Kandlikar, S.G., 1990. A general correlation for saturated two-phase flow boiling heat transfer inside horizontal and vertical tubes, Journal of Heat Transfer 112, 219–228. Lars B.K., 1999. Novel electronic high reliability valve principle for control of direct expansion. In: Proceedings of 20th international congress of refrigeration, Sydney, IIR/IIF, 445. Leducq, D., Guilpart, J., Trystram, G., 2003. Low order dynamic model of a vapor compression cycle for process control design, Journal of Food Process Engineering 26, 67–91. Lenger, M.J., Jacobi, A.M., Hrojak, P.S., 1998. Superheat Stability of an Evaporator and Thermostatic Expansion Valve Air Conditioning and Refrigeration Center. College of Engineering. University of Illinois at Urbana-Champaign. Leonardo, C.S., Christian, J.L.H., Alexandre, T.N., 2010. Assessment of the controlling envelope of a model-based multivariable controller for vapor compression refrigeration systems. Appl. Therm. Eng. 30, 1538–1546. Li, X.Q., Chen, J.P., Chen, Z.J., et al., 2004. A new method for controlling refrigerant flow in automobile air conditioning. Applied Thermal Engineering 24, 1073–1085. Li H., Jeong, S.K., You, S.S., 2009. Feedforward
control of capacity and superheat for a variable speed refrigeration system, Applied Thermal Engineering 29, 1067–1074. Liang, N., Shao, S.Q., Xu, H.B., Tian, C.Q., 2010. Instability of refrigeration system–A review, Energy Conversion and Management 51, 2169–2178. Liang, N., Shao, S.Q., Tian, C.Q., Yan, Y.Y., 2011. Two-phase flow instabilities in horizontal straight tube evaporator, Applied Thermal Engineering 31, 181–187. Mithraratne, P., Wijeysundera N.E., 2001. An experimental and numerical study of the dynamic behaviour of a counter-flow evaporator. Int. J. Refrigeration 24, 554–65. Rasmussen, H., Larsen, L.F.S., 2011. Non-linear and adaptive control of a refrigeration system. IET Control Theory Appl 5, 364–378. Shah, M.M., 1979. A general correlation for heat transfer during film condensation inside pipes. Int. J. Heat Mass Transfer 22(4), 547–556. Steiner, D., Taborek, J., 1992. Flow boiling heat transfer in vertical tubes correlated by an asymptotic model. Heat Transfer Engineering 13 (2), 43–69. Tahat, M.A., Ibrahim, G.A., Probert, S.D., 2001. Performance instability of a refrigerator with its evaporator controlled by a thermostatic expansion-valve. Applied Energy 70, 233-249. Yasuda, H., Toufer, S., 1983. Simulation model of a vapour compression refrigeration system. ASHRAE Trans, 408–24. Zahn, W.R., 1964. A visual study of two-phase flow while evaporating in horizontal tubes. Trans ASME, J Heat Transfer 86, 417–29. Zhang, Q., Canova, M., 2015. Modeling and output feedback control of automotive air conditioning system. Int. J. Refrigeration 58, 207–218.
Page 10 of 17
Superheat region
Two-phase region
Two-phase region
q Superheat region
Superheat region
SV M S> D
>q
(1) evaporator
B N O
Subcooled region
Two-phase region
(2) condenser
Fig. 1. Schematic diagrams of evaporator and condenser.
u
x A -1f
x
Superheat region
Fig. 4. Switching between different heat transfer conditions in the evaporator.
Physical Parameters
Inputs
Two-phase region
Cooled water side
x
EEV
x(0)
Refrigerant side
Condenser Variable-speed compressor
Evaporator
Fig. 2. Solution procedure of the heat exchanger Chilled water side
Chilled water tank (1m 1m 1m)
Variable-speed supply pump
Fig.5. The schematic diagram of the experimental plant. Two-phase region
Superheat region
(1) evaporator in nucleate boiling condition
Two-phase region
Superheat region
(2) evaporator in pure convective boiling condition
Fig. 3. Evaporator in nucleate and pure convective boiling conditions.
Page 11 of 17
Start
Inputs Initial values
Input
Calculate coefficients
Yes
No Nucleate boiling?
Look-up
Look-up
Te , e, g , e, f , H e, g , H e, f , e,i1 (eqn.25)
Te , e, g , e, f , H e, g , H e, f , e,i1 (eqn.24)
Matrix Calculation
Matrix Calculation
dPe , dLe,1 , dTw,e
dPe , dLe,1 , dTw,e
Check for switch to convective boiling (eqn.26)
Check for switch to nucleate boiling (eqn.28)
Integral process
Integral process
Next time step Output
No
Data storage
Appointed time?
Yes End
Fig. 6. Flow diagram for the dynamic simulation program of the evaporator.
Page 12 of 17
Heat flux(kW m-2)
10.2 10 9.8 9.6 9.4 9.2 9 0
100
200
300
400
500
600
700
800
t(s)
Fig. 7. The heat flux response due to the change of compressor speed
Page 13 of 17
a
b
8
Heat flux(kW·m-2)
Degree of superheat(oC)
10
6 4 2 0 0
200
400
600
800
Output heat flux
10 9.5 qONB
9 8.5
1000
Data Model
10.5
0
200
400
t(s)
d
505
495
P c(MPa)
P e(kPa)
c
485
475
0
200
400
600
800
1000
600
800
1000
t(s)
600
800
1.41
1.39
1.37
1.35
1000
0
200
400
t(s)
t(s)
Fig. 8. Validation for the transient response to the change of the EEV opening: a. degree of superheat, b. relationship between output heat flux and qONB, c. evaporating pressure, and d. condensing pressure.
a
b
10
Output heat flux
Heat flux(kW·m-2)
Degree of superheat(oC)
8 6 4 2 0
0
Data Model
10.5
100
200
300
400
10 9.5 9 qONB
8.5
500
0
100
200
t(s)
c
d
520
P c(MPa)
P e(kPa)
400
500
300
400
500
1.41 1.4
510 500 490 480
300
t(s)
1.39 1.38 1.37
0
100
200
300
t(s)
400
500
1.36
0
100
200
t(s)
Fig. 9. Validation for the transient response to the change of the compressor speed: a. degree of superheat, b. relationship between output heat flux and qONB, c. evaporating pressure, and d. condensing pressure.
Page 14 of 17
b
55 50
40
495 485
35 30 0
Model Data
515 505
45
Pe (kPa)
EEV opening(%)
a
100
200
300
400
475
500
0
100
200
1.6
c Pc(MPa)
1.5 1.4 1.3 1.2
0
100
200
300
t(s)
300
400
500
300
400
500
t(s)
400
500
d
9
Degree of superheat(oC)
t(s)
7
5
3 0
100
200
t(s)
Fig. 10. System oscillation analysis: a. EEV-opening, b. evaporating pressure, c. condensing pressure, and d. degree of superheat.
Page 15 of 17
EEV opening(%)
45
40
35 0
100
200
300
400
500
600
700
400
500
600
700
t(s) Degree of superheat(oC)
7 6.5 6 t(s)
5.5 5 0
100
200
300 t(s)
Fig. 11. System oscillation when the system is disturbed by the abrupt change of external parameter
Table 1. Indication expressions of coefficient matrix for the evaporator. Symbol
Indication expression
a e ,11
é d r e, f d r e,g ù ú Ae l e ,1 ê (1 - g e ) +g e dPe dPe ú ëê û
a e ,12
Ae (1 - g e )( r
a e , 21
é ù d r e , f he , f d r e , g he , g Ae l e ,1 ê (1 - g e ) +g e - 1ú êë ú dPe dPe û
a e , 22
Ae (1 - g e )( r
a e ,33
C p , cu r
- r
e, f
e, f
e,g
)
he , f - r
e,g
he , g )
A cu e , cu
Table 2. Calculation expressions of coefficients in heat exchangers. Symbol
Expression
e , i1
S. G. Kandlikar equation (Kandlikar, 1990): i C1Co 2 conv C 3 Bo 4 F fl conv C
c , i1 , c , i 3
C
Dittus-Boelter equation: Nu 0.023 Re i
e ,i 2 , c ,i 2
0.8
0.4
Pr ,
Nu di
Traviss, Baron and Rohsenow equation: Nu
Re
0.9
Pr F1
, i
F2
e,o , c ,o
Nu di
Shah equation(Shah, 1979): o Bw
0.8
(
Di
)
R
0.18
0.2
Do
,
De
B 1496 22T w t
Table 3. Indication expressions of coefficient matrix for the condenser. Symbol
Indication expression
a c ,11
Ac r
a c ,13
d r c ,1 1 é ïì Ac l c ,1 í ( H c ,1 - H c , g ) + êr dPc 2 êë ïî
c ,1
( H c ,1 - H c , g )
a c , 21
Ac r
c, f
a c , 22
Ac r
g c ( H c,g - H c, f ) c, f
a c , 23
ìï Ac l c ,2 í (1 - g c ) r ïî
a c ,31
-
a c ,32
-
a c ,33
Ac L c ,3 (
1 2 1 2
c ,1
+ ( H c ,1 - H c , g )
ù dH c , g ïü c ,1 ú - 1ý dH c úû dPc e ïþ
dr
( H c,g - H c, f )
dH c , f c, f
dPc
é d r c,g + g c ê( H c , f - H c , g ) +r êë dPc
Ac r
c, f
( H c , f - H c ,o )
Ac r
c, f
( H c , f - H c ,o )
1 2
r
dH c , f c, f
c,g
dH c , g ù údPc úû
üï 1ý ïþ
- 1)
dPc
Page 16 of 17
1
a c ,34
2
Ac L c ,3 r
c, f
a c , 41
Ac ( r
a c , 42
Ac g c ( r
a c , 43
ìï Ac í L c ,2 îï
a c ,51
Ac , cu C p , cu r
cu
(T c , w 1 - T c , w 2 )
a c ,55
Ac , cu C p , cu r
cu
L c ,1
a c , 66
Ac , cu C p , cu r
cu
a c , 71
Ac , cu C p , cu r
cu
(T c , w 3 - Tc , w 2 )
a c , 72
Ac , cu C p , cu r
cu
(T c , w 3 - Tc , w 2 )
a c , 77
Ac , cu C p , cu r
cu
L c ,3
c ,1
- r c,g
c, f
- r
) c, f
)
é d r c, f d r c,g ê (1 - g c ) +g c êë dPc dPc
ù æ d r c ,1 1 d r c ,1 dH c , g ú + L c ,1 ç + ç úû 2 dH c dPc è dPc
ö üï ÷ý ÷ ø þï
Table 4. Changes imposed to the valve opening and compressor speed Time(s)
Av[%] fcom[Hz]
t<100
100
200
300
400
40 50
43 45
40 50
44 44
40 50
Page 17 of 17
S0140-7007(16)30231-6 http://dx.doi.org/doi: 10.1016/j.ijrefrig.2016.07.020 JIJR 3393
To appear in:
International Journal of Refrigeration
Received date: Revised date: Accepted date:
18-2-2016 31-5-2016 25-7-2016
Please cite this article as: Yujia Shang, Aiguo Wu, Xing Fang, Yuwen You, Dynamic simulation of electronic expansion valve controlled refrigeration system under different heat transfer conditions, International Journal of Refrigeration (2016), http://dx.doi.org/doi: 10.1016/j.ijrefrig.2016.07.020. 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.
Dynamic simulation of electronic expansion valve controlled refrigeration system under different heat transfer conditions Yujia Shanga,
Aiguo Wua,
Xing Fanga,
Yuwen Youb
College of Electric Engineering and Automation, Tianjin University, No. 92, Weijin Road, Nankai District, Tianjin 300072, China b College of Automation, Tianjin urban construction University, Tianjin 300384, China Correspondence should be addressed to Xing Fang: [email protected] a
Highlights
Based on bubble theory, evaporator is divided into two heat transfer conditions. A switch criteria between two heat transfer conditions is presented. We develop a dynamic model that can predict the instability of system control loop. The oscillation of the refrigeration system is studied using the proposed model.
Abstract: This work presents a new dynamic lumped parameter model that is able to simulate the oscillation of electronic expansion valve (EEV) controlled refrigeration system due to unsuitable control or abrupt change of external parameters. Based on bubble dynamics theory, the heat transfer characteristic is used to divide the evaporator into the stable and unstable conditions, and a matching switch criteria is presented. Together with the model of main components, the dynamic behavior of oscillation of system parameters is simulated successfully. Model validation against experimental data demonstrates the capabilities of the modeling approach in predicting the instability of system control loop. This study provides the basis for further developing advanced control algorithm which will ensure that an EEV-controlled system refrigeration system has better stability. Keywords: Refrigeration system, Simulation, Superheat, Control, Stability
1.Introduction In vapor compression refrigeration systems, the instability of system control loop is the phenomena of oscillation of all parameters such as the refrigerant flow rate, refrigerant pressures and temperatures at different points. The instability will lead lower safety, lower life-span and higher energy consumption (Liang et al., 2010; Mithraratne and Wijeysundera, 2001). Dynamic simulation can be used to study the instability of control loop and provide guidance for design and control strategy of refrigeration system. The simulation of the instability on refrigeration system control characteristics have been noted in some studies for
thermostatic expansion valve (TEV) controlled evaporator refrigeration systems (Yasuda and Toufer, 1983; Ibrahim, 2001; Tahat et al., 2001; Lenger et al., 1998). With the progressively application of EEV, there have been a number of reported studies on the instability in EEV-controlled refrigeration systems (Li et al., 2004; Aprea and Mastrubllo, 2002; Chen et al., 2008; Fallahsohi et al., 2010; Leonardo et al., 2010). Among these researches, a model-driven multivariable controller for vapor compression-EEV refrigeration systems was presented (Leonardo et al., 2010). Fallahsohi et al. (2010) developed a predictive functional control method to control the degree of evaporator superheat (DS) with an
Page 1 of 17
electronic expansion valve. However, the system oscillation can be still observed in these advanced control refrigeration systems when the DS setting is low enough. Therefore, there is a need for a model, with low computational cost, to simulate the dynamic behavior of oscillation for an EEV-controlled refrigeration system and this model should be easily used to develop advanced control algorithms. The key point of modeling is to define the unstable condition, using a macroscopical Nomenclature
parameter of the expansion valve or evaporator, which is the cause of instability. With the frequency-domain analysis method, several characteristic parameters of TEV, such as the thermal capacitance of the bulb and thermal conductance between the bulb and wall were found to be associated with instability (Broersen and Van der Jagt, 1980; Yasuda and Toufer, 1983; Ibrahim, 2001). For the evaporator, there is an evaporator critical time constant, below which the system operates in a stable manner. And
A
state-space matrices
A
exchange surface area[m2]
rcrit
coefficient in Traviss, Baron and Rohsenow equation critical bubble radius
Av
valve opening [%]
void fraction
Cp
specific heat [J(kg K)-1]
com
global efficiency of the compressor
Cv
coefficient of valve discharge
v
compressor volumetric efficiency
C v, p
specific heat at constant pressure [J(kg K)-1]
surface tension [N m-1]
density [kg m-3]
Subscripts pts
two-phase region of evaporator and
K)-1]
C v ,v
specific heat at constant volume [J(kg
d
Diameter [m]
DS
degree of refrigerant superheat Reviewer #3: Response 2
f
compressor frequency [Hz]
H
Enthalpy [kJ kg-1)
H os
outlet enthalpy in isentropic compression
1
superheat region of condenser 2
process [J kg-1]
superheat region of evaporator and two-phase region of condenser
L
length [m]
3
subcooled region of condenser
m
mass flow rate [kg s-1]
avg
average
P
pressure [k Pa)
c
condenser
Pr
Prandtl number
cu
copper
Q
heat transfer rate [W)
chw
chilled water
q
output heat flux [W m-2]
clw
cooled water
q ONB
minimum heat flux for the onset of nucleate
com
compressor
conv
convective
e
evaporator
s
surface exchange by length unit
m-1)
f
liquid refrigerant
Se
heat transfer area of the evaporator [m2]
g
gas refrigerant
i
inner
o
outer
ONB
onset of nucleate boiling
boiling [W Re
m-2)
Reynold number
T
temperature
u
controllable inputs vectors
[m2
[oC]
U
overall heat-transfer coefficient [W
V
volume [m3]
(m2 K) -1]
r sat
refrigerant
Greek symbols
wt
water
w all
Wall metal
heat transfer coefficient [W (m2 K) -1]
this critical time constant is considered to change with the variation of the partial load and two-phase region heat transfer coefficient (Tahat et al., 2001). The dynamic behavior of oscillations was also explained by the strong non-linear course of the
saturation
heat transfer coefficient (Gruhle and Isermann, 1985; Lars, 1999). Chen et al. (2002) demonstrated the influence of heat transfer coefficient of evaporator on the DS. Considering that the variation of heat transfer coefficient may actually
Page 2 of 17
be an integrated result of density-wave oscillation and other oscillations like pressure drop oscillation, Liang et al. (2011) investigated the two-phase flow oscillations of R-22 in a refrigeration system with a horizontal straight tube evaporator. Among these parameters, the heat transfer coefficient impacts the control behavior more Reviewer #1: Response 2 directly and makes the extrapolation of the proposed model to another EEV-controlled refrigeration system more easily. Therefore, this paper developed a dynamic simulation model of the EEV-controlled refrigeration system where heat exchangers adopted the lumped parameter model with the moving-boundary hypothesis. The heat transfer characteristic is used to divide the evaporator into the stable and unstable conditions, and a matching switch criteria is presented. The dynamic model could be used to simulate the oscillation of system parameters by the switch between the stable and unstable condition. This study provides the basis for further developing advanced control algorithm which will ensure that an EEV-controlled system refrigeration system has better stability. The rest of this paper is organized as follows. In Section 2, the component model is presented. Section 3 shows the unstable condition and switching criteria. In Section 4, the experimental test facility is described. Section 5 shows the results and discussion. Finally, in Section 6, the main conclusions of the work are summarized.
2. Component model In this section the modeling of the refrigeration system consisting of the variable speed compressor, EEV, evaporator and condenser is presented. Firstly, the sub-models of four components are developed and then integrating them together. It can be observed that the dynamics of the compressor and EEV are much faster than that of the two heat exchangers. Thus a common assumption is that the time dependency of the compression and expansion processes is not important and the compressor and EEV will be treated as static components. 2.1 Compressor The compressor sub-model provides the refrigerant mass flow rate suctioned from the evaporator and discharged to the condenser, mass flux in the evaporator (Gosney, 1982):
m com f com v vV com
C v ,v
v ( S t - aS t ((
pc
)
Cv , p
- 1))
pe
H com , o H com , i ( H os H com ,i ) / com where
C v , v and C v , p
(2) (3)
are calculated for the saturated
vapor at the evaporating pressure, a and S t are fitted to experimental data as a part of the compressor model validation exercise. 2.2 Expansion valve Liquid refrigerant flowing through an expansion valve can simply be modeled as the following orifice equation: m v Ave C v
c , o ( Pc - Pe )
H v ,i H v ,o
(4) (5)
where Ave Ao Av / 100 is the effective passage flow area, Av is the valve opening[%], Ao is the nominal orifice cross-sectional area, and assuming the enthalpy value is constant in expansion process. 2.3 Evaporator The heat exchanger is divided into control volumes or zones (shown in Fig. 1). According to moving-boundary method, interface boundaries between zones vary over time and are tracked by the model. The governing equations are derived by considering conservation of mass and energy in each control volume. Several assumptions are made in order to reduce these equations: (1) One-dimensional flow for refrigerant in tubes. (2) Pressure drops are assumed to be negligible. (3) Thermodynamic equilibrium is assumed for two-phase fluid. (4) The refrigerants superheated vapour region are assumed as homogeneous flow. (5) The dynamic storage of mass and energy in superheated vapour region is neglected. (6) Exchangers are fully insulated. (7) The water of the heat exchangers is uniform. On the basis of the above, a set of partial derivative equations based on mass energy balance conservation principle can be expressed as:
(1) A
( ) t
(m ) z
0
(6)
Page 3 of 17
A
( h)
t
(m h) z
H e , i , T ch w ] , Eq. (9), (10), (11) can be expressed as s i i (T w T r )
(7)
the matrix equation: 1
(C p A ) w
(T w ) t
x e A e fe (x e , u e ) s o o (T w t T w ) s i i (T w T r )
It can be observed that in evaporator, two-phase refrigerant dominate the major portion of the heat exchangers. Thus a simplified model of the evaporator was proposed by LEDUCQ et al. (2002). These generalized governing equations are applied only on two-phase region, expressed as following Eqs. (9)–(11). Refrigerant mass balance: Ae L e ,1
d e ,1 dt
Ae ( e ,1 e , g )
d L e ,1 dt
m e ,i m e , g
(9)
d e , f he , f dPe
(1 e )
( e , f h e , f e , g h e , g )
d L e ,1 dt
d e , g he , g dPe
e 1)
d Pe dt
Ae (1 e )
m e , i he , i m e , g he , g
dt
e , i1 d e , i (T er 1 T e , w ) e , o 1 d e , o (T ch w T e , w )
a e ,11 Ae a e ,21 0
a e ,12 a e ,22 0
0 0 a e ,33
(17)
(11)
The condenser behaviour is modelled by dividing the heat exchanger into three zones: superheat region, two-phase region and subooled region. Eq. (6), (7) and (8) are applied on all these regions to obtain the condenser model (He, 1996). 1
Static balance equations for the superheated region are written as Eqs. (12)–(14): Refrigerant mass balance: m e ,o m e , g
(12)
Refrigerant energy balance equation in steady state: m co m C p , g d T U e , 2 D e , i d l (T ch w - T e )
(13)
Integrating energy balance equation along the superheated zone, the degree of refrigerant superheat at the evaporator outlet (DS) can be expressed as Eq. (14):
D S (Tchw Te )(1 e
( Le Le ,1 )
U e ,2 d e ,i m com C p , g
)
(16)
2.4 Condenser
Tube wall energy balance: d Te , w
m e ,i m e , o f e ( x e , u e ) m e , i H e , i m e , o H e , g e , i1 d e , i L e ,1 (Te , w Te ) e , i1 d e , i (Ter 1 Te , w ) e , o 1 d e , o (Tchw Te , w )
(10)
e , i 1 d e , i L e ,1 (T e , w T e )
cu C p , cu Ae , cu
Detailed expressions of the evaporator are embodied as Eqs. (16)–(17).
The indication expressions of A e are shown in Table 1. The correlations used in the model for estimating the heat transfer coefficients in the heat exchangers are detailed in Table 2.
Refrigerant energy balance: Ae L e ,1 (
(15)
(8)
(14)
Defining the state variables as x e [ Pe , L e ,1 , T w , e ] and the vector of inputs as u e [ m e , i , m e , o ,
x c A c fc (x c , u c )
(18)
x c [ L c ,1 , L c , 2 , Pc , H c , o , T c , w 1 , T c , w 2 , T c , w 3 ] a c ,11 a c , 21 a c ,31 A c a c , 41 a c ,51 0 a c ,71
0
a c ,13
0
0
0
a c , 22
a c , 23
0
0
0
a c ,32
a c ,33
a c ,34
0
0
a c , 42
a c , 43
0
0
0
0
0
0
a c ,55
0
0
0
0
0
a c ,66
a c ,72
0
0
0
0
(19) 0 0 0 0 0 0 a c ,77
u c [ m c , i , m c , o , H c , i , Tclw ] m c , i H c , i m c , i H c , g c , i 1 d c , i L c ,1 (Tc , w 1 Tc ,1 ) m H m c , o H c , f c , i 2 d c , i L c , 2 (T c , w 2 T c ) c ,o c , g m c , o H c , f m c , o H c , o c , i 3 d c , i L c ,3 (Tc , w 3 Tc ,3 ) m c ,i m c ,o fe (x e , u e ) d (T Tc ,1 ) c , o d c , o (Tclw T c , w 1 ) c ,i1 c ,i c , w1 e , i1 d e , i (Ter 1 Te , w ) c , o d c , o (Tclw Tc , w 2 ) d (T T ) d (T T ) e,w c ,o c ,o clw c ,w 3 e , i1 e , i er 1
(20)
(21)
(22)
Page 4 of 17
rcrit 0 .3 1 0 (1 .5 - Pe / 2 .2 1 0 ) -6
Expressions of all elements in matrix A c are given in Table 3. The condenser and evaporator mean void fraction are fitted to experimental data, being both considered to be equal to 0.8.
Fig.2 portrays the solution procedure of the heat exchanger. The two heat exchangers are combined with the variable speed compressor and the EEV which can provide the time varying inputs. In Fig.2, the time derivative of the state vector is a function ( A -1 f ) of the state vector ( x ), the model inputs, and
6
.
During practical operation, the changes of the EEV opening and compressor speed always increase the refrigerant mass in the evaporator, which causes that the given heat flux is relatively low. When q ONB is higher than the output heat flux of the evaporator ( q ), nucleus boiling will disappears and the vapor is generated by evaporation only on the interface by convective mechanism. Then, it is difficult to keep the dry-out point inside the evaporator because of the disappearance of nucleus boiling and the deterioration of the heat transfer.
numerous physical parameters. The function ( A -1 f ) can be acquired by Eqs. (15) for the evaporator and Eqs. (18) for the condenser. The solution is advanced in time by integrating the derivative from known initial conditions ( x(0) ). The heat exchanger model also provides the refrigerant pressure and outlet refrigerant enthalpy to the compressor and EEV.
3. Unstable condition and switch criteria In saturated flow boiling, the heat transfer coefficient at the nucleate boiling are very high because of the nature of nucleate boiling, and the nucleate and convective components are superimposed by a very complex mechanism. Therefore, the nucleate boiling in evaporator must be activated in order to maintain high heat change efficiency. Due to the surface tension forces, a given amount of heat flux ( q ONB ) is needed for the onset of nucleate boiling. Predictive equations for q ONB , which depends upon a number of operating parameters, are developed from bubble dynamics theory. The following equation was used by Steiner and Taborek, (1992) and is recommended for use with this method: qONB
2 T sat conv rcrit v ( H g - H f )
In order to simulate the dynamic behavior of oscillation, the idea is to divide the evaporator into nucleate and pure convective boiling condition, shown as Fig. 3. Through the switch between two heat transfer conditions, the step change of the heat transfer coefficient can be seen as a strong disturbance to the controller and the cause of control instability. Therefore, the heat transfer coefficient in the two-phase region is expressed as the sum of the convective and the nucleate boiling terms (Kandlikar, 1990): e , i1 C 1C o conv C 3 Bo F fl conv C2
C4
(24)
being C 1 C 4 constant parameters, Co convection number, Bo boiling number, F fl fluid dependent parameter, and conv the liquid convective heat transfer coefficient, computed with the Dittus-Boelter correlation. When the nucleate boiling condition is activated, both of the nucleate and convective boiling term in Eq. (24) can be used, which stated the heat transfer of the evaporator is at stable condition. However, if the pure convective boiling condition is activated, the nucleate boiling term in Eq. (24) is zero and the equation remains valid for convective boiling only: e , i1 C 1C o conv C2
(25)
(23)
where is the surface tension force of R22. It seems that the constant critical radius in Steiner-Taborek model is not applicable for the present condition. Analysis shows that a better fitting to the experimental data can be obtained by multiplying the constant value by a correction factor, and the modified critical radius
The activation of the pure convective boiling condition deteriorates the heat transfer and leads to a sudden descent of the DS, which reduces the safety and efficiency of the refrigeration system. Therefore, the pure convective boiling condition can be seen as unstable condition, and the DS at critical nucleate boiling condition is known as minimal stable value (MSV). Fig. 4 shows the events that trigger the change
Page 5 of 17
of the heat transfer condition. The heat flux discriminant is used when evaporator works at the nucleate boiling condition, and the pure convective boiling condition is activated when Eq. (26) becomes true: qONB q
(26)
where q denotes a threshold that specifies the maximum value of q O N B for the onset of nucleate boiling: q
Q
(27)
Se
being Q the output cooling capacity of the evaporator, S e the heat transfer area of the evaporator ( S e 3 .0 2 m 2 , manufacturer data for the evaporator). It is found that when Eq. (26) is satisfied, the disturbance on q ONB induced by the discontinuous decrease of the heat transfer coefficient leads to the switch data oscillation. In the open-loop simulation (abrupt changes of the compressor frequency or EEV opening), the switching oscillation will eventually vanish. However, the repeating switch between two conditions will not disappear spontaneously during the simulation embedded control algorithm, which causes the pure convective boiling condition cannot be activated normally. In order to expand the application of the proposed model, the DS at the outlet of evaporator is used to judge whether the nucleate boiling condition is activated again when Eq. (28) becomes true: DS MSV
(28)
where MSV is the DS at critical nucleate boiling condition at which q O N B is equal to q and the nucleation is just activated. Therefore, it is reasonable to use the MSV as the dividing point between two different heat transfer conditions. In addition, a dead zone technology has been adopted in order to provide a smooth switching progress during the simulation. For instance, if q O N B q , the pure convective boiling condition is active. Then, the DS discriminant will not be judged in the next few seconds. Fig. 6 shows the flow diagram for dynamic simulation of the evaporator model. Considering that the refrigeration system is slow in response
speed, the simulation step size should be as large as possible without affecting the solution precision. On the other hand, the heat transfer coefficient is discontinuous with the change of heat transfer condition. Discontinuity can be troublesome in the case of large simulation step size, and it may lead to chattering or convergence problems. A dead zone method has been adopted to provide a smooth switching progress during the simulation. Therefore, the variable time-step Runge-Kutta solver is used in order to quicken simulation when ensuring the stability of the proposed model. At the beginning of simulation, the evaporator works at the nucleate boiling condition. For a given time step,#3: the coefficients are first calculated with Reviewer Response 9 the input parameters. According to the heat transfer condition at the preceding time step, the Eq. (26) or Eq. (28) would be used to judge which heat transfer condition in the evaporator is activated at the current time step. When the heat transfer condition is chosen, the thermo-dynamic properties of R22 are looked up from REFPROP7. At the end of the time step, the output parameters of differential form are obtained with matrix calculation and are calculated after integral over time. The output data of parameters are then deposited in storage database. Proceed down the next time step with the same way till the appointed time has been attained.
4. Experimental work 4.1. Experimental facility The experiment to validate the model presented above is carried out in a liquid-to-liquid chiller whose schematic diagram is shown in Fig. 5. The experimental refrigeration system is able to work with the refrigerants R134a and R22, but the refrigerant used for the tests presented here is R22. The compressor is of hermetic variable speed reciprocating type and the expansion device is a pulse-width-modulation (PWM) electric expansion device with variable opening from 0 to 100%. The evaporator and condenser are both double-pipe exchanger heat exchangers, connected to two distinct temperature and flow controlled secondary heat transfer loops. Pure water is used for the evaporator and condenser. The cooling load system regulates the chilled water temperature using immersed electrical resistances, and the nomial output cooling capacity is 35kW. The thermodynamic states of the refrigerant
Page 6 of 17
and secondary fluids at each point are calculated using measurements of pressure and temperature at the inlet and outlet of each element of the facility. These measurements are obtained using Immersion type-T thermocouples (measurement uncertainty of ±0.2 o C ) and pressure transducers (measurement uncertainty of ±0.1%). The refrigerant mass flow rate is measured by a Coriolis-effect mass flow meter with a certified accuracy within ±0.25% of the reading, and the secondary fluids mass flow rates are measured with volumetric flow meters, introducing a maximum error of ± 0.33%. Furthermore, a capacitive sensor is installed to obtain the compressor rotation speed, with a maximum error after calibration of ±1 %. A control and data acquisition system is used for monitoring the experimental variables and for setting the compressor speed and valve opening. 4.2. Identifying the output heat flux model In the refrigeration system, both of the variations of compressor speed ( f ) and EEV opening ( VO ) may effect on the cooling capacity and output heat flux. The disturbance of cooling capacity ( Δ Q ) caused by the variation of the compressor speed ( Δf ) and EEV opening ( ΔVO ) can be represented using a first-order transfer function model (Li et al., 2009). Therefore, the disturbance of output heat flux ( Δq ) is expressed as follows: Δq=G 1 (s) Δf G 2 (s) ΔV O
where
G1
and
G2
(29)
are both first-order transfer
function model: G (s)
K τs+1
G1
and
G2
can be
estimated. Considering that the change of EEV opening has little effect on the output heat flux, G 2 is approximated at zero for simplicity.
G 1 (s)=
Δq
=
0.12
Δ f 32s+1 Δq G 2 (s)= 0 ΔVO
(31) (32)
5. Results and discussion In this section the validation of the model using experimental measurements and correlation analysis are presented. To this end, two sets of experiments have been undertaken: (a) sudden change of EEV opening and compressor frequency, and (b) oscillation of the system when the DS is regulated by P–I controller. Based on the compressor frequency, the output heat flux can be obtained with Eq. (31). The evaporating pressure, condensing pressure and the DS at evaporator outlet are chosen for validation. Finally, the instability of system when the system is affected by external disturbance, is simulated. 5.1. Sudden change of EEV opening and compressor speed The comparison between simulation and experimental results of dynamic process are depicted in Fig. 8 and 9. The test is initiated when the steady-state condition is achieved for the following operating conditions in both tests: EEV opening of 40%, compressor speed of 50 Hz.
(30)
where K represents the static, τ the time constant and s the Laplace operator. In order to identify the parameters of G 1 and G2
Matlab, the parameters of
, two sets of tests have been performed. At the
beginning, fixing the compressor speed at 45Hz and EEV opening at 40% in both tests. Once in steady state, the disturbances of EEV opening (40–44) and compressor speed (45–50) are imposed respectively. The heat flux response due to the sudden reduction of compressor speed is presented in Fig. 7. Using the System Identification Toolbox in
Independent disturbances of EEV opening and compressor frequency are imposed following the sequence shown in Table 4. It is noteworthy that the disturbances of EEV opening in Fig. 8 cause a decrease of the DS at t=200s, and the DS quickly falls to zero at t=600s. The similar results can be seen in Fig. 9 at t=100s and 300s, respectively. During the tests, the change in DS is influenced by two factors. One is the total mass of refrigerant in evaporator, and the other is the heat transfer between the refrigerant and chilled water. The first step increase of the EEV opening and decrease of the compressor speed imply more refrigerant mass
Page 7 of 17
in the evaporator, which results in a longer two-phase region for evaporator, and thus the DS drops. Since the second disturbances are just a little bit stronger than the first, the differences of the the heat transfer between the refrigerant and chilled water. The comparison results show that a good matching between the simulated and experimental data is achieved not only in the nucleate boiling condition but also the pure convective boiling condition. Fig. 8.b and 9.b show the relationship between the actual heat flux and q O N B as the changes of EEV opening and compressor frequency, respectively. The output heat flux is a constant value corresponding to the compressor frequency of 50Hz in Fig. 8.b, and the variation of the output heat flux caused by the change of the compressor frequency is obtained with Eq. (31) in Fig. 9.b. Based on Eq. (23), the given amount of heat flux needed for the onset of nucleate boiling ( q O N B ) can be calculated with the changing DS. At the beginning, the DS is relatively high, and the output heat flux is clearly higher than the value needed. Therefore, the nucleate boiling condition is activated. When the DS reduces at its MSV, the output heat flux is lower than q O N B , which cause that the pure convective boiling condition is activated according to Eq. (26). As the DS is higher than MSV, the nucleate boiling condition is activated again based on Eq. (28) and deterioration of the heat transfer finally disappear. This result provides a further understanding of the stability of control loop. The difference between the output heat flux and q O N B can be seen as the stability margin, which should be positive in order to activate the nucleate boiling condition. It is seen from Fig. 8.b and 9.b that q O N B increases with the increase of EEV opening and decrease of compressor frequency. And the decrease of compressor frequency also causes the decrease of output heat flux. In practical operation, in response to system load changes, compressor speed may be continuously altered to adjust output cooling capacity to match the changes in cooling load. Due to the interaction, the EEV opening should be regulated systematically in order to drive the DS towards its reference. The analysis indicates that the decrease of compressor speed leads to the decrease of the stability margin, which makes the control of DS more difficult. begins. It is worth noting that the steady-state MSV (after the third switch) is almost the reference of
refrigerant mass in the evaporator are little. Therefore, the sudden reductions in the DS in Fig. 8.a and 9.a are due to the sudden deterioration of
5.2. Oscillation of the control loop The oscillations of system parameters when a decreasing step change on the DS is imposed, are depicted in Fig. 10. The compressor speed is 50 Hz, and the EEV opening is derived by P–I algorithm to drive the DS towards the desired references of 8 o C . For the EEV controller, a set of constant P–I values, 1.0 and 0.5, respectively, is adopted. The temperatures of the chilled and cooled water are 12oC and 30oC, respectively. The reference of the DS is step decreased to 5.6 o C at t=50s, and kept on this value. It can be seen that the oscillations of the DS near the reference can be observed after about t=100s, and the system parameters also fluctuate due to the oscillation of the EEV opening. In Fig. 10, a switch from the nucleate boiling condition to pure convective boiling condition can be seen at t=70s. Such an effect is induced by the large magnitude of the imposed change. Due to the over modulation, the DS is lower than MSV during the transient. Consequently, the pure convective boiling condition is activated, which implies the refrigerant in the evaporator cannot absorb enough heat for phase change, and then the DS drops suddenly. However, a further decline of the DS is rejected by the controller which closes the EEV to a smaller position in order to maintain the DS to its reference. As the DS is higher than MSV, the nucleate boiling condition is activated again and deterioration of the heat transfer finally disappears. However, when the DS is close to the reference of 5.6 o C at about t=150s, the pure convective boiling condition is triggered due to the insufficient of the output heat flux to support nucleate boiling, and an abruptly decrease of the DS occurs one more time. The root reason can be concluded that the new balance stability margin is too small, which causes the switch from nucleate to pure convective boiling condition is inevitable with the effect of the controller. Consequently, the periodical change of the heat transfer condition leads to the disturbance on the DS, the controller has to pursue the DS ceaselessly, and the oscillation of the control loop
the DS of 5.6 o C in the example. However, the oscillation of the system can be still emulated
Page 8 of 17
successfully because of the over modulation, as long as the new balance q O N B is close enough to the output heat flux. It is clear that actual fluctuations of the DS control loop and pressures of heat exchangers are consistent with the simulated ones. However, special attention should be paid to the large difference between the tested and simulated condensing pressure at about t=100s. The difference can be concluded that the change of EEV opening in simulation is obviously faster and larger than the actual experimental one within the time interval 50s
The system is first working under the same steady-state condition with the example in section 5.2, and the water temperature goes up 2.5 o C at t=50s. Such variation of the surrounding condition leads to the increase of q O N B and sudden drop of the DS. Albeit the disturbances on the DS have been recognized by the controller, the final stability margin is so small that the periodical switch between two heat transfer conditions cannot be eliminated.
6. Conclusion In this paper, A dynamic lumped parameter model for EEV-controlled refrigeration system has been developed. Two types of heat transfer condition in the evaporator are defined according to different heat transfer mechanisms of refrigerant. Based on the theory of bubble dynamics, a novel technique for the switch between different heat transfer conditions is presented. The difference between the output heat flux and q O N B is seem as Reviewer #1: Response 4 the stability margin. The analysis provides detailed insights to the operational stability due to the change of heat transfer mechanism of refrigerant. In addition, the dynamic behavior of oscillation of system parameters is simulated successfully by the switch between stable and unstable condition. The predicted results in comparison to the corresponding experimental data fall into a reasonable range. The discontinuous changes of the heat transfer coefficient can be used to quantify the system’s inherent characteristic which causes the instability, which is benefit to the design of controller. This study provides the basis for further developing advanced control algorithm which will ensure that an EEV-controlled system refrigeration system has better stability.
Acknowledgements
5.3. Sudden change of the temperature of cooling water
This work is supported by the Natural Science Foundation of China under Grants 61403274.
Finally, the effect of sudden change in the temperature of cooling water on the behavior of the system is simulated by the proposed model, and the transient responses of system parameters when the step changes are imposed on the cooling water temperature, are shown in Fig. 11.
References: Aprea, C., Mastrubllo, R., 2002. Experimental evaluation of electronic and thermostatic expansion valves performance using R22 and R407C. Appl Therm Eng 22, 205–218. Broersen, P.M.T., Van der Jagt M.F.G., 1980. Hunting
Page 9 of 17
of evaporator controlled by a thermostatic expansion valve. Trans ASME, J Dynam Syst Meas Control 102, 130–5. Chen, W., Chen, Z.J., Zhu, R.Q., Wu, Y.Z., 2002. Experimental investigation of a minimum stable superheat control system of an evaporator. Int. J. Refrigeration 25, 1137–1142. Chen, Y., Deng, S., Xu, X., Chan, M., 2008. A study on the operational stability of a refrigeration system having a variable speed compressor. Int. J. Refrigeration 31, 1368–1374. Gosney, W.C., 1982. Principles of Refrigeration, Cambridge. Gruhle, W.D., Isermann, H., 1985. Modeling and control of a refrigerant evaporator. J. Dyn. Syst. Meas. Control 107, 235-240. He, X.D., 1996. Dynamic Modeling and Multivariable Control of Vapor Compression Cycles in Air Conditioning Systems, PhD Thesis, Massachusetts Institute of Technology, America. Ibrahim., G.A., 2001. Effect of sudden changes in evaporator external parameters on a refrigeration system with an evaporator controlled by a thermostatic expansion valve. Int. J. Refrigeration 24, 566–576. Koury, R.N.N., Machado, L., Ismail, K.A.R., 2001. Numerical simulation of a variable speed refrigeration system, Int. J. Refrigeration 24 (2), 192–200. Kandlikar, S.G., 1990. A general correlation for saturated two-phase flow boiling heat transfer inside horizontal and vertical tubes, Journal of Heat Transfer 112, 219–228. Lars B.K., 1999. Novel electronic high reliability valve principle for control of direct expansion. In: Proceedings of 20th international congress of refrigeration, Sydney, IIR/IIF, 445. Leducq, D., Guilpart, J., Trystram, G., 2003. Low order dynamic model of a vapor compression cycle for process control design, Journal of Food Process Engineering 26, 67–91. Lenger, M.J., Jacobi, A.M., Hrojak, P.S., 1998. Superheat Stability of an Evaporator and Thermostatic Expansion Valve Air Conditioning and Refrigeration Center. College of Engineering. University of Illinois at Urbana-Champaign. Leonardo, C.S., Christian, J.L.H., Alexandre, T.N., 2010. Assessment of the controlling envelope of a model-based multivariable controller for vapor compression refrigeration systems. Appl. Therm. Eng. 30, 1538–1546. Li, X.Q., Chen, J.P., Chen, Z.J., et al., 2004. A new method for controlling refrigerant flow in automobile air conditioning. Applied Thermal Engineering 24, 1073–1085. Li H., Jeong, S.K., You, S.S., 2009. Feedforward
control of capacity and superheat for a variable speed refrigeration system, Applied Thermal Engineering 29, 1067–1074. Liang, N., Shao, S.Q., Xu, H.B., Tian, C.Q., 2010. Instability of refrigeration system–A review, Energy Conversion and Management 51, 2169–2178. Liang, N., Shao, S.Q., Tian, C.Q., Yan, Y.Y., 2011. Two-phase flow instabilities in horizontal straight tube evaporator, Applied Thermal Engineering 31, 181–187. Mithraratne, P., Wijeysundera N.E., 2001. An experimental and numerical study of the dynamic behaviour of a counter-flow evaporator. Int. J. Refrigeration 24, 554–65. Rasmussen, H., Larsen, L.F.S., 2011. Non-linear and adaptive control of a refrigeration system. IET Control Theory Appl 5, 364–378. Shah, M.M., 1979. A general correlation for heat transfer during film condensation inside pipes. Int. J. Heat Mass Transfer 22(4), 547–556. Steiner, D., Taborek, J., 1992. Flow boiling heat transfer in vertical tubes correlated by an asymptotic model. Heat Transfer Engineering 13 (2), 43–69. Tahat, M.A., Ibrahim, G.A., Probert, S.D., 2001. Performance instability of a refrigerator with its evaporator controlled by a thermostatic expansion-valve. Applied Energy 70, 233-249. Yasuda, H., Toufer, S., 1983. Simulation model of a vapour compression refrigeration system. ASHRAE Trans, 408–24. Zahn, W.R., 1964. A visual study of two-phase flow while evaporating in horizontal tubes. Trans ASME, J Heat Transfer 86, 417–29. Zhang, Q., Canova, M., 2015. Modeling and output feedback control of automotive air conditioning system. Int. J. Refrigeration 58, 207–218.
Page 10 of 17
Superheat region
Two-phase region
Two-phase region
q Superheat region
Superheat region
SV M S> D
>q
(1) evaporator
B N O
Subcooled region
Two-phase region
(2) condenser
Fig. 1. Schematic diagrams of evaporator and condenser.
u
x A -1f
x
Superheat region
Fig. 4. Switching between different heat transfer conditions in the evaporator.
Physical Parameters
Inputs
Two-phase region
Cooled water side
x
EEV
x(0)
Refrigerant side
Condenser Variable-speed compressor
Evaporator
Fig. 2. Solution procedure of the heat exchanger Chilled water side
Chilled water tank (1m 1m 1m)
Variable-speed supply pump
Fig.5. The schematic diagram of the experimental plant. Two-phase region
Superheat region
(1) evaporator in nucleate boiling condition
Two-phase region
Superheat region
(2) evaporator in pure convective boiling condition
Fig. 3. Evaporator in nucleate and pure convective boiling conditions.
Page 11 of 17
Start
Inputs Initial values
Input
Calculate coefficients
Yes
No Nucleate boiling?
Look-up
Look-up
Te , e, g , e, f , H e, g , H e, f , e,i1 (eqn.25)
Te , e, g , e, f , H e, g , H e, f , e,i1 (eqn.24)
Matrix Calculation
Matrix Calculation
dPe , dLe,1 , dTw,e
dPe , dLe,1 , dTw,e
Check for switch to convective boiling (eqn.26)
Check for switch to nucleate boiling (eqn.28)
Integral process
Integral process
Next time step Output
No
Data storage
Appointed time?
Yes End
Fig. 6. Flow diagram for the dynamic simulation program of the evaporator.
Page 12 of 17
Heat flux(kW m-2)
10.2 10 9.8 9.6 9.4 9.2 9 0
100
200
300
400
500
600
700
800
t(s)
Fig. 7. The heat flux response due to the change of compressor speed
Page 13 of 17
a
b
8
Heat flux(kW·m-2)
Degree of superheat(oC)
10
6 4 2 0 0
200
400
600
800
Output heat flux
10 9.5 qONB
9 8.5
1000
Data Model
10.5
0
200
400
t(s)
d
505
495
P c(MPa)
P e(kPa)
c
485
475
0
200
400
600
800
1000
600
800
1000
t(s)
600
800
1.41
1.39
1.37
1.35
1000
0
200
400
t(s)
t(s)
Fig. 8. Validation for the transient response to the change of the EEV opening: a. degree of superheat, b. relationship between output heat flux and qONB, c. evaporating pressure, and d. condensing pressure.
a
b
10
Output heat flux
Heat flux(kW·m-2)
Degree of superheat(oC)
8 6 4 2 0
0
Data Model
10.5
100
200
300
400
10 9.5 9 qONB
8.5
500
0
100
200
t(s)
c
d
520
P c(MPa)
P e(kPa)
400
500
300
400
500
1.41 1.4
510 500 490 480
300
t(s)
1.39 1.38 1.37
0
100
200
300
t(s)
400
500
1.36
0
100
200
t(s)
Fig. 9. Validation for the transient response to the change of the compressor speed: a. degree of superheat, b. relationship between output heat flux and qONB, c. evaporating pressure, and d. condensing pressure.
Page 14 of 17
b
55 50
40
495 485
35 30 0
Model Data
515 505
45
Pe (kPa)
EEV opening(%)
a
100
200
300
400
475
500
0
100
200
1.6
c Pc(MPa)
1.5 1.4 1.3 1.2
0
100
200
300
t(s)
300
400
500
300
400
500
t(s)
400
500
d
9
Degree of superheat(oC)
t(s)
7
5
3 0
100
200
t(s)
Fig. 10. System oscillation analysis: a. EEV-opening, b. evaporating pressure, c. condensing pressure, and d. degree of superheat.
Page 15 of 17
EEV opening(%)
45
40
35 0
100
200
300
400
500
600
700
400
500
600
700
t(s) Degree of superheat(oC)
7 6.5 6 t(s)
5.5 5 0
100
200
300 t(s)
Fig. 11. System oscillation when the system is disturbed by the abrupt change of external parameter
Table 1. Indication expressions of coefficient matrix for the evaporator. Symbol
Indication expression
a e ,11
é d r e, f d r e,g ù ú Ae l e ,1 ê (1 - g e ) +g e dPe dPe ú ëê û
a e ,12
Ae (1 - g e )( r
a e , 21
é ù d r e , f he , f d r e , g he , g Ae l e ,1 ê (1 - g e ) +g e - 1ú êë ú dPe dPe û
a e , 22
Ae (1 - g e )( r
a e ,33
C p , cu r
- r
e, f
e, f
e,g
)
he , f - r
e,g
he , g )
A cu e , cu
Table 2. Calculation expressions of coefficients in heat exchangers. Symbol
Expression
e , i1
S. G. Kandlikar equation (Kandlikar, 1990): i C1Co 2 conv C 3 Bo 4 F fl conv C
c , i1 , c , i 3
C
Dittus-Boelter equation: Nu 0.023 Re i
e ,i 2 , c ,i 2
0.8
0.4
Pr ,
Nu di
Traviss, Baron and Rohsenow equation: Nu
Re
0.9
Pr F1
, i
F2
e,o , c ,o
Nu di
Shah equation(Shah, 1979): o Bw
0.8
(
Di
)
R
0.18
0.2
Do
,
De
B 1496 22T w t
Table 3. Indication expressions of coefficient matrix for the condenser. Symbol
Indication expression
a c ,11
Ac r
a c ,13
d r c ,1 1 é ïì Ac l c ,1 í ( H c ,1 - H c , g ) + êr dPc 2 êë ïî
c ,1
( H c ,1 - H c , g )
a c , 21
Ac r
c, f
a c , 22
Ac r
g c ( H c,g - H c, f ) c, f
a c , 23
ìï Ac l c ,2 í (1 - g c ) r ïî
a c ,31
-
a c ,32
-
a c ,33
Ac L c ,3 (
1 2 1 2
c ,1
+ ( H c ,1 - H c , g )
ù dH c , g ïü c ,1 ú - 1ý dH c úû dPc e ïþ
dr
( H c,g - H c, f )
dH c , f c, f
dPc
é d r c,g + g c ê( H c , f - H c , g ) +r êë dPc
Ac r
c, f
( H c , f - H c ,o )
Ac r
c, f
( H c , f - H c ,o )
1 2
r
dH c , f c, f
c,g
dH c , g ù údPc úû
üï 1ý ïþ
- 1)
dPc
Page 16 of 17
1
a c ,34
2
Ac L c ,3 r
c, f
a c , 41
Ac ( r
a c , 42
Ac g c ( r
a c , 43
ìï Ac í L c ,2 îï
a c ,51
Ac , cu C p , cu r
cu
(T c , w 1 - T c , w 2 )
a c ,55
Ac , cu C p , cu r
cu
L c ,1
a c , 66
Ac , cu C p , cu r
cu
a c , 71
Ac , cu C p , cu r
cu
(T c , w 3 - Tc , w 2 )
a c , 72
Ac , cu C p , cu r
cu
(T c , w 3 - Tc , w 2 )
a c , 77
Ac , cu C p , cu r
cu
L c ,3
c ,1
- r c,g
c, f
- r
) c, f
)
é d r c, f d r c,g ê (1 - g c ) +g c êë dPc dPc
ù æ d r c ,1 1 d r c ,1 dH c , g ú + L c ,1 ç + ç úû 2 dH c dPc è dPc
ö üï ÷ý ÷ ø þï
Table 4. Changes imposed to the valve opening and compressor speed Time(s)
Av[%] fcom[Hz]
t<100
100
200
300
400
40 50
43 45
40 50
44 44
40 50
Page 17 of 17