Discrete element method (DEM) has become the method of choice for modeling the mechanical response of materials with granular structure. The method consists of representing the material as a collection of particles that interact with each other by transmitting forces through contact between neighboring particles. In DEM, the particles are typically idealized as rigid and the contact forces and associated deformations are modeled through spring-damper systems between the contact particles with various degrees of sophistication and complexity including particle bonding and cohesion necessary for concrete materials. The deformations due to contact include normal displacements as well as tangential displacements. The tangential forces can cause rolling, sliding as well as torsion. For bonded particles, the normal forces can be tensile or compressive. The normal and tangential responses due to applied forces or displacements are controlled by a multitude of user-specified parameters such as the damping coefficients, normal and tangential stiffnesses, tensile and compressive strengths and geometric factors such as the size and shape of the particles, as well as the number and packing arrangement of the particles. This thesis consists of two parts. In the first part, a parametric study of two-particle and three-particle systems was conducted to understand the effect of various DEM parameters on the system response. The DEM component of the commercial finite element software package LS-Dyna was used for the numerical simulations. The parallel-bond contact model was used for capturing the interaction between the particles. The second part of the thesis consists of modeling the uniaxial compressive response of an unreinforced concrete cylinder using a dense packing of spherical particles coupled by parallel-bonds. Guided by the results of the first part, the sensitivity of the predicted results for concrete response to particle size and influence of various parallel-bond parameters on the bulk response was studied to produce a calibrated model that shows to be capable of producing realistic uni-axial compressive behavior as observed in experiments. The results show that the maximum compressive stress in the concrete cylinder largely depends on the strength of parallel bonds and initial particle arrangement or sample void ratio. The axial strain to failure, as predicted by the DEM model considered, is largely dependent on the bond modulus and to a lesser extent on the particle size. Unlike other numerical techniques, particle size does not act as a free parameter that controls resolution; however, for uni-axial compression problem, it affects the overall damage process.