K R Krishnan

David Shallcross
email: davids@research.telcordia.com

Latha Kant
email: lkant@research.telcordia.com

Joint Optimization of Scheduling and Multicast Trees by Column-Generation

Keywords: wireless networks, capacity, routing, scheduling, multicast

In a wireless network, transmissions from the various nodes have to be scheduled so as to avoid mutual interference. The pattern of interference induced by active transmissions depends on the routes along which the link-transmissions have to be scheduled; i.e., the interference that scheduling has to accommodate depends on the routing. Both routing and scheduling are mechanisms to promote the efficient use of network capacity, and in view of their interdependence, it is important to consider their joint optimization. Considering one without the other can create limitations for each function and is non-optimal. We develop an algorithm based on the columngeneration technique of Linear Programming for the joint optimization of routing and scheduling for multicast flows for maximizing network capacity, and demonstrate the benefit of the joint optimization in increased capacity over the case where the routing and scheduling are separately considered.

INTRODUCTION

The intrinsic 'broadcast' nature of wireless transmission makes multicasting a natural application for wireless networks. In this paper, we develop an algorithm based on the columngeneration technique of Linear Programming for the joint optimization of routing and scheduling for multicast flows for maximizing network capacity, and demonstrate the benefit of the joint optimization in increased capacity over the case where the two functions are separately considered.

In a TDMA (Time Division Multiple Access) wireless network, transmissions from the various nodes have to be scheduled so as to avoid mutual interference. Thus, traffic is sent from origin to destination by scheduling the flow on each hop of its route. The pattern of interference created by active transmissions depends on the routes along which the linktransmissions have to be scheduled; i.e., the interference that scheduling has to accommodate depends on the routing. Owing to this connection between routing and scheduling, it is important, for efficient network use, to consider the joint optimization of routing and scheduling. Considering one without the other can create limitations for each function and is non-optimal, as was shown in [START_REF] Krishnan | Joint optimization of transmission schedule and routing for network capacity[END_REF] for the case of unicast flows.

The 'capacity' of a wireless network with respect to multicast flows has been investigated in several recent papers [2 -7], which extend the unicast-capacity results of Gupta and Kumar [START_REF] Gupta | The capacity of wireless networks[END_REF] to derive bounds on multicast flows under certain formulations of what is meant by 'capacity', with assumptions on network topology and loads. In contrast to the notion of capacity formulated in those papers, we work with what might be termed an 'operational' definition of network capacity, stated with reference to the actual traffic demands placed on a given network. According to this 'operational' view, maximizing the capacity of a network means minimizing the fraction of active transmission time needed in the network to carry the specified traffic demands, or, equivalently, maximizing the inactive fraction of transmission time. The inactive fraction of time can then be thought of as the 'excess' capacity available in the network after meeting the traffic demands placed on it. Our interest in this paper is in an algorithm for the joint optimization of paths and linktransmission schedules for multicast flows in a network so as to maximize the excess capacity available in the network, using the conventional implementation of multicast in which an intermediate node replicates and forwards the information it receives to successor nodes, and not the more general case of 'network coding'.

A multicast session is described by the triplet ) , , (

λ D s
, where s is the source-node, D is the set of destination nodes, and λ is the rate of the flow from s to each of the destinations in D . A route for a multicast flow consists of a 'tree' rooted at s , with leaves at the destination nodes D . The routing we consider allows for the traffic of a multicast session to be split among multiple trees. A direct formulation of the joint optimization problem becomes intractable in various situations in networks of realistic size. If a solution is feasible with a single tree for each multicast session, the optimum choice of tree becomes a mixed-integer program. When multiple trees have to be used for flows, the problem can become intractable owing to the large number of trees one may have to consider. Our approach is to formulate the joint optimization problem as a Linear Program and then investigate the possibility of using column-generation [START_REF] Bazaraa | Linear Programming and Network Flows[END_REF] successively to improve on an initial choice of multicast trees for the multicast sessions. In particular, at each step, a solution to the LP allows us to determine, with the use of dual variables, the conditions under which column-generation yields a solution that can improve on the multicast capacity of the initial choice of multicast trees. Thus, the dual variables guide us in the discovery of new candidates for paths for multicast flows in conjunction with the requirements of scheduling, avoiding the need for generating a large number of trees without the assurance that they contribute to improving a solution. We demonstrate the significant capacity gains made by our method in three examples of multicast flows in a 10-node network over the initial case of an optimized transmission schedule for an arbitrary choice of multicast tree.

Previous Work

One finds various approaches in the literature to the problem of finding paths in a wireless network for multicast traffic and for scheduling hop-by-hop transmissions. One approach is to determine [START_REF] Penttinen | Efficient multicast tree algorithm for ad hoc networks[END_REF][START_REF] Penttinen | Minimum cost multicast trees in ad hoc networks[END_REF], for each multicast session, a suitable directed tree in the network, rooted at the source, whose leaf-nodes are the destinations of the multicast. It is possible that nodes other than the source-node and the destination-nodes also appear in the tree, serving as intermediate relay nodes. Thus, the multicast tree is a Steiner tree in the network. Algorithms, exact or heuristic, can be used for determining 'minimum-cost' Steiner trees, where, for example, 'cost' could be the total Euclidean length of links in the tree. However, selecting paths without reference to the need to schedule transmissions to minimize interference could lead to non-optimal choices, as was shown in [START_REF] Krishnan | Joint optimization of transmission schedule and routing for network capacity[END_REF] for the case of unicast flows. Another approach is to select, at regular intervals, the multicast tree to from a set of pre-specified trees [START_REF] Sarkar | A Framework for Routing and Congestion Control for Multicast Information Flows[END_REF], so as to minimize queue lengths in buffers at the time of selection. A third approach is not to bother with constructing a tree at all [START_REF] Toumpis | Performance Bounds for Large Wireless Networks with Mobile Nodes and Multicast Traffic[END_REF], but rely, instead, on network connectivity and a transmission scheme that uses 'minislots' within each time-slot to relay each packet to its various destinations in bounded time. The paper [START_REF] Chaporkar | Wireless Multicast: Theory and Approaches[END_REF] considers a protocol for scheduling hop-by-hop multicast transmission in which a sending node is allowed to transmit whenever the number of its non-interfering receivers exceeds a certain threshold.

After developing the traffic equations for multicast flows in Section 2, we consider in Section 3 the formulation of the constraints introduced by the need for scheduling transmissions to avoid interference. However, there is an important difference between a unicast flow and a multicast flow as to what constitutes interference. Wireless transmission is inherently broadcast in its nature, since a single transmission by a node can be 'heard' by many neighbors, which, in fact, could be a helpful feature in directing multicast flows. This point is discussed in Section 3 and leads to the definition of a 'multicast link'-a set of links belonging to a multicast tree and incident at a common node. Their simultaneous transmissions advance the multicast flow and do not constitute interference with one another, and can be treated as a single 'link' for the purpose of the multicast.

Then, in Section 4, we first examine a Linear Program formulation for optimizing the transmission schedule when the routing is pre-specified, and take up the joint optimization of routing and transmission schedule in Section 5. In Section 6, we show the application of column-generation and dual variables to the joint optimization problem, and state the proposed algorithm in Section 7. Section 8 presents numerical results for several examples to show the effectiveness of the proposed method, and our conclusions are summarized in Section 9.

TRAFFIC FLOW EQUATIONS AND CONSTRAINTS

Suppose there are N nodes, L links, and n multicast sessions in the network. The traffic k λ of multicast session k

uses k τ different trees ] , , [ 1 k k k t t τ L
, with the traffic being divided among those trees in the proportions

} , , { 1 k k k τ α α L , with 1 1 = ∑ α τ = k j kj . Let network. in the used trees multicast of number total 1 = ∑ τ = = n k k T
The mapping of the traffic from the multicast trees onto links of the network is specified by the following expression for load l f on link l :

       = = δ ∑ ∑ α λ = = τ = otherwise , 0 session multicast of tree the of part is link if , 1 where ,..., 1 , 1 1 k j l δ L l f th k lj k lj n k j kj k l k (1).
We put (1) in matrix form with the aid of the following definitions: Let

T n n T k k k T L n k n k f f f ] ,.
τ τ τ τ α α α α α α = α = α α = α = L . For n k ,..., 1 = define the k L τ × matrix ) ( ) ( k lj k k lj k b B δ λ ≡ = . Let ] [ 1 n B B B M L L M =
, an T L × matrix. Then (1) can be rewritten as

α = α ∑ = = B B f k n k k 1 (2).
We represent the set of constraints

n k k j kj ,..., 1 , 1 1 = = ∑ α τ = in the matrix form E H = α (3),
where H is an T n × matrix and E is an n -vector of 1's.

TRANSMISSION SCHEDULE

We let S denote the set of independent link-sets considered for use in the network, where the links of an independent linkset can all be active at the same time without causing mutual interference. Let M S =

. The problem of generating all independent sets is itself NP-hard, and in the numerical examples of a 10-node network presented later, we used the set of all maximal independent link-sets as a 'reasonably complete' set. Let s φ be the fraction of a frame-duration for which independent link-set s is active, with

1 1 ≤ ∑ φ = M s s
as a necessary condition of feasibility of the given loads in the given network under the given routing. With the frame-duration as the unit of time, the fraction of time a given link l is active is the sum of the active durations of each independent link-set to which the link belongs. Thus, if l g denotes the maximum rate at which link l transmits whenever it is active, we have the constraint We note that there is an important difference between a unicast flow and a multicast flow with respect to interference among link-transmissions. For example, in Figure 1, Node A has three incident links. If A is transmitting three separate unicast flows to B, C, and D respectively, these transmissions have to be scheduled for three different time-slots, since, in our model, a node can only handle a single flow in each time slot, and hence only ONE of the links (A→B), (A→C), and (A→D) can be active at any given time. On the other hand, the transmission of a multicast flow by a node on two or more of its incident links is part of the broadcast mechanism by which the multicast flow reaches its destinations, and does NOT constitute interference. Hence, such transmissions can be simultaneous, i.e., for a multicast transmission from node A to nodes B, C, and D, all three links (A→B), (A→C), and (A→D) can be scheduled for the same time-slot. Such a set of links with a common transmitter is defined to be a 'multicast link', and can be treated as a single link for the purpose of a multicast. Thus, interference for multicasts is a less restrictive condition than for unicasts.

   ∈ = ∑ φ ≤ = otherwise , 0 set - link if , where 1 s l g y y f l ls M s s ls l or, in matrix-form, φ ≤ Y f (4).

A. Notion of Interference for Unicast and Multicast Flows

CAPACITY MAXIMIZATION WITH GIVEN MULTICAST TREES

Given a set of multicast sources in the network with traffic rates

JOINT OPTIMIZATION OF TRANSMISSION SCHEDULE

AND MULTICAST TREES When one has to determine the optimum choice of multicast trees for the traffic rather than just choose from a given set of trees, the problem becomes much harder and is now a mixed-integer program. Rather than taking a direct approach to this intractable problem, we seek to use the technique of 'column generation' in which, starting from a small initial set of 'columns' of the constraint matrix (in our case, the columns correspond to multicast trees), we use dual variables to guide a search for a column that would offer improvements to the objective function over the existing columns. The method also has the advantage of making it clear when no such improvement is possible by augmentation of the existing columns.

The dual of the LP ( 5) is constructed in terms of the following dual variables: When the primal and dual LPs are both feasible, they have the same optimum objective value, and the solution of the LP also produces values for the dual variables } , { w v . The next step is to use the dual variables to define a condition that has to be satisfied by any new column (multicast tree) that can offer an improvement in the objective of the primal LP.

COLUMN GENERATION

Suppose that we solve the primal LP ( 5) and its associated dual LP [START_REF] Chaker | Multicast Capacity of Wireless Ad Hoc Networks with Hierarchical Routing[END_REF]. Let the optimum values of the variables be denoted with the superscript *. Then, the 'cost' of the optimum solution of ( 5) is given by

Minimum cost = ∑ = M s s 1 * φ (7).
Suppose that we add another multicast tree t for a specific multicast session

n k k ≤ ≤ 1 ,
. Of course, this would require also that the condition Let the structure of this new multicast tree be denoted by the column-vector

T L b b b ] , , [ 1 L = , where    = otherwise , 0 tree the to belongs link if , 1 t i b i
The new B -matrix corresponding to the addition of tree t is given by [ ]

n k k k B b B B B B L M M M L M M λ 1 2 1 -

, and the corresponding new H matrix has a new column

T n h h h ] , , [ 1 L = in the position         + ∑ τ = 1 1 k j j , with    ≠ = = k i k i h i for , 0 for , 1
From the interpretation of the dual variables [START_REF] Bazaraa | Linear Programming and Network Flows[END_REF] determined by the solution to [START_REF] Chaker | Multicast Capacity of Wireless Ad Hoc Networks with Hierarchical Routing[END_REF], it follows that the addition of tree t can lead to a smaller cost than [START_REF] Chaporkar | Wireless Multicast: Theory and Approaches[END_REF]

only if , i.e., 0 * 1 * 1 * * k k L i i i L i k i i k w b v w b v λ λ < < - ∑ ∑ = = (8).
Condition [START_REF] Gupta | The capacity of wireless networks[END_REF], in fact, implies two requirements:

1. 0 * > k w 2. We are able to construct a tree for multicast session

k with link weights * i v such that k k L i i i w b v λ < ∑ = * 1 *
In the absence of degeneracy [START_REF] Bazaraa | Linear Programming and Network Flows[END_REF], when these two requirements are met, the addition of multicast tree t for multicast session k offers an improvement over the solution [START_REF] Chaporkar | Wireless Multicast: Theory and Approaches[END_REF]. In particular, if, at some point, it is found that no tree can be constructed that meets the criterion (8), we know that an optimal solution has been found (however, the guarantee of optimality does not hold if inexact methods were used to decide that such a tree does not exist).

ALGORITHM FOR CHOICE OF MULTICAST TREES

We now describe the overall procedure followed in applying column-generation in the example to be presented:

1. First, enumerate the set of independent link-sets to be used in the network. For the 10-node network shown below, we, in fact, generated all the maximal independent link-sets.

2. Determine an initial set of multicast trees for the given multicast demands, using unit 'cost' on all links in this initial step. In later iterations, the 'costs' assigned to links will be prescribed by the dual variables on the LP solution of the previous iteration. The 'cost' of a path is the sum of the 'costs' of the links on the path.

Since finding a minimal Steiner tree is itself an NP-hard problem, we use the following simple heuristic algorithm. The same algorithm is used later when constructing trees as part of the column-generation step.

Algorithm for Construction of Steiner Tree

To find a low-total-weight tree connecting a single source to multiple destinations: While there is a destination not yet connected, find a minimum-cost path from the tree constructed so far to a destination not yet connected, and add this path to the tree. b) Extract the dual variables. Assign new link-costs and multicast tree cost-thresholds from the dual variables, in accordance with [START_REF] Gupta | The capacity of wireless networks[END_REF].

