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

Obtaining fluid properties using UDF

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

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   May 17, 2021, 07:31
Default Obtaining fluid properties using UDF
  #1
New Member
 
Join Date: Jan 2018
Posts: 29
Rep Power: 5
UchihaMadara is on a distinguished road
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.



Code:
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
Code:
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
UchihaMadara is offline   Reply With Quote

Old   May 17, 2021, 20:43
Default
  #2
Senior Member
 
Alexander
Join Date: Apr 2013
Posts: 1,723
Rep Power: 28
AlexanderZ will become famous soon enoughAlexanderZ will become famous soon enough
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
__________________
best regards


******************************
press LIKE if this message was helpful
AlexanderZ is offline   Reply With Quote

Old   May 18, 2021, 05:44
Default
  #3
New Member
 
Join Date: Jan 2018
Posts: 29
Rep Power: 5
UchihaMadara is on a distinguished road
Thank you for your response AlexenderZ,


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


Code:
#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.



Code:
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
UchihaMadara is offline   Reply With Quote

Old   May 19, 2021, 20:36
Default
  #4
Senior Member
 
Alexander
Join Date: Apr 2013
Posts: 1,723
Rep Power: 28
AlexanderZ will become famous soon enoughAlexanderZ will become famous soon enough
you are not making global summation for area
you've tried, but didn't use it:
was
Code:
A_glob = PRF_GRSUM1(Amagnitude);
as you are using Amagnitude as a variable to store area you probalbly what to have
Code:
Amagnitude = PRF_GRSUM1(Amagnitude);
and it should be the first line of global summation block, as you are going to use area value below
__________________
best regards


******************************
press LIKE if this message was helpful
AlexanderZ is offline   Reply With Quote

Old   May 24, 2021, 11:13
Default
  #5
New Member
 
Join Date: Jan 2018
Posts: 29
Rep Power: 5
UchihaMadara is on a distinguished road
Tried this and thank you for your suggestion.
But no joy!
UchihaMadara is offline   Reply With Quote

Old   May 24, 2021, 21:47
Default
  #6
Senior Member
 
Alexander
Join Date: Apr 2013
Posts: 1,723
Rep Power: 28
AlexanderZ will become famous soon enoughAlexanderZ will become famous soon enough
show here fixed code and log
__________________
best regards


******************************
press LIKE if this message was helpful
AlexanderZ is offline   Reply With Quote

Reply

Tags
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
Error in Two phase (condensation) modeling adilsyyed CFX 15 June 24, 2015 19:42
Sample UDF for Characterizing the Tube-Side 2-Phase Fluid thodij Fluent UDF and Scheme Programming 0 June 22, 2015 15:09
Changing the properties of a fluid TheForumLord FLUENT 2 June 12, 2014 01:25
energy eqn + constant fluid properties cfx_user Main CFD Forum 1 March 6, 2013 07:13
change particle properties using udf? ljp FLUENT 0 April 1, 2010 15:12


All times are GMT -4. The time now is 23:36.