A numerical investigation has demonstrated that regularizing the Herschel-Bulkley rheological model effectively resolves mathematical discontinuities at the yield stress threshold during boundary layer flow across a stretching sheet. Published by researchers Mohamed Ahmed Bayomy and Karem Mahmoud Ewis in Discover Applied Sciences, the study analyzes the steady, laminar, two-dimensional motion of incompressible micropolar fluids subjected to bioconvection, viscous dissipation, thermal radiation, and chemical reactions within a Darcy porous medium. Traditional non-Newtonian flow formulations involving Bingham plastics and Herschel-Bulkley materials frequently suffer from computational instabilities because vanishing shear rates cause apparent viscosity to approach infinity, generating abrupt boundaries between yielded and unyielded flow zones. By implementing the Papanastasiou regularization approach, the research team rendered the stress jump condition continuous across the entire computational domain, thereby stabilizing numerical convergence and providing detailed solutions for velocity, thermal distribution, mass diffusion, and microrotational behavior.
The underlying physical framework represents complex engineering and biological fluid systems where microscale rotational motions and suspended microstructure play an active dynamic role. Non-Newtonian media exhibiting yield stress followed by shear-dependent flow characteristics appear prominently across diverse industrial domains, including fresh concrete formulations, colloidal suspensions, paints, greases, slurries, starch pastes, and food processing systems. Furthermore, micropolar fluid mechanics account for microscopic phenomena arising from local structure and the intrinsic spinning inertia of fluid elements, making the formulation directly representative of animal blood, liquid crystals containing rigid molecules, and polymeric solutions with dilute suspensions of rigid macromolecules. In these specialized media, individual microscopic elements undergo distinct rotational motion about their centroids, supporting distributed body moments and couple stresses that classical Navier-Stokes models cannot adequately describe. Consequently, combining the Herschel-Bulkley-Papanastasiou framework with micropolar mechanics captures both macroscale non-Newtonian yield dynamics and microscale rotational momentum.
To mathematically formulate this multifaceted transport phenomenon, the investigators constructed a coupled system of nonlinear partial differential equations derived from the fundamental conservation laws governing momentum, microrotation, thermal energy, and mass species diffusion. The momentum formulation accounts for non-Newtonian consistency indices, yield stress parameters, microinertia density, Darcy resistance terms within the permeable substrate, and buoyancy forces driven by temperature and concentration gradients. Thermal energy conservation integrates viscous dissipation alongside Rosseland radiative heat flux approximations, Brownian motion, and thermophoretic particle transport. Similarly, the mass diffusion expression incorporates chemical reaction kinetics alongside coupled thermal-solutal diffusion gradients. By defining appropriate similarity transformations involving non-dimensional stream functions, similarity variables, and dimensionless profiles for axial velocity, microrotation, temperature, and nanoparticles volume fraction, the original partial differential equations were systematically transformed into a coupled set of nonlinear ordinary differential equations defined over a semi-infinite spatial domain.
The physical problem is bounded by realistic interfacial interactions that diverge from standard idealized conditions by incorporating multiple slip boundary conditions at the permeable stretching surface. Classical no-slip assumptions frequently fail when modeling particulate suspensions, emulsions, foams, and polymeric solutions flowing along solid boundaries. Consequently, the research model applies simultaneous slip conditions for axial fluid velocity, microrotation, thermal transport, and mass transfer at the wall where the similarity variable equals zero. The axial velocity condition relates surface slip linearly to the local wall shear stress, while the microrotation boundary parameter regulates how rigid microelements rotate adjacent to the solid interface, spanning configurations from zero wall microrotation to vanishing asymmetric stress tensor states. Meanwhile, thermal and solutal slip parameters govern interfacial heat and mass transfer resistances, interacting directly with wall suction or injection transpiration velocities and far-field ambient conditions where flow disturbances fully decay.
Because analytical solutions for highly coupled nonlinear non-Newtonian boundary layer problems are rarely attainable across general physical parameters, the researchers developed a robust dual-solver numerical framework. The primary solution methodology utilizes an implicit finite difference method applied to a linearized reduction of the transformed differential system. By introducing a reduction variable equating the first derivative of the stream function to an auxiliary velocity profile, the researchers decreased the highest differential order of the momentum equation from three to two, substantially improving numerical stability and iteration convergence. Finite difference approximations were then constructed across an interval divided into uniform mesh increments, with difference schemes selected specifically to yield diagonally dominant coefficient matrices. Iterative updates continued until the root mean square error between successive algorithmic iterations dropped below ten to the negative sixth power, using physically grounded Newtonian profile approximations as initial guesses.
To provide rigorous verification of the finite difference calculations, the authors simultaneously solved the governing equations using the adaptive, residual-controlled fifth-order bvp5c solver implemented within the MATLAB computational environment. Unlike fourth-order routines that evaluate errors indirectly, the fifth-order method directly controls true local truncation error during mesh adaptation. Across comprehensive parametric comparisons spanning yield stress indices, power-law exponents, and Eckert numbers, the relative discrepancy between the finite difference scheme and the fifth-order adaptive solver remained exceptionally small, frequently achieving error bounds between ten to the negative sixth and ten to the negative fourth power. Furthermore, validation against previously published exact, analytical, and numerical benchmarks for stretching sheets, shrinking boundaries, and Newtonian reference cases demonstrated excellent mutual agreement, confirming that the dual computational approach delivers dependable, highly stable predictions for boundary layer structures and interfacial engineering quantities.
Detailed parametric evaluations revealed that variations in the dimensionless yield stress parameter exert a profound influence over momentum, thermal, and concentration boundary layers. As the yield stress parameter increases, the normal velocity and axial velocity decrease within the near-wall region because enhanced yield resistance acts as an effective drag opposing boundary layer expansion. Near the solid boundary, axial velocity drops noticeably with increasing yield stress, although this trend reverses further out in the flow field where secondary transport mechanisms dominate. Conversely, temperature profiles rise with higher yield stress because amplified internal fluid friction and viscous dissipation generate additional thermal energy within the boundary layer. Meanwhile, concentration profiles decrease near the stretching boundary as higher yield stress suppresses mass convection. Notably, the computed boundary layer thickness for temperature measured approximately 2.5 on the similarity scale, compared to nearly 6.0 for velocity, demonstrating that thermal layers can remain significantly thinner than hydrodynamic layers under designated fluid parameters.
The power-law flow index similarly alters transport dynamics as the rheological model transitions between pseudoplastic, Newtonian, and dilatant fluid regimes. For pseudoplastic fluids characterized by power-law exponents below unity, flow resistance near the stretching surface is heightened relative to dilatant fluids possessing flow indices greater than one. Consequently, both normal and axial velocities increase as the power-law index increases, reflecting lower effective resistance in dilatant fluids near the moving wall. In contrast, temperature profiles decrease as the power-law index rises, because pseudoplastic behavior intensifies localized viscous dissipation and interfacial friction, thereby elevating wall-adjacent thermal levels. Concentration fields exhibit an increase near the sheet with larger flow indices due to enhanced convective mass transport. Across these parametric shifts, the hydrodynamic boundary layer thickness varied between 4.0 and 8.0, whereas the thermal and concentration boundary layer thicknesses remained confined near 3.0 and 4.0, respectively, confirming that fluid rheology dictates relative boundary layer growth.
The investigation systematically analyzed the impact of viscous dissipation by evaluating system responses across varying Eckert numbers. Higher Eckert numbers signify that kinetic energy dissipation within the moving fluid contributes substantially to internal thermal generation. Numerical findings demonstrated that increasing the Eckert number from zero to three produces noticeable increases in both normal and axial velocities as well as elevated temperature distributions throughout the boundary layer. This thermal enhancement directly expands the thermal boundary layer thickness toward a value of 4.0. Conversely, microrotation velocity and chemical species concentration decline near the stretching surface under rising Eckert numbers before reversing further into the outer flow field.
Transport characteristics were further modified by chemical and buoyant field variables, including the Schmidt number, concentration expansion parameter, and wall suction. Increasing the Schmidt number from 0.1 to 3.0 sharply reduced normal velocity, axial velocity, and chemical concentration profiles while thinning the concentration boundary layer, as momentum diffusion progressively outpaced mass species diffusion. When the Schmidt number is less than unity, mass diffuses faster than momentum, yielding a thicker concentration boundary layer than the velocity layer, a phenomenon that completely reverses when the Schmidt number exceeds one. Furthermore, augmenting the Grashof number and modified Grashof parameter enhanced natural buoyancy forces, driving upward increases in velocity components, skin friction, couple stress, and Sherwood numbers while lowering the Nusselt number. Similarly, increasing the wall velocity slip parameter dampened near-wall momentum transfer, reducing the skin friction factor, couple stress, and Sherwood number while enhancing the surface Nusselt number.
The authors established clear physical and mathematical boundaries for their computational framework, mapping stable convergence across broad operational domains encompassing power-law exponents up to 2.0, yield stress parameters up to 10, Schmidt numbers up to 2, and material parameters up to 30. The results definitively illustrate that boundary layer thicknesses are fundamentally parameter-dependent, debunking the assumption that hydrodynamic layers are invariably smaller or larger than thermal and solutal layers. While the regularized Papanastasiou method successfully handles the mathematical singularity of yield stress, the researchers identified vital directions for subsequent investigations. Future research priorities include developing exact analytical formulations using the differential transformation method, carrying out physical experimental validations to broaden operational parameter ranges, and extending the model to incorporate applied magnetic fields, full couple stress formulations, and non-Darcian Forchheimer porous medium effects for advanced manufacturing and biological transport designs.
Subject of Research: Numerical investigation of regularized Herschel-Bulkley micropolar boundary layer flow with bioconvection and multiple slip conditions
Article Title: Herschel Bulkley Papanastasiou flow of micropolar fluids under the effect of bioconvection over a stretching sheet through a porous medium with multiple slips
Article References: Bayomy, M. A., & Ewis, K. M. (2026). Herschel Bulkley Papanastasiou flow of micropolar fluids under the effect of bioconvection over a stretching sheet through a porous medium with multiple slips. Discover Applied Sciences, 8(11), Article 1348. https://doi.org/10.1007/s42452-026-09376-5
Image Credits: AI Generated
DOI: 10.1007/s42452-026-09376-5
Keywords: Herschel-Bulkley fluid, Papanastasiou regularization, micropolar fluids, bioconvection, stretching sheet, porous medium, multiple slips, finite difference method, bvp5c solver, boundary layer flow, skin friction factor, Nusselt number, Sherwood number, couple stress, viscous dissipation, non-Newtonian fluids, similarity transformation
News Source: Denise Maddox. (November 13, 2026). Numerical Model Resolves Yield Stress Discontinuities in Complex Micropolar Boundary Flows. Scienmag.
