A finite element method, which directly uses the node electric scalar potentials and node charge densities as variables, is used to predict the electric field and charge density distributions under polarity reversal (PR) voltage. The penalty method is adopted to impose the Dirichlet boundary values and the Crank-Nicolson (C-N) algorithm is then used to solve the transient equation. The surface node charge densities on the perfectly conducting boundaries are adopted to accurately obtain the normal electric field strengths. The method is tested by a two layer coaxial model and applied to analyze the linear and nonlinear transient electric fields and boundary charge densities of a ±500 kV converter transformer under PR voltage.