<p>In this paper, the problem of delaminated doubly-curved composite shells is taken into consideration including variable radii of curvature. The mechanical model was built-up based on the method of four equivalent single layers using the first- and second-order shear deformation shell theories. The governing equations were derived through the principle of virtual work for the delaminated and intact parts. The differential quadrature method was applied to solve the equations for three different delamination scenarios. Besides, the geometry of the shells was varied, too. In this respect, the parabola and ellipse profiles were utilized. In the first stage the deflections and their convergence were taken into account and the comparison with finite element model results was brought to the stage. The second stage was dedicated to the calculation of the J-integral and its mode-II and mode-III components over the four delamination fronts of the embedded rectangular delamination. Based on the fracture mechanical analysis and the maxima of the J-integral distributions, the location of delamination initiation points and the order of initiation were determined. The maximum value of the total energy release rate of the shells with ellipse profile varies between 7.3 and 98.6 % of that of the parabola profile by the first-order shear deformation theory, while the former value is 12.1–100.4% by the second-order shear deformation theory. The results show that the shell geometry (profile curve) and the delamination position influences significantly the mode of fracture and the location the delamination is initiated at.</p>