Marangoni effects on a thin liquid film coating a sphere with axial or   radial thermal gradients

Di Kang; Ali Nadim; Marina Chugunova

arXiv:1703.10471·cond-mat.soft·August 2, 2017

Marangoni effects on a thin liquid film coating a sphere with axial or radial thermal gradients

Di Kang, Ali Nadim, Marina Chugunova

PDF

TL;DR

This paper models the dynamics of a thin liquid film on a sphere under gravity, surface tension, and thermal gradients, analyzing stability and flow patterns influenced by Marangoni effects through numerical simulations.

Contribution

It derives a nonlinear PDE incorporating Marangoni effects for a spherical thin film and analyzes stability under axial and radial thermal gradients with numerical simulations.

Findings

01

Steady states depend on the balance of gravity and Marangoni effects.

02

Axial gradients lead to uniform or localized film accumulation.

03

Radial gradients induce unstable non-axisymmetric modes.

Abstract

We study the time evolution of a thin liquid film coating the outer surface of a sphere in the presence of gravity, surface tension and thermal gradients. We derive the fourth-order nonlinear partial differential equation that models the thin film dynamics, including Marangoni terms arising from the dependence of surface tension on temperature. We consider two different imposed temperature distributions with axial or radial thermal gradients. We analyze the stability of a uniform coating under small perturbations and carry out numerical simulations in COMSOL for a range of parameter values. In the case of an axial temperature gradient, we find steady states with either uniform film thickness, or with the fluid accumulating at the bottom or near the top of the sphere, depending on the total volume of liquid in the film, dictating whether gravity or Marangoni effects dominate. In the case…

Tables2

Table 1. Table 1: Growth rate λ num subscript 𝜆 num \lambda_{\mbox{\tiny num}} calculated for the axisymmetric initial condition ( l = 3 𝑙 3 l=3 , m = 0 𝑚 0 m=0 ) at different times with parameters 𝒮 = 1 𝒮 1 {\cal S}=1 and 𝒩 = 20 𝒩 20 {\cal N}=20 .

$t$	$\max_{x} h (x, t)$	$λ_{num}$
$10^{- 3}$	$1.1277$	$120.18$
$2 \times 10^{- 3}$	$1.2722$	$120.37$
$4 \times 10^{- 3}$	$1.6212$	$121.06$

Table 2. Table 2: Growth rate λ 𝜆 \lambda calculated from the numerical results. The pair of values for each case include the value based on the L ∞ superscript 𝐿 L^{\infty} norm on the left and that based on the L 1 superscript 𝐿 1 L^{1} norm on the right. The numerical tolerance was set to 10 − 7 superscript 10 7 10^{-7} in COMSOL, and the maximum time step was taken to be Δ t = 10 − 6 . Δ 𝑡 superscript 10 6 \Delta t=10^{-6}.

	$t = 10^{- 5}$	$t = 10^{- 4}$	$t = 5 \times 10^{- 4}$
$m = 0$	$120.93$ , $120.16$	$120.47$ , $120.36$	$121.74$ , $121.73$
$m = 1$	$121.81$ , $136, 71$	$133.53$ , $126.20$	$123.20$ , $107.62$
$m = 2$	$139.82$ , $136.07$	$142.50$ , $134.93$	$142.48$ , $134.67$
$m = 3$	$129.90$ , $124.05$	$129.35$ , $124.54$	$132.51$ , $125.33$

Equations123

F (r, θ, ϕ, t) = r - R - h (θ, ϕ, t) = 0 .

F (r, θ, ϕ, t) = r - R - h (θ, ϕ, t) = 0 .

\frac{\partial h}{\partial t} = v_{r} - \frac{v _{θ}}{r} \frac{\partial h}{\partial θ} - \frac{v _{ϕ}}{r sin θ} \frac{\partial h}{\partial ϕ},

\frac{\partial h}{\partial t} = v_{r} - \frac{v _{θ}}{r} \frac{\partial h}{\partial θ} - \frac{v _{ϕ}}{r sin θ} \frac{\partial h}{\partial ϕ},

\frac{1}{r ^{2}} \frac{\partial}{\partial r} (r^{2} v_{r}) + \frac{1}{r sin θ} \frac{\partial}{\partial θ} (sin θ v_{θ}) + \frac{1}{r sin θ} \frac{\partial}{\partial ϕ} v_{ϕ} = 0 .

\frac{1}{r ^{2}} \frac{\partial}{\partial r} (r^{2} v_{r}) + \frac{1}{r sin θ} \frac{\partial}{\partial θ} (sin θ v_{θ}) + \frac{1}{r sin θ} \frac{\partial}{\partial ϕ} v_{ϕ} = 0 .

sin θ (R + h)^{2} v_{r} ∣_{R + h}

sin θ (R + h)^{2} v_{r} ∣_{R + h}

- sin θ (R + h) v_{θ} ∣_{R + h} \frac{\partial h}{\partial θ} - (R + h) v_{ϕ} ∣_{R + h} \frac{\partial h}{\partial ϕ} = 0,

sin θ (R + h)^{2} \frac{\partial h}{\partial t} + \frac{\partial}{\partial θ} (sin θ \int_{R}^{R + h} r v_{θ} d r) + \frac{\partial}{\partial ϕ} (\int_{R}^{R + h} r v_{ϕ} d r) = 0 .

sin θ (R + h)^{2} \frac{\partial h}{\partial t} + \frac{\partial}{\partial θ} (sin θ \int_{R}^{R + h} r v_{θ} d r) + \frac{\partial}{\partial ϕ} (\int_{R}^{R + h} r v_{ϕ} d r) = 0 .

P = p - ρ g \cdot x .

P = p - ρ g \cdot x .

\frac{\partial P}{\partial r}

\frac{\partial P}{\partial r}

\frac{1}{r} \frac{\partial P}{\partial θ}

\frac{1}{r sin θ} \frac{\partial P}{\partial ϕ}

\nabla \cdot \hat{n} \approx \frac{2}{R} - \frac{1}{R ^{2}} (2 h + \frac{1}{sin θ} \frac{\partial}{\partial θ} (sin θ \frac{\partial h}{\partial θ}) + \frac{1}{sin ^{2} θ} \frac{\partial ^{2} h}{\partial ϕ ^{2}}) .

