[Tobi]

I don't think it is a solver-bug but there is probably something missing I don't get until now or there is no appropriate boundary condition (and I don't have time to investigate into that problem). So you have to check it out yourself.

[thiagopl]

Hi everybody,

I'm currently working on a problem of natural convection in a open domain using buoyantPimpleFoam and, as many of you, I faced the same problems on the outlet boundary, say, backflow, accelarating flow, divergence, etc. I tried different combinations of boundary conditions, most of them discussed here in this thread.
Recently I found this work which suggests the combination for p_rgh/U of fixedFluxPressure/zeroGradient for the top boundary and totalPressure/pressureInletOutletVelocity for the sides and bottom boundaries.
I tried those on a test case attached and I've got a good result for the flow field on outlet (see video). In my case the bottom part of the domain is a wall. The boundary conditions were the following:

Pressure: p_rgh

dimensions      [1 -1 -2 0 0 0 0];

internalField   uniform 101325;

boundaryField
{
    outlet
    {
        type            fixedFluxPressure;
        value           $internalField;        
    }
    sides
    {
        type            totalPressure;
        p0              $internalField;
        value           $internalField;
    }
    heat
    {
        type            fixedFluxPressure;
        value           $internalField;
    } 
    base
    {
        type            fixedFluxPressure;
        value           $internalField;
    }
    frontAndBack
    {
        type            empty;
    }
}

Velocity: U

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
    outlet
    {
        type            zeroGradient;      
    }
    sides
    {
        type            pressureInletOutletVelocity;
        inletValue	 uniform (0 0 0);     
        value           uniform (0 0 0);
    }
    heat
    {
        type            noSlip;
    }    
    base
    {
        type            noSlip;
    }
    frontAndBack
    {
        type            empty;
    }
}

Temperature: T

dimensions      [0 0 0 1 0 0 0];

internalField   uniform 293.15;

boundaryField
{
    outlet
    {
        type            inletOutlet;
        inletValue      uniform 293.15;
        value           uniform 293.15;
    }
    sides
    {
        type            inletOutlet;
        inletValue      uniform 293.15;
        value           uniform 293.15;
    }
    heat
    {
        type            fixedValue;
        value           uniform 303.15;        
    } 
    base
    {
        type            zeroGradient;
    }
    frontAndBack
    {
        type            empty;
    }
}

Hope it helps.

[marco.müller]

Just contributing my findings, since I had the same issue, in case someone will find it helpful:

- I was not able to solve my (steady state) problem with buoyantBoussinesqSimpleFoam. It seems to have a current (v2006 for me) implementation that makes it difficult to deal with openings, see:
Wrong implementation of the Boussinesq assumption in heat transfer solvers
- Using buoyantSimpleFoam worked fine for me (without momentumPredictor, which introduced the same artificial convection with middle speed velocities as buoyantBoussinesqSimpleFoam)
- BC setup:
--- lower opening ("inlet"): prghTotalPressure + pressureInletOutletVelocity
--- upper opening ("outlet"): fixedFluxPressure + pressureInletOutletVelocity
- lastly - since I also had forced convection inside my domain -, I successfully used the limitVelocity option to get a stable initialization

[Alczem]

Following the release of OpenFOAM v11, I also want to share some recent tests I did about natural convection and open boundaries since I have been looking into it quite a lot. I think I finally have a decent setup for a steadystate simulation using the fluid solver of OF v11 (similar to buoyantSimpleFoam).

My case is a simple hot box inside a fully open domain, so the air can come and leave frome any side.

I used OpenFOAM v11 because it includes a switch to initialize the hydrostatic pressure in fvSolution so you can use the ph_rgh field and the prghTotalHydrostaticPressure condition in p_rgh. In OpenFOAM v2306, there is a functionObject that does the same initialization but I have yet to try it

Here are the boundary conditions:

p

internalField   uniform 101325;

boundaryField
{
    #includeEtc "caseDicts/setConstraintTypes"

   ".*"
   {
       type            calculated;
       value           uniform 101325;
   }
   
   "(sky|bottom|walls)"
   {
       type            calculated;
       value           uniform 101325;
   }
   
   "vol.*"
    {
        type            empty;
    }
}

p_rgh

internalField   uniform 0;

boundaryField
{
    #includeEtc "caseDicts/setConstraintTypes"

    ".*"
    {
        type            fixedFluxPressure;
        
        value           $internalField;
    }

    bottom
    {
        type            prghTotalHydrostaticPressure;
        meanValue        $internalField;
        p0              $internalField;
        value           $internalField;
    }
    
    sky
    {
        type            prghTotalHydrostaticPressure;
        meanValue        $internalField;
        p0              $internalField;
        value           $internalField;
    }
    
    walls
    {
        type            prghTotalHydrostaticPressure;
        meanValue        $internalField;
        p0              $internalField;
        value           $internalField;
    }
    
    "vol.*"
    {
        type            empty;
    }
}

