Alfredo Pinelli
email: a.pinelli@ciemat.es

I Z Naqavi

Ugo Piomelli

Julien Favier
email: julien.favier@ciemat.es

Immersed-boundary methods for general finite-difference and finite-volume Navier-Stokes solvers

come

Introduction

Over the last few years several authors have turned their attention to immersedboundary methods (IBMs) for their ability to handle moving or deforming bodies with complex surface geometry embedded in a fluid flow. The key feature of IBMs is the fact that the Eulerian grid is not required to conform to the immersed-body geometry, since the no-slip boundary conditions are enforced on the body surface by appropriate boundary forces.

Peskin [START_REF] Peskin | Flow patterns around heart valves: a numerical method[END_REF] presented the first application of this method. The flow motion was solved using a uniform mesh (referred to as the Eulerian grid in the following). The immersed surface was represented by a set of Lagrangian markers, and the boundary forces took the form of singular functions along the immersed surface in the continuous equations.

They were introduced in the discretized equations in conjunction with regularized delta functions that spread (regularize) the force field over the neighboring Eulerian cells. Peskin originally used the IBM to simulate blood flow inside a heart with flexible valves, and the forcing function was computed using Hooke's law [START_REF] Peskin | Flow patterns around heart valves: a numerical method[END_REF][START_REF] Peskin | Numerical analysis of blood flow in the heart[END_REF]. By considering a large value of the spring stiffness, this method could also be applied to rigid boundaries [START_REF] Lai | An immersed boundary method with formal second-order accuracy and reduced numerical viscosity[END_REF][START_REF] Beyer | Analysis of a one-dimensional model for the immersed boundary method[END_REF].

Other authors obtain the singular forces by a feed-back mechanism [START_REF] Goldstein | Modeling a no-slip flow boundary with an external force field[END_REF][START_REF] Höfler | Navier-stokes simulation with constraint forces: finite-difference method for particle-laden flows and complex geometries[END_REF][START_REF] Feng | The immersed boundary-lattice Boltzmann method for solving fluidparticles interaction problems[END_REF][START_REF] Saiki | Numerical simulation of a cylinder in uniform flow: application of a virtual boundary method[END_REF][START_REF] Lee | Stability characteristics of the virtual boundary method in three-dimensional applications[END_REF]: a deviation from the local desired value of velocity (or position) generates a force in the opposite direction, determined by a system of virtual springs and dampers attached to the Lagrangian markers. Undesirable features of these formulations are the introduction of additional free parameters and the compromise that must be taken, when dealing with rigid objects, between severe restrictions on the time step for negligible deformations or larger time steps accepting mild deformation of the embedded body [START_REF] Lai | An immersed boundary method with formal second-order accuracy and reduced numerical viscosity[END_REF][START_REF] Lee | Stability characteristics of the virtual boundary method in three-dimensional applications[END_REF].

To avoid the drawbacks of the original IBM formulation for rigid objects, Fadlun et al. [START_REF] Fadlun | Combined immersed-boundary finitedifference methods for three-dimensional complex flow simulations[END_REF] introduced a direct formulation of the force term. The direct forcing modifies the discretized momentum equation so that the interpolated velocity at the Lagrangian points takes the desired values. Further developments of this technique were proposed, among others, by Kim et al. [START_REF] Kim | An immersed-boundary finite-volume method for simulations of flow in complex geometries[END_REF], Yang and Balaras [START_REF] Yang | An embedded-boundary formulation for large-eddy simulation of turbulent flows interacting with moving boundaries[END_REF], and Taira and Colonius [START_REF] Taira | The immersed boundary method: A projection approach[END_REF]. Kim et al. [START_REF] Kim | An immersed-boundary finite-volume method for simulations of flow in complex geometries[END_REF] used an explicit variant of the direct forcing method that maintains the matrix structure of a standard finite-difference method. Yang and Balaras [START_REF] Yang | An embedded-boundary formulation for large-eddy simulation of turbulent flows interacting with moving boundaries[END_REF] developed a multi-dimensional interpolation for complex boundary shapes. Taira and Colonius [START_REF] Taira | The immersed boundary method: A projection approach[END_REF] proposed the immersed-boundary projection method, in which the pressure and singular forces are treated together as Lagrangian multipliers in the framework of a projection method.

Fadlun et al. [START_REF] Fadlun | Combined immersed-boundary finitedifference methods for three-dimensional complex flow simulations[END_REF] present an example of a flow involving moving boundaries; during the relative motion, however, the boundary force is not smooth. This was recognized by Uhlmann [START_REF] Uhlmann | An immersed boundary method with direct forcing for the simulation of particulate flows[END_REF], who showed that the interpolation procedure used to relate the force at fixed grid nodes and the arbitrarily located boundaries, in simulations with moving bodies, can lead to undesirable force oscillations. This observation led to the formulation of an alternative direct-forcing scheme in which the force is first computed on the Lagrangian markers, then spread onto the neighboring Eulerian nodes. Both interpolation and spreading use the same discrete Dirac kernel. The algorithm conveys very smooth hydrodynamic forces while preserving the global order of the spatial scheme. Recently, Vanella and Balaras [START_REF] Vanella | A moving-least-squares reconstruction for embedded-boundary formulations[END_REF] have presented an extension of the aforementioned method that yields sharp boundary resolution similar to Eulerian direct-forcing schemes and boundary conforming methods, while keeping the simplicity of the original technique.

In most applications of IBMs the underlying Eulerian grid is Cartesian. Curvilinear grids have been used in fewer cases [START_REF] Moin | Advances in large eddy simulation methodology for complex flows[END_REF][START_REF] Roman | An improved immersed boundary method for curvilinear grids[END_REF]. The common argument against using immersed-boundary methods with curvilinear grids is the increased cost of curvilinear codes, compared to Cartesian ones. However, when the grid points must be clustered close to solid boundaries that are not aligned with grid lines, Cartesian grids may be suboptimal, as the refinement is extended to regions of the flow in which it is not needed, or even to the solid region. Furthermore, in wall-bounded turbulent flows it is generally desirable to use grid cells that are longer in the flow direction than in the other two. Here, the use of meshes in which the grid lines are nearly aligned with the streamlines can lead to very significant savings [START_REF] Roman | An improved immersed boundary method for curvilinear grids[END_REF]. Therefore, one of the principal goals of the present contribution is to develop a method that can be applied to both Cartesian (uniform or non-uniform) and curvilinear Eulerian meshes.

One of the reasons that has limited the use of the IBM techniques to Cartesian grids is due to the discrete delta function, which cannot be trivially extended to more complex grid systems without losing some fundamental properties of the regularized force field. In particular, consider the integral conservation of the force (spread onto the Eulerian grid) and of its first moment:

i,j,k f (x i,j,k )∆v i,j,k = m F(X m ) ∆V m , (1)
i,j,k

x i,j,k × f (x i,j,k )∆v i,j,k = m X m × F(X m )∆V m . (2)
In (1-2), x i,j,k are the Eulerian grid nodes (over which the force has been spread), and X m are the Lagrangian markers; ∆v i,j,k is the volume of the (i, j, k) Eulerian cell and ∆V m the volume defined about the m-th Lagrangian marker. The sums are the discrete counterparts of three-dimensional integrals. These properties ensure that the forces and moments exerted by the fluid on the solid can be calculated correctly by integrating the force field.

Equations (1-2) are verified if the regularized delta function δ h (s) has the following properties [START_REF] Roma | An adaptive version of the immersed boundary method[END_REF]:

i,j,k δ h (x i,j,k -X)∆v i,j,k = 1 (3) i,j,k (x i,j,k -X)δ h (x i,j,k -X)∆v i,j,k = 0. (4)
Reproducing conditions [START_REF] Lai | An immersed boundary method with formal second-order accuracy and reduced numerical viscosity[END_REF][START_REF] Beyer | Analysis of a one-dimensional model for the immersed boundary method[END_REF] are easily met on an uniform Cartesian grid by a number of regularized delta functions available in the literature. It will be shown later that correction terms must be introduced in the delta approximant if (1-2) are to be satisfied on an arbitrary underlying mesh.

In this paper we will follow the ideas originated by Liu et al. [START_REF] Liu | Generalized multiple scale reproducing kernel particle methods[END_REF] to build locally regularized functions that verify any number of integral conditions. These local approximants will be used both for interpolating the velocity field and for spreading the singular force field in the framework of a pressure-correction scheme for the incompressible Navier-Stokes equations. We will also demonstrate that the method conserves the order of accuracy of the spatial discretization both in Cartesian (uniform or non-uniform) meshes and in curvilinear ones. It is worth mentioning that, when higher order reproducing conditions are required, the approximant tends toward an exact delta function, and the resolution of the boundary becomes sharper.

The idea of using modified window functions to extend the immersed-boundary method to a general grid system is not new. Indeed it was originally proposed by Zhang et al. [START_REF] Zhang | Immersed finite element method[END_REF] in the context of the finite element method (i.e., the immersed finite element method). Later on, it was extended to a number of applications by several authors [START_REF] Wang | Extended immersed boundary method using FEM and RKPM[END_REF][START_REF] Liu | Immersed electrokinetic finite element method[END_REF][START_REF] Liu | Immersed finite element method and its applications to biological systems[END_REF]. Here we formulate the method in a finite difference/ finite volume context, and most importantly, we clarify the necessary conditions to be met for assembling discrete interpolation and spreading operators that are genuinely dual.

In the following, we will first present the mathematical formulation of the proposed immersed-boundary method and describe the numerical model. We then present the results of the calculation of two-and three-dimensional flows, using both Cartesian (uniform and non-uniform) and curvilinear grids. We will finally draw some conclusions and make recommendations for future work.

Mathematical formulation

The time advancement method

We start by introducing the general framework of a time-discretized formulation based on a fractional step method [START_REF] Chorin | Numerical solution of Navier-Stokes equations[END_REF][START_REF] Temam | Sur l'approximation de la solution des équations de Navier-Stokes par la méthode des pas fractionnaires (I)[END_REF] for the incompressible Navier-Stokes equations. We focus on the pressure-correction method [START_REF] Van Kan | A second-order accurate pressure-correction scheme for viscous incompressible flow[END_REF] with the only aim of introducing the general idea of the proposed immersed-boundary formulation. More details on the actual numerical implementation, including the temporal and spatial discretizations used for each numerical experiment will be given later.

The standard sequence used to solve the incompressible Navier-Stokes equation by a fractional step method is:

