# The Solution of a Coupled Nonlinear System Arising in a Three-Dimensional Rotating Flow Using Spline Method

Authors: Jigisha U. Pandya
Journal: International Journal of Mathematics and Mathematical Sciences
 (2012)

---

## Abstract

The behavior of the non-linear-coupled systems arising in axially symmetric hydromagnetics flow between two horizontal plates in a rotating system is analyzed, where the lower is a stretching sheet and upper is a porous solid plate. The equations of conservation of mass and momentum are transformed to a system of coupled nonlinear ordinary differential equations. These equations for the velocity field are solved numerically by using quintic spline collocation method. To solve the nonlinear equation, quasilinearization technique has been used. The numerical results are presented through graphs, in which the effects of viscosity, through flow, magnetic flux, and rotational velocity on velocity field are discussed.


---

## Body

## 1. Introduction

In fluid mechanics, the problems associated with the flow that occurs due to the rotation of a single disk or that between two rotating disks have been found of interest of many researchers. Flows between finite disks were studied by Dijkstra and van Heijst [1], Adams and Szeri [2] and Szeri et al. [3]. Berker [4] showed that when the two disks are rotating with the same angular speed, there exists a one parameter family of solutions all but one of which is not rotationally symmetric. This result has been extended by Parter and Rajagopal [5], to disks rotating with differing angular speeds; they prove that the rotationally symmetric solutions are never isolated when considered within the full scope of the Navier-Stokes equations. The numerical study of the asymmetric flow has been carried out by Lai et al. [6, 7].Recently, hydromagnetic flow and heat transfer problems have become more important industrially. In view of these, Chakrabarti and Gupta [8] studied the hydromagnetic flow and heat transfer in a fluid, initially at rest and at uniform temperature, over a stretching sheet at a different uniform temperature. Banerjee [9] studied the effect of rotation on the hydromagnetic flow between two parallel plates where the upper plate is porous and solid, and the lower plate is a stretching sheet.In this paper, we analyze the behavior of the solution of the nonlinear coupled systems arising in axially symmetric hydromagnetic flow between two horizontal plates in a rotating system, where the lower plate is a stretching sheet. The governing coupled ordinary differential equations are solved by quintic spline collocation method.In Section2, the mathematical model of the problem given by Vajravelu and Kumar [10] is presented. The quintic spline collocation method is explained in Section 3. The results are displayed in graphical manner in Section 4, and the discussion of results is drawn in Section 5.


## 2. Formulation of the Problem

