I'm having issues in handling the axisymmetric terms (like r*v/y) on nodes set on the symmetry line (y=0). Presently, I extrapolate the r*v/y term from the interior nodes. This provides a stable and satisfactory solution, however apon fine examination of the output, there are some oscillations near the symmetry line. I have also tried to offset the nodes from the symm. line such that there are no nodes at y=0. This eliminates any oscillations, but there are errors associated with the half cells below the symm. line which I don't know how to properly handle (if possible).

Just for reference, I'm modelling the transonic flow in a thrust chamber with an explicit FD Euler solver.

I could really use some more suggestions on how to handle the r*v/y term at y=0.

3D cylindrical polar coordinates and time accurate methods typical of LES codes have substantial problems with centrelines. There are a few papers on the subject and I recall a recent one by Lele in JCP (I think) but do not have the full reference with me. Personally, I have achieved the maintenance of wiggle free convection of smooth profiles by replacing the numerical metrics with their analytical equivalents (i.e. adding small correction involving sines and cosines) but have not written it up. It is an obvious thing to do and I am fairly sure other experienced CFD groups will have done something similar.

2D is not a major problem but to avoid oddities one does need to ensure the numerical representation of the terms behave appropriately as r tends to zero. Again no references to hand but I recall it being discussed in PhD thesis in the 70s and 80s when many of the CFD numerical schemes where first implemented. If you have access to PhDs from this time look particularly at the ones implementing RST turbulence models which had a range of problems with centre lines sufficient to actually change the conserved variables in some instances.

since y->0 and v->0, you can use v/y = d v/d y to compute v/y. AS for rho, axisymmetry implies d rho/ dy =0

> since y->0 and v->0, you can use v/y = d v/d y to compute v/y. AS for rho, axisymmetry implies d rho/ dy =0

This is the correct type of approach for 2D where there are no serious problems.

For 3D the velocity at the centre line is not zero and the gradients in the radial and circumferential velocity components are enormous because of the strong curvature as the singularity is approached. Maintaining reasonable behaviour in this case is far from easy.

The radial component of the velocity on the centerline must be zero for axisymmetric flow (unless there are singular sources or sinks placed on the line). Axisymmetric flow is 2D (x-r) by definition (no variation in circumferential direction).

By 3D you are probably referring to real 3D flow on an axisymmetric grid. That's a different issue (that is more robustly solved by using a different grid topology).