## Introduction

The higher turbine inlet temperature means higher thermal efficiency and power output. Therefore, the ability to operate at higher temperature has been an important factor in improving the performance of aero-engines. Over the past few decades, the turbine inlet temperature has increased from 1,000 K in the first generation to more than 2,000 K Yin and Rao (2020) in the latest generation aero-engines. The increase of temperature and the fact that total amount of radiation rises as the fourth power of the absolute temperature Howell et al. (2020) determine that the radiation and coupled radiation-conduction heat transfer play more and more important roles in hot section of aero-engines.

The radiation and coupled radiation-conduction problems in aero-engines involve two equations: radiative transfer equation (RTE) and energy equation. Considering the cylindrical symmetry of aero-engines, these two governing equations can be expressed in more simple way using cylindrical coordinate which could bring a lot convenience like the implement of boundary conditions. Different from conduction and convection in which only spatial dimension need to be concerned, because of the directional nature of radiation, the angular distribution of radiation also need to be calculated for radiation heat transfer problem. This angular dependence adds to two more dimensions, i.e. polar angle and azimuthal angle. This characteristic of angle-dependent would bring great challenges to the numerical computation. Besides, RTE, the governing equation of radaiton heat transfer, can be considered as a kind of convection-dominated equation from the point of numerical computation. The characteristic of convection-dominated of RTE means that it is very difficult to obtain a stable and accurate numerical solution. The features of radiation heat transfer introduced above determine that solving RTE and obtaining the radiation distribution are the most important and difficult part for both radiation and coupled radiation-conduction heat transfer problems.

In recent years, the radiation and coupled radiation-conduction heat transfer in cylindrical system have evoked wide interests of many researchers. As early as 1982, Fernandes and Francis (1982) gave the rigorous formulations of combined conduction and radiation in concentric cylinders and solved it by Galerkin finite element method. Pandey (1989) employed undetermined parameters method to solve this coupled problem for gray and nongray gases contained between infinitely long concentric cylinders with black surfaces. Aouled-Dlala et al. (2007) investigated coupled radiation-conduction heat transfer in gray hollow spheres and cylinders. They used finite Chebyshev transform (FCT) to improve the performance of discrete ordinates method, and adopted Chebyshev polynomials to approximate the angular derivative term instead of finite difference scheme. The results show that FCT is more accurate than traditional discrete ordinate method. Mishra and Krishna (2011) adopted modified discrete ordinate method (MDOM) and lattice Boltzmann method to analyze coupled radiation-conduction heat transfer in infinite and finite concentric cylinders with absorbing, emitting, and scattering medium.

In this paper, the discontinuous spectral element method (DSEM) is adopted to solve the radiation and coupled radiation-coduction heat transfer in cylindrical coordinate system. For the numerical solution of RTE, the space-angle scheme is used to deal with the lack of angular resolution problem which exists in the results of RTE. Space-angle scheme means both the spatial and angular computational domains of RTE are discretized and solved by DSEM. The results demonstrate that the space-angle DSEM performs much better than the traditional hybrid methods. For coupled radiation-conduction heat transfer problem, spectral element method (SEM) is adopted to solve the energy equation after the radiation distribution is obtained by DSEM and then substituted into the energy equation as radiative source term. This kind of DSEM-SEM scheme could avoid using two sets of grid which would cause the increase of computational cost and the decrease of accuracy. And the results show that this DSEM-SEM scheme is feasible to solve the coupled radiation-conduction heat transfer problem. Then, the effects of various geometric and thermal physical parameters are comprehensively investigated. Finally, these methods are further extended to 2D cylindrical system.

The paper is organized as follow. In section 2, RTE in concentric cylindrical medium and numerical method are introduced. In section 3, the numerical results of radiation and coupled radiation-conduction heat transfer in both 1D and 2D cylindrical system are presented and analyzed. Finally, conclusions are summarized in section 4.

## Methodology

### RTE in concentric cylindrical medium

As shown in Figure 1, the present paper studies the radiation and coupled radiation-conduction heat transfer in cylindrical coordinate system. The θ is polar angle and φ is azimuthal angle. The absorbing, emitting and anisotropic scattering medium is filled in concentric cylinders where the subscripts in and out refer to the inner and outer wall, respectively.

### Physical model of concentric cylinders.

The RTE in one-dimensional cylindrical gray medium can be written as Modest and Mazumder (2021),

