Polycrystalline composites are used in many engineering fields, as they combine significant mechanical and chemical properties. Particularly, Ceramic Matrix Composites (CMCs), comprising ceramic particles or short/long fibers dispersed within a ceramic matrix, exhibit corrosion resistance, high fracture toughness and capability of preserving mechanical properties at extreme temperatures. However, their modeling is still a challenging task, due to the heterogeneous nature of the material, its randomness characteristic and possible activation of nonlinear mechanisms. In this work, a reliable and fast modeling approach to study the nonlinear mechanical response of this category of composites is proposed. This model exploits the advantages of the current Virtual Element Method and the classical Finite Element Method. Each crystal, also referred to as grain, is discretized with a single low order virtual element assuming constant strain interpolation over the element. The interaction between grains is simulated by means of nonlinear interface finite elements based on a damage-friction constitutive law. Therefore, the typical intergranular crack growth is numerically reproduced by avoiding refined finite element grain discretizations, with significant computational cost saving. Reliability of the developed strategy is proved through some numerical tests. Specifically, the response of Alumina/Zirconia composites, a paradigmatic example of CMC, is studied under tensile and shear loads. The investigation is performed at the level of a representative volume element of the material, whose size is properly defined on the basis of the results of a statistical procedure tailored for random composites.