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

Abstract

Numerical solutions of a coupled system of nonlinear partial differential equations modelling the effects of surfactant on the spreading of a thin film on a horizontal substrate are investigated. A CFL condition is obtained from a von Neumann stability analysis of a linearised system of equations. Numerical solutions obtained from a Roe upwind scheme with a third-order TVD Runge-Kutta approximation to the time derivative are compared to solutions obtained with a Roe-Sweby scheme coupled to a minmod limiter and a TVD approximation to the time derivative. Results from both of these schemes are compared to a Roe upwind scheme and a BDF approximation to the time derivative. In all three cases high-order approximations to the spatial derivatives are employed on the interior points of the spatial domain. The Roe-BDF scheme is shown to be an efficient numerical scheme for capturing sharp changes in gradient in the free surface profile and surfactant concentration. Numerical simulations of an initial exponential free surface profile coupled with initial surfactant concentrations for both exogenous and endogenous surfactants are considered.

1. Introduction

In this paper we investigate numerical solutions of a coupled system of hyperbolic/degenerate-parabolic equations [1] modelling the spreading of an insoluble surfactant on the free surface of a thin liquid film. Surfactants are known to decrease the effects of surface tension by creating a spatial variation in the surface tension due to a tangential surface stress. This effect is also known as a Marangoni stress. The coupled system of nonlinear equations is given by [1, 2] ℎ_𝑡+∇⋅(𝐯ℎ)=0,Γ_𝑡+∇⋅(𝐰Γ)=𝛿△Γ, (1) where 1𝐯=2(Γ)−𝛽ℎ∇𝜎3ℎ²𝜅∇ℎ+3ℎ²(Γ)−𝛽∇△ℎ,𝐰=ℎ∇𝜎2ℎ²𝜅∇ℎ+2ℎ²∇△ℎ. (2)

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–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–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 one-dimensional case ℎ=ℎ(𝑥,𝑡) and Γ=Γ(𝑥,𝑡) where surface tension effects are ignored. The resulting coupled system of nonlinear equations is given by 𝜕ℎ+𝜕𝑡𝜕𝑄𝜕𝑥=0,𝜕Γ+𝜕𝑡𝜕𝑃𝜕𝑥=0, (3) where ℎ𝑄=−²2𝜕Γ+ℎ𝜕𝑥3³,𝑃=−ℎΓ𝜕Γ+ℎ𝜕𝑥2²Γ. (4)

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 Γ(𝐱,0) has compact support, then the solution Γ(𝐱,𝑡) 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.

2. 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 𝜕Γ𝜕𝑥_{𝑥=𝑥𝑖}≈Γ^𝑗_𝑖+1−Γ^𝑗_𝑖(),△𝑥+𝒪△𝑥𝜕Γ𝜕𝑥_{𝑥=𝑥𝑖}≈Γ^𝑗_𝑖−Γ^𝑗_𝑖−1(),△𝑥+𝒪△𝑥 (5) the central difference approximation 𝜕Γ𝜕𝑥_{𝑥=𝑥𝑖}≈Γ^𝑗_𝑖+1−Γ^𝑗_𝑖−12△𝑥+𝒪△𝑥², (6) and the high-order central difference approximation 𝜕Γ𝜕𝑥_{𝑥=𝑥𝑖}≈−Γ^𝑗_𝑖+2+8Γ^𝑗_𝑖+1−8Γ^𝑗_𝑖−1+Γ^𝑗_𝑖−212△𝑥+𝒪△𝑥⁴. (7) We therefore obtain the approximations to the fluxes given by 𝑄^𝑗_𝑖=−ℎ^𝑗_𝑖²2Γ^𝑗_𝑖+1−Γ^𝑗_𝑖+ℎ△𝑥3^𝑗𝑖3,−ℎ𝑖=0,_𝑗𝑖22Γ^𝑗_𝑖+1−Γ^𝑗_𝑖−1+ℎ2△𝑥3^𝑗𝑖3,−ℎ1≤𝑖≤2,𝑛−2≤𝑖≤𝑛−1,_𝑗𝑖22−Γ^𝑗_𝑖+2+8Γ^𝑗_𝑖+1−8Γ^𝑗_𝑖−1+Γ^𝑗_𝑖−2+ℎ12△𝑥3^𝑗𝑖3,−ℎ3≤𝑖≤𝑛−3,_𝑗𝑖22Γ^𝑗_𝑖−Γ^𝑗_𝑖−1+ℎ△𝑥3^𝑗𝑖3,𝑃𝑖=𝑛,_𝑗𝑖=−ℎ^𝑗_𝑖Γ^𝑗_𝑖Γ^𝑗_𝑖+1−Γ^𝑗_𝑖+ℎ△𝑥^𝑗𝑖2Γ^𝑗_𝑖2,𝑖=0,−ℎ^𝑗_𝑖Γ^𝑗_𝑖Γ^𝑗_𝑖+1−Γ^𝑗_𝑖−1+ℎ2△𝑥^𝑗𝑖2Γ^𝑗_𝑖2,1≤𝑖≤2,𝑛−2≤𝑖≤𝑛−1,−ℎ^𝑗_𝑖Γ^𝑗_𝑖−Γ^𝑗_𝑖+2+8Γ^𝑗_𝑖+1−8Γ^𝑗_𝑖−1+Γ^𝑗_𝑖−2+ℎ12△𝑥^𝑗𝑖2Γ^𝑗_𝑖2,3≤𝑖≤𝑛−3,−ℎ^𝑗_𝑖Γ^𝑗_𝑖Γ^𝑗_𝑖−Γ^𝑗_𝑖−1+ℎ△𝑥^𝑗𝑖2Γ^𝑗_𝑖2,𝑖=𝑛. (8)