\nabla \cdot \hat{n} \approx \frac{2}{R} - \frac{1}{R ^{2}} (2 h + \frac{1}{sin θ} \frac{\partial}{\partial θ} (sin θ \frac{\partial h}{\partial θ}) + \frac{1}{sin ^{2} θ} \frac{\partial ^{2} h}{\partial ϕ ^{2}}) .

P (θ, ϕ, t) = p_{o} + σ \nabla \cdot \hat{n} + ρ g (R + h) cos (θ) .

P (θ, ϕ, t) = p_{o} + σ \nabla \cdot \hat{n} + ρ g (R + h) cos (θ) .

\hat{n} \cdot (τ_{2} - τ_{1}) \cdot \hat{t} = \nabla_{s} σ,

\hat{n} \cdot (τ_{2} - τ_{1}) \cdot \hat{t} = \nabla_{s} σ,

σ (T) = σ (T_{0}) + \frac{\partial σ}{\partial T}_{T = T_{0}} (T - T_{0}) + \dots = σ_{0} + σ_{1} (T - T_{0}) + \dots

σ (T) = σ (T_{0}) + \frac{\partial σ}{\partial T}_{T = T_{0}} (T - T_{0}) + \dots = σ_{0} + σ_{1} (T - T_{0}) + \dots

μ r \frac{\partial}{\partial r} (\frac{v _{θ}}{r}) = \frac{1}{r} \frac{\partial σ}{\partial θ} \mbox a t r = R + h,

μ r \frac{\partial}{\partial r} (\frac{v _{θ}}{r}) = \frac{1}{r} \frac{\partial σ}{\partial θ} \mbox a t r = R + h,

μ r \frac{\partial}{\partial r} (\frac{v _{ϕ}}{r}) = \frac{1}{r sin θ} \frac{\partial σ}{\partial ϕ} \mbox a t r = R + h .

μ r \frac{\partial}{\partial r} (\frac{v _{ϕ}}{r}) = \frac{1}{r sin θ} \frac{\partial σ}{\partial ϕ} \mbox a t r = R + h .

v_{θ} = \frac{1}{2 μ} \frac{\partial P}{\partial θ} r - \frac{C _{1}}{r} + C_{2},

v_{θ} = \frac{1}{2 μ} \frac{\partial P}{\partial θ} r - \frac{C _{1}}{r} + C_{2},

v_{ϕ} = \frac{1}{2 μ sin θ} \frac{\partial P}{\partial ϕ} r - \frac{C _{3}}{r} + C_{4} .

v_{ϕ} = \frac{1}{2 μ sin θ} \frac{\partial P}{\partial ϕ} r - \frac{C _{3}}{r} + C_{4} .

C_{1}

C_{1}

C_{2}

C_{3}

C_{4}

\frac{\partial h}{\partial t} - \frac{1}{R ^{2} sin θ} \frac{\partial}{\partial θ} (\frac{h ^{3} sin θ}{3 μ} \frac{\partial P}{\partial θ} - \frac{h ^{2} sin θ}{2 μ} \frac{\partial σ}{\partial θ}) - \frac{1}{R ^{2} sin ^{2} θ} \frac{\partial}{\partial ϕ} (\frac{h ^{3}}{3 μ} \frac{\partial P}{\partial ϕ} - \frac{h ^{2}}{2 μ} \frac{\partial σ}{\partial ϕ}) = 0,

\frac{\partial h}{\partial t} - \frac{1}{R ^{2} sin θ} \frac{\partial}{\partial θ} (\frac{h ^{3} sin θ}{3 μ} \frac{\partial P}{\partial θ} - \frac{h ^{2} sin θ}{2 μ} \frac{\partial σ}{\partial θ}) - \frac{1}{R ^{2} sin ^{2} θ} \frac{\partial}{\partial ϕ} (\frac{h ^{3}}{3 μ} \frac{\partial P}{\partial ϕ} - \frac{h ^{2}}{2 μ} \frac{\partial σ}{\partial ϕ}) = 0,

P (θ, ϕ, t) = ρ g R cos θ - \frac{σ}{R ^{2}} (2 h + \frac{1}{sin θ} \frac{\partial}{\partial θ} (sin θ \frac{\partial h}{\partial θ}) + \frac{1}{sin ^{2} θ} \frac{\partial ^{2} h}{\partial ϕ ^{2}}) .

P (θ, ϕ, t) = ρ g R cos θ - \frac{σ}{R ^{2}} (2 h + \frac{1}{sin θ} \frac{\partial}{\partial θ} (sin θ \frac{\partial h}{\partial θ}) + \frac{1}{sin ^{2} θ} \frac{\partial ^{2} h}{\partial ϕ ^{2}}) .

T = T_{0} + k_{1} z = T_{0} + k_{1} r cos θ .

T = T_{0} + k_{1} z = T_{0} + k_{1} r cos θ .

σ = σ_{0} + σ_{1} k_{1} (R + h (θ, ϕ, t)) cos θ \approx σ_{0} + σ_{1} k_{1} R cos θ,

σ = σ_{0} + σ_{1} k_{1} (R + h (θ, ϕ, t)) cos θ \approx σ_{0} + σ_{1} k_{1} R cos θ,

G \equiv \frac{ρ g H ^{2} τ}{3 μ R}, S \equiv \frac{σ _{0} H ^{3} τ}{3 μ R ^{4}}, \mbox an d M \equiv \frac{σ _{1} k _{1} H τ}{2 μ R} .

G \equiv \frac{ρ g H ^{2} τ}{3 μ R}, S \equiv \frac{σ _{0} H ^{3} τ}{3 μ R ^{4}}, \mbox an d M \equiv \frac{σ _{1} k _{1} H τ}{2 μ R} .

σ_{1} k_{1} \sim ϵ^{2} \frac{σ _{0}}{R} \sim ϵ ρ g R .

σ_{1} k_{1} \sim ϵ^{2} \frac{σ _{0}}{R} \sim ϵ ρ g R .

\frac{\partial h}{\partial t} + \frac{1}{sin θ} \frac{\partial}{\partial θ} {sin θ h^{3} [G sin θ + S \frac{\partial}{\partial θ} (2 h + \frac{1}{sin θ} \frac{\partial}{\partial θ} (sin θ \frac{\partial h}{\partial θ}) + \frac{1}{sin ^{2} θ} \frac{\partial ^{2} h}{\partial ϕ ^{2}})] - M sin^{2} θ h^{2}} + \frac{1}{sin ^{2} θ} \frac{\partial}{\partial ϕ} {S h^{3} \frac{\partial}{\partial ϕ} (2 h + \frac{1}{sin θ} \frac{\partial}{\partial θ} (sin θ \frac{\partial h}{\partial θ}) + \frac{1}{sin ^{2} θ} \frac{\partial ^{2} h}{\partial ϕ ^{2}})} = 0.

