[Manuel CFD]

Dear Foamers,

I am new in solving natural ventilation problem by means of OpenFOAM and I encounter some difficulties.

My problem is represented by a building where there are two heat sources at a fixed temperature T=1000 K. I would like to study the convection flows that are generated inside the building through the windows and the doors. This is why I include in the domain also a volume outside the building represented by an hemisphere.

The problem I encounter is that the computation crashes after some iterations, it seems because density lower and upper limits become not physical (the upper limit becomes greater than 6 kg/m3 and the lower limit tends to zero until it becomes negative).

I use the original solver "buoyantSimpleFoam", so I treat the problem as compressible, and I use "realizableKE" as turbulence model.

The boundary conditions I impose are the following:

1. sky (patch representing the upper side of the hemisphere outside the building): totalPressure equal to 101325 Pa and an initial value equal to 101325 Pa too for the "p_rgh" field; "calculated" for the "p" field; "pressureInletOutletVelocity" for "U", imposing an inletValue equal to 0 m/s;
"inletOutlet" for "T", imposing an inletValue equal to 311 K; "inletOutlet" for both "k" and "epsilon", with the inletValue equal to 1e-5 for both of them;

2. heat sources: T equal to 1000 K; for "p_rgh" and "p" fields "buoyantPressure"; for "U" I impose a uniform value equal to 0 m/s; for "k" and "epsilon" I impose the corresponding wall functions;

3. ground and internal walls (these are modelled as baffles): T is set to zeroGradient; U is set equal to 0.0 m/s; p and p_rgh are set as buoyantPressure; and for k and epsilon the corresponding wall functions are set.

As I told you above, the problem is that the computation crashes after few iterations, and it seems that this is due to density divergence: maximum and minimum values are printed on screen at each iteration, and I notice that they continuously diverge till the crash (the minimum value becomes unphysically negative). I already read a post where they suggest to insert a modification inside the "createFields.H" file of buoyantSimpleFoam solver, which consists in inverting the computation of p_rgh by changing the expression "p_rgh = p- rho*gh" with "p = p_rgh + rho*gh"; but nothing seems to change, except the single iteration execution time, which increases a lot.

Does anybody have any suggestion for helping me? I have already tried several different set-ups, but nothing works.

I thank you a lot in advance for your help.

Manuel

[jherb]

Just one idea: Are you sure that the inlet velocity for the sky (inletValue) can be fixed to 0. You are doing a stationary simulation so if mass is flowing out it has to come in somewhere?

[adambarfi]

I think your velocity BC for the sky (U=0) made some problems. if it is very far from the bottom plate, according the pressureInletOutletVelocity you can assume that the velocity gradient is zero.

my suggestion for this boundary condition is p_rgh: buoyantPressure, p: calculated, U:zeroGradient and T:inletOutlet

they are conditions for a boundary that is far from a hot plate

[Manuel CFD]

Dear Mostafa and Joachim,

I thank you a lot for your quick answer and for your suggestions.

Actually I tried imposing the boundary conditions as you said, that is:

on "sky": for "U" zeroGradient; for "p_rgh" buoyantPressure; for "p" calculated; for "T" inletOulet.

For other patches I continued imposing the same bcs as in my previous post.

What I obtained is a result similar to my previous one, that is, the simulation seems to run well at the beginning, but early maximum and minimum rho values start diverging; in particular now the maximum value is at the level of 1.90 Kg/m3.

When I look at the solution by means of paraFoam, what I see is a fluid dynamic field which seems to be smooth everywhere, but in some local isolated points on the sky surface: at those points pressure and velocity values are extremely different from what I expected, that is, ~150000 Pa for pressure and more than 100 m/s for velocity.

What do you think about?

I thank you a lot for your attention and help.

Manuel

[jherb]

Have you checked the mesh quality (checkMesh)? This sounds like problems with the mesh on your sky boundary.

[adambarfi]

mhmm... what is the value of internalField for p and p_rgh? run it with zero value.

can you post the error description?

why are you using buoyantSimpleFoam? try BuoyantBoussinesqSimpleFoam.

[Manuel CFD]

Dear All,

I tried buoyantBoussinesqSimpleFoam but I didn't get any improvements.
I read on other pages of the Forum that a possible problem could be due to the boundary conditions buoyantPressure on the walls. They suggested to change in fixedFluxPressure. I did it but unfortunately still without positive results.
Concerning the mesh, checkMesh does not highlight any problem.
In order to simplify the problem and check the set-up of the simulation I prepared a simplified 2D case. On this geometry the settings give good results. Do you have any suggestions????

I thank you in advance for your help and attention

