Computationally efficient simulation methods for wave energy converters (WECs) are useful in a variety of applications. The simulation task is particularly challenging when nonlinearities are present in the WEC model. Using a Fourier projection of the system inputs and variables, harmonic balance (HB) is a computationally-efficient method to solve for the steady-state motion of a non-linear system, preserving an accurate representation of the non-linear effects. In previous work, HB has been used for the simulation of WECs with one degree of freedom (DoF). Here, HB is presented for WEC systems with an arbitrary number of DoFs. A non-linear, 2-DoF model of the ISWEC wave energy device is used as an example of application. The HB implementation of the ISWEC model is described in detail. Through numerical applications, chosen in both regular and irregular waves, general features of the HB method are exemplified, in particular the exponential convergence rate to the actual mathematical solution, and the sensitivity, in some cases, to the starting point of the HB algoritm