\frac{\partial h}{\partial t} + \frac{1}{sin θ} \frac{\partial}{\partial θ} {sin θ h^{3} [G sin θ + S \frac{\partial}{\partial θ} (2 h + \frac{1}{sin θ} \frac{\partial}{\partial θ} (sin θ \frac{\partial h}{\partial θ}) + \frac{1}{sin ^{2} θ} \frac{\partial ^{2} h}{\partial ϕ ^{2}})] - M sin^{2} θ h^{2}} + \frac{1}{sin ^{2} θ} \frac{\partial}{\partial ϕ} {S h^{3} \frac{\partial}{\partial ϕ} (2 h + \frac{1}{sin θ} \frac{\partial}{\partial θ} (sin θ \frac{\partial h}{\partial θ}) + \frac{1}{sin ^{2} θ} \frac{\partial ^{2} h}{\partial ϕ ^{2}})} = 0.

\frac{\partial h}{\partial t} + \frac{\partial}{\partial x} {(1 - x^{2}) h^{3} [G + S \frac{\partial}{\partial x} (2 h + \frac{\partial}{\partial x} ((1 - x^{2}) \frac{\partial h}{\partial x}) + \frac{1}{1 - x ^{2}} \frac{\partial ^{2} h}{\partial ϕ ^{2}})] - M (1 - x^{2}) h^{2}} + \frac{1}{1 - x ^{2}} \frac{\partial}{\partial ϕ} {S h^{3} \frac{\partial}{\partial ϕ} [2 h + \frac{\partial}{\partial x} ((1 - x^{2}) \frac{\partial h}{\partial x}) + \frac{1}{1 - x ^{2}} \frac{\partial ^{2} h}{\partial ϕ ^{2}}]} = 0,

\frac{\partial h}{\partial t} + \frac{\partial}{\partial x} {(1 - x^{2}) h^{3} [G + S \frac{\partial}{\partial x} (2 h + \frac{\partial}{\partial x} ((1 - x^{2}) \frac{\partial h}{\partial x}) + \frac{1}{1 - x ^{2}} \frac{\partial ^{2} h}{\partial ϕ ^{2}})] - M (1 - x^{2}) h^{2}} + \frac{1}{1 - x ^{2}} \frac{\partial}{\partial ϕ} {S h^{3} \frac{\partial}{\partial ϕ} [2 h + \frac{\partial}{\partial x} ((1 - x^{2}) \frac{\partial h}{\partial x}) + \frac{1}{1 - x ^{2}} \frac{\partial ^{2} h}{\partial ϕ ^{2}}]} = 0,

\frac{\partial ^{2} v _{θ}}{\partial y ^{2}} = \frac{ϵ ^{2} R}{μ} \frac{\partial P}{\partial θ}

\frac{\partial ^{2} v _{θ}}{\partial y ^{2}} = \frac{ϵ ^{2} R}{μ} \frac{\partial P}{\partial θ}

v_{θ} = 0 at y = 0,

v_{θ} = 0 at y = 0,

\frac{\partial v _{θ}}{\partial y} = \frac{ϵ}{μ} \frac{\partial σ}{\partial θ} at y = 1 .

\frac{\partial v _{θ}}{\partial y} = \frac{ϵ}{μ} \frac{\partial σ}{\partial θ} at y = 1 .

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

Marangoni effects on a thin liquid film coating a sphere with axial or radial thermal gradients

Di Kang

[email protected]

Claremont Graduate University

Ali Nadim

[email protected]

Claremont Graduate University

Marina Chugunova

[email protected]

Claremont Graduate University

Abstract

We study the time evolution of a thin liquid film coating the outer surface of a sphere in the presence of gravity, surface tension and thermal gradients. We derive the fourth-order nonlinear partial differential equation that models the thin film dynamics, including Marangoni terms arising from the dependence of surface tension on temperature. We consider two different imposed temperature distributions with axial or radial thermal gradients. We analyze the stability of a uniform coating under small perturbations and carry out numerical simulations in COMSOL for a range of parameter values. In the case of an axial temperature gradient, we find steady states with either uniform film thickness, or with the fluid accumulating at the bottom or near the top of the sphere, depending on the total volume of liquid in the film, dictating whether gravity or Marangoni effects dominate. In the case of a radial temperature gradient, a stability analysis reveals the most unstable non-axisymmetric modes on an initially uniform coating film.

pacs:

TBD

I Introduction

The effect of heating on the dynamics of thin liquid films has been studied by many authors in planar geometry on both smooth and non-smooth substrates. For example an experimental study applying a local heat source to a liquid film falling down a vertical plate was carried out by kabov2002heat They observed the deformation of the film surface into a horizontal bump shape, which became unstable above some critical heating value and eventually gave rise to vertical downstream rivulets. kalliadasis2003marangoni used the integral boundary-layer approximation of the Navier-Stokes system and free-surface boundary conditions to model thin film coating flow on an inclined heated plane and analyzed the linear stability of steady states for different governing dimensionless groups. The stability analysis of a thin liquid film flowing down an inclined heated plane was extended to account for the effects of evaporation or condensation by bankoff . It was also shown that the characteristic time for the development of thin spots decreases as the fourth power of the evaporation rate. Results on the stability of dry patches forming in thin liquid films flowing over a heated planar surface were described by zuber1966stability . It was shown that for liquid metals (high wettability materials), the thermal effects were dominant in determining the stability of dry patches. The three-dimensional evolution of the long-wave Marangoni instability of thin liquid films on a planar surface was studied by oron2000nonlinear . It was seen that the evolving two-dimensional film is unstable to small three-dimensional random perturbations. Large scale surface deformations of a liquid film were studied by bestehorn20033d . Based on a linear stability analysis, they concluded that the shape deformation is unstable due to the Marangoni effect caused by external heating on a smooth and solid substrate. The stability of a heated liquid film on a plane in the case of a temperature dependent fluid viscosity was analyzed by reisfeld1990nonlinear . Through numerical computations, they discovered that when the fluid viscosity is temperature dependent, the rupture time of the film is reduced relative to the constant-viscosity case. jensen1993spreading used lubrication theory to derive evolution equations for a surfactant contaminated liquid film height and cross-sectionally averaged bulk surfactant concentration. They also studied the special case of the governing equations corresponding to the Marangoni flow induced by a locally hot region in the liquid layer. Dynamics of thin liquid films on a heated horizontal wall with parallel grooves on its upper surface were studied experimentally and numerically by alexeev2005marangoni . They found that the heat transfer in liquid thin film depends on the topography of the wall surface.

