[351Cleveland]

Hello forum, have hit a bit of a snag and was wondering if I could get some pointers.

In my solver I have implemented the following equation:

rEQ_ = Cr_*pow(((den_*yliq_)/denliq_,2/15)*(pow(st_,3/5)/(pow(epsilon_,2/5)*pow(denliq_,2/15))))+smv_;

Although this equation does solve when my case is run, the internal field of the results domain only assumes a single value for the whole domain other than the specified boundary conditions. All variables within the equation are scalars with st and denliq as constants.

This equation should only have a solution within the area where yliq > 0 and smv is the small value it should take everywhere else, instead this single value takes up the whole domain.

Am I missing something obvious or is there some issue with the power syntax which I should be aware of?

Cheers.

[adhiraj]

I would write 2/15 as 2.0/15.0 and similarly for the other constants.
Just my 2 cents.

[351Cleveland]

Hey adhiraj

thanks for the idea, unfortunately didn't seem to make a difference. Its worth noting that if I remove the power functions, the equation solves normally with a unique value for each cell but the wrong result. Must confess I'm a little confused.

[wouter]

hello,
If I have to debug such a complicated equation, I try to use a variable for every term/factor and print these so to see which term/factor is the bad one.

just a suggestion

Wouter

[akidess]

Originally Posted by 351Cleveland

rEQ_ = Cr_*pow(((den_*yliq_)/denliq_,2/15)*(pow(st_,3/5)/(pow(epsilon_,2/5)*pow(denliq_,2/15))))+smv_;

Although this equation does solve when my case is run, the internal field of the results domain only assumes a single value for the whole domain other than the specified boundary conditions. All variables within the equation are scalars with st and denliq as constants.

This equation should only have a solution within the area where yliq > 0 and smv is the small value it should take everywhere else, instead this single value takes up the whole domain.

You say all of your variables are scalars, so it is not surprising that your field is constant everywhere. One of your quantities will have to be a volScalarField! E.g. make yliq a volScalarField like rEQ, and then write the equation as:

rEQ_ = Cr_*pow(((den_*yliq_)/denliq_,2/15)*(pow(st_,3/5)/(pow(epsilon_,2/5)*pow(denliq_,2/15)))) * pos(yliq_) + pos(-yliq_) * smv_;

[351Cleveland]

Hi akidess,

I apologise for not being clear enough in my initial question, epsilon is a volScalarField, present in the incompressible RANS k-epsilon model, as is yliq.
The paper which references this particular equation can be found at

http://www.tandfonline.com/doi/abs/1...02200801893796

From your equation:

rEQ_ = Cr_*pow(((den_*yliq_)/denliq_,2/15)*(pow(st_,3/5)/(pow(epsilon_,2/5)*pow(denliq_,2/15)))) * pos(yliq_) + pos(-yliq_) * smv_;

I am wondering about the purpose of the pos(yliq) terms as yliq is bound between 0 and 1, so would it need multiplying by the positive field again?
Also, I thought the addition of the smv_ term on its own would prevent any divide by zero errors? Got a compilation error to do with the multiply by *pos(yliq_) too unfortunately.
Must admit, this is proving a tricky one.

Thanks for the help.

[akidess]

You wrote "This equation should only have a solution within the area where yliq > 0 and smv is the small value it should take everywhere else", so I assumed yliq can be smaller zero. So I guess you just want smv only when yliq = 0:

rEQ_ = Cr_*pow(((den_*yliq_)/denliq_,2/15)*(pow(st_,3/5)/(pow(epsilon_,2/5)*pow(denliq_,2/15)))) * neg(-yliq_) + pos(-yliq_) * smv_;

If yliq = 0, pos(-yliq) should return 1, and neg(-yliq) should return 0. So "pos(-yliq_) * smv_" should return a volScalarField that is smv_ where yliq=0, and 0 everywhere else. Does that make any sense?

[351Cleveland]

Yes, that makes sense, my C++ isn't quite up to scratch yet

I've attached a couple of screenshots to illustrate the problem, the first being at time zero which just shows the boundary conditions, and the next a few steps into the equation (each step = 1e-4s).
The second is the value which rEQ assumes immediately after the start for all time steps, the final is the area which yliq has a positive value at the same time as rEQ was taken.

There may something to do with the constant Cr at he beginning as the assumed value of the red area is identical to that of the constant over the whole domain.

[ngj]

Hi Dominic

I have had problems with the compiler being unable to distinguish between the different types of pow, e.g. std:ow and Foam:ow, so to be clear, I tend always to write

Foam::

in front of such functions as pow, sin, cos, sqrt, etc. If that is the case, then pow(<something>, 2/15) = pow(<something>,0) = 1;

Kind regards,

Niels

[351Cleveland]

Originally Posted by ngj

Hi Dominic

I have had problems with the compiler being unable to distinguish between the different types of pow, e.g. std:ow and Foam:ow, so to be clear, I tend always to write

Foam::

in front of such functions as pow, sin, cos, sqrt, etc. If that is the case, then pow(<something>, 2/15) = pow(<something>,0) = 1;

Kind regards,

Niels

Hi Niels,
might be getting somewhere with that, although compiling correctly it now throws up the following error message when trying to run the case:

--> FOAM FATAL ERROR:
LHS and RHS of + have different dimensions
dimensions : [0 -0.8 1.2 0 0 0 0] + [0 0 0 0 0 0 0]


From function operator+(const dimensionSet& ds1, const dimensionSet& ds2)
in file dimensionSet/dimensionSet.C at line 402.

FOAM aborting

From the dimensions I'm using:

st = dyn/cm = g s^-2
epsilon = cm^2 s^-3
denliq = g cm^-3

den and denliq cancel each other out and yliq is dimensionless

therefore the dimensions should work out as: cm^-1 s^-1
yet, this puts the dimensions as: -0.8 and 1.2.
Does the power affect the dimensions or is there another reason behind this error?

[ngj]

Hi

It suggests that smv_ is dimensionless, which makes sense as long as it is merely a scalar. Making it work you need to define smv_ as a dimensionedScalar with the correct dimensions.

Yes, Foam:ow( dimensionedScalar, exp) does change the dimensions of the dimensionedScalar.

/ Niels

[351Cleveland]

back again,

have been making some progress in the equation - the help was much appreciated. But now I'm stuck with regards to setting the correct units.

In trying to set a dimensionSet for smv in my transport dictionary, I can't find the way to allow the turbulence model to read from the transportProperties file inside the case/constant folder.

Is there a way in which i can set the dimensionSet inside the turbulence model when defining the constant or do I have to read from the external file, and if so, how would I go about this?

Cheers.
Dom.

[351Cleveland]

Problems solved and working nicely now, thanks for the assistance

Turns out it was a combination of issues, the main being the issue with the power syntax, the rEQ field now solves perfectly, and the other problems were compounded through the changing of the units.
Turning off the dimension checking actually helped to work through these issues quicker, details on how to do this can be found on another thread in this forum.

Thanks again.
Dom.