On discretization in time in simulations of particulate flows

We propose a time discretization scheme for a class of ordinary differential equations arising in simulations of fluid/particle flows. The scheme is intended to work robustly in the lubrication regime when the distance between two particles immersed in the fluid or between a particle and the wall tends to zero. The idea consists in introducing a small threshold for the particle-wall distance below which the real trajectory of the particle is replaced by an approximated one where the distance is kept equal to the threshold value. The error of this approximation is estimated both theoretically and by numerical experiments. Our time marching scheme can be easily incorporated into a full simulation method where the velocity of the fluid is obtained by a numerical solution to Stokes or Navier-Stokes equations. We also provide a derivation of the asymptotic expansion for the lubrication force (used in our numerical experiments) acting on a disk immersed in a Newtonian fluid and approaching the wall. The method of this derivation is new and can be easily adapted to other cases.


Introduction
One of the challenges for fluid/particle flow simulations is to provide an accurate resolution of the lubrication regime when the distance between two particles immersed in the fluid or between a particle and the wall becomes very small. Taking aside the problems related to the discretization in space (extremly high gradients of the velocity in the narrow gap between the particle and the wall, for example), we focus our attention in this article on the discretization in time. The lubrication regime is characterized by very high magnitude of the drag force (which can be referred to as the lubrication force in this case). Very small step sizes in time should be thus employed in order to obtain a physically acceptable solution. Our idea is to prohibit the particle from approaching too closely the wall during a simulation. We shall thus choose a threshold q s for the distance q between the particle and the wall and replace the true trajectory of the particle by an approximated one, in which the distance q is kept equal to q s until an eventual rebound of the particle from the wall. The moment of rebound is predicted using an auxiliary quantity (a crude approximation of the velocity) that is computed all along the period of time when the particle is stuck at the distance q s . This approach reminds the gluey particle model of [14,12] where q s is set to zero and the limit of vanishing viscosity is considered. However, our motivations are quite different from that behind the gluey particle model. This model is intended as a simple alternative to the standard governing equations of Navier-Stokes type, eventually corrected by taking into account the roughness of the particle surface. On the other hand, our approach is to take the standard fluid equations for granted (assuming the particle surface to be smooth) and to provide a tool for a robust time discretization of them. It means, in particular, that we would need an accurate enough method valid for any given value of the viscosity, not necessarily small. Note also, that our threshold q s will typically depend on the time step size, so that it is indeed a numerical device and it has no physical meaning.
The plan of the article is as follows: we start by reminding the governing equations in a general setting and by explaining in more detail the difficulties related to the simulations in the lubrication regime in the next section. Section 3 is the core of the paper. The idea of the threshold is rigorously introduced and studied there for the model ordinary differential equation representing the essence of the general setting in the simplest case of a circular (or spherical) particle approaching the wall. The discussion is held on the continuous level in Section 3. The discretization in time is introduced in Section 4 where several implementations of our idea are proposed on the discrete level followed by numerical experiments. We also include an appendix detailing a derivation of the asymptotic expansion for the lubrication force acting on a disk approaching the wall. The method of this derivation is new and can be easily adapted to other cases.
2 Motivations: governing equations for the fluid/particle flows and some difficulties arising in their simulations The general setting of this work is a study of the motion of a rigid particle immersed in a viscous, incompressible fluid with a particular emphasis on situations when the particle approaches the plane. To set the notations, we assume in general that the fluid (with the particle inside) fills a fixed domain Ω ∈ R d with d = 2 or 3, while the region occupied by the particle B t ⊂ Ω varies with time t. We denote the time-dependent fluid domain F t so that Ω = B t ∪ F t at any time t. Supposing that the inertial effects are negligible in the fluid and the no-slip conditions are valid on the boundaries of Ω and B t , the fluid motion is governed by the Stokes equations where u and p are the velocity and the pressure in the fluid, ν and ρ f are the viscosity and the density of the fluid, g is the external force, V = V(t) and ω = ω(t) are the translational and angular velocities of the rigid body B t , r = x−G is the vector pointing from the center of mass of the particle G to a point x on its boundary. The more realistic Navier-Stokes equations may be accomodated into the framework (1) by including the convective term u · ∇u into g. The fluid exerts a net force F and a torque T on the particle given by where D(u) stands for the symmetric gradient of u and n is the unit normal vector on ∂B t directed towards the fluid domain. Note that F and T are indeed functions of only the placement of the particle B t and its translational and angular velocities since the velocity u and the pressure p are uniquely determined in the fluid by the parameters B t , V and ω as the solution to the Stokes equations (1). Moreover, the dependence of F and T on V and ω is linear. Using these notations we write out the equations of motion of the particle as follows where m is the mass of the particle and I t is its inertia tensor, expressed in the fixed Cartesian frame and thus dependent on time. Equations (3) are coupled with the equations describing the propagation of the particle, i.e.Ġ = V for the center of mass andṙ i = ω×r i , i = 1, . . . , d, for the vectors r i fixed in the particle. The force F is a sum of the Archimedes force due to gravity and of the drag force which is purely hydrodynamic, i.e. obtained from (1) by setting g = 0. The particularity of the drag is that it tends very rapidly to ∞ when the particle approaches the wall, thus preventing collisions between them. Indeed, it has been proved in [8] (2D case), [9] (3D case), that a smooth rigid body embedded in a viscous fluid cannot touch the wall in finite time. In the regime of very small distance between the particle and the wall, the drag force is also known as the lubrication force and it is notoriously difficult to take into account in a numerical simulation.
It is noteworthy that these considerations are valid only for smooth surfaces; however, modelling surface roughness of the wall or particle is a much more delicate issue. Several experimental ( [18], [11], [10], [19]) and theoretical ( [10], [17]) works propose in this case to modify the expression for the lubrication drag force by introducing a shift to the distance to the wall, of magnitude strictly lower than the roughness size. This physically means that roughness decreases the dissipation in the system, and that the interaction is similar to that between equivalent smooth surfaces, located at some intermediate position, between the peaks and the valleys of asperities.
To simulate the motion of the particle governed by equation described above, one should be able to solve numerically the Stokes system (1) for any given particle position B t and for any given velocities V and ω. This is a formidable task in itself especially because the position of the particle is not known a priori but changes with time. We do not make precise the choice of the numerical method for the Stokes system. We just assume for the moment that the force F and the torque T can be computed for any B t , V, ω but this computation is in general very expensive. Our primary goal in this article is to devise an efficient discretization in time of (3) using as few as possible solutions of the Stokes system. The simplest idea is to use the following scheme We have introduced here the uniform grid in time t k = k∆t and have denoted the quantities computed at the time step t k by the superscript k. The idea behind the scheme (4)- (7) is to compute first the velocities by (4)-(5) and then to propagate the particle using the last available values of the velocities. The equations (4)-(5) are thus coupled together and also coupled with the solution of the Stokes system on a fixed geometry given by position of the particle B t k−1 on the previous time step t k−1 . The cost of such a computation is normally essentially the same as that of the Stokes system (1) with prescribed V and ω. We need thus one solution of the Stokes system per time step. This approach was successfully used, for example, in [13] in conjunction with a fictitious domain discretization of the Stokes system as in [6]. However, independently from the discretization in space, one would encounter problems when the particle approaches the wall. Indeed, in this lubrication regime the force F explodes and thus the scheme (4)-(5) is no longer valid unless an extremely low value for the time step is used that makes the simulation prohibitively expensive. A commonly used cure for this problem is to introduce short-range repulsion forces between the particle and the wall, as in [6], for example. However, the influence of these (not necessarily realistic) forces on the accuracy of a simulation is not well understood. Another simple idea is just to stop the particle when it tries to penetrate the wall during a numerical simulation. However, it is then not necessarily clear what criterion should be chosen to decide if the particle should eventually bounce off the wall and when should it happen. These questions have a partial answer in the articles [14,12] on the gluey particle model. It is shown there that the particle trajectory satisfies an integro-differential equation in the limit of vanishing viscosity, which is easy to discretize in time using moderate time steps and which predicts the moment of an eventual rebound from the wall. We pursue a similar idea in this article but our aim is to construct an approximated trajectory of the particle in the lubrication regime that would be accurate enough for any given value of the viscosity, not necessarily small.
3 A model ordinary differential equation with lubrication forces

