An extensible lattice Boltzmann method for viscoelastic flows: complex and moving boundaries in Oldroyd-B fluids

Most biological fluids are viscoelastic, meaning that they have elastic properties in addition to the dissipative properties found in Newtonian fluids. Computational models can help us understand viscoelastic flow, but are often limited in how they deal with complex flow geometries and suspended particles. Here, we present a lattice Boltzmann solver for Oldroyd-B fluids that can handle arbitrarily-shaped fixed and moving boundary conditions, which makes it ideally suited for the simulation of confined colloidal suspensions. We validate our method using several standard rheological setups, and additionally study a single sedimenting colloid, also finding good agreement with literature. Our approach can readily be extended to constitutive equations other than Oldroyd-B. This flexibility and the handling of complex boundaries holds promise for the study of microswimmers in viscoelastic fluids.


Introduction
Recent years have seen a surge of interest in the study of viscoelastic fluids, due to increased experimental understanding and several intriguing results that were obtained in these media. In particular, microswimmers in viscoelastic fluids show a richer set of behaviors than possible in simple (Newtonian) fluids, which include: the selfpropulsion of a microswimmer with a single hinge [1,2], which is forbidden in a Newtonian fluid at low Reynolds number by Purcell's scallop theorem [3]; enhanced rotational diffusion of thermophoretic Janus swimmers, due to time-delayed translation-rotation coupling in polymer suspensions [4]; a peak in the motility of Escheria coli bacteria as a function of the polymer concentration and thus complexity of the fluid [5]; and a fundamental change in the way a microorganism propels in response to the rheology of the medium [6]. With the majority of industrially and biologically relevant fluids being viscoelastic [7,8], many more such surprises lie ahead of us.
This has motivated the development of a wide range of theoretical and numerical methods. However, solving the associated hydrodynamic problem remains an open challenge, both in terms of efficiency and in defining the relevant constitutive equations. Much of the numerical work has focused on well-established, albeit basic, models of complex media, such as polymeric fluids described by Oldroyd-B [9] and FENE-P [10,11]. Examples of such solvers applied to microfluidic problems include the finite volume method [12,13], the finite element method [14,15], multi-particle collision dynamics (MPCD) [16,17], dis-sipative particle dynamics [18], the immersed boundary method [6], smoothed-particle hydrodynamics [19], as well as explicit-polymer models based on Stokesian dynamics [20] and MPCD [21]. The open problem is how to simulate a fluid with a well-defined rheological response, while also allowing for the incorporation of colloidal particles. Lattice Boltzmann [22][23][24] (LB) methods hold particular promise to achieve this goal due to their computational efficiency [25] and facile boundary [26] and particle coupling [27][28][29][30], as has been demonstrated in Newtonian media. A wide variety of viscoelastic LB schemes have been conceived over the years [31][32][33][34][35][36][37][38][39][40][41][42][43][44]. However, despite this long history, which we will summarize in sect. 3.2, there remain multiple unresolved issues, especially with regard to boundary conditions.
In this paper, we address the issues of simulation of a viscoelastic fluid using LB with arbitrarily-shaped, moving boundaries. Our method is inspired by the Su et al. [42] algorithm for an Oldroyd-B fluid, which we re-derive as a finite volume scheme similar to that of Oliveira et al. [45]. This ensures momentum conservation and allows us to introduce a boundary coupling that makes no assumptions on the stress at the boundary. Compared to the LB schemes described in the literature, further advantages include low memory usage and the absence of unphysical diffusion terms. After summarizing the relevant theory in sect. 2 and laying out our numerical method in sect. 3, we benchmark our algorithm using several standard rheological tests: time-dependence of the planar Poiseuille flow in sect. 4.1, the instabilities in lid-driven-cavity flow in sect. 4.2, and extensional flow in the four-roll mill in sect. 4.3. Next, we examine the effect of the coupling of translation and rotation on the sedimentation of a sphere in sect. 4.4, showing that we reproduce the shear-induced speed-up. We discuss our findings and conclude with an outlook on future applications in sect. 5.

Theory
In this section, we summarize the equations underlying viscoelastic flow problems. They are commonly split into a Newtonian part and an additional constitutive equation, which describes the stress evolution. In terms of notation, bold symbols denote vectors (Z) i = Z i and bold sans-serif symbols denote tensors (Z) ij = Z ij .