ph_rgh

internalField   uniform 0;

boundaryField
{
    
    ".*"
    {
        type            fixedFluxPressure;
        
        value           $internalField;
    }

    bottom
    {
        type            fixedFluxPressure;
        meanValue        $internalField;
        p0              $internalField;
        value           $internalField;
    }
    
    sky
    {
        type            fixedValue;
        meanValue        $internalField;
        p0              $internalField;
        value           $internalField;
    }
    
    walls
    {
        type            fixedFluxPressure;
        meanValue        $internalField;
        p0              $internalField;
        value           $internalField;
    }
    
    "vol.*"
    {
        type            empty;
    }

}

U

internalField   uniform (0 0 0.0);

boundaryField
{
    #includeEtc "caseDicts/setConstraintTypes"

    ".*"
    {
        type            fixedValue;
        value           $internalField;
    }
    
    walls
    {
        type            slip;
    }

    bottom
    {
        type            zeroGradient;
        //inletValue        uniform (0 0 0);
        value           $internalField;
    }
    
    sky
    {
        type            zeroGradient; //pressureInletOutletVelocity;
        inletValue        uniform (0 0 0);
        value           $internalField;
    }
    
    walls
    {
        type            zeroGradient; //pressureInletOutletVelocity;
        inletValue        uniform (0 0 0);
        value           $internalField;
    }
    
    "vol.*"
    {
        type            empty;
    }

}

T

internalField   uniform 293.15;

boundaryField
{
    ".*"
    {
        type            zeroGradient;
        value           uniform 293.15;
        //value           $internalField;
    }
    
    hot_walls
    {
        type            fixedValue;
        value           uniform 393.15;
    }

    "(sky|bottom|walls)"
    {
        type            inletOutlet;
        inletValue      uniform 293.15;
        value           $internalField;
    }
    
    "vol.*"
    {
        type            empty;
    }
        
}

