Demand forecast and influential factors of cold chain logistics based on a grey model

1.   Introduction

The smoothed particle hydrodynamics (SPH) method is a Lagrangian representation that is based on a particle interpolation to calculate smooth field variables ^[1,2]. Because of its Lagrangian features, the SPH method has fine advantages. The SPH method is used to simulate several industrial applications by simulating highly dynamic and violent flows such as tank sloshing, oil flow in gearboxes, nozzles, and jet impingement. In addition, the SPH method works well with the free surface flows, surface tension in free surface flows, and breaking waves without a special treatment required in this area. Additionally, the simulations of multiphase flows are easily handled by the SPH method. In conventional grid-based approaches, both the Volume-of-Fluid (VOF) method and the Level set method are adept at accurately capturing multiphase interfaces. An additional strength of the SPH method is its ability to handle complex/moving geometries. He et al. ^[3] introduced a new, weakly-compressible, SPH model tailored for multi-phase flows with notable differences in the density, while also allowing for large Courant-Friedrichs-Lewy numbers. Their method involves adjusting the continuity equation to exclude contributions from particles of different phases, solely focusing on neighboring particles of the same phase. Subsequently, they utilize a Shepard interpolation to reset the pressure and density of the particles from other phases. Bao et al. ^[4] presented an upgraded method within the SPH aimed at simulating turbulent flows near walls with medium to high Reynolds numbers. Their approach combined the $\kappa -\varepsilon$ turbulence model and a wall function to enhance the accuracy. To tackle numerical errors and challenges related to the boundary conditions, they split the second-order partial derivative term of the composite function containing the turbulent viscosity coefficient into two components. Ensuring uniform particles, smooth pressure fields, and efficient computation, they employed particle shifting techniques, the $\delta$ -SPH method, and graphics processing unit (GPU) parallel processing in their simulations. Huang et al. ^[5] devised an innovative approach that merged SPH with the Finite Difference Method (SPH-FDM) within an Eulerian framework to analyze the unsteady flow around a pitching airfoil. They introduced an improved Iterative Shifting Particle Technology (IPST) to accurately track the motion of pitching airfoils. The numerical results from the SPH-FDM approach within a Eulerian framework closely aligned with those from either the finite volume method (FVM) or the existing literature, affirming the effectiveness of their proposed method. Further applications of the SPH method and its incompressible scheme, entitled the ISPH method, can be found in the references ^{[6,7,8,9,10,11,12,13,14,15,16,17]}. In recent years, the ISPH-FVM coupling method refers to the integration of the ISPH method with the FVM. This hybrid approach combines the strengths of both techniques to accurately model the fluid flow phenomena, particularly in scenarios involving either complex geometries or multiphase flows. Xu et al. ^[18] introduced the ISPH-FVM coupling method to simulate a two-phase incompressible flow, and combined the advantages of ISPH in interface tracking and FVM in flow field calculations with enhancements in surface tension calculations inspired by VOF and Level Set methods, demonstrating an improved accuracy and stability in complex two-phase flow scenarios. Xu et al. ^[19] developed a 3D ISPH-FVM coupling technique to model a bubble ascent in a viscous stagnant liquid, and employed an exchange of information at overlapping regions between the ISPH and FVM schemes, with an evaluation of surface tension effects using the 3D continuum surface force model (CSF). They validated the accuracy of the interface capture through benchmark tests and compared simulations of a single bubble ascent against experimental and literature-based data; moreover, they explored the method's capabilities in various parameter settings and in studying the dynamics of a two-bubble coalescence, demonstrating the reliability and versatility of the approach in simulating the bubble ascent phenomena. Xu et al. ^[20] proposed an advanced mapping interpolation ISPH-FVM coupling technique tailored to simulate complex two-phase flows characterized by intricate interfaces. The nanofluids were a suspension of nanomaterials in the base fluids. The nanofluids were applied to enhance the thermal features of cooling devices and equipment in electronic cooling systems ^[21,22] and heat exchangers ^[23,24]. In 1995, Choi and Eastman ^[25] achieved an enhanced heat transfer by adding nanoparticles inside a base fluid. The double diffusion of nanofluids (Al₂O₃-H₂O) in an enclosure under the magnetic field impacts was studied by the lattice Boltzmann method ^[26]. There are several uses of coupling magnetic field impacts on nanofluids in thermal engineering experiments ^{[27,28,29,30,31,32]}. The variable magnetic field on ferrofluid flow in a cavity was studied numerically by the ISPH method ^[33]. Aly ^[16] adopted the ISPH method to examine the effects of Soret-Dufour numbers and magnetic fields on the double diffusion from grooves within a nanofluid-filled cavity. To enhance the heat transfer in cavities, the wavy or extended surface or fins were applied. The fins were used in several mechanical engineering applications. Saeid ^[34] studied the existence of different fin shapes on a heating cavity. The numerical attempt based on the finite volume method was introduced for natural convection and thermal radiation in a cavity containing N thin fins ^[35]. Hussain et al. ^[36] examined the effects of fins and inclined magnetic field on mixed convection of a nanofluid inside single/double lid-driven cavities. The recent numerical attempts based on the ISPH method for double diffusion of a nanofluid in a finned cavity, which were introduced by Aly et al. ^[37,38]. This work examines the impacts of Cattaneo-Christov heat flux and magnetic impacts on the double diffusion of a nanofluid-filled cavity. The cavity contains a rotated wavy circular cylinder and four fins fixed on its borders. The ISPH method is utilized to handle the interactions between the rotational speed of an inner wavy cylinder with a nanofluid flow. The established ANN model can accurately predict the values of $\stackrel{-}{Nu}$ and $\stackrel{-}{Sh}$, as demonstrated in this study, in which there is a high agreement between the prediction values generated by the ANN model and the target values.

