Magnetohydrodynamic Stagnation-Point Flow of a Rivlin-Eriksen Nanofluid over a Convectively Heated Cylinder in an Anisotropic Porous Medium

1. Introduction

Magnetohydrodynamics (MHD), which concerns the interaction between electrically conducting fluids and magnetic fields, has attracted sustained research interest due to its extensive applications in engineering and applied sciences. Such applications include electromagnetic casting, cooling of nuclear reactors, MHD generators, plasma propulsion, geothermal energy extraction, astrophysical flows, liquid metal processing, and biomedical engineering. The presence of a magnetic field introduces Lorentz forces that significantly modify momentum transport, thermal behavior, and mass diffusion within the fluid, thereby offering an effective mechanism for controlling flow characteristics in industrial and technological processes (Chandrasekhar 1961; Cramer and Pai 1973; Davidson 2001; Sutton and Sherman 1965).

In recent years, the coupling of MHD effects with non-Newtonian fluid models has gained prominence, as many industrial fluids exhibit viscoelastic behavior that cannot be adequately described by classical Newtonian theory. This has motivated the development of advanced constitutive models capable of capturing both viscous and elastic effects under electromagnetic influence.

Among non-Newtonian fluids, the Rivlin–Ericksen model represents an important class of viscoelastic fluids that accounts for material memory and normal stress effects while retaining mathematical tractability. Rivlin–Ericksen fluids have been successfully applied to model polymer melts, lubricants, paints, and biological fluids (Bird et al. 1987; Rivlin and Ericksen 1955; Tanner 1988). When such fluids are combined with suspended nanoparticles, the resulting Rivlin–Ericksen nanofluid exhibits enhanced thermal conductivity and mass transport properties, making it particularly suitable for advanced heat and mass transfer applications (Buongiorno 2006; Das et al. 2007; Kuznetsov and Nield 2010; Nield and Kuznetsov 2011).

The integration of nanofluid dynamics into MHD viscoelastic flow analysis has opened new avenues for optimizing thermal systems, especially in the presence of thermal radiation, Joule heating, chemical reactions, and porous structures. The flow models incorporating these parameters are typically governed by nonlinear and higher-order partial differential equations, which present substantial analytical and computational complexities in their formulation, analysis, and solution. Stagnation-point flows over cylindrical and stretching surfaces constitute a fundamental class of boundary-layer problems with direct relevance to aerodynamic heating, extrusion processes, wire coating, cooling of cylindrical electronic components, and chemical reactors. The curvature of the surface introduces additional complexity by altering velocity gradients and thermal boundary-layer thicknesses (Ariel 2002; et al. 2014; Hiemenz 1911; Mahapatra and Gupta 2002).

In MHD environments, stagnation-point flows over cylindrical geometries are further influenced by induced magnetic fields, curvature parameters, and slip conditions, which are particularly relevant in microscale and porous systems. Understanding the combined effects of viscoelasticity, nanoparticle transport, magnetic induction, and chemical reaction is therefore essential for the accurate prediction and control of such flows.

Extensive studies have been conducted on Newtonian MHD boundary-layer flows over stretching and stagnation-point surfaces (Pop and Ingham 2010; Schlichting and Gersten 2000; Shercliff 1953). Subsequent investigations extended these models to non-Newtonian fluids, including second-grade, Maxwell, Oldroyd-B, and Rivlin–Ericksen fluids (Abbas and Hayat 2011; Fetecau and Fetecau 2005; Hayat and Abbas 2007; Hayat and Qasim 2012; Rajagopal 1995). The influence of thermal radiation, Joule heating, and viscous dissipation on MHD viscoelastic flows has been examined by several authors (Ali et al. 2011; El-Azeez 2018; Kumari and Nath 2009b; Makinde 2011).

Nanofluid models incorporating Brownian motion and thermophoresis were pioneered by Buongiorno (Buongiorno 2006), leading to numerous studies on MHD nanofluid heat and mass transfer (Ibrahim and Shankar 2017; Sheikholeslami and Ganji 2015; Wang and Mujumdar 2019; Xuan and Li 2000). Chemically reactive nanofluid flows have also received attention due to their relevance in catalytic reactors and biomedical transport (et al. 2001; Noghrehabadi and Pourrajab 2012; Rasheed and Shah 2021).

More recently, studies have addressed induced magnetic field effects and magnetic induction equations in MHD boundary-layer flows (Attia 2009; Davidson 2001; Sutton and Sherman 1965). Semi-analytical techniques such as the Homotopy Perturbation Method (HPM), Homotopy Analysis Method (HAM), Adomian Decomposition Method, and weighted residual approaches have been widely employed to tackle the resulting nonlinear systems (Finlayson 1972; He 1999; He 2000; Liao 2003).

Despite these advances, it is observed that no research has been done on the investigation of viscoelastic memory effects, induced magnetic field, stagnation point flow over cylindrical surfaces,nanofluid heat and mass transfer with chemical reaction, wall slip and suction/injection and the use of hybrid semi-analytical solution strategy such as Homotopy Perturbation Method and Weighted Residual approach. Therefore, the novelty of the present study is centered on the simultaneous consideration of these factors within a unified framework and the successful application of a hybrid Homotopy perturbation Method (HPM) and Galerkin Weighted residual technique to obtain approximate solutions that satisfy all the boundary conditions.

To achieve these objectives, a hybrid semi-analytical approach is employed by combining the Homotopy Perturbation Method (HPM) with the Galerkin weighted residual method. The HPM is used to decompose the highly nonlinear governing equations into a hierarchy of linear subproblems, while a physically motivated exponential trial function is constructed to satisfy the boundary conditions exactly. The remaining unknown coefficients are determined using a Galerkin weighted residual technique, ensuring convergence and near-wall accuracy.

The results demonstrate that magnetic field strength and viscoelastic parameters significantly suppress velocity profiles due to enhanced Lorentz and elastic resistance, while suction stabilizes the boundary-layer. Thermal and concentration fields are strongly influenced by Brownian motion, thermophoresis, radiation, and chemical reaction parameters. The combined methodology yields accurate, physically admissible solutions that provide clear insight into parameter interactions.

2. Formulation of the Problem

This section is premised on the physical concept of problem modeling, leading to governing Partial Differential Equations (PDEs) and the appropriate boundary conditions concerning Rivlin-Ericksen viscoelastic nanofluid flow, heat and mass transfer past an incompressible and laminar, non-inclined solid and hollow cylinders immersed in Rivlin-Eriksen electrically and thermally conducting nanofluid in the presence of chemical reaction of higher order and induced magnetic field. Thermal radiation, higher order heat source or sink, Joule heating and thermal energy dissipation are all considered. By means of extant symmilarity transfomations, the PDEs are reduced to ODEs. We made use of the usual cylindrical coordinate system with velocity $\mathbf{v}=\left(v=u_r,w=u_{\vartheta }=0,u=u_z\right)$ with the fluid located at the point $(r,\vartheta ,z)$ and the assumption of two dimensional flow ensures that $\left(\vartheta ,{\displaystyle \frac{\partial }{\partial \vartheta }}\right)=0$ for axially symmetric physical model. In the model, Prandtl boundary-layer theory is employed i.e., the boundary-layer thickness, $\delta$ has the same order of magnitude as radial coordinate component of the fluid position as well as the velocity component $\mathbf{v}$ as can be seen in Figure 1. The following idealizations are presumed for the formulations:

i

The axial axis is taken as z while the perpendicular to it is r,

ii

The porous medium is anisotropic in nature,

iii

The fluid and porous medium are in perpetual thermal equilibrium without surface slipperiness,

iv

All physical properties are constant except porous-medium permeability and fluid density linear variation in the buoyancy term only in the momentum equations (the so-called Oberbeck-Boussinesq approximation),

v

The cylinder wall is stationary,

vi

The flow is engendered by the stationary wall and/or ambient stream along axial axis $\left(z\right)$ with velocity $U_e\left(z\right)$ respectively,

vii

