{'fraction_non_alphanumeric': 0.025880860845784948, 'fraction_numerical': 0.006063968927113471, 'mean_word_length': 4.882956400093262, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
DEVELOPMENT OF A NUMERICAL TOOL FOR DYNAMIC SIMULATIONS OF TWO-PHASE COOLING SYSTEMS

In this article, the development of a simulation tool for two-phase cooling systems is discussed. The targeted application is modelling of the two-phase pumped loop systems used for Silicon tracking detector cooling at CERN. These systems are similar to vapour compression systems in that the thermal dynamics of such systems are dictated by the two-phase fluid present inside the heat exchangers. To properly account for such dynamics, non-homogenous void fraction based two-phase flow models (used for accurate modeling of vapour compression systems) have been incorporated. The tool has been validated against measurements taken for an R-410A-based residential heat pump unit. Both the heating and cooling mode have been simulated and the results have been compared against measured data. The simulated transients are found to compare well against measured trends. The simulations proceed faster than real-time. The tool shows readiness for use in the design of future detector cooling systems.


INTRODUCTION
In high-energy physics experiments, Silicon tracking detectors are used to uniquely identify fundamental particles created in particle colliders. The Silicon sensors that constitute the tracking detector are sources of heat that require an active cooling system to regulate their temperature during operation. The selected cooling systems must ensure high stability in time and small spatial thermal gradients. Cooling systems based on the Two-Phase Accumulator Controlled Loop (2PACL) concept have been adopted at CERN as the general approach in dealing with these challenges [1]. 2PACL is a specialisation of the more general class of Mechanically Pumped Two-phase Loops, where the saturation temperature in the evaporators is controlled by an Accumulator. This is a large-volume vessel always containing both phases of the refrigerant, whose pressure can be controlled by combined heating and cooling actions. A detailed description of 2PACL systems and their operating principle can be found in Petagna et al. [2], and is omitted here for the sake of brevity. Over the next five years, a large number of such 2PACL systems are required to be developed for installation in several of the experiments at CERN, with an installed cooling power of the order of 1 MW.
1D system-level dynamic simulations can be of great benefit to the designers of thermofluid systems through their ability to predict hypothetical scenarios in a safe and timeefficient manner. In recent years, numerical modelling and simulation have been applied in a wide variety of thermofluid domains and timescales [3][4][5]. The design and commissioning of future tracker cooling systems will require the availability of a similar dynamic simulation tool, and this motivated the present study.
The dynamics of these 2PACL systems share similarities with vapour compression systems commonly used in air-conditioning and refrigeration equipment. The thermal dynamics in both systems are dictated by the two-phase evaporation and condensation of the working fluid inside the heat exchangers. For such systems, the two-phase state must be properly accounted for in modelling exercises.
The simplest method for modelling two-phase flow is the homogeneous model, which adopts the assumptions of thermodynamic equilibrium in the control volume and equal flow velocities of both fluid phases. These assumptions, however, are not valid in all cases. Previous studies [6,7] have demonstrated that the homogeneous flow assumption can lead to faster transients in vapour compression systems, and a prediction of a lower system charge due to over-prediction of the void fraction (the fraction of the channel cross-sectional area occupied by the vapour phase, a key parameter in the analysis of two-phase flows). Thus, models incorporating more sophisticated void fraction models are required.
An alternative method involves the relaxation of the assumption of constant phase velocities adopted in the homogeneous model. In such models, a non-homogeneous voidfraction correlation is used to correct for the fluid density in two-phase flow [7]. This method is attractive since it does not require a significant modification of the form of the governing equations, and has been adopted in this study.
At CERN, the physical modelling platform EcosimPro [8] has previously been used for performing dynamic simulations of cryogenic systems [9]. The successful experience in the cryogenics domain has resulted in EcosimPro being selected at CERN for simulating the silicon tracker cooling systems as well. However, the cryogenics component library adopts the homogeneous assumption for two-phase fluid and, as previously noted, this assumption is inadequate for pumped loop cooling systems.
In this study, a numerical tool is developed in EcosimPro to handle the challenges posed by two-phase cooling systems. The tool has been used to model an R-410A based residential heat pump unit. A detailed set of high-quality measured data is available for the heat pump in two different operating modes. A comparison with these data allows validations of the developed tool against a real system. Furthermore, the simulations are also compared against a model developed on the Modelica/Dymola platform (typically used for physical modelling of energy systems [10]). The Dymola model adopts the homogeneous assumption for twophase flow. Thus, a comparison against it can highlight the impact of using the void-fraction method discussed previously, while also comparing against a well-established tool.

MATHEMATICAL MODELLING
Three partial differential equations (PDEs) dictate the flow of fluids in pipes: the mass conservation (continuity) equation, the momentum conservation equation, and the energy conservation equation. In applying these conservation equations, several assumptions (widely used for the dynamic simulations of thermofluid systems) have been made in the current study:  The working fluid is Newtonian.  The flow is 1-dimensional.  Heat conduction in the flow direction is neglected.  Kinetic and potential energy variations are not taken into account.  Viscous heating is neglected.  Gravitational effects are neglected (these effects can easily be accounted for at a later stage as necessary).
With these assumptions, Eqs. (1), (2) and (3) describe the mass, momentum and energy conservation equations, as described by Levy [11]: In system-level simulations, these equations can be discretized and enforced in a variety of forms such as lumped parameter models, moving boundary models or distributed parameter models. Distributed parameter methods, such as finite difference, finite element or finite volume methods, are normally preferred because they can express a high-level of modelling detail. In particular, the staggered grid scheme [12] was used here to decouple the momentum equation from the continuity and energy conservation equations. In this method, the momentum equation is solved on the momentum grid, which is offset by half a cell width from the thermal grid [13]. The mass and energy balance equations are solved on the thermal grid. Fig. 1 shows the discretization and the indexed grids. It can be seen that the boundaries of the momentum grid align with the centres of the thermal grid and vice versa. This eliminates the need to interpolate the parameters between the cells.

Two-phase flow governing equations
The thermodynamic parameters of pressure and specific enthalpy have been selected as the state variables. Specific enthalpy provides an intuitive understanding of the fluid phase, and the selection of these state variables enables an easier description of the two-phase conservation equations. The discretization procedure is described in detail in Winkler [14], and the discretized mass and energy equations for a given control volume are shown in Eqs. (4) and (5): The partial derivatives of density with respect to pressure and specific enthalpy have been evaluated for the two-phase region using the method described by Thorade and Saadat [15].
Here, H flow refers to the enthalpy flow rate, which has been calculated using the upwind scheme, as shown in Eq. (6): Zhang et al. [16] compared the three forms of the momentum equation for residential heat pump systems: a) no pressure drop, b) steady-state momentum equation, and c) full transient form of the momentum equation. The authors concluded that, in the cases of both fast and slow transients, the results for both the static and transient momentum equations are nearly the same. This indicates that time derivatives do not play a major role in pressure transients. Similar trends are expected in the case of the pumped loop systems since the physics governing the component models are similar. Therefore, the steady state form of the momentum equation has been adopted here, and the pressure drop has been calculated using nominal values provided as inputs to the model, as shown in Eq. (7), and derived from the article by Wetter et al. [17]: where ṁ 0 and dP 0 are the nominal parameters that are derived from the steady state pressure drop conditions. The flow-weighted enthalpy has been calculated as a function of the density-weighted enthalpy using an enthalpy correction term. The first step is the calculation of the slip ratio. In this study, the method of Zivi [18] has been adopted, wherein the slip ratio is calculated as per Eq. (8). Then, the slip ratio is substituted into Eq. (9) to solve for the flow vapour quality. The flow vapour quality is then used to calculate the enthalpy correction term, and finally the flow-weighted enthalpy, as per Eqs. (10) and (11): ĥ = ℎ + ℎ (11) where K = 0.4 is an empirical constant [9], x̂ is the flow-weighted vapour quality, x is the static vapour quality, ĥ represents the flow-weighted specific enthalpy, and h is the densityweighted enthalpy.
The advantage of this method is that, when the flow is single phase, the enthalpy correction term resolves to zero and the flow-weighted and density-weighted enthalpies have the same values. Therefore, the same equations of the mass and energy balance can be used for both cases.

