Finite element method for thermal analysis of concentrating solar receivers
Stanko Shtrakov and Anton Stoilov
South-West University, Blagoevgrad, Bulgaria
Application of finite element method and heat conductivity transfer model for calculation of
temperature distribution in receiver for dish-Stirling concentrating solar system is described. The
method yields discretized equations that are entirely local to the elements and provides complete
geometric flexibility. A computer program solving the finite element method problem is created
and great number of numerical experiments is carried out. Illustrative numerical results are given
for an array of triangular elements in receiver for dish-Stirling system.

1. Introduction.
Cavity receivers for solar concentrating systems absorb a concentrated solar energy, convert it to
heat, and transfer heat to the working gas in power unit (Stirling engine, steam turbine). Heat flux in
such elements is great and material tension is very high. Different shapes and materials were used to
optimize the performance of cavity receivers.
Besides the natural experimental tests, there are possibilities to use serious theoretical studies for
the thermal behavior of the cavity receivers. It would be very useful, because the processes in the
cavity are complicated by high energy intensity and a real measure of temperature distribution is
difficult and expensive.
Because of the special form of the solar receivers, it is not suitable to use ordinary mathematical
technique, such as analytic methods or numeric solutions with finite difference approximation. Finite
elements method is now one of the most popular ways to solve complicated mathematical problems,
especially with irregular definition area. There are many works treating the heat conductivity problems
by using finite elements methods and there is practical experience in this field.
In this paper is presented a mathematical model for heat conductivity processes in cavity receiver
for dish-Stirling system with appliance of finite elements method for solving differential equations.
The studies and numerical tests are made for a typical cylindrical cavity receiver, but algorithm and
computer program created for calculation of temperature distribution in receiver are suited to solving a
different forms of receivers.
The simple construction scheme of typical
Absorber
tube type cavity receiver for dish-Stirling system
is shown in Fig. 1. It comprises absorber plate
(tube system for gas heating), cavity space, walls,
Cavity
aperture and insulation on the external surfaces of
walls.
The absorbing surface is usually placed
Walls
behind the focal point of the concentrator so that
the flux density on the absorbing surface is
reduced. The size of the absorber and cavity
Aperture
walls is typically kept to a minimum to reduce
Insulation
Fig. 1. Cavity receiver for dish-Stirling system heat loss and receiver cost.
Concentrated radiation entering the receiver
aperture diffuses inside the cavity. Most of the energy is directly absorbed by the absorber, and most of
the remainder is reflected or reradiated within the cavity and is eventually absorbed by absorber and
cavity walls. Converted by walls heat is conducted to the absorber plate and is transferred to the
working gas. The major advantage of cavity receiver is that the size of the absorber may be different
from the size of aperture. With a cavity receiver, the concentrator’s focus is placed at the cavity
aperture and the highly concentrated flux spread inside the cavity before encountering the larger

absorbing surface area. This spreading reduces the flux incident on the absorber surface. When incident
flux on the absorbing surface is high, it is difficult to transfer heat through surface without thermally
overstressing materials.
The bottom of the receiver is a heat exchanger for the Stirling engine. Part of solar energy strikes
directly the bottom wall of the receiver and by conduction it delivers heat to exchange pipes.
Cylindrical walls absorb the other part of solar energy and conduct it to the bottom of the receiver. The
external cylindrical surfaces are well insulated for protecting heat losses. Natural convection in cavity
produces air circulation and convective losses. Heat losses are further caused by insensitive radiation
from cavity aperture to the ambient. These losses are not big, because of the small measurements of the
receiver. Nevertheless, all heat losses can be taking into account because they influence the
temperature distribution in the receiver.
2. Theoretical Model
Mathematical model of processes in cavity receiver can be derived on the base of theoretical
treating of heat conductivity in receiver walls, heat transfer to the working gas by absorber plate and
heat losses to the ambient. Because of symmetry, the problem domain can be simply presented as in
fig.2. The main mathematical equation in this problem is the heat conductivity equation for receiver
construction written in cylindrical coordinate system:

