[SHAFKAT91]

Hi,

I am just a beginner with ANSYS. I am using a UDF function to define three properties of a material: Specific heat, density and thermal conductivity. when I use density and thermal conductivity as used defined and specific heat as a constant, fluent runs. But when I try to use specific heat as user defined, following page shows:

Node 0: Process 10516: Received signal SIGSEGV.
----------------------------------------------------------------

Node 1: Process 12556: Received signal SIGSEGV.
---------------------------------------------------------------

MPI Application rank 0 exited before MPI_Finalize() with status 2
--------------------------------------------------------------------------
The fl process could not be started.
--------------------------------------------------

fatal error during one of the computation process.




My code is following:

#include "udf.h"

DEFINE_PROPERTY(cell_conductivity,cell,thread)
{
real f;
real kl;
real ks;
real k;
kl = 0.148;
ks= 0.42;
f= C_LIQF(cell,thread);
k= kl*f+((1-f)*ks);
return k;
}
DEFINE_PROPERTY(density,cell,thread)
{
real f;
real rol;
real ros;
real ro;
rol=789;
ros=840;
f= C_LIQF(cell,thread);
ro= rol*f+((1-f)*ros);
return ro;
}
DEFINE_PROPERTY(spesific_heat,cell,thread)
{
real f;
real cpl;
real cps;
real cp;
cpl=2460;
cps= 1920;
f= C_LIQF(cell,thread);
cp= cpl*f+((1-f)*cps);
return cp;
}

Can anyone please suggest me what this is happening?

[pakk]

Does this happen directly after you initialize, or later in the simulation?

If it happens directly after you initialize, the only thing I can imagine is that in Fluent in liquid fraction (C_LIQF) is calculated after the specific heat is calculated and before conductivity and density are calculated.

In that case, in your calculation of specific heat, you are using a value that is not stored yet, so Fluent does not know what to do.

If that is all true, a solution would be to initialize with a constant specific heat, and change to a user-specified specific heat after that.

[KevinZ09]

Specific heat has its own macro, i.e., DEFINE_SPECIFIC_HEAT. Try using that one instead. I've never had a problem with it.

[pakk]

Originally Posted by KevinZ09

Specific heat has its own macro, i.e., DEFINE_SPECIFIC_HEAT. Try using that one instead. I've never had a problem with it.

I agree. Ignore my earlier advice, and do what KevinZ says.

[SHAFKAT91]

Hello,

Thank you both very much for advice. But I saw that specific heat requires T, Tref, h, yi. I am not working with any of them. Can I modify SPECIFIC_HEAT for calculating the average of solid and liquid portion in a cell.

[KevinZ09]

You should be able to do that. You can modify the DEFINE_SPECIFIC_HEAT marco according to your needs, sort of like you were doing before in the DEFINE_PROPERTY macro that you used. However, you do have to use the specific heat macro, you can't use a regular property macro (as you experienced). So use your mixture formula and just make sure you return the correct cp.

[pakk]

But the DEFINE_SPECIFIC_HEAT macro does not use cell information. How can he calculate C_LIQF(cell,thread) if he does not know "cell"?

[KevinZ09]

Originally Posted by pakk

But the DEFINE_SPECIFIC_HEAT macro does not use cell information. How can he calculate C_LIQF(cell,thread) if he does not know "cell"?

You are right about that, there's no way of accessing cell threads in this macro. You can only have a temperature dependent cp.

What are you precisely trying to do? Do you have a multicomponent flow? Because if so, I think you can use the mixing-law for the mixture cp. That calculates a mass-weighted average if my memory serves me correctly (I used it a while ago, so not entirely sure anymore). But yeah, perhaps give some more details so we can try to help you out further.

[SHAFKAT91]

I am actually trying a melting case. There is no use of h. The melting temperature is 1000K. Once the melting starts, I want to calculate the specific heat in a cell by the following formula:

specific heat= liquid fraction*(Cp)liquid+ (1-liquid fraction)*(Cp)solid

