A specific slip length model for the Maxwell slip boundary conditions in the Navier–Stokes solution of flow around a microparticle in the no-slip and slip flow regimes

In the case of microscopic particles, the momentum exchange between the particle and the gas flow starts to deviate from the standard macroscopic particle case, i.e. the no-slip case, with slip flow occurring in the case of low to moderate particle Knudsen numbers. In order to derive new drag force models that are valid also in the slip flow regime for the case of non-spherical particles of arbitrary shapes using computational fluid dynamics, the no-slip conditions at the particle surface have to be modified in order to account for the velocity slip at the surface, mostly in the form of the Maxwell’s slip model. To allow a continuous transition in the boundary condition at the wall from the no-slip case to the slip cases for various Knudsen (Kn) number value flow regimes, a novel specific slip length model for the use with the Maxwell boundary conditions is proposed. The model is derived based on the data from the published experimental studies on spherical microparticle drag force correlations and Cunningham-based slip correction factors at standard conditions and uses a detailed CFD study on microparticle fluid dynamics to determine the correct values of the specific slip length at selected Kn number conditions. The obtained data on specific slip length are correlated using a polynomial function, resulting in the specific slip length model for the no-slip and slip flow regimes that can be applied to arbitrary convex particle shapes.

of a few micrometres and below can significantly contribute to the spreading of a virus, especially indoors [1]. In order to assess the danger of infection in various cases and flow conditions as well as to study the transport of aerosols in the human respiratory system, computational analysis based on numerical approaches to particle tracking in fluid flows offers a fast, safe and parametric variation-friendly solution [2]. Here, the Langrangian particle tracking in combination with computational fluid dynamics (CFD) solution of the Navier-Stokes equations is the method of choice [3,4]. In numerical simulations of dispersed multiphase flows by means of Lagrangian particle tracking methods, the correct choice of fluid dynamic force models is one of the main challenges. In the case of microparticles, the point particle approach is predominantly applied, especially in tracking large numbers of particles, due to its comparatively low computational demands [5]. With the point particle approach, the shape of a particle is not resolved; hence, the coupling between the flow field and the particle is achieved by means of dedicated force models coupled with empirical correlations of the model parameters. Among the force models, the drag force model is the most important one, as the drag force is the predominant force exerted by a fluid on a particle. In the case of microparticles, the Reynolds numbers are typically well below one and the flow around the particle can be described by Stokes flow [6]. For a sphere with Re p << 1, further denoted as microsphere, in continuum flow the main model parameter, i.e. the drag coefficient is related to the particle Reynolds number [7]. However, as the particle size decreases, molecular effects start to impact the particle-fluid interaction [8], and the additional effect of rarefaction has to be included in the form of a drag model dependence on the Knudsen number Kn. The Knudsen number, defined as the ratio of the gas mean free path λ and a representative physical length scale (such as a gap length over which mass or thermal transport occurs), serves as a means of quantification of the molecular effects in the fluid flow. In the case when the particle size becomes comparable to the mean free path λ of the molecules in the surrounding gas, the no-slip conditions on the particle surface no longer reflect the real conditions. Based on Schaaf and Chambre [9], the flow regimes can be roughly divided into three categories: • Continuum regime (Kn < 0.01), where the continuum assumption holds and the Navier-Stokes (N-S) equations with no-slip boundary conditions are applied in numerical solutions of gas particle flows. • Slip flow regime (0.01 < Kn < 0.1), where the velocity at the particle boundary no longer satisfies the no-slip conditions; hence, the slip flow regime is observed. The fluid flow can still be resolved by solving the Navier-Stokes equations by applying the slip velocity boundary conditions. • Transitional regime (0.1 < Kn < 10), where the continuum theory together with slip conditions begin to break down. • Free molecular regime (Kn > 10), where the continuum theory is no longer valid. When modelling the gas particle flow by the point particle approach, the drag force model requires a correct specification of the drag coefficient. In the low Knudsen number regime, the drag coefficient typically relates to the value of the particle Reynolds number, whereas in the slip flow regime the drag coefficient typically relates to the Knudsen number and a slip correction parameter is introduced. The determination of the slip correction parameter for small particles was covered experimentally by several authors, starting with Knudsen and Weber [10]. Today, there exist several excellent experimental studies on the determination of the slip correction parameter for micro-and nano-spheres, see [11][12][13][14][15][16].
On the other hand, the determination of the slip correction parameter can also be achieved by performing high resolution numerical simulation, using CFD. The CFD solves the Navier-Stokes equations, either with no-slip or slip boundary conditions, leading to a flow field determination around particles of arbitrary shape. The computed pressure and viscous stress values at the particle surface enable straightforward calculation of the resulting drag force, and in turn, when the slip flow conditions are applied, also derivation of the slip correction parameter values. For the continuum flow regime, implementing the no-slip conditions is straightforward [17].
In the case of the slip flow regime, the setting of correct slip conditions is not trivial, as the slip depends on the gas mean free path length as well as on the particle dimensions. Typically, the slip boundary condition is related to the fictitious macroscopic slip velocity used in the context of standard Navier-Stokes equations. A first-order slip boundary condition relates the slip velocity to the magnitude of the shear stress at the wall, i.e. the normal derivative of the tangential velocity. Here, the Maxwell model and the model of Schaaf and Chambre [9] are the most prominent examples. In the latter, the tangential momentum accommodation coefficient σ is used, which is the most commonly employed empirical parameter to illustrate overall gas-surface interaction. Note that the tangential momentum accommodation coefficient is defined in the range of σ = [0; 1], where σ = 1 describes diffuse scattering and σ = 0 specular reflection. With the same thermodynamic conditions, i.e. a constant value of the gas mean free path length, the slip conditions in the case of direct CFD computations vary for different particle sizes. In Sun et al. [18], a constant value of σ = 0.8 for the whole tested Knudsen range (0 < Kn < 0.15) was used (but not substantiated) for the direct computation of drag coefficient, leading to some deviations in the obtained results, compared to the reference experimental data. It was further shown [19] that in the case of curved solid wall surfaces the slip velocity should not be related solely to the normal derivative of the tangential velocity component, but extended and related to the viscous shear stress at the solid wall surface. Also, in order to include the increase in the rarefaction effects in computational models with the increase in the Knudsen values, a slip length model was introduced for the case of the Couette type of flows [20]. Since the latter model is applicable to only planar wall cases, a need arises to derive a velocity slip model that would, in the context of classical Navier-Stokes equations, accurately cover also the case of highly curved geometries, especially for the case of spherical and non-spherical particles.
There exist also numerical approaches that are not based on the classical form of the Navier-Stokes equations. In the context of the extended form of the Navier-Stokes equations, Guo et al. [21] derived a generalised second-order slip boundary condition that allowed to be applied for computation of slip regime flows over a sphere up to a Knudsen number of 0.6. On the other hand, using the extended Navier-Stokes equations form is not practical from a general CFD user point of view, applying vendor CFD codes with no such extensions available to the user. More complex numerical model were also derived, based on molecular dynamics simulations [22], and the solution of the Boltzmann equations that can be used in order to compute the fluid flow around a microparticle [23,24] even in cases of high Knudsen number values. However, these models are not directly applicable in the context of CFD. Note that also thermodynamic theory adds some restrictions on the coefficients in the first-and second-order velocity slip and temperature jump boundary conditions, see Sharipov [25].
The paper is organised as follows. First, the drag force models and the slip correction models for the Lagrangian point particle model and the case of a microsphere are described. This is followed by the description of the CFD-based numerical solution of a 3D flow past a spherical microparticle and its verification for the no-slip microsphere case. The main part of the paper is devoted to the study of the slip boundary conditions in the context of Navier-Stokes solution of flow past a microparticle. Here, the computational procedure for the derivation of a specific slip length (β) model, further denoted as polynomial model, is proposed in the form of the algorithm presented in Fig. 1.
In Sect. 5, the analysis of the obtained results using the Maxwell's and generalised Maxwell's slip models is presented. The paper closes with conclusions.

