# Velocity driven Poiseuille flow using the incompressible lattice Boltzmann method

My goal here is to develop a Poiseuille flow in a channel, using constant pressure and zero velocity as initial conditions. From what I have read, people seem to be mostly using pressure-driven, or force-driven Poiseuille flows. However, I wanted to use prescribed velocity at inlet (left boundary), and, following what I understood, I am using prescribed pressure at outlet (right boundary), both using Zou-He boundary conditions. Hence, I set the fluid density, $$\rho$$, and the $$Y$$ component of the velocity, $$u_y$$ at the outlet, and computed $$u_x$$ as well as the missing distributions accordingly.

However, what I observe is a higher final velocity at the outlet, as shown below (outlet density set to 1).

Is the free-pressure Zou-He boundary condition not suited in this case? My outlet condition is possibly wrong, although I verified it numerous times. Regarding the rest of the code, it gave convincing results on the driven cavity problem, with good accuracy compared to the reference literature results.

• What does the plot show? The horizontal velocity component or its magnitude? So your initial conditions are $\rho, u_x, u_y = \{1, 0, 0\}$, and your boundary conditions at the inlet something like $u_x(y), u_y \approx \{0.05, 0\}$ where you vary $u_x$ over $y$ to obtain a Poiseuille profile and at the outlet $\rho, u_y = \{1, 0\}$? How precisely do you obtain the simulation viscosity? And finally in which order do you apply the Zou/He boundary conditions and the wall? – 2b-t Mar 16 '20 at 15:44
• In addition to @2b-t's questions, is your channel long enough to fully develop? The Poiseuille flow profile is valid in the infinite axial limit, but obviously it can reach that profile before infinity. If you lengthen your domain, does your solution change? – tpg2114 Mar 16 '20 at 15:54
• @2b-t Thanks for answering. - The plot shows the modulus of velocity - The initial condition is (rho, ux, uy) = (1, 0, 0), and the inlet velocity is set with the exact solution of a Poiseuille. At the outlet, I set (rho, uy) = (1, 0) and deduce ux using classical Zou-He method - I apply the boundary conditions bottom/left/right/top + bottom-left/top-left/bottom-right/top-right corners. The B.C. are applied after macroscopic fields computation and before computing equilibrium - Could you elaborate on the viscosity question ? I set u_lbm = 0.05 for the Mach restriction and deduce the rest – Scrimbibete Mar 16 '20 at 16:14
• @tpg2114 when I lengthen the channel, things get worse. I may be wrong, but as I use the exact solution of the Poiseuille as input, I am not sure its length should have an impact on the solution ? – Scrimbibete Mar 16 '20 at 16:19
• @Scrimbibete I mean that your simulation needs a viscosity and so does your Hagen-Poiseuille profile. How do you determine it? Feel also free to add a link to your implementation, I will have a look at it. – 2b-t Mar 16 '20 at 16:20

Seems as if you incorrectly set the characteristic length to the length of the tube $$L$$ and not its diameter $$D$$. What parameter should be chosen as the characteristic length depends on the problem of interest. It is to some degree arbitrary as long as it reflects the physical behaviour. Clearly the domain length $$L$$ should not influence the solution but instead only the thickness of the tube, characterised either by the diameter $$D = 2 R$$ or the radius $$R$$.

Why did this happen?

For laminar, steady, fully developed flow in a pipe one can find the analytical velocity profile (Hagen-Poiseuille equation)

$$u(r) = \frac{\Delta p (R^2 - r^2)}{4 \mu L} \tag{1}\label{1}$$

with the maximum velocity along the center-line of the tube

$$u_{max} = \frac{\Delta p R^2}{4 \mu L} \tag{2}\label{2}$$

and the average velocity

$$\overline{u} = \frac{\Delta p R^2}{8 \mu L} = \frac{u_{max}}{2}. \tag{3}\label{3}$$

The Reynolds number may be introduced with respect to the average velocity $$\overline{u}$$ and the characteristic length $$D = 2 R$$

$$Re = \frac{2 R \rho \overline{u}}{\mu} = \frac{2 R \overline{u}}{\nu} = \frac{D \overline{u}}{\nu}. \tag{4}\label{4}$$

The velocity profile given by equation \eqref{1} is a parabola that can be also expressed as

$$u(r) = \frac{2 \overline{u}}{R^2} (R^2 - r^2) = \frac{u_{max}}{R^2} (R^2 - r^2): \tag{5}\label{5}$$

