1. INTRODUCTION
A numerical approximation of the radiative transfer equation usually
involves macroscopic methods in regions where the medium is opaque,
coupled with methods for solving the photon transport kinetic equation
in a transparent medium. In an opaque medium, the radiative intensity
is close to the radiative equilibrium. The diffusion approximation or
the P1 method, based on the hypothesis of a
directional quasi-isotropy of photon emission distribution provides a
rather accurate approximation (Mihalas &
Mihalas, 1984; Pomraning, 1992). For
the outside diffusion regime, in the case of transparent media, the
macroscopic methods are no longer valid. The microscopic description of
the transfer equation is then preferred even if it is complicated and
time consuming. Also, coupling between diffusion and kinetic solutions
in intermediate regions is very complicated.
The author purposed in previous works a new macroscopic hierarchy of
models that are capable of restituting high radiative nonequilibria
(Dubroca & Feugeas, 1999). The method is
based on the entropy minimization principle applied to the radiative
transfer equation. This idea comes from the approximation of
transitional flows at high altitude introduced by Levermore (1996) for the Maxwell–Boltzmann theory. The
first step of this idea, introduced by Grad (1949), lies on the building of new, continuous
equations taking account additional unknowns. These additional
equations are obtained by considering an adequate closure. However this
approach does not assure the hyperbolic properties of the system.
Levermore (1996) purposed another closure
that ensures the hyperbolicity of the moment equations systems and
dissipate the entropy locally. This approach appears as a general
method to close an infinite hierarchy of macroscopic equations
constructed by taking the moments of a microscopic equation (Feugeas, 1997; Charrier et
al., 1998).
In previous work (Dubroca & Feugeas,
1999), we applied the idea of the entropy minimization principle
to the Bose–Einstein theory to build a hierarchy of models for
the radiative transfer, the M-hierarchy. Each member of this hierarchy
is a hyperbolic system that has a flux limited by the speed of light
and which dissipates entropy locally. The limiting property ensures
that any signal is propagated at a velocity below the speed of light,
contrary to diffusion models that allowed arbitrary and non physical
high velocities. The first member of this hierarchy, the
M1 model, is very attractive because of its
analytic formulation. A first version of this model has been analyzed
by Fort (1997) in the simple case where the
radiative flux is always parallel to a given axis. It has been verified
that the M1 model describes a rather wide class of
situations where anisotropy can be characterized by only one direction.
This model has been used in several physics domains (Audit et al., 2002). In order to enlarge
the capabilities of the M1 model to describe more
complicated anisotropies, several modifications are purposed (Dubroca & Klar, 2002; Klar
& Dubroca, 2002; Klar et al.,
2003). However those are not satisfactory because of losing
interesting properties of the hierarchy. At first, one of the very
attractive characteristics of the M1 model is the
analytic formulation. Besides, every member of the M-hierarchy is built
in the respect of the Galilean invariance of the underlying
distribution. The other members of this hierarchy have no known simple
explicit formulation. However, they conserve other properties of the
M-models: Their flux is limited, they are hyperbolic, and they locally
dissipate the entropy. In the present article we explore the anisotropy
potentialities of the M2 model.
Having recalled the main properties of radiative equilibrium and
entropy in Section 2, we obtain a hierarchy of moment systems as
presented in Section 3. In Section 4, we analyze in detail the first
member of the M-hierarchy, the M1 model. Section 5
addresses the properties of the M2 model. Section 6
is devoted to this anisotropy analysis of the two first members of the
hierarchy. Section 7 concludes the article and proposes some
perspectives.
2. RADIATIVE TRANSFER: EQUILIBRIUM AND
ENTROPY
2.1. Radiative transfer equations
The radiative transfer equation governs the radiative field and its
interaction with matter. It describes the evolution of the
radiative intensity I coupled to the matter energy equation.
For purposes of simplification, we shall take the hypothesis of the
local thermodynamic equilibrium (LTE) and assume that the medium is at
rest. In order to simplify the description, we shall subsequently
assume that the absorption coefficient is independent of the frequency.
This hypothesis is not essential, but it lightens the reasoning pursued
in this article. For the same reason, the scattering effects will also
be neglected. The radiative transfer equations are then written as
follows:
Here, I = I(r,ν,
Ω,t) is a function of the position
r, the frequency ν, the unitary propagation
direction Ω, and of the time t, c
is the speed of light in a vacuum; Tm is
the matter temperature, ρ is the matter specific density, σ is
the absorption coefficient, and Cv is the
specific heat at a constant volume. The radiative collision operator
σ(B(Tm) − I)
establishes a balance between the emission and absorption of the
photons and couples the two equations; S2 is the
unit sphere. The function B =
B(ν,Tm) is the Planck
distribution defined by
where k is the Boltzmann constant and h is the
Planck constant. j =
σB(Tm) represents the matter
emission isotropic source at the temperature of
Tm. The moment of a scalar or
vector-valued function g = g(ν, Ω)
in the space of the frequencies and directions is written
The radiative energy is defined by
ER = 〈I〉, the
radiative flux vector by FR =
〈cΩI〉, and the radiative
pressure tensor by
.
The normalized flux is written as f =
FR /cER.
In the particular case of a Planck radiative intensity, we have
ER(Tm) =
aTm4,
FR(Tm) = 0,
and the pressure tensor reduces to the scalar
PR(Tm) =
aTm4/3 with a =
8π5k4/15h3c3.
In the case of a given radiative intensity, the radiative temperature
TR is defined as from the radiative
energy: ER =
aTR4. In the case of a Planck
radiative intensity, naturally we have TR
= Tm.
2.2. Radiative entropy and equilibrium
If the matter is supposed to behave as a perfect gas, the total
entropy of the system (1) and (2) is the sum of the matter entropy and
of the purely radiative entropy,
HR*(I) =
〈hR*(I)〉, where
hR*(I) is a strictly convex
function (Fort, 1997):
where n is the occupation number density defined as
I =
(2hν3/c2)n .
The total entropy of the system is then
(HR*(I) +
S)/c where S = −ρ
Cv log(Tm) is
the matter entropy. We have the following properties on the entropy of
the system. The total entropy of model (1) and (2) is locally
dissipated (decreases in time):
Moreover, at the thermodynamic equilibrium, T =
TR = Tm, the
Planck distribution (3) or the so-called equilibrium distribution,
produces the minimum radiative entropy under the constraint of a good
reconstruction of the energy:
This problem of minimization is equivalent to solve the Lagrangian
equation:
where α is the Lagrange multiplier vector.
The constraint 〈I〉 =
ER being achieved by definition, the solution
to this minimization problem must satisfy, for all I, the following
condition:
The unique solution to this equation is solution of
in other words
The value of α is determined by the constraint
〈I〉 = aTm4.
Thus, α = 1/Tm c and
Tm = TR at
the Planckian equilibrium.
It is important to remark that relation (7) allows us to show, before
any analytic computation, that α > 0 (because n is
always positive). This result is very important because it insures that
B is strictly positive. We will see further that this property
is true for every member of the M-hierarchy.
3. HIERARCHY OF MOMENT MODELS
Moment models were introduced by Grad (1949) and developed by Levermore (1996; see also Feugeas,
1997; Charrier et al., 1998).
Their construction starts by the choice of a vector-valued subspace,
,
of a finite dimension m of functions of Ω with real values.
is usually chosen as being a space of Legendre polynomials or of the spheric
functions. For a given m, the vector m is composed of the
basic functions of
.
For example, for m = 1, the vector m = (1,
Ω). This is a combination of a scalar and a vector.
Similarly, for m = 2, m = [Ω,
Ω [otimes ] Ω]T.
3.1. Closure of moment equations
Considering the moments of the radiative transfer equation multiplied
by the vector m leads to the construction of radiative moment
equations:
The moment vector is written as ER =
〈mI〉. An additional relation is now needed
to close the chosen moment model. The choice of closure relation is not
unique, and is subject to only one constraint: The radiative intensity
of the model should be nonnegative. As was shown by Levermore (1996) for the Maxwell–Boltzmann theory, the
closure of the radiative moment model (9) can be obtained by finding
the distribution that minimizes the entropy under the constraint of a
good construction of the moments generated in the vectorial space
(Muller & Ruggeri, 1993; Fort, 1997; Dubroca &
Feugeas, 1999):
This minimization problem is equivalent to solve the following
Lagrangian equation:
where α is the vector of the Lagrange multipliers. This
vector is defined such that α·m > 0,
which insures the positiveness of the solution
.
The solution
to this minimization problem must satisfy, for all I, the following
condition:
The solution to this equation is
which leads to
Here, α is determined by the constraint
〈mI〉 = ER.
Besides, the previous relation (11) insures that
α·m > 0, because n > 0.
This result is very important because it guarantees that, for every
member of the hierarchy, B > 0.
The moment models are therefore written in the following form, if
is given by (8):
where
is given by (12).
3.2. Radiative moment models
The hierarchy of the moment models (13) thus generated, offers
interesting structural properties: The normalized flux f is
limited (∥f∥ ≤ 1); the system (13) is hyperbolic;
and the system (13) locally dissipates the entropy.
The first property is extremely important as it guarantees that the
propagation of a perturbation in the radiative medium is limited by the
speed of light. This property is not satisfied in the classical
diffusion, or Eddington, description (P1 model;
Castor et al., 1981; Pomraning, 1992). It lies on the basic assumption that
the angular dependence of the specific intensity is represented by the first
two terms in the spherical harmonic expansion (I =
I0 + Ω·I). This specific
intensity I is then almost isotropic as the angular term is supposed
to be a small perturbation ∥I∥ <<
I0. If this condition is not respected, it leads to
nonphysical negative emission in the flux direction. More generally, the
Pn models fail to respect the property of
limitation of the normalized flux, and require the introduction of a
procedure for the limitation of the radiative flow in the nonopaque
zones.
In turn for the M-models, the second and third properties guarantee
that the system (13) has been properly set mathematically. Limitation
properties on the radiative pressure
tensor can therefore be established. If the normalized flux f
and the Eddington tensor
are moments of the radiative intensity
defined by (12), they satisfy the following constraints:
This results in the flux limitation. For an extreme non-equilibrium,
one finds
.
4. THE FIRST MODEL OF THE HIERARCHY
The M1 model is the first of the hierarchy, and
it corresponds to the vector-valued space,
,
spanned by the vectors 1 and Ω. This is the smallest
vector-valued subspace to restitute an angular anisotropy. The distribution
function has the following form:
The vector A and the constant B are defined by the
conditions
:
The details of this calculation are presented in Dubroca and Feugeas
(1999). The radiative M1
model can now be closed, and the tensor
can be expressed analytically as a function of the moments
ER and FR of
the model:
where
is the Eddington tensor, calculated as a function of the Eddington factor
χ and n = f/∥f∥:
The Eddington factor is a function of ∥f∥:
It is an increasing function of ∥f∥ for which χ
≤ 1 and χ(0) = 1/3. The intrinsic values of the hyperbolic
M1 model at equilibrium (when f = 0) are
equal in absolute value to
,
and have the opposite signs. This situation clearly illustrates the
emission isotropy of the photons that prevails at radiative
equilibrium. On the other hand, in the case of extreme nonequilibrium
(when ∥f∥ tends toward 1), the eigenvalues of the
Jacobian tend towards the speed of light in absolute value and are of
the same sign. In the case of extreme nonequilibrium (the angular
distribution of the radiative intensity moves towards a
δ-function), it is clearly found that the photons all move at the
same speed and in the same direction.
The M1 model, like every member of the
M-hierarchy, has many desirable properties. It has been used in several
domains (Audit et al., 2002).
However, this first model of the M-hierarchy is not able to describe
more than one beam anisotropy. To capture more anisotropic effects,
various modifications of the M1 model are purposed
(Dubroca & Klar, 2002; Klar & Dubroca, 2002; Klar
et al., 2003). However they lose interesting properties
of the hierarchy. For example, the M-hierarchy models respect the
Galilean invariance property of the underlying distribution. Besides,
one of the very attractive characteristic of the first member of the
hierarchy is the explicit analytic formulation. The other members of
this M-hierarchy have no known simple explicit formulation. However,
they conserve all the interesting properties of the M-models
(hyperbolicity, flux is limited, entropic, Galilean invariance). For
this reason, I propose here to explore and describe the potentialities
of the second model of the M-hierarchy.
5. THE M2 MODEL
The M2 model corresponds to the vector-valued space,
,
spanned by the vector Ω and
the tensor Ω [otimes ] Ω. The distribution
function is then expressed as follows:
where the vector A and the tensor
are defined by the conditions :
.
The distribution
achieves the minimum radiative entropy under the constraint of a good
reconstruction of the moments
:
where
is the vector of the Lagrange multipliers. m = [Ω,
Ω [otimes ] Ω]T. As
we have shown, α is defined such that, for any
Ω:
which insure that
.
For any Ω,
That means that C is a symmetric positively defined matrix
for any Ω.
One can readily see that the pressure relation recovers also the
energy relation because Tr(Ω [otimes ]
Ω) = 1. That is the reason why the M2
model only uses those two relations:
where the tensor
is defined as
.
We can note that M2 required the resolution of
nine equations in the three-dimensional case and five equations in the
two-dimensional case. Besides, the M2 model
formally recovers the M1 model. However, one cannot
obtain an analytic formulation of this closure. This is not a big
constraint to solving such models. Kinetic schemes offer very elegant
and effective solutions (Dubroca & Klar,
2002).
6. ANALYSIS OF THE NONEQUILIBRIUM SOLUTIONS
The specific intensity I is a function of the position
r, of the frequency ν, and of the unitary propagation
direction Ω. For a given coordinate system, this vector
can be written as a function of the angle θ: Ω
=T [cos(θ),sin(θ)] . Then, for
a given temperature, for a given point of the space, the normalized
radiative intensity can be represented in cylindrical coordinates (I,
θ, ν). Figure 1b,c presents two
views of this function. For a given θ, the distribution
I(ν) is always a Planck distribution (Fig. 1b) with a temperature depending on θ.
Because the functional form of I(ν) is given, it is
sufficient to plot an angular dependence for an arbitrary given ν.
Therefore, the polar representation provides a useful tool to analyze
the anisotropy of emission (Fig. 1a).
Specific intensity I(θ,ν) for a given temperature:
(a) two-dimensional polar coordinates I(θ) for a given
frequency; (b) I(θ,ν); (c) cylindrical coordinates
(I, θ, ν).
6.1. Analysis of the M1 model
In an isotropic configuration, the underlying intensities, for the
diffusion approximation P1, and for
M1 are identical (Fig.
2a). This rapidly ceases to be the case with a normalized flux
of the order of 0.3, and the difference increases for
∥f∥ = 0.5 (Fig. 2b). When
the normalized flux is of the order of 0.9 (Fig.
2b), that is, in the case of high anisotropies, a non-physical
behavior of the P1 model is evidenced. It would
appear that in this case the radiative intensity is positive in the
direction of the flux, that is, in the direction of highest energy. But
the specific intensity is negative in the opposite direction of the
flux. That indicates that over a certain limit the underlying diffusion
intensity becomes negative. On the contrary, the underlying
M1 intensity remains physically admissible for any
flux.
The M1 model: the underlying distributions of the
M1 model compared with the P1
one: The dashed line correspond to the diffusion approximation
(P1) and the full one corresponds to the first
member of the M-hierarchy (M1). a: Planckian
equilibrium. b: Medium anisotropy. c: Strong anisotropy; the diffusion
approximation distribution can be negative with nonphysical sense. d:
Extreme equilibrium: In the case of the light ray, the diffusion is no
longer available.
When the norm of the normalized flux tends toward 1 (Fig. 2d), the maximum physically acceptable
anisotropy, the difference between the two models increases. The
M1 model tends toward a distribution carried solely
by the radiative flux (a light ray), whereas the model
P1 fails to adopt any special direction. The
P1 model, contrary to the M1
model, does not have any intrinsic flux limitation. For example, for
∥f∥ ≥ 1, the P1 model is still
defined and starts to present a strong region of negative radiative
intensity. This analysis confirms the capacity of the underlying
radiative intensity of the M1 model to capture
strong one-ray radiative nonequilibrium. For a maximum physically
acceptable anisotropy, the M1 model tends toward a
distribution carried solely by the radiative flux (a light ray),
whereas the P1 model does not opt for any direction
in particular and predicts a negative radiative intensity in the
direction of the flux. This model is not able to describe more than a
one-ray regime. The second model of the M-hierarchy offers more
anisotropic potentialities.
6.2. Analysis of the M2 model
At first, as was expected, the underlying distribution of the
M2 model, confirms its ability to recover formally
the M1 description. That is the case for weakly
anisotropic situations (Fig. 3a) and also for
extreme one-ray angular distributions (Fig.
3b). However, this model has much greater potentialities. Figure 4 presents the symmetric bipolar situation
that the M2 model is able to represent. The matrix
has a diagonal asymmetry in the first example (Fig.
4a), and a off-diagonal component for the other one (Fig. 4b). Figure 5 introduces
nonsymmetric bipolarity characteristics associated with the energy flux.
Figure 6 shows various two-ray emission
configurations. Those kind of emissions could be longitudinal (Fig. 6a), oblique (Fig.
6b), and nonsymmetric (Fig. 6c).
The M2 model: The one-ray configuration shows
that the second member. If this hierarchy formally recovers the
M1 model for medium dequilibrium (a) or for a
stronger one (b).
The M2 model: the symmetric bipolarity
characteristic in the main (a) and secondary (b) directions; this could
be the case of the emission seen by a point between two sources, ton
their axe.
The M2 model: symmetric bipolarity
characteristic; this could be the case of the emission of two sources
seen from a point near (a) or far (b) from their axes or of two far
sources (c).
The M2 model can represent various extreme
two-ray emission configurations; (a) symmetric and longitudinal; (b)
symmetric; (c) nonsymmetric.
The analysis of the emission anisotropic potentiality of the
underlying radiative intensity of the M2 model
confirms its capacity to capture strong and various anisotropies. Those
strong potentialities show that this model could be used in the
inertial confinement fusion (ICF) domain.
Let us consider, for example, the case of an emissing sphere (Fig. 7a). Usually, this kind of configuration
cannot be represent by classical diffusive models or by the
M1 model. A punctual source, very far from point M,
could be represent by this first model. However, when M get closer from
the spherical source, the emission becomes bipolar as we can see in the
Figure 7 in the Cartesian description (b) or
in the polar one (c). This configuration can of course be represented
by the second model of the hierarchy (Fig.
7c). This ability to get such bipolar situations can be extended
in the multidimensional cases. Besides, if the source loses its
symmetry, the second member of the hierarchy is able to take it into
account. To illustrate the range of potentiality of this model, we
could also choose the case of a punctual intense source in a gray
middle, or the case of two far and punctual sources. For those cases,
the M2 model can get good solutions.
M2 model: Application of a spherical emissing
source (a); for a point M close from the spherical source, the specific
intensity becomes bipolar as we can see in the Cartesian description
(b) or in the polar one (c).
More generally, this second member seems to be able to get a solution
for the treatment of boundary conditions between gray and transparent
middle.
7. CONCLUSIONS AND PERSPECTIVES
The models of this M-hierarchy offer very attractive properties. The
hyperbolicity allowed us to use effective numerical tools. The
limitation of the flux, and the entropy dissipation insure the physical
credibility of the modelization. The construction of those models
respects particularly the Galilean invariance of the underlying
specific intensity. The first model (M1) is very
effective in the one-ray cases. One very interesting and useful
property of this first member is the analytic formulation of its
closure. However, it can be used in its multigroup formulation (Turpault, 2002) that is not analytic. Besides,
various nonanalytic modifications also have been proposed by different
authors (Dubroca & Klar, 2002; Klar & Dubroca, 2002; Klar
et al., 2003) since the original one (Dubroca & Feugeas, 1999) to get more
anisotropy. The other members of this M-hierarchy have no known simple
analytic formulation of their closure. However, they all respect the
basic properties of the hierarchy and particularly the Galilean
invariance property of the underlying distribution. Besides, the second
model provides enough degrees of freedom to reach a large range of
Des-equilibrium regime. The resolution of this model requires solving
only nine moments equations (3D).
The main difficulty of the moments method is the
realizability of the distribution function: for any moments
set Em, one has to insure that it is
always possible to find a distribution (12). We can show that all
models of this hierarchy are realizable.
We can also note here another very useful characteristic of this
family. All models have the same hyperbolic structure. This property
allows us to use the same numerical tools for each one. Beside, the
boundary conditions treatment is one point very important for those
kinds of radiative regime. For example, we could have to couple
different kind of methods (continuous, kinetic) between gray and
transparent middle. If we want to exchange kinetic informations such as
the specific intensity, we have to respect the kinetic properties.
ACKNOWLEDGMENTS
The author thanks J.-F. Clouet, P. Charrier, B. Dubroca, G. Duffa, G.
Samba, G. Schurtz, and V. Tikhonchuk for numerous discussions that has
been very helpful for this study. He acknowledges particularly V.
Tikhonchuk for his very useful help in writing this article.
