Thin film modelling of magnetic soap films

Abstract

Magnetic soap films under the forcing of an inhomogeneous magnetic field are governed by a wide range of interconnected physics. The study of magnetic soap films requires the development of comprehensive models to support experimental observations. In this study, the thin film approximation is applied to the Navier–Stokes equations to derive a model for the film thickness of magnetic soap films that incorporates the effects of interfacial mobility, surfactant transport and magnetite nanoparticle (NP) transport. This derived model consists of a coupled system of equations for the film thickness, interfacial velocity, interfacial surfactant concentration and magnetite NP concentration. Simulations are performed for both soap films and magnetic soap films by solving the system of equations using the Galerkin finite element method, and results are compared with experiments. Simulation results highlight that interfacial flows can dominate the rate of film thinning and that accounting for the dependence of the magnetisation on the local magnetite NP concentration can influence the predicted speed of magnetically driven flows. Furthermore, simulation results demonstrate that the model is able to predict marginal regeneration in qualitative agreement with the experiments for soap films; the model also predicts the same flow pattern as seen in the experiments for magnetic soap films. Overall, this study advances the state of soap film and magnetic soap film modelling and will contribute to acquiring control over the drainage and stability of magnetic soap films in the long term.

1. Introduction

Soap films are thin films of liquid stabilised by interfaces laden with surface-active agents (surfactants). In the case of free soap films, a gas phase surrounds the soap film on both sides. This is typical of liquid foams, which contain gas dispersed between (surfactant-stabilised) thin films of liquid (Fameau & Fujii 2020). The stability of soap films and liquid foams is of key relevance to several applications, for example, a high foamability and foam stability is desired for shampoos (Mainkar & Jolly 2001) and fire-fighting foams (Magrabi, Dlugogorski & Jameson 2002), whilst foamability and foam stability is typically less desirable for carbonated drinks (Gonzalez Viejo et al. 2019). Soap films thin over time due to drainage (liquid flows) and evaporation, leading to thinner films, which are less stable and more prone to rupture than thicker films (Poulain, Villermaux & Bourouiba 2018). Adding magnetic nanoparticles (NPs) to soap films, to form magnetic soap films, allows for an applied magnetic field to control the drainage and stability of soap films. Lalli et al. (2023) created free magnetic soap films from aqueous-surfactant solutions containing magnetite (Fe$_3$O$_4$) NPs and reported that an inhomogeneous magnetic field can have either a stabilising or destabilising effect on a magnetic soap film, depending on the concentration of magnetite NPs in the film. It is envisaged that further research devoted to acquiring control over the dynamics and stability of magnetic soap films and further understanding the physics of these fluids could lead to the development of novel technologies. One example is the use of magnetic soap bubbles as stimuli-responsive carriers: magnetic fields could be used to control the position of magnetic soap bubbles and trigger the release of a substance, such as a drug, carried by the bubble (Soetanto & Watarai 2000; Tang et al. 2015).

A wide range of physics is involved in the thinning of both soap films and magnetic soap films. In the following, a non-exhaustive introduction to the key phenomena involved in the thinning of soap films is provided. The thickness of a soap film changes with time due to drainage and evaporation. Drainage flows may be driven by gravity and by gradients in capillary pressure, where the local capillary pressure depends on the local surface tension and interface curvature (Manikantan & Squires 2020). Viscous shear from drainage flows of the core liquid (the liquid in the film between the interfaces) or from a surrounding fluid can drive interfacial flows (Couder, Chomaz & Rabaud 1989). The response of the interface depends on the type and surface coverage of surfactants adsorbed at the interface, as adsorbed surfactants resist shear and dilatational deformation. This viscous response is typically characterised through interfacial shear and dilatational viscosities (Sagis 2011; Li & Manikantan 2021). These interfacial viscosities depend on the type of surfactant and generally increase with the surface coverage of surfactant (Fruhner, Wantke & Lunkenheimer 2000; Denkov et al. 2009; Bhamla et al. 2017), which leads to less mobile interfaces. Immobile interfaces, which feature no interfacial flows, arise in the limit of large interfacial viscosities and lead to slow draining and long-lasting films (Mysels 1968; Bruinsma 1995). In the remainder of the text, films with surfactant-laden interfaces that are not immobile are referred to as partially mobile films, and films with interfaces clean of surface-active agents are referred to as fully mobile films. Interfacial flows advect surfactants along the interface, resulting in Marangoni stresses that can act to immobilise the interface or, alternatively, drive interfacial flows (Manikantan & Squires 2020). Diffusion of surfactant along the interface and local stretching or compression of the interface also affect the interfacial surfactant concentration. It is the change in surface tension upon local stretching or compression that results in the elasticity of soap films and allows soap films to survive against disturbances (Couder et al. 1989). When the surfactants stabilising the interface are insoluble, the change in the surface tension resulting from stretching or compression of the interface, which opposes further stretching or compression, is referred to as the Marangoni elasticity (Manikantan & Squires 2020). When the surfactants are soluble, adsorption/desorption of surfactant to/from the interface additionally affects the local interfacial surfactant concentration. If an interface laden with soluble surfactants is stretched or compressed over a long enough time that adsorption/desorption affects the surface tension, the film elasticity is now referred to as the Gibbs elasticity. Many soluble surfactants also aggregate to form micelles once their concentration in the core liquid increases above the critical micelle concentration. Micelles affect adsorption/desorption kinetics (Lucassen 1975) and can lead to stepwise thinning (Krichevsky & Stavans 1997; Ochoa et al. 2021).

For soap films bounded by a Plateau border, a solid frame or a liquid bath, the soap film can be split into three regions (Barigou & Davidson 1994): (i) a thin film region, in which the film has a thickness between several nanometres and several micrometres (Couder et al. 1989); (ii) a meniscus with thickness of the order of a millimetre that connects the thin film to the border, frame or bath (Gennes, Brochard-Wyart & Quéré 2004); and (iii) a transition region between the thin film and the meniscus. Liquid is sucked from the thin film into the meniscus by capillary suction, which arises since the meniscus is at a lower pressure than the thin film due to its curvature (Overbeek 1960). The flow of liquid into the meniscus causes a region of reduced film thickness, known as a marginal pinch, to form in the transition region (Trégouët & Cantat 2021). Gradients in surface tension can destabilise the marginal pinch (Miguet et al. 2021a), resulting in marginal regeneration, an important drainage mechanism for partially mobile soap films (Mysels 1968; Miguet, Rouyer & Rio 2021b). According to the mechanism reported by Nierstrasz & Frens (1999), in marginal regeneration, elements of film with high surface coverage of surfactant are extracted from the meniscus into the thin film. These elements expand due to Marangoni stresses and become thinner as they expand, which creates extra film area that allows for the continual suction of thicker fluid into the meniscus under the constraint that the surface area of the film remains constant over time.

