Michel Bergmann
email: michel.bergmann@inria.fr

Lisl Weynans
email: lisl.weynans@math.u-bordeaux.fr

A Cartesian Method with Second-Order Pressure Resolution for Incompressible Flows with Large Density Ratios

Keywords: incompressible Navier-Stokes equations, projection method, finite differences, Cartesian grid, immersed interfaces, level-set, interface unknowns

An Eulerian method to numerically solve incompressible bifluid problems with high density ratio is presented. This method can be considered as an improvement of the Ghost Fluid method, with the specificity of a sharp second-order numerical scheme for the spatial resolution of the discontinuous elliptic problem for the pressure. The Navier-Stokes equations are integrated in time with a fractional step method based on the Chorin scheme and discretized in space on a Cartesian mesh. The bifluid interface is implicitly represented using a level-set function. The advantage of this method is its simplicity to implement in a standard monofluid Navier-Stokes solver while being more accurate and conservative than other simple classical bifluid methods. The numerical tests highlight the improvements obtained with this sharp method compared to the reference standard first-order methods.

Introduction

Bifluid problems are ubiquitous in nature and in many industrial applications such as combustion in engines, water waves energy converters, and jet printers to cite only a few. In such applications, the density ratio between the two fluids can be large, for instance the ratio is equal to 1000 between water and air. Accurate numerical modeling and numerical simulations of these kind of phenomena are then necessary, in particular to optimize such devices.

In this paper we are thus concerned with the numerical modeling of incompressible bifluid flows with large density ratios, like air and water, and by the accurate description of the phenomena occurring at their interface. We present a sharp Cartesian method for the simulation of incompressible flows with high density and viscosity ratios. This method is an extension of the second-order Cartesian method for elliptic problems with immersed interfaces developed in [START_REF] Cisternino | A parallel second order Cartesian method for elliptic interface problems[END_REF].

Cartesian grids are an attractive alternative to body fitted meshes. Indeed, they avoid complex mesh generation as well as mesh adaptation when unsteady interfaces are considered. Moreover, the numerical resolution of the governing equations can be simplified with an easy parallelization and the use of standard linear algebra libraries. Generally speaking, numerical schemes are easy to implement on a Cartesian mesh because a dimensional splitting is often possible; however, some numerical modeling is necessary near a complex interface that does not fit the background Cartesian grid. This is the case for fluid structure interface and moreover for bifluid interface where the properties of the flow are discontinuous. Indeed, applying naively a numerical scheme originally devised for a flow with constant or continuously varying density will lead to a non-consistent treatment of the interface. Most of the time, it will result in severe stability issues if the density ratio is large as highlighted in [START_REF] Raessi | Consistent mass and momentum transport for simulating incompressible interfacial flows with large density ratios using the level set method[END_REF] and references therein; therefore, as already mentioned, one has to devise specific numerical schemes at the vicinity of the interface. This region is called narrow band and is the set of numerical points that have at least one neighbor on the other side of the interface.

Conservative or non-conservative approaches can both be used to face this issue. Among the non-conservative approaches, one solution is to regularize the properties of the fluids in the vicinity of the interface, so that the density, viscosity, and their derivatives are continuous in the whole computational domain. This idea leads to the well-known "Continuous Surface Force" (CSF) method [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF], where the discontinuous quantities are smoothed near the interface, and in case of a fluid with surface tension, this surface tension is taken into account as a smooth volume force instead of a surface force. This method is widely used (see for instance [START_REF] Sussman | A level-set approach for computing solutions to incompressible two-phase flows[END_REF][START_REF] Galusinski | On stability condition for bifluid flows with surface tension: Application to microfluidics[END_REF]) because it offers a straightforward way to implement the presence of two fluids in an already existing monofluid Navier-Stokes code. However, the exact way that the regularization should be performed is not always clear, and spurious oscillations at the bifluid interface can appear due to errors in the pressure gradient computations. Another non-conservative method introduced by Kang, Fedkiw, and Liu [START_REF] Kang | A Boundary Condition Capturing Method for Multiphase Incompressible Flow[END_REF] after the CSF is the Ghost Fluid Method (GFM). It is based on a firstorder method developed in [START_REF] Liu | A boundary capturing method for Poisson's equation on irregular domains[END_REF] to solve an immersed interface elliptic problem, with a dimensional splitting making the method easy to implement. The resulting linear system is symmetric and has the same structure as the usual matrix to discretize a Poisson equation with variable coefficient on a Cartesian grid. This method has been used successfully in a lot of works, for instance [START_REF] Couderc | Développement d'un code de calcul pour la simulation d'écoulements de fluides non miscibles: Application à la désintégration assistée d'un jet liquide par un courant gazeux[END_REF][START_REF] Tanguy | A Level Set Method for vaporizing two-phase flows[END_REF]. One drawback is that the method is only first-order accurate near the interface [START_REF] Raessi | Consistent mass and momentum transport for simulating incompressible interfacial flows with large density ratios using the level set method[END_REF] and a loss of conservativity of the momentum of each fluid near the interface can occur leading to erroneous velocities.

Non-conservative methods are often associated with a level-set representation of the interface [START_REF] Sussman | A level-set approach for computing solutions to incompressible two-phase flows[END_REF] because the level-set method is itself intrinsically non-conservative at the discrete level, and convenient to use on a Cartesian grid.

The other family of methods is based on the conservative form of the Navier-Stokes equations, where mass and momentum fluxes of each fluid are explicitly computed, see for instance [START_REF] Raessi | Consistent mass and momentum transport for simulating incompressible interfacial flows with large density ratios using the level set method[END_REF][START_REF] Rudman | A volume-tracking method for computing incompressible multifluid flows with large density variations[END_REF][START_REF] Sussman | A coupled level set and volume of fluid method for computing 3D and axisymmetric incompressible two-phase flows[END_REF][START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF][START_REF] Hu | A conservative interface method for compressible flows[END_REF]]. An explicit interface representation is necessary even if the interface do not coincide with grid points. Conservative methods are generally more stable than non-conservative methods. The price for this increased stability is an additional amount of numerical developments due to the interface reconstruction, which can be performed from information carried by Lagrangian markers or by cell quantities such as volume fractions.

Another approach has been developed recently [START_REF] Nangia | A robust incompressible Navier-Stokes solver for high density ratio multiphase flows[END_REF][START_REF] Nangia | A DLM immersed boundary method based wave-structure interaction solver for high density ratio multiphase flows[END_REF] to deal with large density ratio. In this approach, a fully second-order method is obtained at the interface (for both velocity and pressure) with several physical boundary treatments, including velocity and traction boundary conditions.

In this paper we aim to preserve as much as possible the simplicity of the Ghost Fluid method of [START_REF] Kang | A Boundary Condition Capturing Method for Multiphase Incompressible Flow[END_REF], avoiding an explicit identification of the volume fractions near the interface, while improving the accuracy and stability of the pressure computation. We thus propose a method, mainly based on the improvement of the Ghost Fluid method, where the discontinuities across the interface are taken into account in a sharp way with a secondorder scheme inspired from [START_REF] Cisternino | A parallel second order Cartesian method for elliptic interface problems[END_REF]. This second-order treatment improves the conservativity of the method, as it will be proved numerically in the section devoted to numerical validations.

After having described the governing equations for the incompressible bifluid flows that we consider (Section 2), the discretization of these equations in each fluid and at the interface are presented (Sections 3 and 4). The second-order numerical resolution of the elliptic problem arising from the computation of the pressure is introduced (Section 5), and the overall is validated on several two-dimensional numerical test cases (Section 6).

Governing Equations

Flow Equations

We consider a rectangular domain Ω filled with two viscous incompressible fluids with different densities and viscosities. The subdomains Ω -and Ω + corresponding to the two fluids are separated by an interface Γ as depicted in Figure 1.

In this work, these domains are implicitly defined with a scalar function φ, usually called the level-set function, see Section 2.2, that takes different values in each subdomain with a fixed value on the interface. For instance we chose φ = 0 on Γ, φ > 0 in Ω + and φ < 0 in Ω -. The unit normal to the interface is denoted n and the unit tangent vector is denoted η. The density is denoted

ρ = ρ -+ H(φ)(ρ + -ρ -), (1)
and the viscosity is denoted

µ = µ -+ H(φ)(µ + -µ -), ( 2
)
where H is the Heavyside function, i.e., H(φ) = 1 if φ > 0 and H(φ) = 0 if φ < 0. Finally, the two-dimensional velocity vector is denoted u = (u, v).

The flow is modeled in the whole domain with the incompressible Navier-Stokes equations:

ρ(u t + (u • ∇)u) = -∇p + ∇ • τ + ρg -σκ∇H, ∇ • u = 0,
with g the gravitational acceleration vector, τ the viscous stress tensor:

τ = µ(∇u + ∇u T ), (3)
and the term σκ∇H accounting for the surface tension effects, with σ the surface tension itself and κ the local curvature of the interface between the fluids. This formulation of the bifluid incompressible Navier-Stokes equations contains a singular term which is not trivial to handle.

Ω -(φ < 0) ρ -, µ - Ω + (φ > 0) ρ + , µ + Γ (φ = 0) n Figure 1.
Sketch of the computational domain.

Alternatively, the flow can also be modeled in each subdomain with the incompressible Navier-Stokes equations:

ρ(u t + (u • ∇)u) = -∇p + ∇ • τ + ρg, ∇ • u = 0.
The above equations are completed by jump conditions at the interface Γ between the two fluids. In what follows, jumps are defined by the notation [ψ] = ψ +ψ -. The first jump conditions describe the balance between the normal stresses at the interface and the surface tension σ, with κ the local curvature of the interface Γ,

[p -2µ(∇u • n, ∇v • n) • n] = σκ, (4)
[µ(∇u

• n, ∇v • n) • η + (∇u • η, ∇v • η) • n] = 0. ( 5
)
Others jump conditions can be derived from continuity properties across the interface. For instance, for a viscous fluid, the velocity field is continuous across the interface

[u] = 0, (6)
[v] = 0.

Since the material derivative of ( 6) and ( 7) is zero, one can write

0 = ∂[u] ∂t + (u • ∇)[u] = - ∇p ρ + ∇ • τ ρ + g , (8)
which leads to ∇p

ρ = ∇ • τ ρ . (9)
The jump condition for the pressure p can be simplified. We differentiate the jump on the velocity in the tangential direction:

∂u ∂η = 0, ∂v ∂η = 0.
Moreover, because the velocity is divergent-free on each side of the interface,

0 = [∇ • u] = [(∇u • n, ∇v • n) • n + (∇u • η, ∇v • η) • η]. ( 10
)
Combining the two last relationships, we obtain

[(∇u • n, ∇v • n) • n] = 0. (11)
Consequently

