Skip to content

Permeability Calculation Poiseuille

Viewing 10 posts - 1 through 10 (of 10 total)
  • Author
  • #4243

    i`m quite new to openLB and try to implement the calculatiion of the permeability for porous media. Therefore i first want to look at the simple case of the Poiseuille2d Code, since the analytical solution for verification is known.

    The Problem is the following:
    I used the original poiseuille2d Code and only replaced the boundary condition on the front by a pressure boundary condition, so that i can directly use the pressure Gradient as (pin-pout)/lx.
    The mean velocity is also calculated. But if i do just these changes the result doesn`t coincide with the analytical solution…

    The changes in the Code are the following:

    – Parameters for the Simulation:
    const T pin=1000;
    const T pout=500;
    const T density=1e3;
    const T viscosity=1e-6;

    – void prepareLattice()
    sBoundaryCondition.addPressureBoundary(superGeometry, 3, Omega);
    Setting of the inlet pressure etc.

    Unfortunately i am not even sure if the forced-FlowType is correct in this case since it may spoil the boundary condition. Any suggestions how to write down a simple Poiseuille2d Code with pressure boundary conditions just to verify Darcys law. I am quite lost since a couple of days now…

    Thanks in regards,


    Dear Ladasicr,

    the problem may be due to bugs in the GuoZhaoBGKDynamics. These problems are now fixed and will be released soon. Until then you could try the PorousBGKDynamics based on the porous media model of Spaid and Phelan.

    Best regards,


    Dear Fabian,

    thanks for the Reply. Unfortunately i am not that deep into the Topic so far… My Problem may rather sound trivial but it is the following:

    I use pressure BC for the inlet in the poiseuille2d example and after the main Loop i just evaluate the mean velocity in order to calculate the Dary law, e.g. the permeability.
    The Problem occurs now if i set FlowType=unForced, DESCRIPTOR = D2Q9Descriptor and apply the inlet pressure BC as hard Parameters, ie. 100 Pa at the inlet and 10 Pa at the outlet.
    Then the Simulation just breaks because with the given Parameters i get for all Parameters like uMax=inf, avEngergy=-nan, avRho=nan etc.

    Any idea how to fix that simple example?

    Best regards,


    Dear Ladasicr,

    could you post the code of how you set the boundary conditions here? My guess is, that you do not have a smooth start up and this is why the simulation blows up.



    Dear Fabian,

    here is the relevant Code which i included in the Pousoille2d example:

    // In function prepareLattice

    sBoundaryCondition.addPressureBoundary(superGeometry, 3, omega);
    sBoundaryCondition.addPressureBoundary( superGeometry, 4, omega );

    std::vector<T> velocity( 2,T( 0 ) );
    AnalyticalConst2D<T,T> u( velocity );
    AnalyticalConst2D<T,T> rho( 1 );

    // Define inlet and outlet pressure
    AnalyticalConst2D<T,T> rhoIn(converter.getLatticeDensityFromPhysPressure(pin)); // pin in Pa
    AnalyticalConst2D<T,T> rhoOut(converter.getLatticeDensityFromPhysPressure(pout));

    // Initialize all values of distribution functions to their local equilibrium
    sLattice.defineRho( superGeometry, 1, rho);
    sLattice.iniEquilibrium( superGeometry, 1, rho, u );
    sLattice.defineRho( superGeometry, 2, rho);
    sLattice.iniEquilibrium( superGeometry, 2, rho, u );
    sLattice.defineRho( superGeometry, 3, rhoIn);
    sLattice.iniEquilibrium( superGeometry, 3, rho, u );
    sLattice.defineRho( superGeometry, 4, rhoOut);
    sLattice.iniEquilibrium( superGeometry, 4, rho, u );

    // In setBoundaryValues
    if ( iT%iTupdate==0 && iT<= iTmaxStart ) {
    // Smooth start curve, sinus
    // SinusStartScale<T,int> StartScale(iTmaxStart, T(1));

    // Smooth start curve, polynomial
    PolynomialStartScale<T,T> StartScale( iTmaxStart, T( 1 ) );

    // Creates and sets the Poiseuille inflow profile using functors
    T iTvec[1] = {T( iT )};
    T frac[1] = {};
    StartScale( frac,iTvec );
    AnalyticalConst2D<T,T> rho(frac[0]*converter.getLatticeDensityFromPhysPressure(pin));

    sLattice.defineRho( superGeometry, 3, rho );

    I think i did some mistake with the smoothened inlet pressure BC.

    Best regards,


    The code looks mostly fine, but try to smooth pin:

    AnalyticalConst2D<T,T> rho(converter.getLatticeDensityFromPhysPressure(frac[0]*pin));


    Unfortunately the System is still diverging giving something like

    [LatticeStatistics] step=625; t=0.25; uMax=inf; avEnergy=-nan; avRho=nan

    But the Parameters of the System are still the same as in the original example, except the introduction of the inlet and outlet pressure…


    Which Dynamics do you use? If you use the GuoZhaoDynamics, that could be the problem.
    Or maybe your pin is too high, try something like 0.01.


    Dear all,
    In this case which dynamic we should use? Is it ok if I use BGKdynamics?


    In release 1.3 the bug in GuoZhaoBGKdynamics has been fixed, therefore you can use either PorousBGKdynamics or GuoZhaoBGKdynamics.
    See for example “olb-1.3r0/examples/porousMedia/porousPoiseuille3d”, where you can switch between both dynamics.


Viewing 10 posts - 1 through 10 (of 10 total)
  • You must be logged in to reply to this topic.