In order that the phase-fractions sum to 1 it is necessary that the same
diffusivity is used for ALL phases in the implicitPhasePressure option. This is
guaranteed by the new alphaDByAf function which returns a single
surfaceScalarField diffusivity to be used when forming the Laplacian term in the
implicit phase-fraction diffusion correction equation in phaseSystemSolve.
The phase-pressure and turbulent dispersion interface terms are summed over all
phases and interfaces in alphaDByAf to form a single diffusivity.