Stéphane Popinet

Stéphane Zaleski

A front tracking algorithm for the accurate representation of surface tension

A front-tracking algorithm for accurate representation of surface tension

Introduction

A number of methods have been developped in recent years for the solution of problems involving free interfaces. These methods can be divided in two main classes depending on the type of grids used. In the first class, the interface is treated as a boundary between elementary domains. This approach allows a precise discretization of the interfacial terms but requires deformable grids in order to follow the motion of the interface [START_REF] Fyfe | Surface tension and viscosity with lagrangian hydrodynamics on a triangular mesh[END_REF] [START_REF] Glimm | Front tracking applied to Rayleigh-Taylor instability[END_REF]. The second class of methods uses fixed grids to describe the velocity field but requires specific advection schemes in order to preserve the sharpness of the interfacial front. Of these, modern Volume Of Fluid (VOF) advection schemes give good results and ensure an accurate conservation of mass [START_REF] Ashgriz | Flair: Flux line-segment model for advection and interface reconstruction[END_REF][14] [START_REF] Puckett | A high order projection method for tracking fluid interfaces in variable density incompressible flows[END_REF]. Level set type methods are easy to implement but are probably less accurate in term of mass conservation [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF]. Both rely on an implicit description of the interface, given through phase functions (i.e. volume fraction for VOF or distance function for level set), and are usually called front capturing methods. Another group, front tracking algorithms use an explicit discretization of the interfacial discontinuity, and if somewhat more complex to implement, give the precise location and geometry of the interface [START_REF] Tryggvason | Computations of three-dimensional Rayleigh-Taylor instability[END_REF] [START_REF] Unverdi | A front-tracking method for viscous incompressible, multi-fluids flows[END_REF].

All these methods provide good solutions to the problem of interface advection, however accurate representation of surface forces (i.e. surface tension, membrane effects, . . . ) remains a problem when using fixed grids. In particular, most of these methods (including lattice gases [?]) exhibit the so called spurious currents shown in figure 1. These numerical artefacts result Figure 1: Spurious currents around a stationary bubble. The method used for the interface advection is a VOF type scheme. from inconsistent modelling of the surface tension terms and the associated pressure jump. More generally, this poses the problem of the accuracy of the description of the steep gradients occuring at the interface.

Many practical problems involve fluids of very different densities and/or viscosities (water droplets in air, wave breaking for example). Moreover, the surface tension effects are usually large compared to the viscous damping. In these cases, the effect of the unbalanced forces acting on the interface not only reduces the accuracy but can lead to spurious currents which create oscillations strong enough to destroy the interface. The spurious currents are therefore more than simply a numerical inaccuracy but are a real limitation of methods used on fixed grids.

In this paper, we present a front tracking method for the solution of the 2D incompressible Navier-Stokes equations with interfaces and surface effects. We will more particularly investigate the accuracy of the numerical representation of the surface tension and of the associated pressure jump. In the section which follows, we describe the general scheme used to solve the Navier-Stokes equations. The second section gives a detailed description of the front tracking algorithm and in the final section, some test cases investigating the numerical accuracy and convergence properties are presented.

2 General description of the method

Basic equations

We seek to solve the incompressible Navier-Stokes equations with varying density and surface tension. The momentum balance equations are

ρ(∂ t u + u • ∇u) = -∇p + ∇ • µD + σκδ s n + F (1)
∇ • u = 0. ( 2
)
These may be written in conservative form

∂ρu ∂t + ∇ • (ρu ⊗ u) = -∇p + ∇ • µD + σκδ s n + F, (3)
where u = (u, v) is the fluid velocity, ρ = ρ(x, t) is the fluid density, µ = µ(x, t) is the fluid viscosity, D is the viscous stress tensor. The surface tension term is considered to be a force concentrated at the interface, σ is the surface tension, κ the local curvature of the interface, δ s a function equal to one on the interface, zero outside and n the unit outward normal to the interface. Density and viscosity are advected by the velocity field

∂ρ ∂t + u • ∇ρ = 0, (4)
and

∂µ ∂t + u • ∇µ = 0. (5)
In the case of two immiscible fluids, we use the volume fraction C of one of the fluids which verifies the advection equation

∂C ∂t + u • ∇C = 0. ( 6
)
The volume averaged density and viscosity can then be expressed as functions of

C ρ = Cρ 1 + (1 -C)ρ 2 (7)
µ = Cµ 1 + (1 -C)µ 2 . ( 8
)

The Navier-Stokes solver

We use a momentum-conservative version of the Navier-Stokes solver initially developed for the SURFER code [START_REF] Lafaurie | Modelling merging and fragmentation in multiphase flows with surfer[END_REF]. The solution is based on a projection method [START_REF] Peyret | Computational Methods for Fluid Flow[END_REF]. The pressure, volume fraction, momentum and velocity com- [START_REF] Peyret | Computational Methods for Fluid Flow[END_REF] distribution (figure 2). If we associate the control volumes Ω x and Ω y of figure 3, to the component ρu and ρv of the momentum, we can write the integral formulation of (3) without the pressure term as

Ui,j ρ Vi,j+1 ρ Vi,j ρ Ui+1,j ρ Ci,j Pi,j
∂ ∂t Ω ρu + ∂Ω ρu ⊗ u = ∂Ω µD + Ω σκδ s n + ΩF. ( 9
)
The discrete formulation of ( 9) is then

∂ ∂t (ρu Ωij ) + faces l f ρu ⊗ u = faces l f µD + F ij Ωij + Ω ij σκδ s n. ( 10
)
where Ωij is the area of the control volume Ω ij , l f the length of a face and ρu ij is the average value of ρu in Ω ij . For each time step ∆t, the explicit projection method can be subdivided in four steps:

1. An approximate solution ρu ⋆ of the momentum at time (n + 1)∆t is built using an explicit discretization of (10)

1 ∆t (ρu ⋆ ij -ρu n ij ) Ωij = - faces l f ρu⊗u n + faces l f µD n +F n ij Ωij + n Ω ij σκδ s n.
(11) where Φ n is the value of the function Φ at time n∆t.

2. The volume fraction field C is advected using the velocity field at time n∆t. In our case this is done using a Lagrangian advection of the markers and a reconstruction of the volume fraction field.

3. We then need to ensure the divergence-free character of the velocity field u n+1 defined as ρu n+1 /ρ n+1 . The approximate solution u ⋆ = ρu ⋆ /ρ n+1 is corrected with the pressure gradient term

u n+1 = u ⋆ - ∆t ρ n+1 ∇p. (12)
In order to satisfy the continuity equation ( 2), the pressure must obey the Poisson-like equation

∇ • ( ∆t ρ n+1 ∇p) = ∇ • u ⋆ . ( 13
)
This equation is solved efficiently using a multigrid solver [START_REF] Brandt | Guide to multigrid development[END_REF][5][20] [START_REF] Wesseling | An introduction to Multigrid Methods[END_REF].

4. the momentum and velocity fields at time (n+1)∆t are computed using

ρu n+1 = ρu ⋆ + ∆t∇p (14)
and

u n+1 = ρu n+1 /ρ n+1 . (15)

Interface and surface tension

Given a parametric description (x(s), y(s)) of the location of the interface, where s is the curvilinear coordinate, we want to compute the source term due to the surface tension in the momentum equation [START_REF] Glimm | Front tracking applied to Rayleigh-Taylor instability[END_REF]. If Ω is the control volume for one of the components of the momentum and AB, the section of interface inside this control volume (figure 4), the integral source term in ( 9) is

Ω σκδ s n = σ B A κnds, (16)
which can be written using the first Frenet's formula for parametric curves

σ B A κnds = σ B A dt = σ(t B -t A ) ( 17
)
where t is the oriented unit tangent to the curve. The integral source term due to the surface tension is then the sum of the outward unit tangents at the points where the interface enters or exits the control volume.

A

Ω U t tB A B ρ
Figure 4: Method of tensions for computing the integral contribution of the surface tension to the momentum equation.

Surface tension, pressure jump and the pressure gradient

The gradient terms in the momentum equation are discretized using a classical second-order accurate centered scheme. This simple discretization is well suited for the representation of the gradient terms in the bulk of each phase. However, very rapid changes occur near the interface and this type of discretization no longer represents correctly the associated gradients.

As we are interested in phenomena including high surface tension effects, a term of particular concern is the pressure gradient. As first shown by Laplace, in the static case the surface tension is exactly balanced by the jump of the pressure through the interface. The surface tension then yields a pressure jump which can be expressed as

[p] σ = σκ, ( 18
)
where σ is the surface tension coefficient and κ, the curvature. If we look more closely at what happens on our MAC grid, with a simple discretization of the pressure, we can show that we cannot ensure the balance between the pressure gradient and the surface tension (this was first shown by Jie Li [15]). Let us consider the case of a stationary circular bubble. Figure 5 shows a piece of the circular interface and the corresponding pressure cells. The control volume for the ρv component of the momentum is crossed by the interface. As the outward tangents are not colinear, the non zero surface tension contribution to the momentum should be balanced by the pressure gradient ∂p/∂y ≈ h -1 (p1 -p4). However, as the interface does not cross the [p1, p2] segment, there is no pressure jump between p1 and p2 and p1 = p2 (where p1 is the constant pressure inside the bubble). Following the same argument, p2 = p3 = p4 and then p1 = p4. The discretized pressure gradient h -1 (p1-p4) is then zero and cannot balance the surface tension contribution. In practical terms this means that a stationary solution cannot exist using this simple discretization of the pressure gradient at the interface.

In order to avoid spurious currents, we must first find a better way to discretize the pressure gradient near the interface. We must note however that this will only ensure the existence of a stationary solution, with the hypothesis that the exact pressure field is known. The eventual convergence of the method to this stationary solution and its stability are another problem.