The nanofluid is everywhere below the boiling point in the flow regime,

viii

The influence of magnetic polarization are negligible while that of the electric field is significant,

ix

Both nanoparticles and base fluid are in thermal equilibrium,

x

The effects of Hall current and/or ion-slip are negligible.

Under the aforementioned assumptions the boundary layer equations of continuity, momentum, energy, nanoparticle volume fraction for Rivlin-Erikssen nanofluid model can be written as;

Continuity and magnetic solenoidal nature equations in reduced 2D differential forms

$${\displaystyle \frac{\partial v}{\partial r}}+{\displaystyle \frac{\partial ru}{\partial z}}=0,{\displaystyle \frac{\partial r{B}_r}{\partial r}}+{\displaystyle \frac{\partial r{B}_z}{\partial z}}=0$$

(1)

where $\mathbf{B}=(B_r,0,B_z)$ is the induced magnetic field.

Momentum Equation

$${\displaystyle \frac{\partial }{\partial t}}(u-u_e)+u{\displaystyle \frac{\partial u}{\partial z}}-u_e{\displaystyle \frac{\partial u_e}{\partial z}}+v{\displaystyle \frac{\partial u}{\partial r}}={\displaystyle \frac{1}{{\rho }_f}}\left[-{\displaystyle \frac{\partial p}{\partial z}}+\left(\mu +{\mu }^{\prime }{\displaystyle \frac{\partial }{\partial t}}\right){▿}^{2}u+F_z\right]$$

(2)

where the axial component of the resultant body force is

$F_z=\left[-({\rho }_p-{\rho }_f)(\varphi -{\varphi }_{\infty })+(1-{\varphi }_{\infty }){\rho }_f\left({\beta }_T(T-T_{\infty })\right)\right]gsin\alpha +{\displaystyle \frac{1}{{\mu }_m}}\left(B_r{\displaystyle \frac{\partial B_z}{\partial r}}+B_z{\displaystyle \frac{\partial B_z}{\partial z}}\right)$

The dimensional energy balance (see (Adeniyan 2016; Ali et al. 2011; El-Aziz and Afify 2018; Kumari and Nath 2009a; Kuznetsov and Nield 2013; Panigrahi and Panda 2021)) is posited as

$${\left(\rho C\right)}_f\left({\displaystyle \frac{\partial T}{\partial t}}+u{\displaystyle \frac{\partial T}{\partial z}}+v{\displaystyle \frac{\partial T}{\partial r}}\right)=\left[\begin{array}{cc}& k_f{▿}^{2}T-{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial }{\partial r}}\left(r{q}_r\right)+\mu \left\{1+{\displaystyle \frac{{\nu }^{\prime }}{\nu }}{\displaystyle \frac{\partial }{\partial t}}\right\}{\left({\displaystyle \frac{\partial u}{\partial r}}\right)}^2+Q_0{(T-T_{\infty })}^{\varpi }\hfill \\ \hfill & {\left(\rho C\right)}_p\left\{D_B▿\varphi ·▿T+{\displaystyle \frac{D_T}{T_{\infty }}}▿T·▿T\right\}+{\displaystyle \frac{1}{\sigma }}{J}^2\hfill \end{array}\right]$$

(3)

Nanoparticles volume fraction equation (see (Ibrahim and Ul Haq 2016; Rashed et al. 2021)):

$${\displaystyle \frac{\partial \varphi }{\partial t}}+u{\displaystyle \frac{\partial \varphi }{\partial z}}+v{\displaystyle \frac{\partial \varphi }{\partial r}}=D_B{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial }{\partial r}}\left(r{\displaystyle \frac{\partial \varphi }{\partial r}}\right)+{\displaystyle \frac{D_T}{T_{\infty }}}{\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial }{\partial r}}\left(r{\displaystyle \frac{\partial T}{\partial r}}\right)-{\kappa }_r{(\varphi -{\varphi }_{\infty })}^m$$

(4)

Magnetic induction transport equation without Hall current and ion-slip effects, the magnetic induction balance equation and application of boundary layer in $2^D$ takes the form

$${\displaystyle \frac{\partial B_z}{\partial t}}+u{\displaystyle \frac{\partial B_z}{\partial z}}+v{\displaystyle \frac{\partial B_z}{\partial r}}=B_z{\displaystyle \frac{\partial u}{\partial z}}+B_r{\displaystyle \frac{\partial u}{\partial r}}+{\displaystyle \frac{{\eta }_m}{r}}{\displaystyle \frac{\partial }{\partial r}}\left(r{\displaystyle \frac{\partial B_z}{\partial r}}\right)$$

(5)

The Initial and boundary conditions are ((Attia 2009; Kumari and Nath 2009a))

$$t\le 0:u=v=0,T=T_{\infty },\mathbf{B}=\{B_e(t,z),0\},\varphi ={\varphi }_{\infty }\mathrm{for}\mathrm{all}(r,\vartheta ,z)\in {\mathbb{R}}^3,$$

(6)

where $B_z=B_e(t,z)=B_0\left({\displaystyle \frac{z}{\ell \left|1-\lambda t\right|}}\right),(T_f-T_{\infty })=T_0{z}^2,T_0={\displaystyle \frac{b_1}{1-\lambda t}}$. is constant

