USE OF THE MULTIGRID METHODS FOR HEAT RADIATION PROBLEM

We consider the integral equation arising as a result of heat radiation exchange in both convex and nonconvex enclosures of di ﬀ use grey surfaces. For nonconvex geometries, the visibility function must be taken into consideration. Therefore, a geometrical algorithm has been developed to provide an e ﬃ cient detection of the shadow zones. For the numerical realization of the Fredholm integral equation, a boundary element method based on Galerkin-Bubnov discretization scheme is implemented. Consequently, multigrid iteration methods, which are closely related to two-grid methods, are used to solve the system of linear equations. To demonstrate the high e ﬃ ciency of these iterations, we construct some numerical experiments for di ﬀ erent enclosure geometries.


Introduction
Radiative heat exchange plays a very important role in many physical situations.The physical principles of heat radiation are very well understood, and there exists a number of engineering and physics textbooks where a whole hierarchy of different radiation models have been presented (see, e.g., [2,4,7,9]).On the other hand, papers dealing with questions related to heat radiation equation or its numerical realization tend to focus mostly on the simplest possible case of heat radiation exchange in convex enclosures [1,3,8].In fact, one of the most interesting features about transport of heat radiative energy between two points in both convex and nonconvex enclosures of diffuse grey surfaces is its formulation as an integral equation.One of the consequences of this fact is that the pencil of rays emitted at one point can impinge another point only if these two points can "see" each other, that is, the line segment connecting these points does not intersect any surface.The presence of the visibility zones should be taken into consideration in heat radiation analysis whenever the domain, where the heat radiation transfer is taking place, is not convex.
The computation of these visibility zones is not easy, but we were able to develop an efficient algorithm for this purpose and implement it in our computer programme.In dealing with the numerical aspect of this problem, we use the boundary element method based on Galerkin-Bubnov discretization of the boundary integral equation.This leads to a system of linear equations that will be solved iteratively using multigrid methods, which are closely related to two-grid methods.In fact, multigrid methods are among the most efficient methods for solving the linear systems associated with the numerical solution of the heat integral equation.The characteristic feature of the multigrid iteration is its fast convergence.The convergence speed does not deteriorate when the discretization is refined, whereas classical iterative methods slow down for decreasing grid size.As a consequence, one obtains an acceptable approximation of the discrete problem at the expense of the computational work proportional to the number of unknowns, which is also the number of equations of the system.It is not only the complexity which is optimal but also the constant of proportionality is so small that other methods can hardly surpass the multigrid efficiency.Numerical examples are considered here to demonstrate the high performance of these iterations.

Radiation on diffuse grey surfaces
We consider an enclosure Ω ⊂ R 2 with boundary Γ equivalent to the situation of Figure 2.1, where Ω is a conducting body.We assume, for simplicity, that the temperature on Γ is known.
The heat balance on Γ reads as where Q is the heat brought to the surface by conduction, q denotes the radiation emitted by the surface Γ, and J is the energy of incoming irradiation on Γ.For surfaces that are diffuse and grey as emitters and reflectors, the intensity of emitted radiation has the representation (see, e.g., [7]) where ε is the emissivity coefficient (0 < ε < 1), σ is the Stefan-Boltzmann 2) can then be written as In the geometry of Figure 2.1, the irradiation on Γ depends only on the radiation emitted by different parts of Γ itself.Hence, for any point x ∈ Γ, we can write J(x) = Γ F(x, y)q(y)dΓ y . (2.4) The substitution of (2.4) into (2.3)yields the integral equation where the kernel F(x, y) denotes the view factor between the points x and y of Γ.For convex two-dimensional enclosure geometries, F(x, y) has the representation [10] with If the enclosure is nonconvex (Figure 2.1), then we have to take into account the visibility function In this case, the kernel F(x, y) in (2.6) takes the form (2.9) ) is a Fredholm integral equation of the second kind.We introduce the integral operator Some properties of the integral operator (2.11), along with the solvability of (2.5), have been investigated in [10,11].

Construction of the system of equations
The Fredholm integral equation (2.5) can be expressed as where Kq = (1 − ε) Kq and For the numerical simulation of our integral equation, we use the boundary element method.We consider a Galerkin-Bubnov formulation and choose the basis trial functions Φ k (t) with local support Γ k ⊂ Γ.The approximation solution has the general form We let Inserting the ansatz function By introducing the vectors a = (q k ) k=1,...,n and b = g, Φ l,n Γ , l = 1, . . ., n, the matrices M = (M l,k ) l,k=1,...,n , with and S = (S l,k ) l,k=1,...,n , with (3.8) The mass matrix M in (3.8) is symmetric, positive definite, and diagonal dominant.Hence it is invertible.Consequently, (3.8) can always be written in the form To express the fact that the discrete equation (3.9) corresponds to the continuous equation (3.1), we write (3.9) as where q n = a n , g n = M n −1 b n , and

The hierarchy of discrete problem
In general, the discretization parameter n in (3.10) determines the dimension of the matrix system.For the two and multigrid methods, we use the hierarchy in multilevels.Let n 0 be a fixed discretization parameter and h 0 the corresponding step size.By successive halving, we obtain the step sizes More generally, one can consider an arbitrary step size hierarchy The index l is called the level or level number.The step size h l is associated with the dimension parameter n l .Hence, (3.10) at level l is To avoid a double indexing, we write (3.13) as q l = g l + K l q l (l = 0, 1, . . .). (3.14)