pRef (located in constant, and that's why p_rgh and ph_rgh are set at 0 Pa)

dimensions      [1 -1 -2 0 0 0 0];
value           101325;

fvSolution

PIMPLE
{
    momentumPredictor no;
    nNonOrthogonalCorrectors 1;
    //pRefCell        0;
    //pRefValue       0;
    //consistent        yes;
    
    hydrostaticInitialisation yes;
    nHydrostaticCorrectors 10;

}

fvSchemes

ddtSchemes
{
    default     steadyState;
}

gradSchemes
{
    default     cellLimited Gauss linear 1;
}

divSchemes
{
    default         none;
    div(phi,U)      bounded Gauss linearUpwindV default;
    div(phi,K)      Gauss linearUpwind default;
    div(phi,h)      bounded Gauss linearUpwind default;
    div(phi,k)      bounded Gauss linearUpwind default;
    div(phi,epsilon) bounded Gauss linearUpwind default;
    div(phi,omega)     bounded Gauss linearUpwind default;
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear uncorrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         uncorrected;
}

wallDist
{
    method            meshWave;
}

physicalProperties

thermoType
{
    type            heRhoThermo;
    mixture         pureMixture;
    transport       const;
    thermo          hConst;
    equationOfState perfectGas;
    specie          specie;
    energy          sensibleEnthalpy;
}

dpdt            off;

mixture
{
    specie
    {
        molWeight       28.96;
    }
    thermodynamics
    {
        Cp              1004.4;
        Hf              0;
    }
    equationOfState
    {
        rho            1.2;
        pRef        101325;
    }
    transport
    {
        mu              1.831e-05;
        Pr              0.705;
    }
}

A few notes:

• I think the most stable setup for ph_rgh is to have one of the patches at fixedValue and the rest at fixedFluxPressure although it seems to work if fixedFluxPressure is used everywhere but the velocity is not as stable. Using the highest horizontal patch for the fixedValue condition seems ok. I ran the case without any hot patches to see if the air would settle, and the maximume velocity was around 2 cm/s (the domain is 30m tall).
• Setting the U boundary conditions to zeroGradient or pressureInletOutletVelocity did seem to change anything.
• Turbulence needs to be on even if nothing is happening inside the domain. The temperature exploded in my case with turbulence off.
• The main region of interest is too close to the borders of my domain to prevent any perturbations, so I guess better results could be achieved with a larger domain. The vortices are probably due to the narrowness of the domain.
• Watch out when stopping and running the same case again without any modification. The hydrostatic initialization actually modifies the pressure files which means that stopping and restarting the case will initialize the files that are already initialized. It can lead to an unstable simulation (at least for me).

Here is a picture of the temperature field with the surface LIC representation of the velocity. It looks decent to me but I would like to have other opinions on this Cheers!

[thiagopl]

Hi Alczem,

Nice work. Although I posted this and also proposed this thread, which lead to conclusions, I still have problems in setting up natural convection problems in open domains.

I found that OF v10 also includes this switch to initialize the hydrostatic pressure in fvSolution.

I'm testing your setup, but in a transient case (good results by now!). I 'ill post the case here once I have conclusive results.

[Alczem]

Hey!


I remember reading your various posts, I think I read almost every thread related to this kind of simulation


Yep, the switch was already implemented in OF v10, but I mainly use the ESI version and only recently have I been using the foundation version which has some pros (and cons, still figuring it out!).




From some posts (by Tobias Holzmann I think), it sounded as if free convection might only be properly simulated through transient computation, or at least using a transient initial state to initialize the steady run, which is... not encouraging! I tried local time stepping and large CFL numbers but it still so slow, especially if you are simulating fluids and solids in a CHT case hopefully some common setup will emerge from all our tests

[Alczem]

Back again with a similar case to my previous post, only on a much larger domain (which is 2D now and triangles only, I meshed quickly with Salome and had to be fast).


The pressureInletOutletVelocity condition seems actually more stable than zeroGradient, but they achieve very similar results.


A lot of undesired currents remain present near the boundaries and I did not manage to completely mitigate their effects. They are also present when the hot box was disabled to try to reach a perfectly still atmosphere.


Overall I am pleased with the results, but it is slow: it takes more than 5000 iterations to have a satisfactory flow. Maybe the numerics need more tinkering and other boundary conditions can achieve similar results faster. I still take this as a win

[AliVali]

Hi Bodo1993,

It is about three years since your post. I hope you have had found a solution to this problem and are able to share it with me. I am facing exactly the same issue.

Thanks,
Ali

[AliVali]

Originally Posted by Bodo1993

Hi Tobi,

I encounter the same problem in twoLiquidMixingFoam, a solver for two miscible incompressible fluids. In this solver, the mixture density is interpolated linearly (not a constant value). I face the issue of backflow regularly in my simulations. I start to notice the effect once the mixed phase reaches the outlet boundary. I have tried almost all combinations of the outlet and inlet boundaries, but the problem persists. As you mentioned, when I turn the gravity off, I do not observe a back flow. It tells me that there is an issue with the buoyancy term here. I wish there is a solution for this problem and I am wondering if it is a bug in the solver.

I would greatly appreciate your thoughts and suggestions.

Hi Bodo1993,

It is about three years since your post. I hope you have had found a solution to this problem and are able to share it with me. I am facing exactly the same issue.

Thanks,
Ali

[Alczem]

Hey!


I recently have had to work on buoyancy effects again, and I came back to this thread. I managed to slightly improve on my previous solution, so I figured I would share and hopefully get some feedback and see if others have managed to get their open boundaries to work


Basically, I still use the hydrostaticPressure FO (pRef set to 0) from the ESI version to initialize ph_rgh, p_rgh and p. Boundary conditions are as follows:

• ph_rgh : fixedValue for the horizontal open patches, fixedFluxPressure everywhere else, 1bar as initial value
• p_rgh : prghTotalHydrostaticPressure for the openings (both horizontal and vertical), 1bar as initial value
• p : calculated, 1bar as initial value
• U : pressureInletOutletVelocity for the openings (both horizontal and vertical)

The air is set to perfectGas (for incompressiblePerfectGas, totalPressure BCs are good enough).


Now the part I did not understand last time was how to let hot air leave the domain freely instead of "bouncing" on it. To achieve that, hRef has to be set to the height of the horizontal patch that will act as an outlet. Meaning that if you expect your flow to leave by the bottom boundary (in case of cold air), hRef should be set to the coordinate of that boudary.


One more thing that I found useful: once the hydrostaticPressure FO has been run and the fields are initialized, we can disable it and use the readFields FO to read ph_rgh and be able to use the prghTotalHydrostaticPressure for p_rgh without needing to initialize ph_rgh again (because the FO will overwrite the existing pressure fields).



I found this really neat to initialize a transient simulation with a steady-state run and keep using the prghTotalHydrostatic BC but without the need to overwrite the pressure fields at the start of the transient case (which is what would have happened if I used the hydrostaticPressure FO to access ph_rgh).


The behavior of vertical patches is still a bit unclear to me, I intend to run more tests, maybe with a lateral wind pushing the hot plume in the side patches.






I believe that concludes my quest for open boundaries with buoyancy... for now Maybe someone can tell us more on the proper way to treat buoyancy and boundaries in the future, but that is enough for me!