This study concerns one of the most important aspects of Theoretical Fluid Mechanics, namely the regularity of the incompressible Navier-Stokes equation (NSE). The derivation of NSE largely relies on the smoothness of the solution and therefore a careful analysis of its regularity would provide a significant insight into the nature of a fluid motion. In 1989, Foias and Temam showed that the solutions to the 2D NSE on a periodic domain are analytic in time, but in the 3D case this result turned out to be true only for a small interval of time. It means that there is still no guarantee that even if the initial condition v0 is smooth, the solution to the 3D NSE is also smooth. If the amplitude of the periodic v0 is su fficiently small, then unique and smooth solutions are proven to exist for all time. Nevertheless, it is still an open question if smooth initial data with a large magnitude lead to some kind of singularity. There are many different methods of studying this problem, and one of them involves the enstrophy, which is an L2-norm of the vorticity. The usefulness of this quantity relies of the fact that singularity formation strictly depends on the behavior of the enstrophy. Thus, a rigorous analysis of this function may provide an answer to the question concerning the potential blow-up of the solution. Although there is no blow-up in the 2D Navier-Stokes, even on bounded domains, it is still worthwhile to study how much enstrophy can an incompressible flow produce, depending on the initial enstrophy E0. The primary goal of this project is to develop and validate computational tools enabling us to study such problems on bounded domains. The proposed numerical algorithm is based on the vorticity transport equation, while the core mathematical tools allowing us to find the maximum enstrophy growth involve variational optimization. Numerical calculations are performed using a Chebyshev spectral collocation method.