Optimization of heat transfer in star-shaped thermal energy storage systems with NEPCM, and LTNE porous media

1 Introduction

The study of non-Newtonian flow and heat transfer by density driven convection within containers is considered a hot topic in recent decades. These studies become more intriguing when containers feature rotating shapes and the primary fluid consists of innovative substances that effectively enhance heat transfer, such as nanoparticle-enhanced phase change materials (NEPCM). Modern studies on non-Newtonian fluids have shown significant interest in micropolar fluids, especially when enhanced by the addition of certain nanoparticles. This field has become particularly intriguing due to the properties of micropolar fluids, which combine traditional fluid behaviors with micro-effects such as micro-rotation and rotational viscosity. When nanoparticles are added, the thermal and mechanical properties of the fluids are improved. These nanoparticle-enhanced micropolar fluids are utilized in various applications, including industrial coolants and solar energy systems, applications requiring viscosity control, drug delivery systems within the body, in bio-simulations, such as simulating blood flow containing nanoparticles and precision manufacturing processes and 3D printing due to their unique flow properties under varying speeds and directions. Several studies have explored micropolar nanofluids. For instance, Saranya et al. [1] discussed the flow between parallel disks using a micropolar nanofluid suspension. It was noted that the thickness of the boundary layer diminishes as the concentration of nanoparticles increases. The role of chemical reactions in the flow of a magnetized micropolar nanofluid along an extended surface was examined by Shamshuddin et al. [2]. The results revealed that the angular momentum improves as the magnetic parameter increases, while the flow activity decreases. In a related context, considerable attention has been directed toward the diverse behaviors and flow characteristics of micropolar fluids enhanced with various nanoparticles. Prominent investigations in such flow area include analyses of roles of velocity and temperature slip [3], flow over curved stretching surfaces [4], micropolar fluid behavior in the circular domain [5], non-homogeneous micropolar fluid model [6], the application of the Gauss–Lobatto IIIA approach to micropolar flow along an isothermal wedge [7], and the influence of cross-diffusion phenomena [8]. Furthermore, studies on internal flow using NEPCMs have been reported in Refs. [9], [10], [11]. Recent investigations by Ahlawat et al. [12], 13] have provided valuable insights into magneto-micropolar hybrid nanofluid convection and entropy generation within complex cavity geometries. Their studies on porous hexagonal enclosures and cavities with discrete heating configurations demonstrated that the interplay between magnetic field strength, microrotation effects, and nanoparticle concentration can substantially influence heat transfer efficiency and irreversibility. These works highlight the potential of micropolar hybrid nanofluids for optimizing thermal performance in nontraditional geometries, offering a foundation for exploring more intricate configurations such as star-shaped dynamic thermal systems, which represent a promising direction for advanced thermal management and energy storage applications.

Star-shaped thermal energy storage (TES) systems have recently gained attention due to their extended surface geometry, which significantly alters the heat transfer dynamics compared to conventional cylindrical or rectangular storage units. The presence of radial fins or pointed extensions increases the surface-area-to-volume ratio, promoting enhanced thermal interaction with the surrounding phase change or working fluid. When designed to maximize conduction pathways, star-shaped enclosures accelerate the melting or charging process by distributing heat more uniformly across the domain. Conversely, by adjusting the number, thickness, or angular distribution of the arms, the same configuration can be tailored to retard heat propagation, enabling controlled or delayed energy release during the solidification or discharging phase. Thus, the star-like geometry offers a tunable platform for either boosting thermal responsiveness in applications requiring rapid energy uptake or moderating heat release where thermal stability and longer-duration storage are desired.

The study of NEPCM-based suspensions within irregular flow domains has received considerable attention; however, the effects of rotating and stationary internal cylinders remain largely unexplored. For example, Afshar et al. [14] investigated the convection case due to the temperature differences in irregular chamber using a suspension-based on NEPCM. Their study also analyzed the system’s irreversibility, revealing that higher undulation numbers led to increased entropy generation. Nayak et al. [15] examined NEPCM suspensions in cylinders equipped with wavy baffles using a modified Fourier approach, and reported that larger baffle sizes reduced the velocity components. In a similar vein, Bouzidi et al. [16] used a mixture of NEPCM to estimate the rate of the heat transfer in quadrilateral containers, showing that a hybrid concentration of 5 % enhanced the heat transfer rate by approximately 13 % compared to the base water fluid. Hussain et al. [17] explored NEPCM applications in energy storage by studying recto-trapezoidal containers via finite element analysis, concluding that increases in either the Hartmann number or the Péclet number decreased the overall Nusselt number. Additional investigations on the properties and applications of NEPCM, are documented in refs. [18], [19], [20], [21], [22].

