The understanding of damage evolution and fracture propagation in rocks as a brittle material is important for many industrial applications; for instance to investigate mine safety, to analyze material failure in manufacturing and construction and to understand the behaviour of pre-existing fractures on hydraulic fracturing in oil and gas extractions. Currently, numerical studies for fracture propagation in brittle rock specimens with pre-existing fractures show strong mesh dependence. To address the issue of mesh independence, further numerical investigation and improvement of damage models under a fully 3-D assumption is critical.In this thesis, we develop numerical tools to investigate the effect of model parameters, including the heterogeneity of rock by allowing Young's modulus to vary spatially, and of the geometry of a pre-existing fracture on the stress-strain response and damage evolution in rocks. The fracture propagation model uses continuum damage mechanics where the elastic modulus reduces with damage as the rock is weakened through the formation of fractures. Two failure conditions are presented, a micromechanical model [1,2] using the unified strength theory [3,4] and a strain-based form of a modified von-Mises condition [5]. Various damage evolution laws are considered which can simulate a range of behaviour from brittle to quasi-brittle. Four different damage evolution laws are analyzed for a 3-D rectangular single flaw specimen to deliver the most realistic results. The quasi-static equilibrium equation is solved for various cases of brittle failure and damage evolution for specimens under loading using the Finite Element Method (FEM) with tetrahedron elements on unstructured meshes. In contrast to previous implementations based on rectangular grids, the use of unstructured meshes provides higher geometrical flexibility and allows a more accurate way of including geological features, such as flaws, through locally adapted (FEM) meshes.An important factor affecting the mechanical behaviour during the failure process is the heterogeneity of the rock. Heterogeneity is simulated by sampling the initial local Young's modulus from a Weibull distribution over a cubic grid. The values are then interpolated to the computational (unstructured FEM) mesh. This approach introduces an additional length scale for rock heterogeneity represented by the cell size in the sampling grid that is independent of FEM mesh size. This approach is different from previous implementations of rock heterogeneity [6][7][8] where the length scale is determined by the mesh resolution.A well-known problem which arises in "local" damage models is that they suffer from mesh dependency [2] and strain localization [9]. In this work, we apply the non-local implicit gradient damage formulation of Peerlings et al.[5] to address the issue of mesh dependency. The basic concept of the non-local implicit gradient model is that the strain at a given point depends not only on the strain at that point but also on the nearby strain field. The lo...