When it comes to predicting urban airflow, steady Reynolds-averaged Navier-Stokes (SRANS) models that rely on Reynolds stress often face a challenge called the closure problem. This problem involves unresolved structural flaws and uncertainties in the closure coefficients used in the models. Previous attempts to recalibrate coefficients for specific urban flows without breaking the linear constitutive relation have resulted in simulation results constrained by the baseline turbulence model. Therefore, this study aims to enhance the performance of SRANS models by addressing these structural flaws. To achieve this, a novel data-driven framework is proposed. It leverages the deterministic symbolic regression algorithm to discover explicit algebraic expressions for a non-linear Reynolds stress correction model. The robustness of the correction model is ensured by maintaining the linear eddy viscosity model for iterative calculations while keeping the non-linear component frozen. The proposed framework is evaluated using three isolated building cases with varying geometric configurations and inflow boundary conditions. Findings demonstrate that computational fluid dynamics (CFD) predictions incorporating the data-driven non-linear correction model consistently align closer to wind tunnel experimental results compared to both standard and non-linear versions of the k-ε turbulence model. This improvement is reflected in reduced reattachment lengths and more accurate mean velocity distributions in the wake of buildings. However, it should be noted that there is a possibility of overpredicting wind velocity in the windward area. This study introduces valuable insights and additional strategies to enhance the prediction accuracy of SRANS models in urban airflow simulations.