Three-dimensional simulations of Bingham plastic flows with the multiple-relaxation-time lattice Boltzmann model

ABSTRACT This paper presents a three-dimensional (3D) parallel multiple-relaxation-time lattice Boltzmann model (MRT-LBM) for Bingham plastics which overcomes numerical instabilities in the simulation of non-Newtonian fluids for the Bhatnagar–Gross–Krook (BGK) model. The MRT-LBM and several related mathematical models are briefly described. Papanastasiou’s modified model is incorporated for better numerical stability. The impact of the relaxation parameters of the model is studied in detail. The MRT-LBM is then validated through a benchmark problem: a 3D steady Poiseuille flow. The results from the numerical simulations are consistent with those derived analytically which indicates that the MRT-LBM effectively simulates Bingham fluids but with better stability. A parallel MRT-LBM framework is introduced, and the parallel efficiency is tested through a simple case. The MRT-LBM is shown to be appropriate for parallel implementation and to have high efficiency. Finally, a Bingham fluid flowing past a square-based prism with a fixed sphere is simulated. It is found the drag coefficient is a function of both Reynolds number (Re) and Bingham number (Bn). These results reveal the flow behavior of Bingham plastics.


Introduction
Bingham plastics or viscoplastic materials -such as cement mortar, slurries, and suspensions -are widely used in hydraulic and civil engineering. The most important characteristic of Bingham plastics is their yield stress, which defines the point above which flow occurs. Mathematical models for Bingham plastics include the Bingham model, the Herschel-Bulkley model, and the Casson model (Mitsoulis, 2007). The traditional methods of computational fluid dynamics (CFD), such as the finite-volume method (Neofytou, 2005;Xu, Yuan, Repke, & Wozny, 2012) and the finite-element method (Bell & Surana, 1994;Ozmen-Cagatay & Kocaman, 2011), are the most common numerical techniques used to simulate Bingham plastic flows. However, such traditional methods often produce inaccuracies in cases of complex geometries and boundary conditions (I.-B. Lee et al., 2013;Mendoza, Succi, & Herrmann, 2013;Riddle, Carruthers, Sharpe, McHugh, & Stocker, 2004). Despite the adoption of special techniques (Tezduyar, 2001;Tezduyar, Takizawa, Moorman, Wright, & Christopher, 2010) or complex grids (Deng, Mao, Tu, Zhang, & Zhang, 2012) to get round potential inaccuracies, simulations are CONTACT Qi-Cheng Sun qcsun@tsinghua.edu.cn cumbersome, thus prompting the development of a more efficient numerical method for Bingham plastics. The lattice Boltzmann model (LBM) has been a successful substitute for traditional CFD methods in various liquid flows (Haghani, Rahimian, & Taghilou, 2013; C.-H. Lee, Huang, & Chiew, 2015). One of the most popular LBM models is the Bhatnagar-Gross-Krook (BGK) model, in which only one single-relaxation time is used (S. Chen & Doolen, 1998). Given its simplicity, the BGK is an efficient numerical method for non-Newtonian fluids. However, the BGK has a significant drawback concerning low-viscosity fluids where it is often the case that simulations are hampered by numerical instabilities. To overcome this disadvantage, Chai, Shi, Guo, and Rong (2011) introduced the multiple-relaxationtime LBM (MRT-LBM) to model non-Newtonian fluids in 2D. The MRT-LBM achieved a much better performance in terms of numerical stability as different relaxation times can be individually adjusted in accordance with specific problems (Lallemand & Luo, 2000).
The LBM has several advantages, including a simple evolutional equation which is easy to parallelize with a high parallel efficiency . Wang, Zhang, Bengough, and Crawford (2005) successfully simulated an incompressible flow in porous media through a cell-based domain-decomposition on different clusters. The parallel implementation and scaling results of a hybrid lattice Boltzmann (LB)/finite-element code for suspension flow simulations at the Argonne National Laboratory and IBM's Blue Gene/P were presented and discussed by Clausen, Reasor, and Aidun (2010). A flexible patch-based LB parallelization approach integrated into the waLBerla software framework was used to compare the parallel efficiencies of heterogeneous graphics processing unit (GPU) and central processing unit (CPU) clusters (Feichtinger et al., 2011). The LBM is effective in realistic simulations at different parallel efficiencies, which depend on the LBM codes and the architecture of the clusters used. However, in their work, only the BGK model was used in their test and comparison of parallel efficiency. In this study, therefore, we took the next logical step to verify and compare parallel performance using the MRT model. This paper presents a three-dimensional (3D) parallel MRT-LBM for Bingham plastics as an extension of previous research. The 3D MRT-LBM and the Papanastasiou model (Papanastasiou, & Boudouvis, 1997) for Bingham fluids are introduced in section 2 and the impact of relaxation parameters on the MRT model for a Bingham fluid is discussed. Section 3 validates the 3D MRT-LBM for Bingham fluids and the results are presented. Section 4 discusses the parallel performance of the 3D MRT-LBM according to a benchmark cavity flow. Section 5 examines a Bingham fluid flowing past a square-based prism with a fixed sphere and conclusions are drawn in section 6.

