Poiseuille 2D Questions

Palabos General
#1

Hello!

I’ve got some problem with understanding something with Poiseuille 2d flow. I’m using velocity BC on the upper and lower sides of the box, and pressure BC on the inlet and outlet. I don’t understand why when I initialize the system I need to setup velocity based on poiseuilleVelocity function and can’t use some fixed velocity for all particles and let dynamics to stabilize flow like in forcedbgk example? I tried this and for RE number around 10-40 or maybe more u[sub]max[/sub] is ok - it’s converging to about 0.99…, but when I start simulation with RE = 100 or more u[sub]max[/sub] keep on rising over 1.0

Sorry for my poor english.
Thanks in advance for answer.

#2

Hi,

This is about numerical stability. In principle, you can use just any initial condition for this problem, because the channel flow is stationary, at least at modest Reynolds numbers: you’re not interested in the time evolution of the flow. Instead, you’re looking at the steady state solution which is reached after a sufficient number of iterations, no matter what initial condition is used.

However, the LB scheme can become numerically unstable if you violate too many physical principles. If you initialize the velocity to a constant, non-zero value everywhere in the initial condition, the velocity drops from constant to zero from one lattice node to another close to the boundary. This means that (1) the velocity field is discontinuous and (2) the divergence of the velocity is far from being zero.

If you want to avoid cheating by initializing everything to the Poiseuille solution, you can try the following: initialize the velocity to zero everywhere to avoid velocity discontinuities, and the density to a linearly decreasing profile to avoid density discontinuities. This setup should be stable and converge (slowly) to the Poiseuille solution.

#3

I previously tried zero velocity, but i had also density set to 1.

If I understand you correctly i should modify rho in initializing part something like this:
T rho = (T)1 - iX * 0.000001

If i understand correctly, what should be density decreasing rate? It should be computed from the Nx size of the box?

#4

Well, you do have a boundary condition which imposes the value of density at the inlet and at the outlet, right? So, for in-between lattice nodes, simply interpolate between these two values.

#5

Sorry, but I still can’t get it right.

So main part of my function initializing boundries looks like this:


lattice.defineDynamics(0,nx-1, 0,ny-1, &bulkDynamics);
boundaryCondition.addVelocityBoundary1P(1,nx-2,ny-1,ny-1, omega);
boundaryCondition.addVelocityBoundary1N(1,nx-2,   0,   0, omega);
boundaryCondition.addPressureBoundary0N(0,0, 1, ny-2, omega);
boundaryCondition.addPressureBoundary0P(nx-1,nx-1, 1, ny-2, omega);
boundaryCondition.addExternalVelocityCornerNN(0,0, omega);
boundaryCondition.addExternalVelocityCornerNP(0,ny-1, omega);
boundaryCondition.addExternalVelocityCornerPN(nx-1,0, omega);
boundaryCondition.addExternalVelocityCornerPP(nx-1,ny-1, omega);

//density for inlet pressure BC
T rho1 = 1.00001;
//density for outlet pressure BC
T rho2 = 1.0;
T delta_rho = (rho1 - rho2)/nx;
for (int iX=0; iX<nx; ++iX) {
        for (int iY=0; iY<ny; ++iY) {
             //linear decrease of density from inlet to outlet
            T rho = rho1 - iX*delta_rho;
            T u[2] = {0. ,0.};
            lattice.get(iX,iY).defineRhoU(rho, u);
            lattice.get(iX,iY).iniEquilibrium(rho,u);
        }
    }

for (int iY=1; iY<ny-1; ++iY) {
    lattice.get(0,iY).defineRho(rho1);
    }

for (int iY=1; iY<ny-1; ++iY) {
        lattice.get(nx-1,iY).defineRho(rho2);
        }
lattice.initialize();

Results are still very strange, computed u[sub]max[/sub] is different whenever i change the densities on the inlet/outletand the differences are really big. Where I’m doing it wrong?

#6

Spontaneously, I would say that your code looks right. What exactly are you concerned about? It is an expected behavior that uMax depends on the pressure difference. If you push the flow stronger, it flows faster.

#7

I think i need to read more about units conversion to make some sense from these numbers.

For example to the LBunits object i provide: u[sub]max[/sub] = 0.02, Re = 10, lx = 2. , ly = 1. and N = 50 and set my pressure BC to rho1 = 1.001(inlet) and rho2=1.0(outlet).

After 40000 iterations I get steady flow with center velocity u[sub]max[/sub] = 0.52027787, and I don’t know if it is too much, or too slow, or ok. (u[sub]max[/sub] = middleU*dx/dt) Velocity profile generetad in gifs looks ok, also ploted in gnuplot. U[sub]max[/sub] = 3/2u[sub]av[/sub] so maybe everything is ok?

How can I compare my results from pressure gradient simulation with body force simulation - how to calculate body force that i should impose on the system?

#8

You can’t just put the pressure difference to an arbitrary value. In the analytical solution for a Poiseuille flow, the maximum velocity is directly related to the pressure gradient (or to the body force, in a force-driven flow). This formula is quite easy to derive (check out for example this link.

#9

Thank you for your patience in explaining this to me.

If anyone else encountered similar problem here is solution:

  1. Boundry conditions 4 posts up are setup ok.
  2. Pressure is calculated from equation number 3.16 from this page

p - p_0 = c_s^2 \rho

So to get Rho we just multiply pressure by the invCs2 provided by Descriptor.


T rho1 = 1. + (converter.getLatticeU() * 2. * converter.getLatticeNu() * (converter.getNx()-1))/((converter.getNy()-1)/2.)/((converter.getNy()-1)/2) * DESCRIPTOR<T>::invCs2;
T rho2 = 1.0;

Please correct me if I’m wrong, but I think now results are ok.

Thank you Jonas for pointing me to right direction, I was on the page you mentioned before but I wasn’t sure what i should be looking for.