A numerical study of particle migration and sedimentation in viscoelastic couette flow

In this work, a systematic investigation of the migration of sedimenting particles in a viscoelastic Couette flow is presented, using finite element 3D simulations. To this end, a novel computational approach is presented, which allows us to simulate a periodic configuration of rigid spherical particles accurately and efficiently. To study the different contributions to the particle migration, we first investigate the migration of particles sedimenting near the inner wall, without an externally-imposed Couette flow, followed by the migration of non-sedimenting particles in an externally-imposed Couette flow. Then, both flows are combined, i.e., sedimenting particles with an externally-imposed Couette flow, which was found to increase the migration velocity significantly, yielding migration velocities that are higher than the sum of the combined flows. It was also found that the trace of the conformation tensor becomes asymmetric with respect to the particle center when the particle is initially placed close to the inner cylinder. We conclude by investigating the sedimentation velocity with an imposed orthogonal shear flow. It is found that the sedimentation velocity can be both higher or lower then the Newtonian case, depending on the rheology of the suspending fluid. Specifically, a shear-thinning viscosity is shown to play an important role, which is in-line with previously-published results.


Introduction
Suspensions of particles in viscoelastic fluids are encountered in many natural, biological, and industrial processes. Examples of these suspensions are blood, paint, and pharmaceuticals. Therefore, the behavior of particles suspended in a single viscoelastic fluid is of great interest. Many practical applications also show the occurrence of cross-streamline migration of particles. Examples are separation of particles with different sizes, suspensions in microfluidic devices where migration can lead to non-homogenous particle distributions, cells in blood causing clots, which are obstructing the flow, and other biosensing applications [1]. Furthermore, Couette flow devices are often used as a rheometer, to measure the rheology of suspensions. Therefore, the migration behavior of the particles is of great interest, since migration of particles will change the rheology of the fluid.
Understanding the sedimentation behavior of particles in a viscoelastic shear flow can for example be used to optimize the hydraulic fracturing process for obtaining natural gas. In this process, solid particles in polymeric solutions are pumped into hydraulically-induced fractures under high pressure. The particles must keep the fractures open in order to increase the rate of gas recovery. In vertically-oriented fractures, sedimentation causes the particle to settle. Therefore, sedimentation should be limited in order to limit settling and place the particles far into the fracture to increase gas recovery [2]. Migration of a sedimenting particle in a channel near a vertical wall was investigated by Singh and Joseph [3]. They found that a falling particle initially placed close enough to the wall will migrate towards the wall. However, in a Newtonian liquid, the sphere moves away from the vertical wall and attains a steady position between the channel center and the wall.
When the particle was placed in a shear flow, D'Avino et al. [4] found that three different regimes can be recognized in the migration of the particle, depending on the distance of the particle from the wall, where the migration velocity is initially linear with the particle position, but increases when the particle approaches the wall and decreases again when the particle is in very close proximity to the wall.
Migration of a particle in a viscoelastic fluid in a Couette flow was numerically and experimentally investigated by D'Avino et al. [5]. It was found that elasticity promotes migration of the particle towards the outer cylinder for almost every initial position. Only if the particle was initially placed very close to the inner cylinder, it would move towards it. Furthermore, a maximum in the migration velocity was found when the particle moves closer to the outer cylinder. High Newtonian solvent viscosity and shear thinning reduce the migration rate.
The sedimentation of particles in a flow without shear has been studied extensively. A dimensionless number called the sedimentation Weissenberg number relates the polymer relaxation time scale to the flow time scale. It has been shown that at high values for Λ, the deformation of the polymer is dominated by the flow due to the movement of the particle. This can cause significant extension of the polymer around the particle, causing pronounced elastic effects in the wake of the sphere that decrease the particle mobility [6][7][8].
Many research works have also been devoted to the sedimentation of a particle in a shear flow [9]. Van den Brule et al. [10] performed experiments to gain more insight into the effects of fluid elasticity on the settling velocity of spherical particles suspended in viscoelastic fluids. Here, the settling velocity of a spherical particle was measured in three different fluids: a Newtonian fluid, a Boger fluid and a shear-thinning fluid. From this study, it could be concluded that fluid elasticity plays a significant role in reducing the settling velocity of a spherical particle in a viscoelastic fluid. Shear-thinning, however, increased the settling velocity of the particles, since now the viscosity decreases with increasing shear rate. This was also found in another study [11]. Padhy et al. [12] numerically investigated the effect of viscoelasticity on the settling velocity of a particle in a cross-shear flow. Padhy et al. [12] calculated the drag on a particle in their simulations and found that the drag in viscoelastic fluids is increased indirectly through the increase of viscous stresses. The dominant contribution to the viscous stresses becomes asymmetric when sedimentation and shear flow are combined, leading to a decrease in settling velocity. The coupling between the two flows needs to be non-linear for the increase in drag to be substantial. This non-linear coupling is only present in viscoelastic fluids.
Murch et al. [13] recently numerically investigated the growth of viscoelastic wings and how this reduces the particle mobility in a viscoelastic shear flow. They found that for high shear Weissenberg numbers and sedimentation Weissenberg numbers of O(1), the particle mobility is decreased due to the growth of viscoelastic wings around the particle.
Migration is often a problem in flows that are more complex than simple shear, or pure sedimentation. In the current paper, a first step up in complexity is offered, by combining two flows that yield migration. The finite element method is used to investigate numerically the influence of both the sedimentation and shear Weissenberg number on the migration and sedimentation of a particle in a viscoelastic fluid under orthogonal shear in a Couette flow device. The mass of the particle is varied via the downward force on the particle to create different local shear rates around the falling particle, leading to different sedimentation Weissenberg numbers. Furthermore, the orthogonal shear rate is varied to create different shear Weissenberg numbers. The influence of both Weissenberg numbers on the sedimentation and migration of the particle is investigated, using a Giesekus model.
Migration of the particle close to the inner cylinder of the Couette is studied systematically to see how the different flow types influence the migration of the particle and what happens if the two are combined.

