Hi,
I was checking out COMSOL (also implemented in FEM), and they refer to an article in their documentation that provides useful insights for implementing the k-epsilon model.
Additionally, there’s another good article titled “On the implementation of the κ-ε turbulence model in incompressible flow solvers based on a finite element discretization” that you can find on ResearchGate here.
Another resource worth checking out is the paper “A FEniCS-Based Programming Framework for Modeling Turbulent Flow by the Reynolds-Averaged Navier-Stokes Equations” available here. They implemented the equations directly in FEniCS, which could be relevant to your work. see this link for thier implementation in fenics. (CBC.PDESys in Launchpad)
Hope this helps!