The model
Let us consider the simplest setting of the problem described in the previous section: assume Ω ⊂ R 2 is the half-plane {(x, y), y > 0} and the particle is a disk of radius R. Let moreover g(t) = g(t)e 2 and assume the particle is at rest at the initial time. The x-component of the particle velocity and its angular velocity will then vanish at all time. The position of the particle is fully determined by its distance q from the bottom, as in Figure 1. The net force F is the sum of the drag, which is a function of q and V, linear in V, and of the Archimedes force: where n(q) is the drag coefficient computed by the Stokes equations (1). Denoting the y-component of the velocity by v, we are thus led to the following differential equations We are especially interested in the lubrication regime of small q. The asymptotic of n(q) when q → 0 is given by 3 √ 2πν R q 3 2 (see Appendix A) in the 2D case. After eliminating v from the system (8) and going to non-dimensional variables (see Appendix B for the details) we obtain the following equation for q(t) with ε = 3ν ρsR 3 2 2ρs g char (ρs−ρ f ) where ρ s is the density of the solid disk and g char is the characteristic value of g.
In the same way, we can consider the analogous three dimensional problem setting Ω to be the half-space {(x, y, z), z > 0} and the particle to be a ball of radius R. The asymptotic expression is well known in this case (cf. Remark 13) and is given by 6πV ν R 2 q . Performing the same non-dimensionalizations as in the 2D case, we arrive at the equation for q(t) (the distance from the particle to the wall): (10) Table 1: Typical values of the parameter ε in (9) or in (10) taking the viscosity and density of either water or glycerin and different value of the particle radius R and density ρ s .
We remind that equations of the type (9)-(10) are at the basis of the gluey particle model of [14,12]. The model consists in fact in considering the limit ε → 0, which is physically the limit of vanishing viscosity. We see, however, that ε is not necessarily small.

