A Model for Translation and Rotation Resistance Tensors for Superellipsoidal Particles in Stokes Flow

: In this paper, forces and torques on solid, non-spherical, orthotropic particles in Stokes ﬂow are investigated by using a numerical approach on the basis of the Boundary Element Method. Different ﬂow patterns around a particle are considered, taking into account the contributions of uniform, rotational and shear ﬂows, to the force and the torque exerted on the particle. The expressions for the force and the toque are proposed, by introducing translation, rotation and deformation resistance tensors, which capture each of the ﬂow patterns individually. A parametric study is conducted, considering a wide range of non-spherical particles, determined by the parametric superellipsoid surface equation. Using the results of the parametric study, an approximation scheme is derived on the basis of a multivariate polynomial expression. A coefﬁcient matrix for the polynomial model is introduced, which is used as a tunable parameter for a minimization problem, whereby the polynomials are ﬁtted to the data. The developed model is then put to the test by considering a few examples of particles with different shapes, while also being compared to other, currently available solutions. On top of that, the full functionality of the model is demonstrated by considering an example of a pollen grain, as a realistic non-spherical particle. First, a superellipsoid, which best ﬁts the actual particle shape, is found from the considered range. After that, the coefﬁcients of the translation, rotation and deformation resistance tensors are obtained from the present model and compared to the results of other available models. In the conclusion, a superior accuracy of the present model, for the considered range of particles, is established. To the best of the authors knowledge, this is also one of the ﬁrst models to extend the torque prediction capabilities beyond sphere and prolate particles. At the same time, the model was demonstrated to be simple to implement and very conservative with the computational resources. As such, it is suitable for large scale studies of dispersed two-phase ﬂows, with a large number of particles.


Introduction
Multiphase flows with non-spherical particles are being extensively researched in recent years, by authors from numerous fields of engineering and science.Nowadays, with significant development of computational methods, one of the greatest challenges still remains in tackling industrial level simulations.This follows from a fact that they are frequently associated with complex carrier fluid flow fields and a large number of dilutely distributed, non-spherical particles.Problems of this nature are usually treated by using the Euler-Lagrangian approach, based on the point-particle assumption, which is valid for particle sizes much lower than the Kolmogorov length scale of the flow.Such an approach relies on using a model for the prediction of the particle-fluid interaction, instead of using the computationally expensive direct approach.Albeit the simplifications, this approach was proven successful for a number of problems considering disperse flows [1], where small particles are involved.A typical starting point for the formulation of particle-fluid interaction is the Maxey-Riley equation [2] that was originally developed for spheres.It has been repeatedly shown that for sufficiently high particle-to-fluid density ratios (i.e., in solid-gas suspensions), the most important terms are gravity and drag [3][4][5], which dominate the particle translation.The simplest model for drag is that of a sphere in a creeping flow, which was later extended to other flow regimes [6][7][8].Due to its simplicity, accuracy and the fact that the drag formulation for a sphere has been available for the longest, researchers have extensively relied on the assumption of the spherical particle shape.Significant research interest still exists for this approach, as there is a number of applications [9][10][11][12], where the best representation of the particle shape is in fact spherical.On the other hand, it is not uncommon that a reasonable shape specific model does not exist for the actual particle shape and flow regime [13], in which case a spherical drag model can still be used, at least to gain an approximate overview of the observed phenomena.Unlike for the spheres, drag models for non-spherical particles are significantly more complex, not only because of the increased shape formulation effort, but also due to the fact that the particle orientation, with respect to the flow velocity field, needs to be accounted for.One of the first attempts to consider the drag of a simplest non-spherical particle had been carried out by Oberbeck [14], who derived the analytical solution for the drag of an arbitrarily oriented ellipsoid in an uniform creeping flow.Despite the limitation of the linear flow field around the particle, this is still one of the very few models available for non-spherical particles, which is based on the exact expression for the drag prediction.Among the many studies considering the ellipsoidal particles, prolate [15,16] and oblate [17] ellipsoids still seem to be the most popular choices.Nevertheless, for the comprehensive description of the particulate phase involving non-spherical particles, employing only the drag model is not sufficient, as it has been repeatedly demonstrated [18][19][20].It is now well understood that the rotation can significantly affect the particle dynamics and other associated phenomena, such as their deposition.With this in mind, one can no longer neglect the importance of employing an accurate model for the torque, along with those for the drag.Much like for the drag, ellipsoids [21] are one of the very few geometries for which the exact torque expression is available, based on the analytical solutions of the creeping flow field around the particle.Together with the ellipsoid drag model, both have been proven very useful when applied to the problems with fibrous micro-particles, for example in bio-engineering [22,23], medicine [24], paper-making industry [25] and in nuclear reactors processes [26].That said, it is important to remember that non-spherical particles, encountered in a number of natural or industrial flows, are often distinctively irregular, up to the extent where readily available exact models fail to deliver a reasonable solution.To overcome that, drag models based on generalised shape descriptions were proposed as a result of the extensive experimental work throughout the years.In those models, the particle shape is characterised by shape descriptors, i.e., sphericity [27], which is then used as a relation of the particle shape with the drag prediction.Some of the more advanced models also employ additional shape descriptors [28] and even consider the particle orientation [29].Examples that can be found in the literature range from the studies of volcanic ash transport in the atmosphere [30,31], as well as investigations of ice crystal trajectories in an aircraft turbofan compressor [32], to simulating the formation of polymer particles in a downer reactor [33].
Taking into account the facts from the presented literature, we identify the lack of available drag and torque models, based on the exact formulation, for a wider range of non-spherical particle shapes.With the exception of ellipsoids, almost all other nonspherical particles can still only be captured by generalised shape description models, which inevitably comes at a cost of the prediction accuracy [34].For that reason, we propose an extension to the existing models for ellipsoids by introducing a multi-parameter surface description based on the superellipsoidal equation.Some authors estimate [35] that as much as 80% of solids can be represented by superquadratic functions.In comparison with ellipsoids, this expands the parametrization capabilities to a significantly wider range of shapes, such as cylinders, cuboids, rhomboids and even octahedrons.In order to retain the benefits of the exact models, the drag and torque expressions are employed by considering the actual particle shape and its orientation.This allows for the treatment of all non-spherical particles that can be represented by superellipsoids, with essentially the same accuracy as analytical models provide for ellipsoids.Moreover, with this approach, new insights into the complex-shaped particle dynamics can be gained, by exposing the details that most of the generalised shape models would not be able to reconstruct.Hereinafter we present the superellipsoid drag and torque modelling framework, by means of which we obtain the necessary data to derive the required correlations.Following from that, we demonstrate the validity of the numerical approach by conducting a validation study, including the evaluation of mesh resolution and computational domain sizing.In the end, we present the final drag and torque models for superellipsoidal particles, for use in Lagrangian particle tracking applications, together with its accuracy evaluation and the functional demonstration.

