Summary: We consider a finite element approximation of a phase field model for the evolution of voids by surface diffusion in an electrically conducting solid. The phase field equations are given by the nonlinear degenerate parabolic system \[ \gamma\frac{\partial u}{\partial t}-\nabla.(b(u)\nabla[w+\alpha\phi])=0,\qquad w=-\gamma\Delta u+\gamma^{-1}\Psi'(u),\qquad\nabla.(c(u)\nabla\phi)=0 \] subject to an initial condition \(u ^{0}(\cdot )\in [ - 1,1]\) on \(u\) and flux boundary conditions on all three equations. Here \(\gamma \in \mathbb R_{>0}, \alpha \in \mathbb R_{\geq 0}, \Psi \) is a non-smooth double well potential, and \(c(u):=1+u, b(u):=1 - u ^{2}\) are degenerate coefficients. On extending existing results for the simplified two dimensional phase field model, we show stability bounds for our approximation and prove convergence, and hence existence of a solution to this nonlinear degenerate parabolic system in three space dimensions. Furthermore, a new iterative scheme for solving the resulting nonlinear discrete system is introduced and some numerical experiments are presented.