Generalized Stokes equation
The micro-scale flows under consideration take place at low Reynolds numbers, so the hydrodynamics are governed by the time-independent Stokes equations, The first equation corresponds to momentum conservation, and the second equation is the incompressibility condition. d is the number of spatial dimensions, and u and σ are the fluid's flow velocity and stress at position r and time t. F ext is a force applied to the fluid. A Newtonian fluid's stress σ consists of a viscous stress ε and a pressure p: which simplifies eq. (1) to η n is the viscosity of the Newtonian fluid. The more general case of non-Newtonian fluids adds an extra stress τ to eq. (3). τ evolves according to a constitutive equation. Its effect on the flow may be absorbed into eq. (5)'s force via whereê i is the i-th unit vector. The total force F ext = F +F p is a sum of an applied force and the force resulting from viscoelastic stress.

Oldroyd-B fluids
There are many different constitutive equations that describe the wide range of complex fluids encountered in applications. These include Oldroyd-B [9], Jeffreys [46,47], Giesekus [48], FENE-P [10], FENE-CR [49], or Phan-ThienTanner [50]. For simplicity's sake and because it is widely studied, we will focus on Oldroyd-B. We will later indicate how our method can be extended to some of the above more realistic models. Oldroyd-B's τ corresponds to the conformation tensor of the constituent polymers, averaged over a small control volume [51]. It makes several simplifying assumptions about the fluid, including that it is made up of dumbbell polymers with zero equilibrium length and that these are very dilute [52], to arrive at the following constitutive equation: Here, the first term corresponds to advection, the next two terms are due to the polymers being stretched by the velocity gradient, and the final two terms represent the polymer relaxation. λ p is the relaxation time of the polymers, while η p refers to the viscosity added to the fluid by their presence. For use with the finite volume scheme in sect. 3.3, flux and source terms are identified in order to re-cast the equation as a conservation law:

Dimensionless numbers
It is common practice in fluid mechanics to introduce certain dimensionless numbers. Many phenomena do not depend on precise parameter values, but rather on the relative significance of individual physical effects. The Reynolds number gives the ratio of inertial forces to viscous forces: where L is a characteristic length scale of the flow and U a characteristic velocity. Re represents the relative importance of inertia. The Stokes eq. (1) is only valid in the limit of Re 1. The Deborah number is determined by the ratio of the elastic relaxation time to the characteristic time scale of the flow [53]: thus representing the degree of elasticity in response to a deformation. The Weissenberg number relates the elastic relaxation time to the characteristic rate at which the deformation is driven [53]: Finally, it is convenient to introduce the polymer viscosity fraction which can easily be varied while keeping the total viscosity constant.

Numerical methods
Just as the equations in sect. 2 are split into a Newtonian part and a viscoelastic constitutive equation, we employ two separate numerical methods. The former is solved via lattice Boltzmann (LB), while the latter uses the finite volume (FV) method.

Lattice Boltzmann
LB [22,24] constructs solutions to eq. (5) by instead solving the Boltzmann transport equation (BTE), which derives from the same conservation laws. The BTE describes the time evolution of f (r, v, t), which is the probability distribution function of finding a single fluid molecule with velocity v at position r and time t. LB discretizes the BTE on a cubic (square in two dimensions) lattice with grid spacing ∆x and discrete time steps ∆t. Relaxation of f toward its Maxwellian equilibrium is linearized and only a finite set of velocities c i is permitted to allow probability to be exchanged solely between neighboring cells.
The probability distribution is thus replaced by the populations f i (r, t) = f (r, c i , t), with their equilibrium values f eq i (r, t). We use the D3Q19 velocity set in three dimensions and D2Q9 for two-dimensional systems. In the general DdQq notation, d refers to the dimensionality and q to the number of velocity vectors pointing to neighbor cells -here these are the six face and twelve edge neighbors (or four edge and four corner neighbors in two dimensions). The employed two relaxation time (TRT) collision operator relaxes symmetric (+) and antisymmetric (−) linear combinations of f i separately, and only the symmetric relaxation time λ + affects the viscosity of the fluid. λ − can be tuned to improve the accuracy of boundary conditions [54].
The full LB method is given by and −i defined via c −i = −c i . The local fluid density ρ(r, t) appears explicitly because LB does not simulate a perfectly incompressible fluid. For consistency, we did verify in our simulations that the fluid does not compress appreciably. The populations f i and the macroscopic flow fields are connected via ∆ i (r, t) in eq. (16) represents the force F ext applied to the fluid. One possible expression for it is given by Guo et al. [55][56][57]: where w i is the lattice weight factor for c i , · is the scalar/ dot product, ⊗ is the tensor/dyadic product, and Tr is the trace of a tensor. Velocity boundary conditions can be imposed on the fluid by using where r b is a boundary node with velocity u b and r b + c i ∆t is a fluid node. For no-slip conditions u b = 0, this scheme corresponds to a bounce-back of the population.

