How to voxelize and extract blender data?
I just read your question statements and comments.
Generally your approach is fine in my opinion. I try to supply some details here. Hope them help.

External or internal?

Luis B suggested in his comment:


  Have you tried using the Blender Python API to create an internal Python Script(as opposed to an external Python program)?


I think an external script may be more suitable in this case.
The conversion from a geometric representation(which is created by Blender) to computation representation(i.e. something like Eulerian grid) should be coupled with your simulation, not Blender.

For example, the computation grid may be staggered, or you want to use adaptive grid(such as an octree rather than a uniform grid), etc.
All of these requirements depend on your simulation's setting.

How to do the voxelization?

I assume you are using mesh geometry to represent elements and boundaries in your simulation.
Then in this case you need to traverse all meshes and decide which mesh the point belong to. 
For convex mesh you can shoot one ray and check the intersections with the mesh. Two intersections means a point is outside this mesh. Also you need to handle the special case which a point is just right on the surface of mesh.

This process will be slow. But I think it is ok since you only run it once for each simulation configuration.

If you do want to optimize its performance you can consider using some data structures like BVH(bounding volume hierarchy), kd tree or spatial hash.

Things get much simpler if you can use procedural objects. For example you model a sphere using distance field function like $x^2+y^2+z^2-1=0$ rather than polygons. Then it is easier to decide whether a point is in a geometry(evaluate the distance field function and a negative value indicates the point is inside).