<p>Large eddy simulation (LES) achieves a balance between computational efficiency and high-precision simulation of complex flows through accurately analyzing large-scale eddies while simulating small-scale effects, making it particularly valuable for applications in aerospace engineering, marine hydrodynamics and others. This paper studied a hierarchical radial basis functions method for simulating the classical Smagorinsky LES model. The temporal derivatives are first discretized in fully implicit scheme, then nonlinear terms can be linearized using an iterative method on the basis of direct linearization. Then the spatial derivatives are approximated adopting hierarchical radial basis functions (H-RBFs). Numerical experiments with analytical solutions demonstrate that the H-RBFs method achieves high accuracy and remains unaffected by geometric complexity. Moreover, the iterative method achieves progressively higher accuracy than direct linearization with increasing collocation point density. Finally, benchmark tests including classical flow around the cylinder and backward-facing step flow prove that H-RBFs method effectively captures vortex formation.</p>