gmpMaterialPipeFlow.h

Go to the documentation of this file.

/************************************************************************
**
** Copyright (C) 2014 by Carlos Augusto Teixera Mendes
** All rights reserved.
**
** This file is part of the "GeMA" software. It's use should respect
** the terms in the license agreement that can be found together
** with this source code.
** It is provided AS IS, with NO WARRANTY OF ANY KIND,
** INCLUDING THE WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR
** A PARTICULAR PURPOSE.
**
************************************************************************/
 
#ifndef _GEMA_PLUGIN_HYDRAULIC_MATERIAL_PIPEFLOW_H_
#define _GEMA_PLUGIN_HYDRAULIC_MATERIAL_PIPEFLOW_H_
 
#include <gmpFemPhysicsCommonMaterial.h>
 
#include "gmpHydraulicConfig.h"
#include "gmpHydraulicPipe.h"
#include <gmMathUtils.h>
#include <gmTrace.h>
 
 
 
class GMP_HYDRAULIC_PHYSICS_API_EXPORT GmpMaterialPipeFlow: public GmpFemPhysicsCommonMaterial
{
protected:
  enum PipePropertyIds
  {
    Dc_ID,        
    Uf_ID,        
    MFw_ID,       
    GW_ID,                              
    W_ID,                                 
    rhof_ID,                    
  };
 
  enum pipeGaussAttrIds
  {
          FR_GA_ID,          
          FROLD_GA_ID,       
          QFLOW_GA_ID,
          RE_GA_ID,                         
 
          // ------    NO ADDING BELOW THIS LINE
          NUM_GA_IDS,          
  };
 
public:
 
  GmpMaterialPipeFlow(int typeIndex, QString typeName, const GmLogCategory& logger)
    : GmpFemPhysicsCommonMaterial(typeIndex, typeName, logger) {}
 
  static GmpFemPhysicsCommonMaterial* instance(GmSimulationData* simulation, int typeIndex, QString typeName, const GmLogCategory& logger)
  {
    Q_UNUSED(simulation);
    return new GmpMaterialPipeFlow(typeIndex, typeName, logger);
  }
 
  // set pipe material map
  virtual const QVariantMap* materialMetaDataMap()
  {
    S_TRACE();
 
    static QVariantMap m;
 
    if (m.isEmpty())
    {
          // Adds material properties
      m["properties"] = GmpFemPhysicsCommon::ValueList()
        << GmpFemPhysicsCommon::ScalarValue(Dc_ID, "Dp", QObject::tr("Pipe diameter"), "m", true)
        << GmpFemPhysicsCommon::ScalarValue(Uf_ID, "uf", QObject::tr("Dynamic fluid viscosity"), "kPa*s", true)
        << GmpFemPhysicsCommon::ScalarValue(MFw_ID, "Mfw", QObject::tr("Interface biot modulus"), "kPa/m", false)
        << GmpFemPhysicsCommon::ScalarValue(GW_ID, "gw", QObject::tr("Fluid specific weight"), "kN/m3", true)
        << GmpFemPhysicsCommon::ScalarValue(W_ID, "wr", QObject::tr("Wall roughness"), "m", true)
        << GmpFemPhysicsCommon::ScalarValue(rhof_ID, "rhow", QObject::tr("Water density"), "kg/m3", true);
 
 
          // Adds material gauss attributes
          m["gaussAttributes"] = GmpFemPhysicsCommon::ValueList()
                  << GmpFemPhysicsCommon::ScalarValue(FR_GA_ID, "FR", QObject::tr("flow regime"), " ", true, false, 2)
                  << GmpFemPhysicsCommon::HistoryValue(FROLD_GA_ID, FR_GA_ID, 1, true)
                  << GmpFemPhysicsCommon::ScalarValue(QFLOW_GA_ID, "Qf", QObject::tr("Flow in the pipe element"), "m3/s", true, false, 1)
                  << GmpFemPhysicsCommon::ScalarValue(RE_GA_ID, "Re", QObject::tr("Reynolds number of the flow in the pipe"), " ", true, false, 1);
    }
 
    return &m;
  }
 
  virtual double pipeDiameter(const GmElement* e, const GmVector* coord, int ip) const
  {
    S_TRACE();
    assert(e);
    return propertyAc(Dc_ID)->scalarValueAt(e, coord, ip);
  }
 
  virtual double fluidViscosity(const GmElement* e, const GmVector* coord, int ip) const
  {
    S_TRACE();
    assert(e);
    return propertyAc(Uf_ID)->scalarValueAt(e, coord, ip);
  }
 
  virtual double fluidSpecificWeight(const GmElement* e, const GmVector* coord, int ip) const
  {
          S_TRACE();
          assert(e);
          return propertyAc(GW_ID)->scalarValueAt(e, coord, ip);
  }
 
  virtual double pipeWallRoughness(const GmElement* e, const GmVector* coord, int ip) const
  {
          S_TRACE();
          assert(e);
          return propertyAc(W_ID)->scalarValueAt(e, coord, ip);
  }
 
