Abstract: For a nonlinear Poisson-Nernst-Planck mathematical model, in order to improve the stability and efficiency of numerical solution process, the edge-averaged finite element discretization scheme is derived, and a coupled iterative algorithm for numerical solution is given. Under some mild assumptions, the stiffness matrix of edge-averaged finite element discrete scheme is an M-matrix, and the numerical solution is more stable than the standard finite element method. The numerical results show that the L² norm error convergence order of the edge-averaged finite element method is optimal, and the CPU time of the edge-averaged finite element method is about one third of that of the standard finite element method under the same degrees of freedom.