udf for heat generation rate
 

hi every body
I am supposed to write a udf for heat generation rate which is a cosine function
#include"udf.h"
DEFINE_SOURCE(energy_source,c,t,ds,eqn)

{
real x;
ds[eqn]=0.3*cos(0.3*x);
return sin(0.3*x) ;
}
but there are some errors! any one can help me to modidfy it. thanks

Sarah,

Is "x" in your function the x-coordinate of the cell? If so, try the following un-checked code.

The source in this case is actually the sine function, and the derivative (dS) is the cosine function. If your heat generation rate is actually the cosine function as you've listed it, the code is:

ComputerGuy

ComputerGuy
I am really thankful of you for your response , but X is the height of my geometry (of course its better to use y instead of x , I did mistake to type x!) my geometry its here:

Heat generation rate is a cosine function of y (height) :
source =5.958*sin(pi *y/11)+1.866
So
ds= 5.958*pi /11*cos(pi* y/11)
I am grateful if you help me to write a udf for it .
Best regards

Sarah,

I have given you the code you need to make your function work. Instead of x[0], which returns the "x" component of the geometry, you'll need to use x[1], which returns the "y" component.

The only other thing is changing the source=.... and dS[eqn]=.... to the heat generation function and its derivative.

ComputerGuy

Dear ComputerGuy
I am new in udf so there is some thing that I don't know about them , such as ND could you explain me about ND , why you define x[ND_ND] ?
thanks
Sarah

Sarah,

You should read through the Fluent UDF manual. However, from the Fluent 6.3 UDF manual (Fluent, Inc., September 11, 2006), page 3-26:

"The constant ND_ND is dened as 2 for RP 2D (FLUENT 2D) and RP 3D (FLUENT 3D). It can be used when you want to build a 2 x2 matrix in 2D and a 3 x3 matrix in 3D. When you use ND_ND, your UDF will work for both 2D and 3D cases, without requiring any modications."

Thus, when I defined "real x[ND_ND]," I was making x a position vector with dimension ND, where n is the number of dimensions in your simulation.

x[0] would be the x component
x[1] would be the y component
x[2] would be the z component

The manual has extremely thorough examples -- go through it.

ComputerGuy

Hello Computer Guy
 

My problem is also related to heat generation in the energy eqn..

In my model Qg(heat generation) is proportional to exponential function of temperature and time gradient of new variable (dC/dt)
Initial condition of C is specified i.e t=0,C=0
My aim is to calculate a C at different times..
Is this possible in fluent to include the new variable in the fluent UDF...
I have written the UDF for the exponential function of temperature but how to include dc/dt is still a question for me..

Dear ComputerGuy
Yes , you are right , I have read this manual and its so good but there is not any theme about porosity , really in my project I have to define porosity as a function of x and I don’t know how to write a udf for it , everybody who I asked them , say: its impossible or very hard by fluent 6.3.26!
Could you help me about this issue ?of course I am really thankful of you for your great help about heat generation rate .
Best regards
Sarah

Manohar,

Is C a scalar? Is it a property of the fluid which travels with the fluid (like species mass fraction, for example)? Or, is it something about the cell which simply changes with time?

For example, if you wanted to track the change in temperature in a cell such that you approximate dT/dt (or dC/dt in your case), there are macros which allow you to access the previous time step of a given cell. The previous temperature in a cell may be accessed by:
Many of the variables of the previous time step can be accessed in this way. Have a look at the UDF manual, section 3.2.3 (Fluent 6.3) for more information. With this, it should be very simply to approximate a derivative.

If you can't find your answer there, or if Fluent doesn't automatically track what you're looking for, you could always create a few User-Defined Memory Locations to track values from a previous time step.

That is, assume C_UDMI(cell,thread,2) is the current value of C and C_UDMI(cell,thread,1) is the previous value:



If you're looking for a "value" or "property" of the fluid which will move with a fluid parcel (advect, diffuse, etc.), look into User-Defined Scalars.

If you provide a little more information, especially the functional form of your equation and what you're trying to track (i.e., what is C), perhaps you'll be able to get more help.

ComputerGuy

Sarah,

If you can write an equation (or an algorithm) for porosity as a function of geometric position, you can write a UDF for it. It's not that hard. The trick is that as opposed to DEFINE_SOURCE for energy, you must use a DEFINE_PROFILE UDF. Set the region you're interested in to a porous zone, code the DEFINE_PROFILE UDF, compile or interpret the UDF, then load the UDF into the porosity drop-down menu.

Let me know if you need further assistance!

ComputerGuy

Hello Computer Guy
 

