Improved procedure for the computation of Lamb’s coefficients in the physalis method for particle simulation
Abstract
The physalis method is suitable for the simulation of flows with suspended spherical particles. It differs from standard immersed boundary methods due to the use of a local spectral representation of the solution in the neighborhood of each particle, which is used to bridge the gap between the particle surface and the underlying fixed Cartesian grid. This analytic solution involves coefficients which are determined by matching with the finitedifference solution farther away from the particle. In the original implementation of the method this step was executed by solving an overdetermined linear system via the singularvalue decomposition. Here a more efficient method to achieve the same end is described. The basic idea is to use scalar products of the finitedifference solutions with spherical harmonic functions taken over a spherical surface concentric with the particle. The new approach is tested on a number of examples and is found to posses a comparable accuracy to the original one, but to be significantly faster and to require less memory. An unusual test case that we describe demonstrates the accuracy with which the method conserves the fluid angular momentum in the case of a rotating particle.
keywords:
1 Introduction
Sediment transport, settling suspensions, fluidized beds and porous media are but a few examples of fluid flows with fixed or moving particles. As these examples show, such flows are frequently encountered and of great importance in many problems of scientific and technological relevance. From a computational perspective, they pose a major challenge as the particles constitute a very complex and frequently time dependent internal boundary which needs to be adequately resolved for a bona fide simulation.
Until very recently the only practical approach was the EulerianLagrangian pointparticle model, in which the particles are assimilated to mathematical points contributing extra forces to the fluid, while their motion is found by integrating in time a dynamic equation in which the fluid forces are modelled rather than deduced from first principles. The last decade has seen the development of methods capable of dealing with the finite extent of the particles. The most successful such method so far is the one developed by Uhlmann Uhlmann:2005gf (); Uhlmann:2008il (), which has been also used by Lucci et al. Lucci:2010vd ().
An alternative approach, dubbed physalis, first described in Prosperetti:2001fj (); Takagi:2003hk () and, more fully, in Zhang:2005ty (), differs from standard immersed boundary methods in its reliance on a local spectral representation of the solution in the neighborhood of each particle; a brief description is provided in section 2. Advantages of this method are the exact satisfaction of the noslip condition at the particle surface, the great simplification of the calculation of the hydrodynamic forces and couples, and the avoidance of the complex issues arising from the lack of geometrical conformity between the curved particle boundary and the underlying fixed Cartesian grid. Furthermore, the use of a local spectral representation of the solution permits one to describe the effect of each particle with fewer degrees of freedom than conventional finitedifferencebased methods. As a consequence, the grid resolution can be kept relatively low without compromising the accuracy of the solution.
physalis has been thoroughly validated in the original papers and it has proven its value in a number of applications: the sedimentation of 1,024 particles Zhang:2006wk (), the turbulent flow around a fixed particle Naso:2010ti (), the flow induced by a rotating particle Liu:2010tu () and over a porous medium Zhang:2009ky (); LIU:2011hn (). An interesting variant which uses artificial compressibility to generate a timestepping pressure equation has been described and used by Perrin and Hu Perrin:2006cl (); Perrin:2008dh ().
The basic idea of physalis is to use an analytic local solution of the modified (Navier) Stokes equation as a bridge between the particle boundary and the finitedifference solution. The local analytical solution contains undetermined coefficients which are found iteratively by matching the finitedifference solution with the local one.
The significant improvement described in this paper concerns the calculation of the expansion coefficients, which was carried out using the singularvalue decomposition in the original implementation. We show how this step can be executed by taking suitable scalar products with a significant saving in execution time. The purpose of this paper is to describe the implementation of this new method and its verification.
2 General framework
We briefly review here the physalis method for particle simulation, as originally introduced for cylindrical particles Prosperetti:2001fj (); Takagi:2003hk (), and later extended to spherical particles Zhang:2005ty ().
The objective is to simulate, on a fixed Cartesian grid, incompressible fluid flow around spherical particles, while obeying the noslip and nopenetration velocity boundaryconditions on the particle surface. To this end it is necessary to solve the incompressible NavierStokes equations,
(1)  
(2) 
where , , , and respectively denote the velocity, density, pressure, and gravity fields.
As with other immersedboundary methods, since the particles do not conform to the computational grid, the boundary conditions imposed by the particles must be transferred to the adjacent nodes. The physalis method enables this transferral to be done in a spectrally accurate manner.
To proceed, we transform the above equations to a frame in which the particle is at rest. (In the case of more than one particle this step is carried for each particle.) The fluid velocity in the particle restframe is related to , the velocity in the laboratory frame, by
(3) 
where and respectively represent the linear and angular velocities of the particle, and is the position with respect to the particle center. The momentum equation in the particle rest frame is given by
(4) 
and that of continuity is given by
(5) 
Now let
(6) 
where , and represents the particle radius. Then, if we introduce the modified velocity and pressure fields
(7) 
we can write the restframe momentum and continuity equations respectively as
(8)  
(9) 
Since the term within the square brackets vanishes at the particle surface, it will be small in its immediate vicinity. Therefore, in this region, and will approximately satisfy the Stokes equations, given by
(10) 
for which Lamb Lamb:1932wz () developed an exact general solution for flow past a sphere. It should be stressed that the method is intended for application to finiteReynoldsnumber flows. The procedure may be seen as a linearization of the NavierStokes equations about the solidbody motion which the fluid in contact with the particle surface executes at any Reynolds number by virtue of the noslip condition.
The flow field is expanded in terms of vector spherical harmonics (allied to the standard scalar spherical harmonics), which form a complete basis for squareintegrable vector functions on the sphere Prosperetti:2011tu (). The velocity is given by Lamb:1932wz (); Kim:2005wv ()
(11) 
and the pressure by
(12) 
where the , , and are solid harmonics of order :
(13) 
Here are undetermined dimensionless coefficients and the functions are normalized spherical harmonics which are shown in detail in the Appendix. The regular harmonics, corresponding to nonnegative , represent a general incident flow, while the singular harmonics (negative ) represent the disturbance field due to the particle. Although the coefficients are in general complex, the relations (in which the overline denotes the complex conjugate) obviate the need for dealing with complex numbers and the entire calculation can be carried out in the real domain.
The central observation at the root of the physalis method lies in the recognition that, in the immediate neighborhood of the particle, the velocity and pressure fields must be describable as in (11) and (12). Thus, provided the coefficients are known, (11) and (12) can be used to “transfer” the boundary conditions from the particle surface to a cage of adjacent grid nodes enclosing the particle. After this step, the computation is carried out only on the nodes of the fixed Cartesian grid and the complexity deriving from the geometric mismatch between the particle boundary and the grid is eliminated.
With this background, we can now describe the algorithm:

Use these provisional values to generate new velocity boundary conditions at the cage nodes;

Solve again the NavierStokes equations on the finitedifference grid using these boundary conditions;

Repeat to convergence.
It is important to note that the coefficients appearing in (13) embody the physics of the interaction of the particle with the fluid and contain therefore important information. For example, the hydrodynamic force, , and couple, , acting on the particle can be expressed in terms of the coefficients as Zhang:2005ty ():
(14)  
(15) 
where is the volume of the particle, and we have defined , , and , with superscripts and denoting the real and imaginary parts. Thus, unlike other immersed boundary methods, once the coefficients are known the force and couple are known as well with no need for additional calculations.
As discussed in the following section, the specific contribution of the present paper concerns the first step, that of deducing from a given velocity and pressure field the values of the coefficients.
3 Computing the expansion coefficients
Away from the particle(s), the NavierStokes equations (1) are solved on a staggered grid by a standard finitedifference method as described in Zhang:2005ty (). This reference also explains in detail the manner in which the cage of Cartesian nodes enclosing the particle is constructed. Let , with be the cage nodes. In the original algorithm, the coefficients are calculated by requiring that the pressure field given by (12) and the vorticity field obtained from (11) reproduce the finitedifference results. In other words, one forms a linear system involving the coefficients as unknowns, e.g.
(16) 
in which , the lefthand side is as given in (12) and contains the unknown coefficients, while the righthand side is the finitedifference result at the cage nodes. The vorticity field is treated similarly. As explained before, the coefficients thus obtained are then substituted into (11) to generate velocity boundary conditions at the cage nodes.
In this way, four scalar equations are obtained for each cell comprising the cage. Taken together over the cage, this set of equations forms an overdetermined linear system for the unknown coefficients which is solved by means of the singular value decomposition (SVD) algorithm, i.e., in essence, in a leastsquares manner. Although the method is stable, the typically large systems that result (e.g., a matrix of for a gridresolution of 8 cells per radius, and with the analytical expressions truncated at ) prompted us to search for a more efficient procedure.
There is, however, an alternate and perhaps more natural way to obtain the coefficients which relies on the orthogonality property of the spherical harmonics, for example
(17) 
In the spirit of this approach, (16) is replaced by
(18) 
or, upon substituting (12) for and using (17) and recalling that ,
(19) 
It is shown in the Appendix that the coefficients of negative order can be expressed in terms of those of positive order to find
(20) 
This property of the spherical harmonics to single out individual coefficients is analogous to the Fourier transform and, as discussed in section 4, there exists a “Fast Spherical Transform” (FST) Healy:2003wn (), incorporating and scaling as the Fast Fourier Transform. As will be shown below, this procedure proves considerably faster than the SVD approach of Zhang:2005ty (), and comparably accurate.
The other coefficients can be determined by a similar procedure. However there are different possible choices of field variables to use for this purpose, which we now discuss.
In the derivations we will use two orthogonality properties of the vector spherical harmonics similar to (17), namely (see e.g. Ref. Prosperetti:2011tu ())
(21) 
(22) 
It is found that, whatever field one may consider, and appear coupled, while appears separately.
3.1 Computing and
In general it is desirable to choose pressure as one of the field variables used to calculate the coefficients and , since this field is smoother than the velocity field and, additionally, does not vanish near the particle. Thus, (20) is one of the equations that we use to determine the coefficients.
A combination of any two out of three other scalar products, involving the radial component of velocity ; the transversal components of velocity ; or the transversal components of vorticity may be used to solve for and . We will present each of these options below, and, in section 5, demonstrate their respective performance.

