<p>In this paper, we present a cohesive zone model for simulating complex crack propagation, including crack branching and coalescence, in multiphase composites. The adaptive element splitting scheme is utilized in relation with the domain integral to describe arbitrary crack propagation in the matrix as well as crack coalescence and branching phenomena occurring at the interface between the matrix and the inclusion, based on the cohesive surface element approach. To compute the directions of crack growth in the matrix and crack bifurcation at the interface, the maximum strain energy release rate criterion is utilized in conjunction with the virtual grid-based stress recovery scheme. The proposed numerical methodology is validated by solving the fracture process of matrix-inclusion composites for both uniaxial tension and three-point bending configurations. The computational results agree well with the previous numerical and experimental results, including the crack pattern and global load (or stress) versus displacement response, while complex crack evolution, i.e., crack branching and coalescence, is effectively represented using the proposed framework.</p>