On the Analytical Solution of the Kuwabara-Type Particle-in-Cell Model for the Non-Axisymmetric Spheroidal Stokes Flow via the Papkovich–Neuber Representation

: Modern engineering technology often involves the physical application of heat and mass transfer. These processes are associated with the creeping motion of a relatively homogeneous swarm of small particles, where the spheroidal geometry represents the shape of the embedded particles within such aggregates. Here, the steady Stokes ﬂow of an incompressible, viscous ﬂuid through an assemblage of particles, at low Reynolds numbers, is studied by employing a particle-in-cell model. The mathematical formulation adopts the Kuwabara-type assumption, according to which each spheroidal particle is stationary and it is surrounded by a confocal spheroid that creates a ﬂuid envelope, in which the Newtonian ﬂuid moves with a constant velocity of arbitrary orientation. The boundary value problem in the ﬂuid envelope is solved by imposing non-slip conditions on the surface of the spheroid, which is also considered as non-penetrable, while zero vorticity is assumed on the ﬁctitious spheroidal boundary along with a uniform approaching velocity. The three-dimensional ﬂow ﬁelds are calculated analytically for the ﬁrst time, in the spheroidal geometry, by virtue of the Papkovich–Neuber representation. Through this, the velocity and the total pressure ﬁelds are provided in terms of a vector and the scalar spheroidal harmonic potentials, which enables the thorough study of the relevant physical characteristics of the ﬂow ﬁelds. The newly obtained analytical expressions generalize to any direction with the existing results holding for the asymmetrical case, which were obtained with the aid of a stream function. These can be employed for the calculation of quantities of physical or engineering interest. Numerical implementation reveals the ﬂow behavior within the ﬂuid envelope for different geometrical cell characteristics and for the arbitrarily-assumed velocity ﬁeld, thus reﬂecting the different ﬂow/porous media situations. Sample calculations show the excellent agreement of the obtained results with those available for special geometrical cases. All of these ﬁndings demonstrate the usefulness of the proposed method and the powerfulness of the obtained analytical expansions.