Background on viscoelastic LB
As early as 1997, Giraud et al. [31,32] used LB to compute the response of the Jeffreys fluid. This was followed up by Ispolatov and Grant [33], who employed LB to solve a linear Maxwell model, by implementing the elasticstress contribution as a body force onto their fluid. Similar approaches were followed by Li and Fang [34] and Frantziskonis [35]. Later, Frank and Li [36,37] went beyond the body-force coupling and introduced the effect of elastic stress directly into the second moment of the equilibrium distribution, which has recently been revisited by Dellar [38]. Other coupling forms were considered by Onishi et al. [39] and Osmanlic and Körner [40], who employ a Fokker-Planck-like evolution of microscopic dumbbells in a viscous fluid. This type of system is theoretically known to result in a viscoelastic response that resembles Oldroyd-B [58]. More direct approaches to reproducing Oldroyd-B were followed by Karra [41] and Su et al. [42], who solved the stress evolution equation for the corresponding constitutive relation directly using the LB fluid velocity as input to a finite difference scheme. Malaspinas et al. [43] and Su et al. [44] similarly used an LB scheme as a generic differential equation solver and treated the viscoelastic stress tensor component-wise, for both the Oldroyd-B and FENE-P constitutive relations. The LB schemes listed above are not applied to problems with boundaries [31,36], do not require explicit treatment of the stress [33][34][35]37], or use bounce-back rules to impose specific boundary conditions on the stress [32,38]. Some extrapolate stress onto boundaries to allow for cases where no analytic expression exists [43,44], while others can only be applied to systems for which the stress at the boundary is known beforehand [42]. In the following, we build upon this body of knowledge and introduce a general method capable of handling complex and moving boundaries. By doing so, we overcome the limitations of previous viscoelastic LB algorithms.

Finite volume method
FV schemes [59] are well suited for solving problems governed by conservation laws such as eq. (8). Unlike finite difference (FD) schemes, whose first-order versions are identical to first-order FV in the absence of boundary conditions, FV schemes guarantee the conservation of conserved quantities to machine precision. Su et al. [42] originally suggested coupling an FD viscoelastic solver to a Newtonian LB. We initially implemented this FD scheme before developing the method outlined below, but the violation of momentum conservation made moving boundary simulations as presented in sect. 4.4 impossible. Our method is inspired by an LB-coupled FV solver for the electrokinetic equations [60] and has similarities to other FV Oldroyd-B solvers [45].

Discretization
Equation (8) is averaged over one cell's volume V = ∆x d with surface unit normaln to become where Gauß's divergence theorem has been applied and the overbar indicates the volume average. By locatingτ andS at the cell center and J between two cells, the discrete form of this equation is obtained as where we have used the same grid spacing and time step as in sect. 3.1. The neighbor set {c i } does not necessarily need to match the one used in sect. 3.1: we have found D3Q27/D2Q9 to deliver no appreciable advantage over D3Q7/D2Q5 [61] and have thus selected the latter for its lower computational cost. We numerically interpolate Z ∈ {u, τ } as and insert these expressions into eq. (9) to obtain where the projection onto c i and the prefactor account for the case of q > 2d + 1 [60]. We will replace eq. (30) in sect. 3.3.2 with a different expression to improve numerical stability. We need to numerically differentiate u and average over the volume of one cell: Making the central-point approximation and inserting eq. (29) yields the first-order FV discretization which is identical to the corresponding FD scheme. Inserting it into eq. (10) then yieldsS(r, t).
The force eq. (6) is similarly discretized by averaging over the volume of a cell: where τ ij (r + 1 2 c i ∆t) can be obtained via eq. (29). Boundaries across which no stress is transported can be imposed on the FV scheme by using where r b is a boundary node and r b + c i ∆t is a fluid node. τ (r b ) needs to be extrapolated so that the force can continue to be obtained via eq. (36). We found constant extrapolation to be sufficient, but linear or quadratic extrapolation could be employed as needed.

