Unified mean-field modeling of viscous short-fiber suspensions and solid short-fiber reinforced composites

Mean-field homogenization is an established and computationally efficient method estimating the effective linear elastic behavior of composites. In view of short-fiber reinforced materials, it is important to homogenize consistently during process simulation. This paper aims to comprehensively reflect and expand the basics of mean-field homogenization of anisotropic linear viscous properties and to show the parallelism to the anisotropic linear elastic properties. In particular, the Hill–Mandel condition, which is generally independent of a specific material behavior, is revisited in the context of boundary conditions for viscous suspensions. This study is limited to isothermal conditions, linear viscous and incompressible fiber suspensions and to linear elastic solid composites, both of which consisting of isotropic phases with phase-wise constant properties. In the context of homogenization of viscous properties, the fibers are considered as rigid bodies. Based on a chosen fiber orientation state, different mean-field models are compared with each other, all of which are formulated with respect to orientation averaging. Within a consistent mean-field modeling for both fluid suspensions and solid composites, all considered methods can be recommended to be applied for fiber volume fractions up to 10%. With respect to larger, industrial-relevant, fiber volume fractions up to 20%, the (two-step) Mori–Tanaka model and the lower Hashin–Shtrikman bound are well suited.


Motivation
Lightweight components are made of short-fiber reinforced polymers due to their beneficial stiffness-toweight ratio and good formability [1,2] 1 . These components are manufactured by injection molding of fiber suspensions consisting of liquid polymer matrix and suspended fibers. During processing, the local kinematic conditions within the flow influence the evolution of the fiber orientation state which in turn influences the anisotropic elastic properties after the fluid-solid transition of the polymer matrix [3,4]. In view of flow-fiber coupling, the orientation state itself induces anisotropic viscous behavior that affects the flow, which in turn affects the evolution of fiber orientation [3,[5][6][7][8][9]. As a consequence, a physically comprehensive simulation is required to consider the anisotropic viscosity already during mold filling to better represent both the flow kinematics and the evolution of the orientation state [4,10]. The selection of a suitable modeling of the effective anisotropic viscosity is extremely important, since it is hard to be determined experimentally [11,12]. On the one hand, reorientation takes place instantaneously during the experiments, and on the other hand, many independent parameters have to be determined depending on the present material symmetry [11,12]. In this context, homogenization techniques can be used to estimate the effective anisotropic linear viscous behavior of fiber suspensions [4,5,11,12]. In addition, the adjustment to the mean-field methods established in stiffness calculations [13][14][15] enables a micromechanically consistent simulation chain improving the estimation of elastic anisotropy of manufactured parts. Moreover, mean-field models are computationally efficient compared to full-field simulations and, therefore, can be used within industrial mold-filling simulations taking the effects of flow-fiber coupling into account.

Originality of the manuscript and outline
In the framework of this study, mean-field homogenization theory and models concerning the effective linear elastic behavior of solid composites are reviewed. The aim of the manuscript is to elaborate and detail the parallels and differences in the homogenization of anisotropic linear viscous properties of short-fiber suspensions and of anisotropic linear elastic properties of solid short-fiber reinforced composites. Such a novel methodological comparison is important in order to be able to model manufacturing processes of fiber reinforced polymers consistently throughout. Previous studies [4,5,12,16,17] having advanced into this topic are extended by the following details: • The Hill-Mandel condition [18,19] is formulated in detail for fiber suspensions (see Sect The outline of the manuscript is as follows. In Sect. 1, the study is motivated and the originality is described. Furthermore, a detailed literature review is provided and the notation of the manuscript is explained. Section 2 consists of how to describe fibrous microstructures and which statistical assumptions are made for homogenization. The basic relations of mean-field homogenization for linear viscosity and linear elasticity are derived in Sect. 3. A selection of common mean-field models is given in Sect. 4. In Sect. 5, the lower and upper bounds for the linear elastic and linear viscous behavior are given. Upper estimates for the effective viscosity are introduced in Sect. 6. For measured fiber orientation data, the results of all presented mean-field models are shown and discussed in Sect. 7. The outcomes are summarized in Sect. 8, and additional material is given in Appendixes A-D. Particularly, Eshelby's single inclusion problem (SIP) [20,21] is revisited in detail regarding the parallelism between solid and fluid mechanics in the framework of a consistent tensor notation.

State of the art
Mean-field homogenization is an established method in solid mechanics to determine the effective behavior of microstructured heterogeneous materials. In general, mean-field approaches consist of bounding and estimating methods [22][23][24].
According to Hill [25], the first-order bounds based on the works of Voigt [26] and Reuss [27] are physical bounds. As a consequence, only having the constituent's properties and the volume fraction at hand, the effective elastic properties of composite materials must lie between the lower Reuss and the upper Voigt bound. Large phase contrast has an unfavorable effect on the usefulness of these bounds, since the range of the actual properties cannot be narrowed down sufficiently. Another limiting factor of first-order bounds is the restriction to isotropic effective behavior if the constituents show isotropic mechanical behavior. Based on a variational formulation with respect to an auxiliary field, namely stress polarization, Hashin and Shtrikman [28] derived second-order bounds [29][30][31]. These bounds not only account for the constituent's properties and the volume fraction, but also for the geometry of the inclusions in the sense of two-point statistics. Thus, these bounds lead to a smaller range of the effective behavior. In terms of ellipsoidal inclusions, the Hashin-Shtrikman principle is formulated by Willis [31,32] more generally. In the context of texture coefficients of crystal aggregates, the Hashin-Shtrikman bounds are derived and investigated by Lobos and Böhlke [33][34][35].
Based on Eshelby's single inclusion problem [20,21], the dilute distribution model [13,[36][37][38] estimates the effective stiffness neglecting the interaction between the inclusions. This method is limited to small inclusion volume fractions since the physical Reuss bound can be violated [39]. According to Benveniste et al. [37], the dilute distribution approach always ensures the symmetry conditions for the effective stiffness tensor.
Mori and Tanaka [40] use Eshelby's result to model the interaction between the inclusions via the effective strain which is set to the effective matrix strain [41,42]. Brylka [43] extended this concept to consider fiber orientation states using an orientation average scheme [44]. Anisotropic damage modeling of fiber reinforced composites using the Mori-Tanaka approach is done by Schemmann [45]. It should be noted that, as described by Kanaun and Levin [46], the Mori-Tanaka method coincides with the effective field method for the special case of coinciding inclusion shape and shape of the correlation function. As described by Kehrer [24] and Schemmann [45], the advantages of the Mori-Tanaka method lie in its simple and explicit structure. But on the other hand, the symmetries of the effective stiffness tensor based on the Mori-Tanaka method are not ensured generally, as summarized by Kehrer et al. [47]. The symmetries are preserved for two-phase composites [37,48] and for multiphase composites with similar shaped and aligned inclusions [49,50]. The symmetries might be lost for differently shaped inclusions within multiphase composites [49][50][51]. In addition, differently shaped anisotropic inclusions embedded in an isotropic matrix of a two-phase composite may violate the Hashin-Shtrikman bounds [50]. Kanaun and Levin [46] address the drawback of violated index symmetry regarding the Mori-Tanaka method and present a version of the above-mentioned effective field method avoiding this disadvantage.
Based on the concept of Hershey [52] and Kröner [53] for crystal aggregates, Hill [54] and Budiansky [55] developed the so-called self-consistent method for matrix-based composites. The effective elastic properties are estimated by embedding a single inclusion in an infinite large matrix having the unknown effective properties, which may be anisotropic. As described by Kanaun and Levin [46], this procedure refers to the effective medium method. In addition, it is shown that the effective medium method has the drawback of violating the physically implied major tensor symmetry for aligned inclusions. Since the effective properties are unknown, this procedure leads to an implicit formulation which has to be solved iteratively. This method is useful when the inclusion and matrix phases cannot be clearly distinguished [14]. The symmetries of the effective stiffness tensor are preserved for two-phase composites and for multiphase composites with similar shaped and aligned inclusions [37]. In addition, the mixing of the micro-and macroscale within the self-consistent method can be taken as a point of criticism [14]. A further disadvantage of the self-consistent model is the lack of direct interaction between the matrix and the inclusion [23,45]. In contrast, within the generalized self-consistent model [56,57] the inclusions are first embedded in the pure matrix material and secondly in the composite having the effective properties.
A further estimation of the effective properties of a composite is given by the differential scheme. The basic idea is to add the inclusion phase incrementally into the matrix which elastic properties depend on the current inclusion volume fraction. In every step, the dilute distribution model described above is applied. This approach first was introduced by Roscoe [58] using Einstein's estimate [59] of suspensions containing rigid spheres of equal shape in a Newtonian matrix fluid. For macroscopically isotropic solid composites consisting of isotropic phases, the differential scheme was suggested by Roscoe [60] and Boucher [61]. By using tensor calculus, McLaughlin [62] formulated the differential scheme in order to treat two-phase linear elastic anisotropic composites consisting of anisotropic phases and homogeneously distributed inclusions. It is shown that the differential scheme fulfills the Hashin-Shtrikman bounds for an isotropic distribution of spheres and for fiber reinforced composites. Norris [63] formulated a differential approach which includes the differential scheme of Roscoe [60] and Boucher [61] and the self-consistent method of Kröner [53] and Hill [54] as special cases.
In addition to the more classical models described above, there are newer developments for estimating the effective elastic properties of composites. Ponte Castañeda and Willis [64] derived explicit estimates of Hashin-Shtrikman type considering both the inclusion geometries and the spatial distribution. Hori and Nemat-Nasser [65] introduced the double-inclusion model. The basic idea is to embed an ellipsoidal inclusion in an infinite domain of matrix material. The inclusion itself again consists of two phases, one of which takes on the role of the matrix and the other the role of an inclusion. Aboutajeddine and Neale [66] developed a new version of the double-inclusion model in order to generally consider multi-phase composites. Hu and Weng [67] discuss the connection between the double-inclusion model and other mean-field estimates. Zheng and Du [68] and Du and Zheng [69] proposed the explicit interaction direct derivative model to consider multi-phase composites with interactive inclusions of different anisotropic behavior. Two-step homogenization approaches proposed by Pierard et al. [70] divide the composite into domains of which can be homogenized separately. Subsequently, a homogenization over all domains is done to estimate the effective anisotropic behavior. To conclude the overview of the mean-field homogenization of solid composites, reference is made to the literature [13,15,71] which aims to compare different estimates. A comprehensive introduction into homogenization theory can be found in standard works [14,38,39,72,73].
Early works on the topic of estimating the effective scalar viscosity of suspensions consider solid spheres immersed in a Newtonian matrix fluid. Einstein [59] derived an expression for dilute suspensions and Frankel and Acrivos [74] for the concentrated regime, respectively. Graham [75] derived a model which covers both the dilute and the concentrated regime. In order to describe the anisotropic character of fiber suspensions, the effective viscosity has to be based on a tensorial approach. In the context of anisotropic fluids, Ericksen [76] and Hand [77] derived expressions for the stress tensor with respect to non-spherical particles. Based on the transversely isotropic fluid model of Ericksen [78] only valid for dilute suspensions, Lipscomb et al. [79] derived expressions for the unknown material constants. Also, Batchelor [80] considered the effective viscosity of dilute suspensions with rodlike particles and also investigated the non-dilute regime where particle interaction is present [81]. Brenner [82] considered dilute suspensions of Brownian particles and provided expressions for different axisymmetric particle shapes. The effect of Brownian motion in the dilute regime has been studied by Hinch and Leal [83], and an extension to semi-dilute suspensions also has been carried out. Furthermore, Hinch and Leal [84] studied suspension mechanics generally and derived a simplified constitutive model. Based on the model of Hand [77], Hinch an Leal [85] derived constitutive expressions as simple interpolations between weak and strong flow in dilute regimes. For the semi-dilute regime, the model of Dinh and Armstrong [86] is applicable. In contrast, for concentrated suspensions the model of Phan-Thien and Graham [87] and Phan-Thien [88] can be used. The quantitative subdivision of the suspensions into dilute, semi-dilute and concentrated can be taken from Petrie [89]. In order to circumvent numerical difficulties of tensor-based models, Favaloro et al. [90] introduced a model covering anisotropic viscosity via a scalar function.
In the framework of estimating the effective viscosity of fiber suspensions by means of micromechanics, Bertóti and Böhlke [5] have used the Mori-Tanaka model for unidirectional pseudo-domains combined with orientation averaging to cover the fiber orientation state. Furthermore, the results are compared to the phenomenological Dinh-Armstrong model. Bertóti [12] has used the Mori-Tanaka model to determine the three viscosity parameters of Tucker [91] analytically. In addition, results based on fast Fourier transform (full-field simulations) have been used both to evaluate the mean-field approach and to improve the analytical expressions of the viscosity parameters [4,12]. In the context of fast Fourier transform, Schneider [11] has considered primal and dual formulations of cell problems of periodic suspensions and has proved existence and uniqueness theorems. Within fiber orientation evolution, Favaloro [92] considers the Reuss and Mori-Tanaka mean-field approach by using localization based on Eshelby's solution. Periodic suspension consisting of rigid unidirectional fibers in a Newtonian matrix fluid is considered in view of homogenization by Périn and Lévy [93] and governing equations for the special case of fiber-parallel flow are determined. For this configuration, Périn [94] addresses lower and upper bounds for the special case of macroscopic isotropy. Thevenin and Perreux [16] apply the Mori-Tanaka and the self-consistent model estimating the viscosity of fiber suspensions in the framework of two-phase fluid formulations and unidirectional fiber orientation. In addition, the applicability of micromechanical methods from linear elasticity is addressed. Ponte Castañeda [95] applies a Hashin-Shtrikman-type homogenization in order to model the rheological behavior of suspensions with deformable inclusions embedded in Newtonian or viscoplastic matrix fluid. Traxl et al. [17] use nonlinear homogenization in order to determine the effective viscosity of suspensions with various matrix behaviors and with different shapes of embedded inclusions or pores. The analogy to linear elasticity is mentioned, and the dilute distribution, the Mori-Tanaka model and the differential scheme are applied.

Notation
Within this manuscript, a direct tensor notation is preferred. Scalar quantities are given by lower case Latin and Greek letters, e.g., a, b, α, β, whereas vectors are denoted by lower case Latin bold letters, e.g., a, b. Second-order tensors are represented by upper case Latin and lower case Greek bold letters, e.g., A, B, σ , ε and fourth-order tensors by upper case Latin blackboard bold letters, e.g., A, B. The dyadic product is denoted by, for example, a⊗b, A⊗B, and different mappings of vectors and tensors by, e.g., Ab, AB, A[B], AB. The scalar product only being defined between tensors of equal order is indicated by, for example, a · b, A · B. In addition, the box product between second-order tensors is defined by ( A C)[B] = ABC. The second-order identity tensor is denoted by I and I S refers to the identity on symmetric second-order tensors, respectively. The projector on spherical second-order tensors is given by P 1 = I⊗I/3 and P 2 = I S − P 1 indicates the projector on deviatoric second-order tensors. In the context of named quantities and transpositions, sans serif font is used. The transposition of a second-order tensor is denoted by, for example, A T . In case of fourth-order tensors, the major transposition is indicated by, for example, A T H , the transposition of the right index pair by, for example, A T R , and the transposition of the left index pair by, for example, A T L , respectively. The inner transposition is denoted by, for example, A T I . The trace of a tensor is denoted by, for example, tr( A), and the symmetrization operator is indicated by, for example, sym( A). Volume averages are given by, for example, ψ , effective quantities are represented by, for example,ψ, and fluctuations of, for example, ψ are indicated byψ = ψ − ψ . The jump of, for example, ψ over singular surfaces is denoted by ψ = ψ + − ψ − (± describing the volumes separated by ), and the material derivative of, for example, ψ is given byψ. Further details can be found in, for example, Gurtin et al. [96] and Moakher [97].

Description of fibrous microstructures
In this section, the basics of describing fibrous microstructures are discussed to make the manuscript selfcontained. The microstructure of solids is assumed to be constant in time. For simplicity, the microstructure of viscous suspensions is considered constant at each time step in order to apply a uniform procedure. In the following, the specification of the time in the argument list is therefore omitted. This forms the basis for determining the instantaneous effective viscous behavior of suspensions without a changing microstructure during homogenization [95]. Within this study, the ergodicity hypothesis [38,39] is assumed to be valid. As a consequence, the ensemble average of an infinitely large ensemble is equal to the volume average over an infinitely large volume of one single realization [39,46]. Since infinitely large volumes cannot be considered in practice, all following micromechanical considerations refer to a finitely large statistically representative volume element (RVE) without cracks or voids [18,39]. Furthermore, phase-wise constant properties ψ γ = ψ γ are assumed for simplicity with ψ(x) representing an arbitrary microstructural quantity regarding phase γ with the spatial coordinate x. In addition, the considered microstructures are assumed to be statistically homogeneous throughout the manuscript [38,39].
The probability density function (PDF) f (x, n) [44] corresponding to one-point statistics is used to describe fibrous microstructures statistically. This function represents the probability to find a rigid fiber aligned in direction n at the spatial position x and is characterized by the following properties [44]: Furthermore, Brownian motion is not considered and the fibers are modeled neutrally buoyant without direct fiber-fiber interaction [95]. An averaged description of the fiber orientation state can be done with orientation tensors [44,98]. In general, three kinds of these tensors can be distinguished [98], with the second-order tensor N and fourth-order tensor N of the first kind being used in practice [44] The surface of the unit sphere is denoted by S with the surface element dS. Throughout the manuscript, N and N are used within the orientation average procedure discussed in Sect. 4.1. The use of only these two tensors is equivalent to an approximation of f , since all infinitely many orientation tensors are necessary for an exact description of the microstructure [44,98]. Concerning the definition of the orientation tensors of the second and third kind, the reader is referred to Kanatani [98]. These tensors are given in Karl et al. [99] consistent with the notation of this manuscript.

Homogenization of linear viscous and linear elastic properties
In the following, the steps leading to the basic equations of mean-field homogenization for solid composites and viscous suspensions are discussed in parallel. First, the constitutive relations on the micro-and mesoscale are presented and secondly, the macroscale relations are given. In addition, the Hill-Mandel condition for solid composites is reviewed and the corresponding concept is derived in case of viscous suspensions. This section concludes with the basic equations of homogenization, which allows the effective mechanical behavior to be estimated with suitable mean-field models.

Constitutive relations on the micro-and mesoscale
On the micro-and mesoscale, linear elastic behavior of the phases is assumed [18,100] and linear viscous behavior is assumed analogously regarding the phases of the suspension [5,12,16] In Eq. (3), the symmetric Cauchy stress tensor is denoted by σ , C refers to the stiffness tensor with its inverse S, named compliance tensor, and ε represents the symmetric strain tensor regarding small deformations. For viscous suspensions, V stands for the viscosity tensor, with F representing the fluidity as the inverse quantity. The symmetric part of the velocity gradient is named strain rate tensor D and the viscous stress is denoted by σ V . For the special case of incompressibility, tr( D) = 0 holds and therefore, the inversion of V is defined on the symmetric deviatoric subset SymDev [5,12], with P 2 representing the identity on SymDev. In contrast to this is the inversion of C on the symmetric subset Sym with I S as the identity on Sym. Furthermore, the pressure p has no physical meaning in the special case of incompressibility, but rather represents a reaction force or a Lagrange multiplier [5,95]. Another possibility to describe the material behavior on the micro-and mesoscale is given by the elastic strain energy density W and its complement W * for solids [18,39,100], and by the viscous dissipation W V and its complement W * V for suspensions [95] reading on the micro-and mesoscale, C = C T H applies [34]. The major symmetry [3]. Since ε and D are symmetric, C = C T R and V = V T R is usually introduced for convenience but is not required in general. As a consequence, the left minor symmetry C = C T L and V = V T L follows by the major and right minor symmetry.

Constitutive relations on the macroscale and localization tensors
Analogous to the previous section, macroscopic linear elastic behavior of the solid composite is assumed [18,31,100]. For suspensions, the assumption of macroscopic linear viscous behavior [5,12] represents the equivalent leading to the following relations regarding the homogeneous effective fields Based on Eq. (4), the effective strain energy density and its complement for solids [25,29,31,100] can be defined as follows: with the analogy regarding the effective viscous dissipation of suspensions. In order to connect the effective fields with the local fields, the following linear localization operations are used for ε and σ without eigenfields [14] following the work of Hill [18] and Walpole [100]: The localization ofD is used according to Bertóti [12]. The dependence on x is given explicitly in Eq. (7) highlighting local fields. The strain/strain rate localization tensor refers to A S , A V and B S , B V stands for the stress/viscous stress localization tensor, respectively. The index S represents the localization in solid composites and V stands for the localization in viscous suspensions. In general, the localization tensors have both minor symmetries but no major symmetry [18] and they are not invertible. Regarding the volume average, the equations A S,V = I S and B S,V = I S hold [14,18]. As described by Gross and Seelig [14], the localization tensors depend on the complete microstructure and therefore also on the material properties. By using Eqs. (3), (5) and (7), the connection between A S and B S with respect to the effective properties can be derived [18,100]. This procedure can be easily transferred to the viscous localization tensors In the remaining part of Sect. 3, the indices S and V regarding the localization tensors are omitted for convenience. In the context of the equations, it becomes clear whether the solid or viscous localization is used.

Hill-Mandel condition
First, the Hill-Mandel condition [18,19] for solid composites is reviewed to convey the aforementioned parallelism and to make the manuscript self-contained. Secondly, the Hill-Mandel condition for viscous suspensions addressed briefly by Traxl et al. [17] and Bertóti et al. [4] is derived in detail with respect to singular surfaces with a discussion of the associated boundary conditions. To the best of the authors' knowledge, this derivation has not been published yet. The Hill-Mandel condition is known to ensure an energetically correct scale transition [24,101], meaning that the macroscopic strain energy density is equal to the strain energy density of the volume-averaged microscopic quantities [101]. As described by Glüge [101], the Hill-Mandel condition is directly related to the boundary conditions of the RVE and, moreover, it can be seen as a quality criterion of the RVE. Another important point is that no assumptions need to be made concerning the constitutive modeling to fulfill the Hill-Mandel condition [14,17,24]. After Sects. 3.3.1 and 3.3.2, the Hill-Mandel condition is assumed to be fulfilled.

Hill-Mandel condition for solid composites
Based on the decomposition of the local fields σ (x), ε(x) into the volume averages σ , ε and into the local fluctuationsσ (x),ε(x), the Hill-Mandel condition reads or equivalently σ ·ε = 0. Please note that the term σ · ε represents the strain energy density in the linear elastic case. Equation (9) is valid analogously if the strain rate instead of the strain is used. Then, the Hill-Mandel condition is formulated in terms of the stress power. Regarding the derivation of Eq. (9), the reader is referred to further literature, e.g., Nemat-Nasser and Hori [39] and Kachanov and Sevostianov [102]. It can be shown that the fluctuation term σ ·ε depends on the stress vector fluctuationst and the displacement fluctuationsû on the boundary ∂ V \ of the volume V with the singular surfaces . Note that only RVEs without cracks or voids are considered, meaning that t and u are continuous over [100], which is equivalent to t = 0 and u = 0. In summary, Eq. (9) holds if any of the following four cases hold [101]: • Homogeneous displacement on the boundary ∂ V \ :û = 0 • Homogeneous stress on the boundary ∂ V \ :t = 0 • Periodic boundary conditions • Ergodic media (see, e.g., Torquato [38]) in the limit V → ∞: 2 W → σ · ε Motivated by the Hill-Mandel condition and the accessibility of ε and σ over ∂ V \ as shown in Walpole [100], the following volume averages coincide with the effective fields: It should be noted that the Hill-Mandel condition in Eq. (9) holds for arbitrary stress σ and strain ε that satisfy the equilibrium and the compatibility condition, respectively. Both fields are not required to be constitutively coupled.

Hill-Mandel condition for viscous suspensions
The Hill-Mandel condition for viscous suspensions reads or equivalently σ V ·D = 0. In order to derive Eq. (11), the local viscous dissipation W V given in Eq. (4) and its volume average W V are considered The decomposition regarding the local fields σ V (x), D(x) into the volume averages σ V , D and the local fluctuationsσ V (x),D(x) leads to the following volume-averaged viscous dissipation by using σ V = 0 and D = 0: The volume average of the fluctuation terms σ V ·D can be reformulated as follows by applying the product rule, the symmetry of σ V and D = sym grad(v) with the velocity field v To further simplify Eq. (15), the balance of linear momentum for viscous suspensions without body forces is considered in the common form without dimensions [103] ∂v In Eq. (16), fields and operations without dimensions are denoted by (·) * , and Re refers to the Reynolds number. Similar to previous studies [4,11,92,95], Stokes flow (Re 1) is considered in the context of this work for which the pressure and the viscous forces are in equilibrium. As a consequence, local accelerations due to changes in time and convective accelerations are neglected. Note that for the latter, also small spatial gradients of the flow are assumed leading to the simplified linear momentum balance grad( p) = div(σ V ). In this context, the parallels of linear elastic solids and linear viscous suspensions are addressed in the literature [38,95] as long as the kinematic constraints are the same for both suspensions and composites. Next, the splittings p = p +p and σ V = σ V +σ V are used in combination with the linearity of grad(·) and div(·) leading to the linear momentum balance in terms of averages and fluctuations Since p and σ V are constant, Eq. (17) reduces in a trivial way to the linear momentum balance for the fluctuations grad(p) = div(σ V ).
By inserting Eq. (18) into Eq. (15), the volume average of the fluctuation terms read after applying the product rule again For the special case of incompressibility in terms of velocity fluctuations div(v) = 0, Eq. (19) can be expressed as follows In the next step, the Gaussian integral theorem with singular surface [96] is applied to transform the volume integrals into surface integrals To achieve a convenient notation being compatible with the solid mechanics considerations, the stress vectors t V = σ V n and t p = −pn representing viscous and pressure forces are introduced leading to In Eq. (22), the addition of the force vectors t = σ n = (− p I + σ V )n = t p + t V is used to obtain a structure like the one commonly used in solid mechanics. The fluctuation term σ V ·D in Eq. (22)  • Homogeneous velocity on the boundary ∂ V \ :v = 0 • Homogeneous stress on the boundary ∂ V \ :t = 0 • Periodic boundary conditions • Ergodic media in the limit V → ∞: Similar to the previous Sect. 3.3.1 and motivated by the Hill-Mandel condition and the accessibility of σ V and D over ∂ V \ , the effective fields are equal to the volume averages as used by Bertóti [12] The derivation does not imply that the viscous stress and the strain rate are constitutively coupled. Equation (11) corresponds to the Hill-Mandel condition given in Traxl et al. [17] and refers to the macroscopic viscous dissipation addressed by Bertóti et al. [4].

Basic equations of homogenization
Based on the previous sections, the relations representing the effective properties of solid composites and viscous suspensions are derived. For a more comprehensive treatise of solid composites, the reader is referred to the basic literature, e.g., Nemat-Nasser and Hori [39]. By using Eqs. (10) and (23), the following relations hold [18,30,100] In Eq. (24), both the constitutive laws and the localization relations are used. The definitions of the effective stiffnessC, complianceS, viscosityV, and fluidityF are obtained by comparing the coefficients in Eq. (24). By assuming phase-wise constant properties, the following explicit expressions are obtained with c γ representing the volume fraction of phase γ [18,29,100]: The expressionV = VA is used according to Bertóti [12]. It is pointed out that no assumption regarding the material symmetry of the tensors involved has been made. Additional expressions for the effective stiffness and compliance can be derived by considering the macroscopic strain energy density combined with the localization relations as follows [14,18,29,100]: This concept is transferred to the macroscopic dissipation of viscous suspensions leading to the following alternative expressions for the effective viscosity and fluidity: It is easy to show that the definitions of the effective properties given in Eqs. (26) and (27) are equal to those given in Eq. (25) if the Hill-Mandel condition is valid. Note that for deriving the effective viscosity and fluidity, Stokes flow has been assumed. In addition, the algebraic properties of the effective tensors transfer from the corresponding local tensors by comparing the constitutive relations on the micro-and mesoscale in Sect. 3.1 with those on the macroscale in Sect. 3.2 [39].
In the following, the special case of two-phase (·) 1 and (·) 2 solid composites and viscous suspensions are considered. By using the volume averages of the localization tensors [18,100] the expressions A 1 and B 1 in Eq. (25) can be eliminated to receive equations with only one unknown localization tensor A 2 and B 2 . Note that c 1 + c 2 = 1 applies. For clarity, phase 1 is treated as the matrix material (M) and phase 2 represents the fiber material (F) from now on leading to the following basic mean-field equations for solid composites [18,29,100] and viscous suspensions: The equation forV is used according to Bertóti [12] and represents the two-fluid approach of Thevenin and Perreux [16] with the interpretation of V F as an embedded reinforcing fluid. Within the framework of the assumptions made, these equations for the effective properties are exact. For use, however, models for the localization tensors are required [14], a selection of which is presented in Sect. 4. The localization tensors correspond to those from Eq. (7) in the respective context of solid composites or viscous suspensions. In practice, the equation forC is commonly used by modeling the strain localization tensor A F . In the present study, the fibers are considered as rigid bodies with V F → ∞ and, therefore, with a vanishing fluidity F F [12].
For an alternative approach, the reader is referred to the literature [16]. As a consequence, the use of fluidity is suitable in view of numerical treatment [11,12] and Eq. (29) can be rewritten by using Eq. (8): In the last step, phase-wise constant properties VAF F = V F A FF are used. The formulation in Eq. (30) has the advantage that common models for A can be used directly. The last step consists of solving Eq. (30) forF to achieve an explicit equation forV by invertingF on SymDev It is pointed out that the infinitely large fiber viscosity V F is still present in Eq. (31). To get rid of this term leading to a finite expression, a model regarding the strain rate localization has to be used in a first step. Then, V F is pulled into A F and exploited that V −1 F = F F vanishes. This procedure results in the expressions forV given in Sect. 4.

Mean-field homogenization models
In this section, selected mean-field models are presented for solid composites and transferred to viscous suspensions in the context of the basic equations (29). In order to consider fiber orientation states within homogenization, the orientation average scheme is discussed at first. Then, different mean-field models are described briefly and no derivation is given for simplicity. Furthermore, the difference term δC = C F − C M is considered within the averaging operator A S F , since δC A S F = δCA S F holds true for phase-wise homogeneous properties [15,43]. This is done in the same way for viscous suspensions with δV = V F − V M and A V . Furthermore, the phases are assumed to be isotropic. Please note that the following mean-field equations do not hold generally for inhomogeneous and anisotropic phases.

Orientation averaging (OA)
As described above, all models are formulated uniquely with respect to the orientation average procedure [44] accounting for the fiber orientation state. In general, the orientation average of a fourth-order tensor T is defined by [44] which can be expressed explicitly as follows in case of a transversely isotropic tensor T with e 1 -direction as the axis of symmetry [44] T The coefficients b i depend on the components of T [44] The OA operator · OA refers to the averaging operator · F . Throughout this manuscript, the fiber orientation state is assumed to be given and not evolving during the homogenization procedure since, as described before, only the instantaneous effective viscosity is considered (see, e.g., [95]). Obviously, there is a restriction to N and N in Eq. (33), which in the sense of a series expansion of the PDF f represents only part of the microstructure information [44,98]. Thus, the question arises whether homogenization based on the microstructure information of the leading orientation tensors of the first kind N and N is reasonable. According to Brylka [43], the restriction to N for approximating the stiffness of fiber reinforced materials is critical when sharp textures are present. Sharp textures occur in the form of bundled short fibers and affect the reliable applicability of N along with the degree of anisotropy [43]. It is shown that using both N and N lead to smaller stiffness deviations compared to the use of N only. Müller [23] and Müller and Böhlke [104] study the use of N and N in the context of maximum entropy estimates in order to predict anisotropic elastic properties. Based on the study of Hine et al. [105] referring to the maximum entropy character of real short-fiber microstructures, Müller and Böhlke [104] conclude that N leads to reliable stiffness predictions. However, a recommendation is made by Müller and Böhlke [104] to check this on a case-by-case basis and also use N to approximate the stiffness. As a consequence, using the leading fiber orientation tensors to predict the linear viscous and elastic anisotropy in the present study is evaluated as sufficiently accurate.

Dilute distribution model (DD)
As already discussed in Sect. 1.3, the DD method neglects the interaction between the inclusions and, therefore, refers to a single inclusion embedded in an infinitely large matrix. The well-known localization tensor of the SIP is addressed in Appendix A to make the manuscript self-contained and further to formulate the SIP with respect to solid and fluid mechanics in parallel in one publication In Eq. (35), Hill's polarization tensor read P 0,S = EC −1 M and P 0,V = EV −1 M with E representing Eshelby's tensor [16,100]. For the special case of isotropic matrix behavior and spheroidal inclusions, analytical expressions exist for P 0,S [30,43,64]. The expressions for P 0,V are given in Appendix B.1 for the special case of incompressible suspensions to make the manuscript self-contained. Numerical calculation of both P 0,S and P 0,V is necessary for anisotropic matrix behavior as shown in Appendix B.2 for P 0,S . Furthermore, the analogy of bulk modulus/volume viscosity and of shear modulus/shear viscosity is pointed out. By using Eq. (35) and the basic equations (29) for the stiffness [13] and Eqs. (35) and (31) for the viscosity predictions, the equations for the effective behavior read

Mori-Tanaka model (MT)
The concept of the MT model described in Sect. 1.3 leads to the following averaged localization tensor referring to the work of Brylka [43] for both solid composites and viscous suspensions Both versions of the localization tensor A SIP are given in Eq. (35) with respect to viscosity and stiffness estimations. Based on Brylka [43] for stiffness predictions and using Eq. (37) in Eq. (31) for estimating the effective viscosity, the following expressions can be derived: The expression forV corresponds to the orientation-averaged MT model given by Bertóti and Böhlke [5] and Bertóti [12]. In addition, the idea of two-step homogenization [70] is applied to the MT model for stiffness predictions (MT-TS) as shown in Hessman et al. [15]. Then, this concept is transferred to the estimation of the effective viscosity of fiber suspensions. In the first step, the localization tensor in Eq. (37) is simplified for the unidirectional case (UD) as given in Gross and Seelig [14] and in Thevenin and Perreux [16] A MT leading to the following expressions for the effective elastic and viscous behavior [5,14] C Note that the expression forV UD differs from the formulation in Thevenin and Perreux [16] due to the rigid fiber assumption. The second step consists of averaging all UD domains in the context of OĀ In case of viscous suspensions, MT-TS refers to the procedure described by Bertóti and Böhlke [5] leading to an upper bound for the dissipation. The principle drawbacks of MT, which are addressed in Sect. 1.3, should be emphasized again at this point. Note that in the present manuscript, no drawback of MT is present due to the described modeling in Sect. 7.1.

Self-consistent model (SC)
As already discussed in Sect. 1.3, the SC model refers to a single inclusion problem with the matrix having the unknown effective properties. As a consequence, the localization tensor forC [14,15] and forV [16] read with both polarization tensors depending on the anisotropic effective stiffness or viscosity. Therefore, numerical computation ofP 0,S andP 0,V is necessary as described in Appendix B.2. Based on Eq. (42) and the basic equations (29) and (31), the implicit SC estimations read In Eq. (43),C corresponds to the expression given in Hessman et al. [15] after a small manipulation. Furthermore, Müller [23] uses the orientation-averaged SC method forC regarding PDF approximations based on Eq. (32). The equation forV follows directly based on Eq. (36) by applying the rigid fiber assumption and therefore differs from the form given in the literature [16]. It is noted that δC and δV do not depend onC andV. Analogous to the MT model, a two-step homogenization scheme is also considered for the SC model (SC-TS) and written down compactly as follows: The UD expressions in Eq. (44) can be derived by a direct use of Eq. (42) in Eqs. (29) and (31) without applying the OA procedure. Analogously, the expression forV UD differs from the formulation in Thevenin and Perreux [16] due to the rigid fiber assumption. RegardingC UD , the reader is referred to the work of Kanaun and Levin [46] or Gross and Seelig [14]. It should be noted that the drawback of violated index symmetry of this method described by Kanaun and Levin [46] in view of the effective medium method is not present here, since symmetrizing OA procedure is applied.

Differential scheme (DS)
According to the literature review in Sect. 1.3, the basic idea of the DS model is to apply the DD approach and incrementally embed inclusions in a matrix having the unknown effective properties. Using the derivation in the literature [14,62,63] with the localization tensor given in Eq. (42), an adaptation to the OA procedure leads to the following ordinary differential equations for the effective properties: To the best of the authors' knowledge, the DS model in this formulation with respect to OA has not been published yet for both stiffness and viscosity estimates. It is noted that both polarization tensors depend on the anisotropic matrix behavior and has to be computed numerically as described in Appendix B.2. Since the numerical integration of Eq. (45) starts with embedding an inclusion into the pure matrix material, the initial conditionsC(c F = 0) = C M andV(c F = 0) = V M have to be used.

Ponte Castañeda-Willis model (PCW)
As already introduced in Sect. 1.3, the PCW model [64] accounts for different inclusion geometries with different spatial distributions and, therefore, it represents a more general approach than previous models. Throughout this manuscript, all inclusions are assumed to have the same geometry described by P i and the same spatial distribution P d in order to compare the results with those of the other mean-field models discussed above. As a consequence, the localization tensor of the PCW model is given for solid composites [15] and transferred to viscous suspensions as follows: The expression forC follows directly by using Eq. (46) in Eq. (29) leading to the expression given in Hessman et al. [15] after a small algebraic manipulation. The effective viscosity expression is based on the use of Eq. (46) in Eq. (31) The PCW model follows as a special case for a constant spatial distribution of fibers from the interaction direct derivative model described in Sect. 1.3 [15,23]. For the special case of P i = P d = P 0 , meaning that the geometry and the distribution of the fibers have the same constant ellipsoidal characteristics P 0 , the PCW model equals the MT model [45,69]. Note thatV given in Eq. (47) corresponds to the orientation-averaged expression given in Ponte Castañeda [95] for rigid fibers of one single geometry.

Bounds of the effective viscous and elastic properties
In this section, the first-and second-order bounds for the effective viscous and elastic properties are presented restricted to two-phase suspensions and composites with phase-wise constant properties. In the context of viscous suspensions, the fibers are modeled rigid as addressed above. For convenience, the bounds are expressed formally by inequalities, e.g., C 1 ≤ C 2 and V 1 ≤ V 2 , representing the following quadratic forms ∀ε ∈ Sym and ∀ D ∈ SymDev [34,106] As described in Sect. 1.3, the Voigt (V) and Reuss (R) bounds [26,27] are physical bounds and depend on the volumetric composition and properties of the phases. The lower Reuss bound corresponds to the harmonic mean (B = I S ), whereas the upper Voigt bound represents the arithmetic mean (A = I S ) of the phase stiffness or viscosity [14,100]. The advantage is the ease of use, but the bounds are far apart at large phase contrast  [43], which is directly evident in the case of viscosity with a rigid fiber phase leading to macroscopic rigid behavior [14] C As can be seen directly from Eq. (49), isotropic phases always result in isotropic effective properties [14].
To overcome the drawbacks of Voigt and Reuss and to consider both the inclusion geometry and the fiber orientation additionally, the Hashin-Shtrikman two-step scheme (HS) [24,47] is applied in this work and transferred to the homogenization of linear viscous properties. The first step consists of considering an UD domain with the following effective properties for stiffness [24,47] and viscosity [5,12]: For the lower bound, P 0,S depends on C M and P 0,V on V M , whereas C F or V F is used for the upper bound. These transversely isotropic polarization tensors represent the spheroidal inclusion geometry. As already mentioned, considering rigid fibers results in an infinite upper bound for the effective viscosity [14]. In the second step, the bounds of the effective properties are determined by orientation averaging all UD domains. The domains are assumed to be distributed isotropically, leading to the following equations by using Eq. (50) regarding the UD domains for stiffness [24,47] and for viscosity: Note that P • 0 in Eq. (51) describes the isotropic distribution statistics of the domains, not the inclusion geometry. Therefore, P • 0 corresponds to the polarization tensor of a spherical inclusion with its eigenvalues λ 1 , λ 2 given as follows [32,34]: Analogously, P • 0 depends on C M or V M for the lower bound, whereas P • 0 is based on C F or V F for the upper bound. In case of stiffness, (a 1 , a 2 ) = (3K M , 2G M ) is valid with the bulk modulus K M and the shear modulus G M of the matrix. For viscous suspensions, (a 1 , a 2 ) = (3μ v , 2μ s ) holds with the volume viscosity μ v and the shear viscosity μ s of the matrix fluid. The two-step procedure discussed above is shown in Fig. 1 restricted to the finite bounds given in Eq. (51). Note that upper and lower bounds of Hashin-Shtrikman type for the effective dissipation with respect to deformable particles are discussed by Ponte Castañeda [95].

Upper estimates for the effective viscous properties
As discussed before, the upper bounds for the effective viscous behavior are infinitely large for rigid fibers. Therefore, anisotropic upper estimates (UE) based on the DS model and the lower Hashin-Shtrikman bound are suggested. To avoid the severe overestimation of viscosity by the SC and SC-TS models (see Sect. 7.2), the DS model is chosen for the upper estimates. It is emphasized that these are not true upper bounds. The upper estimates are based on the assumption that the DS model represents either the arithmetic (a), or the harmonic (h), or the geometric (g) mean. The corresponding upper estimates UEa, UEh and UEg for the effective viscosity can be determined as follows: Note that the inequality of the meansψ h ≤ψ g ≤ψ a [107] reverses to the inequality for the UĒ V UEh ≥ V UEg ≥V UEa , in a quadratic sense, based on givenV DS andV HS− . After a short calculation, a simple expression for UEg follows: In Eq. (54), exp(y ln(x)) = exp(ln(x y )) = x y has been used in a tensorial sense. Formally, the three upper estimates introduced above can be connected to an unknown homogeneous reference viscosity V 0 in the context of the Hashin-Shtrikman two-step bounding method [24,47]: Here,V UE represents eitherV UEh ,V UEg orV UEa . As described above, P 0,V refers to the inclusion geometry and P • 0,V represents the isotropic distribution of the UD domains, both of which depending on V 0 as the matrix material.

Procedure and settings
In this section, the procedure and settings are described in order to compare all mean-field models, bounds and upper estimates described above. It should be noted that DD is only applicable to dilute concentrations c F 1 and the corresponding results are listed for completeness only to show that DD is not suitable for industrial-relevant concentrations. As stated before, isotropic phases are considered [5,47] For the solid composite, polypropylene (PP) is used as matrix material (E M = 1.6 GPa, ν M = 0.4) and glass fibers are considered (E F = 73 GPa, ν F = 0.22) as given in Schürmann [108] with the Young's modulus E and Poisson's ratio ν. The study refers to a fiber aspect ratio of α = 10 and to the fiber volume fraction of c F = 0.2. This volume fraction corresponds to 41.4 wt-% with the solid mass densities ρ M = 910 kg/m 3 and ρ F = 2540 kg/m 3 [108]. The conversion between mass and volume fraction can be found in Schürmann [108] for given mass densities of the phases. All results given regarding stiffness are normalized by E M . Throughout the calculation, the viscosity is treated dimensionless, meaning that μ s is defined as the characteristic viscosity. As a result, the matrix viscosity tensor without dimension is given by V * M = 3μ v /μ s P 1 + 2P 2 . The polarization tensor P 0,S depending on isotropic matrix behavior (DD, MT, PCW, HS) is implemented for stiffness as given in the literature [30,43,64] for spheroidal inclusions and as given in Eq. (52) for spherical inclusions. Since incompressibility is assumed, the limit μ v /μ s → ∞ has to be considered within P 0,V for homogenization of linear viscous properties [5,16,17,95]. If P 0,V depends on isotropic matrix behavior, the simplified expressions for the incompressible case given in Appendix B.1 are used. If P 0,V depends on anisotropic matrix behavior (SC, DS), numerical integration is used as described for stiffness in Appendix B.2. In this context, a finite value of μ v /μ s has to be chosen to consider incompressibility. Based on the results in Appendix B.1, μ v /μ s = 10 6 is evaluated as sufficiently large. Since the study is limited to incompressible suspensions,V → P 2V P 2 is set after computation [12]. On the one hand, this procedure eliminates the large value of μ v /μ s . On the other hand, this ensures invertibility on SymDev withFV = P 2 and eliminates symmetry errors caused by numerical methods. Furthermore, tr(D) = 0 is valid only numerically in incompressible fluid solvers and in this context P 2V P 2 forces to map only the deviatoric part ofD. In addition, the viscous stress is forced to be deviatoric. The issue of non-deviatoric viscous stress in the context of defining a pressure is addressed in Karl et al. [3] and in Ponte Castañeda [95].
The equations defining the effective material behavior based on SC and DS have to be solved numerically. An explicit fourth-order Runge-Kutta method [109] with step size c F = 0.01 is used for the DS model. The equations of the SC model have been solved iteratively with matlab ® using a nonlinear equation solver based on the Levenberg-Marquardt algorithm [110].
A representative fiber orientation state is chosen as a basis of comparison. This fiber orientation state represents the final anisotropy in view of a solid fiber reinforced composite. For suspensions, on the other hand, the orientation state represents only the instantaneous anisotropy as discussed before. The chosen fiber orientation state for the OA procedure given in the following refers to the measured data N of Hessman et al. [15,111] The orientation tensor N has been determined based on N and the IBOF closure approximation [112]. Note that N given in Eq. (57) has to be transformed to the normalized Voigt notation [113] to be used conveniently in the framework of matrix operations. The elastic anisotropy is visualized via the directional dependent effective Young's modulusĒ(d), whereas the viscous anisotropy is visualized via the effective shear viscosityμ s (d, p) defined as follows: [5,114] 1 WithinĒ, d represents the tensile direction and withinμ s , d refers to the shear direction lying in the shear plane defined by the normal vector p [5,114].

Results and discussion
In Fig. 3, the viscous anisotropy is shown for all mean-field models and bounds introduced above. The effect of the suspended fibers with respect to an increased viscosity and an induced anisotropy is distinct. By comparing the different plots, both the degree of anisotropy and the increase in the viscosity strongly depend on the mean-field model in use. The Reuss bound shown in Fig. 3a predicts an isotropic effective viscosity, whereas DD, HS− and MT/MT-TS predict a weakly anisotropic behavior for the chosen orientation. Optically, MT and MT-TS cannot be distinguished. One can see that the effective viscosity estimation based on DD violates HS− for the chosen parameters which shows that the validity of the dilute suspension is not given. Higher viscosity is predicted by SC-TS and DS, which is due to the fact that the inclusions are embedded in the effective material. This behavior is known from the homogenization of stiffness [13,15,115]. In addition, the anisotropy can be seen more clearly in SC-TS and DS. In Fig. 3b, the PCW model predictions are shown with respect to HS− and the pure matrix fluid. Three different aspect ratios α d = 3, α d = 8, α d = 20 for the ellipsoidal distribution has been chosen arbitrarily to be both below and above the inclusion aspect ratio α i = α = 10 similar to Hessman et al. [15]. The results show that the effective viscosity based on PCW does not change significantly with α d and moreover they almost coincide with the results of the MT model. The results of SC are omitted in Fig. 3a since SC tends to percolate quickly and macroscopically predicts a very high viscosity even at moderate fiber volume fractions [95]. This behavior is known from stiffness homogenization [14] and is even more distinct in the viscous case. By plotting the maximum in-plane viscosity against the fiber volume fraction, Fig. 4a addresses the percolation issue. Typical fiber mass fractions 30 wt-% and 40 wt-% with respect to engineering practice are additionally marked in order to provide reference points for the relevant ranges to practice. The results show that the percolation of SC is distinct. With respect to SC-TS, DS predicts a larger maximum shear viscosity for moderate fiber volume fractions. This behavior is reversed shortly before c F = 0.5 and SC-TS shows percolation, while DS predicts the macroscopic rigidity only in the limit value c F → 1, as described by Gross and Seelig [14]. For comparison, the model of Phan-Thien and Graham [87] (PTG) being valid for concentrated suspensions (c F α > 1) [10] is shown. The results indicate that there is only a small range in which PTG must be used, since it violates HS− and even R for low volume fractions and PTG starts to percolate even before SC-TS for higher volume fractions. In addition, Fig. 4a illustrates that the prediction of anisotropic viscosity is highly model-dependent and can lead to large differences regarding the estimated effective viscosity. Based on the results, mean-field models in which the fibers are embedded in the effective medium (SC, SC-TS, DS) quickly lead to large differences in the anisotropic viscosity regarding the other methods (DD, HS−, MT, PCW). This is due to the fact that these effective medium methods have poles regarding the effective behavior for the case of rigid fibers. To avoid a strong overestimation, the SC method must only be used for c F < 0.1 based on Fig. 4a. It has also been observed that in this range the iterative solution procedure of SC converges well. The deviation of DS and SC-TS from HS−, MT and PCW is moderate for c F < 0.2 and, therefore, all models except SC may be applied in this range without overestimating the effective viscosity due to being too close to the poles shown. As already discussed in Sect. 7.1, DD should not be used since industrial-relevant fiber concentrations do not correspond to dilute suspensions. In addition, R is not capable of predicting anisotropic viscous behavior and, therefore, is evaluated to be not precise enough. Fig. 4a also shows that DS is an appropriate choice for c F < 0.2 in terms of the upper estimates, which are shown in Fig. 4b for the present orientation state. Since DS represents a different mean value, but is constant analogous to HS−, the known order of the mean values is reversed for the upper estimates (see Sect. 6). It should be noted that the upper estimates lying between SC and DS have to be seen as different bounding methods, meaning that the difference in the anisotropic viscosity is desired. Comparing the mean-field predicted anisotropic viscosity with experimental data is not entirely possible as described by Schneider [11]. On the one hand, many independent experiments have to be carried out in order to determine the anisotropic state. On the other hand, any experiment will cause fiber reorientation and thus the orientation state to be measured no longer matches the given orientation state. Schneider [11] also mentions the problem of analytical methods not to be accurate for large fiber volume fractions. This supports the previously given recommendation to use certain mean-field models only in a limited way with respect to the fiber volume fraction. The full-field simulations carried out by Bertóti et al. [4] indicate that a mean-field model based on MT is capable of estimating the effective viscosity of short-fiber suspensions up to c F = 0.2. The special case of spherical inclusions in the context of viscosity homogenization is addressed in Appendix C.
Regarding to the use in flow simulations with fiber-induced anisotropy, all explicit models for viscosity are well suited in view of computational cost. The SC-TS model is only slightly more disadvantageous, as the UD state only has to be solved iteratively once and then the OA procedure is carried out explicitly using the fiber orientation evolving locally within the suspension. On the other hand, SC requires the iterative solution of an implicit equation and DS goes along with numerical integration, both for every local orientation state within the flow. Furthermore, by assuming rigid fibers in viscous matrix fluid, SC may suffer from convergence issues as addressed by Thevenin and Perreux [16].
The results described above are based on a constant matrix viscosity referring to Newtonian behavior. In view of generalized Newtonian matrix behavior [10], the given mean-field methods can also be used estimating the effective viscosity depending on, for example, the shear rate. This topic is discussed briefly in Appendix D in terms of MT.
The effective stiffness regarding different mean-field models and bounds is shown in Fig. 5. Analogous to the viscosity results, the Voigt and Reuss bounds indicate how the stiffness is isotropically increased when isotropic phases are assumed. The Voigt bound corresponds to a circle with radiusĒ/E M = 10.08 and is omitted since a more suitable upper bound is given with HS+. In Fig. 5a, it can be seen that DD almost coincide with HS− for ϕ ≈ π/2 and ϕ ≈ 3π/2. As already discussed, DD may not be used in practice since industrial-relevant fiber concentrations violate the dilute distribution assumption. Furthermore, MT/MT-TS, DS, and SC/SC-TS are close to HS− compared to HS+. The predicted stiffness increases in this order, which is known from the literature [14]. Unlike viscosity, the percolation of SC is not distinct for the chosen fiber volume fraction, which is why both the SC and SC-TS estimations are shown in Fig. 5a. In addition, compared to viscosity, MT and MT-TS can be clearly distinguished visually. It is noticeable that each of the TS models predict lower stiffness than the respective direct model, which coincides with the results of Hessman et al. [15] addressing MT and MT-TS. In Fig. 5b, the results regarding the PCW model for three different spatial distribution aspect ratios α d are shown together with HS+ and HS−. Compared to the viscosity homogenization results, the spatial distribution clearly influences the effective elastic anisotropy. For the considered specification of the composite, the compatibility with HS+ and HS− shows that the chosen spatial distributions give plausible results. As described by Hessman et al. [15], the plausibility of the PCW-estimated elastic anisotropy depends on how the fibers behave in relation to the matrix cell, which altogether depends on the fiber orientation and on the fiber geometry. The special case of spherical inclusions in the context of stiffness homogenization is omitted, since this is well discussed in, for example, Gross and Seelig [14].
The recent publication of Hessman et al. [15] also provides full-field simulation results and experimental data in the context of short-fiber reinforced polymers. Regarding the anisotropy plots therein, the full-field simulation results and the measured direction-dependent stiffness almost coincide with the estimations based on MT-TS and MT. These results lie roughly between HS− and SC, whereas HS− is more suitable for estimating the effective viscosity as described above. It turns out that the PCW does not correctly predict the anisotropy for the chosen parameters and that the difference between HS+ and the full-field and experimental results is large.

Summary and conclusions
Within the scope of this manuscript, the basics of homogenization for linear elastic fiber reinforced composites are reviewed and transferred to linear viscous fiber suspensions. Particularly, the Hill-Mandel condition [18,19] as an important requirement for mean-field modeling is derived in detail for viscous suspensions and, thus, previous studies [4,17] are supplemented. In this context, it is shown that the assumption of Stokes flow [4,11,92,95] is purposeful, leading to four cases for which the viscous Hill-Mandel condition is valid. These four cases can be directly linked to solid composite considerations in the literature [45,101]. In addition, incompressibility in the context of viscous suspensions is assumed and temperature effects are neglected. A wide range of common mean-field models for the effective stiffness [14,15] are reviewed and transferred estimating the effective viscosity, in order to provide a comprehensive basis for a uniform mean-field procedure. All of the discussed mean-field models for both suspensions and composites are formulated in the context of orientation averaging [44] taking into account leading fiber orientation tensors [44,98]. The SIP is revisited in order to formally support the parallelism of homogenization of stiffness and of viscosity. In addition to meanfield modeling, the first-order bounds of Voigt and Reuss [26,27] and the second-order Hashin-Shtrikman two-step scheme [24,47] are reviewed and transferred to viscous suspensions. The analysis shows that the upper bounds lead to rigid behavior on the macroscale for viscous suspensions, since the fiber phase is modeled rigid [14]. In order to provide upper estimates for the effective viscosity, the differential scheme results are seen as different means between the lower Hashin-Shtrikman bound and the upper estimates. For the latter, explicit expressions are derived.
Against the background of an exemplary selected fiber orientation state [15,111], all mean-field models, bounds and upper estimates are computed and compared. The results are limited to a fiber aspect ratio α = 10 and a fiber volume fraction c F = 0.2. To be as general as possible, the results are presented without dimension. For the solid, PP is used as matrix and glass as fiber material [24,47], whereas the viscosity is considered from the beginning without dimension. The results regarding viscosity show that SC tends to overestimate the effective behavior strongly and, therefore, only can be recommended to be used for c F < 0.1. SC-TS and DS also embedding fibers in the effective medium can be used up to c F = 0.2 without large deviations regarding, e.g., MT, HS− and PCW. This is a typical range in industrial applications and corresponds to ≈ 40 wt-% for the chosen materials. Applying DD is not recommended since the dilute distribution assumption is not given for practice-relevant fiber concentrations. By considering the full-field and experimental results of Hessman et al. [15] regarding the effective stiffness of short-fiber reinforced composites, DS, SC and SC-TS overestimate the effective behavior. In view of solid composites, MT, MT-TS and HS− are evaluated positively by considering the results of Hessman et al. [15] and PCW is not recommended for the chosen parameters due to differences in view of the predicted elastic anisotropy.
In conclusion, for c F < 0.1 the deviations between the considered models are small and all can be recommended to be applied in a micromechanically consistent simulation chain. The overestimating property of SC, SC-TS and DS is more distinct for suspensions than for solid composites based on the considered materials and orientation state. As a consequence, for industrially relevant fiber contents up to c F = 0.2 (or ≈ 40 wt-%) applying MT, MT-TS and HS− is recommended for estimating the effective behavior of both viscous suspensions and solid composites consistently. In this range, HS− can be used as a lower bound for both suspensions and composites, whereas UE, DS or SC-TS can be used to estimate the maximum viscosity in this range and HS+ for the upper bound of the solid composite.
In view of flow-fiber coupled generation of fiber orientation data based on fiber-induced anisotropic viscous behavior [3,[6][7][8][9], a uniform mean-field modeling can be applied based on the present work. In case of viscous fiber suspensions, the presented mean-field models can be applied to replace well-known fiber suspension models from the literature, e.g., [79,86,87]. In this context, the fiber-induced anisotropic viscous behavior influences the orientation evolution in a coupled sense and this causes a direct impact on the final anisotropic properties of the solid composite. With respect to injection molding simulations, computation time is a limiting factor, since in every time step homogenization of the linear viscous properties has to be done locally on every grid cell depending on the evolving fiber orientation state. Particularly, explicit mean-field models are computationally efficient and even implicit models can be implemented as two-step procedures saving computation time, since the iterative solution only has to be done once for UD. In this context, mean-field models are superior to full-field approaches, e.g., fast Fourier transform [4,11,116,117], since the latter require a fully resolved microstructure depending on the local orientation state. Especially in injection molding simulation or in highly heterogeneous microstructures where no single averaged orientation state can be considered, full-field approaches are computationally too expensive to be applied in engineering practice. A disadvantage of mean-field methods is that, compared to full-field approaches, only averaged fields can be computed. Therefore, it is not possible to resolve the strongly inhomogeneous fields within the composite or suspension in the area around the fibers. Since full-field approaches are based on a more comprehensive description of the microstructure, the effective viscosity or stiffness is estimated less accurately with mean-field methods. Finally, it should be noted that the drawbacks of some mean-field methods described in Sect. 1.3 are not present in this work, since only two homogeneous, isotropic phases and only rotationally symmetric fibers are considered.
Acknowledgements The partial financial support by the Friedrich und Elisabeth Boysen-Stiftung (Grant BOY-163) is gratefully acknowledged.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Author contribution TK contributed to conceptualization, methodology, formal analysis, software, investigation, data curation, visualization, writing-original draft, and writing-review and editing. TB performed conceptualization, discussion, supervision, funding acquisition, and writing-review and editing. All authors read and approved the manuscript.
Funding Open Access funding enabled and organized by Projekt DEAL.

Conflict of interest
The authors declare that they have no conflict of interest.

Appendix A: Principle of equivalent eigenstrain and extension to the principle of equivalent eigenstrain rate
In this section, the SIP is addressed in view of solid and fluid mechanics considering both material classes in parallel for the first time in a consistent tensor notation. To the best of the authors' knowledge this parallelism has not been published yet. Regarding solid mechanics, the reader is referred Eshelby's work [20,21] and to, e.g., Nemat-Nasser and Hori [39]. The SIP with respect to immiscible fluids refers to the work of Wetzel and Tucker [118] and Bilby et al. [119,120] and is called principle of equivalent eigenstrain rate throughout this manuscript. The following steps formally demonstrate that well-known mean-field models of elasticity can be applied for viscous suspensions. Inhomogeneities (I) with C I embedded in the matrix material C M are replaced by a homogeneous reference material C 0 and a corresponding eigenstrain ε * . The volume of the composite V is considered to fulfill the implication of Eq. (10) and the static balance of linear momentum. As shown in the upper part of Fig. 6,  Fig. 6 Principle of equivalent eigenstrain (upper part) and principle of equivalent eigenstrain rate (lower part) to consider inhomogeneities embedded in a matrix (own sketch based on Gross and Seelig [14] and Brylka [43]) the effective strainε acts on the boundary ∂ V . For viscosity, all quantities are introduced analogously with the interpretations following from the lower part of Fig. 6. The viscosity of the inclusion is not specified in detail. In general, the problem can consist of two immiscible fluids. Then, during the derivation, the result of Eshelby's single inclusion problem [20,21] is used. First, the homogeneous properties C 0 and V 0 are added to the constitutive equations. Note that the index 0 indicates fields within this homogeneous material. The derived fields τ and τ V refer to the stress polarization and the viscous stress polarization [95] representing the difference between the true (viscous) stress and the (viscous) stress induced by the true strain (rate) in the homogeneous reference material In the next step, the true strain (rate) can be expressed as follows by using Eq. (A1): By splitting the (viscous) stress field into the effective and a fluctuating part, the splitting of ε and D can be derived Note that the true fluctuations consist of the fluctuation in the homogeneous reference material and the introduced eigenfields representing the inhomogeneities. By using the expression of the stress polarizations, the eigenfields can be written as follows: The special case of C 0 = C M and V 0 = V M leads to eigenfields only present in the volume of the inclusions In the next step, the fields ε and D are considered inside the inclusion and Eshelby's result is applied to connect the fluctuations inside the inclusions with the corresponding eigenfield via the fourth-order tensor E In Eq. (A6), the nonvanishing field given in Eq. (A5) has been used. It is pointed out that E depends on either the solid or the viscous matrix behavior, depending on whether homogenizing the stiffness or the viscosity. Based on Eq. (A6), the corresponding localization tensors can be derived as follows: Based on Eq. (A6), the fluctuations inside the inclusions read which are directly linked to the stress polarization fields for leading toε In this context, the polarization tensors P 0,S = EC −1 M and P 0,V = EV −1 M are defined mapping the (viscous) stress polarization to the strain (rate) fluctuation in the inclusion. The special case of incompressibility regarding P 0,V is considered in Sec. B.1 to make the manuscript self-contained in view of the methods in use.

Appendix B: Polarization tensor
B.1: Special case of incompressibility in the context of viscosity homogenization As described in Sect. 7.1, the polarization tensor P 0,V for homogenization of incompressible fiber suspensions has to be determined with respect to the limit μ v /μ s → ∞. For the special case of isotropic matrix behavior, this can be done analytically leading to simplified expressions. Regarding spherical inclusions, Eq. (52) reduces to For fiber-like spheroidal inclusions with aspect ratio α > 1 in e 1 -direction, the polarization tensor is used in the following normalized Voigt notation [30,43] with the matrix entries given in Ponte Castañeda and Willis [64] P After applying μ v /μ s → ∞ to the matrix entries above, the simplified expressions given below follow for the case of an spheroidal inclusion in an incompressible matrix fluid: (B13) The incompressible Eshelby tensor E is addressed in the literature [12,17] and also other expressions for the incompressible P 0,V [95,121]. It is pointed out that incompressibility is not considered in view of solid composites. Therefore, the expressions given in Ponte Castañeda and Willis [64] for the matrix entries of P 0,S are used with the elastic constants E M and ν M given in Sect. 7.1. Note that the viscous P 0,V follows directly by substituting the bulk modulus K M with μ v and the shear modulus G M with μ s . Since the expressions in Eq. (B13) cannot be applied for anisotropic matrix behavior, numerical computation of the incompressible viscous polarization tensor is necessary. This procedure is addressed in the following Sect. B.2. In this context, a sufficiently large μ v has to be determined ensuring incompressibility based on the following error measure error 1 (μ v , μ s , α) similar to Hessman et al. [15] error 1 = P 0,V − P 0,V,exact P 0,V,exact .

B.2: Numerical computation for anisotropic matrix behavior
In case of anisotropic elastic matrix behavior, the polarization tensor P 0,S has to be computed numerically. In this context, the following integral over the unit sphere S [31,32] is approximated by a quadrature rule: Note that for homogenization of linear viscous properties, C has to be replaced by V in the equations given in this section representing P 0,V . In Eq. (B15), the quadrature weights are denoted by w i and N int stands for the number of integration points on S. All other tensors are defined by [24,31,32] H = I S K −1 (n⊗n) I S , K = C T I [n⊗n] = C ik jl n k n l , Z = α −1 e 1 ⊗e 1 + e 2 ⊗e 2 + e 3 ⊗e 3 . (B16) In this study, the quadrature points and weights (N int = 5810) are computed according to Lebedev and Laikov [122] with the implementation of Parrish [123]. In order to show that the amount of integration points is sufficiently large, the following error measure error 2 (N int , α) is defined similar to Hessman et al. [15]: error 2 = P 0,S − P 0,S,exact P 0,S,exact .
(B17) In Eq. (B17), P 0,S,exact refers to the analytical expression given in the literature, e.g., [43,64] with C M given in Sect. 7.1, whereas P 0,S represents the numerical integration shown in Eq. (B15). Note that the special case of spheroidal inclusions in an isotropic matrix is considered for simplicity. The error depending on the fiber aspect ratio α and the number of integration points N int is shown in Fig. 7b. The results show that N int = 5810 for α = 10 leads to a well approximated polarization tensor. In general, the error increases with increasing α since large fiber aspect ratios represent a concentration on S to be resolved by the numerical grid.

Appendix D: Generalized Newtonian matrix fluid
Shear rate depending matrix behavior in the context of generalized Newtonian fluids [10,17] is addressed in this section. For simplicity, only the example of a power-law behavior is considered as used in the study of Mezi et al. [6]. By introducing a critical shear rateγ c representing the transition point from Newtonian to non-Newtonian behavior, the viscosity tensor of the matrix fluid is given as follows with respect to the power law [6,10]: It is postulated that the shear rate dependence affects only V M and that the general structure of the basic equations of homogenization remain unchanged. Furthermore, the volume viscosity does not change since incompressibility is assumed. The shear rateγ represents the effective shear rate of the matrix following Thevenin and Perreux [16] The average strain rate of the matrix D M can be expressed in terms ofD by usingD = c M D M + c F D F . By considering the special case of rigid fibers, the average strain rate of the fibers D F vanishes andγ readṡ γ = D ·D Similar to Mezi et al. [6],γ c = 1 s −1 is chosen as the transition point and the shear thinning behavior refers to the exponent b = 0.3 and the parameter a = 1 Pas. For the reference shear rate, the valueγ 0 = 1 s −1 is chosen. In Fig. 9a, the function μ s (γ ) is shown normalized with the constant Newtonian matrix viscosity μ s with respect toγ /γ 0 . After the transition pointγ /γ 0 = 1, the shear thinning behavior is clearly visible. In Fig. 9b, the normalized effective viscosity of the fiber suspension is shown based on MT with respect to selected shear rates marked in Fig. 9a. The parameters and fiber orientation state refer to Sect. 7.1. The results show that the anisotropy does not change while the viscosity magnitude decreases with increasing shear rate. It is noted that μ s (γ ) given in Eq. (D20) is only an example and arbitrarily complex functions can be used also depending on temperature and curing [10,125]. Further information about different non-Newtonian matrix fluids within homogenization of viscous fiber suspensions can be found in, for example, Thevenin and Perreux [16] and Traxl et al. [17].