This paper contributes a scheme for two-dimensional (2D) numerical simulation of hydraulic fracturing processes in geological materials, which consists of a numerical coupling between the finite element method (FEM) and the lattice Boltzmann (LB) model. The Park-Paulino-Roesler potential-based cohesive zone model (PPR) is used to simulate the fracture propagation by means of interface finite elements, such that cohesive forces act on the fracture surface, capturing the softening process. The PPR model is used because it is a generalized fracture model that can represent the fracture process for mode I, mode II and mixed mode I-II, and can be applied to various materials, including heterogeneous materials, such as rock. The FEM and LB are coupled in an iterative process. The paper describes implementation details including procedures for coupling both methods. Examples of hydraulic fracturing process, modeled with the proposed FEM-LB coupling demonstrate the potential of this numerical procedure to model hydraulic fracturing processes in geomaterials of complex geometries.