I think the first step should be to match the single phase results at the dilute limit.
I noticed that in the twoPhaseEulerFoam the near wall G is calculated from turbulent viscosity and set to zero if the y* is less than 11.6:
const scalarField& nuw = nutb.boundaryField()[patchi];
if (yPlus > 11.6)
While in the single phase epsilonWallFunction the G is calculated at all y* values and is calculated from the effective viscosity (nu+nut).
In my test case, the twoPhaseEulerFoam results improve with the effective viscosity based G-calculation and the Spalding turbulent viscosity wall function á la nutSpalartAllmarasWallFunction. The near wall behaviour matches quite well with single phase results and correlations. However, there is still a 1/3 of a difference in the centerline turbulent viscosity.