Modeling the shape dynamics of suspensions of permeable ellipsoidal particles

A dynamic two-scale model is developed for describing the mechanical behavior of suspensions of permeable ellipsoidal particles. The particle dynamics in the proposed model is described in terms of particle positions as well as conformation tensors that capture their size, shape, and orientation. Using non-equilibrium thermodynamics, the macroscopic fluid-dynamics and the particle dynamics on the microstructural level are mutually coupled in a consistent manner. So doing, the link between the macroscopic behavior, e.g. stresses, and the dynamics of the microstructure, e.g. particle shape and size, is established. Finally, the model is cast into a form in which the shape tensor is split into its volumetric and isochoric shape contributions, making it possible to model particles with both shape-preserving size-changes (e.g. swellable particles) and volume-preserving shape-changes (e.g. incompressible yet deformable particles). The size-shape model distinguishes itself in unifying prior knowledge of purely-shape models with that of purely-size models by appropriate choices of the Helmholtz free energy and the generalized mobility.


Introduction
A wide variety of applications nowadays relies on materials where their overall properties can be tailored to meet specific requirements. Soft, permeable particle suspensions provide the versatility required to achieve exactly this purpose, making them particularly useful in paints and inks [1,2] , pharmaceuticals and cosmetics [3,4] , and foods [5] . The fascinating properties of the overall suspension emanate primarily from the properties of the individual particles. On the one hand, the elasticity of the supporting network of the individual particle gives rise to its elastic behavior. The flow of the viscous suspending solvent through this elastic network, on the other hand, results in its viscoelastic behavior.
The rich behavior of permeable-particle suspensions emerges from the fact that permeable particles can undergo size and shape changes in response to different stimuli. For instance, permeable particles in a sufficiently-jammed state undergo rate-dependent volume changes as the viscous background solvent is expelled from the interior of the particle [6,7] . The shape changes in permeable particles are induced by steric effects in concentrated suspensions as a particle impinges against neighboring particles [7] . While elastic shape-changes can be accounted for through soft-interaction potentials [8][9][10] , the effect of the viscous background solvent on both the shape and size dynamics requires accounting for the particle internal degrees of freedom explicitly.
The dynamic two-scale model developed by Hütter et al. [11] presents a new class of models that provide insight about the degrees of freedom of permeable particles. The model accounts for the rate-dependent size change of the particles by treating the particle size as a separate degree of freedom. The developed model, however, focuses on spherical particles for which the particle geometry is described by its radius. Using this model, we have highlighted the effect of the size dynamics on the equilibrium properties [6] , the flow properties [12] , and the stress-relaxation behavior [13] of permeable-particle systems. This paper aims at generalizing the model developed by Hütter et al. [11] towards non-spherical particles. This requires modeling the particle shape explicitly, which naturally includes both the size and shape dynamics.
An essential step in the model development is the suitable choice of a variable that describes the particle shape. Several morphology measures have been introduced in the past, particularly in the field of mixing of immiscible fluids. The interfacial tensor in the Doi-Ohta model [14,15] provides an average description of the morphology of the entire dispersed phase in emulsions and immiscible polymer blends. In the present paper, a tensor is used for each individual particle, due to our interest in a many-particle description, as described in the following. For concentrated permeable-particle suspensions, tracking the exact particle surface can be computationally expensive. Therefore, we consider particles of ellipsoidal shape for simplicity. Ellipsoids cover a wide va-riety of shapes, ranging from platelets to spheres to threads [16,17] . Ellipsoidal droplet shapes are commonly used as an approximation of the microstructure of fluid blends. In fluid mixtures, this coarse-grained description depicts the local features of the microstructure, such as the size, shape, and orientation. An ellipsoid centered at the origin of the coordinate system is described by a second-rank tensor whose eigenvalues are the square of the (inverse) semi-axes lengths of the ellipsoid and whose eigenvectors give the direction of the principal axes of the ellipsoid. In general, the surface of an ellipsoid at a position Q of the coordinate system and whose axes are not necessarily aligned with the coordinate system is described by ( − ) ⋅ ⋅ ( − ) = 1 , with positive definite tensor S . Several efforts have been dedicated to developing models describing the evolution of such tensors in response to deformation fields. For instance, the Maffettone and Minale (MM) model [18] describes the dynamics of ellipsoid droplets in a general flow field, based on a phenomenological description of the driving force and the relaxation mechanism. It has been shown experimentally and numerically how the internal morphology of a system affects its overall properties [19,20] . Rheology can be even a measure for probing the morphology [21] . Iza and Bousmina have highlighted the degree of complexity of the morphologies developed in a fluid-fluid mixture subjected to shear and upon cessation of flow [19] . This elucidates the importance of accurately and consistently describing the stress as a function of the microstructure. In order to achieve this, the evolution of the ellipsoidal tensor in affine deformation is used as the starting point. For purely affine deformation, each point in the ellipsoid is subjected to a velocity gradient = ( ∇ ) T , where is the applied flow field. Consequently, the evolution of the ellipsoid is given by ̇ = − T ⋅ − ⋅ , which is a lower-convected time derivative [17] . The evolution of the ellipsoid in affine deformation can be equivalently described in terms of the inverse tensor = −1 . The evolution of T is upper-convected, that is ̇ = ⋅ + ⋅ T [17] . While the evolution of the particle-related ellipsoidal tensors in this paper is upper convected in nature, the interfacial tensor in the Doi-Ohta model [14,15] has lower-convected characteristics.
In this work, the dynamic two-scale model developed by Hütter et al. [11] for permeable particles is extended to also account for the mechanics and dynamics of the particle shape. Each particle is described with an ellipsoidal tensor. Non-equilibrium thermodynamics, namely the general equation for the non-equilibrium reversible-irreversible coupling (GENERIC) [22][23][24] , is used to ensure that the developed model is thermodynamically consistent. The developed model is expressed in the form of stochastic differential equations, that are suitable for particlebased simulations, i.e. Brownian dynamics simulations. This paper is organized as follows. In Section 2 , the weak formulation of GENERIC is briefly described. This is used in Section 3 to develop a dynamic two-scale model for permeable particles that undergo shape and size changes. In Section 4 , the model is presented in a form suitable for particle-based simulations. In Section 5 , the model is split into purely-size and purely-shape dynamics, and applied to the case of noninteracting ellipsoidal particles. Finally, the paper is concluded with a discussion in Section 6 .

