Arthur Sarthou

Stéphane Vincent

Jean-Paul Caltagirone

Consistant velocity-pressure coupling for second-order L 2 -penalty and direct-forcing methods

Keywords: Navier-Stokes equations, Fictitious domain, Velocity-pressure coupling, Augmented Lagrangian, Pressure-correction methods, Projection methods, Fractionnal-step methods, Penalty method, Penalization method, Immersed boundary method, Incompressible flows

The present work studies the interactions between fictitious-domain methods on structured grids and velocitypressure coupling for the resolution of the Navier-Stokes equations. The pressure-correction approaches are mainly used in this context but the corrector step is generally not modified consistently to take into account the fictitious domain. A consistent modification of the pressure-projection for a high-order penalty (or penalization) method close to the Hikeno-Kajishima modification for the Immersed Boundary Method is presented here. Compared to the first-order correction required for the L 2 -penalty methods, the small values of the penalty parameters do not lead to numerical instabilities in solving the Poisson equation. A comparison of the corrected rotational pressure-correction method with the augmented Lagrangian approach which does not require a correction is carried out.

Introduction

The simulation of real heat and mass transfers often implies interactions between multiphase flows and complex obstacles. Many simulation codes based on structured grids have shown their ability to deal with a large amount of physical phenomena. However, structured grids can not generally match complex interfaces due to their lack of flexibility, and the treatment of problems with complex shapes is unnatural and uneasy with this approach. The fictitious-domain methods have been designed to improve the performances of structured grid codes when complex shapes are necessary. A wide literature has been devoted to the subject during the last decades, especially the last 15 years with the emergence of high-order methods (see for a review [START_REF] Mittal | Immersed boundary methods[END_REF][START_REF] Sarthou | Fictitious domain methods for the elliptic and Navier-Stokes equations[END_REF]).

A first approach to deal with immersed boundaries is the Distributed Lagrangian Multiplier method proposed by Glowinski [START_REF] Glowinski | A distributed lagrange multiplier/fictitious domain method for particulate flows[END_REF]. Used for the Navier-Stokes equations with finite-element methods, the coupling between fluid and solid media is ensured with Lagrange multipliers introduced into the weak formulation of the Navier-Stokes equations. Cartesian grid methods [START_REF] Johansen | A cartesian grid embedded boundary method for poisson's equation on irregular domains[END_REF][START_REF] Mccorquodale | A cartesian grid embedded boundary method for the heat equation on irregular domains[END_REF] propose to use structured grid in the whole domain except near obstacles where unstructured cells are created from structured cells. The implementation of the method is not simple due to the numerous different space configurations of the intersections between cells and objects. The existence of small cells can induce solver trouble and need a special treatment. The Immersed Boundary Method (IBM) was initially presented by Peskin [START_REF] Peskin | Flow patterns around heart valves: A numerical method[END_REF][START_REF] Peskin | The immersed boundary method[END_REF]. Fictitious boundaries are taken into account through a source term activated only near the boundaries. As the source term is weighted with a discrete Dirac function with a non-zero support, the interface influence is spread over some grid cells and a first order of spatial convergence is generally obtained. Another class of IBM, the Direct-forcing (DF) methods, was initially proposed by Mohd-Yusof [START_REF] Mohd-Yusof | Combined immersed boundary/b-spline methods for simulations of flows in complex geometries[END_REF]. The idea here is to impose a no-slip condition directly on the boundary using a mirrored flow over the boundary. In [START_REF] Fadlun | Combined immersed-boundary finitedifference methods for three-dimensional complex flow simulations[END_REF][START_REF] Verzicco | Flow in a impeller stirred tank using an immersed boundary method[END_REF], the correct boundary velocity is obtained by interpolating the solution on the boundary and far from the boundary on grid points in the near vicinity of the interface. In [START_REF] Tseng | A ghost-cell immersed boundary method for flow in complex geometry[END_REF], Tseng and Ferziger use the same principle but extrapolate the solution in ghost cells inside the boundary. This approach can be seen as a generalization of the mirror boundary condition used with the Cartesian staggered grid to impose Dirichlet condition on pressure nodes. As the original conservation equations are solved one node closer to the boundary this approach seems to be more accurate than [START_REF] Fadlun | Combined immersed-boundary finitedifference methods for three-dimensional complex flow simulations[END_REF][START_REF] Verzicco | Flow in a impeller stirred tank using an immersed boundary method[END_REF].

Originally presented in [START_REF] Arquis | Convection mixte dans une couche poreuse verticale non confinée[END_REF] for the conservation equations, the penalty (or from the french designation, penalization) methods for fictitious domains consist in adding specific terms in the conservation equations to play with the order of magnitude of existing physical contributions so as to obtain at the same time and with the same set of equations various physical properties. The Volume Penalty Method (VPM, [START_REF] Khadra | Fictitious domain approach for numerical modelling of navier-stokes equations[END_REF] and the references therein) is a simple way to impose a solution in a part of the numerical domain. The methods imposing the solution are called the L 2 -penalty methods while the H 1 -penalty methods allows a derivative of the solution to be imposed [START_REF] Angot | A penalization method to take into account obstacles in incompressible viscous flows[END_REF]. Classical penalty methods are of first order only since they consider the projected shape of the interface on the Eulerian grid to define the penalty parameters [START_REF] Ramière | A general fictitious domain method with immersed jumps and multilevel nested structured meshes[END_REF]. In [START_REF] Sarthou | Eulerian-Lagrangian grid coupling and penalty methods for the simulation of multiphase flows interacting with complex objects[END_REF][START_REF] Sarthou | Finite Volumes for Complex Applications V, chapter The sub-mesh-penalty method[END_REF], Sarthou et al. have extended penalty methods to higher orders by modifying the expression of the penalty term using implicit interpolations as in [START_REF] Tseng | A ghost-cell immersed boundary method for flow in complex geometry[END_REF]. The method is called the Sub-mesh penalty method (SMPM) and has been applied first to elliptic equations. For the Navier-Stokes equations, the incompressibility of the flows has been ensured with an augmented Lagrangian velocity-pressure coupling [START_REF] Fortin | Méthodes de lagrangien augmenté[END_REF][START_REF] Vincent | Augmented lagrangian and penalty methods for the simulation of two-phase flows interacting with moving solids. application to hydroplaning flows interacting with real tire tread patterns[END_REF]. This method consists in solving the momentum equation with an augmented Lagrangian term which enforces the divergence-free constraint during an iterative process. Hence, the fulfilment of both incompressibility and boundary constraints is obtained with the resolution of a unique equation.

An issue occurs with the pressure-correction methods as the resolution of an additional elliptic equation is performed to rise a pressure and to obtain a solenoidal field. The IBM for the Navier-Stokes equations are generally used with the pressure-correction methods and modify the predictor step only (where the momentum equation is solved). As no modification of the corrector step is performed, the additional boundary constraint is not taken into account in the final velocity and pressure fields and is then no more respected.

In [START_REF] Domenichini | On the consistency of the direct forcing method in the fractional step solution of the navier-stokes equations[END_REF], Domenichini analyzes in details the application of the DF-IBM to the fractional step solution of the Navier-Stokes equations. To focus on the error induced by the non-consistant application of the immersed boundary condition, a spectral solver is used. As can be expected, he notices that the boundary condition is not accurately imposed, even if sub-iterations of the time-splitting can be performed to reduce this error. This problem is not frequently tackled in the literature and fully satisfactory solutions have been found only recently. In [START_REF] Kim | An immersed boundary finite volume method for simulations of flow in complex geometries[END_REF], authors use a mass source and sink term in the pressure equation to preserve the mass balance in the boundary cells but the desired velocity is not exactly imposed on the wall. More recently, Taira and Colonius [START_REF] Taira | The immersed boundary method: A projection approach[END_REF] consider both the boundary forcing of the Peskin IBM and the pressure as Lagrange multipliers. Hence, the time-splitting procedure is applied in the same time and in an equal manner to both quantities. It allows the rigid body and the incompressibility constraints to be satisfied at the same time.

In [START_REF] Ikeno | Finite-difference immersed boundary method consistent with wall conditions for incompressible turbulent flow simulations[END_REF], Ikeno and Kajishima propose a consistant correction for a second-order DF-IBM. The principle is to add the boundary term in the projection step in a consistant way. The update equation of the velocity has to be modified too.

Concerning the L 2 -penalty methods, a solution for the first-order method is recalled in [START_REF] Angot | Finite Volumes for Complex Applications V, chapter Vector Penalty-Projection Methods for the Solution of Unsteady Incompressible Flows[END_REF]. However, applying this modification to a high-order method is quite more challenging. Recently, the correction of [START_REF] Ikeno | Finite-difference immersed boundary method consistent with wall conditions for incompressible turbulent flow simulations[END_REF] has been applied successfully to the SMPM [START_REF] Sarthou | Fictitious domain methods for the elliptic and Navier-Stokes equations[END_REF]. A correction for a direct-forcing penalty method has been introduced by Belliard et al. [START_REF] Belliard | Penalized direct forcing and projection schemes for navierstokes[END_REF].

A consistant correction for a fully implicit high-order L 2 -penalty method is proposed here. The formulation is derived from the penalized momentum equation and thus naturally obtained. Compared to [START_REF] Belliard | Penalized direct forcing and projection schemes for navierstokes[END_REF], all steps of the method take into account the high-order of the penalty term.

The method is applied to the incremental Goda [START_REF] Goda | A multistep technique with implicit difference schemes for calculating two-or threedimensional cavity flows[END_REF] and rotational [START_REF] Timmermans | An approximate projection scheme for incompressible flow using spectral elements[END_REF] pressure-projection methods coupled with the SMPM. These approaches are compared in time and space to the augmented Lagrangian method.

In Section 2, the conservation equations and their discretization are presented. Then, the SMPM is described. The third Section focuses on the consistant correction for the time-splitting methods. In Section 4, numerical tests are performed to study and compare the numerical convergence of the method. The last section concludes the article with a discussion and perspectives are drawn.

Governing equations and base discretization 2.1 Governing equations

We consider the following form of the incompressible Navier-Stokes equations in a domain Ω:

∇ • u = 0 (1) ρ ( ∂u ∂t + (u • ∇)u ) = -∇p + ∇ • [µ(∇u + ∇ T u)] (2)
with u the velocity, ρ the variable fluid density, p the pressure, g the gravity vector, and µ the dynamic viscosity. The Navier-Stokes equations are discretized with implicit finite-volumes on a staggered Cartesian grid. A second-order centered scheme is used to approximate the spatial derivatives while first-order Euler and secondorder Gear schemes are used for the time integration. All the terms are written at time (n + 1)∆t, ∆t being the time-step, except for the non-linear term u n+1 • ∇u n+1 which is linearized as u n • ∇u n+1 for the first-order Gear scheme and as (2u nu n-1 ) • ∇u n+1 for the second-order Gear scheme. The modified semi-discrete form of the original equation ( 2) is then

