
[Sponsors] 
fvOptions limitTemperature crashing in compressibleInterFoam 

LinkBack  Thread Tools  Search this Thread  Display Modes 
October 20, 2020, 04:18 

#21 
Senior Member
Pawan Ghildiyal
Join Date: Nov 2015
Posts: 135
Rep Power: 8 
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; } 

October 20, 2020, 04:22 

#22  
Senior Member
Pawan Ghildiyal
Join Date: Nov 2015
Posts: 135
Rep Power: 8 
Quote:
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; } 

February 22, 2021, 15:28 
Has any one figured this out yet?

#23 
New Member
Adil A
Join Date: Oct 2015
Posts: 3
Rep Power: 8 
I have noticed that this problem doesn't occur for slower flows. Not sure what that means.


December 30, 2021, 03:06 

#24 
New Member
wenjing
Join Date: Apr 2017
Posts: 1
Rep Power: 0 
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 

January 13, 2022, 07:01 

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

February 22, 2022, 16:10 

#26 
Senior Member
Join Date: Sep 2013
Posts: 331
Rep Power: 18 
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 Code:
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()() ) 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 Last edited by Bloerb; March 21, 2022 at 15:21. 

March 22, 2022, 03:34 

#27  
New Member
Join Date: Apr 2021
Posts: 8
Rep Power: 3 
Quote:
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 farfield. However, while the bubble is near the solid wall, the error of negative initial temperature T0 always occurs during the generation of highspeed 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. 

March 22, 2022, 05:00 

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

March 22, 2022, 07:56 

#29 
New Member
Join Date: Apr 2021
Posts: 8
Rep Power: 3 
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 farfield 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 nonreflect 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. 

March 22, 2022, 08:51 

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

March 22, 2022, 10:54 

#31 
New Member
Join Date: Apr 2021
Posts: 8
Rep Power: 3 
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 farfield 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 farfield, 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 farfield 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 farfield. However, it seems like not work in the nearwall, 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 OpenFOAMextended4.1, the compressibleInterFoam do not couple with TEqn, and this solver can fix this problem. 

March 22, 2022, 13:33 

#32 
Senior Member
Join Date: Sep 2013
Posts: 331
Rep Power: 18 
Example:
Code:
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; } ..... } Code:
TEqn.relax(); fvOptions.constrain(TEqn); TEqn.solve(); fvOptions.correct(T); mixture.correctThermo(); < it crashes here for me mixture.correct(); 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) Code:
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_); } } } 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) Code:
laplacianSchemes { default Gauss linear orthogonal; } snGradSchemes { default orthogonal; } 

March 23, 2022, 07:31 

#33  
New Member
Join Date: Apr 2021
Posts: 8
Rep Power: 3 
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 farfield 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: Quote:


March 23, 2022, 12:57 

#34 
Senior Member
Join Date: Sep 2013
Posts: 331
Rep Power: 18 
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.


March 23, 2022, 23:34 

#35 
New Member
Join Date: Apr 2021
Posts: 8
Rep Power: 3 
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 OpenFOAMv9, 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. 

March 30, 2022, 11:21 

#36 
New Member
Sumit
Join Date: Jul 2017
Posts: 5
Rep Power: 7 

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  4  January 14, 2022 07:55 
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 