From figure 5, we can easily see that a better approximation of the pressure gradient can be obtained provided we use a better approximation of the pressure fluxes on the AC segment. Due to the pressure jump at the intersection point B, the pressure flux p AC should be decomposed in two parts p AB and p BC . In the stationary case, the integral contribution of the pressure gradient to the momentum equation would then be discretized as

- Ωy ∂p ∂y = AC p4 -(AB p AB + BC p BC ), (19)
which could be written using ( 18)

AC p4 -(AB p out + BC p in ) = AB (p in -p out ) = AB σκ, (20)
where p out and p in are the pressures outside and inside the bubble respectively. A parametric representation (x(s), y(s)) of the circular interface is

x(s) = x o + κ -1 cos(s) (21) y(s) = y o + κ -1 sin(s). ( 22
)
If we note s B and s D the values of the parameter for the intersection points B and D, we have the relation AB = κ -1 (cos(s B ) -cos(s D )) and then

- Ωy ∂p ∂y = σ(cos(s B ) -cos(s D )). ( 23
)
Using the tension formulation [START_REF] Pawar | Marangoni effects on drop deformation in an extensional flow: The role of surfactant physical chemistry. I. Insoluble surfactants[END_REF], the surface tension contribution to the y component of the momentum can be expressed as

σ D B κn y ds = σ(t y D -t y B ) = σ(κy ′ (s D ) -κy ′ (s B )) = σ(cos(s D ) -cos(s B )). ( 24
) Thus, in the case of a circular interface the exact balance between the pressure gradient term -Ωy ∂p/∂y and the surface tension Ωy σκn y δ s is satisfied provided we know the pressure field and the location of the intersection points B and D exactly.

It is important to note that this approach could be and should be applied to the discretization of all the gradient terms with rapid variations across the interface, in particular the viscous stress tensor for high contrasts of viscosities.

Front tracking algorithm

In this section, we present a detailed description of the front tracking algorithm. The interface is represented using an ordered list of points (x i , y i ), 1 ≤ i ≤ N. A list of connected polynomials (p x i (s), p y i (s)) is constructed using the points and gives a parametric representation of the interface, with s an approximation of the curvilinear coordinate. Both lists are oriented and thus identify the topology of the interface.

Advecting the points

The first step in our algorithm is the advection of the points dataset. We use a simple bilinear interpolation to find the velocity of a point in a given velocity cell (u i,j , u i+1,j , u i,j+1 , u i+1,j+1 )

u(x, y) = u i,j (1 -x -y + xy) + u i+1,j x(1 -y) + (25) u i,j+1 y(1 -x) + u i+1,j+1 xy,
where x and y are the coordinates of the point relative to the (i, j) vertex. The points are then lagrangianly advected using a straightforward first order explicit scheme ,y). Once the points have been advected we need to reconstruct the parametric representation of the interface.

