3D Poiseuille Flow with bounce-back walls


I am trying to solve Poiseuille flow in a 3D tube. Since there is no other possibility, I have implemented a staircase-approximation of the curved tube walls made of bounce-back nodes.
For my project, I need to use velocity boundary conditions, so at the inlet I have a local Poiseuille velocity profile, at the outlet I am using normal outflow conditions. The corner nodes have bounce back dynamics.

So far, for small Reynolds numbers (10-100) everything works fine. But for Re=1000, I need a very high resolution (ny>>100 for the tube diameter) to obtain a proper result, otherwise the simulation is unstable.

For example:

Re=1000, tau =0.6, umax = 0.2 gives ny = Re * (2tau-1)/6 / umax = 167 works (but inaccurate due to large Mach number)
Re=1000, tau =0.56, umax = 0.2 gives ny = Re * (2
tau-1)/6 / umax = 100 is instable (gives NaN after some iterations)

Is there no other possibility to get a stable simulation than increasing resolution? (in fact, for D3Q19 a resolution of nz=ny=167, nx=501, already yields 501167167198 byte = 2.12 GB of memory consumption!)

Moreover, with my current boundary setting, Zou/He or Skordos conditions at inlet/outlet do not work at all (segmentation fault). Is this because of my mixed boundary conditions? (bounce back/velocity)

Thanks for help!

A few comments:

  • Is your simulation turbulent or not? If it is, this means that you are running a direct numerical simulation of fluid turbulence. It is not surprising then that you need a high resolution, because you need to fully resolve all structures in the flow (basically, the size of a lattice cell must me smaller or equal to the Kolmogorov length).
  • Isn’t your Mach number a bit high? Your simulation may be more stable if you use a smaller velocity, as measured in lattice units.
  • The combination of inlet vs. bounce-back nodes currently only works with the regularized boundary condition (“createLocalBoudnaryCondition”). We’re working on extending this to other boundary conditions.

Thank you, jlatt. Currently, I have just implemented a straight tube with Re=1000, referring to diameter and maximal velocity. So the flow should be laminar. In most cases, the simulation first converges towards the Poiseuille solution and then suddenly blows up in a few time steps.

I tried also the outflow-condition instead of normalOutflow. For some parameters, this seems to be stable, but for others it gives same behaviour. I still have to look closer into this.
Has anyone similar experience with high Reynolds numbers in complex (non-straight) geometries?

A related question, posed more generally: Where can I find general information on the stability range of the Lattice Boltzmann method, in 3D? For instance, which is the largest Reynolds number that can be solved stable with a given resolution?

I found the paper Latt&Chopard 2005, which gives stability ranges for cavity flow in 2d. From figure 3 in that paper, I read that for Re=1000, I probably need a resolution of N>150. But is this also valid for Poiseuille flow in 3d?

Hi smeier,

I would like to implement 3D Poiseuille Flow, in a tube.
Can you help me by sharing your experience/code ?

Thank you in advance !

Hi tekicp,

here is my channel setup routine (analagously to the 2d tutorials in Palabos)

My tube is in x-direction. At the inlet, I specify velocity Dirichlet conditions, at the outlet normalOutflow.

A few comments:

  1. I use full bounceback BCs at the walls, since currently there is no other option in Palabos for irregular walls
  2. I use local (velocity) boundary conditions at inlet and outlet. The alternatives are ZouHe or Interp, but they don’t work in combination with bounce-back walls. (at least with my current version of Palabos)
  3. For low resolutions, it makes a large difference where you actually place your wall. I figured out that it is good to have the domain in y- and z-direction one half grid cell larger than the actual pipe. If the pipe has diameter ly=lz in physical units, then one can achieve this with the converter class IncompFlowParameters by calling ny = getNy(true) instead of ny= getNy().
    Finally, I define the wall by the DomainFuntional CylinderShapeDomain with radius (ny-1)/2-0.38 in Lattice units. I played a little bit with this parameter, and found 0.38 to give quite accurate velocity results.
  4. On the choice of parameters … : In my case, I need quite large Reynolds numbers, so I tried to find a good setting for Re=1000. But this seems to work only for high resolutions N>100 (which means a grid of 100x100x300 points if lx=3, ly=lz=1). For lower resolutions, the boundary conditions cause instabilities.

Moreover, I have problems with oscillating pressure waves, even for smaller Reynolds number. Hence, the velocity is very accurate, but the pressure is completely wrong. Since I am a beginner with LBM, I haven’t yet figured out whether this a typical problem or not.

For small Re (Re=10, say) the code works fine and the velocity and pressure is pretty accurate.

It would be nice to hear from your experience!

