
[Sponsors] 
fvOptions limitTemperature crashing in compressibleInterFoam 

LinkBack  Thread Tools  Search this Thread  Display Modes 
April 2, 2019, 13:00 
fvOptions limitTemperature crashing in compressibleInterFoam

#1 
Member
Join Date: Apr 2018
Location: UK
Posts: 65
Rep Power: 5 
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): Code:
limitT { type limitTemperature; active yes; selectionMode all; min 200; max 20000; } Upon trying to run this I get the following error: Code:
Selecting finite volume options model type limitTemperature Source: limitT  selecting all cells  selected 571536 cell(s) with volume 9e07 [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 doublechecked the pEqn.H, UEqn.H and TEqn.H files in compressibleInterFoam allow for the fvOptions functionality. 

April 2, 2019, 13:10 
update: perfectFluid fails too

#2 
Member
Join Date: Apr 2018
Location: UK
Posts: 65
Rep Power: 5 
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" ? 

April 10, 2019, 11:22 
help

#3 
Member
Join Date: Apr 2018
Location: UK
Posts: 65
Rep Power: 5 
Still stuck on this problem  can anyone help, please?


May 4, 2019, 10:00 

#4 
New Member
nahaleh sotoudeh
Join Date: Jun 2017
Posts: 3
Rep Power: 6 

May 4, 2019, 12:08 

#5 
Member
Join Date: Dec 2018
Location: Darmstadt, Germany
Posts: 87
Rep Power: 4 
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). Code:
// 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(); } RP 

July 5, 2019, 07:55 

#6 
Member
Join Date: Apr 2018
Location: UK
Posts: 65
Rep Power: 5 
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 

September 13, 2019, 06:37 

#7 
Member
Arnout
Join Date: Nov 2010
Posts: 46
Rep Power: 12 
Hi,
did you find a solution? Got the same problem... 

January 10, 2020, 04:20 

#8 
New Member
Arne
Join Date: Dec 2018
Posts: 19
Rep Power: 4 
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 

March 10, 2020, 04:22 

#9 
New Member
Join Date: Jan 2015
Location: beyond The Wall
Posts: 2
Rep Power: 0 
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. 

March 10, 2020, 07:29 

#10  
New Member
Arne
Join Date: Dec 2018
Posts: 19
Rep Power: 4 
Quote:
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) 

March 10, 2020, 09:29 

#11 
New Member
Join Date: Jan 2015
Location: beyond The Wall
Posts: 2
Rep Power: 0 
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.
Code:
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; } 

April 28, 2020, 10:46 

#12 
Member
Join Date: Apr 2018
Location: UK
Posts: 65
Rep Power: 5 
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 

April 28, 2020, 12:01 

#13 
Member
Arnout
Join Date: Nov 2010
Posts: 46
Rep Power: 12 
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!) 

April 28, 2020, 12:15 

#14 
Member
Join Date: Apr 2018
Location: UK
Posts: 65
Rep Power: 5 
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): Code:
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 What about the pMin value? 

May 4, 2020, 10:56 

#15 
Member
Join Date: Apr 2018
Location: UK
Posts: 65
Rep Power: 5 
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. Code:
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.1708508826e10 Selecting finite volume options model type limitTemperature Source: limitT_WATER  selecting all cells  selected 1920000 cell(s) with volume 8.1708508826e10 Courant Number mean: 0 max: 0 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 

May 4, 2020, 11:28 

#16 
Member
Join Date: Dec 2018
Location: Darmstadt, Germany
Posts: 87
Rep Power: 4 
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. 

May 4, 2020, 12:36 

#17 
Member
Join Date: Apr 2018
Location: UK
Posts: 65
Rep Power: 5 
Hi Raumpolizei,
First of all thanks for the reply. Here is my checkMesh: Code:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // 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 zipup check OK. Faceface connectivity OK. <<Writing 4 cells with two nonboundary 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 (nonclosed singly connected) (0 0.000218096936827 0.000375) (0.00499524110791 0.000218096936827 0.000375) bottom 3200 6401 ok (nonclosed singly connected) (0 0.000218096936827 0.000375) (0.00499524110791 0.000218096936827 0.000375) right 600 1202 ok (nonclosed singly connected) (0.00499524110791 0.000218096936827 0.000375) (0.00499524110791 0.000218096936827 0.000375) left 0 0 ok (empty) front 1920000 1923801 ok (nonclosed singly connected) (0 0 0.000375) (0.00499524110791 0.000218096936827 0.000375) back 1920000 1923801 ok (nonclosed 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 (nonempty/wedge) directions (1 0 1) Mesh has 3 solution (nonempty) 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 nonempty directions. Boundary openness (1.49847870122e13 1.28518123348e13 9.9480090689e17) OK. Max cell openness = 2.48355969881e16 OK. Max aspect ratio = 31.0351835932 OK. Minimum face area = 6.85510132929e15. Maximum face area = 6.85510130749e10. Face area magnitudes OK. Min volume = 5.35045061858e21. Max volume = 1.06992292132e15. Total volume = 8.1708508826e10. Cell volumes OK. Mesh nonorthogonality Max: 9.69700457025e06 average: 0 Nonorthogonality 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.0290308e08 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 Code:
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 cellcentred values to face values { default linear; } snGradSchemes { default corrected; } // ************************************************************************* // I think the divSchemes are causing problem, specifically div(phi,p) and div(rhoPhi,p) as the linear scheme is essentially centraldifferencing scheme and is not bounded. Perhaps I will try upwind for p and see whether there is any improvement. Any thoughts? 

May 6, 2020, 13:45 

#18 
Member
Join Date: Apr 2018
Location: UK
Posts: 65
Rep Power: 5 
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 timestep. 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 ? 

May 25, 2020, 11:33 

#19 
New Member
Liz
Join Date: Apr 2020
Location: Germany
Posts: 10
Rep Power: 3 
Hi JM27,
I'm having the same problem than you Did you find some solution to it? Cheers! 

May 25, 2020, 11:53 

#20 
Member
Join Date: Apr 2018
Location: UK
Posts: 65
Rep Power: 5 
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? 

Tags 
compressibleinterfoam, foam fatal error, fvoptions, negative temperature, twophasemixturethermo 
Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Can I use fvOptions to couple a solid region and a fluid region?  titanchao  OpenFOAM Running, Solving & CFD  3  July 1, 2021 13:22 
fvOptions: temperatureLimitsConstraint or limitTemperature not working on V3.0+  derekm  OpenFOAM Running, Solving & CFD  6  February 1, 2021 01:16 
topoSet/setSet and fvOptions  pod  OpenFOAM Running, Solving & CFD  5  April 30, 2019 05:41 
twoPhaseEulerFoam + fvOptions limitTemperature  hanness  OpenFOAM Running, Solving & CFD  5  July 19, 2018 08:53 
New output variable for source term in fvoptions  without changing the solver  vincent.clary  OpenFOAM Programming & Development  2  June 26, 2018 05:21 