To predict complex crack propagation behavior in multiphase composite materials, this study introduces a virtual element method (VEM)-based cohesive fracture model. Crack propagation, branching, and coalescence in multiphase composites are effectively captured using the adaptive element splitting scheme within the framework of extrinsic cohesive zone modeling [1]. Due to the flexibility in element shape within VEM, mesh modification processes, such as element splitting, are efficiently handled while maintaining high mesh quality. To determine the crack growth direction in the matrix and the crack branching direction at the interface, the maximum strain energy release rate criterion is employed, while the energy release rate is calculated using the J-integral as a line integral. For the constitutive relationship of the cohesive zone model, the modified PPR potential-based cohesive model is proposed, which provides more flexibility in the shape of the softening curve than the original one [2]. The proposed computational framework is validated through fracture simulations of multiphase composites, including a single fiber-matrix system and a cementitious composite with stiff circular inclusions. The computational results show good agreement with previous numerical and experimental results, including crack patterns and global load-displacement relationships. REFERENCES 1. Choi H., Ju M., Razakamandimby R. D.F.T., Park K., Cohesive zone modeling of crack propagation, branching, and coalescence in multiphase composites, Computational mechanics, 2025, in press 2. Park, K., Paulino, G.H., Roesler, J.R., A unified potential-based cohesive model of mixed-mode fracture, Journal of the mechanics and physics of solids, Vol. 57 (6), 891-908, 2009.