##### (1)
μr[rI(r,Ω)]r1r[ηI(r,Ω)]φ+βI(r,Ω)=S(r,Ω)

with the boundary conditions:

##### (2a)
I(Rin,Ω)=εinIb(Rin)+(1εin)πμ<0I(Rin,Ω)|μ|dμμ>0
##### (2b)
I(Rout,Ω)=εoutIb(Rout)+(1εout)πμ>0I(Rout,Ω)|μ|dμμ<0

where I(r,Ω) is the radiative intensity at spatial position r along angular direction Ω; Ω(θ,φ)=erμ+eψη+ezξ=ersinθsinφ+eψsinθcosφ+ezcosθ is the angular direction described by polar angle θ and azimuthal angle φ; β=κa+κs is extinction coefficient in which κa and κs are absorption coefficient and scattering coefficient, respectively. S(r,Ω) is source term which is defined as,

##### (3)
S(r,Ω)=κaIb+κs4π4πI(r,Ω)Φ(Ω,Ω)dΩ

where Ib is the black body radiative intensity. The scattering phase function Φ(Ω,Ω) describes the probability of photon scattering from incident direction Ω to the direction Ω.

In this paper, the symmetry of one-dimensional cylinder is used. Therefore, the computational domains of 1D case are,

##### (4)
r[Rin,Rout],φ[0,π],θ[0,π2]

### Discontinuous Galerkin method

The discontinuous Galerkin method (DG method) Cockburn (2003) combines features of the finite element and the finite volume framework. Unlike traditional continuous Galerkin method, the DG method works over a trial space of functions that are only piecewise continuous. This means the solution can possibly has a different values on the shared boundary of different elements. An alternative formulation, the so-called weak formulation, is used to get a representation as a piecewise polynomial, in which the polynomial is discontinuous at the element boundaries. The transport equation is taken as an example to introduce the implementation of DG.

##### (5)
(au)=0xRn

First, the solution domain Rn is decomposed into Nel non-overlapping elements Eeln:

##### (6)
EinEjn=ijRn=el=1NelEeln

Then, we need to derive the weak form of the PDE Equation 5. Multiply both sides of Equation 5 by a differentiable test function v, and integrate the resulting equation over the space E to get,

##### (7)
E((au))v=0

using integration by parts,

##### (8)
EuvE(au)v+E(au)nvds=0

where E means the boundary of element E, n denotes the outward normal vector of the boundary. Equation 8 is referred as a weak formulation of Equation 5.

The calculation of DG is carried out on each element, and DG establishes a link between the values of solution in different elements only through flux on the boundary. Therefore, the numerical flux au¯ instead of au need to be introduced to exchange the information between different elements. In this paper, the classical up-winding scheme is considered, in which the numerical flux is modeled as,

##### (9)
au¯=a{u}+|an|[u]

where the operator {} and [] denote the mean value and the jump value of arguments across element boundary, respectively.

##### (10)
{u}=12(u++u)[u]=12(u+u)

Here the superscript operator “+” and “” denote the values at the boundary inside and outside element, respectively.

### Spectral element discretization

In spectral element discretization Patera (1984), the nodal basis functions on each element are constructed by Chebyshev polynomial expansion. For one-dimensional case in standard computing domain [1,1], the nodal basis functions are Lagrange interpolation polynomials through the Chebyshev–Gauss–Lobatto points αCGL,

##### (11)
αjCGL=cos(j1N1π),j=1,2,N

After the construction of basis functions, quantity u to be solved in any position on element can be approximated as,

##### (12)
u=Nei=1uiΓi

where Ne is the number of nodes in element E. the Γi is the the basis function of node i.

### Discontinuous spectral element method discretization of RTE

In order to use DSEM to discretize RTE in cylinderical coordinate system, it is necessary to use the mathematical equivalent form instead of the usual form of the RTE given by Equation 1. However, as shown in Equation 1, there exists a minus sign before the angular differential term /φ in the RTE. This means there are two different mathematical equivalent forms for Equation 1. It is worth noting that it is critical to choose the proper mathematical equivalent form for the successful implement of DSEM. The two different kinds mathematical equivalent forms are:

##### (13a)
Formar(rμI)+φ(ηI)+βrI=rS
##### (13b)
Formbr(rμI)φ(ηI)+βrI=rS

where the dependences of I on r and μ have been omitted to simplify the expression.

#### Physical meaning of the two mathematical equivalent forms

