Scientiﬁc Approaches to Solving the Problem of Joint Processes of Bubble Boiling of Refrigerant and Its Movement in a Heat Pump Heat Exchanger

: The joint processes of bubble boiling and refrigerant movement in heat pump tubes are considered. A coil located on the back side of the photovoltaic panel (PVP) is used as an additional heating surface. A mathematical model of the theoretical determination of the temperature on the front surface of the PVP cooled by freon in the coil from the back side has been created. A feature of the numerical simulation of bubble boiling of refrigerant in a coil was the insufﬁciently detailed study of the process and the lack of references to the results obtained in earlier work. The authors have clariﬁed the boundary conditions and assumptions for numerical simulation of the bubble boiling process of refrigerant R407C. The results were compared with the theoretically found values, and later used in the design of an energy complex consisting of a heat pump and a photovoltaic panel. The mathematical model of theoretical calculation of the PVP temperature and the methodology for constructing a numerical model of bubble boiling of refrigerant as part of the methodology for designing energy technology complexes based on the authors’ plan to present a uniﬁed methodology for energy technology complexes for desalination of seawater based on other types of renewable energy sources.


Introduction
To date, there are no universal methods for calculating the heat transfer coefficient or a universal description of the thermophysical picture of boiling in a liquid stream. However, many empirical methods have been created for various refrigerants under conditions of lowpressure values (up to 2-4 atm) and mass flow rates, and special formulas for calculating heat transfer to small channels (with a diameter of 1 mm and above) and to strictly specified pressure values have been derived [1][2][3][4][5].
Refrigerant R407C was chosen as the working medium for the study. Its thermodynamic properties are close to the common R22, while it does not contain chlorine, and is non-flammable and non-toxic.
Consider energy complexes that combine heat pumps and the use of solar energy. In the scientific paper [6], the authors developed a system based on solar energy for the production of fresh water and electricity. This system consists of a solar module, a solar-powered Rankine cycle, a heat storage subsystem, and a distillation device. Sea water used for desalination is heated by a saturated steam-water mixture from a steam turbine. Using the outlet liquid as a heat source for the distillation device eliminates the need for an external heat exchanger for condensation. The authors of [6] calculated the energy and exergy efficiency, power and fresh water production capacities for each Energies 2023, 16, 4405 2 of 16 element of the system. In addition, a parametric study was carried out to take into account environmental conditions. In research paper [7] the authors analyzed the energy complex for four modules, including a solar heliostat and a receiver, a heat storage device, a Rankine cycle and a desalination system. The steam generator transfers heat to the Rankine cycle. The authors of the study [8] present the balance equations for each element separately and for the entire system as a whole. A mathematical model for the utilization of waste heat from the Rankine cycle has been developed. However, the authors make a number of assumptions: there are no chemical reactions, the potential and kinetic energy practically do not change, the energy costs for the pump are negligible, except for operation in the Rankine cycle, and heat losses occur only through the surface of the reservoir with molten salt and the solar collector.
In [9], the heat loss rate of a solar receiver was calculated. In the research paper [10], after calculating the temperature difference, the authors analyzed the heat losses through the tank with molten salt. The authors of [10] developed a mathematical model of heat accumulation in a tank. This is a rather complicated model, since the volume of the salt melt inside the tank is constantly changing. In [11], a mathematical model of the seawater desalination process was constructed. The study [12] considers the energy and exergy efficiency of the system as a whole. This can be achieved using the main inputs and outputs, while the useful outputs of the system are fresh water and electricity, and the only source of energy is solar energy.
The methodology considered by the authors is a scientific approach to solving the problems of calculating efficiency and optimization in the design of energy technology complexes that combine a heat pump (including a turbo expander for the implementation of the organic Rankine cycle), a photovoltaic panel, a desalination plant, and optionally a solar concentrator and a wind generator. In addition, the scientific approach does not deny the possibility of not only working on the Rankine cycle, since the methodological base does not change. A special feature of the study was the method of heating refrigerant-an additional heater in the form of a coil installed on the back of the photovoltaic panel, Figure 1a,b. Figure 1c also shows the basic scheme of operation of the energy technology complex implemented by the authors. The article considers the first part of the installation-a heat pump.
Let us explain the diagram in Figure 1c. The working substance is refrigerant R407c, which moves along the contour of the system. An evaporator works to convert a gaseous substance into a liquid at low temperature and pressure. The compressor then compresses the liquid refrigerant, increasing the pressure to 19 bar and the temperature to 85 • C. After passing through the compressor, the refrigerant enters the heat exchanger (condenser).
Here it gives up its heat to the water, turning it back into a liquid form under pressure. Water enters the system through the filler neck and exits through the drain cock. After the condensed refrigerant passes through the expansion valve, the pressure drops significantly, causing some of the liquid to evaporate. This mixture of liquid and vapor goes to the evaporator. The liquid turns into vapor, evaporating under the action of the heat of the FEP or an additional heater. Due to this, boiling occurs in the evaporator, which consumes heat from the environment.  1-low-temperature heat source, 2-evaporator, 3-photovoltaic panel, 4-additional (backup) heater, 5-compressor, 6-condenser, 7-expander (capillary tube), 8-heated water, 9-valve, 10-main desalination tank, 11-steam compressor, 12-distillate heat exchanger-condenser, 13-desalinated water outlet, 14-salt water outlet.

