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

Need some support for UDF (my first time)

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

Like Tree4Likes
  • 1 Post By AlexanderZ
  • 3 Post By Tobi

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   June 30, 2021, 15:33
Default Need some support for UDF (my first time)
  #1
Super Moderator
 
Tobi's Avatar
 
Tobias Holzmann
Join Date: Oct 2010
Location: Tussenhausen
Posts: 2,677
Blog Entries: 6
Rep Power: 49
Tobi has a spectacular aura aboutTobi has a spectacular aura aboutTobi has a spectacular aura about
Send a message via ICQ to Tobi Send a message via Skype™ to Tobi
ey everybody, I was in contact with the FLUENT support today regarding the evaporation of H2O2/H2O droplets. As stated in the theory guide, the liquid droplets/phase are assumed to be ideal which means that the activity coefficient is set to one. However, based on the customer, we know that the equilibrium to an ideal assumption is far away and hence, we need to incorporate the activity coefficient in the DPM model. We can use an already defined UDF for manipulating the concentration / partial pressure of each phase. The makro is given here: https://ansyshelp.ansys.com/account/...p_equilib.html
So, I just want to discuss the part I need to modify:


Code:
 /***********************************************************************
   UDF for defining the vapor particle equilibrium
   for multicomponent particles
 ***********************************************************************/
 #include <udf.h>
 DEFINE_DPM_VP_EQUILIB(raoult_vpe,tp,Tp,cvap_surf,Z)
 {
    int is;
    real molwt[MAX_SPE_EQNS];
    Thread *t0 = TP_CELL_THREAD(tp);  /* cell thread of particle location */
    Material *gas_mix = THREAD_MATERIAL(t0); /* gas mixture material */
    Material *cond_mix = TP_MATERIAL(tp); /* particle mixture material */
    int nc = TP_N_COMPONENTS(tp);    /* number of particle components */
    real molwt_cond = 0.;  /* reciprocal molecular weight of the particle */
    for (is = 0; is < nc; is++)
      {
         int gas_index = TP_COMPONENT_INDEX_I(tp,is); /* index of vaporizing
            component in the gas phase */
         if (gas_index >= 0)
           {
            /* the molecular weight of particle material */
            molwt[gas_index] =
           MATERIAL_PROP(MIXTURE_COMPONENT(gas_mix,gas_index),PROP_mwi);
            molwt_cond += TP_COMPONENT_I(tp,is) / molwt[gas_index];
           }
        }
    /* prevent division by zero */
    molwt_cond = MAX(molwt_cond,DPM_SMALL);
    for (is = 0; is < nc; is++)  {
     /* gas species index of vaporization */
     int gas_index = TP_COMPONENT_INDEX_I(tp,is);
     if(gas_index >= 0)
       {
        /* condensed material */
        Material * cond_c = MIXTURE_COMPONENT(cond_mix, is);
        /* condensed component molefraction */
        real xi_cond = TP_COMPONENT_I(tp,is)/(molwt[gas_index]*molwt_cond);
        /* particle saturation pressure */
        real p_saturation = DPM_vapor_pressure(tp, cond_c, Tp);
        if (p_saturation < 0.0)
          p_saturation = 0.0;
        /* vapor pressure over the surface, this is the actual Raoult law */
        cvap_surf[is] = activity * xi_cond * p_saturation / UNIVERSAL_GAS_CONSTANT / Tp;
       }
    }
  /* compressibility for ideal gas */
  *Z = 1.0;
 }
The last line includes a new variable named activity which is not given in the original Makro (so it is the field I need to calculate to incorporate real-gas-effects). The activity variable is a field such as «p_saturation» (in each parcel different). Based on the species we look at the moment, I have to call two functions different functions. I need something like that:


Code:
if (species[i] == "h2o")
{
    activity = activityH2O(temperature, mole-fraction-h2o)
}
else if (species[i] == "h2o2")
{
    activity = activityH2O2(temperature, mole-fraction-h2o)
}
As I never programmed in fluent, it is a bit overwhelming right now (I am programming c++ mostly in QT and OpenFOAM). However, I have no clue how to access the fields of fluent and how to define the two different functions. So I would like to have some support from you. I guess I need something like this:
Code:
double activityH2O(const class temperature, const class mole-fraction)
{
    double temp = 1;
    temp = exp(232*temperature) * moleFraction * some other things
}
However, as I never did something similar in Fluent I have no glue. Mostly, people use pre-defined Makros so I don't even know the Makro name for defining a function such as:
Code:
DEFINE_FUNCTION(activityH2O, ....)
What I get out of the code is that
  • the Temperature of the parcel is Tp
  • the mole fraction of component I is given by xi_cond
So the function should look somehow as:
Code:
double activityH2O(double Tp, double xi_cond)
{
    double tmp;
    // Calc tmp based on the droplet temperature Tp and h2o molefraction xi_cond
    return tmp;
}
But still, I am not sure which class I need to for the objects Tp and xi_cond and if the definition of the function would be fine. I used "double" here in that example, but I am not sure if this is valid? I guess it is «real» as given in the UDF example above (probably it is «real»). However, I am not sure how to compare strings in fluent and how to get the actual name of the species we are looking at. The line
Material * cond_c = MIXTURE_COMPONENT(cond_mix, is);
Gives me the idea that I could use the cond_c pointer to call (maybe) the name of the species such das:
string nameOfSpecies = cond_c->name();
Finally the last question. How can I be sure that this UDF is used later on?


Thank you very much for your support.
Tobi


PS: After re-reading and thinking about the stuff I guess I am close to the guy:


Code:
        /* condensed material */
        Material * cond_c = MIXTURE_COMPONENT(cond_mix, is);
        /* condensed component molefraction */
        real xi_cond = TP_COMPONENT_I(tp,is)/(molwt[gas_index]*molwt_cond);
        /* particle saturation pressure */
        real p_saturation = DPM_vapor_pressure(tp, cond_c, Tp);
        if (p_saturation < 0.0)
          p_saturation = 0.0;

    /* Activity coefficient */
    real activity = activity(tp, Tp, xi_cond);

        /* vapor pressure over the surface, this is the actual Raoult law */
        cvap_surf[is] = activity * xi_cond * p_saturation / UNIVERSAL_GAS_CONSTANT / Tp;
       }
    }
  /* compressibility for ideal gas */
  *Z = 1.0;
 }


 real activity(Thread* tp, real Tp, real x)
{
    real RT = 1.986 * Tp;  
    real B0 = -752 + 0.97 * Tp
    real B1 = 85;
    real B2 = 13


    /* if H2O */  

     return exp(pow(1-x, 2)/RT * (B0 + B1 * (1 - 4*x) + B2 * (1 -2*x) * (1 - 6*x)));


    /* if H2O2 */

     /* return exp(pow(x, 2)/RT * (B0 + B1 * (3 - 4*x) + B2 * (1 - 2*x) * (5 - 6*x))); */

}
__________________
Keep foaming,
Tobias Holzmann
Tobi is offline   Reply With Quote

Old   June 30, 2021, 21:54
Default
  #2
Senior Member
 
Alexander
Join Date: Apr 2013
Posts: 2,249
Rep Power: 32
AlexanderZ will become famous soon enoughAlexanderZ will become famous soon enough
as far as I understand you've almost done with the issue.
you may try to use following code for activity function
Code:
real activity(Thread* t0, real Tp, real x)
{
    real RT = 1.986 * Tp;  
    real B0 = -752 + 0.97 * Tp
    real B1 = 85;
    real B2 = 13

	if (FLUID_THREAD_P(t0) && NNULLP(THREAD_STORAGE(t0, SV_UDM_I)))
        {
		Material *m = THREAD_MATERIAL(t0), *sp;
		int i;
		if (MATERIAL_TYPE(m) == MATERIAL_MIXTURE)
            mixture_species_loop (m,sp,i)
                {
                if ((0 == strcmp(MIXTURE_SPECIE_NAME(m,i),"h2o")) || (0 == strcmp(MIXTURE_SPECIE_NAME(m,i),"H2O")) )
                   {
                    return exp(pow(1-x, 2)/RT * (B0 + B1 * (1 - 4*x) + B2 * (1 -2*x) * (1 - 6*x)));
                   }
				else
					return exp(pow(x, 2)/RT * (B0 + B1 * (3 - 4*x) + B2 * (1 - 2*x) * (5 - 6*x)));
				}
		}
}
real is build in fluent variable type which is real->float, when you run fluent in single precision mode, and real->double in case of double precision
Tobi likes this.
__________________
best regards


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