u * -u n ∆t = -N l (u n , u n-1 ) -Gφ n-1 + 1 Re L(u * , u n ) (5) Lφ = 1 ∆t Du * (6)
u n+1 = u * -∆tGφ n , (7)
where u * is the predicted (non-solenoidal) velocity field, u n is the divergence-free velocity field at time-step n, ∆t is the time step, N l is the discrete non-linear operator, G and D are, respectively, the discrete gradient and divergence operators, L is the discrete Laplacian, φ is a projection variable (related to the pressure field). The operators include coefficients that are specific to the selected time scheme. Note that we assume explicit integration of the advective term, and either explicit or implicit integration of the viscous diffusion, following the technique most often used for the direct or large-eddy simulation of turbulent flows. Following Uhlmann [START_REF] Uhlmann | An immersed boundary method with direct forcing for the simulation of particulate flows[END_REF], the sequence above is modified to impose the desired boundary values on the embedded geometry. The time advancement of the momentum equations is carried out in two stages: First, a fully explicit step analogous to (5) is performed, without any constraint on the embedded boundary:

u * = u n -∆t N l (u n , u n-1 ) -Gφ n-1 + 1 Re L(u n ) (8)
The velocity field obtained from ( 8) is then interpolated onto the embedded geometry Γ, which is discretized through a number of Lagrangian marker points with coordinates

X k : U * (X k , t n ) = I(u * ); (9)
The form of interpolation operator I will be specified later. The values of U * (X k , t n ) are used to determine a distribution of singular forces along Γ, that restore the prescribed boundary values U Γ (X k , t n ) on Γ:

F * (X k , t n ) = U Γ (X k , t n ) -U * (X k , t n ) ∆t . ( 10
)
The singular force field defined over Γ is then transformed by a spreading operator C into a volume force-field defined on the mesh points x i,j,k :

f * (x i,j,k , t n ) = C [F * (X k , t n )] . (11)
The regularized force given by ( 11) is now used as ( 5) is solved again (in explicit or implicit form):

u * -u n ∆t = -N l (u n , u n-1 ) -Gφ n-1 + 1 Re L(u * , u n ) + f * (12)
Finally, the algorithm completes the time step, with the usual solution of the pressure Poisson equation and the consequent projection step [START_REF] Feng | The immersed boundary-lattice Boltzmann method for solving fluidparticles interaction problems[END_REF].

In the next Section we will discuss the key points of the algorithm presented above. They are the definitions of the interpolation operator I, and of the convolution kernel involved in C.

Interpolation and Convolution: Dirac's delta approximants

We propose to use the Reproducing Kernel Particle Method (RKPM) to define interpolation and spreading operators that satisfy the conservation properties (1-2). We begin this subsection by outlining the aspects of the method that have a direct bearing on the definition of the interpolation and spreading operators that we propose. A complete review and analysis of the technique can be found in [START_REF] Liu | Generalized multiple scale reproducing kernel particle methods[END_REF][START_REF] Liu | Reproducing kernel particle methods[END_REF], while applications of the RKPM embedded domain technique in a finite element framework are illustrated in [START_REF] Zhang | Immersed finite element method[END_REF].

The approximation f a (x) of the value of a given smooth function f (s) at point x ∈ Ω can be expressed as a kernel approximation:

f a (x) = Ω w d (x -s)f (s)ds ( 13
)
where w d is the kernel (or weighting) function, and the subscript indicates that the kernel depends on an additional parameter d, the dilation parameter. The non-negative kernel function is assumed to be of compact support (i.e., nonzero in a subdomain Ω I of Ω and zero outside in Ω\Ω I ). The dilation parameter d determines the dimensions of the support Ω I . Note that if the kernel function is the delta function, f a (x) = f (x) and the function is reproduced exactly. Roma et al. [START_REF] Roma | An adaptive version of the immersed boundary method[END_REF] proposed a discrete approximation of the delta function

w d (r) =      1 6 5 -3|r| --3(1 -|r|) 2 + 1 0.5 ≤ |r| ≤ 1.5 1 3 1 + √ -3r 2 + 1 |r| ≤ 0.5 0 otherwise (14)
(where r = (xs)/d) that satisfies the following properties:

1. w d (r) is continuous ∀r ∈ R; 2. w d (r) = 0 if |r| ≥ 1.5; 3. l w d (r -l) = 1, ∀l ∈ N; 4.
l (rl)w d (rl) = 0, ∀r, l;

5.

l [w d (rl)] 2 = 1/2, ∀r, l.

Since the last three properties involve the integers l, the conservation properties (3-4) can be met by a function interpolated using ( 13) and ( 14) only if the nodes are equispaced.

To extend this approach to a non-uniform sequence of nodes, following Liu et al. [START_REF] Liu | Reproducing kernel particle methods[END_REF], we introduce a modified window function

w d (x -s) = n i=0 b i (y, d)(x -s) i w d (x -s), (15)
where the unknown polynomial coefficients b i (y, d) are determined by imposing the reproducing conditions

m 0 (x) = Ω w d (x -s) ds = 1 ( 16
)
m i (x) = Ω (x -s) i w d (x -s) ds = 0 (i = 1, ..., N ), (17)
which are the continuous equivalent of the third and fourth conditions met by w d (xs) on an uniform mesh.

It is interesting to observe that conditions [START_REF] Moin | Advances in large eddy simulation methodology for complex flows[END_REF][START_REF] Roman | An improved immersed boundary method for curvilinear grids[END_REF] imply the exact representation of the elements of the canonical polynomial base 1, x, x 2 , ... . For i = 0, and 1, in fact conditions ( 16) and ( 17) yield:

Ω w d (x -s) ds = 1 ( 18
) Ω s w d (x -s) ds = x; (19)
For i = 2, then, we have

x 2 Ω w d (x -s) ds + Ω s 2 w d (x -s) ds -2x Ω s w d (x -s) ds = 0 (20)
and, therefore,

x 2 = Ω s 2 w d (x -s) ds. (21)
Following this procedure it is possible to show, by induction, that the moment conditions [START_REF] Moin | Advances in large eddy simulation methodology for complex flows[END_REF][START_REF] Roman | An improved immersed boundary method for curvilinear grids[END_REF] are equivalent to the polynomial-reproducing conditions

x i = Ω s i w d (x -s) ds, i = 0, 1, ..., N. (22)
If we insert in [START_REF] Van Kan | A second-order accurate pressure-correction scheme for viscous incompressible flow[END_REF] the definition of the polynomial correction [START_REF] Vanella | A moving-least-squares reconstruction for embedded-boundary formulations[END_REF], it is possible to write the conditions on the modified moments m i (x) in terms of the original moments