Vapor-Liquid Mixture
The condensate is throttled before being fed to the evaporator. At the same time, along with a decrease in pressure, the enthalpy constancy (1) can be set with sufficient accuracy.   The condensate is throttled before being fed to the evaporator. At the same time, along with a decrease in pressure, the enthalpy constancy (1) can be set with sufficient accuracy. The throttling process in the P-i diagram for the selected refrigerant is shown in Figure 2 below in Process 3-4. The throttling process in the P-i diagram for the selected refrigerant is shown in Figure 2 below in Process 3-4. Given the necessary pressure drop and knowing the condensation and evaporation temperatures, it is possible to calculate the degree of dryness of steam at point 4, corresponding to the state of refrigerant at the inlet to the evaporator. The degree of dryness under the specified conditions is equal to x = 0.4. Thus, as a mandatory boundary condition, a fraction of 40% of the mass flow rate for steam is set.

Boiling
The calculation takes into account the flow rate in a pipe with a diameter of 5 mm. Developing the previous suggestion, a mix of bubble and foil boiling is assumed. When the liquid boils and evaporates, bubbles form on the walls, which quickly collapse and form a thin film of liquid when there is enough steam in the channel section.

Flow Mode
This flow mode corresponds to the form of a mixture in which steam forms the core of the flow, while the liquid phase flows in the shape of a foil on the surface of the wall. From the first assumption, it can be concluded that the vapor-liquid blend enters the evaporator in the emulsion mode, which is typical for the initial stage of vaporization. Due to the increase in steam consumption (up to a fraction of 0.4, as can be concluded from the diagram in Figure 2) bubbles tend to accumulate and stick together in the flow into larger ones.
It was assumed that the emulsion regime would prevail throughout the site until the formation of steam shells, that is, large steam plugs ( Figure 3a). As a result of computer modeling, a similar picture was obtained in the initial formulation of the problem. Figure   Figure 2. Heat pump cycle in the P-i diagram. Symbols and units of measurement: red linestemperature, • C, blue lines-entropy, kJ/(kg·K), green lines-specific volume, m 3 /kg, horizontallyenthalpy in kJ/kg, vertically-pressure in bar.
Given the necessary pressure drop and knowing the condensation and evaporation temperatures, it is possible to calculate the degree of dryness of steam at point 4, corresponding to the state of refrigerant at the inlet to the evaporator. The degree of dryness under the specified conditions is equal to x = 0.4. Thus, as a mandatory boundary condition, a fraction of 40% of the mass flow rate for steam is set.

