In this paper the author reviews a version of the global Galerkin that was developed and applied in a series of earlier publications. The method is based on divergence-free basis functions satisfying all the linear and homogeneous boundary conditions. The functions are defined as linear superpositions of the Chebyshev polynomials of the first and second klinds that are combined into divergence free vectors. The description and explanations of treatment of boundary conditions inhomogeneities and singularities are given. Possible implementation for steady state solvers, path-continuation, stability solvers and straight-forward integration in time are discussed. The most important results obtained using the approach are briefly reviewed and possible future applications are discussed.