1. Resolution of the Radiatif Integro-Differential Equation
The fractional step method has been proposed by Marchuk G. I.  for resolution on the “Neutron Transport Equation” and developed by Yanenko N. N.  .
This work adapted this method to solve the radiative thermal equation. This differential equation governs the variation of the luminance in a gray medium with non-reflecting faces. The volume method is used in space  and the fractional step method with time.
The equation is written in the following form:
t is the time, the optical thickness, L the luminance;
With the boundary conditions:
at and (2)
at and (3)
1.2. The Equation’s Models
We have the following general equation, according to those boundary and initial conditions:
for and (5)
for and (6)
and the functions:
・ and depend however of the three variables ;
・ and depend only of the variable , with and as continuous and limited.
Some new variables are also introduced here:
and , for
The first equation then becomes:
With the following boundary conditions:
and the initial condition is, for t = 0:
C = 1 is assumed.
Another change of variables is introduced as well:
The Equations (7.3), (7.4) and (7.2) become:
with the boundary conditions (7.4) and (7.5) being:
The domain of study is .
Consider the space with the following scalar product:
1.3. Evolutionary Writing of the Equation (4)
The Equation (4) with the conditions (5-7), replaced by the Equations (7.6)-(7.10) is presented in the following form:
T: is a period of time.
F: belong to and W0, W(t) belong to ψ and ψ0 respectively.
・ ψ: subspace of functions of L2(D), verifying (AW, W) < ∞.
・ ψ0: subspace of ψ which is formed from the functions which components are continuous in D, and their derivative functions with τ are continuous as well.
where is constant > 0.
2. The Numerical Method
2.1. Approximation According to the Component of Space
The length H is discretized in n nodes τi and the node τi+1/2 is placed between two nodes consecutives τi and τi+1.
The Equation (7.6) is integrated on the ranges , and , and the Equation (7.7) is integrated on the ranges , we obtain the following equations:
* For i = 0:
* For 1 ≤ i ≤ n − 1:
* and for i = n
Introduction of the following approximations:
The components ui and vi+1/2 so check to:
* for i = 0:
Considering to and
* for i=1, n-1:
* and for i = n:
The condition is also taken into account:
Let M (0, 2n) the Hilbert space of vector functions
provided the following scalar product:
Introduce the vectors:
After estimations, the problem is reduced to the following evolutionary system:
where A is an integro-differential operator which is defined according to two operators L and S such as:
We note that A is coercive, ensuring the stability of the process.
2.2. Approximation in the Time
This approximation will be based on the fractional step method, which means making the evolution between two principal steps of time. The process will be described after the decomposition, known as “splitting”, of operator A:
The operators are defined positive on has the following algorithm:
* on : (24)
* on : (26)
* on : (27)
The Fj element is an approximation of F, defined by:
2.3. Final Diagram
With the presented approximations of the three parameters z, μ and t, Equations (15) and (17) can be replaced; for instance, for Equation (15):
Let us consider the following equation:
2.4. Approximation on the Level of μ
The interval [0, 1] is subdivided into the m-1 sub-intervals, which lengths are dμl with .
We substitute by with .
2.5. Definitive Algorithm
The numerical algorithm to solve the integro-differential equation is therefore completed.
After the resolution, we obtain the vector Φl formed by the computed values with principal and the values of ul calculated with the inserted nodes. We will be able to determine and by using the following relation:
The solution of is given by:
2.6. Numerical Results and Interpretation
The results are initially presented with an exponential profile of the temperature. Then, the model’s sensitivities are tested, taking into consideration several physical parameters, such as the albedo or the coefficient of extinction.
The variable C is the celerity of light in the medium. To simplify the study, it is taken equal to 1, because it does not influence the calculation in steady state, which interests us here.
But in the case of an unsteady regime, we return to the evolutionary writing of the problem, which is written with the presence of constant C in the following form:
Which is equivalent to:
It is assumed that C is variable here (and not only equal to 1), a change of variables has to be introduced: and
2) Adaptation of the model
Considering the nature of the boundary conditions related to the modelling of the radiative transfer, the solution of the general equation does not correspond exactly to the one we are expecting. A change of variables has to be introduced.
Where and are the solution of the general model carried out with the change affecting its second member:
3) Interpretation of the results
For an exponential profile of the temperature, as seen in Figure 1, the variation of radiative flux is done in the same direction as that of the temperature. This variation is shown in Figure 2. This agreement in the variations in the temperature
Figure 1. The variation of the temperature profile with Optical thickness T (103 K); example: for 1.2, T = 1200 K.
Figure 2. Evolution of the radiative flux with the Optical thickness.
and radiative flux is justified by approximation of Rosseland  . This approximation is expressed by:
where the radiative transfer is characterized by a radiative conductivity:
k (m−1) is the extinction coefficient per unit volume of the packed bed.
In this study, we were interested in the stationary solution. Figure 3 shows that the solved equation has a good stationary solution.
To see if our code reflects the physical aspect of the radiation phenomena, we determined the effect of different physical parameters of the radiation by using the same profile of temperature.
Influence of albedo:
We present in Figure 4 the effect that the albedo had on the radiative flux for
Figure 3. The variation of radiative flux in accordance with time.
Figure 4. The effect of the albedo on radiative flux with the Optical thickness (1, for ω = 0; 2, for ω = 0.3; 3, for ω = 0.5; 4, for ω = 0.7).
a coefficient of extinction equal to 1000.
It is noticed that when the albedo decreases, the radiative flux increases. This result corresponds to the data of the literature   .
The divergence of the flux vector of radiation is more important when the albedo is low. Physically, for the semi transparent medium not diffusing (no albedo), the effect of the radiative transfer on material is maximum. In contrast, for a medium entirely diffusing (albedo equal to one), only the conductive exchanges influence the establishment of the profile of temperature in the medium.
Consequently, the diffusion tends to decrease the importance of the radiative heat exchanges in the medium.
Influence coefficient of extinction:
We show in Figure 5 the radiative profiles for various coefficients of extinction K (noted beta), for an albedo of 0.5. It is observed that when the coefficient of extinction increases, the radiative flux decreases. This variation is foreseeable, considering the expression of the albedo, according to the coefficient of extinction as well as the preceding effect of the albedo.
Note: these results are obtained with the following parameters:
・ The step of space is equal to 1.
・ The step of time is equal to 10−4.
・ The error is of about to 10−6.
・ l is equal to 10−1.
Figure 5. The effect of the coefficient of extinction (for ω = 0.5) on the evolution of radiative flux profile, (1 for k = 1000 m−1, 2 for k = 1500 m−1, 3 for k = 2000 m−1).
A new numerical method of solving the radiative transfer equation is presented here. It is based on the fractional step method proposed by G. I. Marchuk  and modified by A. Belghit. It seems to be a relatively more simple method compared to the conventional resolution methods (see Case’s normal-mode expansion technique in  , for example). It will be adapted to the resolution of heat and mass transfers in a reactive porous medium.
This work was done with the help of the Physics Department of the Faculty of Sciences (Marrakesh), of Mr. Toufik Cherfi (Civil engineering department, University of La Rochelle) and of Ms. Nour Belghit (of Hyphen Translation).
C Speed of light in the medium, m/s
(C = c0/n with n is the refraction index
and co = 2998.108 m∙s−1)
d Particle diameter, m
E The matrix identity
k Extinction coefficient, m−1
L Luminance, W/(m2Sr1)
n Refraction index
T Temperature, K
σ Boltzmann constant, σ = 1.3805 × 10−23 J∙K−1
δ Scattering coefficient
εp Emissivity of the solid
λ Thermal conductivity of the gas, W/(m∙K)
λr Radiative conductivity
μ Direction cosine between the directed intensity and the positive t axis
σr Stefan-Boltzmann constant, σr = 5.65 × 10−8 Wm2∙K4
h Planck constant, h = 6.6255 × 10−34 J∙s
ρ Density of the gas, kg/m3
ρs Density of the solid, kg/m3
τ Optical thickness, m
Incident radiative flux (solar energy), W
Φ Radiative flux density, , W/m2
ω Albedo, the ratio of the scattering to the Extinction coefficient, ω = δ/k
 Boles, M.A. and Ozisk, M.N. (1972) Simultaneous Ablation and Radiation in an Absorbing, Emitting and Isotropically Scattering Medium. Journal of Quantitative Spectroscopy and Radiative Transfer, 12, 839-847. https://doi.org/10.1016/0022-4073(72)90072-6
 Belghit, A., Daguenet, M. and Reddy, A. (2000) Heat and Mass Transfer in a High Temperature Packed Moving Bed Subject to an External Radiative Source. Chemical Engineering Science, 55, 3967-3978. https://doi.org/10.1016/S0009-2509(99)00575-8
 Belghit, A. (2007) Heat Transfer by Simultaneous Radiation-Conduction and Convection in a High Temperature Packed Bed. Eurotherm 81, Reactive Heat Transfer in Porous Media; Ecole des Mines d’Albi, France.