Enhanced thermal energy storage using micropolar fluids: A numerical study

Abstract

This study develops and solves a micropolar-fluid model for heat storage in a rectangular duct subject to constant wall heat flux. The dimensionless mass, linear and angular momentum, and energy equations are treated with a Runge-Kutta (for the coupled momentum ODEs) and finite-difference scheme (for the energy PDE). Increasing the coupling parameter raises the dimensionless axial velocity while initially reducing and then increasing microrotation; the net effect is a decrease in dimensionless temperature, indicating diminished storage effectiveness. In contrast, higher spin-gradient viscosity reverses these trends and enhances thermal storage. Reynolds and Prandtl numbers exhibit inverse relationships with the temperature profile. Quantitatively, the mean Nusselt number is $\approx 2.53$ at an entrance length of ${\mathrm{x}}^{*}\approx 5.5$ , while the local Nusselt number decays approximately exponentially downstream and approaches zero. For a low-inertia case ( $\mathrm{Re}=5$ ), the peak fluid temperature reached $\sim 95\circ \mathrm{C}$ under the imposed heat flux. These results clarify how micropolar parameters govern heat transfer and storage performance, offering guidance for tuning ducts that use micropolar working fluids.

Introduction

As global energy needs grow, demand for energy services rises, driving a shift to renewable energy, particularly solar, in response to climate change and declining traditional resources. Fossil fuel use has increased emissions, environmental issues, geopolitical tensions, and fuel prices. These challenges pose significant threats to society, making sustainable energy sources essential (Asfar et al., 2023; Khan et al., 2025). Countries are shifting to clean energy, with solar power being the key to its abundance and sustainability (Al-Ghussain et al., 2025). Technological advances and lower costs have made solar energy a vital renewable resource for electricity and water heating. Switching to solar energy benefits the environment, economy, and creates jobs, while driving growth, innovation, and technological progress (Alrbai et al., 2023, 2022).

Solar energy faces challenges such as interruptions caused by sunlight variations and weather changes. These fluctuations affect its reliability, especially for sectors requiring a continuous power supply (Shutaywi et al., 2024). Effective storage systems are therefore essential to store excess energy during peak sunlight and supply it during low-radiation periods (Khan et al., 2024a; Shutaywi et al., 2025). Thermal energy storage (TES) allows short-term retention of thermal energy at varying temperatures, contributing significantly to energy savings and enabling low-cost thermal energy and power generation (Khan et al., 2024b; Raza et al., 2024a). Its importance has grown in addressing supply–demand mismatches in intermittent renewable energy (Khan et al., 2023b; Raza et al., 2024b). Well-designed TES systems minimize thermal losses while maximizing energy extraction efficiency. TES technologies include sensible heat storage, latent heat storage using phase change materials (PCMs), and thermochemical storage (Raza et al., 2024b; Raza and Ashfaq, 2025). TES plays a crucial role in applications such as concentrated solar power and HVAC systems, with recent advancements highlighting its potential to transform energy technologies (Khan et al., 2023c; Rathore et al., 2025). Among solar thermal systems, water heating is widely adopted worldwide, reducing reliance on fossil fuels and greenhouse gas emissions (Khan et al., 2022; Kumar et al., 2025). Evacuated tube collectors (ETCs) achieve higher efficiencies than flat-plate collectors (FPCs) by employing vacuum insulation to minimize conduction and convection heat losses (Anantha Kumar et al., 2023; Kumar et al., 2019). ETCs can absorb up to 80% of solar energy compared to about 60% for FPCs, making them suitable for both heating and cooling applications. Their wide operating temperature range (–30 to 130 °C) enhances durability and reduces maintenance needs (Anantha Kumar, 2024; Anantha Kumar et al., 2025).

In addition to collector design, the choice of working fluid plays a central role in enhancing the efficiency of thermal energy storage systems. Micropolar fluids extend the classical Navier–Stokes framework by incorporating the effects of microrotation and microstructural interactions within the fluid (Khan et al., 2023d; Ramadevi et al., 2020). Unlike Newtonian fluids, which assume constant viscosity and neglect particle-scale effects, micropolar fluids account for the presence of rigid or semi-rigid microelements suspended in the medium (Khan et al., 2023a; Ramadevi et al., 2019a). These particles generate asymmetric stress tensors and additional rotational degrees of freedom, which fundamentally modify the flow and thermal transport characteristics (Ramadevi et al., 2019b). As a result, micropolar fluids can exhibit lower skin friction, altered thermal boundary-layer thickness, and distinctive velocity and temperature profiles compared to conventional fluids (Buruju et al., 2021). In the context of TES, these features are particularly valuable, as they allow for more precise control over energy storage and release within ducts and channels. Parameters such as the coupling constant (which measures the strength of interaction between linear velocity and microrotation) and the spin-gradient viscosity (which characterizes resistance to microrotation) strongly influence the thermal and hydrodynamic response of micropolar fluids. Increasing the coupling parameter, for instance, can enhance velocity but reduce temperature, thereby limiting storage effectiveness, while higher spin-gradient viscosity tends to have the opposite effect by improving thermal retention. Furthermore, dimensionless groups such as the Reynolds and Prandtl numbers interact with micropolar parameters to shape the Nusselt number distribution along the duct, which is a key indicator of heat transfer performance. These distinctive behaviors underscore why micropolar fluids are considered a promising candidate for optimizing solar thermal energy storage, offering pathways to improve reliability, efficiency, and adaptability in systems facing intermittency challenges.

Several studies have explored improving ETC performance. For example, Wang et al. (Wang et al., 2020) showed that transparent-tube collectors (T-TC) achieved higher thermal efficiency (85.0%) than conventional-tube collectors (C-TC, 77.6%), yet the latter produced more useful energy per collector area (641 W vs. 497 W). Gholipour et al. (2020) demonstrated that vacuum tube collectors achieve higher thermal efficiency at greater flow rates (10–40 L/h), with helical coil configurations outperforming spiral and U-tube designs. Using CFD, Shah and Furbo Shah and Furbo (2007) found that shorter tube lengths enhanced efficiency, with optimal inlet flow rates between 0.4 and 1 kg/min.

Research has also focused on optimizing working fluids in ETCs. López-Núñez et al. (2022) reported that a TiO₂-water nanofluid improved outlet temperature, velocity, and thermal efficiency while lowering entropy generation. Khair and Duwairi (2021) studied Al₂O₃–graphite nanofluids in porous media and found improved performance compared to water. Similarly, Maayah et al. (2021) investigated a graphite–CMC non-Newtonian fluid in a circular conduit, showing that reduced porosity improved storage by raising fluid temperature and effectiveness.

Fluids are broadly categorized as Newtonian or non-Newtonian. Newtonian fluids, such as water and air, exhibit constant viscosity regardless of shear rate (Buruju et al., 2021). Non-Newtonian fluids, however, have viscosity dependent on applied shear, producing nonlinear behaviour (Anantha Kumar et al., 2023), which makes them advantageous for specialized applications (Anantha Kumar, 2024). Classical fluid models struggle to capture the behavior of complex fluids like polymers, colloids, and emulsions at micro- and nanoscales (Ramadevi et al., 2019b).

