Received: 30 Aug 2019
Revised: 03 Oct 2019
Accepted: 11 Nov 2019
Published online: 18 Oct 2019
Geometric Optimization of Radiative Enclosures using Sequential Quadratic Programming Algorithm
Shuangcheng Sun^{ }^{1,2*}
^{1 }School of Energy and Power Engineering, Chongqing University, Chongqing, 400044, China
^{2 }Key Laboratory of Lowgrade Energy Utilization Technologies and Systems, Chongqing University, Ministry of Education, Chongqing 400044, China
*Corresponding author: Shuangcheng Sun
School of Energy and Power Engineering, Chongqing University
174, Shazheng Street, Shapingba District, Chongqing, China, 400044
Email: scsun@cqu.edu.cn (S. Sun)
Abstract: The sequential quadratic programming (SQP) algorithm is introduced to optimize the geometry of radiative enclosures filled with participating medium. The design goal is to produce a uniform distribution of radiative heat flux on the prespecified surface of a 2D radiative enclosure. The direct problem involves radiative heat transfer in the participating medium is solved by the discrete ordinate method in a bodyfitted coordinate system. The inverse design task is formulated as an optimization problem. The SQP algorithm is applied to optimize the spatial positions of control points, and the Akima interpolation technique is employed to approximate the geometry of the design surface. The effects of radiative properties of participating medium, the number of control points, and temperature of heating surface on the optimization results are investigated. Retrieval results show that the design goal can be successfully satisfied using the SQP algorithm, and the present methodology is proved to be more efficient than other optimization techniques.
The SQP algorithm was applied to optimize the geometry of radiative enclosures, and hence uniform distributions of radiative heat flux on the design surface were obtained.
Keyword: Inverse radiation problem; Geometric optimization; Radiative heat transfer; Sequential quadratic programming; Inverse geometry design
Radiative enclosures exist in numerous engineering fields, such as solar energy harvesting and conversion, industrial boilers, radiative cooling, and thermal processing.^{15} Geometric optimization of radiative enclosures is of great importance to energy application and environment protection. In recent years, inverse design technique has been successfully applied to optimize the enclosure geometry, in which the design task is formulated as an optimization problem. Compared with conventional forward “trialanderror” method, inverse design approach can achieve an optimal solution in far less computational time.^{68} Quite a few optimization techniques have been developed to solve inverse design problems, through which uniform or desired distributions of temperature and radiative heat flux on the prespecified surface are satisfied.^{939}
In general, the theoretical methodologies for solving geometric optimization problems can be classified into two categories: stochastic algorithms and gradientbased methods. Stochastic algorithms are easy to implement and independent on initial solutions. Meanwhile, many candidates can be generated at each iteration. Lots of inverse design problems have been successfully solved using stochastic algorithms. For instance, Sarvari^{9} applied genetic algorithm (GA) to determine the optimal shape of a 2D radiative enclosure to produce a desired heat flux distribution on the specified design surface. Farahmand et al.^{11} used particle swarm optimization (PSO) algorithm to design the geometry of a 2D radiative enclosure to create the desired temperature and heat flux distributions, and PSO algorithm achieved higher accuracy than GA. Qi et al.^{15, 16} employed PSO and krill herd (KH) algorithms to optimize the geometry of a 2D complex radiative enclosure, through which uniform distributions of radiative heat flux on the design surface were obtained. Lemos et al.^{20} applied generalized extreme optimization (GEO) algorithm to determine the geometrical configuration of filament heaters to produce uniform distributions of temperature and radiative heat flux on the design surface. The stochastic algorithm is proved to be effective in solving geometric optimization problems. However, a large iteration number is required for stochastic algorithms, and thereby the optimization process is timeconsuming.
Gradientbased methods are robust to solve optimization problems and possess the advantages of high efficiency and accuracy. Howell et al.^{1, 7, 12} investigated the inverse design of 2D and 3D radiative enclosures by using KieferWolfowitz, steepest descent, Newton, quasiNewton, and conjugate gradient methods. Retrieval results showed that the inverse design approach was usually more accurate than the forward ‘‘trialanderror’’ design methodology and required far less design time. Sarvari et al.^{26, 27} applied conjugate gradient method (CGM) to optimize the heat source distribution in 2D radiative enclosures filled with gray and nongray participating media, through which the desired temperature and heat flux distributions on the design surface were obtained. Tan et al.^{30, 31} used CGM to design the geometry of radiative enclosures to produce uniform distributions of radiative heat flux on the design surface. Rukolaine^{32} studied the geometric optimization of a 2D radiative enclosure using a combination of random search method and CGM. Daun et al.^{39} conducted a comparison of Tikhonov method, truncated singular value decomposition (TSVD), quasiNewton method, conjugate gradient regularization, and simulated annealing in solving inverse design problems. Retrieval results indicated that conjugate gradient regularization usually provides better optimization accuracy in less computational time and with less storage requirement than other inverse techniques.
Sequential quadratic programming (SQP) algorithm, as a typical gradientbased method and with superlinear convergence ability, is widely used for solving optimization tasks. Qi et al.^{4043} applied SQP algorithm to solve inverse radiation problems, through which the optical property distributions of participating medium were accurately recovered. Kim and Hielscher^{44} proposed a reduced Hessian SQP algorithm to reconstruct the distributions of absorption and scattering coefficients, which can be used for medical imaging. Sun et al.^{45 }employed SQP algorithm to estimate the time and spacedependent heat flux on the boundary of participating medium, and SQP algorithm was proved to be more efficient than CGM, LevenbergMarquardt (LM) method and PSO algorithm. However, to the best knowledge of the author, the application of SQP algorithm for solving geometric optimization problems has not yet been investigated. The main goal of this research is to introduce SQP algorithm to optimize the geometry of radiative enclosures to produce a uniform distribution of radiative heat flux on the design surface. The remainder of this work is organized as follows. The radiative heat transfer in irregular radiative enclosures filled with participating medium is solved in Section 2. The physical model of geometric optimization and the fundamentals of SQP algorithm are described in Section 3. The geometric optimization tasks are solved using the SQP algorithm, and the retrieval results are illustrated in Section 4. Meanwhile, the effects of radiative properties of participating medium, the number of control points and temperature of heating surface on the optimization results are investigated in this section. The main conclusions are provided in Section 5.
For the geometric optimization of radiative enclosures, the direct problem involves radiative heat transfer in enclosures will be resolved at each iteration. Hence, an accurate and efficient direct model is of great significance to the optimization process. The radiative transfer equation (RTE) is applied to describe the heat transfer in participating medium, which can be expressed as:^{46}
(1)
where I(s,W) indicates the radiative intensity at position s and direction W. b_{e}, k_{a} and k_{s} denote the extinction, absorption and scattering coefficients of the participating medium, respectively, and b_{e} = k_{a} + k_{s}. The scattering albedo of the medium can be denoted by $\mathit{\omega}\mathbf{=}{\mathit{\kappa}}_{\mathbf{s}}\mathbf{/}{\mathit{\beta}}_{\mathbf{s}}$. n represents the refractive index of the medium. ${I}_{b}={n}^{2}\sigma {T}^{4}/\pi $ is the radiative intensity of blackbody at temperature T, $\sigma $ is the StefanBoltzmann constant. $\Phi (\Omega \text{'},\Omega )$ represents the scattering phase function, $\Omega \text{'}$and $\Omega $ indicate the incident and scattering directions, respectively. Two kinds of scattering phase functions including linear and HenyeyGreenstein (HG) forms are considered and given as follows:
$\Phi (\Omega \text{'},\Omega )$=1+gcos$(\Omega \text{'},\Omega )$ (2)
(3)
where g is the scattering symmetry factor. g = 0, g > 0 and g < 0 represent isotropically scattering, forward scattering and backward scattering, respectively.
The boundaries of radiative enclosures are assumed to be diffuse gray walls, and thereby the boundary condition for the RTE can be given as:
(4)
The discrete ordinate method (DOM) is employed to solve the direct problem, and the discretized RTE can be expressed as:^{ 46, 47}
(5)
where a and b represent directional cosines. w denotes the quadrature weight. N$\Omega $ is the number of discrete direction. The superscripts m and m' represent discrete directions.
Meanwhile, the discretized boundary condition can be expressed as:
(6)
In the optimization process, the geometric shape of the radiative enclosure is irregular and it is different at each iteration. Thus, the RTE in orthogonal Cartesian coordinates should be transformed into an equation in general bodyfitted coordinates, which can provide a more accurate prediction of the radiative transfer within the participating medium and on the medium boundary.
The Jacobian matrix is used for coordinate transformation, which is defined as:
(7)
Hence, the RTE in bodyfitted coordinate system can be expressed as:
(8)
The control volume method is applied to discrete Eq. (8), and the spatially discretized RTE in a bodyfitted coordinate system can be expressed as:
(9)
where the subscripts e, w, n and s represent the eastern, western, northern and southern boundaries of the control volume, respectively. The subscript P represents the central node of the control volume.
The step scheme is used to close the above RTE, and thereby Eq. (9) can be transformed into the following form:
(10)
where the subscripts E, W, N and S represent the central nodes of the eastern, western, northern and southern control volumes around control volume P, respectively. The coefficients can be given as:
(11)
The conjugate gradients stabilized technique is applied to solve the above discretized RTE, and thereby the distribution of radiative intensity in the medium can be obtained. The radiative heat flux on the boundary of the medium can be calculated by:
(12)
To verify the accuracy and feasibility of the present DOM for solving radiative heat transfer problems in participating medium, two simulation cases involve radiative equilibrium and nonradiative equilibrium are considered. The parameter settings for verification cases are set as same as those in Refs.^{48, 49} as listed in Table 1. As shown in Fig. 1, the retrieval results of the present solution are in good agreement with those of Ref.^{48} solved by collapsed dimension method and Ref.^{49} solved by finite volume method, which illustrates the present DOM in a bodyfitted coordinate system is accurate to solve radiative heat transfer problems in participating medium.
Table 1 Parameter settings of verification cases.
Fig. 1 Verification of the present solution for solving (a) radiative equilibrium and (b) nonradiative equilibrium problems in 2D participating medium.
3.1 Description of geometric optimization
This study aims to produce a uniform distribution of radiative heat flux on the specified surface by optimizing the geometry of radiative enclosure. A 2D radiative enclosure filled with participating medium is considered in the present research. As shown in Fig. 2, the bottom boundary AB is the heating surface with temperature T_{S}, and the other surfaces are cold. The left boundary AC and right boundary BG are adiabatic diffuse walls. The geometric shape of DEF is the surface to be optimized, and a uniform distribution of radiative heat flux on this surface is required.
Fig. 2 Physical model for the geometric optimization of radiative enclosures.
The design task is formulated as an inverse problem which can be solved using optimization techniques. According to the design goal, the objective function can be established as:
(13)
where N is the number of discrete node on the design surface DEF. The subscript av indicates an average value, and Q represents the dimensionless radiative heat flux which can be given as:
(14)
In the optimization process, the design surface DEF is discretized into a series of computational nodes. The geometric shape of the design surface can be fitted by using interpolation techniques based on several control points, which can significantly enhance the optimization efficiency of design tasks. Akima cubic interpolation is an efficient tool for obtaining smooth surface, and it is employed as the interpolation method in this study. The fundamentals of Akima cubic interpolation are available in Refs.^{15, 50} and not repeated here.
To evaluate the optimization results, the relative error is defined as:
(15)
and the average and maximum relative errors are denoted by ${\overline{)\epsilon}}_{rel}$ and ${\epsilon}_{rel,max}$ , respectively.
3.2 Sequential quadratic programming algorithm
The design goal of producing a uniform distribution of radiative heat flux on the surface DEF can be satisfied by minimizing the objective function. SQP algorithm is robust to solve optimization problems, in which the optimization task is transformed into a series of quadratic programming (QP) subproblems. The algorithm will superlinearly converge to the optimum by solving these QP subproblems.
The general nonlinear optimization form for the geometric optimization of radiative enclosures can be expressed as:
(16)
where x indicates the parameters to be optimized (i.e. the spatial positions of control points). c_{i} represents the restriction. The numbers of equality restriction and total restriction are denoted by m_{e} and m, respectively.
The original optimization problem is transformed into the following QP subproblems:^{51}
(17)
where the superscript k indicates the kth iteration. d represents the search direction. H is an approximation of the Hessian matrix of the following Lagrangian function:^{51}
(18)
where λ is the Lagrangian multiplier.
The update of H is critical to the optimization process. Superlinear convergence ability of SQP algorithm can be maintained only when H is a good approximation of Hessian matrix.^{40, 42, 52} To ensure the optimization efficiency of SQP algorithm, the approximation H is generated by:^{52}
(19)
where are determined by:
where δ represents a positive constant and is taken as 0.1≤δ≤0.2.
The spatial positions of control points can be updated by:
^{ }（22）
However, only local optimum can be achieved by the above SQP algorithm, which leads to the optimization results will strongly depend on the initial guessed solution. In order to improve the global optimization ability of SQP algorithm, a penalty function is introduced into the present technique:^{53}
(23)
where r represents the penalty factor and is updated by:
(24)
where ρ indicates a positive constant value.
Meanwhile, the optimization parameter is updated by the following equation instead of Eq. (22):
x^{k+1}=x^{k}+φ^{k}d^{k} (25)
where φ represents the searching step size which satisfies:
F_{p}（x^{k}+φ^{k}d^{k} ）p (x^{k})+θφ^{k}[$\overline{\mathrm{F}}$_{p}（x^{k},d^{k} ）F_{p}(x^{k})] (26)
where 0<φ<1, and $\overline{\mathrm{F}}$_{p} is given as:
(27)
Since the subproblem (17) is an approximation of the original optimization task (16), it is possible that no feasible solution can be achieved by solving the subproblem (17) even though there is feasible solution in the original problem (16). If the approximation is failed, the following QP subproblem is considered:
(28)
where ξ_{i}, η_{i }and η_{i } are additional variables.
In the practical optimization process, the searching step size φ^{k} may be less than 1, which means the SQP algorithm will lose the superlinear convergence ability. To maintain the global optimization and avoid the above Maratos effect, the following modified QP subproblem is introduced [54]:
(29)
where p and a are determined by:^{54}
(30)
(31)
A conic is defined to provide search paths, and the arc is given as:
x^{k+1}=x^{k}+φ^{k}d^{k}+(φ^{k})^{2}+${\hat{d}}^{k}{d}^{k}$
F_{p}（x^{k+1})p (x^{k})+θφ^{k}[$\overline{\mathrm{F}}$_{p}（x^{k},d^{k} ）F_{p}(x^{k})] (33)
The iterative optimization stops until one of the following criteria is satisfied:
(34)
(35)
where ε indicates the specified convergence accuracy. u is a positive constant and it is set as v = 1 in this study.
3.3 Computational procedure of geometric optimization
The main computational procedure of the SQP algorithm for optimizing the geometric shape of radiative enclosures can be summarized as follows:
Step 1: Set the iteration number as k = 0, and define the initial parameters such as H^{0} and r^{0}. Give the initial solution x^{0}.
Step 2: Solve the QP subproblem (17). If there is no feasible solution is obtained, solve the QP subproblem (28) and then go to Step 4.
Step 3: Determine whether Maratos effect occurs. If so, solve the QP subproblem (29).
Step 4: Update the penalty factor r^{k} based on Eq. (24). If the solution of the QP subproblem (29) has been obtained, go to Step 6.
Step 5: Conduct the onedimensional search and calculate the step size j ^{k} which satisfies Eq. (26). Update the spatial positions of control points based on Eq. (25) and then go to Step 7.
Step 6: Conduct the search work on the conic and calculate the step size j ^{k} which satisfies Eq. (33). Then, update the spatial positions of control points based on Eq. (32).
Step 7: Solve the direct problem to get the radiative heat flux distribution on the design surface, and then calculate the penalty function.
Step 8: Determine whether the convergence criterion is satisfied. If so, stop the iterative computation and output the optimization results. Else, go to the next step.
Step 9: Update the approximation matrix H^{k} and set the iteration number k = k + 1, and then go to Step 2.The computation flowchart of SQP algorithm is illustrated in Fig. 3.
Fig. 3 Computational flowchart of SQP algorithm for solving geometric optimization problems.
4.1 Geometric optimization of radiative enclosures
The geometric optimization of 2D radiative enclosures filled with gray participating medium is investigated in this section. The extinction coefficient and scattering albedo of the medium are set as β_{e} = 2.0 m^{}^{1} and ω = 0.5, respectively. The medium is assumed to be isotropically scattering medium (i.e.,Φ = 1.0). The temperature and emissivity of the heating surface are taken as T_{S} = 1000 K and ε_{ω} = 1, respectively. The initial geometric shape of the radiative enclosure is considered as a rectangular enclosure with size of L_{x}×Ly=1.0 × 1.0m^{2}. The computational domain is discretized into 42×42 nodes, and S_{6} quadrature scheme is employed in DOM. The number of control points is set as N_{d} = 4. The SQP algorithm is applied to optimize the spatial positions of control points, and the geometry of the design surface DEF is fitted using Akima cubic interpolation. The convergence accuracy is set as ε = 10^{}^{8}. Each case is implemented using the FORTRAN code, and the developed program is executed on an Intel Core(TM) i58400 PC.
The geometric shape of the design surface and the dimensionless radiative heat flux distribution on this surface are illustrated in Figs. 4 and 5, respectively. As shown, a uniform distribution of radiative heat flux on the design surface is achieved. The average and maximum relative errors are 0.0324 % and 0.0781 %, respectively (see Fig. 6). Thus, the design goal can be excellently satisfied by using the SQP algorithm.
Fig. 4 Initial and optimized geometries of the design surface.
Fig. 5 Initial and final radiative heat flux distributions on the design surface.
Fig. 6 Relative error distribution of radiative heat flux on the design surface.
To illustrate the performance of the present inverse technique for solving the geometric optimization of radiative enclosures, a comparison between SQP algorithm and other optimization techniques is conducted. Five typical inverse methods are considered here, namely CGM, LM method, PSO algorithm, KH algorithm, and GA. The convergence accuracy for all methods is taken as ε = 10^{}^{6} and ε = 10^{}^{8}, respectively. The initial solutions for SQP algorithm, CGM and LM method are set as same values, in which a rectangular enclosure is employed as the initial geometry. Since the swarm diversity is of great importance to intelligent algorithms, the initial spatial positions of control points are randomly generated in the search space^{1.0, 1.1} for PSO algorithm, KH algorithm and GA. Meanwhile, the population size is taken as M = 50 in the above intelligent algorithms. The simulation cases of intelligent algorithms are repeated 20 trials to avoid the randomness of optimization results.
Table 2 lists the retrieval results of different optimization techniques. It can be found that the maximum relative error is less than 0.4%, which indicates that uniform distributions of radiative heat flux on the design surface can be successfully obtained using the abovementioned inverse methods. Moreover, the present SQP algorithm requires much less computational time than other methodologies, especially for highaccuracy requirement. Thus, the SQP algorithm is more efficient for solving the geometric optimization of radiative enclosures than other inverse techniques.
Table 2 Optimization results of different optimization techniques.
4.2 Effect of radiative properties of participating medium
Radiative properties of the participating medium produce direct effects on the radiative heat transfer process, and hence affect the radiative heat flux distribution on the surface of radiative enclosures. Thus, the effect of radiative properties on the geometric optimization results is investigated in this section. The scattering albedo of the medium is taken as ω = 0.5. The geometric optimization results of SQP algorithm for different extinction coefficients are shown in Fig. 7. It can be seen that the larger the extinction coefficient is, the smaller the final radiative heat flux on the design surface will be. Meanwhile, the height of final designed surface decreases with the increment in the extinction coefficient. Uniform distributions of radiative heat flux on the design surface can be achieved using the SQP algorithm for different extinction coefficients. The average relative errors are 0.0330%, 0.0355%, 0.0324% and 0.0176% when the extinction coefficients are taken as β_{e} = 1.0 m^{}^{1}, 1.5 m^{}^{1}, 2.0 m^{}^{1} and 3.0 m^{}^{1} respectively, which demonstrates the robust performance of SQP algorithm for solving geometric optimization problems.
Fig. 7 Geometric optimization results for different extinction coefficients. (a) Initial and final geometries of the design surface, and (b) initial and final distributions of radiative heat flux on the design surface
The extinction coefficient of the medium is set as β_{e} = 1.0 m^{}^{1}, and the scattering albedo is taken as different values. As illustrated in Fig. 8, uniform distributions of radiative heat flux on the design surface are obtained by the SQP algorithm. The final geometric shapes of the design surface and the raditive heat flux distributions on this surface are close for different albedos, which indicates that the scattering albedo has no obvious effects on the geometric optimization results.
Fig. 8 Geometric optimization results for different scattering albedos. (a) Initial and final geometries of the design surface, and (b) initial and final distributions of radiative heat flux on the design surface
The extinction coefficient and scattering albedo of the medium are set as β_{e} = 1.0 m^{}^{1} and ω = 0.5, respectively. Table 3 illustrates the geometric optimization results for different scattering phase functions with different symmetry factors. As shown, good optimization results can be achieved when the radiative enclosure is filled with isotropically scattering, forward scattering and backward scattering media. Meanwhile, the maximum relative error is less than 0.1% for different scattering phase functions and symmetry factors. Thus, for different scattering characteristics of participating medium, the prespecified design goal can be successfully satisfied using the SQP algorithm.
Table 3 Optimization results for different scattering phase functions and symmetry factors
4.3 Effect of the number of control points
The geometry of the design surface is fitted by several control points in the optimization process, and the spatial positions of these control points are optimized by the SQP algorithm. Hence, the number of control points will produce significant effects on the optimization process and final results. The number of control points is taken as N_{d} = 1, 2, 3, 4, 5, 6, 7, 8, and 9, respectively. The spatial positions of control points at the x direction in bodyfitted coordinates are uniformly arranged, and the points symmetric about ξ = 0.5 are taken as same values. The retrieval results for different numbers of control points are shown in Table 4. It can be found that the relative error decreases with the increment in the number of control points from 1 to 8. However, the optimization accuracy decreases when the number of control points is increased to 9, which indicates that too many control points will reduce the sensitivity of radiative heat flux with respect to control points. Moreover, the computational time of oddnumbered control points is less than that of evennumbered (2 < 3, 4 < 5, 6 < 7, 8 < 9), and the computational time increases with the number of control points increases (2 < 4 < 6 < 8, 3 < 5 < 7 < 9). Synthetically considering the optimization accuracy and efficiency, the number of control points is suggested to be N_{d} = 4.
Table 4 Optimization results for different numbers of control points
4.4 Effect of heating source
In the practical engineering applications such as solar energy harvesting and thermal conversation, the energy intensity of external heating source is not invariant. Thus, the effect of the temperature of heating surface on the geometric optimization results is investigated. The retrieval results for different heating temperatures are illustrated in Fig. 9. As shown, the temperature of the heating surface produces few effects on the optimization results in terms of the final geometry and dimensionless radiative heat flux distribution on the designed surface. In other words, a uniform distribution of radiative heat flux on the design surface can be obtained even if the energy intensity of the heating source changes, which is of great significance for the geometric design of energy conversion systems.
Fig. 9 Geometric optimization results for different temperatures on the heating surface. (a) Initial and final geometries of the design surface, and (b) initial and final distributions of radiative heat flux on the design surface
The SQP algorithm is employed for the geometric optimization of radiative enclosures filled with gray participating medium. The DOM in a bodyfitted coordinate system is applied to solve the radiative heat transfer in the irregular enclosure, through which the heat flux distribution on the surface is obtained. The main goal of geometric optimization is to produce a uniform distribution of radiative heat flux on the design surface. The design task is formulated as an optimization problem. The SQP algorithm is applied to optimize the spatial positions of control points, and then the geometric shape of the design surface is fitted by using Akima interpolation technique. Retrieval results show that the specified design requirement can be successfully satisfied by the SQP algorithm. Moreover, the effects of radiative properties of participating medium, the number of control points, and temperature of heating surface on the optimization results are investigated. The following conclusions can be drawn:
(1) The SQP algorithm is proved to be more efficient than other optimization techniques in designing the geometry of radiative enclosures filled with participating medium.
(2) The extinction coefficient of the medium provides a significant effect on the final geometry and radiative heat flux on the designed surface, whereas the scattering albedo has no obvious effects on the optimization results.
(3) The temperature of the heating surface produces few effects on the final geometry and dimensionless radiative heat flux distribution on the design surface.
(4) The number of control points for the geometric optimization of radiative enclosures is suggested to be N_{d} = 4.
Nomenclature
c 
restriction matrix 
d 
search direction 
F 
objective function 
F_{p} 
penalty function 
H 
approximation of Hessian 
I 
radiative intensity, W/(m^{2}sr) 
J 
Jacobian matrix 
L 
size of medium, m 
m 
the total number of restriction 
m_{e} 
the number of equality restriction 
n 
normal vector 
n 
refractive index 
N 
the number of discrete nodes on the design surface 
NΩ 
the number of discrete direction 
q 
radiative heat flux, W/m^{2} 
Q 
dimensionless radiative heat flux 
r 
penalty factor 
s 
spatial position 
T 
temperature, K 
w 
quadrature weight 
x 
optimization parameter or coordinate in orthogonal Cartesian coordinates, m 
y 
coordinate in orthogonal Cartesian coordinates, m 
Greeks symbols


Superscripts/Subscripts 

av 
average value 
b 
blackbody 
e 
eastern boundary of control volume 
E 
central node of eastern control volume 
k 
kth iteration 
m 
mth discrete direction 
n 
northern boundary of control volume 
N 
central node of northern control volume 
P 
central node of control volume 
rel 
relative value 
s 
southern boundary of control volume 
S 
southern surface or central node of southern control volume 
w 
medium wall or western boundary of control volume 
W 
central node of western control volume 
The support of this work by the Natural Science Foundation of Chongqing (CSTC, 2019JCYJMSXMX0441), and the Fundamental Research Funds for the Central Universities (No. 2019CDXYDL0007) are gratefully acknowledged.