Introduction
Flow around particles and fluid-particle interactions are encountered in many processes in physics, engineering, biology and medicine, such as sedimentation, fluidization, flow in granular porous media, flow of colloids of dilute emulsions through porous media and suspensions of living cells. The analytical determination of the flow field in these cases is important, as it enables the calculation of physical quantities of engineering interest, such as the vorticity field, the drag force acting on each particle and the macroscopic pressure gradient. Moreover, they can be used for the modelling of heat and mass transport phenomena, such as adsorption, diffusion and reaction. When the viscous characteristics of self-sufficiency of the model in regard to mechanical energy and it allows the exchange of mechanical energy between the cell and the environment. These differences have to be taken into account in modelling physical and engineering problems.
A more realistic approach to the study of flow through porous media was proposed by Neale and Nader [11], who considered the unit cell to be encapsulated within an unbounded porous medium, and assumed continuous, homogeneous and isotropic permeability that is identical to that of the swarm of spheres. Slip boundary conditions on the unit cell have been proposed in relation to the flow of micropolar fluids. Faltas and Saad [12] investigated Stokes flow past an assemblage of slip eccentric spherical-in-cell models, while Sherief et al. studied the Stokes flow of a micropolar fluid past an assemblage of spheroidal-in-cell model with slip boundary conditions [13] and the interaction between two rigid spheres [14]. Recently, Madasu [15] studied the problem of an incompressible axisymmetric creeping flow caused by a porous spherical particle in a spherical cavity filled with micropolar fluid with both Happel and Kuwabara conditions.
Due to the difficulties encountered in the analytical derivation of Stokes flow for non-spherical geometries, the spherical in-cell model remained the preferred model for almost 150 years, yielding a closed-form analytical expansion for the stream function and the pressure field, while problems in non-spherical geometries were mainly treated numerically. For example, Epstein and Masliyah [16] considered spheroidal-shaped particles and solved the creeping flow problem in an unbounded porous medium numerically, since no analytical expression could be derived for the stream function in the spheroidal coordinates. Dassios et al. [5,17] derived generalized eigenfunctions and introduced complete semiseparable solutions for the stream function in the prolate and the oblate spheroidal geometry expanding the spherical particle-in-cell model to axisymmetric geometries. The solution for the Stokes flow problem through the swarm of particles was obtained by applying either Kuwabara or Happel boundary conditions. Recently, Hadjinicolaou and Protopapas [18], used an inversion method to derive the analytical solution for the Stokes flow within two inverted spheroidal cells and modelled the blood plasma flow through red blood cells. Their solution employed the R-semi-separation method, which they had introduced previously [19].
Furthermore, Vafeas and Dassios [20], by employing the Papkovich-Neuber differential representation as it applies to the Stokes flow and also a Happel-type ellipsoid-in-cell model, obtained the velocity field and total pressure field in terms of harmonic ellipsoidal eigenfunctions. Dassios et al. [21] demonstrated how Papkovich-Neuber representations of the flow fields were correlated to the general solution of Stokes equations in spheroidal geometry. Furthermore, Dassios and Vafeas [22] corroborated the powerfulness and usability of the above 3D Papkovich-Neuber representation by obtaining a closed-form solution for a Happel-type 3D Stokes flow through a swarm of translating and rotating spherical particles.
Numerous non-analytical methods can be found in the literature and provide solutions for the flow problem at hand, including Stokesian dynamics [23] or lattice-Boltzmann simulation [24]. These are beyond the scope of the present research and thus are not discussed further. The analytical treatment of the flow problems is of great significance since it enables the analysis of the geometrical and physical characteristics of the flow with Kuwabara-type boundary conditions [25][26][27] and reveals features of the transport process under consideration, e.g., see for instance [28,29]. Specifically, Djedaidi and Nouri [29] provided theoretical evidence for the existence and uniqueness of flow interactions in heterogeneous porous media without using regularization techniques.
In the present manuscript, we derived for the first time, an analytical solution of a Stokes flow problem assuming a uniform, induced velocity field in a prolate spheroidal fluid cell, arbitrarily orientated, thus overcoming any limitations imposed by the axisymmetric assumptions. Therefore, we used a particle-in-cell model and the Papkovich-Neuber representation instead of the semi-separation method, which provided useful analytical expansions for the velocity and the total pressure. This is a generalization of the results obtained in [17] since we derived the analytical solution representing the flow of a swarm of prolate particles in 3D instead of the corresponding axisymmetric creeping flow, which represents a 2D flow. The special cases of the needle, the disk, the oblate spheroid and the sphere can be derived from the finalized results by choosing the appropriate values of the spheroidal parameters. The prolate spheroidal coordinate system and its analysis [30] are introduced, and fit our problem perfectly, while the results for the oblate geometry can be easily obtained via a simple transformation. The flow fields assume compact closed-form formulae, given as series expansions in terms of prolate spheroidal harmonic eigenfunctions [31] and the analytical results are subjected to the necessary numerical manipulation, in order to provide plots for the implicated fields.
The structure of our manuscript is as follows: in Section 2 we present the fundamental equations and the mathematical development, while in Section 3 we formulate the boundary value problem and derive the analytical solution. In Section 4, we illustrate the obtained results through a variety of graphs of the velocity field, the associate stream lines and the pressure field. In Section 5, a thorough discussion of the research outcomes is given.

Fundamental Equations and Mathematical Development
The governing equations of the steady state Stokes flow of an incompressible, viscous fluid, around particles encapsulated within a smooth, domain Ω R 3 , which is either bounded or unbounded, characterized by dynamic viscosity µ and mass density ρ, are given by the well-known system of partial differential equations [3], which connect the biharmonic velocity field v with the harmonic total pressure field P and assume the form and ∇ · v(r) = 0 for r ∈ Ω R 3 Given the velocity v, the vorticity renders ω = ∇ × v, which reveals that ω is also a harmonic function. Equation (1) states that, in Stokes flow, at any material point of the fluid, there is a balance between the force caused by the pressure gradient and the viscous force. Equation (2) denotes that the velocity field is solenoidal, assuring the incompressibility of the fluid.
Given a fixed positive number c > 0, which specifies the semi-focal distance of our system and by virtue of the variables τ ∈ [1, +∞), ζ ∈ [−1, 1] and ϕ ∈ [0, 2π), we define the transformed prolate spheroidal coordinates [30] as specifying any material point r = x 1x1 + x 2x2 + x 3x3 . The outward unit normal vector on the surface of any spheroid τ = constant is given via the Cartesian basis bŷ where the expression τ 2 − ζ 2 , appearing at the denominator of Equation (4) is always bounded away from zero for all points on or exterior to the surface of the prolate spheroid. Papkovich-Neuber [2] proposed the following three-dimensional differential representation of the solution for Stokes flow, in terms of the harmonic potentials Φ and Ψ, i.e.,