The study of fluid flow induced by the rotation of an inner cylinder within a concentric or eccentric annular domain is a fundamental problem in fluid mechanics with applications in engineering, geophysics, and biological systems. This phenomenon arises when an inner cylinder rotates within a stationary or co-rotating outer cylinder, generating a complex flow field influenced by centrifugal, viscous, and inertial forces. The resulting flow pattern is governed by the rotational speed of the inner cylinder, fluid properties (such as viscosity and density), and geometric parameters, including the radii of the inner and outer cylinders and the gap between them. The significance of studying this type of flow lies in its numerous industrial applications, such as in rotating machinery, journal bearings, and centrifuges, where precise control of fluid flow and heat transfer is essential. Additionally, understanding this flow behavior aids in improving mixing and enhancing convective heat transfer, particularly in heat exchangers. Moreover, it contributes to the modeling of natural processes, including planetary atmospheric dynamics and biological systems such as blood flow in cylindrical vessels [23], [24], [25], [26], [27]. Furthermore, additional related studies on this type of simulation can be found in refs. [28], [29], [30], [31].

Beyond the limitations of earlier works that relied on simplified NEPCM models or canonical geometries, the present study introduces a comprehensive LTNE micropolar framework in a highly irregular star-shaped TES system, supported by a Point-in-Polygon geometric treatment and a wide rotational parametric range, thereby delivering new physical insights unobtainable from previously published configurations. Therefore, the main objective of this study is to simulate the dynamic behavior resulting from the clockwise rotation of an internal cylinder placed within an irregular star-shaped enclosure. The external and internal boundaries are subjected to various dynamic and thermal conditions, making the analysis comprehensive. Several geometries are considered based on the undulation number of the outer contour. The thermal non-equilibrium between the fluid and solid phases is also taken into account. Furthermore, this study introduces a novel approach for analyzing flow and heat transfer within irregular geometries using the finite volume method (FVM). Also, other original contributions can be found in the following points:

1. The first numerical investigation of micropolar NEPCM flow inside a star-shaped porous TES enclosure under LTNE (Local Thermal Non-Equilibrium) conditions.

2. The integration of a dual-energy porous medium model with RBC-based suspension properties, which has not been previously analyzed in such irregular geometries.

3. The examination of a wide range of clockwise angular velocities for an inner rotating cylinder, which extends beyond the narrow ranges considered in the literature.

4. The incorporation of a Point-in-Polygon (PIP) algorithm to accurately resolve complex star-shaped irregular boundaries – a methodological contribution not previously applied in NEPCM-based TES systems.

5. The use of Response Surface Methodology (RSM) to optimize the Nusselt numbers of both fluid and solid phases, which provides quantitative design guidelines that are missing in earlier studies.

2 Problem assumptions

In Figure 1, the geometry of the computational model along with the imposed boundary conditions is presented. The following assumptions are considered:

1. The flow domain is bounded externally by a star-shaped boundary enclosing either a dynamic or stationary core.

2. The irregular outer boundary is defined by r = H + b × sin n ϑ , where H = 1and b = 0.2. The center of the inner cylinder is located at 0,0 with a radius of r ₁ = 0.2.

3. The inner cylinder rotates in the clockwise direction with an angular velocity ω.

4. Several dynamic and thermal cases are examined, namely:

   1. Internally heated, rotating core – the external edge is cold, while the internal core is hot and rotating.

   2. Internally heated, stationary core – the external edge is cold, while the internal obstacle is heated and stationary.

   3. Rotating core, externally heated – the external boundaries are hot, while the internal cylinder is cold and rotating.

   4. Stationary core, externally heated – the external boundaries are hot, while the internal cylinder is cold and stationary.

5. The flow is unsteady and laminar.

6. The host fluid, a micropolar liquid, is modeled as blood containing suspended RBCs, while the NEPCM consists of a nonadecane core and a polyurethane shell. The thermophysical properties of the NEPCM components are provided in Table 1.

7. The active domain is filled with an isotropic, permeable medium.

8. The dual-energy model is adopted, involving separate energy equations for the fluid and solid phases.

9. The porous matrix is composed of aluminum foam, and the entire domain is subjected to an inclined magnetic field.

Figure 1:

Geometry of the computational model with imposed conditions. (a) Internally heated, rotating core. (c) Rotating core, externally heated. (b) Internally heated, stationary core. (d) Stationary core, externally heated.

The eqs. used under the hypothesis are expressed as [32], 33]:

