Degrading convergence when using convective BC, but not with isothermal/adiabatic walls.

[ChristianBauerEng]

Describe the bug
In my transient, compressible simulation the achievable residuals strongly depend on the type of boundary condition that is specified for a wall:

MARKER_ISOTHERMAL / MARKER_HEATFLUX (h=0):
rms_density < 1e-8
rms_energy < 1e-3

MARKER_HEATTRANSFER
rms_density > 1e-5
rms_energy > 1e-1

Here the development of the overall residuals can be seen. (other residuals look similar).

The temperature development along one of the walls can be seen here and is what I'd expect.

A more detailed analysis of local residuals shows, that the energy residuals are the largest directly at the convective wall:

I've tried various numerical schemes, limiters, linear solvers and gradient reconstruction methods, but none of these settings seem to improve convergence with convective BCs.

I am not sure if this is an actual bug or simply a limitation due to the way the convective BC is implemented. So far the solution looks realistic, but since my case is very sensitive to cumulative heat losses this may possibly have an influence on the final outcome.

Case summary

• Solver: RANS, transient
• Domain: 2D Axisymmetric
• Convection: SLAU2, NTS_DUCROS
• Gradients: GREEN_GAUSS
• Turbulence: SST V2003m, Sustaining
• Slope Limiter: Venkat_Wang
• Time: 2nd Order

Full setup, mesh and outputs can be found here.

Desktop (please complete the following information):

• OS: SLES 15.1
• C++ compiler and version: Intel 2021.4
• MPI implementation and version: intelMPI 2019
• SU2 Version: current master branch

[pcarruscag]

You could try turning off or simplifying the contributions to the Jacobian for that type of boundary, see CNSSolver.ccp around line 550. For example, leaving only the Jacobian_i[nDim+1][nDim+1] or none.

[ChristianBauerEng]

I believe this is the section you're talking about:

SU2/SU2_CFD/src/solvers/CNSSolver.cpp

Lines 552 to 559 in 1b08506

	/*--- Total specific energy: e=c_v*T+1/2*v^2 => T=1/c_v(rho*e/rho - 1/2||rho v||^2/rho^2).
	Together with cv=R/(gamma-1) the following Jacobian contributions for the energy equation can be derived. ---*/
	Jacobian_i[nDim+1][0] += Transfer_Coefficient * dTdrho * Area;
	for (unsigned short iDim = 0; iDim < nDim; iDim++)
	Jacobian_i[nDim+1][iDim+1] -= Transfer_Coefficient * dTdrhoe * nodes->GetVelocity(iPoint, iDim) * Area;
	Jacobian_i[nDim+1][nDim+1] += Transfer_Coefficient * dTdrhoe * Area;

Without understanding the structure of the Jacobian (and the overall solver) it's a bit hard for me to understand what's going on here. But if these contributions cause this issue than it should not occur when a Transfer_coefficient of 0.0 is specified, right? I tried that (by setting the BC in the config file accordingly), but the convergence issues remain.

[pcarruscag]

If you restart the heat transfer case (with h=0) from the 0 heat flux case, are the residuals the same or is there a jump?

[ChristianBauerEng]

I did some more testing and can confirm, that also when relaxing the convergence criterion the solution does not change. Therefore I'm positive that I can get good results also with convective BCs.

Nevertheless I tried restarting a simulation with convective BC and h=0 from an isothermal simulation and it diverged immediately after the first timestep. The same thing happened when restarting from an adiabatic simulation:

------------------------------ Begin Solver -----------------------------