For soap films with thicknesses of the order of 100 nm or less, interactions between the interfaces start to become significant. Attractive (conjoining) interactions arise, for example, from van der Waals dispersion forces and act to accelerate film thinning (Bergeron 1999; Stubenrauch & Von Klitzing 2003). On the other hand, repulsive (disjoining) interactions, which act to oppose film thinning, may arise due to electrostatic interactions between electric double layers or steric interactions. The strength of these interactions is quantified by the conjoining/disjoining pressure (Bergeron & Radke 1995). Assuming that the film does not rupture first, film thinning may continue until the disjoining pressure becomes large enough to prevent any further thinning. Soap films in this state are described as metastable (Couder et al. 1989; Bergeron 1999) and are either Newton black films or common black films, depending on which interactions contribute most significantly to the disjoining pressure. Common black films usually result from electric double layer forces and have equilibrium thicknesses between approximately 10 and 100 nm (Overbeek 1960; Lyklema & Mysels 1965; Langevin & Sonin 1994). In contrast, Newton black films are thought to be due to entropic confinement forces (Bergeron 1999) and have smaller equilibrium thicknesses of around 5 nm (Overbeek 1960; Couder et al. 1989).

In addition to drainage, evaporation can have a significant effect on film thinning. The rate of evaporation is dependent on the temperature and relative humidity of the surrounding air (Miguet et al. 2020) and on the temperature of the film, which is affected by evaporative cooling (Boulogne, Restagno & Emmanuelle 2022). Furthermore, the presence of surfactant at the interface is expected to decrease the rate of evaporation with the precise change determined by the surface coverage of surfactant (Rubel & Gentry 1984; Marek & Straub 2001). Similarly, the presence of glycerol in the film is expected to slow the rate of evaporation (Boulogne et al. 2022; Pasquet et al. 2022).

When a magnetic field is applied to a soap film containing magnetic NPs, additional mechanisms need to be considered. The core liquid of a magnetic soap film is a colloidal dispersion of magnetic NPs in a non-magnetic base liquid, or a ferrofluid (Philip 2022). Magnetite (Fe$_3$O$_4$) and maghemite ($\gamma$-${\rm Fe}_2{\rm O}_3$) NPs are widely used as the magnetic NPs in ferrofluids. These NPs usually have a single magnetic domain for the particle sizes in a typical ferrofluid (Li et al. 2017; Kole & Khandekar 2021). When a magnetic field is applied to a ferrofluid, the magnetic dipole moments of the magnetic NPs experience a torque driving them to align with the applied field. The ferrofluid becomes magnetised and the magnetisation will increase as the strength of the applied field is increased until saturation, which occurs when the magnetic dipole moments of all NPs is aligned with the applied field (Rosensweig 2013). If the applied magnetic field is inhomogeneous, the magnetic NPs will be driven in the direction of increasing magnetic field intensity (magnetophoresis) and a magnetically induced flow will result in the magnetic soap film due to the interaction between the NPs and the surrounding liquid molecules. The local magnetisation of a ferrofluid, which determines the strength of magnetic forcing, depends on the local concentration of magnetic NPs and varies nonlinearly with the magnetic field intensity. The local concentration of magnetic NPs is affected by magnetophoresis, flow of the base liquid, diffusion, interactions between the magnetic NPs and interactions between the magnetic NPs and other molecules in the film (Pshenichnikov, Elfimova & Ivanov 2011; Nacev et al. 2012; Ayansiji et al. 2020). Furthermore, the nonlinear relationship between the magnetisation and magnetic field intensity depends on the size distribution of magnetic NPs, magnetic dipole–dipole interactions and the thermal energy of the system, which acts to randomly orient the magnetic dipoles. In addition, the magnetisation is affected by demagnetising fields resulting from the magnetisation of the ferrofluid (Richardi & Pileni 2004), which increase in importance with the concentration of magnetic NPs in the ferrofluid (Pshenichnikov & Burkova 2012). An applied magnetic field can also influence the shape of the meniscus of a magnetic soap film, which may impact the strength of capillary suction (Rosensweig et al. 2005; Singh, Singh & Bhaumik 2024).

The speed of drainage flows is closely related to the viscosity of the draining liquid. The viscosity of a ferrofluid is larger than the viscosity of the base liquid due to the dispersed magnetic NPs (Rosensweig 2013), and the ferrofluid viscosity will increase on application of a stationary magnetic field (Odenbach & Thurm 2002), which is the magnetoviscous effect. Magnetic NPs in ferrofluids are typically coated in a surfactant to prevent them from agglomerating (Rosensweig 2013). Interactions between the surfactant coating and surfactants in the film may affect drainage flows. For example, agglomerates containing surfactant molecules and magnetite NPs formed in the magnetic soap films investigated in Lalli et al. (2023), and some of these agglomerates moved at high speeds in the direction of increasing magnetic field intensity, thus affecting drainage flows.

