2.1. Problem description and key assumptions

The films investigated in Lalli et al. (Reference Lalli, Shen, Dini and Giusti2023) were circular, almost horizontal, formed on a glass boundary, and were created from aqueous-surfactant solutions containing glycerol and dilute concentrations of magnetite NPs. A cylindrical neodymium (Nd$_{2}$Fe$_{14}$B) magnet with diameter 38.5 mm, thickness 10 mm and remanence 1.43 T was used to apply an inhomogeneous magnetic field to each film, as shown in figure 1(a). We refer the reader to Lalli et al. (Reference Lalli, Shen, Dini and Giusti2023) for additional details about the experiments.

Figure 1. A schematic of a magnetic soap film supported by a glass boundary. Panel (a) illustrates the thin film region, which is the region that is modelled and simulated in the present study, in the $x$–$y$ plane. The thin film region is denoted by $\varOmega$; this region has boundary $\partial \varOmega$ and outward-facing unit normal $\boldsymbol {n}_{\rm d}$. A right-handed rectangular Cartesian coordinate system with basis $\{\boldsymbol {e}_{x},\boldsymbol {e}_{y}, \boldsymbol {e}_{z}\}$, coordinates $\{x, y, z\}$ and origin in the bottom left of the thin film is used throughout the derivation. Panel (b) presents a zoomed-in cross-section of a small section of the thin film in (a) in the $x$–$z$ plane. At the dividing surface, the outward-facing unit normal is denoted by $\boldsymbol {n}$ and the tangential vector in this plane by $\boldsymbol {t}_x$. The thin film was assumed to be symmetrical about a centre surface, denoted by $z=C(x,y)$.

The phase interface that separates the core liquid from the gas phase is a three-dimensional region with a thickness of the order of a few angstroms. We modelled the phase interface by defining a two-dimensional dividing surface within the interface region (Slattery Reference Slattery1967; Sagis Reference Sagis2011). The effect of the phase interface on the surrounding phases is accounted for by defining excess quantities on the dividing surface, such as the surface tension, $\gamma$.

A rectangular Cartesian coordinate system, with basis $\{\boldsymbol {e}_{x},\boldsymbol {e}_{y}, \boldsymbol {e}_{z}\}$ and coordinates $\{x, y, z\}$, was used to derive the system of equations governing the thickness evolution of magnetic soap films. The film thickness equation was derived for the half-film thickness, $h$, by assuming symmetry about a centre surface, $C$, with uniform curvature. Since dilute concentrations of magnetite NPs in aqueous-surfactant solutions were used in the experiments to create magnetic soap films, the core liquid was modelled as a water-based ferrofluid containing a dilute concentration of magnetite NPs. Furthermore, the flow of the core liquid was modelled as an incompressible flow of a Newtonian fluid with uniform shear viscosity approximated by that of water, and magnetic dipole–dipole interactions, magnetoviscous effects and demagnetising field effects associated with the magnetisation of the film were neglected. Since the magnetic soap films in the experiments were placed in the symmetry plane of the neodymium magnet such that the applied magnetic field was approximately in the $x$–$y$ plane throughout the films, the gradient $\partial {H}/\partial {z}$ was neglected throughout the films, where $H$ is the magnitude of the magnetic field intensity. The deformation response of the interface was assumed to be governed by the Boussinesq–Scriven surface stress model with uniform surface shear and dilatational viscosities. It was also assumed, for simplicity, that the surfactants stabilising each magnetic soap film were insoluble and that the magnetite NPs in each film were monodisperse.

Figure 1 presents a schematic of the modelling configuration with the coordinate system and several of the key assumptions illustrated. It should be highlighted that it is the thin film, denoted by $\varOmega$ and with boundary $\partial \varOmega$, that is being modelled in the present study. The following notation is used throughout: a $'$ is appended to fluid properties to denote that they describe the gas phase surrounding the film. When no $'$ is present, the fluid properties describe the core liquid (e.g. $\rho '$ is the density of the surrounding gas and $\rho$ is the density of the core liquid). Furthermore, a subscript ‘$\mathrm {s}$’ is used to denote fields defined at the dividing surface. To avoid the derivation becoming too verbose, the centre surface, $C$, is initially discarded and is reintroduced in the final equation for the thickness evolution.

2.2. Governing equations in the core liquid and at the dividing surface

The conservation of mass and balance of momentum demand

(2.1a)$$\begin{gather} \frac{\partial \rho}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\rho \boldsymbol{V}) = 0 , \end{gather}$$

(2.1b)$$\begin{gather}\rho\frac{\partial \boldsymbol{V}}{\partial t} + \rho(\boldsymbol{V} \boldsymbol{\cdot} \boldsymbol{\nabla})\boldsymbol{V} = \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\mathsf{T}} - \boldsymbol{\nabla} \phi + \rho \boldsymbol{g} , \end{gather}$$

where $\rho$ is density, $t$ is time, $\boldsymbol {V}$ is the velocity, $\boldsymbol{\mathsf{T}}$ is the stress tensor, $\phi$ is the conjoining pressure, through which interactions between the surfactant-laden interfaces can be accounted for (Craster & Matar Reference Craster and Matar2009), and $\boldsymbol {g}$ is the gravitational acceleration, which is given by $\boldsymbol {g} = - g \boldsymbol {e}_{z}$. The stress tensor can be expressed as

(2.2)\begin{equation} \boldsymbol{\mathsf{T}} ={-}p\boldsymbol{\mathsf{I}} + \boldsymbol{\mathsf{T}}_{\mathrm{v}} + \boldsymbol{\mathsf{T}}_{\rm m} , \end{equation}

where $p$ is the thermodynamic pressure, $\boldsymbol{\mathsf{I}}$ is the identity tensor, $\boldsymbol{\mathsf{T}}_{\mathrm {v}}$ is the viscous stress tensor and $\boldsymbol{\mathsf{T}}_{\rm m}$ is the magnetic part of the stress tensor. For a Newtonian fluid,

(2.3)\begin{equation} \boldsymbol{\mathsf{T}}_{\mathrm{v}} = \eta (\boldsymbol{\nabla} \boldsymbol{V}+ (\boldsymbol{\nabla}\boldsymbol{V})^{\intercal}) + \lambda(\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{V})\boldsymbol{\mathsf{I}} , \end{equation}

where $\eta$ is the shear viscosity, $\lambda$ is the dilatational viscosity and $^{\intercal }$ is the transpose operator. For a ferrofluid containing a dilute concentration of magnetite NPs and assuming magnetostriction effects are negligible (Rosensweig Reference Rosensweig2013), the magnetic part of the stress tensor given in Cowley & Rosensweig (Reference Cowley and Rosensweig1967) simplifies to

(2.4)\begin{equation} \boldsymbol{\mathsf{T}}_{\rm m} ={-}\tfrac{1}{2}\mu_{0}H^{2}\boldsymbol{\mathsf{I}} + \boldsymbol{B}\boldsymbol{H} , \end{equation}

where $\mu _0$ is the permeability of free space, $\boldsymbol {B}$ is the magnetic flux density, $\boldsymbol {H}$ is the magnetic field intensity with magnitude $H$ and $\boldsymbol {B} \boldsymbol {H}$ is the tensor product of $\boldsymbol {B}$ and $\boldsymbol {H}$. In the absence of demagnetising field effects due to the magnetisation of the film, $\boldsymbol {H}$ is the field due to the neodymium magnet. The magnetic dipoles of the magnetite NPs rotate to align with the applied field by either Brownian relaxation (rotation of the particle) or Néel relaxation (rotation of the magnetic moment within the particle) (Rosensweig Reference Rosensweig1987). The time scale for magnetic relaxation (Rosensweig Reference Rosensweig2013) of the magnetisation, $\boldsymbol {M}$, of the experimentally investigated magnetic soap films is expected to be orders of magnitude less than the characteristic flow time scale (Lalli et al. Reference Lalli, Shen, Dini and Giusti2023). As a result, the magnetisation was approximately parallel with the magnetic field intensity during the lifetime of the investigated magnetic soap films. For an electrically non-conducting fluid containing negligible displacement currents ($\boldsymbol {\nabla } \times \boldsymbol {H} = \boldsymbol {0}$) and with $\boldsymbol {M}$ and $\boldsymbol {H}$ parallel, the magnetic force per unit volume, $\boldsymbol {f}_{\rm m} = \boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol{\mathsf{T}}_{\rm m}$, can be expressed as (Neuringer & Rosensweig Reference Neuringer and Rosensweig1964; Cowley & Rosensweig Reference Cowley and Rosensweig1967)

(2.5)\begin{equation} \boldsymbol{f}_{\rm m} = \mu_{0} M \boldsymbol{\nabla} H , \end{equation}

where $M$ is the magnitude of $\boldsymbol {M}$. When considering the flow of core liquid as an incompressible flow of a Newtonian fluid with uniform shear viscosity, the mass conservation (2.1a) and momentum balance (2.1b) with the stress tensor (2.2) inserted are

(2.6a)$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{V} = 0 , \end{gather}$$

(2.6b)$$\begin{gather}\rho\frac{\partial \boldsymbol{V}}{\partial t} + \rho(\boldsymbol{V} \boldsymbol{\cdot} \boldsymbol{\nabla})\boldsymbol{V} ={-} \boldsymbol{\nabla} (p + \phi) + \eta \nabla^2 \boldsymbol{V} + \mu_0 M \boldsymbol{\nabla} H + \rho \boldsymbol{g} , \end{gather}$$

where $\nabla ^2$ is the Laplace operator.

The dividing surface was defined at the position of 0 surface mass density (Slattery, Sagis & Oh Reference Slattery, Sagis and Oh2007). The jump mass balance, which is a statement of conservation of mass at each point on the dividing surface, then requires (Burelbach, Bankoff & Davis Reference Burelbach, Bankoff and Davis1988; Slattery et al. Reference Slattery, Sagis and Oh2007)

(2.7)\begin{equation} \rho(\boldsymbol{V} - \boldsymbol{V}_{\rm s}) \boldsymbol{\cdot} \boldsymbol{n} = \rho'(\boldsymbol{V}' - \boldsymbol{V}_{\rm s}) \boldsymbol{\cdot} \boldsymbol{n} \quad \text{at } \varSigma , \end{equation}

where $\boldsymbol {V}_{\rm s}$ is the surface velocity (rate of change of position following a particle on the surface), $\boldsymbol {n}$ is the normal vector illustrated in figure 1 and $\varSigma$ denotes the dividing surface. In continuum fluid mechanics, it is assumed that the velocity tangential to the phase interface is continuous across the phase interface (Slattery et al. Reference Slattery, Sagis and Oh2007),

(2.8)\begin{equation} \boldsymbol{\mathsf{P}}_{\rm s} \boldsymbol{\cdot} \boldsymbol{V} = \boldsymbol{\mathsf{P}}_{\rm s} \boldsymbol{\cdot} \boldsymbol{V}_{\rm s} \quad \text{at } \varSigma , \end{equation}