Boiling
The calculation takes into account the flow rate in a pipe with a diameter of 5 mm. Developing the previous suggestion, a mix of bubble and foil boiling is assumed. When the liquid boils and evaporates, bubbles form on the walls, which quickly collapse and form a thin film of liquid when there is enough steam in the channel section.

Flow Mode
This flow mode corresponds to the form of a mixture in which steam forms the core of the flow, while the liquid phase flows in the shape of a foil on the surface of the wall. From the first assumption, it can be concluded that the vapor-liquid blend enters the evaporator in the emulsion mode, which is typical for the initial stage of vaporization. Due to the increase in steam consumption (up to a fraction of 0.4, as can be concluded from the diagram in Figure 2) bubbles tend to accumulate and stick together in the flow into larger ones.
It was assumed that the emulsion regime would prevail throughout the site until the formation of steam shells, that is, large steam plugs ( Figure 3a). As a result of computer modeling, a similar picture was obtained in the initial formulation of the problem.  shows the distribution of the proportion of liquid (red) and steam plugs (gradient to blue) in the tube section. Such a picture was achieved based on the considerations below. 3b shows the distribution of the proportion of liquid (red) and steam plugs (gradient to blue) in the tube section. Such a picture was achieved based on the considerations below. Intensive evaporation of liquid along the boundaries of dry spots (centers of vaporization) is taken as the determining mechanism of bubble boiling. A calculation that takes this factor into account explains the extremely high rate of heat transfer during bubble boiling. When studying a physical model and generalizing experimental data, the analytical Formulas (2) and (3) are adopted for the calculation of heat transfer: The thermal transfer at boiling in the liquid flow is calculated as a combination of convective heat exchange and heat exchange at boiling with the predominance of the latter (4): The final criterion similarity equation takes the form (5): Intensive evaporation of liquid along the boundaries of dry spots (centers of vaporization) is taken as the determining mechanism of bubble boiling. A calculation that takes this factor into account explains the extremely high rate of heat transfer during bubble boiling. When studying a physical model and generalizing experimental data, the analytical Formulas (2) and (3) are adopted for the calculation of heat transfer: The thermal transfer at boiling in the liquid flow is calculated as a combination of convective heat exchange and heat exchange at boiling with the predominance of the latter (4): The final criterion similarity equation takes the form (5): 2. Determination of the refrigerant temperature at the outlet of the evaporator. The first part is presented in the format of a flowchart in Figure 4 according to the Gauss methodology (with input of initial data, forward and reverse, and output of results). As a result, the temperature on the surface of the PVP is calculated with an error of less than 0.1%. The results with an ambient temperature of 20 • C and a wind of 1 m/s are summarized in Table 1. 2. Determination of the refrigerant temperature at the outlet of the evaporator. The first part is presented in the format of a flowchart in Figure 4 according to the Gauss methodology (with input of initial data, forward and reverse, and output of results). As a result, the temperature on the surface of the PVP is calculated with an error of less than 0.1%. The results with an ambient temperature of 20 °C and a wind of 1 m/s are summarized in Table 1.
Losses from above q above Silicon radiation losses Glass radiation losses q gl.rad Heat losses in the system Q losses q above + q sil.rad + q und. + q gl.rad W 11.516 Error rate ∆Q Q−Q losses Q · 100 % 0.067 A 3D diagram is constructed based on the calculation of the temperature of the PVP wall depending on two input arguments: wind speed and ambient temperature. The result is shown in Figure 5: A 3D diagram is constructed based on the calculation of the temperature of the PVP wall depending on two input arguments: wind speed and ambient temperature. The result is shown in Figure 5: The temperature on the evaporator wall is a boundary condition for simulating the refrigerant boiling process. For further theoretical calculation, it is a test value when calculating the thermal transmission ratio.
Theoretically, the approximate value of the thermal transmission ratio is calculated according to the Labuntsov formula: According to the reference books, the saturation temperature is determined to be TH under pressure p at point 4 in the diagram according to Figure 2. The physical properties of saturated liquid refrigerant are searched for at the obtained temperature according to reference data [14]. The Reynolds number is calculated at a given flow velocity ω (8): The result obtained is used to calculate the temperature on the surface of the evaporator tube (6): The temperature on the evaporator wall is a boundary condition for simulating the refrigerant boiling process. For further theoretical calculation, it is a test value when calculating the thermal transmission ratio.
Theoretically, the approximate value of the thermal transmission ratio is calculated according to the Labuntsov formula: According to the reference books, the saturation temperature is determined to be T H under pressure p at point 4 in the diagram according to Figure 2. The physical properties of saturated liquid refrigerant are searched for at the obtained temperature according to reference data [14]. The Reynolds number is calculated at a given flow velocity ω (8): The Nusselt number is calculated according to the Petukhov formula (it implies a developed turbulent regime that will arise due to the presence of turns of the evaporator and high flow velocity of the refrigerant) (9): Next, the convective thermal transmission and boiling heat transfer are compared by comparing the corresponding coefficients. The calculation of the thermal transmission ratio for the forced flow process in the tube is determined as follows (10): The thermal transmission ratio for refrigerant boiling is computed according to the equation corresponding to the geometric and thermophysical similarity. Its approximate calculation is performed by Formula (7) under the condition of a boiling pressure above 1 atm.
IF αboil/αconv. > 3, then the thermal transmission ratio during boiling in the pipe is assumed to be equal to αboil, otherwise-αconv.
As a verification condition for the correctness of the selection of the criterion equation, the surface temperature of the evaporator tube is calculated to be higher. The calculation is correct if the temperature value in this part coincides with the previous one (11): To determine the refrigerant temperature at the evaporator outlet, the amount of heat to bring the steam-water mixture to the state of saturated steam is determined, as well as for its overheating (12): Evaporation section length (13): The remaining section is a refrigerant superheater (14): The refrigerant outlet temperature is (15):  Table 2.