It is evident that the dynamics of soap films and magnetic soap films are controlled by a range of interconnected phenomena, and it is a notorious challenge to decouple and understand the effects of many of these mechanisms from experiments alone (Manikantan & Squires 2020). In this context, modelling can help to reveal the relative importance of different physics and allows for predictions of film behaviour in a range of experimental conditions. In a continuum framework, thin films, bubbles and foams are frequently modelled using interface-tracking or interface-capturing methods (Tezduyar 2006), which require each interface to be followed whilst imposing mass conservation and a momentum balance at each point on the interface. These methods are computationally expensive and require sophisticated techniques to avoid spurious oscillations around each simulated interface (Denner et al. 2017). Thin films can also be modelled using the thin film approximation (lubrication theory), which involves capitalising on the thickness of the film being much smaller than the lateral extent of the film to derive an evolution equation for the film thickness that incorporates essential physics (O'Brien & Schwartz 2002; Craster & Matar 2009). Thin film models, i.e. models derived using the thin film approximation, are computationally less demanding and involve simpler numerics than interface-tracking and interface-capturing methods; therefore, we opted to derive a thin film model to describe the dynamics of soap films and magnetic soap films.

Several thin film models have been developed for the thickness of free soap films. Barigou & Davidson (1994) modelled a circular soap film bounded by copper tubes and accounted for the difference in mobility caused by varying degrees of surfactant adsorption between the thin film, the transition region and the meniscus. However, the thin film was treated simply by assuming fully mobile interfaces with a uniform velocity over the thickness. Aradian, Raphael & De Gennes (2001), too, considered distinct behaviour in each of the three regions, but assumed rather simplified behaviour in the thin film by treating the interfaces as immobile with a uniform surface tension. A thin film model relevant to soap films close to rupture was derived by De Wit, Gallez & Christov (1994). The resulting model consists of a system of three nonlinear equations governing the film thickness, interfacial surfactant concentration, assuming insoluble surfactants, and tangential interfacial velocity. This model has less relevance to soap films with thicknesses greater than approximately 100 nm. Schwartz & Roy (1999) assumed that the flow velocity can be divided into shear and extensional contributions and that the two flows contribute equally to the pressure gradient to derive a model for the thinning of vertical soap films. This model is expected to be valid over a greater range of thicknesses than the model of De Wit et al. (1994) and additionally includes disjoining pressure effects. Nierstrasz & Frens (1999) incorporated interfacial shear and dilatational viscosities, in accordance with the Boussinesq–Scriven surface stress model, into their thin film model for vertical soap films. Soluble surfactants were also considered. In order to model soluble surfactants, equations governing the surfactant concentration in the core liquid and on the interface are required in addition to relations for the rate of surfactant adsorption and desorption (Manikantan & Squires 2020). For thin film models, the equation governing the concentration of surfactant in the core liquid must be integrated over the film thickness to remove the dependence on the depth coordinate, which requires the variation of surfactant concentration over the film depth to be assumed.

The aforementioned studies derive two-dimensional thin film models (film thickness as a function of one spatial dimension); therefore, marginal regeneration can not be simulated, as there can not be an exchange between thinner fluid and thicker fluid at the boundary of a thin film represented by only a single point. Joye, Hirasaki & Miller (1996) and Naire, Braun & Snow (2004) derived three-dimensional thin film models for soap films (film thickness as a function of two spatial dimensions), similar to the model of Nierstrasz & Frens (1999) but under the assumption of insoluble surfactants. Both studies were able to simulate an instability driven by gradients in surface tension that is expected to lead to marginal regeneration. Despite these studies, marginal regeneration in a soap film has yet to be simulated using a thin film model. In addition, it is unclear why several studies have imposed a no-flux boundary condition for soap films supported by a solid wire or frame (Schwartz & Roy 1999; Braun, Snow & Naire 2002) since the thin film approximation is not valid in the meniscus and a net flow of liquid between the thin film and the meniscus is expected. Consequently, a key focus of the present study is to derive a three-dimensional thin film model and solve the resulting system with a non-zero flux boundary condition, with the aim of simulating marginal regeneration.

Several experimental studies have focused on investigating whether magnetic fields can provide an element of control over free magnetic soap films, bubbles and foams: Sudo, Hashimoto & Katagiri (1990) discovered that a magnetic field can be used to stabilise magnetic liquid foams by reducing the rate of foam coarsening. Furthermore, Hutzler et al. (2002) were able to alter the structure of magnetic soap foams in cylindrical tubes using inhomogeneous magnetic fields, and Drenckhan et al. (2003) were able to control the size of the bubbles forming these foams by adjusting the strength of an applied inhomogeneous magnetic field. Homogeneous magnetic fields can also affect film dynamics due to the magnetisation induced by the applied field. Elias et al. (2005) observed that the speed of drainage flows in vertical magnetic soap films was affected by the application of a homogeneous magnetic field, and Mawet et al. (2023) found that the shape of hemispherical magnetic soap bubbles sitting on a solid substrate was altered by an applied homogeneous magnetic field. Lalli et al. (2023) chose to adopt an inhomogeneous magnetic field in their experimental investigation of magnetic soap films, as inhomogeneous magnetic fields provide a more direct control over drainage flows through magnetophoresis. In a different configuration, inhomogeneous magnetic fields are also advantageous in that they can control the position of freely floating magnetic soap bubbles (Rodrigues et al. 2011). The focus of these studies investigating free magnetic soap films, bubbles and foams was predominately experimental, such that only Elias et al. (2005) developed a model for film thinning, which is a simple two-dimensional thin film model accounting for gravity, capillary and demagnetising field effects. Moulton & Pelesko (2010) and Back & Beckham (2012) investigated free magnetic soap films formed on a vertical frame in inhomogeneous magnetic fields. Moulton & Pelesko (2010) derived a two-dimensional thin film model for the film thickness, which was then extended to three dimensions by Back & Beckham (2012) to account for the two-dimensional nature of the applied magnetic field in the plane of the film. Although both studies found some qualitative agreement between simulated films and their experimental results, their models disregard essential physics. For example, their models make the restrictive assumption of immobile interfaces, whilst experimental images in Moulton & Pelesko (2010) and Back & Beckham (2012) suggest that they were investigating partially mobile magnetic soap films (Mysels 1968; Bruinsma 1995). Furthermore, it was assumed that the magnetisation of the core liquid was linearly related to the applied magnetic field intensity through a constant and uniform magnetic susceptibility, which implicitly assumes that the concentration of magnetic NPs is constant and uniform and that a linear relation exists between the magnetisation and applied field intensity. Both studies adopt the no-flux boundary condition at all boundaries, which, as discussed earlier, is inconsistent with the thin film approximation only being valid within the thin film region. The modelling of magnetic soap films has received no further attention, other than the extension of the model in Moulton & Pelesko (2010) to include disjoining pressure effects (Moulton & Lega 2013), such that the state of magnetic soap film modelling is rather limited.

The aim of this study is to advance the state of magnetic soap film modelling, using the experimental configuration and results in Lalli et al. (2023) as a reference. The objectives of this work are then to (i) derive a model for the two-dimensional thickness field of free magnetic soap films under the action of inhomogeneous magnetic fields that includes interfacial mobility, surfactant transport, magnetite NP transport, the dependence of the film magnetisation on the magnetite NP concentration and magnetic field intensity, conjoining/disjoining pressure effects and evaporation; (ii) use the derived model to study the effect of interfacial flows on the rate of film thinning and to understand how accounting for the dependence of the magnetisation on the magnetite NP concentration affects film thickness predictions; (iii) investigate the capability of the model to predict marginal regeneration; and (iv) compare simulated films with the experimental results in Lalli et al. (2023).

The remainder of this paper is organised as follows. Section 2 provides the derivation of the model and the method for solving the system of equations governing the thickness evolution. Section 3.1 examines the contribution of interfacial flows to the rate of film thinning, and § 3.2 investigates how accounting for magnetite NP transport and the dependence of the magnetisation on the magnetite NP concentration affects film thickness predictions. Subsequently, §§ 3.3 and 3.4 present simulations that were performed to allow for a comparison with experiments, and § 4 evaluates the model by comparing these simulations with experiments. Section 4 also provides directions for future work. Concluding remarks close the paper.

2. Method

This section presents the derivation of a three-dimensional thin film model for the free soap films and magnetic soap films investigated in Lalli et al. (2023). The method for solving the resulting system of equations is then introduced and validated.

2.1. Problem description and key assumptions

The films investigated in Lalli et al. (2023) 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. (2023) 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 1967; Sagis 2011). 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 2009), 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 2013), the magnetic part of the stress tensor given in Cowley & Rosensweig (1967) 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 1987). The time scale for magnetic relaxation (Rosensweig 2013) 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. 2023). 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 1964; Cowley & Rosensweig 1967)