void channelSetup( BlockLatticeBase3D<T,DESCRIPTOR>& lattice,
            IncomprFlowParam<T>& parameters,
                   OnLatticeBoundaryCondition3D<T,DESCRIPTOR>& boundaryCondition )

    // define boundary parts
    plint nx = parameters.getNx();
    plint ny = parameters.getNy(true); 
    plint nz = parameters.getNz(true); 
    T dt = parameters.getDeltaT();
    T dx = parameters.getDeltaX();

    Box3D inlet(0, 0, 1, ny-2, 1, nz-2);
    Box3D outlet(nx-1, nx-1, 1, ny-2, 1, nz-2);

    // set inlet boundary condition
    boundaryCondition.setVelocityConditionOnBlockBoundaries(lattice, inlet, boundary::dirichlet);

    // set outlet boundary condition
    boundaryCondition.setVelocityConditionOnBlockBoundaries(lattice, outlet, boundary::normalOutflow);

    // Specify the boundary velocity at inlet
    setBoundaryVelocity (
            lattice, inlet,
            PoiseuilleVelocity3D<T,DESCRIPTOR>(parameters) );

    // set bounce back conditions outside of a cylinder
    defineDynamics(lattice, lattice.getBoundingBox(),
                   new CylinderShapeDomain3D<T>((T)(ny-1)/2, (T)(nz-1)/2, (T)(ny-1)/2-0.38),
                   new BounceBack<T,DESCRIPTOR>);

    // set NoDynamics on larger cylinder
    // for D3Q19, the distance between the radii should be greater than sqrt(2)
    defineDynamics(lattice, lattice.getBoundingBox(),
                   new CylinderShapeDomain3D<T>((T)(ny-1)/2, (T)(nz-1)/2, (T)(ny-1)/2+1.414),
                   new NoDynamics<T,DESCRIPTOR>);

    // Create the initial condition
    initializeAtEquilibrium (lattice, lattice.getBoundingBox(), PoiseuilleVelocity3D<T,DESCRIPTOR>(parameters));


The routine makes use of the classes CylinderShapeDomain3D and PoiseuilleVelocity3D, which are given below.

// A class for defining specific dynamics outside a given cylinder
// defines a cylinder in x-direction, centered at cy, cz (lattice coordinates)
template<typename T>
class CylinderShapeDomain3D : public DomainFunctional3D {
    CylinderShapeDomain3D(T cy_, T cz_, T radius)
        : cy(cy_),
    { }
    virtual bool operator() (plint iX, plint iY, plint iZ) const {
        return util::sqr(iY-cy) + util::sqr(iZ-cz) >= radiusSqr;
    virtual CylinderShapeDomain3D<T>* clone() const {
        return new CylinderShapeDomain3D<T>(*this);
    T cy;
    T cz;
    T radiusSqr;

/// Velocity on the parabolic Poiseuille profile. (in Lattice units)
T poiseuilleVelocity(plint iY, plint iZ, IncomprFlowParam<T> const& parameters) {

        plint ny = parameters.getNy(true);    
        plint nz = parameters.getNz(true);  

        T r_sq = (iY-(T)(ny-1)/2)*(iY-(T)(ny-1)/2) + (iZ-(T)(nz-1)/2)*(iZ-(T)(nz-1)/2);
        T R_sq = (T)(ny-2)*(T)(ny-2)/4; 
        T y_sq = r_sq /R_sq;
        if (y_sq <= 1) {
                return parameters.getLatticeU() * (1 - y_sq);
        } else {
                return 0;

//  used to initialize the velocity for the boundary conditions.
template<typename T, template<typename U> class Descriptor>
class PoiseuilleVelocity3D{
    PoiseuilleVelocity3D(IncomprFlowParam<T> parameters_)
        : parameters(parameters_)
    { }

    // returns the poiseuille velocity; for creating boundary conditions
    void operator()(plint iX, plint iY, plint iZ, Array<T,3>& u) const {
        u[0] = poiseuilleVelocity(iY, iZ, parameters);
        u[1] = T();
        u[2] = T();
    /// This version of the operator returns also a linear decreasing profile for
    ///    the density, to create the initial condition.
    void operator()(plint iX, plint iY, plint iZ, T& rho, Array<T,3>& u) const {
        u[0] = poiseuilleVelocity(iY, iZ, parameters);
        u[1] = T();
        u[2] = T();
        plint nx = parameters.getNx();
        T dt = parameters.getDeltaT();
        T dx = parameters.getDeltaX();
        rho  = 1.0 - DESCRIPTOR<T>::invCs2 * (dt*dt)/(dx*dx)* 48.0/parameters.getRe() * (iX-(nx-1)/2)/(nx-1);

    IncomprFlowParam<T> parameters;


Two quick comments:

  1. It’s true that the Interpolation (Skordos) boundary condition does not work when it touches bounce-back nodes. However, Zou/He should be fine if you use the most recent release (a bug in the Zou/He boundary condition has been corrected).
  2. Pressure waves which bounce back and forth can be very slow to damp out. It’s likely that in your case these waves are due to the initial condition. In this case, you could possibly decrease their effect if the initial condition has a proper density profile, decreasing linearly from inlet to outlet.


I tried the ZouHe boundary conditions with the new Palabos release 0.6.1, and in fact, my first results are much more stable: With the regularizedBGK model and ZouHe velocity conditions and bounce-back walls I have stable and accurate simulations (relative velocity error < 0.005) for Re=1000 with only res=40, tau=0.506, vmax=0.05.

So, many thanks for the update!