Stability improvements
FV and FD schemes are known to exhibit numerical instabilities in certain situations, which result in spatial oscillations or "wiggles" [61]. This is a particularly prominent problem in the context of Oldroyd-B as the model's Péclet number [59], which relates advective transport to diffusive transport, is infinite due to the absence of a diffusive term in eq. (7). We observed stress wiggles when performing the simulations of sects. 4.3 and 4.4 as described in sect. 3.3.1. Solutions proposed for Oldroyd-B include: using higher-order differentiation schemes [42], inserting an artificial diffusion term [43], or storing u and τ on two separate grids shifted relative to each other by half a cell [45,62]. These methods increase computational cost, modify the physics of the system, and make the implementation cumbersome, respectively, so we consider alternative techniques suggested in general FV literature. These include higher-order interpolation [59,63] and differentiation [64] schemes, as well as upwind schemes [59,60]. We resorted to the latter and chose an upwind variant called corner-transport upwind scheme suggested by refs. 60, 65 and employed in our previous work [66,67]. Upwind schemes calculate advective fluxes like eq. (30) not by interpolating quantities to the midpoint between two cells, but by using the quantity from either cell, depending on which way the flow points [59]. Reference 60's method is geometrically motivated by virtually displacing a cell at r by its velocity u(r, t)∆t and calculating the virtual cell's overlap volume with all neighboring cells. This overlap corresponds to the fraction of τ (r, t) to be transferred to the respective neighboring cell. While this in principle results in fluxes in all D3Q27/D2Q9 directions, fluxes beyond the D3Q7/D2Q5 neighbor set are O(u 2 ), making them negligible here.

Moving boundaries
One way of coupling particles to an LB fluid is by the moving boundary method. It was introduced by Ladd [29] and later enhanced by Aidun et al. [30]. This method is applicable for particles much larger than the size of a grid cell and considers the cells inside the particle as no-slip conditions in the particle-co-moving frame. This corresponds to a velocity boundary condition of which can be applied via eq. (26). r, v, and ω are the position, linear, and angular velocity of the particle. Applying the boundary condition to the fluid transfers linear and angular momentum to the particle, corresponding to a force and torque The particle trajectory is obtained by summing these forces and torques, along with any externally applied ones, and integrating numerically with a symplectic Euler integrator. As a particle moves across the lattice, the set of cells it overlaps changes. When a cell at r f is converted from fluid to solid, its fluid populations are deleted. In the reverse case, new fluid populations are created at their equilibrium value, f i (r f , t) = f eq i (r f , t) from eq. (19), whose velocity u b (r f , t) is given by eq. (39). Momentum conservation during creation and destruction of populations is ensured by applying a force to the particle that balances any momentum destroyed or created: The moving boundary method has previously been extended to FV schemes [67,68], but only in the context of ion concentrations propagating according to the electrokinetic equations. In this paper, we take a similar path to apply it to the τ of a viscoelastic medium. Refs. 67, 68 take precautions to ensure that charge is conserved. We do the same here to ensure that stress -whose diagonal elements correspond to stored energy -is not created or destroyed while cells are converted between fluid and solid. Refs. 67, 68 further calculate the fraction of a cell that is overlapped by the particle and use that information to smooth out the conversion process, which they reported to significantly decrease oscillations in the particle's speed. For the simulations in sect. 4.4, we found such smoothing to be unnecessary.
A fluid cell at r f that is destroyed in front of the particle has its stress distributed among the surrounding N f fluid cells as A cell behind the particle that is created with new fluid receives and the corresponding amount is removed from the neighboring cells:

Implementation and extensibility
The methods described above are implemented using the waLBerla C++ framework [25,69]. It allows for efficient and highly parallelized implementation of local algorithms on regular grids and provides several LB implementations and a rigid-body dynamics module. The Python module pystencils [70] can be used to automatically generate code for grid-based algorithms, either for use in Python or for waLBerla. We have extended it with a generator for finite volume discretizations that automatically derives the expressions in sect. 3.3 when provided with the Oldroyd-B eqs. (8) to (10). By instead supplying, for example, the FENE-P constitutive equation [10], we could simulate that model without writing any additional code.
There are several other fluid dynamics software packages that allow the user to provide such equations and automatically derive discretizations for them, e.g. Dedalus [71] or OpenFOAM [72]. The combination of pystencils and waLBerla, however, is unique in that it allows for arbitrarily-shaped boundary conditions that change over time, which can be put to use for the moving boundaries of sect. 3.4. We forgo waLBerla for the two-dimensional simulations, since they do not require rigid-body dynamics or parallelization, and run these simulations completely from Python. In this case, LB is provided by the lbmpy module [73].

Validation and results
In this section, we solve multiple rheological benchmark systems to verify the correctness of our algorithm and implementation by comparing against results from literature. We then simulate a system involving moving boundary conditions and translation-rotation coupling in order to demonstrate the strength of the method.

Time-dependent Poiseuille flow
The planar Poiseuille geometry consists of an infinitely long channel of width L, through which flow is driven by a homogeneous force along the channel, F = F xêx . The channel walls impose a no-slip condition u((x, 0) , t) = u((x, L) , t) = 0, while the infinite length can be achieved via periodic boundary conditions in y-direction. This setup is illustrated in fig. 1 and results in a parabolic steadystate flow profile. Starting this flow in a resting Newtonian fluid causes the steady-state flow to be approached in a monotonous fashion. In a viscoelastic medium, however, the flow velocity can overshoot its steady-state value and then decay to it on a time scale of λ p . Reference 74 provides an analytical expression for the time-dependent  velocity at the center of the channel, u((x, L/2) , t), in a liquid B' model. This model has been shown to be equivalent to Oldroyd-B [75,76].
We choose the channel width L = 28∆x, the applied force F x = 10 −5 ρ∆x 4 /∆t 2 , Newtonian viscosity η n = ρ∆x 2 /∆t − η p , polymer viscosity ratios β ∈ {0.1, 0.3, 0.5, 0.7, 0.9}, and polymer relaxation times λ p /∆t ∈ {1000, 3000, 5000, 7000, 9000} for our test simulations. This corresponds to a Reynolds number of which is well within the low-Reynolds regime we are interested in.  The agreement with the analytics can be improved to around 1% in all cases if L is used as a fit parameter. This is justified by the fact that the boundary position in LB is not guaranteed to be exactly at the edge of the cell [77] and that the extrapolation of eq. (38) introduces an error for the FV method. The resulting L differs from the input parameter by ±0.6 cells, or ±0.3 per boundary, well within the range expected for regular LB.

Lid-driven cavity
The lid-driven cavity consists of a square flow cell of edge length L, with no-slip boundaries on three sides and a constant velocity boundary u((x, L) , t) = (u 0 , 0) on the top side. This is depicted in fig. 4, which also illustrates the shape of the resulting flow: a primary vortex develops near the top center of the flow cell and secondary vortices arise in the lower corners. The exact position of the center of the primary vortex, as well as the position y of the minimum of u x ((L/2, y) , ∞) and the positions x of the minimum and maximum of u y ((x, L/2) , ∞) vary with the flow parameters and have been extensively studied in literature [78][79][80][81][82], making them well-suited for comparison in the following.  We choose width and height L = 194 for the flow cell, Newtonian viscosity η n = ρ∆x 2 /∆t − η p , applied velocity u 0 = 10 −4 ∆x/∆t, polymer viscosity ratio β = 0.5, and polymer relaxation times λ p such that Weissenberg numbers Wi ∈ [0, 1] are obtained. For Wi = 0, β → 0 is also used. The Weissenberg and Deborah numbers coincide as [79] for the system under consideration. The Reynolds number is given by again placing us in the low-Reynolds regime.
For numerical reasons, the velocity boundary condition is not applied as given above. Instead, a regularization is used to remove the infinite flow divergence in the top corners. A common choice is This regularization leaves the qualitative flow features untouched, but thwarts quantitative comparison with the unregularized simulations of ref. 78. The same regularization is employed by refs. 79-82 and shall be used in the comparison below. Figure 5 shows the positions of the primary vortex and the flow velocity extrema in our simulations. Error bars correspond to the size of a cell plus the potential deviation of the true boundary position from the prescribed boundary position. One can see that the general trend from refs. 79-82 is recovered semi-quantitatively, with the exception of the nonlinear deviation of the x-component of the vortex center. Results vary significantly between these references, so that a quantitative comparison is not drawn. However, in view of this, the result in fig. 5 gives confidence in our method's accuracy. The speed with which our results were obtained, as well as the ability to refine these significantly, provide opportunities for future benchmarking.
The flow velocity at the points of interest is shown in fig. 6. Values differ between refs. 79-82 by factors of up to 2, so we only plot the comparison to ref. 82. This reference has matching flow velocities at Wi → 0 and exhibits the same trend of decreasing velocity magnitudes as our results. The vortex is observed to move toward the top left as Wi is increased. The minimum of u x moves down slightly, while both the minimum and the maximum of u y move toward the left. The deviations from the results in literature are expected as the system is very sensitive to resolution, especially at larger Wi. Our resolution was chosen such that the results had sufficiently converged.
We also performed one simulation at β → 0, the Newtonian case, and observe that this yields a different velocity than Wi → 0 at constant β = 0.5. The velocity obtained in the former way agrees with that reported by ref. 79 to within 1%. The latter way corresponds to the case of instantaneous polymer relaxation, but not vanishing viscoelasticity.