  virtual double fluidDensity(const GmElement* e, const GmVector* coord, int ip) const
  {
          S_TRACE();
          assert(e);
 
          return propertyAc(rhof_ID)->scalarValueAt(e, coord, ip);
  }
 
  virtual double frictionFactor(const GmElement* e, const GmVector* coord, int ip, double dp) const
  {
          S_TRACE();
          assert(e);
          // Get properties
          GmGaussAccessor* ReAcc = gaussAttrAc(RE_GA_ID);
          double w = pipeWallRoughness(e, coord, ip);
          double NRe = ReAcc->scalarValueAt(e, ip, coord);
          double fs, ft, fr, fmax, R, dR, deltaft;
          int it = 1, iterMAX = 30;
 
 
          // 1) Compute the friction factor for smooth turbulence, fs
           fs = 0.3164*pow(NRe, -0.25);
 
          // 2) Compute the friction factor with respect to the transition between smooth and rough turbulence, ft
          // Initialize thhe parameters
          ft = fs; 
          R = pow(ft, -0.50) - 1.14 + 2.0*log10(w / dp + 9.35 / NRe * pow(ft, -0.50));
 
          // Newton-Raphson's procedure
          while (abs(R) >= 1e-10)
          {
                  // Derivative of residual function R with respect dgamma
                  dR = -0.50 / pow(ft, 1.50) - 9.350 / (NRe*pow(ft, 1.50)*(w / dp + 9.35 / NRe * pow(ft, -0.50))*log(10.0));
 
                  // Increment of dgamma
                  deltaft = -R / dR;
 
                  // Update
                  ft += deltaft;
 
                  // Check the residual function
                  R = pow(ft, -0.50) - 1.14 + 2.0*log10(w / dp + 9.35 / NRe * pow(ft, -0.50));
 
                  if (it == iterMAX)
                  {
                          gmErrorMsg(logger(), "Error to compute the friction factor. No convergence");
                  }
 
                  it++;
 
          }
 
          // 3) Compute the friction factor for rough turbulence, fr
          fr = pow(1 / (1.14 - 2.0*log10(w / dp)),2.0);
 
          // The friction factor is the maximmum among fs, ft and fr
          fmax = (fs > ft) ? fs : ft;
          fmax = (fmax > fr) ? fmax : fr;
          
          return fmax;
  }
 
  virtual double interfaceBiotModulus(const GmElement* e, const GmVector* coord, int ip) const
  {
    S_TRACE();
    assert(e);
    if (propertyAc(MFw_ID) == NULL)
    {
      return 0;
    }
    return propertyAc(MFw_ID)->scalarValueAt(e, coord, ip);
  }
 
  virtual double pipeConductivity(const GmElement* e, const GmVector* coord, int ip, GmGaussAccessor* DpAcc, GmVector& Pe, double L) const
  {
    S_TRACE();
    assert(e);
        GmGaussAccessor* oldFrAcc = gaussAttrAc(FROLD_GA_ID);        // read from old FR
        double isTurbulent = oldFrAcc->scalarValueAt(e, ip, coord);
    double kl = 0.0;
    double dp = DpAcc->scalarValueAt(e, ip, coord);
    double uf = fluidViscosity(e, coord, ip);
        double gw = fluidSpecificWeight(e, coord, ip);
        double rhoW = fluidDensity(e, coord, ip);
        double w = pipeWallRoughness(e, coord, ip);
        double dH = abs(Pe(0) - Pe(1))/gw;  // head loss in length units
 
        // Compute the gravitational acceleration
        double gAc = 1000.0 * gw / rhoW; // m/s2   
 
    // computes the cubic law
    if (uf > 0.0)
    {
                if (dp > 0.0)
                {
 
                        // Laminar case
                        if (isTurbulent==0)
                        {
                            kl = arma::datum::pi*pow(dp, 4.0) / (128.0 * uf);
                         }
 
                        // Turbulent case
                        // For turbulent the pipe conductivity is computed based on the Darcy-Weisbach and Colebrook equations as presented
                        // by Swammee and Jain (1976). In this expression, the friction factor is explicitly shown and it is independent of 
                        // the Reynolds number.
                        else
                        {
                          
                                double nu = 1000.0*uf / rhoW; // m2/s
                                double A = w / (3.7*dp);
                                double B = 2.51*nu / dp;
                                double alpha = 2.0*gAc*dp*abs(dH) / L;
                                double f1 = -arma::datum::pi * pow(dp, 2.0)*sqrt(2.0*gAc*dp*L / abs(dH)) / (2.0*gw);
                                double f2 = log10(A + B / sqrt(alpha));
                                
                                kl =   f1 * f2;
 
                        }
 
                        // According to Kaufmann and Braun (1999)
                        //else
                        //{
                        //      // Get the friction factor
                        //      double f = frictionFactor(e, coord, ip, dp);
                        //      double alpha = gAc * dp*L / (8.0*f*dH);
                        //      kl = -arma::datum::pi*pow(dp, 2.0)*sqrt(alpha)/ (gw);
                        //}
 
                // Swamee & Swammee (2007)
                        //if (dH < 1e-10)
                        //{
                        //      dH = 1e-10;
                        //}
 
                 // double nu = 1000.0*uf / rhoW; // m2/s
                 // double alpha = gAc * dp*dH / L;
                 // double C1 = pow(128 * nu / (arma::datum::pi*dp*sqrt(alpha)), 4.0);
                 // double C2 = 1.153*pow(pow(415 * nu / (dp*sqrt(alpha)), 8.0) - log(w / (3.7*dp) + 1.775*nu / (dp*sqrt(alpha))), -4.0);
                 // double C = dp * dp / gw * pow(C1 + C2, -0.25);
                 // kl = C * sqrt(gAc*dp*L / (dH));
 
          
      }
    }
    else
    {
      assert(uf != 0 && "Ufw Dynamic viscosity must be different to zero");
    }
 
    return kl;
  }
 