Drag force on a spherical microparticle
In the point particle approach, the motion of the dispersed particles is described in the Lagrangian frame, and a set of ordinary differential equations, which are given by the particle kinematics and Newton's second law, is evaluated along the particle trajectory to obtain the particle location, velocity, angular velocity and orientation. In the case of small particles, typically in the range of 0.1 µm ≤ d p ≤ 5 µm, the major forces are the drag F D , the lift F L , the buoyancy F B and the gravitational force F G . For a spherical particle in a no-shear flow, the lift is zero, and the gravitational and buoyancy force depend only on the volume of the particle and the fluid and particle density, whereas the drag force depends primarily on the fluid flow conditions and particle's shape. In the case of spherical particles, the drag force F D is given by [3]: where m p , V p , ρ p and d p are the mass, volume, density and diameter of the spherical particle. Moreover, u p is the particle velocity and u is the fluid velocity and ρ f denotes the fluid density. Furthermore, C d and C c denote the drag coefficient and the Cunningham correction factor [8]. In the case of Stokes flow, the drag coefficient has the form [26]: with denoting the particle Reynolds number, leading to the equation for the Stokes drag where μ denotes the dynamic viscosity. In Eq. (4), the correction to the drag coefficient, the slip correction factor C c , is applied, which is a standard drag model correction in the case of the slip flow regime. Data on the slip correction factor are available for the case of spherical particles in the form of slip correction models, proposed by Cunningham [8]: where Kn and A denote the Knudsen number [27], and the slip correction parameter. The Knudsen number compares the molecular mean free path λ to a representative physical length scale: As used by various authors, we employ the particle diameter d p as the characteristic length scale, see for example Tao et al. [24]. The mean free path λ is obtained from the relationship of the gas kinetic theory, see Jung et al. [15], as follows: using the mean velocity of the gas moleculesc, see Bird et al. [28], the density of ideal gas ρ = mp kT (9) and a constant φ, which depends on the kinetic theory model. In agreement with Jung et al. [15], who studied atmospherical air conditions, we employ the formulation of Chapman and Enskog [29] where e = 0 when there are repulsive forces between the molecules (φ = 0.491). In the above equations, k denotes the Boltzmann's constant, T the temperature, p the system's pressure and m the molecular mass. The mean free path for standard atmospheric conditions renders λ = 66.8 nm. For common aerosol sizes of d p = 1−10 µm, see Wedel et al. [2,30], we obtain Kn = 0.01−0.067 (slip flow regime). The slip correction parameter A in Eq. (5) is a function of the Knudsen number and three empirical constants denoted in this paper as a, b and c: Moreover, the empirical constants depend on the gas type and particle material [8,18,27]. These slip correction factor expressions are not directly applicable to particles of other shapes, as they were all developed based on experimental studies, involving spherical particles only. In order to develop drag force models for general non-spherical microparticles without the need to perform extensive experimental studies, high-fidelity CFD computational models for numerical simulations of flow around particles of different shapes can be applied. Such studies are frequently used in the development of point particle force models for no-slip conditions [7,17,31,32]. In order to apply a similar computational framework to derive force models for the slip flow regime, specification of physically realistic slip boundary conditions is the key to the success and hence of utmost importance.

