Error subroutine 'cfther'

[Luciano Nagata]

Hello,

I am trying to model a air cooling system in the cs compressible module, but i keep having an error at the cfther routine, with negative values for the internal energy or for density. I've already tried to change the time step parameters lots of times, the initialisation parameters and the mesh (just tried finer meshs, because I can only do tetrahedral meshs with this geometry), but I can't avoid this error.

So any sugestions? I'm quite new to code_saturne, and this is the first compressible model I try to do. I'm sending a part of the listing file and two pictures of my geometry. The flow enters at the big circle and leaves from the two little ones (I'm using a symetric model) and to give you an scale, the diameter of the big circle is 24mm and the mesh size I'm using is 2mm for now (I tried finer meshs, but it takes a lot of extra time, and the results are not good either).

Thanks in advance.

[Brian Angel]

Hello Luciano,

I've looked through your listing file and it appears that the temperature decreases to a value of 0K by iteration 2225. You could try to set min and max values for the temperature in the GUI in the "Clipping" panel under "Numerical parameters -> Equation parameters". This may help to maintain a stable solution. The values do not have to be exact but marginally below and above the minimum and maximum values that you are expecting for this cooling system.

If the temperature differences are reasonably small, you could consider using the default solver with the density coded solely as a function of temperature in the cs_user_physical_properties.f90 routine (Boussinesq approach). This might be easier to get the simulation running and obtain a solution that could give good results. We have used this approach in the past for modelling variable density flows and it has worked well.

If you wish to try this approach, please read the tutorial for flow in a heated square cavity on this web site under "Documentation -> Tutorials".

Regards,

Brian Angel.

[Luciano Nagata]

Hello Brian,

Thanks for your answer. I could "solve" the problem improving my mesh, and using a linear imposition of the inlet pressure (so the pressure goes from 1bar to 2.5 bar in 1e-2 seconds, not instantanily).

But now I'd like to try a non adiabatic case. Is it possible using the compressible module?

And actually, I'm using the compressible model because the inlet air is under a great pressure (2 to 4 bar), not really because the density changes with the temperature, so the defaut solver can still give good results if I use the perfect gas equation for the density (is it possible to specify the density in function of the pressure)?

Best regards,

Luciano Nagata

[Erwan Le Coupanec]

Hello Luciano,

It's good that you managed to run your computation with a progressive increase of the pressure.

Is 2.5 bar the pressure you impose as Ptot (I saw in your listing file that you were using the Ptot, Htot boundary condition) ? Note that this boundary condition is a "reservoir" boundary condition.
What pressure do you impose at your subsonic outlet ?

Otherwise, about the use of the compressible model or not, what range of values for the Mach number do you expect in your case ?
Do you have an idea of the expected density variation () ?

Appart from the compressible module, another choice could be to set the option idilat to 2 (by default it is 1), which allows to take into account the unsteady term in the continuity equation, but regarding your pressure variation this option doesn't seem adapted.

Regards,

Erwan Le Coupanec

[Luciano Nagata]

Hello Erwan,

I impose a total pressure value of 2.6 bar, because the velocity at the inlet is quite low (~ 30m/s), and the outlet pressure is 1atm (1.01325 bar)

About he mach number expected, there is some regions that will have mach ~0.5-0.6 in this model, but for other models I'm going to do I expect some regions with mach 0.9, and for the density variation at the inlet is about 2.8 kg/m3 et at the outlet it should be around 1.6 kg/m3 or even less (minimum 1.2kg/m3).

And about my results, I have a question, even if the outlet imposed pressure is 1atm, the results give me a minimum pressure of 1.75bar isn't it weird? I thought it was a convergence problem, but even if I use more iterations this doesn't change.

Thank you for your help

Best regards

Luciano Nagata

[Luciano Nagata]

And about my results, I have a question, even if the outlet imposed pressure is 1atm, the results give me a minimum pressure of 1.75bar isn't it weird? I thought it was a convergence problem, but even if I use more iterations this doesn't change.

Actually,I found a solution for this problem, I don't exactly why it works, but when I had the problem I was using the gradient calculation method "least square method over extended cell neightborhood", which in a pipe model i've donne before gave me a faster convergence, so I tried changing it to the defaut option (Iterative Handling of non-orthogonalities) and the problem disappeared. Now the results shows a pressure near the outlet that is the same I imposed at the bounduary conditions.

[Erwan Le Coupanec]

Hello Luciano,

The gradient calculation method "iterative handling of non-orthogonalities" is more precised and ensures conservativity, while the "least square method over extended cell neighborhood" is more robust and should be used for example when the quality of the mesh is not good enough.

If the computation runs fine with the default method, i'd say it's better to stick to it.

About one of your previous question, what kind of non adiabatic boundary condition do you want to impose ? With the compressible module, you can't use the gui to impose a heating wall, but you can do it in the user subroutines.

Regards,

Erwan Le Coupanec