Micropolar fluids show promise for thermal energy storage by enhancing heat transfer capacity, heat absorption, and release characteristics. For instance, Duwairi and Chamkha (2005) analyzed boundary-layer flow near an isothermal vertical surface and found that higher vortex viscosity increased friction but reduced heat transfer. Rahman et al. (2009) reported that micropolar fluids reduce both skin friction and heat transfer compared to Newtonian fluids. Qasim et al. (2013) studied micropolar flow over a stretching surface with Newtonian heating, showing that increasing vortex viscosity (K) and Prandtl number (Pr) decreased temperature and thermal boundary-layer thickness. Yadav et al. (2018) investigated micropolar and Newtonian fluids in porous channels, finding that micropolar fluids exhibited lower velocities but could be controlled through permeability and viscosity ratios. Rana et al. (2020) demonstrated that hybrid Cu–Al₂O₃ nanoparticles in micropolar fluids enhanced heat transfer more than Cu alone, although higher vortex and spin-gradient viscosities reduced temperature. Pasha et al. (2022) analytically studied micropolar flow and heat transfer between plates, concluding that higher Peclet numbers improved both heat transfer and concentration. Jalili et al. (2023) analyzed nonlinear 2D micropolar flow between porous barriers, showing that increasing the coupling parameter raised most dimensionless parameters, while Peclet number strongly influenced temperature and concentration. Finally, Al-Sharifi et al. (2023) studied Casson micropolar flow with convective boundary conditions, finding that larger Casson and material parameters increased velocity but reduced temperature and microrotation, whereas higher Prandtl numbers decreased all three quantities.

Despite extensive research on thermal energy storage systems, the behavior of micropolar fluids in vacuum rectangular ducts remains largely unexplored, particularly within the context of solar energy storage applications. Most existing studies focus on conventional Newtonian or simple nanofluids, neglecting the distinctive rotational characteristics of micropolar fluids that critically influence momentum transport, energy dissipation, and heat transfer efficiency. Furthermore, while several numerical models exist for forced convection, few have incorporated the Modified Navier–Stokes framework to capture the coupled effects of microrotation and linear velocity under steady-state conditions.

The present study addresses these gaps by developing a comprehensive numerical model for micropolar fluid flow and heat transfer in a vacuum rectangular duct subjected to constant surface heat flux. By formulating the continuity, linear and angular momentum, and energy equations in dimensionless form and solving them using a hybrid Runge–Kutta and finite difference approach in MATLAB, this work provides a detailed characterization of transitional velocity, microrotation, temperature distribution, and local and mean Nusselt numbers along the duct. Unlike prior studies, this research explicitly examines the combined influence of coupling parameters, spin-gradient viscosity, Reynolds number, and Prandtl number on thermal energy storage performance. The findings reveal previously unreported interactions between micropolar fluid properties and heat transfer dynamics in confined ducts, offering quantitative insights into optimizing energy storage efficiency in solar thermal systems. This study therefore establishes a novel framework for assessing advanced heat transfer fluids in TES applications and provides a foundation for future experimental validation and practical implementation.

Mathematical model

Forced convection heat transfer moves heat from a solid surface to a fluid through forced motion, using pumps, fans, or natural forces like wind. The temperature difference creates a convective heat transfer coefficient, which is higher in forced convection due to enhanced mixing and turbulence. This efficiency is influenced by factors such as fluid velocity, temperature, properties, surface geometry, and roughness. Present paper research focuses on forced convection heat transfer in a duct filled with micropolar fluid, specifically for solar energy applications under constant heat flux in the thermal engineering field. This section mathematically formulates the case of a micropolar fluid under constant solar heat flux in a two-dimensional duct, focusing on convective heat transfer influenced by the fluid's unique properties. These unique characteristics are of great interest in solar energy applications, where efficient heat transfer is essential for optimal energy conversion (Jalili et al., 2023). The spatial structure of the vacuum rectangular duct is two-dimensional, with length and height as the principal directions, forming a rectangular channel. The duct's length, L, is along the x-axis, while its height, 2H, is along the y-axis. The choice of these dimensions depends on the specific requirements of the research and the desired aspect ratio of the duct. The duct can be made of metal or plastic to withstand the study's temperature and pressure conditions. The duct's smooth surface ensures the “no-slip” condition, essential for accurate fluid flow prediction. Dimensions are set by experimental needs, flow velocities (u and v), and heat flux (qs″). Inlet and outlet sections control fluid flow, with no change in cross-section in this study. As Figure 1 illustrates, the spatial structure of the duct in this system is a two-dimensional rectangular-like structure, with defined length and height dimensions, where the micropolar fluid flows and heat transfer processes occur.

Figure 1.

The spatial structure of the duct.

The methodology employed in this study offers several key advantages. By combining the Runge–Kutta fourth-order method for the coupled momentum equations with finite difference methods for the energy equation, the approach ensures high accuracy and stability in solving complex micropolar fluid flows. The dimensionless formulation allows generalization of results and systematic analysis of parameters such as microrotation, coupling constant, spin-gradient viscosity, Reynolds number, and Prandtl number. This comprehensive treatment enables precise prediction of velocity, microrotation, temperature distribution, and Nusselt numbers, providing detailed insights into the thermal energy storage performance of vacuum rectangular ducts.

The present study considers laminar, incompressible, and two-dimensional flow of a micropolar fluid within a vacuum rectangular duct. The fluid is assumed to exhibit microrotation effects and spin-gradient viscosity, while its physical properties, including density, viscosity, thermal conductivity, and specific heat, are considered constant. The duct walls are maintained under a constant surface heat flux, and no-slip and no-spin boundary conditions are applied. For the steady-state analysis, the flow is assumed fully developed, whereas entrance-region effects and daily solar irradiance variations are treated separately, allowing the examination of flow development and transient temperature profiles without altering the momentum equations. Body forces, buoyancy, and radiative effects are neglected, ensuring a focused and consistent analysis of the thermal energy storage performance of micropolar fluids. In the daily variation analysis, a quasi-steady approach is employed, where the steady-state equations are solved repeatedly using time-dependent surface heat flux values to simulate changing solar irradiance throughout the day. While this method does not capture true transient effects, such as fluid inertia or dynamic thermal response, it provides a practical and computationally efficient way to estimate temperature variations and thermal energy storage trends under realistic daily solar conditions. The results highlight periods of peak heat storage and allow evaluation of optimal micropolar fluid parameters for different times of the day.

Navier-Stokes equations can't model micropolar fluids, which have microstructures with molecules that can deform and rotate independently, the micropolar fluid model is a meaningful extension of the traditional Navier-Stokes model. It introduces an angular velocity field for particle rotation, and an additional vector equation for the conservation of angular momentum where the stress tensor lacks symmetry, the conservation law of angular momentum operates independently. This leads to the introduction of an extra equation, the micropolar fluid model includes the micro-rotation equation, enabling analytical solutions based on viscosity coefficients. It distinguishes between dynamic microrotation viscosity, measuring resistance to microrotation, and dynamic Newtonian viscosity, measuring shear resistance, resulting in two types of viscosity: linear and angular. Below, we propose the governing equations for the linear and angular movement and energy of the micropolar fluid in the described system, including the continuity equation for the fully developed region:

The continuity equation (Naterer, 2022):

