Consider the equation for conservation of momentum in an inviscid flow, first in differential form:
$$\displaystyle \frac{\partial}{\partial t} (\rho\vec{u}) = -\nabla\cdot(\rho\vec{u}\vec{u}) – \nabla P $$
and then integrated over an arbitrary control volume $\Omega$, using the divergence theorem on the momentum term and the gradient theorem on the pressure term:
$$\displaystyle \frac{d}{dt} \int_\Omega \rho\vec{u}\,dx = -\int_{\partial\Omega} \rho\vec{u}(\vec{u}\cdot\hat{n})\,dS\, – \int_{\partial\Omega} P\hat{n}\, dS $$
CFD codes typically combine both surface integrals into an “inviscid flux” that leads to discretizations like:
$$\displaystyle \frac{d}{dt} (\rho\vec{u})_c = -\sum_f \frac{A_f}{V_c} (\rho\vec{u}\vec{u} + P\mathbb{I})_f\cdot\hat{n}_{f}$$
for a cell $c$ with volume $V_c$ and faces $f$, each of which as area $A_f$ and a unit normal vector $\hat{n}_f$ pointing outward from cell $c$. $\mathbb{I}$ is the identity matrix, so $P\mathbb{I}\hat{n}_f = P\hat{n}_f$. Conceptually, the term $\rho\vec{u}\vec{u} + P\mathbb{I}$ is evaluated at the face $f$. In practice that requires some approximation, though that won’t matter for present purposes.
Observe that this discretization folds the pressure, which began as a gradient term, into the divergence term arising from the momentum. This causes a problem when repurposing a CFD code designed for 2-D flow for axisymmetric flow because in the cylindrical coordinates used for the latter, the divergence:
$$\displaystyle \begin{align}\nabla\cdot\vec{v} &= \frac{1}{r}\frac{\partial}{\partial r}(rv_r) + \frac{\partial v_z}{\partial z} \\ &= \frac{\partial v_r}{\partial r} + \frac{v_r}{r} + \frac{\partial v_z}{\partial z} \end{align}$$
involves the $v_r/r$ term that has no analog in the gradient:
$$\displaystyle \nabla q = \frac{\partial q}{\partial r}\hat{r} + \frac{\partial q}{\partial z}\hat{z} $$
So starting from the discretization, and replacing the discrete divergence with the axisymmetric divergence the equation that is actually being solved by the code would be:
$$\displaystyle \begin{align}\frac{\partial}{\partial t} (\rho\vec{u}) &= – \frac{1}{r}\frac{\partial}{\partial r}(r (\rho\vec{u}u_r + P\hat{r}) ) – \frac{\partial}{\partial z}(\rho\vec{u}u_z + P\hat{z}) \\ &= -\frac{1}{r}\frac{\partial}{\partial r}(r (\rho\vec{u}u_r ) ) – \frac{P}{r}\hat{r} – \frac{\partial P}{\partial r}\hat{r} – \frac{\partial}{\partial z}(\rho\vec{u}u_z) – \frac{\partial P}{\partial z}\hat{z} \end{align}$$
However, the equation we want to solve is that obtained by replacing the divergence and gradient in the differential form at the top of this post with their axisymmetric forms, which gives:
$$\displaystyle \frac{\partial}{\partial t} (\rho\vec{u}) = – \frac{1}{r}\frac{\partial}{\partial r}(r (\rho\vec{u}u_r )) – \frac{\partial}{\partial z}(\rho\vec{u}u_z) – \frac{\partial P}{\partial r}\hat{r} – \frac{\partial P}{\partial z}\hat{z}$$
The difference is that the actual equation has a $-\frac{P}{r}\hat{r}$ term that should not be present, so it is necessary to add a $+\frac{P}{r}\hat{r}$ source term to the right-hand side to recover the equation that we actually want to solve. So the source of the source is the divergence.