Merge pull request #1383 from su2code/feature_NEMO_Jhs

Addition of species diffusion term to surface heat flux

Author: jtneedels

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -2392,9 +2392,8 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
 
   unsigned long iVertex, iPoint, iPointNormal;
   unsigned short iMarker, iMarker_Monitoring, iDim, jDim;
-  unsigned short T_INDEX = 0, TVE_INDEX = 0, VEL_INDEX = 0;
-  su2double Viscosity = 0.0, WallDist[3] = {0.0}, Area, TauNormal, dTn, dTven,
-            GradTemperature, Density = 0.0, WallDistMod, FrictionVel,
+  unsigned short T_INDEX = 0, TVE_INDEX = 0, VEL_INDEX = 0, RHOS_INDEX = 0;
+  su2double Viscosity = 0.0, WallDist[3] = {0.0}, Area, Density = 0.0, GradTemperature = 0.0,
             UnitNormal[3] = {0.0}, TauElem[3] = {0.0}, TauTangent[3] = {0.0}, Tau[3][3] = {{0.0}}, Cp,
             thermal_conductivity, MaxNorm = 8.0, Grad_Vel[3][3] = {{0.0}}, Grad_Temp[3] = {0.0}, AxiFactor;
   const su2double *Coord = nullptr, *Coord_Normal = nullptr, *Normal = nullptr;
@@ -2527,30 +2526,25 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
       /*--- Compute wall shear stress (using the stress tensor). Compute wall skin friction coefficient, and heat flux
        * on the wall ---*/
 
-      TauNormal = 0.0;
-      for (iDim = 0; iDim < nDim; iDim++) TauNormal += TauElem[iDim] * UnitNormal[iDim];
+      su2double TauNormal = GeometryToolbox::DotProduct(nDim, TauElem, UnitNormal);
 
-      WallShearStress[iMarker][iVertex] = 0.0;
       for (iDim = 0; iDim < nDim; iDim++) {
         TauTangent[iDim] = TauElem[iDim] - TauNormal * UnitNormal[iDim];
         /* --- in case of wall functions, we have computed the skin friction in the turbulence solver --- */
         /* --- Note that in the wall model, we switch off the computation when the computed y+ < 5    --- */
         /* --- We put YPlus to 1.0 so we have to compute skinfriction and the actual y+ in that case as well --- */
         if (!wallfunctions || (wallfunctions && YPlus[iMarker][iVertex] < 5.0))
           CSkinFriction[iMarker](iVertex,iDim) = TauTangent[iDim] * factorFric;
-
-        WallShearStress[iMarker][iVertex] += TauTangent[iDim] * TauTangent[iDim];
       }
-      WallShearStress[iMarker][iVertex] = sqrt(WallShearStress[iMarker][iVertex]);
+      WallShearStress[iMarker][iVertex] = GeometryToolbox::Norm(nDim, TauTangent);
 
       for (iDim = 0; iDim < nDim; iDim++) WallDist[iDim] = (Coord[iDim] - Coord_Normal[iDim]);
-      WallDistMod = 0.0;
-      for (iDim = 0; iDim < nDim; iDim++) WallDistMod += WallDist[iDim] * WallDist[iDim];
-      WallDistMod = sqrt(WallDistMod);
+      
+      su2double WallDistMod = GeometryToolbox::Norm(nDim, WallDist);
 
       /*--- Compute y+ and non-dimensional velocity ---*/
 
-      FrictionVel = sqrt(fabs(WallShearStress[iMarker][iVertex]) / Density);
+      su2double FrictionVel = sqrt(fabs(WallShearStress[iMarker][iVertex]) / Density);
 
       /* --- in case of wall functions, we have computed YPlus in the turbulence class --- */
       /* --- Note that we do not recompute y+ when y+<5 because y+ can become > 5 again --- */