v(r)
and with ∆Φ(r) = 0 and ∆Ψ(r) = 0 for r ∈ Ω R 3 (7) where in prolate spheroidal coordinates, the involved differential operators and stand for the gradient and the Laplacian, respectively, whileτ (which matchesn from Equation (4)),ζ andφ denote the unit coordinate vectors of the system.
In view of the definition of the associated Legendre functions of the first P m n and of the second Q m n kind [31], we introduce the regular on the axis of symmetry (at ζ = ±1) of the prolate spheroid interior and exterior harmonic eigenfunctions of degree n ≥ 0 and of order m = 0, 1, 2, . . . , n u m(q) respectively, wherein for n ≥ 0, denoting the even (e) and the odd (o) parts of the potentials. The functions Y m(q) n play the role of the so-called surface spherical harmonics [31], which are orthogonal with respect to the normalized integral for every n, n ≥ 0, (m, m ) = 0, 1, 2, . . . , (n, n ) and q, q = e, o, where ε 0 = 1 and ε m = 2 otherwise, while the symbol δ denotes the Kronecker delta. Since the functions Y m(q) n form an orthonormal basis according to the inner product of Equation (12), then any smooth and well-behaved function h(ζ, ϕ) can be expanded via where, multiplying Equation (13) by Y m (q ) n and integrating both parts by utilizing Equation (12), the leading coefficients in Equation (13) admit for n ≥ 0, m = 0, 1, 2, . . . , n and q = e, o, the latter being extensively used in the sequel. The implicated harmonic functions to our problem, i.e., the potentials Φ and Ψ that satisfy Equation (7) can be expanded in terms of the harmonic-type eigenfunctions (10) such as n,ext u m(q) n,ext (r) for r ∈ Ω R 3 (15) and n,ext u m(q) n,ext (r) for r ∈ Ω R 3 (16) where the vector constant coefficients that correspond to the vector potential Φ are written both in the Cartesian and the prolate spheroidal geometry via The implicated harmonic functions to our problem, i.e., the potentials Φ that satisfy Equation (7) can be expanded in terms of the harmonic-type eigenfu (10) , m q n p d stands for the scalar constant coefficients that refer to the scalar p Ψ, all sets of constant coefficients being unknown and given for 0 n,p stands for the scalar constant coefficients that refer to the scalar potential Ψ, all sets of constant coefficients being unknown and given for n ≥ 0, m = 0, 1, 2, . . . , n and q = e, o with p = int, ext. Hence, by virtue of Equations (15) and (16) and using classical vector identities [3,31], the flow fields in Equations (5) and (6) (15) and (16) and using classical vector identities [3,31], the flow fields in Equations (5) and (6) or equivalently, as far as the velocity field is concerned, and similarly, giving the total pressure field. Once the velocity field is obtained, the vorticity field is as far as the velocity field is concerned, and similarly, giving the total pressure field. Once the velocity field is obtained, the vorticity field is where, in prolate geometry, the operator ∇ is given by (8), while the compact expressions of the fields are the outcome of a handy redefinition of the eigenfunctions in Equation (10) via the formula (22) in which n ≥ 0, m = 0, 1, 2, . . . , n and q = e, o with p = int, ext. On the other hand, due to Equations (3) and (17), it is easily confirmed that while Equations (8) and (17) are combined to give n,p cos ϕ ∂ ∂ϕ (24) and metry 2022, 14, 170 8 of 21 and 1 cos sin comprising useful formulae, defined for

Boundary Value Problem and Analytical Solution
For the mathematical formulation of the physical problem, two confocal prolate spheroids are considered with semi-focal distance c , whose major axis is parallel to the comprising useful formulae, defined for n ≥ 0, m = 0, 1, 2, . . . , n and q = e, o with p = int, ext.