Methods: Weak formulation of GENERIC
The general equation for the non-equilibrium reversible-irreversible coupling (GENERIC) [22][23][24] is exploited in this paper in order to develop a model that mutually couples mesoscopic degrees of freedom to macroscopic ones in a consistent manner. In this work, the weak formulation of GENERIC formulated in [11,25] is used, summarized in the following. For a closed system, the weak formulation of GENERIC [11,25] imposes the following conditions on the reversible (rev) and irreversible (irr) contributions to the time evolution of the energy E and the entropy S , respectively, The conditions (1) depict the following features of the system. On the one hand, the energy and entropy remain unaffected by the reversible dynamics, which is captured by (1a) and (1b) , respectively. On the other hand, the irreversible dynamics does not affect the total energy and leads to non-negative entropy changes, as given by (1c) and (1d) , respectively. It is noteworthy that, although the weak formulation of GENERIC (1) is less restrictive than its full formulation, it retains many of the essential features. Particularly, the degeneracy conditions in the full GENERIC formulation are reflected in conditions (1b) and (1c) . Using the chain rule, conditions (1) have implications on the evolution of the system variables x . The chain rule for a general functional A is given by where A / x I is a functional derivative of A with respect to x I , z is the integration variable, and the summation runs over all variables in x . In the following, the conditions on the energy and entropy, (1), are used to develop a model for systems of permeable particles. This is achieved by, first choosing a sufficient set of variables describing the system, and second specifying the functionals of energy and entropy in terms of the chosen variables.

