Next Article in Journal
Numerical Simulation Study on the Damage Mechanism of the Combined Perforating Well Testing Tubing in Ultra-Deep Wells
Previous Article in Journal
Study on Numerical Simulation of Formation Deformation Laws Induced by Offshore Shallow Gas Blowout
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Modeling Method for Overheated Zone and Two-Phase Zone of Dry Shell-and-Tube Evaporator in Ship Air Conditioning

1
Marine Engineering College, Dalian Maritime University, Dalian 116026, China
2
Maritime College, Tianjin University of Technology, Tianjin 300384, China
*
Author to whom correspondence should be addressed.
Processes 2024, 12(2), 379; https://doi.org/10.3390/pr12020379
Submission received: 15 December 2023 / Revised: 8 February 2024 / Accepted: 9 February 2024 / Published: 13 February 2024

Abstract

:
This paper researches the heat transfer equation and thermal balance equation of a shell-and-tube evaporator; constructs an accurate mathematical model for the evaporator; and derives equations including detailed and accurate calculation methods for all heat transfer coefficients, such as the refrigerant side heat transfer coefficient, water side heat transfer coefficient, refrigerant kinematic viscosity, density, and specific enthalpy. Adopting this approach involves fitting the relationships between the density, thermal conductivity, kinematic viscosity, and enthalpy of R134a refrigerants in saturated vapor and liquid states. The relationships between superheated gas enthalpy, density, and temperature were also assessed, and heat transfer coefficients were obtained through calculation methods and microelement heat transfer relationships in both the single-phase and two-phase zones, matching empirical formulas concerning the relationship between superheated enthalpy and temperature. Notably, the research utilizes the Simulink approach without relying on M files and S functions to establish the evaporator’s two-phase and superheated zones, as well as an overall simulation model which provides intuitive internal coupling relationships and the coefficient calculation process in the formulas and uses the function “Algebraic Constraint” instead of “memory” or “1/z” to solve algebraic loops, thereby avoiding computation deviations introduced by delays and iterations. Finally, simulation calculations were conducted, and an experimental platform was designed and built for experimental verification which can validate the derived mathematical models. The simulation results, including the evaporator pressure, and chilled water outlet temperature with variation in chilled water mass flow rate, closely matched the experimental outcomes. The simulation model is concise and intuitive. Modifying parameters such as the thermal conductivity of the model material is straightforward, thereby alleviating the workload for researchers. It also facilitates an understanding of model principles for beginners. Moreover, the database generated from the model allows for fault analysis, diagnosis, and decision evaluation.

1. Introduction

