The performance prediction of two counter-rotating vertical axis hydrokinetic turbines is presented in this paper. The flow field is governed by the 3D time-dependent incompressible Navier–Stokes equations. The system of equations is discretized using the Arbitrary Lagrangian-Eulerian Variational Multi-scale formulation for turbulence modeling on moving domains. Sliding interfaces are used to handle the rotor-stator interactions. Weak enforcement of essential boundary conditions is used to relax the requirement of boundary layers resolution. A grid convergence study based on the evaluation of the grid convergence index for the computed torque and the time-averaged axial wake velocity is performed. The grid convergence study shows good convergence with grid refinement in our formulation. Experimental validation of the averaged torque is performed for two different flow conditions and different turbine rotational velocities with good agreement. The finest mesh resolution from the grid convergence study is used for the multiple turbines simulation. The simulation shows almost 25% deficiency in the averaged torque of the downstream turbine. A multi-domain method is introduced to predict the performance of the turbine array at a lower computational cost than the full array simulation. The results from the multi-domain method show good matching on the averaged torque with the full array simulation. The initial study on the effect of struts on the torque generation is presented. The results assure the robustness of the ALE-VMS formulation and how it can be used to simulate multiple turbines at full-scale and full geometric complexity.