Governing Equations
Let us consider the relative flow of a carrier fluid, around a small particle.As discussed in the introduction, we are limiting ourselves to a particle size smaller than the Kolmogorov length scale of the flow.Typically, this also results in a correspondingly small particle Stokes number 2 , where ρ p and ρ f are the particle and fluid densities respectively, d p is the particle diameter and η is the Kolmogorov length scale.It has been repeatedly shown that for a sufficiently small Stokes numbers, the particle almost perfectly traces the fluid streamlines, without any significant detachments from the flow [36].This in turn means that the relative flow velocity around a particle is low and typically falls well into the viscous (Stokes) regime.As a result, we assume a steady state incompressible flow of a Newtonian fluid around a particle, at a very small particle Reynolds number, Re p = || u|| 2 d p /ν 1, where || u|| 2 is the relative flow velocity magnitude and ν is its kinematic viscosity.In this case the advection term in the Navier-Stokes equations can be neglected, which results in the creeping flow (Stokes) equations: Here g is the gravitational acceleration and σ is the Cauchy stress tensor defined as σ = −PI + τ, where P is the pressure, I the identity tensor, and τ the viscous stress tensor.Considering a Newtonian fluid, we model the viscous stresses as , where µ = νρ f , rendering the following form of the Stokes equation Finally, we recognise that gravity is a conservative force, which may be written as a gradient of gravitational potential Φ and introduce a modified pressure as p = P − ρ f Φ, where g = ∇Φ.The final form of the Stokes equation reads The Stokes flow Green's functions satisfy the continuity equation ∇ • u = 0 and are the solutions of the singularly forced Stokes equation.Defining r = r − r 0 and r = |r|, the 3D free-space Green's functions can be written as: Let q = σ • n denote the boundary traction, i.e., the flux of momentum into or out of the boundary.The boundary integral representation for the Stokes problem [37] is: where c( r 0 ) = 2α is twice the solid angle as seen from the point r 0 , i.e., in the interior of the domain c = 8π, at a smooth boundary c = 4π.The normal vector n points into the domain.
The first term on the r.h.s represents the double layer potential of three-dimensional Stokes flow, and the second term is the single-layer potential of three-dimensional Stokes flow.

Boundary Element Solution of Stokes Flow Over a Particle
The boundary integral Equation ( 5) is used as the basis of the Boundary Element Method (BEM) solver, which we developed.The numerical implementation is based on our Laplace BEM solver [38,39].We consider the boundary Γ = ∑ l Γ l to be decomposed into boundary elements Γ l : where n k is the k component of the normal vector pointing from boundary element l into the domain.
Let Φ and Ψ be the interpolation functions used to interpolate the function and flux values within the boundary elements, i.e., , where is the m th nodal value of the function within the l th boundary element.This yields The integrals above depend only on the mesh geometry and are independent of the flow solution.As such, they may be calculated in advance and stored.Since the boundary elements share nodes, the number of integrals that need to be stored is smaller than the number of integrals calculated, since the integrals, which are needed by the same node, can be summed up.To obtain a system of equations from unknown velocities and tractions, we place the source point in all boundary nodes.Storing the integral values in matrices ([T ij ], row corresponding to different source points, columns to different nodes in the mesh) and flow quantities in nodal vectors {u i }, we obtain the following system of equations: Boundary conditions include known values for {u x }, {u y }, {u z }, {q x }, {q y }, {q z }.Collocation points are placed only into nodes, where the value is unknown.A system of linear equations is set up for all unknowns, where in case of unknown {u x } or {q x } Equation ( 8) is used, in case of unknown {u y } or {q y } Equation ( 9) is used and in case of unknown {u z } or {q z } Equation ( 10) is used.Additional details of the BEM employed, such as implementation of integration, can be found in [38,39].
When considering flow over a particle using this method, a no-slip velocity boundary condition is prescribed at the particle surface and the algorithm renders boundary tractions.These can be integrated over the particle's surface to obtain the force and torque exerted on the particle by the fluid.The main advantage of using BEM for this task over domain based methods such as finite volumes is that boundary tractions are obtained by direct solution of the system of equations and not post-processed from the velocity fields.Additionally, only the particle surface must be discretizised, a volume mesh is not necessary, leading to faster computational times.