$${\displaystyle \frac{\mathrm{\partial }\mathrm{u}}{\mathrm{\partial }\mathrm{x}}+{\displaystyle \frac{\mathrm{\partial }\mathrm{v}}{\mathrm{\partial }\mathrm{y}}=0}}$$

(1)

The linear momentum equation (Naterer, 2022):

$$\rho \left(\mathrm{u}{\displaystyle \frac{\mathrm{\partial }\mathrm{u}}{\mathrm{\partial }\mathrm{x}}+\mathrm{v}{\displaystyle \frac{\mathrm{\partial }\mathrm{u}}{\mathrm{\partial }\mathrm{y}}}}\right)=-{\displaystyle \frac{\mathrm{\partial }\mathrm{p}}{\mathrm{\partial }\mathrm{x}}+\mu \left({\displaystyle \frac{{\mathrm{\partial }}^2\mathrm{u}}{\mathrm{\partial }{\mathrm{x}}^2}+{\displaystyle \frac{{\mathrm{\partial }}^2\mathrm{u}}{\mathrm{\partial }{\mathrm{y}}^2}}}\right)+\kappa {\displaystyle \frac{\mathrm{\partial }\mathrm{w}}{\mathrm{\partial }\mathrm{y}}}}$$

(2)

The angular momentum equation (Naterer, 2022):

$$\rho \mathrm{J}\left(\mathrm{u}{\displaystyle \frac{\mathrm{\partial }\mathrm{w}}{\mathrm{\partial }\mathrm{x}}+\mathrm{v}{\displaystyle \frac{\mathrm{\partial }\mathrm{w}}{\mathrm{\partial }\mathrm{y}}}}\right)=\kappa \left({\displaystyle \frac{\mathrm{\partial }\mathrm{v}}{\mathrm{\partial }\mathrm{x}}-{\displaystyle \frac{\mathrm{\partial }\mathrm{u}}{\mathrm{\partial }\mathrm{y}}-2\mathrm{w}}}\right)+\gamma {\displaystyle \frac{{\mathrm{\partial }}^2\mathrm{w}}{\mathrm{\partial }{\mathrm{y}}^2}}$$

(3)

The energy equation (Naterer, 2022):

$$\left(\mathrm{u}{\displaystyle \frac{\mathrm{\partial }\mathrm{T}}{\mathrm{\partial }\mathrm{x}}+\mathrm{v}{\displaystyle \frac{\mathrm{\partial }\mathrm{T}}{\mathrm{\partial }\mathrm{y}}}}\right)={\alpha }_{\mathrm{eff}}\left({\displaystyle \frac{{\mathrm{\partial }}^2\mathrm{T}}{\mathrm{\partial }{\mathrm{x}}^2}+{\displaystyle \frac{{\mathrm{\partial }}^2\mathrm{T}}{\mathrm{\partial }{\mathrm{y}}^2}}}\right)$$

(4)

The following initial and boundary conditions (that are depended on both the hydrodynamic and thermal entrance lengths, considering that these two lengths are not equal in generally).

At $={\mathrm{x}}_{\mathrm{fd},\mathrm{t}}$ ;

$$\mathrm{T}=\mathrm{T}(\mathrm{y},{\mathrm{x}}_{\mathrm{fd},\mathrm{t}})$$

(5)

At $>{\mathrm{x}}_{\mathrm{fd},\mathrm{t}}$ ;

$$\begin{array}{c}{\displaystyle \frac{\mathrm{\partial }\mathrm{T}}{\mathrm{\partial }\mathrm{y}}=0\mathrm{at}\mathrm{y}=0}\\ \mathrm{K}{\displaystyle \frac{\mathrm{\partial }\mathrm{T}}{\mathrm{\partial }\mathrm{y}}={\mathrm{q}}_{\mathrm{s}}^{\mathrm{″}}\mathrm{at}\mathrm{y}=\mathrm{H}}\end{array}$$

(6)

At $\ge {\mathrm{x}}_{\mathrm{fd},\mathrm{h}}$ ;

$$\begin{array}{c}{\displaystyle \frac{\mathrm{\partial }\mathrm{u}}{\mathrm{\partial }\mathrm{y}}=0,{\displaystyle \frac{\mathrm{\partial }\mathrm{w}}{\mathrm{\partial }\mathrm{y}}=0,\mathrm{at}\mathrm{y}=0}}\\ \mathrm{w}=-\mathrm{n}{\displaystyle \frac{\mathrm{\partial }\mathrm{u}}{\mathrm{\partial }\mathrm{y}},\mathrm{u}=0,\mathrm{at}\mathrm{y}=\mathrm{H}}\end{array}$$

(7)

The effects of the micropolar fluid parameters are closely tied to the governing equations. For instance, the coupling parameter directly multiplies the microrotation terms in the momentum equations, influencing the velocity distribution and generating additional rotational stresses that alter energy transport. Similarly, the spin-gradient viscosity appears in the angular momentum equation, enhancing rotational diffusion and thereby increasing the local temperature through improved convective heat transfer. Reynolds and Prandtl numbers scale the convective and diffusive terms in the energy equation, governing the relative influence of advection and conduction. Scaling arguments show that higher coupling parameters increase fluid momentum but reduce thermal penetration depth, while higher spin-gradient viscosity enhances rotational energy dissipation and raises the mean temperature, explaining the trends observed in the numerical results.

The micro-gyration constant n ranges from 0 to 1, relating microrotation to shear stress. n=0 means no rotation; n=0.5 eliminates the asymmetric stress tensor; n=1 indicates turbulent flow. To simplify the governing equations for the micropolar fluid system, we make the following assumptions: steady-state, incompressible fluid with uniform properties in all directions, flow driven by an external force (e.g., a pump), laminar flow, negligible radiation effects, no thermal or mechanical energy generation, fully developed flow, constant surface heat flux, negligible viscous dissipation, and temperature below boiling point, The micropolar fluid pressure depends only on the $\mathrm{x}$ -direction without changing in the direction of $\mathrm{y}$ , i.e., $\frac{\mathrm{\partial }\mathrm{p}}{\mathrm{\partial }\mathrm{y}}\approx 0$ , The micro-rotation vector $\mathrm{w}$ depends only on $\mathrm{y}$ -axis, i.e., $\frac{\mathrm{\partial }\mathrm{w}}{\mathrm{\partial }\mathrm{x}}\approx 0$ , and The negligible net axial conduction in the fully developed region, $\frac{{\mathrm{\partial }}^2\mathrm{T}}{\mathrm{\partial }{\mathrm{x}}^2}\approx 0$ , where $\mathrm{T}$ is the temperature. The assumptions of steady-state, laminar flow, and incompressibility were chosen to simplify the modeling process and focus on the key thermal and fluid dynamic behaviors relevant to thermal energy storage applications. The steady-state assumption allows for the analysis of long-term system behavior, where transient effects are not considered. Laminar flow was assumed due to the relatively low Reynolds number in typical thermal energy storage applications, ensuring that the flow remains smooth and predictable. Incompressibility was assumed based on the low variation in density of the fluid, which is appropriate for micropolar fluids in such systems.

Hence, we have four differential partial equations to deal with. In the following we formulate these governing equations for fully developed regions with different set of boundary conditions, and with the following dimensionless variables definitions (Naterer, 2022):

