CFD and Machine Learning-based Predictive Modeling of Natural Convection in Non-Newtonian Nano-encapsulated Phase Change Material within an Enclosure with a Corrugated Heated Cylinder

✅ Published in Applied Thermal Engineering 2025 (Elsevier, Q1, Impact Factor: 6.9)
📄 Paper: ScienceDirect · DOI: 10.1016/j.applthermaleng.2025.127240

Abstract

This study examines the thermal efficiency of power-law non-Newtonian nano-encapsulated phase change materials (NEPCM) in a square cavity with a corrugated heated cylinder under natural convection. A non-dimensional framework was developed using the Galerkin Finite Element Method (GFEM) via COMSOL Multiphysics, with Polyethylene Glycol (PEG) as the base fluid. Machine learning-based predictive modeling is integrated with CFD. The study numerically examines the effect of Rayleigh number (Ra = 10⁴–10⁶), Hartmann number (Ha = 0–90), power-law index (n = 0.7–1.4), Prandtl number (Pr = 200), Stefan number (Ste = 0.313), and fusion temperature (θ_f = 0.3–0.8). Five ML models — Decision Tree, Random Forest, KNN, XGBoost, and LightGBM — predict Nu̅ and Bejan number (Be̅). At Ha = 90, Nu̅ drops by 65.88% (shear-thinning, n = 0.7). Nu̅ is reduced by 17.35% as the fluid transitions from n = 0.7 to n = 1.4 (Ra = 10⁵). Random Forest achieves R² = 0.9959 (Nu̅) and R² = 0.9999 (Be̅).

Introduction

Nano-encapsulated Phase Change Materials (NEPCMs) are PCMs encapsulated within nanoparticles, enhancing thermal conductivity, stability, and compatibility with carrier fluids. They absorb, store, and release latent heat during phase transitions, making them ideal for thermal energy storage in buildings, solar collectors, electronics cooling, and biomedical applications.

Key challenges: Natural convection in enclosures with complex heated geometries (corrugated cylinders) combined with non-Newtonian rheology and phase-change behavior is highly nonlinear. No prior study had used ML models trained on GFEM simulation data for NEPCM natural convection in enclosures with corrugated heated cylinders.

🎯 Key Contributions

• Non-Newtonian NEPCM Natural Convection — Power-law viscosity model (n < 1 shear-thinning, n > 1 shear-thickening) applied to NEPCM suspension in PEG base fluid under MHD natural convection.
• Corrugated Cylinder Geometry — Sinusoidally corrugated heated cylinder (R=5 corrugations, amplitude B=0.5) induces complex vortex patterns enhancing heat transfer.
• Entropy Generation Analysis — Full entropy production from fluid friction (S_Ft), heat transfer (S_Tt), and magnetic field (S_Mt) computed; Bejan number (Be̅) tracks irreversibility ratio.
• Five ML Surrogate Models — DT, RF, KNN, XGBoost, LightGBM trained on 576 GFEM data points with SHAP and permutation feature importance analysis.
• Polynomial Regression Correlation — 3rd-degree polynomial correlation for Nu̅ developed with R² = 0.944.

Methodology

Physical Setup

Figure 1: Square enclosure (L = H) with corrugated heated cylinder at center. Left/right vertical walls: cold (T = T_c); top/bottom: adiabatic (∂T/∂y = 0); corrugated cylinder: heated (T = T_h). NEPCM particles have PCM core (liquid→solid phase change) + shell structure.

Corrugated wall parametric equations:

x = (p + B × cos(2π·R·q)) × sin(2π·q)
y = (p + B × cos(2π·R·q)) × cos(2π·q)

where: p = 0.2 (cylinder radius), R = 5 (corrugation number), B = 0.5 (wave amplitude), q ∈ [0, 1].

Mesh: 34,684 triangular elements (GFEM via COMSOL Multiphysics). Boundary conditions: u = v = 0 on all walls (no-slip).

NEPCM Bulk Properties

Mixture density: ρ_b = (1 − φ)ρ_f + φρ_p

