[Revolver2006]

Dear all,


In my project I use momentum and energy sources in a dedicated cell zone to simulate the presence of a fan. In order to set the total pressure jump across the fan automatically, I am writing a UDF. The idea is to compute the massflow-averaged total pressure at the inflow and outflow boundary of the cell zone, which I will use to adjust the source terms accordingly. However, I am running into issues, probably with the thread pointing; Here is a piece of the code (in progress):


DEFINE_SOURCE(mass_source, c, t, dS, eqn)
{
real source;

static int i = 0;

static real area;
static real product;
static real massflow;

real x[ND_ND];
real NV_VEC(A);

face_t f;

real u;
real v;
real w;
real V;

real p_operating;
real pressure;

real f_area;
real f_mdot;
real f_product;

real density;
real ptot;

Thread *tb_inflow;
Thread *tb_outflow;

p_operating = RP_Get_Real("operating-pressure");

/*Get the domain*/
Domain *domain = Get_Domain(1);

/* Get inflow boundary thread, face id is 110 */
tb_inflow = Lookup_Thread(domain, 110);

/* Get inflow boundary thread, face id is 111 */
tb_outflow = Lookup_Thread(domain, 111);

/* Loop over all the cell faces in the thread of outflow boundary */
begin_f_loop(f, tb_outflow)
{
if PRINCIPAL_FACE_P(f, tb_outflow)
{
i++;

/* Point to the specific thread */
f_area = F_AREA(A, f, tb_outflow);
f_mdot = F_FLUX(f, tb_outflow);

pressure = F_P(f, tb_outflow);
density = F_R(f, tb_outflow);

u = F_U(f, tb_outflow);
v = F_V(f, tb_outflow);

V = sqrt(pow(u, 2) + pow(v, 2));

ptot = pressure + 1 / 2 * density*pow(V, 2);

f_product = ptot * f_mdot;

massflow += f_mdot;
area += f_area;
product += f_product;
}
}
end_f_loop(f, tb_outflow)

#if !RP_HOST
pt_avg_out = product / massflow;
#endif

source = 1;
dS[eqn] = 0;
return source;
}

The error appears to stem from inserting the tb_outflow thread pointer into the F_U, F_V and F_R macros. Just using F_R(f,t) for example does work, but the numbers do not make any sense. Furthermore, the output generated appears dependent on the amount of parallel processes in use:


real pt_avg_out;
real pt_avg_in;

DEFINE_EXECUTE_AT_END(execute_at_end)
{
#if !RP_HOST
printf("Mass-averaged total pressure at fan exit plane is %g\n", pt_avg_out);
fflush(stdout);
#endif
}


I use in Fluent (18.2) with a coupled compressible pressure-based solver for an axis-symmetric domain. As mentioned, the domain consists of 2 cell zones, one for the fluid domain and one for the fan. Probably I am making some obvious mistakes, but I am running out of ideas (and time for that matter). Any help would be appreciated. thanks!

[obscureed]

Hi Revolver2006,
Some quick comments:

(1) Any time you calculate a sum across cells or faces in parallel, you need to share it across the compute nodes to get the total sum:

#if RP_NODE
  sum1 = PRF_GRSUM(sum1); /* or PRF_GISUM1 for integer */
#endif

(See here for an example of the equivalent calculation of a maximum value: Is it possible? Changing BC from wall to pressure-outlet with an UDF trigger) The host node still won't know the answer, but that does not matter here.

(2) The help files suggest that F_U and F_V are not normally stored on interior face zones -- and I guess that includes face zones between two fluid cell zones. If you really need them, you'll have to take the averages of the two neighbouring cells' values.

(3) It might not matter in a small 2D case, but it is really wasteful to do this whole looping calculation on every visit to DEFINE_SOURCE. Better to put it into a DEFINE_EXECUTE_AT_END and store the result in a static variable.

(4) I accept that this is work in progress, but the current UDF is called "mass_source", which is hopefully not the final plan. I don't think you'd need an energy source for most purposes -- a bit of kinetic energy here and there is not enough to disturb temperature significantly.

(5) Have you checked out the built-in "fan" boundary condition? Another crude way to get something like a fan is to set a fixed value of (for example) x-velocity in the thin fan-like cell zone. (Remember that fixed values work only in the pressure-based solver -- so you're OK in this instance.)

(6) Your compiler will almost certainly help you out, but "1/2" is a classic way to get the integer "0" by mistake (because integer arithmetic rounds fractions down). Your compiler will probably help you out by noticing that these are involved in a real calculation, so it will convert to 1./2. and possibly to 0.5. Much safer to use 0.5 when that's what you mean.

Good luck!
Ed

[Revolver2006]

Dear Ed,


Thank you very much for your comments!

If have managed to get the script working and it is giving me good results. Indeed parallelization of the code was really the key to make it work, together with the fact that (aside from pressure) no flow variables are defined at the boundaries.

In order to reduce the computational load of the simulations, I execute the main code only once every 10 iterations. This also helps to stabilize the flow field when the source term is adjusted.

Regarding the method, I use the source term (momentum and energy) approach because I am dealing with compressible flow for an aerospace application. The built-in fan boundary condition was not giving satisfactory results. Comparison with actual fan data shows good agreement with the current results. Thank you once again!


regards,
Martijn