Old   July 1, 2021, 04:42
Default
  #3
Super Moderator
 
Tobi's Avatar
 
Tobias Holzmann
Join Date: Oct 2010
Location: Tussenhausen
Posts: 2,677
Blog Entries: 6
Rep Power: 49
Tobi has a spectacular aura aboutTobi has a spectacular aura aboutTobi has a spectacular aura about
Send a message via ICQ to Tobi Send a message via Skype™ to Tobi
Hey Alexander,



thank you for your reply. You got everything right. However, I will understand the code you posted. So can you give me some insights into the following lines:


Code:
if (FLUID_THREAD_P(t0) && NNULLP(THREAD_STORAGE(t0, SV_UDM_I)))
So you check if the pointer is valid? By the way, I am checking the particle phase and not the fluid phase. Hence:
Code:
Material *m = THREAD_MATERIAL(t0), *sp;
Should be
Code:
Material *m = TP_MATERIAL(t0);
right? Or is it different? I am sorry, I have no clue regarding the functions and what they do (is there any reference list such as Doxygen? in which I can check the functions)

As both activity coefficients are based on the mole fraction of h2o in the droplet phase I would do it like that (not the best but okay-level):
Code:
    for (is = 0; is < nc; is++)  {
     /* gas species index of vaporization */
     int gas_index = TP_COMPONENT_INDEX_I(tp,is);
     if(gas_index >= 0)
       {
        /* condensed material */
        Material * cond_c = MIXTURE_COMPONENT(cond_mix, is);
        /* condensed component molefraction */
        real xi_cond = TP_COMPONENT_I(tp,is)/(molwt[gas_index]*molwt_cond);
        /* particle saturation pressure */
        real p_saturation = DPM_vapor_pressure(tp, cond_c, Tp);
        if (p_saturation < 0.0)
          p_saturation = 0.0;

    /* ================ ADDED ================== */
    /* Get h2o mole fraction */
    real h2oMole = 0;  
    int indexh2o = -1;
    int gasIndexh2o = -1;
    int ns = 0;
     Material *sp;

    /* Get the species index and gas index of h2o */
    mixture_species_loop(cond_c, sp, ns)
    {
        if ((0 == strcmp(MIXTURE_SPECIE_NAME(cond_c,ns),"h2o") || (0 == strcmp(MIXTURE_SPECIE_NAME(cond_c,ns),"H2O"))
        {
            indexh2o = ns;
            gas_index_h2o = TP_COMPONENT_INDEX_I(tp,ns);
            break;
        }
    }

    /* Check if index is larger than zero - otherwise error */
    if ((indexh2o >= 0) && (gas_index_h2o >= 0))
    {
          /* stop fluent --> show message but how ??? */
    }

    /* get mole fraction of h2o in droplet phase */
     h2oMole = TP_COMPONENT_I(tp,indexh2o)/(molwt[gas_index_h2o]*molwt_cond)

    /* Activity coefficient */
    real activity = activity(tp, Tp, h2oMole);

        /* vapor pressure over the surface, this is the actual Raoult law */
        cvap_surf[is] = activity * xi_cond * p_saturation / UNIVERSAL_GAS_CONSTANT / Tp;
       }
    }
  /* compressibility for ideal gas */
  *Z = 1.0;
 }


 real activity(Thread* tp, real Tp, real x)
{
    real RT = 1.986 * Tp;  
    real B0 = -752 + 0.97 * Tp
    real B1 = 85;
    real B2 = 13

    if (FLUID_THREAD_P(t0) && NNULLP(THREAD_STORAGE(t0, SV_UDM_I)))
    {
        Material *m = TP_MATERIAL(t0), *sp;
        int i;
        if (MATERIAL_TYPE(m) == MATERIAL_MIXTURE)
        mixture_species_loop (m,sp,i)
        {
            if ((0 == strcmp(MIXTURE_SPECIE_NAME(m,i),"h2o")) || (0 == strcmp(MIXTURE_SPECIE_NAME(m,i),"H2O")) )
            {
                return exp(pow(1-x, 2)/RT * (B0 + B1 * (1 - 4*x) + B2 * (1 -2*x) * (1 - 6*x)));
            }
             else if (0 == strcmp(MIXTURE_SPECIE_NAME(m,i),"h2o2")) || (0 == strcmp(MIXTURE_SPECIE_NAME(m,i),"H2O2")) )
            {

                 return exp(pow(x, 2)/RT * (B0 + B1 * (3 - 4*x) + B2 * (1 - 2*x) * (5 - 6*x)));
            }
            else
            {
                   // Stop simulation and show some info that it can only be used for h2o/h2o2 parcels - but how (similar to the indexes)
            }
        }
    }
}
Your feedback would be highly appreciated.
Thank you in advance.
__________________
Keep foaming,
Tobias Holzmann
Tobi is offline   Reply With Quote