Estimates and an approximated solution for the model ODE
Consider the ordinary differential equation q = −n(q)q + g (11) or, equivalently, the system of first-order equations where n(q) is a given differentiable decreasing positive function on (0, ∞) with n(q) = N ′ (q) such that N(q) → −∞ as q → 0. Note that this is valid for n(q) = ε/q 3 2 or n(q) = ε/q, which give the asymptotic of the lubrication force in 2D and 3D respectively. We assume g ∈ L 1 loc (R + ) from now on. Under this hypothesis, Problem (11), completed with appropriate initial conditions, is well-posed, as proved for instance in [14,Prop. 1.1] for the case n(q) = ε/q. The proof is generalized without problem to any n(q) satisfying the above hypotheses, as mentioned in [14, Section 5.1]. (11) with initial conditions: Let us take some threshold value q s . We suppose that q goes below this value after some time t 1 when q(t) hits the threshold q s for the first time. Our aim is to find a suitable approximation of q after t 1 which enables to predict the time t 2 when q goes above the threshold q s without solving (11). To this end, we introducē for t ≥ t 1 and note thatv(t 2 ) =q(t 2 ). Indeed, integrating (11) from t 1 to t > t 1 shows thaṫ Setting here t = t 2 and noting that q(t 2 ) = q s gives the desired result. Our approximation of the trajectory q(t), denoted byq(t), stems from the assumption (verified afterwards in Proposition 3) that the velocityq(t 2 ) at the return point t 2 is small provided the threshold q s is small. Consequently, a good approximation for the time t 2 should be provided by the timet 2 defined as the first time larger than t 1 whenv(t) = 0. The construction of the approximated trajectory is hence the following: we first assume thatq(t) is the same as q(t) until the latter hits q s for the first time at t = t 1 . Next, the trajectoryq(t) is frozen until the time t =t 2 and resumes then again as a solution to (11) starting from q s with zero velocity:q (see Figure 2). In order to study the error committed by introducing the approximated trajectoryq(t), we start by inferring the following estimate: The timet 2 provides a lower bound for t 2 , i.e.t 2 ≤ t 2 . Proof. There are three cases. If t 2 does not exist, there is nothing to prove. If t 2 = t 1 this means t 1 is a local minimum for q so thatq(t 1 ) = 0 and t 1 = t 2 = t 2 . Finally, if t 2 exists and satisfies t 1 < t 2 , since q(t 1 ) = q(t 2 ) = q s the function q(t) has a local minimum t min ∈ (t 1 , t 2 ) whereq(t min ) = 0. Evaluating (15) at We have proved, in particular, that if the true trajectory q(t) ever returns to the values above the threshold q s , i.e. t 2 < ∞, the indicatort 2 will show it, i.e.t 2 < ∞. In Proposition 5, we compute an upper boundt 2 for t 2 . In particular, under suitable assumptions on g the time t 2 exists.
To see moreover thatt 2 may be a good approximation to t 2 we need, as mentioned before, to controlq(t 2 ). The following proposition proves thatq(t 2 ) is not too large with respect to q s . Proposition 3 Letq + (resp. g + ) be the positive part ofq (resp. g):q + = max(q, 0) (resp. g + = max(g, 0)). Then where · denotes the supremum on [t 1 , t 2 ]. In particular,q(t 2 ) ≥ 0, so thatq(t 2 ) ≤ g + n(qs) . Proof. Letq be non-negative on an interval [t min , t max ] ⊂ [t 1 , t 2 ]. Without loss of generality, we can thus assume that t min is a local minimum of q(t) and t max is either a local maximum or t max = t 2 . Then, we introduce If f (t 0 ) > 0 at some t 0 ∈ (t min , t max ) then f (t) rests positive on a small interval around t 0 so thatq < 0 by invoking (11) and thusq(t) is decreasing on this interval. Since q(t) is increasing on [t min , t max ] and n(q) is decreasing, we obtain that f (t) is decreasing locally The results above allow us to fully characterize the error of approximating q(t) byq(t) up to timet 2 and to bound it by certain quantities depending onq alone for time aftert 2 , at least in the special case when g(t) is given by with some positive constants t 0 , g + and a negative function g − (t). Indeed, we prove in the following proposition that the error τ = t 2 −t 2 is bounded by 1/n(q s ) and we remind that n(q s ) → ∞ as q s → 0.
Proposition 4 If g(t) is given by (18), then Moreover, both q(t) andq(t) are non-decreasing for t ≥t 2 and Proof. We first note that t 2 ≥t 2 ≥ t 0 . The first inequality here is already proved in Proposition 2. The second inequality follows fromv(t 2 ) = 0, which can be rewritten as Recalling Proposition 3, we see that which is the first inequality in (19). The second inequality in (19) is We now turn to the study of q(t) andq(t) for time t ≥t 2 , in order to prove that q(t) andq(t) are non-decreasing on this interval. We first observe that q(t) cannot have local maximums at t > t 0 . Indeed, at any local extremum t e > t 0 , equation (11) would implyq(t e ) = g + > 0. Since q(t) is non-decreasing att 2 > t 0 , it should be non-decreasing also everywhere on [t 2 , ∞). The same reasoning applies toq(t). The estimate (20) is now evident in the caset 2 ≤ t ≤ t 2 sinceq(t) ≥ q s and q(t) ≤ q s for such t.
In order to prove (20) in the case t ≥ t 2 we will show that q(t) is squeezed between q τ (t) andq(t) for t ≥ t 2 whereq τ (t) =q(t − τ ). Indeed bothq τ (t) and q(t) satisfy the same equation (11) with the same initial condition at t = t 2 (q(t 2 ) =q τ (t 2 ) = q s ) but possibly different initial velocitiesq(t 2 ) ≥ 0,q τ (t 2 ) = 0. Writing (11) forq τ (t) and q(t), taking the difference thereof and integrating from t 2 to t yieldṡ This shows that if the trajectories of q(t) andq τ (t) intersect at some time t ′ thenq(t ′ ) ≥ q τ (t ′ ) so that q(t) ≥q τ (t) at least for some time after t ′ . Since q(t) ≥q τ (t) holds also for some time after t 2 it should hold for all t ∈ [t 2 , ∞).
It remains to prove that q(t) ≤q(t) for all t ∈ [t 2 , ∞). To this end, we integrate (11) for q(t) from t 2 to t and that forq(t) fromt 2 to t. This yields with the aid of (21) q(t) =q(t 2 ) + N(q s ) − N(q(t)) + We see now that if q(t ′ ) =q(t ′ ) at some time t ′ ≥t 2 then alsoq(t ′ ) =q(t ′ ). It means by the uniqueness of solutions to (11) that the trajectories of q(t) andq(t) either coincide or do not intersect on [t 2 , ∞). Since q(t 2 ) ≤q(t 2 ) this implies q(t) ≤q(t) on [t 2 , ∞).
It is possible to relax the hypotheses of the last Proposition on the particular form of the function g in several ways, if we still remain in the case when q goes below the threshold value q s a finite number of times (in particular, g should be allowed to change sign only a finite number of times). Although we do not have an analogue of these results for a general g(t), we can always provide an easily computable upper bound for t 2 alongside the lower boundt 2 .
Proposition 5 Lett 2 be the first time aftert 2 such that Proof. Considerq(t) defined for t ≥t 2 as the solution tȯ We recall now the property (22) of the original solution and note N(q(t)) < N(q s ) for t ∈ (t 2 , t 2 ) so thatq (t) ≥ t t 2 g(s)ds ≥q(t).
Since q(t 2 ) >q(t 2 ), the trajectory of q(t) lies above that ofq(t) for t ∈ (t 2 , t 2 ) so that q(t) hits the threshold q s beforeq(t). In other words, t 2 <t 2 , wheret 2 is the first time aftert 2 such thatq(t) = q s . This definition oft 2 is equivalent to that in (24).

