Machine learning analysis of thermophysical and thermohydraulic properties in ethylene glycol- and glycerol-based SiO2 nanofluids

NF preparation

SiO₂ NPs powder was purchased from US Research Nanomaterials Inc., USA. EG and G with 99% high purity were obtained from R&M Chemicals. The scanning electron microscopy (SEM) image of the SiO₂ was taken using a Field Emission Scanning Electron Microscope (Supra 55VP FE-SEM, Carl Zeiss). Figure 1a depicts the SEM images. As seen, the NPs are nearly spherical shaped. The particle sizes vary between 8 and 25 nm, with an average particle size estimated as 21 nm²⁹. The EDX analysis of the selected area in the SEM image is displayed in Fig. 1b. The result indicates that the composition of SiO₂ is \(100\%\) consistent with vendor specifications.

NFs are prepared in the respective base liquids, G and EG, following a two-step method¹⁴. Each quantity of NPs required for 0.5%, 0.75%, 1.0%, and 2.0% volume fractions was taken on an electronic balance (TLE 104E, Mettler-Toledo). The formulation was achieved by dispersing SiO₂ in the base liquids with a magnetic agitator to stir the sample in a beaker for 30 min while adjusting the pH value. Following a similar protocol, another set of SiO₂ NFs in a 60:40 by volume G and EG was prepared for comparison.

It is reported that NPs would agglomerate to form clusters and settle over time due to high surface energies³⁰. Kinetic energy is needed to break down the particle clusters into minute sizes, according to Darzi et al.¹⁵. Given that, the samples have been subjected to an ultrasonic homogenization (Labsonic M, Sartorius AG) operating at a 30 kHz frequency for 2 h to improve the stability of the dispersions. The pH control is crucial in colloidal stability as it determines the suspension’s isoelectric point (IEP). Highly acidic NF (low pH) can lead to corrosion with the long-term flow in pipes. The pH value of the NF samples varied between 6 and 7. The pH value was set to 10 by adding NH₃OH to have zeta potential values away from IEP³¹. No dispersant was added during the preparation process to avoid altering the properties of the NF.

The stability of NFs was checked using Zetasizer (Nano ZSP, Malvern), which operates on the principle of dynamic light scattering to measure charge repulsion/attraction between dispersed particles. At pH 10, the dispersions had an average absolute zeta potential of − 33 and − 42 mV indicating that the SiO₂ NFs are stable. Further, a small portion of the NFs is kept under static conditions for months and examined. G-based NFs were stable for over three months, whereas EG-based NFs were stable for up to 1 month without settling in clear storage containers.

Thermophysical property measurement

The thermal property analyzer (KD2 Pro Decagon Devices) based on the transient hot-wire method was used to evaluate the effective TC of NFs within a specified precision of 5.0%. The KS-1 probe with a 60 mm length and 1.28 mm diameter was selected, which provides a transient line heat source. The sample temperature was controlled using an isothermal bath (Vivo-RT2, Julabo), with temperature stabilization better than ± 0.1 K. A rotational rheometer (MCR 302, Anton Paar) was used for effective viscosity measurement. A double-gap concentric cylinder was employed as the measuring geometry. A gap distance of 1.0 mm was allowed between the co-axial cylinders of the system. A Peltier thermostat controlled the cell temperature with a precision of ± 0.1 K. Repeated tests were conducted LVDV-III Ultra Programmable Rheometer. A digital densitometer (DA-645, KEM) which functions on an oscillating U-tube principle, was instrumental in determining the effective density with ± 0.00005 g/cm³ accuracy. Peltier thermoelectric elements enable temperature control within the measuring cell, assuring a precision below 0.03 °C. A differential scanning calorimeter (DSC Q2000, TA Instruments) analyzed the specific heat with an accuracy of 2%. Precise heat measurement was made using the standard test method (ASTM-E1269) under a high-purity nitrogen atmosphere at a 20 °C/min heating rate in the DSC furnace. The device temperature is ± 0.01 °C. A refrigeration cooling system RCS90 was used to conduct specific heat testing at different temperatures.

All devices have been calibrated with either G/EG before the measurements with NFs. The data are collected in the range of temperatures 20–80 °C at atmospheric pressure. Three readings were obtained for each sample at each temperature, and the mean value was stated.

