# Micropolar Couple Stress Nanofluid Flow by Non-Fourier’s-Law Heat Flux Model past a Stretching Sheet

Authors: Gosa Gadisa; Tagay Takele; Shibiru Jabessa
Journal: Journal of Mathematics
 (2021)

---

## Abstract

In this investigation, thermal radiation effect on MHD nonlinear convective micropolar couple stress nanofluid flow by non-Fourier’s-law heat flux model past a stretching sheet with the effects of diffusion-thermo, thermal-diffusion, and first-order chemical reaction rate is reported. The robust numerical method called the Galerkin finite element method is applied to solve the proposed fluid model. We applied grid-invariance test to approve the convergence of the series solution. The effect of the various pertinent variables on velocity, angular velocity, temperature, concentration, local skin friction, local wall couple stress, local Nusselt number, and local Sherwood number is analyzed in both graphical and tabular forms. The range of the major relevant parameters used for analysis of the present study was adopted from different existing literatures by taking into consideration the history of the parameters and is given by0.07≤Pr≤7.0,0.0≤λ,ε≤1.0,0.0≤Rd,Df,Sr,K,≤1.5,0.0≤γE≤0.9,0.9≤Sc≤1.5,0.5≤M≤1.5,0.0≤β≤1.0,0.2≤Nb≤0.4,0.1≤Nt≤0.3. The result obtained in this study is compared with that in the available literatures to confirm the validity of the present numerical method. Our result shows that the heat and mass transfer in the flow region of micropolar couple stress fluid can be enhanced by boosting the radiation parameters.


---

## Body

## 1. Introduction

Numerous mathematical models were proposed to study the rheological properties of non-Newtonian fluids. The fluid model pioneered by Eringen [1] in 1996 revealed the existence of microscopic effects resulting from the local structure and micromotion of the fluid constituents. Moreover, they can sustain couple stresses and comprise the Newtonian models as a special case. The importance of heat and mass transfer and micropolar ﬂuid ﬂow is particularly evident in new and emerging areas of materials processing. Materials such as polymers, alloys, ceramics, composites, semiconductors, and optical materials need thermal energy for fabrication. For instance, temperature control helps to ensure the product quality and consistent production capacity in the polymer extrusion process. It is due to this fact that the heat and mass transfer in the boundary layer flow of non-Newtonian fluid with diverse effects of parameters attracted scholars from all corners [2–6]. These relevant parameters are used to control the heat and mass transfer during the extrusion process. In some studies, non-Fourier’s-law heat flux model is applied to govern the heat and mass transfer in the boundary layer flow region [7, 8].Thermally radiative fluid flows are usually encountered when the difference between the temperature at the surface of the sheet and the ambient temperature is high. In numerous industrial processes, the thermal boundary layer thickness can be altered by the use of the thermal radiation. Examples of such industrial processes include missile technology, nuclear reactors, satellites, power plants, and gas turbines. The effect of radiation parameters in different boundary layer flow regions has been introduced by numerous scholars [9–14]. Mixed convection can be considered as a combination of free and forced convection which occur due to a significant difference in temperature between the surface (wall) and the ambient fluid. Mixed convection has an essential role when the buoyancy force considerably disrupts the flow and thermal fields. Ramzan et al. [15] studied a mixed convection viscoelastic nanofluids past porous media considering Soret-Dufour effect. They employed a homotopy analysis method to solve the proposed problem. Some authors [16–18] analyzed the influence of mixed convection parameters on the boundary layer flow of Oldroyd-B fluid. Most recently, Ibrahim and Gadisa [19] reported the nonlinear convective flow of a couple stress-micropolar nanofluid with the effects of slip and convective boundary conditions.On the other hand, in chemical process engineering, Dufour and Soret have a vital application. Moorthy and Senthilvadivu [20] suggested that when heat and mass transfer processes take place at the same time between the fluxes, the driving potential is of more complex nature, as energy flux can be generated not only by temperature gradients but also by composition gradients. According to the assumption of Fick’s law or Fourier’s, Soret and Dufour effects are typically dilapidated in heat and mass transfer processes. There are, however, exceptions in certain circumstances. For instance, the Soret effect can be utilized for isotope separation and in mixtures between gases with extremely light molecular weight like H2 and He. For average molecular weights like N2 and air, the Dufour effect was found to be of a significant magnitude such that it cannot be neglected [21]. Ahammad and Mollah [22] introduced the concepts of MHD free convection flow and mass transfer over a stretching sheet with Dufour and Soret effects. They solved numerically by applying Runge-Kuta with the shooting technique. Soret and Dufour effect on MHD Casson fluid past a stretching sheet was studied by Hayat et al. [23]. Later on, Ali and Shah [24] reported free-convection MHD micropolar fluid considering Soret and Dufour effects. Following this, different scholars [25–29] scrutinized the effects of Soret and Dufour on micropolar fluid. Most recently, Bhubaneswar et al. [30] forwarded the concept of cross diffusion effects on MHD convection of Casson-Williamson fluid past a stretching surface with radiation and chemical reaction. Ibrahim and Gadisa [19] also analyzed the nonlinear convective flow of a couple stress-micropolar nanofluids with non-Fourier heat flux model past the stretching surface in the presence of slip and convective boundary conditions.The present scrutiny is motivated by the results of the last two papers. As far as we revised and mentioned above, the problem of nonlinear convective flow of micropolar couple stress nanofluid using the Cattaneo-Christov model past the stretching sheet with the effects of thermal radiation and Soret, Dufour, and chemical reaction is still unnoticed. Thus, the main objective of the present study is to fill this gap. We employed the robust numerical technique called GFEM explained in equations (21)–(31). We performed the grid-independence test or grid convergence test to confirm the convergence of the series solution. The effect of the relevant parameters on linear velocity, angular velocity (microrotation), temperature, concentration, local skin friction, local wall couple stress, local Nusselt number, and local Sherwood number is elaborated in both graphical and tabular forms.


## 2. Problem Formulation

In this study, we consider the two-dimensional steady incompressible laminar MHD boundary layer flow of a nonlinear convective micropolar couple stress nanofluid using the Cattaneo-Christov heat flux model in the presence of thermal radiation with the effects of Soret, Dufour, and chemical reactions past the stretching sheet as plotted in Figure1 below. The sheet is stretched linearly with velocity uwx=ax, where a is constant. Applying these suppositions, the governing boundary layer equations with Boussinesq approximations are as follows (Wubshet and Gosa [19]):(1)∂u∂x+∂v∂y=0,(2)u∂u∂x+v∂v∂y=μ+kρf∂2u∂y2+kρf∂N∂y−v′∂4u∂y4−σB02xρfu+g∧1T−T∞+g∧2T−T∞2+g∧3C−C∞+g∧4C−C∞2,(3)u∂N∂x+v∂N∂y=Ωρj∂2N∂y2−kρj2N+∂u∂y,(4)u∂T∂x+v∂T∂y+λEu∂u∂x∂T∂x+v∂v∂y∂T∂y+u∂v∂x∂T∂y+v∂u∂y∂T∂x+2uv∂2T∂x∂y+u2∂2T∂x2+v2∂2T∂y2=kfρfcp∂2T∂y2−1ρfcp∂qr∂y+σDm∂T∂y∂C∂y+DTT∞∂T∂y2+DmkTcscp∂2C∂y2,(5)u∂C∂x+v∂C∂y=Dm∂2C∂y2+DTT∞∂2T∂y2+DmkTTm∂2T∂y2−k1C−C∞,with the following boundary conditions at the surface (i.e., y=0). u=uwx=ax,v=0,N=−0.5∂u/∂y,T=TW,C=CW, and boundary conditions at the far field (i.e., y⟶C∞):(6)u⟶0,N⟶0,T⟶T∞,andC⟶C∞,where uandv are the velocity components in the x and y directions, ν=μ/ρf is the kinematic viscosity, ν′=n/ρf is the couple stress viscosity, n is the couple stress viscosity parameter, μ is the dynamic viscosity, k is the vortex viscosity, kf is the thermal conductivity of the fluid, j is the microinertia density, Dm is the mass diffusivity, and Ω is the spin-gradient viscosity and given by Ω=μ+k/2j=μ1+β/2j, where β=k/μ, uw is the stretching velocity, g is the gravitational acceleration, λ1 and λ2 are the relaxation time and the retardation time, respectively, λEandλC are the Deborah numbers with respect to relaxation time of heat flux and nanoparticles concentration, respectively, Λ1andΛ2 are the linear and nonlinear thermal expansion coefficients due to temperature, Λ3andΛ4 are the linear and nonlinear thermal expansion coefficients due to concentration, ρf is the density of base liquid, cs is the concentration susceptibility, cp is the specific heat capacity of the base fluid, T is the temperature, C is the concentration, σ is the electric conductivity, B02 is the magnetic parameter, k1 is the reaction rate, kT is the thermal-diffusion ratio, Tm is the mean fluid temperature, qr=−4/3σ∗/k∗∂T4/∂y is the radiative heat flux, where σ∗ is the Stefan-Boltzmann constant, and k∗ is the mean absorption coefficient.Figure 1
Geometry of the flow problem.The coupled nonlinear partial differential equations (PDEs)1−5 with boundary conditions (6) above can be reduced to the appropriate coupled nonlinear ordinary differential equations (ODEs) by using the following similarity and dimensionless variables.η=ya/v is the similarity variable and the dimensionless variables f,p,θ,andϕ are defined as follows:(7)ψ=avxfη,N=axavpη,θη=T−T∞Tw−T∞,ϕη=C−C∞CW−C∞.The continuity equation can be satisfied if we define the stream functions as follows:(8)u=∂ψ∂y,v=−∂ψ∂x.Using the introduced similarity transformation above (7), we reduced equations (1)–(5) with the boundary condition (6) to the following coupled highly nonlinear ODEs:(9)1+βf‴+ff″−f′2+βp′−Kfv−Mf′+λθ1+βtθ+λN∗ϕ1+βcϕ=0,1+β2p″−β2p+f″+fp′−f′p=0,1Pr1+43Rdθ″+fθ′−γEff′θ′+f2θ″+Nbθ′ϕ′+DfNtθ′2ϕ″=0,1Scϕ″+fϕ′+Srθ″+1ScNtNbθ″−εϕ=0,with the corresponding boundary conditions:(10)f0=0,f′0=1,p0=−0.5f″0,θ0=1,ϕ0=1atη=0and,f′η⟶0,pη⟶0,θη⟶0,ϕη⟶0asη⟶∞.Here, we affirm the dimensionless parameters as follows:(11)Pr=ρfνfcpkf,γE=λEa,M=σB0aρf,K=ν′aν,β=kμ,λ=GrRex2,βt=Λ2Λ1Tw−T∞,βc=Λ4Λ3Cw−C∞,N∗=Gr∗Gr=Λ3Λ1Cw−C∞Tw−T∞,Gr∗=gΔ3ν2Cw−C∞x3,Sc=vDm,ε=k1a,Gr=gΛ1ν2Tw−T∞x3,Rd=4σ∗T∞3kfk∗,Df=DmKTvcfcpCw−C∞Tw−T∞,Sr=DmKTvTmTw−T∞Cw−C∞,Nb=DBσCw−C∞ν,Nt=DTσTw−T∞T∞ν,where Pr is the Prandtl number, γE is the Deborah number with respect to the relaxation time of the heat flux, M is the magnetic field parameter, λ is the mixed convection parameter/local buoyancy parameter, βt is the nonlinear convection parameter due to temperature, βc is the nonlinear convection parameter due to concentration, N∗ is the ratio of concentration to thermal buoyancy forces, Gr is the Grashof number in terms of temperature, Gr∗ is the Grashof number in terms of concentration, Sc is the Schmidt number, Sr is the Soret, Df is the Dufour number, ε is the chemical reaction rate, β is the material parameter, K is the couple stress parameter, Nb is the Brownian motion parameter, and Nt is the thermophoresis parameter.The engineering physical quantities of interest in this paper are skin friction coefficient, wall couple stress, Nusselt number, and Sherwood number defined as follows:(12)Cf=τwρuw2,Mw=mwρuw2,Nux=xqwkTw−T∞,Shx=xqnpDBCw−C∞,(13)τw=μ+κ∂u∂y+κΝ−n∂3u∂y3y=0,qw=−k∂T∂yy=0,qnp=−DB∂ϕ∂yy=0,mw=Ω∂N∂yy=0.where the wall shear stress τw, wall heat flux qw, and wall mass flux qnp are defined as follows:After substituting the values introduced in the similarity transformation above, we get(14)CfRex=1+12βf″0−Kfiv0,Mw=1+β2jag′0,NuxRex=−1+4Rd3θ′0,ShxRex=−ϕ′0.


## 3. Numerical Simulations

### 3.1. Finite Element Method

We applied the robust numerical technique called the Galerkin finite element method to obtain the solution of coupled nonlinear partial differential equations governing the boundary layer flow. The fundamental steps needed to apply GFEM are dividing the domain into elements/discretization of the domain, the element formulation/derivation of the element equation, the assemblage of the element equation into its global form, and the imposition of boundary conditions and solving a system of linear equations, respectively, (see [19]). The following steps are crucial to apply the finite element method.


#### 3.1.1. Discretization of the Domain

The fundamental concept of the FEM is to divide the domain or region of the problem into small connected parts, called finite elements. The collection of elements is called the finite element mesh. These finite elements are connected in a nonoverlapping manner, such that they completely cover the entire space of the problem.


#### 3.1.2. Generation of the Element Equations

(i)
A typical element is isolated from the mesh and the variational formulation of the given problem is constructed over the typical element.(ii)
Over an element, an approximate solution of the variational problem is supposed, and by substituting this in the system, the element equations are generated.(iii)
The element matrix, which is also known as the stiﬀness matrix, is constructed by using the element interpolation functions.


#### 3.1.3. Assembly of the Element Equations

The algebraic equations obtained are assembled by imposing the interelement continuity conditions. This yields a large number of algebraic equations known as the global finite element model, which govern the whole domain.


#### 3.1.4. Imposition of Boundary Conditions

On the assembled equations, both Dirichlet and Neumann boundary conditions are imposed.


#### 3.1.5. Solution of the Assembled Equations

The assembled equations so obtained can be solved by any of the numerical techniques, namely, LU decomposition method, Gauss elimination method, and so forth.Assuming,(15)f′=g.The system of differential equations above will be reduced to the following equations:(16)1+βg″+fg′−g2+βp′−Kgiv−Mg+λθ1+βtθ+λN∗ϕ1+βcϕ=0,(17)1+β2p″−β2p+g′+fp′−gp=0,(18)1Pr1+43Rdθ″+fθ′−γEfgθ′+f2θ″+Nbθ′ϕ′+DfNtθ′2ϕ″=0,(19)1Scϕ″+fϕ′+Srθ″+1ScNtNbθ″−εϕ=0,with the following corresponding boundary conditions:(20)f0=0,g0=1,p0=−0.5g′0,θ0=1,ϕ0=1atη=0,and,gη⟶0,pη⟶0,θη⟶0,ϕη⟶0asη⟶∞.


### 3.2. Variational Formulation

The variational formulations used in solving differential equations by the finite element method are considered in detail by different scholars [19, 31–34].The variational formulation related to the equations (15)–(19) over a typical element η,ηe+1 is given by(21)∫ηeηe−1w1f′−gdη=0,(22)∫ηeηe+1w21+βg″+fg′−g2+βp′−Kfiv−Mg+λθ1+βtθ+λN∗ϕ1+βcϕdη=0,(23)∫ηeηe+1w31+β2p″−β2p+g′+fp′−gpdη=0,(24)∫ηeη+1w41Pr1+43Rdθ″+fθ′−γEfgθ′+f2θ″+Nbθ′ϕ′+DfNtθ′2ϕ″dη=0,(25)∫ηeηe+1w51Scϕ″+fϕ′+Srθ″+1ScNtNbθ″−εϕdη=0,subjected to the boundary condition (20), where w1w2w3, w4, and w5 are arbitrary weight functions may be regarded as variations in f,g,p,θ,andϕ, respectively.


### 3.3. Finite Element Formulation

At the third step, we look for the approximation solution of the following form:(26)f=∑j=13fjψj,g=∑j=13gjψj,p=∑j=13pjψj,θ=∑j=13θjψj,ϕ=∑j=13ϕjψj,with w1=w2=w3=w4=w5=wii=1,2,3, the quadratic shape functions ψi are defined as(27)ψ1e=ηe+1−ηηe+1+ηe−2ηηe+1−ηe2,ψ2e=4η−ηeηe+1−ηηe+1−ηe2,ψ3e=−η−ηeηe+1+ηe−2ηηe+1−ηe2.where ηe≤η≤ηe+1.At the fourth step, substituting the approximate solution of equation (26) into equations (21)–(25), we acquire the finite element model equation which is given by(28)KeYe=Fe,where Ke denotes the elemental stiffness matrix, Ye is the vector of elemental nodal variables (unknowns), and Fe is the force vector expressed as follow:(29)Ke=K11K12K13K14K15K21K22K23K24K25K31K32K33K34K35K41K42K43K44K45K51K52K53K54K55,Ye=fgpθϕ,Fe=r1r2r3r4r5,where each Kmn is of the order 3x3and rm,m,n=1,2.3,4,5 is of order 3x1. These matrices are defined as(30)Kij11=∫ηeηe+1ψi∂ψj∂ηdη,Kij12=−∫ηeηe+1ψiψjdη,Kij13=0,Kij14=0,Kij15=0,Kij22=−1+β∫ηeηe+1∂ψi∂η∂ψj∂ηdη+∫ηeηe+1ψif¯∂ψj∂ηdη−∫ηeηe+1ψig¯ψjdη−K∫ηeηe+1∂2ψi∂η2∂2ψj∂η2dη−M∫ηeηe+1ψiψjdη,Kij23=β∫ηeηe+1ψi∂ψj∂η,Kij21=0,Kij24=λ∫ηeηe+1ψiψjdη+λβt∫ηeηe+1ψiθ¯ψjdη,Kij25=λN∗∫ηeηe+1ψiψjdη+λN∗βc∫ηeηe+1ψiϕ¯ψjdη,Kij31=Kij34=Kij35=0,Kij32=−β∫ηeηe+1ψi∂ψj∂ηdη,Kij33=−1+β2∫ηeηe+1∂ψi∂η∂ψj∂ηdη−2β∫ηeηe+1ψiψjdη+∫ηeηe+1ψif¯∂ψj∂ηdη−∫ηeηe+1ψig¯ψjdη,Kij44=−1Pr1+4Rd3∫ηeηe+1∂ψi∂η∂ψj∂ηdη+∫ηeηe+1ψif¯∂ψj∂ηdη−γE∫ηeηe+1ψifg¯∂ψj∂ηdη−γE∫ηeηe+1ψiff¯∂2ψj∂η2dη+Nb∫ηeηe+1ψiϕ′¯∂ψj∂ηdη+Nt∫ηeηe+1ψiθ′¯∂ψj∂ηdη,K45=−Df∫ηeηe+1∂ψi∂η∂ψj∂ηdη,K41=K42=K43=K51=K52=K53=0,Kij54=−Sr+NtNb∫ηeηe+1∂ψi∂η∂ψi∂ηdη,Kij55=−1Sc∫ηeηe+1∂ψi∂η∂ψi∂ηdη+∫ηeηe+1ψif¯∂ψi∂ηdη−ε∫ηeηe+1ψiψjdη,ri1=0,ri2=−1+βψi∂g∂ηηeηe+1,ri3=−1+β2ψi∂p∂ηηeηe+1,ri4=−1Pr1+4Rd3ψi∂θ∂ηηeηe+1−Dfψi∂ϕ∂ηηeηe+1,ri5=−1Scψi∂ϕ∂η+Srψi∂θ∂ηηeηe+1,where(31)f¯′=∑j=13f¯j∂ψj∂η,g¯′=∑j=13g¯j∂ψj∂η,p¯′=∑j=13p¯j∂ψj∂η,θ¯′=∑j=13θ¯j∂ψj∂η,ϕ¯′=∑j=13ϕ¯′∂ψj∂η.


### 3.4. Assembly of the System of Equations

At this step, we assemble the elemental systems to obtain the global system, given as follows:(32)KY=F.This results in a large number of algebraic equations that govern the entire domain. And then, the global matrix will be modified by imposing the boundary conditions defined in equation (20). The last step is solving the assembled system of equations by the standard techniques like the Gauss elimination method, LU decomposition, Gauss Jordan method, or any iterative scheme.


## 3.1. Finite Element Method

We applied the robust numerical technique called the Galerkin finite element method to obtain the solution of coupled nonlinear partial differential equations governing the boundary layer flow. The fundamental steps needed to apply GFEM are dividing the domain into elements/discretization of the domain, the element formulation/derivation of the element equation, the assemblage of the element equation into its global form, and the imposition of boundary conditions and solving a system of linear equations, respectively, (see [19]). The following steps are crucial to apply the finite element method.


### 3.1.1. Discretization of the Domain

The fundamental concept of the FEM is to divide the domain or region of the problem into small connected parts, called finite elements. The collection of elements is called the finite element mesh. These finite elements are connected in a nonoverlapping manner, such that they completely cover the entire space of the problem.


### 3.1.2. Generation of the Element Equations

(i)
A typical element is isolated from the mesh and the variational formulation of the given problem is constructed over the typical element.(ii)
Over an element, an approximate solution of the variational problem is supposed, and by substituting this in the system, the element equations are generated.(iii)
The element matrix, which is also known as the stiﬀness matrix, is constructed by using the element interpolation functions.


### 3.1.3. Assembly of the Element Equations

The algebraic equations obtained are assembled by imposing the interelement continuity conditions. This yields a large number of algebraic equations known as the global finite element model, which govern the whole domain.


### 3.1.4. Imposition of Boundary Conditions

On the assembled equations, both Dirichlet and Neumann boundary conditions are imposed.


### 3.1.5. Solution of the Assembled Equations

The assembled equations so obtained can be solved by any of the numerical techniques, namely, LU decomposition method, Gauss elimination method, and so forth.Assuming,(15)f′=g.The system of differential equations above will be reduced to the following equations:(16)1+βg″+fg′−g2+βp′−Kgiv−Mg+λθ1+βtθ+λN∗ϕ1+βcϕ=0,(17)1+β2p″−β2p+g′+fp′−gp=0,(18)1Pr1+43Rdθ″+fθ′−γEfgθ′+f2θ″+Nbθ′ϕ′+DfNtθ′2ϕ″=0,(19)1Scϕ″+fϕ′+Srθ″+1ScNtNbθ″−εϕ=0,with the following corresponding boundary conditions:(20)f0=0,g0=1,p0=−0.5g′0,θ0=1,ϕ0=1atη=0,and,gη⟶0,pη⟶0,θη⟶0,ϕη⟶0asη⟶∞.


## 3.1.1. Discretization of the Domain

The fundamental concept of the FEM is to divide the domain or region of the problem into small connected parts, called finite elements. The collection of elements is called the finite element mesh. These finite elements are connected in a nonoverlapping manner, such that they completely cover the entire space of the problem.


## 3.1.2. Generation of the Element Equations

(i)
A typical element is isolated from the mesh and the variational formulation of the given problem is constructed over the typical element.(ii)
Over an element, an approximate solution of the variational problem is supposed, and by substituting this in the system, the element equations are generated.(iii)
The element matrix, which is also known as the stiﬀness matrix, is constructed by using the element interpolation functions.


## 3.1.3. Assembly of the Element Equations

The algebraic equations obtained are assembled by imposing the interelement continuity conditions. This yields a large number of algebraic equations known as the global finite element model, which govern the whole domain.


## 3.1.4. Imposition of Boundary Conditions

On the assembled equations, both Dirichlet and Neumann boundary conditions are imposed.


## 3.1.5. Solution of the Assembled Equations

The assembled equations so obtained can be solved by any of the numerical techniques, namely, LU decomposition method, Gauss elimination method, and so forth.Assuming,(15)f′=g.The system of differential equations above will be reduced to the following equations:(16)1+βg″+fg′−g2+βp′−Kgiv−Mg+λθ1+βtθ+λN∗ϕ1+βcϕ=0,(17)1+β2p″−β2p+g′+fp′−gp=0,(18)1Pr1+43Rdθ″+fθ′−γEfgθ′+f2θ″+Nbθ′ϕ′+DfNtθ′2ϕ″=0,(19)1Scϕ″+fϕ′+Srθ″+1ScNtNbθ″−εϕ=0,with the following corresponding boundary conditions:(20)f0=0,g0=1,p0=−0.5g′0,θ0=1,ϕ0=1atη=0,and,gη⟶0,pη⟶0,θη⟶0,ϕη⟶0asη⟶∞.


## 3.2. Variational Formulation

The variational formulations used in solving differential equations by the finite element method are considered in detail by different scholars [19, 31–34].The variational formulation related to the equations (15)–(19) over a typical element η,ηe+1 is given by(21)∫ηeηe−1w1f′−gdη=0,(22)∫ηeηe+1w21+βg″+fg′−g2+βp′−Kfiv−Mg+λθ1+βtθ+λN∗ϕ1+βcϕdη=0,(23)∫ηeηe+1w31+β2p″−β2p+g′+fp′−gpdη=0,(24)∫ηeη+1w41Pr1+43Rdθ″+fθ′−γEfgθ′+f2θ″+Nbθ′ϕ′+DfNtθ′2ϕ″dη=0,(25)∫ηeηe+1w51Scϕ″+fϕ′+Srθ″+1ScNtNbθ″−εϕdη=0,subjected to the boundary condition (20), where w1w2w3, w4, and w5 are arbitrary weight functions may be regarded as variations in f,g,p,θ,andϕ, respectively.


## 3.3. Finite Element Formulation

At the third step, we look for the approximation solution of the following form:(26)f=∑j=13fjψj,g=∑j=13gjψj,p=∑j=13pjψj,θ=∑j=13θjψj,ϕ=∑j=13ϕjψj,with w1=w2=w3=w4=w5=wii=1,2,3, the quadratic shape functions ψi are defined as(27)ψ1e=ηe+1−ηηe+1+ηe−2ηηe+1−ηe2,ψ2e=4η−ηeηe+1−ηηe+1−ηe2,ψ3e=−η−ηeηe+1+ηe−2ηηe+1−ηe2.where ηe≤η≤ηe+1.At the fourth step, substituting the approximate solution of equation (26) into equations (21)–(25), we acquire the finite element model equation which is given by(28)KeYe=Fe,where Ke denotes the elemental stiffness matrix, Ye is the vector of elemental nodal variables (unknowns), and Fe is the force vector expressed as follow:(29)Ke=K11K12K13K14K15K21K22K23K24K25K31K32K33K34K35K41K42K43K44K45K51K52K53K54K55,Ye=fgpθϕ,Fe=r1r2r3r4r5,where each Kmn is of the order 3x3and rm,m,n=1,2.3,4,5 is of order 3x1. These matrices are defined as(30)Kij11=∫ηeηe+1ψi∂ψj∂ηdη,Kij12=−∫ηeηe+1ψiψjdη,Kij13=0,Kij14=0,Kij15=0,Kij22=−1+β∫ηeηe+1∂ψi∂η∂ψj∂ηdη+∫ηeηe+1ψif¯∂ψj∂ηdη−∫ηeηe+1ψig¯ψjdη−K∫ηeηe+1∂2ψi∂η2∂2ψj∂η2dη−M∫ηeηe+1ψiψjdη,Kij23=β∫ηeηe+1ψi∂ψj∂η,Kij21=0,Kij24=λ∫ηeηe+1ψiψjdη+λβt∫ηeηe+1ψiθ¯ψjdη,Kij25=λN∗∫ηeηe+1ψiψjdη+λN∗βc∫ηeηe+1ψiϕ¯ψjdη,Kij31=Kij34=Kij35=0,Kij32=−β∫ηeηe+1ψi∂ψj∂ηdη,Kij33=−1+β2∫ηeηe+1∂ψi∂η∂ψj∂ηdη−2β∫ηeηe+1ψiψjdη+∫ηeηe+1ψif¯∂ψj∂ηdη−∫ηeηe+1ψig¯ψjdη,Kij44=−1Pr1+4Rd3∫ηeηe+1∂ψi∂η∂ψj∂ηdη+∫ηeηe+1ψif¯∂ψj∂ηdη−γE∫ηeηe+1ψifg¯∂ψj∂ηdη−γE∫ηeηe+1ψiff¯∂2ψj∂η2dη+Nb∫ηeηe+1ψiϕ′¯∂ψj∂ηdη+Nt∫ηeηe+1ψiθ′¯∂ψj∂ηdη,K45=−Df∫ηeηe+1∂ψi∂η∂ψj∂ηdη,K41=K42=K43=K51=K52=K53=0,Kij54=−Sr+NtNb∫ηeηe+1∂ψi∂η∂ψi∂ηdη,Kij55=−1Sc∫ηeηe+1∂ψi∂η∂ψi∂ηdη+∫ηeηe+1ψif¯∂ψi∂ηdη−ε∫ηeηe+1ψiψjdη,ri1=0,ri2=−1+βψi∂g∂ηηeηe+1,ri3=−1+β2ψi∂p∂ηηeηe+1,ri4=−1Pr1+4Rd3ψi∂θ∂ηηeηe+1−Dfψi∂ϕ∂ηηeηe+1,ri5=−1Scψi∂ϕ∂η+Srψi∂θ∂ηηeηe+1,where(31)f¯′=∑j=13f¯j∂ψj∂η,g¯′=∑j=13g¯j∂ψj∂η,p¯′=∑j=13p¯j∂ψj∂η,θ¯′=∑j=13θ¯j∂ψj∂η,ϕ¯′=∑j=13ϕ¯′∂ψj∂η.


## 3.4. Assembly of the System of Equations

At this step, we assemble the elemental systems to obtain the global system, given as follows:(32)KY=F.This results in a large number of algebraic equations that govern the entire domain. And then, the global matrix will be modified by imposing the boundary conditions defined in equation (20). The last step is solving the assembled system of equations by the standard techniques like the Gauss elimination method, LU decomposition, Gauss Jordan method, or any iterative scheme.


## 4. Results and Discussion

The main target of the present study is to analyze the effects of thermal radiation, diffusion-thermo (Dufour), thermal-diffusion (Soret), chemical reaction, and Cattaneo-Christov model on nonlinear convective MHD micropolar couple stress nanofluids past a linearly stretching surface. The robust numerical method called the Galerkin finite element method (GFEM) is applied to solve the proposed model. We performed grid-invariance test or grid convergence test to confirm the convergence of the series solution. The impact of these pertinent parameters on velocity, angular velocity, temperature, concentration, local skin friction, local wall couple stress, local Nusselt number, and local Sherwood number was analyzed in both graphical and tabular forms. The default values of the present parameters used to plot the graphs are chosen based on the existing literature and parameter history and given below (Wubshet and Gosa [19]):(33)Pr=0.733,λ=0,2,ε=o.2,K=0.3,γE=0.5,Sc=0.9,M=0.5,Sr=o.5,Df=0.4,βt=βc=0.2,N∗=0.3,Rd=0.4,β=1.0,Nb=0.2,Nt=0.1.Figures2 and 3 are plotted to anticipate the control of radiation parameter Rd on linear velocity and temperature distributions. Enlarging radiation parameter in the boundary layer flow of micropolar couple stress nanofluid is to enhance the velocity of the fluid flow as revealed in Figure 2. Radiation in the boundary layer flow region rises as Rd increases; this in turn result enlargement in thermal boundary layer thickness as indicated in Figure 3. That is, to increase the radiation parameter is to initiate the temperature rise of the fluid flow. The Dufour parameter also produces similar effects on the velocity and temperature profiles of the laminar flow as illustrated in Figures 4 and 5. It is shown that the larger values of Df influenced the fluid to flow with faster speed and higher temperature. The enhancement of the Df caused increases in the concentration gradient which caused mass diffusion taking place more rapidly. In this circumstance, the rate of energy transfer associated with the particles became higher. That is why the temperature profile was boosted in the boundary layer flow region as plotted in Figure 5. The Dufour number has influenced insignificantly the concentration species to be lower in the laminar flow as indicated in Figure 6.Figure 2
Velocity distribution for different values of radiation parameter.Figure 3
Temperature distribution for different values of radiation parameter.Figure 4
Velocity distribution for different values of Dufour parameter.Figure 5
Temperature distribution for different values of Dufour parameter.Figure 6
Concentration distribution for different values of Dufour parameter.Figures7 and 8 notice the impacts of Soret number on the concentration profile and temperature profile, respectively. For larger values of Soret parameter, the concentration contour increases significantly whereas the temperature contour decreases insignificantly as noted in Figures 7 and 8, respectively. Figures 9–11 illustrated the impacts of the material parameter β on linear velocity, angular velocity, and temperature of the micropolar-couple stress fluid, respectively. It is revealed that increasing the material parameter β is to increase the linear velocity of the fluid and lower temperature in the flow region. In this circumstance, the angular velocity has not shown consistence as plotted in Figure 10. This result is in good agreement with the study reported by Wubshet Ibrahim and Gosa Gadisa [19]. As plotted in Figures 12 and 13, the couple stress parameter K has quite opposite effect on velocity and temperature of the boundary layer flow. The larger value in the couple stress parameter has the tendency to resist the fluid flow as noted in Figure 12 and caused higher temperature as shown in Figure 13. The control of mixed convection parameter λ on velocity, angular velocity, temperature, and concentration distributions is plotted in Figures 14–17. It is revealed that, in Figure 14, the enhancement in the mixed convection parameter initiates the fluid to flow more rapidly. This is due to the fact that the higher mixed convection parameter associates with the larger thermal buoyancy force which is responsible for the improvement of the linear velocity distribution and decline in temperature and concentration distribution contours (Figures 16 and 17). However, for the increasing values of mixed convection parameter, the angular velocity profile contour has not shown consistency as indicated in Figure 15. It is observed that the microrotation distribution very close to the sheet declines and at some distant from the sheet varies quite opposite with the larger values of mixed convection parameter. Figure 18 demonstrates the influence of first-order chemical reaction rate ε on concentration contour. It shows that concentration diminishes on elevated values of chemical reaction parameters. Figure 19 inspects that the temperature in the boundary layer region is higher in Fourier’s model than in the Cattaneo-Christov heat flux model.Figure 7
Concentration distribution for different values of Soret parameter.Figure 8
Temperature distribution for different values of Soret parameter.Figure 9
Velocity distribution for different values of material parameters.Figure 10
Angular velocity distribution for different values of material parameters.Figure 11
Temperature distribution for different values of material parameters.Figure 12
Velocity distribution for different values of couple stress parameters.Figure 13
Temperature distribution for different values of couple stress parameters.Figure 14
Velocity distribution for different values of mixed convection parameters.Figure 15
Angular velocity distribution for different values of mixed convection parameters.Figure 16
Temperature distribution for different values of mixed convection parameters.Figure 17
Concentration distribution for different values of mixed convection parameters.Figure 18
Concentration distribution for different values of chemical reaction parameters.Figure 19
Temperature distribution for different values ofγE.Figures20–23 are plotted to analyse the effects of radiation parameter Rd versus mixed convection parameter on local skin friction, local wall stress, local Nusselt number, and local Sherwood number, respectively. It is concluded that, for higher values of radiation parameter in the boundary layer region of the flow, both local skin friction and Local wall couple stress decrease insignificantly, while both heat and mass transfer significantly increase. Figure 24 illustrates the effects of Dufour number Df versus γE on heat transfer rate in the laminar flow. It is observed that, with the rise of Dufour number (diffusion-thermo), the heat transfer decrease for micropolar couple stress nanofluid and Soret number (thermal-diffusion) has a decreasing impact on the mass transfer rate as indicated in Figure 25.Figure 20
Local skin friction for different values ofRd versus λ.Figure 21
Local wall couple stress for different values ofRd versus λ.Figure 22
Local Nusselt number for different values ofRd versus λ.Figure 23
Local Sherwood number for different values ofRd versus λ.Figure 24
Local Nusselt number for different values ofDf versus γE.Figure 25
Local Sherwood number for different values ofSr versus ε.The grid-invariance test is performed to maintain the four-decimal point accuracy. It is also called the grid-invariance test or grid convergence test. We used this test to improve the results using a successively smaller step size for the calculations. We started by choosing a coarser mesh with 100 number of elements having step size ofh=0.1. Then, enhancing the number of elements ten times, we obtained a medium mesh with 1000 elements having a step size of h=0.01. Finally, we have a fine mesh of 1500 elements with step size of h=0.0067 and get four-decimal point-accuracy in velocity, angular velocity, temperature, and nanoparticle concentration values. After increasing the number of elements more than 1500, the accuracy is not affected but only to enlarge the compilation time. This is shown in Tables 1–4. Figures 26–29 are plotted to show the course, medium, and fine meshes for every fifth element of mesh.Table 1
Grid-independence test for velocity distribution.
f″ηηCoarse mesh with 100 elements (h=0.1)Medium mesh with 1000 elements (h=0.01)Fine mesh with 1500 elements (h=0.0067)1.51.092511.092541.092692.00.981870.981840.981842.50.931390.931350.931343.00.905800.905750.905743.50.891760.891690.891704.00.883590.883510.883504.50.878680.878570.878565.00.875670.875540.875535.50.873820.873660.873666.00.872700.872510.872516.50.872030.871810.871817.00.871640.871380.871387.50.871130.871130.87113Table 2
Grid-independence test for angular velocity distribution.
p′ηηCoarse mesh with 100 elements (h=0.1)Medium mesh with 1000 elements (h=0.01)Fine mesh with 1500 elements (h=0.0067)1.50.494070.494080.494172.00.412740.412710.412702.50.384520.384470.384473.00.373630.373580.373573.50.369090.369020.369024.00.367070.366970.366974.50.366130.366000.366005.00.365680.365520.365525.50.365460.365280.365286.00.365370.365150.365156.50.365340.365090.365097.00.365350.365050.365057.50.365030.365040.36503Table 3
Grid-independence test for temperature distribution.
θ′ηηCoarse mesh with 100 elements (h=0.4)Medium mesh with 1000 elements (h=0.01)Fine mesh with 1500 elements (h=0.0067)1.50.887700.887700.887702.00.676830.676830.676832.50.555360.555360.555363.00.479650.479650.479653.50.430720.430710.430714.00.398730.398730.398734.50.377890.377890.377895.00.364470.364460.364465.50.355940.355940.355946.00.350610.350610.350616.50.347340.347340.347347.00.345360.345360.345367.50.344180.344180.34418Table 4
Grid-independence test for concentration distribution.
ϕ′ηηCoarse mesh with 100 elements (h=0.1)Medium mesh with 1000 elements (h=0.01)Fine mesh with 1500 elements (h=0.0067)1.50.918680.918670.918672.00.837040.837030.837032.50.807620.807610.807613.00.799260.799240.799243.50.799400.799380.799384.00.802530.802510.802514.50.806220.806190.806195.00.809470.809440.809445.50.812000.811960.811966.00.813830.813780.813786.50.815080.815020.815027.00.815910.815840.815847.50.816360.816360.81636Figure 26
Grid-independence test showing every fifth element of the mesh for velocity profile.Figure 27
Grid-independence test showing every fifth element of the mesh for angular velocity profile.Figure 28
Grid-independence test showing every fifth element of the mesh for temperature profile.Figure 29
Grid-independence test showing every fifth element of the mesh for concentration profile.The assembled system of equations above is nonlinear in nature, and therefore, the iterative scheme is used to find the numerical solution. The system is linearized after incorporating the functionsf¯,g¯,p¯,θ¯,andϕ¯ which are expected to be known at the beginning of the iteration. The iterative process is completed or terminated when the following condition (convergence formula) is satisfied:(34)∑i,jχi,jn∗−χi,jn∗−1≤10−4,whereχdenotes eitherf,g,p,θ,andϕandn∗stands for the iterative step.Excellent convergence has been realized for all results, and it has been confirmed through the tabular and graphical forms for grid-invariance test mentioned above.Table5 indicates that our result is in good agreement with that in the existing literature, whereas Table 6 elaborates the effects with different parameters on local skin friction, local wall couple stress, heat, and mass transfer in the boundary layer flow region of micropolar couple stress fluid.Table 5
Comparison of the values of heat transfer rate−θ′0 forε=β=Κ=M=γE=λ=Sr=Df=Nt=Nb=0.
Pr[35][36][37]Present solution0.070.06630.06560.06560.06670.200.16910.16910.16910.16910.700.45390.45390.53490.45392.000.91130.91140.91140.91137.001.89541.89541.89051.8954Table 6
Numerical values of local skin friction coefficient−f″0, local wall couple stress −p′0, local Nusselt number –θ′0, and local Sherwood number −ϕ′0.
PrλΚγEMRdSrDfScβεNtNb−f″0−p′0–θ′0−ϕ′00.720.20.20.30.50.30.30.40.910.20.10.20.870830.365050.354200.773661.00————————————0.879750.409840.409840.786311.20————————————0.884460.368020.441530.793870.720.5———————————0.709650.284330.395610.80501—1.0———————————0.467980.152660.441750.84252—0.20.5——————————0.870830.365050.354200.77366——1.0——————————0.870830.365050.354200.77366——0.20.5—————————0.873260.365400.363150.77436———1.0—————————0.836170.327360.451940.68164———0.31.0————————1.023720.443990.319130.74550————1.5————————1.083940.449470.330950.65811————0.50.6———————0.863990.363760.404950.76530—————0.9———————0.858770.362850.450560.75959—————0.30.6——————0.870750.364980.349870.79075——————0.9——————0.870680.364910.345750.80706——————0.30.8—————0.860650.361900.222360.72602———————1.2—————0.851120.359020.108040.68622———————0.41.2————0.870740.364890.325510.88681————————1.5————0.870660.364760.300430.98566————————0.91.5———0.811750.316350.355270.82618—————————2.0———0.764020.279270.366530.83395—————————1.00.5——0.870770.364950.320690.90032——————————0.8——0.870700.364860.290191.01546———————————0.2—0.869910.364680.328360.83265———————————0.3—0.869010.364360.315590.84398———————————0.20.30.870120.364800.330360.80542————————————0.40.869420.364570.317420.79865


## 5. Conclusion

In this study, micropolar couple stress nanofluid flow past the stretching surface with the impact of relevant parameters is analyzed. The heat transfer in the boundary layer flow is modeled by the Cattaneo-Christov heat flux model. The robust numerical method called the Galerkin finite element method (GFEM) is applied to solve the proposed model. We performed grid-invariance test or grid convergence test to confirm the convergence of the series solution. The effect of numerous pertinent variables on velocity, angular velocity, temperature, concentration, local skin friction, local wall couple stress, local Nusselt number, and local Sherwood number is analyzed in both graphical and tabular forms, and the following remarks are forwarded:(1)
Both velocity and temperature distributions are increasing functions of radiation parameter and Dufour number.(2)
Chemical reaction and mixed convection parameters have a tendency to retard the concentrations of the species while Soret number revealed quite opposite effect.(3)
Material parameter and couple stress parameter effects are reversed on the velocity and temperature profiles.(4)
Heat and mass transfer in the flow region can be enhanced by boosting the radiation parameters.



---

*Source: 6683711-2021-01-15.xml*