CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > ANSYS > FLUENT > Fluent UDF and Scheme Programming

UDF: How to load data from an external file into a vector

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

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   June 16, 2016, 07:18
Default UDF: How to load data from an external file into a vector
  #1
New Member
 
Join Date: Aug 2015
Posts: 13
Rep Power: 11
EmiS is on a distinguished road
Hello,


I would like to impose inlet velocity values for unsteady flow by a udf. I have 2 files, one containing a single column of time values, and the other containing a column of corresponding velocity values (let's assume that at given time, velocity is uniform across the inlet surface).


I managed to make this work by creating two large arrays with time and velocity values hard-coded directly in the udf source code. Now I'm trying to have the data read from an external file, and I'm stuck. Ideally, I would like to read the values only once at startup, and then keep them in memory. Can this be done?
I'm trying to use the macro DEFINE_EXECUTE_ON_LOADING(name, libname), but unfortunately I have no idea what 'libname' should be.
I aslo thought about using a macro DEFINE_RW_FILE (from UDF docuemntation) but I don't manage to store the data from my external file into a vector "visible" for all the other macros of the UDF (especially the DEFINE_PROFILE)

I pasted a piece of my UDF, which returns to me that the "symbol velocityExp is not defined".

On a more general note, does Fluent and the UDF system support linking to external libraries? Suppose I have my own code written in C/C++ which is fairly large and does some complex computation to obtain boundary values. Is there a way of compiling this separately and then linking it to udf which would only call one high-level function (something like (get_my_bc(time, x, y, z) ....)?


Any working example code to load files on startup would be greatly appreciated.

Thank you very much for your help.


Emilie

/*----------------------------------------------------------------------------*/
/* UDF code - Exctract */
/* */
/*----------------------------------------------------------------------------*/

extern real velocityExp[90];


void read_data_velocity(FILE *ptr_velfile)
{
int i;

ptr_velfile = fopen("ExpData.txt", "r");
printf(":%s:\n", "Read input file for Velocity!");

for (i=0;i<90;i++)
{
fscanf (ptr_velfile, "%lf", &velocityExp[i]);
}

fclose(ptr_velfile);
printf ( "%lu\n", sizeof velocityExp );
printf(":%s:\n", "Close input file for Velocity!");
}


/*----------------------------------------------------------------------------*/

DEFINE_RW_FILE(reader_velocity, ptr_velfile)
{

Message0("Reading UDF data from data file...\n");

read_data_velocity(ptr_velfile);

host_to_node_real(velocityExp, 90);
}

/*----------------------------------------------------------------------------*/
EmiS is offline   Reply With Quote

Old   November 19, 2017, 13:57
Red face Same problem
  #2
FSM
New Member
 
Felipe
Join Date: Nov 2017
Location: Brazil
Posts: 16
Rep Power: 8
FSM is on a distinguished road
Hello, Emis

I am facing the same problem. I am trying to change velocity magnitude and direction of the inlet flow at some specific times. I also built a large code imputing this directly, and I am trying to do this input reading a .txt file with the velocity values. I am using the function fopen to read the file and a for loop to store the values on a vector. I used logical statements to access the values velocity[i] at the specific times, related to the CURRENT_TIME enter.
The fluent interpret the UDF with no error, however, when I try to initialize to run, the fluent faces a fatal error, and closes.
Does anyone know what might wrong?
Emis, did you figure out what happened to your code?

float velocity[4];

FILE *vel_x;
vel_x = fopen("teste.txt","r");

i=1;
for (i=1;i<=4;i++)
{
fscanf (vel_x, "%f", velocity[i]); /* read one value per line */
}
fclose(vel_x); /*close file*/
FSM is offline   Reply With Quote

Old   November 19, 2017, 19:30
Default
  #3
Senior Member
 
Alexander
Join Date: Apr 2013
Posts: 2,363
Rep Power: 34
AlexanderZ will become famous soon enoughAlexanderZ will become famous soon enough
Nice question,
I guess, something wrong with your UDF

Best regards
AlexanderZ is offline   Reply With Quote

Old   November 20, 2017, 12:45
Red face
  #4
FSM
New Member
 
Felipe
Join Date: Nov 2017
Location: Brazil
Posts: 16
Rep Power: 8
FSM is on a distinguished road
Thank you for the answers,

I also guess something is wrong with the code. Probably is on the vector part, But i still don't figuere it out what the problem is. I am posting the code, the test is a .txt file with 5 lines and one different value for each line.
If anyone could help me with that, would be very helpful.

thank you. regards.


#include "udf.h"


/************************************************** */
/*X COMPONENT*/
/************************************************** */

DEFINE_PROFILE(inlet_x_velocity, thread, nv)
{
real x[ND_ND]; /* Position vector*/
real y;
real z;
face_t f;
real current_timestep; /*Atual timestep size*/
real num_timestep;
real current_time;
int n_timestep;
int i, t=10;

//////////////////////////////////////////////////////
extern float velocity[5];

FILE *vel_x; /*Read the velocty file*/
//vel_x = fopen("uvelocity.txt","r"); /*fopen("caminho para o arquivo","modo r (read)")*/
vel_x = fopen("test.txt","r");

/*LOOP TO READ THE FILE AND STORE IN A VECTOR */
i=1; /* initialize before to loop */
for (i=1;i<=5;i++)
{
fscanf (vel_x, "%f", velocity[i]); /* read one value per line */
}
fclose(vel_x); /*close file*/
/////////////////////////////////////////////////////

n_timestep = N_TIME;

current_time = CURRENT_TIME;
current_timestep = CURRENT_TIMESTEP;
num_timestep = N_TIME;



begin_f_loop(f, thread)
{
F_CENTROID(x,f,thread);
y = x[1];

if(num_timestep<t)
{
F_PROFILE(f, thread, nv)=velocity[1];

}
else if(current_time<4.)
{
F_PROFILE(f, thread, nv)=velocity[2];

}
else if(current_time<6.)
{
F_PROFILE(f, thread, nv)=velocity[3];

}
else if(current_time<8.)
{
F_PROFILE(f, thread, nv)=velocity[4];

}
else
{
F_PROFILE(f, thread, nv)=velocity[5];
}
}
end_f_loop(f, thread)
}
FSM is offline   Reply With Quote

Old   November 20, 2017, 21:14
Default
  #5
Senior Member
 
Alexander
Join Date: Apr 2013
Posts: 2,363
Rep Power: 34
AlexanderZ will become famous soon enoughAlexanderZ will become famous soon enough
First check, how your input file is read, successfully?
Code:
/*LOOP TO READ THE FILE AND STORE IN A VECTOR */
i=1; /* initialize before to loop */
for (i=1;i<=5;i++)
{
fscanf (vel_x, "%f", velocity[i]); /* read one value per line */
Message(" velocity[%d] = %f\n",i,velocity[i]);
}
fclose(vel_x); /*close file*/
This code is suitable for single core only.
Also, you may put read_from_file part into other function for instance DEFINE_ON_DEMAND

And why do you use
Code:
 extern float velocity[5];  ?
Code:
real velocity[5];
is enough

On the other hand, if you need dynamic array you may use this code
Code:
real *velocity;
int n;
n = 5; // for your case, or any other value
velocity = (real*)malloc((n)*sizeof(real));
Best regards
AlexanderZ is offline   Reply With Quote

Old   November 21, 2017, 04:14
Default
  #6
Senior Member
 
Join Date: Nov 2013
Posts: 1,965
Rep Power: 26
pakk will become famous soon enough
The problem: "velocity[5]" does not exist! You start counting at 1, but you should start at 0, so use velocity[0] to velocity[4].
pakk is offline   Reply With Quote

Old   November 21, 2017, 08:22
Default
  #7
FSM
New Member
 
Felipe
Join Date: Nov 2017
Location: Brazil
Posts: 16
Rep Power: 8
FSM is on a distinguished road
Thank you all for the help.

I used all of the suggestion, but I still having trouble with this code.
I tested the malloc() to create a dynamic array, (velocity = (real*)malloc((n)*sizeof(real)), however I faced the same problem. The fluent reads the code and interpret without any problem, but when I initialize to run, FLUENT just crash and close.
I am pretty sure that the problem is with the reading part and storage in an array.
AlexanderZ suggested using the DEFINE_ON_DEMAND to read_from_file, could AlexanderZ or someone please clarify me about that.
  • Should I put the DEFINE_ON_DEMAND before the DEFINE_PROFILE?
  • how does it work?
I am new using UDF, and there is a lot to learn. If you guys could give me some ideas would be very helpful.

I am trying to use this code to change the component of the velocity with time. The code below is just a test, my objective is to change the components at more time intervals.

CODE:
Code:
#include "udf.h"

/***********COMPONENT X ****************************/
DEFINE_PROFILE(inlet_x_velocity, thread, nv)
{
	real x[ND_ND];	 /* Position vector*/
	real y;
	real z;
	face_t f;
	real current_time;
	int t=2, i=0;
	int n=5;
	real *velocity;
	
  	FILE *vel;
	vel = fopen("test.txt","r");
		
	velocity = (real*)malloc((n)*sizeof(real));
	
	/*store the data from file .txt to vector*/
	for (i=0;i<=4;i++)
	{
		fscanf (vel, "%f", velocity[i]); /* read one value per line */
		Message(" velocity[%d] = %f\n",i,velocity[i]);
		printf ("\n%f", velocity[i]);
	}
	fclose(vel); /*close file*/

    
	current_time = CURRENT_TIME;	
	
//////////////////////////////////	
	begin_f_loop(f, thread)
	{
	F_CENTROID(x,f,thread);
	y = x[1];

	if(current_time<t)
	{
		F_PROFILE(f, thread, nv)=velocity[1];
	}
	else if(current_time<4.)
	{
		F_PROFILE(f, thread, nv)=velocity[2];
	}
	else if(current_time<6.)
	{
		F_PROFILE(f, thread, nv)=velocity[3];	
	}
	else if(current_time<8.)
	{
		F_PROFILE(f, thread, nv)=velocity[4];
	}
	else
	{
		F_PROFILE(f, thread, nv)=velocity[5];
	}		
}
end_f_loop(f, thread)
}
The Message(" velocity[%d] = %f\n",i,velocity[i]); is also presenting a problem, this error appears on the Fluent console when I try to use it.
ERROR: line 76: function "CX_Message" not found (pc=97).

Thanks in advance.
Best regards,

Felipe.
FSM is offline   Reply With Quote

Old   November 21, 2017, 22:51
Default
  #8
Senior Member
 
Alexander
Join Date: Apr 2013
Posts: 2,363
Rep Power: 34
AlexanderZ will become famous soon enoughAlexanderZ will become famous soon enough
Do you compile this code or interpret it?

Some functions do not supported for interpretation, so you should compile this code.

Best regards
AlexanderZ is offline   Reply With Quote

Old   November 22, 2017, 05:00
Default
  #9
Senior Member
 
Join Date: Nov 2013
Posts: 1,965
Rep Power: 26
pakk will become famous soon enough
The same problem remains!
You write data to
  • velocity[0]
  • velocity[1]
  • velocity[2]
  • velocity[3]
  • velocity[4]
and later you read data from
  • velocity[1]
  • velocity[2]
  • velocity[3]
  • velocity[4]
  • velocity[5]
This is wrong! You try to read velocity[5], which does not exist, and that gives you the fatal error.
pakk is offline   Reply With Quote

Old   December 7, 2017, 07:20
Question
  #10
FSM
New Member
 
Felipe
Join Date: Nov 2017
Location: Brazil
Posts: 16
Rep Power: 8
FSM is on a distinguished road
Thank you all very much for the clarifications.

I made some corrections on the code and now I finally figured out how to read the data file and store the values into a vector. I used a DEFINE_ON_LOAD macro just for the reading process. I am new using UDF and I am still learning to deal with those different UDF macros.

However, I am facing another problem working with the loops. I need to input each value at a specific time according to the timestep. For example in my test case, the interval between each velocity value is 20 sec, which means 200 timesteps (each timestep=0.1s). Furthermore, between one velocity value and the next, I will need to do a linear interpolation. I am trying to use a WHILE loop (I am not sure if it is the better way) to do so, but I am getting a fatal error on fluent.
The code is compiled with no problems and I can create the library libudf with no errors. The error just shows up when I try to Initialize the case to run.
I think this is more of a programming problem. Programming is not one of my strongest ability.
If you guys could help me to fix my code I would really appreciate.

Code:
#include "udf.h"

/*Global variables*/
real *uvelocity;
int *time;
int n=6;  /*number of lines with data on the file*/


DEFINE_EXECUTE_ON_LOADING(on_load,libudf)								/*Read the .txt file and store the velocity values into a vector*/
{
	int i=1;
	float data_vel;
	int data_time;
	FILE *vel;
	FILE *tm;
	vel = fopen("test.txt","r");
	tm = fopen("time.txt","r");
	
	uvelocity = (real*)malloc((n)*sizeof(real));
	time = (int*)malloc((n)*sizeof(int));

/*Places the data into the vector*/
	for (i=1;i<=n;i++)
	{
		data_vel=0;
		data_time=0;
		fscanf (vel, "%f", &data_vel); 	
		fscanf (tm,"%d", &data_time);								/* read one value per line */
		uvelocity[i]=data_vel;
		time[i]=data_time;
		Message("* velocity_U[%d] = %f\n",i,uvelocity[i]);
		Message("* time[%d] = %d\n",i,time[i]);
	}
	fclose(vel); 															
}

DEFINE_ON_DEMAND(list_velocity)										/*check if the values are right on the vector*/
{
	int i=1,j=1;
	for(i=1;i<=n;i++)
	{
		Message("\nVelocity_U[%d]=%f\n",i,uvelocity[i]);
		Message("\nTime[%d]=%d\n",i,time[j]);
		j++;
	}
}

DEFINE_PROFILE(inlet_x_velocity, thread, nv)							/*Define profile to the x velocity component*/
{
	real x[ND_ND];	 
	real y;
	face_t f;
	int n_ts;
	real current_time;
	int  j=1;
	int tL=1, tR=199, deltat=200; 				/* tL= left time;  tR = right Time after deltat timesteps, deltat = number of timesteps between tL and tR */

	current_time = CURRENT_TIME;										
	n_ts = N_TIME;								/*Actual number of timestep*/						
	
//////////////////////////////////	
	begin_f_loop(f, thread)
	{					
	while(n_ts <= 1000)					/*loop to set the velocity values and interpolate between a value and the next one*/ 
	{
		
			F_CENTROID(x,f,thread);
			y = x[1];
			
			if(n_ts==tL)
			{
				F_PROFILE(f, thread, nv)=uvelocity[j];
			}
			else if(n_ts>tL && n_ts<=tR)
			{
				F_PROFILE(f, thread, nv)=((n_ts-tL)*(uvelocity[j+1]-uvelocity[j])/(tR-tL))+uvelocity[j];	
			}
				j=j+1;
				tL=tL+deltat;
				tR=tR+deltat;
		
	}	
	end_f_loop(f, thread)
	}	
}

Data set
Velocity data: File test.txt
Code:
5.0
3.0
6.0
4.0
6.0
0
Thanks in Advance.
FSM is offline   Reply With Quote

Old   December 11, 2017, 03:32
Default
  #11
Senior Member
 
Join Date: Nov 2013
Posts: 1,965
Rep Power: 26
pakk will become famous soon enough
It is the same problem as before.

Code:
uvelocity = (real*)malloc((n)*sizeof(real));
If n=6, this means you have the following memory locations:
  • uvelocity[0]
  • uvelocity[1]
  • uvelocity[2]
  • uvelocity[3]
  • uvelocity[4]
  • uvelocity[5]

But later in your code, you write to them like this:
Code:
for (i=1;i<=n;i++)
{
...
uvelocity[i]=data_vel;
...
}
This means that you write to:

  • uvelocity[1]
  • uvelocity[2]
  • uvelocity[3]
  • uvelocity[4]
  • uvelocity[5]
  • uvelocity[6]

uvelocity[6] does not exist!

In general: start counting at 0, not at 1!
pakk is offline   Reply With Quote

Old   January 1, 2019, 01:02
Default liner interpolation
  #12
New Member
 
vektor
Join Date: Aug 2018
Posts: 23
Rep Power: 8
vek123 is on a distinguished road
so i read the comments, and i am pretty sure that i am bad at these UDF thing, and c programming in general



i created a simple material udf ( only viscosity for now) and, what i did is basically just create two arrays, one for temperature and the other for viscosity, followed by a simple algorithm for the linear interpolation,
it goes like this



Code:
DEFINE_PROPERTY(cell_viscosity, cell, thread)
{
  real mu_lam; 

                    /* local variable  */
  real temp = C_T(cell, thread);   /* local variable  */
 int i; int size_int;

 real T_int{ 4.57E+01, 7.21E+02, 1.94E+03, 3.19E+03, 4.05E+03, 5.35E+03, 6.21E+03, 7.15E+03, 8.90E+03, 9.76E+03,
 1.05E+04, 1.18E+04, 1.31E+04, 1.45E+04, 1.66E+04, 1.83E+04, 2.06E+04, 2.39E+04, 2.68E+04, 3.00E+04};

 real eta_T{8.24E+02, 3.71E+03, 6.94E+03, 9.94E+03, 1.33E+04, 1.68E+04, 1.92E+04, 2.20E+04, 2.50E+04, 
 2.60E+04, 2.58E+04, 2.27E+04, 1.62E+04, 9.65E+03, 4.12E+03, 2.41E+03, 1.82E+03, 1.76E+03, 1.53E+03, 1.00E+03}; /* M.capitelli et al,*/

 size_int = sizeof(T_int);

 for i(i=1;i<= size_int;i++)
{
  if (T_int[1]>temp || T_int[size_int])
    break;
  else if (T_int[i]<temp && T_int[i+1]>temp)
    mu_lam = (eta_T[i]*(T_int[i+1]-temp)+eta_T[i+1]*(temp-T_int[i]))/(T_int[i+1]-T_int[i]);

  return mu_lam;
}
i still don't if this will work, but i wish if someone more familiar with C and udfs to point out any issues, also i want to do the same for,
density, electric and thermal conductivity and the specific heat capacity.
thanks
vek123 is offline   Reply With Quote

Old   January 1, 2019, 20:44
Default
  #13
Senior Member
 
Alexander
Join Date: Apr 2013
Posts: 2,363
Rep Power: 34
AlexanderZ will become famous soon enoughAlexanderZ will become famous soon enough
Code:
 for i(i=1;i<= size_int;i++)
{
  if (T_int[1]>temp || T_int[size_int])
    break;
  else if (T_int[i]<temp && T_int[i+1]>temp)
    mu_lam = (eta_T[i]*(T_int[i+1]-temp)+eta_T[i+1]*(temp-T_int[i]))/(T_int[i+1]-T_int[i]);}
if (T_int[1]>temp || T_int[size_int]) what is it? How it should work from your point of view?
break; is not appropriate in DEFINE_PROPERTY macro, material properties have to be defined at the whole range of simulation (temperature).
else if (T_int[i]<temp && T_int[i+1]>temp), if temperature is lower than T[i] but higher than T[i+1] ?? is it possible?

best regards
AlexanderZ is offline   Reply With Quote

Old   January 1, 2019, 21:58
Default
  #14
New Member
 
vektor
Join Date: Aug 2018
Posts: 23
Rep Power: 8
vek123 is on a distinguished road
Quote:
Originally Posted by AlexanderZ View Post
Code:
 for i(i=1;i<= size_int;i++)
{
  if (T_int[1]>temp || T_int[size_int])
    break;
  else if (T_int[i]<temp && T_int[i+1]>temp)
    mu_lam = (eta_T[i]*(T_int[i+1]-temp)+eta_T[i+1]*(temp-T_int[i]))/(T_int[i+1]-T_int[i]);}
if (T_int[1]>temp || T_int[size_int]) what is it? How it should work from your point of view?
break; is not appropriate in DEFINE_PROPERTY macro, material properties have to be defined at the whole range of simulation (temperature).
else if (T_int[i]<temp && T_int[i+1]>temp), if temperature is lower than T[i] but higher than T[i+1] ?? is it possible?

best regards
hello AlexanderZ,
i really appreciate your help, i really just made that code without trying it, that break i put it in case (i am not sure about this) if the temperature in the model, somehow became smaller than 50[K] or higher that 30000[K] in which my material properties are not defined, so i put break so that the software will immediately exit the loop without wasting time looking for something that is none existent in the first place, i know that is not the way to do it.


EDIT: it should be if (T_int[1]>temp || T_int[size_int]<temp) i forgot to add <temp



in the next statement:


else if (T_int[i]<temp && T_int[i+1]>temp)


this will go through the loop from i=1 to the end and compare temp
which i assume is the temperature solved by fluent ( also not sure about this) with the T_int which is the temperature from material properties table, it will stop at certain i in which

T_int[i]<temp<T_int[i+1] and after that i used the linear interpolation formula. i got the formula from wikipedia
https://en.wikipedia.org/wiki/Linear_interpolation


EDIT2:
Quote:
else if (T_int[i]<temp && T_int[i+1]>temp), if temperature is lower than T[i] but higher than T[i+1] ?? is it possible?
no it means temp is lower than T_int[i+1] and higher than T_int[i]




so what i am thinking/ or what i want to do in the UDF is:


1) read the model temperature from the solver at each cell.
2) compare it the tempature from the material table
3) compute the new material viscosity based on table data.
4) put in the model cell.


