Mathematical Problems in Engineering

Volume 2013 (2013), Article ID 325132, 8 pages

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

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

Received 4 February 2013; Accepted 29 March 2013

Academic Editor: Fazal M. Mahomed

Copyright © 2013 E. Momoniat 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

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] where

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

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 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 which is not relevant to our purposes. Peterson and Shearer [28] make analytical progress in solving (1) by investigating similarity solutions for .

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 (1) is degenerate at [1]. Peterson and Shearer [1] have pointed out that the degeneracy for the case implies that if has compact support, then the solution will have compact support for . This holds true for the case and for a linear equation of state 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 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 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

We implement the multistep method by evaluating (15) for to determine and . We evaluate (16) for , where we use the values for and obtained from (15) to start the scheme. Both (15) and (16) are evaluated for . The values at 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.

The surfactant concentration has two forms. For a preexisting or endogenous surfactant and for a deposited or exogenous surfactant [17, 26, 27] 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 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

#### 3. Stability

In this section we consider the stability of the system (3). We linearise the system around the constant solutions and by making the substitutions where . 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 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 to obtain the system where

In terms of a von Neumann stability analysis we need to show that the amplification factor and . 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 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 and where we have chosen , and . We take our iterations up to a final time of . The mean and maximum differences between the different approximations are defined by the and norms where

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

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

#### References

