2.1. Outline of the modelling procedure

Model development necessitates combined consideration of (i) the fluid dynamics of the liquid layer, (ii) the mass exchange dynamics at the interface and (iii) the conservation equations of surfactant. The derivations, although not particularly difficult, are involved and with rather numerous intermediate steps, so the following outline may be helpful.

The fluid dynamics of the liquid layer is described by the Navier–Stokes equation in the low Reynolds number limit, with free-surface boundary conditions applied along the interface. The dynamic boundary condition contains the local values of surface tension (in the capillary stress term) and its spatial gradients (in the Marangoni stress term). The surface tension is considered a sole function of the local monolayer concentration, $\varGamma _m$ , and thus its values – as well as those of surface elasticity – are derived from an equation of state that follows from Gibbs’ theory, with key novelty the inclusion of the intrinsic compressibility of the adsorbed monolayer.

Surfactant concentration on the interface is described by two terms, $\varGamma _m$ , which gives the surface concentration of the adsorbed monolayer and $\varGamma$ , which, includes both the monolayer surfactant and the surfactant contained in the sublayer reservoir that forms during monolayer collapse. Two conservation equations are thus formulated to follow the spatio-temporal evolution of the above surface concentrations, which include the contributions of interfacial expansion/contraction, of convection and diffusion along the interface and of mass exchange with the environment.

The mass exchange terms are crucial inputs to the surfactant conservation equations, and they consist of the following two contributions. The flux $j_b$ , which describes adsorption/desorption between the monolayer and the bulk, and is modelled by Langmuir kinetics in terms of an adsorption rate and an equilibrium constant. The flux, $j_r$ , which describes internal exchange between the monolayer and the reservoir, and is modelled by a Langmuir-like term that accounts for monolayer replenishment from the reservoir during interface expansion and by a singular term that limits maximum monolayer concentration below $\varGamma _{max}$ during interface contraction. In the following paragraphs, each of the aforementioned dynamic components of the model is considered in detail.

2.2. Dynamics of the liquid layer

The alveolus is modelled as a spherical cap with an opening angle $2 \theta _0$ ( $\theta _0=30^\circ$ for a normal alveolus), as shown in figure 1(a). During oscillation, the opening angle remains constant, i.e. the alveolus retains a self-similar shape (Tsuda, Henry & Butler Reference Tsuda, Henry and Butler2008). Although there are alternative choices for the exact alveolar shape, the truncated sphere remains the most common model geometry in the literature (Kolanjiyil & Kleinstreuer Reference Kolanjiyil and Kleinstreuer2019). The alveolar rim is taken to have a semicircular cross-section of small but finite radius $r_0$ (figure 1 b). Thus, the local septal thickness is ${\sim}2r_0$ . Due to the finite thickness, water and surfactant accumulate on the rim and local values of film thickness and surfactant concentration, which continuously vary during alveolar oscillations, enter the boundary conditions (Bouchoris & Bontozoglou Reference Bouchoris and Bontozoglou2022).

Figure 1. (a) A sketch of the spherical cap with the main problem parameters and (b) magnification of the rim with the local parameter values (red colour marks the hydrophobic end of surfactant molecules).

The liquid layer is Newtonian (Grotberg Reference Grotberg2011), with density $\rho$ and viscosity $\mu$ . The flow is taken to be symmetric in the circumferential direction and is analysed in a spherical coordinate system $(r,\theta,\phi)$ located at the centre of the cap. Thus, the flow field is described by the velocity vector $\boldsymbol{u}=(u_{r}(r,\theta,t),u_{\theta }(r,\theta,t),0)$ and the pressure distribution $p(r,\theta,t)$ . The air–liquid interface is located at $r=R(t)-h(\theta,t)$ , where $h(\theta,t)$ is the liquid film thickness. It is further characterised by its surface tension, $\sigma (\theta,t)$ , which varies in space and time as it is determined by the local concentration of the surfactant monolayer.

Following standard practice in the literature (Podgorski & Gradon Reference Podgorski and Gradon1993; Haber et al. Reference Haber, Butler, Brenner, Emmanuel and Tsuda2000; Zelig & Haber Reference Zelig and Haber2002; Wei et al. Reference Wei, Fujioka, Hirschl and Grotberg2005; Kang et al. Reference Kang, Chugunova, Nadim, Waring and Walther2018), the flow is described by the continuity and the quasi-steady Stokes equation, subject to the no-slip boundary condition on the wall and the kinematic and dynamic boundary conditions along the interface. The dynamic condition includes components normal and tangential to the interface, the former proportional to surface tension and local curvature (capillary stress) and the latter proportional to the local spatial gradient of surface tension (Marangoni stress). The equations are further simplified by a lubrication assumption, which is based on the large scale difference between the $r\hbox{-}$ and $\theta \hbox{-}$ directions ( $h\lt \lt R$ ). The final result is the following evolution equation that describes the spatio-temporal variation of the liquid film thickness, driven by capillary and Marangoni stresses (respectively, the first and second terms in the parenthesis). The details of the mathematical derivation are presented in a previous paper (Bouchoris & Bontozoglou Reference Bouchoris and Bontozoglou2022) and are also discussed by Kang et al. (Reference Kang, Chugunova, Nadim, Waring and Walther2018)

