A Direct Effective Viscosity Approach for Modeling and Simulating Bingham Fluids with the Cumulant Lattice Boltzmann Method

Modeling of fluids with complex rheology in the lattice Boltzmann method (LBM) is typically realized through the introduction of an effective viscosity. For fluids with a yield stress behavior, such as so-called Bingham fluids, the effective viscosity has a singularity for low shear rates and may become negative. This is typically avoided by regularization such as Papanastasiou’s method. Here we argue that the effective viscosity model can be re-interpreted as a generalized equilibrium in which no violation of the stability constraint is observed. We implement a Bingham fluid model in a three-dimensional cumulant lattice Boltzmann framework and compare the direct analytic effective viscosity/generalized equilibrium method to the iterative approach first introduced by Vikhansky which avoids the singularity in viscosity that can arise in the analytic method. We find that both methods obtain similar results at coarse resolutions. However, at higher resolutions the accuracy of the regularized method levels off while the accuracy of the direct method continuously improves. We find that the accuracy of the proposed direct method is not limited by the singularity in viscosity indicating that a regularization is not strictly necessary.


Introduction
A popular way to model the stress of non-Newtonian fluids is by imposing an effective local viscosity. Unlike a Newtonian fluid with constant viscosity, this effective viscosity cannot be drawn in front of the gradient of the stress tensor of the Navier-Stokes equation. Implementing non-Newtonian fluids in Navier-Stokes solvers hence either requires storing the stress tensor explicitly or applying the chain rule to the effective viscosity field which can result in rather complicated differential operators. Being a method derived from kinetic theory, the cumulant LBM intrinsically stores the stress tensor separately from the primitive variables as second-order cumulants. Implementing non-Newtonian behavior through an effective viscosity thus comes naturally in the LBM with all differential operators applied in a consistent order.
Many substances of industrial interest are subject to non-Newtonian stress-strain relationships. One such example is fresh concrete [1]. The development of additive manufacturing techniques in civil engineering construction [2] depends critically on the availability of accurate and efficient numerical models for concrete. Among the many properties of freshly mixed concrete, the yield stress behavior is of particular importance in additive manufacturing processes which highly depends on a controlled transition from the fluid to the solid phase.
In this work we discuss the implementation of the Bingham fluid [3] in the context of the cumulant LBM [4]. The Bingham fluid model is a simple yield stress model with linear stress-strain relationship. Below the yield stress the Bingham fluid essentially behaves like a solid. In an effective viscosity approach this is modeled by assuming an essentially infinite viscosity. The transition between non-yielded and yielded behavior appears as a singularity in the effective viscosity model which apparently requires regularization. The most popular regularization is due to Papanastasiou [5] which relaxes the singularity by an exponential function. Ginzburg and Steiner [6] used Papanastasiou's regularization to implement a Bingham fluid in the LBM. For computing the shear stress they used the non-equilibrium second-order moments which makes the method local and hence computationally efficient. However, for computing the shear stress locally in the LBM the viscosity has to be known. Ginzburg and Steiner used the local viscosity from the previous time step. Vikhansky [7] used an iterative method to recover the correct effective viscosity without further regularization.
Tang et al. [8] introduced the He-Luo pressure based formulation [9] for a Bingham fluid in order to eliminate compressibility effects. For improved stability, Chen et al. implemented the Papanastasiou regularization in the multiple relaxation time (MRT) LBM framework. Further studies on Bingham fluids simulated with the LBM are described in [10] [11] [12]. Most of these works do not explicitly mention how the shear rate is numerically determined. However, the determination of the shear stress in the computation of the effective viscosity is an efficiency concern especially for large scale simulations. The local availability of the shear rate is one of the celebrated advantages of the LBM, but viscosity has to be known a priori to take advantage thereof. The two possibilities found in the above mentioned literature are either the one proposed by Ginzburg and Steiner using the viscosity of the previous time step or the one by Vikhansky solving an implicit relationship for each grid node at each time-step. From these two op-tions the implicit one is more compelling as the use of a previous time step in an Eulerian setting can cause many problems. For example, if no advection equation for the viscosity is solved, the method cannot be Galilean invariant since the viscosity of a moving fluid is stored at a stationary grid. On the other hand, being fully explicit is often considered a computational advantage of the LBM. Being completely explicit implies (among other advantages) that the number of operations per time step is constant such that the wall clock time can be precisely predicted and, most importantly, maintained over the entire time simulated. This is by no means a minor concern, as the load balancing in high performance parallel computation relies heavily on the predictability of the run time of all subroutines. The works discussed above pay little to no attention to the parallel scalability of their approach. Most of them consider only two dimensional cases. In addition, none of these works use state of the art lattice Boltzmann kernels, but instead resort to simple single or multiple relaxation time variants.
In recent years, the LBM saw some significant evolution towards improved accuracy and stability. This improvement is due to the usage of tensor product lattices (i.e. using 27 discrete velocities in 3D), moment matching equilibria incorporating terms of higher than second-order in velocity, multiple relaxation times, Galilean invariant moment transformations (i.e. central moments or cumulants) and a better understanding of statistical independence of the different moments. These novel lattice Boltzmann methods include entropic schemes which typically implement tensor product lattices and tensor product equilibria and are primarily aimed at improving stability [13] [14] as well as central moment methods including so-called cascaded schemes [15] [16] [17] aiming particularly at improving Galilean invariance due to the observation that they are much more stable than moment methods in a static frame. Early multiple relaxation time methods [18] [19] did not pay proper attention to the correct segregation of moments and advocated the use of orthogonal transformation matrices. Even though, some early method used a more appropriate weighted orthogonality condition [20], it is comparatively recent that the requirement of statistical independence between observable quantities relaxing with different relaxation times has been appreciated. The cumulant model is based on the statistical independence concept [4] [21] [22] [23]. It incorporates all of the above mentioned recent developments improving both stability and accuracy of the LBM.
In this paper we present an implementation of the Bingham fluid in the context of the cumulant LBM. Our approach starts similar to the one of Vikhansky [7]. We observe that the implicit problem solved iteratively by Vikhansky actually has an analytic solution. However, this analytic solution has a singularity in viscosity which manifests in a violation of the linear stability constraint of the LBM. Surprisingly, this violation is avoided by the approximate solution of Vikhanski even at arbitrary iteration depth. That is to say, the approximate solution does not converge to the analytic solution at positions where the analytic solution is not physically meaningful. Using the approximate solution can hence be regarded as a special form of regularization. However, to the best of our know-ledge the unregularized analytic solution of the implicit problem has never been systematically compared to the approximate solution. With its violation of the linear stability constraint, the analytic solution could be expected to yield unstable results. This, however, is not observed to be the case. In this paper we show that the effective viscosity ansatz can be reinterpreted as a generalized equilibrium approach for which no violation of stability is indicated. We test our model with two simple planar flows: the flow between two infinite plates and the flow between a rotating and a stationary cylinder. Even though these cases are relatively simple, they both cover the singularity problem of the effective viscosity and include extended areas in which the fluid is expected to behave as a solid.
In the following, we will briefly introduce the cumulant LBM, discuss the implementation of the Bingham fluid and the regularization including a fixed-point iteration scheme which removes the singularity. Next we will discuss an alternative interpretation of the effective viscosity model as a generalized equilibrium model in which the singularity disappears. This is followed by a numerical comparison between the two possibilities and conclusions.