Numerical solution of 3D flow past microparticle
To resolve the fluid flow around a spherical particle and determine the drag force, the Navier-Stokes equations are solved as the governing system of equations, i.e.
The viscous stress tensor τ is related to the velocity field by: In Eqs. (11)(12)(13), u, p and ρ denote the fluid velocity components, the pressure and the fluid density. The equations are solved with the simpleFoam solver of OpenFOAM ® [33, 34], which uses the finite volume method (FVM) to discretise the governing equations.
As the main focus of the present investigation is on the accurate evaluation of the drag force on microparticles, in order to solve the Navier-Stokes equations in the slip flow regime, slip boundary conditions have to be applied.

Slip boundary condition
The implementation of a slip boundary condition is a key challenge to obtain a computational procedure that enables direct evaluation of the drag force acting on a general non-spherical particle. The slip at the boundary depends on the conditions in the fluid phase, which can be assessed based on the Knudsen number Kn [27], see Eq. (6). When the particle diameter d p is in the order of the gas mean free path λ, the fluid molecules at the particle boundary experience on average a certain amount of slip along the wall. The original idea proposed by Maxwell [35] was to relate the slip velocity to the viscous shear stress at the wall, multiplied by a model parameter α, i.e.
Several forms of the slip velocity model have been derived by implementing different approximations to the shear stress, mainly by using first-order velocity derivatives. The most used form of the Maxwell model builds on the variation of the tangential velocity in the direction normal to the wall, denoted here as the conventional Maxwell model. When dealing with curved surfaces, as is the case for spherical and for the vast majority of non-spherical particles, the variation in the wall normal velocity should not be neglected. To take this into account, a generalised version of the Maxwell's slip velocity condition considers also the normal velocity derivative in the tangential direction [19]. In Table 1, several velocity slip models in the conventional and generalised Maxwell forms are referenced. Further, it has to be noted that higher-order slip velocity models have also been developed [21], although their implementation in the context of CFD is more computationally demanding and is therefore not considered in this work. In order to model the slip conditions within the CFD framework, we propose the following generalised slip velocity model: for which we introduce a specific slip length model β = β(Kn), derived from published experimental studies on microparticle drag force correlations and Cunningham-based slip correction factors at standard conditions employing the generalised Maxwell model formulation. By neglecting the term du n /dt in Eq. (15), one can Table 1 Specific slip length (β) models found in the literature obtain the simplified conventional Maxwell model expression. In Eq. (15) (n, τ ) denote the wall normal and tangential direction, respectively, and β describes the ratio of slip length S v to a characteristic length L. For the case of a spherical particle L = d p , we obtain where the slip length S v is defined as [36], employing the viscous slip coefficient σ p . Table 1 gives an overview of the slip length models found in the open literature, expressed in the form of the specific slip length β.
With adopting a fixed tangential momentum accommodation coefficient σ for a certain gas-solid surface, as commonly suggested [25,36,37], the slip length and likewise the slip velocity become linearly proportional to the Knudsen number, see Table 1. However, this assumption, as will be shown, leads to significant errors in computation of the drag force on a spherical microparticle in the slip flow regime, see Fig. 8 for a single σ p value. Moreover, Table 1 displays that the conventional Maxwell model is commonly employed. Nevertheless, when surface curvature is present, e.g. in the case of spherical and most non-spherical particles, the particle boundary normal velocity has an influence on the overall slip behaviour [19], which poses the need for a generalised Maxwell's slip velocity condition. Consequently the challenge of identifying an appropriate nonlinear specific slip length model basing on the generalised Maxwell model remains and is being targeted in the framework of this paper.