where $\boldsymbol{\mathsf{P}}_{\rm s}$ is the surface projection tensor, defined by $\boldsymbol{\mathsf{P}}_{\rm s} = \boldsymbol{\mathsf{I}} - \boldsymbol {n} \boldsymbol {n}$. The mass flux due to evaporation, $J_{\rm e}$, is given by (Burelbach et al. Reference Burelbach, Bankoff and Davis1988; Bai & Gosman Reference Bai and Gosman1996)

(2.9)\begin{equation} J_{\rm e} = \rho' (\boldsymbol{V}' - \boldsymbol{V}_{\rm s}) \boldsymbol{\cdot} \boldsymbol{n} \quad \text{at } \varSigma . \end{equation}

Therefore, $\boldsymbol {V}$ is only equal to $\boldsymbol {V}_{\rm s}$ at the dividing surface when there is no evaporation. The jump momentum balance, which is a statement of momentum balance at each point on the dividing surface, requires (Slattery et al. Reference Slattery, Sagis and Oh2007)

(2.10)\begin{equation} \boldsymbol{\nabla}_{\rm s} \boldsymbol{\cdot} \boldsymbol{\mathsf{T}}_{\rm s} + (\boldsymbol{\mathsf{T}}' - \boldsymbol{\mathsf{T}}) \boldsymbol{\cdot} \boldsymbol{n} - \frac{J_{\rm e}^2}{\rho'}\left(1 - \frac{\rho'}{\rho} \right) \boldsymbol{n} = \boldsymbol{0} \quad \text{at } \varSigma, \end{equation}

where $\boldsymbol {\nabla }_{\rm s}$ and $\boldsymbol {\nabla }_{\rm s} \boldsymbol {\cdot }$ are the surface gradient and surface divergence operators, respectively, and $\boldsymbol{\mathsf{T}}_{\rm s}$ is the surface stress tensor. Appendix A provides additional details about these surface operators. The surface stress tensor $\boldsymbol{\mathsf{T}}_{\rm s}$ can be expressed as the sum of contributions due to the surface tension, $\gamma$, and a surface stress tensor $\boldsymbol{\mathsf{S}}_{\rm s}$ that accounts for surface viscous effects due to the presence of surfactant molecules at the phase interface,

(2.11)\begin{equation} \boldsymbol{\mathsf{T}}_{\rm s} = \gamma \boldsymbol{\mathsf{P}}_{\rm s} + \boldsymbol{\mathsf{S}}_{\rm s} . \end{equation}

Using the surface stress tensor (2.11) and the stress tensor (2.2) in the jump momentum balance (2.10) results in

(2.12)\begin{align} &\boldsymbol{\nabla}_{\rm s} \gamma + 2 \kappa \gamma \boldsymbol{n} + \boldsymbol{\nabla}_{\rm s} \boldsymbol{\cdot} \boldsymbol{\mathsf{S}}_{\rm s} + (p - p')\boldsymbol{n} + (\boldsymbol{\mathsf{T}}_{\mathrm{v}}' - \boldsymbol{\mathsf{T}}_{\mathrm{v}}) \boldsymbol{\cdot} \boldsymbol{n} + (\boldsymbol{\mathsf{T}}_{\rm m}' - \boldsymbol{\mathsf{T}}_{\rm m}) \boldsymbol{\cdot} \boldsymbol{n}\nonumber\\ &\quad = \frac{J_{\rm e}^2}{\rho'}\left(1 - \frac{\rho'}{\rho} \right) \boldsymbol{n} \quad \text{at } \varSigma , \end{align}

where $\kappa$ is the mean curvature at each point on the dividing surface, which is given by

(2.13)\begin{equation} \kappa ={-} \tfrac{1}{2} \boldsymbol{\nabla}_{\rm s} \boldsymbol{\cdot} \boldsymbol{n} . \end{equation}

For a non-magnetisable gas surrounding the magnetic soap film and from the magnetostatic equations (Rosensweig Reference Rosensweig2013),

(2.14)\begin{equation} (\boldsymbol{\mathsf{T}}_{\rm m}' - \boldsymbol{\mathsf{T}}^{}_{\rm m}) \boldsymbol{\cdot} \boldsymbol{n} = \tfrac{1}{2}\mu_{0} (\boldsymbol{M} \boldsymbol{\cdot} \boldsymbol{n})^{2} \boldsymbol{n} . \end{equation}

The Boussinesq–Scriven surface stress model is a surface analogue of the stress deformation relation for a Newtonian fluid. With the Boussinesq–Scriven surface stress model (Slattery et al. Reference Slattery, Sagis and Oh2007),

(2.15)\begin{equation} \boldsymbol{\mathsf{S}}_{\rm s} = 2 \eta_{\rm s} \boldsymbol{\mathsf{D}}_{\rm s} + (\lambda_{\rm s} - \eta_{\rm s}) (\boldsymbol{\nabla}_{\rm s} \boldsymbol{\cdot} \boldsymbol{V}_{\rm s}) \boldsymbol{\mathsf{P}}_{\rm s} , \end{equation}

where $\eta _{\rm s}$ and $\lambda _{\rm s}$ are the surface shear and surface dilatational viscosities, respectively, and $\boldsymbol{\mathsf{D}}_{\rm s}$ is the surface rate of deformation tensor, which is given by

(2.16)\begin{equation} \boldsymbol{\mathsf{D}}_{\rm s} = \tfrac{1}{2} ( \boldsymbol{\mathsf{P}}_{\rm s} \boldsymbol{\cdot} (\boldsymbol{\nabla}_{\rm s} \boldsymbol{V}_{\rm s}) + (\boldsymbol{\nabla}_{\rm s} \boldsymbol{V}_{\rm s})^{\intercal} \boldsymbol{\cdot} \boldsymbol{\mathsf{P}}_{\rm s}) . \end{equation}

Taking the inner product of (2.12) with three linearly independent vectors leads to three scalar jump momentum balance equations, which serve as boundary conditions for the governing equations in the core liquid (2.6). It is convenient to use $\{\boldsymbol {t}_x, \boldsymbol {t}_y, \boldsymbol {n}\}$, where $\boldsymbol {t}_x$ and $\boldsymbol {t}_y$ are two linearly independent vectors tangential to the film surface. In terms of the half-film thickness, $\{\boldsymbol {t}_x, \boldsymbol {t}_y, \boldsymbol {n}\}$ are

(2.17a)$$\begin{gather} \boldsymbol{t}_{x} = \left(1 + \left(\frac{\partial h}{\partial x}\right)^2 \right)^{-{1}/{2}} \left( \boldsymbol{e}_{x} + \frac{\partial{h}}{\partial{x}} \boldsymbol{e}_{z} \right) , \end{gather}$$

(2.17b)$$\begin{gather}\boldsymbol{t}_{y} = \left(1 + \left(\frac{\partial h}{\partial y}\right)^2 \right)^{-{1}/{2}} \left( \boldsymbol{e}_{y} + \frac{\partial{h}}{\partial{y}} \boldsymbol{e}_{z} \right) , \end{gather}$$

(2.17c)$$\begin{gather}\boldsymbol{n} = \left(1 + \left(\frac{\partial h}{\partial x}\right)^{2} + \left(\frac{\partial h}{\partial y}\right)^{2} \right)^{-{1}/{2}} \left( -\frac{\partial{h}}{\partial{x}} \boldsymbol{e}_{x} - \frac{\partial{h}}{\partial{y}} \boldsymbol{e}_{y} + \boldsymbol{e}_{z} \right) . \end{gather}$$

Another condition that must be satisfied at the dividing surface is the kinematic boundary condition, which states that for a function $\varPhi$ that is 0 at each point on the dividing surface (Slattery et al. Reference Slattery, Sagis and Oh2007)

(2.18)\begin{equation} \frac{\partial \varPhi}{\partial t} + \boldsymbol{V}_{{\rm s}} \boldsymbol{\cdot} \boldsymbol{\nabla} \varPhi = 0 . \end{equation}

In this instance, $\varPhi (x,y,z,t) = z - h(x,y,t)$, and the kinematic boundary condition (2.18) is

(2.19)\begin{equation} \frac{\partial h}{\partial t} + u_{\rm s}\frac{\partial h}{\partial x} + v_{\rm s}\frac{\partial h}{\partial y} = w_{\rm s} , \end{equation}

where $\{u_{\rm s}, v_{\rm s}, w_{\rm s}\}$ are the rectangular Cartesian components of $\boldsymbol {V}_{\rm s}$. Furthermore, symmetry boundary conditions are used at the mid-plane of the film,

(2.20)\begin{equation} \frac{\partial u}{\partial z} = 0 , \quad \frac{\partial v}{\partial z} = 0 ,\quad \text{and } w = 0\quad \text{at } z = 0 , \end{equation}

where $\{u, v, w\}$ are the rectangular Cartesian components of $\boldsymbol {V}$.

2.3. Thin film approximation

With a characteristic thickness scale $\mathcal {H}$ and a characteristic length scale $L$, the dimensionless coordinates and thickness are

(2.21a–d)\begin{equation} \tilde{x} = \frac{x}{L} , \quad \tilde{y} = \frac{y}{L}, \quad \tilde{z} = \frac{z}{\mathcal{H}}, \quad \tilde{h} = \frac{h}{\mathcal{H}} , \end{equation}

where $\tilde {\phantom {-}}$ denotes a dimensionless variable. The thin film approximation (or lubrication approximation) involves capitalising on the small size of the thin film parameter, $\epsilon$, which is defined by

(2.22)\begin{equation} \epsilon = \frac{\mathcal{H}}{L} \ll 1 . \end{equation}

All equations are written in dimensionless form and the dominant physics is extracted by expanding dimensionless dependent variables using a series expansion in powers of $\epsilon$, for example,

(2.23a)$$\begin{gather} \tilde{u}(\tilde{x},\tilde{y},\tilde{z},\tilde{t};\epsilon) = \tilde{u}_{0}(\tilde{x},\tilde{y},\tilde{z},\tilde{t}) + \epsilon \tilde{u}_{1}(\tilde{x},\tilde{y},\tilde{z},\tilde{t}) + \epsilon^2 \tilde{u}_{2}(\tilde{x},\tilde{y},\tilde{z},\tilde{t}) + \cdots , \end{gather}$$

(2.23b)$$\begin{gather}\tilde{h}(\tilde{x},\tilde{y},\tilde{t};\epsilon) = \tilde{h}_{0}(\tilde{x},\tilde{y},\tilde{t}) + \epsilon \tilde{h}_{1}(\tilde{x},\tilde{y},\tilde{t}) + \epsilon^2 \tilde{h}_{2}(\tilde{x},\tilde{y},\tilde{t}) + \cdots, \end{gather}$$

and retaining terms at leading order in $\epsilon$. The velocity scale $U$, pressure scale $P$ and scales $M_{{\rm c}}$ and $H_{{\rm c}}$ for the magnetisation and magnetic field intensity are used to introduce the following dimensionless variables:

(2.24)\begin{equation} \left.\begin{gathered} \tilde{u} = \frac{u}{U} , \quad \tilde{v} = \frac{v}{U} , \quad \tilde{w} = \frac{w}{\epsilon U} , \quad \tilde{u}_{\rm s} = \frac{u_{\rm s}}{U} , \quad \tilde{v}_{\mathrm{s}} = \frac{v_{\rm s}}{U} , \quad \tilde{w}_{\rm s} = \frac{w_{\rm s}}{\epsilon U},\\ \tilde{t} = \frac{t U}{L} ,\quad \tilde{p} = \frac{p}{P} , \quad \tilde{\phi} = \frac{\phi}{P} ,\quad \tilde{M}_x= \frac{M_x}{M_{\mathrm{c}}} , \quad \tilde{M}_y = \frac{M_y}{M_{\mathrm{c}}} ,\\ \tilde{M}_z = \frac{M_z}{\epsilon M_{\mathrm{c}}} , \quad \tilde{M} = \frac{M}{M_\mathrm{c}} , \quad \tilde{H} = \frac{H}{H_\mathrm{c}} . \end{gathered}\right\} \end{equation}

Here $\{M_x, M_y, M_z\}$ are the rectangular Cartesian components of the magnetisation $\boldsymbol {M}$. Note that $M_z$ has been scaled by $\epsilon M_{\mathrm{c}}$ as it is expected that $M_z \ll M_x$ and $M_z \ll M_y$ due to the direction of the applied magnetic field. From (2.8), we find at leading order in $\epsilon$ that

(2.25)\begin{equation} u = u_{\rm s} \quad \text{and} \quad v = v_{\rm s} \quad \text{at } \varSigma . \end{equation}

We use (2.25) with the kinematic boundary condition (2.19) and the result of integrating the continuity equation (2.6a) over the film thickness using the Leibniz integral rule and the symmetry boundary condition (2.20) to find that

(2.26)\begin{equation} \frac{\partial{h}}{\partial{t}} + \frac{\partial{Q_x}}{\partial{x}} + \frac{\partial{Q_y}}{\partial{y}} = w_{\rm s} - w(x,y,z=h,t) , \end{equation}

where $Q_x$ and $Q_y$ are the following depth-integrated velocities:

(2.27a,b)\begin{equation} Q_{x} = \int_{0}^{h} u \, \mathrm{d} z , \quad Q_{y} = \int_{0}^{h} v \, \mathrm{d} z . \end{equation}

At leading order in $\epsilon$ and at the dividing surface, $(\boldsymbol {V} - \boldsymbol {V}_{\rm s}) \boldsymbol {\cdot } \boldsymbol {n} = w(x,y,z=h,t) - w_{\rm s}$, which allows the mass flux due to evaporation to be introduced into (2.26),

(2.28)\begin{equation} \frac{\partial{h}}{\partial{t}} + \frac{\partial{Q_x}}{\partial{x}} + \frac{\partial{Q_y}}{\partial{y}} ={-} \frac{J_{\rm e}}{\rho} . \end{equation}

In dimensionless form, (2.28) reads

(2.29)\begin{equation} \frac{\partial \tilde{h}}{\partial \tilde{t}} + \frac{\partial \tilde{Q}_x}{\partial \tilde{x}} + \frac{\partial \tilde{Q}_y}{\partial \tilde{y}} ={-} \tilde{J}_{\rm e} , \end{equation}

where

(2.30a–c)\begin{equation} \tilde{Q}_x = \int_{0}^{\tilde{h}} \tilde{u} \, \mathrm{d} \tilde{z} = \frac{Q_x}{U\mathcal{H}} , \quad \tilde{Q}_y = \int_{0}^{\tilde{h}} \tilde{v} \, \mathrm{d} \tilde{z} = \frac{Q_y}{U\mathcal{H}} , \quad \tilde{J}_{\rm e} = \frac{J_{\rm e}}{\rho U \epsilon} . \end{equation}

The focus now turns to finding $\tilde {u}$ and $\tilde {v}$ so that $\tilde {Q}_x$ and $\tilde {Q}_y$ can be evaluated. The pressure scale $P = {\eta U}/{L}$ is relevant to extensional (plug) flows (Kargupta, Sharma & Khanna Reference Kargupta, Sharma and Khanna2004; Münch, Wagner & Witelski Reference Münch, Wagner and Witelski2005). Since we expect viscous shear from the core liquid to enter the problem at leading order, we instead use the pressure scale $P = {\eta U}/(\epsilon \mathcal {H})$, which allows the derived model to be used in the limit of immobile interfaces. With the pressure scale $P = {\eta U}/(\epsilon \mathcal {H})$, the governing equations for the core liquid (2.6) in dimensionless form are

(2.31a)\begin{gather} \frac{\partial \tilde{u}}{\partial \tilde{x}} + \frac{\partial \tilde{v}}{\partial \tilde{y}} + \frac{\partial \tilde{w}}{\partial \tilde{z}} = 0, \end{gather}

(2.31b)\begin{gather} \epsilon Re \left( \frac{\partial \tilde{u}}{\partial \tilde{t}} + \tilde{u}\frac{\partial \tilde{u}}{\partial \tilde{x}} + \tilde{v}\frac{\partial \tilde{u}}{\partial \tilde{y}} + \tilde{w}\frac{\partial \tilde{u}}{\partial \tilde{z}}\right) =-& \frac{\partial (\kern0.7pt \tilde{p} + \tilde{\phi})}{\partial \tilde{x}} + \epsilon^2 \left(\frac{\partial^{2} \tilde{u}}{\partial \tilde{x}^{2}} + \frac{\partial^{2} \tilde{u}}{\partial \tilde{y}^{2}}\right)\nonumber\\ &+ \,\frac{\partial^{2} \tilde{u}}{\partial \tilde{z}^{2}} + \varPsi \tilde{M}\frac{\partial \tilde{H}}{\partial \tilde{x}} , \end{gather}

(2.31c)\begin{gather} \epsilon Re \left( \frac{\partial \tilde{v}}{\partial \tilde{t}} + \tilde{u}\frac{\partial \tilde{v}}{\partial \tilde{x}} + \tilde{v}\frac{\partial \tilde{v}}{\partial \tilde{y}} + \tilde{w}\frac{\partial \tilde{v}}{\partial \tilde{z}} \right) =-& \frac{\partial (\kern0.7pt \tilde{p} + \tilde{\phi})}{\partial \tilde{y}} + \epsilon^2 \left( \frac{\partial^{2} \tilde{v}}{\partial \tilde{x}^{2}} + \frac{\partial^{2} \tilde{v}}{\partial \tilde{y}^{2}} \right)\nonumber\\ &+\, \frac{\partial^{2} \tilde{v}}{\partial \tilde{z}^{2}} + \varPsi \tilde{M}\frac{\partial \tilde{H}}{\partial \tilde{y}}, \end{gather}

(2.31d)\begin{gather} \epsilon^3 Re \left( \frac{\partial \tilde{w}}{\partial \tilde{t}} + \tilde{u} \frac{\partial \tilde{w}}{\partial \tilde{x}} + \tilde{v} \frac{\partial \tilde{w}}{\partial \tilde{y}} + \tilde{w} \frac{\partial \tilde{w}}{\partial \tilde{z}} \right) =&-\frac{\partial (\kern0.7pt \tilde{p} + \tilde{\phi})}{\partial \tilde{z}} + \epsilon^4 \left( \frac{\partial^{2} \tilde{w}}{\partial \tilde{x}^{2}} + \frac{\partial^{2} \tilde{w}}{\partial \tilde{y}^{2}} \right)\nonumber\\ &+ \,\epsilon^2 \frac{\partial^{2} \tilde{w}}{\partial \tilde{z}^{2}} - \epsilon \mathcal{G} , \end{gather}

where $Re$ is the Reynolds number, $\varPsi$ is a dimensionless magnetic number and $\mathcal {G}$ is a dimensionless gravitational number. These dimensionless numbers are defined by

(2.32a–c)\begin{equation} Re = \frac{\rho U \mathcal{H}}{\eta} , \quad \varPsi = \frac{\mathcal{H} \epsilon \mu_0 M_{\mathrm{c}} H_{\mathrm{c}}}{\eta U} , \quad \mathcal{G} = \frac{\rho g \mathcal{H}^2}{\eta U} . \end{equation}

The Reynolds number can be expressed as $Re = \epsilon Re^{\star }$, where $Re^{\star }$ is of the order of 1 or less. As a result, (2.31) at leading order in $\epsilon$ is

(2.33a)$$\begin{gather} \frac{\partial \tilde{u}_0}{\partial \tilde{x}} + \frac{\partial \tilde{v}_{0}}{\partial \tilde{y}} + \frac{\partial \tilde{w}_{0}}{\partial \tilde{z}} = 0 , \end{gather}$$

(2.33b)$$\begin{gather}0 ={-} \frac{\partial (\kern0.7pt \tilde{p}_0 + \tilde{\phi}_0)}{\partial \tilde{x}} + \frac{\partial^{2} \tilde{u}_0}{\partial \tilde{z}^{2}} + \varPsi \tilde{M}_0 \frac{\partial \tilde{H}}{\partial \tilde{x}} , \end{gather}$$

(2.33c)$$\begin{gather}0 ={-} \frac{\partial (\kern0.7pt \tilde{p}_0 + \tilde{\phi}_0)}{\partial \tilde{y}} + \frac{\partial^{2} \tilde{v}_0}{\partial \tilde{z}^{2}} + \varPsi \tilde{M}_0 \frac{\partial \tilde{H}}{\partial \tilde{y}} , \end{gather}$$

(2.33d)$$\begin{gather}0 ={-} \frac{\partial (\kern0.7pt \tilde{p}_0 + \tilde{\phi}_0)}{\partial \tilde{z}} - \epsilon \mathcal{G} . \end{gather}$$

The subscript 0s in (2.33), which come from the series expansion detailed in (2.23), indicate that leading-order terms have been retained. For brevity, the subscript 0s are omitted in the remainder of the text.

The surface tension was scaled using the maximum surface pressure, $S$,

(2.34)\begin{equation} \tilde{ \gamma} = \frac{\gamma - \gamma_{\textrm{sat}}}{\mathcal{S}} , \quad \mathcal{S} = \gamma_{\mathrm{c}} - \gamma_{\textrm{sat}} , \end{equation}

where $\gamma _{\textrm{sat}}$ is the surface tension of an interface saturated with surfactant and $\gamma _{{\rm c}}$ is the surface tension of a clean interface, i.e. an interface free of surface-active agents. Assuming viscous stresses in the gas phase can be neglected and introducing $\rho ' = \epsilon \rho$, the inner product of the jump momentum balance (2.12) with $\boldsymbol {n}$ when retaining each term at its leading order in $\epsilon$ is

(2.35) \begin{align} \tilde{p} &= \tilde{p}' - \left(\frac{\partial^{2} \tilde{h}}{\partial \tilde{x}^{2}} + \frac{\partial^{2} \tilde{h}}{\partial \tilde{y}^{2}} \right)\left(\epsilon^2 \mathcal{M} \tilde{ \gamma} + \frac{\epsilon^3}{\mathcal{C}}\right) + \tilde{J}_{\rm e}^2 \epsilon^3 Re^{{\star}} (1-\epsilon) \nonumber\\ &\quad- \epsilon^2 \textit{Bq}_{\mathrm{sh}} \left( \frac{\partial \tilde{u}_{{\rm s}}}{\partial \tilde{x}} \left(\frac{\partial^{2} \tilde{h}}{\partial \tilde{x}^{2}} - \frac{\partial^{2} \tilde{h}}{\partial \tilde{y}^{2}}\right) + 2 \frac{\partial^{2} \tilde{h}}{\partial \tilde{x} \partial \tilde{y}} \left( \frac{\partial \tilde{u}_{{\rm s}}}{\partial \tilde{y}} + \frac{\partial \tilde{v}_{{\rm s}}}{\partial \tilde{x}} \right) + \frac{\partial \tilde{v}_{{\rm s}}}{\partial \tilde{y}} \left( \frac{\partial^{2} \tilde{h}}{\partial \tilde{y}^{2}} - \frac{\partial^{2} \tilde{h}}{\partial \tilde{x}^{2}} \right) \right)\nonumber\\ &\quad- \epsilon^2 \textit{Bq}_{\rm d} \left( \frac{\partial^{2} \tilde{h}}{\partial \tilde{x}^{2}} + \frac{\partial^{2} \tilde{h}}{\partial \tilde{y}^{2}} \right) \left( \frac{\partial \tilde{u}_{\rm s}}{\partial \tilde{x}} + \frac{\partial \tilde{v}_{\rm s}}{\partial \tilde{y}} \right) + 2 \epsilon^2 \left(-\frac{\partial \tilde{u}}{\partial \tilde{z}} \frac{\partial \tilde{h}}{\partial \tilde{x}} - \frac{\partial \tilde{v}}{\partial \tilde{z}} \frac{\partial \tilde{h}}{\partial \tilde{y}} + \frac{\partial \tilde{w}}{\partial \tilde{z}} \right)\nonumber\\ &\quad - \frac{1}{2} \frac{M_\mathrm{c}}{H_\mathrm{c}} \epsilon^2 \varPsi \left( - \tilde{M}_x \frac{\partial \tilde{h}}{\partial \tilde{x}} - \tilde{M}_y \frac{\partial \tilde{h}}{\partial \tilde{y}} + \tilde{M}_z \right)^2 \quad \text{at } \varSigma , \end{align}

where $\mathcal {M}$ is the Marangoni number, $\mathcal {C}$ is the capillary number, and $\textit {Bq}_{\mathrm {sh}}$ and $\textit {Bq}_{\rm d}$ are the shear and dilatational Boussinesq numbers, respectively. These dimensionless numbers are defined by

(2.36a–d)\begin{equation} \mathcal{M} = \frac{\mathcal{S} \epsilon}{\eta U} , \quad \mathcal{C} = \frac{\eta U}{\gamma_{\textrm{sat}}} , \quad \textit{Bq}_{\mathrm{sh}} = \frac{\eta_{\rm s} \epsilon}{\eta L} , \quad \textit{Bq}_{\rm d} = \frac{\lambda_{\rm s} \epsilon}{\eta L} . \end{equation}

We used the surrounding gas pressure as the reference pressure, i.e. $\tilde {p}' = 0$. Integrating the leading-order $z$-momentum equation (2.33d) in $z$ and using (2.35) as a boundary condition leads to

(2.37)\begin{equation} \tilde{p} + \tilde{\phi} = \epsilon \mathcal{G} (\tilde{h} - \tilde{z}) - \left(\frac{\partial^{2} \tilde{h}}{\partial \tilde{x}^{2}} + \frac{\partial^{2} \tilde{h}}{\partial \tilde{y}^{2}}\right) \left(\epsilon^2 \mathcal{M} \tilde{\gamma} + \frac{\epsilon^3}{\mathcal{C}}\right) - \tilde{ \varPi} , \end{equation}

where $\tilde { \varPi }$ is the dimensionless disjoining pressure, which is related to the conjoining pressure by $\tilde { \varPi } = - \tilde { \phi }$, and several lower-order terms from (2.35) have been discarded. Having obtained an expression for the pressure, the leading-order $x$-momentum (2.33b) and $y$-momentum (2.33c) equations can be integrated twice in $z$ with the symmetry boundary condition (2.20) and the continuity of tangential velocity across the interface (2.25) to obtain

(2.38a)$$\begin{gather} \tilde{u} = \tilde{u}_{\rm s} + \frac{1}{2}(\tilde{z}^2 - \tilde{h}^2)\left( \frac{\partial (\kern0.7pt \tilde{p} + \tilde{\phi})}{\partial \tilde{x}} - \varPsi \tilde{M} \frac{\partial \tilde{H}}{\partial \tilde{x}} \right) , \end{gather}$$

(2.38b)$$\begin{gather}\tilde{v} = \tilde{v}_{\rm s} + \frac{1}{2}(\tilde{z}^2 - \tilde{h}^2)\left( \frac{\partial (\kern0.7pt \tilde{p} + \tilde{\phi})}{\partial \tilde{y}} - \varPsi \tilde{M} \frac{\partial \tilde{H}}{\partial \tilde{y}} \right) . \end{gather}$$

In deriving (2.38), it was assumed that the magnetisation does not vary over the film depth, as is discussed further in § 2.4. At this point, equations are needed to define the surface velocity components present in (2.38). By taking the inner product of the jump momentum balance (2.12) with $\boldsymbol {t}_x$ and $\boldsymbol {t}_y$, two equations governing $u_{\rm s}$ and $v_{\rm s}$ can be derived. Taking each term in the inner products at leading order in $\epsilon$ results in

(2.39a)$$\begin{gather} \mathcal{M} \frac{\partial \tilde{\gamma}}{\partial \tilde{x}} + \textit{Bq}_{\rm d} \left( \frac{\partial^{2} \tilde{u}_{\rm s}}{\partial \tilde{x}^{2}} + \frac{\partial^{2} \tilde{v}_{\rm s}}{\partial \tilde{x} \partial \tilde{y}} \right) + \textit{Bq}_{\mathrm{sh}} \left( \frac{\partial^{2} \tilde{u}_{\rm s}}{\partial \tilde{x}^{2}} + \frac{\partial^{2} \tilde{u}_{\rm s}}{\partial \tilde{y}^{2}} \right ) = \left.\frac{\partial \tilde{u}}{\partial \tilde{z}}\right\rvert_{\tilde{z} = \tilde{h}}, \end{gather}$$

(2.39b)$$\begin{gather}\mathcal{M} \frac{\partial \tilde{\gamma}}{\partial \tilde{y}} + \textit{Bq}_{\rm d} \left( \frac{\partial^{2} \tilde{u}_{{\rm s}}}{\partial \tilde{y} \partial \tilde{x}} + \frac{\partial^{2} \tilde{v}_{{\rm s}}}{\partial \tilde{y}^{2}} \right) + \textit{Bq}_{\mathrm{sh}} \left( \frac{\partial^{2} \tilde{v}_{\rm s}}{\partial \tilde{x}^{2}} + \frac{\partial^{2} \tilde{v}_{\rm s}}{\partial \tilde{y}^{2}} \right ) = \left. \frac{\partial \tilde{v}}{\partial \tilde{z}} \right\rvert_{\tilde{z} = \tilde{h}} . \end{gather}$$

The velocity gradients in (2.39), which arise from the viscous stress tensor, can be evaluated using the velocities in (2.38). As a result,

(2.40a)$$\begin{gather} \mathcal{M} \frac{\partial \tilde{\gamma}}{\partial \tilde{x}} + \textit{Bq}_{\rm d} \left( \frac{\partial^{2} \tilde{u}_{\rm s}}{\partial \tilde{x}^{2}} + \frac{\partial^{2} \tilde{v}_{\rm s}}{\partial \tilde{x} \partial \tilde{y}} \right) + \textit{Bq}_{\mathrm{sh}} \left( \frac{\partial^{2} \tilde{u}_{\rm s}}{\partial \tilde{x}^{2}} + \frac{\partial^{2} \tilde{u}_{\rm s}}{\partial \tilde{y}^{2}} \right ) = \tilde{h} \left( \frac{\partial (\kern0.7pt \tilde{p} + \tilde{ \phi})}{\partial \tilde{x}} - \varPsi \tilde{M} \frac{\partial \tilde{H}}{\partial \tilde{x}} \right) , \end{gather}$$

(2.40b)$$\begin{gather}\mathcal{M} \frac{\partial \tilde{\gamma}}{\partial \tilde{y}} + \textit{Bq}_{\rm d} \left( \frac{\partial^{2} \tilde{u}_{\rm s}}{\partial \tilde{y} \partial \tilde{x}} + \frac{\partial^{2} \tilde{v}_{\rm s}}{\partial \tilde{y}^{2}} \right) + \textit{Bq}_{\mathrm{sh}} \left( \frac{\partial^{2} \tilde{v}_{\rm s}}{\partial \tilde{x}^{2}} + \frac{\partial^{2} \tilde{v}_{\rm s}}{\partial \tilde{y}^{2}} \right ) = \tilde{h} \left( \frac{\partial (\kern0.7pt \tilde{p} + \tilde{ \phi} )}{\partial \tilde{y}} - \varPsi \tilde{M} \frac{\partial \tilde{H}}{\partial \tilde{y}} \right) . \end{gather}$$

With equations for $\tilde {u}$ (2.38a) and $\tilde {v}$ (2.38b) at hand, the thickness evolution equation can be obtained from (2.29),

(2.41)\begin{align} \frac{\partial \tilde{h}}{\partial \tilde{t}} &={-} \frac{\partial \tilde{u}_{\rm s}\tilde{h}}{\partial \tilde{x}} + \frac{1}{3} \frac{\partial}{\partial \tilde{x}}\left( \tilde{h}^3 \left( \frac{\partial (\kern0.7pt \tilde{p} + \tilde{\phi})}{\partial \tilde{x}} - \varPsi \tilde{M} \frac{\partial \tilde{H}}{\partial \tilde{x}} \right) \right)\nonumber\\ &\quad - \frac{\partial \tilde{v}_{\rm s}\tilde{h}}{\partial \tilde{y}} + \frac{1}{3} \frac{\partial}{\partial \tilde{y}}\left( \tilde{h}^3 \left( \frac{\partial (\kern0.7pt \tilde{p} + \tilde{\phi})}{\partial \tilde{y}} - \varPsi \tilde{M} \frac{\partial \tilde{H}}{\partial \tilde{y}} \right) \right) - \tilde{J}_{\rm e} . \end{align}

There is no explicit $z$ dependence in the derived thickness evolution equation because of the thin film approximation. Consequently, the derived equations can be expressed compactly by introducing a dimensionless planar gradient operator defined by $\tilde {\boldsymbol {\nabla }}_{\rm p} = ({\partial }/{\partial \tilde {x}}) \boldsymbol {e}_{x} + ({\partial }/{\partial \tilde {y}}) \boldsymbol {e}_{y}$, which is a dimensionless two-dimensional version of the three-dimensional gradient operator, $\boldsymbol {\nabla }$. We denote the dimensionless planar divergence operator by $\tilde {\boldsymbol {\nabla }}_{\rm p} \boldsymbol {\cdot }$ and the dimensionless planar Laplace operator (for both scalar and vector fields) by $\tilde {\nabla }_{\rm p}^{2} \,$. At this stage, the centre surface of the film can be introduced such that the dividing surface is located at $z = h(x,y,t) + C(x,y)$. With the dimensionless centre surface given by $\tilde {C} = C / \mathcal {H}$, the equations governing the thickness evolution, written using planar operations, are

(2.42a)$$\begin{gather} \frac{\partial \tilde{h}}{\partial \tilde{t}} ={-} \tilde{\boldsymbol{\nabla}}_{\rm p} \boldsymbol{\cdot} ( \tilde{\boldsymbol{V}}_{\rm s} \tilde{h} ) + \frac{1}{3} \tilde{\boldsymbol{\nabla}}_{\rm p} \boldsymbol{\cdot} ( \tilde{h}^{3} ( \tilde{\boldsymbol{\nabla}}_{\rm p} (\kern0.7pt \tilde{p} + \tilde{\phi}) - \varPsi \tilde{M} \tilde{\boldsymbol{\nabla}}_{\rm p} \tilde{H} ) ) - \tilde{J}_{\rm e} , \end{gather}$$

(2.42b)$$\begin{gather}\tilde{p} + \tilde{\phi} = \epsilon \mathcal{G} (\tilde{h} + \tilde{C} - \tilde{z}) - \left(\epsilon^2 \mathcal{M} \tilde{\gamma} + \frac{\epsilon^3}{\mathcal{C}}\right)( \tilde{\nabla}_{\rm p}^{2} \tilde{h} + \tilde{\nabla}_{\rm p}^{2} \tilde{C} )- \tilde{\varPi} , \end{gather}$$

(2.42c)$$\begin{gather}\mathcal{M} \tilde{\boldsymbol{\nabla}}_{\rm p} \tilde{\gamma} + \textit{Bq}_{\rm d} \tilde{\boldsymbol{\nabla}}_{\rm p} \boldsymbol{\cdot} ((\tilde{\boldsymbol{\nabla}}_{\rm p} \boldsymbol{\cdot} \tilde{\boldsymbol{V}_{\rm s}}) \boldsymbol{\mathsf{I}}) + \textit{Bq}_{\mathrm{sh}} \tilde{\nabla}_{\rm p}^{2} \tilde{\boldsymbol{V}}_{\rm s} = \tilde{h} ( \tilde{\boldsymbol{\nabla}}_{\rm p} (\kern0.7pt \tilde{p} + \tilde{\phi} ) - \varPsi \tilde{M} \tilde{\boldsymbol{\nabla}}_{\rm p} \tilde{H} ) . \end{gather}$$

Equation (2.42c) arises by writing two scalar equations, (2.40a) and (2.40b), as one vector equation, where $\tilde {\boldsymbol {V}}_{\rm s} = \tilde {u}_{\rm s} \boldsymbol {e}_{x} + \tilde {v}_{\rm s} \boldsymbol {e}_{y}$. In order for the thin film approximation to remain valid with a curved centre surface, the curvature of the centre surface, $\kappa _{\textrm{cs}}$, must be everywhere much less than the reciprocal of the characteristic thickness scale, i.e. $\kappa _{\textrm{cs}} \ll 1/\mathcal {H}$ (O'Brien & Schwartz Reference O'Brien and Schwartz2002).

The magnetisation is a function of the molar concentration of magnetite NPs, $c$, and the surface tension depends on the surface excess concentration of surfactant, $\varGamma$. Equations characterising $c$ and $\varGamma$ are derived in §§ 2.4 and 2.5, respectively. These equations need to be solved together with (2.42a), (2.42b) and (2.42c) to determine the time evolution of film thickness.

2.4. Magnetite NP transport

By approximating the molar flux of magnetite NPs as the sum of fluxes due to advective transport, magnetic forcing and diffusion, the transport equation for the molar concentration of magnetite NPs, $c$, in the core liquid of a magnetic soap film can be expressed as (Pshenichnikov et al. Reference Pshenichnikov, Elfimova and Ivanov2011; Nacev et al. Reference Nacev, Komaee, Sarwar, Probst, Kim, Emmert-Buck and Shapiro2012)

(2.43)\begin{equation} \frac{\partial{c}}{\partial{t}} + \boldsymbol{\nabla} \boldsymbol{\cdot} (c \boldsymbol{V}) + \boldsymbol{\nabla} \boldsymbol{\cdot} \left(c b_c \frac{\boldsymbol{f}_{\rm m}}{n}\right) = \boldsymbol{\nabla} \boldsymbol{\cdot} (D_c \boldsymbol{\nabla} c) , \end{equation}

where $b_c$, $n$ and $D_c$ are the mobility, number concentration and gradient diffusion coefficient of the magnetite NPs, respectively. For a dilute dispersion of magnetite NPs in a liquid, the mobility and gradient diffusion coefficient can be approximated by (Batchelor Reference Batchelor1983)

(2.44a,b)\begin{equation} b_c = b_0 (1 - 6.55 \varphi) \quad \text{and} \quad D_c = D_0 (1 + 1.45 \varphi) , \quad \varphi \ll 1 , \end{equation}

where $\varphi$ is the volume concentration of magnetite NPs, and $b_0$ and $D_0$ are the mobility and Brownian diffusion coefficient of an independent sphere, respectively. For a sphere of radius $a$ (Batchelor Reference Batchelor1983),

(2.45)\begin{equation} D_0 = \frac{k_{\rm B} \mathcal{T}}{6{\rm \pi}\eta a} , \end{equation}

where $k_{\rm B}$ is the Boltzmann constant and $\mathcal {T}$ is the absolute temperature. The mobility $b_0$ is related to $D_0$ through $D_0 = b_0 k_{\rm B} \mathcal {T}$ (Peskir Reference Peskir2003). Here, the mobility and gradient diffusion coefficient of magnetite NPs in a thin film have been modelled using equations derived for particles in a bulk liquid, i.e. not thin films of liquid. This is a good approximation when the film thickness dominates the size of the particles dispersed in the film (Vivek & Weeks Reference Vivek and Weeks2015; Kim et al. Reference Kim, Ribbe, Russell and Hoagland2016). For particles with size of the order of the film thickness, Prasad & Weeks (Reference Prasad and Weeks2009) found that diffusion occurred faster in soap films compared with diffusion in a bulk soap film solution. As a result, the mobility and gradient diffusion coefficient are expected to be functions of the film thickness. It is anticipated that the thickness dependence may be derived by considering the Trapeznikov approximation (Prasad & Weeks Reference Prasad and Weeks2009; Vivek & Weeks Reference Vivek and Weeks2015); however, this is not considered here as the film thickness was considerably larger than the particle dimensions throughout most of the film thinning in the experiments.

For a ferrofluid containing a dilute concentration of monodisperse magnetite NPs, the magnetisation of the ferrofluid varies with the magnetic field intensity according to (Rosensweig Reference Rosensweig2013)

(2.46)\begin{equation} \boldsymbol{M} = M_{\textrm{sat}} \mathcal{L} (\xi) \frac{\boldsymbol{H}}{H} , \quad \xi = \frac{\mu_{0} m H}{k_{\rm B} \mathcal{T}}, \end{equation}

where $M_{\textrm{sat}}$ is the saturation magnetisation of the ferrofluid, $\mathcal {L}$ is the Langevin function ($\mathcal {L}(\xi ) = \coth \xi - 1 / \xi$), $\xi$ is the Langevin parameter and $m$ is the magnetic dipole moment of a magnetite NP. Since $M_{\textrm{sat}} = nm$, $M = nm \mathcal {L}(\xi )$.

It is evident that the magnetite NP transport equation (2.43) is three dimensional, whilst the equation governing the film thickness (2.42a) is two dimensional. The three-dimensional transport equation can be integrated over the film thickness by assuming a concentration profile over the film depth (Nierstrasz & Frens Reference Nierstrasz and Frens1999; Huybrechts, Villaret & Hervouet Reference Huybrechts, Villaret and Hervouet2010) to derive a two-dimensional transport equation for the magnetite NPs. We assumed a uniform concentration of magnetite NPs over the film depth since Brownian motion and the surfactant coating around the magnetite NPs act to uniformly disperse the NPs in the core liquid (Rosensweig Reference Rosensweig1987). Upon integrating (2.43) in $z$ and introducing the linear approximations (2.44), the magnetic force (2.5) and the magnetisation relation (2.46), the two-dimensional transport equation for the dimensionless molar concentration of magnetite NPs is

(2.47)\begin{align} &\frac{\partial \tilde{c}}{\partial \tilde{t}} + \frac{\tilde{\boldsymbol{Q}}}{\tilde{h}} \boldsymbol{\cdot} \tilde{\boldsymbol{\nabla}}_{\rm p} \tilde{c} + \frac{\alpha}{Pe_c} \tilde{\boldsymbol{\nabla}}_{\rm p} \boldsymbol{\cdot} (\tilde{c} ( 1 - 6.55 \varphi_{{\rm c}} \tilde{c}) \mathcal{L} (\alpha \tilde{H}) \tilde{\boldsymbol{\nabla}}_{\rm p} \tilde{H} )\nonumber\\ &\quad =\frac{1}{Pe_c} \tilde{\boldsymbol{\nabla}}_{\rm p} \boldsymbol{\cdot} ((1 + 1.45 \varphi_{{\rm c}} \tilde{c}) \tilde{\boldsymbol{\nabla}}_{\rm p} \tilde{c}) , \end{align}

where $\tilde {c} = c / c_{{\rm c}}$, $\tilde {\boldsymbol {Q}} = \tilde {Q}_x \boldsymbol {e}_{x} + \tilde {Q}_y \boldsymbol {e}_{y}$, $\alpha$ is a Langevin-type parameter, $Pe_c$ is a Péclet number for magnetite NP transport and $\varphi _{{\rm c}}$ is the volume concentration for $c = c_{{\rm c}}$. The introduced dimensionless numbers are given by

(2.48a–c)\begin{equation} \alpha = \frac{\mu_0 m H_\mathrm{c}}{k_{\rm B} \mathcal{T}} , \quad \varphi_{{\rm c}} = c_{{\rm c}} \mathcal{N}_{\rm A} \mathcal{V}_{\mathrm{p}} , \quad Pe_c = \frac{U L}{D_0} , \end{equation}

where $\mathcal {N}_{\rm A}$ is Avogadro's constant and $\mathcal {V}_{\mathrm {p}}$ is the volume of each NP including the coating.

2.4.1. Coupling between film thickness and magnetite NP transport

In order to couple the film thickness equation (2.42a) with the magnetite NP transport equation (2.47), a relationship for $\tilde {M}$ as a function of $\tilde {c}$ is needed. From $M = n m \mathcal {L} (\xi )$,

(2.49)\begin{equation} \tilde{M} = \tilde{c} \frac{c_\mathrm{c}}{M_\mathrm{c}} \mathcal{N}_{\rm A} m \mathcal{L} (\alpha \tilde{H}) . \end{equation}

Taking $M_{{\rm c}}$ to be the magnetisation when $c = c_{{\rm c}}$ and $H = H_{{\rm c}}$,

(2.50)\begin{equation} M_\mathrm{c} = c_\mathrm{c} \mathcal{N}_{\rm A} m \mathcal{L} (\alpha) . \end{equation}

Injecting (2.50) into (2.49) results in the required relationship for $\tilde {M}$,

(2.51)\begin{equation} \tilde{M} = \tilde{c} \frac{\mathcal{L} (\alpha \tilde{H})}{\mathcal{L} (\alpha)} . \end{equation}

2.5. Surfactant transport

The dimensionless surface tension in the pressure equation (2.42b) depends on the surfactant dynamics in the film. We assume that the surfactants stabilising each magnetic soap film are insoluble and that they form a Langmuir monolayer at each film interface with the relationship between the surface tension and surface concentration of surfactant given by the Langmuir isotherm (Manikantan & Squires Reference Manikantan and Squires2020),

(2.52)\begin{equation} \gamma = \gamma_{{\rm c}} + \varGamma_{\textrm{sat}} R \mathcal{T} \ln \left(1 - \frac{\varGamma}{\varGamma_{\textrm{sat}}}\right) , \end{equation}

where $R$ is the ideal gas constant, $\varGamma _{\textrm{sat}}$ is the surface surfactant concentration for a saturated interface in units of mol m$^{-2}$ and $\ln$ is the natural logarithm. With $\tilde { \varGamma } = \varGamma / \varGamma _{\textrm{sat}}$, the Langmuir isotherm (2.52) in dimensionless form is

(2.53)\begin{equation} \tilde{ \gamma} = 1 + \varLambda \ln (1 - \tilde{ \varGamma}) , \quad \varLambda = \frac{\varGamma_{\textrm{sat}} R \mathcal{T}}{\mathcal{S}} . \end{equation}

The transport equation for the surface concentration of surfactant along a deforming surface is (Stone Reference Stone1990)

(2.54)\begin{equation} \left(\frac{\partial \varGamma}{\partial t}\right)_{\rm n} + \boldsymbol{\nabla}_{\rm s} \boldsymbol{\cdot} (\varGamma \boldsymbol{\mathsf{P}}_{\rm s} \boldsymbol{\cdot} \boldsymbol{V}_{{\rm s}}) + \varGamma (\boldsymbol{\nabla}_{\rm s} \boldsymbol{\cdot} \boldsymbol{n})(\boldsymbol{V}_{{\rm s}} \boldsymbol{\cdot} \boldsymbol{n}) = D_{\rm s} \nabla_{\rm s}^{2} \varGamma , \end{equation}

where $D_{\rm s}$ is the surface diffusion coefficient for surfactant transport and $\nabla _{\rm s}^{2}$ is the surface Laplace operator. This equation assumes that there is no transport of surfactant between the core liquid and the surface and that the diffusive flux of surfactant along the surface is governed by

(2.55)\begin{equation} \boldsymbol{J}_{\mathrm{D}} ={-} D_{\rm s} \boldsymbol{\nabla}_{\rm s} \varGamma \end{equation}

with $D_{\rm s}$ uniform over the surface. This diffusion relation is not consistent with the Langmuir isotherm; therefore, the following diffusive flux relation that is consistent with the Langmuir isotherm was used (Manikantan & Squires Reference Manikantan and Squires2020):

(2.56)\begin{equation} \boldsymbol{J}_{\mathrm{D}} ={-} \frac{D_{\rm s} \varGamma_{\textrm{sat}}}{\varGamma_{\textrm{sat}} - \varGamma} \boldsymbol{\nabla}_{\rm s} \varGamma .\end{equation}

Here $D_{\rm s}$ was assumed to be uniform over the surface and constant in time. The partial time derivative in the surfactant transport equation (2.54) should be interpreted as the change in $\varGamma$ when following the interface in a direction normal to the original interface (Wong, Rumschitzki & Maldarelli Reference Wong, Rumschitzki and Maldarelli1996; Pereira et al. Reference Pereira, Trevelyan, Thiele and Kalliadasis2007). Each position on the dividing surface can be specified by two parameters known as the surface coordinates (Slattery et al. Reference Slattery, Sagis and Oh2007), and it is easier to work in terms of partial time derivatives when following the evolution of the interface at fixed surface coordinates. We can conveniently use $x$ and $y$ as the surface coordinates since the position vector of any point on the dividing surface is

(2.57)\begin{equation} \boldsymbol{r}_{\rm s} = x \boldsymbol{e}_{x} + y \boldsymbol{e}_{y} + (h(x, y, t) + C(x, y)) \boldsymbol{e}_{z} . \end{equation}

Wong et al. (Reference Wong, Rumschitzki and Maldarelli1996) and Cermelli, Fried & Gurtin (Reference Cermelli, Fried and Gurtin2005) provide a relation for the partial time derivative present in the transport equation (2.54) in terms of the partial time derivative when keeping the surface coordinates fixed, which we denote by $(\frac {\partial \varGamma }{\partial t})_{\circ }$,

(2.58)\begin{equation} \left(\frac{\partial \varGamma}{\partial t}\right)_{\rm n} = \left(\frac{\partial \varGamma}{\partial t}\right)_{{\circ}} - \boldsymbol{U}_{{\rm s}} \boldsymbol{\cdot} \boldsymbol{\nabla}_{\rm s} \varGamma , \end{equation}

where $\boldsymbol {U}_{{\rm s}}$ is the rate of change of spatial position following a point on the surface with fixed surface coordinates, i.e.

(2.59)\begin{equation} \boldsymbol{U}_{\rm s} = \left(\frac{\partial \boldsymbol{r}_{\rm s}}{\partial t}\right)_{{\circ}} . \end{equation}

Following from (2.57), $\boldsymbol {U}_{\rm s} = (\partial h/\partial t) \boldsymbol {e}_{z}$. With (2.56) and (2.58), the surfactant transport equation (2.54) becomes

(2.60) \begin{align} &\left(\frac{\partial \varGamma}{\partial t}\right)_{{\circ}} - \boldsymbol{U}_{\rm s} \boldsymbol{\cdot} \boldsymbol{\nabla}_{\rm s} \varGamma + \boldsymbol{\nabla}_{\rm s} \boldsymbol{\cdot} (\varGamma \boldsymbol{\mathsf{P}}_{\rm s} \boldsymbol{\cdot} \boldsymbol{V}_{\rm s}) + \varGamma (\boldsymbol{\nabla}_{\rm s} \boldsymbol{\cdot} \boldsymbol{n})(\boldsymbol{V}_{\rm s} \boldsymbol{\cdot} \boldsymbol{n}) \nonumber\\ &\quad = D_{\rm s} \boldsymbol{\nabla}_{\rm s} \boldsymbol{\cdot} \left(\frac{\varGamma_{\textrm{sat}}}{\varGamma_{\textrm{sat}}-\varGamma} \boldsymbol{\nabla}_{\rm s} \varGamma \right) , \end{align}

which in dimensionless form and at leading order in $\epsilon$ is

(2.61)\begin{equation} \left(\frac{\partial \tilde{\varGamma}}{\partial \tilde{t}}\right)_{{\circ}} + \tilde{\boldsymbol{\nabla}}_{\rm p} \boldsymbol{\cdot} ( \tilde{\boldsymbol{V}}_{\rm s} \tilde{\varGamma} ) = \frac{1}{Pe_{\rm s}} \tilde{\boldsymbol{\nabla}}_{\rm p} \boldsymbol{\cdot} \left( \frac{1}{1 - \tilde{ \varGamma}} \tilde{\boldsymbol{\nabla}}_{\rm p} \tilde{ \varGamma} \right) , \end{equation}

where $Pe_{\rm s}$ is the Péclet number for surfactant transport along the dividing surface. This Péclet number is given by

(2.62)\begin{equation} Pe_{\rm s} = \frac{U L}{D_{\rm s}} . \end{equation}

2.6. Closures

In order to solve for the film thickness, closures are needed for the evaporation rate and disjoining pressure.

2.6.2. Disjoining pressure

The Derjaguin–Landau–Verwey–Overbeek (DLVO) theory was used for the disjoining pressure in this study (see Ninham Reference Ninham1999 for a list of assumptions inherent in DLVO theory). With the DLVO theory, the disjoining pressure is expressed as the sum of attractive van der Waals forces and repulsive electric double layer forces. Under the assumption of identically charged planar surfaces with a small and constant surface potential, the disjoining pressure is given by the relation (Bergeron Reference Bergeron1999; Matsubara et al. Reference Matsubara, Kimura, Miyao, Shin and Ikeda2021; Agmo Hernández Reference Agmo Hernández2023)

(2.63)\begin{equation} \varPi ={-} \frac{A}{6 {\rm \pi}(2h)^3} + B {\rm e}^{- \zeta 2h} , \end{equation}

where $A$ is the Hamaker constant, $B$ is a parameter that depends on the interface potential and $\zeta ^{-1}$ is the Debye length. In dimensionless form, (2.63) is

(2.64)\begin{equation} \tilde{ \varPi} ={-} \frac{\mathcal{A}}{\tilde{h}^3} + \mathcal{B} \,{\rm e}^{-\mathcal{D} \tilde{h}} , \end{equation}

where the introduced dimensionless numbers are

(2.65a–c)\begin{equation} \mathcal{A} = \frac{A}{48{\rm \pi} \mathcal{H}^3 P} , \quad \mathcal{B} = \frac{B}{P} , \quad \mathcal{D} = 2 \zeta \mathcal{H} . \end{equation}

We used constant values for $\mathcal {A}$, $\mathcal {B}$ and $\mathcal {D}$ throughout the simulations, which does not account for the spatiotemporally varying interfacial surfactant concentration. Note that DLVO theory does not have general applicability, for example, it would not be effective at predicting the equilibrium thickness of Newton black films (Ochoa et al. Reference Ochoa, Gao, Srivastava and Sharma2021). Additional interactions, such as steric forces or solvation forces, may need to be accounted for to more accurately represent the disjoining pressure for films of interest (Bergeron Reference Bergeron1999).

2.7. Solving

The thickness field is governed by the system of coupled partial differential equations (PDEs) defined by (2.42a–c), (2.47) and (2.61), where the variables that must be solved for are $(\tilde {h}, \tilde {p}, \tilde { \boldsymbol {V}}_{{\rm s}}, \tilde {c}, \tilde { \varGamma })$. The Galerkin finite element method was used to solve this coupled system with the addition of suitable initial and boundary conditions and the relations for the magnetisation (2.51), surface tension (2.53) and disjoining pressure (2.64). With (2.42b) in (2.42a), the thickness evolution equation is a nonlinear fourth-order parabolic equation (Liu, Peco & Dolbow Reference Liu, Peco and Dolbow2019). A fourth-order PDE of this type can be solved with mixed finite element methods (Wells, Kuhl & Garikipati Reference Wells, Kuhl and Garikipati2006; Liu et al. Reference Liu, Peco and Dolbow2019; Keita, Beljadid & Bourgault Reference Keita, Beljadid and Bourgault2021), which involve solving the fourth-order equation as a system of two second-order equations. This allows for the use of standard Lagrange finite element basis functions to approximate the unknowns. It was natural to solve for $\tilde {p} + \tilde { \phi }$ (2.42b) as one of the second-order equations.

Uniform initial conditions were used for each of the five variables. In all simulations, the initial pressure was set to 0 and the initial surface velocity was set to $\boldsymbol {0}$. The initial conditions for $\tilde {h}$, $\tilde {c}$ and $\tilde { \varGamma }$, which we denote by $\tilde {h}_{\rm i}$, $\tilde {c}_{\rm i}$ and $\tilde { \varGamma }_{\rm i}$, are given in table 1 for each simulation.

Table 1. The dimensionless numbers used for the governing equations, boundary conditions and initial conditions are provided for the simulations presented in each section. The number of cells, $N_{\textrm{cell}}$, used to discretise the domain to the nearest 100 cells is also provided. Section 3.1(i) refers to the immobile soap film and § 3.1(p) refers to the partially mobile soap films simulated in § 3.1. Section 3.2(u) refers to the simulated magnetic soap films with a uniform concentration of magnetite NPs for two boundary conditions, and § 3.2(t) refers to the same but when considering the transport of magnetite NPs. The $\leftarrow$ symbol indicates that the simulation parameter was the same as the value for the simulation in the column to the left, e.g. the same value was used for the capillary number, $\mathcal {C}$, for all simulations. A hyphen symbolises that a parameter does not affect the results of a simulation, e.g. the value of $\varPsi$ has no effect when there are no magnetite NPs in the film ($\tilde {c}_{\rm i} = 0.0$).

Five boundary conditions are needed on $\partial \varOmega$, one corresponding to each of the five second-order PDEs. As discussed in the introduction, a net flow of liquid is expected between the meniscus and the thin film due to capillary suction and marginal regeneration; therefore, we did not employ a no-flux boundary condition. Instead, since a larger force acts on thicker film due to this suction compared with thinner film (Overbeek Reference Overbeek1960; Mysels Reference Mysels1968), we imposed a boundary condition relating the flux of liquid out of the thin film due to capillary suction to the local film thickness at the boundary,

(2.66a,b)\begin{equation} - \frac{1}{3} \tilde{h}^3 \tilde{\boldsymbol{\nabla}}_{\rm p} \tilde{p}_{{\rm c}} \boldsymbol{\cdot} \boldsymbol{n}_{\rm d} = \mathcal{P} \tilde{h}^3 \quad \text{on } \partial\varOmega , \quad \tilde{p}_{{\rm c}} ={-} \left(\epsilon^2 \mathcal{M} \tilde{\gamma} + \frac{\epsilon^3}{\mathcal{C}}\right) \tilde{\nabla}_{\rm p}^{2} \tilde{h} , \end{equation}

where $\mathcal {P}$ governs the strength of capillary suction and $\tilde {p}_{\rm c}$ is the dimensionless capillary pressure. It is unknown at present how the suction strength depends on the applied magnetic field due to the influence of the applied magnetic field on the meniscus shape (Rosensweig et al. Reference Rosensweig, Elborai, Lee and Zahn2005; Singh et al. Reference Singh, Singh and Bhaumik2024); this could be an interesting area for future research. For some of the simulations in which a magnetic field was applied, we also investigated the boundary condition

(2.67)\begin{equation} -\tfrac 1 3 \tilde{h}^3 (\tilde{\boldsymbol{\nabla}}_{\rm p} \tilde{p} - \varPsi \tilde{M} \tilde{\boldsymbol{\nabla}}_{\rm p} \tilde{H} ) \boldsymbol{\cdot} \boldsymbol{n}_{\rm d} = \mathcal{Q} \tilde{h}^3 \quad \text{on } \partial\varOmega , \end{equation}

which imposes a limit on the drainage flux of core liquid out of the thin film (without disjoining pressure contributions), determined by the parameter $\mathcal {Q}$.

A Neumann boundary condition was used for the film thickness to allow the film thickness to change with time at the boundary of the thin film. The Neumann condition sets $\tilde {\boldsymbol {\nabla }}_{\rm p} \tilde {h} \boldsymbol {\cdot } \boldsymbol {n}_{\rm d}$ on $\partial \varOmega$, which we denote by ${\partial \tilde {h}}/{\partial n_{\rm d}}$.

The divergence of the tensor $\textit {Bq}_{\rm d} (\tilde {\boldsymbol {\nabla }}_{\rm p} \boldsymbol {\cdot } \tilde {\boldsymbol {V}}_{\rm s}) \boldsymbol{\mathsf{I}} + \textit {Bq}_{\mathrm {sh}} \tilde {\boldsymbol {\nabla }}_{\rm p} \tilde {\boldsymbol {V}}_{\rm s}$ appears in (2.42c), and the weak variational form features this tensor operating on $\boldsymbol {n}_{\rm d}$. We chose to adopt the natural boundary condition

(2.68)\begin{equation} \textit{Bq}_{\rm d} (\tilde{\boldsymbol{\nabla}}_{\rm p} \boldsymbol{\cdot} \tilde{\boldsymbol{V}}_{\rm s} ) \boldsymbol{n}_{\rm d} + \textit{Bq}_{\mathrm{sh}} \tilde{\boldsymbol{\nabla}}_{\rm p} \tilde{\boldsymbol{V}}_{\rm s} \boldsymbol{\cdot} \boldsymbol{n}_{\rm d} = \boldsymbol{0} \quad \text{on } \partial \varOmega . \end{equation}

A number of Neumann and Dirichlet boundary conditions were investigated for $\tilde {c}$; however, the most physically consistent behaviour was observed with the boundary condition

(2.69)\begin{equation} \alpha \tilde{c} (1 - 6.55 \varphi_{{\rm c}} \tilde{c}) \mathcal{L} (\alpha \tilde{H}) \tilde{\boldsymbol{\nabla}}_{\rm p} \tilde{H} \boldsymbol{\cdot} \boldsymbol{n}_{\rm d} = (1 + 1.45 \varphi_{{\rm c}} \tilde{c}) \tilde{\boldsymbol{\nabla}}_{\rm p} \tilde{c} \boldsymbol{\cdot} \boldsymbol{n}_{\rm d} \quad \text{on } \partial \varOmega , \end{equation}

which is a balance of the flux of $\tilde {c}$ out of the thin film between magnetic forcing and diffusion. With this boundary condition, the concentration of magnetite NPs averaged over the thin film can change with time due to the advective flux in (2.47).

Nierstrasz & Frens (Reference Nierstrasz and Frens1999) argued that with marginal regeneration, the surfactant that flows out of the thin film during drainage is replaced by an equal amount of surfactant from saturated patches of fluid that are extracted from the film boundary. Consequently, we used a boundary condition that keeps the amount of surfactant covering the thin film constant in time,

(2.70)\begin{equation} \tilde{ \varGamma} \tilde{\boldsymbol{V}}_{\rm s} \boldsymbol{\cdot} \boldsymbol{n}_{\rm d} = \frac{1}{Pe_{\rm s}} \frac{1}{1 - \tilde{ \varGamma}} \tilde{\boldsymbol{\nabla}}_{\rm p} \tilde{ \varGamma} \boldsymbol{\cdot} \boldsymbol{n}_{\rm d} \quad \text{on } \partial\varOmega , \end{equation}

which is a no-flux surfactant boundary condition.

We use the following compact notation to express the problem of finding the time evolution of film thickness in weak form (Langtangen & Logg Reference Langtangen and Logg2017):

(2.71a,b)\begin{equation} \langle \boldsymbol{\nabla} h, \boldsymbol{\nabla} k \rangle = \int_{\varOmega} \boldsymbol{\nabla} h \boldsymbol{\cdot} \boldsymbol{\nabla} k \, \mathrm{d}\varOmega , \quad \langle h, k \rangle_{\partial\varOmega} = \int_{\partial \varOmega} h k \, \mathrm{d}s . \end{equation}

The mixed weak form of the problem, with the boundary conditions inserted into the surface integral terms and the initial conditions not repeated here for brevity, is to find $(\tilde {h}, \tilde {p}, \tilde { \boldsymbol {V}}_{{\rm s}}, \tilde {c}, \tilde { \varGamma }) \in K^{(0)} \times {K^{(1)}} \times {{ \boldsymbol {K}}^{(2)}} \times {K^{(3)}} \times {K^{(4)}}$ such that

(2.72a)\begin{align} &\left\langle \frac{\partial \tilde{h}}{\partial \tilde{t}}, {k^{(0)}} \right\rangle + \langle \tilde{\boldsymbol{\nabla}}_{\rm p} \boldsymbol{\cdot} ( \tilde{\boldsymbol{V}}_{\rm s} \tilde{h} ), {k^{(0)}} \rangle + \frac 13 \langle \tilde{h}^3 \tilde{\boldsymbol{\nabla}}_{\rm p} (\kern0.7pt \tilde{p} + \tilde{ \phi}), \tilde{\boldsymbol{\nabla}}_{\rm p} {k^{(0)}} \rangle + \langle \tilde{J}_{\rm e}, {k^{(0)}} \rangle \nonumber\\ &\quad + \frac 13 \varPsi \langle \tilde{\boldsymbol{\nabla}}_{\rm p} \boldsymbol{\cdot} (\tilde{h}^3 \tilde{M} \tilde{\boldsymbol{\nabla}}_{\rm p} \tilde{H}), {k^{(0)}} \rangle - \frac 13 \epsilon \mathcal{G} \langle \tilde{h}^3 \tilde{\boldsymbol{\nabla}}_{\rm p} (\tilde{h} + \tilde{C}) \boldsymbol{\cdot} \boldsymbol{n}_{\rm d}, {k^{(0)}} \rangle_{\partial\varOmega} \nonumber\\ &\quad + \langle \mathcal{P} \tilde{h}^3, {k^{(0)}} \rangle_{\partial\varOmega} + \frac{1}{3} \langle \tilde{h}^3 \tilde{\boldsymbol{\nabla}}_{\rm p} \tilde{\varPi} \boldsymbol{\cdot} \boldsymbol{n}_{\rm d}, {k^{(0)}} \rangle_{\partial\varOmega} = 0 \quad \forall {k^{(0)}} \in {K^{(0)}} , \end{align}

(2.72b)\begin{align} &\langle \tilde{p} + \tilde{ \phi}, {k^{(1)}} \rangle - \epsilon \mathcal{G} \langle \tilde{h} + \tilde{C}, {k^{(1)}} \rangle - \left\langle \left(\epsilon^2 \mathcal{M} \tilde{ \gamma} + \frac{\epsilon^3}{\mathcal{C}}\right) \tilde{\boldsymbol{\nabla}}_{\rm p} \tilde{h} , \tilde{\boldsymbol{\nabla}}_{\rm p} {k^{(1)}} \right\rangle\nonumber\\ &\quad + \langle \tilde{ \varPi}, {k^{(1)}} \rangle + \left\langle \left(\epsilon^2 \mathcal{M} \tilde{ \gamma} + \frac{\epsilon^3}{\mathcal{C}}\right) \frac{\partial \tilde{h}}{\partial n_{\rm d}}, {k^{(1)}} \right\rangle_{\partial\varOmega} = 0 \quad \forall {k^{(1)}} \in {K^{(1)}} , \end{align}

(2.72c)\begin{align} &\mathcal{M} \langle \tilde{\boldsymbol{\nabla}}_{\rm p} \tilde{ \gamma}, {{ \boldsymbol {k}}^{(2)}} \rangle - \int_{\varOmega} (\textit{Bq}_{\rm d} (\tilde{\boldsymbol{\nabla}}_{\rm p} \boldsymbol{\cdot} \tilde{\boldsymbol{V}}_{\rm s}) \boldsymbol{\mathsf{I}} + \textit{Bq}_{\mathrm{sh}} \tilde{\boldsymbol{\nabla}}_{\rm p} \tilde{\boldsymbol{V}}_{\rm s}) \boldsymbol{:} \tilde{\boldsymbol{\nabla}}_{\rm p} {{ \boldsymbol {k}}^{(2)}} \, \mathrm{d}\varOmega\nonumber\\ &\quad- \langle \tilde{h} (\tilde{\boldsymbol{\nabla}}_{\rm p} (\kern0.7pt \tilde{p} + \tilde{ \phi}) - \varPsi \tilde{M} \tilde{\boldsymbol{\nabla}}_{\rm p} \tilde{H}), {{ \boldsymbol {k}}^{(2)}} \rangle = 0 \quad \forall {{ \boldsymbol {k}}^{(2)}} \in {{ \boldsymbol {K}}^{2}} , \end{align}

(2.72d)\begin{align} &\left\langle \frac{\partial \tilde{c}}{\partial \tilde{t}}, {k^{(3)}} \right\rangle + \left\langle \frac{\tilde{\boldsymbol{Q}}}{\tilde{h}} \boldsymbol{\cdot} \tilde{\boldsymbol{\nabla}}_{\rm p} \tilde{c}, {k^{(3)}} \right\rangle + \frac{1}{Pe_c} \langle (1 + 1.45 \varphi_{{\rm c}} \tilde{c}) \tilde{\boldsymbol{\nabla}}_{\rm p} \tilde{c}, \tilde{\boldsymbol{\nabla}}_{\rm p} {k^{(3)}} \rangle \nonumber\\ &\quad - \frac{\alpha}{Pe_c} \langle \tilde{c} ( 1 - 6.55 \varphi_{{\rm c}} \tilde{c}) \mathcal{L} (\alpha \tilde{H}) \tilde{\boldsymbol{\nabla}}_{\rm p} \tilde{H} , \tilde{\boldsymbol{\nabla}}_{\rm p} {k^{(3)}} \rangle = 0 \quad \forall {k^{(3)}} \in {K^{(3)}} , \end{align}

(2.72e)\begin{align} &\left\langle \left(\frac{\partial \tilde{\varGamma}}{\partial \tilde{t}}\right)_{{\circ}}, {k^{(4)}} \right\rangle - \langle \tilde{\varGamma} \tilde{\boldsymbol{V}}_{\rm s}, \tilde{\boldsymbol{\nabla}}_{\rm p}{k^{(4)}} \rangle + \frac{1}{Pe_{\rm s}} \left\langle \frac{1}{1 - \tilde{ \varGamma}} \tilde{\boldsymbol{\nabla}}_{\rm p} \tilde{ \varGamma}, \tilde{\boldsymbol{\nabla}}_{\rm p} {k^{(4)}} \right\rangle = 0 \quad \forall {k^{(4)}} \in {K^{(4)}} , \end{align}

where $\boldsymbol {:}$ is the Frobenius inner product, $k^{(i)}$ are test functions and $K^{(i)}$ are suitable infinite-dimensional function spaces containing the exact solution to the problem. Since the centre surface was assumed to be of uniform curvature, $\tilde {\nabla }_{\rm p}^{2} \tilde {C}$ was set to 0 in (2.72b).

The FEniCSx project, which comprises Basix (Scroggs et al. Reference Scroggs, Baratta, Richardson and Wells2022a,Reference Scroggs, Dokken, Richardson and Wellsb), UFL (Alnæs et al. Reference Alnæs, Logg, Ølgaard, Rognes and Wells2014), FFCx and DOLFINx (Baratta et al. Reference Baratta, Dean, Dokken, Habera, Hale, Richardson, Rognes, Scroggs, Sime and Wells2023), was used to find approximate solutions to the variational problem (2.72). All time derivatives were discretised using the backward Euler scheme with a time step of $\Delta \tilde {t} = 0.001$, except where otherwise noted, and the finite element method was used at each time step to solve the time-discrete system of PDEs on a triangulation of the domain with linear Lagrange basis functions used to approximate all variables being solved for. The triangulation was created using Gmsh (Geuzaine & Remacle Reference Geuzaine and Remacle2009), with increased refinement from the centre to the circumference of the circular domain to better resolve large gradients in $\tilde {h}$ and $\tilde { \varGamma }$ near the domain boundary. Table 1 presents the number of cells, $N_{\textrm{cell}}$, used in each simulation. The number of cells used for each simulation was chosen by investigating increasing numbers of cells until the boundary conditions were enforced with sufficient accuracy and the average thickness over time had converged with the number of cells. Newton's method, through an interface from DOLFINx to PETSc's Newton solver (Dalcin et al. Reference Dalcin, Paz, Kler and Cosimo2011), was used to solve the system of nonlinear algebraic equations arising after implicit temporal discretisation and finite element spatial discretisation. Convergence criteria based on the residual were employed with tolerances for the absolute and relative residuals set to $10^{-10}$. The multifrontal massively parallel sparse direct solver (MUMPS) with LU factorisation (Amestoy et al. Reference Amestoy, Duff, L'Excellent and Koster2000) was used to solve the linear system in each Newton iteration. All simulations were run in parallel on 16 processors (Intel Xeon Silver 4216 processor $-$2.10 GHz base frequency) using the message passing interface (Dalcin & Fang Reference Dalcin and Fang2021). For the simulations presented in § 3.4 ($N_{\textrm{cell}} = 195\,400$, $\Delta \tilde {t} = 0.001$, $\tilde {t} \in [0, 0.5]$), the simulation time was approximately 2 h.

The dimensionless numbers that were used for each simulation are provided in table 1. The properties and scales in table 2 of Appendix B, or values similar to those provided, were used to calculate many of the dimensionless numbers in table 1. The characteristic length scale was chosen as 2 mm on the basis that the strength of capillary suction is determined by the meniscus curvature and the characteristic size of the meniscus is given by the capillary length (Gennes et al. Reference Gennes, Brochard-Wyart and Quéré2004; Berg, Adelizzi & Troian Reference Berg, Adelizzi and Troian2005), which is $\sqrt {{\gamma }/(\rho g)} \approx 2$ mm. Furthermore, the velocity scale $U = 4\times 10^{-6}$ m s$^{-1}$ was found by setting the capillary number to $\mathcal {C} = 10^{-7}$ so that capillary effects enter the problem at leading order, i.e. $\epsilon ^3 / \mathcal {C} \sim 1$.

2.8. Validation

Appendix C presents the model derived in this study adapted to free vertical magnetic soap films. Given that the model in Appendix C simplifies to the model in Moulton & Pelesko (Reference Moulton and Pelesko2010) under several assumptions, we can validate our numerical implementation of the model derived in the present study by comparison with the simulation results in Moulton & Pelesko (Reference Moulton and Pelesko2010). With the assumptions used in Moulton & Pelesko (Reference Moulton and Pelesko2010), which are immobile interfaces, uniform surface tension, a linear relationship between the film magnetisation and the magnetic field intensity ($\tilde {M} = \tilde {H}$), no disjoining pressure effects and no evaporation, the model detailed in Appendix C simplifies to

(2.73a)$$\begin{gather} \frac{\partial \tilde{h}}{\partial \tilde{t}} = \frac{1}{3} \tilde{\boldsymbol{\nabla}}_{\rm p} \boldsymbol{\cdot} ( \tilde{h}^{3} ( \tilde{\boldsymbol{\nabla}}_{\rm p} \tilde{p} - \varPsi \tilde{H} \tilde{\boldsymbol{\nabla}}_{\rm p} \tilde{H} - \mathcal{G} \boldsymbol{e}_{x} ) ) , \end{gather}$$

(2.73b)$$\begin{gather}\tilde{p} ={-} \frac{\epsilon^3}{\mathcal{C}} \tilde{\nabla}_{\rm p}^{2} \tilde{h} . \end{gather}$$

This system was solved on the same domain and with the same simulation parameters, magnetic field, initial conditions and boundary conditions as in Moulton & Pelesko (Reference Moulton and Pelesko2010). Figure 2 presents a comparison between our simulation results and the results in figure 7 of Moulton & Pelesko (Reference Moulton and Pelesko2010). Since Moulton & Pelesko (Reference Moulton and Pelesko2010) adopt the thin film approximation in the meniscus region, where it is not valid, $x = 0$ mm and $x = 45$ mm correspond to the very borders of the frame in these simulations. The agreement is almost exact when no magnetic field was applied. There is a slight deviation between the results in figure 2(b–d); however, we believe that this is because the Dirichlet boundary condition deviates from $h = 0.5$ mm at the top boundary ($x = 0$ mm) in the simulations performed by Moulton & Pelesko (Reference Moulton and Pelesko2010) when a magnetic field was applied. As a result, we consider this sufficient validation for the numerical implementation of the derived model. In § 3 the derived model is used to simulate horizontal magnetic soap films.

Figure 2. The grey line represents the initial condition. The black line and blue circles in each figure indicate the thickness after 60 s for the simulations in this study and the simulations performed in Moulton & Pelesko (Reference Moulton and Pelesko2010), respectively. The blue circles were obtained by sampling the data from figure 7 of Moulton & Pelesko (Reference Moulton and Pelesko2010) using WebPlotDigitizer (Rohatgi, Rehberg & Stanojevic Reference Rohatgi, Rehberg and Stanojevic2018). The four cases are (a) no magnetic field applied, (b) 14.5 mm between film and magnet, (c) 8.5 mm between film and magnet, and (d) 2.5 mm between film and magnet.