Axisymmetric near-axis fixes using L’Hôpital-consistent formulations

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -2985,12 +2985,20 @@ void CFVMFlowSolverBase<V, FlowRegime>::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){
       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 {
+      /*--- 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

SU2_CFD/src/numerics/flow/flow_sources.cpp

@@ -63,27 +63,65 @@
   su2double Pressure_i, Enthalpy_i, Velocity_i, sq_vel;
   unsigned short iDim, iVar, jVar;
 
-  if (Coord_i[1] > EPS) {
-
-    yinv = 1.0/Coord_i[1];
+  /*--- Common calculations for both branches ---*/
+  su2double rho = U_i[0];                     // density
+  su2double u = U_i[1]/U_i[0];                // u-velocity
+  su2double v = U_i[2]/U_i[0];                // v-velocity
+  su2double r = Coord_i[1];                   // radial coordinate
+  su2double dv_dr = PrimVar_Grad_i[2][1];     // ∂v/∂r (radial velocity gradient)
+
+  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 +145,18 @@
       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 +250,34 @@
   }
 
   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);