Julien Sorel Djeumegni

Myriam Lazard

Hervé Vital Le Dez

Thierry Tagne

Julien Sorel Djeumegni
email: julien.sorel.djeumegni@univ-poitiers.fr

Hervé Thierry

Tagne Kamdem

Radiative heat transfer in a

Keywords: semi-transparent, semi-analytical, cavity, Bickley-Naylor functions

Nomenclature

Gaussian quadrature abscissa 𝜑 𝑚𝑖𝑛 , 𝜑 𝑚𝑎𝑥 minimum and maximum boundary integrals to be discretized 𝜑 𝑙 azimuthal angle at position 𝑙

Introduction

For several decades, prediction of radiative heat transfer at high temperature in semitransparent media has been a major challenge, giving its relevant applications in several engineering fields [START_REF] Siegel | Thermal Radiation Heat Transfer[END_REF][START_REF] Steward | The calculation of radiative heat flux in a cylindrical furnace using the Monte Carlo Method[END_REF][START_REF] Le Dez | Radiative heat transfer in a semi-transparent medium enclosed in a two-dimensional square cavity[END_REF][START_REF] Salah | Numerical solution of the radiative transfer equation in a cylindrical enclosure with finite volume method[END_REF][START_REF] Wang | Meshless method for solving radiative transfer problems in complex two-dimensional and three-dimensional geometries[END_REF]. The cases of complex or irregular geometries remain particularly interesting in this paper, because they are useful for computational fluid dynamics [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF], while attention is drawn on a two-dimensional geometry with square cavity. For this purpose, providing a solution to this engineering problems firstly involves a deep predictive analysis on heat transfer within the semi-transparent medium concerned.

Early, Sànchez and Smith [START_REF] Sànchez | Surface Radiation Exchange for Two-Dimensional Rectangular Enclosures Using the Discrete -Ordinates Method[END_REF] worked on surface radiation exchange within rectangular enclosures using discrete ordinates method (DOM), in order to better predict heat fluxes. The enclosure concerned was designed with many opaque protrusions and obstructions surfaces, wich scatters, and emits radiation. Constant radiative properties, radiant exchange between surfaces separated by a transparent medium has been computed, to come out with an accomodate arrangement of surfaces, for an optimal prediction. At the end, by overcoming shadowing and obstruction effect, heat flux was calculated and compared with those of irradiation analysis [START_REF] Siegel | Thermal Radiation Heat Transfer[END_REF] with good agreement. Adams and Smith [START_REF] Adams | Three-dimensional discrete-ordinates modeling of radiative transfer in a geometrically complex furnace[END_REF], used a similar method to apply in the case of a three dimensional furnace with internal cooling-pipes behaving as obstacles during the energy transfer by radiation process. The aim was to model directional shadowing effect caused by cooling pipes in radiation intensity field. Hence, predicted values for incident wall flux and radiative transfer were computed and compared with experimental data. In the same order, Chai and coworkers [START_REF] Chai | Treatment of irregular geometries using a Cartesiancoordinates-based discrete-ordinates method[END_REF] did a convenient analysis of protrusions, obstructions and curved or inclined surfaces, with a simple procedure dealing with irregular cartesian coordinates-based geometries. The work was applied to many two-dimensional problems, by taking in account both discrete ordinates method [START_REF] Chai | Treatment of irregular geometries using a Cartesiancoordinates-based discrete-ordinates method[END_REF]and finite volume method (FVM) [START_REF] Chai | Treatment of irregular geometries using a Cartesiancoordinates finite-volume radiation heat transfer procedure[END_REF]. Further, Coelho and coworkers [START_REF] Coelho | Modelling of radiative heat transfer in enclosure with obstacles[END_REF] develop a modelling of radiative heat transfer in enclosures carrying very small thickness of the obstacle. Whereas, in some applications like power station boilers, and in most others relevant engineering problems, radiation is not the only physical phenomenon involved. Hence, they realized a suitable model of radiation heat transfer using discrete ordinates and finite volume methods, in order to couple it with fluid flow codes. At the end, heat flux was compared with those one obtained through zonal methods calculations. The summary is that, predicting by DOM remain the most economical method, despite the fact that, it can be applied to complex geometries, it suffers from shortcoming known as ray effect [START_REF] Chai | Ray effect and false scattering in the discrete ordinates method[END_REF], due to angular discretization in radiative intensity. To overcome this difficulty, an iterative scheme based to generate by stochastic issue, the sets of ordinates directions was implemented [START_REF] Pasini | Numerical experiments on the application of the diamond scheme to sets of discrete directions obtained from a random numbers generator[END_REF]. Another technic which consists to separate the radiative intensity into two parts using modified differential approximations and spherical hamonics method was also adapted [START_REF] Olfe | Radiative equilibrium of a gray medium bounded by nonisothermal walls[END_REF][START_REF] Modest | The modified differential approximation for radiative transfer in general three dimensional media[END_REF], and applied to a two-dimensional and three dimensional enclosures with non-uniform boundary conditions. The idea was also extended by Sakami and Charette [START_REF] Sakami | Application of a modified discrete ordinates method to twodimensional enclosures of irregular geometry[END_REF], through the modified discrete ordinates method (MDO) to solve the problem of ray effect into DOM, in an absorbing, emitting and scattering medium enclosed in a cavity of irregular geometry. They applied it for a prediction of radiative heat transfer in complex geometry with obstrusion/obstacles [START_REF] Sakami | Analysis of radiative heat transfer in complex twodimensional enclosures with obstacles using the modified discrete ordinate method[END_REF]; the aim was to counter the ray effect inherent in DOM, when obstructions are emphasized in enclosure and the results proposed showed smooth agreement with liretature. Kim and coworkers [START_REF] Kim | Unstructured finite-volume method for radiative heat transfer in a complex two-dimensional geometry with obstacles[END_REF] studied radiative transfer in twodimensional enclosure with obstacles such as bafles, pins, or pipes, in order to bring also their contribution to practical engineering applications. They extended the work of [START_REF] Baek | Nonorthogonal finite-volume solutions of radiative heat transfer in a three -dimensional enclosure[END_REF], with unstructured grid system and finite-volume method was considered to predict radiative heat transfer by adopting a triangular mesh. The main geometries investigated were, a square enclosure with finned internal cylinder, a furnace with embedded cooling pipes and a semicircular enclosure with cylinder hole. All the results achieved were in good agreement with the one obtained by zone method. Wilson and Tan [START_REF] Wilson | Radiative transfer in a two-dimensional rectangular annulus medium[END_REF] solve numerically the radiative equation in a two-dimensional rectangular absorbing-emitting, and anisotropically scattering medium with rectangular annulus. Taking in consideration diffusely emitting and reflecting boundaries surfaces, a finite difference method combined with a product of discrete-ordinates quadrature scheme were developped. Their method proposed can be easily modified, to deduce results proposed by previous authors in condition of none existing internal annulus or cavity. Unfortunately, there was also a lack of precision on results proposed, althought they are in good agreement with literature.

Recently, Le Dez and Sadat [START_REF] Le Dez | Radiative transfer in a semi-transparent medium enclosed in a cylindrical annulus[END_REF] did an exact analytical description of radiative transfer in a semi-transparent medium enclosed in a cylindrical annulus, bounded by hot diffusely reflecting surfaces. The authors shown that radiative field can be described in this geomety by purely geometric weighting coefficients. Since, the annulus is a non participating medium, playing the role of obstruction in the geometry, temperature field is computed for any combined heat transfer involving radiation. Their results, applied to some examples was in good agreement with literature. The most recents works dealing with two-dimensional enclosures carrying inner square body, are the one studying coupling between radiation, convection and/or radiation [START_REF] Mezrhab | Computation of combined naturalconvection and radiation heat transfer in a cavity having a square body at its center[END_REF][START_REF] Lari | Combined heat transfer of radiation and natural convection in a square cavity containing participating gases[END_REF][START_REF] Saravanan | Coupled thermal radiation and natural convection heat transfer in a cavity with a heated plate inside[END_REF][START_REF] Chen | Numerical analysis on the radiation-convection coupled heat transfer in an open-cell foam filled annulus[END_REF]. Mostly, the participating medium is filled by a fluid; air or a gaz flowing through enclosure, whereas inner square represents an obstruction. all this works involve numerical methods, such as discrete ordinates and finite volume methods used to solve the governing equations, in order to evaluate temperature field and heat flux inside the participating medium.

Therfore, importance of this present paper is to provide an exact semi-analytical method to modelize radiation within complex geometries, in order to solve with accuracy some thermal engineering problems. Otherwise stated, attention is mainly focussed on heat transfer by radiation and intend to extend the work of Djeumegni et al [START_REF] Djeumegni | Modeling of radiative heat transfer in a gray semi-transparent medium with internal fluid cavity limited by black boundary surfaces[END_REF]. The aim there, is to evaluate through a relevant analysis, the radiative quantities within a two dimensional semi-transparent medium, carrying a non-participating inner square body, for a further combining with convective and/or conductive flow. The mathematical issue is the exact method which keeps an hybrid formulation combining angular and space integrals, to avoid mistakes of inherent effects in computations.

The paper is structured as follows: the description of the geometry and governing equations of the problem, followed by establishment of exact expressions of incident radiation and radiative flux vector. Then, numerical results are presented under tabulated data, with various simulations, to access the correctness of the present method.

Geometry and governing equations

Geometry

Analysis of energy transfer by radiation, inside a two-dimensional semi-transparent medium containing a centered inner square cavity, enhances a deep examination on how radiative quantities should be evaluated.

𝐷

Considering a (𝑥, 𝑦) reference coordinate system Fig. 1a-b, the semi-transparent medium is composed of two sets: the participating one, which is a solid bounded by the enclosure, and a non-participating one, which is represented by the square cavity. The participating medium is designed also by a square having length sides of 𝐻 and the inner square cavity of length ℎ. Propagation of radiation is done through (𝑒 ⃗ 𝑥 , 𝑒 ⃗ 𝑦 ) plane, since z-direction remains infinite, while there is no variation of parameters along this axis. Consequently, because the inner square cavity is present, the geometrical analysis of heat transfer by radiation within the participating medium 𝐷 , becomes more complex. This is illustrated by the fact, when the boundary surfaces of the enclosure are at imposed temperature as in the present case, the geometry of the problem shows existence of five (05) different cases to model ray propagation. They are shown on figures Fig. 234, with all of them depending of the ratio

