Numerical Determination of Truncation Orders in the Correction Method for Stokes Equations

In [15], the authors propose an accurate method, namely the correction method, for computing hydrodynamic interaction between very closed spherical particles in a Stokes uid. The accuracy of this method depends on two truncation parameters for approximating the Neumann to Dirichlet matrix and the velocity correction respectively. In this paper, we establish a numerical determination to estimate these parameters. We perform some numerical experiments to present our method.


Introduction
In recent papers, many mathematicians have studied the hydro dynamic interactions between close particles immersed in a Stokes uid to simulate the motions of nano scale swimmers robots which is designed from nano-sized medical devices (see [2], [3] and references therein).One of the main diculties is to compute the total forces and torques exerted by the particles.Our work is motivated by this problem.
We consider N non intersecting particles immersed in a viscous uid.For simplicity, the particles are identical balls B 1 , B 2 ,..., B N with radius 1 and centers z 1 , z 2 ,..., z N ∈ R 3 , respectively.We assume that the closed balls B i do not intersect and that the uid lls the rest of the space.The uid occupies the domain Ω := R 3 \ N i=1 B i .We assume moreover that the uid inertia eects are negligible compared to the viscosity (i.e. the Reynolds number is very small Re 1) so that the velocity u and the pressure p solve the stationary Stokes equations in the uid domain, where σ = ∇u+∇u T −pId is the stress tensor in the uid and Id is the identity matrix.On the surfaces of the particles, we consider a no-slip boundary condition, u = u i on ∂B i , i = 1, 2, ..., N.
where the velocity u i corresponds to a rigid displacement.It is characterized by the velocity U i at the center z i of the ball B i and by the angular velocity ω i (U i , ω i ∈ R 3 ), We are interested in solutions which decay at innity, i.e., which fullls u(x) → 0, p(x) → 0, as |x| → ∞.The existence and uniqueness of a solution to (1) is classical in the Hilbert space (see [10] and [17]).
The surface density of force exerted on the uid at some point x of the surface ∂B i is given by f i (x) = ∇u + ∇u T − pId .ni , (3) where n i denotes the exterior normal on the surface of the particle B i .The total force and total torque exerted by the particle B i on the uid are given by the following formulas, We recall that the solution to the transmission problem u is given by the convolution of the surface force density f i with the Green tensor for the Stokes equation as follows This explicit formula gives the velocity eld everywhere as soon as the force densities f i are known.However, the data of the problem are the velocity eld u i , not these force densities.
We are led to consider the following Dirichlet to Neumann operator This operator is positive and symmetric, its inverse is the corresponding Neumann to Dirichlet operator" N D. In the initial problem, we only need to compute approximations of DN when (u | ∂B i ) 1≤i≤N is a nite sequence of rigid motion.Moreover, we do not need a complete description of (f i ) 1≤i≤N but only the projections of these densities.In short, we only need a projection of the operator DN on a nite dimensional space of dimensions 6N .This projection is called the friction operator as follows Unfortunately, we do not have a nice explicit expression for DN such as (4).To compute accurate approximations of F starting from (4), the direct method consists in 1/ approximating N D by a Galerkin method, 2/ inverse this approximate operator, 3/ project this inverse on the space of rigid motions.For a xed position of the particles, the direct method has a very good behavior as we send the truncating order to innity.On the contrary, if we consider a sequence of congurations with at least two particles getting closer and closer and with different prescribed velocities, the distributions of forces concentrate near the contact points.In this case, the convergence of the direct method degenerates.
A large amount of research has been studied to develop numerical tools to approximate the friction operator, such as [5], [6], [8], [9], [13].Recently, in [15], the authors developed a method which is called the correction method for computing very accurate numerical solutions.This method depends on two truncating parameters for approximating the Dirichlet to Neumann operator DN and computing the correction velocities respectively.In the correction method, these parameters are chosen to obtain a good approximation of friction operator.The main goal of this paper is to present a numerical method to determine truncating parameters for the correction method.
The rest of this paper is organized as follows.
In the next section, we describe the main idea of the correction method which is based on the singular-regular decomposition.We also present the interpolation method for computing the singular elds and the computation of correction velocities for two close particles.In Section 3. , we propose some numerical computations to estimate the truncation orders using in the correction method.In Section 4. , numerical experiments are performed for four particles problem.
Finally, we give several conclusions and perspectives in the last section.

