<p>In this paper, we aim to design two energy-stable and efficient finite element schemes for simulating the ferrofluid flows based on the well-known Shliomis model. The model is a highly nonlinear, coupled, multi-physics system, consisting of the Navier–Stokes equations, magnetostatic equation, and magnetization field equation. We propose two reliable numerical algorithms with the following desired features: linearity and unconditional energy stability. Several key techniques are used to achieve the required features, including the auxiliary variable method, consistent terms method, prediction-correction method, and semi-implicit stabilization method. The first scheme is based on a hybrid continuous/discontinuous finite elements spatial approximation, and the second utilizes decoupled continuous finite element spatial discretization. We have rigorously demonstrated that the proposed schemes are unconditionally energy stable and carried out extensive numerical simulations to illustrate the accuracy and stability of the developed schemes, as well as some interesting controllable characteristics of the ferrofluid flows.</p>