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
and
, are continuously variable during system operation and are represented by
and
, 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
,
, and
, and the temperatures are
,
, and
, respectively. The water temperature is represented as
,
, and
, denoted as
,
, and
, respectively. The tube wall temperatures at points 5 and 6 are denoted as
and
, with the average temperature,
, represented by
.
and
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, , the heat transfer area, , and the average temperature difference, , i.e., . “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]:
is the refrigerant mass flow rate (kg/s); is the chilled water mass flow rate (kg/s); and are the inlet and outlet enthalpy values of the refrigerant (kJ/kg); is the overall heat transfer coefficient of the condenser (kW/(m2·K)); is the total surface area of the tube (m2); and 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
is based on the relationship between saturated vapor enthalpy and temperature for R134a:
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:
where,
,
,
,
,
,
,
,
,
,
The formula for the heat exchanger’s heat transfer area and heat transfer coefficient in Equation (2) is as follows:
is the heat transfer coefficient on the refrigerant side of the superheated zone (kW/(m2·K)); is the heat transfer coefficient on the water side of the superheated zone (kW/(m2·K)); is the flow velocity of chilled water inside the tube (m/s); is the outer tube diameter (m); is the inner tube diameter (m); is the thermal conductivity of the tube wall (kW/(m·K)); and is the number of tubes.
The determination of the heat transfer coefficient in the single-phase zone can be obtained as follows:
where,
,
is the heat transfer coefficient on the refrigerant side of the superheated zone, (kW/(m2·K)); is the refrigerant flow velocity (m/s); is the kinematic viscosity of refrigerant (m2/s); is the density of refrigerant (kg/m3); is the specific heat capacity of refrigerant at constant pressure (kJ/(kg·°C)); and 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:
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:
- (2)
Calculation method for the overall heat transfer coefficient in the formula:
- (3)
Calculation method for the surface heat transfer coefficient on the chilled water side:
- (4)
Calculation method for the surface heat transfer coefficient on the refrigerant side [
14,
15]:
where,
,
,
,
,
.
is the average density of the refrigerant in the two-phase zone (kg/m3); is the average thermal conductivity of the refrigerant in the two-phase zone (kW/(m·K)); is the average kinematic viscosity of the refrigerant in the two-phase zone (m2/s); is the latent heat of vaporization of the refrigerant (kJ/kg); 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:
where,
,
,
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:
where,
,
,
.
- (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:
- (6)
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,
, mass flow rate,
, refrigerant mass flow rate,
, and evaporation temperature,
. The relationship expressions for the two unknowns, namely, the outlet temperature of chilled water,
, and refrigerant outlet temperature,
, in the superheated zone, are solved through the function
. Then, in
Figure 6,
is taken as a known parameter input, allowing the solution of these two unknowns. Subsequently, in the two-phase zone,
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,
.
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, , 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, , of 12 mm and an outer diameter, , of 16 mm, while the chilled water mass flow rate is set to 0.5 kg/s. For the condenser, the designed heat exchange capacity is 35 kW, with an inner diameter, , of 14 mm and an outer diameter, , of 18 mm, while the cooling water mass flow rate, , 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 and 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.