  virtual double pipeTangentConductivity(const GmElement* e, const GmVector* coord, int ip, GmGaussAccessor* DpAcc, GmVector& Pe, double L) const
  {
          S_TRACE();
          assert(e);
          GmGaussAccessor* oldFrAcc = gaussAttrAc(FROLD_GA_ID); // read from old FR
          double isTurbulent = oldFrAcc->scalarValueAt(e, ip, coord);
          double kt = 0.0;
          double dp = DpAcc->scalarValueAt(e, ip, coord);
          double uf = fluidViscosity(e, coord, ip);
          double gw = fluidSpecificWeight(e, coord, ip);
          double rhoW = fluidDensity(e, coord, ip);
          double w = pipeWallRoughness(e, coord, ip);
          double dH = (Pe(0) - Pe(1))/ gw;  // head loss in length units
          double signaldH;
 
          if (dH >= 0)
          {
                  signaldH = 1.0;
 
          }
          else
          {
                  signaldH = -1.0;
          }
 
          // Compute the gravitational acceleration
          double gAc = 1000.0*gw / rhoW; // m2/s
 
          // computes the cubic law
          if (uf > 0.0)
          {
                  if (dp > 0.0)
                  {
                          // Laminar case
                          if (isTurbulent == 0)
                          {
 
                                  kt = arma::datum::pi * pow(dp, 4.0) / (128.0 * uf);
                          }
 
                          // Turbulent case
                          // For turbulent the pipe conductivity is computed based on the Darcy-Weisbach and Colebrook equations as presented
                          // by Swammee and Jain (1976). In this expression, the friction factor is explicitly shown and it is independent of 
                          // the Reynolds number.
                          else
                          {
                              double nu = 1000.0*uf / rhoW; // m2/s
                                  double A = w / (3.7*dp);
                                  double B = 2.51*nu / dp;
                                  double alpha = 2.0*gAc*dp*abs(dH) / L;
                                  double f1 = -arma::datum::pi * pow(dp, 2.0)*sqrt(2.0*gAc*dp*L / abs(dH)) / (2.0*gw);
                                  double f2 = log10(A + B / sqrt(alpha));
                                  double df1h1 = signaldH * arma::datum::pi*pow(dp, 3.0)*gAc*L / (4.0*gw*pow(dH, 2.0))*sqrt(2.0*abs(dH) / (gAc*dp*L));
                                  double df2h1 = -1.0*signaldH*B / (2.0*abs(dH)*(A*sqrt(alpha) + B)*log(10));
 
                                  kt = dH * (df1h1*f2 + f1 * df2h1) + f1 * f2;
                          }
 
                  }
          }
          else
          {
                  assert(uf != 0 && "Ufw Dynamic viscosity must be different to zero");
          }
 
          return kt;
  }
 
  virtual void ckeckFlowRegime(const GmElement* e, const GmVector* coord, int ip, GmGaussAccessor* DpAcc, GmVector& Fi) const
  {
          S_TRACE();
          GmGaussAccessor* FrAcc    = gaussAttrAc(FR_GA_ID); // save no new FR
          double Qflow = abs(Fi(0));
          double dp = DpAcc->scalarValueAt(e, ip, coord);
          double uf = fluidViscosity(e, coord, ip);
          double rhoW = fluidDensity(e, coord, ip);
          double nu = 1000.0*uf / rhoW; // m2/s
          double Re = 4.0*Qflow / (arma::datum::pi*nu*dp); // Compute the Reynolds number
 
 
          // Check if the flow is turbulent (Re > 2200)
          if (Re > 2200)
          {
                  FrAcc->setScalarValue(e, ip, 1);  // turbulent
          }
          else
          {
                  FrAcc->setScalarValue(e, ip, 0);  // laminar
          }
 
          GmGaussAccessor* QfAcc = gaussAttrAc(QFLOW_GA_ID);
          QfAcc->setScalarValue(e, ip, Qflow);
 
          // Set the Reynolds number 
          GmGaussAccessor* ReAcc = gaussAttrAc(RE_GA_ID);
          ReAcc->setScalarValue(e, ip, Re);
  }
 
};
 
#endif