m i (x) = Ω (x -s) i w d (x -s) ds = δ i0 (where δ ij is Kronecker's delta): m 0 (x) = Ω w d (x -s) ds = n i=0 b i (x)m i (x) = 1 m 1 (x) = Ω (x -s) w d (x -s) ds = n i=0 b i (x)m i+1 (x) = 0 • • • = • • • m j (x) = Ω (x -s) j w d (x -s) ds = n i=0 b i (x)m i+j (x) = 0 • • • = • • • m N (x) = Ω (x -s) N w d (x -s) ds = n i=0 b i (x)m i+N (x) = 0.
which results in a symmetric linear system

      m 0 m 1 • • • m N • • • • • • • • • • • • m j m j+1 • • • m N +j • • • • • • • • • • • • m N m N +1 • • • m 2N             b 0 • b j • b N       =       1 • 0 • 0       (23)
for the unknown polynomial coefficients b i , i = 0, 1, ..., N .

Modified window functions that verify the reproducing conditions in higher dimensions can be obtained by imposing the exact representation of a complete polynomial basis. Note that in 1D to ensure the integral conservation of the force and of its moments, one needs all the linear combinations of {1, x} to be exactly represented; in 2D, the linear combinations of {1, x, y} and, in 3D, the ones related with {1, x, y, z}. The exact representation of higher degree polynomials provides a sharper definition of the boundary [START_REF] Vanella | A moving-least-squares reconstruction for embedded-boundary formulations[END_REF]. Of course, the exact representation of higher order polynomials should be also considered if the underlying order of the Eulerian scheme is increased.

In 2D or 3D problems the mother window function can be given as a Cartesian product of ( 14) with itself:

w δ,η (x -s, y -t) = w δ (x -s) w η (y -t) (24)
or

w δ,η,σ (x -s, y -t, z -v) = w δ (x -s) w η (y -t) w σ (z -v) (25)
Here, δ, η and σ are again the dilatation parameters in the three coordinate directions.

Next we look for corrected window functions

w δ,η (x -s, y -t) = b 0 + (x -s)b 1 + (y -t)b 2 + (x -s)(y -t)b 3 +(x -s) 2 b 4 + (y -t) 2 b 5 w δ,η (x -s, y -t) (26)
(where b i = b i (δ, η, x, y) for i = 0, ..., 5) in 2D, and similarly,

w δ,η,σ (x -s, y -t, z -v) = b 0 + (x -s)b 1 + (y -t)b 2 + (z -v)b 3 + (x -s)(y -t)b 4 +(y -t)(z -v)b 5 + (z -v)(x -s)b 6 + (x -s) 2 b 7 +(y -t) 2 b 8 + (z -v) 2 b 9 w δ,η,σ (x -s, y -t, z -v) (27)
in 3D (where b i = b i (δ, η, σ, x, y, z) and i = 0, . . . , 10). By imposing the exact representation of the members of the polynomial basis, a symmetric linear system M 2D/3D b = e 1 is obtained. Where

M 2D =    m 0,0 m 1,0 m 0,1 m 1,1 m 2,0 m 0,2 m 1,0 m 2,0 m 1,1 m 2,1 m 3,0 m 0,2 m 0,1 m 1,1 m 0,2 m 1,2 m 2,1 m 0,3 m 1,1 m 2,1 m 1,2 m 2,2 m 3,1 m 2,3 m 2,0 m 3,0 m 2,1 m 3,1 m 4,0 m 2,2 m 0,2 m 1,2 m 0,3 m 1,3 m 2,2 m 0,4    (28)
and

M 3D =        m 0,0,0 m 1,0,0 m 0,1,0 m 0,0,1 m 1,1,0 m 0,1,1 m 1,0,1 m 2,0,0 m 0,2,0 m 0,0,2 m 1,0,0 m 2,0,0 m 1,1,0 m 1,0,1 m 2,1,0 m 1,1,1 m 2,0,1 m 3,0,0 m 1,2,0 m 1,0,2 m 0,1,0 m 1,1,0 m 0,2,0 m 0,1,1 m 1,2,0 m 0,2,1 m 1,1,1 m 2,1,0 m 0,3,0 m 0,1,2 m 0,0,1 m 1,0,1 m 0,1,1 m 0,0,2 m 1,1,1 m 0,1,2 m 1,0,2 m 2,0,1 m 0,2,1 m 0,0,3 m 1,1,0 m 2,1,0 m 1,2,0 m 1,1,1 m 2,2,0 m 1,2,1 m 2,1,1 m 3,1,0 m 1,3,0 m 1,1,2 m 0,1,1 m 1,1,1 m 0,2,1 m 0,1,2 m 1,2,1 m 0,2,2 m 1,1,2 m 2,1,1 m 0,3,1 m 0,1,3 m 1,0,1 m 2,0,1 m 1,1,1 m 1,0,2 m 2,1,1 m 1,1,2 m 2,0,2 m 3,0,1 m 1,2,1 m 1,0,3 m 2,0,0 m 3,0,0 m 2,1,0 m 2,0,1 m 3,1,0 m 2,1,1 m 3,0,1 m 4,0,0 m 2,2,0 m 2,0,2 m 0,2,0 m 1,2,0 m 0,3,0 m 0,2,1 m 1,3,0 m 0,3,1 m 1,2,1 m 2,2,0 m 0,4,0 m 0,2,2 m 0,0,2 m 1,0,2 m 0,1,2 m 0,0,3 m 1,1,2 m 0,1,3 m 1,0,3 m 2,0,2 m 0,2,2 m 0,0,4        (29)
Here, b is the array of unkown polynomial coefficients and e 1 is the first unit vector of the canonical basis in R 6 in 2D ( R 10 in 3D). The elements of the matrices are:

m i,j = Ω I (x -s) i (y -t) j w δ,η (x -s, y -t) ds dt (30)
in 2D, and

m i,j,k = Ω I (x -s) i (y -t) j (z -v) k w δ,η,σ (x -s, y -t, z -v) ds dt dv (31)
in 3D. The choice of the support Ω I and the solution procedure required to ensure that the system [START_REF] Liu | Reproducing kernel particle methods[END_REF] is not singular will be described in the following Subsection.

Interpolation and convolution: discrete approach

We now discuss the implementation of the IB technique based on RKPM in a discrete, generalized-coordinate finite-difference or finite-volume code. We limit ourselves to the 2D case with a second-order correction polynomial, since the 3D extension is straightforward.

First, the embedded-boundary curve must be discretized into a number of nodes X I , I = 1, . . . , N e . Around each node X I we define a rectangular cage Ω I that contains at least three nodes of the underlying mesh in each direction, while at the same time minimizing the number of grid nodes contained in the cage. The cage and underlying mesh are sketched in Figure 1. Following the definition [START_REF] Uhlmann | An immersed boundary method with direct forcing for the simulation of particulate flows[END_REF], the edges of the rectangle measure 3δ in x and 3η in y (recall that δ and η are the dilation parameters in the coordinate directions). Including at least three nodes in each direction in the cage is required to avoid a singular moment matrix when considering second-order polynomials for the correction of the original window function.

To determine the dimensions of the rectangular support centred in X I , first we look for the grid node closest to X I , x î, = (x î, , y î, ); next, we consider (in a structured mesh) the set of nodes neighboring x î, , N I = {x î+k,+l }, for k and l = -1, 0, 1; we then evaluate h ± x and h ± y as:

       h + x (X I ) = max {|x i,j -x i-1,j | : x i,j , x i-1,j ∈ N I } , h - x (X I ) = min {|x i,j -x i-1,j | : x i,j , x i-1,j ∈ N I } , h + y (X I ) = max {|y i,j -y i,j-1 | : y i,j , y i,j-1 ∈ N I } , h - y (X I ) = min {|y i,j -y i,j-1 | : y i,j , y i,j-1 ∈ N I } . (32)
Based on these values, we define the length of the edges of the rectangle (3δ I and 3η I ) through the local dilation factors:

δ I = 5 6 h + x (X I ) + 1 6 h - x (X I ) + ǫ x (X I ) (33)
η I = 5 6 h + y (X I ) + 1 6 h - y (X I ) + ǫ y (X I ) (34)
where ǫ x (X I ) and ǫ y (X I ) are small fractions of the local mesh spacing, and are added to avoid the support boundary touching some of the support nodes (to be defined later); in this case the window function would be zero at those nodes, making the discretized moment matrix singular. Finally, a set of mesh nodes that fall within the cage is sought:

S I = x i,j : |x î, -x i,j | < 3 2 δ I and |y î, -y i,j | < 3 2 η I . (35)
We have verified that with this particular choice, when the underlying mesh is reasonably smooth, the set of nodes within the support is at least nine almost everywhere (27 in 3D).

The elements of the moment matrix (30) must be evaluated numerically to assemble the local window function centered in X I . We approximate the entries in the moment matrix relative to node X I = (X I , Y I ) using the mid-point quadrature rule:

m I i,j = k,l∈S I (x k,l -X I ) i (y k,l -Y I ) j w δ I ,η I (x k,l -X I , y k,l -Y I ) ∆A k,l , (36)
where ∆A k,l is the area of the cell centered at x k,l . The extension of such an approximation to the 3D case ( 31) is trivial. Once the discrete moment matrix is assembled for each Lagrangian node X I , the coefficients of the correction polynomials are found by solving, at each Lagrangian point, the symmetric linear system: M I b I = e 1 , for I = 1, ..., N e . M I are the discrete equivalent of ( 28) or (29) (i.e., 6 × 6 in a 2D case and 10 × 10 in 3D, with a second-order correction polynomial). Due to the very low values that the window function may take at the nodes close to the boundary of the support cage, the moment matrix may become ill-conditioned. This problem can be avoided by rescaling the linear system, solving the equivalent one

H I M I (H I ) -1 b I = e 1 ,
where

H I = diag 1, 1 δ I , 1 η I , 1 δ I η I , 1 δ 2 I , 1 η 2 I (37)
(the 3D form of H I is a trivial extension of the one above) in two stages:

H I M I c I = e 1 , and
H I c I = b I . ( 38
)
The coefficient matrix of the first linear system can be equivalently obtained by normalizing the distances (x k,l -X I ) and (y k,l -Y I ) appearing in [START_REF] Tritton | Experiments on the flow past a circular cylinder at low Reynolds numbers[END_REF] with the dilation parameters (δ I and η I ). The matrix product that follows in [START_REF] Lu | Calculation of the timing of vortex formation from an oscillating cylinder[END_REF] is needed to undo the scaling.

The methodology developed so far allows the definition of a localized window function w δ I ,η I (x -X I ) to be used in the convolution integrals symbolically introduced in ( 9) and [START_REF] Kim | An immersed-boundary finite-volume method for simulations of flow in complex geometries[END_REF]. In particular, given a component of the velocity field u i (x, y) known at the mesh nodes x k,l ∈ S I , the interpolated value at node X I on the embedded line can be approximated by:

U i (X I ) = I(u i ) = k,l∈S I u i (x l,k ) w δ I ,η I (x k,l -X I , y k,l -Y I ) ∆A k,l , (39)
having used the same quadrature rule as in [START_REF] Tritton | Experiments on the flow past a circular cylinder at low Reynolds numbers[END_REF]. Once the force component F i (X I ) is found from [START_REF] Fadlun | Combined immersed-boundary finitedifference methods for three-dimensional complex flow simulations[END_REF], the distribution of the singular forces over the mesh nodes can be obtained using a discrete counterpart of (11) as convolution operator. This operator can be determined by using a quadrature formula over a strip surrounding Γ:

f i (x k,l ) = C(F i ) = Ne I=1 F i (X I ) w δ I ,η I (x k,l -X I , y k,l -Y I ) ε I ∆s I ( 40
)
where ∆s I is length of the arc joining X I+1/2 to X I-1/2 , and ε I is a characteristic strip-width related to the local dilation coefficients of the window function w δ I ,η I .

To determine the correct value of ε I we start by considering the value of the force on the markers obtained by interpolation from the nodes of the underlying grid:

F i (X I ) = k,l∈S I f i (x l,k ) w δ I ,η I (x k,l -X I , y k,l -Y I ) ∆A k,l (41)
Next we replace the values of f i (x l,k ) with those obtained from the spreading step (40):

F i (X I ) = k,l∈S I ∆A k,l w δ I ,η I (x k,l -X I ) × Ne K=1 F i (X K ) w δ K ,η K (x k,l -X K ) ε K ∆s K .
(42) Rearranging ( 42) the following conditions are obtained:

F i (X I ) = Ne K=1 a I,K ε K F i (X K ), I = 1, ..., N e (43)
where a I,K is the (discrete) integral of the product of the I th and K th window functions over the support of the former one multiplied by ∆s K :

a I,K = ∆s K k,l∈S I w δ I ,η I (x k,l -X I ) w δ K ,η K (x k,l -X K ) ∆A k,l (44)
In matrix notation, the system (43) can be written as:

[A diag(ε 1 , . . . , ε Ne )] F i = F i ; (45)
By requiring that the local width ε i is independent of the actual force distribution F i , the condition det [A diag(ε 1 , . . . , ε Ne ) -I] = 0 (46) (were I is the identity matrix) is obtained. The array ε that verifies the constraint (46) can be found by solving the linear system:

A ε = 1 (47) with 1 = (1, 1, . . . , 1) T .
It is worth noting that, by virtue of ( 41) and ( 44), a I,K /∆s K can be interpreted as the value of w δ K ,η K (x -X K ) evaluated at x = X I . Therefore, each I th equation in (47) is a condition imposing that the sum of all the K window functions evaluated at node X I , weighted by ε K ∆s K , is one. In other words, the weighted window functions ε K ∆s K w δ K ,η K (x -X K ) used to carry out the spreading step, must satisfy a partitionof-unity condition [START_REF] Babuška | The partition of unity method[END_REF].

The conditioning of matrix A depends on the ratios between the distances of the Lagrangian nodes ∆s K and the dilation factors δ I and η I , or equivalently the ratio between ∆s K and the local Eulerian grid size. To determine an optimal criterion for the selection of the number of Lagrangian nodes to be used to discretize an embedded contour ). The first column contains the number of Lagrangian markers uniformly distributed along Γ. The second and third columns show the values of the two extreme eigenvalues of A: note that the condition number (λmax/λ min ) grows from 6.67 up to 7.01 × 10 5 in the worst case. The fourth column shows the error after the spreading-interpolation procedure. The last two columns give the average value of the solution of system (47) and its root-meansquare. As a reference value, the interpolation error (only from Eulerian to Lagrangian mesh) for the case ∆s/∆x = 0.9873 is approximately 0.0273.

on a given Eulerian mesh we have carried out numerical experiments by considering a circle Γ embedded in a uniform mesh (with the same spacings in the x and y directions) and in a non-uniform one. The circle is discretized with equispaced Lagrangian nodes in the former case and non-equispaced nodes in the latter. We have defined a test function g(x, y) = sin πx cos πy to be used to measure the error induced by the spreading operation. We then evaluate the infinity norm of the difference between g(X I , Y I ) and g IC (X I , Y I ), i.e., the values of the test function at the same nodes obtained by first spreading g(X I , Y I ) onto the Eulerian mesh and then re-interpolating the discrete values over Γ using [START_REF] Williamson | Defining a universal and continuous Strouhal-Reynolds number relationship for the laminar vortex shedding of a circular cylinder[END_REF]:

g IC (X I , Y i ) = I {C [g(X I , Y I )]} . ( 48
)
Table 1 presents a summary of the results obtained for the case of the circle embedded onto an uniform Cartesian Eulerian mesh. The maximum distance between the equispaced Lagrangian nodes is small enough that the function can be represented accurately. The error increases when ∆s becomes smaller than the Eulerian mesh spacing. Furthermore, when the Lagrangian nodes become too close to each other, the solution of (47) becomes oscillatory (as shown by min(λ) → 0, while max(λ), which depends on the dimension of the support, remains approximately constant). The deterioration of the solution observed as the number of Lagrangian nodes is increased depends on the fact that we are re-interpolating a function that had been sampled on Γ at a higher resolution, resulting in a Gibbs-like phenomenon. A similar behavior was observed when considering a non-uniform Cartesian mesh. We conclude that the spacing of Lagrangian nodes should be approximately equal to the local grid size. If this requirement is satisfied, the linear system (47) is well-conditioned and the resulting solution ε(s) turns out to be smooth and positive at all Lagrangian nodes. Note that the positiveness of ε(s) is consistent with the interpretation of the discrete spreading [START_REF] Taneda | Experimental investigation of the wake behind a sphere at low Reynolds numbers[END_REF] as an integral over a finite stripe.

When the underlying grid spacing is coarser than the Lagrangian spacing between neighboring markers, it is still possible to obtain an acceptable solution for ε(s) either using a singular-value-decomposition when solving (47) (by filtering our the smallest singular values, or by solving a regularized version of the given linear system). We will not discuss further this issue, since the criterion for the distribution of the Lagrangian markers given above appears sufficient to impose the desired boundary value on the immersed boundary accurately. Finally, it is worth noting that the conservation properties ( 1) and ( 2) are verified independently of the value assigned to ε, whose role is only to enforce a partition-of-unity condition. Furthermore, the actual entries in matrix A do not need to be computed explicitly, since the action of A over a given vector ε p is simply a sequence of a spreading action of the unity function over Γ, followed by an interpolation according to (40) (using ε p ) and [START_REF] Williamson | Defining a universal and continuous Strouhal-Reynolds number relationship for the laminar vortex shedding of a circular cylinder[END_REF], respectively. As shown above, if the Lagrangian spacing is comparable to the local Eulerian grid size, any Krylov-type iterative method, applied to (47), converges in few iterations (typically 4-6) without any preconditioning.

The Navier-Stokes solver

The immersed-boundary method discussed so far has been implemented in a curvilinear finite-volume solver for the incompressible Navier-Stokes equations [START_REF] Silva Lopes | Simulations of isotropic turbulence using a non-orthogonal grid system[END_REF]. The timeadvancement follows the procedure roughly described above in Equations (5-7).

In particular, the equations are discretized on a non-staggered grid system using a curvilinear finite volume code. The method of Rhie and Chow [START_REF] Rhie | Numerical study of the turbulent flow past an airfoil with trailing edge separation[END_REF] is used to avoid pressure oscillations. Both convective and diffusive fluxes are approximated by secondorder central differences. A second-order semi-implicit fractional-step procedure [START_REF] Kim | Application of a fractional step method to incompressible Navier-Stokes equations[END_REF] is used for the temporal discretization. The Crank-Nicolson scheme is used for the temporal discretization of wall-normal diffusive terms, and the second-order Adams-Bashforth scheme for all the other terms. Fourier transforms are used to reduce the three-dimensional Poisson equation into a series of two-dimensional Helmoltz equations in wavenumber space, which are then solved iteratively using the biconjugate gradient stabilized (BiCGStab) method. The code is parallelized using the MPI message-passing library and the domain-decomposition technique, and has been widely tested [START_REF] Silva Lopes | Simulations of isotropic turbulence using a non-orthogonal grid system[END_REF][START_REF] Silva Lopes | Large-eddy simulation of the flow in an S-duct[END_REF][START_REF] Radhakrishnan | Reynolds-averaged and large-eddy simulations of turbulent non-equilibrium flows[END_REF][START_REF] Radhakrishnan | Wall-modeled large-eddy simulations of flows with curvature and mild separation[END_REF] in simulations of turbulent flows using curvilinear, body-fitted grids.

Results

In this section the accuracy, convergence and robustness of the RKPM-based immersedboundary method proposed here will be investigated. We consider first the laminar flow around a circular cylinder; we then examine the flow over a two-dimensional hill to show the implementation of the scheme on non-uniform and non-orthogonal grids. Then we implement and test the scheme in a three-dimensional case: the flow around the sphere is simulated in a range of Reynolds numbers that covers both steady and unsteady cases. Finally, we consider a case in which the body is moving, by examining the impulsive start of a two-dimensional flat plate orthogonal to its direction of motion.

Flow around a circular cylinder

The steady flow around a circular cylinder is considered at two Reynolds numbers (based on the free-stream velocity U ∞ and the cylinder diameter D), Re D = 30 and Re D = 185. At the low Reynolds number the flow remains steady, while at Re D = 185 periodic vortex shedding is expected; in this case the results depend critically on the accurate reproduction of the vorticity field in the vicinity of the cylinder, which, in turn, depends on a correct prescription of pressure and shear forces at the solid boundary.

The flow around the cylinder is simulated using a domain of dimensions [-9D, 40D] in the streamwise (x) direction and [-17D, 17D] in the vertical (y) direction; the center of the cylinder is at (0, 0). The dimensions of the domain are comparable to those used in other numerical studies [START_REF] Uhlmann | An immersed boundary method with direct forcing for the simulation of particulate flows[END_REF][START_REF] Vanella | A moving-least-squares reconstruction for embedded-boundary formulations[END_REF]. A steady uniform velocity is assigned at the inlet plane, while at the outlet plane convective boundary condition are specified. Free-slip conditions are applied on the top and the bottom of the domain.

The grid is uniform near the cylinder, in the region -0.6D ≤ x ≤ D and -D ≤ y ≤ D. Outside this region, it is stretched with a stretching ratio of less than 1.02. Four different grids (with ∆x = ∆y = 0.0025D, 0.005D, 0.01D and 0.02D in the uniform region) are considered with 1260, 629, 316 and 153 uniformly spaced Lagrangian points on the boundary of the cylinder. The results for the finest grid are considered accurate, and the L 2 and L ∞ norms of the error obtained on the coarser grids are calculated and shown in Figure 2. The results demonstrate the second-order accuracy of the method. The simulation results are also compared with the experimental data. The most important physical feature of this flow is the presence of a recirculating region in the wake of the cylinder. The important parameters associated with the wake are shown in Figure 3; l is the length of the wake, a is the distance from the cylinder to the center of the wake vortex, b is the distance between the centres of the wake vortices, and θ is the angle of separation measured from x-axis. Another important flow parameter is the drag coefficient, C D = 2D/ρU 2 ∞ D 2 (where D is the drag force). Since the interpolation and spreading operators conserve the force, the drag force D can be calculated directly from the summation of the forces at all Lagrangian points:

D = Ne l=1 F * (X l ) • e 1 ∆s l , ( 49
)
where e 1 is the unit vector in the flow direction. The streamlines and vorticity contours for this flow are shown in Figure 4. The wake parameters computed from the simulation are in good agreement with the experimental data, as shown in Table 2. To illustrate the effect of the forcing in the vicinity of the immersed body, in Figure 5 we show profiles of pressure p and of the u and v velocity components near the cylinder along three lines: a horizontal one, a vertical one, and one inclined by 45 o from the horizontal (and oriented upstream). Note that the region from [-0.5, 0.5] along the x and y-directions is within the cylinder. We observe some fluctuations in the velocity component normal to the boundary [u at the stagnation point, Figure 5(a), and v at the top of the cylinder, Figure 5(e)], while the velocity component parallel to the boundary goes to zero smoothly. The oscillations, however, are confined to the three points of the support, where the flow is expected to be unphysical, and are not corrupting the rest of the flow field. Also, notice that the interpolation of the velocities to the Lagrangian point results in an exact enforcement of the no-slip condition even when oscillations are present. The pressure shows a large jump, again restricted to the support points, and varies smoothly away from the body. Figures 5 (g), (h) and (i) show the velocities u, v and pressure p variation along a line at 45 o from the stagnation point.

There are no fluctuations in the velocity components, and pressure and velocities go to zero at the cylinder wall smoothly. To further quantify the behavior of the flow in the vicinity of the cylinder boundary a third-order spline extrapolation of the velocity near the boundary is also plotted on the curves (as a dash-dot-dot line) in Figure 5(a) and (e). The computation gives zero velocity at the Lagrangian points; however, the extrapolated velocity is not zero. The absolute difference between the extrapolated velocity and the Lagrangian-point value (i.e., the no-slip condition) reduces at second-order rate as the grid is refined, as shown in Figure 6. Also notice that there is motion inside the cylinder (since we do not enforce a zero-velocity condition there); the inner flow is, however, decoupled from the outer one.

We then examined a higher Re D case, to verify that the unsteadiness of the flow is captured correctly. The domain and the boundary conditions for the higher Reynolds number case were similar to the low-Re D case. Only two grids were considered, with ∆x = ∆y = 0.01D and 0.005D, respectively, in the region of interest (-0.6D ≤ x ≤ D and D ≤ y ≤ D). 3. The results obtained with the present method are well within the range of values obtained by other researchers.

Flow over a two-dimensional hill

An important property of the current forcing scheme is its applicability to non-uniform and non-orthogonal grids without loss of accuracy. To demonstrate this property a two dimensional flow over a small hill in a channel is considered here with a Reynolds number of Re = 600 based on the channel height H. The domain size is 35H × H. The shape of the hill is defined by y h = 0.15H sin 2 [π(x -2.0)]; thus, the maximum height is 15% of the channel height, its peak is at x = 2.5. The no-slip boundary condition is enforced on the top and bottom walls of the channel and a periodic boundary condition is assumed in the streamwise direction. The domain, however, is long enough that the flow returns to a fully developed state, with a parabolic profile, before the end of the computational domain.

Four different grid topologies, shown in Figure 8, are considered. First a nonorthogonal body-fitted grid is used, which conforms to the hill. In this case hill becomes

