Discrete fracture network has significant effects on the thermal-hydraulic-mechanical-chemical coupling analysis in a highly fractured medium. Due to the complex configuration of interior fracture network, generation of a high resolving mesh is crucial in achieving accurate numerical results. An adaptive mesh refinement method is developed for generating adaptive meshes in a domain with discrete fracture networks. An optimized size function is embedded into Persson's mesh generator. The constrained Delaunay triangulation algorithm is integrated into the analogous force equilibrium method. An assessment method of mesh size quality is proposed. Fractures are pre-meshed with a user-specified length in advance. The whole domain is then adaptively meshed using a constrained Delaunay method in a locally relaxation manner. Varying meshes can be generated with fine meshes around fractures and coarse meshes in the domain. Two case studies are carried out and compared with those obtained from conventional methods. Results show that the present method is advantageous in generating a high resolving mesh. Finer meshes are adaptively generated around fractures and gradually changed to sparser meshes away from fractures. Smoothly and gradually varying meshes are generated with little sacrifice of the size and angle quality. The adaptive mesh refinement method can be used for a wide range of thermal, hydraulic, mechanical, chemical simulations and their numerical coupling system in a fractured medium. © 2014 Elsevier B.V.