Fewer results are available for heated thin liquid films in non-planar geometries. Experimental results on film breakdown and heat flux for sub-cooled water films flowing downwards on the outside of a vertical uniformly-heated tube can be found in the work of fujita1978heat . Some numerical simulations of the flow field and associated heat transfer coefficient for a thin liquid film on a horizontal rotating disk were presented by rahman1992numerical . wilson1994onset studied the onset of steady Marangoni convection for a thin liquid film coating a rigid sphere for cases where the free surface was either deformable or non-deformable. It was seen that the liquid layer is always stable when heated from the outside and may be unstable when heated from the inside if the magnitude of the non-dimensional Marangoni number is sufficiently large. In this work, we consider the same geometry as wilson1994onset but rather than solve for the temperature distribution as part of the solution after linearization, we simply assume that the temperature field is externally imposed. The lubrication approximation in the spherical geometry, as used in our own recent work kang2016dynamics but without rotation and with added Marangoni effects, allows us to derive a nonlinear initial-boundary-value problem in the form of a single fourth-order degenerate parabolic partial differential equation for the film thickness. We use this as a basis to study the case where the temperature varies in the vertical direction in the presence of gravity. We then examine the case of radial temperature gradients as well and compare the results from numerical simulations of the nonlinear system with the predictions from a linearized stability analysis.

Other related areas of research where Marangoni effects at a free surface play an important role include thermo-capillary migration of attached droplets along a surface with an imposed temperature gradient ford1994 ; sui2014moving ; evaporative instabilities in thin liquid films and climbing films hosoi2001evaporative ; sultan2005evaporation ; modulating surface stresses for microfluidic actuation darhuber2005principles ; and non-isothermal spreading of liquid drops or thin fluid ridges on solid substrates ehrhard1991non ; dunn2009quasi . The mathematical and computational challenges of solving higher order partial differential equations such as the thin-film equation on complex curved surfacesgreer2006fourth ; azencot2015functional are also topics of active research.

The article has the following structure: in Section II, we describe the physical model and derive the governing nonlinear partial differential equation for the evolution of the film profile under the lubrication approximation. The driving forces in the model include gravity, surface tension and Marangoni effects. In Section III, we consider the case of an axial temperature variation and show that it is possible to find steady states where the downward draining flow due to gravity is exactly balanced by an upward Marangoni flow. Depending on the amount of fluid present, one can obtain a uniform film, or have the accumulation of fluid near the top or bottom of the sphere. The streamfunction for the internal flow is obtained for the case of uniform film thickness, illustrating that while the profile is steady, there is continual flow occurring within the film. In Section IV, we consider the case with a radially varying temperature field. In the absence of gravity, we complete a linear stability analysis starting with a uniform film and identify the set of unstable modes that include both axisymmetric and non-axisymmetric shapes. We also carry out numerical simulations of the evolution of a slightly perturbed initial profile and compare the numerical growth rates of the unstable modes with the predictions from the linearized analysis. For the axisymmetric case, even in the presence of gravity and with Marangoni effects, we find an energy functional for the system with a dissipation property, i.e., with the energy monotonically decreasing during the time evolution.

II Model formulation and governing equations

In this section we derive the model for the time evolution of a thin viscous liquid film on a solid spherical object of radius $R$ in the presence of gravity $\bm{g}$ and with an imposed steady temperature field which will be described later. The derivation parallels the one in our recent paperkang2016dynamics which focused on the effect of centrifugal forces but without Marangoni effects. Here, we do not consider rotational effects but do allow for surface tension variation with temperature and non-axisymmetric shapes. Let $r$ denote the radial distance from the origin, $\theta$ the polar angle ( $\theta=0$ corresponds to the positive $z$ -direction) and $\phi$ the azimuthal angle in spherical coordinates. The imposed temperature field is a function of position given by $T(r,\theta,\phi)$ . The function $h(\theta,\phi,t)$ denotes the thickness of the thin film coating at time $t$ at each point on the spherical substrate. We assume that the model satisfies the standard lubrication approximation constraints. The interface of the thin film $r=R+h(\theta,\phi,t)$ is the zero-level set of the function

[TABLE]

The kinematic boundary condition ${D\mathcal{F}}/{Dt}=0$ results in

[TABLE]

at $r=R+h(\theta,\phi,t).$

The continuity equation $\bm{\nabla\cdot v}=0$ takes the form

[TABLE]

Multiplying this equation by $r^{2}\sin\theta$ , integrating along the $r$ -direction from $r=R$ to $r=R+h$ , and applying the Leibniz rule to bring the derivatives outside the integrals, we obtain

[TABLE]

where the no-penetration boundary condition $v_{r}=0$ at the solid surface $r=R$ has also been applied. When combined with the kinematic boundary condition, this yields the time evolution equation for the thin film thickness $h(\theta,\phi,t)$ as

[TABLE]

In order to get the partial differential equation which models the time evolution of the coating thickness $h(\theta,\phi,t)$ , we need to relate the velocity components $v_{\theta}$ and $v_{\phi}$ to $h(\theta,\phi,t)$ and its derivatives. This can be done by integrating the components of the momentum equation under the Reynolds lubrication approximation.

Let us introduce the modified pressure $P$ as follows:

[TABLE]

Here $\bm{x}$ is the position vector and $\bm{g}$ is the gravitational acceleration which is taken to act in the negative $z$ -direction. In terms of the modified pressure, under the lubrication approximation, the $r$ -, $\theta$ - and $\phi$ -components of the momentum equation simplify tokang2016dynamics

[TABLE]

Eq. (7) implies that the modified pressure is independent of $r$ and may only depend on $\theta$ , $\phi$ and $t$ . The normal stress balance at $r=R+h$ , requires the pressure in the film to be equal to the uniform pressure $p_{o}$ in the external air plus the capillary contribution associated with the surface tension $\sigma$ of the interface and its curvature $\nabla\cdot\hat{{\bf n}}$ , i.e., $p|_{r=R+h}=p_{o}+\sigma\nabla\cdot\hat{{\bf n}}$ . Here, $\hat{{\bf n}}=\nabla{\cal F}/|\nabla{\cal F}|$ (with ${\cal F}$ given by Eq. (1)) is the unit normal at the interface pointing toward the exterior. When the film thickness $h(\theta,\phi,t)$ is much smaller than the sphere radius $R$ , the curvature can be approximated in a series in $h$ , with the first two terms given by:

[TABLE]

In terms of the modified pressure, the normal stress balance becomes:

[TABLE]

To incorporate Marangoni effects into the analysis, we need to consider the tangential stress condition at the interface. The Marangoni effect can be thought of as an effective tangential stress along an interface due to presence of surface tension gradients. When the dependence of surface tension $\sigma$ on the temperature gives rise to such stresses resulting in a flow, this phenomenon is called thermo-capillary convection. The tangential stress balance at the interface is

[TABLE]

where $\mathbf{\hat{n}}$ and $\mathbf{\hat{t}}$ are the normal and tangential unit vectors at the interface and $\bm{\tau}_{1}$ and $\bm{\tau}_{2}$ are the viscous stresses in the liquid and air phases respectively, the latter being negligible. For small variations in temperature, the dependence of surface tension on temperature can be linearized about a base state in the form

[TABLE]

where $\sigma_{0}\equiv\sigma(T_{0})$ and $\sigma_{1}\equiv\left.\frac{\partial\sigma}{\partial T}\right|_{T_{0}}$ . Typically, surface tension is a decreasing function of temperature, meaning that $\sigma_{1}<0$ .

Under the lubrication approximation, the $\theta$ - and $\phi$ -components of the tangential stress balance at the interface take the forms:

[TABLE]

Noting that the modified pressure $P$ is independent of $r$ , the momentum equations (8) and (9) can be integrated to obtain the following general solutions for $v_{\theta}$ and $v_{\phi}$ :

[TABLE]

To solve for the integration “constants” (these are independent of $r$ but may depend on the other variables), we apply the no-slip boundary condition on the solid surface and the tangential stress balances (13) and (14), which include the Marangoni effects, on the liquid-air interface. These result in

[TABLE]

Having found the relationships between $v_{\theta}$ , $v_{\phi}$ and the film thickness $h(\theta,\phi,t)$ , the integrals within the evolution equation (5) can be evaluated and simplified. In the lubrication limit, the film thickness $h$ is everywhere small compared to the sphere radius $R$ , so that all factors $R+h$ and $R-h$ that appear in the resulting equation can be approximated by $R$ to leading order. This results in the following self-contained fourth-order partial differential equation for the evolution of the film thickness:

[TABLE]

where

[TABLE]

The Marangoni terms in the evolution equation (15) are the contributions to the flux components proportional to ${\partial\sigma}/{\partial\theta}$ and ${\partial\sigma}/{\partial\phi}$ . However, there is a subtle issue with this derivation. In particular, since surface tension $\sigma$ appears in the expression for the pressure $P$ , it may seem that differentiating $P$ with respect to $\theta$ or $\phi$ , as needed in the above evolution equation, would introduce additional Marangoni terms involving the derivatives of $\sigma$ . However, scaling arguments given below will show that in the normal stress balance which leads to the dependence of $P$ on $\sigma$ , one can simply use the leading constant term $\sigma_{0}$ for the surface tension, whereas in the actual Marangoni terms, it will be the next term $\sigma_{1}$ that plays the main role.

For the purpose of scaling the equation in the following sections, we can define the characteristic film thickness $H$ to be the average thickness of the thin film based on its initial profile. One can then define the dimensionless film thickness to be $\hat{h}=h/H$ . The lubrication approximation requires $\epsilon={H}/{R}\ll 1$ . Also, one can choose an arbitrary time scale $\tau$ whose value can be determined later by scaling arguments and define a dimensionless time $\hat{t}=t/\tau$ . Within the thin film equation (15), we have driving forces associated with gravity, surface tension and Marangoni effects. The scaling we shall choose for each of the cases below will be one in which all three of these effects are of similar orders of magnitude, balancing one another. This will be done for two different imposed temperature fields (vertical and radial distributions) in the next two sections.

III Vertical temperature field

III.1 The leading order equation

Let us suppose that an externally imposed temperature field exists which varies linearly in the vertical coordinate $z$ with a vertical temperature gradient $k_{1}$ ; that is,

[TABLE]

From the assumed linear dependence of surface tension on temperature, we thus find

[TABLE]

at the liquid-air interface.

Using this result in the evolution equation (15), scaling $h$ with $H$ and $t$ with $\tau$ as suggested earlier, and simplifying to get the equation for $\partial\hat{h}/\partial\hat{t}$ , results in a dimensionless equation in which the following three dimensionless groups appear respectively in front of the contributions to the flux due to gravity, surface tension and Marangoni effects (in the $\theta$ -derivative):

[TABLE]

Note that the Marangoni parameter $\cal M$ may be negative since $\sigma_{1}$ and $k_{1}$ can have either sign. We seek a model in which all the above mentioned effects are present and of the same orders of magnitude. To balance these three terms, we thus need

[TABLE]

Equivalently, the constant part of surface tension $\sigma_{0}$ must be fairly large compared to gravitational effects ( $\sigma_{0}\sim\rho gR^{2}/\epsilon$ ) in order for capillarity and gravity to both be important. Also the typical variation in surface tension across the length of the film ( $\sigma_{1}k_{1}R$ ) due to temperature gradients, must be quite small compared to the constant part of the surface tension $\sigma_{0}$ , i.e., $\sigma_{1}k_{1}R\sim\epsilon^{2}\sigma_{0}$ . In this limit, all three effects (gravity, capillarity and Marangoni) contribute equally to the film evolution and the dimensionless groups $\cal G$ , $\cal S$ and $\cal M$ introduced above are all of order unity. Time scale $\tau$ can be chosen to make any one of the three parameters equal to one exactly, if desired. Also, in the normal stress balance at the interface, $p|_{r=R+h}=p_{o}+\sigma\nabla\cdot\hat{{\bf n}}$ , one can replace $\sigma$ with $\sigma_{0}$ without affecting the subsequent evolution equation at leading order.

Under this scaling, in the $\phi$ -derivative term of the evolution equation, only the capillarity term due to the normal stress balance contributes at leading order. The scaled evolution equation (from which the hats have been dropped for clarity), takes the form:

[TABLE]

