Modeling of macroscopic stresses in a dilute suspension of small weakly inertial particles

In this paper we derive asymptotically the macroscopic bulk stress of a suspension of small inertial particles in an incompressible Newtonian fluid. We apply the general asymptotic framework to the special case of ellipsoidal particles and show the resulting modification due to inertia on the well-known particle-stresses based on the theory by Batchelor and Jeffery.


Introduction
Suspensions of small, arbitrarily shaped particles in a Newtonian fluid are of great interest in physical, biological and engineering sciences, see among others [22,5,16].Here one of the important questions, which often arises, is how the carrier flow is influenced by the suspended solids.This effect strongly depends on the particle size and mass, the characteristics of the flow but also on the scale of interest.In general, the initial configuration of the particles in the fluid domain involves some kind of stochastic nature, additionally the particles may interact with each other and also turbulent flow fluctuations may strongly alter the deterministic path of each particle.Thus looking on a cut-out of the domain on the microscale, where each particle has a significant size and that way a strong impact on the streamlines, the problem is dominated by stochastic effects.On the macroscale on the other hand, where the size of the particles is small compared to the dimensions of the flow, the suspension can be described as a homogeneous, non-Newtonian fluid and the task consists of modeling the corresponding bulk stresses.
One step towards the derivation of the macroscopic suspension behavior is the study of the two-way coupled particle-fluid problem.Oberbeck described the disturbance flow generated by an ellipsoidal particle moving with a constant translational velocity through a stationary viscous fluid in [17], Edwardes extended these results to the case of a rotating ellipsoid in [6].In the case of a dilute suspension of small rigid spheres in a stationary Newtonian fluid, Einstein gave a derivation of the change in viscosity of the carrier fluid by first considering the motion of a single particle in a sufficient small region, such that the corresponding flow may be treated as linear, and then solving the stationary Stokes equations for the disturbance flow.The contribution of many particles followed by superposition of every single disturbance flow field [7].An analogous procedure was applied to ellipsoidal inertia-free particles by Jeffery [10].The model equation that describes the evolution of the principal axis of a prolate ellipsoid of revolution in a viscous fluid is known as Jeffery's equation.Batchelor presented a general framework of modeling stresses in a suspension of rigid particles in an incompressible Navier-Stokes flow [1].In the case of a dilute suspension, he applied his theory to ellipsoidal inertia-free particles using the results of Jeffery and derived an analytical term for the particle-stresses in the macroscopic description of the suspension.The study of non-dilute suspensions involves the modeling of the interactions between particles, see for example [2,14,4,8].Another general approach for suspension modeling was presented in [25], Date: January 29, 2022 ⋆ Corresponding author, email: vibe@math.fau.de,phone: +49 9131 85 67215 1 FAU Erlangen-Nürnberg, Lehrstuhl Angewandte Mathematik I, Cauerstr.11, D-91058 Erlangen, Germany .
1 where Cauchy's stress principle was applied to the suspension, reproducing classical results as well as showing additional effects induced by spatial non-uniformities of the dispersed phase.The influence of particle concentration was also considered in [20,21].Rigorous asymptotical homogenization strategies for suspensions were used for example in [12,3].The work by Junk & Illner [11] deals with the strict asymptotic derivation of Jeffery's equation, using expansions in the small size ratio (ratio between the characteristic length scale associated to the particle and the one associated to the fluid).It yields a general strategy for constructing a correction to the undisturbed Navier-Stokes solution to account for the presence of a particle.
In this paper we extend the asymptotic approach by Junk & Illner by taking into account inertia.We particularly classify three different types of inertia and use the asymptotic results for a single particle in the general framework of suspension modeling according to Batchelor.We study the case of a dilute suspension and consider particle-fluid interactions as well as gravitational forces.The regard of inertia requires the introduction of a special splitting for the forces associated with the bulk stresses.This way we deduce an extended influence of the solid phase on the macroscopic behavior of the suspension induced by small deviations of the particle movement from the streamlines of the carrier fluid (Theorem 11).Apart from a modification of the particle-stress tensor we obtain an additional effective particle force term.We demonstrate the inertia related effects for the wellstudied example of a dilute suspension of (prolate) ellipsoidal particles (Corollary 14).
We organize this paper as follows: In Section 2 we set up the model for a single particle in a fluid flow with respect to different types of inertia.We give a brief overview of the asymptotic derivation including the essential definitions and assumptions.In Section 3 we sketch the basic ideas of Batchelor's framework and use the asymptotic results to formulate corresponding consistent suspension models which account for weakly inertial particles.We illustrate the results for the special case of ellipsoidal particles in Section 4 and conclude with a summary in Section 5. Technical details to the example are provided in the Appendix.
Throughout this paper we use the following notation: Notation 1. Scalar values are typeset as ordinary characters a ∈ R. Vectors are indicated by small bold characters, v ∈ R n , n > 1, their components are denoted by (v) i = v i ∈ R, i = 1, . . ., n.In the following we will omit the ranges of the indices if they are obvious.Matrices are written as large bold letters, M ∈ R n×m with the components (M) ij = M ij , the identity matrix is denoted by I.
Tensors of higher order are treated as linear mappings between the corresponding spaces and written as ordinary large characters, e.g.B : R 3 → R 3×3 , or as the corresponding dyadic product.We use a tensor calculus notation, viz.
analogously for higher order tensors: and so on.For v, w ∈ R 3 the cross-product is given by (v×w) k = i,j ǫ ijk v i w j , with the Levi-Civita symbol for (i, j, k) an even permutation of (1, 2, 3), −1, for (i, j, k) an odd permutation of (1, 2, 3), 0, else.
We use the mapping B(v) ij = k ǫ ikj v k which assigns a vector to a skew-symmetric matrix such that B(v) As for the differential operators, the Jacobian is denoted by (∂ y f) ij = ∂ yj f i and the gradient by its transpose (∇ y f) ij = ∂ yi f j .For scalar-valued functions these notations are certainly interchangeable.As usual, the nabla operator is formally treated as a vector and for example used for the divergence and rotation operators: div(u) = ∇ • u, rot(u) = ∇ × u.Higher order derivatives are noted as ∂ y...y u for a sufficiently smooth u.