C D C rms L St Present ∆x = 0.005D
1.430 0.423 0.196 Present ∆x = 0.01D 1.509 0.428 0.199 Vanella and Balaras [START_REF] Vanella | A moving-least-squares reconstruction for embedded-boundary formulations[END_REF] 1.377 0.461 -Guilmineau and Queutey [START_REF] Guilmineau | A numerical simulation of vortex shedding from an oscillating circular cylinder[END_REF] 1.280 0.443 0.195 Lu and Dalton [START_REF] Lu | Calculation of the timing of vortex formation from an oscillating cylinder[END_REF] 1.310 0.422 0.195 Williamson [START_REF] Williamson | Defining a universal and continuous Strouhal-Reynolds number relationship for the laminar vortex shedding of a circular cylinder[END_REF] --0.193 an integral part of the domain boundary and direct no-slip boundary condition is applied. This is considered the reference case, with which the others are compared. The maximum grid size in the region 2H ≤ x ≤ 4H and 0 ≤ y ≤ 0.25H is 0.01H. This is the region of interest for this flow, since it covers the hill and the recirculation zone.

The second grid is also a non-orthogonal grid, but the hill does not coincide with a grid line, and the no-slip conditions on the hill itself are imposed by the immersed-boundary method. The maximum grid size in the region of interest for this case is also 0.005H. The third grid is a non-uniform orthogonal grid, clustered in y near the top and bottom of the hill, and in x near its leading and trailing edges. In the central part of the hill the grid is stretched in both directions. Here the grid size in the region of interest varies between 0.0025H and 0.0075H. The fourth grid is a uniform orthogonal grid, also with a grid size ∆x = ∆y = 0.005H in the region of interest. Except for the first case, the hill is defined as an immersed boundary with 650 uniformly distributed Lagrangian points along the geometry, and the no-slip condition on the hill is applied through the forcing described previously. The vorticity contours on the hill and in the recirculation region for the four grids are shown in Figure 9. The non-orthogonality and non-uniformity of the grid has no significant effect on the large-scale flow features. The vorticity contours for all the three grids with the IBM compare well with the body-fitted non-orthogonal grid. A quantitative comparison of the u and v velocity profiles is shown in Figure 10. The velocity profiles are compared at three locations, x = 2.5 (the peak), x = 2.75 (50% along the downward slope) and x = 3.2 (in the recirculation region). The u and v velocity profiles for all the four grids are in very good agreement with each other. The maximum error appears to be on the v velocity component, at the hill crest. Even there, the maximum difference between the velocity calculated using the boundary-fitted grid and the methods with the immersed boundary and the uniform grid is less than 7%, decreasing to 3% when non-uniform mesh is used, and the flow gradients are resolved better. Taneda [START_REF] Coutanceau | Experimental determination of the main features of the viscous flow in the wake of a circular cylinder in uniform translation. Part 1. Steady flow[END_REF] Present Calculations Fadlun et. al. [START_REF] Fadlun | Combined immersed-boundary finitedifference methods for three-dimensional complex flow simulations[END_REF] Zhang and Zheng [START_REF] Guilmineau | A numerical simulation of vortex shedding from an oscillating circular cylinder[END_REF] Johnson and Patel [START_REF] Lu | Calculation of the timing of vortex formation from an oscillating cylinder[END_REF] Figure 11: Length of the separation bubble at different Reynolds number.

Flow around a sphere

The examples discussed above show that the current scheme is able to simulate twodimensional flows with complex immersed geometries very well. The application of the scheme to three-dimensional immersed objects is tested next by simulating the flow around a sphere at low Reynolds numbers. A Cartesian domain [-4D, 10D]×[-4D, 4D]× [-4D, 4D] is used. The sphere is placed with its center at (0, 0, 0). A uniform grid of ∆x = ∆y = ∆z = 0.05D is used around the sphere; the grid is stretched away from the sphere. The sphere is defined using 1258 equally spaced Lagrangian points on its surface. A range of Reynolds numbers Re D from 30 to 100 is considered. The length of the separation bubble for these Reynolds numbers is compared with the data by Taneda [START_REF] Taneda | Experimental investigation of the wake behind a sphere at low Reynolds numbers[END_REF] presented in Batchelor [START_REF] Batchelor | An Introduction to Fluid Mechanics[END_REF] and with two other immersed-boundary techniques (by Fadlun et al. [START_REF] Fadlun | Combined immersed-boundary finitedifference methods for three-dimensional complex flow simulations[END_REF] and Zhang and Zheng [START_REF] Zhang | An improved direct-forcing immersed-boundary method for finite difference applications[END_REF]) in Figure 11, which shows that the prediction of the length of the separation bubble compares well with other immersed-boundary methods. Another important parameter associated with this flow is the coefficient of pressure C p = 2(pp ∞ )/ρU 2 ∞ , where p is the pressure at the sphere wall and p ∞ is pressure in undisturbed flow field. The C p from the present simulation is in very good agreement with the results of Fadlun et al. [START_REF] Fadlun | Combined immersed-boundary finitedifference methods for three-dimensional complex flow simulations[END_REF] and Zhang and Zheng [START_REF] Zhang | An improved direct-forcing immersed-boundary method for finite difference applications[END_REF] at Re D = 100, as shown in Figure 12.

The flow around a sphere at low Reynolds numbers remains symmetric and no shedding occurs. However, at a slightly higher Reynolds number the flow become asymmetric; at Re D = 280 vortex shedding is started (Johnson and Patel [START_REF] Johnson | Flow past a sphere up to a reynolds number of 300[END_REF]). We performed a calculation in the shedding region, at Re D = 300, using a Cartesian domain with dimensions [-4D, 11D] × [-4.5D, 4.5D] × [-4.5D, 4.5D]. The grid size in the vicinity of the sphere is ∆x = ∆y = ∆z = 0.018D, and is stretched away from the sphere. The sphere is defined using 10,000 equally spaced Lagrangian points on its surface.

The mean streamwise velocity along the x-axis, U avg , in the wake region is compared 22 Present Calculations Fadlun et. al [START_REF] Fadlun | Combined immersed-boundary finitedifference methods for three-dimensional complex flow simulations[END_REF] Zhang and Zheng [START_REF] Guilmineau | A numerical simulation of vortex shedding from an oscillating circular cylinder[END_REF] Figure 12: Coefficient of pressure Cp at Reynolds number Re D = 100.

in Figure 13 with the data obtained by Johnson and Patel [START_REF] Johnson | Flow past a sphere up to a reynolds number of 300[END_REF] using a body-fitted grid.

There is good agreement between the present calculation and the reference results. The shedding from the sphere results in some unique coherent structures shown in Figure 14.

The vortices are visualized as isosurfaces of the second invariant of the velocity gradient tensor, Q = (u i,j u j,i )/2 and the frames are 20 time units apart. The vortices appear like hairpin vortices. The heads of the shed vortex form very close to the sphere in the wake region and its legs stretch as it moves away from the sphere. The legs are attached to the head of the next vortex. The head of the alternate shedded vortices are inclined in vertically opposite directions with respect to the x-axis.

Suddenly accelerated normal flat plate

As a last example, we consider the case of a moving body, an infinitesimally thin finite flat plate of height h, suddenly accelerated from rest to a constant velocity U o in the direction normal to its surface in a fluid at rest. The Reynolds number is Re h = U o h/ν = 1000. This value matches two simulations of this case by Mittal et al. [START_REF] Mittal | A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries[END_REF], who used of a direct-forcing method, and Koumoutsakos and Shiels [START_REF] Koumoutsakos | Simulations of the viscous flow normal to an impulsively started and uniformly accelerated flat plate[END_REF] who employed a vortex-particle method. The latter is particularly well suited for this problem, since, at least at early times, the vorticity is confined to a small subregion of the domain, thereby easing the computational requirements for the vortex simulations.

In our calculation, the grid is uniform in the direction of the plate motion, while in the normal direction a stretched grid has been employed (with a stretching ratio of 1.045) to cluster the nodes near the ends of the plate. The computational domain, 12h × 9h (in the streamwise and normal direction), has been discretized using 2550 × 260 points, giving a minimum grid spacing ∆x min = ∆y min = 0.0012h. The Lagrangian markers are clustered according to the local vertical mesh spacing. Figure 15 shows the evolution of the wake behind the plate at three times. Very good agreement with the results in the literature (see Figure 18 of [START_REF] Mittal | A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries[END_REF] and Figure 5 of [START_REF] Koumoutsakos | Simulations of the viscous flow normal to an impulsively started and uniformly accelerated flat plate[END_REF]) is found.

Figure 16 shows the temporal variation of the computed non dimensional bubble length (i.e., the length of the region of reverse flow s/h measured on the centerline behind the plate, normalized by plate height) obtained from the current simulations. Again, an excellent agreement is found with the results of Koumotsakos and Sheils [START_REF] Koumoutsakos | Simulations of the viscous flow normal to an impulsively started and uniformly accelerated flat plate[END_REF]. It is worth noting that this case demonstrates the ability of the forcing algorithm to handle infinitesimally thin bodies including the singular points at the two extremities, with no need to handle this configuration through the use of auxiliary ghost-cell.

Conclusions

We have presented a novel interpolation-spreading procedure in the context of immersedboundary type methods. The global algorithm follows the ideas originally introduced by Peskin [START_REF] Peskin | Flow patterns around heart valves: a numerical method[END_REF], and lately extended to rigid moving bodies by Uhlmann [START_REF] Uhlmann | An immersed boundary method with direct forcing for the simulation of particulate flows[END_REF] and Vanella and Balaras [START_REF] Vanella | A moving-least-squares reconstruction for embedded-boundary formulations[END_REF], among others. These approaches differ from the direct-forcing methodology introduced by [START_REF] Fadlun | Combined immersed-boundary finitedifference methods for three-dimensional complex flow simulations[END_REF], since the evaluation of the forcing function is done on the Lagrangian markers distributed on the embedded geometry, rather than on the Eulerian grid points surrounding it. In principle, the main advantages of this strategy compared to existing direct-forcing schemes is the ease and versatility of its implementation with almost any available solver for parabolic equations. This feature allows decoupling the computation of the required forcing from the computational grid itself. The forcing procedures proposed so far in the literature, however, do not fully exploit the theoretical advantages of the general methodology. The present contribution aims to develop a numerical method leading to a genuine plug-and-play module that can be applied to existing solvers based on any spatial discretization.

The method proposed here draws heavily from ideas developed for meshless methods, particularly from the Reproducing Kernel Particle Method (RKPM) [START_REF] Liu | Generalized multiple scale reproducing kernel particle methods[END_REF], its application as an immersed boundary tool in the finite element context [START_REF] Zhang | Immersed finite element method[END_REF], and from the partitionof-unity concept [START_REF] Babuška | The partition of unity method[END_REF], for the design of kernel functions that preserve the order of accuracy of the underlying spatial scheme. In particular, the main improvements with respect to other techniques are:

1. It can be applied to codes relying upon any spatial discretization (Cartesian nonuniform, structured body fitted and unstructured grids). 2. One can preserve the underlying order of the spatial scheme while improving the sharpness of the immersed boundaries, in principle, by selecting the order of the reproducing conditions. 3. The spreading and interpolation operators are dual of each other. 4. It can be applied to moving solid objects 5. The computational cost of the forcing procedure is a small fraction of the overall cost required by the original solver. 6. The forces and moments caused by the fluid are calculated in a straightforward manner by integrating the forcing field.

These qualities of the proposed algorithm have been demonstrated numerically applying the proposed forcing method in a curvilinear, finite-volume incompressible Navier-Stokes solver for steady and unsteady flow over two-dimensional and three-dimensional objects, both at rest and in motion. Cartesian uniform and non-uniform grids have been used, as well as curvilinear non-orthogonal meshes. In all cases the numerical results were in excellent agreement with data in the literature, showing the robustness, accuracy, and range of applicability of the proposed methodology.

No tests using unstructured grids have been performed in this study. Nevertheless, the method is inherently general and the only issues that should be faced for this class of discretizations would concern a strategy to determine the support of the window function and a quadrature rule consistent with a cell centered or cell vertex approach.

Figure 1 :

 1 Figure 1: Definition of the support cage. The dashed line is the embedded curve and the dashed rectangle is the support cage Ω I (centered about X I ). • Lagrangian markers on the embedded curve, grid points, grid points within the support.

Figure 2 :

 2 Figure 2: Grid convergence study. (a) L 2 and (b) L∞ norms of the error.

Figure 3 :

 3 Figure 3: Definition of different parameters of the wake.

Figure 4 :

 4 Figure 4: Streamlines and vorticity contours for the steady-state flow around a cylinder at Re D = 30. l/D a/D b/D θ C D Present 1.70 0.56 0.52 48.05 o 1.80 Coutanceau and Bouard [35] 1.55 0.54 0.54 50.00 o -Tritton [36] ----1.74

Figure 5 : 2 Figure 6 :

 526 Figure 5: Behavior of u, v, and p near the body. Stagnation point region: (a) u, (b) v, (c) p. Top of the cylinder: (d) u, (e) v and (f ) p. Along a line at 45 o from point of stagnation: (g) u, (h) v and (i) p.

Figure 7 Figure 7 :

 77 Figure 7: Flow around a cylinder at Re D = 185. (a) Time-series of lift and drag coefficients; (b) instantaneous vorticity contours and streamlines.

Figure 8 :Figure 9 :

 89 Figure 8: Grids for two-dimensional hill test. (a) Body-fitted non-orthogonal grid. (b) Non-orthogonal grid with immersed boundary. (c) Non-uniform orthogonal grid with immersed body. (d) Uniform orthogonal grid with immersed boundary. Every fourth grid line is shown except in (a), where every line is shown. Only part of the domain is shown.

Figure 10 :

 10 Figure 10: Comparison of velocity profiles at x = 2.5, x = 2.75 and x = 3.2 over the hill. (a) u, (b) v.Body fitted non-orthogonal grid; Non-orthogonal grid with IB; Non-uniform orthogonal grid with IB;Uniform orthogonal grid with IB.

Figure 13 :

 13 Figure 13: Mean streamwise velocity in the wake region.

Figure 14 :

 14 Figure 14: Visualization of the coherent structures shedding from the sphere at Re D = 300. Frames (a), (b) and (c) are 20D/U∞ time units apart.

Figure 15 :

 15 Figure 15: Spanwise vorticity contours at three dimensionless times tUo/h: (a) 0.5; (b) 1 and (c) 2. Positive vorticity corresponds to solid lines, negative values to dashed lines.

Figure 16 :

 16 Figure 16: Dimensionless length of the reversed-flow region behind the plate. Solid line and circles corresponds to the present results, × are values obtained by Koumotsakos and Shiels [45].