We consider the steady flow of an electrically conducting fluid between two horizontal parallel plates when the fluid and the plates rotate in unison about an axis normal to the plates with an angular velocityΩ. A Cartesian coordinate system is considered in such a way that the x-axis is along the plate, the y-axis is perpendicular to it, and the z-axis is normal to the xy-plane as shown in Figure 1.Figure 1
Flow configuration.The origin is located at the centre of the channel, and the plates are located aty=-h and h. The lower plate is is stretched by introducing two equal and opposite forces so that the position of the point (0,-h,0) remains unchanged. A uniform magnetic flux with density Bo is acting along y-axis about which the system is rotating. The upper plate is subjected to a constant wall injection with a velocity Vo. The equations of motion in a rotating frame of reference are
(2.1)u∂u∂x+υ∂u∂y+2Ωω=-1ρ∂p*∂x+ν[∂2u∂x2+∂2u∂y2]-σBo2ρu,(2.2)u∂υ∂y=-1ρ∂p*∂y+ν[∂2υ∂x2+∂2υ∂y2],(2.3)u∂ω∂x+υ∂ω∂y-2Ωu=ν[∂2ω∂x2+∂2ω∂y2]-σBo2ρω,(2.4)∂u∂x+∂υ∂y=0,
where u,υ, and ω denote the fluid velocity components along the x,y, and z directions, ν is the kinematics coefficient of viscosity, ρ is the fluid density, and p* is the modified fluid pressure. The velocity is independent of z and that all derivatives towards z therefore do not appear in the equations of motion. This leads to the absence of ∂p*/∂z in (2.3) which implies that there is a net cross-flow along the z-axis.The boundary conditions are(2.5)u=Ex,υ=0,ω=0aty=-h,u=0,υ=-υ0,ω=0aty=h.We introduce nondimensional variables(2.6)η=yh,u=Exf′(η),υ=-Ehf(η),ω=Exg(η),
where a prime denotes differentiation with respect to η.Substituting (2.6) in (2.1) to (2.4), we have
(2.7)-1ρ∂p*∂x=E2x[f′2-ff′′-f′′′R+M2Rf′+2K2Rg],(2.8)-1ρh∂p*∂η=E2h[ff′+1Rf′′],(2.9)g′′-R[f′g-fg′]+2K2f′-M2g=0,
where
(2.10)R=Eh2ν,viscosityparameter,M2=σB02h2ρν,magneticparameter,K2=Ωh2ν,therotationparameter.
Equation (2.7) with the help of (2.8) can be written as
(2.11)f′′′-R[f′2-ff′′]-2K2g-M2f′=A,
where A is a constant.Differentiation of (2.11) with respect to η gives
(2.12)f′′′′-R(f′f′′-ff′′′)-2K2g′-M2f′′=0.
Thus we solve the following nonlinear system numerically for several sets of values of the parameters:
(2.13)f′′′′-R(f′f′′-ff′′′)-2K2g′-M2f′′=0,g′′-R(f′g-fg′)+2k2f′-M2g=0,
subject to the boundary conditions
(2.14)f=0,f′=1,g=0atη=-1,f=λ,f′=0,g=0atη=1,
where λ=υ0/Eh a parameter depending on the y component of velocity at the upper plate.


## 3. Quintic Spline Collocation Method