Four-roll mill
The four-roll mill consists of a square cell with length L and periodic boundary conditions. A force field of is applied to it, resulting in four counter-rotating rolls as illustrated in fig. 7. Reference 83 provides an analytical prediction for the steady-state stress in the vicinity of the We choose cell size L = 214∆x/ √ 2, Newtonian viscosity η n = 1.5ρ∆x 2 /∆t, polymer viscosity ratio β = 1 3 , maximum velocity u 0 = 10 −3 ∆x/∆t, and polymer relaxation times λ p /∆t ∈ [1000, 24000]. The simulation is run until sufficiently converged, which we find to be the case at t = 20λ p . The Weissenberg number is given by [83,84] and the Reynolds number is low at We found that our simulations lead to a decoupling of the stress at the center point from the rest of the domain due to the upwind scheme from sect. 3.3.2. To avoid this, we rotated the lattice by 45 • relative to the system as indicated in fig. 7, while ensuring that the periodic continuation of the system remains intact. We would like to stress that this is a rather unusual situation, which only appears here due to the high level of symmetry and the divergence at the central point. Such behavior will not commonly appear in soft matter systems, but when it does, it is easily identified in the stress profiles. This gives users a means to eliminate potentially problematic simulation runs.  by less than 1% between the two fits, yet the latter fit is significantly better. This is because fitting an exponent is very sensitive to small deviations. For Wi eff < 1/4, the structure of the stress profile is not captured well by the fit. This is due to the lack of a singularity, as eq. (53) was constructed with a singularity in mind [83]. Beyond this value, three regimes of solutions are recovered: continuous and differentiable at the center (Wi eff < 1/3), continuous but not differentiable at the center (1/3 ≤ Wi eff < 1/2), and diverging at the center (Wi eff > 1/2). We reproduce the expected regimes, albeit with the caveat that divergences in our scheme are not present, due to the smoothing of solutions that its discretization imposes.   Fig. 9. Geometry of the sedimenting sphere system. A sphere of radius R sediments under velocity v due to an applied force F in a periodic cubic box of length L. A torque M is applied to the sphere to rotate it with velocity ω.

