Foam::pow function error
Hello again forum,
I've been having trouble with a power function I've implemented in a modified compressible k-epsilon turbulence model. The equation is as follows: rEQ = Cr * Foam::pow((den_*yliq)/denliq_,0.133)*((st_/Foam::pow(epslion_,0.4)*63.1) There are no errors during the compiling process and this identical equation has worked in a code based on the incompressible k-epsilon RAS model. Cr is a constant and den, yliq , denliq and st are all volScalarFields; when attempting to solve, the following error is thrown up: Starting time loop Time = 1e-05 Code:
smoothSolver: Solving for Ux, Initial residual = 1, Final residual = 0.0637079, No Iterations 4 as when removing the yliq term, the equation solves normally with no errors, however yliq does not have any values at zero or negative as a filter function has been included to prevent this? Can anybody shed some light on the problem? Thanks. Dom. |
And what about denliq, as you are dividing by that?
Also, I would expect your solver to complain about epslion, instead of epsilon. |
denliq is set at a constant dimensioned value of 1000 - but is implemented as a field rather than a constant due to a number of other tests I was running, unfortunately substituting it with an integer or removing it altogether still results in the same problem.
I attribute epslion to my inconsistent typing :) It is epsilon in the equation itself. |
Can you show you complete correct function? Especially your yvap equation, as it is not solving anything, according to your logfile.
|
Certainly, I've included the yliq equation - yvap currently isn't fully implemented so won't be solving anything at this stage.
Code:
rEQ_ = Cr_*(Foam::pow((den_*yliq_)/denliq_,0.1333))*(st_/(Foam::pow(epsilon_,0.4)*63.1)); |
How do you initialize epsilon?
|
Fixed value (very turbulent flow) at twin air inlets and the compressible epsilon wall function at the boundaries with a value of 1e-12 and zeroGradient at the outlet - its my first compressible case so I'm adapting as best I can from the incompressible version.
|
I mean the internalField. If that is zero, it would explain the division by zero you are probably encountering.
|
What is the sign of
Code:
(den_*yliq_)/denliq As you receive error from libm I guess you're out of the function domain. |
Ah, my apologies for missing that out. That is also set at a fixed value of 0.3 - identical to one of the air inlet values. I've managed to get a modified form of the equation working as follows:
Code:
|
Quote:
|
It is really difficult to tell what is going wrong, as you provide incomplete information. Clearly from the piece of code you are sharing, the problem arises after the solving (and construction) of yliqEqn. Please, share your whole code, or just try to old school debug with Info statements.
|
After taking your advice, I built a debug version of openfoam 2.2.1 and ran my code using gdb as the debugger and found something very interesting from the error log:
Code:
Since the s term term comes out negative despite the filter function attached to yliq - can I assume that it is not working correctly and producing -ve values for yliq? Or does the problem lie elsewhere in the term? Thanks for the assistance. |
I'm not quite sure you've shown right piece of code. According to solver output yliq equation is solved. Can you show the code after solution of yvap equation?
|
I've included the yliq equation and the correction routine along with the equations implemented immediately after in the .C file. Although the yliq field seems to solves without issue, it is only when being used to calculate rEQ_ that the problem seems to occur.
I'm beginning to wonder if the correction routine I put in (based on the number of cells in the mesh) only takes effect after everything is solved - as the output file for yliq contains no negative values at all.....but seems strange that it worked in exactly the same form for the incompressible version of the same model I made. Code:
P { margin-bottom: 0.21cm; } tmp<fvScalarMatrix> yliqEqn |
Why not use either
Code:
forAll(yliq_, cellI) or Code:
yliqMin = dimensionedScalar("yliqMin", ...dimensions..., 0.0); Code:
yliqMin = dimensionedScalar("yliqMin", ...dimensions..., 0.0); |
It appears that's solved the problem nicely - looks like the routine I used wasn't working correctly and causing negative values to slip through into the following equations.
I would like to thank you all very much for your help, hopefully this is the last time the problem will show up. Thanks again. Dominic. |
Dear Dominic, I want to emphasize that the problem was solved within 20 minutes after posting the relevant lines. Please keep in mind to carefully prepare you post when asking a question here. If you reread the thread, you will see that I already asked you to post the relevant piece of the code two days ago, which could have saved you a lot of time.
|
All times are GMT -4. The time now is 22:04. |