The study of radiative transfer within participating gaseous and particulate media has become increasingly important in the prediction of the combustion process of hydrocarbons for various scientific and industrial applications. The radiative transfer equation (RTE) is an integro-differential equation in five independent variables describing the physical process of radiative transfer. The angular dependency of the RTE makes it exceedingly difficult to solve by deterministic methods. Several approximate deterministic methods for the RTE have been developed over time. Two most promising candidates, the discrete ordinates method (DOM) and the spherical harmonics PN method, are often used to solve the RTE even though both of them have their limitations. The DOM discretizes the entire solid angle by a finite number of ordinate directions and integrals over direction are replaced by numerical quadrature. DOM is relatively simple to implement but suffers from ray effects and false scattering and requires an iterative solution for scattering media or reflecting surfaces. On the other hand, the spherical harmonics PN method is a spectral method that solves the RTE by approximating the angular distribution of the intensity by a truncated series of spherical harmonics. Despite the popularity of the lowest order of the PN method, i.e., the P1 method, the potential of high-order PN methods has never been fully explored. This is partly due to cumbersome mathematics, and to lack of research in this area compared with the effort and progress made in its most popular counterpart, the DOM. Increasing of the order of PN is expected to overcome the difficulty of optically thin and optically intermediate conditions or domains with optically thin and optically intermediate regions, which is the motivation for this research. The Photon Monte Carlo (PMC) method is so far the most accurate method; unlike the DOM/FVM and PN methods, the stochastic PMC method gives an exact solution to the RTE. However, the PMC method can be computationally expensive since a large number of rays must be traced, which prevents it from wider applications in evaluating radiative transfer within combustion simulations.This study focuses on a recently-developed general PN formulation consisting of N(N+1)/2 second-order elliptic PDEs and their Marshak's boundary conditions for arbitrary 3-D geometries. The number of equations and unknowns can be further reduced to (N+1)^2/4 for two-dimensional geometries by taking advantage of the geometric characteristics of spherical harmonics. Special boundary conditions, including symmetry/specular reflection boundaries, walls with specified radiative flux, cyclic boundaries and mixed diffuse-specular surfaces have also been developed for high-order PN methods. The high-order PN methods (up to the order of 7) have been implemented within the finite volume-based OpenFOAM open-source libraries. The performance of high-order PN methods is demonstrated by solving a number of examples covering a wide range of different geometries and varying radiative properties including coupled simulations of a turbulent jet flame and a frozen snapshot study of a high-temperature oxy-natural gas burner. The goal of these examples is to test the performances of the high-order PN methods with respect to all kinds of factors, e.g., order of PN, overall optical thickness, geometry, homogeneity of radiative properties, etc., as well as to verify the finite volume implementations of the high-order PN method on OpenFOAM.