Model development
Similar to the model developed earlier in [11] , a two-scale model is developed in this work, where both scales are mutually coupled. For instance, a deformation applied on the macroscopic level distorts the microstructure, this results in unbalanced interactions between the particles, which in turn give rise to macroscopic stresses. In this section, we derive a model that consistently couples both scales, particularly by providing a constitutive relation for the stress in terms of mesoscopic variables.

Choice of variables
For the two-scale model described for spherical particles [11] , the macroscopic level, i.e. the solvent-particle system, is treated as a nonisothermal fluid. The macroscopic variables are, hence, the mass density ( r ), the momentum density ( ) = ( )∕ , where v is the macroscopic velocity field, and the temperature field ϑ( r ). In all these variables, r is the macroscopic position. On the mesoscopic level, overdamped particle dynamics is considered. That is the particle velocities relax to the equilibrium distribution much faster than the time required for the applied deformation to cause a significant change in velocity. The reader is referred to [11] for more detail. Mesoscopically, each particle i is, hence, described with the position of its center Q i measured relative to r , and also by a tensor that captures the shape of the particle T i . The latter is introduced in place of the particle radius in [11] , in order to describe the shape and size of the particle.
For practical reasons, a distribution function p of the mesoscopic states of all particle positions { } =1 , …, and shape tensors { } =1 , …, is used a dynamic variable for the mesoscopic level of description. To account for inhomogeneous situations, the distribution function is made dependent on the macroscopic position, leading to = ( , { } , { }) , where {…} denotes the collection of variables for all particles, i.e. including all terms from = 1 until = . Averages over mesoscopic states, denoted by ⟨ · ⟩, can be conveniently described in terms of p as where n denotes the number of N -particle systems per unit volume, and d Q d T denotes the volume element in the high-dimensional space of all particle position vectors and shape tensors. In summary, the full set of independent variables describing both levels of description, macroscopic and mesoscopic, respectively, is given as

Energy and entropy functionals
The model development involves the specification of the energy and entropy functionals in terms of the dynamic variables x . This step, however, is not essential at this point. Once the model has been developed in a form that is as general as possible, the energy and entropy functionals can be specified for the particles of interest. The only assumption that is made at this point is that the momentum density only affects the kinetic energy, which is part of the total energy of the system. The internal energy int and the entropy S , that both remain unspecified at this point, do not depend on the momentum density. In other words, the energy and entropy functionals are given by

General form of the evolution equations
In order to derive the full evolution equations of each variable in x given in (5) , one can make use of prior knowledge of some characteristics of these evolution equations, in particular, in terms of kinematics. In contrast, the constitutive relation for the stress tensor and the irreversible contributions to the evolution equations are not known a priori. These contributions shall be worked out later.
Concerning notation, throughout this paper, particle indices are denoted by lower-case Latin symbols i, j, k , … , while Cartesian indices are given in lower-case Greek symbols, e.g. , , , … The Einstein notation is used for the Cartesian components, i.e. the summation notation is suppressed for the repeated Greek subscripts. The general form of the evolution equations of the dynamic variables (5) can be written as where is the total macroscopic stress tensor, Θ is the temperature change, and ̇ and ̇ stand for the rates of change of the particle position and shape, respectively. All these terms are unknown at this point and will be discussed and specified in the following sections. It is to be noted that Θ, ̇ , and ̇ have both reversible and irreversible contributions. It is also pointed out that, in deriving the above-equations and in the following derivations, the boundary terms in the space of microstructure states are assumed to vanish. In addition, there are no macroscopic fluxes through the boundaries of the system, as it is a closed system.

