Parametric analysis of magnetic field-dependent viscosity and advection–diffusion between rotating discs

Abstract

The constitutive expressions of unsteady Newtonian fluid are employed in the mathematical formulation to model the flow between the circular space of porous and contracting discs. The flow behavior is investigated for magnetic field-dependent (MFD) viscosity and heat/mass transfers under the influence of a variable magnetic field. The equation for conservation of mass, modified Navier–Stokes, Maxwell, advection diffusion and transport equations are coupled as a system of ordinary differential equations. The expressions for torques and magnetohydrodynamic pressure gradient equation are derived. The MFD viscosity $\vartheta$ , magnetic Reynolds number ${\aleph }_{em}$ , squeezing Reynolds number ${\aleph }_b$ , rotational Reynolds number ${\aleph }_a$ , magnetic field components ${\aleph }_c$ , ${\aleph }_d$ , pressure $F_{\text{pres}}$ and the torques ${{\varrho }^{\prime }}_0$ , ${\varrho }_1$ which the fluid exerts on discs are discussed through numerical results and graphical aids. It is concluded that magnetic Reynolds number causes an increase in magnetic field distributions and decrease in tangential velocity of flow field, also the fluid temperature is decreasing with increase in magnetic Reynolds number. The azimuthal and axial components of magnetic field have opposite behavior with increase in MFD viscosity.

Introduction

Study of squeezing flow between parallel discs has attracted engineers and scientists due to its vast engineering applications such as compression and injection shaping, braking devices, pumping of heart, rotating machinery, bearing with liquid metal lubrication, hydraulic shock absorbers, modeling of metal and plastic sheets, thin fiber, paper sheets formations, and squeezed film in power transformations. According to Lenz’s law of motion of a conductor into a magnetic field, electric current is induced in the conductor and creates its own magnetic field. When currents are induced by motion of a conducting fluid through a magnetic field, a Lorentz force acts on the fluid and modifies its motion. In magnetohydrodynamics (MHD), the motion modifies the field and vice versa. This makes the theory highly nonlinear. ^1,2 Hughes and Elco ³ studied the squeezed flow of an electrically conducting fluid between two discs in the presence of a magnetic field. It has been shown that the load capacity of the normal force which the fluid exerts on the upper disc is dependent on the MHD interactions in the fluid. Verma ⁴ also studied the squeeze film lubrication of a magnetic fluid between two approaching surfaces in the presence of an externally applied magnetic field. Krieger et al. ⁵ studied the MHD lubrication flow between parallel stationary discs in an axial magnetic field. They also obtained excellent agreement between theory and experimental results until the transition to turbulent flow occurred. Rahimi-Gorji et al. ⁶ investigated the unsteady squeezing nanofluid flow and heat transfer using Galerkin method for the solution of governing nonlinear differential equations in the presence of variable magnetic field. Acharya et al. ⁷ studied the squeezing flow of two types of nanofluids such as Cu-water and Cu-kerosene between parallel plates in the presence of variable magnetic field. The nonlinear differential equations are solved numerically by RK4 method with shooting technique and analytically using differential transformation method. Saidi and Tamim ⁸ investigated the unsteady three-dimensional nanofluid flow, heat and mass transfer in a rotating system in the presence of an externally applied magnetic field. Khan et al. ⁹ studied the MHD flow of a viscous incompressible fluid between parallel discs. The governing system of equations is solved using homotopy analysis method (HAM). Khan and Shah ¹⁰ investigated the magnetic squeeze film flow between rotating discs with induced magnetic field effects taking the radial, azimuthal, and axial components of magnetic field as

$B_r=\frac{\delta r}{2D(1-\delta t)}m\prime (\eta ),B_{\theta }=\frac{r{N}_o}{\sqrt{1-\delta t}}n(\eta ),B_z=\frac{-\delta M_o}{\sqrt{1-\delta t}}m(\eta )$

with velocity components

$u=\frac{\delta r}{1-\delta t}f\prime (\eta ),v=\frac{\Omega r}{1-\delta t}g,w=\frac{\delta r}{1-\delta t}f(\eta )$

computing the response of radial and azimuthal magnetic fields to squeezing rates and relative disc rotation. The flow of a fluid film squeezed between two rotating parallel plane surfaces in the presence of a magnetic field applied perpendicular to the surface was studied by Hamza, ¹¹ and it showed that the magnetic forces and the centrifugal inertial forces have opposite effects. Shah et al. ¹² studied the effects of magnetic field on lubrication flow between two discs and showed that the magnetic field increases the fluid load on upper disc. Kuzma et al. ¹³ studied the MHD squeeze film by taking the fluid-inertial effect and buoyant forces. Their results showed an excellent agreement between theory and experiment. Hayat et al. ¹⁴ studied the squeezing flow of nanofluid between two parallel plates in the presence of magnetic field. They neglected the induced magnetic field for small magnetic Reynolds number and concluded that temperature of fluid is decreasing for large values of squeezing parameter. Nabhani and El Khlifi ¹⁵ numerically investigated the squeezing film between two porous circular discs in the presence of an externally applied magnetic field. An implicit finite difference scheme is used to discretize the governing unsteady nonlinear equation and solve the new system of equations by Gauss–Seidel method. The radiation effect on squeezing flow between parallel discs is studied by Mohyud-Din and Khan. ¹⁶ HAM is employed to obtain the expression for velocity and temperature profiles. They observed an opposite behavior of the velocity profile for suction and injection of fluid for all involved parameters. Hatami et al. ¹⁷ investigated the heat transfer of nanofluid flow between parallel plates in the presence of variable magnetic field. Using homotopy perturbation method, they concluded that Nusselt number has direct relationship with Brownian motion parameter. Heat transfer analysis of squeezing flow between parallel discs is studied by Azimi and Riazi ¹⁸ for GO-water nanofluid flowthe squeeze film lubrication of a magnetic fluid between. Verma ⁴ studied the squeeze film lubrication of a magnetic fluid between two approaching surfaces in the presence of an externally applied magnetic field. He assumed that the applied magnetic field M has components of the form

$M_x=M(x)\text{cos}\theta , M_y=M(x)\text{sin}\theta , M_z=0 \hspace{0.17em}\text{where \hspace{0.05em}} \hspace{0.17em}\theta =\theta (x,y)$

The squeezing flow with rotating discs under the influence of variable magnetic field is also studied in detail by Ibrahim ¹⁹ and Kumari et al. ²⁰ Sunil et al. ²¹ and Domairry and Aziz ²² conducted a theoretical investigation of the combined effect of magnetic field-dependent (MFD) viscosity and rotation on ferroconvection in the presence of dust particles subjected to a transverse uniform magnetic field and provided the approximate analytic solution for the squeezing flow of viscous fluid between parallel discs with suction or blowing. There is great interest in the rotation of fluids in a rotating magnetic field with viscosity, where large rotational velocities can be achieved. This phenomenon can be used for the construction of various kinds of pumps, automotive magnetorheological shock absorbers, novel aircraft landing gear systems, biological prosthetics, and also for nonmechanical mixing of fluids, for example, spacecraft fuels. Because of its great scientific and applied interest, it is attempted to discuss the influence of rotation and how MFD viscosity affects the magnetization in fluid.

Nomenclature