where, the velocity vector is v ̄ * = u * i ̄ + v * j ̄ , the fluid temperature is T s e * , the solid temperate is T m f * and B ̄ is the strength of the applied magnetic field. The subjected boundary conditions are:

1. Internally heated, rotating core

   1. On the star-shaped boundary:

      u * = v * = N * = 0 , T s e * = T m f * = T c * .

   2. On the cylindrical core boundary:

      u * = ω y * − c 2 * , V = − ω x * − c 1 * , N * = 0 , T s e * = T m f * = T h * .

2. Internally heated, stationary core

   1. On the star-shaped boundary:

      u * = v * = N * = 0 , T s e * = T m f * = T c *

   2. On the cylindrical core boundary:

      u * = v * = N * = 0 , T s e * = T m f * = T h *

3. Rotating core, externally heated

   1. On the star-shaped boundary:

      u * = v * = N * = 0 , T s e * = T m f * = T h *

   2. On the cylindrical core boundary:

      u * = ω y * − c 2 * , V = − ω x * − c 1 * , N * = 0 , T s e * = T m f * = T c *

4. Stationary core, externally heated

   1. On the star-shaped boundary:

      u * = v * = N * = 0 , T s e * = T m f * = T h *

   2. On the cylindrical core boundary:

      (6) u * = v * = N * = 0 , T s e * = T m f * = T c *

Here, it should be mentioned that these boundary conditions represent four operational configurations combining inner/outer heating with rotating or stationary cores. In the internally heated cases, the cylindrical core is maintained at the hot temperature while the star-shaped boundary is kept cold, with the core either rotating or stationary. In the externally heated cases, the heating is applied at the star-shaped boundary while the core remains cold, again with either a rotating or stationary cylinder. In all cases, no-slip and zero microrotation conditions are enforced on the solid boundaries, while the rotating-core conditions follow the standard solid-body velocity formulation. Also, the suspension’s properties are expressed as:

In eq. (7), T _f refers to a fixed fusion temperature, ΔT* is a slight variation in temperature leads to a phase change, ρ s e * represent the mixture density, ρ c * refers to the core density, ρ s * is the density of the shell, and ρ n p * and ρ b f * indicate the density of NEPCM nanoparticles, and the base fluid, respectively. The weight ratio between the core and shell, denoted by Ω, is about 0.447. Moreover, C p , n p , e ff * ​ refers to the effective heat capacity of the nanoparticles, Nc is variable thermal conductivity parameter, Nv is variable viscosity parameter and k s e * stands for the thermal. The previous system is converted from dimensional to nondimensional form using the following dimensionless variables.

Applying eq. (8) gives:

where, the Rayleigh number R a = g ∗ β b f ∗ T h ∗ − T c ∗ H * 3 ν b f ∗ α b f * , the Hartmann Number H a = B ∗ H ∗ σ b f ∗ μ b f * , the Darcy Number D a = K ∗ H * 2 , and the Prandtl number Pr = ν b f ∗ α b f ∗ = C p , b f ∗ μ b f ∗ k b f * . Additionally, H s = h ∗ k b f ∗ H * 2 is the fluid-solid heat transfer coefficient and Δ = k 1 μ b f * is the vortex viscosity parameter.

The converted boundary conditions are:

1. Internally heated, rotating core

   1. On the star-shaped boundary:

      U = V = N = 0 , T m f = T s e = 0 .

   2. On the cylindrical core boundary:

      U = ω Y − c 2 , V = − ω X − c 1 , N = 0 , T m f = T s e = 1 .

2. Internally heated, stationary core

   1. On the star-shaped boundary:

      U = V = N = 0 , T m f = T s e = 0

   2. On the cylindrical core boundary:

      U = V = N = 0 , T m f = T s e = 1 .

3. Rotating core, externally heated

   1. On the star-shaped boundary:

      U = V = N = 0 , T m f = T s e = 1

   2. On the cylindrical core boundary:

      U = ω Y − c 2 , V = − ω X − c 1 , N = 0 , T m f = T s e = 0

4. Stationary core, externally heated

   1. On the star-shaped boundary:

      U = V = N = 0 , T m f = T s e = 1

   2. On the cylindrical core boundary:

      (15) U = V = N = 0 , T m f = T s e = 0

In the previous equation, ω is the rotational velocity. Furthermore, the dimensionless heat capacity ratio Cr is expressed as:

where, λ = ρ s * ρ c * C p , c * + ρ s * ρ c * Ω C p , s * ρ b f * ρ s * C p , b f * + ρ b f * ρ c * Ω C p , b f * , S t e = ρ b f * C p , b f * ( T h * − T c * ) ( ρ s * + Ω ρ c * ) L s f * ρ s * ρ c * , T f = T f * T h * − T c * , Δ T = Δ T * T h * − T c * ,