- E. R. Peterson and M. Shearer, “Simulation of spreading surfactant on a thin liquid film,”
*Applied Mathematics and Computation*, vol. 218, no. 9, pp. 5157–5167, 2012. View at Publisher · View at Google Scholar · View at MathSciNet - R. Levy and M. Shearer, “The motion of a thin liquid film driven by surfactant and gravity,”
*SIAM Journal on Applied Mathematics*, vol. 66, no. 5, pp. 1588–1609, 2006. View at Publisher · View at Google Scholar · View at Scopus - R. V. Craster and O. K. Matar, “Dynamics and stability of thin liquid films,”
*Reviews of Modern Physics*, vol. 81, no. 3, pp. 1131–1198, 2009. View at Publisher · View at Google Scholar · View at Scopus - A. B. Afsar-Siddiqui, P. F. Luckham, and O. K. Matar, “The spreading of surfactant solutions on thin liquid films,”
*Advances in Colloid and Interface Science*, vol. 106, no. 1–3, pp. 183–236, 2003. View at Publisher · View at Google Scholar · View at Scopus - O. E. Jensen, “Self-similar, surfactant-driven flows,”
*Physics of Fluids*, vol. 6, no. 3, pp. 1084–1094, 1994. View at Google Scholar · View at Scopus - M. Renardy, “On an equation describing the spreading of surfactants on thin films,”
*Nonlinear Analysis, Theory, Methods and Applications*, vol. 26, no. 7, pp. 1207–1219, 1996. View at Publisher · View at Google Scholar · View at Scopus - M. Renardy, “A singularly perturbed problem related to surfactant spreading on thin films,”
*Nonlinear Analysis, Theory, Methods and Applications*, vol. 27, no. 3, pp. 287–296, 1996. View at Publisher · View at Google Scholar · View at Scopus - B. J. Fischer and S. M. Troian, “Thinning and disturbance growth in liquid films mobilized by continuous surfactant delivery,”
*Physics of Fluids*, vol. 15, no. 12, pp. 3837–3845, 2003. View at Publisher · View at Google Scholar · View at Scopus - O. K. Matar and S. M. Troian, “Spreading of a surfactant monolayer on a thin liquid film: onset and evolution of digitated structures,”
*Chaos*, vol. 9, no. 1, pp. 141–153, 1999. View at Google Scholar · View at Scopus - M. R. E. Warner, R. V. Craster, and O. K. Matar, “Dewetting of ultrathin surfactant-covered films,”
*Physics of Fluids*, vol. 14, no. 11, pp. 4040–4054, 2002. View at Publisher · View at Google Scholar · View at Scopus - K. Kargupta and A. Sharma, “Morphological self-organization by dewetting in thin films on chemically patterned substrates,”
*Journal of Chemical Physics*, vol. 116, no. 7, pp. 3042–3051, 2002. View at Publisher · View at Google Scholar · View at Scopus - K. Kargupta and A. Sharma, “Creation of ordered patterns by dewetting of thin films on homogeneous and heterogeneous substrates,”
*Journal of Colloid and Interface Science*, vol. 245, no. 1, pp. 99–115, 2002. View at Publisher · View at Google Scholar · View at Scopus - D. Halpern, O. E. Jensen, and J. B. Grotberg, “A theoretical study of surfactant and liquid delivery into the lung,”
*Journal of Applied Physiology*, vol. 85, no. 1, pp. 333–352, 1998. View at Google Scholar · View at Scopus - M. S. Borgas and J. B. Grotberg, “Monolayer flow on a thin film,”
*Journal of Fluid Mechanics*, vol. 193, pp. 151–170, 1988. View at Publisher · View at Google Scholar - D. P. Gaver and J. B. Grotberg, “The dynamics of a localized surfactant on a thin film,”
*Journal of Fluid Mechanics*, vol. 213, pp. 127–148, 1990. View at Google Scholar · View at Scopus - D. P. Gaver III and J. B. Grotberg, “Droplet spreading on a thin viscous film,”
*Journal of Fluid Mechanics*, vol. 235, pp. 399–414, 1992. View at Google Scholar - O. E. Jensen and J. B. Grotberg, “Insoluble surfactant spreading on a thin viscous film: shock evolution and film rupture,”
*Journal of Fluid Mechanics*, vol. 240, pp. 259–288, 1992. View at Google Scholar · View at Scopus - O. E. Jensen and J. B. Grotberg, “The spreading of heat or soluble surfactant along a thin liquid film,”
*Physics of Fluids A*, vol. 5, no. 1, pp. 58–68, 1993. View at Google Scholar · View at Scopus - O. E. Jensen, D. Halpern, and J. B. Grotberg, “Transport of a passive solute by surfactant-driven flows,”
*Chemical Engineering Science*, vol. 49, no. 8, pp. 1107–1117, 1994. View at Google Scholar · View at Scopus - D. Halpern and J. B. Grotberg, “Dynamics and transport of a localised soluble surfactant on a thin film,”
*Journal of Fluid Mechanics*, vol. 237, pp. 1–11, 1992. View at Google Scholar - Y. L. Zhang, O. K. Matar, and R. V. Craster, “A theoretical study of chemical delivery within the lung using exogenous surfactant,”
*Medical Engineering and Physics*, vol. 25, no. 2, pp. 115–132, 2003. View at Publisher · View at Google Scholar · View at Scopus - R. V. Craster and O. K. Matar, “Surfactant transport on mucus films,”
*Journal of Fluid Mechanics*, vol. 425, pp. 235–258, 2000. View at Google Scholar - O. K. Matar, R. V. Craster, and M. R. E. Warner, “Surfactant transport on highly viscous surface films,”
*Journal of Fluid Mechanics*, vol. 466, pp. 85–111, 2002. View at Publisher · View at Google Scholar · View at Scopus - Y. L. Zhang, O. K. Matar, and R. V. Craster, “Analysis of tear film rupture: effect of non-Newtonian rheology,”
*Journal of Colloid and Interface Science*, vol. 262, no. 1, pp. 130–148, 2003. View at Publisher · View at Google Scholar · View at Scopus - L. W. Schwartz, R. V. Roy, R. R. Eley, and H. M. Princen, “Surfactant-driven motion and splitting of droplets on a substrate,”
*Journal of Engineering Mathematics*, vol. 50, no. 2-3, pp. 157–175, 2004. View at Publisher · View at Google Scholar · View at Scopus - B. D. Edmonstone and O. K. Matar, “Flow of surfactant-laden thin films down an inclined plane,”
*Journal of Engineering Mathematics*, vol. 50, pp. 141–156, 2004. View at Google Scholar - H. A. Stone, “A simple derivation of the time-dependent convective-diffusion equation for surfactant transport along a deforming interface,”
*Physics of Fluids A*, vol. 2, no. 1, pp. 111–112, 1990. View at Google Scholar · View at Scopus - E. R. Peterson and M. Shearer, “Radial spreading of a surfactant on a thin liquid film,”
*Applied Mathematics Research eXpress*, vol. 2011, no. 1, pp. 1–22, 2011. View at Publisher · View at Google Scholar · View at Scopus - P. L. Roe, “The use of the Riemann problem in finite-difference schemes,” in
*Proceedings of the 7th International Conference on Numerical Methods in Fluid Dynamics*, vol. 141 of*Lecture Notes in Physics*, pp. 354–359, Springer, New York, NY, USA, 1980. - H. C. Yee, “Upwind and symmetric shock-capturing schemes,”
*NASA Technical Memorandum*89464, 1987. View at Google Scholar - S. Jaisankar and S. V. Raghurama Rao, “A central Rankine-Hugoniot solver for hyperbolic conservation laws,”
*Journal of Computational Physics*, vol. 228, no. 3, pp. 770–798, 2009. View at Publisher · View at Google Scholar · View at Scopus - H. Q. Yang and A. J. Przekwas, “A comparative study of advanced shock-capturing schemes applied to Burgers' equation,”
*Journal of Computational Physics*, vol. 102, no. 1, pp. 139–159, 1992. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - C.-W. Shu, “Total-variation-diminishing time discretizations,”
*Journal on Scientific and Statistical Computing*, vol. 9, no. 6, pp. 1073–1084, 1988. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - S. Gottlieb and C. W. Shu, “Total variation diminishing Runge-Kutta schemes,”
*Mathematics of Computation*, vol. 67, no. 221, pp. 73–85, 1998. View at Google Scholar · View at Scopus - A. Iserles,
*A First Course in the Numerical Analysis of Differential Equations*, Cambridge University Press, Cambridge, UK, 1996. View at MathSciNet - T. G. Myers, “Thin films with high surface tension,”
*SIAM Review*, vol. 40, no. 3, pp. 441–462, 1998. View at Google Scholar · View at Scopus - J. A. Diez, L. Kondic, and A. Bertozzi, “Global models for moving contact lines,”
*Physical Review E*, vol. 63, no. 1, Article ID 011208, 13 pages, 2000. View at Publisher · View at Google Scholar · View at Scopus