Mathematical model of a single particle in a fluid
In this section we provide the ingredients for the suspension model.For this purpose, we consider a single small oriented particle suspended in a Newtonian fluid and study its influence on the carrier flow under the assumption that the particle is small compared with the dimensions of the flow.Proceeding from the three-dimensional interface problem we derive an asymptotic model for the particle in the fluid flow by help of expansions in the size parameter ǫ (size ratio between particle and flow domain).The asymptotic analysis follows the procedure that was established for inertiafree particles in [11].The novelty is the extension to inertial particles which we classify with respect to different types.

Three-dimensional interface problem.
Dimensional problem.A particle is a rigid body that we model as an open, bounded and timedependent domain E(t) ⊂ R 3 with a smooth boundary ∂E(t) and time t ∈ R + 0 .The particle is characterized by its center of mass c : R + xdx and an orthonormal righthanded director triad d i : R + 0 → S 2 1 , i = 1, 2, 3, where S 2 1 denotes the unit sphere in R 3 .Additionally we define a time-invariant reference state E ⊂ R 3 with the center of mass lying in the origin and the directors identified with the Cartesian basis vectors {e 1 , e 2 , e 3 }.The function which maps every point of the reference state to the actual time-dependent state is completely presented by rigid body motions (translation and rotation), i.e. x : 2.1), with R ∈ SO(3) being the rotation matrix associated with the director triad, R ij = e i • d j .We assume that the small oriented particle is suspended in an incompressible Newtonian fluid.The regular domain Ω contains the fluid and the particle.The acting forces are due to particle-fluid interaction and gravity.The model of first principles for the particle and the flow consists for all t > 0 of the incompressible Navier-Stokes equations in the fluid domain Ω \ E(t), the no-slip condition on the interface as well as the kinematics and dynamics of the particle: Here, u and p denote the fluid velocity and pressure, S[u] = −pI + µ f ∇u + ∇u T the Newtonian stress tensor with µ f being the fluid viscosity and ρ f the fluid density.The particle related quantities v and ω describe the linear and angular velocities of the rigid body.The inertia tensor is given by T with the time-invariant part Ĵ = |E| −1 E y 2 I − y ⊗ y dy and ρ p is the particle density.The gravitational acceleration is given by g with normalized direction e g .The system is completed by appropriate initial conditions for the flow and the particle as well as boundary conditions for the flow on ∂Ω.
Dimensionless formulation.We introduce two characteristic lengths associated with the fluid system x and with the particle y as well as a characteristic time t.We scale the velocities as u = v = x/t and the pressure with p = uµ f /x.For the sake of a simple notation, we keep the same symbols for all variables as in the dimensional formulation and get the dimensionless system: ) ) ) The model (2.2) is characterized by four dimensionless parameters: the size ratio ǫ = y/x, the density ratio ρ = ρ p /ρ f , the Reynolds number Re = u xρ f /µ f (ratio of inertial and viscous forces in the fluid) and the Froude number Fr = u/ √ xg (ratio of inertial and gravitational forces in the fluid).The Newtonian stress in dimensionless form * reads as S[u] = −pI+2E[u] with E[u] = 0.5(∇u+∇u T ).The size parameter enters the bijection between the reference and time-dependent state: x(y, t, ǫ) = ǫR(t) • y + c(t), (2.3a) The density ratio ρ is now the key parameter that allows us to define different inertial regimes.