. Energy Equations
To apply solutions by the numerical method, the formulas are used:

The Transfer Equation
Generally used in models with one or two differential equations, the transfer equation can be presented in this form: The interpretation of the parameters for each type of equations is presented in Table 3. Table 3. Explanation of parameters.
The turbulence model is used in the calculations of transition k-kl-ω.

Modeling in a Software Environment ANSYS Fluent
Standard k-ε (+2E) is a two-parameter RANS model, which is very economical in terms of computational resources and is used for a wide class of flows. Wall functions are used and the model is suitable only for steady turbulent flows. The disadvantage is that in many problems there is a strong eddy diffusion. It should be noted that this model is able to take into account the effect of compressibility and natural convection, which is useful in the calculation of industrial cooling systems. However, this model is not suitable for solving the problem of refrigerant boiling, since it gives a large error and low convergence of results during boiling.
The authors came to the conclusion that transitional flow models should be used (the models only model the transition from laminar to turbulent flow; the wall function is used): Transition k-kl-ω (+3E) The k-kl-ω (+3E) transition turbulence model is a widely used turbulence model that combines the k-ε and k-ω models together with the Baldwin-Lomax model. This combination of models provides a reliable method for predicting flow behavior in transients (which is useful in calculating boiling due to its erratic nature and difficult prediction).
The k-kl-ω (+3E) transition model has a number of advantages over other turbulence models. It is able to accurately predict the onset of the transition from laminar to turbulent flow, which is a challenge for many turbulence models.
Transition k-kl-ω (+3E) is a three-parameter model, suitable for solving heat transfer problems in elements that create turbulence (for example, on a turbine blade).
However, the k-kl-ω (+3E) transition model is not without its drawbacks. The model is relatively complex and requires significant computational resources to solve. In addition, it is more prone to numerical instability than other simpler turbulence models, potentially making it difficult to use in some problems. Finally, the model relies on the use of a transition factor to another model, which can lead to inaccuracies and errors if the factor is set incorrectly.
Border conditions

