diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index e44ad6a10a7..eed445193a5 100644 --- a/Common/include/CConfig.hpp +++ b/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. diff --git a/Common/include/parallelization/mpi_structure.cpp b/Common/include/parallelization/mpi_structure.cpp index 03152e709e4..9d9e2ac7c61 100644 --- a/Common/include/parallelization/mpi_structure.cpp +++ b/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 void CBaseMPIWrapper::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); } diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index c168ced0830..ff8d9aa7285 100644 --- a/Common/src/CConfig.cpp +++ b/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::max()); /*--- Options related to the species solver. ---*/ diff --git a/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp b/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp index 6b9c4517298..f4b85d2096d 100644 --- a/SU2_CFD/include/solvers/CSpeciesFlameletSolver.hpp +++ b/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. diff --git a/SU2_CFD/include/solvers/CSpeciesSolver.hpp b/SU2_CFD/include/solvers/CSpeciesSolver.hpp index 3a657d97c01..43ab3f5c537 100644 --- a/SU2_CFD/include/solvers/CSpeciesSolver.hpp +++ b/SU2_CFD/include/solvers/CSpeciesSolver.hpp @@ -146,7 +146,7 @@ class CSpeciesSolver : public CScalarSolver { * \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). diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index 04757eddef6..f3d9c4383f8 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/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(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, diff --git a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp index a6d2930a904..710a5981078 100644 --- a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp +++ b/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) { diff --git a/config_template.cfg b/config_template.cfg index 6f6d5ee9bcc..4bda172f890 100644 --- a/config_template.cfg +++ b/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.