2.2.
Tracer particles of different inertial type.As we will show our asymptotic approach results in particles which follow the streamlines of the surrounding fluid.Hence, we refer to them as tracer particles.To distinguish between different types of inertial particles, we introduce an ǫ-dependent mass function α mass (ǫ).With the particle density ρ p = ρ p (ǫ) = α mass (ǫ)m p /V p , the density ratio and the size parameter are then related according to Here m p , V p and m f , V f are the mass and the volume associated with the particle and the fluid, respectively, m = m p /m f .The freely selectable mass function has no physical meaning, but allows to balance the inertial terms in the particle momentum balance in different ways, yielding several models of inertial particles.In particular, we set α mass (ǫ) = ǫ k , k ∈ N, and define three inertial regimes with respect to the behavior of the density ratio in the limit ǫ → 0, i.e. heavy (k = 2), normal (k = 3) and light weighted (k > 3) tracer particles, see ratio satisfies ρ ≡ const leading to the inertia-free particle model that was investigated in [11].The cases k ≤ 1 are not covered by the asymptotic simplification as we will comment on in Remark 7.
2.3.Asymptotic analysis.In this subsection we generalize the asymptotic approach of [11] to the introduced inertial type classification.We particularly point out the differences in the asymptotic analysis that arise for the heavy tracer particles.
Asymptotic expansion.We assume that under the considered inertia types the particle influences the carrier flow locally but not globally.This way we follow [11] and make the following assumptions.
1) The fluid is essentially undisturbed by the particle.
2) The fluid motion induces a rotation of the particle of O(1).
3) The particle is far away from the boundary of the fluid domain ∂Ω.
For the quantities of the particle we take a regular expansion in powers of ǫ, while the fields of the fluid flow are separated in a global and a local part: The local fields of the fluid quantities are expressed in the local coordinates of the particle reference state, The scaling is chosen in such a way that the gradient of the pressure balances the Laplacian of the velocity.It is a consequence of the fact that the bijective mapping between the reference and time-dependent state is ǫ-dependent, see (2.3b).(This way the chain rule generates a factor of ǫ −1 every time u loc is being differentiated with respect to x or t.) Asymptotic solution.Our goal is to formulate a solution of the complete system (2.2) up to an error of O(ǫ).Therefore, we consider only a finite number of asymptotic coefficients in the expansions.
To address them we set ) Inserting the expansion coefficients into (2.3)yields the following form of the bijective mapping between the reference and the time-dependent state: x e (y, t, ǫ) = ǫR e (t, ǫ) • y + c e (t, ǫ), It is shown in ([11], Lemma 5) that (R e ) −1 = (R e ) T + O(ǫ 2 ) and (R e ) −1 = O(1), thus y e is well-defined.The asymptotic coefficients of the mapping are particularly given by Moreover, we use Taylor expansions for functions of the form f (x(y, t, ǫ)) in c 0 in the following, i.e.
with the differential operators To emphasize the general structure of the asymptotic result in the subsequent Lemma 4 and to stress the relevance of certain terms for the suspension model in Section 3 we introduce abbreviations for some expressions. , where the derivatives of u 0 and the Newtonian stresses are evaluated in c 0 .Moreover, abbreviate the linear accelerations by and the angular accelerations by Lemma 4 (Asymptotic one-particle model).Let the following four requirements be fulfilled, then u e , p e , c e , v e , ω e and R e defined in (2.5) are a solution of the complete system (2.2) up to an order of O(ǫ) for ǫ ↓ 0. (R1) Let the global flow velocity u 0 and pressure p 0 be the solutions of the incompressible Navier-Stokes equations in with the boundary condition u 0 = u on ∂Ω, and the initial condition u 0 (x, 0) = u(x, 0) for x ∈ Ω \ E(0), and (R2) Let the particle related coefficients satisfy the following set of conditions: Let with the initial conditions c 0 (0) = c(0), c i (0) = 0, i = 1, 2 and let the zero-order velocity suffice the so-called 'tracer condition' Let additionally the matrices R 0 , R 1 solve the differential equations with the initial conditions R 0 (0) = R(0) and R 1 (0) = 0. (R3) For every t > 0 let the local fields solve the stationary Stokes equations on the unbounded exterior of E with the Dirichlet and integral conditions for i = 1, 2 and the functions h i , f i , g i given in Abbreviation 3. Additionally let u loc,i and p loc,i fulfill the following decay properties: for some c i , d i , f i , i = 1, 2 independent of y. (R4) Let the mass function (2.4) be given as α mass (ǫ) = ǫ k with k ≥ 2.
Proof.The proof results from a straight forward computation where the asymptotic coefficients are inserted into the full system (2.2) and Taylor expansions are applied for the global flow fields, whenever they are evaluated at the particle boundary (see [11] for normal tracer particles).The introduction of the different inertial types does not change the general mathematical structure of the problem but only affects the linear and angular momentum balances of the particle, i.e. the functions f i , g i in (R3).We sketch the steps of the proof for the sake of completeness.
Inserting the asymptotic expansions for c, v, ω and R in (2.2c) and applying (2.7a), (2.7c) shows that the particle kinematics (2.2c) is fulfilled up to an order of O(ǫ 2 ).
The conservation of mass in (2.2a) is fulfilled exactly since each local velocity field is assumed to be divergence-free.
Considering the no-slip equation (2.2b), the use of the asymptotic expansions and the Taylor series of u 0 in c 0 results in the following conditions that are exactly the Dirichlet conditions of u loc,i (2.8b) and the tracer condition (2.7b), Inserting the asymptotic expansion in the momentum equation of the fluid and applying the chain rule for the local fields yields Since u 0 satisfies the Navier-Stokes equations and the local fields satisfy the stationary Stokes equations according to (R1) and (R3), the momentum balance holds up to an error of O(ǫ).
As for the momentum balances of the particle, we insert the asymptotic expansion on both sides of the equations.On the right-hand side we use a Taylor series for the term S[u 0 ] in c 0 and keep the integrals for the local fields, this leads to for the linear momentum and for the angular momentum.These are exactly the integral conditions (2.8c) and (2.8d).
The requirement (R4) is needed to avoid contradictory conditions as we explain in Remark 7.
In Lemma 4 we assume that (R3) is fulfilled despite the fact that the problem seems to be overdetermined at first glance by the presence of Dirichlet as well as integral conditions on ∂E.This is however not true as the following lemma shows.
Lemma 5 (Solvability conditions of Stokes problem).Let w q , p q , q = 1, . . ., 6, be the solutions of the following six stationary Stokes problems related to the boundary conditions of the three elementary translations and rotations: w q = e q , w q+3 = y × e q , q = 1, 2, 3, y ∈ ∂E, w q ≤ c q / y , ∇ y w q , p q ≤ d q / y 2 , y ≥ f q > 0 (2.10) for constants c q , d q , f q .Set the following surface moments associated with w q , q = 1, . . ., 6 to be 1) Let (R3) of Lemma 4 be fulfilled, then the following solvability conditions hold: ) Proof.The key ingredient for this proof is the existence of the Green formula for solutions of the Stokes problem ( [11], Lemma 6), which connects the stresses and boundary conditions of two Stokes solutions and which needs the decay properties (2.8e) as well as those in (2.10).
For 1): The equations (2.8c) and (2.8d) are equivalent to for q = 1, 2, 3.With the help of the Green formula the roles of w q and u loc,i in the right-hand sides of (2.12) can be interchanged.Thus (2.12) is equivalent to We note that for i = 1, 2 the terms h i of Abbreviation 3 have the form , which holds for any M ∈ SO(3) and any ω ∈ R 3 , and where Inserting h i and sorting the resultant terms with respect to ω i−1 and v i leads to which directly results in (2.11).
For 2) we consider (2.11) which is equivalent to (2.14) and thus to (2.13), since h i is the Dirichlet condition of u loc,i by assumption.We use the Green formula and get (2.12).
For 3) we analogously arrive at (2.13).The assumption that the integral conditions are fulfilled is equivalent to (2.12).Combining both results gives the equalities Since u loc,i and w q fulfill the decay property the above equations imply u loc,i = h i on ∂E by means of the Green formula.
1) The solvability conditions (2.11) together with the particle related equations (2.7) build a system of nonlinear ordinary differential equations (ODEs) for c i and R i−1 , i = 1, 2, where the ODE for the ith asymptotic coefficient depends on lower order coefficients, while the surface moments depend only on the geometry of the particle.The order of the corresponding ODEs depends on the choice of the mass function α mass .For example in the case α mass (ǫ) = ǫ, f 1 depends on dv 1 /dt and thus (2.11a) becomes second order for c 1 with an additional initial condition v 1 (0) = 0.However, under suitable regularity assumptions on u 0 , the corresponding right-hand sides of the ODE systems are at least continuous and thus at least local solutions exist by Theorem of Peano.2) Solving the solvability conditions of the Stokes problems (2.11) determine the Dirichlet conditions (2.8b).However, it is not our aim to examine the existence, uniqueness and regularity of the Stokes solutions given by (2.8b), (2.8e) in general in this paper.Instead we will discuss some special problems in Section 4 where analytical solutions are available.3) Just like in [11], Lemma 4 and 5 imply a procedure for constructing a O(ǫ)-solution of (2.2): First the Navier-Stokes equations (2.6) are solved to get u 0 and p 0 , then the path c 0 of the particle is obtained from the tracer condition (2.7b).The solution of the solvability conditions (2.11a) together with the corresponding differential equations (2.7) for i = 1 provides the coefficients c 1 , v 1 and ω 0 , R 0 that define the Dirichlet condition (2.8b).The stationary Stokes problem (2.8) for i = 1 is solved.Finally, the last two steps are repeated to get u loc,2 which completes the procedure.
Remark 7 (Density ratio scaling).We want to point out the effect of the density ratio scaling introduced in (2.4).The momentum balances of the particle (2.9) can be written as i+1 , m ω i+1 denoting the according integral terms appearing on the right-hand sides of (2.9a) and (2.9b).The terms k i , ℓ i describe the linear and angular accelerations of the particle, while the integral terms connect them to the Dirichlet conditions of the corresponding local flow fields via the solvability conditions (2.11).Different choices of the mass function lead to the following pattern for the balancing of the accelerations with the integral terms: (2.15)Here each row shows the O(ǫ ℓ )-correction of the momentum equations and each column corresponds to the choice of power of ǫ of the mass function.
For k ≤ 1 the condition dv 0 /dt = Fr −2 e g arises (see (2.15)).In combination with the tracer condition (2.7b), which is independent of the mass function α mass , this leads to the implication that du 0 (c 0 (t), t)/dt = Fr −2 e g in contradiction to the claim, u 0 solving the Navier-Stokes equations in Ω, according to (R1) of Lemma 4. This fact indicates that in this regime Assumption 2, 1) does not hold anymore, since the inertial effects of the particle are too strong and the perturbation of the surrounding fluid is of the same order as u 0 .Thus, in consistence to the asymptotic, we restrict our classification of inertial types in Table 2.1 to k ≥ 2.
2.4.Extension of Stokes solutions to the particle domain.The formulation of the suspension model requires the definition of velocities and stresses on the whole domain Ω, including the interior of the particle.Hence, we need to introduce a meaningful extension of the quantities to the particle domain E. A possibility is to consider E as a part of the fluid domain and use the so-called singularity solutions which express the Stokes solutions by means of the Green's dyadic (Oseen-Burgers tensor) [13].Such solutions have the disadvantage of not being bounded in the neighborhood of the origin, which is an important requirement in our modeling as we will see in Section 3.Alternatively, one could also think of a continuous extension of the Dirichlet conditions (2.8b) to E and model the stresses as Newtonian.This leads to a discontinuity in the stresses at the particle boundary which contradicts with the desire for a smooth macroscopic stress on Ω.Therefore, we follow here the approach of [18] and treat the particle as a fluid with a rigidity constraint.Proceeding from a dimensional description (analogously to Section 2.1), we obtain a consistent model for the dimensionless stresses and velocities in the particle domain by help of the presented asymptotic techniques.
The fluid in the particle domain is characterized by the velocity z : E(t) × R + 0 → R 3 of a material point and the non-Newtonian stress T : for all t > 0 with the flow velocity u solving (2.1a) and appropriate initial conditions in consistency with (2.1).The stress T acts as Lagrange multiplier to the rigidity constraint in (2.17a) and can be expressed as symmetric gradient field of the three-dimensional unknown λ : . † The velocity in the rigid body domain can be explicitly stated in terms of the particle quantities (center of mass c, linear and angular velocities v, ω of (2.1)) as z = B(ω)•(x−c)+v.Using this expression and the ODE for the director triad (2.1c) simplifies (2.17) to a boundary value problem for T, resp.λ: To transform the problem to the particle reference state we set In dimensionless form it is given by We model T as composition of the Newtonian stress of u 0 and some disturbance, T(y, As asymptotic expansion we take a regular power series in the size parameter ǫ for T loc and use a Taylor series in c 0 for S[u 0 ](x(y, t, ǫ), t).With we can formulate the asymptotic model for the stresses in Lemma 8. † The following considerations can equivalently be formulated for λ, but since we are only interested in the stresses themselves, we do not use the inner structure of T explicitly.It is however important to note that T have only three degrees of freedom such that the subsequent boundary value problems are well-posed.
Lemma 8 (Asymptotic model of stresses in particle domain).Let the requirements of Lemma 4 be fulfilled.Let T loc,i = ∇ y λ loc,i + ∇ y λ T loc,i , i = 1, 2 be the solutions of the following boundary value problems for t > 0 where the force terms are given by Then Te , R e , c e solve the system (2.18) at least up to an error of O(ǫ).
Proof.Similarly to the proof of Lemma 4, the approximation order follows from a straight forward computation.Inserting Te on the left side of (2.18b), using the expansion described in the proof of Lemma 4 for S[u] on the right-hand side and taking (2.19b) into account yields that (2.18b) holds up to O(ǫ 2 ).
When taking the divergence of Te with respect to y, the terms of the form R i • S[u 0 ] • R j vanish, since the Newtonian stresses are evaluated at c 0 .With the asymptotic expansions of R and c, (2.18a) is fulfilled up to O(ǫ) provided that (2.19a) holds.

Remark 9 (Local velocity fields and forces in particle domain).
1) By construction, the stress T induces a rigid body velocity z = ǫ(dR/dt) • y + v in E, while its approximation T e = R e • Te • R eT results in the approximation of z, namely Since the functions v i and R i do not depend on y, we can rewrite z e as z e = u 0 (x(y, t, ǫ), t) + ǫR e (h 1 (y, t) + ǫh 2 (y, t)) + O(ǫ 3 ), using Abbreviation 3, where h i (y, t) is the continuous extension of the Dirichlet conditions of u loc,i from Lemma 4 to E. This way we can extend the local disturbance velocity fields to E by setting u loc,i = h loc,i .The so defined velocity is bounded by the assumptions of Lemma 4.

