CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM

Issues with power syntax

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   November 4, 2011, 12:04
Default Issues with power syntax
  #1
New Member
 
Dominic
Join Date: Jan 2011
Location: Leeds, UK
Posts: 25
Rep Power: 6
351Cleveland is on a distinguished road
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.
351Cleveland is offline   Reply With Quote

Old   November 4, 2011, 12:55
Default
  #2
Senior Member
 
Adhiraj
Join Date: Sep 2010
Location: Pennsylvania, United States
Posts: 101
Rep Power: 6
adhiraj is on a distinguished road
I would write 2/15 as 2.0/15.0 and similarly for the other constants.
Just my 2 cents.
adhiraj is offline   Reply With Quote

Old   November 9, 2011, 06:23
Default
  #3
New Member
 
Dominic
Join Date: Jan 2011
Location: Leeds, UK
Posts: 25
Rep Power: 6
351Cleveland is on a distinguished road
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.
351Cleveland is offline   Reply With Quote

Old   November 9, 2011, 18:11
Default just a suggestion
  #4
Senior Member
 
Wouter van der Meer
Join Date: May 2009
Location: Elahuizen, Netherlands
Posts: 132
Rep Power: 8
wouter is on a distinguished road
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
wouter is offline   Reply With Quote

Old   November 10, 2011, 04:13
Default
  #5
Senior Member
 
akidess's Avatar
 
Anton Kidess
Join Date: May 2009
Location: Delft, Netherlands
Posts: 919
Rep Power: 17
akidess will become famous soon enough
Quote:
Originally Posted by 351Cleveland View Post
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_;
akidess is offline   Reply With Quote

Old   November 10, 2011, 07:43
Default
  #6
New Member
 
Dominic
Join Date: Jan 2011
Location: Leeds, UK
Posts: 25
Rep Power: 6
351Cleveland is on a distinguished road
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.
351Cleveland is offline   Reply With Quote

Old   November 10, 2011, 07:56
Default
  #7
Senior Member
 
akidess's Avatar
 
Anton Kidess
Join Date: May 2009
Location: Delft, Netherlands
Posts: 919
Rep Power: 17
akidess will become famous soon enough
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?
akidess is offline   Reply With Quote

Old   November 10, 2011, 08:10
Default
  #8
New Member
 
Dominic
Join Date: Jan 2011
Location: Leeds, UK
Posts: 25
Rep Power: 6
351Cleveland is on a distinguished road
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.
Attached Images
File Type: jpg Screenshot-3.jpg (48.6 KB, 7 views)
File Type: jpg Screenshot-2.jpg (47.3 KB, 6 views)
File Type: jpg Screenshot-4.jpg (47.5 KB, 6 views)
351Cleveland is offline   Reply With Quote

Old   November 10, 2011, 08:47
Default
  #9
ngj
Senior Member
 
Niels Gjoel Jacobsen
Join Date: Mar 2009
Location: Deltares, Delft, The Netherlands
Posts: 1,620
Rep Power: 25
ngj will become famous soon enoughngj will become famous soon enough
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
Code:
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
ngj is offline   Reply With Quote

Old   November 10, 2011, 09:29
Default
  #10
New Member
 
Dominic
Join Date: Jan 2011
Location: Leeds, UK
Posts: 25
Rep Power: 6
351Cleveland is on a distinguished road
Quote:
Originally Posted by ngj View Post
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
Code:
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?
351Cleveland is offline   Reply With Quote

Old   November 10, 2011, 09:39
Default
  #11
ngj
Senior Member
 
Niels Gjoel Jacobsen
Join Date: Mar 2009
Location: Deltares, Delft, The Netherlands
Posts: 1,620
Rep Power: 25
ngj will become famous soon enoughngj will become famous soon enough
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
ngj is offline   Reply With Quote

Old   November 17, 2011, 07:15
Default
  #12
New Member
 
Dominic
Join Date: Jan 2011
Location: Leeds, UK
Posts: 25
Rep Power: 6
351Cleveland is on a distinguished road
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 is offline   Reply With Quote

Old   November 18, 2011, 07:51
Default
  #13
New Member
 
Dominic
Join Date: Jan 2011
Location: Leeds, UK
Posts: 25
Rep Power: 6
351Cleveland is on a distinguished road
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.
351Cleveland is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
what is syntax error : missing ')' before ';' aleisia Fluent UDF and Scheme Programming 8 March 10, 2015 16:42
Using Perl script and power Syntax in CFX-Pre omidiut CFX 3 October 18, 2011 06:10
[Other] Turbogrid power syntax la7low ANSYS Meshing & Geometry 0 January 15, 2011 23:20
Power Syntax Command BGK CFX 4 June 20, 2010 06:25
error while compiling the USER Sub routine CFD user CFX 3 November 25, 2002 16:16


All times are GMT -4. The time now is 17:44.