sandbox/Antoonvh/law-of-the-wall.h

    Law of the Wall

    Here we implement a closure for turbulent flow over a rough surface according to the so-called “log-law of the wall”. The idea is that fully resolving a turbulent boundary layer near a rough surface within a larger domain may be unfeasable for very high Reynodlsnumber flows. This is expected to be the case when the grid resolution is much larger than the typical size the roughness elements. Rather ad hoc, the effect of a wall can be summarized as that is exerts a drag force (momentum flux) on the flow.

    We use the well established result for a profile of a horizontal velocity profile (u(y)) in a wall-bounded turbulent flow in the so-called log layer:

    \displaystyle u(y) = \frac{u_{\tau}}{k}\text{ln}\left( \frac{y}{y_0}\right),

    where \text{ln}\left( x \right) is the natural logarithm of a dummy variable x, u_\tau is the friction velocity, k = 0.4 the Von Karmann constant, y is the height above the roughness elements that are characterezed by a so-called ‘roughness length’ y_0, i.e the height at which the velocity vanishes.

    Using this, we write for the discretized horizontal velocity in the cell (U_n = \sqrt{u_x^2 + u_z^2}) adjacent to the rough and no-slip boundary,

    \displaystyle U_n = \frac{\int_{y_0}^{y_0+\Delta} u \text{d}y}{\Delta} = \frac{u_\tau}{k\Delta}\left[ y\left( \text{ln}\left( \frac{y}{y_0}\right) -1\right) \right]_{y_0}^{y_0+\Delta}.

    When we assume \Delta \gg y_0, we can rewrite this result as,

    \displaystyle u_\tau\approx \frac{kU_n}{\text{ln}\left( \frac{\Delta }{y_0} \right)-1 }.

    Using Pythagoras,

    \displaystyle u_{\tau,x}^2 + u_{\tau,z}^2 = u_\tau^2,

    we may write,

    \displaystyle u_{\tau,i}^2 = \left( \frac{u_i}{U_n} \right) ^2 u_\tau^2

    which closes the problem, as we have expressed the stress as a function of the known values and constants u_x, u_z, y_0, \Delta and k

    implementation

    We assume there is no heterogeneity for y_0 and in order to prevent an excessive number of calculations of the logaritmic function, we make use of a look-up table (lut) containing values for \left( \frac{k}{\text{ln} \left( \frac{\Delta}{y_0} \right) - 1} \right) ^ 2 that can be used for adaptive and non-adaptive simulations using upto 19 levels of refinement. Noting that it is the user’s responsibility not to apply the closure when the assumptions do not make sense. I.e. e.g. in the Roughness sublayer.

    #ifndef Y_0
#define Y_0  (0.01)
#endif
double lut[20];

event init (t = 0){
  double k = 0.4;
  lut[0] = 0.; //level > 0
  for (int m = 1; m <= 19; m++)
    lut[m] = sq(k/(log(L0/(((double)(1<<m))*Y_0)) - 1.));
}

    Now we use the midpoint method to advance in time as we take into account the surface drag at the bottom surface.

    #define dvdt(u) (sign(u)*lut[level]*sq(u)/Delta)
event law_of_the_wall(i++){
  double ui; //scratch for the mid-point-value estimate 
  foreach_boundary(bottom){
    ui = u.x[] - (dt*dvdt(u.x[])/2.); 
    u.x[] -= dt*dvdt(ui);
    ui = u.z[] - (dt*dvdt(u.z[])/2.);
    u.z[] -= dt*dvdt(ui);
  }
}

    Tests