CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   Fluent UDF and Scheme Programming (http://www.cfd-online.com/Forums/fluent-udf/)
-   -   How to make a UDF to have a velocity profile in a square channel (http://www.cfd-online.com/Forums/fluent-udf/64463-how-make-udf-have-velocity-profile-square-channel.html)

Gigis May 12, 2009 10:22

How to make a UDF to have a velocity profile in a square channel
 
i have a 3D domain as a square channel. I had set the inlet as velocity inlet in Gambit..And now come into Fluent, i have to use a UDF programming to set a parabolic flow veloctiy inlet and difference pressure at the inlet.. But i dont know how to write a programme about this 3D domain... Anyone can help? thanks


Here is what I tried:

/************************************************** ********************
vprofile.c
UDF for specifying steady-state velocity profile boundary condition
************************************************** ********************/
#include "udf.h"
#include "iostream"


DEFINE_PROFILE(inlet_x_velocity, thread, position)
{
real x[ND_ND]; /* this will hold the position vector */
real y, z;
face_t f;
x[0]=0.07
begin_f_loop(f, thread)
for(i=1:15;i=i+2)
{
F_CENTROID(x, f, thread); /*fluent finds an error in this line, why??*/
y = x[0];
z = x[0];
F_PROFILE(f, thread, position)=(4.*0.0005*0.0005*0.1/(3.14*3.14*3.14*0.001*0.006))*sigma((1/(i*i*i))*(1-((cosh(i*3.14*y/0.0005))/(cosh(i*3.14/2))))*(sin(i*3.14*z/0.0005));
}
end_f_loop(f, thread)
}

Daniel Tanner May 12, 2009 17:03

The UDF looks ok, though I think you are not understanding the UDF or "position vector" fully.
begin_f_loop(f, thread) is used to loop over every element (face) on the inflow boundary.

F_CENTROID(x, f, thread) returns the "position vector" (x) of the centroid of each face for every iteration of the loop. The "position vector" is x[] where x[0] is the x value, x[1] is the y value and x[2] is the z value of each centroid. Remember there is one centroid for each face element of the inlet surface. NOTE: you have used x[0] = 0.07 before the loop, this is meaningless as it will be overwritten by F_PROFILE at each iteration.

F_PROFILE(f, thread, position) allows you to prescribe the value of a variable at this inlet surface, e.g., the velocity at each wall centroid. Hence, y = x[0] and z = x[0] just means that all are X position. Note: ND_ND tells the UDF if you are running a 2D or 3D solver, hence, x holds three values.

Your loop looks strange, check the syntax, i.e., "for(i=1:15;i=i+2)" might be something like "for(i = 0; i < 10; i++)". You must also define i, i.e., "int i;"

Does this answer your problem? If not submit the actual error or the nature of the error (e.g., non-sensible values, compilation problems ...).

Gigis May 12, 2009 18:43

Thanks for reply

This is my new UDF below, But I still do not understand everything into.

It does not works anymore. Fluent says:" ligne 17 invalid lvalue in assignment".

Could you help me a little bit again?

Thank you

Gigis

/************************************************** ********************
vprofile.c
UDF for specifying steady-state velocity profile boundary condition
************************************************** ********************/
#include "udf.h"
#include "iostream"


DEFINE_PROFILE(inlet_x_velocity, thread, position)
{
real x[ND_ND]; /* this will hold the position vector */
real y, z;
face_t f;

x=x[0]

(L.17) begin_f_loop(f, thread)
int i
for(i = 0; i < 10; i+=2)
{
F_CENTROID(x, f, thread);
y = x[1];
z = x[2];
F_PROFILE(f, thread, position)=(4.*0.0005*0.0005*0.1/(3.14*3.14*3.14*0.001*0.006))*sigma((1/(i*i*i))*(1-((cosh(i*3.14*y/0.0005))/(cosh(i*3.14/2))))*(sin(i*3.14*z/0.0005));
}
end_f_loop(f, thread)
}

Daniel Tanner May 13, 2009 11:55

Some quick observations:

1) the "int i;" command must come at the start with the rest of the declared variables (i.e., x, y, z ...)
2) all lines must finish with ";"
3) better to use x[0], x[1] and x[2] directly (not x, y, z)

#include "udf.h"
#include "iostream" /* Don't know why you need this and should probably be "iostream.h"*/

DEFINE_PROFILE(inlet_x_velocity, thread, position)
{
real x[ND_ND]; /* this will hold the position vector, these are real numbers */
int i; /* i is an integer and statement must finish with ";"*/
face_t f;

x=x[0]; /* you had no semi-colon and this means nothing, it will be overwritten in the next loop! So best to remove it!*/

begin_f_loop(f, thread) /* this is the loop that iterates through every element on the inlet boundary!!!! This is the same as a for(i = 0; i < Number_of_elements_on_inlet_face; i++)*/

for(i = 0; i < 10; i+=2) /* what does this do??????? This will stop the boundary condition from being applied to every part of the inlet */
{
F_CENTROID(x, f, thread); /* this statement returns the centroid location to position vector x, hence above x = x[0]; means nothing. */
/* The next line prescribes the calculated value as a boundary condition (e.g., velocity) for the current inlet element. Hence, the loop Begin_f_loop iterates through each element and prescribes the boundary condition based on its location, which is determined by F_Centroid */
F_PROFILE(f, thread, position)=(4.*0.0005*0.0005*0.1/(3.14*3.14*3.14*0.001*0.006))*sigma((1/(i*i*i))*(1-((cosh(i*3.14*x[1]/0.0005))/(cosh(i*3.14/2))))*(sin(i*3.14*x[2]/0.0005));
}
end_f_loop(f, thread)
}

Daniel Tanner May 13, 2009 12:01

See the example below! This will write a profile to the inlet face. For a velocity inlet boundary condition this is effectively setting a velocity at each inlet element as a function of the x*y*z location of each face centroid!

I hope this helps!

#include "udf.h"

DEFINE_PROFILE(inlet_x_velocity, thread, position)
{
real x[ND_ND];
face_t f;

begin_f_loop(f, thread)
{
F_CENTROID(x, f, thread);
F_PROFILE(f, thread, position)= x[0]*x[1]*x[2]; /* same as x*y*z */
}
end_f_loop(f, thread)
}

Gigis May 14, 2009 16:39

OK thanks for the answer,

So, I tried to arrange my UDF and now it seems to work! I just made another one more simple. without the sum (=sigma) et without "for(i = 0; i < 10; i+=2)" so.

But, now, the matter is there is not the velocity value I want. The maximum value is about 0.00032m/s, While in the UDF i defined x[0]=0.07m/s.

It is one thing I cannot understand, could you please explain me or give me a solution?

Thanks a lot

Regis

Here my UDF:
/************************************************** ********************
vprofile.c
UDF for specifying steady-state velocity profile boundary condition
************************************************** ********************/
#include "udf.h"
#include "iostream.h"


DEFINE_PROFILE(inlet_x_velocity, thread, position)
{
real x[ND_ND]; /* this will hold the position vector */
face_t f;
x[0]=0.07;

begin_f_loop(f, thread)

{
F_CENTROID(x, f, thread);
F_PROFILE(f, thread, position)=(4.*0.0005*0.0005*0.1/(3.14*3.14*3.14*0.001*0.006))*(((1/(1*1*1))*(1-((cosh(1*3.14*x[1]/0.0005))/(cosh(1*3.14/2))))*(sin(1*3.14*x[2]/0.0005)))+((1/(3*3*3))*(1-((cosh(3*3.14*x[1]/0.0005))/(cosh(3*3.14/2))))*(sin(3*3.14*x[2]/0.0005)))+((1/(5*5*5))*(1-((cosh(5*3.14*x[1]/0.0005))/(cosh(5*3.14/2))))*(sin(5*3.14*x[2]/0.0005)))+((1/(7*7*7))*(1-((cosh(7*3.14*x[1]/0.0005))/(cosh(7*3.14/2))))*(sin(7*3.14*x[2]/0.0005))));
}
end_f_loop(f, thread)
}

Daniel Tanner May 14, 2009 17:40

Sorry if I was unclear!

Firstly, you set x[0] = 0.07 but it does not appear in your equation for the profile! If you wish to insert it then you can just insert 0.07 directly, i.e., (10*x[1])+1 is the same as (10*0.07)+1.

Secondly, if you wish to use x[0] = 0.07 then you must do so after the F_CENTROID command as this sets x[0] to the x value of the centroid of the current face! i.e.:
begin_f_loop(f, thread)
{
F_CENTROID(x, f, thread);
x[0]=0.07;
F_PROFILE(f, thread, position)=.....

Thirdly, I am not an expert in C but I would suggest that you do not use integers in your calculations. In some programming languages this can lead to the processor performing an integer calculation rather that single or double precision, i.e. replace (1/(1*1*1)) with (1.0/(1.0*1.0*1.0)). Not sure if it makes a difference in Fluent but no harm to include the first decimal place.

Finally, if the profile is correct but the peak value is not as desired then just multiply everything to scale to the correct peak value, i.e. if you want a peak of 1 m/s then multiply the profile by 1/0.00032:
F_PROFILE(f, thread, position)=(1.0/0.00032)*.....

n7310889 September 3, 2012 07:59

Hi Dani,
I'm writing an udf to assign heat generation profile in the outer wall of (symmetric right hand side half section) of a 3D annular fluid domain of 33mm radius. This fluid is flowing through a 35 mm radius and 8 m long horizontal tube. The heat generation (q) profile varies along the circumference (horizontal x & vertical y direction) of the tube, and assumed constant & steady state along the length (horizontal z direction). Coordinates of the axis of rotation of the tube are (0,0,0) and (0,0,8.0).
My earnest request to you, please check the udf written for a Fluent 3D model. I'll remain grateful to you. Also please check the logic of assigning radius value in the macro. Since my model contains only fluid domain of 33 mm, then what should be assigned value for radius among 33 mm of fluid body or 35 mm including tube? Please suggest me.

#include "udf.h"

#define RADIUS 0.033
#define CONVEC_LOS_COFICENT 0.05
#define TUBE_EMISIVITY 0.1378
#define GLAS_EMISIVITY 0.05
#define SIGMA 5.67e-08
#define DIRECT_NORMAL_INSOLATION 937.9
#define T_KELVIN 273.0


DEFINE_PROFILE(wallheatfluxprofile,thread,index)
{
real x[ND_ND];
real glas_tube_emsivity,T_amb,T_sky,T_glass,T_gas,rad_l oss,convec_loss;
double a,z,q;


face_t f;

begin_f_loop(f,thread)
{
F_CENTROID(x,f,thread);

glas_tube_emsivity=(TUBE_EMISIVITY*GLAS_EMISIVITY)/((TUBE_EMISIVITY+GLAS_EMISIVITY-(TUBE_EMISIVITY*GLAS_EMISIVITY));

T_amb =T_KELVIN+28.5;
T_sky =T_KELVIN+14.0;
T_glass =T_KELVIN+121.0;
T_gas =T_KELVIN+30.0;

rad_loss =TUBE_EMISIVITY*SIGMA*((2*pow(WALL_TEMP_INNER (f,thread),4)-pow(T_amb,4)-pow(T_sky,4))+glas_tube_emsivity*SIGMA*(pow(WALL_T EMP_INNER (f,thread),4)-pow(T_glass,4));
convec_loss =CONVEC_LOS_COFICENT*(WALL_TEMP_INNER (f,thread)- T_gas);

z=x[2];
for(z=0.0; z<=1.5; z++)
{
if((z>=0.0) && (z<=1.2))
{
for(a=0.0; a<=180.0; a++)

{
x[0] =RADIUS*sin(a);
x[1] =-RADIUS*cos(a);


if ((a>=0.0) && (a<=15.9))
{
q =4.1643*pow(a,3)-60.818*pow(a,2)+485.19*a+32320;
}
if ((a>15.9) && (a<=47.5))
{
q =0.1288*pow(a,3)-6.7478*pow(a,2)+293.03*a+37992;
}
if ((a>47.5) && (a<=86.5))
{
q =1.2582*pow(a,3)-256.73*pow(a,2)+15690*a-250299;
}
if ((a>86.5) && (a<=180.0))
{
q =-0.0004*pow(a,3)-0.0716*pow(a,2)+47.947*a-3026.1;
}


F_PROFILE(f,thread,index) =(q*DIRECT_NORMAL_INSOLATION/937.9)-rad_loss-convec_loss;
}
else
{
F_PROFILE(f,thread,index) =0.0;
}
}
end_f_loop(f,thread)
}

mactech001 January 13, 2013 23:20

output value in DOS window
 
Quote:

Originally Posted by Daniel Tanner (Post 215970)
Some quick observations:

#include "iostream" /* Don't know why you need this and should probably be "iostream.h"*/


}

Dear Daniel,

hello.

i would like to output the value i apply to F_PROFILE to the DOS window by using 'cout'. the following is my code:

#include "udf.h"
#include "iostream.h"
DEFINE_PROFILE(wallheatgenerate,thread,i)
{
real source=0.001;
face_t f;
begin_f_loop(f,thread)
F_PROFILE(f,thread,i)=source+0.001;
cout<<"F_PROFILE output:"<<F_PROFILE(f,thread,i)<<endl;
end_f_loop(f,thread)
}

but, in Fluent, the compiling claims that it can't open the include file "iostream.h".

can you help me to see what's wrong with the code please?

my compiler is VS C++ 2008 express, and it compiles the 'Hello world' code fine with no problems.


All times are GMT -4. The time now is 15:26.