Without opening your part here are my suggestions. Crank up the mesh density. As a rule of thumb you need at least 3 elements through the thickness to account for bending in first order elements even more with tets like cosmos uses. Ignore the values near the constraints for comparison even though they may be real. Beam theory does not account for stress concentrations near boundary conditions. I hope this helps.
Thanks for your reply and comments. When I refine the mesh around supports the deviation becomes worse. I'm wondering why this deviation can only be observed for round bars or tubes. With a solid square section there is no difference between using Solid or Beam mesh. How can I find the average maximum bending stress apart form those stress concentration areas? The results shouldn't be that different from beam flexural formula.