$$\begin{array}{rl}& {\mathrm{u}}^{*}={\displaystyle \frac{\mathrm{u}}{{\mathrm{u}}_{\mathrm{m}}},{\mathrm{x}}^{*}={\displaystyle \frac{\mathrm{x}}{\mathrm{H}},{\mathrm{y}}^{*}={\displaystyle \frac{\mathrm{y}}{\mathrm{H}},{\mathrm{w}}^{*}={\displaystyle \frac{\mathrm{w}}{{\mathrm{w}}_{\mathrm{m}}},\theta =\left[{\displaystyle \frac{\mathrm{T}-\mathrm{Te}}{{\mathrm{q}}_{\mathrm{s}}^{\mathrm{″}}.\mathrm{H}{\displaystyle \frac{1}{\mathrm{K}}}}}\right]}}}}\\ & {\mathrm{u}}_{\mathrm{m}}={\displaystyle \frac{{\mathrm{H}}^2.\mathrm{p}}{\mu },{\mathrm{w}}_{\mathrm{m}}={\displaystyle \frac{{\mathrm{u}}_{\mathrm{m}}}{\mathrm{H}},\mathrm{p}=-{\displaystyle \frac{\mathrm{\partial }\mathrm{p}}{\mathrm{\partial }\mathrm{x}}=\mathrm{constant},{\hat{\mathrm{k}}}^{*}={\displaystyle \frac{\kappa }{\mu },\xi ={\displaystyle \frac{\gamma }{\mu \cdot {\mathrm{H}}^2}}}}}}\end{array}$$

(8)

Where ${\mathrm{u}}^{*},{\mathrm{x}}^{*}$ , ${\mathrm{y}}^{*},{\mathrm{w}}^{*},\theta$ are dimensionless parameters velocity, x-y coordinate x y, microrotation, and temperature.

Where ${\mathrm{u}}_{\mathrm{m}}$ is the mean value of the velocity component of $\mathrm{u}$ , $\mathrm{H}$ is the high of the duct considered from its center, ${\mathrm{w}}_{\mathrm{m}}$ is the mean value of the micro-rotation angular velocity, ${\mathrm{q}}_{\mathrm{s}}^{″}$ is the value of the constant heat flux, and is the uniform initial temperature at the beginning of the duct, and ${\mathrm{T}}_{\mathrm{e}}=\mathrm{T}(0,\mathrm{y})$ is the initial uniform temperature at $\mathrm{x}=0$ , ${\hat{\mathrm{k}}}^{*}$ is coupling parameter, $\xi$ is the spin gradient viscosity parameter.

The important feature of hydrodynamic conditions in the fully developed region is that the radial velocity component $\mathrm{v}$ is zero ( $\frac{\mathrm{\partial }\mathrm{v}}{\mathrm{\partial }\mathrm{y}})$ , and hence according to the equation (2), the gradient of the axial velocity component ( $\frac{\mathrm{\partial }\mathrm{u}}{\mathrm{\partial }\mathrm{x}}$ ) become also zero. The equations (9, 10, 11) represent the final governing dimensionless equations that will be solved together in the next section, and this will be done numerically regarding the following dimensionless initial and boundary conditions (that are depended on both the hydrodynamic and thermal entrance lengths, considering that these two lengths are not equal in generally):

$$\frac{{\mathrm{\partial }}^2{\mathrm{u}}^{*}}{\mathrm{\partial }{{\mathrm{y}}^{*}}^2}+{\hat{\mathrm{k}}}^{*}\frac{\mathrm{\partial }{\mathrm{w}}^{*}}{\mathrm{\partial }{\mathrm{y}}^{*}}+1=0$$

(9)

$${\hat{\mathrm{k}}}^{*}.\frac{\mathrm{\partial }{\mathrm{u}}^{*}}{\mathrm{\partial }{\mathrm{y}}^{*}}+2{\hat{\mathrm{k}}}^{*}.{\mathrm{w}}^{*}=\zeta .\frac{{\mathrm{\partial }}^2{\mathrm{w}}^{*}}{\mathrm{\partial }{{\mathrm{y}}^{*}}^2}$$

(10)

$$\frac{{\mathrm{\partial }}^2\theta }{\mathrm{\partial }{{\mathrm{y}}^{*}}^2}=\mathrm{Re}.{\mathrm{p}}_{\mathrm{r}}.{\mathrm{u}}^{*}.\frac{\mathrm{\partial }\theta }{\mathrm{\partial }{\mathrm{x}}^{*}}$$

(11)

At ${\mathrm{x}}^{*}=\frac{{\mathrm{x}}_{\mathrm{fd},\mathrm{t}}}{\mathrm{H}}$ ;

$$\theta =\theta \left({\mathrm{y}}^{*},{\mathrm{x}}^{*}={\displaystyle \frac{{\mathrm{x}}_{\mathrm{fd},\mathrm{t}}}{\mathrm{H}}}\right)$$

(12)

At ${\mathrm{x}}^{*}>\frac{{\mathrm{x}}_{\mathrm{fd},\mathrm{t}}}{\mathrm{H}}$ ;

$$\begin{array}{c}{\displaystyle \frac{\mathrm{\partial }\theta }{\mathrm{\partial }{\mathrm{y}}^{*}}=0\mathrm{at}{\mathrm{y}}^{*}=0}\\ {\displaystyle \frac{\mathrm{\partial }\theta }{\mathrm{\partial }{\mathrm{y}}^{*}}=1\mathrm{at}{\mathrm{y}}^{*}=1}\end{array}$$

(13)

At ${\mathrm{x}}^{*}\ge \frac{{\mathrm{x}}_{\mathrm{fd},\mathrm{h}}}{\mathrm{H}}$ ;

$$\begin{array}{c}{\displaystyle \frac{\mathrm{\partial }{\mathrm{u}}^{*}}{\mathrm{\partial }{\mathrm{y}}^{*}}=0,{\displaystyle \frac{\mathrm{\partial }{\mathrm{w}}^{*}}{\mathrm{\partial }{\mathrm{y}}^{*}}=0,\mathrm{at}{\mathrm{y}}^{*}=0}}\\ {\mathrm{w}}^{*}=-\mathrm{n}.{\displaystyle \frac{\mathrm{\partial }{\mathrm{u}}^{*}}{\mathrm{\partial }{\mathrm{y}}^{*}},{\mathrm{u}}^{*}=0,\mathrm{at}{\mathrm{y}}^{*}=1}\end{array}$$

(14)

Numerical solution

For clarity and to demonstrate the realism of the nondimensional parameters, representative dimensional values are provided. The rectangular duct is assumed to have a height of 0.05 m and a width of 0.1 m. The micropolar fluid has a density of 1000 kg/m³, dynamic viscosity of 0.001 Pa·s, and thermal conductivity of 0.6 W/m·K. The surface heat flux is taken as 5000 W/m², consistent with typical solar thermal applications. Using these values, the dimensionless parameters—Reynolds number, Prandtl number, coupling parameter, and spin-gradient viscosity—fall within physically realistic ranges relevant to solar thermal energy storage scenarios.

This study employed numerical methods to solve the governing equations for steady state. The Runge-Kutta fourth order (RK4) method was used for coupled ordinary differential equations to find dimensionless velocity and microrotation, while the finite difference method determined dimensionless temperature based on boundary and initial conditions. All solutions were derived under constant surface heat flux conditions.