c) For each multicast demand, search for low-cost multicast trees (using the algorithm above), and add to the current set of multicast trees any new trees found that have total costs less than the costthreshold in [START_REF] Gupta | The capacity of wireless networks[END_REF] for the corresponding multicast demand.

d) If no trees are found below the required costthreshold (or if the number of iterations has reached the pre-set limit), exit the loop, and take the solution to be completed. In each case, the process of column-generation was started by solving the LP (5) with an arbitrary tree for each multicast, and then repeating the LP after adding trees that were found admissible by the criterion (8), making use of the dual variables corresponding to the LP solution. The iterative procedure stops either when no new trees can be found that meet [START_REF] Gupta | The capacity of wireless networks[END_REF] or when the number of iterations has reached a preset limit (30, in our example).

RESULTS

As noted earlier, a measure of the spare capacity in the network is given by 

       -∑ = M s s 1 1
φ , which is the metric that we seek to maximize. The results obtained for the three cases are shown in Table 1. The effectiveness of the method is borne out by the difference between the initial spare capacity (corresponding to the initial choice of a multicast tree) and the final spare capacity that produced at the end of the columngeneration iterations. Table 1 also shows, for each case, the number of iterations taken to reach the final solution, as well as the number of trees and the number of independent link-sets that are used in that solution.

SUMMARY

The need to schedule link transmissions in wireless networks to avoid interference among links leads to limitations on the capacity of the network that are quite separate from the inherent capacity of the links themselves. The choice of routing also determines how well the capacities of various paths are put to use in carrying traffic. Scheduling and routing interact in their effect on capacity -links are scheduled so as to avoid interference with other links, while the interference pattern that scheduling has to accommodate depends on traffic routing. As a result, for efficient use of the network, routing and scheduling should be considered together: considering one without the other is non-optimal. This was shown in [START_REF] Krishnan | Joint optimization of transmission schedule and routing for network capacity[END_REF] for the case of unicast flows.

In this paper, we have dealt with the case of multicast flows. Here, the joint optimization involved is that of scheduling and the determination of multicast routes (Steiner trees, with root at multicast source and leaf-nodes at destinations). A direct formulation of this problem leads to an intractable mixed-integer program. We have proposed an iterative method, based on LP duality theory and columngeneration, in which, in each iteration, we derive conditions for higher capacity to be achievable by a new choice of tree. If such a tree cannot be found, the current solution is guaranteed to be optimal. This guarantee of optimality holds when exact methods are used to check for the existence of trees that meet the required conditions. If non-exact methods have to be used for tree-generation (since tree-generation is itself a hard problem), optimality at termination is not guaranteed. In practice, one also uses a pre-set bound on the number of iterations. In any case, at each step, we know whether the next Link-bandwidth = 10 Mb/s for all links Traffic Rate for Multicast Sessions = 2 Mb/s step can yield a better solution or not. Our examples show the efficacy of the method and the significant capacity-gain that results from joint consideration of routing and scheduling.

Figure 1 .

 1 Figure 1. Interference for Unicast and Multicast Flows

  the problem of capacitymaximization as the problem of supporting the traffic of the given sources while minimizing the active time of the network the links are inactive, as a measure of spare capacity in the network beyond that required to carry the given the operational definition of capacity that we adopt in this paper. Thus, for a given choice of multicast

3 .

 3 Repeat the following steps (a), (b), (c), and (d) for a specified maximum number of iterations: a) Solve the LP (5) with the given independent linksets and the current choice of multicast trees.

Figure 2 . 10 -

 210 Figure 2. 10-Node Network used in Example The algorithm was applied to the 10-node network shown below in Figure 2, for three different cases of multicast traffic demands. The bandwidths of all the links were taken to be 10 Mb/sec and the traffic rate of each multicast session in the network was taken to be 2Mb/sec. The multicast traffic sessions used in the three cases are as follows: Case I Multicast T1 = [Node 10 → Nodes (6, 9)] Case II Multicast T2 = [Node 7 → Nodes (2, 6)] Case III Multicast T3 = T1 + T2 = [Node 10 → Nodes (6, 9)] + [Node 7 → Nodes (2, 6)] In each case, the process of column-generation was started by solving the LP (5) with an arbitrary tree for each multicast, and then repeating the LP after adding trees that were found admissible by the criterion[START_REF] Gupta | The capacity of wireless networks[END_REF], making use of the dual variables corresponding to the LP solution. The iterative procedure stops either when no new trees can be found that meet[START_REF] Gupta | The capacity of wireless networks[END_REF] or when the number of iterations has reached a preset limit (30, in our example).As noted earlier, a measure of the spare capacity in the

TABLE 1

 1

		.	MULTICAST CAPACITY WITH COLUMN-GENERATION
		Spare Capacity =
	Case		    -∑ = M s 1 1	φ	s	   	# Iterations	# Independent Link-Sets	# Trees
		Initial Final Increase
	I	0.4	0.60		0.20	18	4	2
					(50%)
	II	0.2	0.63		0.43	18	7	3
					(215%)
	III	0.0	0.36		0.36	30	13	10
					(∞)

WNC3 2010

ACKNOWLEDGEMENT

This research was carried out by collaborative participation in the Communications and Networks Consortium sponsored by the U. S. Army Research Laboratory under the Collaborative Technology Alliance Program, Cooperative Agreement DAAD19-01-2-0011. The U. S. Government is authorized to reproduce and distribute reprints for Government purposes notwithstanding any copyright notation thereon.


Claire Bost

Georges-Henri Cottet

G.-H Cottet

E Maitre

Convergence analysis of a penalization method for the three-dimensional motion of a rigid body in an incompressible viscous fluid

come

Convergence analysis of a penalization method for the three-dimensional motion of a rigid body in an incompressible viscous fluid

Introduction

In this paper we are concerned with the numerical analysis of a penalization method for the twoway interaction of a rigid body with an incompressible fluid in three dimensions. The traditional numerical approach to deal with fluid-structure problems is the so-called ALE (for Arbitrary Lagrangian Eulerian) method where fluids (resp solids) are described in an Eulerian (resp Lagrangian) framework. Fluids are computed on a moving mesh fitting the solids and stress and velocity continuity are used to derive the appropriate boundary conditions on the fluid/solid interface [START_REF] Maury | Direct simulations of 2D fluid-particle flows in biperiodic domains[END_REF].

Another way is to use a fixed mesh and rely on a Lagrange-Galerkin approach where the fluid is solved using a characteristics-finite-element method in a space where the velocity field is rigid in the solid, and the solid itself moved by a Lagrangian equation. In [START_REF] San Martín | Convergence of the Lagrange-Galerkin method for the equations modelling the motion of a fluid-rigid system[END_REF] the convergence of a scheme following this approach is proved in the case where the densities of fluid and structure are equal.