1.
From the experience of modeling the emulsion flow, the condition of a two-phase flow at the inlet is set (previously, only one liquid phase was set at the inlet, and not a mixture); 2.
The extreme pressures in the cycle were determined (3.5 and 19 bar, respectively); 3.
The saturation temperature is set from the theoretical calculation;

4.
Based on points 2 and 3, a certain degree of dryness of the vapor-liquid mixture at the inlet to the evaporator is set.
Computer modeling of the flow process of refrigerant R407C in the channel of a brass evaporator with an external diameter of 10 mm and a wall width of 1 mm with a dispersed-ring flow is carried out. Figure 6a shows a geometric model. The simulation includes the condition of fluid flow in the gravity field with the formation of vapor bubbles. The group of problems of the effect of the heating wall on the evaporated liquid is solved with satisfactory accuracy by a non-stationary formulation using the Transition k-kl-ω turbulence model. This is a three-parameter model of the RANS group, suitable for a narrow range of tasks, which allows modeling the transformation from laminar to turbulent flow, in which wall functions are used. The main solvable equations embedded in the mathematical model are the energy equation and the turbulence and turbulence-dissipation model equations.
The lifting force (causing the mixing of the flow) occurs due to convection, described by the Boussinesq approximation. This approximation consists in solving the system of Navier-Stokes equations and the energy equation using a special dependence of density on temperature only when calculating the volumetric force.
The approximation allows us to take into account the gravity field, that is, the aspi- The search for unknowns in the equations is carried out using structured and bodyadapted prismatic grids (Figure 6b), which allows solving the problem when describing complicated geometry (sweep method for cylindrical bodies; at least 2 million prismatic and tetrahedral elements). A partition in the wall region (a layer of at least three cells) was also performed, which allows using the boundary layer theory for correct calculation of the flow at the walls of the tube, as well as correct calculation of the redistribution of the vapor and liquid phase.
The simulation includes the condition of fluid flow in the gravity field with the formation of vapor bubbles. The group of problems of the effect of the heating wall on the evaporated liquid is solved with satisfactory accuracy by a non-stationary formulation using the Transition k-kl-ω turbulence model. This is a three-parameter model of the RANS group, suitable for a narrow range of tasks, which allows modeling the transformation from laminar to turbulent flow, in which wall functions are used. The main solvable equations embedded in the mathematical model are the energy equation and the turbulence and turbulence-dissipation model equations.
The lifting force (causing the mixing of the flow) occurs due to convection, described by the Boussinesq approximation. This approximation consists in solving the system of Navier-Stokes equations and the energy equation using a special dependence of density on temperature only when calculating the volumetric force.
The approximation allows us to take into account the gravity field, that is, the aspiration of steam bubbles in the direction opposite to gravity, respectively, to take into account the orientation of the tube in space.  Figure 7a shows that due to the instability of the process, convergence proceeds slowly. However, Figure 7b-d shows the totals for temperatures and fractions of the liquid and vapor phases of refrigerant. The results are indeed close to theoretical and experimental data in accuracy.

Results
It should be noted that there are few research papers devoted to the joint processes of heat transfer, motion and boiling, both in the Russian Federation and in the world today. Basically, scientific laboratories use experimental data obtained earlier on stands and on-site objects. The complexity of research in the field of bubble boiling consists in the transience of processes (high intensity), and a large volume of newly obtained mixtures (for example, the same R407c is a non-azeotropic mixture of three refrigerants). In other words, there are a lot of new artificially produced mixtures, and full-scale work on numerical modeling and verification of at least heat transfer coefficients according to previously obtained formulas and ratios has not been carried out recently. The authors in their works try to partially compensate for the lack of "fresh" data that could be used by both design engineers and research engineers.

