Granular materials, which can range from nanometers to centimeters, are widely encountered as particles or powders in nature and in engineered systems. Some examples of granular materials include sands, gravels, biomass feedstocks, rice grains, pharmaceutical pills, glass beads, and lunar or martian regolith simulants. The attributes of and the interactions between constituent particles dictate the bulk mechanical properties and behavior of granular materials. The discrete element method (DEM), a particle-based modeling technique, has become a useful numerical tool for modeling the bulk behavior (such as mechanical, flow, and breakage behavior) of granular materials, where interactions between the granular materials' individual constituent particles are explicitly modeled. Particle attributes such as size, shape, and flexibility play a fundamental role in affecting the bulk mechanical properties and behavior of granular materials, and the DEM should be able to effectively capture those attributes. Therefore, the main objective of this dissertation is to develop discrete element models for stiff and flexible granular materials with increased complexity in particle shape representation and breakage behavior and to apply the DEM models to laboratory and industrial scale processes. Addressing the main objective, four studies have been carried out. The first study focuses on experimentally characterizing the LHS-1 lunar regolith simulant and developing spherical particle model with rolling resistance in DEM to model its mechanical behavior and understand the regolith-drill interaction. Laboratory experiments have been performed to characterize the physical and mechanical properties of LHS-1, including particle size, specific gravity, and shear strength. Results from characterization tests are used in the development and calibration of the DEM model. A direct shear test model is developed in DEM. In this study, the LHS-1 particles are modeled as rigid, spherical particles in DEM. To compensate for the shape of the particles, rolling resistance model in DEM is used. The direct shear test model is calibrated and validated using the direct shear test results. Furthermore, the DEM model of a regolith-drill system is developed to demonstrate its application for more advanced applications. The second study focuses on developing an FT4 DEM model for pine residues using multi-sphere particles and understanding the rheological behavior of anatomical fractions of pine residues. The irregular-shaped particles of the pine residue are modeled using the multi-sphere method, which is also known as the clumped-sphere method in DEM. Particle size, shape, stiffness, and density are modeled explicitly for the anatomical fractions of pine residues, i.e., needle, stem, and bark. Efforts have been made to understand the rheological properties of the anatomical fractions. The results show that stems recorded the highest force and torque on the impeller of the FT4 rheometer, which corresponds to the highest flow energy, whereas needles have the lowest flow energy. The results also show that in a whole sample (with a mixture of stem, needle and bark), increasing the stem percentage increases the flow energy while increasing the needle percentage decreases the flow energy. These observations are consistent with experimental data. The analyses of tip speed and particle sizes show that these two factors do not have a clear impact on the recorded force and torque responses. The third study focuses on developing a bonded-sphere DEM-based regression model for bulk density prediction of switchgrass given key biomass characteristics, i.e., moisture content, particle size and size distribution parameter. The predicted bulk density values of switchgrass from the regression model can be used to estimate the mass flow rate of switchgrass. A bulk density model is developed in DEM to generate data points for the regression model. Because there is no breakage in bulk density test, bonded-sphere model without breakage is used to model the deformable particles of switchgrass in this study. The regression model is validated using the experimental results on switchgrass. The cultivars of switchgrass that are used for this work are Cave-in-Rock and Alamo. Therefore, the regression model may provide better results for these two cultivars. Comparison between the experimental-based regression and DEM-based regression show that the DEM-based regression model performs better. The fourth study focuses on developing a bonded-sphere DEM model for predicting switchgrass size reduction in an industrial-scale hammer mill. As part of the mechanical preprocessing for bioenergy conversion, size reduction is an important pretreatment of biomass feedstock. At Idaho National Lab (INL), a proposed feeding system is demonstrated through the process demonstration unit (PDU) and the second stage grinder G2 is a hammer mill in the system. A DEM model is developed for the hammer mill in this study. The particles of the switchgrass are modeled using the bonded-sphere model in DEM with particle breakage enabled. Efforts have been made to better handle the computational expenses while maintaining simulation accuracy. Therefore, the coarse-graining method is applied to the hammer mill model to decrease the number of particles. The results show that coarse-grained DEM model with practical considerations can model the original model (without coarse-graining) effectively. This model can be used to understand the grinding behavior of switchgrass feedstock.