A Fully Resolved Computational Fluid Dynamics Study of the Boundary Layer Flow of an Aqueous Nanoliquid Comprising Gyrotactic Microorganisms over a Stretching Sheet: The Validity of Conventional Similarity Models

Abstract

When materials are processed in the form of sheets that are stretched, cooling is often required. Coolants have been developed to maximize the rate of heat transfer away from the sheet, including by adding nanoparticles and microorganisms to control the physical properties of the fluid. Such coolants perform well, but the interaction between them and the sheet is not yet fully understood. Most of the articles found in the literature have used similarity models to solve the set of governing equations. In this method, the governing equations can be mapped into a set of 1-D equations and solved easily. However, care should be taken when using this method as the validity of this method is ensured only in the fully developed region, far away enough from the extrusion slit. The present study, therefore, aims to explore the reliability of a similarity model by comparing it with a full computational fluid dynamics (CFD) approach. In this work, the boundary layer flow of a nanoliquid comprising gyrotactic microorganisms in both the developed and undeveloped regions of a stretching sheet is studied using computational fluid dynamics with the finite difference approach, implemented using FORTRAN. The results of the CFD method are compared against the similarity analysis results for the length of the developed and undeveloped regions. This study, for the first time, distinguishes between the undeveloped and fully developed regions and finds the region in which the similarity analysis is valid. The numerical results show that the critical Reynolds numbers for the boundary layers of the concentration of the nano-additives and of density of the microorganisms are equal. To achieve an agreement between the CFD and the similarity model within 5%, the Grashof number for the hydrodynamic boundary layer must be <4 × 10⁴. Increasing the bioconvection Rayleigh number leads to a decrease in the skin friction coefficient. The length of the region in which the microorganism’s density is not fully developed remains approximately constant for 10³ < Gr < 10⁵. Nonetheless, this length reduces significantly when the Grashof number increases from 10⁵ to 10⁶. The reduced Nusselt number, Nu_r, increases when the density difference of the microorganisms increases.

1. Introduction

The large number of applications of stretching sheets in different industrial processes has motivated many researchers to focus on different physical phenomena related to such sheets. Several industrial and engineering implementations involve the cooling of metal sheets [1], paper and glass production [2,3], polymer processing [4], and the production of plastic sheets [5]. In such processes, the heat transfer rate between the stretching sheet and the working fluid plays a vital role in the quality of the manufactured products. To minimize process times and maximize productivity, it is often desirable to maximize the rate of heat transfer. Accordingly, different methods of cooling have been purposed such as injection (blowing), using a discontinuous stretching surface, and linear and nonlinear stretching velocities [6,7,8,9,10,11]. In addition, another approach could be altering the properties of the coolant fluid using single or hybrid nanoparticles, nano-encapsulated phase change materials, or by imposing magnetic (or electric) fields to magnetically (or electrically) conducting fluids [12,13,14,15,16].

Nanoparticles are widely used to enhance the thermal properties of working fluids with low thermal conductivities such as engine oil, ethylene glycol, and water [17,18] to improve the heat transfer rates in different processes such as engine cooling [19], thermal energy transport in heat exchangers [20,21], and the cooling of electronic devices [22]. Studies have been published with models for the thermophysical properties of various nanoliquids, many of which determine the effective dynamic viscosity and thermal conductivity from an empirical analysis or theory [23,24,25,26,27]. Publications have also reported that the use of nanoparticles with a volume fraction of $\le 5\%$ can increase the effective thermal conductivity of the nanoliquid by $>20\%$ [28,29,30,31,32].

Bioconvection is fluid motion caused by microscopic density gradients generated by the communal motion of swimming microorganisms [33,34,35,36]. The increase in the base fluid density due to the aggregation of motile microorganisms in a certain region causes bioconvection. Swimming microorganisms within the nanoliquid flow can amplify the mass transfer rate and the mixing, and prevent the instability of nanoparticles [37,38].

The study of the flow of nanoliquids containing microscopic organisms is not very old; the first study, published in 2010, surveyed the onset of bioconvection in a vertical flow arising from free convection [39]. Since then, other studies have investigated the boundary layer flows of the forced and free convections of nanoliquids containing microorganisms over stretching and static sheets [40,41,42,43]. Javad et al. [41] studied the flow of a fluid comprising nanoparticles and microorganisms over a porous stretching sheet. They found that the density of the microorganism is enhanced by raising the thermophoresis parameter and reduced by increasing the Schmidt number. Shahid et al. [42] numerically analyzed the characteristics of the flow of a nanofluid-containing gyrotactic over a stretching sheet in a porous medium. The results showed that the permeability parameter retards the motion of microorganisms. Sulaiman [43] reported that an increment in bioconvection Lewis and Peclet numbers in an Oldroyd-B nanofluid with gyrotactic microorganisms increases the concentration of the microorganisms. In all of these studies, a similarity analysis has been used for obtaining the flow behavior. However, no discussion on the undeveloped regions where this method is not valid could be found in literature [40,41,42,43].

