An extensible lattice Boltzmann method for viscoelastic flows: complex and moving boundaries in Oldroyd-B fluids
Michael Kuron, Cameron Stewart, Joost de Graaf, Christian Holm
JJournal manuscript No. (will be inserted by the editor)
An extensible lattice Boltzmann method for viscoelastic flows:complex and moving boundaries in Oldroyd-B fluids
Michael Kuron , Cameron Stewart , Joost de Graaf , and Christian Holm Institute for Computational Physics, University of Stuttgart, Allmandring 3, 70569 Stuttgart, Germany, e-mail: [email protected] Institute for Theoretical Physics, Center for Extreme Matter and Emergent Phenomena, Utrecht University, Princetonplein 5,3584 CC Utrecht, The NetherlandsNovember 24, 2020
Abstract.
Most biological fluids are viscoelastic, meaning that they have elastic properties in additionto the dissipative properties found in Newtonian fluids. Computational models can help us understandviscoelastic flow, but are often limited in how they deal with complex flow geometries and suspendedparticles. 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 confinedcolloidal suspensions. We validate our method using several standard rheological setups, and additionallystudy a single sedimenting colloid, also finding good agreement with literature. Our approach can readilybe extended to constitutive equations other than Oldroyd-B. This flexibility and the handling of complexboundaries holds promise for the study of microswimmers in viscoelastic fluids.
Recent years have seen a surge of interest in the studyof viscoelastic fluids, due to increased experimental un-derstanding and several intriguing results that were ob-tained in these media. In particular, microswimmers inviscoelastic fluids show a richer set of behaviors than pos-sible in simple (Newtonian) fluids, which include: the self-propulsion of a microswimmer with a single hinge [1, 2],which is forbidden in a Newtonian fluid at low Reynoldsnumber by Purcell’s scallop theorem [3]; enhanced rota-tional diffusion of thermophoretic Janus swimmers, dueto time-delayed translation-rotation coupling in polymersuspensions [4]; a peak in the motility of
Escheria coli bac-teria as a function of the polymer concentration and thuscomplexity of the fluid [5]; and a fundamental change inthe way a microorganism propels in response to the rheol-ogy of the medium [6]. With the majority of industriallyand biologically relevant fluids being viscoelastic [7, 8],many more such surprises lie ahead of us.This has motivated the development of a wide rangeof theoretical and numerical methods. However, solvingthe associated hydrodynamic problem remains an openchallenge, both in terms of efficiency and in defining therelevant constitutive equations. Much of the numericalwork has focused on well-established, albeit basic, mod-els of complex media, such as polymeric fluids describedby Oldroyd-B [9] and FENE-P [10, 11]. Examples of suchsolvers applied to microfluidic problems include the finitevolume method [12, 13], the finite element method [14, 15],multi-particle collision dynamics (MPCD) [16, 17], dis- sipative particle dynamics [18], the immersed boundarymethod [6], smoothed-particle hydrodynamics [19, 20], aswell as explicit-polymer models based on Stokesian dy-namics [21] and MPCD [22]. The open problem is how tosimulate a fluid with a well-defined rheological response,while also allowing for the incorporation of colloidal par-ticles. Lattice Boltzmann [23–25] (LB) methods hold par-ticular promise to achieve this goal due to their computa-tional efficiency [26] and facile boundary [27] and particlecoupling [28–31], as has been demonstrated in Newtonianmedia. A wide variety of viscoelastic LB schemes havebeen conceived over the years [32–45]. However, despitethis long history, which we will summarize in sect. 3.2,there remain multiple unresolved issues, especially withregard to boundary conditions.In this paper, we address the issues of simulation of aviscoelastic fluid using LB with arbitrarily-shaped, mov-ing boundaries. Our method is inspired by the Su et al. [43]algorithm for an Oldroyd-B fluid, which we re-derive as afinite volume scheme similar to that of Oliveira et al. [46].This ensures momentum conservation and allows us tointroduce a boundary coupling that makes no assump-tions on the stress at the boundary. Compared to the LBschemes described in the literature, further advantages in-clude low memory usage and the absence of unphysicaldiffusion terms. After summarizing the relevant theory insect. 2 and laying out our numerical method in sect. 3, webenchmark our algorithm using several standard rheolog-ical tests: time-dependence of the planar Poiseuille flow insect. 4.1, steady shear flow in sect. 4.2, the instabilities in a r X i v : . [ phy s i c s . f l u - dyn ] N ov Michael Kuron et al.: An extensible lattice Boltzmann method for viscoelastic flows lid-driven-cavity flow in sect. 4.3, and extensional flow inthe four-roll mill in sect. 4.4. Next, we examine the effectof the coupling of translation and rotation on the sedimen-tation of a sphere in sect. 4.5, showing that we reproducethe shear-induced speed-up. We discuss our findings andconclude with an outlook on future applications in sect. 5.
In this section, we summarize the equations underlyingviscoelastic flow problems. They are commonly split intoa 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-serifsymbols denote tensors ( Z ) ij = Z ij . The micro-scale flows under consideration take place atlow Reynolds numbers, so the hydrodynamics are gov-erned by the time-independent Stokes equations, d (cid:88) k =1 ∂∂r k σ ki ( r , t ) = − F ext i ( r , t ) , (1) d (cid:88) k =1 ∂∂r k u k ( r , t ) = 0 . (2)The first equation corresponds to momentum conserva-tion, and the second equation is the incompressibility con-dition. d is the number of spatial dimensions, and u and σ are the fluid’s flow velocity and stress at position r andtime t . F ext is a force applied to the fluid.A Newtonian fluid’s stress σ consists of a viscous stress ε and a pressure p : σ ij ( r , t ) = ε ij ( r , t ) − p ( r , t ) δ ij , (3) ε ij ( r , t ) = η n (cid:18) ∂∂r i u j ( r , t ) + ∂∂r j u i ( r , t ) (cid:19) , (4)which simplifies eq. (1) to η n d (cid:88) k =1 ∂ ∂r k u i ( r , t ) = ∂∂r i p ( r ) − F ext i ( r , t ) . (5) η n is the viscosity of the Newtonian fluid.The more general case of non-Newtonian fluids addsan extra stress τ to eq. (3). τ evolves according to a con-stitutive equation. Its effect on the flow may be absorbedinto eq. (5)’s force via F p ( r , t ) = d (cid:88) j =1 ˆ e j d (cid:88) i =1 ∂∂r i τ ij ( r , t ) , (6)where ˆ e 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 resultingfrom viscoelastic stress. There are many different constitutive equations that de-scribe the wide range of complex fluids encountered inapplications. These include Oldroyd-B [9], Jeffreys [47,48], Giesekus [49], FENE-P [10], FENE-CR [50], or Phan-Thien-Tanner [51]. For simplicity’s sake and because it iswidely studied, we will focus on Oldroyd-B. We will laterindicate how our method can be extended to some of theabove more realistic models. Oldroyd-B’s τ corresponds tothe conformation tensor of the constituent polymers, av-eraged over a small control volume [52]. It makes severalsimplifying assumptions about the fluid, including that itis made up of dumbbell polymers with zero equilibriumlength and that these are very dilute [53], to arrive at thefollowing constitutive equation: ∂∂t τ ij ( r , t ) = − d (cid:88) k =1 u k ( r , t ) ∂∂r k τ ij ( r , t )+ d (cid:88) k =1 τ ik ( r , t ) ∂∂r k u j ( r , t )+ d (cid:88) k =1 τ kj ( r , t ) ∂∂r k u i ( r , t )+ η p λ p (cid:18) ∂∂r i u j ( r , t ) + ∂∂r j u i ( r , t ) (cid:19) − λ p τ ij ( r , t ) . (7)Here, the first term corresponds to advection, the nexttwo terms are due to the polymers being stretched bythe velocity gradient, and the final two terms representthe polymer relaxation. λ p is the relaxation time of thepolymers, while η p refers to the viscosity added to the fluidby their presence. For use with the finite volume schemein sect. 3.3, flux and source terms are identified in orderto re-cast the equation as a conservation law: ∂∂t τ ij ( r , t ) = − ∂∂r k J ijk ( r , t ) + S ij ( r , t ) (8) J ijk ( r , t ) = u k ( r , t ) τ ij ( r , t ) (9) S ij ( r , t ) = τ ij ( r , t ) d (cid:88) k =1 ∂∂r k u k ( r , t ) (cid:124) (cid:123)(cid:122) (cid:125) =0 per eq. (2) + d (cid:88) k =1 τ ik ( r , t ) ∂∂r k u j ( r , t )+ d (cid:88) k =1 τ kj ( r , t ) ∂∂r k u i ( r , t )+ η p λ p (cid:18) ∂∂r i u j ( r , t ) + ∂∂r j u i ( r , t ) (cid:19) − λ p τ ij ( r , t ) . (10) ichael Kuron et al.: An extensible lattice Boltzmann method for viscoelastic flows 3 It is common practice in fluid mechanics to introduce cer-tain dimensionless numbers. Many phenomena do not de-pend on precise parameter values, but rather on the rel-ative significance of individual physical effects. The Rey-nolds number gives the ratio of inertial forces to viscousforces: Re = ρU Lη , (11)where L is a characteristic length scale of the flow and U a characteristic velocity. Re represents the relative im-portance of inertia. The Stokes eq. (1) is only valid in thelimit of Re (cid:28)
1. The Deborah number is determined bythe ratio of the elastic relaxation time to the characteristictime scale of the flow [54]:De = λ p UL , (12)thus representing the degree of elasticity in response toa deformation. The Weissenberg number relates the elas-tic relaxation time to the characteristic rate at which thedeformation is driven [54]:Wi = λ p ˙ γ. (13)Finally, it is convenient to introduce the polymer viscosityfraction β = η p η , (14)which can easily be varied while keeping the total viscosity η = η n + η p (15)constant. Just as the equations in sect. 2 are split into a Newtonianpart and a viscoelastic constitutive equation, we employtwo separate numerical methods. The former is solved vialattice Boltzmann (LB), while the latter uses the finitevolume (FV) method.
LB [23, 25] constructs solutions to eq. (5) from the Boltz-mann transport equation (BTE), which derives from thesame conservation laws. The BTE describes the time evo-lution of f ( r , v , t ), which is the probability distributionfunction of finding a single fluid molecule with velocity v at position r and time t . LB discretizes the BTE ona lattice, typically a square/cubic lattice in two/three di-mensions, with grid spacing ∆x and discrete time steps ∆t . Relaxation of f toward its Maxwellian equilibrium islinearized and only a finite set of velocities c i is permittedto allow probability to be exchanged solely between neigh-boring cells. The probability distribution is thus replaced by the populations f i ( r , t ) = f ( r , c i , t ), with their equi-librium values f eq i ( r , t ). We use the D3Q19 velocity set inthree dimensions and D2Q9 for two-dimensional systems.In the general D d Q q notation, d refers to the dimension-ality and q to the number of velocity vectors pointing toneighbor cells — here these are the six face and twelveedge neighbors (or four edge and four corner neighborsin two dimensions). The employed two relaxation time(TRT) collision operator relaxes symmetric (+) and anti-symmetric ( − ) linear combinations of f i separately, andonly the symmetric relaxation time λ + affects the viscos-ity of the fluid. λ − can be tuned to improve the accuracyof boundary conditions [55].The full LB method is given by f i ( r + c i ∆t, t + ∆t ) = f i ( r , t ) − λ + (cid:0) f + i ( r , t ) − f eq+ i ( r , t ) (cid:1) − λ − (cid:0) f − i ( r , t ) − f eq − i ( r , t ) (cid:1) + ∆ i ( r , t ) (16)with f ± i ( r , t ) = 12 ( f i ( r , t ) ± f − i ( r , t )) , (17) f eq ± i ( r , t ) = 12 (cid:0) f eq i ( r , t ) ± f eq − i ( r , t ) (cid:1) , (18) f eq i ( r , t ) = w i ρ ( r , t ) (cid:18) c i · u ( r , t ) (19)+ 16 ( c i · u ( r , t )) − u ( r , t ) (cid:19) ,η n = ρ ( r , t ) (cid:18) λ + − (cid:19) , (20) λ − = 316 λ + , (21)and − i defined via c − i = − c i . The local fluid density ρ ( r , t ) appears explicitly because LB does not simulatea perfectly incompressible fluid. For consistency, we didverify in our simulations that the fluid does not compressappreciably. The populations f i and the macroscopic flowfields are connected via ρ ( r , t ) = q (cid:88) i =1 f i ( r , t ) , (22) u ( r , t ) = 1 ρ ( r , t ) q (cid:88) i =1 f i ( r , t ) c i + 12 F ext ( r , t ) ∆t. (23) ∆ i ( r , t ) in eq. (16) represents the force F ext appliedto the fluid. One possible expression for it is given by Guoet al. [56–58]: ∆ i ( r , t ) = 3 w i ∆x ∆t ρ (cid:20) F ext ( r , t ) · c i + 32 Tr (cid:18) G ( c i ⊗ c i ) ∆t ∆x − G (cid:19)(cid:21) (24) Michael Kuron et al.: An extensible lattice Boltzmann method for viscoelastic flows with G = 2 − λ e (cid:0) u ( r , t ) ⊗ F ext ( r , t )+ F ext ( r , t ) ⊗ u ( r , t ) (cid:1) , (25)where w i is the lattice weight factor for c i , · is the scalar/dot product, ⊗ is the tensor/dyadic product, and Tr is thetrace of a tensor.Velocity boundary conditions can be imposed on thefluid by using f i ( r b + c i ∆t, t + ∆t ) ≡ f − i ( r b , t ) + 6 ρw i ∆t ∆x c i · u b , (26)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, thisscheme corresponds to a bounce-back of the population. As early as 1997, Giraud et al. [32, 33] used LB to com-pute the response of the Jeffreys fluid. This was followedup by Ispolatov and Grant [34], who employed LB tosolve a linear Maxwell model, by implementing the elastic-stress contribution as a body force onto their fluid. Sim-ilar approaches were followed by Li and Fang [35] andFrantziskonis [36]. Later, Frank and Li [37, 38] went be-yond the body-force coupling and introduced the effect ofelastic stress directly into the second moment of the equi-librium distribution, which has recently been revisited byDellar [39]. Other coupling forms were considered by On-ishi et al. [40] and Osmanlic and K¨orner [41], who employ aFokker-Planck-like evolution of microscopic dumbbells ina viscous fluid. This type of system is theoretically knownto result in a viscoelastic response that resembles Oldroyd-B [59]. More direct approaches to reproducing Oldroyd-Bwere followed by Karra [42] and Su et al. [43], who solvedthe stress evolution equation for the corresponding con-stitutive relation directly using the LB fluid velocity asinput to a finite difference scheme. Malaspinas et al. [44]and Su et al. [45] similarly used an LB scheme as a genericdifferential equation solver and treated the viscoelasticstress tensor component-wise, for both the Oldroyd-B andFENE-P constitutive relations. Phillips and Roberts [60]provide a more in-depth review of the cited methods forviscoelastic fluids, as well as methods for generalized New-tonian fluids; for a discussion of LB methods that dealwith viscoelastic behavior of active fluids, see the reviewby Carenza et al. [61].The LB schemes listed above are not applied to prob-lems with boundaries [32, 37], do not require explicit treat-ment of the stress [34–36, 38], or use bounce-back rules toimpose specific boundary conditions on the stress [33, 39].Some extrapolate stress onto boundaries to allow for caseswhere no analytic expression exists [44, 45], while otherscan only be applied to systems for which the stress atthe boundary is known beforehand [43]. In the following, we build upon this body of knowledge and introduce ageneral method capable of handling complex and movingboundaries. By doing so, we overcome the limitations ofprevious viscoelastic LB algorithms.
Our method is inspired by an LB-coupled finite volume(FV) solver for the electrokinetic equations [62] and hassimilarities to other FV Oldroyd-B solvers [46]. FV meth-ods [63] are well suited for solving problems governed byconservation laws such as eq. (8) since they guarantee theconservation of, e.g., momentum and energy to machineprecision. We found that the hybrid scheme of finite dif-ferences (FD) and LB originally suggested by Su et al. [43]led to violation of energy conservation in the constitutiveequation in the presence of boundaries, which translatedinto a violation of the conservation of momentum. Mov-ing boundary simulations as presented in sect. 4.5 weretherefore impossible.
Equation (8) is averaged over one cell’s volume V = ∆x d with surface unit normal ˆ n to become ∂∂t ¯ τ ij ( r , t ) = − V (cid:90) V d (cid:88) k =1 ∂∂r k J ijk ( r , t )d V + ¯ S ij ( r , t )= − V (cid:90) ∂V d (cid:88) k =1 J ijk ( r , t ) n k d S + ¯ S ij ( r , t ) , (27)where Gauß’s divergence theorem has been applied andthe overbar indicates the volume average. By locating ¯ τ and ¯ S at the cell center and J between two cells, the dis-crete form of this equation is obtained as¯ τ ij ( r , t + ∆t ) ≈ − V q (cid:88) (cid:96) =1 d (cid:88) k =1 J ijk ( r + 12 c (cid:96) ∆t, t ) c (cid:96)k + ¯ S ij ( r , t ) + ¯ τ ij ( r , t ) , (28)where we have used the same grid spacing and time stepas in sect. 3.1. The neighbor set { c i } does not necessarilyneed to match the one used in sect. 3.1: we have foundD3Q27/D2Q9 to deliver no appreciable advantage overD3Q7/D2Q5 [64] and have thus selected the latter for itslower computational cost.We numerically interpolate Z ∈ { u , τ } as Z ( r + 12 c i ∆t, t ) ≈
12 ( Z ( r , t ) + Z ( r + c i ∆t, t )) (29)and insert these expressions into eq. (9) to obtain J i ( r + 12 c i ∆t, t ) =1 | c i | A τ ( r + 12 c i ∆t, t ) (cid:18) u ( r + 12 c i ∆t, t ) · c i (cid:19) , (30) ichael Kuron et al.: An extensible lattice Boltzmann method for viscoelastic flows 5 where the projection onto c i and the prefactor A = 12 d q (cid:88) (cid:96) =1 | c (cid:96) | . (31)account for the case of q > d + 1 [62]. We will replaceeq. (30) in sect. 3.3.2 with a different expression to improvenumerical stability.We need to numerically differentiate u and averageover the volume of one cell: (cid:90) V ∂∂r i u ( r , t )d V = ∆x/ (cid:90) − ∆x/ · · · ∆x/ (cid:90) − ∆x/ (cid:18) u ( r + ∆x e i , t ) − u ( r − ∆x e i , t ) (cid:19) d r · · · d r d d r i . (32)Making the central-point approximation ∆x/ (cid:90) − ∆x/ · · · ∆x/ (cid:90) − ∆x/ u ( r ± ∆x e i , t ) d r · · · d r d d r i ≈ ∆x d − u ( r ± ∆x e i , t ) (33)and inserting eq. (29) yields the first-order FV discretiza-tion ∂∂r i u ( r , t ) ≈ ∆x ( u ( r + ˆ e i , t ) − u ( r − ˆ e i , t )) , (34)which is identical to the corresponding FD scheme. Insert-ing it into eq. (10) then yields ¯ S ( r , t ).The force eq. (6) is similarly discretized by averagingover the volume of a cell: F p j ( r , t ) = 1 V (cid:90) V d (cid:88) i =1 ∂∂r i τ ij ( r , t )d V = 1 V (cid:90) ∂V d (cid:88) i =1 τ ij ( r , t ) n i d S (35) ≈ A q (cid:88) (cid:96) =1 | c (cid:96) | d (cid:88) i =1 τ ij ( r + 12 c (cid:96) ∆t, t ) c (cid:96)i , (36)where τ ij ( r + c i ∆t ) can be obtained via eq. (29).Boundaries across which no stress is transported canbe imposed on the FV scheme by using J ( r b + 12 c i ∆t ) ≡ , (37)where r b is a boundary node and r b + c i ∆t is a fluidnode. τ ( r b ) needs to be extrapolated so that the force cancontinue to be obtained via eq. (36). We found constantextrapolation τ ( r b ) ≡ τ ( r b + c i ∆t ) (38)to be sufficient, but linear or quadratic extrapolation couldbe employed as needed. FV and FD schemes are known to exhibit numerical in-stabilities in certain situations, which result in spatial os-cillations or “wiggles” [64]. This is a particularly promi-nent problem in the context of Oldroyd-B as the model’sP´eclet number [63], which relates advective transport todiffusive transport, is infinite due to the absence of a diffu-sive term in eq. (7). We observed stress wiggles when per-forming the simulations of sects. 4.4 and 4.5 as describedin sect. 3.3.1. Solutions proposed for Oldroyd-B include:using higher-order differentiation schemes [43], insertingan artificial diffusion term [44], or storing u and τ ontwo separate grids shifted relative to each other by half acell [46, 65]. These methods increase computational cost,modify the physics of the system, and make the implemen-tation cumbersome, respectively, so we consider alterna-tive techniques suggested in general FV literature. Theseinclude higher-order interpolation [63, 66] and differentia-tion [67] schemes, as well as upwind schemes [62, 63].We resorted to the latter and chose an upwind vari-ant called corner-transport upwind scheme suggested byrefs. 62, 68 and employed in our previous work [69, 70].Upwind schemes calculate advective fluxes like eq. (30) notby interpolating quantities to the midpoint between twocells, but by using the quantity from either cell, dependingon which way the flow points [63]. Reference 62’s methodis geometrically motivated by virtually displacing a cellat r by its velocity u ( r , t ) ∆t and calculating the virtualcell’s overlap volume with all neighboring cells. This over-lap corresponds to the fraction of τ ( r , t ) to be transferredto the respective neighboring cell. While this in princi-ple results in fluxes in all D3Q27/D2Q9 directions, fluxesbeyond the D3Q7/D2Q5 neighbor set are O ( u ), makingthem negligible here. One way of coupling particles to an LB fluid is by the mov-ing boundary method. It was introduced by Ladd [30] andlater enhanced by Aidun et al. [31]. This method is appli-cable for particles much larger than the size of a grid celland considers the cells inside the particle as no-slip con-ditions in the particle-co-moving frame. This correspondsto a velocity boundary condition of u b ( r b , t ) = v ( t ) + ω ( t ) × ( r b − r ( t )) , (39)which can be applied via eq. (26). r , v , and ω are theposition, linear, and angular velocity of the particle. Ap-plying the boundary condition to the fluid transfers linearand angular momentum to the particle, corresponding toa force and torque F ( t ) = V q (cid:88) i =1 c i ( f i ( r b , t ) + f − i ( r b + c i ∆t, t )) , (40) T ( t ) = V q (cid:88) i =1 ( r b − r ) × c i ( f i ( r b , t ) + f − i ( r b + c i ∆t, t )) . (41) Michael Kuron et al.: An extensible lattice Boltzmann method for viscoelastic flows
The particle trajectory is obtained by summing these for-ces and torques, along with any externally applied ones,and integrating numerically with a symplectic Euler inte-grator.As a particle moves across the lattice, the set of cellsit overlaps changes. When a cell at r f is converted fromfluid to solid, its fluid populations are deleted. In the re-verse case, new fluid populations are created at their equi-librium value, f i ( r f , t ) = f eq i ( r f , t ) from eq. (19), whosevelocity u b ( r f , t ) is given by eq. (39). Momentum conser-vation during creation and destruction of populations isensured by applying a force to the particle that balancesany momentum destroyed or created: F ( t ) = ± V∆t ρ ( r f , t ) u ( r f , t ) . (42)The moving boundary method has previously been ex-tended to FV schemes [70, 71], but only in the context ofion concentrations propagating according to the electroki-netic equations. In this paper, we take a similar path toapply it to the τ of a viscoelastic medium. Refs. 70, 71take precautions to ensure that charge is conserved. Wedo the same here to ensure that stress — whose diagonalelements correspond to stored energy — is not created ordestroyed while cells are converted between fluid and solid.Refs. 70, 71 further calculate the fraction of a cell that isoverlapped by the particle and use that information tosmooth out the conversion process, which they reportedto significantly decrease oscillations in the particle’s speed.For the simulations in sect. 4.5, we found such smoothingto be unnecessary.A fluid cell at r f that is destroyed in front of the par-ticle has its stress distributed among the surrounding N f fluid cells as τ ( r f + c i ∆t, t + ∆t ) = τ ( r f + c i ∆t, t ) + 1 N f τ ( r f , t ) . (43)A cell behind the particle that is created with new fluidreceives τ ( r f , t + ∆t ) = 1 N f + 1 q (cid:88) i =1 τ ( r f + c i ∆t, t ) , (44)and the corresponding amount is removed from the neigh-boring cells: τ ( r f + c i ∆t, t + ∆t ) = τ ( r f + c i ∆t, t ) − τ ( r f , t + ∆t ) N f . (45) The methods described above are implemented using the waLBerla
C++ framework [26, 72]. It allows for efficientand highly parallelized implementation of local algorithmson regular grids and provides several LB implementationsand a rigid-body dynamics module. The Python module pystencils [73] can be used to automatically generate codefor grid-based algorithms, either for use in Python or for
L F y x u x yL Fig. 1.
Geometry of the planar Poiseuille flow system. Aforce F is applied to a fluid in a periodic channel of width L ,which leads to a parabolic profile across the channel for theflow velocity u x along the channel. The dashed line indicateswhere the flow velocity is measured for further analysis. waLBerla . We have extended it with a generator for fi-nite volume discretizations that automatically derives theexpressions in sect. 3.3 when provided with the Oldroyd-B eqs. (8) to (10). By instead supplying, for example, theFENE-P constitutive equation [10], we could simulate thatmodel without writing any additional code.There are several other fluid dynamics software pack-ages that allow the user to provide such equations andautomatically derive discretizations for them, e.g. Deda-lus [74] or
OpenFOAM [75]. The combination of pysten-cils and waLBerla , however, is unique in that it allows forarbitrarily-shaped boundary conditions that change overtime, which can be put to use for the moving boundaries ofsect. 3.4. We forgo waLBerla for the two-dimensional sim-ulations, since they do not require rigid-body dynamics orparallelization, and run these simulations completely fromPython. In this case, LB is provided by the lbmpy module[76].
In this section, we solve multiple rheological benchmarksystems to verify the correctness of our algorithm and im-plementation by comparing against results from literature.We then simulate a system involving moving boundaryconditions and translation-rotation coupling in order todemonstrate the strength of the method.
The planar Poiseuille geometry consists of an infinitelylong channel of width L , through which flow is driven bya homogeneous force along the channel, F = F x ˆ e x . Thechannel walls impose a no-slip condition u (( x, (cid:124) , t ) = u (( x, L ) (cid:124) , t ) = 0, while the infinite length can be achievedvia periodic boundary conditions in y -direction. This set-up is illustrated in fig. 1 and results in a parabolic steady-state flow profile. Starting this flow in a resting Newtonianfluid causes the steady-state flow to be approached in amonotonous fashion. In a viscoelastic medium, however,the flow velocity can overshoot its steady-state value andthen decay to it on a time scale of λ p . This is due to thestorage of energy in the elastic medium, which is released ichael Kuron et al.: An extensible lattice Boltzmann method for viscoelastic flows 7 -0.50.00.51.01.52.02.53.03.54.04.5 u x ( x , L / , ∞ ) / u x , N e w t o n i a n ( x , L / , ∞ ) t/∆t β = 0 . β = 0 . β = 0 . β = 0 . β = 0 . Fig. 2.
Velocity at the center of a planar Poiseuille channelover time for varying viscosity ratios β and polymer relaxationtime λ p = 3000 ∆t . Symbols are numerical calculations, andlines show the analytical prediction by ref. 77. The solid linesuse L from the input parameters, whereas the dashed linesallow it to be a free fit parameter. u x ( x , L / , ∞ ) / u x , N e w t o n i a n ( x , L / , ∞ ) t/∆tλ p = 1000 ∆tλ p = 3000 ∆tλ p = 5000 ∆t λ p = 7000 ∆tλ p = 9000 ∆t Fig. 3.
Velocity at the center of a planar Poiseuille channelover time for viscosity ratio β = 0 . λ p . Symbols are numerical calculations, and linesshow the analytical prediction by ref. 77. The solid lines use L from the input parameters, whereas the dashed lines allow itto be a free fit parameter. back into the fluid on the time scale λ p . Reference 77provides an analytical expression for the time-dependentvelocity at the center of the channel, u (( x, L/ (cid:124) , t ), in a liquid B’ model. This model has been shown to be equiv-alent to Oldroyd-B [78, 79].We choose the channel width L = 28 ∆x , the appliedforce F x = 10 − ρ∆x /∆t , Newtonian viscosity η n = ρ∆x /∆t − η p , polymer viscosity ratios β ∈ { . , . , . , . , . } , and polymer relaxation times λ p /∆t ∈ { , , , , } for our test simulations. This cor-responds to a Reynolds number ofRe = ρu x (( x, L/ (cid:124) , ∞ ) Lη = ρF x L η = 0 . , (46) u L y x u u x yL Fig. 4.
Geometry of the planar Couette flow system. A veloc-ity boundary condition of u is applied to one side of a fluid ina periodic channel of width L , which leads to a linear profileacross the channel for the flow velocity u x along the channel. which is well within the low-Reynolds regime we are in-terested in.Figure 2 shows the flow velocity u x (( x, L/ (cid:124) , t ) overtime for various polymer viscosity ratios β at constantpolymer relaxation time λ p = 3000 ∆t . One can see thatthe magnitude of the overshoot increases with β . For thelargest values of β , the flow can even decay to its finalspeed in an oscillatory fashion. Figure 3 keeps β = 0 . λ p . Here it is clear that the magni-tude of the overshoot increases with λ p , which is also thecharacteristic decay time of the overshoot.Figures 2 and 3 additionally show the analytical resultfrom ref. 77 for comparison. The agreement with the an-alytics can be improved to around 1% in all cases if L isused as a fit parameter. This is justified by the fact thatthe boundary position in LB is not guaranteed to be ex-actly at the edge of the cell [80] and that the extrapolationof eq. (38) introduces an error for the FV method. The re-sulting L differs from the input parameter by ± . ± . The planar Couette geometry is similar to that of sect. 4.1,but replaces the applied force with a velocity boundarycondition of u (( x, L ) (cid:124) , t ) = u ˆ e x on one of the planes, asillustrated in fig. 4. This relative motion leads to a linearsteady-state velocity profile across the channel. The firstnormal stress difference, N = τ xx − τ yy = 2 η p λ p u L , (47)is expected to be constant over the entire channel, asobtained by solving eq. (7) with the given velocity pro-file. We choose the channel width L = 28 ∆x , the ap-plied velocity u = 10 − ∆x/∆t , Newtonian viscosity η n = ρ∆x /∆t − η p , polymer viscosity ratios β ∈ { . , . , . , . } , and polymer relaxation times λ p /∆t ∈ [1000 , t = 10 λ p . Wefind that N agrees with the prediction to within 0 . ∼ Michael Kuron et al.: An extensible lattice Boltzmann method for viscoelastic flows u L x Ly u
Fig. 5.
Geometry of the lid-driven cavity system. A squareflow cell of size L has no-slip boundaries on three sides and aconstant-velocity boundary condition of u along the fourth.The resulting flow develops a primary vortex near the top mid-dle of the cell. Along the dashed lines, flow velocity minima andmaxima are found at the red crosses. The blue arrow indicateshow the vortex center moves as Wi is increased from 0 to 1. The lid-driven cavity consists of a square flow cell of edgelength L , with no-slip boundaries on three sides and aconstant velocity boundary u (( x, L ) (cid:124) , t ) = u ˆ e x on thetop side. This is depicted in fig. 5, which also illustratesthe shape of the resulting flow: a primary vortex developsnear the top center of the flow cell and secondary vorticesarise in the lower corners. The exact position of the cen-ter of the primary vortex, as well as the position y of theminimum of u x (( L/ , y ) (cid:124) , ∞ ) and the positions x of theminimum and maximum of u y (( x, L/ (cid:124) , ∞ ) vary with theflow parameters and have been extensively studied in lit-erature [81–85], making them well-suited for comparisonin the following.We choose width and height L = 194 for the flow cell,Newtonian viscosity η n = ρ∆x /∆t − η p , applied velocity u = 10 − ∆x/∆t , polymer viscosity ratio β = 0 .
5, andpolymer relaxation times λ p such that Weissenberg num-bers Wi ∈ [0 ,
1] are obtained. For Wi = 0, β → λ p u L (48)for the system under consideration. The Reynolds numberis given by Re = ρu Lη = 0 . , (49)again placing us in the low-Reynolds regime.For numerical reasons, the velocity boundary conditionis not applied as given above. Instead, a regularizationis used to remove the infinite flow divergence in the topcorners. A common choice is u (( x, L ) (cid:124) , t ) = 16 u (cid:16) xL (cid:17) (cid:16) − xL (cid:17) ˆ e x . (50) . . . . . . . . . . . . . x / L o r y / L Wi u x minimum u y maximum u y minimum Vortex center x Vortex center y refs. [82–85] Fig. 6.
Positions of the primary vortex and flow extremain the lid-driven cavity. Colors refer to the different points.Symbols are our results, while the hatched areas indicates therange covered by the numerical results from refs. 82–85. -0.2-0.10.00.10.2 . . . . u x / u o r u y / u Wi u x minimum u y maximum u y minimum Newtonianref. [85] Fig. 7.
Velocity of the flow extrema in the lid-driven cav-ity. Colors refer to the different points. Symbols are our re-sults, while the line refers to numerical results from ref. 85.The square symbols indicate Newtonian simulations ( β → β = 0 . This regularization leaves the qualitative flow features un-touched, but thwarts quantitative comparison with theunregularized simulations of ref. 81. The same regulariza-tion is employed by refs. 82–85 and shall be used in thecomparison below.Figure 6 shows the positions of the primary vortex andthe flow velocity extrema in our simulations. Error barscorrespond to the size of a cell plus the potential devi-ation of the true boundary position from the prescribedboundary position. One can see that the general trendfrom refs. 82–85 is recovered semi-quantitatively, with theexception of the nonlinear deviation of the x -component ofthe vortex center. Results vary significantly between thesereferences, so that a quantitative comparison is not drawn.However, in view of this, the result in fig. 6 gives confi-dence in our method’s accuracy. The speed with which our ichael Kuron et al.: An extensible lattice Boltzmann method for viscoelastic flows 9 F FF FL L y x
Fig. 8.
Geometry of the four-roll mill. Four counter-rotatingforces F are applied to a periodic square flow cell of size L .This leads to a pure extensional flow at the center of the cell.Velocity and stress will be measured along the dashed line.The dotted square indicates the actual simulation domain used,which still obeys the periodic boundary conditions. results were obtained, as well as the ability to refine thesesignificantly, provide opportunities for future benchmark-ing.The flow velocity at the points of interest is shown infig. 7. Values differ between refs. 82–85 by factors of upto 2, so we only plot the comparison to ref. 85. This refer-ence has matching flow velocities at Wi → u x moves downslightly, while both the minimum and the maximum of u y move toward the left. The deviations from the resultsin literature are expected as the system is very sensitiveto resolution, especially at larger Wi. Our resolution waschosen such that the results had sufficiently converged.We also performed one simulation at β →
0, the New-tonian case, and observe that this yields a different ve-locity than Wi → β = 0 .
5. The velocityobtained in the former way agrees with that reported byref. 82 to within 1%. The latter way corresponds to thecase of instantaneous polymer relaxation, but not vanish-ing viscoelasticity.
The four-roll mill consists of a square cell with length L and periodic boundary conditions. A force field of F ( r , t ) = 8 π η n u L sin (cid:0) πL x (cid:1) cos (cid:0) πL y (cid:1) cos (cid:0) πL x (cid:1) sin (cid:0) πL y (cid:1) (51)is applied to it, resulting in four counter-rotating rolls asillustrated in fig. 8. Reference 86 provides an analyticalprediction for the steady-state stress in the vicinity of thecentral point, where the flow is purely extensional, i.e., u (( L/ , L/ (cid:124) , t ) = α (ˆ e x − ˆ e y ).We choose cell size L = 214 √ ∆x , Newtonian viscosity η n = 1 . ρ∆x /∆t , polymer viscosity ratio β = , maxi- mum velocity u = 10 − ∆x/∆t , and polymer relaxationtimes λ p /∆t ∈ [1000 , t = 20 λ p . The Weissenberg number is given by [86, 87]Wi = 4 πλ p u L , (52)and the Reynolds number is low atRe = ρu Lη = 0 . . (53)We found that our simulations lead to a decoupling ofthe stress at the center point from the rest of the domaindue to the upwind scheme from sect. 3.3.2. To avoid this,we rotated the lattice by 45 ◦ relative to the system as in-dicated in fig. 8, while ensuring that the periodic continu-ation of the system remains intact. We would like to stressthat this is a rather unusual situation, which only appearshere due to the high level of symmetry and the divergenceat the central point. Such behavior will not commonly ap-pear in soft matter systems, but when it does, it is easilyidentified in the stress profiles. This gives users a meansto eliminate potentially problematic simulation runs.Figure 9a-d shows the stress component τ xx along avertical line through the center of the cell, which is markedwith a red cross in fig. 8. A comparison with ref. 86 isdrawn by fitting with its local solution, τ xx (( L/ , y ) (cid:124) , ∞ ) = 2 η p Wi eff λ p (1 − eff ) + C (cid:12)(cid:12)(cid:12)(cid:12) yL − (cid:12)(cid:12)(cid:12)(cid:12) − . (54)We fit via C while keeping Wi eff = λ p α constant, as wellas via both C and Wi eff . We find that Wi eff only differsby less than 1% between the two fits, yet the latter fit issignificantly better. This is because fitting an exponent isvery sensitive to small deviations. For Wi eff < /
4, thestructure of the stress profile is not captured well by thefit. This is due to the lack of a singularity, as eq. (54) wasconstructed with a singularity in mind [86]. Beyond thisvalue, three regimes of solutions are recovered: continuousand differentiable at the center (Wi eff < / / ≤ Wi eff < / eff > / eff that we obtained from the fits via Wi. Compar-ison with the corresponding plot from ref. 86 is excellentup to Wi eff ≈ . ≈ . So far, all systems investigated were two-dimensional andhad constant boundary conditions. To demonstrate ouralgorithm’s capabilities beyond this, we simulate the sedi-mentation of a rotating sphere. A sphere of radius R is − .
04 0 0 . (a) (b)(c) (d)(e) τ xx (( L / , y ) | , ∞ ) / η p W i e ff λ p ( − W i e ff ) y/L − . . − .
04 0 0 . τ xx (( L / , y ) | , ∞ ) / η p W i e ff λ p ( − W i e ff ) y/L − . . − .
04 0 0 . τ xx (( L / , y ) | , ∞ ) / η p W i e ff λ p ( − W i e ff ) y/L − . . − .
04 0 0 . τ xx (( L / , y ) | , ∞ ) / η p W i e ff λ p ( − W i e ff ) y/L − . . . . . . W i e ff Wi simulationeq. (54) via C , Wi eff eq. (54) via C ref. [86] Fig. 9. (a-d) Stress τ xx (( L/ , y ) (cid:124) , ∞ ) near the center of thefour-roll mill for different polymer relaxation times λ p . Symbolsare our results, with their connecting line coming from fittingeq. (54) with an added offset. The other lines are fits witheq. (54) via one or two parameters. (e) Wi eff plotted over theWeissenberg number Wi. Symbols are our results, while theline comes from ref. 86. placed in a cubic box of size L with periodic bound-ary conditions. A constant force F = F z ˆ e z is applied tothe sphere and the counterforce − F is distributed evenlyamong all fluid cells so that the net momentum of the sys-tem remains zero. Furthermore, a constant torque M = M z ˆ e z is applied to the sphere to rotate it around the z -axis; a counter-torque on the fluid is not needed [88]. Thegeometry is illustrated in fig. 10.We choose our parameters as R = 8 ∆x , L/R ∈ [7 . , F z = 0 . ρ∆x /∆t , η n = ρ∆x /∆t , η p /η n ∈ { , , , } and λ p = 6000 ∆t . The simulation is run until thevelocity v of the sphere has converged, for which t = 10 λ p tends to suffice. We can assume M z = 0 since it doesnot change the order of magnitude of the sedimentation FML − F R z x
Fig. 10.
Geometry of the sedimenting sphere system. A sphereof radius R sediments under velocity v due to an applied force F in a periodic cubic box of length L . A torque M is appliedto the sphere to rotate it with velocity ω . . . v z / v S t o k e s Wi β = 0 . , De = 0 . β = 0 . , De = 0 . β = 0 . , De = 0 . β = 0 . , De = 0 . Fig. 11.
Sedimentation velocity of a rotating sphere in anOldroyd-B fluid. Symbols with error bars are our results, whilethe lines are theoretical predictions from ref. 89 for the sameparameters (ending at Wi max as per eq. (60)). The inset dis-plays an enlargement of the lower left region of the large graphindicated by the dotted box. velocity v [89] and employ Stokes’ law, v Stokes = F πηR , (55)in order to estimate the Reynolds number for our param-eter range asRe = 2 ρv Stokes Rη = ρF πη ∈ [0 . , . , (56)which lies well in the low-Reynolds regime. The Weissen-berg and Deborah numbers of the system are given by[89] Wi = λ p ω z (57)De = λ p v R , (58)where ω = ω z ˆ e z is the measured angular velocity of thesphere. v is the sedimentation velocity measured for ω z =0, with all other parameters kept equal. ω z can be varied ichael Kuron et al.: An extensible lattice Boltzmann method for viscoelastic flows 11 by changing the applied torque M z . M z is chosen such thatwe cover a range of Weissenberg numbers while stayingbelow a certain value of the tangential velocity v t = ω z R in order to not jeopardize the LB’s stability. To achievethis, we define a maximum surface Reynolds numberRe t,max = 2 ρv t,max Rη ≡ . , (59)which can be used to obtain a maximum allowed Weissen-berg number asWi max = λ p ω max = λ p Re t,max η ρR . (60)The parameters provided above correspond to four setsof simulations with different polymer viscosity fractions β .Within each set, the variation of ω z or M z corresponds toa change in Wi, which makes the horizontal axis of fig. 11.To obtain the value on the vertical axis, first an exponen-tial decay is fitted to v ( t ) to extrapolate to t → ∞ , andthen simulations at different L are used to extrapolate itto L → ∞ . The fit error of these two processes is used toobtain the plot error bars. In fig. 11, we also compare toan analytical solution by Housiadas [89], who expanded v/v Stokes in terms of De for arbitrary β and χ = Wi / De.Agreement is mostly within error bars up to Wi ≈
1. De-viations beyond that are comparable to those found byref. 89’s own comparison to numerical results from ref. 90for similar parameters. This shows that our method repro-duces the analytical solution in its range of validity, whilebehaving similar to other methods beyond that realm.
We have introduced a method to simulate Oldroyd-B flu-ids with lattice Boltzmann. It uses moving boundaries toallow for the simulation of suspended colloids. We vali-dated our method against several rheological benchmarkproblems and determined it to correspond well with lit-erature for Weissenberg and Deborah numbers and vis-cosity fractions between zero and one, a regime relevantfor many colloidal systems. We also validated our methodfor a specific colloidal problem, a sphere sedimenting un-der an applied torque, where analytical predictions arerecovered in their regime of validity. Computational ef-fort scales linearly with the number of fluid cells, whilethe computational cost of adding particles is negligiblecompared to that of simulating the fluid. Published dataon the benchmarks we considered for this work coveredonly a small parameter space, i.e. the few most relevantpoints, therefore we will make our full data set availableto serve as a reliable reference for future investigations.The simulation code will also be provided to enable oth-ers to study similar systems at parameters and resolutionsof their choosing. Finally, thanks to the use of automaticcode generation, our model and implementation are eas-ily extensible to other viscoelastic models. Incorporatingthermal fluctuations [91] is also conceivable. Our viscoelastic, moving-boundary LB facilitates fu-ture study of dense colloidal suspensions in viscoelasticfluids. This might include the collective sedimentation ofcolloids [92], which goes beyond the single-body effectsdiscussed in sect. 4.5. The field of self-propelled colloidsis of particular interest to us. Previous reports of vis-coelastic enhancement of rotational diffusivity [4], for ex-ample, have spurred interest in the community. Simula-tion studies [22] however could not discern whether thiswas an effect of viscoelasticity or merely of an inhomo-geneous polymer concentration. Our method does awaywith the explicit consideration of polymers and might set-tle such questions. Besides effective propulsion models [22,93], fully-resolved propulsion models [70] might also beused, which would permit investigating complex phenom-ena arising from the interplay of hydrodynamics, viscoe-lasticity, electrostatics and phoretic interactions, such asthose experimentally studied in ref. 94. Our new and ex-tensively validated method provides a first stepping stonetoward such future physical modeling.
We are grateful to Martin Bauer for help with pystencils andthank Ashreya Jayaram, Alexander Morozov, Becca Thomases,and Rudolf Weeber for useful discussions and Fabian H¨auslfor helpful comments on the manuscript. We acknowledge theDeutsche Forschungsgemeinschaft (DFG) for funding throughthe SPP 1726 “Microswimmers: from single particle motion tocollective behavior” (HO1108/24-2) and through the EXC 2075“SimTech” (390740016). JdG further acknowledges funding byan NWO START-UP grant (740.018.013).
Author contributions
Conceptualization: MK, JdG; Calculations and analysis:MK, CS; Writing: MK, CS, JdG; Supervision: MK, JdG,CH; Funding acquisition: JdG, CH; Resources: CH.
Research data
The numerical code and analysis scripts used to obtainthe data presented in this publication are available at https://doi.org/10.24416/UU01-2AFZSW , along with arepresentative subset of the data.
References
1. T. Qiu, T.C. Lee, A.G. Mark, K.I. Morozov,R. M¨unster, O. Mierka, S. Turek, A.M. Leshansky,P. Fischer, Nat. Commun. , 5119 (2014)2. T. Normand, E. Lauga, Phys. Rev. E , 061907(2008)3. E. Purcell, Am. J. Phys. , 3 (1977)4. J.R. Gomez-Solano, A. Blokhuis, C. Bechinger, Phys.Rev. Lett. , 138301 (2016)5. V.A. Martinez, J. Schwarz-Linek, M. Reufer, L.G.Wilson, A.N. Morozov, W.C. Poon, Proc. NationalAcad. Sci. , 17771 (2014)
6. C. Li, B. Qin, A. Gopinath, P.E. Arratia,B. Thomases, R.D. Guy, J. Royal Soc. Interface , 20170289 (2017)7. R.G. Larson, The Structure and Rheology of ComplexFluids (Oxford University Press, New York, 1999),ISBN 97801951219718. S. Peker, S. Helvaci,
Solid-Liquid Two Phase Flow (Elsevier Science, 2008), ISBN 97804445223759. J. Oldroyd, Proc. Royal Soc. Lond. A , 523 (1950)10. A. Peterlin, J. Polym. Sci. Part B: Polym. Lett. , 287(1966)11. R. Bird, P. Dotson, N. Johnson, J. Non-Newton. FluidMech. , 213 (1980)12. Y. Zhang, G. Li, A.M. Ardekani, Phys. Rev. Fluids ,023101 (2018)13. S. De, J. Kuipers, E. Peters, J. Padding, Phys. Rev.Fluids , 053303 (2017)14. L. Zhu, E. Lauga, L. Brandt, Phys. Fluids , 051902(2012)15. G. Li, A.M. Ardekani, Phys. Rev. Lett. , 118001(2016)16. D. Toneian, G. Kahl, G. Gompper, R.G. Winkler, J.Chem. Phys. , 194110 (2019)17. S. Sahoo, S.P. Singh, S. Thakur, Soft Matter , 2170(2019)18. B. ten Bosch, J. Non-Newton. Fluid Mech. , 231(1999)19. A. V´azquez-Quesada, P. Espa˜nol, R.I. Tanner,M. Ellero, J. Fluid Mech. , 1070 (2019)20. J. King, S. Lind, arXiv preprint arXiv:2009.12245(2020)21. A.K. Townsend, H.J. Wilson, J. Non-Newton. FluidMech. , 136 (2018)22. K. Qi, E. Westphal, G. Gompper, R.G. Winkler, Phys.Rev. Lett. , 068001 (2020)23. G.R. McNamara, G. Zanetti, Phys. Rev. Lett. ,2332 (1988)24. F. Higuera, S. Succi, R. Benzi, Europhys. Lett. , 345(1989)25. T. Kr¨uger, H. Kusumaatmaja, A. Kuzmin, O. Shardt,G. Silva, E.M. Viggen, The Lattice BoltzmannMethod: Principles and Practice (Springer, Cham,2017), ISBN 978331944647926. M. Bauer, S. Eibl, C. Godenschwager, N. Kohl,M. Kuron, C. Rettinger, F. Schornbaum,C. Schwarzmeier, D. Th¨onnes, H. K¨ostler et al.,Comput. Math. Appl. (2020)27. Q. Zou, X. He, Phys. Fluids , 1591 (1997)28. B. D¨unweg, A.J.C. Ladd, in Advanced ComputerSimulation Approaches for Soft Matter Sciences III (Springer-Verlag Berlin, Berlin, Germany, 2009), Vol.221 of
Advances in Polymer Science , pp. 89–16629. P. Ahlrichs, B. D¨unweg, Int. J. Mod. Phys. C , 1429(1998)30. A.J.C. Ladd, J. Fluid Mech. , 285 (1994)31. C.K. Aidun, Y. Lu, E.J. Ding, J. Fluid Mech. ,287 (1998)32. L. Giraud, D. d’Humi`eres, P. Lallemand, Int. J. Mod.Phys. C , 805 (1997) 33. L. Giraud, D. d’Humi`eres, P. Lallemand, Europhys.Lett. , 625 (1998)34. I. Ispolatov, M. Grant, Phys. Rev. E , 056704(2002)35. H.B. Li, H.P. Fang, Chin. Phys. , 2087 (2004)36. G.N. Frantziskonis, Phys. Rev. E , 066703 (2011)37. X. Frank, H.Z. Li, Phys. Rev. E , 036309 (2005)38. X. Frank, H.Z. Li, Phys. Rev. E , 056307 (2006)39. P.J. Dellar, SIAM J. Sci. Comput. , A2507 (2014)40. J. Onishi, Y. Chen, H. Ohashi, Prog. Comput. FluidDyn. , 75 (2005)41. F. Osmanlic, C. K¨orner, Comput. Fluids , 190(2016)42. S. Karra, Master’s thesis, Texas A & M University(2007)43. J. Su, J. Ouyang, X. Wang, B. Yang, Phys. Rev. E , 053304 (2013)44. O. Malaspinas, N. Fi´etier, M. Deville, J. Non-Newton.Fluid Mech. , 1637 (2010)45. J. Su, J. Ouyang, X. Wang, B. Yang, W. Zhou, J.Non-Newton. Fluid Mech. , 42 (2013)46. P.J. Oliveira, F. Pinho, G. Pinto, J. Non-Newton.Fluid Mech. , 1 (1998)47. H. Jeffreys, The Earth: Its Origin, History and Physi-cal Constitution (Cambridge University Press, 1976),ISBN 978052120648848. R. Bird, R. Armstrong, O. Hassager,
Dynamics ofPolymeric Liquids, Volume 1: Fluid Mechanics (Wi-ley, 1987), ISBN 978047180245749. H. Giesekus, J. Non-Newton. Fluid Mech. , 69(1982)50. M. Chilcott, J. Rallison, J. Non-Newton. Fluid Mech. , 381 (1988)51. N. Phan-Thien, R.I. Tanner, J. Non-Newton. FluidMech. , 353 (1977)52. M. Deville, T.B. Gatski, Mathematical Modeling forComplex Fluids and Flows (Springer Science & Busi-ness Media, 2012), ISBN 978364243560753. R. Bird, C. Curtiss, R. Armstrong, O. Hassager,
Dy-namics of Polymeric Liquids, Volume 2: Kinetic The-ory (Wiley-Interscience, 1987), ISBN 978047180244054. J. Dealy, Rheol. Bull. , 14 (2010)55. I. Ginzburg, F. Verhaeghe, D. d’Humi`eres, Commun.Comput. Phys. , 427 (2008)56. Z. Guo, C. Zheng, B. Shi, Phys. Rev. E , 046308(2002)57. U.D. Schiller, Ph.D. thesis, Johannes Gutenberg-Universit¨at Mainz (2008)58. U.D. Schiller, Comput. Phys. Commun. , 2586(2014)59. R.B. Bird, J.M. Wiest, Annu. Rev. Fluid Mech. ,169 (1995)60. T.N. Phillips, G.W. Roberts, IMA J. Appl. Math. ,790 (2011)61. L.N. Carenza, G. Gonnella, A. Lamura, G. Negro,A. Tiribocchi, Eur. Phys. J. E , 81 (2019)62. F. Capuani, I. Pagonabarraga, D. Frenkel, J. Chem.Phys. , 973 (2004) ichael Kuron et al.: An extensible lattice Boltzmann method for viscoelastic flows 13
63. H. Versteeg, W. Malalasekera,
An Introduction toComputational Fluid Dynamics: The Finite Vol-ume Method (Pearson Education, 2007), ISBN978013127498364. S. Mazumder,
Numerical Methods for Partial Differ-ential Equations: Finite Difference and Finite VolumeMethods (Elsevier, 2015), ISBN 978012803504765. M. Alves, F. Pinho, P. Oliveira, J. Non-Newton. FluidMech. , 207 (2001)66. C. Rhie, W.L. Chow, AIAA J. , 1525 (1983)67. ˇZ. Lilek, M. Peri´c, Comput. Fluids , 239 (1995)68. P. Colella, J. Comput. Phys. , 171 (1990)69. G. Rempfer, G.B. Davies, C. Holm, J. de Graaf, J.Chem. Phys. , 044901 (2016)70. M. Kuron, G. Rempfer, F. Schornbaum, M. Bauer,C. Godenschwager, C. Holm, J. de Graaf, J. Chem.Phys. , 214102 (2016)71. N. Rivas, S. Frijters, I. Pagonabarraga, J. Harting, J.Chem. Phys. , 144101 (2018)72. C. Godenschwager, F. Schornbaum, M. Bauer,H. K¨ostler, U. R¨ude, A framework for hybrid paral-lel flow simulations with a trillion cells in complexgeometries , in
Proceedings of the International Con-ference on High Performance Computing, Networking,Storage and Analysis (Association for Computing Ma-chinery, 2013), p. 3573. M. Bauer, J. H¨otzer, D. Ernst, J. Hammer, M. Seiz,H. Hierl, J. H¨onig, H. K¨ostler, G. Wellein, B. Nestleret al.,
Code Generation for Massively Parallel Phase-Field Simulations , in
Proceedings of the InternationalConference for High Performance Computing, Net-working, Storage and Analysis (Association for Com-puting Machinery, 2019), pp. 1–3274. K.J. Burns, G.M. Vasil, J.S. Oishi, D. Lecoanet, B.P.Brown, Phys. Rev. Res. , 023068 (2020)75. H. Jasak, A. Jemcov, Z. Tukovic et al., OpenFOAM: AC++ library for complex physics simulations , in
Inter-national Workshop on Coupled Methods in NumericalDynamics (IUC Dubrovnik Croatia, 2007), Vol. 1000,pp. 1–2076. M. Bauer, H. K¨ostler, U. R¨ude, arXiv preprintarXiv:2001.11806 (2020)77. N. Waters, M. King, Rheol. Acta , 345 (1970)78. S.C. Xue, R. Tanner, N. Phan-Thien, J. Non-Newton.Fluid Mech. , 33 (2004)79. K.S. Park, Y. don Kwon, Korea-Australia Rheol. J. , 47 (2009)80. I. Ginzbourg, P. Adler, J. Phys. II , 191 (1994)81. K. Yapici, B. Karasozen, Y. Uludag, J. Non-Newton.Fluid Mech. , 51 (2009)82. R. Sousa, R. Poole, A. Afonso, F. Pinho, P. Oliveira,A. Morozov, M. Alves, J. Non-Newton. Fluid Mech. , 129 (2016)83. S. Dalal, G. Tomar, P. Dutta, J. Non-Newton. FluidMech. , 59 (2016)84. T.W. Pan, J. Hao, R. Glowinski, Int. J. Numer. Meth-ods Fluids , 791 (2009)85. F. Habla, M.W. Tan, J. Haßlberger, O. Hinrichsen, J.Non-Newton. Fluid Mech. , 47 (2014) 86. B. Thomases, M. Shelley, Phys. Fluids , 103103(2007)87. F. Pimenta, M. Alves, J. Non-Newton. Fluid Mech. , 85 (2017)88. L.P. Fischer, T. Peter, C. Holm, J. de Graaf, J. Chem.Phys. , 084107 (2015)89. K.D. Housiadas, Phys. Rev. Fluids , 103301 (2019)90. A. Castillo, W.L. Murch, J. Einarsson, B. Mena, E.S.Shaqfeh, R. Zenit, Phys. Rev. Fluids , 063302 (2019)91. M. H¨utter, P.D. Olmsted, D.J. Read, arXiv preprintarXiv:2011.02797 (2020)92. W.L. Murch, E.S.G. Shaqfeh, Phys. Rev. Fluids ,073301 (2020)93. J.P. Binagia, A. Phoa, K.D. Housiadas, E.S.G.Shaqfeh, J. Fluid Mech. , A4 (2020)94. S. Saad, G. Natale, Soft Matter15