Thermal conductivity: k_b/k_f = 1 + N_c·φ, where N_c = N_v = 3.0

Power-law viscosity: μ_b/μ_f = (1 + Nv·φ)|γ̇|^{n-1}

NEPCM core specific heat (with phase change): includes latent heat term via sine-smoothed fusion function f to ensure continuity at melting boundary (T_f ± ΔT_Mr/2).

Heat capacity ratio: Cr = Cp_b·ρ_b / (Cp_f·ρ_f) = (1−φ) + φλ + φ·f/(δ·Ste)

Stefan number: Ste = (ρCp)(T_w − T_c)(ρ_s + lρ_c) / (h_sf·ρ_c·ρ_s)

Dimensionless Parameters

📊 Grid Independence (Ra=10⁶, Pr=200, Ha=90, Ste=0.313)

Selected: 34,684 elements. Middle mesh shows negligible variation (<0.1%) vs finer mesh — optimal accuracy–cost balance.

Code Validation

• Qualitative (Fig. 2): Isotherms, streamlines, and heat capacity ratio compared vs Ghalambaz et al. (2019) for Ra=10⁵, Ste=0.313, θ_f=0.3, Pr=6.2 — excellent visual agreement.
• Quantitative (Fig. 3): Vertical velocity (v at y=0.5) and Nu̅ profiles compared vs Turan et al. (2011) for non-Newtonian power-law fluids (n=0.6–1.8; Ra=10³–10⁶; Pr=100) — adequate agreement confirmed.

📊 Results — Parametric Study

Effect of Rayleigh Number (Ra) on Nu̅

• Ra↑: Buoyancy forces dominate viscous forces → stronger vortices → wavy crowded isotherms → NEPCM melting zone closer to corrugated cylinder.
• At Ra=10⁴: 4 corner vortices, smooth isotherms (conduction-dominated).
• At Ra=10⁵: Bottom vortices split, flow more chaotic.
• At Ra=10⁶: Large-scale circulations suppress small eddies; isotherms become highly curved.

Nu̅ increases 350% as Ra rises from 10⁴ to 10⁶ at Ha=0.

Effect of Hartmann Number (Ha) on Nu̅ and Be̅

Higher Ha → stronger Lorentz force → suppresses convection → smoother isotherms, widening slender zone at cylinder top. Phase transition region narrows and moves away from cylinder.

Effect of Power-law Index (n) on Nu̅

Local Nu̅ reduced by 30.52% at arc length 0.2 as n rises from 0.7 to 1.4.

• n=0.7 (shear-thinning): 6 vortices, tightly packed streamlines near heated cylinder, highly curved isotherms → maximum Nu̅.
• n=1.0 (Newtonian): 4 stable vortices, uniform temperature distribution.
• n=1.4 (shear-thickening): Narrow vortices, smoother isotherms, diffusive transitions → minimum Nu̅.

Effect of Fusion Temperature (θ_f)

θ_f has no discernible effect on overall isotherm contours or streamlines. It shifts the NEPCM phase transition region: higher θ_f brings the melting region closer to the corrugated cylinder, enlarging it. At θ_f = 0.8, the melting region is most prominent at the enclosure upper side.

📊 Table 2 — Entropy Profile (Ra=10⁵, Pr=200, Ste=0.313)

Key entropy findings:

• S_Ft and S_Tt decrease as n rises; Be̅ increases by ~120% (Ha=0, θ_f=0.3) as n goes 0.7→1.4.
• For n=0.7, Ha=0→90 with θ_f=0.3: S_Ft drops 97.25%, S_Tt drops 68.83%, S_St drops 81.74%.
• S_Mt decreases for shear-thinning fluids (n<1) but increases for shear-thickening (n>1) as Ha rises.

📊 Table 4 — ML Model Performance Metrics (Test Set)

Random Forest best overall: lowest RMSE for Nu̅ (0.3845) and highest R² (0.9959). DT has lowest MAPE for Nu̅ (1.0522%). All models achieve R² > 0.97 for Nu̅; RF/DT/LightGBM all hit R² = 0.9999 for Be̅.

