Numerical biaxial tests of cohesionless granular materials were conducted using a flexible boundary model to investigate the evolution of the quantity and directional probability of force chains with varied axial strain, under different confining pressures. It has found that the quantity of high-stress particles which don't form force chains decrease as the deviatoric stress increases, indicating that more loads are born by force chains. Statistically, the maximum length of force chains is 9 particles. As axial strain increases, the probability of force chains oriented in the direction [60 degrees, 120 degrees] evolves alike with deviatoric stress. While the probability of force chains oriented in the two directions [40 degrees, 60 degrees] and [120 degrees, 140 degrees] evolves contrary to deviatoric stress. Two types of shear failure were found in biaxial tests: single shear band and two conjugated shear bands. According to the comparison of the probability of force chains oriented in [40 degrees, 60 degrees] and [120 degrees, 140 degrees], a single shear band occurs within the angle range with obviously less probability, while two conjugated shear bands occur when the probabilities of force chains in these two directions have slight difference.