As mentioned before, the computational domain of 1D case is: r[Rin,Rout],φ[0,π],θ[0,π2], which means η=sinθcosφ0 in the whole computational domain. After applying the DSEM to Form a (equation 13a) and Form b (equation. 13b), respectively, one can obtain:

##### (14a)
FormaEI[μrΓr(η)Γφ+βrΓ]drdφdθ+Enr(μr)IΓ¯dl+Enφ(η)IΓ¯dl=ErSΓdrdφdθ
##### (14b)
FormbEI[μrΓr+(η)Γφ+βrΓ]drdφdθ+Enr(μr)IΓ¯dlEnφ(η)IΓ¯dl=ErSΓdrdφdθ

It is obvious that the only difference between the two mathematical equivalent forms lies in the boundary integral term ±Enφ(η)IΓ¯dl. Therefore, the influence of this boundary integral term ±Enφ(η)IΓ¯dl need to be studied.

First, we consider Form a. For the case when nφ=1 (Figure 2a), using up-winding scheme and substituting Equation 9, Equation 10 and nφ=1 into boundary integral term, we can obtain,

### Diagram of φ boundaries. (a) boundary of nφ=−1. (b) boundary of nφ=1.

##### (15)
Enφ(η)IΓ¯dl=E(1)(η){I}Γdl+E|(1)(η)|[I]IΓdl=E(1)(η)I+Γdl=Enφ(η)I+Γdl

The equation above means that I+ inside the element E as shown in Figure 2a is the upstream physical quantity. That is to say that radiative intensity with larger φ is the upstream and radiation with smaller φ is the downstream. The radiation propagates from larger φ to smaller φ.

For the case when nφ=1 (Figure 2b), Substituting Equation 9, Equation 10 and nφ=1 into boundary integral term, we can obtain,

##### (16)
Enφ(η)IΓ¯dl=E(1)(η){I}Γdl+E|(1)(η)|[I]IΓdl=E(1)(η)IΓdl=Enφ(η)IΓdl

The equation above means that I outside the element E as shown in Figure 2b is the upstream physical quantity. In other words, radiation transfers from larger φ to smaller φ. In conclusion, the physical meaning of Form a is that radiation propagates from larger φ to smaller φ in φ direction.

Similarly, the physical meaning of Form b can be analyzed in the same way. For the case when nφ=1, we can get,

##### (17)
Enφ(η)IΓ¯dl=[E(1)(η){I}Γdl+E|(1)(η)|[I]IΓdl]=E(1)(η)IΓdl=Enφ(η)IΓdl

The equation above means that I outside the element is the upstream physical quantity. That is to say that the radiative intensity with smaller φ is upstream. Radiation propagates from smaller φ to larger φ.

When nφ=1, we will get,

##### (18)
Enφ(η)IΓ¯dl=[E(1)(η){I}Γdl+E|(1)(η)|[I]IΓdl]=E(1)(η)I+Γdl=Enφ(η)I+Γdl

The equation above means that I+ inside the element is upstream physical quantity in φ direction and radiation propagates from smaller φ to larger φ.

In conclusion, the physical meaning of Form b is that in φ direction radiation propagates from smaller φ to larger φ. This conclusion means that the physical meanings of Form a and Form b are just the opposite. Therefore, we need to reveal the real physical scene of radiation propagation in cylindrical coordinate system, and then determine which mathematical equivalent form should be chosen in this paper to implement DSEM.

#### Radiation propagation in cylindrical coordinate system

In angular direction φ, there exists an implicit relationship between the “upstream” and “downstream”. The angular computational domain in the cylindrical symmetric RTE is φ[0,π]. Figure 3a and 3b present two typical propagations of radiation to clarify the implicit relationship of the “upstream” and “downstream” in angular φ direction. It is obvious that the position with larger φ is the “upstream” while smaller φ is the “downstream”, radiation propagates from larger φ to smaller φ.

### The two typical propagations of the radiation in r−φ plane. (a) typical propagation path 1 (b) typical propagation path 2.

It can be seen that the physical meaning of Form a is consistent with the physical reality of radiation propagation in the cylindrical coordinate system while the physical meaning of Form b is just the opposite. Therefore, the Form a (equation 14a) is the mathematical equivalent which we should choose in this paper while the Form b is not allowed to be used.

### Spectral element method discretization of energy equation