My model is related due to absorption(Where gas is getting absorbed by solid).
I am solving energy equation(Unsteady) for the cylindrical model
In my model
Qg = exp(-const/T)*dc/dt ----(1)
where dc/dt = (pg - Peq)/peq*(some constant) ----(2)
where pg =constant= 3 bar
peq is a function of Temperature(T).
C is concentration of gas in the solid
My aim is to find to find T and C for different time.
I can substitute (2) eqn in (1) eqn and total heat source will be in function of temperature but finding C by means of temperature will be a long procedure.
Is there any way to find both C and T parallel by writing UDF??

Dear ComputerGuy
Thank you It is very kind of you . well I have written an equation for porosity , Its here :
Porosity= 0.39*(1+1.05*exp(-100*(y-1)/0.06)) 1 < y<1.425
Porosity= 0.39*(1+1.05*exp(-100*(1.85-y)/0.06)) 1.425<y<1.85
y=1 , y=1.85 porosity =0.39
And I have written a udf for it :


#include "udf.h"

DEFINE_PROFILE(porosity_profile,t,i)
{
real x[ND_ND];
real y;
cell_t c;

begin_c_loop(c,t)
{

y=x[1];
if(1<y<1.425)
F_PROFILE(c,t,i)=0.39*(1+1.05*exp(-100*(y-1)/0.06));
else if(1.425<y<1.85)
F_PROFILE(c,t,i)=0.39*(1+1.05*exp(-100*(1.85-y)/0.06));
}
end_c_loop(c,t)
}

is it corresponding to the equation of porosity? And when we use udf for porosity so we have to write a udf for viscous resistance and inertial resistance as well , so should we interpereted and hook udfs respectively ?, It means we should first hook porosity_udf then viscous resistance_udf and inertial resistance_udf , isn’t it ?
I should also thank you since your advices help me in my project so far .
Best regards
Sarah

Manohar,

If I understand correctly:

1. Heat is being generated by the absorbtion of a gas into a solid.
2. The heat generation rate depends on the temperature, as well as the difference in pressure of the gas from its equilibrium value (vapor pressure).
3. The vapor pressure is a function of temperature, and is therefore time dependent.
4. You want an expression for, or an output of, both T and C

My suggestions would be the following:

1. Have Fluent monitor and output the temperature as a function of flow time (not iteration); this can be volume averaged, mass averaged, etc. You need to decide where you want to monitor the temperature, because there appears to be a transient, spatially dependent thermal behavior in your system. Thus, your choice of location may change your expression for T as a function of time.
2. Have Fluent monitor and output Peq. Better yet, set up a custom field fuction which is From there, you could very easily numerically integrate the result (which is dC/dt) and come up with an answer for C. If that doesn't suit you, regress the data (i.e., fit a curve through it), develop an expression which approximates dC/dt, and analytically integrate that.

Many times problems of this sort don't have a nice analytical solution. If they did, we probably wouldn't be using Fluent to solve them.

ComputerGuy

Sarah,

Your UDF was basically correct, although didn't address the case where y=1 or y=1.85. You also need to assign a vector to x, which I did with the C_CENTROID command. I've changed it below:

See also: http://my.fit.edu/itresources/manual...udf/node42.htm

Examples 5 and 6 show the usage of both F_PROFILE and C_PROFILE for defining porosity. You can do the same for viscous and inertial resistances, using basically the same code as above, but modified with the appropriate values.

ComputerGuy

Dear ComputerGuy
Thank you for modification of my code , but now I have faced a funny problem! While I interpreting a udf , previous udf would be removed , for example , in my case ,first I interpreted and hooked udf of energy _source and then udf of porosity but I saw energy_source udf has been removed , what should I do till fluent accepts many udfs not only one . thanks
Best regards
Sarah

Sarah,

Put all the UDF's in the same .c file. For instance:

#include "udf.h"

DEFINE_SOURCE()
{}

DEFINE_PROFILE()
{}

etc..

All of the appropriate functions will be available in the drop-down lists.

ComputerGuy

Dear ComputerGuy
Thank you so much! , however there is still a question for me , you declared that I use the same code which written for porosity for viscous and inertial resistances , but It is too hard since in my project these are defined as:
C1= (Dp^2 *porosity^3)/(150(1-porosity)^3)
C2=(3.5*(1-porosity))/(Dp*porosity^3))



So can I write a udf for them like this:
#include "udf.h"
DEFINE_PROFILE(inertial_res,t,i)

