[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:

 /***********************************************************************
   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:

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:

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:

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:

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:

        /* 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))); */

}

[AlexanderZ]

as far as I understand you've almost done with the issue.
you may try to use following code for activity function

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]

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:

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:

Material *m = THREAD_MATERIAL(t0), *sp;

Should be

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):

    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.

[AlexanderZ]

you are right, you need

Material *m = TP_MATERIAL(t0);

as you are working with particles

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

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

[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.,

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!!!

[AlexanderZ]

from Ansys FLuent Customizatiion manual

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)

[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.

[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.