[diego_angeli]

Dear FOAMers,

I don't know if this issue has ever been discussed before, but nevertheless I feel that it's important to raise it up, hoping that core OF developers can have a look at the matter.

According to the release notes of Version 7 (https://openfoam.org/release/7/) which I only read today, "The Boussinesq equation of state can now be applied to any buoyancy solver, deprecating the specialized buoyantBoussinesq[SP]impleFoam solvers".

As a convinced (and stubborn) user of the buoyantBoussinesq[SP]impleFoam solvers (and of the conjugateHeatFoam in the foam-extend version), and as a researcher and DNS code developer dealing with buoyancy-induced flows since 15 years, I am quite disappointed by the news, since, in my opinion, this corresponds to the end of the possibility of simulating buoyant flows with a correct implementation of the Boussinesq assumption in OpenFOAM.

To put it shortly:
- the Boussinesq assumption is not an equation of state
- in the Boussinesq assumption, density varies linearly with temperature only in the buoyancy term
- the usefulness of the Boussinesq assumption from a numerical point of view is the possiblility to exploit the advantages of the incompressible framework to simulate what, as a matter of fact, is a low-Mach number compressibility effect, using a gauge pressure to avoid numerical instability

In the OpenFOAM implementation (please do correct me if I am wrong):
- rho = rho0*(1 - beta*(T - T0)) is treated as an equation of state
- this is applied to all terms in the conservation equations
- now all heat transfer solvers are based on a compressible framework, and pressure is absolute. This could lead to ill-conditioned linear systems when, e.g. free convection problems, especially with open boundaries or with big, closed domains mimicking open spaces (think of, say, flow originating from an unconfined heated wire), are studied.

I am aware of the fact that there was a need to rationalize the heat transfer compartment in OpenFOAM, but I suggest finding a way to make the Boussinesq assumption a variant of the rhoConst equationOfState, only affecting the buoyancy term in the momentum equation, rather than an equationOfState as it is now. This would not only be more correct and coherent with the theory, but it also would restore the exact implementation of the buoyantBoussinesq[SP]impleFoam solvers.

One more remark: the "BernardCells" of the tutorial are misspelt. The phenomenon is known worldwide as Rayleigh-Bénard convection, the second surname being that of the French physicist Henri Bénard. I believe that a tutorial of such a fantastic code like OpenFOAM should have correct spelling of scientists' names.

Thanks a lot by now for anyone's attention, and replies/remarks

Diego

[kerim]

Dear All,
I was slightly surprised by the name of folder BernardCells when testing buoyantPimpleFoam. After visualizing the calculation results in Paraview, I realized that this is Bénard – Rayleigh convection in Bénard cells!

Kerim

[Diro7]

Hello Diego!

I'm starting to work on solvers based on the legacy buoyantBoussinesqPimpleFoam and I thought that updating them to the newest version of OF at the beginning of my work could be a good idea.

After several days I realized too that the new OF solvers were based on different physics. Also, the new buoyant[SP]impleFoam tutorials present some strange convergence behaviour.

I'm not new to OF, but I'm not a great expert of compressible solvers.
I'm definitely supposed to work on large, possibly closed domains and heavy liquids.
As an expert, do you have some specific advice or would you just recommend to switch back to older releases?

Thank you very much.

Andrea

[HPE]

To my observation, buoyantBoussinesq*Foams are still present in v1912. Please correct me if I am wrong.

[diego_angeli]

Originally Posted by Diro7

Hello Diego!

I'm starting to work on solvers based on the legacy buoyantBoussinesqPimpleFoam and I thought that updating them to the newest version of OF at the beginning of my work could be a good idea.

After several days I realized too that the new OF solvers were based on different physics. Also, the new buoyant[SP]impleFoam tutorials present some strange convergence behaviour.

I'm not new to OF, but I'm not a great expert of compressible solvers.
I'm definitely supposed to work on large, possibly closed domains and heavy liquids.
As an expert, do you have some specific advice or would you just recommend to switch back to older releases?

Thank you very much.

Andrea

Dear Andrea,

which cases are you targeting?

In the Boussinesq domain I would stick with the legacy solvers if you have to do fluid-only simulations, while I'd use the foam-extend conjugateHeat(Simple)Foam solvers for conjugate heat transfer.

Best regards

Diego

[diego_angeli]

Originally Posted by HPE

To my observation, buoyantBoussinesq*Foams are still present in v1912. Please correct me if I am wrong.

HPE, thanks for the info!

This could surely be useful to "patch back" the Boussinesq solvers in the .org version. However, my point is more conceptual about what is and what is not the Boussinesq approximation... and this seems to be ignored by the developers of the biggest open source continuum mechanics library.

[HPE]

Hi,

Might the org developers dont know the issue or this thread? I would suggest to submit an issue ticket in .org-Mantis, otherwise.

And/or in .com-GitLab even though the solvers exist there.

Thank you for your contributions.

[damu1414]

Dear Andrea,

I have been trying to get a converged solution for natural convection in square cavity case using buoyantSimpleFoam. This should be one of the simplest problems however, as you said, I am unsure about the default convergence criteria set in the solver. I would be grateful to have a discussion here. Apart from the convergence criteria, I would like to know if using buoyantPimpleFoam is a better choice for steady state problems. True that piso and pimple are transient solvers but, how bad is using that for steady state? Are there any other issues you have noticed in v7 of OpenFOAM specifically in buoyancy related solvers?

FYI, I am new to OF and have been using it for a year

Dear Diego and Kerim, request you too to join.

Thank you

[Diro7]

Dear damu1414,

I'm not sure I'm the best one around here to give general advice but I'll try to give you some hints that could be helpful with the discussion.
Also, I'm not so sure this thread is the most appropriate one for your questions, anyway...

Originally Posted by damu1414

I have been trying to get a converged solution for natural convection in square cavity case using buoyantSimpleFoam. This should be one of the simplest problems however, as you said, I am unsure about the default convergence criteria set in the solver. I would be grateful to have a discussion here.

I do not have great experience with density-based solvers, as after my initial troubles I switched back to OF6 to be able to use the Boussinesq-type ones (since my problem at hand hardly needed compressible fluid modelling).
Generally speaking, since you are talking about convergence it should be useful to include your case setup, and somewhat more precise info on results (maybe residuals etc).
If everything is set "correctly", or without further information, the only advice I can give you is to try to play with relaxation factors and see what appens.

Originally Posted by damu1414

Apart from the convergence criteria, I would like to know if using buoyantPimpleFoam is a better choice for steady state problems. True that piso and pimple are transient solvers but, how bad is using that for steady state?

I can't see any particular reason why you shouldn't at all use transient solvers for steady state calculations, apart from computational efficiency and some slightly more complex case setup.
As a matter of fact, in my (very limited) experience some problems can be quite tough to be brought to convergence with steady state algorithms.
Transient solvers may be helpful in such cases as stability can be enforced on more "physical" grounds (by means of Courant control etc) instead of a bare iterative procedure and somewhat arbitrary relaxation factors.

Originally Posted by damu1414

Are there any other issues you have noticed in v7 of OpenFOAM specifically in buoyancy related solvers?

As I said before, my adventures with OF7 and buoyant solvers came to an end quite quickly, so I'm afraid I'm of no help with regard to this

Best,

Andrea

[damu1414]

Dear Andrea

Many thanks for the response. Please see attached the case files. This is for a square cavity with left and right walls kept at hot and cold temperatures respectively and top and bottom being adiabatic. I have used the numericals to match Ra=10^4.
Physical properties of air taken at 1atm and 100degC. I have also included the mesh file(done in gambit).

FYI, I am trying to validate my results(Nusselt number to be specific) with the literature NATURAL CONVECTION OF AIR IN A SQUARE CAVITY : A BENCH MARK NUMERICAL SOLUTION by G. DE VAHL DAVIS

I would request you to have a look and please let me know if any further info is required.

Thank you

[Diro7]

Two quick preliminary considerations:

1) Why didn't you use blockMesh for mesh generation? Your geometry is so simple and already available in lots of OpenFOAM tutorials