Runge Kutta Method (Solution of coupled linear & angular momentum):

Solving the coupled non-dimensional equations of linear and angular momentum, represented in equations (9) and (10), using the Runge-Kutta method, specifically the fourth-order (RK4) method. This approach helps us determine the dimensionless transitional velocity (u*) and the dimensionless microrotation (w*) which are depend on y*.

The fourth-order Runge-Kutta (RK4) method is a popular numerical technique for solving ordinary differential equations (ODEs). It efficiently approximates solutions with high accuracy, making it a preferred choice for many applications., It uses four slope estimates (K) at each step to compute the next value, with a step interval of (h = 0.01).

Here are the steps detailing the use of this method:

Introducing a new variables to transform the higher-order system equation (second order) into a first-order equation, for equation (9) and (10).

$$\frac{{\mathrm{d}}^2{\mathrm{u}}^{*}}{\mathrm{d}{{\mathrm{y}}^{*}}^2}+{\hat{\mathrm{k}}}^{*}\frac{\mathrm{d}{\mathrm{w}}^{*}}{\mathrm{d}{\mathrm{y}}^{*}}+1=0$$

(15)

$${\hat{\mathrm{k}}}^{*}.\frac{\mathrm{\partial }{\mathrm{u}}^{*}}{\mathrm{\partial }{\mathrm{y}}^{*}}+2{\hat{\mathrm{k}}}^{*}.{\mathrm{w}}^{*}=\zeta .\frac{{\mathrm{\partial }}^2{\mathrm{w}}^{*}}{\mathrm{\partial }{{\mathrm{y}}^{*}}^2}$$

(16)

$$U_1=u*,U_2={\displaystyle \frac{{\mathrm{du}}^{*}}{{\mathrm{dy}}^{*}},W_1=w*,W_2={\displaystyle \frac{{\mathrm{dw}}^{*}}{{\mathrm{dy}}^{*}}}}$$

(17)

The system becomes:

$$\frac{\mathrm{d}{\mathrm{U}}_1}{\mathrm{d}{\mathrm{y}}^{*}}=U_{2,}\frac{\mathrm{d}{\mathrm{U}}_2}{\mathrm{d}{\mathrm{y}}^{*}}={\frac{{\mathrm{d}}^2{\mathrm{u}}^{*}}{\mathrm{d}{{\mathrm{y}}^{*}}^2}}_{,}{W}_2=\frac{\mathrm{d}{\mathrm{W}}_1}{\mathrm{d}{\mathrm{y}}^{*}},\frac{\mathrm{d}{\mathrm{W}}_2}{\mathrm{d}{\mathrm{y}}^{*}}=\frac{{\mathrm{d}}^2{\mathrm{w}}^{*}}{\mathrm{d}{{\mathrm{y}}^{*}}^2}$$

(18)

• Once the equations have been rearranged and the new variables substituted, implement the Runge-Kutta method combined with the shooting technique. Define the system of equations as illustrated below:

  $$f_1={\displaystyle \frac{{\mathrm{dU}}_2}{{\mathrm{dy}}^{*}}=-1-{\hat{\mathrm{k}}}^{*}{\mathrm{W}}_2}$$

  (19)
  $$f_2={\displaystyle \frac{{\mathrm{dW}}_2}{{\mathrm{dy}}^{*}}={\displaystyle \frac{{\hat{\mathrm{k}}}^{*}.{\mathrm{U}}_2+2{\hat{\mathrm{k}}}^{*}.{\mathrm{W}}_1}{\zeta }}}$$

  (20)

• Establishing the boundary condition and providing the initial guess for the system

  $$\begin{array}{rl}& @y*=0,U_2=0,W_2=0\\ & @y*=1,W_1=-0.5U_{2,}{U}_1=0\end{array}$$

  (21)

And the initial guess for the unknowns at y* = 0 are:

$$\begin{array}{rl}& u*(0)=U_1(0)=0.6\\ & w*(0)=W_1(0)=0.07\end{array}$$

• Sixteen values of slopes (K) are calculated for each step, with four (4) values for each function, where the value of (K₁) at the beginning of interval, (K₂, K₃) at the midpoint of interval, and (K₄) at the end of interval.

• Next, calculate the values of U₁, U_2, W_1, W₂ after applying the following equation:

  $$U_{1,n}{}_{+}{}_1=U_{1,n+}{\displaystyle \frac{{\mathrm{K}}_{1,1}+2{\mathrm{K}}_{2,1}+2{\mathrm{K}}_{3,1+}{\mathrm{K}}_{4,1}}{6}}$$

  (22)
  $$U_{2,n}{}_{+}{}_1=U_{2,n+}{\displaystyle \frac{{\mathrm{K}}_{1,2}+2{\mathrm{K}}_{2,2}+2{\mathrm{K}}_{3,2+}{\mathrm{K}}_{4,2}}{6}}$$

  (23)
  $$W_{1,n}{}_{+}{}_1=W_{1,n+}{\displaystyle \frac{{\mathrm{K}}_{1,3}+2{\mathrm{K}}_{2,3}+2{\mathrm{K}}_{3,3+}{\mathrm{K}}_{4,3}}{6}}$$

  (24)
  $$W_{2,n}{}_{+}{}_1=W_{2,n+}{\displaystyle \frac{{\mathrm{K}}_{1,4}+2{\mathrm{K}}_{2,4}+2{\mathrm{K}}_{3,4+}{\mathrm{K}}_{4,4}}{6}}$$

  (25)

Continue iterating until the steps reach y* = 1 and compare the value with the boundary condition at this point if it is matching with it then the initial guess is correct otherwise adjust the initial guess value to achieve boundary condition.

In MATLAB, this method utilized the built-in routine bvp4c, designed for solving boundary value problems. MATLAB requires converting these equations into a system of coupled nonlinear differential equations, represented in the matrix format below:

$$\left[\begin{array}{c}{\mathrm{U}}_2\\ -1-{\hat{\mathrm{k}}}^{*}{\mathrm{W}}_2\\ {\mathrm{W}}_2\\ {\displaystyle \frac{{\hat{\mathrm{k}}}^{*}.{\mathrm{U}}_2+2{\hat{\mathrm{k}}}^{*}.{\mathrm{W}}_1}{\zeta }}\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right]$$

The boundary conditions required for the solutions are as follows:

$$\left[\begin{array}{c}{\mathrm{U}}_2(0)\\ {\mathrm{W}}_2(0)\\ {\mathrm{U}}_1(1)\\ {\mathrm{W}}_1(1)\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\\ -0.5{\mathrm{U}}_2(1)\end{array}\right]$$

Next, the system of equations is solved using MATLAB's bvp4c function. This function employs a three-stage Lobatto IIIa formula, ensuring accuracy and continuity within the specified interval. Finally, the outcomes are presented in the form of tables and graphs which illustrate in next section.

Numerical solution of constant surface heat flux in fully developed region for energy equation

In this section shown the solving of the linear second order partial differential of the dimensionless energy equation (11), according to the initial and boundary conditions in the equations (13) by finite difference method using MATLAB. The partial differential equation was discretized to transform it into an algebraic for each node by approximation by Taylor expansion as shown below:

$$\frac{\mathrm{\partial }\theta }{\mathrm{\partial }{\mathrm{x}}^{*}}=\frac{\left({\theta }_{\mathrm{i},\mathrm{j}}-{\theta }_{\mathrm{i}-1,\mathrm{j}}\right)}{\mathrm{\Delta }{\mathrm{x}}^{*}}$$

(26)

where the first order backward difference approximation.

$$\frac{{\mathrm{\partial }}^2\theta }{\mathrm{\partial }{{\mathrm{y}}^{*}}^2}=\left(\frac{{\theta }_{\mathrm{i},\mathrm{j}+1}-2{\theta }_{\mathrm{i},\mathrm{j}}+{\theta }_{\mathrm{i},\mathrm{j}-1}}{{(\mathrm{\Delta }{\mathrm{y}}^{*})}^2}\right)$$

(27)

while the second order is central difference approximation.

The finite difference approach to solve equation (3.20) is:

$${\mathrm{u}}^{*}{\displaystyle \frac{({\theta }_{\mathrm{i},\mathrm{j}}-{\theta }_{\mathrm{i}-1,\mathrm{j}})}{\mathrm{\Delta }{\mathrm{x}}^{*}}=\mathrm{Re}.\mathrm{Pr}\left({\displaystyle \frac{{\theta }_{\mathrm{i},\mathrm{j}+1}-2{\theta }_{\mathrm{i},\mathrm{j}}+{\theta }_{\mathrm{i},\mathrm{j}-1}}{{(\mathrm{\Delta }\mathrm{y}*)}^2}}\right)}$$

(28)

The finite difference method addresses the equation at specific, discrete points within the domain, known as grid points, as depicted in Figure 3. Here, Δy* indicates the step difference along the y-axis, while Δx* represents the step difference along the x-axis. The variables i and j correspond to the number of nodes along the x* and y* axes which were (100) for each, respectively, based on this number, Δx* and Δy* were set to (0.35, 0.01) respectively in order to obtain the solution of equation (11) in matrix form.

Assuming that Re and Pr are constant and u* is a function of y* after executive the Runge kutta code as shown previously, it was applied in the finite difference code. The discretized equation is used to form a linear system of algebraic equation in the following manner:

$$A.\theta =b$$

(29)

A matrix (N*N) was named (A) which is the N is the number of the nodes in the space and it was equal (10⁴ * 10⁴), and each row in this matrix represents the coefficient of (θ) for each node, while the matrix (b) shows the initial and boundary condition and its size was (10⁴ * 1) and considering that the initial guess value of (θ) for all nodes on y* were equal (1), finally the system which aforementioned above equation (29) was solved by using inverse matrix:

$$\theta =A^{-1}.b$$

(30)

The results of dimensionless temperature and dimensionless (mean and wall) temperature are presented in the form of graphs which illustrate in next section.

Numerical solution testing and results validation

To ensure the reliability of the numerical results, a comprehensive analysis of convergence, stability, and grid independence was conducted. The stability of the finite difference method (FDM) applied to the energy equation was assessed using the von Neumann stability criterion. As shown in Figure 2 (left), the numerical error decreases logarithmically as the time step size is refined, indicating a stable and convergent solution. Similarly, the convergence behavior of the Runge-Kutta fourth-order (RK4) method used for solving the momentum equations was analyzed. Figure 1 (right) demonstrates a consistent reduction in numerical error with decreasing step size, confirming the accuracy and stability of the RK4 method. Also, a grid independence study was performed by comparing numerical results at different grid resolutions (50 × 50, 100 × 100, and 200 × 200). Figure 3 presents the variation of the mean Nusselt number across these grid sizes. It was observed that the Nusselt number remains nearly constant between the medium (100 × 100) and fine (200 × 200) grids, with a relative difference of less than 1%. This confirms that the numerical solution is independent of further grid refinement, ensuring computational efficiency without compromising accuracy. Based on this analysis, the smaller grid was selected for the simulations, balancing precision and computational cost.

Figure 2.

Convergence and stability of numerical solution for both FDM (a) and RK4 (b) methods.

Figure 3.

Grid independence study for the numerical solution.

The verification of the current codes for the numerical solution, using the Runge-Kutta fourth order (RK4) method and Finite Difference Method (FDM) were examined using previous research papers to ensure that the solution was correct. Considered the cases from the previous studies along with their boundary conditions and implemented them in the current codes used in this study. The papers dealt with two parallel permeable porous walls and a two-dimensional micropolar fluid flow between them (Jalili et al., 2023), and (Mirgolbabaee et al., 2017). Considering the following differences from the current study:

The previous studies utilized the stream function $\psi (\mathrm{x},\mathrm{y})$ to express vertical and axial velocities ( $\mathrm{u},\mathrm{v}$ ).

The dimensionless stream function and temperature varied only along the vertical coordinate.

The applied boundary conditions were as follows:

$$\begin{array}{c}f={\mathrm{f}}^{\mathrm{\prime }}=g=0,\theta =\varphi =1,\mathrm{when}\xi =-1\\ f=\theta =\varphi =0,{\mathrm{f}}^{\mathrm{\prime }}=-1,g=1\mathrm{when}\xi =1\end{array}$$

(31)

After implementing the RK4 and FDM numerical schemes to solve the linear and angular momentum, as well as energy equations, in MATLAB, the obtained results for the dimensionless stream function (f) and dimensionless temperature (θ) as a function of the y-coordinate were found to be in good agreement with those from the benchmark studies, as illustrated in Figure 4 and Table 1. This verification confirms the accuracy and reliability of the present numerical approach. The results demonstrate a high degree of convergence between the current and previous values for both dimensionless stream function and temperature, with an error rate not exceeding 0.86% indicated a strong agreement and compatibility between the previous and current results, validating the accuracy and reliability of the codes used in this study.

Figure 4.

The results validation with Jalili et al. and Mirgolbabaee et al. for the dimensionless stream function f(y).

Table 1.

The results validation with Jalili et al. and Mirgolbabaee et al. for the dimensionless temperature.

ζ (Y-axis)	θ (ζ)
	Modified AGM (Jalili et al., 2023)	(Mirgolbabaee et al., 2017)	Present work
−1	1	1	1
−0.8	1.006453	1.012256	1.021
−0.6	0.999566	1.010345	1.026
−0.4	0.96911	0.983647	1.003
−0.2	0.908586	0.925411	0.9453
0	0.8155	0.832841	0.8511
0.2	0.691186	0.707079	0.7222
0.4	0.540247	0.552896	0.5642
0.6	0.369597	0.377888	0.3852
0.8	0.187126	0.190972	0.1945
1	0	0	0

Results and discussion

The findings of this study have important implications for practical thermal energy storage applications, particularly in solar energy systems and industrial heat exchangers. The ability of micropolar fluids to enhance heat transfer rates, especially in systems that require efficient thermal storage, could lead to improved performance in such applications. This research opens up possibilities for designing more efficient thermal energy storage systems using micropolar fluids.

The results were obtained for the fully developed region in the duct, with equations solved and presented in figures for transitional velocity, microrotation, Nusselt numbers (the ratio of convective to conductive heat transfer), and temperature profiles based on various parameters like coupling, spin gradient viscosity, Prandtl, and Reynolds numbers. Heat flux data were tailored for Amman using experimental data from the PVGIS website.