Three schemes for ODE (11)
We first describe a straightforward Euler discretization in time of (11) rewritten as the system (12). Introducing the time step ∆t and the discrete times t k = k∆t, k = 1, 2, . . . we thus consider the following Algorithm 1 is inspired by the real fluid-particle simulations in which one first finds the new velocity at each time step and then moves the particle with this velocity. An immediately evident drawback of the scheme in Algorithm 1 is that it does not necessarily provide a positive approximation q k . We remind that negative values of q are unphysical. Moreover, even if approximations q k remain positive but become tiny at some time steps, one would require very small stepsize to obtain an accurate solution.
Step 1 . For k = 1, 2 Remark 6 One may argue that this can be cured by resorting to a fully implicit scheme that couples the evaluation of the velocity with the displacement of the particle: Indeed, this scheme provides a positive solution for q, for example, in the important case n(q) = ε/q. To see this we eliminate v k from (26), which gives an equation for q k only The function of q k in the left-hand side is increasing form −∞ to ∞ when q k goes from 0 to ∞. This means that there is the unique positive solution q k for any given v k−1 and q k−1 . However, scheme (26) would be too expensive in a real fluid/particle simulation where n(q k ) should be recomputed after any change in q k through a numerical approximation of the Stokes equations.
Remark 7 Difficulties related to the discretization of these differential equations have already been pointed out in [14] (see Remark 3.1): it is shown that very small values of q can be reached (below the smallest real number which can be stored by standard numerical softwares and also below intermolecular distances). Different approaches have been proposed in previous related works to deal with similar problems: in [14, Section 5.2], the author suggests the introduction of a cut-off function for the microscopic distance γ = ε ln(q), defined in terms of the roughness of the surface (see also [12, Section 2.5]); in [12, Section 2.4.3], for the case of fluid/particle simulation, a constraint for the distance q is defined in terms of the mesh size of the fluid domain.
In order to overcome these difficulties, in this work we propose to introduce a small threshold value q s > 0, to replace the exact solution q(t) by the approximated oneq(t), as summarized in (16) on the continuous level. Discretization in time of this idea introducing also an approximation ofv(t) defined in (14) is detailed in Algorithm 2. Note that the exact threshold is replaced in Step 2 of this algorithm by the last available value of q k before Algorithm 2 (with an a priori fixed threshold) Step 0. Initialize (v 0 , q 0 ) and choose q s > 0.
Step 3. For k = k 2 + 1, k 2 + 2, . . . update (v k , q k ) as in Step 1. If q k goes again under q s at some time step, then reintroducev and keep switching between Steps 1 and 2 back and forth.
passing below q s , i.e. q k 1 −1 , so that the passage time t 1 is legitimately approximated by t k 1 −1 . As suggested by Proposition 4 a good choice for q s would be such that 1/n(q s ) = C∆t with some constant C so that the criterion for switching from Step 1 to Step 2 can be rewritten as n(q k ) ≥ 1/(C∆t). Indeed, extrapolating the result of Proposition 4 to a general g(t), we see then that the error of approximation of the exact return time t 2 bȳ t 2 ≈ k 2 ∆t is of order ∆t. Moreover, the error |q(t) − q(t)| is dominated by |q(t) −q(t − τ )| with τ of order ∆t. Thus the error caused by the introduction of the threshold should be of of the same order as that of Euler scheme itself. However, an optimal choice of the constant C above would require some a posteriori error indicators to control the error of the Euler discretization.
Generally speaking, a discretization of an ODE using the constant time step is often not optimal concerning the CPU time needed to achieve a desired accuracy. This is especially true for our ODE (11) since the coefficient n(q(t)) can change enormously during a simulation. One could try therefore to modify our algorithms by introducing a non constant stepsize chosen at each step using an error indicator. A general recipe for an optimal adaptation of the stepsize, according to [1] is to choose the stepsizes so that the discretization error per time step remains approximately constant during the whole simulation. As is well known, the error per time step in a Euler scheme like that of Algorithm 1, is given to the leading order by max(|q(t k )|, |v(t k )|)∆t 2 k /2 where ∆t k = t k − t k−1 is the stepsize on step k which can be varying. Assuming that (q k−1 , v k−1 ) and (q k , v k ) are Euler approximations at times t k−1 and t k respectively, the second derivatives of q and v can be computed approximately by deriving the equations in (12) with respect to time and replacing the Algorithm 3 (with adapted variable stepsize and automatically chosen threshold ) Step 0. Initialize (v 0 , q 0 ) for t 0 = 0, choose the first time step ∆t 1 > 0, the minimal tolerated time step ∆t min > 0 and the error tolerance per time step tol > 0. Step and calculate the approximation for the error committed on this time step by e k = max(|q(t k )|, |v(t k )|)∆t 2 k /2 whereq(t k ),v(t k ) are computed using (29)-(30). If e k ≤ tol and q k > 0 then we accept the just computed values (v k , q k ) and proceed to the next time step increasing the time step to ∆t k+1 = tol e k ∆t k . Otherwise, if e k > tol or q k ≤ 0, we reject the approximation (v k , q k ) and try to recompute it by (28) with a smaller time step ∆t k := min( tol e k ∆t k , ∆t k /2). If the new approximation (v k , q k ) is still not sufficiently accurate, we try to diminish the time step again and again until an acceptable approximation (v k , q k ) is computed and proceed only then to the next time step taking ∆t k+1 equal to the smallest value of ∆t k used on step k. However, if in the process of reducing ∆t k we come to a time step smaller than ∆t min at some k = k 1 , we abandon the calculation of (v k 1 , q k 1 ) and switch to Step 2 taking the initial valuesv k 1 −1 = v k 1 −1 , q s = q k 1 −1 and setting ∆t to the current value of ∆t k .
Step 3. For k = k 2 + 1, k 2 + 2, . . . update (v k , q k ) as in Step 1. If q k goes again under q s at some time step, then reintroducev and keep switching between Steps 1 and 2 back and forth.
missing derivatives by finite differences. This gives Thus, denoting E k = max(|q(t k )|, |v(t k )|)/2 where the second derivatives are replaced by the formulas above, the strategy to choose the time step would be to require E k ∆t 2 k = tol = const where tol is prescribed tolerance. As this recipe can lead to extremely small ∆t k we combine it with the threshold approximation by switching toq(t) only when ∆t k becomes smaller than some minimal admissible time step size. This idea is implemented in Algorithm 3.

