[peob]

I am encountering a peculiar problem in a number of my simulations with pimpleFoam and rhoPimpleFoam. When I run my cases at reasonably large Courant numbers (maxCo from 50 to 10000, say) they are nicely stable, but I don't have the time-accuracy I would like. So I lower the maxCo down to something on the order of 1 or smaller (e.g. 0.3) to improve my time-accuracy, however, my solution then goes unstable. This is the exact opposite behavior from what I would expect: stable at small maxCo and unstable at large maxCo.

Attached is an example pimpleFoam case that demonstrates the behavior. Unfortunately I cannot include the grid, because it exceeds the file size limitations. However, if you happen to have cfMesh on your system, all the files are set up for you to simply run "cartesianMesh" and you should then be able to generate the grid. NOTE... you will also have to scale the grid from millimeters to meters before running the case. Simply execute the command: transformPoints -scale '(0.001 0.001 0.001)'. Boundary conditions and all of the other settings are already there; no need to change them.

If you set the maxCo to a large value (say 100) then the solution converges just fine. However, if you set the maxCo to a small value (try 1.0, or 0.5 or smaller) the solution goes unstable. It may take a few hundred iterations to see the instability grow, but as the solution continues to run, the instability worsens. I have defined a couple of animation slices (y and z-direction) where you can plot the progression of the instability in the U-field.

I'm running OpenFOAM version 3.0.1 on Ubuntu 14.04 LTS. In the attached case, I used cfMesh to generate the cartesian mesh. CheckMesh indicates that I have a good quality mesh.

Has anyone else encountered this problem? I welcome any thoughts or suggestions.

Thanks.
Phil

[tas38]

Originally Posted by peob

I am encountering a peculiar problem in a number of my simulations with pimpleFoam and rhoPimpleFoam. When I run my cases at reasonably large Courant numbers (maxCo from 50 to 10000, say) they are nicely stable, but I don't have the time-accuracy I would like. So I lower the maxCo down to something on the order of 1 or smaller (e.g. 0.3) to improve my time-accuracy, however, my solution then goes unstable. This is the exact opposite behavior from what I would expect: stable at small maxCo and unstable at large maxCo.

Attached is an example pimpleFoam case that demonstrates the behavior. Unfortunately I cannot include the grid, because it exceeds the file size limitations. However, if you happen to have cfMesh on your system, all the files are set up for you to simply run "cartesianMesh" and you should then be able to generate the grid. NOTE... you will also have to scale the grid from millimeters to meters before running the case. Simply execute the command: transformPoints -scale '(0.001 0.001 0.001)'. Boundary conditions and all of the other settings are already there; no need to change them.

If you set the maxCo to a large value (say 100) then the solution converges just fine. However, if you set the maxCo to a small value (try 1.0, or 0.5 or smaller) the solution goes unstable. It may take a few hundred iterations to see the instability grow, but as the solution continues to run, the instability worsens. I have defined a couple of animation slices (y and z-direction) where you can plot the progression of the instability in the U-field.

I'm running OpenFOAM version 3.0.1 on Ubuntu 14.04 LTS. In the attached case, I used cfMesh to generate the cartesian mesh. CheckMesh indicates that I have a good quality mesh.

Has anyone else encountered this problem? I welcome any thoughts or suggestions.

Thanks.
Phil

I agree with you, this behavior is counter-intuitive.

A couple of things...
1. Do you see the same behavior if you turn off the adaptive time stepping and simply increase dt?
2. nOuterCorrectors = 2, so you likely do not meet the residual values you set to terminate your PIMPLE loop, i.e. the code is running 'nearly' in PISO mode. You may want to change this to a larger value so your residuals are met within each time step's sub-loop to take advantage of the solver's PIMPLEness.
3. Since you are running in nearly PISO mode, do you see any improvement by increasing the nCorrectors?

[peob]

Troy,

Thanks for the reply. I tried your suggestions and here is what I found...

1) I do see the same behavior when I use a fixed time-step. I ran time-step sizes that yielded maxCo values in the range of 0.5 to 1.0
2) I ran nOuterCorrectors in the range of 2 to 40 and saw no change in behavior
3) I ran nCorrectors in the range of 1 to 20 and saw no change in behavior

I also ran in just PISO mode (nOuterCorrectors=1) with nCorrectors in the range of 2 to 20 and saw no change in behavior.

I should also mention that I took this case and converted it over to run in Fluent. I ran in Fluent with time-step sizes that yielded maxCo values in the range of 0.05 to 10, and the solution was stable in all cases. So, it looks like there may be something peculiar going on with OpenFOAM.

I welcome any further suggestions. I think that I'll try a few more ideas, and if nothing improves, then I'll submit this to the mantis bug report.

Thanks again for your suggestions.

[tas38]

Could you provide an external link to the mesh? I would like to try and run the code.

[peob]

I'll see what I can do. I won't be able to get you the grid until at least Tuesday, May 31st. Thanks for taking the time to look into this.

[peob]

Troy,