Boundary Value Problem and Analytical Solution
For the mathematical formulation of the physical problem, two confocal prolate spheroids are considered with semi-focal distance c, whose major axis is parallel to the x 3 -axis of a Cartesian coordinate system. The surface of the inner spheroid is solid and stationary, indicated by S α at τ = τ α = a 3 /c, while the fluid layer is confined by the outer fictitious spheroidal surface, indicated by S β at τ = τ β = b 3 /c, where a 3 and b 3 are the corresponding major semi-axes of the two confocal spheroids. Obviously a 1 = a 2 and b 1 = b 2 are the equal minor semi-axes of the inner and the outer spheroid, respectively, which are lying on the x 1 x 2 -plane. An arbitrarily orientated and uniformly approaching velocity generates the non-axisymmetric flow in the fluid layer between the two prolate spheroidal surfaces. Therefore, the domain of interest is now defined by The applied boundary conditions are similar to those in the Kuwabara-type particlein-cell model [10], which yield (29) and Equation (28) expresses the non-slip flow condition, while the non-homogeneous relationship in Equation (29) implies that there is a flow across the boundary of the fluid envelope S β , the second equality being the outcome of the inner product of Equation (4) with Equation (26), resulting in the known function, g(ζ, ϕ). Furthermore, in compliance to Kuwabara's model, the vorticity on the external spheroid is assumed to vanish, as it is shown in Equation (30). This completes the formulation of a well-defined boundary value problem.
At this point, we must to make an important note. As is easily observed from Equations (19)- (21), in order to identify them, we have to calculate four sets of scalar coefficients from the application of the boundary condition on the interior spheroidal surface and four sets of scalar coefficients from the application of the boundary condition on the exterior spheroidal surface, which makes eight sets of unknowns in total. However, the boundary conditions from Equations (28)-(30) offer us three, plus one, plus three, that is, seven scalar relations, meaning that we need one more condition so that the problem is wellposed. This disadvantage is cancelled due to the flexibility of the differential representation theory. In detail, since the Papkovich-Neuber general solution uses four potentials (one vector and one scalar) to describe the three-component vector velocity field, we are left with a certain degree of freedom in selecting the unknown constant coefficients that correspond to the potentials, which will provide us with the eighth missing condition, conveniently chosen according to physics or/and mathematics (e.g., see [24] for the applicability of this aspect).
Analogous results for the oblate spheroidal geometry [31] can be easily obtained through the simple transformation τ → iλ and c → −ic (31) where 0 ≤ λ ≡ sinhη < +∞ and c > 0 are the characteristic variables of the oblate spheroidal system. On the other hand, the limiting case of the needle can be recovered by a prolate spheroid, where 0 < a 2 = a 1 a 3 < +∞ and 0 < b 2 = b 1 b 3 < +∞, while in the case where 0 < a 2 a 1 = a 3 < +∞ and 0 < b 2 b 1 = b 3 < +∞, the oblate spheroid becomes a circular disk. Finally, the sphere corresponds to the values a 1 = a 2 = a 3 = a and where a, b > 0 are the radii of the concentric spheres that correspond to the Kuwabara spherical particle-in-cell model.
In order to mathematically manipulate the boundary conditions from Equations (28)- (30) in an easy manner, we project the vector velocity and vorticity fields from Equations (18) or (19) and (21) to the prolate spheroidal basis via and respectively, while the pressure field is given by the expansion of Equation (20), all the above relations are valid in view of Equations (17) and (23)- (25). Thus, taking into account the results of Equations (32) and (33), we proceed to the application of the pre-mentioned boundary conditions, which provides us with seven independent scalar conditions, particularly, three from Equation (28) for k = 1, 2, 3, one from Equation (29) for k = 4 and three from Equation (30) for k = 5, 6, 7 that have to be supplemented by the missing eighth fictitious condition for k = 8, according to the flexibility of the Papkovich-Neuber differential representation, as discussed earlier. However, this condition must be appropriately chosen with respect to mathematical or/and physical explanation, in order to obtain accurate and validated results. To this end and for reasons of clarity and writing convenience, we primarily proceed to an alternative convenient definition of the unknown constant coefficients, such as for every n ≥ 0, m = 0, 1, 2, . . . , n and q = e, o Consequently, we obtain the compact form of the seven sets of boundary conditions from Equations (28)- (30), plus the extra eighth condition (i.e., for k = 1, 2, . . . , 8), those being where δ k,4 denotes the Kronecker delta (δ k,4 = 1 when k = 4, otherwise zero), which comprise standard homogeneous conditions, except the non-homogeneous fourth one that corresponds to Equation (29) and the eighth one that implies the appropriately imposed relationship between the unknown constant coefficients. The functions f m(q),k n,j p within Equation (35) for n ≥ 0, m = 0, 1, 2, . . . , n and q = e, o when k, j p = 1, 2, . . . , 8, also according to Equation (22), yield f m(q),1 n,j p f m(q),1 n,j p for k = 3, as far as Equation (28) is concerned, while for k = 4, matching the non-homogeneous condition of Equation (29). On the other hand, the rest of the leading functions imply f m(q),5 n,j p =1,2 (ζ, ϕ) = τ 2 β − 1R m n,p τ β ζP m n (ζ) cos ϕ f for k = 6, f m(q),7 n,j p =1,2 (ζ, ϕ) = τ β R m n,p τ β 1 − ζ 2 P m n (ζ) − τ 2 β − 1 R m n,p τ β ζP m n (ζ) cos ϕ f n,j p =5,6 (ζ, ϕ) = τ 2 β − 1 1 − ζ 2 R m n,p τ β ζP m n (ζ) + τ β R m n,p τ β P m n (ζ) f for k = 7, which concern the vanishing vorticity condition of Equation (30). Finally, the key method to the solution is the arbitrary choice of the final functions for k = 8, those being whenever k, j p = 1, 2, . . . , 8 and where the leading constants F m/m (q/q ),k n/n ,j p and G m (q ) n of the summations for n, n ≥ 0, (m, m ) = 0, 1, 2, . . . , (n, n ) and q, q = e, o with k, j p = 1, 2, . . . , 8, by virtue of Equation (14), enjoy the expressions F m/m (q/q ),k n/n ,j p respectively. In doing so, we manage to transfer the difficulty of solving Equation ( where Equation (66) stands for eight separate conditions for each k = 1, 2, . . . , 8. Then, Equation (66) represents a system of linear algebraic equations, which can be solved using cut-off techniques by imposing an indispensable common upper limit for both of the twodegree indexes, let us assume n = n = 0, 1, 2, . . . , N, in order to obtain quadratic forms, such as where for m = m = 0, 1, 2, . . . , n and q = q = e, o, we obtain

Numerical Implementation
In this section, we apply the obtained results to a particular spheroidal cell in order to demonstrate the three-dimensional vector velocity field, its two-dimensional intersectional mapping on the coordinate planes and the pressure field on every axis. The numerical implementation is based upon the use of the Cartesian coordinate system (x 1 , x 2 , x 3 ), so that the velocity field of Equation (19) is translated in Cartesian geometry as x 3 ) and the total pressure field in Equation (20) is assumed to yield P = P(x 1 , x 2 , x 3 ), all these expressions use the coefficients derived from Equation (66) and by virtue of the definitions from Equation (34).
Additionally, in what follows, each field or physical property is provided in SI units, which we do not show for reasons of writing convenience and clarity.
In order to represent the swarm of particles, we choose the solid volume fraction γ of the cell to be equal to that of the swarm. We recall that both of the surfaces S α and S β are confocal and c is their common semi-focal distance. As stated previously, we denote a 3 as the long and a 1 = 1 as the short semi-axis of the solid internal prolate spheroid, and b 3 as the long and b 1 as the short semi-axis of the outer prolate spheroid [5], which satisfies Therefore, solving Equation (70) with respect to b 3 , we arrive at where, which provides the necessary prerequisites for the proceeding numerical elaboration of the produced results.
To plot the graphs of the fields, we have to calculate the constant coefficients, given in Equation (34), of the series expansion of the velocity and the pressure field. We solve the system (Equation (67)) using cut-off techniques assuming N = 1, which leads to a linear quadratic system of 48 unknowns and equations. Furthermore, without loss of generality and for computational convenience, we assume a simple linear form for φ m(q) n,j p (ζ, ϕ) = ϕ with j p = 1, 2, . . . , 8, where ϕ ∈ [0, 2π] and U = (1, 2, 3) are the arbitrarily orientated and uniformly approaching velocity, while µ = 1 stands for the dynamic viscosity. On that account, each graph that follows (see for instance Figures 1-11 in the sequel), is drawn by taking a 3 = 3 and γ = 0.05, so that τ α = 1.1, τ β = 1.6 and c = 2.8.
Therein, the velocity field and the streamlines are plotted in the bounded region, wherein x 1 , x 2 , x 3 ∈ [−5, 5]. In particular, in Figure 1, we demonstrate the velocity field in a 3D frame and in Figures 2-10, we present pairs of graphs (the vector plot and the corresponding streamlines) for all the combinations of any two of the velocity components v i , v j for i, j = 1, 2, 3 and i < j on each coordinate plane x 2 x 3 (for x 1 = 0), x 1 x 3 (for x 2 = 0) and x 1 x 2 (for x 3 = 0). It is worth noting that the length of each vector in the components of the velocity field indicates the strength of the field in the particular position.                      , v v in the elliptical ring on the coordinate plane 2 3 x x for 1 0 x = .        , v v in the elliptical ring on the coordinate plane 1 3 x x for 2 0 x = . Figure 10. Vector plot (a) and streamlines (b) of velocity components 2 3 , v v in the circular ring on the coordinate plane 1 2 x x for 3 0 x = . Figure 11. Total pressure field Ρ on 1 x -axis (a), on 2 x -axis (b) and on 3 x -axis (c). , v v in the elliptical ring on the coordinate plane 1 3 x x for 2 0 x = . Figure 10. Vector plot (a) and streamlines (b) of velocity components 2 3 , v v in the circular ring on the coordinate plane 1 2 x x for 3 0 x = . Figure 11. Total pressure field Ρ on 1 x -axis (a), on 2 x -axis (b) and on 3 x -axis (c). , v v in the elliptical ring on the coordinate plane 1 3 x x for 2 0 x = . Figure 10. Vector plot (a) and streamlines (b) of velocity components 2 3 , v v in the circular ring on the coordinate plane 1 2 x x for 3 0 x = . Figure 11. Total pressure field Ρ on 1 x -axis (a), on 2 x -axis (b) and on 3 x -axis (c). Initially, taking into account that the approaching velocity that we used was U = (1, 2, 3), we can observe that the three-dimensional velocity field in Figure 1 is affected by the imposed main velocity and it is orientated in a similar direction, which is also displayed in the two-dimensional intersections of Figures 2-10. In the sequel, by analyzing the vector plots given in Figures 2a, 3a, 4a, 5a, 6a, 7a, 8a, 9a and 10a, we can verify the conclusion obtained from Figure 1 that the choice of the approaching velocity affects the velocity field, since in most of these figures in which x 3 -axis is employed, the velocity vectors are directed towards the x 3 -axis. More precisely, from Figures 2a, 3a and 4a, we see that the vector field created by the vectors v 1 , v 2 is symmetric in the fluid envelope area because the velocity component v 3 is bigger than the others, so it dominates in the velocity field. On the other hand, comparing Figures 5a and 6a, we see that the vector components v 1 , v 3 create a field, which is strong for the positive values of x 3 , a result that is valid for the vector components v 2 , v 3 , given in Figure 8a. On the contrary, the velocity field for v 2 , v 3 in Figure 9a shows that the vectors are of almost equal measure in the plane x 2 = 0. Moreover, from Figures 2a, 5a and 8a, we can observe the velocity fields that occur for two of the velocity's components in the plane x 2 x 3 and we deduce that the velocity component v 2 is the one that contributes the least. The velocity components v i , v j for i, j = 1, 2, 3 and i < j in the plane x 1 x 3 are presented in Figures 3a, 6a and 9a, and show that the vectors of equal measure in most cases, except for the negative values of x 3 for v 1 , v 3 in Figure 6a, where the measures are smaller. In the plane x 3 = 0, we see that the field created by v 1 , v 2 is axisymmetric with respect to the x 2 -axis (Figure 4a), while the field for the components v 2 , v 3 and v 1 , v 3 , drawn in Figures 7a and 10a are symmetric with respect to (0, 0). Finally, Figure 11 depicts the variation in the total pressure field with respect to one of the variables x 1 , x 2 , x 3 , setting the other two equal to zero, inside the envelope's region, where x 1 , x 2 ∈ (−3.49, −1.29) ∪ (1.29, 3.49) and x 3 ∈ (−4.48, −3.1) ∪ (3.1, 4.48). Particularly, from Figure 11a,b we can conclude that the absolute value of the pressure decreases when we are moving away from the solid prolate spheroid into the fluid envelope along the x 1 -axis or x 2 -axis, where the pressure fields are axisymmetric. Moreover, the pressure of the fluid envelope is negative on x 1 -axis and it is approximately zero on the boundary of the fictitious prolate, while it is positive on the x 2 -axis and it is approximately one on the boundary of the fictitious prolate. In Figure 11c, we see that the pressure in the x 3 -axis in not axisymmetric and it is increasing in the lowest area of the fluid envelope (from the exterior to the interior prolate spheroid), while in the upper area, the pressure is increasing from the interior to the exterior prolate spheroid. In both cases, the pressure increases from negative values.