p	Pressure (N m−2)	C	Dimensional concentration
$r,\theta ,z$	Cylindrical polar coordinates	$\overrightarrow{r}$	Radius vector of the disc
ur	Radial velocity (m s−1)	$\overrightarrow{V}$	Velocity vector
$u_{\theta }$	Azimuthal velocity (m s−1)	Greek symbols
uz	Axial velocity (m s−1)	$\omega$	Rotation vector
t	Time (s)	${\Omega }_l$	Lower disc angular velocity
Tu	Temperature at upper disc (K)	$\kappa$	Thermal conductivity (W m−1 K−1)
Tl	Temperature at lower disc (K)	$\mu$	Dynamic viscosity $(\text{Pa}\hspace{0.17em}\text{s})$
Cu	Concentration at upper disc	$\nu$	Kinematic viscosity (kg m−1 s−1)
Cl	Concentration at lower disc	$\sigma$	Relative angular velocity
$Pr$	Prandtl number (v/k)	$\rho$	Fluid density (kg m−3)
${\aleph }_c$	Strength of magnetic field	α	Positive constant
${\aleph }_b$	Squeezing parameter	$\eta$	Similarity variable
${\aleph }_d$	Strength of magnetic field	$\Gamma$	Transformed fluid temperature
$D(t)$	Distance between two discs (m)	${\sigma }^s$	Stefan–Boltzmann constant
B	Induced magnetic field	No	Constant number
${\delta }_o$	Soret number	$\Psi$	Transformed fluid concentration
$\lambda$	Dufour number	${\kappa }^a$	Mean absorption coefficient
cp	Specific heat of fluid (J kg−1 K−1)	Subscript
D	Molecular diffusion coefficient	u	Fluid condition on upper disc
kT	Thermal diffusion ratio	l	Fluid condition on lower disc
Tm	Mean fluid temperature (K)	Superscript
qr	Radiative heat flux (W m−2)	$\ast$	Dimensionless variable
${\Re }_d$	Radiative parameter	$\text{'}$	Derivative with respect to $\eta$
${\delta }_c$	Schmidt number

Existing information on the topic witnessed that the MFD viscosity of a squeezing flow between porous discs for viscous fluid under the influence of heat/mass transfers and externally applied variable magnetic field in polar coordinates has never been reported, and this is the very first study in the literature. In the following sections, the problem is formulated, analyzed, and discussed through graphs and tables.

Formulation of the problem

Consider an unsteady, incompressible, laminar, and electrically conducting viscous fluid flow caused by the rotation of two squeezing discs, parallel and smooth as shown in Figure 1. The discs are rotating with different angular velocities proportional to $\frac{{\Omega }_l}{1-\alpha t}$ and $\frac{{\Omega }_u}{1-\alpha t}$ , where ${\Omega }_l$ , ${\Omega }_u$ are the representative angular velocities of lower and upper discs, respectively, with dimension $t^{-1}$ . ²³ The distance between discs at time t is $D(t)=l{(1-\alpha t)}^{0.5}$ , $\alpha t\le 1$ apart, where l is representative length equivalent to the discs separation at $t=0$ . ¹⁹ The lower porous disc is prohibited from moving in the axial direction, where upper disc is rotating and moving toward or away from constrained lower disc. The fluid is assumed to be incompressible having a variable viscosity ²⁴ given by $\mu ={\mu }_o(1+{\delta }_v\mathrm{.}B)$ , where ${\mu }_o$ is the viscosity of the fluid when there is no magnetic field. The variation coefficient of viscosity ${\delta }_v$ has been taken to be isotropic, that is, ${\delta }_{vr}={\delta }_{v\theta }={\delta }_{vz}$ . Hence, $\mu$ in the component form can be written as ${\mu }_r={\mu }_o(1+{\delta }_v{B}_r)$ , ${\mu }_{\theta }={\mu }_o(1+{\delta }_v{B}_{\theta })$ , and ${\mu }_z={\mu }_o(1+{\delta }_v{B}_z)$ . Further, B, M, and H are related by $B={\mu }_p(M+H)$ where M is the magnetization when the magnetic field is H and ${\mu }_p$ is the magnetic permeability of vacuum. The effect of shear dependence on viscosity is not considered since it has negligible effect for a mono-dispersive system of large rotation and high field. As a first approximation for small field variation, linear variation of magneto-viscosity has been used, hence $B={\mu }_p(M_o+H_o)$ where H_o is uniform magnetic field and M_o is the magnetization when the magnetic field is H_o . The movable disc is affected by the magnetic field H defined as

$H_r=\frac{\alpha r{M}_o}{{\mu }_2(1-\alpha t)}, H_{\theta }=\frac{r{N}_o}{{\mu }_2(1-\alpha t)}, H_z=\frac{-r{M}_o}{{\mu }_1{(1-\alpha t)}^{0.5}}$

where M_o and N_o are used to dimensionless H_r , $H_{\theta }$ , H_z and ${\mu }_1$ , ${\mu }_2$ are the magnetic permeability of outside and inside media between the two discs, respectively. Following experimental study, the above parameters of magnetic field are zero on the lower disc. ²⁴ The induced magnetic field $(B_r,B_{\theta },B_z)$ in the fluid is generated by the magnetic field H. Under these assumptions, the governing equations in vector form are ^25,26 :

Equation of conservation of mass is

$\nabla \cdot U=0$

1

• Equation of conservation of momentum is

$\rho \left[\frac{\partial U}{\partial t}+(U\cdot \nabla )U\right]=-\nabla p+\mu {\nabla }^{2}U+\frac{1}{{\mu }_2}\left[(\nabla \times B)\times B\right]$

2

• Maxwell’s equations are

$\nabla \cdot B=0,\frac{\partial B}{\partial t}=\nabla \times (V\times B)+\frac{1}{\varrho {\mu }_2}{\nabla }^{2}B$

3

• Equations of heat and mass transfer are

Figure 1.

Geometry of the problem.

$\rho (V\cdot \nabla )T=\frac{1}{C_p}\nabla \cdot (k\nabla )T+\frac{1}{C_p}\nabla q_r+\frac{D{k}_T}{C_s{c}_p}\nabla \cdot (\nabla C)$

4

$\rho (V\cdot \nabla )C=\nabla \cdot (D\nabla )C+\frac{D{k}_T}{C_s{T}_m}\nabla \cdot (\nabla T)$

5

Boundary conditions

The boundary conditions are chosen as

$\begin{array}{l}{u}_r=0,u_{\theta }=\frac{{\Omega }_{l}r}{1-\alpha t},u_z=-\frac{w_o}{{(1-\alpha t)}^{0.5}},B_r=0,B_{\theta }=0,B_z=0,T=T_l,C=C_l\hspace{1em}\text{at}\hspace{1em}z=0\hfill \\ u_r=0,\hspace{1em}{u}_{\theta }=\frac{{\Omega }_{u}r}{1-\alpha t},\hspace{1em}{u}_z=\frac{-0.5\alpha l}{{(1-\alpha t)}^{0.5}},\hspace{1em}{B}_r=0,\hspace{1em}{B}_{\theta }=\frac{r{N}_o}{1-\alpha t},\hspace{1em}{B}_z=\frac{-r{M}_o}{{(1-\alpha t)}^{0.5}}\hfill \\ T=T_u,\hspace{1em}C=C_u\hspace{1em}\hspace{1em}\text{at}\hspace{1em}z=D(t)\hfill \end{array}$

6

where u_r , $u_{\theta }$ , and u_z are respectively the radial, azimuthal, and axial components of velocity; p is the pressure; T is temperature; C is concentration; ${\alpha }_t$ is thermal diffusivity; $\rho$ is fluid density; D is diffusion coefficient; T_m is mean fluid temperature; T_l and C_l denote the temperature and concentration at the lower disc while T_u and C_u are the temperature and concentration at the upper disc, respectively; $\nu$ is kinematic viscosity; $\varrho$ is electrical conductivity; c_p is specific heat at constant pressure; $\kappa$ is conductivity; ${\omega }_0$ is velocity of upper disc; ${\Omega }_u,{\Omega }_l$ are the angular velocities of upper and lower discs; ${\kappa }_T$ is thermal diffusion ratio; ${\sigma }^s$ is Stefan–Boltzmann constant; ${\kappa }^a$ is the mean absorption coefficient; and q_r is the radiative heat flux such that $\nabla q_{{}_r}=\frac{16{\sigma }^s{T}_u^3}{3{\kappa }^a}\frac{\partial T^2}{\partial z^2}$ (Khan and Shah ¹⁰ ).

The following similarity transformations ²⁷ are chosen for reducing the partial differential equations (1) to (6) from cylindrical coordinates to a system of ordinary differential equations

$\begin{array}{l}{u}_r=r\frac{\partial }{\partial z}F(z,t)=0.5\frac{\alpha r}{1-\alpha t}f\prime (\eta ),\hspace{1em}{u}_{\theta }=rG(z,t)=\frac{{\Omega }_{l}r}{1-\alpha t}g(\eta )\hspace{1em}\hfill \\ u_z=-2F(z,t)=-\frac{\alpha l}{{(1-\alpha t)}^{0.5}}f(\eta ),B_r=r\frac{\partial }{\partial z}M(z,t)=0.5\frac{\alpha r{M}_o}{l^2(1-\alpha t)}m\prime (\eta )\hspace{1em}\hfill \\ B_{\theta }=rN(z,t)=\frac{2\nu r{N}_o}{\alpha l(1-\alpha t)}n(\eta ),\hspace{1em}{B}_z=-2M(z,t)=-\frac{\alpha M_o}{l{(1-\alpha t)}^{0.5}}m(\eta )\hspace{1em}\hfill \\ \Gamma =\frac{T-T_u}{T_l-T_u},\hspace{1em}\Psi =\frac{C-C_u}{C_l-C_u}\hspace{1em}\text{where}\hspace{1em}\eta =\frac{z}{l\sqrt{1-\alpha t}}\hfill \end{array}$

7

Equation of continuity is identically satisfied and the momentum, Maxwell’s, heat and mass transfer equations take the following form

$\vartheta {\aleph }_b^{2}f″″-{\aleph }_b^3\left[3f″+\eta f″\prime -2ff″\prime +2{\aleph }_c^2{\aleph }_{em}\left(\eta mm″+mm\prime +2m^{2}f″-2fmm″\right)\right]+2{\aleph }_a^2(gg\prime -{\aleph }_d^{2}nn\prime )=0$

8

$\vartheta g″-{\aleph }_b(2g+\eta g\prime +2gf\prime -2fg\prime )+2{\aleph }_c{\aleph }_d(mn\prime -nm\prime )=0$

9

$m″-{\aleph }_{em}(m+\eta m\prime -2fm\prime +2mf\prime )=0$

10

${\aleph }_{d}n″-{\aleph }_d{\aleph }_{em}(2n+\eta n\prime -2fn\prime )-2{\aleph }_{c}mg\prime =0$

11

$(3{\Re }_d+4)\Gamma ″+3{\Re }_d\lambda \Psi ″-6{\Re }_d{P}_r{\aleph }_{b}f\Gamma \prime =0$

12

$\Psi ″+{\delta }_c{\delta }_o\Gamma ″+2{\delta }_c{\aleph }_{b}f\Psi \prime =0$

13

and the boundary conditions are reduced to

$\begin{array}{l}f(0)=\delta ,\hspace{1em}f\prime (0)=0,\hspace{1em}g(0)=1,\hspace{1em}\Gamma (0)=1,\hspace{1em}\Psi (0)=1,\hspace{1em}m(0)=0,\hspace{1em}n(0)=0\hfill \\ f(1)=0.5,\hspace{1em}f\prime (1)=0,\hspace{1em}g(1)=\frac{{\Omega }_u}{{\Omega }_l}=\sigma ,\hspace{1em}\Gamma (1)=0,\hspace{1em}\Psi (1)=0,\hspace{1em}m(1)=1,\hspace{1em}n(1)=1\hfill \end{array}$

14

where $\sigma =\frac{{\Omega }_u}{{\Omega }_l}$ is the relative angular velocity of discs, ${\aleph }_b=\frac{\alpha l^2}{2\nu }$ is squeeze Reynolds number, ${\aleph }_a=\frac{{\Omega }_l{l}^2}{\nu }$ is rotational Reynolds number, $B_t=\varrho {\mu }_2\nu$ is Bachelor number, ${\aleph }_{em}={\aleph }_b{B}_t$ is magnetic Reynolds number, $\vartheta =1+{\delta }_v{\mu }_p(M_o+H_o)$ is MFD viscosity parameter, ${\aleph }_c=\frac{M_o}{l\sqrt{{\mu }_2\rho }}$ is strength of magnetic field in z direction, ${\aleph }_d=\frac{N_o}{\Omega \sqrt{{\mu }_2\rho }}$ is strength of magnetic field in $\theta$ direction, $\delta =\frac{w_0}{\alpha l}$ is suction/injection parameter, $P_r=\frac{\nu }{{\kappa }_T}$ is Prandtl number, ${\delta }_c=\frac{\nu }{D}$ is Schmidt number, ${\Re }_d=\frac{{\kappa }^a\kappa }{4{\sigma }^s{T}_u^3}$ is radiation parameter, ${\delta }_o=\frac{D(T_l-T_u)}{\nu T_m(C_l-C_u)}$ is Soret number, and $\lambda =\frac{D(C_l-C_u){\kappa }_T}{C_s{C}_p\nu (T_l-T_u)}$ is the Dufour number.

Approximate analytical solution

The analytic method HAM ^24,28 is used to solve system of equations (8) to (13). Due to HAM, the functions $f(\eta )$ , $g(\eta )$ , $m(\eta )$ , $n(\eta )$ , $\Gamma (\eta )$ , and $\Psi (\eta )$ can be expressed by a set of base functions ${\eta }^c, c\ge 0$ as

$f_m(\eta )={\displaystyle \sum _{\zeta =0}^{\infty }}{a}_{\zeta }{\eta }^{\zeta }$

15

$g_m(\eta )={\displaystyle \sum _{\zeta =0}^{\infty }}{b}_{\zeta }{\eta }^{\zeta }$

16

$m_m(\eta )={\displaystyle \sum _{\zeta =0}^{\infty }}{c}_{\zeta }{\eta }^{\zeta }$

17

$n_m(\eta )={\displaystyle \sum _{\zeta =0}^{\infty }}{d}_{\zeta }{\eta }^{\zeta }$

18

${\Gamma }_m(\eta )={\displaystyle \sum _{\zeta =0}^{\infty }}{e}_{\zeta }{\eta }^{\zeta }$

19

${\Psi }_m(\eta )={\displaystyle \sum _{\zeta =0}^{\infty }}{f}_{\zeta }{\eta }^{\zeta }$

20

where $a_{\zeta }$ , $b_{\zeta }$ , $c_{\zeta }$ , $d_{\zeta }$ , $e_{\zeta }$ , and $f_{\zeta }$ are the constant coefficients to be determined. Initial approximations are chosen as follows

$f_0(\eta )=(2A-1){\eta }^3-\frac{3}{2}(2A-1){\eta }^2+A$

21

$g_0(\eta )=(\sigma -1)\eta +1$

22

$m_0(\eta )=\eta$

23

$n_0(\eta )=\eta$

24

${\Gamma }_0(\eta )=1-\eta$

25

${\Psi }_0(\eta )=1-\eta$

26

The auxiliary operators are chosen as

$$\begin{gathered} {\ell }_f=\frac{{\partial }^4}{\partial {\eta }^4},\hspace{1em}{\ell }_g=\frac{{\partial }^2}{\partial {\eta }^2},\hspace{1em}{\ell }_m=\frac{{\partial }^2}{\partial {\eta }^2},\hspace{1em}{\ell }_n=\frac{{\partial }^2}{\partial {\eta }^2},\hspace{1em} \\ {\ell }_{\Gamma }=\frac{{\partial }^2}{\partial {\eta }^2},\hspace{1em}{\ell }_{\Psi }=\frac{{\partial }^2}{\partial {\eta }^2} \end{gathered}$$

27

with the following properties

${\ell }_f({\zeta }_1{\eta }^3+{\zeta }_2{\eta }^2+{\zeta }_3\eta +{\zeta }_4)=0$

28

${\ell }_g({\zeta }_5\eta +{\zeta }_6)=0$

29

${\ell }_m({\zeta }_7\eta +{\zeta }_8)=0$

30

${\ell }_n({\zeta }_9\eta +{\zeta }_{10})=0$

31

${\ell }_{\Gamma }({\zeta }_{11}\eta +{\zeta }_{12})=0$

32

${\ell }_{\Psi }({\zeta }_{13}\eta +{\zeta }_{14})=0$

33

where ${\zeta }_1$ , ${\zeta }_2$ , ${\zeta }_3$ , ${\zeta }_4$ , ${\zeta }_5$ , ${\zeta }_6$ , ${\zeta }_7$ , ${\zeta }_8$ , ${\zeta }_9$ , ${\zeta }_{10}$ , ${\zeta }_{11}$ , ${\zeta }_{12}$ , ${\zeta }_{13}$ , and ${\zeta }_{14}$ are arbitrary constants.

The zeroth order deformation problems can be obtained as

$(1;\beta ){\ell }_f\left[\overline{f}(\eta ;\beta )-f_0(\eta )\right]=q{\hslash }_f{N}_f\left[\overline{f}(\eta ;\beta ),\overline{g}(\eta ;\beta ),\overline{m}(\eta ;\beta ),\overline{n}(\eta ;\beta )\right]$

34

$(1;\beta ){\ell }_g\left[\overline{g}(\eta ;\beta )-g_0(\eta )\right]=q{\hslash }_g{N}_g\left[\overline{f}(\eta ;\beta ),\overline{g}(\eta ;\beta ),\overline{m}(\eta ;\beta ),\overline{n}(\eta ;\beta )\right]$

35

$(1;\beta ){\ell }_m\left[\overline{m}(\eta ;\beta )-m_0(\eta )\right]=q{\hslash }_m{N}_m\left[\overline{f}(\eta ;\beta ),\overline{m}(\eta ;\beta ),\overline{n}(\eta ;\beta )\right]$

36

$(1;\beta ){\ell }_n\left[\overline{n}(\eta ;\beta )-n_0(\eta )\right]=q{\hslash }_n{N}_n\left[\overline{f}(\eta ;\beta ),\overline{m}(\eta ;\beta ),\overline{n}(\eta ;\beta )\right]$

37

$(1;\beta ){\ell }_{\Gamma }\left[\overline{\Gamma }(\eta ;\beta )-{\Gamma }_0(\eta )\right]=q{\hslash }_{\Gamma }{N}_{\Gamma }\left[\overline{f}(\eta ;\beta ),\overline{\Gamma }(\eta ;\beta ),\overline{\Psi }(\eta ;\beta )\right]$

38

$(1;\beta ){\ell }_{\Psi }\left[\overline{\Psi }(\eta ;\beta )-{\Psi }_0(\eta )\right]=q{\hslash }_{\Psi }{N}_{\Psi }\left[\overline{f}(\eta ;\beta ),\overline{\Gamma }(\eta ;q),\overline{\Psi }(\eta ;\beta )\right]$

39

The nonlinear operators of equations (9) to (13) are defined as

$$\begin{gathered} N_f\left[\overline{f}(\eta ;\beta ),\overline{g}(\eta ;\beta ),\overline{m}(\eta ;\beta ),\overline{n}(\eta ;\beta )\right] \\ =\vartheta {\aleph }_b^2\frac{{\partial }^4\overline{f}(\eta ;\beta )}{\partial {\eta }^4}-{\aleph }_b^3[\eta \frac{{\partial }^3\overline{f}(\eta ;\beta )}{\partial {\eta }^3}+3\frac{{\partial }^2\overline{f}(\eta ;\beta )}{\partial {\eta }^2}-2f\frac{{\partial }^3\overline{f}(\eta ;\beta )}{\partial {\eta }^3}+2{\aleph }_c^2{\aleph }_{em}(\eta \overline{m}(\eta ;\beta )\frac{{\partial }^2\overline{m}(\eta ;\beta )}{\partial {\eta }^2}+\overline{m}(\eta ;\beta )\frac{\partial \overline{m}(\eta ;\beta )}{\partial \eta } \\ +2{\overline{m}}^2(\eta ;\beta )\frac{{\partial }^2\overline{f}(\eta ;\beta )}{\partial {\eta }^2}-2\overline{m}(\eta ;\beta )\overline{f}(\eta ;\beta )\frac{{\partial }^2\overline{m}(\eta ;\beta )}{\partial {\eta }^2})]+2{\aleph }_d^2\left(\overline{g}(\eta ;\beta )\frac{\partial \overline{g}(\eta ;\beta )}{\partial \eta }-{\aleph }_d^2\overline{n}(\eta ;\beta )\frac{\partial \overline{n}(\eta ;\beta )}{\partial \eta }\right) \end{gathered}$$

40

$$\begin{gathered} N_g\left[\overline{f}(\eta ;\beta ),\overline{g}(\eta ;\beta ),\overline{m}(\eta ;\beta ),\overline{n}(\eta ;\beta )\right]=\vartheta \frac{{\partial }^2\overline{g}(\eta ;\beta )}{\partial {\eta }^2}-{\aleph }_b^2[2\overline{g}(\eta ;\beta )+\eta \frac{\partial \overline{g}(\eta ;\beta )}{\partial \eta }+2\overline{g}(\eta ;\beta )\frac{\partial \overline{f}(\eta ;\beta )}{\partial \eta } \\ -2\overline{f}(\eta ;\beta )\frac{\partial \overline{g}(\eta ;\beta )}{\partial \eta }+2{\aleph }_c{\aleph }_d\left(m\frac{\partial \overline{n}(\eta ;\beta )}{\partial \eta }-\overline{n}(\eta ;\beta )\frac{\partial \overline{m}(\eta ;\beta )}{\partial \eta }\right)] \end{gathered}$$

41

$N_m\left[\overline{f}(\eta ;\beta ),\overline{m}(\eta ;\beta ),\overline{n}(\eta ;\beta )\right]=\frac{{\partial }^2\overline{m}(\eta ;\beta )}{\partial {\eta }^2}-{\aleph }_{em}(\overline{m}(\eta ;\beta )+\eta \frac{\partial \overline{m}(\eta ;\beta )}{\partial \eta }-2\overline{f}(\eta ;\beta )\frac{\partial \overline{m}(\eta ;\beta )}{\partial \eta }+2\overline{m}(\eta ;\beta )\frac{\partial \overline{f}(\eta ;\beta )}{\partial \eta })$

42

$N_n\left[\overline{f}(\eta ;\beta ),\overline{m}(\eta ;\beta ),\overline{n}(\eta ;\beta )\right]={\aleph }_c\frac{{\partial }^2\overline{n}(\eta ;\beta )}{\partial {\eta }^2}-{\aleph }_d{\aleph }_{em}(2\overline{n}(\eta ;\beta )+\eta \frac{\partial \overline{n}(\eta ;\beta )}{\partial \eta }-2f\frac{\partial \overline{n}(\eta ;\beta )}{\partial \eta })-2{\aleph }_c\overline{m}(\eta ;\beta )\frac{\partial \overline{g}(\eta ;\beta )}{\partial \eta }$

43

$N_{\Gamma }\left[\overline{f}(\eta ;\beta ),\overline{\Gamma }(\eta ;\beta ),\overline{\Psi }(\eta ;\beta )\right]=(3{\Re }_d+4)\frac{{\partial }^2\overline{\Gamma }(\eta ;\beta )}{\partial {\eta }^2}+3{\Re }_d\lambda \frac{{\partial }^2\overline{\Psi }(\eta ;\beta )}{\partial {\eta }^2}-6{\Re }_d{P}_r{\aleph }_b\overline{f}(\eta ;\beta )\frac{\partial \overline{\Gamma }(\eta ;\beta )}{\partial \eta }$

44

$N_{\Psi }\left[\overline{f}(\eta ;\beta ),\overline{\Gamma }(\eta ;\beta ),\overline{\Psi }(\eta ;\beta )\right]=\frac{{\partial }^2\overline{\Psi }(\eta ;\beta )}{\partial {\eta }^2}+{\delta }_c{\delta }_o\frac{{\partial }^2\overline{\Gamma }(\eta ;\beta )}{\partial {\eta }^2}+2{\delta }_c{\aleph }_b\overline{f}(\eta ;\beta )\frac{\partial \overline{\Psi }(\eta ;\beta )}{\partial \eta }$

45

where $\beta$ is an embedding parameter; ${\hslash }_f$ , ${\hslash }_g$ , ${\hslash }_m$ , ${\hslash }_n$ , ${\hslash }_{\Gamma }$ , and ${\hslash }_{\Psi }$ are the nonzero auxiliary parameters; and N_f , N_g , N_m , N_n , $N_{\Gamma }$ , and $N_{\Psi }$ are the nonlinear parameters.

For $\beta =0 \text{and} 1$ , we have

$\begin{array}{l}\overline{f}(\eta \mathrm{,0})=f_o(\eta ),\hspace{1em}\hspace{1em}\overline{f}(\eta \mathrm{,1})=f(\eta )\hfill \\ \overline{g}(\eta \mathrm{,0})=g_o(\eta ),\hspace{1em}\hspace{1em}\overline{g}(\eta \mathrm{,1})=g(\eta )\hfill \\ \overline{m}(\eta \mathrm{,0})=m_o(\eta ),\hspace{1em}\hspace{1em}\overline{m}(\eta \mathrm{,1})=m(\eta )\hfill \\ \overline{n}(\eta \mathrm{,0})=n_o(\eta ),\hspace{1em}\hspace{1em}\overline{n}(\eta \mathrm{,1})=n(\eta )\hfill \\ \overline{\Gamma }(\eta \mathrm{,0})={\Gamma }_o(\eta ),\hspace{1em}\hspace{1em}\overline{\Gamma }(\eta \mathrm{,1})=\Gamma (\eta )\hfill \\ \overline{\Psi }(\eta \mathrm{,0})={\Psi }_o(\eta ),\hspace{1em}\overline{\Psi }(\eta \mathrm{,1})=\Psi (\eta )\hfill \end{array}$

46

so we can say that as $\beta$ varies from 0 to 1, $\overline{f}(\eta \mathrm{,0})$ , $\overline{g}(\eta \mathrm{,0})$ , $\overline{m}(\eta \mathrm{,0})$ , $\overline{n}(\eta \mathrm{,0})$ , $\overline{\Gamma }(\eta \mathrm{,0})$ , and $\overline{\Psi }(\eta \mathrm{,0})$ vary from initial guesses $f_0(\eta )$ , $g_0(\eta )$ , $m_0(\eta )$ , $n_0(\eta )$ , ${\Gamma }_0(\eta )$ , and ${\Psi }_0(\eta )$ to exact solution $f(\eta )$ , $g(\eta )$ , $m(\eta )$ , $n(\eta )$ , $\Gamma (\eta )$ , and $\Psi (\eta )$ , respectively.

Taylor’s series expansion of these functions yields

$f(\eta ;\beta )=f_0(\eta )+{\displaystyle \sum _{m=1}^{\infty }}{\beta }^m{f}_m(\eta )$

47

$g(\eta ;\beta )=g_0(\eta )+{\displaystyle \sum _{m=1}^{\infty }}{\beta }^m{g}_m(\eta )$

48

$m(\eta ;\beta )=m_0(\eta )+{\displaystyle \sum _{m=1}^{\infty }}{\beta }^m{m}_m(\eta )$

49

$n(\eta ;\beta )=n_0(\eta )+{\displaystyle \sum _{m=1}^{\infty }}{\beta }^m{n}_m(\eta )$

50

$\Gamma (\eta ;\beta )={\Gamma }_0(\eta )+{\displaystyle \sum _{m=1}^{\infty }}{\beta }^m{\Gamma }_m(\eta )$

51

$\Psi (\eta ;\beta )={\Psi }_0(\eta )+{\displaystyle \sum _{m=1}^{\infty }}{\beta }^m{\Psi }_m(\eta )$

52

$\begin{array}{l}{f}_m(\eta )={\frac{1}{m!}\frac{{\partial }^{m}f(\eta ;\beta )}{\partial {\eta }^m}|}_{\beta =0},g_m(\eta )={\frac{1}{m!}\frac{{\partial }^{m}g(\eta ;\beta )}{\partial {\eta }^m}|}_{\beta =0},m_m(\eta )={\frac{1}{m!}\frac{{\partial }^{m}m(\eta ;\beta )}{\partial {\eta }^m}|}_{\beta =0}\hfill \\ n_m(\eta )={\frac{1}{m!}\frac{{\partial }^{m}n(\eta ;\beta )}{\partial {\eta }^m}|}_{\beta =0},{\Gamma }_m(\eta )={\frac{1}{m!}\frac{{\partial }^m\Gamma (\eta ;\beta )}{\partial {\eta }^m}|}_{\beta =0},{\Psi }_m(\eta )={\frac{1}{m!}\frac{{\partial }^m\Psi (\eta ;\beta )}{\partial {\eta }^m}|}_{\beta =0}\hfill \end{array}$

53

It should be noted that the convergence of above series strongly depends upon ${\hslash }_f$ , ${\hslash }_g$ , ${\hslash }_m$ , ${\hslash }_n$ , ${\hslash }_{\Gamma }$ , and ${\hslash }_{\Psi }$ .

Assuming that these nonzero auxiliary parameters are chosen so that equations (34) to (39) converge at $\beta =1$ . Therefore one can obtain

$f(\eta )=f_0(\eta )+{\displaystyle \sum _{m=1}^{\infty }}{f}_m(\eta )$

54

$g(\eta )=g_0(\eta )+{\displaystyle \sum _{m=1}^{\infty }}{g}_m(\eta )$

55

$\Gamma (\eta )={\Gamma }_0(\eta )+{\displaystyle \sum _{m=1}^{\infty }}{\Gamma }_m(\eta )$

56

$\Psi (\eta )={\Psi }_0(\eta )+{\displaystyle \sum _{m=1}^{\infty }}{\Psi }_m(\eta )$

57

Differentiating the deformation equations (34) to (39) $m \text{times}$ with respect to $\beta$ and putting $\beta =0$ , we have

${\ell }_f\left[f_m(\eta )-{\chi }_m{f}_{m-1}(\eta )\right]={\hslash }_f{R}_{f,m}(\eta )$

58

${\ell }_g\left[g_m(\eta )-{\chi }_m{g}_{m-1}(\eta )\right]={\hslash }_g{R}_{g,m}(\eta )$

59

${\ell }_m\left[m_m(\eta )-{\chi }_m{m}_{m-1}(\eta )\right]={\hslash }_m{R}_{m,m}(\eta )$

60

${\ell }_n\left[n_m(\eta )-{\chi }_m{n}_{m-1}(\eta )\right]={\hslash }_n{R}_{n,m}(\eta )$

61

${\ell }_{\Gamma }\left[{\Gamma }_m(\eta )-{\chi }_m{\Gamma }_{m-1}(\eta )\right]={\hslash }_{\Gamma }{R}_{\Gamma ,m}(\eta )$

62

${\ell }_{\Psi }\left[{\Psi }_m(\eta )-{\chi }_m{\Psi }_{m-1}(\eta )\right]={\hslash }_{\Psi }{R}_{\Psi ,m}(\eta )$

63

subject to the boundary conditions

$\begin{array}{l}{f}_m(0)=A,\hspace{1em}{f\prime }_m(0)=0,\hspace{1em}{g}_m(0)=0,\hspace{1em}{m}_m(0)=0,\hspace{1em}{n}_m(0)=0,\hspace{1em}{\Gamma }_m(0)=0,\hspace{1em}{\Psi }_m(0)=1\hfill \\ f_m(1)=0.5,\hspace{1em}{f\prime }_m(1)=0,\hspace{1em}{g}_m(1)=\sigma ,\hspace{1em}{m}_m(1)=1,\hspace{1em}{n}_m(1)=1,\hspace{1em}{\Gamma }_m(1)=0,\hspace{1em}{\Psi }_m(1)=0\hfill \end{array}$

64

where

$\begin{array}{c}{R}_{f,m}(\eta )=\vartheta {\aleph }_b^2{{f″}^{″}}_{m-1}(\eta )-{\aleph }_b^3[3{f″}_{m-1}(\eta )+\eta {{f″}^{\prime }}_{m-1}(\eta )-2{\displaystyle \sum _{j=0}^{m-1}}{f}_j(\eta ){{f″}^{\prime }}_{m-j-1}\\ +2{\aleph }_c^2{\aleph }_{em}{\displaystyle \sum _{j=0}^{m-1}}{m}_j(\eta )\left(\eta {m″}_{m-j-1}+{m\prime }_{m-j-1}+2m_{m-j-1}{f″}_{m-j-1}-2f_{m-j-1}{m″}_{m-j-1}\right)]\\ +2{\aleph }_d^2\left({\displaystyle \sum _{j=0}^{m-1}}{g}_j(\eta ){g\prime }_{m-j-1}(\eta )-{\aleph }_d^2{\displaystyle \sum _{j=0}^{m-1}}{n}_j(\eta ){n\prime }_{m-j-1}(\eta )\right)-8\Omega {g\prime }_{m-1}(\eta )-M^2{f″}_{m-1}(\eta )\end{array}$

65

$\begin{array}{c}{R}_{g,m}(\eta )=\vartheta {g″}_{m-1}(\eta )-{\aleph }_b[\eta {g\prime }_{m-1}(\eta )+2g_{m-1}+2{\displaystyle \sum _{j=0}^{m-1}}{g}_j(\eta ){f\prime }_{m-j-1}(\eta )-2{\displaystyle \sum _{j=0}^{m-1}}{f}_j(\eta ){g\prime }_{m-j-1}(\eta )\\ +2{\aleph }_c{\aleph }_d\left({\displaystyle \sum _{j=0}^{m-1}}{m}_j(\eta ){n\prime }_{m-j-1}(\eta )-{\displaystyle \sum _{j=0}^{m-1}}{n}_j(\eta ){m\prime }_{m-j-1}(\eta )\right)]\end{array}$

66

$\begin{array}{c}{R}_{m,m}(\eta )={m″}_{m-1}(\eta )-{\aleph }_{em}[\eta {m\prime }_{m-1}(\eta )+m_{m-1}-2{\displaystyle \sum _{j=0}^{m-1}}{f}_j(\eta ){m\prime }_{m-j-1}(\eta )+2{\displaystyle \sum _{j=0}^{m-1}}{m}_j(\eta ){f\prime }_{m-j-1}(\eta )\\ -{\displaystyle \sum _{j=0}^{m-1}}{n}_j(\eta ){m\prime }_{m-j-1}(\eta ))]\end{array}$