Geometric model and computational domain
Discretisation of the governing equations and numerical computation is performed with the open-source code OpenFOAM ® , which uses a FVM framework. The OpenFOAM ® solver simpleFoam and the S I M P L EC algorithm are adopted to solve the Navier-Stokes equations. A spherical mesh is generated by employing block Mesh and snappy H ex Mesh. This generation of the computational domain is an important part of the set-up, as the problem under investigation is of elliptical nature and care has to be taken in order to exclude the influence of the outer domain boundaries on the flow field around the particle. A stationary sphere (u p = 0) with diameter d p is placed at the centre of the spherical mesh to ensure equal distance to the outer boundaries. In the first step, the no-slip boundary condition is selected at the particle surface and the Stokes flow regime is computed to validate the numerical set-up. The computational domain and boundary conditions are sketched in Fig. 2.
Once the CFD solution is obtained, the fluid force exerted on the surface of the sphere can directly be obtained by integrating the pressure and viscous stresses over the particle's surface,   Fig. 2 Due to the symmetry of the sphere and the Stokes flow conditions, the obtained force has only one nonzero component, the drag force F D in the direction of the fluid flow. Using Eq. (1), one can then readily calculate the drag coefficient.
In the following, first the mesh resolution is investigated. Therefore, seven different meshes are generated by gradually refining the inner base cube, see Fig. 2. Simultaneously, the cell growth in the remaining spherical domain is adjusted to ensure a gradual growth towards the outer boundaries. The obtained mesh statistics are provided in Table 2. Figure 3 displays the resulting drag coefficient values, which are normalised to the Stokes drag coefficient. The classic Stokes coefficient (denoted as C d,St ) is given in Eq. (2) and is valid for Re p << 1 [3]. Several extensions to higher Re p values can be found in the literature [3]. For completeness, we include two extended drag laws [38]: which is valid up to Re p ≤ 5 and [39]: which is applicable up to Re p ≤ 800 [3]. In the following, we employ these three analytical solutions and empirical equations as references for verification purposes. For completeness, these extended drag laws are displayed in Fig. 3 alongside the present drag coefficients computed with the OpenFOAM ® simulations. Figure 3 indicates that from mesh M4 onward (see Table 2) mesh resolution independence is achieved. In general, a reasonable agreement with all drag law values is found for all the presented meshes. In the case of Re p = 0.1, increasing the number of base cube elements leads to drag coefficients that approach the numerical prediction of C d,SN . Besides, for Re p = 0.001 only a slight over-prediction is obtained as mesh convergence is reached. However, this over-prediction phenomenon is attributed to the domain size influences, rather than the mesh resolution, which is investigated in the following chapter. Consequently, we conjecture insufficient mesh resolution with drag force under-prediction. Based on the presented results, mesh M4 is considered as appropriate for further investigations.   [39], C d,O S [38], C d,St [26] Moreover, special attention has to be given to properly size the computational domain as a considerable challenge exists when employing computational methods for the prediction of extremely low-Re flows. In such cases, viscous stresses are predominant over convective stresses. Besides, the momentum transport is mainly governed by diffusion that affects a substantial spatial range when approaching the Stokes flow limit. Therefore, the computational domain span is considered a key factor in the process of establishing a grid independent solution. Multiple reports already exist where the significance of the domain and grid configuration is highlighted [31,40].
In the following, a domain sizes study for six different grid sizes (D1−D6) is conducted and presented in Fig. 4.
In Fig. 4, domain size independence is observable for both considered Reynolds numbers with increasing grid span width. For Re p = 0.1, the impact of the numerical framework is negligible from mesh D4 onward, where the drag coefficient approaches the reference results of Oseen [38], who extended the Stokes analysis by considering first-order inertia terms [3]. Besides, for Re p = 0.001, domain independence is evident for the mesh D6 where the present drag coefficients as well as both considered extended drag coefficient models approach the Stokes solution. Accordingly, we conjecture an insufficient domain size with drag force overprediction. Therefore, special attention has to be given to the combination of mesh resolution and domain size as both errors could cancel out in unfavourable situations. By using the computational domain D6, the influence of the numerical framework is considered negligible as mesh and domain independence is already established. Accordingly, the configuration of D6 with n cells ≈ 1.4 · 10 6 is chosen for further investigations.
In the following, cell count reduction is targeted to minimise the computational time. Hence, we reduce the amount of cells in the spherical far field and simultaneously increase the related cell aspect ratio. The resolution of the inner base cube and therefore the mesh in the proximity of the spherical particle remain unchanged in comparison with the mesh D6. With this procedure a cell count reduction of 73% is achieved (n cells = 0.4 · 10 6 ) while maintaining equal accuracy.