📊 Table 5 — CFD vs Random Forest Prediction Comparison

RF predictions closely align with CFD results across all parameter combinations.

📊 Feature Importance (Permutation & SHAP Analysis)

Consistent across all 5 models on both training and test sets:

• Ra drives both Nu̅ and Be̅ (dominant buoyancy effect)
• n strongly influences Nu̅ (viscosity changes convection intensity)
• Ha moderately affects both
• θ_f has minimal influence on either target — consistent across all models, SHAP, and permutation methods

📊 Polynomial Regression Correlation for Nu̅ (R² = 0.944)

Nu̅ = 2.478×10⁻⁶ + 0.0005·Ra + 4.964×10⁻⁸·n + 1.406×10⁻⁶·Ha + 2.598×10⁻⁸·θ_f
     − 3.636×10⁻⁹·Ra² − 0.0001·Ra·n − 7.534×10⁻⁷·Ra·Ha
     + 4.949×10⁻⁸·n² + 2.435×10⁻⁸·n·θ_f + 4.037×10⁻¹⁵·Ha²
     + 6.733×10⁻⁷·Ha·θ_f + 1.433×10⁻⁸·θ_f²
     + 3.228×10⁻¹⁵·Ra³ + 2.792×10⁻¹¹·Ra²·n + 3.774×10⁻¹³·Ra²·Ha
     + 2.381×10⁻⁵·Ra·n² + 2.50×10⁻⁷·Ra·n·Ha + 3.977×10⁻¹⁰·Ra·Ha²
     + 5.37×10⁻⁸·n³ + 2.468×10⁻⁸·n²·θ_f + 1.35×10⁻⁸·n·θ_f²
     + 3.62×10⁻⁷·Ha·θ_f² + 8.67×10⁻⁹·θ_f³

Key physics embedded: Ra positive effect (dominant), n negative cross-terms (viscous suppression), Ha negative (magnetic damping), θ_f negligible.

Key Conclusions

CFD findings:

• Nu̅ increases 350% as Ra rises from 10⁴ to 10⁶ at Ha=0 — buoyancy-driven convection dominates
• Nu̅ drops 65.88% as Ha rises from 0 to 90 (n=0.7) — Lorentz force suppresses convection
• Nu̅ reduced 17.35% as fluid shifts from n=0.7 to n=1.4 at Ra=10⁵ — shear-thickening lowers heat transfer
• Nu̅ falls 61.95% at Ra=10⁶ as n goes from 0.7 to 1.4
• Total entropy S_St decreases 81.74% as Ha rises (n=0.7, θ_f=0.3)
• Be̅ increases ~120% as n rises (Ha=0) — thermal irreversibility grows with shear-thickening
• Fusion temperature θ_f shifts NEPCM melting region toward the cylinder but has negligible effect on Nu̅

ML findings:

• Random Forest most reliable: R²=0.9959 (Nu̅), R²=0.9999 (Be̅), RMSE=0.3845 (Nu̅), 0.0036 (Be̅)
• Decision Tree best MAPE: 1.0522% (Nu̅), 1.1948% (Be̅)
• Ensemble methods (RF, XGBoost, LightGBM) substantially outperform KNN for complex nonlinear relationships
• Ra is universally the most important feature; θ_f is consistently negligible

📚 Citation

@article{chowdhury2025cfd,
  title={CFD and Machine learning-based predictive modeling of natural convection in non-Newtonian nano-encapsulated phase change material within an Enclosure with a corrugated heated cylinder},
  author={Chowdhury, Md. Mejbah Ullah and Ahad, Jawad Ibn and Molla, Md. Mamun and Nag, Preetom and Rahman, Azad},
  journal={Applied Thermal Engineering},
  volume={278},
  pages={127240},
  year={2025},
  publisher={Elsevier},
  doi={10.1016/j.applthermaleng.2025.127240}
}