Table 1 shows the measured properties of G and EG. Comparisons have been made with the values reported in the literature^32,33,34,35. The Hewitt³³ data correlates well with thermal conductivity values within a maximum of 0.25 and 0.02% deviations for glycerol and ethylene glycol, respectively. Hewitt and Cabaleiro’s specific heat data showed a 0.9 and 5.6% variation, while the Lide thermal conductivity data deviated from the measured values by 1.3 and 0.6%. A maximum deviation in viscosity of 5.9 and 1.1% for G and EG was observed when compared with Lide³⁵ and Quijada-Maldonado³², respectively, in the temperature range of 25 to 80 °C. The overall deviation of the calibration results in all experiments was better than 6% from the reference values. Table 2 presents the uncertainty of the measured thermophysical properties of glycerol (G) and ethylene glycol (EG) in percentage.

A comparison of the thermophysical properties of SiO₂ NFs with the base liquid is shown in Fig. 2a–d. It was established that the TC enhancement was highest at 1.0% concentration with values of 4.2% and 10.7%, respectively, for G and EG-based NFs. The NF exhibited Newtonian behavior with viscosity independent of the shear rate in a similar manner to that evidenced by Tadjorodi et al.³⁶ and Żyła and Jacek³⁷. The viscosity of SiO₂ NFs increased by 27% and 33% compared to base liquids G and EG. The suspension density increased by nearly 2% approximately. SiO₂ NFs exhibited lower values of effective specific heat than their base fluids. Meanwhile, the specific heat decreased by nearly 2.7% and 1.5%, over a temperature range of 25 to 80 °C. Further, the measured density values of NFs are consistent with the calculated values using the mixing theory relation within 1.0% deviation in each case. The specific heat data deviated by 3% from the classical thermal equilibrium model. Escher et al.³⁸ have reported a good agreement for density and specific heat values with mixture relations for water-based SiO₂ NFs for a concentration of up to 31%.

The experimental data for the measured SiO₂ NFs thermal conductivity, viscosity, density and specific heat was regressed to generate an empirical Eqs. (1)–(16) as a function of particle volume concertation (φ) and temperature (T) for further analysis of heat transfer as follows:

SiO₂-EG nanofluid

SiO₂-G nanofluid

With the parameter units \(\varphi \) (%) and \(T\) (^oC).

Heat transfer test setup

The experimental system has a closed-loop design composed of five essential components: the test tube, power supply, cooling arrangement, measurement system, and data acquisition. The experimental setup schematic diagram is depicted in Fig. 3. The test section involves a single stainless-steel tube with a bell-mouth entry, thermocouples, and split heaters. The test section consists of a tube 1.0 m long with outer and inner diameters of 6.35 and 4.57 mm, respectively. The outer tube surface is embedded with two-cylinder split-body electric heaters. Each heater, with a 750 W maximum power rating, is wrapped with ceramic fiber insulation and connected to a variable transformer. Eight K-type thermocouples have been used to quantify and record the temperatures at different places. Two of the thermocouples are spot welded to the test section, each at an axial distance of 150 mm from either end of the tube, to measure the wall temperature; four thermocouples at equidistant positions of 100 mm from the tube end are located to measure the surface temperatures and two thermocouples at the inlet and outlet of the test section. All thermocouples were calibrated within a normal of ± 0.1 °C.

The cooling unit comprises a chiller, circulating pump, water tank, and temperature control system. The chiller, rated at 0.74 kW, was connected to the plate-type heat exchanger to regulate the fluid temperature at the test section inlet. A 3.0 hp horizontal multistage pump (AB, Teral) circulated the fluid in the test section. A digital vortex flow meter (SV4200, IMF) with a range of 1–20 LPM is utilized to quantify the fluid flow. A volumetric cylinder made of acrylic of 1.0 L capacity with scale graduation was connected at the test-section exit as a working fluid reservoir to check the flow rate of the fluid visually. Two absolute pressure transducers (GP-100, Keyence) were installed at the pipe inlet and outlet to measure the pressure drop (∆P) across the test section. The calibration range of the sensors from 0 to 10 MPa is ± 1%. All measuring instruments in the circuit were connected to the data logger for recording output signals.

The experiments were undertaken between 6 and 12 LPM flow rates corresponding to flow Reynold numbers of 1200–22,000, while the working fluid temperature is maintained at 80 °C. All the readings are recorded at steady-state conditions. The circuit is cleaned with water and air-dried between successive experiments.

Data analysis

The heat energy \({Q}_{i}\) provided to the working fluid in the heat section is a function of electric current \(I\) and voltage:

Simultaneously, the rate of heat transfer \({Q}_{a}\) was evaluated from the mass flow rate \(\dot{m}\) and the fluid temperatures at the inlet and outlet of the tube:

Under a steady state, the energy available with the hot fluid exiting the test section should equal the heat removed by the cooling liquid in the chiller. Newton’s law of cooling relation is used to evaluate the convective HTC as follows:

where the A, \({T}_{w}\) and \({T}_{b}\) are surface area, wall temperature and the fluid bulk or average temperature, respectively computed as:

The average Nusselt number was estimated based on the convective HTC \(h\), tube diameter \(D\), and TC of the fluid \(k\) as

Furthermore, the turbulent HTC can be considered from Dittus–Boelter equation in the following form:

where c, m, and n represent the coefficients suitable for NF experimental data. The \(Re\), and Prandtl number \((Pr)\) terms are defined as follows:

The NF properties used for heat transfer analysis are determined at the bulk temperature. Further details on the derived equations for the thermophysical properties of the NF are given in Ref.²⁹. The nondimensional FF \(f\) was calculated from Darcy–Weisbach equation, which relates the ∆P, pipe length \(L\), hydraulic diameter \(D\), fluid density \(\rho \), and average velocity \(u\) as follows:

Uncertainty analysis

Analysis of the experimental uncertainty was undertaken to validate the precision of the measurements. The uncertainties in heat transfer characteristics were estimated based on the error approach presented by Beckwith et al.³⁹, following the protocol described in^17,40. The instrument and uncertainties estimated from the measured parameters are presented in Table 2.

Theoretical correlations for Heat transfer

Model 1

A correlation for Nusselt number of single-phase liquids under fully developed transition and turbulent flows is given by Gnielinski⁴¹ as:

where the fanning friction factor is given as \({f}_{F}={\left(1.58\text{ln}Re-3.28\right)}^{-2}\). Equation (25) is valid in \(the 2300\ge Re\ge {10}^{6}\) range.

Model 2

Other correlations are considered for a developing flow in a circular tube with a small velocity boundary layer thickness. Del Giudice⁴² developed a model for the developing flow heat transfer in a pipe exposed to uniform wall heat flux with the consideration of temperature dependence of viscosity and thermal conductivity as:

where, \({X}^{*}=\frac{L}{{D}_{h}RePr}\) \(n=0.761{\left(RePr\right)}^{0.0224}-0.000109 RePr\)

Where Pn_µ is viscosity Pearson number = (β \({q}_{w}^{“}\) D/k_e); β =  − (dµ/dt)/µ; \({q}_{w}^{“}\) is the heat flux at the tube surface W/m²; D is the tube inner diameter; k_e is the TC at tube entry temperature. Equation (26a) is valid for \(5.0 \le \text{ Pr }\le 100;\) 10^–4 ≤ X^* ≤ X^*_max. The value of \({X}_{max}^{*}\frac{1}{ 4{Pn}_{\mu }}ln\left[\frac{{\left({Re}_{b}\right)}_{max}}{{Re}_{e}}\right]\) is estimated for the experimental conditions to be 0.08.

Model 3

Muzychka and Yovanovich⁴³ presented a model for predicting Nu in the combined entrance region of a tube valid for uniform wall flux boundary conditions given by

where, \(f\left({R}_{e}\right)=\frac{12}{ \sqrt{\epsilon \left(1+\epsilon \right)} \left[1-\frac{192\epsilon }{{\pi }^{5}}tanh\left(\frac{\pi }{2\epsilon }\right)\right]}\); \(f\left({P}_{r}\right)=\frac{0.866}{ \sqrt{\epsilon \left(1+\epsilon \right)} \left[1-\frac{192\epsilon }{{\pi }^{5}}tanh\left(\frac{\pi }{2\epsilon }\right)\right]}\); \(m=2.27+1.65{P}_{r}^{1/3}\); \({z}^{*}={L/{D}_{h}}{R}_{e}{P}_{r}\). For a circular tube, \(\epsilon =1, {C}_{1}=3.86, {C}_{2}=1.5, {C}_{3}=0.501,\) and \({C}_{4}=2\)

Equation (26b) is valid for 0 < Z^* < ∞ and 0.01 < Pr < ∞.

Machine learning

The experimental data collected in the last section was employed to develop a comprehensive set of models to prognosticate the thermohydraulic behavior of Ethylene Glycol and Glycerol based non-porous SiO₂ nanofluids. A battery of Python-based open-source libraries was used in the Jupytor environment. A total of five ML techniques were employed for the development of prediction-models in this case. The LR, RF, DT, XGBoost, and AdaBoost models were chosen for their various strengths: LR for simplicity and interpretability, RF, and XGBoost for excellent prediction accuracy, DT for simple decision-making, and AdaBoost for improving weak learners.

The LR was used to prepare the baseline model while the other four namely Random Forest (RF), Extreme gradient boosting (XGBoost), Adaptive boosting (AdaBoost), and Decision tree (DT) were compared against it. A brief description for each is provided as follows:

Linear regression

