The present work deals with the numerical resolution of 3D-1D coupled problems arising from geometrical model reduction. A domain decomposition approach is used to de-couple the problem on the 1D domains from the problem in the 3D bulk domain. The de-coupling is based on a multi-field domain decomposition method, with new variables introduced at the internal interfaces, whereas the global solution is recovered through the minimization of a functional, expressing the error in the matching of the solutions at the interfaces. One complexity in 3D-1D simulations is related to the necessity of refining the mesh of the 3D domain near the inclusions to correctly reproduce the irregular behaviour of the solution in these regions. This partially offsets the advantage of the geometrical reduction. Here the use of the Extended Finite Element method is proposed to overcome this issue. The discrete function space for the 3D solution is enriched with additional basis functions and thus avoiding mesh refinement. A quadrature rule is proposed for the numerical integration of the enrichment functions that have a log-type singularity around the inclusions. The method is applied to the simulation of blood flow in a tissue with an embedded vascular network.