An Improved Multi-Relaxation Time Lattice Boltzmann Method for the Non-Newtonian Influence of the Yielding Fluid Flow in Cement-3D Printing

The multi-relaxation time lattice Boltzmann method (MRT-LBM) has an excellent performance in dealing with the complex flow in many different areas. According to the specific behavior of the fluids, it also has some shortcomings when applied to some special flow like as the non-Newtonian flow. In Cement-3D printing, the fluids always exhibit according to the yielding behavior. When using the standard MRT-LBM, the simulation maybe divergent. In order to solve the problem, this work presents an improved MRT-LBM considering the non-Newtonian effect as a special forcing term to ensure the stable and accurate simulation. Finally, the Poiseuille flow was used to validate the feasibility of the proposed method.


Introduction
3-D printing is an advanced technology to model parts with complex structures [1,2]. It has been widely applied in the mechanical engineering, art, bioengineering, and construction fields [3][4][5]. When it is used in the construction field, the materials are always inorganic mixtures. The fluids are typically non-Newtonian fluids [6][7][8][9]. The flow of the materials in the pipes cannot be measured easily, so simulations are required.
In recent years, lattice Boltzmann method (LBM) has developed into a more popular method for fluid simulation. It can be applied in flow simulation with complex fluid types or various structures [10][11][12]. Although LBM has been widely applied in many areas, it is not stable or accurate in certain cases such as the non-Newtonian flow, fluids with lower or higher Reynolds number, and so on. For Newtonian fluids, the viscosity is a constant, while the viscosity of non-Newtonian fluids will change with the shearing rate. When the viscosity is close to 0 or bigger than 1, the simulation process will be unstable or inaccurate. Power-law fluids are the most common non-Newtonian fluids [13][14][15][16]. Many researchers have modified LBM for the power-law flow to achieve stable and accurate simulations. Niu et al. showed that if the relaxation time can be kept in (1/2, 1), a better stability and higher accuracy will be guaranteed [17]. Based on the above theory, Gabbanelli applied a truncated LBM in a power-law flow, which shows a better simulation [18]. Boyd conducted a local method for the power-law tunnel flow, and the results show a higher accuracy than the truncated method [19]. Nejat et al. applied a second-order lattice Boltzmann method in a power-law flow simulation [20]. Orestis et al. presented a local lattice Boltzmann method for power-law fluids to avoid divergence. It should be noted that the power-index should be kept in a certain range to

Rheological Equation of Yielding Fluids
Common yielding fluids mainly include Bingham fluids and Herschel-Bulkley fluids. They can be expressed using a universal equation: .
where . γ is the shearing rate, τ a is the shearing stress, K is the viscosity coefficient, τ 0 is the initial yielding stress, and n is the power-law index. When n = 1, the fluids will be Bingham fluids; otherwise, it corresponds to Herschel-Bulkley fluids. The kinematic viscosity ν can be calculated using: In 1997, Papanastasiou and Boudouvis proposed a modified equation to simplify the constitutive model [26]: where ρ is the density, and m is used to control the increase of stress. When m approaches zero, the equation will become the constitutive model of power-law fluids.

