A novel semi-analytical (SA) model has been proposed to analyse mode I delamination in bonded composite structures including fasteners used for crack arrest. The SA model uses hierarchical shape functions describe the displacement field over the DCB specimen. The specimen is divided in three domains with each a set of hierarchical shape functions to accurately describe the displacement field around the crack tip, where the middle domain moves along with the crack tip. The Principle of Minimal Potential Energy combined with the Rayleigh-Ritz Method are used to determine the set of equilibrium equations. The total strain energy of the virtual DCB specimen is defined in terms of; the bending strain energy of the laminate, the strain energy of the bonding area between the sub-laminates modelled by an elastic foundation, and the fastener strain energy. The hole and the fastener position in the DCB specimen were accounted for in each of these terms. The overall displacement field obtained by solving the esultant set of equilibrium equations is used to define the specimen's compliance function, from which the Strain Energy Release Rate (SERR) can be determined for different load levels and crack lengths. An iterative scheme based on the Newton Raphson method is used to compute the crack lengths associated with different load levels during the delamination propagation. The results of the SA model in terms of load-displacement curves were compared with Cohesive Zone based Nonlinear FE predictions and experimental results. A very good agreement between results obtained using the SA model, FE model and experimental results was found. The proposed SA model allows large parametric studies at very low computational cost, being a powerful tool for preliminary analysis, design, and optimization of bolted-bonded composite parts.