sandbox/lopez/oscillation.c

Oscillating walls in a two-layer sandwich of fluid and solid.

This is the first test case reported in Sugiyama et al. 2011 who proposed a numerical scheme for fluid-structure coupled problems involving interaction of (incompressible) hyperlastic materials with (incompressible) fluids by modifying slightly available classic numerical schemes for incompressible fluid dynamics.

#include "navier-stokes/centered.h"
#include "src/hyperelastic.h"
#include "vof.h"

#define VW 1    //Amplitude of the wall velocity
#define OM M_PI //Frequency of the oscillation
#define C1 2.5  //Solid is a Neo-Hookean of G = 5
#define MU_F 1  //Fluid viscosity
#define LS 0.5  //Width of the solid layer
#define LF 0.5  //Width of the fluid layer 

face vector visc[];
scalar beta1v[];
scalar f[];
scalar * interfaces = {f};

u.t[top] = dirichlet (VW*sin(OM*t));

We will consider the upper half of the sandwich so the bottom is at rest.

u.t[bottom] = dirichlet (0.);

int MAXLEVEL = 6;

int main() {
  L0 = LS + LF ;
  origin (-L0/2, 0);
  periodic (right);
  DT = .005;
  mu = visc;
  beta1 = beta1v;
  stokes = true;
  init_grid (1 << MAXLEVEL);
  //  vector v = B.x;
  // v.n[bottom] = dirichlet(0.);
  run();
}

event init (i = 0) {
  vertex scalar psi[];

If a planar interface coincides with cell faces, i.e. the VOF scalar jumps from 0 to 1 in adjacent cells across the interfaces, fractions() complains. Displacing a little bit the interface avoid the problem.

  foreach_vertex ()
    psi[] = -y + (LS + 0.0001); 
  fractions (psi, f);

We follow the idea of Sugiyama el al (2011). In their work, it is defined \hat{\mathbf{B}} = \phi^q \mathbf{B} being \phi the solid volume fraction and q an elegible exponent. Note that \hat{\mathbf{B}} is, in contrast to \mathbf{B}, defined in all the computational domain. The equation governing its temporal evolution is still UCD(\hat{\mathbf{B}}) = 0.

  foreach() {
    u.x[] = 0.;
    foreach_dimension()
      B.x.x[] = sqrt(clamp(f[],0, 1.));
  }
}

event properties (i++) {
  foreach_face() {
    double T = (f[]+f[-1])/2.;
    visc.x[] = MU_F*(1-T);
  }

As in Sugiyama et al (2011) we will assume that the hyperelastic stresses in half-full cells will be proportional to the volume fraction of solid. So \beta_1 will be proportional to \sqrt{f}.

  foreach()
    beta1v[] = 2.*C1*sqrt(clamp(f[],0,1.));
}

In an imcompressible solid, the determinant of the Cauchy left deformation is an invariant, (det(\mathbf{B}))_t = 0.

event logfile (i++) {
  scalar invariant[];
  foreach () 
    invariant[] = B.x.x[]*B.y.y[]-B.y.x[]*B.x.y[];

  printf("t = %g sum = %g\n", t, statsf(invariant).sum);
}

Velocity profiles at some particular instants are saved.

event vprofile (t = {2, 2.8, 4, 4.8}) {
  char name[80];
  sprintf (name, "vprof-%g", t);
  FILE * fp = fopen (name, "w");
  foreach () 
    fprintf(fp,"%g %g \n", y, u.x[]);
  fclose(fp);
}

#if 0
#if QUADTREE
event gfsview (i += 4) {
  static FILE * fpg = popen("gfsview2D -s vista.gfv","w");
  output_gfs (fpg, t= t);
}
#endif
#endif

set xlabel 'y'
set ylabel 'v_x'
unset key
plot 'oscillation.sugiyama' u 1:2 w l, 'oscillation.sugiyama' u 1:3 w l, \
     'vprof-2' u 1:2, 'vprof-4' u 1:2,					\
     'vprof-2.8' u 1:2, 'vprof-4.8' u 1:2

Velocity field distribution at instant t=0 and t=0.8. (script)

References