A mixed time-domain/frequency-domain method is proposed for modelling dense wave energy converter (WEC) arrays with non-linear power take-off (PTO). The model is based on a harmonic balance method which describes the system response in the frequency domain, while evaluating the non-linear PTO force and solving the system equations of motion in the time domain. The non-linear PTO force is computed with Lagrange multipliers. In order to apply the proposed method for WEC array responses in real sea states, the time series is split into time windows and the simulation is carried out individually per window. The method is demonstrated by investigating the dynamics of the Ocean Grazer WEC array (OG-WEC) with an adaptable piston pumping system. The key parameters thought to possibly influence model accuracy, including the number of harmonic components, the length of the time window and overlay, are discussed. It is shown that the proposed model can significantly reduce the computational cost with an acceptable accuracy penalty.