I am doing the same thing using PROPERTY macros for density and thermal conductivity. But for specific heat, I have hit a wall.

[KevinZ09]

The enthalpy you just define it as cp(T-Tref). Look in the manual for the DEFINE_SPECIFIC_HEAT macro and you'll see how it works a bit.

And as said earlier, you can't use DEFINE_PROPERTY for specific heat, you have to use the aforementioned macro, because enthalpy and specific heat have to be related. So if you want to use an UDF for specific heat, you need to use that one.

But again, if you have a multicomponent model, you can set the cp in the mixture panel such that it uses a mass-weighted average I believe. So the question is, are you using a multicomponent model?

[ujwal rajan]

Hey guys!

Can someone help me??? [I am very new to this!]
I am facing the same problem. I am performing a FSI over an airfoil using a UDF. I was first asked to do this for a mock case with mesh already biult. My task now is to apply the same UDF for another mesh and airfoil.
The response file will give the pitch response, plunge response, lift and pitching moment

The code is given below:

#include "matrixOperations.h"
#define alpha RP_Get_Real("angle-of-attack")
#define alphas RP_Get_Real("pitch-freeplay")
real mass = 10.0;
real Salpha = 9.0;
real Ialpha = 17.4;
real Kh = 45000.0;
real Kalpha = 200000.0;
real hs = 0.0;
//real alphas = 0.0349066;
real Kh_hat = 0.0;
real Kalpha_hat = 0.0;
square_matrix *mmatrix = NULL;
square_matrix *minvmatrix = NULL;
vector *fvector = NULL;
vector *qvector = NULL;
vector *qdotvector = NULL;
FILE *response = NULL;
const size_t dof = 2;
real current_time;
void initmatrix()
{
current_time = -1.0;
mmatrix = square_matrix_alloc_zeros(dof);
mmatrix->data[0] = mass;
mmatrix->data[1] = Salpha;
mmatrix->data[2] = Salpha;
mmatrix->data[3] = Ialpha;
minvmatrix = square_matrix_alloc(dof);
square_matrix_inverse(mmatrix, minvmatrix);
fvector = vector_alloc_zeros(dof);
qvector = vector_alloc_zeros(dof);
qdotvector = vector_alloc_zeros(dof);
}
void f(real time, const vector * y, const vector * z, vector * fnResult)
{
square_matrix *kamatrix = square_matrix_alloc_zeros(dof);
square_matrix *kbmatrix = square_matrix_alloc_zeros(dof);
square_matrix *kcmatrix = square_matrix_alloc_zeros(dof);
vector *qamatrix = vector_alloc_zeros(dof);
vector *qbmatrix = vector_alloc_zeros(dof);
vector *qcmatrix = vector_alloc_zeros(dof);
qamatrix->data[0] = y->data[0];
qamatrix->data[1] = y->data[1];
qbmatrix->data[0] = hs;
qbmatrix->data[1] = alphas;
if (y->data[0] > hs)
{
kamatrix->data[0] = Kh;
kbmatrix->data[0] = -Kh;
kcmatrix->data[0] = Kh_hat;
qcmatrix->data[0] = pow((y->data[0] - hs), 3.0);
}
else if (y->data[0] < -hs)
{
kamatrix->data[0] = Kh;
kbmatrix->data[0] = Kh;
kcmatrix->data[0] = Kh_hat;
qcmatrix->data[0] = pow((y->data[0] + hs), 3.0);
}
else
{
kamatrix->data[0] = 0.0;
kbmatrix->data[0] = 0.0;
kcmatrix->data[0] = 0.0;
}
if (y->data[1] > alphas)
{
kamatrix->data[3] = Kalpha;
kbmatrix->data[3] = -Kalpha;
kcmatrix->data[3] = Kalpha_hat;
Message0("In one!\n");
qcmatrix->data[1] = pow((y->data[1] - alphas), 3.0);
}
else if (y->data[1] < -alphas)
{
kamatrix->data[3] = Kalpha;
kbmatrix->data[3] = Kalpha;
kcmatrix->data[3] = Kalpha_hat;
Message0("In two!\n");
qcmatrix->data[1] = pow((y->data[1] + alphas), 3.0);
}
else
{
kamatrix->data[3] = 0.0;
kbmatrix->data[3] = 0.0;
kcmatrix->data[3] = 0.0;
}
vector *fnterm0 = vector_alloc_zeros(dof);
vector *fnterm1a = vector_alloc_zeros(dof);
vector *fnterm1b = vector_alloc_zeros(dof);
vector *fnterm1c = vector_alloc_zeros(dof);
vector *fnResulta = vector_alloc_zeros(dof);
vector *fnResultb = vector_alloc_zeros(dof);
square_matrix *fnterm2a = square_matrix_alloc_zeros(dof);
square_matrix *fnterm2b = square_matrix_alloc_zeros(dof);
square_matrix *fnterm2c = square_matrix_alloc_zeros(dof);
square_matrix_vector_multi(minvmatrix, fvector, fnterm0);
square_matrix_matrix_multi(minvmatrix, kamatrix, fnterm2a);
square_matrix_matrix_multi(minvmatrix, kbmatrix, fnterm2b);
square_matrix_matrix_multi(minvmatrix, kcmatrix, fnterm2c);
square_matrix_vector_multi(fnterm2a, qamatrix, fnterm1a);
square_matrix_vector_multi(fnterm2b, qbmatrix, fnterm1b);
square_matrix_vector_multi(fnterm2c, qcmatrix, fnterm1c);
vector_sub(fnterm0, fnterm1a, fnResulta);
vector_sub(fnResulta, fnterm1b, fnResultb);
vector_sub(fnResultb, fnterm1c, fnResult);
square_matrix_free(kamatrix);
square_matrix_free(kbmatrix);
square_matrix_free(kcmatrix);
vector_free(qamatrix);
vector_free(qbmatrix);
vector_free(qcmatrix);
vector_free(fnterm0);
vector_free(fnterm1a);
vector_free(fnterm1b);
vector_free(fnterm1c);
vector_free(fnResulta);
vector_free(fnResultb);
square_matrix_free(fnterm2a);
square_matrix_free(fnterm2b);
square_matrix_free(fnterm2c);
}