67

$\begin{array}{c}{R}_{n,m}(\eta )={\aleph }_d{n″}_{m-1}(\eta )\\ -{\aleph }_{em}{\aleph }_d\left[\eta {n\prime }_{m-1}(\eta )+2n_{m-1}-2{\displaystyle \sum _{j=0}^{m-1}}{f}_j(\eta ){n\prime }_{m-j-1}(\eta )\right]\\ -2{\aleph }_c{\displaystyle \sum _{j=0}^{m-1}}{m}_j(\eta ){g\prime }_{m-j-1}(\eta )\end{array}$

68

$$\begin{gathered} R_{\Gamma ,m}(\eta )=(3{\Re }_d+4){\Gamma ″}_{m-1}(\eta )+3{\Re }_d\lambda {\Psi ″}_{m-1}(\eta ) \\ -6{\Re }_d{P}_r{\aleph }_b{\displaystyle \sum _{j=0}^{m-1}}{f}_j(\eta ){\Gamma \prime }_{m-j-1}(\eta ) \end{gathered}$$

69

$$\begin{gathered} R_{\Psi ,m}(\eta )={\Psi ″}_{m-1}(\eta )+{\delta }_c{\delta }_o{\Gamma ″}_{m-1}(\eta ) \\ +2{\aleph }_b{\delta }_c{\displaystyle \sum _{j=0}^{m-1}}{f}_j(\eta ){\Psi \prime }_{m-j-1}(\eta ) \end{gathered}$$

