AbstractA numerical framework for simulations of wake interactions associated with a wind turbine column is presented. A Reynolds‐averaged Navier‐Stokes (RANS) solver is developed for axisymmetric wake flows using parabolic and boundary‐layer approximations to reduce computational cost while capturing the essential wake physics. Turbulence effects on downstream evolution of the time‐averaged wake velocity field are taken into account through Boussinesq hypothesis and a mixing length model, which is only a function of the streamwise location. The calibration of the turbulence closure model is performed through wake turbulence statistics obtained from large‐eddy simulations of wind turbine wakes. This strategy ensures capturing the proper wake mixing level for a given incoming turbulence and turbine operating condition and, thus, accurately estimating the wake velocity field. The power capture from turbines is mimicked as a forcing in the RANS equations through the actuator disk model with rotation. The RANS simulations of the wake velocity field associated with an isolated 5‐MW NREL wind turbine operating with different tip speed ratios and turbulence intensity of the incoming wind agree well with the analogous velocity data obtained through high‐fidelity large‐eddy simulations. Furthermore, different cases of columns of wind turbines operating with different tip speed ratios and downstream spacing are also simulated with great accuracy. Therefore, the proposed RANS solver is a powerful tool for simulations of wind turbine wakes tailored for optimization problems, where a good trade‐off between accuracy and low‐computational cost is desirable.