[dhaya400]

Hello,

I am working on a unsteady flow problem with multiphase mixture model and run the simulation in Parallel mode.

For my problem, I have written a udf to update the pressure value of face boundary at every time step using the mass flow of the previous time step.

Code:
#include"udf.h"

real Q12; /* current time step flow rate */
#define BCA_ID 63 /* domain id*/
DEFINE_EXECUTE_AT_END(flow_rate_bca)
{
#if !RP_HOST
Domain *d;
real m;
cell_t c;
Thread *t;
face_t f;
d = Get_Domain(1);
t= Lookup_Thread(d, BCA_ID);
/* loop to calculate mass flow rate at outlet */
begin_f_loop(f,t)
{
if (PRINCIPAL_FACE_P(f,t))
{
m += F_FLUX(f,t);
}
}
end_f_loop(f,t)
Q12 = m / 1060.0; /* convert mass flow rate to volume flow rate */
#if RP_NODE
/* Perform node synchronized actions here. Does nothing in Serial */
Q12 = PRF_GRSUM1(Q12);
#endif /* RP_NODE */
m = 0.0;
#endif /* !RP_HOST */
}

DEFINE_PROFILE(pressure_out_bca,t,i)
{
#if !RP_HOST
real Z = 5.1918E+07;
real C = 8.697E-10;
real R = 10.6080E+08;
real ts = CURRENT_TIMESTEP;
face_t f;
begin_f_loop(f,t)
{
if (CURRENT_TIME == 0)
{
F_UDMI(f,t,0) = 13332.2;
F_UDMI(f,t,1) = F_UDMI(f,t,0);
F_UDMI(f,t,2) = 1.22E-05;
}

/* F_UDMI(f,t,0) - previous time step pressure
F_UDMI(f,t,1) - current time step pressure
F_UDMI(f,t,2) - previous time step flow rate
Q12 - current time step flow rate
*/

else
{
F_UDMI(f,t,3) = Q12;
F_UDMI(f,t,1) = ((F_UDMI(f,t,0)*R*C)+(F_UDMI(f,t,3)*R*ts)+(F_UDMI( f,t,3)*Z*ts)+(Z*R*C*F_UDMI(f,t,3))-(Z*R*C*F_UDMI(f,t,2))) / (ts +(R*C));
F_UDMI(f,t,0) = F_UDMI(f,t,1);
F_UDMI(f,t,2) = F_UDMI(f,t,3);
}
}
end_f_loop(f,t)
begin_f_loop(f,t)
{
F_PROFILE(f,t,i) = F_UDMI(f,t,1); /* Calculate pressure and apply to outlet*/
}
end_f_loop(f,t)
#endif /* !RP_HOST */
}


Question:
When I physically calculate the pressure values and check with pressure output files from simulation, it is not matching. Is there any error in the code which is causing this problem?

[pakk]

It could also be an error in your own calculation...

How big is the difference? Which values do you get, which values did you calculate? Which equation do you think you are using?

[dhaya400]

Hi Pakk,

Thank you for your reply.

I am using the following equation to calculate the pressure values for the outlet.

P(n+1) = {[P(n)*R*C] + [Q(n+1)*R*ts] + [Q(n+1)*Z*ts] + [Z*R*C*Q(n+1)] - [Z*R*C*Q(n)]} / (ts +(R*C))

F_UDMI(f,t,0) corresponds to P(n)
F_UDMI(f,t,1) corresponds to P(n+1)
F_UDMI(f,t,2) corresponds to Q(n)
F_UDMI(f,t,3) corresponds to Q(n+1)

Initializing, Q(n) = 1.22E-05 , P(n) = 13332.4 for CURRENT_TIME = 0

Result from simulation:
Q (n) Q(n+1) P(n) P(n+1)
1.22E-05 3.266E-05 13332.2 133324.4
3.266E-05 8.902E-05 133324.4 10423.9

Result from physical calculation:
Q (n) Q(n+1) P(n) P(n+1)
1.22E-05 3.266E-05 13332.2 14629.7
3.266E-05 8.902E-05 14629.7 18429.7

[AlexanderZ]

you may write to console each term of your equation using Message0("text"); macro

code seems to be correct.
the only thing I would do -> split time step definition with applying the value:
was
real ts = CURRENT_TIMESTEP;
to be
real ts;
ts = CURRENT_TIMESTEP;

[pakk]