Cumulant Lattice Boltzmann Model
The effective viscosity of non-Newtonian fluids can vary substantially such that it is imperative to use a model base with a large range of attainable viscosities.
Early lattice Boltzmann models based on single relaxation time collision operators [24] had a rather limited range of viscosity due to stability issues. In the last two decades substantial progress has been made, in particular regarding the extension of LBM to lower viscosity limits. It is of note that by proper scaling of the grid spacing and the time step any consistent numerical method with limited stability can, in principle, model any viscosity. However, if viscosity is not constant throughout the domain, the ratio between the smallest and the largest attainable viscosity determines which viscosity range can effectively be simulated. Therefore, it is important for the underlying fluid solver for non-Newtonian fluids to have good accuracy and stability properties over a large range of viscosities. In this paper we use the cumulant LBM [4] which is among the LBM with the largest stability range as well as with the highest accuracy for the given D3Q27 stencil used. The cumulant LBM discretizes the Boltzmann equation in the hydrodynamic limit. It is a discrete velocity model using 27 discrete speeds on a tensor product lattice. The discrete velocity distribution function ijk f streams along the discrete velocity directions x y z ie je ke + + in discrete time steps t ∆ from lattice node to lattice node. At lattice nodes the distribution undergoes a local collision which is performed in cumulant space. Cumulants are the statistically independent observable quantities of a statistical distribution.
Previous lattice Boltzmann models either applied the collision operator directly in distribution space [24] or in moment space [18]. Moments can be defined in a large variety of ways [25]. Important properties of a moment base are the chosen frame of reference (e.g. Eulerian leading to raw moments or Lagrangian implying central moments), their orthogonality properties ((non-)orthogonal, weighted orthogonal (i.e. Hermite moments)), and their grouping according to symmetry properties (rotational invariance). Performing collisions in moment or cumulant space allows for a substantially improved control of the properties of the method as different quantities can evolve on different time scales when using different relaxation times. However, for this to be consistent, the different moments have to be statistically independent. Traditionally this is attempted to be enforced by orthogonalization. However, orthogonalization has no unique definition for moments. Moments orthogonal in one frame of reference are not orthogonal in another frame of reference. Even in a Lagrangian frame of reference, moments which are orthogonal in equilibrium are not orthogonal when out of equilibrium.
Cumulants overcome this inconsistency as they are derived from the definition of statistical independence.
The cumulant collision operator assigns an individual relaxation rate abc ω to each cumulant abc C such that the post collision cumulants abc C * can be computed as: ( ) eq abc abc abs abc abc where the equilibrium cumulant eq abc C is typically zero for non-conserved cumulants but can also be chosen otherwise in certain circumstances, for example if energy conservation is not considered as is usually the case for the incompressible Navier-Stokes equation. Asymptotic analysis [26] of the lattice Boltzmann equation reveals two relationships of interest for the modeling of non-Newtonian fluids.
The viscosity is related to the relaxation rate of second-order cumulants (i.e. ω ω ω ω = = = and the corresponding diagonal terms of second-order of the tensor abc ω with the exception of the trace): The shear stress in the LBM can be written as: For a known viscosity the shear rate can be locally computed from the second-order cumulants: Here we omitted the other permutations of Equations (5) and (6). The LBM does not require any additional finite differencing for the calcula-tion of the shear rate. However, it is of note that the relaxation time, which is a function of viscosity, appears in the equations for the shear rate such that an implicit problem arises if the viscosity depends on the shear rate.