Problem Statement
The shear flow is created in a Couette flow device, by rotating the outer cylinder with a velocity U w , while keeping the inner cylinder stationary. When studying sedimentation, a force is applied to the particle in the z-direction. To minimize the computational costs, only a part of the Couette flow device is modeled, as is shown in Figure 1. Periodicity is assumed on surfaces S 2 and S 4 . The top and bottom of the domain are assumed to be far enough away from the particle, to not feel the presence of the particle. The angle φ, together with the inner radius R i , the outer radius R o , and the height H define the size and shape of the domain. The spherical particle has radius a, and its boundary is denoted by ∂P. A cylindrical coordinate system will be used throughout this paper, with components, using the convention [r, θ, z]. The computational domain of the Couette is denoted by Ω. Note that the positive z-direction is directed in the gravitation direction. Figure 1. Problem description of a rigid, spherical particle suspended in a viscoelastic Couette flow. The Couette flow is set up by moving the outer cylinder in the tangential direction, with velocity U w , while the inner cylinder remains stationary. Simultaneously, a force acts on the particle in the z-direction, causing it to sediment. The inner cylinder boundary is denoted by S 1 , the outer cylinder boundary is denoted by S 3 . The top and bottom of the Couette are denoted by S 5 and S 6 , respectively. The periodic side boundaries are denoted by S 2 and S 4 , where S 4 is rotated over an angle φ with respect to surface S 2 .

Balance Equations
It is assumed that the fluid is incompressible, that inertia can be neglected, and that there are no external body forces acting on the fluid. This leaves the following equations for the balance of momentum and balance of mass: where u is the fluid velocity and σ is the Cauchy stress tensor, which consists of three contributions: where p is the pressure, I the unit tensor, and η s the solvent or Newtonian viscosity. Furthermore, D = (∇u + ∇u T )/2 is the rate of deformation tensor, and τ p represents the viscoelastic stress tensor, given by: where c is the conformation tensor and G the polymer modulus, which equals η p /λ. Here, η p is the polymer viscosity and λ the polymer relaxation time. The evolution of the conformation tensor c can be described by the following equation using the Giesekus model [14]: where α is the mobility parameter and () denotes the upper-convected derivative, which is defined as: where D()/Dt is the material derivative denoted by D()/Dt = ∂()/∂t + u · ∇().