Settling sphere
So far, all systems investigated were two-dimensional and had constant boundary conditions. To demonstrate our algorithm's capabilities beyond this, we simulate the sedimentation of a rotating sphere. A sphere of radius R is placed in a cubic box of size L 3 with periodic boundary conditions. A constant force F = F zêz is applied to the sphere and the counterforce −F is distributed evenly among all fluid cells so that the net momentum of the system remains zero. Furthermore, a constant torque M = M zêz is applied to the sphere to rotate it around the zaxis; a counter-torque on the fluid is not needed [85]. The geometry is illustrated in fig. 9. We choose our parameters as R = 8∆x, L/R ∈ [7.5, 30], F z = 0.008ρ∆x 4 /∆t 2 , η n = 1 6 ρ∆x 2 /∆t, η p /η n ∈ {0, 1 2 , 1, 2} and λ p = 6000∆t. The simulation is run until the velocity v of the sphere has converged, for which t = 10λ p tends to suffice. We can assume M z = 0 since it does not change the order of magnitude of the sedimentation velocity v [86] and employ Stokes' law, in order to estimate the Reynolds number for our parameter range as which lies well in the low-Reynolds regime. The Weissenberg and Deborah numbers of the system are given by [86] Wi = λ p ω z where ω = ω zêz is the measured angular velocity of the sphere. v 0 is the sedimentation velocity measured for ω z = 0, with all other parameters kept equal. ω z can be varied by changing the applied torque M z . M z is chosen such that we cover a range of Weissenberg numbers while staying below a certain value of the tangential velocity v t = ω z R in order to not jeopardize the LB's stability. To achieve this, we define a maximum surface Reynolds number which can be used to obtain a maximum allowed Weissenberg number as The parameters provided above correspond to four sets of simulations with different polymer viscosity fractions β. Within each set, the variation of ω z or M z corresponds to a change in Wi, which makes the horizontal axis of fig. 10. To obtain the value on the vertical axis, first an exponential decay is fitted to v(t) to extrapolate to t → ∞, and then simulations at different L are used to extrapolate it to L → ∞. The fit error of these two processes is used to obtain the plot error bars. In fig. 10, we also compare to an analytical solution by Housiadas [86], who expanded v/v Stokes in terms of De for arbitrary β and χ = Wi/De. Agreement is mostly within error bars up to Wi ≈ 1. Deviations beyond that are comparable to those found by ref. 86's own comparison to numerical results from ref. 87 for similar parameters. This shows that our method reproduces the analytical solution in its range of validity, while behaving similar to other methods beyond that realm.

Summary and outlook
We have introduced a method to simulate Oldroyd-B fluids with lattice Boltzmann. It uses moving boundaries to allow for the simulation of suspended colloids. We validated our method against several rheological benchmark problems and determined it to correspond well with literature for Weissenberg and Deborah numbers and viscosity fractions between zero and one, a regime relevant for many colloidal systems. We also validated our method for a specific colloidal problem, a sphere sedimenting under an applied torque, where analytical predictions are recovered in their regime of validity. Computational effort scales linearly with the number of fluid cells, while the computational cost of adding particles is negligible compared to that of simulating the fluid. Published data on the benchmarks we considered for this work covered only a small parameter space, i.e. the few most relevant points, therefore we will make our full data set available to serve as a reliable reference for future investigations. The simulation code will also be provided to enable others to study similar systems at parameters and resolutions of their choosing. Finally, thanks to the use of automatic code generation, our model and implementation are easily extensible to other viscoelastic models.
Our viscoelastic, moving-boundary LB facilitates future study of dense colloidal suspensions in viscoelastic fluids. This might include the collective sedimentation of colloids [88], which goes beyond the single-body effects discussed in sect. 4.4. The field of self-propelled colloids is of particular interest to us. Previous reports of viscoelastic enhancement of rotational diffusivity [4], for example, have spurred interest in the community. Simulation studies [21] however could not discern whether this was an effect of viscoelasticity or merely of an inhomogeneous polymer concentration. Our method does away with the explicit consideration of polymers and might settle such questions. Besides effective propulsion models [21,89], fully-resolved propulsion models [67] might also be used, which would permit investigating complex phenomena arising from the interplay of hydrodynamics, viscoelasticity, electrostatics and phoretic interactions, such as those experimentally studied in ref. 90. Our new and extensively validated method provides a first stepping stone toward such future physical modeling.
We are grateful to Martin Bauer for help with pystencils and thank Ashreya Jayaram, Alexander Morozov, Becca Thomases, and Rudolf Weeber for useful discussions. We acknowledge the Deutsche Forschungsgemeinschaft (DFG) for funding through the SPP 1726 "Microswimmers: from single particle motion to collective behavior" (HO1108/24-2). JdG further acknowledges funding by an NWO START-UP grant (740.018.013).

Research data
The numerical code and analysis scripts used to obtain the data presented in this publication are available at https: //doi.org/10.24416/INSERTED_LATER, along with a representative subset of the data.