2.   Mathematical analysis

The initial physical model and its particle generation are represented in Figure 1. The physical model contains a square cavity filled with a nanofluid and containing four fins and a wavy circular cylinder. The embedded wavy circular cylinder rotates around its center by a rotational velocity $\mathrm{V} = \omega \left(r-{r}_{0}\right)$ with adiabatic thermal/solutal conditions. The four fins are fixed at the center of the cavity borders with a low temperature/concentration ${T}_{c}\&{C}_{c}$. The source is located on the left wall of a cavity with ${T}_{h}\&{C}_{h}$, the plane walls are adiabatic, and the right wall is maintained at ${T}_{c}\&{C}_{c}$. It is assumed that the nanofluid flow is laminar, unsteady, and incompressible. Figure 2 depicts the architecture and fundamental configuration structure of the created multilayer perceptron (MLP) network. A model with 14 neurons in the hidden layer of the MLP network model was identified so that the number of neurons were included in the hidden layer ^[39,40]. To build ANN models, the data-collecting process must be idealistically structured. Three methods were used to arrange the 1, 020, 888 data points that were used to construct the ANN model. In the created MLP network model, the $\stackrel{-}{Nu}$ and $\stackrel{-}{Sh}$ values were estimated in the output layer, while the $\tau$ and $Sr$ values were defined as input parameters in the input layer. A total of 70% of the data was set aside for model training, 15% for validation, and 15% for testing. The recommended training algorithm in the network model was the Levenberg-Marquardt (LM) algorithm. The LM algorithm effectively minimizes nonlinear least squares problems by iteratively adjusting parameters based on an objective function, utilizing the gradient descent and Gauss-Newton methods, while dynamically controlling the step size to prevent overshooting; this achieves convergence to the optimal parameter values, making it widely used in various fields for an efficient optimization. Below are the transfer functions that are employed in the MLP network's hidden and output layers:

The mathematical equations utilized to compute the mean squared error (MSE), coefficient of determination (R), and margin of deviation (MoD) values, which were selected as the performance criteria, are outlined below ^[41]:

Figure 1.  Primary model of rotated wavy cylinder and four fins in the cavity's center borders.

DownLoad: Full-Size Img PowerPoint

Figure 2.  The fundamental design of the intended MLP network setup.

DownLoad: Full-Size Img PowerPoint

The regulating equations in the Lagrangian description ^[38,42,43] are as follows:

The dimensionless quantities ^[38,42,43] are as follows:

After substituting Eq (6) into the dimension form of governing Eqs (1)-(5), then the dimensionless equations ^[38,42,43] are as follows:

The Soret number is $Sr = \frac{{D}_{2}}{{\alpha }_{f}}\left(\frac{\left({T}_{h}-{T}_{c}\right)}{\left({C}_{h}-{C}_{c}\right)}\right)$, ${\delta }_{C} = \frac{{\nu }_{f}{\delta }_{1}}{{L}^{2}}$ is the Cattaneo Christov heat flux parameter, and $Du = \frac{{D}_{1}}{{\alpha }_{f}}\left(\frac{\left({C}_{h}-{C}_{c}\right)}{\left({T}_{h}-{T}_{c}\right)}\right)$ is Dufour number. The Hartmann number is $Ha = \sqrt{\frac{{\sigma }_{\mathrm{f}}}{{\mu }_{\mathrm{f}}}}{B}_{0}L$, the buoyancy ratio is $N = \frac{{\beta }_{C}\left({C}_{h}-{C}_{c}\right)}{{\beta }_{T}\left({T}_{h}-{T}_{c}\right)}$, Rayleigh's number is $Ra = \frac{\mathrm{g}{\beta }_{T}\left({T}_{h}-{T}_{c}\right){L}^{3}}{{\nu }_{f}{\alpha }_{\mathrm{f}}}$, the Prandtl number is $Pr = \frac{{\nu }_{f}}{{\alpha }_{f}},$ and Lewis's number is $Le = \frac{{\alpha }_{f}}{{D}_{m}}$.

2.1.   Dimensionless boundary conditions

Left wall, $\theta = \mathrm{\Phi } = 1, U = 0 = V,$ (17)

Right wall, $\theta = \mathrm{\Phi } = 0, U = 0 = V,$

Plane walls, $\frac{\partial \mathrm{\Phi }}{\partial Y} = \frac{\partial \theta }{\partial Y} = 0, U = 0 = V,$

Four fins, $\theta = \mathrm{\Phi } = 0, U = 0 = V,$

Wavy circular cylinder $\frac{\partial \mathrm{\Phi }}{\partial Y} = \frac{\partial \theta }{\partial Y} = 0, U = {U}_{rot}, V = {V}_{rot},$

The mean Nusselt/Sherwood numbers are as follows:

2.2.   Nanofluid properties

The nanofluid properties are as follows (^[44,45,46]):

3.   ISPH method

The solving steps in the ISPH method ^{[42,43,47,48]} are as follows:

First:

Second:

Third:

Fourth:

Fifth:

The shifting technique ^[49,50] is as follows:

The ISPH method's flowchart encompasses initiating the computational domain and boundary conditions, advancing through the time steps to estimate particle densities, solving the pressure fields, computing the heat and mass transfer, updating the particle velocities and positions, and ultimately analyzing results, all to accurately simulate the fluid dynamics. The flowchart of the ISPH method is shown in Figure 3.

Figure 3.  The flowchart illustrating the ISPH method.

DownLoad: Full-Size Img PowerPoint

4.   Validation tests

This section introduces several numerical tests of the ISPH approach, beginning with an examination of the natural convection around a circular cylinder within a square cavity. Figure 4 illustrates a comparison of isotherms during the natural convection (NC) at different Rayleigh numbers ($Ra$) between the ISPH approach and the numerical findings by Kim et al. ^[51]. The ISPH approach yielded results consistent with those obtained by Kim et al. ^[51].

Figure 4.  Comparison of the isotherms during NC between Kim et al. ^[51] (ᅳ), and the ISPH method (- - -) at $Ra = {10}^{3}, {10}^{4}, {10}^{5}$, and ${10}^{6}$.

DownLoad: Full-Size Img PowerPoint

Second, the ISPH method's validation for the natural convection in a cavity with a heated rectangle source is presented. Figure 5 compares the ISPH method with numerical simulations (Fluent software) and experimental data from Paroncini and Corvaro ^[51] for isotherms at various Rayleigh numbers and hot source heights. The ISPH method shows good agreement with both experimental and numerical results across different scenarios. Additionally, Figure 6 illustrates comparisons between the experimental, numerical, and ISPH method results for streamlines at different Rayleigh numbers and hot source lengths, confirming the accuracy of the ISPH method simulations without redundancy.

Figure 5.  The comparison of isothermal lines involves assessing (i) the present ISPH method, (ii) numerical simulation (Fluent software), and (iii) experimental data according to Paroncini and Corvaro. ^[51].

DownLoad: Full-Size Img PowerPoint

Figure 6.  Streamline comparisons involve evaluating (i) the current ISPH method, (ii) numerical simulations, and (iii) experimental data from the cited source ^[52].