[p] = 2[µ](∇u • n, ∇v • n) • n + σκ. ( 12
)
Finally, we will use Equations ( 9) and [START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF] to compute the pressure jump at the interface.

Interface Description

The interface between the two fluid subdomains is implicitly defined by a scalar function φ. Local geometrical information on the function φ is needed to obtain an accurate discretization of the interface. To this purpose we use the level set method, introduced by Osher and Sethian [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamiltonâ Jacobi formulations[END_REF] and described in [START_REF] Sethian | Level Set Methods and Fast Marching Methods[END_REF][START_REF] Sethian | Implementation, and Application of Level Set and Fast Marching Methods for Advancing Fronts[END_REF][START_REF] Osher | Level Set Methods and Dynamic Implicit Surfaces[END_REF]. A common choice for the level-set function φ is the signed distance function to the interface:

φ(x) =    dist Γ (x) if x ∈ Ω + , -dist Γ (x) if x ∈ Ω -, 0 if x ∈ Γ. (13)
The zero isoline of φ thus represents implicitly the interface Γ immersed in the computational domain.

We assume that the interface is smooth enough, so that the derivatives of the level-set function in the vicinity of the interface are well-defined. A useful property of the level-set function is a straightforward computation of its normal with the formula

n(x) = ∇φ(x) |∇φ(x)| . ( 14
)
In the same way, the curvature of the interface can be computed with the formula

κ = ∇ • n. (15)
We impose a curvature threshold 1/h, with h the grid spacing, corresponding to the minimum size of a bubble for a given spatial discretization.

For a moving interface as it is usually the case for bifluid problems, the flow density and viscosity are updated with φ tracking the interface thanks to the transport equation

φ t + û • ∇φ = 0, (16)
where the velocity fields û coincides with the flow velocity field u on the interface Γ.

Different choices for the value of û in Ω + and Ω -can be a priori used. A natural choice we have considered is û = u in the whole domain:

φ t + u • ∇φ = 0. ( 17
)
Another possible choice is the extension velocity introduced in [START_REF] Adalsteinsson | The Fast Construction of Extension Velocities in Level Set Methods[END_REF].

Navier-Stokes Monofluid Solver and Numerical Method for Interface Evolution

The computational domain is discretized on a two-dimensional uniform Cartesian grid with a grid spacing ∆x = ∆y = h. However, the following approach stands for nonuniform Cartesian meshes. The points on the Cartesian grid are defined as M i,j = (x i , y j ). Similarly, we denote by u ij the approximation of u at the point (x i , y j ). In what follows, all the unknowns are collocated in space on Cartesian meshes. An odd-even coupling can sometimes be observed when one uses collocated unknowns on a Cartesian grid. This problem can be fixed using some corrections such as [START_REF] Rhie | Numerical study of the turbulent flow past an airfoil with trailing edge separation[END_REF][START_REF] Mittal | A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries[END_REF]. However, in all the applications considered in this paper we have not observed any odd-even coupling.

In this section we present the numerical solver devised for monofluid incompressible Navier-Stokes equations that is used in each fluid when the interface is not taken into account. We also provide information about the numerical methods used to compute the evolution of the level-set function.

Flow Computation

We use a classical projection method [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF][START_REF] Temam | Sur l'approximation de la solution des equations de Navier-Stokes par la méthode des pas fractionnaires II[END_REF] to solve the Navier-Stokes equations in each fluid outside the narrow band. In what follows, a non-incremental projection method is used, i.e., the guess value for the pressure in the prediction step is zero.

We thus compute successively:

u * -u n ∆t = -(u n • ∇)u n + 1 ρ (∇ • τ) n + g (prediction step), (18)
u n+1 -u * ∆t = - ∇p ρ (correction step). ( 19
)
The convective terms are computed with a fifth-order WENO scheme, and the viscous terms with an explicit second-order centered finite-difference scheme. The time integration is performed with a first-order explicit Euler scheme, which is consistant with the use of a non-incremental version of the projection scheme.

The pressure appearing in Equation ( 19) is computed through the resolution of a Poisson equation in order to enforce the divergence-free condition. At each point in one subdomain, the following relationship is satisfied:

∇ • 1 ρ ∇p = ∇ • u * ∆t . ( 20
)
We will provide additional details about the jump conditions that have to be satisfied across the interface for this problem in Section 4. On the exterior boundary of the domain Ω, Neumann boundary conditions are satisfied:

∇p ρ = u b -u * ∆t , (21)
where u b is the value of the velocity to be imposed on the external boundaries. We compute at each iteration an adaptive time step taking into account the restrictions due to convection, viscosity, and surface tension. The convective time step restriction is given by

∆t |u| max ∆x + |v| max ∆y ≤ 1. ( 22
)
with |u| max and |v| max the maximum magnitudes of the horizontal and vertical velocities. The viscous time step restriction is given by

∆t max µ - ρ -, µ + ρ + 2 ∆x 2 + 2 ∆y 2 ≤ 1 (23)
We also apply a time step restriction associated with the surface tension evaluated only in the narrow band. This time step restriction is similar to the one in [START_REF] Kang | A Boundary Condition Capturing Method for Multiphase Incompressible Flow[END_REF] and in [START_REF] Couderc | Développement d'un code de calcul pour la simulation d'écoulements de fluides non miscibles: Application à la désintégration assistée d'un jet liquide par un courant gazeux[END_REF], and is in this context usually the most restrictive:

∆t σ|κ| min(ρ + , ρ -) min(∆x 2 , ∆y 2 ) ≤ 1. ( 24
)
Finally, at each time step, the overall algorithm is the following:

• Prediction: evaluate convective and diffusive fluxes and compute u * , •

Interface evolution: convect the level-set with velocity u and re-initialize if necessary,

•

Construction and resolution of the linear system for the pressure, •

Correction step: update velocity with pressure gradient.

Numerical Method for the Level-Set Evolution

The computation of the level-set function should be performed very accurately when one deals with moving interfaces. Indeed, as the level-set method is not intrinsically conservative, a lack of accuracy in the computation of the level set evolution results often in a substantial loss of mass for one of the fluids. It can also increase the problem of transfer of momentum between both fluids and generate spurious velocity oscillations. Moreover, if one wants to compute the curvature of the interface from ( 15), the level-set function needs to be accurate enough (at least third-order) so that the finite difference formulas used to discretize (15) are consistent.

Unfortunately, the property of the signed distance function is usually lost when the interface evolves with the flow velocity. The norm of the level set gradient can be far from unity. These gradients variations of the level-set are harmful to the accuracy of the numerical evaluation of the normal to the interface and the curvature. To circumvent the problem, Sussman et al. [START_REF] Sussman | A level-set approach for computing solutions to incompressible two-phase flows[END_REF] introduced a reinitialization algorithm to recover the signed distance function through the resolution of the eikonal equation

|∇φ| = 1. ( 25
)
Several methods have been developed over the years to perform this reinitialization step, either by using a relaxation method and searching a stationary solution to a time dependent Hamilton-Jacobi equation [START_REF] Sussman | A level-set approach for computing solutions to incompressible two-phase flows[END_REF] or using a Gauss-Seidel-based method as in fast marching methods [START_REF] Sethian | A fast marching level set method for monotonically advancing fronts[END_REF][START_REF] Rouy | A viscosity solutions approach to shape-from-shading[END_REF] or fast sweeping methods [START_REF] Tsai | Fast sweeping algorithms for a class of Hamilton-Jacobi equations[END_REF].

The usual level set strategy for an evolving interface is thus the following:

• Use a transport equation to update φ.

•

From time to time, reinitialize φ with the signed distance function.

One of the most widespread option in the literature is to combine a fifth-order WENO scheme [START_REF] Jiang | Efficient Implementation of Weighted ENO Schemes[END_REF] with a RK3 scheme, for the transport and for the reinitialization through a relaxation equation. It provides a high-order yet stable resolution; however, although the reinitialization procedure, performed with such a numerical scheme, may improve mass conservation, it also introduces some error by slightly moving the interface as shown in [START_REF] Russo | A remark on computing distance functions[END_REF]. For this reason, Russo and Smereka [START_REF] Russo | A remark on computing distance functions[END_REF] introduced a subcell fix taking into account the interface location in the reinitialization procedure. This technique was extended to a higher order accuracy in [START_REF] Duchene | Second-order accurate computation of curvatures in a level set framework using novel high-order reinitialization schemes[END_REF] through the use of third-order ENO schemes near the interface.

Usually, this reinitialization steps are performed uniformly every each n iterations. A recent study [START_REF] Luddens | Enablers for high-order level set methods in fluid mechanics[END_REF] proposes a strategy to sample the reinitialization steps based on interface deformation criteria. This should also be coupled with the high-order decentered reinitialization scheme of [START_REF] Duchene | Second-order accurate computation of curvatures in a level set framework using novel high-order reinitialization schemes[END_REF] near the interface.

In what follows, we use the classical option of the fifth-order WENO scheme for the spatial discretization. Since it is the most commonly used technique in the literature, it will allow us to distinguish the effects of the new scheme for the pressure computation from the effects of the reinitialization technique.

Navier-Stokes Solver near the Interface

Notations

Let us introduce some definitions and notations. A grid point is defined to be irregular if at least one of its neighbors is on the other side of the interface, i.e., if the sign of φ changes between this point and at least one of its neighbors, see Figure 2. The set of irregular points is called the narrow band. All the other points are called regular grid points.

We define the interface point I i,j,E = ( xi,j,E , y j ) as the intersection of the interface Γ and the East segment [M ij M i+1j ], if it exists. Similarly, the interface points I i,j,W = ( xi,j,W , y j ), I i,j,N = (x i , ỹi,j,N ) and I i,j,S = (x i , ỹi,j,S ) are defined as the intersection of the interface and the

West [M i-1j M ij ], North [M ij M ij+1 ] and South [M ij-1 M ij
] segments, respectively. With this notation the same interface point can be described in two different ways

I i,j,S = I i,j-1,N or I i,j,E = I i+1,j,W . ( 26
)
The set of interface points is denoted Γ h , see Figure 2 for an illustration.

• • • • • • • • • • • • • • • • • • • • • • • • • • I i,j,E = I i+1,j,W I i+2,j+1,N = I i+2,j+2,S i-1 i i+1 i+2 i+3 i+4 i+5 j-1 j j+1 j+2 j+3 • • • • • • • • • Figure 2.
Example of geometrical configuration, with regular grid points (white circles), irregular points (black circles), and interface points (diamonds) with the two possible notations.

Modeling Choices for the Discontinuities across the Interface

In this study, the values of the viscosity and the densities are discontinuous across the interface; therefore, if the numerical scheme for solving incompressible Navier-Stokes equations described previously in Section 3.1 was applied on the irregular grid points of the narrow band, the approximations of the following terms would not be consistent:

•

Prediction step: viscous terms.

•

Correction step: divergence of the predicted velocity, elliptic operator, and gradient of the pressure.

This lack of consistency could eventually lead to stability problems. The computation of the convective terms is not mentioned in the above enumeration because it is performed with a fifth-order WENO scheme, which automatically provides spatial adaptivity. Therefore, we assume that the gradients computed with the WENO scheme are decentered near the interface, and consequently, consistent. Moreover, the levelset function is classically evolved with such a scheme, because it is crucial to have a good accuracy in the computation of the interface evolution; thus, it seems coherent to have the same numerical scheme for the convection of the interface and the convection of the fluids.

In this work, we use two different strategies to handle the lack of consistency near the interface, one for the viscous terms and another for the pressure computation. For the viscous term in the prediction step, we follow a continuous approach and regularize the quantities used for the computation of the viscous terms. It has been proven in [START_REF] Desjardins | A spectrally refined interface approach for simulating multiphase flows[END_REF][START_REF] Herrmann | The influence of density ratio on the primary atomization of a turbulent jet in crossflow[END_REF] that this continuous approach provides correct accuracy for high Reynolds number flows. It has also been used successfully in [START_REF] Raessi | Consistent mass and momentum transport for simulating incompressible interfacial flows with large density ratios using the level set method[END_REF][START_REF] Brackbill | A continuum method for modeling surface tension[END_REF]. A sharp approach for the viscous terms could probably improve the accuracy of the simulations; however, the complexity of the computations would be increased due to the treatment of the jump conditions for the viscous terms (5) implying derivatives of the velocity components in both normal and tangential directions. Moreover, if one needs to use an implicit treatment of the viscous terms, such a sharp treatment would become more complex to handle.

The discretization of the prediction step becomes:

u * -u n ∆t = -(u n • ∇)u n + 1 ρ ∇ • μ(∇u + ∇u T ) n + g. (27)
In practice the viscosity and the inverse of the density were regularized by a discrete convolution [START_REF] Raessi | Consistent mass and momentum transport for simulating incompressible interfacial flows with large density ratios using the level set method[END_REF]:

16 μi,j = 4 µ i,j + 2 µ i+1,j + 2 µ i-1,j + 2 µ i,j+1 + 2 µ i,j-1 +µ i+1,j+1 + µ i+1,j-1 + µ i-1,j+1 + µ i-1,j-1 , 16 ρi,j = 4 ρ i,j + 2 ρ i+1,j + 2 ρ i-1,j + 2 ρ i,j+1 + 2 ρ i,j-1 + 1 ρ i+1,j+1 + 1 ρ i+1,j-1 + 1 ρ i-1,j+1 + 1 ρ i-1,j-1 .
Then, we discretize the viscous terms with a classical second-order centered scheme as in Section 3.1.

In the correction step, according to the jump conditions ( 4)-( 12) presented in Section 2, the pressure satisfies an elliptic problem with discontinuous values of the solution and its derivative across the interface:

∇ • 1 ρ ∇p = ∇ • u * ∆t in Ω + ∪ Ω -, [p] = σ κ + 2[µ] (u n , v n ) • n on Γ, ∇p ρ = ∇ • τ ρ on Γ.
Because the viscous terms in the prediction step are handled with a regularization approach, we have [µ] = 0 and ∇•τ ρ = 0. Therefore, the pressure computed for the correction step satisfies rather:

∇ • 1 ρ ∇p = ∇ • u * ∆t in Ω + ∪ Ω -, (28)
[p] = σ κ on Γ,

∇p ρ = 0 on Γ. ( (29)
) 30
The details of the resolution of this elliptic problem will be provided in Section 5. Let us remark however that as the values of p are discontinuous across the interface, two values of p will be created at each point of the interface considered in the numerical scheme.

Gradient and Divergence for Correction Step

The predicted velocity u * obtained after the prediction step ( 18) is defined only on grid points. We need to compute the divergence of this predicted velocity in order to solve the elliptic Equation [START_REF] Adalsteinsson | The Fast Construction of Extension Velocities in Level Set Methods[END_REF]. However, since the two fluids have different properties across the interface and the derivatives of the velocity are not necessarily continuous, we need to use a decentered stencil on each side of the interface. Consequently, we have to compute two values for u * on each interface point, one for each side of the interface. In practice, as jump conditions for u * are not available, we perform simply linear extrapolations from the grid values on the interface points. Then, to compute the divergence of u * on an irregular grid point M i,j , we use a standard five point stencil, see Figure 3. Formally this is equivalent to a standard first-order decentered scheme.

More precisely, we denote u * S the value of the predicted velocity u * on the nearest point to M i,j in the south direction (possibly an interface point), with coordinates (x S , y S ). Similarly, we define u * N , u * W and u * E and the associated coordinates (x N , y N ), (x W , y W ) and (x E , y E ). The discretization reads

∇ • u * i,j = u * E -u * W x E -x W + v * N -v * S y N -y S .
Similarly, in order to keep a consistent discretization, the gradient of the pressure p appearing in Equation ( 19) is also computed with an adapted decentered stencil near the interface, see Figure 3. More precisely, with the same notations S, N, W and E, as before, the discretization reads

(∇p) i,j =    p E -p W x E -x W p N -p S y N -y S   . ( 31
)
If one of the discretization point is an interface point, we consider the value of the pressure on this point corresponding to the same subdomain than point M i,j . Indeed we recall that, since the pressure is discontinuous across the interface, two pressure unknowns (one for each subdomain) are computed at an interface point. If no interface point is involved, the numerical scheme reduces to the classical second-order central finite differences scheme.

• • • • • • M i+2,j+2 • • • • • M i,j i-1 i i+1 i+2 i+3 i+4 i+5 j-1 j j+1 j+2 j+3 Figure 3.
Example of geometrical configuration, with points involved in the discretization of the divergence of the predicted velocity and pressure gradient at grid points M i,j and M i+2,j+2 in black.

Numerical Resolution of Elliptic Problems with Immersed Interfaces

The elliptic problem with discontinuous values across an interface ( 28)-( 30) is solved with the second-order method developed in [START_REF] Cisternino | A parallel second order Cartesian method for elliptic interface problems[END_REF]. The interface points are used to impose the jump conditions across the interface in the numerical scheme. Since the pressure is discontinuous across the interface, two unknowns are created, one for each side of the interface. The accuracy of this method is based on the use of unknowns located at the interface. The size of this linear system is thus augmented with two unknowns for each interface point. These interface unknowns are used to discretize the flux jump conditions and the elliptic operator accurately enough to get a second-order convergence in maximum norm. Actually, to this purpose, near the interface the elliptic operator needs to be discretized with a first-order truncation error, and the fluxes with a second-order truncation error. For a visual explanation of the discretization we refer to Figure 4. The advantage of using this method, compared to the reference work of [START_REF] Kang | A Boundary Condition Capturing Method for Multiphase Incompressible Flow[END_REF] is that the jump conditions in the correction step are solved with second-order accuracy instead of firstorder. The drawback is that the linear system is not symmetric anymore and it is solved with the preconditioned GMRES method.

Discrete Elliptic Operator

We use a standard five-point stencil including the grid point M i,j and its nearest neighbors in each direction: interface or grid points. More precisely, we denote p S the value of the solution on the nearest point in the south direction, with coordinates (x S , y S ).

Similarly, we define p N , p W , and p E and the associated coordinates (x N , y N ), (x W , y W ), and (x E , y E ). Since the density is piecewise constant, the discretization reads

∇.( 1
ρ ± ∇p) i,j = 1 ρ ± ∆p = 1 ρ ± p E -p ij x E -x i - p ij -p W x i -x W x E -x W 2 + 1 ρ ± p N -p ij y N -y j - p ij -p S y j -y S y N -y S 2 , ( 32
)
where ρ ± stands for ρ + or ρ -. Example of stencils for the discretization of jump conditions. Points involved in the discretization of the x-derivative of the pressure at interface point I i,j-1,E and of the y-derivative of the pressure at interface point I i+1,j+1,N in black. For I i,j-1,E both derivatives are expressed with second-order accuracy while for I i+1,j+1,N the left derivative is expressed with second-order and the right derivative with first-order accuracy.

• • • • M i+2,j+2 • M i+3,j • • • • • i-1 i i+1 i+2 i+3 i+4 i+5 j-1 j j+1 j+2 j+3 • • • • • I i+1,j+1,N • •• • I i,j-1,E i-1 i i+1 i+2 i+3 j-1 j j+1 j+2

Discrete Flux Transmission Conditions

As written before, at each interface point we create two additional unknowns, called interface unknowns. We denote them by p ± i,j,γ with γ = E, W, N or S. The interface unknowns carry the values of pressure on each side of the interface.

Contrarily to [START_REF] Cisternino | A parallel second order Cartesian method for elliptic interface problems[END_REF], we do not have a jump condition on the normal derivative, but on the whole gradient, as expressed in formula [START_REF] Duchene | Second-order accurate computation of curvatures in a level set framework using novel high-order reinitialization schemes[END_REF]. To obtain the same number of equations and unknowns we have to chose in which direction we want to project this gradient equality. As we use a Cartesian grid, it is easier to discretize the xand y-derivatives than derivatives in other directions. Consequently, we discretize the following jump conditions at each interface point I i,j,γ , with γ = N, S, W, E.

p + i,j,γ -p - i,j,γ = σ κ i,j,γ , (33)
1 ρ + (∂ x p + ) i,j,γ - 1 ρ -(∂ x p -) i,j,γ = 0 if γ = E, W. (34)
1

ρ + (∂ y p + ) i,j,γ - 1 ρ -(∂ y p -) i,j,γ = 0 if γ = N, S. (35)
We want the truncation error of the discretization of flux equality [START_REF] Chenadec | A monotonicity preserving conservative sharp interface flow solver for high density ratio two-phase flows[END_REF] and ( 35) to be second-order accurate in order to solve the problem with a second-order accuracy. A possible configuration of the interface is illustrated in Figure 4. In the x-direction, it is straightforward to compute a second-order approximation of the x-derivative with three a priori nonequidistant points. For example we approximate the flux on the left side of interface point I i,j,E , if it exists, with the values of p on the points M i-1,j , M i,j and I i,j,E with the formula:

(∂ x p ± ) i,j,E ≈ (p i-1,j -p ± i,j,E )(x i -xi,j,E ) ∆x(x i-1 -xi,j,E ) - (p i,j -p ± i,,Ej )(x i-1 -xi,j,E ) ∆x(x i -xi,j,E ) . ( 36
)
The right x-derivative is approximated in the same way.

(∂ x p ± ) i,j,E ≈ - (p i+2,j -p ± i,j,E )(x i+1 -xi,j,E ) ∆x(x i+2 -xi,j,E ) + (p i+1,j -p ± i,,Ej )(x i+2 -xi,j,E ) ∆x(x i+1 -xi,j,E ) . ( 37
)
The same discretization holds for the y-derivative. The formulas ( 36) and ( 37) are consistent if both grid points involved in the formula, for instance M i-1,j and M i,j , belong to the same domain. If on one side of the interface the two closest grid points aligned with the intersection point do not belong to the same subdomain, then the second-order discretization is not possible anymore. In this case, we use a first-order discretization involving only two points: the interface point and the closest grid point on the same side of the interface. Such a case is illustrated on Figure 4. In fact, this first-order discretization is equivalent to the Ghost Fluid method [START_REF] Kang | A Boundary Condition Capturing Method for Multiphase Incompressible Flow[END_REF].

Let us notice that, because we use a dimensional splitting for the jump conditions across the interface, it is quite straightforward to eliminate the interface unknowns from the linear system. We simply inject expressions [START_REF] Grenier | An Hamiltonian interface SPH formulation for multi-fluid and free-surface flows[END_REF] and [START_REF] Khedkar | The inertial sea wave energy converter (ISWEC) technology: Device-physics, multiphase modeling and simulations[END_REF] in the jump condition [START_REF] Chenadec | A monotonicity preserving conservative sharp interface flow solver for high density ratio two-phase flows[END_REF], and use the resulting equality to express p ± i,j,E as a function of p i-1,j , p i,j , p i+1,j and p i+2,j . This expression for p ± i,j,E can then be used in the discretization of the elliptic operator (32). The local curvature κ i,j,γ at the interface point I i,j,γ is computed in the following way. We first compute on all irregular grid points the value κ = φ 2

x φ yy + φ 2 y φ xx -2φ x φ y φ xy (φ 2

x + φ 2 y ) 3/2 [START_REF] Bhalla | Simulating water-entry/exit problems using Eulerian-Lagrangian and fully-Eulerian fictitious domain methods within the open-source IBAMR library[END_REF] with centered second-order finite-difference formulas. Then we perform a one-dimensional linear interpolation of these values on the interface points.

Numerical Results and Validations

This section is devoted to the numerical simulations and validations.

After having studied the well-balanced behavior of the proposed approach, i.e., the capacity of the method to preserve equilibria state for a bubble, in Section 6.1, we perform a quantitative validation for the well-known dam break problem from Section 6.2. We then present other numerical simulations for the rising of bubbles with different sizes that can qualitatively been compared to some reference results in Section 6.3.

Equilibria Preservation for a Bubble: the Parasitic Oscillations

This first test case aims to assess the influence of the interface curvature error on the stability of the numerical scheme. A bubble is located at the center of the computational domain in Sections 6.1.1 and 6.1.2. Due to the Laplace law and the concavity of the interface, the pressure inside the bubble is larger than the pressure outside. If the curvature of the interface is computed numerically, the errors due to the numerical approximation in the right-hand side of Equation [START_REF] Herrmann | The influence of density ratio on the primary atomization of a turbulent jet in crossflow[END_REF] cause small errors in the resolution of the pressure equation and the system is thus not well balanced. These errors create artificial values of the velocity near the interface which should theoretically be zero. These artificial velocities are often called parasitic currents. The amplitude of the parasitic currents is an indication of the stability and the accuracy of the numerical method, and especially of the pressure computational step. Indeed, they are the only source of numerical errors.

In what follows several comparisons with other references methods, the Ghost Fluid and the CSF methods in Section 6.1.1 and the Volume of Fluid method in Section 6.1.2, are presented on slightly different test cases.

Comparison with the Ghost Fluid and the CSF Methods

We use the same parameters as in [START_REF] Couderc | Développement d'un code de calcul pour la simulation d'écoulements de fluides non miscibles: Application à la désintégration assistée d'un jet liquide par un courant gazeux[END_REF], where a Ghost Fluid and a CSF method were implemented. The amplitude of the parasitic currents, compared to the results in [START_REF] Couderc | Développement d'un code de calcul pour la simulation d'écoulements de fluides non miscibles: Application à la désintégration assistée d'un jet liquide par un courant gazeux[END_REF], are reported in Table 1. The L ∞ and L 2 norms are computed over the whole domain Ω. The initial configuration is described in Figure 5. As it can be observed, the amplitude of the parasitic currents generated by our method is several orders of magnitude smaller than those of the CSF method, and significantly lower than those of the Ghost Fluid method when the grid is refined. Table 1. Comparison between the new method and the numerical results obtained in [START_REF] Couderc | Développement d'un code de calcul pour la simulation d'écoulements de fluides non miscibles: Application à la désintégration assistée d'un jet liquide par un courant gazeux[END_REF] for the Ghost Fluid method and the CSF method for parasitic oscillations, at time t = 1.

                     L = 2 cm, R = 1 cm, ρ int = 1000 kg • m -3 , µ int = 10 -3 Pa • s, ρ ext = 1 kg • m -3 , µ ext = 10 -5 Pa • s, σ = 0.1 N • m -1 (39) 2 L

Ghost Fluid Method

CSF New Method N L ∞ Error L 2 Error L ∞ Error L 2 Error L ∞ Error L 2 Error
16 8.08 ×10 -3 1.88 ×10 -3 3.55 ×10 -2 1.94 ×10 -2 5.21 ×10 -3 7.31 ×10 -5 32 3.42 ×10 -4 7.50 ×10 -5 3.12 ×10 -2 1.18 ×10 -2 9.26 ×10 -5 1.42 ×10 -6 64 5.13 ×10 -5 7.97×10 -6 2.12 ×10 -2 5.44 ×10 -3 1.36 ×10 -5 1.47 ×10 -7 128 2.79 ×10 -5 4.74 ×10 -6 6.44 ×10 -3 1.38 ×10 -3 2.22 ×10 -6 1.92 ×10 -8

Comparison with a Volume of Fluid Method

We now compare the behavior of our method to the Volume of Fluid method developed in [START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF]. The density and viscosity ratio are both chosen to be one for this test case.

The coefficient σ is chosen so as to obtain an Ohnesorge number Oh = µ σρD satisfying

Oh 2 = 1 12000
. The maximum velocity is computed for varying grids at non-dimensional

time t * = t T = 250, with T = Dµ σ .                      L = 1.25 m, R = 1 m, ρ int = 1 kg • m -3 , µ int = 10 -3 Pa • s, ρ ext = 1 kg • m -3 , µ ext = 10 -3 Pa • s, σ = 0.00012 N • m -1 (40)
A comparison between our method and the Volume of Fluid method is presented in Table 2. The new method provides a better accuracy than the Volume of Fluid method for the coarsest grid. As expected, the Volume of Fluid method outperforms our new approach for finer meshes due to more sophisticated schemes near the interface. Nonetheless, Table 2 show a second-order accuracy for our new method. This test case is studied in [START_REF] Raessi | Consistent mass and momentum transport for simulating incompressible interfacial flows with large density ratios using the level set method[END_REF][START_REF] Chenadec | A monotonicity preserving conservative sharp interface flow solver for high density ratio two-phase flows[END_REF], and based on experiments conducted in [START_REF] Martin | An experimental study of the collapse of liquid columns on a rigid horizontal plane[END_REF]. The initial configuration is a water column at rest in air. The initial height and width of the column are both 5.715 cm. The domain size is 40 cm×10 cm. The physical constants are the same than for the rising bubble (Section 6.1.1). For more details, we refer the reader to [START_REF] Raessi | Consistent mass and momentum transport for simulating incompressible interfacial flows with large density ratios using the level set method[END_REF]. We present in Figure 6 the interface evolution at non-dimensional times T = t g/H = 0, 1, 2, 3, 4, with H the initial height of the water column. The computations are performed with 256 × 64 points.

Figure 7 present the temporal evolution of the water front, compared to the experimental results [START_REF] Martin | An experimental study of the collapse of liquid columns on a rigid horizontal plane[END_REF], to results with the Ghost Fluid method used for pressure resolution, and to the conservative method of Raessi and Pitsch [START_REF] Raessi | Consistent mass and momentum transport for simulating incompressible interfacial flows with large density ratios using the level set method[END_REF]. We observe that the front propagation is in agreement with the experimental results and the results of the conservative method [START_REF] Raessi | Consistent mass and momentum transport for simulating incompressible interfacial flows with large density ratios using the level set method[END_REF]. This means that, though the method is not strictly conservative, the numerical errors due to momentum transfer across the interface are not large enough to slow down the propagation of the front. It is not the case for instance for the Ghost Fluid method, as it can be noticed in Figure 7 and has been reported in [START_REF] Raessi | Consistent mass and momentum transport for simulating incompressible interfacial flows with large density ratios using the level set method[END_REF]. In Figure 13, the interface evolution and the vorticity values for the Ghost Fluid method (used for the pressure computation) and the new method are plotted. One observes small oscillations of the vorticity near the interface between two fluids for the Ghost Fluid method, but not for the new method. The conservation of the partial mass of air and momentum of air in each dimension for the new method, and the Ghost Fluid method are depicted in Figures 141516respectively. As air is the less dense fluid, it is more prone to big oscillations if erroneous mass and momentum transfers between the two fluids occur. One can observe that with the new method and its second-order pressure resolution, the oscillations of these quantities are notably less present than for the Ghost Fluid. The partial mass is better conserved too.

Conclusions

We have developed a new method on Cartesian grids for the simulation of incompressible flows with large density ratios. This method relies on a sharp resolution of the pressure term across the interface defined by a level-set function. The advantage of the proposed approach is its simplicity to implement in an existing Cartesian monofluid Navier-Stokes solver as for instance the Ghost Fluid or CSF methods. It is only necessary to modify the stencil for the pressure equation at some irregular grid points by adding one additional point. This Cartesian scheme uses thus additional unknowns located on the interface to discretize with second-order accuracy the jump conditions across the interface. The viscous term is treated with a regularizing approach that allows us to eliminate terms in the jump conditions. Indeed, it has been shown [START_REF] Desjardins | A spectrally refined interface approach for simulating multiphase flows[END_REF][START_REF] Herrmann | The influence of density ratio on the primary atomization of a turbulent jet in crossflow[END_REF] that the regularization has no significant impact on the accuracy of the results. Numerical results show that this method leads to more accurate and stable results than some reference methods as for instance the well known Ghost Fluid or CSF methods. The conservation of volume of each phase is also better conserved. We thus take advantage of both the regularity of the interface defined by the level-set function (curvature is numerically consistent), and the mass conservation of each phase. The CLSVOF [START_REF] Sussman | A coupled level set and volume of fluid method for computing 3D and axisymmetric incompressible two-phase flows[END_REF] can also provide good mass conservation and precision of the interface but with a more complex numerical implementation.

Future works include an extension of the method to three-dimensional problems with interactions with solids including the numerical simulation of wave energy converters [START_REF] Khedkar | The inertial sea wave energy converter (ISWEC) technology: Device-physics, multiphase modeling and simulations[END_REF][START_REF] Bhalla | Simulating water-entry/exit problems using Eulerian-Lagrangian and fully-Eulerian fictitious domain methods within the open-source IBAMR library[END_REF]. We also aim to study the sensibility of the reinitialization procedure of the level-set function defining the bifluid interface, as for instance the one developed in [START_REF] Luddens | Enablers for high-order level set methods in fluid mechanics[END_REF].

Figure 4 .

 4 Figure 4. Left: points involved in the discretization of the elliptic operator at grid nodes M i+2,j+2 and M i+3,j in black, right:Example of stencils for the discretization of jump conditions. Points involved in the discretization of the x-derivative of the pressure at interface point I i,j-1,E and of the y-derivative of the pressure at interface point I i+1,j+1,N in black. For I i,j-1,E both derivatives are expressed with second-order accuracy while for I i+1,j+1,N the left derivative is expressed with second-order and the right derivative with first-order accuracy.

Figure 5 .

 5 Figure 5. Test case of the static bubble with parasitic oscillations.

Figure 6 .

 6 Figure 6. Evolution of the interface for the dam break problem at non-dimensional times T = t g/H = 0, 1, 2, 3, 4.

Figure 7 .

 7 Figure 7. Evolution of the front of propagation: comparison between experimental data and several numerical methods: the Ghost Fluid method (non-conservative method), the conservative method of Raessi and Pitsch and our new method, The dimensionless location of the front z a is plotted as a function of the dimensionless time T = t g/H.

Figure 9 .

 9 Figure 9. Evolution of the interface and vorticity values for the small bubble test case, at times 0, 0.02, 0.035, 0.05, resolution: 40 × 60 (left); 80 × 120 (middle); 160 × 240 (right).

Figure 10 .

 10 Figure 10. Evolution of the interface and vorticity values for the large bubble test case, at times 0, 0.2, 0.35, 0.5, resolution: 40 × 60 (left); 80 × 120 (middle); 160 × 240 (right).

Figure 12 .

 12 Figure 12. Initial fluid domain for the test case of the two rising bubbles in water.

Figure 13 .

 13 Figure 13. Cont.

Figure 13 .

 13 Figure13. Evolution of the interface and vorticity values for the two bubbles test case at times t = 0., 0.04, 0.08, 0.12, resolution 120 × 180, (left): new method for pressure, (right): Ghost Fluid method for pressure.

Figure 14 .Figure 15 .

 1415 Figure 14. Evolution in time of partial mass of air with N = 60, 90, 120 points in x-direction with (a) Ghost Fluid method, (b) new method.

Figure 16 .

 16 Figure 16. Evolution in time of partial y-momentum of air N = 60, 90, 120 points in x-direction with (a) Ghost Fluid method, (b) new method.

Table 2 .

 2 Numerical results for parasitic oscillations at non-dimensional time t * = 250 for[START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF] and our method.

	∆x	Error L ∞ for [12]	Error L ∞ for Our Method
	2.5/16	7.3 ×10 -4	7.48 ×10 -5
	2.5/32	4.5 ×10 -6	4.7 ×10 -6
	2.5/64	5.5 ×10 -8	1.26 ×10 -6
	6.2. Collapse of a Water Column: the Dam Break Problem

Funding: This study has been carried out with financial support from the French State, managed by the French National Research Agency (ANR) in the frame of the Investments for the future Programme IdEx Bordeaux (ANR-10-IDEX-03-02).

Rising of Air Bubble in Water

We study the evolution of fluid bubbles rising in an heavier fluid, and compare our results to several methods in the literature. The initial configuration is described in Figure 8.

Comparison with the Ghost Fluid Method

We consider air bubbles rising in water, as in the test case proposed for the Ghost Fluid method in [START_REF] Kang | A Boundary Condition Capturing Method for Multiphase Incompressible Flow[END_REF]. The value of the physical parameters are

ρ water = 1000 kg/m 3 , µ water = 1.137 × 10 -3 kg/ms, ρ air = 1.226 kg/m 3 , µ air = 1.78 × 10 -5 kg/ms, σ = 0.0728 kg/s 2 , g = -9.8 m/s 2 .

(41)

We consider two cases: a small bubble with R = 1/300 m and a large one R = 1/3 m. In the first case, the surface tension plays an important role in the evolution of the interface because of the high bubble curvature. In the second case, the surface tension has less influence, and larger deformations occur. The interface of the small bubble and the vorticity values are plotted at times t= 0, 0.02, 0.035, 0.05 in Figure 9. The interface of the large bubble and the vorticity values are plotted at times t= 0, 0.2, 0.35, 0.5 in Figure 10. Our numerical results are in good agreement with [START_REF] Kang | A Boundary Condition Capturing Method for Multiphase Incompressible Flow[END_REF]. 6 R 9 R in [START_REF] Kang | A Boundary Condition Capturing Method for Multiphase Incompressible Flow[END_REF] or 10 R in [START_REF] Grenier | An Hamiltonian interface SPH formulation for multi-fluid and free-surface flows[END_REF] 2 R in [START_REF] Grenier | An Hamiltonian interface SPH formulation for multi-fluid and free-surface flows[END_REF] 3 R in [START_REF] Kang | A Boundary Condition Capturing Method for Multiphase Incompressible Flow[END_REF] gaz liquid R

Comparison with SPH and the Level-Set Method

This test case is taken from [START_REF] Grenier | An Hamiltonian interface SPH formulation for multi-fluid and free-surface flows[END_REF], and inspired from a test case presented in [START_REF] Sussman | A level-set approach for computing solutions to incompressible two-phase flows[END_REF]. It gives us the opportunity to compare our method to another class of methods, based on the SPH formulation. The initial configuration is described on Figure 8. The values of the physical parameters are

025 m, ρ water = 1000 kg/m 3 , µ water = 1.137 × 10 -3 kg/ms, ρ air = 1.226 kg/m 3 , µ air = 1.78 × 10 -5 kg/ms, σ = 0.0728 kg/s 2 , g = -9.8m/s 2 .

(42)

The evolution of the interface and the vorticity values are plotted on Figure 11, for 120 × 200 grid points. We observe that the interface deforms in a way similar to the results in [START_REF] Grenier | An Hamiltonian interface SPH formulation for multi-fluid and free-surface flows[END_REF].

Two Air Bubbles in Water

Two air bubble are initially at rest in water, see Figure 12. At the same time that they are rising, their interaction produces larger deformations than for a single bubble. This test case is meant to assessing the conservativity of the new method, as the increase in pression resolution is meant to increase this conservativity.

The value of the physical parameters are

Conflicts of Interest:

The authors declare no conflict of interest.


Federico Zaccaria
email: federico.zaccaria3@unibo.it

Edoardo Idá
email: edoardo.ida2@unibo.it

Sébastien Briot
email: sebastien.briot@ls2n.fr

Marco Carricato
email: marco.carricato@unibo.it

Workspace Computation of Planar Continuum Parallel Robots

Keywords: Continuum Parallel Robots, Workspace Evaluation, Kantorovich Theorem

Continuum parallel robots (CPRs) comprise several flexible beams connected in parallel to an end-effector. They combine the inherent compliance of continuum robots with the high payload capacity of parallel robots. Workspace characterization is a crucial point in the performance evaluation of CPRs. In this paper, we propose a methodology for the workspace evaluation of planar continuum parallel robots (PCPRs), with focus on the constant-orientation workspace. An explorative algorithm, based on the iterative solution of the inverse geometrico-static problem is proposed for the workspace computation of a generic PCPR. Thanks to an energy-based modelling strategy, and derivative approximation by finite differences, we are able to apply the Kantorovich theorem to certify the existence, uniqueness, and convergence of the solution of the inverse geometrico-static problem at each step of the procedure. Three case studies are shown to demonstrate the effectiveness of the proposed approach.

I. INTRODUCTION

C ONTINUUM robots (CRs) [1] are manipulators that achieve movement through controlled displacement and deformation of slender elastic links. CRs are well suited for maneuvers in complex curvilinear pathways and when intrinsic flexibility is important. In particular, continuum parallel robots (CPRs) employ multiple slender links connected in parallel to an end-effector (EE) and, compared with their serial counterparts, they may exhibit higher precision, stability, and payload capacity [START_REF] Black | Parallel continuum robots: Modeling, analysis, and actuation-based force sensing[END_REF], [START_REF] Chen | Kinetostatics modeling and analysis of parallel continuum manipulators[END_REF]. As recently highlighted [START_REF] Mauzé | Nanometer precision with a planar parallel continuum robot[END_REF], CPRs are good candidates for high precision manipulation tasks, or applications where intrinsic compliance is required and serial CRs suffer from payload limitations. While the basic architecture of a CPR is similar to a Gough-Stewart platform, different designs may lead to an increase in performance [START_REF] Childs | Leveraging geometry to enable highstrength continuum robots[END_REF]. In order to assess which design better suits a specific operational task, performance quantification is required.

Workspace evaluation, i.e. the identification of all poses where the robot may lie in a stable equilibrium, is a crucial performance assessment tool. Though several geometrical, discretization and numerical methods are available in the literature for rigid-link manipulators [START_REF] Merlet | Parallel robots[END_REF], workspace computation algorithms for CRs are at a preliminary stage. Joint-space sampling approaches are used in [START_REF] Burgner-Kahrs | Workspace characterization for concentric tube continuum robots[END_REF], [START_REF] Trivedi | Dexterity and workspace analysis of two soft robotic manipulators[END_REF], for the computation of the CR reachable workspace. However, these approaches are computationally expensive to achieve accuracy. To overcome this disadvantage, an approximated approach is proposed in [START_REF] Cao | Workspace analysis of tendon-driven continuum robots based on mechanical interference identification[END_REF]. However, because this method can detect only points on the external surface of the workspace, it seems difficult to extend it to the identification of inner boundaries. An optimization-based approach is proposed in [START_REF] Amehri | Discrete cosserat method for soft manipulators workspace estimation: an optimizationbased approach[END_REF] for the workspace computation under a general modelling strategy that can encompass various CRs. However, even in this case, inner boundaries are not easy to detect.

As far as CPRs are concerned, the workspace of a 3-leg CPR is calculated in [START_REF] Singh | Optimal work space of parallel continuum manipulator consisting of compact bionic handling arms[END_REF] by intersecting the workspace of each leg, with the latter being independently obtained by sampling the actuation space and computing the resulting EE poses. This approach, also used in [START_REF] Nuelle | Modeling, calibration, and evaluation of a tendon-actuated planar parallel continuum robot[END_REF] for a planar CPRs (PCPRs), allows the results to be obtained in a straightforward way, but it may lack accuracy since no interaction force is modelled between flexible legs and the rigid EE. A different strategy is employed in [START_REF] Orekhov | Analysis and validation of a teleoperated surgical parallel continuum manipulator[END_REF] for a 6-DoF CPR: the actuation space is sampled uniformly and, by iteratively solving the forward geometricostatic problem, the reachable workspace is obtained. Strain limits are considered by excluding configurations that overcome a certain threshold. This approach can be employed for a generic CPR, but singularity and stability analyses are not performed. Equilibrium stability is instead analyzed in [START_REF] Till | Elastic stability of cosserat rods and parallel continuum robots[END_REF], from an optimal-control point of view. The above limitations are also encountered in [START_REF] Mauzé | Nanometer precision with a planar parallel continuum robot[END_REF], where an inverse geometricostatic discretization approach is successfully used for the evaluation of the constant orientation workspace of PCPRs, but the modelling strategy does not include singularity and stability analysis. Similarly, an inverse kinematics discretization strategy is used in [START_REF] Lilge | Tendon Actuated Continuous Structures in Planar Parallel Robots: A Kinematic Analysis[END_REF], and singular configurations are additionally detected by employing a simplified mathematical model; however, the latter cannot account for the action of external loads on the robot. The workspace of a two-leg CPR is computed in [START_REF] Altuzarra | Position analysis in planar parallel continuum mechanisms[END_REF], by including the possibility of analyzing different working modes (i.e multiple solutions of the inverse problem [START_REF] Bonev | Working and assembly modes of the agile eye[END_REF]). However, the modelling strategy, based on the solution of elliptic integrals and limited to planar robots, does not seem to be easily applicable to a generic CPR.

A fundamental part of each workspace evaluation algorithm is the computation of the robot pose. In CPRs, this task is not trivial, since the geometry of the manipulator is not sufficient 0000-0000/00$00.00 c 2021 IEEE to describe the pose of the whole robot, and its configuration is defined by the elastic deformation of the links. This problem is called geometrico-static and, in CPRs, forward and inverse problems do not generally admit an analytical solution. Consequently, the geometrico-static model is usually simplified by introducing numerical approximations [START_REF] Simo | A three-dimensional finite-strain rod model. Part II: Computational aspects[END_REF] with the approximate solution being computed by numerical schemes (e.g. Newton-based methods). However, it may happen that no or multiple solutions exist for the same problem, and in this case solution certification becomes a crucial point. Interval analysis (IA) applies a branch-and-prune approach to solve a set of equations: it is mainly employed (i) when it is needed to certify that all solutions are found in a bounded region [START_REF] Altuzarra | Certified kinematics solution of 2-dof planar parallel continuum mechanisms[END_REF], and (ii) when dealing with bounded model parameter uncertainties [START_REF] Tan | Pose characterization and analysis of soft continuum robots with modeling uncertainties based on interval arithmetic[END_REF], [START_REF] Iqbal | A guaranteed approach for kinematic analysis of continuum robot based catheter[END_REF]. As a drawback, computational time of IA methods exponentially depends on the number of variables [START_REF] Merlet | Interval analysis for certified numerical solution of problems in robotics[END_REF]. On the other side, the Kantorovich theorem [START_REF] Kantorovich | On Newton's method for functional equations[END_REF] is a useful tool when the solution certification is important. This theorem establishes sufficient conditions on the initial guess of the Newton iteration to certify the existence, uniqueness, and convergence of the solution on a defined region [START_REF] Goldsztejn | A three-step methodology for dimensional tolerance synthesis of parallel manipulators[END_REF]. While the application of the Kantorovich theorem is computationally efficient, certification of multiple solutions, and dealing with uncertainties are not possible, in contrast to IA approaches.

Our goal is to develop a methodology for the workspace evaluation of PCPRs, a class of CPRs which recently obtained higher scientific interest [START_REF] Mauzé | Nanometer precision with a planar parallel continuum robot[END_REF], [START_REF] Lilge | Tendon Actuated Continuous Structures in Planar Parallel Robots: A Kinematic Analysis[END_REF]. We want to include several aspects usually of interest for workspace evaluation that stateof-the-art approaches cannot include simultaneously, namely equilibrium stability, singularity identification, strain limits, and external load influence. For this purpose, we decided to adopt the modelling strategy proposed in [START_REF] Zaccaria | An Analytical Formulation for the Geometrico-static Problem of Continuum Planar Parallel Robots[END_REF], since it allows us to apply external loads to any robot component, to enforce generic geometric constraints, to assess equilibrium stability, and to identify singular configurations (by the procedure introduced in [START_REF] Briot | Singularity conditions for continuum parallel robots[END_REF]).

Our workspace-evaluation algorithm is based on the iterative solution of the inverse geometrico-static problem (IGSP) [START_REF] Briot | Singularity conditions for continuum parallel robots[END_REF]. In order to certify the existence of a unique solution of the IGSP at each step of the workspace computation algorithm, and thus certify that all the computed configurations are characterized by the same working mode, we propose to apply the Kantorovich theorem. With this theorem, it is possible to ensure that a computed solution lies in the neighborhood of a previously computed configuration, and to certify the convergence of the numerical method used to find the new IGSP solution. By employing an energy-based modelling strategy combined with finite differences [START_REF] Zaccaria | An Analytical Formulation for the Geometrico-static Problem of Continuum Planar Parallel Robots[END_REF], we can obtain an analytical expression of the constants required for the solution certification through the Kantorovich theorem that may be not trivial to obtain with other approaches (i.e. the continuous formulation in [START_REF] Black | Parallel continuum robots: Modeling, analysis, and actuation-based force sensing[END_REF], which may require the solution of an optimization problem). As it will also be shown, this methodology for the workspace computation is time-efficient 1 . It allows obtaining certified results in a reduced computational time, which is beneficial since fast workspace computation algorithms may be used for workspace-guided design tools.

We additionally propose the use of a preconditioning matrix for the Newton iteration, to alleviate the numerical illconditioning of the IGSP problem and, consequently, to facilitate the IGSP solution certification.

Finally, we propose an adaptive-grid workspace-exploration algorithm. At each step of the workspace algorithm, the grid is refined in case the solution of the IGSP is not certifiable, so that a refined grid is used only where necessary (in contrast to the fixed grid flooding algorithm of [START_REF] Briot | Singularity conditions for continuum parallel robots[END_REF]): the overall computation time is therefore reduced. The algorithm is then employed to identify stable and unstable workspace regions, and to detect singularities that determine workspace boundaries as in [START_REF] Briot | Singularity conditions for continuum parallel robots[END_REF].

The paper is structured as follows. Section II introduces the PCPR modelling strategy, analyzes configuration stability and singularities, and discusses the Kantorovich theorem application. Section III describes the novel workspace exploration algorithm in detail. Section IV proposes three case studies to show the effectiveness of our approach. Section V draws conclusions and outlines future work directions.

II. MODELLING

A. Modelling Approach

Consider a slender beam of length L as represented in Fig. 1a. Shear and extensibility are neglected. Each rod is discretized into N elements of equal length L e = L/N and the orientation of the i-th element (i = 1, ..., N ), with respect to (w.r.t.) a fixed frame, is θ i . The deformation energy of the beam, assumed as initially straight, can be written as:

V b (θ 1 , ..., θ N ) = V e + V l (1)
where V e and V l represent beam's deformation energy and the potential energy of concentrated and distributed loads acting on the beam, computed by approximating derivatives with finite differences (see [START_REF] Zaccaria | An Analytical Formulation for the Geometrico-static Problem of Continuum Planar Parallel Robots[END_REF], [START_REF] Briot | Singularity conditions for continuum parallel robots[END_REF] for their expressions). We assume that the PCPR comprises n beams, as in Fig. 1b, with the generic beam being identified by the index k = 1, ..., n. We can identify n controlled variables 2 , collected in the vector q A = [q a1 , ..., q an ]. Moreover, we define

θ k = [θ 2k , ..., θ N k ] and θ = [θ 1 , ..., θ n ].
The position and orientation of the EE are denoted as p P and φ, respectively, and grouped in the array q P = [p P , φ]. The array x = [θ, q P ] contains all non-actuated variables. Finally, the total potential energy of the PCPR can be written as:

V t (q A , x) = n k=1 V b,k (q ak , θ k ) -f T (p P -p * P ) -m(φ -φ * )
(2) where f , m are, respectively, a EE constant force and a constant moment (which ensures that the overall system is conservative), and (.) * denotes the undeformed state.

Due to the parallel architecture, the configuration variables are related by geometric constraints, such as closure equations. If beams are connected to the platform by revolute joints, closure equations can be written as:

p Bk (q Ak , θ k ) -p Bk (p P , φ) = 0 (3)
where p Bk (q Ak , θ k ) is function of the leg variables [START_REF] Zaccaria | An Analytical Formulation for the Geometrico-static Problem of Continuum Planar Parallel Robots[END_REF] and p Bk (p P , φ) is function of platform variables. Other beam connections to the EE can be modelled (such as fixed constraints [START_REF] Mauzé | Nanometer precision with a planar parallel continuum robot[END_REF]): Φ is the array stacking m φ general geometric constraints.

A static robot configuration is feasible, for an assigned value of q A , if and only if the manipulator is in a stable equilibrium, that is, the total potential energy is at a minimum. However, variables are related by constraints Φ, and local extrema of V t are characterized by Lagrange conditions [START_REF] Nocedal | Numerical Optimization[END_REF]. Assuming that ∇ x Φ has full rank, x is a local extrema of V t if Lagrange multipliers λ exist such as:

∇ x L = 0 , ∇ λ L = 0 (4)
where the Lagrangian function L is:

L(q A , x, λ) = V t (q a , x) + Φ(q a , x) T λ (5)
In order to assess if a configuration belongs to the static workspace, we are interested in determining if a value of q A exists so that an assigned EE pose q d P is an equilibrium, that is, we seek the solution of the IGSP, defined as: 2 In the case of a revolute joint, q ak = θ 1k whereas, in the case of a prismatic joint, the orientation θ 1k of the first element is a constant value and q ak defines the position of the first node of the beam.

F(y) =      0 = ∇ x L = ∇ x V t + ∇ x Φ T λ 0 = ∇ λ L = Φ 0 = q P -q d P (6)
Equations ( 6) form a square system of N eq = n(N + 2) + 3 nonlinear equations in y = (q A , θ, q P , λ), which can be solved by using root-finding techniques, such as the Newton method.

B. Equilibrium Stability and Singularity Conditions

Once a solution of Eq. ( 6) is identified, its stability, and its singularities can be evaluated by the method presented in [START_REF] Briot | Singularity conditions for continuum parallel robots[END_REF]. The Jacobian matrix J of Eq. ( 6), can be computed as:

J(y) = ∂F(y) ∂y =   A L U L P L Λ L A Φ U Φ P Φ 0 0 0 I 0   (7)
where:

• A L = ∇ q A (∇ x L) , U L = ∇ θ (∇ x L) • P L = ∇ q P (∇ x L) , Λ L = ∇ λ (∇ x L) • A Φ = ∇ q A Φ , U Φ = ∇ θ Φ , P Φ = ∇ q P Φ
and I is the identity matrix. The analytical expression of all terms in J is obtainable, but not reported here for brevity. Equilibrium stability can be assessed by verifying that the reduced hessian matrix H r of the total potential energy is positive definite [START_REF] Nocedal | Numerical Optimization[END_REF]. H r can be obtained as:

H r = Z T ∂ 2 V t ∂x 2 Z = Z T U L P L Z (8)
where Z denotes the left nullspace of Λ L . Moreover, configuration singularities can be assessed if the following matrices degenerate [START_REF] Briot | Singularity conditions for continuum parallel robots[END_REF]:

T 1 = Z T A L Z T U L A Φ U Φ , T 2 = Z T P L Z T U L P Φ U Φ (9)
Type 1 (or serial) singularities are related to the degeneracy of T 1 and define the limits of the robot workspace. Type 2 (or parallel) singularities, related to T 2 , delimit stable-to-unstable regions.

C. Kantorovich Theorem for Newton Iteration and Preconditioning

If a Newton method is employed for the solution of the IGSP in Eq. ( 6), the Kantorovich theorem establishes sufficient conditions for the existence and uniqueness of the solution [START_REF] Kantorovich | On Newton's method for functional equations[END_REF]. Let χ, δ, λ ∈ R + , y 0 be the initial guess given to the numerical solver, and B(y 0 , 2δ) = {y : yy 0 ≤ 2δ} a ball of radius 2δ centered in y 0 . For convenience, let the infinite norm being used. Let J(y) be the Jacobian matrix of F(y) w.r.t. y, and χ, δ, λ are chosen such as:

χ ≥ J -1 (y 0 ) (10) δ ≥ J -1 (y 0 )F(y 0 ) (11) λ ≥ max y∈B   max h∈[1,Neq]   i,j ∂ 2 F h (y) ∂y i ∂y j     (12)
Then, if χ, δ, λ exist such that 2χδλ ≤ 1, the Kantorovich theorem states that a solution y * ∈ B(y 0 , t * ) exists, where

t * = 1- √ 1-2χδλ λχ ∈ [δ, 2δ].
Moreover, the solution is unique inside B(y 0 , 2δ), and the Newton iteration, starting from y 0 , converges to y * . Kantorovich constants have a numerical meaning: χ is related to the absolute conditioning of the problem, δ represents the closeness to the linearized solution, and λ is influenced by the non-linearity of the problem.

However, if a discretization approach is used for modelling PCPRs, the Jacobian matrix J is usually large and ill-conditioned, resulting in high values of χ. Therefore, the certification of the existence, uniqueness, and convergence of the IGSP solution holds only for y 0 sufficiently close to y * . To overcome this difficulty, we can preconditionate the IGSP to reduce the value of χ. We propose to solve an equivalent and modified inverse geometrico-static problem (MIGSP), where each IGSP equation is multiplied by a constant that represents the physical dimensions of the equation terms:

F M (y) =                    0 = 1 EI (∇ θ V t + ∇ θ Φ T λ) 0 = L EI (∇ p P V t + ∇ p P Φ T λ) 0 = 1 EI (∇ φ V t + ∇ φ Φ T λ) 0 = 1 L Φ 0 = 1 L (p P -p d P ) 0 = φ -φ d (13)
or equivalently:

F M (y) = MF(y) (14)
where M is a diagonal matrix collecting the coefficients that multiply each equation of Eq. ( 13). The Jacobian matrix of the MIGSP w.r.t. y is:

J M (y) = ∂F M (y) ∂y = M ∂F(y) ∂y = MJ(y) (15)
The Kantorovich constants of the MIGSP can be obtained from Eqs. [START_REF] Amehri | Discrete cosserat method for soft manipulators workspace estimation: an optimizationbased approach[END_REF], [START_REF] Singh | Optimal work space of parallel continuum manipulator consisting of compact bionic handling arms[END_REF], [START_REF] Nuelle | Modeling, calibration, and evaluation of a tendon-actuated planar parallel continuum robot[END_REF] by employing F M , J M instead of F, J. While the computation of χ, δ is straightforward, the computation of λ requires to analytically compute all the second derivative of MIGSP equations, which is long and involved. According to Eq. ( 12), we have to identify where the sum in absolute value of the second derivatives assumes the maximum value inside the ball for each equation. However, this sum is not constant w.r.t. y and we may identify the maximum by solving a constrained optimization problem, at the cost of high computational time. We preferred to approximate λ with absolute value inequalities [START_REF] Nocedal | Numerical Optimization[END_REF]: in this way we consider the worst-case scenario and we obtain an expression of λ that depends on y 0 , δ. The details of this computation are not provided here due to space limitations.

III. WORKSPACE EVALUATION

This Section describes in detail the innovative workspacecomputation strategy, which we call Adaptive Flooding Algorithm (Alg. 1 lines 1-19). In particular, we modified the Flooding Algorithm of [START_REF] Briot | Singularity conditions for continuum parallel robots[END_REF] by introducing a grid-adaptation procedure: during each iteration, the grid size is adjusted when the certification is not feasible with the initial grid size. This is done in order to certify as much workspace as possible, since the Kantorovich constant δ depends on the distance of the initial guess from the solution. We focus on the constant-orientation workspace of PCPRs, i.e. the set of all Algorithm 1: Adaptive flooding algorithm. possible locations of the robot EE that can be reached with a given orientation, though our approach can be easily extended to other types of workspace. A discretization approach is employed and the MIGSP is solved iteratively over a grid of EE locations. The first step of the algorithm requires the initialization of some entities: a 2-dimensional uniform grid is generated accordingly to an initial stepsize (s i ). The grid discretizes a user-defined box (assumed to fully include the workspace) and in the middle of each square of the grid, a point representing the EE position is placed (Fig. 2a). An initial point (named p P ) where the explorative algorithm starts is selected within the grid. For the first solution of the MIGSP, an initial guess obtained through constant-curvature modelling approach is employed. As shown in [START_REF] Lilge | Tendon Actuated Continuous Structures in Planar Parallel Robots: A Kinematic Analysis[END_REF], the inverse problem under constant-curvature approach admits a finite number of solutions (but usually low solution accuracy), and we can decide which solution to use as initial guess, in order to let the algorithm converge to the desired working mode [START_REF] Bonev | Working and assembly modes of the agile eye[END_REF]. The resulting output configuration is stacked in the Config list to be employed as a future guess. Then, neighboring EE positions (Fig. 2a) are computed and saved in the toDo list, which contains points to be processed. The toDoEnd list, that contains points that are neighbors to Type 2 singular configurations, is initialized as empty. Then, the iterative algorithm starts, a new p P is obtained from the first element of the toDo list and employed as the EE location to be reached. The initial guess y 0 for the MIGSP, the Kantorovich flag associated to y 0 and neighbors of p P are obtained through the FindGuess procedure. The latter, which plays an important role in the certification of the MIGSP solution, is detailed in subsection III-A.

If the current step s c (i.e the edge size of the square centered in p P ) is greater to the minimum stepsize allowed for the computation (s m ), the value of the Kantorovich constants is checked to certify the solution: if 2χδλ ≤ 1, the MIGSP is solved and new points to be processed are created in the Compute procedure (detailed in subsection III-B). If the manipulator is not cuspidal [START_REF] Wenger | Cuspidal and noncuspidal robot manipulators[END_REF], and the solution is certified, the working mode of the initial guess is preserved also on the configuration resulting from the MIGSP solution. In the case of the solution is not certifiable, the grid is refined: the square centered in p P is divided into new four equal squares (Fig. 2b), and new four points p new are placed in the middle of each new square. Then, the original EE location p P is replaced with p new and added to the toDo list. In the case the minimum stepsize is reached (s c ≤ s m ), the Compute procedure is executed even if the MIGSP solution is not certified, in order to fully compute the workspace. Finally, the algorithm restarts until some elements are present in the lists toDo and toDoEnd. If toDo is empty, then it is refilled with toDoEnd (how toDo and toDoEnd are managed is explained in subsection III-B).

A. Choice of the Initial Guess

The choice of the initial guess at each iteration plays an important role in the MIGSP solution certification. The routine of the initial guess selection is described in lines 20-32 of Alg. 1. Given a desired EE location p P , we want to identify an initial-guess configuration y 0 to be used for the MIGSP solution. Initially, the distance r N for which two EE locations are considered to be neighbors is set as √ 2s c (Fig. 2a). This way, neighbors of p P are identified in the grid and stacked in the array n. Neighbors in the workspace are extracted from n and collected in the array n W K . However, caused by the grid refinement process, it can happen that no neighbors are in the workspace (Fig. 2c): in this case, the radius r N for which points are considered to be neighbors is multiplied by two, and the selection of n is repeated until workspace points are found.

In order to increase the possibility to certify the MIGSP solution, the Kantorovich constants for the configurations associated to n W K EE locations are computed. Then, the EE location that ensures the lowest value of 2χδλ (named n best ) is identified, and the robot configuration y 0 associated to n best is extracted from Config. If the solution certification is not required, n best can be chosen as the one that ensures the best inverse conditioning of J to speed up the computation.

B. Computation Process

This subsection describes the routine for the computation of the MIGSP solution and for the creation of new points to be processed (lines 33-46 of Alg. 1). Starting from a given initial guess y 0 and a set of neighboring EE locations n, the MIGSP is solved by a Newton scheme. Then, T 1 , T 2 are computed according to Eqs. [START_REF] Cao | Workspace analysis of tendon-driven continuum robots based on mechanical interference identification[END_REF]. In particular, we compute the inverse condition number of T 1 , T 2 in order to identify their degeneracy. Then, we check if the Newton solver converges and the resulting configuration is not Type-1 singular. Moreover, mechanical constraints are verified in mechconstr: these include strain limits on the legs, as well as joint limits. If the check succeeded, outputs associated with the resulting configuration y (e.g. internal energy of the robot, equilibrium stability, number of inflection points) are computed in Values. These results, as well as singularity flags and Kantorovich constants, are saved in Results and y is stored in Config as a future initial guess. Subsequently, neighbors not in the workspace and not in toDo,toDoEnd, are stored in n 1 . If the actual configuration y is not T 2 -singular, n 1 is added to the toDo list, else in toDoEnd.

In this way, Type 1 singularities, associated with boundaries of the workspace, are not crossed but only approached. Type 2 singularities, which delimit stable from unstable regions, are crossed in a second stage of the algorithm (only when toDo is empty) in order to discover possible stable regions separated by unstable transitions.

IV. CASE STUDIES

In this Section, three case studies, focusing on different PCPRs (Fig. 3), are reported. Workspace evaluation is discussed with a focus on the identification of different regions (e.g. stable, unstable, regions where stress limits are exceeded) and certification of the IGSP solution. Our algorithm is compared to the flooding algorithm of [START_REF] Briot | Singularity conditions for continuum parallel robots[END_REF], to demonstrate the benefit of the grid-adaptation routine in terms of computational time. For all case studies, beams are made of harmonic steel with Young modulus E = 210 GPa, maximum stress σ max = 1800 MPa, density ρ = 7800 kg/m3 , length L = 1m , circular cross-section of radius r = 1mm. Simulations are performed in the Matlab enviroment.

A. RF RF R robot

This subsection investigates the workspace of a RF RF R robot. The aim of this case study is twofold: on the one hand, it shows the capability of our algorithm to detect singularities and unstable regions, as well as to include external loads and strain limits in the model; on the other hand, it investigates the influence of the stepsize and the preconditioner on the IGSP solution certification.

This manipulator, borrowed from [START_REF] Altuzarra | Position analysis in planar parallel continuum mechanisms[END_REF], [START_REF] Zaccaria | An Analytical Formulation for the Geometrico-static Problem of Continuum Planar Parallel Robots[END_REF], has two actuated revolute joints in A 1 , A 2 (R) and two flexible links (F ) connected by a passive revolute joint centered in B (Fig. 3a). The distance between the actuators is L A1A2 = 0.4m. An external force of 1.5N is applied on the EE, and legs are subjected to gravity. Simulations are performed with N = 50, ensuring sufficient MIGSP solution accuracy.

The workspace of the RF RF R, computed by the algorithm presented in Section III, is shown in Figs. 4a,4b,4c. Configurations are marked as singular when the inverse condition number of matrices reported in Eq. ( 9) is lower than a certain threshold T OL. We experienced that T OL = 10 -6 correctly identify singularities. Stress limits are considered by evaluating whether a first-order approximation of stress [START_REF] Antman | Nonlinear Problems of Elasticity[END_REF] on each element does not overcome σ max . By considering stress limits, the workspace is considerably reduced, in particular in the working mode showed in Fig. 4c (0.96 m 2 without considering stress limits to 0.36 m 2 with stress limit included).

We also tested the influence of the minimum stepsize s m and the influence of the preconditioner on the MIGSP certification of the workspace, with a focus on the stable and feasible workspace (i.e where stress limits or joint bounds are not exceeded). To quantify how many configurations are computed in a certified way, we introduce the certified percentage of the workspace C % as:

C % = 100 • A C /A W ( 16
)
where A C is the certified area and A W is the workspace area. These areas are obtained by summing the area of each square of the grid that lies in the workspace (for A W ) and by summing the area of each square that is computed in a certified way and belongs to the workspace for A C . Results are reported in Table I. With s i = 4mm, by passing from s m = 2mm (Fig. 4a) to s m = 1mm (Fig. 4b), C % grows from 80.4% to 92.1%. As expected, the computational time 3 increases from 32 min to 84 min. To achieve C % = 92.1% the flooding algorithm of [START_REF] Briot | Singularity conditions for continuum parallel robots[END_REF] required 192 min with 1 mm stepsize, which is considerably higher. If the preconditioner is not used, C % drops to 35.1% (with s m = 1mm) and at the same time, the computational time reaches 383 min. The increase of the computational time is due to the higher number of processed points (18.9 • 10 5 compared to 7.1 • 10 5 in the preconditioned case) for the adaptation of the grid, and not for a considerable increase in the MIGSP solution time. This is confirmed by the data relative to the case of s m = s i = 4mm (where no grid refinement is possible), with the same grid being employed with and without the preconditioner, and the resulting computational times are comparable. Also in the case of a different working mode (Fig. 4c), we certified a significant amount of the workspace (70.2%) in a reduced time (44 min) with s m = 1mm.

B. 3 -RF R robot

This subsection investigates the workspace of the 3 -RF R robot. This case study shows the possibility of identifying the workspace and certifying the IGSP solution with different EE orientations. This manipulator borrowed from [START_REF] Lilge | Tendon Actuated Continuous Structures in Planar Parallel Robots: A Kinematic Analysis[END_REF] has three actuated revolute joints A 1 , A 2 , A 3 (R) and three flexible links (F ) connected by passive revolute joints B 1 , B 2 , B 3 at a rigid EE. Actuators are equally placed along a circumference of radius r b = 0.6m, whereas passive revolute joints are placed on a circumference of radius r p = 0.15m (Fig. 3b). Simulations are performed with N = 30 and no external loads are included.

The workspace of the 3-RF R robot is illustrated in Fig. 4d and 4e by fixing φ = -π 3 , φ = -π 6 , respectively. In both cases, no unstable regions and Type-2 singularities are detected. Stress limits are included, but no point exceeds σ max . As before, s m = 1mm is chosen to guarantee a sufficient value of C % (81% and 84.6%) in a reasonable computational time (11 and 10 min). Again, the flooding algorithm of [START_REF] Briot | Singularity conditions for continuum parallel robots[END_REF] required higher computational time to obtain the same C % (20 and 18.5 min with 1 mm stepsize).

C. 3 -P F R robot

This subsection studies the workspace of the 3-P F R robot. This case study shows the possibility of analyzing PCPRs with different actuators and including joint limits. This manipulator is similar to the one proposed in [START_REF] Mauzé | Nanometer precision with a planar parallel continuum robot[END_REF], except for the connection of the flexible links with the platform (passive revolute joints in our case, in contrast with fixed connections in [START_REF] Mauzé | Nanometer precision with a planar parallel continuum robot[END_REF]). The 3 -P F R robot has three actuated prismatic joints A 1 , A 2 , A 3 (P ) and three flexible links (F ) connected by passive revolute joints B 1 , B 2 , B 3 at a rigid EE. Actuators are equally spaced along a circumference of radius r b = 0.65m, and passive revolute joints are placed on the platform of radius r p = 0.15m (Fig. 3c). Simulations are performed with N = 30 and no external loads are included.. The workspace of the 3-P F R robot is illustrated in Fig. 4f, where the EE orientation is φ = π 6 and 1.4m-long rails are symmetrically placed around a circle of radius r b . As for the 3 -RF R robot, maximum strain limits are included, but not exceeded. With s m = 1 mm, we reached C % = 94.5% in 5.2 min of computational time whereas the algorithm of [START_REF] Briot | Singularity conditions for continuum parallel robots[END_REF] required 17 min with 1 mm stepsize. TABLE I: Influence of the preconditioner and the stepsize on the certification of the workspace. Simulations, relative to the workspace displayed in Fig. 4a, 4b, are performed with s i = 4mm, and N = 50.

V. CONCLUSIONS

This paper presented an adaptive flooding algorithm for the workspace computation of PCPRs. The algorithm may identify unstable regions, singularity loci, incorporate external loads, set maximum stress limits and joint bounds. Thanks to an energy-based modelling strategy approximated through finite differences for derivatives, the IGSP solution was certified in terms of existence, uniqueness, and convergence of the solution, by verifying Kantorovich conditions during the Newton-based problem-solving procedure. With this approach, we certified the IGSP solution over a large percentage of the workspace in a reduced computational time in comparison with previous algorithms. However, with large workspaces and/or small stepsizes, the flooding approach may require the computation of a large number of points, which may be not computationally efficient. Future work will be directed on the integration of the adaptive flooding algorithm on a workspace-guided design tool for PCPRs. Moreover, the authors will investigate the workspace evaluation of CPRs and the certification of the solution of the IGSP problem in the spatial case.

Fig. 1 :

 1 Fig. 1: (a) Discretization of a flexible link and (b) representation of a generic PCPR.

Fig. 2 :

 2 Fig. 2: Workspace algorithm: (a) representation of the grid, points along the grid and condition for neighborhood, (b) grid refinement strategy; (c) situation when neighbor radius should be enlarged.

Fig. 3 :

 3 Fig. 3: Three PCPRs object of our case studies: the RF RF R (a), the 3 -RF R (b), and the 3 -P F R (c). Relevant design dimensions are displayed.

Fig. 4 :

 4 Fig. 4: Workspaces of PCPRs. Type 1 and Type 2 singularities are drawn in red and black, respectively. Certified workspace is depicted in blue, non certified workspace in light blue and non certified unstable regions in yellow. Non certified stable regions where stress limits or joint limits are exceeded are represented in green and grey, respectively.

1

  Initialize grid, toDo, toDoEnd, Config, Results; 2 while toDo = ∅ or toDoEnd = ∅ do

	3	if toDo= ∅ then
	4	toDo ← toDoEnd; toDoEnd ← ∅ ;
	5	end
	6	p P = toDo(1), toDo ← toDo\p P ;
	7	[flag k ,y0,n] = FindGuess(p P , sc);
	8	if sc > sm then
	9	if flag k ≤ 1 then
	10	Compute(y0,n);
	11	else
	12	pnew = Generate points(p P ,sc);
	13	Replace p P with pnew;
	14	toDo ← toDo ∪ pnew
	15	end
	16	else
	17	Compute(y0,n);
	18	end
	19 end
	20 Function FindGuess(p P , sc): 21 Set neighbors radius r N = √	2sc;
	22	do
	23	n = Find neighbors EE positions to p P ∈ grid;
	24	n W K = n ∈WK;
	32	return [flag k , y0, n];
	33 Function Compute(y0,n):
	34	y = Solve IGSP starting from y0;
	35	[T1, T2] = Singularity(y);
	36	if (Solver Converged & T1 < T OL & mechconstr(y)) then
	37	Values = CalculateOutputs(y);
	42	else
	43	toDoEnd ← toDoEnd ∪ n1;
	44	end
	45	end
	46	return;

25 if (n W K = ∅) then 26 r N ← 2r N ; 27 end 28 while (n W K = ∅); 29 Compute Kantorovich flag for each n W K ; 30 n best = n W K with best Kantorovich flag; 31 y0 = Config(n best ); 38 Save Values in Results, Save y in Config; 39 n1 = select n / ∈(WK,toDo,toDoEnd); 40 if T2 < T OL then 41 toDo ← toDo ∪ n1;

Since the model proposed in[START_REF] Zaccaria | An Analytical Formulation for the Geometrico-static Problem of Continuum Planar Parallel Robots[END_REF] uses discretized continuum-robot equations, a large number of independent variables is needed for the accurate solution of the IGSP, and high computational time is thus expected for IAbased approaches[START_REF] Merlet | Interval analysis for certified numerical solution of problems in robotics[END_REF] (a) (b)

Results are obtained by a CPU Intel Core i7-8700K,3.7GHz,32Gb RAM