2.
The correction method Our method consists in taking advantage of the linearity of the Stokes equations for rewriting the elds (u, p) as a superposition where each couple (u 0 , p 0 ) solves the Stokes equations in Ω and (u c , p c ) solves the Stokes equations in the ctitious uid domain: The couple (u c , p c ) handle the large variations of (u, p) localized in the small gap between B i and B j which are due to the dierence between the prescribed velocities on ∂B i and ∂B j .Precisely, for c = (i, j) ∈ P, we introduce the velocity eld which vanishes if and only if the solid B i ∪B j follows a rigid motion.The singular eld (u c , p c ) are dened as the unique solution of the problem By linearity, the remaining part (u 0 , p 0 ) solves the Stokes problem in Ω.The boundary conditions u 0 for this problem are set so that the total velocity eld satises the boundary conditions u i specied in the original problem where, At the end we aggregate the dierent contributions.With obvious notation, for k = 1, 2, ..., N , Notice that the singular solution (u c , p c ) associated to a pair of close particles c = (i, j) ∈ P do not contribute to the forces and torques exerted by the surface of a third particle B k , k / ∈ {i, j} : Similarly, using the Levi-Civita antisymmetric symbol ε αβγ and Einstein summation convention on greek indices, we compute, As a consequence, the total force and torque exerted by the particle B k on the uid are given by The advantage of decomposing the solution resides in the possibility of using dierent methods for solving problems ( 5) and (6).The singular parts are solution of the Stokes equations (5) around only two solid particles.We will approximate these singular parts by interpolating in pre-computed tables.The remaining parts solves the Stokes equations ( 6) in the original domain but with modied boundary conditions which do not necessarily correspond to rigid motions of the particles.The remaining regular part may be approximated by using any standard numerical method.
Let us rst consider problem (5).For c = (i, j) ∈ P, by changing coordinates, we may assume that z i = −(1 + d c /2)e z and z j = (1 + d c /2)e z .In the new coordinates, the velocity w c uniquely decomposes as where e 1 and e 2 are two unit vectors orthogonal to e z .Hence, the solution of ( 5) can be decomposed as where, for Z = A, B, B , or C, the couple (u Z , p Z ) solves the Stokes equations in the domain where B dc ± denotes the solid sphere with unit radius and center ±(1 + d c /2)e z .The dierence between these problems comes from the specic boundary conditions, where w Z are dened as follows, for x ∈ ∂B dc When solving independently the second or the third problem, we may rotate the frame so that e 1 or e 2 coincide with e x .We end with four family of problems only depending on the distance d c .More precisely, in view of (8), we need approximations of T Z (d c ) := Using the symmetries of the problems, the corresponding total forces and torques on ∂B dc − are deduced from the former.
For the computation of the boundary conditions (7) satised by the remaining regular part (u 0 , p 0 ), we also need approximations of v Z (x, d c ) := u Z (x), for x ∈ Ω dc and Z = A, B, B , C.
In the next section, we describe a procedure for computing these quantities.The method is based on known asymptotic as d c → 0, direct computations and interpolation in the parameter d c .
Let us now consider problem (6).
It is of the same nature as the original problem: solve the Stokes equations in the uid domain surrounding the particles.The new problem looks even more complex since we have substituted the function w 0 for the simple rigid motions u i that can be described with 6N parameters.However, by construction, (u 0 , p 0 ) is a very regular vector eld, even in the limit of touching particles.
As a consequence, applying standard numerical methods to problem (6), we can compute approximations of (u 0 , p 0 ) with an accuracy that does not depend on the distance d c between close particles.
In the next section, we describe a procedure for computing the approximations of the singular part and the boundary condition for the regular part.The main idea is to interpolate the needed quantities into a grid of known values which has been computed once for all during a pre-processing step.

The interpolation method for computing the singular elds
As explained in the discussions at the end of Section 2.1., for each c = (i, j) ∈ P, the singular part (u c , p c ) can be decomposed as a combination of four parts (u Z , p Z ) which are solutions of four family of problems only depending on the where Ω dc and w Z are given by ( 10) and ( 11) respectively.Recall that the uid domain Ω dc only c 2017 Journal of Advanced Engineering and Computation (JAEC) depends on the distance d c .We need to compute approximations of f Z and T Z given by ( 12), (13).
Our method is based on asymptotic formulas for the total force and torque at small distance, direct computations and interpolation in the parameter d c .
In a pre-processing step, we decompose ( fZ ) k and (w Z ) k , for k = B dc − , B dc + , in the basis of vector spherical harmonics (see [16]) as follows, By truncating the above series up to order M max , with M max large, the discrete Neumann to Dirichlet matrix N D Z,dis. is computed as described in [18].Then we compute accurate approximations of the surface force density f dis.
Z by solving the linear problem, where W Z = w Z i,α k=i,j,α=1,...,Mmax . By this direct method, we may compute f dis.For instance, let us consider the rst problem Z = A. We are interested in the total force and torque exerted by the rst particles B dc − .In this case, from the symmetries and the asymptotic formulas given by (see [7]) We guess that F A (d c ) expands as The constants C 1 , C These values are obtained by the direct method with a very large L max .
This ends the preprocessing step.
Eventually, when needed, we use the cubic spline interpolation method to estimate R A (d c ) for any non-tabulated distance d c ∈ (0, δ) from the tabulated values.In Table 1, we show the result of some numerical tests realized in order to estimate the error due to the interpolation method.