DownLoad: Full-Size Img PowerPoint

5.   Results and discussion

The performed numerical simulations are discussed in this section. The range of pertinent parameters is the Cattaneo Christov heat flux parameter $\left({0\le \delta }_{c}\le 0.001\right)$, the Dufour number $\left(0\le Du\le 2\right)$, the nanoparticle parameter $\left(0\le \phi \le 0.1\right)$, the Soret number $\left(0\le Sr\le 2\right)$, the Hartmann number $\left(0\le Ha\le 80\right)$, the Rayleigh number $\left({10}^{3}\le Ra\le {10}^{5}\right)$, Fin's length $\left({0.05\le L}_{Fin}\le 0.2\right)$, and the radius of a wavy circular cylinder $\left(0.05\le {R}_{Cyld}\le 0.3\right)$. Figure 7 shows the effects of ${\delta }_{c}$ on the temperature $\theta$, concentration $\mathrm{\Phi },$ and velocity field $\mathrm{V}$. Increasing ${\delta }_{c}$ slightly enhances the temperature, concentration, and velocity field within an annulus. Figure 8 shows the effects of $Du$ on $\theta$, $\mathrm{\Phi }$, and $\mathrm{V}$. Physically, the Dufour number indicates the influence of concentration gradients on thermal energy flux. An increase in $Du$ has a minor influence on $\theta$ and $\mathrm{\Phi }$. Increasing $Du$ growths the velocity field. These results are relevant to the presence of a diabatic wavy cylinder inside a square cavity, which shrinks the contributions of $Du$. Figure 9 represents the effects of $Ha$ on $\theta$, $\mathrm{\Phi },$ and $\mathrm{V}$. The Hartmann number $\left(Ha\right)$ is a proportion of an electromagnetic force to a viscous force. Physically, the Lorentz forces represent the magnetic field that shrinks the flow velocity. Consequently, an increment in $Ha$ shrinks the distributions of θ and Φ. The maximum of the velocity field declines by $87.09\%$ as $Ha$ powers from 0 to 80. The study of the inner geometric parameters including the length of inner fins ${L}_{Fin}$ and the radius of a wavy circular cylinder ${R}_{Cyld}$ on the features of temperature, concentration, and velocity field is represented in Figures 10 and 11. In Figure 10, an expanded of ${L}_{Fin}$ supports the cooling area, which reduces the strength of $\theta$ and $\mathrm{\Phi }$ within a cavity. As the inner fins represent blockages for the nanofluid flow in a square cavity, the maximum of the velocity field reduces by $48.65\%$ as the ${L}_{Fin}$ boosts from $0.05$ to $0.2$. In Figure 11, regardless of the adiabatic condition of an inner wavy cylinder, an increment in ${R}_{Cyld}$ reduces the distributions of temperature and concentration across the cavity. Furthermore, due to the low rotating speed of the inner wavy cylinder, an expanded wavy cylinder reduces the velocity field. It is seen that the maximum of the velocity field shrinks by $55.42\%$ according to an increase in ${R}_{Cyld}$ from $0.05$ to $0.3$. Figure 12 represents the effects of $\phi$ on $\theta$, $\mathrm{\Phi }$, and $\mathrm{V}$. Here, due to the presence of an inner wavy cylinder, the contribution of $\phi$ is minor in enhancing the cooling area and reducing the strength of the concentration inside a cavity, whilst the physical contribution of $\phi$ in increasing the viscosity of a nanofluid appears well by majoring the velocity field. The velocity's maximum reduces by $36.52\%$ according to an increasing concentration of nanoparticles from 0 to 10%. Figure 13 establishes the effects of $Ra$ on $\theta$, $\mathrm{\Phi }$, and $\mathrm{V}$. As $Ra$ represents the ratio of buoyancy to thermal diffusivity and it is considering a main factor in controlling the natural convection. Accordingly, there are beneficial enhancements in the strength of the temperature and concentration in a cavity according to an increase in $Ra$. Furthermore, the Rayleigh number effectively works to accelerate the velocity field in a cavity. Figure 14 indicates the effects of $Sr$ on $\theta$, $\mathrm{\Phi },$ and $\mathrm{V}$. Physically, the Soret influence expresses diffusive movement that is created from a temperature gradient. Increasing $Sr$ enhances the strength of the temperature; this increment strongly enhances the strength of the concentration. Increasing $Sr$ from 0 to 2 supports the velocity's maximum by $134.3\%$. Figure 15 indicates the average $\stackrel{-}{Nu}$ and $\stackrel{-}{Sh}$ under the variations of $Sr$, ${R}_{Cyld}$, $Ha$, and ${\delta }_{C}$. The average $\stackrel{-}{Nu}$ is enhanced according to an increment in $Sr$, whilst it decreases alongside an increase in ${R}_{Cyld}$ and $Ha$. The average $\stackrel{-}{Sh}$ decreases alongside an increase of $Sr$, ${R}_{Cyld}$, and $Ha$. It is found that the variations of ${\delta }_{C}$ slightly changes the values of $\stackrel{-}{Nu}$ and $\stackrel{-}{Sh}$.