70

and ${\chi }_m=\{1, \text{if} m>1, \text{and} 0, m=1$ .

Finally, the general solution of equations (58) to (63) can be written as

$$\begin{gathered} f_m(\eta )={\displaystyle {\int }_0^{\eta }}{\displaystyle {\int }_0^{\eta }}{\displaystyle {\int }_0^{\eta }}{\displaystyle {\int }_0^{\eta }}{\hslash }_f{R}_{f,m}(z)dzdzdzdzdz+{\chi }_m{f}_{m-1}+{\zeta }_1{\eta }^3 \\ +{\zeta }_2{\eta }^2+{\zeta }_3\eta +{\zeta }_4 \end{gathered}$$

71

$g_m(\eta )={\displaystyle {\int }_0^{\eta }}{\displaystyle {\int }_0^{\eta }}{\hslash }_g{R}_{g,m}(z)dzdz+{\chi }_m{g}_{m-1}+{\zeta }_5\eta +{\zeta }_6$

72

$m_m(\eta )={\displaystyle {\int }_0^{\eta }}{\displaystyle {\int }_0^{\eta }}{\hslash }_m{R}_{m,m}(z)dzdz+{\chi }_m{m}_{m-1}+{\zeta }_7\eta +{\zeta }_8$

73

$n_m(\eta )={\displaystyle {\int }_0^{\eta }}{\displaystyle {\int }_0^{\eta }}{\hslash }_n{R}_{n,m}(z)dzdz+{\chi }_m{n}_{m-1}+{\zeta }_9\eta +{\zeta }_{10}$

74

${\Gamma }_m(\eta )={\displaystyle {\int }_0^{\eta }}{\displaystyle {\int }_0^{\eta }}{\hslash }_{\Gamma }{R}_{\Gamma ,m}(z)dzdz+{\chi }_m{\Gamma }_{m-1}+{\zeta }_{11}\eta +{\zeta }_{12}$

75

${\Psi }_m(\eta )={\displaystyle {\int }_0^{\eta }}{\displaystyle {\int }_0^{\eta }}{\hslash }_{\Psi }{R}_{\Psi ,m}(z)dzdz+{\chi }_m{\Psi }_{m-1}+{\zeta }_{13}\eta +{\zeta }_{14}$

76

and so the exact solution $f(\eta )$ , $g(\eta )$ , $m(\eta )$ , $n(\eta )$ , $\Gamma (\eta )$ , and $\Psi (\eta )$ becomes

$\begin{array}{l}f(\eta )\approx {\displaystyle \sum _{n=0}^m}{f}_n(\eta ),\hspace{1em}g(\eta )\approx {\displaystyle \sum _{n=0}^m}{g}_n(\eta ),\hspace{1em}m(\eta )\approx {\displaystyle \sum _{n=0}^m}{m}_n(\eta )\hfill \\ n(\eta )\approx {\displaystyle \sum _{n=0}^m}{n}_n(\eta ),\hspace{1em}\Gamma (\eta )\approx {\displaystyle \sum _{n=0}^m}{\Gamma }_n(\eta ),\hspace{1em}\Psi (\eta )\approx {\displaystyle \sum _{n=0}^m}{\Psi }_n(\eta )\hfill \end{array}$

77

Optimal convergence control parameters

It must be remarked that the series solutions (58) to (63) contain the nonzero auxiliary parameters ${\hslash }_f$ , ${\hslash }_g$ , ${\hslash }_m$ , ${\hslash }_n$ , ${\hslash }_{\Gamma }$ , and ${\hslash }_{\Psi }$ which determine the convergence region and also rate of the homotopy series solutions. To obtain the optimal values of ${\hslash }_f$ , ${\hslash }_g$ , ${\hslash }_m$ , ${\hslash }_n$ , ${\hslash }_{\Gamma }$ , and ${\hslash }_{\Psi }$ , here the so-called average residual error defined by Liao ²⁸ was used as

${\epsilon }_m^f=\frac{1}{\zeta +1}{\displaystyle \sum _{j=0}^{\zeta }}{\left[N_f{\left({\displaystyle \sum _{i=0}^m}\overline{f}(\eta ),{\displaystyle \sum _{i=0}^m}\overline{g}(\eta ),{\displaystyle \sum _{i=0}^m}\overline{m}(\eta ),{\displaystyle \sum _{i=0}^m}\overline{n}(\eta )\right)}_{n=j\delta n}\right]}^{2}d\eta$

78

${\epsilon }_m^g=\frac{1}{\zeta +1}{\displaystyle \sum _{j=0}^{\zeta }}{\left[N_g{\left({\displaystyle \sum _{i=0}^m}\overline{f}(\eta ),{\displaystyle \sum _{i=0}^m}\overline{g}(\eta ),{\displaystyle \sum _{i=0}^m}\overline{m}(\eta ),{\displaystyle \sum _{i=0}^m}\overline{n}(\eta )\right)}_{n=j\delta n}\right]}^{2}d\eta$

79

${\epsilon }_m^m=\frac{1}{\zeta +1}{\displaystyle \sum _{j=0}^{\zeta }}{\left[N_m{\left({\displaystyle \sum _{i=0}^m}\overline{f}(\eta ),{\displaystyle \sum _{i=0}^m}\overline{m}(\eta ),{\displaystyle \sum _{i=0}^m}\overline{n}(\eta )\right)}_{n=j\delta n}\right]}^{2}d\eta$

80

${\epsilon }_m^n=\frac{1}{\zeta +1}{\displaystyle \sum _{j=0}^{\zeta }}{\left[N_n{\left({\displaystyle \sum _{i=0}^m}\overline{f}(\eta ),{\displaystyle \sum _{i=0}^m}\overline{m}(\eta ),{\displaystyle \sum _{i=0}^m}\overline{n}(\eta )\right)}_{n=j\delta n}\right]}^{2}d\eta$

81

${\epsilon }_m^{\Gamma }=\frac{1}{\zeta +1}{\displaystyle \sum _{j=0}^{\zeta }}{\left[N_{\Gamma }{\left({\displaystyle \sum _{i=0}^m}\overline{f}(\eta ),{\displaystyle \sum _{i=0}^m}\overline{\Gamma }(\eta ),{\displaystyle \sum _{i=0}^m}\overline{\Psi }(\eta )\right)}_{n=j\delta n}\right]}^{2}d\eta$

82

${\epsilon }_m^{\Psi }=\frac{1}{\zeta +1}{\displaystyle \sum _{j=0}^{\zeta }}{\left[N_{\Psi }{\left({\displaystyle \sum _{i=0}^m}\overline{f}(\eta ),{\displaystyle \sum _{i=0}^m}\overline{\Gamma }(\eta ),{\displaystyle \sum _{i=0}^m}\overline{\Psi }(\eta )\right)}_{n=j\delta n}\right]}^{2}d\eta$

83

Due to ²⁸

${\epsilon }_m^t={\epsilon }_m^f+{\epsilon }_m^g+{\epsilon }_m^m+{\epsilon }_m^n+{\epsilon }_m^{\Gamma }+{\epsilon }_m^{\Psi }$

84

where ${\epsilon }_m^t$ is the total squared residual error. The total average squared residual error is minimized by employing Mathematica package BVPh $2.0$ . ²⁸

The torques exerted on the discs

The frictional moment or torque which the fluid exerts on the upper disc is given ²⁴

${\tau }_u=2\pi \mu {\displaystyle {\int }_0^1}{r}^2{\left(\frac{\partial u_{\theta }}{\partial z}\right)}_{z=D}dr$

but $u_{\theta }={\Omega }_{u}r{(1-\alpha t)}^{-1}g(\eta )$ so the above equation becomes

$\begin{array}{l}{\tau }_u=\frac{\pi \mu {\Omega }_u{c}^4}{2l{(1-\alpha t)}^{1.5}}g\prime (1)\\ \text{or \hspace{0.05em}} {\tau }_u^{*}=g\prime (1) \text{where} {\tau }_u^{*}=\frac{2l{(1-\alpha t)}^{1.5}}{\pi \mu {\Omega }_u{c}^4}\end{array}$

85

For the lower disc, the corresponding result is

${\tau }_l^{*}=g\prime (0)$

86

which are the dimensionless exerted torques of the fluid on upper and lower discs.

The normal force or pressure of fluid on discs

According to Shah et al., ²⁴ the pressure or the normal force which the fluid exerts on the upper disc is given as

$F=2\pi \left[{\displaystyle {\int }_0^1}rP(r\mathrm{,1},t)dr-{\displaystyle {\int }_0^a}r{P}^{+}(r\mathrm{,1},t)dr\right]$

87

where $P^{+}(r\mathrm{,1},t)$ in for conditions on the side of the disc and $p(r\mathrm{,1},t)$ denotes the pressure at the edge of the disc at time t. Let us assume that $\frac{\partial P^{+}(r\mathrm{,1},t)}{\partial r}=0$ and using equation (8), we get

$$\begin{gathered} \frac{\partial p}{\partial r}=\frac{\rho r}{{(1-\alpha t)}^2}[\frac{\alpha \nu }{2l^2}\vartheta f″\prime +{\Omega }_u^2{g}^2-\frac{{\alpha }^2}{4}(2f\prime +\eta f″+f″-2ff″)-\frac{2}{\rho {\mu }_2}\left(\frac{{\alpha }^2{M}_o^2}{4l^2}mm″+N_o^2{n}^2\right)] \\ \Upsilon (t,\eta )=\frac{1}{{(1-\alpha t)}^2}[\frac{\alpha \nu }{2l^2}\vartheta f″\prime +{\Omega }_u^2{g}^2-\frac{{\alpha }^2}{4}(2f\prime +\eta f″+f″-2ff″)-\frac{2}{\rho {\mu }_2}\left(\frac{{\alpha }^2{M}_o^2}{4l^2}mm″+N_o^2{n}^2\right)] \end{gathered}$$

88

where $\Upsilon (t,\eta )=\frac{1}{\rho r}\frac{\partial p}{\partial r}$ . Using equations (87) and (88), we have

$\begin{array}{l}F=\frac{\pi \rho {\alpha }^2{a}^4}{16(1-\alpha t{)}^2}\left[2{\aleph }_c^{2}m″(1)-{\aleph }_b^{-1}\vartheta f″\prime (1)-{\left(\frac{{\aleph }_a}{{\aleph }_b}\right)}^2({\sigma }^2-2{\aleph }_d^2)\right]\\ F_{\text{pres}}=\left[2{\aleph }_c^{2}m″(1)-{\aleph }_b^{-1}\vartheta f″\prime (1)-{\left(\frac{{\aleph }_a}{{\aleph }_b}\right)}^2({\sigma }^2-2{\aleph }_d^2)\right]\end{array}$

89

where $F_{\text{pres}}=\frac{16(1-\alpha t{)}^2}{\pi \rho {\alpha }^2{a}^4}$ , which is the dimensionless pressure on the upper disc. The positive or negative numerical values of F _pres will be according to the force which the fluid acting on the upper disc is in the positive or negative direction of the z-axis, respectively.

Error analysis

The problem under consideration is solved by $\text{Mathematica}$ package BVPh 2.0 for a maximum residual error ${10}^{-40}$ using a computing machine with installed memory of $4.00$ GB and processor of Intel(R) core(TM) $\text{i}5-2450\text{M}$ CPU $2.50$ GHz. Analyses are carried out using $40\text{th}$ order approximations. Error analysis performed in Figure 2 and tabulated results given in Tables 1 to 9 support the authentification of results for different involved physical parameters.

Figure 2.

Total residual error profile of (a) $f(\eta )$ , (b) $g(\eta )$ , (c) $m(\eta )$ , (d) $n(\eta )$ , (e) $\Gamma (\eta )$ , and (f) $\Psi (\eta )$ for fixed values of $\lambda =2$ , $\delta =1$ , $P_r=0.1$ , ${\delta }_c=1$ , ${\delta }_o=0.2$ , ${\Re }_d=1$ , ${\aleph }_{em}=0.2$ , ${\aleph }_b=0.1$ , ${\aleph }_c=0.1$ , ${\aleph }_d=0.5$ , ${\aleph }_a=0.05$ , $\vartheta =2$ , and $\sigma =1$ .

Table 1.

Total residual error for different order of approximations taking fixed values of $\lambda =2$ , $\delta =1$ , $P_r=0.1$ , ${\delta }_c=1$ , ${\delta }_o=0.2$ , ${\Re }_d=1$ , ${\aleph }_{em}=0.2$ , ${\aleph }_b=0.1$ , ${\aleph }_c=0.1$ , ${\aleph }_d=0.5$ , ${\aleph }_a=0.05$ , $\vartheta =2$ , and $\sigma =1$ .

m	${\epsilon }^f{}_m$	${\epsilon }^g{}_m$	${\epsilon }^m{}_m$	${\epsilon }^n{}_m$	${\epsilon }^{\Gamma }{}_m$	${\epsilon }^{\Psi }{}_m$	CPU time (s)
2	$6.7\times {10}^{-7}$	$5.0\times {10}^{-7}$	$7.7\times {10}^{-7}$	$6.1\times {10}^{-7}$	$6.8\times {10}^{-3}$	$9.6\times {10}^{-3}$	0.69
6	$2.3\times {10}^{-11}$	$7.2\times {10}^{-14}$	$1.2\times {10}^{-13}$	$1.7\times {10}^{-14}$	$8.8\times {10}^{-6}$	$1.5\times {10}^{-6}$	5.2
11	$6.4\times {10}^{-17}$	$9.7\times {10}^{-20}$	$2.0\times {10}^{-19}$	$5.3\times {10}^{-21}$	$6.5\times {10}^{-8}$	$1.7\times {10}^{-10}$	11.7
16	$1.7\times {10}^{-22}$	$1.6\times {10}^{-25}$	$3.3\times {10}^{-25}$	$9.2\times {10}^{-27}$	$5.8\times {10}^{-12}$	$1.7\times {10}^{-13}$	24.3
21	$4.4\times {10}^{-28}$	$2.9\times {10}^{-31}$	$5.0\times {10}^{-31}$	$1.7\times {10}^{-32}$	$7.1\times {10}^{-15}$	$3.3\times {10}^{-17}$	44.1
26	$1.1\times {10}^{-33}$	$4.8\times {10}^{-34}$	$3.4\times {10}^{-33}$	$1.3\times {10}^{-34}$	$1.4\times {10}^{-18}$	$1.8\times {10}^{-20}$	74.3
31	$7.8\times {10}^{-36}$	$4.6\times {10}^{-34}$	$3.5\times {10}^{-33}$	$1.2\times {10}^{-34}$	$7.6\times {10}^{-22}$	$5.3\times {10}^{-24}$	119.1
36	$8.1\times {10}^{-36}$	$4.6\times {10}^{-34}$	$3.5\times {10}^{-33}$	$1.2\times {10}^{-34}$	$2.6\times {10}^{-25}$	$1.8\times {10}^{-27}$	184.5
40	$8.1\times {10}^{-36}$	$4.6\times {10}^{-34}$	$3.5\times {10}^{-33}$	$1.2\times {10}^{-34}$	$5.2\times {10}^{-28}$	$2.8\times {10}^{-30}$	256.5

Table 2.

Computations for $f(\eta )$ , $g(\eta )$ , $m(\eta )$ , $n(\eta )$ , $\theta (\eta )$ , and $\Psi (\eta )$ with $\lambda =2$ , $\delta =1$ , $P_r=0.1$ , ${\delta }_c=1$ , ${\delta }_o=0.2$ , ${\Re }_d=1$ , ${\aleph }_{em}=0.2$ , ${\aleph }_b=0.1$ , ${\aleph }_c=0.1$ , ${\aleph }_d=0.5$ , ${\aleph }_a=0.05$ , $\vartheta =2$ and $\sigma =1$ and various values of $\eta$.

$\eta$	HAM result					Numerical result
	$f(\eta )$	$g(\eta )$	$m(\eta )$	$\Gamma (\eta )$	$\Psi (\eta )$	$f(\eta )$	$g(\eta )$	$m(\eta )$	$\Gamma (\eta )$	$\Psi (\eta )$
0.1001	0.0140	0.0993	0.0929	0.9015	0.8981	0.0140	0.0993	0.0929	0.9015	0.8981
0.2002	0.0521	0.9876	0.1862	0.8030	0.7962	0.0522	0.9876	0.1862	0.8030	0.7962
0.3003	0.1082	0.9834	0.2805	0.7044	0.6944	0.1083	0.9834	0.2805	0.7044	0.6944
0.4004	0.1763	0.9808	0.3716	0.6055	0.5929	0.1763	0.9808	0.3716	0.6055	0.5930
0.5005	0.2504	0.9800	0.4737	0.5063	0.4919	0.2504	0.9800	0.4737	0.5063	0.4919
0.6006	0.3244	0.9809	0.5737	0.4065	0.3915	0.3244	0.9809	0.5737	0.4065	0.3915
0.7007	0.3924	0.9835	0.6763	0.3059	0.2918	0.3924	0.9835	0.6763	0.3060	0.2918
0.8008	0.4483	0.9876	0.7818	0.2046	0.1931	0.4483	0.9876	0.7818	0.2046	0.1931
0.9009	0.4862	0.9932	0.8902	0.1023	0.0955	0.4862	0.9932	0.8902	0.1023	0.0955

HAM: homotopy analysis method.

Table 3.

Convergence of HAM solution for different orders of approximation for $f″(0)$ , $-g^{\prime }(0)$ , $-m^{\prime }(0)$ , $-n^{\prime }(0)$ , $-\Gamma \prime (0)$ , and $-\Psi \prime (0)$ when $\lambda =2$ , $\delta =1$ , $P_r=0.1$ , ${\delta }_c=1$ , ${\delta }_o=0.2$ , ${\Re }_d=1$ , ${\aleph }_{em}=0.2$ , ${\aleph }_b=0.1$ , ${\aleph }_c=0.1$ , ${\aleph }_d=0.5$ , ${\aleph }_a=0.05$ , $\vartheta =2$ , and $\sigma =1$ .

m	$f″(0)$	$-g^{\prime }(0)$	$-m^{\prime }(0)$	$-n^{\prime }(0)$	$-\Gamma \prime (0)$	$-\Psi \prime (0)$
1	3.00595917	0.07676626	−0.92589583	−0.93460354	0.99932281	1.01406377
5	3.00701050	0.07387897	−0.92738629	−0.93331303	0.98418369	1.01812071
10	3.00701099	0.07387896	−0.92738629	−0.93331304	0.98371264	1.01825226
15	3.00701099	0.07387896	−0.92738629	−0.93331304	0.98370961	1.01825520
20	3.00701099	0.07387896	−0.92738629	−0.93331304	0.98370953	1.01825522
25	3.00701099	0.07387896	−0.92738629	−0.93331304	0.98370953	1.01825522
30	3.00701099	0.07387896	−0.92738629	−0.93331304	0.98370953	1.01825522
35	3.00701099	0.07387896	−0.92738629	−0.93331304	0.98370953	1.01825522
40	3.00701099	0.07387896	−0.92738629	−0.93331304	0.98370953	1.01825522

HAM: homotopy analysis method.

Table 4.

Optimal values of convergence control parameters versus different orders of approximation with fixed values of $\lambda =2$ , $\delta =1$ , $P_r=0.1$ , ${\delta }_c=1$ , ${\delta }_o=0.2$ , ${\Re }_d=1$ , ${\aleph }_{em}=0.2$ , ${\aleph }_b=0.1$ , ${\aleph }_c=0.1$ , ${\aleph }_d=0.5$ , ${\aleph }_a=0.5$ , $\vartheta =2$ , and $\sigma =1$ .

Order	hf	hg	hm	hn	$h_{\Gamma }$	$h_{\Psi }$	${\epsilon }_m^t$	CPU time (min)
2	−102.065	−1.0029	−1.0396	−2.0345	−3.1942	−1.0741	$5.0684\times {10}^{-2}$	0.1
3	−302.151	−1.6707	−1.3943	−2.4366	−1.7222	−0.5789	$4.7062\times {10}^{-2}$	1.2
4	−375.804	−1.2747	−1.3993	−2.6746	−0.1852	−0.8075	$4.508\times {10}^{-3}$	7.7
5	−586.382	−1.7362	−1.6181	−2.7142	−0.1694	−0.7751	$3.7891\times {10}^{-4}$	53.9
6	−781.235	−2.2512	−1.8266	−2.9008	−0.1522	−0.5291	$4.2399\times {10}^{-5}$	431.1

Table 5.

Computations for $f″(0)$ , $-g^{\prime }(0)$ , $-m^{\prime }(0)$ , and $-\Gamma \prime (0)$ with $\lambda =2$ , $\delta =1$ , $P_r=0.1$ , ${\delta }_c=0.1$ , ${\delta }_o=2$ , ${\Re }_d=1$ , ${\aleph }_{em}=0.2$ , ${\aleph }_b=0.01$ , ${\aleph }_c=0.1$ , ${\aleph }_d=0.5$ , ${\aleph }_a=0.5$ , $\sigma =1$ and various values of $\vartheta$ .

$\vartheta$	HAM result					Numerical result
	$f″(0)$	$-g^{\prime }(0)$	$-m^{\prime }(0)$	$-n^{\prime }(0)$	$-\Gamma \prime (0)$	$f″(0)$	$-g^{\prime }(0)$	$-m^{\prime }(0)$	$-n^{\prime }(0)$	$-\Gamma \prime (0)$
0.1	3.0396	0.1456	0.9274	0.9334	0.9998	3.0657	0.1458	0.9276	0.9335	0.9998
0.5	3.0111	0.0298	0.9274	0.9334	0.9998	3.0111	0.0298	0.9274	0.9333	0.9998
1	3.0054	0.0149	0.9274	0.9334	0.9998	3.0054	0.0150	0.9274	0.9333	0.9998
1.5	3.0036	0.0099	0.9274	0.9224	0.9998	3.0036	0.0100	0.9274	0.9223	0.9998
2	3.0027	0.0075	0.9274	0.9224	0.9999	3.0027	0.0075	0.9274	0.9223	0.9998

HAM: homotopy analysis method.

Table 6.

Computations for $f″(0)$ , $-g^{\prime }(0)$ , $-m^{\prime }(0)$ , and $-\Gamma \prime (0)$ with $\lambda =2$ , $\delta =1$ , $P_r=0.1$ , ${\delta }_c=0.1$ , ${\delta }_o=2$ , ${\Re }_d=1$ , ${\aleph }_b=0.01$ , ${\aleph }_c=0.2$ , ${\aleph }_d=0.5$ , ${\aleph }_a=0.05$ , $\vartheta =0.5$ and $\sigma =1$ and various values of ${\aleph }_{em}$ .

${\aleph }_{em}$	HAM result					Numerical result
	$f″(0)$	$-g^{\prime }(0)$	$-m^{\prime }(0)$	$-n^{\prime }(0)$	$-\Gamma \prime (0)$	$f″(0)$	$-g^{\prime }(0)$	$-m^{\prime }(0)$	$-n^{\prime }(0)$	$-\Gamma \prime (0)$
0.5	3.0105	0.0298	0.8319	0.8445	0.9998	3.0105	0.0298	0.8320	0.8445	0.9998
1	3.0013	0.0298	0.7012	0.7206	0.9998	3.0096	0.0298	0.7013	0.7207	0.9998
1.5	3.0101	0.0298	0.5976	0.6204	0.9998	3.0091	0.0298	0.5976	0.6204	0.9998
2	3.0085	0.0298	0.5139	0.5381	0.9998	3.0086	0.0298	0.5139	0.5381	0.9998
2.5	3.0082	0.0298	0.4453	0.4698	0.9998	3.0082	0.0298	0.4453	0.4698	0.9998

HAM: homotopy analysis method.

Table 7.

Computations for $f″(0)$ , $g^{\prime }(0)$ , $m^{\prime }(0)$ , and ${\vartheta }^{\prime }(0)$ with $\lambda =2$ , $\delta =1$ , $P_r=1$ , ${\delta }_c=2$ , ${\delta }_o=2$ , ${\Re }_d=1$ , ${\aleph }_{em}=2$ , ${\aleph }_b=0.1$ , ${\aleph }_d=0.5$ , ${\aleph }_a=1$ , $\vartheta =2$ , $\sigma =1$ and various values of ${\aleph }_c$ .

${\aleph }_c$	HAM result					Numerical result
	$f″(0)$	$-g^{\prime }(0)$	$-m^{\prime }(0)$	$-n^{\prime }(0)$	$-\Gamma \prime (0)$	$f″(0)$	$-g^{\prime }(0)$	$-m^{\prime }(0)$	$-n^{\prime }(0)$	$-\Gamma \prime (0)$
0.1	3.0734	0.0739	0.5150	0.5393	1.0133	3.0734	0.0740	0.5150	0.5394	1.0133
0.5	3.0735	0.0739	0.5150	0.5392	1.0133	3.0733	0.0739	0.5150	0.5392	1.0133
1	3.0743	0.0738	0.5150	0.5391	1.0133	3.0748	0.0738	0.5150	0.5390	1.0133
1.5	3.0760	0.0737	0.5150	0.5391	1.0133	3.0750	0.0738	0.5150	0.5388	1.0133
2	3.0774	0.0736	0.5150	0.5390	1.0133	3.0770	0.0737	0.5150	0.5386	1.0133

HAM: homotopy analysis method.

Table 8.

Fluid pressure and torques on upper and lower discs with $\lambda =2$ , $\delta =1$ , $P_r=0.1$ , ${\delta }_c=1$ , ${\delta }_o=0.2$ , ${\Re }_d=1$ , ${\aleph }_{em}=0.2$ , ${\aleph }_b=0.1$ , ${\aleph }_c=0.1$ , ${\aleph }_d=0.5$ , $\sigma =1$ and various values of $\vartheta$ .

$\vartheta$	HAM solution			Numerical solution
	$F_{\text{pres}}$	${\rho }^{\prime }(1)$	${\rho }^{\prime }(0)$	$F_{\text{pres}}$	${\rho }^{\prime }(1)$	${\rho }^{\prime }(0)$
0.1	−90.6705	0.0335	−0.0338	−90.6706	0.0334	−0.0339
0.5	−90.5757	0.0343	−0.0340	−90.5757	0.0343	−0.0341
1	−90.2796	0.0354	−0.0343	−90.2796	0.0355	−0.0343
1.5	−89.7862	0.0365	−0.3457	−89.7862	0.0365	−0.0345

HAM: homotopy analysis method.

Table 9.

Fluid pressure and torques on upper and lower discs with $\lambda =0.2$ , $\delta =1$ , $P_r=1$ , ${\delta }_c=0.1$ , ${\delta }_o=2$ , ${\Re }_d=5$ , ${\aleph }_{em}=0.2$ , ${\aleph }_c=1.5$ , ${\aleph }_d=0.5$ , ${\aleph }_a=0.5$ , $\sigma =0.3$ , $\vartheta =1.5$ and various values of ${\aleph }_b$ .

${\aleph }_b$	HAM solution			Numerical solution
	$F_{pres}$	${\rho }^{\prime }(1)$	${\rho }^{\prime }(0)$	$F_{pres}$	${\rho }^{\prime }(1)$	${\rho }^{\prime }(0)$
0.1	−192.47	−0.7583	−0.6939	−192.47	−0.7583	−0.6939
0.5	−18.01	−0.9779	−0.6655	−18.01	−0.9779	−0.6655
1	−9.49	−1.2845	−0.6185	−9.49	−1.2845	−0.6185
1.5	−9.21	−1.5982	−0.5983	−9.21	−1.5982	−0.5983

HAM: homotopy analysis method.

Figure 2 illustrates the maximum average residual error at different orders of approximation for $f(\eta )$ , $g(\eta )$ , $m(\eta )$ , $n(\eta )$ , $\Gamma (\eta )$ , and $\Psi (\eta )$ . It is clear from subfigures that error is almost continuously reduced up to $28\text{th}$ order of approximation. Table 1 presents the total residual error for different order of approximations taking fixed values of $\lambda =2$ , $\delta =1$ , $P_r=0.1$ , ${\delta }_c=1$ , ${\delta }_o=0.2$ , ${\Re }_d=1$ , ${\aleph }_{em}=0.2$ , ${\aleph }_b=0.1$ , ${\aleph }_c=0.1$ , ${\aleph }_d=0.5$ , ${\aleph }_a=0.05$ , $\vartheta =2$ , and $\sigma =1$ . Table 2 shows comparison of HAM results to numerical values of $f(\eta )$ , $g(\eta )$ , $m(\eta )$ , $n(\eta )$ , $\Gamma (\eta )$ , and $\Psi (\eta )$ with $\lambda =2$ , $\delta =1$ , $P_r=0.1$ , ${\delta }_c=1$ , ${\delta }_o=0.2$ , ${\Re }_d=1$ , ${\aleph }_{em}=0.2$ , ${\aleph }_b=0.1$ , ${\aleph }_c=0.1$ , ${\aleph }_d=0.5$ , ${\aleph }_a=0.05$ , $\vartheta =2$ and $\sigma =1$ and various values of $\eta$ . Table 4 presents optimal values of convergence control parameters as well as the maximum values of total average squared residual error versus different order of approximation. Here, it is noticed that the solution obtained from momentum, energy and transport equations converges to exact solution as we increase the order of approximation. More justification of our accurate solution is supported with the help of Table 3 for numerical values of $f″(0)$ , $g\prime (0)$ , $m\prime (0)$ , $n\prime (0)$ , $\Gamma \prime (0)$ , and $\Psi \prime (0)$ . Here, it is clear that the solutions almost converge at $10\text{th}$ order of approximations. Tables 5 to 9 compare the numerical values of HAM and BVP $4\text{c}$ for different values of $\vartheta ,{\aleph }_a,{\aleph }_b,{\aleph }_{em}$ and fixed values of remaining involved parameters.

Results and discussions

The mathematical formulation for the constitutive expressions of unsteady Newtonian fluid is employed to model the flow between the circular space of porous and squeezing discs in the form of equations (8) to (13) subject to the boundary conditions given in equation (14). These equations are solved and compared for numerical investigations through HAM and BVP $4\text{c}$ . HAM is an analytical method that gives solution in series form while BVP4c package is a numerical solver which is an adoptive mesh method that adopts nonuniform meshes for setting error in each step size. BVP $4\text{c}$ is described as collocation codes which solves BVPs by computing a cubic spline on each subinterval of a mesh of given interval, so the method can be viewed as a collocation method or a finite difference method with continuous extension. Parametric analyses are carried out for the dimensionless parameters such as MFD viscosity parameter $\vartheta$ , magnetic Reynolds number ${\aleph }_{em}$ , rotational Reynolds number ${\aleph }_a$ , squeezing Reynolds number ${\aleph }_b$ , axial magnetic force parameter ${\aleph }_c$ , and tangential magnetic force parameter ${\aleph }_d$ . The influence of these flow parameters is depicted both graphically and numerically for velocity components $f\prime (\eta )$ , $g(\eta )$ , $f(\eta )$ ; magnetic field components $m(\eta )$ , $n(\eta )$ ; temperature variation $\Gamma (\eta )$ ; and mass transport variation $\Psi (\eta )$ . Representative values are used to simulate physically realistic flows.

Figure 3(a) and (b) illustrates the effect of rotational Reynolds number ${\aleph }_a$ on $f\prime (\eta )$ and $g(\eta )$ . An increase in ${\aleph }_a$ means increase in angular velocity of lower disc or decrease in fluid viscosity. As expected, increase in angular velocity of lower disc is pushing the fluid in radial direction due to which $f\prime (\eta )$ is increasing as shown in Figure 3(a). It is also noticed that $f\prime (\eta )$ is gradually decreasing as it crosses the central region. The same but opposite behavior of the fluid is seen for azimuthal velocity $g(\eta )$ . The effect of rotational Reynolds number is also illustrated for the magnetic field components $m(\eta )$ and $n(\eta )$ in Figure 4(a) and (b). It is seen from this figure that increasing angular velocity of lower disc gives strength to both azimuthal and axial components of magnetic field. Figure 5 displays the influence of magnetic Reynolds number on the axial and azimuthal induced magnetic field distribution. It is obvious from the figure that increasing the ratio of the fluid flux to the magnetic diffusivity is increasing azimuthal induced magnetic field. However, a slight increase is seen for $n(\eta )$ which decreases near the upper disc.

Figure 3.

Profile of (a) $f^{\prime }(\eta )$ and (b) $g(\eta )$ for different values of ${\aleph }_a$ and fixed values of $\delta =2$ , $\sigma =1.3$ , $\lambda =1$ , $P_r=2$ , ${\delta }_c=0.5$ , ${\delta }_o=0.5$ , ${\Re }_d=5$ , ${\aleph }_{em}=0.1$ , ${\aleph }_b=0.01$ , ${\aleph }_c=7$ , ${\aleph }_d=3$ , and $\vartheta =1$ .

Figure 4.

Profile of (a) $m(\eta )$ and (b) $m(\eta )$ for different values of ${\aleph }_a$ and fixed values of $\delta =2$ , $\sigma =1.3$ , $\lambda =1$ , $P_r=2$ , ${\delta }_c=0.5$ , ${\delta }_o=0.5$ , ${\Re }_d=5$ , ${\aleph }_{em}=0.1$ , ${\aleph }_b=0.01$ , ${\aleph }_c=7$ , ${\aleph }_d=3$ , and $\vartheta =1$ .

Figure 5.

Profile of (a) $m(\eta )$ and (b) $n(\eta )$ for different values of ${\aleph }_{em}$ and fixed values of $\delta =1$ , $\sigma =1$ , $\lambda =0.1$ , $P_r=2$ , ${\delta }_c=2$ , ${\delta }_o=0.5$ , ${\Re }_d=5$ , ${\aleph }_b=0.5$ , ${\aleph }_c=0.1$ , ${\aleph }_d=2$ , ${\aleph }_a=3$ , and $\vartheta =1$ .

The effect of ${\aleph }_{em}$ is illustrated for the $g(\eta )$ and $\Gamma (\eta )$ in Figure 6. It is clear from this figure that as fluid flux is getting higher, fluid has no effect on azimuthal velocity distribution, but as fluid moves toward upper disc its velocity decreases. Similarly, fluid temperature is decreasing with increase of fluid flux or decrease of magnetic diffusivity. Figures 7 and 8 investigate the effect of $\sigma$ on $g(\eta )$ , $f\prime (\eta )$ , and $m(\eta )$ , respectively. The relative angular velocity parameter $\sigma =\frac{{\Omega }_u}{{\Omega }_l}$ represents the rotations of discs in same direction $(\sigma >0)$ and in opposite direction $(\sigma <0)$ . By choosing suitable values of $\sigma$ in the range $-1\le \sigma \le 1$ , it is possible to investigate the different flow regimes in fluid domain. The influence of $\sigma$ on $g(\eta )$ is shown in Figure 7(a). It is noticed that rotation of discs in same direction is increasing $g(\eta )$ near upper disc. It is also noticed that as fluid moves toward upper disc $g(\eta )$ is gradually increasing. Rotation of discs in opposite direction has opposite behavior on $g(\eta )$ as shown in Figure 7(b); similarly Figure 8 shows that axial magnetic field is increasing by rotation of discs in same direction and vice versa during opposite rotations.

Figure 6.

Profile of (a) $g(\eta )$ and (b) $\Gamma (\eta )$ for different values of ${\aleph }_{em}$ and fixed values of $\delta =1$ , $\sigma =1$ , $\lambda =0.1$ , $P_r=2$ , ${\delta }_c=2$ , ${\delta }_o=0.5$ , ${\Re }_d=5$ , ${\aleph }_b=0.5$ , ${\aleph }_c=0.1$ , ${\aleph }_d=2$ , ${\aleph }_a=3$ , and $\vartheta =1$ .

Figure 7.

(a and b) Profile of $g(\eta )$ for different values of $\sigma$ and fixed values of $\delta =1$ , $\lambda =1$ , $P_r=1$ , ${\delta }_c=1$ , ${\delta }_o=0.5$ , ${\Re }_d=0.5$ , ${\aleph }_{em}=3$ , ${\aleph }_b=1$ , ${\aleph }_c=2$ , ${\aleph }_d=7$ , ${\aleph }_a=1$ , and $\vartheta =1$ .

Figure 8.

(a and b) Profile of $m(\eta )$ for different values of $\sigma$ and fixed values of $\delta =1$ , $\lambda =1$ , $P_r=1$ , ${\delta }_c=1$ , ${\delta }_o=0.5$ , ${\Re }_d=0.5$ , ${\aleph }_{em}=3$ , ${\aleph }_b=1$ , ${\aleph }_c=2$ , ${\aleph }_d=7$ , ${\aleph }_a=1$ , and $\vartheta =1$ .

Figures 9 and 10 study the influence of MFD viscosity parameter $\vartheta$ on $f\prime \eta )$ , $g(\eta )$ , $m(\eta )$ , and $n(\eta )$ . As fluid viscosity is getting higher, the radial velocity starts increasing, but when fluid crosses the mid-domain it starts decreasing. It is also noticed that increase in viscosity is decreasing azimuthal velocity. Maximum decrease is seen in mid-domain of fluid. Figure 11 depicts the influence of dimensionless azimuthal magnetic force parameter ${\aleph }_c$ on radial velocity and axial induced magnetic field distribution. The influence of ${\aleph }_c$ will be dominant in the vicinity of upper disc where its angular velocity ${\Omega }_u$ is less than the angular velocity of the lower disc ${\Omega }_l$ . The radial velocity is accelerating up to mid-domain with acceleration of ${\aleph }_c$ as shown in Figure 11(a), whereas $m(\eta )$ is decelerating. Similarly, Figure 12(a) indicates that $f(\eta )$ is getting higher with increase of ${\aleph }_d$ . Hence, ${\aleph }_d$ should be used to increase the axial velocity. Figure 13 is made for the H-curves and Figure 14 shows the effect of ${\aleph }_c$ on $m\prime (\mathrm{\eta })$ . Tables 5 to 9 numerically investigate the current problem for different values of involved parameters. All results have excellent agreement with the results obtained via BVP $4\text{c}$ . The variation of load or pressure and torques that the fluid exerts on upper and lower discs are shown in Tables 8 and 9. The normal motion of upper disc and rotational motion are set to be the same, that is, ${\aleph }_a={\aleph }_b=1$ . We noticed that the load or pressure $F_{\text{pres}}$ is always negative and decreases for both $\vartheta$ and ${\aleph }_b$ . Also, increase in $\vartheta$ is increasing the torque on upper disc. Table 5 numerically investigates the influence of $\vartheta$ for $f″(0)$ and $-g\prime (0)$ . It is clear from the table that increase in viscosity decreases both $f″(0)$ and $-g\prime (0)$ . Similarly, Table 6 depicts the influence of magnetic Reynolds number ${\aleph }_{em}$ which decrease both $-m\prime (0)$ and $-n\prime (0)$ .

