I don't know for sure, but the analytical formula may be limited to bending deformation only and may not include shear deformation. Shear deformation will be included in the simulation, so there will not be agreement for that reason.
But your question seems to be about obtaining a mesh independent result. You are comparing hex dominant with tetrahedron. You have not specified the element order: linear or quadratic, which is more important that hex vs tetrahedron. Program controlled is quadratic, so that may be assumed, but it is better to set what you want.
I suggest you create a shell model using a surface of the same dimensions. When meshed with quadratic quad elements, that will converge rapidly to a mesh independent value as the size of elements are reduced.
When you have a thin walled solid mesh undergoing bending, the most important parameter is the number of elements through the thickness. This is very difficult to do with tetrahedrons. One way to do that is to go into SpaceClaim and split the thickness in half, then use the Share button to merge the nodes at the center plane.
However, with a mesh method of Sweep with the source face being the large area, you can freely enter the number of divisions in the sweep direction, which is through the thickness. Make a series of elements through the thickness of 1, 2, 4, 8. Try that with linear and then repeat with quadratic elements. Do all of this with the element size found using the shell elements.
You can save some time if you use 1/4 symmetry.