[UchihaMadara]

Hello ,


I am trying to run one of those cases, where the properties from outlet from few surfaces are assigned as inlet to another.


In this setup, I intend to calculate total pressure and velocity on the surface, and assigned to another surface with the change in friction factor value. However I am not sure if I have set the the total pressure and velocity with the correct formula.


The massflow rates have been giving values in the correct neighborhood.
however the other two are currently problemsome.


TLR- Could anyone point out, how can I obtain the total pressure and velocity at a surface using UDF?



As I am using a constant density assumption, I am not obtaining this with F_R(f,t). Dia, mu, L are also constants.

DEFINE_PROFILE(pressure_outlet, thread, position)
{
    /* Averages static pressure on downstream interface inflow BC and applies it at the upstream outflow interface 
     */

    int ID;
    face_t f;
    real NV_VEC(area);

    real Amagnitude, pressure, press_avg;
    real A_sum, A_glob, velocity_glob, density_sum, density_glob ;
    real mdot_sum, mdot_glob, Flux, mflux_glob , reynolds, velocity, velocity_sum, friction_factor;
    real PA_glob;

    face_t f_upstream;

    Domain* domain_upstream;
    domain_upstream = Get_Domain(1);
    real density = 0.9462;
    int zone_ID = THREAD_ID(thread);     /* set the ID of the downstream interface */

        if (zone_ID == out1)
    {
        ID = in1;
    }
    else if (zone_ID == out2)
    {
        ID = in2;
    }
    else if (zone_ID == out3)
    {
        ID = in3;
    }
    else if (zone_ID == out4)
    {
        ID = in4;
    }
    else
    {
        ID = 10000;
    }

    if (ID!=10000)
    {        
        Thread* thread_upstream = Lookup_Thread(domain_upstream, ID);
        mdot_sum = 0.0;
        A_sum = 0.0; 
        velocity_sum = 0.0;
        reynolds = 0.0;
        
        begin_f_loop(f_upstream, thread_upstream)
            if PRINCIPAL_FACE_P(f_upstream, thread_upstream)
            {
                F_AREA(area, f_upstream, thread_upstream);
                Amagnitude += NV_MAG(area);
                Flux += F_FLUX(f_upstream, thread_upstream);
                pressure +=  F_P(f_upstream, thread_upstream); 
                 velocity  += F_FLUX(f_upstream, thread_upstream)/(density*NV_MAG(area));
                /*pressure = F_P(f_upstream, thread_upstream) +  0.5 * density * pow(velocity,2);*/
                /*PA_sum += pressure * NV_MAG(area);*/
                /*velocity_sum += velocity;*/
            }
        end_f_loop(f_upstream, thread_upstream)

        mdot_glob = PRF_GRSUM1(Flux    * Amagnitude);
        A_glob = PRF_GRSUM1(Amagnitude);
        mflux_glob = mdot_glob / A_glob; 

        PA_glob = PRF_GRSUM1(pressure*Amagnitude);

        velocity_sum = PRF_GRSUM1(velocity*Amagnitude);
        velocity_glob = velocity_sum / A_glob;  /* AWA velocity */

        reynolds = density * fabs(velocity_glob) * Dia / mu;

        friction_factor =  1/pow((1.8*log10(reynolds)-1.5),2);

        press_avg = (PA_glob/A_glob) + (0.5*density*pow(velocity_glob,2)) - (friction_factor*(L/Dia)*(0.5*density*pow(velocity_glob,2)))*frelax;
        
    }
    else
    {
        press_avg =  0.0;
    }
        begin_f_loop(f, thread)
        {            
            F_PROFILE(f, thread, position) = press_avg;
        }
        end_f_loop(f, thread)    
    Message0("reynolds at %d is %g\n", zone_ID, reynolds); 
    Message0("ff at at %d is %g\n", zone_ID, friction_factor);
    Message0("velocity %d is %g\n", zone_ID, velocity_glob);
    Message0("calculated_p at %d is %g\n", zone_ID, press_avg);
    Message0("\n");
}

Outputs looks as follows

reynolds at 51 is -7.5611e+07
ff at at 51 is -nan
velocity 51 is -290.34
calculated_p at 51 is -nan

reynolds at 50 is -2.11107e+08
ff at at 50 is -nan
velocity 50 is -810.635
calculated_p at 50 is -nan

reynolds at 49 is -1.34733e+08
ff at at 49 is -nan
velocity 49 is -517.363
calculated_p at 49 is -nan

