Numerical computation of the critical energy constant for two-dimensional Boussinesq equations