[WIP] Isothermal streamwise periodicity

Common/include/CConfig.hpp

@@ -1044,6 +1044,7 @@ class CConfig {
   su2double Streamwise_Periodic_PressureDrop;        /*!< \brief Value of prescribed pressure drop [Pa] which results in an artificial body force vector. */
   su2double Streamwise_Periodic_TargetMassFlow;      /*!< \brief Value of prescribed massflow [kg/s] which results in an delta p and therefore an artificial body force vector. */
   su2double Streamwise_Periodic_OutletHeat;          /*!< /brief Heatflux boundary [W/m^2] imposed at streamwise periodic outlet. */
+  su2double Streamwise_Periodic_LambdaL;            /*!< /brief Exp coefficient for iso-thermal BCs Streamwise Periodic. */
 
   su2double *FreeStreamTurboNormal;     /*!< \brief Direction to initialize the flow in turbomachinery computation */
   su2double Restart_Bandwidth_Agg;      /*!< \brief The aggregate of the bandwidth for writing binary restarts (to be averaged later). */
@@ -2997,6 +2998,12 @@ class CConfig {
    */
   unsigned short GetnMarker_HeatFlux(void) const { return nMarker_HeatFlux; }
 
+  /*!
+   * \brief Get the total (local) number of isothermal markers.
+   * \return Total number of isothermal markers.
+   */
+  unsigned short GetnMarker_Isothermal(void) const { return nMarker_Isothermal; }
+
   /*!
    * \brief Get the total number of rough markers.
    * \return Total number of heat flux markers.
@@ -9605,4 +9612,16 @@ class CConfig {
    */
   SA_ParsedOptions GetSAParsedOptions() const { return saParsedOptions; }
 
+  /*!
+   * \brief Set Lambda L for Streamwise Periodic
+   * \return bool
+   */
+  void SetStreamwise_Periodic_LamdaL(su2double value) { Streamwise_Periodic_LambdaL = value; }
+
+  /*!
+   * \brief Set Lambda L for Streamwise Periodic
+   * \return bool
+   */
+  su2double GetStreamwise_Periodic_LamdaL(void) const { return Streamwise_Periodic_LambdaL; }
+
 };

Common/include/option_structure.hpp

@@ -1778,6 +1778,7 @@ enum ENUM_OBJECTIVE {
   TOPOL_DISCRETENESS = 63,      /*!< \brief Measure of the discreteness of the current topology. */
   TOPOL_COMPLIANCE = 64,        /*!< \brief Measure of the discreteness of the current topology. */
   STRESS_PENALTY = 65,          /*!< \brief Penalty function of VM stresses above a maximum value. */
+  STREAMWISE_PERIODIC_LAMBDAL= 71,
 };
 static const MapType<std::string, ENUM_OBJECTIVE> Objective_Map = {
   MakePair("DRAG", DRAG_COEFFICIENT)
@@ -1813,6 +1814,7 @@ static const MapType<std::string, ENUM_OBJECTIVE> Objective_Map = {
   MakePair("SURFACE_PRESSURE_DROP", SURFACE_PRESSURE_DROP)
   MakePair("SURFACE_SPECIES_0", SURFACE_SPECIES_0)
   MakePair("SURFACE_SPECIES_VARIANCE", SURFACE_SPECIES_VARIANCE)
+  MakePair("STREAMWISE_PERIODIC_LAMBDAL", STREAMWISE_PERIODIC_LAMBDAL)
   MakePair("CUSTOM_OBJFUNC", CUSTOM_OBJFUNC)
   MakePair("FLOW_ANGLE_OUT", FLOW_ANGLE_OUT)
   MakePair("MASS_FLOW_IN", MASS_FLOW_IN)
@@ -2457,6 +2459,8 @@ struct StreamwisePeriodicValues {
   su2double Streamwise_Periodic_InletTemperature;   /*!< \brief Area avg static Temp [K] at the periodic inlet. Used for adaptive outlet heatsink. */
   su2double Streamwise_Periodic_BoundaryArea;       /*!< \brief Global Surface area of the streamwise periodic interface. */
   su2double Streamwise_Periodic_AvgDensity;         /*!< \brief Area avg density on the periodic interface. */
+  su2double Streamwise_Periodic_LambdaL;            /*!< \brief Temperature Gradient in iso-thermal BCs. */
+  su2double Streamwise_Periodic_ThetaScaling;       /*!< \brief Scaled Temperature for iso-thermal BCs. */
 };
 
 /*!

Common/src/CConfig.cpp

@@ -4957,8 +4957,8 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
       SU2_MPI::Error("Streamwise Periodic Flow + Incompressible Euler: Not tested yet.", CURRENT_FUNCTION);
     if (nMarker_PerBound == 0)
       SU2_MPI::Error("A MARKER_PERIODIC pair has to be set with KIND_STREAMWISE_PERIODIC != NONE.", CURRENT_FUNCTION);
-    if (Energy_Equation && Streamwise_Periodic_Temperature && nMarker_Isothermal != 0)
-      SU2_MPI::Error("No MARKER_ISOTHERMAL marker allowed with STREAMWISE_PERIODIC_TEMPERATURE= YES, only MARKER_HEATFLUX & MARKER_SYM.", CURRENT_FUNCTION);
+    // if (Energy_Equation && Streamwise_Periodic_Temperature && nMarker_Isothermal != 1)
+    //   SU2_MPI::Error("No MARKER_ISOTHERMAL marker allowed with STREAMWISE_PERIODIC_TEMPERATURE= YES, only MARKER_HEATFLUX & MARKER_SYM.", CURRENT_FUNCTION);
     if (Ref_Inc_NonDim != DIMENSIONAL)
       SU2_MPI::Error("Streamwise Periodicity only works with \"INC_NONDIM= DIMENSIONAL\", the nondimensionalization with source terms doesn;t work in general.", CURRENT_FUNCTION);
     if (Axisymmetric)

SU2_CFD/include/numerics/flow/flow_sources.hpp

@@ -343,6 +343,8 @@ class CSourceIncStreamwise_Periodic final : public CSourceBase_Flow {
   bool turbulent; /*!< \brief Turbulence model used. */
   bool energy;    /*!< \brief Energy equation on. */
   bool streamwisePeriodic_temperature; /*!< \brief Periodicity in energy equation */
+  bool bool_heat_flux_bc; /*!< \brief HEAT Flux BC boolean. */
+  bool bool_isotherml_bc; /*!< \brief ISO THERMAL BC boolean. */
   su2double Streamwise_Coord_Vector[MAXNDIM] = {0.0}; /*!< \brief Translation vector between streamwise periodic surfaces. */
 
   su2double norm2_translation, /*!< \brief Square of distance between the 2 periodic surfaces. */

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -2980,6 +2980,8 @@ su2double CFVMFlowSolverBase<V,R>::EvaluateCommonObjFunc(const CConfig& config)
     case SURFACE_SPECIES_VARIANCE:
       objFun += weight * config.GetSurface_Species_Variance(0);
       break;
+    case STREAMWISE_PERIODIC_LAMBDAL:
+    objFun += weight * config.GetStreamwise_Periodic_LamdaL();
     case CUSTOM_OBJFUNC:
       objFun += weight * Total_Custom_ObjFunc;
       break;

SU2_CFD/src/numerics/flow/flow_sources.cpp

@@ -703,6 +703,8 @@ CSourceIncStreamwise_Periodic::CSourceIncStreamwise_Periodic(unsigned short val_
   turbulent = (config->GetKind_Turb_Model() != TURB_MODEL::NONE);
   energy    = config->GetEnergy_Equation();
   streamwisePeriodic_temperature = config->GetStreamwise_Periodic_Temperature();
+  bool_heat_flux_bc = (config->GetnMarker_HeatFlux() > 0);
+  bool_isotherml_bc = (config->GetnMarker_Isothermal() > 0);
 
   for (unsigned short iDim = 0; iDim < nDim; iDim++)
     Streamwise_Coord_Vector[iDim] = config->GetPeriodic_Translation(0)[iDim];
@@ -729,21 +731,42 @@ CNumerics::ResidualType<> CSourceIncStreamwise_Periodic::ComputeResidual(const C
   /*--- Compute the periodic temperature contribution to the energy equation, if energy equation is considered ---*/
   if (energy && streamwisePeriodic_temperature) {
 
-    scalar_factor = SPvals.Streamwise_Periodic_IntegratedHeatFlow * DensityInc_i / (SPvals.Streamwise_Periodic_MassFlow * norm2_translation);
+    /*--- Compute scalar factors, if there is a HEAT FLUX BC ---*/
+    if (bool_heat_flux_bc) {
+      scalar_factor = SPvals.Streamwise_Periodic_IntegratedHeatFlow * DensityInc_i / (SPvals.Streamwise_Periodic_MassFlow * norm2_translation);
 
-    /*--- Compute scalar-product dot_prod(v*t) ---*/
-    dot_product = GeometryToolbox::DotProduct(nDim, Streamwise_Coord_Vector, &V_i[1]);
+      /*--- Compute scalar-product dot_prod(v*t) ---*/
+      dot_product = GeometryToolbox::DotProduct(nDim, Streamwise_Coord_Vector, &V_i[1]);
+    }
+
+    /*--- Compute scalar factors, if there is an ISOTHERMAL BC ---*/
+    if (bool_isotherml_bc) {
+
+      /*--- Compute three terms associated with the source terms for iso-thermal BCs ---*/
+      // Reference Eq(20) Stalio et. al, doi:10.1115/1.2717235
+
+      dot_product = GeometryToolbox::DotProduct(nDim, Streamwise_Coord_Vector, PrimVar_Grad_i[3]);
+      su2double term1 = 2 * Thermal_Conductivity_i * dot_product;
+      su2double term2 = - SPvals.Streamwise_Periodic_LambdaL * Thermal_Conductivity_i * V_i[3];
+      su2double term3 = - V_i[3] * V_i[1] * DensityInc_i *  config->GetSpecific_Heat_Cp();
+      scalar_factor = SPvals.Streamwise_Periodic_LambdaL * (term1 + term2 + term3);
+      dot_product = 1.0;
+    }
 
     residual[nDim+1] = Volume * scalar_factor * dot_product;
 
     /*--- If a RANS turbulence model is used, an additional source term, based on the eddy viscosity gradient is added. ---*/
     if(turbulent) {
+      if (bool_heat_flux_bc) {
+        /*--- Compute a scalar factor ---*/
+        scalar_factor = SPvals.Streamwise_Periodic_IntegratedHeatFlow / (SPvals.Streamwise_Periodic_MassFlow * sqrt(norm2_translation) * Prandtl_Turb);
+      }
+
+      if (bool_isotherml_bc)
+        scalar_factor = (-V_i[3] * SPvals.Streamwise_Periodic_LambdaL) * config->GetSpecific_Heat_Cp() / Prandtl_Turb;
 
-      /*--- Compute a scalar factor ---*/
-      scalar_factor = SPvals.Streamwise_Periodic_IntegratedHeatFlow / (SPvals.Streamwise_Periodic_MassFlow * sqrt(norm2_translation) * Prandtl_Turb);
-
-      /*--- Compute scalar product between periodic translation vector and eddy viscosity gradient. ---*/
-      dot_product = GeometryToolbox::DotProduct(nDim, Streamwise_Coord_Vector, AuxVar_Grad_i[0]);
+        /*--- Compute scalar product between periodic translation vector and eddy viscosity gradient. ---*/
+        dot_product = GeometryToolbox::DotProduct(nDim, Streamwise_Coord_Vector, AuxVar_Grad_i[0]);
 
       residual[nDim+1] -= Volume * scalar_factor * dot_product;
     } // if turbulent

SU2_CFD/src/output/CFlowIncOutput.cpp

@@ -178,6 +178,8 @@ void CFlowIncOutput::SetHistoryOutputFields(CConfig *config){
     AddHistoryOutput("STREAMWISE_MASSFLOW", "SWMassflow", ScreenOutputFormat::FIXED, "STREAMWISE_PERIODIC", "Massflow in streamwise periodic flow");
     AddHistoryOutput("STREAMWISE_DP",       "SWDeltaP",   ScreenOutputFormat::FIXED, "STREAMWISE_PERIODIC", "Pressure drop in streamwise periodic flow");
     AddHistoryOutput("STREAMWISE_HEAT",     "SWHeat",     ScreenOutputFormat::FIXED, "STREAMWISE_PERIODIC", "Integrated heat for streamwise periodic flow");
+    if(config->GetnMarker_Isothermal() > 0)
+        AddHistoryOutput("STREAMWISE_LAMBDAL",     "SWLambdaL",     ScreenOutputFormat::FIXED, "STREAMWISE_PERIODIC", "Exponential factor for iso-thermal temperature gradient");
   }
 
   AddAnalyzeSurfaceOutput(config);
@@ -262,6 +264,8 @@ void CFlowIncOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSolv
     SetHistoryOutputValue("STREAMWISE_MASSFLOW", flow_solver->GetStreamwisePeriodicValues().Streamwise_Periodic_MassFlow);
     SetHistoryOutputValue("STREAMWISE_DP", flow_solver->GetStreamwisePeriodicValues().Streamwise_Periodic_PressureDrop);
     SetHistoryOutputValue("STREAMWISE_HEAT", flow_solver->GetStreamwisePeriodicValues().Streamwise_Periodic_IntegratedHeatFlow);
+    if(config->GetnMarker_Isothermal() > 0)
+        SetHistoryOutputValue("STREAMWISE_LAMBDAL", flow_solver->GetStreamwisePeriodicValues().Streamwise_Periodic_LambdaL);
   }
 
   /*--- Set the analyse surface history values --- */

SU2_CFD/src/output/CFlowOutput.cpp

@@ -451,7 +451,7 @@ void CFlowOutput::SetAnalyzeSurface(const CSolver* const*solver, const CGeometry
     Tot_Surface_Temperature += Temperature;
     config->SetSurface_Temperature(iMarker_Analyze, Temperature);
 
-    su2double Pressure = Surface_Pressure_Total[iMarker_Analyze] * config->GetPressure_Ref();
+    su2double Pressure = Surface_Pressure[iMarker_Analyze] * config->GetPressure_Ref();
     SetHistoryOutputPerSurfaceValue("SURFACE_STATIC_PRESSURE", Pressure, iMarker_Analyze);
     Tot_Surface_Pressure += Pressure;
     config->SetSurface_Pressure(iMarker_Analyze, Pressure);

SU2_CFD/src/solvers/CIncEulerSolver.cpp

@@ -1643,6 +1643,9 @@ void CIncEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_cont
       if(streamwise_periodic_temperature && turbulent)
         numerics->SetAuxVarGrad(nodes->GetAuxVarGradient(iPoint), nullptr);
 
+      /*--- Load the Prim variable gradient that we already computed. ---*/
+      numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), nullptr);
+
       /*--- Compute the streamwise periodic source residual and add to the total ---*/
       auto residual = numerics->ComputeResidual(config);
       LinSysRes.AddBlock(iPoint, residual);

SU2_CFD/src/solvers/CIncNSSolver.cpp

@@ -120,7 +120,7 @@ void CIncNSSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container
 void CIncNSSolver::GetStreamwise_Periodic_Properties(const CGeometry *geometry,
                                                      CConfig *config,
                                                      const unsigned short iMesh) {
-
+    bool turbulent = (config->GetKind_Turb_Model() != TURB_MODEL::NONE);
   /*---------------------------------------------------------------------------------------------*/
   // 1. Evaluate massflow, area avg density & Temperature and Area at streamwise periodic outlet.
   // 2. Update delta_p is target massflow is chosen.
@@ -199,6 +199,13 @@ void CIncNSSolver::GetStreamwise_Periodic_Properties(const CGeometry *geometry,
     /*---------------------------------------------------------------------------------------------*/
 
     su2double HeatFlow_Local = 0.0, HeatFlow_Global = 0.0;
+    su2double dTdn_Local = 0.0, dTdn_Global = 0.0;
+    su2double Volume_Temp_Local = 0.0, Volume_Temp_Global = 0.0;
+    su2double Volume_TempS_Local = 0.0, Volume_TempS_Global = 0.0;
+    su2double Volume_Local = 0.0, Volume_Global = 0.0;
+    su2double Volume_VTemp_Local = 0.0, Volume_VTemp_Global = 0.0;
+    su2double turb_b1_coeff_Local = 0.0, turb_b1_coeff_Global = 0.0;
+    su2double turb_b2_coeff_Local = 0.0, turb_b2_coeff_Global = 0.0;
 
     /*--- Loop over all heatflux Markers ---*/
     for (auto iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
@@ -222,13 +229,90 @@ void CIncNSSolver::GetStreamwise_Periodic_Properties(const CGeometry *geometry,
           HeatFlow_Local += FaceArea * (-1.0) * Wall_HeatFlux/config->GetHeat_Flux_Ref();;
         } // loop Vertices
       } // loop Heatflux marker
+
+      if (config->GetMarker_All_KindBC(iMarker) == ISOTHERMAL) {
+        /*--- Identify the boundary by string name and retrive ISOTHERMAL from config ---*/
+        su2double GradT[3] = {0.0,0.0,0.0};
+
+        for (auto iVertex = 0ul; iVertex < geometry->nVertex[iMarker]; iVertex++) {
+
+          auto iPoint = geometry->vertex[iMarker][iVertex]->GetNode();
+
+          if (!geometry->nodes->GetDomain(iPoint)) continue;
+
+          const auto AreaNormal = geometry->vertex[iMarker][iVertex]->GetNormal();
+
+          /*Compute wall heat flux (normal to the wall) based on computed temperature gradient*/
+          for(auto iDim=0; iDim < nDim; iDim++) {
+
+            GradT[iDim] = nodes->GetGradient_Primitive(iPoint, 3, iDim);
+
+            dTdn_Local += GradT[iDim]*AreaNormal[iDim]*nodes->GetThermalConductivity(iPoint);
+
+          } // loop dim
+        } // loop Vertices
+      } // loop Isothermal Marker
     } // loop AllMarker
 
     /*--- MPI Communication sum up integrated Heatflux from all processes ---*/
     SU2_MPI::Allreduce(&HeatFlow_Local, &HeatFlow_Global, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
+    SU2_MPI::Allreduce(&dTdn_Local, &dTdn_Global, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
+
+    /*--- Loop over all heatflux Markers ---*/
+    for (unsigned long iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) {
+
+      if (!geometry->nodes->GetDomain(iPoint)) continue;
+
+      su2double volume = geometry->nodes->GetVolume(iPoint);
+
+      su2double Temp = nodes->GetTemperature(iPoint);
+
+      Volume_Local += volume;
+
+      Volume_TempS_Local += volume * Temp;
+
+      Volume_Temp_Local += volume * Temp * nodes->GetThermalConductivity(iPoint);
+
+      // coeff_b1 for turbulence
+      if (turbulent && (config->GetnMarker_Isothermal() != 0))
+        turb_b1_coeff_Local += Temp * nodes->GetAuxVarGradient(iPoint, 0, 0) * config->GetSpecific_Heat_Cp() * volume  / config->GetPrandtl_Turb();
+
+      Volume_VTemp_Local += volume * Temp * nodes->GetVelocity(iPoint, 0) * nodes->GetDensity(iPoint);
+
+    } // points
+
+    /*--- MPI Communication sum up integrated Heatflux from all processes ---*/
+    SU2_MPI::Allreduce(&Volume_Temp_Local, &Volume_Temp_Global, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
+    SU2_MPI::Allreduce(&Volume_VTemp_Local, &Volume_VTemp_Global, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
+    SU2_MPI::Allreduce(&Volume_Local, &Volume_Global, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
+    SU2_MPI::Allreduce(&Volume_TempS_Local, &Volume_TempS_Global, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
+    if (turbulent && (config->GetnMarker_Isothermal() != 0))
+      SU2_MPI::Allreduce(&turb_b1_coeff_Local, &turb_b1_coeff_Global, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
 
     /*--- Set the solver variable Integrated Heatflux ---*/
-    SPvals.Streamwise_Periodic_IntegratedHeatFlow = HeatFlow_Global;
+    if (config->GetnMarker_HeatFlux() > 0)
+      SPvals.Streamwise_Periodic_IntegratedHeatFlow = HeatFlow_Global;
+
+    if (config->GetnMarker_Isothermal() > 0) {
+      SPvals.Streamwise_Periodic_ThetaScaling = Volume_TempS_Global/Volume_Global;
+
+      for (unsigned long iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++)
+        nodes->SetStreamwise_Periodic_RecoveredTemperature(iPoint, nodes->GetTemperature(iPoint)/SPvals.Streamwise_Periodic_ThetaScaling);
+
+      /*--- Set the solver variable Lambda_L for iso-thermal BCs ---*/
+      su2double b0_coeff =  Volume_Temp_Global; 
+      su2double b1_coeff = Volume_VTemp_Global * config->GetSpecific_Heat_Cp() + turb_b1_coeff_Global;
+      su2double b2_coeff = -dTdn_Global;
+
+      /*--- Find the value of Lambda L by solving the quadratic equation ---*/
+      su2double pred_lambda = (- b1_coeff + sqrt(b1_coeff * b1_coeff - 4 * b0_coeff * b2_coeff))/(2 * b0_coeff);
+      if (!config->GetDiscrete_Adjoint())
+        SPvals.Streamwise_Periodic_LambdaL -= 0.01 * (SPvals.Streamwise_Periodic_LambdaL - pred_lambda);
+      else
+        SPvals.Streamwise_Periodic_LambdaL = pred_lambda;
+    //   cout<<"Lambda :: "<<SPvals.StrLambdaeamwise_Periodic_LambdaL<<endl;
+      config->SetStreamwise_Periodic_LamdaL(SPvals.Streamwise_Periodic_LambdaL);
+    } // if isothermal
   } // if energy
 }