In this study, heat transfer is evaluated for the fluid and solid phases individually, along with their overall contribution, and the associated coefficients are defined as:

The average heat transfer rates are:

The overall heat transfer rate is given as:

The corresponding formulations for the time-averaged heat transfer rates are:

where n′ is the outward normal vector, and L corresponds to the characteristic length of the rotating core.

Using the above values, the Prandtl number is Pr  = 22.86, and the suspension’s electrical conductivity is given by:

3 Solution technique

The use of the Finite Volume Method (FVM) in irregular geometries poses significant challenges for researchers, primarily due to the difficulty of accurately distinguishing points located inside the computational domain from those lying on its boundaries. To address this issue, a novel strategy has been adopted by combining the FVM with MATLAB’s inpolygon function. This function allows the classification of computational points into interior and boundary categories with high accuracy, thereby simplifying the setup of the numerical scheme in irregular domains.

The implementation begins by specifying the boundary of the irregular star-shaped cavity as two coordinate arrays (x, y). Subsequently, an array of candidate grid points (numerical nodes) is generated to represent the computational mesh. Using the inpolygon function, each point is tested to determine whether it lies inside the domain or on the boundary. This classification ensures proper enforcement of the governing equations within the interior and correct application of boundary conditions at the edges.

Once the mesh of the domain is obatined, the governing equations (9)– (14) are discretized via the finite volume approachs. A first-order upwind scheme is employed for convective derivatives, while the central difference scheme is applied to discretize the Laplacian terms. The resulting algebraic system is solved iteratively using the Successive Over-Relaxation (SOR) method.

The convergence criterion is set to the order of 10⁻⁶. To determine the appropriate grid, a grid-size independence test was performed at n = 3, Ra = 10⁵, Da = 10⁻², ω = −10and T _f  = 0.6 for the case of an internally heated and rotating core, as presented in Table 2. Noticeably, the grid size of 241 × 241 is suitable for lower Rayleigh numbers (Ra = 10³ − 10⁴), whereas higher values of Ra require grid sizes of 251 × 251 and 261 × 261. In this context, Figure 2 shows the structured mesh within the active domain 251 × 251 for n = 3 and b = 0.2. Furthermore, Table 3 presents a verification of time-step independence at n = 3, Ra = 10⁵, Da = 10⁻², ω = −10, and T _f  = 0.6 for the internally heated and rotating core case. It was found that a time step of Δt = 10⁻⁴is suitable for all computations.

Figure 2:

Structured mesh within the active domain (251 × 251) with n = 3 and b = 0.2.

The reliability of this approach is confirmed by strong agreement with the previous outcomes obtained by Kim et al. [37], which focused on natural convection in a regular domain containing a stationary cylinder, as can be observed from Table 4. Furthermore, Krane and Jessee [38] presented an experimental study on density-driven buoyancy convection within a square chamber. Here, a comparison with the results of Krane and Jessee [38] is performed and presented in Figure 3. Very good agreement is found between the results.

Figure 3:

Experimental data comparison (Krane and Jessee [38], at Ra = 1.89 × 10⁵).

However, while the FVM provides accurate solutions, direct numerical optimization in such complex domains would require an impractically large number of simulations. To overcome this limitation, the Response Surface Methodology (RSM) is employed as a complementary optimization tool.

4 Optimization of crucial physical variables

In this part, the heat transfer rates for the fluid and solid phases N u a v g , f , N u a v g , s together with the total heat transfer Qt are optimized using an effective technique that is response surface methodology (RSM). This approach is a statistical and mathematical technique designed to model and optimize responses influenced by multiple factors. In a related context, RSM builds a surrogate regression model that efficiently approximates the relationship between input variables and system responses instead of conducting an exhaustive simulation for all possible parameter combinations. The following three stages can describe the RSM technique:

1. Design of experiments (DoE): A statistical design such as the Central Composite Design (CCD) is used to select parameter combinations for numerical simulations, ensuring efficient exploration of the parameter space.

2. Surrogate model development: Simulation outputs from the FVM are fitted to a second-order polynomial model:

   R = β 0 + ∑ i = 1 k β i X i + ∑ i = 1 k β i i X i 2 + ∑ i < j β i j X i X j + ε

   where R is the predicted response, X _i are the design variables, β are regression coefficients, and ɛ represents residual error.

3. Optimization and validation: The surrogate model is employed to identify optimal conditions using desirability functions. Predictions are validated against additional FVM simulations not included in the regression fit to ensure robustness.

In this work, the primary input factors are:

1. Undulation number (n) of the star-shaped cavity that has range 3 ≤ n ≤ 7.

2. Heat transfer coefficient (Hs) with a range 0.1 ≤ Hs ≤ 100.

3. Fusion temperature (T _f ) of the NEPCM with a range 0.3 ≤ T _f  ≤ 0.9.

The referenced case of this process is set as: ω = −10, Ha = 10, Ra = 10⁴, Da = 10⁻² and Δ = 1. Also, the number of observations is 20 with error degrees of freedom 10. The following estimated coefficients for the second-order polynomial model for Nu_avg,f presented in Table 5.

Also, the estimated coefficients for the optimization of the average Nusselt number for the solid phase Nu_avg,s and the total heat transfer rate (Qt) are presented in Tables 6 and 7. In these tables:

1. Estimate represents the regression coefficients used in the model equation.

2. SE denotes the standard error, which quantifies the uncertainty of each estimated coefficient in which a smaller SE indicates higher reliability.

3. tStat refers to the t-statistic that it calculated as the ratio of the coefficient estimate to its standard error, indicating the strength of the effect.

4. pValue indicates the significance level that show whether the corresponding effect is statistically meaningful (typically, values below 0.05 are considered significant).

From these investigations, the following outcomes can be summarized:

1. The optimal values of n, Hs, and T _f for maximizing the fluid-phase heat transfer are 6.9781, 0.10014, and 0.38151, respectively, with a maximum predicted fluid-Nusselt number of 14.2699.

2. The optimal values of n, Hs, and T _f for optimizing the solid-phase heat transfer are 7.0000, 0.10002, and 0.3000, respectively, with a maximum predicted solid-Nusselt number of 5,189.86.

3. The optimal values of n, Hs, and T _f for enhancing the overall heat transfer are 7.0000, 0.10005, and 0.3000, respectively, with a maximum predicted total heat transfer rate of 1,329.06.

From the physical view, these outcomes confirm that increasing surface irregularity enhances the effective contact area between the heated wall and the phase change material. In this case, the melting is accelerated and a stronger thermal mixing is obtained. Also, the high predicted values of the Nusselt coefficients for both the solid and fluid phases disclosed that the NEPCM suspension effectively improves thermal conductivity while maintaining sufficient latent heat storage capacity.

From the practical view, this optimization implies that careful adjustment of wall temperature, rotational effects, and shape modulation can yield a balance between rapid energy charging and stable thermal response. Hence, the RSM findings validate the physical reasoning observed in prior simulations as well as offer a quantitative framework for designing compact and energy-efficient TES systems with complex geometries.

5 Results and discussion

This part presents a discussion of the obtained outcomes, highlighting the physical interpretations and underlying mechanisms governing the observed trends. The influences of the governing parameters on the flow field, thermal behavior, and transport characteristics are analyzed systematically. Here the key-factors are examined in wide ranges such as, the rotational velocity ω where −25 ≤ ω ≤ 0, the undulation number n where 3 ≤ n ≤ 9 with amplitude b = 0.2, the heat transfer coefficient Hs where 0.1 ≤ Hs ≤ 1,000, the Rayleigh number Ra where 10³ ≤ Ra ≤ 2 × 10⁵, the fusion temperature 0.1 ≤ T _f  ≤ 0.9, the Darcy number 10⁻⁴ ≤ Da ≤ 10⁻², the vortex viscosity parameter 0 ≤ Δ ≤ 6 and the dimensionless time where 0 < t ≤ 0.5. the referenced case for the illustrations are ϕ = 0.04, Nv = Nc = 3, Ha = 10, Hs = 10 and Φ = π/3.

In Figures 4 and 5, the variations of average Nusselt coefficient for the fluid phase Nu_avg,f and for the solid phase Nu_avg,s with the rotational velocity ω for different values of the undulation number n at Da = 10⁻², Ra = 10⁴, Δ = 1, T _f  = 0.6 are illustrated. It can be observed that the influence of ω on Nu_avg,f is more pronounced than on Nu_avg,s , primarily due to the dominance of conduction within the solid isothermal region. Increasing the magnitude of ω in the clockwise direction leads to a slight reduction in Nu_avg,f , indicating that stronger rotational motion tends to suppress convective circulation near the heated boundaries. In contrast, increasing the undulation number n enhances both Nu_avg,f and Nu_avg,s , as the wavy wall geometry promotes stronger mixing and increases the effective heat transfer area. From the physical view, the clockwise rotation of the cylinder causes that the induced flow alters the natural convection pattern. At higher rotational speeds, the centrifugal effect disrupts buoyancy-driven circulation, leading to weaker convective heat transfer in the fluid region and as a results a decrease in Nu_avg,f is given. Also, the number of surface waves on the boundary is raising as n is altered expanding the heat transfer surface area and generating localized vortices near the crests and troughs. These vortices enhance fluid mixing and promote stronger convection, leading to an overall rise in both fluid and solid Nusselt numbers.