Simulation Run using the Single-zone Driver
The simulation will run for 2000 time steps.
+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|   Time_Iter|  Inner_Iter|    Cur_Time|   Time(sec)|Linear_Solve|     Max CFL|    rms[Rho]|   rms[RhoU]|   rms[RhoV]|   rms[RhoE]|      rms[k]|      rms[w]|   Avg_Press|Avg_Press(NS|Avg_Press(NS|Avg_Press(NS|Avg_Press(NS|HF(NS_Resona|HF(NS_Resona|    tavg[HF]|tavg[Avg_Pre|tavg[Avg_Tem|
+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|        2000|           0|  5.0000e-08|  4.2625e-02|           3|  1.0000e+02|   -2.736211|   -0.228378|   -0.533349|    2.754171|    0.915535|    7.481489|  3.6909e+05|  0.0000e+00|  1.0000e+05|  1.0000e+05|  1.6909e+05|  4.5645e-11|  9.7518e-12|  0.0000e+00|  0.0000e+00|  0.0000e+00|
|        2000|           5|  5.0000e-08|  1.8793e-02|           9|  1.0000e+02|   -6.912682|   -4.751492|   -4.808573|   -1.347847|   -4.059524|    4.964360|  3.6910e+05|  0.0000e+00|  1.0000e+05|  1.0000e+05|  1.6910e+05|  1.3436e-12|  1.5626e-12|  0.0000e+00|  0.0000e+00|  0.0000e+00|
|        2000|          10|  5.0000e-08|  1.6455e-02|           9|  1.0000e+02|   -8.116638|   -6.253576|   -6.088205|   -2.591315|   -6.098708|    2.464595|  3.6910e+05|  0.0000e+00|  1.0000e+05|  1.0000e+05|  1.6910e+05| -1.7511e-12|  7.4888e-13|  0.0000e+00|  0.0000e+00|  0.0000e+00|
|        2000|          12|  5.0000e-08|  1.6035e-02|           9|  1.0000e+02|   -8.700803|   -6.595298|   -6.652705|   -3.158904|   -6.897678|    1.470020|  3.6910e+05|  0.0000e+00|  1.0000e+05|  1.0000e+05|  1.6910e+05| -3.5221e-12|  1.5395e-13|  0.0000e+00|  0.0000e+00|  0.0000e+00|
+-----------------------------------------------------------------------+
|        File Writing Summary       |              Filename             |
+-----------------------------------------------------------------------+
|SU2 binary restart                 |Flow_Restart_Convection_ZeroRestart_02000.dat|
|Paraview                           |Volume_Output_Convection_ZeroRestart_02000.vtu|
|CSV file                           |Surface_Output_Convection_ZeroRestart_02000.csv|
+-----------------------------------------------------------------------+
+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|   Time_Iter|  Inner_Iter|    Cur_Time|   Time(sec)|Linear_Solve|     Max CFL|    rms[Rho]|   rms[RhoU]|   rms[RhoV]|   rms[RhoE]|      rms[k]|      rms[w]|   Avg_Press|Avg_Press(NS|Avg_Press(NS|Avg_Press(NS|Avg_Press(NS|HF(NS_Resona|HF(NS_Resona|    tavg[HF]|tavg[Avg_Pre|tavg[Avg_Tem|
+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+


Error in "void CSolver::SetResidual_RMS(const CGeometry *, const CConfig *)":
-------------------------------------------------------------------------
SU2 has diverged (NaN detected).
------------------------------ Error Exit -------------------------------

[ChristianBauerEng]

It turns out using convective BCs drastically changes the global flow pattern.

Here you can see a comparison of pressures between a simulation with isothermal wall (top, 300 K) and convective BC (bottom, 300K, 20). Only the BCs in the conical part on the right side are different. Numerical settings and further BCs are identical.
https://github.com/su2code/SU2/assets/42602495/ba378660-ddd3-417c-a75e-f5ab54542697

Here you can see the difference in temperature.
https://github.com/su2code/SU2/assets/42602495/a6fcd09d-a065-46a5-b04a-60052fe35ed4

The difference means, that in one simulation (isothermal wall, top) I get strong, periodic oscillations, while in the other one (with convection) the oscillations are dampened out almost immediately.

Is there a SU2 unit test case for the validation of the convective boundary condition I can run?

[ChristianBauerEng]

I'm still not clear on how the convection BC is supposed to work.

If I understand correctly up to this point HEAT_FLUX and HEAT_TRANSFER are largely identical:

SU2/SU2_CFD/src/solvers/CNSSolver.cpp

Lines 476 to 479 in 1b08506

	else if (kind_boundary == HEAT_TRANSFER) {
	const su2double Twall = nodes->GetTemperature(iPoint);
	Wall_HeatFlux = Transfer_Coefficient * (Tinfinity - Twall);
	}

Essentially in both cases the heat flux is calculated and saved to Wall_HeatFlux

For both BCs this is taken into account via the residuals:

SU2/SU2_CFD/src/solvers/CNSSolver.cpp

Lines 491 to 496 in 1b08506

	/*--- Apply a weak boundary condition for the energy equation.
	Compute the residual due to the prescribed heat flux.
	The convective part will be zero if the grid is not moving. ---*/
	su2double Res_Conv = 0.0;
	su2double Res_Visc = Wall_HeatFlux * Area;

and the residual contributions is added:

SU2/SU2_CFD/src/solvers/CNSSolver.cpp

Lines 530 to 532 in 1b08506

	/*--- Convective and viscous contributions to the residual at the wall ---*/
	LinSysRes(iPoint, nDim+1) += Res_Conv - Res_Visc;

But I'm a bit lost why HEAT_TRANSFER is treated differently here:

SU2/SU2_CFD/src/solvers/CNSSolver.cpp

Lines 538 to 559 in 1b08506

	if (implicit) {
	if (kind_boundary == HEAT_TRANSFER){
	/*--- It is necessary to zero the jacobian entries of the energy equation. ---*/
	if (!dynamic_grid)
	for (auto iVar = 0u; iVar < nVar; ++iVar)
	Jacobian_i[nDim+1][iVar] = 0.0;
	const su2double oneOnRho = 1.0 / nodes->GetDensity(iPoint);
	const su2double oneOnCv = (Gamma - 1.0) / config->GetGas_ConstantND();
	const su2double Vel2 = nodes->GetVelocity2(iPoint);
	const su2double dTdrho = oneOnRho * ( -Tinfinity + oneOnCv * 0.5 * Vel2);
	const su2double dTdrhoe = oneOnCv * oneOnRho;
	/*--- Total specific energy: e=c_v*T+1/2*v^2 => T=1/c_v(rho*e/rho - 1/2||rho v||^2/rho^2).
	Together with cv=R/(gamma-1) the following Jacobian contributions for the energy equation can be derived. ---*/
	Jacobian_i[nDim+1][0] += Transfer_Coefficient * dTdrho * Area;
	for (unsigned short iDim = 0; iDim < nDim; iDim++)
	Jacobian_i[nDim+1][iDim+1] -= Transfer_Coefficient * dTdrhoe * nodes->GetVelocity(iPoint, iDim) * Area;
	Jacobian_i[nDim+1][nDim+1] += Transfer_Coefficient * dTdrhoe * Area;

Why is this special treatment necessary if HEAT_TRANSFER essentially is just a different way to calculate Wall_HeatFlux? If I disable the additional Jacobian contributions, shouldn't the solution still be accurate?

Regards,
Christian

[pcarruscag]

There is a contribution to the Jacobian because the heat flux depends on the temperature at the node. Does it work better if you remove that contribution?

[ChristianBauerEng]

On it, I've disabled the Jacobian contributions in my branch.

Initial tests look primising, but I'll need a longer-running test to validate it on my case. Will update once the case has finished.

[ChristianBauerEng]

After some testing I can report the following:

• With Jacobian contributions disabled convergence is greatly increased and on-par with the isothermal/heat flux configurations
• After around 12800 timesteps the solution diverges and SU2 crashes.
• For some reason the global flowfield is drastically different between convective and isothermal/adiabatic cases. As the convective case is somewhere between adiabatic and isothermal conditions this should not be the case.

Here you can see a comparison of pressure between convective case (top) and isothermal case (bottom). The only difference is the type of boundary of the conical walls.

The same comparison by temperature:

I'm a bit lost as to what's happaning here...