Old   July 1, 2021, 21:52
Default
  #4
Senior Member
 
Alexander
Join Date: Apr 2013
Posts: 2,249
Rep Power: 32
AlexanderZ will become famous soon enoughAlexanderZ will become famous soon enough
you are right, you need
Code:
Material *m = TP_MATERIAL(t0);
as you are working with particles


Code:
if (FLUID_THREAD_P(t0) && NNULLP(THREAD_STORAGE(t0, SV_UDM_I)))
as far as I understand, you are right here too, we are checking if pointer is valid, if memory was allocated.

Unfortunately, I don't know about the list of available UDF function similar to Doxygen. Some documentations could be found in Ansys Fluent Customization manual.

To stop the fluent you may use macro
Code:
Error("text");
Which stops fluent iterations and return the message, which is defined.

But I'm not sure if you can rerun simulation after this error, which is not convenient
__________________
best regards


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

Old   July 3, 2021, 08:34
Default
  #5
Super Moderator
 
Tobi's Avatar
 
Tobias Holzmann
Join Date: Oct 2010
Location: Tussenhausen
Posts: 2,677
Blog Entries: 6
Rep Power: 49
Tobi has a spectacular aura aboutTobi has a spectacular aura aboutTobi has a spectacular aura about
Send a message via ICQ to Tobi Send a message via Skype™ to Tobi
Hi, thank you for your hint regarding the ERROR function. Do you know if I can throw a message to the console somehow. E.g.,

Code:
msg("Calculate activity");
Just to get the information that this UDF is really used? No idea how to tell Fluent to use the UDF. The only possible way is to add the UDF to the Fluent Case and Compile it. But not sure if Fluent knows to use this function directly or if one has to activate it.


EDIT: Regarding the message, I guess that I will try this one printing values from udf!!!
__________________
Keep foaming,
Tobias Holzmann
Tobi is offline   Reply With Quote

Old   July 8, 2021, 20:21
Default
  #6
Senior Member
 
