diff --git a/Apps/Common/CFDenums.h b/Apps/Common/CFDenums.h index 162a3316..8afe4774 100644 --- a/Apps/Common/CFDenums.h +++ b/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 }; } diff --git a/Apps/Common/CMakeLists.txt b/Apps/Common/CMakeLists.txt index fc051d16..3a46a2df 100644 --- a/Apps/Common/CMakeLists.txt +++ b/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}) diff --git a/Apps/Common/SAWallLaw.C b/Apps/Common/SAWallLaw.C new file mode 100644 index 00000000..3758c4b9 --- /dev/null +++ b/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; +} diff --git a/Apps/Common/SAWallLaw.h b/Apps/Common/SAWallLaw.h new file mode 100644 index 00000000..b18271ea --- /dev/null +++ b/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 +#include + +/*! + \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