Have you had a chance to look at my sample case yet?

Thanks.
Phil

[peob]

If anyone still happens to be following this thread....

I submitted this problem to the mantis bug reporting site, and the answer I received was that it was a "user support issue" and not a bug, because my grid wasn't the best.

Well, I haven't given up on this yet, because being able to run reasonable time-steps for time accuracy on "industrial-type" grids is pretty important to my livelihood. My grids are typically on the order of 10-50 million cells. The one I submitted to the bug reporting site was a mere 9978 cells, because any bigger and it would be too big to download. As a consequence of the low cell count, it wasn't the best of grids, even though it did pass checkMesh. I understand this.

I have since generated some bigger and better quality meshes (using cfMesh's cartesianMesh generator). These grids range from 87172 cells to about 2.6 million cells. I have found the same troubling behavior (unstable at small time-step sizes, but stable at large time-step sizes) when running pimpleFoam using OpenFOAM version 3.0.1.

A couple of interesting points...

1) when I run the same case using interFoam (version 3.0.1) at small time-steps (maxCo ~0.1) the solution is stable (unlike when I run pimpleFoam). My geometry has one inlet and one outlet. I initialize the inlet region with liquid and let it run at a large maxCo until the liquid fills the domain (stable). Then I switch over to maxCo =0.1 and run it another 20,000 time steps, and the solution remains stable.

2) I also tried running pimpleFoam using foam-extend version 3.2, with the same exact case as OpenFOAM version 3.0.1, and the solution ran stable at maxCo = 0.1, whereas OpenFOAM version 3.0.1 diverged.

I have tried different "gradSchemes", "snGradSchemes", "divSchemes", and "lapacianSchemes". I have tried running in PISO mode and PIMPLE mode with nCorrectors from 1 to 20, nOuterCorrectors from 1 to 40, nonOrthogonalCorrectors from 1 to 10. I have also adjusted the relaxation factors. Still I encounter the instability.

Enough for now. I'll keep updating as I find out anything new.

If anyone else has experienced a similar phenomenon, please let me know. Maybe I'm the only one who has encountered this problem?

Thanks for reading.

[hoemmaria]

Hi, Phil.

I am experiencing the same problem when using pimpleDyMFoam. I am sorry to say that I do not know anything about why this happens, so I am wondering if you have gained any more knowledge?

Kind regards,
Maria

[peob]

Maria,

I have not been able to gain any more knowledge, because I had to move on to other analyses, which do not require time-accuracy, so I've been running cases at the larger Courant Numbers (Co>20).

I have been debating whether I should resubmit this problem to the mantis bug report, but I don't want to use up any goodwill I may have left, since the reply to my first posting was that this wasn't a bug, but a user error. I do, however, find it interesting that the solution is unstable with pimpleFoam, yet stable if I run interFoam with only a single (liquid) phase. I would think that this might point to something. Unfortunately, I am an old-timer from the days of Fortran-77, and I am struggling to learn enough about C++ for me to be able to track down the problem myself. I'll keep working at it, but am not too optimistic that I'll get there, which is why I have asked for help.

Regarding your problem using pimpleDyMFoam, can you run your case with pimpleFoam (no DyM)? If so, does it go unstable at very small time-steps (very small Co, e.g. Co<0.5)? If it does, can you try setting up your case to run with interFoam, but with just a single (liquid?) phase, so see if that is stable, like my case was?

BTW, at what Courant Number (and time-step size) does your case start to go unstable? The trouble is that once the instability has progressed, the code automatically changes the time-step size to try to accommodate the Courant Number, and the time-step size can get very, very small.

[hoemmaria]

Phil,

I made my case work by playing around with different values of the solver correctors - nNonOrthogonalCorrectors, nCorrectors and nOuterCorrectors. Thus, I am not really sure I had the same problem as you. I am not very familiar with C++ myself, so I do not really know everything that is going on with the code.

I have a max Co=0.5 in my simulation and only one phase.

Kind regards,
Maria

[peob]

