## Abstract

We use a fully coupled, three-dimensional, finite-element method to study the evolution of the surface-tension-driven instabilities of a liquid layer that lines an elastic tube, a simple model for pulmonary airway closure. The equations of large-displacement shell theory are used to describe the deformations of the tube and are coupled to the Navier–Stokes equations, describing the motion of the liquid.

The liquid layer is susceptible to a capillary instability, whereby an initially uniform layer can develop a series of axisymmetric peaks and troughs, analogous to the classical instability that causes liquid jets to break up into droplets. For sufficiently high values of the liquid's surface tension, relative to the bending stiffness of the tube, the additional compressive load induced by the development of the axisymmetric instability can induce non-axisymmetric buckling of the tube wall. Once the tube has buckled, a strong destabilizing feedback between the fluid and solid mechanics leads to an extremely rapid further collapse and occlusion of the gas-conveying core of the tube by the liquid. We find that such occlusion is possible even when the volume of the liquid is too small to form an occluding liquid bridge in the axisymmetric tube.

## 1. Introduction

In healthy lungs, the airway walls are lined by a thin layer of fluid, which ensures that the pulmonary epithelial cells remain hydrated; in addition, the liquid lining forms part of the mechanism for the removal of foreign particles. Airway closure occurs when the gas-conveying core of the airway, the lumen, becomes blocked by the lung-lining fluid. In such cases, gas is ‘trapped’ behind the blockage, limiting the volume of gas that can be expired and leading to hysteresis in the pressure–volume curves of air-filled lungs. Airway closure tends to occur in the smaller airways towards the end of expiration, and has been demonstrated both experimentally (Burger & Macklem 1968; Engel *et al*. 1975) and histologically (Hughes *et al*. 1970).

Under normal physiological conditions, closed airways reopen during the early stages of inspiration and the rupture of the liquid bridges can generate lung sounds, e.g. pops and crackles (Suki *et al*. 1994), but does not affect gas exchange. Under pathological conditions, however, closure can occur in the larger airways and the vessels often remain closed for substantial portions of the breathing cycle, impairing the exchange of gases and reducing the efficiency of the lungs. Pulmonary pathologies that render the airways more susceptible to closure include respiratory distress syndrome, in which the surface tension at the air–liquid interface is elevated; emphysema, in which the tissue supporting the airway walls is destroyed, so that the walls are, in effect, more deformable; and pulmonary oedema, where the volume of the lung-lining fluid is increased. The details of how such changes in the physical properties of the lungs affect closure remain unclear; in part owing to the complex nature of the system dynamics, which are largely determined by the strong interaction between the elastic airway wall and the viscous lung-lining fluid.

In mechanical terms, airway closure is believed to be initiated by an axisymmetric, surface-tension-driven instability of the liquid film, analogous to the classical capillary instability that causes liquid jets to break up into droplets (Plateau 1873; Rayleigh 1879). In the undeformed airway, an approximately circular, cylindrical tube is lined by a uniform film of viscous fluid (liquid) surrounding a less viscous fluid (air) in the core. This configuration is unstable to axisymmetric, finite-wavelength perturbations when the length of the tube is greater than the inner circumference of the liquid film. The uniform film evolves into a series of periodically spaced, finite-amplitude lobes connected by thin films. If there is sufficient fluid present, the lobes can grow large enough to occlude the tube completely by forming axisymmetric liquid bridges (Gauglitz & Radke 1988; Kamm & Schroter 1989; Johnson *et al*. 1991).

The initial development of the instability is independent of the properties of the airway wall, but changes in fluid pressure during the evolution will alter the mechanical load on the wall, which can result in substantial wall deformation. Halpern & Grotberg (1992) studied the effects of wall elasticity on the instability and found that the reduced fluid pressure in the growing lobes causes the airway wall to move inwards, reducing the volume of fluid required to occlude the airway, compared with the case of a rigid wall. Halpern & Grotberg's (1992) study considered only axisymmetric deformations, however, and they did not investigate the possibility of non-axisymmetric buckling instabilities, known to affect elastic tubes under large compressive loads (Flaherty *et al*. 1972; Yamaki 1984).

Heil (1999) found that elastic tubes occluded by axisymmetric liquid bridges are indeed unstable to non-axisymmetric perturbations, provided that the surface tension of the lung-lining fluid is sufficiently large, relative to the bending stiffness of the tube. Furthermore, the statically stable, non-axisymmetric equilibrium configurations were found to require much smaller volumes of fluid than the corresponding axisymmetric liquid bridges. Nonetheless, the existence of such non-axisymmetric liquid bridges does not imply that they can ever be realized in a continuous evolution from the initial, axisymmetric, thin-film configuration. For example, airway closure at small fluid volumes would not be possible if the compressive load required to induce buckling were only reached once an axisymmetric liquid bridge has formed.

