AbstractThe hybrid finite element‐peridynamic (FEM‐PD) models have been evidenced for their exceptional ability to address hydro‐mechanical coupled problems involving cracks. Nevertheless, the non‐local characteristics of the PD equations and the required inversion operations when solving fluid equations result in prohibitively high computational costs. In this paper, a fast explicit solution scheme for FEM‐PD models based on matrix operation is introduced, where the graphics processing units (GPUs) are used to accelerate the computation. An in‐house software is developed in MATLAB in both CPU and GPU versions. We first solve a problem related to pore pressure distribution in a single crack, demonstrating the accuracy of the proposed method by a comparison of FEM‐PD solutions with those of PD‐only models and analytical solutions. Subsequently, several examples are solved, including a one‐dimensional dynamic consolidation problem and the fluid‐driven hydraulic fracture propagation problems in both 2D and 3D cases, to comprehensively validate the effectiveness of the proposed methods in simulating deformation and fracture in saturated porous media under the influence of hydro‐mechanical coupling. In the presented numerical results, some typical strong dynamic phenomena, such as stepwise crack advancement, crack branching, and pressure oscillations, are observed. In addition, we compare the wall times of all the cases executed on both the GPU and CPU, highlighting the substantial acceleration performance of the GPU, particularly when tackling problems with a significant computational workload.