The velocity at the center-line $$r=0$$ it will be maximum $$u_{max}$$ and at the walls it will be zero.

As you can now see, if you wrongly choose $$L = c D >> D$$ as the characteristic length while imposing the velocity profile according to equation \eqref{5}, the initial velocity profile (set at the inlet boundary) will not be the one that would actually develop in a flow with the chosen Reynolds number. Instead the Reynolds number of the developed flow will be $$c$$ times higher! As a consequence the velocity profile on the left visually corresponds to a lower Reynolds number than the developed profile on the right.

How to set-up the simulation correctly

Simply use the definition of the Reynolds number with the characteristic diameter $$D_{lb}$$ and not the tube length $$L_{lb}$$ to set the viscosity:

From the definition of the Reynolds number in lattice-Boltzmann ($$lb$$) units

$$Re = \frac{D_{lb} \overline{u}_{lb}}{\nu_{lb}}$$

one can determine the simulation viscosity $$\nu_{lb}$$

$$\nu_{lb} = \frac{D_{lb} \overline{u}_{lb}}{Re}$$

and thus enforce a certain kinetic relaxation time $$\tau_{lb}$$

$$\tau_{lb} = \frac{\nu_{lb}}{c_{S_{lb}}^2} + \frac{\Delta t_{lb}}{2}$$

which for most classic LBM schemes with $$\Delta x_{lb} = \Delta t_{lb} = 1$$ and a lattice speed of sound of $$c_{S_{lb}} = \frac{1}{\sqrt{3}}$$ takes the form

$$\tau_{lb} = \frac{6 \nu_{lb} + 1}{2}.$$

More appropriate: Pressure-drop periodic boundaries

Imposing the developed analytical profile is somehow of limited usage. It is basically just a unit test that the solver works correctly as it should preserve the solution given as the boundary condition. Pressure-drop boundaries on the other hand have the advantage that you can control the simulation just using a pressure drop, the developed profile emerges naturally over time.

If you are interested in implementing pressure-driven flow, I would recommend you either the comparably simple method by Zhang/Kwok (but it also rescales the non-equilibrium distribution) or the more sophisticated approach by Kim/Pitsch. I implemented them myself on $$D2Q9$$- and $$D3Q19$$-lattices and may be able to answer your questions. As far as I can remember, pressure-driven flows where a forcing term is supposed to mimic the pressure drop suffer from artifacts near the walls in combination with simple bounce-back no-slip walls as such a combination does not account for the impact of the forcing term at the wall.

• That is a very detailed answer, thank you a lot for that. I have taken into account your advices on the proper simulation set-up, but I still experience problems. What seems to generate this large velocity on the outlet of the domain is to explicitely set the outlet density to 1. When I don't do it, things get much better (although not perfect). My original intent was to make the final code public, so I will switch it to public now anyway, if you feel like giving it a look, I will be having some free time in the coming weeks anyway... github.com/jviquerat/lbm – Scrimbibete Mar 18 '20 at 10:39
• @Scrimbibete As the entire Italian territory is on quarantine I have lots of time as well the up-coming weeks. ;) I will have a look at your code later today! – 2b-t Mar 18 '20 at 11:02
• That's nice of you, feel free to open issues/request PRs :) – Scrimbibete Mar 18 '20 at 11:09
• @Scrimbibete I just had a quick glance at your code. Is there a reason you used $\overline{u}_{you} = \frac{2}{3} u$ instead of $\overline{u}_{me} = \frac{1}{2} u$ as suggested in my post? If I enforce the pressure at the outlet, use my $\overline{u}_{me}$ and increase the time step then the results look fine to me. Be aware that you need a lot of time steps in order to reach the steady solution as your initial condition is quite unphysical and will create shock wave perturbations that will be reflected by the outlet pressure boundary and thus the first few time steps look very unphysical. – 2b-t Mar 18 '20 at 13:05
• @Scrimbibete If you want to improve this behaviour you could try Latt's regulararised boundary conditions that try to estimate the non-equilibrium distribution $f^{(neq)}$ by $f^{(1)}$ from the Chapman-Enskog expansion in combination with non-equilibrium bounce-back journals.aps.org/pre/abstract/10.1103/PhysRevE.77.056703 or implement the more stable Guo's extrapolation boundaries iopscience.iop.org/article/10.1088/1009-1963/11/4/310/meta. They both improve pressure wave reflection behaviour and stability for higher Reynolds numbers significantly as well. – 2b-t Mar 18 '20 at 13:15