An Euler forward approximation to the time derivative is given by 𝜕ℎ𝜕𝑡_{𝑡=𝑡𝑗}≈ℎ_𝑖^𝑗+1−ℎ^𝑗_𝑖(),△𝑡+𝒪△𝑡𝜕Γ𝜕𝑡_{𝑡=𝑡𝑗}≈Γ_𝑖^𝑗+1−Γ^𝑗_𝑖().△𝑡+𝒪△𝑡 (9) A BDF approximation to the time derivative is given by 𝜕ℎ𝜕𝑡_{𝑡=𝑡𝑗}≈3ℎ_𝑖^𝑗+1−4ℎ^𝑗_𝑖+ℎ_𝑖^𝑗−12△𝑡+𝒪△𝑡²,𝜕Γ𝜕𝑡_{𝑡=𝑡𝑗}≈3Γ_𝑖^𝑗+1−4Γ^𝑗_𝑖+Γ_𝑖^𝑗−12△𝑡+𝒪△𝑡². (10) A third-order TVD Runge-Kutta approximation to the time derivative to the free surface ℎ(𝑥,𝑡) yields ℎ_𝑖^(𝑗+1/3)=ℎ^𝑗_𝑖−△𝑡𝑄△𝑥^{𝑗()𝑖+1/2}−𝑄^𝑗_()𝑖−1/2,ℎ_𝑖^(𝑗+2/3)=34ℎ^𝑗_𝑖+14ℎ_𝑖^(𝑗+1/3)−14△𝑡𝑄△𝑥^{()(𝑗+1/3)𝑖+1/2}−𝑄₍₎^(𝑗+1/3)_𝑖−1/2,ℎ_𝑖^𝑗+1=13ℎ^𝑗_𝑖+23ℎ_𝑖^(𝑗+2/3)−23△𝑡𝑄△𝑥^{()(𝑗+2/3)𝑖+1/2}−𝑄₍₎^(𝑗+2/3)_𝑖−1/2. (11)

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 𝜕ℎ𝜕𝑡_{𝑡=𝑡𝑗}+1𝑄△𝑥^{𝑗()𝑖+1/2}−𝑄^𝑗_()𝑖−1/2=0,𝜕Γ𝜕𝑡_{𝑡=𝑡𝑗}+1𝑃△𝑥^{𝑗()𝑖+1/2}−𝑃^𝑗_()𝑖−1/2=0. (12)

The fluxes are approximated by the three-point central difference schemes with a numerical viscosity term of Roe [29, 30] given by 𝑄^𝑗_()𝑖+1/2=12𝑄^𝑗_𝑖+𝑄^𝑗_𝑖+1−𝑎^𝑗_()𝑖+1/2ℎ^𝑗_𝑖+1−ℎ^𝑗_𝑖,𝑃^𝑗_()𝑖+1/2=12𝑃^𝑗_𝑖+𝑃^𝑗_𝑖+1−𝑏^𝑗_()𝑖+1/2Γ^𝑗_𝑖+1−Γ^𝑗_𝑖, (13) where 𝑎^𝑗_()𝑖+1/2=𝑄^𝑗_𝑖+1−𝑄^𝑗_𝑖ℎ^𝑗_𝑖+1−ℎ^𝑗_𝑖,ℎ^𝑗_𝑖+1≠ℎ^𝑗_𝑖,ℎ_𝑖Γ_𝑖+1−Γ_𝑖−ℎ△𝑟^𝑗𝑖2,ℎ^𝑗_𝑖+1=ℎ^𝑗_𝑖,𝑏^𝑗_()𝑖+1/2=𝑃^𝑗_𝑖+1−𝑃^𝑗_𝑖Γ^𝑗_𝑖+1−Γ^𝑗_𝑖,Γ^𝑗_𝑖+1≠Γ^𝑗_𝑖,−ℎ^𝑗_𝑖²2,Γ^𝑗_𝑖+1=Γ^𝑗_𝑖. (14) Using a finite volume approximation to the spatial derivatives in (3) with an Euler approximation to the time derivatives we obtain ℎ_𝑖^𝑗+1=ℎ^𝑗_𝑖+𝜆𝑄^𝑗_()𝑖+1/2−𝑄^𝑗_()𝑖−1/2,Γ_𝑖^𝑗+1=Γ^𝑗_𝑖+𝜆𝑃^𝑗_()𝑖+1/2−𝑃^𝑗_()𝑖−1/2. (15) 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 ℎ_𝑖^𝑗+1=43ℎ^𝑗_𝑖−13ℎ_𝑖^𝑗−1+2𝜆3𝑄^𝑗_()𝑖+1/2−𝑄^𝑗_()𝑖−1/2,Γ_𝑖^𝑗+1=43Γ^𝑗_𝑖−13Γ_𝑖^𝑗−1+2𝜆3𝑃^𝑗_()𝑖+1/2−𝑃^𝑗_()𝑖−1/2. (16)

We implement the multistep method by evaluating (15) for 𝑗=0 to determine ℎ¹_𝑖 and Γ¹_𝑖. We evaluate (16) for 𝑗=1,2,…, where we use the values for ℎ¹_𝑖 and Γ¹_𝑖 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 𝑄^𝑗₍₎^𝐿𝑊_𝑖+1/2=12𝑄^𝑗_𝑖+𝑄^𝑗_𝑖+1−△𝑡𝑎△𝑥^{𝑗()𝑖+1/22}ℎ^𝑗_𝑖+1−ℎ^𝑗_𝑖,𝑃^𝑗₍₎^𝐿𝑊_𝑖+1/2=12𝑃^𝑗_𝑖+𝑃^𝑗_𝑖+1−△𝑡𝑏△𝑥^{𝑗()𝑖+1/22}Γ^𝑗_𝑖+1−Γ^𝑗_𝑖, (17) such that 𝑄^∗𝑗_()𝑖+1/2=𝑄^𝑗_()𝑖+1/2+Φ(𝑟)𝑄^𝑗₍₎^𝐿𝑊_𝑖+1/2−𝑄^𝑗_()𝑖+1/2,𝑃^∗𝑗_()𝑖+1/2=𝑃^𝑗_()𝑖+1/2+Φ(𝑤)𝑃^𝑗₍₎^𝐿𝑊_𝑖+1/2−𝑃^𝑗_()𝑖+1/2,ℎ𝑟=^𝑖−ℎ_𝑖−1ℎ_𝑖+1−ℎ_𝑖Γ,𝑤=^𝑖−Γ_𝑖−1Γ_𝑖+1−Γ_𝑖. (18) After some numerical experimentation we find that the limiter that works best is the minmod limiter given by Φ(𝑟)=max[0,min(1,𝑟)]. (19)

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.

The surfactant concentration has two forms. For a preexisting or endogenous surfactant Γ(𝑥,0)=1, (20) and for a deposited or exogenous surfactant [17, 26, 27] Γ(𝑥,0)=exp−𝑥². (21) We consider an initial exponential film profile given by ℎ(𝑥,0)=exp−𝑥². (22) 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. (23) The continuity boundary conditions 𝜕ℎ𝜕𝑥_{𝑥=0,𝑥=∞}=0,𝜕Γ𝜕𝑥_{𝑥=0,𝑥=∞}=0 (24) close the problem. The boundary conditions (23) coupled with (24) give ℎ^𝑗₀=1,ℎ^𝑗₁Γ=1,_𝑗0=1,Γ^𝑗₁=1. (25) 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 ℎ_𝑛^𝑗+1=ℎ^𝑗_𝑛,ℎ^𝑗+1_𝑛−1=ℎ^𝑗_𝑛−1,Γ_𝑛^𝑗+1=Γ^𝑗_𝑛,Γ^𝑗+1_𝑛−1=Γ^𝑗_𝑛−1. (26)

3. 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+𝜖ℎ₀(𝑥,𝑡),Γ=1+𝜖Γ₀(𝑥,𝑡), (27) where 𝜖≪1. Substituting (27) into (3) and separating to leading order coefficients in 𝜖 we obtain the linear system 𝜕ℎ₀+𝜕𝑡𝜕ℎ₀−1𝜕𝑥2𝜕²𝑔₀𝜕𝑥²=0,𝜕𝑔₀+1𝜕𝑡2𝜕𝑔₀+𝜕𝑥𝜕ℎ₀−𝜕𝜕𝑥²𝑔₀𝜕𝑥²=0. (28) 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 𝜕ℎ₀𝜕𝑥_{𝑥=𝑥𝑖}≈ℎ₀^𝑗_𝑖+1−ℎ₀^𝑗_𝑖−12△𝑥+𝒪△𝑥²,𝜕𝑔₀𝜕𝑥_{𝑥=𝑥𝑖}≈𝑔₀^𝑗_𝑖+1−𝑔₀^𝑗_𝑖−12△𝑥+𝒪△𝑥²,𝜕²𝑔₀𝜕𝑥²_{𝑥=𝑥𝑖}≈𝑔₀^𝑗_𝑖+1−2𝑔₀^𝑗_𝑖+𝑔₀^𝑗_𝑖−1△𝑥²+𝒪△𝑥². (29) We approximate the time derivatives by the forward difference approximations 𝜕ℎ₀𝜕𝑡_{𝑡=𝑡𝑗}≈ℎ_0𝑖^𝑗+1−ℎ₀^𝑗_𝑖(),△𝑡+𝒪△𝑡𝜕𝑔₀𝜕𝑡_{𝑡=𝑡𝑗}≈𝑔_0𝑖^𝑗+1−𝑔₀^𝑗_𝑖().△𝑡+𝒪△𝑡 (30) We obtain the forward-time central-space (FTCS) approximation to (28) given by ℎ_0𝑖^𝑗+1=ℎ₀^𝑗_𝑖−△𝑡ℎ2△𝑥^0𝑗𝑖+1−ℎ₀^𝑗_𝑖−1+△𝑡2△𝑥²𝑔₀^𝑗_𝑖+1−2𝑔₀^𝑗_𝑖+𝑔₀^𝑗_𝑖−1,𝑔_0𝑖^𝑗+1=𝑔₀^𝑗_𝑖−△𝑡𝑔4△𝑥^0𝑗𝑖+1−𝑔₀^𝑗_𝑖−1−△𝑡ℎ2△𝑥^0𝑗𝑖+1−ℎ₀^𝑗_𝑖−1+△𝑡△𝑥²𝑔₀^𝑗_𝑖+1−2𝑔₀^𝑗_𝑖+𝑔₀^𝑗_𝑖−1. (31)(32) To perform a von Neumann stability analysis we substitute ℎ₀^𝑗_𝑖=𝐻^𝑗𝑒^{𝐼𝜔𝑖△𝑥},𝑔₀^𝑗_𝑖=𝐺^𝑗𝑒^{𝐼𝜔𝑖△𝑥}, (33) into (32) where 𝐼²=−1 to obtain the system 𝐻^𝑗+1𝐺^𝑗+1𝐻=𝐴^𝑗𝐺^𝑗, (34) where𝐴=1−𝐼△𝑡()△𝑥sin𝜔△𝑥△𝑡△𝑥²(())−cos𝜔△𝑥−1𝐼△𝑡()△𝑥sin𝜔△𝑥1−𝐼△𝑡()+2△𝑥sin𝜔△𝑥2△𝑡△𝑥²(cos(𝜔△𝑥)−1). (35)

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 𝜌(𝐴)<1. (36) The matrix 𝐴 admits the eigenvalues 𝜆₁=1−𝐼△𝑡(),𝜆2△𝑥sin𝜔△𝑥₂=1+2△𝑡△𝑥²(cos(𝜔△𝑥)−1)−𝐼△𝑡().△𝑥sin𝜔△𝑥 (37) The spectral radius, 𝜌(𝐴), is defined as 𝜌(𝐴)=max𝜆_𝑖,𝑖=1,2. (38) We find that 𝜆₁=1+△𝑡4△𝑥²sin²(),𝜆𝜔△𝑥₂=1+2△𝑡△𝑥²(cos(𝜔△𝑥)−1)²+△𝑡△𝑥²sin²(𝜔△𝑥). (39) Condition (36) therefore gives 1+2△𝑡△𝑥²(cos(𝜔△𝑥)−1)²+△𝑡△𝑥²sin²(𝜔△𝑥)<1. (40) The nonlinear equation (40) can be simplified to give the Courant-Friedrichs-Lewy (CFL) condition △𝑡△𝑥²<12, (41)

where we have chosen 𝜔=(𝜋+2𝑧𝜋)/△𝑥 for 𝑧∈𝕀.

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.

4. 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⁻⁶ where we have chosen 𝑥₀=0, 𝑥_𝑛=4 and ℎ=(𝑥_𝑛−𝑥₀)/𝑛. We take our iterations up to a final time of 𝑡=1.0. The mean and maximum differences between the different approximations are defined by the 𝐿₂ and 𝐿_∞ norms where 𝐿₂=𝑢𝑥_𝑖,𝑡_𝑗−𝑢^∗𝑥_𝑖,𝑡_𝑗2=△𝑥^𝑛_𝑖=0𝑢𝑥_𝑖,𝑡_𝑗−𝑢^∗𝑥_𝑖,𝑡_𝑗^21/2,𝐿_∞=𝑢𝑥_𝑖,𝑡_𝑗−𝑢^∗𝑥_𝑖,𝑡_𝑗∞=max_𝑖𝑢𝑥_𝑖,𝑡_𝑗−𝑢^∗𝑥_𝑖,𝑡_𝑗. (42)

The approximations 𝑢(𝑥_𝑖,𝑡_𝑗) and 𝑢^∗(𝑥_𝑖,𝑡_𝑗) correspond to numerical solutions from two different schemes. We tabulate our results in Table 1.

Table 1 

Comparing 𝐿₂ and 𝐿_∞ norms for Roe-TVD, Roe-Sweby-TVD, and Roe-BDF schemes where we have taken 𝑛=500 and Δ𝑡=2.5×10⁻⁶.

From the results in Table 1 we note that the Roe-BDF approach produces the smallest error for the 𝐿₂ 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⁻⁵. This implies that △𝑡/△𝑥²=0.390625 and △𝑡/△𝑥²=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.

Table 2 

Comparing run times in seconds for an endogenous surfactant using OpenMP where we have taken 𝑛=500 and Δ𝑡=2.5×10⁻⁵.

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 method as implemented by Levy and Shearer [2] as well as Peterson and Shearer [1].

(a)

(b)

Figure 1 

Numerical simulation of the evolution of the initial free surface (22) and an endogenous surfactant (20) where we have taken 𝑛=500 and △𝑡=2.5×10−5.

(a)

(b)

Figure 2 

Numerical simulation of the evolution of the initial free surface (22) and an exogenous surfactant (21) where we have taken 𝑛=500 and △𝑡=2.5×10−5.

5. 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 high-order 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 △𝑡/△𝑥²<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.

Acknowledgments

E. Momoniat acknowledges support from the National Research Foundation of South Africa. M. M. Rashidi thanks the Centre for Differential Equations, Continuum Mechanics and Applications at the University of the Witwatersrand, Johannesburg, for their hospitality during his research visit.