Iteration schemes for (3.14)
In this section, we try to show how multigrid methods can be implemented to solve the linear system (3.14).

Picard iteration
One can solve the system of equations (3.14) iteratively.The simplest iteration is the Picard iteration.The (i + 1)st iterate q i+1 l is obtained by inserting the ith iterate q i l into the right-hand side of (3.14):

.16)
A sufficient convergence condition is the matrix norm estimate K l < 1. (3.17)

Two-grid iteration
This method consists of two major steps.The Picard step (the so-called smoothing step since the error is smoothed) and the coarse-grid correction ql → ql − Pδ l−1 .In fact the usual procedure of the two-grid iteration of level l for one iteration step q i l → q i+1 l can be illustrated in Algorithm 3.1.

.19)
Stopping criterion: Coarse grid correction: Here r is n l × n l−1 restriction matrix and P is n l−1 × n l prolongation matrix.The indices l − 1 and l are used for the coarse grid and fine grid, respectively.

Convergence of the two-grid method
Since the mapping q i l → q i+1 l of the two-grid algorithm is affine, it has the representation [5,6]

.24)
where A TGM l is the two-grid iteration matrix that has the representation

.25)
A sufficient condition for the convergence of the two-grid method is the matrix norm estimate (3.26)

Multigrid iteration
Even though the two-grid iteration reduces the amount of computational work drastically, the solution of the coarse-grid equation (3.22) still takes the major part of the work.The problem to be solved in (3.22) reads respectively.Obviously, (3.27) has the same form as the original equation (3.15) which is solved by the two-grid method at levels l − 1 and l − 2.
Then it becomes necessary to solve an auxiliary equation Again, the two-grid algorithm can be applied to levels l − 2, l − 3, and so forth .The resulting algorithm known as the multigrid iteration uses all discretization levels.Such iteration is given by Algorithm 3.2.

Convergence of the multigrid method
By analogy with (3.24), we write one step q i l → q i+1 l of the multigrid algorithm in the form where A MGM l is the multigrid iteration matrix that is recursively defined by (see [5,6]) An alternative representation of (3.34) is (3.35)A sufficient convergence condition is the matrix norm estimate (3.36)

Realization of the visibility function v(x, y)
We represent here an efficient algorithm for the computation of the visibility function v(x, y) given in (2.8).To satisfy the first requirement, we set x = X(t 0 ) and define Next, determine t 1 = t 0 + .When z Γ = X(t 1 ), the assertion follows immediately.

Numerical experiments and results
Since the convergence requirements (regularity, consistency, and stability) for the two-grid and multigrid iterations are satisfied [10], we can now apply these algorithms to solve the linear system of equations (3.14).For the numerical application, we choose the emissivity coefficient as ε = 0.2, the Stefan-Boltzmann constant has the value σ = 5.6696 × 10 −8 (W/m 2 K 4 ), and the surface temperature will be given by the function with T 1 = 1000K and T 2 = 1800K.The mass matrix M = (M l,k ) l,k=1,...,n and the right-hand side b n = g, Φ l,n Γ , with g(t) = εσT 4 (t), either can be calculated analytically exact for special geometries or numerical integration is applied.The computation of the stiffness matrix S = (S l,k ) l,k=1,...,n = KΦ k,n , Φ l,n Γ has been performed numerically using Gauss quadrature.Theoretically, Galerkin method requires a time-consuming double integral over Γ for the calculation of every element of this stiffness matrix.Due to this fact, we have used the Gauss quadrature with respect to fast computation, that is, by evaluating the kernel of the integral equation as seldom as possible.Theoretical and numerical error estimates for this problem have been presented in [10].
Moreover, in the case of a nonconvex enclosures, the main problem is the efficient detection of the shadow zones to calculate the visibility function v(x, y) appearing as a part of the stiffness matrix S. Thus, a geometrical algorithm was developed (see Section 4) to determine the visibility function in two dimensions.To this end, we consider the following examples.seconds.Note that the step size h l is associated with the dimension parameter n l , where h l = 1/n l with n l = 2 l and l is called the level number.
One sees clearly in Table 5.1 the efficient performance of the two-grid and multigrid techniques for solving this problem in regards to the number of iteration steps and CPU time required to achieve fast convergence.

Nonconvex enclosure
Example 5.2.As an example of the nonconvex enclosure, we consider the curve shown in Figure 5.1.In this case, the visibility function v(x, y) must be taken into consideration and consequently its geometrical algorithm is implemented in our computer programme.Computation of this visibility function has been presented in Section 4. Table 5.2 shows the numerical results for this nonconvex case.In fact, one concludes similar remarks to those reported in Table 5.1.

Figure 2 . 1 constant
Figure 2.1 ) then (3.5) can be rewritten asM n − S n a n = b n .

5. 1 .
Convex enclosure Example 5.1.Let Ω be the domain of a unit square.Table 5.1 shows the numerical results for solving (3.14) by the two-grid and multigrid methods.It contains both the number of iteration steps and the CPU time in t = 0 t = 0.25 t = 0.5 t = 0.75