**3. Skin Friction and Local Nusselt Number**

Shear stresses at lower disk in radial and tangential directions are τ*zr* and τ*z*<sup>θ</sup>

$$\tau\_{zr} = \mu\_{nf} u\_{z}|\_{z=0} = \frac{\mu\_{f} r \Omega\_{1} f^{\prime\prime}(0)}{h(1-\phi)^{2.5}},\\ \tau\_{z0} = \mu\_{nf} u\_{z}|\_{z=0} = \frac{\mu\_{f} r \Omega\_{1} g^{\prime}(0)}{h(1-\phi)^{2.5}},\tag{37}$$

The total shear stress is

$$
\tau\_w = \left(\tau\_{zr}^2 + \tau\_{z0}^2\right)^{1/2},\tag{38}
$$

Coefficients of drag force at *z* = 0, and *z* = *h* for the disk are

$$\begin{split} \mathbb{C}\_{f\_1} &= \frac{\tau\_w \underset{\rho\_f(r\Omega\_1)^2}{\operatorname{Re}\_f(r\Omega\_1)^2}}{\rho\_f(r\Omega\_1)^2} = \frac{1}{\operatorname{Re}\_r(1-\phi)^{25}} \Big[ \left(f''(0)\right)^2 + \left(g'(0)\right)^2 \Big]^{1/2}, \\ \mathbb{C}\_{f\_2} &= \frac{\tau\_w \underset{\rho\_f(r\Omega\_2)^2}{\operatorname{Re}\_f(r\Omega\_2)^2}}{\rho\_f(r\Omega\_2)^2} = \frac{1}{\operatorname{Re}\_r(1-\phi)^{25}} \Big[ \left(f''(1)\right)^2 + \left(g'(1)\right)^2 \Big]^{1/2}, \end{split} \tag{39}$$

Here, *Rer* represents local Reynolds number.

The dimensional form of *Nu* (the local Nusselt number) is

$$Nu = \frac{k\_{nf}(\rho c\_p)\_f}{\rho\_f k\_f},\tag{40}$$

By using transformation given in Equations (25), Equation (40) becomes

$$(1 - ct)^{1/2} \text{Nu}\_1 = -\frac{k\_{nf}}{k\_f} \theta'(0)\_\prime (1 - ct)^{1/2} \text{Nu}\_2 = -\frac{k\_{nf}}{k\_f} \theta'(1)\_\prime \tag{41}$$

### **4. Numerical Method**

In current model, MATLAB built-in-function bvp4c is used to solve coupled ordinary differential equations (ODE's) (Equations (26–36)) with mentioned boundary conditions (32). The computational purpose of the infinite domain is restricted to η = 4 which is enough to indicate the asymptotic behavior of the solution. The theme numerical scheme needs initial approximation with tolerance 10<sup>−</sup>6. The initial taken estimation must meet the boundary conditions without interrupting the solution technique. We obtain a system comprising three first-order differential equations given below:

$$\begin{aligned} f' &= y\_2, \\ f'' &= y\_3, \\ f''' &= y\_4 \\\ g' &= \frac{y}{A}(A\_1(\frac{1}{2}y\_3 + \frac{\eta}{2}y\_4) + \text{Re}(2y\_1y\_4 - 2y\_5y\_6) - \frac{\sigma\_{\eta f}}{\sigma\_f}\frac{\text{MRc}(y\_5 - my\_6)}{\text{B}(1 + \text{m}^2)}), \\ g' &= y\_5 \\ g' &= y\_6, \\ y\_2 &= \frac{B}{R}R[\left(y\_5 + \frac{1}{2}\eta y\_6\right)A\_1 + 2\left(y\_2y\_5 - y\_1y\_6\right)] - \frac{\sigma\_{\eta f}}{\sigma\_f}\frac{\text{MRc}(y\_5 + \eta y\_2)}{A(1 + \text{m}^2)}, \\ \theta &= y\_7, \\ \theta &= y\_8 \\ y\_3 &= \frac{1}{\frac{\eta}{R}\frac{1}{\eta\eta\eta} - \frac{1}{A\_1}(y\_1)^2 - 2\eta y\_1}(A\_1\{s+y\_7+\frac{1}{2}\eta y\_8\} + (s+y\_7)y\_2 - 2y\_1y\_8 + \\ \gamma\|\{s+y+\frac{\eta}{2}\eta y\_8\} + y\_2\{y\_2+\frac{1}{2}\eta y\_3\}(s+y\_7) + 2y\_2\{s+y\_7+\frac{1}{2}\eta y\_8\} - \\ 6y\_1y\_8 + (y\_1+\eta y\_2)\theta' + (y\_2)^2(s+y\_7) - \frac{\Lambda}{A\_1}y\_2y\_8 - 2y\_1y\_3(s+y\_7) - \frac{\eta}{\Omega\tau}(s+y\_7)), \end{aligned} \tag{42}$$

With suitable boundary condition

$$\begin{aligned} y\_1(0) &= 0, y\_2(0) = \gamma\_1, y\_5(0) = 1, y\_7(0) = 1 - \mathbf{s}, \\ y\_1(1) &= 0, y\_2(1) = \gamma\_2, y\_5(1) = \Omega, y\_7(1) = 0 \end{aligned} \tag{43}$$