A ship’s air conditioning system is a crucial piece of auxiliary equipment, providing a comfortable rest environment and serving as a key guarantee for the safety and efficient production of the crew. The evaporator, one of the main four components of the refrigeration unit, significantly influences the efficiency of the entire refrigeration system due to its heat exchange characteristics. Although there is a wealth of literature and various methods for researching ship air conditioning systems, the current mature approach to studying mathematical models of ship air conditioning systems typically involves initially independently constructing the heat transfer and exchange equations for the four major components of the air conditioning system. Subsequently, these equations are simultaneously solved to obtain the overall performance of the system. However, this method often lacks clarity in conveying the relationships between parameters such as the pressure, temperature, and enthalpy of refrigerants and cooling mediums as they transfer through the four components. In contrast, the use of Simulink to construct mathematical models offers an intuitive, concise, and easily understandable approach. For researchers, this method is very user-friendly and popular, significantly reducing the workload of research tasks and making it easy for beginners to grasp the principles of the model.
In 2021, Pritam Dey and colleagues evaluated the overall performance of a shell-and-tube evaporator using numerical methods. They analyzed the flow boiling heat transfer characteristics of nanorefrigerants inside the tubes using correlations for nanofluids. The numerical solution of the boiling and pressure drop model for the nanorefrigerants in the evaporator was determined using EES software [1]. In 2022, Feihu Chen et al. investigated the heat transfer characteristics of the evaporator terminal, cold distribution unit, and refrigerant flow distribution in a water-cooled multiple jet heat pipe system for data centers [2]. In 2022, Moharana et al. analysed the effects of different parameters affecting heat transfer in tube bundles on the shell side of two-phase shell-and-tube heat exchangers based on the available literature. The studies included were subdivided into effect of performance parameters for flow boiling and for pool boiling [3]. In 2022, Santa et al. proposed a suitable mathematical model for the physical process in the tubular evaporator of the heat pump. Theoretical results were obtained by dividing the evaporator into control volumes and solving the corresponding system of equations of the proposed mathematical model using the Runge–Kutta and Adams–Moulton predictor–corrector methods [4]. While these scholars have established models for refrigeration units, none of them have developed a mathematical model that distinguishes the superheated zone and two-phase zone. Additionally, they have not constructed Simulink models and simulations for these specific components.
In 2021, Xu Guowen proposed a parallel collaborative simulation method to accurately predict the dynamic changes in the thermal environment of air-conditioned rooms. This method aimed to enhance the efficiency of testing in the air conditioner thermal comfort laboratory, achieving efficient coupling and collaborative simulation between the air conditioning system model and the room thermal environment model [5]. In 2022, Zhiyi Wang et al. established a Simulink system simulation and Simulink-M file optimization model to study the performance of the unit in variable outdoor air conditions by combining experiments, and the optimal performance of coefficient and the corresponding optimal bypass fresh air flow rate of the unit in different outdoor temperature was obtained [6]. In 2023, H.H. El-Ghetany and others used Simulink for the numerical analysis of a refrigeration system using silicone gel and water. They conducted a case study modeling various cooling load capacities for adsorption systems [7]. While these studies introduce the establishment of refrigeration system mathematical models based on Simulink, a thorough examination of the literature did not reveal detailed information on the calculation methods for all coefficients in the formulas. Moreover, there is no indication of precise coefficient solutions. Additionally, a complete construction of a refrigeration system mathematical model entirely based on the Simulink approach, along with a detailed explanation of the coefficient calculation process in the formulas, was not found in these sources.
In the exploration of control decisions for ship air conditioning systems, the authors encountered the challenge of a lack of intuitive and easily modifiable platforms for studying the different structures and performance parameters of the air conditioning system. This limitation made it difficult to obtain diverse data and observe various phenomena. Consequently, the authors aim to construct a Simulink simulation platform for the refrigeration system. Considering that the heat exchange and transfer models for the evaporator and condenser are the most complex components of the refrigeration system, particularly due to the phase change in the refrigerant, this paper focuses on studying the evaporator as a key component.
Before establishing the evaporator model, several assumptions need to be made:
(1)
The refrigerant flows in a one-dimensional manner within the tube, with only axial movement considered and no consideration of radial flow;
(2)
The cross-sectional area inside and outside the heat exchanger tube remains constant, neglecting minor variations in structure at the interface;
(3)
The refrigerant outside the tube flows in a one-dimensional manner, engaging in pure counter-current heat exchange with the cold water;
(4)
The thermal resistance of the tube wall is ignored;
(5)
Pressure drop is neglected;
(6)
It is assumed that the gas–liquid mixture in the two-phase region is uniformly mixed, and the mass flow rate remains constant.
Based on the assumptions mentioned above, this paper initially analyzed the heat exchange processes in the superheated and two-phase zones of a shell-and-tube evaporator. Subsequently, thermal balance equations and heat transfer equations for the refrigerant and chilled water on both sides of the shell were established. Then, by employing calculation methods for surface heat transfer coefficients on the chilled water side and the refrigerant side, as well as the differential heat transfer calculation method, formulas for calculating the overall heat transfer coefficient in the thermal balance and heat transfer equations were provided. The specific procedures are as follows:
Initially, by establishing heat exchange and transfer formulas for the two zones, along with the calculation methods for all coefficients, such as the refrigerant side heat transfer coefficient, water side heat transfer coefficient, refrigerant kinematic viscosity, density, and specific enthalpy within these formulas, a precise mathematical model for the evaporator is developed. Subsequently, utilizing the Simulink approach without relying on any script functions, the paper constructs a detailed model for both the two-phase and superheated zones of the evaporator, followed by comprehensive simulation calculations [8]. Finally, an experimental setup for an air conditioning system is designed and implemented for validation. The calculation of formula coefficients involves matching empirical formulas or fitting relationships for the density, thermal conductivity, kinematic viscosity, and enthalpy of R134a refrigerants in saturated vapor and liquid states. This fitting process is combined with the fitting formulas for superheated gas from refs. [9,10]. Through this approach, all coefficients related to heat transfer are derived for the equations in both the single-phase and two-phase zones. The schematic diagram of the paper’s structure is shown in Figure 1:
When solving linear equations in the constructed simulation model, the function “Algebraic Constraint” is employed instead of “memory” or “1/z” functions. This choice mitigates computational inaccuracies arising from delay and iterative calculations, ensuring precise results.
The simulation model can be utilized to investigate parameters such as the structural performance of the evaporator and refrigerant flow control. It facilitates the selection and optimization of the evaporator and enables precise control of the system. Moreover, the process and approach used to construct the evaporator model can also be applied to develop a Simulink model for the condenser.
The author’s next plan involves establishing a Simulink model for the three zones of the condenser and subsequently integrating the compressor, expansion valve, and environmental load models to create a comprehensive Simulink simulation platform for the entire air conditioning system. This platform can be used to continuously monitor changes in variables such as evaporator pressure, evaporator temperature, condenser pressure, and condenser temperature, as well as temperature parameters within the refrigerated compartments. Ultimately, this platform will support the establishment of a comprehensive fault diagnosis model, control predictions, and decision libraries for the entire refrigeration unit [11,12].

2. Evaporator Heat Transfer Process