Figure 7.  The effects of ${\delta }_{c}$ on $\theta$, $\mathrm{\Phi }$ and $\mathrm{V}$ at $N = 2, \gamma = 45°,$ $Du = 0.6$, $Sr = 1.2$, $\phi = 0.05$, $Ha = 10$, ${L}_{Fin} = 0.1, {R}_{Cyld} = 0.25$ and $Ra = {10}^{4}$.

DownLoad: Full-Size Img PowerPoint

Figure 8.  The effects of $Du$ on $\theta$, $\mathrm{\Phi }$ and $\mathbf{V}$ at $N = 2, \gamma = 45°,$ ${\delta }_{c} = 0.0001$, $Sr = 1.2$, $\phi = 0.05$, $Ha = 10$, ${L}_{Fin} = 0.1, {R}_{Cyld} = 0.25$ and $Ra = {10}^{4}$.

DownLoad: Full-Size Img PowerPoint

Figure 9.  The effects of $Ha$ on $\theta$, $\mathrm{\Phi }$ and $\mathbf{V}$ at $N = 2, \gamma = 45°,$ ${\delta }_{c} = 0.0001$, $Sr = 1.2$, $\phi = 0.05$, $Du = 0.6$, ${L}_{Fin} = 0.1, {R}_{Cyld} = 0.25$ and $Ra = {10}^{4}$.

DownLoad: Full-Size Img PowerPoint

Figure 10.  The effects of ${L}_{Fin}$ on $\theta$, $\mathrm{\Phi }$ and $\mathbf{V}$ at $N = 2, \gamma = 45°,$ ${\delta }_{c} = 0.0001$, $Sr = 1.2$, $\phi = 0.05$, $Du = 0.6$, $Ha = 10, {R}_{Cyld} = 0.25$ and $Ra = {10}^{4}$.

DownLoad: Full-Size Img PowerPoint

Figure 11.  The effects of ${R}_{Cyld}$ on $\theta$, $\mathrm{\Phi }$ and $\mathbf{V}$ at $N = 2, \gamma = 45°,$ ${\delta }_{c} = 0.0001$, $Sr = 1.2$, $Ha = 10$, $Du = 0.6$, ${L}_{Fin} = 0.1, \phi = 0.05,$ and $Ra = {10}^{4}$.

DownLoad: Full-Size Img PowerPoint

Figure 12.  The effects of $\phi$ on $\theta$, $\mathrm{\Phi }$ and $\mathbf{V}$ at $N = 2, \gamma = 45°,$ ${\delta }_{c} = 0.0001$, $Sr = 1.2$, $Ha = 10$, $Du = 0.6$, ${L}_{Fin} = 0.1, {R}_{Cyld} = 0.25$ and $Ra = {10}^{4}$.

DownLoad: Full-Size Img PowerPoint

Figure 13.  The effects of $Ra$ on $\theta$, $\mathrm{\Phi }$, and $\mathbf{V}$ at $N = 2, \gamma = 45°,$ ${\delta }_{c} = 0.0001$, $Sr = 1.2$, $Ha = 10$, $Du = 0.6$, ${L}_{Fin} = 0.1, \phi = 0.05,$ and ${R}_{Cyld} = 0.25$.

DownLoad: Full-Size Img PowerPoint

Figure 14.  The effects of $Sr$ on $\theta$, $\mathrm{\Phi }$, and $\mathbf{V}$ at $N = 2, \gamma = 45°,$ ${\delta }_{c} = 0.0001$, $Ha = 10$, $Du = 0.6$, ${L}_{Fin} = 0.1, \phi = 0.05,$ $Ra = {10}^{4},$ and ${R}_{Cyld} = 0.25$.