$$t>0:\left\{\begin{array}{cc}\hfill u& =U_e,v=V_w(t,z),-k{\displaystyle \frac{\partial T}{\partial r}}=h_s(T_f-T_w),D_B{\displaystyle \frac{\partial \varphi }{\partial r}}+{\displaystyle \frac{D_T}{T_{\infty }}}{\displaystyle \frac{\partial T}{\partial r}}=0,\hfill \\ \hfill \varphi & ={\varphi }_w(t,z)={\varphi }_{\infty }+{\displaystyle \frac{b_0}{[1-\lambda t]}}{z}^2,B_z={\displaystyle \frac{\partial B_x}{\partial r}}=0\mathrm{on}r=a,\hfill \\ \hfill u_e(t,z)& =U_e\left({\displaystyle \frac{z}{\ell [1-\lambda t]}}\right),{\displaystyle \frac{\partial u}{\partial r}}\to 0,T\to T_{\infty },B_z\to B_e(t,z)=B_0\left({\displaystyle \frac{z}{\ell [1-\lambda t]}}\right),\varphi \to {\varphi }_{\infty }\mathrm{as}r\to \infty .\hfill \end{array}\right.$$

(7)

Introducing the following similarity transformations;

$$\left\{\begin{array}{cc}\hfill \eta & ={\displaystyle \frac{r^2-a^2}{2a\lambda }}\sqrt{{\displaystyle \frac{U_e}{Vz}}}+\epsilon ,u_z\equiv u=U_w{f}^{\prime }\left(\eta \right),\theta ={\displaystyle \frac{T-T_{\infty }}{T_f-T_{\infty }}},\phi ={\displaystyle \frac{\varphi -{\varphi }_{\infty }}{{\varphi }_f-{\varphi }_{\infty }}},u_r=v=-{\displaystyle \frac{a}{r}}\sqrt{{\displaystyle \frac{\nu U_e}{z}}}f\left(\eta \right),\hfill \\ \hfill B_z& =B_0\left({\displaystyle \frac{z}{a}}\right)h^{\prime }\left(\eta \right),B_r=-{\displaystyle \frac{B_0}{r}}\sqrt{{\displaystyle \frac{Vz}{U_w}}}h\left(\eta \right),a\ne 0\&U_w=U\left({\displaystyle \frac{z}{\ell [1-\lambda t]}}\right).\hfill \end{array}\right.$$

(8)

The similarity transformation given in Equation (8) above, derived using Lie group analysis see (Adeniyan 2016; Adigun et al. 2021) satisfies the conservation of mass flow (continuity equation) and the soleinoidal nature law of magnetic field. The dimensionless momentum, energy, chemically reactive nanoparticle balance, and the magnetic induction balance equations are derived as

$$\begin{array}{cc}& (1+2\gamma \eta )f^{\prime \prime \prime }+2\gamma f^{\prime \prime }+f{f}^{\prime \prime }-f^{\prime 2}-A\left(f^{\prime }+{\displaystyle \frac{1}{2}}\eta f^{\prime \prime }\right)+Ri(\theta -N_r\phi )+Q{Z}^2(h^{\prime 2}-h{h}^{\prime \prime }-1)\hfill \\ \hfill & A{\beta }_p\left({\displaystyle \frac{1}{2}}\eta (1+2\gamma \eta )f^{\prime \prime \prime \prime }+(2+6\gamma \eta )f^{\prime \prime \prime }+4\gamma f^{\prime \prime }\right)+A\Omega -{\Omega }^2=0\hfill \end{array}$$

(9)

$$\begin{array}{cc}& {\displaystyle \frac{1+R_d}{P_r}}\left((1+2\gamma \eta ){\theta }^{\prime \prime }+2\gamma {\theta }^{\prime }\right)+E_c(1+2\gamma \eta )f^{\prime \prime 2}+f{\theta }^{\prime }-n{f}^{\prime }\theta +G{\theta }^{\varpi }\hfill \\ \hfill & +(1+2\gamma \eta )\left(N_b{\phi }^{\prime }{\theta }^{\prime }+N_t{\theta }^{\prime 2}\right)+Q\Lambda E_c{h}^{\prime \prime 2}=0\hfill \end{array}$$

(10)

$${\displaystyle \frac{1}{S_c}}\left((1+2\gamma \eta ){\phi }^{\prime \prime }+2\gamma {\phi }^{\prime }\right)+{\displaystyle \frac{1}{L_e{P}_r}}{\displaystyle \frac{N_t}{N_b}}\left((1+2\gamma \eta ){\theta }^{\prime \prime }+2\gamma {\theta }^{\prime }\right)-{\kappa }^2{\phi }^m+f{\phi }^{\prime }-n{f}^{\prime }\phi =0$$

(11)

$$\Lambda \left((1+2\gamma \eta )h^{\prime \prime \prime }+2\gamma h^{\prime \prime }\right)-h{f}^{\prime \prime }+f{h}^{\prime \prime }-A_t\left(h^{\prime }+{\displaystyle \frac{1}{2}}\eta h^{\prime \prime }\right)=0$$

(12)

Their corresponding boundary conditions are

$$\begin{array}{cc}& f\left(ϵ\right)=-f_w,f^{\prime }\left(ϵ\right)=\lambda f^{\prime \prime }\left(ϵ\right),{\theta }^{\prime }\left(ϵ\right)=-Bi(1-\theta \left(ϵ\right)),N_b{\phi }^{\prime }\left(ϵ\right)+N_t{\theta }^{\prime }\left(ϵ\right),h\left(ϵ\right)=h^{\prime \prime }\left(ϵ\right)=0\mathrm{at}\eta =ϵ\hfill \\ \hfill & f^{\prime }\left(\eta \right)\to 1,f^{\prime \prime }\left(\eta \right)\to 0,\theta \left(\eta \right)=\phi \left(\eta \right)\to 0,h^{\prime }\left(\eta \right)\to 1,\mathrm{as}\eta \to \infty \hfill \end{array}$$

(13)

where $ϵ\ll 1$ and $\lambda ={\displaystyle \frac{3}{2}}\varsigma {\beta }_p$, $N_r={\displaystyle \frac{({\rho }_p-{\rho }_f)}{{\rho }_f(1-{\varphi }_{\infty }){\beta }_T(T_f-T_{\infty })}}$ is the buoyancy ratio, $Ri=(1-{\varphi }_{\infty }){\beta }_T(T_w-T_{\infty }){\displaystyle \frac{gz}{U_w^2}}$ is the modified local Richardson number, $A_t={\displaystyle \frac{\lambda a}{U|1-\lambda t|}}$ is the local transient (unsteadiness) parameter with the constant $l=a$, $Q={\displaystyle \frac{B_0^2}{\rho {\mu }_m{U}^2}}$ is the Chandrasekhar number, $\Lambda ={\displaystyle \frac{(1-\lambda t){\mu }_m}{lU}}$ is the local inverse modified magnetic Prandtl number, $\gamma ={\displaystyle \frac{1}{a}}\sqrt{{\displaystyle \frac{\nu z}{U_w}}}$ is the curvature parameter, $Ec={\displaystyle \frac{U_w^2}{C_f\Delta T}}$ is the Eckert number, $Pr={\displaystyle \frac{\nu }{{\alpha }_f}}$ is the (Thermal) Prandtl number, ${\alpha }_f$ is the thermal diffusivity, $R_d={\displaystyle \frac{16{\sigma }^{*}}{3k_f{k}^{*}}}{T}_{\infty }^3$ is the Rosseland radiation parameter, $N_b={\tau }_0{\displaystyle \frac{D_B\Delta \varphi }{\nu }},N_t={\tau }_0{\displaystyle \frac{D_T\Delta T}{\nu T_{\infty }}}$ are the Brownian and thermophoretic parameters, $Le={\displaystyle \frac{{\alpha }_f}{D_B}}$ is the modified Lewis number, ${\kappa }^2={\displaystyle \frac{z{\kappa }_r}{U_w}}{(\Delta \varphi )}^{m-1}$ is the chemical reaction parameter of the nanofluid, $\varsigma =N\sqrt{{\displaystyle \frac{U_w}{\nu U_w}}}$ is the slip parameter, $f_w=V_w\sqrt{{\displaystyle \frac{z}{\nu U_w}}}$ is the surface mass flux (suction/injection) at the wall surface. $\Omega ={\displaystyle \frac{U_e}{U}}$ is the velocity ratio of the free stream.

2.1. Parameters of Engineering Interest

The important physical quantities of interest are the local skin-friction coefficient $C_f$, the local Nusselt number $N{u}_z$, and the local Sherwood number $NS{h}_z$ are defined as

$$C_f={\displaystyle \frac{{\tau }_w}{{\rho }_n{U}_{\infty }^2}},N{u}_z={\displaystyle \frac{z{q}_w}{k_f\Delta T}},\mathrm{and}NS{h}_z={\displaystyle \frac{z{j}_w}{D_B\Delta \varphi }}$$

(14)

where ${\tau }_w$ is the skin-friction or the shear stress along the stretching surface and $q_w$ is the wall heat flux incorporating the impact of thermal radiation and $j_w$ is the nanofluid volume-fraction (flux) transfer at the wall of the cylinder are given respectively as follows

$$\left\{\begin{array}{cc}& {\tau }_w={\displaystyle \frac{{\mu }^{\prime }}{{\rho }_f}}{\displaystyle \frac{\partial }{\partial t}}{\left({\displaystyle \frac{\partial u}{\partial r}}\right)}_{r=\epsilon }\hfill \\ \hfill & q_w={\displaystyle \frac{z}{\Delta T}}{k}_f\left(1+{\displaystyle \frac{16{\sigma }^{*}{T}_{\infty }^3}{3k_f{k}^{*}}}\right){\left({\displaystyle \frac{\partial T}{\partial r}}\right)}_{r=\epsilon },\mathrm{and}\hfill \\ \hfill & j_w=D_B{\left({\displaystyle \frac{\partial \varphi }{\partial r}}\right)}_{r=\epsilon }\hfill \end{array}\right.$$

(15)

substituting the non-dimensionless parameters, we obtain

$$\left\{\begin{array}{cc}& R{e}_z^{{\displaystyle \frac{1}{2}}}{C}_f={\displaystyle \frac{3}{2}}A{\beta }_p{f}^{\prime \prime }\left(\epsilon \right)\hfill \\ \hfill & R{e}_z^{-{\displaystyle \frac{1}{2}}}N{u}_z=(1+R_d)Bi(1-\theta \left(\epsilon \right))\hfill \\ \hfill & R{e}_z^{-{\displaystyle \frac{1}{2}}}NS{h}_z={\phi }^{\prime }\left(\epsilon \right)\hfill \end{array}\right.$$

(16)

where $R{e}_z={\displaystyle \frac{U_{w}z}{\nu }}$ is the local Reynolds number.

3. Method of Solution

The transformed momentum equation obtained in the preceding section constitutes a highly nonlinear boundary-value problem due to the combined effects of viscoelasticity, magnetic field, wall slip, and suction/injection. The presence of higher-order derivatives arising from the Rivlin–Ericksen fluid model further precludes the possibility of obtaining an exact closed-form solution. Consequently, a semi-analytical approach is adopted to obtain an approximate solution that preserves the essential physics of the flow while satisfying all boundary conditions exactly. The solution strategy involves the construction of a physically motivated exponential trial function, followed by the determination of the unknown coefficients using the imposed boundary conditions and a weighted residual technique. The details of the assumed solution structure and the procedure for evaluating the unknown constants are presented in the subsequent subsections.

3.1. Applications of Homotopy Perturbation Method to the Problem

$$\begin{array}{cc}& {\mathcal{H}}_1(\eta ,f,\theta ,\phi ,h;p)\equiv (1-p){\displaystyle \frac{1}{2}}A{\beta }_p\eta (1+2\gamma \eta )f^{\prime v}\hfill \\ \hfill & +p\left[\begin{array}{cc}& 1+2\gamma \eta )f^{\prime \prime \prime }+2\gamma f^{\prime \prime }+f{f}^{\prime \prime }-f^{\prime 2}-A\left(f^{\prime }+{\displaystyle \frac{1}{2}}\eta f^{\prime \prime }\right)+Ri(\theta -N_r\phi )\hfill \\ \hfill & A{\beta }_p\left({\displaystyle \frac{1}{2}}\eta (1+2\gamma \eta )f^{\prime v}+(2+6\gamma \eta )f^{\prime \prime \prime }+4\gamma f^{\prime \prime }\right)+Q{Z}^2(h^{\prime 2}-h{h}^{\prime \prime }-1)+A\Omega -{\Omega }^2\hfill \end{array}\right]\hfill \end{array}$$

(17)

$$\begin{array}{cc}& {\mathcal{H}}_2(\eta ,f,\theta ,\phi ,h;p)\equiv (1-p){\displaystyle \frac{1+R_d}{P_r}}{\theta }^{\prime \prime }\hfill \\ \hfill & +p\left[\begin{array}{cc}& {\displaystyle \frac{1+R_d}{P_r}}\left((1+2\gamma \eta ){\theta }^{\prime \prime }+2\gamma {\theta }^{\prime }\right)+E_c(1+2\gamma \eta )f^{\prime \prime 2}+f{\theta }^{\prime }\hfill \\ \hfill & -n{f}^{\prime }\theta +G{\theta }^{\varpi }+(1+2\gamma \eta )\left(N_b{\phi }^{\prime }{\theta }^{\prime }+N_t{\theta }^{\prime 2}\right)+Q\Lambda E_c{h}^{\prime \prime 2}\hfill \end{array}\right]\hfill \end{array}$$

(18)

$${\mathcal{H}}_3(\eta ,f,\theta ,\phi ,h;p)\equiv (1-p){\displaystyle \frac{1}{Sc}}{\phi }^{\prime \prime }+p\left[\begin{array}{cc}& {\displaystyle \frac{1}{S_c}}\left((1+2\gamma \eta ){\phi }^{\prime \prime }+2\gamma {\phi }^{\prime }\right)+{\displaystyle \frac{1}{L_e{P}_r}}{\displaystyle \frac{N_t}{N_b}}\left((1+2\gamma \eta ){\theta }^{\prime \prime }+2\gamma {\theta }^{\prime }\right)\hfill \\ \hfill & -{\kappa }^2{\phi }^m+f{\phi }^{\prime }-n{f}^{\prime }\phi \hfill \end{array}\right]$$

(19)

$${\mathcal{H}}_4(\eta ,f,\theta ,\phi ,h;p)\equiv (1-p)\Lambda h^{\prime \prime \prime }+p\left[\Lambda \left((1+2\gamma \eta )h^{\prime \prime \prime }+2\gamma h^{\prime \prime }\right)-h{f}^{\prime \prime }+f{h}^{\prime \prime }-A_t\left(h^{\prime }+{\displaystyle \frac{1}{2}}\eta h^{\prime \prime }\right)\right]$$

(20)

The solutions of the equations above are assumed to be written as a power series in p, as given in Equation (21)

$$\begin{array}{c}\hfill f=f_0+p{f}_1+p^2{f}_2+p^3{f}_3+\dots \\ \hfill \theta ={\theta }_0+p{\theta }_1+p^2{\theta }_2+p^3{\theta }_3+\dots \\ \hfill \phi ={\phi }_0+p{\phi }_1+p^2{\phi }_2+p^3{\phi }_3+\dots \\ \hfill h=h_0+p{h}_1+p^2{h}_2+p^3{h}_3+\dots \end{array}$$

(21)

Substituting (21) and its respective derivatives into Equations (18)–(20), together with the boundary conditions in Equation (13) and arranging the equations according to the powers of the embedding parameter p

$$p^{\left(0\right)}:{\displaystyle \frac{1}{2}}A{\beta }_p\eta f_0^{\prime v}=0,{\displaystyle \frac{1+R_d}{P_r}}{\theta }_0^{\prime \prime }=0,{\displaystyle \frac{1}{Sc}}{\phi }_0^{\prime \prime }=0,\Lambda h_0^{\prime \prime \prime }=0$$

(22a)

$$f_0\left(ϵ\right)=-f_w,f_0^{\prime }\left(ϵ\right)=\lambda f_0^{\prime \prime }\left(ϵ\right),{\theta }_0^{\prime }\left(ϵ\right)=-Bi(1-{\theta }_0^{\prime }\left(ϵ\right)),N_b{\phi }_0^{\prime }\left(ϵ\right)+N_t{\theta }_0^{\prime }\left(ϵ\right)=0,h_0\left(ϵ\right)=h_0^{\prime \prime }\left(ϵ\right)=0$$

(22b)

$$f_0^{\prime }(\infty )\to 1,f_0^{\prime \prime }(\infty )\to 0,{\theta }_0(\infty )\to 0,{\phi }_0(\infty )\to 0,h_0^{\prime }(\infty )\to 1$$

(22c)

$$\begin{array}{c}\hfill p^{\left(1\right)}:\begin{array}{cc}& {\displaystyle \frac{1}{2}}A{\beta }_p\eta {\displaystyle \frac{d^4{f}_1}{d{\eta }^4}}+(2\eta y+1){\displaystyle \frac{d^3{f}_0}{d{\eta }^3}}+2y{\displaystyle \frac{d^2{f}_0}{d{\eta }^2}}+f_0{\displaystyle \frac{d^2{f}_0}{d{\eta }^2}}-{\left({\displaystyle \frac{d{f}_0}{d\eta }}\right)}^2-A\left({\displaystyle \frac{d{f}_0}{d\eta }}+{\displaystyle \frac{1}{2}}\eta {\displaystyle \frac{d^2{f}_0}{d{\eta }^2}}\right)\hfill \\ \hfill & +\mathrm{Ri}\left({\theta }_0-N_r{\phi }_0\right)+Q{Z}^2\left[{\left({\displaystyle \frac{d{h}_0}{d\eta }}\right)}^2-h_0{\displaystyle \frac{d^2{h}_0}{d{\eta }^2}}-1\right]+A{\beta }_p\left[(6\eta y+2){\displaystyle \frac{d^3{f}_0}{d{\eta }^3}}+4y{\displaystyle \frac{d^2{f}_0}{d{\eta }^2}}\right]\hfill \\ \hfill & +A\Omega -{\Omega }^2\hfill \end{array}=0\end{array}$$

(23a)

$$\begin{array}{cc}\hfill & f_1\left(ϵ\right)=0,f_1^{\prime }\left(ϵ\right)=\lambda f^{\prime \prime }\left(ϵ\right),f_1^{\prime }(\infty )\to 0,f_1^{\prime \prime }(\infty )\to 0\hfill \end{array}$$

(23b)

$$\begin{array}{c}\hfill \begin{array}{cc}& {\displaystyle \frac{1+R_d}{Pr}}{\displaystyle \frac{d^2{\theta }_1}{d{\eta }^2}}+{\displaystyle \frac{1+R_d}{Pr}}\left(2\eta \gamma {\displaystyle \frac{d^2{\theta }_0}{d{\eta }^2}}+2\gamma {\displaystyle \frac{d{\theta }_0}{d\eta }}\right)+E_c(2\eta \gamma +1){\left({\displaystyle \frac{d^2{f}_0}{d{\eta }^2}}\right)}^2+f_0{\displaystyle \frac{d{\theta }_0}{d\eta }}\hfill \\ \hfill & -n{\displaystyle \frac{d{f}_0}{d\eta }}{\theta }_0+G{\theta }_0+(2\eta \gamma +1)\left[N_b{\displaystyle \frac{d{\phi }_0}{d\eta }}{\displaystyle \frac{d{\theta }_0}{d\eta }}+N_t{\left({\displaystyle \frac{d{\theta }_0}{d\eta }}\right)}^2\right]+Q\Lambda Ec{\left({\displaystyle \frac{d^2{h}_0}{d{\eta }^2}}\right)}^2\hfill \end{array}=0\end{array}$$

(24a)

$$\begin{array}{cc}\hfill & {\theta }_1^{\prime }\left(ϵ\right)=Bi{\theta }_1\left(ϵ\right),{\theta }_1(\infty )=0\hfill \end{array}$$

(24b)

$$\begin{array}{c}\hfill \begin{array}{cc}& {\displaystyle \frac{1}{\mathrm{Sc}}}{\displaystyle \frac{d^2{\phi }_1}{d{\eta }^2}}+{\displaystyle \frac{1}{\mathrm{Sc}}}\left(2\eta \gamma {\displaystyle \frac{d^2{\phi }_0}{d{\eta }^2}}+2\gamma {\displaystyle \frac{d{\phi }_0}{d\eta }}\right)+{\displaystyle \frac{N_t}{\mathrm{Le}Pr{N}_b}}\left[(2\eta \gamma +1){\displaystyle \frac{d^2{\theta }_0}{d{\eta }^2}}+2\gamma {\displaystyle \frac{d{\theta }_0}{d\eta }}\right]\hfill \\ \hfill & -{\kappa }^2{\phi }_0+f_0{\displaystyle \frac{d{\phi }_0}{d\eta }}-n{\displaystyle \frac{d{f}_0}{d\eta }}{\phi }_0\hfill \end{array}=0\end{array}$$

(25a)

$$\begin{array}{cc}\hfill & N_b{\phi }_1^{\prime }\left(ϵ\right)+N_t{\theta }_1^{\prime }\left(ϵ\right)=0,{\phi }_1(\infty )\to 0\hfill \end{array}$$

(25b)

$$\begin{array}{cc}\hfill & \Lambda {\displaystyle \frac{d^3{h}_1}{d{\eta }^3}}+\Lambda \left(2\eta \gamma {\displaystyle \frac{d^3{h}_0}{d{\eta }^3}}+2\gamma {\displaystyle \frac{d^2{h}_0}{d{\eta }^2}}\right)-h_0{\displaystyle \frac{d^2{f}_0}{d{\eta }^2}}+f_0{\displaystyle \frac{d^2{h}_0}{d{\eta }^2}}-A_t\left({\displaystyle \frac{d{h}_0}{d\eta }}+{\displaystyle \frac{1}{2}}\eta {\displaystyle \frac{d^2{h}_0}{d{\eta }^2}}\right)=0\hfill \end{array}$$

(26a)

$$\begin{array}{cc}\hfill & h_1\left(ϵ\right)=h_1^{\prime \prime }\left(ϵ\right)=0,h_1^{\prime }(\infty )\to 0\hfill \end{array}$$

(26b)

$$\begin{array}{c}\hfill p^{\left(2\right)}:\begin{array}{cc}& {\displaystyle \frac{1}{2}}A{\beta }_p\eta {\displaystyle \frac{d^4{f}_2}{d{\eta }^4}}+(2\eta \gamma +1){\displaystyle \frac{d^3{f}_1}{d{\eta }^3}}+2\gamma {\displaystyle \frac{d^2{f}_1}{d{\eta }^2}}+f_0{\displaystyle \frac{d^2{f}_1}{d{\eta }^2}}+f_1{\displaystyle \frac{d^2{f}_0}{d{\eta }^2}}-2{\displaystyle \frac{d{f}_0}{d\eta }}{\displaystyle \frac{d{f}_1}{d\eta }}\hfill \\ \hfill & -A\left({\displaystyle \frac{d{f}_1}{d\eta }}+{\displaystyle \frac{1}{2}}\eta {\displaystyle \frac{d^2{f}_1}{d{\eta }^2}}\right)+R_i\left({\theta }_1-N_r{\phi }_1\right)+Q{Z}^2\left[2{\displaystyle \frac{d{h}_0}{d\eta }}{\displaystyle \frac{d{h}_1}{d\eta }}-h_0{\displaystyle \frac{d^2{h}_1}{d{\eta }^2}}-h_1{\displaystyle \frac{d^2{h}_0}{d{\eta }^2}}\right]\hfill \\ \hfill & +A{\beta }_p\left[{\eta }^2\gamma {\displaystyle \frac{d^4{f}_1}{d{\eta }^4}}+(6\eta \gamma +2){\displaystyle \frac{d^3{f}_1}{d{\eta }^3}}+4\gamma {\displaystyle \frac{d^2{f}_1}{d{\eta }^2}}\right]\hfill \end{array}=0\end{array}$$

(27a)

$$\begin{array}{cc}\hfill & f_2\left(ϵ\right)=0,f_2^{\prime }\left(ϵ\right)=\lambda f_2^{\prime \prime }\left(ϵ\right),f_2^{\prime }(\infty )\to 0,f_2^{\prime \prime }(\infty )\to 0\hfill \end{array}$$

(27b)

$$\begin{array}{c}\hfill \begin{array}{cc}& {\displaystyle \frac{1+R_d}{Pr}}{\displaystyle \frac{d^2{\theta }_2}{d{\eta }^2}}+{\displaystyle \frac{1+R_d}{Pr}}\left(2\eta \gamma {\displaystyle \frac{d^2{\theta }_1}{d{\eta }^2}}+2\gamma {\displaystyle \frac{d{\theta }_1}{d\eta }}\right)+2E_c(2\eta \gamma +1){\displaystyle \frac{d^2{f}_0}{d{\eta }^2}}{\displaystyle \frac{d^2{f}_1}{d{\eta }^2}}+f_0{\displaystyle \frac{d{\theta }_1}{d\eta }}+f_1{\displaystyle \frac{d{\theta }_0}{d\eta }}\hfill \\ \hfill & -n{\displaystyle \frac{d{f}_0}{d\eta }}{\theta }_1-n{\displaystyle \frac{d{f}_1}{d\eta }}{\theta }_0+G{\theta }_1+(2\eta \gamma +1)\left[N_b{\displaystyle \frac{d{\phi }_0}{d\eta }}{\displaystyle \frac{d{\theta }_1}{d\eta }}+N_b{\displaystyle \frac{d{\phi }_1}{d\eta }}{\displaystyle \frac{d{\theta }_0}{d\eta }}+2N_t{\displaystyle \frac{d{\theta }_0}{d\eta }}{\displaystyle \frac{d{\theta }_1}{d\eta }}\right]\hfill \\ \hfill & +2Q\Lambda E_c{\displaystyle \frac{d^2{h}_0}{d{\eta }^2}}{\displaystyle \frac{d^2{h}_1}{d{\eta }^2}}\hfill \end{array}=0\end{array}$$

(28a)

$$\begin{array}{cc}\hfill & {\theta }_2^{\prime }\left(ϵ\right)=Bi{\theta }_2\left(ϵ\right),{\theta }_2(\infty )=0\hfill \end{array}$$

(28b)

$$\begin{array}{c}\hfill \begin{array}{cc}& {\displaystyle \frac{1}{S_c}}{\displaystyle \frac{d^2{\phi }_2}{d{\eta }^2}}+{\displaystyle \frac{1}{S_c}}\left(2\eta \gamma {\displaystyle \frac{d^2{\phi }_1}{d{\eta }^2}}+2\gamma {\displaystyle \frac{d{\phi }_1}{d\eta }}\right)+{\displaystyle \frac{N_t}{L_{e}Pr{N}_b}}\left[(2\eta \gamma +1){\displaystyle \frac{d^2{\theta }_1}{d{\eta }^2}}+2\gamma {\displaystyle \frac{d{\theta }_1}{d\eta }}\right]-{\kappa }^2{\phi }_1\hfill \\ \hfill & +f_0{\displaystyle \frac{d{\phi }_1}{d\eta }}+f_1{\displaystyle \frac{d{\phi }_0}{d\eta }}-n{\displaystyle \frac{d{f}_0}{d\eta }}{\phi }_1-n{\displaystyle \frac{d{f}_1}{d\eta }}{\phi }_0\hfill \end{array}=0\end{array}$$

(29a)

$$\begin{array}{cc}\hfill & N_b{\phi }_2^{\prime }\left(ϵ\right)+N_t{\theta }_2^{\prime }\left(ϵ\right)=0,{\phi }_2(\infty )\to 0\hfill \end{array}$$

(29b)

$$\begin{array}{cc}& \Lambda {\displaystyle \frac{d^3{h}_2}{d{\eta }^3}}+\Lambda \left(2\eta \gamma {\displaystyle \frac{d^3{h}_1}{d{\eta }^3}}+2\gamma {\displaystyle \frac{d^2{h}_1}{d{\eta }^2}}\right)-h_0{\displaystyle \frac{d^2{f}_1}{d{\eta }^2}}-h_1{\displaystyle \frac{d^2{f}_0}{d{\eta }^2}}+f_0{\displaystyle \frac{d^2{h}_1}{d{\eta }^2}}+f_1{\displaystyle \frac{d^2{h}_0}{d{\eta }^2}}\hfill \end{array}$$

(30a)

$$\begin{array}{cc}\hfill & -A_t\left({\displaystyle \frac{d{h}_1}{d\eta }}+{\displaystyle \frac{1}{2}}\eta {\displaystyle \frac{d^2{h}_1}{d{\eta }^2}}\right)=0\hfill \end{array}$$

(30b)

$$\begin{array}{cc}\hfill & h_2\left(ϵ\right)=h_2^{\prime \prime }\left(ϵ\right)=0,h_2^{\prime }(\infty )\to 0\hfill \end{array}$$

(30c)

To obtain an approximate analytical solution of the transformed momentum balance, energy, nanoparticle volume-fraction, and the magnetic equations governing the MHD Rivlin–Ericksen fluid, physically consistent trial function are constructed for the dimensionless stream function $f\left(\eta \right)$, thermal function $\theta \left(\eta \right)$, nanoparticle function $\phi \left(\eta \right)$, and magnetic function $h\left(\eta \right)$. This approach has been widely adopted in boundary-layer studies involving non-Newtonian fluids and MHD effects, as it guarantees the correct asymptotic behaviour while retaining mathematical tractability.

In MHD boundary-layer flows, the velocity profile must satisfy the following fundamental requirements:

$\left(i\right).$ finite velocity at the wall due to suction/injection and slip effects; $(ii.)$ monotonic approach to the free-stream velocity as $\eta \to 0$, and $(iii.)$ exponential decay of shear stress far from the surface, reflecting boundary-layer confinement under magnetic damping.

To satisfy these conditions, the following composite exponential forms are assumed:

$$f_0\left(\eta \right)=a_1+b_1\eta +c_1{e}^{-k\eta }+d_1\eta e^{-k\eta }$$

(31)

$${\theta }_0\left(\eta \right)=a_2+b_2{e}^{-k\eta }$$

(32)

$${\phi }_0\left(\eta \right)=a_3+b_3{e}^{-k\eta }$$

(33)

$$h_0\left(\eta \right)=a_4+b_4\eta +c_4{e}^{-k\eta }$$

(34)

where $a_i,b_i,c_i,\mathrm{and}{d}_i$, $i=1,2,3,4$ are arbitrary constants to be determined and $k>0$ is a known decay parameter related to magnetic-field strength, fluid elasticity, or effective viscosity.

Solving Equation (31) subject to the boundary conditions in (22) to obtain

$$f_0\left(\eta \right)=-f_w-c_1+\eta +c_1{e}^{-k\eta }+d_1\eta e^{-k\eta },\mathbf{where}{d}_1={\displaystyle \frac{ck(\lambda k+1)-1}{1+2\lambda k}}$$

(35)

$${\theta }_0\left(\eta \right)={\displaystyle \frac{Bi}{k+Bi}}{e}^{-k\eta },$$

(36)

$${\phi }_0\left(\eta \right)=-{\displaystyle \frac{N_t}{N_b}}{\displaystyle \frac{Bi}{k+Bi}}{e}^{-k\eta },\mathbf{and}$$

(37)

$$h_0\left(\eta \right)=\eta$$

(38)

The constant $c_1$ remains free and can be established using the weighted residual method.

3.2. Weighted Residual (Galerkin) Method

The use of a Galerkin weighted residual technique to determine the remaining constant ensures that the governing momentum equation is satisfied in an averaged sense over the boundary-layer. Since the dominant near-wall correction is governed by the exponential decay, we select

$$W\left(\eta \right)=e^{k\eta }$$

(39)

The choice of an exponential weighting function:

• guarantees convergence of the integrals,
• emphasizes the near-wall region where MHD and viscoelastic effects are dominant,
• is mathematically consistent with the assumed solution structure.

The weighted residual condition is imposed as

$${\int }_0^{\infty }\mathcal{R}\left(\eta \right)W\left(\eta \right)d\eta =0.$$

(40)

For steady incompressible MHD Rivlin–Ericksen flow after HPM transformation, the momentum equation typically assumes the form

$${\displaystyle \frac{1}{2}}A{\beta }_p\eta f_0^{\prime v}\left(\eta \right)=0,$$

(41)

Let the residual of the momentum equation be defined as

$$\mathcal{R}\left(\eta \right)={\displaystyle \frac{1}{2}}A{\beta }_p\eta f_0^{\prime v}\left(\eta \right).$$

(42)

substituting Equations (39) and (42) into Equation (41)

$${\int }_0^{\infty }\left[{\displaystyle \frac{1}{2}}A{\beta }_p\eta f_0^{\prime v}\left(\eta \right)\right]e^{-k\eta }d\eta =0,$$

(43)

Solving Equation (43) by substituting the fourth derivatives of $f_0^{\prime v}$ for the free variable $c_1$ as

$$c_1={\displaystyle \frac{3}{2k+\lambda k^2}}$$

(44)

Thus, the initial approximate solution of the momentum equation is given as

$$f_0\left(\eta \right)=-f_w+\eta -{\displaystyle \frac{3}{2k+\lambda k^2}}\left(1-e^{-k\eta }\right)+{\displaystyle \frac{1}{2+\lambda k}}\eta e^{-k\eta }$$

(45)

substituting $c_1$ into (35) to get

$$f_0\left(\eta \right)=-f_w+\eta -{\displaystyle \frac{3}{2k+\lambda k^2}}\left(1-e^{-k\eta }\right)+{\displaystyle \frac{1}{2+\lambda k}}\eta e^{-k\eta },$$

(46)

$${\theta }_0\left(\eta \right)={\displaystyle \frac{Bi}{k+Bi}}{e}^{-k\eta },$$

(47)

$${\phi }_0\left(\eta \right)=-{\displaystyle \frac{N_t}{N_b}}{\displaystyle \frac{Bi}{k+Bi}}{e}^{-k\eta },\mathbf{and}$$

(48)

$$h_0\left(\eta \right)=\eta$$

(49)

Equation (46) facilitate the analytical computation of $f_i\left(\eta \right),{\theta }_i\left(\eta \right),{\phi }_i\left(\eta \right),\mathrm{and}{h}_i\left(\eta \right)$, $i=1,2,\dots$ together with their respective boundary conditions. However, the exact closed form for $f_i\left(\eta \right),{\theta }_i\left(\eta \right),{\phi }_i\left(\eta \right),\mathrm{and}{h}_i\left(\eta \right)$, $i=0,1,2,\dots$ are exceedingly lengthy due to

• numerous parametric values
• exponential integrals $E_i(-k\eta )$
• incomplete Gamma functions $\Gamma (\eta ,k\eta )$
• polynomial coefficients from multiple integrations

For practical applications, numerical values are substituted for all parameters, solve the boundary value problem in Maple 2016 using dsolve, and present results for $f_i\left(\eta \right),{\theta }_i\left(\eta \right),{\phi }_i\left(\eta \right),\mathrm{and}{h}_i\left(\eta \right)$, $i=1,2,\dots ,N$. Then the final results take the form

$$f\left(\eta \right)=f_0\left(\eta \right)+f_1\left(\eta \right)+f_2\left(\eta \right)+\dots$$

(50)

$$\theta \left(\eta \right)={\theta }_0\left(\eta \right)+{\theta }_1\left(\eta \right)+{\theta }_2\left(\eta \right)+\dots$$

(51)

$$\phi \left(\eta \right)={\phi }_0\left(\eta \right)+{\phi }_1\left(\eta \right)+{\phi }_2\left(\eta \right)+\dots$$

(52)

$$h\left(\eta \right)=h_0\left(\eta \right)+h_1\left(\eta \right)+h_2\left(\eta \right)+\dots$$

(53)

The adopted semi-analytical approach ensures exact satisfaction of all boundary conditions while capturing the essential nonlinear, magnetic, and viscoelastic effects inherent in Rivlin–Ericksen fluid flow.

4. Discussion of Results

This section examines the effects of key governing parameters on the flow, thermal, and mass-transport characteristics of a chemically reactive MHD Rivlin–Ericksen nanofluid over a cylindrical surface. The ensuing analysis from the research is substantiated by nine composite figures, each with subfigures and two tabulated datasets reporting the computed values of the skin-friction coefficient, Nusselt number, and Sherwood number. Throughout all graphical representations, directional arrows superimposed on the solution profiles denote the monotonic response of the respective field variables to successive increments in the governing parameter under consideration, thereby providing unambiguous visual indication of parametric sensitivity and directions.

4.1. Velocity Field and Momentum Transport

Figure 2 and Figure 3 illustrate the influence of the suction/injection parameter $f_w$, unsteadiness parameter A, viscoelastic parameter ${\beta }_p$, curvature parameter $\gamma$, Richardson number $Ri$, Chandrasekhar number Q, scaled axial position Z, and velocity ratio parameter $\Omega$ on the dimensionless velocity profile $f^{\prime }\left(\eta \right)$. Figure 2a demonstrates that increasing the suction parameter $(f_w>0)$ suppresses the velocity profile and reduces the momentum boundary-layer thickness. Suction extracts low-momentum fluid from the vicinity of the wall, thereby stabilizing the boundary layer and enhancing wall shear stress. Conversely, fluid injection $(f_w<0)$ supplies additional fluid into the boundary layer, resulting in thicker momentum layers and higher velocities away from the wall. The effect of the unsteadiness parameter A is displayed in Figure 2b. Larger values of A significantly increase the velocity throughout the boundary layer. Physically, increasing unsteadiness intensifies the temporal acceleration of the flow, promoting stronger momentum transport and delaying the decay of the velocity field toward the free-stream state. Figure 2c shows that the Rivlin–Ericksen viscoelastic parameter ${\beta }_p$ reduces the velocity profile. The viscoelastic stresses generated by fluid memory effects introduce an additional resistance to fluid motion, thereby suppressing momentum transport. The reduction is most pronounced near the cylinder surface where elastic stresses attain their highest magnitude. The influence of the curvature parameter $\gamma$ is illustrated in Figure 2d. Increasing $\gamma$ decreases the velocity profile and compresses the momentum boundary layer. Stronger curvature enhances geometric confinement of the flow around the cylinder, producing steeper velocity gradients and greater wall shear. Figure 3a depicts the effect of the Richardson number $Ri$. Increasing $Ri$ causes a slight reduction in velocity. Since $Ri$ measures the ratio of buoyancy forces to inertial forces, stronger buoyancy modifies the stagnation-point structure and increases flow resistance within the boundary layer. The Chandrasekhar number Q shown in Figure 3b significantly suppresses the velocity field. Larger values of Q correspond to stronger magnetic fields, which generate Lorentz forces opposing fluid motion. The resulting electromagnetic braking effect converts kinetic energy into internal energy and consequently reduces momentum transport. Figure 3c indicates that increasing the scaled axial position Z decreases the velocity profile. Since Z appears quadratically in the magnetic contribution of the momentum equation, larger axial positions strengthen the effective magnetic resistance and accelerate the decay of the velocity field. The influence of the velocity ratio parameter $\Omega$ is presented in Figure 3d. Increasing $\Omega$ enhances the velocity distribution throughout the boundary layer. Physically, a larger free-stream-to-wall velocity ratio increases the driving force of the external stagnation flow, resulting in stronger momentum transport and thicker velocity profiles.

Overall, suction, magnetic effects, viscoelasticity, curvature, Richardson number, and axial position suppress the velocity field, whereas larger unsteadiness and velocity ratio parameters promote fluid motion and enhance momentum transport.

4.2. Temperature Distribution and Heat Transfer

Figure 4 and Figure 5 display the effects of the unsteadiness parameter A, viscoelastic parameter ${\beta }_p$, Prandtl number $Pr$, Eckert number $Ec$, temperature exponent n, and curvature parameter $\gamma$ on the dimensionless temperature distribution $\theta$($\eta$).

Figure 4a reveals that increasing the unsteadiness parameter A enhances the temperature field and thickens the thermal boundary layer. Greater unsteadiness strengthens temporal energy accumulation within the fluid, allowing thermal energy to diffuse farther from the wall. Similarly, Figure 4b shows that larger values of the viscoelastic parameter ${\beta }_p$ increase the temperature profile. The suppression of fluid motion caused by viscoelastic stresses reduces convective cooling and promotes the retention of thermal energy within the boundary layer. The influence of the Prandtl number is presented in Figure 5a. Increasing $Pr$ markedly reduces the temperature distribution. Since $Pr$ is inversely proportional to thermal diffusivity, larger values weaken thermal diffusion and consequently produce thinner thermal boundary layers and higher heat-transfer rates. Figure 5b demonstrates that the Eckert number $Ec$ substantially elevates the temperature profile. The Eckert number measures the conversion of kinetic energy into thermal energy through viscous dissipation. As $Ec$ increases, additional heat is generated within the fluid, resulting in higher temperatures and thicker thermal layers. The effect of the temperature exponent n is shown in Figure 5c. Increasing n reduces the temperature distribution throughout the boundary layer. Larger values of n strengthen the temperature-dependent convective cooling mechanism represented in the energy equation, leading to lower thermal energy retention. Figure 5d indicates that increasing the curvature parameter $\gamma$ enhances the temperature field. Stronger curvature alters the thermal diffusion path around the cylindrical geometry and promotes thermal energy accumulation close to the wall, thereby thickening the thermal boundary layer.

From a heat-transfer perspective, parameters that increase temperature $(A,{\beta }_p,Ec,\mathrm{and}\gamma )$ generally reduce the wall temperature gradient and consequently diminish the Nusselt number, whereas increasing $Pr$ enhances heat transfer due to stronger thermal gradients at the wall.

4.3. Concentration Field and Mass Transport

Figure 6, Figure 7, Figure 8 and Figure 9 illustrate the influence of suction/injection $f_w$, unsteadiness A, viscoelastic parameter ${\beta }_p$, curvature parameter $\gamma$, Richardson number $Ri$, Chandrasekhar number Q, scaled axial position Z, velocity ratio $\Omega$, Prandtl number $Pr$, Eckert number $Ec$, temperature exponent n, Brownian motion parameter $N_b$, thermophoresis parameter $N_t$, Schmidt number $Sc$, and Lewis number $L_e$ on the nanoparticle concentration profile $\varphi \left(\eta \right)$.

Figure 6a shows that increasing suction decreases nanoparticle concentration and thins the concentration boundary layer. Suction removes nanoparticle-laden fluid from the vicinity of the wall, thereby enhancing mass-transfer rates. Injection produces the opposite effect. Figure 6b,c indicate that increasing either the viscoelastic parameter ${\beta }_p$ or the unsteadiness parameter A suppresses nanoparticle concentration. Both parameters strengthen the effective transport resistance and reduce species diffusion within the boundary layer. Figure 6d demonstrates that increasing the Richardson number $Ri$ decreases concentration. Stronger buoyancy effects modify the flow structure and enhance species transport away from the wall, resulting in lower concentration levels. The effects of the Chandrasekhar number Q, scaled axial position Z, curvature parameter $\gamma$, and velocity ratio $\Omega$ are shown in Figure 7. Increasing Q reduces concentration due to magnetic damping of fluid motion, which limits nanoparticle transport. Similarly, increasing Z lowers concentration because of enhanced magnetic influence along the axial direction. Increasing curvature $\gamma$ also suppresses concentration owing to geometric confinement effects. In contrast, larger values of $\Omega$ increase concentration by strengthening the external flow and promoting nanoparticle transport. Figure 8 illustrates the thermal parameters. Increasing the radiation parameter $R_d$ slightly elevates concentration because radiative heating enhances nanoparticle mobility. Increasing the Prandtl number $Pr$ reduces concentration due to weaker thermal diffusion and reduced nanoparticle transport. Conversely, higher Eckert numbers $Ec$ increase concentration through viscous-dissipation-induced heating. The temperature exponent n decreases concentration by strengthening thermal decay effects. The effects of Brownian motion $N_b$, thermophoresis $N_t$, Schmidt number $Sc$, and Lewis number $L_e$ are presented in Figure 9. Increasing $N_b$ substantially enhances concentration because random nanoparticle motion promotes species dispersion throughout the boundary layer. Similarly, larger thermophoretic parameters $N_t$ increase concentration by driving nanoparticles away from the heated wall into the fluid domain.

Figure 9c shows that increasing the Schmidt number $Sc$ decreases concentration. Since $Sc$ is inversely proportional to mass diffusivity, larger values reduce species diffusion and thin the concentration boundary layer. Likewise, Figure 9d indicates that increasing the Lewis number $L_e$ suppresses concentration because thermal diffusion becomes increasingly dominant over mass diffusion.

Consequently, suction, viscoelasticity, unsteadiness, magnetic effects, curvature, Richardson number, Schmidt number, and Lewis number reduce nanoparticle concentration, whereas thermophoresis, Brownian motion, velocity ratio, radiation, and viscous dissipation enhance species transport within the boundary layer.

4.4. Magnetic Profile

Figure 10 illustrates the effects of the local transient unsteadiness parameter $A_t$ and the inverse modified magnetic Prandtl number $\Lambda$ on the induced magnetic field profile $h^{\prime }\left(\eta \right)$. Figure 10a shows that increasing the local transient unsteadiness parameter $A_t$ suppresses the induced magnetic field. Larger values of $A_t$ strengthen temporal magnetic diffusion effects, thereby reducing magnetic induction generated by fluid motion. Consequently, the magnetic boundary layer becomes thinner and the induced field approaches its free-stream value more rapidly. Figure 10b demonstrates that increasing the inverse modified magnetic Prandtl number $\Lambda$ also reduces the magnetic profile. Physically, larger $\Lambda$ corresponds to stronger magnetic diffusivity relative to momentum diffusivity. Enhanced magnetic diffusion weakens the coupling between the fluid flow and the induced magnetic field, resulting in lower magnetic induction throughout the boundary layer. The results confirm the strong coupling between hydrodynamic and electromagnetic transport mechanisms. Parameters that enhance magnetic diffusion $(A_t\mathrm{and}\Lambda )$ weaken magnetic induction and reduce the thickness of the magnetic boundary layer, whereas lower values promote stronger magnetic field generation within the flow domain.

4.5. Skin Friction, Nusselt, and Sherwood Numbers

Table 1 and Table 2, together with Figure 9, summarize the effects of key parameters on the engineering quantities of interest. Skin friction increases with magnetic parameter, viscoelastic parameter, curvature, and suction due to enhanced wall shear, while slip and injection reduce it by weakening near-wall velocity gradients.

The Nusselt number decreases with increasing radiation and viscous dissipation as elevated temperatures reduce wall-temperature gradients. In contrast, a higher Prandtl number and suction enhance heat transfer by thinning the thermal boundary-layer.

The Sherwood number increases with Schmidt number and chemical reaction parameter owing to steeper concentration gradients at the wall. Brownian motion reduces mass-transfer, whereas thermophoresis enhances it, consistent with the observed concentration profiles.

4.6. Overall Physical Interpretation

Across all results, suction consistently enhances momentum, heat, and mass transfer by stabilizing the boundary-layer, whereas injection has a destabilizing effect. Magnetic and viscoelastic parameters exhibit similar suppressive trends on velocity but differ in their physical origins. The combined HPM–Galerkin approach yields physically consistent solutions that clearly capture these interactions while satisfying all boundary conditions.

5. Conclusions

In this study, a semi-analytical investigation of chemically reactive magnetohydrodynamic stagnation-point flow of a Rivlin–Ericksen nanofluid over a cylindrical surface has been carried out. The governing nonlinear boundary-layer equations were solved using a hybrid Homotopy Perturbation–Galerkin weighted residual approach, enabling accurate satisfaction of all boundary conditions and clear interpretation of parameter effects. Based on the presented results, the following conclusions are drawn:

1.

The applied magnetic field and Rivlin–Ericksen viscoelastic parameter both suppress the fluid velocity and enhance wall shear stress; however, the underlying mechanisms differ, with Lorentz forces dominating electromagnetic resistance and elastic stresses arising from fluid memory effects.

2.

Suction consistently stabilizes the boundary-layer by reducing its thickness, leading to increased skin friction, heat transfer, and mass-transfer rates, whereas injection produces the opposite effect by thickening the boundary-layers.

3.

Surface curvature intensifies near-wall velocity gradients and increases skin friction, while velocity slip weakens wall shear and reduces surface drag.

4.

Thermal radiation and viscous dissipation significantly elevate fluid temperature and thicken the thermal boundary layer, resulting in reduced Nusselt numbers due to diminished wall-temperature gradients.

5.

Brownian motion and thermophoresis both influence thermal and concentration fields, with thermophoresis exerting a stronger effect on boundary-layer thickness due to directed nanoparticle migration.

6.

Chemical reaction and Schmidt number strongly suppress nanoparticle concentration and enhance mass transfer rates, whereas Brownian motion reduces the Sherwood number by smoothing concentration gradients.

7.

The induced magnetic field is strongly coupled with the velocity field, increasing with magnetic Prandtl number and injection, and diminishing under suction.

Overall, the results demonstrate that the combined effects of viscoelasticity, nanoparticle transport, magnetic induction, and surface curvature play a crucial role in controlling momentum, heat, and mass transfer characteristics. The adopted hybrid semi-analytical methodology provides an efficient and reliable framework for analyzing complex MHD non-Newtonian nanofluid flows and can be readily extended to other geometries and transport mechanisms relevant to advanced thermal and industrial applications.