In my earlier post, I was using OpenFOAM v3.0.1, but have continued to encounter this problem with versions 4.0, 4.1, 5.0, 6.0. (I haven't compiled 7.0 yet, but I suspect the problem will persist)

I'm pretty confident that there is something not quite right with the algorithm implementation. I find it peculiar that I can run exactly the same problem on exactly the same grid and with pimpleFoam the solution is unstable, but with interFoam the solution is stable. Seems like someone should be able to find out what is different between the implementation of these two algorithms, and identify the problem. I'm an old-timer FORTRAN guy, and don't know how to read C++ code, so I need some help in this area.

[Tobermory]

Folks - I have not reviewed your cases in detail, but let me offer the following observation to see if it helps.

Reducing the time step significantly will reduce your temporal differencing errors ... if this is first order temporal differencing then you are removing some of the numerical stabilising (damping) that is present in your high Co runs. That is probably why your low Co runs are much more lively.

You see the same effect when refining mesh spacing - what runs nice and stably with a coarse mesh can go haywire as you refine the mesh and remove the numerical diffusion that is stabilising the solution.

[jdw135]

Sorry to necro a dead thread here, but I'm currently experiencing this same problem. Running a simple shear layer domain, 2 inlets (high and low speed), one outlet. Case runs smoothly at Co 0.5-1 with a uniform grid and cell size = 5mm. If I drop to 1mm or less cell size, I get negative initial temperature errors in large regions of the grid, and using an fvOptions to clamp temperature to a reasonable range results in non-physical results. Anyone have more info?

I'm on 20.06 using rhoPimpleFoam.

[Tobermory]

Difficult, Jonathan to diagnose with so little info, but from what you've described, I wonder whether your initial fields are specified correctly?

It seems that you need the added numerical diffusion from a coarser mesh and time step to smear out some problems in the initial time steps. This diffusion disappears of course when you make the mesh 5x smaller (I am assuming you drop the time step by 5x as well to preserve the same Co?).

Take a look at your initial fields - make sure you specify everything consistently, initial fields and boundaries; don't set zero initial fields for quantities like nut; try find out why your temperature field is going nuts - wrong thermoProps? wrong BC? wrong initial field?

Good luck - it can be really frustrating to track these issues down, but work through it methodically - and when you find the answer, write it down so you know next time! Don't rely on memory - trust me!

[jdw135]

I will admit that I am extremely new to openfoam, so I'm sure there are issues with my initial fields. For instance, is there a resource that discusses how to properly set initial fields for things like nut? Also, can I set an initial total temperature field?

I did change my time step as I changed the grid spacing to ensure I kept the same Co, but I still ran into issues. I'm almost certain my problem lies in lack of knowledge, but resources for openfoam are spread so vastly across the internet, and across so many versions, that finding reliable documentation and tutorials has been quite difficult.

[Tobermory]

Agreed - it can be pretty daunting, starting out - and no, I haven't ever found a good primer to walk through the basics. ... sorry! The good news is that it becomes easier with time ... true of anything I guess.

For nut, you can calculate this from your turbulence model. For example, if it's a k-epsilon derivative then nut is defined as nut = Cmu * k^2 / epsilon, with Cmu=0.09. So, assuming that you are initialising the domain with some initial values of k0, epsilon0, you can calculate a value of nut0 and use that to initialise the nut field. This provides some turbulent diffusion in the first few iterations, and stabilises things a little. As for k0 & epsilon0, assume a low turbulence intensity eg 1%, and a suitable lengthscale, e.g. the diameter of your inlet. You also want to initialise alphat - use nut0 and the Prandtl number to calculate that.

I am not sure about total temperature - I assume that you are solving an equation for T, so you need to supply an initial T field that is in balance with your U and p fields ... you can specify total T on the boundaries, but again make sure that it is in balance with everything else.

Hope that helps - and good luck. If you are still having problems, you can always try posting the case to see if someone has the time to debug it for you.

[NablaDyn]

In my experience, negative temperature might also, i.a., occur in case of bad mesh, nonsensical mesh dimensions (e.g. defined in mm instead of m) or inappropriate thermophysicalProperties settings, boundary conditions.

[jdw135]

Originally Posted by Tobermory

Agreed - it can be pretty daunting, starting out - and no, I haven't ever found a good primer to walk through the basics. ... sorry! The good news is that it becomes easier with time ... true of anything I guess.

For nut, you can calculate this from your turbulence model. For example, if it's a k-epsilon derivative then nut is defined as nut = Cmu * k^2 / epsilon, with Cmu=0.09. So, assuming that you are initialising the domain with some initial values of k0, epsilon0, you can calculate a value of nut0 and use that to initialise the nut field. This provides some turbulent diffusion in the first few iterations, and stabilises things a little. As for k0 & epsilon0, assume a low turbulence intensity eg 1%, and a suitable lengthscale, e.g. the diameter of your inlet. You also want to initialise alphat - use nut0 and the Prandtl number to calculate that.

I am not sure about total temperature - I assume that you are solving an equation for T, so you need to supply an initial T field that is in balance with your U and p fields ... you can specify total T on the boundaries, but again make sure that it is in balance with everything else.

Hope that helps - and good luck. If you are still having problems, you can always try posting the case to see if someone has the time to debug it for you.

Thanks so much of the information and help. I will work through these things today.

In my experience, negative temperature might also, i.a., occur in case of bad mesh, nonsensical mesh dimensions (e.g. defined in mm instead of m) or inappropriate thermophysicalProperties settings, boundary conditions.

My grid is square, with a uniformly-spaced structured grid and an aspect ratio of 1 everywhere. I am fairly certain that mesh quality is at the bottom of the list of suspects in this particular case.

[jdw135]

I have to say though, these bugs are getting frustrating. I'm now being told 'mixingLength' is not found in my epsilon BC file, but it's clearly there. Snippet posted below.

    lsinlet
    {
        type            turbulentMixingLengthDissipationRateInlet;
        mixingLength    0.7;    // 0.7 inlet diameter
        value           $internalField;
    }