Verification of the numerical model: the continuum flow over a sphere
This section targets the no-slip flow over an unbounded sphere, which has been extensively studied by many researchers [24]. The main attention is drawn to the drag coefficient C d of the sphere. As displayed in Fig. 5, good agreement can be observed between the present and literature results for Re p < 0.1. For Re p > 0.1, the best agreement is achieved to Schiller and Naumann [39], which indicates the applicability of the set-up beyond the Stokes regime. Figure 6a-d presents the final overall spherical domain, the inner cubic element, the boundary layers and the resolution of the spherical particle.

Discussion on conventional and generalised Maxwell model
In this section, we discuss the use of conventional and generalised Maxwell models. The analytical total drag coefficient of an unconfined sphere in the slip flow regime is given by [19] and expressed in terms of the correction to the Stokes formula for the stationary sphere: with • Conventional Maxwell model: C c,conv = 1 + 2σ p λ/d / 1 + 4σ p λ/d , • Generalised Maxwell model: C c,gen = 1 + 4σ p λ/d / 1 + 6σ p λ/d .
In Fig. 7, comparison of the analytical descriptions of [19] to the present Stokes drag corrections, obtained by implementing β B Maxwell slip conditions, see Table 2, in the CFD framework, is presented. One can note an excellent agreement of the numerical and analytical values for both types of Maxwell conditions which validates the implementation of the slip velocity boundary condition in OpenFOAM as described in Sect. 9.   [15] However, Fig. 7 also clearly indicates that the consequence of misapplying the conventional slip boundary condition on curved surfaces results in a notable altering of the drag coefficient. The conventional form of the slip boundary condition leads to a significant underestimation of the drag compared to the generalised version, which is exacerbating with increasing rarefaction effects. This underestimation of the drag is even more pronounced for increasing σ p values as displayed in Fig. 7d. This states the importance of imposing the generalised boundary condition on spherical particles in order to properly capture the physics of the present application. On the other hand, the question arises whether the linear dependence of β(Kn) in the slip velocity boundary conditions can, in the framework of the CFD, lead to accurate results of the drag force computation in the case of objects with curved surfaces. As will be seen, this is not the case; however, a higher-order dependence of β on Kn does lead to significantly improved accuracy of the CFD computations.