Bingham Fluid
The Bingham fluid [3] extends a Newtonian fluid by a constant yield stress 0 τ : Our aim is to impose Equation (7) onto the LBM. This is done through equating Equations (3) and (7): According to Equation (5) γ is a function of ω . To make this relationship explicit, we define: With this Equation (9) can be stated as: This can be solved for the relaxation rate to give explicitly: We will call the relaxation rate in Equation (13) . According to Equation (2) viscosity goes to infinity when 0 ω = and viscosity will be negative for 0 ω < . Regularization procedures like the one by Vikhansky [7] are imposed to avoid the singularity in viscosity.

Regularization
From Equation (8) it is not obvious why the effective viscosity ν should become negative before the shear rate γ reaches zero. The reason for the problem encountered in the last subsection is in the implicit dependence of γ on the relaxation rate when computed from the second-order cumulants. If we ignore this implicit dependence and solve Equation (9) for ω we obtain: This solution is always positive, but it requires γ which in lattice Boltzmann is only known through Instead of using the apparently unfeasible analytic solution, we can solve Equation (14) through a fixed-point iteration as follows [7]: Since all terms on the right hand side of Equation (16) are positive and since we start the iteration from a positive value 1 n ω + is always guaranteed to be positive.
In order to accelerate the execution of the fixed-point iteration we rewrite Equation (16): We note that Equation (17) has a simple analytic solution for 1 n n ω ω + = which can easily be verified to be identical with Equation (13) which we are trying to avoid. In fact a solution of Equation (17) exists only for 0 n ω ≠ which is the mathematical reason for the analytic solution appearing to be unfeasible. The obtained iteration can be easily unrolled for a fixed iteration depth. For example, the iteration of depth five reads: In his original paper [7] Vikhansky claims that two to three iterations are sufficient to reach convergence. However, it is unclear to us whether he started from the relaxation rate of the previous time step which would certainly reduce the number of iterations but would come at additional cost in memory consumption and data transfer and would not be Galilean invariant. It is of note that the fixed-point iteration is not expected to converge quickly. The convergence order can be computed directly by comparing two subsequent iterations: The rate of change is constant implying linear convergence. However, we should remind ourselves that we deliberately replaced the exact solution by this approximation for the purpose of regularization and avoiding the singularity of the exact solution. The use of faster converging methods like Aitken extrapolation or Newton iteration would reintroduce the possibility for the relaxation rate to become negative and could hence eliminate the advantage of regularization. That is, of course, if regularization is really as advantageous as implied through the avoidance of the singularity in viscosity.
Naively, one could assume that a singularity in viscosity and the violation of the linear stability regime should result in unstable simulations. However, in the current case, as demonstrated below, this problem did not manifest in our simulations. A simple reason why the violation of the linear stability constraint might not be relevant in the current case is seen from the fact that according to Equation (13) a problem occurs only if the shear rate is small enough. Instability would necessarily manifest itself in a locally increased shear rate which would force the relaxation rate back into its stable range. Thus, even though the linear stability constrained is violated, this violation is self-limiting and consequently non-linearly stable.
We also have to take the nature of the LBM into consideration and how it compares to a classical Navier-Stokes solver. In a classical Navier-Stokes solver the stress tensor is constructed by computing the shear rate and multiplying it with the viscosity. A singularity in viscosity is numerically unmanageable in this context and obviously has to be avoided. However, in the LBM, viscosity is an emergent property of the relaxation rate. The singularity occurs when 0 ω = which is not necessarily a problem for the lattice Boltzmann algorithm itself as it does not divide by ω . The problem of a singular viscosity is hence a very serious numerical difficulty for the Navier-Stokes equation, but it is not necessarily relevant for the LBM.
It is instructive to recall that the effective viscosity ansatz is only an algebraic trick to express the yield stress through the relaxation rate. Alternatively we could specify the yield stress in the equilibrium second-order cumulants and leave the relaxation rates untouched. A similar dual approach exists in the correction for cubic error terms in viscosity of the LBM. Dellar [27] decided to incorporate this correction into a modified relaxation rate while we implemented it by a modified equilibrium function of second-order cumulants [4]. Even though these approaches appear to be rather different in concept, a detailed investigation shows that both methods are algebraic transformations of each other and should hence yield identical results up to round-off errors. In the following we demonstrate the equivalence of the modified relaxation rate to a generalized equilibrium.