(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 2007). The jump mass balance, which is a statement of conservation of mass at each point on the dividing surface, then requires (Burelbach, Bankoff & Davis 1988; Slattery et al. 2007)

(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. 2007),

(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. 1988; Bai & Gosman 1996)

(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. 2007)

(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 2013),

(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. 2007),

(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. 2007)

(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 2004; Münch, Wagner & Witelski 2005). 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 2002).

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. 2011; Nacev et al. 2012)

(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 1983)

(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 1983),

(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 2003). 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 2015; Kim et al. 2016). For particles with size of the order of the film thickness, Prasad & Weeks (2009) 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 2009; Vivek & Weeks 2015); 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 2013)

(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 1999; Huybrechts, Villaret & Hervouet 2010) 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 1987). 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 2020),

(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 1990)

(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 2020):

(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 1996; Pereira et al. 2007). Each position on the dividing surface can be specified by two parameters known as the surface coordinates (Slattery et al. 2007), 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. (1996) and Cermelli, Fried & Gurtin (2005) 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.1. Evaporation rate

In the simplest case, soap films can be simulated assuming a constant and uniform rate of evaporation (Huang et al. 2020; Ishida et al. 2020; Pasquet et al. 2022). A more accurate treatment of evaporation involves solving the energy balance equation in the core liquid and the jump energy balance equation at the dividing surface and introducing a constitutive relation between the mass flux due to evaporation and the interface temperature (Burelbach et al. 1988; Sultan, Boudaoud & Amar 2005). To isolate the effects of drainage phenomena, all simulations were conducted without evaporation, i.e. $\tilde {J}_{{\rm e}} = 0$. We leave it as an area of future work to couple the system of equations introduced in the present contribution with the energy balance equations.

2.6.2. Disjoining pressure

The Derjaguin–Landau–Verwey–Overbeek (DLVO) theory was used for the disjoining pressure in this study (see Ninham 1999 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 1999; Matsubara et al. 2021; Agmo Hernández 2023)

(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. 2021). 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 1999).

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 2019). A fourth-order PDE of this type can be solved with mixed finite element methods (Wells, Kuhl & Garikipati 2006; Liu et al. 2019; Keita, Beljadid & Bourgault 2021), 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 1960; Mysels 1968), 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. 2005; Singh et al. 2024); 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 (1999) 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 2017):

(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. 2022a,b), UFL (Alnæs et al. 2014), FFCx and DOLFINx (Baratta et al. 2023), 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 2009), 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. 2011), 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. 2000) 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 2021). 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. 2004; Berg, Adelizzi & Troian 2005), 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 (2010) 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 (2010). With the assumptions used in Moulton & Pelesko (2010), 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 (2010). Figure 2 presents a comparison between our simulation results and the results in figure 7 of Moulton & Pelesko (2010). Since Moulton & Pelesko (2010) 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 (2010) 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 (2010), respectively. The blue circles were obtained by sampling the data from figure 7 of Moulton & Pelesko (2010) using WebPlotDigitizer (Rohatgi, Rehberg & Stanojevic 2018). 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.

3. Results

In § 3.1, simulation results are used to investigate the significance of interfacial flows on the rate of film thinning, and in § 3.2, simulations results are used to study how accounting for magnetite NP transport and the dependence of the magnetisation on the magnetite NP concentration affects film thickness predictions. The simulations presented in §§ 3.1 and 3.2 were performed on a dimensionless domain defined by $\tilde { \varOmega } = \{(\tilde {x}, \tilde {y}) \mid (\tilde {x} - 0.5)^2 + (\tilde {y} - 0.5)^2 \leq 0.25 \}$; the simulation parameters were chosen to investigate the phenomena of interest and are not necessarily physically realistic. In §§ 3.3 and 3.4, more realistic simulation parameters were chosen to align with the experiments in Lalli et al. (2023). The length scale $L$ was held fixed at the capillary length, and the simulations were conducted on a larger dimensionless domain defined by $\tilde { \varOmega } = \{(\tilde {x}, \tilde {y}) \mid (\tilde {x} - 9.125)^2 + (\tilde {y} - 9.125)^2 \leq 18.25 \}$. Table 1 provides the dimensionless numbers that were used in each simulation.

3.1. Effects of surface mobility and surface flows – soap films

