su2code
Symmetry BC not working (perfectly)
Hi all,
I am simulating a 2D channel with two triangular bumps in a Mach = 2.0 inviscid flow. When slicing the solution along x-normal planes, I can see that the symmetry condition does not lead to zero-gradients in the normal direction at the euler wall. Only the pressure variable has a zero-gradient at the wall. The situations doesn't change if I try with Roe, AUSMPLUSUP or JST schemes, nor if I change the way gradients are computed. Also refining the mesh does not improve the solution. I tried running the same case reverting SU2 to the commit previous to the implementation of the new symmetry BC and it works fine. The euler walls are the top and the bottom ones.
Here is the general flow solution colored by the Mach number
Here the solution with the new symmetry BC
Here the solution with the old symmetry BC
I have also tried with a subsonic flow and the result is the same. Notice how before the triangular bump the gradients are zero at the euler walls (red line).
All the simulations have converged to a RMS of the residual density of $10^{-13}$.
Here you can find the mesh and the config file. https://polimi365-my.sharepoint.com/:f:/g/personal/10507725_polimi_it/Em_r77HJAJpOiceuw5YSeAMB2CDPDoI5RJT_4CTe9fj7Aw?e=zz2jBN
I think that the problem is that when setting the normal gradients to 0 in the code then we are unable to reconstruct the primitive variables at the mid-point of the wall normal edges.
Desktop (please complete the following information):
- OS: Ubuntu
- C++ compiler and version: gcc 9.4.0
- MPI implementation and version: 4.0.3
- SU2 Version: v8, commit 77fed34
It looks like the energy also has nonzero normal gradients at the wall... Gradients are corrected in the correctGradientsSymmetry routine, called from green-gauss or least-squares. We get the local normal and then make a call to correctGradient, where we remove the normal component. All solver variables are corrected here, so not sure what is going on here. Any ideas, @pcarruscag ?
If the gradient is set to 0 then it is not possible to correctly reconstruct the variable at the mid point of the edge. What if we also set the variable to be equal to the one of the normal neighbor (for scalar quantities)?
I am talking about the ComputeFlux function in numerics_simd/flow/convection/roe.hpp, where the primitives are reconstructed.
Maybe I got it. The problem seems to be the fact that not all of the euler walls are aligned. Indeed, when I try to simulate a channel flow with only the upper bumps, the bottom euler wall (Y = 0) has gradients equal to 0.
Maybe it has to do with normals and we get back to #2383 .
This may be due to some of the dissipation terms that were kept for stability
/*--- Compute the residual using an upwind scheme. ---*/
auto residual = conv_numerics->ComputeResidual(config);
/*--- We include an update of the continuity and energy here, this is important for stability since
* these fluxes include numerical diffusion. ---*/
for (auto iVar = 0u; iVar < nVar; iVar++) {
if (iVar < iVel || iVar >= iVel + nDim) LinSysRes(iPoint, iVar) += residual.residual[iVar];
}
I have already tried removing those but the gradients were still not 0. I think the problem might be with the computations of the normals.
We do not strongly enforce the gradients at the walls. I think the normal computations are fine, it is I think the more complex boundary layer after the bump that causes strong gradients close to the surface. It is then more difficult to weakly enforce the gradients to be zero when the initial gradients are large. I wouldn't directly modify the values at the nodes, this is bad for convergence.
What do you mean that you do not strongly enforce the gradients at the wall? Isn't the correctGradientsSymmetry function doing exactly that? Plus, it is an euler simulation, thus no boundary layer is present.
Maybe something strange happens when a symmetry boundary condition is not aligned with the principal axes. I say this because when I tried removing the bumps on the lower boundary then the symmetry condition worked just fine on that boundary. On the other hand, the top boundary was still affected by the problem.
The gradients are corrected at the level of the green-gauss routine, and that is fed to the Euler equations. So we have the correct gradients in a field called gradients, and then we modify the linear solver residuals based on a normal correction. But we do not directly force the gradients of the solution of the Euler equations to zero like what you suggest, where we force the first points in the domain normal to the wall to have the wall value. With boundary layer I mean the first cell(s) above the wall . Where the shock hits the wall and then reflects, you have large gradients in e.g. energy, which then have to modified. btw, do you have literature results for this case?
I see what you mean.
I've found this case in the Nishikawa's limiters paper (https://arc.aiaa.org/doi/10.2514/6.2022-1374) and he doesn't seem to have this problem (Figure 14b, page 34)
Thanks for the paper. In the paper, the inlet density is 1, we have 2, is this just a nondimensionalization? If I simply scale the density by 2 then the results looks pretty similar to the result in the paper. It is difficult to see if Nishikawa has nonzero normal gradients of density close to the wall.
I also have density equal to 1 at the inlet in my simulations. Our simulations are non-dimensionalized. Also re-scaling the colorbar and the density to align with the paper's results, as you said, the problem is barely visible. But my point is that I don't understand why this happens with the new implementation of the symmetry BC whereas it was not happening before. Plus, it seems wrong as the gradients at the boundaries are not zero.
Thank you for testing the symmetry branch. What do you get for this case if you simulate it with the diamond shape in the middle and flat symmetry on top or bottom? If that is hard to setup, maybe top and bottom periodicity could be used with a small CFL, running explicit if necessary.
Thank you for testing the symmetry branch. What do you get for this case if you simulate it with the diamond shape in the middle and flat symmetry on top or bottom? If that is hard to setup, maybe top and bottom periodicity could be used with a small CFL, running explicit if necessary.
Yes, I tried it and it works better.
However, I have noticed that I have used Euler for the upper and lower lines instead of symmetry. When I try to use symmetry on the boundaries before and after the diamond then I get some sort of wake behind the diamond shape, which is kind of present also when the diamond is in mid air although much finer. Changing Symmetry BCs to Periodic BCs the solution improves, but I have asymmetries on the rear of the diamond shape between upper an lower wall. Theoretically speaking, the simulations with periodic and symmetry BCs should be equivalent to the one with the mid-air diamond shape, but they are not.
Original mesh
Original mesh with symmetry BC before and after diamond (it's the one that took more to converge)
With diamond shape mid-air
Original mesh with symmetry BCs as Periodic BCs
All simulations are at deep convergence (I have also tried the central diamond with periodic BCs and solution is the same):
No changes seen if I try to converge it further.
I think some of those differences may be due to new rules that Nijso worked on to merge intersecting Euler walls, which effectively creates slightly different normals near the leading and trailing edges of the diamond depending on the setup. You see a little normal shock on the Sym-Euler-Sym configuration that is not there with the Euler-Euler-Euler.
I don't know what else we can do to make the solution more symmetric. If we use those degrees of freedom to put conditions on gradients of variables it turns into some kind of finite difference and we no longer control the flux that is going through those edges.
Contour plot of density for the finest mesh (dx=0.0025, coarse == original mesh=0.01)
At the interface of the diamonds there's a visible wiggle, similar to what we see at the wall. You also see in the 2d plot above that the density gradient at the interface of the diamonds increases with increasing mesh size. Does the density gradient at the wall have the same root cause as the density gradient at the diamond interface?
For y-momentum, at the symmetry line at x=0 we see that the radial velocity becomes zero with mesh refinement. At the Euler wall, the gradient has also converged to a constant value
first order (MUSCL=OFF) with crazy large mesh:
No wiggles at the diamond interfaces, so that seems to be a higher order disturbance.
first order (MUSCL=OFF) with crazy large mesh:
No wiggles at the diamond interfaces, so that seems to be a higher order disturbance.
The isolines are more curved on the symmetry plane behind the diamond than the other symmetry plane on the opposite side. This may be associated to the change in boundary normal due to the diamond shape.
First order and symmetry BCs also for Fluent? What happens with second order?
First order and symmetry BCs also for Fluent? What happens with second order?
This was with a coarse mesh and second order in Fluent, using 'symmetry' for the bottom and 'inviscid wall' for the top, but switching to symmetry does not change the results. Fluent uses ROE-FDS. I also tried AUSM, but results are very similar. Also third order is similar.
What if you use symmetry also for the top? Do you obtain the same "wake" that I got?
No, switching between euler walls and symmetry does not seem to do anything.
Then maybe there is a problem with the symmetry plane. Can you please share the cas file that you are using?
public link to the fluent testcase: https://tue.data.surfsara.nl/index.php/s/kvHgjHxLJfJJyjn