The pressure Poisson equation depends of the divergence (mass continuity in incompressible flows). If the divergence is set zero before calculating the pressure, the pressure field will not exist.

You must use equations based on mass continuity to guaranty div u = 0 at all flow, even in interior.  If you need to set div u=0 at the boundaries, you need to do the same for the interior field.

Actually, the pressure Poisson equation is derived from the original set of governing equations where div u=0 is used. However, this Poisson equation does not guaranty the continuity equation if we do not set div u=0 at the boundaries. Mathematically, when div u=0 is imposed at the boundary, the equations are identical if the initial velocity field is divergence free. However, numerically, I found that, div u is not zero in the interior of the fluid, in some region div u>0, in some region, div u

You are correct, isn't the Poisson equation that guaranty the continuity equation. It's perfect normal that you got some regions that div u .neq. 0, it is normal even for staggered mesh. There are two important things there, the sum of div must be zero (just like you getting), and the div u .neq. 0 can't be very high. The higher div u .neq. 0 that you got are the higher velocities and pressure points where your solution is weaker (in numerical stability point of view), and if your simulation diverges, probably the motive is this point, this high div u .neq. 0.

Well, this is what I think about this, your situation is just fine.

But why do the original governing equations are singular, i.e the mass continuity equation+the momentum equation.

I know the reason now.