Computation of correction velocities
In this section, we present the interpolation method to compute the coecients of the correction velocities.
We consider again the problem (14).Let B R be the ball of radius R = 3 centered at the origin of the coordinate system.This ball contains the two balls B dc − and B dc + .We want to determine the velocity U Z (r, d c ) for r ∈ R 3 \ B R .
We rst compute the force densities on the boundary of B dc + and B dc − using the direct method with a large truncating order.Then, we can deduce the velocity eld U Z everywhere using the explicit formula (4) and vector spherical harmonics (see [18]).On the other hand, we know that the velocity eld in R 3 \ B R reads We then only have to tabulate the coecients g T Z,l,m , g I Z,l,m , g N Z,l,m .These coecients are ob- tained by projecting U Z (•, d c ) on the basis of rescaled vector spherical harmonics on ∂B R .In a last step, we use (15) to obtain the corresponding coecients in the vector spherical harmonic basis on ∂B(0, 1).
In practice, the series ( 15) is truncated at some order Lmax .We call Lmax the correction truncation order.Notice that this truncation order may be dierent than L max dened in [18].The choice of Lmax will be discussed in the next section.

Numerical determination of the truncation orders
In the correction method, when we approximate the correction w 0 determined by (7) and the Neumann to Dirichlet matrix DN , we have to choose two truncating parameters: L max for approximating the Neumann to Dirichlet matrix and Lmax for approximating the velocity cor- rections.These quantities prescribe the number of vector spherical harmonics used for the discretization.The natural question is how can we choose these parameters such that the solution has a given accuracy?How do they depend on the distances between the particles?In this section, we present a numerical estimation of these parameters.

Correction truncation order
Let us consider the problem (1) with three unit balls.We assume that their centers lie on the vertical axis with corresponding coordinates z 1 = (0, 0, 0), z 2 = (0, 0, 2 + d), and z 3 = (0, 0, 4 + d + D).We assume moreover that the two rst balls translate along the vertical axis with opposite velocities and that the third c 2017 Journal of Advanced Engineering and Computation (JAEC)

Truncation order for solving the problem
We now consider the same three-sphere conguration as in the previous section.For computational time problem, we could not calculate f dis.L, Lmax for very large values of L. The error on the surface force density is estimated by the dierence between two consecutive values of L with Lmax determined in the previous section.
For every L ≥ 1, we dene The truncation order L max is chosen as follows, for a given small real number ε > 0, In fact, the truncation order L max can be also estimated with another denition of the density error, where L∞ is very large.The two errors are very close in the numerical computation.Hence, it is more convenient to use the rst denition (16).
• The rst case: D = D 0 > δ, In this case, the truncation order mainly depends on the distance d.We conclude that for isolated pairs of particles D ≥ δ, we see that the critical truncation level is a monotonic increasing function of d.
We can choose  We made these tests with d varying from 0.01 to 0.5.Even for d = 0.01, the truncation order L max = 15 lead to an error smaller than ε = 10 −6 (see Figure 5).
• The second case: D = D 0 < δ, In this case, the optimal truncation order depends on both d and D. This truncation order tends to innity as both d and D go to 0.
In practice, we see on the graphic that we can

Conclusions and perspectives
In conclusion, we have presented an accurate method for the computations of hydrodynamic forces between spherical particles suspended in a Stokes uid.The main improvement of this new method compared with the Stokesian Dynamics is that the inuence of the singular force densities between two closed particles on the neighboring particles is also computed.For this reason, the computational cost for this method is larger.The main part of the computational time is due to the computation of the correction velocities and their projection on the vector spherical harmonics basis.On the other hand, these computations are independent from one sphere to another and could be easily parallelized.This should solve the main drawback of the method.Moreover, we proposed some numerical determination for some parameters using in this method.
In our research, we only consider spherical particles.The main advantage of this shape is that the computation can be based on the vector spherical harmonics basis.The methods generalize to arbitrary smooth particles.In this case, we should use a boundary nite element method instead of the decomposition in vector spherical harmonics.

Z
as a function of d c for a nite number of distances, say d c ∈ D dis := {d 0 , λd 0 , λ 2 d 0 , ...} for some small d 0 and some λ > 1. Combining the explicit asymptotic formula of the force density with this discrete set of accurately computed values, we obtain approximations of f dis.Z (d c ) by interpolation for every 0 < d c < δ.

Finally, using a polynomial 10 Fig. 2 :
Fig. 2: The absolute errors of interpolation correspond to the coecient g N Z,5,0 (dc) (right) in function of dc in the case Z = A.
choose L max as an ane function of log 10 D and log 10 d in the region of [L max,opt ≥ 40] (see Figure 6).c 2017 Journal of Advanced Engineering and Computation (JAEC) very fast and requires a small level of truncation order to get an accurate result.
B i and B j , d ij = |z i − z j | − 2.The set of pairs of close particles is dened as ticles c .The Figure1shows the behavior of the rest term R A (d c ). Fig. 1: The term R A (dc) in a function of dc.
2 , C 3 and C 4 are then determined by using a least square approximation based on highly accurate numerical simulations performed for a small number of small values of d Tab.1: The absolute errors of interpolation.R A (d c ) for d c ranging in a nite subset of (0, δ).