Numerical treatment
Dermont et al. [19] highlight the fact that simulations involving reversible or zero mass flow rates are challenging for numerical solvers to handle since the solver is required to take very small time steps to ensure convergence. The authors recommend a series of techniques for more robust modelling of such cases. Some of their recommendations have been adopted in the current library, namely: 1. Pressure drop regularization: Since the mass flow rate and the pressure drop are correlated with a quadratic equation, the mass flow rate vs. pressure drop curve has an infinite slope at the origin, as shown in Fig. 2 a ('Default' continuous line). This causes oscillations to appear in the mass flow calculation, due to the fast dynamics near the zero region. The small time constants in this low flow region cause the solver to reduce the step size drastically. In order to alleviate this issue, a monotonically increasing third-order polynomial [20] has been used to approximate the curve in the small flow region, as illustrated in Fig. 2 a ('Regularized' dashed line). This ensures that the curve has a finite slope near the zero region, and the integrator is not forced to take extremely small time steps to maintain convergence. The error introduced is of second-order importance and does not affect the quality of the predictions. 2. Heat transfer coefficient smoothing: The heat transfer coefficients vary significantly between the single-phase and two-phase regions. This discontinuous change of the heat transfer coefficient can lead to numerical issues when a control volume transitions from one phase to the other. This issue can be circumvented through the application of a smoothing technique between the different phases, shown in Fig. 2 b. In this study, in particular, a spline interpolation method has been used to smooth the heat transfer coefficient. The liquid, twophase and vapour heat transfer coefficients are specified as user inputs. Then, in the region around the single-phase to two-phase flow transition (0 < x < 0.1, and 0.9 < x < 1), the spline interpolation is used between the coefficients, as schematically shown in Fig. 2 b. This ensures a smooth and gradual variation of the heat transfer coefficient around the transition between single-phase and two-phase, while again the error introduced is of second-order importance and does not affect the quality of the predictions. The individual component models for the various components are derived using the modelling techniques described in this section, and the individual details are omitted here for the sake of brevity. Further details of the models may be found in the paper by Ling et al. [21].

