sandbox/Antoonvh/vof-tracer-particles.h

    Particles and interfaces

    Tracer-particle paths may cross a tracer interface (represented by volume fractions), leading to inconsistencies. A hack to circumvent this complicated problem is to modify estimated particle positions based on the position of the reconstructed interface. Such “corrections” could be small, and may have a huge inpact on the resulting particle paths.

    Using the code in this header file, the user can choose to restrict particles to a specific side of the interface, to stay on the interface or to be unbounded (regular tracer particles). This is specified by an argument passed to new_vof_tracer_particles(). For ten particles use:

    ...
    Particles Pin  = new_vof_tracer_particles (10, 1);  // f[] == 1 phase
    Particles Pout = new_vof_tracer_particles (10, 0);  // f[] == 0 phase
    Particles Pon  = new_vof_tracer_particles (10, 3);  // On the interface
    Particles Pun  = new_vof_tracer_particles (10, -1); // Unbounded
    ...

    These particles should then be positioned accordingly.

    Implementation

    The tracer-particles code is extended to also include an array (Pphase) that marks the phase of a particle.

    The reposition() function takes the vof field and a list of particles (referenced by Particles P) as the input. It returns the number of particles that it repositioned. Note that only the particles inside partial cells are repositioned.

    extern scalar f; // the vof field
    #define u uf     // Also used for VOF advection
    #include "tracer-particles.h"
    #undef u
    #include "fractions.h"
    
    int * Pphase = NULL;
    
    Particles new_vof_tracer_particles (long unsigned int n, int phase) {
      Particles P = new_tracer_particles (n);
      Pphase = realloc (Pphase, sizeof(int) * (P + 1));
      Pphase[P] = phase;
      return P;
    }
    
    event free_vof_tracers (t = end) {
      free (Pphase);
      Pphase = NULL;
    }
    
    int reposition (scalar f, Particles P) {
      int rep = 0;
      if  (Pphase[P] >= 0) { //Do something
        double mir = Pphase[P] < 2 ? 2. : 1.; //Mirror *or* on interface
        foreach_particle_in(P) {
          Point point = locate (x, y, z); 
          coord cc = {x,y,z}; //Cell centre
          if (f[] > 1e-6 && f[] < 1 - 1e-6) {
    	coord n = interface_normal (point, f); //pointing outwards.
    	double ff = Pphase[P] == 1 ? f[] : 1 - f[];
    	if (Pphase[P] != 1) {
    	  foreach_dimension()
    	    n.x *= -1;
    	}
    	normalize (&n);
    	double alpha = plane_alpha (ff, n);
    	double ALP = 0;
    	foreach_dimension()
    	  ALP += n.x*(p().x - cc.x);
    	ALP /= Delta;
    	if (ALP > alpha || mir != 2) { //Reposition particle
    	  foreach_dimension()
    	    p().x -= mir*(ALP - alpha)*n.x*Delta; 
    	  rep++;
    	}
          }
        }
      }
      return rep;
    }
    
    event defaults (i = 0) 
      P_RK2 = true;        //Match the VOF advection order of accuracy
    
    event tracer_particles_step1 (i += 2, last) {//After RK step 3 (i.e 2)
      foreach_P_in_list (tracer_particles) {
        
        particle_boundary (P);
        reposition (f, P);
      }
    }

    Todo

    • Allow unbounded particles along side bounded counterparts.
    • Bind particles to the interface itself

    Test

    Usage