myembed-moving.h

Advection of a discrete rigid body

Advection of a discrete rigid body

We detail here the first-order in time coupling algorithm we use to advect the discrete rigid body \Gamma_{\Delta}, of boundary \delta \Gamma_{\Delta}, compatible only with the centered Navier-Stokes solver. We assume that there is only one rigid body and consider two advection scenarios:

imposed embedded boundary motion;
explicit weak fluid-solid coupling.

Setup

In each scenario, the discrete rigid body \Gamma_{\Delta} is characterized by: the position of the center of mass \mathbf{x}_{\Gamma} stored in p_p, the translation velocity of the center of mass \mathbf{u}_{\Gamma} stored in p_u, the angular velocity about the center of mass \mathbf{\omega}_{\Gamma} stored in p_w and their corresponding accelerations stored in p_au, p_aw.

coord p_p;        // Position
coord p_u, p_w;   // Velocity
coord p_au, p_aw; // Acceleration

The remaining unknown is the location and shape of the discrete rigid boundary \delta \Gamma_{\Delta}, which we describe using a user-defined distance function \Phi given by the function p_shape(). Note that this function must contain a call to the fractions_cleanup() function.

extern void p_shape (scalar c, face vector f, coord p);

Fluid-solid coupling

The following no-slip Dirichlet boundary condition for velocity and Neumann boundary condition for pressure on the discrete rigid boundary \delta \Gamma_{\Delta} allow us to couple the motion of the fluid and the discrete rigid body \Gamma_{\Delta}:

\displaystyle \left\{ \begin{aligned} & \mathbf{u} = \mathbf{u}_{\Gamma} = \mathbf{u_{\Gamma}} + \omega_{\Gamma} \times \left(\mathbf{x} - \mathbf{x}_{\Gamma}\right) \\ & {\nabla}_{\Gamma} p = -\rho \frac{\mathrm{D} \mathbf{u}}{\mathrm{D}t} \cdot \mathbf{n}_{\Gamma} + \left(\nabla \cdot \left(\mu \nabla \mathbf{u} \right)\right) \cdot \mathbf{n}_{\Gamma} = -\rho \left( \frac{\mathrm{d} \mathbf{u_{\Gamma}}}{\mathrm{d}t} + \frac{\mathrm{d} \mathbf{\omega_{\Gamma}}}{\mathrm{d}t} \times \left(\mathbf{x} - \mathbf{x}_{\Gamma}\right) + \mathbf{\omega_{\Gamma}} \times \mathbf{\omega_{\Gamma}} \times \left(\mathbf{x} - \mathbf{x}_{\Gamma}\right) \right) \cdot \mathbf{n}_{\Gamma} + \left(\nabla \cdot \left(\mu \nabla \mathbf{u} \right)\right) \cdot \mathbf{n}_{\Gamma}. \end{aligned} \right.

In practice, we however use the following homogeneous Neumann boundary condition for pressure: \displaystyle {\nabla}_{\Gamma} p = 0 This condition is suitable for a fixed rigid body and we have found that using it with a moving rigid body does not significantly affect the computed solution.

No-slip Dirichlet boundary condition for velocity

The function velocity_noslip_x() computes the previously defined no-slip boundary condition for the x-component of the velocity \mathbf{u}.

foreach_dimension()
static inline double velocity_noslip_x (Point point,
					coord up, coord wp, coord pp,
					double xc, double yc, double zc)
{
  assert (cs[] > 0. && cs[] < 1.);

We first compute the relative position \mathbf{r} = \mathbf{x} - \mathbf{x}_{\Gamma}.

  // The coordinate x,y,z are not permuted with foreach_dimension()
  coord r = {xc, yc, zc};
  foreach_dimension() {
    r.x -= pp.x;
    if (Period.x) {
      if (fabs (r.x) > fabs (r.x + (L0)))
	r.x += (L0);
      if (fabs (r.x) > fabs (r.x - (L0)))
	r.x -= (L0);
    }
  }

We then compute the veolcity (translation + rotation).

#if dimension == 2
    coord sgn = {-1, 1};
    return (up.x) + sgn.x*wp.x*(r.y);
#else // dimension == 3
    return (up.x) + wp.y*(r.z) - wp.z*(r.y);
#endif // dimension
}

Neumann boundary condition for pressure

The function pressure_gradient() returns in da the contribution of all the components of the particle acceleration and viscous stresses to the pressure gradient.

static inline void pressure_acceleration (Point point,
					  coord dup, coord dwp,
					  coord wp, coord pp,
					  double xc, double yc, double zc,
					  coord * da)
{
  assert (cs[] > 0. && cs[] < 1.);

We first compute the relative position \mathbf{r} = \mathbf{x} - \mathbf{x}_{\Gamma}.

  // The coordinate x,y,z are not permuted with foreach_dimension()
  coord r = {xc,yc,zc};
  foreach_dimension() {
    r.x -= pp.x;
    if (Period.x) {
      if (fabs (r.x) > fabs (r.x + (L0)))
	r.x += (L0);
      if (fabs (r.x) > fabs (r.x - (L0)))
	r.x -= (L0);
    }
  }

We first compute the acceleration of the particle gu.

  coord gu;
  foreach_dimension()
    gu.x = 0.;
  
  // Rotational acceleration
  coord dw1, dw2;
#if dimension == 2
  coord sgn = {-1, 1};
  foreach_dimension() {
    dw1.x = sgn.x*dwp.x*(r.y);
    dw2.x = - sq(wp.x)*(r.x);
  } 
#else // dimension = 3
  double w1 = 0., w2 = 0.;
  foreach_dimension() {
    w1 += wp.x*(r.x);
    w2 += sq(wp.x);
  }
  foreach_dimension() {
    dw1.x = dwp.y*(r.z) - dwp.z*(r.y);
    dw2.x = w1*wp.x - w2*(r.x);
  }
#endif // dimension

We include the translational and rotational accelerations.

  foreach_dimension()
    gu.x = dup.x + dw1.x + dw2.x;

We then compute the viscous contribution.

  coord gmu;
  foreach_dimension()
    gmu.x = 0.;
  
  /* foreach_dimension() { */
  /*   scalar s = u.x; */
  /*   double a = 0.; */
  /*   foreach_dimension() */
  /* 	a += (mu.x[1]*face_gradient_x (s, 1) - mu.x[0]*face_gradient_x (s, 0))/Delta; */
  /*   double b, c = embed_flux (point, u.x, mu, &b); */
  /*   a += (b + c*u.x[]); */
  /*   gmu.x = -a; */
  /* } */

Finally, we combine both contributions.

  foreach_dimension()
    da->x = (gu.x + gmu.x); 
}

foreach_dimension()
static inline double pressure_acceleration_x (Point point,
					      coord dup, coord dwp,
					      coord wp, coord pp,
					      double xc, double yc, double zc)
{
  coord da;
  pressure_acceleration (point, dup, dwp, wp, pp, xc, yc, zc, &da);
  return da.x;
}

The following function finally computes the Neumann boundary condition for pressure, where the inward unit normal \mathbf{n}_{\Gamma} is pointing from fluid to solid. To switch to a homogenous Neumann boundary condition, the user can set the scalar attribute neumann_zero to true. The default is false.

attribute {
  bool neumann_zero;
}

static inline double pressure_neumann (Point point,
				       coord dup, coord dwp,
				       coord wp, coord pp,
				       double xc, double yc, double zc)
{
  coord da = {0., 0., 0.};
  pressure_acceleration (point, dup, dwp, wp, pp, xc, yc, zc, &da);

  coord b, n;
  embed_geometry (point, &b, &n);
  
  double dpdn = 0.;
  foreach_dimension()
    dpdn += (-da.x)*n.x;
  dpdn *= rho[]/(cs[] + SEPS);
  return dpdn;
}

Dump and restore a particle

typedef struct {
  coord c;      // Center of mass
  coord u, w;   // Velocity
  coord au, aw; // Acceleration
} particle;

struct p_Dump {
  char * file; // File name
  particle * list; // List of particles
  FILE * fp; // File pointer
  bool unbuffered;
};

void p_dump (struct p_Dump p)
{
  FILE * fp = p.fp;
  char def[] = "p_dump", * file = p.file ? p.file : p.fp ? NULL : def;

  char * name = NULL;
  if (file) {
    name = (char *) malloc (strlen(file) + 2);
    strcpy (name, file);
    if (!p.unbuffered)
      strcat (name, "~");
    if ((fp = fopen (name, "w")) == NULL) {
      perror (name);
      exit (1);
    }
  }
  assert (fp);

  // Get particle data (only 1 particle for now)
  int p_n = 1;
  particle * p_list = p.list;

  // Dump particle data
  fwrite(p_list, sizeof(*p_list), p_n, fp);
  
  /* free (p_list); */
  if (file) {
    fclose (fp);
    if (!p.unbuffered)
      rename (name, file);
    free (name);
  }
}

bool p_restore (struct p_Dump p)
{
  FILE * fp = p.fp;
  char * file = p.file;
  if (file && (fp = fopen (file, "r")) == NULL)
    return 0;
  assert (fp);

  // Read particle data (only 1 particle for now)
  int p_n = 1;
  particle * p_list = p.list;

  if (fread (p_list, sizeof(*p_list), (p_n), fp) < 1) {
    fprintf (ferr, "#p_restore(): error reading particle data\n");
    exit (1);
  }

  foreach_dimension() {
    p_p.x  = p_list->c.x;
    p_u.x  = p_list->u.x;
    p_w.x  = p_list->w.x;
    p_au.x = p_list->au.x;
    p_aw.x = p_list->aw.x;
  }
  
  /* free (p_list); */
  if (file)
    fclose (fp);
  
  return true;
}

Initialization

event defaults (i = 0)
{  
  foreach_dimension() {
    p_p.x  = 0.;
    p_u.x  = 0.;
    p_w.x  = 0.;
    p_au.x = 0.;
    p_aw.x = 0.;
  }
}

event init (i = 0)
{

We decrease the value of the CFL (same value as the one used with the VOF algorithm).

  CFL = 0.5;

  if (!restore (file = "restart")) { // No restart

We initialize the embedded boundary in the test case file. We also initialize the velocity in the test case file.

  }
  else { // Restart

We first restore the particle properties.

    particle pp_restore;
    bool p_restart = p_restore ("p_restart", &pp_restore);
    assert (p_restart == true);

We then need to initialize the face fraction fs since it is not dumped.

    p_shape (cs, fs, p_p);
  }

We then initialize the volume fraction at the previous timestep csm1. This needs to be done even when restarting the simulation as csm1 is not dumped.

  trash ({csm1});
  foreach()
    csm1[] = cs[];
  boundary    ({csm1});
  restriction ({csm1}); // Since restriction/prolongation depend on csm1

Finally, we define the boundary conditions for the velocity, the pressure gradient g and the presssure p on the embedded boundaries.

  p[embed]  = neumann (p.neumann_zero ?  0. : pressure_neumann (point, (p_au), (p_aw), (p_w), (p_p), x, y, z));
  pf[embed] = neumann (pf.neumann_zero ? 0. : pressure_neumann (point, (p_au), (p_aw), (p_w), (p_p), x, y, z));

#if dimension == 2
  u.n[embed]   = dirichlet (velocity_noslip_x (point, (p_u), (p_w), (p_p), x, y, z));
  u.t[embed]   = dirichlet (velocity_noslip_y (point, (p_u), (p_w), (p_p), x, y, z));
  uf.n[embed]  = dirichlet (velocity_noslip_x (point, (p_u), (p_w), (p_p), x, y, z));
  uf.t[embed]  = dirichlet (velocity_noslip_y (point, (p_u), (p_w), (p_p), x, y, z));

  g.n[embed] = dirichlet (p.neumann_zero ? 0. : pressure_acceleration_x (point, (p_au), (p_aw), (p_w), (p_p), x, y, z)); 
  g.t[embed] = dirichlet (p.neumann_zero ? 0. : pressure_acceleration_y (point, (p_au), (p_aw), (p_w), (p_p), x, y, z));
#else // dimension == 3
  u.n[embed]   = dirichlet (velocity_noslip_x (point, (p_u), (p_w), (p_p), x, y, z));
  u.t[embed]   = dirichlet (velocity_noslip_y (point, (p_u), (p_w), (p_p), x, y, z));
  u.r[embed]   = dirichlet (velocity_noslip_z (point, (p_u), (p_w), (p_p), x, y, z));
  uf.n[embed]  = dirichlet (velocity_noslip_x (point, (p_u), (p_w), (p_p), x, y, z));
  uf.t[embed]  = dirichlet (velocity_noslip_y (point, (p_u), (p_w), (p_p), x, y, z));
  uf.r[embed]  = dirichlet (velocity_noslip_z (point, (p_u), (p_w), (p_p), x, y, z));
  
  g.n[embed] = dirichlet (p.neumann_zero ? 0. : pressure_acceleration_x (point, (p_au), (p_aw), (p_w), (p_p), x, y, z)); 
  g.t[embed] = dirichlet (p.neumann_zero ? 0. : pressure_acceleration_y (point, (p_au), (p_aw), (p_w), (p_p), x, y, z));
  g.r[embed] = dirichlet (p.neumann_zero ? 0. : pressure_acceleration_z (point, (p_au), (p_aw), (p_w), (p_p), x, y, z));
#endif // dimension

As rho is used in the Neumann boundary condition for pressure, we need rho on all levels of the grid.

  restriction ({rho});

  boundary ({p, u, g, pf, uf});
}

Timestep

We modify the maximun timestep dtmax to account for the velocity of the discrete rigid body \Gamma_{\Delta}. Note that this event occurs before moving the embedded boundaries to their t^{n+1} position. We therefore use the position and boundary conditions at time t^n.

event stability (i++)
{
  foreach(reduction(min:dtmax)) {
    if (cs[] > 0. && cs[] < 1.) {

      // Barycenter and normal of the embedded fragment
      coord b, n;
      embed_geometry (point, &b, &n);

      // Local maximum velocity, in the direction of the normal
      double umax = 0.;
      foreach_dimension() {

We use here the boundary condition on the embedded boundary.

	bool dirichlet = true;
	double ub = (u.x.boundary[embed] (point, point,
					  u.x, &dirichlet));
	assert (dirichlet);	
	umax += (ub*n.x);
      }
      
      // Non-restrictive timestep (independent of *cs* and *fs*)
      double dte = Delta/(fabs (umax) + SEPS);
      if (dte < dtmax)
	dtmax = dte;
    }
  }
}

Prediction

event advection_term (i++)
{

In case of a periodic domain, we shift the coordinates of the center of mass p_p.

  coord p_o = {(X0), (Y0), (Z0)};
  foreach_dimension() {
    if (Period.x) {
      if (p_p.x < p_o.x)
	p_p.x += (L0);
      if (p_p.x > p_o.x + (L0))
	p_p.x -= (L0);
    }
  }

Step 1

We store the volume fraction defined at time t.

  trash ({csm1});
  foreach()
    csm1[] = cs[];
  boundary    ({csm1});
  restriction ({csm1}); // Since restriction/prolongation depend on csm1

Step 2 (prediction)

We advance the embedded boundary to time t+dt. This step requires the user to define the quantities p_p, p_u, p_w, p_au and p_aw at time t+dt.

  p_shape (cs, fs, p_p);

We then make sure not to use any values in newly emerged cells by setting the flag emerged to false and update all boundaries as the boundary conditions depend on the embedded boundaries.

  emerged = false;
  boundary (all);

We update the fluid properties to account for changes in the metric.

  event ("properties");

Step 3 (emerged cells)

Since the advection event is explict, we define here the values of the centered velocity u and centered pressure gradient g in emerged cells. We also define the pressures p and pf in emerged cells to provide an improved initial guess to the multigrid projection solver.

In the solid cells, we set all variables to 0. This is necessary to avoid mesh adaptation inside the solid boundaries. This might however lead to inaccuracies when using the default restriction_average operator.

  for (scalar s in {p, u, g, pf})
    foreach()
      if (cs[] <= 0.)
	s[] = 0.;

In the emerged cells, we use an extrapolation along the normal, using the function embed_extrapolate, to update the velocity u and the pressure gradient g (both use Dirichlet boundary conditions).

  for (scalar s in {u, g})
    foreach() {
      if (csm1[] <= 0. && cs[] > 0.) {

  	// Normal emerged cell
  	if (cs[] < 1.) {

  	  // Cell centroid, barycenter and normal of the embedded fragment
	  coord c, b, n;
	  embed_geometry (point, &b, &n);
	  double alpha = plane_alpha (cs[], n);
	  plane_center (n, alpha, cs[], &c); // Different from line_area_center
	
  	  // Dirichlet boundary condition on the embedded boundary
  	  bool dirichlet = true;
	  double sb = (s.boundary[embed] (point, point, s, &dirichlet));
	  assert (dirichlet);
	  
  	  // Emerged cell value
  	  s[] = embed_extrapolate (point, s, cs, n, c, sb);
  	}
	
  	// Pathological emerged cell (cs = 1)
  	else {
	  int navg = 0;
	  double savg = 0.;
	  foreach_neighbor(1)
	    if (cs[] > 0. && (emerged || csm1[] > 0.)) {
	      navg += 1;
	      savg += s[];
	    }
	  s[] = savg/(navg + SEPS);
  	}
      }
    }

As the pressure uses a Neumann boundary condition, we update the value of the pressures p and pf in emerged cells using their average over neighboring cells. Note that these values of pressure are used only as the initial condition for the Poisson solver.

  for (scalar s in {p, pf})
    foreach() {
      if (csm1[] <= 0. && cs[] > 0.) {
	int navg = 0;
	double savg = 0.;
	foreach_neighbor(1)
	  if (cs[] > 0. && (emerged || csm1[] > 0.)) {
	    navg += 1;
	    savg += s[];
	  }
	s[] = savg/(navg + SEPS);
      }
    }

Before using the boundary function, we set the emerged flag to true to indicate that all emerged cells have been updated and can now be used.

  emerged = true;
  boundary (all);

Step 4 (boundary conditions)

As rho is used in the Neumann boundary condition for pressure, we need rho on all levels of the grid.

  restriction ({rho});
}

TO DO

Multi-particles
Improve fluid-solid coupling algorithm (to remove added mass effect for density ratios close to 1)