(2.1) \begin{equation} \frac {\partial h}{\partial t} + \frac {2h\dot {R}}{R} - \frac {1}{R^2\sin \theta }\frac {\partial }{\partial \theta }\left (\frac {h^3\sin \theta }{3\mu }\frac {\partial p}{\partial \theta }-\frac {h^2\sin \theta }{2\mu }\frac {\partial \sigma }{\partial \theta }\right) = 0. \end{equation}

It is noted that, in the lubrication approximation, pressure is uniform in the radial direction and can be described by the following expression, using the normal boundary condition on the interface:

(2.2) \begin{equation} p(\theta,t) = p_{air} - \sigma \left [\frac {2}{R}+\frac {2\,h}{R^2}+\frac {1}{R^2\sin \theta }\frac {\partial }{\partial \theta }\left (\sin \theta \,\frac {\partial h}{\partial \theta }\right)\right]. \end{equation}

The term in parenthesis in (2.1) is equal to $-Q(\theta,t)/2\pi$ , where

(2.3) \begin{equation} Q(\theta,t)=2\pi R\,\sin \theta \int _{R-h}^{R} u_{\theta }\,{\rm d}r=2\pi \sin \theta \left (-\frac {h^3}{3\mu }\frac {\partial p}{\partial \theta }+\frac {h^2}{2\mu }\frac {\partial \sigma }{\partial \theta }\right)\,, \end{equation}

is, in the lubrication limit, the volumetric flow rate along the entire $\phi \hbox{-}$ circumference. Thus, (2.1) may be rewritten as

(2.4) \begin{equation} \frac {\partial }{\partial t}\big (2\pi R^2\sin \theta \,h\big) + \frac {\partial Q}{\partial \theta } = 0, \end{equation}

which is intuitively more evident as the mass balance of liquid in a differential control volume along the entire $\phi \hbox{-}$ periphery at the location $\theta$ .

The derivations leading to (2.1) – as well as the entire paper – neglect the effect of the dynamics of airflow inside the alveolus during the breathing cycle. The following simple scaling argument, based on the continuity of shear stress and velocity at the air–liquid interface, shows that the effect is indeed negligible. More specifically, if $u_{s,a}$ is the liquid interfacial velocity tentatively imposed by gas shear

(2.5) \begin{equation} \mu _a\left (\frac {\partial u}{\partial y}\right)_a\sim \mu \frac {\partial u}{\partial y}\Rightarrow \mu _\alpha \frac {u_a-u_{s,a}}{\bar {R}}\sim \mu \frac {u_{s,a}-0}{\bar {H}}\Rightarrow u_{s,a}\sim u_a\frac {\mu _a}{\mu }\frac {\bar {H}}{\bar {R}}. \end{equation}

Given that $\mu _a/\mu \sim 1:50$ and $\bar {H}/\bar {R}\sim 0.3:100$ , the liquid velocity imposed by air shear is four orders of magnitude smaller than the typical air velocity inside the alveolus, i.e. two to three orders of magnitude smaller than the shearing velocity on the liquid film triggered by Marangoni stresses (as shown by Bouchoris & Bontozoglou (Reference Bouchoris and Bontozoglou2022), the shearing velocity is one to two orders of magnitude smaller than the radial velocity of the alveolar wall).

An assumption, implicit in the above formulation, is the neglect of the potential role of surfactant aggregates in the liquid film dynamics. The dynamics of the aggregates themselves is probably related to the kinematics of surfactant exchange with the interface and is not in direct interest to the present study. However, given the relative sizes of the aggregates ( ${\sim} 100\,\unicode{x03BC}\rm m$ ) and the liquid film ( ${\sim} 200{-}1000\,\unicode{x03BC}\rm m$ ), it is reasonable to expect that, with increasing number density, the aggregates will have a non-negligible effect on the liquid dynamics. Possibilities for including this effect in more refined simulations are mentioned in the Concluding Remarks.

2.3. Dynamics of the surfactant

The pulmonary surfactant is a mixture of lipids and proteins, which are practically insoluble in water. It is indicative that the critical micelle concentration (CMC) of the major phospholipids (i.e. the maximum concentration of monomer beyond which monomers combine to form micelles) is $10^{-9}{-}10^{-10}\,\textrm{M}$ (Marsh & King Reference Marsh and King1986). As a result, a portion of the surfactant populates the interface with local surface concentration $\varGamma (\theta,t)$ , while the excess amount resides in the bulk in the form of aggregates (Zuo et al. Reference Zuo, Veldhuizen, Neumann, Petersen and Possmayer2008; Bykov et al. Reference Bykov, Milyaeva, Isakov, Michailov, Loglio, Miller and Noskov2021). There is ample experimental evidence (Baoukina & Tieleman Reference Baoukina and Tieleman2011; Parra & Perez-Gil Reference Parra and Perez-Gil2015; Bouchoris & Bontozoglou Reference Bouchoris and Bontozoglou2021) indicating that the adsorption/desorption of surfactant to/from the interface occurs through its direct interaction with these aggregates. The pulmonary surfactant mixture is currently modelled by a single generic surfactant, which mimics the dynamic behaviour of the actual preparation (Wüstneck et al. Reference Wüstneck, Perez-Gil, Wüstneck, Cruz, Fainerman and Pison2005).

