In this study, a three-dimensional, steady-state, non-isothermal, two-phase-flow computational fluid dynamics model is developed for polymer electrolyte fuel cells. It is implemented into an open-source library. All major transport processes, heat and mass transfer, electron/proton transfer, electrochemical reactions are taken into account and solved simultaneously. A two-phase Eulerian-Eulerian algorithm is applied to describe flow in the gas channels on the cathode side of the cell. The present model is validated experimentally, by means of polarization curves and local current density distributions. Fairly good agreement is obtained for a variety of current densities. It is observed that the local current density distributions change as a function of gas humidity.