The assumption of immobile interfaces is commonly made in thin film models for soap films and magnetic soap films. In this section we compare immobile soap films ($\textit {Bq}_{\textrm{sh}} \rightarrow \infty$ and $\textit {Bq}_{\rm d} \rightarrow \infty$, or sufficient Marangoni stresses) with partially mobile soap films to investigate the importance of including surface mobility and surface flows in thin film models, where the surface mobility can be thought of as the inverse of the Boussinesq numbers (Stone et al. 2002; Dame et al. 2005), although the surface mobility also depends on Marangoni effects (Stebe, Lin & Maldarelli 1991). Note that the fully mobile limit occurs when the interfaces are free of surface-active agents ($\textit {Bq}_{\rm sh} = \textit {Bq}_{\rm d} = 0$, $\gamma = \gamma _{{\rm c}}$), which is not investigated here as surfactants are a prerequisite for films to survive for any reasonable amount of time. Gravitational effects were set to 0 ($\mathcal {G} = 0$) to isolate the effects of interfacial flows that arise in partially mobile films. An immobile soap film was simulated by assuming that the interface was saturated with surfactant, i.e. $\gamma = \gamma _{\textrm{sat}}$ and $\tilde { \gamma } = 0$, which was imposed in the simulations by setting the Marangoni number to $\mathcal {M} = 0$. Marangoni stresses arising due to interfacial flows driven by viscous shear from the core liquid are expected to have an immobilising effect on the interface. Consequently, the Marangoni number was varied alone for the simulated partially mobile soap films to investigate interfacial flows of varying magnitudes.

In all simulations, the thickness field remained axisymmetric over time. Figure 3(a) presents the thickness profiles of an immobile soap film over time. The immobile film thins due to the capillary suction boundary condition with all drainage within the film governed by spatial variations in the capillary pressure. The thinning slows with time due to the $\tilde {h}^3$ dependence in the thickness evolution equation (2.42a) and the decrease in capillary suction force at the boundary with time (2.66). Furthermore, as time progresses, a region of increased film thickness forms in the centre surrounded by a thinner ring. This behaviour has been reported for circular thin films and occurs when the rate of thinning due to capillary suction exceeds the curvature-related thinning at the centre of the film (Joye, Hirasaki & Miller 1992; Pugh 2016).

Figure 3. (a) Film thickness profiles for $\tilde {t} \in [0, 1]$ along a line passing through the centre of the film at $\tilde {y} = 0.5$ for an immobile soap film and a partially mobile soap film with $\mathcal {M} = 100$. The time between the thickness profiles is $\Delta \tilde {t} = 0.04$, and time increases in the vertically downwards direction, indicated by darker colours. (b) The average film thickness over the thin film as a function of time for an immobile film and three partially mobile films.

The thinning progresses considerably faster for the partially mobile soap film shown in figure 3(a) compared with the immobile film. This occurs because a surface velocity field develops that is radially outward from the centre due to viscous shear forcing from the core liquid that is also moving radially outwards from the centre. Surfactants adsorbed at the interface are advected radially outwards from the film centre, leading to Marangoni stresses that act to decelerate the surface velocity.

Faster surface flows arise with smaller Marangoni numbers, causing the average thickness over the film to fall faster with time, as shown in figure 3(b). The film thickness reaches an equilibrium when the disjoining pressure is sufficient to prevent any further drainage or evaporation. The precise equilibrium thickness obtained in this set of simulations is dependent upon the disjoining pressure parameters used in (2.63) and the strength of the capillary suction boundary condition. For the partially mobile film with $\mathcal {M} = 10$, we were able to simulate an equilibrium film by reducing the time step and slowly ramping down the Neumann slope condition from 30$^\circ$ to 0$^\circ$ once disjoining pressure effects became significant to avoid numerical instabilities caused by large gradients in the disjoining pressure. The simulated equilibrium film thickness ($2h$) was 45.6 nm ($\tilde {h} = 2.28\times 10^{-3}$), which is a typical thickness of common black films (Bergeron 1999).

Figure 3(b) presents how the rate of film thinning approaches the rate of thinning for the immobile soap film as the Marangoni number is increased, highlighting how Marangoni stresses can immobilise the interface. Overall, the results presented in this section highlight the importance of including interfacial flows in models for the thinning of soap films, as interfacial flows can significantly affect the rate of film thinning.

3.2. Magnetite NP transport – magnetic soap films

The drainage flows in magnetic soap films induced by an inhomogeneous magnetic field are caused by the transport of magnetic NPs within the film; this transport can affect the local magnetic NP concentration and the film magnetisation. Despite this, previous models for the thinning of magnetic soap films under the forcing of inhomogeneous magnetic fields calculate the local film magnetisation by implicitly assuming that the concentration of magnetite NPs remains uniform throughout the film (Moulton & Pelesko 2010; Back & Beckham 2012). It is therefore of interest to understand how accounting for variations in the magnetite NP concentration and the dependence of the film magnetisation on the magnetite NP concentration affects film thickness predictions. Consequently, simulations were performed with magnetite NP transport by solving for the local molar concentration of NPs using (2.47) and without magnetite NP transport by imposing a constant and uniform concentration of magnetite NPs. The present results evaluate the possible effects and importance of accounting for magnetic NP transport when simulating magnetic soap films in inhomogeneous magnetic fields. The magnetic field used in the simulations is presented in figure 12(a).

With the capillary suction boundary condition (2.66), the thickness initially increases from right to left due to magnetically induced flows in the direction of increasing magnetic field intensity, as shown in figure 4(d). The thickness then falls to a minimum on the left side of the domain, as illustrated by figure 4(a,d), due to the magnetic drainage flux out of the thin film, which is given by $\frac {1}{3} \tilde {h}^3 \varPsi \tilde {M}(\tilde {c}, \tilde {H}) \tilde {\boldsymbol {\nabla }}_{\rm p} \tilde {H} \boldsymbol {\cdot } \boldsymbol {n}_{\rm d}$. The thickness on the left side of the domain falls faster when accounting for magnetite NP transport since the concentration of magnetite NPs, and therefore, the magnetisation, is greater on the left side of the film with magnetite NP transport compared with the case of uniform concentration. In contrast, the magnetic soap films in the experiments generally increased in thickness in the direction of increasing magnetic field intensity (Lalli et al. 2023). As a result, we also investigated boundary condition (2.67), which limits the drainage flux of core liquid out of the thin film. With this boundary condition, we see a continual increase in thickness from right to left in figure 4(b,e). The difference between the thickness profiles for uniform $\tilde {c}$ and with the transport of $\tilde {c}$ is now smaller. Nevertheless, due to the build up of concentration of magnetite NPs on the left side of the film and the depletion of magnetite NPs on the right side of the film, as illustrated in figure 4(c,f), there is a larger thickness on both the left and right sides of the film when considering magnetite NP transport. This is because magnetically driven flows are stronger on the left side of the film and weaker on the right side of the film when accounting for magnetite NP transport.