2) Did you notice that under the buoyantSimpleFoam tutorials folder there is a buoyantCavity tutorial which is very similar to your case? It seems to run fine

Andrea

[damu1414]

Dear Andrea,

I was looking for an update before responding to you. And yes, there is.

I switched to OF5(I had to downgrade to Ubuntu 18.04 though) and applied the buoyantBoussinesqPimpleFoam solver for my case. To my surprise, it gave me results matching with the existing literature. As discussed earlier in this thread, it seems to be an issue with OF7. Am not sure.

To your question of why am not using blockMesh. I have been using gambit for quite sometime and some old habits die hard. . Also, am new to OF and learning in-house meshing tools.

Thanks a ton for you and Diego once again.

[kwardle]

I am interested in the original intent of this thread. Any update on reverting this change for Boussinesq solvers? I have certainly observed challenges with instability using the new solver that I never saw previously.

[Tobermory]

Sounds like the answer is "no" - no change. Possible workarounds are:
1. use rhoConst EOS in the solver (with rho = rho0) , and then add your own momentum source term for the Boussinesq gravity term (rho0*beta*(T - T0))
2. or, build you own Boussinesq solver from [sp]impleFoam, adding in the temperature equation and the Boussinesq gravity term.

Option 2 would be cleaner, since you avoid having to deal with p_rgh, but is obviously a little bit more effort. You can reuse much of the coding from buoyant[SP]impleFoam however, for the T equation. Good luck.