Improved MRT-LBM
Standard MRT-LBM in two-dimensional space (D2Q9 model) with a forcing term can be expressed as [27]: where r is the displacement vector, e i is the discrete velocity, F is an external force, and δt is the time step size, which is always set as δt = 1. The second item on the right side relates to the external force. f is a distribution function and f eq is an equilibrium distribution function, which is defined as: where c s represents the lattice sound speed, where c 2 s = c 2 /3; u is the velocity vector; ω i is the weight coefficient, given by ω 0 = 4/9, ω i = 1/9 for i = 1-4, and ω i = 1/36 for i = 5-8. M is a transformation matrix: The discrete velocity e i is expressed as: 2c(cos 2i−9 4 π, sin 2i−9 4 π), i = 5, 6, 7, 8 The lattice speed c = δx/δt, δx represents the lattice step size, whose value is always set as 1.
The principal diagonal symmetric matrix S mainly has an important effect on the relaxation process. The specific expression is given by: where s 0 , s 3 , s 5 , s 7 , and s 8 are common relaxation factors, s 3 and s 5 have the same meaning, and s 7 and s 8 have the same meaning. s 1 is related to the viscosity, and s 2 can be used to adjust the stability of numerical simulation. s 4 and s 6 can be used to improve the accuracy of MRT-LBM. The factors have the following relations: s 3 = s 5 , where they have no influence on macroscopic process; s 4 = s 6 , where they affect the accuracy of MRT model; and s 7 = s 8 , where they represent relaxation factors. Then, s 0 is an important related factor to the density, and it also has no effect on the macroscopic process. s 1 is usually related to viscosity. The stability can also be improved by adjusting it. The factor s 2 also has an impact on the stability. Here, s 1 and s 2 are set to be 0.8. According to the above description, s 0 , s 3 , and s 5 are selected arbitrarily as 0. s 4 and s 6 are equal to 1.9. s 7 and s 8 will be discussed further below. With MRT-LBM, kinematic viscosity is defined as: For a standard single relaxation time LBM (LBM mentioned below refers to a standard single relaxation time LBM), the strain rate tensor is described as: where ρ is the density, and τ is the relaxation time, where is τ = 1/s 8 . However, the strain rate tensor in standard MRT-LBM is different from LBM, which is expressed as: Then, the second invariant of strain rate tensor can be acquired using Equation (11) as: With Equation (12), the shearing rate can be calculated as: . γ = 2D I I (13) According to the equality of two different expressions in Reference [21], then substituting viscosity into the related equation, the stress tensor σ αβ can be described as: where δ αβ is a Kronecker delta, K is the viscosity coefficient, and P is the pressure. In the following section, the forcing term will be used to explain the non-Newtonian effect. For standard MRT-LBM with external force, the forcing term is expressed as [23]: Further, the item F in Equation (15) can be calculated using: The Navier-Stokes equation would be expanded by Chapman-Enskog at incompressible limit as: (17) where u α and u β are the components of u, and µ is the dynamic viscosity. The relation between kinematic viscosity and dynamic viscosity is ν = µ/ρ. Also, the above equation can also be described using the following equation, which includes the forcing term: Combining Equations (3) and (14), with Equation (17), the force related term F in Equation (18) can be obtained as: Thus, the non-Newtonian effect described as a special external force will be achieved. Combining the shearing rate obtained from Equation (13) with Equation (3), the parameter s 8 for next iteration could be calculated. The specific programming process based on software MATLAB 2014a (MathWorks, Natick, MA, USA) can be conducted according to the following steps: (1) Conduct the physical transformation according to the dimension theory. The Reynolds number is taken as the key criteria, and then the other parameters used in simulation can be obtained. (2) Set the simulation domain and initial conditions, especially the initial value of the factor s 8 .
(3) Calculate the equilibrium distribution function according to Equation (5). (4) Conduct the collision step, where the specific function is shown as: where (r,t) = Mf(r,t), eq (r,t) = Mf eq (r,t), and the force term is mainly used to explain the non-Newtonian effect of yielding fluids. The calculations can be seen in Equations (15), (16), and (19). (5) Calculate the strain rate tensor and shearing rate by using Equations (11)- (13), then the kinetic viscosity can be obtained using Equation (3) and the relaxation factor s 8 will be updated. (6) Conduct the streaming step, where the streaming function is shown as: (7) Process the boundary condition, where the non-equilibrium bounce-back scheme is taken as the boundary method. (8) Conduct the next calculation from Step (3). (9) Calculate the density and velocity, where the corresponding equations are:

Rheological Equation of Fluids
Many experiments have been done to determine the detailed types and proportions of materials. The specific components are ordinary Portland cement, water, polycarboxylate water reducing agent, modified rod soil, and cellulose ether. The corresponding proportion were 1:0.4:1.5 :0.4 :1.2 . The rheological equation for mixture fluids can be obtained using a viscometer, which is described as: With Equation (23), it can be found that the mixture fluids were Herschel-Bulkley fluids. It belongs to typical yielding non-Newtonian fluids.

Structure of Extruding Head
In the Cement-3D printing experiment, the structure of the printing head is the key technology. According to repeated testing, the single screw extruder was taken as the printing head. The detailed structure is shown as Figure 1a. When expanded along a spiral, a cavity can be obtained, as seen in Figure 1b.
(4) Conduct the collision step, where the specific function is shown as: where (r,t) = Mf(r,t), eq (r,t) = Mf eq (r,t), and the force term is mainly used to explain the non-Newtonian effect of yielding fluids. The calculations can be seen in Equations (15), (16) and (19). (5) Calculate the strain rate tensor and shearing rate by using Equations (11)- (13), then the kinetic viscosity can be obtained using Equation (3) and the relaxation factor s8 will be updated.
(6) Conduct the streaming step, where the streaming function is shown as: (7) Process the boundary condition, where the non-equilibrium bounce-back scheme is taken as the boundary method.
(8) Conduct the next calculation from Step (3). (9) Calculate the density and velocity, where the corresponding equations are:

Rheological Equation of Fluids
Many experiments have been done to determine the detailed types and proportions of materials. The specific components are ordinary Portland cement, water, polycarboxylate water reducing agent, modified rod soil, and cellulose ether. The corresponding proportion were 1:0.4:1.5‰:0.4‰:1.2‰. The rheological equation for mixture fluids can be obtained using a viscometer, which is described as: With Equation (23), it can be found that the mixture fluids were Herschel-Bulkley fluids. It belongs to typical yielding non-Newtonian fluids.

Structure of Extruding Head
In the Cement-3D printing experiment, the structure of the printing head is the key technology. According to repeated testing, the single screw extruder was taken as the printing head. The detailed structure is shown as Figure 1a. When expanded along a spiral, a cavity can be obtained, as seen in Figure 1b. The specific meanings and values of the key parameters have been listed in Table 1.

Flow Simulation using Improved MRT-LBM
The lattice number was set as 480 × 160. The rotation speed of the screw was 40 r/min. The diameter of the screw was designed to be 50 mm. The velocity was only set at the upper surface in Figure 1b and the velocities at other three surfaces were 0. With the above description, flow simulations were conducted. A streamlines figure could be obtained as seen in Figure 2. The simulation of the paste flow using standard MRT-LBM was also conducted, while the result was divergent. It was mainly caused by the non-Newtonian effect. The specific meanings and values of the key parameters have been listed in Table 1.

Flow Simulation using Improved MRT-LBM
The lattice number was set as 480 × 160. The rotation speed of the screw was 40 r/min. The diameter of the screw was designed to be 50 mm. The velocity was only set at the upper surface in Figure 1b and the velocities at other three surfaces were 0. With the above description, flow simulations were conducted. A streamlines figure could be obtained as seen in Figure 2. The simulation of the paste flow using standard MRT-LBM was also conducted, while the result was divergent. It was mainly caused by the non-Newtonian effect. In addition, velocity components were also analyzed. The velocity distributions of the component u at different directions are given in Figure 3. In Figure 3a, it can be found that the biggest velocity u was near the inner wall of the barrel and the smallest velocity u was near the outer wall of the screw. In Figure 3b, it shows that the biggest velocity u was near the center line of the channel and the smallest velocity u was at the wall of swirling leaves. In addition, velocity components were also analyzed. The velocity distributions of the component u at different directions are given in Figure 3. In Figure 3a, it can be found that the biggest velocity u was near the inner wall of the barrel and the smallest velocity u was near the outer wall of the screw. In Figure 3b, it shows that the biggest velocity u was near the center line of the channel and the smallest velocity u was at the wall of swirling leaves. The specific meanings and values of the key parameters have been listed in Table 1.

Flow Simulation using Improved MRT-LBM
The lattice number was set as 480 × 160. The rotation speed of the screw was 40 r/min. The diameter of the screw was designed to be 50 mm. The velocity was only set at the upper surface in Figure 1b and the velocities at other three surfaces were 0. With the above description, flow simulations were conducted. A streamlines figure could be obtained as seen in Figure 2. The simulation of the paste flow using standard MRT-LBM was also conducted, while the result was divergent. It was mainly caused by the non-Newtonian effect. In addition, velocity components were also analyzed. The velocity distributions of the component u at different directions are given in Figure 3. In Figure 3a, it can be found that the biggest velocity u was near the inner wall of the barrel and the smallest velocity u was near the outer wall of the screw. In Figure 3b, it shows that the biggest velocity u was near the center line of the channel and the smallest velocity u was at the wall of swirling leaves. It can be found from Figure 4a that the biggest velocity v was at the midway between the barrel and the screw. Then the velocity v was always 0 along the center line of channel. In Figure 4b, the velocity v near the outer wall of the screw was approximately 0 and it fluctuated near the wall of the swirling leaves. It was bigger when the position was closer to the inner wall of the barrel.  It can be found from Figure 4a that the biggest velocity v was at the midway between the barrel and the screw. Then the velocity v was always 0 along the center line of channel. In Figure 4b, the velocity v near the outer wall of the screw was approximately 0 and it fluctuated near the wall of the swirling leaves. It was bigger when the position was closer to the inner wall of the barrel.