Figure 4. The top row presents $\tilde {h}$ and $\tilde {c}$ at $\tilde {t} = 0.10$. Panels (a,b) present the thickness field with magnetite NP transport for boundary conditions (2.66) and (2.67), respectively. Panel (c) presents the concentration field for magnetite NPs for boundary condition (2.67). The profiles in (d), (e) and (f) are plots over the dashed centrelines marked in (a), (b) and (c) over time for $\tilde {t} \in [0.0, 0.2]$. The time between profiles is $\Delta \tilde {t} = 0.02$, and the progression of time is indicated by arrows and progressively darker colours. In (d) and (e), the dashed profiles are for a uniform concentration of magnetite NPs and the solid profiles are for magnetite NP transport governed by (2.47).

With boundary condition (2.67), the thickness profiles are more realistic; however, this boundary condition results in $\tilde {\boldsymbol {\nabla }}_{\rm p} \tilde {p} _{{\rm c}} \boldsymbol {\cdot } \boldsymbol {n}_{\rm d}$ varying around the boundary. Although it seems reasonable that there would be a limit on the drainage flux of core liquid out of the thin film, it is unclear, at present, whether there is any physical justification for $\tilde {\boldsymbol {\nabla }}_{\rm p} \tilde {p} _{{\rm c}} \boldsymbol {\cdot } \boldsymbol {n}_{\rm d}$ varying around the boundary of the thin film, but it could possibly be attributed to the effect of the magnetic field on the meniscus shape (Rosensweig et al. 2005; Singh et al. 2024). Overall, accounting for the transport of magnetite NPs in thin film models can significantly affect film thickness computations since the magnetisation depends on the local NP concentration. The size of this effect is strongly dependent on the dimensionless numbers $\alpha$, $\varphi _{{\rm c}}$ and $Pe_{c}$, which control the strength of magnetophoresis and diffusion relative to advective transport, and the choice of the boundary condition for $\tilde {c}$. Further research is needed to understand the behaviour of magnetite NPs at the thin film boundary, $\partial \varOmega$, to justify the boundary condition used for $\tilde {c}$ (2.69) or to suggest an alternative boundary condition that may be more realistic.

3.3. Simulating experimental soap films

There are a number of unknowns from the experiments, such as the surface shear and dilatational viscosities, the initial surface excess concentration of surfactant and the initial film thickness (the films were initially too thick to exhibit interference colours (Lalli & Giusti 2023)) such that it is only possible for a qualitative comparison between the simulations and the experiments.

The length scale $L$ was held fixed at the capillary length and the size of the dimensionless simulation domain was increased to $\tilde { \varOmega } = \{(\tilde {x}, \tilde {y}) \mid (\tilde {x} - 9.125)^2 + (\tilde {y} - 9.125)^2 \leq 18.25 \}$ to align with the experiments. With a domain of this size, capillary pressure effects are no longer significant over the entire domain, with negligible capillary pressure gradients in the centre of the domain, in contrast to the soap films simulated in § 3.1. Based on the method used to create the soap films in the experiments, we hypothesise that the pressure was slightly larger under the soap films studied. As a result, we imposed an axisymmetric centre surface, $C$, of uniform curvature in the simulations, created by slicing a spherical bubble of radius 100 mm at a level such that $\tilde {C}$ was 168.0 in the centre of the film and 0.0 at the thin film boundary.

Figure 5 presents the dimensionless thickness field over a segment of the simulation time for two initial surface concentrations of surfactant (top row is $\tilde { \varGamma }_{\rm i} = 0.05$ and bottom row is $\tilde { \varGamma }_{\rm i} = 0.20$). With an initially uniform thickness field, a minimum in the thickness field arises in the centre of the film for both initial concentrations of surfactant, primarily due to gravitational drainage and shear-driven surface flows. The surface flows advect surfactant outwards from the film centre, resulting in larger Marangoni stresses with increasing distance from the film centre, as illustrated in figure 6(a) for $\tilde {t} = 0.01$. The film thickness and surfactant distribution is axisymmetric at this time. Furthermore, the Marangoni stresses are similar with $\tilde { \varGamma }_{\rm i} = 0.05$ and $\tilde { \varGamma }_{\rm i} = 0.20$ since the Marangoni number, which features in (2.42c), is the same in each instance. After some thinning has taken place, the axisymmetry in the simulations is broken by the large Marangoni stresses at the boundary of the thin film combined with the suction boundary condition (2.66). Figure 5 presents how film that is thinner than the surrounding film is then frequently extracted from the boundary of the film before being driven towards the film centre. These thin patches of fluid initially have a higher surfactant concentration than the surrounding film, and their motion towards the film centre is driven by Marangoni stresses and the difference in gravitational forcing on thin patches of fluid compared with surrounding thicker fluid (Huang et al. 2020). This simulated behaviour appears to be in agreement with the mechanism for marginal regeneration outlined by Nierstrasz & Frens (1999).

Figure 5. Simulated thickness fields for $\tilde { \varGamma }_{\rm i} = 0.05$ (a) and $\tilde { \varGamma }_{\rm i} = 0.20$ (b) for $\tilde {t} \in [0.03, 0.165]$. The time between each thickness field is $\Delta \tilde {t} = 0.015$.

Figure 6. (a) Strength of Marangoni stresses relative to viscous shear at $\tilde {t} = 0.01$ along a line passing through the centre of the film at $\tilde {y} = 9.125$ for $\tilde {\varGamma }_{\rm i} = 0.05$ and $\tilde {\varGamma }_{\rm i} = 0.20$. The grey line represents the initial condition. (b) The average thickness over the entire film (solid line) and over a central circular section of film defined by $\tilde {\varOmega }_{{\rm c}} = \{(\tilde {x}, \tilde {y}) \mid (\tilde {x} - 9.125)^2 + (\tilde {y} - 9.125)^2 \leq 9.125 \}$ (dashed line) as a function of time for $\tilde {\varGamma }_{\rm i} = 0.05$ and $\tilde {\varGamma }_{\rm i} = 0.20$.

Figure 6(b) presents the variation of the average thickness over a central section of film and over the entire film as a function of time. There is initially a sharp change in the thickness averaged over the centre of the film for both initial surfactant concentrations. This is caused by large surface velocities arising in the initial absence of surface tension gradients, which are driven by viscous shear from the core liquid. Marangoni stresses acting to immobilise the interface arise quicker with the larger initial concentration of surfactant, $\tilde { \varGamma }_{\rm i} = 0.20$. Consequently, the film with a smaller initial concentration of surfactant drains considerably faster over time.