The description of the interfacial surfactant composition becomes more involved when the alveolar oscillations are strong enough to include the stages of monolayer collapse and subsequent replenishment. The surface concentration, $\varGamma$ , now consists of the sum of the monolayer and reservoir concentrations, $\varGamma =\varGamma _{m}+\varGamma _{r}$ . During compression, the monolayer concentration does not increase beyond the maximum value that corresponds to the minimum attainable surface tension. Thus, $\varGamma _{m}\rightarrow \varGamma _{max}$ , where $\sigma (\varGamma _{max})=\sigma _{min}$ , and any further increase of $\varGamma$ results in an equal increase of $\varGamma _{r}$ . It is important to re-iterate that, capillary and Marangoni stresses vary with $\varGamma _m$ and not with $\varGamma$ . Figure 2 provides a sketch of the interfacial region and the surfactant exchange processes during compression (figure 2 a) and during expansion (figure 2 b).

Figure 2. Schematic of the interface with the monolayer and the reservoir protrusions, and the bulk with surfactant aggregates. (a) Compression that leads to monolayer collapse and growth of the reservoir. (b) Re-expansion, when the monolayer concentration falls below $\varGamma _{max}$ and is replenished by surfactant from the reservoir. Red colour marks the hydrophobic end of surfactant molecules and bars in (b) (not to scale) indicate the sizes of the monomer, the surface layer and the aggregates as measured by Xu et al. (Reference Xu, Yang and Zuo2020).

The local value of surface tension, $\sigma$ , along the interface is a function of the local monolayer concentration, i.e. $\sigma =\sigma (\varGamma _m)$ , and is described by the equation of state to be developed shortly. Considering a space coordinate, $s$ , along the interface, the variation of surface tension and monolayer concentration are related by the chain rule, i.e. $\partial \sigma /\partial s=({\rm d}\sigma /{\rm d}\varGamma _m)(\partial \varGamma _m/\partial s)$ . The sensitivity of $\sigma$ to $\varGamma _m$ is expressed by the Gibbs elasticity, $E$ , where

(2.6) \begin{equation} E = -\frac {{\rm d}\sigma }{{\rm d}\ln \varGamma _m} = -\varGamma _m \frac {{\rm d}\sigma }{{\rm d}\varGamma _m}. \end{equation}

A novel feature of the surfactant model is the inclusion of an intrinsic compressibility, $\alpha$ , of the adsorbed monolayer, as defined and justified by Fainerman, Miller & Kovalchuk (Reference Fainerman, Miller and Kovalchuk2002) and Kovalchuk et al. (Reference Kovalchuk, Loglio, Fainerman and Miller2004, Reference Kovalchuk, Miller, Fainerman and Loglio2005). More specifically, the molar surface area of the monolayer, $\varOmega$ , is taken to vary linearly with surface pressure, $\varPi =\sigma _0-\sigma$ , according to the relation

(2.7) \begin{equation} \varOmega = \varOmega _{0} (1-\alpha \varPi), \end{equation}

where $\sigma _0$ is the surface tension of pure water and $\varOmega _0$ is the molar area at zero surface pressure. This correction is particularly significant for dense monolayers and was recently shown to offer quantitative agreement of the model dynamics with laboratory measurements using actual pulmonary preparations (Saad et al. Reference Saad, Neumann and Acosta2010; Xu et al. Reference Xu, Yang and Zuo2020).

In terms of $\varOmega$ , the monolayer coverage, $\gamma$ , is $\gamma =\varGamma _m \varOmega$ , and thus the surface concentration at interfacial saturation, $\varGamma _{m,\infty }$ (which corresponds to $\gamma =1$ ), varies with surface pressure, and is given by

(2.8) \begin{equation} \varGamma _{m,\infty } = \frac {1}{\varOmega } = \frac {1}{\varOmega _{0}(1-\alpha \varPi)}. \end{equation}

Using Gibbs’ theory, which is valid for a Gibbs dividing surface and an ideal bulk phase, the following equation of state is derived (Kovalchuk et al. Reference Kovalchuk, Loglio, Fainerman and Miller2004), connecting surface pressure to interfacial coverage:

(2.9) \begin{equation} \frac {\varPi \varOmega _0}{\mathcal{R}\mathcal{T}}\left (1-\alpha \frac {\varPi }{2}\right) = -\ln (1-\gamma),\end{equation}

where $\mathcal{R}$ is the universal gas constant and $\mathcal{T}$ the absolute temperature.

The surface concentration of surfactant varies continuously during breathing as a result of the periodic expansion and contraction of the interface. The variation of $\varGamma _m$ from the equilibrium value, $\varGamma _{eq}$ (corresponding to a non-oscillating alveolus for which $\varGamma _{eq}\equiv \varGamma _{m,eq}$ ) causes a mass flux, $j_b$ , of surfactant between the monolayer and the bulk. Assuming a Langmuir adsorption isotherm, equilibrium is given by the expression

(2.10) \begin{equation} KC_{10}=\frac {\varGamma _{eq}}{\varGamma _{eq,\infty }-\varGamma _{eq}},\end{equation}

where $\varGamma _{eq,\infty }$ is the saturation value at equilibrium conditions, $K$ is the equilibrium constant and $C_{10}$ the CMC of the monomer in the bulk. Thus, mass exchange with the bulk is expressed as follows in terms of an adsorption rate, $k_{ads}$ :

(2.11) \begin{equation} j_b = k_{ads}C_{10}\left [(\varGamma _{m,\infty }-\varGamma _m)-\frac {\varGamma _m}{K\,C_{10}} \right], \end{equation}

