sandbox/Antoonvh/SGS.h

Apply the Vreman Sub-Grid Scale Model

This header file can be used to run with a sub-grid scale model. By default, the eddy viscosity is based on the Vreman Eddy viscosity model. This file should included after the Navier-Stokes/centered.h module is included but before grid adaptation is carried out. This module also applies diffusion to the list of tracers.

#include "tracer.h"
#include "diffusion.h"
#include "vreman.h"

We initialize global variables: K_h and K_m are the diffusivity for heat and momentum, respectively. Pr is the turbulent Prandlt number, defined as Pr=\frac{K_h-\nu}{K_m-\nu}, with \nu being the molecular viscosity (molvis). Csmag is the classical Smagorinsky constant and tracers is a list of diffusive flow tracers, e.g. the Buoyancy field.

face vector Km[],Kh[];
(const) face vector Pr;
scalar Evis[]; //Cell Centered diffusivity
double molvis;  
double Csmag;
scalar * tracers;

Since K_m\propto K_h \propto \Delta^2, proper care is required for evaluating corresponding diffusivities at the different levels of resolution boundaries.

static inline void Evisprol(Point point,scalar s){
  foreach_child()
    Evis[]=bilinear(point,Evis)/4.; 
}
static inline void Evisres(Point point,scalar s){
  double sum = 0.;
  foreach_child()
    sum += s[];
  s[] = sum/2.;
}

We set some default values and for the parameters that may be overwritten by the users’ init event.

event defaults(i=0){
  if (dimension!=3) //Allow to run, but give a warning
    fprintf(stdout,"Warning %dD grid. The used formulations only make sense for 3D turbulence simulations\n",dimension);
  mu=Kh;
  Pr=unityf;
  molvis=0.;
  Csmag=0.12;
  Evis.prolongation=Evisprol;
  Evis.restriction=Evisres;

On tree grids we do not care about the diffusivities on refined and/or coarsend cells and faces. These should be properly reevaluated before they apear in any computation

#if TREE  
  Evis.refine=no_restriction;
  Evis.coarsen=no_restriction;
  foreach_dimension(){
    Kh.x.refine=no_restriction;
    Km.x.coarsen=no_restriction;
  }
#endif
}

We calculate the centered eddyviscosity and translate it into the required face-field diffusivity.

event Eddyvis(i++){
  eddyviscosity(Csmag,u,molvis,Evis); 
  boundary({Evis});
  foreach_face(){
    Km.x[]=(Evis[]+Evis[-1])/2;
    Kh.x[]=(Pr.x[]*(Km.x[]-molvis))+molvis;
  }

In 3D, there are 4 finer faces per coarser face. So consistency with K_m,K_h\propto \Delta^2 is conviniently achieved by applying the Boundary_flux function for these face fields

  boundary_flux({Km,Kh}); 
}

For some user-friendliness, the SGS-mixing of the tracers is hard-coded into this header file.

event tracer_diffusion(i++){
  for (scalar s in tracers)
    diffusion(s,dt,Kh);
}