What is $V_\mathrm{eff}$?
The effective potential from S&T isn't the gravitational potential, $V$, that appears in places like $\vec{F}_g = -m \nabla V$. That effective potential comes from a particular way of writing down the radial part of the geodesic equation for a point mass in Schwarzschild spacetime,
$$\left(\frac{dr}{d\tau}\right)^2 = \tilde{E}^2 - V_\mathrm{eff},$$
where $\tau$ is the proper time along the geodesic path, $\tilde{E}$ and $\tilde{l}$ are the energy and angular momentum per unit mass of the orbiting object, finally S&T (and most authors) define the effective potential without the square root,
$$V_\mathrm{eff} = \left(1-\frac{2M}{r} \right)\left(1+\frac{\tilde{l}^2}{r^2} \right).$$
This effective potential has more going on than just the Newtonian-like gravitational potential. Expanding it out there are four terms:
$$V_\mathrm{eff} = 1 - \frac{2M}{r} + \frac{\tilde{l}^2}{r^2} - \frac{2M\tilde{l}^2}{r^3}$$
- The $1$ comes from the rest energy of the orbiting body. It's kind of like an "$\tilde{m}$" term, the mass per unit mass of the orbiting body. It might be worth remembering that the energy in the relativistic equation includes the rest energy, which does not directly come into play in the Newtonian version of the two body problem.
- The $\frac{2M}{r}$ is the regular Newtonian potential.
- The $\frac{\tilde{l}^2}{r^2}$ term is the effect of the conserved angular momentum. It is sometimes called the "centrifugal" part of the potential. It also appears in the Newtonian two body problem.
- The cross term $\frac{2M\tilde{l}^2}{r^3}$ is new to the GR case, and mixes gravitational and centrifugal effects.
This is exact for the case it describes, but the question is what do you want to do with the GR corrections?
What do we want to model?
The papers you point to want to model accretion disks. This is not a simple point mass moving along a geodesic. For accretion, one might worry about viscous forces within the accreting matter. In which case the matter doesn't follow a geodesic of the spacetime, so we would need to start over from a new differential equation. Someone might also worry about gravitational interactions between particles in the accretion disk. It would be way easier to model those interactions as Newtonian rather than trying to solve the N-body problem in GR. There could be other forces as well, like effects from magnetic fields.
There is one big issue that jumps out at me in trying to connect the two body problem to the accretion disc problem: the energy and angular momentum of a particle in the accretion disk aren't conserved. Even in strictly Newtonian accretion dynamics, particles must lose angular momentum somehow in order to fall in. So pulling the centrifugal, $\tilde{l}$ term into an effective potential (which is common in the Newtonian two body problem, too) isn't helpful for an accretion disk. If the method of shedding angular momentum is viscous forces, then the particle's energy isn't conserved either. The viscosity provides friction which converts mechanical energy to heat, removing it from the dynamics. The method of simplifying the two body differential equation in terms of conserved quantities doesn't work. In the accretion problem $\tilde{E}$ and $\tilde{l}$ are functions of time, not constants.
For this modeling task, it seems like we want to solve Newton's second law:
$$m \frac{d^2\vec{r}}{dt^2} = \vec{F}_{g,\mathrm{BH}} + \vec{F}_\mathrm{visc} + ...$$
The effective potential for two body geodesics isn't going to tell us $\vec{F}_{g,\mathrm{BH}}$. What Paczynski and Wiita wanted was a modification on $V = -\frac{2M}{r}$ that would tell them how the gravitational field of the central BH is modified by GR, without any of the orbit specific, centrifugal stuff. With that in hand they can approximate $\vec{F}_{g,\mathrm{BH}}$ on the disk and easily add in the viscous forces, or whatever else too.
Interestingly, Abramowicz (2009) showed that the Paczynski-Wiita potential, which seems to have come from pure intuition, is the correct psuedo-Newtonian potential. It can be obtained by isolating the gravitational part from the centrifugal part of the GR two-body effective potential.