Particle Motion
It is assumed that the inertia of the particle can be neglected and that a force F ext = Fe z is applied to the particle. The balance of forces and the balance of torques that act on the particle can be expressed as follows: where n p is the outwardly-directed unit normal vector on the boundary of the particle ∂P, X is the position of the center point of the particle, and σ is evaluated on the fluid side of the particle boundary ∂P. The translational velocity of the particle is denoted by U, which is an unknown in the problem, but will be such that Equation (7) is satisfied. The relation between the particle position and velocity is obtained by the following kinematic equation:

Boundary Conditions
On the particle boundary and rigid walls S 1 and S 3 , a no-slip condition is imposed. On S 1 , the velocity is set to zero, whereas on S 3 , the velocity is tangential to the curved surface. Periodicity of the surfaces S 2 and S 4 is assumed. Surfaces S 5 and S 6 are assumed to be far enough from the particle to not feel the presence of the particle. Symmetry is prescribed in surfaces S 5 and S 6 , i.e., the velocities in the rand θ-direction are left free, while the velocity in the z-direction is prescribed to be zero. The velocities on the side curves of surfaces S 5 and S 6 are periodic in the θ-and r-direction. This yields the following boundary conditions: t · e r = t · e θ = 0 on S 5 and S 6 , where ω is the particle angular velocity and the subscripts containing a number denoting the velocities on the surfaces S 1 up to S 6 , respectively. Equation (16) imposes a periodic boundary condition for the velocities on surfaces S 2 and S 4 . The traction vector on the surface with outwardly-directed normal n is denoted by t = σ · n. The rotation φ of the velocities and tractions from surfaces S 2 to S 4 is performed by the rotation tensor R. More details on the rotation tensor can be found in Appendices A and B. On surface S 6 , the solution of the conformation of c S of the problem without a particle is prescribed since this surface is an inflow surface due to the rigid body movement of the whole mesh with the particle, as will be explained in Section 3.3. This conformation c S is obtained by solving a separate 1D problem of a Couette flow without the particle and using the steady state value of c S . This means that for the conformation tensor, the following boundary condition is applied: Furthermore, a periodic boundary condition is also applied for the conformation tensor on surfaces S 2 and S 4 . This periodic boundary condition is imposed by the following equation in the same way as is done for the velocities, but now for a tensor instead of a vector: An initial condition is needed for the conformation tensor c. In order to reduce the computation time, the steady-state solution for the problem without a particle is prescribed as the initial condition for the conformation:

Arbitrary Lagrange-Euler Formulation
The domain is described with a boundary-fitted mesh that is moved in time in such a way that the mesh moves with the particle, but not necessarily with the fluid. Therefore, the governing equations have to be rewritten in the Arbitrary Lagrange-Euler (ALE) formulation [15]. The movement of the grid has to be compensated. This only applies to the evolution equation of the conformation tensor, since this is the only equation that contains a convective term. The material derivative in Equation (6) is now written as: where ∂c/∂t| ξ denotes the time derivative in a fixed grid point and u m is the mesh velocity. Details about the mesh velocity will be explained in Section 3.3.

Numerical Method
The finite element method is used to solve the problem. For the velocity and the pressure, isoparametric, tetrahedral P 2 P 1 (Taylor-Hood) elements are used, whereas for the conformation, tetrahedral P 1 elements are used. In order to solve the constitutive equation, the log-conformation representation [16] and DEVSS-G are used for stability [17].

