<p>Most existing research on two-phase electrohydrodynamic (EHD) flows is based on Taylor’s leaky dielectric model. A critical issue with this classical model is its failure to ensure thermodynamic consistency, which is essential for accurately representing the underlying physics. As a macroscopic method, the lattice Boltzmann (LB) method has shown promise in simulating two-phase EHD flows, but the thermodynamically consistent LB method for such complex flows is rarely reported or discussed. In this paper we first present a theoretical analysis of the two-phase EHD flows by using Onsager’s variational principle, which is an extension of Rayleigh’s principle of least energy dissipation and, naturally, guarantees thermodynamic consistency. It shows that the governing equations of the model include the hydrodynamic equations, the Cahn-Hilliard equation coupled with additional electrical effect, and the full Poisson-Nernst-Planck electrokinetic equations. After that, a coupled LB scheme is constructed for simulating two-phase EHD flows. In particular, in order to handle two-phase EHD flows with a relatively larger electric permittivity ratio, we also introduce a delicately designed discrete forcing term into the LB equation for electrostatic field. The current LB method is validated through simulations of electro-osmotic flow in a flat microchannel, droplet deformation in a uniform electric field both with and without shear flow, the equilibrium configuration of a droplet in electric fields, and droplet detachment during reversed electrowetting. The numerical results align well with analytical solutions and experimental data, demonstrating the feasibility of the proposed method.</p>