RESIDENTIAL HEAT PUMP
The test unit for this study is a 3-ton, refrigerant R-410A heat pump for air-conditioning and space-heating applications. It was used in a previous study [22] to investigate the impact of replacing R-410A with lower global warming potential refrigerants (so-called 'drop-in' replacement). The indoor and outdoor units are round-tube, plate-fin heat exchangers. A single-stage scroll compressor circulates the refrigerant around the heat pump. In the cooling mode, a Thermostatic Expansion Valve (TXV) is used as the expansion device while in the heating mode a short-tube orifice with a fixed opening diameter is used. A reversing valve switches the flow direction. In this study, only results involving R-410A as the working fluid are considered.
A schematic of the system is shown in Fig. 3. The outdoor unit was kept inside an environmental chamber with controlled temperature and humidity. Long pipes transported the refrigerant between the indoor and outdoor units, and the indoor unit was placed inside a closed air loop. The cyclic test conditions are derived from the ASHRAE Standard 116 [23], and are listed in Table I. The D test represents operation in the cooling mode, while the High Temperature Cyclic test represents heating mode. Further details of the test setup can be found in Alabdulkarem et al. [22].

Model implementation
The heat pump model has been implemented in the EcosimPro physical modelling platform. The platform incorporates object-oriented modelling constructs as well as acausal modelling which facilitate the task of the modeller. In contrast to causal modelling techniques, acausal modelling allows for equations to be stated without consideration of computational order. This has the benefit of allowing the programmed equations to more closely resemble their mathematical form when written on paper. The object-oriented approach allows for the development of components as aggregates of subcomponents. This enables easier addition or removal of modelling capabilities, such as accounting for gravitational head or including detailed pressure drop predictions. The models have been implemented in EcosimPro version 5.6.0. The calculations are performed on a 64-bit Windows 7 workstation with an Intel i7-6700 CPU and 8 GB of RAM. A screenshot of the EcosimPro model is shown in Fig. 4.