Comparison of the Received Data
In the course of the research, the authors were able to compare the temperatures on the back surface of the PVP: obtained experimentally, theoretically, and by numerical modeling. The discrepancy of the data turned out to be acceptable for engineering techniques. The closest results were obtained with theoretical calculation and numerical modeling. The difference was 0.1 • C. At the same time, when compared with experimental data on thermal sensors, the discrepancy is already more significant, approximately 0.3 • C, if we take into account the error of measuring instruments-0.5 • C. However, even with such values of temperature deviations, it is necessary to understand the order of the measured values-from −10 to + 90 • C. That is, with a measurement range of 100 • C, the error turns out to be insignificant.
Thus, it can be said with confidence that the method of determining the temperature of the PVP and the data introduced into the boundary conditions during numerical modeling open up new opportunities for designers of energy technology complexes based on renewable energy sources, as new scientific approaches to temperature prediction appear, and consequently, more accurate determination of the design features of equipment, in particular, heat exchangers.  It should be noted that there are few research papers devoted to the joint processes of heat transfer, motion and boiling, both in the Russian Federation and in the world today. Basically, scientific laboratories use experimental data obtained earlier on stands and on-site objects. The complexity of research in the field of bubble boiling consists in the transience of processes (high intensity), and a large volume of newly obtained mixtures (for example, the same R407c is a non-azeotropic mixture of three refrigerants). In Numerical modeling and Formula (15) compared the obtained data. The discrepancy was about 0.1-0.5 • C.
If we compare the results in Tables 2 and 4, we can conclude that numerical simulation gives results that are closer to the theoretical calculation data, although, for the analytical calculation, the convergence turned out to be quite good.

Limitations of the Study
Regarding the limitations of the flow regime under study, it should be noted that after optimizing the mathematical model of heat transfer to the evaporator tube and a number of numerical calculations, it was found that the flow regime in the evaporator is dispersed-annular. Steam in the emulsion mixture is redistributed so that the main phase (liquid with 60% of the flow rate) moves along the pipe surface and in the form of small drops in steam, that is, in an annular mode.
The authors cannot guarantee the applicability of the obtained results for other ratios of the vapor and liquid phases.

Prospects for the Use of a Heat Pump in an Energy Complex for Seawater Desalination
The goal of creating plants with a thermal process for obtaining fresh water from sea water has led to the development of innovative concepts and design schemes with high efficiency. Particular attention is paid to the issues of optimizing operating modes, replacing expensive structural materials, pipe surfaces with other heater configurations, and reducing scale formation by an increase in the initial temperature of heating demineralized water in order to reduce the consumption of thermal and electrical energy. This goal can be achieved through the utilization of low-grade and waste heat. The cost of water produced in the desalination process is largely determined by the energy source used in the construction of large distillation plants [6][7][8][9][10][11][12].
In the future, the authors will develop a methodology for calculating the exergy fluxes of a desalination plant. An organic Rankine cycle with a turbo expander is used instead of a heat pump for the advanced scheme of the energy complex as a whole, and additional components for the use of renewable energy sources, such as a wind turbine and a solar concentrator, can be implemented. This approach allows for an integrated approach to energy conversion, using multiple sources for maximum efficiency.
First of all, new achievements in the field of low-temperature technology will be considered; for example, the authors of articles [15,16] explore the properties of refrigerant R407c. At the same time, such installations, which are energy technological complexes, must be environmentally friendly. The authors currently have scientific developments in this area. For example, in the article [17], the basics of constructing mathematical models for trigeneration and multigeneration complexes are considered.

Conclusions
When optimizing the mathematical model of heat transfer to the evaporator tube and a number of numerical calculations, it was found that the flow regime in the evaporator is dispersed-annular. The steam in the emulsion mixture is redistributed so that the main phase (liquid with a fraction of 60% of the flow rate) moves along the surface of the tube in the form of small droplets in steam-in an annular mode ( Figure 2).
As the two-phase flow moves, there is a change to the opposite extreme case, when steam flows along the wall, and the liquid moves in the form of small droplets-in the dispersed mode ( Figure 2). The wet steam reaches the saturation state, and then the dry steam overheats (process 4-1 in Figure 1).
It is also revealed that in order to obtain the correct refrigerant heat transfer coefficient for a coil type evaporator, it is necessary to derive a criterion similarity equation based on experimental data.