Computational procedure to determine the β = β(Kn) relationship
As already stated, the main aim of the present investigation is the derivation of a specific slip length relationship that would enable a scalable set-up of the slip boundary conditions on microparticles according to their size. To obtain such relationship, empirical data on drag force measurements on microparticles in the form of oil drops and polymeric hard spheres (PSL) are used. In all the studies, a general relationship for the slip correction factor employs a slip correction parameter A, which is a function of the Knudsen number, see Eq. (5). The experimentally obtained expressions for the slip correction parameter A for different particle material, but of standard spherical shape, at standard conditions for air, are listed in Table 3.
Based on the known data from Table 3, the computational algorithm for the determination of the β = β(Kn) relationship is set as follows:

Algorithm 1: Determination of the β = β(Kn)
1. Choose the model case No. X from Table 3  With the obtained data on the β = β (Kn) relationship, the derivation of an approximation function to fit the data of the step 3 in Algorithm 1 is performed with the following approximation choice: which is further denoted as the polynomial model. Note that b 0 = 0 for all the considered models.

Results and analysis
In this section, the computational method and results are presented. First, to describe the computational approach to determine the β(Kn) relationship, one of the available experimental slip correction studies is selected and the results of performing the distinct steps in Algorithm 1 are reported. In the second part of this section, we report on the attempt to derive a more general β(Kn) model by incorporating the data from the considered experimental studies into the approximation procedure, and on the analysis of the resulting errors in the quality of the approximations. Finally, the results obtained for the derived models are compared to existing data and guidelines for their use are given.

