rough-embed.h

Embedded law of the wall
- Implementation
- References

Embedded law of the wall

The so-called “law of the wall” is a major experimental result on wall-bounded turbulent flows. It may be exploited to include the main effect of a rough surface in simulations without computing the turbulent transfer directly.

The variation of the “averaged” value of a velocity component u_i at a distance d from a wall, inside a wall-bounded turbulent flow may exibit the logaritmic law of a no-slip wall:

\displaystyle U(d) = \frac{{u}_{*}}{\kappa} \mathrm{ln}\left(\frac{d}{z_0}\right),

with \kappa the Von Karman constant and z_0 the roughness length (d \gg z_0). We can compute u_* from,

\displaystyle u_* = \frac{\kappa}{\mathrm{ln}\left(\frac{d}{z_0}\right)}U(d).

This relates to the flux (F_\phi) of a dummy variable \phi:

\displaystyle F_\phi = \frac{u_* \kappa}{\mathrm{ln}(\frac{d}{z_{0,\phi}})} \left( \phi (d) - \phi(0) \right),

where z_{0,\phi} is an \phi-specific roughness length (by default, it is equal to z_0. These scalars must be added to a list called rough_list, their surface value and roughness length are defined as scalar-field attributes. *

scalar * rough_list = NULL;

attribute {
  double rl;
  double s_val;
}

#ifndef Z_0 //z0 for momentum -> default for scalars 
#define Z_0 (L0/(10*(1 << grid->maxdepth)))
#endif

event defaults (i = 0) {
  foreach_dimension() {
    scalar s = u.x;
    rough_list = list_add (rough_list, u.x);
    s.rl = Z_0;
    s.s_val = 0; //no-slip;
  }
}

// A user interface for adding scalar fields
struct rough_s {
  scalar s;
  double s_val;
  double z0;
};

void add_rough_scalar (struct rough_s p) {
  scalar s = p.s;
  rough_list = list_add (rough_list, s);
  s.rl = p.z0 ? p.z0 : Z_0;
  s.s_val = p.s_val; 
}

The code in this header file implements the closure, and applies the resulting tendencies. It is supposed to be used with embedded boundaries and the centered solver.

Implementation

#if (dimension == 2)
#define VOL (sq(Delta))
#elif (dimension == 3)
#define VOL (cube(Delta))
#endif
double karman = 0.3679; //= exp(-1)

Next, the “work horse”: a function that computes F_\phi for all fields in rough_list. Much of the used method is not very different from using dirichlet() conditions with embed.h and we follow the wisdom presented there. Note that the interpolator mixes the finite-volume representation of the scalar-field values with the vaulues at cell-centered locations. The written implementation considers a surface whose normal component is mostly in the x direction. The function returns u_*.

foreach_dimension() { 
  static inline double u_tau_x (Point point, scalar cs, coord n, coord p,
				scalar * list, double * flx) {// List and fluxes 
    foreach_dimension()
      n.x = - n.x;     // inward not outward
    double d[2] = {nodata, nodata}, v[2][list_len(list)], u_star = 0; //Distances and values;
    bool defined = true;
    foreach_dimension()
      if (defined && !fs.x[(n.x > 0.)])
	defined = false;
    if (defined) {
      for (int l = 0; l < 2; l++) {
	int i = (l + 1)*sign(n.x);
	d[l] = (i - p.x)/n.x;
	double y1 = p.y + d[l]*n.y;
	int j = y1 > 0.5 ? 1 : y1 < -0.5 ? -1 : 0;
	y1 -= j;
#if dimension == 2
	if (fs.x[i + (i < 0),j] && fs.y[i,j] && fs.y[i,j+1] &&
	    cs[i,j-1] && cs[i,j] && cs[i,j+1]) {
	  int c = 0;
	  for (scalar s in list) 
	    v[l][c++] = quadratic (y1, (s[i,j-1]), (s[i,j]), (s[i,j+1]));
	}
#else // dimension == 3
	double z = p.z + d[l]*n.z;
	int k = z > 0.5 ? 1 : z < -0.5 ? -1 : 0;
	z -= k;
	bool defined = fs.x[i + (i < 0),j,k];
	for (int m = -1; m <= 1 && defined; m++)
	  if (!fs.y[i,j,k+m] || !fs.y[i,j+1,k+m] ||
	      !fs.z[i,j+m,k] || !fs.z[i,j+m,k+1] ||
	      !cs[i,j+m,k-1] || !cs[i,j+m,k] || !cs[i,j+m,k+1])
	    defined = false;
	if (defined) {
	  // bi-quadratic interpolation
	  int c = 0;
	  for (scalar s in list) {
	    v[l][c++] =
	      quadratic (z,
			 quadratic (y1,
				    (s[i,j-1,k-1]), (s[i,j,k-1]), (s[i,j+1,k-1])),
			 quadratic (y1,
				    (s[i,j-1,k]),   (s[i,j,k]),   (s[i,j+1,k])),
			 quadratic (y1,
				    (s[i,j-1,k+1]), (s[i,j,k+1]), (s[i,j+1,k+1])));
	  }
	}
#endif // dimension == 3
	else {
	  d[l] = nodata; 
	  break;
	}
      }
      if (d[0] != nodata && d[1] != nodata) { // Two points
	double U[2] = {0., 0.};
	for (int l = 0; l < 2; l++) {
	  int dim = 0;
	  foreach_dimension() {
	    U[l] += sq(v[l][dim]);
	    dim++;
	  }
	}
	for (int l = 0; l < 2; l++) 
	  U[l] = U[l] > 0 ? sqrt(U[l]) : 0;
	u_star = karman*(U[0]/(log(Delta*d[0]/Z_0)) +
			 U[1]/(log(Delta*d[1]/Z_0)))/2.; //mean
	int c = 0;
	for (scalar s in list) {
	  flx[c] = -u_star*((v[0][c] - s.s_val)/(log(Delta*d[0]/s.rl)) +
			    (v[1][c] - s.s_val)/(log(Delta*d[1]/s.rl)))/2.; //mean
	  c++;
	}
      }
      else if (d[0] != nodata || d[1] != nodata) { //Single point interpolation
	int l;
	if (d[0] != nodata)
	  l = 0;
	else
	  l = 1;
	double U = 0.;
	int dim = 0;
	foreach_dimension() {
	  U += sq(v[l][dim]);
	  dim++;
	}
	U = U > 0 ? sqrt(U) : 0; 
	u_star = karman*U/(log(Delta*d[l]/Z_0));
	int c = 0;
	for (scalar s in list) {
	  flx[c] = -u_star*(v[l][c] - s.s_val)/(log(Delta*d[l]/s.rl));
	  c++;
	}
      }
    } else { //else Poor estimate based on the embedded cell itself :(
      d[0] = max (10*Z_0, fabs (p.x/n.x));
      double U = 0;
      foreach_dimension()
	U += sq(u.x[]);
      U = U > 0 ? sqrt(U) : 0.;
      u_star = karman*U/(log(Delta*d[0]/Z_0));
      int c = 0;
      for (scalar s in list) 
	flx[c++] = -u_star*(s[] - s.s_val)/(log(Delta*d[0]/s.rl));
    }
    return u_star;
  }
}

The function rough_wall() calles u_tau_i() depending on the dominant orientation of the boundary. The function returns u_*.

static inline double u_tau (Point point, scalar cs, coord n, coord p,
			    scalar * list, double * flx) {
#if dimension == 2
  foreach_dimension()
    if (fabs(n.x) >= fabs(n.y))
      return u_tau_x (point, cs, n, p, list, flx);
#else // dimension == 3
  if (fabs(n.x) >= fabs(n.y)) {
    if (fabs(n.x) >= fabs(n.z))
      return u_tau_x (point, cs, n, p, list, flx);
  }
  else if (fabs(n.y) >= fabs(n.z))
    return u_tau_y (point, cs, n, p, list, flx);
  return u_tau_z (point, cs, n, p, list, flx);
#endif // dimension == 3
  return nodata;
}

The drag term is included just before the diffusive step. The total change (m^{1+\mathrm{dimension}}s^{-1}= u_{*}) for each boundary fragment is stored in the flxdt vector.

event viscous_term (i++) {
  scalar * dsl = list_clone (rough_list); //change for each field
  foreach() { //Boundary cells.
    for (scalar ds in dsl)
      ds[] = 0.;
    double scs = 0;
    if (cs[] > 0. && cs[] < 1.) {
      foreach_neighbor(1)
	scs += cs[];
      double sf = 1/Delta; //Scale factor
      foreach_dimension()
	sf *= Delta;
      coord n = facet_normal (point, cs, fs), p;
      double alpha = plane_alpha (cs[], n);
      double area = sf*plane_area_center (n, alpha, &p);
      double flxs[list_len(rough_list)];
      u_tau (point, cs, n , p, rough_list, flxs);
      int c = 0;
      for (scalar ds in dsl)
	ds[] = dt*area*flxs[c++]/scs;
    }
  }
  boundary (dsl);

The drag is added to the cells near the boundary (smearing). For additional stability and physical consistency, the new velocity must be between the current value and zero…

Note that VOL = dv()/cs[].

  foreach() {
    scalar s, ds;
    for (s, ds in rough_list, dsl) {
      double se = 0;
      foreach_neighbor(1)
	se += ds[];
      if (cs[] > 0) {
	double new = s[] > s.s_val ? max(s[] + se/VOL, s.s_val) : min (s[] + se/VOL , s.s_val);
	s[] = s[] > s.s_val ? min (new, s[]) : max (new, s[]);
      }
    }
  }
  boundary (rough_list);
}

References

Juan C. Bergmann, 1998, A physical interpretation of von Karman’s constant based on asymptotic considerations—A new value. Journal of the Atmospheric Sciences. doi