Thermochemically consistent backflow treatment for FGM

Common/include/CConfig.hpp

@@ -635,6 +635,8 @@ class CConfig {
   unsigned short nInc_Outlet;      /*!< \brief Number of inlet boundary treatment types listed. */
   su2double Inc_Inlet_Damping;     /*!< \brief Damping factor applied to the iterative updates to the velocity at a pressure inlet in incompressible flow. */
   su2double Inc_Outlet_Damping;    /*!< \brief Damping factor applied to the iterative updates to the pressure at a mass flow outlet in incompressible flow. */
+  bool Inc_Outlet_BackflowPrevention;          /*!< \brief Enable removal of the reversed normal velocity component at outlet faces where backflow is detected. */
+  unsigned long Inc_Outlet_BackflowPrevention_Iter; /*!< \brief Number of outer iterations for which backflow prevention is active; after this the BC reverts to pure Neumann. Default ULONG_MAX (always active). */
   bool InletUseNormal;             /*!< \brief Flag for whether to use the local normal as the flow direction for a pressure inlet. */
   su2double Linear_Solver_Error;   /*!< \brief Min error of the linear solver for the implicit formulation. */
   su2double Deform_Linear_Solver_Error;          /*!< \brief Min error of the linear solver for the implicit formulation. */
@@ -5212,6 +5214,19 @@ class CConfig {
    */
   su2double GetInc_Outlet_Damping(void) const { return Inc_Outlet_Damping; }
 
+  /*!
+   * \brief Check whether backflow prevention is enabled at incompressible outlets.
+   * \return True if the reversed normal velocity component should be zeroed on backflow faces.
+   */
+  bool GetInc_Outlet_BackflowPrevention(void) const { return Inc_Outlet_BackflowPrevention; }
+
+  /*!
+   * \brief Get the outer-iteration limit for backflow prevention at incompressible outlets.
+   * Backflow prevention is active while GetOuterIter() < this value.
+   * Returns ULONG_MAX (always active) when INC_OUTLET_BACKFLOW_PREVENTION_ITER is not set.
+   */
+  unsigned long GetInc_Outlet_BackflowPrevention_Iter(void) const { return Inc_Outlet_BackflowPrevention_Iter; }
+
   /*!
    * \brief Get the kind of mixing process for averaging quantities at the boundaries.
    * \return Kind of mixing process.

Common/include/parallelization/mpi_structure.cpp

@@ -101,12 +101,12 @@ void CBaseMPIWrapper::Error(const std::string& ErrorMsg, const std::string& Func
 
   /* Check if this rank must write the error message and do so. */
   if (Rank == MinRankError) {
-    std::cout << std::endl << std::endl;
-    std::cout << "Error in \"" << FunctionName << "\": " << std::endl;
-    std::cout << "-------------------------------------------------------------------------" << std::endl;
-    std::cout << ErrorMsg << std::endl;
-    std::cout << "------------------------------ Error Exit -------------------------------" << std::endl;
-    std::cout << std::endl << std::endl;
+    std::cerr << std::endl << std::endl;
+    std::cerr << "Error in \"" << FunctionName << "\": " << std::endl;
+    std::cerr << "-------------------------------------------------------------------------" << std::endl;
+    std::cerr << ErrorMsg << std::endl;
+    std::cerr << "------------------------------ Error Exit -------------------------------" << std::endl;
+    std::cerr << std::endl << std::endl;
   }
   Abort(currentComm, EXIT_FAILURE);
 }
@@ -130,12 +130,12 @@ void CBaseMPIWrapper::CopyData(const void* sendbuf, void* recvbuf, int size, Dat
 template <typename ScalarType>
 void CBaseMPIWrapper<ScalarType>::Error(const std::string& ErrorMsg, const std::string& FunctionName) {
   if (Rank == 0) {
-    std::cout << std::endl << std::endl;
-    std::cout << "Error in \"" << FunctionName << "\": " << std::endl;
-    std::cout << "-------------------------------------------------------------------------" << std::endl;
-    std::cout << ErrorMsg << std::endl;
-    std::cout << "------------------------------ Error Exit -------------------------------" << std::endl;
-    std::cout << std::endl << std::endl;
+    std::cerr << std::endl << std::endl;
+    std::cerr << "Error in \"" << FunctionName << "\": " << std::endl;
+    std::cerr << "-------------------------------------------------------------------------" << std::endl;
+    std::cerr << ErrorMsg << std::endl;
+    std::cerr << "------------------------------ Error Exit -------------------------------" << std::endl;
+    std::cerr << std::endl << std::endl;
   }
   Abort(currentComm, 0);
 }

Common/src/CConfig.cpp

@@ -1396,6 +1396,10 @@ void CConfig::SetConfig_Options() {
   addDoubleOption("INC_INLET_DAMPING", Inc_Inlet_Damping, 0.1);
   /*!\brief INC_OUTLET_DAMPING \n DESCRIPTION: Damping factor applied to the iterative updates to the pressure at a mass flow outlet in incompressible flow (0.1 by default). \ingroup Config*/
   addDoubleOption("INC_OUTLET_DAMPING", Inc_Outlet_Damping, 0.1);
+  /*!\brief INC_OUTLET_BACKFLOW_PREVENTION \n DESCRIPTION: Remove the reversed normal velocity component from the outlet ghost state when backflow is detected, preventing spurious inflow through the outlet. \ingroup Config*/
+  addBoolOption("INC_OUTLET_BACKFLOW_PREVENTION", Inc_Outlet_BackflowPrevention, false);
+  /*!\brief INC_OUTLET_BACKFLOW_PREVENTION_ITER \n DESCRIPTION: Number of outer iterations for which backflow prevention is active. After this iteration count the outlet BC reverts to pure Neumann, allowing physically real backflow to develop. Defaults to ULONG_MAX (always active). \ingroup Config*/
+  addUnsignedLongOption("INC_OUTLET_BACKFLOW_PREVENTION_ITER", Inc_Outlet_BackflowPrevention_Iter, numeric_limits<unsigned long>::max());
 
   /*--- Options related to the species solver. ---*/
 

SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp

@@ -189,6 +189,22 @@ class CSpeciesFlameletSolver final : public CSpeciesSolver {
   void BC_Inlet(CGeometry* geometry, CSolver** solver_container, CNumerics* conv_numerics, CNumerics* visc_numerics,
                 CConfig* config, unsigned short val_marker) override;
 
+  /*!
+   * \brief Impose the outlet boundary condition.
+   * When backflow prevention is enabled and reversed flow is detected at a face,
+   * the scalar ghost state is set to the inlet prescribed values (fresh mixture)
+   * instead of the Neumann interior values, preventing burned-gas scalars from
+   * re-entering the domain and pushing FGM look-ups outside the manifold.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] solver_container - Container vector with all the solutions.
+   * \param[in] conv_numerics - Description of the numerical method.
+   * \param[in] visc_numerics - Description of the numerical method.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] val_marker - Surface marker where the boundary condition is applied.
+   */
+  void BC_Outlet(CGeometry* geometry, CSolver** solver_container, CNumerics* conv_numerics, CNumerics* visc_numerics,
+                 CConfig* config, unsigned short val_marker) override;
+
   /*!
    * \brief Impose the (received) conjugate heat variables.
    * \param[in] geometry - Geometrical definition of the problem.

SU2_CFD/include/solvers/CSpeciesSolver.hpp

@@ -146,7 +146,7 @@ class CSpeciesSolver : public CScalarSolver<CSpeciesVariable> {
    * \param[in] val_marker - Surface marker where the boundary condition is applied.
    */
   void BC_Outlet(CGeometry* geometry, CSolver** solver_container, CNumerics* conv_numerics, CNumerics* visc_numerics,
-                 CConfig* config, unsigned short val_marker) final;
+                 CConfig* config, unsigned short val_marker) override;
 
   /*!
    * \brief Impose the isothermal wall Dirichlet boundary condition (value).

SU2_CFD/src/solvers/CIncEulerSolver.cpp

@@ -2639,7 +2639,10 @@ void CIncEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container,
   string Marker_Tag  = config->GetMarker_All_TagBound(val_marker);
 
   su2double Normal[MAXNDIM] = {0.0};
+  int nBackflow_loc = 0;
 
+  const bool backflow_prevention = config->GetInc_Outlet_BackflowPrevention() &&
+                                   (config->GetOuterIter() < config->GetInc_Outlet_BackflowPrevention_Iter());
   INC_OUTLET_TYPE Kind_Outlet = config->GetKind_Inc_Outlet(Marker_Tag);
 
   /*--- Loop over all the vertices on this boundary marker ---*/
@@ -2674,6 +2677,16 @@ void CIncEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container,
 
     P_domain = nodes->GetPressure(iPoint);
 
+    /*--- Compute the face area and normal velocity (positive = outflow, negative = backflow). ---*/
+
+    const su2double Area = GeometryToolbox::Norm(nDim, Normal);
+    const su2double Vn   = GeometryToolbox::DotProduct(nDim, &V_domain[prim_idx.Velocity()], Normal) / Area;
+
+    if (Vn < 0.0) {
+      SU2_OMP_ATOMIC
+      nBackflow_loc += 1;
+    }
+
     /*--- Compute a boundary value for the pressure depending on whether
      we are prescribing a back pressure or a mass flow target. ---*/
 
@@ -2744,6 +2757,26 @@ void CIncEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container,
 
     }
 
+    /*--- Backflow prevention via velocity reflection.
+     *
+     * When reverse flow is detected at the outlet (Vn < 0),
+     * we reflect the normal velocity component to prevent 
+     * backflow and stabilize the solution.
+     * 
+     * Factor of 2 is used to ensure that the reflected velocity
+     * is sufficient to counteract the backflow, by applying a
+     * restorative force that is proportional to the detected backflow velocity.
+     *
+     * Dynamic pressure is intentionally NOT applied here.  For
+     * variable-density / FGM cases the penalty is too small to
+     * overcome expansion-driven backflow and creates a thermodynamically
+     * inconsistent ghost state that pushes FGM controlling variables outside the manifold. ---*/
+
+    if (Vn < 0.0 && backflow_prevention) {
+      for (iDim = 0; iDim < nDim; iDim++)
+        V_outlet[iDim + prim_idx.Velocity()] -= 2.0 * Vn * Normal[iDim] / Area;
+    }
+
     /*--- Neumann condition for the temperature. ---*/
 
     V_outlet[prim_idx.Temperature()] = nodes->GetTemperature(iPoint);
@@ -2830,6 +2863,16 @@ void CIncEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container,
   }
   END_SU2_OMP_FOR
 
+  /*--- Print a warning if backflow was detected on this outlet marker. ---*/
+  BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
+    if (nBackflow_loc > 0) {
+      cout << "WARNING [Rank " << rank << "]: Backflow detected at outlet marker \""
+           << Marker_Tag << "\": " << static_cast<unsigned long>(nBackflow_loc)
+           << " face(s) have reversed normal velocity." << endl;
+    }
+  }
+  END_SU2_OMP_SAFE_GLOBAL_ACCESS
+
 }
 
 void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver_container, CConfig *config,

SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp

@@ -496,6 +496,142 @@ void CSpeciesFlameletSolver::BC_Inlet(CGeometry* geometry, CSolver** solver_cont
   CSpeciesSolver::BC_Inlet(geometry, solver_container, conv_numerics, visc_numerics, config, val_marker);
 }
 
+void CSpeciesFlameletSolver::BC_Outlet(CGeometry* geometry, CSolver** solver_container, CNumerics* conv_numerics,
+                                       CNumerics* visc_numerics, CConfig* config, unsigned short val_marker) {
+  SU2_ZONE_SCOPED
+
+  /*--- When backflow prevention is disabled, or the outer-iteration window has elapsed,
+   *    fall back to the base Neumann outlet so that physically real backflow (e.g.
+   *    recirculation zones at a converged solution) is not artificially suppressed. ---*/
+  if (!config->GetInc_Outlet_BackflowPrevention() ||
+      config->GetOuterIter() >= config->GetInc_Outlet_BackflowPrevention_Iter()) {
+    CSpeciesSolver::BC_Outlet(geometry, solver_container, conv_numerics, visc_numerics, config, val_marker);
+    return;
+  }
+
+  const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
+  const string Marker_Tag = config->GetMarker_All_TagBound(val_marker);
+
+  int nBackflow_loc    = 0;  /*--- Number of faces with backflow seen this call. ---*/
+  su2double Vn_min_loc = 0.0;  /*--- Most negative Vn seen this call. ---*/
+
+  // Compute a manifold-consistent enthalpy state for backflow faces.
+  // V_ref is the inlet bulk velocity; used to normalise the blending weight so that
+  // alpha -> 0 for marginal backflow and alpha -> 1 for strong backflow.
+  string Inlet_Tag;
+  bool found_inlet = false;
+  for (auto iMarker = 0u; iMarker < config->GetnMarker_CfgFile(); iMarker++) {
+    const string tag = config->GetMarker_CfgFile_TagBound(iMarker);
+    if (config->GetMarker_CfgFile_KindBC(tag) == INLET_FLOW) {
+      Inlet_Tag    = tag;
+      found_inlet  = true;
+      break;
+    }
+  }
+  if (!found_inlet) {
+    /*--- No inlet boundary defined — nothing to blend against; fall back to Neumann. ---*/
+    CSpeciesSolver::BC_Outlet(geometry, solver_container, conv_numerics, visc_numerics, config, val_marker);
+    return;
+  }
+
+  const su2double T_inlet = config->GetInletTtotal(Inlet_Tag);
+  const su2double V_ref   = config->GetInletPtotal(Inlet_Tag) / config->GetVelocity_Ref();
+
+  /*--- Determine thermodynamically consistent inlet enthalpy via Newton method. ---*/
+  su2double H_inlet;
+  GetEnthFromTemp(solver_container[FLOW_SOL]->GetFluidModel(),
+                  T_inlet, config->GetInlet_SpeciesVal(Inlet_Tag), &H_inlet);
+
+  SU2_OMP_FOR_STAT(OMP_MIN_SIZE)
+  for (auto iVertex = 0u; iVertex < geometry->nVertex[val_marker]; iVertex++) {
+    const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode();
+
+    if (!geometry->nodes->GetDomain(iPoint)) continue;
+
+    /*--- Outward face normal and normal velocity — needed for both strong and weak paths. ---*/
+    su2double Normal[MAXNDIM] = {0.0};
+    for (auto iDim = 0u; iDim < nDim; iDim++) Normal[iDim] = -geometry->vertex[val_marker][iVertex]->GetNormal(iDim);
+
+    auto V_domain = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint);
+    const su2double Area = GeometryToolbox::Norm(nDim, Normal);
+    const su2double Vn   = GeometryToolbox::DotProduct(nDim, &V_domain[prim_idx.Velocity()], Normal) / Area;
+
+    if (Vn < 0.0) {
+      SU2_OMP_ATOMIC
+      nBackflow_loc += 1;
+    }
+
+    /*--- Strong BC path: prescribe inlet scalars on backflow faces; otherwise copy
+     from the interior neighbour as usual. ---*/
+    if (config->GetMarker_StrongBC(Marker_Tag)) {
+      if (Vn < 0.0) {
+        nodes->SetSolution_Old(iPoint, Inlet_SpeciesVars[val_marker][iVertex]);
+      } else {
+        auto Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor();
+        nodes->SetSolution_Old(iPoint, nodes->GetSolution(Point_Normal));
+      }
+      LinSysRes.SetBlock_Zero(iPoint);
+      for (auto iVar = 0u; iVar < nVar; iVar++) Jacobian.DeleteValsRowi(iPoint, iVar);
+      continue;
+    }
+
+    /*--- Weak BC: flow ghost state set by CIncEulerSolver::BC_Outlet. ---*/
+    auto V_outlet = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, iVertex);
+
+    /*--- Choose scalar ghost state:
+     *   Outflow  (Vn >= 0): Neumann — interior values propagate outward (standard).
+     *   Backflow (Vn <  0): Dirichlet — Use blended approach to create thermodynamically consistent state. ---*/
+    su2double backflow_scalars[MAXNVAR];
+    const su2double* scalar_ghost;
+
+    if (Vn < 0.0) {
+      const su2double alpha      = min(-Vn / V_ref, 1.0);
+      const su2double H_interior = nodes->GetSolution(iPoint, I_ENTH);
+
+      /*--- Start from the interior solution, then blend the enthalpy toward the
+       *   thermodynamically consistent inlet value and clamp PV to zero. ---*/
+      for (auto iVar = 0u; iVar < nVar; iVar++)
+        backflow_scalars[iVar] = nodes->GetSolution(iPoint, iVar);
+
+      backflow_scalars[I_ENTH]    = (1.0 - alpha) * H_interior + alpha * H_inlet;
+      backflow_scalars[I_PROGVAR] = config->GetSpecies_Clipping_Min(I_PROGVAR);
+
+      scalar_ghost = backflow_scalars;
+    } else {
+      scalar_ghost = nodes->GetSolution(iPoint);
+    }
+
+    conv_numerics->SetPrimitive(V_domain, V_outlet);
+    conv_numerics->SetScalarVar(nodes->GetSolution(iPoint), scalar_ghost);
+    conv_numerics->SetNormal(Normal);
+
+    if (dynamic_grid)
+      conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint));
+
+    if (conv_numerics->GetBoundedScalar()) {
+      const su2double* velocity = &V_outlet[prim_idx.Velocity()];
+      const su2double density = solver_container[FLOW_SOL]->GetNodes()->GetDensity(iPoint);
+      conv_numerics->SetMassFlux(BoundedScalarBCFlux(iPoint, implicit, density, velocity, Normal));
+    }
+
+    auto residual = conv_numerics->ComputeResidual(config);
+    LinSysRes.AddBlock(iPoint, residual);
+    if (implicit) Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i);
+  }
+  END_SU2_OMP_FOR
+
+  /*--- Print backflow diagnostics once per BC call (no MPI reduction — safe). ---*/
+  BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
+    if (nBackflow_loc > 0 && rank == MASTER_NODE) {
+      cout << "  Flamelet BC_Outlet [" << Marker_Tag << "]: "
+           << nBackflow_loc
+           << " backflow face(s), Vn_min = " << Vn_min_loc << " m/s"
+           << "  --> applying inlet scalars (PV=0, H=H_inlet) on backflow faces." << endl;
+    }
+  }
+  END_SU2_OMP_SAFE_GLOBAL_ACCESS
+}
+
 void CSpeciesFlameletSolver::BC_Isothermal_Wall_Generic(CGeometry* geometry, CSolver** solver_container,
                                                         CNumerics* conv_numerics, CNumerics* visc_numerics,
                                                         CConfig* config, unsigned short val_marker, bool cht_mode) {

config_template.cfg

@@ -343,6 +343,16 @@ INC_OUTLET_TYPE= PRESSURE_OUTLET
 % Damping coefficient for iterative updates at mass flow outlets. (0.1 by default)
 INC_OUTLET_DAMPING= 0.1
 %
+% Activate backflow prevention at outlet markers for incompressible flows (NO, YES).
+% If YES, the normal velocity at outlet faces will be set to zero if backflow is detected.
+% This is only active for the first INC_OUTLET_BACKFLOW_PREVENTION_ITER iterations,
+% after which the BC will switch to the specified type in INC_OUTLET_TYPE.
+INC_OUTLET_BACKFLOW_PREVENTION= NO
+%
+% Number of iterations to apply backflow prevention at outlet markers for incompressible flows
+% If no value is specified, backflow prevention will be applied for the entire simulation.
+INC_OUTLET_BACKFLOW_PREVENTION_ITER= 1000
+%
 % Bulk Modulus for computing the Mach number
 BULK_MODULUS= 1.42E5
 % Epsilon^2 multipier in Beta calculation for incompressible preconditioner.