The characteristics of R134a are similar to R12, which is colorless, odorless, non-toxic, non-flammable, and non-explosive. It has a higher latent heat of vaporization than R12, is incompatible with mineral lubricating oil, does not contain chlorine atoms, and does not deplete the ozone layer. It possesses excellent safety performance and is considered an excellent long-term alternative refrigerant. Therefore, in this paper and the experimental apparatus, R134a is chosen as the refrigerant.
The presence of a subcooled liquid zone is generally not considered unless it is specific to certain refrigeration systems. Therefore, based on the different states of the refrigerant, the evaporator can be divided into two zones: the superheated zone and the two-phase zone.
The enthalpy–pressure relationship of the refrigerant during the circulation in the refrigeration unit is depicted in Figure 2. Point 1 represents the superheated gaseous state of refrigerant at the compressor outlet, points 2 and 3 respectively represent the saturated gaseous and saturated liquid states of refrigerant in the condenser, point 4 represents the subcooled liquid state of refrigerant at the condenser outlet, point 5 represents the two-phase state of refrigerant at the evaporator inlet, point 6 represents the saturated gaseous state of refrigerant in the evaporator, and point 7 represents the superheated gaseous state of refrigerant at the evaporator outlet. In Figure 3, points 5, 6, and 7 have the same meanings as described above. This paper focuses on modeling and simulating the evaporator, which corresponds to points 5–7 in the diagram. The process is divided into the two-phase zone (5–6) and the superheated zone (6–7). The lengths of the two-phase and superheated zones, denoted as l e , 1 and l e , 2 , are continuously variable during system operation and are represented by z 41 and z 42 , respectively.
The variations in refrigerant enthalpy and the process of chilled water temperature change are illustrated in Figure 3. The refrigerant enters at point 5 and exits at point 7, while the chilled water flows in the opposite direction, entering at point 7 and exiting at point 5. At point 5, the refrigerant is in a low-temperature liquid–vapor mixture state, while point 6 represents saturated gas and point 7 corresponds to superheated gas. The specific enthalpies of the refrigerants at locations 5, 6, and 7 are h r , e , 1 , h r , e , 2 , and h r , e , 3 , and the temperatures are z 11 , z 12 , and z 13 , respectively. The water temperature is represented as t w , e , 3 , t w , e , 2 , and t w , e , 1 , denoted as z 23 , z 22 , and z 21 , respectively. The tube wall temperatures at points 5 and 6 are denoted as t p , e , 1 and t p , e , 2 , with the average temperature, t p , e , represented by z 31 . t e and t c represent the evaporator and condenser temperatures, respectively.
If the heat flux of the refrigerant transferred to the chilled water is referred to as the heat transfer heat flux of the heat exchanger and if the heat dissipation heat flux of the heat exchanger to the external environment is neglected, then the heat transfer heat flux through the heat exchanger, the heat release heat flux of the refrigerant, and the heat absorption heat flux of the chilled water are equal. The heat transferred during the heat transfer process, i.e., the heat transfer heat flux, ϕ , of the heat exchanger, can be expressed using the average heat transfer coefficient, K , the heat transfer area, A , and the average temperature difference, Δ t m , i.e., ϕ = K A Δ t m . “The heat transfer process” refers to the process in which heat is transferred from one side of the fluid at the wall to the other side through the wall. The key factors are the heat transfer coefficient and the average temperature difference, which will be discussed in detail in the equations for the superheated zone and two-phase zone that follow.

3. Evaporator Heat Transfer Equations

3.1. Superheated Zone

Due to the equality of heat absorption by the refrigerant on the inner side of the shell and the heat release by the chilled water on the outer side, the following thermal balance and heat transfer equations are obtained [13]:
m r , e h r , e , 3 h r , e , 2 = m w , e c p , w , e z 21 z 22
m r , e h r , e , 3 h r , e , 2 = K e , 1 A e , 1 z 21 + z 22 z 12 t r , e , 3 2
m r , e is the refrigerant mass flow rate (kg/s); m w , e is the chilled water mass flow rate (kg/s); h r , e , 2 and h r , e , 3 are the inlet and outlet enthalpy values of the refrigerant (kJ/kg); K e , 1 is the overall heat transfer coefficient of the condenser (kW/(m2·K)); A e , 1 is the total surface area of the tube (m2); and c p , w , e is the specific heat of chilled water at constant pressure (kJ/(kg·°C)).
The coefficients in the equations are calculated as follows:
The calculation of h r , e , 2 is based on the relationship between saturated vapor enthalpy and temperature for R134a:
h r , e , 2 = 1.145 × 10 5 × z 12 3 0.001375 × z 12 2 + 0.5846 × z 12 + 398.8
According to reference [9], the relationship between the refrigerant’s superheated temperature, saturated enthalpy, and superheated enthalpy allows for the calculation of the superheated enthalpy. The relationship can be expressed as follows:
B h r , e , 3 3 [ b 3 ( z 13 + 273 ) 3 + b 2 ( z 13 + 273 ) 2 + b 1 ( z 13 + 273 ) + 1 ] h r , e , 3 2 [ b 6 ( z 13 + 273 ) 3 + b 5 ( z 13 + 273 ) 2 + b 4 ( z 13 + 273 ) + 1 ] h r , e , 3 [ b 9 ( z 13 + 273 ) 3 + b 8 ( z 13 + 273 ) 2 + b 7 ( z 13 + 273 ) + 1 ] = 0
where, b 1 = 0.3501551 , b 2 = 1.967904 × 10 3 , b 3 = 3.227027 × 10 6 , b 4 = 7.485942 , b 5 = 1.734667 × 10 1 , b 6 = 3.091114 × 10 4 , b 7 = 625.2481 , b 8 = 6.143028 , b 9 = 9.275674 × 10 3 ,
B = [ 1 + b 1 ( z 12 + 273 ) + b 2 ( z 12 + 273 ) 2 + b 3 ( z 12 + 273 ) 3 ] / h r , e , 2 + [ 1 + b 4 ( z 12 + 273 ) + b 5 ( z 12 + 273 ) 2 + b 6 ( z 12 + 273 ) 3 ] / h r , e , 2 2 + [ 1 + b 7 ( z 12 + 273 ) + b 8 ( z 12 + 273 ) 2 + b 9 ( z 12 + 273 ) 3 ] / h r , e , 2 3
The formula for the heat exchanger’s heat transfer area and heat transfer coefficient in Equation (2) is as follows:
K e , 1 = [ 1 a r , e , gs + d o , e d i , e 2 λ p , e × d o , e d i , e + 1 a w , e , gs × d o , e d i , e ] 1 × 0.001
A e , 1 = n r , e π d o , e + d i , e 2 z 41
a w , e , gs = ( 1395.6 + 23.26 × ( z 21 + z 22 ) 2 ) × u w , e 0.8 d i , e 0.2
a r , e , gs is the heat transfer coefficient on the refrigerant side of the superheated zone (kW/(m2·K)); a w , e , gs is the heat transfer coefficient on the water side of the superheated zone (kW/(m2·K)); u w , e is the flow velocity of chilled water inside the tube (m/s); d o , e is the outer tube diameter (m); d i , e is the inner tube diameter (m); λ p , e is the thermal conductivity of the tube wall (kW/(m·K)); and n r , e is the number of tubes.
The determination of the heat transfer coefficient in the single-phase zone can be obtained as follows:
α r , e , gs = λ r , e , gs d o , e × 0.023 × ( u r , e , gs × d o , e ν r , e , gs ) 0.8 × ( ν r , e , gs × ρ r , e , gs × c p - r , e , gs λ r , e , gs ) 0.3
where, λ r , e , gs = 0.01 , u r , e , gs = 4 × m r , e ρ r , e , gs × n r , e × π × d i , e 2
u w , e = m w , e ρ m , e × ( S e n r , e × π × d o , e 2 4 )
ν r , e , gs = 5.268 × 10 12 × ( z 13 + z 11 2 ) 3 + 6.457 × 10 10 × ( z 13 + z 11 2 ) 2 2.85 × 10 8 × ( z 13 + z 11 2 ) + 7.157 × 10 7
c p - r , e , gs = 9.313 × 10 7 × ( z 13 + z 11 2 ) 3 + 2.496 × 10 5 × ( z 13 + z 11 2 ) 2 + 0.003617 × ( z 13 + z 11 2 ) + 0.8995
λ r , e , gs is the heat transfer coefficient on the refrigerant side of the superheated zone, (kW/(m2·K)); u r , e , gs is the refrigerant flow velocity (m/s); ν r , e , gs is the kinematic viscosity of refrigerant (m2/s); ρ r , e , gs is the density of refrigerant (kg/m3); c p - r , e , gs is the specific heat capacity of refrigerant at constant pressure (kJ/(kg·°C)); and S e is the overall cross-sectional area of the evaporator (m2).
Based on the fitting formula for superheated density from reference [9], the superheated density is determined. The saturated gas density in the equation is obtained using the following fitted relationship:
ρ r , e , g = 9.972 × 10 5 × z 11 3 + 0.007164 × z 11 2 + 0.4437 × z 11 + 14.52