2)
The force terms fi are not equal to f i of Abbreviation 3, but E fi dy = f i holds true in consistency with (2.19b).

Kinetic model for a particle suspension
In this section we present our main result, the asymptotical (macroscale) description for a particle suspension under consideration of inertial effects.For this purpose we set up mass and momentum balances by considering the particle suspension as a stochastic homogeneous fluid whose behavior is characterized by non-Newtonian stresses.To obtain a model for the qualitative behavior of this fluid, which is independent of the specific realization, we analyze the averaged equations.With the help of the ergodicity assumption we are able to derive an analytical expression for the bulk stresses.This procedure goes originally back to Batchelor in [1], the novelty of our work is the use of rigorous asymptotical results as model ingredients as well as the regard of inertial effects.These effects involve additional net volume forces and arise from the small deviations of the particles motions from the streamlines of the surrounding fluid.

Macroscopic stresses and forces.
The core of the suspension model (3.2) are the macroscopic stresses and forces that we deduce from the asymptotic one-particle model.We realize the fluctuation related quantities that are marked with the index ′ by superposing the local disturbance fields .loc of the particles (Assumption 10, (A6)).
Surface stresses.The derivation of the model for Σ s goes along [1].It is based on the idea of applying the ergodicity hypothesis for E[Σ s ] and thus integrating the stresses over a suitable averaging volume V(x, t) which contains N (x, t) particles E k (t), k = 1, . . ., N .The underlying assumption on where the stress fluctuations outside the particle are treated as Newtonian and inside the particles with respect to an appropriate extension.We particularly model the velocity and stress fluctuations by means of the local disturbance fields given in Lemma 4 for u ′ and in Lemma 8 for T ′ .With (A2) we may restrict to the influence of the local fields generated by the kth particle on E k and neglect particle interactions.Using the two identities and the divergence theorem on (3.3) results in Hence, we obtain for the averaged surface stresses where S p denotes the particle-induced stress tensor.The arising pressure term is incorporated here in the Newtonian stresses (as Lagrange multiplier to the incompressibility constraint).To express the integrals in S p (3.4b) with respect to the particle reference state we use the identities and for any scalar-valued smooth function f .This implies Using Lemma 8, the terms associated with the center of mass of each particle c k cancel each other out and, since f1 is independent of y and E ydy = 0 by definition of center of mass in reference state, we get Body forces.For the random body force we use a similar approach.We assume that the overall fluctuating force in the averaging volume V(x, t) is generated by contributions of each single particle in the corresponding particle domain where I denotes the indicator function, i.e.I A (x) = 1 for x ∈ A and zero otherwise.Additionally, since the particles do not interact with each other, the random force generated by a single particle is modeled as the divergence of the local stresses in E k .The average E[b] V (x, t) then follows as The last equality holds by Lemma 8 and Remark 9, 2).
Reynolds-kind stress tensor.The last step is the treatment of the Reynolds-kind stress term appearing in (3.2).The key here is the fact that the local disturbance velocity is scaled with ǫ and decreases fast enough at distance from the particle, while in its vicinity as well as in the particle domain it is bounded: By Lemma 4 the local fields fulfill u loc,i ∼ O(r −1 ) for r = y ≫ 0 and, since y = ǫ −1 x − c , it holds u loc ∼ O(ǫ 2 ) for x − c ≫ 0. By (A2) and (A5) we find around every particle in V a ball B r k of radius r k , containing only the particle E k and it holds: where Br k denotes the ball centered at the origin with radius rk .The first integral in the last expression is bounded by Lemma 4 and the last by Remark 9. Altogether we get E[u ′ ⊗u ′ ] V ∼ O(ǫ 4 ).