Reversible dynamics
The dynamics of each variable in x contains both reversible and irreversible contributions. Reversible contributions arise from dynamics that are invariant to time reversal [11] . This is analogous to affine deformation in continuum mechanics. The reversible rate of change of the position and shape of each particle can be obtained in the case of affine deformation under an imposed flow field v .
Under affine deformation, the reversible rates of change in the position and shape of a particle i are given by respectively, where = . The former, Eq. (8a) , is identical to the reversible dynamics of the particle position given in [11] . The reversible dynamics of T i , (8b) , is based on the discussion in the introduction ( Section 1 ). The only difference is that (8b) is written in a form that, whatever the symmetry properties of the actual state T i , the change of T i is manifestly symmetric. In other words, the model is set-up for arbitrary T i , but that, when starting with symmetric T i , the symmetry of T i will be preserved during the dynamics. This means that, the subset of symmetric tensors { T i } is invariant with respect to the dynamics. In this way, i.e. by not enforcing the symmetry of T i as part of their definition, complications when taking derivatives with respect to T i are avoided.
In order to model rigid particles, it is noted that these particles do not deform in size or shape, but they can change in orientation due to, for instance, the applied flow. To account for this case, one requires mixed (Gordon-Schowalter) derivatives instead of only upper convected (8b) , or lower convected, as described in [26][27][28] for the flow behavior of rigid axisymmetric particles.
Using condition (1b) together with (8), the reversible contribution to the temperature change can be obtained as where s is entropy density per unit volume. In deriving (9) , we have used the following relation with = and = .
In order to derive the stress tensor in terms of the particle dynamics, the condition for the conservation of energy in reversible dynamics (1a) together with (9) is used. The stress tensor is obtained as with int the internal energy density of the system. The expressions (12b) and (12c) resemble derivatives of the Helmholtz free energy; for more detail, the reader is referred to [11] .
In deriving (11) , we have used relation (10) with = int and = int . In addition, we have assumed that the consistency relation int = holds.

Irreversible dynamics
As stated in Section 3.3 , the temperature change Θ as well as the rate of change of the mesoscopic variables { Q i } and { T i } contain contributions due to irreversible dynamics. In this section, we derive these contributions making use of the thermodynamic conditions (1c) and (1d) . In fact, also effects such as heat conduction and viscous flow can be considered. However, in the present work, we neglect these effects for clarity, in order to focus on the particle dynamics. Interested readers are referred to [11] , where the effect of heat conduction and viscous flow have been discussed.
In this work, we consider the overdamped particle dynamics as the main effect for irreversible dynamics. Forces acting on the particles result directly in a change in the particle positions and shapes, since inertia effects are neglected [11] . This dissipative effect results in a temperature change, in addition to a change in the particle distribution. In other words, it results in an irreversible contribution to Θ, { ̇ } , and { ̇ } .
Using condition (1c) , the irreversible contribution to the temperature change is obtained as Using condition (1d) together with (14) , A quasi-linear solution to Eq. (15) for n > 0 and ϑ > 0 is given by where is a positive semi-definite and symmetric N × N × 3 × 3 tensor with the property The tensor is a positive semi-definite N × N × 3 × 3 × 3 × 3 tensor that has the symmetry properties Eq. (18b) ensures that irreversible changes in T i are symmetric. It is noted that the ansatz (16) implies that the thermodynamic forces and fluxes related to position and shape are mutually decoupled. This is done for simplicity rather than for necessity, and the model formulation can be extended in this direction in a straightforward way. In order to keep the model as general as possible, the tensors and remain unspecified in deriving the model. Particularly, they may depend on the mesoscopic degrees of freedom.

