In this paper, an accurate analytical subdomain model to calculate electromagnetic performances of bearingless flux-switching permanent magnet (BFSPM) machines is developed. The entire region of the machine is divided into six subdomains, including slots, air gap, magnets, and so on. According to the Neumann boundary conditions between the subdomains and the iron and the continuity conditions at the interface between the adjacent subdomains, the magnetic field distributions of subdomains are obtained by solving the Maxwell equations. In addition, a method to compensate the saturation effect is proposed, and the magnetic saturation of the stator and rotor is considered by the actual B-H characteristics. Based on the proposed analytical method, the air-gap flux density distributions of the machine with different iron material under linear and nonlinear conditions are predicted. Further, the electromagnetic performances including the cogging torque, electromagnetic torque, flux linkage, back-EMF, inductance, and suspension force are calculated. These analytical results are compared with those obtained by the Finite Element Analysis (FEA) to verify the accuracy of the proposed analytical method. INDEX TERMS Bearingless flux-switching machines, subdomain model, magnetic field distribution, magnetic saturation.