Alexander
Join Date: Apr 2013
Posts: 2,249
Rep Power: 32
AlexanderZ will become famous soon enoughAlexanderZ will become famous soon enough
from Ansys FLuent Customizatiion manual
Code:
After the UDF that you have defined using DEFINE_DPM_VP_EQUILIB is interpreted (Interpreting UDFs (p. 313)) or compiled (Compiling UDFs (p. 319)), the name of the argument that you supplied as the first DEFINE macro argument will become visible and selectable in the Create/Edit Materials dialog box in ANSYS Fluent. Note that before you hook the UDF, you’ll need to create particle injections in the Injections dialog box with the type Multicomponent chosen.
So you need to hook UDF function to execute it through simulation.
As you've mentioned it is possible to write information to console using Message("text"); macro (works same as printf in C)
__________________
best regards


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

Old   July 9, 2021, 14:57
Default
  #7
Super Moderator
 
Tobi's Avatar
 
Tobias Holzmann
Join Date: Oct 2010
Location: Tussenhausen
Posts: 2,677
Blog Entries: 6
Rep Power: 49
Tobi has a spectacular aura aboutTobi has a spectacular aura aboutTobi has a spectacular aura about
Send a message via ICQ to Tobi Send a message via Skype™ to Tobi
Hey Alex,

thank you very much for your support. The library is ready and works as expected. The last week I was using the manuals of fluent intensively and hence, I could solve a lot of problems. In addition, your feedback is highly valuable and hence, I want to thank you for taking your private time for your replys.
__________________
Keep foaming,
Tobias Holzmann
Tobi is offline   Reply With Quote

Old   August 2, 2021, 07:14
Default
  #8
Super Moderator
 
Tobi's Avatar
 
Tobias Holzmann
Join Date: Oct 2010
Location: Tussenhausen
Posts: 2,677
Blog Entries: 6
Rep Power: 49
Tobi has a spectacular aura aboutTobi has a spectacular aura aboutTobi has a spectacular aura about
Send a message via ICQ to Tobi Send a message via Skype™ to Tobi
Hey everybody,

I just want to give all of you (who are doing similar things) the procedure you need to make a correct H2O/H2O2 evaporation analysis:
  • You need to create a new UDF that handles the Equilibrium calculation (as given above)
  • Add the activity factor as described above
  • Modify the saturation pressure for H2O and H2O2 according to the air concentration. Here you must do the following:
    • Get cell in which the parcel is located
    • Calculate the mole fraction of H2O and H2O2 inside the cell
    • Get the cell pressure
    • Based on the three data, you need to calculate the real saturation pressure of H2O and H2O2 (iteratively)
  • Use the real saturation pressure for the concentration equilibrium

Done. The main challange was the iterative process of the H2O and H2O2 saturation pressure as these guys are not a function of temperature only as they are a function of the other species as well.

Doing so, you will have a very accurate analysis.
Some speedup possibilities:
  • One can use a function that calculates the real saturation pressures based on the cell temperature and save it to a UDMI
  • Use the UDMI stored real saturation pressures for the calculation inside the equilibrium function

This will save a lot of computational power as the vapor equilibrium calculation is called several times. However, the simplification lead to less accuracy as the saturation pressures are (a) not based on the parcel temperature (b) not updated during the DPM iteration.
AlexanderZ, shauwai and Hello_Wang like this.
__________________
Keep foaming,
Tobias Holzmann
Tobi is offline   Reply With Quote

Reply

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
[solidMechanics] Support thread for "Solid Mechanics Solvers added to OpenFOAM Extend" bigphil OpenFOAM CC Toolkits for Fluid-Structure Interaction 686 December 22, 2022 10:10
Elastic or Plastic Deformations in OpenFOAM NickolasPl OpenFOAM Pre-Processing 31 November 22, 2021 10:26
Extrusion with OpenFoam problem No. Iterations 0 Lord Kelvin OpenFOAM Running, Solving & CFD 8 March 28, 2016 12:08
simpleFoam error - "Floating point exception" mbcx4jc2 OpenFOAM Running, Solving & CFD 12 August 4, 2015 03:20
dynamic Mesh is faster than MRF???? sharonyue OpenFOAM Running, Solving & CFD 14 August 26, 2013 08:47


All times are GMT -4. The time now is 07:55.