# Sudden Expansion flow

I am a new user of palabos.
I try to solve sudden Expansion flow problem.
In order to solve it I took the tutorial code of Poiseuille flow and changed the geometry.
However when running the simulation it seems that the walls are not solid and I still recieve a tube tunnel in the result.
The code is attached and I will be greatfull for any asistance.
Thank you.

[code=“cpp”]
#include “palabos2D.h”
#include “palabos2D.hh”
#include
#include
#include

/* Code 1.5 in the Palabos tutorial
*/

using namespace plb;
using namespace std;

typedef double T;
#define DESCRIPTOR plb::descriptors::D2Q9Descriptor

const int nx = 1200;
const int ny = 400;

/// Velocity on the parabolic Poiseuille profile
T poiseuilleVelocity(plint iY, IncomprFlowParam const& parameters) {
T y = (T)iY / parameters.getResolution();
return 4.parameters.getLatticeU() * (y-yy);
}

/// A functional, used to initialize the velocity for the boundary conditions
template
class PoiseuilleVelocity {
public:
PoiseuilleVelocity(IncomprFlowParam parameters_)
: parameters(parameters_)
{ }
/// This version of the operator returns the velocity only,
/// to instantiate the boundary condition.
void operator()(plint iX, plint iY, Array<T,2>& u) const {
u[0] = poiseuilleVelocity(iY, parameters);
u[1] = T();
}
/// This version of the operator returns also a constant value for
/// the density, to create the initial condition.
void operator()(plint iX, plint iY, T& rho, Array<T,2>& u) const {
u[0] = poiseuilleVelocity(iY, parameters);
u[1] = T();
rho = (T)1;
}
private:
IncomprFlowParam parameters;
};

void channelSetup (
MultiBlockLattice2D<T,DESCRIPTOR>& lattice,
IncomprFlowParam const& parameters,
OnLatticeBoundaryCondition2D<T,DESCRIPTOR>& boundaryCondition )
{

``````Box2D topLow(0, nx/2-1, ny*5/8-1, ny*5/8-1);
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, topLow, boundary::freeslip );
Box2D topSide(nx/2-1, nx/2-1, ny*5/8-1 , ny-1);
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, topSide, boundary::freeslip );
Box2D topWall(nx/2-1, nx-1,  ny-1, ny-1);
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, topWall, boundary::freeslip );
Box2D bottomLow(0, nx/2-1, ny*3/8-1, ny*3/8-1);
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, bottomLow, boundary::freeslip );
Box2D bottomSide(nx/2-1, nx/2-1, ny*3/8-1, 0);
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, bottomSide, boundary::freeslip );
Box2D bottomWall(nx/2-1, nx-1, 0, 0);
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, bottomWall, boundary::freeslip );
Box2D inlet(0, 0, ny*3/8+2, ny*5/8+2);
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, inlet );
Box2D outlet(nx-1, nx-1, 2, ny-2);
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, outlet); //, boundary::outflow );
boundaryCondition.addInternalVelocityCornerPN( nx/2, ny*5/8, lattice, boundary::freeslip );
boundaryCondition.addInternalVelocityCornerPP( nx/2, ny*3/8, lattice, boundary::freeslip );
boundaryCondition.addExternalVelocityCornerPN( nx/2, ny, lattice, boundary::freeslip );
boundaryCondition.addExternalVelocityCornerPP( nx/2, 0, lattice, boundary::freeslip );

// Specify the boundary velocity.
setBoundaryVelocity (
lattice, lattice.getBoundingBox(),
PoiseuilleVelocity<T>(parameters) );

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

lattice.initialize();
``````

}

void writeGifs(MultiBlockLattice2D<T,DESCRIPTOR>& lattice, plint iter)
{
const plint imSize = 600;
ImageWriter imageWriter(“leeloo”);
imageWriter.writeScaledGif(createFileName(“u”, iter, 6),
*computeVelocityNorm(lattice),
imSize, imSize );
}

int main(int argc, char* argv[]) {
plbInit(&argc, &argv);

``````global::directories().setOutputDir("./tmp/");

// Use the class IncomprFlowParam to convert from
//   dimensionless variables to lattice units, in the
//   context of incompressible flows.
IncomprFlowParam<T> parameters(
(T) 1e-1,  // Reference velocity (the maximum velocity
//   in the Poiseuille profile) in lattice units.
(T) 100.,  // Reynolds number
300,       // Resolution of the reference length (channel height).
3.,        // Channel length in dimensionless variables
1.         // Channel height in dimensionless variables
);
const T imSave   = (T)0.1;  // Time intervals at which to save GIF
//   images, in dimensionless time units.
const T maxT     = (T)3.1;  // Total simulation time, in dimensionless
//   time units.

writeLogFile(parameters, "Poiseuille flow");

MultiBlockLattice2D<T, DESCRIPTOR> lattice (
parameters.getNx(), parameters.getNy(),
new BGKdynamics<T,DESCRIPTOR>(parameters.getOmega()) );

OnLatticeBoundaryCondition2D<T,DESCRIPTOR>*
boundaryCondition = createLocalBoundaryCondition2D<T,DESCRIPTOR>();

channelSetup(lattice, parameters, *boundaryCondition);

// Main loop over time iterations.
for (plint iT=0; iT*parameters.getDeltaT()<maxT; ++iT) {
if (iT%parameters.nStep(imSave)==0 && iT>0) {
pcout << "Saving Gif at time step " << iT << endl;
writeGifs(lattice, iT);
}
// Execute lattice Boltzmann iteration.
lattice.collideAndStream();
}

delete boundaryCondition;
``````

}

I just want to change the geometry from box to those shown in the picture.

And I don’t understand why I dont get “solid walls” with the code that I wrote.

http://www.palabos.org/documentation/userguide/_images/boundaries.gif

[code=“cpp”]
// Velocity boundary condition on bottom wall.
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, Box2D(1, nx/8-1, ny3/8, ny3/8) );
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, Box2D(nx/8-1, nx/8-1, ny3/8, 1) );
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, Box2D(nx/8-1, nx-1, 1, 1) );
// Velocity boundary condition on top wall.
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, Box2D(1, nx/8-1, ny
5/8, ny5/8) );
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, Box2D(nx/8-1, nx/8-1, ny
5/8, ny-1) );
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, Box2D(nx/8-1, nx-1, ny-1, ny-1) );

``````boundaryCondition.addInternalVelocityCornerPN( nx/8-1, ny*3/8, lattice );