Journal of Applied Mathematics

Volume 2013 (2013), Article ID 181939, 6 pages

http://dx.doi.org/10.1155/2013/181939

Research Article

Numerical Investigation of the Steady State of a Driven Thin Film Equation

A. J. Hutchinson, C. Harley, and E. Momoniat

Centre for Differential Equations, Continuum Mechanics and Applications, School of Computational and Applied Mathematics, University of the Witwatersrand, Private Bag 3, Wits Johannesburg 2050, South Africa

Received 6 November 2012; Accepted 29 December 2012

Academic Editor: Mehmet Pakdemirli

Copyright © 2013 A. J. Hutchinson 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.

Abstract.

A third-order ordinary differential equation with application in the flow of a thin liquid film is considered. The boundary conditions come from Tanner's problem for the surface tension driven flow of a thin film. Symmetric and nonsymmetric finite difference schemes are implemented in order to obtain steady state solutions. We show that a central difference approximation to the third derivative in the model equation produces a solution curve with oscillations. A difference scheme based on a combination of forward and backward differences produces a smooth accurate solution curve. The stability of these schemes is analysed through the use of a von Neumann stability analysis.

1. Introduction

In this paper we investigate numerical solutions of the third-order ordinary differential equation (ODE)

where is a constant to be determined. Equation (1) is the steady state case of the partial differential equation

where the flux is given by

This equation describes the flow of a thin liquid film, where denotes the film thickness. The flux terms represent surface shear and gravity, where the forces act in opposing directions, and the diffusion term on the right-hand side represents surface tension. The surface shear term may arise due to temperature or concentration gradients or to an external shear force. Equation (2) with (3) is a special case of the more general partial differential equation

where , , and are nonnegative parameters. Equation (4) models the flow of a thin film of fluid with thickness above an inclined plane. The variable is the distance down the plane, and the variable represents time. The parameters , , and describe the effects of surface tension, gravity, and the angle of inclination of the plane, respectively. This model is derived in the work by Bertozzi et al. [1], and numerical solutions are obtained and investigated. In our model, we consider a vertical plane, and thus . With and by considering the steady state solution, that is, , (4) reduces to

Dividing (5) by and, for simplicity, replacing and with and , respectively, (5) becomes

which can be rewritten as follows:

In (7), the parameter is proportional to Marangon's number which is a ratio of the surface tension to viscous forces, and the parameter is proportional to the Bond number which is a ratio of gravity to viscous forces [2]. The more general equation,

was solved numerically for the case where the arbitrary constant of integration was set equal to , and the asymptotic solution was investigated in [2].

We can nondimensionalize the equation using the following transformations:

and after integrating once with respect to , (8) becomes

The case where has been studied in [2]. We will consider (10) for the case , where the value of will be determined through the application of one of the boundary conditions. In this paper, we use the conditions provided by Tanner's problem [3] in the spreading of silicon oil drops on a horizontal surface modelled as a thin film. These conditions are given by

In his work Tanner [3] investigated surface tension-driven spreading of an oil drop modelled as a thin film. The problem is closed by the boundary condition

where is a positive constant. The autonomous nature of (10) makes the specification of arbitrary. The boundary condition (12) describes the physical phenomenon of touchdown [4] or, more precisely, the point at which the free surface of the thin film makes contact with the solid substrate. Given that (1) is singular at , the boundary condition (12) is imposed as

where is the height of the precursor film. A discussion of this approximation can be found in Bertozzi [5, 6].

In this paper, we use a symmetric and nonsymmetric finite difference scheme to solve (1). This method was applied in [4] to the following equation:

In [7] a shooting method was applied to (14), and the asymptotic solution was also studied. In this research, we will conduct a von Neumann stability analysis as a means of investigating the stability of the methods employed.

This paper is set out as follows: in Section 2.1 we implement a symmetric finite difference approximation to the third derivative with truncation error . We see that this scheme is zero stable and a spectral analysis gives inconclusive results. Thus, a von Neumann stability analysis is implemented. The solution obtained via the symmetric difference scheme contains oscillations. In Section 2.2 we implement a non-symmetric difference approximation to the third derivative. This scheme has truncation error . It is shown to be zero stable; however, a von Neumann stability analysis is also implemented. The solution of the non-symmetric difference scheme does not contain any oscillations.

