Comparison of Different Analytic Solutions to Axisymmetric Squeezing Fluid Flow between Two Infinite Parallel Plates with Slip Boundary Conditions

Abstract

We investigate squeezing flow between two large parallel plates by transforming the basic governing equations of the first grade fluid to an ordinary nonlinear differential equation using the stream functions 𝑢_𝑟(𝑟,𝑧,𝑡)=(1/𝑟)(𝜕𝜓/𝜕𝑧) and 𝑢_𝑧(𝑟,𝑧,𝑡)=−(1/𝑟)(𝜕𝜓/𝜕𝑟) and a transformation 𝜓(𝑟,𝑧)=𝑟²𝐹(𝑧). The velocity profiles are investigated through various analytical techniques like Adomian decomposition method, new iterative method, homotopy perturbation, optimal homotopy asymptotic method, and differential transform method.

1. Introduction

The study of squeezing flows has widespread applications in chemical engineering, industrial engineering, mechanical engineering, biomechanics, and food industry. Valves and diarthrodial joints are also the examples for squeeze flows relevant in biology and bioengineering. The first application to squeeze flow problem was made by Stefan in 1874 [1]. The motion of a thin film of lubricant, squeezed flow between two stationary parallel plane surfaces were reported by Tichy and Winner [2] and Wang and Watson [3]. The theoretical and experimental studies of squeezing flows have been conducted by many researchers [4–11]. The mathematical studies of these flows are concerned primarily with the nonlinear partial differential equations which arise from the Navier-Stokes equations. These equations have no general solutions, and only a few exact solutions have been attained by confining some physical aspects of the original problem [12]. To solve these nonlinear differential equations, different perturbation and analytical techniques have been extensively used in fluid mechanics and engineering [13].

In the literature only a few papers deal with the comparison of different analytical methods. In this paper we study the squeezed flow between two large parallel plates with slip boundary conditions. The velocity profile is obtained using various analytical techniques like Adomian decomposition method (ADM), new iterative method (NIM), homotopy perturbation (HPM), optimal homotopy asymptotic method (OHAM), and differential transform method (DTM) [14–23]. The residual of each technique is computed and a comparison is made to assess the efficiency of the above techniques. We select DTM for analyzing the velocity profile under different flow parameters.

Squeezing flows are produced by vertical movements of boundaries or by applying external normal forces. Commonly two types of boundary conditions are employed. For a viscous fluid at a solid wall, it is generally accepted that the fluid velocity matches the velocity of the solid boundary, and it is known as no slip boundary condition. While the no-slip condition is experimentally proven to be accurate for a number of macroscopic flows. Navier [24] proposed a general boundary condition that incorporates the possibility of fluid slip at a solid boundary. He assumed that the velocity𝑢_𝑥 at a solid surface is proportional to the shear rate at the surface, that is, 𝑢_𝑥=𝛽𝜕𝑢_𝑥/𝜕𝑦,where 𝛽is the slip length or slip coefficient. If𝛽=0, the generally assumed no-slip boundary condition is obtained, and if 𝛽 a finite constant, fluid is slip occurs at the wall. Its effect depends upon the length scale of the flow [25–28].

2. Basic Equation

We consider a steady axisymmetric flow where the velocity vector, 𝑢̃ is represented by𝑢̃=𝑢_𝑟(𝑟,𝑧),0,𝑢_𝑧(𝑟,𝑧).(2.1) In the absence of body forces, the Navier-Stokes equations are obtained for the first grade fluid by using equations of continuity and momentum,−𝜌𝑢̃×𝑤+∇∇⋅𝜌𝑢̃=0,2||𝑢||̃²+𝑝=−𝜂∇×𝑤,(2.2) where 𝜌is the constant density, 𝑝 is the pressure, 𝜂 is the viscosity, and 𝑤=∇×𝑢̃ is the vorticity vector.

Substituting (2.1) into (2.2), we get the following Continuity equation:𝜕𝑢_𝑟𝜕𝑟+𝑢𝑟^𝑟+𝜕𝜕𝑢𝑧^𝑧=0.(2.3)𝑟-component of N.S equation:𝜕𝜕𝑟𝜌2𝑢^2𝑟+𝑢^2𝑧+𝑝−𝜌𝑢_𝑧Ω(𝑟,𝑧)=−𝜂𝜕Ω(𝑟,𝑧)𝜕𝑧.(2.4)𝑧-component of N.S equation:𝜕𝜕𝑧𝜌2𝑢^2𝑟+𝑢^2𝑧+𝑝−𝜌𝑢_𝑧Ω(𝑟,𝑧)=𝜂1𝑟𝜕(𝑟Ω𝜕(𝑟𝑟,𝑧)),(2.5) where Ω(𝑟,𝑧)=𝜕𝑢_𝑟/𝜕𝑧−𝜕𝑢_𝑧/𝜕𝑟 is the vorticity function.

Now we define a function which is known as generalized pressure,𝑝̂=𝜌2𝑢^2𝑟+𝑢^2𝑧+𝑝.(2.6) Using (2.6), (2.4) and (2.5) take the following form:𝜕𝜕𝑟𝑝̂−𝜌𝑢^𝑧Ω(𝑟,𝑧)=−𝜂𝜕Ω(𝑟,𝑧)𝜕𝜕𝑧,𝜕𝑧𝑝̂−𝜌𝑢^𝑧Ω(𝑟,𝑧)=𝜂1𝑟𝜕(𝑟Ω𝜕(𝑟𝑟,𝑧)).(2.7) We now introduce the stream functions, 𝑢_𝑟(𝑟,𝑧,𝑡)=(1/𝑟)(𝜕𝜓/𝜕𝑧),𝑢_𝑧(𝑟,𝑧,𝑡)=(−1/𝑟)×(𝜕𝜓/𝜕𝑟),and obtain the following results: Ω(𝑟,𝑧)=−1𝑟𝐸²𝜓,(2.8) where 𝐸²=𝜕²/𝜕𝑟²−(1/𝑟)(𝜕/𝜕𝑟)+𝜕²/𝜕𝑧²,𝜕𝜕𝑟𝑝̂−𝜌𝑟1²𝜕𝜓𝜕𝑟𝐸²𝜓=−𝜂𝑟𝜕𝜕𝑧𝐸²𝜕𝜓,𝜕𝑧𝑝̂−𝜌𝑟1²𝜕𝜓𝜕𝑧𝐸²𝜓=−𝜂𝑟𝜕𝜕𝑟𝐸²𝜓.(2.9) Eliminating 𝑝̂, from (2.9), we obtain−𝜌𝜕𝜓,𝐸²𝜓/𝑟²𝜕(𝑟,𝑧)=𝜂𝑟𝐸⁴𝜓.(2.10) We consider viscous incompressible fluid, squeezed between two large planar and parallel plates, separated by a distance 2𝑑. The plates are moving towards each other with velocity 𝑈. The surfaces of both plates are covered by special material with slip length (slip coefficient)𝛽. For small values of 𝑈 the gape distance 2𝑑 between the plates varies slowly with the time 𝑡, so the flow can be taken as quasisteady [29–32]. See Figures 1, 2, 3, and 4.

Figure 1 

Plots of the residuals, blue: ADM, green: NIM, red: HPM, black: OHAM and dashed: DTM. In this figure, a comparison of the residuals is shown graphically, and it reveals that OHAM is more suitable technique for this problem. One can easily observe that the curve of the OHAM residual is approximately normally distributed while the curves of the residuals of other methods are 𝐽-shaped and highly skewed.

Figure 2 

𝑧-component 𝑢𝑧(𝑟,𝜃,𝑧)=(−1/𝑟)(𝜕𝜓/𝜕𝑟) of the velocity profile given by OHAM, when 𝑅=1,𝛾=1.

Figure 3 

𝑟-component 𝑢𝑟(𝑟,𝜃,𝑧)=(1/𝑟)(𝜕𝜓/𝜕𝑧)of the velocity profile given by OHAM, when 𝑅=1, 𝛾=1.

Figure 4 

Vorticity function Ω(𝑟,𝑧)=𝜕𝑢𝑟/𝜕𝑧−𝜕𝑢𝑧/𝜕𝑟 where the velocity components have been obtained by OHAM, for 𝑅=1, 𝛾=1.

The boundary conditions are as follows:

   (i)At𝑧=𝑑,𝑢_𝑟=𝛽𝜕𝑢^𝑟𝜕𝑧,𝑢^𝑧(ii)At𝑧=0,𝑢_𝑧=−𝑈.=0,𝜕𝑢^𝑟𝜕𝑧=0.(2.11) Now we use the transformation 𝜓(𝑟,𝑧)=𝑟²𝑓(𝑧).(2.12)

By virtue of 𝐸²=𝜕²/𝜕𝑟²−(1/𝑟)(𝜕/𝜕𝑟)+𝜕²/𝜕𝑧² and (2.12), the compatibility equation (2.10) and the boundary conditions equation (2.11) become𝐹^𝑖𝑣(𝑧)+2𝜌𝜂𝐹(𝑧)𝐹^′𝐹(0)=0,𝐹_(𝑧)=0,(0)=0,𝐹(𝑑)=𝑈2,𝐹^(𝑑)=𝛽𝐹^(𝑑).(2.13) Now by introducing dimensionless parameters𝐹^∗=𝐹𝑈/2,𝑧^∗=𝑧𝑑,𝛾=𝑑𝛽,𝑅=𝜌𝑑𝜂𝑈,(2.14) and dropping “*” for simplicity, the boundary value problem (2.13) become,𝑑⁴𝐹𝑑𝑧⁴+𝑅𝐹𝑑³𝐹𝑑𝑧³=0,(2.15) with boundary conditions𝐹(0)=0,𝐹^(0)=0,𝐹(1)=1,𝐹^(1)=𝛾𝐹^(1).(2.16)

3. Analytical Techniques

In this section, we give the basic idea of various analytical techniques and evaluate the velocity profile of our problem by considering 𝑅=1,𝛾=1. The residuals of all the techniques are computed and the results are displayed in Table 2.

3.1. Adomian Decomposition Method

According to [13, 14], we consider the differential equation𝐿(𝐹(𝑧))+𝑀(𝐹(𝑧))+𝑁(𝐹(𝑧))=𝑔(𝑧),(3.1) where 𝐿is the operator of the highest order derivative with respect to𝑧, 𝐿=𝑑⁴/𝑑𝑧⁴, 𝑅is the reminder of the linear term, and the nonlinear term is represented by𝑁(𝐹(𝑧)). Operating 𝐿⁻¹ on both sides of (3.1) we get the following: 𝐹(𝑧)=𝛼₀+𝛼₁𝑧+𝛼₂𝑧²2!+𝛼³𝑧³3!+𝐿⁻¹(𝑔(𝑧))−𝐿⁻¹𝑀(𝐹(𝑧))−𝐿⁻¹𝑁(𝐹(𝑧)),(3.2) where𝛼_𝑖:𝑖=0,1,2,3, the constants can be determined by using initial or boundary conditions.

The unknown function 𝐹(𝑧)can be expressed by an infinite series of the form𝐹(𝑧)=^∞_𝑛=0𝐹_𝑛(𝑧),(3.3) where𝐹₀(𝑧)=𝛼₀+𝛼₁𝑧+𝛼₂𝑧²/2!+𝛼₃𝑧³/3!+𝐿⁻¹(𝑔(𝑧)) and 𝐹_𝑛+1=−𝐿⁻¹𝑀(𝐹_𝑛(𝑧))−𝐿⁻¹(𝐴_𝑛), 𝑛=0,1,2,…. The nonlinear term 𝑁(𝐹(𝑧)) is decomposed by an infinite series of polynomial given by𝑁(𝐹(𝑧))=^∞_𝑛=0𝐴_𝑛,(3.4) where 𝐴_𝑛 are the so-called Adomian polynomials that can be determined by the formula𝐴_𝑛=1𝑛!𝑑𝑑𝜆^𝑛𝑛𝑁^𝑛_𝑖=0𝜆^𝑖𝐹_𝑖_𝜆=0.(3.5)

It has been observed that these polynomials can be constructed for a wide class of nonlinear functions.

The solution𝐹(𝑧)=∑^∞𝑛=0𝐹_𝑛(𝑧) is approximated by the truncated series of order𝐾, that is,𝐹(𝑧)=^𝐾_𝑛=0𝐹_𝑛(𝑧).(3.6) In our case𝑔(𝑧)=0, 𝑀(𝐹(𝑧))=0.

Now appling ADM on (2.15) and (2.16), for 𝑛=0,1,…,5.

We obtain 𝐹₀(𝑧)=𝐴𝑧³𝐹6+𝐵𝑧₁(𝑧)=−1120𝐴𝐵𝑧⁵+𝐴²𝑧⁷𝐹5040₂(𝑧)=𝐴𝐵²𝑧⁷1680+𝐴²𝐵𝑧⁹22680+𝐴³𝑧¹¹⋮.1108800(3.7)

Considering the Adomian 5th-order solution,𝐹(𝑧)=𝐹₀(𝑧)+𝐹₁(𝑧)+⋅⋅⋅+𝐹₅(𝑧)+𝑂𝑧²⁰.(3.8)

The boundary conditions at 𝑧=1 are used to get the following values of 𝐴and 𝐵.𝐴=1.74911,𝐵=0.71896.(3.9)

By substituting these values our solution is𝐹(𝑧)=0.718965𝑧+0.291518𝑧³−0.0104795𝑧⁵−0.0000688452𝑧⁷+0.0000701133𝑧⁹−3.49247×10⁻⁶𝑧¹¹−4.73949×10⁻⁷𝑧¹³+5.95015×10⁻⁸𝑧¹⁵−8.90909×10⁻¹⁰𝑧¹⁷−1.65991×10⁻⁹𝑧¹⁹+𝑂𝑧²⁰.(3.10)

3.2. New Iterative Method

The basic idea of new iterative method [15, 16]. Consider the following nonlinear general differential equation:𝐿(𝐹(𝑧))+𝑀(𝐹(𝑧))+𝑁(𝐹(𝑧))=𝑔(𝑧),(3.11) where 𝐿is the operator of the highest order derivative with respect to𝑧, 𝐿=𝑑⁴/𝑑𝑧⁴, 𝑀 is the reminder of the linear term, and the nonlinear term is represented by𝑁(𝐹(𝑧)). Operating 𝐿⁻¹ on both sides of (3.11) we get,𝐹(𝑧)=𝛼₀+𝛼₁𝑧+𝛼₂𝑧²2!+𝛼³𝑧³3!+𝐿⁻¹(𝑔(𝑧))−𝐿⁻¹𝑀(𝐹(𝑧))−𝐿⁻¹𝑁(𝐹(𝑧)),(3.12) where 𝛼_𝑖:𝑖=0,1,2,3, are constants to be determined by using initial or boundary conditions.

The unknown function 𝐹(𝑧)can be expressed by an infinite series of the form𝐹(𝑧)=^∞_𝑛=0𝐹_𝑛(𝑧),(3.13) where 𝐹₀(𝑧)=𝛼₀+𝛼₁𝑧+𝛼₂𝑧²/2!+𝛼₃𝑧³/3!+𝐿⁻¹𝑔(𝑧) and 𝐹_𝑛+1=−𝐿⁻¹𝑀(𝐹_𝑛(𝑧))−𝐿⁻¹(𝐺_𝑛), 𝑛=0,1,2,… The nonlinear term 𝑁(𝐹(𝑧)) is decomposed by an infinite series of polynomials given by𝑁(𝐹(𝑧))=^∞_𝑛=0𝐺_𝑛,(3.14) where 𝐺_𝑛=𝑁(∑^𝑛𝑘=0𝐹_𝑘)−𝑁(∑^{𝑛−1𝑘=0}𝐹_𝑘) and 𝐺₀=𝑁(𝐹₀(𝑧)).

The solution𝐹(𝑧)=∑^∞𝑛=0𝐹_𝑛(𝑧), is approximated by the truncated series of order 𝐾, that is,𝐹(𝑧)=^𝐾_𝑛=0𝐹_𝑛(𝑧).(3.15)

In our case 𝑔(𝑧)=0 and 𝑀(𝐹(𝑧))=0.

Now using NIM on (2.15) and (2.16), for 𝑛=0,1,…,5.

We obtain the following: 𝐹₀(𝑧)=𝐴𝑧³𝐹6+𝐵𝑧₁(𝑧)=−1120𝐴𝐵𝑧⁵+𝐴²𝑧⁷𝐹5040₂(𝑧)=−𝐴𝐵²𝑧⁷1680+𝐴²𝐵𝑧⁹22680+𝐴³𝑧¹¹1108800−𝐴²𝐵²𝑧¹¹1900800−𝐴³𝐵𝑧¹³38438400−𝐴⁴𝑧¹⁵⋮.3962649600(3.16) Considering the NIM 5th-order solution,𝐹(𝑧)=𝐹₀(𝑧)+𝐹₁(𝑧)+𝐹₃(𝑧)+𝐹₄(𝑧)+𝐹₅(𝑧)+𝑂𝑧²⁰,(3.17) we have the following:𝐹(𝑧)=𝐵𝑧+𝐴𝑧³6−1120𝐴𝐵𝑧⁵+−𝐴²5040+1𝐴6𝐵80²𝑧⁷+𝐴²𝐵22680−2𝐴41𝐵9³2𝑧⁹+𝐴³1108800−241𝐴²𝐵²39916800+38𝐴0𝐵16⁴0𝑧¹¹+−71𝐴³𝐵239500800+131𝐴²𝐵³207567360−65𝐴89𝐵4⁵40𝑧¹³+−1051𝐴⁴217945728000+1357𝐴³𝐵²25147584000−4759𝐴²𝐵⁴87178291200𝑧¹⁵+2243𝐴⁴𝐵1140023808000−359𝐴³𝐵³50523782400+179𝐴²𝐵⁵130288435200𝑧¹⁷+21919𝐴⁵844757641728000−71999𝐴⁴𝐵²167094918144000+104977𝐴³𝐵⁴300358272614400−31𝐴²𝐵⁶649813570560𝑧¹⁹+𝑂𝑧²⁰.(3.18) Using the boundary conditions at 𝑧=1, we get the following value of 𝐴 and 𝐵𝐴=1.74911,𝐵=0.71897.(3.19) The approximate solution is as follows:𝐹(𝑧)=0.718965𝑧+0.291518𝑧³−0.0104795𝑧⁵−0.0000688448𝑧⁷+0.0000701132𝑧⁹−3.49247×10⁻⁶𝑧¹¹−4.73948×10⁻⁷𝑧¹³+5.95015×10⁻⁸𝑧¹⁵−8.34647×10⁻¹¹𝑧¹⁷−1.18034×10⁻⁹𝑧¹⁹+𝑂𝑧²⁰.(3.20)

3.3. HPM

To illustrate the basic idea of homotopy perturbation method [17–21], we consider the following nonlinear differential equation:𝐴(𝐹)=𝑓(𝑟).(3.21)𝐴 is a general differential operator, the operator 𝐴 can usually be divided into two parts 𝐿and 𝑁, where 𝐿 is linear, and 𝑁 is nonlinear:𝐴=𝐿+𝑁.(3.22) so𝐿(𝐹)+𝑁(𝐹)−𝑓(𝑟)=0,𝑟∈Ω.(3.23)𝑓(𝑟)is a known analytic function.

With the boundary condition 𝐵(𝐹,𝜕𝐹/𝜕𝑛)=0,𝑟∈Γ.

𝐵 is a boundary operator, and Γis the boundary of the domain Ω.

Now we construct the following homotopy:𝐻(𝑣,𝑝)=(1−𝑝)𝐿(𝑣)−𝐿𝐹₀+𝑝[𝐴𝑣−𝑓(𝑟)],(3.24) where 𝑝∈[0,1]is an embedding parameter and 𝐹₀ is the first approximation that satisfied the boundary condition. To get an approximate solution, we expand 𝐹(𝑟,𝑝)in Taylor’s series about𝑝 in the following manner:𝐹(𝑟)=𝑣₀(𝑟)+^∞_𝑚₌₁𝑣_𝑚(𝑟)𝑝^𝑚.(3.25) Plugging (3.25) into (3.24) and then equating the coefficient of like powers of𝑝, we get the following problems which are directly integrable.

Zeroth-order problem:𝐹𝐹^0(𝑖𝑣)(𝑧)=0,₀(0)=0,𝐹^′0(0)=0,𝐹₀(1)=1,𝐹^0(1)=𝐹^0(1).(3.26)

First-order problem:𝐹𝐹^1𝑖𝑣(𝑧)=−𝐹₀(𝑧)𝐹^0(𝑧),₁(0)=0,𝐹^1(0)=0,𝐹₁(1)=0,𝐹^1(1)=𝐹^1(1).(3.27)

Second-order problem:𝐹𝐹^2𝑖𝑣(𝑧)=−𝐹₁(𝑧)𝐹^0(𝑧)−𝐹₀(𝑧)𝐹^1(𝑧),₂(0)=0,𝐹^2(0)=0,𝐹⋮₂.(1)=0,𝐹^2(1)=𝐹^2(1)(3.28) We consider the following 5th-order solution,𝐹𝐹((𝑧𝑧))==𝐹23₀7(𝑧9)6+40𝐹2₁17(𝑧7)8+09𝐹3₃91(𝑧6)4+04𝐹9₄𝑧(𝑧)+𝐹₅(𝑧)+𝑂𝑧²⁰,−3310695821245667419664353068210670𝑧0₇000+722438625316510010035703830881612280808090𝑧0³−44131018577416052800667574162308030𝑧0⁵−23493804016674220877𝑧28₁₃3200000+17132518713686749235965010𝑧0⁹0−883719021770748070𝑧0¹¹00−90921367129818598690𝑧0₁₉00+1227349972352173336𝑧00¹⁵000+15423904885192𝑧6¹⁷720018688959184896000+𝑂𝑧²⁰.(3.29)

3.4. OHAM

According to [22–24], we consider the following differential equation:𝐴(𝐹(𝑧))=𝑓(𝑟),𝑟∈Ω,(3.30) with the boundary conditions𝐵𝐹,𝜕𝐹(𝑧)𝜕𝑛=0,𝑟∈Γ,(3.31) where 𝐴is a differential operator, 𝐵is a boundary operator, and 𝑓(𝑟)is a known function of 𝑟: 𝑟∈Ω. The operator 𝐴 can be written as 𝐴=𝐿+𝑁, where𝐿 is linear and 𝑁is a nonlinear operator. In OHAM we first construct a homotopy equation,(1−𝑝)𝐿(𝐹(𝑟,𝑝))−ℎ(𝑝)[𝐿(𝐹(𝑟,𝑝))+𝑁(𝐹(𝑟,𝑝))−𝑓(𝑟)]=0,(3.32) where 𝑝∈[0,1]is an embedding parameter, ℎ(𝑝)is a nonzero auxiliary function for 𝑝≠0 andℎ(0)=0𝐿(𝐹(𝑧))=0,for𝑝=0𝐴(𝐹(𝑧))=𝑓(𝑟),for𝑝=1.(3.33) The solution𝐹(𝑟,0)=𝑣₀(𝑟) of 𝐿(𝐹(𝑧))=0 traces the solution curve 𝑣(𝑟) continuously as 𝑝 approaches to 1, where 𝑣₀ is the solution of the zeroth-order problem that will come in the next few lines. We next choose the auxiliary function ℎ(𝑝) in the following form:ℎ(𝑝)=^𝑚_𝑖=1𝑝^𝑖𝐶_𝑖,(3.34) where 𝐶₁,𝐶₂,…are the convergence controlling constants which are to be determined. To get an approximate solution, we expand 𝐹(𝑟,𝑝)in Taylor’s series about𝑝 in the following manner:𝐹𝑟,𝑝,𝐶_𝑖=𝑣₀(𝑟)+^∞_𝑚₌₁𝑟,𝐶₁,𝐶₁,…,𝐶_𝑚𝑝^𝑚.(3.35) Now after substituting the auxiliary function ℎ(𝑝) and 𝐹(𝑟,𝑝,𝐶_𝑖)in homotopy equation we compare the coefficient of like powers of𝑝, to obtain the following linear equations.

Zeroth-order problem:𝐿𝑣₀=0,𝐵𝑣₀,𝜕𝑣⁰𝜕𝑛=0.(3.36)

First-order problem:𝐿𝑣₁−1+𝐶₁𝑣₀𝐿+𝐶₁𝑓(𝑟)=𝐶₁𝑁₀𝑣₀,𝐵𝑣₁,𝜕𝑣¹𝜕𝑛=0.(3.37)

Second-order problem:𝐿𝑣₂−1+𝐶₁𝐿𝑣₁−𝐶₂𝐿𝑣₀−𝐶₂𝑓(𝑟)=𝐶₂𝑁₀𝑣₀+𝐶₁𝑁₁𝑣₁,𝐵𝑣₂,𝜕𝑣²𝜕𝑛=0,(3.38) and so on.

If the series is convergent at𝑝=1 for suitable auxiliary constants 𝐶₁,𝐶₂,…, then𝑣(𝑟)=𝐹𝑟,𝐶_𝑖=𝑣₀(𝑟)+^∞_𝑚₌₁𝑣_𝑚𝑟,𝐶₁,𝐶₁,…,𝐶_𝑚.(3.39)

The result of the mth-order approximations are given by𝑣̆(𝑟)=𝑣₀(𝑟)+^𝑚_𝑖=1𝑣_𝑖𝑟,𝐶₁,𝐶₁,…,𝐶_𝑖.(3.40)

Residual of the solution is𝑅𝑟,𝐶₁,𝐶₂,…,𝐶_𝑚=𝐴𝑣̆−𝑓(𝑟).(3.41) If𝑅=0, 𝑣̆ will be the exact solution, but it does not happen specially in nonlinear problems. To find the optimal values of𝐶_𝑖, many methods can be applied. We follow the method of least squares. According to the method of least squares, we first construct the functional𝐽𝑟,𝐶₁,𝐶₂,…,𝐶_𝑚=^𝑏𝑎𝑅²𝑑𝑟,(3.42) and then minimizing it, we have𝜕𝐽𝜕𝐶₁=𝜕𝐽𝜕𝐶₂=⋯=𝜕𝐽𝜕𝐶_𝑚=0,(3.43) where 𝑎and 𝑏are in the domain of the problem. With these constants known, the approximate solution (of order𝑚) is well determined.

Now applying OHAM to (2.15) and (2.16), we obtain the following problems which are directly integrable.

Zeroth-order problem:𝐹𝐹^0(𝑖𝑣)(𝑧)=0,₀(0)=0,𝐹^′0(0)=0,𝐹₀(1)=1,𝐹^0(1)=𝐹^0(1).(3.44)

First-order problem:𝐹𝐹^1𝑖𝑣(𝑧)=𝐶₁𝐹₀(𝑧)𝐹^0(𝑧)+𝐹^0𝑖𝑣(𝑧)+𝐶₁𝐹^0𝑖𝑣(𝑧),₁(0)=0,𝐹^1(0)=0,𝐹₁(1)=0,𝐹^1(1)=𝐹^1(1).(3.45)

Second-order problem:𝐹^2𝑖𝑣(𝐹𝑧)=𝐶₁𝐹₁(𝑧)𝐹^0(𝑧)+𝐶₁𝐹₀(𝑧)𝐹^1(𝑧)+𝐹^1𝑖𝑣(𝑧)+𝐶₁𝐹^1𝑖𝑣(𝑧),₂(0)=0,𝐹^2(0)=0,⋮𝐹.₂(1)=0,𝐹^2(1)=𝐹^2(1)(3.46) Considering the OHAM 5th-order solution and using the method of least squares, we obtain 𝐶₁=−0.93281. Hence the solution is 𝐹(𝑧)=0.718965𝑧+0.291517𝑧³−0.0104795𝑧⁵−0.0000688407𝑧⁷+0.0000700828𝑧⁹−3.4536×10⁻⁶𝑧¹¹−4.87073×10⁻⁷𝑧¹³+5.80616×10⁻⁸𝑧¹⁵+3.31269×10⁻⁹𝑧¹⁷−5.61766×10⁻¹⁰𝑧¹⁹+𝑂𝑧²⁰.(3.47)

3.5. NDSolve

NDSolve is a mathematica code, utilized for solution of ordinary and partial differential equations. This code is also used for differential-algebraic equations and system of ordinary differential equation. NDSolve gives solution on discrete points rather than for the function 𝐹 itself. List interpolation is used for the construction of approximating polynomial.

Apply NDSolve to (2.15) and (2.16), the following approximate solution is obtained:𝐹(𝑧)=0.718965𝑧−2.5833×10⁻⁶𝑧²+0.291529𝑧³−0.000025575𝑧⁴−0.0104737𝑧⁵+0.0000970056𝑧⁶−0.000294814𝑧⁷+0.000231011𝑧⁸−0.0000378459𝑧⁹+0.0000122062𝑧¹⁰.(3.48)

3.6. DTM

According to [33, 34], the basic idea of differential transforms method (DTM) starts from the following definition.

If 𝐹(𝑧) is a given function, its differential transform is defined as follows:𝐹(𝑟)=1𝑟!𝑑^𝑟𝐹(𝑧)𝑑𝑧^𝑟||||_𝑧=0.(3.49) The inverse transform of 𝐹(𝑟) is defined by𝐹(𝑧)=^∞_𝑟=0𝑧^𝑟𝐹(𝑟).(3.50) In actual application, the function 𝐹(𝑧) is expressed by a finite series𝐹(𝑧)=^𝑁_𝑟=0𝑧^𝑟𝐹(𝑟).(3.51) Equation (3.51) implies that 𝐹(𝑧)=∑^{∞𝑟=𝑁+1}𝑧^𝑟𝐹(𝑟). is negligibly small.

The fundamental operations of the DTM are given in Table 1.

Table 1 

Table 2 

Comparison of residuals. In this table residuals of all the techniques at the different domain points are displayed. We observe that OHAM provides more stable solution.

3.6.1. Analysis of the Method

Consider a fourth-order boundary value problem𝐹⁽⁴⁾(𝑧)=𝐺(𝑧,𝐹),0<𝑧<𝐻,(3.52) with the boundary conditions:𝐹(0)=𝛼₀,𝐹(𝐻)=𝛼₁,𝐹⁽¹⁾(𝐻)=𝛼₂,𝐹⁽²⁾(0)=𝛼₃,(3.53) where 𝛼_𝑖:𝑖=0,1,2,3 are given values.

The differential transform of (3.52) is as follows:𝐹(𝑟+4)=𝐺(𝑟)∏^4𝑖=1(𝑟+𝑖),(3.54) where 𝐺(𝑟) is the differential transform of𝐺(𝑧,𝐹).

The transformed boundary conditions (3.53) are given by𝐹(0)=𝛼₀,^𝑁_𝑟=0𝐻^𝑟𝐹(𝑟)=𝛼₁,^𝑁_𝑟=0𝑟𝐻^𝑟𝐹(𝑟)=𝛼₂,𝐹(2)=𝛼³2.(3.55) Using (3.54) and (3.55) values of 𝐹(𝑖)∶𝑖=4,5,…are obtained which give the following series solution up to 𝑂(𝑧^𝑁+1),𝐹(𝑧)=^𝑁_𝑟=0𝑧^𝑟𝐹(𝑟)+𝑂𝑧^𝑁+1,(3.56) applying DTM on (2.15) and (2.16), the transformed boundary conditions and differential are as follows: 𝐹(0)=0,𝐹(1)=𝑎,𝐹(2)=0,𝐹(3)=𝑏,𝐹(𝑟+4)=𝑟!(𝑟+4)−𝑅^𝑟_𝑘=0(𝑘+1)(𝑘+2)(𝑘+3)𝐹(𝑘+3)𝐹(𝑟−𝑘).(3.57) Using (3.57), we obtain the following values of 𝐹(𝑖)∶𝑖=1,2,3…,19.

For these values, the unknowns 𝑎and 𝑏 are determined by the following system:¹⁹_𝑟=0𝐹(𝑟)=1,¹⁹_𝑟=0𝑟𝐹(𝑟)=¹⁹_𝑟=0𝑟(𝑟−1)𝐹(𝑟).(3.58) We find, 𝑎=1.53266and𝑏=−0.570353.

Now using the inverse differential transform, the following approximate solution of 𝑂(𝑧²⁰)is obtained:𝐹(𝑧)=1.53266𝑧−0.570353𝑧³+0.043708𝑧⁵−0.00710856𝑧⁷+0.00130068𝑧⁹−0.000251907𝑧¹¹+0.0000492109𝑧¹³−9.6216×10⁻⁶𝑧¹⁵+1.86893×10⁻⁶𝑍¹⁷−3.60367×10⁻⁷𝑧¹⁹+𝑂𝑧²⁰.(3.59)

4. Conclusion

In this paper we have used Adomian’s decomposition method, new iterative method, homotopy perturbation method, optimal homotopy asymptotic method, and differential transform method to an axisymmetric squeezing flow problem. Though all the methods are based on Taylor’s series expansion, they produce different results because each one has its own environment. Homotopy methods combine homotopy from topology and perturbation method. It has been observed that homotopy perturbation method may suffer convergence in many problems while OHAM controls the convergence region by employing the auxiliary function. Differential transform method, Adomian decomposition method, and new iterative method are the straightforward application of Taylor’s series, and they suffer from divergence in general and particularly in initial value problems. Besides all these facts, ADM, HPM, and DTM can lead easily to closed-form solutions of many problems.