while for an insoluble surfactant $j_b=0$ . It is noted that the above expression for mass exchange is based on kinetic resistance at the interface (rather than diffusion in the bulk), an assumption strongly supported by the literature (Ingenito et al. Reference Ingenito, Mark, Morris, Espinosa, Kamm and Johnson1999; Saad et al. Reference Saad, Neumann and Acosta2010).

Mass exchange between the monolayer and the reservoir, $j_r$ , changes direction during the breathing cycle, because the reservoir stores surfactant during compression (so that $\varGamma _m$ does not exceed $\varGamma _{max}$ ) and releases it back to the monolayer during re-expansion of the alveolus (when $\varGamma _m$ falls below $\varGamma _{max}$ ). This dual role is described by the following expression:

(2.12) \begin{equation} j_{r} = k_r ( \varGamma _{m,\infty } - \varGamma _{m})\left (\varGamma - \varGamma _m\right) - \frac {k_{rc}}{\left (\varGamma _{max}-\varGamma _m\right)^2}. \end{equation}

The first term in (2.12) gives the rate of replenishment of the monolayer during expansion and is proportional to the potential for replenishment, $ ( \varGamma _{m,\infty } - \varGamma _{m})$ , and to the available surfactant concentration in the local reservoir, $ (\varGamma -\varGamma _m)$ . The second term limits the monolayer concentration to be below the maximum value $\varGamma _{max}$ during compression. Parameter $k_{rc}$ may be viewed as a measure of the stiffness of the dense monolayer, in the sense that the smaller $k_{rc}$ the closer the monolayer concentration is to its limiting value without collapsing and the more abrupt the collapse. By choosing a very small value $k_{rc}=10^{-5}$ , the contribution of the second term in (2.12) becomes significant only when $\varGamma _m$ approaches close enough to $\varGamma _{max}$ .

2.4. Mass balances of surfactant

The aforementioned variables, $\varGamma$ and $\varGamma _m$ , are functions of location and time, and their evolution is determined by appropriate mass conservation equations. Such equations must include convection and diffusion along the interface and mass exchange with the bulk and the reservoir, and must also take into account the effect of expansion/contraction of the interfacial area on surface concentration. It is noted that there is no need for a mass balance in the bulk, because the typical surfactant loading is much higher than the CMC of the monomer. Thus, the bulk is always saturated with monomer and the effect of bulk diffusion is negligible.

First, the case of weak alveolar oscillations – corresponding to shallow breathing manoeuvres – is considered. In this case, the surface concentration never approaches $\varGamma _{max}$ , so $\varGamma _m\equiv \varGamma$ . Thus, a single conservation equation is required, which – when simplified by a lubrication approximation – leads to the following result (Kang et al. Reference Kang, Chugunova, Nadim, Waring and Walther2018; Bouchoris & Bontozoglou Reference Bouchoris and Bontozoglou2022):

(2.13) \begin{equation} \begin{aligned} &\frac {\partial \varGamma }{\partial t} + \frac {2\varGamma \dot {R}}{R} - \frac {1}{R^2\sin \theta }\frac {\partial }{\partial \theta }\left (\frac {\varGamma h^2\sin \theta }{2\mu }\frac {\partial p}{\partial \theta }-\frac {\varGamma h\sin \theta }{\mu }\frac {\partial \sigma }{\partial \theta }\right)\\ &\qquad - \frac {\mathcal{D}_s}{R^2\sin \theta }\frac {\partial }{\partial \theta }\left (\sin \theta \,\frac {\partial \varGamma }{\partial \theta }\right) = j_b. \end{aligned} \end{equation}

Here, $\mathcal{D}_s$ is the surface diffusivity of surfactant and the terms in the two parentheses represent, respectively, the contributions of convection and surface diffusion. Eqaution (2.13) may take the following condensed form:

(2.14) \begin{equation} \frac {\partial }{\partial t} \big(2\pi R^2\sin \theta \;\varGamma \big) + \frac {\partial Q_{\varGamma }}{\partial \theta } = 2\pi R^2 \sin \theta \, j_b, \end{equation}

where $Q_{\varGamma }(\theta,t)$ is the interfacial flow rate of surfactant along the entire $\phi -$ circumference

(2.15) \begin{equation} Q_{\varGamma }(\theta,t)=2\pi R\,\sin \theta \left (u_s\varGamma -\frac {\mathcal{D}_s}{R} \frac {\partial \varGamma }{\partial \theta }\right),\end{equation}

and $u_s$ is the interfacial water velocity

(2.16) \begin{equation} u_s\equiv u_{\theta }(R-h,\theta,t)=-\frac {1}{2\mu }\frac {\partial p}{\partial \theta }\frac {h^2}{R}+\frac {1}{\mu }\frac {\partial \sigma }{\partial \theta }\frac {h}{R}. \end{equation}

In order to describe deep-breathing manoeuvres, i.e. alveolar oscillations that include monolayer collapse and replenishment, the following modifications are necessary. As noted previously, $\sigma =\sigma (\varGamma _m)$ and $E=E(\varGamma _m)$ , i.e. capillary and Marangoni forcing are both driven by the monolayer concentration and not the total surface value, $\varGamma$ . Thus, convective and diffusive fluxes have as a driving force the monolayer concentration $\varGamma _m$ . However, the reservoir extensions are attached to the monolayer (Parra & Perez-Gil Reference Parra and Perez-Gil2015), and consequently are transported with it. Thus, for each mol of monolayer surfactant transported by convection or diffusion $\varGamma _r/\varGamma _m$ also follow, or equivalently a total of $(\varGamma _r+\varGamma _m)/\varGamma _m=\varGamma /\varGamma _m$ is transported. As a result, the interfacial flow rate, $Q_{\varGamma }$ , now becomes