Figure 9.

Profile of (a) $f^{\prime }(\eta )$ and (b) $g(\eta )$ for different values of $\vartheta$ and fixed values of $\delta =2$ , $\sigma =1$ , $\lambda =0.1$ , $P_r=2$ , ${\delta }_c=0.5$ , ${\delta }_o=0.5$ , ${\Re }_d=5$ , ${\aleph }_{em}=1$ , ${\aleph }_b=0.02$ , ${\aleph }_c=1$ , ${\aleph }_d=2$ , and ${\aleph }_a=0.2$ .

Figure 10.

Profile of (a) $m(\eta )$ and (b) $n(\eta )$ for different values of $\vartheta$ and fixed values of $\delta =2$ , $\sigma =1$ , $\lambda =0.1$ , $P_r=2$ , ${\delta }_c=0.5$ , ${\delta }_o=0.5$ , ${\Re }_d=5$ , ${\aleph }_{em}=1$ , ${\aleph }_b=0.02$ , ${\aleph }_c=1$ , ${\aleph }_d=2$ , and ${\aleph }_a=0.2$ .

Figure 11.

Profile of (a) $f^{\prime }(\eta )$ and (b) $m(\eta )$ for different values of ${\aleph }_c$ and fixed values of $\delta =1$ , $\sigma =1$ , $\lambda =0.1$ , $P_r=0.2$ , ${\delta }_c=0.2$ , ${\delta }_o=0.2$ , ${\Re }_d=5$ , ${\aleph }_{em}=3$ , ${\aleph }_b=1.2$ , ${\aleph }_d=0.02$ , ${\aleph }_a=0.2$ , and $\vartheta =1$ .

