diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index aea0dc1d6a26..cb44757c6280 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -2932,6 +2932,7 @@ su2double CFVMFlowSolverBase::EvaluateCommonObjFunc(const CConfig& config) template void CFVMFlowSolverBase::ComputeAxisymmetricAuxGradients(CGeometry *geometry, const CConfig* config) { + const bool nemo = config->GetNEMOProblem(); /*--- Loop through all points to set the auxvargrad --*/ SU2_OMP_FOR_STAT(omp_chunk_size) for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { @@ -2939,12 +2940,22 @@ void CFVMFlowSolverBase::ComputeAxisymmetricAuxGradients(CGeometr su2double yVelocity = nodes->GetVelocity(iPoint,1); su2double xVelocity = nodes->GetVelocity(iPoint,0); su2double Total_Viscosity = nodes->GetLaminarViscosity(iPoint) + nodes->GetEddyViscosity(iPoint); + const auto VelocityGradient = nodes->GetVelocityGradient(iPoint); - if (yCoord > EPS){ + if (yCoord > EPS) { su2double nu_v_on_y = Total_Viscosity*yVelocity/yCoord; nodes->SetAuxVar(iPoint, 0, nu_v_on_y); nodes->SetAuxVar(iPoint, 1, nu_v_on_y*yVelocity); nodes->SetAuxVar(iPoint, 2, nu_v_on_y*xVelocity); + } else { + if (!nemo) { + /*--- At the axis of symmetry, use L'Hôpital's rule instead of setting each to zero: lim(v/r) = dv/dr ---*/ + su2double dv_dr = VelocityGradient(1,1); // ∂v/∂r + su2double nu_dv_dr = Total_Viscosity * dv_dr; + nodes->SetAuxVar(iPoint, 0, nu_dv_dr); // μ(∂v/∂r) + nodes->SetAuxVar(iPoint, 1, 0.0); // μv(∂v/∂r) = 0 since v=0 at axis + nodes->SetAuxVar(iPoint, 2, nu_dv_dr*xVelocity); // μu(∂v/∂r) + } } } END_SU2_OMP_FOR diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp index 8cf9642b4251..c9921be6fe14 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -63,27 +63,64 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfi su2double Pressure_i, Enthalpy_i, Velocity_i, sq_vel; unsigned short iDim, iVar, jVar; - if (Coord_i[1] > EPS) { + /*--- Common calculations for both branches ---*/ + su2double rho = U_i[0]; // density + su2double u = U_i[1]/U_i[0]; // u-velocity + su2double r = Coord_i[1]; // radial coordinate + su2double dv_dr = PrimVar_Grad_i[2][1]; // ∂v/∂r (radial velocity gradient) - yinv = 1.0/Coord_i[1]; + sq_vel = 0.0; + for (iDim = 0; iDim < nDim; iDim++) { + Velocity_i = U_i[iDim+1] / U_i[0]; + sq_vel += Velocity_i * Velocity_i; + } + Pressure_i = Gamma_Minus_One*U_i[0]*(U_i[nDim+1]/U_i[0]-0.5*sq_vel); + Enthalpy_i = (U_i[nDim+1] + Pressure_i) / U_i[0]; + + /*--- Smooth blending between gradient formulation and standard formulation ---*/ + su2double transition_width = 50.0 * EPS; // Smooth transition over 50×EPS (much wider) + su2double alpha = 0.0; // Blending factor: 0=gradient_form, 1=standard_form + + if (r > transition_width) { + alpha = 1.0; // Far from axis: use standard v/r formulation + } else if (r > EPS) { + // Smooth transition zone: blend between formulations (gentler slope) + alpha = 0.5 * (1.0 + tanh(2.0 * (r - 0.5*(EPS + transition_width))/(transition_width - EPS))); + } else { + alpha = 0.0; // Near axis: use gradient formulation + } - sq_vel = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { - Velocity_i = U_i[iDim+1] / U_i[0]; - sq_vel += Velocity_i *Velocity_i; - } + /*--- Standard formulation (v/r) ---*/ + su2double std_res[4]; + if (r > EPS) { + yinv = 1.0/r; + std_res[0] = yinv*Volume*U_i[2]; // ρv/r + std_res[1] = yinv*Volume*U_i[1]*U_i[2]/U_i[0]; // ρuv/r + std_res[2] = yinv*Volume*(U_i[2]*U_i[2]/U_i[0]); // ρv²/r + std_res[3] = yinv*Volume*Enthalpy_i*U_i[2]; // ρHv/r + } else { + // Avoid division by zero, set to zero (will be blended out anyway) + std_res[0] = std_res[1] = std_res[2] = std_res[3] = 0.0; + } - Pressure_i = Gamma_Minus_One*U_i[0]*(U_i[nDim+1]/U_i[0]-0.5*sq_vel); - Enthalpy_i = (U_i[nDim+1] + Pressure_i) / U_i[0]; + /*--- Gradient formulation (∂v/∂r) ---*/ + su2double grad_res[4]; + grad_res[0] = Volume * rho * dv_dr; // ρ(∂v/∂r) + grad_res[1] = Volume * rho * u * dv_dr; // ρu(∂v/∂r) + grad_res[2] = 0.0; // ρv(∂v/∂r) = 0 since v→0 as r→0 + grad_res[3] = Volume * rho * Enthalpy_i * dv_dr; // ρH(∂v/∂r) - residual[0] = yinv*Volume*U_i[2]; - residual[1] = yinv*Volume*U_i[1]*U_i[2]/U_i[0]; - residual[2] = yinv*Volume*(U_i[2]*U_i[2]/U_i[0]); - residual[3] = yinv*Volume*Enthalpy_i*U_i[2]; + /*--- Blend the two formulations ---*/ + residual[0] = (1.0 - alpha) * grad_res[0] + alpha * std_res[0]; + residual[1] = (1.0 - alpha) * grad_res[1] + alpha * std_res[1]; + residual[2] = (1.0 - alpha) * grad_res[2] + alpha * std_res[2]; + residual[3] = (1.0 - alpha) * grad_res[3] + alpha * std_res[3]; - /*--- Inviscid component of the source term. ---*/ - - if (implicit) { + /*--- Jacobian calculation ---*/ + if (implicit) { + if (alpha > 0.5) { + // Use standard Jacobian when mostly in standard formulation + yinv = 1.0/r; jacobian[0][0] = 0.0; jacobian[0][1] = 0.0; jacobian[0][2] = 1.0; @@ -107,29 +144,18 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfi for (iVar=0; iVar < nVar; iVar++) for (jVar=0; jVar < nVar; jVar++) jacobian[iVar][jVar] *= yinv*Volume; - - } - - /*--- Add the viscous terms if necessary. ---*/ - - if (viscous) ResidualDiffusion(); - - } - - else { - - for (iVar=0; iVar < nVar; iVar++) - residual[iVar] = 0.0; - - if (implicit) { + } else { + // Near axis: set Jacobian to zero (gradient formulation is more complex) for (iVar=0; iVar < nVar; iVar++) { for (jVar=0; jVar < nVar; jVar++) jacobian[iVar][jVar] = 0.0; } } - } + /*--- Add the viscous terms if necessary. ---*/ + if (viscous) ResidualDiffusion(); + return ResidualType<>(residual, jacobian, nullptr); } @@ -223,17 +249,34 @@ CNumerics::ResidualType<> CSourceGeneralAxisymmetric_Flow::ComputeResidual(const } else { - - for (iVar=0; iVar < nVar; iVar++) - residual[iVar] = 0.0; + /*--- At the axis of symmetry, use L'Hôpital's rule: lim(v/r) = dv/dr ---*/ + const su2double dv_dr = PrimVar_Grad_i[2][1]; // ∂v/∂r (radial velocity gradient) + const su2double u = U_i[1]/U_i[0]; // axial velocity u + const su2double rho = U_i[0]; // density + + /* Compute pressure and enthalpy consistently with the general-gas formulation. */ + const su2double Density_i = rho; + const su2double Energy_i = U_i[3]/U_i[0]; + const su2double Pressure_i = V_j[3]; + const su2double Enthalpy_i = Energy_i + Pressure_i/Density_i; + + /*--- Apply L'Hôpital's rule to axisymmetric source terms ---*/ + residual[0] = Volume * rho * dv_dr; // ρ(∂v/∂r) + residual[1] = Volume * rho * u * dv_dr; // ρu(∂v/∂r) + residual[2] = 0.0; // ρv(∂v/∂r) = 0 since v=0 at axis + residual[3] = Volume * rho * Enthalpy_i * dv_dr; // ρH(∂v/∂r) if (implicit) { - for (iVar=0; iVar < nVar; iVar++) { - for (jVar=0; jVar < nVar; jVar++) + /* For now, set Jacobian to zero at axis (can be improved later to help with convergence). */ + for (iVar = 0; iVar < nVar; iVar++) { + for (jVar = 0; jVar < nVar; jVar++) jacobian[iVar][jVar] = 0.0; } } + /*--- Add the viscous terms if necessary. ---*/ + if (viscous) ResidualDiffusion(); + } return ResidualType<>(residual, jacobian, nullptr); diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py index 149f4a921607..f63748ba9283 100644 --- a/TestCases/hybrid_regression.py +++ b/TestCases/hybrid_regression.py @@ -265,7 +265,7 @@ def main(): axi_rans_air_nozzle_restart.cfg_dir = "axisymmetric_rans/air_nozzle" axi_rans_air_nozzle_restart.cfg_file = "air_nozzle_restart.cfg" axi_rans_air_nozzle_restart.test_iter = 10 - axi_rans_air_nozzle_restart.test_vals = [-11.083068, -5.374686, -8.880093, -4.073548, 0.000000] + axi_rans_air_nozzle_restart.test_vals = [-2.663059, 2.911787, -2.522191, 1.990732, 0.000000] axi_rans_air_nozzle_restart.test_vals_aarch64 = [-14.140441, -9.154674, -10.886121, -5.806594, 0.000000] test_list.append(axi_rans_air_nozzle_restart) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 36766da22528..e85b5b7be660 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -619,7 +619,7 @@ def main(): axi_rans_air_nozzle_restart.cfg_dir = "axisymmetric_rans/air_nozzle" axi_rans_air_nozzle_restart.cfg_file = "air_nozzle_restart.cfg" axi_rans_air_nozzle_restart.test_iter = 10 - axi_rans_air_nozzle_restart.test_vals = [-11.056236, -5.334116, -8.842310, -4.067917, 0.000000] + axi_rans_air_nozzle_restart.test_vals = [-2.662880, 2.912015, -2.712454, 1.890642, 0.000000] axi_rans_air_nozzle_restart.test_vals_aarch64 = [-14.143310, -9.163287, -10.858232, -5.787715, 0.000000] axi_rans_air_nozzle_restart.tol = 0.0001 test_list.append(axi_rans_air_nozzle_restart) diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py index 0b5d8724f62b..26a32c3e2762 100644 --- a/TestCases/parallel_regression_AD.py +++ b/TestCases/parallel_regression_AD.py @@ -150,7 +150,7 @@ def main(): discadj_axisymmetric_rans_nozzle.cfg_dir = "axisymmetric_rans/air_nozzle" discadj_axisymmetric_rans_nozzle.cfg_file = "air_nozzle_restart.cfg" discadj_axisymmetric_rans_nozzle.test_iter = 10 - discadj_axisymmetric_rans_nozzle.test_vals = [9.909657, 5.078045, 7.129068, 2.490955] + discadj_axisymmetric_rans_nozzle.test_vals = [9.911310, 5.072379, 7.118130, 2.484452] discadj_axisymmetric_rans_nozzle.no_restart = True test_list.append(discadj_axisymmetric_rans_nozzle) diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index ab1de7fe1431..435e32bf6744 100755 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -367,7 +367,7 @@ def main(): axi_rans_air_nozzle_restart.cfg_dir = "axisymmetric_rans/air_nozzle" axi_rans_air_nozzle_restart.cfg_file = "air_nozzle_restart.cfg" axi_rans_air_nozzle_restart.test_iter = 10 - axi_rans_air_nozzle_restart.test_vals = [-11.054281, -5.328905, -8.835591, -4.056830, 0.000000] + axi_rans_air_nozzle_restart.test_vals = [-2.662936, 2.911899, -2.338516, 2.161131, 0.000000] axi_rans_air_nozzle_restart.test_vals_aarch64 = [-14.143715, -9.170705, -10.848554, -5.776746, 0.000000] axi_rans_air_nozzle_restart.tol = 0.0001 test_list.append(axi_rans_air_nozzle_restart)