The motions of ice floes in linear waves were simulated by coupling Discrete Element Method (DEM) and Computational Fluid Dynamics (CFD). The interactions between ice floes are investigated by DEM. The hydrodynamics of ice floes, mainly including the drag force and the buoyancy, are calculated by CFD. In the simulation, the ice floes are treated as discrete elements, and the contact forces between ice floes are determined by the Hertz-Mindlin (no-slip) contact model. The shape of ice floes is an approximate square composed of a number of spherical faces, which can reduce the computation cost. The waves are treated as First Order Airy wave, which is linear in nature and applied to small amplitude waves in shallow liquid depth ranges. The volume of fluid (VOF) method was adopted to capture the free surface. The simulation results are in agreement with the actual situation to a certain extent.