Linear regression (LR) is considered most basic form of supervised ML algorithm. In this case, a linear equation is fitted to the actual data for modeling the correlation with independent variables (features) and a dependent variable (target). The objective in this case is to represent their mutual relationship. It can be expressed as follows:

Herein, y is the target (dependent variable), \({x}_{1}, {x}_{2}, {x}_{3}\dots .\) are features (independent variable), \({\beta }_{0}, {\beta }_{1}, {\beta }_{2}, \dots .\) are coefficients and \(\varepsilon \) denotes the error.

LR algorithm is employed to locate the line which may provide the best fit to that minimize the sum of squared differences between the actual and predicted values.

Random Forest

Random Forest (RF) is a type of ensemble learning system. It is developed for regression type complex problems. In the training phase, RF generates a large number of decision trees, each of these is trained on a random portion of the training data and characteristics.

Let we denote, a training data set denoted as ‘X’ having n samples and m features. In this case y is the target variable, T denotes total count of decision trees in the test forest, X_i denotes a random subset taken out from training data X, sampled with replacement in case when i ranges from 1 to T. Similarly in case of features F_j denotes the random subset wherein, the j ranges from 1 to T.

In case of each decision tree i, a sample for training denotes as (X_i, y_i) is randomly selected from (X,y). Then a decision tree Di is trained on (X_i, y_i) employing a split criterion on the basis of least MSE.

The final prediction using a RF-based model is expressed as:

In this case \(\widehat{y}\) is the forecasted output D_i(x) denotes the forecast from i-th decision tree for input x.

To summarize, RF regression employs the aggregate of several DTs trained on random subsets of data to generate robust and precise forecasts for regression problems.

Decision Tree

Decision Tree (DT) is a fundamental and flexible ML method which can be used for regression of data. DT creates a hierarchical tree framework having core nodes denoting feature-based decisions while leaf nodes represent the predicted values.

Let us denote a training data set as ‘X’ with n samples and m features. In this scenario, y is the goal variable, D is the DT model.

In the training phase, the DT splits the feature space recursively into subgroups on the basis of feature values. At each node, the DT selects the feature and split threshold to keep the MSE as low as possible. In the regression work, the prediction is done by traversing the tree from root to leaf node and assigning the mean value of the target variable inside the leaf node to the input sample.

Mathematically, the DT-based forecast can be expressed in simple terms as:

It can be summarized that DT-based regression separates the feature space recursively and then makes forecast using the mean target variable value inside each zone, resulting in interpretable and simple regression models.

Extreme gradient boosting

Extreme gradient boosting (XGBoost) belongs to the gradient boosting family of ML. XGBoost is known for its exception prediction performance. In the gradient boosting process, it sequentially combines weak learners so as to form a robust prognostic model. XGBoost make use gradient descent approach for minimizing the loss function ‘L’ such that

In this case, the \(\theta \) represents the parameters of model, whereas \({y}_{i} and \widehat{{y}_{i}}\) denotes the actual and predicted values. K denotes the number of weak learners in form of trees, and \(\Omega ({\text{f}}_{k})\) denotes the regularization term applied to each tree.

XGBoost uses the additive approach for model building as:

Herein, \({\text{f}}_{k}({\text{x}}_{i})\) represents the prediction of the k-th tree in case of i-th sample.

Regularization techniques like L1 and L2 are used for controlling the complexity of the individual trees.

Herein, T is number of leaves in tree, \(\omega \) denotes the leaf weights, and \(\gamma \) and \(\lambda \) are regularization parameters.

XGBoost also provides insights into feature importance by calculating the gain, which assesses the importance of each feature to the model’s performance.

To summarize, XGBoost improves the model’s accuracy by successively minimizing a loss function with gradient descent and incorporating regularization approaches to control model complexity. Its ability to offer feature importance and handle missing values render it an effective tool for regression problems across multiple domains.

Adaptive boosting

Adaptive boosting (AdaBoost) combines multiple weak leaners h(x) to form a strong predictor f(x). In mathematical terms the final predictor F(x) is a weighted sum of learners of weaker learners:

Herein, T is the total number of weal learners, \({\alpha }_{1}\) is the weight allotted to the t-th weak learner, and \({h}_{1}(x)\) is the prediction of t-th weak learner.

In the training phase, AdaBoost allots weight to each instance of training \(({x}_{i}, {y}_{i})\), in case when \({x}_{i}\) denotes the input while \({y}_{i}\) denotes true label. In starting phase all weights are set equally as:

\({{w}_{i}}^{(1)}= \frac{1}{N}\);where N is the total count of instances.

AdaBoost fits a weak learner to the training data \({w}_{i}\) in each iteration t, it subsequently computes the weighted error \({\varepsilon }_{t}\) of each weak learner.