Validation and Discussion
The obtained results demonstrate the stability of the new, improved MRT-LBM. To ensure accuracy, we took the theoretical results of Poiseuille flow to compare with numerical results of the proposed method. The theoretical solution of yielding fluids is shown as: where H represents the distance between the two parallel plates, and yτ is a critical point, which is equal to: We consider that the pivotal factors mainly included initial yielding stress, power-index, and lattice number. The accuracy will be discussed using the result of the relative error. The equation of the error er is: The power index n was set as 0.5, 1, 2, where these values correspond to shear-thinning fluids, Bingham fluids, and shear-thickening fluids, respectively. For different power indexes, the effect of lattice numbers and initial yielding stresses were explored. To avoid the effect of other factors, we set the pressure gradient / P x ∂ ∂ = −1.0 × 10 −5 , viscosity coefficient 0 μ = 0.01, density ρ = 1, and distance H = 1. Lattice numbers were set in the range of 50-250 and the initial yielding stresses were selected as 1.0 × 10 −6 , 2.5 × 10 −6 , and 3.5 × 10 −6 .
In Figure 5, it shows the numerical velocities and theoretical velocities for different initial yielding stresses when the power index n = 0.5 and lattice number was 200. It can be seen that the numerical velocities were consistent with the theoretical velocities at different initial yielding

Validation and Discussion
The obtained results demonstrate the stability of the new, improved MRT-LBM. To ensure accuracy, we took the theoretical results of Poiseuille flow to compare with numerical results of the proposed method. The theoretical solution of yielding fluids is shown as: where H represents the distance between the two parallel plates, and y τ is a critical point, which is equal to: We consider that the pivotal factors mainly included initial yielding stress, power-index, and lattice number. The accuracy will be discussed using the result of the relative error. The equation of the error e r is: The power index n was set as 0.5, 1, 2, where these values correspond to shear-thinning fluids, Bingham fluids, and shear-thickening fluids, respectively. For different power indexes, the effect of lattice numbers and initial yielding stresses were explored. To avoid the effect of other factors, we set the pressure gradient ∂P/∂x = −1.0 × 10 −5 , viscosity coefficient µ 0 = 0.01, density ρ = 1, and distance H = 1. Lattice numbers were set in the range of 50-250 and the initial yielding stresses were selected as 1.0 × 10 −6 , 2.5 × 10 −6 , and 3.5 × 10 −6 .
In Figure 5, it shows the numerical velocities and theoretical velocities for different initial yielding stresses when the power index n = 0.5 and lattice number was 200. It can be seen that the numerical velocities were consistent with the theoretical velocities at different initial yielding stresses. In Figure 6, it gives the relative error for each case of three initial yielding stresses with different lattice numbers. The relative error decreased with increasing lattice number. stresses. In Figure 6, it gives the relative error for each case of three initial yielding stresses with different lattice numbers. The relative error decreased with increasing lattice number.  When the power index n = 1, it corresponds to a special kind of yielding fluids named Bingham fluids, and the velocities of the simulation and theory at three different initial yielding stresses are shown in Figure 7. It shows that the numerical results also had a good consistency with the theoretical results for all cases. The figure of relative errors are given in Figure 8 for three situations. It suggests that the bigger the lattice number, the smaller the relative error will be. stresses. In Figure 6, it gives the relative error for each case of three initial yielding stresses with different lattice numbers. The relative error decreased with increasing lattice number.  When the power index n = 1, it corresponds to a special kind of yielding fluids named Bingham fluids, and the velocities of the simulation and theory at three different initial yielding stresses are shown in Figure 7. It shows that the numerical results also had a good consistency with the theoretical results for all cases. The figure of relative errors are given in Figure 8 for three situations. It suggests that the bigger the lattice number, the smaller the relative error will be. When the power index n = 1, it corresponds to a special kind of yielding fluids named Bingham fluids, and the velocities of the simulation and theory at three different initial yielding stresses are shown in Figure 7. It shows that the numerical results also had a good consistency with the theoretical results for all cases. The figure of relative errors are given in Figure 8 for three situations. It suggests that the bigger the lattice number, the smaller the relative error will be. Materials 2018, 11, x FOR PEER REVIEW 9 of 12  When the power index n = 2, it belongs to shear-thickening fluids. The velocities of simulation and theory at three different initial yielding stresses are shown in Figure 9. In Figure 10, it gives the tendency between relative error and lattice number. Both the good consistencies and small relative errors demonstrate that the proposed method is effective.  When the power index n = 2, it belongs to shear-thickening fluids. The velocities of simulation and theory at three different initial yielding stresses are shown in Figure 9. In Figure 10, it gives the tendency between relative error and lattice number. Both the good consistencies and small relative errors demonstrate that the proposed method is effective. When the power index n = 2, it belongs to shear-thickening fluids. The velocities of simulation and theory at three different initial yielding stresses are shown in Figure 9. In Figure 10, it gives the tendency between relative error and lattice number. Both the good consistencies and small relative errors demonstrate that the proposed method is effective.
In Figures 5, 7, and 9, the only difference of conditions is the power index. It was found that the power index had an important effect on velocity. With the increasing of the power index, the velocity increased rapidly. Also, the initial yielding stress affected the velocity significantly. The smaller initial yielding stress led to a bigger velocity. For different initial yielding stresses and power indexes, all results showed a good consistency. This indicates that the proposed method is stable and valid.
In Figures 6, 8, and 10, for different power indexes and initial yielding stresses, the similarity was that the relative error decreased with increasing lattice number. The slope of each line was very close to −1. The smaller errors indicated that the improved method is accurate.  In Figures 5, 7, and 9, the only difference of conditions is the power index. It was found that the power index had an important effect on velocity. With the increasing of the power index, the velocity increased rapidly. Also, the initial yielding stress affected the velocity significantly. The smaller initial yielding stress led to a bigger velocity. For different initial yielding stresses and power indexes, all results showed a good consistency. This indicates that the proposed method is stable and valid.
In Figures 6, 8, and 10, for different power indexes and initial yielding stresses, the similarity was that the relative error decreased with increasing lattice number. The slope of each line was very close to −1. The smaller errors indicated that the improved method is accurate.