3.2. Two-Phase Zone

The symbols in the two-phase zone have the same meanings as in the superheated zone, with “g” representing saturated vapor, “gs” representing the superheated zone, “l” representing saturated liquid, and “lg” representing the two-phase zone. The derivation of the equations in the two-phase zone follows a similar process to the superheated zone:
(1)
Thermal balance equations and heat transfer equations for refrigerant and chilled water on both sides of the shell tube:
m r , e h r , e , 2 h r , e , 1 = m w , e c p , w , e z 22 z 23
m r , e h r , e , 2 h r , e , 1 = K e , 2 A e , 2 z 22 + z 23 2 z 11 2
h r , e , 1 = 6.752 × 10 6 × t c 3 + 0.001531 × t c 2 + 1.341 × t c + 199.4
(2)
Calculation method for the overall heat transfer coefficient in the formula:
K e , 2 = [ 1 a r , e , lg + d o , e d i , e 2 λ p , e × d o , e d i , e + 1 α w , e , lg × d o , e d i , e ] 1 × 0.001
A e , 2 = d o , e + d i , e 2 n r , e x 42
(3)
Calculation method for the surface heat transfer coefficient on the chilled water side:
a w , e , lg = ( 1395.6 + 23.26 × ( z 22 + z 23 ) 2 ) × u w , e 0.8 d i , e 0.2
(4)
Calculation method for the surface heat transfer coefficient on the refrigerant side [14,15]:
a r , e , lg = 0.725 n m 9.81 ρ r , e , lg λ r , e , lg 3 v r , e , lg 1 4 r e , lg 1 4 ( z 31 z 11 ) 1 4 d o , e 1 4
where, ρ r , e , lg = α ρ r , e , g + ( 1 a ) ρ r , e , l , λ r , e , lg = a λ r , e , g + ( 1 a ) λ r , e , l , ν r , e , lg = μ r , e , lg / ρ r , e , lg , μ r , e , lg = χ μ r , e , g + ( 1 χ ) μ r , e , l , r r , e , lg = h r , e , g h r , e , 1 .
ρ r , e , lg is the average density of the refrigerant in the two-phase zone (kg/m3); λ r , e , lg is the average thermal conductivity of the refrigerant in the two-phase zone (kW/(m·K)); ν r , e , lg is the average kinematic viscosity of the refrigerant in the two-phase zone (m2/s); r r , e , lg is the latent heat of vaporization of the refrigerant (kJ/kg); μ r , e , lg is the average dynamic viscosity of the refrigerant in the two-phase zone (Pa·s); α is the void fraction; and χ is the mass dryness fraction.
The relationships between the density, thermal conductivity, kinematic viscosity, and specific enthalpy of saturated vapor and temperature are fitted to obtain the relationships as follows:
λ r , e , g = 0.01
μ r , e , g = c r , e , g e x p m r , e , g θ T
where, c r , e , g = 14.27 × 10 6 , m r , e , g = 0.902 ,
θ T = l n ( z 11 + 273 ) + 0.5 [ 1 1 / ( z 11 + 273 ) ] 2 × [ 1 1 / ( z 11 + 273 ) l n ( z 11 + 273 ) ] × { 1 0.1 [ 1 1 / ( z 11 + 273 ) ] 4 } .
Subsequently, based on the numerical relationships between the density, thermal conductivity, kinematic viscosity, and specific enthalpy of saturated liquid refrigerant with temperature, the corresponding empirical correlation is derived:
ρ r , e , l = 9.914 × 10 5 × z 11 3 0.007573 × z 11 2 3.24 × z 11 + 1294.0
λ r , e , l = 0.07
μ r , e , l = c r , e , l e x p { b θ ( T ) + m r , e , l θ ( T ) 1 / 3 }
where, c r , e , l = 41.06 × 10 6 , m r , e , l = 1.962 , b = 1.8 .
(5)
The calculation method for infinitesimal heat transfer:
Due to the introduction of a new unknown parameter, namely, the tube wall temperature, when calculating the surface heat transfer coefficient on the refrigerant side in the two-phase zone, additional infinitesimal heat transfer equations on both sides of the tube wall need to be included accordingly. The concept of infinitesimal heat transfer involves dividing the heat exchanger into several computational elements along the flow direction of the refrigerant and working fluid. Each computational element can be considered as a finite-length small heat exchanger. For each computational element, the inlet states of the refrigerant and working fluid are known. By solving the energy conservation equations established on each element, the outlet states of the refrigerant and working fluid for each element can be obtained:
a r , e , lg ( z 31 z 11 ) = [ ( z 22 + z 23 ) 2 z 31 ] / [ ( 1 a w , e , lg + d o , e d i , e 2 λ p , e ) d o , e d i , e ]
(6)
Evaporator tube length:
l = z 41 + z 42