D3Q19 MRT-LBM
D' Humieres (1994) first proposed the MRT-LBM, also known as a generalized LBM. A few years later, a 3D MRT-LBM was developed (d 'Humieres, Ginzburg, Krafczyk, Lallemand, & Luo, 2002), which demonstrated better numerical stability compared to a BGK-LBM for a Newtonian fluid. In this paper, a 3D MRT-LBM with a D3Q19 model is used. The evolution equation is: ( where is the density distribution function of a fluid particle at position x and time t with velocity c i , and f i eq is the equilibrium distribution function of f i (x, t), which is given as The weighting factors ω i in the D3Q19 model are given as .
( 3 ) The corresponding particle velocities are . (4) The only difference between MRT and BGK is the collision matrix = M −1 SM, where M is the transformation matrix. The role of M is to change f and f eq from the velocity space to the moment space with m = M · f and m eq = M · f eq , where m denotes the moment vector: and the equilibrium value of m eq is given as follows: Details of M, mand m eq can be found in d 'Humieres et al. (2002).

Papanastasiou model for Bingham fluids
Different constitutive equations have been proposed to model the stress-deformation behavior of viscoplastic materials (Bird, Dai, & Yarusso, 1983), among which the Bingham model is most popular due to its simplicity and can be shown as: where τ and τ 0 are the shear stress and yield stress, respectively,γ is the shear rate, and μ p represents the plastic viscosity. Papanastasiou and Boudouvis (1997) modified the equation as follows: where m is introduced to overcome the discontinuity at τ 0 . In this paper, m is set to 10 8 (S.-G. Chen, Sun, Jin, & Liu, 2014). From Equation (11), the apparent viscosity can be obtained: where |γ | = 2 γ = 2 3 α,β=1 Chai et al. (2011) and Chai and Zhao (2012) proved that S αβ is a function of the density distribution function: The relaxation factor s 9 , s 13 can be obtained using Equation (7) once the apparent viscosity is known. The criterion used to track the flow (yielded) of Bingham fluids is when the magnitude of the extra stress |τ | exceeds the yield stress τ 0 .

The effect of relaxation parameters in the MRT-LBM for Bingham fluid
For Newtonian fluid, it is understood that the MRT-LBM can overcome numerical instabilities at very low viscosities. However, the potential advantages of applying the MRT-LBM to Bingham fluid have not, so far, been studied. In this section, we examine the effect of relaxation parameters in MRT-LBM for the Bingham fluid using a numerical approach. For the sake of simplicity, we first consider a D2Q9 MRT-LBM. The results can easily be extended to 3D.
In the MRT-LBM, the density distribution functions, which have no physical meaning, can be decoupled through transformation. The new moment representations have an obvious physical significance, such as hydrodynamic quantities and their fluxes, etc. Thus, each relaxation parameter can be controlled independently using this mechanism.
In the D2Q9 MRT-LBM there are nine relaxation parameters: S = (s 0 , s 1 , s 2 , s 3 , s 4 , s 5 , s 6 , s 7 , s 8 ). We consider a classical benchmark -Poiseuille flow in dimensionless form. The boundary conditions are the same as those in S.-G. Chen et al. (2014). The equivalent pressure gradient is −∇P = 1.667 × 10 −5 , the plastic viscosity μ p = 0.24, and the yield stress τ 0 = 0.004. The lattice spacing δ x = 0.001 m, the time step δ t = 0.001 s, and the density of the Bingham fluid ρ 0 = 1000 kg/m 3 . Firstly, the following two cases were compared (Figure 1(a)). The relaxation parameters were the same as in our previous work for case 1, where S = (0, 1.1, 1.0, 0, 1.2, 0, 1.2, 1/ω, 1/ω), whereas all relaxation parameters were equal to 1/ω for case 2. It can be shown that the BGK-LBM cannot converge to an analytical solution, especially in the unyielded region where the shear rate is close to zero. On the contrary, the MRT-LBM solution agrees very well with the analytical solution, with an average error of only 0.25%.
Considering that s 0 , s 3 and s 5 are related to the conserved quantities of density ρ and linear momentum j x , j y , collisions in the moment space do not change these conserved quantities. Thus, in case 3 s 0 , s 3 and s 5 are set to 0, while other relaxation parameters are set to 1/ω. It can be seen from Figure 1(b) that the deviations in the unyielded region are smaller than in case 2. Therefore, some higher-order moments have important effects for a Bingham fluid, which are often insignificant in a Newtonian fluid simulation.
unyielded region can be found in Figure 1(c). However, the maximum velocity is not in good agreement with the analytical solution.
The relaxation parameter s 2 , which is related to the fourth-order moment ε, is considered in case 5 and set to 1.0. It is obvious from Figure 1(d) that the result (the average error is 1.8%) is much better than in the other cases. The relaxation parameter s 1 is determined by the bulk viscosity and s 7 and s 8 are determined by shear viscosity. If s 1 is equal to 1.1 and s 7 and s 8 are set to 1/ω, the standard MRT-LBM is recovered and a more accurate solution is obtained.
From the numerical analysis of the relaxation parameters it was found that the higher-order moments have a significant effect on the numerical stability in the unyielded region where the shear rate is close to 0. By adjusting the relaxation parameters carefully, the MRT-LBM can overcome this numerical instability in the unyielded region. We conclude, therefore, that the MRT-LBM is more suitable than the BGK-LBM for the simulation of Bingham fluids.
Through the same method as above, the relaxation parameters in 3D can also be determined when the error between numerical and analytical solutions is small enough. In this paper, s 1 = 1.13, s 2 = s 10 = 1.4, s 4 = 1.2, and s 16 = 1.85.

Dimensionless numbers
In this work, the Reynolds number Re and the Bingham number Bn are used to present the simulation results for standard comparisons: where ρ is the density of the Bingham fluid, L c is the characteristic length, and U c is the characteristic velocity. For spherical particles the drag coefficient C d , which is a dimensionless parameter, is usually employed: where |F x | is the drag force in the flow direction and r is the radius of the spherical solid particle. Previous studies (Beris, Tsamopoulos, Armstrong, & Brown, 1985;Blackery & Mitsoulis, 1997;Yoshino, Hotta, Hirozane, & Endo, 2007) have used Stokes' drag coefficient C s to measure the relationship between this dimensionless parameter and the yield stress. The current study also uses C s , which is expressed as

3D Poiseuille flow through the cross-section of a circular tube
A 3D steady Poiseuille flow through a circular tube was used to validate the 3D MRT-LBM. All the conditions for this simulation are in dimensionless form. The simulation domain has 48 × 60 × 60 uniform lattices (Z × 2R × 2R). The pressure boundary conditions put forward by Zou and He (1997) are used at the inlet and outlet. The pressure gradient G = −8.33 × 10 −6 in lattice units. The halfway bounce-back boundary condition is applied to other walls. The lattices outside the tube remain unchanged to save on computational time. The plastic viscosity μ p is set to 0.2, and Re is fixed at 0.95. The Bn varies from 0 to 80 as the yield stress τ 0 increases from 0 to 0.00010. The lattice spacing δ x = 0.001 m, the time step δ t = 0.001 s and the density of the Bingham fluid ρ 0 = 1000 kg/m 3 . The analytical solution is given by Chatzimina, Georgiou, Argyropaidas, Mitsoulis, and Huilgol (2005): where r 0 is the yield point given by r 0 = − 2τ 0 G , and u max = − G 4μ (R − r 0 ) 2 . Figure 2 shows the analytical solutions (solid lines) and numerical results (dots). The velocity profile is a parabola at first, when yield stress is small. As the yield stress increases the profile changes to a flat plateau (i.e., an unyielded region). The numerical results match the analytical solutions well, which proves that the 3D MRT-LBM is an effective tool for the simulation of Bingham fluids.
The accuracy of the MRT-LBM approach has been studied by Chai et al. (2011) and it is second-order accurate in space, which is affected by fluid compressibility. Thus, the maximum Ma (Ma = u max C s ) in their simulation is 0.016 when τ 0 is 0, which is much smaller than 1.0 in order to ensure a reasonably accurate solution.

3D Poiseuille flow through the cross-section of a square tube
The 3D MRT-LBM approach was further employed to simulate the Poiseuille flow through a square tube. The simulation domain has 60 × 60 × 60 uniform lattices; all other computational conditions are the same as those in section 3.1, including the pressure gradient and boundary conditions. As is shown by Figure 3, our results are consistent with those of Vikhansky (2008). The unyielded zones include two types: dead zones at the corners and a plug zone at the center. A single parameter Bi, which is similar to the Bingham number Bn, governs the distribution of the yielded and unyielded zones. The rates of the radius of the unyielded region to the side length at different values of Bi are 0.51, 0.74, and 0.95.

The parallel performance of the MRT-LBM
Equations (1), (8), and (9) indicate that the calculation of the LBM is local and that only information regarding the nearest neighbor cells is used during each time step (Kandhai et al., 1998). This spatial locality makes the LBM very well suited to parallel computing.
An equal-subdomain partitioning technique was used to divide the domain into subdomains and the message-passing interface (MPI) library (Figures 4 and  5) was adopted to transfer data between subdomains (Feichtinger et al., 2011;Wang et al., 2005). A Bingham fluid flowing in a cubic cavity was used as a benchmark to verify the parallel performance of the 3D MRT-LBM. Numerical experiments were conducted on the High-Performance Computing Architecture of the Tsinghua National Laboratory for Information Science and Technology, China. The details of this architecture are as follows: (1) Peak performance of 104 Teraflops with 740 computing nodes (37 racks of 20 nodes). (2) Each node contains two Intel Xeon X5670 (2.93 GHz, 12 MB cache) CPUs, each with six cores (a total of 8880 processor cores; 1 rack = 240 cores).  (3) 32 or 48 GB RAM per node (28.9 TB in total).
(4) Bandwidth of up to 5.1 GB/s.
Generally, the turnaround of one iteration T total is the sum of the times for calculation T cal and communication T com . Thus, where a and b are parameters related to the architecture, U is the computation speed, V is the communication speed, N is the lattice resolution, and n is the number of processors (Wang et al., 2005). Weak scaling and strong scaling (Kandhai et al., 1998) are two of the most common methods for the evaluation of scaling performance. Concerning weak scaling, the subdomain size for one processor remains fixed while the domain size is enlarged. As for strong scaling, the processing core number is increased, whereas the domain size is fixed. Two parameters, the speedup factor Note: The subdomain size on each process is fixed to 20×20×20, 30×30×30, and 40×40×40. The whole domain sizes vary from 20×20×20 to 360×360×360. S and efficiency E, are calculated to assess the parallel computation: where T 1 and T n are the running times when a single processor and n processors are used, respectively. Equation (19) indicates that in weak scaling, N 3 n is constant; thus, T cal remains fixed whereas T com increases as the lattice resolution N increases and T total is a linear function of N 2 . By contrast, in strong scaling, T com and N remain unchanged whereas T cal decreases with an increasing processor number; thus, T total is proportional to 1/n. The weak scaling results are shown in Figure 6. The larger the subdomain size, the longer the time required to run 4000 iterations. The total run time T total , and particularly the communication time T com , increase slowly with an increasing number of cores; T total and N 2 also hold a linear relationship (Equation (19)).
The strong scaling results are shown in Figure 7. For one domain size, the total computing time keeps a perfect linear relationship with the inverse of n, consistent with the theoretical result from Equation (19). The speedup first increases and then saturates as more processors are used for a fixed domain size. When the same number of processors is used, a higher speedup and efficiency can be obtained with a larger computational domain. Both speedup and efficiency decrease as the number of processors increases because more time is spent on communications. Besides, the efficiency E with different processor numbers also relates to the domain size: a larger domain size implies a higher efficiency. Sometimes, the Note: The whole domain size is fixed to 100×100×100, 300×300×300, and 500×500×500. The subdomain sizes vary from 10×10×12 to 500×500×500.
speedup can be larger than the number of processors, and the efficiency is higher than that which results from the optimized use of cache memory.
In conclusion, the proposed 3D MRT-LBM is appropriate for parallel implementation, especially for large domain sizes. The parallel efficiency is around 80% when the domain size is 500×500×500 and 792 processors are used. From Equation (8), the evolution equation of the MRT model is more complex than that of the BGK model, and a matrix needs to be solved during each step. Thus the computational cost of the MRT model is larger; however, the time spent on communication is almost the same in both models. As a result, the parallel speedup or efficiency of the MRT is better than that of the BGK since the percentage of communication time is smaller. The parallel performance of the BGK and MRT models are also dependent on the domain size and the number of processors used. Generally, with a larger domain size and more processors, a better parallel efficiency can be achieved with the MRT model. For example, with a domain size of 500×500×500, the total computational time of the MRT is twice that of the BGK when 24 cores are used, but this decreases to 1.2 times when 792 cores are used. When 24 cores and 792 cores are used, the speedup increases from 24 to 475 and the efficiency decreases from 1.00 to 0.60 for the BGK model, contrasting with a speedup increase from 24 to 607 and a decrease in efficiency from 1.00 to 0.78 for the MRT model. Therefore, the parallel performance of the MRT model is better than that of the BGK model, although the total computational time of the MRT is greater. Concerning parallel efficiency, it is certain that the total computational time of the MRT model can be further decreased through the use of more processors.

Bingham fluid flows around a fixed sphere
A Bingham fluid flowing around a sphere is an important benchmark problem in fluid mechanics. Because no analytical solution yet exists, many researchers (Beris et al., 1985;Blackery & Mitsoulis, 1997;Prashant & Derksen, 2011) have tried to solve this problem using different numerical methods, such as the finite element method and the LBM. In this paper, a solution for the problem is investigated based on the proposed MRT-LBM. The calculated domain is a square-based prism with a size of 4×4×6 (L/d×L/d×H/d). The spatial resolution is such that the spherical radius spans 6 lattices. The sphere is assumed to be stationary, and the walls of the prism and the fluid move at a constant velocity V (Figure 8). Such a resolution is sufficient for achieving a reasonable accuracy (Derksen, 2008;Derksen & Sundaresan, 2007). The circular boundary is approximated by a series of stairs, and the moment exchange method is applied to calculate the force imposed on the sphere.
To ensure a creeping flow, Re is set to 0.001, while Bn is varied from 0.108 to 544.5 to draw comparisons with the results of Blackery and Mitsoulis (1997). In all cases, C s is calculated; the yielded and unyielded regions are shown in Figure 9. These qualitative results are consistent with those reported by Blackery and Mitsoulis (1997).
Aside from the qualitative results, a quantitative term C s is compared with that of Blackery and Mitsoulis (1997), who determined an approximate relationship    Blackery and Mitsoulis (1997).
(21) As shown in Figure 10 and Table 1, the numerical results using the MRT-LBM (square dots) are in good agreement with those derived from Equation (21) (solid lines). The average deviation between the present result and the reference data is only 1.22%. As yield stress increases, the resistance effect of the fixed sphere becomes larger. As a result, C s increases as Bn increases.
As with the 2D case (S.-G. Chen et al., 2014), the drag coefficient C d is composed of two parts: In Equation (22), C dB and C dN are the drag coefficients from Bingham fluid and Newtonian fluid, respectively.
A relation for C dN is provided according to the work of Owen, Leonardi, and Feng (2011): Concerning a spherical particle in a Bingham plastic, for settlement to occur the diameter of the spherical particle has to be larger than a critical diameter d c . Bethmont et al. (2003) have shown that a general form of this criterion can be written as where τ 0 is the yield stress, ρ s is the density of the particle, ρ f is the density of the fluid, and K is a constant. However, the existing data suggest significantly dispersed values for K (between 11 and 25; see Roussel, 2006). For example, Beris et al. (1985) determined that K = 20.97 while Ansley and Smith's (1967) determined that K = 16.5. In this paper, the value of K is presented according to the results of the MRT-LBM. According to Equation (24), the drag force on a spherical particle to overcome the yield stress can be expressed as follows: Introducing Equations (25) and (14) into Equation (15), C dB can be derived theoretically: Combining Equation (26) with Equation (23), the drag coefficient can be expressed as follows:  Table 2. With the data in Table 2, a value of K = 21.61 was determined using the least squares method with an associated minimum error of 7.4%. This result is similar to that of Beris et al. (1985), which was obtained using the finite-element method. Thus, C d = 28.81Bn+21.12 Re + 6.3 Re 0.5 + 0.25. Figure 11 compares the numerical results (unfilled dots) with Equation (28) (colored lines), which further proves that Equation (28) can be applied to predict the drag coefficient of the sphere particle in Bingham flows.  Figure 11. Correlations between C d , Bn, and Re.
The present results indicate that increasing Bn substantially increases the drag coefficient in the system when Re is constant, while increasing Re decreases the drag coefficient for the same values of Bn; Newtonian fluids are also characterized by this phenomenon.

Conclusions
A 3D MRT-LBM has been proposed for Bingham fluids in this paper. Numerical evidence supporting the assertion that the MRT-LBM is better than the BGK-LBM in terms of stability and accuracy has been presented. The MRT-LBM was validated using the Poiseuille flow, which served as a benchmark. The flat plateau of the velocity profile widens as the yield stress increases. The consistency between the numerical and analytical solutions validates the 3D MRT-LBM. The parallel computing experiments show that the proposed model is appropriate for parallel implementation, with a high level of parallel efficiency. The application of the MRT-LBM to a two-phase system, which comprised the interaction between a Bingham fluid and a fixed sphere, indicated that the highly viscous region (unyielded region) becomes larger as Bn increases. The qualitative relationship between C S and Bn also confirms the effectiveness of the 3D MRT-LBM. Finally, concerning the drag coefficient on a fixed sphere, it was found that, in contrast to the Newtonian fluid, the drag coefficient is impacted by both Re and Bn.
This work demonstrated some characteristic flow behavior of Bingham plastics but used only a relatively simple scenario for the verification of the 3D parallel MRT-LBM because the current main focus was on the development of the numerical tool.
Simulating realistic flow problems such as selfcompacted concrete flowing through porous material is our ultimate aim because such flows are of principal significance in hydraulic engineering. Therefore, our future work will focus on the 3D fluid-particle interaction for Bingham fluids.