In concrete structures, opened cracks contribute significantly to the transfer of shear and normal stresses through the contact forces occurring between fractured surfaces. Such contact forces are due to protruding asperities, engaged by interlocking and friction. This physical phenomena, termed as Aggregate Interlock, contributes significantly to the shear-strength of cracked concrete structures. The complex phenomenon depends on the topography of cracked surfaces, where material protruding from one side may engage with the opposite one. Various analytical models have been proposed to determine and quantify the contact stresses developing due to aggregate interlock. However, these analytical models do not include the full aspect of the problem such as roughness of crack surfaces at different length scales, elasto-plastic deformation of crack surfaces and process of fracture and interlock simultaneously occurring along the crack surfaces. Therefore,the aim of this thesis, is to investigate the role of such parameters in determining the contact stresses by numerically modeling of contact between concrete surfaces using numerical framework such as Boundary Integral method and Finite Element method. This research proposes to statistically analyze the influence these parameters have on the overall shear-strength of concrete and to use the information thus gained for improving upon the existing analytical-models. The first part of this objective is to better understand the consistency of the widely used analytical approaches such as Two-Phase model and particularly of the physical interpretation of the fitting parameters usually employed to match the experimental results. This thesis suggests an extension of the model by introducing the induced surface alterations (elasto-plastic deformation, degradation of matrix material), which allows accurate semi-analytic predictions of shear stress resistance. The proposed model is further strengthened by a number of tailored experiments, comprising concrete with various aggregate sizes, and various loading kinematics. The second part of the objective is to investigate the role small-scale roughness, at micro-scale, of a crack in concrete plays. In this thesis, we propose a methodology to numerically generate representative rough surfaces in concrete while prescribing characteristic signatures of real concrete surfaces. With help of numerical modeling of contact between such representative surfaces, it is shown that inclusion of small-scale roughness eliminates the need of any fitting parameter, making it a reliable and physically-based method for predicting shear transfer phenomena in concrete. An empirical power-law predicting the shear resistance in concrete is a direct outcome, which accounts for micro-scale roughness and aggregate distribution. The complex problem of aggregate interlock is generally simplified as a problem of two completely fractured surfaces coming in contact. However, the contact between surfaces occurs simultaneously along the fracture process. Therefore, the third part of this thesis is to account for the simultaneous process of fracture and contact between fracturing surfaces. A new numerical framework is developed that couples cohesive zone modeling, representing fracture, and a Node-to-segment contact algorithm to handle normal and frictional contact conditions. The proposed model allows to handle contact between the newly cracked surfaces.