Static building blocks
In Section 3 , the model is developed for a general form of the system internal energy int and entropy S . The internal energy and entropy of the system are defined in this section, in order to make the model system specific.
Looking at the constituents of the overall system, it is clear that the suspending solvent as well as the N -particles contribute to the energy and entropy of the system. In particular, the free energy density of only the solvent is given in terms of the densities of internal energy and entropy of the solvent as where and depend on the density and temperature of the system. In the presence of particles, an additional contribution emerges due to the effective interaction energy Φ of the particles. This effective interaction energy has the character of a Helmholtz free energy and can hence be decomposed into energetic and entropic contributions [11,22] , Each of these contributions is, in general, density-and temperaturedependent in addition to their dependence on the microstructure. Based on (19) and (20) , the internal energy of the system consists of a contribution due to the internal energy of the solvent, , and a contribution due to the particles, Φ E . Similarly, the entropy of the system is composed of a contribution due to the entropy of the solvent, , and an entropic contribution due to the particles, Φ S . Therefore, where B is the Boltzmann constant and p 0 is a normalization constant for the particle distribution. In (21) , it is assumed that the densities of both energy and entropy depend on the fields x only locally, i.e., they do not depend on their derivatives, e.g. ∇ r x . The last term in the entropy expression represents the configurational entropy, which eventually ensures that the Boltzmann distribution is recovered at equilibrium. It is to be noted that f in (12a) denotes the free energy density of the entire system, containing three contributions: (i) the free energy density of the suspending solvent , (19) , and two contributions related to the suspended particles, namely (ii) due to the effective interaction Φ and (iii) due to the configurational entropy, i.e. the last contribution to the entropy (21b) . The expressions (21) are used in the following for giving explicit forms of the constitutive relation for the stress tensor as a function of the microstructure, and of the particle dynamics given by the Fokker-Planck equation for the evolution of p .

Constitutive expression for the stress tensor
Using the form of internal energy and entropy given in (21) , one obtains explicit expressions for F ; p and F ; , The (non-dissipative contribution to the) stress tensor (11) , in turn, becomes where the solvent pressure is defined as

Particle dynamics
By virtue of the specific forms of the internal energy and entropy (21) , a concrete realization of the evolution equation of the distribution function p ( Eq. (7d) ) can be obtained. Making use of the force-flux relation (16) for the irreversible particle-dynamics and Eq. (22b) , the dynamics of p is described by .
This Fokker-Planck equation can be expressed in the form of stochastic differential equations in the Itô interpretation [29] as with ∑ , , representative of the fluctuation-dissipation theorem [22,30] . The symmetry property of the increment of T i in time leads to Using the fluctuation-dissipation relations (27) and (28) , as well as the symmetry property (29) , it can be shown that the symmetry properties of in Eq. (17) and in Eq. (18) follow. Each of the increments in the Wiener processes, namely d w for the position dynamics and d W for the shape dynamics, should be uncorrelated in time. This implies and It is also pointed out that fluctuations in the particle position are not correlated with fluctuations in the particle shape, ⟨ , ,

Split of shape and size changes
The model developed in Section 3 and applied in Section 4 employs the shape tensor T i for the quantification of the microstructure. The shape tensor T i , in its general form, contains information about both the shape and the size of particle i . For practical reasons, it is sometimes useful to separate the effects of size and shape. Considering for instance the case of incompressible particles, these are particles that undergo shape changes while keeping their volume constant. This situation is rather common in treating dilute fluid mixtures, in which fluid particles are dispersed in another fluid matrix [18,31] .
This section aims at differentiating between changes in size and shape, making it suitable to recover limiting cases such as incompressible particles and particles that are rigid in terms of shape change. The unconstrained shape tensor T i is split into (i) a volumetric part expressed in terms of the particle volume,  = √ det , and (ii) an isochoric part that accounts only for the particle shape expressed in terms of ̂ = ∕ 3 √ det . It is noted that the actual volume of an ellipsoid is given by (4 ∕3)  , but for convenience the prefactor 4 /3 is not included in the definition of  . In this section, we show the procedure of expressing the model obtained in Sections 3 and 4 in terms of constrained shape tensors { ̂ } and particle volumes {  } . In addition, it is shown how the model presented in [11] can be recovered as a special case.

