Okay so I have imposed limitations on variables that make the source term non-physical. The limitations constrain the force term similar to this graph I developed in MATLAB.
The simulation diverged once again. I changed the gradient volume fraction to the magnitude of the gradient... I am not sure about this.
Also, do you think that the gradient of the volume fraction is respective of dimension?