The efficient solution of large sparse linear systems is essential in computational fluid dynamics (CFD) applications, where it serves as a foundation for computing numerical solutions of partial differential equations (PDEs). These numerical solutions are required for simulating fluid flows, heat transfer and combustion in real-world scenarios. In CFD applications, the accurate model of physical phenomena relies on the discretization of continuous PDEs into algebraic systems, typically in the form of large sparse linear systems. The size of these systems is derived from the spatial discretization in which the number of finite volumes is the same as the dimension of the coefficient matrix. As the computational domain becomes more detailed and has more elements, the resulting linear systems of equations show sparsness patterns that can significantly impact the efficiency of iterative linear solvers [1]. The algebraic multigrid (AMG) method [2] has proven to be an efficient scalable parallel solver [3], but recent high-performance computers present new challenges for its parallel usage. This study presents an algorithm which uses polynomial smoothing as an easy and effi- cient way to improve the efficiency and scalability of the AMG solver. As the literature suggests, there is no single ”best” polynomial; the proper choice depends on the eigen- value distribution, which is rarely known a priori [4]. Specifically, we demonstrate how the Chebyshev polynomial is attractive for usage on parallel architectures, due to the lack of memory-intensive inner product calculations [5]. The newly implemented algorithm will be tested for single-phase, incompressible, laminar and turbulent flows. The proposed algorithm’s expected outcome is a faster solution phase. We expect that our findings will help address the challenges that arise when using recent high-performance computers and that the proposed algorithm can serve as an effective solution procedure for solving large sparse linear systems in CFD applications.