Static building blocks
Let us first consider the static building blocks, particularly derivatives of the Helmholtz free energy Φ with respect to the unconstrained tensors { T i }. These terms appear in the constitutive relation for the stress tensor (23) , as well as in the driving forces for irreversible dynamics in (25) and (26) . The partial derivatives can be re-written as where the subscript (uc) denotes the unconstrained derivative, i.e. the derivative where the constraint on ̂ is not taken into account. The projection tensor P i is defined as with The reader is referred to [25,32,33] for more detail on the issue of constrained derivatives. The relation (33) can be used, e.g., in (23) to obtain an expression for the stress tensor as a function of the constrained shape tensors { ̂ } and sizes {  } . For symmetric T i and given that the symmetry of the increments of T i is preserved ( Eqs. (26b) and (29) ), the term involving derivatives with respect to T i in the stress tensor (23) assumes the form (36) where the symmetric and deviatoric parts of a second-rank tensor are denoted by It is to be noted that, in deriving (36) , the symmetry of T i has been only used in explicit occurrences of T i , but not when calculating partial derivatives with respect to ̂ . Eq. (36) shows that the split of the shape tensor naturally decomposes the stress tensor into a stress contribution that is deviatoric and related to derivatives of Φ with respect to ̂ , and an isotropic contribution that is related to derivatives of Φ with respect to  . A similar procedure can be used to express the driving forces for the irreversible dynamics.

Particle shape and size dynamics
According to Itô's formula [29,34] for stochastic calculus, the evolution of the constrained tensor ̂ and size  can be expressed, in terms of the evolution of T i , respectively. To make the following steps more transparent, the evolution of T i , given by (26b) , is expressed in the form The evolution of the constrained shape tensors ̂ can be written in the form, Based on (39a) and (38) , the quantities ̂ and ̂ can be expressed in terms of A and B as , , where the projection tensor P i is defined by Eq. (34) . The increment of the constrained tensor ̂ should respect the conservation of particle volume, that is ( det ̂ ) = 0 . In the following, it is shown that (39b) with (40) does not lead to any volume changes. Based on Itô's formula, using (39b) and (40) , and with and with (35a) , one finds which implies that the volume of the particle remains constant over time. The above derivation highlights the importance of using proper stochastic calculus in terms of Itô's formula, specifically the secondorder derivative terms. If the second-order terms were not considered, the increments of ̂ given by (39b) with (40) would not be volume preserving.
Similarly to the relation between (38) and (39b) , the quantities  and  in the evolution equation for the volume  , are related to the quantities A and B as , , One obtains explicit forms for ( ̂ , ̂ ) and (  ,  ) by inserting explicit expressions for A and B in (40) and (44) , respectively. It is to be noted that the evolution of ̂ and  are mutually coupled through the dependence of ( ̂ , ̂ ) and (  ,  ) on both ̂ and  (see (40) and (44) ).

Limit of spherical particles
The interpretation of the evolution of the shape tensor in terms of the evolution of the constrained tensor ̂ and the volume  facilitates describing the dynamics of spherical particles. These are particles that undergo no or only insignificant shape changes. A spherical particle is described by an isotropic shape tensor = 2 , or equivalently a constrained shape tensor ̂ = and  = 3 , where R i is the instantaneous particle radius. To obtain a model for spherical particles, on the one hand, one can explicitly describe the particles with an isotropic shape tensor, by enforcing ̂ = as a mathematical constraint. On the other hand, one can eliminate (significant) shape changes via the (Helmholtz free) energy Φ, by strongly penalizing deviations from the spherical shape. In this case, the shape dynamics is governed by rapid relaxation towards equilibrium, the latter being the subset of spherical particles. Moreover, the shape changes introduced by imposed deformation are considered insignificant compared to the shape relaxation. This justifies the use of an isotropic tensor, in particular ̂ = , in the evolution equation of the particle volume (see (43b) and (44) ). Doing so, i.e., using ̂ = in (43b) with (44) , for example the reversible deformation of the particle volume is found to be  | rev =  , which leads to for the reversible dynamics of the particle radius. Eq. (45) is in fact identical to what has been used in the derivation of the model by Hütter et al. [11] , in the limit of affine deformation of the particle surface. With respect to the transition from particle volume to particle radius, the following general comment is in place. In order to recover the distribution of the mesoscopic states in terms of particle radii ( p R ) instead of volumes (  ), the transformation of the measure must be taken into account for each particle, which implies It is well known that a transformation of variables leads to a modification of the entropy (e.g., see [22] ). Since the entropy enters the stress tensor (23) and the irreversible dynamics (25) , special care must be taken when relating R -based and -based models.