where parameters $\cal G$ , $\cal S$ and $\cal M$ are given by equation (17).

With the change of variable $x=-\cos\theta$ , the equation becomes

[TABLE]

over the range $-1<x<1$ and $t>0$ . It is this form of the equation that is used in the numerical simulations presented below.

III.2 Uniform films and their stability

Equation (18) admits a steady solution with a uniform film thickness $h$ provided that ${\cal G}h^{3}-{\cal M}h^{2}=0$ , or when $h={\cal M}/{\cal G}$ . In terms of dimensional parameters, this corresponds to a film of thickness $3\sigma_{1}k_{1}/(2\rho g)$ provided that $\sigma_{1}k_{1}>0$ . For a given slope $\sigma_{1}$ of the surface tension versus temperature curve and a given vertical temperature gradient $k_{1}$ , only a film of that particular thickness can remain uniform along the entire spherical surface. In that state, the draining of the film due to gravity perfectly balances the upward flow caused by the temperature gradient (and Marangoni effect) to keep the film steady and perfectly uniform. A natural question to ask is what happens if an initial film of that thickness is perturbed slightly in $\theta$ and $\phi$ , or if the volume of the film is such that its initially uniform thickness is either smaller or larger than that particular value.

For convenience, if we focus on the case ${\cal G}-{\cal M}=0$ , the uniform film thickness would be $h=1$ . To check the stability of that state by numerical simulations, in Figure 1 we show the evolution of a slightly perturbed film computed by COMSOL and observe that, without perturbing the overall volume of the film, such a film relaxes back to its original uniform state where the gravity and Marangoni effects are in balance. The uniform steady state is thus stable with respect to small perturbations which do not change the total mass (or volume) of the liquid.

However, that uniform state is unstable to perturbations that change the mass (or average initial thickness) of the liquid film away from its equilibrium state. If the the perturbation increases the mass slightly, the gravitational term dominates and and the uniform shape develops a droplet profile at the bottom of the sphere, as seen in Figure 2, but if the perturbation decreases the mass or the average film thickness, the Marangoni effect wins and the uniform profile develops a bump near the top of the sphere, as seen in Figure 3.

If ${\cal G}-{\cal M}\neq 0$ , then the uniform thickness $h=1$ is not a stationary solution. Starting with a film of average initial thickness 1, if ${\cal G}-{\cal M}>0$ , the combined gravity and Marangoni effects will drive the fluid downward and the steady state shape will be like a hanging drop near the bottom of the sphere, i.e., gravity wins (see Figure 4); whereas if ${\cal G}-{\cal M}<0$ , the combined gravity and Marangoni effect will drive the fluid upward to form a bump near the top of the sphere, i.e., Marangoni wins (see Figure 5).

The smaller the initial perturbation in volume, the smaller would be the deviation of the final stationary shape from the uniform one. This can be seen in Figure 6 where for the case ${\cal G}={\cal M}=10$ that still corresponds to equilibrium thickness $h=1$ , the initial profile and the final steady state profile are plotted when the initial thickness has values 0.99, 0.999, 1.001 and 1.01. As expected, initial films that are thicker than the uniform equilibrium value develop maxima near the bottom of the drop (at $x=1$ ), whereas the initially thinner films develop minima near the bottom. Interestingly, the steady state film thickness is not a monotonically increasing or decreasing function in either case.

III.3 Streamlines for the steady surface flow

For the case $\mathcal{G}=\mathcal{M}$ , The uniform film $h=1$ is a steady state solution. To obtain the streamlines of the flow taking place within this uniform film, it is convenient to transform the radial coordinate to a new variable $y$ (scaled and dimensionless) that ranges from 0 to 1 across the film thickness. Let us choose $y={(r-R)}/{(\epsilon R)}$ . Recall that $\epsilon=H/R$ where $H$ would be the dimensional uniform film thickness. Rewrite Eq. (8) and Eq. (13) in terms of variable $y$ and take the leading order to get

[TABLE]

with boundary conditions

[TABLE]

Note that only variable $y$ is dimensionless, with all other variables and parameters still being dimensional. Solve for $v_{\theta}$ to get

[TABLE]

While it may appear that the two terms in this equation are of different orders, because the Marangoni term is an order of $\epsilon$ smaller than the gravity term, the two terms in the equation are of the same order of magnitude, as will be seen below. When the film is of uniform thickness, from the expression (16) for the pressure field we find that

[TABLE]

Also, from (17), when $\cal G=\cal M$ , we have that:

[TABLE]

Thus the expression for $v_{\theta}$ simplifies to:

[TABLE]

This is a quadratic profile which changes sign at $y=2/3$ , with the fluid flowing downward due to gravity in the bulk region $0<y<2/3$ , but flowing back up due to the Marangoni effect in $2/3<y<1$ near the interface, with zero net flux across the whole cross section. Substituting $v_{\theta}$ into the continuity equation (3), we can solve for $v_{r}$ :

[TABLE]

The streamfunction corresponding to this flow is given by

[TABLE]

In this case, the velocity components are related to the streamfunction by $v_{r}=\frac{\epsilon}{\sin\theta}\frac{\partial\psi}{\partial\theta}$ and $v_{\theta}=\frac{-1}{\sin\theta}\frac{\partial\psi}{\partial y}$ . The streamlines are illustrated in Figure 7.

IV Radial temperature field

IV.1 The leading order equation

In this section, we consider the specific case where the externally imposed temperature field is a linear function of $r$ only:

[TABLE]

Temperature $T_{0}$ is a uniform temperature along the surface of the solid sphere and depending on whether the gradient $k_{2}$ is positive or negative, temperature increases or decreases outward from the surface. As such, the surface tension $\sigma$ of the liquid-air interface at coordinate $(\theta,\phi)$ depends directly on the film thickness $h(\theta,\phi,t)$ . In this sense, we can treat $\sigma$ as a function of $h$ :

[TABLE]

The sign of the product $\sigma_{1}k_{2}$ determines whether surface tension is higher or lower in regions where the film is thinner or thicker.

Using this result for the surface tension in the evolution equation (15), scaling $h$ with $H$ and $t$ with $\tau$ as before and simplifying results in a dimensionless evolution equation in which the following three dimensionless groups appear respectively in front of the contributions to the flux due to gravity, surface tension and Marangoni effects (in the $\theta$ -derivative):

[TABLE]

The first two parameter are the same as before, while the third, measuring Marangoni effects, is a little different. To balance these terms so that all effects are comparable and present at leading order we now need