Asymptotical suspension model. By means of the macroscopic stress and force descriptions
we formulate an asymptotic suspension model to (3.2) that is consistent to the one-particle model and valid up to order O(ǫ 4 ) in the size parameter ǫ.
Theorem 11 (Asymptotical model of a suspension with weakly inertial tracer particles).Let Assumption 10 be fulfilled, and let the local behavior of any particle in the suspension be determined by Lemma 4 and Lemma 8. Then the suspension is macroscopically described up to an error of O(ǫ 4 ) by supplemented with appropriate initial and boundary conditions.The respective particle-induced stress and force are with N (x, t) being the number of particles in the averaging volume V(x, t) and with index k marking the asymptotic coefficients of the kth particle.
The form of the particle-stress Σ p in (3.5b) is known from the work by Batchelor [1], while our asymptotical approach for inertial particles (Section 2.2) give rise to an additional body force b p that is generated by small deviations of the particles' center of mass from the streamlines of the surrounding fluid.Being O(ǫ 2 ), this extra-force dominates the particle contribution to the momentum of the fluid in the case of heavy tracer particles.One should further notice that the suspension behaves like a Newtonian fluid up to an error of O(ǫ 2 ) in consistency to the underlying Assumption 10, (A7).

Special case of a suspension with ellipsoidal particles
The suspension of particles that have the same density as the surrounding viscous carrier fluid and the shape of prolate ellipsoids is often treated in literature; the structure of the particle-induced stresses is well-known in this case, see [1] and [15].In this section we illustrate the inertial particles' effects by comparing our asymptotical suspension model of Theorem 11 with those classical results.The determination of the respective analytical forms for Σ p and b p requires knowledge of the Newtonian stresses of the local disturbance velocity S[u loc,1 ] at ∂E, an explicit form of the Dirichlet conditions for u loc,1 and also of the differential equation for c 1 , which is hidden in the solvability conditions in Lemma 5. Therefore, we first study the solvability conditions, deriving the ODE for c 1 and an expression for ω 0 , also characterizing the Dirichlet conditions for u loc,1 and then use the corresponding local stresses to evaluate the integrals in (3.5).The necessary calculations are quite technical and lengthy as they are strongly related to the geometrical properties of the ellipsoids.For the sake of completeness they are summerized in the appendix.4.1.General properties.We start with introducing the quantities that characterize the geometry of arbitrary ellipsoids and discuss then the implications arising for the asymptotical framework.
An ellipsoid is given by E = DB 1 , where B 1 denotes the unit ball in R 3 and D = diag(d 1 , d 2 , d 3 ) is the diagonal matrix with the lengths of the semi axes d i > 0, i = 1, 2, 3.The surface moments appearing in the solvability conditions (2.11) in Lemma 5 can be provided as with the geometry dependent constants ζ q , q = 1, 2, 3. Appendix A.2 is dedicated to their derivation, the explicit expressions for ζ q are stated at its end.From (2.11) we see that the structure of the surface moments implies a decoupling of the conditions for the linear and angular velocities, hence we get . In our set up, g 1 ≡ 0 always holds (cf.Abbreviation 3).Consequently, the Dirichlet condition (2.8b) is presented by where the matrix-valued function A 1 = A 1 (t) is independent of y and the time-dependent vector h const 1 (t) becomes where δ = d 1 d 2 d 3 .For the definition of the geometry dependent scalars α o i and matrix M we refer to Appendix A.1.