Publications associated with the boundary layer flow of nanoliquids containing motile microorganisms on either stretching or static sheets have employed a similarity analysis to model the governing physics. There is no study investigating the flow physics of a nanoliquid containing mobile microorganisms in the entrance region (near the extrusion slit in the stretching sheet). The present study aims to explore the reliability of the mentioned method by comparing it with fully resolved CFD simulations. This study, for the first time, distinguishes between the undeveloped and fully developed regions and identifies the region in which the similarity analysis is valid. For this purpose, critical Reynolds and Grashof numbers are obtained for different locations along the sheet (different X values), defined as values for which the disagreement between the CFD and the similarity analysis is no more than 5%. The critical values define a range of applicability for the similarity solution method. In this paper, we present the Materials and Methods (including the problem formulation, numerical approach, mesh sensitivity analysis, and validity) in Section 2. The results and corresponding discussion are presented in detail in Section 3, and the conclusions of the study are presented in the final section.

2. Materials and Methods

2.1. Problem Formulation

A simple schematic view of the problem with applicable boundary conditions shows a sheet linearly stretched along the x-axis (Figure 1, where u and v are the components of velocity in the x- and y-direction, respectively; T is temperature, c is the concentration of nanoparticles and can take values between 0 and 1, n is the concentration of microorganisms, the subscript _w refers to the surface of the sheet, and the subscript $\infty$ refers to the quantities far from the sheet). In this case, the velocity is at the sheet, u_w = ax, where a is a constant. The temperature, T_w, the concentration of nanoparticles, c_w, and the density of the microorganisms on the sheet, n_w, are assumed to be constant. Far from the sheet, the velocity is zero and the temperature, nanoparticle concentration, and microorganism density are T_∞, c_∞, and n_∞, respectively.

To model the problem efficiently, simplifying assumptions are applied:

• The flow is incompressible, viscous, steady, and laminar.
• The sheet is impermeable.
• There is no slip on the sheet surface.
• The nanoliquid is dilute so the nanoparticles do not interact with each other or affect the motion of the microorganisms.
• The Boussinesq approximation may be applied to any change in density.
• The Brownian motion and thermophoresis are considered.

Applying the assumptions above leads to a continuity equation (Equation (1)) [38]; the conservation of momentum in the x-direction (Equation (2)) and y-direction (Equation (3)) [38], the conservation of energy (Equation (4)), the conservation of the nanoparticles (Equation (5)) [4,38], and the conservation of the microorganisms (Equation (6)) [38,44]. In Equations (2) and (3), ρ_f is the density of the coolant fluid, p is the pressure in the fluid, u and v are the components of velocity in the x- and y-direction, μ_f is the viscosity of the fluid, g is the acceleration due to gravity, β_f is the volumetric thermal expansion coefficient of the fluid, γ is the average volume of a microorganism, and ρ_p is the density of the nanoparticles. In Equation (4), T is the temperature, α_f is the thermal diffusivity of the fluid, D_B is the Brownian diffusion coefficient, and D_T is the thermophoresis diffusion coefficient. In Equation (5), c is the concentration of the nanoparticles and in Equation (6), n is the density of the microorganisms and D_m is the diffusion coefficient of the microorganisms.

$$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0.$$

(1)

$${\rho }_f\left(u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}\right)=-\frac{\partial p}{\partial x}+{\mu }_f\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right).$$

(2)

$${\rho }_f\left(u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}\right)=-\frac{\partial p}{\partial y}+{\mu }_f\left(\frac{{\partial }^{2}v}{\partial x^2}+\frac{{\partial }^{2}v}{\partial y^2}\right)+g\left\{\begin{array}{l}{\rho }_f{\beta }_f\left(1-c_{\infty }\right)\left(T-T_{\infty }\right)\\ -\left({\rho }_p-{\rho }_f\right)\left(c-c_{\infty }\right)\\ -\gamma \left({\rho }_m-{\rho }_f\right)\left(n-n_{\infty }\right)\end{array}\right\}.$$

(3)

$$u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}={\alpha }_f\left(\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}\right)+\tau \left\{\begin{array}{l}{D}_B\left(\frac{\partial c}{\partial x}\frac{\partial T}{\partial x}+\frac{\partial c}{\partial y}\frac{\partial T}{\partial y}\right)\\ +\frac{D_T}{T_{\infty }}\left({\left(\frac{\partial T}{\partial x}\right)}^2+{\left(\frac{\partial T}{\partial y}\right)}^2\right)\end{array}\right\}$$

