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:
- I use full bounceback BCs at the walls, since currently there is no other option in Palabos for irregular walls
- 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)
- 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.
- 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_),
radiusSqr(util::sqr(radius))
{ }
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;
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{
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;
};