repeat int the next time step.
this is the original UDF that i modified on:

/************************************************** *******************
UDF that simulates solidification by specifying a temperature-
dependent viscosity property
************************************************** ********************/

DEFINE_PROPERTY(cell_viscosity, c, t)
{
real mu_lam;
real temp = C_T(c, t);

if (temp > 288.)
mu_lam = 5.5e-3;
else if (temp > 286.)
mu_lam = 143.2135 - 0.49725 * temp;
else
mu_lam = 1.;

return mu_lam;
}
vek123 is offline   Reply With Quote

Old   January 1, 2019, 22:36
Default
  #15
New Member
 
vektor
Join Date: Aug 2018
Posts: 23
Rep Power: 8
vek123 is on a distinguished road
i removed the break thing and simply put it this way, i added the possibility if temp is exactly the temperature in the table.


DEFINE_PROPERTY(cell_viscosity, cell, thread)
{
real mu_lam;

/* local variable */
real temp = C_T(cell, thread); /* local variable */
int i; int size_int;

real T_int{ 4.57E+01, 7.21E+02, 1.94E+03, 3.19E+03, 4.05E+03, 5.35E+03, 6.21E+03, 7.15E+03, 8.90E+03, 9.76E+03,
1.05E+04, 1.18E+04, 1.31E+04, 1.45E+04, 1.66E+04, 1.83E+04, 2.06E+04, 2.39E+04, 2.68E+04, 3.00E+04};

real eta_T{8.24E+02, 3.71E+03, 6.94E+03, 9.94E+03, 1.33E+04, 1.68E+04, 1.92E+04, 2.20E+04, 2.50E+04,
2.60E+04, 2.58E+04, 2.27E+04, 1.62E+04, 9.65E+03, 4.12E+03, 2.41E+03, 1.82E+03, 1.76E+03, 1.53E+03, 1.00E+03}; /* M.capitelli et all,*/

size_int = sizeof(T_int);

for i(i=1;i<=size_int;i++)
{

if(temp == T_int[i])
mu_lam = eta_T[i];
else if (T_int[i]<temp && T_int[i+1]>temp)
mu_lam = (eta_T[i]*(T_int[i+1]-temp)+eta_T[i+1]*(temp-T_int[i]))/(T_int[i+1]-T_int[i]);

return mu_lam;
}
vek123 is offline   Reply With Quote

Old   January 2, 2019, 02:21
Default
  #16
Senior Member
 
Alexander
Join Date: Apr 2013
Posts: 2,363
Rep Power: 34
AlexanderZ will become famous soon enoughAlexanderZ will become famous soon enough
Code:
if(temp == T_int[i])
mu_lam = eta_T[i];
this temp == T_int[i] approach is bad, because what is the chance you will have exactly T_int[i] temperature?

Code:
DEFINE_PROPERTY(cell_viscosity, cell, thread)
{
real mu_lam; 

/* local variable */
real temp = C_T(cell, thread); /* local variable */
int i; int size_int;

real T_int{ 4.57E+01, 7.21E+02, 1.94E+03, 3.19E+03, 4.05E+03, 5.35E+03, 6.21E+03, 7.15E+03, 8.90E+03, 9.76E+03,
1.05E+04, 1.18E+04, 1.31E+04, 1.45E+04, 1.66E+04, 1.83E+04, 2.06E+04, 2.39E+04, 2.68E+04, 3.00E+04};

real eta_T{8.24E+02, 3.71E+03, 6.94E+03, 9.94E+03, 1.33E+04, 1.68E+04, 1.92E+04, 2.20E+04, 2.50E+04, 
2.60E+04, 2.58E+04, 2.27E+04, 1.62E+04, 9.65E+03, 4.12E+03, 2.41E+03, 1.82E+03, 1.76E+03, 1.53E+03, 1.00E+03}; /* M.capitelli et all,*/

size_int = sizeof(T_int);

for i(i=1;i<=size_int;i++)
{

if(temp > T_int[size_int])
mu_lam = eta_T[size_int];
else if (T_int[i]<temp && T_int[i+1]>temp)
mu_lam = (eta_T[i]*(T_int[i+1]-temp)+eta_T[i+1]*(temp-T_int[i]))/(T_int[i+1]-T_int[i]);
else
mu_lam = eta_T[0]; // true if temperature is equal or lower than T_int[0]

return mu_lam;
}
I didn't test this code

best regards
AlexanderZ is offline   Reply With Quote

Old   January 2, 2019, 02:53
Default
  #17
New Member
 
vektor
Join Date: Aug 2018
Posts: 23
Rep Power: 8
vek123 is on a distinguished road
Quote:
Originally Posted by AlexanderZ View Post
Code:
if(temp == T_int[i])
mu_lam = eta_T[i];
this temp == T_int[i] approach is bad, because what is the chance you will have exactly T_int[i] temperature?

Code:
DEFINE_PROPERTY(cell_viscosity, cell, thread)
{
real mu_lam; 

/* local variable */
real temp = C_T(cell, thread); /* local variable */
int i; int size_int;

real T_int{ 4.57E+01, 7.21E+02, 1.94E+03, 3.19E+03, 4.05E+03, 5.35E+03, 6.21E+03, 7.15E+03, 8.90E+03, 9.76E+03,
1.05E+04, 1.18E+04, 1.31E+04, 1.45E+04, 1.66E+04, 1.83E+04, 2.06E+04, 2.39E+04, 2.68E+04, 3.00E+04};

real eta_T{8.24E+02, 3.71E+03, 6.94E+03, 9.94E+03, 1.33E+04, 1.68E+04, 1.92E+04, 2.20E+04, 2.50E+04, 
2.60E+04, 2.58E+04, 2.27E+04, 1.62E+04, 9.65E+03, 4.12E+03, 2.41E+03, 1.82E+03, 1.76E+03, 1.53E+03, 1.00E+03}; /* M.capitelli et all,*/

size_int = sizeof(T_int);

for i(i=1;i<=size_int;i++)
{

if(temp > T_int[size_int])
mu_lam = eta_T[size_int];
else if (T_int[i]<temp && T_int[i+1]>temp)
mu_lam = (eta_T[i]*(T_int[i+1]-temp)+eta_T[i+1]*(temp-T_int[i]))/(T_int[i+1]-T_int[i]);
else
mu_lam = eta_T[0]; // true if temperature is equal or lower than T_int[0]

return mu_lam;
}
I didn't test this code

best regards
That's true actually there is a very small chance of it happening. And your edit looks spot on. I will try that one.
Thanks
I am still far away from running the final simulation " electric arc simulation". Since it is a bit complicated but i may come back with some more question about udfs during the process.

Thank you so much.
vek123 is offline   Reply With Quote

Old   January 14, 2019, 01:55
Default
  #18
New Member
 
vektor
Join Date: Aug 2018
Posts: 23
Rep Power: 8
vek123 is on a distinguished road
Quote:
Originally Posted by AlexanderZ View Post
Code:
if(temp == T_int[i])
mu_lam = eta_T[i];
this temp == T_int[i] approach is bad, because what is the chance you will have exactly T_int[i] temperature?

Code:
DEFINE_PROPERTY(cell_viscosity, cell, thread)
{
real mu_lam; 

/* local variable */
real temp = C_T(cell, thread); /* local variable */
int i; int size_T_int;

real T_int{ 4.57E+01, 7.21E+02, 1.94E+03, 3.19E+03, 4.05E+03, 5.35E+03, 6.21E+03, 7.15E+03, 8.90E+03, 9.76E+03,
1.05E+04, 1.18E+04, 1.31E+04, 1.45E+04, 1.66E+04, 1.83E+04, 2.06E+04, 2.39E+04, 2.68E+04, 3.00E+04};

real eta_T{8.24E+02, 3.71E+03, 6.94E+03, 9.94E+03, 1.33E+04, 1.68E+04, 1.92E+04, 2.20E+04, 2.50E+04, 
2.60E+04, 2.58E+04, 2.27E+04, 1.62E+04, 9.65E+03, 4.12E+03, 2.41E+03, 1.82E+03, 1.76E+03, 1.53E+03, 1.00E+03}; /* M.capitelli et all,*/

size_int = sizeof(T_int);

for i(i=1;i<=size_int;i++)
{

if(temp > T_int[size_int])
mu_lam = eta_T[size_int];
else if (T_int[i]<temp && T_int[i+1]>temp)
mu_lam = (eta_T[i]*(T_int[i+1]-temp)+eta_T[i+1]*(temp-T_int[i]))/(T_int[i+1]-T_int[i]);
else
mu_lam = eta_T[0]; // true if temperature is equal or lower than T_int[0]

return mu_lam;
}
I didn't test this code

best regards

it had some small problems, i fixed it and it works


Code:
#include "udf.h"
DEFINE_PROPERTY(cell_viscosity, cell, thread)
{
real mu_lam; 

/* local variable */
real temp = C_T(cell, thread); /* local variable */
int i; int size_T_int;

real T_int[] = { 4.57E+01, 7.21E+02, 1.94E+03, 3.19E+03, 4.05E+03, 5.35E+03, 6.21E+03, 7.15E+03, 8.90E+03, 9.76E+03,1.05E+04, 1.18E+04, 1.31E+04, 1.45E+04, 1.66E+04, 1.83E+04, 2.06E+04, 2.39E+04, 2.68E+04, 3.00E+04};

real eta_T[] = {8.24E+02, 3.71E+03, 6.94E+03, 9.94E+03, 1.33E+04, 1.68E+04, 1.92E+04, 2.20E+04, 2.50E+04, 2.60E+04, 2.58E+04, 2.27E+04, 1.62E+04, 9.65E+03, 4.12E+03, 2.41E+03, 1.82E+03, 1.76E+03, 1.53E+03, 1.00E+03}; /* M.capitelli et all,*/

int size_T_int = sizeof(T_int)/sizeof(int);

for (i=0;i<=size_T_int-1;i++)
{

if(temp > T_int[size_T_int-1])
mu_lam = eta_T[size_T_int-1];
else if (T_int[i]<temp && T_int[i+1]>temp)
mu_lam = (eta_T[i]*(T_int[i+1]-temp)+eta_T[i+1]*(temp-T_int[i]))/(T_int[i+1]-T_int[i]);
else
mu_lam = eta_T[0]; // true if temperature is equal or lower than T_int[0]
}
return mu_lam;
}
vek123 is offline   Reply With Quote

Reply

Tags
boundary condition u, global vector, load text file, udf

Thread Tools Search this Thread
Search this Thread:

Advanced Search
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 Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
[Other] How to use finite area method in official OpenFOAM 2.2.0? Detian Liu OpenFOAM Meshing & Mesh Conversion 4 November 3, 2015 03:04
[swak4Foam] Problem installing swak_2.x for OpenFoam-2.4.0 towanda OpenFOAM Community Contributions 6 September 5, 2015 21:03
Trouble compiling utilities using source-built OpenFOAM Artur OpenFOAM Programming & Development 14 October 29, 2013 10:59
[swak4Foam] Error bulding swak4Foam sfigato OpenFOAM Community Contributions 18 August 22, 2013 12:41
[blockMesh] Axisymmetrical mesh Rasmus Gjesing (Gjesing) OpenFOAM Meshing & Mesh Conversion 10 April 2, 2007 14:00


All times are GMT -4. The time now is 04:32.