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

Authors: E. Momoniat; M. M. Rashidi; R. S. Herbst
Journal: Mathematical Problems in Engineering
 (2013)

---

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


---

## Body

## 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]
(1)ht+∇·(vh)=0,Γt+∇·(wΓ)=δ△Γ,
where
(2)v=12h∇σ(Γ)-β3h2∇h+κ3h2∇△h,w=h∇σ(Γ)-β2h2∇h+κ2h2∇△h.The free surface of the thin film is given byh=h(x,y,t) and the surfactant concentration by Γ=Γ(x,y,t). The nondimensional constants β and κ balance gravity, capillarity, and Marangoni stress, respectively. The nondimensional constant δ=1/Pe where Pe 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 caseh=h(x,t) and Γ=Γ(x,t) where surface tension effects are ignored. The resulting coupled system of nonlinear equations is given by
(3)∂h∂t+∂Q∂x=0,∂Γ∂t+∂P∂x=0,
where
(4)Q=-h22∂Γ∂x+h33,P=-hΓ∂Γ∂x+h22Γ.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,t) will have compact support for t>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 Section2 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 definehij=h(xi,tj) and Γij=Γ(xi,tj) where the spatial domain is discretized into n+1 equidistant intervals of width △x and xi=i△x. The time tj is defined by tj=j△t where △t is the time step length. We approximate the spatial derivatives in (4) by the low-order forward and backward difference approximations
(5)∂Γ∂x|x=xi≈Γi+1j-Γij△x+𝒪(△x),∂Γ∂x|x=xi≈Γij-Γi-1j△x+𝒪(△x),
the central difference approximation
(6)∂Γ∂x|x=xi≈Γi+1j-Γi-1j2△x+𝒪(△x2),
and the high-order central difference approximation
(7)∂Γ∂x|x=xi≈-Γi+2j+8Γi+1j-8Γi-1j+Γi-2j12△x+𝒪(△x4).
We therefore obtain the approximations to the fluxes given by
(8)Qij={-(hij)22(Γi+1j-Γij△x)+(hij)33,i=0,-(hij)22(Γi+1j-Γi-1j2△x)+(hij)33,1≤i≤2,n-2≤i≤n-1,-(hij)22(-Γi+2j+8Γi+1j-8Γi-1j+Γi-2j12△x)+(hij)33,3≤i≤n-3,-(hij)22(Γij-Γi-1j△x)+(hij)33,i=n,Pij={-hijΓij(Γi+1j-Γij△x)+(hij)2Γij2,i=0,-hijΓij(Γi+1j-Γi-1j2△x)+(hij)2Γij2,1≤i≤2,n-2≤i≤n-1,-hijΓij(-Γi+2j+8Γi+1j-8Γi-1j+Γi-2j12△x)+(hij)2Γij2,3≤i≤n-3,-hijΓij(Γij-Γi-1j△x)+(hij)2Γij2,i=n.An Euler forward approximation to the time derivative is given by(9)∂h∂t|t=tj≈hij+1-hij△t+𝒪(△t),∂Γ∂t|t=tj≈Γij+1-Γij△t+𝒪(△t).
A BDF approximation to the time derivative is given by
(10)∂h∂t|t=tj≈3hij+1-4hij+hij-12△t+𝒪(△t2),∂Γ∂t|t=tj≈3Γij+1-4Γij+Γij-12△t+𝒪(△t2).
A third-order TVD Runge-Kutta approximation to the time derivative to the free surface h(x,t) yields
(11)hi(j+1/3)=hij-△t△x(Qi+(1/2)j-Qi-(1/2)j),hi(j+2/3)=34hij+14hi(j+1/3)-14(△t△x)(Qi+(1/2)(j+1/3)-Qi-(1/2)(j+1/3)),hij+1=13hij+23hi(j+2/3)-23(△t△x)(Qi+(1/2)(j+2/3)-Qi-(1/2)(j+2/3)).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
(12)∂h∂t|t=tj+1△x(Qi+(1/2)j-Qi-(1/2)j)=0,∂Γ∂t|t=tj+1△x(Pi+(1/2)j-Pi-(1/2)j)=0.The fluxes are approximated by the three-point central difference schemes with a numerical viscosity term of Roe [29, 30] given by
(13)Qi+(1/2)j=12(Qij+Qi+1j-|ai+(1/2)j|(hi+1j-hij)),Pi+(1/2)j=12(Pij+Pi+1j-|bi+(1/2)j|(Γi+1j-Γij)),
where
(14)ai+(1/2)j={Qi+1j-Qijhi+1j-hij,hi+1j≠hij,hi(Γi+1-Γi△r)-(hij)2,hi+1j=hij,bi+(1/2)j={Pi+1j-PijΓi+1j-Γij,Γi+1j≠Γij,-(hij)22,Γi+1j=Γij.
Using a finite volume approximation to the spatial derivatives in (3) with an Euler approximation to the time derivatives we obtain
(15)hij+1=hij+λ(Qi+(1/2)j-Qi-(1/2)j),Γij+1=Γij+λ(Pi+(1/2)j-Pi-(1/2)j).
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
(16)hij+1=43hij-13hij-1+2λ3(Qi+(1/2)j-Qi-(1/2)j),Γij+1=43Γij-13Γij-1+2λ3(Pi+(1/2)j-Pi-(1/2)j).We implement the multistep method by evaluating (15) for j=0 to determine hi1 and Γi1. We evaluate (16) for j=1,2,…, where we use the values for hi1 and Γi1 obtained from (15) to start the scheme. Both (15) and (16) are evaluated for i=1,2,…,n-1. The values at i=0 and i=n 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
(17)QLWi+(1/2)j=12(Qij+Qi+1j-△t△x(ai+(1/2)j)2(hi+1j-hij)),PLWi+(1/2)j=12(Pij+Pi+1j-△t△x(bi+(1/2)j)2(Γi+1j-Γij)),
such that
(18)Q*i+(1/2)j=Qi+(1/2)j+Φ(r)(QLWi+(1/2)j-Qi+(1/2)j),P*i+(1/2)j=Pi+(1/2)j+Φ(w)(PLWi+(1/2)j-Pi+(1/2)j),r=hi-hi-1hi+1-hi,w=Γi-Γi-1Γi+1-Γi.
After some numerical experimentation we find that the limiter that works best is the minmod limiter given by
(19)Φ(r)=max[0,min(1,r)].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(20)Γ(x,0)=1,
and for a deposited or exogenous surfactant [17, 26, 27]
(21)Γ(x,0)=exp(-x2).
We consider an initial exponential film profile given by
(22)h(x,0)=exp(-x2).
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
(23)h(0,t)=1,Γ(0,t)=1.
The continuity boundary conditions
(24)∂h∂x|x=0,x=∞=0,∂Γ∂x|x=0,x=∞=0
close the problem. The boundary conditions (23) coupled with (24) give
(25)h0j=1,h1j=1,Γ0j=1,Γ1j=1.
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 x=∞ 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 x=∞ gives
(26)hnj+1=hnj,hn-1j+1=hn-1j,Γnj+1=Γnj,Γn-1j+1=Γn-1j.


## 3. Stability

In this section we consider the stability of the system (3). We linearise the system around the constant solutions h=1 and Γ=1 by making the substitutions
(27)h=1+ϵh0(x,t),Γ=1+ϵΓ0(x,t),
where ϵ≪1. Substituting (27) into (3) and separating to leading order coefficients in ϵ we obtain the linear system
(28)∂h0∂t+∂h0∂x-12∂2g0∂x2=0,∂g0∂t+12∂g0∂x+∂h0∂x-∂2g0∂x2=0.
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
(29)∂h0∂x|x=xi≈h0i+1j-h0i-1j2△x+𝒪(△x2),∂g0∂x|x=xi≈g0i+1j-g0i-1j2△x+𝒪(△x2),∂2g0∂x2|x=xi≈g0i+1j-2g0ij+g0i-1j△x2+𝒪(△x2).
We approximate the time derivatives by the forward difference approximations
(30)∂h0∂t|t=tj≈h0ij+1-h0ij△t+𝒪(△t),∂g0∂t|t=tj≈g0ij+1-g0ij△t+𝒪(△t).
We obtain the forward-time central-space (FTCS) approximation to (28) given by
(31)h0ij+1=h0ij-△t2△x(h0i+1j-h0i-1j)+△t2△x2(g0i+1j-2g0ij+g0i-1j),(32)g0ij+1=g0ij-△t4△x(g0i+1j-g0i-1j)-△t2△x(h0i+1j-h0i-1j)+△t△x2(g0i+1j-2g0ij+g0i-1j).
To perform a von Neumann stability analysis we substitute
(33)h0ij=HjeIωi△x,g0ij=GjeIωi△x,
into (32) where I2=-1 to obtain the system
(34)(Hj+1Gj+1)=A(HjGj),
where(35)A=(1-I△t△xsin(ω△x)△t△x2(cos(ω△x)-1)-I△t△xsin(ω△x)1-I△t2△xsin(ω△x)+2△t△x2(cos(ω△x)-1)).In terms of a von Neumann stability analysis we need to show that the amplification factor|Hj+1/Hj|<1 and |Gj+1/Gj|<1. Alternatively, we need to show that the iterative system (34) is bounded. We can show this by showing that the spectral radius ρ(A) satisfies the condition
(36)ρ(A)<1.
The matrix A admits the eigenvalues
(37)λ1=1-I△t2△xsin(ω△x),λ2=1+2△t△x2(cos(ω△x)-1)-I△t△xsin(ω△x).
The spectral radius, ρ(A), is defined as
(38)ρ(A)=max(|λi|),i=1,2.
We find that
(39)|λ1|=1+(△t4△x)2sin2(ω△x),|λ2|=(1+2△t△x2(cos(ω△x)-1))2+(△t△x)2sin2(ω△x).
Condition (36) therefore gives
(40)(1+2△t△x2(cos(ω△x)-1))2+(△t△x)2sin2(ω△x)<1.
The nonlinear equation (40) can be simplified to give the Courant-Friedrichs-Lewy (CFL) condition
(41)△t△x2<12,where we have chosenω=(π+2zπ)/△x for z∈𝕀.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 theMinGW version of C++. We first consider a simulation of all three numerical schemes taking n=250 and △t=2.5×10-6 where we have chosen x0=0, xn=4 and h=(xn-x0)/n. We take our iterations up to a final time of t=1.0. The mean and maximum differences between the different approximations are defined by the L2 and L∞ norms where
(42)L2=∥u(xi,tj)-u*(xi,tj)∥2=(△x∑i=0n‍|u(xi,tj)-u*(xi,tj)|2)1/2,L∞=∥u(xi,tj)-u*(xi,tj)∥∞=maxi|u(xi,tj)-u*(xi,tj)|.The approximationsu(xi,tj) and u*(xi,tj) correspond to numerical solutions from two different schemes. We tabulate our results in Table 1.Table 1
ComparingL2 and L∞ norms for  Roe-TVD,  Roe-Sweby-TVD, and  Roe-BDF schemes where we have taken n=500 and Δt=2.5×10-6.
Method
h
(
x
,
t
)
g
(
x
,
t
)
L
2
L
∞
L
2
L
∞
Endogenous surfactant
Roe-TVD versus  Roe-Sweby-TVD
0.00694248
0.01257900
0.0015539400
0.002015310
Roe-Sweby-TVD versus  Roe-BDF
0.00168707
0.00780482
0.0000845757
0.000414136
Roe-TVD  versus  Roe-BDF
0.00721239
0.01252690
0.0015349300
0.002011100
Exogenous surfactant
Roe-TVD  versus  Roe-Sweby-TVD
0.00919323
0.01143190
0.0111612
0.03605900
Roe-Sweby-TVD  versus  Roe-BDF
0.00636653
0.01057280
0.0021225
0.00227495
Roe-TVD  versus  Roe-BDF
0.00571950
0.00672977
0.0105552
0.03629150From the results in Table1 we note that the Roe-BDF approach produces the smallest error for the L2 and L∞ 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 taken=500 and △t=2.5×10-5. This implies that △t/△x2=0.390625 and △t/△x2=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 usingOpenMP where we have taken n=500 and Δt=2.5×10-5.
Final time
Nonparallel
Parallel
1
48 s
43 s
2
96 s
86 s
5
237 s
227 s
10
457 s
431 sFrom Table2 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 n values) and long-time solutions, a significant improvement is noted.In Figures1 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].Numerical simulation of the evolution of the initial free surface (22) and an endogenous surfactant (20) where we have taken n=500 and △t=2.5×10-5.
(a)
(b)Numerical simulation of the evolution of the initial free surface (22) and an exogenous surfactant (21) where we have taken n=500 and △t=2.5×10-5.
(a)
(b)


## 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 △t/△x2<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.



---

*Source: 325132-2013-04-24.xml*