(2.17) \begin{equation} \begin{aligned} Q_{\varGamma }(\theta,t)=2\pi R\,\sin \theta \left [u_s(\varGamma _m)\,\varGamma _m-\frac {\mathcal{D}_s}{R} \frac {\partial \varGamma _m}{\partial \theta }\right]\left (\frac {\varGamma }{\varGamma _m}\right) \\ =2\pi R\,\sin \theta \left [u_s(\varGamma _m)\,\varGamma -\frac {\mathcal{D}_s}{R} \frac {\partial \varGamma _m}{\partial \theta }\left (\frac {\varGamma }{\varGamma _m}\right)\right]. \end{aligned} \end{equation}

Implementation of the above modifications leads to the following conservation equation for the total surface concentration of surfactant:

(2.18) \begin{equation} \begin{aligned} \frac {\partial \varGamma }{\partial t} + \frac {2\varGamma \dot {R}}{R} - \frac {1}{R\sin \theta }\frac {\partial }{\partial \theta }\left [\sin \theta \, u_s(\varGamma _m)\,\varGamma -\sin \theta \,\frac {\mathcal{D}_s}{R} \frac {\partial \varGamma _m}{\partial \theta }\left (\frac {\varGamma }{\varGamma _m}\right)\right] = j_b, \end{aligned} \end{equation}

with the exchange rate with the bulk given by (2.11). It is noted that the internal exchange rate $j_r$ , given by (2.12), results only in a redistribution of surfactant between the monolayer and the reservoir, and therefore does not enter the conservation equation for the total surface concentration. However, it is critical for the correct modelling of the evolution of the monolayer concentration, $\varGamma _m$ , which is described by the following expression:

(2.19) \begin{equation} \begin{aligned} \frac {\partial \varGamma _m}{\partial t} + \frac {2\varGamma _m\dot {R}}{R} - \frac {1}{R\sin \theta }\frac {\partial }{\partial \theta }\left [\sin \theta \,u_s(\varGamma _m)\,\varGamma _m-\sin \theta \,\frac {\mathcal{D}_s}{R} \frac {\partial \varGamma _m}{\partial \theta }\right] = j_b + j_r. \end{aligned} \end{equation}

2.5. Initial and boundary conditions

The evolution of liquid film thickness and surfactant surface concentration as described by (2.1), (2.18) and (2.19) – or (2.1) and (2.13) for the simpler case of weakly nonlinear oscillations – is supplemented by appropriate initial and boundary conditions. A convenient starting point for all the computations to be presented next is the equilibrium structure of a motionless alveolus. At equilibrium, both the interface surfactant concentration and the capillary pressure are uniform in $\theta$ , and the interface attains a spherical shape with surface tension $\sigma _{eq}=\sigma _{eq}(\varGamma _{eq})$ . Then, (2.1) has a trivial solution for the liquid film thickness, $H(\theta)$ , which is a linear function of $\cos \theta$ determined by two parameters, the mean thickness, $\bar {H}$ , and the thickness over the rim, $H(\theta _0)\equiv H_0$ (Bouchoris & Bontozoglou Reference Bouchoris and Bontozoglou2022). The computations with this initial condition indicate that in most cases a single oscillation cycle is sufficient to reach the periodic solution.

The existence of a continuous, and rather thick, liquid layer over the sharp rim ( $H_0=0.14\;\unicode{x03BC} \textrm {m}$ as estimated from the data of Bastacky et al. (Reference Bastacky, Lee, Goerke, Koushafar, Yager, Kenaga, Speed, Chen and Clements1995)) is attributed to the existence of repulsive disjoining pressure forces (Starov & Ivanov Reference Starov and Ivanov2004), which develop between the interface, the epithelium and surfactant aggregates squeezed in between. This picture is further supported by recent measurements (Xu et al. Reference Xu, Yang and Zuo2020) that estimate the combined aggregate and monolayer thickness to be of the order of $120{-}150\; \textrm {nm}$ , i.e. in quantitative agreement with the data of Bastacky et al. (Reference Bastacky, Lee, Goerke, Koushafar, Yager, Kenaga, Speed, Chen and Clements1995). Apart from determining $H_0$ , the disjoining pressure has no role in the dynamics, because its characteristic time scale is very large, even compared with the breathing frequency (Oron, Davis & Bankoff Reference Oron, Davis and Bankoff1997; Zelig & Haber Reference Zelig and Haber2002). The non-zero value of film thickness over the rim was shown by Bouchoris & Bontozoglou (Reference Bouchoris and Bontozoglou2022) to play a key role in the development of tangential motion during the oscillation cycle.

Boundary conditions are applied at the apex, $\theta =\pi /2$ , and along the rim of the alveolus at $\theta =\theta _0$ . Symmetry is assumed at the alveolar apex, and thus

(2.20) \begin{equation} \frac {\partial h}{\partial \theta }=\frac {\partial \varGamma }{\partial \theta }=\frac {\partial \varGamma _m}{\partial \theta }=0 \quad \textrm {at} \quad \theta =\pi /2. \end{equation}