3.4. Simulating experimental magnetic soap films

Two magnetic soap films with different concentrations of magnetite NPs ($\tilde {c}_{\rm i} = 0.5$ and $\tilde {c}_{\rm i} = 1.0$) were simulated with magnetic field intensity given by figure 12(c). As $Pe_{c} \rightarrow \infty$ in the magnetite NP transport equation (2.47), the transport of $\tilde {c}$ is dominated by advection. In this limit, the solution for $\tilde {c}$ is uniform for a uniform initial concentration of magnetite NPs. Preliminary simulations showed that the solution for $\tilde {c}$ may faithfully be taken as uniform for $Pe_c = 2000$ and the other dimensionless number given in table 1. Consequently, a uniform concentration of magnetite NPs was adopted for the simulations performed in this section for computational efficiency. A centre surface of uniform curvature was imposed with smaller height than in § 3.3 to satisfy the assumption $\partial {H}/\partial {z} = 0$ over the film depth. This surface was obtained by slicing a bubble of radius 1 m such that $\tilde {C}$ was 16.7 in the centre of the film and 0.0 at the thin film boundary.

Figure 7 presents the variation of the average thickness over the left side of the thin film, right side of the thin film and entire thin film as a function of time for the two magnetic soap films studied. Magnetically induced flows in the direction of increasing magnetic field intensity drove the average thickness to be larger on the left side of the film compared with the right side. The difference in average thickness between the left and right sides of the film is greater for the more concentrated magnetic soap film since a larger concentration of magnetite NPs leads to stronger magnetically driven flows. The average thickness over the entire thin film falls faster with $\tilde {c}_{\rm i} = 1.0$ compared with $\tilde {c}_{\rm i} = 0.5$, as the drainage flux out of the thin film due to magnetic forcing, $\frac {1}{3} \tilde {h}^3 \varPsi \tilde {M}(\tilde {c}, \tilde {H}) \tilde {\boldsymbol {\nabla }}_{\rm p} \tilde {H} \boldsymbol {\cdot } \boldsymbol {n}_{\rm d}$, is larger for $\tilde {c}_{\rm i} = 1.0$.

Figure 7. Average thickness over the left-half ($\tilde {\varOmega }_{\mathrm {l}}$), right-half ($\tilde {\varOmega }_{\mathrm {r}}$) and entire thin film ($\tilde {\varOmega }$) as a function of time for $\tilde {c}_{\rm i} = 0.5$ and $\tilde {c}_{\rm i} = 1.0$.

Thickness fields over a subset of the simulation time are presented in figure 8 for the two magnetic soap films studied. Figure 8 shows that film that is thinner than the surrounding fluid is sporadically extracted from the film boundary before being driven away from the film boundary. Thin patches of fluid extracted from the left side of the film move approximately in the $\boldsymbol {e}_{x}$ direction, while thin patches extracted from the right side of the film are deflected away from moving radially towards the centre of the film.

Figure 8. Simulated thickness fields for $\tilde {c}_{\rm i} = 0.5$ (a) and $\tilde {c}_{\rm i} = 1.0$ (b) for $\tilde {t} \in [0.20, 0.38]$. The time between each thickness field is $\Delta \tilde {t} = 0.02$.

The simulated films feature fast flows in the direction of increasing magnetic field intensity that drive fast flows around the circumference of the film in the reverse direction. This behaviour is evident from the surface velocity fields in figure 9(a,c). The flows around the boundary eventually separate from the boundary and orient towards the interior of the film. The separation can be seen in figure 8 for $\tilde {t} = 0.20$, with thinner fluid being driven from the boundary at the point of separation according to the surface flows in figure 9(a,c). The point of separation occurs further to the right for the magnetic soap film with a higher concentration of magnetite NPs, as the magnetically driven surface velocities are larger for this film. In figure 9(b,d) the surface velocities are smaller due to Marangoni stresses, and the separation of the flows from the boundary occur earlier (further to the left) for each of the two magnetic soap films. Again, the separation occurs earlier for $\tilde {c}_{\rm i} = 0.5$ compared with $\tilde {c}_{\rm i} = 1.0$.

Figure 9. A subset of surface velocity vectors for $\tilde {c}_{\rm i} = 0.5$ at (a) $\tilde {t} = 0.20$ and (b) $\tilde {t} = 0.38$ and for $\tilde {c}_{\rm i} = 1.0$ at (c) $\tilde {t} = 0.20$ and (d) $\tilde {t} = 0.38$, where $\tilde {V}_{\rm s}$ is the magnitude of $\tilde {\boldsymbol {V}}_{\rm s}$.

4. Comparison with experiments

This section evaluates the performance of the model by comparing the results in §§ 3.3 and 3.4 with experiments. Recommendations for future work are also discussed.

4.1. Soap films

To the best of our knowledge, marginal regeneration has not been simulated before using a thin film model. The computations performed using the thin film model derived in the present study were in qualitative agreement with the drainage due to marginal regeneration observed in the experiments, namely, thin patches of fluid were extracted from the boundary of the film before moving towards the centre of the film whilst thicker fluid was sucked into the meniscus. The trajectories of the thin patches in figure 5 are in qualitative agreement with the trajectories illustrated in figure 9 of Lalli et al. (2023).