reynolds at 48 is -8.88475e+07
ff at at 48 is -nan
velocity 48 is -341.167
calculated_p at 48 is -nan

[AlexanderZ]

to get ff you are using nothing else but reynolds number, which is defined.
so the problem here could be in log10()
try add #include <math.h> in header

also you are using
Amagnitude
Flux
pressure
velocity
which are not defined before summation, fix it. they are not zeros by default.

compile your code

[UchihaMadara]

Thank you for your response AlexenderZ,


Yes, I do have all that you have mentioned, was just showing the DEFINE_Profile bit

#include "udf.h"
#include <math.h>

#define in1 36
#define in2 35
#define in3 34
#define in4 33
#define out1 51
#define out2 50
#define out3 49
#define out4 48

#define frelax 0.0001

#define L 0.672
#define Dia 0.006
#define mu 2.18e-5
#define PI 3.14159
#define density 0.9462
/*
Step 1: integrate the mass flux at every outlet and set value on corresponding BC
Step 2: integrate the static pressure at every inlet and set value on corresponding BC
*/

DEFINE_PROFILE(massflux_inlet, thread, position)
{
    /* Integrates mass flux through upstream outflow interface
       Sets the mass flux profile on the downstream interface
     */

    int ID; /* the ID of zone with upstream boundary */
    face_t f;
    real area[ND_ND];

    real Amagnitude, Flux, mflux_glob;
    real mdot_sum, mdot_glob;
    real A_sum, A_glob;

    face_t f_upstream;
    Domain* domain_upstream;          /* domain is declared as a variable   */
    domain_upstream = Get_Domain(1);  /* returns fluid domain pointer       */

    /* Find the mass flux through the upstream interface outlet */
    /* set the ID of the upstream outlet interface surface */
    int zone_ID = THREAD_ID(thread);
        if (zone_ID == in1)
    {
        ID = out1;
    }
    else if (zone_ID == in2)
    {
        ID = out2;
    }
    else if (zone_ID == in3)
    {
        ID = out3;
    }
    else if (zone_ID == in4)
    {
        ID = out4;
    }
    else
    {
        ID = 10000;
    }
    if (ID != 10000)
    {
        Thread* thread_upstream = Lookup_Thread(domain_upstream, ID);
        Flux  = 0.0;
        Amagnitude = 0.0;
        begin_f_loop(f_upstream, thread_upstream)
            if PRINCIPAL_FACE_P(f_upstream, thread_upstream)
            {
                F_AREA(area, f_upstream, thread_upstream);
                Amagnitude += NV_MAG(area);
                Flux += F_FLUX(f_upstream, thread_upstream);
            }
        end_f_loop(f_upstream, thread_upstream)
        
            mdot_glob = PRF_GRSUM1(Flux*Amagnitude);
            A_glob = PRF_GRSUM1(Amagnitude);
            mflux_glob = mdot_glob / A_glob; 
    
    }
    else
    {
        mflux_glob = 0.0;
    }

    begin_f_loop(f, thread)
    {    
        F_PROFILE(f, thread, position) = mflux_glob;
    }
    end_f_loop(f, thread)
    Message0("Amagnitude is %f\n", Amagnitude);
    Message0("massflowrate at %d is %f\n", zone_ID, mflux_glob);    
}


