Numerical Investigation of Thin Film Spreading Driven by Surfactant Using Upwind Schemes

1 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 2Mechanical Engineering Department, Engineering Faculty of Bu-Ali Sina University, Hamedan, Iran 3University of Michigan-Shanghai Jiao Tong University Joint Institute, Shanghai Jiao Tong University, Shanghai, China

The free surface of the thin film is given by ℎ = ℎ( , , ) and the surfactant concentration by Γ = Γ( , , ). The nondimensional constants and balance gravity, capillarity, and Marangoni stress, respectively. The nondimensional constant = 1/ where is the Peclet number. Surfactants have important applications in both industrial and biological applications.
The reviews of Craster and Matar [3] and Afsar-Siddiqui et al. [4] contain applications of surfactant and thin film flow problems to industry. These include thin film spreading driven by surfactant [5][6][7] and the use of surfactants in industrial coating [8,9]. Warner et al. [10] study the effects of surfactants on the dewetting of ultrathin films. This is important for templating of films in microelectronics [11,12]. Halpern et al. [13] perform a theoretical study of the delivery of surfactant into the lung. This is part of a much larger investigation dealing with the delivery of treatment for neonate and adult respiratory distress syndrome [14][15][16][17][18][19][20][21]. Extensions to this works have been undertaken by Craster and Matar [22] and Matar et al. [23] by modelling the effect of surfactants on a layer of mucus modelled as a non-Newtonian fluid. This is a significant improvement of the Newtonian models considered previously. Other biological models in which the use of surfactants is relevant are applications to the human eye [24]. Schwartz et al. [25] model cell division and motility as a consequence of the interaction of surfactants with the free surface of the cell membrane.
In this paper we consider numerical solutions of the onedimensional case ℎ = ℎ( , ) and Γ = Γ( , ) where surface tension effects are ignored. The resulting coupled system of nonlinear equations is given by The model equations (3) have been derived in [17,26,27] where surface tension effects have been ignored. Levy and Shearer [2] have considered numerical solutions of the coupled system (3) by implementing an implicit scheme and solving the resulting equations using a Newton's method.
Of particular interest in the numerical investigation are the shock type structures that develop in the travelling wave solutions. Peterson and Shearer [1] consider numerical solutions of a one-dimensional case of (3) for = = 0 and a linear equation of state for the surface tension. They consider a front tracking numerical scheme and an implicit numerical scheme where the resulting equations are solved using a Newton's method. A front capturing method is also introduced based on Godunov's method which is very effective in capturing the moving front. Peterson and Shearer then go on to consider two-dimensional spreading for + > 0 which is not relevant to our purposes. Peterson and Shearer [28] make analytical progress in solving (1) by investigating similarity solutions for = = = 0. In this paper we use an explicit upwind numerical scheme of Roe [29,30] to solve the coupled system (3). Upwind schemes are typically implemented to solve hyperbolic partial differential equations while the coupled system (3) is hyperbolic/degenerate-parabolic [1]. For = 0 (1) is degenerate at Γ = 0 [1]. Peterson and Shearer [1] have pointed out that the degeneracy for the case = 0 implies that if Γ(x, 0) has compact support, then the solution Γ(x, ) will have compact support for > 0. This holds true for the case = 1 and = 0 for a linear equation of state (Γ) = 1 − Γ considered in this paper. It is this compact support that makes it possible to implement an upwind scheme to solve the hyperbolic/degenerate-parabolic coupled system. Another advantage of using the Roe scheme is, as Jaisankar and Raghurama Rao [31] have shown, that the approximation to the flux gradient used in the Roe scheme is related to the speed of the shock through the Rankine-Hugoniot condition. This ties in with the front tracking and front capturing methods implemented by Peterson and Shearer [1]. The standard Roe formulation is compared to a Roe-Sweby scheme [32] with a minmod limiter. An advantage of an explicit scheme over an implicit scheme is the ease with which we can implement the scheme in parallel. We implement the explicit scheme using OpenMP with C++.
A further improvement we make to the numerical scheme is to improve the order of the time integration. TVD (total variation diminishing) Runge-Kutta schemes [33,34] are popular high-order time integration schemes. At each time step more than one iteration of the numerical scheme is required for an improvement in the accuracy of the scheme. This can prove to be computationally expensive for a coupled system of nonlinear equations. We compare a TVD Runge-Kutta approximation to the time derivative to a second-order approximation to the time derivative given by a BDF (backward difference formula) leading to an A-stable multistep method [35]. Unlike the TVD Runge-Kutta scheme the BDF scheme requires only one Roe scheme evaluation at each time step.
The paper is divided up as follows. In Section 2 we derive the numerical scheme to solve the coupled system of nonlinear equations (3). In Section 3 we consider the stability of an FTCS scheme which gives the stability of the equivalent upwind scheme. In Section 4 we consider simulations of the numerical scheme. Concluding remarks are made in Section 5.