Some numerical tests
As an illustration, we consider the ODE (11) with n(q) = ε/q 3/2 corresponding to the lubrication force in 2D. We report several numerical results setting for t > 2, and the initial conditions q(0) = 1,q(0) = 0. In all the cases, we have at our disposal a very accurate reference solution, which we call "exact" in what follows. We first compare Algorithm 1 (without threshold) vs. Algorithm 2 with the threshold such that n(q s ) = 1/(20∆t). The choice of the coefficient 20 in the last formula is somewhat arbitrary, but it works fine in our test cases. The first series of numerical experiments is performed taking ε = 0.1. The results are reported in Figure 3. Note that, although the solution obtained with Algorithm 1 is positive at ∆t = 0.01 and smaller, the introduction of the threshold in Algorithm 2 seems to enhance the quality of the solution at large time. Let us now turn to another series of numerical experiments taking ε = 10 −3 . The results are reported in Figure 4. Algorithm 1 is now unable to produce a physically acceptable solution even with very small ∆t = 0.0001. Algorithm 2, on the other hand, works fine. In order to observe better the quasi-contact region and the effect of introducing the threshold, we zoom on small distances by passing to a log-scale in Figure 5 (left). We report also some numerical experiments aiming at the determination of an optimal value of the threshold q s . As said already, our strategy is to choose it so that n(q s ) = 1/(C∆t) taking C = 20 in all the preceding experiments.   Figure 4, i.e. varying ∆t and choosing q s so that n(q s ) = 1/(20∆t). Right: same test with ∆t = 0.001 fixed and different values of q s calculated so that n(q s ) = 1/(C∆t) and C varying in the range [1,2000]. and taking three different values for q s by varying C. These results confirm two general observations: taking q s too small may deteriorate the accuracy of the solution (indeed, the original goal of the introduction of q s was to avoid the extremely small values of q). On the other hand, taking q s too large may unnecessarily perturb the solution in the regions where a straight discretization could give more precise results. Thus, there should be an optimal choice of q s which seems to be not too far from the formula mentioned above with C = 20. We conclude by a series of experiments using Algorithm 3 with adapted stepsize. We take again the same governing equation with ε = 10 −3 and run Algorithm 3 setting the tolerance to tol = 10 −5 and taking a number of values for ∆t min . The results are reported in Figure 6 with the evolution of q(t k ) on the left and that of ∆t k on the right. We give in particular the results with ∆t min = 0 so that Algorithm 3 is reduced to a standard time marching scheme with automatically adapted step sizes. It is interesting to note that step sizes ∆t of the order 10 −4 are actually sufficient to satisfy our error tolerance criterion almost everywhere apart from a very small range of time around 1 where the Algorithm chooses ∆t of the order 10 −8 . Setting ∆t min to 10 −3 or 10 −4 avoid such small time steps and gives fairly good results.