{
cell_t c;
begin_c_loop(c,t)
{
F_PROFILE(c,t,i) = 3.5*(1 - C_POR(c,t))
/(d_p*pow(C_POR(c,t),3));
}
end_c_loop(c,t)
}
C_POR(c,t) is porosity , so I have to modify porosity code of course I don’t know how do it! Could you help me?
Thanks
Sarah

Sarah,

You need to use a similar format to the porosity function I wrote, not the same one.


Short of writing all the UDF's for you, here's what I'd suggest:

1. Enable a user-defined memory location
2. In your DEFINE_PROFILE function for porosity (already written below in this thread), also write the porosity to the C_UDMI(c,t,0)
3. Create two new functions: one for C1 and one for C2 (DEFINE_PROFILEs)
4. Using the definitions you've just written, define C1 and C2, but instead of using C_POR, which I didn't know existed as a built-in definition, you could simply use C_UDMI(c,t,0), because C_UDMI(c,t,0) now represents the porosity in the cell.

Are you trying to simulate a packed bed? If so, have you read the Fluent manual (specificially, look at equation 7.19-15 and 7.19-16) for the Ergun and Blake-Kozeny equations.

Regards,
ComputerGuy

Dear computerGuy
Yes , I am supposed to simulate a packed bed and I’ve read Fluent manual which related to porous media modeling , I ‘ve taken your advices and wrote this code:
#define C_UDMI(c,t,0) F_PROFILE(c,t,i)
DEFINE_PROFILE(porosity_profile,t,i)
{
real x[ND_ND];
real y;
cell_t c;
begin_c_loop(c,t)
{
C_CENTROID(x,c,t);
y=x[1];
if((y==1.) || (y==1.85))
{
F_PROFILE(c,t,i)=0.39;
}
if ((y>1.) && (y<=1.425))
{
F_PROFILE(c,t,i)=0.39*(1.+1.05*exp(-100.*(y-1.)/0.06));
}
if ((y>1.425) && (y<1.85))
{
F_PROFILE(c,t,i)=0.39*(1.+1.05*exp(-100.*(1.85-y)/0.06));
}

}

end_c_loop(c,t)
}
// inertial resistance udf .


DEFINE_PROFILE(inertial_res,t,i)

{
cell_t c;
begin_c_loop(c,t)
{
F_PROFILE(c,t,i) = 3.5*(1 - C_UDMI(c,t,0))
/(0.06*pow(C_UDMI(c,t,0),3));
}
end_c_loop(c,t)
}
//viscous resistance udf .

DEFINE_PROFILE(viscous_res,t,i)
{
cell_t c;
begin_c_loop(c,t)
{
F_PROFILE(c,t,i) = 150*pow((1 - C_UDMI(c,t,0)),2)
/(0.0036*pow(C_UDMI(c,t,0),3));
}
end_c_loop(c,t)
}


Is it correct? Of course I could hook it to fluent like porosity code which I had written my self but you modified it , really you suggested to enable a user-defined memory location, why?

Best regards
Sarah

Sarah,

I think everything looks right, except for the first line. Try this:

By doing it this way, you'll write the value of porosity into a user defined memory location (which you also have to enable in Fluent, by the way). I have never seen C_POR(c,t), and cannot find it in the documentation, and thus I figured this is a safe way of referencing the cell porosity. I don't do a whole lot of UDF coding in porous regions, so perhaps you've discovered a keyword I don't know.

Let us know if this works for you!

ComputerGuy

Thank Computer guy
 

i gone through u r suggestions..The second point was not clear for me..can u please describe in detail??

Manohar:

1) Numerical integration: http://en.wikipedia.org/wiki/Numerical_integration
2) Custom field functions in Fluent: http://my.fit.edu/itresources/manual...g/node1221.htm

