Since natural fractures are often non-equidimensional, the circular disc model still has great limitations. By contrast, the elliptical disc model is more applicable to representing natural fractures, especially for slender ones. This paper developed a universal elliptical disc (UED) model by incorporating the center point, size, and azimuth of fractures as variables. Specifically, with respect to the azimuth of elliptical fractures in three-dimensional (3D) space, we proposed a paradigm to construct its probability density function (PDF) by coupling the orientation and rotation angle of long axis based on three coordinate transformations. To illustrate the construction process of the PDF of the fracture azimuth, we took the orientation following the Fisher distribution and the rotation angle following Von Mises distribution as an example. A rock slope is used to show the use of the developed UED model, and the 3D DFNs for the slope rock mass are generated by Monte Carlo simulation. In addition, the DFNs for the rock mass are also generated based on the existing circular disc model and non-universal elliptical disc model. The comparison results from the three models clearly illustrate the superiority of the UED model over the existing circular and non-universal elliptical disc models.