Table 1 :

 1 Spreading and interpolation error for a uniform, 2D Cartesian Eulerian grid covering the surface [0, 5] × [0, 5] and a circle of unitary radius with center at (2.5, 2.5

	Ne	∆s/∆x	min(λ)	max(λ) ||error||∞	ε	rms(ε)
	105	1.1818	1.48	9.85	0.0131	0.103	0.103
	115	1.0791	0.714	9.83	0.0131	0.103	0.103
	125	0.9921	0.0648	9.81	0.0175	0.102	0.103
	135	0.9209	7.85×10 -3	9.80	0.0235	0.103	0.116
	145	0.8577	2.72×10 -3	9.80	0.0349	0.103	0.197
	155	0.8024	1.11×10 -3	9.81	0.103	0.103	0.401
	165	0.7510	4.27×10 -4	9.81	0.175	0.104	1.30
	175	0.7115	1.61×10 -4	9.80	0.313	0.105	2.90
	185	0.6719	4.94×10 -5	9.80	0.486	0.106	4.46
	195	0.6364	1.39×10 -5	9.81	1.84	0.119	21.1

Table 2 :

 2

Comparison of wake parameters and drag coefficient for steady-state flow around a cylinder at Re D = 30 with experimental data.

Table 3 :

 3 Comparison of coefficients of drag and lift, and Strouhal number for the flow around a cylinder at Re D = 185.

Acknowledgements

This research was began while AP was visiting the Department of Mechanical and Materials Engineering at Queen's University. The computational support of the High Performance Computing and Virtual Laboratory, Queen's University site, is gratefully acknowledged. AP and JF acknowledge the financial support of the Spanish Ministry of Innovation and Science through the Consolider grant Supercomputación y e-Ciencia. UP acknowledges the support of the Natural Science and Engineering Research Council of Canada (NSERC) under the Discovery Grant program.


Zhuhong Shao

Huazhong Shu

Jiasong Wu

Zhifang Dong

Gouenou Coatrieux

Jean-Louis Louis Coatrieux

Double color image encryption using iterative phase retrieval algorithm in quaternion gyrator domain

published or not. The documents may come

Double color image encryption using iterative phase retrieval algorithm in quaternion gyrator domain

Introduction

The image encryption technique, by its capability to ensure the security of information transmission and communication, has attracted a growing attention since the double random phase encoding technique was proposed by Refregier and Javidi [START_REF] Refregier | Optical image encryption based on input plane and Fourier plane random encoding[END_REF] in 1995. In general, the input image is encrypted into an undistinguishable image and decrypted when received with the proper granted keys. In the past decades, various optical or digital encryption algorithms for gray images have been developed [START_REF] Refregier | Optical image encryption based on input plane and Fourier plane random encoding[END_REF][START_REF] Liu | Optical image encryption by cascaded fractional Fourier transforms with random phase filtering[END_REF][START_REF] Hennelly | Optical image encryption by random shifting in fractional Fourier domains[END_REF][START_REF] Nishchal | Securing information using fractional Fourier transform in digital holography[END_REF]. Some advanced techniques of optical encryption have emerged lately, such as the spatially incoherent image encryption system [START_REF] Zang | Optical image encryption with spatially incoherent illumination[END_REF], ptychography based encryption [START_REF] Shi | Optical image encryption via ptychography[END_REF] as well as asymmetric cryptosystem using random binary phase modulations [START_REF] Liu | Asymmetric cryptosystem using random binary phase modulation based on mixture retrieval type of Yang-Gu algorithm[END_REF]. Furthermore, a new and effective optical system using phase modulation with sparse representation was reported in [START_REF] Chen | Phase-modulated optical system with sparse representation for information encoding and authentication[END_REF], which achieves higher security by employing optical authentication method.

Color images remain however of major importance in many applications [START_REF] Swain | Color indexing[END_REF][START_REF] Gevers | Color-based object recognition[END_REF]. Zhang and Karim [START_REF] Zhang | Color image encryption using double random phase encoding[END_REF] were the first to introduce a single-channel encryption algorithm for color image using double random phase encoding in the Fourier domain. Multiple channels encryption algorithms for color image were meanwhile investigated by applying the existing encryption algorithms for gray image on the red, green and blue channels independently. They make use, for example, of Fresnel transform [START_REF] Chen | Optical color image encryption by wavelength multiplexing and lensless Fresnel transform holograms[END_REF], fractional Fourier transform [START_REF] Chen | Color information processing (coding and synthesis) with fractional Fourier transforms and digital holography[END_REF][START_REF] Sui | Single-channel color image encryption using phase retrieve algorithm in fractional Fourier domain[END_REF] or gyrator transform [START_REF] Abuturab | Color image security system using double random-structured phase encoding in gyrator transform domain[END_REF][START_REF] Sui | Color image encryption based on gyrator transform and Arnold transform[END_REF]. Based on the quaternion theory [START_REF] Hamilton | Elements of Quaternions[END_REF], recent works introduced color image encryption algorithms in a fully vectorial form. They rely on the discrete quaternion Fourier transform [START_REF] Wang | Double random-phase encryption based on discrete quaternion Fourier-transforms[END_REF] and the quaternion gyrator transform [START_REF] Shao | Quaternion gyrator transform and its application to color image encryption[END_REF], respectively. The interest of the quaternion gyrator transform is that the presence of rotation angle enlarges the keys space and thus enhances the system security.

In addition to the aforementioned encryption algorithms for single gray or color image, multiple-image encryption algorithm capable to fuse more than one target image into a single illegible image has attracted great interest [START_REF] Situ | Multiple-image encryption by wavelength multiplexing[END_REF][START_REF] Shi | Multiple-image hiding in the Fresnel domain[END_REF][START_REF] Hwang | Multiple-image encryption and multiplexing using a modified Gerchberg-Saxton algorithm and phase modulation in Fresnel-transform domain[END_REF][START_REF] Chang | Position multiplexing multiple-image encryption using cascaded phase-only masks in Fresnel transform domain[END_REF][START_REF] Chang | Wavelength multiplexing multiple-image encryption using cascaded phase-only masks in the Fresnel transform domain[END_REF][START_REF] Alfalou | Double random phase encryption scheme to multiplex and simultaneous encode multiple images[END_REF][START_REF] Deng | Multiple-image encryption using phase retrieve algorithm and intermodulation in Fourier domain[END_REF][START_REF] Wang | Fully phase multiple-image encryption based on superposition principle and the digital holographic technique[END_REF][START_REF] Sui | Multiple-image encryption based on phase mask multiplexing in fractional Fourier transform domain[END_REF][START_REF] Kong | Multiple-image encryption scheme based on cascaded fractional Fourier transform[END_REF][START_REF] Tao | Double image encryption based on random phase encoding in the fractional Fourier domain[END_REF][START_REF] Liu | Double image encryption based on iterative fractional Fourier transform[END_REF][START_REF] Liu | Triple image encryption scheme in fractional Fourier transform domains[END_REF][START_REF] Wu | Four-image encryption method based on spectrum truncation, chaos and the MODFrFT[END_REF][START_REF] Joshi | Color image encryption and decryption for twin images in fractional Fourier domain[END_REF][START_REF] Joshi | Fractional Fourier transform based image multiplexing and encryption technique for four-color images using input images as keys[END_REF][START_REF] Liu | Double image encryption by using iterative random binary encoding in gyrator domains[END_REF][START_REF] Li | Double-image encryption based on iterative gyrator transform[END_REF]. The superiority of multiple-image encryption relies on the minimization of storage requirements for secret data and the improvement in transmission efficiency. Although the previous encryption algorithms for single gray or color image can be directly used for successive encryption of multiple images, the cost and the complexity of such systems will be expected to significantly increase. The drawbacks of such approach have been discussed for the double random phase (DRP) encoding technique in [START_REF] Alfalou | Double random phase encryption scheme to multiplex and simultaneous encode multiple images[END_REF]. To surmount these deficiencies, advanced encryption algorithms for multiple images have been proposed these last years. By combining the iterative Fourier transform algorithm with the classical double random phase encryption system, Alfalou and Mansour [START_REF] Alfalou | Double random phase encryption scheme to multiplex and simultaneous encode multiple images[END_REF] proposed a two level encryption/decryption system to encode multiple images. Meanwhile, the iterative phase retrieval algorithm for gray or binary images associated or not with other techniques in Fresnel transform domain [START_REF] Shi | Multiple-image hiding in the Fresnel domain[END_REF][START_REF] Hwang | Multiple-image encryption and multiplexing using a modified Gerchberg-Saxton algorithm and phase modulation in Fresnel-transform domain[END_REF][START_REF] Chang | Position multiplexing multiple-image encryption using cascaded phase-only masks in Fresnel transform domain[END_REF][START_REF] Chang | Wavelength multiplexing multiple-image encryption using cascaded phase-only masks in the Fresnel transform domain[END_REF], fractional Fourier transform domain [START_REF] Sui | Multiple-image encryption based on phase mask multiplexing in fractional Fourier transform domain[END_REF][START_REF] Kong | Multiple-image encryption scheme based on cascaded fractional Fourier transform[END_REF][START_REF] Tao | Double image encryption based on random phase encoding in the fractional Fourier domain[END_REF][START_REF] Liu | Double image encryption based on iterative fractional Fourier transform[END_REF][START_REF] Liu | Triple image encryption scheme in fractional Fourier transform domains[END_REF][START_REF] Wu | Four-image encryption method based on spectrum truncation, chaos and the MODFrFT[END_REF][START_REF] Joshi | Color image encryption and decryption for twin images in fractional Fourier domain[END_REF][START_REF] Joshi | Fractional Fourier transform based image multiplexing and encryption technique for four-color images using input images as keys[END_REF] as well as in gyrator transform domain [START_REF] Liu | Double image encryption by using iterative random binary encoding in gyrator domains[END_REF][START_REF] Li | Double-image encryption based on iterative gyrator transform[END_REF] were also intensively investigated. The fractional Fourier transform is widely used due to the fact that its parameter, the fractional order, enlarges the key space. In contrast, the gyrator transform has not received so far much attention and in particular for image encryption.

Throughout the literature, however, there are few studies focused on the encryption algorithm for double or multiple color images. An encryption algorithm for double color images implemented in the fractional Fourier domain was reported in [START_REF] Joshi | Color image encryption and decryption for twin images in fractional Fourier domain[END_REF], where two color images were in a first step converted into their indexed formats, a data matrix and a colormap matrix. Afterwards, the encryption is implemented on both data matrices. In most cases, such method will degrade the quality of the restituted images due to the approximation made in the transformation between the RGB color space and its indexed format. The same pre-processing strategy was used for four color images encryption scheme in [START_REF] Joshi | Fractional Fourier transform based image multiplexing and encryption technique for four-color images using input images as keys[END_REF], this encryption being carried out through cascaded fractional Fourier transforms.

In this paper, our main objective is the design of a double color image encryption algorithm based on quaternion algebra. The rest of the paper is organized as follows. Section 2 briefly reviews the main features of quaternions and quaternion gyrator transform (QGT) for describing color images. In section 3, the algorithm for encryption and decryption of double color images in quaternion gyrator domain is described. Numerical simulations are reported in section 4 to demonstrate the performance of the proposed solution. Some conclusions are drawn in section 5.

Quaternions and Quaternion gyrator transform

Quaternions

A quaternion number can be written as: q a ib jc kd = + + + , where a, b, c, d are real numbers and i, j, k are orthogonal imaginary parts obeying the following rules:

2 2 2 1 i j k = = = ! , ij ji k = ! = , jk kj i = ! = , ki
ik j = ! = . Any quaternion q can also be rewritten as:

(

) ( ) q a ib c id j = + + +

. It must be noticed that the quaternion multiplication is non-exchangeable. The conjugate of q is: q a ib jc kd = ! ! ! . When 0 a = , q becomes a pure quaternion. The quaternion q can also be expressed in polar form as:

I q q e ! =
, where the norm is

2 2 2 2 q a b c d = + + + , with
+ + = + + and 2 2 2 1 tan b c d a ! " # $ + + % & = % & '
( representing the eigenaxis and the eigenangle, respectively.

An RGB image f(x, y) can be represented by a pure quaternion matrix as

( , ) ( , ) ( , ) ( , ). R G B f x y if x y jf x y kf x y = + + (1)
where the subscripts {R, G, B} represent the color components of the image. Based on the quaternion representation, each color pixel is encoded as a pure quaternion, thus the ! correlations between them can be effectively preserved.

Quaternion gyrator transform

The gyrator transform, with rotation angle !, has been introduced by Rodrigo et al. in the field of optics [START_REF] Rodrigo | Gyrator transform: properties and applications[END_REF]. Some applications have been already reported for image encryption and specific noise reduction [START_REF] Abuturab | Color image security system using double random-structured phase encoding in gyrator transform domain[END_REF][START_REF] Sui | Color image encryption based on gyrator transform and Arnold transform[END_REF][START_REF] Liu | Double image encryption by using iterative random binary encoding in gyrator domains[END_REF][START_REF] Li | Double-image encryption based on iterative gyrator transform[END_REF][START_REF] Rodrigo | Gyrator transform: properties and applications[END_REF][START_REF] Rodrigo | Applications of gyrator transform for image processing[END_REF]. Like the fractional Fourier transform, the gyrator transform belongs to a class of linear canonical integral transform and can be viewed as an extension of the Fourier analysis. Its advantage is also to enlarge the key space and consequently strengthen the system security. The traditional gyrator transform for gray image can be extended to quaternion domain for describing color image. Due to the non-commutative property of quaternion multiplication, the right-side quaternion gyrator transform (QGT) with rotation angle ! is defined as [START_REF] Shao | Quaternion gyrator transform and its application to color image encryption[END_REF] 1

( ) cos ( ) [ ( , )]( , ) ( , ) exp 2 sin sin uv xy uy vx QG f x y u v f x y dxdy ! ! µ " ! ! + # + $ % = & ' ( ) ** (2)
where µ is a pure unit quaternion and

( ) 3
i j k µ = + + is chosen in this paper. For ! = "/2,
the right-side of Eq. ( 2) reduces to the right-side quaternion Fourier transform with the rotation of coordinates at "/2; for ! = 3"/2, to the inverse quaternion Fourier transform with the rotation of coordinates at "/2. The inverse quaternion gyrator transform corresponds to the QGT at angle ! " . The implementation of QGT is achieved through two left-side quaternion Fourier transforms (for more details about quaternion Fourier transform and quaternion gyrator transform, refer to [START_REF] Shao | Quaternion gyrator transform and its application to color image encryption[END_REF][START_REF] Said | Fast complexied quaternion Fourier transform[END_REF]).

Principle of encryption and decryption

Encryption algorithms using iterative phase retrieval technique for double gray images have been reported in [START_REF] Sui | Multiple-image encryption based on phase mask multiplexing in fractional Fourier transform domain[END_REF][START_REF] Li | Double-image encryption based on iterative gyrator transform[END_REF]. They consist of a series of forward and backward iterations. To deal with double color images, three groups of gray images (the color image is decomposed into red, green and blue channels and each component is viewed as a gray image) will be respectively encrypted. However, this approach does not make use of the inter-relationships between color components. It also leads to three times larger iteration loops and phase functions (e.g. the phase functions will go up to 12). Obviously, this can be a serious handicap for modern applications. Here, a novel encryption algorithm for double color images is proposed by combining the iterative phase retrieval algorithm with quaternion gyrator transform.

Description of the proposed algorithm

Without loss of generality, suppose that f(x, y) and g(x, y) be two color images to be encrypted, and h(x, y) the ciphertext image. The phase functions 1

! , 2 ! , 1 ! , 2 ! and ! are distributed in the interval [0, 2 ]
! , and assume that the following relationship between the images and the phases holds

2 1 2 1 11 1 12 2 21 1 22 2 exp( ) { { exp( )}exp( )} { { exp( )}exp( )} h G G f G G g ! ! " " #$ µ % µ % µ & µ & = = (3)
where µ 11, µ 12 , µ 21 , and µ 22 are the unit pure quaternions and the five phase functions are unknown at the beginning. With the expression in Eq. ( 3), the information of both color images can be hidden into the single image h. Hereinafter, an iterative phase retrieval algorithm is used to obtain the phase functions 1 ! , 2 ! , 1 ! and 2 ! . According to Eq. ( 3), the following formula can be derived

1 2 2 1 21 1 11 1 12 2 22 2 exp( ) { { { { exp( )}exp( )}}exp( )}. g G G G G f ! ! " " µ # µ $ µ $ µ # % % = % (4) ! At the initialization stage, 1 1 ! , 1 2 ! , 1 2 ! distributed in [0, 2 ]
! are randomly generated.

After n-th iteration, with the known of the phase functions

n n n n g G G G G f ! ! " " µ # µ # µ $ % % = % (5)
Then, the phase 1 n ! can be determined by

1 ârg{ } n n g g ! = " ( 6
)
where g is the conjugate of g, and arg{!} denotes the phase extraction operator.

Subsequently substituting the phase function 1 n ! into Eq. ( 4), the functions

1 2 n ! + , 1 2 n ! + , 1 1 n ! + in the next iteration are updated as belows, 2 2 1 1 1 2 11 1 12 2 21 1 1 arg { { { exp( )}exp( )}} { exp( )} n n n n G G G f G g ! " " ! # µ $ µ $ µ # % + & ' = ( ) * + (7) 2 2 1 1 1 1 2 21 1 22 2 11 1 1 arg { { { exp( )}exp( )}} { exp( )} n n n n G G G g G f ! " " ! # µ $ µ $ µ # % + + & ' = ( ) * + (8)
{ }

1 2 2 1 1 1 1 1 21 1 22 2 12 2 arg { { { { exp( )} exp( )}} exp( )} n n n n f G G G G g ! ! " " # µ $ µ $ µ # % % + + + % (9)
To evaluate the quality of the recovered color image with respect to the reference one and also to end the iteration, the correlation coefficient (CC) in color space is used

2 2 { , , } {[ ( )][ ( )]} . 3 {[ ( )] } {[ ( )] } n n c c c c n n c R G B c c c c E g E g g E g CC E g E g E g E g ! " " = " " # (10)
where n c g denotes the R, G, B channels of the recovered color image at the n-th iteration step, and E{•} the expected value operator. Once the CC value, selected as the convergence criterion, reaches a predefined threshold, then the iteration process stops.

Finally, the optimized four phases are obtained as follows

1 1 1 2 2 1 1 2 2 , , , . n n n n ! ! ! ! " " " " # = = = = (11)
According to Eq. ( 3), the final encrypted image h(x, y) is obtained by the cascaded QGT performed on color image f(x, y) with rotation angles ! 1 , ! 2 and phase functions 1

! , 2 ! as 2 1 11 1 12 2 { { exp( )}exp( )} h G G f ! ! µ " µ " = (12)
Alternatively, it can be achieved by the cascaded QGT performed on color image g(x, y) with rotation angles # 1 , # 2 and phase functions $ 1 , $ 2 as

2 1 21 1 22 2 { { exp( )}exp( )} h G G g ! ! µ " µ " = (13)
As pointed in [START_REF] Alfalou | Double random phase encryption scheme to multiplex and simultaneous encode multiple images[END_REF], a perfect encryption algorithm should parallel a decryption algorithm to allow us retrieving the secret information by means of the granted keys as a receiver. In other words, the encryption algorithm should be reversible. In the proposed encryption algorithm, the accompanied phase function % can be considered as the public key. The rotation angles ! 1 , ! 2 and the phase masks 1 ! , 2 ! are private keys for recovering the color image f(x, y), while the angles # 1 , # 2 and the phase masks $ 1 , $ 2 are private keys for recovering the color image g(x, y). The parameters µ, µ 11 , µ 12 , µ 21 , µ 22 can be viewed as auxiliary keys. With the above known essential keys, the two target color images f(x, y) and g(x, y) can be restituted by implementing the inverse encryption algorithm process as follows

1 2 12 2 11 1 { { exp( )}exp( )}exp( ). f G G h ! ! "# µ $ µ $ % % = % % (14) 1 2 22 2 21 1 { { exp( )}exp( )}exp( ). g G G h ! ! "# µ $ µ $ % % = % % (15)
!

Complexity analysis

The computational complexity of the proposed encryption algorithm will be analyzed next. It depends on two main factors. First, to extract each phase function in the n-th loop, the quaternion gyrator transform and its inverse transform must be computed four times. Here, the left-side quaternion Fourier transform is obtained via the 2D fast complex-valued Fourier transform [START_REF] Said | Fast complexied quaternion Fourier transform[END_REF] and introduced into the calculation of quaternion gyrator transform. For an image with size M " N, the computation complexity of the 2D fast Fourier transform has been shown to be

2 ( log ) MN MN
! [START_REF] Chang | Position multiplexing multiple-image encryption using cascaded phase-only masks in Fresnel transform domain[END_REF]. The second factor is the number of iterations. It has been set to 1000 in all experiments, a value considered as a good balance between time computation and accuracy in image recovery. In short, the computation load depends largely upon the image size and the iteration number.

Numerical simulation results

To corroborate the validity and the effectiveness of the proposed encryption algorithm, numerical simulations were carried out on a computer with Intel Core2 Duo CPU E8400@3.00GHz and 3G DRAM and with the MATLAB 2011a. The tested color images selected from USC-SIPI image database [41] consist of four groups as shown in Fig. 1, all with size 256"256. The system parameters are ( )

2 lm i k µ = + (l, m = 1, 2), ! 1 =0.16, ! 2 =0.25, # 1 =0.21, # 2 =0
.45 (other values can be selected). The normalized mean square error (NMSE) [START_REF] Chen | Optical color image encryption by wavelength multiplexing and lensless Fresnel transform holograms[END_REF][START_REF] Liu | Double image encryption by using iterative random binary encoding in gyrator domains[END_REF] were used to evaluate the similarity between the decrypted image ˆ( , ) f x y and its reference f(x, y), which is defined by

1 1 2 0 0 1 1 2 0 0 ( , ) ( , ) NMSE . ( ,
) N M x y N M x y f x y f x y f x y ! ! = = ! ! = = ! = " " " " (16)
Simulations have been done first to test the security and the reversibility of the algorithm. The four phase functions for each tested group were obtained after 1000 iteration loops. The evolution of CC values is displayed in Fig. 2, where the images in the upper row and the lower row for each tested group shown in Fig. 1 are respectively used as the reference image to terminate the iteration in the phase retrieval process. Table 1 lists the NMSE values of decrypted images, where all the images in the lower row are used as reference image to obtain the four phases. It should be noted that for any tested group, four sets of decrypted results can be obtained due to the fact that each of them can be used as the reference to determine the convergence or for encryption. From the results of table 1, it can be observed that for each tested group, the error of the first decrypted image is larger due to the fact that the encypted data h obtained from the second tested image g(x, y) is used as the public keys. When using Eq. ( 12), the error of the second decrypted image will be larger which is caused by the smaller residual error when obtaining phase functions by means of the phase retrieval algorithm. Figure 3 shows the results of encryption and granted decryption of the first tested group (Lena and Pepper). It can be seen that with the knowledge of the granted keys, the original color images can be restituted with a high quality. In contrast, without knowing the rotation angles or phase functions, almost no useful information can be recovered from the decrypted results, as demonstrated in Fig. 4. In order now to test the sensitivity of rotation angles to small changes, the decryption was performed by fixing one angle value and by varying the other, where the deviation # ranges from -0.10 to 0.10 with step 0.005. Others keys used for encryption and decryption in this simulation are set to the correct values. From the NMSE values shown in Fig. 5, it can be concluded that when the deviation is equal or greater than 0.01, the decrypted images (for brevity not displayed here) cannot be recognized clearly. Therefore, the rotation angles of QGT are sensitive keys. During data transmission, information loss or noise contamination may occur. Another simulation was carried out so as to evaluate the performance of the proposed encryption algorithm against loss of data. To demonstrate the interest of the present encryption algorithm for double color images, a comparison with the encryption algorithm based on fractional Fourier transforms (FrFT) [START_REF] Joshi | Color image encryption and decryption for twin images in fractional Fourier domain[END_REF] is provided. To do so, the information loss severity p (varying from 10% to 60% with an increment 10%) around the center is regarded as an occlusion attack to the encrypted data. From the results illustrated in Fig. 6, it can be seen that the NMSE values of the FrFT-based encryption algorithm are more stable. However, the NMSE values of the present algorithm evolve linearly with the percentage p. When increasing the occlusion percentage p, the similarity between recovered image and its original one becomes obviously smaller. This behavior can be explained by the fact that a significant part of the QGT information is distributed around the image center considered as the origin. Although the NMSE values of the FrFT-based encryption algorithm are smaller in some limited intervals, we found that no useful information can be recovered from the decrypted image even when the smallest occlusion percentage 10% is used in reduplicative experiments. Conversely, as displayed in Fig. 7, relevant decrypted information can be obtained when using the proposed algorithm for p values up to 50%. In other words, the proposed encryption ! algorithm shows more robustness against loss of encrypted data. In order to test the robustness of our approach against noise attack, a white Gaussian noise with zero mean and varying standard deviations & was added to the four tested group images. The same noise was added to the encrypted images obtained with fractional Fourier transforms [START_REF] Joshi | Color image encryption and decryption for twin images in fractional Fourier domain[END_REF] for comparison purpose. Fig. 8 displays the NMSE values of the decrypted images from these noisy encrypted images with the FrFT-based and the QGT-based encryption algorithm. The NMSE values corresponding to the proposed algorithm are much smaller. Two examples of the resulting decrypted images are shown in Fig. 9, where the noise standard deviation is 0.05. It can be concluded that the proposed encryption algorithm shows here also more robustness against noise.

Conclusion

A novel double color image encryption algorithm was presented in this paper. The originality stands on the holistic representation for color image by using quaternion algebra and iterative phase retrieval algorithm in quaternion gyrator domain. The phase masks used for encryption are obtained through iterative phase retrieval algorithm in quaternion gyrator domain, and the encrypted image is obtained via cascaded quaternion gyrator transforms. Knowing all the granted keys, as the parameters of quaternion gyrator and phases, two color images can be restituted with a high quality. The results of the numerical simulations have demonstrated the feasibility and the effectiveness of the proposed encryption system. When compared with the double color images encryption algorithm in the fractional Fourier domain, the proposed encryption algorithm has shown a better robustness to occlusions and additive Gaussian noise. It must be also emphasized that, when both color images need to be encrypted for the FrFT-based encryption algorithm, our encryption approach only requires to encrypt one image, and thus provides savings in data storage while improving the transmission efficiency. To further reinforce the encryption system, other techniques such as sparse representation will be taken into consideration in the near future. In addition, the proposed encryption algorithm can also be extended to encrypt multiple gray or colored images.

Fig. 1 .

 1 Fig. 1. Tested color images: each column constitutes a tested group (Primary color images are selected from USC-SIPI image database).

Table 1 . 30 Fig. 2 .Fig. 3 .

 13023 Fig. 1 (b)

Fig. 4 .

 4 Fig. 4. Decrypted results with randomly guessed keys: (a) and (b) correspond to Lena, (c) and (d) correspond to Pepper: (a) and (c) correspond to phases randomly guessed, (b) and (d) correspond to angles randomly guessed.

Fig. 5 .

 5 Fig. 5. The NMSE values of decrypted images when varying one rotation angle while the other

Fig. 6 .

 6 Fig. 6. The NMSE values of decrypted images with different information loss severities.

Fig. 7 .

 7 Fig. 7. Decrypted results of encrypted data with 50% information loss.

Fig. 8 .

 8 Fig. 8. The NMSE values of decrypted images with different Gaussian white noise using different encryption algorithms.

Fig. 9 .

 9 Fig. 9. Examples of decrypted images from noisy encrypted data with standard deviation 0.05.

& & & &

Acknowledgments

This work was supported by the National Basic Research Program of China under Grant 2011CB707904, by the National Natural Science Foundation of China under Grants 61073138, 61271312, 61201344, and 11301074, by the Ministry of Education of China under Grants 20110092110023 and 20120092120036, the Key Laboratory of Computer Network and Information Integration (Southeast University), Ministry of Education, and by Natural Science Foundation of Jiangsu Province under Grants SBK200910055 and BK2012329. The authors would like to thank the anonymous reviewers for their constructive comments and suggestions.