The accurate simulation of steel structures requires a precise model of the joint behaviour. The methods proposed by the steel codes are based on either rotating springs or involved models of springs and rigid bars. In this article, a precise method to model the stiffness of 2D bolted steel connections is presented. First, the joint is accurately modelled using finite elements (FE). Then, the FE model is condensed to a cruciform element of 4 nodes (12 degrees of freedom) by constraining each side cross-section to a node located at its centre of gravity. Subsequently, forces are applied to each node to compute the flexibility matrix, which is then used to construct the stiffness matrix that is finally decomposed through singular value factorization. Following this procedure, a parametric study is conducted to build the training and validation sets of the metamodel. Kriging and Radial Basis Functions are chosen to metamodel and predict the stiffness matrices of the cases not included in the parametric study. Finally, steel structures are analysed with both complete finite elements and surrogate models, and the results are used to confirm the accuracy of the proposed method.