Manuel

[adambarfi]

I need some Information about your project, I did a simple model of three dimensional with the info you gave in your first post, but Rayliegh=1e5, It worked!

1- when you used buoyantBoussinesqSimpleFoam, what was the error?
2- are you using of turbulence models?
3- what schemes you used for discretization?

[Manuel CFD]

Dear Mostafa,

1) It crashes afeter having solved h equation; the error reported by the solver is "maximum iteration number exceeded" during the computation of T by the thermo object.

2) we are using realizableKE imposing on sky (atmosphere) inletOutlet condition for both k (inlet value equal to 1e-5) and epsilon (1e-6) and on corresponding wall functions on all the walls. We have two kind of walls: one with an assigned temperaure (T=1000 K) and the other one adiabatic.

3) We have used the following fvScheme file:

gradSchemes
{
default Gauss linear;
}

divSchemes
{
default bounded Gauss limitedLinear 1.0;
div(phi,U) bounded Gauss limitedLinearV 1.0;
div(phi,K) bounded Gauss limitedLinear 1.0;
div(phi,h) bounded Gauss limitedLinear 1.0;
div(phi,k) bounded Gauss limitedLinear 1.0;
div(phi,epsilon) bounded Gauss limitedLinear 1.0;
div(phi,omega) bounded Gauss limitedLinear 1.0;
div((muEff*dev2(T(grad(U))))) Gauss linear corrected;
}

laplacianSchemes
{
default Gauss linear corrected;
laplacian(muEff,U) Gauss linear corrected;
laplacian(Dp,p_rgh) Gauss linear corrected;
laplacian(alphaEff,h) Gauss linear corrected;
laplacian(DkEff,k) Gauss linear corrected;
laplacian(DepsilonEff,epsilon) Gauss linear corrected;
laplacian(DomegaEff,omega) Gauss linear corrected;
}

interpolationSchemes
{
default linear;
}

We have also tried to use for all divergence schemes upwind instead of the limitedLinear but nothing has changed.

Thank you,

Manuel

[Manuel CFD]

I miss to tell you that on sky patch I imposed buoyantPressure for the p_rgh field.

Manuel

[adambarfi]

according to my experiences "maximum iteration number exceeded" error occurs when the convergence is very weak. this may happen for various reasons: mesh quality, mesh size, your time step ...

try the following:
1- change the geometry scale by 0.1 and then change other parameters so that your dimensionless numbers stay the same as before.
2- increase mesh number
3- decrease the deltaT
4- try Gauss upwind for Divergence scheme
5- download some papers that are related to your project and see what is their BCs, mesh sizes and the discretization schemes.

I hope it helps you

[wc34071209]

You may need to try to turn off the gravity to see if the problem still exists or not.

[Manuel CFD]

Dear All,

After a lot of attempts, I managed to solve the problem inserting bounding values for all the fields (temperature, velocity, pressure and density). In this way I obtained a converged solution which shows realistic flow fields. This computation was performed by imposing a 1000 K fixed temperature on heat source surfaces.
Starting from this solution, I tried to perform another computation, replacing this temperature boundary condition with the compressible::turbulentHeatFluxTemperature one, in order to specify heat fluxes on these surfaces (in the range between 500-30000 W/m2). Unfortunately this boundary condition is causing a lot of problems, because the temperature on these surfaces is exceeding the 2000 K upper bounding value. I use this bounding value because the highest temperature that I am going to expect is about 1600-1700 K. I do not understand why it happens.
I tried also to impose lower heat fluxes values, limited to the 5000 W/m2 upper value. Also in this case the temperature diverge.
Furthermore using the wallHeatFlux post-processing utility, I am checking the wall heat fluxes on these surfaces and they do not match with the imposed ones.
Even adopting the simplified 2D geometry mentioned in the previous posts, the problem still remains.
Have you any suggestions??

Thank you,

Manuel

[olivierdambron]

Dear Manuel,

I am experiencing the very same issues that you describe and I have nto been able to solve the 'maximum number of interation exceeded'.

I am not very confident with my BCs and Schemes.

Any chance you could share the case that worked for you? I'm trying to analyze the airflow induced by buoyancy and radiation inside a solar chimney using buoyantSimpleFoam.

Kind regards,
Olivier

[Hip2BL7]

I was experiencing the same issue as you OP, but I managed to cure my issue by under-relaxing my rho values to <0.1 (my solution then converged quickly). Think I'm going to try buoyantPimpleFoam instead.


My case is turbulent flow through a heated pipe, Re = 5e5.

[MariMari]

Hi
Your post was good. thank you very much.
Good luck

excel square-root