Suspension model of weakly inertial ellipsoidal particles.
We combine the results from the one-particle asymptotics and deduce the macroscopic suspension description for arbitrarily shaped tracer ellipsoids (TE) that we even specify for prolate ellipsoids in Corollary 14.
Arbitrarily shaped ellipsoids.The particle-induced stress tensor for arbitrarily shaped ellipsoids becomes by means of (3.5b) and (4.7) The force term b p (3.5c) depends via the functions f k i (Abbreviation 3) on the choice of the mass function (2.4).With respect to the three different inertial types of tracer ellipsoids we obtain where φ = |E|N/|V| denotes the volume fraction.We use here the asymptotic approximation Theorem 11) and Assumption 2, (A8) that allows us to evaluate the gradient of the velocity at x instead of the center of mass of the corresponding particle.
Remark 12 (Macroscopic definition of the local quantities).In the derivation of b p we used the law of large numbers to replace the discrete arithmetic mean of N particles with the expectation E[.], which is in consistency with Assumption 10.This transformation was especially done to achieve a description that is easily comparable with the results presented in literature, see e.g.[1,15,19].However, it involves a small technical issue: While the ensemble average needs only the well-defined velocity fields of each particle v k 1 in V(x, t), the expression E[dv 1 /dt](x, t) presupposes the definition of an underlying random velocity field solves almost surely the corresponding ODE (4.6a).Having the meaning of the averages of the particle related quantities in mind we stick to the presented notation for reasons of readability.
Prolate ellipsoids.For a prolate ellipsoid the lengths of the semi axes satisfy It is convenient to introduce the aspect ratio a r = d 1 /d 2 > 1 and the parameter ν = (a 2 r − 1)(a 2 r + 1) −1 .Its rotational behavior can be expressed in terms of the main director p 1 , (R 0 ) ij = e i • p j as the angular velocity ω 0 (4.6b) becomes due to the orthonormality of {p 1 , p 2 , p 3 }.In combination with (4.9) is often called Jeffery's equation in remembrance of [10].In the asymptotical context the equation was derived for normal tracer ellipsoids in the work of Junk & Illner [11].As for the particle-stresses Σ p (4.8), it can be shown by a technical but straight forward calculation [23] that with δ = d 1 d 2 2 and the geometry dependent M and a i as given in Abbreviation 13 (cf.Appendix A.1).The expression holds for every particle k.Using Assumption 2 and shifting the isotropic part ))I into the pressure of the Newtonian stresses, we can state the particlestresses in the well-known form with volume fraction φ and deformation gradient tensor E[u] (see e.g.[15]), Abbreviation 13 (Geometrical parameters of prolate ellipsoids for particle-stresses).The geometrical parameters a 1 , a 2 , a 3 of (4.10) are (cf.[9,19]) In addition, where the coefficients are given by Corollary 14 (Asymptotical model of a suspension with prolate weakly inertial tracer ellipsoids).Let the assumptions of Theorem 11 be valid and let the volume fraction φ = |E|N/|V| be independent of x, then the suspension of prolate ellipsoidal particles is macroscopically described up to an error of O(ǫ 4 ) by supplemented with appropriate initial and boundary conditions.It is Σ p according to (4.10) and Remark 15 (Properties of additional inertia related forces in suspension model).1) Compared to the classical results in literature, e.g.[15], our model includes a change in the overall Newtonian stress of the fluid and additionally, in the case of heavy tracer ellipsoids, an extra body force.The source of these effects can be identified in (3.5c) to originate from the one-particle associated functions f k i for i = 1, 2, which in turn have the role of source terms in the ODEs for c k i according to the solvability conditions (2.11), clearly seen in the case of ellipsoidal particles (4.6a).This implies that the additional macroscopical effects have their origin in the deviation of the particle center of mass from the streamlines of the undisturbed fluid, which is given by This dependence is also clearly observed for ellipsoidal particles.In the general case the solvability conditions for the linear and angular velocity may not decouple, but still the f i terms are associated to the ODEsystem for the correction of particle movement.Thus throughout this text we refer to these terms in Theorem 11 as originating from the small disturbance.
2) In the special case ρ ≡ 1, we see by (2.4) that m/|E| ≡ 1 and the contribution due to the relative motion vanishes, similarly to the classical results in [1].
Remark 16 (Splitting of the inner force).As mentioned in Section 3.1, the splitting approach (3.1) for the overall inner force into a force generated by surface stresses and a volume force is only needed, since the linearity of the mean value with respect to differentiation is not generally valid for volume-based averages.Here, we want to illustrate this fact based on our results: As we see from Theorem 11, the volume average for the surface stresses yields (apart from the Newtonian stresses) only Σ p that can analytically be simplified to the classical particle-stress term of Corollary 14.This is a contribution of O(ǫ 3 ) which is dependent neither on the local coordinates y of a single particle -and thus cannot generate an effective force of O(ǫ 2 )-, nor on the mass function α mass -and thus cannot change despite different inertia models.In contrast, the average of the volume force E[b] V generates a contribution of O(ǫ 2 ) and changes accordingly to the choice of the inertia model.