Transversal components of vorticity: Upon taking the curl of the velocity given by (11) and projecting on the transversal direction we find
(25) where for brevity we have defined the operator as, e.g.,
(26) Upon taking the scalar product of this equation with and using (22) we obtain
(27) where again (48) and (49) have been used to eliminate the coefficients with negative index.
Having in hand equation (20), along with either (23), (24), or (27), we can now solve for the two unknowns, and .
3.2 Computing
Three scalar products can be formed that involve :
4 Implementation
The secondorder accurate discretization of the NavierStokes equations, along with the construction of the cage and the manner in which the particle boundary conditions are applied at the cage surface, are described in detail in Zhang:2005ty (). We will therefore focus our discussion here on those aspects of the present methodology that diverge from the original version, namely, the calculation of the expansion coefficients.
The scalar products appearing in (18) and the following equations amount to integrations carried out over a spherical surface concentric with and enclosing the particle. In general, therefore, the quadrature nodes do not coincide with the finitedifference grid points. We thus require a procedure for approximating the values of , , and at the quadrature nodes on the basis of those on the finitedifference grid. To this end we adopt two separate approaches: interpolation, and local linearization. For the reasons explained below, we limit the latter to those scalar products involving and in relations (20), (24), and (28).
4.1 Interpolation
Since the spatial and temporal discretization is secondorder accurate, an interpolation of the same order or higher is appropriate. We experimented with both cubic splines and quadratic Lagrangepolynomials, both implemented in the GNU Scientific Library Anonymous:wz (). These proved comparably accurate, but the Lagrange interpolation is considerably more expedient and was hence adopted.
Given the point at which the value of e.g. the pressure is required, the interpolation proceeds as follows:

Let be the cell containing ;

Centered on, and including this cell, construct a cube of 27 cells having coordinates , with , at which is known;

Arbitrarily choosing to interpolate first on the  planes of the cube, we interpolate in for each and , and evaluate the results at , yielding ;

Interpolate in for each , and evaluate at , yielding ;

A final interpolation in is performed and evaluated at , resulting in the desired approximation of .
This threedimensional, secondorder accurate interpolation therefore requires 13 onedimensional interpolations in total.
4.2 Local linearization
For pressure, we write
(32) 
where is defined as above, and is the cellcenter which, given the staggeredgrid formulation, is where resides. The gradient is approximated with secondorder accurate (at ) central differences, using adjoining cells.
For velocity, the corresponding expression is
(33) 
Since the velocity is defined at face centers, we approximate , to secondorder accuracy, with the average across the cell. The gradients are calculated with central differences. The diagonal components reside at the cellcenter. As the remaining components reside at edgecenters, we approximate them to secondorder accuracy with the average across the cell.
While the quadratic interpolation has higher accuracy than the linearization, the latter has some advantages that make it worth exploring. Firstly, note that the estimate is divergencefree to the same accuracy as the velocity field at the nodes, since
(34) 
while the same is not necessarily true of the interpolation. Secondly, the operation count involved with the linearization is lower than for the interpolation. Lastly, the spatial extent of the region involved with the computation is reduced, as only the six adjacent, facesharing cells need be considered.
The same approach applied to the vorticity would only give a zeroorder accuracy. An effort to improve this low accuracy, on the other hand, would require a comparatively much larger set of cells. This is the reason why we did not pursue this method for the vorticity options.
With these procedures for sampling the flow field at the quadrature points on the sphere, the next task is to evaluate the scalar products. The present integrands (e.g., in equation (17)) are oscillatory and increasingly so with higher . Particular care must therefore be exercised in distributing the quadrature points in order to maintain stability. Given the homogeneity of the azimuthal () direction, a uniform grid , is appropriate, where and is the number of quadrature points in . The natural distribution in the polar () direction, spanned by the associated Legendre polynomials (cf. equations (46) and (47)), is , where has the same range as .
The scalar product on the sphere can essentially be considered as the product of discrete Fourier and Legendre transforms in the azimuthal and polar directions, respectively. Using the quadrature points described above, the cost of a direct numerical evaluation of these transforms scales as . Faster alternatives exist, however. The Fast Fourier Transform can naturally be used in the azimuthal direction, and there additionally exists a Fast Legendre Transform, due to Healy:2003wn (). The combination of these fast algorithms reduces the cost scaling to . The authors of Healy:2003wn () have released a Cimplementation of their method, the SpharmonicKit^{1}^{1}1See http://www.cs.dartmouth.edu/geelong/sphere/, which we embed in our code.
Lastly, we note that the Fast Spherical Transform of Healy:2003wn () includes only the spherical harmonics , and not gradients thereof. To compute the scalar products involving (e.g., equation (28)), we express these in terms of the :
(35)  
(36) 
It should be noted that, in spite of the explicit appearance of the imaginary unit, the actual computations are carried out in the real domain only.
5 Results
In this section we contrast the various aspects of the singularvalue decomposition (SVD) and the scalar product (SP) version of physalis. In particular we examine relative speed and accuracy, and address the question of which SP is best suited.
5.1 Comparing the speed of the SVD and SP
As mentioned in section 3, Zhang:2005ty () compute the expansion coefficients by matching the vorticity and pressure, based respectively on the analytical expansions in equations (11) and (12), to those from the finitedifference field, resulting in a linear system, . is a rectangular matrix of size , where and , where is the degree at which the summations in equations (11) and (12) are truncated; is the vector containing the unknown expansion coefficients; and contains the values from the finitedifference field. and respectively represent the number of vorticity and pressure nodes comprising the cage.
The commonly used SVD algorithm can be shown Lawson:1974wg () to scale asymptotically as . In contrast, the cost of the FST scales as , where is the number of samples. If the flow variables , , and are bandlimited so that , then will suffice to calculate the coefficients without error (this is analogous to the NyquistShannon sampling theorem Shannon:1949fb ()), assuming that the samples are errorfree Healy:2003wn (). In the general case we do not know and, therefore, to prevent aliasing, take . In terms of , the cost of the FST therefore scales similarly.
As shown in Zhang:2005ty (), suffices to adequately describe flows with particle Reynolds numbers up to about 100, with a grid resolution of , where is the grid spacing. As such, the asymptotic scalings above are not as relevant as empirical measurements of the cost associated with a single calculation of the coefficients. Such measurements are shown in figure 1, both as a function of (with fixed ) and of the grid resolution (with fixed ). The SP approach proves over two orders of magnitude faster than the SVD at all shown. We also note that the cost of both the interpolation and the linearization, described in the previous section, is independent of the grid resolution as the number of quadrature points is only determined by . This is apparent from figure 1, right, while the cost of the SVD grows strongly with .
While the SP has been shown to be significantly faster than the SVD, there are some caveats of note. Firstly, the most expensive portion of the overall calculation is the solution of the Poisson equation for pressure. The speedup in the calculation of the coefficients, while impressive, will therefore only affect the overall calculation fractionally, and not proportionally. Notwithstanding, the overall gains can be significant, particularly in the presence of many particles. Secondly, we note that the SP is only faster for moving particles. Since the cage of stationary particles is fixed, the SVD needs only to be computed once, and the resulting decomposition can be reused, resulting in comparable speed with the SP. If the particle is moving, on the other hand, the cage moves with it (once the particle moves a distance of from the center of the cage, a new cage is constructed Zhang:2005ty ()), requiring a call to the SVD. Although this does not happen at every time step, the effect can be nonetheless significant as shown in section 5.5.
5.2 Recovery of assigned coefficients in Stokes flow
The simplest test of the accuracy with which the coefficients are calculated can be carried out on a flow field generated analytically by using the expressions (11) and (12). We use these relations with arbitrarily chosen coefficients to assign the fields at the grid points. We choose (for a total of coefficients assigned), and a spatial resolution of . The SPs are computed at . Given that the entire field in this case is Stokes, the solution is, in theory, insensitive to the choice of . However, the further out one is from the surface of the particle, the finer the Cartesian grid becomes relative to the scale () on which the vary. Later we will demonstrate how the accuracy varies with for flows in which the Stokes region has a finite extent.
To quantify the accuracy with which the assigned coefficients are recovered from the assigned fields, we define an error metric as
(37) 
where and respectively denote the assigned and recovered values; and are defined in the same way.
Figure 2 shows the results so obtained. Generally, the recovery by all methods is reasonable, with a maximum error level of , achieved by the SVD. In fact, all of the SPs prove more accurate than the SVD under the error metric. In contrast the SVD incurs a larger error at higher values of . This instability, presumably associated with the numerical implementation of the SVD, is particularly apparent in , but can also be noted in and . In terms of and , the SP involving the interpolated proves most accurate, followed by those involving and , and the linearization estimate of (recall, from section 4, that linearization is only used with and ). The interpolated also performs best in terms of , but now followed by , , and the linearized .
Although the reasons for the differing performance, in terms , , and , are somewhat unclear, we can make some observations. Firstly, we note that the interpolated is likely more accurate than due to its vector nature: the relevant SPs (equations (24) and (28)) involve two velocity components, presumably resulting in some cancellation of error. While the same is true of , it should be remembered that this is a derived quantity and some error will be associated with the differentiation of the velocity field. The results involving the linearized approximations in equations (32) and (33) are likely inferior to the interpolated ones due to firstorder accuracy of the former. Stokes flow with contains relatively rapid variations, which are better captured by the secondorder method. This is evidenced by setting , at which degree the two methods produce nearly the same results (not shown). A further confirmation is provided by the results of the next sections where the two methods show comparable performance in smoother flows in which only the lowest order coefficients () prove significant.
5.3 Drag on a sphere in a triplyperiodic array
In this section we examine the accuracy of the SVD and SPs in a pressuredriven flow through a triplyperiodic array of spheres. As in Zhang:2005ty (), we write the pressure gradient as
(38) 
in which is periodic, is the unit vector in the flow direction , and is a constant. By applying a momentum balance in the flow direction at steady state, it can be shown that the force on the particle is given by
(39) 
where . We use a periodic domain of size , with the particle fixed at the center. The flow is initially stationary, at which time a pressure gradient is applied. This accelerates the flow until it reaches a steady state where the drag force on the particle matches that applied by the pressure gradient, resulting in a Reynolds number of , where is the mean flow rate. The accuracy of the force calculation in equation (14) can then be gauged by comparing with the exact result in equation (39).
Using a resolution of , we calculate the error as a function of , with the SPs being computed at . Figure 3 show the results so obtained. Good results are obtained for all methods, with the error being less than for , as was the case in Zhang:2005ty (). Notable is the performance of SP involving , suffering the lowest error out of all the methods, including the SVD. The two SPs involving now prove comparably accurate, with the interpolated one only having a slight advantage. This is likely due to the smoothness of the flow field, as the coefficients are dominant. Incidentally, these are the coefficients involved with the force calculation, as can be seen from equation (14). In contrast, the SP involving suffers the largest error of all. While all velocities vanish at the surface of the particle, the radial component does so more rapidly than the transversal ones. This causes larger errors in the interpolation, as the signal to noise ratio is diminished. This explains why proves less accurate here than in the previous section, where the SPs were computed at .
Since the SPs are computed at a given distance from the particle surface, this test case provides an opportunity to explore the extent of the Stokes region, at least for this flow configuration. As one moves closer to the surface, SPs involving velocity will at some point be increasingly in error as the velocity vanishes. The SPs involving vorticity and pressure will also suffer error in this region, as neither is continuous through the surface of the cage where the boundary conditions are applied. Thus, upon moving towards the cage, one will eventually interpolate through this discontinuity, resulting in error. The error will also grow as we move further from the particle and out of the Stokes region.
Figure 4 shows error calculations analogous to those in figure 3, but now as a function of , with fixed at 3. The behavior of the different SPs with is not uniform, with having a minimum in error at , while the minimum for those involving velocities occurs close to .
Lastly, we note that, as this test case is nearly axisymmetric in the vicinity of the particle, the coefficients are very small and as such we cannot draw conclusions as to which SP serves best to compute this. This is addressed in the next subsection.
5.4 Couple on a sphere in a triplyperiodic array
In the previous section we presented results demonstrating the accuracy with which the coefficientcalculation of the SVD and different SPs conserves the momentum of flow driven by a pressure gradient. In this section we perform an analogous set of calculations, but this time involving the conservation of angular momentum. This is a stringent test of which few examples, if any, can be found in the literature.
A particle rotated at a constant rate will impart angular momentum to the surrounding fluid. In turn, in a finite or periodic domain, at steady state the couple at the surface of the particle must be balanced by the couple on the fluid on the outer boundary. The calculation proceeds as follows.
We take the crossproduct of , the position vector relative to the particle center, with the momentum equation (1) at steady state and in the absence of gravity, and integrate over the domain
(40) 
where
(41) 
By the symmetry of the stress tensor we have
(42)  
An application of the divergence theorem transforms equation (40) to
(43) 
where is the surface of the domain and that of the particle; is the normal vector directed out of the fluid. We have dropped the integral over the particle surface in the lefthand side as the radial component of velocity vanishes there. For periodic boundary conditions, which we assume, the remaining integral in the lefthandside also vanishes. Furthermore
(44) 
is just the hydrodynamic couple on the particle (with the minus sign due to the direction of the normal) and equation (43) therefore becomes
(45) 
We take so that and . We use a domain of size , a grid resolution of , and a rotation rate of . We compute the error , as a function of , analogously with the previous section. Figure 5 shows the results. With the exception of , all SPs prove more accurate than the SVD for . At first glance, this seems at odds with the results of the previous section, where was the most accurate SP. However, with , and, therefore, . Hence, the second component of this SP is lost. Further, a significant portion of is associated with the solidbody rotation of the particle, which gets subtracted out in the frame of the particle, decreasing the signaltonoise ratio in the calculation of , which in turn causes higher error levels. In contrast, the SP associated with now performs best out of all candidates, which is also at odds with the results of the previous section. As in the previous section, the SPs involving the interpolated, and linearized, estimates of provide practically the same results.
We also compute the above error for fixed , but varying the at which the SPs are computed. Figure 6 shows the results so obtained. As was the case in figure 4, the error for different SPs does not behave uniformly in . Generally, though, it seems that there is some leeway in the placement of , with providing reasonable results.
The difference between the performance of the various SPs in this section and the previous section is instructive. Clearly, which SP will perform best depends on the flow conditions. If, for example, is very small, it will not deliver a reliable estimate of the (via equation (30)). On the whole, it seems that the combination of and provides consistently good estimates across widely varying flow conditions. One would further suggest using the linearization over the interpolation, given the simplicity of the former and its other advantages described in section 4.
5.5 A falling particle
The test cases in the previous sections were all for nontranslating particles. We now describe some tests on a freely falling particle in a parallelepipedal container of length and cross section . The boundary conditions in the lateral and vertical directions are periodic and noslip, respectively. The ratio of the particle density to the fluid density is , the kinematic viscosity is m/s, and the acceleration of gravity has the standard value =9.81 m/s. The simulation is carried out with and 8, respectively using and 3. The timestep in all cases was , resulting in a CFL number of 0.45 at . In light of the previous results, we narrow our scope of SPs to those involving the linearization of and , and those involving the quadratically interpolated and .
The particle is released from rest at . If no provision is made, at the weight of the particle would be suddenly imposed on the fluid, thus generating strong force oscillations. To avoid this artifact we have ramped up the value of the acceleration of gravity in an exponential manner over a time .
Figure 7 shows the evolution of the Reynolds number . In an unbounded fluid the steady state Reynolds number predicted by the SchillerNaumann correlation Schiller:1933vh () is 25, but, due to the lateral constraint, here it is close to 22.5 in all cases. Essentially, the effect of the constraint can be viewed as a blockage effect due to the image particles. While the steadystate Reynolds number varies slightly between the different methods and at the two different resolutions, the largest difference is less than 3%.
Since the cage changes with time, for the reasons explained in section 5.1, this is also a good test case to examine the CPUtime used by particlerelated operations. Figure 8 shows the cumulative time in seconds thus expended, as measured using an Intel 2.66GHz CPU. As expected, both SPs provide a significant speedup, particularly for . For , the SPs have cumulatively cost roughly 14 s at , while the SVD costs roughly 109 s, larger by a factor of 8. The cost associated with the SVD increases by a factor of 6 upon increasing the resolution to , while that of the SPs only increases by a factor of 2 (recall that and 3 respectively for and 8). This is consistent with the time measurements demonstrated in figure 1. The cumulative savings in the cost of particlerelated operations would be proportionally increased in the presence of multiple particles. Lastly, the cost associated with the SP based on the linearized velocity and pressure is roughly 50% less than that based on the interpolated transversal vorticity and pressure.
It is generally the case with immersedboundary methods that force oscillations are produced as the particle translates with respect to the underlying grid. It was shown in Zhang:2005ty () that physalis is not immune from this problem, and it is therefore of interest to examine the performance of the SPs also in this respect.
Generally speaking, the origin of the force oscillations must be sought in the minor discontinuities that the calculated force undergoes as the particle moves. The same mechanism is responsible for the force oscillations encountered upon a sudden release of the particle as mentioned before. A discontinuity in the hydrodynamic force on the particle, by the actionreaction principle, causes an impulsive force on the fluid, which responds with a pressure impulse. This is the reason why the artificial compressibility approach of Perrin:2006cl (); Perrin:2008dh () suffers from this difficulty much less than the original physalis.
In our method, at the end of each timestep, the position of the particle center is compared with that of the cage center and, if the difference is greater than , the cage is displaced Zhang:2005ty (). As a consequence, a different set of nodes is used to compute the coefficients at the next time step, and this causes an unavoidable small discontinuity in the force calculation even though the displacement of the cage is small enough that both the old and the new cages are entirely in the Stokes region. This mechanism leads one to expect that the force oscillations decrease with a refinement of the grid, as it was indeed observed in Zhang:2005ty () and as will be shown presently. Another contributing factor may be the fact that the location at which the boundary conditions are applied on the fluid is different when the cage is moved. Furthermore, some of the grid points that were previously inside the cage become fluid points and need to be reinitialized as such, as described by Zhang:2005ty (). However, this is likely not a significant source of force oscillations due to the iterations in the inner loop (steps 14 in section 2).
Figure 9 shows the time evolution of the normalized force on the particle for and 8. It should be noted that, for clarity, the curves for the original SVD method and the SP with have been translated by . The expected decrease of the amplitude of oscillations with increasing grid refinement is evident for all three methods, although the SVD curves show significant spikes even with the finer grid. The amplitude of the oscillations with both SP methods is smaller, does not present any spikes and exhibits a greater reduction on the finer grid. The greatest reduction is obtained by the SP with the linearized which thus is seen to perform well also from this point of view.
5.6 Taking the scalar products closer to the particlesurface
In the above sections, the SPs were taken over a sphere with a radius with with , where is the radius of the particle. While this does not arise concerns in singleparticle simulations, it might under circumstances in which particles are expected to be in close proximity (such as for flow through a tightlypacked porous medium) and even to collide. In this section we show that is not an intrinsic limitation of our method and that even is permissible, allowing interparticle contact.
As noted before, vorticity and pressure are not continuous through the cage. This is the reason why, for a given cage size and grid resolution, there exists a lower limit on below which the calculation of vorticity and pressure becomes in error. While this is the case for both the interpolation and local linearization, the latter is less affected due the smaller number of cells involved, thus allowing a smaller before the cage is breached.
The constraint of operating outside the cage can be met either by using appreciably greater than 1, as done before, or by decreasing as we now show. The latter option is permissible since the analytic expressions (11) and (12) embody a smooth analytic continuation of the pressure and velocity fields inside the particle. Thus, the boundary conditions on the cage may be applied on the inside as well.
As shown in figure 10, by taking , the SP involving gives excellent results even with . As the velocities vanish at , the relevant SPs are increasingly in error close to the surface, even with . In situations where is required, one would thus be advised to use . Actually, in all the examples considered, except that of section 5.4, this option gives excellent results and the associated computational overhead is only modestly larger than for the linearization method.
6 Summary and conclusions
In this work we have presented a new approach to calculating the unknown expansion coefficients in Lamb’s solution for the nearparticle flow field in the physalis method. In the previous approach of Zhang:2005ty (), the coefficients were obtained via the SVD solution of a set of linear equations formed by equating the analytical expansion, evaluated at the nodes of the finitedifference field around the particle, to the finitedifference field itself (cf. equation (16)). In the new approach we take advantage of the orthogonality of the set of vector harmonics by taking suitable scalar products thereof with values obtained from the finitedifference field.
We have described how the finitedifference field can be estimated on the spherical surface over which the scalar products are taken and further how the relevant integrals are computed efficiently via a combination of discrete Fourier and Legendre transforms on the sphere.
There are several options for which fields should be used to calculate the coefficients. Our results suggest that the local linearization method based on the pressure and the transverse velocity (section 4.2) is efficient and accurate provided the scalar product can be computed on a surface with a radius exceeding that of the particle by 12 mesh lengths. For simulations where the particles are very close or can collide, a more accurate option would be the use of the transverse vorticity in addition to the pressure (sections 3.1 and 3.2).
Via various examples involving stationary, rotating, and sedimenting particles, we have demonstrated that the new approach is more accurate, more stable in terms of forceoscillations, and significantly faster than the SVD approach of Zhang:2005ty (). The speedup is particularly significant at higher Reynolds numbers, as . As the number of points comprising the cage scales as , the cost of the SVD scales as , (cf. section 5.1). While the same resolution is of course required when using the scalar product, the associated cost scales only as . The speedup reported in this paper are therefore likely to be further enhanced at higher Reynolds numbers.
Acknowledgements
This work was supported by the IMPACT Institute, the Netherlands.
Appendix
The , in equation (13) and others, are the orthogonal spherical harmonics,
(46) 
where and are respectively the polar and azimuthal angles, and . For , the associated Legendre polynomials are given by
(47) 
The constant is included in order to render the orthonormal.
The regular and singular solid harmonics can be related through the noslip boundary condition on the particle surface. The resulting expressions are
(48)  
(49)  
(50) 
References
 (1) M. Uhlmann, An immersed boundary method with direct forcing for the simulation of particulate flows, J. Comput. Phys. 209 (2) (2005) 448–476.
 (2) M. Uhlmann, Interfaceresolved direct numerical simulation of vertical particulate channel flow in the turbulent regime, Phys. Fluids 20 (5) (2008) 053305.
 (3) F. Lucci, A. Ferrante, S. Elghobashi, Modulation of isotropic turbulence by particles of Taylor lengthscale size, J. Fluid Mech. 650 (2010) 5–55.
 (4) A. Prosperetti, H. N. Oguz, Physalis: A New o(N) Method for the Numerical Simulation of Disperse Systems: Potential Flow of Spheres, J. Comput. Phys. 167 (1) (2001) 196–216.
 (5) S. Takagi, H. Oguz, Z. Zhang, A. Prosperetti, PHYSALIS: a new method for particle simulation  Part II: Twodimensional NavierStokes flow around cylinders, J. Comput. Phys. 187 (2) (2003) 371–390.
 (6) Z. Zhang, A. Prosperetti, A secondorder method for threedimensional particle simulation, J. Comput. Phys. 210 (2005) 292–324.
 (7) Z. Zhang, L. Botto, A. Prosperetti, Microstructural Effects in a FullyResolved Simulation of 1,024 Sedimenting Spheres, IUTAM Symposium on Computational Approaches to Multiphase Flow (2006) 1–10.
 (8) A. Naso, A. Prosperetti, The interaction between a solid particle and a turbulent flow, New Journal of Physics 12.
 (9) Q. Liu, A. Prosperetti, Wall effects on a rotating sphere, J. Fluid Mech. 657 (2010) 1–21.
 (10) Q. Zhang, A. Prosperetti, Pressuredriven flow in a twodimensional channel with porous walls, J. Fluid Mech. 631 (2009) 1–21.
 (11) Q. Liu, A. Prosperetti, Pressuredriven flow in a channel with porous walls, J. Fluid Mech. 679 (2011) 77–100.
 (12) A. Perrin, H. H. Hu, An explicit finitedifference scheme for simulation of moving particles, J. Comput. Phys. 212 (1) (2006) 166–187.
 (13) A. Perrin, H. H. Hu, An explicit finite difference scheme with spectral boundary conditions for particulate flows, J. Comput. Phys. 227 (20) (2008) 8776–8791.
 (14) H. Lamb, Hydrodynamics, Cambridge University Press, 1932.
 (15) A. Prosperetti, Advanced Mathematics for Applications, Cambridge University Press, 2011.
 (16) S. Kim, S. J. Karrila, Microhydrodynamics, principles and selected applications, Dover Publications, 2005.
 (17) D. Healy, D. N. Rockmore, P. Kostelec, S. Moore, FFTs for the 2sphereimprovements and variations, J. Fourier anal. appl. 9 (4) (2003) 341–385.
 (18) B. Gough (Ed.), GNU Scientific Library Reference Manual, 3rd Edition, Network Theory Ltd., 2009.
 (19) C. L. Lawson, R. J. Hanson, Solving least squares problems, Prentice Hall, 1974.
 (20) C. E. Shannon, Communication in the Presence of Noise, Proc. IRE 37 (1) (1949) 10–21.
 (21) L. Schiller, A. Naumann, Fundamental calculations in gravitational processing, Z. Ver Dtsch. Ing. 77 (1933) 318–320.