Application of the shape-size model
In this section, the model split in terms of the shape and size is further concretized. In all of that, it is assumed for simplicity that the model is formulated in three spatial dimensions, which is relevant particularly when using the Cayley-Hamilton theorem. Furthermore, we restrict our attention to non-interacting ellipsoids, i.e. the (Helmholtz free) energy Φ is a sum of single-particle contributions, called i . In general, i is an arbitrary function of the invariants of T i . However, without loss of generality, it proves to be convenient to write i as a function of the first ( ̂ 1 , ) and second ( ̂ 2 , ) invariants of the constrained shape tensor ̂ , and the third invariant ( I 3, i ) of the unconstrained tensor T i , It is noted that, in order to develop a model for the unconstrained tensors { T i }, the energy (49) is considered as a function of the unconstrained tensors, rather than the corresponding isochoric and volumetric parts. Based on (49) and making use of the Cayley-Hamilton theorem, the expression (36) related to the stress tensor (23) can be written in the form where i, K denotes the derivative of i with respect to its K th argument. The first two terms on the right-hand side (r.h.s.) are deviatoric in nature, by virtue of t r ( ̂ −1 ) = ̂ 2 , , and related to the dependence of i on ̂ 1 , and ̂ 2 , . In contrast, the third contribution on the r.h.s. is isotropic, i.e. volumetric, and originates from the dependence of i on I 3, i . In the above and in the sequel, it is used that T i is symmetric, given that the initial condition is symmetric, and that this symmetry is not violated by the dynamics.
The particle mobility tensor is chosen as follows. First, it is assumed that the ellipsoids also do not interact dynamically, more precisely that Λ ij, should vanish for i ≠ j . And second, for the further specification, the approach in [17] is inspiring, where also the concept of an ellipsoidal shape tensor is employed. For the example in this paper, the following choice is made, where Λ is a scalar prefactor that does not depend on the shape tensor. For completeness, it is mentioned that the restriction to constant Λ is merely used to keep this illustrative example simple; relaxation of this restriction (e.g. to exactly match cases studied in [17] ) is straightforward and does not pose any conceptual difficulties. The mobility tensor (51a) satisfies the symmetry properties (18). By expressing the symmetric and positive T i as = ⋅ T , one can show that indeed satisfies not only the symmetry condition (29) , but also the fluctuation-dissipation theorem (28) . Based on (51a) , the divergence of the mobility tensor can be calculated. A careful calculation results in (see [35] for details) Based on the Helmholtz free energy (49) and the relaxation-tensor related terms (51) , and using the Cayley-Hamilton theorem, the quantities ̂ and ̂ given by (40) in the evolution of the constrained tensor (39b) are obtained as The explicit form of  and  given by (44) in the evolution of the particle volume (43b) are obtained as Similar to the discussion with respect to the stress-tensor related terms in (50) , the following can be noted about the dynamics of shape and volume described by (52) and (53) , respectively. The relaxation of shape contained in (52a) is driven by the change in i with respect to ̂ 1 , and ̂ 2 , . In contrast, the relaxation of volume contained in (53a) is driven by the change in i with respect to I 3, i .
The thermal fluctuations in shape, described by (51b) , make use of the multiplicative decomposition of the shape tensor, and similarly for (52b) . In order to implement that in a numerical simulation, the following procedure can be used. Omitting the particle index i to simplify notation, it can be shown that a choice for C that satisfies = ⋅ T has the following components in a Cartesian coordinate system, with 2 = 11 22 − 2 12 and ̄ 2 = 11 23 − 12 13 . It is pointed out that none of the arguments in the square-roots are negative for a positive definite T . What could happen though is that T 11 and/or J 2 approach zero (e.g. when imposing severe uniaxial elongation in the 3-direction). In order to avoid, in a practical way, the division by small numbers in (54) as much as possible, one can rotate the coordinate system in such a way that for the rotated ̃ one obeys ̃ 11 ≥ ̃ 2 ≥ ̃ 3 .