Equivalence of Effective Relaxation Rate and Generalized Equilibrium Approaches
Asymptotic analysis as, for example, demonstrated in appendix G of [4] can be used to calculate the relationship between second-order cumulants and the velocity derivatives: Note that Equation (20) and (21) are identical to Equation (5) and (6), respectively (up to an unspecified equilibrium cumulant each). The equilibrium here is understood as a generalized equilibrium as introduced by Asinari [28] for ana-lyzing the cascaded LBM. The generalized equilibrium is an attractor for the collision operator and does not necessarily have a direct thermodynamic interpretation. For the respective second-order cumulants the thermodynamic equilibria would be zero which is why they are omitted from Equations (5) and (6). Here we include them for the purpose of showing that an effective relaxation rate can be equivalently understood as a generalized equilibrium. For this purpose we consider Equation (5) with the equivalent relaxation rate a ω resulting from the rheology model and the fixed Newtonian relaxation rate ω ∞ : It is hence seen that the equivalence between the modified relaxation rate in Equation (22) and the generalized equilibrium at fixed relaxation rate in Equation (22) is established by choosing the specific equilibrium Equation (24). It is of particular note here that this relationship is independent of how a ω was obtained in the first place, implying that any method using an effective relaxation rate can be rewritten as a method with fixed relaxation rate and generalized equilibrium. This hence applies to Vikhansky's model for Bingham fluids [7] in the same way as for Dellar's Galilean correction to the viscosity [27]. The inverse, however, is not true as we see that the generalized equilibrium depends on the cumulant itself implying that only generalized equilibria of a very specific form can be turned into an effective relaxation rate.
In the generalized equilibrium form, the linear stability constraint is not necessarily violated if 0 a ω < . It hence seems to be admissible to use the analytic solution of the effective relaxation rate even if it might become negative.