Conclusion
In this paper we presented a model for a suspension of small particles, which involves inertial effects.The latter were identified as the ability of particles to deviate from the streamlines of the surrounding fluid as a consequence of different scaling of the particle momentum balance.This was achieved by an appropriate scaling of the particle mass, respectively the density ratio.To keep the results of [11], we retained the basic asymptotic approach and restricted our choice of the mass function accordingly.Having characterized the microscale behavior of the particles, we modeled the particle suspension following [1], supplementing a strategy to take additional forces into account, which have their origin in the relative motion of the particles.Afterwards we gave different models for the corresponding inertial particle regimes.These models are composed of incompressible Navier-Stokes-like equations with modified, non-Newtonian stress and force.Besides the classical part in the stress, we found a modification entering the overall viscosity of the fluid.To illustrate the general results, we applied the theory to the classical example of symmetrical ellipsoidal particles.
One practical restriction of our approach is the fact that we cannot give a recipe on how to choose an inertial regime for a given particle and fluid, since our model of inertia was formulated by means of an asymptotic behavior of the particle.For a concrete situation, one has to observe the particle motion and then, based on the intensity of deviation of the particle motion from the fluid streamlines, decide which regime is applicable.

Appendix A. Analytical statements for ellipsoidal geometry
In the appendix we summarize some fundamental results of different works for ellipsoidal particles, the individual results are from Oberbeck [17], Edwardes [6], Jeffery [10] and Junk & Illner [11].Since we use especially the solutions of Oberbeck and Jeffery in the derivation of the suspension model for ellipsoidal particles in Section 4, we provide here the relevant arguments.In Appendix A.1 From (A.2b), (A.2c) and (A.2f) it follows ) and lim s→∞ δ −1 = 0.This is an important feature of the functions, which will be used in Lemma A.2. Next, we present some asymptotic properties of these functions for y → ∞.
Proof.From the definition of λ(y), it holds and thus 0 Set d = max i d i , then the behavior of the algebraic functions δ and µ results as For the integral functions the following statement is used: Let g : (0, ∞) → R + be a continuous function with g.
Since Q was arbitrary, this is still true for Q → ∞.Another statement also needed is: Let α, a, s > 0 then 1 − (a/s + 1) −α ≤ αa/s.This can be directly concluded from the fact that the function f (s) = αa/s + (a/s + 1) −α is strictly decreasing and lim s→∞ f = 1.Consequently, The same steps lead to We summarize some results of [17,6,10] for the disturbance flow of an ellipsoid in a Stokes flow in Lemma A.2.
Lemma A.2 (Oberbeck and Jeffery solutions).Consider an ellipsoid with its geometry associated functions.Let v ∈ R 3 , A ∈ R 3×3 be a constant vector, respectively matrix, where tr(A) = 0. Define Proof.The proof of this lemma can be found in [17] for the translational boundary condition and in [10] for the linear one.For the sake of completeness, we show the relevant steps here.First, using 0 = ∆ y χ = ∆ y ω = ∆ y ψ j from (A.The last result simplifies further, if A is skew-symmetric, i.e.A = B(v) for some v ∈ R 3 .Then, with the constants c i = (d 2 j α o j + d 2 k α o k ) −1 , where (i, j, k) is a permutation of (1, 2, 3).To evaluate the integrals over the surface of an ellipsoid we transform them to the unit sphere.According to [11]  for an integrable function F : R 3 → R 3 × R 3 .Now consider the six Stokes problems in Lemma 5.According to Lemma A.2 the analytical solutions (w q , p q ) are given by (A.4) with v = e q for q = 1, 2, 3 and by (A.5) with A = B(−e q−3 ) for q = 4, 5, 6.To calculate the surface moments s q , t q , V q and W q , we first apply the transformation to the unit sphere (A.7) and afterwards use (A.6a) or (A.6b) for the corresponding moments: with v q = O • e q and v q+3 = −c q B(e q ) • D • z for q = 1, 2, 3. Since the integral of a homogeneous polynomial p(z) fulfills S 2 1 p(z)ds(z) = 0 if p has odd degree, it follows s q = 0, V q = 0, for q = 1, 2, 3, t q = 0, W q = 0, for q = 4, 5, 6.

Figure 2 . 1 .
Figure 2.1.Lagrangian description, bijective mapping between reference state and the actual time-dependent state.

Table 2 .
1.For k = 3, the density In the following S[.] always designates the Newtonian stress tensor of the corresponding solutions of the Navier-Stokes or the Stokes equations and E[.] denotes the associated deformation gradient tensor.

2 )
Let (2.11)hold.Let (u loc,i , p loc,i ) be the solution of the stationary Stokes problem (2.8a) with the Dirichlet condition (2.8b) and with the decay property (2.8e), then the integral conditions (2.8c) and (2.8d) are fulfilled.3) Let (2.11) hold.Let (u loc,i , p loc,i ) be the solution of the stationary Stokes problem (2.8a) with the integral conditions (2.8c) and (2.8d) and with the decay property (2.8e), then the Dirichlet condition (2.8b) is fulfilled.4) The linear systems (2.11) are invertible.
exclusively determined by the source term in (4.6a).The last ingredient for the suspension model are the Newtonian stresses of u loc,1 at ∂E.Since the Stokes problems in Lemma 4 are linear, the local velocity field u loc,1 is given by the linear combination of the Oberbeck u loc,1,Ob and Jeffery solutions u loc,1,Je (see Lemma A.2 stated in Appendix A.1).This yields the following stress terms that can be computed by means of the techniques provided in Appendix A.2, ∂ES[u loc,1,Ob