Alternate methods can be devised where the whole fluid-solid system is seen as a multiphase flow and the fluid/solid interface is captured implicitly rather than explicitly. Likewise, the interface continuity conditions are recovered in an implicit fashion. The rigid motion inside the solid Without loss of generality we assume that M = 1. Then x G (t) = Ωs(t) ρ s x dx, J(t) = Ωs(t) ρ s (r 2 Ir ⊗ r) dx where r(x, t) = xx G (t). The system is subject to a body density force g (usually gravity).

Weak formulation

The basic formulation of this fluid-solid coupling is the following : given initial conditions,

x ′ G (0) = v 0 g , ω u (0) = ω 0 u , u = u 0 , Ω s (0) = Ω 0 s ( 1
)
supplemented with

x G (0) =

Ω 0 s x dx, X s (x, 0) = x, (2)
find t → Ω s (t) and (x, t) → (u(x, t), p(x, t)) solution for t > 0 of

ρ f (u t + (u • ∇)u) -2µ div(D(u)) + ∇p = ρ f g on Ω f (t), (3)
div u = 0 on Ω f (t), ( 4)

u = 0 on ∂Ω, (5)
u = x ′ G + ω u × r on ∂Ω s (t), (6)
x ′′ G (t) = g + ∂Ωs(t)

(Σ n) ds, (7)
J(t)ω ′ u (t) = -ω u (t) × (J(t) ω u (t)) + Ωs(t) ρ s (r × g) dx + ∂Ωs(t) r × (Σ n) ds, (8)
Ω s (t) = X s (t, Ω 0 s ), ( 9)

∂X s ∂t = x ′ G (t) + ω u (t) × r(X s (t), t), (10)
where n denotes the unit outward normal on ∂Ω s (t), and Σ is the fluid stress tensor.

In this formulation the last two equations describe the rigid motion of Ω s (t). In order to give a weak formulation of this problem, let us introduce some function spaces. From now on, u will denote the velocity field on the whole computational domain Ω. We define V = {u ∈ H 1 0 (Ω), div u = 0}, H = {u ∈ L 2 (Ω), div u = 0, u • n = 0 on ∂Ω} and, with the notations of [START_REF] San Martin | Global Weak Solutions for the Two-Dimensional Motion of Several Rigid Bodies in an Incompressible Viscous Fluid[END_REF] extended to the three dimensional case,

