Finite element modeling of charged species transport has enabled analysis, design and optimization of a diverse array of electrochemical and electrokinetic devices. These systems are represented by the Poisson-Nernst-Plank equations coupled with the Navier-Stokes equation, with a key quantity of interest being the current at the system boundaries. Accurately computing the current flux is challenging due to the small critical dimension of the boundary layers (small Debye layer) that require fine mesh resolution at the boundaries. We resolve this challenge by using the Dirichlet-to-Neumann transformation to weakly impose the Dirichlet conditions for the Poisson-Nernst-Plank equations. The results obtained with weakly imposed Dirichlet boundary conditions showed excellent agreement with those obtained when conventional boundary conditions with highly resolved mesh were employed. Furthermore, the calculated current flux showed faster mesh convergence using weakly imposed conditions compared to the conventionally imposed Dirichlet boundary conditions. This approach substantially reduces the computational cost of modeling electrochemical systems.