Law of the wall consistent with SA turbulence model

git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@2476 e10b68d5-8a6e-419e-a041-bce267b0401d

Apps/Common/CFDenums.h

@@ -27,7 +27,9 @@ namespace CFD //! CFD scope
     BOUSSINESQ = 8,
     ALE        = 16,
     SUPG       = 32,
-    VMS        = 64
+    VMS        = 64,
+    SPALDING   = 128,
+    SA         = 256
   };
 }
 

Apps/Common/CMakeLists.txt

@@ -39,6 +39,7 @@ INCLUDE(cmake/UseMultiArch.cmake)
 # Common Navier-Stokes application sources
 ADD_LIBRARY(CommonIFEM BDF.C
                        Spalding.C
+		       SAWallLaw.C
                        StabilizationUtils.C
                        TimeIntUtils.C)
 TARGET_LINK_LIBRARIES(CommonIFEM ${IFEM_LIBRARIES})

Apps/Common/SAWallLaw.C

@@ -0,0 +1,99 @@
+//==============================================================================
+//!
+//! \file SAWallLaw.C
+//!
+//! \date Jun 19 2013
+//!
+//! \author Runar Holdahl / SINTEF
+//!
+//! \brief Implementation of Spalding parametrization of a turbulent
+//!        boundary layer. Computes the mean velocity uplus parallel
+//!        to a solid wall given the distance y and the tangential
+//!        velocity component ut.
+//!
+//==============================================================================
+
+#include "SAWallLaw.h"
+
+SAWallLaw::f_SA_wall_law::f_SA_wall_law() : B(5.03339087905055799), a1(8.148221580024245),
+					    a2(-6.9287093849022945), b1(7.4600876082527945), b2(7.468145790401841),
+					    c1(2.5496773539754747), c2(1.3301651588535228), c3(3.599459109332379),
+					    c4(3.6397531868684494)
+{}
+
+
+bool SAWallLaw::computeYplus(double y, double nu, double utan, double& yplus) const
+{
+  // Initial guess
+  yplus = y;
+
+  // yplus is zero for zero tangential velocity
+  if (utan < 1.0e-12)
+    return true;
+  
+  // Residual
+  double uplus = y*utan/(nu*yplus); 
+  double r = uplus - SA.up(yplus);
+  
+  int it = 0;
+  double drdyplus;
+  while ((fabs(r) > rtol) && (it < maxit)) {
+    drdyplus = -y*utan/(nu*yplus*yplus) - SA.dupdyp(yplus);
+    yplus -= r/drdyplus;
+    uplus = y*utan/(nu*yplus); 
+    r = uplus - SA.up(yplus);
+    it++;
+  }
+
+  return (it < maxit);
+}
+
+
+bool SAWallLaw::computeUstar(double y, double nu, double utan, double& ustar) const
+{
+  double yplus;
+  if (!this->computeYplus(y,nu,utan,yplus))
+    return false;
+
+  double uplus = SA.up(yplus);
+  ustar = utan/uplus;
+  return true;
+}
+
+
+bool SAWallLaw::computeTauB(double y, double nu, double utan, double& tauB) const
+{
+  // Constant
+  const double CbI = 4.0;
+
+  // Initial guess 
+  tauB = CbI*nu/y;
+
+  if (utan < 1.0e-12)
+    return true;
+
+  double yplus;
+  if (!this->computeYplus(y,nu,utan,yplus))
+    return false;
+
+  double uplus = SA.up(yplus);
+  double ustar = utan/uplus;
+  tauB = ustar*ustar/utan;
+
+  return true;
+}
+
+
+bool SAWallLaw::computeNu(double y, double nu, double utan, double& nuwall) const
+{  
+  double yplus;
+  if (!this->computeYplus(y,nu,utan,yplus))
+    return false;
+
+  //! Wall BC for Spalart-Allmaras
+  double uplus = SA.up(yplus);
+  double ustar = utan/uplus;;
+  nuwall = K*ustar*y;
+  
+  return true;
+}

Apps/Common/SAWallLaw.h

@@ -0,0 +1,106 @@
+//==============================================================================
+//!
+//! \file SAWallLaw.h
+//!
+//! \date Jun 19 2013
+//!
+//! \author Runar Holdahl / SINTEF
+//!
+//! \brief Implementation of Spalding parametrization of a turbulent
+//!        boundary layer. Computes the mean velocity uplus parallel
+//!        to a solid wall given the distance y and the tangential
+//!        velocity component ut.
+//!
+//==============================================================================
+
+#ifndef SA_WALL_LAW_H
+#define SA_WALL_LAW_H
+
+#include <cmath>
+#include <iostream>
+
+/*!
+  \brief Class representing Spalart-Allmaras parametrization of a turbulent
+         boundary layer
+*/
+
+class SAWallLaw 
+{
+ protected:
+  // Function defining Spalart-Allmaras parametrization
+  class f_SA_wall_law
+  {
+  protected:
+    //! Model constants
+    const double B;
+    const double a1;
+    const double a2;
+    const double b1;
+    const double b2;
+    const double c1;
+    const double c2;
+    const double c3;
+    const double c4;
+    
+  public:
+    //! \brief Constructor defining constants
+    f_SA_wall_law(); 
+    //! \brief Empty destructur
+    virtual ~f_SA_wall_law() {};
+    
+    //! \brief Spalart-Allmaras function definition
+    //! \param[in] yp y plus parameter
+    double up(double yp) const
+    {
+      double yppa1 = yp + a1;
+      double yppa2 = yp + a2;
+      
+      double up =  B + c1*log(pow(yppa1,2.0) + pow(b1,2.0));
+      up -= c2*log(pow(yppa2,2.0) + pow(b2,2.0));
+      up -= c3*atan2(b1,yppa1);
+      up -= c4*atan2(b2,yppa2);  
+      
+      return up;
+    }
+    
+    double dupdyp(double yp) const
+    {
+      double yppa1 = yp + a1;
+      double yppa2 = yp + a2;
+
+      double value  = 2.0*c1*yppa1/(pow(yppa1,2.0) + pow(b1,2.0));
+      value -= 2.0*c2*yppa2/(pow(yppa2,2.0) + pow(b2,2.0));
+      value -= c3/(yppa1*(1.0+pow(b1/yppa1,2.0)));
+      value -= c4/(yppa2*(1.0+pow(b2/yppa2,2.0)));
+      
+      return value;
+    }
+  };
+  
+ protected:
+  // Newton-Raphson parameters
+  double rtol;                // Residual tolerance
+  int    maxit;               // Maximal number of iterations
+  
+  // Parameters
+  const double K;             // von Karman constant
+  
+  // Function definition
+  f_SA_wall_law SA;           // Spalart-Allmaras parametrization
+
+ public:
+  //! \brief The default constructor initializes parameters.
+ SAWallLaw(double eps = 1.0e-10, int mit = 400) : rtol(eps), maxit(mit), K(0.41) {}
+  
+  //! \brief Empty destructor
+  virtual ~SAWallLaw() {}
+
+  bool computeYplus(double y, double nu, double utan, double& yplus) const;
+  bool computeUstar(double y, double nu, double utan, double& ustar) const;
+  bool computeTauB(double y, double nu, double utan, double& tauB) const;
+  bool computeNu(double y, double nu, double utan, double& nuwall) const;
+
+ 
+};    
+
+#endif