We present the first three-dimensional, fully compressible gas-dynamics simulations in 4π geometry of He-shell flash convection with proton-rich fuel entrainment at the upper boundary. This work is motivated by the insufficiently understood observed consequences of the H-ingestion flash in post-asymptotic giant branch (post-AGB) stars (Sakurai's object) and metal-poor AGB stars. Our investigation is focused on the entrainment process at the top convection boundary and on the subsequent advection of H-rich material into deeper layers, and we therefore ignore the burning of the proton-rich fuel in this study. We find that for our deep convection zone, coherent convective motions of near global scale appear to dominate the flow. At the top boundary convective shear flows are stable against Kelvin-Helmholtz instabilities. However, such shear instabilities are induced by the boundary-layer separation in large-scale, opposing flows. This links the global nature of thick shell convection with the entrainment process. We establish the quantitative dependence of the entrainment rate on grid resolution. With our numerical technique, simulations with 10243 cells or more are required to reach a numerical fidelity appropriate for this problem. However, only the result from the 15363 simulation provides a clear indication that we approach convergence with regard to the entrainment rate. Our results demonstrate that our method, which is described in detail, can provide quantitative results related to entrainment and convective boundary mixing in deep stellar interior environments with very stiff convective boundaries. For the representative case we study in detail, we find an entrainment rate of 4.38 ± 1.48 × 10-13 M ⊙ s-1.