buoyantSimpleFoam and watertank
 

Dear Foamers,

I am working (since 2 weeks) on a very simple simulation.

What I want to simulate:

Something like that: http://www.wiga-energietechnik.de/bi...mage/lwsp2.gif


What I did:

- I meshed the whole geometry with a corse and very fine mesh
- I build polynoms for water thermodynamics (30°C-70°C)
- I changed the thermodynamics for water
- Simulation is LAMINAR

- Inlet 4e-5 m³/s

- At the inlet I have a very simple pipe installation but the solver blow up every time so I just set an Inlet + Outlet (thats all - see pictures).

Now my problem:

Every BC I set make problems.

I am not 100% sure how I should set the p_rgh BC for inlet/outlet/wall.

p is calculated.


For U and T its clear.

The solver is working just for 1 or none iterations.
If I set of the gravity the simulation is working.

It seems that the solver is calculating my water with a compressibility because after the first time step I get extrem huge velocity fields in the big domain.


I tried a lot of BC for U + p_rgh - fixed at the outlet - pressureInletOutletVelocty etc.

Does someone can give me a hint how to set these BC right?

Relaxationfactors are decreased to 0.1.
linearUpwind + limitedLinear schemes are used etc...


Interesting fact:

without gravitation the simulation is working.
With gravitation the mass flux cant be calculated:
For gravity I have to set the outlet p_rgh to fixedValue that I can calculate the first timestep but after that the solver blow up again. A picture of p, p_rgh and U is included in the attachment.

The error message:
Compared with the velocity field its clear why its blow up but I can not fix it.


Any suggestions would be appreciated.

Regards Tobi

Hi Tobi,

seems to me as your problem might arise due to the initial conditions of p. In your final solution the p field should be close to the hydrostatic field and p_rgh close to constant. However, if you assume p to be constant in the beginning then you get the p_rgh field shown in your plot which leads to high velocities and possibly even to a crash. The problem is that in the very first step the solver uses the p-field to calculate the p_rgh field although your boundary conditions are set for p_rgh and p is set to calculated. So what you can do is either adapt the solver (which is just a change of one line) or probably a little bit easier use funkySetFields to set the initial pressure field to the hydrostatic pressure.

Hope that helps
Hannes

Hi Tobi,

As a by-pass you could run initially with gravity turned off. Than after say 200 iterations turn gravity on. Unfortunately it is not run-time modifiable so you have to stop (and save) the run, change gravity and rerun from latestTime.

Regards,
Tom

Hi,

thanks for your replays.


I tried Toms hints.
The pictures show p, p_rgh and U after 500 iterations.

Then i turned gravity on and the 501 step is shown in the 2nd replay.

@hannes: what lines had to be modified?
@all: could it be possible that my polynomes are wrong?

Here are the pictures at 501 iterations.
After that my solver blow up:
One hint for p_rgh:

Normally I should have such a pressure field shown in the picture, am I not?
Because my tank is 2m high and the fluid is water (1000kg/m³).

So as I see in this formula the higher I get in the tank the lower p0 should be because:

Is that correct?
Therefor p and p_rgh should not be the same?

Hi all,

just one question. Why is in interFoam the p0 calculation for nonclosedVolumes like that:
and in buoyant***Foam like that:
Regards
Tobi

Hi all,

I made a test with chtMultiRegionSimpleFoam with only one region.
Just a very simple case - box with inlet and outlet at the top of the box.

If I use the thermodynamics out of the liquidHeater tutorial:
Its working.
After switching to my own thermodynamic with the polynoms its blow up after the 4th iteration:

Therefor I made a test just with the first coefficient like:
for all other polynoms too. But the result is the same - after 3 iterations the solver blow up due to maximum Iteration reached:

Any hints?
I am out of mind and have no further ideas at the moment.

Regards
Tobi

Hi all,

just notice:

I am stupid! :p


All values of mu are wrong. Instead of writing:

503,89e-6

I wrote

0,50389e-6


After changing this, the easy test case is working. Now I am going to check if my bigger project is working too.

Okay.

I checked it out with the other case.
Without gravity its working (bouyantSimpleFoam and chtMultiRegionSimpleFoam).