ρ ( γ 1 u n+1 + γ 2 u n + γ 3 u n-1 ∆t + ( (γ 4 u n + γ 5 u n-1 ) • ∇ ) u n+1 ) = -∇(γ 4 p n + γ 5 p n-1 ) + ∇ • [µ(∇u n+1 + ∇ T u n+1 )] (3)
with the additional constraint ∇ • u n+1 = 0. The values of γ i depend on the temporal scheme as

• γ 1 = 1, γ 2 = -1, γ 3 = 0, γ 4 = 1, γ 5 = 0 for the Gear 1 or Euler scheme • γ 1 = 3 2 , γ 2 = -2, γ 3 = 1 2 , γ 4 = 2, γ 5 = -1
for the Gear 2 scheme. In the next parts, the Euler scheme is generally written for the sake of simplicity.

The linear system resulting from the discretization is solved with a BiCG-Stab II solver [START_REF] Gustafsson | On First and Second Order Symmetric Factorization Methods for the Solution of Elliptic Difference Equations[END_REF], preconditioned by a Modified and Incomplete LU method [START_REF] Van Der | Bi-cgstab: a fast and smoothly converging variant of bi-cg for the solution of nonsymmetric linear systems[END_REF].

The velocity-pressure coupling 2.2.1 Pressure-correction methods

Most of the finite-volume CFD codes on Eulerian grids use pressure-correction (or fractional-step) methods. The idea is to obtain first a predicted velocity u * satisfying the momentum equation only. This field is not solenoidal as nothing constrains this condition. In a second step, the projection, the pressure is risen with respect to the divergence of u * . The third step consists in updating the velocity according to the pressure gradient obtained with the second step. We consider RHS the sum of the convective and diffusive terms of the equation (2). The half discretization in time gives:

ρ ( u n+1 -u n ∆t ) = RHS n+1 -∇p n+1 . ( 4
)
This equation is solved, but as here ∇•u n+1 ̸ = 0, the solution is denoted u * . We define u ′ such as u n+1 = u ′ +u * and p ′ such as p n+1 = p ′ + p * . Hence, the predictor step solves:

ρ ( u * -u n ∆t ) = RHS * -∇p n . ( 5
)
To obtain the final velocity, the following equation is used :

ρ ( u n+1 -u * ∆t ) = -∇p ′ . ( 6
)
This equation can be constructed with two points of view. For the first one, we consider (4)-(5) while RHS ′ is neglected (implementations which keep RHS ′ are more difficult to perform) introducing an additional error as the convective and diffusive terms are only considered at step * . The second point of view uses the Hodge-Helmholtz orthogonal decomposition of the space L 2 (Ω)

d = H⊕H ⊥ where H = { u ∈ L 2 (Ω) d , ∇ • u = 0, u • n = 0 on ∂Ω } and H ⊥ = { ∇ϕ, ϕ ∈ H 1 (Ω) } .
Hence, the predicted velocity field can be corrected by a pressure gradient to obtain a solenoidal field. Equation ( 6) is considered as the second step of a time-splitting where the part of the solution deriving from a potential is added to the predicted field to obtain the solenoidal solution. The pressure increment is practically obtained by solving the divergence of (6) :

∇ • u * = ∇ • ∆t ρ ∇p ′ . ( 7
)
Once the pressure increment is obtained, velocity and pressure are updated:

p n+1 = p ′ + p n (8) u n+1 = u * - ∆t ρ ∇p ′ . ( 9
)
In [START_REF] Timmermans | An approximate projection scheme for incompressible flow using spectral elements[END_REF], Timmermans et al. proposes a correction of this last step replacing (8) by

p n+1 = p ′ + p n -µ∇ • u * . ( 10
)
This correction gives a consistent pressure boundary condition while the standard incremental algorithms gives an artificial Neumann boundary condition for the pressure. An overview of the different projection methods is performed in [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF]. Concerning the fictitious domains, the IBM applied to the NS equations are generally designed for the projector step only. As no modification of the corrector step is performed, the additional boundary constraint is not taken into account and is then not respected.

The augmented Lagrangian method

The augmented Lagrangian (AL) method [START_REF] Fortin | Méthodes de lagrangien augmenté[END_REF] consists in adding a term ∇(dr∇ • u) to the momentum equation of the NS equations so as to enforce the divergence free constraint. The pressure is updated with the Uzawa method [START_REF] Arrow | Studies in linear and Nonlinear Programming -Iterative method for concave programming[END_REF]. The parameter dr sets the magnitude of the constraint and must be chosen according to the magnitude of the other terms of the equation to avoid low numerical performances of the solver and to obtain a suitable physical solution. Iterative solvers can be very sensitive to the magnitude of dr (the condition number of the matrix varies linearly with respect to dr [START_REF] Févriére | On the penalty-projection method for the navierstokes equations with the mac mesh[END_REF] and the direct solvers allow higher values of dr to be taken) and a high parameter implies an increase of the number of internal iterations of the solver. A too high parameter penalizes the initial equation and leads to a strictly incompressible velocity field with no respect to the initial momentum equation. Choosing a suitable parameter is not trivial. Furthermore, multiphase flows [VSC] can induce strong variations of the densities and viscosities, and require dr to vary accordingly. To tackle this issue, the parameter dr can be determined according to the physical quantities [START_REF] Vincent | An adaptative augmented lagrangian method for three-dimensional multimaterial flows[END_REF] or the coefficients of the discretization matrix [VSC].

The base algorithm of the augmented Lagrangian method is now described. Starting with u * ,0 = u n and p * ,0 = p n , while

||∇ • u * ,m || > ϵ , solve (u * ,0 , p * ,0 ) = (u n , p n ) ρ ( u * ,m -u * ,0 ∆t + u * ,m-1 • ∇u * ,m ) -∇(dr∇ • u * ,m ) = -∇p * ,m-1 + ∇ • [µ(∇u * ,m + ∇ T u * ,m )] p * ,m = p * ,m-1 -dr∇ • u * ,m (11)
Although the AL method is an iterative procedure, one iteration is generally acceptable to reach a sufficiently small divergence. To enforce an immersed Dirichlet BC, the penalty term χ ε (Πu n+1 -u D ) (described in details in section 3) is added to the momentum equation. We obtain the following simplified formulation :

ρ ( u n+1 -u n ∆t + u n • ∇u n+1 ) -∇(dr∇ • u n+1 ) = -∇p n + ∇ • [µ(∇u n+1 + ∇ T u n+1 )] + χ ε (P u n+1 -u D ) (12)
p n+1 = p n -dr∇ • u n+1 . ( 13
)
This last equation is not a splitting in time. By taking the divergence of (13), on can see that ∇ • p n+1 is still present in (12) and dr∇ • u n+1 can be seen as the implicit pressure increment. Hence, the AL methods allows large time steps to be used. Furthermore, no boundary conditions are required for the pressure. A small number of AL iterations could be required to obtain an acceptable divergence [VSC]. If a machine accuracy fulfilment of the divergence-free constraint is desired, the number of required iterations can be prohibitive (especially with iterative solvers). A simple solution is to use a penalty-projection method as presented in the next section.

The penalty-projection methods

The penalty-projection methods [START_REF] Shen | On error estimates of some higher order projection and penalty-projection methods for navier-stokes equations[END_REF][START_REF] Févriére | On the penalty-projection method for the navierstokes equations with the mac mesh[END_REF] are a combination of the augmented Lagrangian and pressurecorrection methods. An augmented Lagrangian term (called in this case the penalty term) is added to the momentum equation which produces a velocity field u * with a moderate divergence. A projection step is then performed to obtain a solenoidal velocity field u n+1 . The final pressure is obtained with p n+1 = p n + p ′ -dr∇ • u * for the standard incremental scheme. For the rotational scheme, dr is replaced by dr + µ in the pressure update. Compared to the augmented Lagrangian method taken alone, the weaker requirement on the divergence for the first step allows to use a smaller dr parameter (which decreases the conditioning number of the linear system) and less AL sub-iterations. Compared to the pressure-correction methods, the intermediate predicted velocity field is closer to the solenoidal final field. Another penalty-projection has been proposed by Breil and Caltagirone in [START_REF] Caltagirone | A vectorial projection method for solving the navier-stokes equations[END_REF]. The difference lies in the projection step which solves a vector equation on the velocity:

∇(∇ • u ′ ) = -∇(∇ • u * ) (14)
where u ′ is a velocity increment such as u n+1 = u * + u ′ . The authors demonstrates that a very small amount of solver iterations (less than ten for a BiCG-Stab II with ILU preconditioning) are required to reach the divergencefree constraint. An extension of this method is presented in [START_REF] Angot | Finite Volumes for Complex Applications V, chapter Vector Penalty-Projection Methods for the Solution of Unsteady Incompressible Flows[END_REF].

3 The L 2 -penalty methods

The L 2 -penalty methods are a class of fictitious domain method used to impose a Dirichlet or Neumann boundary condition on a complex interface. We specify to avoid confusion that the penalty methods and the penaltyprojection methods are not related except that both add a term in the conservation equation to enforce a specific behavior of the solution.

Base principle

Let us consider the original domain of interest denoted by Ω 0 , typically the fluid domain, which is embedded inside a simple computational domain Ω ⊂ R d . The auxiliary domain Ω 1 , typically a solid particle or an obstacle, is then such that Ω = Ω 0 ∪ Σ ∪ Ω 1 where Σ is an immersed interface (see Fig. 1 left). Let n be the unit outward normal vector to Ω 0 on Σ. Our objective is to numerically impose the adequate boundary or interface conditions on the interface Σ. The continuous L 2 -penalty method for the incompressible Navier-Stokes equations consists in adding a term χ ε (u -u D ) into the momentum equation

ρ ( ∂u ∂t + (u • ∇)u ) = -∇p + ρg + ∇ • [µ(∇u + ∇ T u)] + χ ε (u -u D ) (15)
where 0 < ε ≪ 1 denotes the penalty parameter and χ is the Heaviside function such as