It has been argued (Bouchoris & Bontozoglou Reference Bouchoris and Bontozoglou2022) that the boundary conditions along the rim are critical for correct modelling of the dynamics, as they dictate the gradients with the interior that drive the flow and transport phenomena. In the present investigation, the boundary conditions proposed by Bouchoris & Bontozoglou (Reference Bouchoris and Bontozoglou2022) are used, which amount to simple mass balances that equate water and surfactant accumulation over the rim to the respective fluxes from the interior at $\theta _0$ .

In terms of the volumetric flow rate of water around the circumference of the rim, $Q(\theta _0,t)$ , the mass balance of water takes the form

(2.21) \begin{equation} \begin{gathered} -2\,Q(\theta _0,t) = \frac {d}{{\rm d}t}\left [2\pi R \sin \theta _0 \left (\pi \frac {(r_0+h_0)^2}{2}-\pi \frac {r_0^2}{2}\right)\right]. \end{gathered} \end{equation}

Similar mass balances for the total and monolayer surfactant, taking into account convection and diffusion at the interface and exchange with the bulk and the reservoir, lead to the following expressions:

(2.22) \begin{equation} -2 Q_{\varGamma }(\theta _0,t) + 2\pi R \sin \theta _0 \left.\pi (r_0+h_0) j_b\right |_{\theta _0} = \frac {d}{{\rm d}t}\big [2\pi R\sin \theta _0\,\pi \left (r_0+h_0\right)\varGamma _0\big], \end{equation}

(2.23) \begin{equation} -2 Q_{\varGamma m}(\theta _0,t) + 2\pi R \sin \theta _0 \left.\pi (r_0+h_0) (j_b+j_r)\right |_{\theta _0} = \frac {d}{{\rm d}t}\big [2\pi R\sin \theta _0\,\pi \left (r_0+h_0\right)\varGamma _{m0}\big]. \end{equation}

The multiplier two in. (2.21), (2.22) and (2.23) accounts for flux coming to the rim from both sides of the mid-alveolar wall and is appropriate for alveoli in close packing (type A alveoli), as they occur in the lowest respiratory bronchioles. The minus sign indicates flow in the negative $\theta \hbox{-}$ direction, i.e. towards the rim.

The application of mass balances only – and not flow equations – over the rim is based on the assumption that the rim dynamics does not influence the problem, or, in other words, that it is totally enslaved to the dynamics of the alveolus. This assumption has been justified by an order-of-magnitude analysis, which revealed that Marangoni forcing is much faster than capillary drainage (Bouchoris & Bontozoglou Reference Bouchoris and Bontozoglou2022). Thus, the flow induced inside the alveolus by Marangoni stresses dominates the system dynamics.

3.1. The final equations in dimensionless form

The characteristic scales used to non-dimensionalise the problem variables are taken from the equilibrium conditions. Thus, we consider a motionless alveolar cap of radius $\bar {R}$ , coated by a liquid film whose mean thickness is $\bar {H}$ . With these choices, the lubrication parameter is formally defined as

(3.1) \begin{equation} \epsilon =\frac {\bar {H}}{\bar {R}}. \end{equation}

The liquid is loaded by surfactant aggregates and equilibrates with an adsorbed monolayer of surface concentration $\varGamma _{\textit{eq}}$ , which results in surface tension $\sigma _{eq}$ and surface elasticity $E_{eq}=-\varGamma _{eq} \left. ({\rm d}\sigma /{\rm d}\varGamma) \right |_{\varGamma =\varGamma _{eq}}$ . Finally, the characteristic time is the breathing period, $T$ .

With the above scales, the following dimensionless variables are defined: $h^*=h/\bar {H}$ , $H^*=H/\bar {H}$ , $R^*=R/\bar {R}$ , $r_0^*=r_0/\bar {R}$ , $t^*=t/T$ , $\varGamma ^*=\varGamma /\varGamma _{eq}$ , $\sigma ^*=\sigma /\sigma _{eq}$ , $E^*=E/E_{eq}$ , $p^*=p\,\bar {R}^2/(\sigma _{eq}\bar {H})$ , $j_b^*=j_b T/\varGamma _{eq}$ and $j_r^*=j_r T/\varGamma _{eq}$ . The governing evolution equations are transformed as follows in terms of the dimensionless variables:

(3.2) \begin{equation} \frac {\partial h^*}{\partial t^*} + \frac {2h^*\dot {R^*}}{R^*} + \frac {1}{R^{*2}\sin \theta } \frac {\partial Q^*}{\partial \theta } = 0,\end{equation}

(3.3) \begin{equation} \frac {\partial \varGamma ^*}{\partial t^*} + \frac {2\varGamma ^*\dot {R^*}}{R^*}+ \frac {1}{R^{*2}\sin \theta }\frac {\partial Q_{\varGamma }^*}{\partial \theta }-j_b^* = 0,\end{equation}

and

(3.4) \begin{equation} \frac {\partial \varGamma _m^*}{\partial t^*} + \frac {2\varGamma _m^*\dot {R^*}}{R^*}+ \frac {1}{R^{*2}\sin \theta }\frac {\partial Q_{\varGamma m}^*}{\partial \theta }-(j_b^*+j_r^*) = 0. \end{equation}

