# Full text of "NASA Technical Reports Server (NTRS) 19940019177: Thermal analysis of combinatorial solid geometry models using SINDA"

## See other formats

THERMAL ANALYSIS OF N94-23650 COMBINATORIAL SOLID GEOMETRY MODELS USING SINDA Capt Diane Gerencser Mr George Radke Capt Rob Introne Phillips Laboratory, Kirtland AFB, NM Mr John Klosterman Mr Dave Miklosovic Battelle, Columbus, OH SUMMARY 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. INTRODUCTION 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. 223 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 named AutoSINDA. 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. ALGORITHM DEVELOPMENT Overview 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. Assumptions "^ 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 224 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 Model Discretization 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 by, 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. 225 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) ry The center of the sphere is located at the average coordinate. *c - 2 + y. (3) - ^ ^TTrm The cross-sectional area through the center of the sphere is given by, A — - (4) 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. -£ i-i (5) 1 - 1 The volume of the element is determined by, LA, ( 6 ) The surface area of the element is determined by, 1 ^ A *J " n, FtTfS (7) 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 226 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 2 A n ray* ( 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. A-ctmdJ-k (11) The average length, l^v.,!-*. from the "center" of element I that is normal to the conductive area, Acood.i-k is found by. 2 n E 4 amdjt f-1 ( 12 ) The non-linear radiative conductance (inverse of resistance) from element I to element k is contained in E,_*. 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. 227 ( 13 ) w i - croafj ^ — 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, ! sstL Ay (14) 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 resistance. 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 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 directional distribution. 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) where 228 0 - aumfa? + 2 a ( 18 ) and 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, * section 13 ( 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 229 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, 9 (24) 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 spot size. 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. SAMPLE CASES 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 round-off errors. 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. 230 Table I Comparison of the temperature at the center of the box between the analytical solution and the numerical approximation time (s) anal (5cm) anal (5.1cm) CO tl c n=6 n=9 n= 1 1 n= 17 1 69.8 72.0 88.2 79.7 74.2 71.6 73.8 2 26.9 29.0 33.2 21.7 29.8 28.8 30.4 3 9.8 11.0 0.0 m 11.2 10.9 11.7 4 3.5 4.1 0.0 2.3 4.2 4.1 4.5 5 1.3 1.6 0.0 0.7 1.6 1.5 1.7 6 0.5 0.6 0.0 0.2 0.6 0.6 0.7 7 0.2 0.2 0.0 0.1 0.2 0.2 0.3 8 0.1 0.1 0.0 0.0 0.1 0.1 0.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 transfer coefficients. 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 231 Table II Comparison of the temperature at the center of a sphere between the analytical solution and the numerical solution time (s) anal (5.0) anal (5.1) Ti- ll c n-6 n= 11 n= 13 1 99.6 99.7 93.3 97.1 99.4 99.5 2 88.7 90.0 78.5 84.7 91.1 91.2 3 68.9 71.2 62.8 68.5 75.1 74.9 4 50.6 53.2 49.2 53.2 58.3 57.9 5 36.5 38.9 38.0 40.5 44.2 43.7 6 26.1 28.2 29.2 30.6 33.1 32.6 7 18.6 20.4 22.4 22.9 24.6 24.2 8 13.3 14.7 17.1 17.2 18.3 17.9 9 9.5 10.6 13.1 12.9 13.6 13.3 10 6.7 7.7 10.0 9.6 10.1 9.8 11 4.8 5.5 mm mm mm mm 12 3.4 4.0 5.8 mm 5.6 5.4 13 mm 2.9 4.4 4.0 4.1 4.0 14 2.0 BEM 3.1 3.1 3.0 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. 232 iiillid I till Table III Physical Model of the 3x4x5 cm Box for Validation Box Face Outward Normal Face Area, cm z Material Broadband Emissivity 1 -X 20.0 Cast Aluminum Alloy, Type 360 0.200 2 +X 20.0 Aluminized Kapton 1, Kapton side out 0.740 3 -Y 15.0 Carbon Black Tedlar Film 0.800 4 + Y 15.0 Highly Reflective Aluminum Alloy, 2024-A1 0.030 5 -Z 12.0 1 White Paint, Chemglaze A276 0.881 6 + Z 12.0 Black Paint, Chemglaze A382 0.975 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 visualization. 233 Table IV Results of the Radiative Transfer Validation Box Faces AutoSINDA Results [W/K 4 ] Analytic Results 6 Surface Elements 600 Surface Elements Krad %e Krad %t 1-2 7.06995xia 12 7.22104xia 12 2.14 7. 16098x1 0 12 1.29 1-3 5. 04044x1 O' 12 5.10710xia 12 1.32 5.10193xl0 12 1.22 1-4 1.65488x1 O’ 13 1.68235xia 13 1.66 1.73681x10-" 4.96 1-5 4.520l9xia 12 4.55183xia 12 0.70 4.56767x1 O' 12 1.05 1-6 4.87479X10- 12 4.47 5.06076xl0 12 3.82 2-3 2.1I906xia" 2.18800xia 11 3.26 2.17028x10" 2.42 2-4 7.23058xia 13 7.20755xia 13 -0.32 7.32650xl0 13 1.33 2-5 1.901 17xia M 1.9501 lxia 11 2.58 1.93365X10 11 1.71 2-6 2. 09295x1 a 11 2. 18181x1a 11 4.25 2.14510x1a 11 2.49 3-4 5.39 125x1 a 13 5.66541xia 13 5.09 5.40254xl0 13 0.21 3-5 1.47630x1a 11 1.56077x1a 11 5.72 1.53137x1a 11 3.73 3-6 1.68516xia n 1.74621xia u 3.63 1.69953xlO u 0.85 4-5 5.11 129x la 13 5.14137xia 13 0.59 5.19854xl0 13 1.71 4-6 5.64379xia 13 5.75225xia 13 1.92 5.75496x1 O' 13 1.97 5-6 1.10285xia u 1.36538xia n 23.8 1.12036x10" 1.5 234 AUTOSINDA FLOW CHART Reformat SINDA output for postprocessing Display Temperatures on Solid Model Figure 1. Flow Chart of the AutoSINDA Process. 235