RESULTS AND DISCUSSION
The simulated models were compared against measured data, as well as against a component library developed previously in Dymola and used for simulating the heat pump under investigation. The Dymola simulations were performed using similar component models, although a homogeneous assumption was adopted for the two-phase region. Thus, a comparison of the results from the two different models will highlight the impact of accounting for the two-phase flow using the void-fraction based method discussed previously.

Heating mode results
The heating mode results are shown in Fig. 5. During the off-period before startup, the refrigerant accumulates in the coldest parts of the system. In the heating case, the outdoor ambient conditions are colder than indoors, thus the refrigerant collects in the outdoor components. Specifically, a significant amount is stored in the accumulator. Just before startup, the refrigerant is at 10.3 bar, corresponding to the saturation pressure of R-410A at 8.3°C. Upon startup, the compressor motor ramps up to its maximum rotational speed of 3500 RPM. The refrigerant rapidly migrates from the low pressure side to the high pressure side of the system. This causes a sharp decrease in the suction pressure and a corresponding increase in the discharge pressure for the first two minutes. The compressor mass flow rate quickly peaks, as evident in Fig. 5 e. As refrigerant continues to migrate to the high pressure side, the density of the refrigerant in the low-pressure side decreases. Thus, the density of the vapour entering the compressor drops, causing a decrease in the mass flow rate after the initial maxima. The refrigerant flow through the valve takes longer to establish in comparison to the compressor. Two phenomena dictate the valve flow rate: the pressure difference across the valve and the inlet liquid density. Initially, the indoor coil does not fully condense the refrigerant to subcooled state, and thus the inlet density is low. Furthermore, the pressure difference across the valve is small, and thus the flow through the valve is correspondingly small. Once the mass flow rate through the valve starts to increase, however, the rate of pressure change slows down. It is seen in Fig. 5 a that, in six minutes, the system does not attain steady state and the pressures are still evolving. A comparison of the temperatures is given in Fig. 5 b. Upon startup, the discharge temperature increases sharply to around 60°C and then decreases slightly. After shutdown, the discharge temperature in the test setup is seen to decrease instantaneously down to the saturation temperature before increasing. This is attributed to the static head created by the indoor unit (condenser) being placed at a higher level than the compressor. After shutdown, the refrigerant flow through the compressor ceases immediately. Thus, the liquid within the condenser migrates in the reverse direction as the compressor driven flow is due to gravity. The models do not account for gravity and thus cannot capture this phenomenon. Instead, the simulations show a gradual decay of temperature approaching the ambient conditions as the compressor continues to lose heat through natural convection and radiation. The air outlet temperatures are shown in Fig. 5 c for both the indoor and the outdoor coils and match the experimental trends well during the on-period. During the off-period, the outdoor unit fan is turned off, while the indoor unit air (circulated using the closed air loop blower) stays on. The off-cycle indoor coil air temperature predictions for both the platforms evolve faster than measured data. This is due to the models not incorporating the thermal inertia of the closed air loop, which continues to retain heat. The heating up of the outdoor coil air is also not as significant in the simulations as in the experimental data, and this is due to the simulations not accounting for the fact that heat is transferred from the hot compressor located inside the outdoor unit.
The compressor power consumption and the indoor unit air side capacities are shown in Fig. 5 d. These parameters help to quantify the performance of the heat pump, and can be used to investigate the cycling losses in the system. The cyclic capacity is predicted fairly well by EcosimPro, although the Dymola model overpredicts the air side capacity. This overprediction is a direct consequence of the faster startup transients caused due to the homogenous model assumption in the Dymola model.
The refrigerant charge distribution is shown in Fig. 5 f. Note that only the simulated data are shown, since measuring the charge in individual components is not trivially accomplished. Only the overall system charge can be measured while filling the plant. Indeed, a significant benefit of conducting such simulations is the ability to predict hard-to-measure parameters such as charge distribution. The plot shows that during the on-period, nearly a quarter of the system charge resides in the piping (primarily in the long liquid lines), indicating the importance of accounting for the piping volumes. After the off period, the refrigerant initially migrates towards the outdoor unit and eventually, only vapour is present in the indoor unit. After about 8 minutes, the charge shifts into the accumulator and the migration stabilises at around ten minutes.
The off-cycle trends are less well predicted than the on-cycle trend. Specifically, the refrigerant pressures equalize quicker in the simulations than in the test setup. Once the system power is turned off, the refrigerant charge migrates to the accumulator. In the real plant, however, the accumulator resides inside the Outdoor unit, and is placed alongside the hot compressor. This is seen as the rise in the temperature of the outdoor coil in Fig. 5 c after shutdown, caused by the heat radiation from the hot compressor. Thus, the migration to the accumulator is slower. Further, the indoor unit is placed at a much higher elevation than the outdoor unit but gravitational effects are neglected in the current model. These phenomena cause some differences in the off-period trends. Overall, however, EcosimPro and Dymola predictions compare well, with the Dymola models showing faster-than-real transients due to the adoption of the homogeneous model.