A Particle in Stokes Flow
We aim to investigate the fluid flow induced forces and torques on non-spherical particles.Our main goal is to extend the currently available modelling options for non-spherical particles, as highlighted in the introduction.Since we want to consider a wider range of particle shapes, it is a good idea to employ a parametric surface definition approach.In order to enable shape manipulation options, we introduce a superellipsoidal particle [40][41][42] with four parameters λ 1 , λ 2 , e 1 and e 2 and fix its volume to V = π/6, which is the volume of a sphere with diameter d p = 1.The three semiaxes are where B(•) is the beta function.Placing the frame of reference in the centre of the superellipsoid and aligning the coordinate system axes with the semiaxes of the superellipsoid, we can define the surface of the superellipsoid as In Figure 1 a couple of representative superellipsoidal particles are shown highlighting the wide variety of shapes, which can be modelled by the superellipsoidal parametrization.
The particle frame of reference (PFR) is also shown.We consider a particle moving in a fluid, where u is the fluid velocity in the PFR.The force on the particle can be calculated by where q = σ • n is the boundary traction and K is the translation resistance tensor expressed in the PFR.The boundary traction represents the sink of momentum at the particle surface and thus its integral over the surface yields the force.When we consider a spinning particle in a stationary fluid, the torque is T = Γ r × qdΓ = −µπc 3 Ω • ω where Ω is the rotation resistance tensor and ω is the angular velocity of the particle in the PFR.The flow shear also affects the rotation dynamics of the particle.Taking both contributions into account, the torque can be expressed as where Π is the deformation resistance tensor.Here f , g, h are the elements of the deformation rate tensor and ξ, η and χ are elements of the spin tensor, defined as Making use of the Lorentz reciprocal theorem [43], Happel and Brenner [8] showed that the translation resistance tensor and the rotation resistance tensors are symmetric.As an immediate consequence of the properties of symmetric tensors, every particle must posses principal axes of translation, i.e., three perpendicular directions such that translating without rotation to one of them, the particle experiences a force only in that direction.In the case, when the particle is orthotropic (possess three mutually perpendicular symmetry planes), it follows from symmetry considerations that the principal axes of translation are normal to the symmetry planes.Moreover, similar can be applied for the particle rotation in a uniformly rotating and shearing flow.If we employ a two-dimensional uniform flow, with pure rotation or pure shear, blank with the particle symmetry planes, the resulting torque is limited to a single principal axis, normal to the symmetry plane where the flow is employed.Thus, considering orthotropic particles and choosing the particle reference frame blank with symmetry planes, the translation K, rotatation Ω and deformation Π resistance tensors are diagonal.That said, one must determine nine components to completely describe the force (13) and torque (14) experienced by a particle in an arbitrary creeping flow field.The K, Ω and Π tensors are expressed with respect to the particle frame of reference.Using a general tensor rotation formalism, we can transform them into the inertial (Eulerian) frame by writing K = R T KR, where K is the translation resistance tensor in the inertial frame of reference and R is the rotation matrix.The transformation of Ω and Π tensors is analogous to that of the K tensor.

Numerical Setup
For the purpose of this study, we employ a three-dimensional numerical model of a superellipsoidal particle in creeping flow.The numerical domain, in the shape of a sphere (Figure 2), is generated with the particle under consideration positioned in its centre.
Since Stokes flow is diffusive in nature and the simulations are steady state, we expect the perturbation of the fluid due to the particle to be non-negligible in a very large region.The Stokes's analytical solution for flow over a sphere reveals that the perturbation in flow velocity diminishes as 1/r measuring away from the particle.This means that the boundary conditions must be placed far away from the particle in order not to have a significant effect on the force and torque calculations.Therefore, in order to correctly assess the flow field and the resistance and rotation tensors, we take into account the findings of some recent studies listed below.Chadil et al. [44] performed a DNS study of the flow around the sphere, where they used an rectangular domain of 16d p × 8d p .Likewise, Zastawny et al. [45] also employed an rectangular domain, sized to 20d p × 20d p × 10d p , in their immersed boundary method (IBM) setup for non-spherical particles.An extensive study of low-Re flow around non-spherical particles was conducted by Andersson and Jiang [46], who selected a volume equivalent cube edge length as a comparative domain measure.Edge lengths from 20d p to 170d p were examined, where the length of 135d p was claimed to produce a sufficiently domain independent solution.In a similar fashion, Štrakl et al. [47] employed a cube-shaped domain, with the base edge length of 160d p , for their parametric study of drag and lift forces acting on non-spherical particles.At least in the latter of the four studies, the authors had to accept a compromise regarding the domain size selected for the study, as the mesh element count was required to be kept within reasonable limits.This is typically an inevitable consequence of the available computational resources and time, especially in the case of parametric studies, where a high number of simulations is required.Overcoming this limitation comes as one of the main advantages of using BEM, where meshing the volumetric domain is not required.With this in mind, if we proceed to the question of the domain size in the present work, we are able to chose the diameter of the outer sphere almost arbitrarily, which we finally set to a value of 1024d p .In the case of creeping flow around a sphere, using the domain of this size would reduce the effects of the domain boundary to ∼0.1%, compared to ∼1%, as achieved by the domains used in [46,47].The boundary condition, imposed at the surface of the outer sphere (domain boundary), is a known velocity field.We decide to process the contributions of different resistance tensors separately.Following from that, we employ nine different velocity fields, to take in to account all three principal axes for the translation, rotation and shear flow.In the case of estimating the translation resistance tensor K we prescribe the velocity corresponding to the translation of the particle in a stationary fluid, along the three axes of symmetry, namely u = − î (Figure 3, left).The unit vectors ( î, ĵ, k) define the standard basis of the coordinate vector space.To obtain the rotation resistance tensor Ω, we prescribe the velocity field corresponding to a uniform rotation of the particle, around the three axes of symmetry, namely u = −z ĵ + y k (Figure 3, centre).Analogously, for the deformation resistance tensor Π we prescribe the two-dimensional pure shear velocity field, normal to the three symmetry axes, namely u = z ĵ + y k (Figure 3, right).At the same time, we prescribe a no-slip condition at the surface of the particle, to close the set of the required boundary conditions.According to ( 8)- (10), one can see that the result of such simulation is a traction at the particle surface q, which is displayed in Figure 3.When integrated according to Equations ( 13) and ( 14), this yields the force and the torque exerted on the particle by the fluid.

