Correction of symmetry plane implementation

[bigfooted]

Proposed Changes

The implementation of the symmetry plane is incomplete. We follow here the book of Blazek, Computational Fluid Dynamics: Principles and Applications . According to Blazek, (chapter 8.6) 4 conditions have to be met on a symmetry plane:

1. No flux across the boundary
2. $\overline n \cdot \nabla \phi = 0$ (gradient of scalars are zero)
3. $\overline n \cdot \nabla (\overline U \cdot \overline t) = 0$ (gradient of tangential velocity normal to boundary)
4. $\overline t \cdot \nabla (\overline U \cdot \overline n) = 0$ (gradient of normal velocity along the boundary)

Points 2-4 all deal with gradients and can be dealt with in the gradient computation, i.e. Green-Gauss or Least Squares. According to Blazek, 2 approaches can be followed. "One possibility is to construct the missing half of the control volume by mirroring the grid on the boundary. The fluxes and gradients are then evaluated like in the interior using reflected flow variables." This approach can be implemented in an easy way when computing the Green-Gauss gradients. In SU2, routines are already in place that deal with GG gradients on boundaries. Here, we just have to identify the symmetry planes and mirror the flux through the faces. Blazek continues: "The second methodology computes the fluxes for the halved control volume (but not accross the boundary). The components of the residual normal to the symmetry plane are then zeroed out. It is also necessary to correct normal vectors of those faces of the control volume, which touch the boundary. The modification consists of removing all components of the face vector, which are normal to the symmetry plane. The gradients also have to be corrected according to Eq. (8.40) [points 2,3,4 above]"

1. Modify Green-Gauss and Least-Squares gradient computation for symmetry planes -> This seems to work fine
2. Fix normal components of edges connected to symmetry plane -> This seems to work fine
3. Fix residual components of velocity on the symmetry plane -> This seems to work fine
4. Multiple symmetry planes, including Euler walls -> This seems to work fine
5. multigrid -> Works fine now that we use the agglomeration as in the paper of Diskin.
6. grid movement
7. -> Works fine with the new fix to the residuals (energy update)

Related Work

#1168

PR Checklist

Put an X by all that apply. You can fill this out after submitting the PR. If you have any questions, don't hesitate to ask! We want to help. These are a guide for you to know what the reviewers will be looking for in your contribution.

• [x] I am submitting my contribution to the develop branch.
• [x] My contribution generates no new compiler warnings (try with --warnlevel=3 when using meson).
• [x] My contribution is commented and consistent with SU2 style (https://su2code.github.io/docs_v7/Style-Guide/).
• [x] I used the pre-commit hook to prevent dirty commits and used pre-commit run --all to format old commits.
• [ ] I have added a test case that demonstrates my contribution, if necessary.
• [ ] I have updated appropriate documentation (Tutorials, Docs Page, config_template.cpp), if necessary.

[bigfooted]

So are things exploding? If you feed in a solution and look at the gradients after the first iteration does it look like the modifications are doing the job?

I am using a very coarse 2D channel flow with the top a wall and bottom a symmetry, left inlet, right outlet (laminar incompressible navier stokes). The line figure is the normal velocity components at the symmetry axis, which should be zero. The black line develop. The red line is only the Green-Gauss correction. The peak at 0.5 is due to the mixed-mesh problem that the GG method has. The blue line is the current PR.

[bigfooted]

The edge correction in PhysicalGeometry is correctly modifying the normals, all edges between ipoint and jpoint that lie on the symmetry plane are identified, and the normal of the dual edge for those points is then corrected:

...
iPoint 329 and jPoint 338 are in the symmetry plane 
normal=-0.0284601 0
...

But only activating this edge correction is immediately causing issues. I would expect that this correction would have as an effect that the fluxes through those edges would not have components normal to the symmetry plane anymore.

[pcarruscag]

So modifying the gradient + residual is an improvement and it's just the normals that cause issues? Or the change to the residual is also problematic?

[bigfooted]

So modifying the gradient + residual is an improvement and it's just the normals that cause issues? Or the change to the residual is also problematic?

The modification of the Green-Gauss gradients only is working correctly. If I take the normal velocity on the symmetry plane as a performance indicator, then I observe that there is a mean (or maximum) discrepancy $\Delta V$ for develop, which reduces to less than $0.5 \Delta V$ for the new GG. [edit: additionally, the convergence behavior does not deteriorate for the GG modification]

[bigfooted]

Below is a new testcase, laminar flow around a cylinder. In the line figure I show the normal velocity on the axis y=0, with the solid red line the results from develop. The dashed line is with only the new Green-Gauss correction. It always leads to an improvement. These 2 results were obtained with an undisturbed grid. I then distort the grid a bit (as seen in the contour plot) and rerun the GG correction and I get the green line. Still an improvement, but a bit more bumps on the line. I then correct the edge normals as done in CPhysicalGeometry. I then get the jagged line shown. The correction using the residuals in BC_Sym_Plane seems to be ill-defined because the final converged solution depends on the initial solution [so I do not show tests with it here, I kept the original bc_sym_plane].

[bigfooted]

Symmetry plane is correct now. Below the figure of the super coarse diverging channel flow. Top figure: corrected symmetry, Bottom figure: original. The symmetry plane is in the middle of the figure. The top figure now has exactly V=0 on the symmetry plane.

[bigfooted]

Here is a picture of the diverging channel flow with the pressure contours. This is not related to the current implementation, but notice that at the corner where the outlet and the wall meet, there is a small pressure peak.

[bigfooted]

Update: Everything seems to work correctly now.

[bigfooted]

Results for laminar flow around a cylinder. right-top figure shows normal velocity at the symmetry axis, which should be zero. old (develop) can show very high values for normal velocity. new implementation: normal velocity < 1e-8 right-bottom figure: normal velocity along a vertical line in front of the cylinder. The normal velocity should be zero at the symmetry plane. Convergence is the same as develop.

[bigfooted]

Flow around a sphere, taking a 5 degree slice with 2 symmetry planes. Top:new implementation. bottom: old (develop). This picture shows axial velocity, for both converged to (pressure)R=-10.

[bigfooted]

TO DO: add testcase with 2 symmetry planes (laminar flow around sphere, slice)