The single slip correction case
The methodology for fitting β(Kn) to experimental results is exemplarily described for case No. 2 [16]. First, the Maxwell model proposed by Schaaf and Chambre [9] is implemented in OpenFOAM ® . In the following, several CFD computations are performed by varying the Knudsen number and σ p = β/Kn values. Note that a variation of σ p for a given Knudsen number leads to a change of β. As stated by Sharipov [25], basing on the CL scattering law the upper limit for is σ p = 2.845. In the case of diffuse scattering, σ p values can reach higher values; however, Sharipov [25] recommends employing σ p = 1.0. To obtain an appropriate fit for the mentioned range, we, therefore, choose to investigate a σ p,max = 3.0 to enclose the mentioned range of Sharipov [25]. In this regard, a σ p range of 1.0 ≤ σ p ≤ 3.0 for different cases of the Kn values is examined. Note that the conventional Maxwell model was derived for planar and non-rotating surfaces, and thus is not intended for use on curved surfaces [19]. Although in the present case we are dealing with the surface of a sphere, it is nevertheless interesting to investigate the implications of using the conventional (inappropriate) form of the Maxwell slip boundary condition on the spherical particle, as it is often misapplied [19]. Figure 8 shows the computational results for C d obtained from the CFD model results normalised by the Stokes model solution C d,St , using the conventional and generalised Maxwell implementation. In addition, Fig. 8 includes the reference results based on the experimentally derived equation of No. 2 [16]. It is evident that the application of a fixed σ p value, which leads to a linear β(Kn) relation, cannot reproduce the experimental results of Rader [16]. The deviations of this linear β(Kn) model to the experimental results are especially  Moreover, Fig. 8 is highlighting the upper Kn limit for which the fitting procedure of general and conventional model are applied for the reference case No. 2 [16], i.e Kn ≤ 0.22 for conventional and Kn ≤ 0.12 for generalised Maxwell model. With this limitation we ensure a σ p,max ≤ 3.0 to represent the whole slip coefficient range mentioned by Sharipov [25].
In the next step, the Maxwell model, see Eq. (15), is adapted to the reference case No. 2 [16], by deriving an appropriate β = β(Kn) relation. In this regard, appropriate β estimations for varying Knudsen numbers, which are further denoted as β * , need to be obtained. The β * values are estimated by employing the following equation: In the following, the present specific slip length model, see Eq. (23), is fitted to the identified β * values for case No. 2 (β * N o.2 ), which is achieved by using the MATLAB ® toolbox for curve fitting. The fitting parameters obtained are listed in Table 4. Also, the approximation functions are visualised in Fig. 9 along with the obtained β * N o.2 values. As indicated in Fig. 10, the polynomial fit achieves excellent agreement with the reference results until the upper fitting limit of Kn ≤ 0.12. The same methodology can be used to generate fits for different literature references. The obtained fitting parameters for the presented approximation models for six selected experiments are presented in Table 4. Furthermore, the resulting approximation functions obtained by the generalised Maxwell model are displayed in Fig. 11.

The β analysis
Finally, after the successful derivation of the β = β(Kn) relationships for individual reference data, an attempt is made to derive a more general β model by simultaneously incorporating the data from all experimental   [20] studies under investigation into the approximation procedure. In summary, data from experimental results for oil droplets [11,13,16], and for polymeric hard spheres [11,14,15], are considered here with the aim of deriving three different general relations for β: oil droplets in air (Model 1), polymeric hard spheres in air (Model 2) and ensemble data of all six reference cases (Model 3). Figure 12a, b displays the resulting normalised drag coefficients of the investigated experiments for oil droplets in air as well as PSL spheres in air and indicates that only small deviations occur between the considered experimental curves over the investigated Knudsen range, which provides a solid basis for the derivation of a more generally valid β model.
In the following, we present the results of the derivations of approximation models for β for the generalised Maxwell model. First, β values are collected from the analysis described in Sect. 5.1. Then, the polynomial model approximation functions are derived for the oil droplets in air, polymeric hard spheres (PSL) in air and for the ensemble data. The resulting models, together with the corresponding experimental values, are presented in Fig. 13a-c, while a graphical comparison of the generated approximation functions is given in Fig. 13d.  Table 3; [11], [13], [16], [12], [14], [15] (a) (b) (c) (d)

