Adding squared penalty terms to get rid of constraints is a simple approach giving an accuracy of order 1/penalty factor only. Hence it is not recommended for high accuracy unless you let the penalty go to infinity during the computation. But a high penalty factor makes the Hessian very ill-conditioned,
which limits the total accuracy achievable without taking into account the constraints explicitly.
Note that bound constraints are much easier to handle than general constraints, whence they are virtually never converted to penalties.
The solver L-BFGS-B (used with an about 5-dimensional history) usually solves bound constrained problems very reliably and fast in both low high dimensions.
Exceptions are misconvergence on problems that can become very flat far off the solutions, where it is easy to get stuck with a descent method.
We made lots of experiments on very diverse functions in many different dimensions, with many different solvers available, as we needed a very robust bound-constrained solver as part of our global optimization software.
L-BFGS-B clearly stands out as general purpose method, though of course on sme problems other solvers perform significantly better. Thus I'd recommend L-BFGS-B as a first choice, and would try alternative techniques just in case L-BFGS-B handles your particular class of problems poorly.