But after switching on the gravity its still the same. :(

Regards Tobi

Perhaps your problem is the same as described in this thread:
http://www.cfd-online.com/Forums/ope...roperties.html

You could also try the transient solver (buoyantPimpleFoam)

Dear Tobi,

the pictures you showed match exactly the problem we had. I'm pretty sure that your results will improve when you start the simulations with hydrostatic pressure distribution for p. What happens is when you turn on gravity the pressure field rapidly changes in order to be consistent so that's why you get those strong gradients in the p_rgh field which again results in the high velocities.
If you don't want to set your p field with funkySetFields then you should adapt your solver. Copy the solver and the only thing you have to change in createFields.H is the line reading
For the buoyantSimpleFoam solver it is line 71
Replace it by solving for p:
This way you initialise using the constant p_rgh field instead of the constant p field which is much closer to the real solution.

Regards
Hannes

Hi all,

thanks for the hints but it is still not working.
On my easy case the chtMultiRegionSimpleFoam solver (only with one domain) is working. But the buoyantSimpleFoam is not working. Additionally I tryed the hints hannes said but not with success.


In the attachment you find a picture of the chtMultiSImpleFoam solver and the solution.
I have no further idea at the moment.

For me it seems that the buoyantSimpleFoam is not a good solver for the thing I want to do. Furthermore the chtMultiSimpleFoam solver is not working in my big case too.

Thanks for your help but I think I have to give up on that project.
Maybe "WATER" is not very common for the solvers.

At least I had a look into the cht solver and the createFields.H file.
There is the same calculation as in the buoyantSimpleFoam.
So I have no idea why this solver is working and the other one not :/




Regards Tobi

Hi all,

now the pressure fields are the same (buoyantSimpleFoam and chtMultiRegionSimpleFoam) .

I had to change the schemes to get the same results :)


Now I am going to check if its working wit a bigger domain.

Summary with a bigger domain:


- chtMultiRegionSimpleFoam isn 't working anymore
- buoyantSimpleFoam isn 't working anymore
- myBuoyantSimpleFoam with the modification in the createFields.H is working.

Hannes I think the work is done now.
I will check my official geometry now.

I keep you posted.

Hi Hannes - your hint is working fine :)
Thanks a lot.


Just one question to that.
Why isn 't it implemented as you wrote?

I think there is any reason for that?

Additionally with the PIMPLE algorithm it is not working.
Do you have any experience with that?

Maybe I will initialize it with simple and then switch to pimple.

Regards Tobi

Hi all,

with the modified buoyantSimpleFoam solver my case is working and the steady state result is very nice.

After initialize this solution with the buoyantPimpleFoam solver I get crazy p_rgh and p fields again :(

Hi Tobi,

I can't really tell you why it is implemented that way, I'm not aware of any restrictions at that point. Maybe the thinking is, that it is easier to initialise a pressure field which is a little more intuitiv then to initialise p_rgh when not starting with constant fields.
However, concerning the problem with buoyantPimpleFoam I could only guess. First of all, one thing that might become neccessary is to increase your writePrecision in controlDict, the standard six (or eight?) digits are by far not sufficient when small fluctuations in p_rgh are concerned, so that might be a reason why a restart might fail. Otherwise it should be possible to run the simulation with the buoyantPimpleFoam starting from a constant field when the same change to the solver is performed (change in createFields.H).
Could you provide some more information on the crash if the above does not help?

Hannes

Dear Hannes,

thanks for your replay.

I also tried to start the simulation with the change in the createField.H.
But without success.

Additionally I changed the time precision like you said but the same - crash.
Here is the output:

I added a picture of the initial solution.

Regards Tobi

If I read your log file correctly, you don't use the nOuterCorrectors loop. Have you tried nOuterCorrectors > 1 with underrelaxation?

Hi,

PIMPLE and underrelaxation?
Is underrelaxation not changing the real solution (accuracy in time)?

Well I used nOuterCorrectors > 1 but without success anyway. The error is different. I reach the maximum iteration in the temperature calculation:

Hi Tobi,

what I meant was to increase the write-precision and not the time-precision. At least that helps me a lot when restarting simulations.
Perhaps you could also try with larger time-steps. Your Courant-number is extremely low. Maybe it's also worth trying with nNonOrthogonalCorrectors > 0 (I reckon your mesh is not 100% orthogonal?).

Hannes

Hi Hannes,

I increased the "writePrecision" and not the "timePrecision" :)
At the moment I increased both to 15, set my dT lower and the nonOrthogonal to 1.

