I made a mistake in my comment to your question, but here is the correct version:
Also, as you have found out yourself, it is the second global_matrix that needs to be transformed.
The correct formula is:
matrix_world = global_matrix * obj.matrix_world * global_matrix.transpose()
However this only works if global_matrix is just a rotation (which it is in your case). In general what you need to calculate is not the transpose of the matrix, but the inverse. In the case of rotation matrices the inverse is equal to the transpose.
So the general case looks like this:
matrix_world = global_matrix * obj.matrix_world * global_matrix.inverse()