[obscureed]

Hi Mel B,

Three phases are more difficult than two phases. So, it would be a good idea to start with two phases: build up the model, test the settings, get some results. You only need to extend it to include air if the air is actually relevant; on the other hand, if the air/water interface is horizontal, and negligible amounts of air are entrained into the mixture, then a flat zero-shear surface on top of a two-phase simulation is likely to be good enough.

You mention packing limits etc, which makes me think that you are (or should be) using the granular Eulerian model for the solid phase. This is a whole different topic, not requiring UDFs but requiring some tinkering and experimentation. The default settings are not very useful, so some research into good settings will be beneficial.

Good luck again!
Ed

[Mel B]

thank you for your answers
i will take your advice

[sd_athena]

Hi there,
I would like to do exactly the same thing as DF15:

Originally Posted by DF15

In my results from a simulation I have an unphysical layer of air in the first cell layer on a wall which I know should be completely wet, i.e. the volume fraction of air should be zero. I would therefore like to create a UDF which loops through all cells adjacent to this wall, and if the volume fraction of air is below a certain limit I want to set the volume fraction to 0.

But before changing the volume fraction of bottom layer, I want to first output the bottom volume fraction to see what it's like, but I ran into the same issue as DF15:

Originally Posted by DF15

I see that the water volume fractions are not 1-volume_fraction air, and they are sometimes negative....

To start, I'm running a serial simulation for a jet impinging on a moving wall problem. I have modified DF15's original code to output cell coordinates along with the volume fraction, and deleted the parallel processing bit. Here is my UDF:

#include "udf.h"

DEFINE_ADJUST(my_remove_air, domain)
{
  /* Create the thread pointers for the wall */
  Thread *t;              /* Mixture level thread */
  Thread *t_air;          /* Phase level thread */
  Thread *t_water;        /* Phase level thread */
  face_t f;               /* Face thread */
  cell_t c;               /* Cell thread */
  domain = Get_Domain(1); /* Get the domain */

  real volume_fraction_air = 1.0;
  real volume_fraction_water = 1.0;
  real x[ND_ND];

  /* Find mixture level thread and phase level threads */
  t = Lookup_Thread(domain, 15);     /* Surface thread ID from Fluent */
  t_air = THREAD_SUB_THREAD(t, 0);   /* Primary phase */
  t_water = THREAD_SUB_THREAD(t, 1); /* Secondary phase */

  /* Loop through all faces on surface */
  begin_f_loop(f, t)
  {
    /* Save coordinate to "x" */
    F_CENTROID(x, f, t);
    /* Get thread for adjacent cell */
    c = F_C0(f, t);

    /* Get volume fractions of air and water */
    volume_fraction_air = C_VOF(c, t_air);
    volume_fraction_water = C_VOF(c, t_water);

    /* Update volume fractions */
    if (volume_fraction_air > 0.01 && volume_fraction_air < 0.3)
    {
      Message("AIR: Point( %f , %f):\n", x[0], x[1]);
      Message("Volume fraction (air, water):( %f, %f)\n\n", volume_fraction_air, volume_fraction_water);
    }

    if (volume_fraction_water > 0.7 && volume_fraction_water < 1.0)
    {
      Message("WATER: Point( %f , %f):\n", x[0], x[1]);
      Message("Volume fraction (air, water):( %f, %f)\n\n", volume_fraction_air, volume_fraction_water);
    }
  }
  end_f_loop(f, t)
}

At the bottom layer of cells, I'm expecting to see some liquid cells around the impingement point (should be within -0.0006<x<0.0006), and air cells elsewhere, but instead, I get <0.3 air fraction cells at wrong places, where I wouldn't see "liquid" if I visualize the fraction in FLUENT contour. Also, the air fraction and liquid fraction never add up to 1. And water fraction can sometimes be a huge number. Is my code wrong somewhere? Please shed some light on this!