Terms $Q^*,\;Q_{\varGamma }^*, \;Q_{\varGamma m}^*$ are respectively the dimensionless volumetric flow rates of water, and total and monolayer surfactant, the former scaled by $\bar {Q}=2\pi \bar {R}^2\bar {H}/T$ and the two latter by $\bar {Q}_{\varGamma,eq}=2\pi \bar {R}^2 \varGamma _{eq}/T$ .

Three trivial algebraic manipulations are undertaken in order to derive the final evolution equations from ((3.2)–(3.4)). First, the independent spatial variable $\theta$ is replaced by $x=-\cos \theta$ . Second, the fluxes of water and surfactant are expanded according to the expressions in (2.1), (2.15) and (2.17), using also (2.2) and (2.16) for the pressure and the surface velocity. Third, to avoid symbol overload, the dimensionless variables are written from now on without the asterisk, e.g. $h$ signifies $h^*$ , $\varGamma$ signifies $\varGamma ^*$ , $R$ signifies $R^*$ etc. There should be no confusion with this last change, since only dimensionless variables appear in the rest of the paper.

Thus, the final evolution equations are

(3.5) \begin{equation} \frac {\partial h}{\partial t} + \frac {2h\dot {R}}{R} + \frac {1}{R^2} \frac {\partial Q}{\partial x} = 0,\end{equation}

(3.6) \begin{equation} \frac {\partial \varGamma }{\partial t} + \frac {2\varGamma \dot {R}}{R}+ \frac {1}{R^2}\frac {\partial Q_{\varGamma }}{\partial x}-j_b = 0,\end{equation}

and

(3.7) \begin{equation} \frac {\partial \varGamma _m}{\partial t} + \frac {2\varGamma _m\dot {R}}{R}+ \frac {1}{R^2}\frac {\partial Q_{\varGamma m}}{\partial x}-(j_b+j_r) = 0,\end{equation}

which are solved for initial conditions

(3.8) \begin{equation} h(x,0)=H(x), \quad \varGamma (x,0)=\varGamma _m(x,0)=1,\end{equation}

and subject to the following boundary conditions at $x_0=-\cos \theta _0$ , with subscript $0$ signifying value at the rim:

(3.9) \begin{equation} -2\,Q(x_0,t)=\sqrt {1-x_0^2}\;\pi R\left [\left (r_0+\epsilon h_0\right)\frac {{\rm d}h_0}{{\rm d}t}+h_0\left (r_0+\epsilon \frac {h_0}{2}\right)\frac {\dot {R}}{R}\right],\end{equation}

(3.10) \begin{equation} -2\,Q_{\varGamma }(x_0,t)=\sqrt {1-x_0^2}\;\pi R\left [(r_0+\epsilon h_0)\left (\frac {{\rm d}\varGamma _0}{{\rm d}t}-j_{b0}\right)+\epsilon \varGamma _0\frac {{\rm d}h_0}{{\rm d}t}+\varGamma _0(r_0+\epsilon h_0)\frac {\dot {R}}{R}\right],\end{equation}

(3.11) \begin{align}\notag &-2\,Q_{\varGamma m}(x_0,t)\\ &\quad =\sqrt {1-x_0^2}\;\pi R\left [(r_0+\epsilon h_0)\left (\frac {{\rm d}\varGamma _{m0}}{{\rm d}t}-j_{b0}-j_{r0}\right)+\epsilon \varGamma _{m0}\frac {{\rm d}h_0}{{\rm d}t}+\varGamma _{m0}(r_0+\epsilon h_0)\frac {\dot {R}}{R}\right]. \end{align}

The symmetry boundary conditions at $\theta =\pi$ are now trivially satisfied, provided the derivatives $\partial h/\partial x,\;\partial \varGamma /\partial x,\;\partial \varGamma _m/\partial x$ are finite at $x=1$ .

The fluxes that appear in ((3.5)–(3.7)) are given in expanded form by the following expressions:

(3.12) \begin{equation} Q=\epsilon ^3\frac {Ca^{-1}}{3R^4}\sigma (\varGamma _m) h^3 \big(1-x^2 \big) \frac {\partial }{\partial x}(2h+g)-\epsilon \frac {Ma}{2R^2}h^2 \big(1-x^2 \big) \frac {E(\varGamma _m)}{\varGamma _m}\frac {\partial \varGamma _m}{\partial x},\end{equation}

(3.13) \begin{equation} \begin{aligned} Q_{\varGamma }&=\epsilon ^3\frac {Ca^{-1}}{2R^4}\varGamma \sigma (\varGamma _m) h^2 \big(1-x^2 \big) \frac {\partial }{\partial x}(2h+g) - \epsilon \frac {Ma}{R^2}\varGamma h \big(1-x^2 \big) \frac {E(\varGamma _m)}{\varGamma _m}\frac {\partial \varGamma _m}{\partial x} \\ &\quad-\frac {Pe_s^{-1}}{R^2} \big(1-x^2 \big)\frac {\partial \varGamma _m}{\partial x}\left (\frac {\varGamma }{\varGamma _m}\right),\end{aligned} \end{equation}

(3.14) \begin{equation} \begin{aligned} Q_{\varGamma m}&=\epsilon ^3\frac {Ca^{-1}}{2R^4}\varGamma _m\,\sigma (\varGamma _m) h^2 \big(1-x^2 \big) \frac {\partial }{\partial x}(2h+g) \\ &\quad -\epsilon \frac {Ma}{R^2}\varGamma _m\, h \big(1-x^2 \big) \frac {E(\varGamma _m)}{\varGamma _m}\frac {\partial \varGamma _m}{\partial x}- \frac {Pe_s^{-1}}{R^2}\big(1-x^2 \big)\frac {\partial \varGamma _m}{\partial x},\end{aligned} \end{equation}

