an sigFpe error on Turbulence
 

Hi Ehsan,

As I've told you in the past, the turbulence fields cannot be initialized with 0 (zero), even if they're calculated in the next iteration. In this case, you have this:
The line "value uniform 0;" should start with the same initial value you've used for other BCs and internal fields, namely ".00001943998".

Best regards,
Bruno

oh,what a careless doing!
let me test it.

Hi
the new error is:

Hi Ehsan,

I saw that error yesterday and forgot to tell you about it: it's probably due to a division by zero.
Problem is that it does not indicate where exactly this occurs. And since you have 2 patches using "groovyBC", debugging this in a single step is extremely complicated.

My suggestion is to reduce the complexity by diving the problem into smaller problems:

1. Do not run the case in parallel.
2. Use only one patch with groovyBC.
3. If this is not enough, remove the complexity of the expressions.
4. The error message does give the hint that something went wrong at
   Code:

   ---

   Foam::compressible::RASModels::kOmegaSST::correct()

   ---

   This means that one of the turbulence related fields has the problem, so focus on those fields. But keep in mind that "U" and "p" could possibly also affect this!

GroovyBC should have an option for debugging:
Problem is that I did not get any information from it :(


Good luck!
Bruno

Hi dear Bruno
you told me about crashing.
I have used different initial values(internalField and on the walls)
I had set k without mag(U) in below part but it didn't complait about that and i noticed it accidently,what may be the cause?
I'm testing on realizableKE to see if the error is related to turbulence model.
in laminar it works fine as I told you.
please let me know if you find the cause.

Hi
I did 1 and 2 without success.
maybe the problem is on mut and alphat: in zeroGradients?

Hi Bruno
I found out its because of epsilon(or omega)expression I have put there.
when I change it to a constant value run begins.this sensitivity doesn't exist on k expression(it runs when k is a formula return of a constant)
could anyone give me an advice why this error occurs on epsilon(only)?
----------------------------------------------------------------------------
I also replaced k formula but it's not because of k variable in epsilon file.
this is the error:
I don't want a constant epsilon I need it be variable.whats wrong with it?(variable k works with constant epsilon)

Hi Ehsan,

I got your email and I've taken a look at the files you had sent me before.

OK, so the problem is due to a "divide" gone wrong. The expression in question is:
The only division here is the division by "l".

If we look at the formulas in questions for calculating "l", we find these:
So, for "D_H" to be zero, most likely it's due to "a" or "b".
"a" is constant and "b" is variable... therefore, the problem is "b".
Which leads us to "pos()", which is known for easily having positions with components equal to 0 (zero). Or because "max(p.y)==p.y". This is the problem you're having!


Now, my question is: why are you using a variable "b" and not a constant "b"?

Because from what I know (which might not be 100% correct), the "hydraulic diameter" is usually a fixed value, because it's the diameter of the whole entrance.
It's like the Reynolds number: it's usually the same for the whole inlet... at least if the flow is the same for the whole inlet! ;)

Best regards,
Bruno

Oh,what an accurate man you are!exactly its for that!I 'll test some minutes later.
as you see in my case,the entrance opens to outside gradually(like a wall moves and opens to an environment) so I thought I have to use a variable hydraulic diameter and since the opening is 0 at first and becomes full opened during movement of the wall D_H should be variable.so i have to use full D_H in your opinion or my thought is correct?:)

I forgot about the "dynamic opening" detail... but... then why aren't you using "min" as well?

then it becomes whole the channel height again.I wonder why pts doesn't work?!

:confused: But "pos()" only refers to the positions of the face centers on the patch. It does not refer to points on the whole domain.
Therefore, if the points on your patch are dynamically moving, then the "pos()" values should adapt accordingly as well! Therefore, the min/max values refer only to the size of the current patch.

As I wrote above, "pos()" refers to the patch only, therefore there will be at least a few points where this is true:
Therefore, when the values are equal, the subtraction of the two values equals to zero... which leads to the division by zero on the very last expression.

I thought you already knew that "pos()" referred to the position of the face centres of the patch only.

yes I knew.the patch left has height equal to .004m(4mm) and min(pos().y) is equal to 0 and max(pos().y) is equal to .004 then the difference expression becomes .004 that is the height totally.
imagine there are 3 cells on the left patch.points refer to inlet opened and vertical line indicate that the face of cell is a wall(fixed velocity and zeroGradient p and T) and it opens gradually like this:

