[NightWing]

Hi All

I was searching some old OpenFOAM version and came across the solver named 'BoussinesqBuoyantFoam' which is available in OpenFoam 1.6 version.
This solver looks like the base code of a simple natural convection which contains the basic equations of U and T as shown below:

fvVectorMatrix UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
- fvm::laplacian(nu, U)
==
-beta*(T - T0)*g
);

solve(UEqn == -fvc::grad(p)); <-- This term is used when momentum predictor is used. (is it right?)

solve
(
fvm::ddt(T)
+ fvm::div(phi, T)
- fvm::laplacian(DT, T)
);

rho = rho0*(scalar(1) - beta*(T - T0));


But this solver is taken out in the after versions of OpenFOAM i.e. from 1.7 onwards and this solver got replaced (i feel so) with 'buoyantBoussinesqPimpleFoam' and buoyantBoussinesqSimpleFoam' and so on.

But I want to ask this.

-beta*(T - T0)*g <--- This term in the UEqn.H in the 'buoyantBoussinesqFoam' solver mentioned in OF 1.6 (i.e. the buoyancy term) is not seen in the R.H.S of UEqn.H in 'buoyantBoussinesqPimpleFoam or SimpleFoam.

fvVectorMatrix UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
+ turbulence->divDevReff(U)
);

UEqn.relax();

if (momentumPredictor)
{
solve
(
UEqn
==
fvc::reconstruct
(
(
fvc::interpolate(rhok)*(g & mesh.Sf())
- fvc::snGrad(p)*mesh.magSf()
)
)
);
}

can somone explain where the buoyancy term in this equation of UEqn.H is accounted?approxiamtion.

Or is it like the solver 'buoyantBoussineqFoam' mentioned in OF 1.6 is the base code and generalized code for a simple problem dealing with Boussineq approximation and it has nothing to with buoyantBoussinesqPimpleFoam?

Looking forward for suggestions
Best Regards!

[Tobi]

Hi,

a) please use code-tags
b) The newer version have rhok included and is defined as:

    // Kinematic density for buoyancy force
    volScalarField rhok
    (
        IOobject
        (
            "rhok",
            runTime.timeName(),
            mesh
        ),
        1.0 - beta*(T - TRef)
    );

I think it should be clear now, isn't it?

c) The older version used the PISO algorithm (I think), hence you are limited to small time-steps and there was no solver for steady-state solutions. Due to that they introduced PIMPLE and SIMPLE. Note that PIMPLE can work as PISO.

[NightWing]

 
// Kinematic density for buoyancy force
    volScalarField rhok
    (
        IOobject
        (
            "rhok",
            runTime.timeName(),
            mesh
        ),
        1.0 - beta*(T - TRef)
    );

yes, indeed. But as i went through the files of UEqn.H and PEqn.H, the reference of "rhok" can be found like this:

In UEqn.H

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

        fvOptions.correct(U);
    }

This bit of code suggest that, if I put my momentumpredictor 'yes' in 'fvSolutions' directory, this loops gets read. Then we have the equation which i suppose to be in the form

fvVectorMatrix UEqn
    (
        fvm::ddt(U)
      + fvm::div(phi, U)
      - fvm::laplacian(nuw, U)
     ==
       -g*beta*(T-Tref)
    );

    UEqn.relax();

is that what you have mentioned in the previous post???

and also in PEqn.H

surfaceScalarField phig(-rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());

which i Guess 'phig' talks about the turbulence part.

As i am doing a problem which is laminar, i have removed the turbulent files from the base solver.

In that case, the UEqn.H of the base solver 'buoyantBoussinesqPimpleFoam'

 fvVectorMatrix UEqn
    (
        fvm::ddt(U)
      + fvm::div(phi, U)
      + turbulence->divDevReff(U)
     ==
        fvOptions(U)
    );

    UEqn.relax();

the divDevReff(U) term, i have replaced with

fvVectorMatrix UEqn
    (
        fvm::ddt(U)
      + fvm::div(phi, U)
      - fvm::laplacian(nuw, U)
     ==
        fvOptions(U)
    );

    UEqn.relax();

    // along with the pimple loop as the base solver

Does this make sense at all? Am getting little by little of information everyday and keeps adding up my understanding of OF.

Thanks for your valuable time and suggestions

Best Regards

[Tobi]

Hi,

first of all. In Foam we are talking always about FLUXES. Hence it is not necessary to calculate the first guess of the velocity field (momentumPredictor on). Unfortunately I had never time to go deep into the pressure coupling (pEqn.H) but the phig are just some fluxes due to the fact that phi always represents fluxes. By the way, why you think that phig takes turbulence into account. Finally it will - yes, because we construct the UEqn.H (does not matter if we have momentumPredictor on or off) and using the UEqn matrix, we build the rAU field (that contains the turbulence etc).