Cooling mode results
The cooling mode results are shown in Fig. 6. It can be observed that some of the system characteristics, such as the refrigerant temperatures and pressures, largely resemble the heating mode. The differences in transients between the two modes arise mainly due to the presence of a Thermostatic Expansion Valve (TXV) in the cooling mode, which actively controls the evaporator superheating. The mass flow rate through the TXV, shown in Fig. 6 e, reaches a maximum in the initial stages and then decreases. This is caused by the thermal inertia of the sensor bulb of the valve. In the initial stages, the refrigerant rapidly migrates from the evaporator, reducing the suction pressure. However, the sensor bulb pressure (dictated by the bulb temperature), takes time to adapt due to the thermal resistances of the piping and the sensor bulb as well as the contact between them. Thus, a pressure imbalance is created on the diaphragm that forces the valve to remain open. Once the bulb begins to respond, however, the valve begins to control the superheat, and a subsequent reduction in mass flow rate is seen. The differences in discharge temperature shown in Fig. 6 b occur due to differences in the compressor model between the libraries. The Dymola library incorporates a more simplified description for the suction and discharge chamber models for the compressor, and is, therefore, unable to predict the discharge temperature to the same accuracy of the EcosimPro model. The EcosimPro model simulates the Cooling mode cycle in 1324 seconds. The simulations thus achieve the objective of proceeding faster than real time. Finally, the system charge (mass of refrigerant) in the actual unit was 5.2 kg. The charge inventory in the simulated EcosimPro models was 3.69 kg, which is 71 % of the system charge. By contrast, the Dymola models, based on the homogeneous assumption, could only account for 52 % of the system charge.
It is, thus, seen that using the slip-ratio based void-fraction modelling method enables the accounting of a greater fraction of actual system charge, as well as a more accurate prediction of system transients, without gravely affecting simulation speed.

CONCLUSION
In this article, a simulation tool has been developed in the physical modelling platform EcosimPro for the purposes of modelling two-phase cooling systems with intended use in high energy physics applications. Slip-ratio based void-fraction correlations have been adopted to account for two-phase fluid flow. The method used requires minimal modification of the single-phase flow governing equations. The tool has been used to simulate a residential heat pump unit, which has similar thermal dynamics as the pumped loop systems used in high energy physics. Both heating and cooling modes of the heat pump are simulated and the simulation results are compared with measured data, as well as against a different library previously developed in Dymola. The Dymola model adopts the simpler homogeneous assumption for the two-phase region. The EcosimPro model is able to account for a much greater portion of the system refrigerant charge compared to Dymola model, highlighting the benefit of using void-fraction correlations for the two-phase regime.
The numerical tool shows readiness for tackling the challenges of simulating the dynamics of two-phase cooling systems in future high energy physics applications. It will next be used for assisting in the design and commissioning phases of the next generation of the pumped loop cooling systems used at CERN.