In this paper the author reviews methodology of a version of the global Galerkin that was developed and applied in a series of his 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 types that are combined in 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 deliberated.