Since direct numerical simulations of natural convection flows cannot be performed at high Ra-numbers, a dynamically less complex mathematical formulation is sought. In the quest for such a formulation, we consider regularizations (smooth approximations) of the nonlinearity. The regularization method basically alters the convective terms to reduce the production of small scales of motion by means of vortex stretching. In doing so, we propose to preserve the symmetry and conservation properties of the convective terms exactly. This requirement yields a novel class of regularizations that restrain the convective production of smaller and smaller scales of motion by means of vortex stretching in an unconditional stable manner, meaning that the velocity cannot blow up in the energy-norm (in 2D also: enstrophy-norm). The numerical algorithm used to solve the governing equations preserves the symmetry and conservation properties too. The regularization model is successfully tested for a 3D natural convection flow in air-filled (Pr = 0.71) differentially heated cavity of height aspect ratio 4 at Ra = 10 10 and 10 11 . Moreover, a method to dynamically determine the regularization parameter (local filter length) is also proposed and tested.