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.
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.
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.