Curved Boundary Conditions
Being a Cartesian grid based method, the LBM is applicable to real world industrial problems only because effective methods to incorporate grid refinement and curved boundary conditions have been incorporated. Most of these methods are dependent on the local viscosity and their application to non-Newtonian fluids requires some adjustment. We leave the important topic of grid refinement to future work and discuss the implementation of an off-grid velocity boundary conditions applicable to moving curved walls. In lattice Boltzmann, velocity boundary conditions are conveniently implemented via a bounce back scheme with modified velocity. Bounce back is applied on links entering the fluid domain during the streaming step. In its simplest implementation, bounce back returns the population leaving the domain as the population entering the domain in opposite direction. This results in a zero velocity boundary condition approximately halfway between grid nodes. The Open Journal of Fluid Dynamics no-slip (i.e. bounce back) boundary is turned into a velocity boundary condition by adding twice the momentum of the moving wall [29]: Here ( ) x qe x t q t qf x t f x t q t q q q * + ∆ + ∆ + ∆ + = + + + (28) It is observed that the population arrives at the correct location. However it appears to arrive at the wrong time instance ( ) . This is of minor concern as the LBM recovers the Navier-Stokes equation only in diffusive scaling The non-Newtonian fluid model enters the boundary condition only in the determination of ω . This can be implemented for the iterative as well as for the analytic model of the Bingham fluid's viscosity.

Implementation Details
The Bingham fluid model with unrolled iterative regularization according to Equation (18) and the method based on the analytic solution are implemented in the massively parallel lattice Boltzmann solver Virtual Fluids [34]. The base model is the cumulant lattice Boltzmann solver with 27 velocities and optimized relaxation rates as introduced in [35] and analyzed in [36] [37]. In this method, the relaxation rates of cumulants of order three and higher are selected depending on the relaxation rates of second-order cumulants and this optimization is only valid in a specific range of relaxation rates associated with low viscosity. For the high viscosity encountered below the yield threshold the optimization would introduce additional singularities in the higher cumulants. This is avoided here by selecting the relaxation rates of the higher-order cumulants according to the fixed asymptotic relaxation rate ω ∞ which means that we are applying the generalized equilibrium instead of the effective viscosity ansatz. A caveat to this choice is that ω ∞ is somewhat arbitrary such that the optimization obtained in [35] for Newtonian fluids might not be optimal for non-Newtonian fluids.
However, the complex optimization procedure in [35] heavily relies on the linearization of the LBM and requires substantial further research to be fully applied to the non-linearity of the Bingham fluid. We will hence leave a Bingham specific optimization of the relaxation rates to future research.
The lattice Boltzmann algorithm including the computation of the effective relaxation rate for boundary conditions and the collision operator is given in

Numerical Validation
The main purpose of the current paper is to demonstrate that the viscosity singularity which has to be avoided in Navier-Stokes implementations of the effective viscosity model does not arise in the LBM. For demonstration purpose we chose two canonical test cases which are both popular and relevant for Bingham fluids: Poiseuille flow and Taylor-Couette flow. Yet, even for these comparatively simple setups proper convergence studies for yield stress fluids are rarely found in literature.

Poiseuille Flow
A classical and simple test for the implementation of a Bingham fluid is flow between two infinite plates driven by a constant force. This Poiseuille type flow results in a parabolic flow near the boundaries and a plug flow in the center of the channel for a Bingham fluid as described in [38]: where x F is the force amplitude driving the flow, 0 z is the yield point given , H is the half height of the channel and µ ρν = is the dynamic viscosity. Driving the flow by a force x F is typically preferred in benchmark simulations over the more physically correct pressure gradient as the latter cannot be implemented with periodic boundaries and would hence introduce boundary effects from the pressure boundary conditions. We used a simulation domain of size 2 2 2 H × × . The viscosity at the lowest resolution µ was set to 0.005, the yield stress  Bingham numbers defined respectively as: The speed of sound is a constant in the LBM such that consistency between the dimensionless numbers is established by scaling . This scaling implies that the time step t ∆ is scaled proportionally to the grid spacing x ∆ which also keeps the velocity scale U constant.
The disadvantage of this so-called acoustic scaling is that a finite error in Mach number persists such that absolute convergence is not expected. In order for this to be small we chose rather small values for the force and the yield stress.
To quantify the error we compute the following L2 norm ( ) ( )

Taylor-Couette Flow
Our next test is a planar Taylor-Couette flow depicted in Figure 3. We investigate the flow between a rotating outer and a resting inner cylinder. Perpendicular to the plane periodic boundary conditions are used such that the setup is quasi two dimensional. Interpolated bounce back with and without velocity is used for the outer and the inner cylinder, respectively. In this test case, flow attached to the outer cylinder will be below the yield threshold and hence move as a solid with the same angular velocity as the outer cylinder. This plug flow domain attached to the outer cylinder is reduced in size with higher angular momentum.