We do not rigorously prove that the schemes investigated in this paper converge, but rather show that the schemes are stable. Furthermore, we show that convergence does take place by evaluating the absolute error at each iteration. We show, in a similar fashion to Momoniat [4], that even if a numerical scheme is stable, the boundary conditions produce results which are physically unstable. In this instance, this instability manifests as oscillations in the solution. Furthermore, stability and accuracy should not necessarily be accessed based upon the truncation error of the scheme: in this case the most accurate and stable scheme is the non-symmetric scheme with a truncation error of and a stability criteria of , instead of and as for the symmetric scheme where the solution curve exhibited oscillations.

2. Finite Difference Methods

In the sections that follow, we let and . The choice of comes from a consideration of the equation at , where

We may claim; however, that which means that since . Doing so makes sense given the implementation of the physical condition that the volume flux at vanishes, where is the height of the precursor film. Thus, we may choose such that (10) becomes

2.1. Symmetric Difference Approximation

We approximate (1) subject to the boundary conditions given by (11) and (13) by implementing a symmetric difference approximation. We define to be the approximate value of at a point , where , , , and , and the step length is given by , where is a positive integer. The first derivative at a point is approximated as

and the second derivative at a point is approximated as

The third derivative is approximated by a central difference approximation,

and using (18) we obtain

Simplifying this we get

Approximating (1) using (21) and after simplifying, we obtain a nonlinear difference equation,

with a truncation error of . We now define to be the th iteration value to . Using this notation, (22) can be evaluated iteratively in the following manner:

The boundary conditions are incorporated by requiring that satisfies (11) and (13). The boundary conditions and are implemented by

When evaluating (23) at the point , a value for is required. In thin film theory, it is reasonable to suppose that the height of the fluid at position is the same at [4]. In other words, we can obtain a value for by setting

Also, in order to calculate using (25), a value for is required. A central difference approximation is used to approximate the boundary condition ,

and setting (26) equal to zero, we have

The difference equation (23), together with the boundary conditions (24), (25), and (27), can be written more compactly in matrix form as follows:

where

An initial guess, , to that satisfies the boundary conditions (11) and (13) is given by

The corresponding homogeneous equation of (23) is given by

This equation has been studied extensively in [4]. This difference scheme is zero stable, and the results of a spectral analysis were found to be inconclusive with regard to the convergence of the difference scheme. Plots of (31) were found to contain oscillations.

Given that we need to consider an alternative form of analysis we choose to implement the von Neumann stability analysis. In order to do so we first need to linearize our scheme (23). We notice that the difference scheme is linear in except for the term. Rewriting (23) we obtain

Let

and by making the substitution (33) in our difference scheme (32) we obtain

which is now linear in .

We now implement a Von Neumann stability analysis by making the substitution,

into (34), where and is a constant. The simplification that follows here can be found in [4], where is given by

whereas, in this paper, is given by (33). In [4] ; however, the sign of defined by (33) is not as straightforward to determine. Following the same procedure as described in [4] we obtain the following expression for the amplification factor:

This can be further simplified, resulting in an upper bound for given by

In order to obtain the maximum and minimum values of , consider the plots of as a function of given in Figures 1 and 2. We see that the minimum value for is obtained at and is zero. In order to obtain the value at which is a maximum, we find the second derivative of and set it equal to zero. Solving this in Mathematica, we find that attains its maximum value at . The value of is . In other words, we have deduced the following results:

Now always satisfies (38). Set to obtain

The von Neumann stability criteria is satisfied if (40) is satisfied. We can, therefore, conclude that the iteration (32) is conditionally stable, provided the value of the step length remains below the maximum value of . In Table 1 we show the absolute error at each iteration. After five iterations the absolute error max is . We plot the results after five iterations in Figure 3. In both Table 1 and Figure 3 we have chosen , where , and thus, the choice of satisfies the von Neumann stability criteria.

Table 1: Table showing the absolute error after each iteration of (28), where .

Error

max

1

2

3

4

5

Figure 1: Plot of for .

Figure 2: Plot of for .

Figure 3: Solution to the symmetric difference scheme for .