void g(real time, const vector * y, const vector * z, vector * gnResult)
{
vector_scale(z, 1.0, gnResult);
}
void runge4(real t, vector * y, vector * z, real delta_t)
{
vector *r4k1 = vector_alloc_zeros(dof);
vector *r4k2 = vector_alloc_zeros(dof);
vector *r4k3 = vector_alloc_zeros(dof);
vector *r4k4 = vector_alloc_zeros(dof);
vector *r4l1 = vector_alloc_zeros(dof);
vector *r4l2 = vector_alloc_zeros(dof);
vector *r4l3 = vector_alloc_zeros(dof);
vector *r4l4 = vector_alloc_zeros(dof);
vector *r4temp0 = vector_alloc_zeros(dof);
vector *r4temp1 = vector_alloc_zeros(dof);
vector *r4temp2 = vector_alloc_zeros(dof);
g(t, y, z, r4temp0);
vector_scale(r4temp0, delta_t, r4k1);
f(t, y, z, r4temp0);
vector_scale(r4temp0, delta_t, r4l1);
vector_scale(r4k1, 0.5, r4temp0);
vector_add(y, r4temp0, r4temp1);
vector_scale(r4l1, 0.5, r4temp0);
vector_add(z, r4temp0, r4temp2);
g(t + delta_t * 0.5, r4temp1, r4temp2, r4temp0);
vector_scale(r4temp0, delta_t, r4k2);
f(t + delta_t * 0.5, r4temp1, r4temp2, r4temp0);
vector_scale(r4temp0, delta_t, r4l2);
vector_scale(r4k2, 0.5, r4temp0);
vector_add(y, r4temp0, r4temp1);
vector_scale(r4l2, 0.5, r4temp0);
vector_add(z, r4temp0, r4temp2);
g(t + delta_t * 0.5, r4temp1, r4temp2, r4temp0);
vector_scale(r4temp0, delta_t, r4k3);
f(t + delta_t * 0.5, r4temp1, r4temp2, r4temp0);
vector_scale(r4temp0, delta_t, r4l3);
vector_add(y, r4k3, r4temp0);
vector_add(z, r4l3, r4temp1);
g(t + delta_t * 0.5, r4temp0, r4temp1, r4temp2);
vector_scale(r4temp2, delta_t, r4k4);
f(t + delta_t * 0.5, r4temp0, r4temp1, r4temp2);
vector_scale(r4temp2, delta_t, r4l4);
vector_add(r4k2, r4k3, r4temp0);
vector_scale(r4temp0, 2.0, r4temp1);
vector_add3(r4k1, r4temp1, r4k4, r4temp0);
vector_scale(r4temp0, 1.0 / 6.0, r4temp1);
vector_add(y, r4temp1, r4temp0);
vector_scale(r4temp0, 1.0, y);
vector_add(r4l2, r4l3, r4temp0);
vector_scale(r4temp0, 2.0, r4temp1);
vector_add3(r4l1, r4temp1, r4l4, r4temp0);
vector_scale(r4temp0, 1.0 / 6.0, r4temp1);
vector_add(z, r4temp1, r4temp0);
vector_scale(r4temp0, 1.0, z);
vector_free(r4k1);
vector_free(r4k2);
vector_free(r4k3);
vector_free(r4k4);
vector_free(r4l1);
vector_free(r4l2);
vector_free(r4l3);
vector_free(r4l4);
vector_free(r4temp0);
vector_free(r4temp1);
vector_free(r4temp2);
}
DEFINE_ON_DEMAND(open_close)
{
if (response == NULL)
{
response = fopen("response.txt", "w");
}
else
{
fclose(response);
}
}
void update_force_vector()
{
Domain *domain = Get_Domain(1);
Thread *thread = Lookup_Thread(domain, 3);
face_t f;
real r0[ND_3];
r0[0] = 0.5*qvector->data[0]*sin(alpha*M_PI / 180.0);
r0[1] = qvector->data[0]*cos(alpha*M_PI / 180.0);
r0[2] = 0.0;
real area[ND_3], rf[ND_3], r[ND_3];
real force[ND_3], forceV[ND_3], forceP[ND_3],
forceTot[ND_3], moment, momentTot;
N3V_S(force, = , 0.0);
N3V_S(forceTot, = , 0.0);
N3V_S(forceV, = , 0.0);
N3V_S(forceP, = , 0.0);
moment = 0.0;
momentTot = 0.0;
real Lift = 0.0;
real Momt = 0.0;
N3V_S(r, = , 0.0);
N3V_S(rf, = , 0.0);
N3V_S(area, = , 0.0);
#if RP_NODE || !PARALLEL
begin_f_loop_int(f, thread)
{
F_AREA(area, f, thread);
F_CENTROID(rf, f, thread);
N3V_VV(r, = , rf, -, r0);
N3V_VS(forceV, = , F_STORAGE_R_N3V(f, thread, SV_WALL_SHEAR), *, -1.0);
N3V_VS(forceP, = , area, *, F_P(f, thread));
N3V_VV(force, = , forceV, +, forceP);
moment = r[0] * force[1] - r[1] * force[0];
N3V_V(forceTot, += , force);
momentTot += moment;
}
end_f_loop_int(f, thread)
int i;
for (i = 0; i < 3; i++)
{
forceTot[i] = PRF_GRSUM1(forceTot[i]);
}
momentTot = PRF_GRSUM1(momentTot);
Lift = cos(alpha*M_PI / 180.0)*forceTot[1] - sin(alpha*M_PI / 180.0)*forceTot[0];
Momt = momentTot;
#endif
node_to_host_real_2(Lift, Momt);
fvector->data[0] = Lift;
fvector->data[1] = Momt;
}
DEFINE_ON_DEMAND(storeCoord)
{
Domain *domain = Get_Domain(1);
Thread *t;
int n;
face_t f;
Node *v;
thread_loop_f(t, domain)
{
begin_f_loop(f, t)
{
f_node_loop(f, t, n)
{
v = F_NODE(f, t, n);
N_UDMI(v, 0) = NODE_X(v);
N_UDMI(v, 1) = NODE_Y(v);
}
}
end_f_loop(f, t)
}
}
DEFINE_GRID_MOTION(eqofmotion, domain, dynamic_thread, time, dtime)
{
if (mmatrix == NULL)
{
initmatrix();
}
if (time != current_time)
{
if (time < 0.0)
{
qvector->data[0] = 0.0;
qvector->data[1] = 0.01745329252*sin(100 * time);
qdotvector->data[0] = 0.0;
qdotvector->data[1] = 1.745329252*cos(100 * time);
}
else
{
update_force_vector();
runge4(time, qvector, qdotvector, dtime);
}
current_time = time;
}
if (response != NULL)
{
fprintf(response, "%14.6e %14.6e %14.6e %14.6e %14.6e\n", time, qvector->data[0], qvector->data[1], fvector->data[0], fvector->data[1]);
}
Thread *tf = DT_THREAD(dynamic_thread);
face_t f;
Node *v;
real dx[ND_ND];
int n;
SET_DEFORMING_THREAD_FLAG(THREAD_T0(tf));
Clear_Node_Flags(domain, NODE_VISITED_FLAG);
begin_f_loop(f, tf)
{
f_node_loop(f, tf, n)
{
v = F_NODE(f, tf, n);
if (NODE_POS_NEED_UPDATE(v))
{
NODE_POS_UPDATED(v);
dx[0] = ((N_UDMI(v, 0) + 0.5)*cos(qvector->data[1]) - N_UDMI(v, 1)*sin(qvector->data[1]) - 0.5) - (qvector->data[0] * sin(alpha*M_PI / 180.0));
dx[1] = ((N_UDMI(v, 0) + 0.5)*sin(qvector->data[1]) + N_UDMI(v, 1)*cos(qvector->data[1])) + (qvector->data[0] * cos(alpha*M_PI / 180.0));
NV_V(NODE_COORD(v), = , dx);
}
}
}
end_f_loop(f, tf)
}