∂ ⎛ ∂T ⎞ ∂ ⎛ ∂T ⎞ 1 ∂T
= ∆2T = 0
⎟+ λ
⎜λ
⎟ + ⎜λ
∂r ⎝ ∂r ⎠ ∂z ⎝ ∂z ⎠ r ∂r

(1)

where T is the temperature in the walls, r and z – space variables in radial and high direction (fig.2)
and λ – heat conductivity coefficient [W/m K]. ∇ 2 is Poason’s operator for vector form description of
equation.

z

E
C
D

The boundary conditions, needed for the solution
of eqn (1), can be obtained from consideration of energy
balance for different surfaces of the cavity receiver.
External surfaces of the receiver (D) are insulated and the
heat transfer to the ambient can be presented by the next
equation:

λ

A

∂Τ
= K D (Τ − Τa )
∂r

(2)

Here Ta is the ambient temperature, T – temperature
of external surface and Kd – overall heat transfer
coefficient from the receiver external surface to the
ambient air (including insulation).
r
The inner surfaces of receiver (A and C)
B
exchange heat energy with air in cavity and absorb solar
Fig.2 Theoretical scheme of receiver radiation entering the space in receiver. Boundary
condition for these surfaces can be written as:

λ

(

)

∂Τ
= α B Τ − Τf + q
∂r

(3)

where αB is a convective transfer coefficient in receiver cavity [W/m2 K], Tf – air temperature in cavity,
T – temperature of inner receiver surface and q – solar energy flux [W/m2].
Surface E can be modeled as a heat exchanger for heating the working gas in Stirling motor.
Boundary condition in this case is:

λ

∂Τ
= K W (Τ − Τwg )
∂r

(4)

Heat transfer coefficient Kw depends on parameters of convective heat transfer to the working fluid.
Temperature of working fluid Twg is defined from the thermodynamic treating of Stirling process. In
general, boundary conditions can be written in common form as follow:

λ

∂T
= h eT + C e ,
∂n

(5)

where he refers to the convection coefficient and Ce describes the convection coefficient with the
external temperature and solar flux. Superscription e refers to the surface index (e). Derivative ∂T/∂n is
taken to the outward normal direction n for every boundary surface.
Differential problem (1) – (4) can be solved by numerical methods. The special form of problem
domain (fig.2) leads to creating an unstructured grid for descretization of calculation area. This is the
reason for using the finite element method for solving the problem.
3. Finite Element Method
In the finite element method (FEM), the problem domain is discretized and represented by an
assembly of finite elements. The method yields discretized equations that are entirely local to the
element. As a result, the discrete equations are developed in isolation and are independent of the mesh
configuration. In this way, finite elements can readily accommodate unstructured and complex grids,
unlike other numerical methods (finite difference method) requiring structured (rows and columns)
format. Finite element method provides complete geometric flexibility.
Various approaches to the finite element method formulation are used, but the most prevalent is the
Galerkin method [1,2]. It provides the correct number of basis function and the same resulting
equations as methods based on other (more complicated) types of formulations.
Galerkin’s method selects the weight functions equal to the basis functions (shape functions) of the
approximate solution. It will be demonstrated on the example of solving the two-dimensional thermal
conductivity problem in cavity receiver for dish-Stirling system – equation (1) with boundary condition
(2) – (4).
Commonly encountered elements in two-dimensional configuration include a linear triangle,
bilinear rectangle, and bilinear quadrilateral. Triangular elements are well suited to irregular
boundaries, and spatial and scalar interpretation can be accommodated in terms of a linear polynomial.
Nodes are assigned to location in the element at
local node 1
which unknown function, such as temperature, is to be
(r1,z1)
determined. Nodes are often placed only at the corners
of the elements (Fig.3), but additional nodes can be
placed internally or along the element boundaries.
The main problem in finite elements method is to
local node 3
choose shape or integration function to provide
(r3,z3)
approximate variation of the depended variable within
local node 2
each element between the values of the nodes. The
(r2,z2)
simplest shape function is linear. Each shape function
is a local interpretation function that is defined only
Fig. 3. Linear triangle
within elements containing a particular node. Linear
triangle element (Fig.3) refers to linear interpolation along a side or within the element. For an
interpolation involving a scalar shape form function, N(r,z) within a linear triangle can be described as:
N(r,z)=a+br+cz.
(6)
where the unknown coefficients, a, b, c, can be determined on the basis of the substitution of nodal
values. For example, at node 1, the position is (r1, z1) and the scalar shape function is N1.
By using the shape functions, an approximate solution Tˆ ( r , z ) for T(r,z) is assumed in the
form