The solution obtained via a symmetric difference scheme given by (23) oscillates. These oscillations diminish as the value of the step length is reduced. However, even for there are still oscillations present. This seems to indicate that the oscillations occur in the homogeneous part of the difference scheme. Similar to [4] we believe these oscillations occur as a result of imposing the additional boundary condition necessary to iterate the numerical scheme. In the next section we show that a non-symmetric difference approximation provides smooth accurate solutions even though the truncation error is higher than the truncation error for the symmetric difference scheme.

2.2. Nonsymmetric Difference Approximation

An alternative method to approximate the third-order ODE (1) is considered below. Forward and backward difference approximations are used to approximate the first- and second-order derivatives instead of the central (symmetric) approximations considered above. The approximations are given by

Substituting (45) into (1) results in the following difference equation:

which has a truncation error of . After multiplying by we obtain

Equation (47) can be evaluated iteratively as

To check for 0 stability, we consider the following equation:

which is the homogeneous part of (47). This equation has been studied extensively in [4]. This scheme was proven to be zero stable, and no oscillations were present in the solution curve. We can write the above system in the form

where

An initial guess satisfying the boundary conditions (11) and (13) (as before) is

A von Neumann stability analysis was used to determine the stability of the system (48). To perform the analysis equation (48) was linearized as follows:

where

In [4] the value assigned to is . The analysis follows the same procedure with the value for defined by (54) replacing the value for in [4]. The amplification factor, , is given by

An upper bound for is given by

Using the method for determining the maximum and minimum values for from the previous section, we have that a similar computation for yields

The von Neumann stability criteria give that

which shows that the scheme is conditionally stable. In Table 2 we show the absolute error max at each iteration for six iterations—after six iterations the absolute error is . As done previously, in both Table 2 and Figure 3 we have chosen , where such that the choice of   satisfies the von Neumann stability criteria. A plot of the nonsymmetric difference scheme is given in Figure 4 where it can be seen that there are no oscillations present.

Table 2: Table showing the absolute error after each iteration of (50), where .

Error

max

1

2

3

4

5

6

Figure 4: Solution to the nonsymmetric scheme for .

3. Concluding Remarks

In order to investigate the steady state solutions of the partial differential equation (2) we integrate with respect to . In doing so an arbitrary constant arises, denoted by , the value of which was deduced by using the physical condition that the volume flux at vanishes. Symmetric and non-symmetric finite difference schemes were implemented in order to obtain steady state solutions. In this paper, we use 0 stability to show that the finite difference schemes used to approximate the third-order derivative in the model equation are stable. A von Neumann stability analysis is implemented to obtain criteria on the step length which makes the iterative evaluation of the nonlinear difference scheme stable.

The symmetric difference scheme has a truncation error and is shown to be 0-stable. The nonlinear difference scheme is solved iteratively and displays unstable results even though the step length maintains the von Neumann stability criteria. These oscillations persist even as is reduced. The non-symmetric difference scheme, which has a truncation error of and is 0 stable, in turn provided a stable solution. The results obtained in this paper reiterate the fact that finite difference approximations with a small truncation error should not thoughtlessly be used to solve problems in which the boundary conditions affect the behaviour of the solution.

Acknowledgment

C. Harley acknowledges support from the National Research Foundation, South Africa, under Grant no. 79184.

References

A. L. Bertozzi, A. Münch, and M. Shearer, “Undercompressive shocks in thin film flows,” Physica D, vol. 134, no. 4, pp. 431–464, 1999.

E. Momoniat, “A nonlinear eigenvalue problem from thin-film flow,” Journal of Engineering Mathematics, 2012.

L. H. Tanner, “The spreading of silicone oil drops on horizontal surfaces,” Journal of Physics D, vol. 12, pp. 1473–1484, 1979.

E. Momoniat, “Numerical investigation of a third-order ODE from thin film flow,” Meccanica, vol. 46, no. 2, pp. 313–323, 2011.

A. L. Bertozzi, “Symmetric singularity formation in lubrication-type equations for interface motion,” SIAM Journal on Applied Mathematics, vol. 56, no. 3, pp. 681–714, 1996.

A. L. Bertozzi, “The mathematics of moving contact lines in thin liquid films,” Notices of the American Mathematical Society, vol. 45, no. 6, pp. 689–697, 1998.

E. Momoniat, “On the determination of the steady film profile for a non-Newtonian thin droplet,” Computers & Mathematics with Applications, vol. 62, no. 1, pp. 383–391, 2011.