The continuous evolution of the system from the initial axisymmetric configuration, allowing for the possibility of non-axisymmetric deformation, was studied by White & Heil (2005), who used a modified lubrication-theory model to describe the behaviour of the fluid. In particular, White & Heil (2005) performed a linear stability analysis of the evolving (primary) axisymmetric instability to non-axisymmetric perturbations, and found that, for sufficiently high values of the surface tension, the system was linearly unstable to non-axisymmetric perturbations before the formation of any liquid bridges. If the initial transmural pressure is large enough, the tube buckles uniformly in a ‘ring’ mode, and the most unstable azimuthal mode is two-lobed (Heil & White 2002). If, on the other hand, the initial transmural pressure is too small to cause uniform buckling, the most unstable azimuthal mode is three-lobed. White & Heil (2005) also presented numerical simulations of the nonlinear evolution of the system.

Although a lubrication-theory description of the fluid mechanics is appropriate in the initial phases of the instability, it becomes invalid during the latter stages of the evolution, when the underlying assumptions are violated: the film becomes very thick and the interface becomes highly curved. In this paper, we complete the investigation of White & Heil (2005) by examining the time-evolution of the three-dimensional system using the three-dimensional Navier–Stokes equations to describe the behaviour of the lung-lining liquid.

## 2. The model

We model the airway as an elastic tube of undeformed radius *R*, length *L**, wall thickness *h*, Young's modulus *E* and Poisson's ratio *ν*. The tube is lined by an incompressible, Newtonian fluid of density *ρ*, viscosity *μ* that surrounds an air core, and the surface tension at the air–liquid interface is assumed to be a constant, *σ**. In reality, the airway epithelial cells secrete a natural surfactant, which dynamically alters the surface tension at the air–liquid interface, but this effect is neglected in the present model.

The problem is formulated in Cartesian coordinates, , where the *x*_{3}-direction is chosen to coincide with the axis of the tube and *x*_{1} and *x*_{2} are the transverse coordinates, see figure 1.

In the initial configuration, the liquid forms a film of uniform thickness *H*=*H**/*R*; throughout this paper asterisks are used to distinguish dimensional quantities from their dimensionless equivalents. The tube is subject to an external pressure *P**^{(ext)} relative to the air pressure in the lumen, which we choose to be the reference pressure and set to zero.

### (a) Wall equations

We use geometrically nonlinear Kirchhoff–Love shell theory to model the deformation of the elastic tube (e.g. Green & Zerna 1954; Wempner 1973). The assumptions underlying the theory are that (i) the dimensionless wall thickness *h*/*R* is small; (ii) material lines normal to the undeformed midplane remain unstretched and normal to the midplane during the deformation; and (iii) the ratio of the wall thickness, *h*, to the minimum radius of curvature of the deformed shell remains small. Assumption (ii) implies that the deformation of the shell may be completely specified by the displacements of its midplane, * v*=

**/*

**v***R*. The position of the undeformed midplane is given by , where , and are dimensionless Lagrangian coordinates; and the deformed midplane position is thus(2.1)Throughout this paper, Greek indices take the values

*α*=1, 2, and Latin indices take the values

*i*=1, 2, 3. We also use the convention that repeated indices are summed over all possible values of the index.

In order to characterize the geometry of the shell's undeformed midplane, we determine the base vectors **a**_{α}=**r**_{,α}, where the comma denotes partial differentiation with respect to *ζ*^{α}, and construct the covariant midplane metric tensor , with determinant *a*=*a*_{11}*a*_{22}−*a*_{12}*a*_{21}=1. We also define a curvature tensor, where is the (inward) unit normal to the midplane.

Uppercase letters are used to denote shell variables associated with the deformed midplane and we define the deformed midplane base vectors, **A**_{α}=**R**_{,α}; deformed covariant midplane metric tensor, , with determinant *A*; and deformed curvature tensor, , where is the (inward) unit normal to the deformed midplane.

The finite deformations of a general elastic body of undeformed volume, , loaded by tractions, * t*, acting on the deformed surface, , may be deduced from the principle of virtual displacements. In physical terms, the principle states that the external work during a virtual displacement must equal the internal work stored in the deformed body by the corresponding virtual change in its strain energy, and can be written as(2.2)where

*σ*

^{ij}is the second Piola–Kirchhoff stress tensor,

*ϵ*

_{ij}is the Green–Lagrange strain tensor and

*is the displacement field of the body. The symbol*

**v***δ*is used to indicate variations in a quantity.

The large bending deformations of the thin-walled elastic tube that occur in the present system only generate small strains, allowing us to employ a linear constitutive equation (Hooke's law). We determine the governing equations for the shell from the general expression (2.2) by assuming a stress-free initial state, plane stress in the midplane and integrating through the thickness of the shell analytically, which gives(2.3)Here, * f*=

**/*

