3D Poiseuille Flow with bounce-back walls

Hi,

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!

• 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 ?

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.

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));

lattice.initialize();
}

``````

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 {
public:
CylinderShapeDomain3D(T cy_, T cz_, T radius)
: cy(cy_),
cz(cz_),
{ }
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);
}
private:
T cy;
T cz;
};

``````
``````
/// 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{
public:
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);
}

private:
IncomprFlowParam<T> parameters;
};

``````

Hi,