(3.15) \begin{equation} j_b = St_b \left [ (\varGamma _{m,\infty }-\varGamma _m)-\frac {\varGamma _m}{KC_{10}}\right],\end{equation}

(3.16) \begin{equation} j_r = St_r (\varGamma _{m,\infty }-\varGamma _m) (\varGamma -\varGamma _m) - St_{rc}\frac {1}{\left (\varGamma _{max}-\varGamma _m\right)^2},\end{equation}

where, as a reminder, the properties $\sigma$ and $E$ in (3.12)–(3.14) are written with the dependence on monolayer concentration made explicit. It is noted that the function $g(x,t)$ , which appears in these expansions, is a substitution for the expression

(3.17) \begin{equation} g(x,t)=\frac {\partial }{\partial x}\left (\big(1-x^2 \big)\frac {\partial h}{\partial x}\right). \end{equation}

This substitution is necessary for the finite-element solution of the problem, as the introduction of the new unknown function $g(x,t)$ eliminates the higher than second-order derivatives in $h$ . Finally, the dimensionless numbers that appear in (3.12)–(3.16) are defined as follows:

(3.18) \begin{equation} \begin{aligned} Ca^{-1}=\frac {\sigma _{eq}}{\mu (\bar {R}/T)}\;, \quad Ma=\frac {E_{eq}}{\mu (\bar {R}/T)}\;, \quad Pe_s^{-1}=\frac {\mathcal{D}_s T}{\bar {R}^2}\;,\\ St_b=\frac {T}{1/(k_{ads}C_{10})}\;, \quad St_r=\frac {T}{1/(k_r\varGamma _{eq})}\;, \quad St_{rc}=\frac {T}{\varGamma _{eq}^3/k_{rc}}, \end{aligned} \end{equation}

and may be viewed as ratios of the breathing period to the characteristic times of capillary, Marangoni, diffusion, adsorption and replenishment phenomena.

3.2. Numerical solution

The evolution equations ((3.5)–(3.7)) are discretised in space with a Galerkin finite-element method. The nodes of the independent spatial variable $x \in [x_0,1]$ are clustered close to $x_0$ , where the solution changes faster. Within each element, the primary unknowns, $h,\,g,\,\varGamma,\,\varGamma _m$ , are expanded by quadratic Lagrangian basis functions, $\phi _i$ , and the evolution equations are weighted integrally with the same basis functions. The final weak forms are as follows, after applying integration by parts to reduce the order of differentiation:

(3.19) \begin{equation} \int _{x_0}^{1} \frac {\partial h}{\partial t}\phi _i \,{\rm d}x + \int _{x_0}^{1} \frac {2h\dot {R}}{R}\phi _i \,{\rm d}x + \frac {1}{R^2}\left [Q\phi _i\right]_{x_0}^{1} -\frac {1}{R^2}\int _{x_0}^{1} Q\frac {d\phi _i}{{\rm d}x}\,{\rm d}x = 0,\end{equation}

(3.20) \begin{equation} \int _{x_0}^{1} g\,\phi _i \,{\rm d}x - \left [\big(1-x^2 \big)\frac {\partial h}{\partial x}\phi _i\right]_{x_0}^{1} + \int _{x_0}^{1} \big(1-x^2 \big)\frac {\partial h}{\partial x}\frac {d\phi _i}{{\rm d}x}\,{\rm d}x = 0,\end{equation}

(3.21) \begin{equation} \int _{x_0}^{1} \frac {\partial \varGamma }{\partial t}\phi _i \,{\rm d}x + \int _{x_0}^{1} \frac {2\varGamma \dot {R}}{R}\phi _i \,{\rm d}x + \frac {1}{R^2}\left [Q_{\varGamma }\phi _i\right]_{x_0}^{1} -\frac {1}{R^2}\int _{x_0}^{1} Q_{\varGamma }\frac {d\phi _i}{{\rm d}x}\,{\rm d}x-\int _{x_0}^{1}j_b\phi _i \,{\rm d}x= 0,\end{equation}

and

(3.22) \begin{equation} \begin{aligned} \int _{x_0}^{1} \frac {\partial \varGamma _m}{\partial t}\phi _i \,{\rm d}x + \int _{x_0}^{1} \frac {2\varGamma _m\dot {R}}{R}\phi _i \,{\rm d}x + \frac {1}{R^2}\left [Q_{\varGamma m}\phi _i\right]_{x_0}^{1} -\frac {1}{R^2}\int _{x_0}^{1} Q_{\varGamma m}\frac {d\phi _i}{{\rm d}x}\,{\rm d}x \\ -\int _{x_0}^{1}(j_b+j_r)\phi _i \,{\rm d}x = 0.\end{aligned} \end{equation}

The fluxes in the integrated terms in (3.19)–(3.22) are by symmetry equal to zero at $x=1$ , and are evaluated at $x=x_0$ from the boundary conditions, (3.9)–(3.11). Time derivatives are treated with an implicit Euler scheme and the spatial integrals are evaluated by three-point Gaussian quadrature. Finally, the resulting system of nonlinear algebraic equations is solved with a Newton–Raphson iterative scheme.