**f***K*is the traction per unit area of the deformed midplane, non-dimensionalized by the bending modulus of the shell, ; and are the strain and bending tensors, respectively, which characterize the deformation of the shell's midplane and may be related to the Green–Lagrange strain tensor; andis the plane stress stiffness tensor, non-dimensionalized by Young's modulus, where

*a*

^{αβ}is the contravariant metric tensor of the undeformed midplane.

*A*is the determinant of the metric tensor of the deformed midplane, see above; its presence in the load terms reflects the fact that the load acts on the

*deformed*surface.

We take the variations in equation (2.3) with respect to the displacements of the midplane and their derivatives, and write the displacement vector in Cartesian components. The domain is decomposed into *N*_{S} two-dimensional, isoparametric, Hermite finite elements, in which the nodal displacements and gradients are independent degrees of freedom (Bogner *et al*. 1967). Equation (2.3) becomes(2.4)Here, *ψ*^{(S)} are the two-dimensional Hermite shape functions that interpolate the displacements and their derivatives at the nodal points.

### (b) Fluid equations

The initial motion of the fluid is driven by a balance between surface tension and viscosity and so we non-dimensionalize the fluid velocity and pressure by * u*=

**/(*

**u***σ**/

*μ*) and

*p*=

*p**/(

*σ**/

*R*), respectively. The dimensionless position vector is

*=*

**x****/*

**x***R*and

*t*=

*t**/(

*μR*/

*σ**) is the natural time-scale for the system. The dynamics of the liquid lining are governed by the dimensionless, unsteady Navier–Stokes equations(2.5a)and conservation of mass(2.5b)

In our non-dimensionalization, the Reynolds number, *Re*=*ρσ***R*/*μ*^{2} is a material constant.

We decompose the fluid domain into *N*_{F} (three-dimensional) isoparametric, Taylor–Hood finite elements (Taylor & Hood 1973), in which the pressure is approximated by tri-linear polynomials, *ψ*^{(P)}, and the global coordinates and velocity components by tri-quadratic polynomials, *ψ*^{(F)}. The free-surface is parametrized by the distance from the elastic tube, where *ζ*^{α} are the Lagrangian coordinates that parametrize the tube's midplane. The distance is measured in the direction * S*(

*ζ*

^{1},

*ζ*

^{2}): a continuous field of pre-determined unit vectors directed from fixed points

*(*

**R***ζ*

^{1},

*ζ*

^{2}) on the deformed tube wall towards the tube axis (0, 0,

*R*

_{3}(

*ζ*

^{1},

*ζ*

^{2})). The location of the free surface is, therefore,(2.6)and we use isoparametric, bi-quadratic elements to interpolate both the free-surface height and the Lagrangian coordinates

*ζ*

^{α}; a generalization of Kistler & Scriven's (1983) method of spines.

The momentum equations (2.5*a*) are weighted by the tri-quadratic velocity basis functions, *ψ*^{(F)}, and integrated over the fluid domain, *V*, to give the weak form(2.7a)Here, **u**^{(m)} is the velocity of the moving mesh, found by taking the time-derivative of the nodal positions; *S*_{fs} is the free surface; *S* is the union of all other surfaces bounding the fluid domain, namely, the tube wall and regions in the planes *x*_{3}=0, *L*; * n* is the unit normal to the bounding surface (directed out of the fluid); and

*C*is the line bounding the free surface, parametrized by the coordinate

*s*; see figure 1.

In equation (2.7*a*), the dynamic boundary condition at the free surfacehas been used to replace the integral of normal stresses over *S*_{fs} by an integral involving the surface curvature, *κ*=*κ***R* (Ruschak 1980). The surface-curvature integral is further modified by using the surface divergence theorem, Weatherburn (1955, p. 240, equation 26), which leads to a representation in terms of the surface coordinates *ζ*^{α}. The covariant base vectors of the free surface, the metric tensor and its determinant are given byThe notation [**g**_{β}]_{i} indicates the *i*th component of the base vector **g**_{β} in the global Cartesian basis. Finally, the unit vector * m* is tangent to the free surface and normal to the curve

*C*, and is directed out of the computational domain.

The weak form of the continuity equation follows from weighting equation (2.5*b*) by the tri-linear pressure basis functions, *ψ*^{(P)} and integrating over *V*(2.7b)

The system is completed by the kinematic boundary condition at the free surfacewhich is weighted by the bi-quadratic basis functions for the free-surface height, *ψ*^{(H)}, and integrated over the free surface(2.7c)providing a set of equations that determine the unknown film thickness .

### (c) Boundary conditions

#### (i) Fluid–solid coupling

The motion of the elastic tube is coupled to the fluid motion through the no-slip boundary condition(2.8)which requires that the fluid velocity at the wall is equal to the wall velocity. The fluid exerts a traction on the tube wall, and the load terms in equation (2.4) are given by(2.9)where *σ*=*σ**/*RK* represents the dimensionless surface tension and *P*^{(ext)}=*P**^{(ext)}/*K* is the dimensionless external pressure.

#### (ii) Symmetry conditions

We restrict the tube length, *L*=*L**/*R* to the range *π*(1−*H*)<*L*<2*π*(1−*H*) and assume symmetry in the axial direction at the ends of the tube. The restriction ensures that the system is unstable to the primary axisymmetric instability, but that only one lobe will form.