(4)

where $\tau ={\left(\rho c_{\mathrm{p}}\right)}_{\mathrm{p}}/{\left(\rho c_{\mathrm{p}}\right)}_{\mathrm{f}}$.

$$u\frac{\partial c}{\partial x}+v\frac{\partial c}{\partial y}=D_B\left(\frac{{\partial }^{2}c}{\partial x^2}+\frac{{\partial }^{2}c}{\partial y^2}\right)+\frac{D_T}{T_{\infty }}\left(\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}\right).$$

(5)

The macroscopic motion of the microorganism can be expressed as:

$$\overrightarrow{\nabla }\cdot \left(n\overrightarrow{u}+n\overrightarrow{q}-D_m\nabla n\right)=0.$$

(6)

The first term in the parentheses shows the flux of the microorganisms due to macroscopic fluid motion, the second one represents the directional swimming of the microorganisms up the nanoparticle gradients, and the last term is the diffusive process that models all random motions of the microorganisms, respectively. Note that $\overrightarrow{u}=(u,v)$ is the velocity vector of the nanofluid and $\overrightarrow{q}$ is the average swimming velocity of the microorganism and can be approximated as [38]:

$$\overrightarrow{q}=b{W}_c\nabla c$$

(7)

where b has the dimensions of length and W_c is the maximum cell swimming velocity (the product bW_c is assumed to be constant). Substituting Equation (7) for Equation (6) and expanding them results in the following equation for the motion of the microorganism:

$$u\frac{\partial n}{\partial x}+v\frac{\partial n}{\partial y}+\frac{b{W}_c}{\left(c_w-c_{\infty }\right)}\left[\frac{\partial }{\partial y}\left(n\frac{\partial c}{\partial y}\right)+\frac{\partial }{\partial x}\left(n\frac{\partial c}{\partial x}\right)\right]=D_m\left(\frac{{\partial }^{2}n}{\partial x^2}+\frac{{\partial }^{2}n}{\partial y^2}+2\frac{{\partial }^{2}n}{\partial x\partial y}\right).$$

(8)

The boundary conditions applied in the dimensional coordinates lead to further equations that describe the flow (Equations (9)–(12)) [4] where H is the height of the computational domain and L is the total length of the sheet.

$$u=v=0,T=T_{\infty },c=c_{\infty },n=n_{\infty }atx=0 \& 0\le y\le H.$$

(9)

$$u=u_w=ax,v=0,T=T_w,c=c_w,n=n_{w}aty=0 \& 0\le x\le L.$$

(10)

$$u=0,\frac{\partial v}{\partial y}=0,T=T_{\infty },c=c_{\infty },n=n_{\infty }aty\to \infty \& 0\le x\le L.$$

(11)

$$\frac{\partial u}{\partial x}=0,\frac{\partial v}{\partial x}=0,\frac{\partial c}{\partial x}=0,\frac{\partial n}{\partial x}=0atx=L \& 0\le y\le H.$$

(12)

The dimensionless parameters, θ, C, N, X, Y, U, V, and P are introduced to convert Equations (1)–(6) and all boundary conditions to dimensionless coordinates (Equation (13)). The pressure terms in Equations (2) and (3) can be replaced by the stream function, $\Psi$, and the vorticity function, $\omega$ (Equations (14) and (15)).

$$\begin{array}{l}\theta =\frac{T-T_{\infty }}{T_w-T_{\infty }},C=\frac{c-c_{\infty }}{c_w-c_{\infty }},N=\frac{n-n_{\infty }}{n_w-n_{\infty }}\\ X=\frac{x}{L},Y=\frac{y}{L},U=\frac{u}{u_L},V=\frac{v}{u_L},P=\frac{p}{\rho u_L{}^2}\end{array}$$

(13)

$$U=\frac{\partial \Psi }{\partial Y},V=-\frac{\partial \Psi }{\partial X}.$$

(14)

$$\omega =\frac{\partial V}{\partial X}-\frac{\partial U}{\partial Y}.$$

(15)

