The problem of constraint preservation for discretizations of non-linear PDEs is addressed on the example of the hyperbolic Yang-Mills equations in temporal gauge. These equations preserve a nonlinear divergence field analogous to the electric charge for Maxwell's equations. We introduce and discuss several discretizations of these equations on finite element spaces of Lie algebra valued differential forms. Numerical experiments indicate that simply restricting the variational formulation to the Galerkin spaces yields substantial drift in the charge, contrary to the linear Maxwell case. We then propose a fully discrete method constrained with Lagrange multipliers, for which we prove discrete charge conservation and observe excellent energy conservation.