Acknowledgements.
We would like to thank the organizers of the conference ECM'09 (International Conference of Microfluidics and Complex flow) Lassaad El Asmi and Mourad Ismail for giving us the opportunity to present this work. We had a number of useful discussions during this conference, especially with Aline Lefebvre-Lepot and François Feuillebois, to whom we are very grateful.

A An estimate for the drag (lubrication) force
This part is devoted to computing the first order expansion of the viscous drag exerted by a fluid on a disk approaching a plane wall. Let Ω ⊂ R 2 be the half plane Ω = {(x, y) with y > 0}, the boundary of Ω is ∂Ω = {(x, 0), x ∈ R} (see Figure 1). We assume that the fluid fills the domain F q := Ω \ B q where B q := B((0, R + q), R) is the domain occupied by the solid disk of radius R > 0 and placed at the distance q > 0 from the boundary is ∂Ω. Our notations are lightly different here from that of Section 2. Indeed, domains F q and B q are now indexed by q rather than t since the eventual dependence of q on time does not interest us in this Appendix. Assume moreover that the disk has a prescribed velocity V = V e 2 and zero angular velocity. The motion of the fluid is governed by Stokes equations (1) with g = 0 and zero boundary conditions at the infinity for u. It is classical that (1) is well-posed in a suitable framework and has a unique solution (u, p) (the pressure p is defined up to a constant). In particular, the viscous drag force F exerted by a fluid flow (u, p) on B q and computed by (2) is well-defined as a function of V and q. Our aim is to prove the following result: with ε(q) → 0 when q → 0.