When i try implementing the same UDF, it gives the following error:
MPI Application rank 0 exited before MPI_Finalize() with status 2
the fl process could not be started.

i tried changing the compiler directives to only :
#if !PARALLEL
{

}

When i do this, it seems to get past the error but only the time gets updated on the response file. The pitch response, plunge response, lift and pitching moment are always displayed as zero as time progresses.

[ujwal rajan]

Hey guys, I found the problem, I had to match the zone ID for the airfoil surface.

Sent from my ONEPLUS A3003 using CFD Online Forum mobile app

[Stac93]

Hi guys,

i am facing the same Problem. In my case, i have to define a temperature and pressure dependend viscosity and density
I interpreted the following udf:






/************************************************** *******************
UDF for viscosity and density
************************************************** ********************/
#include "udf.h"


DEFINE_PROPERTY(HLP46_density, c, t)
{
real rho;

real temp = C_T(c,t);
real pska, A0, A1, B, C, CP0, D, E1, E2, EO, G, RhoO, Toff;
real p = C_P(c,t);

pska = 10 / ((p / 100000)); /*[kg m^-1 s^-2])/(1 [kg m^-1 s^-2])*/
Toff = temp-273.15;
A0 = 0.640157 * 100;
A1 = 0.443242 * 10000; /*[K^-1]*/
B = 0.171801 * 10000;
C = 0.215003;
CP0 = (10^5)*(A0+A1*temp); /**[J kg^-1 K^-1]**/
D = -0.618199; /** [kg m^-3] **/
E1 = 0.115299 * 10000;
E2 = -0.134291 * 10000;
EO = 0.269774; /** [kg m^-3] **/
G = 0.001; /** [Pa s] **/
RhoO = 0.886860 * 1000; /**[kg m^-3] *10^3**/

rho = Toff*(EO*(1-(exp(-pska/E1)+E2*pska))+D)+RhoO/(1-C*log10((B+pska)/B));
return rho;
}
DEFINE_PROPERTY(HLP46_viscosity,c,t)
{
real mu_lam, Toff, G, a0, a1, a2, a3, a4, a5, a6, a7, a8, pska;
real T = C_T(c,t);
real p = C_P(c,t);

Toff = T - 273.15 ;/*(1[K])-1*/
mu_lam = G * exp((a0 - a2 * Toff) / (a1 + Toff)) * exp(((a3 + a4 * pska) * exp(-Toff / (a5 + a6 * pska)) + (a7 + a8 * pska)) * pska);
G = 0.001;
a0 = 0.716789 * 1000;
a1 = 0.116482 * 1000;
a2 = 0.361897 * 10;
a3 = 0.822870 / 100;
a4 = -0.311457 * 100000;
a5 = 0.570805 * 1000;
a6 = -0.228366;
a7 = -0.562193 / 100;
a8 = 0.300900 / 100000;
pska = 1 / ((p / 100000));

return mu_lam;
}



