In this study, considering the elastoplastic behavior of contacting asperities on micro-scale and the contact pressure distribution of bolted joint interfaces on macro-scale, an improved contact stiffness model is presented. Firstly, the bolted T-stub joint is chosen as the object of study, and the distribution of normal contact load over the interface under preload conditions is studied. Then, combining the continuity, smoothness and monotonicity of the constitutive relations of normal contact load-deformation in three different regimes of plastic- elastoplastic-elastic of contacting asperities, the contact stiffness of individual asperity is analytically derived. Meanwhile, the relation between the distribution of contact stiffness and the distribution of normal contact load for bolted joint interfaces is determined based on fractal theory. The presented model overcomes the shortcomings of the contact stiffness without elastoplastic behavior obtained from the MB model and reveals the inhomogeneous distribution features of interfacial contact stiffness in bolted joints. Finally, the finite element (FE) analysis for flange-bolted joints incorporating the proposed model and the relevant experimental verification are performed. The results show that the average error of the first 4-order natural frequencies between the FE results and the experimental ones is only 0.72%, and the predicted results in frequency-domain are consistent with the measured ones, which illustrates that the proposed model can be used in interface modeling of jointed structures.