This paper provides a variational approximation by finite difference energies of functionals related to variational models in fracture mechanics for linearly--elastic materials in then framework of the Griffith theory of brittle fracture. The weak formulation of such a problem leads to functionals of the type \[ \mu\int_\Omega|{\mathcal E}u(x)|^2 dx+ \frac{\lambda}{2}\int_\Omega|\text{div}\,u(x)|^2 dx+ \int_\Omega\Phi([u],\nu_u)d{\mathcal H}^{n-1}\tag{*} \] defined on the space \(SBD(\Omega)\) of integrable functions \(u\) whose symmetrized distributional derivative \(Eu\) is a bounded Radon measure with density \({\mathcal E}u\) with respect to the Lebesgue measure and with a singular part concentrated on an \((n-1)\)-dimensional set \(J_u\), on which it is possible to define a normal \(\nu_u\) in a weak sense and one-sided traces. Here as usual, \({\mathcal H}^{n-1}\) is the \((n-1)\)-dimensional Hausdorff measure. The authors construct both the discrete and continuous approximation of \((*)\) by means of techniques different from those used by \textit{A. Chambolle} [M2AN, Math. Model. Numer. Anal. 33, No.2, 261-288 (1999; Zbl 0947.65076)], \textit{A. Braides} and \textit{M.S. Gelli} [J. Convex Anal. 9, No.2, 363-399 (2002; Zbl 1031.49022)], and \textit{M. Gobbino} [Commun. Pure Appl. Math. 51, No.2, 197-228 (1998; Zbl 0888.49013)], where it is reduced to the 1--dimensional case by an integral-geometric approach. In the present paper, such an approach is not possible due to the presence of the divergence term. Both approximations are performed in dimension 2, while in the discrete scheme boundary value problems are also treated, and an extension to dimension 3 is given.