The fifth degree spline is used to find numerical solutions to the boundary value problems discussed in (2.13) together with (2.14). A detailed description of spline functions generated by subdivision is given by de Boor [11].Consider equally spaced knots of a partitionπ: a=xo<x1<x2<⋯<xn=b on [a,b]. Let S5[π] be the space of continuously differentiable, piecewise, Quintic polynomials on π. That is, S5[π] is the space of Quintic polynomials on π. The Quintic spline is given by Bickley [12] and by G. Micula and S. Micula [13]
(3.1)s(x)=a0+b0(x-x0)+12c0(x-x0)2+16d0(x-x0)3+124e0(x-x0)4+1120∑k=0n-1Fk(x-xk)+5,
where the power function (x-xk)+ is defined as
(3.2)(x-xk)+={x-xk,ifx>xk,0,ifx≤xk,Consider a fourth-order linear boundary value problem of the form(3.3)yiv(x)+p(x)y′′′(x)+q(x)y′′(x)+r(x)y′(x)+t(x)y(x)=m(x),a≤x≤b,
subject to the boundary conditions
(3.4)α0y0+β0yn′+γ0yn′′+δ0yn′′′=η0,α1y0′+β1yn+γ1yn′′+δ1yn′′′=η1,α2y0′′+β2yn+γ2yn′′+δ2yn′′′=η2,α3y0′′′+β3yn+γ3yn′′+δ3yn′′′=η3,
where y(x),p(x),q(x),r(x),t(x),  and m(x) are continuous functions defined in the interval x∈[a,b];η0,η1,η2,η3 are finite real constants.Let (3.1) be an approximate solution of (3.3), where a0, b0, c0, d0, e0, F0, F1,…,Fn-1 are real coefficients to be determined.Letx0,x1,…,xn be n+1 grid points in the interval [a,b], so that
(3.5)xi=a+ih,i=0,1,…n;x0=a,xn=b,h=b-an.
It is required that the approximate solution (3.1) satisfies the differential equation at the point x=xi. Putting (3.1) with its successive derivatives in (3.3), we obtain the collocation equations as follows:
(3.6)∑k=0n-1Fk{(xi-xk)++12p(xi)(xi-xk)+2+16q(xi)(xi-xk)+3+124r(xi)(xi-xk)+4+1120t(xi)(xi-xk)+5}+e0{1+p(xi)(xi-x0)+12q(xi)(xi-x0)2+16r(xi)(xi-x0)3+124t(xi)(xi-x0)4}+d0{p(xi)+q(xi)(xi-x0)+12r(xi)(xi-x0)2+16t(xi)(xi-x0)3}+c0{q(xi)+r(xi)(xi-x0)+12t(xi)(xi-x0)2}+b0{r(xi)+t(xi)(xi-x0)}+a0{t(xi)}=m(xi),i=0,1,2,…,n.
From boundary conditions,
(3.7)∑k=0n-1Fk(δ02(b-xk)+2+γ06(b-xk)+3)+e0(δ0(b-a)+γ02(b-a)2)+d0(δ0+γ0(b-a))+c0(γ0)+b0(β0)+a0(α0)=η0,∑k=0n-1Fk(δ12(b-xk)+2+γ16(b-xk)+3+β1120(b-xk)+5)+e0(δ1(b-a)+γ12(b-a)2+β124(b-a)4)+d0(δ1+γ1(b-a)+β16(b-a)3)+c0(γ1+β12(b-a)2)+b0(β1(b-a)+α1)+a0(β1)=η1,∑k=0n-1Fk(δ22(b-xk)+2+γ26(b-xk)+3+β2120(b-xk)+5)+e0(δ2(b-a)+γ22(b-a)2+β224(b-a)4)+d0(δ2+γ2(b-a)+β26(b-a)3)+c0(γ2+β22(b-a)2+γ2(b-a)+α2)+b0(β2(b-a))+a0(β2)=η2,∑k=0n-1Fk(δ32(b-xk)+2+γ36(b-xk)+3+β3120(b-xk)+5)+e0(δ3(b-a)+γ32(b-a)2+β324(b-a)4)+d0(δ3+γ3(b-a)+β36(b-a)3+α3)+c0(γ3+β32(b-a)2)+b0(β3(b-a))+a0(β3)=η3.
Using the power function (x-xk)+ in the above equations, a system of n+5 linear equations in n+5 unknowns a0, b0, c0, d0, e0, F0, F1,…,Fn-1 is thus obtained. This system can be written in matrix-vector form as follows:
(3.8)AX=B,
where X=[Fn-1, Fn-2,…,F2, F1, F0, e0, d0, c0, b0, a0]T,  B=[η3, η2, η1, η0, m(xn), m(xn-1),…,m(x1), m(x0)]T.The coefficient matrixA is an upper triangular Hessenberg matrix with a single lower subdiagonal, principal and upper diagonal having nonzero elements. Because of this nature of matrix A, the determination of the required quantities becomes simple and consumes less time. The values of these constants ultimately yield the quintic spline s(x) in (3.1).In case of nonlinear boundary value problem, the equations can be converted into linear form using quasilinearization method (Bellman and Kalaba [14]), and hence this method can be used as iterative method. The procedure to obtain a spline approximation of yi(i=0,1,2,…,j, where j denotes the number of iteration) by an iterative method starts with fitting a curve satisfying the end conditions and this curve is designated as yi. We obtain the successive iterations yi’s with the help of an algorithm described as above till desired accuracy.


## 4. Quintic Spline Solution

We solve following nonlinear system numerically for several sets of values of the parameters:(4.1)f′′′′-R(f′f′′-ff′′′)-2K2g′-M2f′′=0,(4.2)g′′-R(f′g-fg′)+2K2f′-M2g=0.
subject to the boundary conditions
(4.3)f=0,f′=1,g=0atη=-1,f=λ,f′=0,g=0atη=1.The spline collocation method is used to solve the differentiation system (4.1) to (4.3). Equation (4.2) is a linear equation of order two, whereas (4.1) is a nonlinear equation of order four. For solving nonlinear equation by spline collocation method, we require a linear form of differentiation equation. The quasilinearization technique transforms (4.1) into linearized form as(4.4)fi+1iv+(Rfi)fi+1′′′-(Rfi′+M2)fi+1′′-(Rfi′′)fi+1′+(Rfi′′′)fi+1=2K2gi′+R(fifi′′′-fi′fi′′),
with boundary conditions
(4.5)fi+1(-1)=0,fi+1′(-1)=1,fi+1(1)=λ,fi+1′(1)=0.
The linear equation (4.2) can be written as
(4.6)g′′+Rfg′-(Rf′+M2)g=-2K2f′
with boundary conditions
(4.7)g(-1)=0,g(1)=0.
We use the quintic spline method to find solutions of (4.4) of fourth order with the conditions (4.5), whereas the cubic spline method is used for finding solutions of second-order equation (4.6) and (4.7).The quintic spline given by(4.8)s(ηi)=a0+b0(ηi-η0)+12c0(ηi-η0)2+16d0(ηi-η0)3+124e0(ηi-η0)4+1120∑k=0n-1Fk(ηi-ηk)+5
is an approximate solution of the problem given by (4.4) and (4.5). Substituting s(η) in both of above equations, we obtain the collocation as
(4.9)∑k=0N-1Fk[(ηi-ηk)++Rfi2(ηi-ηk)+2-(Rfi′+M2)6(ηi-ηk)+3-Rfi′′2(ηi-ηk)+4+Rf′′′120(ηi-ηk)+5]+e0[1+Rfi(ηi-η0)-(Rfi′+M2)2(ηi-η0)2-Rfi′′6(ηi-η0)3+Rfi′′′24(ηi-η0)4]+d0[Rfi-(Rfi′+M2)(ηi-η0)-Rfi′′2(ηi-η0)2+Rfi′′′6(ηi-η0)3]+c0[-(Rfi′+M2)-Rfi′′(ηi-η0)+Rfi′′′2(ηi-η0)2]+b0[Rfi′′′(ηi-η0)-Rfi′′]+a0[Rfi′′′]=2K2gi′+R(fifi′′′+fifi′′),i=0(1)N.
First two boundary conditions in (4.5) immediately give
(4.10)a0=0,b0=1,
and other two boundary conditions are
(4.11)1120∑k=0N-1Fk(ηN-ηk)+5+124e0(ηN-η0)4+16d0(ηN-η0)3+12c0(ηN-η0)2+b0(ηN-η0)+a0=λ,124∑k=0N-1Fk(ηN-ηk)+5+16e0(ηN-η0)3+12d0(ηN-η0)2+c0(ηN-η0)+b0=0.
First of all, initial assumptions regarding fi,fi′,fi′′,fi′′′ are necessary. Let a curve aη3+bη2+cη+d be fitted through the points η=-1 and η=1 satisfying the boundary conditions (4.5). This requirement yields a=(1-λ)/4, b=-1/4, c=(3λ-1)/4 and d=(1+2λ)/4 so that f(η)=((1-λ)/4)η3-(1/4)η2+((3λ-1)/4)η+((1+2λ)/4).Further, using the cubic spline(4.12)s(ηi)=a0+b0(ηi-η0)+12c0(ηi-η0)2+16∑k=0n-1dk(ηi-η0)+3.
In (4.6) and (4.7), the following set of equations is obtained:
(4.13)∑k=0N-1dk[(ηi-ηk)++0.1fi2(ηi-ηk)+2-(0.1fi′+16)(ηi-ηk)+3]+c0[1+0.1fi(ηi-η0)-(0.1fi′+12)(ni-η0)2]+b0[0.1fi-(0.1fi′+1)(ηi-η0)]+a0[-(0.1fi′+1)]=-2K2fi′,i=0(1)N,a0=0
and the last equation is
(4.14)16∑k=0N-1dk(ηN-ηk)+3+12c0(ηN-η0)2+b0(ηN-η0)+a0=0.A straight lineg(η)=aη+b can be fitted to the boundary conditions for g. A line g(η)=0 satisfied the boundary conditions (4.7). Therefore, initial values of g and g′ are also zero.ForN=10, R=0.1 and several sets of values of the parameters λ,M2 and K2, the velocity components f,f′, and g are obtained in graphical form. The behavior of velocity components are shown in Figures 2, 3, and 4. The results obtained by the spline collocation method are compared with the graphical solutions obtained by Vajravelu and Kumar [10], which are tested by comparing the analytical solutions (for small value of R) with the exact solutions f and g. The numerical results thus obtained by the spline method to find f,f′, and g are in close agreement with the available data. In Figures 2, 3, and 4, η equals a dimensionless injection parameter, M2 is a magnetic parameter, and K2 is a rotational parameter.Figure 2
Velocity profile for(R=0.1,N=10).Figure 3
Velocity profile forf′(R=0.1,N=10).Figure 4
Velocity profile forg (R=0.1, N=10).With improper initial guess, convergence is not observed every time. As discussed by Vajravelu and Kumar [10], this is due to inherent instability of boundary-value problem. Because of this, integration from the starting point of the domain may produce rapidly increasing solutions, which may occasionally lead to overflow before the end-point is reached.


## 5. Discussion of the Results

The behavior of the nondimensional velocity componentf for R=0.1 is shown in Figure 2. From Figure 2, it is evident that the Lorenz force decreases the velocity component f (see curves II and III), but the value of f increases with an increase in the value of the parameter λ as seen in curves III and IV. Further, in Figure 2, curves IV and V show that f decreases steadily for small K2, whereas for large value of K2, fincreases near the lower plate and decreases near the upper plate.Figure3 describes the behavior of f′(η), which is proportional to velocity along parallel plates, for several sets of values of the parameters only λ,M2 and K2 (the rotation parameter) are varied while R is kept constant at 0.1. From Figure 3, curves III and IV, it is evident that the velocity component f′ increases with the parameter λ(=υ0/Eh, y component of velocity at the upper plate), and the increase is maximized near the stretching plate. However, the effect of the Lorenz force (magnetic parameter M2) of f′ is to increase it near the upper plate and to decrease it near the lower (stretching) plate that is displayed from the curves II and III of Figure 3. Further, the effect of small K2 is to increase f′ at the upper plate and to decrease it near the lower plate is evident from curves I and II. But the effect of large K2 is to increase f′ near the plates and to decrease it at the centre of the channel. Also, the rotation of the channel brings humps near the plates, indicating the occurrence of a boundary layer near the plates. Furthermore, for large K2 the coriolis force and the magnetic field that act against the pressure gradient cause reversal of the flow.The transverse velocityg increases as the parameter λ increases, which is shown in curves III and IV of Figure 4. But this is quite opposite to the phenomenon with the magnetic parameter M2. However, the rotation parameter K2 increases g and the maximum of g occurs near the stretching sheet for large K2.From above figures, it is observed that for largeK2 convergence in the results is not guaranteed. This fact is also observed by Vajravelu and Kumar [10].We conclude from the above numerical experience that the spline collocation method can be successfully applied to solve the fourth-order nonlinear coupled equation, which governs the hydromagnetic fluid flow. This widens the field of applicability to the higher-order coupled differential equations. This encourages the applications to other types of problems.Another useful conclusion is that the selection of the domain is not restricted to positive interval only. That is we can successfully apply the above method for the negative intervals as the domain of the problem.



---

*Source: 702458-2012-08-15.xml*