If you have an analytical representation of Peq as a function of temperature, you will be able to (without a UDF) have Fluent calculate it for you (#2). You can, under the contours of the solution, examine every custom field function. You can also output the value of the custom field function using a surface or volume monitor. Then, using #1, you can derive a value for C as a function of time.

Regards and Happy New Year,
ComputerGuy

ِDear ComputerGuy
I tried it but there are some errors when I want to initialize ,of course when I hook another udfs except udfs of viscous and inertial resistances , there is not any error and its done , I think maybe its error is related to C_UDMI(c,t,0), in fluent ,in fluid panel there are 2 directions of viscous and inertial resistances of course I give the same value to both of them while I assume porosity is constant , I saw C_POR in a literature which related to porous media !
Can you some advice about it
Happy New Year
cheers

This is unchecked. It doesn't use C_UDMI, and is not as efficient as it could be. However, I think it does what you want. I have looked through the Fluent manual and can't find C_POR, so if you must use it, I'm afraid you're on your own. Let us all know where you find it and if it works.

ComputerGuy

Dear ComputerGuy
I owe thanks to you whose advice and guidance have helped me so much in my project .
Really I had to change your code a little , I just eliminated : C_UDMI(c,t,0)=F_PROFILE(c,t,i) , F_PROFILE(c,t,i)=c_2, F_PROFILE(c,t,i)=c_1; from your code andafter that it ran.
I don’t know why with these while in initialization , fluent reports error But without these there is not any problem! By the way I don’t have any bias to use C_POR ! I just saw it in a literature .
Thank you and wish the best for you .
Sarah

Sarah,

I'm glad it worked. I left C_UDMI(c,t,0)=F_PROFILE(c,t,i); in the code accidentally, and you're right to take it out. However, in the C1 and C2 code, if you take out the line F_PROFILE(c,t,i)=c_2; and F_PROFILE(c,t,i)=c_1;, Fluent will not change your viscous and inertial resistances, and the result will probably not be correct according to the physics you're trying to impose. Make sure you double-check your results. Also, make sure you change Dp; I have assumed a particle diameter of 100 microns, but surely your particles are of a different size.

Regards,
ComputerGuy

Dear ComputerGuy
You are right but in your code , inside loops you indicated c_1 and c_2 instead of F_PROFILE , in inverse of this statement yes you are right . really how can I chek the result of udfs ? and If I want to plot temp of particles of porous medium how can I specify its zone surface? Yes, I changed the value of Dp to 0.06 .
Best regards
Sarah

Hi everyone,

I'm trying simulation porous media in rectangular channel, but the result isn't suitable with any research.

So, would you help me. I wish someone can check my simulation and give some reports if there is something wrong.

Thank you for your help.
Please send your e-mail, than i will send you my works to to your email.
my e-mail: oky.andytya.net@gmail.com

Regrads,
OKY Andytya P

note:
I use ANSYS Fluent 6.3 [CFD]

Hi Computer Guy,
I have a problem about my UDF code.
You mentioned that "Fluent doesn't automatically track what you're looking for, you could always create a few User-Defined Memory Locations to track values from a previous time step." I checked the UDF MANUAL and found that F_UDMI and C_UDMI store the face value and cell value respectively. If I want to store a calculation result from the previous step, can I use this Macros?
Here's part of my code. Should i use C_UDMI or F_UDMI?
And for the C_T_M1(c,t), it can only be used with UDS, how can I define a scalar?

Best Regards

Molixu,

First: You have posted this question in several threads. My suggestion, to help those who are trying to help you, is to please stop posting the same question to multiple threads.

Now, to answer your question about cumulative heat. Try the code below. Before you start your simulation, run the On Demand macro UDMInit. You need to have 1 user defined memory location activated in Fluent. I haven't tried this code, so there may be syntax errors, but it should work with minor modifications.

ComputerGuy

Really sorry for posting the same question..I was too worried about my problem...:confused: Thank you for your reply.
What I need is the total heat obtained in the domain of each step, acummulating them until a certain value and than change the boundary condition. I wrote the codes as follows, write Q_tot to a txt file and read it before using . But I found that in the txt file, the figure is always zero. and the BC didn't act in the way I want. Could you please give me some advice? Thanks a lot!

Inlet mass flow dependent on pressure at previous time step
 

Hello,

I am trying to simulate a flow in which the inlet mass flow rate depends on the pressure at previous time step in the domain. I tried implementing this in a udf using C_P_M1(c, t). But interpretation of udf always leads to segmentation violation.

Can someone suggest a way to implement this using udf.

Thanks,

Varun.

Just my 2 ct. to confirm that C_POR(c,t) does exist in order to retrieve the porosity for cells. You can find its definition in mem.h.

Hope this could help.

CFD Simulation of Gas turbine Engine_ Combustion Chamber
 

Hello everyone,

Can anyone suggest me a good document/book to start with simulation in gas turbine engine (combustion chamber).

I am planning to do the post processing in FLUENT. Is Openfoam software better than Fluent for combustion purpose?

Also, I intend to use C language to develop the code. Does Fluent take code written in Fortran. And which one is better? Fortran or C?

Please help me as I'm a beginner in this field.

I am planning to code advance combustion model and then link it to fluent/openfoam to simulate combustion chamber in a gas turbine.

I have certain parts of the code written in Fortran. Is there anyway I can directly use the fortran code in FLUENT?

Do I need to convert all the FORTRAn code into C language to be able to use it in FLUENT? Any method?

Please help.

2D position dependent heat source
 

hi,
I need to write an UDF for 2D heat source (on a face).The source has a parabolic position dependent & time independent function (source=x*x+2*x+1 ). Slould I use C-CENTROID or F-CENTROID in my programm?What is wrong with my programm?
With thanks & regards.
#include "udf.h"
DEFINE_SOURCE(energy_source,c,t,dS,eqn)
{
real x[ND_ND];
real source;
cell_t c;
begin_c_loop(c,t)
{
C_CENTROID(x,c,t);
source = x[0]*x[0]+2*x[0]+1;
dS[eqn] = 2*x[0]+2;
return source;
}
end_c_loop(c,t)
}

gas turbine combustor
 

Dear juzer

i think you could able simulate this model, in OpenFoam, more easy than fluent, and OF is much better than fluent.
your question is wrong, that you say witch one is better c or fortran
both are language for writing program,
if you want to use fluent you should write udf function, basically it use c.
but Openfoam is much simpler than all.

i myself now working of modeling one dimensional gas turbine combustor and i must use FORTRAN, have you ever seen ATEC code?

many regards.
Reza khodadadi.

problem with implementing porosity variation as a function of distance from the wall
 

Hi all,

It looks like you were able to solve a similar problem earlier. Now I am also trying to implement porosity variation as a function of distance from the wall for a packed bed reactor. My UDF has separate functions to define profiles for interfacial heat transfer coefficient, porosity and the specific surface area. I've pasted my code herewith.

I am using this .c file to load for the profiles of porosity, heat transfer coefficient and the specific surface area in FLUENT 14 with a thermal non-equilibrium model selected for the porous medium. I seem to be getting an error message when I load the porosity profile function.

It says - error - invalid argument [1] wrong type. Again, the .c file compiles and builds without any issues but pops this up when I try to use the porosity function..

Would you by any chance know what the bug is? Any help would be much appreciated.

/************************************************** ***
This is a UDF to calculate heat transfer coefficient and interfacial area
for the packed bed reactor example provided in FLUENT so that the 2 temperature
non-equilibrium model can be put to use

************************************************** ******/
#include "udf.h"

#define d_b 0.008 /* mean diameter of particles */
real A_V_sphere = 6.0/d_b; /* area-to-volume ratio of a sphere*/


DEFINE_PROFILE(Heat_trans_coeff,t,i)
{
cell_t c;
real Nu,Re,Pr;
real dens, visc; /*Fluid*/
real cond, cp; /* Fluid*/
real eps_b; /*Porosity of the bed*/

begin_c_loop(c,t)
{
eps_b = C_POR(c,t); /*Porosity fluid*/
dens = C_R(c,t); /*Density of fluid*/
visc = C_MU_L(c,t); /*Viscosity fluid*/
cond = C_K_L(c,t); /*Conductivity fluid*/
cp = C_CP(c,t); /*Specific heat fluid*/
Re = ND_MAG(C_U(c,t),C_V(c,t),C_W(c,t))*d_b*dens*eps_b/visc ; /*Reynolds number has been calculated based on pore velocity*/
Pr = cp*visc/cond ;
Nu = 2.+ 1.1 * pow(Re,0.6) * pow(Pr,1./3.);
F_PROFILE(c,t,i) = Nu*cond/d_b ;
}
end_c_loop(c,t)
}



DEFINE_PROFILE(Interfacial_Area_Density,t,i)
{
cell_t c;
real eps_b;
real A_b_V; /* area-to-volume ratio or the required specific surface area */
begin_c_loop(c,t)
{
eps_b = C_POR(c,t);
F_PROFILE(c,t,i) = (1.-eps_b)*A_V_sphere;
}
end_c_loop(c,t)
}


DEFINE_PROFILE(porosity_function,t,i)
{
cell_t c;
real x[ND_ND]; /*This will hold the position vectors*/
real y;
real a1;
real a2=6; /*to specifiy the porosity variation function*/
real eps_inf=0.37;
real eps;
a1 = (1./eps_inf)-1;


begin_c_loop(c,t)
{
C_CENTROID(x,c,t);
y = (H - x[1])/ d_b ;
eps = eps_inf*(1. + a1*exp(-1*a2*y));
F_PROFILE(c,t,i) = eps;
}
end_c_loop(c,t)
}

no value given for variable "x".
DEFINE_SOURCE(energy_source,c,t,ds,eqn)

{
real x=1.0;
ds[eqn]=0.3*cos(0.3*x);
return sin(0.3*x) ;
}

Hi,
I'm writing an udf to assign heat flux profile on 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 flux (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 z,glas_tube_emsivity,T_amb,T_sky,T_glass,T_gas,rad _loss,convec_loss;
double a,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)
}