**3. Entropy Optimization and Bejan Number**

There are three causes of entropy Optimization in the current problem. Heat transport, viscous dissipation, and mass diffusion all generate entropy. The entropy equation is written as follows [18,21]:

$$\begin{split} S\_{\mathcal{S}} = \frac{k\_{nf}}{T\_{\infty}^{-2}} \Big( 1 + \frac{16\sigma^{\*}T^{3}}{3\overline{k}^{4}k\_{nf}} \Big) \left( \frac{\partial T}{\partial \overline{r}} \right)^{2} + \frac{1}{T\_{\infty}} \Big( \mu\_{nf} \left( \frac{\partial u}{\partial \overline{r}} \right)^{2} + \frac{1}{\overline{\rho}\overline{c}} \left( \frac{\partial u}{\partial \overline{r}} \right)^{2} - \frac{1}{6\overline{\rho}\overline{c}^{3}} \left( \frac{\partial u}{\partial \overline{r}} \right)^{4} \Big) \\ &+ \frac{RD}{T\_{\infty}} \Big( \frac{\partial T}{\partial \overline{r}} \frac{\partial C}{\partial \overline{r}} \Big) + \frac{RD}{T\_{\infty}} \left( \frac{\partial C}{\partial \overline{r}} \right)^{2} \end{split} \tag{16}$$

By using the typical entropy generation rate to convert Equation (16) to dimensionless form, we obtain:

