Can we extract contact normals in a scene that undergoes in rigid body simulation?
Don't know much about how rigid bodies work behind the scenes, and haven't really scripted rigid body sims before. 

But if the contact normals are the normals of the faces that were touched during collision, you might be able to get them with a little hack: set up dynamic paint to paint the faces that are touched by the colliding object, then use a script to find all the faces that were painted (via the vertex color map created by the dynamic paint simulation).

Here's a little script that finds these normals if you use a red color for painting the touched object (assumes that the painted object is the active object, with the dynamic paint modifier applied and the default vertex color map name was used for the dynamic paint sim):

import bpy, bmesh
from collections import defaultdict

def get_vcolors( obj ):
    ''' 
        Vert colors are arranged in an unintuitive way, according to loop vertices, go figure.
        Also, every vertex has as many vertex colors as the number of faces its included in.
        This function outputs a dictionary where every vertex index has a list of vertex colors 
        associated with that vertex.
    '''
    vcolor_dict = defaultdict(list)
    mesh        = obj.data
    color_layer = obj.data.vertex_colors['dp_paintmap']

    i = 0
    for poly in mesh.polygons:
        for idx in poly.loop_indices:
            loop  = mesh.loops[idx]
            color = color_layer.data[i].color
            vcolor_dict[loop.vertex_index].append(color)
            i += 1

    return vcolor_dict

o = bpy.context.object
vcolors = get_vcolors( o )

# Find all the verts that have a red value above 0.9
red_verts = []
for v in vcolors:
    for c in vcolors[v]:
        if c.r &gt; 0.5:
            red_verts.append( v )
            continue

# Find the faces that include a red vert and their normals
red_faces = []
for i,f in enumerate(o.data.polygons):
    for rv in red_verts:
        if rv in f.vertices:
            red_faces.append( i )

red_faces = list(set(red_faces) ) # Remove doubles
red_faces_normals = [ o.data.polygons[i].normal for i in red_faces ]


You can find the setup I used for the simulations here (blend file):
http://www.pasteall.org/blend/37568