DEFINE_PROFILE(pressure_outlet, thread, position)
{
    /* Averages static pressure on downstream interface (in:*) inflow BC and applies it at the upstream outflow interface  (out:*)
     */

    int ID; /* the ID of zone with downstream boundary */
    face_t f;
    real area[ND_ND];

    real Amagnitude, pressure, press_avg;
    real PA_sum, PA_glob;
    real A_sum, A_glob, velocity_glob, density_sum, density_glob ;
    real mdot_sum, mdot_glob, Flux, mflux_glob , reynolds, velocity, velocity_sum, friction_factor;

    face_t f_downstream;
    Domain* domain_downstream;         
    domain_downstream = Get_Domain(1); 
    
    /* set the ID of the downstream interface */
    int zone_ID = THREAD_ID(thread); /* applied to */

        if (zone_ID == out1)
    {
        ID = in1;
    }
    else if (zone_ID == out2)
    {
        ID = in2;
    }
    else if (zone_ID == out3)
    {
        ID = in3;
    }
    else if (zone_ID == out4)
    {
        ID = in4;
    }
    else
    {
        ID = 10000;
    }
    if (ID!=10000)
    {        
        Thread* thread_downstream = Lookup_Thread(domain_downstream, ID);
        Amagnitude = 0.0; 
        Flux = 0.0;
        pressure = 0.0;
        
        begin_f_loop(f_downstream, thread_downstream)
            if PRINCIPAL_FACE_P(f_downstream, thread_downstream)
            {
                F_AREA(area, f_downstream, thread_downstream);
                Amagnitude += NV_MAG(area);
                Flux += F_FLUX(f_downstream, thread_downstream);
                pressure +=  F_P(f_downstream, thread_downstream); 
            }
         end_f_loop(f_downstream, thread_downstream)

mdot_glob = PRF_GRSUM1(Flux*Amagnitude);
        A_glob  = PI * 0.25 * pow(Dia,2);
        mflux_glob = mdot_glob / A_glob;
        PA_glob = PRF_GRSUM1(pressure*Amagnitude);
        velocity_glob = mflux_glob / (density*A_glob);
        reynolds = density * fabs(velocity_glob) * Dia / mu; 
        friction_factor = 1.0/pow((1.8*log10(reynolds)-1.5),2);
        press_avg = (PA_glob/A_glob) + (0.5*density*pow(velocity_glob,2)) - (friction_factor*(L/Dia)*(0.5*density*pow(velocity_glob,2)));
    }
    else
    {
        press_avg =  0.0;
    }
        begin_f_loop(f, thread)
        {            
            /*velocity = mflux_glob/ (density*(PI*0.25*pow(Dia,2)));*/
            F_PROFILE(f, thread, position) =  press_avg;
        }
        end_f_loop(f, thread)    
    Message0("\n");
    Message0("Amagnitude is: %0.12f\n", Amagnitude);
    Message0("A_glob is %f\n", A_glob);
    Message0("PA_glob is: %f\n", PA_glob);
    Message0("massflowrate at in%d: is %f\n", zone_ID, mflux_glob);
    Message0("reynolds at %d is %f\n", zone_ID, reynolds);
    Message0("ff at at %d is %f\n", zone_ID, friction_factor);
    Message0("velocity %d is %f\n", zone_ID, velocity_glob);
    Message0("calculated_p at %d is %f\n", zone_ID, press_avg);
    Message0("\n");
 }

the output for which is something like. Here the main concern is that despite following the manual example for obtaining total area, I still get a zero, and perhaps due to it, the other values suffers. Would you have a suggesstion to correct Amagnitude? Idon't see why this should yeild 0.

Amagnitude is 0.000000000000
massflowrate at 36 is 0.000223
Amagnitude is 0.000000000000
massflowrate at 35 is 0.000283
Amagnitude is 0.000000000000
massflowrate at 34 is 0.000173
Amagnitude is 0.000000000000
massflowrate at 33 is 0.000176

Amagnitude is:            0
A_glob is 0.000028
PA_glob is: 0.538942
massflowrate at in51: is -0.000707
reynolds at 51 is 6883.657060
ff at at 51 is 0.034191
velocity 51 is -26.432700
calculated_p at 51 is 18125.940007


Amagnitude is:            0
A_glob is 0.000028
PA_glob is: 0.738375
massflowrate at in50: is -0.001392
reynolds at 50 is 13549.274851
ff at at 50 is 0.028366
velocity 50 is -52.028146
calculated_p at 50 is 23326.701371


Amagnitude is:            0
A_glob is 0.000028
PA_glob is: 0.009078
massflowrate at in49: is -0.001392
reynolds at 49 is 13549.275805
ff at at 49 is 0.028366
velocity 49 is -52.028150
calculated_p at 49 is -2466.904770


Amagnitude is:            0
A_glob is 0.000028
PA_glob is: -0.401181
massflowrate at in48: is -0.000925
reynolds at 48 is 9008.356102
ff at at 48 is 0.031680
velocity 48 is -34.591377
calculated_p at 48 is -15631.361372

[AlexanderZ]

you are not making global summation for area
you've tried, but didn't use it:
was

A_glob = PRF_GRSUM1(Amagnitude);

as you are using Amagnitude as a variable to store area you probalbly what to have

Amagnitude = PRF_GRSUM1(Amagnitude);

and it should be the first line of global summation block, as you are going to use area value below

[UchihaMadara]

Tried this and thank you for your suggestion.
But no joy!

[AlexanderZ]

show here fixed code and log