[Diro7]

Originally Posted by kwardle

I am interested in the original intent of this thread. Any update on reverting this change for Boussinesq solvers? I have certainly observed challenges with instability using the new solver that I never saw previously.

You can rather easily port the boussinesq solvers from OF6 to newer code releases. Some parts of the OF libraries needed by the legacy boussinesq solvers have been moved/renamed, but AFAIK all the functionality is still present somewhere and the exact behaviour can be reproduced.

Andrea

[kerim]

Dear damu1414,
I am sorry for late joining to this discussing.
First of all there is inconsistency between blockMeshDict (3D) and the files (2D) from 0 folder in your caseFiles.

I would like to repeat your OF5 results. Could you share your right caseFiles?

In the meantime, if you're interested, please see papre, where we've made a few corrections to buoyantCavity tutorial mentioned by Duro:
https://aip.scitation.org/doi/abs/10.1063/5.0071571

[fnohpat]

Dear FOAMers,

If someone could please advise me on the following: I am implementing buoyancy force with the Boussinesq approximation in the momentum equations as a source term in the "fvModels" file.

buoyancyBoussinesqForce
{
type coded;
select all;

field U;

codeInclude
#{
#};

codeAddSup
#{
Pout<< "**codeAddSup**" << endl;
#};

codeAddRhoSup
#{
Pout<< "**codeAddRhoSup**" << endl;
const scalar T0 = 293.5;
const scalar rho0 = 1.20329;
const scalar beta = 0.003407155;

const volScalarField& T = mesh().lookupObject<volScalarField>("T"); //T field

//double buoyancy = -0.001;
double buoyancy= -9.81*rho0*(1.0-beta*(T-T0));
//Pout<< "T = "<<T<< endl;

vectorField& USource = eqn.source();
USource += vector(0, buoyancy, 0); //,
#};

codeAddAlphaRhoSup
#{
Pout<< "**codeAddAlphaRhoSup**" << endl;
#};
}


But, I have the following error:

error: cannot convert ‘Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> >’ to ‘double’ in initialization
make: *** [/opt/openfoam11/wmake/rules/General/transform:26: Make/linux64GccDPInt32Opt/codedFvModelTemplate.o] Error 1

The error is in: double buoyancy= -9.81*rho0*(1.0-beta*(T-T0));

What is the correct variable type for buoyancy?

I'm solving natural convection in a square cavity with OpenFOAM 11. For verification purposes, I want to solve it using the Boussinesq approximation.

Thansk a lot
Felipe Noh

[Diro7]

Originally Posted by fnohpat

Dear FOAMers,

