[pgh]

Will limit temperature at end of each time step

limitFieldsTminmax
{
// Mandatory entries (unmodifiable)
type limitFields;
libs (fieldFunctionObjects);

// Mandatory entries (runtime modifiable)
fields (T);
limit both;
min 200;
max 1400;

// Optional (inherited) entries
writeControl outputTime;
}

[pgh]

Originally Posted by JM27

Hi redTo,

Unfortunately not, I am still struggling with this issue.

Perhaps you can elaborate more on your problem and maybe we can help each other out?

Use this functionObject:
https://www.openfoam.com/documentati...mitFields.html

limit temperature at end of each time step. Note , this will not go in fvOption but goes as functionObject



limitFieldsTminmax
{
// Mandatory entries (unmodifiable)
type limitFields;
libs (fieldFunctionObjects);

// Mandatory entries (runtime modifiable)
fields (T);
limit both;
min 200;
max 1400;

// Optional (inherited) entries
writeControl outputTime;
}

[aadil]

I have noticed that this problem doesn't occur for slower flows. Not sure what that means.

[lyu]

Hi JM27,

I'm having the same problem as you, Did you find some solution to it?
if you know how to solve this, pls let me know. thanks in advance

[JM27]

Hi,

Unfortunately the fvOptions method did not work out for me. The way I resolved this was by suppressing the solution of the temperature equation in the solver and solving for an isothermal flow, but of course whether this applies to your case depends on the physics of the problem you are trying to model. In my case it was fine to do this. Hope this helps.

[Bloerb]

This error often happens (atleast for me) with high speed flows and sudden changes in pressure. For example a water jet hitting a wall.


The main problem is the stiffness of the system caused by using an incompressible equationOfState. E.g systems like water (rhoConst) and air (perfectGas).
This is not an error in temperature but pressure, because pressure waves in incompressible media are instantaneous. You basically get cavitation, but the solver can't handle it.

Quick example...you have a pipe. Maybe even 1D. You initialize half of it with alpha=1 the rest with 0. Now you start the simulation with a high speed flow. Which yields a shockwave travelling in front of your fluid. This wave will reflect at the outlet with twice the speed and once it hits your fluid again you get pressure fluctuations. You also get shocks inside the initialized fluid region. Add a non perfect mesh to the mix and you basically can't solve this with compressibleInterFoam. Atleast not without serious work.

Hence you should always print the minimal pressure of your domain with a function object. Damping these waves with waveTransmissive boundary conditions can help. Changing rhoConst to something like adiabaticPerfectFluid with constants chosen in such a way, that the water becomes slightly compressible helps.
The same applies to the air phase. Adiabatic compression can be more stable. What also helps is uncorrected instead of corrected in fvSchemes (which removes non orthogonal corrections and hence reduces accuracy).

