[JM27]

Hi there,


I am running some simulations in compressibleInterFoam and I want to make the solver more robust by preventing the temperature from becoming negative.



After some research on how this can be done, I have tried using the elegant fvOptions function to limit the temperature in the domain. I am using OpenFOAM v5.



The fvOptions file I set is as follows (values are arbitrary, just want to test functionality with compressibleInterFoam):

 limitT
    {
        type                  limitTemperature;
        active               yes;     
        selectionMode   all;   
        min                   200;            
        max                  20000;   
    }

Upon trying to run this I get the following error:

Selecting finite volume options model type limitTemperature
    Source: limitT
    - selecting all cells
    - selected 571536 cell(s) with volume 9e-07
[0] 
[0] 
[0] --> FOAM FATAL ERROR: 
[0] Not implemented
[0] 
[0]     From function const Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> &Foam::twoPhaseMixtureThermo::he() const
[0]     in file twoPhaseMixtureThermo.H at line 138.
[0] 
FOAM parallel run aborting
[0] 
[2] 
[2] 
[2] --> FOAM FATAL ERROR: 
[2] Not implemented
[2] 
[2]     From function const Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> &Foam::twoPhaseMixtureThermo::he() const
[2]     in file twoPhaseMixtureThermo.H at line 138.
[2] 
FOAM parallel run aborting
[2] 
[1] 
[1] 
[1] --> FOAM FATAL ERROR: 
[1] Not implemented
[1] 
[1]     From function [0] #0  Foam::error::printStack(Foam::Ostream&)const Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> &Foam::twoPhaseMixtureThermo::he() const
[1]     in file twoPhaseMixtureThermo.H at line 138.
[1] 
FOAM parallel run aborting
[1] 
[1] #0  Foam::error::printStack(Foam::Ostream&)[2] #0  Foam::error::printStack(Foam::Ostream&)[3] 
[3] 
[3] --> FOAM FATAL ERROR: 
[3] Not implemented
[3] 
[3]     From function const Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> &Foam::twoPhaseMixtureThermo::he() const
[3]     in file twoPhaseMixtureThermo.H at line 138.
[3] 
FOAM parallel run aborting
[3] 
[3] #0  Foam::error::printStack(Foam::Ostream&) at ??:?

I understand that the error is in the twoPhaseMixtureThermo however I cannot seem identify the cause of the problem yet. The phases involved are air and water. I use perfectGas (compressible) to model air and rhoConst for water (incompressible)...



Could this be the cause of the problem ? Does limitTemperature function not work with incompressible phases?


I will keep on looking into this and post the solution if found, however any help will be greatly appreciated!



Thanks

**************************************

Note:

1. I tried running the case with fvOptions specified in /constant or in /system, however this does not seem to make a difference (as expected).



2. I have been through some tutorials to understand the fvOptions functionality however there are none available on limitTemperature.


3. I have double-checked the pEqn.H, UEqn.H and TEqn.H files in compressibleInterFoam allow for the fvOptions functionality.

[JM27]

Update:

Tried running the same case but setting the equation of state for the liquid to perfectFluid instead of rhoConst. The aim here was to test whether the error arises due to the incompressibility of the liquid.

This fails, and the same error posted above persists.

Should I include the following in the make/options file of compressibleInterFoam.C?

"-I$(LIB_SRC)/fvOptions/lnInclude" ?

[JM27]

Still stuck on this problem - can anyone help, please?

[nhlstd]

Originally Posted by JM27

Still stuck on this problem - can anyone help, please?

Hi,

Did you find the solution?
I still stuck on this problem as you.

[raumpolizei]

Hello there,