To your question to the UEqn... first we construct the matrix:

    fvVectorMatrix UEqn
    (
        fvm::ddt(U)
      + fvm::div(phi, U)
      + turbulence->divDevReff(U)
     ==
        fvOptions(U)
    );

and if we have the MomentumPredictor turned on, we also take into account the buoyancy force:

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

That is finally that one:

        solve
        (
             fvm::ddt(U)
          + fvm::div(phi, U)
          + turbulence->divDevReff(U)
          ==
             fvOptions(U)
          + fvc::reconstruct
             (
                 (
                   - ghf*fvc::snGrad(rhok)
                   - fvc::snGrad(p_rgh)
                 )*mesh.magSf()
             )
        );

Finally to your turbulence modification. Yes it is correct. You also can check out pisoFoam.C

        fvVectorMatrix UEqn
        (
            fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)
          + g*rho*beta*(T-Tinf)
        );

[NightWing]

Originally Posted by Tobi

Finally to your turbulence modification. Yes it is correct. You also can check out pisoFoam.C

        fvVectorMatrix UEqn
        (
            fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)
          + g*rho*beta*(T-Tinf)
        );

would it be like this then for my laminar case (without rho in my buoyancy term as it is an in-compressible case)

fvVectorMatrix UEqn
        (
            fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)
          + g*beta*(T-Tinf)
        );

Moreover, if i use the term g*beta*(T-Tinf) as such to construct the UEqn.H, then i need not require to use the momentum predictor loop again (suppose it is turned off or on)

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

        fvOptions.correct(U);
    }

Hence the only terms in my UEqn.H would look like this

fvVectorMatrix UEqn
        (
            fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)
          + g*beta*(T-Tinf)
        );

       // Done without any momentum.predictor loop

Am i right Tobi?

Thank you for your valuable suggestions

Best Regards

[Tobi]

Hi,

okay,.. you want to solve boussinesq approximation in a laminar flow field. Then use boussinesqBouyantPimple/SimpleFaom. There is no need to modify anything. Keep your turbulence->divDevReff(U) and set turbulence or RANSModel to laminar. Thats it. From the version you mentioned (the old stuff), note that there (if I remember correctly) we use the pressure field p instead of the p_rgh pressure field. Therefore you get a more complex equation or have to rearrange them. Finally OpenFOAM changes the code from readability to "why are theses terms there and what does they mean". This is simple because of numerical reasons... if I remember the old code, it was much easiert to follow the equations but in any case, the new code is more stable, and that is the main point.

If you really want to understand why the solvers in 2.3.x or 3.x.x looks different to the old versions, derive the equations yourself, think about numerics and then you get the point. Good luck.

[NightWing]

Originally Posted by Tobi

Hi,

okay,.. you want to solve boussinesq approximation in a laminar flow field. Then use boussinesqBouyantPimple/SimpleFaom. There is no need to modify anything. Keep your turbulence->divDevReff(U) and set turbulence or RANSModel to laminar. Thats it. From the version you mentioned (the old stuff), note that there (if I remember correctly) we use the pressure field p instead of the p_rgh pressure field. Therefore you get a more complex equation or have to rearrange them.

Indeed, that' true! old solvers works around with p whereas recent solvers uses p_rgh

Why I tried to modify the solver 'buoyantBoussinesqPimpleFoam'?

if you look at this solver, it is a transient solver for turbulent flows. moreover this solver uses a radiation model.

As i wanted to work around with Boussinesq approximation and my flow laminar, i removed the dependency files and header files which calls the turbulence part and radiation models.

Because, if you look at the 0/folder, it helps to remove the fields of alphat, k, epsilon, nut etc, which is irrelavant for my problem. Moreover, my ultimate aim is to non-dimensionalize these equations so that i can study a simple natural convection problem which makes use of Boussinesq approximation.

Hence, by avoiding the complexities of turbulence and all, it would help to quickly use some non-dimensional parameters to convert these simple equations (without any turbulence and all).

So it seems like if I want to follow the basic equations as mentioned in old solvers like OF 1.6 (which talks about p rather than p_rgh), i want to change my PEqn.H of OF 2.x.x in the same way as old solvers which talks about p.

Or rather use OF 1.6 itself to solve my problem.

Right now, As like you suggested me in another thread, i am trying to incorporate the non-dimensional parameters in the post-processing stage.


Thank you for your time and valuable suggestions

Best regards

[Tobi]

If you are interested in steady-state, just use buossinesqBouyantSimpleFoam. I can not figure out the problem. In addition, if the solver uses radiation too, just set radiation to false (thats it). As I told you, the solver is for laminar and turbulent flow fields!