Mesh Study and Validation
To validate the proposed algorithm and to estimate its accuracy, we first consider a particle in the shape of a prolate ellipsoid, by varying the axial ratio λ 1 and keeping the rest of the superellipsoidal parameters at unity (λ 2 = 1, e 1 = 1, e 2 = 1).For such a particle, Jeffery [21] analytically derived the non-zero components of the three resistance tensors in creeping flow.Ravnik et al. [19] established that their validity is limited to particles, whose size is much smaller than the Kolmogorov length scale.The expressions for the tensor components, as given by Jeffery [21], are and The nondimensional coefficients α 0 and γ 0 were defined by Gallily and Cohen [48] as Starting from the special case of λ 1 = 1 (where the particle is in fact a sphere) the two nondimensional coefficients are lim λ 1 →1 α 0 = lim λ 1 →1 λ 2 1 γ 0 = 2 3 , while the resulting nonzero tensor components are K ii = 6 and Ω ii = 8.We continue with two additional prolate particles, by setting λ 1 = 2 and λ 1 = 5, for which we compute the tensor components in accordance with ( 17)- (19).Together with the spherical particle, they are used as a reference, to confirm the validity of the proposed numerical approach.The domain discretization is employed in two steps.First a domain boundary (the outer sphere) is meshed separately with approximately 4.1 • 10 3 degrees of freedom, which is kept constant for all cases.After that, the particle meshing is employed, where in total four different meshes are considered.We begin with mesh A, which has 7.4 • 10 3 degrees of freedom and continue with mesh B with 15.9 • 10 3 , mesh C with 31.1 • 10 3 and mesh D with 48.4 • 10 3 degrees of freedom.As we are trying to capture different particle shapes, some discrepancy between the meshes was unavoidable, so the above given DOFs are the average values for all particles.The results, in terms of the diagonal components of the translation, rotation and deformation resistance tensors, obtained for all meshes, are given in Table 1.As evident, the numerical results compare well with the analytical values.The average relative errors of all three particles, for all non-zero tensor components, are 2.1% for mesh A, 0.21% for mesh B, 0.15% for mesh C and 0.14% for mesh D.Moreover, to also demonstrate the adequacy of the mesh, when capturing real superellipsoids, we select two additional particles by employing all four superellipsoidal parameters.The first, a rectangular disc particle, is obtained by using λ 1 = 5, λ 2 = 5, e 1 = 0.2, e 2 = 0.2 and the second, a pyramid-like particle, by using λ 1 = 10, λ 2 = 10, e 1 = 1.5, e 2 = 1.5.We repeat the above process, where all non-zero tensor components are computed, for the same four mesh densities.In Figure 4, the meshes A-D of the rectangular disc particle are visualised.Due to the lack of analytical expressions and the lack of availability of experimental results, we have to base our verdict regarding the mesh accuracy on the convergence of the values of interest alone.This is displayed in Figure 5, where the resulting tensor components are plotted against the number of DOF for each mesh.The resistance tensor components K xx and K zz are shown in Figure 5a, Ω xx and Ω zz in Figure 5b and Π xx in Figure 5c.We notice that not all tensor components have been presented, which is supported by the fact that the two considered particles possess both axial ratios of the same value (λ 1 = λ 2 ), so the following simplifications are observed in the resulting tensor components: K xx = K yy , Ω xx = Ω yy , Π xx = −Π yy and Π zz = 0.In general, the convergence behaviour is similar to the previous study, where spherical and prolate particles have been considered.Albeit some indications already appeared in the results from Table 1, the plots in Figure 5 clearly show the lesser sensitivity of the K tensor components to the mesh density, compared to the components of the Ω and Π tensors.For example, in case of the rectangular-disc particle, the average change in the K tensor components, from mesh A to mesh D is around 1.8%, while the same for Ω and Π yields around 5.1%.Nevertheless, if we take mesh B for our starting point, the average tensor component changes reduce to approximately 0.4% for K, 1.0% for Ω and 1.1% for Π.We can observe that all of the tensor components start to settle at the latest around 20 • 10 3 DOF, which is consistent with what we have already established when considering the sphere and the prolate ellipsoids.After the extensive mesh evaluation, we are choosing mesh B as appropriate for the main parametric study.This mesh is also considered to have the optimal trade-off between accuracy and the required computational resources.

Data Inspection
The models of shape dependant translation, rotation and deformation resistance tensor components are based on the conducted parametric study, where we investigate approximately 5.4 • 10 3 superellipsoidal particles, which are obtained from the following range for the superellipsoidal parameters: ] and e 2 = [0.2,1.8].Additionally, we employ a minor simplification by limiting the parameter λ 2 ≤ λ 1 , as we recognise that, for example, the particle with axial ratios of λ 1 = 5, λ 2 = 3 is in fact geometrically identical to the particle with axial ratios of λ 1 = 3, λ 2 = 5 rotated by 90 • .This in turn allows us to reduce the final count of the particles in the study from 9.8 • 10 3 to 5.4 • 10 3 , without loosing significant data in the set.That said, we compute each of the nine tensor components, required to construct the K, Ω and Π tensors, for the remaining range of particles in the study.At this point, we are left with a decent amount of raw data, obtained by the direct numerical computation with BEM.In order to simplify the process of putting the obtained results into use, we aim to find a suitable approximation scheme.
We have to realise that the mathematical space in consideration is essentially fourdimensional (λ 1 ,λ 2 ,e 1 ,e 2 ).As a result, we cannot expect to find any trivial solution for producing an accurate fit expression, to capture the complex nature of the data in question.Instead, we decide to make a number of data projections, for all nine tensor components, where we fix some of the parameters to constant values, which, in essence, reduces the dimensionality of the observation points.As an example of that, we present the tensor components K zz , Ω zz and Π zz on a 3D plot, where we fix two dimensions to a constant.In Figure 6 the tensor components are shown with respect to λ 1 and λ 2 , while e 1 and e 2 are fixed to unity.Likewise, in Figure 7 the data is displayed with respect to e 1 and e 2 , while λ 1 and λ 2 are 5 and 3 respectively.By observing the Figures 6 and 7, one can immediately establish that the data shows distinctively non-linear characteristics.We can maybe identify a slightly less pronounced non-linear behaviour in the case of K zz , whereas the behaviour of the Ω zz and Π zz components is relatively similar.With this in mind, we proceed by finding approximation functions, to reconstruct the resulting tensor components with respect to the particle geometric parameters.Initially, we limit ourselves to fit only a single variable function, while fixing other parameters to constant, due to the enormous complexity involved, when considering more than one parameter at once.By making a few educated guesses, we employ a general polynomial as a one-dimensional non-linear fit expression, which we fit to data observations, presented on Figures 6 and 7. Four functions are fitted to the results of each tensor component, namely two by fixing λ 1 and two by fixing λ 2 in Figure 6 and analogously with e 1 and e 2 in Figure 7.After comparing the fit accuracy for different order polynomials, we establish a third order polynomial as appropriate for fitting the parameters λ 1 and λ 2 and a second order polynomial for e 1 and e 2 .To confirm that, we evaluate the corresponding coefficient of determination R 2 for each fitted function, where we obtain R 2 = 0.98 as an average for all fitted functions.

Approximation Scheme Derivation
We are deriving a multivariate approximation scheme to predict each of the nine tensor components individually.We begin by first writing a general m − th order univariate polynomial, using vector notation, as p m (x) = [ 1 x x 2 ... x m ] T .At the moment we assume its coefficients to be equal to unity.We set up the initial expression for each geometric parameter, employing the polynomial order as just established, which yields p 3 (λ 1 ), p 3 (λ 2 ), p 2 (e 1 ) and p 2 (e 2 ).Prior to introducing the momentarily missing polynomial coefficients, we combine the first two and the second two polynomials together, to obtain two bi-variate polynomials.We write their tensor product and vectorize [49] its transpose back to the column vector, namely p 3,3 (λ 1 , λ 2 ) = vec p 3 (λ 1 )p 3 (λ 2 ) T T and p 2,2 (e 1 , e 2 ) = vec p 2 (e 1 )p 2 (e 2 ) T T .At this point we define the matrix of polynomial coefficients A, with dimensions 16 × 9, where 16 is the dimension of a vector p 3,3 (λ 1 , λ 2 ) and 9 is the dimension of a vector p 2,2 (e 1 , e 2 ).With everything set up so far, we write where f (λ 1 , λ 2 , e 1 , e 2 ) is the resulting approximation scheme for an individual tensor component.To obtain the coefficients of A we reformulate ( 22) into an optimization problem, again, for each tensor component separately.First, let us define the sum of squared errors, which we try to minimise, namely min where f i (λ 1 , λ 2 , e 1 , e 2 ) is the approximation scheme result (for the individual tensor component), evaluated for the i − th particle, y i is the tensor component of the i − th particle shape and n is the number of considered particle shapes.Considering the full range of the numerical results, we detect an increased standard deviation for Ω and Π results, which in turn increases the prediction difficulty for the approximation scheme.This increase of standard deviation in the results is mostly contributed by the exponentially rising magnitudes of the Ω and Π tensor components, with respect to λ 1 .To reduce this drawback as far as possible, we split the numerical results into two sub-ranges, to group together the values of similar magnitude and thus reduce the approximation error.By introducing a simple selection criteria, based on a single parameter λ 1 we obtain the two ranges, namely R 1 for λ 1 ≤ 5 and R 2 for λ 1 > 5.For solving the presented optimization problem, we use a multivariate non-linear regression method, from the software package GEKKO optimization suite [50].
As a result, we obtain the polynomial coefficient matrices with best-fitting terms for the individual tensor components in both ranges, namely The resulting coefficient matrices are given in the Appendix A. With that, we are able to determine all tensor components, for an arbitrary particle k within the range of the study, by simply evaluating (22) for λ 1 , λ 2 , e 1 , e 2 of the k − th particle, using A R (α) as the coefficient matrix, where index α denotes the tensor component and R the range under consideration.

Validation of the Approximation Scheme
Let us investigate the accuracy of the derived approximation schemes for the tensor components.Please note that there are particles in the parametric study, for which some of the tensor components are zero (or very close to zero, due to numerical errors).This in turn raises a question of representing the relative errors, as in those cases a relatively small absolute discrepancy would cause erratically high errors.For that purpose, we decide to use a modified relative error expression, which, on the example of K zz component, is given by where ε is a modified relative error, K * zz is the tensor component of the particle, as predicted by the approximation function, K zz is a numerical result of the tensor component and ||K || F is the Frobenius norm of the numerically computed K tensor.In Figures 8 and 9 we present the ε errors for K zz , Ω zz and Π zz tensor components, plotted in terms of λ 1 , λ 2 , and e 1 , e 2 respectively.If we first focus on Figure 8, we see that the distribution of K zz errors is fairly flat, with relatively low average error, indicated as ε(K zz ) in the right hand corner of the plot.On the contrary, the distributions for Ω zz and Π zz are showing a slightly increased dispersion between the values on the displayed range, with comparatively increased average errors.This is an expected behaviour, as we already identified that the standard deviation in the results of the latter two tensor components is much higher than in the results of the first.On the other hand, the plots in terms of e 1 , e 2 (Figure 9) are showing a very flat distribution for all considered tensor components.Following from that, we establish that the parameters λ 1 and λ 2 are significantly contributing to the standard deviation of the results and consequently to the increase of the approximation difficulty.This is also the main reason for the selection of the third order polynomial as the approximation function for λ 1 and λ 2 and the second order polynomial for e 1 and e 2 .Finally, to draw a conclusion from the above, we give the final values of modified relative errors, in terms of average for each tensor, namely ε(K) = 0.11%, ε(Ω) = 0.28% and ε(Π) = 0.31%.

Demonstration of the Force and Torque Model Using Parametrically Defined Particles
To evaluate the performance of the derived models on some common particle shapes, we put them to the test by comparing the results to various existing methods that have been extensively employed so far.For testing, we choose five particle shapes, namely the sphere (I), a prolate with λ 1 = 5 (I I), an ellipsoid with λ 1 = 5 and λ 2 = 3 (I I I), a rectangular disc superellipsoid with λ 1 = λ 2 = 5 and e 1 = e 2 = 0.2 (IV) and a pyramid-like superellipsoid with λ 1 = λ 2 = 10 and e 1 = e 2 = 1.5 (V).We start by predicting the translation resistance tensor K components.As already highlighted in the introduction, a simplification of the particle shape to a sphere is still commonly in use, so we first compare to the analytical model for the drag on a sphere.Another readily available analytical model that we include for comparison, is the drag on a prolate, first derived by Oberbeck.For more irregular, non-spherical particles, a number of empirical models, based on experimental results, are available.We use the models of Haider and Levenspiel [27], Leith [28] and Holzer and Sommerfeld [29] for our comparison.The results in terms of predicted K tensor components are given in Table 2.The shaded area in the last row of each particle section is the average of all component errors, computed by the modified relative error expression (24).For particles, where analytical results are available, we use those as a reference, while for other particles we use the results as obtained by the direct simulation.Obviously, particles I and I I are best captured by the corresponding analytical models, where the error is zero.On the other hand, as soon the particle geometry starts to deviate spherical or prolate shape, the errors of the analytical models start to increase significantly.In those cases, all three empirical models perform much better, where the most accurate off all three appears to be the one from Holzer and Sommerfeld.However, we can see that for all particles under consideration, the present model significantly outperforms all of the considered empirical models, while also performing very close to the accuracy of the analytical models, for the particles I and I I.The average errors obtained by the present model are of negligible magnitude, staying well below the order of 1%.
Table 2. K tensor coefficients are determined for particles I-V, using different methods.First, the analytical result is presented for the considered particle (applies only for particles I and I I).In the next column, a direct result from the present BEM code is presented.This is followed by the results obtained analytically with two simplified particle representations, namely the sphere and the prolate (considering λ 1 of the particle).After that, the empirical models from Haider and Levenspiel [27], Leith [28] and Holzer and Sommerfeld [29] follow.In the last column, the result of the derived approximation scheme is shown.

Particle a
Coeff.K

An.
Res.We continue with the prediction of the rotation resistance tensor Ω components.By analogy with the K tensor, we compare the results to other available models.One of the very few models available for the rotational dynamics of non-spherical particles is Jefferys analytical model for prolate ellipsoids.Despite our greatest efforts, we were not able to find any empirical models based on experiments that would be appropriate for modelling general non-spherical particles.Therefore, we use Jefferys model to produce the results for the sphere and for the prolate ellipsoid (using the particle λ 1 as the axial ratio).These are then used as a comparison for the effectiveness of the present approximation scheme, which is listed in Table 3.The fact that the availability of rotational dynamics models for non-spherical particles (other than prolate ellipsoids) is practically non-existent only further highlights the importance of the present model.Up until now, authors researching distinctively non-spherical particles, such as for example particles I I I, IV or V, have often neglected the rotational dynamics.The only practical alternative so far has been to simplify the particle shape to prolate and model it with the Jefferys model.Although this is still better than not considering the rotations at all, it is evident from Table 3 that the accuracy of such assumption quickly degrades when the actual particle shape migrates away from the prolate ellipsoid.For exactly this purpose, the advantage of the present model shows its true potential.One can observe that much like for the K tensor, the average prediction errors for all considered particles are well below the order of 1%.This also applies to predicting the sphere or the prolate particle, where the present model is able to almost perfectly reconstruct the analytical results.Table 3. Ω tensor coefficients are determined for particles I-V, using different methods.First, the analytical result is presented for the considered particle (applies only for particles I and I I).In the next column, a direct result from the present BEM code is presented.This is followed by the results obtained analytically with two simplified particle representations, namely the sphere and the prolate (considering λ 1 of the particle).In the last column, the result of the derived approximation scheme is shown.

Particle a
Coeff  As already mentioned, flow shear also affects the rotational dynamics of the particle, which we capture with the deformation resistance tensor Π.We predict its components in the same way as for the Ω tensor, where we also use the same comparison method.The results are presented in Table 4, from which we can draw a similar conclusion to Ω and K tensors, as far as the present model performance is concerned.Unlike for the Ω and K, in case of a Π tensor, some of its components are essentially zero, where the absolute error is few times greater than the target value.However, if we compare the magnitude of this error to the magnitude of the other non-zero tensor components, the ratio is actually negligible.From that we can assume that the practical consequences of this, in terms of particle rotational dynamics in the flow, would be essentially unnoticeable for any realistic engineering application.Table 4. Π tensor coefficients are determined for particles I-V, using different methods.First, the analytical result is presented for the considered particle (applies only for particles I and I I).In the next column, a direct result from the present BEM code is presented.This is followed by the results obtained analytically with two simplified particle representations, namely the sphere and the prolate (considering λ 1 of the particle).In the last column, the result of the derived approximation scheme is shown.

Demonstration of the Force and Torque Model Using a Realistic Particle
To show the applicability of the present model to a realistic engineering problem, we decide to demonstrate the process of predicting the translation, rotation and deformation resistance tensor components on an example of the realistic pollen particle [51].In Figure 10a we present a microscope image of the pollen grain under consideration.We begin by creating a 3D reconstruction of the particle, which was employed to the best of our abilities, using all of the available information from [51], the result of which is the surface shown in Figure 10b.To generate reference results, we employ the presented numerical model on the geometry of the reconstructed particle.We use the same set-up as established in the section of mesh study and model validation.Prior to obtaining the tensor components in question from our approximation scheme, we need to determine the four superellipsoidal parameters for the considered shape.We determine the bounding box of the particle, by finding its most-outer points.To obtain the optimal position in space, we rotate the particle so that we get the smallest possible bounding box.In the next step, we use bounding box edge lengths to identify the approximate axial ratios λ 1 and λ 2 and assume e 1 = e 2 = 1 to first obtain the ellipsoid shape that approximately fits the particle in question.The surface of an superellipsoid is given by (12), which yields unity for all points that lie on the surface.From that, we can form a numerical optimization problem, by writing an objective function [S(x, y, z) − 1] 2 , which we try to for the given set of parameters λ 1 , λ 2 , e 1 , e 2 .The resulting superellipsoidal particle, as obtained from the GEKKO optimization suite [50], has the parameters λ 1 = 1.96, λ 2 = 1.83, e 1 = 0.564, e 2 = 0.472 and is displayed in Figure 10c.We repeat the same process of determining the tensor coefficients as above, including the comparison with all available correlations.The results are presented in Table 5.
Table 5. K, Ω and Π tensor coefficients are determined for an example of the pollen grain.First the present BEM code result is presented, obtained for the actual reconstructed particle geometry, as illustrated in Figure 10b.This is followed by the results obtained analytically with two simplified particle representations, namely the sphere and the prolate (considering λ 1 of the particle).After that, the empirical models from Haider and Levenspiel [27], Leith [28] and Holzer and Sommerfeld [29] follow (applies only for the K tensor components).In the last column, the result of the approximation scheme is shown, for the superellipsoid (λ 1 = 1.96, λ 2 = 1.83, e 1 = 0.564, e 2 = 0.472), that best matches the pollen grain.

Coeff
If we focus initially on the translation resistance, we see again that the analytical models for the sphere or prolate produce quite poor results, which is unsurprising, since the particle shape is very different from the sphere or the prolate.As expected, the empirical correlations from [27][28][29] definitely provide some improvement to that.Much like with the results for particles I-V, the present approximation scheme offers the best replicate of the numerically obtained results, with almost perfect accuracy.In case of the rotation and deformation resistance tensors, the gap between the accuracy of currently available models and the present model is even more pronounced.We already established that the analytical model of Jeffery was essentially the only available approximation so far, which performs poorly for particles very different to prolates.Especially in this cases, the results of the present model, offer by far the greatest improvement.Overall we see that the average errors range from less than 1% for the translation, to ≈2% and ≈3% for the rotation and deformation resistance tensors respectively.

