The Finite Element Method and Thermal Desktop

ABSTRACT

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:

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.

INTRODUCTION

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 addressed.

FIRST LAW FEM

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 unrealistic.

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 a computer.

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 identical temperatures.

Figure 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 drastically.

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 manner.

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 region).

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 step function.

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 N_i for the apportioning function for node i. We can also use Green’s theorem [2] to simplify the Laplacian in EQ (5) into a gradient. The energy into a node by conduction is therefore also given by EQ (6).

Note that the gradient of the apportioning function, N_i, 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 by k).

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, G_ij is the same quantity as G_ji. 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 unrealistic.

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 [1]. 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:

EQ(14)

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":

EQ(15)

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 (14).

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 results.

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 node.

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 [3]. 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 geometries.

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 temperatures.

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 [4] 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 [4] 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 FEM element

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 problem.

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 be missed.

COMPATIBILITY 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, N_i, 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 [5].

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 FD approach.

INTEGRATION WITH SINDA/FLUINT

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 automatically.

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 centroids.

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 [6] or TSS [7].

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.

INTEGRATION 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 [8]. 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 [9].

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 [10]. 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 [11]. 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™ [12], that incorporates these features is currently underway.

A 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 structural analysis.

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 elements.

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® [13]. 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 based system.

SUMMARY

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 capabilities.

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.

Table 1: Summary

ACKNOWLEDGEMENTS

Thermal Desktop development has been funded by the NASA Marshall Space Flight Center, with Mr. Bill Till serving as technical monitor.

REFERENCES

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 1996

For more information on this subject, please contact the author:

Tim Panczak
C&R Technologies, Inc.
9 Red Fox Lane
Littleton, Colorado 80127-5710

About Us | Products | Services | Support | What's New | Resources
Home | Request a Quote | Site Map | Feedback

Copyright ® 1999-2006 Cullimore & Ring Technologies, Inc. All rights reserved.