Substituting Equations (13)–(15) into Equations (1)–(6) leads to dimensionless equations in terms of the stream and vorticity functions, Ψ and ω (Equations (16)–(20), where $Nb$ is the Brownian motion parameter, $Nt$ is the thermophoresis parameter, $Gr$ is the Grashof number (Equation (21)), $Nr$ is the buoyancy ratio parameter, $Re$ is the Rayleigh number, $Pr$ is the Prandtl number, $Le$ is the traditional Lewis number, $Lb$ is the bioconvection Lewis number (Equation (22)), $Pe$ is the Peclet number, $\Omega$ is the difference in the concentration of the microorganisms, and $Rb$ is the bioconvection Rayleigh number (Equation (23)).

$$\frac{{\partial }^2\Psi }{\partial X^2}+\frac{{\partial }^2\Psi }{\partial Y^2}=-\omega$$

(16)

$$U\frac{\partial \omega }{\partial X}+V\frac{\partial \omega }{\partial Y}=\frac{1}{Re}\left(\frac{{\partial }^2\omega }{\partial X^2}+\frac{{\partial }^2\omega }{\partial Y^2}\right)+\frac{Gr}{R{e}^2}\left(\frac{\partial \theta }{\partial X}-Nr\frac{\partial C}{\partial X}-Rb\frac{\partial N}{\partial X}\right).$$

(17)

$$\frac{\partial \Psi }{\partial Y}\frac{\partial \theta }{\partial X}-\frac{\partial \Psi }{\partial X}\frac{\partial \theta }{\partial Y}=\frac{1}{RePr}\left(\frac{{\partial }^2\theta }{\partial X^2}+\frac{{\partial }^2\theta }{\partial Y^2}\right)+\frac{Nb}{Re}\left(\begin{array}{l}\frac{\partial C}{\partial X}\frac{\partial \theta }{\partial X}\\ +\frac{\partial C}{\partial Y}\frac{\partial \theta }{\partial Y}\end{array}\right)+\frac{Nt}{Re}\left(\begin{array}{l}{\left(\frac{\partial \theta }{\partial X}\right)}^2\\ +{\left(\frac{\partial \theta }{\partial Y}\right)}^2\end{array}\right).$$

(18)

$$\frac{\partial \Psi }{\partial Y}\frac{\partial C}{\partial X}-\frac{\partial \Psi }{\partial X}\frac{\partial C}{\partial Y}=\frac{1}{ReLe}\left(\frac{{\partial }^{2}C}{\partial X^2}+\frac{{\partial }^{2}C}{\partial Y^2}\right)+\frac{1}{ReLe}\left(\frac{Nt}{Nb}\right)\left(\frac{{\partial }^2\theta }{\partial X^2}+\frac{{\partial }^2\theta }{\partial Y^2}\right).$$

(19)

$$\frac{\partial \Psi }{\partial Y}\frac{\partial N}{\partial X}-\frac{\partial \Psi }{\partial X}\frac{\partial N}{\partial Y}+\left(\frac{Pe}{ReLb}\right)\left[\begin{array}{l}\frac{\partial }{\partial X}\left(\left(N+\Omega \right)\frac{\partial C}{\partial X}\right)\\ +\frac{\partial }{\partial Y}\left(\left(N+\Omega \right)\frac{\partial C}{\partial Y}\right)\end{array}\right]=\frac{1}{ReLb}\left(\frac{{\partial }^{2}N}{\partial X^2}+\frac{{\partial }^{2}N}{\partial Y^2}+\frac{{\partial }^{2}N}{\partial X\partial Y}\right).$$

(20)

$$Nb=\frac{\tau D_B\left(c_w-c_{\infty }\right)}{\nu },Nt=\frac{\tau D_T\left(T_w-T_{\infty }\right)}{\nu T_{\infty }},Gr=\frac{g{\beta }_f\left(1-c_{\infty }\right)\left(T_w-T_{\infty }\right)L^3}{{\nu }^2}.$$

(21)

$$Nr=\frac{\left({\rho }_p-{\rho }_f\right)\left(c_w-c_{\infty }\right)}{{\rho }_f{\beta }_f\left(1-c_{\infty }\right)\left(T_w-T_{\infty }\right)},Re=\frac{u_{L}L}{\nu },Pr=\frac{\nu }{\alpha },Le=\frac{\nu }{D_B},Lb=\frac{\nu }{D_m}.$$

(22)

$$Pe=\frac{b{W}_c}{D_m}, \Omega =\frac{n_{\infty }}{n_w-n_{\infty }},Rb=\frac{\gamma \left({\rho }_m-{\rho }_f\right)\left(n_w-n_{\infty }\right)}{{\rho }_f{\beta }_f\left(1-c_{\infty }\right)\left(T_w-T_{\infty }\right)}.$$

(23)

The boundary conditions can also be rewritten in terms of Ψ and ω (Equations (24)–(27)).

$$\Psi =0,\omega =-\frac{{\partial }^2\Psi }{\partial X^2},\theta =0,C=0,N=0atX=0.$$

(24)

$$\frac{\partial \Psi }{\partial Y}=X,\omega =-\frac{{\partial }^2\Psi }{\partial X^2},\theta =1,C=1,N=1atY=0.$$

(25)

$$\frac{\partial \Psi }{\partial Y}=0,\omega =-\frac{{\partial }^2\Psi }{\partial X^2},\theta =0,C=0,N=0atY\to \infty .$$

(26)

$$\frac{{\partial }^2\Psi }{\partial X^2}=0,\frac{\partial \omega }{\partial X}=0,\frac{\partial \theta }{\partial X}=0,\frac{\partial C}{\partial X}=0,\frac{\partial N}{\partial X}=0atX=1.$$

(27)

The partial differential equations can be changed to the similarity coordinates by defining the necessary parameters (Equation (28), in which η is the similarity variable and f is the dimensionless stream function) [4]. The parameters defined in Equation (28) satisfy the continuity equations. Using the similarity parameters, the governing equations may be rewritten (Equations (29)–(32)) as may the boundary conditions (Equation (33)). The prime (′) symbol refers to differentiation with respect to the position.

$$\eta =YR{e}^{\frac{1}{2}}, f^{\prime }\left(\eta \right)=\frac{U\left(X,Y\right)}{X}, \theta (\eta )=\theta \left(X,Y\right), C(\eta )=C\left(X,Y\right),N(\eta )=N\left(X,Y\right).$$

(28)

$$f^{‴}+f{f}^{\prime }-{f^{\prime }}^2=0.$$

(29)

$${\theta }^{″}+Pr\left(f+Nb{C}^{\prime }\right)+PrNt{{\theta }^{\prime }}^2=0.$$

(30)

$$C^{″}+Lef{C}^{\prime }+\left(\frac{Nt}{Nb}\right){\theta }^{″}=0.$$

(31)

$$N^{″}+\left(Lbf-Pe{C}^{\prime }\right)N^{\prime }-Pe{C}^{″}N-Pe\Omega C^{″}=0.$$

(32)

$$\begin{array}{l}{f}^{\prime }\left(0\right)=1, \theta \left(0\right)=1, C\left(0\right)=1,{\rm N}\left(0\right)=1\\ f^{\prime }\left(\infty \right)=0, \theta \left(\infty \right)=0, C\left(\infty \right)=0,{\rm N}\left(\infty \right)=0\end{array}$$

(33)

The wall shear stress, τ_∞, the heat flux, q_w, the mass flux of the nano-additives, q_m, and the flux of the microorganisms, q_n, on the sheet surface are defined in Equation (34). The essential physical parameters in this investigation are the skin friction coefficient, C_f, the Nusselt number, Nu, the Sherwood number, Sh, and the concentration of the microorganisms, Nn (Equation (35)).

$${\tau }_w=\mu {\left.\frac{\partial u}{\partial y}\right|}_{y=0}, q_w=-k{\left.\frac{\partial T}{\partial y}\right|}_{y=0}, q_m=-D_B{\left.\frac{\partial c}{\partial y}\right|}_{y=0}, q_n=-D_n{\left.\frac{\partial n}{\partial y}\right|}_{y=0}$$

(34)

$$C_f=\frac{{\tau }_w}{\frac{1}{2}{\rho }_f{u}_w^2},Nu=\frac{x{q}_w}{k_f\left(T_w-T_{\infty }\right)},Sh=\frac{x{q}_m}{D_B\left(c_w-c_{\infty }\right)},Nn=\frac{x{q}_n}{D_n\left(n_w-n_{\infty }\right)}.$$

(35)

Combining the dimensionless variables (Equation (13)) and the physical parameters (Equation (34)), new relations can be obtained (Equation (36)). Additional reduced physical quantities for the skin friction coefficient, C_fr, the Nusselt number, Nu_r, the Sherwood number, Sh_r, and the density number of the microorganisms, Nn_r, may also be defined (Equation (37)).

$$C_f=\frac{2}{X^{2}Re}{\left.\frac{\partial U}{\partial Y}\right|}_{Y=0},Nu=-X{\left.\frac{\partial T}{\partial Y}\right|}_{Y=0},Sh=-X{\left.\frac{\partial C}{\partial Y}\right|}_{Y=0},Nn=-X{\left.\frac{\partial N}{\partial Y}\right|}_{Y=0}.$$

(36)

$$\begin{array}{l}{C}_{fr}=\frac{2}{XR{e}^{\frac{1}{2}}}{\left.\frac{\partial U}{\partial Y}\right|}_{Y=0},N{u}_r=-R{e}^{-\frac{1}{2}}{\left.\frac{\partial T}{\partial Y}\right|}_{Y=0},S{h}_r=-R{e}^{-\frac{1}{2}}{\left.\frac{\partial C}{\partial Y}\right|}_{Y=0,}\\ N{n}_r=-R{e}^{-\frac{1}{2}}{\left.\frac{\partial N}{\partial Y}\right|}_{Y=0}\end{array}$$

(37)

2.2. Numerical Approach, Mesh Sensitivity Analysis, and Validity

The computational domain needs to be limited to enable the solution to the CFD model to be derived. Therefore, the upper edge should be located at an appropriate distance from the sheet, far away enough from the sheet that the boundary layer fits entirely within the domain but not so far away as to prevent the computer from reaching a converged solution. The height of the computational domain can be approximated by using Equation (38) [4].

$$\frac{H}{L}={\eta }_{\infty }R{e}^{-\frac{1}{2}}.$$

(38)

The finite difference method is employed to discretize the partial and ordinary differential equations (PDEs and ODEs). The advection and diffusion terms in the PDEs are discretized by utilizing the upwind and central schemes, respectively. The point successive over-relaxation method is applied to solve the discretized equation of the stream function whereas the remaining PDEs employ the alternating direction implicit algorithm. The point successive over-relaxation method is used for solving a linear system of the equations, leading to a faster convergence. The tridiagonal matrix algorithm is employed to solve the controlling ODEs. FORTRAN is used for the simulations.

Water (Pr = 6.2) is used as the base fluid because motile microorganisms only survive in water [40]. The nanofluid boundary layers over stretching sheets have previously been studied with Nb and Nt in the range of 0.1–0.5 and Le in the range of 1–25 [4,8], Nr between 0.1 and 0.5 [39], and Re between 20 and 1000 [4]. Considering the parameters related to bioconvection and the microorganisms, Pe and Lb have both been considered in the range of 0–3 [44], Rb in the range 0–2.4 [43], and Ω has been studied in the range of 0–1 [39].

As there are drastic gradients of velocity, temperature, nanoparticle concentrations, and density of microorganisms near the solid boundaries, the mesh in these regions should be denser than in the other regions. A grid independence test is carried out for Re = 10³, Gr = 10⁴, Le = 10, Nb = Nt = Nr = Rb = Ω = 0.1, and Pe = Lb = 1. These values are selected to be close to the validation case (Figure 2). The solution process is run until the convergence criteria reach less than 10⁻⁶. The error percentage of the reduced values is less than 1% when the grid size is 135 × 135 (

Table 1. Grid independence test with errors calculated with respect to the finest grid tested.

Grid Size	Nnr	Error %	Shr	Error %	Nur	Error %	Cfr	Error %
60 × 60	2.3153	13.43	1.7885	15.28	1.6153	71.20	−1.781	11.72
90 × 90	2.6747	0.40	2.1110	0.56	0.9435	3.37	−2.018	1.67
135 × 135	2.6639	0.02	2.0992	0.11	0.9764	0.78	−1.985	0.52
200 × 200	2.6644	-	2.1015	-	0.9841	-	−1.995	-

). Accordingly, a mesh of size 135 × 135 is chosen for subsequent computations.

For validation, the results concluded from this study were evaluated with those reported by Noghrehabadi et al. [4]. The excellent conformity between the present and published results proves the accuracy and reliability of the current code (Figure 2).

3. Results and Discussion

Figure 3 demonstrates the velocity distribution of the different distances from the beginning of the sheet (extrusion slit) for the various values of Reynolds number, Re. When the Re is low (Re = 25), the velocity exhibits a maximum close to the extrusion slit.

This velocity peak is attributed to the fact that, because of temperature gradient, the buoyancy force at the beginning of the sheet is high and decreases when moving to the other end of the sheet. On the other hand, the shear force (applied on the fluid) imposed by the stretching sheet is a function of the sheet velocity and increases along the sheet. As a result, when the velocity of the stretching sheet is low, the buoyancy force (natural convection) is dominant at the beginning of the sheet and a velocity peak can be observed near the sheet. However, for higher shear forces, the velocity field is dominated by the velocity at the surface of the sheet and, therefore, the maximum velocity can be observed there.

At X ≥ 0.5 (i.e., more than half way along the sheet), the velocity distribution decays exponentially from a maximum value at the sheet surface to zero far from the sheet. This trend is consistent with the physics of fluid flow. Furthermore, the length of the sheet at which a boundary layer is fully developed increases with the Reynolds number. When Re = 200, there is no velocity overshoot. When Re > 200, the velocity profiles at all X values match the similarity analysis. Hence, Re = 200 can be considered to be the critical Reynolds number to develop the desired flow.

The influence of the Reynolds number, Re, on the dimensionless temperature distribution shows that the dimensionless temperature decreases exponentially from a maximum value on the sheet surface to a value of zero at a far distance from the sheet for all values of the Reynolds number and at all locations (Figure 4). The growth rate of the thermal boundary layer increases with the Reynolds number, Re. As a result, the thermal boundary layer can be fully developed in lower X values.

The length of the thermal region that is not fully developed increases so that at Re = 200 and for X > 0.1, the temperature profiles of the CFD and the similarity solutions match. The thermal boundary layer approximated by the similarity analysis also does not match any Reynolds number for X < 0.05.

In a similar way to the temperature boundary layer, at the lowest value of Reynolds number, Re = 25, no boundary layer for the concentration of nanoparticles ever fully develops (Figure 5). The similarity analysis does not match the CFD results for X < 0.05. The entrance length (the length of the undeveloped region) for the concentration of the microorganisms is longer than that for the other boundary layers (Figure 6). As the Reynolds number increases, the developed regions for the concentration of the nanoparticles and the density of the microorganism form and grow. As a result, the concentration boundary layers of the nanoparticles and microorganisms can be fully developed in lower X values.

As shown in Figure 7, the changes of the reduced friction coefficient, C_fr, along the x-axis for the different Reynolds numbers show that at low Reynolds numbers (25 ≤ Re ≤ 100), C_fr is very high near the extrusion slit due to the drastic velocity gradients. The velocity gradient reduces with an increasing Reynolds number and so does C_fr.

To distinguish the developed and undeveloped regions, the relative errors between the outcomes of the CFD (obtained by solving the PDEs) and the similarity analysis (obtained by solving ODEs) may be calculated (Equation (39) where χ represents a reduced skin friction coefficient, C_fr, a reduced Nusselt number, Nu_r, a reduced Sherwood number, Sh_r, or a reduced concentration of microorganisms, Nn_r).

$$\Delta =\left|\frac{{\chi }_{CFD}-{\chi }_{sim}}{{\chi }_{sim}}\right|\times 100.$$

(39)

Critical Reynolds and Grashof numbers can be obtained for the different locations along the sheet (different X values) based on Equation (39) by selecting a threshold below which the CFD and the similarity analysis are accepted to agree. For example, Figure 8 shows that the critical Reynolds numbers for the boundary layers for the concentrations of nanoparticles and microorganisms are equal as a function of X when the maximum permitted discrepancy between the results of the CFD calculations and the similarity analysis is Δ < 5%. As seen, the critical Reynolds numbers reduce as the distance from the slit increases; i.e., for a larger X. When X < 0.065, the boundary layers approximated by the similarity analysis are also not valid for any Reynolds number below 1000. For X ≥ 0.1, the critical Reynolds numbers are 200, 180, 400, and 400 for the reduced skin friction coefficient, C_fr, the reduced Nusselt number, Nu_r, the reduced Sherwood number, Sh_r, and the reduced concentration of nanoparticles, Nn_r, respectively.

The average error versus the Reynolds number is shown in Figure 9. When Re < 80, the average error in the reduced skin friction coefficient, C_fr, is 34% whereas the average error of the other characteristics for all Reynolds numbers is less than 21%. The results also show that for a maximum permitted discrepancy between the results of the CFD calculations and the similarity analysis, Δ ≤ 5%; the approximations of the similarity analysis are acceptable when Re ≥ 180, 200, 400, and 400 for the reduced skin friction coefficient, C_fr, the reduced Nusselt number, Nu_r, the reduced Sherwood number, Sh_r, and the reduced concentration of microorganisms, Nn_r, respectively.

A high overshoot in the velocity appears when the Grashof number Gr = 10⁶ (Figure 10). This can be explained by the fact that the flow field near the slit is governed by the buoyancy force when the Grashof number is high. Hence, the maximum velocity cannot be observed on the sheet [38]. Increasing the Grashof number leads to a slower development of the flow, which means that more of the sheet is in the region where the boundary layer is not fully developed. The Grashof number can be thought of as resisting the development of the flow and boundary layer.

As the Grashof number, Gr, increases, so does the reduced skin friction coefficient, C_fr, along the entire sheet surface. This is because the buoyancy is more significant as Grashof number increases. This causes the velocity gradient to increase and so C_fr increases in the region before the flow is fully developed. The length of the region in which the boundary layer is not fully developed also increases when the Grashof number increases (Figure 11). In a similarity analysis, the boundary layer is a constant thickness everywhere across the sheet and so the CFD solution diverges from the similarity analysis as the Grashof number increases.

Figure 12a shows the critical Grashof number corresponding with the distance from the extrusion slit for Re = 400, Pr = 6.2, Le = 10, Nb = Nt = Nr = Rb = Ω = 0.1, Pe = Lb = 1, and Δ < 0.05. As seen, for X ≤ 0.1 and Δ < 0.05, the velocity boundary layer approximated by the similarity analysis is valid when Gr ≤ 8 × 10⁴. The average error in the reduced skin friction coefficient, $\overline{\Delta C_{\mathrm{fr}}}$, as a function of the Grashof number varies linearly with Gr (Figure 12b). The critical Grashof number for approximating the hydrodynamic boundary layer is $4\times {10}^4$ when $\overline{\Delta C_{\mathrm{fr}}}$ ≤ 5%.

The effect of the Grashof number on the reduced density number of microorganisms, Nn_r, along the sheet surface for Re = 400, Pr = 6.2, Le = 10, Nb = Nt = Nr = Rb = Ω = 0.1, and Pe = Lb = 1 is small when the Grashof number increases from 10³ to 10⁵; however, as the Grashof number increases further to 10⁶, there is a strong effect (Figure 13). The gradient of the concentration of microorganisms reduces in the undeveloped region as the Grashof number increases, which results in a thicker boundary layer. It can be observed that, due to the coupling between the flow field and the concentration field of the microorganisms, the graph for the reduced density of microorganisms (Nn_r) is completely consistent with the velocity field (f’(η)) observed in high Grashof numbers (Figure 10). Increasing the velocity at the beginning of the sheet (equivalent to increasing Gr) boosts Nn_r. However, increasing the Grashof number too much results in the fluid velocity becoming negative (viz. the fluid flows in the negative x-direction), as shown in Figure 10 at distances far from the slit. This leads to a reduction in Nn_r.

The effects of the bioconvection Rayleigh number, Rb, on the reduced skin friction coefficient, C_fr, the reduced Nusselt number, Nu_r, the reduced Sherwood number, Sh_r, and the reduced concentration of microorganisms, Nn_r, are depicted in Figure 14. Increasing Rb causes C_fr to become smaller because, in the buoyancy term of the momentum equation, Rb creates a negative buoyancy, causing a resistant flow against the forced convection. When Rb becomes more significant, this resistance to the flow is stronger. Consequently, the velocity diminishes within the hydrodynamic boundary layer, leading to a decrease in C_fr. An increase in Rb can also cause a slight decrease in Nu_r, Sh_r, and Nn_r in the undeveloped region. The reduction in the net flow velocity causes the rate of heat, mass, and transfer of the microorganisms to diminish.

Increasing the difference in the density of the microorganisms from the surface to the far-field, Ω, causes a slight increase in the reduced Sherwood number, Sh_r, and a larger increase in the reduced Nusselt number, Nn_r, along the sheet surface (Figure 15). The density of the swimming microorganisms diminishes as Ω increases, so the mobility of the nanoliquid increases, as does the reduced Sherwood number, Sh_r (Figure 15a). Increasing Ω leads to an increase in the gradient in the density of the microorganisms between the sheet surface and the far-field. Therefore, the rate of transfer of the microorganisms, Nu_r, increases as Ω increases (Figure 15b).

4. Conclusions

This work studied the boundary layer flow of a nanoliquid, which included swimming microorganisms, in the regions of developed and undeveloped flow across a linearly stretching sheet. To distinguish between the developed and the underdeveloped regions, the computational fluid dynamics results were compared with the results of a similarity analysis. The findings are summarized below:

• At low values of Reynolds numbers (Re = 25 and Re = 100), the velocity experienced an overshoot in the vicinity of the extrusion slit. However, the maximum velocity occurred at the sheet surface when the Reynolds number was high.
• The similarity analysis could not approximate the flow characteristics near the extrusion slit at low values of Reynolds numbers.
• Having a mismatch less than 5% between the CFD and the similarity analysis at a position of 6.5% of the way along the sheet, the critical Reynolds number above which the similarity analysis matched the computational fluid dynamics solution was more than 1000 for Gr = 10⁴, Pr = 6.2, Le = 10, Nb = Nt = Nr = Rb = Ω = 0.1, and Pe = Lb = 1.
• An increase in the Grashof number, Gr, increased the length of the region in which the flow was undeveloped, leading to a decrease in the reliability of the similarity analysis. It was found that the similarity analysis was within 5% of the computational fluid dynamics solution at a position of 10% of the way along the sheet and beyond this for Gr < 8 × 10⁴, Re = 400, Pr = 6.2, Le = 10, Nb = Nt = Nr = Rb = Ω = 0.1, and Pe = Lb = 1.
• Increasing the Grashof number from 10³ to 10⁵ did not affect the length of the region in which there was no fully developed boundary layer for the density of the microorganisms. However, this region diminished strongly as the Grashof number increased from 10⁵ to 10⁶.
• Increasing the bioconvection Rayleigh number, Rb, caused the reduced skin friction coefficient, C_fr, to decrease, but had a weaker effect on the reduced Nusselt number, Nu_r, the reduced Sherwood number, Sh_r, and the reduced concentration of microorganisms, Nn_r.
• Increasing the difference in the density of the microorganisms from the surface to the far-field, $\Omega$, led to an increase in the reduced Nusselt number, Nu_r.

This study provides evidence that caution should be taken when a similarity solution is employed for solving the flow characteristics over a stretching sheet.