*3.1. Procedure*

This section is devoted to obtaining the numerical solutions of resulting Equations (23) and (24) by means of the Runge–Kutta method with a shooting technique [38] using MATLAB software. This technique is preferred for two following reasons; firstly, the thermal energy equation is nonlinearly coupled and, secondly, due to missing of condition. For missing condition, first of all, we are supposed to convert the given system of differential equations into first-order initial value problems in order to carry out systematic guessing of missing initial conditions which will continue until the desired accuracy and convergence are not achieved. The following equations explain in detail the numerical process, which is the prerequisite for the adopted numerical technique.

$$
\mu\_f = \lg\_1 \tag{31}
$$

in which of *uf* is the velocity of the fluid phase. As it is an iterative scheme, in which each step has a possible error that can be successively reduced by changing higher order derivatives of *uf* , in terms of first-order ordinary differential equations as follows:

$$\frac{du\_f}{dy} = \mathbf{g}\_1' = \mathbf{g}\_2 \tag{32}$$

$$\frac{d^2 u\_f}{dy^2} = \mathcal{g}\_2' = \mathcal{g}\_3 \tag{33}$$

$$\frac{d^3 u\_f}{dy^3} = \mathbf{g}'\_3 = \mathbf{g}\_4 \tag{34}$$

where prime denoted the differentiation with respect to *y*. In view of Equations (31)–(34) the transformed form of Equation (23) is obtained as:

$$\mathbf{g}'\_{4} = \boldsymbol{\gamma}^{2}\mathbf{g}\mathbf{s} - \left(\frac{\boldsymbol{\gamma}^{2}M^{2}}{1-\mathbf{C}}\right)\mathbf{g}\_{1} - \left(\frac{\boldsymbol{\gamma}^{2}}{1-\mathbf{C}}\right)\boldsymbol{P} \tag{35}$$

Similar to the preceding pattern, one finds no absurdness to convert thermal differential Equation (24) into first-order system by making the following suppositions:

$$
\Theta = \lg\_5 \tag{36}
$$

$$\frac{d\Theta}{dy} = \mathbf{g}'\_5 = \mathbf{g}\_6 \tag{37}$$

By using Equations (36) and (37) in Equation (24), we have:

$$\mathcal{g}'\_{6} = \frac{B\_{r}}{m} \left(\mathcal{g}\_{2}\right) \left(\mathcal{g}\_{4}\right) - B\_{r} \left(\mathcal{g}\_{2}\right)^{2} \tag{38}$$

In view of Equations (35) and (38), the associated boundary conditions given in Equations (25)–(30) at the lower and upper plate can be obtained as:

$$\begin{aligned} (i) \ \text{g}\_1 &= \beta\_1 \Big(\text{g}\_2 - \frac{\text{g}\_4}{\text{y}^2}\Big) \\ (ii) \ \text{g}\_2 &= c\_1 \\ (iii) \ \text{g}\_3 &= 0 \\ (iv) \ \text{g}\_4 &= c\_2 \\ (v) \ \text{g}\_5 &= 0 \\ (vi) \ \text{g}\_6 &= c\_3 \end{aligned} \qquad \text{where } y = -1 \tag{39}$$

In the same way given and missing conditions at the upper wall are:

$$\begin{aligned} (i) \quad &g\_1 = -\beta\_1 \Big(g\_2 - \frac{g\_4}{\gamma^2}\Big) \\ (ii) \quad &g\_2 = c\_4 \\ (iii) \quad &g\_3 = 0 \\ (iv) \quad &g\_4 = c\_5 \\ (v) \quad &g\_5 = 1 \\ (vi) \quad &g\_6 = c\_6 \end{aligned} \qquad \text{where } y = 1 \tag{40}$$

where *c*1, *c*2, *c*3, *c*4, *c*5, and *c*<sup>6</sup> are the missing conditions which can be easily determined during the routine calculation.

#### *3.2. Validation*

The numerical results for both phases are presented in Tables 1–3. Table 1 offers the variation of velocity for both phases against the slip parameter. Table 2 shows the simultaneous variations in the velocities for single- and two-phase flows at different points. Table 3 displays the thermal variation at different points within the given domain when *M* = 0.5, γ = 2.0, and *Br* = 2.0. It is found that the results extracted numerically are compatible with the physical expectations and satisfy all the subjected conditions as shown graphically. This provides a useful check that the presented solutions are correct.


**Table 1.** Variation in the velocity of both phases against slip parameter.

**Table 2.** Variation in the velocities for single- and two-phase flows at different points.



**Table 3.** Thermal variation at the different points of given domain.

#### *3.3. Discussion*

In this section, a concise study of pertinent parameters is graphically presented in Figures 2–10. Figures 2 and 3 are plotted to examine the influence of the magnetic parameter on the motion of couple stress fluid and metallic Hafnium particles. In both graphs, a clear decline in velocities for higher values of the magnetic parameter is observed. Nevertheless, the theory of Hannes Alfven explains the same phenomenon involving the interaction of magnetic fields being induced into an electrically conducting fluid system. This phenomenon produces Alfven waves which result in clear retardation of fluid's speed. However, in Figures 4 and 5, the density of the Hafnium particles brings out a different result as compared to magnetic fields. The major push of pressure on the fluid on slippery walls, the hydro motion in both phases is supported by the addition of extra metallic particles. Consequently, velocity increases by increasing the number of particles. Such factors can be regarded as to attenuate the interaction of fluid particle or interparticle collision allowing the particles to move with least resistance. The most significant parameter which constitutes the existence of the present fluid flow is couple stress parameter γ. It is observed that the fluid particle additives, contribute to expediting the movement. This may cause obscurity and vagueness in the mind of a reader, but Equations (23), (28), and (31) provide enough clues about the inverse influence of couple stress parameter on the flow that attenuates the force of friction/drag arising from the effect base fluid's accumulation. This constitutes a size-dependent effect in the base fluid, in addition to minimizing the rotational field of the fluid particles. Hence, rapid fluid flow is observed in both Figures 6 and 7. However, the contribution of slippery walls is not negligible, as they assist the metallic particles to frisk freely in the liquid. The role of the slip parameter that supports the velocity of both phases is spotted in Figures 8 and 9. Generally, it is believed that slippery walls only snag the flow because of their behavior as a retarding force. Against all such expectations, in the present study, slip effects bring about unprecedented change by increasing the velocity of the fluid, as shown in Figure 8. This change is due to the inverse influence of γ, given in Equations (28) and (31) which rebuffs all such perception that slip parameter merely hampers the flow. The change in temperature through Brinkman number *Br* is sketched in Figure 10. It is revealed that the higher values of Brinkman heat up the fluid temperature.

**Figure 2.** Effects of magnetic fields on the flow.

**Figure 3.** Effects of magnetic fields on the particles.

**Figure 4.** Effects of *C* on the flow.

**Figure 5.** Effects of *C* on the particles.

**Figure 6.** Couples stress parameter affecting the flow.

**Figure 7.** Couples stress parameter affecting the motion of particles.

**Figure 8.** Effects of slip parameter on the flow.

**Figure 9.** Effects of slip parameter on the motion of particles.

**Figure 10.** Role of Brinkman number on the temperature.