Right after initialisation, i get the error: The fl process could not be started. Therfore i Initialized it before i did set the user defined viscosity and density. But when i want to start the calculation then i get the folliwing error and fluent collapses:
999999: mpt_accept: Error: accept failed: Invalid argument

Please give me some tips why it doesn't work.
Thank you in advance

[AlexanderZ]

initialize first, then set material properties

/************************************************** *******************
UDF for viscosity and density
************************************************** ********************/
#include "udf.h"


DEFINE_PROPERTY(HLP46_density, c, t)
{
real rho;
real pska, A0, A1, B, C, CP0, D, E1, E2, EO, G, RhoO, Toff,p,temp ;

temp = C_T(c,t);
p = C_P(c,t);

pska = 10 / ((p / 100000)); /*[kg m^-1 s^-2])/(1 [kg m^-1 s^-2])*/
Toff = temp-273.15;
A0 = 0.640157 * 100; 
A1 = 0.443242 * 10000; /*[K^-1]*/
B = 0.171801 * 10000;
C = 0.215003;
CP0 = (10^5)*(A0+A1*temp); /**[J kg^-1 K^-1]**/
D = -0.618199; /** [kg m^-3] **/
E1 = 0.115299 * 10000;
E2 = -0.134291 * 10000;
EO = 0.269774; /** [kg m^-3] **/
G = 0.001; /** [Pa s] **/
RhoO = 0.886860 * 1000; /**[kg m^-3] *10^3**/

rho = Toff*(EO*(1-(exp(-pska/E1)+E2*pska))+D)+RhoO/(1-C*log10((B+pska)/B));
return rho;
}
DEFINE_PROPERTY(HLP46_viscosity,c,t)
{
real mu_lam, Toff, G, a0, a1, a2, a3, a4, a5, a6, a7, a8, pska,T, p;
T = C_T(c,t);
p = C_P(c,t);
Toff = T - 273.15 ;/*(1[K])-1*/ 

G = 0.001; 
a0 = 0.716789 * 1000;
a1 = 0.116482 * 1000;
a2 = 0.361897 * 10;
a3 = 0.822870 / 100;
a4 = -0.311457 * 100000;
a5 = 0.570805 * 1000;
a6 = -0.228366;
a7 = -0.562193 / 100;
a8 = 0.300900 / 100000;
pska = 1 / ((p / 100000));
mu_lam = G * exp((a0 - a2 * Toff) / (a1 + Toff)) * exp(((a3 + a4 * pska) * exp(-Toff / (a5 + a6 * pska)) + (a7 + a8 * pska)) * pska);
return mu_lam;
}

pressure must be different from 0

best regards

[pakk]

1. Always try to compile, don't interpret.
2. Find mistakes by looking at compile errors.
3. Find mistakes by making your code simpler
4. You calculate mu_lam before you calculate the constants that it depends on.
5. 10^5 does not mean 100000 in c.