A controlled surface wrinkling pattern has been widely used in diverse applications, such as stretchable electronics, smart windows, and haptics. Here, we focus on hexagonal wrinkling patterns because of their great potentials in realizing anisotropic and tunable friction and serving as a dynamical template for making non-flat thin films through self-assembling processes. We employ large-scale finite element simulations of a bilayer neo-Hookean solid (e.g., a film bonded on a substrate) to explore mechanical principles that govern the formation of hexagonal wrinkling patterns and strategies for making nearly perfect hexagonal patterns. In our model, the wrinkling instabilities are driven by the confined film expansion. Our results indicate robust hexagonal patterns exist at a relatively small modulus mismatch (on the order of 10) between the film and substrate. Besides, the film expansion should not exceed the onset of wrinkling value too much to avoid post-buckling patterns. By harnessing the imperfection insensitivity of one-dimension sinusoidal wrinkles, we apply a sequential loading to the bilayer structure to produce the nearly perfect hexagonal patterns. Lastly, we discuss the connection between the simple bilayer model and the gradient structures commonly existed in experiments.