For copuled radiation-conduction heat transfer problem in cylindrical coordinate system, different from the RTE which involving three-dimensional rφθ domain, the one-dimensional energy equation with radiative source term is,

##### (19)
d2Tdr2+1rdTdr=4πκak(Ib(T(r))14πG(r))

where k is thermal conductivity, T is the temperature value of the latest iteration step.

Since the coupling between the energy equation and the equation of radiation transfer is highly nonlinear due to the relationship between Ib and the fourth power of temperature T, it is necessary to employ iterative method to solve the copuled radiation-conduction problem.

The energy equation Equation 19 is a classical elliptic equation. Therefore, a good result can be obtained by directly using the spectral element method (SEM) without introducing discontinuous Galerkin scheme. The mesh of SEM adopts the same number of elements and nodes as DSEM to avoid using two sets of computational grid which would cause the increase of computational cost and the decrease of accuracy. Applying the SEM to Equation 19, we can obtain:

##### (20)
EdTdrdΓdrdr+EdTdrΓ(nr1)dl+E1rdTdrΓdr=E4πκak(Ib(T(r))14πG(r))Γdr

### The extension to 2D cylindrical system

Figure 4 depicts the 2D cylindrical system, the z direction need to be considered compared with 1D case. Therefore, from the mathematical point of view, the radiation heat transfer problem in 2D cylindrical system considered here is actually 4D case, the polar angle θ, azimuthal angle φ, spatial r and z.

### Physical model of 2D concentric cylinders.

#### Radiation propagation in the 2D cylindrical system

The RTE in 2D cylindrical gray medium is,

##### (21)
μr[rI(r,z,Ω)]r1r[ηI(r,z,Ω)]φ+ξI(r,z,Ω)z+βI(r,z,Ω)=S(r,z,Ω)

The boundary conditions of gray walls we considered in this paper are:

##### (22a)
I(Rin,z,Ω)=εinIb(Rin)+(1εin)πμ<0I(Rin,z,Ω)|μ|dμμ>0
##### (22b)
I(Rout,z,Ω)=εoutIb(Rout)+(1εout)πμ>0I(Rout,z,Ω)|μ|dμμ<0
##### (22c)
I(r,Zdown,Ω)=εdownIb(Zdown)+(1εdown)πξ<0I(r,Zdown,Ω)|ξ|dξξ>0
##### (22d)
I(r,Zup,Ω)=εupIb(Zup)+(1εup)πξ>0I(r,Zup,Ω)|ξ|dξξ<0

There also exists the symmetry for 2D cylinder and the domains we need to considered are,

##### (23)
r[Rin,Rout],z[Zdown,Zup],φ[0,π],θ[0,π]

Similar to the analysis in 1D cylinder case, for the successful implement of DSEM, the mathematical equivalent form of Equation 21 we need to choose is.

##### (24)
r(rμI)+φ(ηI)+z(ξrI)+βrI=rS

And the corresponding formula after applying DSEM to Equation 24 is,

##### (25)
EI[μrΓr(η)ΓφξrΓz+βrΓ]drdφdθdz+Enr(μr)IΓ¯dl+Enφ(η)IΓ¯dl+Enz(ξr)IΓ¯dl=ErSΓdrdφdθdz

#### Coupled radiation-conduction problem in the 2D cylindrical system

The copuled radiation-conduction heat transfer problem in 2D cylindrical system involves two-dimensional rz domain, the 2D energy equation with radiative source term is,

##### (26)
2Tr2+1rTr+2Tz2=4πκak(Ib(T(r,z))14πG(r,z))

After applying the SEM to Equation 26, we would obtain:

##### (27)
EdTdrdΓdrdrdz+EdTdrΓ(nr[1,0])dl+E1rdTdrΓdrdzEdTdzdΓdzdrdz+EdTdzΓ(nz[0,1])dl=E4πκak(Ib(T(r,z))14πG(r,z))Γdrdz

## Result and discussion

### Grid independence tests

