Numerical investigation of the generalized ... - Semantic Scholar
Applied Mathematics and Computation 217 (2010) 2631–2638
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Numerical investigation of the generalized lubrication equation E. Momoniat *, C. Harley, E. Adlem Centre for Differential Equations, Continuum Mechanics and Applications, School of Computational and Applied Mathematics, University of the Witwatersrand, Johannesburg, Private Bag 3, Wits 2050, South Africa
a r t i c l e
i n f o
a b s t r a c t Explicit, implicit–explicit and Crank–Nicolson implicit–explicit numerical schemes for solving the generalized lubrication equation are derived. We prove that the implicit– explicit and Crank–Nicolson implicit–explicit numerical schemes are unconditionally stable. Numerical solutions obtained from both schemes are compared. Initial curves with both zero and finite contact angles are considered. Ó 2010 Elsevier Inc. All rights reserved.
Keywords: Thin film Generalized lubrication Crank–Nicolson
1. Introduction In this paper we consider numerical solutions of the generalized lubrication equation [14,20]
! 3 @h @ n@ h ; ¼ h @t @x @x3
h P 0;
ð1:1Þ
where the non-negative function h(x, t) denotes the height of the free surface above the solid substrate, x is the horizontal coordinate and t the time. The generalized lubrication Eq. (1.1) models the spreading of a thin film driven by surface tension. The constant n denotes the kind of flow. The case n = 3 corresponds to the surface tension driven spreading of a thin Newtonian fluid [14]. The case n = 1 models the thickness of a thin film in a Hele-Shaw cell [10]. The fourth-order nonlinear partial differential Eq. (1.1) is degenerate because the coefficient of the highest derivative tends to zero as h tends to zero. Non-negative solutions admitted by (1.1) have been investigated by Beretta et al. [4]. Smyth and Hill [24] investigate similarity solutions of (1.1) which lead to waiting time solutions. A generalisation of this similarity solution has been obtained in Momoniat et al. [19] using a Lie group approach and is given by
jðx þ a1 Þ4 hðx; tÞ ¼ t þ a2
!1=3 ;
ð1:2Þ
where a1 and a2 are constants and j = 81/56. They also investigate asymptotic solutions valid for small 0 < n 1. Further similarity solutions have been investigated by Witelski [28] and Bernoff and Witelski [6]. Investigations into the behaviour of analytical solutions of (1.1) have been undertaken by Bernis et al. [5] and Bowen and King [9]. Oron et al. [22] discuss applications of thin film flow. Myers [20] presents a review of analytical results relating to various investigations of the generalized lubrication equation. Asymptotic solutions to (1.1) for dip coating, draining flows, contact lens flows and spreading drops are discussed. In this paper we derive an explicit, implicit–explicit and a Crank–Nicolson implicit–explicit numerical scheme to determine numerical solutions of (1.1). The degenerate nature of (1.1) makes imposing boundary conditions at the contact line
* Corresponding author. E-mail addresses: [email protected], [email protected] (E. Momoniat). 0096-3003/$ - see front matter Ó 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2010.08.001
2632
E. Momoniat et al. / Applied Mathematics and Computation 217 (2010) 2631–2638
h = 0 very challenging. Bertozzi [7,8] focuses specifically on the contact line singularity and other singularities that can develop in solutions of (1.1). Existence and the behaviour of solutions of (1.1) based on entropy estimates are investigated by Passo et al. [23]. Numerical schemes based on these entropy estimates in higher spatial dimensions have been considered by Grün [15]. Finite element approximations to (1.1) are investigated by Barrett et al. [2] and An et al. [1]. Zhornitskaya and Bertozzi [29] develop finite difference schemes that ensure h P 0. The finite difference schemes satisfy conservation of mass, conservation of surface energy dissipation and a nonlinear energy dissipation. Ha et al. [16] have investigated numerical solutions of (1.1) coupled with a nonlinear source term @/@x(u2u3). Ha et al. [16] solve the nonlinear diffusion term given by (1.1) using a Crank–Nicolson scheme. In this paper we show that an implicit–explicit scheme is unconditionally stable and is more efficient for solving the nonlinear diffusion Eq. (1.1). In this paper we firstly derive an explicit numerical scheme for solving (1.1). We show that this explicit numerical scheme is conditionally stable. We then derive an implicit–explicit scheme. Implicit–explicit schemes for solving hyperbolic partial differential equations have been investigated by McGuire and Morris [17]. Taha and Ablowitz [26] in particular show how these implicit–explicit schemes can be used to obtain numerical solutions of the nonlinear Schrodinger equation. We prove the implicit–explicit scheme is unconditionally stable. We then derive a Crank–Nicolson [11] implicit–explicit scheme. We prove that this Crank–Nicolson implicit–explicit scheme is unconditionally stable. Numerical solutions obtained from the implicit-scheme and the Crank–Nicolson implicit–explicit scheme are applied to initial curves with both zero and finite contact angles. The advantage of the approach taken in this paper is that we obtain numerical schemes that are unconditionally stable. The numerical schemes are evaluated by solving a linear system of equations. This makes the numerical schemes we have derived easy to implement. The accuracy of our numerical solutions depend on the steplengths in the time and spatial directions. The paper is divided up as follows: In Section 2 we derive an explicit, implicit–explicit and a Crank–Nicolson implicit– explicit numerical scheme. We show that the explicit numerical scheme is conditionally stable while the implicit–explicit and Crank–Nicolson implicit–explicit numerical schemes are unconditionally stable. In Section 3 we implement the implicit–explicit and Crank–Nicolson implicit–explicit numerical schemes by solving a linear system of equations. Concluding remarks are made in Section 4. 2. Stability analysis In this section we derive and investigate the stability of numerical schemes to solve (1.1). Derivatives in (1.1) are approximated using finite differences. We approximate the time derivative by a forward difference approximation
@hi;j hi;jþ1 hi;j þ OðDt Þ: @t Dt
ð2:1Þ
A central difference approximation to the first and second derivatives is given by
@hi;j hiþ1;j hi1;j þ O Dx2 ; @x 2Dx @ 2 hi;j hiþ1;j 2hi;j þ hi1;j þ O Dx2 ; @x2 Dx2
ð2:2Þ ð2:3Þ
where hi,j = h(xi, tj). The domain x 2 [a, a] is divided into m + 1 intervals of equal length. We define xi+1 = a + iDx where Dx = 2a/m. We approximate the third derivative as a first derivative of the second derivative to obtain the approximation
@ 3 hi;j hiþ2;j 2hiþ1;j þ 2hi1;j hi2;j þ O D x2 : @x3 2 Dx3
ð2:4Þ
To overcome the contact line singularity an appropriate boundary condition is given by
hða; tÞ ¼ ;
1;
ð2:5Þ
where is the height of the precursor film [18,13,21,25,8]. To close the problem we include the additional boundary conditions
@hða; tÞ ¼ 0: @x
ð2:6Þ
The boundary condition (2.6) implies that the contact angle, the angle made by the film with the substrate, is zero. This implies that complete wetting takes place. We impose the normalizing condition [27,12]
hð0; 0Þ ¼ 1;
ð2:7Þ
i.e. the height of the film at the center is one. An appropriate initial profile that satisfies the boundary conditions (2.5), (2.6) and (2.7) on a large enough finite domain is given by
hðx; 0Þ ¼ exp x2 : If the domain under consideration is given by [a, a] then the height of the precursor film is given by
ð2:8Þ
= exp(a2).
E. Momoniat et al. / Applied Mathematics and Computation 217 (2010) 2631–2638
2633
2.1. Explicit approximation Approximating the outer spatial derivative in (1.1) first we obtain
@hi;j 1 @ 3 hiþ1;j @ 3 hi1;j n n h ¼ hiþ1;j i1;j 2Dx @t @x3 @x3
!
þ O Dx2 :
ð2:9Þ
Note that both the nonlinear term hn and the linear derivative @ 3h/@x3 are approximated explicitly. Substituting (2.1) and (2.9) into (1.1) and dropping the truncation error term we obtain the difference equation
hi;jþ1 ¼ hi;j
i Dt h n n hiþ1;j hiþ3;j 2hiþ2;j þ 2hi;j hi1;j hi1;j hiþ1;j 2hi;j þ 2hi2;j hi3;j : 4 4 Dx
ð2:10Þ
Setting
ai;j ¼ khni;j ; k ¼
Dt 4 D x4
ð2:11Þ
and dropping the truncation errors we write (2.10) as:
hi;jþ1 ¼ hi;j aiþ1;j hiþ3;j 2hiþ2;j þ 2hi;j hi1;j þ ai1;j hiþ1;j 2hi;j þ 2hi2;j hi3;j :
ð2:12Þ
To perform a von Neumann stability analysis we set
hi;j ¼ U j eIxiDx
ð2:13Þ
2
in (2.12) where I = 1 and x a constant. Substituting (2.13) into (2.12) we obtain
U jþ1
¼ 1 baiþ1;j þ cai1;j ;
ð2:14Þ
b ¼ e3IxDx 2e2IxDx þ 2 eIxDx
ð2:15Þ
c ¼ eIxDx 2 þ 2e2IxDx e3IxDx :
ð2:16Þ
Uj where
and
The stability criteria is usually specified in terms of the amplification factor g where
U jþ1 jg j ¼ j : U
ð2:17Þ
Taking absolute values on the left and right hand side of (2.14) we get
U jþ1 jg j ¼ j ¼ 1 baiþ1;j þ cai1;j : U
ð2:18Þ
Imposing the stability criteria
jg j < 1
ð2:19Þ
on (2.18) we obtain the expression
1 baiþ1;j þ cai1;j < 1:
ð2:20Þ
The boundary condition (2.7) at the center of the film places an upper bound on the value of ai1,j 6 k and ai+1,j 6 k where hi1,j = hi+1,j = 1. The assumption that hi1,j = hi+1,j = 1 is a valid one as there is no source term to increase the height of the film. Since surface tension is the only driving force the height of the film will decrease until a steady profile is reached [27]. The inequality (2.20) simplifies to
j1 þ kðc bÞj < 1:
ð2:21Þ
From definitions (2.15) and (2.16) we find that 2
2
b c ¼ 16 sin ðxDxÞ sin ðxDx=2Þ:
ð2:22Þ
The inequality (2.21) simplifies to
2 2 1 16k sin ðxDxÞ sin ðxDx=2Þ < 1:
ð2:23Þ
2634
E. Momoniat et al. / Applied Mathematics and Computation 217 (2010) 2631–2638 2
2
Given the bound j sin ðxDxÞ sin ðxDx=2Þj 6 1 on the trigonometric function we simplify (2.23) to
0 < k < 1=8:
ð2:24Þ
From (2.11) we express the inequality (2.24) as:
0 < Dt < Dx4 =2:
ð2:25Þ
We can thus conclude that the explicit scheme (2.12) is conditionally stable as the stability criteria (2.19) holds provided (2.25) is true. The stability criteria (2.25) for the explicit scheme suggests that it is impractical to use the explicit scheme (2.12) to determine numerical solutions of (1.1) as a large number of iterations are required to obtain solutions for reasonable times. As an example, if we consider Dx = 0.1 then we require Dt < 0.0005 for stability. To obtain solutions at t = 1.0 more than 20000 iterations are required. A spatial step of Dx = 0.1 gives a spatial truncation error of Oð102 Þ which is too big to give reliable solutions. 2.2. Implicit–explicit scheme An implicit–explicit scheme for solving (1.1) can be obtained by writing (2.9) as:
@hi;j 1 @ 3 hiþ1;jþ1 @ 3 hi1;jþ1 n n hi1;j ¼ hiþ1;j 3 2Dx @t @x @x3
!
þ O Dx2 :
ð2:26Þ
Substituting (2.1) and (2.4) into (2.26), re-arranging terms and dropping the truncation error term we obtain
hi;jþ1 þ aiþ1;j hiþ3;jþ1 2hiþ2;jþ1 þ 2hi;jþ1 hi1;jþ1 ai1;j hiþ1;jþ1 2hi;jþ1 þ 2hi2;jþ1 hi3;jþ1 ¼ hi;j :
ð2:27Þ
We perform a linear stability analysis by substituting (2.13) into (2.27) to obtain
U jþ1 1 þ baiþ1;j cai1;j ¼ U j :
ð2:28Þ
We again impose the conditions ai+1,j 6 k and ai1,j 6 k, where hi+1,j = hi1,j = 1. Then from (2.28)
U jþ1 Uj
!
½1 þ kðb cÞ ¼ 1:
ð2:29Þ
The stability criteria (2.19) combined with (2.22) simplifies (2.29) to
1 < 1: 1 þ 16k sin2 ðxDxÞ sin2 ðxDx=2Þ
ð2:30Þ
The condition (2.30) is satisfied for k > 0 and hence the implicit–explicit scheme (2.27) is unconditionally stable as the criteria (2.19) is always satisfied. 2.3. Crank–Nicolson implicit–explicit scheme A Crank–Nicolson [11] improvement on the scheme (2.26) is given by
@hi;j 1 @ 3 hiþ1;jþ1 @ 3 hiþ1;j n þ ¼ hiþ1;j 4Dx @t @x3 @x3
!
n hi1;j
@ 3 hi1;jþ1 @ 3 hi1;j þ @x3 @x3
!!
þ O Dx2 :
ð2:31Þ
Substituting (2.1) and (2.4) into (2.31), re-arranging terms and dropping the truncation error term we obtain
1 1 hi;jþ1 þ aiþ1;j hiþ3;jþ1 2hiþ2;jþ1 þ 2hi;jþ1 hi1;jþ1 ai1;j hiþ1;jþ1 2hi;jþ1 þ 2hi2;jþ1 hi3;jþ1 2 2 1 1 ¼ hi;j aiþ1;j hiþ3;j 2hiþ2;j þ 2hi;j hi1;j þ ai1;j hiþ1;j 2hi;j þ 2hi2;j hi3;j : 2 2
ð2:32Þ
We perform a linear stability analysis by substituting (2.13) into (2.32) to obtain
1 1 U jþ1 1 þ baiþ1;j cai1;j ¼ U j 1 baiþ1;j cai1;j : 2 2
ð2:33Þ
As before we impose the conditions ai+1,j 6 k and ai1,j 6 k where hi+1,j = hi1,j = 1. The stability criteria (2.19) combined with (2.22) simplifies (2.33) to
1 8k sin2 ðxDxÞ sin2 ðxDx=2Þ < 1: 1 þ 8k sin2 ðxDxÞ sin2 ðxDx=2Þ
ð2:34Þ
The condition (2.34) is satisfied for k > 0 and we can therefore conclude that the Crank–Nicolson implicit–explicit scheme (2.32) is unconditionally stable as the criteria (2.19) is always satisfied.
2635
E. Momoniat et al. / Applied Mathematics and Computation 217 (2010) 2631–2638
In the next section we compare numerical solutions to (1.1) obtained from the implicit–explicit scheme (2.27) and the Crank–Nicolson implicit–explicit scheme (2.32). 3. Comparison of numerical schemes In the previous section we have indicated that the explicit scheme is impractical to implement due to the large number of iterations that would be required to obtain accurate results. In this section we show how the implicit–explicit scheme (2.27) and the Crank–Nicolson implicit–explicit scheme (2.32) can be implemented using linear iterative schemes. The implicit–explicit scheme (2.27) is written as the linear system
Ahðjþ1Þ ¼ hðjÞ ;
ð3:1Þ
where 2
1 0 0 0 0 0 0 6 0 2a2 a2 0 0 6 a2 ð1þ2a2 Þ 6 6 A ¼ 6 2a1 a1 a3 ð1þ2a1 þ2a3 Þ a1 2a3 a3 0 6 2a2 a4 ð1þ2a2 þ2a4 Þ a2 2a4 a4 4 a2 0 a3 2a3 a5 ð1þ2a3 þ a5 Þ a3 2a5
0 0 0 0 a5
3
an5 2an5 an3 ð1þ2an5 þ2an4 Þ an5 2an3 an3 0 0 an4 2an4 an2 ð1þ2an4 þ2an2 Þ an4 2an2 an2 7 7 7 0 0 an3 2an3 an1 ð1þ2an3 þ2an1 Þ an1 an3 2an1 7 7 7 0 0 0 an2 2an2 0 ð1þ2an2 Þ an2 5 0
0
0
0
0
0
0
1
ð3:2Þ and
h iT ðjÞ ðjÞ ðjÞ ðjÞ hðjÞ ¼ h0 ; h1 ; . . . ; hn1 ; hn :
ð3:3Þ
We have imposed the boundary conditions (2.5) and (2.6) as: jþ1
h0
j
hn
ðjÞ
h2 ¼ h2 ;
¼ h0 ;
ðjþ1Þ
ðjÞ
¼ hn
ð3:4Þ
and ðjÞ
h1 ¼ h1 ;
ðjÞ
ðjÞ
ðjÞ
ðjÞ
hnþ1 ¼ hn1 ;
ðjÞ
ðjÞ
hnþ2 ¼ hn2 ;
ð3:5Þ
respectively. The Crank–Nicolson implicit–explicit scheme (2.32) is implemented as:
Ahðjþ1Þ ¼ BhðjÞ ;
ð3:6Þ
where A is given by (3.2) and 2
1 0 0 0 0 0 0 0 6 0 2a2 a 2 0 0 0 6 a2 ð1 2a2 Þ 6 B¼6 a1 2a3 a3 0 0 6 2a1 a1 þ a3 ð1 2a1 2a3 Þ 6 2a2 a4 ð1 2a2 2a4 Þ a2 2a4 a4 0 4 a 2 0 a3 2a3 a5 ð1 2a3 a5 Þ a3 2a5 a5
an5 2an5 an3 ð1 2an5 2an4 Þ an5 0 an4 2an4 an2 ð1 2an4 2an2 Þ
0 0 0
0 0 0
an3 0 0
2an3 an2 0
an1 2an2 0
3 0 7 an2 7 7 7: ð1 2an3 2an1 Þ an1 þ an3 2an1 7 7 0 ð1 2an2 Þ an2 5 0 0 1 2an3
an4
an3 2an2
ð3:7Þ We use MATHEMATICA to solve the systems (3.1) and (3.6). In the case of (3.1) the system is evaluated as:
hðjþ1Þ ¼ A1 hðjÞ ;
ð3:8Þ
while for (3.6) the system is evaluated as:
hðjþ1Þ ¼ A1 BhðjÞ :
ð3:9Þ
Both the implicit–explicit scheme and the Crank–Nicolson implicit–explicit scheme converge to the same solution. To test which approach is the most efficient we consider steady solutions to (1.1). For steady solutions we have hi,j+1 = hi,j. The implicit–explicit scheme (3.1) is then evaluated as:
zþ1 ðjÞ z A hðjÞ ¼ h ;
ð3:10Þ
where z is the iteration number and ðhðjÞ Þ0 is an initial guess to the steady profile. Similarly, the Crank–Nicolson implicit– explicit scheme (3.6) is written as:
zþ1 z A hðjþ1Þ ¼ B hðjÞ :
ð3:11Þ
Note here that both (3.10) and (3.11) need to be iterated until the convergence criteria
zþ1 z max hðjÞ hðjÞ < d;
ð3:12Þ
where d is the tolerance is satisfied. We use the zero-contact exponential curve (2.8) as our initial guess. Choosing d = 0.01 we tabulate the number of iterations to convergence to the tolerance d in Table 1. From Table 1 we note that the Crank–Nicolson scheme is more costly as it takes longer to converge.
2636
E. Momoniat et al. / Applied Mathematics and Computation 217 (2010) 2631–2638 Table 1 Table comparing number of iterations required to reach a tolerance of d = 0.01. We have chosen m = 500, Dt = 0.1, a = 4 and Dx = 2a/m = 0.016. The initial curve is given by h(x, 0) = exp(x2). n
Implicit–explicit
Crank–Nicolson implicit–explicit
2 3
10 8
14 12
Fig. 1. Numerical solutions admitted by (1.1) obtained from the implicit–explicit scheme (3.1) for n = 1, n = 2 and n = 3 at t = 0.5. We have chosen m = 500, Dt = 0.1, a = 4 and Dx = 2a/m = 0.016. The initial curve is given by h(x, 0) = exp(x2).
As both the implicit–explicit scheme and the Crank–Nicolson implicit–explicit scheme converge to the same solution it is not necessary to plot the output of both schemes. In Fig. 1 we plot the numerical solution obtained from the implicit–explicit scheme (3.1) for n = 1, n = 2 and n = 3. We note from Fig. 1 that as n increases the film height increases at a fixed time. We next consider initial profiles with finite contact angle. For the arbitrary initial profile
hðx; 0Þ ¼ f ðxÞ
ð3:13Þ
the appropriate boundary conditions are given by [20]
@hða; tÞ df ðaÞ ¼ ; @x dx
3
@ 3 hða; t Þ d f ðaÞ ¼ : 3 @x3 dx
ð3:14Þ
The precursor film height is imposed as:
hða; tÞ ¼ f ðaÞ:
ð3:15Þ
Choosing a = 1 we have
f ð1Þ ¼ ;
df ð1Þ ¼ tan /; dx
f ð0Þ ¼ 1;
ð3:16Þ
where / is the contact angle. A fourth-order polynomial that satisfies (3.16) is given by
1 f ðxÞ ¼ 1 þ 2ð 1Þx2 ð 1Þx4 x2 x2 1 tan /: 2
ð3:17Þ
Imposing the boundary conditions (3.14) as:
hð1; tÞ ¼ ;
@hð1; tÞ df ð1Þ ¼ ; @x dx
3
@ 3 hð1; tÞ d f ð1Þ ¼ ; 3 @x3 dx
ð3:18Þ
we end up solving the system
Ahðjþ1Þ ¼ ChðjÞ ;
ð3:19Þ
E. Momoniat et al. / Applied Mathematics and Computation 217 (2010) 2631–2638
Fig. 2. Numerical solutions admitted by (1.1) obtained from the implicit–explicit scheme (3.19) for n = 3 at t = 0.1. We have chosen m = 250, Dt = 0.1, a = 1 and Dx = 2a/m = 0.008. The initial curve is given by hðx; 0Þ ¼ 1 þ 2ð 1Þx2 ð 1Þx4 12 x2 ðx2 1Þ tan /.
2637
= 0.01,
where
2
1 6 6 1 6 6 1 6 6 6 0 6 6 0 6 6 6 0 6 6 . C ¼ 6 .. 6 6 0 6 6 6 0 6 6 0 6 6 6 0 6 6 4 0 0
0 0
0
0
0 0
0
0
1 0 0 0 0 0 2 0 2 1 0 0
0 0
0 0
0 0
1
0
0 0
0
0
0 0
0
1 0 0
0
0
0 0 .. .. . .
0 .. .
0 1 0 .. .. . . . ..
0 .. .
0 .. .
0 0
0
0
0 1
0
0
0 0
0
0
0 0 1 0
0 0
0
0
0 0
0 1
0 0 0 0
0 0
0 0
0 0 0 0
0 1 0 0
0 0
0
0
0 0
0
0
0 0
0
0 0 0 0
0 0
0
3
7 07 7 07 7 7 0 0 0 07 7 0 0 0 07 7 7 0 0 0 07 7 .. .. .. 7 7: . . . 7 0 0 0 07 7 7 0 0 0 07 7 0 0 0 07 7 7 2 0 2 1 7 7 7 0 0 1 1 5 0 0 0 1
ð3:20Þ
The system (3.19) is the implicit–explicit scheme for finite contact angle /. In Fig. 2 we plot the numerical solution obtained from (3.19) for different values of /. Concluding remarks are made in the next section. 4. Concluding remarks In this paper we have derived unconditionally stable numerical schemes for solving the lubrication Eq. (1.1). The numerical schemes are based on implicit–explicit and Crank–Nicolson implicit–explicit finite difference approximations. The numerical schemes have a truncation error OðDt þ Dx2 Þ. Both numerical schemes are easily implemented using linear iteration. The correctness of the results obtained in this paper are verified by the fact that both the implicit–explicit and Crank–
2638
E. Momoniat et al. / Applied Mathematics and Computation 217 (2010) 2631–2638
Nicolson implicit–explicit finite difference schemes give the same solutions. We have shown that the implicit–explicit scheme is more efficient than the Crank–Nicolson scheme. This is an improvement on the work of Ha et al. [16]. The finite difference schemes presented in this work are easier to implement than the schemes derived by Zhornitskaya and Bertozzi [29]. The accuracy of the numerical schemes can be improved by considering a higher order approximation to the third-order derivative given by
@ 3 hi;j hiþ3;j þ 8hiþ2;j 13hiþ1;j þ 13hi1;j 8hi2;j þ hi3;j þ O Dx4 : @x3 8Dx3
ð4:21Þ
This improves the truncation error to OðDt þ Dx4 Þ. The resulting matrices are however more complicated and would require more computation to solve. Another possible approach that reduces a nonlinear partial differential equation to a linear system is the quasi-linearisation approach introduced by Bellman and Kalaba [3]. This approach does however need iterations to take place at each time step. References [1] R. An, Y. Li, K. Li, Finite element approximation for fourth-order nonlinear problem in the plane, Appl. Math. Comput. 194 (2007) 143–155. [2] W. Barrett, J.F. Blowey, H. Garcke, Finite element approximation of a fourth order nonlinear degenerate parabolic equation, Numer. Math. 80 (1998) 525–556. [3] R.E. Bellman, R.E. Kalaba, Quasilinearization and Non-Linear Boundary Value Problems, American Elsevier, New York, 1965. [4] E. Beretta, M. Bertsch, R. Dal Passo, Nonnegative solutions of a fourth-order nonlinear degenerate parabolic equation, Arch. Ration. Mech. Anal. 129 (1995) 175–200. [5] F. Bernis, J. Hulshof, J.R. King, Dipoles and similarity solutions of the thin film equation in the half-line, Nonlinearity 13 (2000) 413–439. [6] A.J. Bernoff, T.P. Witelski, Linear stability of source-type similarity solutions of the thin film equation, Appl. Math. Lett. 15 (2002) 599–606. [7] A.L. Bertozzi, Symmetric singularity formation in lubrication-type equations for interface motion, SIAM J. Appl. Math. 56 (1996) 681–714. [8] A.L. Bertozzi, The mathematics of moving contact lines in thin liquid films, Not. Am. Math. Soc. 45 (1998) 689–697. [9] M. Bowen, J.R. King, Asymptotic behaviour of the thin film equation in bounded domains, Eur. J. Appl. Math. 12 (2001) 135–157. [10] P. Constantin, T.F. Dupont, R.E. Goldstein, L.P. Kadanoff, M.J. Shelley, S.-M. Zhou, Droplet breakup in a model of the Hele-Shaw cell, Phys. Rev. E 47 (1993) 4169–4181. [11] J. Crank, E. Nicolson, A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type, P. Camb. Philos. Soc. 43 (1947) 50–67. [12] B.R. Duffy, S.K. Wilson, A third-order differential equation arising in thin-film flows and relevant to Tanner’s Law, Appl. Math. Lett. 10 (1997) 63–68. [13] E.B. Dussan, On the spreading of liquids on solid surfaces: static and dynamic contact lines, Annu. Rev. Fluid Mech. 11 (1979) 371–400. [14] H.P. Greenspan, On the motion of a small viscous droplet that wets a surface, J. Fluid Mech. 84 (1978) 125–143. [15] G. Grün, On the convergence of entropy consistent schemes for lubrication type equations in multiple space dimensions, Math. Comput. 72 (2003) 1251–1279. [16] Y. Ha, Y.-J. Kim, T.G. Myers, On the numerical solution of a driven thin film equation, J. Comput. Phys. 227 (2008) 7246–7263. [17] G.R. McGuire, J.Ll. Morris, Explicit–implicit schemes for the numerical solution of nonlinear hyperbolic systems, Math. Comput. 29 (1975) 407–424. [18] S. Middleman, The effect of induced air flow on the spin coating of viscous liquids, J. Appl. Phys. 62 (1987) 2530–2532. [19] E. Momoniat, T.G. Myers, S. Abelman, New solutions for surface tension driven spreading of a thin film, Int. J. Nonlinear Mech. 40 (2005) 523–529. [20] T.G. Myers, Thin films with high surface tension, SIAM Rev. 40 (1998) 441–462. [21] J.A. Moriarty, L.W. Schwartz, E.O. Tuck, Unsteady spreading of thin liquid films with small surface tension, Phys. Fluids A3 (5) (1991) 733–742. [22] A. Oron, S.H. Davis, S.G. Bankoff, Long-scale evolution of thin liquid films, Rev. Mod. Phys. 69 (1997) 931–980. [23] R.D. Passo, H. Garcke, G. Grün, On a fourth-order degenerate parabolic equation: global entropy estimates, existence, and qualitative behavior of solutions, SIAM J. Math. Anal. 29 (1998) 321–342. [24] N.F. Smyth, J.M. Hill, High-order nonlinear diffusion, IMA J. Appl. Math. 40 (1988) 73–86. [25] M.A. Spaid, G.M. Homsy, Stability of Newtonian and viscoelastic dynamic contact lines, Phys. Fluids 8 (1996) 460–478. [26] T.R. Taha, M.J. Ablowitz, Analytical and numerical aspects of certain nonlinear evolution equations: II. Numerical, nonlinear Schrodinger equation, J. Comput. Phys. 55 (1984) 203–230. [27] L.H. Tanner, The spreading of silicone oil drops on horizontal surfaces, J. Phys. D. Appl. Phys. 12 (1979) 1473–1484. [28] T.P. Witelski, Similarity solutions of the lubrication equation, Appl. Math. Lett. 10 (1997) 107–113. [29] L. Zhornitskaya, A.L. Bertozzi, Positivity-preserving numerical schemes for lubrication-type equations, SIAM J. Numer. Anal. 37 (1999) 523–555.
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Numerical investigation of the generalized lubrication equation E. Momoniat *, C. Harley, E. Adlem Centre for Differential Equations, Continuum Mechanics and Applications, School of Computational and Applied Mathematics, University of the Witwatersrand, Johannesburg, Private Bag 3, Wits 2050, South Africa
a r t i c l e
i n f o
a b s t r a c t Explicit, implicit–explicit and Crank–Nicolson implicit–explicit numerical schemes for solving the generalized lubrication equation are derived. We prove that the implicit– explicit and Crank–Nicolson implicit–explicit numerical schemes are unconditionally stable. Numerical solutions obtained from both schemes are compared. Initial curves with both zero and finite contact angles are considered. Ó 2010 Elsevier Inc. All rights reserved.
Keywords: Thin film Generalized lubrication Crank–Nicolson
1. Introduction In this paper we consider numerical solutions of the generalized lubrication equation [14,20]
! 3 @h @ n@ h ; ¼ h @t @x @x3
h P 0;
ð1:1Þ
where the non-negative function h(x, t) denotes the height of the free surface above the solid substrate, x is the horizontal coordinate and t the time. The generalized lubrication Eq. (1.1) models the spreading of a thin film driven by surface tension. The constant n denotes the kind of flow. The case n = 3 corresponds to the surface tension driven spreading of a thin Newtonian fluid [14]. The case n = 1 models the thickness of a thin film in a Hele-Shaw cell [10]. The fourth-order nonlinear partial differential Eq. (1.1) is degenerate because the coefficient of the highest derivative tends to zero as h tends to zero. Non-negative solutions admitted by (1.1) have been investigated by Beretta et al. [4]. Smyth and Hill [24] investigate similarity solutions of (1.1) which lead to waiting time solutions. A generalisation of this similarity solution has been obtained in Momoniat et al. [19] using a Lie group approach and is given by
jðx þ a1 Þ4 hðx; tÞ ¼ t þ a2
!1=3 ;
ð1:2Þ
where a1 and a2 are constants and j = 81/56. They also investigate asymptotic solutions valid for small 0 < n 1. Further similarity solutions have been investigated by Witelski [28] and Bernoff and Witelski [6]. Investigations into the behaviour of analytical solutions of (1.1) have been undertaken by Bernis et al. [5] and Bowen and King [9]. Oron et al. [22] discuss applications of thin film flow. Myers [20] presents a review of analytical results relating to various investigations of the generalized lubrication equation. Asymptotic solutions to (1.1) for dip coating, draining flows, contact lens flows and spreading drops are discussed. In this paper we derive an explicit, implicit–explicit and a Crank–Nicolson implicit–explicit numerical scheme to determine numerical solutions of (1.1). The degenerate nature of (1.1) makes imposing boundary conditions at the contact line
* Corresponding author. E-mail addresses: [email protected], [email protected] (E. Momoniat). 0096-3003/$ - see front matter Ó 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2010.08.001
2632
E. Momoniat et al. / Applied Mathematics and Computation 217 (2010) 2631–2638
h = 0 very challenging. Bertozzi [7,8] focuses specifically on the contact line singularity and other singularities that can develop in solutions of (1.1). Existence and the behaviour of solutions of (1.1) based on entropy estimates are investigated by Passo et al. [23]. Numerical schemes based on these entropy estimates in higher spatial dimensions have been considered by Grün [15]. Finite element approximations to (1.1) are investigated by Barrett et al. [2] and An et al. [1]. Zhornitskaya and Bertozzi [29] develop finite difference schemes that ensure h P 0. The finite difference schemes satisfy conservation of mass, conservation of surface energy dissipation and a nonlinear energy dissipation. Ha et al. [16] have investigated numerical solutions of (1.1) coupled with a nonlinear source term @/@x(u2u3). Ha et al. [16] solve the nonlinear diffusion term given by (1.1) using a Crank–Nicolson scheme. In this paper we show that an implicit–explicit scheme is unconditionally stable and is more efficient for solving the nonlinear diffusion Eq. (1.1). In this paper we firstly derive an explicit numerical scheme for solving (1.1). We show that this explicit numerical scheme is conditionally stable. We then derive an implicit–explicit scheme. Implicit–explicit schemes for solving hyperbolic partial differential equations have been investigated by McGuire and Morris [17]. Taha and Ablowitz [26] in particular show how these implicit–explicit schemes can be used to obtain numerical solutions of the nonlinear Schrodinger equation. We prove the implicit–explicit scheme is unconditionally stable. We then derive a Crank–Nicolson [11] implicit–explicit scheme. We prove that this Crank–Nicolson implicit–explicit scheme is unconditionally stable. Numerical solutions obtained from the implicit-scheme and the Crank–Nicolson implicit–explicit scheme are applied to initial curves with both zero and finite contact angles. The advantage of the approach taken in this paper is that we obtain numerical schemes that are unconditionally stable. The numerical schemes are evaluated by solving a linear system of equations. This makes the numerical schemes we have derived easy to implement. The accuracy of our numerical solutions depend on the steplengths in the time and spatial directions. The paper is divided up as follows: In Section 2 we derive an explicit, implicit–explicit and a Crank–Nicolson implicit– explicit numerical scheme. We show that the explicit numerical scheme is conditionally stable while the implicit–explicit and Crank–Nicolson implicit–explicit numerical schemes are unconditionally stable. In Section 3 we implement the implicit–explicit and Crank–Nicolson implicit–explicit numerical schemes by solving a linear system of equations. Concluding remarks are made in Section 4. 2. Stability analysis In this section we derive and investigate the stability of numerical schemes to solve (1.1). Derivatives in (1.1) are approximated using finite differences. We approximate the time derivative by a forward difference approximation
@hi;j hi;jþ1 hi;j þ OðDt Þ: @t Dt
ð2:1Þ
A central difference approximation to the first and second derivatives is given by
@hi;j hiþ1;j hi1;j þ O Dx2 ; @x 2Dx @ 2 hi;j hiþ1;j 2hi;j þ hi1;j þ O Dx2 ; @x2 Dx2
ð2:2Þ ð2:3Þ
where hi,j = h(xi, tj). The domain x 2 [a, a] is divided into m + 1 intervals of equal length. We define xi+1 = a + iDx where Dx = 2a/m. We approximate the third derivative as a first derivative of the second derivative to obtain the approximation
@ 3 hi;j hiþ2;j 2hiþ1;j þ 2hi1;j hi2;j þ O D x2 : @x3 2 Dx3
ð2:4Þ
To overcome the contact line singularity an appropriate boundary condition is given by
hða; tÞ ¼ ;
1;
ð2:5Þ
where is the height of the precursor film [18,13,21,25,8]. To close the problem we include the additional boundary conditions
@hða; tÞ ¼ 0: @x
ð2:6Þ
The boundary condition (2.6) implies that the contact angle, the angle made by the film with the substrate, is zero. This implies that complete wetting takes place. We impose the normalizing condition [27,12]
hð0; 0Þ ¼ 1;
ð2:7Þ
i.e. the height of the film at the center is one. An appropriate initial profile that satisfies the boundary conditions (2.5), (2.6) and (2.7) on a large enough finite domain is given by
hðx; 0Þ ¼ exp x2 : If the domain under consideration is given by [a, a] then the height of the precursor film is given by
ð2:8Þ
= exp(a2).
E. Momoniat et al. / Applied Mathematics and Computation 217 (2010) 2631–2638
2633
2.1. Explicit approximation Approximating the outer spatial derivative in (1.1) first we obtain
@hi;j 1 @ 3 hiþ1;j @ 3 hi1;j n n h ¼ hiþ1;j i1;j 2Dx @t @x3 @x3
!
þ O Dx2 :
ð2:9Þ
Note that both the nonlinear term hn and the linear derivative @ 3h/@x3 are approximated explicitly. Substituting (2.1) and (2.9) into (1.1) and dropping the truncation error term we obtain the difference equation
hi;jþ1 ¼ hi;j
i Dt h n n hiþ1;j hiþ3;j 2hiþ2;j þ 2hi;j hi1;j hi1;j hiþ1;j 2hi;j þ 2hi2;j hi3;j : 4 4 Dx
ð2:10Þ
Setting
ai;j ¼ khni;j ; k ¼
Dt 4 D x4
ð2:11Þ
and dropping the truncation errors we write (2.10) as:
hi;jþ1 ¼ hi;j aiþ1;j hiþ3;j 2hiþ2;j þ 2hi;j hi1;j þ ai1;j hiþ1;j 2hi;j þ 2hi2;j hi3;j :
ð2:12Þ
To perform a von Neumann stability analysis we set
hi;j ¼ U j eIxiDx
ð2:13Þ
2
in (2.12) where I = 1 and x a constant. Substituting (2.13) into (2.12) we obtain
U jþ1
¼ 1 baiþ1;j þ cai1;j ;
ð2:14Þ
b ¼ e3IxDx 2e2IxDx þ 2 eIxDx
ð2:15Þ
c ¼ eIxDx 2 þ 2e2IxDx e3IxDx :
ð2:16Þ
Uj where
and
The stability criteria is usually specified in terms of the amplification factor g where
U jþ1 jg j ¼ j : U
ð2:17Þ
Taking absolute values on the left and right hand side of (2.14) we get
U jþ1 jg j ¼ j ¼ 1 baiþ1;j þ cai1;j : U
ð2:18Þ
Imposing the stability criteria
jg j < 1
ð2:19Þ
on (2.18) we obtain the expression
1 baiþ1;j þ cai1;j < 1:
ð2:20Þ
The boundary condition (2.7) at the center of the film places an upper bound on the value of ai1,j 6 k and ai+1,j 6 k where hi1,j = hi+1,j = 1. The assumption that hi1,j = hi+1,j = 1 is a valid one as there is no source term to increase the height of the film. Since surface tension is the only driving force the height of the film will decrease until a steady profile is reached [27]. The inequality (2.20) simplifies to
j1 þ kðc bÞj < 1:
ð2:21Þ
From definitions (2.15) and (2.16) we find that 2
2
b c ¼ 16 sin ðxDxÞ sin ðxDx=2Þ:
ð2:22Þ
The inequality (2.21) simplifies to
2 2 1 16k sin ðxDxÞ sin ðxDx=2Þ < 1:
ð2:23Þ
2634
E. Momoniat et al. / Applied Mathematics and Computation 217 (2010) 2631–2638 2
2
Given the bound j sin ðxDxÞ sin ðxDx=2Þj 6 1 on the trigonometric function we simplify (2.23) to
0 < k < 1=8:
ð2:24Þ
From (2.11) we express the inequality (2.24) as:
0 < Dt < Dx4 =2:
ð2:25Þ
We can thus conclude that the explicit scheme (2.12) is conditionally stable as the stability criteria (2.19) holds provided (2.25) is true. The stability criteria (2.25) for the explicit scheme suggests that it is impractical to use the explicit scheme (2.12) to determine numerical solutions of (1.1) as a large number of iterations are required to obtain solutions for reasonable times. As an example, if we consider Dx = 0.1 then we require Dt < 0.0005 for stability. To obtain solutions at t = 1.0 more than 20000 iterations are required. A spatial step of Dx = 0.1 gives a spatial truncation error of Oð102 Þ which is too big to give reliable solutions. 2.2. Implicit–explicit scheme An implicit–explicit scheme for solving (1.1) can be obtained by writing (2.9) as:
@hi;j 1 @ 3 hiþ1;jþ1 @ 3 hi1;jþ1 n n hi1;j ¼ hiþ1;j 3 2Dx @t @x @x3
!
þ O Dx2 :
ð2:26Þ
Substituting (2.1) and (2.4) into (2.26), re-arranging terms and dropping the truncation error term we obtain
hi;jþ1 þ aiþ1;j hiþ3;jþ1 2hiþ2;jþ1 þ 2hi;jþ1 hi1;jþ1 ai1;j hiþ1;jþ1 2hi;jþ1 þ 2hi2;jþ1 hi3;jþ1 ¼ hi;j :
ð2:27Þ
We perform a linear stability analysis by substituting (2.13) into (2.27) to obtain
U jþ1 1 þ baiþ1;j cai1;j ¼ U j :
ð2:28Þ
We again impose the conditions ai+1,j 6 k and ai1,j 6 k, where hi+1,j = hi1,j = 1. Then from (2.28)
U jþ1 Uj
!
½1 þ kðb cÞ ¼ 1:
ð2:29Þ
The stability criteria (2.19) combined with (2.22) simplifies (2.29) to
1 < 1: 1 þ 16k sin2 ðxDxÞ sin2 ðxDx=2Þ
ð2:30Þ
The condition (2.30) is satisfied for k > 0 and hence the implicit–explicit scheme (2.27) is unconditionally stable as the criteria (2.19) is always satisfied. 2.3. Crank–Nicolson implicit–explicit scheme A Crank–Nicolson [11] improvement on the scheme (2.26) is given by
@hi;j 1 @ 3 hiþ1;jþ1 @ 3 hiþ1;j n þ ¼ hiþ1;j 4Dx @t @x3 @x3
!
n hi1;j
@ 3 hi1;jþ1 @ 3 hi1;j þ @x3 @x3
!!
þ O Dx2 :
ð2:31Þ
Substituting (2.1) and (2.4) into (2.31), re-arranging terms and dropping the truncation error term we obtain
1 1 hi;jþ1 þ aiþ1;j hiþ3;jþ1 2hiþ2;jþ1 þ 2hi;jþ1 hi1;jþ1 ai1;j hiþ1;jþ1 2hi;jþ1 þ 2hi2;jþ1 hi3;jþ1 2 2 1 1 ¼ hi;j aiþ1;j hiþ3;j 2hiþ2;j þ 2hi;j hi1;j þ ai1;j hiþ1;j 2hi;j þ 2hi2;j hi3;j : 2 2
ð2:32Þ
We perform a linear stability analysis by substituting (2.13) into (2.32) to obtain
1 1 U jþ1 1 þ baiþ1;j cai1;j ¼ U j 1 baiþ1;j cai1;j : 2 2
ð2:33Þ
As before we impose the conditions ai+1,j 6 k and ai1,j 6 k where hi+1,j = hi1,j = 1. The stability criteria (2.19) combined with (2.22) simplifies (2.33) to
1 8k sin2 ðxDxÞ sin2 ðxDx=2Þ < 1: 1 þ 8k sin2 ðxDxÞ sin2 ðxDx=2Þ
ð2:34Þ
The condition (2.34) is satisfied for k > 0 and we can therefore conclude that the Crank–Nicolson implicit–explicit scheme (2.32) is unconditionally stable as the criteria (2.19) is always satisfied.
2635
E. Momoniat et al. / Applied Mathematics and Computation 217 (2010) 2631–2638
In the next section we compare numerical solutions to (1.1) obtained from the implicit–explicit scheme (2.27) and the Crank–Nicolson implicit–explicit scheme (2.32). 3. Comparison of numerical schemes In the previous section we have indicated that the explicit scheme is impractical to implement due to the large number of iterations that would be required to obtain accurate results. In this section we show how the implicit–explicit scheme (2.27) and the Crank–Nicolson implicit–explicit scheme (2.32) can be implemented using linear iterative schemes. The implicit–explicit scheme (2.27) is written as the linear system
Ahðjþ1Þ ¼ hðjÞ ;
ð3:1Þ
where 2
1 0 0 0 0 0 0 6 0 2a2 a2 0 0 6 a2 ð1þ2a2 Þ 6 6 A ¼ 6 2a1 a1 a3 ð1þ2a1 þ2a3 Þ a1 2a3 a3 0 6 2a2 a4 ð1þ2a2 þ2a4 Þ a2 2a4 a4 4 a2 0 a3 2a3 a5 ð1þ2a3 þ a5 Þ a3 2a5
0 0 0 0 a5
3
an5 2an5 an3 ð1þ2an5 þ2an4 Þ an5 2an3 an3 0 0 an4 2an4 an2 ð1þ2an4 þ2an2 Þ an4 2an2 an2 7 7 7 0 0 an3 2an3 an1 ð1þ2an3 þ2an1 Þ an1 an3 2an1 7 7 7 0 0 0 an2 2an2 0 ð1þ2an2 Þ an2 5 0
0
0
0
0
0
0
1
ð3:2Þ and
h iT ðjÞ ðjÞ ðjÞ ðjÞ hðjÞ ¼ h0 ; h1 ; . . . ; hn1 ; hn :
ð3:3Þ
We have imposed the boundary conditions (2.5) and (2.6) as: jþ1
h0
j
hn
ðjÞ
h2 ¼ h2 ;
¼ h0 ;
ðjþ1Þ
ðjÞ
¼ hn
ð3:4Þ
and ðjÞ
h1 ¼ h1 ;
ðjÞ
ðjÞ
ðjÞ
ðjÞ
hnþ1 ¼ hn1 ;
ðjÞ
ðjÞ
hnþ2 ¼ hn2 ;
ð3:5Þ
respectively. The Crank–Nicolson implicit–explicit scheme (2.32) is implemented as:
Ahðjþ1Þ ¼ BhðjÞ ;
ð3:6Þ
where A is given by (3.2) and 2
1 0 0 0 0 0 0 0 6 0 2a2 a 2 0 0 0 6 a2 ð1 2a2 Þ 6 B¼6 a1 2a3 a3 0 0 6 2a1 a1 þ a3 ð1 2a1 2a3 Þ 6 2a2 a4 ð1 2a2 2a4 Þ a2 2a4 a4 0 4 a 2 0 a3 2a3 a5 ð1 2a3 a5 Þ a3 2a5 a5
an5 2an5 an3 ð1 2an5 2an4 Þ an5 0 an4 2an4 an2 ð1 2an4 2an2 Þ
0 0 0
0 0 0
an3 0 0
2an3 an2 0
an1 2an2 0
3 0 7 an2 7 7 7: ð1 2an3 2an1 Þ an1 þ an3 2an1 7 7 0 ð1 2an2 Þ an2 5 0 0 1 2an3
an4
an3 2an2
ð3:7Þ We use MATHEMATICA to solve the systems (3.1) and (3.6). In the case of (3.1) the system is evaluated as:
hðjþ1Þ ¼ A1 hðjÞ ;
ð3:8Þ
while for (3.6) the system is evaluated as:
hðjþ1Þ ¼ A1 BhðjÞ :
ð3:9Þ
Both the implicit–explicit scheme and the Crank–Nicolson implicit–explicit scheme converge to the same solution. To test which approach is the most efficient we consider steady solutions to (1.1). For steady solutions we have hi,j+1 = hi,j. The implicit–explicit scheme (3.1) is then evaluated as:
zþ1 ðjÞ z A hðjÞ ¼ h ;
ð3:10Þ
where z is the iteration number and ðhðjÞ Þ0 is an initial guess to the steady profile. Similarly, the Crank–Nicolson implicit– explicit scheme (3.6) is written as:
zþ1 z A hðjþ1Þ ¼ B hðjÞ :
ð3:11Þ
Note here that both (3.10) and (3.11) need to be iterated until the convergence criteria
zþ1 z max hðjÞ hðjÞ < d;
ð3:12Þ
where d is the tolerance is satisfied. We use the zero-contact exponential curve (2.8) as our initial guess. Choosing d = 0.01 we tabulate the number of iterations to convergence to the tolerance d in Table 1. From Table 1 we note that the Crank–Nicolson scheme is more costly as it takes longer to converge.
2636
E. Momoniat et al. / Applied Mathematics and Computation 217 (2010) 2631–2638 Table 1 Table comparing number of iterations required to reach a tolerance of d = 0.01. We have chosen m = 500, Dt = 0.1, a = 4 and Dx = 2a/m = 0.016. The initial curve is given by h(x, 0) = exp(x2). n
Implicit–explicit
Crank–Nicolson implicit–explicit
2 3
10 8
14 12
Fig. 1. Numerical solutions admitted by (1.1) obtained from the implicit–explicit scheme (3.1) for n = 1, n = 2 and n = 3 at t = 0.5. We have chosen m = 500, Dt = 0.1, a = 4 and Dx = 2a/m = 0.016. The initial curve is given by h(x, 0) = exp(x2).
As both the implicit–explicit scheme and the Crank–Nicolson implicit–explicit scheme converge to the same solution it is not necessary to plot the output of both schemes. In Fig. 1 we plot the numerical solution obtained from the implicit–explicit scheme (3.1) for n = 1, n = 2 and n = 3. We note from Fig. 1 that as n increases the film height increases at a fixed time. We next consider initial profiles with finite contact angle. For the arbitrary initial profile
hðx; 0Þ ¼ f ðxÞ
ð3:13Þ
the appropriate boundary conditions are given by [20]
@hða; tÞ df ðaÞ ¼ ; @x dx
3
@ 3 hða; t Þ d f ðaÞ ¼ : 3 @x3 dx
ð3:14Þ
The precursor film height is imposed as:
hða; tÞ ¼ f ðaÞ:
ð3:15Þ
Choosing a = 1 we have
f ð1Þ ¼ ;
df ð1Þ ¼ tan /; dx
f ð0Þ ¼ 1;
ð3:16Þ
where / is the contact angle. A fourth-order polynomial that satisfies (3.16) is given by
1 f ðxÞ ¼ 1 þ 2ð 1Þx2 ð 1Þx4 x2 x2 1 tan /: 2
ð3:17Þ
Imposing the boundary conditions (3.14) as:
hð1; tÞ ¼ ;
@hð1; tÞ df ð1Þ ¼ ; @x dx
3
@ 3 hð1; tÞ d f ð1Þ ¼ ; 3 @x3 dx
ð3:18Þ
we end up solving the system
Ahðjþ1Þ ¼ ChðjÞ ;
ð3:19Þ
E. Momoniat et al. / Applied Mathematics and Computation 217 (2010) 2631–2638
Fig. 2. Numerical solutions admitted by (1.1) obtained from the implicit–explicit scheme (3.19) for n = 3 at t = 0.1. We have chosen m = 250, Dt = 0.1, a = 1 and Dx = 2a/m = 0.008. The initial curve is given by hðx; 0Þ ¼ 1 þ 2ð 1Þx2 ð 1Þx4 12 x2 ðx2 1Þ tan /.
2637
= 0.01,
where
2
1 6 6 1 6 6 1 6 6 6 0 6 6 0 6 6 6 0 6 6 . C ¼ 6 .. 6 6 0 6 6 6 0 6 6 0 6 6 6 0 6 6 4 0 0
0 0
0
0
0 0
0
0
1 0 0 0 0 0 2 0 2 1 0 0
0 0
0 0
0 0
1
0
0 0
0
0
0 0
0
1 0 0
0
0
0 0 .. .. . .
0 .. .
0 1 0 .. .. . . . ..
0 .. .
0 .. .
0 0
0
0
0 1
0
0
0 0
0
0
0 0 1 0
0 0
0
0
0 0
0 1
0 0 0 0
0 0
0 0
0 0 0 0
0 1 0 0
0 0
0
0
0 0
0
0
0 0
0
0 0 0 0
0 0
0
3
7 07 7 07 7 7 0 0 0 07 7 0 0 0 07 7 7 0 0 0 07 7 .. .. .. 7 7: . . . 7 0 0 0 07 7 7 0 0 0 07 7 0 0 0 07 7 7 2 0 2 1 7 7 7 0 0 1 1 5 0 0 0 1
ð3:20Þ
The system (3.19) is the implicit–explicit scheme for finite contact angle /. In Fig. 2 we plot the numerical solution obtained from (3.19) for different values of /. Concluding remarks are made in the next section. 4. Concluding remarks In this paper we have derived unconditionally stable numerical schemes for solving the lubrication Eq. (1.1). The numerical schemes are based on implicit–explicit and Crank–Nicolson implicit–explicit finite difference approximations. The numerical schemes have a truncation error OðDt þ Dx2 Þ. Both numerical schemes are easily implemented using linear iteration. The correctness of the results obtained in this paper are verified by the fact that both the implicit–explicit and Crank–
2638
E. Momoniat et al. / Applied Mathematics and Computation 217 (2010) 2631–2638
Nicolson implicit–explicit finite difference schemes give the same solutions. We have shown that the implicit–explicit scheme is more efficient than the Crank–Nicolson scheme. This is an improvement on the work of Ha et al. [16]. The finite difference schemes presented in this work are easier to implement than the schemes derived by Zhornitskaya and Bertozzi [29]. The accuracy of the numerical schemes can be improved by considering a higher order approximation to the third-order derivative given by
@ 3 hi;j hiþ3;j þ 8hiþ2;j 13hiþ1;j þ 13hi1;j 8hi2;j þ hi3;j þ O Dx4 : @x3 8Dx3
ð4:21Þ
This improves the truncation error to OðDt þ Dx4 Þ. The resulting matrices are however more complicated and would require more computation to solve. Another possible approach that reduces a nonlinear partial differential equation to a linear system is the quasi-linearisation approach introduced by Bellman and Kalaba [3]. This approach does however need iterations to take place at each time step. References [1] R. An, Y. Li, K. Li, Finite element approximation for fourth-order nonlinear problem in the plane, Appl. Math. Comput. 194 (2007) 143–155. [2] W. Barrett, J.F. Blowey, H. Garcke, Finite element approximation of a fourth order nonlinear degenerate parabolic equation, Numer. Math. 80 (1998) 525–556. [3] R.E. Bellman, R.E. Kalaba, Quasilinearization and Non-Linear Boundary Value Problems, American Elsevier, New York, 1965. [4] E. Beretta, M. Bertsch, R. Dal Passo, Nonnegative solutions of a fourth-order nonlinear degenerate parabolic equation, Arch. Ration. Mech. Anal. 129 (1995) 175–200. [5] F. Bernis, J. Hulshof, J.R. King, Dipoles and similarity solutions of the thin film equation in the half-line, Nonlinearity 13 (2000) 413–439. [6] A.J. Bernoff, T.P. Witelski, Linear stability of source-type similarity solutions of the thin film equation, Appl. Math. Lett. 15 (2002) 599–606. [7] A.L. Bertozzi, Symmetric singularity formation in lubrication-type equations for interface motion, SIAM J. Appl. Math. 56 (1996) 681–714. [8] A.L. Bertozzi, The mathematics of moving contact lines in thin liquid films, Not. Am. Math. Soc. 45 (1998) 689–697. [9] M. Bowen, J.R. King, Asymptotic behaviour of the thin film equation in bounded domains, Eur. J. Appl. Math. 12 (2001) 135–157. [10] P. Constantin, T.F. Dupont, R.E. Goldstein, L.P. Kadanoff, M.J. Shelley, S.-M. Zhou, Droplet breakup in a model of the Hele-Shaw cell, Phys. Rev. E 47 (1993) 4169–4181. [11] J. Crank, E. Nicolson, A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type, P. Camb. Philos. Soc. 43 (1947) 50–67. [12] B.R. Duffy, S.K. Wilson, A third-order differential equation arising in thin-film flows and relevant to Tanner’s Law, Appl. Math. Lett. 10 (1997) 63–68. [13] E.B. Dussan, On the spreading of liquids on solid surfaces: static and dynamic contact lines, Annu. Rev. Fluid Mech. 11 (1979) 371–400. [14] H.P. Greenspan, On the motion of a small viscous droplet that wets a surface, J. Fluid Mech. 84 (1978) 125–143. [15] G. Grün, On the convergence of entropy consistent schemes for lubrication type equations in multiple space dimensions, Math. Comput. 72 (2003) 1251–1279. [16] Y. Ha, Y.-J. Kim, T.G. Myers, On the numerical solution of a driven thin film equation, J. Comput. Phys. 227 (2008) 7246–7263. [17] G.R. McGuire, J.Ll. Morris, Explicit–implicit schemes for the numerical solution of nonlinear hyperbolic systems, Math. Comput. 29 (1975) 407–424. [18] S. Middleman, The effect of induced air flow on the spin coating of viscous liquids, J. Appl. Phys. 62 (1987) 2530–2532. [19] E. Momoniat, T.G. Myers, S. Abelman, New solutions for surface tension driven spreading of a thin film, Int. J. Nonlinear Mech. 40 (2005) 523–529. [20] T.G. Myers, Thin films with high surface tension, SIAM Rev. 40 (1998) 441–462. [21] J.A. Moriarty, L.W. Schwartz, E.O. Tuck, Unsteady spreading of thin liquid films with small surface tension, Phys. Fluids A3 (5) (1991) 733–742. [22] A. Oron, S.H. Davis, S.G. Bankoff, Long-scale evolution of thin liquid films, Rev. Mod. Phys. 69 (1997) 931–980. [23] R.D. Passo, H. Garcke, G. Grün, On a fourth-order degenerate parabolic equation: global entropy estimates, existence, and qualitative behavior of solutions, SIAM J. Math. Anal. 29 (1998) 321–342. [24] N.F. Smyth, J.M. Hill, High-order nonlinear diffusion, IMA J. Appl. Math. 40 (1988) 73–86. [25] M.A. Spaid, G.M. Homsy, Stability of Newtonian and viscoelastic dynamic contact lines, Phys. Fluids 8 (1996) 460–478. [26] T.R. Taha, M.J. Ablowitz, Analytical and numerical aspects of certain nonlinear evolution equations: II. Numerical, nonlinear Schrodinger equation, J. Comput. Phys. 55 (1984) 203–230. [27] L.H. Tanner, The spreading of silicone oil drops on horizontal surfaces, J. Phys. D. Appl. Phys. 12 (1979) 1473–1484. [28] T.P. Witelski, Similarity solutions of the lubrication equation, Appl. Math. Lett. 10 (1997) 107–113. [29] L. Zhornitskaya, A.L. Bertozzi, Positivity-preserving numerical schemes for lubrication-type equations, SIAM J. Numer. Anal. 37 (1999) 523–555.