My mesh is good:
At least I tried to switch of the gravity again but without success.

Regards
Tobi

The idea of PIMPLE is to combine the PISO and the SIMPLE algorithms. You iterate several times over one time step and use underrelaxation to get a stable solution. Of course, if you only did one iteration per time step + underrelaxation your solution would not be correct. Therefore, you do multiple iterations (check your residuals).

So I would suggest, set nOuterCorrectors to e. g. 20 and and for the relaxation factors something like that:
You probably have to adjust these factors. See also the tutorials for pimpleFoam (.../OpenFOAM/OpenFOAM-2.2.x/tutorials/incompressible/pimpleFoam/...)

Hi Joachim,

I know that the PIMPLE solver is a combination of PISO + SIMPLE. But your hint makes sence.

To your suggestion to have a look at the tutorials.
In the incompressible solver there are given the relaxation factors but always no underrelaxation. The factors always set to 1. In all other tutorials (used PIMPLE) its the same.

But as I said befor your hint makes sence to make nOuterCorrectors = 20 or bigger to resolve the underrelaxation.


I will give it a try.

Hi all,

there is something strange in my simulation.
After changing the fvSolution with underrelaxation and using nOuterCorrections you can see that p_rgh is divergating:

If someone what to have a look at the case you can download it here:

www.holzmann-cfd.de/upload/case_schichtspeicher.tar.gz


Regards Tobi

Add the following line to your thermophysicalProperties:
see: http://www.openfoam.org/version2.2.0/thermophysical.php
Then your case is running (with the "vanilla" buoyantPimpleFoam solver). Here is my fvSolution file which I also changed (nNonOrthogonalCorrectors, the values of the relaxationFactors) but which might also run with your settings:

Hi Joachim,


