Numerical simulation of rigid particles in Stokes flow: lubrication correction for general shapes of particles.

We address the problem of numerical simulation of suspensions of rigid particles in a Stokes ﬂow. We focus on the inclusion of the singular short range interaction eﬀects (lubrication eﬀects) in the simulations when the particles come close one to another. The problem is solved without introducing new hypothesis nor model. As in [11], the key idea is to decompose the velocity and pressure ﬂows in a sum of a singular and a regular part. In this article, the singular part is computed using an explicit asymptotic expansion of the solution when the distance goes to zero. This expansion is similar to the asymptotic expansion proposed in [6] but is more appropriate for numerical simulations of suspensions. It can be computed for any locally convex (that is the particles have to be convex close to the contact point) and regular shape of particles. Using [6] as an intermediate result, we prove that the remaining part is regular in the sense that it is bounded independently of the distance. As a consequence, only a small number of degrees of freedom are necessary to obtain accurate results. The method is tested in dimension 2 for clusters of two or three aligned particles with general rigid velocities. We show that, as expected, the convergence is independent of the distance.


Introduction
We consider in this article suspensions of rigid particles in a Stokes flow.Various domains of application can be modeled and studied using such systems: concretes or reinforced plastics in industry, silting up in natural environment, micro-swimmers or blood tests in biology, wastewater treatment in ecology...The applications we have in mind are for example migration in couette flows, dynamics of aggregates, or the study of the interphase stress contribution to the total stress.Unfortunately, the behavior of such systems is still partially misunderstood and precise numerical simulations are needed to understand the mechanisms behind their rheological behavior as well as their origin.A large amount of research has been developed to design numerical tools in order to study the motion of rigid bodies embedded in a Stokes flow.One can see for example [13] for a finite-element method, [1,16] for multipole decomposition methods, [4] for finite-difference methods or [7,8] for lattice Boltzmann simulations in case of a Navier-Stokes fluid.
In view of the previously cited applications, we have to consider dense systems for which the particles are close to the maximal packing and therefore close to contact.It consists in a challenging domain in which very few results are available.From a computational point of view, we have to consider suspensions with densities from 30% to 50%, involving several thousands of particles.It is well known, that, for such physical studies, one needs to compute precisely the velocity and pressure fields in the fluid domain.In particular, one has to take into account carefully close interactions between particles (such as lubrication or contact) and their feedback on the flow.The objective of this paper is to propose a method allowing to compute a precise approximation of the velocity and pressure fields using coarse meshes, even for small distances between the particles.Moreover, we design a method allowing to deal with non-spherical particles.
Equations of the problem.
We consider Ω = R 3 and denote by (B i ) i=1,...,N N rigid particles in Ω.We assume that the fluid domain F = Ω \ ∪B i is filled with a Newtonian fluid governed by the Stokes equations.
We consider in the following the corresponding Dirichlet problem.We suppose that the velocity on the boundary of each particle is given: the velocity of particle i is denoted by u i ∈ H 1/2 (∂B i ) and is supposed to satisfy the incompressibility condition ∂Bi u i • n = 0, ∀i = 1, . . ., N. ( Note that this condition is fulfilled if particle i undergoes a rigid motion: The problem is written as follows, Find (u, p) ∈ H 1 (F)×L 2 0 (F) such that −∆u + ∇p = 0, in F; div u = 0, in F; u = u i , on ∂B i , i = 1, . . ., N ; ( where L 2 0 (F) = {q ∈ L 2 (F) : F q = 0}.Since we suppose that Ω = R 3 , we add the following vanishing condition at infinity lim |(x,y,z)|→∞ u(x, y, z) = 0.
This, together with assumption (1) on u i , ensures that system (2) admits a unique solution.We denote by St(u 1 , . . ., u N ) this solution.
Remark 1. 1 We focus here on the Dirichlet problem: the velocities of the particles are supposed to be known and the objective is to compute precisely the corresponding pressure and velocity fields.This is a first step towards the simulation of suspensions, in which the velocity of the particles is part of the unknowns of the Stokes problem.Indeed, to compute these unknown velocities, some methods are based on the computation of the inverse of the resistance matrix (which provides the forces for given velocities).The computation of this resistance matrix is achieved solving several Dirichlet problems.Moreover, the study of dense suspensions in three dimensions leading to huge numerical problems, lot of other methods are based on iterative algorithms.To implement these algorithms, one has to define the matrix-product vector which, again, is typically based on a Dirichlet problem.
A singular problem due to the lubrication effect.
The main difficulty to solve problem (2) in case of dense suspensions is to be able to deal with the so called lubrication effect: when two particles get close one to another, their relative velocities, together with the incompressibility constraint, induce large variations in the velocity fields and high values for the pressure field, located in the thin gap between the two particles.For example, if d is the distance between the particles and if x // denotes the coordinates in a plane orthogonal to the line of the centers of the spheres with origin at the intersection with this line, the density of the pressure asymptotically behaves as p(x // / √ d)/d 2 when d goes to zero.We provide on Figure 1 an example of velocity and pressure fields for two particles moving towards each other (the detail of the corresponding configuration will be given in Subsection 3.2.1).
This behavior is said to be singular in the sense that the norm of the solutions are not bounded when the distance goes to zero.As a consequence, for a given number of degrees of freedom, the error when computing a numerical solution to the problem is increasing when the distance between the particles goes to zero.It is clear that, to properly capture the relevant phenomenon in numerical simulations, one will need to use a high number of degrees of freedom.For example, in the finite-element simulations presented Section 3.2.1,we can observe that we need about 8 elements in the gap between the particles to capture the force exerted on the particles in the direction along the line of their centers.This becomes very expensive for 3-dimensional dense suspensions (both for CPU time and memory reasons).However, in such situation, most of the particles are close and the effect of lubrication dominates the global behavior of the system.The challenge is then to find a suitable numerical strategy that takes into account the lubrication effects, without increasing the size of the mesh.
A singular/regular decomposition of the forces.
Suppose you want to compute the forces (F i ) i=1...N exerted on the particles: for example to compute the resistance matrix associated to the problem.A widely used strategy has been proposed in [3] in the context of Stokesian Dynamics simulations.The Stokesian Dynamics method is based on the computation of the resistance matrix at each time step.The authors propose a method to compute accurately the forces exerted on the particles, in order to be able to compute a precise resistance matrix, taking lubrication into account.
The founding remark is the following: if two particles are embedded in an infinite fluid, have the same radius and have opposite velocities with unit norm, the forces exerted on these particles only depend on the distance d and can be written as a function of one variable: d → F 2sph i (d).Since they only depend on a unique parameter, they can be easily estimated using both an asymptotic development when the distance goes to zero (see e.g.[2]) and a mean square approach for larger distances (very precise computations are achieved off-line for a discrete set of distances to construct an approximation for any value of the distance).Now, consider a general configuration of N particles.The singular part of the force F sing i acting on particle i is estimated as the sum of the contributions of each pair of close particles including particle i, each of these contributions being estimated using the previous off-line computations.The singular part of the forces being estimated, one then has to compute the remainder of the force F reg i , due to multi-particle interactions, to the physical boundaries of the domain, etc...This regular part can be approximated using coarse meshes.We refer to [3] for more details.Finally, the method can be rewritten as computing where F sing i is computed precisely from off-line computations and F reg i is computed on-line, using a small number of degrees of freedom.This method has been used more recently to include lubrication in the computation of the resistance matrix using direct solvers [4].
Such a method allows to take lubrication into account when computing the forces exerted on the particles.It converges when the number of degrees of freedom used to compute the regular part of the force increases.However, due to the off-line computations, it is restricted to spherical mono-disperse particles.Indeed, for polydisperse spherical particles, computing the singular solution would necessitate an approximation depending on 2 parameters: the distance and the ratio of the two radii.Concerning non-spherical particles, the singular part depends on the distance, the geometry of the particles in the gap and the orientation of the particles...Such estimations have never been achieved because too costly... Another kind of correction of the forces has been proposed more recently in [9], where an explicit asymptotic development of the local distribution of the force on the boundary of the particles has been used.This allows to free the authors from the off-line estimation step so that they can deal with non-spherical particles.However, the lubrication correction is added to the force computed by the direct numerical solver and therefore, when the mesh step goes to zero, numerical parameter tuning is necessary to avoid computing twice the force.
Note also that, although these methods compute precisely the total force exerted on the particles, they do not correct the velocity and pressure flows.
A new singular/regular decomposition of the flow.
The key point to take lubrication into account on the whole flow is to decompose the velocity and pressure fields into a regular and a singular part: u = u reg + u sing and p = p reg + p sing .Using the linearity of Stokes problem and assuming that the singular field is known, it remains to approximate the solution to the following regular problem: The aim is to construct a suitable singular field u sing such that: • u sing is known analytically; • the corresponding solution u reg is bounded regardless of the distance, that is there exists C > 0 independent of d such that Doing so, the right-hand side of (4) can be computed precisely and its solution can be approximated using a small number of degrees of freedom.A first method to compute such a singular field was proposed in [11].The authors followed the idea proposed for the Stokesian Dynamics: they considered the singular solution (u sing , p sing ) to the Stokes problem for two spherical particles with the same radius and having opposite velocities with unit norm.This solution only depends on the distance d between the particles: d → (u sing (d), p sing (d)).They estimate off-line the coefficients of the solution on a suitable basis.The remaining regular part of the solution is then approximated using a small number of degrees of freedom.Doing so, the global flow u = u sing + u reg is computed precisely, including lubrication effects, even for low numbers of degrees of freedom.However, the method is strongly dependent on the estimation of the singular field which can be achieved only if it depends on a single parameter.This reduces the use of the method to mono-disperse spherical suspensions.
A way to avoid off-line computations would be to derive an explicit analytical expression for the singular field.One could think to use the exact expressions of the field solution to the problem with the two particles approaching one to another.Such expressions do not exist for any kind of particles but are available for pairs of spheres (see e.g.[12] for two equal spheres or more recently [5] for two unequal spheres).They are expressed as infinite series which have to be truncated to be used in a numerical code.Unfortunately, for a given precision, the order of truncation increases when the distance between the particles goes to zero.This implies parameter tuning to choose the order of truncation and makes such a method potentially costly for suspensions involving many pairs of close particles.
To obtain a method numerically tractable, one can recall that, since (u reg , p reg ) is solution to the Stokes problem (4) with (u sing , p sing ) in the right-hand side, it is sufficient to find (u sing , p sing ) for which the following quantities are bounded when the distance goes to zero: together with u i − u sing in H −1/2 (∂B i ) for i = 1 . . .N .As a consequence, there is no need to derive the expression of an exact singular flow in the whole domain for any distance: it is sufficient to use an asymptotic expression when the distance goes to zero.A first step towards such a choice for the singular field was proposed in [10].In that article, the singular field is based on an explicit asymptotic expansion of the flow (u in , q in ) in the gap between the particles.The expansion is simply truncated to be used as the singular field.Doing so, the author shows that, using a finiteelement discretization for the regular problem, she obtains a good approximation of the whole flow in case of two spherical particles.However, the numerical results still depend on the distance.
In the present article, we propose a new explicit expression for the singular field which is based on the same asymptotic flow (u in , q in ).Our asymptotic field is similar to the one proposed in [6], which is modified to be more tractable from a numerical point of view: the singular field we propose presents less oscillations and vanishes outside the gap between the particles.Here, we use estimations proved in [6] as an intermediate result to show that our remaining regular field is bounded independently of the distance.Doing so, we obtain numerical results for which the error is independent of the distance which allows to use coarse meshes, independently of the size of the gap.Moreover, tests for non-spherical particles are achieved, showing that the method can be extended to more general suspensions.
Note that a decomposition of the solution has also been proposed in [14] in the context of boundary integral equations: the force density is decomposed in a singular and a regular part.The method proposed needs parameter tuning to make it converge when the number of degrees of freedom increases.Here, we propose a decomposition of the velocity and pressure fields so that the method is adapted to volumic numerical solvers (finite-element, finite-differences...).Moreover, no parameter tuning is necessary to make the method converge.
Outline of the paper.
The outline of the article is the following: after this introduction, the second section is dedicated to the case of two single particles.We explain how we construct the singular field (u sing , p sing ) and how the results of [6] are used to bound the remaining regular field.Numerical simulations are presented in Section 3 to validate the method.In Section 4, we investigate the case of three aligned particles, to show how the method can be extended to more than one singularity.To finish, we discuss in Section 5 the extension to more general suspensions.

The two particles case
We consider here the two particles case.Using the linearity of Stokes problem, one can decompose the boundary condition and write As a consequence, it is sufficient to solve accurately the problem involving an obstacle and a moving particle, which is the case treated in [6].Following this reference, let us rewrite this elementary problem and precise the notations.
We denote by B the obstacle and B the moving particle with velocity u .

Figure 2: The two particles case
The problem is the following Without loss of generality we define the origin of coordinates as the unique point of ∂B that achieves the distance d between ∂B and ∂B (see Figure 2).Let e ⊥ be the unit normal outward to ∂B at the origin and we assume that e ⊥ = e z .All along the paper we will use the notations ⊥ and // to distinguish the normal and the tangential components: • for any vector field v ∈ R 3 we define its normal and tangential components as follows, • we note ∇ // (respectively div // ) the tangential component (that is the component corresponding to e x and e y ) of operator ∇ (respectively div ) and ∇ ⊥ its normal component.
We also suppose there exists a parametrization of the boundaries ∂B and ∂B close to the origin: there exist L > 0 and smooth functions γ b and γ t on the two-dimensionnal ball centered at 0 with radius L denoted by B(0, L) such that for any (x, y) ∈ B(0, L), Classically, we suppose that there is a "single non-degenerate contact" in B(0, L), that is there exists a unique contact point in B(0, L) × R when the distance d goes to 0 (we refer the reader to [6] for the corresponding assumptions on the functions γ t and γ b ).
Finally, the gap between the particles is denoted by G L : The assumption of "single non-degenerate contact" implies that one needs the geometry of the particles to be locally convex close to the contact point.

Explicit asymptotic field proposed in [6]
The expressions of the singular field we use in the paper and defined in Section 2.2 are based on an asymptotic expansion of the flow when the distance goes to zero.Formal computations can be found in [2] where the author estimates the forces exerted on two particles almost in contact and undergoing a rigid motion.One can extend the corresponding computations to obtain a formal development of the flow in the narrow gap between the particles.These computations have been rigorously derived and generalized to any boundary condition on the particles in [6].We recall here the expression of the asymptotic field proposed in [6] and the estimations provided in the article.Indeed, this asymptotic field, together with the corresponding estimations, will be used in the Section 2.3 as an intermediate field in the computations in order to estimate u reg .

A local asymptotic development of the flow
Following [2,6], we provide here an explicit expression for the singular field in the gap G L between the particles.
• Step 1. Singular part of the boundary condition.Using the linearity of Stokes equations, the first step is to decompose the boundary condition u (rigid or not) into a sum of two boundary conditions.The first part, denoted by u as , is supposed to be the singular one, leading to singular behaviors of the flow when the distance goes to zero.It has to include the whole singularities of the boundary condition in the sense that the remaining part has to generate flows with a bounded norm when the distance goes to zero.By [2,6] the singular part can be chosen as • Step 2. Local singular pressure field.Denoting the distance between the boundaries of the particles by the asymptotic pressure field q in does not depend on z and is defined in B(0, L) as the solution to the following 2-dimensional problem: • Step 3. Local singular velocity field.From the 2-dimensional scalar field q in defined in B(0, L), the corresponding asymptotic velocity field is where, for any 2-dimensional scalar field q defined in B(0, L) the 3-dimensional vector field w[q] is given in G L by Note that these formula provide a divergence free velocity field.

Truncation and singular field choice in [6]
To construct a singular field in F, a regular truncation function χ L is used in [6], which does not depend on z, is equal to 1 in B(0, L/2) and vanishes outside B(0, L).
• The singular pressure field in G L .The singular pressure field is simply chosen as the truncation of the local singular field q in : Note that since the pressure field p sing does not depend on z, one has ∇ ⊥ p sing = 0.
• The singular velocity field in G L .In order to choose the singular velocity field, one has to recall that, to bound (u reg , p reg ) independently of the distance, we need to bound ∆u sing − ∇p sing in H −1 (see ( 5)).Since ∇ ⊥ p sing = 0, the tangential and normal components are studied separately.Following expression (9) for the local singular velocity field from the local pressure field, one can try to compute a singular velocity field with the following form: u sing = w[q] for a given q to be chosen.In that case, the normal component of the first term to be bounded in (5) rewrites: and its tangential component, using equation (10a): In view of these computations, the authors in [6] choose to use the singular pressure field p sing both for the tangential and normal part of the velocity field in G L : Indeed, this choice cancels the gradient of the singular pressure field in the right-hand side of ( 12) and ( 13), this term being the most singular one.
• Extension to F. Then the velocity and pressure fields are naturally extended by zero on the whole domain F. Note that v in does not match the boundary condition u on the moving particle: we have • Modification of the velocity field oustide G L/4 to match the boundary condition.This velocity field is then modified outside G L/4 : it is written as where v ext vanishes in G L/4 and is designed to obtain a velocity field v as in F which is divergence free and satisfies the same boundary conditions as u.This correction v ext has to be bounded independently of the distance.Here is how it is constructed: -The authors first define an outer domain Ω independent of d (see Figure 3) such that: -In order to match the boundary condition on ∂B , the authors define v ext as the solution to the following Stokes problem in Ω , with the correcting boundary condition which is extended by zero to ∂Ω .The solution v ext is then also extended by zero to the whole domain F.
Remark 2.2 Note that the asymptotic field v as does not vanish outside the gap G L .

Asymptotic results from [6]
Let us now recall three estimations proved in [6] that will be useful to estimate u reg .
First, the authors provide the following estimation for the inner field [6, proof of lemma 20]: there exists C > 0 independent of d such that, Note that, this estimation is useful to bound in G L the right-hand side of the Stokes problem solved by their regular field (see (14) together with ( 12) and ( 13)).
The authors also provide the following estimation for the external field [6, Lemma 18]: there exists C > 0 independent of d such that, Finally, it is proved in [6, Lemma 4] that the field v as is a good approximation of u when the distance goes to zero: there exists C > 0 independent of d, but depending on γ b , γ t , ∂B and u , such that Moreover, if the function γ is radial one has,

A new explicit expression for the singular field
In this section, we explain how we use the results of [6] to obtain an expression for the singular field that is tractable for numerical simulations of multi-particle suspensions.From the local asymptotic field (u in , q in ) in G L , we now want to define the singular field (u sing , p sing ) ∈ H 1 (F) × L 2 0 (F) in the whole fluid domain F. First of all, note that, to make it tractable for numerical simulations of multi-particle suspensions, we need that the singular velocity field vanishes far from the gap G L .Indeed, the right-hand side and the boundary condition for u reg in (4) depend on u sing .Computing the right-hand side of the regular problem everywhere in F together with the value of the singularity generated by a pair of particles on the boundary of all the other particles in the suspension would be too costly.This prevents using the asymptotic field v as = v in + v ext of [6] for general suspensions involving many couples of close particles.Note also that a singular field vanishing outside the gaps would allow the parallelization of each pair of particles' contribution.
The local singular part of the asymptotic field v in vanishes outside G L so that one could think to use it as a singular field.However, its normal component v in,⊥ involves many derivatives of the truncation function, which increases the number of computations and induces (controlled) oscillations in the singular flow.
We propose in the following a new choice for the singular field.As in [6], it is based on q in defined in (8).The choice for the singular pressure field is the same as the one proposed in the previous section and given by (11).The difference relies in the choice for the singular velocity field.
• The new singular velocity field in G L .According to (12) and (13) we see that to cancel the gradient of the singular pressure field, it is in fact sufficient to choose the singular velocity field consistently with p sing only for the tangential part: u sing // = w[p sing ] // .Concerning u sing ⊥ , to avoid the computation of too many derivatives of χ L , contrary to [6] we simply consider the truncation of w[q in ] // .To sum up, our choice for the singular field in G L is: • Extension to F. As in the previous section the velocity and pressure fields are naturally extended to F by zero.
We plot on Figure 4 the singular field corresponding to the reference situation with two particles moving with opposite velocities, for which the solution is plotted on Figure 1.
From this, we can plot on Figure 5 the corresponding regular field (u reg , p reg ) = (u, p) − (u sing , p sing ).
We observe that, within the gap, the regular velocity flow is almost constant and the maximum of pressure flow is much smaller than for the reference solution.This suggests that it will be much more easy to compute numerically this field, with not too many degrees of freedom.Using asymptotic results from [6] (recalled in Subsection 2.1.3),we will prove in Subsection 2.3 that this regular flow is indeed bounded independently of the size of the gap.This will allow to obtain convergence results independent of the distance in the numerical tests presented Section 3.
Remark 2.3 Note that contrary to the fields u in = w[q in ] defined in (9) or v in = w[p sing ] in (14), since we do not use the same pressure field for the tangential and normal terms, the velocity field u sing is no longer divergence free.Remark 2.5 A straightforward choice for the singular velocity field could have been, as for the pressure, to truncate u in for both the tangential and normal components: use u in = w[q in ]χ L in G L , extended by zero in F. This choice is of course different from the actual one since w[p sing ] // = w[q in χ L ] // = w[q in ] // χ L and unfortunately, it does not make the gradient of the singular pressure cancel in (5).As a consequence, the corresponding regular field (u reg , p reg ) is not bounded independently of the distance.The corresponding numerical error based on the computation of the remaining field (u reg , p reg ) provides better results than when computing directly (u, p).The error is less dependent on the distance but still, it is increasing when the distance decreases (see [10]).
Remark 2.6 Let us again emphasize the differences between the singular field u sing used in this paper and the asymptotic expansion v as = v in + v ext defined in [6].
• The field u sing , as well as v in , is designed to represent the singularity in the gap between the particles.They both vanish outside G L , their tangential components are equal while their normal components are different: • In [6] the inner velocity is corrected and the authors finally consider the asymptotic field v as = v in + v ext which fits the boundary condition on the moving particle but does not vanish anymore outside G L .In this paper we use u sing as singular velocity field, without any modification: it does not fit the boundary condition on the moving particle but vanishes outside G L (see again the singular field plot on Figure 4 for example).

Estimation for u reg
The regular field u reg is defined as u reg = u − u sing .Since v as = v in + v ext and u sing // = v in,// , we obtain The first three terms on the right-hand side are bounded from [6] (see equations ( 17), ( 18) and ( 19)).It remains to bound the last term.To do so, χ L being regular, we know that there exists C L > 0 such that, To bound these two terms, we first remark that these fields vanish outside G L and that, concerning the first term of the sum, we have As a consequence, in order to bound u sing ⊥ , it is sufficient to bound ∇w[q in ] ⊥ together with ∇ // w[q in ] // in L 2 (G L ).To do so, we recall that v in = w[p sing ] so that w[q in ] is constructed from q in as v in is constructed from p sing in [6].As a consequence, the computations of [6] to obtain (17) can be carried out the same way to obtain a similar upper bound, which provides a constant C independent of d such that . Moreover, if the function γ is radial one has, u reg ≤ C.

Numerical validation for pairs of particles
Due to the computational cost of three dimensional direct simulations, we propose here to validate the method using numerical simulations in dimension two.The results of the previous section can straightforwardly be extended to this situation.Note that, concerning short-range interaction and lubrication, dimension two is more challenging than dimension three.Indeed, the situation is more constrained since the fluid has one dimension less to escape the gap between the particles: the singularity induced by lubrication effects is more important.
After detailing the numerical scheme we use, we test the method for different configurations for the case of two (spherical or not) particles with general rigid velocities.Note that, since a finite-element discretization is used, the fluid domain is bounded, which is not the case in the theoretical results of the previous section.In fact, if the particles are apart from the boundaries, the previous results still hold.In case the particles are close to the boundaries, one should add the corresponding singularity to the computations.

Numerical scheme
We propose to solve (4) using a finite-element method.We consider an inf-sup stable pair of finite-element spaces for the Stokes problem V h ⊂ H 1 (F) and M h ⊂ L 2 (F).The singular field is detailed in Section 2.2 using a truncated asymptotic expansion (see equations (11) and (20)).Thus, since (u sing , p sing ) are known, we are looking to compute an approximation (u reg h , p reg h ) of the regular field solution to the following Stokes problem: where u sing h is the projection of u sing onto the space made of the traces on ∂F of functions of V h .Note that the terms ∆u sing , ∇p sing and div u sing have to be computed as precisely as possible and cannot be computed using u sing h .Then, the discrete approximation of the field (u, p) solution to problem (2) is given by, It is well-known that non-homogeneous Stokes problems as (4) or ( 22) have a unique solution if and only if a compatibility condition is satisfied by the right-hand side of the problem.These compatibility conditions are respectively, Thanks to divergence theorem and assumption (1) on u i the first compatibility condition is satisfied, and the continuous problem (4) has a unique solution.However, the second compatibility condition does not hold.Indeed, the volumic integral involves u sing while the surface integral involves u sing h .Moreover, note that at the implementation level, the integrals will be computed numerically, using quadrature formula, which will give rise to a second source of incompatibility.As a consequence, even if all these formula are consistent when h goes to zero, the compatibility condition is not verified at the discrete level and problem (22) may not have a solution for h > 0.
Thus, to remedy this problem, for a given ε > 0 we solve the following penalized Stokes problem, Integrating the divergence equation we now obtain, In the classical Stokes equation, the pressure field is defined up to an additive constant.The penalized Stokes problem allows us to fix the constant so that the mean value of the pressure satisfies the previous equation to compensate the fact that, from a numerical point of view, the compatibility condition is not satisfied.Thus, the problem (23) admits a unique solution.
Remark 3.1 Parameter ε has to be chosen sufficiently small so that the error induced in the divergence equation is smaller than the one due to the finite-element discretization.On the other hand, the finite-element discretization being fixed, letting ε go to zero in the numerical tests shall induce numerical errors since the discrete problem for ε = 0 has no solution.

Numerical simulations
We run in this section two tests involving a pair of particles.In the first one, we consider two spherical particles moving along the line of their centers (see Figure 6a) and we check that the error for the approximation of (u, p) in that case does not depend on the distance.Then, we show how the explicit expansion computed for a normal relative displacement of two spherical particles can be used to solve more general cases such as non-spherical particles undergoing general rigid motions, that is u i = V i + w i (x − x i ) ⊥ for i = 1, 2 (see Figure 6b).The numerical tests have been implemented with FreeFEM [15] and we choose V h = P 1 isoP 2 and M h = P 1 .Since we do not know the exact solution for this type of problem, we compare our results with a reference solution (u ref , p ref ) computed with a direct solver on a very fine mesh.In each case we choose the penalization parameter ε = 10 −8 .To compute the right-hand side of the problem, u sing is projected on the very fine mesh used to compute the reference solution.The linear system associated to problem ( 23) is solved with a direct solver.The boundary condition is projected on the space V h and is imposed by penalization.
For each configuration, the tests are run for five different distances between the particles: d = r 2 /10, d = r 2 /20, d = r 2 /30, d = r 2 /50 and d = r 2 /100.We plot the absolute error between the approximate solution and the reference solution (in H 1 0 -norm for each component of the velocity and in L 2 -norm for the pressure) for the first three distances.For the last two distances, the reference solution is not very accurate so that the errors are not relevant.The behaviour of the computed forces and torques for all the distances are also studied by comparing the values computed on the reference mesh and on coarser meshes.

Two spherical particles moving along the line of their centers
The configuration is the one given in Figure 6a with r 1 = 0.07, r 2 = 0.1, u 1 = e x , u 2 = −e x .In that case one has u as = u .
We study the convergence of the method when the size of the mesh h goes to zero.We compare the results with the ones obtained without any correction.We display on Figure 7 the fine mesh used to compute the reference solution (which is displayed on Figure 1), together with three of the meshes used for the computations.We observe in Figure 8 that without correction (dashed line) the error increases as the distance between the two particles decreases while, as expected, the error no longer depends on the distance when using the new method with correction (solid line).Moreover, the error between the approximate solution and the reference solution is better with correction than without correction.
Note that we plotted the absolute error which is proved to be bounded from our theoretical results.The relative error is much smaller and, using the correction, it decreases when the particles approach.Indeed, the norm of the solution increases when the distance goes to zero.We provide in Table 1 the values for these norms, which gives an idea of the orders of magnitudes.
In Figure 9, we plot the local error with and without correction.One can observe that this error is located in the gap between the particles and that the correction improves the results.Note that for each figure the Table 1: Norms of the reference solution color scale is different.Indeed, it is necessary to use the color scale corresponding to each figure rather than a common color scale to better observe, in each case, where the error is localized as well as its order of magnitude.Let us now consider the forces (see equation ( 3)) exerted on the particles in the longitudinal direction, that is We compare in Figure 10 their numerical approximation with the corresponding reference value (in black).We only consider the components according to e x since the component according to e y and the torque are null, due to the symmetry of the configuration.Note that due to the high variation of the normal forces for these small distances, we plot the results in logarithmic scale.Thus, since the force F x 1 is negative we actually plot −F x 1 .As expected, a symmetrical behaviour between the two particles is observed.In addition, we observe that the forces obtained with correction (solid lines) converge faster towards the reference forces than those obtained without correction (dashed lines).Note also that, the smaller is the distance, the greater is the improvement due to the new method.

Two non-spherical particles with general velocities: approximation of the singular field and the geometry
We consider now the case of non-spherical particles with a general rigid motion , w 1 = 3, w 2 = 5, r 2 = 1/20 and the radius of the osculating circle to particle 1 at the contact point is r 1 = 1/15.
To apply the theoretical results given in Section 2.2, we must be able to explicitly compute the singular field (u sing , p sing ) defined by (11) and (20).In Section 3.2.1 since we considered two spheres moving along the line of their centers, the singular field could easily be explicitly computed.In more general configurations, this computation could be trickier.Here, in order to go back to the previous computations, we choose to make approximations both on the geometry of the particles and on the singular boundary condition u as .The first approximation is a geometrical approximation: we approximate the non-spherical boundaries (∂B 1 in our case) by the osculating circle at the contact point (dashed lines in Figure 6b).This amounts to approaching the parametrization of the boundaries γ t and γ b in (11) and (20).
Then, we approximate the boundary conditions: we consider the velocity of the boundary at the contact point and use its projection along the normal direction instead of u as : Doing so, we take into account the part of u as leading to the most singular behavior of the flow.
To sum up, the singular field is the one that would be computed for two spheres approximating the particles (osculating circle at the contact point) and moving along the line of their centers with velocities equal to the projection along this direction of their initial velocities at the contact point.Then, explicit computations can be achieved as for the previous case.
Due to these approximations, the norm of (u reg , p reg ) now depends on the distance but in a less singular way that the original field (u, p).We check here that these simplifications are sufficient to catch the lubrication effects.As for the previous case we begin by plotting in Figure 11 the error between the approximate solution and the reference solution.We observe that even if the whole singularity is not taken into account in the singular field, the behavior for the velocity field (Figures 11b and 11c) is very similar to what was observed in the previous section (see Figure 8).Indeed, here again, the error with correction does not depend on the distance and the error is smaller than the error without correction.However, we notice that the error for the pressure in Figure 11a depends a little bit on the distance in the case with correction, but still much less than without correction.Note that the pressure being the most singular part of the solution, it is not surprising that the dependence on the distance is more perceptible for this term.
As for the spherical case, we compare in Figure 12 the forces numerically computed with the reference forces (in black).Since we consider general rigid motions we also plot the component according to e y and the torque, that is for i = 1, 2 We observe that for the component according to e x the force with correction approaches the reference force much better than the force without correction (see Figures 12a and 12d).As expected, these normal forces along e x are well corrected since the new boundary condition ũ as takes into account the normal motion of the contact points.For the component according to e y as well as for the computation of the torque, this is no longer the case but the behaviours with and without correction are similar.

Extension to 3 aligned particles and 2 singularities
We come back to the 3-dimensional problem of multi-particle systems (2) for which the solution is written St(u 1 , . . ., u i , . . ., u N ).Once again, the linearity of Stokes equation allows us to write this problem as a sum • Case 2 (see Figure 13b): In that case, two singularities have to be taken into account, due to the relative velocities between B 2 and B 1 and between B 2 and B 3 .With evident notations, we define L as the union of the gaps between the two pairs of particles.The singular field is defined as (u sing , p sing ) = (u sing,2,1 , p sing,2,1 ) + (u sing,2,3 , p sing,2,3 ).
Again, (u sing,2,i , p sing,2,i ) is the singular field computed in G 2,i L from ( 11) and (20) and extended to zero in F (of course, one has to use an appropriate change of coordinate before using formula (11) and (20): for each singularity, the origin of the coordinates is the point on ∂B 2 which achieves the minimum between B i and B = B 2 and e ⊥ is the unit vector from the origin to particle i).Again, the singular field vanishes outside G L .

Construction of the intermediate field v as
In order to prove that, using this singular field, the corresponding regular field u reg = u − u sing is bounded, we proceed as in Section 2. We construct an intermediate asymptotic field v as = v in + v ext , which is an extension to the one constructed in [6] for two particles.
The inner field v in is defined from the inner field for two particles as u sing has been defined from the singular field for two particles: with evident notations, Case 1 (see Figure 13a): in .Case 2 (see Figure 13b): L/2 in case 2. It now remains to define v ext so that the couple (v in , v ext ) shares the same properties as in Section 2. To do so, we recall that we need to construct v ext , bounded independently of the distance and such that v as = v in + v ext satisfies the boundary conditions of the initial problem.Following Section 2 and [6], v ext is defined as the solution to a Stokes problem in a suitable outer domain Ω .The main difficulty is to be able to choose Ω such that the boundary condition is corrected and v ext is bounded independently of the distance.To this end, as in Section 2, the domain Ω has to satisfy: • Ω is a smooth subdomain of F \ G L/4 (with evident notation); • Ω is independent of d; • Ω ⊂ F for any d ∈ (0, 1]; • ∂Ω contains the part of the boundary of ∂B on which v in = u .Supposing, as in Section 2, that for each singularity, the contact is "single, non-degenerate", such a domain exists for the two configurations studied.One of the choices for Ω is represented on Figure 13. From this construction, the computations and estimations in Section 2 can be extended to obtain the asymptotic result (19) and one gets the required estimate (21) on u reg .

Case of 3 non-spherical aligned particles with general velocities: numerical results
In this section, we consider three non-spherical particles with general rigid motions u i = V i + w i (x − x i ) ⊥ .As in Section 3.2.2,we approximate the geometry of the particles at the contact point using the osculating circles and the singular velocity is computed using the more singular part of the boundary condition: We choose V 1 = (1, 5), V 2 = (−1, −2), V 3 = (1.5, −5), w 1 = 3, w 2 = 5, w 3 = 10, r 2 = 1/20 and the radius of the osculating circles at the contact points are r 1 = 1/15 and r 3 = (1.3) 2 /20.The corresponding configuration is plotted on Figure 14.
The convergence of the method is illustrated on Figures 15 and 16 where we plot respectively the error versus h on the flow and the convergence versus h of the forces and torques on the three particles.The results are similar to the ones obtained in Section 3.2.2 for two particles: the error on the flow is better when using the singular/regular decomposition and do not depend on the distance for the velocity.It slightly depends on the distance for the pressure field but less than without correction, which again is due to the fact that the whole singularity is not taken into account.The forces along the horizontal direction are better approximated using the new method and the other components behave similarly with or without correction.To sum up, from a theoretical point of view, we propose in this article a new explicit expression to compute a singular field (u sing , p sing ) in three dimensions.We prove that the remaining field (u reg , p reg ) = (u, p) − (u sing , p sing ) is regular in the sense that it is bounded independently of the distance.We achieve this study for two particles or three aligned particles.The particles are not necessarily spherical (but need to be locally convex close to the contact point, see Remark 2.1) and undergo general rigid velocities.The method is easily transposed to the 2-dimensional case and the numerical behavior of the method for two or three aligned particles is investigated in 2 dimensions.We show that, with the new method, the error made on the computed fields does not depend on the distance anymore and that the forces are better approximated.
The extension of the method to more general configurations involving any number of non aligned particles is straightforward: mimicking the case of three aligned particles, one would like to sum up the different contributions of the singularities induced by the pairs of close particles.To illustrate this, let us consider three particles in a L-shape configuration (see Figure 17).Suppose that the radius of the particles are r 1 = 2, r 2 = 1, and r 3 = 1.5.Particles 1 and 2 are aligned on the horizontal axis and the third particle is above the second one.Suppose that particle 2 is at rest and that both particles 1 and 3 move towards it.In that case, one has two singularities to sum: between particles 1 and 2 and between particles 2 and 3.As in Section 4.1 (dealing with three aligned particles), to prove that the remaining field is regular, it is sufficient to construct an adapted Ω for the two following cases: B = B 2 and B = B 1 .We show in Figure 17 the corresponding constructions, together with the reference solution computed on a very fine mesh for distances between particles 1 and 2 and particles 2 and 3 equal to r 1 /10.The demonstration of the regularity of the remaining field can then be extended which is confirmed by the fact that the reference solution in Figure 17c behaves as the sum of two singularities concentrated in the gaps between the pairs of particles 1 and 2 and particles 2 and 3.As a consequence, the method can be extended in that case and convergence results should be similar to the ones obtained in the case of aligned particles.As long as it is possible to construct a domain Ω verifying the assumptions recalled in Subsection 4.1.2,one can choose the singular field as the sum of the pairwise singular fields and prove that the remaining part of the field is regular.We provide on Figure 18 an example in dimension 2 for which this construction is not possible.This example exhibits a closed subdomain in which the fluid is confined when the distance between the particles goes to zero.It is made of three particles located at the apexes of an equilateral triangle.All the particles have the same radius r and the distance between the particles is equal to r/10.In Figure 18a we observe that it is not possible to construct a domain Ω (independent of the distance and included in F for any distance) whose boundary takes into account the blue part of ∂B 2 when the distance goes to zero whereas we should correct the boundary condition on this part of ∂B 2 .We plot in Figure 18b the pressure field when one of the particles (bottom right) is at rest and the two other ones move towards it.We can observe that the solution does not behave as the sum of the singularities in the gaps: a new kind of singularity appears in the confined zone.Note that such a situation exhibiting a confined zone of fluid does not exist in dimension 3.As a consequence, one can always construct an adapted Ω and the method should work for any configuration of particles.The numerical behavior of the method for more general configurations in dimension 3 will be investigated in a forthcoming work.
The behavior of the method for higher numbers of particles will also be investigated.Note that the contribution to the source term of each pair of close particles in (4) is independent so that its implementation can be parallelized.Concerning the behavior of the error, note that we proved that adding all the singularities in the source term leads to a regular remaining field: the solution to the corresponding regular problem is bounded independently of the distance.As a consequence, the number of particles being given, the numerical error itself is bounded independently of the distance, allowing to use coarse meshes when the particles get close one to another.
To finish, let us comment the fact that, as already mentioned, the particles have to be locally convex close to the contact point.The method being local, no hypothesis is made on the global convexity of the particles

Figure 1 :
Figure 1: Reference flow for two particles with opposite velocities along the horizontal direction.

Remark 2 . 4
The resulting velocity field u sing ∈ C ∞ (F) has support in G L .It vanishes on the boundary of the obstacle and u sing = u on ∂B ∩ G L/2 .Note that u sing does not satisfy the boundary condition on the whole boundary ∂B of the moving particle.
(a) Singular pressure flow (b) Singular horizontal velocity (c) Singular vertical velocity

Figure 4 :
Figure 4: Singular flow for two particles with opposite velocities along the horizontal direction.

( a )
Regular pressure flow (b) Regular horizontal velocity (c) Regular vertical velocity

Figure 5 :
Figure 5: Regular flow for two particles with opposite velocities along the horizontal direction.

Figure 7 :
Figure 7: Different meshes used for the computation

for uy d = r2 10 with correction d = r2 20 with correction d = r2 30 with correction d = r2 10 without correction d = r2 20 without correction d = r2 30 without correctionFigure 8 :
Figure 8: Absolute error in function of h for 3 different distances for two spherical particles

Figure 9 :
Figure 9: Local error for two particles with opposite velocities along the horizontal direction, d = r/16, h = 0.012.

Figure 10 :
Figure 10: Comparison of the forces in function of h for two spherical particles

− p ref L 2 (a) L 2 -error for p 10 Figure 11 :
Figure 11: Absolute error in function of h for 3 different distances in the general case

Figure 14 :
Figure 14: Configuration for the 3 non-spherical particles.Particles, together with their osculating circles at the contact point.

Figure 15 :
Figure 15: Absolute error in function of h for 3 different distances for 3 aligned particles Figure 17: Three particles in a L-shape configuration.
Evolution of the configuration when the distance goes to zero (b) Pressure field solution to the problem

Figure 18 :
Figure 18: A limitation in the two-dimensional case: three particles in a T-shape configuration, existence of a confined zone of fluid.