THERMAL ANALYSIS OF N94-23650
COMBINATORIAL SOLID GEOMETRY MODELS
Capt Diane Gerencser
Mr George Radke
Capt Rob Introne
Phillips Laboratory, Kirtland AFB, NM
Mr John Klosterman
Mr Dave Miklosovic
Battelle, Columbus, OH
Algorithms have been developed using Monte Carlo techniques to determine the thermal
network parameters necessary to perform a finite difference analysis on Combinatorial Solid
Geometry (CSG) models. Orbital and laser fluxes as well as internal heat generation are modeled to
facilitate satellite modeling. The results of the thermal calculations are used to model die infra-red
(IR) images of targets and assess target vulnerability. Sample analyses and validation are
presented which demonstrate code products.
CSG models of targets have been used for many years in performing various phenomenologic
analyses such as nuclear particle transport, kinetic energy weapon effects, survivability, susceptibility,
and lethality studies. Many CSG vehicle target description databases are under construction and
include highly detailed, three dimensional solid geometry models of red and blue, high value strategic
and tactical targets. Ideally, these data sets include all of the information necessary to perform
detailed thermal analysis for thermal performance simulation. Applications of thermal analysis on
these data sets include thermal design analysis, laser weapon effects survivability analysis, and IR
image synthesis for sensor and seeker performance of target susceptibility, acquisition and tracking,
simulation of system IR imaging, and automatic target recognition algorithm training.
SINDA and other network methods for thermal analysis algorithms have been successfully
employed to model many complex thermophysical systems. Their success may be partially credited
to the simplicity of modeling very complicated thermal control, radiative transfer and nonlinear
thermal problems. However, these codes are not usually linked directly to geometric solid geometry
models, and therefore the burden of constructing the input data has been placed on the code user.
Often, the only cross correlation for temperature predictions is a sketch or drawing and a node
number penciled in by hand.
At the Phillips Laboratory and Battelle Memorial Institute, algorithms have recently been
developed to marry the two technologies, allowing high resolution thermal analysis of these CSG
target sets for such studies. The algorithms include three dimensional geometry discretization,
thermal network parameter computation including interelemental radiative transfer computation.
orbital and external thermal environment simulation, transient source computation including internal
heat generations and thermal control. The computed data is then reformatted to a SINDA input file
format for transient or steady state thermal predictions. The resultant thermal profiles are explicitly
linked to the geometric models. This correlation increases the usefulness of the data by allowing IR
image and signature predictions from the total set of geometry, thermophysical, thermoradiative and
temperature data. Although originally created for imagery analysis, the algorithms were modified to
include time and space resolved laser fluxes. The thermal results of the laser loads on a target can be
used in vulnerability analysis. All of the algorithms are combined in a suite of codes collectively
THERMAL SOLVER CHOSEN
SINDA was chosen as the thermal solver for AutoSINDA for a variety of reasons. Its
application is flexible, able to model time and temperature dependent parameters, convection,
conduction, internal and external radiation, heat pipes, thermal control systems, and phase changes.
SINDA also allows programmable access through user defined variables and has many solution
schemes which allow both transient and steady state solutions. In addition, SINDA has been used and
tested for many years and has a heritage as a thermal analysis code in the satellite community. The
SINDA code has a broad base of users within the community from which to draw expertise.
An overview of the AutoSINDA codes is presented in Figure 1. The model is first discretized.
The discretized model has the same format as the original model and has lost no geometric fidelity.
Next the thermal network parameters are calculated using a Monte Carlo technique. Orbital loadings
are applied using only first incident absorption of solar and earth fluxes. Diffuse fluxes such as earth
shine is modeled as a distributed impressed flux. The result is an element/load table as a function of
time. AH thermal network parameters and loads are reformatted for SINDA input. SINDA is then
used for transient thermal analysis. There is a one to one correspondence between the discretized
e ements in the solid model and the SINDA nodes. The temperatures correlated to the discretized
geometry are used to determine laser effects or the resulting IR image of the target.
"^ e AutoSINDA model is limited by the following assumptions: 1) all materials are opaque
2) all materials are perfect diffuse or perfect specular, and 3) only first incident external radiative flux
absorption is modeled.
CSG Geometric Description
A combinatorial solid geometry model is described by Boolean combinations of solid
primitives. The lowest hierarchy unit of these models is called the "element" (not related to element
as in Finite Element). Each element should be composed of a contiguous solid material and have a
single surface property or coating. An example of a single element would be a base plate for an
antenna or a panel for a structural housing.
Each element is a built by a combination of solid primitives. This entails complex intersections
and unions of void and solid ellipsoids, elliptical cones and polyhedra. Note that in general these
combinatorial intersections of volumes allows very complex solid descriptions in a compact format.
The target model is comprised of many elements assembled in a series of hierarchial sets to
describe the complete system model. Many of these models are quite detailed, modeled down to the
electronic card level.
Algorithm Basic Processes
In order to perform thermal analyses on the CSG descriptions, five basic processes are
executed on the geometry model. These are:
Process 1: Model Discretization
Process 2: Thermal Network Parameter Computation
Process 3: Environmental Simulation
Process 4: Application of Environmental Flux & Internal Heat Generations
Process 5: Data Reformatting to SINDA Input Deck
For each element in the CSG model, a bounding box is found using the minimum and
maximum x, y, and z coordinates for the element. The number of divisions in each direction which
will produce a length less than or equal to the maximum mesh size, for that element is determined
n r “ + 1 (1)
where r is the axis being divided and "int" takes the truncated integer of the expression. Each axis is
divided into n, equal segments. A plane going through each division is cut perpendicular to the axis.
Planes which cut through the x, y, and z axis intersect to form new boxes whose sides are at most
long. Each box created that overlaps with any part of the original element is "ANDed" with the
original element using boolean algebra. This creates a new element in the discretized model. This
algorithm discretizes all elements in this "rectangular" manner and may create "strange" or very small
new elements. If an element from the original model has a convex surface, this discretization scheme
is undesirably susceptible to creating elements that are in two discontiguous pieces. In this case,
increasing the discrete resolution will reduce this anomaly.
Thermal Network Parameter Computation
A Monte Carlo technique is employed to determine volumes, surface areas, conductive
resistances, and radiative resistances. For each element, I, a sphere is found which completely
surrounds the element. The sphere is found by determining the minimum and maximum coordinates
for each primitive in the element (an element can be composed of several primitives which are
"ANDed" and "ORed" together). The maximum coordinate for the element is determined by the
maximum coordinate of all element primitives that are "Ored" and the minimum maximum coordinate
for all the primitives that are "ANDed". The minimum coordinate for the element is found similarly.
Note that this technique is very simply and efficiently evaluated but does not necessarily find the
smallest sphere containing the element, which may be extremely complex.
The radius of this sphere is determined by,
U _ ~ + fomc ~ * few ~ ^mfa) (2)
The center of the sphere is located at the average coordinate.
The cross-sectional area through the center of the sphere is given by,
A — -
Rays are shot in an inward random direction from a random location on the surface of the
surrounding sphere. A total of n^ are shot. The number of times a ray hits the surface of the
element, n^,, is recorded. Every time a ray hits the surface of the element, there is one entrance and
one exit point. For every hit, the distance traveled through the element from the entrance to the exit
point, dL, is recorded. The total length traveled through the element is determined by summing the
distances traveled through the element from each hit.
1 - 1
The volume of the element is determined by,
( 6 )
The surface area of the element is determined by,
A *J "
For every entrance and exit point that a ray found, a ray is shot from that location on the
surface of the element in the outward normal direction. The ray is tracked until it either hits another
element or escapes to space. If the ray hits another element, k, the distance traveled to the element k
is determined. If the distance traveled is less than a tolerance distance for conduction, the point
where the ray originated is located on a conductive surface. The number of times a point on a
conductive area is found from the element I to element k, is collected in n^,^. If the origination
point is located on a conductive surface, a ray is shot from that location in the inward normal
direction until it exits the element. The distance traveled to the exit of element I is recorded in the
array d k . When the ray exits element I, a check is made to see if it conducts to the same element, k,
at that location. If it does, d k is divided by two. The distance to the conductive center from element
I to k will be determined by half the average value of d*. The reason d k is divided by two if the ray
conducts to element k on both sides of element I, is that instead of heat flowing through element I,
heat will flow in the inward/outward direction from the element.
For every entrance and exit point found that is not located on a conductive area, the energy
that is. emitted from the element to all the other elements and space is found. A ray is shot in a
random outward direction. The ray is given an initial energy weight, w, of:
(O - Otf
( 8 )
The ray is tracked until it either hits another element or escapes to space. If the ray escapes to space,
the remaining weight, «, is added to and the ray is terminated. If the ray hits an element, k,
the energy deposited in E,.* is determined by,
E lk " E i* + e *“ <9)
The weight of the ray is reduced by the energy deposited in the element it hit,
o ■ u - e k a
( 10 )
The ray continues in a random direction from the point it hit on element k until an exit criteria is met.
Exit criteria are in the form of a maximum number of reflections made by the ray and fraction of
original weight left in the ray.
The thermal network from element I to all other elements and space is determined. For every
element k that element I conducts to, the conductive area between I and k is found by.
The average length, l^v.,!-*. from the "center" of element I that is normal to the conductive area,
Acood.i-k is found by.
( 12 )
The non-linear radiative conductance (inverse of resistance) from element I to element k is contained
Upon completion of this algorithm for every element, the arrays from i to j and j to i are
combined. Since the conductive area between element i and j is computed separately for element i
and element j, they must be combined. A weighted average between the two areas is taken. They
are weighted by the cross-sectional area of the other’s surrounding cylinder.
( 13 )
w i -
^ — t , + ^-crottj
W 2 "
■^cradHi + ^cmtj
Ay - ^ee njhj * yv ‘A*mdJM
The conductive resistance from element i to element j can then be determined by,
where is the thermal conductivity. The "resistive" conductor between element i and element j is
found by adding the two resistors in series from the elements to the conductive area,
+ Rjm (^
The conductive heat flow between element i and j is then found by taking the inverse of the
To find the radiative heat flow between element i and j, a weighted average between and
Ej^j is taken. They are weighted the same way as the areas,
E H “ w i e h + ( 16 )
Environmental simulation is accomplished using a standard orbital code which reads 'the orbital
elements and epoch of the satellite and writes the satellite ephemeris to a data file. This information
includes the direct solar flux and direction, the earth shine flux, the solar albedo flux and the
orientation of the satellite with respect to the earth.
Application of External Fluxes & Internal Heat Generations
The direct sun, laser flux, and earth albedo are modeled as collimated loads. Although the
earth albedo is by no means collimated, this model is implemented for computational simplicity. This
model does conserve the energy of the earth albedo, primarily making the gross approximation in the
Earthshine is modeled in a distributed fashion. The distributed earthshine is modeled as
originating from a flat earth disk, broken into 13 sections of equal area. Each area radiates one
thirteenth of the total earth flux on the space object. The disk radius is defined as,
r d - sinQ (17)
0 - aumfa? + 2 a
( 18 )
altitude of satellite
( 19 )
and r e is the IR apparent radius of the earth,
r t - 6378.2 fan + 30 hn(IR apparent atmosphere) (20)
The disk is broken into an inner disk, a four-piece middle annulus, and an eight-piece outer annulus.
The area of each of the thirteen sections is,
( 21 )
The radius of the inner circle is.
( 22 )
The outer radius of the middle annulus is,
The middle annulus is sectioned at 0, 90, 180, and 270 degrees. The outer annulus is sectioned at 0,
45, 90, 135, 180, 225, 270, and 315 degrees.
The effective directions for each section to impinge the satellite is modeled deterministically as
coming from a point at the center of each sector. The center is defined by the area average radius for
the annulus sectors and the center of the circle for the inner circle sector. Each direction is assigned
one thirteenth of the total earthshine flux.
The flux and direction of all external loads are applied to the solid model by ray tracing.
Currently only first-incidence actions of externally applied loads are modeled. The absorbed flux
going into each element is numerically evaluated by ray tracing the entire model in a semi-systematic
uniform grid. The grid is defined by the square which surrounds the circle defined by the model-
enveloping sphere. The enveloping sphere is defined as the sphere centered on the models bounding
rectangular box with a radius extending to the bounding box’s comer. Ray tracing is accomplished
through a user-specified mesh size.
The radiant heat transfer into each element from loads other than laser loads is calculated by
multiplying the number of hits on the element, the area each ray represents, the spectral (or banded)
absorptivity of each element, and the uniform incident flux of the load. Since laser loads do not have
uniform flux, the laser model employs rays with an energy equal to its representative area times the
radial variation of flux according to the gaussian laser profile,
where q,, is the peak flux (central intensity) of the laser, r is the radial distance from the ray being
traced to the aim point of the target, and a is the Gaussian deviation parameter for characterizing the
Internal Heat Generations
Internal loadings can be specified as a constant heat generation or an on/off thermal control
heat generation. The internal heats are assigned to the geometry model’s original elements before
discretization and are then distributed to the discrete elements based on mass fraction.
Data Reformat to SINDA Input Deck
Upon completion of the four previous processes, all necessary data to perform thermal analysis
now exists. The data is formatted to a SINDA compatible input deck and subsequently submitted to
the SINDA process. This reformatting operation allows the use of both temperature dependent or
independent material properties. The option to cyclically repeat a single orbit is also allowed. Any
of the SINDA transient solution schemes may be used as well as the steady state solvers.
Conductive Model Validation
Case 1: Solid Cube with a Zero Temperature Boundary Condition
Two samples cases were created to test the validity of the conduction method, the first sample
case is a box with sides length 5 cm. The temperature of the box is initially 100 K when a zero
temperature is imposed at the boundary for times greater than 0. In order to apply the above
algorithms, a thin shell .1 cm thick is modeled around the box. The analytical solution for a box with
sides of 5 and 5.1 cm is compared to the solution predicted by the algor ithms f or differe nt mesh sizes
in Table 1. In the table, n is the number of divisions in the x, y, and z directions.
As the mesh size gets finer (n gets larger), the algorithms approach the analytical solution with
the exception of the finest discretization. The final column in Table 1 shows a divergent behavior as
the number of divisions was increased to 17 in each direction. This is probably caused by numerical
Case 2: Sphere with a Zero Temperature Boundary Condition
The second sample case is a sphere with radius 5 cm initially set to 100 K when a zero
temperature is imposed at the boundary for times greater than 0. Again, in order to apply the above
algorithms, a thin shell 0.1 cm thick is modeled around the sphere. The analytical solution for the
sphere with radii of 5 and 5.1 cm is compared to the solution predicted by the algorithms for different
mesh sizes in Table 2.
Table I Comparison of the temperature at the center of the box between the analytical solution
and the numerical approximation
n= 1 1
Radiative Transfer Validation
In order to validate the radiative transfer model, a simple test geometry was constructed using
the solid geometry modeling software. A 6 sided hollow box was created 3 x 4 x 5 cm in internal
dimensions x, y and z respectively. Note that the radiative transfer mechanism accounts for both
internal and external radiative heat transfer, so it was necessary to ensure that this box was
constructed with no external view factors. Two methods were used to compute the radiative heat
The first method involved using a 6 facet analytical model for 6 simple view factors. This
method assumes uniform distribution of all irradiance and radiosity. The method of Gebhart was used
to find the radiative interchange factors and the matrix.
The second method involved discretizing each face of the box and numerically computing view
factors for each small element within the faces. Upon determining the hundreds of view factors, the
method of Gebhart was used to compute the radiative transfer factors. The large matrix of radiative
interchange factors then had to be summed over each face using the appropriate summation principles.
Note that the AutoSINDA algorithms inherently model the second case more accurately as it does not
assume the uniform irradiance and radiosity distribution.
The two distinct methods were chosen to allow complete comparison, knowing that the second
method models the radiative transfer more accurately. Table 3 summarizes the model of the 6 sided
box and its 300K broad band radiative properties.
Case 1: 6 Node Analytical Radiative Transfer
The internal radiative interchange factors for this geometry, K^, were computed using a
simple 6 node view factor analysis of the 6 sides of the box. This was done analytically using the
view factor formula for two right plates with a common edge. The method of Gebhart was used to
Table II Comparison of the temperature at the center of a sphere between the analytical
solution and the numerical solution
mathematically compute the radiative interchange factors and the resultant K ra(1 matrix.
Case 2: 6 Node, Discretized Radiative Transfer
The internal radiative transfer was modeled by discretizing each face of the box. Ideally this
will increase the accuracy of the radiative interchange. However, since the AutoSINDA radiative
transfer algorithms only tabulate this data in the 6x6 matrix, this very large set of data had to be
summed over the subelements within each face to back out the net radiative interchange for each ot
the six faces.
Radiative Transfer Validation Results
AutoSINDA was executed on the 6 sided solid geometry model using 10 thousand rays per
element, allowing a maximum of 30 reflections per ray. A mixed mode reflective model was chosen
allowing both photon analog radiative transfer propagation as well as the continuous ray progression.
The execution time for this case is on the order of 30 seconds wall clock time for a SunSPARC 10.
iiillid I till
Table III Physical Model of the 3x4x5 cm Box for Validation
Cast Aluminum Alloy, Type 360
Aluminized Kapton 1, Kapton side out
Carbon Black Tedlar Film
Highly Reflective Aluminum Alloy,
White Paint, Chemglaze A276
Black Paint, Chemglaze A382
Table IV shows the results of AutoSINDA in comparison to both of the methods described
above. Overall, the mean deviation for the six node view factor method is about 4% with a
maximum deviation of 24%, so the comparison against the simple 6 node model is quite unacceptable.
However, using the more accurate many node model, the mean deviation falls to only 2% with
a maximum deviation on the order of 5%, an increase of almost 5 times the accuracy. Hence it is
concluded that the AutoSINDA algorithms compare very well with the more accurate method of
computing the radiative transfer.
SAMPLE OF COMPLICATED GEOMETRY
As a sample of a complicated geometry, Figure 2 shows a spiral antenna attached to a flat plate
which has been subjected to orbital and laser loadings. This example is not presented to compare to
analytical or measured data, but merely to show the range of applications and the usefulness of data
Table IV Results of the Radiative Transfer Validation
[W/K 4 ]
6 Surface Elements
600 Surface Elements
7. 16098x1 0 12
5. 04044x1 O' 12
1.65488x1 O’ 13
4.56767x1 O' 12
1.901 17xia M
1.9501 lxia 11
2. 09295x1 a 11
2. 18181x1a 11
5.39 125x1 a 13
5.11 129x la 13
5.75496x1 O' 13
AUTOSINDA FLOW CHART
Figure 1. Flow Chart of the AutoSINDA Process.