Dimensionless (transitional velocity and microrotation)

Figures 5(a) and 5(b) and Figures 6(a) and 6(b) illustrate the dimensionless translational velocity $u^{*}$ and microrotation $w^{*}$ across the dimensionless duct coordinate $y^{*}$ for various values of the coupling parameter k^̀* $(0.1,0.2,0.3,0.4)$ and spin-gradient viscosity parameter $\zeta (0.001,0.01,0.1,1)$ . The results show that increasing k^̀* enhances $u^{*}$ , as higher vortex viscosity increases the fluid's resistance to rotation, necessitating higher translational velocity to satisfy conservation of linear and angular momentum. At the duct center ( $y^{*}=0$ ), $w^{*}$ initially decreases with higher k^̀* but reverses beyond $y^{*}\approx 0.35$ , where velocity gradients are larger and angular shear stresses dominate, leading to increased microrotation. Conversely, lower $\zeta$ values result in higher translational velocity but initially lower microrotation at the center; $w^{*}$ then rises, intersecting other curves at several points before reaching the duct wall, where it attains its maximum value. These trends indicate a strong interplay between rotational viscosity effects and flow development, highlighting how micropolar parameters can be tuned to optimize momentum and energy transfer for enhanced thermal storage. Furthermore, the observed behavior emphasizes that both $u^{*}$ and $w^{*}$ are sensitive to microstructural fluid properties, which have direct implications for designing efficient TES systems using micropolar fluids.

Figure 5.

(a). Dimensionless transitional velocity with y* at different ${\hat{\mathrm{k}}}^{*}$ . (b). Dimensionless microrotation with y* at different ${\hat{\mathrm{k}}}^{*}$ .

Dimensionless temperature for different coupling parameters

Figure 7 (a) and 7 (b) illustrate the variation of dimensionless temperature $\theta$ across the dimensionless duct coordinate $y^{*}$ for different values of the coupling parameter k^̀* at two axial positions, $x^{*}=8.75$ and $x^{*}=$ 35. The results indicate that increasing k^̀* leads to a decrease in $\theta$ at each axial position. This trend can be attributed to the enhanced translational velocity associated with higher coupling parameters, which accelerates the fluid flow and reduces the residence time of the fluid within the heated region, thereby limiting its ability to absorb and store thermal energy effectively. Additionally, the interaction between linear and angular momentum in micropolar fluids redistributes energy within the fluid layers, contributing to the observed decrease in dimensionless temperature near the duct center. Along the flow direction, $\theta$ gradually increases as the fluid absorbs heat from the constant surface heat flux, demonstrating the combined effects of conduction and convection in the duct. The rate of temperature increase is more pronounced for lower k^̀* values, indicating that the interplay between rotational viscosity and velocity gradients strongly influences the thermal energy storage capability. These observations highlight the critical role of the coupling parameter in controlling both the flow structure and thermal transport within micropolar fluids, providing valuable guidance for optimizing duct geometries and micropolar fluid properties for enhanced solar thermal energy storage applications.

Figure 6.

(a). Dimensionless transitional velocity with y* at different ζ. (b). Dimensionless microrotation with y* at different ζ.

Dimensionless temperature for different spin gradient viscosity parameters

Figures 8 (a) and 8(b) illustrate the effect of the spin-gradient viscosity parameter $\zeta$ on the dimensionless temperature $\theta$ across the dimensionless duct coordinate $y^{*}$ at different axial positions $x^{*}$ . The results show that increasing $\zeta$ leads to higher and more uniform microrotation throughout most of the duct crosssection. This enhanced microrotation promotes greater mixing of fluid layers, which in turn improves convective heat transfer and stabilizes the temperature distribution. At the same time, the translational velocity is relatively lower for higher $\zeta$ , which increases the residence time of the fluid in the heated region, allowing more thermal energy to be absorbed and stored. Consequently, the overall thermal energy storage capacity of the micropolar fluid is enhanced. These findings demonstrate that $\zeta$ plays a dual role by simultaneously controlling rotational motion and optimizing heat retention, emphasizing its significance for designing efficient thermal energy storage systems using micropolar fluids.

Dimensionless temperature for different Reynolds number

The effect of Reynold number on dimensionless tempreture with y* for several points x* was represented in Figure 9 (a,b), the Reynold number is the ratio between inertia force and viscouss force so when increasing reynold number means that the viscouss force decreased so that the mean translational velocity will increase which may reduce the effectiveness capacity of energy stored in fluid layers.

Figure 7.

(a). Dimensionless temperature with y* at x* = 8.75 for different ${\hat{\mathrm{k}}}^{*}$ . (b). Dimensionless temperature with y* at x* = 35 for different ${\hat{\mathrm{k}}}^{*}$ .

Dimensionless temperature for different Prandtl number

Figures 10(a) and 10(b) depicts the effect of varying Prandtl numbers (Pr) on the dimensionless temperature $\theta$ at different axial positions $x^{*}$ . The results show that lower Prandtl numbers correspond to higher $\theta$ values across the duct. This behavior can be explained by the definition of the Prandtl number as the ratio of momentum diffusivity to thermal diffusivity. A lower $Pr$ indicates higher thermal diffusivity relative to momentum diffusivity, resulting in a thicker thermal boundary layer compared to the hydrodynamic boundary layer. Consequently, heat diffuses more rapidly into the fluid, increasing the overall temperature. In contrast, higher Prandtl numbers limit thermal diffusion, leading to thinner thermal boundary layers and lower fluid temperatures. These observations highlight the critical role of $Pr$ in governing energy transport and storage in micropolar fluids, offering insight into selecting suitable fluid properties for efficient thermal energy storage systems.

Dimensionless mean temperature

Figure 11 (a,b,c,d) depicted dimensionless mean temperature along the duct x* with selective values of coupling parameter, spring gradient viscosity parameter, Reynold number, and Prandtl number. It is clear from the previous figures when the Reynold number and Prandtl number rise up the dimensionless temperature decreased, because there are not enough time to capture more heat from solar source due to higher velocity and lower thermal conductivity respectively, thus the dimensionless mean temperature (average temperature) is going to be decreased on the other hand the impact of coupling parameter is inversely with mean temperature, but the spin gradient viscosity parameter is proportional with it, due to being coupled with velocity, and noticed that the dimensionless mean temperature increased linearly in case of constant surface heat flux.

Figure 8.

(a). Dimensionless temperature with y* at x* = 8.75 for different ζ. (b). Dimensionless temperature with y* at x* = 35 for different ζ.

Dimensionless wall temperature

Figures 12(a)-(d) illustrate the effect of the previously discussed parameters on the dimensionless wall temperature ${\theta }_w$ at $y^{*}=1$ along the axial coordinate $x^{*}$ . In the entrance region, the difference between the wall and mean temperatures ( $\mathrm{\Delta }\theta ={\theta }_w-{\theta }_m$ ) is relatively small, as the convective heat transfer coefficient h is high. This high heat transfer at the inlet requires only a small temperature difference to maintain the prescribed surface heat flux, resulting in the characteristic parabolic behavior of ${\theta }_w$ in the entrance region. As the flow develops downstream, the local heat transfer coefficient decreases and stabilizes, which increases the required driving force ( $\mathrm{\Delta }\theta$ ) to sustain the same surface heat flux. Consequently, the wall temperature rises almost linearly in the fully developed region, reflecting the steady-state balance between heat addition by the wall and the convective transport by the micropolar fluid. These results clearly demonstrate the influence of both flow development and micropolar parameters on wall heat transfer and thermal energy storage, providing insights for optimizing duct designs in solar thermal applications.