Semi-Analytic Solution
A semi-analytic solution for the laminar Taylor-Couette flow of a Bingham fluid is usually derived in dimensionless form scaled with the Reynolds number of the inner cylinder [39] [40]. This unfortunately precludes the case where the inner cylinder is at rest. We hence present here a semi-analytic solution for this particular case. Following Landry et al. [40] we start from the condition for stationary momentum in cylindrical coordinates: The stress in the Bingham fluid is given in cylindrical coordinates as: where ( ) V r is the tangential velocity at radius r. For the purpose of this derivation we assume 0 rθ τ τ > and plug Equation (36) into Equation (35). This gives rise to a 2nd order ordinary differential equation which we supplement Open Journal of Fluid Dynamics With ( ) The yield radius y R is unknown and it is recovered from the condition of continuity for Equation (37). Basically we need to solve for y R . This is done here numerically such that the final solution is only semi-analytic. As there are multiple solutions we have to pick the one for which y i R R > . Figure 4 shows the profiles of the angular velocity in between the two cylinders for four different angular velocities of the outer cylinder, for four different resolutions (

Results
) and the two methods under investigation. The yield stress is We observe that a plug flow attached to the outer cylinder develops for sufficiently low angular velocities and that the different resolutions and the analytic and iterative methods agree well with respect to the velocity profiles. Figure 5 depicts the convergence behavior of the L2 norm of the angular velocity profiles of the Taylor-Couette solution when compared to the semi-analytic result. It is observed that both methods under investigation show very similar results at coarse resolution. However, at higher resolutions the convergence of the iterative method slows down substantially. It is of note that the convergence of the LBM using the analytic relaxation rate also deteriorates at the highest resolution, especially for the cases with a high proportion of non-yielded fluid in the domain. However, this deterioration is much less pronounced than in the case of the iterative method.

Conclusions
In this work we presented a cumulant LBM for the simulation of Bingham fluids. The ansatz used is similar to the one proposed by Vikhansky, however, in contrast to his work we solve the implicit problem analytically. The resulting model implies the existence of negative viscosities at shear rates below the yield threshold. At higher shear rates the viscosity returns to a positive value such that the method remains stable through self-limiting even though it appears to be linearly unstable. The singularity in viscosity occurring between the yielded and non-yielded state does not affect the performance of the lattice Boltzmann method much since viscosity is never explicitly used in the method and appears only as an emerging property derived from the relaxation rate. This is in stark contrast to an effective viscosity model implemented in a Navier-Stokes solver where the viscosity has to be explicitly specified such that the singularity would be fatal to the simulation.
We showed that an effective viscosity model can always be translated to a generalized equilibrium model with fixed viscosity. This result can also be applied in other contexts, most notably for Dellar's Galilean correction of viscosity using a modified relaxation rate [27].
Using the analytic solution for the relaxation rate has two important advantages over Vikhansky's method. First, it is simpler and more efficient since no iteration is required and the relaxation rate can be obtained by evaluating a simple analytic expression. Second, the analytic method is more accurate. The latter could appear obvious as we are evidently comparing the exact solution to a slowly converging approximation. However we have to recall that the necessity of some sort of regularization for the singular viscosity is almost universally accepted in literature. While the loss of convergence due to regularization in the case of high resolution is not surprising, we could show that regularization can be avoided altogether due to the way of how the LBM deals with viscosity and that avoiding regularization improves convergence and therefore overall accuracy. Finally, we note that more complex fluids will still require iterative function solvers to determine the effective viscosity or the generalized equilibrium. Even using the analytic relaxation rate, the convergence properties of our me-thod are below the expectation of a second order accurate method. Simulating a yield stress fluid includes the modeling of a quasi-solid phase which, for an explicit method on a Cartesian mesh, is not a well posed problem. Convergence studies are rarely shown for the simulation of such complex fluids. As the solid phase remains to be an approximation even with the analytic relaxation rate, it is essential for the underlying numerical method to support a high viscosity contrast. The cumulant LBM has been established for fluids with very low viscosity.
In this paper we showed that the simulation of the solid domain is not limited by the requirement of a finite viscosity in this method.