Computer Physics Reports 351 (1987) 351-369

North-Holland, Amsterdam

FINITE ELEMENT METHODS

351

FOR THE TRANSPORT

EQUATION

Pierre LESAINT
Faculth des Sciences, Universitk de Besanqon, 25030 Besanqon Cedex, France

Continuous
and discontinuous
methods are described
Monodimensional
spherical geometry is then considered.
rectangular or triangular meshes.

for solving the transport

equation in x-y geometry.
We end up with the description of nodal methods on

0167-7977/87/$07.00
0 Elsevier Science Publishers B.V.
(North-Holland Physics Publishing Division)

352

P. Lesaint / Fimte element methods for the transport equution

Contents
..........
1. introduction
....................................
2. Position of the problem and discrete ordinate method. ......
3. Continuous and discontinuous finite element methods for the (S. ,;) &iki~
..........
geometry ......................................
..........
3.1. Continuous methods. ..........................
..........
3.2. Discontinuous methods. ........................
4. Continuous methods for the one-dimensional spherical problem ..........
..........
4.1. Continuous method ...........................
..........
4.2. Discontinuous method .........................
5. Nodal methods for the (x, y) planar geometry. ...........
..........
..........
5.1. Rectangular meshes ...........................
..........
5.2. Triangular meshes ............................
..........
References .......................................

353
353
357
358
361
362
363
364
365
365
366
369

P. iksaint

353

/ Finite element methods for the transport equation

1. Introduction
To solve the neutron transport equation, a commonly used procedure consists, in a first, step
to discretize with respect to the angles. We get a set of discrete ordinate equations consisting for
each angular direction in a first-order hyperbolic problem. This problem has for a long time been
solved on orthogonal
meshes by finite difference methods. For non-orthogonal
meshes, characteristic methods, finite elements or nodal methods may be used. The aim of this paper is to
present some finite element and nodal methods for the spatial discretization
of the transport
equation.
In the first section, we introduce the general problem and show how one can get the discrete
ordinates equations.
In the second section, we define continuous and discontinuous
methods for the planar x, y
transport equation. We then show how these methods apply to the monodimensional
problem in
spherical geometry.
In the last section, we define nodal methods on rectangular and on triangular meshes.

2. Position of the problem and discrete ordinate method

The time-dependent

neutron

transport

equation

is written

as follows,

for x E D c R3 and

UEVCR3

au

at

+ u grad,u(x,

k(x,
/V

u, t) + (J(x, u)u(x,

u, u)u(x,

u, t>

u, t) du+f(x,

(2.1)

u, t),

(24
where U(X, u, t) represents the particle flux at the geometrical point x and streaming with the
velocity u.
The integral at the right hand side of eq. (2.1) takes into account the scattering and the fission
properties. The kernel k(x, u, u) is non-negative
and, in the isotropic case, depends only on x,
on the modulus 1u 1 and 1u 1 and on the angle p = arccos (u, u)/ 1u 1 1u 1. The source term is
represented by f( x, u, t). In most of the cases, we have
I/=
Boundary
u(x,

{uER3;
conditions

a< Iu(

sb}.

may frequently

(2.3)
be written

as follows

u, t) = given function,

for all x E CID such that

domain D.

u. n -C 0, where

(2.4)
n denotes

the outer

normal

on the boundary

aD of

P. Lesomt

354

Eq. (2.1) written

in condensed

/ Finite element methods for the transport equation

manner,

reads as follows

l tLu=Ku+f,

(2.5)

where the operator L represents the spatial streaming of particules and takes into account the
total cross section a; the integral operator of the right hand side of eq. (2.1) is denoted K.
Assume that f is independent
of t. A necessary condition for the convergence as t + + CYZ
of
the solution u of problem (2.5) to the solution ii of the stationary problem
Lii=Ku+f

(2.6)

is that the eigenvalues

problems:

CIof the operator

i) To find the eigenvalue

K - L satisfy

Re (Y< 0. We have to solve the following

(Ywith bigger real part, solution

of

Lu + au = Ku.

(2.7)

ii) To solve, when Re LY< 0 for all the eigenvalues

(2.6).

CI of problem

To solve problem i) we first determine,

the spectral radius
values of (Y.This is done by a power inverse iterative method

(2.7), the stationnary

problem

X(a) = p((L + a)- K) for given

LU n++aun+*=Ku.

(2.8)

The value of (Ysuch that A( a) = 1 may then be found by a procedure

[8] that this eigenvalue cx is of maximal real part.
A usual method for solving problem (6) is the iterative algorithm
Lu r~+ =Ku"+f

of dichotomy.

It is shown

(2.9)

which converges when p( L- K) c 1. It is shown [8] that this last condition is satisfied when
Re (Y< 0 for all the eigenvalues (Y,solution of (2.7). To treat the dependence
of the solution ii of
problem (2.6) with respect to the velocity u, we first discretize the space V by considering
different modulus of the velocity of the particles (energy groups)
a I u1 <

v2 < . . . < vI _<h.

The flux u,( x, w) corresponding

to each energy group v, is then a function of the geometrical
position x and the two angular cosines represented
by w (it is sufficient to represent a vector
w E Iw with modulus equal to one). Denote by S the unit sphere in Iw3, we may write
w . grad,u;

+ u;u, = c
J

Js

k,, ( x, w, w)u/(x,

w) dw+f,(x,

w)

(2.10)

P. Lmaint

for 1 s i I I, which can also be written

L,ui=

To solve the integro-differential

1

as
(2.11)

1 Sill.

CK,jz4j+f,,

Lu!+l)=

35s

/ Finite element methods for the transport equation

system (2.11), we may use the following

CK,j~J+l+

iterative

algorithm

1lilI.

CKjjuS+fi,

(2.12)

j>i

j<i

This procedure is usually called external iterations.

For each iteration and for each energy group, we have to solve the following
(internal iterations)

(Lj - K;)uj+l) = Sj,

where S! is known
follows jS]

1I i

from previous

w. grad,u @+l)(x,

type of problems

(2.13)

I,

calculations.

w) + ~(x)u+r(x,

Dropping

w) = /k(x,
s

away the index

w, ~)u@+r)(x,

i, eq. (2.13) reads as

w) dw+f(x,

0)
(2.14)

which leads for each iteration

we grad,u(x,

to the following

w) + a(x)u(x,

For two-dimensional

planar

p$ +vg+

geometry,

A discretization
with respect
compute the term:
/_l, V) = F(x,

+ v$

aY

eq. (2.15) becomes

w) =

to the angular

Js

k(x,

A truncated
expansion of Legendre
discrete ordinate formulation
pg

(2.15)

CJ).

F(x, y, p, v).

(724 =

F(x,

w) = F(x,

problem

+ a+,

Y, P,?

direction

o, w)u()(x,
polynomials

%J =F(x,

(p, V) (with p2 + y2 I 1) is then needed

w) dw+f(x,
[3] is often

to

0).
used,

leading

to the following

Y, PL,, v,),

1 I m I M = N( N + 2)/2 for the S, approximation

32).

(for even values of N generally

not exceding

356

P. Lesaint / Finite clement methods for the transport equation

It remains

to treat the spatial variables

j.~$j+v*+f~~=I:
dY
u(x,

for

y) = given value,

for the generic problem

(x, .y)ED,

on

a-D=

(2.16)

{(x,

y)EaD;

where n = (n,, nv) denotes the outer normal on aD.

In the case of the one-dimensional
spherical geometry,

$+_-

1 -/AZ au

+au=F,

O<r<R,

-1

co}.

pn,+~n,

(2.17)

we get

<p-c<.

(2.18)

a/J
u (R , p) = given function

for

with p==.._c/l_~l.
Eq. (2.18) may also be written

p$(ru)+r&((l- p)u)

j_~cO,

as follows

+ aru

= r*F.

(2.19)

Before starting with the discretization

of problems (2.16) and (2.18), we give some fundamental
properties of those problems. Integrating eq. (2.16) over the domain D, we get

Jn
where a + D = dD - (a - D), this relation

[31 (property
Let d/dy
+ v au/ay),
*+ou=F,
dv

showing

the conservation

F dx d-v,

of particles

as expressed

in

1).
denote the derivative
we get

along the characteristic

direction

(p,

v) (i.e. du/dy

= p au/ax

(2.20)

which means that the solution u at the point (x,,, y,,) depends only on the values of u and F at
line going through the point (x,, y,,) and lying
the points (x, y) belonging to the characteristic
upstream this point (property 2). The numerical methods described later on shall satisfy these
two properties, at least in a discrete sense.
For the one-dimensional
spherical problem, we define

P. Lesaint / Finite element methods for the transport equation

Fig. 1. Spherical

geometry

351

in one dimension.

Integrating eq. (2.19) over the domain [R,,R 2]X [- 1, + l] we get the following relation for the
conservation of particules

R;m,( R2) -

R:m,(

R,)

+ ~Rzc+%z,,(r)
dr = ~~*/_;rS

dr dp.

RI

Let d/dy denote a derivative along the characteristic lines defined by v2(1 - ,u*) = cste (i.e.
(d/dy) 24= /I, au/% t- ((1- jJ2)/r) ~~/~~) we still get relation (2,20) and the solution ECat any
point MO= (rO, pO) depends only on the values of U, cf and F along the characteristic line going
through M,, upstream from M, (fig. 1).

3. Continous and discontinuous finite element methods for the (x, y) planar geometry
In this section, we assume that the domain D is convex. Assume that the domain
subdivided in triangles and (or) convex quadrilaterals, denoted in the sequel by T.
Let a - T (resp. a + T) = {(x, y) E i3T, pn,,,+ vnY,T<: 0 (resp. > O)}, where
(n X,T, ny,T) denotes the outer normal on G?

D is
nT=

P. Lesaint / Finite element methods for the transport equation

358

Fig. 2. Order of sweeping

To calculate the approximate

solution uh on the domain D, we set up the following principle,
which allows to satisfy property 2 in a discrete sense.
Calculations
may be done on an element T c D if we already know the
value of the
approximate
solution uh along the lighted part a - T of the boundary L3Tof the element T.
Such a procedure is always possible [14], the ordering of the elements T being non-unique (fig.
2) and algorithms for computing the order of sweeping a general mesh are available [16,17]. This
procedure is quite obvious for an orthogonal mesh, but for a general mesh the algorithm must be
run, and the order of sweeping must be stored, for each discrete angular direction prior to
computing
the solution. In the sequel, we assume that the approximate
solution u,~ to be
computed on the element T is already known on the lighted part 8 - T of the boundary of T.
We let P(k) and Q(k) be, respectively, the spaces of polynomials
p and 4 defined by

P(%

Y) =

a&Y-

It/Sk

4(x,

Y)

arJxY.

3.1. Continuous methods

Let P, be the space of shape functions,
au/ay + UU.
T. Let Hu=p
au/&x+v

representing

the approximate

solution

on the element

P. &saint

/ Finite element methods for the transport equation

359

;J---j;;rja,
as

a4

(It, v)

Fig. 3. Finite difference

The approximate

solution

J(Hu,-F)udxdY=O

uh E P, is defined
forall

UEQ~,

schemes.

by

(3.1)

u,, is known at the degrees of freedom 2 of the space P, associated with the part a - T of the
boundary. The dimension of the space of test functions Q, is equal to dimension P, - dimension
_E. To satisfy property 1 (i.e. conservation of the particles) we assume that the space Qr contains
the constant function.
Example 1 (fig. 3): Let T be a rectangle and P,= Q(1).
The degrees of freedom are the values of p at the corners of the rectangle. If the characteristic
directions (p, V) is not parallel to any side of the rectangle, we have dim 2 = 3, dim Q, = 1 so
that we write the equation
((Hu,-F)dxdY=O,

u,EP,

JT

which is exactly

the standard

(3.4

SN scheme [3]

Example 2 (fig. 3): Let T be a rectangle and P, = P(1).

The degrees of freedom are the values at the midpoints a, 1 I i I 4 of the sides of the rectangle
and the diamond relation is satisfied, i.e. p(q) +~(a,) =~(a,)
+~(a,).
We have dim .Z= 2,
dim QT = 1 and we still write eq. (3.2), ui + u3 = uz + u4 = 224, which is the well-known
diamond scheme,

360

P. Lesuint / Finite element methods

for the

transport equation

Fig. 4.

For these two examples, the error is of order h. with h = max( a.~, h-v).
Exumple 3 (fig. 4): Let T be a triangle A,A,A,.
Depending
on the characteristic
direction (p, v), we may have one or two lighted sides. Let
P, = P(l), the degrees of freedom being the values of the flux at the vertices of the triangle. If
one side of the triangle is lighted, then dim 2 = 2 and we write one equation on the triangle.
If two sides are lighted, then dim 1: = 3, which means that the solution u,~ would already be
defined on the triangle T and would not depend on the physical parameters associated with this
triangle, which is not admissible.
Now let P, = P(2) the degrees of freedom being the values of the flux at the vertices and at
the midpoints of the sides of the triangle. When one side of the triangle is lighted, we have to
write three equations, for example
(~~,~-~)~dxdy=O

forall

vEP(l).

/T

When two sides of the triangle

J7(Flu, -

are lighted. we have to write only one equation

F) dx dy = 0.

We will see further that continuous

freedom (nodal methods).

methods

on triangles

may be defined

with less degrees

of

Example 4: Consider the four nodes isoparametric

quadrilateral.
We may have one, two or three
lighted sides. In the last case the approximate
solution uh is already known on the quadrilateral

P. Lesaint / Finite element methods for the transport equation

361

and the physical parameters associated with this element have no influence
will see in the next section how to avoid such a situation.
Example 5: Consider a rectangular
uh E Q(k) by writing

JT

(Hu,-F)v
uh

dx dy=O

domain

forall

given on a - T f rom previous

D and rectangular

elements

on the solution.

We

T, we define the solution

u~Q(k--1),
calculations.

(3.3)

In the case of rectangular elements we show [lo] that the method is stable for a discrete L, norm
defined by approximating
the integral Jr.w2 dx dy by the quadrature
formula using the k*
points with coordinates
equal to the Gauss-Legendre
abscissa. If the exact solution is sufficiently smooth, the error measured by the norm defined above is of order hk+.
Remark 3.1: Assume D = [O,l] X [O,l], p > 0, v > 0, u = g on 8 - D, where g is continuous along
a - D and smooth on each side of a - D. Then the exact solution u is continuous on D, but the
first derivative of u has a jump when crossing the characteristic
line going through the origin.
The exact solution is smooth on the two parts of the domain D, and the order to approximation
is no longer valid.
Remark 3.2: The approximate
solution can also be defined
Gauss-Legendre
points, instead of writing eq. (2.3).
3.2. Discontinuous

by assuming

collocation

at the

methods

Let P, be a space of shape functions

defined by

on the element

T. The approximate

JT(HU~-F)V~X~~-__.(~~,+V~)(U~-~~)U~S=O,

forall

solution

UEP~,

uh E P, is

(3.4)

where n = (n,, ny) denotes the outer normal on aT, and &, is the value of the approximate
solution on a - T, known from previous calculations.
To calculate the solution u,, we must solve a linear system of equations on each triangle.
Assume that the space P, contains the constant function, then the conservation
of particles is
achieved in the following sense

/ a+T

(/HZ, + vnY)uh ds +

Let [u] = u, - u_ denote

JT

(JU/, dx dy =

the jump

J a-T

of a piecewise

l~~,+vn,l<,
defined

ds+
function

JT

Fdx

u, across

dy.

(3.5)

the face S of a

P. Lesaint / Finite element methods for the transport equation

362

triangle T, where u, and u_ are written,

value of u. We have the following stability

sa_f2w2

I un,y+

respectively,
for the down-stream
estimate [6]

vy,,I ds.

(3.6)

In the case when T is a triangle and P, = P(k) or when the element

quadrilateral,
we get the following error estimate [6]
)I u -

Assume
distance

Uh

11L2(a)

uh

II

T is a Q(k)

o(hk+12).

now that the element

between the midpoints

I] u -

and the upstream

isoparametric

(3.7)

T is a rectangle, or a slightly distorted

quadrilateral,
i.e. the
of the diagonals of T is of order h2, then we may write [13]
(3.8)

L2(f2)= @@+'b

Remark 3.3: For rectangular

meshes, both continuous
and discontinuous
methods can be run
and the order of approximation
is the same in both cases. However the discontinuous
method is
more stable, thanks to the jumps appearing along the sides of the elements. A good illustration of
this extra stability is found when applying both continuous
and discontinuous
methods for
solving the differential
equation du/dx
= uu = 0 on [0, h], u(O) = g. Using polynomials
with
degree I 1, we get

Uh(h)

1 - h/2
1 +X/25?

for the continuous

uh(h-)

(3.9)

method,

with h = ah, and

1 - h/3
1 + 2h/3

(3.10)

+ A2/6 g

for the discontinuous

method.
Eqs. (3.9) and (3.10) define approximations
of exp( -X) lying, respectively, on and below the
diagonal of the Pade table and it is known that the approximation
(3.10) is well suited for stiff
differential equations.

4. Continuous

and discontinuous

methods

for the one-dimensional

spherical

problem

Let R = IAr, 1 = JAp, r, = iAr, p, =jAp, K,, = [I;, ;+,I X [p,, P,+~]. The problem (2.19) is
solved by sweeping the rectangles following the caracteristic
direction as shown in fig. 5.

P. Lesaint / Finite element methods for the transport equation

Fig. 5. Nodal methods

4.1. Continuous

363

on rectangles.

method

The approximate
solution belongs to P(1) on each element and the continuity at the interface
of two elements is insured at only one point of this interface. The numerical scheme reads as
follows

(4.1)

for all the elements

K,,. The integrals

are approximated

by quadrature

formulas:

(4.2)

(4.3)

The weight wl+ 1,2 and the abscissa Y,+1,2 are chosen so that the quadrature
be exact for the
space of polynomials
spanned by r and r* [18]. The weight ~~+r,~ and the coordinate
~,+r,~
may be defined so that the quadrature formula be exact for the space of polynomials spanned by
p and 1 - 3~~ [12]. For other choices of the weights, we refer to refs. [18,12]. The quantities
P = 0, we discretize the reflexion
u,(& ~.,+r,~) are given, for -J <j I - 1. At the boundary
condition ~(0, p) = ~(0, -p), which gives

P. Lesuint / Finite element methous for the transport equution

364

To define the values of uh for p = - 1, we solve the differential

dq
-__
dr +ocp=F(r,

-1).

equation

O<r-=zR,

cp( R ) = given value.

in the following

way

dVh

I:, . ,

Ji

-dy+q,,-F(r.

-l))rdr=O,

O<i<I-

1,

(4.5)

r,

the integral

being computed

with the quadrature

We can show [12] that conservation

of particles
Let 11. 11y denote the weighted norm defined

formula

is achieved
by

(4.3) and we write

in a discrete

sense.

(4.7)
and consider the discrete norm
II . I/ y,h derived by computing
the integrals appearing
in
expression (4.7) with the quadrature
formulas (4.2) and (4.3). Let h = max( AI-, APL), we have the
asymptotic error estimates [12], for a smooth solution u
I/ u - u/, 112,3.h = 0( Iz~ ]log h I ,j2),
11u -

uh

11*.,, =

O( h6 llog h

II u

uh

II 0.h

0(h]log

1 /2).

hl2).

Remurk 4.1: The scheme (4.1) to (4.6) is the analog for the spherical geometry of the diamond
scheme for planar geometry and we could expect an error estimate of order h2. In fact since the
error
scheme is not really centered on the mesh K,,, one can show [12,18] that the truncation
with respect to the Y variable is of order Ar2/r, + 1,,2 on each element.
4.2. Discontinuous

method

The approximate
solution belongs to Q(l) on each element, the degrees of freedom being the
values of the solution at the corners of the element. Continuity at the interface of two elements is
not required. The numerical scheme reads as follows
+r$((l
Jhj+r2uh)

-p)u,,)

pr2nr + r(1 - /+z,)(

+or2u,T-r2Fjuh

Uh - &)Ch ds = 0,

drdy

(4.8)

365

P. Lesaint / Finite element methods for the transport equation

for all u,, E Q(l), for each rectangle K,,, where &, denote the value of the approximate
solution
outside K,,; the scalars n, and nP denote the components
of the outer normal
along a - I$,
with length equal to one on a - K,,. Let

Taking

v,, equal to one in relation

r ,:lmt,J2(ri+l)

(4.8), we get the following

I,h(ri) + /_:1/rz+r2(
ouh
- F)
rt2m
r,

The scheme (4.8) is stable for the norm

estimate

Remark

I] . I( 2 defined

conservation

property:

dr dp = 0.

by relation

(4.7) and we have the error

4.2:

More accuracy might be needed on the boundary r = 0. It is then possible to use a

discretization
with polynomials with degree 2 or 3 with respect to the r variable. The formulation
of the scheme remains the same. The cost of computation
is increased on each mesh but the
error is respectively of order Ap2 + Ar3 and Ap2 + Ar4 for the ]I - 112 norm.
5. Nodal methods for the (x, y) planar geometry
The origin of the nodal methods are to be found in ref. [19] and the bibliography
of this
article. In this section we shall describe this method for rectangular meshes as is done in ref. [19]
and for triangular meshes [ll]. With the same arguments like in ref. [ll], we can also derive a
method for regular hexagonal meshes.
5. I. Rectangular

meshes

Consider a rectangle AiA2A3A4,

with lighted faces A,A,
take the mean value of eq. (1.16) with respect to x. We get

and A,A,

(fig. 5). For any y, we

b
U(.G

Y)

dx

i
= +-lhFdx
a

- &(u(b,

4.~

Y>

dx

Ja
y) - ~(a,

Y>>.

This differential equation is then integrated exactly

value of eq. (2.17) with respect to y. We get

= &ldy

- &(u(x,

d) - u(x,

c)).

(5.1)
from c to d. For any x, we take the mean

(54

366

P. &saint

/ Finite element methods for the trunsport equation

This differential
equation is then integrated
exactly from a to b. We have got two relations
between the fluxes and their mean values on each side of the rectangle. Approximating
the fluxes
on each side by their mean values, we get the following scheme
+ 1 -e-p

u
T

+ l-e-

u
R

uT

1 -epDuL+
ff

_puB+

u R=e

-aUL

1 --;-^

1 -e-BF
u

(5.3)

us + 1 - eeaF,

(5.4)

where the subscripts T, R, B and L are respectively, written for top, right, bottom and left, and
where (Y= u Ax/p, p = u A Y/Y. The mean value of the source term F is denoted on each mesh
by F. The average ii of the flux on each mesh is then defined by the balance equation

uu=F-~(u,-u,~)-~(uT-u,).

(5.5)

AY

For more sophisticated

schemes involving higher order moments we refer to ref. [19]. This
method is shown to be stable for small values of a and p [ll]. Denoting by ii the mean value of
the exact solution u on each rectangle T belonging to the domain D, we get the following error
estimate
[Ax Ay c

(ii - ii)2j

= 0(h).

(5.6)

TcD

Numerical experiments
(5.6) is not optimal.

[19] show an error of order h2, which means that the theoretical

estimate

5.2. Triangular meshes

Let T be the triangle A,A,A,

Y =

t&i -

V)Y,

with two lighted faces CA and CB. We define the mapping

M-r

+ t(i + 77)Y, + (1 - 5)_v3.

The triangle T is the image of the reference triangle ? by the mapping MT, as shown in fig. 6.
transformaTo any function v defined on T, we associate i? = 110 M, . Using this coordinates
tion, equation (1.16) may be written as
a-

aa
a<

+b-

aa
317

-+uJG=J@,

a = iA,A,(p~n,

where

+ in_,.),, > 0,

(5.7)
a + b = -A,A,(~H.,

+ YPZ,.),~> 0,

P. &saint

367

/ Finite element methods for the transport equation

Fig. 6. Nodal method

on a triangle with two lighted sides.

and
J=-

area T
area T

= area T.

For a fixed value of <, we integrate

eq. (5.7) from TJ= -5

to 17= +[

and we use the identity

and we get

Approximating
the fluxes a( 5, - 5) and i2(& 5) by their already known mean values on the
eq. (5.8) from [ = 0 to 6 = 1
lighted faces u3u2 and u3u1, we integrate exactly the differential
and we get an expression of the mean value of the flux on the face u1u2.
Now let T be a triangle A,A,A,
with only one lighted face A,A,. We define the mapping M,
by
x = Ex, + 77x3 + (I-

C - rl)x,,

Y = 5Y, + 77Y3 + (1 - 5 - 77)Y,.

The triangle

T is the image of the reference

triangle

? by the mapping

M,

as shown in fig. 7.

368

P. Lesaint / Finite element methodsfor the transport equatron

Using this coordinates

transformation

aa

ai-i

CY~ +,8zi

we get the equation

+aJC=JF,

where

a = -A,A,(w,

+ 1.>13
<0

,l3= -A,A&n,

+ YH.& < 0.

(5.9)

For a fixed value of E, 0 < 6 < 1, we integrate

and we use the identity

eq. (5.9) with respect

to 77 from 17= 0 to n = 1 - .$

We integrate exactly this differential

equation from < = 0 to 5 = 1. Appproximating
the fluxes
along the faces b,b, and h,b, by their mean values, we get a first relation between the mean
values of the flux on the three faces of the triangle. A second relation may be derived by
integrating a similar differential
equation with respect to 77. The method is shown to be stable
[ll] and the theoretical
rate of convergence
is O(h). Remark 5.1: Another approach
to the
problem is done in ref. [15], where a nodal method is used in the case of two lighted faces and a
characteristic
method in the case of one lighted face.

Fig. 7. Nodal method

on a triangle with one lighted side

P. Lmaint / Finite element methods for the transport equation

369

References
[l] R.E. Alcouffe and E.W. Larsen, A review of characteristic
methods used to solve the linear transport equation.
Proc. ANS-ENS
Intern. Topical Meeting on advances in mathematical
methods for the solution of nuclear
engineering problems, Munich, Fed. Rep. Germany (27-29 April 1981).
[2] G.I. Bell and S. Glasstone, Nuclear reactor theory (Van Nostrand Reinhold, New York, 1970).
[3] B.G. Carlson and K.D. Lathrop, Transport theory - the method of discrete ordinates. Computing
methods in
reactor physics, eds. H. Greenspan,
Kelber C.N. and D. Okrent (Gordon and Breach Science Publishers, New
York, 1968).
[4] P.G. Ciarlet, The finite element method for elliptic problems (North-Holland,
Amsterdam,
1978).
[5] J.J. Duderstadt
and W.R. Martin, Transport theory, (John Wiley, New York, 1979).
[6] C. Johnson, Finite element methods for convection
diffusion problems.
Fifth Intern. Symp. on Computing
Methods in Engineering and Applied Sciences, INRIA, Versailles (dec. 1981).
[7] E.W. Larsen and R.E. Alcouffe, The linear characteristic
method for spatially discretizing the discrete ordinates
equations in (x, y) geometry. Proc. ANS-ENS Intern. Topical Meeting on Advances in Mathematical
Methods
for the Solution of Nuclear Engineering Problems, Munich, Fed. Rep. Germany (27-29 April 1981,.
[S] P. Lascaux, Equation
du transport.
Analyse Mathematique
et Calcul NumCrique
pour les Sciences et les
Techniques, eds. R. Dautray and J.L. Lions, tome 3 (Masson, Paris, 1985) pp. 982-1016.
[9] K.D. Lathrop, Spatial Differencing of the Transport Equation: Positivity vs. Accuracy. J. Comput. Phys. 4 (1969)
475.
[lo] P. Lesaint, Continuous and Discontinuous
Finite Element Methods for Solving the Neutron Transport Equation.
MAFELAP 1975 Brunel University, ed. J. Whiteman (Academic Press, New York, 1976).
(111 P. Lesaint, Nodal Methods for the Transport Equation, MAFELAP 1984 ed. J. Whiteman (Academic Press, New
York, 1985).
[12] P. Lesaint, Approximation
de Problemes de Transport en Neutronique.Computing
methods in applied science,
second Intern Symp, dec. 1975 IRIA. Lecture Note in Economics and Mathematical
Systems, No 134 (Springer,
Berlin, 1976).
[13] P. Lesaint and P.A. Raviart, On a Finite Element Method for Solving the Neutron
Transport
Equation.
Mathematical
aspects of finite elements in partial differential equations, ed. C. de Boor (Academic Press, New
York, 1974) pp. 89-123.
[14] P. Lesaint and J. Germ-Roze, Isoparametric
Finite Element Methods for the Neutron Transport Equation. Intern.
J. Numer. Meth. Eng. 10 (1976) 171.
[15] R. Paternoster,
A Linear Characteristic
Nodal Transport
Method for the Two-Dimensional
(x, y) Geometry
Multigroup Discrete Ordinates Equations over an Arbitrary Triangle Mesh.
[16] W.H. Reed, Triangular Mesh Difference Schemes for the Transport Equation. LA-4769, Los Alamos Scientific
Laboratory (1971).
[17] W.H. Reed and T.R. Hill, Triangular Mesh Methods for the Neutron Transport
Equation. Proc. 1973 conf. on
Mathematical
Models and Computational
Techniques for Analysis of Nuclear Systems, Ann Arbor, Michigan
(9-11 april 1973).
[18] W.H. Reed and K.D. Lathrop, Truncation Error Analysis of Finite Difference Approximations
to The Transport
Equation. Nucl. SC. Eng. 41 (1970) 237.
[19] W.F. Walters and R.D. ODell, Nodal Methods for the Discrete Ordinates
Transport
Problems in (x, y)
geometry. Proc. ANS-ENS Intern. Topical Meeting on Advances in Mathematical
Methods for the Solution of
Nuclear Engineering Problems, Munich, Fed. Rep. Germany (27-29 Aril 1981).