Remark 9
The equivalent result in the three-dimensional setting is well-known but it is new in two dimensions. Whatever the novelty of this result, the main contribution in this section is that we provide a new way for computing the first order expansion of the viscous drag. This new method is more robust than the known computations [2,3,15,16]. In particular, we claim it extends to the three-dimensional setting with small changes and can be adapted to other boundary conditions. As an illustration, our method is a welldesigned tool for estimating the distance between solid bodies in solutions to the fluidstructure interaction system with the full newtonian Navier Stokes equation on the fluid domain [5,8,9].
The proof of Theorem 8 is divided into three steps. First, we recall the variational formulation for the solution (u, p) to (1) and apply it to compute the associated drag coefficient n(q) in the expression for viscous drag F = −n(q)V. Secondly, we deduce from the variational formulation the lower and upper bounds for n(q). We conclude the proof by computing an asymptotic expansion of the bounds of this range.
1. The variational formulation for the Stokes system (1). As (u, p) in (1) and the formula (2) for the drag depend linearly on V, we set V = e 2 in what follows. It is classical to extend the velocity u to the whole domain Ω by setting it to e 2 inside B q and to look for u in the following space D 1,2 0 (Ω) := w ∈ L 1 loc (Ω) with ∇w ∈ L 2 (Ω), div w = 0 in Ω and w = 0 on ∂Ω . The solution of (1) is associated with the minimization of Dirichlet integral D(w) = Ω |∇(w)(x, y)| 2 dxdy over the subset Y q of D 1,2 0 (Ω): Let us denote this minimum by E D (q), i.e.
As the Dirichlet integral D is a strictly convex functional on Y q it has a unique minimizer u q for any q > 0. It is easy to see that u q gives the solution of (1). Indeed, we have for any w ∈ C ∞ c (F q ) such that div w = 0, Ω ∇u(x, y) : ∇w(x, y)dxdy = 0, which is the variational formulation of (1). The pressure p q ∈ L 2 loc (Ω) can be then recovered as the Lagrange multiplier corresponding to the constraint div u = 0. We notice that p q is unique up to a constant and the ellipticity of the Stokes problem implies that (u q , p q ) is smooth on F q so that it furnishes a classical solution to (1). We refer to [4] for more details and also for a proof of the converse implication.
Given w ∈ C ∞ c (Ω) such that div w = 0 and w = e 2 on B q , a straightforward integration by parts yields (recall that n is the unit normal looking into F q ): Taking a suitable family of approximation of u q (see [4] for details) we obtain in the limit: The drag F is parallel to the velocity V for symmetry reasons, i.e. F = (F · e 2 )e 2 . We have thus the following result.
Lemma 10 For any given q > 0 and V parallel to e 2 , the drag force is given by F = −n(q)V with n(q) = νE D (q).
2. Upper and lower bounds for E D (q). We apply the above variational formulation to bound E D (q). To this end, given q > 0, we first prove that the Dirichlet integral of any u ∈ Y q is greater than some m D (q) depending only on q. This gives us a lower bound for E D (q). Then, we construct a suitableũ q ∈ Y q and compute its Dirichlet integral M D (q). This gives an upper bound for E D (q). We finally compare m D and M D for small values of q to prove Theorem 8.
On the other hand, we compute the Dirichlet integral D(u) with respect to ψ. This yields: We denote I(ψ) the integral on the right-hand side of the above inequality. Thus, D(u) is greater than the minimum of I(ψ) over smooth ψ satisfying the boundary conditions (32). This minimum is computed in the following lemma: Lemma 11 Given q > 0 and ψ ∈ C ∞ (Ω q ), satisfying boundary conditions (32) there holds I(ψ) ≥ I(ψ q ) where : , ∀ (x, y) ∈ Ω q .
3. Asymptotic expansion of m D (q) and M D (q).
Proof. As already seen in the proof of the preceding Lemma, m D (q) is given by the integral (34) with C = 0 and η(t) = t 2 (3 − 2t). Integrating with respect to t and performing a change of variables x → √ qx yields m D (q) = 12 Expanding to the first order, we have: We remark moreover that standard geometric arguments imply that δ q (x) ≥ q + x 2 /2 or δ q ( √ qx)/q ≥ 1 + x 2 /2 so that we can apply the Lebesgue theorem to obtain: One can also check that ∂ yy ψ q is the second derivative of ψ q which diverges the fastest when q goes to 0. Thus, M D (q) − m D (q) = o(m D (q)) as q → 0 so that (35) holds also for M D (q), which proves Theorem 8.

