Abstract: To address the challenges posed by complex inclusion geometries, difficult solution of coupled multi-field responses, and the high computational cost of conventional finite element methods for particulate magnetoelectric composites with randomly distributed inclusions, a multiphysics-coupled Voronoi cell finite element method (MC-VCFEM) is developed for microscale field analysis and macroscopic effective-property prediction. In the proposed method, Voronoi tessellation is used to partition the inclusions and their neighboring matrix regions, while the mechanical equilibrium equation and the Gauss’s laws for electricity and magnetism are formulated within a unified framework based on a generalized modified complementary energy functional. Within each element, self-equilibrated stress, electric displacement, and magnetic flux density fields satisfying the corresponding governing equations are independently constructed, while displacement, electric potential, and magnetic potential are interpolated along the element boundaries to enforce interelement compatibility. Taking BaTiO₃-CoFe₂O₄ particulate magnetoelectric composites as a representative material system, microstructural models with circular inclusions, elliptical inclusions, elliptical inclusions with different inclination angles, and multiple randomly distributed inclusions are established. Their local coupled mechanical-electrical-magnetic responses are investigated and validated against reference finite element solutions obtained using refined meshes. On this basis, the effects of inclusion volume fraction, shape, and orientation on the effective elastic coefficients, dielectric permittivities, magnetic permeabilities, and magnetoelectric coupling coefficients are examined. The results show that the proposed method significantly reduces the number of elements and total degrees of freedom while maintaining close agreement with conventional finite element solutions. Inclusion geometry has a pronounced influence on local field concentrations and macroscopic effective properties. In particular, changes in the orientation of elliptical inclusions further cause directional differences in the effective material parameters.