Figure 4:

Variation of the average Nusselt number of the fluid phase Nu_avg,f with the rotational velocity (ω) for different undulation numbers (n) at Da = 10⁻², Ra ⁼ 10⁴, Δ = 1, T _f  = 0.6.

Figure 5:

Variation of the average Nusselt number of the solid phase Nu_avg,s with the rotational velocity (ω) for different undulation numbers (n) at Da = 10⁻², Ra = 10⁴, Δ = 1, T _f  = 0.6.

Figures 6 and 7 illustrate the variations of the mean Nusselt coefficient for the fluid phase (Nu_avg,f ) with time (t) for different values of the fluid–solid heat transfer coefficient (Hs) and the Rayleigh number (Ra) at ω = −10, n = 5, Da = 10⁻², Δ = 1, T _f  = 0.6. At early times, significant fluctuations in Nu_avg,f are observed, while as time progresses, Nu_avg,f gradually approaches a constant value, indicating the attainment of a steady-state condition. Increasing Hs enhances the thermal coupling between the fluid and solid phases, which facilitates energy exchange and reduces the local temperature gradients. Conversely, higher Ra values lead to stronger buoyancy-driven convection, resulting in higher Nu_avg,f and improved heat transfer performance. Physically, the sharp variations in Nu_avg,f at the lower t occur due to transient temperature adjustments between the fluid and solid regions while as the system evolves, thermal equilibrium is gradually established, leading to a stable heat transfer rate. Also, faster energy exchange across the interface is obtained at the higher Hs which causes smoothing temperature gradients and moderating variations in Nu_avg,f . Furthermore, a raising in Ra intensifies buoyancy forces, strengthens convective circulation, and enhances thermal transport, which manifests as higher Nu_avg,f values.

Figure 6:

Variation of the average Nusselt number of the fluid phase Nu_avg,f with the time parameter (t) for different heat transfer coefficient (Hs) at ω = −10, n = 5, Da = 10⁻², Ra = 10⁵, Δ = 1, T _f  = 0.6.

Figure 7:

Variation of the average Nusselt number of the fluid phase Nu_avg,f with the time parameter (t) for different Rayleigh number (Ra) at ω = −10, n = 5, Da = 10⁻², Ra = 10⁵, Δ = 1, T _f  = 0.6.

The variations of the average Nusselt number for the suspension, Nu_avg,f , together with the overall heat transfer, Qt, with respect to changes in the Darcy number (Da) for different values of the fusion temperature T _f are illustrated in Figures 8 and 9. The other parameters are fixed at ω = −10, n = 5, Ra = 10⁵. Remarkably, an increase in Da results in an enhancement of both Nu_avg,f and Qt. Physically, this behavior is attributed to the higher permeability of the medium, which intensifies the convective transport and consequently increases the temperature gradients between the solid and fluid phases. On the other hand, increasing the fusion temperature leads to a reduction in both Nu_avg,f and Qt. These profiles owing to that as the fusion temperature increases, the temperature difference between the hot and cold regions becomes smaller. This weakens the thermal driving force for convection, leading to lower heat transfer rates.

Figure 8:

Variation of the average Nusselt number of the fluid phase Nu_avg,f with the Darcy number (Da) for different fusion temperature (T _f ) at ω = −10, n = 5, Ra = 10⁵.

Figure 9:

Variation of the overall heat transfer rate Qt with the Darcy number (Da) for different fusion temperature (T _f ) at ω = −10, n = 5, Ra = 10⁵.

In Figures 10–12, 3D illustrations and contour plots for the fluid Nusselt number Nu_avg,f with the raising in the some governing parameters, namely, vortex viscosity coefficient Δ, the rotational velocity ω, the fusion temperature T _f , the Darcy number Da, heat transfer coefficient Hs, and the Rayleigh number Ra are illustrated. The value of the Rayleigh-Darcy number is fixed at Ra _D = Ra∗Da = 10. Here, it should be mentioned that a fixed cylinder with The lower values of Δ means less viscous suspension. Therefore the increase in Δ causes a reduction in values of Nu_avg,f due to the lower convective situation. Furthermore, the raising in Ra causes an increase in density-driven convection and hence the gradients of the temperature are enhanced while a reduction in Nu_avg,f is obtained as Hs is altered due the temperature differences eithin the flow area.

Figure 10:

3D surface and contour plot of average Nusselt coefficient for the fluid phase Nu_avg,f for the alterations of Δ and ω at n = 5, Da = 10⁻², Ra = 10³, Ha = 10, T _f  = 0.6.

Figure 11:

3D surface and contour plot of average Nusselt coefficient for the fluid phase Nu_avg,f for the alterations of T _f and Da at ω = 0, n = 5, Ra = 10³, Ha = 10.

Figure 12:

3D surface and contour plot of average Nusselt coefficient for the fluid phase Nu_avg,f for the alterations of Hs and Ra at ω = 0, n = 5, Da = 10⁻², Ha = 1, T _f  = 0.6.

Figure 13 show the vartiations of the streamlines (ψ), fluid and solid isotherms (T _se, T _mf), heat capacity ratio (Cr) and angular velocity (N) for the aformentioned cases at n = 5, b = 0.2, Ra = 10³, Da = 10⁻², Ha = Hs = 10, T _f  = 0.6, ω = −20. In the first column, a rotating and internally heated core is considered. This case reveals that the clockwise rotation of the inner cylinder causes the fluid to concentrate in the left region of the irregular shape, with a maximum stream function of 0.25752. Moreover, the isotherms for both the fluid and solid phases indicate that the most heated zone occurs around the internal surface. The heat capacity distribution shows an active melting/solidification zone around the cylinder, with a largely inactive region near the outer edge. The microrotation profiles follow the stream function pattern due to the velocity source terms, exhibiting a small active area on the right-hand side. In the second column, where a fixed and heated inner cylinder is considered, a symmetrical pattern is observed with 47.15 % lower activity compared to the previous case. The thermal boundary layers around the inner cylinder exhibit a slight enhancement compared to the rotating case. The case where the inner cylinder is rotated in the clockwise direction and the outer edge is heated is placed in the 3rd column. The flow concentrates in the right hand side with ψ max = 0.30066 and the micropolar angular velocity get similar behavior while the isotherms contours revealed a major active area compared to the inner heated case. Here the lower values of the average heat transfer coefficient are owing to that it is computed around the inner cold cylinder. The active zone regarding the melting/solidification is higher than that of inner case. In the fourth column, where the outer edge is heated and the inner cylinder remains stationary, the flow activity decreases, and the fluid motion becomes more symmetrical. The overall flow strength is reduced by 44.73 % compared to the previous case, while the temperature gradients are intensified.

Figure 13:

Vartiations of the streamlines (ψ), fluid and solid temperature (T _se, T _mf), heat capacity ratio (Cr) and angular velocity (N) for the different considered cases at n = 5, b = 0.2, Ra = 10³, Da = 10⁻², Ha = Hs = 10, T _f  = 0.6, ω = −20.

Using the case of an internal rotating and hot cylinder with a clockwise angular velocity of ω = −20, the variations of the streamlines (ψ), fluid isotherms (T _se), heat capacity ratio (Cr), and angular velocity (N) with respect to changes in the undulation number (n) at an amplitude of R a = 1 0 3 , T f = 0.6 , Φ = π 3 , D a = 1 0 − 2 are displayed in Figure 14. When the value of n is small, the geometry exhibits lower complexity, which promotes stronger convective motion and fluid mixing. In this case, the fluid temperature and heat capacity ratio show nearly uniform distributions. As n increases, the flow mixing is suppressed, leading to a weaker convective regime. Conversely, the isotherm contours reveal an increase in the thickness of the thermal boundary layers as n varies. In a related context, Figure 15 depicts the vartiations of the streamlines (ψ), fluid temperature (T _se), heat capacity ratio (Cr) and angular velocity (N) for the raising of ω at Ra = 10³, Da = 10⁻², n = 3, b = 0.2, T _f  = 0.6, Φ = π/3. In this case, the rotating cylinder is heated while the outer edge is cold. Noticeably, the features of ψ and N are concentrated mainly in the left region of the flow domain, indicating a significant enhancement in flow activity as the clockwise rotation increases. The value of ψ max rises from 0.17614 at ω = −5 to 0.26909 at ω = −20, whereas the temperature gradients become weaker. From the physical view, the fluid circulation within the domain is enhanced as the rotational motion is rising. This stronger rotation generates greater shear forces and induces more vigorous mixing in the fluid region adjacent to the cylinder. Consequently, the convective transport of momentum becomes more dominant, resulting in a higher magnitude of the stream function ( ψ max ) and more active flow structure. Furthermore, vartiations of the streamlines (ψ), fluid temperature (T _se), heat capacity ratio (Cr) and angular velocity (N) for the alterations of Ra at n = 5, ω = −10, Da = 10⁻³, T _f  = 0.6, Φ = π/3 are illsutrtaed in Figure 16. It should be mentioned here that the impact of clockwise rotation is more significant at lower values of Ra. However, as Ra increases, the density-driven convection becomes dominant compared to the fluid circulation induced by internal rotation. Various flow configurations, temperature fields, and heat capacity ratio distributions are observed as Ra increases. The flow activity is enhanced from ψ max = 0.056804 at Ra = 10³ to, ψ max = 14.6768 at Ra = 2 × 10⁵. Also, the heat transfer rate around the cylinder gets an augmentation by 5.83 % as Ra is varied from 10³ to 2 × 10⁵. Further, the melting/solidification process is enhanced as Ra increase. These behaviors are due to that the buoyancy forces arising from temperature-induced density differences become stronger. These buoyancy-driven motions enhance natural convection within the domain, leading to more vigorous fluid circulation and greater thermal energy transport.