If someone could please advise me on the following: I am implementing buoyancy force with the Boussinesq approximation in the momentum equations as a source term in the "fvModels" file.

buoyancyBoussinesqForce
{
type coded;
select all;

field U;

codeInclude
#{
#};

codeAddSup
#{
Pout<< "**codeAddSup**" << endl;
#};

codeAddRhoSup
#{
Pout<< "**codeAddRhoSup**" << endl;
const scalar T0 = 293.5;
const scalar rho0 = 1.20329;
const scalar beta = 0.003407155;

const volScalarField& T = mesh().lookupObject<volScalarField>("T"); //T field

//double buoyancy = -0.001;
double buoyancy= -9.81*rho0*(1.0-beta*(T-T0));
//Pout<< "T = "<<T<< endl;

vectorField& USource = eqn.source();
USource += vector(0, buoyancy, 0); //,
#};

codeAddAlphaRhoSup
#{
Pout<< "**codeAddAlphaRhoSup**" << endl;
#};
}


But, I have the following error:

error: cannot convert ‘Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> >’ to ‘double’ in initialization
make: *** [/opt/openfoam11/wmake/rules/General/transform:26: Make/linux64GccDPInt32Opt/codedFvModelTemplate.o] Error 1

The error is in: double buoyancy= -9.81*rho0*(1.0-beta*(T-T0));

What is the correct variable type for buoyancy?

I'm solving natural convection in a square cavity with OpenFOAM 11. For verification purposes, I want to solve it using the Boussinesq approximation.

Thansk a lot
Felipe Noh

After a very quick glance:

- Since T is a volScalarField, your best chance is to cast buoyancy as a volScalarField as well. This is exactly what the error message is suggesting since, from what I remember, volScalarField is shorthand for <Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>>. For sure, casting it as a mere double cannot work.

- Also T - T0 might be problematic since T has physical dimension and T0 is a scalar (i.e. a mere double again) with no dimensions. But I'm not sure about this, I'm not writing openfoam since quite a long time now

[fnohpat]

Originally Posted by Diro7

After a very quick glance:

- Since T is a volScalarField, your best chance is to cast buoyancy as a volScalarField as well. This is exactly what the error message is suggesting since, from what I remember, volScalarField is shorthand for <Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>>. For sure, casting it as a mere double cannot work.

- Also T - T0 might be problematic since T has physical dimension and T0 is a scalar (i.e. a mere double again) with no dimensions. But I'm not sure about this, I'm not writing openfoam since quite a long time now

Thanks Andrea:


I have made some changes but I have more errors.

buoyancyBoussinesqForce
{
type coded;
select all;

field U;

codeInclude
#{
#include "createTime.H"
//#include "createMesh.H"
#};

codeAddSup
#{
Pout<< "**codeAddSup**" << endl;
#};

codeAddRhoSup
#{
Pout<< "**codeAddRhoSup**" << endl;
const scalar T0 = 293.5;
const scalar rho0 = 1.20329;
const scalar beta = 0.003407155;

const volScalarField& T = mesh().lookupObject<volScalarField>("T"); //T field

volScalarField buoyancy
(
IOobject
(
"buoyancy",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar ("buoyancy",dimensionSet(1,-2,-2,0,0,0,0),0)
);

vectorField& USource = eqn.source();

//double buoyancy = -0.001;
buoyancy= -9.81*rho0*(1.0-beta*(T-T0));
//Pout<< "T = "<<T<< endl;


USource += vector(0, buoyancy, 0); //,
#};

codeAddAlphaRhoSup
#{
Pout<< "**codeAddAlphaRhoSup**" << endl;
#};
}



Could I send you the complete case for you to review?

Thanks

[Diro7]

No, of course you can't. I don't have time to review your entire case and, as I said, I'm not working with openfoam anymore.
If you need someone to do your job for you, I'm sure you'll find plenty of people here who do that for a living.

Otherwise, you'll find plenty of free information in the forum or elsewhere to help you tackle your problem. Posting each iteration of your case every few hours or asking direct work from other people is not the proper way to ask for help here.

Good luck,

Andrea