The grid independence tests are conducted in this section. Radiative equilibrium situation is considered firstly, the radius ratio of concentric cylinder Rin/Rout=0.1, the outer wall remains cold, i.e. Tout0. The cylinder is filled with isotropically scattering media, because of the isotropic scattering, the scattering albedo does not change the distribution of radiative heat flux. Figure 5a gives the change of dimensionalless radiative heat flux qr=qr/σTin4 at the inner wall with extinction coefficient β when the number of nodes in the element Mr×Mφ×Mθ changes and the number of elements is fixed as Nr×Nφ×Nθ=8×8×8. The results in Figure 5a show that when Mr×Mφ×Mθ is larger than 8×8×8, the dimensionalless radiative heat flux qr at the inner wall remains almost unchanged and in consistent with the result of the modified discrete coordinate method (MDOM) Mishra and Krishna (2011). Figure 5b presents the changing trend of qr at the inner wall with extinction coefficient under the condition that the number of nodes in the element Mr×Mφ×Mθ=8×8×8 remains constant and the number of elements Nr×Nφ×Nθ changes. It can be seen that, the results are not in good agreement with the those in reference when Nr×Nφ×Nθ=2×2×2. With the number of elements increases from 3×3×3 to 6×6×6, the results agree well with reference and remain basically unchanged.

### Radiation heat transfer in 1D cylindrical system

Note that the radiative equilibrium problem is already considered in the grid independence tests above, then another typical radiation heat transfer problem, non-radiative equilibrium case would be numerically analyzed in the following section.

For non-radiative equilibrium problem, the radiative participating medium is the radiation heat source. The concentric cylinder boundaries remain cold, i.e. Tin=Tout=0. The change of dimensionalless radiative heat flux qr at both inner and outer wall with extinction coefficient β are depicted in Figure 6. In Figure 6a, the scattering albedo ω=0. And it can be seen that for a given value of β, the dimensionalless radiative heat flux qr increases with the increase of radius ratio Rin/Rout in both inner and outer wall. In Figure 6b qr decreases all the time when scattering albedo ω becomes larger. This is because the larger ω means more energy the medium would be scattered, and the radiative heat flux would decrease consequently. It is obvious that the numerical results of DSEM are in good agreement with the results in reference Mishra and Krishna (2011) for the non-radiative equilibrium problem.

### Coupled radiation-conduction heat transfer in 1D cylindrical system

Table 1 shows the dimensionless total heat flux in inner wall qt(Rin) and outer wall qt(Rout) under conditions of different conduction-radiation parameter Ncr and different scattering albedo ω for the case Rin/Rout=0.5, εin = εout = 1, Tout/Tin=0.5 and β=1. The mesh used in this case is Nr×Nφ×Nθ=16×8×8 and Mr×Mφ×Mθ=5×5×5. It can be seen that the numerical results of DSEM-SEM used in this paper are basically consistent with the results of FCT and SCM, Thus, it is verified that the DSEM-SEM can correctly solve the coupled radiation-conduction problem in one-dimensional cylindrical system.

### Values of the dimensionless total heat flux in inner wall qt∗(Rin) and outer wall qt∗(Rout) obtained by DSEM-SEM, FCT Aouled-Dlala et al. (2007) and SCM Sun et al. (2021).

Ncrωqt(Rin)qt(Rout)
FCTSCMDSEM-SEMFCTSCMDSEM-SEM
10.91.64361.64211.63900.82180.82100.8197
0.51.64881.64681.64470.82440.82340.8231
0.11.65371.65121.64990.82680.82560.8260
0.10.93.45233.43633.40631.72611.71831.7053
0.53.50453.48403.46271.75221.74221.7380
0.13.55293.52713.51311.77631.76381.7660
0.010.921.540321.380721.079310.770010.692110.5583
0.521.990721.793721.548010.995310.898810.8178
0.122.317222.088921.867511.158611.046510.9918

After verifying the feasibility of DSEM-SEM scheme, we study the effects of conduction-radiation parameter Ncr on dimensionless temperature for coupled radiation-conduction problem. Figure 7 gives the dimensionless temperature distribution for different Ncr and the case of pure heat conduction. It can be seen that, when Ncr is relatively small, the dimensionless temperature has an obvious nonlinear distribution. However, with the further increase of Ncr, the dimensionless temperature distribution tends to be linear and gets closer to the profile of pure heat conduction.

### Effects of conduction-radiation parameter Ncr on dimensionless temperature distribution.

The conduction-radiation parameter Ncr=kβ/4σTref3 is defined as the ratio of conductive heat transfer and radiation heat transfer. It is obvious that, for a large value of Ncr, conduction plays a dominant role in coupled heat transfer problem. On the contrary, a small value of Ncr means radiation becomes much more pronounced. For the conduction-dominated problem, the dimensionless temperature tends to be linear while the dimensionless temperature tends to be nonlinear for the radiation-dominated problem. Thus, with the increasing of conduction-radiation parameter, the distribution of dimensionless temperature changes from nonlinear to linear.