Conclusions
The non-Newtonian effect will result in a poor stability and accuracy when using standard LBM in the numerical simulation of yielding fluids. To improve the situation, a modified multi-relaxation-time lattice Boltzmann method was proposed. The non-Newtonian part was modeled with a special forcing term. Then the improved method was applied in a simulation of mixed paste in the extruder of Cement-3D printing. Stable streamlines and velocity distributions  In Figures 5, 7, and 9, the only difference of conditions is the power index. It was found that the power index had an important effect on velocity. With the increasing of the power index, the velocity increased rapidly. Also, the initial yielding stress affected the velocity significantly. The smaller initial yielding stress led to a bigger velocity. For different initial yielding stresses and power indexes, all results showed a good consistency. This indicates that the proposed method is stable and valid.
In Figures 6, 8, and 10, for different power indexes and initial yielding stresses, the similarity was that the relative error decreased with increasing lattice number. The slope of each line was very close to −1. The smaller errors indicated that the improved method is accurate.

Conclusions
The non-Newtonian effect will result in a poor stability and accuracy when using standard LBM in the numerical simulation of yielding fluids. To improve the situation, a modified multi-relaxation-time lattice Boltzmann method was proposed. The non-Newtonian part was modeled with a special forcing term. Then the improved method was applied in a simulation of mixed paste in the extruder of Cement-3D printing. Stable streamlines and velocity distributions

Conclusions
The non-Newtonian effect will result in a poor stability and accuracy when using standard LBM in the numerical simulation of yielding fluids. To improve the situation, a modified multi-relaxation-time lattice Boltzmann method was proposed. The non-Newtonian part was modeled with a special forcing term. Then the improved method was applied in a simulation of mixed paste in the extruder of Cement-3D printing. Stable streamlines and velocity distributions were obtained. For the velocity component u, the biggest velocity in two directions occurred at the position of the inner wall of the barrel and the center of the channel, respectively. The smallest velocity in two directions occurs at the position of the outer wall of screw and the wall of swirling leaves, respectively. For the velocity component v, it was 0 when the position was at the center of the channel and the outer wall of the screw. In order to prove the stability and accuracy of the method, numerical solutions of the Poiseuille flow were compared with theoretical solutions. The numerical results were consistent with the theoretical results. For the cases of different initial yielding stresses and power-indexes, the relative errors were small enough to verify the feasibility of the improved method.