Figure 14:

Vartiations of the streamlines (ψ), fluid temperature (T _se), heat capacity ratio (Cr) and angular velocity (N) for the alterations of n at  b = 0.2 , R a = 1 0 3 , T f = 0.6 , Φ = π 3 , D a = 1 0 − 2 , ω = − 20 (internally heated, rotating core).

Figure 15:

Vartiations of the streamlines (ψ), fluid temperature (T _se), heat capacity ratio (Cr) and angular velocity (N) for the raising of ω at Ra = 10³, Da = 10⁻², n = 3, b = 0.2, T _f  = 0.6, Φ = π/3 (internally heated, rotating core).

Figure 16:

Vartiations of the streamlines (ψ), fluid temperature (T _se), heat capacity ratio (Cr) and angular velocity (N) for the alterations of Ra at n = 5, ω = −10, Da = 10⁻³, T _f  = 0.6, Φ = π/3 (internally heated, rotating core).

6 Conclusions

A numerical simulation was performed to examine the flow dynamics and heat transfer characteristics inside a star-shaped flow domain filled with a dual-energy porous medium. The working suspension consists of a red blood cell (RBC) mixture containing NEPCMs with nonadecane cores and polyurethane shells, while aluminum foam serves as the porous matrix. The study focuses on the flow structures and thermal fields, and the heat transfer rates for the fluid and solid phases are evaluated at the inner rotating cylinder. The key parameters considered include the clockwise rotational velocity, fusion temperature, Darcy number, vortex viscosity parameter, and Rayleigh number. Comparisons between inner/outer heating and rotating/stationary inner cylinders were also performed. The following major outcomes can be highlighted:

1. Clockwise rotation of the inner cylinder significantly influences fluid Nusselt number (Nu_avg,f ​) but has a minor effect on the solid phase due to dominant conduction.

2. Stronger rotation enhances shear forces, promotes mixing, and increases ψ max .

3. Higher undulation numbers enhance both Nu_avg,f and Nu_avg,s ​ by increasing surface area and generating localized vortices.

4. Increasing Ra strengthens density-driven convection, intensifies flow circulation, and raises ψ max from 0.0568 to 14.6768.

5. Small n leads to uniform temperature and heat capacity distributions, while larger n thickens thermal boundary layers and suppresses convective mixing.

6. Heat transfer around the cylinder increases by ∼5.83 % as Ra rises from 10³ to 2 × 10⁵.

7. Higher Da enhances permeability, boosting convective transport and overall heat transfer (Qt).

8. Higher Hs improves thermal coupling between fluid and solid phases, leading to faster attainment of steady-state conditions and smoother temperature distributions.

9. Symmetric flow patterns emerge when the cylinder is stationary or the outer edge is heated.

10. Optimizing these parameters can enhance convective transport and improve thermal performance in dual-energy porous systems.

The present investigation provides valuable insights for developing efficient thermal energy storage and waste-heat recovery systems based on micropolar NEPCM suspensions in porous media. The optimized parameters identified herein can guide the design of compact, high-performance units for industrial applications such as power generation, metallurgy, and process heat recovery, enabling improved energy utilization and reduced environmental impact.

Future research can build upon this work by conducting experimental validation of the proposed numerical model under realistic thermal and flow conditions to confirm its scalability. Further studies should also explore different nanoparticle types, encapsulation materials, and porous matrix configurations to enhance the thermophysical stability and storage capacity of NEPCMs. In addition, integrating machine learning or optimization algorithms could help predict optimal operating conditions in real time for large-scale systems. Extending the model to three-dimensional, transient, and multi-phase flow scenarios would also provide deeper understanding and facilitate the design of next-generation thermal management and energy recovery technologies