4. Building the Simulink Model

According to the above relationships, construct the following Simulink simulation diagrams. Figure 4 provides an overall simulation diagram of the evaporator for both the superheated and two-phase zones based on the input–output relationships. The calculation models for the chilled water, refrigerant, and evaporator performance parameters are shown in Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11 and Figure 12.
Figure 5 presents the Simulink solution model for Equation (2), established for the superheated zone. In this model, the input parameters are the temperature of the chilled water entering the evaporator, z 21 , mass flow rate, m w , e , refrigerant mass flow rate, m r , e , and evaporation temperature, z 12 . The relationship expressions for the two unknowns, namely, the outlet temperature of chilled water, z 22 , and refrigerant outlet temperature, z 13 , in the superheated zone, are solved through the function f ( z ) = 0 . Then, in Figure 6, z 22 is taken as a known parameter input, allowing the solution of these two unknowns. Subsequently, in the two-phase zone, z 22 is used as a known input again for the calculation of the two-phase zone equations.
In Figure 7, the coefficient B is determined based on the evaporation temperature and refrigerant saturated gas specific enthalpy. Subsequently, this value of B is incorporated into the empirical formula relating superheated specific enthalpy and temperature, allowing the calculation of the superheated specific enthalpy, h r , e , 3 .
In Figure 8, the heat transfer coefficient is calculated based on the evaporator heat transfer formula. The water-side surface heat transfer coefficient and refrigerant-side surface heat transfer coefficient involved in the heat transfer formula are calculated in Figure 9 and Figure 10, respectively. The refrigerant-side surface heat transfer coefficient is determined based on the flow velocity, kinematic viscosity, specific heat at constant pressure, and empirical formula fitted to the refrigerant’s superheated density in the superheated zone. This is established in Figure 11 and Figure 12, and its value is computed. Then, this value is input into Figure 10 to obtain the refrigerant-side surface heat transfer coefficient.
The Simulink diagram for the two-phase zone allows for iterative solutions to determine the evaporator temperature and chilled water outlet temperature in the two-phase zone. The evaporator temperature is taken as a known condition and fed into the superheated zone as a known parameter input. The chilled water outlet temperature in the two-phase zone is then determined, serving as the evaporator outlet temperature. The construction of the Simulink model diagram for the two-phase zone follows a similar approach to the superheated zone, and detailed description is omitted here.

5. Simulation Results