There are some key differences between the marginal regeneration behaviour observed in the simulations (see figure 5) and in the experimentally investigated soap films (see figure 10a–c), which are now discussed. In the experiments, marginal regeneration started almost instantly after film formation, whilst the breaking of axisymmetry occurred later after around 20 s in the simulations. The mechanism that leads to pinch destabilisation and marginal regeneration, which is yet to be fully understood (Trégouët & Cantat 2021), is not incorporated into the model since the pinch lies outside of the simulated region; therefore, it may be more realistic to trigger marginal regeneration in the simulations by perturbing the film thickness or surface surfactant concentration. As evident from figures 5 and 10(a–c), the frequency of thin patches was much larger in the experiments compared with the simulations, suggesting that marginal regeneration contributed more significantly to the rate of thinning in the experiments. The thin patches were also typically larger in the simulations compared with the experiments. Since the frequency and size of thin patches is expected to be closely tied to the pinch destabilisation mechanism, a more accurate treatment of this destabilisation mechanism may lead to a better agreement between the simulations and the experiments. Furthermore, the thin patches of fluid in the experiments had well-defined, although deformable, boundaries, while the simulated thin patches had more diffusive boundaries. We hypothesise that the behaviour observed in the experiments is governed by the surfactant monolayers covering the thin patches of fluid, preventing surrounding thicker fluid from flowing into the thin patches and equilibrating their thickness with the surrounding fluid. It would be interesting to conduct further research to explore how thin patches of fluid interact with surrounding fluid, which could inform how the model could be improved to more accurately simulate marginal regeneration. The diffusive nature of the thin patches of fluid is also related to the choice of simulation parameters, such as the Péclet number for surface surfactant transport. Additionally, it may be important to account for the variation of the interfacial viscosities with the surface concentration of surfactant. Spatial gradients in the surface shear and dilatational viscosities would then enter the equation governing the surface velocity (2.42c) from the Boussinesq–Scriven surface stress model (2.15).

Figure 10. The thinning of a soap film (a–c) and magnetic soap film in an inhomogeneous magnetic field (d–f) (Lalli et al. 2023), where $t$ is the time from film formation. The relationship between film thickness and interference colours was computed by applying an interference relation derived for monochromatic waves at a number of discrete wavelengths (Lalli & Giusti 2023).

A black film formed in the centre of experimentally investigated soap films, as seen in figure 10(a), within around 10 s from film formation. The black centre then grew with time, as evident from figure 10(b,c), due to drainage and evaporation. Although a minimum in thickness was observed in the centre of the simulated soap films, black film did not form. This may be attributed to the aforementioned differences in drainage due to marginal regeneration and since evaporation was switched off in all simulations.

As many surfactants are soluble, future work should focus on extending the model from insoluble surfactants to soluble surfactants. This will allow the effects of surfactant adsorption and desorption on film thinning to be explored.

4.2. Magnetic soap films

The simulated magnetic soap films shown in figure 8 feature an increase in thickness in the direction of increasing magnetic field intensity (from right to left), in alignment with the magnetic soap film shown in figure 10(d–f). Furthermore, the simulations predicted faster magnetically driven drainage flows with a larger initial concentration of magnetite NPs in the magnetic soap film, as was also reported by Lalli et al. (2023).

The arrows in figure 10(d) illustrate the direction of the fastest drainage flows: there are flows towards the magnet in the centre of the film with reverse flows around the film boundary that then separate from the film boundary and move towards the interior of the film. Our model predictions in figures 8 and 9 predict this same flow pattern.

The trajectories of the thin patches of fluid in figure 8 are in good agreement with the trajectories in figure 9 of Lalli et al. (2023); however, in similarity with soap films, the size and frequency of thin patches differs significantly between simulation and experiment. For example, many thin patches of fluid are present in the reverse flows around the boundary in figure 10(d–f), whilst only a few thin patches of fluid were observed in these flows in the simulated films in figure 8. As previously mentioned, a more accurate treatment of the pinch destabilisation mechanism could lead to a closer agreement between simulation and experiment. Furthermore, the frequency of thin patches of fluid depends on the proposed capillary suction boundary condition (2.66) and the suction strength, $\mathcal {P}$. This boundary condition allows the model to reproduce behaviour in qualitative agreement with marginal regeneration. Nevertheless, further research is needed to validate this boundary condition for film borders subject to capillary suction. In addition, further research is needed to justify a suitable value for the Neumann boundary condition for the film thickness, based on the behaviour in the transition region and the meniscus.

The derived model assumes a dilute concentration of magnetite NPs in the core liquid. The validity of the model could be extended to higher concentrations of magnetite NPs by considering interactions between the NPs, magnetoviscous effects and demagnetising field effects. Furthermore, the Langevin relation (2.46) describing the nonlinear magnetisation of the ferrofluid assumes that the magnetite NPs in the ferrofluid are monodisperse. Since most ferrofluids are polydisperse, it may be valuable in future work to investigate magnetisation relations derived for polydisperse NPs (Ivanov et al. 2007).

With the simulation parameters used in this study, convergence is achieved with a relatively large time step, $\Delta \tilde {t} = 0.001$ (Liu et al. 2019). It is anticipated that some form of stabilisation, such as streamline-upwind Petrov–Galerkin stabilisation or a subgrid-scale method (Codina 1998; John & Knobloch 2007; Liu et al. 2019), will be needed when forcing terms are particularly large, such as large values of $\varPsi$ or $\mathcal {G}$, or when the Péclet numbers for magnetite NP and surfactant transport are particularly large. Finally, future work should focus on defining the range of variables and dimensionless numbers over which the pressure scale $P = {\eta U}/(\epsilon \mathcal {H})$ remains appropriate. For films close to the fully mobile limit, the pressure scale suitable for extensional flows, $P = {\eta U}/{L}$, is expected to be more relevant (Kargupta et al. 2004; Münch et al. 2005), whilst for films close to rupture, inertial terms become significant such that the pressure scale $P = \rho U^2$ would be more relevant (Erneux & Davis 1993; De Wit et al. 1994).

5. Concluding remarks

A system of equations governing the thickness field of free magnetic soap films has been derived using the thin film approximation. This system of equations incorporates several physical phenomena that were not accounted for in previous magnetic soap film models, such as interfacial flows and the dependence of the magnetisation on the local concentration of magnetite NPs. This model was used to simulate both soap films and magnetic soap films. It was found that interfacial flows can dominate the rate of film thinning, depending on the strength of Marangoni stresses created by the interfacial flows. Accounting for variations in the magnetite NP concentration can also have a notable effect on film thickness predictions; however, this is less important for high Péclet number ($Pe_c$) flows. A capillary suction boundary condition was proposed with the aim of simulating marginal regeneration. An exchange between thin and thick fluid at the film boundary was observed in the simulations for soap films, and the simulated thin patches of fluid flowed from the film boundary towards the film centre, in agreement with the marginal regeneration phenomena observed in the experiments. The derived model was also able to reproduce several features observed in experimentally investigated magnetic soap films, such as faster magnetically driven drainage flows in the direction of increasing magnetic field intensity with larger concentrations of magnetite NPs and reverse flows that separate from the film boundary and flow into the interior of the film. This work advances the state of magnetic soap film modelling and is an important step towards more accurately modelling soap films in which marginal regeneration is occurring.