(why it doen't keep its order?!)

I found out why its zero.in most top cell max(pos().y) is equal to pos().y
I thought it may be solved by using pts(so that it give us the height of first cell from top of patch) but this error occurred :

Hi Bruno
I foun a way.I need the cell height.how to obtain a cell height?
maybe: (.003 is the width of cell.
is it true?
--------------------------------
Bruno it works nooooooooooooooow!
thank you for helping a looooooooooooooooooot!

a divergence occurs during run:
I even have set k and epsilon as constants:
k:
epsilon:

Hi Ehsan,

As for the error you showed in that post, there was a square root of "0" or of a negative value, which is why it crashed.

There are only 2 "sqrt" in "Foam::compressible::RASModels::realizableKE::corr ect()", at "src/turbulenceModels/compressible/RAS/realizableKE/realizableKE.C":

• Here:
  Code:

  ---

  volTensorField gradU(fvc::grad(U_));
volScalarField S2(2*magSqr(dev(symm(gradU))));
volScalarField magS(sqrt(S2));

  ---

• And here:
  Code:

  ---

  sqrt((mu()/rho_)*epsilon_)

  ---

So either "gradU" has zeros in it; or the second "sqrt" has got "mu()" or "epsilon_" equal to zero in at least one entry of the respective fields.

Best regards,
Bruno

Hi Bruno
You think its related to initial conditions?
I'll test on the backward cases i sent to you before.
Could you please do some tests which of cases you want tomorrow night and let me know your opinion to reach the answer more rapidly?
Thank you.

Hi Bruno.I attached a test case.
and this error occurs when run yPlus220:

Hi Ehsan,

I haven't looked at the test case yet, but from the error you've shown, it looks like the problem is that you forgot about using "-compressible".

I'll give some more feedback, once I manage to look at the case in about 30min.

yes.it should be for that reason.
could you please have a look into the main case also?

Uhm... I have the case running in serial (sequential mode) for over 600 seconds already and it hasn't crashed yet.

At which time step is it meant to crash?

this simple case doesn't crash,although the values of turbulence variales increase constantly.which value is suitable to be used in internal field and enterance?are the current values proper?
And crashing occurs in the case with groovyBC.

Hi Ehsan,

As tired as I am, I'm trying to figure out what to suggest. This is unfortunately waaaay beyond my level of experience.

From what I can figure out, "omega" shoots up to very high values at the walls, due to the sudden pressure/mass flow intake. I think that injecting higher values of k-omega won't affect the results... the pressure difference is simply too high. Either that, or the initial pressure is already too high inside the domain, leading "omega" to increase so much!?

OK, there are several layers of complexity at work here and you'll have to isolate one at a time:

1. Pressure: it's too high at the inlet. Start the simulation with a very small difference in pressure intake. Examine how things simulate.
   • Use atmospheric pressure levels as reference, instead of using the high pressure levels right away.
2. As you gradually increase pressure, you should begin to notice two effects:
   1. The mesh resolution starts to seem insufficient (if it was coarse to begin with).
   2. k-omega or k-epsilon start to become less functional.
   The mesh should be adapted, if it seems to need to be. And perhaps looking for other turbulence models would be a good idea, if these models start to loose their effectiveness.
3. And as I said in the past in other threads: you need simple cases where you can analytically calculate what the end result should be. I think you've already found several papers/thesis that should be able to give you some examples.
4. Pressure driven simulations can be harder to simulate. Sometimes it's easier to first do some approximate tests by injecting speed, rather than pressure. This strategy can help figure out if there is something wrong with the chosen pressure boundary conditions... or at least, this usually works on incompressible cases.

And that's pretty much all I can figure out to suggest :(

Good luck!
Bruno

thanks Bruno for your time and effort
Did you do tests on simple case with constant entrance height?
Please have a look into the case with groovyBC to see wich values are appropriate.
Feel free to change turbulence and other values to make it work.
you have a more well ordered way to investigate things than me.
Thanks

Hi Bruno,

At the moment I have no zero parameters in my boundary conditions nor initial conditions. I am trying to implement the SST k-omega solver for the conjugate heat transfer problem and I get this error which I don't quite understand. I was able to get the case working for a k-epsilon; however, for some strange reason when I change the solver and add values for kappat, nut, and omega the solver crashes. I was wondering if you could please help me understand what is this message exactly saying. Where is the equation being divided by zero and what is causing it?

Cheers!

Dear Lucas,

I guess there is something wrong with your mesh.
Please run below code and post the result.
Regards,
CFDUser_

Greetings Lucas,

As shown above, the problem is in the method "F2()": https://github.com/OpenFOAM/OpenFOAM...OmegaSST.C#L72
To know better which exact division is giving the problem, you would have to build the Debug version: http://openfoamwiki.net/index.php/HowTo_debugging

But my guess is that the new fields you've added, there were some that you initialized with 0, which as stated in post #2, is wrong.

Best regards,
Bruno

PS: Please wrap code output with the "[CODE]" markers, as explained in my signature link: Posting code and output with [CODE]

Hey Bruno,

Thanks for the help! I got it working. The F2() parameter depends on the inverse of the velocity field. Initially I left it set as (0 0 0) but I made a minor change and it worked. Sorry about the code, I was not aware about the link you posted. Thanks for letting me know!

Cheers!