This paper conducted Simulink simulations and established an experimental system to validate the simulation results. The experimental setup consists of components such as a compressor, dry shell-and-tube evaporator, shell-and-tube condenser, electronic expansion valve, water pipes, and water valves. The experimental platform connections are arranged as shown in Figure 13: the sequence of pipeline connections involves the compression of gas refrigerant from the compressor to the evaporator, where it becomes a high-temperature, high-pressure superheated gas. The gas then goes to the condenser, where it condenses into a high-temperature subcooled liquid. After passing through the expansion valve for throttling and pressure reduction, it flows into the evaporator. There, it absorbs heat, turns into a low-temperature gas, and finally returns to the compressor inlet. Cooling water for the condenser is provided by a centrifugal pump, and the refrigerant in the evaporator releases heat, then enters the air conditioner to absorb heat. The air, cooled and dehumidified by the air conditioner, is sent into the cabin by a fan. Temperature sensors are installed at the inlet and outlet of the refrigerant in the evaporator pipeline and the inlet and outlet of the cooling water in the condenser pipeline to measure the temperatures. Pressure gauges are installed at the inlet and outlet of the compressor to measure suction and discharge pressures.
In the experimental unit room, with an outdoor temperature of 34 °C, the external heat infiltration is 2.5 kW, lighting and equipment contribute 3.5 kW, and the human load is 12.5 kW. The compressor has a rated power of 9 kW, and the initial assumed value for the refrigerant mass flow rate, m r , e , is 0.5 kg/s, set as a global variable. The evaporator is designed for a heat exchange capacity of 23 kW, with an inner diameter, d i , e , of 12 mm and an outer diameter, d o , e , of 16 mm, while the chilled water mass flow rate m w , e is set to 0.5 kg/s. For the condenser, the designed heat exchange capacity is 35 kW, with an inner diameter, d i , c , of 14 mm and an outer diameter, d o , c , of 18 mm, while the cooling water mass flow rate, m w , c , is set to 2.67 kg/s.
After the experimental system is operational, the air conditioning system initiates the cooling process. The indoor temperature gradually decreases, and the chilled water temperature also decreases. The evaporator pressure and temperature decrease accordingly. During the process of evaporator pressure decreasing from 2.43 bar to 0.97 bar, the evaporator temperature decreases from −5 °C to −28 °C. The simulated and measured values of the chilled water inlet and outlet temperatures during this process are shown in Figure 14 and Figure 15.
Figure 14 shows the temperature of the chilled water entering the superheated zone of the evaporator, which is input parameter for simulation, and Figure 15 depicts the outlet temperature of the chilled water in the two-phase zone of the evaporator. From the figures, it can be observed that, initially, the temperature of the chilled water entering the superheated zone of the evaporator is 32.4 °C. As the refrigerant absorbs heat and the chilled water releases heat in the evaporator, the outlet temperature decreases to 20.6 °C. Subsequently, the chilled water enters the air conditioner, absorbing heat as the indoor temperature drops. With the gradual decrease in indoor temperature, the inlet and outlet temperatures of the chilled water also decrease, reducing the heat load. This, in turn, leads to a gradual decrease in evaporator pressure, reaching 19.1 °C and 10.1 °C in the later stages.
Additionally, from the graphs, it can be observed that the initial temperature difference between the inlet and outlet of the chilled water is 11.8 °C, which later decreases to 9 °C. This phenomenon is attributed to the fact that, during the initial stages of system operation, the indoor temperature is high, resulting in a significant heat load and higher heat absorption by the refrigerant in the evaporator. As the indoor temperature decreases over time, the heat load diminishes, and the refrigerant absorbs less heat in the evaporator, leading to a reduction in the temperature difference between the inlet and outlet of the chilled water.
Figure 16 is a three-dimensional coordinate graph with the chilled water inlet temperature on the x-axis, evaporator pressure on the y-axis, and chilled water outlet temperature on the z-axis. The chilled water inlet temperature and evaporator pressure are input parameters, while the chilled water outlet temperature is the unknown parameter being solved. It is evident from the graph that the simulated values closely match the measured values.
In the mathematical model of the evaporator, by considering the chilled water inlet and outlet temperatures as known conditions and treating the evaporator pressure as an unknown variable, mathematical equations can be established following the modeling principles outlined in the text, particularly those for the superheated and two-phase zones. Consequently, the evaporator pressure can be solved. In this experimental system, the temperatures measured by temperature sensors in the chilled water inlet and outlet pipelines of the evaporator are used as known conditions. These temperatures are input into the heat balance and heat transfer equations for the evaporator. This allows the inverse solution for the evaporator pressure. The simulated and measured values of the evaporator pressure during the variation in inlet and outlet temperatures are presented in Figure 17, indicating a close match between the simulated and measured values.
To further validate the reliability of the simulation model, the mass flow rate of the chilled water was increased from 0.5 kg/s to 0.6 kg/s. Subsequent simulation and experimental tests were conducted, and the actual and theoretical values of the chilled water outlet temperature are depicted in Figure 18. In the figure, the green “o” and red “*” lines represent the simulation and measured values of the chilled water outlet temperature when the mass flow rate is 0.5 kg/s. The blue “+” and pink “∆” lines, on the other hand, depict the scenario when the mass flow rate of chilled water is 0.6 kg/s.
It can be observed that after increasing the mass flow rate of the chilled water, the simulation and measured values of the chilled water outlet temperature closely align, and the trend of change remains essentially consistent with the period before the alteration in mass flow rate. Additionally, from the graph, it is evident that compared to the period before the change in mass flow rate, the chilled water outlet temperature initially increases and then steadily decreases. This is because, under constant heat absorption conditions, the change in temperature of the chilled water is inversely proportional to its mass flow rate. In the initial phase of flow rate increase, the temperature rise is mitigated; as the refrigerant’s cooling capacity increases, more heat is absorbed, leading to a faster temperature decrease in the later stages.
During the refrigeration process, the simulated values of evaporator pressure, chilled water inlet temperature, and chilled water outlet temperature obtained from the simulation model closely match the measured data. This alignment suggests that the heat exchange and heat transfer model constructed in this study is accurate, and the formula coefficient calculation method is precise. Moreover, the model can also compute the unknown lengths l e , 1 and l e , 2 in the evaporator two-phase zone and superheated zone, providing a deeper understanding of the refrigerant phase change process.

6. Conclusions