Weak Formulation
In this section, the weak formulation of the balance equations is presented. The motion of the particle is enforced as a constraint, using the combined equation of motion approach [18]. The periodic boundary condition for the velocity is also enforced as a constraint, using Lagrange multipliers. The weak form of this problem can now be formulated as follows: find u, p, U, c, ω and λ such that: for all admissible test functions v, H, q, d, χ, V , µ. Furthermore, D v = (∇v + (∇v) T )/2, and (·, ·) denotes an inner product on domain Ω, whereas ·, · ∂P is an inner product on ∂P; τ and ν are parameters for SUPG and DEVSS-G stabilization, respectively, and s = log c. More information on log-conformation stabilization and the function g can be found in [16], whereas more information on the DEVSS-G method and the projected velocity gradient G can be found in [17]. A finite element method on moving meshes aligned with the particle boundary is used to solve the equations. Since the elements are aligned with the particle boundary, the inner product on ∂P can be calculated in the nodes. This is done using a collocation method [18]: where n coll is the number of collocation points, which corresponds to the nodes on the particle boundary. Furthermore, µ k and x k are the Lagrange multiplier and coordinate of the kth collocation point, respectively. The same approach can be used for the weak formulation of the periodic boundary condition on the velocity. Note that the periodic surfaces have to be meshed such that the meshes are also periodic.

Implementation Periodic Boundary Conditions
In order to impose periodic boundary conditions for the velocity, the velocity has to be prescribed such that the e r , e θ and e z components on S 2 and S 4 are equal. This means that u 4 is equal to u 2 , but rotated. The actual implementation is done in Cartesian coordinates. The same approach is used for the periodic boundary condition for the conformation tensor, only now, a tensor rotation is used instead of a vector rotation. More details on the implementation can be found in Appendices A and B.

Mesh Movement
Since the mesh will move in time, it is necessary to move the mesh to track the particle boundary. This is done using the ALE method. The starting mesh in the ALE method is used until it becomes too distorted. When this happens, a new mesh, covering the same domain as the old one, is generated using Gmsh [19], and the solution on the old mesh is projected onto the new one. The projection problem is solved to obtain the solution variables on the new mesh. This remeshing and solving the projection problem is done in the same way as was done by Jaensson et al. [20]. However, the projection leads to inaccuracies, and the remeshing procedure is very time consuming. Therefore, the number of remeshing steps should be reduced as much as possible during the simulations.
The mesh velocity u m can be set to an arbitrary value, and this can be exploited to reduce the remeshing steps, or even avoid remeshing [15]. Therefore, the whole mesh moves as a rigid body with the particle in the θ-and z-direction, while the mesh still deforms in the r-direction due to the movement of the particle. The mesh displacement in the r-direction due to the migration of the particle can be obtained by solving a Laplace equation, in order to guarantee a smooth mesh motion [15].

Prescribing the Conformation Tensor
Since the conformation tensor obeys a hyperbolic partial differential transport equation, the conformation has to be prescribed at inflow boundaries [21,22]. The bottom surface S 6 of the domain can be seen as an inflow surface, since the whole mesh is moved as a rigid body with the particle in the z-direction. Therefore, a solution of the conformation tensor from the same problem without a particle is prescribed at surface S 6 .

Time Integration
The time stepping approach now contains the following steps: Step 1: Compute the particle position at the beginning of every time step. The new particle position can be found using an explicit Euler scheme for Equation (9): For subsequent time steps, however, a second-order Adams-Bashforth scheme is used: Step 2: Update the mesh, using the new particle position as explained in Section 3.3.
Step 3: The mesh velocities can now be obtained by numerically differentiating the mesh coordinates. For the first time step, the mesh velocity u m is obtained in each node using a backward differencing scheme: where x are the nodal coordinates of the mesh. For subsequent time steps, a second-order backward differencing scheme is used: The zand θ-component of the mesh velocity are set equal to the particle velocity in the zand θ-direction. In this way, the whole mesh moves with the particle.
Step 4: Compute u n+1 and p n+1 . The method of D'Avino and Hulsen [23] for decoupling the momentum balance from the constitutive equation is applied, using DEVSS-G for stability [17].
Step 5: After solving for the new velocities and pressures, the actual conformation tensor c n+1 is found using a second-order, semi-implicit Gear scheme with conformation prediction for Equation (5) using SUPG and the log-conformation representation for stability. The conformation tensor is solved component wise, meaning that not all conformation tensor components are known at the same time. Because of this, it is not possible to enforce periodicity of c on surfaces S 2 and S 4 via a constraint. Therefore, an iteration is performed until the periodic solution for c is found. During this iteration, the conformation tensor has to be described on the inflow-nodes of surfaces S 2 and S 4 . See Section 3.3.1. Therefore, it is first computed which of the nodes on surface S 4 are inflow nodes by computing whether (u − u m ) · n is positive. Here, n is the normal vector of surface S 4 . Due to periodicity, outflow nodes on S 2 have to be inflow nodes on S 4 .