Upwind Numerical Scheme
We define ℎ = ℎ( , ) and Γ = Γ( , ) where the spatial domain is discretized into + 1 equidistant intervals of width △ and = △ . The time is defined by = △ where △ is the time step length. We approximate the spatial derivatives in (4) by the low-order forward and backward difference approximations the central difference approximation and the high-order central difference approximation Mathematical Problems in Engineering 3 We therefore obtain the approximations to the fluxes given by An Euler forward approximation to the time derivative is given by A BDF approximation to the time derivative is given by A third-order TVD Runge-Kutta approximation to the time derivative to the free surface ℎ( , ) yields The third-order TVD Runge-Kutta approximation to the surfactant concentration will have a similar form. A finite volume approximation to (3) is given by The fluxes are approximated by the three-point central difference schemes with a numerical viscosity term of Roe [29,30] given by where Using a finite volume approximation to the spatial derivatives in (3) with an Euler approximation to the time derivatives we obtain A finite volume approximation to the spatial derivatives in (3) combined with the BDF approximation to the time derivative given by (10) leads to the numerical scheme 4

Mathematical Problems in Engineering
We implement the multistep method by evaluating (15) for = 0 to determine ℎ 1 and Γ 1 . We evaluate (16) for = 1, 2, . . ., where we use the values for ℎ 1 and Γ 1 obtained from (15) to start the scheme. Both (15) and (16) are evaluated for = 1, 2, . . . , − 1. The values at = 0 and = are determined from the boundary conditions. The initial conditions come from the physics of the problem being considered.
The Roe approximations to the flux given by (13) coupled with the third-order Runge-Kutta approximation to the time derivative lead to a Roe-TVD numerical scheme. A Roe-Sweby scheme [32] combines the high-order Roe flux approximation (13) with the low-order Lax-Wendroff approximation such that * After some numerical experimentation we find that the limiter that works best is the minmod limiter given by The Roe-Sweby scheme is coupled with the third-order TVD approximation to the time derivative leading to a Roe-Sweby-TVD numerical scheme. The final scheme we consider is the Roe-BDF scheme given by (16) where the approximations to the fluxes are given by (13). The Roe-TVD, Roe-Sweby-TVD, and Roe-BDF numerical schemes are simulated in the next section.
We consider an initial exponential film profile given by We are interested in simulating the behaviour of a moving front at the leading edge of a thin film. We thus fix the height of the film and surfactant concentration at the origin such that ℎ (0, ) = 1, Γ (0, ) = 1.
The continuity boundary conditions close the problem. The boundary conditions (23) coupled with (24) give The model equation (3) is degenerate since the coefficient of the highest derivative tends to zero as the film height tends to zero [36]. To avoid numerical difficulties which arise out of this degeneracy for the boundary conditions at = ∞ we assume the height of the film is not zero but rather the height of a precursor film [37]. In this paper we fix the height of the precursor film to be the height of the initial film profile at the end point. We make a similar assumption with the surfactant concentration. This precursor boundary condition coupled with the derivative condition at = ∞ gives