This paper utilizes Simulink without relying on M-files and S-functions to construct a simulation model for the evaporator, enabling the calculation of the lengths of the superheated and two-phase zones under specific evaporating pressure and temperature conditions, as well as the chilled water outlet temperature. Conversely, when the chilled water outlet temperature is known, the evaporating temperature and pressure can also be calculated.
To build the Simulink model for the entire refrigeration system, it is necessary to establish the coupling relationships among the evaporator, condenser, compressor, and expansion valve. The refrigerant state in the condenser includes the superheated zone, two-phase zone, and subcooled zone, where the superheated and subcooled zones belong to the single-phase zone. This is similar to the refrigerant state in the evaporator. Therefore, the modeling process and approach used for the evaporator model can be referenced for constructing a Simulink model for the condenser. Next, by inputting the evaporating pressure and refrigerant outlet temperature from the evaporator into the compressor and setting the inlet and outlet pressures of the expansion valve as the condensing pressure and evaporating pressure, respectively, the entire refrigeration unit’s Simulink model can be established. This model can be used to determine the continuous variation parameters, such as evaporating pressure, evaporating temperature, condensing pressure, condensing temperature, and the refrigerated cabin temperature. This facilitates the establishment of a comprehensive unit database for data analysis and decision evaluation.
The constructed model can be utilized for system fault analysis and diagnosis. By comparing historical data, if the chilled water outlet temperature decreases slowly, it may be indicative of issues within the evaporator system. Possible reasons could include blockages in the heat exchanger-chilled water pipes leading to reduced flow or scaling on the evaporator, in turn resulting in a decrease in the heat transfer coefficient. Further investigation can be conducted by integrating pipeline pressure with the simulated heat transfer coefficient to pinpoint the specific cause. If the pipeline pressure is normal, it could indicate a problem caused by scaling.
The model can also support energy-saving analysis and control of the system. Utilizing the established refrigeration system model, power consumption models for cooling water pump units, chilled water pump units, and refrigeration units can be created. By adjusting the frequency of pumps or compressors based on intelligent algorithms such as neural networks, the system’s power consumption can be minimized while satisfying requirements for refrigerant superheat, chilled water outlet temperature, and room temperature.
The modeling approach presented in this paper not only provides modeling ideas and a research platform for researchers, significantly facilitating their work and reducing their workload, but also makes it easy for beginners to grasp the principles of refrigeration system modeling, enabling them to quickly delve into more in-depth research. Moreover, the model can be applied for fault analysis, diagnosis, and decision evaluation.

Author Contributions

Conceptualization, Z.H. and Q.Z.; methodology, Q.Z. and Y.T.; software, X.L.; validation, Q.Z., X.L. and Z.W.; formal analysis, Z.H.; investigation, Z.W.; resources, X.W.; data curation, Z.W. and X.W.; writing—original draft preparation, Z.H.; writing—review and editing, Z.H.; visualization, Y.T.; supervision, J.Z.; project administration, J.Z.; funding acquisition, Z.H. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by [Ministry of Industry and Information Technology Project: Innovation Project of the Offshore LNG Equipment Industry Chain] grant number [CBG3N21-2-7] and the APC was funded by [Ministry of Industry and Information Technology Project: Innovation Project of the Offshore LNG Equipment Industry Chain] grant number [CBG3N21-2-7].

Data Availability Statement

Data are contained within the article.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Dey, P.; Mandal, B.K. Performance Enhancement of a Shell-and-tube Evaporator Using Al2O3/R600a Nanorefrigerant. Int. J. Heat Mass Transf. 2021, 170, 121015. [Google Scholar] [CrossRef]
  2. Chen, F.; Liao, S.; Gong, G. Finite Time Thermodynamic Analysis and Optimization of Water Cooled Multi-spilit Heat Pipe System (MSHPS). Energy Built Environ. 2022, 3, 373–383. [Google Scholar] [CrossRef]
  3. Moharana, S.; Bhattacharya, A.; Das, M.K. A critical review of parameters governing the boiling characteristics of tube bundle on shell side of two-phase shell and tube heat exchangers. Therm. Sci. Eng. Prog. 2022, 29, 101220. [Google Scholar] [CrossRef]
  4. Santa, R.; Bosnjakovic, M.; Cikic, A. Experimental and Numerical Testing of Heat Pump Evaporator. Appl. Sci. 2022, 12, 11973. [Google Scholar] [CrossRef]
  5. Xu, G. Simulation of Air Conditioning Systems and Room Thermal Environments Based on SimulinkFluent Co-Simulation. Master’s Thesis, Huazhong University of Science and Technology, Wuhan, China, 2021. [Google Scholar]
  6. Wang, Z.; Lu, X.; Wang, G. Experimental and simulation performance optimization of an exhaust air heat recovery heat pump unit. Int. J. Low-Carbon Technol. 2022, 17, 686–695. [Google Scholar] [CrossRef]
  7. El-Ghetany, H.H.; Omara, M.A.; Abdelaziz, G.B.; Abdelhady, R.G. Design of Silica Gel/Water Adsorption Chiller Powered by Solar Energy for Air Conditioning Applications. J. Energy Storage 2023, 63, 107055. [Google Scholar] [CrossRef]
  8. Sanama, C.; Xia, X. Control-oriented Modelling for Multivariable Control of Refrigerant Dynamics in Vapor Compression Systems. In Proceedings of the 2023 14th International Conference on Mechanical and Intelligent Manufacturing Technologies (ICMIMT), Cape Town, South Africa, 26–28 May 2023; pp. 332–344. [Google Scholar]
  9. Zhao, D.; Wu, Z.; Ding, G. Fast Calculation Method for Supercritical Refrigerant Thermodynamic Properties. J. Eng. Thermophys. 2008, V29, 1645–1648. [Google Scholar]
  10. Wu, Z. Rapid calculation of thermodynamic properties of refrigerants. J. Shanghai Jiao Tong Univ. 2006, 40. [Google Scholar] [CrossRef]
  11. Khaled, N.; Mathur, H. Modeling and Optimal Control of Vehicle Air Conditioning System. In Intelligent Manufacturing and Energy Sustainability: Proceedings of ICIMES 2020; Springer: Singapore, 2021; pp. 275–283. [Google Scholar]
  12. Zhang, Y. Parameter identification of evaporator models based on social particle swarm algorithm. Comput. Era 2022, 53–57. [Google Scholar] [CrossRef]
  13. Prajapati, M.; Shah, M.; Soni, B. Performance evaluation of a novel geothermal energy integrated single effect evaporator desalination with software simulation. J. Clean. Prod. 2023, 407, 137087. [Google Scholar] [CrossRef]
  14. Wallis, G.B. One-Dimensional Two-Phase Flow; McGraw-Hill Book Company: New York, NY, USA, 1969. [Google Scholar]
  15. Zhou, B.; Zhang, B. The Parameter Calculation of the Thermophysical Properties of HFC-134a. J. South. Yangtze Univ. (Nat. Sci. Ed.) 2005, 4, 390–393. [Google Scholar]
