marginal_embed_2.c

This code demonstrates a bug where the no-slip condition is violated at the embedded boundary in pressure-driven viscous microchannel flow. It is similar to Antoon van Hooft’s code for scalar mixing in a tube flow where the key difference is here I numerically compute the flow field instead of prescribing the exact solution.

This plot shows the spatial distribution of the solute at the final simulated time step.

This plot shows the spatial distribution of the x-component of the velocity at the final simulated time step. I expect a fully developed parabolic flow with zero slip at the embedded boundaries.

This plot shows the spatial distribution of the pressure field at the final simulated time step.

This plot shows the spatial distribution of the y-component of the velocity at the final simulated time step.

#include "embed.h"
#include "navier-stokes/centered.h"
#include "diffusion.h" 
#include "run.h"

#define tfinal 1
#define time_step 0.1
#define length 5
#define aspect_ratio 10

scalar c[], * tracers = {c};

int main() {

The goal is to simulate viscous Stokes flow, although the convective term of the Navier-Stokes equations can be turned on (leaving the default stokes = false) for debugging purposes.

  # if 1
  stokes = true;
  //TOLERANCE = HUGE;
  NITERMIN = 2;
  DT = time_step;
  # endif

The suggested workaround for this already-documented bug to adjust the domain size or origin did not do the trick.

  L0 = length; origin (0, 0);
  run();
}

event init (t = 0) {

You can define microchannel geometry here. As suggested in a previous post, moving the boundary (try adjusting the multiplying factor in yoffset) did not help.

  double width = L0/aspect_ratio; double yoffset = length*0.5;
  solid (cs, fs, -(width/2 - y + yoffset)*(-width/2 - y + yoffset));
  refine (cs[] && level < 10);
  fractions_cleanup(cs, fs);

  // boundary conditions
  c[left] = dirichlet(0); c[right] = dirichlet(0); c[embed] = neumann(0);

Both BC’s #1 and #2 led to some instabilities near the channel inlet on the left, and have very different flow fields despite similar pressure fields.

  // BC #1: assign a pressure gradient
  //p[left] = dirichlet(1)*cs[]; p[right] = dirichlet(0);
  
  // BC #2: uniform inflow and outflow
  //u.n[left] = dirichlet(1)*cs[]; u.n[right] = dirichlet(1)*cs[];

The no-slip condition is implemented here but seems to be ignored by the solver.

  // no-slip at embedded boundaries
  u.n[embed] = dirichlet(0); u.t[embed] = dirichlet(0);

Alternatively to either BC #1 or #2, you can drive flow through this acceleration term (making sure to select BC #1). In this case, the no-slip condition is observed, however if the advective term of the solute is turned on, the concentration c goes to zero everywhere.

  # if 1
  // using acceleration term to assign dpdx
  const face vector g[] = {-1.,0.};
  a = g;
  mu = fm;
  # endif

Initialize the solute as a thin plug.

  // intialize solute cloud
  double startpoint = length/4;
  foreach() 
    c[] = (startpoint < x)*(x < startpoint*1.1)*cs[];
}

The advection and diffusion of the scalar are enforced here, although while debugging the flow field I comment out the advection call.

// simple advection diffusion of passive scalar c
event advection_and_diffusion (i++) {
  face vector kappa[];
  foreach_face()
    kappa.x[] = fs.x[]*0.01;
  diffusion (c, dt, kappa);
  //advection({c}, u, DT);
}

// Generate images at the final timestep.
event stop (t = tfinal) {
  output_ppm (c, file = "c.png", n=400);
  output_ppm (u.x, file = "ux.png", n=400);
  output_ppm (u.y, file = "uy.png", n=400);
  output_ppm (p, file = "p.png", n=400);
}