# Numerical Treatment for Darcy-Forchheimer Flow of Sisko Nanomaterial with Nonlinear Thermal Radiation by Lobatto IIIA Technique

Authors: Iftikhar Uddin; Rizwan Akhtar; Zhu Zhiyu; Saeed Islam; Muhammad Shoaib; Muhammad Asif Zahoor Raja
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

In this study, the competency of numerical computational framework based on Lobatto IIIA technique is utilized for dynamical analysis of the Darcy-Forchheimer flow of Sisko nanomaterial with nonlinear thermal radiation. The resultant PDEs of the Sisko fluid model expressions are transformed into system of nonlinear ODEs by exploiting the similarity variables. Graphical representations and numerical illustrations are used to envisage the characteristics of various physical parameters of interest on velocity profile, nanoparticles concentration, and temperature distribution of Sisko fluidic system. In addition, skin friction and Nusselt number are numerically examined with observation that the material parameter of Sisko fluid increases the velocity profile as well as Nusselt number while decreasing temperature, concentration profiles, and skin friction coefficient.


---

## Body

## 1. Introduction

Material of the type having pores is termed as porous medium, in fluid dynamics; usually it is considered that the porous medium is filled with nanofluid. The study of nanofluid flow through porous medium has wide range of applications due to the fact that it has become a field of high interest for researchers, mathematicians, engineers, geologists, and architectures. Some of the considerable applications are flow of liquid in reservoirs, catalytic reactors, production of oil, heat exchanger, etc. Transfer of heat in porous medium has extensively been used in the heat pipe technology, nonwoven ingredients, production of papers, electronic based technology, storage of energy, etc. [1–4]. Henry Darc, in 1856, a French engineer, pioneered the flow of homogeneous liquids through porous medium during his valuable work on liquid flow through the sand beds. The law of Darcy is not effective for the situation when boundary and inertial effects appears at higher rate of fluid flow. After that the work was extended by Dutch scientist named Philippes Forchheimer in 1901. He extended the velocity equation by adding the square of velocity term in the equation of momentum of the problem to predict the behavior of boundary layer flow and inertia [5]. The term Forchheimer was later introduced by Muskat [6]. Pal and Mondal [7] constructed a fluid model with Darcy-Forchheimer effects over stretching sheet and he examined that with the decrease of the nanoparticle concentration profile, the value of electric field parameter increased. Ganesh et al. [8] studied the Darcy-Forchheimer MHD nanoliquid.Thermal radiation, either linear or nonlinear, has a very significant role in different high temperature processes in the field of fluid dynamics. It has diverse uses in several types of power plants and industries. For final possessions of desired qualities, the control of heat transfer is very significant and necessary in industries. Modern system linked to astrophysical flow, reactors cooling, electrical power generators, vehicles in the space, and plasmas are working by useful applications of radiation either linear or nonlinear. Latif et al. [9] studied the nonuniform viscosity fluid model with nonlinear radiation. Unsteady flow in existence of source/sink and radiation effects is modeled by Hayat et al. [10]. Pal [11] studied thermal radiation in a channel the peristaltic transport by using convective conditions. Pandey and kumar [12] studied the thermal radiation as well as mass flux in porous medium on natural convection. The effect of thermal radiation natural convection has been investigated by Huang [13]. Kothandapani and Prakash [14] worked out upon the influence of radiation over the peristalsis of magneto hydrodynamics in a channel. Some renowned researchers also worked on thermal radiation [15–23]. In these studies, the radiation with nonlinearity was introduced by the researchers with in fluid flow models. In these investigations extra variable contains difference between temperature of wall and ambient temperature. Due to high temperature variation, the transfer of radiation amongst the two objects is fortified. More details regarding nonlinear radiation are given [24–33]. Some recent work related to Non-Newtonian fluids, radiation, porous medium, and nanofluid is provided here [34–37]. These are the motivating factors for the authors to investigate in the model and analysis of nonlinear radiated nanofluid flow of Sisko fluid. The noticeable features of the study are listed as follows.(i)
An innovative design of Sisko nanofluid with nonlinear radiation and Darcy-Forchheimer flow involving porous medium has been modeled and PDEs of the model are transformed into an equivalent system of ODEs using suitable similarity transformations.(ii)
The concept of boundary layer is incorporated to develop the fluidic problem. The energy expression of the problem is developed by adding nonlinear radiation term.(iii)
Numerical analysis of the Sisko fluid model involving nanoparticles is presented by invoking the Lobatto IIIA method in order to investigate the effects of involved physical variables on velocity, concentration, and temperature by variation of physical quantities.(iv)
The application of Lobatto IIIA method to this model is an innovative work with these effects. Lobatto IIIA is a version of bvp4c method, which is based on finite difference method. Other methods can also solve this problem, i.e., NDSolve, RK-4, HAM, etc. this method has a strength to solve highly nonlinear systems of ODEs.(v)
Results presented in numerical and graphical representations endorsed the worth of the methodology as viable and effective computing solver for visualizing the dynamics of Sisko nanofluidic model.Outline of the rest of the paper is given as follows: mathematical model of the problem is represented in Section2, in Section 3, results and discussion of numerical experiments are provided, while the conclusion and future research work are presented in the last section.


## 2. Mathematical Development

Consider two-dimensional Sisko nanofluid flow through a porous medium and nonlinear stretching sheet as shown in Figure1. The coordinatey is taken as normal to the sheet while x-axis is along the sheet. A hot fluid is taken for heating or cooling the surface of sheet connectively. It is considered that uniformly concentrated fluid with nanoparticles is taken at the surface of the stretching sheet. The two dimensional governing equations with boundary layer are given below [38, 39]:(1)∂u∂x+∂v∂y=0,(2)u∂u∂x+υ∂u∂y=aρ∂2u∂y2-bρ∂∂y-∂u∂yn-vKu-Fu2,(3)u∂T∂x+υ∂T∂y=α∂2T∂y2+τDB∂C∂y∂T∂y+DTT∝∂T∂y2-1ρcp∂qr∂y,(4)u∂C∂x+υ∂C∂y=DB∂2C∂y2+DTT∝∂2T∂y2,with the boundary conditions(5)y=0,u=uwx=cxs,υ=0,-k∂T∂y=hfTf-T,C=Cw,y→∞:u→0,T→T∞,C→C∞,where the inertia coefficientF is defined by(6)F=CbxKRosseland radiation approximation is given as(7)qr=-4σ∗3k∗∂T4∂yTerm σ∗ in the above equation is Stefan-Boltzmann parameter while k∗ is mean absorption parameter. After simplifying (7), we have(8)∂qr∂y=-16σ∗3k∗∂T∞3∂y∂2T∂y2.In order to obtain ODEs from the above PDEs, we use the similarity transformations(9)u=cxsf′,υ=-UReb-1/n+11n+1s2n-1+1f+s2-n-1ηf′,η=yxReb1/n+1,ψ=xUReb-1/n+1fη,ϕη=C-C∞Cw-C∞,θη=T-T∞Tf-T∞.Continuity equation is identically satisfied, while (2)–(4) yield(10)Af′′′+n-f′′n-1f′′′+s2n+1+1n+1ff′′-sf′2-k1f′-Frf′2=0,(11)θ′′+Prs2n+1+1n+1fθ′+1+43Rdθ′′+43Rdθw-133θ′2θ2+θ3θ′′+3θw-122θ′2θ+θ2θ′′+3θw-1θ′2+θθ′′+Prλθ+Nbθ′ϕ′+Ntθ′2=0,(12)ϕ′′+PrLes2n-1+1n+1fϕ+NtNbθ′′=0,with boundary conditions(13)fη=0,f′η=1,θ′η=-γ1-θη,ϕη=1atη=0,f′η→0,θη→0,ϕη→0asη→∞,where parameters appearing in (9)-(16) are defined as(14)Rea=ρUxa,Reb=ρxnU2-nb,A=Reb2/n+1Rea,k1=avReaKρU2,λ=aQ0ReacpρU2,Le=αDB,Pr=aReaαρReb-2/n+1,γ=hkvc,Nb=τDBαCw-C∞,Fr=CbK1/2,Nt=DTT∞ταTf-T∞,α=kρcpfτ=ρCpsρCpf,Rd=4σT∞3k∗k∞,whereA is material parameter of Sisko fluid. The local Nusselt number is (15)Nux=xqwkTf-T∞,NuxReb-1/n+1=-θ′0,and qw is the heat flux, given by(16)qw=-k∂T∂yy=0Figure 1
Geometry representing the problem.


## 3. Results and Discussion

The transformed fluid flow problem defined in (10)–(13) has been solved numerically by exploiting the strength of Labotto IIIA method and invoking MATLAB. routine “bvp4c” for solution of boundary value problem of differential equation. We refer the interested reader to the article [40] for the necessary additional details of the Labotto III solver in [40] and citations therein. The solutions of (10)–(13) are also calculated with explicit Runge-Kutta (RK) method by using the Mathematica routine “NDSolve” for numerical solution of differential equations. The results of Lobatto III A are compared with explicit RK method and error analysis shows the matching of the order 10−08 to 10−09 exist between both methodologies which validate the correctness of the proposed numerical procedure. Therefore, in the rest of the study numerical procedure of Labotto IIIA is used for numerical solution of the problem.Characteristics of some influential parameters like heat Prandtl numberPr, the Froude number Fr, Brownian motion Nb parameter, thermophoresis number Nt, Lewis number Le, fitted rate constant n, the material parameter of Sisko fluid A, the radiation parameter Rd, the porosity permeability parameter k1, the heat generation/absorption variable λ, the Biot number γ, and θw the temperature ratio parameter on the velocity f′(η), temperature θ(η), and concentration ϕ(η) are shown. These numerical outcomes have been illustrated in Figures 2–10 and Tables 1–7. Numerical computation has been performed by choosing the default estimations of involved variables in the problem as Fr=n=A=λ=0.5, Pr=1.2, Le=1.5, Nb=Nt=0.2, s=1.5, k1=0.1, γ=1, Rd=0.9, and θw=0.7. Figures 2(a), 2(b), and 2(c) show the influence of Fr on f(η), f′(η), and ϕ(η). For increasing values of Fr, f(η) and f′(η) decrease while ϕ(η) increases. Figure 2(d) discloses the effect ofn on f(η). By increasing n, f(η) shows decreasing behavior. Figures 3(a), 3(b), and 3(c) depict the decreasing behavior of f′(η), θ(η), and ϕ(η), respectively, for increasing values of n. In Figure 3(d) increasing quality is illustrated for θ(η) for higher values of Prandtl number Pr. The concentration profile ϕ(η) shows decreasing characteristic for increasing Prandtl number Pr, Lewis number Le, and Brownian number Nb, which are shown in Figures 4(a), 4(b), and 4(d) respectively. Increase in Lewis Number Le has decreasing effect because the larger Lewis number makes diffusivity of molecules smaller and due to this reason, the concentration field is decreased. The temperature profile θ(η) is enhanced in Figure 4(c) for higher Nb. In Figures 5(a) and 5(b), it is clearly seen that increasing Nt causes increasing temperature θ(η) and concentration ϕ(η) profiles. Enhancement in temperature profile due to increase in thermophoresis number is because of nanoparticles. The nanoparticles which are closer to hot fluid shift towards cold fluid due to increase in thermophoresis parameter, which makes the temperature boundary layer thick. It is obvious in Figures 5(c) and 5(d) that increase in A causes increase in f(η) and f′(η), respectively. From Figures 6(a) and 6(b), it is clear that temperature θ(η) and concentration ϕ(η), respectively, both are decreasing for higher values of A. Also, Figures 6(c) and 6(d) both represent the decreasing behavior of f(η) and f′(η), respectively, for enhanced values of parameter s. Enhancement is seen in θ(η) and ϕ(η) for increase in s, which is shown in Figures 7(a) and 7(b), respectively. An increase in k1 results in decrease in f(η) and f′(η), which is visualized in Figures 7(c) and 7(d), respectively, and increasing effect of k1 is shown for θ(η) and ϕ(η) for increasing values of k1 in Figures 8(a) and 8(b), respectively. It is found in Figures 8(c) and 8(d) that θ(η) and ϕ(η) are enhanced, respectively, for higher estimations of γ. Figures 9(a), 9(b), 9(c), and 9(d) visualize the effect of λ on f(η), f′(η), θ(η), and ϕ(η), respectively. Higher estimations of λ cause decrease in f(η) and f′(η), while increase in θ(η) and dual behavior in concentration profile ϕ(η). In the domain from 0 to 2, ϕ(η) decreases while in the domain from 2 to 7 ϕ(η) it shows opposite behavior. Figure 10(a) shows that increase in Rd causes decrease in θ(η) and dual behavior of ϕ(η) is seen for higher estimations of Rd in Figure 10(b). In the domain from 0 to 2, ϕ(η) increases while in the domain from 2 to 7, ϕ(η) decreases. Figure 10(c) discloses the increasing behavior of θ(η) for higher values of the temperature ratio parameter θw on temperature distribution θ(η). This is because greater ratio in the temperature increases the temperature boundary layer. Figure 10(d) implies the dual characteristic of ϕ(η) for higher values of θw. In the domain from 0 to 2, ϕη diminishes while in domain from 2 to 7, ϕ(η) enhances. Table 1 shows the variation of values of involved parameters in different scenarios. Scenario 1 shows the variation of Fr; scenario 2, variation ofn; scenario 3, variation of Pr; scenario 4, variation of Le; scenario 5, variation of Nb; scenario 6, variation of Nt; scenario 7, variation of A; Scenario 8, variation of s; scenario 9, variation of k1; scenario 10, variation of γ; scenario 11, variation of λ; scenario 12, variation of Rd; and scenario 13, variation of θw. Table 2 shows skin friction -f′′(0) coefficients for various scenarios. Table 3 shows numerical values for surface heat transfer rate –θ′0 in different scenarios.Table 1
Parameters of physical quantities in the Sisko nano fluid model.
SCENARIOS Case-1 Case-2 Case-3 Case-4 1 Fr = 0.1 Fr = 0.3 Fr = 0.6 Fr = 0.9 2 n=0.5 n=0.6 n=0.7 n=0.8 3 Pr=1.0 Pr=1.2 Pr=1.4 Pr=1.6 4 Le =1.2 Le =1.4 Le =1.6 Le =1.8 5 Nb=0.7 Nb =0.9 Nb =1.3 Nb =1.8 6 Nt =0.2 Nt =0.4 Nt =0.6 Nt =0.8 7 A=0.2 A=0.4 A=0.6 A=0.8 8 s=1.1 s=1.3 s=1.5 s=1.7 9 k 1 = 0.1 k 1 = 0.3 k 1 = 0.5 k 1 = 0.7 10 γ = 1.0 γ = 1.2 γ = 1.4 γ = 1.6 11 λ = 0.2 λ = 0.4 λ = 0.6 λ = 0.8 12 R d = 0.9 R d = 1.2 R d = 1.5 R d = 1.8 13 θ w = 0.5 θ w = 0.6 θ w = 0.7 θ w = 0.8Table 2
Numerical data of skin friction-f′′(0) for the Sisko model.
SCENARIOS Case-1 Case-2 Case-3 Case-4 1 1.0902 1.1504 1.2367 1.3189 2 1.2085 1.17529 1.1501 1.1310 3 1.2085 1.2085 1.2085 1.2085 4 1.2085 1.2085 1.2085 1.2085 5 1.2085 1.2085 1.2085 1.2085 6 1.2085 1.2085 1.2085 1.2085 7 1.4636 1.2768 1.1508 1.0579 8 1.0902 1.1504 1.2085 1.2646 9 1.2085 1.2936 1.3747 1.4522 10 1.2085 1.2085 1.2085 1.2085 11 1.0902 1.1504 1.2085 1.2646 12 1.2085 1.2085 1.2085 1.2085 13 1.2085 1.2085 1.2085 1.2085Table 3
Computations for surface heat transfer rate (-θ′0) where: Fr=n=A=λ=0.5, Pr=1.2, Le=1.5, Nb=Nt=0.2, s=1.5, k1=0.1, γ=1, Rd=0.9 and θw=0.7.
SCENARIOS Case-1 Case-2 Case-3 Case-4 1 0.2984 0.2976 0.2965 0.3812 2 0.2968 0.3059 0.3134 0.3197 3 0.3100 0.2968 0.2824 0.2666 4 0.2976 0.2971 0.2966 0.2961 5 0.2968 0.2924 0.2837 0.2751 6 0.2968 0.2937 0.2906 0.2875 7 0.2940 0.2960 0.2976 0.2989 8 0.2984 0.2976 0.2968 0.2961 9 0.2968 0.2950 0.2934 0.2920 10 0.2968 0.3111 0.3221 0.3308 11 0.3571 0.3204 0.2687 0.1791 12 0.2968 0.3396 0.3672 0.3868 13 0.3818 0.3462 0.2968 0.2174Table 4
Analysis on the basis of number of meshes of the solver for Sisko model.
SCENARIOS Case-1 Case-2 Case-3 Case-4 1 866 889 911 1003 2 906 935 1092 1189 3 906 906 906 910 4 905 906 906 906 5 906 906 907 911 6 906 914 909 969 7 774 924 896 867 8 866 889 906 943 9 906 1012 1115 1223 10 906 906 906 906 11 866 889 906 998 12 906 906 906 906 13 906 906 906 909Table 5
Analysis on the basis of magnitude of relative errors of the solver for Sisko model.
SCENARIOS Case-1 Case-2 Case-3 Case-4 1 1.451e-06 2.990e-07 2.487e-07 2.247e-07 2 2.658e-07 4.473e-07 4.152e-07 7.453e-07 3 1.341e-07 2.658e-07 4.897e-07 2.028e-07 4 1.228e-07 2.017e-07 3.300e-07 5.433e-07 5 2.658e-07 3.305e-07 2.358e-07 1.002e-07 6 2.658e-07 1.228e-07 1.206e-07 1.728e-07 7 2.768e-07 2.266e-07 2.256e-06 6.363e-07 8 1.451e-06 2.990e-07 2.658e-07 2.325e-07 9 2.658e-07 2.110e-07 1.789e-07 2.631e-07 10 2.658e-07 2.464e-07 2.343e-07 2.304e-07 11 1.501e-06 2.761e-07 2.731e-07 1.425e-07 12 2.658e-07 2.469e-07 2.386e-07 2.304e-07 13 2.604e-07 2.431e-07 2.658e-07 9.887e-08Table 6
Analysis on the basis of number of ODEs evaluated of the solver for Sisko model.
SCENARIOS Case-1 Case-2 Case-3 Case-4 1 32260 34336 34666 36045 2 34591 35025 37379 38835 3 34591 34591 34590 34649 4 34576 34591 34590 34590 5 34591 34591 34606 34665 6 34591 34711 34635 35535 7 35267 34861 32650 32273 8 32260 34336 34591 35146 9 34591 36181 37726 39346 10 34591 34590 34590 34590 11 32260 34337 34591 35970 12 34591 34591 34591 34591 13 34590 34591 34591 34636Table 7
Analysis on the basis of number of boundary conditions evaluated of the solver for Sisko nano fluid model.
SCENARIOS Case-1 Case-2 Case-3 Case-4 1 97 98 98 98 2 98 98 98 98 3 98 98 98 98 4 98 98 98 98 5 98 98 98 98 6 98 98 98 98 7 118 98 97 97 8 97 98 98 98 9 98 98 98 98 10 98 98 98 98 11 97 98 98 98 12 98 98 98 98 13 98 98 98 98Figure 2
Variation of Fr onf, f′andϕ, while n on f.
(a)
Effect onf(η) by varying Fr (b)
Effect onf′(η) by varying Fr (c)
Effect onϕ(η) by varying Fr (d)
Effect onf(η) by varying nFigure 3
Variation of n onf′, θandϕ, while Pr on θ.
(a)
Effect onf′(η) by varying n (b)
Effect onθ(η) by varying n (c)
Effect onϕ(η) by varying n (d)
Effect onθ(η) by varying PrFigure 4
Variation ofPr on ϕ, Le on ϕ, Nb on θ, and Nb on ϕ.
(a)
Effect onϕ(η) by varying Pr (b)
Effect onϕ(η) by varying Le (c)
Effect onθ(η) by varying Nb (d)
Effect onϕ(η) by varying NbFigure 5
Variation ofNt on θ and ϕ, while that ofA on f and f′.
(a)
Effect onθ(η) by varying Nt (b)
Effect onϕ(η) by varying Nt (c)
Effect onf(η) by varying A (d)
Effect onf′(η) by varying AFigure 6
Variation of A onθ(η) and φ(η), while that of s on f(η) and f′(η).
(a)
Effect onθ(η) by varying A (b)
Effect onϕ(η) by varying A (c)
Effect onf(η) by varying s (d)
Effect onf′(η) by varying sFigure 7
Variation of s onθ(η) and φ(η), while that of k1 on f(η) and f′(η).
(a)
Effect onθ(η) by varying s (b)
Effect onϕ(η) by varying s (c)
Effect onf(η) by varying k1 (d)
Effect onf′(η) by varying k1Figure 8
Variation ofk1 on θ(η) and ϕ(η), while that of γ on θ(η) and ϕ(η).
(a)
Effect onθ(η) by varying k1 (b)
Effect onϕ(η) by varying k1 (c)
Effect onθ(η) by varying γ (d)
Effect onϕ(η) by varying γFigure 9
Variation ofλ on f(η), f′(η), θ(η), and ϕ(η).
(a)
Effect onf(η) by varying λ (b)
Effect onf′(η) by varying λ (c)
Effect onθ(η) by varying λ (d)
Effect onϕ(η) by varying λFigure 10
Variation of Rd onθandϕ while that of θw on θ(η) and φ(η).
(a)
Effect onθ(η) by varying Rd (b)
Effect onϕ(η) by varying Rd (c)
Effect onθ(η) by varying θw (d)
Effect onϕ(η) by varying θwError and complexity analysis of Lobatto IIIA solver is presented in Tables4–7. Table 4 explains the number of meshes for each case of different scenarios. Table 5 provides maximum error for each case of different scenarios, which show that reasonable accuracy attained by the scheme. Table 6 depicts the number of ODEs evaluated for each case of various scenarios. Table 7 explains the number of boundary conditions (BCs) evaluated during each case for different scenarios. The presented complexity indices can be used for comparison of the given scheme with future similar studies for the Sisko nanofluidic model.


## 4. Conclusions

Characteristics of radiated Sisko nanofluid are studied. Main features of the present analysis are listed below:(i)
The parameter A increasesf(η) and f′(η) while decreasing θ(η) and ϕ(η).(ii)
Velocity profilef′η decreases while concentration profile ϕ(η) increases with enhancement in s.(iii)
An increase ink1 shows decrease in f(η) and f′(η) while showing increase in θ(η) and ϕ(η).(iv)
Enhancement is observed inθ(η) and ϕ(η) for higher γ.(v)
A decrease inf(η) and f′(η) while increasing behavior is observed for higher values of λ.(vi)
Skin friction enhances for higher values ofFr, s, k1, and λ while diminishing for higher values of n and A.(vii)
Higher values ofFr, Pr, Le, Nb, Nt, s, k1, λ, and θw decay the heat transfer rate while reverse behavior is observed forn, A,γ, and Rd.



---

*Source: 8974572-2019-05-13.xml*