### The radiation and coupled radiation-conduction heat transfer in 2D cylindrical system

#### The radiation heat transfer in 2D cylindrical system

The non-radiative equilibrium is considered in 2D cylindrical system. The inner, outer, up and down walls of the 2D cylinder are all cold, i.e. Twall0. And all of the four walls are black diffusive boundaries. The height of cylinder is H=ZupZdown=2, Rout=1, radius ratio Rin/Rout=0.5. Figure 8 shows the dimensionless radiative heat flux along the axial z direction at the outer wall of a two-dimensional cylinder with various radiative heat transfer parameters.

### The influence of different radiative heat transfer parameters on dimensionless radiative heat flux qr∗. (a) extinction coefficient β. (b) scattering albedo ω.

Figure 8 shows the results of DSEM in the 2D cylinder case, the results are compared with those of modified discrete coordinate method (MDOM) in ref. Mishra et al. (2011). The comparison shows that in various cases, the results of DSEM are all in good agreement with those in reference, which verifies the feasibility of DSEM in two-dimensional cylindrical system. The effects of different radiative heat transfer parameters on dimensionless radiative heat flux qr will be further analyzed.

Figure 8a shows the influence of extinction coefficient β on qr when the scattering albedo ω=0 and radius ratio is set to Rin/Rout=0.5. It can be seen from the figure that with the increase of extinction coefficient β, the qr increases significantly. This is because when the scattering albedo ω is fixed, the larger the extinction coefficient β means there would be more interaction between photon and medium.Therefore, the radianve heat flux qr would be larger with the increase of β.

The influence of scattering albedo ω on dimensionless radianve heat flux qr is given in Figure 8b, in which the extinction coefficient β and radius ratio Rin/Rout remain unchanged. The scattering albedo ω represents the degree of scattering. As scattering albedo approaches 1, the scattering becomes more stronger, and the medium would absorb less radiative energy. Therefore, when other parameters remain unchanged, the increase of scattering albedo ω would lead to a decrease trend of radiative heat flux. And this trend is reflected in Figure 8b clearly.

#### Coupled radiation-conduction heat transfer in 2D cylindrical system

In this chapter DSEM-SEM scheme would be extended from 1D cylinder case to solve the coupled radiation-conduction problem in 2D cylindrical system. The results solved by Finite Difference Method-Finite Volume Method (FDM-FVM) in ref. Mishra et al. (2011) are also presented for comparison, in which the radiative intensity are obtained by FVM while the energy equation is solved by FDM.

In this case, the height of 2D cylinder H=ZupZdown=1, Rout=1, the medium inside cylinder is isotropic scattering. Radiative heat transfer parameters are: extinction coefficient β=1, scattering albedo ω=0. The outer and up wall are cold wall Tout=Tup=0 while Tdown=100=Tref, the inner wall is adiabatic boundary, which means T/r=0.

Figure 9 shows the influence of radius ratios Rin/Rout on dimensionless temperature Θ of 2D concentric cylinder when the conduction-radiation parameter Ncr=0.1 remains unchanged. The temperature distribution at axial central plane z=(Zdown+Zup)/2 is given. It can be seen that as Rin/Rout decreases, the dimensionless temperature at the axial central plane gradually increases. This is because the space inside the cylinder becomes larger when Rin/Rout decreases, and more heat can be transferred from the down wall Tdown with high temperature to the axial central plane, which causes the increase of dimensionless temperature in this plane. And it can be expected that the temperature distribution would reach the highest when the Rin/Rout=0 (i.e. a solid 2D cylinder), and the corresponding Θ at r=0 in axial central plane would be 0.5.

## Conclusion

In this paper, two different kinds mathematical equivalent forms of RTE in cylindrical coordinate system are studied and compared. Through the analysis from the perspective of physical meaning, the correct form of RTE in cylindrical coordinate system used for DSEM is given. Then, both the spatial and angular domain are discretized by DSEM to solve the radiation heat transfer problem in cylinderical system. Besides, considering the numerical characteristics of both RTE and energy equation, DSEM-SEM scheme is proposed to solve the coupled radiation-conduction problem. Then, several cases are given and analyzed to prove the feasibility for solving radiation and coupled radiation-conduction heat transfer in cylindrical system, Finally, these methods are further extended to 2D cylindrical system which is closer to the real case of aero-engines.