Having a perfect mesh helps a lot. You can basically solve everything on pure hex mesh. compressibleInterIsoFoam also helps

    fvScalarMatrix TEqn
    (
        fvm::ddt(rho, T) + fvm::div(rhoPhi, T) - fvm::Sp(contErr, T)
      - fvm::laplacian(turbulence.alphaEff(), T)
      + (
             mixture.totalInternalEnergy()
          ?
            fvc::div(fvc::absolute(phi, U), p)()() // - contErr/rho*p
          + (fvc::ddt(rho, K) + fvc::div(rhoPhi, K))()()
          - (U()&(fvModels.source(rho, U)&U)()) - contErr*K
          :
            p*fvc::div(fvc::absolute(phi, U))()()
        )
       *(
           alpha1()/mixture.thermo1().Cv()()
         + alpha2()/mixture.thermo2().Cv()()
        )

Different formulations of the temperature equation as in the foundation version, or different treatment of the continuity error did help slightly.


But...

Atleast from my perspective compressible multiphase simulations at high speeds with one incompressible medium is extremely difficult to keep stable and i haven't found any settings that do not eventually crash. As many have stated they are pretty accurate and converged until they fail. My suspicion is that a sudden pressure shocks cause this. Even with thousands of PIMPLE loops to converge every time step this fails. I have only managed to solve these so far by making the incompressible side more and more compressible. And hence loosing a lot of accuracy. And by creating absolutely perfect hex meshes.




The exact crash is in mixture.correctThermo() from twoPhaseMixtureThermo

[Gary Smith]

Originally Posted by Bloerb

This error often happens (atleast for me) with high speed flows and sudden changes in pressure. For example a water jet hitting a wall.


The main problem is the stiffness of the system caused by using an incompressible equationOfState. E.g systems like water (rhoConst) and air (perfectGas).
This is not an error in temperature but pressure, because pressure waves in incompressible media are instantaneous. You basically get cavitation, but the solver can't handle it.

Quick example...you have a pipe. Maybe even 1D. You initialize half of it with alpha=1 the rest with 0. Now you start the simulation with a high speed flow. Which yields a shockwave travelling in front of your fluid. This wave will reflect at the outlet with twice the speed and once it hits your fluid again you get pressure fluctuations. You also get shocks inside the initialized fluid region. Add a non perfect mesh to the mix and you basically can't solve this with compressibleInterFoam. Atleast not without serious work.

Hence you should always print the minimal pressure of your domain with a function object. Damping these waves with waveTransmissive boundary conditions can help. Changing rhoConst to something like adiabaticPerfectFluid with constants chosen in such a way, that the water becomes slightly compressible helps.
The same applies to the air phase. Adiabatic compression can be more stable. What also helps is uncorrected instead of corrected in fvSchemes (which removes non orthogonal corrections and hence reduces accuracy).

Having a perfect mesh helps a lot. You can basically solve everything on pure hex mesh. compressibleInterIsoFoam also helps

    fvScalarMatrix TEqn
    (
        fvm::ddt(rho, T) + fvm::div(rhoPhi, T) - fvm::Sp(contErr, T)
      - fvm::laplacian(turbulence.alphaEff(), T)
      + (
             mixture.totalInternalEnergy()
          ?
            fvc::div(fvc::absolute(phi, U), p)()() // - contErr/rho*p
          + (fvc::ddt(rho, K) + fvc::div(rhoPhi, K))()()
          - (U()&(fvModels.source(rho, U)&U)()) - contErr*K
          :
            p*fvc::div(fvc::absolute(phi, U))()()
        )
       *(
           alpha1()/mixture.thermo1().Cv()()
         + alpha2()/mixture.thermo2().Cv()()
        )

Different formulations of the temperature equation as in the foundation version, or different treatment of the continuity error did help slightly.


But...

Atleast from my perspective compressible multiphase simulations at high speeds with one incompressible medium is extremely difficult to keep stable and i haven't found any settings that do not eventually crash. As many have stated they are pretty accurate and converged until they fail. My suspicion is that a sudden pressure shocks cause this. Even with thousands of PIMPLE loops to converge every time step this fails. I have only managed to solve these so far by making the incompressible side more and more compressible. And hence loosing a lot of accuracy. And by creating absolutely perfect hex meshes.




The exact crash is in mixture.correctThermo() from twoPhaseMixtureThermo

hi, Bloerb
I think your answer is quietly good and gave me a lot of inspiration. but, still, I gat some issues with this compressibleInterFoam solvers, could you please help me?
I am doing a cavitation bubble collapsing simulation. I need to use the axisymmetric grid to reduce computation costs. Judging from the checkMesh file, my meshes are fine except for the max aspect ratio is relatively large, which is about 300. The simulation can be carried out successfully in the far-field. However, while the bubble is near the solid wall, the error of negative initial temperature T0 always occurs during the generation of high-speed flow or pressure sudden change. This is my basic case and here is my questions.
Q1:For the EOS of fluid, I have already used the adiabaticPerfectFluid, but I didn't find a proper EOS for air, I'm still using perfectGas. would you please tell me which EOS for air can become adiabatic and compression? that will help me a lot.
Q2:For p, I have used the waveTransmissive boundary condition, but for p_rgh I am not clear what boundary condition to use, especially for walls. And in many papers, only the boundary of P are mentioned. I've tried a lot of different boundary conditions, but this only changes when the negative temperature error occurs.
Hope you can help me, I was really suffered by this error. Thanks in advance.

[JM27]

Hi Gary,

Q1: You can use adiabaticPerfectFluid for the gas phase also, simply set the constant B to zero in thermophysicalProperties to obtain the polytropic EoS for the gas. Be careful with your terminology, both gas and liquid are fluids. I think you meant to say you are using adiabiaticPerfectFluid as the EoS for the liquid.

Q2: p is calculated from p_rgh in compressibleInterFoam and you can use the calculated BC for p if gravity is not important in your simulation (i.e. g=0, which is normally the case for cavitation simulation). Please refer to the OpenFOAM userguide to learn more about this BC. I wouldn't use waveTransmissive unless you can really understand what it does, I think this unnecessarily complicates things.

It would be nice to know how you get on with these changes. Goodluck!

[Gary Smith]

hi, JM27

Thanks for your reply. Actually, as you can see, I am a rookie of OpenFOAM, and I did mean to say is that for liquid I used adiabaticPerfectFluid as the EoS. Thanks for pointing this out.

And I've tried that you told me using the adiabaticPerfectFluid for gas as well, here are my parameters values for gas: p0=3540; rho0=0.025601; B=0; gamma=1.335;
for liquid: p0=1e+05; rho0=996; B=3.31e+8; gamma=7.15; is that OK?

I do the simulation both in the far-field and near the solid wall, unfortunately, they both failed because of the negative initial temperature error. But there is something different from my previous error. Since I monitor the max and min values for P and T in the computation domain, this time the Pmin didn't suddenly change and I can see the Tmin gradually decrease to a negative number over dozens of time steps.

In my previous error report, the minimum pressure, which I set to 100, appeared, and then after a few steps, the error was reported. This phenomenon often occurs near the collapsing time.

And I use waveTransmissive because is a non-reflect BC, and I also need to monitor the energy change, especially the wave energy, which is related to pressure, and results using this BC look pretty good.

By the way, I have been browsing a lot of your posts lately and benefited a lot, thanks to you for that.

[JM27]

Originally Posted by Gary Smith

Thanks for your reply. Actually, as you can see, I am a rookie of OpenFOAM, and I did mean to say is that for liquid I used adiabaticPerfectFluid as the EoS. Thanks for pointing this out.

The terminology you use has nothing to do with OpenFOAM but with fluid mechanics knowledge.

Originally Posted by Gary Smith

And I've tried that you told me using the adiabaticPerfectFluid for gas as well, here are my parameters values for gas: p0=3540; rho0=0.025601; B=0; gamma=1.335;
for liquid: p0=1e+05; rho0=996; B=3.31e+8; gamma=7.15; is that OK?

How are you determining your p0 and rho0 values inside the gas?

Originally Posted by Gary Smith

I do the simulation both in the far-field and near the solid wall, unfortunately, they both failed because of the negative initial temperature error. But there is something different from my previous error. Since I monitor the max and min values for P and T in the computation domain, this time the Pmin didn't suddenly change and I can see the Tmin gradually decrease to a negative number over dozens of time steps.

In my previous error report, the minimum pressure, which I set to 100, appeared, and then after a few steps, the error was reported. This phenomenon often occurs near the collapsing time.

And I use waveTransmissive because is a non-reflect BC, and I also need to monitor the energy change, especially the wave energy, which is related to pressure, and results using this BC look pretty good.

By the way, I have been browsing a lot of your posts lately and benefited a lot, thanks to you for that.

I assume that by far-field you mean in an infinite volume of liquid without any walls in the vicinity of the bubble?

The computation of p and T is related in the solver, so while the errors you see occuring are slightly different, it might be the same issue causing them. Some food for thought, since you are using an adiabatic EoS, then you no longer need to solve for temperature, unless you are interested in temperature changes. Worth thinking about the initial and boundary conditions you are setting for temperature and whether you need the temperature equation at all...

It is up to you if you want to use waveTransmissive for pressure, but if the solution is diverging, it might be worth trying a more "basic" BC for pressure, such as fixedValue of zeroGradient. The simpler the setup, the easier the detective work you need to do to figure out what is throwing the error will be

[Gary Smith]

The p0 and rho0 are what our research group has been using, and the gamma of the gas is from the paper, which I'm not sure is correct.

And yes, the far-field means that the bubble is in an unbounded flow field with no walls, just imagine a cube filled with water, and the cavitation bubble is in the middle. And near the solid wall is like a literal meaning.

Maybe I didn't mention it, which is my bad, but the settings of my work, like the choice of Eos: perfectGas for gas, adiabaticPerfectFluid for liquid and the BCs for 0 file, do not cause negative temperature errors in the far-field, the results of bubble radius with time had a great agreement with experimental results. But when the cavitation bubble is placed near the solid wall, there will be a negative temperature error.

Here is my BCs for 0 file: the BC for pressure field is chosen to be waveTransmissive with the atmospheric pressure p=101325Pa and a relaxation length scale lInf=0.0034. For the velocity U the BC pressureInletOutletVelocity is specified, for alha.water and temperature T, the zeroGradient BC is required. As for p_rgh, I've tried many different BCs, and now I'm using waveTransmissive, which I change the field to p_rgh. The above set of far-field BCs calculates pretty well. At the wall the pressure, temperature and volume fraction alpha.water were supplemented with zeroGradient, whereas the velocity used noslip. The p_rgh used fixedFluxPressure.

Of course, in the beginning, I used some basic BCs, such as fixedValue, and then changed it step by step to the above. It turns out that these settings performed well in the far-field. However, it seems like not work in the near-wall, so I want to know what is wrong and fix it.

And my work really needs the temperature equation, I know that some versions of OpenFOAM like OpenFOAM-extended-4.1, the compressibleInterFoam do not couple with TEqn, and this solver can fix this problem.

[Bloerb]

Example:

 
thermoType
{
    type            heRhoThermo;
    mixture         pureMixture;
    transport       const;
    thermo          eConst;
    equationOfState adiabaticPerfectFluid;
    specie          specie;
    energy          sensibleInternalEnergy;
}

mixture
{
    equationOfState
    {
        p0          1.0e5;
        rho0        1.0;
        gamma       1.4;
        B           0.0;
    }    
    .....
}

use a high B to make it basically incompressible. I should also note here, that the limitTemperature fvOption and limitFields functionObject do not stop this to my knowledge. It crashes inside TEqn:

TEqn.relax();      

fvOptions.constrain(TEqn);      

TEqn.solve();      

fvOptions.correct(T);      

mixture.correctThermo(); <-- it crashes here for me     

mixture.correct();

The negative temperature crashes are happening inside mixture.correctThermo() from twoPhaseMixtureThermo. limitTemperature fvOption / limitFields functionObject correct either fields T.alu / T.air instead of T or they limit the enthalpy/internal energy. Neither of these seem to work. Because they are not limiting the T value used in correctThermo.

I basically created my own fvOption for correct limiting. (The foundation version is probably doing it correctly, but you might have to remove the phase keyword)

    scalarField& Tc = he.primitiveFieldRef();
    
    forAll(cells_, i)
    {
        const label celli = cells_[i];
        Tc[celli] = max(min(Tc[celli], Tmax_), Tmin_);
    }
    
    volScalarField::Boundary& Tbf = he.boundaryFieldRef();

        forAll(Tbf, patchi)
        {
            fvPatchScalarField& Tp = Tbf[patchi];

            if (!Tp.fixesValue())
            {
                    forAll(Tp, facei)
                    {
                        Tp[facei] = max(min(Tp[facei], Tmax_), Tmin_);
                    }
            }
        }

which is nothing but a loop over all cells and boundaries to limit between min and max on the T field.
Now will this make your simulation run fine? No it simply prolongs the inevitable. Making a simulation run with limiters is always just spraying the dead cat with febreze. Might stop it stinking for a bit, but does not get rid of the underlying problem.

compressibleInterFoam/interFoam are not the correct solver for highly viscous flows (problems due the viscosity gradient at the interface), high density ratios and high speed shock wave like flows.


In many cases you get really high and unphysical velocities before the crash, your pressure is just bouncing like crazy and not like you'd expect etc.




@Garry


Try this(if you mesh has no non Orthogonality according to checkMesh)

laplacianSchemes
{
    default         Gauss linear orthogonal;
}
snGradSchemes
{
    default         orthogonal;
}

And increase pMin in thermoPhysicalProperties (i know counter productive for cavitation bubbles, but this has a huge impact). waveTransmissive is used to make outlets not reflect pressure waves, which they shouldn't. Hence in your case probably not what you primarily need if your boundary is very far away. If you see pressure waves reflecting from there it obviously helps, but this was merely mentioned by me to illustrate the point that this happens because the pressure solution is shit, which in turn leads to shitty temperatures, which when used to update the thermophysical properties crashes the simulation. And what i am describing isn't what crashes most peoples simulations. Most often the mesh is just pure garbage or the initialization/boundary conditions are shit.

[Gary Smith]

hi, Bloerb

Thanks for your reply and suggestions. After reading your reply, I have a better understanding of the reason for this negative temperature error.

As you said about the downside of the limiter, I wouldn't use the limiters unless I have to. And the limiter you created looks like it should work just fine.

And I replaced corrected with orthogonal in the above two. The results in the far-field are similar. Unfortunately, the error is still reported in the near wall.

To the accuracy of the simulation, I have to set pMin to a small value. As for my quality of the mesh, I guess my mesh isn't too bad, here is the log of checkMesh:

Overall domain bounding box (-0.03 0 -0.0013098283) (0 0.03 0.0013098283)
Mesh has 2 geometric (non-empty/wedge) directions (1 1 0)
Mesh has 3 solution (non-empty) directions (1 1 1)
Wedge front with angle 2.5 degrees
Wedge back with angle 2.5 degrees
All edges aligned with or perpendicular to non-empty directions.
Boundary openness (1.2927415e-17 1.0615266e-13 5.3768975e-13) OK.
Max cell openness = 2.259683e-16 OK.
Max aspect ratio = 301.6458 OK.
Minimum face area = 1.9404864e-12. Maximum face area = 2.7837102e-06. Face area magnitudes OK.
Min volume = 9.7556593e-18. Max volume = 2.7989816e-09. Total volume = 1.1788455e-06. Cell volumes OK.
Mesh non-orthogonality Max: 0 average: 0
Non-orthogonality check OK.
Face pyramids OK.
Max skewness = 0.33079647 OK.
Coupled point location match (average 0) OK.

Mesh OK.

[Bloerb]

interFoam like solvers are extremely susceptible to aspect ratios. The reason being that resolving a change in the interface across those aspect ratios is extremely difficult. I wouldn't recommend using aspect ratios like that. Especially if it crosses the interface. Otherwise the mesh is perfect. And orthogonal schemes are definitely the right choice, higher accuracy and better stability. Corrected is used to correct out high non orthogonality. Uncorrected or orthogonal have a big impact on stability. If it's just a blockMesh you could upload your case and i can check it for you.

[Gary Smith]

That would be wonderful. Here is my blockMeshDict. And the mesh is generated with the local refinement, so the high max aspect ratio seems like inevitable.

Another thing that you might need to note is I am using OpenFOAM-v9, and there are some differences compared with what you are using. In blockMeshDict you may need change [calc""] to [eval{}]. (just in case it may bother you). Anyway, I appreciate your kindful help.

[sumitzanje]

Originally Posted by JM27

Hi redTo,

Unfortunately not, I am still struggling with this issue.

Perhaps you can elaborate more on your problem and maybe we can help each other out?

Have you solved this issue? I'm also facing the same issue.

[CorbinMG]

Originally Posted by sumitzanje

Have you solved this issue? I'm also facing the same issue.

To place limits on the temperature field in compressibleInterFoam you can add the following to the end of the controlDict file to enable the limitFields function object:

functions
{
limitT
{
    type            limitFields;
    libs            (fieldFunctionObjects);

    fields          (T);
    limit           both;
    min             299;
    max             5000;
    log             false;
    writeControl    none;
}
}

This method worked for me, while the fvOptions method did not.
Note: I am using OpenFOAM v2206

Good luck!
Best regards,
Corbin

[mohammedtareq]

Hi JM27
did you figure out the solution of this problem becouse i facees same problem with nagetive temprature

[shock77]

Depending on your solver, you might have to implement it yourself when fvOptions is not included. Some years ago I did that for rhoCentralFoam. Its not too much work: Sudden Increase of residuals