[TABLE]

which is different from the scaling in the case of a vertical temperature gradient with respect to the Marangoni effect.

Under this scaling, surface tension and Marangoni effects will also appear in the $\phi$ -derivative terms. The leading order evolution equation assumes the form

[TABLE]

After the change of variable $x=-\cos\theta$ and upon neglecting gravity for the time being (i.e., setting ${\cal G}=0$ ) the equation becomes

[TABLE]

In the absence of gravity and with a radially symmetric temperature field, it is expected that a uniform film thickness (e.g., $h(x,\phi,t)=1$ ) should be a solution to the above equation, which is obviously the case. We are interested in the stability of this solution. In order to investigate this, consider a slight perturbation of that solution in the form:

[TABLE]

Substituting into the evolution equation and taking the limit as $\varepsilon\rightarrow 0$ , we obtain the linearized equation (dropping the tilde for convenience):

[TABLE]

Upon assuming a modal decomposition for the small perturbation $h$ in the form

[TABLE]

where $P_{l}^{m}(x)$ is the associated Legendre function, which is a solution of

[TABLE]

the growth or decay rates $\lambda$ are found to be

[TABLE]

Interestingly, the rates are independent of the index $m$ of the azimuthal modes.

The stability of a mode with index $l$ depends upon the sign of ${\cal N}$ . (Note that the dimensionless group $\cal S$ is always positive.) If ${\cal N}<0$ , $\lambda$ is negative for all positive $l$ , which mean that the uniform film thickness is a stable solution. If ${\cal N}>0$ , $\lambda$ may be positive for some small values of $l$ , but will become negative with increasing $l$ ; therefore a certain range of mode numbers may be unstable. An intuitive explanation of the instability is as follows. For ${\cal N}>0$ , surface tension is higher where the film is thicker and lower where the film is thinner. As such, if an initially uniform film is perturbed, the Marangoni effect which drives the surface flow from low to high surface tension regions causes the bulk fluid to be “pumped” away from the thin regions (with low surface tension) towards the thick regions (with high surface tension), thereby causing the thin regions to become thinner and the thick regions to become thicker, amplifying the initial perturbation.

IV.2 Numerical simulations

To check the stability predictions numerically, we choose ${\cal S}=1$ and ${\cal N}=20$ . Also we choose the amplitude of the initial perturbation to be $\varepsilon=10^{-5}$ . The relationship between the eigenvalues $\lambda$ and modes $l$ is plotted in Figure 8 to illustrate that the higher modes are stabilized ( $\lambda<0$ ) due to surface tension $\cal S$ . The most unstable mode for these parameter values corresponds to $l=3$ .

For an axisymmetric initial condition, we choose $l=3$ and $m=0$ . We can simplify the evolution equation when the film thickness $h$ is independent of $\phi$ . We carry out the numerical simulations of the full nonlinear equation in COMSOL Multiphysics. The initial condition is

[TABLE]

For $l=3$ , the predicted growth rate is $\lambda=120$ . We can estimate the growth rate from the numerical simulations using

[TABLE]

The results are shown in Table 1. We find that the nonlinear time evolution of the film thickness is captured well by the linear approximation for small times. The analytically computed profile from the linearized model and the numerical solution of the full nonlinear equation are plotted together for comparison in Figure 9, together with the $L^{2}$ -distance between them as a function of time. While for small times there is good agreement between the linear and nonlinear dynamics, as time goes by, the two solutions deviate significantly.

If we consider both $\theta$ and $\phi$ dependence and carry out the simulations using the two-dimensional equation, with $l=3$ and different choices of $m$ , we can verify the analytical stability result with varying degrees of accuracy depending on the value of $m$ . Since the three-dimensional results are difficult to compare visually, we also provide a comparison of the height profiles along the a spiral line on the sphere: $\phi=\pi(1-\cos\theta)$ , which corresponds to the diagonal line of the rectangular domain $-1<x<1$ and $0<\phi<2\pi$ .

The result for $m=0$ through $m=3$ are shown in Figures 10–13 respectively. For the axisymmetric $m=0$ case, good agreement remains through time $t=0.01$ while for the $m=2$ and $m=3$ cases, similar agreement can only be obtained until a much shorter time of $t=0.001$ . The case $m=1$ yields the worst agreement, which is associated with the Legendre function $P_{3}^{1}(x)$ being the least smooth among $P_{3}^{m}(x)$ . The $L^{2}$ -difference between computational and analytical results for different $m$ are shown in Figure 14, illustrating that the linear model approximates the short time non-linear dynamics sufficiently well.

Using a similar procedure to the earlier axisymmetric case, the growth rate derived from the simulation is shown in Table 2. For each case, the growth rates are estimated using both $L^{\infty}$ and $L^{1}$ norms of the numerical results. Whereas the analytical growth rate is predicted to be 120 for all cases, the numerical results show varying degrees of agreement, with differences as high as 18% in some cases.

IV.3 Energy functional for the axisymmetric case

When the film profile is independent of $\phi$ , we have been able to find a corresponding energy functional whose time-derivative is negative under the governing dynamics, showing that the system is dissipative. We still take the temperature field to to be purely in the radial direction and gravity to be zero. The corresponding evolution equation is given by

[TABLE]

where parameter $\cal S$ relates to surface tension and $\cal N$ to the Marangoni effect.

Multiply $h$ by the multiplier

[TABLE]

and integrate over $x$ in $[-1,1]$ to obtain the functional $E(h)$ :

[TABLE]

Here, subscripts on $h$ refer to corresponding partial derivatives.

Take the derivative of $E$ with respect to time $t$ and simplify using integration-by-parts and by invoking the evolution equation to relate the time-derivative $h_{t}$ to the spatial derivative term, resulting in

[TABLE]

As $h$ and $(1-x^{2})$ are non-negative on $[-1,1]$ , we conclude that $\frac{d}{dt}E(h)\leq 0$ .

Mathematically, the energy functional (21) allows one to prove the existence of the minimizer for the axisymmetric case. However, in addition to the energy decay property, one must also show that this functional is bounded from below. To establish a lower bound for $E(h)$ , we can use the result from kang2016weak that $h(x,t)$ belongs to a Sobolev space $H^{1}$ with weight $1-x^{2}$ . This establishes that the integral $\int_{-1}^{1}(1-x^{2})h_{x}^{2}dx$ is bounded from below. The other two terms contributing to $E(h)$ are easily seen to be bounded from below and therefore so is $E(h)$ itself.