Results
In this section, the convergence of the migration problem will be investigated. After that, the results of the migration of a sedimenting particle without Couette flow, migration of a particle in a Couette flow without sedimentation, and migration of a particle with the two flows combined will be discussed. Finally, the results of the sedimentation of a particle in a Couette flow for two different sets of fluid parameters will be discussed. In order to increase the sedimentation Weissenberg number, the applied downward force on the particle is increased. Since a force is applied on the particle, the sedimentation velocity of the particle is not known a priori. Therefore, the sedimentation Weissenberg number is defined as follows: where U N is the Stokes velocity defined as U N = F/(6πη 0 a). The shear Weissenberg number is defined as: The ratio between the solvent viscosity and the total viscosity can be expressed as the ratio between the solvent viscosity and the polymer viscosity and can be described as β = η s /(η p + η s ). The viscosity ratio β is chosen to be 0.35, together with a non-linear parameter α of 0.2. The particle is positioned initially at a distance 3a from the inner cylinder of the Couette.
The outer radius of the outer cylinder has a value of 90-times the particle radius, while the radius of the inner cylinder has a value of 50-times the particle radius, yielding a gap size of 40-times the particle radius. The height of the Couette flow device is 80-times the particle radius. To save computational costs, only a part of the Couette flow device is modeled, where surface S 4 is rotated over an angle of θ = 45 • with respect to surface S 2 . By using θ = 45 • , we actually model a Couette flow with eight particles in it, because of the periodic boundary conditions. For this angle, the distance between the particles is larger than 40-times the particle radius, which is large enough to exclude the effect of particle-particle interactions [12].

Convergence
In this section, the convergence of the migration problem will be investigated. Convergence was the hardest to obtain for a high shear Weissenberg number combined with a high sedimentation Weissenberg number. Therefore, a convergence study was performed for a shear Weissenberg number of Wi = 2 and a sedimentation Weissenberg number of Λ = 2. In order to choose a mesh with elements small enough that convergence is reached, the migration velocity of a particle, made dimensionless by dividing by a/λ, was plotted for different element sizes h p around the particle. For these simulations, meshes were constructed with six elements across the gap of the Couette flow device and h p = 20, h p = 30, and h p = 40 elements on the equator of the particle. This element size was extended over a distance a from the particle and then slowly increased to the element size of the surrounding mesh over a distance 9a. The mesh was refined near the inner radius of the Couette, where 15 elements were placed along the inner radius. The mesh for h p = 30 is shown in Figure 2. The results for mesh convergence are shown in Figure 3a. Time convergence was studied by plotting again the dimensionless migration velocity of the particle, for a mesh with h p = 30, but now with different time step sizes ∆t. The result is shown in Figure 3b. Note that we use quadratic elements for the velocity, meaning that h p = 30 means 61 nodes on the equator of the particle, and six elements across the gap means 13 nodes across the gap.    Figure 3. Migration velocity of a particle placed at a distance 3a from the inner cylinder of a Couette for different element sizes around the particle h p and different time step sizes ∆t.
These figures show that for small enough element sizes and time step sizes, the lines superimpose. The convergence is studied here for the most complicated scenario of high Weissenberg numbers and the particle placed originally close to the inner cylinder. Since Figure 3 shows relatively good convergence, we expected that convergence would be better for the less difficult scenarios. In order to have a good balance between accuracy and computational costs, a element size of h p = 30 and a time step size of ∆t = 0.05 were chosen for the simulations.

