It is not an easy thing to do actually and Blender can´t do it out of the box.

If you got a quad which triangles do not share the same surface normal, you got one problem.

You can calculate the average normal out of those two triangles, but what then?

you can either move one of the vertices that is off, or two, you can move all 4 to flatten it, but do you move all in the same direction, or do you average over their position.

Also changing the vertices of one quad, means you create a triangular fold in the adjecant quad, then you got to flatten this one.

In the end you´ll have a recursive algorithm, rushing over the whole mesh, flattening one quad, and then flattening the 8 quads next to it, which got wonky for flattening one, those 8 again, effect 8 each, also the one you initially flattened.

In the end you got an exponentional algorithm, running recursive towards a termination condition like, stop if all the quads are “almost” flat.

And on top of that you´ll have a completely changed geometry and your mesh will look nothing like what it did.