In order to examine the stability of the system to different axial wavenumbers, we add a small perturbation of the form *P*_{c} cos(*M* *ζ*^{2}) to the external pressure *P*^{(ext)} and impose azimuthal symmetry in the planes *ζ*^{2}=0 and *π*/*M*.

Exploiting the symmetry of the system restricts the computational domain to *ζ*^{1}∈[0,*L*/2] and *ζ*^{2}∈[0,*π*/*M*]. At the symmetry planes, the boundary conditions for the solid are(2.10a)where *x*_{n} is the coordinate normal to the symmetry plane; *v*_{n} is the displacement component normal to the symmetry plane and , are displacement components tangential to the plane, and normal to each other. The boundary conditions for the fluid are(2.10b)where *u*_{n} is the fluid velocity in the direction normal to the symmetry boundary and , are the fluid tractions tangential to the symmetry plane, and normal to one another.

If the normal to the symmetry plane is one of the global Cartesian coordinate directions, the implementation of equations (2.10*a*,*b*) is straightforward. For example, imposing axial symmetry at the ends of the tube yields essential boundary conditions for the solid displacement and derivatives(2.11a)where *v*_{i} is the component of the displacement in the *x*_{i}th coordinate direction. There is one essential boundary condition for the fluid velocity(2.11b)and transverse fluid tractions are set to zero by omitting the surface integrals over the planes *x*_{3}=0, *L* in equation (2.7*a*)—the natural boundary conditions.

Imposing symmetry at *ζ*^{2}=0 (in the plane *x*_{2}=0) is similar, we omit the surface integral over *x*_{2}=0 in equation (2.7*a*) and apply the essential conditions(2.12a)(2.12b)If *M*≠2, then the normal to the plane *ζ*^{2}=*π*/*M* does not coincide with a coordinate direction and the implementation of the boundary conditions (2.10*a*,*b*) must be altered. Resolving the normal and tangential directions into the global Cartesian coordinate system, the solid boundary conditions become(2.13a)Similarly, the boundary condition on the fluid velocity normal to the symmetry plane is(2.13b)

The boundary conditions on the fluid tractions imply that the traction has a component only in the normal direction, *t*_{i}=*t*_{k}*n*_{k}*n*_{i}, which implies that(2.13c)and this expression is used in the surface integral over the symmetry plane in equation (2.7*a*). The equations (2.13*a*,*b*) are treated as additional constraints and are imposed at all nodes on the symmetry plane via Lagrange multipliers.

### (d) Numerical implementation

The problem involves a number of disparate time-scales and, in the interests of efficiency, Gresho & Sani's (2000, pp. 805–6) adaptive, ‘predictor–corrector’, linear, multistep method was used. The predictor is a second-order, variable-step (explicit) leapfrog scheme and the corrector is a second-order, variable-step (implicit) BDF2 scheme. We denote the time after the *n*th time-step by *t*_{n+1}=*t*_{n}+Δ*t*_{n}, where Δ*t*_{n} need not remain constant. The relative change in nodal positions between time-steps is used to control the error and *d*_{ijn} is defined to be the local truncation error in the *x*_{i}-direction at node *j* after the *n*th time-step(2.14)where *R*_{ijn} is the *i*th component of the position vector to node *j* after the *n*th time-step and is the corresponding predicted value. We use a root mean square error norm(2.15)where is the total number of nodes. The size of the next time-step is then chosen according to the formula(2.16)where *ϵ* is an error control parameter, usually 10^{−4}≤*ϵ*≤10^{−8}. Unless stated otherwise, runs were conducted with *ϵ*=10^{−6}.

At *t*=0, the fluid velocities are set to zero and an initial perturbation of the form *H*(1+0.01 cos(*x*_{3}*π*/*L*)) is applied to the film height in order to initiate the Rayleigh instability.1 The first time-step was taken using a second-order trapezoid (Crank–Nicolson) scheme. The next three time-steps were taken using the BDF2 scheme, but with a constant Δ*t*. All subsequent time-steps were allowed to vary according to the scheme described above. If the scale factor *S*≡(*ϵ*/‖*d*_{n}‖)^{3}<0.8, the predicted value is sufficiently different from the true solution that we reject the solution and repeat the step with a smaller Δ*t*. If *S*>4, we set *S*=4, so that Δ*t* is not increased by a factor of more than four between any two consecutive steps.

The coupled system of nonlinear, algebraic equations that arises at each time-step was solved numerically using a Newton–Raphson method, and a frontal solver (Duff & Scott 1996) was used to assemble and invert the Jacobian matrices. The entries in the Jacobian matrices that arise from the fluid equations were calculated analytically, but all others were determined by finite differencing the governing equations.

The fluid–structure interaction solver has been previously validated in a number of steady problems (Hazel & Heil 2003*a*,*b*). The time-dependent portions of the code were validated by comparison with the lubrication-theory-based computations of White & Heil (2005), and the initial growth rates of the axisymmetric instability. In order to validate the application of the symmetry boundary conditions in planes that are not coordinate axes, contact pressures for the axially uniform deformation of an elastic tube (ring) were compared with those of Flaherty *et al*. (1972) for axial wavenumbers of two, three and four. In addition, it was confirmed that the evolution of the primary axisymmetric instability remained the same for the fluid domains of different azimuthal extents: *π*/2, *π*/3 and *π*/4. Final code validation was obtained by repeating selected calculations with finer spatial and temporal discretizations; see figures 3 and 10 below.

## 3. Results

### (a) Parameter values

In all cases, we used a Poisson's ratio of *ν*=0.49 and set the wall thickness to *h*/*R*=1/20; in reality the wall thickness of the distal airways is slightly higher, but this value is close to the upper limit for which thin-shell theory is applicable.

The initial load on the elastic tube is the combined load due to the external pressure *P*^{(ext)} and the internal load caused by the jump in pressure across the air–liquid interface, *σ*/(1−*H*). Following White & Heil (2005), we shall refer to this combined load as the initial pressure(3.1)

Heil & White's (2002) linear stability analysis for the case of axially uniform (ring-mode) tube deformations demonstrated that the tube becomes unstable to non-axisymmetric perturbations when *P*^{(init)}>3. In the present context, we wish to examine cases where the system is stable to such perturbations and so we set *P*^{(init)}=2.9. Throughout the parameter variations, *P*^{(init)} will remain fixed so that we may determine the effects of variations in *σ* without the (most obvious) effect of increasing the initial load on the tube. Physically, keeping *P*^{(init)} constant means that any increase in surface tension is balanced by a compensatory decrease in external pressure.2

Unless otherwise stated, we set the length of the tube to be the half-length of the fastest-growing wavelength of the initial axisymmetric instability, (Goren 1962; Hammond 1983; Halpern & Grotberg 1992). The effects of variations in *L* are considered in §3*c*.

We found, as in the two-dimensional problem (Heil & White 2002), that fluid inertia has a negligible effect on the results (see §3*e* for further details) and we shall present results only for the case *Re*=0. Finally, we observed that the time-scale *μR*/*σ** used in the non-dimensionalization of the governing equations does not, in fact, reflect the time-scale of the surface-tension-driven collapse. Following White & Heil (2005), we present results on the more appropriate time-scale .

### (b) Variations in *σ*

Figure 2 illustrates the surface-tension-driven, non-axisymmetric collapse of a liquid-lined elastic tube by plotting the instantaneous fluid volumes at four different times during the system's evolution for an initial film thickness of *H*=0.1. The surface tension is *σ*=70>*σ*_{(min)}≈63, where *σ*_{min} represents White & Heil's (2005) theoretical prediction for the minimum value of the surface tension at which the system will become unstable to non-axisymmetric perturbations. The initial, uniform-film configuration shown in figure 2*a* is susceptible to the axisymmetric Rayleigh–Plateau instability, which causes the fluid to drain from one end of the tube to the other, and a sizeable axisymmetric fluid lobe has developed by *t*=5.169 (figure 2*b*). The fluid pressure decreases in the developing lobe, leading to a locally increased compression on the elastic wall. Eventually, the compressive load becomes large enough to cause the tube to buckle non-axisymmetrically and a notable change in the tube-wall shape is observed by *t*=7.860 (figure 2*c*). The film thickness in the buckling region is large and the moderate azimuthal pressure gradients induced by weak non-axisymmetries of the air–liquid interface are sufficient to redistribute the fluid such that the air–liquid interface remains approximately axisymmetric throughout the evolution. Once the film thickness has become large enough, the non-axisymmetric collapse accelerates, and the radius of the air–liquid interface tends to zero very rapidly; the tube is completely occluded by the liquid at *t*=7.861.

A rapid decrease in the radius of the air–liquid interface was also observed in the axisymmetric computations of Halpern & Grotberg (1992), the lubrication-theory model of White & Heil (2005), and the two-dimensional calculations of Heil & White (2002). It is a consequence of a strong destabilizing feedback between the fluid and solid mechanics: when the tube buckles, the radius of the air–liquid interface decreases, causing an increase in the curvature of the interface and, accordingly, an increase in the pressure drop across the interface. The increased pressure drop results in a lower fluid pressure and hence a higher compressive load on the tube wall, which causes the tube to collapse even further.

We now describe the temporal evolution of the system by plotting the radius of the air–liquid interface at the four positions shown in the inset in figure 3. The solid lines in figure 3 show the evolution of the four radii for *σ*=70. At the initial time, *t*=0, the film is uniform and the air–liquid interface has a constant radius of 0.9. The development of the axisymmetric lobe is indicated by the decrease in radii *R*_{1} and *R*_{2}, and the corresponding increase in the radii *R*_{3} and *R*_{4} at the other end of the tube where the film thins. Following the initial exponential growth (for *t*≲0.75), the axisymmetric instability saturates and the system continues to evolve very slowly towards an axisymmetric equilibrium state in which all the fluid is contained in the lobe; see Gauglitz & Radke (1988) or White & Heil (2005) for a more detailed description of this phase of the evolution. As fluid continues to drain into the main lobe, the curvature of the air–liquid interface increases, leading to a slow increase in the compressive load on the tube wall. Eventually, the compression becomes large enough to cause the tube to buckle non-axisymmetrically, and ultimately, this leads to an extremely rapid decrease in radius in the lobe region, at *t*≈7.86.

For values of *σ* slightly above *σ*_{min}, as in the case considered here, the axisymmetric lobe must be almost fully developed before the compressive load is sufficient to cause the tube to buckle non-axisymmetrically. In this regime, the time at which the ultimate non-axisymmetric collapse occurs is extremely sensitive to small changes in the governing parameters. To illustrate this, the dashed line in figure 3 shows the system's evolution at a slightly larger surface tension, *σ*=71. Before the onset of the non-axisymmetric instability, the system's evolution differs very little from that when *σ*=70; for instance, at *t*=6.3, just before the non-axisymmetric collapse for *σ*=71, the radii of the air–liquid interfaces differ by less than 0.3%. However, the extremely slow evolution of the axisymmetric state causes the system with *σ*=71 to reach the critical compressive load much earlier than for *σ*=70.

The extreme sensitivity of the system in this regime is particularly challenging for any numerical time-stepping scheme. The results for *σ*=70 using a smaller error-control parameter, *ϵ*=10^{−8} are shown as square markers in figure 3 and are in good agreement with the case when *ϵ*=10^{−6}.

In previous lubrication-theory-based computations of the system, it was not possible to follow the evolution much beyond the onset of the final catastrophic collapse. In contrast, solving the Navier–Stokes equations permits us to follow the evolution all the way to the point of closure, at which the radius of the air–liquid interface reaches zero and the airway is completely occluded by the fluid.

Figure 4 shows the temporal evolution of the four control radii for three different values of *σ*. For *σ*=50<*σ*_{min}, the system is unstable to the axisymmetric instability, but the compressible load never becomes large enough to cause non-axisymmetric buckling and the tube remains open. At *σ*=70, just above the critical value, there is a relatively long period in which the primary instability appears to saturate, followed by a rapid collapse, as discussed above. At *σ*=100, the surface tension is large enough to provoke catastrophic non-axisymmetric collapse of the tube before the axisymmetric instability reaches saturation. The diamond-shaped markers in figure 4 indicate the times at which the lubrication-theory-based, linear stability analysis predicts the system to become unstable to non-axisymmetric perturbations.

### (c) Variations in *L*

If the tube length differs from that of the half-length of the fastest-growing wavelength of the axisymmetric instability, , the initial (linear) growth-rate of the primary instability is, of course, reduced. The subsequent nonlinear development of the system depends on whether the tube length is greater or smaller than . Figure 5 shows the evolution of the four control radii for three different tube lengths, , and 1.95*π*(1−*H*), when *H*=0.1 and *σ*=70.

When , a greater volume of liquid is available for recruitment into the developing lobe, and the liquid film remains thicker throughout, permitting a faster axial redistribution of the fluid. Hence, the nonlinear development proceeds more rapidly than when , and the ultimate catastrophic collapse occurs at significantly earlier times. The axial extent of the collapsed region appears to remain localized and constant, however, taking place over approximately two tube diameters, see figure 6.

When , less fluid is available for recruitment into the lobe, and thin films develop that retard the axial redistribution of the liquid. The nonlinear development is, therefore, always slower than that when , delaying the ultimate collapse of the tube. For sufficiently short tubes, the lack of fluid completely suppresses the collapse; the radius of curvature of the air–liquid interface in the lobe remains relatively high, even when the lobe is fully developed, and the pressure jump across the interface is not large enough to induce buckling of the tube.

### (d) Variations in *H*

White & Heil's (2005) linear stability analysis predicted that variations in the film thickness, *H*, would have two main effects: (i) a reduction in *H* reduces the growth rate of the primary axisymmetric instability; (ii) at smaller film thicknesses, there is an increase in the minimum surface tension, *σ*_{min}, required to make the system unstable to non-axisymmetric perturbations. The second effect is a consequence of the smaller curvature of the developing lobe and, hence, smaller pressure jump across the air–liquid interface. Both these predictions are confirmed by our Navier–Stokes simulations.

In the context of the pulmonary airway closure problem, one of the key parameters is the minimum volume of fluid required to form occluding liquid bridges. In an attempt to determine this value, Gauglitz & Radke (1988) performed numerical simulations of the axisymmetric Rayleigh–Plateau instability in rigid tubes. The simulations were started from a uniform-film configuration, perturbed by the fastest-growing unstable mode. A film thickness of at least *H*≈0.12 was required to form an occluding liquid bridge during the long (but finite) period of the simulations. It is possible that thinner films could evolve towards occluding liquid bridges over even longer periods of time. The occlusion of a rigid tube is impossible, however, if the volume of fluid contained in the liquid lining is less than the volume required to form a ‘minimal liquid bridge’—a liquid bridge in which the opposite air–liquid interfaces come into contact, as sketched in figure 7. Assuming that the liquid bridge develops from the fastest growing perturbation to the axially uniform initial state, conservation of fluid volume implies that airway closure in rigid tubes is impossible for ; see White & Heil (2005).

In order to assess whether non-axisymmetric instabilities allow the formation of occluding liquid bridges at smaller fluid volumes, we performed simulations with a film thickness of *H*=0.0395 and a surface tension of *σ*=320, which is larger than White's (2003) prediction of *σ*_{min}≈290. Figure 8 illustrates the system's evolution with plots of four instantaneous fluid volumes, while figure 9 shows the evolution of the four control radii as functions of time.

The initial dynamics of the system are the same as for the thicker film: the axisymmetric instability develops and saturates. Eventually, the compressive loading in the lobe region is high enough to cause notable non-axisymmetric buckling (*t*=0.843; figure 8*a*). In this case, however, the film in the lobe is significantly thinner than when *H*=0.1 and the elevated flow resistance weakens the transverse flows that try to return the air–liquid interface to an axisymmetric shape. This is illustrated in the inset in figure 9, which shows that the radii *R*_{1} and *R*_{2} differ over the range 0.846<*t*<0.854. The strongly elevated fluid pressures in the thin-film regions are strong enough to cause an outward buckling of the tube in these regions, as shown in figure 8*c*,*d*. The tube shapes in the final stages of the evolution are, therefore, quite different to the case when *H*=0.1. The development of the transverse draining flows and high fluid pressures oppose the collapse, but cannot inhibit the final catastrophic collapse towards an occluded state. In this case, the final collapse becomes so rapid that the adaptive time-stepping scheme has to take excessively small time-steps, which forced us to halt the computations. We believe that these were the most technically challenging computations, and so we used them as the basis for validation of the spatial discretization, see appendix A.

### (e) The effect of fluid inertia

The inclusion of fluid inertia did not significantly affect the results, even though the nominal Reynolds number, *Re*=*ρσ***R*/*μ*^{2}, used in the non-dimensionalization of the Navier–Stokes equations can be very large: e.g. *Re*=5000 in the terminal bronchioli, based on Halpern & Grotberg's (1992) estimates for the physical parameter values, *R*=250 μm, *σ**=20 dynes cm^{−1}, *ρ*_{w}=*ρ*=1000 kg m^{−3}, *μ*=10^{−3} kg (m s)^{−1}. This is because the fluid is initially at rest and the non-dimensional fluid velocities (on the *σ**/*μ*-scale) remain very small until the final stages of the collapse. During these final stages, the presence of fluid inertia was found to ‘slow down’ the collapse very slightly, but did not otherwise alter the system's dynamics.

## 4. Summary and discussion

This paper complements and completes White & Heil's (2005) investigation of three-dimensional instabilities of surface-tension-driven non-axisymmetric instabilities of the pulmonary airways. Using the Navier–Stokes equations to model the dynamics of the lung-lining fluid has allowed us to confirm that for values of *σ* greater than the critical *σ*_{min} predicted by White & Heil's (2005) linear stability theory, the system undergoes a non-axisymmetric collapse that eventually leads to complete occlusion of the airway by the fluid. In particular, we found that non-axisymmetric instabilities allow airway closure to occur in cases where the initial film thickness is too small for axisymmetric liquid bridges to form. The parameter values used in our simulations are in a range that is appropriate for the conditions in the small airways. For instance, Halpern & Grotberg's (1992) estimates for the physical parameters in the terminal bronchioli correspond to a non-dimensional surface tension of *σ*=120. Furthermore, for a non-dimensional film thickness of *H*=0.1, the time-scale *μ*/(*KH* ^{3})=1.5 s, indicating that the instabilities develop over a time-scale that is comparable to the duration of the breathing cycle (see White & Heil 2005).

Our computations were performed with periodic boundary conditions, but we found that the non-axisymmetric collapse remains localized over approximately two tube diameters, provided that the length of the tube is greater than the fastest-growing axial wavelength of the axisymmetric instability, . Therefore, we do not expect the additional structural support that may be provided by airway bifurcations in the lungs to significantly affect our results.

In the interests of simplicity, the present study has also ignored a number of other potentially important physical effects. For instance, we have not represented the effects of pulmonary surfactant in our model. In previous studies (Halpern & Grotberg 1993; Otis *et al*. 1993; Cassidy *et al*. 1999) it was found that surfactant acts to delay airway collapse. In part, this effect is due to a reduction in mean surface tension. Additionally, surfactant is advected towards the developing fluid lobe, leading to a higher local surfactant concentration and, therefore, a lower surface tension. The ensuing surface-tension gradient between the fluid film and the developing lobe acts in the opposite direction to the pressure gradient and opposes the further development of the lobe. Therefore, we expect that the inclusion of pulmonary surfactant in the present model would serve to increase the time until airway collapse, but would not be able to suppress the collapse altogether.

Intermolecular forces, such as van der Waals forces, were also ignored though they could become important in regions of small film thickness where they might cause the film to rupture. Although intermolecular forces could be included in our model, their inclusion would not necessarily lead to a better description of the actual physics of airway closure. We believe that local surface irregularities of the airway wall, which were also neglected in our model, would begin to affect the behaviour of the fluid long before short-range molecular forces become important.

The use of a linear constitutive equation for the airway wall is only justified provided that the strain remains small. We monitored the components of the strain tensor, *γ*_{αβ}, during the evolution of the system. The greatest strains were induced by circumferential compression of the shell in the buckling region. In all cases, *γ*_{22} remained less than 5%, with values greater than 1% attained only during the final stages of collapse.

The model of the airway wall is necessarily simplified, neglecting the external tethering, the wall's multilayer structure, viscoelastic effects and any material nonlinearities. However, we do not believe that the inclusion of any of these effects would be able to suppress airway closure by the mechanism analysed in this paper. Once the radius of the air–liquid interface has become sufficiently small, the instability is driven by the strongly destabilizing feedback between fluid and solid mechanics, which is independent of the precise details of the wall mechanics. Moreover, the effects of variations in the non-dimensional surface tension, *σ*, the tube length, *L*, and the film thickness, *H*, on the system's dynamics are expected to remain the same. It is entirely likely, however, that the inclusion of external tethering and/or wall damping would require higher values of *σ* to induce collapse on the same time-scales as in our simplified model. Previous two-dimensional studies (Wang *et al*. 1983; Wiggs *et al*. 1997) have demonstrated that external tethering and inclusion of the multilayer structure of the airway wall increase the most unstable azimuthal buckling wavenumber, which is in accordance with physiological observations (Lambert *et al*. 1994), indicating that the wavenumber of collapsed airways is considerably higher than that predicted by our model.

## Acknowledgments

We gratefully acknowledge financial support from the EPSRC. The HSL library routine MA42 and the cfortran.h header file were used in the computations.

## Footnotes

As this paper exceeds the maximum length normally permitted, the authors have agreed to contribute to production costs.

↵In fact, such a perturbation is not required. After sufficiently long periods, the accumulation of random numerical round-off errors is sufficient to initiate the Rayleigh instability. Nevertheless, in order to provide meaningful comparisons, each parameter study is started with same initial perturbation.

↵We note that this process is similar to the therapeutic procedure known as PEEP (positive end expiratory pressure) ventilation. The procedure is used for the mechanical ventilation of patients whose lung liquid lining has an abnormally high surface tension, e.g. in respiratory distress syndrome. The internal pressure of the airways is elevated by a constant amount so that the mechanical load on the airways is reduced and the airways do not collapse.

- Received September 6, 2004.
- Accepted December 1, 2004.

- © 2005 The Royal Society