Conclusions and Discussion
An analytical methodology for dealing with the problem of Stokes flow through a swarm of spheroidal particles with Kuwabara-type boundary conditions was developed. For the problem at hand, the prolate spheroidal configuration was considered, while the results for the oblate spheroidal case can be readily obtained via a simple complex transformation. Moreover, interesting degenerate cases of the needle, disk and the limiting case of the sphere can be derived just by setting suitable values to the geometrical parameters that enter into the flow fields' expressions. For example, the case of the sphere is taken as the semifocal distance c tends to zero and the semi-axis ratio α 3 /α 1 goes towards unity, while the needle case is taken as the semi-axis ratio becomes much bigger than one. The conducted research was based on the reliable particle-in-cell technique by which the creeping flow through the swarm of the spheroidal particles is modelled as flow within a fluid envelope that surrounds a single solid particle and attains a fictitious spheroidal boundary, which expresses the disturbance to the flow that causes the rest of the particles of the assemblage and the volume of the fluid envelope to match the fluid volume fraction of the swarm. Hence, in view of the Kuwabara model, which assumes an arbitrarily moving fluid around a motionless particle, the non-slip flow condition on the surface of the inner spheroid was supplemented by the boundary conditions on the outer spheroidal surface, i.e., the zero normal velocity component and vorticity.
Under the low Reynolds hydrodynamic consideration, the fundamental differential representation of Papkovich-Neuber for Stokes flow was used and the non-axisymmetric flow fields were obtained in terms of easy-to-handle harmonic functions. The difficulty of this problem emerges in the determination of the unknown coefficients, which reveal inherent characteristics of the corresponding potentials. This difficulty was caused by the appearance of particular indeterminacies depicting the complexity of the spheroidal geometry. On that account, the flexibility of the particular general solution was demonstrated by imposing appropriate additional conditions, which made the calculations feasible. In particular, the aforementioned boundary conditions, enriched by the additional set of equations that were imposed due the degrees of freedom that the Papkovich-Neuber representation offered us, led to a set of equations, defined on the surfaces of the two spheroids. Therein, the application of standard orthogonality arguments created a linear system for the unknown constant coefficients, which were solved incrementally via cut-off techniques. Evidently, the application of the representation theory, combined with the imposed conditions, offered closed-form expressions for the flow fields under consideration as series expansions of prolate spheroidal harmonic eigenfunctions. The analytical formulae of the velocity and the pressure field were numerically implemented using cut-off techniques. For computational convenience, N was set equal to 1 and all the arbitrary functions φ m(q) n,j p (ζ, ϕ) with j p = 1, 2, . . . , 8 given in (61), were set equal to ϕ.
Consequently, it is obvious that advanced analytical or semi-analytical solutions in fluid dynamics are quite valuable tools for the assessment of the validity of numerical methods and solutions. In addition, bearing in mind that very important physical laws are expressed through analytical formulation, the necessity of tackling problems with a credible mathematical basis is obvious. In our view, there is always a need for these types of analytical methods, which can coexist with pure numerical techniques that aim to provide solutions to problems arising in the aforementioned physical applications. In this respect, work in progress involves the utilization of the differential representation theory in solving such non-axisymmetric boundary value problems in more complicated geometries, such as ellipsoidal, bi-spherical, toroidal, etc., wherein randomly-shaped particles in porous media can be depicted more effectively.