We have numerically investigated an initial data construction of the black hole universe model which consists of cubic lattice of identical rotating black holes. It is the solution of vacuum Einstein constraint equations. To realize our numerical construction of the model, a cubic domain having periodic boundary conditions on each sides is introduced. Our numerical solutions are two parameter families of solutions characterized by the ADM angular momentum of the black hole and the coordinate size of its cubic domain. Our results show that, as the separation of rotating black holes increases, the relation between the Hubble expansion rate and the mass density of black holes becomes similar to the case of the Einstein de Sitter universe although the rotation slows down their agreement compared to the case of non-rotating black hole universe. The behaviors of scale factors defined in various ways are quite different from those in the non-rotating case, but become indentical as the distance between neighboring black holes becomes large. It is also shown that the presence of rotation makes the geometry nearby cubic surfaces anisotropic with negative curvature. In particular, we point out that the relation between the Hubble expansion rate and the mass density of black holes is well described by that of the Kasner universe when the rotating black holes become close.