Despite recent advances in computer aided design (CAD) based tools,
spacecraft thermal analysis remains outside the realm of finite
element method (FEM) based analysis. The primary complaints against
FEM often cited are:
is not based on physical principles.
codes do not provide procedural modeling for heaters, heat
pipes, or other abstract thermal control components.
codes generate inappropriately large thermal models.
However, a failure on the part of existing FEM based codes
does not invalidate the advantages of the Finite Element Method.
Properly implemented, FEM based systems can have significant advantages.
A simple first law interpretation of FEM is presented, and shows
that finite difference (FD) and FEM meshes may co-exist in the same
thermal model, and solved using traditional analyzers such as SINDA/FLUINT.
A description of an integrated FD/FEM based system that efficiently
satisfies all areas of spacecraft thermal analysis, including thermal
radiation, is also presented.
Despite the advantages that FEM has in integrating with CAD generated
design geometry and analysis models from other analysis disciplines,
its use for thermal analysis has seen very little application in
the aerospace community. The primary reason is that existing FEM
codes do not currently satisfy some key needs of the aerospace thermal
analyst. Many are based on structural codes and have solvers that
do not allow procedural modeling of thermal control components such
as heaters, and have weak or no support for thermal radiation.
The situation is understandable in that the thermal aerospace market
is small compared to most FEM code suppliers major customer base.
Structural analysis tools are applicable for most engineering markets.
In addition, there is virtually no demand from the thermal aerospace
community for FEM based tools.
The lack of demand is due in part to the fact that existing tools
do not meet all aerospace analysis requirements, but also due to
many misconceptions about FEM. This paper has two objectives. One
is to present a new view of FEM, based on familiar first law principles.
Many of the complaints against FEM, such as negative conductors
are explained when viewed from a conservation of energy principle.
When viewed in these terms, FEM can be seen to be fully compatible
with existing analyzers and modeling approaches.
The second objective is to present a description of a FEM system
that does satisfy the requirements of the aerospace thermal analyst.
In particular, the problems associated with radiation analysis are
One reason that FEM is not used more widely in the aerospace industry
is the belief that the method is purely an abstract mathematical
approach, rather than a first-law approach. The fact that it produces
"negative" conductors reinforces the belief that the method is physically
Part of the beauty of FEM is that it can be applied to a wide variety
of problems including stress, dynamics, fluid flow, electromagnetism,
and heat transfer. The generality of the method is highlighted in
most texts on FEM and the method is usually presented as an abstract
technique for solving any governing differential equation.
A weighted residual approach whereby the error in a trial solution
is minimized by the application of successive integration of orthogonal
weighting functions is often used to reduce the governing differential
equation into a set of algebraic equations that can be solved on
The mathematics involved in the derivation are typically beyond
what is found in the undergraduate engineering curriculum. Such
a difficult abstract approach does not inspire the confidence in
the method that one can easily obtain from a simple energy balance
approach. However, FEM does have a simple first-law interpretation
as applied to thermal problems.
POINT SOURCES - Consider the four nodes shown in Figure 1 along
with a point source (for instance a laser beam) of heat being applied
exactly at the corner of the nodes. From the symmetry of the problem,
it is expected that all four nodes will be at the same temperature.
A reasonable approach is to divide the heat from the point source
equally among the four nodes so that the model will also predict
1: Point heat source being deposited at the corner of four nodes
If the point source is displaced an infinitesimally small distance
from the corner and into the interior of one of the nodes, then
the usual approach is to lump all of the energy into that node.
Because of the infinitesimally small displacement away from the
corner, the temperatures of the nodes in reality will only change
an infinitesimal amount. However, predicted temperatures will show
a much greater change between the two nearly identical cases, since
the way in which energy has been distributed to the nodes has changed
If the predicted temperatures change dramatically from dividing
the energy equally among the nodes, compared to lumping it all into
one node, then the problem needs a finer nodal resolution. However,
no matter how fine the model is nodalized, the node that receives
all of the energy will still have a predicted temperature higher
than in reality.
Rather than increasing the amount of nodes in the thermal model
to obtain better accuracy, a different scheme for distributing the
energy may be employed. When the point source is at the corner of
the four nodes, the most accurate results are obtained by equally
dividing the energy among the four nodes. On the other hand, when
the point source is at the exact center of a control volume, the
most accurate results are obtained when all of the energy is lumped
into that node. The transition between these two extreme conditions
can be handled in a smooth manner, rather than in the traditional
step change fashion. As the point source moves closer to a control
volume center, more energy is apportioned to that node, and less
to the neighboring nodes.
The apportionment of energy can be handled by a simple bilinear
interpolation scheme, carried out over an imaginary rectangle connecting
the centers of the four nodes (see Figure 2). The energy is divided
such that the sum of the energy distributed to all of the nodes
equals the total energy being deposited.
Figure 2: Interpolation region used to apportion
energy between neighboring nodes
The position inside the interpolation region can be denoted by
two scalar variables, u and v, that vary between zero
and unity. For example, the center of node 3 has u,v coordinates
(0,0) and the center of node 2 has coordinates (1,1). The amount
of energy each node receives by a point source of energy inside
the interpolation region is given by multiplying the total energy
of the point source by the factors listed in EQ (1).
The sum of all the factors is unity for all u,v positions,
ensuring that the deposited energy is conserved. The distribution
of energy into a node based on its proximity to the control volume
center can be viewed as the fundamental technique of FEM for heat
transfer. In fact, all heat sources are treated in this consistent
DISTRIBUTED SOURCES - Suppose that the nodes in Figure 2 are being
illuminated by a distributed, rather than a point source. A distributed
source is handled in the same manner by dividing it up into many
small areas, each of which is considered as a point source. For
example, consider a flux incident on the interpolation region that
can be expressed as a function of the u,v coordinates. The
amount of heat that is assigned to node 2 is given by:
This equation is the integral over the interpolating region of
the heat flux per unit area multiplied by the apportioning function
for node 2 given in EQ (1) (assuming a rectangular interpolation
The traditional method of apportioning energy into a control volume
may be viewed as a special case where the apportioning function
is a step function that is unity over the control volume, and zero
everywhere else. Note also that for a heat load that is constant
over the interpolation region, the results obtained using the apportioning
functions in EQ (1) gives the same results as a using the traditional
CAPACITANCE - The thermal capacitance of a node can also be derived
by using a consistent apportioning of energy viewpoint. The differential
form of the basic heat conduction equation relates the rate of internal
energy being stored at a point to the net rate of energy being deposited
by conduction and volumetric heat sources:
It can be seen that each of the terms has the same units - heat
flux per unit volume. The differential form of the heat conduction
equation treats all terms as point sources. We can use the smoothly
varying energy apportioning function to determine the contribution
due to stored energy to a node's energy balance equation.
For example, using the interpolation region shown in Figure 2,
the contribution to node 2's energy balance equation due to sensible
heating or cooling is given by:
The term enclosed in the square brackets is the thermal capacitance
associated with node 2 that is contributed by the material spanned
by the interpolation region. If the density and specific heat are
constant over the interpolation region, the integral assigns one
fourth of the thermal mass to node 2, the same as using step function
apportioning. The total thermal mass of node 2 is determined by
the sum of the contributions for all interpolation regions that
contain node 2.
CONDUCTION - The net energy accumulating at a point due to conduction
can also be apportioned to each node. Using the example in Figure
2 again, the net energy into node 2 by conduction is given by:
To simplify the discussion, let us use the symbol Ni
for the apportioning function for node i. We can also use Green's
theorem  to simplify the Laplacian in EQ (5) into a gradient.
The energy into a node by conduction is therefore also given by
Note that the gradient of the apportioning function,Ni,
points towards the node center. The apportioning function can
be visualized as a "tent" with a value of unity at the node center,
and a value of zero at the far boundaries of the interpolating region.
The gradient gives the direction of steepest ascent. When negated
and dotted with the temperature gradient, it gives the heat flowing
per unit volume towards the node by conduction (when multiplied
The contributions to the node's energy balance equation for heat
sources and thermal capacitance can be solved for directly. However,
the equation for the conduction term contains the temperature gradient.
We must therefore supply an approximation to the temperature gradient
as a function of the unknown temperatures so that the integration
may be carried out.
We can approximate the temperature gradient over the interpolation
interval by first approximating the temperature:
One can see that by substituting the apportioning functions given
by EQ (1) into EQ (7) that the above formulation is a simple bilinear
interpolation of the temperatures at the node centers over the interpolation
region. The gradient of the temperature is now given by:
Note that the same functions are used to interpolate the temperatures
across the region as are used to apportion energy into a node. In
FEM parlance, using the same apportioning functions (or "weighting"
functions) as the interpolating functions (or "shape" functions)
is called the Galerkin method.
The Galerkin method has a special significance in that it guarantees
a symmetric set of terms in the set of nodal energy balance equations.
Substituting EQ (8) into EQ (6) gives:
The net heat flowing into a node by conduction is given by an algebraic
equation of the four unknown temperatures. Furthermore, because
of symmetry, Gij is the same quantity as Gji.
Because the interpolating functions sum to unity at every point
in the interpolating region, the following relation is also true:
The terms to node i's energy balance equation due to conduction
can therefore be expressed in the familiar form:
NEGATIVE CONDUCTORS - One of the hardest conceptual aspects of
FEM to overcome for thermal engineers that have become accustomed
to a lumped parameter electrical analogy view of thermal modeling
is that for some types of interpolation regions, namely triangles
with interior angles greater than 90 degrees, negative conductors
can result. This is often the basis for eliminating FEM as an acceptable
modeling approach, since these negative terms are regarded as physically
It must be recognized that a conduction term between two nodes
produced by FEM does not represent the heat transfer between
the two nodes. The conduction terms in EQ (11) considered together
represent the total energy balance for a thermal node. The net heat
into a node is given by a linear combination of temperatures, which
can be rearranged into a form that is compatible with SINDA/FLUINT
. This form looks like "conductors", but they are not.
Let us consider a simple example to see the difference between
the global energy balance for a node and individual heat flows between
nodes. The traditional centroid approach uses an algebraic approximation
to the integral form of the heat conduction equation:
The above equation is put into algebraic form by using summations
as approximations for the integrals, and differences as approximations
for the differential terms:
The summation on the left is usually carried out with only one
integration interval, the entire control volume. The mean value
theorem tells us that the best place to evaluate the integral, in
lieu of any other knowledge of the solution, is at the centroid
of the volumetric element. The centroid method could also be called
a "finite-sum-finite-difference" method.
For rectangular nodes, the summation on the right hand side of
EQ (13) is broken down into four intervals corresponding to each
face. The temperature gradient normal to the control volume surface
is approximated by the temperature difference between the centroids
of adjacent nodes divided by the distance between them as shown
in Figure 3. Note that the line connecting the nodes must be perpendicular
to the control volume surface for an accurate approximation. Using
a cosine projection does not always yield accurate results.
Figure 3: Centroid method approximations to the
heat conduction equation
We may extend this approach and formulate a more accurate approximation
to the heat flowing across the control surfaces by using more than
one interval for each face in the summation on the right hand side
of EQ (13). In this case, the right hand face is divided into two
intervals (see Figure 4) and the temperature gradient normal to
the control surface at the centers of each interval is approximated.
This approach could be considered as a second-order centroid method
with respect to approximating the conduction terms.
Figure 4: Centroid method with two intervals used
to approximate conductive heat flow across right face
We can use interpolation regions again to assist in evaluating
the temperatures at the points a through d in order to compute the
temperature gradient normal to the control surface interval. Consider
a square geometry with unit thickness and unit lengths for the sides
of the control volumes. The algebraic approximation to the heat
flow through the right hand face is given by:
The heat flow from node 2 to node 1 is now expressed as a function
of six nodal temperatures, rather than just nodes 1 and 2. This
expression can be algebraically rearranged to look like "conductors":
We have negative "conductors" connecting node 1 to nodes 5 and
3. This does not mean that negative energy flows between
these nodes! The heat flow between nodes 1 and 2 is given by a collection
of "conductors", not just a single term. The algebraic equations
are rearranged only to make the input compatible with thermal analyzers
such as SINDA/FLUINT. Once this data is read in by SINDA/FLUINT,
the numeric computations that are performed are equivalent to EQ
Notice that in EQ (14) only the terms for T2 and T1 can be rearranged
to look like a typical conductor. The terms for T3 to T6 were made
to look like conductors by artificially adding terms for T1 in EQ
(15). The four artificially added terms for T1 all sum to zero,
yielding the original equation.
To complete the energy balance equation for node 1, the heat flow
through the other faces must also be considered. When all of the
conductors generated from considering the heat flow through all
of the faces are summed together, the network shown in Figure 5
Figure 5: Conductors obtained using CV-FEM
In this simple example, heat flows between two nodes are represented
by negative conductors. The conductors representing the heat flows
across each face could be input to SINDA/FLUINT as is, or common
node-to-node terms can be summed together to yield a composite conductor
which represents components of heat flows between more than one
This method can be extended to arbitrary quadrilateral interpolation
regions, by computing the gradient of the temperature at the centers
of the control surface intervals using EQ (8). This method is called
CV-FEM . In the example of square control volumes, the composite
conductor is positive, since the negative conductors are offset
by larger positive conductors. This may not be so for skewed triangular
Even when the net algebraic expression yields only positive conductors,
they cannot be viewed individually as representing the heat transfer
by conduction between the two nodes that they connect. Positive
conductors exist between the center node in Figure 5 and the diagonal
nodes, however, no heat flows directly from the center node to these
nodes. If the value of the heat flow between two nodes is desired,
the original expression given by EQ (14) must be used with the predicted
For well formed interpolation regions, the generated conductors
will be positive. When regions become highly skewed, negative conductors
can result. Recall that the heat flow by conduction is expressed
by interpolating the temperatures across the element. When
elements become skewed, extrapolation results. For example,
consider the problem of determining a surface temperature for a
simple 1D finite difference problem as shown in Figure 6 by extrapolating
the temperatures at nodes 1 and 2. The expression for the surface
temperature contains negative terms. Extrapolation, rather than
interpolation, also occurs for 2D regions with large interior angles,
resulting in negative terms to represent the heat flow.
Figure 6: Extrapolation yields negative terms
The presence of large or highly skewed elements may generate negative
conductors, but this does not automatically indicate a high degree
of error. What is important is how well the actual temperatures
match with the assumptions of how the temperature varies over the
interpolating region. For quadrilateral regions, the temperature
is assumed to vary bi-linearly. For triangular regions, the temperature
is assumed to vary linearly.
Costello  has performed a series of "thermal patch" tests, where
an arbitrary shaped triangular element is imbedded inside a larger
rectangular region. The net heat flow across the patch is computed
by setting one side to zero and the other side to unity and solving
the network using the FEM generated conductors. This heat flow was
compared to the exact solution and was found to match exactly for
all interior angles chosen for the triangle, even those that exceeded
90 degrees and caused negative conductors.
In this case, the actual solution varied linearly across the patch
(zero on one side to unity on the other). The linear approximation
used for the triangle is sufficient to match this case exactly,
even with highly skewed elements that generate negative conductors.
Costello  also points out an apparent paradox with respect to
FEM. He considers a single triangular element as shown in Figure
7. A negative conductor will exist between nodes 1 and 3, since
the triangle has an interior angle greater than 90 degrees. Suppose
that node 1 is a sink at zero degrees and node 2 is a sink at 1
degree. It is desired to maintain node 3 at a temperature of 1 degree
by applying heat as needed. Since the temperature of node 3 matches
node 2, it is assumed that no heat is transferred between nodes
3 and 2. Since the conductor between node 1 and 3 is negative, heat
must be paradoxically removed, rather than added to node 3 to maintain
its temperature at 1 degree.
Figure 7: An inappropriate use of a single triangular
The problems arising with the application of this element to this
thermal problem can be explained by first considering the assumptions
imposed by the element. A triangular element assumes a linear variation
of temperature between the nodes. In this example, the temperatures
of all three nodes have been specified, and therefore so has the
temperature gradient within the element. For these particular temperatures,
the gradient is perpendicular to the face formed by nodes 2 and
3. For a triangular element, the gradient is constant everywhere
over the interpolation region. Isotherms for the assumed temperature
distribution are shown.
Had this triangle been imbedded in a larger network of elements,
the gradient shown might be reasonable. However, in this particular
thermal problem, the sides of the triangle are adiabatic, and the
actual isotherms must be perpendicular to the boundaries. The true
temperature solution differs significantly from the assumptions
used by the triangular element. A better approach would be to use
more than one element (bisecting the element into to two triangles
would yield a good result), or use the knowledge of the expected
form of the solution to generate a conduction term by hand.
One final observation should be made about the example in Figure
7. Because of the negative conduction term between nodes 1 and 3,
it is assumed that negative energy must flow between them, a physically
unreal situation. Again we must consider all of the conductors together
in order to compute the heat flow between nodes on an element. Consider
the control volumes for this example:
Figure 8: Heat flow through control volumes for
a triangular element
Given the prescribed temperatures, the gradient is perpendicular
to the right hand face and is constant everywhere in the element.
The temperature gradient is shown at the centers of the control
surfaces that separate the nodes. Qualitatively, it can be seen
that heat does flow from node 3 to node 1. It can also be seen that
heat flows from node 2 to node 3, even though they are at the same
temperature. If this element were in a larger network and this gradient
actually existed, then heat actually would flow across this surface.
In anticipating the solution to our example problem, however, we
would expect that the temperature gradient point from node 1 to
node 3, almost perpendicular to the gradient imposed by our choice
of element and temperatures. Heat must be removed from node 3 not
because negative energy is flowing from node 1, but because positive
energy is flowing from node 2.
While potentially useful in aiding understanding of the heat transfer
process, the electrical network analogy is not valid for all types
of thermal problems. Some thermal problems can be posed such
that they have an electrical network equivalent. However, not all
thermal problems are subjected to the constraints of the lumped
parameter electrical network analogy. It is not that negative
"conductors" are unrealistic, but that it may be physically impossible
to construct a lumped parameter electrical analogue of a given thermal
With respect to FEM, conductors should be viewed as a by product
of algebraically rearranging a valid first law expression for the
heat flow. They are conduction terms rearranged into "conductor"
form. Approaching thermal modeling with the fundamental first law
heat transfer equations in mind, and an understanding of how differential
and integral quantities are approximated by algebraic means can
avoid many pitfalls and provide modeling approaches that might otherwise
WITH FD METHODS
The derivation of the energy balance equation terms, (heat sources,
capacitance, and conduction) has been presented here for a simple
rectangular interpolating region. For arbitrarily shaped regions,
the Jacobian must be included in the integration to relate a dudv
element to the actual area. The gradients of the interpolating functions,
Ni, must also be converted from u,v coordinates to x,y,z
space. However, the purpose here is not to present the detailed
mechanics of generating energy balance equation coefficients using
FEM or CV-FEM. That is a task that is better suited for a computer
than by hand, and is covered in detail by many references .
The purpose here is to present the method in terms of an energy
balance approach so that it can become apparent how FEM can integrate
with FD methods, traditional solution techniques, and with radiation
analysis. FEM can be viewed as a technique where energy from a source
is apportioned to a node based on the relative closeness to the
center of the control volume. The apportioning functions are in
general arbitrary, but a simple form results from using the simplest
representation, linear interpolation. External sources, thermal
mass, and conduction are all viewed in the same manner.
The "element" in FEM is nothing more than a region of the problem
demarcated to use a particular set of apportioning functions. A
nodal point can still be viewed as representing a location in a
control volume that is used to represent its average temperature.
The centers of the elements surrounding a node may be connected
together and viewed as the control volume for a thermal node.
Furthermore, the option exists when formulating the energy balance
equation for a node to use either the step function apportioning
functions (FD), or the smoothly varying apportioning functions (FEM)
for any of the types of terms (heating source, capacitance, or conduction).
For example, consider the geometry shown in Figure 9.
Figure 9: Energy into a node may be assigned with
either smooth or step apportioning functions
A small strip heater is placed near node 6. The control volume
for node 6 can be viewed as the boundary connecting the centers
of the elements that surround the node. The thermal mass could be
computed using this control volume, or by using EQ (4). The energy
dissipated by the strip heater could be lumped completely into node
6, or distributed to the nodes that share the element using the
smoothly varying apportioning functions.
Using a FEM approach, the strip heater could be subdivided into
smaller regions, and then heat from each region distributed according
to the apportioning functions. But this is not a necessity. If this
part of the model was being constructed by hand, then it may be
that the improved accuracy is not worth the trouble to calculate
the distributed heat loads into all of the nodes near the strip
heater. In this case, lumping the all of the generated heat into
node 6 is the best use of engineering judgement.
Once an element is defined, a computer can generate the conduction
and capacitance terms, saving much tedious and error prone work
by hand. If part of the model is not represented by geometry, however,
then abstract nodes using FD generated conductors can be connected
to any of the nodes in the FEM mesh.
For example, suppose the mesh in Figure 9 represents a plate and
the connection from node 12 to node 11 represents a simplified view
of a heat leak from a mounting point on the plate. The heat leak
might be due to conduction through bracketry, insulating washers,
and bolts. Using a simple abstraction of the mounting hardware rather
than a detailed FEM model of each bracket component makes sense
for this model (although a detailed FEM model may have been used
to deduce the simple abstract representation).
The use of FEM does not exclude any other modeling approach. A
consistent view of FEM as an energy balance approach allows a thermal
model to consist of FEM generated and FD generated terms.
Where amenable to a geometric representation, FEM can provide significant
savings in the modeling effort. Where geometry is not available
or inappropriate, FD methods should be used. FEM should be viewed
as a modeling tool that complements, not excludes, the existing
As shown by EQ (2), EQ (4), and EQ (11), nodal energy balance terms
generated using FEM are directly compatible with SINDA/FLUINT. Any
1D, 2D or 3D element that uses linear interpolation will produce
compatible conductors and capacitances.
It is often thought that SINDA/FLUINT is a finite-difference electrical
analogy thermal analyzer. While this is an acceptable view, a more
accurate description is that it is a general purpose equation solver,
and can accept FEM generated coefficients without modification.
Thus the solution to many of the complaints about solvers bundled
with FEM codes is to simply not use them.
As discussed earlier, additional information must be input to SINDA/FLUINT
if particular node-to-node heat flows are to be calculated. A finite
element pre-processor could automatically generate the data files
and subroutines needed for these calculations. Rather than manually
computing the heat flow by obtaining the temperatures and conductors
between nodes in user logic, special routines could be supplied
as part of the standard SINDA/FLUINT library to perform this function
Another complaint about FEM is that the network is harder to adjust
in order to correlate with test data. This is true if individual
network elements are being modified by hand, since individual conductors
do not represent heat flows. Correlation is best handled by using
the heat flow routines to gain information about the model, and
then adjusting the thickness, conductivity, specific heat, or shape
of the original FEM geometry.
Other approaches have been formulated to integrate FEM meshes with
thermal analyzers such as SINDA/FLUINT, usually based on an attempt
to convert the FEM mesh to a FD representation using the element
Part of the motivation for the FD conversion is to generate a network
using a familiar approach, since many analysts have an adverse view
of FEM, but mostly to be compatible with existing radiation analysis
tools. Formulating the network based on the elements makes it easier
to interface the mesh to traditional tools such as TRASYS  or
Converting a FEM mesh to some sort of FD representation does allow
compatibility with existing thermal radiation analyzers, but does
have other drawbacks. As EQ (13) shows, the accuracy of centroid
methods depends on having the lines connecting the centroids perpendicular
to the sides of the elements. Significant errors can occur when
nodes are "skewed." Using a cosine scaling does not always produce
acceptable results as shown in Figure 10.
The temperature gradient in any direction is determined by the
projection of the temperature gradient vector. In the example shown,
the cosine of the gradient approximated by using the two nodal points
is significantly less than the actual gradient perpendicular to
the control volume boundary. Projecting the control surface area
instead of the gradient yields the same error. If the actual temperature
gradient lies in the same direction as the line connecting nodes
1 and 2, no error will result. Error increases as the temperature
gradient becomes perpendicular to the centroid-to-centroid line.
Figure 10: Error in heat flow when line connecting
nodes is not perpendicular
Post processing and model-to-model data mapping are also hampered
when using a centroid approach in that element values must be converted
back to nodal values. Nodal values are usually taken as the average
of the element temperatures that are connected to the node. This
produces an artificial smoothing of the data.
Another approach is to place an additional node at the center of
each element, and then connect this node to the original nodes with
smaller FEM elements. The node at the center of the element is used
for radiation heat transfer. This is a valid approach, but generates
many additional nodes. Using the mesh directly avoids these complications.
The most direct and beneficial method of interfacing FEM generated
meshes with traditional thermal analyzers is to use them exactly
as they are. The problem does not lie with incompatibilities
of the FEM mesh, but with the radiation analyzer. The best approach
is therefore to fix the problems with the radiation analyzer, rather
than deal with the problems and drawbacks of FEM-to-FD mesh translation.
WITH RADIATION ANALYSIS
The aerospace thermal analyst is not only concerned with conduction,
but more importantly, the effects from orbital heating and from
radiation exchange among the spacecraft surfaces. These calculations
often consume more computer resources than the solution step to
predict temperatures. The lack of integration with suitable radiation
analysis tools is probably the biggest reason for the lack of widespread
use of FEM techniques.
As discussed earlier, nodes on a FEM mesh may be viewed in the
traditional control volume fashion, or by using a smoothly varying
apportioning of energy. Radiation exchange can be performed with
codes such as TRASYS or TSS by using the traditional control volume
view of nodes on a FEM mesh. Control surfaces for each node are
generated by subdividing the element into smaller polygons as shown
in Figure 11.
Figure 11: Control volume view of nodes used to
generate radiation exchange surfaces in a FEM mesh
Node 1 in Figure 11 would have four separate polygons generated
for use by TRASYS or TSS. This approach will produce accurate results,
but efficiency will suffer and large database files will result
due to the many surfaces produced. Another approach is possible
that yields better accuracy and is consistent with FEM energy apportioning.
Using EQ (2), the effects of direct incident heating can be computed
directly for nodes on a FEM mesh. The method works well with both
area integration and Monte Carlo methods. For the area integration
technique, the element is subdivided into smaller regions, each
of which is considered as a point source. The method fits particularly
well with Monte Carlo methods, since each ray is considered as a
point source bundle of energy. The energy deposited by the point
sources are distributed to the nodes using the energy apportioning
functions associated with the element.
Radiation exchange factors can also be computed directly for a
FEM mesh, using a consistent energy apportioning approach . Just
as the interpolation functions are used to interpolate the temperature
across an element, EQ (7), they can be used to interpolate the radiosity
across an element (see Figure 12). Rather than the constraint of
an isothermal emitting surface, a FEM surface can emit non-isothermally,
yielding better accuracy for the same size nodal regions.
In this approach, a ray leaving an element represents an energy
bundle consisting of contributions from each of the nodes used by
the element. The interpolation functions are used to determine the
relative contribution of energy to the ray from each node. This
technique generates global node-to-node radiation exchange factors
directly and avoids the element to node conversion problems. This
method has also been incorporated into ESATAN/ESARAD .
Figure 12: Raytracing with finite elements
Alternatively, instead of considering a ray to be composed of energy
from each node, calculations can be performed for a single node
at a time if the interpolation function is used to weight the generation
of the emitting location on the element. Using this approach, each
ray will have a energy value for one node only, but more rays will
be fired closer to the nodal point.
Modeling radiation leaving a surface uniformly, as with traditional
isothermal approaches, is equivalent to a viewing a geometric model
using a "flat" shading model. Each facet appears as a single shade
of color. Using the non-isothermal FEM mesh, radiation leaving a
surface is equivalent to using a smooth "Phong" type of shading,
where a curved surface appears curved, rather than faceted.
In fact, the use of thermal aerospace radiosity techniques have
found recent application in the graphics world as a technique for
generating extremely realistic images . The problem in graphics
image generation is an almost exact analogue to the thermal radiation
problem. The only difference is that the quantity of interest for
image generation is the total energy leaving a surface (which strikes
your eye), rather than the total energy being absorbed. Many methods
have been formulated including one that incorporates a traditional
double area integration scheme using finite elements . This
is an area of active research and will likely spawn useful algorithms
that are applicable to thermal radiation modeling.
It should be reiterated that a radiation analyzer does not need
to be exclusively based on FEM. A program that accepts traditional
planar and conic "uniform-radiosity" surfaces and arbitrary shaped
"non-uniform-radiosity" FEM meshes can be constructed. The use of
FEM for radiation modeling does not exclude traditional techniques.
The development of a new radiation analyzer, RadCAD™ ,
that incorporates these features is currently underway.
BETTER MESH GENERATOR
Once a radiation analyzer is created that can analyze FEM meshes
directly and produce nodal heating rates and radiation exchange
factors, then the conduction and capacitance network created using
FEM can be used directly with SINDA/FLUINT.
However, one more drawback still remains. Using traditional radiation
analyzers, a cylinder can be modeled using relatively few circumferential
nodes. Using existing FEM based meshing utilities, a cylinder must
be modeled using enough planar elements to geometrically approximate
the cylinder. The nodal density is driven by geometric fidelity,
rather than thermal accuracy concerns. One usually ends up with
a thermal model with many more nodes that what one could obtain
using TRASYS or TSS.
This is a problem with existing FEM codes, however, not with
the finite element method. Existing FEM codes have been developed
with a library of element types that are primarily designed for
Elements specific to thermal aerospace applications can
be formulated. In the example shown in Figure 2, a simple bilinear
scheme was used to approximate the temperature variation over the
flat rectangular element. The same interpolation scheme could be
used for a rectangular region wrapped around a cylinder as shown
in Figure 13.
Figure 13: A flat element yields the same capacitance
and conduction terms as a curved cylindrical element
One can image a rectangular plate subjected to a certain set of
boundary conditions to produce a certain temperature distribution.
This distribution is unchanged if the plate is folded or curved
in any manner that does not "stretch" the surface. Topologically
the curved surface is the same as the flat surface with respect
to computing capacitance and conduction terms.
The finite element method may be applied to any region on
a curved surface (even "stretched" regions), as long as a suitable
method of interpolating the temperature distribution over the region
can be formulated. The set of familiar primitives currently in use,
such as cones, spheres, toroids, etc., may be used with the finite
element method. If such element types were available, thermal FEM
models could be formulated with no penalty associated with excessive
In addition, elements do not have to be restricted to be parallel
to the principle parametric directions of the surface. Current analyzers
allow surfaces such as cylinders to be subdivided along constant
height and angular lines only. Improved modeling flexibility could
be obtained by first specifying the general orientation and dimensions
of a cylinder, and then demarcating regions using a set of three
or four (height, angle) pairs. Three and four sided surface elements
are mapped to curved surfaces using parametric coordinates, rather
than x,y,z values. The radiation analyzer first performs
an intersection test with the high level graphics primitive, then
determines which element is hit, and then distributes the appropriate
amount energy into the nodes on the element.
Curved elements are simple to implement for thermal calculations,
but complicate things considerably when used for stress or dynamics
calculations. A folded piece of paper behaves much differently under
stress than does a flat piece of paper. This is the reason that
only flat elements are typically available in existing FEM codes.
A FEM based meshing program specific to thermal aerospace needs
could contain other types of specialized elements to aid in spacecraft
modeling. Just as there are specialized (and somewhat abstract)
stress elements such as dampers, springs, and point masses, specific
thermal elements could be created such as heaters, heatpipes, contact
conductance, and MLI blankets. An element can be used to approximate
a given differential equation over a given domain, and it can also
be used to denote parts of the model that are more abstractly represented,
or represented by procedures (user logic). Models are built using
thermal modeling "objects", some that generate data and SINDA/FLUINT
logic, some that use FD techniques to generate thermal mass and
conductors, and some that use FEM techniques to generate capacitances
and conduction terms.
These new FEM techniques are being integrated with traditional
modeling methods in a program called Thermal Desktop® .
Thermal Desktop will be the first system to integrate a finite element
method specifically tailored to aerospace thermal analysis, traditional
FD methods, and non-geometric modeling methods into a single CAD
The use of FEM for aerospace applications has been limited for
a number of reasons, foremost being the inability of currently available
FEM based codes to produce efficiently sized models, supply necessary
procedural modeling capabilities, and provide adequate radiation
analysis capabilities. Current suppliers of FEM based systems have
little motivation to expand the capabilities of existing codes due
to low demand and a generally adverse view of the finite element
method by the aerospace community.
Previous approaches to integrating FEM meshes with existing tools
have been hampered by the incompatibility with existing radiation
analyzers, not by a fundamental incompatibility with existing thermal
analyzers. A better solution is the fix the incompatibilities with
current radiation analyzers, rather than transforming an otherwise
accurate and useful FEM mesh.
Viewing FEM from first law principles adds an existing modeling
technique to the thermal analyst's toolbox, in conjunction with
current procedural and FD based modeling practices. A closer look
at the physics involved and the algebraic approximations used explains
the phenomenon of negative "conductors." Conductors generated using
FEM and higher order FD methods do not represent the heat flow between
the two nodes that they connect. The conductors represent components
of heat flows between many nodes.
Adding the capability for curved elements on conic surfaces to
existing FEM based mesh generators will allow the creation of combined
conduction/radiation models with nodal densities comparable with
existing approaches. Allowing arbitrary orientation of elements
on curved surfaces and additional modeling elements such as MLI
blankets and heaters will provide new and powerful model building
Table 1 summarizes the theses of this paper, identifying the problems
often cited in the aerospace community regarding FEM codes, and
their resolution. Once aerospace thermal engineers are able to separate
the negative aspects of the current structural oriented FEM meshers
and solvers from the positive aspects of FEM, and are provided with
tools that are truly appropriate to the unique aspects of thermal
control, they will find that they have gained CAD integration and
concurrent engineering, and have lost only the task of generating
models by hand.
(lack of first law basis, or violation of second law)
lack of procedural
continue to use
codes such as SINDA/FLUINT with FEM generated networks
models too large
for thermal problems
elements (Thermal Desktop®)
continue to use
codes such as SINDA/FLUINT
Thermal Desktop development has been funded by the NASA Marshall
Space Flight Center, with Mr. Bill Till serving as technical monitor.
1. Cullimore, B.A., "Advances in SINDA/FLUINT
and SINAPS," 26th ICES Conference, 1996
2. Kreyszig, E., Advanced Engineering Mathematics,
Fourth Edition, John Wiley & Sons, New York, 1979
3. Schneider, G. E., and Raw, M. J., "A New Control-Volume-Based
Finite Element Procedure for the Numerical Solution of the Fluid
Flow and Scalar Transport Equations," Final Report to Atomic Energy
of Canada Limited, Chalk River, Ontario, January 1986
4. Costello, F. A., and Costello C. F., "Node
Geometries and Conductances in Spacecraft Thermal Models," SAE paper
951698, 25th ICES Conference, 1995
5. Zienkiewicz, O. C., and Morgan, K., Finite
Elements and Approximation, John Wiley and Sons, New York, 1984
6. Lockheed Engineering and Management Services
Company, "Thermal Radiation Analysis System Users Manual," JSC 22964,
Contract NAS 917900, October 1991
7. Panczak, T., et al., "Thermal Synthesizer System:
An Integrated Approach to Spacecraft Thermal Analysis," SAE paper
911582, 22nd ICES Conference, July 1992.
8. Chin, J. H., Panczak, T. D., and Fried, L.,
"Finite Element and Raytracing in Coupled Thermal Problems," Proceedings
of the Sixth International Conference on Numerical Methods in Thermal
Problems, pp. 683-701, Pineridge Press, Swansea, U.K., 1989
9. Cook, G. M., Hodgetts, C., Stock, N. J., Planas
Almazan, P., "Using Finite Element Discretisation for Spacecraft
Thermal Analysis," SAE paper 951608, 25th ICES Conference, 1995
10. Sillion, F. X., and Puech, C., Radiosity
and Global Illumination, Morgan Kaufmann Publishers, Inc., San
Francisco, California, 1994
11. Zatz, H. R., Galerkin Radiosity: "A Higher
Order Solution Method for Global Illumination," Computer Graphics
Proceedings, ACM SIGGRAPH, 1993
12. Panczak, T., and Ring, S., "RadCAD: Integrating
Radiation Analysis with Modern CAD systems," Presented at the 26th
ICES conference, July 1996.
13. Panczak, T., and Ring, S., "True Computer
Aided Thermal Engineering: A Geometric and Nongeometric Simulation
Environment," Final Report for NASA/MSFC, Contract NAS8-40666, July