Fig. 13
Generalised Maxwell model: Derived group models for oil droplets, polymeric hard spheres (PSL) in air and ensemble model. Experimental data for oil droplets in air: [11], [13], [16] and for PSL spheres in air: [15], [14], [12]. Derived models: Model 1, Model 2, Model 3. Reference linear slip length model by Pan and Liu [20] [11,13,16] b [12,14,15] As visible in Fig. 13a-c, the obtained model functions are a sufficiently accurate representation of each experimental data set over the Knudsen range studied. Moreover, Fig. 13d indicates a good agreement between the derived models for Kn < 0.1. However, towards higher Knudsen values, the discrepancies between the approximation functions become more pronounced, leading to comparatively higher β values for Model 1 (oil droplets in air) and lower β values for Model 2 (PSL spheres in air) ( Table 5).  [16]. Figure 14 clearly presents that the considered slip models found in the literature that base on a linear β(Kn) function are unable to capture the experimental results of Rader [16] for conventional ( [19] (conv.), [20,25]) as well as generalised Maxwell formulation ( [19] (gen.)). This clearly highlights the need for a nonlinear β(Kn) relation in order to overcome the strongly rising deviations of linear models for Kn > 0.025. Solely in the region Kn < 0.025, the discrepencies between the considered literature models and the present polynomial model can be neglected, as underlined in Fig. 14a. For Kn > 0.025, the advantages of a nonlinear β(Kn) relation, as employed in the presented polynomial model, outweigh the standard linear models. This is indicated by an excellent agreement of the present polynomial model (Model 3) to Rader [16] up to Kn ≈ 0.15 with a maximum deviation of 1.7%, whereas the considered linear models already reached a deviation of more than 9.1%.

Conclusions
As the size range of particles under consideration in computational multiphase studies reaches the micrometre and sub-micrometre range, it is important to implement particle fluid dynamic force models able to accurately cover also the slip flow regime. With the point particle approach and the case of a sphere the Cunningham slip correction models are commonly used. When non-spherical particles are encountered, slip correction models have to be modified in order to account for the particle's shape. To achieve this, a particle resolved fluid flow simulation by means of CFD, presented in the paper, needs to implement slip flow boundary conditions at the particle surface. The most common choice is the Maxwell slip velocity condition, with the tangential momentum accommodation coefficient as the main model parameter. Since the latter is recognised as a constant value that does not depend on the Knudsen number, the Maxwell type boundary condition becomes a linear function of the Knudsen number. When used in the framework of a CFD simulation of a flow past a sphere, significant errors in the computed drag force values in the slip flow regime can therefore occur, as highlighted in the paper. In order to avoid this effect, the Maxwell slip conditions are rewritten with the specific slip length as the main modelling parameter. Since particle shapes are seldom composed of planar surfaces, for which the conventional Maxwell slip model is valid, the generalised version of the Maxwell's slip velocity condition is applied. The available results of the experimental drag studies on spherical microparticles in the form of oils drops and polymeric hard spheres are used in a novel dedicated CFD-based computational algorithm to obtain the correct specific slip length values for the Knudsen range of 0 < Kn < 0.12. This is followed by derivation of a polynomial function-based specific slip length dependency on the Knudsen number, which can be used in a CFD study of microparticle dynamics in the slip flow as well as in the no-slip flow regimes. Furthermore, it is shown that implementing the conventional form of the Maxwell slip boundary condition leads to a significant underestimation of the drag compared to the generalised Maxwell form, which is exacerbating with increasing rarefaction effects. Finally, we conclude that the present investigation introduced a computational procedure for the determination of the specific slip length model based on experimental slip correction data of spherical microparticles and estimate that the developed model is applicable also to general convex non-spherical particle shapes, delivering superior results with respect to imposing classical form of the Maxwell slip conditions, derived for the planar geometries.