K(t) = {u ∈ V, D(u) = 0 in Ω s (t)} = {u ∈ V, ∃(V u , ω u ) ∈ R 3 × R, u = V u + ω u × r in Ω s (t)}.
Next we define a density on the whole domain by setting ρ = ρ s χ Ωs(t) + ρ f χ Ω f (t) , where χ A denotes the characteristic function of set A, which takes value 1 inside A and 0 outside. Let us note Q = Ω×]0, T [. Then the weak formulation is the following [START_REF] Hillairet | Lack of collision between solid bodies in a 2D constant density incompressible viscous flow[END_REF]: given initial conditions

H 0 = χ Ω 0 s , ρ 0 = ρ s H 0 + ρ f (1 -H 0
) and u = u 0 ∈ K(0), find (x, t) → (ρ(x, t), u(x, t), H(x, t)) such that

                                 u ∈ L ∞ (0, T, H) ∩ L 2 (0, T, V), H, ρ ∈ C(0, T ; L q (Ω)) ∀q ≥ 1, u(t) ∈ K(t) for a.e. t ∈]0, T [, with Ω s (t) = {x ∈ Ω, H(x, t) = 1}, ∀ξ ∈ H 1 (Q) ∩ L 2 (0, T ; K(t)), Ω [ρu • ∂ t ξ + (ρ(u • ∇)u -2µD(u)) : D(ξ) + ρg • ξ] dx = d dt Ω ρu • ξ dx, ∀ψ ∈ C 1 (Q), ψ(T ) = 0, T 0 Ω H ∂ψ ∂t + Hu • ∇ψ dxdt + Ω H 0 ψ(0) dx = 0, T 0 Ω ρ ∂ψ ∂t + ρu • ∇ψ dxdt + Ω ρ 0 ψ(0) dx = 0. ( 11
)
Note that we could equivalently have defined ρ = ρ s H + ρ f (1 -H), as ρ 0 is piecewise constant, and transported by the same velocity field than H.

Penalized problem

For η > 0, we consider the following penalized problem:

given (ρ η (0) = ρ 0 , u η (0) = u 0 η , H η (0) = χ Ω 0 s ), to find (ρ η , u η , p η , H η ), with ρ η , H η ∈ L ∞ (]0, T [×Ω), u η ∈ L ∞ (0, T ; H)∩L 2 (0, T ; V), p η ∈ L 2 (Q), ∂ ∂t (ρ η u η ) ∈ L 2 (0, T ; V ′ ) solution in the sense of distributions in Q of ∂(ρ η u η ) ∂t + div(ρ η u η ⊗ u η ) -2µ div(D(u η )) + ∇p η + 1 η ρ η H η (u η -u η,s ) = ρ η g (12)
div u η = 0 (13)

u η,s = 1 M η Ω ρ η u η H η dx + J -1 η Ω ρ η (r η × u η )H η dx × r η ( 14
)
ρ η t + u η .∇ρ η = 0 ( 15
)
H η t + u η,s .∇H η = 0 (16)
We set Ω η s (t) = {x ∈ Ω, H η (x, t) = 1}. In equation [START_REF] Hillairet | Lack of collision between solid bodies in a 2D constant density incompressible viscous flow[END_REF] we divided the first term by M η = Ω ρ η H η dx, which is not constant in time in general. On contrary we have

|Ω η s (t)| = Ω H η dx = |Ω 0 s | since u η,s
is divergence free and H η vanishes on ∂Ω (we assumed no contact of the solid with ∂Ω). The inertia tensor is defined as

J η = Ω ρ η H η (r 2 η I -r η ⊗ r η ) dx = Ω η s (t) ρ η (r 2 η I -r η ⊗ r η ) dx. with r η = x -x Gη = x -Ω ρ η H η x dx. For a ∈ R 3 \ {0}, a T J η a = Ω η s (t) ρ η |r η × a| 2 dx ≥ min(ρ s , ρ f ) Ω η s (t)
|r η × a| 2 dx (see estimate [START_REF] Patankar | A formulation for fast computations of rigid particulate flows[END_REF]). This last quantity being strictly positive for an open nonempty integration set, J η is nonsingular

(we recall that |Ω η s (t)| = |Ω 0 s | > 0).
Before stating our convergence result, a few remarks are in order.

First one may wonder about the well-posedness of the above problem. However it will directly result from the a priori estimates and convergence arguments given in the following that this problem does have at least a weak solution. Indeed, these arguments could easily be used to show the convergence of the solutions to a linearized version -or finite-dimension approximation -of ( 12)-( 16). Next we can observe that in this model we penalize the difference between u η and the projection of u η onto velocity fields rigid in the solid domain, namely u η,s (see lemma 3.1 below).

The density is transported with the original velocity field so that estimates on the Navier-Stokes equations are easier to obtain. The characteristic function is transported by the rigid velocity so that the shape of Ω η s (t) remains undeformed (this is exactly the Eulerian counterpart of (9-10)). As observed in [START_REF] Coquerelle | A vortex level-set method for the two-way coupling of an incompressible fluid with colliding rigid bodies[END_REF] this has a practical importance (in particular it means that the rigid solid can be recovered exactly through simple algebra from its initial shape). As far as numerical analysis is concerned, it also provides "for free" regularity properties on the computed rigid body, as soon as the initial body is smooth. The price to pay is that the level sets of ρ η and H η do not coincide, i.e. in general we do not have ρ η H η = ρ s as in the non penalized formulation. Note also that in principle we should prescribe a boundary value for H η on ∂Ω when u η,s is inward. Since our analysis is restricted to times when the solid body does not approach the boundary of the computational box, we can take this boundary value to be zero, which amounts to solve (16) on R n and take its restriction to Ω.

Let us now comment on the differences between the present method and the methods in [START_REF] San Martin | Global Weak Solutions for the Two-Dimensional Motion of Several Rigid Bodies in an Incompressible Viscous Fluid[END_REF][START_REF] Patankar | A formulation for fast computations of rigid particulate flows[END_REF][START_REF] Sharma | A fast computation technique for the direct numerical simulation of rigid particulate flows[END_REF]. The method used in [START_REF] San Martin | Global Weak Solutions for the Two-Dimensional Motion of Several Rigid Bodies in an Incompressible Viscous Fluid[END_REF] to prove the existence of the fluid-body problem uses a penalization acting on the deformation. It does not require to compute the rigid body velocity u η,s . However it does not guarantees a strictly rigid transport of the body, which is clearly a drawback form a practical point of view. The projection method in [START_REF] Patankar | A formulation for fast computations of rigid particulate flows[END_REF][START_REF] Sharma | A fast computation technique for the direct numerical simulation of rigid particulate flows[END_REF] can be seen as a particular discretization of the present penalization method, with an explicit time-discretization of the penalization term. In practice this result in a limitation of the penalization coefficient which has some We will come back to this point in our numerical illustrations in section 4. Numerical comparisons for a benchmark challenging 3D case are also given in [START_REF] Coquerelle | A vortex level-set method for the two-way coupling of an incompressible fluid with colliding rigid bodies[END_REF].

In the following sections we will prove the convergence of at least a subsequence of (ρ η , u η , p η , H η )

to the weak solution defined above. Next section starts with some a priori estimates which will provide weak convergence of subsequences. In section 3.3 we will have to use more sophisticated tools adapted from [START_REF] San Martin | Global Weak Solutions for the Two-Dimensional Motion of Several Rigid Bodies in an Incompressible Viscous Fluid[END_REF] to get some strong convergence in u η which will allow us to pass to the limit in nonlinear terms of (P η ). More precisely we prove the following result.

Theorem 2.1. Under the regularity assumptions of section 2, let (ρ η , u η , p η , H η ) a solution of (P η ). Then there exists a subsequence of (ρ η , u η , H η ) and functions (ρ, u, H) such that ρ η → ρ, H η → H strongly in C(0, T ; L q (Ω)) for all q ≥ 1,

u η → u strongly in L 2 (Q) and weakly in L 2 (0, T ; H 1 0 (Ω)) ∩ L ∞ (0, T ; L 2 (Ω))
and such that (ρ, u, H), is a solution of [START_REF] Glowinski | A Fictitious Domain Approach to the Direct Numerical Simulation of Incompressible Viscous Flow past Moving Rigid Bodies: Application to Particulate Flow[END_REF].

Before proceeding to the proof, let us point out a few remarks. For a sake of simplicity in the notations we have stated our penalization method and theorem for a single rigid body. It will be apparent from the proof below that it readily extends to the case of several bodies. Furthermore, the time to which the convergence result is restricted, is essentially the time for which contact of the rigid body do not touch the boundary of Ω (in the case of several bodies it would be the time on which we can ensure that contact between bodies do not happen). As a result if we consider periodic boundary conditions and a single body convergence holds for all times.

3 Proof of theorem 2.1

The following lemma states that u η,s , as defined in (P η ), is the projection of u η onto velocity fields which are rigid on Ω η s (t).

Lemma 3.1. Let ξ be a rigid velocity field, i.e. such that ξ(x) = V ξ + ω ξ × r(x) for some constant vectors V ξ ∈ R 3 and ω ξ ∈ R 3 . Then if u η,s is defined by ( 14) there holds

Ω ρ η H η (u η -u η,s ) • ξ dx = 0. ( 17
)
Moreover, the result holds if ξ is a time dependent velocity field rigid in Ω η s (t) at time t.

Proof. Let the mean translation and angular velocities be defined as

V u = 1 M η Ω ρ η H η u η dx ω u = J -1 η Ω ρ η H η (r η × u η ) dx then Ω ρ η H η (u η -u η,s ) • ξ dx = Ω ρ η H η [u η -(V u + ω u × r η )] • [V ξ + ω ξ × r η ] dx = V ξ • Ω ρ η H η u η dx + ω ξ • Ω ρ η H η (r η × u η ) dx -V u • V ξ Ω ρ η H η dx -V u • ω ξ × Ω ρ η H η r η dx -V ξ • ω u × Ω ρ η H η r η dx - Ω ρ η H η (ω u × r η ) • (ω ξ × r η ) dx = V ξ • (M η V u ) + ω ξ • (J η ω u ) -V u • (M η V ξ ) -V u • ω ξ × Ω ρ η H η r η dx -V ξ • ω u × Ω ρ η H η r η dx - Ω ρ η H η (ω u × r η ) • (ω ξ × r η ) dx. As (ω u × r η ) • (ω ξ × r η ) = (ω ξ • ω u )r 2 η -(r η • ω ξ )(r η • ω u ), we have Ω ρ η H η (ω u × r η ) • (ω ξ × r η ) dx = ω ξ • (J η ω u ).
Finally, by definition of r η , Ω ρ η H η r η dx = 0, and we get

Ω ρ η H η (u η -u η,s ) • ξ dx = ω ξ • (J η ω u ) -ω ξ • (J η ω u ) = 0.

Estimates for transport and Navier-Stokes equations

In all the sequel, C denotes a positive constant. At this stage, we consider a given time interval [0, T ]. The value to which T must be restricted will be given later in this section.

Standard estimates for transport equations ( 15) and [START_REF] Lions | Quelques méthodes de résolutions des problèmes aux limites non linéaires[END_REF] show that ρ η and H η are bounded in L ∞ (0, T ; L ∞ (Ω)). More precisely, for all time t ∈ [0, T ],

ρ min := min(ρ s , ρ f ) ≤ ρ η (x, t) ≤ max(ρ s , ρ f ) H η (x, t) ∈ {0, 1} a.e. x ∈ Ω. ( 18
)
Thus, up to extracting a subsequence, we can assume that

ρ η ⇀ ρ in L ∞ (0, T, L ∞ (Ω)) weak*, (19)
and

H η ⇀ H in L ∞ (0, T, L ∞ (Ω)) weak*, (20)
where H and ρ satisfy the bounds [START_REF] Patankar | A formulation for fast computations of rigid particulate flows[END_REF]. We set Ω s (t) = {x ∈ Ω, H(x, t) = 1}. Concerning the Navier-Stokes equations, using u η as a test function in the weak formulation of (12) and using [START_REF] Janela | A Penalty Method for the simulation of fluid-body rigid interaction[END_REF] we classically obtain:

1 2 Ω ∂(ρ η |u η | 2 ) ∂t dx + 2µ Ω |D(u η )| 2 dx + 1 η Ω ρ η H η (u η -u η,s ) • u η dx = Ω ρ η g • u η dx.
From Lemma 3.1, since u η,s is a rigid velocity field, we get

Ω ρ η H η (u η -u η,s ) • u η dx = Ω ρ η H η (u η -u η,s ) 2 dx.
Collecting terms we get, since from ( 18)

H η = H η , 1 2
d dt √ ρ η u η 2 L 2 (Ω) +µ D(u η ) 2 L 2 (Ω) + 1 η √ ρ η H η (u η -u η,s ) 2 L 2 (Ω) ≤ √ ρ η u η L 2 (Ω) g L ∞ (Q) ρ η L 2 (Q)
which upon time integration on [0, T ] gives

√ ρ η (t)u η (t) 2 L 2 (Ω) + 2µ D(u η ) 2 L 2 (Q) + 2 η √ ρ η H η (u η -u η,s ) 2 L 2 (Q) ≤ √ ρ η0 u η0 2 L 2 (Ω) + C T 0 √ ρ η (s)u η (s) L 2 (Ω) ds.
Applying Gronwall Lemma, Poincaré inequality and bounds from [START_REF] Patankar | A formulation for fast computations of rigid particulate flows[END_REF] gives the following estimates :

u η bounded in L 2 (0, T, H 1 0 (Ω)), (21)
√ ρ η u η and u η bounded in L ∞ (0, T, L 2 (Ω)), ( 22
) 1 √ η √ ρ η H η (u η -u η,s ) and 1 √ η H η (u η -u η,s ) bounded in L 2 (0, T, L 2 (Ω)). ( 23
)
Thus we can extract subsequences from ρ η , u η and H η , still denoted by ρ η , u η and H η , such that

u η ⇀ u in L 2 (0, T, H 1 0 (Ω)) weak, ( 24
) √ ρ η u η ⇀ χ and u η ⇀ u in L ∞ (0, T, L 2 (Ω)) weak*, (25)
√ ρ η H η u η - √ ρ η H η u η,s → 0 and H η u η -H η u η,s → 0 in L 2 (0, T, L 2 (Ω)) strong. ( 26
)
The identification of χ with √ ρu results from strong convergence results proved by Lions and DiPerna on transport equations. [START_REF] Di Perna | Ordinary differential equations, transport theory and Sobolev spaces[END_REF] theorem II.4, (24) and incompressibility imply

ρ η → ρ in C(0, T, L q (Ω)) strong ∀q ∈ [1, +∞[ (27)
with ρ solution of

ρ t + u • ∇ρ = 0 on Ω×]0, T [, ρ = ρ 0 on Ω × {0}.
From this strong convergence we can pass to the limit in the product √ ρ η u η : given v ∈ L q (0, T ; L r (Ω)) with q > 2 and r > 6 5 , we write

T 0 Ω ( √ ρ η u η - √ ρu)v dxdt = T 0 Ω (u η -u) √ ρv dxdt + T 0 Ω ( √ ρ η - √ ρ)u η v dxdt.
From the injection of H 1 into L 6 in dimension less or equal to 3 the first integral converges toward 0. For the second integral we use the strong convergence (following (27)) of √ ρ η in L s for a s such that u η v is in L s ′ where s ′ is the conjugate exponent of s. Thus we have

√ ρ η u η ⇀ √ ρu in L q (0, T, L r (Ω))
weak, for all q < 2, r < 6.

(28)

Setting T and passing to the limit in the rigid velocity

This rigid velocity is defined by

u η,s (x, t) = u η,G (t) + ω η (t) × r η (x, t), with u η,G (t) = 1 M η Ω ρ η u η H η dx and ω η (t) = J -1 η Ω ρ η (r η × u η )H η dx.
First we note that M η = Ω ρ η H η dx is bounded from below independently of η since ρ η is bounded from below and Ω H η dx = |Ω 0 s | > 0 does not depend on η or t. From the bounds on ρ η , H η and

u η it is straightforward to show that u η,G (t) bounded in L ∞ (0, T ).
Likewise, from the definition of J η we observe that for a ∈ R 3 \ {0}

a T J η a ≥ min(ρ s , ρ f ) Ωs(t) |r η × a| 2 dx > 0.
Moreover, the initial solid is regular and transported by a rigid velocity. We thus know that there is a ball of radius R > 0 centered on the center of gravity x Gη included into Ω η s (t). Then the above estimates implies

a T J η a ≥ min(ρ s , ρ f ) B(x Gη ,R) |r η × a| 2 dx = min(ρ s , ρ f ) B(0,R) |x × a| 2 dx = C(R)|a| 2 with C(R) = 2R 5 π 15 > 0. Taking a = J -1 2 η b (J η is symmetric) we get for all b ∈ R 3 \ {0}, b T J -1 η b = |J -1 2 η b| 2 ≤ 1 C(R) |b| 2 ,
which proves that each coefficient of J -1 η is bounded independently of η and t. From the bounds on u η , H η and ρ η this implies that

ω η (t) is bounded in L ∞ (0, T ).
In particular this implies that the solid velocity u η,s is bounded in L ∞ by some constant M independent of η and time. We can now define the maximum time for which the convergence result will be proved. If we denote by d 0 the initial distance between solid and the boundary ∂Ω then choosing for instance T = d 0 /2M ensures that the body will not touch the boundary for

t ∈ [0, T ].
In all the sequel we will assume this value of T .

From the above estimates we can ensure that there exists u G (t) and ω(t) in L ∞ (0, T ) such that, up to the extraction of subsequences,

u η,s ⇀ u s := u G + ω × r in L ∞ (0, T, L ∞ (Ω)) weak*.
Now we point out that taking the gradient of the rigid velocity field u η,s gives

∇u η,s =   0 -ω 3 η ω 2 η ω 3 η 0 -ω 1 η -ω 2 η ω 1 η 0  
so that the ∇u η,s (and all subsequent space derivatives) is also bounded in L ∞ (0, T ; L ∞ (Ω)). In particular

u η,s ⇀ u s in L 2 (0, T, W 1,∞ (Ω)) weak*. ( 29
)
We now wish to prove that u s (or equivalently u G and ω) has a similar structure as u η,s , that is, to pass to the limit in the expression defining u η,s . Using (29), the already mentioned compactness results of [START_REF] Di Perna | Ordinary differential equations, transport theory and Sobolev spaces[END_REF], applied to the transport equation on H η now gives

H η → H a.e. and in C(0, T, L p (Ω)) strong ∀p ∈ [1, +∞[ (30)
with H verifying

H t + u s .∇H = 0 on R n × (0, T ) H = H 0 on R n × {0}.
Note that this Cauchy problem has been set in R n because u s does not vanish on ∂Ω, but H vanishes outside Ω. However we can prove by passing to the limit in (26) that Hu = Hu s . Thus div(uH) = div(u s H) and H verifies a transport equation with velocity field u on Ω (note that no boundary conditions are needed on ∂Ω since u is zero on the boundary).

This convergence gives us the strong convergence of r η in C(0, T, L p (Ω)), ∀p ≥ 1, and from ( 27),( 30) and ( 24) we can easily pass to the limit in the expression of u η,G and ω η to get

u G (t) = Ω ρuH dx Ω ρH dx and ω(t) = Ω ρ(r 2 I -r ⊗ r)H dx -1 Ω ρ(r × u)H dx.

Strong convergence of u η

The remaining part of the proof is more technical since we aim to prove the strong convergence of at least a subsequence of u η in order to be able to pass to the limit in the inertial term of Navier-Stokes equations. Classically, this is obtained thanks to a Fourier transform in time which provides an estimate on some fractional time derivative of u η which brings compactness [START_REF] Lions | Quelques méthodes de résolutions des problèmes aux limites non linéaires[END_REF][START_REF] Temam | Navier-Stokes equations and numerical analysis[END_REF].

Here these technics can not be used since the solid is moving. We instead rely on tools developed in [START_REF] San Martin | Global Weak Solutions for the Two-Dimensional Motion of Several Rigid Bodies in an Incompressible Viscous Fluid[END_REF].

Thereafter we will use, for σ > 0 and r ∈ [1/2, 1], the following notations

• Ω s,σ (t) = {x ∈ Ω, d(x, Ω s (t)) < σ}, • V 0 = {v ∈ L 2 (Ω), div v = 0, v • n = 0 on ∂Ω}, • V r = {v ∈ H r (Ω), div v = 0, v = 0 on ∂Ω}, • K r σ (t) = {v(t) ∈ V r , D(v(t)) = 0 in D ′ (Ω s,σ (t))} (which is a closed subset of V r ),
• P r σ (t) the orthogonal projection in the sense of the H r norm of V r on K r σ (t).

To prove the strong convergence of a subsequence of u η in L 2 (Q) we write

T 0 Ω |u η -u| 2 dxdt ≤ 1 ρ min T 0 Ω |ρ(u 2 η -u 2 )| dxdt + T 0 Ω |2ρu • (u -u η )| dxdt .
From (25) the second integral on the right side converges to 0, thus

T 0 Ω |u η -u| 2 dxdt ≤ 1 ρ min T 0 Ω |ρ η u 2 η -ρu 2 | dxdt + T 0 Ω |(ρ η -ρ)u 2 η | dxdt + l η ,
where l η → 0 when η → 0. Moreover, by ( 21) and ( 27) the second integral on the right hand side converges to 0, and

T 0 Ω |u η -u| 2 dxdt ≤ 1 ρ min T 0 Ω |ρ η u η • P r σ (u η ) -ρu • P r σ (u)| dx dt + T 0 Ω |ρ η u η • (u η -P r σ (u η ))| dx dt + T 0 Ω |ρu • (P r σ (u) -u)| dx dt + l η ≤ 1 ρ min ( ρ η u η • P r σ (u η ) -ρu • P r σ (u) L 1 (Q) + ρ η L ∞ (Q) u η L 2 (Q) P r σ (u η ) -u η L 2 (Q) + ρ L ∞ (Q) u L 2 (Q) P r σ (u) -u L 2 (Q) + l η
Finally, as ρ η is bounded in L ∞ (0, T, L ∞ (Ω)) and using [START_REF] Sarthou | Eulerian-Lagrangian grid coupling and penalty methods for the simulation of multiphase flows interacting with complex objects[END_REF], we get

T 0 Ω |u η -u| 2 dxdt ≤ 1 ρ min ( ρ η u η • P r σ (u η ) -ρu • P r σ (u) L 1 (Q) + C P r σ (u η ) -u η L 2 (Q) + C P r σ (u) -u L 2 (Q) + l η (31)
This decomposition shows that the sought convergence essentially amounts to prove that (up to the extraction of subsequences)

lim σ→0 P r σ (u) -u L 2 (Q) = 0 (32) lim σ→0 lim η→0 P r σ (u η ) -u η L 2 (Q) = 0 (33) lim σ→0 lim η→0 ρ η u η • P r σ (u η ) -ρu • P r σ (u) L 1 (Q) = 0 (34)
To prove (32-34) we will make use of some Lemma that we state now.

Lemma 3.2. Let (f n ) be a sequence of functions bounded in L p (0, T ) for some p > 2 and converging to 0 almost everywhere on [0, T ]. Then (f n ) converges strongly to 0 in L 2 (0, T ).

Proof. Let ε > 0. From Egorov theorem ([5], p.75), the almost everywhere convergence implies that there exists

A ε ⊂ [0, T ] such that |[0, T ]\A ε | < ε f n → 0 uniformly on A ε which means ∃N ∈ N, ∀n ≥ N, ∀t ∈ A ε , |f n (t)| 2 < ε. Therefore Aε (f n (t)) 2 dt ≤ ε|A ε | ≤ εT.
Since, from the bound in the L p norm

[0,T ]\Aε (f n (t)) 2 dt ≤ [0,T ]\Aε 1 q dt 2/q [0,T ] (f n (t)) p dt 2/p ≤ ε 2/q C with 1 q + 1 p = 1 2 , we get T 0 (f n (t)) 2 dt ≤ εT + ε 2/q C
which proves the L 2 convergence.

         -∆w(t) + ∇p(t) = 0 on Ω\Ω s,σ (t), div w(t) = 0 on Ω\Ω s,σ (t), w(t) = g(t)
on ∂Ω s,σ (t), w(t) = 0 on ∂Ω.

Then there exists σ 0 > 0 and C > 0, such that for all σ < σ 0 , we have the following estimate:

w(t) 2 L 2 (Ω\Ωs,σ(t)) ≤ C u(t) L 2 (Ωs,σ(t)) ∇u(t) L 2 (Ωs,σ(t))
The proof of this result is postponed to the appendix.

Lemma 3.4. The limits H, u and u s defined in (30), ( 24) and (29) verify

Hu = Hu s . (35)
Proof. First we claim that H η u η ⇀ Hu in L q (0, T, L r (Ω)) weak, with q < 2 and r < 6 (36)

Indeed, let us introduce v ∈ L q (0, T, L r (Ω)) with q > 2 and r > 6 5 . We have

T 0 Ω (H η u η -Hu) • v dxdt = T 0 Ω H(u η -u) • v dxdt + T 0 Ω (H η -H)u η • v dxdt
From the injection of H 1 into L 6 in dimension less or equal to 3 and (24) we get

u η ⇀ u in L 2 (0, T, L 6 (Ω)) weak,
and

, since H is bounded in L ∞ (Q), lim η→0 T 0 Ω H(u η -u) • v dx dt = 0.
Moreover from (30) we easily get

lim η→0 T 0 Ω (H η -H)u η • v dx dt = 0.
We next show that

H η u η,s ⇀ Hu s in L p (0, T, L p (Ω)) weak, ∀p ∈ [1, +∞[. (37)
Let v ∈ L p (0, T, L p (Ω)) with p > 1. We write

T 0 Ω (H η u η,s -Hu s ) • v dxdt = T 0 Ω H(u η,s -u s ) • v dxdt + T 0 Ω (H η -H)u η,s • v dx dt.
As H is bounded in L ∞ (Q), with (29) we get:

lim η→0 T 0 Ω H(u η,s -u s ) • v dx dt = 0
In addition, from (30) we get

lim η→0 T 0 Ω (H η -H)u η,s • v dx dt = 0.
By (36) and (37) we thus deduce H η u η -H η u η,s ⇀ Hu -Hu s in L q (0, T, L r (Ω)) weak, with q < 2 and r < 6 (38)

Finally we recall that

H η u η -H η u η,s → 0 in L 2 (Q) strong, (39)
and the desired result follows. The result is finally obtained by identifying the limits in (38) et (39).

Lemma 3.5.

∀σ > 0, ∃η 0 > 0, ∀η < η 0 , ∀t ∈ [0, T ], Ω η s (t) ⊂ Ω s,σ (t) and Ω s (t) ⊂ Ω η s,σ (t). ( 40
)
Proof. From (30) with p = 1 we have

∀ε > 0, ∃η 0 > 0, ∀η < η 0 , ∀t ∈ [0, T ], Ω |H η (x, t) -H(x, t)| dx < ε, (41)
which means

∀ε > 0, ∃η 0 > 0, ∀η < η 0 , ∀t ∈ [0, T ], |Ω η s (t) \ Ω s (t)| + |Ω s (t) \ Ω η s (t)| < ε. ( 42
)
By contradiction we suppose that we can find σ 0 > 0 such as ∀η 0 > 0, there exists η < η 0 and t ∈ [0, T ] for which at least one of the inclusions of (40) is false. Assume the first inclusion is false.This means that we can find x η (t) ∈ Ω η s (t) such as d(x η (t), Ω s (t)) > σ 0 . Ω η s (t) is a rigid deformation of Ω η s (0) so its boundary is C 2 . Thus, there exists a sufficiently small ρ independent of η, such as for each point of Ω η s (t) there exists a ball of radius ρ > 0 containing this point and included in Ω η s (t). Then there exists also a ball of radius ρ := min(ρ, σ 0 /3) containing the point and included in Ω η s (t). This latter ball is included in Ω η s (t) \ Ω s (t). Indeed it contains a point at distance more than σ 0 from Ω s (t) and its diameter is less than 2σ 0 /3. We thus obtained that

∃σ 0 , ∀η 0 > 0, ∃η > 0, ∃t ∈ [0, T ], |Ω η s (t) \ Ω s (t)| > π ρ2
with ρ independent of η and t. This contradicts (42). A similar argument shows that the second inclusion cannot hold either.

A key point in the convergence is the following result.

Lemma 3.6.

lim σ→0 P r σ (u) -u L 2 (0,T,V r ) = 0, ∀r ∈ [1/2, 1[ (43)
The proof of this result is postponed to the appendix.

The next lemma, which is also proved in the appendix, is essentially a rephrasing of the previous one with u η instead of u. The difference is that we do not have anymore u ηu η,s = 0 ins Ω η s (t), but we do have an estimate on it, from [START_REF] Temam | Navier-Stokes equations and numerical analysis[END_REF], which allows to pass to the limit as η goes to 0.

Lemma 3.7.

lim σ→0 lim η→0 P r σ (u η ) -u η L 2 (0,T,V r ) = 0, ∀r ∈]1/2, 1[ (44)
Lemma 3.8.

lim σ→0 lim η→0 ρ η u η • P r σ (u η ) -ρu • P r σ (u) L 1 ([0,T ]×Ω) = 0, ∀r ∈]1/2, 1[. ( 45
)
Proof. Let r ∈]1/2, 1[ and σ > 0. From Lemma 3.5 there exists η 0 > 0 such that ∀η < η 0 ,

Ω η s (t) ⊂ Ω s,σ/3 (t), ∀t ∈ [0, T ].
Let η < η 0 . Arguing as in [START_REF] San Martin | Global Weak Solutions for the Two-Dimensional Motion of Several Rigid Bodies in an Incompressible Viscous Fluid[END_REF] we split [0, T ] in N T subintervals

I k = [(k -1)τ, kτ ], τ = T /N T , k = 1, .., N T . We choose N T large enough (depending on σ) such that Ω s,σ/2 (kτ ) ⊂ Ω s,σ (t) and Ω s,σ/3 (t) ⊂ Ω s,σ/2 (kτ ), ∀t ∈ I k , ∀k = 1, .., N T . (46)
This is possible because Ω s,σ (t) is moving with a rigid velocity field, with L 2 regularity in time :

the flow X(t, x) generated by this velocity field is thus continuous in time. Let x ∈ Ω s,σ/2 (kτ ).

There exits y ∈ Ω s (kτ ) such that |y -x| ≤ σ/2. But y = X(kτ, z) for some z ∈ Ω s (0) and therefore

|y -X(t, z)| ≤ σ/2 for t ∈ I k if τ is small enough. Therefore |x -X(t, z)| ≤ σ.
This proves the first inclusion in (46). The second inclusion is proved by a similar continuity argument. On each subinterval I k , k = 1, .., N T , we consider the momentum equation

ρ η ∂u η ∂t + ρ η (u η • ∇)u η -2µ div(D(u η )) + ∇p η + 1 η ρ η H η (u η -u η,s ) -ρ η g = 0.
Let us consider a test function ξ vanishing outside I k and such that ξ(., t) ∈ K 1 σ/2 (kτ ) for t ∈ I k .

Since ξ(., t) is rigid on Ω s,σ/2 (kτ ) ⊃ Ω s,σ/3 (t) ⊃ Ω η s (t), Lemma 3.1 yields

I k Ω [ρ η u η • ξ t + (ρ η u η ⊗ u η -2µD(u η )) : D(ξ) + ρ η g • ξ] dxdt = 0.
From bounds given by ( 18), ( 21) and ( 22) we derive the following estimates:

I k Ω D(u η ) : D(ξ) dxdt ≤ D(u η ) L 2 (I k ,L 2 (Ω)) D(ξ) L 2 (I k ,L 2 (Ω)) ≤ C ξ L 2 (I k ,H 1 0 (Ω)) ≤ C ξ L 2 (I k ,K 1 σ/2 (kτ )) ≤ C ξ L 4 (I k ,K 1 σ/2 (kτ )) , I k Ω (ρ η u η ⊗ u η ) : D(ξ) dxdt ≤ I k ρ η L ∞ (Ω) u η ⊗ u η L 2 (Ω) D(ξ) L 2 (Ω) dt ≤ ρ η L ∞ (I k ,L ∞ (Ω)) I k u η 2 L 4 (Ω) ξ H 1 0 (Ω) dt ≤ C I k u η 1 2 L 2 (Ω) ∇u η 3 2 L 2 (Ω) ξ H 1 0 (Ω) dt ≤ C u η 1 2 L ∞ (I k ,L 2 (Ω)) I k ∇u η 3 2 L 2 (Ω) ξ H 1 0 (Ω) dt ≤ C ∇u η 3 2 L 2 (I k ,L 2 (Ω)) ξ L 4 (I k ,H 1 0 (Ω)) ≤ C ξ L 4 (I k ,H 1 0 (Ω)) ≤ C ξ L 4 (I k ,K 1 σ/2 (kτ )) ,
and

I k Ω ρ η g • ξ dxdt ≤ ρ η L ∞ (I k ,L ∞ (Ω)) ξ L 4 (I k ,K 1 σ/2 (kτ )) ≤ C ξ L 4 (I k ,K 1 σ/2 (kτ )) .
Collecting terms we get

I k Ω ρ η u η • ξ t dxdt ≤ C ξ L 4 (I k ,K 1 σ/2 (kτ )) .
As ξ(., t) ∈ K 1 σ/2 (kτ ), ξ t (., t) ∈ K 1 σ/2 (kτ ) ⊂ K 0 σ/2 (kτ ) and we have

| ρ η u η , ξ t | = | ρ η u η , P 0 σ/2 (kτ )ξ t | = | P 0 σ/2 (kτ )(ρ η u η ), ξ t | = | d dt P 0 σ/2 (kτ )(ρ η u η ), ξ |.
Therefore

I k Ω d dt P 0 σ/2 (kτ )(ρ η u η ) • ξ dxdt ≤ C ξ L 4 (I k ,K 1 σ/2 (kτ )) ,
which means that

d dt P 0 σ/2 (kτ )(ρ η u η ) bounded in L 4 3 (I k , (K 1 σ/2 (kτ )) * ). ( 47
) Moreover ρ η u η is bounded in L 2 (I k , L 2 (Ω)), P 0 σ/2 (kτ )(ρ η u η ) is bounded in L 2 (I k , K 0 σ/2 (kτ )). ( 48
)
Since K 0 σ/2 (kτ ) ⊂ (K r σ/2 (kτ )) * compactly for r > 0, and (K r σ/2 (kτ )) * ⊂ (K 1 σ/2 (kτ )) * continuously for r < 1, by the Aubin-Simon Lemma (see e.g. [START_REF] Boyer | Eléments d'analyse pour l'étude de quelques modèles d'écoulements de fluides visqueux incompressibles[END_REF], p. 98) with ( 47) and (48), we obtain the relative compactness of the sequence

P 0 σ/2 (kτ )(ρ η u η ) in L 2 (I k , (K r σ/2 (kτ )) * ) for all r ∈]1/2, 1[.
From (24) we deduce

lim η→0 P 0 σ/2 (kτ )(ρ η u η ) = P 0 σ/2 (kτ )ρu in L 2 (I k , (K r σ/2 (kτ )) * ) strong, ∀r ∈]1/2, 1[. ( 49
)
Since from (46), we have

P 0 σ/2 (kτ )P r σ (t) = P r σ (t) ∀t ∈ I k ∀r ∈]1/2, 1[, (50)
and we can write

I k ρ η u η , P r σ (t)(u η ) L 2 (Ω) dt = I k ρ η u η , P 0 σ/2 (kτ )P r σ (t)u η L 2 (Ω) dt = I k P 0 σ/2 (kτ )(ρ η u η ), P r σ (t)u η L 2 (Ω) dt = I k P 0 σ/2 (kτ )(ρ η u η ), P r σ (t)u η (K r σ/2 ) * ,K r σ/2 dt. The sequence (u η ) is bounded in L 2 (0, T, V r ) for all r ∈]1/2, 1[, therefore (P r σ (t)u η ) is bounded in L 2 (0, T, K r σ/2 ) for all r ∈]1/2, 1[.
Therefore there exists a subsequence of P r σ (t)u η still denoted P r σ (t)u η such that

P r σ (t)u η ⇀ P r σ (t)u in L 2 (0, T, K r σ/2 ) weak. ( 51
)
Passing to the limit in η yields

lim η→0 I k ρ η u η , P r σ (t)(u η ) L 2 (Ω) dt = I k P 0 σ/2 (kτ )ρu, P r σ (t)u L 2 (Ω) dt = I k ρu, P 0 σ/2 (kτ )P r σ (t)u L 2 (Ω) dt = I k ρu, P r σ (t)u L 2 (Ω) dt.
Summing over k = 1, .., N T , we finally obtain

lim η→0 ρ η u η • P r σ (u η ) -ρu • P r σ (u) L 1 (Q) = 0, which implies lim σ→0 lim η→0 ρ η u η • P r σ (u η ) -ρu • P r σ (u) L 1 (Q) = 0.
We can now conclude to the strong convergence of u η .

Let ε > 0. From Lemma 3.6,

∃σ 0 > 0, ∀σ < σ 0 , P r σ (u) -u L 2 (Q) < ε.
From Lemma 3.7,

∃σ 0 > 0, ∀σ < σ 0 , ∃η 0 > 0, ∀η < η 0 , P r σ (u η ) -u η L 2 (Q) < ε,
and by Lemma 3.8,

∃σ 0 > 0, ∀σ < σ 0 , ∃η 0 > 0, ∀η < η 0 , ρ η u η P r σ (u η ) -ρuP r σ (u) L 1 (Q) < ε.
We therefore get from (31) (up to the extraction of a subsequence)

∃η 0 > 0, ∀η < η 0 , T 0 Ω |u η -u| 2 dxdt < Cε
which means that (still up to a subsequence)

u η → u in L 2 (Q) strong. ( 52
)
Classically, we also obtain from ( 19)

ρ η u η ⇀ ρu in L 2 (Q) weak. ( 53
)

Passing to the limit

Let us now prove that as η goes to zero, a subsequence of (u η , ρ η ) converges toward (u, ρ) solution of the weak formulation [START_REF] Glowinski | A Fictitious Domain Approach to the Direct Numerical Simulation of Incompressible Viscous Flow past Moving Rigid Bodies: Application to Particulate Flow[END_REF]. Indeed : We have proved that

ρ η ⇀ ρ in L ∞ (0, T, L ∞ (Ω)) weak * . Therefore ρ ∈ L ∞ (Q).
We have proved that

u η ⇀ u in L 2 (0, T, V ) weak, that √ ρ η u η is bounded in L ∞ (0, T, L 2 (Ω)) and
ρ η bounded from above and below in L ∞ (0, T, L ∞ (Ω)). This implies that u η bounded in L ∞ (0, T, L 2 (Ω)), thus its weak limit belongs to

L ∞ (0, T, L 2 (Ω)) ∩ L 2 (0, T, V ): u ∈ L ∞ (0, T, L 2 (Ω)) ∩ L 2 (0, T, V ).
From Lemma 3.4, we have Hu = Hu s = H(u G + ω × r) where H is the characteristic function of Ω s (t). Thus

u(t) ∈ K(t).
Using compactness results of DiPerna-Lions we already obtained that ρ and H are solutions of transport equations with u and u s as velocities. For H this means that for all ψ ∈ C 1 (Q) with

ψ(T ) = 0, T 0 Ω H ∂ψ ∂t + Hu s • ∇ψ dxdt + Ω H 0 ψ(0) dx = 0.
As from Lemma 3.4, Hu s = Hu, H is also solution of

T 0 Ω H ∂ψ ∂t + Hu • ∇ψ dxdt + Ω H 0 ψ(0) dx = 0.
In other terms H, like ρ satisfies a transport equation with velocity u.

Let us finally check that u satisfies the momentum equation.

Let 12) and ( 15) we get

σ > 0. If ξ σ ∈ H 1 (Q) ∩ L 2 (0, T ; K 1 σ (t)), from (
Ω ∂(ρ η u η ) ∂t + div(ρ η u η ⊗ u η ) -2µ div(D(u η )) + ∇p η + 1 η H η ρ η (u η -u η,s ) -ρ η g • ξ σ dx = 0.
From Lemma 3.5, there exists η 0 such that η < η 0 implies:

Ω H η ρ η (u η -u η,s ) • ξ σ dx = 0.
By integration by parts

Ω -2µ div(D(u η )) • ξ σ dx = Ω 2µD(u η ) : D(ξ σ ) dx, Ω div(ρ η u η ⊗ u η ) • ξ σ dx = Ω -(ρ η u η ⊗ u η ) : D(ξ σ ) dx, Ω ∂(ρ η u η ) ∂t • ξ σ dx = d dt Ω ρ η u η • ξ σ dx - Ω ρ η u η • ∂ξ σ ∂t dx.
As a result

Ω ρ η u η • ∂ξ σ ∂t + (ρ η u η ⊗ u η -2µD(u η )) : D(ξ σ ) + ρ η g • ξ σ dx = d dt Ω ρ η u η • ξ σ dx.
We have already established that

u η ⇀ u in L 2 (0, T, H 1 0 (Ω)) weak, u η → u in L 2 (0, T, L 2 (Ω)) strong, ρ η u η ⇀ ρu in L 2 (0, T, L 2 (Ω)
) weak, and ρ η → ρ in L 2 (0, T, L 2 (Ω)) strong. Letting η goes to zero, we thus obtain

Ω ρu • ∂ξ σ ∂t + (ρu ⊗ u -2µD(u)) : D(ξ σ ) + ρg • ξ σ dx = d dt Ω ρu • ξ σ dx.
which corresponds to the weak formulation [START_REF] Glowinski | A Fictitious Domain Approach to the Direct Numerical Simulation of Incompressible Viscous Flow past Moving Rigid Bodies: Application to Particulate Flow[END_REF]. This holds for any

ξ σ ∈ H 1 (Q) ∩ L 2 (0, T ; K 1 σ (t)),
for arbitrary σ > 0, and, since the time interval has been chosen to guarantee that there is no contact with the boundary, by Proposition 4.3 of [START_REF] San Martin | Global Weak Solutions for the Two-Dimensional Motion of Several Rigid Bodies in an Incompressible Viscous Fluid[END_REF], for all ξ ∈ H 1 (Q) ∩ L 2 (0, T ; K(t)). This ends the proof of theorem 2.1.

Numerical simulations

We give here a few numerical illustrations of the penalization method in 2D. We only sketch the numerical discretization and we refer to [START_REF] Bost | Méthodes Level-Set et pénalisation pour le calcul d'interactions fluide-structure[END_REF] for a more detailed description and further numerical results.

Numerical 2D and 3D validations against experimental and other numerical results for a discretization of the present penalization method in a vorticity formulation of the Navier-Stokes equations are given in [START_REF] Coquerelle | A vortex level-set method for the two-way coupling of an incompressible fluid with colliding rigid bodies[END_REF].

We choose a time-step ∆t and denote by a superscript n discretization of all quantities at time t n = n∆t. For each time integration, we split the penalization model ( 12)-( 16) as follows

• We solve the following variable density flow problem and obtain:

u = u n -∆t(u n .∇)u n + ∆t ρ n µ∆ u - ∆t ρ n ∇p n+1 + ∆tg with ρ n = ρ s H n + ρ f (1 -H n )
• We compute the rigid velocity the rigid velocity corresponding to ũ:

u s = u G + w × r n , u G = Ω ρ n u H n dx Ω ρ n H n dx , ω = J -1 Ω ρ n (r n × u) H n dx
• We penalize the flow with this rigid velocity inside the solid body:

u n+1 -u ∆t = 1 η H n (u s -u n+1 ) ⇔ u n+1 = u + ∆t η H n u s 1 + ∆t η H n (54)
• We finally advect the solid with the rigid velocity:

H n+1 = H n -∆t u s • ∇H n
The variable-density flow equations are solved by a projection method using a MAC staggered grid: p and H share the same nodes, while the horizontal and vertical components of velocity are located respectively between a pressure node and its right or upper neighbor. This arrangement of the degrees of freedom ensures an algebraic divergence free vector field during the projection step, and prevents parasitic pressure modes. The advection equation is discretized by a fifth order WENO scheme, using the real velocity for the density and the rigid velocity field for H, ensuring a constant shape of the rigid body.

Note that we have used an implicit time discretization of the velocity penalization, which allows to use a very small penalization parameter η. Using an explicit method would require this value no to be smaller than ∆t. It can indeed be checked that a explicit method with η = ∆t essentially amounts to the projection method [START_REF] Patankar | A formulation for fast computations of rigid particulate flows[END_REF]. We will see below that using smaller values of η together with an implicit scheme has a significant effect on the accuracy of the method.

In order to numerically validate the penalization method, we consider the case of the sedimentation of a rigid cylinder in two dimensions (see [START_REF] Glowinski | A Fictitious Domain Approach to the Direct Numerical Simulation of Incompressible Viscous Flow past Moving Rigid Bodies: Application to Particulate Flow[END_REF], [START_REF] Coquerelle | A vortex level-set method for the two-way coupling of an incompressible fluid with colliding rigid bodies[END_REF]). The domain Ω = [0, 2] × [0, 6] is filled with an incompressible viscous fluid initially at rest, of density ρ f = 1 and viscosity µ f = 0.01.

The rigid cylinder of radius R = 0.125 and density ρ s = 1.5 is initially centered in (1, 4), and we apply the gravity force g = -980. In this experiment, the Reynolds number based on the cylinder diameter varies from 0 to about 250.

In order to verify how the rigid constraint is satisfied in the solid, we monitor at time t = 0.1 the L 2 -norm of the discrete deformation tensor defined by:

D(u) 2 L 2 (Ωs(t)) := ij H ij D 2 11 (u ij ) + 2D 2 12 (u ij ) + D 2 22 (u ij ) (∆x) 2
We fix ∆x = 1/256 and ∆t = 10 -4 , and compute this norm for values of η from 10 -4 to 10 -12 , at

t = 0.1.
The results presented in table 1 indicate a convergence of the penalization method, as far as the deformation is concerned, with first order in η. Note that in our decomposition of the penalization model, we have chosen to verify the divergence free constraint in the stage that precedes the computation of the rigid velocity. As a result, the solution at the end of a complete time step exhibits a boundary layer in the divergence of the velocity, proportional to ∆t η (u S -u)•n) L 2 (∂Ω S ) .

If we had chosen to apply the pressure correction after the equation (54) the divergence of the velocity would be algebraically zero at the end of the time step but the deformation inside the body would be significantly higher.

In figure 1 we show the profiles of the vertical velocity for several values of η, corresponding to a cross section at the center of the cylinder. We can observe that below η = 10 -8 one may consider that we obtained converged velocity results. The projection method of [START_REF] Patankar | A formulation for fast computations of rigid particulate flows[END_REF][START_REF] Sharma | A fast computation technique for the direct numerical simulation of rigid particulate flows[END_REF] of replacing the step (54) by u n+1 = u s . This corresponds to the following explicit discretization of the penalization term:

u n+1 -u ∆t = 1 η H n (u s -u)
with ∆t = η. We show the results obtained with these parameters and this discretization of the penalization term. As far as precision is concerned, one can notice the benefit of using larger penalization parameters combined with an implicit time discretization of the penalization term.

This confirms the results reported in [START_REF] Coquerelle | A vortex level-set method for the two-way coupling of an incompressible fluid with colliding rigid bodies[END_REF].

Conclusion

We have presented and analyzed a penalization method that extends the method of [START_REF] Ph | A penalization method to take into account obstacles in incompressible viscous flows[END_REF] to the case of a rigid body moving freely in an incompressible fluid. The proof is based on compactness arguments. Numerical illustrations have been provided to illustrate our convergence result. The benefit of using very large penalization parameters combined with an implicit time discretization of the penalization term, compared to the projection method [START_REF] Sharma | A fast computation technique for the direct numerical simulation of rigid particulate flows[END_REF] which corresponds to a particular explicit time discretization for the penalized equation, has been demonstrated. While this was not our primary goal, an outcome of our convergence study is an existence result for a weak formulation of the coupling between a rigid solid and a fluid. Let us shortly discuss how this result compares with existing ones [START_REF] Grandmont | Existence de solutions d'un problème de couplage fluidestructure bidimensionel instationnaire[END_REF][START_REF] Desjardins | Existence of Weak Solutions for the Motion of Rigid Bodies in a Viscous Fluid[END_REF][START_REF] Conca | Existence of solutions for the equations modelling the motion of a rigid body in a viscous fluid[END_REF][START_REF] San Martin | Global Weak Solutions for the Two-Dimensional Motion of Several Rigid Bodies in an Incompressible Viscous Fluid[END_REF]. In [START_REF] Grandmont | Existence de solutions d'un problème de couplage fluidestructure bidimensionel instationnaire[END_REF] local in time existence and uniqueness of strong solutions was proved. The Eulerian approach was developed in [START_REF] Desjardins | Existence of Weak Solutions for the Motion of Rigid Bodies in a Viscous Fluid[END_REF] where global in time existence of weak solutions was proved in dimension 2, without collisions. In the three-dimensional case, to our knowledge only local in time existence of weak solutions was obtained, since L 2 regularity of the time derivative of velocity was required (and therefore global existence would imply global existence of strong solutions). In [START_REF] Conca | Existence of solutions for the equations modelling the motion of a rigid body in a viscous fluid[END_REF][START_REF] Gunzburger | Global Existence of Weak Solutions for Viscous Incompressible Flows around a Moving Rigid Body in Three Dimensions[END_REF] the existence of global weak solution in three dimensions for one ball shaped solid, with possible collision with the boundary, was proved. In [START_REF] San Martin | Global Weak Solutions for the Two-Dimensional Motion of Several Rigid Bodies in an Incompressible Viscous Fluid[END_REF] the existence of global weak solutions for several rigid bodies with collisions was proved in dimension 2. Our results prove the existence of global in time weak solutions in three dimensions, before collision. By contrast with [START_REF] Conca | Existence of solutions for the equations modelling the motion of a rigid body in a viscous fluid[END_REF][START_REF] Gunzburger | Global Existence of Weak Solutions for Viscous Incompressible Flows around a Moving Rigid Body in Three Dimensions[END_REF], this result can easily be generalized to the case of several bodies by introducing indicator functions, rigid velocities and penalization terms corresponding to each body. To our knowledge, the existence of global in time weak solutions for several bodies with collisions is an open problem in three dimensions.

Appendix

This section si devoted to the proof of some technical lemmas that were used in section 3.

Proof of Lemma 3.3 Let φ(t) ∈ L 2 (Ω\Ω s,σ (t)), and (v(t), q(t)) ∈ H 1 (Ω\Ω s,σ (t))×L 2 (Ω\Ω s,σ (t))

solution of the Stokes problem      -∆v(t) + ∇q(t) = φ(t) on Ω\Ω s,σ (t), div v(t) = 0 on Ω\Ω s,σ (t), v(t) = 0 on ∂ Ω\Ω s,σ (t) .
Since we assumed a C 2 regularity on Ω\Ω 0 s this regularity is conserved through rigid motion and, for σ small enough (say σ < σ 0 for some σ 0 > 0) to Ω\Ω s,σ (t). The regularity results of Agmon-Douglis-Nirenberg on the linear Stokes problem (see [START_REF] Temam | Navier-Stokes equations and numerical analysis[END_REF], prop. 2.3. p. 35) give

(v(t), q(t)) ∈ H 2 (Ω\Ω s,σ (t)) × H 1 (Ω\Ω s,σ (t))
and there exists C > 0 such that

v(t) H 2 (Ω\Ωs,σ(t)) + q(t) H 1 (Ω\Ωs,σ(t)) ≤ C φ(t) L 2 (Ω\Ωs,σ(t)) .
Note that, with our definition of T , the constant C, which depends on the geometry of the domain boundary, can be taken independent of t ∈ [0, T ] and σ, provided σ 0 is taken small enough. We can then write

g(t)q(t) • nds ≤ g(t) L 2 (∂Ωs,σ(t)) ∇v(t) L 2 (∂Ωs,σ(t)) + q(t) L 2 (∂Ωs,σ(t)) ≤ C g(t) L 2 (∂Ωs,σ(t)) v(t) H 2 (Ω\Ωs,σ(t)) + q(t) H 1 (Ω\Ωs,σ(t)) ≤ C g(t) L 2 (∂Ωs,σ(t)) φ(t) L 2 (Ω\Ωs,σ(t)) ≤ C u(t) 1 2 L 2 (Ωs,σ(t)) ∇u(t) 1 2 L 2 (Ωs,σ(t)) φ(t) L 2 (Ω\Ωs,σ(t)) .
This proves the assertion.

Proof of Lemma 3.6

Step 1: We first show how to construct for a.e. t ∈]0, T [ a function v σ (., t) ∈ K r σ (t) such that lim σ→0 v σ (., t)u(., t) V r = 0 a.e. on ]0, T [. Let σ > 0 and v σ (., t) such that

         -∆v σ (., t) + ∇p(., t) = -∆u(., t) on Ω\Ω s,σ (t), div v σ (., t) = 0 on Ω\Ω s,σ (t), v σ (., t) = u s (., t) on ∂Ω s,σ (t), v σ (., t) = 0 on ∂Ω,
where

u s = 1 M Ω ρu H dx + J -1 Ω ρ(r × u)H dx × r.
By lemma 3.4, u(., t) = u s (., t) on Ω s (t). Extending v σ (., t) by u s (., t) in Ω s,σ (t), we have v σ (., t) ∈ K r σ (t). We set e σ (., t) = v σ (., t)u(., t). It satisfies

        
-∆e σ (., t) + ∇p(., t) = 0 on Ω\Ω s,σ (t), div e σ (., t) = 0 on Ω\Ω s,σ (t), e σ (., t) = u s (., t)u(., t) on ∂Ω s,σ (t), e σ (., t) = 0 on ∂Ω.

We extend e σ (., t) by u s (., t)u(., t) in Ω s,σ (t), so that e σ (., t) = 0 in Ω s (t).

We claim that lim σ→0 e σ (., t) L 2 (Ω) = 0 a.e. on ]0, T [.

In Ω s (t) e σ (., t) = 0, thus e σ (., t) Step 2: By definition of P r σ , P r σ u(., t)u(., t) V r ≤ v σ (., t)u(., t) V r thus the pointwise convergence on v σ we just obtained implies lim σ→0 P r σ u(., t)u(., t) V r = 0 a.e. on ]0, T [. (59)

Step 3: f σ : t → P r σ u(., t)u(., t) V r is measurable on [0, T ] and since 0 ∈ K r σ , f σ Therefore, thanks to lemma 3.2, lim σ→0 f σ L 2 (0,T ) = 0, which means lim σ→0 P r σ uu L 2 (0,T,V r ) = 0.

Proof of Lemma 3.7

Step 1: We construct for a.e. fixed t ∈ [0, T ] a function v ησ (., t) ∈ K r σ (t) such that lim σ→0 lim η→0 v ησ (., t)u η (., t) V r = 0 a.e. on ]0, T [. Let σ > 0 and v ησ (., t) solution of the following Stokes problem outside Ω η s,σ (t):          -∆v ησ (., t) + ∇p(., t) = -∆u η (., t) on Ω\Ω η s,σ (t), div v ησ (., t) = 0 on Ω\Ω η s,σ (t), v ησ (., t) = u η,s (., t) on ∂Ω η s,σ (t), v ησ (., t) = 0 on ∂Ω.

Extending v ησ (., t) by u η,s (., t) in Ω η s,σ (t), we have v ησ (., t) ∈ K r σ (t). We then introduce e ησ (., t) = v ησ (., t)u η (., t). It verifies          -∆e ησ (., t) + ∇p(., t) = 0 on Ω\Ω η s,σ (t), div e ησ (., t) = 0 on Ω\Ω η s,σ (t), e ησ (., t) = u η,s (., t)u η (., t) on ∂Ω η s,σ (t), e ησ (., t) = 0 on ∂Ω, and we extend it by u η,s (., t)u η (., t) in Ω η s,σ (t).

We claim that lim σ→0 lim η→0 e ησ (., t) L 2 (Ω) = 0 a.e. on ]0, T [.

From lemma 3.5, for a given σ > 0, there exists η 0 > 0 such that ∀η < η 0 , Ω η s (t) ⊂ Ω s,σ (t) and Ω s (t) ⊂ Ω η s,σ (t).

Let η < η 0 . We write e ησ (., t) ≤ C(η + η

1 2 + σ 2 ) 1 2 .
Therefore, for a fixed value of σ we can pass to the limit in η, and then pass to the limit in σ, to

This strong convergence can be turned into an almost everywhere in t convergence up to the extraction of a subsequence. The rest of the proof is adapted in a straightforward way from that of Lemma 3.6.

Figure 1 :

 1 Figure 1: Sedimentation of a two-dimensional cylinder, for ∆x = 1/256 and ∆t = 10 -4 . Vertical velocity in an horizontal cross-section through the center of the cylinder at t = 0.1 for several values of the penalization parameter.

e

  σ (., t) V r = 0 ∀r ∈ [1/2, 1[ a.e. on ]0, T [. (58)

L 2 r

 2 t) = 0 a.e. on [0, T ], f σ is measurable on [0, T ], f σ (0,T )≤ C with r < 1.

2 L 2 (

 22 Ω) dt = 0.

Table 1 :

 1 Sedimentation of a two dimensional cylinder. Errors on D(u) L 2 (Ωs(t)) and convergence orders at t = 0.1 for ∆x = 1/256 and ∆t = 10 -4 .

	essentially consists

  The integral of v • ∇p vanishes since v is divergence-free and vanishes on ∂ Ω\Ω s,σ (t) . Then using classical trace theorems in Sobolev spaces we get

	Ω\Ωs,σ(t)	w(t) • φ(t) dx = -	∂Ωs,σ(t)	g(t) •	∂v(t) ∂n	ds +	∂Ωs,σ(t)
	Ω\Ωs,σ(t)	w(t) • φ(t) dx = -	Ω\Ωs,σ(t)	w(t) • ∆v(t) dx +	Ω\Ωs,σ(t)	w(t) • ∇q(t) dx
			= -	∂(Ω\Ωs,σ)(t)	w(t) •	∂v(t) ∂n	ds +	Ω\Ωs,σ(t)	∇w(t) • ∇v(t) dx
				+	∂(Ω\Ωs,σ)(t)	w(t)q(t) • nds -	Ω\Ωs,σ(t)	div w(t)q(t) dx
			= -	∂Ωs,σ(t)	w(t) •	∂v(t) ∂n	ds +	∂(Ω\Ωs,σ)(t)	∂w(t) ∂n	• v(t)ds
				-	Ω\Ωs,σ(t)	∆w(t) • v(t) dx +	∂Ωs,σ(t)	w(t)q(t) • nds
			= -	∂Ωs,σ(t)	w(t) •	∂v(t) ∂n	ds +	Ω\Ωs,σ(t)	v(t) • ∇p(t) dx
				+	∂Ωs,σ(t)	w(t)q(t) • nds

  ≤ C e σ (., t) L 2 (Ωs(t)) ∇e σ (., t) L 2 (Ωs(t)) + σ 2 ∇e σ (., t)Next, as e σ (., t) = u s (., t)-u(., t) in Ω s,σ (t) and u(., t) and u s (., t) are in H 1 0 (Ω), we get ∇e σ (., t) L 2 (Ωs,σ(t)\Ωs(t)) ≤ ∇e σ (., t) L 2 (Ωs,σ(t)) ≤ ∇e σ (., t) L 2 (Ω) ≤ C, where C is independent of σ. This gives e σ (., t) L 2 (Ωs,σ(t)) = e σ (., t) L 2 (Ωs,σ(t)\Ωs(t)) ≤ Cσ.In order to prove that this convergence also holds in V r we first note thate σ (., t) H 1 (Ω) ≤ C a.e. on ]0, T [(56)as is readily seen from estimates on the Stokes problem verified by e σ . By interpolation (see e.g.[1], p. 135), we obtain e σ (., t) V r ≤ e σ (., t)

			2 L 2 (Ω) = e σ (., t)	2 L 2 (Ωs,σ(t)\Ωs(t)) + e σ (., t)	2 L 2 (Ω\Ωs,σ(t)) .
	Since Ω s,σ (t)\Ω s (t) has width 2σ, from the proof of lemma 5.10 of [10] we have a.e. on ]0, T [,
	e σ (., t)	2 L 2 (Ωs,σ(t)\Ωs(t))	≤ C e σ (., t)	2 L 2 (∂Ωs(t)) + σ 2 ∇e σ (., t)	2 L 2 (Ωs,σ(t)\Ωs(t))
						2
						L 2 (Ωs,σ(t)\Ωs(t))
			= Cσ 2 ∇e σ (., t)	2 L 2 (Ωs,σ(t)\Ωs(t))
	By Lemma 3.3 we thus get
		e σ (., t)
					1-r L 2 (Ω) e σ (., t)	r H 1 (Ω)

2

L 2 (Ω\Ωs,σ(t)) ≤C e σ (., t) L 2 (Ωs,σ(t)) ∇e σ (., t) L 2 (Ωs,σ(t)) ≤ Cσ.

Collecting the above estimates, we conclude that lim σ→0 e σ (., t) 2 L 2 (Ω) = 0.

2

  L 2 (Ω) = e ησ (., t) Since e ησ (., t) = u η,s (., t)u η (., t) in Ω s,σ (t) and u η , u η,s are bounded in L 2 (0, T ; H 1 0 (Ω)), we have

			0	T	∇e ησ (., t)	2 L 2 (Ωs,σ(t)) dt ≤ C.	(67)
	With (62) and (67) we are now in position to estimate the integral over [0, T ] of (64-66). By
	Cauchy-Schwarz inequality:
	0	T	e ησ (., t)	2 L 2 (Ωs,σ(t)\Ω η s (t)) dt ≤ C(η	1 2 + σ 2 ),

T 0 e ησ (., t)

2 L 2 (Ω\Ωs,σ(t)) dt ≤ C T 0 e ησ (., t) 2 L 2 (Ωs,σ(t)) dt 1 2