You update F_UDMI(f,t,1) before you update F_UDMI(f,t,0), which seems wrong. Does not explain your current problem, however.

[dhaya400]

Hi Alexander and Pakk,

Thank you for your reply.
I have incorporated your suggestions and made few changes in the code for updating the previous time step values.
The following code is moved out of face loop and made separate face loop for the code.
F_UDMI(f,t,0) = F_UDMI(f,t,1);
F_UDMI(f,t,2) = F_UDMI(f,t,3);

Using Message command, I have checked the values of F_UDMI variables with physical calculation, the pressure values are matching. But if I check F_UDMI values with pressure values from the report files it is not matching. Why?
For example, F_UDMI = 16733 Pa, Pressure from report file = 16084 Pa

I want to update the pressure values for mixture domain only. But while printing the values, I come to know that the udf is updating for all the domains. So how to make the define profile macro to update mixture domain alone?

Revised code:

#include"udf.h"
real Q12; /* current time step flow rate */
#define BCA_ID 63 /* domain id*/
DEFINE_EXECUTE_AT_END(flow_rate_bca)
{
#if !RP_HOST
Domain *d;
real m;
cell_t c;
Thread *t;
face_t f;
d = Get_Domain(1);
t= Lookup_Thread(d, BCA_ID);
/* loop to calculate mass flow rate at outlet */
begin_f_loop(f,t)
{
if (PRINCIPAL_FACE_P(f,t))
{
m += F_FLUX(f,t);
}
}
end_f_loop(f,t)
Q12 = m / 1060.0; /* convert mass flow rate to volume flow rate */
#if RP_NODE
/* Perform node synchronized actions here. Does nothing in Serial */
Q12 = PRF_GRSUM1(Q12);
#endif /* RP_NODE */
m = 0.0;
#endif /* !RP_HOST */
}

DEFINE_INIT(init_pressure_bca,d)
{
Thread *t1;
cell_t c;
t1 = Lookup_Thread(d,BCA_ID); // for BCA
real P = 13332.2;
thread_loop_c(t1, d)
{
begin_c_loop_all(c,t1)
{
C_P(c,t1) = P;
}
end_c_loop_all(c, t1)
}
#if RP_HOST
Message("BCA Pressure Initalised\n");
#endif
}

DEFINE_PROFILE(pressure_out_bca,t,i)
{
#if !RP_HOST
real Z = 5.1918E+07;
real C = 8.697E-10;
real R = 10.6080E+08;
real ts;
ts = CURRENT_TIMESTEP;
face_t f;
begin_f_loop(f,t)
{
if (CURRENT_TIME == 0)
{
F_UDMI(f,t,0) = 13332.2;
F_UDMI(f,t,1) = F_UDMI(f,t,0);
F_UDMI(f,t,3) = 1.22E-05;
}
/* F_UDMI(f,t,0) - previous time step pressure
F_UDMI(f,t,1) - current time step pressure
F_UDMI(f,t,2) - prevoius time step flow rate
Q12 - current time step flow rate
*/
else
{
F_UDMI(f,t,3) = fabs (Q12);
F_UDMI(f,t,1) = ((F_UDMI(f,t,0)*R*C)+(F_UDMI(f,t,3)*R*ts)+(F_UDMI( f,t,3)*Z*ts)+(Z*R*C*F_UDMI(f,t,3))-(Z*R*C*F_UDMI(f,t,2))) / (ts +(R*C));
}
}
end_f_loop(f,t)
begin_f_loop(f,t)
{
F_UDMI(f,t,0) = F_UDMI(f,t,1);
F_UDMI(f,t,2) = F_UDMI(f,t,3);
}
end_f_loop(f,t)
begin_f_loop(f,t)
{
F_PROFILE(f,t,i) = F_UDMI(f,t,1); /* Calculate pressure and apply to outlet*/
}
end_f_loop(f,t)
#endif /* !RP_HOST */
}

[dhaya400]

Hi,

While searching for multiphase -Define profile macro, I cam across the below line from the udf manual:

In multiphase cases a DEFINE_PROFILE UDF may be called more than once (particularly if the profile is used in a mixture domain thread). If this needs to be avoided, then add the prefix MP_ to the UDF name. The function will then be called only once even if it is used for more than one profile.

What does this line covey? Can anyone explain with example?
I think.. if I got the answer to this question. I will solve the error in the above post.