Upon minimizing $E(h)$ subject to the mass conservation constraint

[TABLE]

we obtain the Euler-Lagrange equation with a Lagrange multiplier $\lambda$ (note that $\lambda$ in this section is not the same as the growth rate from the linearized stability analysis):

[TABLE]

Numerically we find that with a slightly perturbed initial condition $h_{0}(x)=1+0.001\sin(2\pi x)$ and parameters ${\cal G}=0$ , ${\cal S}=1$ and ${\cal N}=-2$ , the thin film will become uniform and during this process, energy will approach to an equilibrium corresponding to $\lambda=2$ (see Figure 15).

It turns out to be possible to include the effects of gravity, characterized by parameter $\cal G$ , in this derivation. In that case, one can show that the energy functional becomes:

[TABLE]

with the corresponding Euler-Lagrange equation:

[TABLE]

These nonlinear differential equations that include the logarithmic term $\ln(h)$ can be solved to obtain the long-time, energy-minimizing, equilibrium shape of the axisymmetric film in the presence of a radial temperature gradient. Numerically we find that with a slightly perturbed initial condition $h_{0}(x)=1+0.001\sin(2\pi x)$ and parameters ${\cal G}=2$ , ${\cal S}=1$ and ${\cal N}=-3$ , gravity will drive the fluid to accumulate near the bottom of the sphere and during this process, the energy approaches an equilibrium corresponding to the value $\lambda=2.21$ (see Figure 16).

V Discussion

For the case of a thin viscous liquid film coating a sphere with an imposed axial (i.e. vertical) or radial temperature gradient, we have been able to derive a closed-form nonlinear fourth-order partial differential equation that governs the evolution of the film thickness accounting for gravitational, capillary and Marangoni effects.

Assuming a linear variation of surface tension with temperature and a linearly varying axial temperature field, we found that there is a particular uniform film thickness (i.e., fluid volume) where, despite a continual internal flow, the film profile remains constant and steady. In that state, the downward draining flow due to gravity is balanced exactly by the upward flow due to Marangoni effects (provided that the surface tension near the north pole is higher than that near the south pole). Moreover, our analysis showed that if the fluid volume is slightly increased from that state, the fluid tends to accumulate near the bottom of the sphere, whereas if the volume is slightly decreased, the film thins out near the bottom and bulges out above, though not exactly at the north pole. This analysis suggests a potential experimental method for measuring, indirectly, the slope of the surface tension versus temperature curve. For instance, one can slowly add volume to a liquid coating to see at what point the film thickness become exactly uniform. Conversely, for a fluid that might slowly evaporate over time, one can start with excess fluid, which produces a bulge near the bottom of the sphere and observe, as the volume decreases slowly by evaporation, the film first become uniform in thickness and then thin out near the bottom. If the (dimensional) steady uniform film thickness $h_{s}$ can be measured in this manner, the parameter $\sigma_{1}$ can be obtained from our theoretical prediction: $\sigma_{1}=2\rho gh_{s}/3k_{1}$ . Typically, both $k_{1}$ and $\sigma_{1}$ would be negative numbers in this case.

For a radially varying temperature field, we were able to show that an initially uniform film may become unstable provided that surface tension becomes lower in regions where the film thins out and higher where it bulges out. So, if surface tension decreases with increasing temperature, the imposed temperature field must be higher at the sphere surface and decrease radially outward. This is consistent with the finding by wilson1994onset . In the axisymmetric case, we could also find an energy functional with the dissipation property and bounded from below, so that the evolution of the film thickness takes the energy level toward a minimum. To show that the Euler-Lagrange equation that describes the minimum of the energy functional does have a proper mathematical structure, we can apply local analysis to show that it admits zero-contact-angle non-negative solutions near the transition between dry and wet zones on the surface. Without loss of generality, let us assume that the touch-down point is at $x_{0}=0$ . Hence, locally in the neighborhood of $x_{0}$ the differential equation (22) simplifies to

[TABLE]

where $\alpha={\cal N}/{\cal S}$ and $\beta=\lambda/{\cal S}$ . By reduction of order, i.e., multiplying the equation by $h^{\prime}$ and integrating we obtain

[TABLE]

where $c$ is the integration constant. A zero contact angle $h^{\prime}(0)=0$ implies that the integration constant is $c=0$ resulting in

[TABLE]

To get a local non-negative solution we would choose the plus sign. Local existence follows from the continuity of the right side of the equation at $h(0)=0$ at $0^{+}$ .

Acknowledgements.

This work was partially supported by a grant from the Simons Foundation (#275088 to MC). DK thanks the Institute of Mathematical Sciences at Claremont Graduate University for a Director’s Fellowship that supported his doctoral studies.

Bibliography23

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1(1) Kabov, O. A., Scheid, B., Sharina, I. A., & Legros, J. C. (2002). Heat transfer and rivulet structures formation in a falling thin liquid film locally heated. International Journal of Thermal Sciences, 41(7), 664-672.
2(2) Kalliadasis, S., Kiyashko, A., & Demekhin, E. A. (2003). Marangoni instability of a thin liquid film heated from below by a local heat source. Journal of Fluid Mechanics, 475, 377-408.
3(3) Bankoff, S. G. (1971). Stability of liquid flow down a heated inclined plane. International Journal of Heat and Mass Transfer, 14(3), 377-385.
4(4) Zuber, N., & Staub, F. W. (1966). Stability of dry patches forming in liquid films flowing over heated surfaces. International Journal of Heat and Mass Transfer, 9(9), 897-905.
5(5) Oron, A. (2000). Nonlinear dynamics of three-dimensional long-wave Marangoni instability in thin liquid films. Physics of Fluids, 12(7), 1633-1645.
6(6) Bestehorn, M., Pototsky, A., & Thiele, U. (2003). 3D large scale Marangoni convection in liquid films. The European Physical Journal B-Condensed Matter and Complex Systems, 33(4), 457-467.
7(7) Reisfeld, B., & Bankoff, S. G. (1990). Nonlinear stability of a heated thin liquid film with variable viscosity. Physics of Fluids A: Fluid Dynamics, 2(11), 2066-2067.
8(8) Jensen, O. E., & Grotberg, J. B. (1993). The spreading of heat or soluble surfactant along a thin liquid film. Physics of Fluids A: Fluid Dynamics, 5(1), 58-68.