x n+1 i = x n i + ∆t u(x, y) (26)
y n+1 i = y n i + ∆t v(x

Constructing the polynomials

We have chosen to use connected cubic polynomials with continuous first and second order derivatives. This type of curve is usually known as cubic splines [START_REF] Walsh | Best approximation properties of the spline fit[END_REF][11] [START_REF] Press | Numerical Recipes, The Art of Scientific Computing[END_REF]. The parametric representation must be periodic as the interfaces are mostly self-connected (drops, bubbles, periodic wave trains . . . ).

We then need to choose a parameter s in order to interpolate the two sets of points (x i , s i ) and (y i , s i ). A simple choice is an approximation of the curvilinear coordinate

s i = i-1 j=1 (x j+1 -x j ) 2 + (y j+1 -y j ) 2 . ( 27
)
The connection conditions for the interpolating polynomials lead to two pseudo tridiagonal systems Ba = c, one for each coordinate of the parametric curve, where B is a N 2 matrix of the form

B =          b 1 c 1 0 0 α a 2 b 2 c 2 0 0 0 . . . . . . . . . 0 0 0 a N -1 b N -1 c N -1 α 0 0 a N b N          (28)
where the α coefficients arise from the periodicity condition [START_REF] Hämmerlin | Numerical Mathematics[END_REF]. The solution of this type of system can be reduced to the solution of two tridiagonal systems which are easily solved using classical techniques [START_REF] Press | Numerical Recipes, The Art of Scientific Computing[END_REF].

The construction of the interpolating parametric spline curve from the (x i , y i ) set of points then requires the solution of four tridiagonal systems of size N 2 which can be done in O(N) operations. All the other operations of our marker algorithm deal with local computations along the interface and are consequently also of order N. The ratio between the time spent in the marker algorithm and in the computations done on the bulk of the fluid (the Navier-Stokes solver) is then of order 1/N. As the domain size increases the proportion of computational time needed by the marker representation decreases. In practice, the marker computation accounts for less than ten percent of the time for a 64 2 mesh.

Redistribution

As the interface evolves, the markers may drift due to non zero tangential velocities and we may need more markers if the interface is stretched by the flow. We then need to redistribute the markers in order to ensure an homogeneous distribution of points along the interface. This is done at each time step using the interpolating curve (x(s), y(s)). As s is an approximation of the curvilinear coordinate, if we choose a redistribution length l, the new number of markers is N new = s N /l and the points are redistributed as (x new i , y new i ) = (x(il), y(il)). l is usually chosen as ∆x, which yields an average number of one marker per computational cell. Decreasing this length does not apparently improve the accuracy and in some cases leads to instabilities.

Computing the volume fraction

In the case of flows with varying density and/or viscosity between the phases, we need to solve the advection equations ( 4) and ( 5). This can be reduced to the advection of the volume fraction field C. We then have to create the new volume fraction field corresponding to the new parametric interpolation of the interface. In the cells crossed by the interface, we want to compute the area of the enclosed domain Ω (figure 6). This problem can be reduced to the computation of a circulation along a parametric curve as illustrated below. Following Stokes theorem, if P and Q are two functions of the space

if we choose P = 0, Q = x we get

∂Ω xdy = Ω dxdy, (30)
which can be written using the parametric description of ∂Ω x(s)y ′ (s)ds = Ω.

For our third order polynomial parametric function (x(s), y(s)),

x(s) = a x s 3 + b x s 2 + c x s + d x y(s) = a y s 3 + b y s 2 + c y s + d y , (32)
this integral is given by

x(s)y ′ (s)ds = 1 2 a x a y s 6 + 1 5 (3b x a y + 2a x b y )s 5 + 1 4 (a x c y + 2b x b y + 3c x a y )s 4 + 1 3 (b x c y + 2c x b y + 3d x a y )s 3 + 1 2 (c x c y + 2d x b y )s 2 + d x c y s. ( 33
)
We must then add the circulation along the faces of the cell. In the case of our cartesian grid, three useful observations can be made:

1. As for the horizontal faces dy = 0, only the vertical faces can lead to a non-zero contribution to the circulation.

2. If a vertical face is not crossed by the interface its contribution to the circulation is xdy = x dy = ±x, which is an integer number.

3. We impose a maximum CFL number of 0.5. In this case the variation of the volume fraction during any time step can not be larger than 0.5 in any cell. Consequently, we can write

c n -0.5 ≤ c n+1 ≤ c n + 0.5 (34)
If we compute only the fractional part of the circulation, cn+1 , c n+1 is defined as cn+1 + i, where i is an unknown integer. However i is uniquely defined by (34) and can take only four distinct values {-1, 0, 1, 2}.

It is now easy to compute the volume fraction, we just follow the interface, as defined by the list of polynomials. We test whether the interface intersects the horizontal and/or vertical grid lines. The point of intersection is computed using a combination of bisection and Newton-Raphson methods [START_REF] Press | Numerical Recipes, The Art of Scientific Computing[END_REF]. If there is no intersection (eg. segment BC in figure 6), we just add the circulation along the whole polynomial segment, computed using (33), to the cell volume fraction c. If there is one (or more) intersection, we compute the circulation for each part of the polynomial segment and add the circulation to the corresponding cell (segments AB and CD). When the intersection is on a vertical segment, we add the vertical contribution to the circulation to one cell and subtract it from the neighboring cell (figure 7). Once all the interfaces contributions have been added, we use (34) to compute the value of the integer constant i to be added to get c n+1 . With this algorithm, the cases where an interface (or several interfaces) crosses more than once the same cell are treated implicitly.

Cells which were crossed by the interface at the previous time step have their volume fraction field updated to zero or one according to whether they now lie outside or inside the domain respectively. The inequality (34) gives the correct value.

Surface tension contribution to the momentum equation

Following the tension formulation [START_REF] Pawar | Marangoni effects on drop deformation in an extensional flow: The role of surfactant physical chemistry. I. Insoluble surfactants[END_REF], if the interface crosses a momentum control volume boundary, we compute the value of the parameter s (using the same solver as above) for the intersection point. We then add the surface tension contribution σt where t is the unit tangent to the curve pointing outward for the considered cell: t = ±(x ′ (s)i+y ′ (s)j)/ x ′2 (s) + y ′2 (s), where i and j are the unit vectors in the x and y direction. As this computation is done once for the two neighboring cells, the momentum contributions of the surface tension cancel exactly. As the control volumes for each component of the momentum do not coincide due to our staggered grid (see figure 3), we need to repeat this calculation twice, shifting the position of the markers respective to the C control volume, by (-1/2, 0) for the ρu component and (0, -1/2) for the ρv component.

The pressure gradient correction

In the case of our staggered grid, the x component of the integral contribution of the pressure gradient is usually given by Ωx ∂p/∂x = h(p i,j -p i-1,j ). If the interface is crossing the vertical face BC of the control volume, we have a pressure jump either between p i,j and p i,j-1 (case (a): EB < h/2), or between p i,j and p i,j+1 (case (b): EB > h/2). We can have a better approximation of the flux hp BC , using a linear correction δhp i,j to the pressure fluxes on either sides of the interface:

• if EB < h/2: hp BC = EBp i,j-1 + ECp i,j = hp i,j + δhp i,j = hp i,j + EB(p i,j-1 -p i,j ).

• if EB ≥ h/2: hp BC = EBp i,j + ECp i,j+1 = hp i,j + δhp i,j = hp i,j + (h - EB)(p i,j+1 -p i,j ).
A correction δhp i,j -δhp i-1,j can then be applied to the integral term Ω ∂p/∂x. This is done practically when computing the surface tension contribution expressed at the same intersection point E. Ideally, this pressure gradient correction should be applied when solving the pressure equation ( 13), thus getting the right pressure field (where hp BC = hp ij ). We have chosen however a simpler approach where the pressure gradient correction is considered as a source term in the right hand side of the momentum balance equation ( 10). This term is then computed using the pressure field at time n∆t, which is less accurate and probably less stable. We will see however that the method gives good results for difficult test cases. This pressure gradient correction is probably the simplest choice available. However several other solutions exist and it would be interesting to assess their accuracy and convergence properties. A simple extension is the choice of higher order interpolation schemes on both sides of the interface, in order to obtain a more accurate representation of the pressure jump. The mean square interpolation proposed by Shyy et al. [START_REF] Shyy | Computational Fluid Dynamics with Moving Boundaries[END_REF] would be interesting too and could be used to obtain any term with rapid variations across the interface (e.g. viscous terms for high contrasts of viscosity).

Results

In this section, we present some simple test cases intended to illustrate the ability of our method to cope with high tension surface flows without loss of accuracy. All the computations have been done on an IBM RS6000/370 workstation. Typically one time step on a 128 2 grid requires 0.7 second, the most time consuming procedure being the multigrid solver.

Stationary bubble and spurious currents

An interesting test case is the verification of the stationary Laplace solution for a circular bubble or droplet. As shown previously, in principle a stationary solution exists for our method and consequently Laplace law should be verified exactly. However, one cannot interpolate a circle (which is C ∞ ) exactly using a parametric spline curve (C 2 in our case) and small differences will subsist. Two questions follow, the stability of this approximate solution and its numerical convergence to the theoretical solution.

Using dimensional arguments one can show that the amplitude of the spurious currents must be proportional to σ/µ. A relevant measure of the amplitude of the spurious currents is then the coefficient of proportionality K = U max µ/σ which is characteristic of the accuracy of the method. Using this velocity scale, we can define a Reynolds number as σρD/µ 2 , where D is the diameter of the bubble. This number is actually 1/Oh 2 , where Oh is known as the Ohnesorge number [START_REF] Ohnesorge | Formation of drops by nozzels and the break-up of liquid jets[END_REF]. Table 1.a illustrates the constant character of K over a broad range of Ohnesorge numbers. We used a 32 2 Cartesian mesh with periodic boundary conditions in the x direction and reflecting conditions on the horizontal walls, the density and viscosity ratios are one, the ratio between the diameter and the box width is 0.4. We measured the maximum amplitude of the spurious currents after 250 characteristic time scales (t = t phys σ/D/µ = 250).

1/Oh

Table 1: Amplitude of the spurious currents around a circular bubble. (a) Independence of the non-dimensional maximum velocity with the Ohnesorge number on a 32 2 mesh. (b) Convergence of the non-dimensional maximum velocity with the spatial resolution.

A second series of tests was performed with the same geometry, for a value of 1/Oh 2 of 12,000 and increasing spatial resolutions. Table 1.b shows a convergence of the method to the theoretical solution which is faster than second order in the spatial resolution (the exponent is approximately 2.75). For the maximum resolution of 256 2 the absolute error is of the order of the machine precision.

The same test was done using the Piecewise Linear Interface Calculation (PLIC/CIAM) code [START_REF] Li | Calcul d'interface affine par morceaux[END_REF]. Figure 9 shows that without the correction to the pressure gradient, the amplitude of the spurious currents is independent of the spatial resolution for PLIC, which results in spurious currents 10 5 times stronger than with the markers method for the 256 2 case.

Capillary wave

Another simple but important test is the solution of the damped oscillation of a capillary wave. The linear theory for the small amplitude oscillation of an interface between two inviscid fluids of equal density ρ in a non-bounded domain, gives the dispersion relation

ω 2 0 = σk 3 2ρ , ( 35
)
where k is the wavenumber [START_REF] Landau | Fluid Mechanics[END_REF].

Normal mode analysis

However, in the case of a viscous fluid, this relation must be corrected. To derive this new dispersion relation, we can use the formalism first introduced by Meunier [?], which relates the time derivative Ḣk (ω) of the amplitude to the forcing F k (ω) applied to the interface via a transfer function R k (ω). For small amplitude perturbations, we have

Ḣk (ω) = R k (ω)F k (ω), (36)
with

R k (ω) = ik(1 -k/q) 2ρω , q = k 2 -iρω/µ. (37)
We look for solutions of the form H k (ω) = H 0 e i(kx-ωt) , we then have the relation Ḣk (ω) = -iωH k (ω). The forcing F k (ω) due to the surface tension is

F k (ω) = -σk 2 H k (ω). ( 38
)
Substituting this relations in (36) yields the dispersion relation

ω 2 = σk 3 2ρ (1 -k/q) = ω 2 0 (1 -k/q), (39)
which gives the viscous correction to the inviscid case (35) and the damping ratio of the amplitude. We can write (39) with q as unknown, this yields a polynomial equation

q(q -k)(q + k) 2 = - σk 3 ρ 2µ 2 , (40)
which can be solved numerically.

Initial value problem

However, as will be shown from the numerical simulations, this normal mode analysis is not suitable for the study of the related initial-value problem.

As shown by Prosperetti et al. [START_REF] Cortelezzi | Small amplitude waves on the surface of a layer of a viscous liquid[END_REF][21], an analytical solution for the initial value problem exists in the case of the small-amplitude waves on the interface between two superposed viscous fluids, provided the fluids have the same dynamic viscosity. We consider here the case of an initial sinusoidal perturbation of amplitude H 0 , the two fluids being at rest. If we take as fundamental time ω -1 0 , with ω 0 given by ω 2 0 = σk 3 /(ρ 1 +ρ 2 ), and as fundamental length k -1 , we can set τ = ω 0 t, ǫ = νk 2 /ω 0 .

Using these non-dimensional time and viscosity, the analytical solution for the non-dimensional amplitude a = H/H 0 is given in compact form by

a(τ ) = 4(1 -4β)ǫ 2 8(1 -4β)ǫ 2 + 1 erfc(ǫ 1/2 τ 1/2 ) (42) + 4 i=1 z i Z i ω 2 0 z 2 i -ǫω 0 ×exp[(z 2 i -ǫω 0 )τ /ω 0 ]erfc(z i τ 1/2 /ω 1/2 0 ),
where the z i 's are the four roots of the algebraic equation

z 4 -4β(ǫω 0 ) 1/2 z 3 + 2(1 -6β)ǫω 0 z 2 + (43) 4(1 -3β)(ǫω 0 ) 3/2 z + (1 -4β)(ǫω 0 ) 2 + ω 2 0 = 0,
and Z 1 = (z 2 -z 1 )(z 3 -z 1 )(z 4 -z 1 ) with Z 2 , Z 3 , Z 4 obtained by circular permutation of the indices. The dimensionless parameter β is given by β = ρ 1 ρ 2 /(ρ 1 + ρ 2 ) 2 .

Comparison with the numerical simulations

The test case for the capillary wave is a square box divided in two equal parts by a sinusoidal perturbation. The wave length is equal to the box width. The boundary conditions are free slip on the top and bottom walls, periodic along the horizontal axis. The ratio between the initial interface perturbation H 0 and the box height is 0.01. The Ohnesorge number is 1/ √ 3000, the non dimensional viscosity ǫ ≈ 6.472 × 10 -2 , the two fluids densities are the same and thus β = 1/4. Figure 10 shows the evolution of the amplitude with time for the numerical solution and both the normal mode and initial value analytical solutions. The normal mode curve is obtained using H 0 e λτ /ω 0 cos(ωτ /ω 0 ) where ω and λ are given by the numerical solution of (39).

Grid

= 10ρ 2 , β ≈ 0.826, ǫ ≈ 4.799 × 10 -2 , Oh 1 = 1/ √ 3000, Oh 2 = 1/ √ 30000.
In order to study the spatial convergence of the method, we have repeated this test for resolutions of 8 2 , 16 2 , 32 2 and 64 2 . Table 2 summarizes the results. We obtain a good convergence and a relative error of order 10 a spatial resolution of 32 2 compared to the initial value analytical solution of Prosperetti. Moreover, we treat here the case of a bounded domain of aspect ratio hλ = 1/2. The analytical solution of Prosperetti is given in the infinite depth case and should be corrected with a term of order cotanhkh = cotanhπ ≈ 1.0037.

The previous test was performed with a density ratio of one (β = 1/4). Table 3 illustrates the results for a density ratio of ten (β ≈ 0.826).

Rayleigh-Taylor instability

In order to illustrate the capability of our method to deal with more complex cases, we present here a classical test. The original version of this computation has been performed by Puckett et al. [START_REF] Puckett | A second-order projection method for tracking fluid interfaces in variable density incompressible flows[END_REF][24] using a VOF type method. A one meter wide, four meters high rectangular domain is discretized using a 64 × 256 grid. The fluid densities are 1.225 and 0.1694 kg/m 3 . The fluid viscosities are 0.00313 kg.m -1 .s -1 . The interface between the fluids is an initially sinusoidal perturbation of amplitude 0.05 meters. Figure 12 shows the evolution of the interface at times 0, 0.7, 0.8 and 0.9 seconds. The maximum mass fluctuation of our method is approximately 0.14%, which is larger than the observed variation for a VOF type method (about 0.01%). The interface evolution compares well with the results of Puckett et al. [START_REF] Puckett | A second-order projection method for tracking fluid interfaces in variable density incompressible flows[END_REF] and with the simulation done using PLIC/CIAM on a 128 × 512 grid (figure 13).

We can note however, that we have not yet implemented a reconnection mechanism to deal with topology changes. Consequently, our method conserves the small filamentary structures even for low resolutions. We are studying a reconnection mechanism based on a physical representation of the short range interactions between interfaces. Whereas VOF type methods introduce an arbitrary reconnection length (usually the mesh spacing), this method should provide a physical criterium for the reconnections. The marker representation gives the precise location of the interface and is therefore well-suited for the study of small scale effects.

Conclusions

In summary, we have presented a front tracking numerical algorithm able to deal accurately with surface tension terms. For the case of a stationary t=0.9 t=0.7 t=0.0 t=0.8

Figure 12: Rayleigh-Taylor instability on a 64 × 256 grid using the fronttracking algorithm.

t=0.9 t=0.7 t=0.0 t=0.8

Figure 13: Rayleigh-Taylor instability on a 128 × 512 grid using the PLIC/CIAM algorithm.

bubble, the method shows a convergence to the theoretical solution which is faster than second order in the spatial resolution. For reasonable mesh sizes, the spurious currents usually present when using fixed grids can then be reduced to the machine precision. Good results are obtained for the test case of a capillary wave with viscous damping. We obtain a very good agreement with Prosperetti's linearized solution to the initial value problem. The relative error for the time evolution of the amplitude is smaller than 10 -2 when using 32 or more points per wavelength. The method is robust and remains accurate even for very small Ohnesorge numbers (10 -3 or less). As we use a fixed Eulerian grid, the Navier-Stokes algorithm is efficient. Moreover, due to its one dimensional nature, the surface tracking part is negligible in term of computational time. The ideas developed in this paper can be generalized to different types of Navier-Stokes solvers and are not limited to projection methods on cartesian grids. In particular, we emphasized the necessity of a precise and consistent discretization of the surface tension and of the associated pressure jump. This probably requires the precise knowledge of the interface position relative to the underlying grid. Various applications require an accurate representation of the surface effects: motion of drops of water in air [START_REF] Zaleski | Two-dimensional navier-stokes simulation of deformation and break-up of liquid patches[END_REF][22], wave breaking [START_REF] Banner | Wave breaking in deep water[END_REF] or sonoluminescence [?]. Moreover, the exact description of the interface position allow an easy implementation of the membrane models found in biological mechanics [START_REF] Barthes-Biesel | Microrheological models of red blood cell mechanics[END_REF][10] and of the effect of surfactants [START_REF] Pawar | Marangoni effects on drop deformation in an extensional flow: The role of surfactant physical chemistry. I. Insoluble surfactants[END_REF]. An interesting study could also be a detailed sub-grid modelling of the reconnections between interfaces [START_REF] Eggers | Universal pinching of 3d axisymmetric free-surface flow[END_REF].

We are working on an axisymmetric version of the code. In the future, we would like to develop a fully three-dimensional front tracking algorithm.

Figure 2 :Figure 3 :

 23 Figure 2: MAC discretization of the pressure, volume fraction, momentum and velocity components

Figure 5 :

 5 Figure 5: The discretization of the pressure gradient. If a simple discretization is used near the interface p1 = p2 = p3 = p4 and the pressure gradient h -1 (p1 -p4) cannot balance the surface tension.

Figure 6 :

 6 Figure 6: computation of the volume fraction

Figure 7 :

 7 Figure 7: Intersection of an interface with a vertical face. Two cases arise depending on the topology of the domain.

Figure 8 :

 8 Figure 8: The pressure gradient correction for a MAC grid. The location of the interface relative to the pressure nodes yields two cases: (a): EB < h/2, (b): EB > h/2.

Figure 9 :

 9 Figure 9: Amplitude of the spurious currents versus spatial resolution for the PLIC and markers methods.

Figure 10 :

 10 Figure10: Time evolution of the amplitude of a capillary wave in the viscous case. Our front-tracking method is used on a 128 2 square grid. The theoretical curves are obtained from a normal mode analysis and from the exact solution to the initial value problem in the linearized viscous case (Prosperetti).

Figure 11 :

 11 Figure 11: Close up of the previous figure illustrating the differences between the normal mode an initial value solutions.

Table 2 :

 2 Evolution of the relative error between numerical computations and the normal mode and initial value analytical solutions. The relative error is the r.m.s. of the differences between the solutions divided by the amplitude of the initial perturbation. We compare the solutions for the first 25 nondimensional time units. ρ 1 = ρ 2 , β = 1/4, ǫ ≈ 6.472 × 10 -2 , Oh = 1/

	size Error/Initial value Error/Normal mode
	8 2	0.2972	0.3077
	16 2	0.0778	0.0959
	32 2	0.0131	0.0332
	64 2	0.0098	0.0307
	128 2	0.0065	0.0280
			√	3000.
	Grid size Error/Initial value
	8 2		0.3593
	16 2		0.1397
	32 2		0.0566
	64 2		0.0264
	128 2	0.0148

Table 3 :

 3 Evolution of the relative error between numerical computations and the initial value analytical solution. ρ 1


E Ory

H Yuan

A Prosperetti

S Popinet

S Zaleski

Growth and Collapse of a Vapor Bubble in a Narrow Tube

Keywords: Nos.: 47.55.Dz, 47.60.+i, 47.11.+j

The uid mechanical aspects of the axisymmetric growth and collapse of a bubble in a narrow tube lled with a viscous liquid are studied numerically.

The tube is open at both ends and connects two liquid reservoirs at constant pressure. The bubble is initially a small sphere and growth is triggered by a large internal pressure applied for a short time. After this initial phase, the motion proceeds by inertia. This model simulates the e ect of an intense, localized, brief heating of the liquid which leads to the nucleation and growth of a bubble. The dimensionless parameters governing the problem are discussed and their e ects illustrated with several examples. It is also shown that, when

I. INTRODUCTION

In an ink{jet printer drops are produced and ejected due to the rapid growth of bubbles in a narrow ow passage just upstream of the nozzle exit (see e.g. Ref. 2,1,5]). This is but one example of the possible application of vapor bubbles in small devices, where they o er the interesting potential of mechanical actuation, pumping, and other ow e ects without mechanical moving parts. In these situations bubbles are produced by small heaters in the tube wall to which a brief current pulse is applied. A very small mass of liquid reaches a temperature comparable with the critical temperature, and a vapor bubble nucleates and rapidly grows. In a recent paper we have studied the thermal and attendant phase{change processes that take place in such a situation 18]. We found that, due to the work of expansion and vapor condensation on the cold liquid surface, the initial high vapor pressure in the bubble is very quickly dissipated. The rest of the bubble evolution essentially proceeds by inertia. In that study the uid mechanical aspects of the problem were treated by means of a simpli ed model based on potential ow and an approximate treatment of viscous e ects. The purpose of this article is to simulate the actual ow behavior by solving numerically the full Navier{Stokes equations. We shall however simplify the thermal problem by assuming that the bubble internal pressure at time 0 is brie y raised to a level p above the saturation level p 0 ' 0 corresponding to the undisturbed liquid temperature. This high internal overpressure is maintained for a time , after which it falls to zero again. In another recent paper 19] we have compared the results of this simple model with the one of Ref. 18] nding a very good agreement.

The advancing of a gas bubble in a tube is a classic problem treated by many others since the early work of Bretherton 3]; Ref. 10] gives a good review and some more recent work can be found in Refs. [START_REF] Chikanawa | Generation, growth and collapse of a bubble in conductive liquid heated directly by electric current[END_REF]7]. The situation considered in this paper is however di erent as, rst, the bubble undergoes a signi cant volume change and, second, the ow is highly transient as opposed to the steady or quasi{steady situations considered before. In addition, inertial e ects are signi cant and therefore we deal with the full Navier{Stokes equations, rather than the Stokes equations as many of the previous authors. The situation considered in the study of boiling in porous media (see e.g. Ref. 12,15,11,[START_REF] Scardovelli | Direct numerical simulation of free{surface and interfacial ow[END_REF]) is also quite di erent from the present one as, in that case, heat is added to the system continuously from the entire pore wall, rather than locally and impulsively as here. Among the earlier work on the speci c problem studied here, the closest one is that of Asai et al. 2] who, however, developed a quasi{one{dimensional model and were chie y concerned with drop formation in the ink{jet system rather than with speci c uid mechanical phenomena.

An interesting and unexpected result that we nd is that a single bubble cyclically growing and collapsing away from the midpoint of the channel is capable of inducing a mean uni{directional ow capable of pumping liquid from one reservoir to the other.

II. PRELIMINARY CONSIDERATIONS

We consider a tube of length L and inner diameter D connecting two liquid reservoirs at the same pressure P 1 (Fig. 1). We simplify the problem by assuming axial symmetry and, at the initial instant, place a small spherical bubble of diameter D i on the axis at the center of the tube. The pressure in the bubble is maintained at the level p 0 + p for 0 t < , after which it falls to p 0 for the remainder of the growth and collapse cycle; p 0 is the saturated vapor pressure at the undisturbed liquid temperature.

On the basis of a simple conceptual model according to which, in its growth, the bubble pushes in opposite directions two equal liquid columns of length 1 2 L with an acceleration ( p + p 0 P 1 )= 1 2 L , we nd that a characteristic velocity scale for the problem is

V = p + p 0 P 1 1 2 L ; (1)
where is the duration of the high{pressure phase. In terms of this characteristic velocity we de ne Reynolds, Strouhal, and Weber numbers in the usual way by

Re = V D ; St = D V ; W e = D V 2 :
(2)

Here , , and are the density, viscosity, and surface tension coe cients of the liquid.

It should be noted that the velocity V de ned by (1) is of the order of the peak velocity during the process and, therefore, the value of the Reynolds number de ned here tends to underestimate the magnitude of viscous e ects. In the literature the capillary number Ca = W e=Re is also frequently encountered, but its signi cance here is somewhat obscured by the highly transient nature of the phenomenon investigated.

Other parameters characterizing the problem are the ratio of the applied overpressure to the pressure di erence P 1 p 0 that tends to collapse the bubble: P = p P 1 p 0 ;

(3) the aspect ratio of the channel:

A = L D ; (4)
and the ratio of the initial bubble diameter D i to the channel diameter:

D = D i D : (5)
Other quantities of interest can be expressed in terms of the previous parameters. For example, the ratio of the viscous di usion time to the capillary wave period is

T d = D 2 s D 3 = Re
p W e : [START_REF] Chikanawa | Generation, growth and collapse of a bubble in conductive liquid heated directly by electric current[END_REF] This quantity (which equals the square root of the Laplace number or the inverse of the Ohnesorge number), may also be interpreted as the order of the number of periods necessary for viscosity to damp out the longest capillary waves that the system can support.

We shall compare the results with those of a simpli ed model developed in earlier work in which the bubble is approximated as a cylinder occupying the entire cross section of the tube, with some provisions made for surface tension and viscous e ects 18]. A full description and a complete numerical solution are presented in 19]. If the model is further simpli ed neglecting surface tension e ects, interfacial curvature, and viscosity, it lends itself to an approximate analysis and a consideration of the results found in this way is useful to interpret the numerical ones presented below. For simplicity we only consider the simpli ed solution valid provided the bubble surfaces remain su ciently far from the tube ends. With this assumption it can be shown that the maximum horizontal extension X m of the bubble is: X m ' p + p 0 P 1 P 1 p 0 p 2 1 2 L : (7) In the context of this simple model the bubble is a cylinder and therefore its maximum

volume is V max = 1 4 D 2 X m or, in dimensionless form, V max 1 4 D 2 L = P A St : (8)
On the basis of the same simple model one nds that the processes of growth and collapse are nearly the mirror image of each other and that the bubble lifetime t B is given by t B ' 2 P ; (9) or

V t B L = 2 P A St : (10)
With the aid of these relations we have the following expression for the ratio of the viscous di usion length to the tube diameter:

T v = p t B D = s 2P
St Re : (11) This parameter can also be interpreted as the ratio of the viscous di usion time to the bubble lifetime. Another signi cant parameter is the ratio of the period of capillary waves with wavelength of the order of D to the bubble lifetime:

T = s D 3 t 2 B = St 2P
p W e : (12) This parameter gives an idea of the number of capillary wave periods during the bubble life.

III. THE MODEL

We solve numerically, by a method brie y described later, the standard full, incompressible, viscous Navier{Stokes equations in the interior of the tube. Since only a minute amount of liquid is heated, and only for a brief fraction of the duration of the process, we treat the uid properties constant. The usual no{slip condition is applied at the walls. At the bubble surface the normal stress is taken to equal the internal bubble pressure corrected for surface tension e ects while the tangential stress vanishes. The boundary conditions at the tube ends are less straightforward and require some discussion.

The simplest option is to impose that the pressure at the tube ends is the undisturbed reservoir pressure P 1 . This choice can be justi ed during the growth phase of the bubble as, due to ow separation, the uid issues from the tube ends similarly to a jet. Approximating the ow as parallel, and neglecting viscous e ects in the ow induced by the jet in its surroundings, we then nd that the pressure at the tube outlets is indeed approximately equal to P 1 . The accuracy of this prescription during the collapse phase, however, is harder to justify as the liquid now enters the tube ends essentially omnidirectionally in a fashion reminiscent of a sink ow. If we approximate the ow as potential, an application of the Bernoulli integral shows that the pressure at the tube's ends (index e) is given by @ e @t + 1 2 u 2 e + P e = P 1 ; (13) where e is the velocity potential and u e the local velocity, assumed uniform over the cross section. In the context of potential ow, it can be proven that the velocity potential e satis es 9,8] e 1 2 Du e = 0 ; (14) where the upper sign is for the left end and the lower sign for the right end of the tube. With this relation we then nd P e = P 1 1 2 u 2 e 1 2 D @u e @t : (15) The last term can be shown to be negligible on the basis of the following argument. Integrate the momentum equation along the axis of the tube from the bubble surface to the tube end assuming a uniform velocity. This procedure gives `(t) du e dt

= P e P b ; [START_REF] Satik | A pore{network study of bubble growth in porous media driven by heat transfer[END_REF] where `is the length of the liquid column between the free surface and the tube end; the pressure in the liquid at the bubble surface, P b , di ers from the bubble internal pressure p 0 because of surface tension. Upon using this result to eliminate the liquid acceleration in (15) we nd

P e P 1 = 1 2 u 2 e 1 2 D ` (t)
(P e P b ) : [START_REF] Scardovelli | Direct numerical simulation of free{surface and interfacial ow[END_REF] In this paper we only consider situations such that D `, and we therefore feel justi ed in dropping the last term. In principle, at the end of the collapse, if opposite sides of the bubble were to collide, P b would suddenly become very large and this approximation could be inaccurate. We do not include this nal stage of the collapse in the present simulations. In summary, then, we set P e = P 1 during bubble growth while, during collapse, we impose

P e = P 1 1 2 u 2 e : (18)
The switching between the two conditions is done at the moment of maximum bubble expansion when the velocity is essentially zero, which maintains a desirable degree of continuity.

To estimate the relative importance of the second term in [START_REF] Wilson | The unsteady expansion and contraction of a long two{dimensional vapour bubble between superheated or subcooled parallel plates[END_REF] it is useful to rewrite this equation in dimensionless form: (20

)
Since V is of the order of the maximum velocity and, for the cases considered here, the combination of dimensionless parameters appearing in this equation is not large, one would only expect a small contribution from this end correction, which is con rmed by the numerical results presented later in connection with Fig. 4.

All the results presented here are computed using a free{surface version of the code described in [START_REF] Yuan | The pumping e ect of growing and collapsing bubbles in a tube[END_REF][START_REF] Satik | A pore{network study of bubble growth in porous media driven by heat transfer[END_REF]14]. The full incompressible axisymmetric Navier{Stokes equations are solved by a projection method. The interface is discretized using a set of marker points linked by cubic splines. This description supplies the accurate knowledge of the position and curvature of the interface necessary to include the surface tension contribution. The marker points are advected by a simple bilinear interpolation scheme and are redistributed at every time step to insure a uniform distribution as the bubble deforms. Axial symmetry is enforced in the calculation. In the radial direction the integration domain extends between the axis and the wall of the tube, while in the axial direction it includes the entire tube length.

The interface velocity is calculated by extrapolation from the liquid region, which requires that a su cient number of grid points exist between the bubble surface and the tube wall.

To avoid an excessive grid re nement as this liquid layer gets thinner, we arbitrarily set the outward radial velocity of the interface to zero when it has advanced as far as a few grid points from the tube wall. By varying the number of grid points between 1 and 3, and by using di erent grids, we have veri ed that this procedure leads only to minute di erences that are undetectable within the thickness of the lines used in the graphs. When the bubble is as close as this to the tube wall, the radial velocity of the interface is very small anyway. We have also conducted the standard grid re nement tests to verify the grid independence of the solutions that will be presented in the next section. Typically we used square cells with 32 grid points in the radial direction.

IV. RESULTS

Figure 2 shows the bubble shape at di erent instants of time during growth (upper panel) and collapse. In this as well as in all the examples that follow the bubble originates at the midpoint of the tube and therefore we only show the right half of the bubble. Note that the tube extends well beyond the right boundary of the gures all the way to z=D = 16. The dimensionless parameters have the following values: Re = 306.25, St = 12.76, W e = 21.44, P = 50, A = 32, D = 0.1 (see Table ). Here and in the following gures (except the last two) the duration of the internal overpressure is V =L = 2.45 10 3 .

It may also be useful to have in mind a concrete physical realization of this system in terms of dimensional quantities. With the physical properties of water at 20 C, the previous parameters would correspond e.g. to a 2 mm{long tube, with a diameter of 62.5 m. The reservoir pressure is P 1 = 1 atm, the overpressure p is 50 bars, and the overpressure duration is = 1 s; the characteristic velocity V is 4.9 m/s and the initial bubble diameter 6.25 m.

As mentioned before, in the practical applications we envisage the initial temperature of the bubble surface would be of the order of 80% the critical temperature (see e.g. 1]). The corresponding saturation pressure would be much higher than the value of 50 bars used here. However, in that case, the bubble would start from a small nucleus and a large part of the initial overpressure would be expended in balancing surface tension and in expanding the nucleus. Since here we start already from a macroscopic bubble, and we keep the pressure overpressure constant until it falls to zero, rather than allowing it to gradually decay, the value that we use is realistic. This point was examined in Ref. 19] where it was shown that, assuming the product p equal to the time integral of the bubble internal overpressure resulted, in a good comparison with the results of a thermally more sophisticated model in which the initial bubble temperature was higher, but the cooling of the vapor was accounted for 18].

In Fig. 2 successive pro les are separated by V t=L = 3.92 10 3 dimensionless units (corresponding, for the dimensional values given before, to 1.6 s); hence the internal pressure has already fallen to p 0 before the second pro le shown.

The bubble remains approximately spherical until it grows to occupy about 2/3 of the cross section. At that point preferential growth in the axial direction begins, with a relatively high curvature that gradually becomes more blunt as the ow slows down. The curvature of the end surfaces of the bubble during collapse remains relatively small and the motion is very nearly one{dimensional until the very end. The rapid reversal from an oblate to a prolate shape between the next to the last and the last contours shown is reminiscent of a similar phenomenon encountered in the collapse of a free, initially oblate bubble 4]. The phenomenon is caused by the fact that, when inertia dominates, the liquid velocity tends to be inversely correlated with the local radius of curvature of the interface.

The maximum dimensionless volume of the bubble predicted by ( 8) is in this case 0.122, to be compared with a calculated dimensionless volume of 0.076. In addition to the approximations inherent in (8), this di erence is due to the fact that the precise bubble shape is rounded at the ends, which further decreases the volume with respect to the cylindrical shape assumed in (8). The bubble dimensionless lifetime estimated from (10) is 0.245, while the computed one is 0.187. The di erence between the two is enhanced by the fact that the bigger bubble predicted by the simple model takes longer to collapse and has, therefore, a longer lifetime. The parameter T v de ned in (11) has the value 0.16, implying that viscous e ects are only signi cant over a small fraction of the channel diameter. The thickness of the residual lm near the wall is limited by the numerical considerations described in the previous section and therefore is not quantitatively signi cant here.

The capillary number evaluated from the formal relation Ca = W e=Re has the value 0.07. However, by using directly the de nition Ca = V= with the computed time{dependent velocity, one nds values ranging approximately between 0.19 and 0, with most of the time Ca < 0.03{0.04. Thus surface tension is considerably less important than viscosity although, due to the transient nature of the ow, the results applicable to the steady case (see e.g. Ref. [START_REF] Chikanawa | Generation, growth and collapse of a bubble in conductive liquid heated directly by electric current[END_REF]7]) are not directly relevant here. The parameter T de ned in (12) equals 0.59, so that the longest capillary waves have a period comparable with the bubble characteristic time, which explains the relative bluntness of the bubble ends. The parameter T d of ( 6) is about 66, implying that capillary waves are only very lightly damped by viscosity.

A sixfold decrease in the Reynolds number from 306.25 to 51.04 (Fig. 3) produces a bubble about 20% shorter with a dimensionless lifetime of 0.159. Neither e ect is contained in (7) and ( 9) which are evidently inapplicable to this relatively high{viscosity case. In qualitative agreement with the increased value of T v , the liquid lm now is not limited by the numerics and is thicker than before, but the other general features of the process are not signi cantly a ected. The increased damping of capillary waves (T d is now about 11) permits a greater curvature of the bubble end surfaces. With a larger Reynolds number, one nds that the bubble maximum volume increases but little else changes with respect to the case of Fig. 2. The bubble volume vs. time is shown in Fig. 4 for di erent values of the Reynolds number. Here, as well as in Figs. 5, 10, and 13, time is made dimentionless according to

T = V t L = 1 A St t :
(21)

The dashed line shows the e ect of modifying the boundary condition at the tube ends discussed in section III by simply imposing P = P 1 during both growth and collapse. The e ect on the bubble dynamics is evidently rather small. The maximum volume and the bubble lifetime are shown vs. the Reynolds number in Fig. 5. The bubble shape is markedly a ected by an increase of the Weber number. Figure 6 is for an increase to W e = 120. In the rst place the bubble grows longer as the pressure drop across the interface necessary to balance the e ect of surface tension is less. Secondly, the liquid lm clearly has a non{uniform thickness as the reduced surface tension decreases the pressure gradient in the liquid due to curvature. Finally, the larger value of T allows for the presence of capillary waves that are particularly evident during the collapse phase. If the Weber number is reduced to 10.72 (Fig. 7) the picture is similar to that of Fig. 2 for W e = 21.44 except that, due to the capillary pressure drop, the bubble does not grow quite as long.

All other things being equal, an increase in the Strouhal number means a shorter duration of the overpressure phase. The bubble grows therefore less, as implied by (8) and as shown in Fig. 8 for St = 76.53. The simple expression (8) predicts for this case a bubble dimensionless volume of 0.02, against the exact value of 0.017. For this case T = 3.54, and therefore capillary waves are prominent on the bubble surface. If St is reduced to 2.39 (Fig. 9) the bubble grows so large that the hypothesis under which the result ( 8) is derived is violated so that the predicted dimensionless volume, 0.65, is quite di erent from the actual computed value of 0.22. Figure 10 shows the bubble volume vs. time for several values of St. According to the relation (10), the dimensionless bubble lifetime should be inversely correlated with St, which is qualitatively veri ed. Quantitatively, (10) is only approximately correct for the larger Strouhal numbers due to the neglect of end e ects.

Figures 11 and12 illustrate the e ect of increasing and decreasing the pressure parameter P whose in uence on the volume vs. time is shown in Fig. 13. The proportionality of the bubble lifetime to P predicted by ( 9) is well veri ed provided the bubble does not grow too long.

Finally, the in uence of the aspect ratio is illustrated in Figs. 14 (A = 8) and 15 (A = 64). When plotted versus z=D there is remarkably little di erence between the two results in spite of the large variation in the aspect ratio. This nding is in approximate agreement with the result (7) that would predict X max =D ' P=St independent of A. The di erence is however in the expected direction when the maximum volume is normalized as in ( 8): the bubble grows to occupy about 22% of the channel volume for A = 8, but only about 4% for A = 64. The bubble lifetime increases by about 30% going from A = 8 to A = 64.

The e ect of the bubble initial diameter is small. It is found that increasing D i =D from 0.1 to 0.8 only leads to a 20% increase in the maximum volume and a 10% increase in the bubble lifetime.

V. SINGLE{BUBBLE PUMPING

A very interesting phenomenon is encountered when the bubble grows and collapses away from the tube midpoint. An example is shown in Fig. 16 where the initial bubble nucleus is positioned 3/8 of the tube length away from the left end of the tube. Here all the dimensional parameters have the same values as in Fig. 2 except that = 2 s. As in the previous cases, the area shown in the gure is a small fraction of the computational domain that extends for 0 z 2000 m.

It is evident from the gure that the axial position where the collapse is completed { 859 m from the tube left end { is di erent from that of the original bubble nucleus { 750 m away from the tube left end. Thus, the bubble moves 109 m (i.e., 1.74 tube diameters) away from the closest end of the tube. Since the bubble essentially occludes the tube during most of its life, this result suggests that the growth and collapse cycle is capable of imparting a net displacement to the liquid. This fact is con rmed by Fig. 17 which shows the liquid velocities on the tube axis at the two ends of the tube vs. time for this case. It is seen here that at the end of the collapse there is a residual positive net velocity, i.e. directed toward the end of the tube farther from the bubble.

Thus, by generating successive bubbles in a periodic fashion, one would be able to pump the liquid from one reservoir to the other one along the tube. Further considerations on this interesting and potentially useful e ect can be found in Ref . 19].

Figure 16 presents also other elements of interest. In the rst place, it is clear that the direction of motion does not reverse at the same time at the two ends of the bubble. Due to the greater inertia of the liquid column on the right, outward motion in this direction continues somewhat longer than that on the left. Secondly, there is a marked di erence in shape between the two ends of the bubbles during growth and, especially, during collapse. Due to the di erent inertia of the liquid on the two sides of the bubble, the velocity with which the two interfaces displace is also di erent. The velocity of the left interface is relatively large, and the shape here has several features reminiscent of those encountered in Fig. 6 where the large value of the Weber number, W e = 114.33, also caused a prevailing of inertia over surface tension. The motion of the right bubble interface is instead relatively slow due to the large inertia of the longer right liquid column. Here, particularly at the beginning of the collapse, one encounters a tendency for the bubble surface to detach from the tube wall which is reminiscent of other slow collapse cases found before, e.g. in Fig. 8 for St = 76.53 or Fig. 12 for P = 12.5.

VI. CONCLUSIONS

We have studied the uid dynamics of the growth and collapse of a bubble in a small tube under the action of a short{lived internal overpressure. In a practical application, one would be interested in using the bubble growth to achieve mechanical actuation, possibly with a high repetition rate. The growth would be triggered by powering an electric heater for a short time to bring the temperature of a small liquid volume in the vicinity of the critical temperature.

In this note we have avoided the complexity of phase change and heat di usion by applying a constant internal overpressure for a short time of the order of 1% of the total bubble lifetime. In spite of this brief duration, the expansion is mechanically very energetic without requiring a large amount of thermal energy. It should therefore be relatively easy to avoid a progressive increase of the system's temperature by dissipating heat by conduction through the tube walls, which is an essential requirement for achieving a fast repetition rate.

As shown in Ref. 18], in an actual application the cooling of the bubble occurs on a scale of microseconds because of the work of expansion and the conduction of heat into the surrounding cold liquid. After this time the temperature of the bubble remains essentially constant at the ambient value, the vapor pressure becomes negligible, and the entire dynamics are governed by inertia. This behavior is a consequence of the fact that, for a liquid such as water, the vapor density at normal temperature is so low that the latent heat required to maintain the bubble lled with saturated vapor can be supplied without altering appreciably the temperature of the liquid at the liquid-vapor interface. The phenomenon is analogous in normal ow cavitation and the reader may consult e.g. Ref. 13] for more details and a quantitative argument.

Finally, a very interesting and somewhat surprising nding of this paper has been the pumping action described at the end of the previous section. This e ect { which we have in the meantime demonstrated experimentally { could have useful practical applications and will be pursued in future work.

  p 0 P 1 the rst term in the right-hand side is very close to 1. The second term can be written

5 .

 5 FIGURES

TABLE I .

 I Values of various dimensionless parameters used in the gures of this article.

				TABLES
	Figure #	P	Re	We	St	A	T v	T d	T
	2	50	306.25	21.44	12.76	32	0.16	66.14	0.59
	3	50	51.04	21.44	12.76	32	0.39	11.02	0.59
	6	50	306.25 114.33	12.76	32	0.16	28.64	1.37
	7	50	306.25	10.72	12.76	32	0.16	93.54	0.42
	8	50	306.25	21.44	76.53	32	0.07	66.14	3.54
	9	50	306.25	21.44	2.39	32	0.37	66.14	0.11
	11	100	306.25	21.44	12.76	32	0.23	66.14	0.30
	12	12.5 306.25	21.44	12.76	32	0.08	66.14	2.36
	14	50	306.25	21.44	12.76	8	0.16	66.14	0.59
	15	50	306.25	21.44	12.76	64	0.16	66.14	0.59
	16 & 17	50	612.5	85.75	3.19	32	0.23	66.14	0.30

TABLE CAPTION I

 CAPTION . Values of various dimensionless parameters used in the gures of this article. Sketch of the situation modelled in this paper. The computational domain is limited to the interior of the tube and is bounded by the dashed lines at the tube ends. 2. Bubble shape at di erent instants during growth (upper panel) and collapse (lower panel); here and in the following the last computed shape in the upper panel is repeated in the lower one. The dimensionless parameters have values A = 32, Re = 306.25, St = 12.76, W e = 21.44 (see Table). Because of symmetry, only the right half of the bubble is shown. Note that the actual computational domain extends between z=D = 16 and z=D = 16 and therefore considerably beyond the left and right frame boundaries. Successive bubble shapes are separated by 3.92 10 3 dimensionless time units.3. As in Fig.2except that Re = 51.04. Successive bubble shapes are separated approximately by 4.54 10 3 dimensionless time units.4. Bubble volume vs. time for di erent values of the Reynolds number; all other parameter values as in Fig.2. The dashed line di ers from the adjacent solid line only in that the condition P e = P 1 is applied during both the growth and collapse phases.

	FIGURE CAPTIONS
	1.

Liquid velocity induced by the bubble growth and collapse on the axis of symmetry at the left (solid line) and right (dashed line) ends of the tube for the case of the previous gure. Note that at the end of the collapse there is a residual velocity away from the left end of the tube.

ACKNOWLEDGMENT

The portion of this work carried out at the Johns Hopkins University was supported by AFOSR under grant F49620{96{1{0386.


David Dureisseix

Folding optimal polygons from squares

What is the largest regular n-gon that fits in a unit square? Can it be folded from a square piece of paper using standard moves from origami? Answering the first question is relatively easy, using simple ideas from geometry. The second is more interesting; our answer illustrates the difference between origami and the standard compass-and-straightedge constructions of the Greeks, where, for instance, the 7-gon cannot be constructed. Not only can we fold a 7-gon, but we can fold the largest one possible from a given square piece of paper.

Origami (from japanese oru, to fold and kami, paper), is the ancient art of paperfolding. When we fold a paper in half, we create a line segment and bisect a length. These simple moves can be combined to reproduce any compassand-straightedge construction [START_REF] Alperin | A mathematical theory of origami constructions and numbers[END_REF][START_REF] Row | Geometric Exercises in Paper Folding[END_REF]. Thus, by origami, as with an unmarked straightedge and compass, we can construct roots of any second-order polynomial from a given unit length.

However, many constructions known to be impossible under the standard Greek rules, such as trisecting a given angle, become possible with origami. For instance, using a construction technique due to Lill and first used for origami by M. P. Beloch [START_REF] Beloch | Sul metodo del ripiegamento della carta per la risoluzione dei problemi geometrici[END_REF], we can construct roots of cubic polynomials by folding [START_REF] Alperin | A mathematical theory of origami constructions and numbers[END_REF][START_REF] Edwards | Folding quartic roots[END_REF][START_REF] Geretschläger | Euclidean constructions and the geometry of origami[END_REF][START_REF] Justin | Résolution par le pliage de l'équation du troisième degré et applications géométriques[END_REF]. Origami also simplifies certain constructions that are possible, but cumbersome, with compass and straightedge.

Since origami often begins with a square piece of paper, we propose not only to fold a regular n-gon, but to fold the one with the largest area that fits in the square. Such polygons will be called optimal polygons. For instance, the side of the largest equilateral triangle that fits in a unit square is known to have length √ 6 -√ 2 ≈ 1.035 (shown in Figure 4). Wetzel [START_REF] Wetzel | [END_REF] takes this as the starting point for his article "Fits and Covers," which gives many answers to similar problems, but does not address our question.

Our first step is to determine the proper orientation of an optimal polygon with respect to the square. We do this in complete generality and then consider how to construct them by folding. We show how to fold the optimal hexagon and pentagon, which can also be constructed with the standard compass and straightedge. Moving into the realm of techniques that break the Greek rules, we trisect an angle and show how to fold the optimal 7-gon and 9-gon, neither of which can be constructed with straightedge and compass alone. It turns out that in each case we fold a stellated polygon as an intermediate step.

Are you eager to fold the optimal 11-gon? If so, you will have to invent a folding technique that permits you to construct roots of a quintic polynomial! Figure 1: The width of a strip that contains a tilted polygon

Facts about optimal polygons

The goal is to find the largest regular n-gon that can be folded from a square piece of paper, for n ≥ 3. Of course, the case n = 4 is trivial, with no folding required. For the general case, let us review some facts about the regular n-gon.

Let R be the radius of the circumscribed circle and r the radius of the inscribed circle. (Another name for r is the apothem of the polygon.) The reader may wish to confirm that r = R cos(π/n). The diameter of the ngon, denoted L, is the maximum distance between any two of its points. A contrasting quantity is the altitude, meaning the shortest perpendicular distance from a vertex to an opposite side. When n is even, these quantities are simple, the diameter is just 2R and the altitude is 2r. When n is odd, the altitude is rather easily seen to be R + r, which is 2R cos 2 (π/(2n)). The diameter can then be found to be 2R cos(π/(2n)). (Figure 1 will help with this.)

Another useful quantity is the side of a stellated polygon, l. Figure 1 shows it to be 2R sin(2π/n). The side of the polygon is h = 2R sin(π/n).

Before fitting our n-gon into a square, we first fit it into a strip. It turns out to be simplest to consider the n-gon to have a fixed radius R and find out how wide the strip must be to contain it. Depending on the orientation of the n-gon, the necessary width will fall somewhere between the diameter and the altitude. Since the altitude is smaller, we might decide that it is best to orient the polygon with its altitude along one dimension of the square. Unfortunately, the n-gon is always fatter in the perpendicular direction.

Therefore, let us find the narrowest strip of paper that can contain a given polygon when it is tilted with respect to the strip at an angle θ. For a fixed rotation angle θ, the minimum strip width is denoted a(θ). The odd case is shown in Figure 1.

Aided by Figure 1, which shows the odd case, the reader can verify the following formula for a(θ), given in terms of the number n of edges (or vertices) of the polygon, and the radius R of the circumscribed circle:

• if n is odd, a(θ) is (π/n)-periodic, and a(θ) = L cos(θ -π/(2n)) for θ ∈ [0, π/n].
Here, L = 2R cos(π/(2n)) is the diameter of the polygon; Now that we know the narrowest strip that contains a tilted regular polygon, we add a second strip of paper, orthogonal to the first one, and require that it too must contain the polygon. Observe from Figure 2 that the width of this second strip will be a(θ + π/2). If the polygon is to fit into a square, each strip must have width A(θ) = max{a(θ), a(θ + π/2)} as shown in Figure 3. Thus, if we minimize A(θ), we find the smallest square, which is equivalent to find the largest regular polygon within a given square. From the previous expressions for a(θ) (depicted in Figure 3), we derive the following values for the side of the smallest square:

• if n is even, a(θ) is (2π/n)-periodic, and a(θ) = L cos(θ -π/n) for θ ∈ a(θ) θ a(θ) θ a ( θ + π / 2 ) Figure 2: Two perpendicular strips 0 θ 3π/2n a(θ) a(θ + π/2) n even n odd a(θ + π/2) θ 0 π/4n π/2n A(θ) a(θ) A(θ) 3π/4n
• θ opt = π/(4n) (modulo π/(2n)) if n is odd, • θ opt = π/(2n) (modulo π/n) if n is even.
Note that in each case, a(θ opt ) = a(θ opt + π/2) and so one can conclude that each side of the square touches at least one vertex of the optimal polygon. Moreover, using the formulas for these angles θ opt , the reader may prove that each optimal polygon has at least one diagonal of the square as an axis of symmetry. Figure 4 shows the optimal polygon placements up to the octagon (n = 8 edges).

We remark that this pattern also gives us the optimal polygons that fit any rectangular piece of paper, and not only square ones. Folding the optimal hexagon and pentagon

The case of the optimal hexagon (n = 6 edges) involves the angle 2π/6, and is very easy to construct. Given a square of paper with side 1, we use the previous formula for θ opt to deduce that R = 1/(2cos(π/12)). This gives us the edge length,

h = R = √ 2( √ 3 -1)/2. The side of the stellated hexagon is l = 2h cos(π/6) = √ 2(3 - √ 3)/2.
Our construction, which first produces the stellated hexagon, is shown in Figure 5.

• Step 1: Fold corner B onto the central vertical line to create line AE, which meets diagonal BD at F . Since sin(∠B AD) = 1/2, this is a simple technique to get an angle ∠B AD = π/6. Then, ∠BAE = ∠EAB = ∠BAB /2 = (π/2 -∠B AD)/2 = π/6 again. (We have trisected the angle ∠BAD = π/2; this is easy for this particular angle, but more difficult for general case, as we will see for the nonagon. Of course, if have only straightedge and compass we cannot do trisect a general angle at all.) Notice that ∠EAC = π/4 -∠BAE = π/12, so,

DF = DO + OF = OA + OA tan(∠EAC) = √ 2(1 + tan(π/12))/2 = √ 2(3 - √ 3)/2
. This is the length l of the desired polygon side. The following steps are needed to move the crease of length l in a correct position to obtain an edge of the optimal hexagon.

y z F C B G D F H A B G { H A B E I | B E I J B E J K } ~ K L B x C A E D B G Figure 7
: Folding sequence of the optimal pentagon the diagonal, GH is an edge of the optimal stellated hexagon.

• Step 3: Unfold. Bring H onto the main diagonal at I, with G as an end point of the crease.

• Step 4: Fold GI and HI.

• Step 5: Complete the stellated polygon, using symmetries.

Folding a pentagon (n = 5 edges) requires the angle 2π/5, and is more difficult than folding the hexagon. As early as 1989, Roberto Morassi [START_REF] Morassi | The elusive pentagon[END_REF] designed an origami construction of the optimal pentagon. The technique shown in Figure 7 has been developed independently and seems much simpler. As before, the stellated version of the polygon is used. With an initial square of unitary edge length, we get l = 1/ cos(π/20) in Figure 6.

• Step 1: With D as the middle of the edge, bisect the angle ∠BAD. The crease is AC, tan(∠DAE) = 1/2 and ∠BAC = ∠CAD = (1/2)∠BAD = (1/2)(π/2 -∠DAE). We can compute tan(∠BAC) = ( √ 5 -1)/2 = BC; this is the so-called "golden ratio."

• Step 2: Fold C on the central horizontal line DG at F , with B as an end point of the crease. Since cos(π/5) = 1/( √ 5 -1) and BG = 1/2, we get BF = BC = ( √ 5-1)/2 = 1/(2 cos(π/5)) = BG/ cos(π/5). This allows us to conclude that cos(∠F BG) = BG/BF = cos(π/5) and ∠F BG = π/5.

• Step 3: Bisect ∠F BG to get ∠HBA = π/10. Unfold.

• Step 4: Bisect ∠ABH (fold behind) to get ∠ABI = π/20 and BI = 1/ cos(π/20). This is the length l of the optimal pentagon edge. As before, the following steps are needed to move the crease in a correct position.

•

Step 5: Bring I on BE at J.

• Step 6: Fold in half BJ to get K. Unfold.

• Step 7: KL = BJ = BI = l is the correct edge.

• Step 8: Complete the polygon.

Folding other optimal polygons

If the optimal square, triangle, and octagon, are the easiest regular polygons to design (the reader may begin to try to fold them, except for the square...), the optimal hexagon and pentagon are the next easiest. As soon as a regular polygon can be constructed by folding, the corresponding optimal polygons can be also folded with a technique similar to the one used in this paper [START_REF] Dureisseix | Searching for optimal polygon -remarks about a general construction and application to heptagon and nonagon[END_REF]. In a publication from 1837, P. L. Wantzel [START_REF] Wantzel | Recherches sur les moyens de reconnaître si un problème de géométrie peut se résoudre à la règle et au compas[END_REF] demonstrated which regular polygons are constructible with straight-edge and compass (see also [START_REF] Carrega | Théorie des corps : la règle et le compas[END_REF]). A necessary condition concerns the number of edges of these polygons: they must have n = 2 p f 1 f 2 ...f s edges, with p as an integer, and where the f i are all different primes of the form 2 m +1, where m is also an integer. This result can be further simplified because a necessary (but not sufficient) condition for these f i to be primes, is to be Fermat numbers, that is, numbers of the form 2 2 m + 1, m still being an integer. Up to the author knowledge, the only known Fermat primes to day are 3, 5, 17, 257 and 65537 [START_REF] Weisstein | The CRC Concise Encyclopedia of Mathematics[END_REF][START_REF] Guy | Mersenne Primes. Repunits. Fermat Numbers. Primes of Shape[END_REF]. Therefore, one can expect that all the previously folded polygons (n = 3, n = 2 2 = 4, n = 5, n = 2 1 × 3 = 6, n = 2 3 = 8) can also be built with straight-edge and compass. This won't be the case for heptagon (n = 7 edges) and nonagon (n = 9 edges), for instance. Using the technique mentioned in the introduction to solve third-order equations by folding, the previous set of constructible polygons can be extended to the set of polygons with n = 2 m 3 q g 1 g 2 ...g s edges, where the g i are all different primes, of the form 2 p 3 r + 1, m, p, q, and r being integers (see for instance [START_REF] Scimeni | Draw of a regular heptagon by the folding[END_REF][START_REF] Justin | Résolution par le pliage de l'équation du troisième degré et applications géométriques[END_REF][START_REF] Alperin | A mathematical theory of origami constructions and numbers[END_REF]). Such a construction will be necessary to fold the optimal heptagon and nonagon.

To introduce this new construction, let us recall one of its applications: the trisection of an arbitrary angle θ [START_REF] Fusimi | Trisection of angle by Hisashi Abe[END_REF], known to be impossible with Euclidean constructions (for a nice discussion on this subject, the reader may refer to the web page http://www.jimloy.com/geometry/trisect.htm). To trisect an arbitrary angle, consider the construction of Figure 8: an isoceles triangle AA B has γ as half main angle. Build the perpendicular AD to AB to get ∠A AD = γ. The main idea is then to look at the axis of symmetry ∆ that The trisection of an angle requires to solve a third-order polynomial equation, while intersecting circles and lines (basically, the Euclidean constructions) requires only second-order polynomial equations. Equipped with the operation of Figure 9, the paperfolder is now able to construct additional polygons.

Since folding a nonagon (n = 9 edges) is a little easier than folding the heptagon, let us begin with it. Figure 10 describes shortly the corresponding sequence. Folding it requires precision (and a large square of paper), and proving that an exact optimal nonagon is obtained is a not so easy task [START_REF] Dureisseix | Searching for optimal polygon -remarks about a general construction and application to heptagon and nonagon[END_REF]. Both are left as challenges to the reader. Here are some guidelines: Here, we have a 9/2 regular star polygon).

• Step 4: Report this distance at JK.

•

Step 5: Complete the nonagon (quite challenging, again).

An even more challenging construction is the one of the heptagon (n = 7 edges). This time, the angle 2π/7 is involved, still anattainable with Euclidean constructions. Therefore, it requires also the previous technique, designed to solve any third-order polynomial root [START_REF] Justin | Résolution par le pliage de l'équation du troisième degré et applications géométriques[END_REF] (let us just mention that 2 cos(2π/7) is a root of t 3 + t 2 -2t -1 = 0, and the solving technique is the same as in [START_REF] Huzita | Drawing the regular heptagon and the regular nonagon by origami (paper folding)[END_REF][START_REF] Geretschläger | Folding the regular heptagon[END_REF]). Figure 11 describes shortly the corresponding sequence:

• Step 1: Precrease diagonal and central lines. Fold in half and unfold.

• Step 2: Use previous technique, to get a fold HI of slope 2 cos(2π/7) = AI/AH (G is located at the middle of the previous two horizontal folds).

• Step 3: Fold in half and unfold.

•

Prospects

Several regular polygons can be folded with Euclidean constructions. Their optimal versions can be folded also, though practically, they become less and less easy to obtain. With a special basic fold that is recalled in this paper, and that is now widely spread through the community of paperfolders, more can be done. But even with it, not all of the polygons can be folded: for instance, the first unreachable regular polygon is the hendecagon (n = 11 edges). A construction, simple enough to enter the standard repertoire, and allowing the construction of the regular hendecagon, would require the solution of a 5th-order polynomial equation (the trigonometric functions of 2π/11 are roots of such a polynomial equation). It is still to come.

Figure 3 :

 3 Figure 3: Strip problem solving procedure

Figure 4 :Figure 5 :

 45 Figure 4: Placement of optimal polygons x y z A B

• Step 2 :

 2 Turn over the model. Split DF in two by folding, and create the fold GH. Its length is GH = DF and due to symmetry with respect to

Figure 6 :

 6 Figure 6: Optimal pentagon and its stellated version (or pentagram)

Figure 8 :Figure 9 :

 89 Figure 8: Illustration of trisection problem y x

• Step 1 :• Step 2 :

 12 Precrease diagonal and central lines. Fold π/3 = ∠BAE; D should lie on the central horizontal line. Unfold. With the trisection method, fold 2π/9 = ∠BAF (then, ∠F AC = π/4-2π/9 = π/36). The point G is the intersection of AF and the central vertical line. Since AH = 1/2, AG = 1/(2 cos(2π/9)).

Figure 10 :

 10 Figure 10: Folding sequence of the optimal nonagon

Step 4 :• Step 5 :Figure 11 :

 4511 Figure 11: Folding sequence of the optimal heptagon

Acknowledgements. The author wishes to thank the generous referees for their help and fruitful discussions about the topic of this note.