Figure 12.

Profile of (a) $f(\eta )$ and (b) $g(\eta )$ for different values of ${\aleph }_d$ and fixed values of $\delta =2$ , $\sigma =1.3$ , $\lambda =1$ , $P_r=2$ , ${\delta }_c=0.5$ , ${\delta }_o=0.5$ , ${\Re }_d=5$ , ${\aleph }_{em}=0.01$ , ${\aleph }_b=0.01$ , ${\aleph }_c=0.7$ , ${\aleph }_a=1$ , and $\vartheta =1$ .

Figure 13.

H-curves profile of $f″(0)$ , $g^{\prime }(0)$ , $m^{\prime }(0)$ , $n^{\prime }(0)$ , $\Gamma \prime (0)$ , and $\Psi \prime (0)$ for fixed values of $\lambda =0.01$ , $\delta =0.01$ , $P_r=0.05$ , ${\delta }_c=0.4$ , ${\delta }_o=0.4$ , ${\Re }_d=0.2$ , ${\aleph }_{em}=0.01$ , ${\aleph }_b=0.05$ , ${\aleph }_c=0.03$ , ${\aleph }_d=0.01$ , ${\aleph }_a=0.03$ , $\vartheta =0.05$ , and $\sigma =0.1$ .

Figure 14.