thanks for having a look at the case.
Its working ? On my computer its not working :(

I changed the fvSolution to yours and added the additionally line in the thermodynamics to switch the pressure-work off.

I did not pass the first iteration

I did not change the standard OpenFOAM solver!

The only other change I made was that I run the case on 6 processors in parallel (as specified in your decompseParDict). Here are all changes:
I checked again and it also runs as a single process.
My OpenFOAM version is the commit 9f19321c01ca1e30d339b48eef3f70f5534ca7f4 from the 2.2.x branch at github.

And here the log (single process) until I stopped the run with CTRL+C:

Hi,

my fault was the entry in the thermodynamics...
Thanks for your help!
Regards Tobi

Now I am really confused:
off or no shouldn't matter.

But it should read dpdt not dtdp.

Hi,

that was a mistake I made. I corrected it in the last post.
That "no" "false" "off" should be the same is correct. I made any other strange mistake that I can not say at the moment. Maybe I wrote dtdp instead of dpdt in the thermodynamics ?!


Thanks for your help. Now everything is working.
Regards
Tobi

Output dictionary files into log
 

There is a useful debug switch in the global controlDict:
This global controlDict has to be copied from $WM_PROJECT_DIR/etc/ to $HOME/.OpenFOAM/2.2.x (or what ever version you use). Then make the changed marked in bold font above. All dictionary files willl be printed to STDOUT whenever they are (re)read. This is very useful if you want to know which settings you really used for a simulation.

Hi Tobi,

I have a similar application and want to study the degradation of the stratification in the storage tank. I have initialized a temperature profile and want to see how the heat loss and the convective mixing currents at the wall affect the stratification under no-flow conditions. My case works well with buoyantBoussinesqPimpleFoam. Now I want to run it with buoyantPimpleFoam because I want use chtMultiRegion later. The simulation stops after a couple of iterations due to pressure problems.

Could you please summarize your steps?
- Are you using now the solver with the changed createFields.H and have you initialized p or p_rgh?
- What boundary conditions are using for p and p_rgh?
- I have seen that you are using GAMG for the pressure, while this solver works quite good for my incompressible case, it does not work with buoyantPimpleFoam.
- Could you post your fvSolution and fvScheme?

I have also disabled dpdt and use polynomials in the thermophysical properties.

I tried to download your zip file but it is not available anymore.

Thanks and best regards

Bernhard

Dear Bernhard,

I think in your case it is possible to start with your steady-state solution, arent you?

Can you give me a hint why you have to use buoyantPimpleFoam befor using cht?

How are your Settings of the fvSolution / pimple algorithm and your Relaxation factors?

What about your schemes?

Is is a forced or buoyancy driven flow?

Regards Tobi

It is a buoyancy-driven case, the geometry is a cylinder as wedge with upper, lower and outer wall. There is no inlet or outlet. I use setFields and funkySetFields to initialize a temperature profile so that the upper region is 823 K and the lower 563 K and a transition zone in between (my fluid is molten salt). I use groovyBC for the outer wall to apply a variable gradient. The aim is to investigate the influence of the wall currents caused by the heat loss on the stratification. The result for the incompressible case is looking like that.
Attachment 27918

Attachment 27919

If I switch now to buoyantPimpleFoam and make the necessary changes for the BC's the simulation is not really running (only for some iterations). I am using fixedFluxPressure for p_rgh (for the walls) and calculated for the p. I tried to initialize a hydrostatic pressure for p with funkySetField but it does not help. While I was using a linearUpwind Scheme in the incompressible case for U and T, I have set all div schemes to Gauss upwind for the buoyantPimpleFoam case. I am using further a lowRe turbulence model. Increasing the number of correctors did not really help.

This is my fvSolution file:
and these are the thermophysical properties
So far, to keep it as simple as possible, only the density depends on the temperature.
The fvSchemes are quite simple
I tried a lot of settings but the pressure blows up everytime.
However, if I change the patch type for the upper boundary (the upper face area of the cylinder) from wall to patch and use an outletInlet BC for U instead of a no-slip condition and a fixedValue or totalPressure for the p_rgh the simulation is running. I dont know if these BCs are suitable to treat the upper boundary as an opening representing a free surface if the storage is not completely full.

While the GAMG solver for the pressure was not working at all for the case with only walls as boundaries, it is faster than the PCG solver with the top patch as opening.

I am using buoyantPimpleFoam to check the required settings before I switch to cht.

I would be thankful for some hints.

Hi,

1. your controlDict
2. fvSolution - was that all? Did you not use underrelaxation?

Regards Tobi

Hi Tobi,

thats my ControlDict
Even adding relaxation factors and increasing the number of correctors does not help
I have the feeling that the p_rgh BCs are not suitable. I am using fixedFluxPressure for all the wall boundaries similar to the hotRoom tutorial. I think I have to change the BC or even the patch type for the top face of the cylinder.

I've attached my case, if you want to have a look, it contains an Allrun script which does more or less everything such as initializing pressure and temperature by funkySetFiels. The parallel configuration uses 4 processors.

Maybe you can have a look.

Thanks

Bernhard

Hi Bernhard,

1. your mesh is not okay for a rotation 2d case:

This should be the first point to clear.


2. your pressure gradient is 100000 Pa to 89.4 Pa ? That makes no sence to me :)

3. why do you use a turbulence model :) ?

4. what fluid do you have?

Regards Tobi

Thanks for the hint, I actually thought that I have checked the mesh but I found a little mistake in the blockMeshDict, the z-points for all the arcs were wrong. After the change the mesh passed checkMesh, however, after executing ./Allrun and checking the mesh again I had the same error, so I removed the renumberMesh -overwrite and now it also passes checkMesh after Allrun. Is seems the renumberMesh is mixing something up.
Anyway that is the new mesh file. Attachment 27949
I have initialized p by calculating the hydrostatic pressure with the average density of my fluid (1822 kg m-3). The density is 1740 kg m-3 in the upper region of the storage at 550 °C and 1905 kg m-3 at 290°C in the lower region.
With a storage height of 2 m, I would expect a pressure of approx. 35750 Pa, what is shown in ParaView if I have a look at cell center values, if I switch to interpolated values, the values I have specified for the BCs of p are mixing that up. So I have manually set the values for top (0 Pa) and bottom (35747.64 Pa) and the value for the side wall to 0??? What should be the inital value for the side wall? I probably have to restructure my Allrun script and use $internalField.

My fluid is molten salt, it is a heat transfer fluid and heat storage medium used in solar thermal power plants, withstanding temperatures upto 600 °C.

I have found some literature for water storages at stand-by conditions and there was a good agreement between numerical simulation (with low-Re model) and experimental data. Thats more or less the reason for the turbulence model.

I will check if the changes in p have some effect...