𝐻 ℎ ∈ [1,3[ ∪ [3, 2 + √5[ ∪ [2 + √5, +∞[. Physically,
it means radiation that reaches each numbered sub-surfaces denoted by 𝜁 𝑛 , with 𝑛 ∈ [START_REF] Siegel | Thermal Radiation Heat Transfer[END_REF][START_REF] Wilson | Radiative transfer in a two-dimensional rectangular annulus medium[END_REF] for Fig. 2a, Fig. 2b and Fig. 4a, 𝑛 ∈ [START_REF] Siegel | Thermal Radiation Heat Transfer[END_REF][START_REF] Sakami | Application of a modified discrete ordinates method to twodimensional enclosures of irregular geometry[END_REF] for Fig. 3.a and 𝑛 ∈ [START_REF] Siegel | Thermal Radiation Heat Transfer[END_REF][START_REF] Talukdar | Discrete transfer method with the concept of blocked-off region for irregular geometries[END_REF] for Fig. 3.b, does not follow the same expression, because the inner square cavity behaves like an obstruction for some incoming rays pathways passing through it. However, propagation of radiation within the participating medium cannot obey to the same law of prediction at each point studied. (𝑏)

The participating medium studied is gray, with constant radiative characteristics and simultaneously emits, absorbs but doesn't scatter radiation. The medium is divided in equal isothermal mesh grid surfaces of ∆𝑥 × ∆𝑦, which also depends of surfaces length as

𝐻 𝑁 𝑥 -1 × 𝐻 𝑁 𝑦 -1
, where 𝑁 𝑥 and 𝑁 𝑦 represent respectively the mesh positions following 𝑥 and 𝑦 directions. Therefore, each point 𝑀 𝑖𝑗 is located at the center of each cell, having coordinates 𝑀 𝑖𝑗 (𝑥̅ 𝑖 , 𝑦 ̅ 𝑗 ), and obeying to the position

(𝑖, 𝑗) ∈ [2, 𝑁 𝑥 -1] × [2, 𝑁 𝑦 -1].
Therefore, radiative parameters desired are calculated for all point 𝑀 𝑖𝑗 belonging the participating medium, such that:

𝑀 𝑖𝑗 ( (𝑖 - 3 2
)

𝐻 𝑥 𝑁 𝑥 -2 (𝑗 - 3 2
)

𝐻 𝑦 𝑁 𝑦 -2 ) (1)
The radiation incoming at point 𝑀 𝑖𝑗 appearing on Fig. 4b, originates from west, north and partially south and east boundary surfaces of the participating medium, added to the one of west and north boundary surfaces of the inner square cavity. "When there is a spectral dependence of radiative properties, like in gaseous participating medium, absorption coefficients can be reordered in a set of absorbing bands, and coupled to the radiative transfer equation. Therefore, semi-analytical method used to solve the problem becomes not appropriate because of the present irregular geometry compared to numerical methods [START_REF] Sakami | Analysis of radiative heat transfer in complex twodimensional enclosures with obstacles using the modified discrete ordinate method[END_REF][START_REF] Chai | Treatment of irregular geometries using a Cartesian coordinates finite volume radiation heat transfer procedure[END_REF][START_REF] Talukdar | Discrete transfer method with the concept of blocked-off region for irregular geometries[END_REF], which can be directly involved."

Governing equations

Energy transfer by radiation in a gray semi-transparent medium, which absorbs, emits, but does not scatter radiation, at steady state condition is represented by a differential equation in the plane (𝑥, 𝑦), as: , 𝜋] could be deduced by simple rotation and translation.

𝑠𝑖𝑛𝜃𝑐𝑜𝑠𝜑 𝜕𝐼(𝑥,
Hence, incident radiation originating from all the directions of propagation following Ω ⃗⃗⃗ , is set as :

𝐺(𝑥, 𝑦) = 2 ∫ ∫ 𝐼(𝑥, 𝑦, 𝜃, 𝜑) 𝜑=2𝜋 𝜑=0 𝜃= 𝜋 2 𝜃=0 𝑠𝑖𝑛𝜃𝑑𝜃𝑑𝜑 (4)
The radiative flux vector , is also classically given by :

𝑞 ⃗ 𝑟 (𝑥, 𝑦) = 2 ∫ ∫ 𝐼(𝑥, 𝑦, 𝜃, 𝜑) 𝜑=2𝜋 𝜑=0 𝜃= 𝜋 2 𝜃=0 𝑠𝑖𝑛𝜃 2 ( 𝑐𝑜𝑠𝜑 𝑠𝑖𝑛𝜑 0 ) 𝑑𝜃𝑑𝜑 (5)
Consequently, temperature field is iteratively deduced from the divergence relation of radiative flux vector at equilbrium state conditions by :

∇ ⃗ ⃗⃗ . 𝑞 ⃗ 𝑟 (𝑥, 𝑦) = 𝑘 𝑎 (4𝜋𝐼 𝑏 (𝑇(𝑥, 𝑦)) -𝐺(𝑥, 𝑦)) (6
𝐸 𝑆 𝑊 𝑁 1 𝜑 𝑁 1 𝜑 𝑁 2 𝜑 𝑊 2 𝜑 𝑊 1 𝜑 𝑆 1 𝜑 𝑆 2 𝜑 𝐸 1 𝜑 𝐸 2

Incident radiation

Let evaluate quantitatively the incident radiation, which exists within the domain 𝐷 . The first task consists to establish the respective equations of lines which delimit each sub-surface where incoming radiation 𝐼(𝑥, 𝑦, 𝜃, 𝜑) is find; after integrating it from all directions of Ω ⃗⃗⃗ , in order to obtain the incident radiation 𝐺(𝑥, 𝑦) in the domain 𝐷 , the following relations are obtained and shown as:

 Sub-surface 7 : (𝑥, 𝑦) ∈ ]0, 𝐻-ℎ 2 [ × ] 𝐻+ℎ 2 , 𝐻[
Incoming radiation viewed from north surface of the participating medium is projected along 𝑒 ⃗ 𝑧 on the (𝑥, 𝑦) plane following the curvilinear abscissa 𝛿(𝑥, 𝑦, 𝜃, 𝜑) is obtained from Eq.( 3) by: 𝐼 𝑁 1 (7) (𝑥, 𝑦, 𝜃, 𝜑) =

𝜎 𝐵 𝑇 𝑁 1 4 𝜋 . 𝑒 -𝑘 𝑎 { 𝐻-𝑦 𝑠𝑖𝑛𝜑𝑠𝑖𝑛𝜃 } + 𝑘 𝑎. 𝜎 𝐵 𝜋 ∫ 𝑇 4 (𝑥 ′ , 𝑦 ′ ). 𝑒 -𝑘 𝑎 { 𝑦 ′ -𝑦 𝑠𝑖𝑛𝜑𝑠𝑖𝑛𝜃 } 𝑑𝑦 ′ 𝑠𝑖𝑛𝜑𝑠𝑖𝑛𝜃 𝑦 ′ =𝐻 𝑦 ′ =𝑦 (7)
with, 𝑥 ′ = 𝑥 ± (𝑦 ′ -𝑦)𝑡𝑎𝑛𝜑.

From all directions of vector Ω ⃗⃗⃗ , radiation covers 𝜑 𝑁 1 = 𝑚𝑒𝑠(𝑀 𝑖𝑗 𝑁 𝐸 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ , 𝑀 𝑖𝑗 𝑁 1 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ ), and 𝜑 𝑁 2 = 𝑚𝑒𝑠(𝑀 𝑖𝑗 𝑁 1 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ , 𝑀 𝑖𝑗 𝑁 𝑊 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ ), therefore we can write:

𝐺 𝑁 1 (7) (𝑥, 𝑦) = 2 ∫ ∫ 𝐼 𝑁 1 (𝑥, 𝑦, 𝜃, 𝜑)𝑠𝑖𝑛𝜃𝑑𝜃𝑑𝜑 𝜃= 𝜋 2 𝜃=0 𝜑= 𝜋 2 +𝜑 𝑁 2 𝜑= 𝜋 2 -𝜑 𝑁 1 (8)
with,

𝜑 𝑁 1 = 𝑡𝑎𝑛 -1 { 𝐻-𝑥 𝐻-𝑦 } and 𝜑 𝑁 2 = 𝑡𝑎𝑛 -1 { 𝑥 𝐻-𝑦 } .
Thanks to Eq. ( 7), and a change of variables in term of 𝜑 , with two successive angular integrations similar to the one developed in [START_REF] Chen | Numerical analysis on the radiation-convection coupled heat transfer in an open-cell foam filled annulus[END_REF], solution of the incident radiation is now given by:

𝐺 𝑁 1 (7) (𝑥, 𝑦) =

2𝜎 𝐵 𝑇 𝑁 1 4 𝜋 {𝐵 𝑖𝑠 2 (𝑘 𝑎 (𝐻 -𝑦), 𝜑 𝑁 1 ) + 𝐵 𝑖𝑠 2 (𝑘 𝑎 (𝐻 -𝑦), 𝜑 𝑁 2 )} + 2𝑘 𝑎. 𝜎 𝐵 𝜋 ∫ ∫ 𝑇 4 (𝑥 ′ = 𝑥 + (𝑦 ′ -𝑦)𝑡𝑎𝑛𝜑, 𝑦 ′ ) 𝑦 ′ =𝐻 𝑦 ′ =𝑦 𝜑=𝜑 𝑁 1 𝜑=0 𝐾 𝑖 1 { 𝑘 𝑎 (𝑦 ′ -𝑦) 𝑐𝑜𝑠𝜑 } 𝑑𝑦 ′ 𝑑𝜑 𝑐𝑜𝑠𝜑 + 2𝑘 𝑎. 𝜎 𝐵 𝜋 ∫ ∫ 𝑇 4 (𝑥 ′ = 𝑥 -(𝑦 ′ -𝑦)𝑡𝑎𝑛𝜑, 𝑦 ′ ) 𝑦 ′ =𝐻 𝑦 ′ =𝑦 𝜑=𝜑 𝑁 2 𝜑=0 𝐾 𝑖 1 { 𝑘 𝑎 (𝑦 ′ -𝑦) 𝑐𝑜𝑠𝜑 } 𝑑𝑦 ′ 𝑑𝜑 𝑐𝑜𝑠𝜑 (9)
where,

𝐾 𝑖 𝑛 (𝑢) = ∫ 𝑒 -𝑢 𝑠𝑖𝑛𝜃( 𝑠𝑖𝑛𝜃) 𝑛-1 𝑑𝜃 𝜃= 𝜋 2 𝜃=0
, 𝑛 ∈ ℕ, 𝑢 ∈ ℝ + (10) and,

𝐵 𝑖𝑠 𝑛 (𝑢, 𝜃) = ∫ 𝐾 𝑖 𝑛 ( 𝑢 𝑐𝑜𝑠𝜑 ) (𝑐𝑜𝑠𝜑) 𝑛-2 𝜑=𝜃 𝜑=0 𝑑𝜑, 𝑛 ∈ ℕ, 𝑢 ∈ ℝ + (11)
represent respectively a set of Bickley's Naylor functions, and Altaç modified Bessel functions [START_REF] Altaç | Integrals Involving Bickley and Bessel Function in radiative transfer and generalized exponential integral function[END_REF]. For numerical computation, Eq. ( 9) is considered on its discretized form. Hence, at each center cell, this incident radiation propagates following some conditions:

{ 𝑦 𝑗 ′ = 𝑦 𝑗 + (𝐻 -𝑦 𝑗 )𝛿 𝑚 𝑥 𝑖 ′ = 𝑥 𝑖 ± (𝐻 -𝑦 𝑗 )𝛿 𝑚 𝑡𝑎𝑛𝜑 (12)
where, 𝛿 𝑚 ∈ [0,1], represents the quadrature abscissa, 𝑚 the abscissa number, such 𝑚 ∈ [1, 𝑀] and 𝑀 is total number of quadrature set. For instance, during the rays pathway, there exist particular positions of points 𝑀 𝑖𝑗 to follow, and that are directly related to the couple point

(𝑝, 𝑞) ∈ ℕ 2 , such that [𝑥 𝑝 - ∆𝑥 2 , 𝑥 𝑝 + ∆𝑥 2 ] and 𝑦 𝑗 ′ ∈ [𝑦 𝑞 - ∆𝑦 2 , 𝑦 𝑞 + ∆𝑦 2 ].
It helps to deduce easily their position in the reference coordinate, by [(𝑝 -

) ∆𝑥 , (𝑞 -

) ∆𝑦],

where :

{ 𝑝 ≤ 𝑖 + 1 2 ± (𝑁 𝑦 -𝑗 - 1 2 ) 𝛿 𝑚 𝑡𝑎𝑛𝜑 𝑞 ≤ 𝑗 + 1 2 + (𝑁 𝑦 -𝑗 - 1 2
) 𝛿 𝑚 [START_REF] Pasini | Numerical experiments on the application of the diamond scheme to sets of discrete directions obtained from a random numbers generator[END_REF] In the same order, azimuth angle 𝜑 is also discretized by the use of Gauss quadrature such that:

𝑑𝜑 = (𝜑 𝑚𝑎𝑥 -𝜑 𝑚𝑖𝑛 )𝛽 𝑙 , 𝑙 ∈ {1,2, … , 𝑁 𝜑 } , 𝜑 𝑙 ∈ 0,1 (14)
where, 𝜑 𝑚𝑎𝑥 and 𝜑 𝑚𝑖𝑛 are respectively the maximun and minimum boundary integrals to discretize. 𝛽 𝑙 is the angular abscissa, 𝑁 𝜑 is the number of Gauss quadrature used to approximate respective integrals in Eq. ( 9). Let define a constant, 𝛼 = 2𝜎 𝐵 𝜋 in the problem, finally, the discretized form of incident radiation from north location having boundary surface 𝛤 𝑁 1 becomes in terms of center cell coordinates (𝑥̅ 𝑖 , 𝑦 ̅ 𝑗 ), working with the square cells, ∆𝑥 = ∆𝑦 :

∀(𝑖, 𝑗) ∈ [2, 𝐸 { 𝐻-ℎ 2∆𝑥 }] × [𝐸 { 𝐻+ℎ 2∆𝑦 } , 𝑁 𝑦 ],
𝐺 𝑁 1 (7) (𝑖, 𝑗) = 𝛼𝑇 𝑁 1 4

) (15) with, 𝜑 ̅ 𝑁 1 (𝑖,𝑗) = 𝑡𝑎𝑛 -1 { 𝑢 𝑖 𝑢 𝑗 }, 𝜑 ̅ 𝑁 2 (𝑖,𝑗) = 𝑡𝑎𝑛 -1 { 𝑣 𝑖 𝑢 𝑗 } , 𝜑 𝑙 1 = 𝛽 𝑙 × 𝜑 ̅ 𝑁 1 (𝑖,𝑗) and 𝜑 𝑙 2 = 𝛽 𝑙 × 𝜑 ̅ 𝑁 2 (𝑖,𝑗) , 𝑢 𝑗 = (𝑁 𝑦 -𝑗 + 3 2 ) ∆𝑦 , 𝑢 𝑖 = (𝑁 𝑥 -𝑖 + 3 2 ) ∆𝑥 , 𝑣 𝑖 = (𝑖 - 3 2
) ∆𝑥, and 𝐸 𝑢 denoted an integer part of a real 𝑢 ∈ ℝ.

The same calculation is performed for incident radiation 𝐺 𝑊 1 (7) from west location delimited by the angle 𝑚𝑒𝑠(𝑀 𝑖𝑗 𝑁 𝑊 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ , 𝑀 𝑖𝑗 𝑆 𝑊 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ ) = 𝜑 𝑊 1 + 𝜑 𝑊 2 , where the boundary surface is 𝛤 𝑊 1 , then:

𝐺 𝑊 1 (7) (𝑖, 𝑗) = 𝛼𝑇 𝑊 1 4

with 𝑗) and 𝜑 𝑙 4 = 𝛽 𝑙 × 𝜑 ̅ 𝑊 2 (𝑖,𝑗) , and

, 𝜑 ̅ 𝑊 1 (𝑖,𝑗) = 𝑡𝑎𝑛 -1 { 𝑣 𝑗 𝑣 𝑖 } , 𝜑 ̅ 𝑊 2 (𝑖,𝑗) = 𝑡𝑎𝑛 -1 { 𝑢 𝑖 𝑣 𝑖 }, 𝜑 𝑙 3 = 𝛽 𝑙 × 𝜑 ̅ 𝑊 1 (𝑖,
𝑣 𝑗 = (𝑗 - 3 2
) ∆𝑦.

The exact semi-analytical expressions of incident radiation denoted 𝐺 𝑆 1 (7) are calculated under the aperture 𝜑 𝑆 such that

𝜑 𝑆 = 𝑚𝑒𝑠(𝑀 𝑖𝑗 𝑆 𝑊 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ , 𝑀 𝑖𝑗 𝑆 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ ) = 𝜑 𝑆 1 + 𝜑 𝑆 2 .
similarly those of 𝐺 𝐸 1 (7) , 𝐺 𝑊 2 (7) and 𝐺 𝑁 2 (7) are respectively performed under the apertures represented in the annex 1. Finally, within sub-surface 7, radiation intensity 𝐺 (7) (𝑥̅ 𝑖 , 𝑦 ̅ 𝑗 ) is the contribution of rays following all apertures illustrated before, hence: 𝐺 (7) (𝑖, 𝑗) = 𝐺 𝑁 1 (7) (𝑖, 𝑗) + 𝐺 𝑊 1 (7) (𝑖, 𝑗) + 𝐺 𝑆 1 (7) (𝑖, 𝑗) + 𝐺 𝐸 1 (7) (𝑖, 𝑗) + 𝐺 𝑊 2 (7) (𝑖, 𝑗) + 𝐺 𝑁 2 (7) (𝑖, 𝑗)

In the other sub-surfaces 2,12,17 shown on Fig. 4a, incident radiation are deduced by translation of indices using expressions already established, then:

𝐺 (2) (𝑖, 𝑗) = 𝐺 (7) ((𝑁 𝑥 -𝑖 +

) , (𝑗 -

)) ; 𝐺 (12) (𝑖, 𝑗) = 𝐺 (7) ((𝑖 - )) and 𝐺 (17) (𝑖, 𝑗) = 𝐺 (7) ((𝑁 𝑥 -𝑖 + )).

Therefore, the other exact expressions of incident radiations corresponding to all sub-surfaces 𝜁 𝑛 , 𝑛 ∈ [1,20]\ 7 follow a very similar method of evaluation, but computed and presented in the tables. Finally the process end, regardeless of :

(𝑖, 𝑗) ∈ (]2, 𝑁 𝑥 -2[ × ]2, 𝑁 𝑦 -2[)\ (]𝐸 { (𝐻-ℎ) 2∆𝑥 } , 𝐸 { (𝐻+ℎ) 2∆𝑥 }[ × ]𝐸 { (𝐻-ℎ) 2∆𝑦 } , 𝐸 { (𝐻+ℎ) 2∆𝑦 }[),
where 𝐸(𝑢) denotes the integer part of a real 𝑢 ∈ ℝ, the exact incident radiation expected at Eq. ( 4) is the given by:

𝐺(𝑖, 𝑗) = ∑ 𝐺 (𝑘) (𝑖, 𝑗) 20 𝑘=1 (18)
The boundary conditions related to the participating medium are as follow:

for 𝑥̅ 𝑖 = 0, and 𝑦 ̅ 𝑗 ∈ ]0, 𝐻[ , we have 𝐺(𝑖, 𝑗) = 𝜎 𝐵 𝑇 𝑊 1 4 ; for 𝑥̅ 𝑖 = 𝐻 and 𝑦 ̅ 𝑗 ∈ [0, 𝐻], 𝐺(𝑖, 𝑗) = 𝜎 𝐵 𝑇 𝐸 1 4 ; for 𝑥̅ 𝑖 ∈ ]0, 𝐻[ and 𝑦 ̅ 𝑗 = 0, we have 𝐺(𝑖, 𝑗) = 𝜎 𝐵 𝑇 𝑆 1 4 ; for 𝑥̅ 𝑖 ∈ ]0, 𝐻[and 𝑦 ̅ 𝑗 = 𝐻, 𝐺(𝑖, 𝑗) = 𝜎 𝐵 𝑇 𝑁 1 4 .

At the surfaces of inner square cavity, they are set such that : , the result of 𝐺(𝑖, 𝑗) = 𝜎 𝐵 𝑇 𝑁 2 4 .

The results of exact semi-analytical expressions of incident radiation are then computed numerically for better proof.

Temperature field

Temperature field 𝑇(𝑖, 𝑗) is deduced From Eq. ( 6) at radiative equilibrium, using the relation ∇ ⃗ ⃗⃗ . 𝑞 ⃗ 𝑟 (𝑖, 𝑗) = 0. Thanks to Gauss quadrature used to discretize the space-angular integrals, which define the incident radiation within the participating medium 𝐷 . Once, solved by iteration process the resulting equation is:

𝑘 𝑎 (4𝜋𝐼 𝑏 (𝑇(𝑖, 𝑗)) -∑ 𝐺 (𝑘) (𝑖, 𝑗) 20 𝑘=1
) = 0 , 𝑘 𝑎 ≠ 0

Consequently, the resulting temperature field is used in the process of determination the radiative flux field.

Radiative flux field

Following a similar procedure as shown in section 2, radiative flux vector calculated previously under apertures shown on Fig. 5 and delimited by northern boundary surface 𝛤 𝑁 1 belonging the enclosure, is:

𝑞 ⃗ 𝑟 𝑁 1 = 2 ∫ ∫ 𝐼 𝑁 1 (𝑥, 𝑦, 𝜃, 𝜑)(𝑠𝑖𝑛𝜃) 2 ( 𝑐𝑜𝑠𝜑 𝑠𝑖𝑛𝜑 ) 𝑑𝜃𝑑𝜑 𝜃= 𝜋 2 𝜃=0 𝜑= 𝜋 2 +𝜑 𝑁 2 𝜑= 𝜋 2 -𝜑 𝑁 1 (20)
For several sub-surfaces of incoming radiation 𝐼 𝑁 1 (𝑥, 𝑦, 𝜃, 𝜑) is calculated, in order to deduce components of radiative flux vector. Then, within sub-surface 7, the relations below are obtained using similar calculations as the one presented in the reference [START_REF] Djeumegni | Modeling of radiative heat transfer in a gray semi-transparent medium with internal fluid cavity limited by black boundary surfaces[END_REF].

 ∀(𝑖, 𝑗) ∈ [2, 𝐸 { 𝐻-ℎ 2∆𝑥 }] × [𝐸 { 𝐻+ℎ 2∆𝑦 } , 𝑁 𝑦 ],
Following aperture 𝜑 𝑁 = 𝑚𝑒𝑠(𝑀 𝑖𝑗 𝑁 𝐸 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ , 𝑀 𝑖𝑗 𝑁 1 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ ) and delimited by boundary surface 𝛤 𝑁 1 , radiative flux vector is set by: 𝑞 ⃗ 𝑟,𝑁 1 (7) = 𝛼𝑇 𝑁 1 4 ( 𝐶 𝑖𝑠 3 (𝑘 𝑎 𝑢 𝑗 , 𝜑 ̅ 𝑁 1 (𝑖,𝑗) ) -𝐶 𝑖𝑠

where, 𝐶 𝑖𝑠 𝑛 , 𝑛 ∈ ℕ is a modified Bessel function [START_REF] Talukdar | Discrete transfer method with the concept of blocked-off region for irregular geometries[END_REF] set by:

𝐶 𝑖𝑠 𝑛 (𝑢, 𝜃) = ∫ 𝐾 𝑖 𝑛 ( 𝑢 𝑐𝑜𝑠𝜑 ) (𝑐𝑜𝑠𝜑) 𝑛-3 𝑠𝑖𝑛𝜑 𝜑=𝜃 𝜑=0 𝑑𝜑 (22)
Radiative flux within sub-surface 7, following aperture 𝜑 𝑊 = 𝑚𝑒𝑠(𝑀 𝑖𝑗 𝑁 𝑊 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ , 𝑀 𝑖𝑗 𝑆 𝑊 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ ) and delimited by boundary surface 𝛤 𝑊 1 is also deduced as:

𝑞 ⃗ 𝑟,𝑊 1 (7) = 𝛼𝑇 𝑊 1 4 ( -𝐵 𝑖𝑠

In Eq. [START_REF] Wilson | Radiative transfer in a two-dimensional rectangular annulus medium[END_REF][START_REF] Le Dez | Radiative transfer in a semi-transparent medium enclosed in a cylindrical annulus[END_REF] Consequently, as from the one of incident radiation, the rest of exact expressions of radiative flux (𝑞 ⃗ 𝑟,𝑆 1 (7) , 𝑞 ⃗ 𝑟,𝐸 1 (7) , 𝑞 ⃗ 𝑟,𝑊 2 (7) , and 𝑞 ⃗ 𝑟,𝑁 2 (7) ) that evaluate quantitatively what is find within sub-surface 7 are shown in the annex 2. The others exact semi-analytical equations are not given there for conciseness, since the method of modelling is similarly the same. The sum of contribution of radiative flux within sub-surface 7 yields to:

𝑞 ⃗ 𝑟 (7) (𝑖, 𝑗) = 𝑞 ⃗ 𝑟,𝑁 1 (7) (𝑖, 𝑗) + 𝑞 ⃗ 𝑟,𝑊 1 (7) (𝑖, 𝑗) + 𝑞 ⃗ 𝑟,𝑆 1 (7) (𝑖, 𝑗) + 𝑞 ⃗ 𝑟,𝐸 1 (7) (𝑖, 𝑗) + 𝑞 ⃗ 𝑟,𝑊 2 (7) (𝑖, 𝑗) + 𝑞 ⃗ 𝑟,𝑁 2 (7) (𝑖, 𝑗) (

Once, Eq. ( 23) enables to deduce some of radiative flux components from sub-surfaces 7, by translation of indices 𝑖 and 𝑗. However, it generates:

𝑞 ⃗ 𝑟 (2) (𝑖, 𝑗) = 𝑞 ⃗ 𝑟 (7) ((𝑁 𝑥 -𝑖 +

) , (𝑗 -

)); 𝑞 ⃗ 𝑟 (12) (𝑖, 𝑗) = 𝑞 ⃗ 𝑟 (7) ((𝑖 -

) , (𝑁 𝑗 -𝑗 +

)) and 𝑞 ⃗ 𝑟 (17) (𝑖, 𝑗) = 𝑞 ⃗ 𝑟 (7) ((𝑁 𝑥 -𝑖 + )).

while,

(𝑖, 𝑗) ∈ (]2, 𝑁 𝑥 -2[ × ]2, 𝑁 𝑦 -2[)\ (]𝐸 { (𝐻-ℎ) 2∆𝑥 } , 𝐸 { (𝐻+ℎ) 2∆𝑥 }[ × ]𝐸 { (𝐻-ℎ) 2∆𝑦 } , 𝐸 { (𝐻+ℎ) 2∆𝑦 }[),
the vector form of radiative flux is obtained by:

𝑞 ⃗ 𝑟 (𝑖, 𝑗) = ∑ 𝑞 ⃗ 𝑟 (𝑘) (𝑖, 𝑗) 20 𝑘=1 (25)
In the same way, radiative flux at boundary surfaces of enclosure are given from the following relations:

if

𝑥̅ 𝑖 = 0,

Numerical results

This section is devoted to the presentation of numerical results obtained from exact semianalytical solutions, previously presented. The results of radiative quantities correspond to the ratio modified Bessel functions are computed by using a specific FORTRAN algorithm similar to [START_REF]Algorithm 609, a portable FORTRAN subroutine for the Bickley functions Kin (x)[END_REF] and designed for the specific problem. Therefore, accuracy of the present model based on ray tracing method is firstly checked with a reference data [START_REF] Crosbie | Radiative transfer in a two-dimensional rectangular medium exposed to diffuse radiation[END_REF] for a simple geometry, by comparing values of temperature field at line = 𝐻/2 , when the size of inner square cavity is set to zero. Secondly, the local radiative flux is also plotted for complex geometry and validated with some benchmark results [START_REF] Sànchez | Surface Radiation Exchange for Two-Dimensional Rectangular Enclosures Using the Discrete -Ordinates Method[END_REF][START_REF]Algorithm 609, a portable FORTRAN subroutine for the Bickley functions Kin (x)[END_REF]. Thirdly, tabulated data and various simulations are also performed for several sizes of the inner square cavity ℎ , and for different values of the optical thickness 𝜏. South boundary surface of the participating medium is hot (𝑇 𝐸 1 = 𝑇 𝑊 1 = 𝑇 𝑁 1 = 0, 𝑇 𝑆 1 ≠ 0) and (𝑇 𝐸 2 = 𝑇 𝑊 2 = 𝑇 𝑁 2 = 𝑇 𝑆 2 = 0). The result is compared when ℎ = 0, with the reference data established by Crosbie and Schrenker [START_REF] Altaç | Integrals Involving Bickley and Bessel Function in radiative transfer and generalized exponential integral function[END_REF], when the optical thickness is 𝜏 = 0.25.

The choice of many cells is used for better accuracy in the results expected. While, at least (251 × 251 ) cells are implemented, and show exact results for small and average optical thickness (𝜏 = 𝑘 𝑎 𝐻), but is not really the same in terms of large optical thickness [START_REF]Algorithm 609, a portable FORTRAN subroutine for the Bickley functions Kin (x)[END_REF]. The Fig. 6 , where 𝑇 * = 0.25. This exact value is also obtained on the reference benchmark data [START_REF] Crosbie | Radiative transfer in a two-dimensional rectangular medium exposed to diffuse radiation[END_REF]. Practically, instead in the case of solid (foam, glass, etc..), fluid or gas participating media attention is focussed on small and average optical thickness; consequently the second integral part of Eq. ( 3) disappears . However, it makes the semi-transparent medium behaves like a none internal emitting source of radiation, through which rays propagate from boundary surfaces. From these approximations, previous analytical equations could be simplified and nondimensional temperature 𝑇 * becomes proportional to radiation intensity 𝐺 * . When the optical thickness becomes large, methods like Rosseland approximation are adopted by solving 𝑑𝑖𝑣(𝑘 𝑅 𝑔𝑟𝑎𝑑 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ 𝑇) = 0, where 𝑘 𝑅 denotes the Rosseland radiative conductivity, shown by the relation: 𝑘 𝑅 = (16𝜎 𝐵 /3𝑘 𝑎 )𝑇 3 . Eq. ( 36) is computed at radiative equilibrium after several iterations and shows a smooth correlation with literature, but carrying also a little errors

|∆𝜀| = |𝑇 𝐶𝑟𝑜𝑠𝑏𝑖𝑒 4 -𝑇 𝑝𝑟𝑒𝑠𝑒𝑛𝑡 4
|, no longer to |∆𝜀| ≤ 0.2%, and allows us to have a good confidence in the proposed method. Although most of innovative engineering problems dealing with heat transfer within semi-transparent media, applied to complex geometries like the present, convection mode dealt with [START_REF] Mezrhab | Computation of combined naturalconvection and radiation heat transfer in a cavity having a square body at its center[END_REF][START_REF] Lari | Combined heat transfer of radiation and natural convection in a square cavity containing participating gases[END_REF][START_REF] Saravanan | Coupled thermal radiation and natural convection heat transfer in a cavity with a heated plate inside[END_REF]. Such being the case, it is firstly very important to control with accuracy the behavior of radiation heat transfer mode, before coupling it to Navier's Stokes equations; this is the reason why the actual work is proposed.

Complex geometry

The complex geometry is Configured quite similar to Fig. 1a, and has already been investigated by Sànchez and Smith [START_REF] Sànchez | Surface Radiation Exchange for Two-Dimensional Rectangular Enclosures Using the Discrete -Ordinates Method[END_REF], followed by Chai et al [START_REF] Chai | Treatment of irregular geometries using a Cartesiancoordinates-based discrete-ordinates method[END_REF], and Talukdar [START_REF] Talukdar | Discrete transfer method with the concept of blocked-off region for irregular geometries[END_REF]. In this case, all surfaces of the participating medium are black with length 𝐻 = 1.0 𝑚. Following Fig. 1a, 𝑇 𝑊 1 = 310𝐾 and the rest of boundary surfaces 𝛤 𝐸 1 , 𝛤 𝑁 1 , 𝛤 𝑆 1 are maintained at isothermal temperature 𝑇 = 300𝐾 for the participating medium. Only the south boundary surface of the inner square cavity 𝛤 𝑆 2 is imposed at temperature 𝑇 𝑆 2 = 300𝐾 and the other boundary surfaces are cold. Therefore, the result of local radiative flux in term of non-dimensional length 𝜍 is displayed on Fig. 7 bellow and validated with the radiation/irradiation method (RIM) from Siegel and Howell [START_REF] Siegel | thermal radiation heat transfer[END_REF], and DOM proposed by Sànchez and Smith [START_REF] Sànchez | Surface Radiation Exchange for Two-Dimensional Rectangular Enclosures Using the Discrete -Ordinates Method[END_REF]. 𝑞 𝑝𝑟𝑒𝑠𝑒𝑛𝑡 )/𝑞 𝑟𝑒𝑓 ] gives a value less than 0.1 percent with RIM, whereas one obtain a value near to 0.2 percent compared with DOM. It has already been proven that, more 𝑀 is higher, more the result is accurate when using DOM, and one can suitably deduce from the work proposed by Sànchez and Smith [7] between RIM and DOM, a similar conclusion between the present method and DOM. Furthermore, computation of heat flux using the present method remains also in good agreement with the one of discrete transfer method (DTM) already proposed by Talukdar [START_REF] Talukdar | Discrete transfer method with the concept of blocked-off region for irregular geometries[END_REF], but is not plotted there for concisness. Accuracy on the present method with RIM derived on the exact method of calculation the integral form of radiative transfer equation.

Tabulated data of radiative quantities

Tables 1 and5 show for optical thickness 𝜏 = 0.25, the discrete values of the non-dimensional temperature 𝑇 * at the southwest side of the participating medium, because the geometry is symmetric. It is shown that, there is no heat transfer with the inner square, that means (𝑥, 𝑦) ∈ at the southwest side of the participating medium is deduced from exact semi-analytical

Numerical simulation of radiative quantities

Numerical simulations of non-dimensional temperature and radiative flux fields are displayed for various values of ratio 𝐻/ℎ that characterize the geometry of the semi-transparent medium and performed from lower to higher values of optical thickness. It is noticeable that, solve the RTE for this complex geometry, is more easy to deal with skewed non-orthogonal meshes grid, but in this work, because of the choice of exact method of resolution based on ray-tracing method, orthogonal meshes are adapted for it. All the simulations below are performed with 251 × 251 cells. Fig. 7(a-c) represent a non-dimensional temperature fields 𝑇 * within the semi-transparent medium when there is no inner cavity (ℎ = 0), and only south boundary surface 𝛤 𝑆 1 is hot. Optical thickness used, varies from lower to higher values (𝜏 = 0.01, 𝜏 = 1, 𝜏 = 5). Simulation results remains in good agreement with literature [START_REF] Le Dez | Radiative heat transfer in a semi-transparent medium enclosed in a two-dimensional square cavity[END_REF].When 𝜏 = 0.01, it corresponds to a very low absorption coefficients 𝑘 𝑎 , and temperature distribution reaches quickly a north boundary surface of the participating medium with a value of 𝑇 * evaluated between 0.1 and 0.2. When 𝜏 = 1, 𝑇 * is deduced between 0.0 and 0.1, but when 𝜏 = 5, temperature in the semi-transparent medium is almost zero, except near the south boundary surface. Therefore, higher the optical thickness, more the absorption in the medium. Those of Fig. 8(a-c) show the corresponding ycomponents of non-dimensional radiative flux 𝑞 𝑦 * = 𝑞/𝜎 𝐵 𝑇 4 performed in conditions similar to Fig. 7 (a-c), and a radiative flux decreases suddenly when optical thickness increases. Simulations of radiative flux on Fig. 8(a-c) displayed above, show radiation propagates following parabolic and hyperbolic shapes. On Fig. 7(a-c), when temperature decreases with ray path length, radiative flux does the opposite. Fig. 9(a-c) and Fig. 10(a-c) describe also respectively 𝑇 * and 𝑞 𝑦 * when both south and north boundary surfaces 𝛤 𝑆 1 , 𝛤 𝑁 1 are hot. These particular simulations are useful when there is a need to study a coupling with natural convection of participating gases as show in [START_REF] Mezrhab | Computation of combined naturalconvection and radiation heat transfer in a cavity having a square body at its center[END_REF]. On Fig. 7(a), symmetry on ray propagation is shown and at the medium line 𝑥 = 𝐻/2, one have 𝑇 * = 0.25 already proven in the literature. On Fig. 8(a-c), radiative flux is ineffective, such 𝑞 𝑦 * = 0 at the same line 𝑥 = 𝐻/2 and there is no radiative flux following x-axis (𝑞 𝑥 * = 0). On Fig. 11(a-c) and Fig. 12(a-c), one have depicted the simulations of radiation when the respective south and east boundary surfaces of the participating medium are hot, with a broad range of optical thicnesses from 𝜏 = 0.01,1,5 . While temperature propagates in the semitransparent medium, the inner square cavity blocks rays from south to north boundary surfaces for Fig. 11(a-c) surfaces and from west to east boundary surfaces for Fig. 12(a-c). In subsurfaces 5 already presented in Fig. 4a, there is no ray that reaches from south boundary surface to this place, then radiation which passes through results only from emission of the medium. Consequently the effect of blockage of radiation deviates the trajectory of hot particles and The Fig. 13(a-c) and Fig. 14(a-c) present simulations of non-dimensional temperature and radiative flux fields when the south and north are simultaneously hot. One can remark that, temperature field near the east, north boundary surfaces on theses Fig. 13 (a-c), and Fig. 14 (ac) is really symmetric, but non-uniform as compared from the one displayed on Fig. 9(a-c) and Fig. 10 (a-c). This is because the inner square cavity does not allow rays propagate from bottom to top of the semi-transparent medium. A similar phenomenon can be observed when 𝛤 𝑊 1 and 𝛤 𝐸 1 are also simultaneously hot and the other surfaces are set to be cold. When 𝐻 = 3ℎ, the triangular sub-surfaces 5,6,15,16 from Fig. 2a are well observed on the suitable simulations shown on Fig. 15(a-c). The results of non-dimensional temperature presented there, justify and strengthen the geometrical analysis developed at section 2.1. One also argue that, more the size of the inner cavity, higher the blockage of rays in the semitransparent medium. Similarly, to the case of simple geometry, if absorption coefficient increases within the participating medium optical thickness directly increases, then radiation will quickly be attenuated; it will be an important loss of energy transfer; it is the case of macroporous media (fibers, foam, etc…). The non-participating medium is well represented through an empty white square, which does not participates; practically it can carry a fluid or a gas at steady conditions, or by another solid medium. There is symmetry in temperature distribution, but also some strange points where radiation changes suddenly from one neighboring point to another; this is because the inner square cavity blocks the rays on both sides of the middle lines When, the size of the inner square cavity becomes very small like 𝐻 = 8ℎ compared to the cases of 𝐻 = 3ℎ or 𝐻 = (2 + √5)ℎ, obsctruction of rays becomes ineffective. Then, when optical thickness is small, radiation propagates more easily in the medium and there is not a significant loss of energy. Moreover, the influence of optical thickness remains unchanged, it is shown on Fig. 17(a-b) and Fig. 18(a-b). The main difference when the size of the inner square cavity becomes very small as compare to the one of large size is firstly based on their respective geometries. Another approach used on simulations is to make participate the boundary surfaces of the inner square cavity 𝛤 𝑊 2 , 𝛤 𝐸 2 and 𝛤 𝑆 2 , 𝛤 𝑁 2 , and the rest boundary surfaces of the participating medium are maintained cold. This is to show how internal emission within the semi-transparent medium can be observed (Fig. 21(a-c) to Fig. 24(a-c)), by variating the size of inner cavity concerned for an range of optical thicknesses. One can mention that, when the inner square cavity becomes very small, the values of temperature within the semi-transparent medium are almost equal to zero, but radiation propagates following the same shape.

as from the figures illustrated on the geometrical analysis, simulations proposed confirm the idea described in sub-section 2.1, because the radiative quantities have different analytical expressions within each sub-surfaces of the participating medium. The smaller the size of the internal cavity, the more internal sub-surfaces are created.This is a new method used to describe radiative transfer based of exact method using ray tracing and that will be performed later by adding another conditions and hypothesis.

Conclusion and remarks

A deep semi-analytical study, followed by numerical simulations within a gray semitransparent medium having a centered inner square cavity has been presented. Temperature, incident radiation, and radiative flux vector fields have been evaluated. The method used there, was based on ray tracing method adapted for simple and complex geometries, and followed by the used of specific functions which calculate integrals solutions of RTE. After a rigourous analysis of the five possible and different cases of predicting radiative transfer in the participating medium concerned, exact discretized forms of equations governing heat transfer have been established. Thanks to Gauss quadrature used for this purpose. Numerical results of temperature have been compared with the one of literature when the size of the inner square is reduced to zero, and when it is different to zero and there is good agreement. Some data of nondimensional temperature and incident radiation have been tabulated and proposed as benchmark results, then several simulations have also been performed for the same radiative quantities at different broad ranges of optical thicknesses. The method produces smooth temperature and flux fields without oscillatory errors, realized in a satisfactory computation time. The main difficulties encountered were based on how to analyze the irregular geometry, how to establish an adequate mesh grids and how to implement the numerical computation of specific functions used. Hence, after overcome with a solution of the problem stated at introduction, an extension of this present work will be to couple radiation with convection, where the participating medium will be a fluid flow (liquid or gas), and introduce the condition of diffusely reflecting boundary surfaces, to make it more practical and relevant to solve engineering problems.

Appendix A. Annex1

Incident radiation following aperture of 𝜑 𝑆 = 𝑚𝑒𝑠(𝑀 𝑖𝑗 𝑆 𝑊 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ , 𝑀 𝑖𝑗 𝑆 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ ) delimited by boundary surface 𝛤 𝑆 1 :

𝐺 𝑆 1 (7) (𝑖,𝑗) and 𝜑 𝑙 12 = 𝛽 𝑙 × 𝜑 ̅ 𝑁 22 (𝑖,𝑗) .

Appendix A. Annex2

Radiative flux following apperture 𝜑 𝑆 = 𝑚𝑒𝑠(𝑀 𝑖𝑗 𝑆 𝑊 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ , 𝑀 𝑖𝑗 𝑆 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ ) delimited by boundary surface 𝛤 𝑆 1 is given by: 𝑞 ⃗ 𝑟,𝑆 1 (7) = 𝛼𝑇 𝑆

Fig. 1 .

 1 Fig.1. Geometry of the semi-transparent medium in (𝑥, 𝑦) plane for (a) and (b) in (𝑧, 𝑦) plane.

Fig. 2 .

 2 Fig.2. Sub-surfaces of incoming ray: case of 𝐻 = 3ℎ for (a) and 𝐻 𝑥,𝑦 < 3ℎ for (b).

Fig. 3 .

 3 Fig.3. Sub-surfaces of incoming ray: 𝐻 ∈ ]3ℎ, (2 + √5)ℎ[ for (a) and 𝐻 > (2 + √5)ℎ for (b).

Fig. 4 .

 4 Fig.4. Sub-surfaces of incoming ray: case of 𝐻 = (2 + √5)ℎ for (a), and ray propagation within sub-surface 7 for (b).

) 2 . 3 .Fig. 5 .

 235 Fig.5. Meshing of the semi-transparent medium

4 𝜋 4 𝜋 4 𝜋 4 𝜋𝐻 2 [ 2 ,

 444422 and 𝑦 ̅ 𝑗 ∈ ]0, 𝐻[, 𝑞 ⃗ 𝑟 (𝑖, 𝑗) = 𝜎 𝐵 𝑇 𝑊 1 𝑒 ⃗ 𝑦 ; when, 𝑥̅ 𝑖 = 𝐻 and 𝑦 ̅ 𝑗 ∈ ]0, 𝐻[, therefore 𝑞 ⃗ 𝑟 (𝑖, 𝑗) = 𝜎 𝐵 𝑇 𝐸 1 𝑒 ⃗ 𝑦 ; if 𝑥̅ 𝑖 ∈ ]0, 𝐻[ and 𝑦 ̅ 𝑗 = 0, then 𝑞 ⃗ 𝑟 (𝑖, 𝑗) = 𝜎 𝐵 𝑇 𝑆 1 𝑒 ⃗ 𝑥 ; for 𝑥̅ 𝑖 ∈ ]0, 𝐻[ and 𝑦 ̅ 𝑗 = 𝐻, 𝑞 ⃗ 𝑟 (𝑖, 𝑗) = 𝜎 𝐵 𝑇 𝑁 1 𝑒 ⃗ 𝑥 . At boundary surfaces of inner square cavity, the similar expressions are also: 𝑥̅ 𝑖 = and 𝑦 ̅ 𝑗 = 𝐻+ℎ the radiative flux is 𝑞 ⃗ 𝑟 (𝑖, 𝑗) = 𝜎 𝐵 𝑇 𝑁 2 4 𝜋 𝑒 ⃗ 𝑥 .

2 +

 2 √5[ ∪ [2 + √5, +∞[, where 𝐻 = 1𝑚. Bickley-Naylor and Altaç

3. 1 .

 1 Numerical verification 3.1.1. Simple geometry Numerical result of non-dimensional temperature field denoted as 𝑇 * = ( 𝑇(𝑥̅ 𝑖 ,𝑦 ̅ 𝑗 )

Fig. 6 .

 6 Fig.6. Non-dimensional temperature field 𝑇 * at ℎ = 0 , 𝑥 = 𝐻/2 , 𝑁 𝑥 × 𝑁 𝑦 = 251 × 251 grids

Fig. 7 .

 7 Fig.7. local heat flux field in the participating medium with inner square cavity

3. 3 . 1 .Fig. 7 .Fig. 8 .Fig. 9 .

 31789 Fig.7. Non-dimensional temperature field for 𝑇 𝑆 1 ≠ 0, ℎ = 0: (a) 𝜏 = 0.01, (b) 𝜏 = 1, (c) 𝜏 = 5.

Fig. 10 .

 10 Fig.10. Y-component of non-dimensional radiative flux at 𝑇 𝑆 1 ≠ 0 and 𝑇 𝑁 1 ≠ 0, ℎ = 0: (a) 𝜏 = 0.01, (b) 𝜏 = 1, (c) 𝜏 = 5.

3. 3 . 1 .

 31 Fig.11. Non-dimensional temperature field for 𝑇 𝑆 1 ≠ 0 , 𝐻 = (2 + √5)ℎ: (a) 𝜏 = 0.01, (b) 𝜏 = 1, (c) 𝜏 = 5.

Fig. 12 .

 12 Fig.12. Non-dimensional temperature field for 𝑇 𝑊 1 ≠ 0 , 𝐻 = (2 + √5)ℎ: (a) 𝜏 = 0.01, (b) 𝜏 = 1, (c) 𝜏 = 5.

  modifies the magnitude of the radiative quantity expected. When the value of optical thickness increases, the medium cools very quickly.

Fig. 13 .

 13 Fig.13. Non-dimensional temperature field for 𝑇 𝑆 1 ≠ 0 and 𝑇 𝑁 1 ≠ 0, 𝐻 = (2 + √5)ℎ: (a) 𝜏 = 0.01, (b) 𝜏 = 1, (c) 𝜏 = 5.

Fig. 14 .

 14 Fig.14. Y-component of non-dimensional radiative flux at 𝑇 𝑆 1 ≠ 0, 𝐻 = (2 + √5)ℎ: (a) 𝜏 = 0.01, (b) 𝜏 = 1, (c) 𝜏 = 5.

2 .Fig. 15 .

 215 Fig.15. Non-dimensional temperature field for 𝑇 𝑆 1 ≠ 0 and 𝑇 𝑁 1 ≠ 0, 𝐻 = 3ℎ, : (a) 𝜏 = 0.01, (b) 𝜏 = 1, (c) 𝜏 = 5.

Fig. 16 .

 16 Fig.16. Y-component of non-dimensional radiative flux for 𝑇 𝑆 1 ≠ 0 and 𝑇 𝑁 1 ≠ 0, 𝐻 = 3ℎ : (a) 𝜏 = 0.01, (b) 𝜏 = 1, (c) 𝜏 = 5.

Fig. 17

 17 Fig.17. Non-dimensional temperature field for 𝑇 𝑆 1 ≠ 0 and 𝑇 𝑁 1 ≠ 0, 𝐻 = 8ℎ: (a) 𝜏 = 0.01, (b) 𝜏 = 1, (c) 𝜏 = 5.

  Fig.17. Non-dimensional temperature field for 𝑇 𝑆 1 ≠ 0 and 𝑇 𝑁 1 ≠ 0, 𝐻 = 8ℎ: (a) 𝜏 = 0.01, (b) 𝜏 = 1, (c) 𝜏 = 5.

Fig. 18 .

 18 Fig.18. Non-dimensional temperature field for 𝑇 𝐸 1 ≠ 0 and 𝑇 𝑊 1 ≠ 0, 𝐻 = 8ℎ: (a) 𝜏 = 0.01, (b) 𝜏 = 1, (c) 𝜏 = 5.

4 .

 4 Fig.19. Non-dimensional temperature field for 𝑇 𝑆 2 ≠ 0, 𝑇 𝑁 2 ≠ 0, 𝑇 𝑊 2 ≠ 0, 𝑇 𝐸 2 ≠ 0, 𝐻 = 2ℎ: (a) 𝜏 = 0.01, (b) 𝜏 = 1, (c) 𝜏 = 5.

Fig. 20 .Fig. 22 .

 2022 Fig.20. X-component of Non-dimensional radiative flux field for 𝑇 𝑆 2 ≠ 0, 𝑇 𝑁 2 ≠ 0, 𝑇 𝑊 2 ≠ 0, 𝑇 𝐸 2 ≠ 0, 𝐻 = 2ℎ: (a) 𝜏 = 0.01, (b) 𝜏 = 1, (c) 𝜏 = 5.

Fig. 23 .

 23 Fig.23. Non-dimensional temperature field for 𝑇 𝑆 2 ≠ 0, 𝑇 𝑁 2 ≠ 0, 𝑇 𝑊 2 ≠ 0, 𝑇 𝐸 2 ≠ 0, 𝐻 = 8ℎ: (a) 𝜏 = 0.01, (b) 𝜏 = 1, (c) 𝜏 = 5.

Fig. 24 .

 24 Fig.24. Y-component of non-dimensional radiative flux field for 𝑇 𝑆 2 ≠ 0, 𝑇 𝑁 2 ≠ 0, 𝑇 𝑊 2 ≠ 0, 𝑇 𝐸 2 ≠ 0, 𝐻 = 8ℎ.

  non-dimensional radiative flux in the medium (𝑊𝑚 -2 ) 𝑥-component of the radiative flux (𝑊𝑚 -2 ) 𝑦-component of the radiative flux (𝑊𝑚 -2 ) 𝑇 radiation temperature in the medium(𝐾) 𝑇 * non-dimensional temperature in the medium 𝑢 real number at which 𝐵 𝑖𝑠 𝑛 , 𝐶 𝑖𝑠 𝑛 are evaluated

	𝑁 𝜑 , 𝑁 Radiative heat transfer in a two-dimensional semi-transparent medium number of angular quadrature, and cells with a centered inner square cavity 𝑞 𝑥 𝑞 ⃗ 𝑟 * 𝑞 ⃗ 𝑟 vector radiative flux (𝑊𝑚 -2 )
	𝑞 𝑦
	Greek symbols
	∆𝑥, ∆𝑦	length of the cell following 𝑥 and 𝑦 respectively (𝑚)
	𝛽 𝑙 𝜎 𝐵	quadrature angular abscissa Stefan-Boltzmann constant (5.67 10 -8 𝑊. 𝑚 -2 . 𝐾 -4 )
	𝜑, 𝜃	azimuthal and zenith angle of unit vector Ω ⃗⃗⃗
	Ω ⃗⃗⃗	unit radiation propagation vector
	𝜕𝐷 𝑒𝑥𝑡 , 𝜕𝐷 𝑖𝑛𝑡 boundary domains for enclosure and inner square cavity respectively ∇. 𝑞 ⃗ 𝑟 radiative flux divergence(𝑊𝑚 -3 )
	𝛿	radiative path length from boundary surface to attenuated point in the
		medium
	𝛿 𝑚
	𝑘 𝑎	absorption coefficient (𝑚 -1 )
	(𝑒 ⃗ 𝑥 , 𝑒 ⃗ 𝑦 )	unit vectors of 𝑥, 𝑦 directions
	𝐵 𝑖𝑠 𝑛 , 𝐶 𝑖𝑠 𝑛 𝐺	Altaç modified Bessel functions volumic incident radiation (𝑊𝑚 -3 )
	𝐺 *	non-dimensional incident radiation in the medium
	𝐻	length of external cavity along 𝑥 and 𝑦 direction (𝑚)
	ℎ	length of the inner square body (𝑚)
	(𝑖, 𝑗, 𝑝, 𝑞)	cells numbering
	𝐼	one directional incoming radiation intensity (𝑊𝑚 -2 𝑆𝑟)
	𝐼 0	black body radiation intensity (𝑊𝑚 -2 𝑆𝑟)
	𝐾 𝑖 𝑛	Bickley-Naylor functions
	𝑙, 𝑚	angular and spatial numbering quadrature

  The boundary domain, either for enclosure denoted by 𝜕𝐷 𝑒𝑥𝑡 or inner square cavity denoted by 𝜕𝐷 𝑖𝑛𝑡 are at imposed constant temperature, and boundary surfaces are supposed to be black (𝜀 = 1). 𝑇 𝐸 1 , 𝑇 𝑁 1 ,𝑇 𝑊 1 , 𝑇 𝑆 1 denote a constant temperature for east, north, west and south boundary surfaces of the participating medium, and 𝑇 𝐸 2 ,𝑇 𝑁 2 , 𝑇 𝑊 2 , 𝑇 𝑆 2 for east, north, west and south boundary surfaces of the inner square cavity respectively.

  𝑘 𝑎 the absorption coefficient, and 𝐼 𝑏 the Planck's black body radiation intensity depending of temperature medium. Heat transfer in this kind of medium is subjected to boundary conditions: 𝐼(𝑥, 𝑦, 𝜃, 𝜑) = 𝐼 𝑏 1 [𝑇(𝑥, 𝑦)], for (𝑥, 𝑦) ∈ 𝜕𝐷 𝑒𝑥𝑡 , and 𝐼(𝑥, 𝑦, 𝜃, 𝜑) = 𝐼 𝑏 2 [𝑇(𝑥, 𝑦)], for (𝑥, 𝑦) ∈ 𝜕𝐷 𝑖𝑛𝑡 . In this case, 𝐼 𝑏 1 and 𝐼 𝑏 2 are respectively the Planck's source functions at boundary surfaces of the enclosure and at boundary surfaces of the inner square cavity. Let define by, 𝛿 = 𝛿(𝑥, 𝑦, 𝜃, 𝜑) the curvilinear abscissa of the ray pathway; the solution of Eq. (1) for any point along the participating medium enables to obtain 𝐼(𝑥, 𝑦, 𝜃, 𝜑) by the relation:

			𝑦,𝜃,𝜑) 𝜕𝑥	+ 𝑠𝑖𝑛𝜃𝑠𝑖𝑛𝜑	𝜕𝐼(𝑥,𝑦,𝜃,𝜑) 𝜕𝑦	= 𝑘 𝑎 𝐼 𝑏 [𝑇(𝑥, 𝑦)] -𝐼(𝑥, 𝑦, 𝜃, 𝜑)	(2)
	where, 𝐼(𝑥, 𝑦, 𝜃, 𝜑) is the one directional incoming intensity, 𝑇(𝑥, 𝑦) a temperature field in the
	participating medium, 𝐼(𝑥, 𝑦, 𝜃, 𝜑) = 𝐼 0 (𝑥, 𝑦)𝑒 -𝑘 𝑎. 𝛿(𝑥,𝑦,𝜃,𝜑) +	𝑘 𝑎. 𝜎 𝐵 𝜋 ∫ 0 𝛿(𝑥,𝑦,𝜃,𝜑)	𝑇 4	(𝑥 ′ , 𝑦 ′ )𝑒 -𝑘 𝑎. 𝛿 ′ 𝑑𝛿′	(3)
	with, 𝛿 ′ = 𝛿(𝑥 ′ , 𝑦 ′ , 𝜃, 𝜑), for each (𝑥 ′ , 𝑦 ′ ) ∈ [0, 𝑥] × [0, 𝑦]. 𝜎 𝐵 is the Stephan-Boltzman
	constant and 𝐼 0 (𝑥, 𝑦) represents incoming intensity leaving each boundary surface. Since, the
	symmetry around the respective geometries already illustrated allows to study the problem for
	𝜃 ∈ [0,	𝜋 2	], whereas others quantities following 𝜃 ∈ ]	𝜋 2

  {𝐵 𝑖𝑠 2 (𝑘 𝑎 𝑢 𝑗 , 𝜑 ̅ 𝑁 1 (𝑖,𝑗) ) + 𝐵 𝑖𝑠 2 (𝑘 𝑎 𝑢 𝑗 , 𝜑 ̅ 𝑁 2

							(𝑖,𝑗) )}
	+ 𝛼𝑘 𝑎 (𝜑 ̅ 𝑁 1 (𝑖,𝑗) )𝑢 𝑗 ∑	𝑁 𝜑 𝑙=1 ∑	𝑀 𝑚=1	𝜔 𝑚 𝜔 𝑙 𝑐𝑜𝑠𝜑 𝑙 1	𝑇 4 (𝑝, 𝑞)𝐾 𝑖 1 ( 𝑘 𝑎 𝑢 𝑗 𝛿 𝑚 𝑐𝑜𝑠𝜑 𝑙 1	)
	+ 𝛼𝑘 𝑎. (𝜑 ̅ 𝑁 2 (𝑖,𝑗) )𝑢 𝑗 ∑	𝑁 𝜑 𝑙=1 ∑	𝑀 𝑚=1	𝜔 𝑚 𝜔 𝑙 𝑐𝑜𝑠𝜑 𝑙 2	𝑇 4 (𝑝, 𝑞)𝐾 𝑖 1 (	𝑘 𝑎 𝑢 𝑗 𝛿 𝑚 𝑐𝑜𝑠𝜑 𝑙 2

  {𝐵 𝑖𝑠 2 (𝑘 𝑎 𝑣 𝑖 , 𝜑 ̅ 𝑊 1 (𝑖,𝑗) ) + 𝐵 𝑖𝑠 2 (𝑘 𝑎 𝑣 𝑖 , 𝜑 ̅ 𝑊 2 (𝑖,𝑗) )}

	+𝛼𝑘 𝑎. (𝜑 ̅ 𝑊 1 (𝑖,𝑗) )𝑣 𝑖 ∑	𝑁 𝜑 𝑙=1 ∑	𝑀 𝑚=1	𝜔 𝑚 𝜔 𝑙 𝑐𝑜𝑠𝜑 𝑙 3	𝑇 4 (𝑝, 𝑞)𝐾 𝑖 1 (	𝑘 𝑎 (1-𝛿 𝑚 )𝑣 𝑖 𝑐𝑜𝑠𝜑 𝑙 3	)
	+𝛼𝑘 𝑎. (𝜑 ̅ 𝑊 2 (𝑖,𝑗) )𝑣 𝑖 ∑	𝑁 𝜑 𝑙=1 ∑	𝑀 𝑚=1	𝜔 𝑚 𝜔 𝑙 𝑐𝑜𝑠𝜑 𝑙 4	𝑇 4 (𝑝, 𝑞)𝐾 𝑖 1 (	𝑘 𝑎 (1-𝛿 𝑚 )𝑣 𝑖 𝑐𝑜𝑠𝜑 𝑙 4	)

  𝑚𝑒𝑠(𝑀 𝑖𝑗 𝑠 𝑤 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ , 𝑀 𝑖𝑗 𝑛 𝑤 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ ), 𝑚𝑒𝑠(𝑀 𝑖𝑗 𝑛 𝑤 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ , 𝑀 𝑖𝑗 𝑛 𝑒 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ ), 𝜑 𝐸 = 𝑚𝑒𝑠(𝑀 𝑖𝑗 𝐸 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ , 𝑀 𝑖𝑗 𝑁 𝐸 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ ) = 𝜑 𝐸 1 + 𝜑 𝐸 2 and are

  3 (𝑘 𝑎 𝑢 𝑗 , 𝜑 ̅ 𝑁 2 (𝑖,𝑗) ) 𝐵 𝑖𝑠 3 (𝑘 𝑎 𝑢 𝑗 , 𝜑 ̅ 𝑁 1 (𝑖,𝑗) ) + 𝐵 𝑖𝑠 3 (𝑘 𝑎 𝑢 𝑗 , 𝜑 ̅ 𝑁 2

	+ 𝛼 𝑘 𝑎. (𝜑 ̅ 𝑁 2 (𝑖,𝑗) )𝑢 𝑗 ∑	𝑁 𝜑 𝑙=1 ∑	𝑀 𝑚=1	𝜔 𝑚 𝜔 𝑙 𝑐𝑜𝑠𝜑 𝑙 2	𝑇 4 (𝑝, 𝑞)𝐾 𝑖 2 ( 𝑘 𝑎 𝑢 𝑗 𝛿 𝑚 𝑐𝑜𝑠𝜑 𝑙 2	) ( -𝑠𝑖𝑛𝜑 𝑙 2 𝑐𝑜𝑠𝜑 𝑙 2	)
								(𝑖,𝑗) ) )
	+ 𝛼𝑘 𝑎. (𝜑 ̅ 𝑁 1 (𝑖,𝑗) )𝑢 𝑗 ∑	𝑁 𝜑 𝑙=1 ∑	𝑀 𝑚=1	𝜔 𝑚 𝜔 𝑙 𝑐𝑜𝑠𝜑 𝑙 1	𝑇 4 (𝑝, 𝑞)𝐾 𝑖 2 (	𝑘 𝑎 𝑢 𝑗 𝛿 𝑚 𝑐𝑜𝑠𝜑 𝑙 1	) (	𝑠𝑖𝑛𝜑 𝑙 1 𝑐𝑜𝑠𝜑 𝑙 1	)

  3 (𝑘 𝑎 𝑣 𝑖 , 𝜑 ̅ 𝑊 1 (𝑖,𝑗) )-𝐵 𝑖𝑠 3 (𝑘 𝑎 𝑣 𝑖 , 𝜑 ̅ 𝑊 2 (𝑖,𝑗) ) -𝐶 𝑖𝑠 3 (𝑘 𝑎 𝑣 𝑖 , 𝜑 ̅ 𝑊 1 (𝑖,𝑗) ) + 𝐶 𝑖𝑠 3 (𝑘 𝑎 𝑣 𝑖 , 𝜑 ̅ 𝑊 2

							(𝑖,𝑗) ) )
	+ 𝛼𝑘 𝑎 (𝜑 ̅ 𝑊 1 (𝑖,𝑗) )𝑣 𝑖 ∑	𝑁 𝜑 𝑙=1 ∑	𝑀 𝑚=1	𝜔 𝑚 𝜔 𝑙 𝑐𝑜𝑠𝜑 𝑙 4	𝑇 4 (𝑝, 𝑞)𝐾 𝑖 2 ( 𝑘 𝑎 (1-𝛿 𝑚 )𝑣 𝑖 𝑐𝑜𝑠𝜑 𝑙 4	) ( -𝑐𝑜𝑠𝜑 𝑙 3 -𝑠𝑖𝑛𝜑 𝑙 3	)
	+ 𝛼𝑘 𝑎 (𝜑 ̅ 𝑊 2 (𝑖,𝑗) )𝑣 𝑖 ∑	𝑁 𝜑 𝑙=1 ∑	𝑀 𝑚=1	𝜔 𝑚 𝜔 𝑙 𝑐𝑜𝑠𝜑 𝑙 3	𝑇 4 (𝑝, 𝑞)𝐾 𝑖 2 ( 𝑘 𝑎 (1-𝛿 𝑚 )𝑣 𝑖 𝑐𝑜𝑠𝜑 𝑙 3	) ( -𝑐𝑜𝑠𝜑 𝑙 4 𝑠𝑖𝑛𝜑 𝑙 4	)

  the components of radiative flux vector, 𝑞 𝑟,𝑁 1 𝑥 , 𝑞 𝑟,𝑊 1 𝑥 , 𝑞 𝑟,𝑊 1 𝑦 and 𝑞 𝑟,𝑁 1 𝑦 following 𝑥 and 𝑦 directions are carried respectively by specific functions 𝐶 𝑖𝑠 𝑛 and 𝐵 𝑖𝑠 𝑛 , and 𝐾 𝑖 𝑛 .

  is parabolic decreasing, because the values of temperature variate within the participating medium following parabolic shapes, from hot surface 𝛤 𝑆 1 to cold surface 𝛤 𝑁 1 . Verification on calculations is more precise at 𝑥 =

	𝐻 2	and 𝑦 =	𝐻 2

Table 1

 1 Table1and 5 present it, on lines Non-dimensional temperature field 𝑇 * within the semi-transparent medium, for 𝐻 = 3ℎ, 𝑘 𝑎 = 0.25 𝑚 -1 , 𝑁 𝑥 × 𝑁 𝑦 = 251 × 251,𝑇 𝑆 1 ≠ 0, the rest of boundary surfaces are cold.

	𝑥 𝐻	= 0.5 and	𝑥 𝐻	=

Table 2

 2 Non-dimensional incident radiation 𝐺 * inside the semi-transparent medium, for 𝐻 = 3ℎ, 𝑘 𝑎 = 0.25 𝑚 -1 , 𝑁 𝑥 × 𝑁 𝑦 = 251 × 251,𝑇 𝑆 1 ≠ 0, the rest of boundary surfaces are cold. will let more rays pass from the heat surface 𝛤 𝑆 1 to 𝛤 𝑁 1 compared as the one of 𝐻 = 3ℎ. It is important to mention that, prediction is well done when the size of the inner square cavity variates. On the Tables 1 and Tables 5, it is noticeable that, temperature decreases in the semi-transparent medium, when one moves away from the south to the north boundary surface. At the same time, the choice of quadrature used in computation affects also the desired values, but with a slight difference; it is more visible on Table4. A similar analysis is done for the non-dimensional incident radiation given by Table2 and 3. The respective values of non-dimensional temperature field 𝑇 * and incident radiation 𝐺 * =

	𝐺 *

Table 3

 3 Non-dimensional incident radiation 𝐺 * inside the semi-transparent medium, for 𝐻 = (2 + √5)ℎ, 𝑘 𝑎 = 0.25 𝑚 -1 , 𝑁 𝑥 × 𝑁 𝑦 = 301 × 301, 𝑇 𝑆 1 ≠ 0, the rest of boundary surfaces are cold.

	𝐺 *

Table 4

 4 Influence of cell numbers on temperature field 𝑇 * at 𝑥 = 𝐻/2 , 𝑘 𝑎 = 0.25 𝑚 -1 .

	𝐻 = (2 + √5)ℎ

Table 5

 5 Non-dimensional temperature field 𝑇 * inside the semi-transparent medium, for 𝐻 = (2 + √5)ℎ, 𝑘 𝑎 = 0.25 𝑚 -1 , 𝑁 𝑥 × 𝑁 𝑦 = 301 × 301, 𝑇 𝑆 1 ≠ 0, the rest of boundary surfaces are cold. sub-sections 2.3.1 and 2.3.4. When the number of cells increases, it allows obtaining accurate results; nevertheless, if the number of cells choiced becomes extremely important, the specific functions 𝐾 𝑖𝑛 , 𝐵 𝑖𝑠 𝑛 , 𝐶 𝑖𝑠 𝑛 , can make diverge the results.

	𝑻 *

  𝜑 𝑙 5 = 𝛽 𝑙 × 𝜑 ̅ 𝑆 1(𝑖,𝑗) and 𝜑 𝑙 6 = 𝛽 𝑙 × 𝜑 ̅ 𝑆 2(𝑖,𝑗) .Following aperture of 𝜑 𝐸 = 𝑚𝑒𝑠(𝑀 𝑖𝑗 𝐸 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ , 𝑀 𝑖𝑗 𝑁 𝐸 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ ) and delimited by boundary surface 𝛤 𝐸 1 is : 𝜑 𝑙 9 = 𝛽 𝑙 × 𝜑 ̅ 𝑊 21 (𝑖,𝑗) and 𝜑 𝑙 10 = 𝛽 𝑙 × 𝜑 ̅ 𝑊 22 (𝑖,𝑗) . Following aperture of 𝜑 𝑁 2 = 𝑚𝑒𝑠(𝑀 𝑖𝑗 𝑛 𝑤 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ , 𝑀 𝑖𝑗 𝑛 𝑒 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ ) and delimited by boundary surface 𝛤 𝑁 2 is: 𝜑 𝑙 11 = 𝛽 𝑙 × 𝜑 ̅ 𝑁 21

	+𝛼𝑘 𝑎 ( 𝜑 ̅ 𝑁 21 (𝑖,𝑗) ) ( ∆𝑦 𝑣 𝑗	-	(𝐻+ℎ) 2∆𝑦 ) ∆𝑦 ∑	𝑁 𝜑 𝑙=1 ∑	𝑀 𝑚=1	𝜔 𝑚 𝜔 𝑙 𝑐𝑜𝑠𝜑 𝑙 11	𝑇 4 (𝑝, 𝑞)𝐾 𝑖 1 (	𝑘 𝑎 ( ∆𝑦 𝑣 𝑗 𝑐𝑜𝑠𝜑 𝑙 11 -(𝐻+ℎ) 2∆𝑦 )∆𝑦𝛿 𝑚	)
	-𝛼𝑘 𝑎 (𝜑 ̅ 𝑁 22 (𝑖,𝑗) ) ( ∆𝑦 𝑣 𝑗	-	(𝐻+ℎ) 2∆𝑦 ) ∆𝑦 ∑	𝑁 𝜑 𝑙=1 ∑	𝑀 𝑚=1	𝜔 𝑚 𝜔 𝑙 𝑐𝑜𝑠𝜑 𝑙 12	𝑇 4 (𝑝, 𝑞)𝐾 𝑖 1 (	𝑘 𝑎 ( ∆𝑦 𝑣 𝑗 𝑐𝑜𝑠𝜑 𝑙 12 -(𝐻+ℎ) 2∆𝑦 )∆𝑦𝛿 𝑚	)
	(𝑖, 𝑗) = 𝛼𝑇 𝑆 1 4 {𝐵 𝑖𝑠 2 (𝑘 𝑎 𝑣 𝑗 , 𝜑 ̅ 𝑆 1 (𝑖,𝑗) ) + 𝐵 𝑖𝑠 2 (𝑘 𝑎 𝑣 𝑗 , 𝜑 ̅ 𝑆 2 (𝑖,𝑗) )} + 𝛼𝑘 𝑎. (𝜑 ̅ 𝑆 1 (𝑖,𝑗) )𝑣 𝑗 ∑ 𝑁 𝜑 𝑙=1 ∑ 𝜔 𝑚 𝜔 𝑙 𝑐𝑜𝑠𝜑 𝑙 5 𝑀 𝑚=1 𝑇 4 (𝑝, 𝑞)𝐾 𝑖 1 ( 𝑘 𝑎 (1-𝛿 𝑚 )𝑣 𝑗 𝑐𝑜𝑠𝜑 𝑙 5 with, 𝜑 ̅ 𝑁 21 (𝑖,𝑗) = 𝑡𝑎𝑛 -1 { (𝐻+ℎ) 2∆𝑥 -𝑣 𝑖 ∆𝑥 𝑣 𝑗 ∆𝑦 -(𝐻+ℎ) 2∆𝑦 }, 𝜑 ̅ 𝑁 22 (𝑖,𝑗) = 𝑡𝑎𝑛 -1 { (𝐻-ℎ) 𝑣 𝑖 2∆𝑥 -∆𝑥 𝑣 𝑗 ∆𝑦 2∆𝑦 -(𝐻+ℎ) },	)
	+ 𝛼𝑘 𝑎. (𝜑 ̅ 𝑆 2 (𝑖,𝑗) )𝑣 𝑗 ∑	𝑁 𝜑 𝑙=1 ∑	𝑀 𝑚=1	𝜔 𝑚 𝜔 𝑙 𝑐𝑜𝑠𝜑 𝑙 6	𝑇 4 (𝑝, 𝑞)𝐾 𝑖 1 ( 𝑘 𝑎 (1-𝛿 𝑚 )𝑣 𝑗 𝑐𝑜𝑠𝜑 𝑙 6	)
	with, 𝜑 ̅ 𝑆 1 (𝑖,𝑗) = 𝑡𝑎𝑛 -1 { 𝑣 𝑗 } , 𝜑 ̅ 𝑆 2 𝑣 𝑖 (𝑖,𝑗) = 𝑡𝑎𝑛 -1 { } , 𝐺 𝐸 1 (𝐻-ℎ) 2∆𝑥 𝑣 𝑖 -∆𝑥 𝑣 𝑗 ∆𝑦 -(𝐻-ℎ) 2∆𝑦 (7) (𝑖, 𝑗) = 𝛼𝑇 𝐸 1 4 {𝐵 𝑖𝑠 2 (𝑘 𝑎 𝑢 𝑖 , 𝜑 ̅ 𝐸 1 (𝑖,𝑗) ) + 𝐵 𝑖𝑠 2 (𝑘 𝑎 𝑢 𝑖 , 𝜑 ̅ 𝐸 2 (𝑖,𝑗) )}
	+𝛼𝑘 𝑎. (𝜑 ̅ 𝐸 1 (𝑖,𝑗) )𝑢 𝑖 ∑	𝑁 𝜑 𝑙=1 ∑	𝑀 𝑚=1	𝜔 𝑚 𝜔 𝑙 𝑐𝑜𝑠𝜑 𝑙 7	𝑇 4 (𝑝, 𝑞)𝐾 𝑖 1 ( 𝑘 𝑎 𝑢 𝑖 𝛿 𝑚 𝑐𝑜𝑠𝜑 𝑙 7	)
	+𝛼𝑘 𝑎. (𝜑 ̅ 𝐸 2 (𝑖,𝑗) )𝑢 𝑖 ∑	𝑁 𝜑 𝑙=1 ∑	𝑀 𝑚=1	𝜔 𝑚 𝜔 𝑙 𝑐𝑜𝑠𝜑 𝑙 8	𝑇 4 (𝑝, 𝑞)𝐾 𝑖 1 ( 𝑘 𝑎 𝑢 𝑖 𝛿 𝑚 𝑐𝑜𝑠𝜑 𝑙 8	)
	with, 𝜑 ̅ 𝐸 1 (𝑖,𝑗) = 𝑡𝑎𝑛 -1 { 𝑢 𝑗 𝑢 𝑖 } , 𝜑 ̅ 𝐸 2 (𝑖,𝑗) = 𝑡𝑎𝑛 -1 {	(𝐻+ℎ) 2∆𝑦 -(𝐻+ℎ) 2∆𝑥 -∆𝑦 𝑣 𝑗 ∆𝑥 𝑣 𝑖	}, 𝜑 𝑙 7 = 𝛽 𝑙 × 𝜑 ̅ 𝐸 1 (𝑖,𝑗) and 𝜑 𝑙 8 = 𝛽 𝑙 × 𝜑 ̅ 𝐸 2 (𝑖,𝑗) .
	Following aperture of 𝜑 𝑊 2 = 𝑚𝑒𝑠(𝑀 𝑖𝑗 𝑠 𝑤 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ , 𝑀 𝑖𝑗 𝑛 𝑤 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ ) and delimited by boundary surface 𝛤 𝑊 2 is :
	𝐺 𝑊 2 (7) (𝑖, 𝑗) = 𝛼𝑇 𝑊 2 4 {𝐵 𝑖𝑠 2 (𝑘 𝑎 ( (𝐻-ℎ) 2∆𝑥 -	𝑣 𝑖 ∆𝑥 ) ∆𝑥, 𝜑 ̅ 𝑊 21 (𝑖,𝑗) ) -𝐵 𝑖𝑠 2 (𝑘 𝑎 ( (𝐻-ℎ) 2∆𝑥 -	𝑣 𝑖 ∆𝑥 ) ∆𝑥, 𝜑 ̅ 𝑊 22 (𝑖,𝑗) )}
	+ 𝛼𝑘 𝑎. (𝜑 ̅ 𝑊 21 (𝑖,𝑗) ) ( (𝐻-ℎ) 2∆𝑥	-	𝑣 𝑖 ∆𝑥	) ∆𝑥 ∑	𝑁 𝜑 𝑙=1 ∑	𝑀 𝑚=1	𝜔 𝑚 𝜔 𝑙 𝑐𝑜𝑠𝜑 𝑙 9	𝑇 4 (𝑝, 𝑞)𝐾 𝑖 1 (	𝑘 𝑎 ( (𝐻-ℎ) 2∆𝑥 -∆𝑥 𝑣 𝑖 𝑐𝑜𝑠𝜑 𝑙 9 )∆𝑥𝛿 𝑚	)
	-𝛼𝑘 𝑎. (𝜑 ̅ 𝑊 22 (𝑖,𝑗) ) ( (𝐻-ℎ) 2∆𝑥	-	𝑣 𝑖 ∆𝑥	) ∆𝑥 ∑	𝑁 𝜑 𝑙=1 ∑	𝑀 𝑚=1	𝜔 𝑚 𝜔 𝑙 𝑐𝑜𝑠𝜑 𝑙 10	𝑇 4 (𝑝, 𝑞)𝐾 𝑖 1 (	𝑘 𝑎 ( (𝐻-ℎ) 2∆𝑥 𝑐𝑜𝑠𝜑 𝑙 10 -𝑣 𝑖 ∆𝑥 )∆𝑥𝛿 𝑚	)
	with, 𝜑 ̅ 𝑊 21 (𝑖,𝑗) = 𝑡𝑎𝑛 -1 { }, 𝐺 𝑁 2 𝑣 𝑗 ∆𝑦 -(𝐻-ℎ) 2∆𝑦 (𝐻-ℎ) 2∆𝑥 -𝑣 𝑖 ∆𝑥 }, 𝜑 ̅ 𝑊 22 (𝑖,𝑗) = 𝑡𝑎𝑛 -1 { 𝑣 𝑗 ∆𝑦 (𝐻+ℎ) -2∆𝑦 (𝐻-ℎ) 2∆𝑥 -𝑣 𝑖 ∆𝑥 (7) (𝑖, 𝑗) = 𝛼𝑇 𝑁 2 4 {𝐵 𝑖𝑠 2 (𝑘 𝑎 ( 𝑣 𝑗 ∆𝑦 -(𝐻+ℎ) 2∆𝑦 ) ∆𝑦, 𝜑 ̅ 𝑁 21 (𝑖,𝑗) ) -𝐵 𝑖𝑠 2 (𝑘 𝑎 ( ∆𝑦 -𝑣 𝑗	(𝐻+ℎ) 2∆𝑦 ) ∆𝑦, 𝜑 ̅ 𝑁 22 (𝑖,𝑗) )}

  1 4 ( -𝐶 𝑖𝑠 3 (𝑘 𝑎 𝑣 𝑗 , 𝜑 ̅ 𝑆 1 (𝑖,𝑗) ) + 𝐶 𝑖𝑠 3 (𝑘 𝑎 𝑣 𝑗 , 𝜑 ̅ 𝑆 2 (𝑖,𝑗) ) -𝐵 𝑖𝑠 3 (𝑘 𝑎 𝑣 𝑗 , 𝜑 ̅ 𝑆 1 (𝑖,𝑗) ) -𝐵 𝑖𝑠 3 (𝑘 𝑎 𝑣 𝑗 , 𝜑 ̅ 𝑆 2 (𝑖,𝑗) ) ) + 𝛼𝑘 𝑎 (𝜑 ̅ 𝑆 1 (𝑖,𝑗) )𝑣 𝑗 ∑ Radiative flux following apperture 𝜑 𝐸 = 𝑚𝑒𝑠(𝑀 𝑖𝑗 𝐸 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ , 𝑀 𝑖𝑗 𝑁 𝐸 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ ) delimited by boundary surface 𝛤 𝐸 1 is given by: 𝑞 ⃗ 𝑟,𝐸 1 (7) = 𝛼𝑇 𝐸 1 4 ( 𝐵 𝑖𝑠 3 (𝑘 𝑎 𝑢 𝑖 , 𝜑 ̅ 𝐸 1 (𝑖,𝑗) ) + 𝐵 𝑖𝑠 3 (𝑘 𝑎 𝑢 𝑖 , 𝜑 ̅ 𝐸 2 (𝑖,𝑗) ) 𝐶 𝑖𝑠 3 (𝑘 𝑎 𝑢 𝑖 , 𝜑 ̅ 𝐸 1 (𝑖,𝑗) ) -𝐶 𝑖𝑠 3 (𝑘 𝑎 𝑢 𝑖 , 𝜑 ̅ 𝐸 2 Radiative flux following apperture 𝜑 𝑊 2 = 𝑚𝑒𝑠(𝑀 𝑖𝑗 𝑠 𝑤 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ , 𝑀 𝑖𝑗 𝑛 𝑤 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ ) delimited by boundary surface 𝛤 𝑊 2 is shown as: Radiative flux following apperture 𝜑 𝑁 2 = 𝑚𝑒𝑠(𝑀 𝑖𝑗 𝑛 𝑤 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ , 𝑀 𝑖𝑗 𝑛 𝑒 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ ) and delimited by boundary surface 𝛤 𝑁 2 is set as :

	𝑞 ⃗ 𝑟,𝑊 2 (7) = 𝛼𝑇 𝑊 2 4 (	-𝐵 𝑖𝑠 3 (𝑘 𝑎 ( 𝐶 𝑖𝑠 3 (𝑘 𝑎 ( (𝐻-ℎ) (𝐻-ℎ) 2∆𝑥 2∆𝑥 --𝑣 𝑖 𝑣 𝑖 ∆𝑥 ∆𝑥 ) ∆𝑥, 𝜑 ̅ 𝑊 21 ) ∆𝑥, 𝜑 ̅ 𝑊 21 (𝑖,𝑗) ) + 𝐵 𝑖𝑠 3 (𝑘 𝑎 ( (𝐻-ℎ) 2∆𝑥 (𝑖,𝑗) ) -𝐶 𝑖𝑠 3 (𝑘 𝑎 ( (𝐻-ℎ) 2∆𝑥 --∆𝑥 𝑣 𝑖 ∆𝑥 ) ∆𝑥, 𝜑 ̅ 𝑊 22 (𝑖,𝑗) ) ) ∆𝑥, 𝜑 ̅ 𝑊 22 (𝑖,𝑗) ) 𝑣 𝑖 )
	-𝛼𝑘 𝑎. (𝜑 ̅ 𝑊 21 (𝑖,𝑗) ) ( (𝐻-ℎ) 2∆𝑥	-	𝑣 𝑖 ∆𝑥	) ∆𝑥 ∑	𝑁 𝜑 𝑙=1 ∑	𝑀 𝑚=1	𝜔 𝑚 𝜔 𝑙 𝑐𝑜𝑠𝜑 𝑙 10	𝑇 4 (𝑝, 𝑞)𝐾 𝑖 2 (	𝑘 𝑎 ( (𝐻-ℎ) 2∆𝑥 𝑐𝑜𝑠𝜑 𝑙 10 -𝑣 𝑖 ∆𝑥 )∆𝑥𝛿 𝑚	) ( -𝑠𝑖𝑛𝜑 𝑙 10 𝑐𝑜𝑠𝜑 𝑙 10	)
	+ 𝛼𝑘 𝑎. (𝜑 ̅ 𝑊 22 (𝑖,𝑗) ) ( (𝐻-ℎ) 2∆𝑥	-	𝑣 𝑖 ∆𝑥	) ∆𝑥 ∑	𝑁 𝜑 𝑙=1 ∑	𝑀 𝑚=1	𝜔 𝑚 𝜔 𝑙 𝑐𝑜𝑠𝜑 𝑙 9	𝑇 4 (𝑝, 𝑞)𝐾 𝑖 2 (	𝑘 𝑎 ( (𝐻-ℎ) 2∆𝑥 -∆𝑥 𝑣 𝑖 𝑐𝑜𝑠𝜑 𝑙 9 )∆𝑥𝛿 𝑚	) ( -𝑠𝑖𝑛𝜑 𝑙 9 𝑐𝑜𝑠𝜑 𝑙 9	)
	𝑞 ⃗ 𝑟,𝑁 2 (7) = 𝛼𝑇 𝑁 2 4 ( -𝐵 𝑖𝑠 3 (𝑘 𝑎 ( 𝐶 𝑖𝑠 3 (𝑘 𝑎 ( 𝑣 𝑗 ∆𝑦 𝑣 𝑗 -∆𝑦 -(𝐻+ℎ) 2∆𝑦 ) ∆𝑦, 𝜑 ̅ 𝑁 21 (𝑖,𝑗) ) -𝐶 𝑖𝑠 3 (𝑘 𝑎 ( ∆𝑦 𝑣 𝑗 (𝐻+ℎ) 2∆𝑦 ) ∆𝑦, 𝜑 ̅ 𝑁 21 (𝑖,𝑗) ) + 𝐵 𝑖𝑠 3 (𝑘 𝑎 ( 𝑣 𝑗 -∆𝑦 -(𝐻+ℎ) 2∆𝑦 ) ∆𝑦, 𝜑 ̅ 𝑁 22 (𝑖,𝑗) ) 2∆𝑦 ) ∆𝑦, 𝜑 ̅ 𝑁 22 (𝑖,𝑗) ) (𝐻+ℎ) )
	𝑁 𝜑 𝑙=1 ∑ 𝑁 𝜑 𝑙=1 ∑ 𝑀 𝑀 𝑚=1 𝑐𝑜𝑠𝜑 𝑙 11 𝜔 𝑚 𝜔 𝑙 𝑐𝑜𝑠𝜑 𝑙 5 𝜔 𝑚 𝜔 𝑙 𝑚=1 (𝑖,𝑗) )𝑣 𝑗 ∑ (𝐻+ℎ) 2∆𝑦 ) ∆𝑦 ∑ + 𝛼𝑘 𝑎 (𝜑 ̅ 𝑆 2 (𝑖,𝑗) ) ( +𝛼𝑘 𝑎. (𝜑 ̅ 𝑁 21 𝑣 𝑗 ∆𝑦 -𝑁 𝜑 𝑙=1 ∑ 𝜔 𝑚 𝜔 𝑙 𝑐𝑜𝑠𝜑 𝑙 6 𝑀 𝑚=1 -𝛼𝑘 𝑎. (𝜑 ̅ 𝑁 22 (𝑖,𝑗) ) ( 𝑣 𝑗 ∆𝑦 -(𝐻+ℎ) 2∆𝑦 ) ∆𝑦 ∑ 𝑁 𝜑 𝑙=1 ∑ 𝜔 𝑚 𝜔 𝑙 𝑐𝑜𝑠𝜑 𝑙 12 𝑇 4 (𝑝, 𝑞)𝐾 𝑖 2 ( 𝑇 4 (𝑝, 𝑞)𝐾 𝑖 2 ( 𝑘 𝑎 (1-𝛿 𝑚 )𝑣 𝑗 𝑐𝑜𝑠𝜑 𝑙 5 𝑘 𝑎 ( 𝑣 𝑗 ∆𝑦 -(𝐻+ℎ) ) ( 2∆𝑦 )∆𝑦𝛿 𝑚 -𝑠𝑖𝑛𝜑 𝑙 5 𝑐𝑜𝑠𝜑 𝑙 11 ) ( -𝑠𝑖𝑛𝜑 𝑙 11 𝑐𝑜𝑠𝜑 𝑙 11 ) -𝑐𝑜𝑠𝜑 𝑙 5 𝑇 4 (𝑝, 𝑞)𝐾 𝑖 2 ( 𝑘 𝑎 (1-𝛿 𝑚 )𝑣 𝑗 𝑐𝑜𝑠𝜑 𝑙 6 ) ( 𝑠𝑖𝑛𝜑 𝑙 6 -𝑐𝑜𝑠𝜑 𝑙 6 ) 𝑀 𝑚=1 𝑇 4 (𝑝, 𝑞)𝐾 𝑖 2 ( 𝑘 𝑎 ( 𝑣 𝑗 ∆𝑦 -(𝐻+ℎ) 2∆𝑦 )∆𝑦𝛿 𝑚 𝑐𝑜𝑠𝜑 𝑙 12 ) ( -𝑠𝑖𝑛𝜑 𝑙 12 ) 𝑐𝑜𝑠𝜑 𝑙 12	)
												(𝑖,𝑗) )	)
	+ 𝛼𝑘 𝑎. (𝜑 ̅ 𝐸 1 (𝑖,𝑗) )𝑢 𝑖 ∑	𝑁 𝜑 𝑙=1 ∑	𝑀 𝑚=1	𝜔 𝑚 𝜔 𝑙 𝑐𝑜𝑠𝜑 𝑙 7	𝑇 4 (𝑝, 𝑞)𝐾 𝑖 2 (	𝑘 𝑎 𝑢 𝑖 𝛿 𝑚 𝑐𝑜𝑠𝜑 𝑙 7	) (	𝑐𝑜𝑠𝜑 𝑙 7 𝑠𝑖𝑛𝜑 𝑙 7	)
	+ 𝛼𝑘 𝑎. (𝜑 ̅ 𝐸 2 (𝑖,𝑗) )𝑢 𝑖 ∑	𝑁 𝜑 𝑙=1 ∑	𝑀 𝑚=1	𝜔 𝑚 𝜔 𝑙 𝑐𝑜𝑠𝜑 𝑙 8	𝑇 4 (𝑝, 𝑞)𝐾 𝑖 2 (	𝑘 𝑎 𝑢 𝑖 𝛿 𝑚 𝑐𝑜𝑠𝜑 𝑙 8	) ( -𝑠𝑖𝑛𝜑 𝑙 8 𝑐𝑜𝑠𝜑 𝑙 8	)

Acknowledgment

Institut Pprime, University of Poitiers, CNRS, and ENSMA UPR 3346 in France full supported this work. The authors would like to extend their thanks to Pr. Karl Joulain, the director of Institut Pprime, who provided us an adequate logistic and digital equipment.