(And further down the road: if I want to modify the volume fraction, is it correct just do C_VOF(c, t_air) = 0; ? This line doesn't look grammatically correct to me...)

Attached here is my console output:

WATER: Point( -0.024553 , 0.000000):
Volume fraction (air, water):( 1.225000, 0.808232)

WATER: Point( -0.024488 , 0.000000):
Volume fraction (air, water):( 1.225000, 0.808232)

WATER: Point( -0.022731 , 0.000000):
Volume fraction (air, water):( 0.000000, 0.891896)

WATER: Point( -0.021105 , 0.000000):
Volume fraction (air, water):( 0.000000, 0.700325)

AIR: Point( -0.016941 , 0.000000):
Volume fraction (air, water):( 0.018332, 0.000000)

AIR: Point( -0.015363 , 0.000000):
Volume fraction (air, water):( 0.050668, 0.000000)

AIR: Point( -0.015379 , 0.000000):
Volume fraction (air, water):( 0.083759, 0.000000)

AIR: Point( -0.015233 , 0.000000):
Volume fraction (air, water):( 0.042202, 0.000000)

AIR: Point( -0.015249 , 0.000000):
Volume fraction (air, water):( 0.037195, 0.000000)

AIR: Point( -0.015200 , 0.000000):
Volume fraction (air, water):( 0.255415, 8406014382993345560005773999945249098457649811319658221261105377105996984900204045989904783493320367167334726754715575354461540949770328919402126586252368113313832418337812290118622493609911685882041727903820973078597410572206080.000000)

AIR: Point( -0.015216 , 0.000000):
Volume fraction (air, water):( 0.258801, 0.000000)

AIR: Point( -0.013606 , 0.000000):
Volume fraction (air, water):( 0.046318, 0.000000)

AIR: Point( -0.013622 , 0.000000):
Volume fraction (air, water):( 0.036467, 0.000000)

AIR: Point( -0.013574 , 0.000000):
Volume fraction (air, water):( 0.248774, 0.000000)

AIR: Point( -0.013590 , 0.000000):
Volume fraction (air, water):( 0.255933, 1415904190692875361038252955304465630115559044936289027645842231490168993143107507857414795208457459272622434946884331084204803670362613215973251472386990428993215477566345280859382482624590118912.000000)

AIR: Point( -0.011150 , 0.000000):
Volume fraction (air, water):( 0.159145, 34395815283264944020442328314845921807297559782693342197075083912706231639984653042658805569755460645848384277899533588381756014738832315539673205335912877559294546561139212288.000000)

AIR: Point( -0.011117 , 0.000000):
Volume fraction (air, water):( 0.042907, 0.000000)

AIR: Point( -0.011085 , 0.000000):
Volume fraction (air, water):( 0.042907, 0.000000)

AIR: Point( -0.011052 , 0.000000):
Volume fraction (air, water):( 0.042907, 22825143118965125085500476286251066076191524716354110593553745103214464278772278459609359938378446835604477012485365668489508495658533348179786451996859930680765668031536108144393517769924292957154915121263924511497599702671229452288.000000)

AIR: Point( -0.011020 , 0.000000):
Volume fraction (air, water):( 0.042907, 0.000000)

AIR: Point( -0.010987 , 0.000000):
Volume fraction (air, water):( 0.021453, 0.000000)

AIR: Point( -0.010955 , 0.000000):
Volume fraction (air, water):( 0.021453, 0.000000)

AIR: Point( -0.010922 , 0.000000):
Volume fraction (air, water):( 0.021453, 0.000000)

AIR: Point( -0.008043 , 0.000000):
Volume fraction (air, water):( 0.100000, 0.000000)

AIR: Point( -0.005115 , 0.000000):
Volume fraction (air, water):( 0.083000, 0.000000)

AIR: Point( -0.004530 , 0.000000):
Volume fraction (air, water):( 0.200000, 0.000000)

AIR: Point( -0.001114 , 0.000000):
Volume fraction (air, water):( 0.050000, 0.000000)

AIR: Point( -0.000023 , 0.000000):
Volume fraction (air, water):( 0.100000, 0.000000)

AIR: Point( 0.000096 , 0.000000):
Volume fraction (air, water):( 0.050000, 0.000000)


 Reversed flow on 122 faces of pressure-outlet 13.

 Reversed flow on 76 faces of pressure-outlet 14.

 Reversed flow on 525 faces of pressure-outlet 16.

 Reversed flow on 142 faces of pressure-outlet 17.
!  128 solution is converged
   128  1.6489e-06  1.9287e-04  1.7494e-04  9.4805e-04  0:00:18   14
/display/object/display v-fraction

> 
Flow time = 0.00014s, time step = 14

Calculation complete.

Thank you!

[vinerm]

The code should actually but you are trying to debug two things at a time; the code and Fluent solution. This is making it more complicated. Here are the suggestions

1. Remove the line containing Get_Domain since domain is being provided by Fluent from its argument. The value is already there

2. Use THREAD_SUB_THREADS or an mp_thread_loop_f instead of general purpose loop.

3. Both of your logical conditions are almost same. So, you only need one of those because both volume fractions are connected.

And for debugging the code, test it on a case where you already know what the result should be or compare against cell values reported in Contour plots.

[sd_athena]

Originally Posted by vinerm

The code should actually but you are trying to debug two things at a time; the code and Fluent solution. This is making it more complicated. Here are the suggestions

1. Remove the line containing Get_Domain since domain is being provided by Fluent from its argument. The value is already there

2. Use THREAD_SUB_THREADS or an mp_thread_loop_f instead of general purpose loop.

3. Both of your logical conditions are almost same. So, you only need one of those because both volume fractions are connected.

And for debugging the code, test it on a case where you already know what the result should be or compare against cell values reported in Contour plots.

Hi Vinerm,
Thank you so much for your suggestions!
I have modified my code accordingly, but the output is still not correct.
I am now hooking the UDF to a case/data file that I already know the result. From the reported contour plot, I am expecting to see that on the bottom wall (y=0): water_fraction > 10% on (-0.003 < x < 0.033); water_fraction < 10% on (-0.033 < x < -0.003). However, the output from UDF is like this:

====Thread id 15 ======
WATER: Point( -0.024391 , 0.000000):
Volume fraction (air, water):( 0.000000, 0.293361)

WATER: Point( -0.024293 , 0.000000):
Volume fraction (air, water):( 0.000000, 0.808232)

WATER: Point( -0.024228 , 0.000000):
Volume fraction (air, water):( 0.000000, 0.808232)

WATER: Point( -0.024098 , 0.000000):
Volume fraction (air, water):( 0.000000, 0.404115)

WATER: Point( -0.022666 , 0.000000):
Volume fraction (air, water):( 0.000000, 0.127668)

WATER: Point( -0.021446 , 0.000000):
Volume fraction (air, water):( 0.000000, 1.000000)

WATER: Point( -0.018144 , 0.000000):
Volume fraction (air, water):( 0.000000, 0.159145)

WATER: Point( -0.015168 , 0.000000):
Volume fraction (air, water):( 0.000000, 1.000000)

WATER: Point( -0.015184 , 0.000000):
Volume fraction (air, water):( 0.000000, 1.000000)

WATER: Point( -0.015070 , 0.000000):
Volume fraction (air, water):( 0.000000, 1.000000)

WATER: Point( -0.015086 , 0.000000):
Volume fraction (air, water):( 0.000000, 1.000000)

WATER: Point( 0.015197 , 0.000000):
Volume fraction (air, water):( 0.000000, 0.355090)

WATER: Point( 0.015135 , 0.000000):
Volume fraction (air, water):( 0.000000, 0.317810)

First off, air_fraction and water_fraction still don't add up to 1. (I checked, they add up to one for other face thread though, just not on the bottom wall thread.)
Second, the air_fraction is all zero.
Third, the locations of these mostly-water cells do not match with the reported locations from the Contour plots.

Could you please give me some advice on how to debug this?

FYI, my modified UDF is like this:

#include "udf.h"

DEFINE_ADJUST(new_remove_air, domain)
{
  cell_t c;
  face_t f;
  Thread *mixture_thread;
  Thread **pt = THREAD_SUB_THREADS(mixture_thread);
  int t_id = 1;

  real volume_fraction_air = 1.0;
  real volume_fraction_water = 1.0;
  real x[ND_ND];

  /* Loop through phase face threads */
  mp_thread_loop_f(mixture_thread, domain, pt)
  {
    t_id = THREAD_ID(mixture_thread);
    if (t_id == 15) /* Wall face has id of 15 */
    {
      Message("====Thread id %d ======\n", t_id);
      begin_f_loop(f, mixture_thread)
      {
        F_CENTROID(x, f, mixture_thread);
        /* Get thread for adjacent cell */
        c = F_C0(f, mixture_thread);

        /* Get volume fractions of air and water */
        volume_fraction_air = C_VOF(c, pt[0]);
        volume_fraction_water = C_VOF(c, pt[1]);

        /* Output volume fractions */
        if (volume_fraction_water > 0.1 && volume_fraction_water < 1.1)
        {
          Message("WATER: Point( %f , %f):\n", x[0], x[1]);
          Message("Volume fraction (air, water):( %f, %f)\n\n", volume_fraction_air, volume_fraction_water);
        }
      }
      end_f_loop(f, mixture_thread)
    }
  }
}

Thank you!

[sd_athena]

Originally Posted by MehdiAmin

I have this same problem! But I can't be able to fix it yet and follow your repliesand guides to anyone know how to solve this?
Thank you so much.
lectin

May I ask, is your problem with the wrong output of air/water fractions, or with non-physical air layers at the bottom wall where it should be completely wet?
(I'm having both problems and I'm hoping to solve them one by one. Just curious what your problem is )

[vinerm]

Though there are few things that are incorrect, such as, Thread **pt = THREAD_SUB_THREADS(mixture_thread), because at that stage mixture_thread does not have a useful value, it is being reset later. You just need to mention Thread **pt;

The bug though is with pt being fetched. mp_thread_loop_f provides pt, which is a point to face thread for phase i. In the code, this is being used for C_VOF. C_VOF requires pt pointing to cell thread. I am not very sure if these point to same but if you read the last two sections at the following link, it appears that they point to different threads. You have to fetch pt for cell thread for use in C_VOF. Try using this pt for something that is face value, such as, F_U for both phases and then compare against contours. If it matches for velocity, then it is confirmed that you need different pt. All you can try using a cell loop instead to confirm if C_VOF returns correct values when cell loop is being used.

https://www.afs.enea.it/project/nept...df/node105.htm

[sd_athena]

Originally Posted by vinerm

C_VOF requires pt pointing to cell thread...
All you can try using a cell loop instead to confirm if C_VOF returns correct values when cell loop is being used.

You are exactly right, Vinerm! Thank you so much for pointing these out.

I have modified my UDF to use cell loop and identify the bottom layers of cells by y-coordinate. Now the output volume fractions from UDF are consistent with FLUENT contour plots.

Also, I tried to modify the volume fractions to perfectly wet by:

C_VOF(c, pt[0]) = 0.;
C_VOF(c, pt[1]) = 1.;

I am surprised that, these lines actually worked, because I find it weird to assign a numerical value to the result of some function call... Maybe the UDF library somehow overloaded "=" for the C_VOF() function?
(I'm glad that my code is working, just a bit curious about the weird syntax.)

For anyone who want to do the same thing, here is my working UDF:

#include "udf.h"

DEFINE_ADJUST(remove_air, domain)
{
  cell_t c;
  Thread *cell_threads;
  Thread **pt; /* point to cell thread for phase i. */

  real volume_fraction_air, volume_fraction_water;
  real x[ND_ND];

  /* Loop through phase cell threads */
  mp_thread_loop_c(cell_threads, domain, pt)
  {
    begin_c_loop(c, cell_threads)
    {
      C_CENTROID(x, c, cell_threads);

      if (x[1] < 0.000026) /* Identify bottom two layers of cells */
      {
        /* Get volume fractions of air and water */
        volume_fraction_air = C_VOF(c, pt[0]);
        volume_fraction_water = C_VOF(c, pt[1]);

        /* Find cells with air-liquid mixtures */
        if (volume_fraction_water > 0.7)
        {
          /* Output coordinates and vf to check */
          Message("Volume fraction (air, water):( %f, %f)\n\n", volume_fraction_air, volume_fraction_water);
          Message("WATER: Point( %f , %f):\n", x[0], x[1]);

          /* Reset volume fraction to complete liquid */
          C_VOF(c, pt[0]) = 0.;
          C_VOF(c, pt[1]) = 1.;
        }
      }
    }
    end_c_loop(c, cell_threads)
  }
}

[vinerm]

Good to know that it works. There is nothing weird in the code. C_ and F_ based variables can be both lvalue and rvalue. You can use these to fetch or assign data.

[sd_athena]

Originally Posted by vinerm

Good to know that it works. There is nothing weird in the code. C_ and F_ based variables can be both lvalue and rvalue. You can use these to fetch or assign data.

Thank you for the explanation. I understand now!

[MohammedGhazi]

Hi everyone,

I want to write a UDF for the surface tension force(source term in the momentum equation of the VOF model ) as the attached file.
Can anyone help me with that?

[rhrh]

Originally Posted by sd_athena

Hi Vinerm,
Thank you so much for your suggestions!
I have modified my code accordingly, but the output is still not correct.
I am now hooking the UDF to a case/data file that I already know the result. From the reported contour plot, I am expecting to see that on the bottom wall (y=0): water_fraction > 10% on (-0.003 < x < 0.033); water_fraction < 10% on (-0.033 < x < -0.003). However, the output from UDF is like this:

====Thread id 15 ======
WATER: Point( -0.024391 , 0.000000):
Volume fraction (air, water):( 0.000000, 0.293361)

WATER: Point( -0.024293 , 0.000000):
Volume fraction (air, water):( 0.000000, 0.808232)

WATER: Point( -0.024228 , 0.000000):
Volume fraction (air, water):( 0.000000, 0.808232)

WATER: Point( -0.024098 , 0.000000):
Volume fraction (air, water):( 0.000000, 0.404115)

WATER: Point( -0.022666 , 0.000000):
Volume fraction (air, water):( 0.000000, 0.127668)

WATER: Point( -0.021446 , 0.000000):
Volume fraction (air, water):( 0.000000, 1.000000)

WATER: Point( -0.018144 , 0.000000):
Volume fraction (air, water):( 0.000000, 0.159145)

WATER: Point( -0.015168 , 0.000000):
Volume fraction (air, water):( 0.000000, 1.000000)

WATER: Point( -0.015184 , 0.000000):
Volume fraction (air, water):( 0.000000, 1.000000)

WATER: Point( -0.015070 , 0.000000):
Volume fraction (air, water):( 0.000000, 1.000000)

WATER: Point( -0.015086 , 0.000000):
Volume fraction (air, water):( 0.000000, 1.000000)

WATER: Point( 0.015197 , 0.000000):
Volume fraction (air, water):( 0.000000, 0.355090)

WATER: Point( 0.015135 , 0.000000):
Volume fraction (air, water):( 0.000000, 0.317810)

First off, air_fraction and water_fraction still don't add up to 1. (I checked, they add up to one for other face thread though, just not on the bottom wall thread.)
Second, the air_fraction is all zero.
Third, the locations of these mostly-water cells do not match with the reported locations from the Contour plots.

Could you please give me some advice on how to debug this?

FYI, my modified UDF is like this:

#include "udf.h"

DEFINE_ADJUST(new_remove_air, domain)
{
  cell_t c;
  face_t f;
  Thread *mixture_thread;
  Thread **pt = THREAD_SUB_THREADS(mixture_thread);
  int t_id = 1;

  real volume_fraction_air = 1.0;
  real volume_fraction_water = 1.0;
  real x[ND_ND];

  /* Loop through phase face threads */
  mp_thread_loop_f(mixture_thread, domain, pt)
  {
    t_id = THREAD_ID(mixture_thread);
    if (t_id == 15) /* Wall face has id of 15 */
    {
      Message("====Thread id %d ======\n", t_id);
      begin_f_loop(f, mixture_thread)
      {
        F_CENTROID(x, f, mixture_thread);
        /* Get thread for adjacent cell */
        c = F_C0(f, mixture_thread);

        /* Get volume fractions of air and water */
        volume_fraction_air = C_VOF(c, pt[0]);
        volume_fraction_water = C_VOF(c, pt[1]);

        /* Output volume fractions */
        if (volume_fraction_water > 0.1 && volume_fraction_water < 1.1)
        {
          Message("WATER: Point( %f , %f):\n", x[0], x[1]);
          Message("Volume fraction (air, water):( %f, %f)\n\n", volume_fraction_air, volume_fraction_water);
        }
      }
      end_f_loop(f, mixture_thread)
    }
  }
}

Thank you!

Hi, I am also doing eulerian model and want to set the volume fraction of secondary phase to zero. That code really works for me. Thanks!
May I ask you that does the code set the volume fraction to zero for each iteration or it just set the volume fraction to zero in the post processing. Thanks a lot!

[sd_athena]

Originally Posted by rhrh

Hi, I am also doing eulerian model and want to set the volume fraction of secondary phase to zero. That code really works for me. Thanks!
May I ask you that does the code set the volume fraction to zero for each iteration or it just set the volume fraction to zero in the post processing. Thanks a lot!

Hi rhrh,
I'm glad my code was helpful for you!
I believe my code sets the volume fraction to zero after each iteration, which is a feature of any `DEFINE_ADJUST` macro. If you think setting volume fraction after each iteration is too frequent, you may use something like this to let it only runs after 5 time steps, for example:

#define EXECUTE_TIME_STEPS 5

DEFINE_ADJUST(remove_air, domain)
{
  if (N_TIME % EXECUTE_TIME_STEPS == EXECUTE_TIME_STEPS)
  {
   /* Code for setting volume fractions */
  }
}

[rhrh]

Thanks for your help! I got it.