N

Tˆ (r , z ) = ∑ T j N j (r , z )

(7)

j =1

where Tj are the temperatures at the nodes desired for the solution, and Nj are the shape functions that
are equal to 1 at each node. When this approximation is substituted into the energy equation (1), there
is a residual that depends on r and z:

∂ ⎛ ∂Tˆ ⎞ ∂ ⎛ ∂Tˆ ⎞ 1 ∂Tˆ
⎟+ λ
⎟ + ⎜λ
⎜λ
= ∆2Tˆ = Re s (r , z )
⎟
⎜
⎟
⎜
∂r ⎝ ∂r ⎠ ∂z ⎝ ∂z ⎠ r ∂r

(8)

It is desired to be obtained a solution that, in average sense over the entire volume, is as close
as possible to the exact solution each time. Variational principles are applied to minimize the residual.
A set Wi(r,z) of independent weighting function is applied, and the residual is made orthogonal with
respect to each of the weighting functions. This provides the following integral, which is evaluated for
each of the sets of independent weighting functions,

∫∫ Re s(r , z )W (r , z )dS

(9)

i

S

The integration is over the whole volume in which the solution is being obtained. In the
Galerkin method the weighting functions are chosen to be the same function set as the shape functions.
Equation (9) provides the Galerkin form of energy equation for each of the weighting functions,

⎡ N
⎤
T j ∇ 2 N j (r , z )⎥ ⋅ N i (r , z ) ⋅ dS
∫∫S ⎢⎣λ ⋅ ∑
j =1
⎦

(10)

Boundary conditions can be incorporated through the boundary integral term in above
equation. Considering the general form of the boundary conditions (5) and using integration by parts,
the following form of equation (10) can be received [1]:

⎡ N ∂N i ∂N j ∂N i ∂N j 1 ∂N j
⎤
+
+
(
N i )T j ⎥ ⋅ dS = ∫ [ N i (r , z )(h eTi + C e )dn
∫∫S ⎢⎣λ ∑
∂r ∂r
∂z ∂z
r ∂r
j =1
n
⎦

(11)

Evaluating equation (11) for each
i provides N simultaneous algebraic
equations for the Tˆ ( r , z ) . In the matrix
form, this is:
[Kij] [Tj] = [Fj]
(12)
Since each shape function Nj(r,z)
(and hence each weighting function in
Galerkin method) is zero except within an
element containing Tj, the resulting matrix
7
8
9
10
11
12
[Kij] for solving the simultaneous
6
4
10
2
8
equations for Tj is banded and sparse. It is
3
1
9
5
7
usually
banded along the diagonal (banded
3
4
5
6
1
2
matrix).
The finite element method yields
Fig. 4. Discretization of domain
discretized equations that are entirely local
to the element and hence the global matrix [Kij] is a simple combination (sum) from local matrixes of
each node.
Solving the system of algebraic equations (12), the unknown temperatures Tj can be received
as a result. Because of the special form of [Kij] matrix, different numerical techniques have been
developed [1,2] for solving the system (12).

4. Numerical solution techniques
The mentioned above technique is demonstrated on the example of cavity receiver (fig.1).
Using triangle elements the problem domain can be discretized as it is shown in fig.4. Each element is
numbered and nodes coordinates (ri,zi) are specified.
Interpolation of the temperature in triangle element can be written in terms of the shape
functions, Ni, Nj, Nk, as follow:
T (e ) (r , z ) = Ti N i(e ) (r , z ) + T j N (je ) (r , z ) + Tk N k(e ) (r , z ) where function N ie = a ie + bie r + cie z ( i =1, 2, 3
) is determined by considering the next conditions:

N ie (ri , z i ) = 1

N ie (r j , z j ) = N ie (rk , z k ) = 0 .

Superscription e refers to the element number (e). Coefficients a, b, and c in shape form functions have
the following form:

⎧ai = (r j z k − rk z j ) / 2 S
⎪
⎨ bi = (z j − z k ) / 2 S
⎪ c = (r − r ) / 2 S
k
j
⎩ i

⎧a j = (rk z i − rk z i ) / 2 S
⎪
⎨b j = ( z k − z i ) / 2S
⎪
⎩c j = (ri − z k ) / 2S

⎧a k = (ri z j − r j z i ) / 2S
⎪
(13)
⎨bk = (z i − z j ) / 2S
⎪
⎩c k = (r j − ri ) / 2S

where S is the area of triangle element.
Finite element method is well-established numerical technique for heat conduction problem in
orthogonal coordinate system. Heat conductive problem in cylindrical coordinate system (equation (1))
differs from the orthogonal coordinate problem with next integral in equation (11):

λ ∫∫
S

1 ∂N j
N i drd
r ∂r

(14)

The shape function derivatives in equation (11) are constants: ∂Ni/∂r = bi and ∂Ni/∂z = bi. This
gives a simple form of the local matrix for element n (part of global matrix [Kij]) in orthogonal
coordinate system (without contribution of integral (14)). With regarding the above values of
derivatives and not considering the boundary conditions the matrix can be written in the following
form [1]:

k (n)

⎡ (bi n ) 2 + (ci n ) 2
⎢ n n
n
n
= λ ⎢bi b j + ci c j
⎢b n b n + c n c n
i
k
⎣ i k

bi b j + ci c j
n
n
(b j ) 2 + (c j ) 2
n
n
n
n
b j bk + c j c k
n

n

n

n

n
n
n
n
bi bk + ci c k ⎤
n
n
n
n⎥
b j bk + c j c k ⎥
n
n
(bk ) 2 + (c k ) 2 ⎥⎦

(15)

Integral (14) is more complicated because there is 1/r member and Ni is in form (6) with
coordinates r and z. The integral (14) can be divided in three parts:

I = λ ⋅ b j ∫∫

Ni
drdz = λ ⋅ b j [
r

ai

∫∫ r drdz + ∫∫ b drdz + ∫∫ c
i

i

z
drdz ] = λ ⋅ b j [ I1 + I 2 + I 3 ]
r

ntegral
I2 is easy to be calculated because bi is constant and: I2 = bi S, where S is the area of triangle element.
Other integrals can be calculated by integrating over triangle area (fig. 5). By using the triangle
elements with particular disposition as it is shown in fig.4 and dimensions presented on the fig.5 the
integrals I1 and I3 can be presented as:

I1 =

α i r2
[1 +ν (lnν − 1)]
µ
2

r
3 2 r
I 3 = [ r1 ln ν ( z 2 + 1 ) − z 2 ∇ r ( r1r2 − r1 − 1 )]
µ
4
4
2µ
ci

(17)

I

µ=

where

r
∇r
и ν = 2
r1
∇z
z

3

3

3

2

r2= r3

z
1

2

r

1

r2

r1= r3

r

2
r1
(б )

(a)

r2
1

r1= r3
(в )

Fig.5 Triangle elements
These values must be added to the corresponding members in matrix [Kij].
The cylindrical heat conductivity problem can be solved by approximate methods. It is necessary
in case of common disposition and form of triangle elements. One of these methods is to make
integration in (14) assuming coordinates r and z as constants and equal to the coordinates of the mass
center of the element - rm and zm. In this case approximate solution can be received, and it will be so
close to the real solution as much as small are elements. Integral (16) can be presented as:

⎡a
z ⎤
I = b j S e ⎢ i + bi + ci m ⎥
rm ⎦
⎣ rm

(18)

Other approximate method is determined by rearranging the original equation (1) in form:

∂ ⎛ ∂T ⎞ ∂ ⎛ ∂T ⎞
⎜ λr
⎟ + ⎜ λr
⎟=0
∂r ⎝ ∂r ⎠ ∂z ⎝ ∂z ⎠

(19)

This form resembles to the form of conductivity equation written in Decart’s coordinates if
product λr is interpreted as modified conductivity coefficient. It will vary for different elements and will
depend on radius r. For integration in (11) it can be assumed approximate approach, for example to use
constant radius rm of the mass center of element. In this case the matrix (15) can be used, but the
coefficient λ will be different for each element of the matrix.
Boundary conditions (5) contribute the matrix elements by adding values for triangle sides that
touch the boundaries of the receiver: surfaces A, B, C, D, E (fig.2). Contribution is determined by
integrating the last members in (11). It differs for elements with equal indexes and with different
indexes:
gii = gjj = gkk = he∆r/3 ; gij = gjk = gik = he∆r/6 ;
(20)
Values in right column of matrix equation are determined for triangle sides that contact the
ambient or exchange elements of receiver. These are formed by integrating the members not containing
temperature in last part of (11):

Fi =

1 e
1
⋅ C ⋅ ∇r or Fi = ⋅ C e ⋅ ∇z
2
2

(21)

5. Programming and result
A computer program solving the finite element problem has been created and great number of
numerical experiments has been carried out. The program comprises two main parts. The first part
includes program modules for automatic discretization and numbering the elements and nodes in mesh
configuration. These modules are specific for different problems and configuration of domain. The

85
4

85 0

80 0

73 5

80 0

85 0

73 5

73
2

2
73

4
85

main result of performance of program modules is forming the matrix with numbers and coordinates of
triangle element nodes.
The second part of the program modules uses created in the first module data to form the basic
matrix Kij – equation (12). Special procedure for solving the band matrix is used to calculate the
requested parameters (temperatures, heat flux etc.). This part is standard and can be used for solving
different tasks.
The program and mathematical treating of finite element method were verified by using the
exact analytical solution of thermal conductivity problem. It is known, that for simple configuration of
domain, the heat conductivity problem has analytical solution. Such configuration is a cylindrical
domain. For such a simple configuration the finite element method gives results, which a very close to
the analytical solution (difference is below 0.1%).
Fig. 7 presents the
73
0
0
73
temperature distribution in
receiver for dish-Stirling
85
8
8
85
system,
calculated
by
computer program. Some of
the used data for boundary
conditions
are
assessed
710
710
approximately. This means
that
the
temperature
distribution in receiver must
be interpreted only in rough
figures. When the conditions
for receiver performance are
700
700
defined more precise, the
results
for
temperature
distribution will be more
accurately.
Fig.7 Temperature distribution in Receiver
6. Conclusion
Finite element method has been applied for heat transfer problem in receiver for dish-Stirling
system. It uses unstructured and complex grids, which are suited to the special forms of absorber
element. Finite element method provides complete geometric flexibility.
Presented mathematical model and algorithm for program is universal and can be used for
different tasks. In this paper are discussed only principle aspects of using the finite element method for
heat transfer processes in receiver for solar receivers in concentrating solar systems. Numerical
experiments carried out in this research show that this technique gives very good approaches to the real
thermal processes. It can be used successfully to investigate different conditions and forms of receivers
for concentrating solar systems and other thermal equipments.
References:
1. Zienkiewich O.C., K. Morgan, Finite Elements and Approximation, John Wiley & Sons,
New York, 1983.
2. Zienkiewich O.C., The Finite Element Method in Engineering Science, McGraw-Hill,
London, 1971.
3. Naterer G.F, Heat Transfer in Single and Multiphase Systems, Taylor & Francis, New York,
1997
4. Siegal R, J.R.Howell, Thermal Radiation Heat Transfer, Taylor & Francis, New York, 1998

5. Strang G, G.J. Fox, An Analysis of the Finite Element Method, Prentice-Hall,Inc,
Engeleword Cliffs, 1973.
6. Norrie D.H., G de Vires, An Introduction to Finite Element Analysis, Calgary, Alberta,
Canada.
7. Kreith F., W.B. Black, Basic Heat Transfer, Harper and Row, Publishers, New York, 1980.