Discussion
This paper presents a dynamic two-scale model that describes the mechanics of suspensions of permeable ellipsoidal particles. Following the principles of non-equilibrium thermodynamics, the macroscopic fluid dynamics is consistently coupled with the particle dynamics. In addition to the particle positions, the particle dynamics is captured by conformation tensors that describe the particle size, shape, and orientation of each particle. The general form of the model highlights the connection between the macroscopic response, i.e. the constitutive relation for the stress tensor (23) , and the mesoscopic particle dynamics, i.e. the Fokker-Planck equation for the distribution of particle positions, shapes, and sizes (25) . The particle dynamics is equivalently expressed in the form of stochastic differential equations (26) , that are suited for Brownian dynamics simulations. The fluctuationdissipation theorem holds not only for the particle-position dynamics, but also for the particle-shape dynamics. A specific system realization is achieved by the appropriate choice of the generalized mobility tensor and the Helmholtz free energy. For example, the effective particle permeability can be captured by the former, while particle interactions and their elastic properties can be described through the latter. Specifically, different material properties for each individual particle can be implemented.
The model is conveniently expressed in a form that differentiates between changes in size and shape, respectively. This interpretation facilitates the application of constraints on the shape or the size of the particles. Modeling particles with volume-preserving shape changes and shape-preserving size-changes can be practically achieved by strongly penalizing size and shape changes, respectively, through the Helmholtz free energy. Furthermore, the significance of the formulation using the size-shape split lies in the fact that it provides a natural unification of well-known models for purely-shape dynamics, e.g. [17,18] , and the work based on purely-size dynamics [6,[11][12][13] .
The concrete realization of the split size-shape model, presented in this work for non-interacting particles, allows for direct comparison with, for instance, the Maffettone Minale (MM) model [18] . In light of the recent work of Mwasame et al. [17] , the MM model is recovered using a Helmholtz free energy function MM that depends solely on the second invariant of the constrained shape tensor, i.e. MM = MM ( ̂ 2 , ) . The particle relaxation in the MM model is driven by the particle surface tension which allows the particle to recover its equilibrium spherical shape. Going beyond only-shape dynamics, one can modify the Helmholtz free energy to accommodate size changes as well. In particular, this can be achieved, e.g., via an additive contribution to the Helmholtz free energy function,  , that depends only on the particle volume, that is = MM ( ̂ 2 , ) +  ( 3 , ) . The volume contribution to the energy could take a form that depends on the particle elastic properties in a fashion similar to the energy used in [6] .
Although non-interacting particles have been considered in the application of the size-shape model, this restriction was used merely to minimize the complexity of the derivation. The application of the general approach can be readily extended to interacting particles, if so desired. An essential step in accounting for interacting ellipsoids is the calculation of the distance between ellipsoids. Numerical algorithms have been developed for the calculation of (i) the surface-surface distance between non-overlapping ellipsoids [36] and the signed distance between overlapping ellipsoids [37] , and (ii) the distance of closest approach [38] . Alternatively, a contact function has been formulated [39] which distinguishes between overlapping, non-overlapping, and externally-tangent ellipsoids. One can in principle employ the measure of choice for accounting for the distance between ellipsoids in an interaction potential, e.g. [40][41][42][43] . Doing so, the tools presented in this paper can be employed to elucidate the effect of microstructure on the mechanical and dynamical behavior of suspensions of permeable ellipsoidal particles, by numerical studies.
In order to model complex particles that are not necessarily symmetric and/or convex, one needs to think about how to account for the effect of the imposed flow and for the particle interactions. For both of these issues, lack of symmetry and/or lack of convexity (e.g. red blood cells) lead to intricacies. In this, the choice of structural variables for describing the particles plays a crucial role. For instance, Janus particles can be described with a vector that points towards the direction of preferred properties [44] , while liquid-crystal banana-shaped structures can be approximated by v-shaped particles described by two unit vectors and a bend angle [45,46] . Applying the model presented in the current paper to non-symmetric/non-convex particles requires approximating their shape by an ellipsoid. This approximation could result in inaccurate descriptions of the deformation behavior and of the calculation of particle interactions.