Flat shading for non-planar polygons
Without knowing the internals of Blender, I would say it uses shared normals for the quads, splitting the normals only on edges between the quads, not between the triangles. So your idea of uploading a normal per quad is probably close to the result you are looking for.

If you load the mesh as quads, you are set. If you are using an old version of OpenGL where quads still exist, you may even be able to use solution 3).

If you load the mesh as triangles, you will have to figure out which triangles originated in quads. You can likely do this simply by enforcing split normals beyond a fixed angle between adjacent triangles.

I am not sure there is a simple dynamic solution, where you just render any mesh. You may be able to use the geometry shader to look at the adjacent vertices, but the easiest version is to smooth the normals once before upload.