[chymalaia]

Hello everyone,

Since there are no predefined solvers with those charcteristics, I'm trying to construct one wich uses VOF (it's a 2 liquid phases problem), solves for temperature and also considers buoyancy through the Boussinesq approximation (wich causes less processing effort than considering a compressible flow). For that I decided to modify the interFoam solver, having already added (to a copy of it) the energy equation by following some tutorials and threads in this forum.

My problem now is with the Boussinesq aproximation. For those who don't know it, the approximation basically uses one constant property of the liquid (Beta), gravity (g) the temperature (T) and one reference temperature (T0), as follows:

[...other terms of momentum equation] + (1 - Beta*(T - T0))*g

So if I'm correct I'll have to add that to the momentum equation, at UEqn.H. The problem is I don't know how to do that, mainly how to to define Beta to use it in that file and the notation to write that in the momentum equation.

Does anyone know where I could find information that would help me comprehend that? Until now, I have not found anything that suited...

Thanks in advance for any help.

Rosivaldo.

[Fabio1893]

Hello Rosivaldo,

if you need beta as a constant, you can simply define it in your UEqn.H file. It might be not the best way to do that, since you will have to change the value in your source code and recompile, if you want to use a new fluid, but it works. It is also possible to interpolate beta from a table, dependent on the temperature.

Fabian

[Diro7]

Hi,

assuming you are talking about the Foundation OpenFOAM versions, Boussinesq-based solvers where removed starting from version 7.
However you can still use version 6, or at least browse the source code.
Look for the buoyantBoussinesqPimpleFoam solver.

From UEqn.H

    if (pimple.momentumPredictor())
    {
        solve
        (
            UEqn
         ==
            fvc::reconstruct
            (
                (
                  - ghf*fvc::snGrad(rhok)
                  - fvc::snGrad(p_rgh)
                )*mesh.magSf()
            )
        );

        fvOptions.correct(U);
    }

The red line is the one you are looking for. The rhok field stands for the relative linearised density variation. It is updated in TEqn.H:

rhok = 1.0 - beta*(T - TRef);

Beta and TRef are simply read from the transportProperties file.
I'm not familiar with VOF solvers, but this should be a good starting point.
To understand the implementation, I suggest this nice blog entry:
https://caefn.com/openfoam/solvers-b...nesqpimplefoam

Hope it helps!

Andrea