Migration
In this section, it is investigated how the two complex flows (shear and sedimentation) influence the migration of a particle placed close to the inner cylinder of a Couette flow device. In order to do this, first the effect of the two separate complex flows is investigated. This is done for different initial particle positions in the Couette flow device, close to the inner cylinder. The dimensionless particle radial position is defined as: where x p is the radial position of the center of the particle at time t/λ = 10 and a is the particle radius. Therefore, Γ indicates how many times the particle radius is away from the inner cylinder.
The migration velocity at time t/λ = 10 is used in the figures in this Results Section. For a time that is ten-times the relaxation time of the fluid, it was assumed that the viscoelastic stresses were fully developed, and therefore, all transient effects were due to the migration of the particle. For the representation of the results, the cross-section indicated in red in Figure 2 was used. From here on, this cross-section is referred to as S CS .

Migration for a Falling Particle without Shear in Couette Flow
The migration of a falling particle close to the inner cylinder of the Couette was first investigated without applying a shear flow. This means that the velocity of the outer cylinder U w = 0. The result is shown in Figure 4. Here, the migration velocity U r of the particle is made dimensionless with a/λ.  This figure shows that the movement of the particle towards the inner cylinder slows down significantly when the particle approaches the inner cylinder. For small values of Λ, it even seems to become positive, i.e., the particle is migrating away from the inner cylinder. This effect of the wall is much more severe for the case of a migrating particle in a Couette flow without sedimentation than for a migrating particle without shear, but with sedimentation, as we will see in Section 4.2.2. The figure also shows that the migration velocity goes to zero when the particle position is further away from the inner cylinder. The movement of the particle towards the inner cylinder speeds up for increasing Λ (i.e., higher negative U r ). The trace of the conformation tensor around the particle for initial particle positions Γ i = 2 and Γ i = 5 and Λ = 1 without shear is shown using S CS in Figure 5 at time t/λ = 10. This figure shows that the trace of the conformation tensor around the particle becomes slightly asymmetric and larger at the side of the particle close to the cylinder, when the particle is positioned close to the inner cylinder.

Migration of a Particle in Couette Flow without Sedimentation
The migration of the particle is now investigated for the case where there is no sedimentation, but the shear rate is varied. This means that the downward force on the particle is set to zero and the velocity of the outer cylinder U w is varied. The result is shown in Figure 6.
This figure shows that the movement of the particle towards the inner cylinder speeds up with increasing shear Weissenberg number. For low values of Wi, the migration velocity is nearly zero, independent of the position of the particle. For particles close to the inner cylinder, the migration velocity is negative, indicating that the particle migrates towards the inner cylinder. For larger distances from the inner cylinder, the migration velocity is positive for high values of Wi, i.e., the particle seems to migrate away from the inner cylinder. Figure 7 shows the trace of the conformation tensor around the migrating particle with Wi = 1 and no sedimentation, for an initial position Γ i = 2 and Γ i = 5 at time t/λ = 10 using S CS . This figure shows that the presence of the inner cylinder causes the trace of the conformation to by asymmetric around the particle when the particle is placed close to it. The conformation is higher at the side of the particle close to the inner cylinder, indicating that the polymers are more stretched in the small gap between the particle and the inner cylinder. (a) tr(c) for Γ i = 2 (b) tr(c) for Γ i = 5 Figure 7. The trace of the conformation tensor around the particle for different initial positions from the inner cylinder at time t/λ = 10 for Λ = 0 and Wi = 1, using S CS .