DownLoad: Full-Size Img PowerPoint

Figure 15.  Average $\stackrel{-}{Nu}$ and $\stackrel{-}{Sh}$ under the variations of $Sr$, ${R}_{Cyld}$, $Ha$, and ${\delta }_{C}$.

DownLoad: Full-Size Img PowerPoint

Figure 16 displays the performance chart that was made to examine the produced ANN model's training phase. The graph illustrates how the MSE values, which were large at the start of the training phase, got smaller as each epoch progressed. By achieving the optimal validation value for each of the three data sets, the MLP model's training phase was finished. The error histogram derived from the training phase's data is displayed in Figure 17. Upon scrutinizing the error histogram, one observes that the error values are predominantly centered around the zero-error line. Additionally, relatively low numerical values for the errors were achieved. In Figure 18, the regression curve for the MLP model is displayed, indicating that R = 1 is consistently employed for training, validation, and testing. According to the results of the performance, regression, and error histograms, the built ANN model's training phase was successfully finished. Following the training phase's validation, the ANN model's prediction performance was examined. First, the goal values and the $\stackrel{-}{Nu}$ and $\stackrel{-}{Sh}$ values derived by the ANN model were compared in this context. The goal values for every data point are displayed in Figure 19, along with the $\stackrel{-}{Nu}$ and $\stackrel{-}{Sh}$ values that the ANN model produced. Upon a closer inspection of the graphs, it becomes evident that the $\stackrel{-}{Nu}$ and $\stackrel{-}{Sh}$ values derived by the ANN model perfectly align with the desired values.

Figure 16.  The performance of MLP model during training processes.

DownLoad: Full-Size Img PowerPoint

Figure 17.  The error histogram for the MLP model.

DownLoad: Full-Size Img PowerPoint

Figure 18.  The regression outline.

DownLoad: Full-Size Img PowerPoint

Figure 19.  The ANN model and target values yielded the values of $\stackrel{-}{Nu}$ and $\stackrel{-}{Sh}$ for every data point.

DownLoad: Full-Size Img PowerPoint

6.   Conclusions

The current work numerically examined the effects of a Cattaneo-Christov heat flux on the thermosolutal convection of a nanofluid-filled square cavity. The square cavity had four fins located on the cavity's boarders and a wavy circular cylinder, which was located in the cavity's center, and rotated in a circular form. This study checked the impacts of the Cattaneo Christov heat flux parameter, the Dufour number, a nanoparticle parameter, the Soret number, the Hartmann number, the Rayleigh number, Fin's length, and the radius of a wavy circular cylinder on the temperature, concentration, and velocity field. The results are highlighted within the following points:

• The parameters of the embedded blockages, including the fins length and radius of a wavy circular cylinder, worked well to improve the cooling area and reduce the strength of concentration as well as the velocity field.

• Adding an extra concentration of nanoparticles up to 10% enhanced the cooling area and reduced the maximum velocity field by 36.52%.

• The Rayleigh number beneficially contributed to power the double diffusion within an annulus and was effective in accelerating the nanofluid velocity.

• The Soret number effectively contributed to change the distributions of the temperature, the concentration, and the velocity field within an annulus as compared to the Dufour number.

• The excellent accord between the prediction values produced by the ANN model and the goal values indicated that the developed ANN model could forecast the $\stackrel{-}{Nu}$ and $\stackrel{-}{Sh}$ values with a remarkable precision.

Author contributions

Munirah Alotaibi: Conceptualization, Funding acquisition, Methodology, Writing – original draft; Abdelraheem M. Aly: Formal Analysis, Investigation, Validation, Visualization, Writing – review & editing. All authors have read and approved the final version of the manuscript for publication.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors extend their appreciation to the Deanship of Scientific Research at King Khalid University, Abha, Saudi Arabia, for funding this work through the Research Group Project under Grant Number (RGP. 2/38/45). This research was funded by the Princess Nourah bint Abdulrahman University Researchers Supporting Project number (PNURSP2024R522), Princess Nourah bint Abdulrahman University, Riyadh, Saudi Arabia.

Conflict of interest

The authors declare that they have NO affiliations with or involvement in any organization or entity with any financial interest.