Profile of $m^{\prime }(\eta )$ for different values of ${\aleph }_c$ and fixed values of $\delta =1$ , $\sigma =0.1$ , $\lambda =0.1$ , $P_r=0.2$ , ${\delta }_c=0.2$ , ${\delta }_o=0.2$ , ${\Re }_d=5$ , ${\aleph }_{em}=3$ , ${\aleph }_b=1.2$ , ${\aleph }_d=0.02$ , ${\aleph }_a=0.2$ , and $\vartheta =1$ .

Concluding remarks

In this study, the MFD viscosity of a squeezing flow for viscous fluid under the influence of heat/mass transfers and externally applied variable magnetic field has been investigated. HAM is used to determine the series solution of fluid velocity components $f(\eta )$ , $f\prime (\eta )$ , $g(\eta )$ ; magnetic field components $m(\eta )$ , $n(\eta )$ ; and mass and temperature distribution $\Gamma (\eta )$ , $\Psi (\eta )$ . Exceptional stability and convergence characteristics have been demonstrated with HAM. The physical key parameters emerging have been investigated graphically in detail including dimensionless $\sigma$ , ${\aleph }_d$ , ${\aleph }_{em}$ , ${\aleph }_c$ , ${\aleph }_a$ , and $\vartheta$ .

Main upshots of this article are presented as below:

The induced magnetic field distribution $m(\eta )$ and $n(\eta )$ are increasing functions of rotational Reynolds number. Further, the radial and azimuthal velocities are directly affected by an increase in ${\aleph }_a$ .

Magnetic Reynolds number causes an increase in magnetic field distributions and decrease in tangential velocity.

Fluid temperature decreases by increase in the magnetic Reynolds number.

The azimuthal and axial components of magnetic field have opposite behavior with increase in MFD viscosity.

It is observed from Table 5 that increase in MFD viscosity decreases $-f″(0)$ and $-g\prime (0)$ .

It is clear from Table 6 that increase in magnetic Reynolds number decreases both $-m\prime (0)$ and $-n\prime (0)$ .