Abstract: To address the challenges posed by complex inclusion geometries, strongly coupled multiphysics 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 electric and magnetic Gauss equations are formulated within a unified framework based on a generalized modified complementary energy functional. Within each element, 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. This formulation enables the generalized internal fields to be determined directly and reduces the dependence of local-field accuracy on extensive mesh refinement near material interfaces. Taking BaTiO₃-CoFe₂O₄ particulate magnetoelectric composites as a representative material system, microstructural models with circular, elliptical, inclined elliptical, and multiple randomly distributed inclusions are established. Their coupled mechanical-electrical-magnetic responses are investigated and validated against reference finite element solutions obtained using refined meshes. Comparisons of contour distributions and path profiles show that the proposed method accurately captures the local distributions of stress, electric displacement, and magnetic flux density, including pronounced variations near matrix-inclusion interfaces. 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. Increasing the inclusion volume fraction significantly affects the effective mechanical, electrical, magnetic, and magnetoelectric responses, while variations in inclusion shape and orientation alter the directional dependence of the coupled response. In particular, elliptical inclusions exhibit a clear orientation dependence, leading to marked directional variations in the effective material properties.