Stability
In this section we consider the stability of the system (3). We linearise the system around the constant solutions ℎ = 1 and Γ = 1 by making the substitutions ℎ = 1 + ℎ 0 ( , ) , where ≪ 1. Substituting (27) into (3) and separating to leading order coefficients in we obtain the linear system We now perform a von Neumann stability analysis on the linear system (28) to determine the Courant-Friedrichs-Lewy (CFL) condition for numerical stability. We approximate the spatial derivatives in (28) by central difference approximations Mathematical Problems in Engineering 5 We approximate the time derivatives by the forward difference approximations We obtain the forward-time central-space (FTCS) approximation to (28) given by To perform a von Neumann stability analysis we substitute into (32) where 2 = −1 to obtain the system ) .
In terms of a von Neumann stability analysis we need to show that the amplification factor | +1 / | < 1 and | +1 / | < 1. Alternatively, we need to show that the iterative system (34) is bounded. We can show this by showing that the spectral radius ( ) satisfies the condition The matrix admits the eigenvalues The spectral radius, ( ), is defined as We find that Condition (36) therefore gives The nonlinear equation (40) can be simplified to give the Courant-Friedrichs-Lewy (CFL) condition where we have chosen = ( + 2 )/ △ for ∈ I. In this section we have shown that the FTCS numerical scheme given by (31) and (32) will give stable results provided the CFL condition (41) is satisfied. This condition holds true for the Roe upwind scheme derived in Section 2. In the next section we consider the evolution of the initial exponential thin film profile for both endogenous and exogenous surfactants.

Simulation of Numerical Scheme
We simulate the numerical schemes on a 2.66 GHz Windows machine with 6.00 GB of RAM with 2 Intel QUAD core CPUs using the MinGW version of C++. We first consider a simulation of all three numerical schemes taking = 250 and △ = 2.5 × 10 −6 where we have chosen 0 = 0, = 4 and ℎ = ( − 0 )/ . We take our iterations up to a final time of = 1.0.
The approximations ( , ) and * ( , ) correspond to numerical solutions from two different schemes. We tabulate our results in Table 1. From the results in Table 1 we note that the Roe-BDF approach produces the smallest error for the 2 and ∞ norms for both an endogenous surfactant as well as for an exogenous surfactant. We interpret these results to mean that the Roe-BDF approach will produce the most stable results with the least variation for long-time simulations.
We next consider long-time simulations at high resolution where we take = 500 and △ = 2.5 × 10 −5 . This implies that △ / △ 2 = 0.390625 and △ / △ 2 = 0.003125. We simulate this case in parallel using OpenMP. In Table 2 we show the advantage of using OpenMP over nonparallel computations for long-time simulations of an endogenous surfactant. From Table 2 we note that implementing the numerical scheme on OpenMP is faster than using the nonparallel formulation. While the improvement is not significant for the times considered here, when one considers an increase in the resolution (larger values) and long-time solutions, a significant improvement is noted.
In Figures 1 and 2 we obtain the propagating front solutions observed by Levy and Shearer [2]. We have obtained these solutions without recourse to an implicit numerical Mathematical Problems in Engineering

Concluding Remarks
In this paper we have considered numerical solutions of a coupled system of nonlinear equations modelling the effects of surfactant on the free surface of a thin film on a horizontal substrate. The original system has been investigated by Levy and Shearer [2] and solved by implementing an implicit numerical scheme coupled with Newton's method. We have shown that an explicit upwind scheme coupled with highorder approximations to the spatial derivatives and a BDF scheme for the time integration in which the time step and space step satisfy the CFL condition △ / △ 2 < 1/2 are able to capture the sharp changes in gradient that occur in both the free surface profile and surfactant concentration. The upwind schemes are easier to implement than implicit schemes and can be implemented in parallel.