Conclusions
In the present work, we first describe the numerical framework for the prediction of forces and torques on non-spherical particles in Stokes flow.Due to the nature of the problem, we employ the BEM approach, which is proven to be superior in terms of computational cost and accuracy for the proposed problem set-up, as it enables a direct evaluation of the boundary tractions, leading to excellent accuracy of computed forces and torques acting on a particle.Additionally we also show the advantages of BEM for the problems where a large computational domain size is required, by including the comparison with similar approaches, employed by other authors.The numerical setup is described along with boundary conditions and the results of the mesh study and model validation are presented, where the final parameters for the parametric study are selected.This is followed by the presentation of the numerical study results, where a few examples of 3D data visualizations are given, on a limited range of parameters.The methods of data simplification are explained, where the initial, approximate relationships are first found on a limited data range, by fixing some of the parameters to constant, and thus temporarily reducing the dimensionality of the data.After that, the reconstruction of the simplified relationships into a functional multivariate approximation scheme is presented.Next, we demonstrate the reconstructed model functionality on various different particle examples, where we examine the model performance in terms of predicting the translation, rotation and deformation resistance tensors.The results of the present model are compared with different readily available models, where it is established that the present model delivers significantly improved results as compared to existing models.Furthermore, we identify that the existing models were only able to predict the translation resistance of the non-spherical particles with fairly acceptable accuracy.As it turns out, there is a significant lack of models that are able to account for the rotational dynamics of the non-spherical particles (other than prolates).As a result, researchers were often forced to neglect the rotational dynamics, or they had to rely on the simplification of the particle shape, which we demonstrate to be very inaccurate for most of complex particles.That said, we prove that the present model offers significant advantages, especially in modelling the rotational dynamics, as it can predict the rotation and deformation tensor components with essentially the same accuracy as for the translation resistance.As a final result, we have established that the average modified errors remain well below the order of 1% for all particles used for comparison.Even when trying to approximate a complex, nonsymmetrical pollen particle, with a best-fitting superellipsoid, we are able to reconstruct the tensor components with the accuracy in the order of 1% for the translation, 2% for the rotation and 3% for the deformation resistance.To the best of our knowledge, the present model is outperforming currently available models in terms of accuracy for all non-spherical particles (except for the sphere and the prolate, where the analytical models are obviously superior), within the considered parameter range.Moreover, it is also one of the first models to account for both, the translational and rotational dynamics, on such a wide range of complex, non-spherical particles.With the derivation of the approximation scheme it is also very simple to implement in practically any numerical Lagrangian particle tracking framework and it comes with a very moderate computational cost, as compared to some of the empirical model approaches, which require complicated projection area computation with every position change, relative to the flow.With this in mind, we are convinced that this is a significant step towards closing the gap in research, when it comes to modelling distinctively non-spherical particles.R 1 parameter limits:

Figure 1 .
Figure 1.Examples of superellipsoids with particle frame of reference defined.From the left a triaxial ellipsoid, oval, cylinder, cuboid, rhomboid and octahedron are shown.All shapes are formed from the same axes (a, b, c), by only altering the exponential parameters (e 1 , e 2 ).

Figure 2 .
Figure 2. Numerical domain in the shape of a sphere, with the diameter of the 1024d p .The particle under consideration is positioned in the centre of the domain, with its volume equivalent diameter d p .The boundary condition, on the surface of the outer sphere is a constant velocity field u, while a no slip condition u| p = 0 is imposed on the particle surface.

Figure 3 .
Figure 3. Boundary traction magnitude || q|| 2 at the surface of a rectangular disc particle, shown for the uniform flow (left), rotational flow (centre) and the shear flow (right).

Figure 4 .Figure 5 .
Figure 4. Surface meshes A, B, C and D on the example of rectangular disc particle.Mesh elements in the bottom right corner of the particle are displayed in magnified detail, to show the quality of the captured rounded edge.

Figure 6 .Figure 7 .
Figure 6.Tensor component values are presented with respect to λ 1 and λ 2 , for the examples of tensor components K zz , Ω zz and Π zz , shown at (a-c) sections respectively.Parameters e 1 and e 2 are fixed to unity.Note the triangular shape of the plot surface, due to the imposed limitation λ 2 ≤ λ 1 .Black and blue dotted lines are displaying the one-dimensional fit functions, which are positioned at the fixed values of λ 1 and λ 2 as indicated in the plot legend.

Figure 8 .Figure 9 .
Figure 8. Modified relative errors ε are presented, with respect to λ 1 and λ 2 , for the examples of tensor components K zz , Ω zz and Π zz , shown at (a-c) sections respectively.Parameters e 1 and e 2 are fixed to unity.The average error for the plot is displayed in the top right corner.

Figure 10 .
Figure10.grain in consideration for testing the performance of the present model.The picture of the original particle[51] is shown in (a), reprinted with permission from[51], copyright 2016 John Wiley and Sons.In (b) the reconstructed 3D surface of the particle is displayed and in (c) the best matching superellipsoid (λ 1 = 1.96, λ 2 = 1.83, e 1 = 0.564, e 2 = 0.472) is given.

Table 1 .
Comparison of the analytical and numerical values of the translation and rotation tensors in the case of a sphere and two prolate ellipsoids.