Figure 1. The structure of the paper.
Figure 1. The structure of the paper.
Processes 12 00379 g001
Figure 2. Pressure–enthalpy diagram of refrigeration cycle.
Figure 2. Pressure–enthalpy diagram of refrigeration cycle.
Processes 12 00379 g002
Figure 3. Changes in refrigerant specific enthalpy and chilled water temperature within the evaporator.
Figure 3. Changes in refrigerant specific enthalpy and chilled water temperature within the evaporator.
Processes 12 00379 g003
Figure 4. Overall simulation diagram for the superheated and two-phase zones.
Figure 4. Overall simulation diagram for the superheated and two-phase zones.
Processes 12 00379 g004
Figure 5. Model for calculating the chilled water outlet temperature, t w , e , 2 , of the superheated zone.
Figure 5. Model for calculating the chilled water outlet temperature, t w , e , 2 , of the superheated zone.
Processes 12 00379 g005
Figure 6. Model for the calculating refrigerant outlet temperature, t r , e , 3 , of the superheated zone.
Figure 6. Model for the calculating refrigerant outlet temperature, t r , e , 3 , of the superheated zone.
Processes 12 00379 g006
Figure 7. Calculation of B and the refrigerant outlet enthalpy, h r , e , 3 , of the superheated zone.
Figure 7. Calculation of B and the refrigerant outlet enthalpy, h r , e , 3 , of the superheated zone.
Processes 12 00379 g007
Figure 8. Model for calculating the evaporator heat transfer coefficient, K e , 1 , of the superheated zone.
Figure 8. Model for calculating the evaporator heat transfer coefficient, K e , 1 , of the superheated zone.
Processes 12 00379 g008
Figure 9. Model for calculating the chilled water heat transfer coefficient, a w , e , gs , of the superheated zone.
Figure 9. Model for calculating the chilled water heat transfer coefficient, a w , e , gs , of the superheated zone.
Processes 12 00379 g009
Figure 10. Model for calculating the refrigerant heat transfer coefficient, a r , e , gs , of the superheated zone.
Figure 10. Model for calculating the refrigerant heat transfer coefficient, a r , e , gs , of the superheated zone.
Processes 12 00379 g010
Figure 11. Model for calculating the refrigerant flow rate, u r , e , gs , kinematic viscosity, ν r , e , gs , and specific heat at constant pressure, c p - r , e , gs , of the superheated zone.
Figure 11. Model for calculating the refrigerant flow rate, u r , e , gs , kinematic viscosity, ν r , e , gs , and specific heat at constant pressure, c p - r , e , gs , of the superheated zone.
Processes 12 00379 g011
Figure 12. Model for calculating the refrigerant density, ρ r , e , gs , of the superheated zone.
Figure 12. Model for calculating the refrigerant density, ρ r , e , gs , of the superheated zone.
Processes 12 00379 g012
Figure 13. Experimental device connection diagram.
Figure 13. Experimental device connection diagram.
Processes 12 00379 g013
Figure 14. Chilled water inlet temperature.
Figure 14. Chilled water inlet temperature.
Processes 12 00379 g014
Figure 15. Chilled water outlet temperature.
Figure 15. Chilled water outlet temperature.
Processes 12 00379 g015
Figure 16. The relationship between the outlet temperature of chilled water and the inlet temperature and evaporation pressure.
Figure 16. The relationship between the outlet temperature of chilled water and the inlet temperature and evaporation pressure.
Processes 12 00379 g016
Figure 17. The relationship between evaporation pressure and the inlet and outlet temperatures of chilled water.
Figure 17. The relationship between evaporation pressure and the inlet and outlet temperatures of chilled water.
Processes 12 00379 g017
Figure 18. Temperature of the chilled water outlet with variation in the chilled water mass flow rate.
Figure 18. Temperature of the chilled water outlet with variation in the chilled water mass flow rate.
Processes 12 00379 g018
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

He, Z.; Zhang, Q.; Wei, Z.; Liao, X.; Wu, X.; Zhang, J.; Tan, Y. Modeling Method for Overheated Zone and Two-Phase Zone of Dry Shell-and-Tube Evaporator in Ship Air Conditioning. Processes 2024, 12, 379. https://doi.org/10.3390/pr12020379

AMA Style

He Z, Zhang Q, Wei Z, Liao X, Wu X, Zhang J, Tan Y. Modeling Method for Overheated Zone and Two-Phase Zone of Dry Shell-and-Tube Evaporator in Ship Air Conditioning. Processes. 2024; 12(2):379. https://doi.org/10.3390/pr12020379

Chicago/Turabian Style

He, Zhibin, Qi Zhang, Zhenghao Wei, Xingzhe Liao, Xiaoyu Wu, Jundong Zhang, and Yanghui Tan. 2024. "Modeling Method for Overheated Zone and Two-Phase Zone of Dry Shell-and-Tube Evaporator in Ship Air Conditioning" Processes 12, no. 2: 379. https://doi.org/10.3390/pr12020379

APA Style

He, Z., Zhang, Q., Wei, Z., Liao, X., Wu, X., Zhang, J., & Tan, Y. (2024). Modeling Method for Overheated Zone and Two-Phase Zone of Dry Shell-and-Tube Evaporator in Ship Air Conditioning. Processes, 12(2), 379. https://doi.org/10.3390/pr12020379

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop