Table of Contents Author Guidelines Submit a Manuscript
Discrete Dynamics in Nature and Society
Volume 2011, Article ID 285809, 19 pages
Research Article

Spectral Approximation of an Oldroyd Liquid Draining down a Porous Vertical Surface

1Department of Mathematics, University of Gaziantep, 27310 Gaziantep, Turkey
2Department of Elementary Mathematics Education, Mevlana University, 42003 Konya, Turkey
3Department of Mathematics, Nigde University, 51240 Nigde, Turkey

Received 29 June 2011; Revised 5 October 2011; Accepted 10 October 2011

Academic Editor: Carlo Piccardi

Copyright © 2011 F. Talay Akyildiz et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


Consideration is given to the free drainage of an Oldroyd four-constant liquid from a vertical porous surface. The governing systems of quasilinear partial differential equations are solved by the Fourier-Galerkin spectral method. It is shown that Fourier-Galerkin approximations are convergent with spectral accuracy. An efficient and accurate algorithm based on the Fourier-Galerkin approximations for the governing system of quasilinear partial differential equations is developed and implemented. Numerical results indicating the high accuracy and effectiveness of this algorithm are presented. The effect of the material parameters, elasticity, and porous medium constant on the centerline velocity and drainage rate is discussed.

1. Introduction

Thin-film drainage down porous vertical surfaces is important in industry. Draining films occur in processes as diverse as dip coating, electroplating, enameling, emptying storage vessels, and oil recovery mechanisms [1, 2]. Spectral projection and corresponding error analysis of the system of nonlinear partial differential equations arising in the free drainage start-up flow of Oldroyd four constant liquids over a porous vertical surface is considered.

Literature review reveals that this problem is not considered. But for the case of impermeable wall, Goshawk and Waters [3] and Pennington and Waters [4] investigated the drainage of an Oldroyd four-constant liquid from a vertical surface via a finite difference method. But the problem they consider is a special case of the expended investigation in this paper, and error analysis is not explored in their work. Again, for case of steady flow (or start-up phase neglected), the literature more richer, in this case, Keeley et al. [5] investigate the drainage of thin films of non-Newtonian liquids from vertical surface, and the behavior of the Phan Thain-Tanner models are investigated in detail [6, 7]. In the present study, Galerkin’s method of the system of quasilinear partial differential equations governing the free drainage problem is investigated for a porous vertical surface. It is shown that method converges and that the convergence is not at all dependent on whether or not the physical parameters of the problem assume special values. The paper is organized as follows. The problem is defined in Section 2, and some basic results on Fourier approximations are given. A suitable Fourier-Galerkin approximation for the problem under consideration is proposed in Section 3 and error analysis given following [811]. Efficient and robust algorithms for the problem under consideration are constructed and numerical results presented in Section 4.

2. Mathematical Formulation and Preliminary Results on Fourier Approximation

Consider a thin liquid film draining down a flat porous vertical surface defined by Cartesian coordinates (𝑥,𝑦,𝑧). The 𝑥 axis points vertically downwards, the solid surface lies in the plane 𝑦 = 0 with the thickness of the liquid film measured in the positive 𝑦 direction, and the 𝑧 axis is positioned perpendicular to the gravitational force completing a set of right-handed axes. The nondimensionalized equations of motion and the dimensionless Oldroyd four constant constitutive model form a quasilinear system of PDEs, where (details can be found in [3, 4, 12, 13] for the interested reader)𝜕𝑆𝑥𝑦=𝜕𝑦𝜕𝑢𝜕𝑡1+𝛼2𝑆𝑢,(2.1)𝑥𝑦+𝑆1𝜕𝑆𝑥𝑦+1𝜕𝑡2𝜇1𝑆𝑥𝑥𝜕𝑢=𝜕𝑦𝜕𝑢𝜕𝑦+𝑆2𝜕2𝑢𝑆𝜕𝑦𝜕𝑡,(2.2)𝑥𝑥+𝑆1𝜕𝑆𝑥𝑥𝜕𝑡2𝑆1𝑆𝑥𝑦𝜕𝑢𝜕𝑦=2𝑆2𝜕𝑢𝜕𝑦2,(2.3) where above we used the following dimensionless parameters as in [3]𝑔𝑦=𝑌𝜈21/3𝑔,𝑥=𝑥𝜈21/3𝑔,=𝐻(𝑥,𝑡)𝑦=𝑌𝜈21/3𝑔,𝑡=𝑇𝑦=𝑌2𝜈1/3,𝑆1=𝜆1𝑔2𝜈1/3,𝑆2=𝜆2𝑔2𝜈1/3𝑈,𝑢=(𝜈𝑔)1/3𝑉,𝑣=(𝜈𝑔)1/3,𝑆𝑖𝑗=𝑃𝑖𝑗𝜇𝑔𝜈1/3,𝜇1=𝜇0𝑔2𝜈1/3,𝛼2=(𝜈𝑔)1/3.𝜌𝑔(2.4) And here 𝑢=𝑢(𝑦,𝑡),𝑆𝑥𝑦=𝑆𝑥𝑦(𝑦,𝑡),𝑆𝑥𝑥=𝑆𝑥𝑥(𝑦,𝑡) are the dimensionless velocity and the dimensionless deviatoric stress tensor; 𝛼2,𝑆1,𝑆2, and 𝜇1 represent the porous medium constant, dimensionless relaxation and retardation time constants, and a dimensionless material parameter, respectively. No slip at the wall 𝑦=0 and zero shear rate on the free surface of the liquid are assumed,𝑢(0,𝑡)=0,𝑡0,𝜕𝑢𝜕𝑦=0on𝑦=.(2.5) The liquid is at rest at 𝑡=0, therefore initial conditions are𝑢(𝑦,0)=𝑆𝑥𝑦(𝑦,0)=𝑆𝑥𝑥(𝑦,0)=0.(2.6) To calculate the shape of the film profile at a given time 𝑡, the thickness is allowed to vary with 𝑥 while assuming the flow is still locally parallel. Combining the material derivative at the free surface𝑣(𝑥,,𝑡)=𝜕𝜕𝑡+𝑢(𝑥,,𝑡)𝜕,𝜕𝑥(2.7) with the equation of continuity yields a differential equation in ,𝜕𝜕𝑥0𝜕𝑢=𝜕𝑑𝑦+𝑢𝜕.𝜕𝑡(2.8) Introducing the flow rate 𝑄(,𝑡) per unit width across the film thickness ,𝑄(,𝑡)=0𝑢(𝑥,𝑦,𝑡)𝑑𝑦,(2.9) differentiating 𝑄 with respect to , and substituting the result into (2.8) and integrating give𝑥(,𝑡)𝑥0()=𝑡0𝜕𝑄𝜕𝑑𝜏,(2.10) where 𝑥0()=𝑥(,0) is the initial profile and can be chosen to represent any suitable initial shape. Equation (2.10) effectively determines the position of the free surface 𝑥 as a function of and 𝑡.

Next some mathematical notation is introduced. Denote the inner product in 𝕃2(0,) by(𝑓,𝑔)=0𝑓(𝑦)𝑔(𝑦)𝑑𝑦.(2.11) If 𝑓𝕃2(0,𝐻), then Fourier sine series is defined as𝑓(𝑦)=𝑘=1𝑓𝑠(𝑘)sin(2𝑘1)𝜋𝑦,2(2.12) where𝑓𝑠2(𝑘)=𝜋0sin(2𝑘1)𝜋𝑦2𝑓(𝑦)𝑑𝑦𝑘=1,2,.(2.13) Similarly, Fourier cosine series is defined as𝑓𝑓(𝑦)=𝑐(0)2+𝑘=1𝑓𝑐(𝑘)cos(2𝑘1)𝜋𝑦,2(2.14) where𝑓𝑐2(𝑘)=𝐻𝐻0cos(2𝑘1)𝜋𝑦2𝑓(𝑦)𝑑𝑦𝑘=1,2,.(2.15) Denote by 𝐻𝑚 the Sobolev norm, given by𝑓2𝐻𝑚=𝑘=1|||1+2𝑘12|||2𝑚||𝑓𝑠(||𝑘)2,||𝑓𝑐||(0)22+𝑘=1|||1+2𝑘12|||2𝑚||𝑓𝑐(||𝑘)2.(2.16)

The space of periodic Sobolev functions on the interval [0,]is defined as the closure of the space of smooth periodic functions with respect to the 𝐻𝑚-norm and will be simply denoted by 𝐻𝑚. In particular, the space 𝕃2(0,)with norm denoted by 𝕃2 is recovered for 𝑚=0. We now define subspaces of 𝕃2(0,)spanned by the set𝐷𝑁=2sin((2𝑘1)𝜋𝑦/2),𝐷,1𝑘𝑁𝑁=2cos((2𝑘1)𝜋𝑦/2).,1𝑘𝑁(2.17) The operators 𝑃𝑁 and 𝑃𝑁 denote the orthogonal, self-adjoint projection of 𝕃2 onto 𝐷𝑁 and 𝐷𝑁defined, respectively, by𝑃𝑁𝑓(𝑦)=𝑁𝑘=1sin(2𝑘1)𝜋𝑦𝑓2𝑠𝑃(𝑘),𝑁𝑓𝑓(𝑦)=𝑐(0)2+𝑁𝑘=1(cos2𝑘1)𝜋𝑦𝑓2𝑐(𝑘).(2.18) For𝑓𝐻𝑚, the estimates:𝑓𝑃𝑁𝑓𝐿2𝐶𝑝𝑁𝑚𝜕𝑚𝑦𝑓𝐿2,𝑓𝑃𝑁𝑓𝐻𝑛𝐶𝑝𝑁𝑛𝑚𝜕𝑚𝑦𝑓𝐿2,(2.19) hold for an appropriate constant 𝐶𝑝 and a positive integer 𝑛. The reader is referred to [8] for the proof of these inequalities.

The space of continuous functions from the interval [0,𝑇] into the space 𝐻𝑛is denoted by 𝐶([0,𝑇],𝐻𝑛). Similarly, we also consider the space 𝐶([0,𝑇],𝐷𝑁), where the topology on the finite-dimensional space 𝐷𝑁 can be given by any norm. Finally, note the inverse inequality𝜕𝑚𝑦𝜑𝐿2𝑁𝑚𝜑𝐿2,(2.20) which holds for integers 𝑚>0 and 𝜑𝐷𝑁. A proof of this estimate can also be found in [11]. We will make use of the Sobolev lemma, which guarantees the existence of a constant 𝑐 such thatsup𝑦||||𝑓(𝑦)𝑐𝑓𝐻1.(2.21)

We now note that exactly the same estimates hold for 𝑃𝑁. In the following, it will always be assumed that a solution of our problem (2.11)–(2.15) exists on some time interval[0,𝑇] with a certain amount of spatial regularity. In particular, we suppose that a solution exists in the (𝐶([0,𝑇],𝐻1))3 space for some𝑇>0. With these preliminaries in place, we are now set to tackle the problem of defining a suitable spectral projection of (2.11)–(2.15) and proving the convergence of such a projection. First, the Fourier-Galerkin method is presented and a proof of convergence given.

3. The Fourier-Galerkin Method

{𝑒𝑘(𝑦),𝑘𝜖}={2sin((2𝑘1)𝜋𝑦/2)/,𝑘} and {𝑓𝑘(𝑦),𝑘𝜖}={2cos((2𝑘1)𝜋𝑦/2)/,𝑘} are chosen to be an orthonormal basis of the Hilbert space𝕃20[0,] and 𝕃2[0,], respectively. Then, the subspace of these Hilbert spaces spanned by the 𝐷𝑁={2sin((2𝑘1)𝜋𝑦/2)/,1𝑘𝑁} and𝐷𝑁={2cos((2𝑘1)𝜋𝑦/2)/,0𝑘𝑁},respectively. Fourier-Galerkin approximation of (2.1)–(2.5) are find the functions 𝑢𝑁(𝑡)𝐷𝑁, 𝑆𝑁𝑥𝑦(𝑡), and 𝑆𝑁𝑥𝑥𝐷𝑁 for all 0𝑡𝑇, such that𝜕𝑡𝑢𝑁1𝜕𝑌𝑆𝑁𝑥𝑦,𝜔1+𝛽2𝑢,𝜔1[],𝑆=0,𝑡0,𝑇(3.1)𝑁𝑥𝑦+𝑆1𝜕𝑡𝑆𝑁𝑥𝑦+12𝜇1𝑆𝑁𝑥𝑥𝜕𝑦𝑢𝑁𝜕𝑦𝑢𝑁𝑆2𝜕2𝑡𝑦𝑢𝑁,𝜔2[]𝑆=0,𝑡0,𝑇,(3.2)𝑁𝑥𝑥+𝑆1𝜕𝑡𝑆𝑁𝑥𝑥2𝑆1𝑆𝑁𝑥𝑦𝜕𝑦𝑢𝑁2𝑆2𝜕𝑦𝑢𝑁2,𝜔2[]𝑢=0,𝑡0,𝑇,(3.3)𝑁(0)=0,𝑆𝑁𝑥𝑦(0)=0,𝑆𝑁𝑥𝑥(0)=0,(3.4) for all 𝜔1𝐷𝑁 and for all 𝜔2𝐷𝑁. Since for each 𝑡, 𝑢𝑁(,𝑡), 𝑆𝑁𝑥𝑦(,𝑡), and 𝑆𝑁𝑥𝑦(,𝑡) have the form𝑢𝑁(𝑦,𝑡)=𝑁𝑘=1̂𝑢𝑁(𝑘,𝑡)2sin((2𝑘1)𝜋𝑦/2),𝑆𝑁𝑥𝑦(𝑦,𝑡)=𝑁𝑘=1̂𝑠𝑁𝑥𝑦(𝑘,𝑡)2cos((2𝑘1)𝜋𝑦/2),𝑆𝑁𝑥𝑥(𝑦,𝑡)=𝑁𝑘=1̂𝑠𝑁𝑥𝑥(𝑘,𝑡)2cos((2𝑘1)𝜋𝑦/2).(3.5)

Taking 𝜔1=2/𝐻sin((2𝑘1)𝜋𝑦/2),𝜔2=2/𝐻cos((2𝑘1)𝜋𝑦/2),1𝑘𝑁in (3.1)–(3.3) yields the following system of equations for the Fourier coefficients of 𝑢𝑁, 𝑆𝑁𝑥𝑦, and 𝑆𝑁𝑥𝑥:𝑑𝑑𝑡̂𝑢𝑁(𝑘,𝑡)=(2𝑘1)𝜋2̂𝑠𝑁𝑥𝑦(𝑘,𝑡)+22𝜋(2𝑘1),(3.6)̂𝑠𝑁𝑥𝑦(𝑘,𝑡)+𝑆1𝑑𝑑𝑡̂𝑠𝑁𝑥𝑦(𝑘,𝑡)=(2𝑘1)𝜋2̂𝑢𝑁(𝑘,𝑡)+𝑆2(2𝑘1)𝜋𝑑2𝑑𝑡̂𝑢𝑁(𝑘,𝑡)+𝜇1(2𝑘1)𝜋42𝑁𝑖,𝑗=1𝑐𝑖𝑗𝑁̂𝑠𝑁𝑥𝑥(𝑖,𝑡)̂𝑢𝑁(𝑘,𝑡),(3.7)̂𝑠𝑁𝑥𝑥(𝑘,𝑡)+𝑆1𝑑𝑑𝑡̂𝑠𝑁𝑥𝑥(𝑘,𝑡)(2𝑘1)𝜋2̂𝑢𝑁(𝑘,𝑡)2𝑆1(2𝑘1)𝜋42𝑁𝑖,𝑗=1𝑐𝑖𝑗𝑁̂𝑠𝑁𝑥𝑦(𝑖,𝑡)̂𝑢𝑁(𝑘,𝑡)=2𝑆2(2𝑘1)𝜋422𝑁𝑖,𝑗=1𝑑𝑖𝑗𝑁̂𝑢𝑁(𝑖,𝑡)̂𝑢𝑁(𝑗,𝑡),(3.8)̂𝑢𝑁(𝑘,𝑡)=0,̂𝑠𝑁𝑥𝑦(𝑘,𝑡)=0,̂𝑠𝑁𝑥𝑥𝑐(𝑘,𝑡)=0,(3.9)𝑖𝑗𝑁=0cos(2𝑖1)𝜋𝑦2sin(2𝑗1)𝜋𝑦2cos(2𝑁1)𝜋𝑦𝑑2𝑑𝑦,𝑖𝑗𝑁=0sin(2𝑖1)𝜋𝑦2sin(2𝑗1)𝜋𝑦2cos(2𝑁1)𝜋𝑦2𝑑𝑦.(3.10)

This is a nonlinear system of ordinary differential equations for the functions {̂𝑢𝑁(𝑘,𝑡),̂𝑠𝑁𝑥𝑦(𝑘,𝑡),̂𝑠𝑁𝑥𝑥(𝑘,𝑡)}𝑁𝑘=1; by standard existence theory, there is a unique solution which exists on some time interval[0,𝑇𝑁), where 𝑇𝑁 possibly may be equal to𝑇. Since the argument is standard, the proof is omitted here. The main result of this paper is the fact that the Galerkin approximation {𝑢𝑁,𝑆𝑁𝑥𝑦,𝑆𝑁𝑥𝑥} converges to the exact solution {𝑢,𝑆𝑥𝑦,𝑆𝑥𝑥} when 𝑢 is smooth enough. This is stated in the next theorem.

Theorem 3.1. Suppose that a solution {𝑢,𝑆𝑥𝑦,𝑆𝑥𝑥} of (2.1)–(2.5) exists in the space (𝐶([0,𝑇],𝐻𝑚))3 for 𝑚1 and for some time 𝑇>0. If 𝑃𝑁𝑢(,𝑡)𝑢𝑁(𝑢,𝑡)𝑐𝑁1𝑚 and 𝑃𝑁𝑆𝑥𝑦𝑆𝑁𝑥𝑦𝑐𝑁1𝑚, then, for large enough 𝑁, there exists a unique solution {𝑢𝑁,𝑆𝑁𝑥𝑦,𝑆𝑁𝑥𝑥} of the finite dimensional problem (3.1)–(3.4). Moreover, there exist constants Γ1,Γ2, and Γ3 such that sup[]𝑡0,𝑇𝑢𝑢𝑁𝐿2Γ1𝑁1𝑚,sup[]𝑡0,𝑇𝑆𝑥𝑦𝑆𝑁𝑥𝑦𝐿2Γ2𝑁1𝑚,sup[]𝑡0,𝑇𝑆𝑥𝑥𝑆𝑁𝑥𝑥𝐿2Γ3𝑁1𝑚.(3.11)

Before the proof is given, note that the assumptions of the theorem encompass the existence of constants 𝜅, 𝜅1, and 𝜅2 such thatsup[]𝑡0,𝑇(𝑢𝑦,𝑡)𝐻𝑚𝜅,sup[]𝑡0,𝑇𝑆𝑥𝑥(𝑦,𝑡)𝐻𝑚𝜅1,sup[]𝑡0,𝑇𝑆𝑥𝑦(𝑦,𝑡)𝐻𝑚𝜅2.(3.12) In particular, it follows that there are other constants 𝜓, 𝜓1, and 𝜓2 such thatsup[]𝑡0,𝑇(𝑢𝑦,𝑡)𝐻𝑚𝜓,sup[]𝑡0,𝑇𝑆𝑥𝑥(𝑦,𝑡)𝐻𝑚𝜓1,sup[]𝑡0,𝑇𝑆𝑥𝑦(𝑦,𝑡)𝐻𝑚𝜓2.(3.13) The main ingredient in the proof of the theorem is a local error estimate which will be established by the following lemma.

Lemma 3.2. Suppose that the solution {𝑢𝑁,𝑆𝑁𝑥𝑦,𝑆𝑁𝑥𝑥} of (3.1)–(3.4) exists on the time interval [0,𝑡𝑁] and that sup𝑡[0,𝑡𝑁]𝑢𝑁(𝑦,𝑡)𝐻22𝜓, sup𝑡[0,𝑡𝑁]𝑆𝑁𝑥𝑦(𝑦,𝑡)𝐻22𝜓1, sup𝑡[0,𝑡𝑁]𝑆𝑁𝑥𝑥(𝑦,𝑡)𝐻22𝜓2, sup𝑡[0,𝑡𝑁]|𝑃𝑁𝑢(,𝑡)𝑢𝑁(,𝑡)|𝛽𝑁1𝑚, and sup𝑡[0,𝑡𝑁]|𝑃𝑁𝑆𝑥𝑦(,𝑡)𝑆𝑁𝑥𝑦(,𝑡)|𝛽1𝑁1𝑚, then the error estimate: sup𝑡0,𝑡𝑁𝑢𝑢𝑁𝐿2Γ1𝑁1𝑚,sup𝑡0,𝑡𝑁𝑆𝑥𝑦𝑆𝑁𝑥𝑦𝐿2Γ2𝑁1𝑚,sup𝑡0,𝑡𝑁𝑆𝑥𝑥𝑆𝑁𝑥𝑥𝐿2Γ3𝑁1𝑚,(3.14) holds for constant Γ1 which is the function of 𝑇,𝛼2,𝐶𝑝,𝜅,𝛽, and 𝑐, constantΓ2 which is the function of 𝑇,𝑆1,𝑆2,𝜅1,𝜓,𝜓1,𝜓2,𝑐1,𝛽1, and𝑐, and constant Γ3 which is the function of 𝑇,𝑆1,𝜅1,𝜓,𝜓1,𝑐1,𝜅, and 𝑐.

Proof. Let𝜎1=𝑃𝑁𝑢𝑢𝑁,𝜎2=𝑃𝑁𝑆𝑥𝑦𝑆𝑁𝑥𝑦,𝜎3=𝑃𝑁𝑆𝑥𝑥𝑆𝑁𝑥𝑥. Also, from the definition of 𝑃𝑁 and 𝑃𝑁,we have 𝜕𝑦𝜎1=𝑃𝑁𝜕𝑦𝑢𝜕𝑦𝑢𝑁,𝜕𝑦𝜎2=𝑃𝑁𝜕𝑦𝑆𝑥𝑦𝜕𝑦𝑆𝑁𝑥𝑦,𝜕𝑦𝜎3=𝑃𝑁𝜕𝑦𝑆𝑥𝑥𝜕𝑦𝑆𝑁𝑥𝑥. We apply𝑃𝑁,𝑃𝑁 and 𝑃𝑁 to both sides of (2.1)–(2.3), respectively. Since 𝑃𝑁, 𝑃𝑁 commute with derivation, we obtain 𝜕𝑡𝑃𝑁𝑢+𝛼2𝑢=1+𝜕𝑦𝑃𝑁𝑆𝑥𝑦,𝑃𝑁𝑆𝑥𝑦+𝑆1𝜕𝑡𝑃𝑁𝑆𝑥𝑦+12𝜇1𝑃𝑁𝑆𝑥𝑥𝜕𝑦𝑢=𝑃𝑁𝜕𝑦𝑢+𝑆2𝜕𝑡𝑃𝑁𝜕𝑦𝑃𝑢,𝑁𝑆𝑥𝑥+𝑆1𝜕𝑡𝑃𝑁𝑆𝑥𝑥2𝑆1𝑃𝑁𝑆𝑥𝑦𝜕𝑦𝑢=2𝑆2𝑃𝑁𝜕𝑦𝑢2.(3.15) We multiply these equations with test functions 𝜎1𝐷𝑁,𝜎2𝐷𝑁, and 𝜎3𝐷𝑁, respectively, integrate over [0,], and subtract the resulting expressions from (3.1), (3.2), and (3.3) to get 12𝑑𝜎𝑑𝑡12𝐿2+𝛼2𝜎12𝐿2=𝑃𝑁11𝑁,𝜎1+𝑃𝑁𝜕𝑦𝑆𝑥𝑦𝑆𝑁𝑥𝑦,𝜎1,𝜎22𝐿2+𝑆112𝑑𝜎𝑑𝑡22𝐿2+12𝜇1𝑃𝑁𝑆𝑥𝑥𝜕𝑦𝑢𝑆𝑁𝑥𝑥𝜕𝑦𝑢𝑁,𝜎2=𝑃𝑁𝜕𝑦𝑢𝑢𝑁,𝜎2+𝑆2𝑃𝑁𝜕2𝑡𝑦𝑢𝑢𝑁,𝜎2,𝜎22𝐿2+𝑆112𝑑𝜎𝑑𝑡22𝐿2+12𝜇1𝑃𝑁𝑆𝑥𝑥𝜕𝑦𝑢𝑆𝑁𝑥𝑥𝜕𝑦𝑢𝑁,𝜎2=𝑃𝑁𝜕𝑦𝑢𝑢𝑁,𝜎2+𝑆2𝜕𝑡𝑃𝑁𝜕𝑦𝑢𝑢𝑁,𝜎2,𝜎32𝐿2+𝑆112𝑑𝜎𝑑𝑡32𝐿22𝑆1𝑃𝑁𝑆𝑥𝑦𝜕𝑦𝑢𝑆𝑁𝑥𝑦𝜕𝑦𝑢𝑁,𝜎3=2𝑆2𝑃𝑁𝜕𝑦𝑢2𝜕𝑦𝑢𝑁2,𝜎3.(3.16) Since 𝜎1𝐷𝑁,(𝜎2,𝜎3)𝐷𝑁, 12𝑑𝜎𝑑𝑡12𝐿2+𝛼2𝜎12𝐿2=𝜕𝑦𝜎2,𝜎1𝜕𝑦𝜎2𝐿2𝜎1𝐿2,(3.17) then we have 12𝑑𝜎𝑑𝑡12𝐿2+𝛼2𝜎12𝐿2𝜕𝑦𝜎2𝐿2||𝜎2||.(𝐻,𝑡)(3.18) Consequently, from hypothesis and Gronwall’s inequality, we obtain sup𝑡0,𝑡𝑁𝜎1𝐿2Γ1𝑁1𝑚,(3.19) where Γ1 is a function of 𝑇,𝛽,𝛼2, and𝑐: 𝜎22𝐿2+𝑆112𝑑𝜎𝑑𝑡22𝐿2+12𝜇1𝑆𝑥𝑥𝜕𝑦𝑢𝑆𝑁𝑥𝑥𝜕𝑦𝑢𝑁,𝜎2=𝜕𝑦𝜎1,𝜎2+𝑆2𝜕2𝑡𝑦𝜎1,𝜎2.(3.20) Hence, we get 𝜎22𝐿2+𝑆112𝑑𝜎𝑑𝑡22𝐿2+12𝜇1𝑆𝑥𝑥𝜕𝑦𝑢𝑆𝑁𝑥𝑥𝜕𝑦𝑢𝑁,𝜎2𝜕𝑦𝜎1𝐿2𝜎2𝐿2+𝑆2𝜕2𝑡𝑦𝜎1𝐿2𝜎2𝐿2,𝜎(3.21)32𝐿2+𝑆112𝑑𝜎𝑑𝑡32𝐿22𝑆1𝑆𝑥𝑦𝜕𝑦𝑢𝑆𝑁𝑥𝑦𝜕𝑦𝑢𝑁,𝜎3=2𝑆2𝜕𝑦𝑢2𝜕𝑦𝑢𝑁2,𝜎3.(3.22) Let us estimate third term on the left-hand side of (3.20) in the time interval [0,𝑇𝑁): 𝑆𝑥𝑥𝜕𝑦𝑢𝑆𝑁𝑥𝑥𝜕𝑦𝑢𝑁,𝜎2=𝑆𝑥𝑥𝑆𝑁𝑥𝑥𝜕𝑦𝑢+𝑢𝑁+𝑆𝑁𝑥𝑥𝜕𝑦𝑢𝑢𝑁+𝜕𝑦𝑢𝑁𝑆𝑥𝑥𝑆𝑁𝑥𝑥,𝜎2.(3.23) Thus, ||𝑆𝑥𝑥𝜕𝑦𝑢𝑆𝑁𝑥𝑥𝜕𝑦𝑢𝑁,𝜎2||sup𝑦||𝜕𝑦𝑢+𝑢𝑁||𝑆𝑥𝑥𝑆𝑁𝑥𝑥𝐿2𝜎2𝐿2+sup𝑦||𝑆𝑁𝑥𝑥||𝜕𝑦𝑢𝑢𝑁𝐿2𝜎2𝐿2+sup𝑦||𝜕𝑦𝑢𝑁||𝑆𝑥𝑥𝑆𝑁𝑥𝑥𝐿2𝜎2𝐿2𝑐𝑢+𝑢𝑁𝐻2𝑢+𝑐𝑁𝐻2𝑆𝑥𝑥𝑃𝑁𝑆𝑥𝑥𝐿2+𝑃𝑁𝑆𝑥𝑥𝑆𝑁𝑥𝑥𝐿2𝜎2𝐿2𝑆+𝑐𝑁𝑥𝑥𝐻1𝑢𝑃𝑁𝑢𝐻1+𝑃𝑁𝑢𝑢𝑁𝐻1𝜎2𝐿2×𝐶5𝑐𝜓𝑝𝑁𝑚𝑆𝑥𝑥𝐻𝑚+𝜎2𝐿2𝜎2𝐿2+2𝑐1𝐶𝑝𝑁1𝑚𝑢𝐻𝑚+𝜎2𝐿2𝜎2𝐿2.(3.24) Hence, 𝜎22𝐿2+𝑆1𝑑𝜎𝑑𝑡2𝐿2𝜓2𝜎2𝐿2+5𝑐𝜓𝜅𝑁𝑚+2𝑐𝜅1𝜓1𝐶𝑝𝑁1𝑚+𝜕𝑦𝜎1𝐿2𝜓2𝜎12𝐿2+5𝑐𝜓𝐶𝑝𝜅𝑁𝑚+2𝑐𝜅1𝜓1𝐶𝑝𝑁1𝑚+𝜕𝑦𝜎1𝐿2𝑐1𝜎2𝐿2+𝑐𝑁1𝑚+𝜕𝑦𝜎1𝐿2+𝑆2𝜕2𝑡𝑦𝜎1𝐿2.(3.25) Using the hypothesis and Gronwall’s inequality, sup𝑡0,𝑡𝑁𝜎2𝐿2Γ2𝑁1𝑚,(3.26)Γ1 is a function of 𝑇,𝑆1,𝑆2,𝜅1,𝜓,𝜓1,𝜓2,𝑐1,𝛽1, and𝑐. Estimating the RHS of (3.22) in the time interval [0,𝑇𝑁), 𝜕𝑦𝑢2𝜕𝑦𝑢𝑁2,𝜎3=𝜕𝑦𝑢+𝑢𝑁𝑢𝑢𝑁,𝜎3=𝜕𝑦𝑢+𝑢𝑁𝑢𝑢𝑁,𝜎3+𝑢+𝑢𝑁𝜕𝑦𝑢𝑢𝑁,𝜎3=𝜕𝑦𝑢+𝑢𝑁𝑢𝑢𝑁,𝜎3+𝑢+𝑢𝑁𝜕𝑦𝑢𝑃𝑁𝑢,𝜎3+𝑢+𝑢𝑁𝜕𝑦𝑃𝑁𝑢𝑢𝑁,𝜎3,|||𝜕(3.27)𝑦𝑢2𝜕𝑦𝑢𝑁2,𝜎3|||sup𝑦||𝜕𝑦𝑢+𝑢𝑁||𝑢𝑢𝑁𝐿2𝜎3𝐿2+sup𝑦||𝜕𝑦𝑢+𝑢𝑁||𝜕𝑦𝑢𝑃𝑁𝑢𝐿2𝜎3𝐿2+||||𝐻0𝑢+𝑢𝑁𝜕𝑦𝜎3𝜎3||||𝑑𝑦𝑐𝑢+𝑢𝑁𝐻2𝑢𝑃𝑁𝑢𝐿2+𝑃𝑁𝑢𝑢𝑁𝐿2𝜎3𝐿2+𝑐𝑢+𝑢𝑁𝐻1𝑢𝑃𝑁𝑢𝐻1𝜎3𝐿2+12𝐻0𝜎23||𝜕𝑦𝑢+𝑢𝑁||𝐶𝑑𝑦3𝑐𝜓𝑝𝑁𝑚𝑢𝐻𝑚+𝜎3𝐿2𝜎3𝐿2+3𝑐𝜓𝐶𝑝𝑁1𝑚𝑢𝐻𝑚𝜎3𝐿2+12sup𝑦||𝜕𝑦𝑢+𝑢𝑁||𝐻0𝜎23𝑑𝑦.(3.28) Noting that the last integral is bounded by(1/2)3𝑐𝜓, the estimate is |||𝜕𝑦𝑢2𝜕𝑦𝑢𝑁2,𝜎3|||𝜎3𝑐𝜓3𝐿232𝜎3𝐿2+𝐶𝑝𝑢𝐻𝑚𝑁𝑚+𝑁1𝑚.(3.29)(𝑆𝑥𝑦𝜕𝑦𝑢𝑆𝑁𝑥𝑦𝜕𝑦𝑢𝑁,𝜎3)can be estimated in exactly the same way as (𝑆𝑥𝑥𝜕𝑦𝑢𝑆𝑁𝑥𝑥𝜕𝑦𝑢𝑁,𝜎2). Then, (3.22) as a whole is estimated as 𝜎3𝐿2+𝑆1𝑑𝜎𝑑𝑡3𝐿2𝜓2𝜎3𝐿2+5𝑐𝜓𝐶𝑝𝜅𝑁𝑚+2𝑐𝜅1𝜓1𝐶𝑝𝑁1𝑚3+3𝑐𝜓2𝜎3𝐿2+𝐶𝑝𝜅𝑁𝑚+𝑁1𝑚.(3.30) Therefore, we get 𝑐1𝜎3𝐿2+𝑆1𝑑𝜎𝑑𝑡3𝐿2𝑐𝑁1𝑚.(3.31) Then, using the Gronwall’s inequality, sup𝑡0,𝑡𝑁𝜎3𝐿2Γ3𝑁1𝑚,(3.32) where Γ3 is a function of 𝑇,𝑆1,𝜅1,𝜓,𝜓1,𝑐1, and𝑐. Since 𝑢𝑢𝑁𝐿2=𝑢+𝑃𝑁𝑢𝑃𝑁𝑢𝑢𝑁𝐿2𝑢𝑃𝑁𝑢𝐿2+𝑃𝑁𝑢𝑢𝑁𝐿2,𝑆(3.33)𝑥𝑦𝑆𝑁𝑥𝑦𝐿2=𝑆𝑥𝑦+𝑃𝑁𝑆𝑥𝑦𝑃𝑁𝑆𝑥𝑦𝑆𝑁𝑥𝑦𝐿2𝑆𝑥𝑦𝑃𝑁𝑆𝑥𝑦𝐿2+𝑃𝑁𝑆𝑥𝑦𝑆𝑁𝑥𝑦𝐿2,𝑆(3.34)𝑥𝑥𝑆𝑁𝑥𝑥𝐿2=𝑆𝑥𝑥+𝑃𝑁𝑆𝑥𝑥𝑃𝑁𝑆𝑥𝑥𝑆𝑁𝑥𝑥𝐿2𝑆𝑥𝑥𝑃𝑁𝑆𝑥𝑥𝐿2+𝑃𝑁𝑆𝑥𝑥𝑆𝑁𝑥𝑥𝐿2.(3.35) Using (2.18) and (3.19) in (3.33), (2.18) and (3.26) in (3.34), (2.18) and (3.32) in (3.35), respectively, we obtain sup𝑡0,𝑡𝑁𝑢𝑢𝑁𝐿2Γ1𝑁1𝑚,sup𝑡0,𝑡𝑁𝑆𝑥𝑦𝑆𝑁𝑥𝑦𝐿2Γ2𝑁1𝑚,sup𝑡0,𝑡𝑁𝑆𝑥𝑥𝑆𝑁𝑥𝑥𝐿2Γ3𝑁1𝑚,(3.36) where Γ1,Γ2, and Γ3 are constants, functions of (𝑇,𝛼2,𝐶𝑝,𝜅,𝛽, and 𝑐), (𝑇,𝑆1,𝑆2,𝜅1,𝜓,𝜓1,𝜓2,𝑐1,𝛽1, and 𝑐), and of (𝑇,𝑆1,𝜅1,𝜓,𝜓1,𝑐1,𝜅, and𝑐), respectively.

Lemma 3.3. Suppose that the solution {𝑢𝑁,𝑆𝑁𝑥𝑦,𝑆𝑁𝑥𝑥} of (3.1)–(3.3) exists on the time interval [0,𝑡𝑁] and that sup𝑡[0,𝑡𝑁]𝑢𝑁(𝑦,𝑡)𝐻22𝜓, sup𝑡[0,𝑡𝑁]𝑆𝑁𝑥𝑦(𝑦,𝑡)𝐻22𝜓1, sup𝑡[0,𝑡𝑁]𝑆𝑁𝑥𝑥(𝑦,𝑡)𝐻22𝜓2, sup𝑡[0,𝑡𝑁]|𝑃𝑁𝑢(,𝑡)𝑢𝑁(,𝑡)|𝛽𝑁1𝑚, and sup𝑡[0,𝑡𝑁]|𝑃𝑁𝑆𝑥𝑦(,𝑡)𝑆𝑁𝑥𝑦(,𝑡)|𝛽1𝑁1𝑚, then the error estimate: sup𝑡0,𝑡𝑁𝑢𝑢𝑁𝐻2Γ1𝑁3𝑚,sup𝑡0,𝑡𝑁𝑆𝑥𝑦𝑆𝑁𝑥𝑦𝐻2Γ2𝑁3𝑚,sup𝑡0,𝑡𝑁𝑆𝑥𝑥𝑆𝑁𝑥𝑥𝐻2Γ3𝑁3𝑚,(3.37) holds for the constants Γ1,Γ2, and Γ3. The proof of the Lemma follows from (3.19), (3.26), and (3.32) after application of the triangle inequality and the inverse inequality (2.19).

Proof of Theorem 3.1. To extend the estimate of the first inequality in (3.14) to the time interval[0,𝑇]𝑡𝑁 unspecified in Lemma 3.2 is defined as 𝑡𝑁[]=sup𝑡0,𝑇forall𝑡𝑢𝑡,𝑁𝑦,𝑡𝐻22𝜓.(3.38) Thus, the time 𝑡𝑁 corresponds to the largest time in [0,𝑇] for which the 𝐻2-norm of 𝑢𝑁 is uniformly bounded by 2𝜓. Since 𝑢𝑁(𝑦,0)𝐻2=𝑃𝑁(𝑦,0)𝐻2, 𝑢𝑁(𝑦,0)𝐻2𝑢(𝑦,0)𝐻2𝜓,(3.39) therefore 𝑡𝑁>0 for all 𝑁. Note that 𝑡𝑁 is smaller than the maximum time of existence of the solution 𝑇𝑁. Now, we need to show that there exists 𝑁𝐿 such that 𝑡𝑁=𝑇𝑁𝑁𝐿,(3.40) and therefore the supremum in (3.14) holds on [0,𝑇]. From the definition (3.38), we either have 𝑡𝑁=𝑇 or 𝑡𝑁<𝑇 in which case𝑢𝑁(𝑦,𝑡)𝐻2=2𝜓. Now assume that 𝑡𝑁<𝑇, then 𝑢2𝜓=𝑁𝑦,𝑡𝑁𝐻2𝑢𝑁𝑦,𝑡𝑁𝑢(𝑦,𝑡)𝐻2+𝑢(𝑦,𝑡)𝐻2=𝑢𝑁𝑦,𝑡𝑁𝑢(𝑦,𝑡)𝐻2+𝜓.(3.41) Hence, we obtain 𝑢𝜓𝑁𝑦,𝑡𝑁𝑢(𝑦,𝑡)𝐻2.(3.42) On the other hand, Lemma 3.3 implies 𝜓Γ1𝑁3𝑚(3.43) or Γ𝑁1𝜓1/𝑚3.(3.44) In conclusion, for𝑁𝐿>(Γ1/𝜓)1/𝑚3, we cannot have 𝑡𝑁<𝑇 and claim (3.40) holds. It follows that 𝑁𝑁𝐿 the solution 𝑢𝑁 of (3.1) is defined on[0,𝑇], since as noted before 𝑡𝑁<𝑇𝑁, and, from (3.14), sup[]𝑡0,𝑇𝑢(𝑦,𝑡)𝑢𝑁(𝑦,𝑡)𝐿2Γ1𝑁1𝑚.(3.45)
In exactly the same way, we can extend the estimate of the second and third inequalities in (3.14) to the time interval[0,𝑇] and show that sup[]𝑡0,𝑇𝑆𝑥𝑦(𝑦,𝑡)𝑆𝑁𝑥𝑦(𝑦,𝑡)𝐿2Γ2𝑁1𝑚,sup[]𝑡0,𝑇𝑆𝑥𝑥(𝑦,𝑡)𝑆𝑁𝑥𝑥(𝑦,𝑡)𝐿2Γ3𝑁1𝑚.(3.46)

4. Numerical Results and Discussion

The system of differential equations (3.6)–(3.9) is of the following form𝑑𝑑𝑡̂𝑢𝑁(𝑘,𝑡)=𝐺1̂𝑠𝑁𝑥𝑦,𝑑(𝑘,𝑡)𝑑𝑡̂𝑠𝑁𝑥𝑦(𝑘,𝑡)=𝐺2̂𝑠𝑁𝑥𝑥(𝑖,𝑡),̂𝑠𝑁𝑥𝑦(𝑘,𝑡),̂𝑢𝑁(𝑘,𝑡),𝑆1,𝑆2,𝑑𝑑𝑡̂𝑠𝑁𝑥𝑥(𝑘,𝑡)=𝐺3̂𝑠𝑁𝑥𝑥(𝑖,𝑡),̂𝑠𝑁𝑥𝑦(𝑘,𝑡),̂𝑢𝑁(𝑘,𝑡),𝑆1,𝑆2,̂𝑢𝑁(𝑘,𝑡)=0,̂𝑠𝑁𝑥𝑦(𝑘,𝑡)=0,̂𝑠𝑁𝑥𝑥(𝑘,𝑡)=0.(4.1) Runge-Kutta method is applied to this system. The integrals in equations (2.8) and (2.9) are calculated analytically and numerically, respectively, with 𝜕𝑄/𝜕approximated by a central difference formula. To illustrate the spectral accuracy, the time step is chosen to be sufficiently small so that the error is dominated by the spatial discretization. The free drainage of the Oldroyd-B liquid (𝜇1=0) for which an exact analytical solution is possible is considered first [14]. Figure 1 compares the exact analytical solution in [14] with the approximate solution with 𝑁=5 nodes only for both permeable and impermeable wall. The exact and approximate solutions are indistinguishable in the figure. The error log10(𝑢𝑁𝑢𝐿[0,])at 𝑡=1 of the Fourier-Galerkin approximations with increasing number of nodes for the drainage of Oldroyd-B liquid is listed in Table 1.

Table 1
Figure 1: Comparison of exact analytical and numerical solutions of centerline velocity field for 𝜇1=0, 𝑆1=2, 𝑆2=1, and =1 for both permeable wall (𝛼=1) and impermeable wall.

This shows that numerical results are at least accurate up to the seventh decimal for 𝑁=20. The aim of this paper is to elaborate the effects of the nonlinear parameter 𝜇1 and porous medium parameter 𝛼2 on the centerline velocity and drainage rate. The effect of these parameters on the velocity field is shown in Figures 2 and 3. Figure 2 displays the effect of the nonlinear parameter 𝜇1 on the centerline velocity when the wall is impermeable 𝛼=0, and 𝑆1=2, 𝑆2=1, and =1. Clearly, the overshoot gradually disappears as numerical values of the nonlinear parameter 𝜇1 increase. In addition, the steady centerline velocity increases with increasing values of 𝜇1. Figure 3 explores the effect of porosity on the centerline velocity for 𝜇1=10. Increasing the porosity parameter triggers a decrease in the value of the centerline velocity. The difference between the relaxation and retardation times,𝑆1𝑆2, is a measure of the elasticity of the Oldroyd four-constant liquid, the greater the difference is the more elastic the liquid is. The effect of elasticity on the centerline velocity of constant viscosity Oldroyd four-constant liquids is shown in Figures 4 and 5 for permeable and impermeable walls, respectively. In either case, the centerline velocity increases with an increase in elasticity. The effect of the nonlinear parameter 𝜇1 and porous medium parameter 𝛼2 on the drainage rate is examined in Figures 6 and 7, respectively. Since in all cases the nonlinear parameter 𝜇1 has an increasing effect on the velocity, we expect increasing 𝜇1 will lead to a thinner film over either type of wall, impermeable or permeable. That is evident in Figure 6, which shows that liquid drain more rapidly as𝜇1 is increased from zero. Since the porous medium parameter 𝛼2 has a decreasing effect on the velocity in all cases (both permeable and impermeable wall), we expect increasing 𝛼2 will lead to a thicker film over either impermeable or permeable wall, Figure 7. The effect of the elasticity on the drainage rate is shown in Figure 8 for three liquids of which liquid 1 is the most elastic and liquid 3 is the least elastic. Liquid 1 drains more rapidly than liquid 2, which in turn drains more rapidly than liquid 3.

Figure 2: The effect of the nonlinear parameter 𝜇1 on the centerline velocity profile for 𝑆1=2, 𝑆2=1, 𝛼=0, and =1.
Figure 3: The effect of the porous medium parameter 𝛼2 on the centerline velocity profile for 𝑆1=2, 𝑆2=1, 𝜇1=10, and =1.
Figure 4: The effect of the elasticityon the centerline velocity profile for 𝜇1=10 and =1 for permeable wall (𝛼=1).
Figure 5: The effect of the elasticityon the centerline velocity profile for 𝜇1=10 and =1 for impermeable wall.
Figure 6: The effect of the material parameter 𝜇1 on the drainage rate for 𝑆1=2, 𝑆2=1 for impermeable wall.
Figure 7: The effect of porous medium parameter 𝛼2on the drainage rate for 𝑆1=2, 𝑆2=1, and 𝜇1=10.
Figure 8: The effect of elasticity on the drainage rate for 𝛼=1 and 𝜇1=1.


  1. J. Q. Xu, “Modeling unsteady-state gravity-driven flow in porous media,” Journal of Petroleum Science and Engineering, vol. 62, no. 3-4, pp. 80–86, 2008. View at Publisher · View at Google Scholar · View at Scopus
  2. J. A. Goshawk, N. D. Waters, G. K. Rennie, and E. J. Staples, “Enhancement of the drainage of non-Newtonian liquid films by oscillation,” Journal of Non-Newtonian Fluid Mechanics, vol. 51, no. 1, pp. 21–60, 1994. View at Publisher · View at Google Scholar · View at Scopus
  3. J. A. Goshawk and N. D. Waters, “The effect of oscillation on the drainage of an elastico-viscous liquid,” Journal of Non-Newtonian Fluid Mechanics, vol. 54, pp. 449–464, 1994. View at Publisher · View at Google Scholar · View at Scopus
  4. S. V. Pennington, N. D. Waters, and E. W. Williams, “The numerical simulation of an oldroyd liquid draining down a vertical surface,” Journal of Non-Newtonian Fluid Mechanics, vol. 34, no. 2, pp. 221–239, 1990. View at Publisher · View at Google Scholar · View at Scopus
  5. A. M. Keeley, G. K. Rennie, and N. D. Waters, “Draining thin films-part 1,” Journal of Non-Newtonian Fluid Mechanics, vol. 28, no. 2, pp. 213–226, 1988. View at Publisher · View at Google Scholar · View at Scopus
  6. Y. Dimakopoulos and J. Tsamopoulos, “On the gas-penetration in straight tubes completely filled with a viscoelastic fluid,” Journal of Non-Newtonian Fluid Mechanics, vol. 117, no. 2-3, pp. 117–139, 2004. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  7. M. Pavlidis, Y. Dimakopoulos, and J. Tsamopoulos, “Fully developed flow of a viscoelastic film down a vertical cylindrical or planar wall,” Rheologica Acta, vol. 48, no. 9, pp. 1031–1048, 2009. View at Publisher · View at Google Scholar · View at Scopus
  8. Y. Maday and A. Quarteroni, “Error analysis for spectral approximation of the Korteweg-de Vries equation,” RAIRO Modélisation Mathématique et Analyse Numérique, vol. 22, no. 3, pp. 499–529, 1988. View at Google Scholar · View at Zentralblatt MATH
  9. H. Kalisch and X. Raynaud, “Convergence of a spectral projection of the Camassa-Holm equation,” Numerical Methods for Partial Differential Equations, vol. 22, no. 5, pp. 1197–1215, 2006. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  10. A. Quarteroni and A. Valli, Numerical Approximation of Partial Differential Equations, Springer, Berlin, Germany, 1997.
  11. C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang, Spectral Methods in Fluid Dynamics, Springer, New York, NY, USA, 1988.
  12. F. T. Akyildiz, K. Vajravelu, and H. Ozekes, “Fourier-Galerkin domain truncation method for Stokes' first problem with Oldroyd four-constant liquid,” Computers and Mathematics with Applications, vol. 55, no. 11, pp. 2452–2457, 2008. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  13. N. D. Waters and A. M. Keeley, “Start-up of an elastico-viscous liquid draining from a vertical surface,” Journal of Non-Newtonian Fluid Mechanics, vol. 22, no. 3, pp. 325–334, 1987. View at Publisher · View at Google Scholar · View at Scopus
  14. F. T. Akyildiz and D. A. Siginer, “Start-up of an elastico-viscous liquid draining from a porous vertical surface,” Journal of Non-Newtonian Fluid Mechanics. In press. View at Publisher · View at Google Scholar