Multi-point axisymmetric boundary condition for Euler equations

Question

I'm solving 2D axisymmetrical Euler's equations in conservative form: $$ \frac{\partial U}{\partial t} + \frac{\partial F(U)}{\partial x} + \frac{\partial G(U)}{\partial r} = H(U) $$ where $$ U = \left( \begin{array}{c} r\rho \\ r\rho u \\ r\rho v \\ re\end{array} \right), \; F(U) = \left( \begin{array}{c} r\rho u \\ r(\rho u^2 + p) \\ r\rho uv \\ r(e+p)u\end{array} \right), \; G(U) = \left( \begin{array}{c} r\rho v \\ r\rho uv \\ r(\rho v^2 + p) \\ r(e+p)v\end{array} \right), \; H(U) = \left( \begin{array}{c} 0 \\ 0 \\ p \\ 0 \end{array} \right), $$

using finite-difference WENO5 method.

How to correctly impose discrete axisymmetric boundary condition at $r=0$? One paper regarding boundary condition (for full cylindrical coordinates) just mentions that symmetry condition should be used for axisymmetric flows, but without any details.

When I was using MacCormack's method (2nd order, for each grid point one more neighboring point in each direction is needed), the procedure was pretty simple (C-syntax):

//internal points
for (k = 1; k <= k_max - 1; k++)
{
    for (l = 1; l <= l_max - 1; l++)
    {
        R[k][l] = ... (calculated by MacCormack's method); 
        U[k][l] = ...;
        V[k][l] = ...;
        P[k][l] = ...;
    }
}
//boundary at r = 0
for (k = 0; k <= k_max; k++)
{
    R[k][0] = R[k][1]; 
    U[k][0] = U[k][1];
    V[k][0] = 0;
    P[k][0] = P[k][1];
}
//other boundaries...

where R,U,V,P are arrays for $\rho, u, v, p$, first array index k is for $x$, second l is for $r$ (uniform square grid) and l=0 correspond to $r=0$. This boundary condition seems to work well.

With WENO5, the problem is that for each grid point, three more points in each direction are used in the stencil, so specifying one point at l=0 isn't enough.

My current ideas are:

1. Explicitly set three near-axis points:

//internal points
for (k = 3; k <= k_max - 3; k++)
{
    for (l = 3; l <= l_max - 3; l++)
    {
        R[k][l] = ... (calculated by WENO5 method); 
        U[k][l] = ...;
        V[k][l] = ...;
        P[k][l] = ...;
    }
}
//boundary at r = 0
for (k = 0; k <= k_max; k++)
{
    for (l = 0; l <= 2; l++)
    {
        R[k][l] = R[k][l]; 
        U[k][l] = U[k][l];
        V[k][l] = 0;
        P[k][l] = P[k][l];
    }
}
//other boundaries...

2. Add three ghost points for $r<0$ (so that l=3 at $r=0$) and set them with mirrored values:

//internal points
for (k = 3; k <= k_max - 3; k++)
{
    for (l = 3; l <= l_max - 3; l++)
    {
        R[k][l] = ... (calculated by WENO5 method); 
        U[k][l] = ...;
        V[k][l] = ...;
        P[k][l] = ...;
    }
}
//boundary at r = 0
for (k = 0; k <= k_max; k++)
{
    /values at r < 0
    for (l = 0; l <= 2; l++)
    {
        R[k][l] = R[k][6 - l]; 
        U[k][l] = U[k][6 - l];
        V[k][l] = -V[k][6 - l];
        P[k][l] = P[k][6 - l];
    }
    //values at r = 0
    R[k][3] = R[k][4]; 
    U[k][3] = U[k][4];
    V[k][3] = 0;
    P[k][3] = P[k][4];
}
//other boundaries...

or mirror with reversed sign (since values at ghost points are multiplied by negative $r$ in equations):

    ...
    for (l = 0; l <= 2; l++)
    {
        R[k][l] = -R[k][6 - l]; 
        U[k][l] = -U[k][6 - l];
        V[k][l] = V[k][6 - l];
        P[k][l] = -P[k][6 - l];
    }
    ...

Methods 1 and 2 generate noticable artifacts near the symmetry axis.

3. Shift the grid by half-step so that no point is placed exactly at $r=0$ and use ghost-points. This method leads to computation instability, though.

What would be the correct way?

Answer 1

As you mentioned above, the symmetry axis is treated in the same way that a general symmetry plane setting the ghost point values to their “mirror image” counterparts (negating the velocity). However, you have to be careful because you are using the geometric flux form of the Euler equations. You are "redefining" your governing equations multiplying them by the radial coordinate $r$. Your variables and fluxes at the ghost nodes where $r<0$ are being affected by this change in sign of your radial coordinate and so you might be imposing unphysical negative densities in those nodes. You can find more details in this paper:

https://www.sciencedirect.com/science/article/pii/S004579301500119X

Also, I want to point out to you that should have $\alpha r^{\alpha-1}p$ instead of only $p$ in your pressure nozzling term on the right-hand side of your governing equations. Here is another paper that shows the derivation of the form of the Euler equations that you are currently using:

https://link.springer.com/article/10.1007/s00193-017-0784-y