@@ -2562,35 +2556,41 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
       /// TODO: Move these ifs to specialized functions.
       if (!nemo){
 
-        GradTemperature = 0.0;
-
         if (FlowRegime == ENUM_REGIME::COMPRESSIBLE) {
-          for (iDim = 0; iDim < nDim; iDim++) GradTemperature -= Grad_Temp[iDim] * UnitNormal[iDim];
+          GradTemperature = -GeometryToolbox::DotProduct(nDim, Grad_Temp, UnitNormal);
 
           Cp = (Gamma / Gamma_Minus_One) * Gas_Constant;
           thermal_conductivity = Cp * Viscosity / Prandtl_Lam;
         }
         if (FlowRegime == ENUM_REGIME::INCOMPRESSIBLE) {
           if (energy)
-            for (iDim = 0; iDim < nDim; iDim++) GradTemperature -= Grad_Temp[iDim] * UnitNormal[iDim];
+            GradTemperature = -GeometryToolbox::DotProduct(nDim, Grad_Temp, UnitNormal);
 
           thermal_conductivity = nodes->GetThermalConductivity(iPoint);
         }
-
         HeatFlux[iMarker][iVertex] = -thermal_conductivity * GradTemperature * RefHeatFlux;
 
       } else {
 
-        const auto thermal_conductivity_tr = nodes->GetThermalConductivity(iPoint);
-        const auto thermal_conductivity_ve = nodes->GetThermalConductivity_ve(iPoint);
-        const auto Grad_PrimVar            = nodes->GetGradient_Primitive(iPoint);
-
-        dTn = 0.0; dTven = 0.0;
-        for (iDim = 0; iDim < nDim; iDim++) {
-          dTn   += Grad_PrimVar[T_INDEX][iDim]*UnitNormal[iDim];
-          dTven += Grad_PrimVar[TVE_INDEX][iDim]*UnitNormal[iDim];
-        }
-        HeatFlux[iMarker][iVertex] = thermal_conductivity_tr*dTn + thermal_conductivity_ve*dTven;
+        unsigned short iSpecies, nSpecies = config->GetnSpecies();
+
+        const auto& thermal_conductivity_tr = nodes->GetThermalConductivity(iPoint);
+        const auto& thermal_conductivity_ve = nodes->GetThermalConductivity_ve(iPoint);
+        const auto& Grad_PrimVar            = nodes->GetGradient_Primitive(iPoint);
+        const auto& Ds                      = nodes->GetDiffusionCoeff(iPoint);
+        const auto& hs                      = nodes->GetEnthalpys(iPoint);
+
+        /*--- Compute enthalpy transport to surface due to mass diffusion ---*/ 
+        su2double sumJhs = 0.0;
+        for (iSpecies = 0; iSpecies < nSpecies; iSpecies++)
+          sumJhs += Ds[iSpecies]*hs[iSpecies]*GeometryToolbox::DotProduct(nDim, Grad_PrimVar[RHOS_INDEX+iSpecies], UnitNormal);
+        
+        su2double dTn   = GeometryToolbox::DotProduct(nDim, Grad_PrimVar[T_INDEX], UnitNormal);
+        su2double dTven = GeometryToolbox::DotProduct(nDim, Grad_PrimVar[TVE_INDEX], UnitNormal);
+        
+        /*--- Surface energy balance: trans-rot heat flux, vib-el heat flux,
+        enthalpy transport due to mass diffusion ---*/ 
+        HeatFlux[iMarker][iVertex] = thermal_conductivity_tr*dTn + thermal_conductivity_ve*dTven + sumJhs;
       }
 
       /*--- Note that y+, and heat are computed at the

SU2_CFD/include/variables/CNEMONSVariable.hpp

@@ -42,12 +42,11 @@ class CNEMONSVariable final : public CNEMOEulerVariable {
   VectorType Viscosity_Ref;     /*!< \brief Reference viscosity of the fluid. */
   VectorType Viscosity_Inf;     /*!< \brief Viscosity of the fluid at the infinity. */
   MatrixType DiffusionCoeff;    /*!< \brief Diffusion coefficient of the mixture. */
-  CVectorOfMatrix Dij;            /*!< \brief Binary diffusion coefficients. */
+  CVectorOfMatrix Dij;          /*!< \brief Binary diffusion coefficients. */
   VectorType LaminarViscosity;  /*!< \brief Viscosity of the fluid. */
   VectorType ThermalCond;       /*!< \brief T-R thermal conductivity of the gas mixture. */
   VectorType ThermalCond_ve;    /*!< \brief V-E thermal conductivity of the gas mixture. */
