In poorly developed fractured rocks, the contribution of individual fracture on rock conductivity should be considered. However, due to the lack of data, a deterministic approach cannot be used. The conventional way to model discrete fractures is to use a Poisson process, with prescribed distribution, for fracture size and orientation. Recently, a stochastic approach, based on the idea that the elastic energy due to fractures follows a Boltzmann distribution, has been used to generate realizations of correlated fractures in two dimensions. The elastic energy function has been derived by applying the appropriate physical laws in an elastic medium. The resulting energy function has been used in the simulated annealing algorithm to generate the realizations of two dimensional fracture networks. The main contribution of this work is to extend this technique to 3D, and to better incorporate geological field observations. In 3D, the method has adjusted the orientation of fractures to three orthogonal sets. Moreover, we investigate the effects of boundary condition, fracture size distribution, and the anisotropy of the medium. We have observed that far field stress can control the orientation of fractures. As a result, this fracture modeling technique can be used to stochastically generate sub-seismic fractures.