χ(x) = { 1 if x ∈ Ω 1 0 if x ∈ Ω 0 .
In Ω 0 , the penalty term vanishes and the original momentum equation is retrieved. In Ω 1 , the equation tends to u = u D when ε -→ 0.

Discretization

For the sake of simplicity the method is described in 2D for a scalar equation. The computational domain Ω is approximated with a curvilinear mesh T h composed of N × M (×L in 3D) cell-centered finite volumes (V I ) for I ∈ E, E being the set of index of the Eulerian orthogonal curvilinear structured mesh. Let x I be the vector coordinates of the center of each volume V I . The local characteristic space step h I of the volume V I is defined as the maximum length of V I in each direction, whereas h denotes the Eulerian mesh step: h = sup I∈E h I . This grid is used to discretize the conservation equations. A dual grid is introduced for the management of the penalty method (in finite-difference discretization, the primal mesh is used). The grid lines of this dual cell-vertex mesh are defined by the network of the cell centers x I . The volumes of the dual mesh are denoted by (K I ). The Eulerian unknowns are noted ϕ I which are the approximated values of ϕ(x I ), i.e. the solution at the cell centers x I .

The discrete interface Σ h , hereafter called the Lagrangian mesh, is given by a discretization of the original interface Σ. It is described by a piecewise linear approximation of Σ:

Σ h = {σ l ∈ P d-1 1 , l ∈ L f }, K being the cardinal of L f
and L f being the set of index of the Lagrangian mesh. Typically, σ l are segments in 2D and triangles in 3D. The vertices of each face σ l are denoted by x l,i for i = 1, d and the set of all vertices is : {x l , l ∈ L v }. The intersection points between the grid lines of the Eulerian dual mesh and the faces σ l of the Lagrangian mesh are denoted by {x i , i ∈ I} (see Fig. 1 right). The cell centers x I are sorted according to their location inside Ω 0 or Ω 1 with the discrete Heaviside function χ I = χ(x I ). This function is computed from Σ h with a Thread Ray-casting method [START_REF] Sarthou | A second-order curvilinear to cartesian transformation of immersed interfaces and boundaries. application to fictitious domains and multiphase flows[END_REF]. The principle is to cast a ray from each Eulerian nodes. If the number of intersections between Σ h and the ray is odd, the node is inside the object, otherwise outside. The algorithm needs LM N K/ max(L, M, N ) intersection tests and is faster than the classical Ray-casting method [START_REF] Ogayar | Point in solid strategies[END_REF] which requires LM N K intersection tests. New sets of Eulerian points x I are defined near the interface such as one neighbor x J verifying χ J ̸ = χ I exists, i. e. the segment [x I ; x J ] is cut by Σ h . These Eulerian "interface" points are also sorted according to their location inside Ω 0 or Ω 1 . Two sets {x I , I ∈ N 0 } and {x I , I ∈ N 1 } are so obtained, where

N 0 = {I, x I ∈ Ω 0 , ∃ neighb(x I ) ∈ Ω 1 } and N 1 = {I, x I ∈ Ω 1 , ∃ neighb(x I ) ∈ Ω 0 }.
The spatial order of the L 2 -penalty method directly depends on the discretization of χ ε (ϕ -ϕ D ). The VPM [START_REF] Angot | Contribution à l'étude des transferts thermiques dans des systèmes complexes aux composants électroniques[END_REF] discretizes the term for a node ϕ i as χi ε (ϕ i -ϕ D ) with ϕ D the boundary value at the vicinity of x i . As the solution is considered as constant in each penalized cells, the shape of the immersed boundary is perceived by the fluid as being stair-stepped, or rasterized. This inaccurate description of the interface implies a first-order of spatial convergence [START_REF] Ramière | Convergence analysis of the q 1 -finite element method for elliptic problems with non-boundary-fitted meshes[END_REF].

To reach a second-order of spatial accuracy, the sub-mesh penalty method (SMPM) [SVCA08, SVAC08] discretizes the penalty term with χi ε (Π i ϕ -ϕ D (x l )) where Π i is a polynomial interpolator. We consider a point x I with I ∈ N 1 and only one neighbor x J in Ω 0 . The Lagrangian point x l is the intersection between [x I ; x J ] and Σ h (Fig. 1 right). The solution ϕ = ϕ l has to be satisfied at x l which implies in a discrete point of view Π I ϕ(x l ) = ϕ D (x l ) with Π the P 1 1 -interpolator (one-dimensional linear polynomial) between the Eulerian unknowns u I and u J :

ϕ l = α I ϕ I + α J ϕ J with 0 < α I , α J < 1 and α I + α J = 1 (16)
The coefficients α I and α J are determined such as

Π I ϕ(x I ) = ϕ I and Π I ϕ(x J ) = ϕ J . If now x I has a second neighbor x K in Ω 0 , the intersection x m between [x I ;
x K ] and Σ h is considered. We choose x p , a new point of Σ h between x l and x m (see Fig. 1 left). Practically, the barycenter between x m and x l is used. The resulting point x p is not necessarily on Σ h but it does not spoil the second-order precision of the method since the distance d(x p , Σ h ) between x p and Σ h is varying like O(h 2 ). The solution ϕ p (x p ) is then approximated using a P 2 1 -interpolation (two-dimensional linear polynomial)of the values ϕ I , ϕ J and ϕ K :

ϕ p = α I ϕ I + α J ϕ J + α K ϕ K , 0 < α I , α J , α K < 1 , α I + α J + α K = 1 (17)
We can also use a Q 2 1 -interpolation of ϕ I , ϕ J , ϕ K and ϕ L , by extending the interpolation stencil with the point x L which is the fourth point of the cell of the dual mesh defined by x I , x J and x K (see Fig. 1

.right).

If Σ h is regular enough, x I has almost never a third neighbor in Ω 0 . However, if it is the case, the first-order L 2 -penalty term χI ε (ϕ I -ϕ D (x l )) is used. In any case, by decreasing the Eulerian mesh step h, we also decrease the number of points x I having more than two neighbors in Ω 0 . It has to be noticed that the penalty term is only required for nodes x I with I ∈ N 1 , i.e. for the nodes in Ω 1 with at least one neighbor in Ω 0 . What happens inside Ω 1 and far from Σ has no impact on the flow in Ω 0 and so is of secondary interest (practically, the parts of the discretization matrix corresponding to these nodes are removed before the resolution of the liner system).

Concerning the application of the method to a Curvilinear grid, the computational domain can be "unfold" into a Cartesian domain where most of the computations are performed. Hence, the penalty constraints are build from the Cartesian grid with the standard Cartesian routines. The curvilinear to Cartesian projection is described in [START_REF] Sarthou | A second-order curvilinear to cartesian transformation of immersed interfaces and boundaries. application to fictitious domains and multiphase flows[END_REF] .

4 Penalty correction for the pressure correction methods

First-order correction

For a first-order penalty term, such a correction is easy to performed [START_REF] Angot | Finite Volumes for Complex Applications V, chapter Vector Penalty-Projection Methods for the Solution of Unsteady Incompressible Flows[END_REF]. Let us now write the momentum Navier-Stokes equation with a first-order penalty term :

ρ ∂u ∂t = +RHS -∇p + χ ε (u -u D ) (18)
In the standard method, (4)-( 5) is considered to obtain the equation ( 6). The same operation is performed with the penalty term and the equation (6) becomes:

ρ ( u n+1 -u * ∆t ) = -∇p ′ + χ ε (u n+1 -u * ) (19)
and the equation (7) becomes:

∇ • u * = ∇ • ( ∆t ρ - χ ε ) -1 ∇p ′ . ( 20
)
The velocity is then updated using (19) :

u n+1 = u * -( ∆t ρ - χ ε ) -1 ∇p ′ . ( 21
)
It is important to notice that the correction terms due to the fictitious domain method appears naturally with this walkthrough. The standard pressure update is not modified.

Remark 4.1 A classical problem with the elliptic equation (20) in the penalized approach is that the diffusion coefficient varies in O(ε -1 ) in Ω 1 . Hence, a strong imposition of the penalty constraint produces a very high diffusion coefficient which leads to numerical instabilities. In [START_REF] Belliard | Penalized direct forcing and projection schemes for navierstokes[END_REF], with the assumption that ρ is constant over the domain, the authors propose the following correction to avoid numerical instabilities :

ρ ∆t ∇ • u * = ∇ • ε ε + χ ∇p n+1 + 1 ε (p n+1 -p 0 ) ( 22
)
where p 0 is a prescribed pressure. The last term produces a L 2 -penalization of the pressure equation.

Remark 4.2 In [START_REF] Belliard | Penalized direct forcing and projection schemes for navierstokes[END_REF], authors use the non-incremental version of the scheme and a first-order correction for the pressure equation. However, the boundary condition is explicitly chosen as being the linear extrapolation of the solution and a second-order of spatial convergence is reached for a stationary case. In our approach, all the terms due to the penalization are of second-order and implicit.

Higher-order correction

A fully implicit second-order correction is now proposed. A consistant correction following the precedent walkthrough with a penalty term of higher order is much more delicate. The first-order term χi ε

(u i -u D (x l )) is replaced by χi ε (Π i u -u D (x l )) with Π i u = ∑ j/xj ∈N eighb(xi)
α j u j . The pressure equation (20) becomes, for each node x i :

∇ • u * i = ∇ • ( ρ ∆t - χ i ε Π i ) -1 ∇p ′ (23)
which requires the calculation of the matrix corresponding to ( ρ ∆t -χi ε Π i ) -1 before the resolution of the initial system. A more simple method is chosen here even if it would be interesting to evaluate this first method. We start back with the equation (6) with a penalty term by following the natural walkthrough of the first-order correction :

u n+1 i -u * i = - ∆t ρ ∇p ′ + ∆tχ i ρε Π i ( u n+1 -u * ) (24)
developed as

(1 - ∆tχ i ρε Π i )(u * -u n+1 ) = ∆t ρ ∇p ′ . ( 25
)
The key of this construction is to introduce Π ′ i , a new interpolator such as:

Π ′ i u = ∑ j/xj ∈N eighb(xi),i̸ =j α j u j . ( 26
)
Hence,

Π i u = α i u i + Π ′ i u and we have ( 1 - ∆tχ i ρε (α i + Π ′ i ) ) (u * -u n+1 ) = ∆t ρ ∇p ′ . ( 27
)
By construction, the nodes involved by Π ′ i are always in Ω 0 where χ = 0 so the original pressure-correction method occurs. Hence, using (9) we have

( 1 - ∆tχ i α i ρε ) (u * i -u n+1 i ) - ( ∆tχ i ρε Π ′ i ∆t ρ ∇p ′ ) = ∆t ρ ∇p ′ (28)
and we obtain

u n+1 i = u * i - ( ∆tχ i ρε -∆tχ i α i Π ′ i ∆t ρ ∇p ′ + ∆tε i ρε -∆tχ i α i ∇p ′ ) . ( 29
)
the equation of the velocity update. Using the divergence of (29) and ∇ • u n+1 i = 0 we obtain the final correction equation

∇ • u * i = ∇ • ( ∆tχ i ρε -∆tχ i α i Π ′ i ∆t ρ ∇p ′ + ∆tε i ρε -∆tχ i α i ∇p ′ ) . ( 30
)
The parameter ε ≪ 1 and we obtain at the limit :

∆tχ i ρε -∆tχ i α i -→ 0 and ∆tε ρε -∆tχ i α i -→ ∆t ρ in Ω 0 when ε -→ 0 (31)
as χ(x i ) = 0 for x i ∈ Ω 0 and

∆tχ i ρε -∆tχ i α i -→ - 1 α i and ∆tε ρε -∆tχ i α i -→ 0 in Ω 1 when ε -→ 0. ( 32
)
The pressure update is still p n+1 = p n + p ′ (-µ∇ • u * ). Using (24) to build the velocity update, the standard equation ( 9) is recovered in Ω 0 .

A major advantage of this formulation is that the diffusion coefficient in the pressure equation (30) has generally an absolute value of magnitude ∆t/ρ which avoid the numerical instability of the first-order method. Critical values appears when the interface is very close to the penalized node or a neighbor of the penalized node. If we consider the case where a node x I ∈ Ω 1 has one neighbor x J ∈ Ω 0 , the penalty constraint is α I u I + (1 -α I )u J = u D (x l ), and here Π ′ I u = (1 -α I )u J . Hence, the diffusion coefficient is (α I -1)∆t/(α I ρ) and is critical if the interface Σ h is close to the penalized node. In this case the coefficient is very small and the same problem as with the firstorder correction occurs. However we noticed that a correct solution can be obtained with a diffusion coefficient of magnitude 10 -10 . If the interface position leads to a smaller value, one can slightly move the interface to decrease α I . The case α I ≪ 1 provides a large diffusion coefficient and produces a H 1 -penalty term [START_REF] Angot | A penalization method to take into account obstacles in incompressible viscous flows[END_REF], that is to say a Neumann boundary condition in the considered cell.

Comparison to the Hikeno-Kajishima correction

This approach has been proposed by [START_REF] Ikeno | Finite-difference immersed boundary method consistent with wall conditions for incompressible turbulent flow simulations[END_REF] for a DF-IBM method and is presented here for our interpolator Π. They propose the following velocity update

u n+1 i = u * i - ( (1 -χ i ) ∆t ρ ∇p ′ -χ i Π ′ i ∆t ρ ∇p ′ α i ) (33)
which is directly introduced and not obtained from the solved momentum equation (the consistency of the method is demonstrated afterwards). In our formulation, all the method is derived from the original momentum equation and the prediction equation so the consistency is naturally deduced. From (33) and the incompressibility constraint, the following equation is obtained :

∇ • u * i = ∇ • ( (1 -χ i ) ∆t ρ ∇p ′ -χ i Π ′ i ∆t ρ ∇p ′ α i ) . ( 34
)
and the pressure is updated as in the standard method

p n+1 = p ′ + p n . ( 35
)
The method is applied too to the non-incremental fractional step method.

As can be seen, the equations obtained are different than with the present approach. Obviously, it is due to the way the boundary condition is imposed (penalty term with a coefficient 1/ε of dominant magnitude in our case or IBM direct-forcing term with Dirac function in [START_REF] Ikeno | Finite-difference immersed boundary method consistent with wall conditions for incompressible turbulent flow simulations[END_REF]). However, the final pressure correction equation with the penalized correction (30) tends to (34) when ε tends to zero. The same result is obtained with the corresponding velocity updates.

Remark 4.3 For this correction (and the penalty correction presented here when ε → 0), the velocity in Ω 1 is updated as

u n+1 i = u * i + Π ′ i ∆t ρ ∇p ′ α i ( 36
)
By construction of the interpolator

Π i u = α i u i + Π ′ i u, no node of Ω 1 is involved in Π ′ i .
Hence, as the pressure correction in Ω 0 is

u n+1 i = u * i - ∆t ρ ∇p ′ , ( 37
)
one can replace ∆t ρ ∇p ′ by (u * i -u n+1 i ) to obtain u n+1 i = u * i - Π ′ i (u * i -u n+1 i ) α i . ( 38
)
Considering the initial interpolator Π i , we obtain

Π i u n+1 i = Π i u * i ( 39
)
and the boundary constraint obtained in the predictor step is conserved. It induces as well that Π i u ′ i = 0. The update of the velocity near the velocity is not necessarily null contrary to the linear combination of these velocities.

Temporal accuracy for the incremental pressure-correction

We study here the temporal accuracy of the pressure for the penalized incremental pressure-correction method with a first-order Euler scheme. The temporal accuracy for the base method is O(∆t) [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF]. The ideal equation system to solve is:

[ A G D 0 ] [ u n+1 p n+1 ] = [ r 0 ] ( 40
)
with G the gradient operator, D the divergence operator, A the following sub-matrix

A = { ρ ∆t + N -V - χ ε Π } (41)
with N the linearized discretization of the inertial term and V the discretization of the viscous term. The second member r is defined as

r i = ∆t ρ u * i - χ i ε Π i u * . ( 42
)
The velocity update (29) allows to write

[ I BG 0 I ] [ u n+1 p n+1 ] = [ u * p n+1 ] ( 43
)
with I the identity matrix and B the matrix such as

B = { ∆tχ i ρε -∆tχ i α i Π ′ i ∆t ρ + ∆tε ρε -∆tχ i α i } . ( 44
)
At last, the pressure elliptic equation yields

[ A 0 D -DBG ] [ u * p n+1 ] = [ r 0 ] . ( 45
)
As first noticed by Perot [START_REF] Perot | An analysis of the fractional step method[END_REF], the two systems ( 43) and ( 45) are a block LU decomposition of the system

[ A ABG D 0 ] [ u n+1 p n+1 ] = [ r 0 ] ( 46
)
which allows us to write the error of the scheme using (40) :

AB = { ( ρ ∆t + N -V - χ ε Π ) ( ∆tχ ρε -∆tχα Π ′ ∆t ρ + ∆tε ρε -∆tχα )} . ( 47
)
We study the error in Ω 0 . We consider the matrices

N 0 , N 1 , V 0 and V 1 such as V = V 0 + V 1 and N = N 0 + N 1 .
The matrice N 0 is such that {N 0 } i,j = {N } i,j if x j ∈ Ω 0 , else {N 0 } i,j = 0. The same occurs for V 0 . It is necessary to split these contributions as the value of χ varies in B. Hence, we have

AB = { ( ρ ∆t + N 0 -V 0 ) ∆t ρ + (N 1 -V 1 ) ( ∆t ρε -∆tα Π ′ ∆t ρ + ∆tε ρε -∆tα )} ( 48
)
which tends to term in {1 + O(∆t)} when ε tends to zero. The first order is retrieved. We consider now the error in Ω 1 :

AB = { ( N 0 -V 0 -1 ε Π ′ ) ∆t ρ + ( ρ ∆t + N 1 -V 1 -1 ε α ) ( ∆t ρε-∆tα Π ′ ∆t ρ + ∆tε ρε-∆tα )} = {(N 0 -V 0 ) ∆t ρ + ( -1 ε + ρ ρε-∆tα - ∆tα ε(ρε-∆tα) ) Π ′ ∆t ρ +(N 1 -V 1 ) ( ∆t ρε-∆tα Π ′ ∆t ρ + ∆tε ρε-∆tα ) + ρε ρε-∆tα -∆tα ρε-∆tα } = { (N 0 -V 0 ) ∆t ρ + (N 1 -V 1 ) ( ∆t ρε-∆tα Π ′ ∆t ρ + ∆tε ρε-∆tα ) + 1 } (49)
which also tends to a term in {1 + O(∆t)} when ε tends to zero. Hence, the order of the original method is retrieved. The same result with a quite similar analysis is obtained by [START_REF] Ikeno | Finite-difference immersed boundary method consistent with wall conditions for incompressible turbulent flow simulations[END_REF].

Value of the penalty parameter

The value of the penalty parameter ε has an influence on the solution as the solution u * obtained with the momentum equations converges toward the desired solution for the L 2 -norm with an order 3/4 in ε [START_REF] Angot | A penalization method to take into account obstacles in incompressible viscous flows[END_REF].

For the first-order method, the value of the parameter is more critical. In the present approach, the empirical value ε ≈ 10 -10 is used. This value can varies depending on the linear solver. Nonetheless, taking ε ≈ 10 -20 in (21) will ensure the desired velocity inside the obstacle. For the second-order correction, ε can be taken sufficiently small to ensure Πu * = u D at machine accuracy when the high-order penalty term is used. The converged correction can be also directly implemented.

Numerical experiments

The performances of the modified pressure-correction are evaluated here. The method is compared to the augmented Lagrangian method which does not require a correction. It is also an opportunity to compare the approaches in a more general point of view. When no analytical solution is available, a Richardson extrapolation is used to compute a reference solution [START_REF] Roache | Verification and validation in computational science and engineering[END_REF]. We consider three values h 1 , h 2 and h 3 of a numerical parameter verifying consecutive ratios of two. The convergence rate θ and the reference solution f ext are given by:

θ = ln ( f h 1 -f h 2 f h 2 -f h 3 ) ln ( h1 h2 ) ( 50
)
f ext = ( h2 h3 ) θ f h3 -f h2 ( h2 h3 ) θ -1 (51)
The asymptotical convergence zone has to be reached to obtain a relevant extrapolation.

Cylindrical Couette flow

We consider a Couette flow between two cylinders of radius R 1 = 0.5m and R 2 = 3m. Their angular velocities are ω 1 = 0 rad.s -1 and ω 2 = 2 rad.s -1 . The solution is

v θ (r) = ω 2 R 2 2 -ω 1 R 2 1 R 2 2 -R 2 1 r + (ω 1 -ω 2 )R 2 2 R 2 1 R 2 2 -R 2 1 1 r (52)
for the velocity and

p(r) = ρa 2 r 2 2 - b 2 2r 2 + 2ab log(r) (53)
for the pressure with

a = ω 2 R 2 2 -ω 1 R 2 1 R 2 2 -R 2 1 and b = (ω 1 -ω 2 ) R 2 1 R 2 2 R 2 2 -R 2 1 .
The NS equations are solved in a domain Ω = [ -2 ; 2 + 0.1

√ 2 ] × [ -2 ; 2 + 0.2 √ 2 ].
The analytical solution is imposed on ∂Ω. The SMP method is used to impose a Dirichlet BC on the inner circle. For the augmented Lagrangian method, a value of the parameter dr = 10 is chosen. 1: L 2 relative errors in space and corresponding orders for the velocity and the pressure -augmented Lagrangian an incremental pressure-correction methods

The convergence study of the L 2 relative spatial error is given in table 1 and plotted in figure 2. These results are obtained with the AL method. Similar results are obtained with the rotational pressure-correction with a negligible differential on the L 2 error demonstrating the spatial accuracy of the modification. As expected for such a case [START_REF] Sarthou | Eulerian-Lagrangian grid coupling and penalty methods for the simulation of multiphase flows interacting with complex objects[END_REF], a second order is obtain for the velocity. The convergence rate for the pressure error is around 1.5.

The table 2 gives the same convergence study for the rotational method without immersed-boundary correction for a time step ∆t = 1s. As can be seen in figure 2, the lack of correction has almost no influence on the pressure. The convergence rate for the velocity is lower but acceptable. A factor ten is obtained between the solution with and without the correction for the finest mesh. The time evolution of the solution is now evaluated. For a case with Dirichlet boundary conditions, [START_REF] Guendelman | Coupling water and smoke to thin deformable and rigid shells[END_REF][START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF] give for the rotational method a rate of O(∆t 2 ) for the velocity and O(∆t 3 2 ) for the pressure. Simulations with a 64 × 64 mesh are conducted with different time steps and velocity-pressure coupling methods. The instant T p when the L 2 error on the pressure reaches 1.5 × 10 -3 and the instant T v when the L 2 error on the velocity reaches 1.5 × 10 -4 are considered to study the convergence. The table 4 shows the convergence of these values for the augmented Lagrangian and rotational pressure-correction methods, and for the Euler temporal schemes. The reference values are computed with the Richardson extrapolation using the three more refined values. A clear first order of convergence is obtained for the velocity and the pressure for both velocity-pressure coupling methods. Except for the larger time-steps, both methods gives the same results.

The figure 3 shows the evolution of the spatial error with respect to the number of time iterations for ∆t = 0.1s, ∆t = 1s and ∆t = 10s while the table 3 gives the values of T v and T p for these time steps. Due to its strong implicitation, the AL method is always converged faster than the pressure-corrections. One can notice that the temporal evolution of the solution for the rotational method without correction is quite similar to the evolution of the corrected methods up to an error of 10 -3 on the velocity. For this case, the interest of the correction seems to be minor. The figure 3 shows that the convergence of the incremental pressure-correction is much slower than with the other methods. For the present case, reaching the same level of error than with the two others methods is prohibitive.

The curves of convergence for the second-order Gear scheme are show in figure 4 (the error convergence for the velocity with the augmented Lagrangian method and for the first-order temporal scheme is given for comparison). The values are given in table 5. An irregular convergence is observed for the velocity. For time steps inferior to 0.01s, the error is under the error for the Euler scheme for all values and the asymptotical convergence zone seems to be reached. A first-order of convergence rate is obtained for the pressure. An order of about 0.8 is obtained for the velocity with the pressure-correction which is far from the theoretical order of 2. The augmented Lagrangian method has lower errors but its convergence rate is about 0.5 in the asymptotical zone. Studies conducted in [START_REF] Févriére | On the penalty-projection method for the navierstokes equations with the mac mesh[END_REF] suggest that the augmented Lagrangian should reach a first-order for the velocity and the pressure. The same behavior is obtained in [AJL]. As the convergence rates are better for the first-order temporal scheme, a saturation effect can be involved. Compared to the classical studies, we have the immersed boundary correction here which could cause this saturation. The value of the parameter ε is not involved here as the converged equations in ε gives the same results. In this case the theoretical error study suggests a first-order of convergence for the pressure. The use of a relatively coarse mesh seems not to be involved as we use the Richardson extrapolation to compute the solution so one can suppose that the spatial error does not mix with the temporal error. This point has to be investigated further. The figure 4 shows the convergence for the augmented Lagrangian method with dr = 100. The error on the pressure is close to the other methods. For the smallest time steps, the error on the velocity oscillates around a value so we cannot use the Richardson extrapolation. This values is different from the value obtained with the other methods, so another saturation effect seems to be involved. For this reason, the reference value of the pressure-correction and the augmented Lagrangian with dr = 10 is taken. Even if the convergence is stopped for the smaller time steps, an excellent error (compared to the other methods) is obtained.

To finish with this case, the table 6 shows the spatial errors for various time steps with the rotational method without correction. Contrary to the corrected methods, the error at the stationary state depends on the time step even if its influence is small here. A quite surprising results is that the error decreases when the time step increases while Domenichini [START_REF] Domenichini | On the consistency of the direct forcing method in the fractional step solution of the navier-stokes equations[END_REF] has noticed the contrary (but for different cases and with a spectral solver).

Flow past a cylinder

The instationary flow past a cylinder of unit diameter is now simulated to study the temporal order of the method for an instationary case. We The vorticity and the pressure are shown in figure (5). On can see that the vorticity in the periodic Bénard-von Kármán vortex street is strongly decaying with the X direction compared to the standard solution of the litteraure. This difference is due to the coarseness of the mesh. However, the aim here is not to compare our results with the literacy so the size of the computational mesh is relatively moderate.

The table 8 gives the values of a period of oscillation (adimentionalised by the minimum time step) for different time steps with the augmented Lagragian an rotational methods with the first and second-order Gear schemes for the time derivatives. The convergence order is determined with the Richardson extrapolation performed with the three more refined time steps. The results in term of relative error are given in figure 7.

For both time schemes, the differences between AL and the rotational methods are non-negligible for the larger time steps, with a greater accuracy for the AL. As for the precedent case, it shows the advantage of the AL to deal with larger time steps. For the other time steps, both methodology reach an almost similar accuracy. For the smaller time steps (where the asymptotic convergence seems to be reached) the convergence orders on the velocity and the pressure is about 1 for the Euler scheme and about 1.2 for the second-order Gear scheme. From [START_REF] Guendelman | Coupling water and smoke to thin deformable and rigid shells[END_REF], the rate of error for the L 2 -norm of the velocity of O(∆t 3 2 ) is expected.

Conclusion

The correction of the first-order L 2 -penalty method for the pressure-correction methods has been extended to a second-order with the Sub-mesh penalty method. The correction converges toward the Hikeno and Kajishima [START_REF] Ikeno | Finite-difference immersed boundary method consistent with wall conditions for incompressible turbulent flow simulations[END_REF] correction which is designed for a direct-forcing IBM method. The consistency of the method is directly deduced from its construction.

A brief theoretical analysis has proven that the temporal error of the pressure correction method with a firstorder Gear scheme was not altered. Again, this point is similar to the Hikeno-Kajishima correction. A study with higher integration schemes is now desirable.

Numerical experiments have been carried out. The correction has been compared to the augmented Lagrangian method and the same results have been obtained in space for the cylindrical Couette flow. For this first case, convergence rates of 2 and 1.3 in the L 2 -error norm for the velocity and the pressure have been obtained. It corresponds to the known performances of the rotational method with Dirichlet boundary conditions. Concerning the flow past a cylinder at Re = 100, the study has shown a maximum convergence order for both augmented Lagrangian and rotational methods of 1.20. Those results are close to the literature where a convergence rate between O(∆t) and O(∆t 3 2 ) is expected [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF]. A combination with the second-order open boundary conditions of [PGA11] could be investigated.

As for a small enough penalty parameter ε the present methodology is equivalent to a corrected direct-forcing IBM, especially the method of [START_REF] Tseng | A ghost-cell immersed boundary method for flow in complex geometry[END_REF], the conclusions of this study can be extended to this method.

This work is also a more general comparison between the rotational pressure-correction and the augmented Lagrangian methods (and this last method can a priori be applied to any DF-IBM method). The results show that the spatial accuracy is the same for both methods. Concerning the time accuracy, the AL approach seems to be more efficient with very-high time steps. For moderate and low time-steps, one cannot conclude. Almost no differences have been obtained with the case of the flow past a cylinder. For the Couette flow, quite similar results are obtained for the pressure. For the velocity, the results depends on the time step. For the smaller time steps, the convergence of the AL method decreases while the absolute temporal accuracy is still better than with the rotational method. It has been shown that the value of the penalty parameter dr had an influence on the convergence rate. The influence of the number of sub-iterations could be evaluated too. Increasing the number of sub-iteration generally enhances the convergency at the expense of the computational cost. All these considerations show the complexity of a comparative study between those methods (with and without immersed boundary modification). A future work devoted to an extended comparison would be of high interest, especially if simulations of multiphase flows are performed.

Mesh Relative L²-error

  consider a cylinder of diameter D in a domain Ω = [ -10R ; 15R ] × [ -10R ; 10R ]. The inlet velocity V and the fluid properties are set such as the Reynolds number is equal to 100. The computational mesh is composed of 175 × 150 cells with an inner zone of dimensions [ -D ; 2D ] × [ -D ; D ] with a constant space step covered by 75 × 50 cells. The figure 6 shows the mesh and the position of the cylinder. An Orlanski open boundary condition [Orl76] is imposed for the outflow.

Figure 2 :Figure 4 :Figure 5 :Figure 7 :

 2457 Figure 2: Evolution of the spatial error for the velocity u and the pressure p with (COR) and without (NC) the correction for the Couette flow

  Mesh L 2 rel. error. u Order L 2 rel. error.

				p Order
	16	3.409 × 10 -3		5.898 × 10 -3
	32	4.708 × 10 -4	2.86	4.222 × 10 -3	0.48
	64	1.307 × 10 -4	1.85	1.269 × 10 -3	1.73
	128	3.314 × 10 -5	1.98	3.979 × 10 -4	1.67
	256	8.032 × 10 -6	2.04	1.339 × 10 -4	1.57
	512	2.111 × 10 -6	1.93	5.040 × 10 -5	1.41
	1024	5.281 × 10 -7	2.00	2.040 × 10 -5	1.30
	Table

Table 2 :

 2 Mesh L 2 rel. error. u Order L 2 rel. error. L 2 relative errors in space and corresponding orders for the velocity and the pressure -Timmermans method without immersed boundary correction

	p Order

Table 3 :

 3 ∆t T v AL T p AL T v SPT T p SPT Temporal convergence of the error on the velocity and the pressure for the augmented Lagrangian and the rotational methods for large time steps for the Couette flow

			10	42.27	47.76	439.8	538.2
			1	5.801	5.811	24.18	29.99
			0.1	2.506	2.266	2.573 2.4214
	∆t	T u AL Order		T p AL Order	T u SPT Order	T p SPT Order
	Ref	2.107			1.863	2.107	1.863
	128 × 10 -3 5.168 × 10 -1		5.172 × 10 -1	6.402 × 10 -1	1.034 × 10 +0
	64 × 10 -3 2.591 × 10 -1	1.00 2.584 × 10 -1	1.00 2.835 × 10 -1	1.18 2.968 × 10 -1	1.80
	32 × 10 -3 1.315 × 10 -1	0.98 1.304 × 10 -1	0.99 1.353 × 10 -1	1.07 1.381 × 10 -1	1.10
	16 × 10 -3 6.630 × 10 -2	0.99 6.539 × 10 -2	1.00 6.651 × 10 -2	1.02 6.654 × 10 -2	1.05
	8 × 10 -3 3.328 × 10 -2	0.99 3.268 × 10 -2	1.00 3.299 × 10 -2	1.01 3.259 × 10 -2	1.03
	4 × 10 -3 1.665 × 10 -2	1.00 1.636 × 10 -2	1.00 1.640 × 10 -2	1.01 1.614 × 10 -2	1.01
	2 × 10 -3 8.339 × 10 -3	1.00 8.185 × 10 -3	1.00 8.182 × 10 -3	1.00 8.019 × 10 -3	1.01
	1 × 10 -3 4.177 × 10 -3	1.00 4.096 × 10 -3	1.00 4.082 × 10 -3	1.00 3.984 × 10 -3	1.01

Table 4 :

 4 Temporal convergence of the error on the velocity and the pressure for the augmented Lagrangian and the rotational methods -Gear 1 scheme for the Couette flow

	∆t	T u AL Order	T p AL Order	T u SPT Order	T p SPT Order
	Ref	2.107		1.863		2.107	1.863
	128 × 10 -3 1.918 × 10 -1		8.743 × 10 -1		3.852 × 10 -2	6.724 × 10 -1
	64 × 10 -3 2.531 × 10 -2	2.92 4.095 × 10 -2	4.42 6.041 × 10 -3	2.67 3.887 × 10 -2	4.11
	32 × 10 -3 3.569 × 10 -3	2.83 2.356 × 10 -2	0.80 1.157 × 10 -3	2.38 1.846 × 10 -2	1.07
	16 × 10 -3 2.769 × 10 -4	3.69 1.290 × 10 -2	0.87 1.802 × 10 -3 -0.64 9.048 × 10 -3	1.03
	8 × 10 -3 6.918 × 10 -4 -1.32 6.687 × 10 -3	0.95 1.251 × 10 -3	0.53 4.632 × 10 -3	0.97
	4 × 10 -3 5.336 × 10 -4	0.37 3.428 × 10 -3	0.96 7.212 × 10 -4	0.80 2.326 × 10 -3	0.99
	2 × 10 -3 3.756 × 10 -4	0.51 1.753 × 10 -3	0.97 4.120 × 10 -4	0.81 1.168 × 10 -3	0.99
	1 × 10 -3 2.645 × 10 -4	0.51 8.976 × 10 -4	0.97 2.357 × 10 -4	0.81 5.862 × 10 -4	0.99

Table 5 :

 5 Temporal convergence of the error on the velocity and the pressure for the augmented Lagrangian and the rotational methods -Gear 2 scheme for the Couette flow ∆t L 2 rel. error. u L 2 rel. error. p

	10 -2	7.421 × 10 -4	5.954 × 10 -3
	10 -1	5.593 × 10 -4	1.929 × 10 -3
	1	5.242 × 10 -4	1.322 × 10 -3
	10	5.197 × 10 -4	1.288 × 10 -3
	10 2	5.196 × 10 -4	1.287 × 10 -3

Table 6 :

 6 L 2 relative errors in space for the non-corrected rotational method for various time steps for the Couette flow ∆t/∆t min Period/∆t min (RPC) Order Period/∆t min (AL) Order

	Ref	4.667 × 10 2		4.694 × 10 2
	128	1.137 × 10 3		1.043 × 10 3
	64	7.973 × 10 2	1.02	7.506 × 10 2	1.02
	32	6.282 × 10 2	1.03	6.241 × 10 2	0.86
	16	5.493 × 10 2	0.96	5.530 × 10 2	0.88
	8	5.107 × 10 2	0.90	5.134 × 10 2	0.92
	4	4.896 × 10 2	0.94	4.920 × 10 2	0.96
	2	4.785 × 10 2	0.96	4.806 × 10 2	1.00
	1	4.726 × 10 2	0.96	4.750 × 10 2	1.00

Table 7 :

 7 Period of the oscillations with the rotational pressure correction (RPC) and the augmented Lagrangian (AL) with a Gear 1 temporal scheme for the flow past a cylinder at Re = 100

Table 8 :

 8

	∆t/∆t min Period/∆t min (RPC) Order Period/∆t min (AL) Order
	Ref	4.673 × 10 2		4.691 × 10 2
	128	1.083 × 10 3		1.016 × 10 3
	64	7.386 × 10 2	1.18	7.278 × 10 2	1.08
	32	5.671 × 10 2	1.44	5.762 × 10 2	1.27
	16	5.051 × 10 2	1.40	5.102 × 10 2	1.39
	8	4.832 × 10 2	1.25	4.853 × 10 2	1.36
	4	4.740 × 10 2	1.23	4.760 × 10 2	1.27
	2	4.701 × 10 2	1.18	4.722 × 10 2	1.21
	1	4.685 × 10 2	1.18	4.706 × 10 2	1.21

Period of the oscillations with the rotational pressure correction (RPC) and the augmented Lagrangian (AL) with a Gear 2 temporal scheme for the flow past a cylinder at Re = 100

arthur.sarthou@gmail.com, Institut de Mécanique des Fluides de Toulouse, UMR CNRS 5502, INPT, UPS; 1 Allée Camille Soula, 31400 Toulouse, France

Université de Bordeaux, UMR CNRS 5295, IPB, 16, avenue Pey-Berland,

Pessac Cedex, France

Iterations


On the origins of Conflict of Interest in the Mutual Fund Industry

Keywords: JEL classification : G23, G24, G28 Mutual funds, Broker Dealer, Regulation, conflict of interest Classification JEL : G23, G24, G28 OPCVM, Conseiller financier, Réglementation, Conflit d'intérêt

This article discusses conflicts of interest existing between investors and their financial advisors.

Several contributions in this field treated this question in relation either with invested amount or with expected holding period of the mutual fund. Our paper considers an approach globalizing these two parameters. The first part describes the organization of mutual fund fees in the US. Thanks to a numerical simulation, the second part emphasizes the conflict of interest focusing on the expected holding period. The third one analyses the changes of variables and their consequences on characteristics of investors concerned by this conflict of interest.

On the origins of a Conflict of Interest in the Mutual Fund

Industry

Sébastien M. Lemeunier 1 1 Professor, Department of finance European Business School Paris and Adjunct Professor, Department of Economics, ESSEC Business School Paris Email : slemeunier@gmail.com, Telephone : 0033 620 23 22 10

Introduction

The mutual fund industry in the United States has known a tremendous growth over the last 30 years.

By 2009, half of American households held shares in one or several mutual funds (ICI 2010).

Moreover 8,000 mutual funds manage 11,121 billion dollars in the United-States (ICI (2010)). This success has been supported by substantial marketing efforts and generous commission schemes to induce sales brokers.

To purchase mutual fund shares indeed, investors can turn to many types of intermediaries and have choices as how to pay for their services. They may pay a "sales load" at the time they purchase shares, or a deferred sales load when they redeem shares, and/or they may have to pay annually 12b-1 fees that are deducted from the fund's assets. Combined together, these fees entailed the creation of different share classes typically labeled class A, B and C shares. Investors make a choice between these share classes according to their characteristics in terms of invested amounts and expected holding period. However it may occur that the related compensation scheme of their financial advisor doesn't match their particular needs. In doing so, the existence of this conflict of interest drove to numerous litigation cases (Senator et al (2006) and Krawczyk (2004).

So what are the characteristics of investors exposing them to this conflict of interest? Are investors more or less exposed as determining parameters them to this conflict of interest vary?

Most of studies explain this conflict of interest either with the Expect Holding period O'Neal (1999b), Livingston and O'Neal (1998), Davis (1995) or with the invested amounts Foster (2009), Senator et al 2006, SEC 2004. The latter has been the most frequently treated especially because it entailed many litigation cases. Our approach integrates these two causes to explain this conflict of interest.

To answer these questions, a preliminary part introduces an overview of mutual fund fees in the US.

Thanks to numerical simulations in a second part, it is determined characteristics of investors concerned by a conflicting situation with their financial advisor. It is made the restrictive assumption that invested amounts are lower than $50,000. However, this assumption will be relaxed in the last part. In a third part, varying variables of the model, it is emphasized the evolution of these characteristics.

Mutual fund fees

Investors willing to buy mutual fund shares do not all call for brokers' services and it implies a different pricing structure. It is common practice to distinguish load from no-load mutual funds.

Load funds use brokerage firms to distribute their shares or any other intermediary between them and investors. These intermediaries advise investors in order to guide their choices and to provide additional services2 . Investors care about these guidance and services and it justifies the payment of load fees.

A no-load fund is a mutual fund whose shares are sold without a sales commission and with limited distribution fees3 . No-load funds adopt a more direct relationship with their clients in order to promote their shares. It relies on relatively cheaper distribution methods (advertisement, direct mails…).

In doing so, non-sophisticated investors tend to invest in mutual funds with load fees, under the influence of marketing and their brokers' recommendations ICI (1997), Del Guercio (2002) whereas sophisticated investors choose no-load mutual funds. Thereafter this paper focuses on nonsophisticated investors as it is considered conflicts of interests between investors and their brokers.

According to their invested amount and their expected holding period, investors have to make a choice among different classes of shares which result from combinations of three types of fees. Three main share classes are usually mentioned and typically denoted by A, B, C as described in table 1 (Annex 1):

-Class A shares charge front-end loads declining according to the level of the invested amount and 12b-1 fees relatively lower in comparison to the other share classes. Class A shares are the only one to propose front-end loads. In addition to their profitable aspect, front-end loads induce investors to stay longer in the fund in order to amortize them. Compared to the Class B and C shares described hereafter, funds with Class A shares are more suitable for investors with a relatively longer expected holding period.

-Class B shares are an alternative to Class A shares for investors with a long-term investment horizon. They adopt a coercive approach to prevent from an early withdrawal of the investor.

Instead of paying front-end loads, the investor pays deferred loads in case of redemption of shares.

These deferred loads decrease with each year the share is held. Deferred loads decline to zero over 7 years. Additionally 12b-1 fees for Class B and C shares are of the same order and then relatively higher to those applied to Class A shares. Class B shares are usually converted into Class A shares after the eighth year of investment to avoid that investors bear higher 12b-1 fees too long.

-Class C shares charge relatively higher 12b-1 fees. They also include deferred loads of 1 % if the investor redeems his or her shares within the first year and zero the subsequent years. Mutual funds using Class C Shares are also called "level load funds". Considering the relative combinations of fees they use, Class C shares are particularly suitable for investors with a short term expected holding period. Moreover, they allow a more reactive management of the investor's portfolio changing more frequently of mutual fund at a lower cost.

The conflict of interest analysis

The conflict of interest results from the inadequacy between interests of the investor in terms of fees charged and the related compensation of the broker. In such a situation a broker may be induced to advise a class of share not suitable to an investor according to his or her characteristics in terms of expected holding period or invested amount. In doing so, what are the characteristics of investors concerned by conflicting situations with their financial advisor?

In this part, it is only considered situations where invested amounts are lower than $50,000 to avoid any questions concerning discount policy and to isolate characteristics related to the expected holding period of the investor. Characteristics concerning the invested amount will be treated further.

The model

To demonstrate the conflict of interest, we compare4 the annual Holding-Period Return (HPR) of investors with the Present Value Commission (PVC) of brokers for each share class. A numerical simulation is implemented and results are graphically represented. A first step consists in presenting formulas related to the annual Holding-Period Return and the Present Value Commission. Next simulations are applied according to different expected holding-period and then preferences of investors and brokers are compared.

The annual holding-period return of investors (HPR)

The Holding Period Return is the total net return of a mutual share over the period during which it was held. To compare this HPR between different investment horizons and between share classes we annualize the HPR to obtain a percentage per year.

For each class of share, the related net expected annualized holding-period return is compared in order to estimate the most suitable share class according to the expected holding period of the investor. Each class A, B, C share has a specific combination of fees as described in Table 1 according to the corresponding formulas hereafter:

( )( ) ( ) [ ] { } 1 1 1 1 / 1 250 , , - - × + - = t t A daily daily A t A E R FL HPR
(1)

FL i denotes the front-end loads charged by the fund i, R daily stands for its daily expected returns and E daily is the daily expense ratio5 . t is the holding period, in years.

( ) ( ) [ ] { } ( ) ( ) [ ] ( ) ( ) [ ] { } z z > - - × + × - × + = ≤ - - - × + = t , 1 E 1 R 1 E 1 R 1 HPR t , 1 DL E 1 R 1 HPR 1/t z) - 250(t A daily, daily 250z B daily, daily t B, 1/t t B, 250t B daily, daily t B, (2)
For the Class B share, the formula takes into account the deferred loads (DL) for the corresponding year and the conversion into Class A share after z = 8 years.

( ) ( ) [ ] { } 1 1 1 / 1 , 250 , , - - - × + = t t C t C daily daily t C DL E R HPR (3)
The annual net Holding Period Return is graphically represented using figures appearing in Table 2.

An average gross return of 12% per year has been assumed for this numerical simulation. The chart below displays the evolution of the annual Holding Period Return according to the expected holding period of investors for each share class. Each point on one of these curves corresponds to a date at which an investor withdraws his or her money from a mutual fund. The corresponding coordinate on the Y axis is the average net return an investor could benefit each year until the date of withdrawal. Assuming an average gross return of 12%, an investor that held a fund with a Class A share for 4 years can expect to obtain each year an average net return of 9.2%.

Because HPR A and HPR B are very similar, the comparison will mainly focus on HPR A and HPR C for the relevance of the analysis. The intersection point of HPR A and HPR C means that for the corresponding holding period, an investor is indifferent between having a Class A share and a Class C share. t* inv denotes this intersection point.

Chart 1 shows that an investor with an investment horizon longer than t* inv = 7.9 years would prefer to acquire a Class A share. Conversely, an investor with a shorter investment horizon would prefer to have a Class C share. t* inv , the intersection point of HPR A and HPR C , is solution of the system: * * , , ,,,* 1 1

inv inv t C t A HPR HPR = (4) ( )( ) ( ) [ ] { } ( ) ( ) [ ] { } 1 1 1 1 1 1 1 * * * * / 1 , 250 , , / 1 250 , , - - - × + = - - × + - in inv inv inv t t C t C daily C daily t t A daily A daily A DL E R E R FL ( ) ( ) ( ) ( ) ( ) ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - × + - × + - - - = C daily C daily A daily A daily A t C inv E R E R FL DL t , ,
1 1 ln 250 1 ln (5)
The Present Value Commission of brokers (PVC)

Financial advisors are rewarded for getting a new client through two types of commissions:

-An initial sales commission (IC i ) stemming from front-end loads and paid when the new client enters into the fund and -A trailing sales commission (TC i ) that is paid quarterly to the financial advisor as long as his or her client keeps his or her shares.

As described in Table 1 the financial advisor's earnings depend of the share class sold to the investor.

The Present Value Commission represent the total compensation that the broker receives for the entire duration the investor keeps his or her shares.

The trailing commission for the fund i (TC i ) is paid to the broker on each quarter and deducted from 12b-1 fees. TC i is computed on a daily basis n. Assuming that a quarter includes 63 working days, the following expression gives the trailing commission for the class i share on the quarter m:

[ ] n m m daily i daily daily i m i E R TC y TCquarterl + - = ∑ - + = )) 1 ( 63 63 1 , , , ) 1 )( 1 ( (6)
Then, the expression of the Present Value Commission can be written as the sum of the Initial Commission and the discounted Trailing Commissions. T is the whole length of time in quarters the investor holds his or her shares and k is the discount rate.

Apart from front-end loads applied to the Class A share, the formula used to obtain the Present Value Commission is the same for the three classes of share i.

∑ = + - + = T m m m i i i T i k y TCquarterl FL IC PVC 1 , , ) 1 ( ) 1 ( (7)
To solve this geometric sequence (5), our developments are the following:

∑ = + - - - + = T m m i i m i i i i i i T i k FL IC PVC 1 63 63 63 , ) 1 ( ) 1 ( 1 ) 1 ( β β β β β α which equals ∑ = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + + = T m m i i i T i k IC PVC 1 63 , 1 β γ (8) with [ ] ) 1 )( 1 ( ,daily i daily i E R - + = β , α i = TCquarterly i and ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - - - = 1 ) 1 ( 64 63 i i i i i i i FL β β β β α γ
The present Value Commission of brokers is graphically represented according to the expected holding-period of investors for each share class6 as below. According to these equations and using

∑ ∑ = = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + + = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + + * * 1 63 1 63 1 1 Brok Brok t m m C C C t m m A A A k IC k IC β γ β γ (9)

Divergent / convergent interest areas and characteristics of investors

Comparing the broker's PVC with the investor's HPR according to various investment horizons, the conflict of interest arises indeed.

• A broker could prefer to advise a Class C share to a long term investor, while a Class A share would be more suitable.

• A broker could prefer to advise a Class A share to a short term investor, while a Class C share would be more suitable.

To be more precise a comparison of indifference points of investors and brokers reveals a potential conflict of interest but not for every expected holding period. It may be delimited areas of holding period with a conflict of interest and areas of holding period where interests of brokers and investors converge.

Points of indifference t* inv and t* Brok between classes A and C shares are respectively set after 7.9 years corresponding to a HPR of 10.11% and after 4.4 years corresponding to a PVC of 6.21%.

Based on these key figures, these areas of divergent and convergent interests are shown in the Table 3. As can be seen t* inv and t* Brok may vary with variables of the model. Therefore the convergent interest is unstable and its dynamics should be analyzed.

Dynamics of the convergent interest area

Assuming that investors are uniformly distributed according to their expected holding period, variations of the convergent interest area may give information about the proportion of investors concerned by the conflict of interests. This study may help whether to suggest solutions to enlarge the convergent interest area or to stress situations in which more or less investors are exposed to conflicts of interests.

Convergent and divergent interest areas vary according to exogenous factors as the discount rate, expectations about the market trend, the specific return of a fund, the level of fees and commissions. It is then analyzed effects of these variations on the convergent interest area. The assumption of an invested amount lower than $50,000 will be relaxed at the end of this part.

To keep consistency with the philosophy of share classes, the analysis of the conflict of interest respects the following assumptions:

• Let E A , E C , FL A , IC A , IC C , TC A , TC C , R A , R C ≥ 0 and FL c = 0 Then, one can infer that brokers' risk apprehension represented by an increase of k enlarges the conflict of interest areas at the expense of the convergent interest area. It entails that more unstable is the investment horizon of the client or more uncertain is the advisor's professional situation and more the convergent interest area shrinks.

• Let E A ≤ E C , IC A ≥ IC C and TC A ≤ TC C .
Additionally, the conflict of interest disappears between brokers with a professional time horizon shorter than t* brok and investors with a long term horizon, because he or she will prefer to recommend Class A shares.

-effect of market trend expectations Proposition 2 : If market trend forecasts are optimistic, they imply a higher expected gross return and the convergent interest area increases.

Proof 2 : for an identical variation of R daily,A and R daily,C the expression ,, 1 1 1 1 ln stay the same and t * inv does not change:

( ) ( ) ( ) ( ) ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - × + - × + C daily C daily A daily A daily E R E R , ,
0 * 2 = ∂ ∂ ∂ C A inv R R t . Since
TC quarterly,C is higher than TC quaterly,A and according to equation ( 7), an identical rise of R dailyA and R daily,C will imply a higher increase of the slope for PVC C than for PVC A . t * Brok decreases such as

C A brok R R t ∂ ∂ ∂ * 2 <0 and ∆(t * inv -t * brok ) > 0.
So, optimistic market trend forecasts do not only imply an additional demand for mutual funds (Karceski ( 2002)) but they also imply an enlargement of the convergent interest area. The corollary is that pessimistic market trend forecasts decrease the convergent interest area. Nevertheless, this potential increase of the conflict of interest area is "compensated" by the decrease of the demand that characterizes an expected bear market period.

-effect of a variation of fees and commissions

Effect of a FL A increase

Proposition 3 : If FL A increases the convergent interest area increases. ,,1

Proof 3: equation (5) shows that any increase of FL A implies that -ln(1 -FL

( ) ( ) ( ) ( ) ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - × + - × + C daily C daily A daily A daily E R E R , ,
0 ) ( * * > ∂ - ∂ A brok inv E t t and 0 ) ( * * < ∂ - ∂ C brok inv E t t .
An increase of the daily expense ratio E A and/or front-end loads FL A implies an increase of the convergent interest area. From the conflict of interest perspective, one can say that an increase of Class

A share fees has a virtuous effect9 , all other things being equal. On the contrary, an increase of fees related to Class C share entails a decrease of the convergent interest area.

Effect of brokers' commission increase

Any variation of commissions has no impact on the net average return of investors. Therefore, it is only considered their effect on t* brok .

Proposition 7: an increase of commissions from Class A shares or Class C shares decreases or increases the convergent interest area, respectively. In doing so, considering proposition (6) and reminding that 12b-1 fees are included in the daily expense ratio E i , if an increase of brokers' commissions is funded by an increase of fees in the same proportion, the effect on the convergent interest area is cancelled.

-Evolution of the conflict of interest according to invested amounts The delimitation of these conflict of interest areas has been done assuming that invested amounts are lower than $50,000. Nevertheless, depending on the importance of the amount invested, Class A shares offer discounted loads. It is then examined if investors with larger amounts will be more or less exposed to the conflict of interest. It is applied to front-end loads the sliding scale according to the invested amount and breakpoints as described in Table 1.

Additionally, while it is observed the evolution of the convergent interest area for class A and C shares, it is included the Class B share. However, if the invested amount is higher than the breakpoint $100,000 class A shares offer a discount to investors. For this level of invested amount, investors with a shorter expected holding period prefer to buy Class A Shares. On the other hand, brokers prefer advising a Class B Share to have higher From $50,000, brokers are indifferent between class A and B share and from $100,000 they systematically are more rewarded when they sell a class B share rather than selling a class A share.

From this breakpoint, higher invested amount are and narrower is the convergent interest area.

Additionally, including Class B Share, the convergent interest disappears for lower invested amount $500,000 in our example).

This last diagram stresses that investors the most exposed to the conflict of interest exposure are those with the largest amounts. This conclusion takes into account a comparison between class A, B and C shares and explains the conflict of interest with the expected holding period and the invested amount.

Conclusion

It has been demonstrated the existence of a conflict of interest between investors and their financial advisor. To be concerned by these potentially conflicting situations, we have determined characteristics of investors in terms of expected holding period for their mutual fund shares and invested amount.

Isolating investors with an invested amount lower than $50,000 , it has been stressed with numerical simulations the existence of a convergent interest area and areas of conflict of interest according to their expected holding period. For these amounts lower than $50,000 it has been focused on a comparison between Class A and C shares. Varying variables of the model, we show that characteristics of investors exposed to the conflict of interest change.

Any increase of market trend expectations or Front end loads or daily expense ratio of class A share implies an increase of the convergent interest area.

An increase of the discount rate of brokers or daily expense ratio of Class C share shrinks the convergent interest area.

Integrating invested amounts and comparisons with class B shares, it appears that the convergent interest area shrinks for amounts higher than $100,000 and disappears for amounts higher than $500,000. From this breakpoint, any kind of investor is exposed to a potential conflict of interest with his or her financial advisor. 12b-1 fees are an annual marketing or distribution fee on a mutual fund. They are paid every year as long as the investor holds his or her shares.

Bibliography

Deferred loads : also known as a back-end sales charge, they are imposed when an investor redeems shares. The percentage charged declines the longer shares are held.

Other expenses: annual fees non related to distribution expenses.

Front-end loads are sales charge applied at the time of the initial purchase for a mutual fund. It is deducted from the investment amount.

Annex 2

Proof 7 :

 7 according to equation(8), an increase of IC A or IC C leads to a parallel upward shift of PVC A or PVC C and t * brok increases or diminishes, respectively: (7) an increase of TC A or an increase of TC C entails a steeper PVC A curve or a steeper PVC C curve, respectively.

  As a result of a numerical simulation, this chart shows the evolution of indifference points between class A and C shares and the conflict of interest according to invested amounts. The dotted line represents the expanse of the convergent interest area between brokers and investors. It steadily decreases to equal zero for amounts exceeding 1 million dollars. Surprisingly, investors with the largest amounts are also the most concerned by this conflict of interest. In other words, it is shown that the most important clients are and the most they are exposed to the conflict of interest with their financial advisor.

  commissions 10 . Hereafter, it is shown the evolution of points of indifference Chart 4 : Evolution of points of indifference between class A

10

  However the SEC have prohibited this type of practice mentioning it into an amendment of the Rule 12b-1 in 2004. Securities and Exchange Commission [2004] « Prohibition on the Use of Brokerage Commissions to Finance Distribution."; SEC [release No. IC-26591; File No. S7-09-04] October 14, 2004

•

  Alexander, G., Jones, D., Nigro P., 1998. Mutual fund shareholders: characteristics, investor knowledge, and sources of information. Financial Services Review 7, 301-316 • Capon N., Fitzsimons G., Prince R., 1996. An individual level analysis of the mutual fund investment decision. Journal of Financial Services Research 10, 59-82. • Collins S., 2004. The effect of 12b-1 plans on mutual fund investors. Working paper, ICI. • Davis L., 1995. A New Look at SEC Rule 12b-1. Securities Regulation Law Journal, Vol 23. • Del Guercio D., Tkac P., 2002. The determinants of the flows of funds of managed portfolios: Mutual funds vs. pension funds. Journal of Financial and Quantitative Analysis Vol. 37 N°4 pages 523-557. • Foster P., 2009. Mutual Fund Multi-Class Offerings: Addressing Conflicts of Interest Through Meaningful Disclosure and Robust Sales Practice Protocols. The investment Lawyer Vol. 16, No. 6 III, ed. Financial Services Perspectives and Challenges (HBS Press, Boston, Mass). • Walsh L., 2004. The costs and benefits to fund shareholders of 12b-1 plans: an examination of fund flows, expenses and returns. ; SEC; May 2004 • Wilcox R. T., [2003]; "Bargain hunting or star gazing? Investor's preferences for stock mutual funds"; Journal of Business Vol 76 n°4; 2003 • Zhao Xinge [2003]; "The role of brokers and financial advisors behind investments into load funds"; Working paper; November 2003 Annex 1

Present Value Commission (PVC) of brokers with a risk free rate

  figures in Table 1, we conduct a numerical simulation as shown in Table 2. Calculations are based on a gross annual return of 12% and a risk free rate of 3% has been applied to discount commissions of brokers. In this Chart, an intersection point is also noticeable between the Present Value Commissions related to classes A and C shares. This point is denoted by t* Brok representing the holding period for which brokers are indifferent between advising a class A share and a Class C share. t* Brok appears after 4.4 years corresponding to a PVC of 6,21%. t* Brok is solution of the system 7 :

									PVC	A	,	t	* Brok	=	C PVC	,	t	* Brok
	5% 40% 45% Chart 2 : 0% 10% 15% 20% 25% 30% 35% Present Value Commission																	PVC A PVC B PVC C
	0	5	10	15	20	25	30	35	40	45		50	55	60	65	70	75	80
						Holding period in quarters

Table 3 :

 3 Conflict of interest between

brokers and advisors Area with a conflict of interest Area with convergent interests Area with a conflict of interest Points of indifference

  For the period of time bounded by the segment [t* inv ; t* Brok ] interests of brokers ans nvestors converge. This area with convergent interests is spread over a period of t* inv -t* Brok = 3.5 years.In doing so, investors characterized by an expected holding period included in [t* inv ; t* Brok ] are not concerned by this conflict of interest. Out of this area, investors are all concerned.

	0	t *	Brok = 4.4 years	t *	inv = 7.9 years
	Brokers prefer to advise	A	C		C
	Is preferable for Investors	C	C		A

  An increase of the daily expense ratio of the class A or an increase of the daily expense ratio of the Class C share increase or reduces the convergent interest area, respectively.

	diminishes. It entails a reduction of t *	brok . Therefore	* brok A FL ∂ ∂ t	<	0	and	* inv A FL ∂ ∂ t	>	0	, hence for an
	increase of FL A ,	∆	( t	* inv	-	t	* brok	)	>	0
	Effect of a daily expense ratio increase 8
	Proposition 4: Proof 4:
											A ) raises and then t *	inv
	increases. If FL A increases, γ A decreases and according to expression (8) the slope of PVC A

Table 1 Typical Class structures for Load Mutual Funds Typical class structure for Load Mutual Fund Broker compensation arrangement Front-end Share class 12b-1 fees Other expenses Deferred load Initial investment Load Initial commission Trailing commission

 1

					< $50,000	5.75%	5%
					$50,000 -$100,000 4.50%	4%
	A	0.25%	0.70%	0%	$100,000 -$250,000 3.50% $250,000 -$500,000 2.50%	3.20% 2.25%	0.25%
					$500 000-$1,000,000	2%	1.70%
					>1 000 000$	0%	1%
				5% the 1st
	B	1%	0,70%	year 4% the 2 nd year 3% the 3rd year	0% Conversion
				3% the 4th year
				2% the 5th year
				1% the 6th year

into Class A share after the 8th year of investment.

  Neal E. S.[1999];"Mutual fund share classes and conflict of interest between brokers and investors"; Working paper; 1999 O'Neal provides his article with data obtained from an average sample composed of 20 equity mutual funds among the largest.

						4%	0.25%
				1% the first
	C	1%	0.70%	year 0% after the first	0%	1%	1%
				year
	Cf: O'

Table 2 Numerical simulations of investor Annual Holding-Period Returns and Broker Present Value Commissions Annual Holding-Period Return Present Value Commission Investment Duration in

 2

	years	HPR A,t	HPR B,t	HPR C,t	PVC A,t	PVC B,t	PVC C,t
	0	0%	0%	0%	5,00%	4,00%	1,00%
	1	4.56%	5.11%	9.11%	5,25%	4,26%	2,03%
	2	7.70%	8.28%	10.11%	5,51%	4,54%	3,14%
	3	8.77%	9.28%	10.11%	5,80%	4,83%	4,32%
	4	9.31%	9.55%	10.11%	6,10%	5,15%	5,59%
	5	9.63%	9.84%	10.11%	6,43%	5,49%	6,95%
	6	9.85%	10.01%	10.11%	6,79%	5,86%	8,39%
	7	10.01%	10.11%	10.11%	7,18%	6,24%	9,95%
	8	10.12%	10.11%	10.11%	7,59%	6,69%	11,60%
	9	10.21%	10.20%	10.11%	8,04%	7,16%	13,38%
	10	10.29%	10.28%	10.11%	8,52%	7,67%	15,28%
	11	10.35%	10.34%	10.11%	9,04%	8,22%	17,31%
	12	10.39%	10.39%	10.11%	9,60%	8,82%	19,48%
	13	10.44%	10.43%	10.11%	10,21%	9,46%	21,81%
	14	10.47%	10.47%	10.11%	10,86%	10,15%	24,30%
	15	10.50%	10.50%	10.11%	11,56%	10,90%	26,96%
	16	10.53%	10.53%	10.11%	12,32%	11,70%	29,81%
	17	10.56%	10.55%	10.11%	13,14%	12,57%	32,85%
	18	10.58%	10.57%	10.11%	14,02%	13,50%	36,11%
	19	10.60%	10.59%	10.11%	14,96%	14,51%	39,60%
	20	10.61%	10.61%	10.11%	15,99%	15,59%	43,33%

This advice includes the fund selection, the asset allocation of the investor's portfolio and taxation. Services usually include the opening of an account, financial orientated press, an availability to answer to investors' questions, an overview of the account by web…

Distribution fees, also called Rule 12b-1 fee can not exceed 0.25 percent per year. (ICI 2009 Mutual fund fact book)

This approach comparing share classes has been used byCollins (2004),O'Neal (1999a and 1999b), Livingston and O'Neal (1998) andDavis (1995).

The daily expense ratio results from the annual expense ratio comprising 12b-1 fees and management fees charged to the investor every year and as long as he or she holds his or her shares of mutual fund.

Class A and B shares are also very similar from brokers' point of view. It is then focused on the comparison between A and C share classes.

The existing conditions for the intersection point t * brok are IC A > IC C and γ C > γ A .

It is reminded that E daily , i is the daily expense ratio that includes "other expenses" and 12b-1 fees.

It can also be demonstrated that E A has a larger impact than FL A on the convergent interest area.