Dynamic viscosity of low GWP refrigerants in the liquid phase: An empirical equation and an artificial neural network

This study presents a simple correlation for describing the temperature and pressure dependence of the liquid dynamic viscosity of low GWP refrigerants, namely HydroFluoroOlefins (HFOs) and HydroChloroFluoroOlefins (HCFOs). The model has 3 input parameters (i.e., reduced temperature, reduced pressure, and acentric factor) and 6 coefficients which were regressed on 794 experimental data collated from the literature for 7 alternative refrigerants (i.e., R1233zd(E), R1234yf, R1234ze(E), R1234ze(Z), R1224yd(Z), R1336mzz(E), and R1336mzz (Z)). Moreover, a multi-layer perceptron neural network for the liquid dynamic viscosity of the studied fluids was developed from the selected dataset. The artificial network has the same 3 input parameters of the correlation and one hidden layer with 19 neurons. The results of the proposed correlation proved that it is an accurate model for calculating the dynamic viscosity of the studied liquid refrigerants, despite its simplicity. It ensured an average absolute relative deviation of the liquid dynamic viscosity (AARD( η )) of 2.88 %, lower than that given by other literature correlations. As expected, the multi-layer perceptron neural network provided the best results for all the selected refrigerants (AARD( η ) = 0.86 % for the complete dataset), proving that it can be considered a reference for the development of other models.