Local Nusselt number

Figure 13 (a, (b), (c), (d)) discussed the local Nusselt number along the duct for selective values of parameters ( ${\hat{\mathrm{k}}}^{*}$ , ζ, Re, Pr) where the local Nusselt number is calculated with respect to Reciprocal of the fraction for the difference between the wall and the entrance temperature, as shown in the graphs the local Nusselt number is approximately decays exponentially until it reaches zero because of decreasing convective heat transfer coefficients along the flow direction.

Figure 9.

(a). Dimensionless temperature with y* at x* = 8.75 for different Re. (b). Dimensionless temperature with y* at x* = 35 for different Re.

Mean Nusselt number

Figures 14(a)-(d) presents the variation of the mean Nusselt number ( $N{u}_m$ ) along the axial coordinate $x^{*}$ for different values of the coupling parameter (k^̀*), spin-gradient viscosity ( $\zeta$ ), Reynolds number ( $Re$ ), and Prandtl number $(Pr)$ . The mean Nusselt number is calculated based on the reciprocal of the dimensionless temperature difference between the wall and the fluid bulk. As expected, $N{u}_m$ is highest in the entrance region due to the strong thermal gradients and high convective heat transfer at the duct inlet. As the flow develops downstream, $N{u}_m$ gradually decays and stabilizes, reaching approximately $99\text{\%}$ of its fully developed value, $N{u}_m\approx 2.53$ , at $x^{*}=5.5$ , marking the onset of the fully developed thermal region. Increasing the coupling parameter, Reynolds number, or Prandtl number extends the entrance region, as stronger velocity gradients or higher adjective transport delay the establishment of a fully developed thermal profile. In contrast, higher values of the spin-gradient viscosity parameter $\zeta$ enhance rotational diffusion and mixing, which accelerates thermal equilibration and shortens the entrance length. These results emphasize how micropolar fluid parameters, flow inertia, and thermal properties collectively govern the development of convective heat transfer, providing crucial guidance for optimizing duct geometries and fluid characteristics to maximize thermal energy storage efficiency in solar energy applications.

Temperature variation during the day at different heat flux

In Figure 15 (a), (b), (c), (d) as shown below represented the temperature difference between outlet and inlet mean temperature with hourly time for several values of heat flux in Amman for three months (January, June, and September), the data which was taken from PVGIS website according to equation 4.26. Two values of the variables ( ${\hat{\mathrm{k}}}^{*}$ , ζ, Re, Pr) were taken for each month to demonstrate their effect on temperature differences, it was clear that the best heating storage in micropolar fluid occurred in September due to high solar irradiance, Additionally, the time that recorded the highest temperature variation for all months was at solar noon.it was observed that the highest value of variation temperature for almost parameters was achieved (60–70) Celsius in September, except for when the Reynolds number was set to 5, as it had the highest impact on temperature differences, reaching approximately 95 degrees Celsius in September, This indicates the optimal values for the parameters to achieve the maximum heat storage within the micropolar fluid was ( ${\hat{\mathrm{k}}}^{*}$ =0.4, ζ = 0.01,Re = 5, Pr = 5) were determined for September. The observed temperature variation of 95°C in September was calculated based on the fully developed temperature profiles along the length of the duct, considering the impact of the Reynolds number at that time of year. The high temperature variation observed during this period is a result of the combination of increased fluid velocity and the specific thermophysical properties of the micropolar fluid, leading to more effective thermal storage in the system. This provides valuable insight into the potential seasonal impacts on system performance.

Figure 10.

(a). Dimensionless temperature with y* at x* = 8.75 for different Pr. (b). Dimensionless temperature with y* at x* = 35 for different Pr.

Figure 11.

(a). Dimensionless mean temperature with x* for different ${\hat{\mathrm{k}}}^{*}$ . (b). Dimensionless mean temperature with x* for different ζ. (c). Dimensionless mean temperature with x* for different Re. (d). Dimensionless mean temperature with x* for different Rr.

Figure 12.

(a). Dimensionless wall temperature with x* for different ${\hat{\mathrm{k}}}^{*}$ . (b). Dimensionless wall temperature with x* for different ζ. (c). Dimensionless wall temperature with x* for different Re. (d). Dimensionless wall temperature with x* for different pr.

Figure 13.

(a). Local Nusselt number with x* for different ${\hat{\mathrm{k}}}^{*}$ . (b). Local Nusselt number with x* for different ζ. (c). Local Nusselt number with x* for different Re. (d). Local Nusselt number with x* for different Rr.

Figure 14.

(a). Mean Nusselt number with x* for different ${\hat{\mathrm{k}}}^{*}$ . (b). Mean Nusselt number with x* for different ζ. (c). Mean Nusselt number with x* for different Re. (d). Mean Nusselt number with x* for different Rr.

Figure 15.

(a). Variation temperature with during the day for different heat flux and ${\hat{\mathrm{k}}}^{*}$ . (b). Variation temperature with during the day for different heat flux and ζ. (c). Variation temperature with during the day for different heat flux and Re. (d). Variation temperature with during the day for different heat flux and Pr.

Regarding the temperature variation during the day, the September temperature variation of 95°C was calculated using experimental data from the PVGIS website, which provides solar radiation and temperature data specific to the region. This data was used to estimate the temperature profile during a typical day in September, which was then used in the current model to simulate the temperature distribution and assess the heat transfer behavior.

Conclusions

This study presents a detailed numerical investigation of micropolar fluid flow and heat transfer in a steadystate vacuum rectangular duct, emphasizing their potential for thermal energy storage in solar applications. The results show that micropolar fluid parameters strongly influence flow and thermal behavior. Increasing the coupling parameter raised the dimensionless transitional velocity, while initially decreasing microrotation before it increased, leading to a reduction in the dimensionless temperature. In contrast, higher spring gradient viscosity enhanced temperature retention, improving the fluid's energy storage capacity. Reynolds and Prandtl numbers were also influential, with lower values supporting higher temperature variation and more effective heat storage.

Analysis of Nusselt numbers revealed that the mean Nusselt number stabilized around 2.53 at an entrance length of ${\mathrm{x}}^{*}=5.5$ , while the local Nusselt number decreased exponentially along the duct. Temperature variations throughout the day indicated that September provided the highest potential for heat storage, with maximum temperature differences reaching approximately $95\circ \mathrm{C}$ at $\mathrm{Re}=5$ . By optimizing micropolar fluid parameters ( ${\mathrm{k}}^{*}=0.4,\zeta =0.01,\mathrm{Re}=5,\mathrm{Pr}=5$ ), the system achieved its best thermal storage performance, demonstrating the significant role of microrotation and spin-gradient effects in enhancing energy retention. The findings demonstrate how micropolar fluid parameters influence heat transfer and energy storage, providing new theoretical insights into microrotation and spin-gradient effects. Practically, they guide the design of optimized TES systems and suggest future studies on experimental validation and advanced micropolar fluid configurations.