$$\begin{split} \mathcal{N}\_{\mathbb{G}} &= \left( \frac{k\_{\mathbb{M}}}{k\_f} + \mathcal{R}d(\theta(\theta\_{\omega}-1)+1)^3 \right) \big( 1+2\eta\gamma\rangle a\_1 \theta'^2 + L\_1(1+2\eta\gamma)\theta'\phi' \\ &+ L\_1(1+2\eta\gamma)\frac{a\_2}{a\_1}\phi'^2 + \mathcal{B}r(1+2\eta\gamma) \left( \left( \frac{1}{\left(1-\phi\right)^{2.5}} + a \right) f''^2 - \frac{1}{3} a\lambda f''^4 \right) \end{split} \tag{17}$$

where *NG* <sup>=</sup> *Sg <sup>S</sup>*<sup>0</sup> <sup>=</sup> *Sg <sup>k</sup> <sup>f</sup> <sup>U</sup>*0(*Tω*−*T*∞)/*T*∞*μ<sup>f</sup> <sup>L</sup>* denotes the total entropy production, *<sup>α</sup>*<sup>1</sup> <sup>=</sup> *<sup>T</sup>ω*−*T*<sup>∞</sup> *T*∞ and *α*<sup>2</sup> = *<sup>C</sup> <sup>ω</sup>*−*C*<sup>∞</sup> *<sup>C</sup>*<sup>∞</sup> are the temperature ratio and concentration ratio variables, respectively, and *Br* <sup>=</sup> *<sup>μ</sup><sup>f</sup> <sup>U</sup>*<sup>0</sup> <sup>2</sup>*z*<sup>2</sup> *<sup>k</sup> <sup>f</sup> <sup>L</sup>*2(*Tω*−*T*∞) and *<sup>L</sup>*<sup>1</sup> <sup>=</sup> *RD*(*Cω*−*C*∞) *k f* are Brinkman number and diffusion parameter, respectively.

Equation (17) can be written using the pattern shown below:

$$N\_{\rm G} = N\_{\rm H} + N\_{\rm f} + N\_{\rm m} \tag{18}$$

Here, *NH* represents entropy production because of heat transmission, *Nf* represents entropy production because of fluid friction, and *Nm* represents entropy production. The Bejan number is defined as:

$$Be = \frac{N\_H + N\_m}{N\_G} \tag{19}$$

$$B\varepsilon = \frac{\begin{pmatrix} \frac{k\_{nf}}{k\_f} + R(\theta(\theta\_\omega - 1) + 1)^3 \Big) (1 + 2\gamma\eta) \theta'^2 a\_1 + L(1 + 2\gamma\eta) \theta' \phi' + L(1 + 2\gamma\eta) \frac{a\_1}{a\_1} \phi'^2\\ \left(\frac{k\_{nf}}{k\_f} + R d(\theta(\theta\_\omega - 1) + 1)^3 \right) (1 + 2\eta\gamma) a\_1 \theta'^2 + R\gamma(1 + 2\eta\gamma) \left(\left(\frac{1}{(1 - \phi\_1)^{2.5}} + a\right) f^{\eta 2} - \frac{1}{3} a\lambda f^{\rho 4}\right) \end{pmatrix}}{+ L\_1(1 + 2\eta\gamma) \theta' \phi' + L\_1(1 + 2\eta\gamma) \frac{a\_1}{a\_1} \phi'^2} \tag{20}$$

#### **4. Numerical Method**

The solution mechanism for the currently constructed model is computed in this part manipulating the bvp4c technique (shooting scheme). The bvp4c technique (shooting scheme) in the MATLAB tool is used to solve the ODEs (7)–(11) via (12) numerically. First, we transform a higher-order system to a first-order system for this technique. To accomplish this, we follow the steps below:

$$\begin{array}{l} f = y\_{1\prime} f' = y\_{2\prime} f'' = y\_{3\prime} f''' = y\_3' \\ \theta = y\_{4\prime} \theta' = y\_{5\prime} \theta'' = y\_5' \\ \phi = y\_6, \phi' = y\_{7\prime} \phi'' = y\_7' \end{array} \tag{21}$$

$$\begin{split} y\_{3}^{\prime} &= \frac{\left( \begin{array}{c} \frac{4}{3}a\lambda\gamma(1+2\eta\gamma)y\_{3}^{-3} - 2\eta\left(\frac{1}{(1-\Phi\_{1})^{15}} + a\right)y\_{3} \\ + \frac{6s}{(1-\Phi\_{1})^{25}}y\_{2} + F\_{\eta}y\_{2}^{-2} - \left( (1-\Phi\_{1}) + \Phi\_{1}\frac{\rho\_{\star}}{\rho\_{f}} \right) (y\_{1}y\_{3} - y\_{2})^{2} \end{array} \right)}{\left( (1+2\eta\gamma)\left(\frac{1}{(1-\Phi\_{1})^{25}} + a\right) - \lambda a (1+2\eta\gamma)^{2}y\_{3} \right)^{2}} \\ & \begin{cases} Ec(1+2\eta\gamma) \left( \frac{1}{3}\lambda a (1+2\eta\gamma)y\_{3}^{-4} - \left( \frac{1}{(1-\Phi\_{1})^{25}} + a \right) y\_{3}^{2} \right) - \frac{k\_{\eta}/k\_{f}}{P}\gamma y\_{5} \\ - \left( (1-\Phi\_{1}) + \Phi\_{1}\frac{\left(\rho\_{\star}\right)\_{\star}}{\left(\rho\_{\star}\right)\_{f}} \right) \left( N(1+2\eta\gamma)y\_{5}^{2} + Nb(1+2\eta\gamma)y\_{5}y\_{7} + y\_{1}y\_{5} \right) \end{cases} \\ y\_{5}^{\prime} &= \frac{-\frac{3\eta R}{P\kappa\_{\mu}/k\_{f}} (1+2\eta\gamma) \left( y\_{4}(\theta\_{\omega}-1) + 1 \right)^{2} (\theta\_{\omega}-1)y\_{5}^{2}}{\frac{k\_{\eta}/k\_{f}}{P\eta} \left( 1+2\eta\gamma\right) + \frac{R}{P\kappa\_{\mu}/$$

with

$$\begin{aligned} y\_1(0) &= 0, y\_2(0) = 1, y\_4(0) = 1, y\_6(0) = 1 \\ y\_2(\infty) &\to 0, y\_4(\infty) \to 0, y\_6(\infty) \to 0 \end{aligned} \tag{25}$$