Remark 13
The same method can be easily adapted to the three dimensional case. More precisely, the result of Lemma 10 is still valid : the viscous drag exerted by a fluid on a sphere can be computed as the result of a minimization problem of the same type as the one given by equation (31). Next, since the problem is axisymmetric, one can use cylindrical coordinates and perform the same asymptotical analysis as in the 2D case, in order to obtain that with ε(q) → 0 when q → 0.
Numerical evaluation of the drag force.
In order to compare the theoretical estimate given by Proposition 8 to a numerical evaluation of the force F, we carry out the following simulations: we begin with solving the Stokes problem in a fluid with viscosity ν = 1 and then compute the corresponding force exerted by the fluid on a particle of radius R = 0.1 situated above a plane at a distance q. The computations are performed by using the Finite-Element solver FreeFem++ [7] and the results are plotted in Figure 7 (dashed line with squares). They agree with the asymptotic expansion F ∼ 3 √ 2πV ν R q 3 2 (solid line).  Let us first consider the 2D case. Eliminating the velocity v from (8) and specifying m = ρ s πR 2 with ρ s the density of the solid, we arrive at In order to understand the typical values of parameters in this ODE, we should pass to non-dimensional variables, which can be introduced as follows: where the radius R is used as the length scale, the time scale is denoted by T and g char is a typical value of the external force density, so that g char ∼ 10m/s 2 . We choose then the time scale T so that the non-dimensional external force becomes of order 1, i.e.