Introduction
Over the years, to decrease global greenhouse gas (GHG) emissions, various international environmental agreements and regulations (e.g., Regulation (EU) No. 517/2014No. 517/ (2014) ) and UNEP ( 2016)) have been approved to limit the use and production of several conventional long-lived GHG refrigerants, such as hydrofluorocarbons (HFCs).Consequently, a search for low global warming potential (GWP) alternatives, known as "fourth-generation" refrigerants (McLinden and Huber, 2020), is currently underway at a worldwide level.In particular, McLinden and his co-workers (Domanski et al., 2017;McLinden et al., 2017;McLinden and Huber, 2020) systematically searched for potential low GWP refrigerants by thermodynamic, environmental, and safety selection criteria.As a result, they found a limited amount of synthetic alternative working fluids, mostly hydrofluoroolefins (HFOs) and hydrochlorofluoroolefins (HCFOs), having suitable characteristics for various HVAC&R applications.However, a trade-off between suitable thermophysical properties and flammability is required for some of them.Moreover, few experimental data for the thermodynamic and transport properties of the selected potential alternatives are available in the literature (Bobbo et al., 2018).Consequently, models to accurately describe their properties are necessary.
Among the thermophysical properties required to design and simulate refrigeration systems, the dynamic viscosity (η) of refrigerants is essential to accurately evaluate the pressure drop and heat and mass transfer in system components, such as compressors and heat exchangers.But comprehensive theoretical-based models cannot provide accurate results due to the complex nature of the liquid η.For this reason, empirical estimation methods and semi-empirical correlations for the liquid η of fluids have been proposed in the literature (Poling et al., 2001;Viswanath et al., 2007).Several models based on various approaches and characterized by different levels of complexity have been specifically developed to describe the viscosity behavior of refrigerants and their mixtures, especially in the liquid phase.Many models have coefficients specific to each refrigerant that should be regressed from experimental data.Some of the main estimation methods for refrigerants and their mixtures are as follows: extended corresponding states (ECS) models used in REFPROP 10.0 (Huber, 2018;Lemmon et al., 2018); models taking into account the residual entropy scaling (Bell et al., 2016;Kang et al., 2022;Liu et al., 2020;Yang et al., 2021); semi-empirical correlations based on various theories (Assael et al., 1995;Latini et al., 2002;Liu et al., 2022;Teja et al., 1999;Wang et al., 2007;Yousefi et al., 2019); models based on equations of state (EoSs) (Khosharay et al., 2018); artificial neural network (ANN) models (Di Nicola et al., 2022;Ghaderi et al., 2013;Wang et al., 2020;Zhi et al., 2018).Some of the simplest models requiring readily available and well-known input parameters are briefly described below.
Latini and his co-authors (Latini et al., 2002(Latini et al., , 1996(Latini et al., , 1990) developed simple semi-empirical correlations to estimate the liquid dynamic viscosity dependence on the temperature at saturated conditions for various refrigerant families and their mixtures.The proposed models were obtained by the combination of the Batschinski equation (Batschinski, 1913) modified by Hildebrand (Hildebrand, 1977) with the rule of Mathias (Reid et al., 1987).Their latest form is the following: with where η is in mPa s, T r = T T c − 1 is the reduced temperature, T is the temperature in K, T c is the critical temperature in K, M is the molecular mass in kg kmol − 1 , T Br = T B T c − 1 is the reduced normal boiling point temperature, T B is the boiling point temperature in K, and C h, α, β are regressed coefficients.
In the original works, Eqs. ( 1) and (2) coefficients were calculated for the methane and ethane refrigerant families, yielding deviations of η generally lower than 6 %.A modified version of Eq. ( 1) has been recently proposed by Di Nicola et al. (2022) to model the dynamic viscosity dependence on temperature and pressure for low GWP refrigerants in the liquid phase.The modified model has the following expression: where p r = p p c − 1 is the reduced pressure, p is the pressure in MPa, p c is the critical pressure in MPa, D and E are regressed coefficients, and A has the expression of Eq. (2).By minimizing the deviations between the experimental and calculated data, the following values of the coefficients were obtained for the analyzed low GWP refrigerants: C = 1.376, h = 1.387, α = − 0.05275, β = − 4.9, D = 0.9758, and E = 1.139.
They ensured an average absolute relative deviation of 4.43 % for the studied dataset.Recently, Liu et al. (2022) presented a semi-empirical model based on the Modified Enskog Theory to describe the dynamic viscosity of liquid refrigerants, among which different HFOs and HCFOs, and their mixtures.The proposed equation is the following: where η 0 is the dilute gas viscosity in Pa s, ρ is the number density in m − 3 , b' is the covolume in m 3 and g(σ) is the radial distribution function.Instead, B corresponds to a fluid-specific empirical parameter that depends on temperature and pressure.In the original work, the Chapman-Enskog solution of the Boltzmann equation was used to calculate η 0 by considering that the molecular interactions can be roughly represented by the Lennard-Jones model.Instead, the Peng-Robinson EoS was used to calculate b ' and g(σ).As highlighted by the authors, the comparison with other literature models proved the high accuracy of their simple correlation for the selected refrigerants.In recent years, neural networks and machine learning have become increasingly popular in calculating the thermophysical properties of fluids.These techniques are particularly well suited for addressing complex problems with many interacting variables.For example, recent applications of these approaches have included the development of predictive models for several thermophysical properties, including surface tension (Gharagheizi et al., 2011;Pazuki et al., 2011;Pierantozzi et al., 2021), thermal conductivity (Eslamloueyan and Khademi, 2009;Pierantozzi and Petrucci, 2018) and volumetric properties (Taghizadehfard et al., 2019;Yousefi and Karimi, 2013).Different neural networks were proposed to model η of refrigerants and their mixtures (Di Nicola et al., 2022;Ghaderi et al., 2013;Wang et al., 2020;Zhi et al., 2018).Apart from the ANN proposed by Ghaderi et al. (2013), they can also calculate η for various refrigerants with low GWP.This work presents a general correlation specifically developed to describe the liquid dynamic viscosity dependence on temperature and pressure for low GWP refrigerants.One of the main features of this model is its greater simplicity compared to the literature correlations described above.In fact, it has a simple expression with only 3 input parameters (i.e., reduced temperature, reduced pressure, and acentric factor) and 6 coefficients that can be easily implemented with a programming language or even a spreadsheet.The coefficients were regressed on all the measurements collated from the literature.In addition, an ANN model having the same 3 input parameters has been proposed to predict the viscosity of the studied refrigerants.The results provided by the correlation and ANN were compared with those obtained from literature models.Finally, to ensure higher usability of the proposed correlation and ANN, an executable file that is easy to use and free for research and teaching is presented.

Data selection and analysis
Through a literature survey, 794 experimental data of η were selected for the following seven low GWP alkenes refrigerants: R1233zd (E), R1234yf, R1234ze(E), R1234ze(Z), R1224yd(Z), R1336mzz(E), and R1336mzz(Z).The selection of the η data was carried out by following a fluid-by-fluid and critical analysis.In particular, we discarded experimental points with deviations greater than three sigma with respect to the mean values, clearly beyond the common trend, and without details on the measurement conditions (i.e., pressures and densities).Measurements at T r > 0.9 were also neglected because simple estimation models, such as the ones studied here, generally cannot accurately describe the η behavior close to the critical point.
Table 1 provides the number of points and data references for the studied fluids, together with the temperature, pressure, and dynamic viscosity ranges.Table 2 summarizes the physical properties of the studied refrigerants and their references.Most of the properties reported in Table 2 were collated from the studies presenting the EoSs used in REFPROP 10.0 (Lemmon et al., 2018).However, the critical properties of R1224yd(Z) (Sakoda and Higashi, 2019) and R1336mzz(Z) (Alam et al., 2017) were selected from other reliable sources.Instead, since this refrigerant is not available in the list of pure fluids included in REFPROP 10.0, the R1336mzz(E)'s properties reported by Tanaka et al. (2017) Table 1 Liquid dynamic viscosity measurements for the selected low GWP refrigerants together with the relative combined (u(η)) and expanded (U(η)) uncertainties reported in the original works.were used in this study.
Fig. 1 shows the behaviors of the η measurements selected for the seven analyzed refrigerants as a function of T r and p r .It is evident from Fig. 1 that their liquid viscosity increases as the temperature decreases and the pressure increases.However, the effect of pressure on η is limited up to moderate pressures; in fact, the values of η are almost constant at up to p r = 1 and fixed T r .

Models for the liquid viscosity of low GWP refrigerants
In this section, details about the proposed correlation and ANN for the liquid dynamic viscosity of the low GWP refrigerants are reported.It is worth pointing out that an executable file for calculating η of the studied fluids with both the proposed correlation and ANN is provided in the supplementary material, together with detailed instructions for its use.The executable file is free for research and teaching.

Proposed correlation
Starting from the correlation proposed by Latini and his co-authors (Eq.( 1)) and its modified version (Eq.( 3)), a new correlation that combined simplicity and reliability was designed for the η of the studied fluids.The development of the model was based on the idea that it should have a limited number of coefficients and input physical properties well-known for the refrigerants under study.To this end, we analyzed various equations developed by simplifying Eq. ( 3) and characterized by only three input properties and six coefficients.In addition to T r and p r , which allow to describe the temperature and pressure dependence of the liquid dynamic viscosity, the other properties shown in Table 2 were tested.It was found that the correlations with the acentric factor provided lower deviations than the others; in particular, the following equation ensured the best results: where η is in mPa s, ω is the acentric factor, and A, B, C, D, E, and F are regressed coefficients.
The greater simplicity of Eq. ( 5) compared to Eq. ( 3) is proved by the fact that it involves only three input properties.Although the acentric factor was not considered in the models developed by Latini and his coauthors, this property was used in the correlations to calculate the viscosity of many fluids.In this regard, the acentric factor appears in the model proposed by Lucas (1981) to describe the effect of pressure on liquid viscosity.It was also used to calculate the intrinsic molal volume in the model for the viscosity of liquids proposed by Przedziecki and Sridhar (1985).Moreover, it was adopted to obtain the saturated liquid density used in the correlation proposed by Dutt (1990) to estimate the kinematic viscosity of petroleum crude oil fractions.
To take into account the uncertainty of the experimental data in the regression of the coefficients, the values of A, B, C, D, E, and F of Eq. ( 5) of the seven refrigerants were obtained by the minimization of the following objective function OF for the complete dataset (794 data points): where η exp is the experimental liquid dynamic viscosity, η calc is the calculated liquid dynamic viscosity, U(η) is the expanded uncertainty, and N is the number of experimental data.The values of U(η) reported in Table 1 were used in the regression.The Levenberg-Marquardt optimization algorithm (Moré, 2006) was used for regression.Moreover, fluid-specific coefficients were regressed for each refrigerant by minimizing OF.Table 3 presents the obtained values.

Artificial neural network
Artificial neural networks (ANNs) are mathematical models considered a key component of machine learning and artificial intelligence because they can recognize complex patterns and adapt their behavior on the basis of new information (Hochreiter and Schmidhuber, 1997).In particular, they can predict a given dataset using a series of mathematical relationships known as artificial neurons contained in "layers" (Rumelhart et al., 1986).Once the model has been trained to identify patterns within the data used for its development, the ANN can predict the outcome for a completely new dataset.In fact, it can process multiple inputs through several hidden layers to finally arrive at an output layer.This method is particularly useful in cases where a large amount of data is available or where performing analysis is made difficult by the lack of data (Maren et al., 2014).
The steps involved in building a neural network include collecting the relevant data and then developing the model using training, validation, and test datasets (Faúndez et al., 2020).The training dataset allows the calculation of the network's parameters so that it can predict new values as best as possible.The validation set is used to evaluate the network's internal performance and modify the hyperparameters during the tuning process.The test set allows the evaluation of the accuracy of the network's tuning process.In fact, this dataset is not part of the training process, so the network is unrelated to this data.In selecting the best ANN configuration, it is crucial to analyze the results provided by the test set during the training process.
The neurons' and layers' numbers are two of the most important elements in choosing the right architecture.They determine the complexity of the model and the amount of information the network can process.It is essential to select a proper number of neurons that ensures accurate results without having a network that is too sensitive to data.In fact, having a network with too many neurons may improve its overall performance, but then expose the network to the problem of overfitting (Oyedotun et al., 2017).In this case, the network fits the training data very well, but it cannot describe the underlying pattern of the data, so it does not make reliable predictions on unseen data.The layers' number determines how the inputs are connected to the output nodes.Despite multiple hidden layers can improve the ANN accuracy, it was shown that an increase in their number can enhance the network's susceptibility to convergence on suboptimal local minima, which can detrimentally affect the generalization capabilities of the model (Karsoliya, 2012;Liu et al., 2007).Another key element is the activation function which determines the level of interaction between the neurons in the network.This function is used to define the relationship between the input and output variables of a network.In other words, it decides how the output node responds to a particular input.Some common activation functions used in ANNs are the logistic sigmoid, rectified linear unit, hyperbolic tangent, and tanh.
Once the aspects described above have been defined, the next step is to build and train the model.The ANN training usually involves feeding the neural network with input data, observing how the system responds, and making adjustments to ensure the output is as close as possible to the desired value.This process is repeated several times until the output produced by the network matches the desired result.
To ensure a good trade-off between computational complexity and model performance, an ANN with only one hidden layer, three input parameters, and one output parameter (i.e., the liquid dynamic viscosity) was chosen as the network configuration for calculating η for the studied liquid refrigerants.The same input properties of Eq. ( 5) were selected as input parameters: the reduced temperature (T r ), the reduced pressure (p r ), and the acentric factor (ω). Since viscosity is a non-linear phenomenon with respect to temperature, the sigmoid was chosen as the activation function, which can guarantee an accurate response in nonlinear problems.To have low deviations without the problem of overfitting, it was decided to have a network with one hidden layer and the fewer number of neurons that could guarantee good performance.As can be seen from Fig. 2, the ANN with 19 neurons in the hidden layer provided the best performance on all the datasets during the training process.The dataset was divided using random allocation into three subsets for training, validation and testing purposes.The percentage of instances in each subset is 80 % for training, 10 % for validation, and 10 % for testing.
Therefore, as shown in Fig. 3, the ANN with a hidden layer with 19 neurons, three input parameters (T r , p r , ω), and one output layer (η) was selected as the network architecture to calculate η for the studied fluids.

Results and discussions
Table 3 shows the results obtained with the proposed correlation (Eq.( 5)) using the coefficients regressed for the complete dataset and for each refrigerant.The coefficients obtained from all the selected data gave an average absolute relative deviation of η (AARD(η)) equal to 2.88 %.Instead, as expected, the fluid-specific coefficients ensured more accurate results, yielding AARD(η) always lower than 2 %.The AARD (η) is defined as: where η exp is the experimental liquid dynamic viscosity, η calc is the calculated liquid dynamic viscosity, and N is the number of experimental data.
In general, the results reported in Table 3 are satisfactory and show uniformity in terms of coefficients and deviations between the various fluids under investigation, even considering the different number of experimental data points collated for each fluid.Fig. 4 shows the consistency between η calculated from Eq. ( 5) with the coefficients obtained for the complete dataset and the experimental data for the studied refrigerants.
Going into more detail with the result analysis, the AARD(η) and maximum absolute relative deviations (MARD(η)) for each fluid provided by the proposed correlation (with the coefficients regressed for the overall dataset) are reported in Table 4.This table also presents the results of Eq. (3), Eq. ( 4), the proposed ANN, and REFPROP 10.0.For Eq.
(3), the coefficients for the alternative refrigerants proposed by Di Nicola et al. (2022) were used in the calculation.For Eq. ( 4), the fluid-specific B coefficients available in the original work (Liu et al., 2022) were considered here.For R1336mzz(E), it was impossible to perform calculations for Eq. ( 4) because of the lack of its B coefficients.Moreover, Table 4 does not show the REFPROP 10.0 results for R1336mzz(E) because this fluid is not yet available on the software routines.Therefore, it is worth noting that, unlike the proposed correlation, REFPROP 10.0 and Eq. ( 4) do not allow calculating η of all the studied refrigerants.
Comparing Table 3 and Table 4, it is worth noting that, even when the complete dataset is regressed, deviations for Eq. ( 5) are similar to the ones obtained when the fluids were regressed individually, with the only exception of R1234ze(Z).This outcome confirms the validity of the model.
Analyzing in detail the results achieved with the presented correlation, Figs. 5 and 6 show deviations between the selected measurements and the calculations obtained from Eq. ( 5) using the coefficients obtained from the complete dataset versus reduced temperatures and pressures, respectively.From these figures, higher deviations for R1234ze(Z) appear for the proposed equation, both when the analysis is performed regarding reduced temperatures and reduced pressures.Generally, experimental data for R1234ze(Z) were obtained at low  reduced pressures and high reduced temperatures.However, excluding R1234ze(Z), all fluids' deviations are well within 10 %.
Comparing the different correlation performances, however, it is evident that the proposed model (Eq.( 5)) gives, on average, better results than the ones already present in the open literature (Eqs.( 3) and ( 4)).On the other hand, REFPROP 10.0 provided an accurate viscosity description for the refrigerants available in the software.This result was expected since it uses models derived from their experimental data.
Regarding the results obtained with the ANN, Table 4 shows that this model is undoubtedly the most accurate and can therefore be considered a useful reference for the other methods and even a lower limit for the deviations provided by a model.In fact, its results are promising, especially considering the dataset size and the limited amount of input parameters used to train the model.The relatively low values of MARD (η) for all the fluids prove that the network works well on the data used for training and the data used to test its prediction capability.Furthermore, by comparing the results of the ANN with that reported in Table 3, it can be seen that the network ensured even better results than those obtained from Eq. ( 5) with fluid-specific coefficients.this study, an empirical correlation and an ANN were proposed to represent the liquid dynamic viscosity dependence on temperature and pressure for low GWP refrigerants.Both models have only three input parameters (i.e., reduce temperature, reduce pressure, and acentric factor) and were developed using 794 experimental data collated from the literature for 7 refrigerants.
An accurate description of the dynamic viscosity for the studied liquid refrigerants was ensured by the proposed correlation using the coefficients regressed from all the selected experimental data (AARD (η) of 2.88 % for the complete dataset).In this regard, the proposed model provided lower deviations than other literature correlations, such as a modified version of the correlation proposed by Latini and his coauthors (Latini et al., 2002(Latini et al., , 1996(Latini et al., , 1990) (AARD (η) = 4.43 %) and the correlation proposed by Liu et al. (2022) (AARD (η) = 3.30 %).Instead, REFPROP 10.0 gave better results than the correlation (AARD (η) = 1.46 %).However, REFPROP 10.0 did not allow calculating η for R1336mzz (E) because no models for this refrigerant are available on the software routines.Therefore, the proposed correlation can be considered a reliable tool for calculating the liquid dynamic viscosity of the analyzed low GWP refrigerants, especially for the less-known refrigerants, such as R1336mzz(E).In fact, unlike other literature models, it ensures a good trade-off between accuracy and simplicity for these refrigerants.
From the training process, it was found that a single-layer ANN with  5) using the coefficients regressed for the complete dataset versus the experimental data for the studied refrigerants.

Table 4
AARD(η)% and MARD(η)% given by Eq. (3), Eq. ( 4), Eq. ( 5), the artificial neural network (ANN), and REFPROP 10.0 for each fluid and the complete dataset.19 neurons provided the best results.In particular, this network architecture gave the lowest deviations for all the studied fluids with respect to the other studied models (AARD (η) = 0.86 % for the whole dataset).This outcome proved that the ANN is undoubtedly the most accurate model and can be considered a good reference for developing other models.Finally, it is worth remarking that, to allow researchers and scholars to use the models presented in this work, an executable file for calculating the dynamic viscosity of the studied low GWP refrigerants with both the proposed correlation and ANN is provided in the supplementary material.

Fig. 1 .
Fig. 1. 3D behaviors of the liquid dynamic viscosity experimental points (η) for the low GWP refrigerants as a function of the reduced temperature (T r ) and reduced pressure (p r ) (a) and their projections to the η -T r plane (b) and η -p r plane (c).The filled and open symbols refer to the measurements at p r < 1 and p r ≥1, respectively.

Fig. 2 .
Fig. 2. AARD(η) for the complete, training, validation, and test datasets as a function of the neurons' number.

Fig. 3 .
Fig. 3. Schematic diagram of the proposed ANN model.T r is the reduced temperature, p r is the reduced pressure, ω is the acentric factor, and η (mPa s) is the liquid dynamic viscosity.

Fig. 4 .
Fig. 4. Liquid dynamic viscosities (η) obtained from Eq. (5) using the coefficients regressed for the complete dataset versus the experimental data for the studied refrigerants.

Fig. 5 .
Fig. 5. Deviations between the experimental data and calculations of Eq. (5) as a function of the reduced temperature.

Fig. 6 .
Fig.6.Deviations between the experimental data and calculations of Eq. (5) as a function of the reduced pressure.

Table 2
Physical properties of the selected refrigerants.