-  vector<su2double> thermalconductivities;
-  vector<su2double> Ds;
+  MatrixType Enthalpys;         /*!< \brief Species enthalpies of the mixture. */
 
   su2double inv_TimeScale;      /*!< \brief Inverse of the reference time scale. */
 
@@ -100,6 +99,12 @@ class CNEMONSVariable final : public CNEMOEulerVariable {
    */
   inline su2double* GetDiffusionCoeff(unsigned long iPoint) override { return DiffusionCoeff[iPoint]; }
 
+  /*!
+   * \brief Get the species enthalpy.
+   * \return Value of the species enthalpy.
+   */
+  inline su2double* GetEnthalpys(unsigned long iPoint) override { return Enthalpys[iPoint]; }
+
   /*!
    * \brief Get the laminar viscosity of the flow.
    * \return Value of the laminar viscosity of the flow.

SU2_CFD/include/variables/CVariable.hpp

@@ -937,6 +937,12 @@ class CVariable {
    */
   inline virtual su2double GetMassFraction(unsigned long iPoint, unsigned long val_Species) const { return 0.0; }
 
+  /*!
+   * \brief Get the species enthalpy.
+   * \return Value of the species enthalpy.
+   */
+  inline virtual su2double* GetEnthalpys(unsigned long iPoint) { return nullptr; }
+
   /*!
    * \brief A virtual member.
    * \param[in] iPoint - Point index.

SU2_CFD/src/variables/CNEMONSVariable.cpp

@@ -61,6 +61,7 @@ CNEMONSVariable::CNEMONSVariable(su2double val_pressure,
   LaminarViscosity.resize(nPoint)          = su2double(0.0);
   ThermalCond.resize(nPoint)               = su2double(0.0);
   ThermalCond_ve.resize(nPoint)            = su2double(0.0);
+  Enthalpys.resize(nPoint, nSpecies)       = su2double(0.0);
 
   Max_Lambda_Visc.resize(nPoint) = su2double(0.0);
   inv_TimeScale = config->GetModVel_FreeStream() / config->GetRefLength();
@@ -77,7 +78,7 @@ CNEMONSVariable::CNEMONSVariable(su2double val_pressure,
 
 bool CNEMONSVariable::SetPrimVar(unsigned long iPoint, CFluidModel *FluidModel) {
 
-  unsigned short iVar, iSpecies;
+  unsigned long iVar, iSpecies;
 
   fluidmodel = static_cast<CNEMOGas*>(FluidModel);
 
@@ -98,13 +99,21 @@ bool CNEMONSVariable::SetPrimVar(unsigned long iPoint, CFluidModel *FluidModel)
 
   SetVelocity2(iPoint);
 
-  Ds = fluidmodel->GetDiffusionCoeff();
+  const auto& Ds = fluidmodel->GetDiffusionCoeff();
   for (iSpecies = 0; iSpecies < nSpecies; iSpecies++)
     DiffusionCoeff(iPoint, iSpecies) = Ds[iSpecies];
 
+  su2double T   =  Primitive(iPoint,nSpecies);
+  su2double Tve =  Primitive(iPoint,nSpecies+1);
+
+  su2double* val_eves = GetEve(iPoint);
+  const auto& hs = fluidmodel->ComputeSpeciesEnthalpy(T, Tve, val_eves);
+  for (iSpecies = 0; iSpecies < nSpecies; iSpecies++)
+    Enthalpys(iPoint, iSpecies) = hs[iSpecies];
+  
   LaminarViscosity(iPoint) = fluidmodel->GetViscosity();
 
-  thermalconductivities    = fluidmodel->GetThermalConductivities();
+  const auto& thermalconductivities = fluidmodel->GetThermalConductivities();
   ThermalCond(iPoint)      = thermalconductivities[0];
   ThermalCond_ve(iPoint)   = thermalconductivities[1];