path: root/src/boundary/wallFunctionBoundaryPostProcessors3D.hh

/*  This file is part of the OpenLB library
 *
 *  Copyright (C) 2018 Marc Haussmann
 *  E-mail contact: info@openlb.net
 *  The most recent release of OpenLB can be downloaded at
 *  <http://www.openlb.net/>
 *
 *  This program is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU General Public License
 *  as published by the Free Software Foundation; either version 2
 *  of the License, or (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public
 *  License along with this program; if not, write to the Free
 *  Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
 *  Boston, MA  02110-1301, USA.
*/

#ifndef WALLFUNCTION_BOUNDARY_POST_PROCESSORS_3D_HH
#define WALLFUNCTION_BOUNDARY_POST_PROCESSORS_3D_HH

#include "wallFunctionBoundaryPostProcessors3D.h"
#include "core/finiteDifference3D.h"
#include "core/blockLattice3D.h"
#include "dynamics/firstOrderLbHelpers.h"
#include "core/util.h"
#include "utilities/vectorHelpers.h"

namespace olb {

template <typename T, typename S>
Musker<T,S>::Musker(T nu, T y, T rho) : AnalyticalF1D<T,S>(1), _nu(nu), _y(y),_rho(rho)
{
  this->getName() = "Musker";
}

template <typename T, typename S>
bool Musker<T,S>::operator()(T output[], const S tau_w[])
{
  T y_plus = _y*sqrt(tau_w[0]/_rho)/_nu;

  T a = 5.424;
  T b = 0.119760479041916168;
  T c = 0.488023952095808383;
  T d = 0.434;
  T e = 3.50727901936264842;

  output[0] = sqrt(tau_w[0]/_rho)*(a*atan(b*y_plus - c) +
                                   d*log(pow(y_plus+10.6, 9.6)/pow(pow(y_plus, 2) - 8.15*y_plus + 86, 2)) - e);

  // Condition for the sub-viscous layer : TODO MARC H
  if (output[0] < 0) {
    output[0] = y_plus * sqrt(tau_w[0]/_rho);
  }

  return true;
}

template <typename T, typename S>
PowerLawProfile<T,S>::PowerLawProfile(T nu, T u2, T y2, T y1, T rho) : AnalyticalF1D<T,S>(2), _nu(nu), _u2(u2), _y2(y2), _y1(y1), _rho(rho)
{
  this->getName() = "PowerLawProfile";
}

template <typename T, typename S>
bool PowerLawProfile<T,S>::operator()(T output[], const S input[])
{
  T tau_w = 0.0246384 * _rho * pow(_nu, 0.25) * pow(_u2, 1.75) / pow(_y2, 0.25);
  T u_tau = sqrt(tau_w/_rho);
  T y_plus = _y1 * u_tau / _nu;

  if (y_plus > 30.0) {
    output[0] = tau_w;
    output[1] = u_tau * 8.3 * pow(y_plus, 1./7.);
  } else if (y_plus < 30.0  && y_plus > 5.) {
    output[0] = tau_w;
    output[1] = u_tau * (-2.81742 + 4.85723 * log(y_plus));
  } else {
    output[0] = 2.*_u2 * _rho * _nu / _y2;
    if (output[0] < 0.) {
      output[0]*=-1.;
    }
    u_tau = sqrt(output[0]/_rho);
    y_plus = _y1 * u_tau / _nu;
    output[1] = u_tau * y_plus;
  }
  return true;
}


////////  WallFunctionBoundaryProcessor3D ////////////////////////////////
template<typename T, typename DESCRIPTOR>
WallFunctionBoundaryProcessor3D<T,DESCRIPTOR>::WallFunctionBoundaryProcessor3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_,
                                                                            BlockGeometryStructure3D<T>& blockGeometryStructure,
                                                                            std::vector<int> discreteNormal, std::vector<int> missingIndices,
                                                                            UnitConverter<T, DESCRIPTOR> const& converter, wallFunctionParam<T> const& wallFunctionParam,
                                                                            IndicatorF3D<T>* geoIndicator)
  : x0(x0_), x1(x1_), y0(y0_), y1(y1_), z0(z0_), z1(z1_),
    _blockGeometryStructure(blockGeometryStructure),
    _discreteNormal(discreteNormal), _missingIndices(missingIndices),
    _converter(converter), _wallFunctionParam(wallFunctionParam)
{
  // Define Direction and orientatation
  discreteNormalX = _discreteNormal[0];
  discreteNormalY = _discreteNormal[1];
  discreteNormalZ = _discreteNormal[2];

  int Type_BC = discreteNormalX*discreteNormalX + discreteNormalY*discreteNormalY + discreteNormalZ*discreteNormalZ;
  T normal_norm = sqrt(Type_BC); // l2 norm : magnitude of the vector
  if (Type_BC == 1) { // Straight plane</