the reason for that is (as the ouput states) that the function he() of the class twoPhaseMixtureThermo is not implemented. If you take a look at the documentation (of18.12, https://www.openfoam.com/documentati...ce.html#l00131) and check the function he(), you will understand why (also, see the code below). The temperature limitation seem to be achieved by changing the transported energy variable (the variable for which a transport equation is solved in your solver). Limiting the temperature therefore relies on this function to be available. For multiphase mixtures this is not working as the energy variable is different for each phase and no "generalized" energy variable is available. Even if it was available, it would be impossible to change this energy variable yielding same limits in T for both phases without further information. This is due to for instance different heat capacities. A solution could be to implement your own version of the limiting temperature feature which handles twophasemixturethermo by accessing the energy variable via thermo1().he() and thermo2().he(). If however you are only interested in one phase, you could just outcomment the NotImplemented macro in the function (which I do not recommend).

// from the twoPhaseMixtureThermo class
virtual volScalarField& he()
{
   NotImplemented; // this is causing the error, someone put it there, reason why before commenting it.
   return thermo1_->he();
}

Good luck!
RP

[JM27]

Hi RP,

thank you for your reply and sorry for not answering before, but I just noticed this message for some reason.

I'm not sure I fully understand why the temperature goes negative in the first place. I use an oscillatory pressure boundary, which seems to be fine if I choose an amplitude that does not go below atmospheric pressure i.e. p_oscillating < p_atmospheric.

However, if I do chose p_oscillating to be >= p_atmospheric, it causes problems. The 'negative' pressure will cause the temperature in the domain to go negative accordingly. I sort of manage to bypass this by altering the value of pMin in thermophysicalProperties dict, but I am not sure whether there is a better way of doing this and I do wish to make the solver more robust against this issue.

You seem to be quite experienced in compressibleInterFoam, so if you have any suggestions they would be much appreciated.

Thanks

JM

[The King]

Hi,

did you find a solution? Got the same problem...

[arsimons]

Hello everyone

I think I might help a little with this problem.


First of, I think that a negative pressure should never be possible since it is physically impossible (just like a negative temperature). As such, the oscillation of the pressure has no limitation as long as the amplitude is smaller than the mean pressure at this boundary.


Second, I found a way to use the limitTemperature in multiphase. I have simply added 'phase' in the 'limitT' part of fvOptions. It looks like this now (my present phases are 'air' and 'water'):


limitT
{
type limitTemperature;
active yes;

min 101;
max 1000;
selectionMode all;
phase air;
}


HOWEVER
I do still have a problem when I try to run this in parallel. In that case, it seems to ignore the fvOptions (it errors at the same time with the same negative value for the temperature as when fvOptions was not present) even though it is certainly read.


Does anyone have experience with limitTemperature in parallel run?


Thanks in advance.



Best regards
Arne

[EngMec]

Hi JM27:
I also used compressibleInterFoam and experienced the same problem with temperature limiter. My case worked fine with v6 and v7 but this error just keeps poping up on v5. so I think it is some kind of bug that they later fixed.

[arsimons]

Originally Posted by EngMec

Hi JM27:
I also used compressibleInterFoam and experienced the same problem with temperature limiter. My case worked fine with v6 and v7 but this error just keeps poping up on v5. so I think it is some kind of bug that they later fixed.

Are you sure it was fixed? Since someone told me that temperature limiter can never work in compressibleInterFoam since it works with temperature rather than enthalpy and the temperature limiter apparently works on enthalpy? (I am sorry if I am wrong but I have not checked this myself)

[EngMec]

I noticed the post by [raumpolizei]. I don't known how they did it but this error disappears and my case stops exploding. So, I guess it works, somehow. Here is the limiter I used on a cluster with v7.

limitTAir
{
    type		limitTemperature;
    active		yes;
	
    selectionMode all;
    min			50;
    max			350;
    phase		air;
}


limitTWater
{
    type		limitTemperature;
    active		yes;
	
    selectionMode all;
    min			295;
    max			305;
    phase		water;

}

[JM27]

Hi EngMec,

I also tried using two separate limitTemperature fvOptions, one for air and one for water.

I also monitor the maximum and minimum values of temperature in the solution and I observe some very strange behaviour. Although fvOptions is read at the beginning of the simulation, the maximum temperature in the domain still exceeds the maximum I set in fvOptions (the simulation almost reaches the end).

The solution eventually diverges due to a very low negative initial temperature, which also falls out of bounds of the limitTemperature...

I am starting to think that this is a bug? Or fvOptions does not works with multiphase? The OF version I use now is v6.

Someone please help

[The King]

Hi JM27,

I have been struggling with this negative temperature for some time. An improved mesh (courser) where I increased the thickness of the first wall layer solved it.

I started to save small time steps (0.01s) and could see that the solution start to be strange at my outlet with some back flow and than the temperature start dropping. If you post your settings and BC's, I can take a look as well but I don't know if I will find something.

Good luck (and keep distance!)

[JM27]

Thanks for the reply.

May I ask what boundary conditions you are using at the outlet?

I am simulating internal flow and I suspect that this could be an effect of the boundary condition at the outlet. I read that fixedValue may cause some flow distortion at outlet and was thinking of trying fixedMeanValue instead (or similar)...

Here are my boundary conditions. The flow is an internal flow in a gap with two parallel walls and an outlet (right BC):

Valid fields:
    volScalarField	alpha.water
    volVectorField	U
    volScalarField	p_rgh
    volScalarField	p
    volScalarField	T

group	: wall
    scalar		alpha.water zeroGradient
    scalar		p_rgh	  fixedFluxPressure
    scalar		p		  calculated
    scalar		T		  zeroGradient
    vector		U		  noSlip

patch	: right //outlet
    scalar		alpha.water zeroGradient
    scalar		p_rgh	  fixedValue
    scalar		p		  calculated
    scalar		T		  zeroGradient
    vector		U		  zeroGradient

empty	: left // axis
    scalar		alpha.water empty
    scalar		p_rgh	  empty
    scalar		p		  empty
    scalar		T	 	  empty
    vector		U		  empty

group	: wedge
    scalar		alpha.water wedge
    scalar		p_rgh	  wedge
    scalar		p		  wedge
    scalar		T		  wedge
    vector		U		  wedge

With regards to changing my mesh, this is not really an option as I am trying to replicate an existing study, especially the first layer thickness.

What about the pMin value?

[JM27]

Hello everyone,

Unfortunately I am still having issues with this case

The fvOptions limitTemperature function improves the solution as it runs almost entirely. However, the solution eventually crashes due to an initial negative temperature even though the fvOptions is being used!!

I've checked the log file, and can see that the fvOptions function is being read and created.

No MRF models present

Creating finite volume options from "constant/fvOptions"

Selecting finite volume options model type limitTemperature
    Source: limitT_AIR
    - selecting all cells
    - selected 1920000 cell(s) with volume 8.1708508826e-10
Selecting finite volume options model type limitTemperature
    Source: limitT_WATER
    - selecting all cells
    - selected 1920000 cell(s) with volume 8.1708508826e-10
Courant Number mean: 0 max: 0

However, almost half way through the simulation I start to note that the maximum velocity in the domain increases unrealistically, sometimes even reaching 1000m/s.
Similarly, when this happens the temperature starts to decrease below the minimum value set in fvOptions. I set minimum value to be min T = 285K in fvOptions, but when the velocity starts to increase the minimum T in the domain decreases rapidly

The strange thing is, that I am comparing my results to existing cases in literature and the agreement of the alpha field is excellent for 90% of the simulation, after which the simulation crashes.

I have tested several boundary conditions and compared my setup with tutorials for compressibleInterFoam, but I am still unable to identify what is causing this divergence problem.

I am now testing various discretisation schemes as I suspect that this might be causing the solution to lose its boundedness....

Please help

[raumpolizei]

Indeed, this does not sound healthy. How is the overall mesh quality? What schemes are you using, is the simulation converging at each timestep? If yes what are your convergence criteria? These are the questions that could lead you to a solution of your problem. Good luck!
RP

PS: Just my personal opinion (I do not want to harm anybody here): I think limiters are a bad solution most of the time since they tend to overshadow other problems (crappy mesh or wrong schemes). If you are not simulating extreme physics like supercritical flows or something alike, you should be trying to work without them. Bounding scalars goes in the ultima ratio direction.

[JM27]

Hi Raumpolizei,


First of all thanks for the reply.


Here is my checkMesh:

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create polyMesh for time = 0

Enabling all (cell, face, edge, point) topology checks.

Enabling all geometry checks.

Time = 0

Mesh stats
    points:           3847001
    internal points:  0
    edges:            9610200
    internal edges:   1916201
    internal edges using one boundary point:   0
    internal edges using two boundary points:  1916201
    faces:            7683200
    internal faces:   3836200
    cells:            1920000
    faces per cell:   5.9996875
    boundary patches: 6
    point zones:      0
    face zones:       0
    cell zones:       0

Overall number of cells of each type:
    hexahedra:     1919400
    prisms:        600
    wedges:        0
    pyramids:      0
    tet wedges:    0
    tetrahedra:    0
    polyhedra:     0

Checking topology...
    Boundary definition OK.
    Cell to face addressing OK.
    Point usage OK.
    Upper triangular ordering OK.
    Face vertices OK.
    Topological cell zip-up check OK.
    Face-face connectivity OK.
  <<Writing 4 cells with two non-boundary faces to set twoInternalFacesCells
    Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces...
    Patch               Faces    Points   Surface topology                   Bounding box
    top                 3200     6401     ok (non-closed singly connected)   (0 -0.000218096936827 0.000375) (0.00499524110791 0.000218096936827 0.000375)
    bottom              3200     6401     ok (non-closed singly connected)   (0 -0.000218096936827 -0.000375) (0.00499524110791 0.000218096936827 -0.000375)
    right               600      1202     ok (non-closed singly connected)   (0.00499524110791 -0.000218096936827 -0.000375) (0.00499524110791 0.000218096936827 0.000375)
    left                0        0        ok (empty)                        
    front               1920000  1923801  ok (non-closed singly connected)   (0 0 -0.000375) (0.00499524110791 0.000218096936827 0.000375)
    back                1920000  1923801  ok (non-closed singly connected)   (0 -0.000218096936827 -0.000375) (0.00499524110791 0 0.000375)

Checking geometry...
    Overall domain bounding box (0 -0.000218096936827 -0.000375) (0.00499524110791 0.000218096936827 0.000375)
    Mesh has 2 geometric (non-empty/wedge) directions (1 0 1)
    Mesh has 3 solution (non-empty) directions (1 1 1)
    Wedge front with angle 2.49999996192 degrees
    Wedge back with angle 2.49999996192 degrees
    All edges aligned with or perpendicular to non-empty directions.
    Boundary openness (-1.49847870122e-13 1.28518123348e-13 -9.9480090689e-17) OK.
    Max cell openness = 2.48355969881e-16 OK.
    Max aspect ratio = 31.0351835932 OK.
    Minimum face area = 6.85510132929e-15. Maximum face area = 6.85510130749e-10.  Face area magnitudes OK.
    Min volume = 5.35045061858e-21. Max volume = 1.06992292132e-15.  Total volume = 8.1708508826e-10.  Cell volumes OK.
    Mesh non-orthogonality Max: 9.69700457025e-06 average: 0
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 0.330796465394 OK.
    Coupled point location match (average 0) OK.
    Face tets OK.
    Min/max edge length = 5.0290308e-08 0.000436193873654 OK.
    All angles in faces OK.
    Face flatness (1 = flat, 0 = butterfly) : min = 1  average = 1
    All face flatness OK.
    Cell determinant (wellposedness) : minimum: 0.0106786000355 average: 3.31472668332
    Cell determinant check OK.
    Concave cell check OK.
    Face interpolation weight : minimum: 0.375 average: 0.49844202787
    Face interpolation weight check OK.
    Face volume ratio : minimum: 0.333333333332 average: 0.993108806431
    Face volume ratio check OK.

Mesh OK.

End

and my discretisation schemes:

ddtSchemes 
{
   default                  Euler;                           // first order in time
} 

gradSchemes
{
    default             Gauss linear;                
}

divSchemes // Convection Schemes Settings
{
 
   default                                            none;    
   div(phi,alpha)                                  Gauss vanLeer; 
   div(rhoPhi,U)                                   Gauss linearUpwind grad(U);  // UEqn
   div(((rho*nuEff)*dev2(T(grad(U)))))    Gauss linear; // vanLeer;
   div(rhoPhi,T)                                    Gauss upwind;    // TEqn
   div(rhoPhi,p)                                    Gauss linear;
   div(phi,p)                                        Gauss linear;
   div(rhoPhi,K)                                    Gauss linear;    // TEqn     
   div(phirb,alpha)                                Gauss interfaceCompression 1; // alphaEqn
}

laplacianSchemes
{
   default              Gauss linear corrected;
}

interpolationSchemes
// interpolating cell-centred values to face values
{
   default                      linear;                              
}
snGradSchemes
{
     default                  corrected;
}


// ************************************************************************* //

The overall checkMesh is OK, but as highlighted in bold, the aspect ratio of some cells is quite high ~31. These are cells in the near-wall region were I refine the mesh to resolve boundary layer flow.


I think the divSchemes are causing problem, specifically div(phi,p) and div(rhoPhi,p) as the linear scheme is essentially central-differencing scheme and is not bounded.

Perhaps I will try upwind for p and see whether there is any improvement.

Any thoughts?

[JM27]

I am still struggling with this problem, can someone help me please?

I have since tried to avoid using the limitTemperature in fvOptions by suppressing solution of the temperature equation in compressibleInterFoam. This makes the solution stable (well at least it's working so far, just not sure if it's a healthy hack to introduce), so I guess the problem is arising due to the temperature not converging (~1000 iterations) in every time-step.
I'm not really interested in the solution of temperature so I guess this is OK, or rather better than limiting the temperature to some range...

Does anyone have any suggestions with regards to suitable schemes for this problem ?

[redTo]

Hi JM27,

I'm having the same problem than you Did you find some solution to it?

Cheers!

[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?