Migration of a Falling Particle in a Shear Flow
Now, the two complex flows are combined by first fixing the shear Weissenberg number (Wi = 2) and varying the sedimentation Weissenberg number and vice versa. The results are shown in Figure 8.
Both figures show that when the particle approaches the inner cylinder Γ = 2, this slows down the movement of the particle towards it for high shear and sedimentation Weissenberg numbers. Here, the particle nearly touches the inner cylinder at t/λ = 10, drastically reducing the mobility of the particle in the r-direction. Furthermore, it can be observed that combining the two complex flows (solid lines) leads to significantly larger migration velocities than the superposition of the migration velocities of the two separate flows (dashed lines). This effect was more pronounced for large values of both Weissenberg numbers, indicating that for large values of Λ and high Wi, the coupling of the flows becomes non-linear, and therefore, the migration velocity for the combination of the two flows will no longer be a simple superposition of the migration velocities of the two separate flows.  Figure 9 shows the trace of the conformation tensor at time t/λ = 10 for initial particle positions Γ i = 2 and Γ i = 5 for Λ = 1 and Wi = 1, using S CS . This figure again shows that the trace of the conformation tensor is larger at the side of the particle close to the cylinder, when the particle approaches it. For Γ = 5 the trace of the conformation tensor is symmetric and high at both sides of the particle due to the shear flow.
(a) tr(c) for Γ i = 2 (b) tr(c) for Γ i = 5 Figure 9. The trace of the conformation tensor around the particle for different initial positions from the inner cylinder at time t/λ = 10 for Λ = 1 and Wi = 1, using S CS .
These contour plots show, that it is possible to influence the migration velocity of a particle by combining shear and sedimentation flows and changing the shear and sedimentation Weissenberg numbers. Figure 10b for example shows that the migration velocity can be accelerated, when both Weissenberg numbers are high. The migration can be slowed down by lowering one, or both of these Weissenberg numbers. Markers x indicate the data points; between these points, the data are interpolated using the MATLAB R function contourf.

Sedimentation
To investigate the influence of shear and sedimentation Weissenberg number on the sedimentation of a particle, the two complex flows are combined. The particle is now placed initially in the middle of the Couette flow device. Simulations are performed for two different sets of material parameters shown in Table 1. Here, Fluid 1 is used to mimic a Boger fluid with a large value for β and small value for α [24], so the effect of shear thinning can be excluded from the results. Fluid 2 has a much higher value for α, meaning that the fluid is more shear thinning, and a slightly lower value for β. These parameters are the same as for the migration simulations. The zero-shear viscosity η 0 = η s + η p is the same for both fluids. The sedimentation velocity for both fluids, scaled with the sedimentation velocity in a Newtonian fluid with viscosity η 0 = η + η p , is plotted as a function of the shear Weissenberg number for different sedimentation Weissenberg numbers in Figure 11. This figure shows that the sedimentation velocity is indeed decreasing with increasing shear Weissenberg number. This effect is much more pronounced in Fluid 1 than in Fluid 2. Furthermore, the sedimentation velocity is higher in Fluid 2 than in Fluid 1, and for small Wi, the sedimentation velocity is almost equal to the sedimentation velocity of a particle in a Newtonian fluid in Fluid 1, whereas it is higher than the Newtonian velocity in Fluid 2. This confirms that shear thinning is increasing the sedimentation velocity of a falling particle, also decreasing the effect of elasticity that causes the sedimentation velocity to decrease for higher Wi, as was also found in the literature [10,11]. The contour plots of the sedimentation velocity for both fluids are shown in Figure 12. The contour plots also clearly show different behavior for the two fluids. For Fluid 1, the sedimentation velocity clearly decreases with increasing Wi, while for Fluid 2, the sedimentation velocity seems almost independent of Wi for small Λ, whereas the effect of decreasing velocity for increasing Wi for higher Λ is much less pronounced as for Fluid 1. For Fluid 2, the effect of increasing Λ is also bigger than for Fluid 1.

