In concrete structures, shear failure is one of the failure mechanisms that should be estimated carefully due to its brittle nature. As the technology advances, the capacity of shear in concrete structures can be calculated with the analytical formulation in codes and standards and simulated with NLFEA. The shear capacity is affected by the size which is a phenomenon observed from experiments where the change in the structure size does not have a linear relation to the structural strength. This can be referred to as the experimental size effect. When estimating the shear capacity with NLFEA, this effect must be included in order to avoid over-estimation of capacity. However, another size effect can be found from the result of NLFEA. The change in structural size has an influence on the global response of concrete structures with different numerical configurations in a simulation with NLFEA. A numerical configuration consists of several numerical parameters. This effect is referred to as the numerical size effect. Therefore, it is necessary to incorporate the experimental size effect and reduce the influence of the numerical size effect in a numerical model in order to increase the accuracy and precision of the simulation with NLFEA. The aim of this study is to investigate how the numerical size effect influences the NLFEA results, which will help understand how to improve the accuracy of the simulation. The investigation is done by studying numerical models with a variation on several numerical parameters, such as load increment, error tolerance, the maximum number of iterations and mesh size. The study is focused on identifying the influence of these numerical parameters for cases with different structural sizes and comparable shear slenderness. This is done by observing the global behaviour of the structure and comparing NLFEA results to their respective experimental results. Any difference in the simulated global response is analysed by observing the difference between the experimental and the simulated crack pattern and the convergence state of the respective step. The study is done in three stages. The study is initiated by investigating how the flexural shear failure mode can be obtained with NLFEA. Two specimens with comparable shear slenderness (a/d), 3.70 and 3.91 for beam A123A1 and H123A respectively, are used as a reference for the simulation model. It was found that bad convergence was achieved at the later stage of the analysis and by not accepting non-converged steps, the results became dependent on the convergence state of the analysis where different peak loads were achieved in different numerical configurations for the same case study. The study continues by identifying the correlation of numerical parameters with the results of NLFEA by creating simulation models with different iterative incremental solutions, load increment, error tolerance, and the maximum number of iterations. All load steps are accepted in this stage, regardless of their convergence state, and the rotating crack model is used in these numerical models. The NLFEA results are post-processed by setting up an upper limit in the error tolerance in order to minimize the variation of the peak load and attain a higher consistency of the results. At the last stage, another study is performed to identify the correlation between mesh sizes with the results of NLFEA by simulating several models with different mesh sizes. All load steps are also accepted in this stage, regardless of their convergence state, and the rotating crack model is used in these numerical models The NLFEA results are also post-processed with the same criteria attained in the previous stage. It is concluded that the results of NLFEA are dependent on the following numerical parameters, the load increment, error tolerance, and the maximum number of iterations, due to different initiations and propagation of dowel crack in the non-converged steps, which influences the convergence condition in the later stage of the NLFEA. The formation and rapid propagation of the dowel crack results in the excessive change in the principal strain direction and magnitude, which induces high relative energy variation and out-of-balance force. These numerical parameters indirectly contribute to the error found in each step in the analysis. An additional study is done by replacing the crack model with a fixed crack model with a damage-based shear retention model, but it was found that this did not solve the problem due to the excessive change in the shear retention factor on the elements that lead to premature failure of the beam. The observation related to the numerical size effect reveals that the mesh size controls the propagation rate of the dowel crack. The use of a smaller mesh size will reduce the propagation rate of the dowel crack after its initiation due to a smaller increment of the crack width between each load step. As a result, the increment of the crack width becomes smaller as the mesh size decrease. This becomes more pronounced due to the effect of excessive change of the principal stress-strain direction and magnitude. Full Newton-Raphson method has been proven to give the smallest variation of the peak load ratio (Pnorm), ranging from 0.887-1.140 for specimen A123A1 and 1.055-1.246 for specimen H123A, compared to the other iterative-incremental method with the use of load increment ratio (load increment/experimental peak load) within 0.03-0.006 in combination with evaluation of acceptable analysis step with error tolerance, where the predefined error tolerance must be set to Gerr = 0.0001 and Ferr = 0.01, and non-converged steps can still be accepted if the relative energy variation and out-of-balance force of the corresponding step are below 0.03 and 0.58 respectively. An additional criterion related to the shear displacement on the major flexural shear crack (Δ) should be added to the evaluation of the acceptable analysis step in order to prevent an excessive crack opening on the major flexural shear crack at the peak load. The analysis step can still be considered valid if the respective step satisfies the criterion of the error tolerance and the vertical crack opening on the major flexural crack on the respective step is below the ultimate tensile crack width (CWt,u), which can be calculated from the Guidelines for Nonlinear Finite Element Analysis of Concrete Structures (Hendriks & Roosen, 2020). A mesh size ratio which larger than h/15 is required to have the correct simulated failure.