In this paper, bulk drag coefficient in rigid dense vegetation is investigated mainly by using a three dimensional numerical simulation model CADMAS-SURF/3D by incorporating Immersed Boundary Method to calculate flow around the vertical cylinder in the Cartesian grid. Large Eddy Simulation is also incorporated as a turbulence model. Firstly, validation of the developed model is conducted with a single cylinder in the flow field based on literature. All the results obtained here (Re=300, 3,900 and 8,000) show good agreement with the reference data in literature. After the validation, multiple cylinders are allotted in three different densities (S/D=2.8, 2.0, 1.4) in a numerical wave tank and numerical simulations are conducted to investigate bulk drag coefficient. The result shows that the ratio of bulk drag coefficient to drag coefficient, which represents a reduction, is not just a function of density but a function of parameter 2a/S, in which 2a is stroke of the motion and S is cylinder distance. 2a is less than S, the effect of the density is neglected because the wake does not reach the other cylinders even when the density is high. On the contrary, it might affect the ratio of bulk drag coefficient to drag coefficient when the stroke of the motion is larger than the cylinder distance even when the density is low. In general, the ratio of bulk drag coefficient to drag coefficient decreases when 2a/S increases.