Conclusions
A numerical model was developed in order to investigate numerically the influence of shear and sedimentation Weissenberg number on the migration and sedimentation of a particle in a Couette flow device. DEVSS-G and the log-conformation representation are used for numerical stability. To reduce computational costs, only a part of the Couette flow device is modeled. Periodic boundary conditions are applied at the symmetry planes, using vector and tensor rotations. We believe that this domain excludes the effect of particle-particle interactions. In order to reduce the amount of remeshing steps, the mesh moves as a rigid body with the particle in the θ-and z-direction, while it still deforms due to the movement of the particle in the r-direction. This is done using the ALE method.
The particle was initially placed at different positions close to the inner cylinder of the device. First, the influence of the two separate complex flows was investigated by either applying no velocity U w on the outer cylinder, and therefore Wi = 0, or applying no downward force on the particle, and therefore Λ = 0. The suspending fluid is a Giesekus fluid.
First, the migration velocity of a falling particle without applying a shear flow was studied. It was found that the migration velocity was almost zero when the particle was far away from the inner cylinder. When the particle approached the inner cylinder, the particle migrated towards it, and this migration velocity was increased by increasing Λ. At a certain position, the particle got so close to the inner cylinder that the mobility of the particle in the r-direction was drastically decreased, decreasing the migration velocity. Secondly, the migration of a particle in a Couette flow, without sedimentation, was studied. Here, it was found that for small values of Wi, there was no migration. When Wi was increased, the particle started to migrate towards the inner cylinder when it was initially placed close to it. When the particle was initially positioned further away from the cylinder, the particle migrated away from it. When the two flows were combined, it was found that combining the two complex flows leads to a higher migration velocity than the sum of the migration velocities of the two separate flows. This effect was more pronounced for higher Weissenberg numbers, indicating that for high values of Λ and Wi, the coupling of the two flows becomes non-linear. The trace of the conformation tensor was found to be asymmetric around the particle when the particle was in close proximity of the inner cylinder. The trace of the conformation is higher at the side of the particle close to the inner cylinder, indicating that the polymers are more stretched here.
To study the influence of shear and sedimentation Weissenberg number on the sedimentation of a particle, two different sets of material parameters were used, where Fluid 2 is much more shear thinning than Fluid 1 and fluid 1 has a large β, mimicking a Boger fluid. It was confirmed that shear thinning increases the sedimentation velocity of the particle. Furthermore, it was found that the sedimentation velocity indeed decreases for increasing Wi in Fluid 1, whereas this effect was less pronounced for Fluid 2. For small values of Λ, this effect was not seen at all for Fluid 2, indicating that shear thinning opposes the effect of elasticity. Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Periodic Boundary Condition Velocity
A schematic overview of the domain in 2D is shown in Figure A1. The domain is shown in 2D for clarity, because the velocity in the z-direction does not undergo a rotation. The z-component is now the vector pointing into the paper. To obtain the velocity components on surface S 2 equal to the velocity components on S 4 , u 2 has to be rotated over a counterclockwise angle φ, since surface S 4 is rotated over an angle φ with respect to surface S 2 . The rotation of a vector can be performed by taking the dot product with a tensor R, which rotates the vector in a constant coordinate system. After applying the rotation, the velocity on surface S 2 equals the velocity on S 4 . In a Cartesian coordinate system (e x , e y , e z ), the rotation is expressed as follows: Using Equation (A1), U x , U y , and U z on surface S 4 can be expressed as follows: U y = U x sin φ + U y cos φ, (A3) The same approach can be used to calculate the traction components on surfaces S 4 and S 2 , in order to obey the periodic boundary conditions.

Appendix B. Periodic Boundary Condition Conformation Tensor
where R is a tensor that performs a rotation in a constant coordinate system. Furthermore, the inner product of a rotated tensor A with a rotated vector x gives a rotated vector b, and the same thing holds for the unrotated tensor and vectors: If these equations are combined with Equation (A5), the following can be obtained: Equation (A8) can be multiplied with R −1 , (R −1 · R = I), to obtain: Equation (A9) combined with Equation (A7) states that A = R −1 · A · R. If A = R −1 · A · R is pre-multiplied with R and post-multiplied with R −1 , we obtain: R · A · R −1 = R · R −1 · A · R · R −1 . Invoking that R is an orthogonal tensor (R T = R −1 ), the rotation of tensor A to tensor A can be performed as follows: