a 

a 

NASA Contractor Report 1851 17 

Thermal Evaluation of Advanced Solar 

Dynamic Heat Receiver Performance 

Roger A. Crane 

The University of South Florida 

Tampa, Florida 

a 

March 1989 

a 

9 

e 

a 

Prepared for 

Lewis Research Center 

Under Grant NAG3-851 

NIA 

National Aeronautics and 

Space Administration 

(NaSA-Ca-laSl17) TEEBfiAL EVALUATIOI OE N 89- 272 56 

AWbNCBD SOLAR DYNAHIC HEBT RECEIVER 

PERPORHANCE Final Report, Jan - Dec. (University of South Plorida) 131 pCSCL lOA Unclas 

1988 

63/44 0224747 
e 

ii 

Table of Contents 

e 

1 Introduction ........................................... 1 

1.1 Advanced Solar Dynamic Systems ..................... 2 

1.1.2 Advanced Solar Dynamic Stirling Cycle ......... 3 

1.2 Advanced Solar Dynamic Receiver Design ............. 4 

1.3 Project Scope ...................................... 6 

1.1.1 Advanced Solar Dynamic Brayton Cycle .......... 2 

2 Analysis Approach ...................................... 11 

2.1 Basic Cycles ....................................... 11 

2.2 Concentrators ...................................... 12 

2.3 PCM Enhancement .................................... 14 

2.4 Receiver Concepts .................................. 15 

2.4.1 Base Receiver Concept ......................... 15 

2.4.2 Heat Pipe Concept ............................. 16 

2.4.3 Pebble Bed Concept ............................ 17 

2.4.4 Cavity Heat Pipe .............................. 17 

2.4.5 Brayton Cycle Cavity Receiver ................. 19 

2.5 Hybrid Elements .................................... 20 

2.6 Model .............................................. 22 

2.6.1 Physical model ................................ 24 

2.6.2 Void model ..................................... 24 

2.6.3 Energy equation ............................... 26 

2.6.4 Computer Program .............................. 27 

3 Parametric Study ....................................... 28 

3.1 Base Line Receiver ................................. 28 

3.2 Heat Pipe Receiver ................................. 34 

3.3 Pebble Bed Receiver ................................ 39 

3.4 Cavity Heat Pipe Receiver .......................... 41 

3.5 Cavity Receiver .................................... 45 

4.1 Attainable Performance ............................. 50 

4 Discussion ............................................. 50 

4.2 PCM Thermal Enhancement ............................ 53 

4.3 Parametric Studies ................................ 54 

5 Conclusions ............................................ 56 

6 References ............................................. 57 
Table of Figures 

a 

Advanced Solar Dynamic Brayton Cycle ..................... 3 

Advanced Solar Dynamic Stirling Cycle .................... 4 

Weight Distribution for Solar Dynamic System ............. 6 

Incident Solar Flux Distribution ......................... 12 

Receiver Energy Distribution ............................. 13 

Finite Difference Element Arrangement .................... 24 

Base Line Orbital Turbine Inlet Temperature .............. 29 

Effects of PCM Enhancement on Base Line Turbine Inlet 

Temp ..................................................... 30 

Effects of Element Number on Base Line Turbine Inlet Temp . .......................................................... 30 

Effects of Base Line Receiver L/D on Turbine Inlet Temp . . 31 

Effects of Base Line Element P/D on Turbine Inlet Temp ... 32 

Effects of Base Line Element Number on Canister Temp ..... 32 

Effects of Base Line Receiver L/D on Canister Temp ....... 33 

Effects of Base Line P/D on Canister Temp ................ 33 

Effects of PCM Selection on Heat Pipe Orbital Vapor Temp . .......................................................... 34 

Effects of PCM Enhancement on Heat Pipe Vapor Temp ....... 35 

Effects of Element Number on Heat Pipe Vapor Temp ........ 36 

Effects of Receiver L/D on Heat Pipe Vapor Temp .......... 36 

Effects of Element Number on Heat Pipe Canister Temp ..... 37 

Effect of Receiver L/D on Heat Pipe Canister Temp ........ 38 

Effect of PCM Enrichment on Heat Pipe Canister Temp ...... 38 

Effect of Element P/D on Pebble Bed Turbine Inlet Temp ... 39 

Effect of Pebble Bed Receiver L/D on Turbine Inlet Temp . . 40 

Effect of Pebble Bed Element P/D on Canister Temp ........ 40 

Effect of Pebble Bed Receiver L/D on Canister Temp ....... 41 

Effects of Cavity Heat Pipe Element Number on Canister 

Temp ..................................................... 42 

Effect of Cavity Heat Pipe Apex Angle on Cavity Temp ..... 43 

Effects of Cavity Heat Pipe Receiver L/D on Cavity Temp . . 44 

Effects of PCM Enhancement on Cavity Heat Pipe Cavity 

Temp ..................................................... 45 

Cavity Receiver Orbital Turbine Inlet Temperature ........ 46 

Effects of Cavity Receiver Element P/D on Canister Temp . . 47 

Effects of Cavity Receiver PCM Absorptivity on Cansiter 

Temp ..................................................... 47 

Effects of Cavity Receiver L/D on Canister Temp .......... 48 

Effects of Cavity Receiver Element P/D on Canister Temp . . 49 

a 

e 
e 

1 

1 Introduction. 

A study has been conducted by the University of South Florida 

to evaluate elements of Advanced Solar Dynamic power generation. 

This effort has been designed to determine the thermal performance 

of basic receiver concepts and the effects of various phase change 

materials. The study was funded by the NASA-Lewis Research Center 

under contract NAG3-851 and conducted during the period of January 

1988 through December 1988 at the University of South Florida, 

Tampa, Florida. 

Associated with a significant increase in space exploration, 

scientific research, commercial development and potential military 

weapons systems, comes the need for increased capacity, space 

operable, electrical power sources. Photovoltaic power generation 

is nowa welldevelopedtechnology successfully utilized in numerous 

commercial and military satellites. The technology is especially 

attractive for space applications in terms of mechanical simplicity, 

utilizing no moving parts and easing probldms of deployment and 

maintenance. Their major limitation is relatively low energy 

efficiency, typically in the range of 8 percent. This is well below 

the levels encountered in many relatively simple terrestrial systems 

which regularly achieve efficiencies of 20 to 25 percent. 

As power requirements in space increase, the size penalties 

placed on design by the photovoltaic technology become less 

acceptable. Whereas current satellites have generally required 

less than 10 kWe, missions planned within the next 10 years [l] 

require one to two orders of magnitude increase. Advanced phases 

of the space station incorporate design plans for a number of solar 

dynamic cycle electrical generators providing a total power level 

of between 300 and1000 kWe. Otherunclassifiedlarge scalemissions 

include a GEO communications platform, an air/ocean traffic control 

system, a GEO payload delivery system and a micro g materials 

processing system. With a clear need established for an advanced, 

light weight, increased capacity power generation system, NASA has 

undertaken the development of a space solar dynamic technology. 

It is clear that with the development of such advanced technologies, 

substantial improvements can be made in system size and weight. 

Concepts currently under consideration are for applications in 

low earth orbit (LEO). These include solar powered Brayton and 

Stirling cycles operation at peak cycle temperatures of between 

1000 and 14.00 K. To ensure continuous power during the eclipse 

portion of the LEO, thermal energy storage is accomplished through 

the use of a suitable phase change material (PCM) . With proper 

selection of materials the use of latent energy storage will provide 

relatively compact, low weight storage ata uniformhigh temperature. 

Both fluoride salts and metallic alloys are being evaluated for 

these applications. Typically the fluoride salts provide la,rge 

e 

e 

e 

e 

e 

e 
2 

a 

0 

a 

a 

a 

a 

heats of fusion with suitable containment chemical compatibility; 

they are limited by poor thermal transport properties and large 

volume changes. Metallic alloys provide outstanding thermal 

transport properties and restricted volume changes but do not have 

demonstrated containment materials. 

Theworkdescribedherein addresses a portionofthat development 

effort. The goal is first to examine a number of alternative 

receiver designs to determine their basic thermal performance. 

Designs providing large thermal gradients or large temperature 

swings during orbit are susceptible to early mechanical failure. 

In addition significant temperature variations in the working fluid 

have an adverse impact on the thermal efficiency of the cycle. By 

reviewing the thermal performance of each basic design, the relative 

merits of the basic concepts are to be compared. In addition, the 

effect of thermal enhancement and metal utilization provides a 

partial characterization of the performance improvements to be 

achieved by developing these technologies. The results are then 

thought to be applicable in determining the justification for 

development of new technologies. Finally the analysis includes a 

basic geometrical parametric study of each of the designs. Near 

the optimal design point it appears that there is generally a trade 

off between thermal performance and receiver weight. The results 

of this study will show how design parameters thermally affect 

overall system weight. 

1.1 Advanced Solar Dynamic Systems 

Initially the ASD program gave consideration to three basic 

power cycles: the Rankine, Brayton and Stirling engines. An 

organic Rankine cycle, using toluene as 3 working fluid, was 

considered but determined to be of limited potential for advanced 

applications. Organics tend to break down chemically at elevated 

temperatures, limiting the cycle to about 700 K. Since 

improvements in thermodynamic efficiency can most easily be 

achieved by increases in peak cycle temperature, consideration 

of the Rankine cycle then turned to either using liquid metals 

a as working fluids or to finding an alternative stable, high 

temperature fluid. Preliminary design studies have shown that 

such systems are of significantly greater specific mass than 

their alternatives so that attention is currently concentrating 

on future development of the Brayton and Ltirling cycles. 

a 

1.1.1 Advanced Solar Dynamic Brayton Cycle 

A schematic of the Brayton cycle is shown in Figure 1. 

The cycle utilizes a mixture of xenon and helium as the working 

fluid. Earlier studies have indicated that a suitable 

combination of high and low molecular weight gases provides 

an optimummixturewith relativelyhighdensity and high thermal 

a 
a 

ORIGINAL PASE SS 

OF POOR QUALITY 

3 

0 

a 

a 

e 

0 

conductivity [2]. The gas is compressed in a single stage 

mixed flow compressor, heated as it ir passed successively 

througha recuperator and a receiver, expandedthrough aturbine 

and finally cooled by returning through the recuperator and 

passing through the heat exchanger. Energy is supplied by 

focusing solar energy from a parabolic collector into the 

receiver; waste heat is radia.ted into space. 

neat pipes 

heat ekchonqer -. 

Concentrator 

! 

i 

\ 

li I 

a 

a 

a 

a 

Figure 1. Advanced Solar Dynamic Brayton Cycle. 

1.1.2 Advanced Solar Dynamic Stirling Cycle 

The Stirling cycle is represented schematically in Figure 

2. The engine includes a piston and cylinder arrangement in 

which the working fluid moves alternatcly to chambers above 

and below the piston via a set of bypass ports in the cylinder 

walls. Theupperchamber isthermallyconnectedtotheradiator. 

The piston is attached directly to a linear alternator to 

provide foraverysimple systemof electrical powergeneration. 

The Stirling engine represents a technological step beyond 

the Brayton engine. Whereas the Brayton cycle is widely used 

in terrestrial applications the Stirling engine is still in 

its infancy. The advantage lies primarily in its high specific 

power in larger generator sizes. Results of an earlier study 

[4] show that a 7 kWe Stirling system will weigh about 20% 

less than a comparable Brayton system, 40% less than a planer 

photovoltaic array, and 45% less than a photovoltaic array 

with a focused concentrator; at 35 kWe the advantage increases 
4 

a 

e 

0 

e 

e 

e 

e 

a 

a 

to about 25% below the Brayton cycle, 50% less than a photovoltaic planer array, and 85% less than a photovoltaic array 

with a focused concentrator. NASA is currently testing a 

prototype Stirling engine to obtain additional data for 

evaluating its suitability for long term space applications. 

I - 

heat exchanger 1 1 1 1 - heat I )pes 

1 

Figure 2. Advanced Solar Dynamic Stirling Cycle. 

1.2 Advanced Solar Dynamic Receiver Design 

Both the Brayton and Stirling cycles are adaptable to a wide 

range of high power applications. In low earth orbit provisions 

are made for uninterrupted power generation during the eclipse 

portion of the orbit by storing excess energy in periods of solar 

incidence as a sensible and/or latent heat. A receiver containing 

the storage media and designed for use with the Space Station is 

shown in Figure 3. On the right side of the receiver is an 

aperture placed at the focal point of the concentrator. Energy 

passingthroughtheaperture is absorbed in theelementscontaining 

the PCM located on the inside perimeter of the receiver; a portion 

is stored for later use while the remainder is transferred to 

the working fluid flowing through a concentric passage inside 

the PCM. 

a 

a 
e 

5 

e 

e 

0 

GRAPH IT E 

PCWORKING FI UID TUB 

INLET MANIFOLD 

(ALUMINUM) 

e Figure 3. Space Station Receiver Design. [3] 

e 

0 

Design studies [4] initiated to determine the overall weight 

distribution within a typical solar dynamic system have clearly 

established design targets for substantial performance 

improvements. Results, shown in Figure 4, indicate that the 

receiver would likely be the single heaviest component with about 

45% of the overall system weight. This is followed by the 

concentrator, with about 20% of the system weight, and the 

structure and radiator, each constituting about 10%. In an 

attempt to reduce overall system weight, it is natural that 

attention focus first on the receiver. In so doing it is also 

important to note that in most receiver designs, the PCM constitutes only a small fraction of the actual receiver weight. 

Receiver weight distributions for the base case design are shown 

in Table I. Here the PCM constitutes only 70 kg out of an overall 

receiver weight of 461 kg. This is a meer 15% and clearly suggests 

substantial opportunity for systemweightsavingsthrough improved 

receiver design. 

0 

e 

0 

e 
6 

Figure 4 

_.-- --. '. 

/' 

ENGINE 

CONDITIONING 

Weight Distribution for Space Station Solar Dynamic 

System. 

1.3 Project Scope. 

T h i s s t u d y i s d e s i g n e d t o c o m p a r e t h e b a s i c t h e r m a l performance 

of five receiverconcepts. Thegoal is todetermine theefficiency 

of each in utilization of the PCM, the magnitude of the thermal 

gradients resulting, the extent of the resultant peak cycle 

temperature swings and their impact on overall cycle efficiency. 

Concepts were selected for analysis based partially on their 

respective merit, but also with the goal of including the widest 

possible range of distinct design approaches. Those selected 

include (a) a base case taken from the Space Station design 

effort, (b) a heat pipe configuration, (c) a pebble bed concept, 

(d) a Sanders Brayton cycle design and (e) a cavity heat pipe 

arrangement. The base case design is shown in Figure 3. The 

remaining concepts are shown in Figures 5 through 8 respectively. 
7 

TABLE I. BASIC CASE RECEIVER WEIGHT DISTRIBUTION 

Component I 

Phase Change Material 

Working Fluid Tube 

Containment Canister 

Manif old 

Formed Insulation Mandrel 

Backwall support Plate 

Multifoil Insulation 

Aluminum Outer Shell 

Baffels and Supports 

Aperture Plate 

Aperture Shield 

TOTAL 

Weight, kg 

70 

17 

132 

21 

32 

15 

47 

59 

23 

30 

l5 

461 

Figure 5. Heat Pipe Receiver. [3] 
0 

8 

0 

I 

I 

I 

0 

0 

0 

INLET 

N~~~~~~~ 

. APERTURE 

PACKED BED TUB 

FL CONTl O"TLET ENCAPSULATION 

MATERIAL, PCM / 

.ow 

UNMENT 

JBE 

TUBE DETAIL 

a 

~~ 

Figure 6. Pebble Bed Receiver. [3] 

0 

0' 

Each concept has been proposed by one or more outside contractors to NASA. However, the approach undertaken here is not 

to attempt to compare specific designs. These designs have been 

used to identify basic configurations. Each configuration has 

then been normalized to utilize the same group of PCMs and 

evaluated under a similar set of operation conditions irrespective 

of those originally proposed. The results are not seen as a 

basis for comparison of the respective proposals but instead are 

intended as a thermal comparison of the overall configurations. 

Inaddition to considerations of basic rnceiver configuration, 

optimal selection of those dimensional parameters which affect 

thermal performance will be of concern. Specifically this study 

will include a parametric study of each of the basic receiver 

designs to establish the range of performance achievable with 

each. While no direct trade-off will be performed between 

improvements in thermal performance and possible adverse effects 

on receiver weight, it is envisioned that the results might be 

useful in such a comparison. 
9 

e 

a 

0 

e 

e 

TES MODULE 

I36 EA) 

MODULE 

BRANDS ___ --- . -.-_ 

STIRLING 

HEATER 

HEAD 

Na VAPOR 

Na LlOUlD 

EVAPORATION 

DOME 

MlILTIFOIL 

INSUUTION 

Figure 7. Cavity Heat Pipe Receiver. [5] 

Thethirdconsideration in thisstudy is that of PCMselection. 

Thermodynamic considerations would normally lead to a design in 

which all heating would occur at a temperature as high as practical 

below that imposed by the physical characteristics of the high 

temperature structure. This has led to the search for a material 

having a melting point near that imposed by the strength limitations of existing high temperature alloys and simultaneously 

possessing a high heat of fusion. 

A preliminary list of candidate phase change materials 

indicates numerous possibilities. The choice is complicated by 

the fact that many of these materials do not have well defined 

thermo-physical properties [6]. The characteristics sought here 

are not generally duplicated in other applications so that there 

has not heretofore been an intensive effort to establish properties. interest makes 

measurement particularly difficult. At such high temperatures 

materials are frequentlyquite corrosivemaking containment itself 

a problem. This analysis has been undertaken to determine the 

relative merits in utilization of either thermal conductivity 

Moreover the high temperature range of 
10 

e 

enhanced salts or of metallic PCMs. By including direct evaluations of the thermal performance of the respective receiver 

concepts incorporating such alternative materials is possible to 

ascertain the potential for improved performance from each. 

Ultimately such results may offer insight or lend weight to the 

decision to 

I 

pursue development ofthemore promisingtechnologies. 

MULTlFOlL 

INSULATION 

\ \ 

RADIATION WINmWS 

SHIELD 

CORRUGATED 

iEAT EXCHANGER 

WALL 

L I 

Figure 8. Brayton Cycle Cavity Receiver. [5] 

e 
11 

2 Analysis Approach 

0 

0 

a 

2.1 Basic Cycles 

Performance parameters for the cycle have been established 

from typical Brayton Cycle conditions. The selection of the 

Brayton cycle appears natural in that the receivers selected for 

comparison were, with the exception of the Cavity Heat Pipe and 

the Stirling Heat Pipe designs, originally designed for that 

cycle. It is not thought that this assumption will significantly 

detract from the generality of the results. There appears little 

reason to believe that operating conditions would change dramatically with either cycle so that results are thought to be 

generally applicable. 

The basic cycle performance characteristics are shown in 

Table 11. The base design for comparison is a Brayton cycle 

engine producing 7 kwe from a continuous 25.9 kwt input. This 

corresponds to an overall cycle efficiency of 27%. Variations 

in the turbine inlet temperaturewill affect both cycle efficiency 

and power output; these effects may probahly be best handled by 

controlling the working fluid flow rate or bypassing a portion 

of the flow around the receiver. In any case these control 

mechanisms are not used for purposes of this analysis and both 

receiver inlet temperature and flow rate are maintained at a 

constant level. The coolant composition, inlet temperature and 

flow rate are taken as the nominal values used in the base design. 

This selection might be considered to favor the base case but in 

preliminary studies, designs have not shown a significant sensitivity to the selection. In any case it was desired to compare 

the respective designs under similar operating conditions so that 

utilization of base case parameters appear4 more reasonable than 

selecting conditions proposed with one of the advanced designs. 

TABLE 11. BASIC CYCLE PERFORMANCE CHARACTERISTICS 

Parameter 

Engine Power 

Engine Thermal Input 

Receiver Losses 

Solar Incidence Period 

Eclipse Period 

Working Fluid 

Coolant Flow Rate 

Coolant Inlet Temperature 

Specification 

7 kWe 

25.9 kwt 

3.1 kwt 

. 54 min. 

36 min. 

714 He 29% Xe 

0.25 kg/s 

800 K 
12 

Parameter 

Dish diameter 

Rim angle 

Slope error 

Aperture radius 

Aperture offset 

Dish misallignment 

2.2 Concentrators 

The performance of any receiver in which the PCM is directly 

irradiated by the incident solar flux will be strongly a function 

of the incident flux distribution. In the study reporte'd here 

all flux distribution have been determined from the CAV2 [7] 

computer program. The input parameters used in each of the cases 

are as given in Table 111. The receiver diameter and depth were 

selected to accommodate the PCM elements in each case without 

additional clearance. The elements were located at a distance 

from the aperture so that the incident flux would begin to become 

significant at the first PCM cannister and was extended to exactly 

the element length. The inside diameter of the receiver was set 

at the outer diameter formed by the PCM elements. 

TABLE 111. BASIC CONCENTRATOR CHARACTERISTICS - 

Specification 

3.59 m 

4 50 

1% 

0.091 m 

Om 

0 mr 

* 

a 

e 

a 

a 

e 

rn I 

Figure 9. Base Case Incident Solar Flux Distribution. 
a 

a 

0 

0 

e 

a 

0 

0 

13 

The resultant flux distribution calculated for the base case 

configuration is shown in Figure 9.' The distribution is shown 

for receiver arrangements with L/D ratios varing from 1 to 2 as 

used in the parametric evaluation. Clearly the flux is quite 

high near the aperture end of the receiver leading to difficulties 

in obtaining uniform PCM utilization. The problem is somewhat 

more severewith thehigher L/D ratios indicating a higher incident 

flux to compensate for the lower reflectance from the bottom of 

the receiver. Receivers utilizing a truncated cone configuration 

can provide a more uniform incident flux distribution and a more 

appropriate PCM allocation. 

simple in form, are somewhat difficult to represent numerically without introducing certain round off error. In order 

to smooth the function it has been decided to integrate the 

area under the curve allowing representation of the cumulative 

heat input as a function of axial position. 

has been performed using standard numerical techniques and 

resultant curve is shown in Figure 10. Also shown is a third 

order curve fit of the integral used in the receiver program. 

The technique allows for relatively accurate representation of 

the CAV2 results but also provides a somewhat clearer picture 

of the tendency for poor flux distribution. 

that 48% of the incident energy enters the first 20% of the 

receiver length. 

The shape of the incident flux curves, while relatively 

The integration 

It may be noted 

,I- ./ Je 

Figure 10. Receiver Energy Distribution. 
14 0 

~ 

'0 

0 

0 

0 

0 

0 

a 

2.3 PCM Enhancement 

While fluoride salts have beenconsidered as prime candidates 

for selection as the phase changematerial, they are characterized 

by a low thermal conductivity. In order to achieve a high level 

of thermal utilization and low temperature gradients designers 

have relied on avariety ofthermalenhancementschemes. Proposed 

schemes have included high conductivity metal fins extending 

through the PCM, wire meshes within the PCM, metallic foams in 

which the PCM can be placed and a graphite matrix in which 

individual PCM pellets are to be arrayed. Manufacturing of 

certain of these concepts represents an unproven technology. 

Also unproven is the long term integrity .:f the concepts during 

normal orbital cycles. Nevertheless there remains the need to 

develop techniques of improving the thermal performance of such 

materials and it would seem prudent at this time to evaluate the 

impact of specific material enhancement schemes so as to perhaps 

limit development to only the most promising concepts. 

The approach used here has been to develop a model based 

largely around a llrandomll metal mesh. It is probably easiest to 

envision such a mesh as being constructed from layers of wires 

lying in a plane. Stacked on the first layer will be a second 

layer of wires arranged so that the wires in adjacent lamina 

appear to form a square grid. The model is designed around a 

large number of such planes. If one of these laminas is viewed 

from the edge so as to be seen from the end of the wires, the 

wires and PCM may be considered to act to transfer heat in 

parallel. The thermal resistance associated with this plane is: 

2 

( 1 - E ) * h. pc,,, + c. * A. I; v,L, R, = 

Adjacent planes, when viewed from the same perspective, will 

appear as layers of PCM acting in series with the wire segments. 

The effective thermal resistance of these laminae are given as: 

I P-2.r +- 

kPcAi-kE,,, 0.5.kpcnr In( -)' 

R, = P 

The composite system will appear as a large number of such planes 

all acting in parallel with one another: 

11 - I 

---+- 

R, K,) IC", 

When this system is viewed from a direction normal to both layers 

of wires, each lamina will appear as a cylinder in parallel with 

PCM. The equivalent resistance of this arrangement is given as: 
0 15 

a 

0 

0 

0 

0 

0 

0 

It is assumed that these two orientations are randomly dispersed 

within the elements so that the overall effect is to produce an 

isotropic effective conductivity.' It is further assumed that 

these resistances will combine in the same manner as structural 

properties in a composite material. From random orientations of 

fibers in a composite materials: 

This approach has be utilized to evaluate the effect of adding 

the various enhancement materials. The results are shown in 

Table IV. In the listing, the density is given as the density 

of PCM per unit volume. Since voids are assumed to remain 

distributed within the solid phase the same value is reported 

for liquid and solid. Likewise, the specific heat reported is 

based solely on the PCM so that the enhailcement material does 

not add to sensible energy storage capacity. The thermal 

conductivity values include the effects of distributed voids 

within the solid and may be taken as the effective enhanced 

conductivity. 

Table IV. Thermal Transport Properties of Enhanced PCM 

Material 

LiF 22 CaF2 

10% BN 

18% BN 

10% Grap. 

18% Grap. 

Ge . 

Density 

kJ/m3 

2097.0 

1887.3 

1719.4 

1887.3 

1719.4 

5010.2 

Specific 

Heat 

kJ/kg K 

(Sol/liq) 

1.84/1.97 

1.84/1.97 

1.84/1.97 

1.84/1.97 

1.84/1.97 

0.38/0.39 

Thermal 

Conductivity 

kW/m K 

(Sol/ 1 i.r) 

.00413/.00173 

.00492/. 00254 

.00554/.00320 

.00643/. 00387 

.00825/. 00559 

.0405/.0424 

Heat of 

Fusion 

kJ/W 

753 

753 

753 

753 

753 

481 

Density 

gJ/m3 

~~ 

1.597 

1.421 

1.295 

1.421 

1.295 

2.410 

2.4 Receiver Concepts 

2.4.1 Base Receiver Concept 

that shown in Figure 3 for the space station. 

The base case receiver design is essentially similar to 

In this design 

a 
0 

Parameter 

Tube Length, m 

Number of Tubes 

Sidewall/PCM Thickness 

Coolant Tube ID 

Canister ID 

Canister OD 

Canister Wall Thickness 

E 1 emen t P/ D 

PCM 

PCM Mass 

16 

Specification 

1.8 

75 

12.5% 

0.00841 m 

0.01095 m 

0.015 m 

0.00071 m 

1.66 

LiF 22 CaF2 

96 kg 

0 

0 

a 

0 

a 

2.4.2 Heat Pipe Concept 

The heat pipe receiver design is shown in Fisure < 5. The 

receiver is constructed in three axial sections. The first 

of these sections is located nearest the receiver aperture and 

includes no PCM canisters. Solar energy entering through the 

aperture strikes the central tubes directly. Inside the tube 

is a heat pipe working fluid. The evaporated fluid then 

transports the energy into the remaining two sections. In the 

central section the tube is surrounded by PCM canisters, with 

a configuration much like those in the base case design. A 

portion of the energy from the heat pipe fluid is conducted 

into the PCM and stored as latent energy during periods of 

solar incidence. The remaining energy from the heat pipe is 

transported to the last section of the receiver where a heat 

exchanger couples the heat pipe to the engine. During the 

solar eclipse the heat pipe drops in pressure and temperature 

and begins to remove energy from the PCM and transports it to 

the heat exchanger in the last section of the receiver. 

Dimensions on the base case receiver are shown in Table V. 

a 
a 

Parameter 

Tube Length, m 

Number of Tubes 

Packing Fraction 

Coolant Tube ID 

Canister Wall Thickness 

PCM 

PCM Mass 

PCM Sphere OD 

0 

Specification 

0.91 m 

15 

70.0% 

0.82 m 

0.00127 m 

LiF 22 CaF2 

96 kg 

0.014 m 

0 

a 

a 

0 

TABLE VI. HEAT PIPE RECEIVER CHARACTERISTICS 

Parameter 

Tube Length, m 

Number of Tubes 

Sidewall/PCM Thickness 

Coolant Tube ID 

Canister ID 

Canister OD 

Canister Wall Thickness 

PCM 

PCM Mass 

Specification 

0.7116 m 

20 

12.5% 

0.0482 m 

0.0522 m 

0.0820 m 

0.00127 m 

LiF 22 CaF2 

91.5 kq 

2.4.3 Pebble Bed Concept 

The pebble bed receiver is shown in Figure 6. In 

17 

this 
0 18 

a 

0 

e 

0 

0 

Brayton engine with the inclusion of a suitable heat exchanger 

within the receiver. In this design the receiver is configured 

as a right cylinder. The solar incidence from a parabolic 

collector arrives at the receiverandpassesthrough an aperture 

at the bottom. Inside the aperture the flux strikes the 

evaporator section of the cavity. The vapor then transfers 

heat to condensing surfaces onboth the thermal storagematerial 

(TSM) and the engine high temperature heat exchanger. 

The thermal storage material is contained in a number of 

canisters arranged around the periphery of the receiver. The 

side walls of these structures are tapered inward toward the 

apex, so that as the canisters are placed side by side they 

bend into a closed cylinder configuration. This cylinder is 

then banded together and placed inside the receiver. The apex 

of the canisters are covered with 3. condenser/evaporator 

material, generally a wire mesh. During the period of solar 

incidence, this area receives a continuous heat input, gradually 

melting the TSM. During the eclipse portion of the orbit, 

energy is removed from the TSM by evaporating the heat pipe 

fluid from this area. 

The shape of individual canisters has been established as 

a method of thermal enhancement. In effect, if one were to 

consider the region between roofs, this region acts as a 

triangular shaped fin with the heat pipe fluid replacing the 

fin material. The arrangement provides for fins of extremely 

high effective conductivity and of relatively low equivalent 

mass. Increasing the number of fins will necessarily imply 

an increased number of canisters, so that careful trade-offs 

are necessary between thermal enhancement and system weight. 

The high apex shape may also serve to alleviate the problem 

of void management, but this has not yet been confirmed. 

Certainly the narrowing walls are reminiscent of geometries 

proposed for fuel management within zero ?g environments. 

However Marangoni forces may well counter these effects during 

the heating phase and a careful analysis will be required to 

ascertain the final deposition. 

From a thermal aspect, the problem is not so severe as in 

those designs which allow direct solar incidence onto the 

canisters. In those cases the heat flux at any point in the 

receiver was fixed by the collector and receiver geometries. 

If there were inadequate material to absorb the incident 

radiation, surface temperatures would rise and a thermally 

induced failure might occur. In the cavity heat pipe design 

it is not the flux on the canister surface which is fixed, but 

rather the temperature. The temperature throughout the cavity 

will be nearly uniform and adequate to drive the energy into 
e 

Parameter 

Tube Length, m 

Number of Tubes 

Packing Fraction 

Coolant Tube ID 

Canister Wall Thickness 

PCM 

PCM Mass 

19 

Specification 

0.91 m 

15 

70.0% 

0.82 m 

0.00127 m 

LiF 22 CaF2 

96 kg 

0 

t L,le TSM, neither more nor less. At any point where low thermal 

storage capacity might occur, the energy input would simply 

decrease locally and the corresponding energy would be stored 

elsewhere. Only inthe casewhere thevoid covered a significant 

fraction of the overall heat transfer'surface would a problem - 

arise. In this analysis it is assu:,lcd that massive 

coverage of the heated surface will be avoided. 

TABLE VII. CAVITY HEAT PIPE CHARACTERISTICS 

0 

2.4.5 Brayton Cycle Cavity Receiver 

a 

0 

e 

0 

void 

The Brayton cycle cavity receiver concept is shown in 

Figure 8. A major feature of this particular design are a 

number of innovative concepts to decrease reciever weight. 

Most notable is the method used for PCM containment. The 

overall configuration of the receiver is that of a truncated 

cone with corregations running from the base toward the apex. 

The PCM is held into the troughs of the coLL-egations by capillary 

action while being contained within the receiver by a suitable 

window across the receiver aperture. A second smooth conical 

section is situated outside the first, with the interstitial 

space then serving as a passage for coolant. This particular 

arrangement is interesting in that it provides for a relatively 

light structure and includes a relatively limited number of 

fabrication steps. 

Solar radiation which enters the receiver inpinges directly 

on the PCM. Since the PCM is normally highly transparent to 

the incident radiation, suitable addatives are to be included 

to ensure adequate absorption. While this approach is as yet 

unproven, it is interesting from a heat transfer viewpoint in 

that it allows the designer a degree of control over the 

absorbed energy distribution. An indirect consequence of this 

is that high element P/D ratios are not required to achieve 

effective distribution of the stored energy. The net result 

is an extremely compact receiver with assocated weight savings 

in structure and insulation. 
20 

a 

e 

a 

a 

e 

This concept is the only one of the designs considered 

which utilizes'the truncated cone geometry to distribute PCM 

with higher fractions near the aperture and a decreasing 

fraction at the more distant locations. By situating the PCM 

where it can be most effectively utilized a somewhat decreased 

mass is required. 

TABLE IX. BRAYTON CYCLE RECEIVER CHARACTERISTICS 

I Parameter I Specification 

Corrugation Length 

Number of Corrugations 

Height of Corrugation 

Corrugation Wall Thickness 

Base Radius of Truncated Cone 

Top Radius of Truncated Cone 

PCM 

PCM Mass 

0.654 m 

100 

.064127 m 

.001 m 

.386 m 

.154 m 

Lif 22 Caf 

96 Ks 

2.5 Hybrid Elements 

Withintheparametricstudiesperformedoneachofthe receiver 

concepts, four PCM concepts have been evaluated. These have 

included (a) a fluoride salt, (b) a fluoride salt enhanced with 

Boron Nitride, (c) a fluoride salt enhanced with ATJ graphite 

and (d) a metallic PCM. This grouping, while considering only 

a limited number of alternatives, is thought to bracket the range 

of thermal transport properties of most materials under consideration for Advanced Solar Dynamic appl ;cations at this time. 

As a part of this study several additional materials were 

considered and, while not included in the parametric evaluation, 

were evaluated for thermal properties. In this case NaF was 

selected as the fluoride salt with Germanium being selected as 

the metallic alternative. Enhanced combinations each included 

NaF with the following enhancement schemes: (a) 12% of volume 

consisting of Niobium fins, (b) 12% Nickel fins, (c) 18% niobium 

wire mesh, (d) 18% nickel wire mesh and (e) 18% pyrolytic graphite 

mesh. 

An alternative arrangement considered was to encapsulate 

Germanium into spheres and to fill the interstices in the packing 

with either a conductivity enhancement material, i.e. graphite, 

or with a nonencapsulated PCM, i.e. NaF. Using uniform sized 

spheres it is possible to achieve packings of spheres of about 

60%. By using 2 significantly different sphere sizes packing 

densities of up to about 80% are achievable. 

e 
21 

a 

a 

0 

a 

a 

a 

0 

A final arrangement considered included packings of uniform 

sized spheres of encapsulated Germanium with pyrolytic graphite 

fibers traveling along the major axis of the packing to enhance 

thermal conductivity. The interstitial spaces were then filled 

with NaF to completely utilized the available volume. 

The results of an analysis of each of these schemes to 

determine the effective thermal transport properties is shown in 

Table X. In each of these comparisons the ;cost important thermal 

properties are the effective thermal conductivity and the energy 

density. The effective conductivity is of course that hypothetical 

thermal conductivity which would produce the same temperature 

difference in a large sample as would the nonhomogeneous material 

under consideration. It may be viewed as indicative of the 

temperature gradients which each sample would produce in an 

operating system. The energy density is the latent energy stored 

per unit volume by the nonhomogeneous element. It may be viewed 

as an indication of the required size of receiver which might 

incorporate that particular PCM. 

What is not indicated within this comparison is the ease of 

development of each of the respective concepts. Whereas Germanium 

appears to be the choice candidate on the basis of energy density 

and thermal conductivity, problems of containment may suggest an 

alternative selection. However almost any other selection of 

Germanium in combination with other materials considered here 

tends to result in an improved energy density. The potential 

combination with NaF, if it could be fabricated, certainly 

indicates a potential for compact design. However this particular 

combination results in an unusually low thermal conductivity and 

it not thought to merit further consideration. 

The enhancement of NaF with either of the two graphites, ATJ 

or pyrolytic graphite, will serve to significantly improve the 

effective thermal conductivity. The use of pyrolytic graphite 

offers a significant advantage over nickel additions. However 

any addition of a nonmelting material will necessarily decrease 

the energy density. The choice of a fluoride salt with any of 

the thermal enhancement additions is a low risk choice in terms 

of technological development while offering satisfactory thermal 

conductivity improvement. Itwill, of course, require a relatively 

large and heavy receiver. 
e 

e 

e 

0 

e 

a 

e 

a 

e 

22 

Table X. Thermal Transport Properties of Advanced Enhancement 

PCM 

Material 

NaF 

Ge 

Fins : 

NaF/Nio. 

NaF/Ni. 

Mesh: 

NaF/Nio. 

NaF/Ni. 

NaF/Pyro. 

Matrix: 

.4NaF/.6Ge 

.2NaF/. 8Ge 

. 6Ge/Pyro. 

.8Ge/Pyro. 

.8Ge/ATJ 

Hybrid: 

NaF/Ge/ Pyro 

~ ~~ 

Density 

k~/m3 

1948.6 

5010.2 

2744.4 

2744.4 

3142.3 

3142.3 

3785.6 

4397.9 

3685.3 

4347.7 

4347.7 

3864.9 

Specific 

Heat 

kJ/kg K 

(Sol/liq) 

1.12/1.24 

.381/.396 

1.01/ 

1.01 

.964/1.06 

.967/1.07 

.680/. 730 

.530/. 560 

.640/. 652 

.640/. 652 

.go/. 91 

.499/. 516 

Thermal 

Conductivity 

kW/m K 

(Sol/liq) 

.003/.00013 

.0405/ 0.3 24 

.0087/.0069 

.0105/.0081 

.0044/.0035 

.0052/.0042 

.0045/. 0021 

.0053/.0026 

.0559/.0577 

.0340/.0365 

.0188/.0203 

.0204/.0179 

Heat of 

Fusion 

kJ/kg 

798 

481 

493 

493 

401 

401 

401 

545 

508 

392 

443 

443 

455 

Energy 

Density 

gJ/m3 

1.537 

2.401 

1.360 

1.360 

1.261 

1.261 

1.261 

2.061 

2.235 

1.446 

1.926 

1.926 

1.757 

By a process of elimination it would appear that the 

aggressive, high paybackchoicewouldbe a c;,;ibinationofGermanium 

encapsulated or contained in graphite. If high density, i.e. 

80%, packing can be achievedthis approachcombines b o t h u n u s u a l l y 

high thermal conductivity and high energy density. 

2.6 Model 

In this study effort, a number of semi-two dimensional, 

transient, finite difference code were developed to evaluate each 

of the receiver systems. The codes were then used to perform 

several pertinent parametric studies. These studies were directed 

not only as establishing the performance ofthe proposed systems, 

but also toward optimizing such performance and establishing 

directions of future development. 

e 

e 
e 23 

e 

e 

e 

e 

e 

e 

e 

While describing the model, reference will be made to the 

base line configuration. In certain of the other concepts this 

approach has beenmodifiedto represent the specific arrangements, 

but in general the overall methods remain uniform throughout. 

The. approach has been to utilize a two dimensional, transient 

analysis with an implicit formulation. All conduction is assumed 

to occur in the axial (R-Theta or X-Y, depending on the coordinate 

system selected) plane. A series of such axial elements are 

stacked between the coolant inlet and outlet. For each axial 

element the solar input is determined based on the results of 

the CAV2 output. Heat input into the cooltlnt is determined based 

on the calculated convective coefficient, coolant inlet temperature and the average wall temperature. Successive axial 

elements are connected thermally through the change in coolant 

temperature as it flows downstream. This approach utilizing the 

implicit calculational scheme provides for stable solutions for 

relatively long time steps and therefore is thought to be most 

applicable to computations involving medium to small computers. 

The two dimensional finite difference program is designed to 

describe the system shown in Figure' 11. Symmetry is assumed so 

that only one half of the element is modeled. Within the 

semi-annulus an arbitrary number of elements is used in the R 

and Theta directions. Generally this will require annular 

segimental elements where the dR segiments are uniformly spaced 

and, in general, unequal to dTheta. Within the metal containment 

representing the first and last radial elements, a somewhat 

different radial length is used corresponding to the actual 

thickness of the containment wall. 

All radiant heat transfer is assumed to occur through the 

outside circumference. In general this flux will be nonuniformly 

distributed azimuthally. In the model the total incident flux 

is constituted from a uniformly distributEd portion representing 

the indirect radiation and a frontal portion representing a direct 

radiation component. Because of symmetry both the 0 and 180 

degree walls are insulated. All convective losses are assumed 

to occur through the inside circumference. Convective coefficients are determined from conventional pipe flow correlations 

and assume a developing thermal boundary layer based on the 

distance from the inlet to the first element. 

There are, within the literature, a number of methods of 

presenting the energy equation for solution of problems involving 

phase change. The method employed here is not taken directly 

from any of these but is a logical extension of the methods 

proposed by Pantankar [8] for the general solution of conduction 

problems. Emphasis is placed on writing relatively short, simple 

programs using methods that are broadly understood within the 

heat transfer community. While it might be argued that certain 
24 

0 

e 

a 

of the advanced procedures provide for improved accuracy or 

stability, the advantages of presenting a scheme generally 

understood within the potential range of users or reviewers and 

ease of technical support are considered to be compelling. 

1 

insula t t' cl, .. 

1 t 1 .> 11 I ( I 1 c' rl 

solar- Lnfl { LX 

I i 

Figure 11. Finite Difference Element Arra:iqement. 

2.6.1 Physical model. 

The physical model is somewhat complicated by a number of 

factors including zero g effects, inadequate physical property 

data and the need for generality where various thermal 

enhancement devices might be considered. Particular caution 

must be employed in such choices in that the selection of 

differing models for alternative conditions may, in fact, 

overshadow all other parameters under consideration. This is 

particularly true in terms of the void model and while alternate 

enhancement schemes may actually justiiy differing models a 

uniform model has been incorporated here in an attempt at 

generality. 

2.6.2 Void model. 

Considerable uncertainty exists as to void location during 

the heating/cooling cycle. A characteristic of the fluoride 

salts currently being considered is that they tend to undergo 

significant density changes with phase change. Normally this 

will involve a density increase of between 10 and 25% as the 

material solidifies. Thus upon solidification voids will tend 

a 
25 e 

0 

0 

a 

0 

to form so that the question remains as to how such voids 

react. If they fail to coalesce then such voids will likely 

remain dispersed within the solid phase. If they do coalesce 

then the question remains as to how sucl; voids migrate and to 

what position. 

At this time an experiment is planned for a space flight 

test in an attempt to gain actual data on void migration. 

Prior to acquiring such data any model will be, at best, only 

speculative. In those cases where a high conductivity foam 

is to be used for enhancement of the thermal conductivity, it 

might reasonably be anticipated that migration will be 

restricted and that voids will remain within the cells where 

they are formed. It is thought possible that viscous effects 

might produce similar distributions in the more general case. 

In the absence of experimental data it is assumed that all 

voids form microspheres which are uniformly distributed within 

the solid. The remaining material will then exhibit no 

convective component to heat transfer and will retain a uniform 

effective density throughout. Note that the thermal conductivity will be affected. In this case the granular material 

model proposed by Meredith and Tobias [9 J is utilized to predict 

the effective thermal conductivity of the overall system. For 

the case in which radiation is neglected, the voids may be 

assumed to have zero thermal conductivity and the equation 

reduces to the form: 

e 

where [ = I~C VOICI tI-actlol1. 

An alternate approach is used in the case of the cavity 

receiver. In this case if voids are assumed to coalesce they 

may be assumed to migrate to the liquid-cavity interface and 

to leave a 100% dense solid phase. This is in contrast to 

those enclosed PCM systems where the final location of void 

migration remains unanswered. To remain consistent with the 

analysis of the original contractor proposing this concept it 

was decided to assume such migration in this one case. This 

approach will tend to result in higher calculated thermal 

conductivities between the melt front and the gas coolant. 

The effect will clearly be most significant when the PCM is 

largely solidified. This occurs early in the heating cycle 

and late in the cooling cycle. It should be noted that the 

time early in the heating cycle does not correspond to a time 

of temperature extremes so that the consequences on engine 

0 
a 26 

d 

a 

e 

performance and thermal stresses are not thought to be significant. However, the period late in the cooling cycle is a 

time of temperature extremes. The minimum predicted turbine 

inlet temperatures may be somewhat nonconservative. Thermal 

stresses at this same point in time are not thought to be as 

severe as those occurring late in the heating cycle so that 

these results remain uneffected. 

Computationally the void migration is accomplished by 

allowing a mass interchange between any element undergoing a 

phase change and the outermost element within the same radial 

segment. At each time step the movement of the melt front is 

calculated within the element undergoing phase change. The 

change in mass is then calculated. It is further assumed that 

the transferred mass leaves and enters the respective elements 

at the melting temperature as a method of ensuring that 

appropriate energy balances are maintained. 

2.6.3 Energy equation. 

The general transient energy equation for a ?stationary 

system may be written in terms of sensible and latent energy 

accumulation and conductive transport as: 

In finite difference form this equation may be written as: 

The major feature in descretizing this formula is in the 

treatment of the latent accumulation term. The resultant 

equation is nonlinear so that a trial and error procedure is 

required. Since it is anticipated that the temperature of any 

element with a partial phase change should be near the melting 

temperature, the velocity may be adjusted in each iteration 

so as to cause the associated node temperature to vanish. 

where 11 = a convergence factor. 
a 27 

e 

a 

In the cases reviewed here the convergence factor could be 

taken as between 0.25 and 2.0. The velocity term obtained in 

this manner is cumulative. That is, it is underevaluated in 

each iteration so that summing the value for subsequent 

iterations will eventually lead to a convergent value where 

the local temperature approaches the melting value. A second 

convergencecriteria is also includedbasedonanenergybalance 

for the overall axial segiment. In general the goal is to 

seek energy convergence to within .1%. 

2.6.4 Computer Program. 

Selectionof an implicit scheme tosolve the nodal equations 

is designed to enhance numerical stability and to permit the 

use of reasonably long time increments. This is viewed as a 

necessity for such programs when utilizing relatively small 

computational facilities. Amajor advantage to such an approach 

is that the resultant program ishighlyportable, being designed 

primarily for PC usage. 

In order to simplify the programming and to reduce the 

number of parameters being stored a hybrid computational scheme 

is employedwhichincluded features of bothGaussian elimination 

and Gauss Seidel. Considering the figure of the finite 

difference scheme, the approach is to successively solve for 

the temperatures along each of the vertical columns, using a 

one dimensional (tridiagonal) method. While evaluating the 

temperatures along any one column, temperatures along any other 

column are assumed to be known based on the last iteration. 

Utilizing successive sweeps, this method provides for a rapid 

convergence while limiting memory requirements. 
0 28 

* 

e 

0 

a 

3 Parametric Study 

In the parametric studies which follow, the effects of the basic 

geometric parameters influencing overall thermal performance of the 

canisters are examined. These include the number of canisters, 

element geometry or spacing and the use of alternative thermal 

enhancement materials. During the transient, temperatures are 

evaluated at two times, 54 and 90 minutes. These times correspond 

to the end of the solar incidence, where maximum temperatures are 

achieved, and to the end of the solar eclipse, where minimum 

temperatures are encountered. The implications of suchtemperatures, 

while largely qualitative, directly influence both receiver life 

and cycle efficiency. 

One important aspect ofthe study is to underline the differences 

in thermal performance between metals and salts as PCM's. Because 

of the inherently low thermal conductivity of the salts, their 

parametric range has been extended with a view of establishing 

potential for improved performance with improved geometry. In the 

case of metals this same extension ha3 appeared unnecessary. 

3.1 Base Line Receiver. 

The base line receiver utilizes convective transport directly 

to the working fluid rather than utilizing a heat pipe as the 

otherdesigns consideredhere. Thenatureofthedesignconstrains 

the parametric study to some extent so that certain combinations 

of enhancement and configuration perturbation are not possible. 

This is notalimitation ofthedesign, but rather oftheparameters 

selected. In order to take advantage of the full range of 

enhancement the system would require a different inlet temperature 

and/or flow rate. The values chosen for the study basically 

correspond to those initially proposed in the design. This 

decision has meant that considerations OL enhancement are thus 

addressed on a more limited scale. 

It should also be noted that any enhancement effect due to 

the canister side walls is neglected. If foam enhancement is 

added it is assumed that the foam will provide for the void 

management so that a large number of individual canisters is not 

required. Since exclusion of the canister walls is thought 

necessary for evaluation of enhancement effects it has also been 

omitted fromthebasecaseto fully indicate relative improvements. 

The turbine inlet temperature for the base line case during 

normal orbit is shown in Figure 12. Thermcihynamically the outlet 

temperature must average 1000 K (LiF 22 CaF2 melts at 1039 K) 

and variations around this value are due to alternate melting 

and sensible heating or solidifying and sensible cooling of the 

PCM. The temperature is seen to increase during the solar 

a 
* 29 

incidence reaching a value about 30 degrees above the average 

and cooling rapidly during eclipse to about 30 degrees below the 

0 average. 

e 

0 

0 

Figure 12. Base Line Orbital Turbine Inlet Temperature. 

The effects of thermal enhancement are seen in Figure 13. 

Note that the plotted values do not extend to include the 18% 

enhancement value. A projection of the partial low temperature 

lines would indicate that this case would fall above the average 

gas temperature violating the first law. However it is seen that overthe effective range the lower temperature is increased 

an average of 3.4O and 4.2O for each percent addition of BN and 

pyrolytic graphite, respectively and that the enhanced materials 

can be utilized effectively to control temperature swings. Full 

utilization of the enhancement would require either an increase 

in coolant flow rate or reoptimization of the inlet temperature. 

Increasing the number of elements serves to reduce the radius 

of the individual canisters and subsequently reduces the thermal 

resistance within each. As indicated in F;.gure 14 the effect is 

fairly significant up to between 60 and 75 elements producing 

about a 3O improvement per element. In view of the weight 

advantage of ‘reducing the number of elements it would appear that 

the optimal condition would occur at a number near 60 or 65 

elements. 
30 

e 

Ease Line Receiver 

P/D 166 

L/D I 75 

-'JO~T, 

I 

, 4 . , b , , a , , , , ~ , , , IFet :.elit Enhilrir:ement Z 

Figure 13. Effects of PCM Enhancement on Turbine Inlet 

Temperature. 

c 

1 L 

Figure 14. Effects of Element Number on Turbine 

Temperature. 

Inlet 

The effects of increasing the receiver L/D ratio are shown 

in Figure 15. Again increasing the length of each element tends 

to reduce the radius of individual canisters again reducing the 

effective thermal resistance. The magnitude of improvement is 

only about 1.3O per percent change in L/D. While the effect 

shown is relatively small it may, in fact, be something of an 

overstatement. In the cases demonstrated here the' receiver flux 

a 
31 

a 

a 

e 

e 

distribution was assumed to remain constant with each of the 

configurations. In fact as the L/D increases the heat flux will 

tend to become somewhat more maldistributed. In this case the 

improvements seen at high L/D ratios will be somewhat reduced. 

I 

Figure 15. Effects of Receiver L/D on Turbine Inlet Temperature. 

The effect of increasing the element circumferential pitch 

to diameter ratio would be to decrease the azimuthal variation 

in the incident heat flux. The improvement in the utilization 

of the PCM will then lead to an improvement in the uniformity of 

the inlet turbine gas temperature. This effect is seen in Figure 

16. It is noted that the reduction in variation appears nearly 

linear over the P/D range of 1.25 to 2.0 with a corresponding 

improvement of 30°. 

In Figures 17 through 19 the maximum and minimum canister 

temperatures are shown at 54 minutes in the first section of the 

element. These values are obtained from the same runs in which 

the turbine inlet temperatures were ca1c::lated. What will be 

generally noted is that the PCM temperatures show much the same 

trends but are much more sensitive to the geometrical parameters. 

Where poor PCM utilization is encountered there is a shift from 

storing energy in its latent form to storage as sensible heat. 

The inevitable results include large variations in turbine inlet 

temperature accompanied by high peak PCM temperatures. 

The PCM temperature differential as a result of increasing 

the element number is shown in Figure 17. The design is subjected 

to considerable peaking of the incident heat flux and shows 

considerable sensible heating until the number of elements nears 
0 ORIGINAL PAGE 9s 

OF POOR QUALITY 

32 

0 

90. The temperature difference across the canister is also 

decreased with increasing number of elements. The improvement 

is somewhat compressed by the scale of the figure. In this case 

the temperature differential decreases from about 60° to 25O as 

the number of elements is increased from 45 to 90. Note that 

improvement in canister temperature differences is over 50%. 

a 

a 

e 

c 

Figure 16. Effects of Element P/D on Turbine Inlet 

ture. 

TemperaFigure 17. Effects of Element Number on Canister Temperature 

Difference. 
0 33 

a 

0 

0 

0 

0 

The variation in PCM temperature with receiver length to 

diameter ratio is seen in Figure 18. The improved utilization 

of PCM serves to decrease the peak temperature from about 305O 

to about llOo above melting as the L/D increases from 1 to 2. 

The canister temperature difference is seen to decrease from 

about 40° to nearly 20° over the same range. 

320 

Base Line Receiver 

Lip 22 CaF 

P/D 175 

r 

1 on 4 

r 

?ILL11 

Figure 18. Effects of Receiver L/D on Canister Temperature 

Difference. 

Figure 19. Effects of Element P/D on Canister Temperature 

Difference. 
34 

e 

The effect of element pitch to diameter ratio on PCM temperatures is shown for the base line receiver in Figure 19. Peak 

temperatures decrease from about 270°to about125O above melting 

as the P/D increases from 1.25 to 2. Canister temperature 

differences decrease from about 40° to about 20° over the same 

range. 

3.2 Heat Pipe Receiver. 

I 

I 

Figure 20. 

Temperature. 

In Figure 20 the heat pipe vapor temperature is indicated 

through the steady orbit. It can be seen that temperatures change 

markedly between eclipse and solar incidence. The heat pipe 

fluid must operate at above melting to store energy during the 

incidenceperiodandbelowmeltingtoremoveenergyduringeclipse. 

As either period progresses the melt front advances through the 

PCM so that a continually increasing temperature difference is 

required to sustain the heat transfer rate. One can also see 

the effects of thermal enhancement of the PCM. Additions of 18% 

by volume of BN or ATJ graphite substantially reduce the indicated 

temperature variations with overall variations decreasing from 

about llOo to about 75O and 60° respectively. The use of a metal 

PCM is seen to virtually eliminate all such variations providing 

only about a 25O change. Note that virtually all of the variation 

in the last case occurs in the very last time step. This is an 

indication of excessive solidification and use of only sensible 

energy. With some adjustment of the operating conditions this 

drop can be eliminated reducing the overall variation to only 

about 5O. 

Effects of PCM Selection on Heat Pipe Vapor 

a 
35 

a 

a 

Q 

The maximum and minimum heat pipe vapor temperatures corresponding to the end of solar incidence (54 min.) and the end 

of solar eclipse (90 min.) as a function of the percent of thermal 

enhancement are shown in Figure 21. Also shown are the corresponding swings encountered with the metal PCM (Germanium). 

Enhancement materials serve to improve the radial conduction but 

theeffectispartiallycounteredbydisplacingPCM. This produces 

an improved conductivity but requires that heat be conducted over 

a larger distance. While the effect is significant in either of 

the enhancement materials, both are substantially less effective 

than the use of metals. From the heat pipe vapor temperatures 

produced itwouldappearthebetween 8 and15% of either enhancement 

material is required to produce a significant improvement. 

Figure 21. Effect of Enhancement Fraction on Heat Pipe Vapor 

Temperature. 

An increased number of elements has the effect, as seen in 

Figure 22, of decreasing the radial thickness of PCM in each 

canister. The result should be a decrease in heat pipe vapor 

temperature variation as the thermal resistance between the heat 

pipe fluid and the melt front is decreased. With LiF 22 CaF2 

the variation in heat pipe vapor temperatures is seen to be 

substantially reduced as the number of elements is increased to 

between about 20 and 25. Beyond this point the impact of further 

additions appears limited. The effect can be achieved with far 

fewer elements with some degree of thermal enhancement. The 

metal PCM is shown only for the case of it, elements. It would 

be expected that the effect of increasing the canister number 

for the Ge PCM would result in a similar curve for the salts. 
36 

a 

Figure 22. Effect of Element Number on Heat Pipe Vapor 

Temperatures. 

Q 

Figure 23. Effect of Receiver L/D Ratio on Heat Pipe Vapor 

Temperatures. 

As the receiver is made longer the radial dimension of the 

individual canisters is decreased. Thermally this results in a 

decrease in both the heat flux and the radial resistance through 

which energy must flow. As seen in Figure 23 the effect appears 

to be relatively large for ratios up to 1. For longer receivers 
ORlGENAL PAGE IS 

OF POOR QUALITY 

37 

a 

a 

e 

a 

there appears to be only a small improvement in the variation. 

Again the effect of the solar flux profile would be to moderate 

the improvement at the longer L/D ratios. 

Both the absolute PCM temperature and the temperature difference across the canister are seen to decrease as the number 

of elements is increased. Results shown in Figure 24 indicate 

that the LiF 22 CaF2, having a relatively lowthermalconductivity, 

will require considerable sensible energy storage to ensure 

adequate potential to transfer heat to the heat pipe vapor during 

the eclipse phase. The increase in the number of canisters from 

15 to 30 serves to drop the peak temperature from about 120 above 

melting to about 30°. This is seen to occ?!r in an exponentially 

decaying fashion. The minimum temperature is not shown here but 

will, by design, be quite close to melting. Thus the peak 

temperature above melting is essentially the same as the temperature difference. 

T------ 

'-0 * r 

1 

, 

I 

i 

1 1 

- ";y i I 

: 51' I-. 

c" 2 =u 1 

ir 

.. '. 

, e-- J 0 BPI '. 

. -- I 

'. 

I ,Jplllt? ...- 

I -7 -IT--, I , 7T--r4 

? ili c I !

.- I 

-. ._ 

1-1 , ~ 

-- - -----%--- -_. I 

I1 - --- - .- --. 

- , 'j 

j --- 77---- 

1- II -1 -J 

1c - - 

Piumbrr t,t tlem-nts 

Figure 24. Effects of Element Number on Canister Temperature 

Difference. 

As seen in Figure 25 the effect of increasing the receiver 

length to diameter ratio is also to decrease both PCM temperature 

and temperature difference. In this case the temperature difference for LiF 22 CaF2 is decreased from about 100 to about 40 

degrees as the L/D ranges from 0.66 to 1.66. 

a 
ORIG!F:AL PAGE IS 

OF POOR QUALtTY 

3% 

LIF 22 cot ! 

Figure 25. Effects of Receiver L/D on Canister Temperature 

Difference. 

e 

Heat Pipe Receive 

L,’D 1 07 

Ei-?-rnentc, LO i !,., 

i.8 

Figure 26. Effects of PCM Enrichment on Canister Temperature 

Difference. 

The PCM maximum and minimum temperatures are plotted against 

the element enrichment in Figure 26. As noted earlier, the effect 

of enrichment is to improve thermal conductivity butthis occurs 

at the cost of increasing the length of the path through which 

heat conduction must occur. The utilization of 18% Boron Nitride 

or Graphite is to decrease the peak Element temperature from 

about 75 degrees above melting to around 50 and 30 respectively. 
0RiGCNAL PAGE IS 

OF POOR QUALITY 

39 

e 

a 

a 

e 

0 

In this case the minimum temperature remains, by design, 

essentially at the melt temperature. The element temperature 

differences will, therefore, have the same value. 

3.3 Pebble Bed Receiver. 

Themodeling ofthe pebblebed assumedthatthe pebble stacking 

would provide for efficient gas mixing so that the working fluid 

provides for a degree of PCM enhancement. This effect reduces 

the radial thermal resistance in the elements so that the impact 

of the maldistribution of input heat flux associated with low 

element circumferential pitch to diameter ratios is significantly 

reduced. As seen in Figure 27 the effect of varying the P/D 

between 1.4 and 2 is negligible on the minimum turbine inlet 

temperature, typically about 990 K throughout the range. The 

maximum turbine inlet temperature, corresponding to the end of 

solar incidence is seen to vary only between about 1050 down to 

1070 K over this same range. The effect of graphite enhancement 

is seen to be small but typically reduces the maximum temperature 

by 20 to 30 degrees over this same range. 

Figure 27. Effects of Element P/D on Turbine Inlet 

ture. 

TemperaAs seen in Figure 28, the effect of varying the receiver 

length to diameter ratio is also small. Varying its value from 

0.7 to 1.5 appears to decrease the turbine inlet temperature by 

about 15 degrees. This same variation occurs at both the peak 

and the minimum temperatures. 
40 

Figure 28. Effects of Receiver L/D on Turbine Inlet Temperature. 

Figure 29. Effects of Element P/D on Canister Temperature 

Difference. 

The effect of element pitch to diameter ratio on the PCM 

temperature is shown in Figure 29. The effect is somewhat greater 

than that seen on the turbine inlet temperature. Here the P/D 

is varied from 1.4 to 2 the maximum PCM temperature at 54 minutes 

is seen to decrease from about 1122 to 1095 K. At the same time 

the minimum PCM temperature increases from about 1017 to 1028 K. 

The overall variation across the element is then seen to decrease 
0 

e 

a 

0 

0 

41 

by about 25%. Both effects would be thought to be beneficial 

toward reducing high temperature, stress related structural 

problems. 

The PCM maximum temperature and temperature difference at 54 

minutes is also seen to be a weak function of the receiver length 

to diameter ratio as seen in Figure 30. As L/D varies from 0.75 

to 1.6 the peak PCM temperature is decreased from 1110 to 1100 

K. The temperature difference over this same range decreases 

from about 95 to 70 degrees. 

-1 Pebble 8erl Rrceivrrl 

Figure 30. Effects of Receiver L/D on Canister Temperature 

Difference. 

3.4 Cavity Heat Pipe Receiver. 

As noted earlier each canister has associated with it a heat 

pipe wedge which serves effectively as a fin for thermal 

enhancement. Clearly the effect of additional fins will be to 

reduce temperature differentials within the element. Because 

the distances through which the energy is being driven is reduced 

with increased canisters, the temperature difference between the 

heat pipe fluid and the PCM melt temperature is reduced. The 

net effect is that cavity temperature variations between the end 

of the periods of solar incidence and solar eclipse are reduced. 

The requirement for small high temperature source variations is 

thought to be more severe for the Stirling engine than for the 

Brayton; it is noted that a Stirling receiver would therefore 

normally incorporate a higher number of elements than a corresponding Brayton receiver. 
ORIGINAL PAGE E 

OF POOR QUALITY 

42 

e 

a 

0 

* 

e 

e 

The results of the parametric study indicating the effect of 

number of canisters is shown in the Figure 31. Results indicate 

that temperature reductions average about 3.5 K per canister for 

the LiF 22 CaF2 without enhancement and the BN enhanced PCM. The 

pyrolytic graphite demonstrates a somewhat reduced effect with 

only about a 2.5 K reduction per canister. For the 100% PCM case 

this would imply that doubling the number of canisters from 18 

to 36 would decrease the cavity temperature variation from about 

175 to about 115 K. While this trend is somewhat reduced with 

the enhanced materials, 18% pyrolytic graphite will still undergo 

an improvement of from about 125 to about 65 K. The Ge, having 

an extremely high thermal diffusivity, shows a quite low temperature variation. In spite of the low value overall trends 

appear the same as for the remaining materials. Clearly the 

number of canisters is seen to be a primary factor in determining 

the magnitude of cavity temperature variation. 

Figure 31. Effects of Element Number on Canister Temperature 

Difference. 

A small apex angle will also tend to reduce temperature 

variations by reducing the maximum distance that heat must be 

transported into the solid. So long as a portion of the PCM 

remains partially melted the entire system is then coupled to 

the melting temperature. Short thermal distances are then seen 

to be close coupled whereas longer distances permit greater 

thermal swings. It is further noted that the small apex angles 

which tend to enhance heat transfer, involve configurations with 

large surface to volume ratios so that system weight will tend 

to suffer. Clearly additional trade-offs are required to optimize 

the system in this regard. 
43 

e 

a 

0 

0 

The effect of apex angle on cavity temperature oscillations 

is shown in the Figure 32. The relationship appears to be nearly 

linear over the range shown, about 3 K per degree. Earlier 

results, not shown here, indicate that the enhancement effects 

begin to decrease above angles of about 22 degrees so that cavity 

temperature swings approach an asymptotic value. In effect 

decreasing the apex angle from the asymptotic limit to about 15 

degrees will only decrease temperature variations from about 125 

to a little over 100 K. A further reduction to a very sharp 

angle of 10 degrees will only drop the variation to 90 K. It is 

concluded that apex angle provides only iL;odest improvement in 

cavity temperature variation and generally will be of secondary 

importance. 

Figure 32. Effects of Element Apex Angle on Cavity Temperature Difference. 

Increasing element length will have the effect of decreasing 

the other element dimensions, specifically the depth and width. 

Since the latter two dimensions are the most significant in terms 

of the heat transfer problem, the net effect will be to decrease 

cavity temperature variations. 

As seen in Figure 33 for the unenhanced PCM, increasing the 

receiver length to diameter ratio from 0.5 to 2.0 produces a 

corresponding decrease in cavity temperature variations from 

about 160 to 80 K. For the two foam enhanced schemes the effect 

is similar; temperature variations are about halved. As would 

be expected the temperature variations for the metals, having 

very high thermal diffusivities, are quite small. As with the 
44 

0 

0 

0 

0 

0 

0 

salts, the overall trend is toward decreased variations with 

increased L/D. In all cases the cavity temperature variation is 

seen to be highly dependent on receiver L/D. 

Figure 33. Effects of Receiver L/D on Cavity Temperature 

Difference. 

The effects of thermal enhancement are shown in Figure 34. 

Theeffectofenhancementtendstobenonlinearand small additions 

tend to produce limited results. This tendency is particularly 

seen for the BN case where additions of less than 10% show no 

significant overall improvement in therm;:: performance. A 20% 

addition has the overall effect of reducing cavity temperature 

swings from about 115 to about 70 K. The pyrolytic graphite 

would appear to offer more substantial thermal improvement; cavity 

temperatures are reduced to 90 and 50 K for foam additions of 10 

and 20%, respectively. Thermal enhancement is thereby seen to 

provide an effective means of further temperature control. 
'* 45 

0 

0 

0 

0 

a 

a 

a 

e 

+- 

-2u 1 Pvrolvtic 

Figure 34. Effects of PCM Enhancement on Cavity Temperature 

Difference. 

3.5 Cavity Receiver. 

The variation in the turbine inlet temperature associated 

with the cavity receiver is shown in Figure 35. The range is 

seen to be relatively low with a minimum of 975 and a maximum of 

1015. This is the least variability of any of the designs 

considered. The effect is thought to be closely related to the 

large number of corrugations considered. Each of the other 

designs showed a tendency toward decreasinr, variations with large 

number of elements. One advantage of this concept is that it 

seems to make the utilization of large numbers of PCM subdivisions 

without the weight and mechanical complexity of certain of the 

other designs. 

shadow the effects of other parameters on the turbine inlet 

temperature. Results showed little discernable variation in 

these temperatures with L/D, P/D enhancement, etc. Again, if 

the unit were designed with a number of corrugations more consistent with the number of elements in other receiver designs 

it is thought that this receiver would follow the established 

patterns. 

ters described the associated drawings have been omitted. 

The effect of the large number of corrugations is to overBecause of the lack of variability in the paramea 
e 46 

a 

0 

a 

a 

0 

0 

0 

I) 

Figure 35. Turbine Inlet Temperatures Through Orbit. 

The effects of PCM enhancement on the PCM maximum and minimum 

temperatures are seen in Figure 36. Effects are seen to follow 

a pattern similar to that of other receivers. Clearly there is 

a tendency to reduce peak temperatures as enhancement is added; 

minimum temperatures appear to remain dirtually unchanged. 

However the effect is seen to be relatively small. The receiver 

is thought to be designed with sufficient corrugations that the 

dimensions on individual PCM ltelementstt are small. The peak 

temperature is seen to decrease by about 20 degrees with 18% BN, 

about 35 degrees with a similar amount of graphite addition. 

Improvements in the temperature differentials are also relatively 

small, on the order of 20 to 30%. It should be recognized that 

with the transparency of the PCM much of the thermal energy 

bypasses the PCM reducing the required differential to produce 

melting in the PCM. 

A unique feature of this design is the ability to vary the 

energy distribution in the PCM by varying the absorptivity. This 

is done by controlling the quantity of absorbing material added. 

It would be anticipated that an optimal value of absorptivity 

could be found. A large value would tend to cause excessive 

heating at the PCM surface and result in increased temperature 

differentials. A low value would produce excessive absorption 

at themetalcorrugation-PCM interface. Under an extreme condition 

virtually all energy would be added to the gas during periods of 

solar incidence with little storage occurring. Under these 

conditions no temperature difference woullt occur across the PCM. 

While the results shown in Figure 37 suggest that in the range 

indicated, the lower the absorptivity the better, this result is 

0 
0 

2 toso- EN Cr I 

ORlGlMAL PAGE 1s 

OF POOR QtlliEIfY 

47 

misleading. Preliminary results suggest that a value of between 

0.8 and 0.95 would provide a suitable balance between these two 

considerations. 

:r- L,/D I I \A -- Corruqations I :5~,,,, 

',; Enhancement 

Figure 37. Effects of PCM Absorptivity on Canister 

ture Difference. 

TemperaThe effect of increasing the element length to diameter ratio 

on the PCM temperature and temperature difference are shown in 

Figure 38. The result may initially appear counter to those 

developed for other concepts in that temperatures and temperature 

0 
e 

48 

e 

e 

e 

e 

e 

e 

e 

e 

e 

differences increase with increasing L/D This is an artifact 

of the method used in varying the L/D in the parameter study. 

The decision was made to hold the minimum receiver diameter 

constant while allowing the larger diameter to vary with length 

(The corrugation angles were held fixed). The net result was to 

increase the PCM radial dimensions as length was added. 

It may be observed that the minimum PCM temperature is held 

slightly below the melting point at all times so that this 

parameter shows little variability. The maximum temperature 

increases relatively rapidly as the radial PCM dimension 

increases. (Note that the added fraction of absorbing material 

must be continually adjusted to maintain a cunstant absorptivity. ) 

Under these conditions the temperature differential is seen to 

nearly double, from 80 to 145 degrees as the L/D is increased 

from 0.66 to 1.33. 

Figure 38. Effects of Receiver L/D on Canister Temperature 

Difference. 

As seen in Figure 39 the effect of increasing the number of 

corrugations is similar to that of increasing the number of 

elements for other concepts. Here the range extends from 75 to 

150. Under these conditions the minimum temperature shows little 

variation, remaining slightly below the melting point. The 

maximum temperature decreases from about 1165 to about 1130 for 

the unenhanced PCM. Additions of 18% BN decreased this peak 

temperature from 1135 to 1110 over this same range. Graphite 

provides for a further reduction ranging from 1110 to 1090 K. 

From the increasing slope of the curve it would appear that 

e 
49 e 

a 

e 

OWICllNAL PAGE 13 

OF POOR QUALITY 

significant reductions in the number of corrugations below the 

indicated range would lead to rather large increases in these 

values. 

----7--- , It1 I I) +--- I 

YU 110 150 I50 70 

rlurriner ot Corrugations 

Figure 39. Effects of Element P/D on Canister Temperature 

Difference. 

e 
e 50 

4 Discussion 

It should be emphasized that in a evaluation of this type, 

includingonlythermalperformance ofthe relative concepts, relative 

comparisons are somewhat limited and should not be used as a sole 

basis for receiver selection. Moreover there remain a certain 

arbitrariness in establishing certain of the configurations. 

Specifically those' designs including canisters were not credited 

with the thermal enhancement properly attributed to that design. 

This was done in order to illustrate a continuity in the effects 

of thermal enhancement where it was assumed that PCM distribution 

would be accommodated by the foam filler itself. Further it is not 

clear to the authors that the effects of potential liquid PCM 

circulation and thermal insulation due to possible void coalescence 

and migration have been equitably evaluated between competing 

concepts. Finally the use of the CAV2 generated incident heat 

fluxes may unfairly penalize the non heat pipe concepts by failing 

to account for radiant redistribution of the incident heat flux 

within the receiver. Nevertheless it is believed that the results 

presented are useful in a number of ways. In a general sense they 

may serve to provide reasonable estimates of the performance 

achievable. They also provide an indication of the benefits 

attributable to various enhancement alternatives. Finally they 

serve to indicate that range of individual parameters which tend 

to optimize thermal performance within each of the respective 

designs. 

e 

e 

e 

e 

4.1 Attainable Performance 

Attainable thermal performance is evaluated here in two 

distinct ways. The first involves the variation in receiver 

outlet temperature, generally termed either turbine inlet temperature, heat pipe temperature, or engine head temperature. The 

second involves PCM peak temperatures or PCM thermal differences. 

These parameters relate to the thermal stress and thermal cycling 

within the PCM containment and have a direct bearing on the 

ultimate reliability of the receiver. 

The cyclic thermal variation for the base line receiver is 

shown in Figure 12 where the overall variation is seen to be of 

the order of 60 K. The basic problem of optimizing performance 

with this design is the maldistribution cf incident energy in 

e the axial direction. Excess incident energy is provided at the 

aperture end ofthe receiver tending to provide excessive sensible 

heating of the PCM during the latter stages of the heating cycle. 

Insufficient heating occurs at the opposite end of the receiver 

so that complete solidification and sensible heat removal during 

the eclipse phase remains a problem. The design problem then 

becomes one of making a trade off between excess temperature near 
0 

51 

0 

a 

e 

a 

e 

e 

the receiver inlet versus increased engine hot side temperature 

variations. Themost direct solution ofthe problem is to increase 

the amount of PCM in the receiver. 

The corresponding plot for the heat pipe receiver is shown 

in Figure 20. It is important to note here that the number of 

elements in the heat pipe system is reduced to 20 from the original 

75 used in the base line receiver. It is also important to note 

that the calculated results for the heat pipe system do not 

include any enhancement effects from the canister side walls. 

As a consequence the LiF 22 CaF2 results indicate somewhat larger 

temperature swings than would be observed in such a system. Based 

on prior runs a reasonable approximation for such a system would 

be to assume that the effects of adding 12% metal fins would be 

to produce about 2/3 of the improvement of adding 18% BN. Thus 

the anticipated temperature variation for the system is thought 

to be about 94 K. Increasing the number of elements from 20 to 

30 with the heat pipe system, as seen in Figure 22, would serve 

to decrease the temperature difference to about 50 K. The two 

systems are strikingly similar in terms of individual components 

and would seem to offer the most viable direct comparison of any 

two receivers under Consideration. The improvements associated 

withtheheatpipeconfigurationmaybeattributedto a combination 

of factors. First the heat pipe serves to smooth the incident 

heat flux. Second, energy utilized by the engine during the 

period of solar incidence is not forced through the PCM, but 

rather is transported directly from the aperture section of the 

heat pipe directly to the engine. The resultant reduction in 

heat flux through the PCM can only serve to reduce the required 

temperature differential. 

The engine hot side temperature varidtion for the cavity 

receiver is shown in Figure 35. This design also utilized a 

direct solar incidence on the PCM with the working fluid flowing 

past the PCM. As with the base case this design has the problem 

of excessive heat addition near the aperture and insufficient 

heat flux at the outlet end. It does benefit from a somewhat 

higher number of canisters than the base line receiver and from 

a truncated cone shape which redistributes the incident heat flux 

to an extent. The result of these two improvements is seen in 

reducing the working fluid temperature variation to about 47 K. 

It is concluded that an overall terrLerature variation of 

between 50 and 100 K is reasonably achievable with an unenhanced 

LiF 22 CaF2 PCM receiver. As will be noted in the following 

sections, the lower range of values will typically be associated 

with arrangements optimized for thermal performance rather than 

weight. 
e 

52 

e 

0 

0 

e 

d 

The range of radial temperature variation across the individual elements is also of importance, not to engine performance, 

but rather as an indication of the thermal stresses established 

within the PCM containment. This variation is normally a maximum 

at the end of solar incidence and is therefore evaluated at that 

time. Shown in Figure 19 is the maximum and minimum temperatures 

within that range of canisters nearest the aperture of the base 

line case. Since the reference condition is at a P/D of about 

1.66 it can be seen that the maximum and minimum temperatures 

correspond to about 172 and 148 degrees above melting or 1211 

and 1187 respectively. Note that thevariation across the element 

is only 24 degrees but the absolute values are relatively large. 

This is due to the high axial incident flux near the aperture. 

The cavity receiver, as shown in Figure 39, tends to have 

slightly improved PCM temperature differentials over the base 

line case. With 125 corrugations the PCM maximum is 1132 and 

minimum is 1022 K. The reduction in the maximum temperature can 

be attributed partially tothe improved incident flux distribution 

associated with the truncated cone cdnfiguration. It is also 

important to note that the coolant passages in this design are 

suitable for conforming to axial needs. in this analysis it has 

been assumed that the convection coefficient is axially adjusted 

to provide an axial NTU distribution which matches the fraction 

of heat removal tothe fraction of incidentenergy. Upon reviewing 

the final results it was realized that this effect could be 

carried to a greater extent. That is the matching of the rate 

of heat removal to the incident heat flux leads to a situation 

in which the orbital temperature variation tends to produce axial 

PCM temperature patterns in which the minimum temperature appears 

as something of a mirror image of the maximum temperature with 

an axis about the melting point. A further refinement of the 

method would increase the rate of heat removal at axial locations 

corresponding to the highest PCM temperatures. This would distort 

the mirror image decreasing the peak temperatures and producing 

significant sensible heat removal during the heat removal phase 

at these same locations. 

In comparison both the heat pipe and the cavity heat pipe 

receivers operate at maximum and minimum PCM temperatures of 

about 60 and 0 degrees above melting or 1099 and 1039 K, 

respectively. The magnitude of the temperatures is significantly 

reduced due to smoothing of the axial flux in the heat pipe. The 

temperature differential is, however, considerably larger. This 

may be attributed to the reduced number of elements and the 

subsequent increase in individual canister radius. 

a 

e 
53 

a 

a 

0 

a 

a 

a 

a 

The pebble bed arrangement is by it nature somewhat different 

in terms of element temperature variations. The values shown in 

Figure 30 indicate that the temperature difference across the 

PCM mass ranges from 1025 to 1110 K. Note that this difference 

is not across a single sphere but represents the temperature 

difference between spheres at the center and edge of the bed. 

Forthis reasoncomparisons of metal stress between these elements 

and the alternative designs are difficult. In this design the 

PCM sphere size is essentially 1/3 of the element containment 

ID. A simple analysis would then suggestatemperature difference 

across each sphere of about 30 to 35 degrees. 

4.2 PCM Thermal Enhancement 

The effects of thermal enhancement on the cycle hot side 

temperature are shown for the base line receiver in Figure 13. 

For the heat pipe receiver these same results are shown in Figure 

20. Neither case provides a clear comparison of the various 

enhancement schemes. As noted earlier, the set of conditions 

chosen for the base liner receiver do not lend themselves to the 

use of Ge or 18% additions of BN or pyrolytic graphite. As a 

consequence consideration has been limited to 10% additions of 

the latter two enhancement materials. Moi-eover both cases are 

somewhat complicated by the question of whether and how side wall 

should be considered. The approach proposed here is simply to 

compare the results neglecting the side walls but to weight 

certain of the other cases more heavily in making the final 

comparison. 

In the base line case the cyclic variations are reduced from 

about 77 K to 43 and 37 K for a 10 % addition of BN and graphite, 

respectively. This corresponds to a 44% and a 52% reduction in 

variation from the base line. The enhancement of the heat pipe 

receiver is seen to reduce the variation from about 110 K to 

about 75, 60 and 5 K for the 18% BN 18% graphite and the Ge, 

respectively. This results in a 32% improvement for the BN, and 

a 45% improvement for the graphite. Note that while the while 

the enhancement fraction is much lower for the base line case 

the improvement achieved is greater. The introduction of Ge as 

a phase change material reduces the overall temperature variation 

to near zero and would clearly be the most advantageous choice 

from thermal considerations. 

A somewhat simpler comparison is made in the various 

enhancement schemes in the cavity heat pipe receiver. Referring 

to Figure 31, it is seen that using the Leference 36 elements 

the introduction of thermal enhancement reduces the variation 

from 115 to 90 and 62 K, respectively. The improvement is seen 

to be a reduction of 22% and 45% respectively. While the Ge PCM 
a 54 

a 

0 

0 

0 

a 

a 

has not been evaluated for this same number of elements, certainly 

the variation will be less than for the 2.1 element case or less 

than 15 K. The corresponding reduction is then greater than 87%. 

The cavity receiver provides a second relatively direct 

comparison. As seen in Figure 34, the BN serves to reduce the 

variation from about 117 to 95 K, or about 19%. Similarly the 

graphite reduces the variation to about 70 K or about 40%. 

Summarizing the cases above it is concluded that thermal 

enhancement with BN would normally be expected to decrease engine 

hot side temperature variations to between 20 and 25%, graphite 

to between 40 and 45% and Ge to between 85 and 95%. 

Canister temperature differences are shown for the cavity 

heat pipe receiver can be taken from Figure 32 recognizing that 

at 54 minutes the minimum temperature will be at the back wall 

of the element where only a residual solid mass of PCM remains. 

As a consequence the element temperature difference will be 

between the melting temperature and the peak temperature. It 

would appear then that the addition of 18% BN would serve to 

reduce temperature variation across the element from 65 to 48 K 

or by 26%. Graphite enhancement reduces the differential to 32 

K or by 51%. Germanium PCM would encounter only about a 5 K 

difference at this same time representing a reduction of 92%. 

The cavity receiver indicates a much larger temperature 

difference across the elements as indicated in Figure 36. Here 

the unenhanced material operates with a peak temperature difference of about 120 K. Eighteen percent additions of BN reduce 

this to about 100 K or about 17% and similar additions of graphite 

to about 79 K or by 34%. 

Summarizing these conclusions it appears that BN enhancement 

would potentially reduce peak PCM temperature differences by 

between 15 and 30%. Similar additions of graphite would reduce 

these same differences by between 35 and 50%. Germanium PCM 

would virtually eliminate such differences, reducing them by 

about 90% as compared to the unenhanced salt. 

4.3 Parametric Studies 

In reviewing the thermal performance of each of the design 

concepts, it becomes clear that most of the corresponding performance curves appear similar in general form. That is the 

cycle hot side temperature, whether it be termed the turbine 

inlet temperature or heat pipe temperature, tends to become more 

uniform as the effective thickness of the PCM is decreased. 

Several methods remain for controlling t!id effective thickness 

in each of these designs. 
55 

e 

0 

Generally, increasing the receiver length will tend to 

redistribute PCM so as to provide a smaller radial path for energy 

storage. In the study shown here this cifect is shown as the 

receiver length to diameter (L/D) ratio. An apparent exception 

to decreased temperature variations with increased L/D is the 

case of the cavity receiver. This case involves the use of an 

absorbing phase in the PCM. It is thought that the higher 

temperatures are associated with the steeper energy gradient 

associated with the increased concentration of absorbing material 

utilized to maintain a constant average absorptivity. 

In each of the designs it was found that increasing the number 

of elements, or corrugations in the case of the cavity receiver, 

tended to reduce variations in the cycle hot side temperature. 

This again is attributable to the reduction in the effective 

thickness of PCM through which the heat flux must travel. In 

those cases which did not utilize the heat pipe concept to 

redistribute the incident heat flux, it was found that increasing 

the element pitch to diameter ratio would tend to redistribute 

incident heat flux in a more uniform method also reducing temperature differences. 

a 

a 
a 56 

5 Conclusions. 

While the base line concept represents a relatively high degree 

of development in thermal receiver design, substantial improvements 

in thermal performance may still be achievec,. Generally these are 

accomplished through'one of several methods including (1) the use 

of heat pipes to redistribute incident heat flux in a uniform manner, 

(2) better utilization of the PCM material with less sensible heating 

and cooling, (3) introduction of concepts which completely or 

partially bypass energy flow around the PCM and directly to the 

engine, (4) introduction of concepts which can practically accommodate larger number of elements or introduce other effective methods 

of reducing the thickness of PCM through which heat transport occurs 

or (5) introduction of enhanced thermal conductivity materials. 

Within this range of design features it is possible to reduce 

peak PCM temperature from about 1210 to around 1090 K without the 

benefit of thermal enhancement. The base line receiver is found 

to provide relatively uniform engine hot side temperatures with 

variations of about 60 K. While certain of the concepts may in 

fact increase this value, reductions to about 45 K are possible 

while retaining the LiF 22 CaF2 PCM. 

The introduction of BN enhancement is found to marginally 

beneficial in terms of improvedthermal performance. The development 

of a graphite enhancement material would provide somewhat higher 

levels of performance and would receive a much higher priority. 

The development of Germanium or other metallic PCM offers a quite 

high level of thermal performance is, if containment or other 

potential development problems can be overcome, offer substantial 

improvements over all other enhancement schemes. 

a 

a 

0 
6 References 

a 

1. 

2. 

3. 

4. 

5. 

6. 

7. 

8. 

9. 

Wallin, W. E. & Dustin, M. O., "Advanced Space Solar Dynamic 

Power Systems Beyond IOC", 22nd IECEC, Philadelphia, Penn, Aug. 

1987. 

Vanco, M. R., "Analytical Comparison of kelative Heat Transfer 

Coefficients and Pressure Drops of Inert Gases and Their Binary 

Mixtures", Internal Report, NASA Lewis Research Center. 

Strumpf, H. J. & Coombs, M.G., "Advanced Heat Receiver Conceptual 

Design Studyt1, NASA CR-180901, December 1987. 

Juhasz, J. J., Coles-Hamilton, C. E. and Lacy, D. E., "Impact 

of Thermal Energy Storaged Properties on Solar Dynamic Space 

Power Conversion System Mass", 22nd IECEC, Philadelphia, August 

Kesseli, J. & Saunders, R., "Advanced Heat Receiver Conceptual 

Design Study", Presented at NASA Lewis Research Center, Oct. 

1986. 

Lacy, D. E. & Coles-Hamilton, C., "Selection of High Temperature 

Thermal Energy Storage Materials for Advanced Solar Dynamic 

Space Power SystemsI1, 22nd IECEC, Philadelphia, Penn, Aug. 1987. 

Diver, R. B. & Andraka, C. E., I1CAV2-A PC Based Computer Program 

for Predicting Incident Solar Flux Distributions Inside Dish 

Cavity Receivers" Sandia Report SAND86-0927, Feb. 1987. 

Patankar, S. V., Humeri cal Heat Tr ansfer and F1 uid Flow, 

Hemisphere Publishing Corp., New York, 1.980. 

Meredith, R. E. &I Tobias, C. W., llConductivities in Emulsions, 

Electro chemical Society Journ al, Vol. 108, 286, 1961. 

10-14, 1987. 

a 
APPENDIX A 

PROGRAM FOR THE THERMAL ANALYSIS 

OF THE 

BASE LINE RECEIVER 

AND THE 

HEAT PIPE RECEIVER 
0 

c*************************************************************** 

c*************************************************************** 

c*************************************************************** 

C * 

C PROGRAM SOLDTESS * 

0 

0 

a 

0 

0 

c * 

C SOLAR DYNAMIC TRANSIENT ENERGY STORAGE SIMULATION * 

C * 

C This program is the main routine for the TESS code. * 

C Its purpose is to read & print data and to call * 

C analysis. Checks for convergence are ..i;lde here. * 

C calculational subroutines so as to complete the thermal * 

C * c*************************************************************** 

C * 

C*************************************************************** 

C * 

C NOMENCLATURE 

C 

C System dimensions 

C 

C CS = Coefficient for Simpson 1/3 Integration 

C DR = Diff. element dimension in radial direction, m 

C DPHI = Diff. element dim. in azimuth.Ll direction, rad 

C DMAX = Maximum move in melt position in one iteration 

C DMELT = Fractional movement of melt position in element 

C DTIME = Time increment, sec 

C DZ = Diff. element dimension in axial directrion, m 

C GASMAS = Coolant flow rate, kg/sec 

C LENGTH = Axial length of stack, m 

C OVRRLX = Over relaxiation factor on melt movement 

C POS = Axial position of current node, m 

C R = Radius of element (I,J), m 

C RI = Outside radius of inner can ring, m 

C RICAN = Inside radius of inner can ring, m 

C RMELT = Fraction of PCM in element wlted 

C RN = Radius at outside edge of element (I,J), m 

C RO = Inside radius of outer can ring, m 

C ROCAN = Outside radius of outer can ring, m 

C RS = Radius at inside edge of node (1,J) , m 

C ROMELT = Previous fraction of PCM melted 

C T = Projected relative temperature of node 

C TcIJ = Absolute coolant temperature of gas at pos., K 

C TGAS = Relative coolant temperature at axial position 

C TGASIN = Absolute coolant temp. entering axial position,K 

C TGASP = Absolute coolant temperature at mid DZ, K 

C THGAS = Absolute coolant temperature at inlet, K 

C THETA = Absolute temperature of node (I,J), K 

C THETAG = Absolute temperature of gas, K 

C THMELT = Absolute melting temperature of storage media, K 

II) 
e 

A2 

0 

e 

0 

a 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

0 

D 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

TIME = Time into orbit run, min 

TIMER = Time of orbit, min 

TO = Previous relative temperature of node (1,J) 

TP = Projected relative temperature of node (1,J) 

WALLS = Ratio of wall conduction area to storage area 

Material 

CPGAS = 

CVLIQ = 

CVMET = 

CVPURE = 

CVSOL = 

DENLIQ = 

DENMET = 

DENPURE = 

DENSOL = 

properties 

Coolant specific heat, kJ/kg K 

Liquid specific heat, kJ/kg K 

Ring specific heat, kJ/kg K 

Specific heat of PCM (solid & liquid) kJ/kg K 

Solid specific heat, kJ/kg F 

Liquid density, kg/m**3 

Ring density, kg/m**3 

Density of PCM (solid t liquid) kg/m**3 

Solid density, kg/m**3 

F = Smooth tube adiabatic friction factor 

GASMIX = Fraction of gas 1 in coolant mix 

H = Convection coefficient of gas W/m**2 K 

HFL = Heat of fusion, kJ/kg 

K = Eff. thermal conductivity of node (I,J), W/m K 

KE = Thermal conductivity of node (I,J+l), W/m K 

KN = Thermal conductivity of node (I+l,J), W/m K 

KP = Thermal conductivity solid with void, W/m K 

KS = Thermal conductivity of node (I-l,J), W/m K 

KW = Thermal conductivity of node (1,J-1), W/m K 

K - 1 = Thermal conductivity constant (See GASPROP) 

K-2 = Thermal conductivity constant (See GASPROP) 

K - 3 = Thermal conductivity constant (See GASPROP) 

K-4 = Thermal conductivity constant (See GASPROP) 

KLIQ = Liquid thermal conductivity, W/m K 

KMET = Ring thermal conductivity, W/m K 

KPURE = Eff. thermal cond. of PCM (solid 61 liq.) W/m K 

KSOL = Solid thermal conductivity, W/m K 

MW = Molecular weight of gas component 

VISc-1 = Viscosity constant (see GASFPOP-)- 

VISC-2 = Viscosity constant (see GASPROP) 

VISC-3 = Viscosity constant (see GASPROP) 

VISC 4 = Viscosity constant (see GASPROP) 

VOID- = Void fraction in element 

Indicies 

I = Counter in radial direction 

ICOUNT = Loop counter 

IFLAG = Flag to indicate insufficient convergence 

IMAX = Maximum number of loops 

J =.Counter in azimuthal direction 

L = Counter in axial direction 
A3 

m 

a 

a 

a 

0 

0 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

LIO = Counter on output file 

NORBIT = Number of orbit cycles 

NI = Number of radial nodes 

NJ = Number of azimuthal nodes 

NL = Number of axial nodes 

NTIME = Number of time increments per orbit 

Matrix coefficients 

A = Conduction matrix coefficient 

B = Conduction matrix coefficient 

C = Conduction matrix coefficient 

D = Conduction matrix coefficient 

on 1-1 node 

on I node 

on 1+1 node 

on constant vector 

Thermal resistances & capacitances 

AE = Thermal conductance of J-1 node, W/K 

AN = Thermal conductance of 1+1 node,.W/K 

AS = Thermal conductance of 1-1 node, W/K 

AW = Thermal conductance of J+l node, W/K 

APO = Thermal capacitance of node (1,J) , kJ/m K 

Heat fluxes 

QAZ = Outside surface heat flux, W/m**2 

QDIR = Direct solar heat flux, W/m**2 

QGASP = Total heat added to gas at position DZ/2, kJ 

QGASDZ = Total heat added to gas at postion DZ, kJ 

QIJ = Radiation to element (I,J), W/m**2 

QIDIR = Indirect solar heat flux, W/m**2 

QINTR = Heat input on internal surface, W/m**2 

QZ = Integrated avg. axial multiplier on heat flux 

QZZ = Local axial multiplier on heat flux 

QZ- = Constant on fit of QZ * 

C*************************************************************** 

C 

C 

PROGRAM SOLDTESS 

REAL KLIQ,KMET,KSOL,LENGTH,MIN(ll) ,NTU 

DIMENSION FLOGAS(lO),QDIR(lO),QINDIR(lO),QINT(lO) 

DIMENSION TGAS(lO),THETA(lO),THGAS(lO,lO,lO),TP(l5,lO,lO) 

COMMON/BLOKA/CVLIQ,CVMET,CVSOL,DENLIQ,DENLIQ,DENMET,DENSOL,HFL, 

C 

, KLIQ, KMET, KSOL, LENGTH, PI, RI, RICAN, RO, ROCAN, 

,RMELT(15,10,10),ROMELT(l5,lO,lO),VOID,WALLS 

COMMON/BLOKB/F, GASMIX 

COMMON/BLOKC/DR, DTIME, DPHI, DZ ,HEAT,QAZ,QZ ,TO (15,10,10) 

COMMON/BLOKD/CPGAS,DELTDZ,DELTP,GASMAS,H,NTU,TGASIN,TGASP, 

' ,THMELT 
A4 

0 

a 

0 

e 

0 

0 

0 

0 

COMMON/BLOKE/A(15) ,B(15) ,C(15) ,D(15) ,T(15,10,10) 

COMMON/BLOKF/I,J,L,NI,NJ,NL 

COMMON/BLOKG/POS,QZl,QZ2,QZ3,QZ4 

C 

................................................................ 

C * 

C OPEN INPUT AND OUTPUT FILES * 

C * 

C*************************************************************** 

C 

OPEN (8, FILE=' INPUT. TXT' ) 

OPEN(7,FILE='OUTPUT.TXT') 

II=1 

PI=3.14159 

READ(8,2100,ERR=93)IMAX,NORBIT,NI,NJ,NL,NTIME 

READ(8,2200,ERR=94)RICAN,RI,RO,ROCAN,T..~NGTH,DUMl 

READ(8,2200,ERR=95) CVLIQ, DENLIQ,KLIQ,WALLS,THMELT,DUM2 

READ(8,2200,ERR=96)CVSOL,DENSOL,KSOL,ERRX,HFL,DUM3 

READ(8,2200,ERR=97)CVMET,DENMET,DENMET,~ET,GASMIX,TIMER,BORL 

READ(8,2200,ERR=98)QZl,QZ2,QZ3,QZ4,SORL 

1 READ (8,2200, ERR=99, END=5) MIN (11) , TGAS (11) , FLOGAS (11) , ,QDIR(II) ,QINDIR(II) ,QINT(II) 

II=II+l 

MIN (11) =l. E9 

GO TO 1 

C 

C***************************************~*********************** 

c * 

C INITIALIZE LOOP COUNTERS AND ERROR FLAGS * 

C * 

C*************************************************************** 

C 

5 DELTP=O.O 

ITEMP=O 

DR= (RO-RI ) /FLOAT (NI -2 ) 

DPHI=PI/FLOAT (NJ) 

DZ=LENGTH/FLOAT (NL) 

DTIME=60.O*TIMER/FLOAT(NTIME) 

DO 10 II=l,NI 

DO 10 JJ=l,NJ 

DO 10 LL=l,NL 

C 

T(I1, JJ, LL)=-50.0 

TP( 11, JJ, LL) =T (11, JJ, LL) 

C 

IF(I1 .LT. N1)THEN 

ELSE 

ENDIF 

RMELT(I1, JJ, LL) =O. 0 

RMELT ( I I, J J , LL) = 1.0 

8 
A5 

0 

e 

0 

a 

0 

a 

a 

a 

0 

I) 

C 

c*************************************************************** 

C * 

C * 

C*************************************************************** 

C 

10 CONTINUE 

C BEGIN CALCULATIONAL LOOPS ON LENGTH, TIME, ANGLE * 

DO 100 NO=l,NORBIT 

NPUT=O 

TIME=O. 0 

C 

DO 100 N=l,NTIME 

IF (MIN (NPUT+l) . LE. TIME+O .1E-6) NPUT=NPUT+l 

GASMAS=FLOGAS(NPUT) 

THGAS (N, NO, 1) =TGAS (NPUT) 

TIME=FLOAT (N) *DTIME/60. 0 

C 

DO 100 L=l,NL 

POS=DZ*FLOAT (L) 

TGAS IN=THGAS (N , NO, L) 

DO 20 II=l,NI 

DO 20 JJ=l,NJ 

C 

ROMELT (11, JJ, L) =RMELT (11, JJ, L) 

TO (11, JJ, L) =T(II, JJ, L) 

20 CONTINUE 

C 

ICOUNT=O 

ICOUNT=ICOUNT+l 

TGASP=THGAS (N, NO, L) +DELTP 

DO 40 J=l,NJ 

30 IFLAG=O 

C 

CALL GASPROP(RICAN,T(l,J,L) ,DZ) 

CALL SOLAR (LENGTH, QDIR (NPUT) , QINDIR (NPUT) , ROCAN) 

CALL CONDUCT (BORL, SORL, QINT (NPUT) ) 

CALL SOLVER 

CALL HEATGAS (T ( 1, J , L) ) 

40 CONTINUE 

C 

DO 50 II=l,NI 

DO 50 JJ=l,NJ 

ERROR=ABS (T( 11, JJ, L) -TP( 11, JJ, L) ) /THMELT 

IF(ERR0R .GT. ERRMAX)IFLAG=2 

IF(ERR0R .GT. ERRMAX)ITEMP=ITEMP+l 

TP ( 11, JJ , L) =T ( 11, JJ , L) 

50 CONTINUE 

C 

IF(IC0UNT .GT. 1MAX)THEN 

a 
A6 

C 

a 

C 

0 

a 

0 

8 

C 

C 

C 

C 

WRITE (7,1000) ITEMP 

IFLAG=O 

ENDIF 

IF(1FLAG .GE. l)GO TO 30 

THGAS (N,NO, L+1) =THGAS (N,NO, L) +DELTDZ 

IF(N0 .LT. NORB1T)GO TO 100 

WRITE(7,1050) 

WRITE(7,1250) 

WRITE(7,llOO) POS,TIME,NO 

WRITE (7,1150) QZ, HEAT 

DO 70 I=l,NI 

DO 60 J=l,NJ 

THETA (J) =T (I, J, L) +THMELT 

WRITE (7,1200) (THETA(J) , J=l,NJ) 

60 CONTINUE 

70 CONTINUE 

WRITE (7,13 00) TGASP, THGAS (N, NO, L+1) , H 

WRITE (7,1050) 

WRITE(7,1250) 

DO 80 I=l,NI-l 

WRITE(7,1350) (RMELT(1, J, L) , J=l,NJ) 

80 CONTINUE 

WRITE (7,1350) (RMELT ( 1, J, L) , J=l, NJ) 

100 CONTINUE 

c****************************************** .................... C * 

C INPUT AND OUTPUT FORMAT STATEMENTS * 

C * 

C*************************************************************** 

C * 

1000 FORMAT(lX,'FAILED TO CONVERGE AFTER ',16,' ITERATIONS.') 

1050 FORMAT(lX,' ') 

1100 FORMAT(lX,'Axial Position = ',1F5.3,' m. Time = I, 

1150 FORMAT(lX,'Heat Input Frac. = ',1F6.4,' Heat = ',1F9.4) 

1200 FORMAT(lOF7.1) 

& 1F6.1,' min. Cycle = ',Illt.') 

1250 FORMAT(lX,?------------------------------------------------ ___----------------------------------___________ 

I 1 &==? ?--------------------I _____--------------- 

CF6.1,'K H = ',F7.4,'KW/m2-Kt) 

1300 FORMAT(lX,'Avg. Gas Temp. = ',F6.1,'K Out. Gas temp.', 

1 3 5 0 FORMAT ( 10 F7 .4 ) 

1400 FORMAT(lX,'INPUT FORMAT ERROR IN 1ST LINE') 

1401 FORMAT(lX,'INPUT FORMAT ERROR IN 2ND LINE') 

1402 FORMAT(lX,'INPUT FORMAT ERROR IN 3RD LINE') 
A7 

e 

a 

a 

e 

a 

1403 FORMAT(lX,’INPUT FORMAT ERROR IN 4TH LINE’) 

1404 FORMAT(lX,’INPUT FORMAT ERROR IN 5TH LINE‘) 

1405 FORMAT(lX,’INPUT FORMAT ERROR IN 6TH LINE‘) 

1406 FORMAT(lX,’INPUT FORMAT ERROR IN TRANSIENT SECTION’) 

2100 FORMAT(16,5110,2F10.0) 

2200 FORMAT (2F10.4,4E10.4) 

STOP 

STOP 

STOP 

STOP 

STOP 

STOP 

STOP 

STOP 

END 

93 WRITE(*,1400) 

94 WRITE (*, 1401) 

95 WRITE (*, 1402) 

96 WRITE(*,1403) 

97 WRITE(*,1404) 

98 WRITE (*, 1405) 

99 WRITE(*,1406) 
A8 

0 

0 

0 

0 

C*******************.******************************************* 

C************************************************************** 

C * 

C SUBROUTINE SOLAR * 

C * 

C This program calculates the solar incidence on * 

C the outside perimeter of the containment ring * 

C q1! = q" (direct) *cos (angle) + q" (indirect) * 

C * 

C * 

C*****************X******************************************** 

C************************************************************** 

C 

C 

SUBROUTINE SOLAR(LENGTH,QDIR,QINDIR,ROCAN) 

DIMENSION CS (21) 

REAL LENGTH 

COMMON/BLOKC/DR, DTIME, DPHI, DZ, HEAT, QAZ, QZ ,TO (15,10,10) 

COMMON/BLOKF/I,J,L,NI,NJ,NL 

COMMON/BLOKG/POS,QZl,Q22,QZ3,QZ4 

DATA CS(1) ,CS(3) ,CS(5) ,CS(7) ,CS(9) ,CS(11) ,CS(13)/ 

/l. 0,6*2.0/ 

DATA CS(15) ,CS(17) ,CS(19) ,CS(2) ,CS(4) ,CS(6) ,CS(8)/ 

/3*2.0,4*4.0/ 

DATA CS(10) ,CS(12) ,CS(14) ,CS(16) ,CS(18) ,CS(20) ,CS(21)/ 

/6*4., 1 ./ 

DATA PI/3.14159/ 

C 

C*************************************************************** 

C * 

C QZ USES SIMPSON'S 1/3 INTEGRATION TO FIND THE * 

C INTETRATED AVERAGE VALUE OF THE AXIAL POWER CURVE * 

C OVER THE INTERVAL DZ * 

C * 

C*************************************************************** 

C 

XL= (POS-DZ)/LENGTH 

QZ=O. 0 

DO 15 II=1,21 

C 

IF(I1 .GT. l)XL=XL+O.OS*DZ/LENGTH 

QZZ=QZl*XL*EXP(l.O-(XL/QZ2)**QZ3)+QZ4*(XL)**O.5 

QZ=QZ+CS(II) *(QZZ)/60.0 

15 CONTINUE 

C 

C******************************************~******************* 

C * 

C FIND THE PROJECTION OF THE DIRECT SOLAR INCIDENCE * 

C * 

c************************************************************** 
e 

A9 

C 

* 

0 

a 

e 

ANGLEl=PI*FLOAT(J)/FLOAT(NJ) 

ANGLE2=PI*FLOAT(J-l)/FLOAT(NJ) 

IF (J+J-2 .GT. NJ) THEN 

ELSEIF (J+3 .GT. NJ)THEN 

ELSE 

ENDIF 

C 

. SINANGLE=O. 0 

SINANGLE=l.O-SIN(ANGLE2) 

SINANGLE=SIN(ANGLEl)-SIN(ANGLE2) 

C 

C************************************************************** 

C * 

C ADD DIRECT AND INDIRECT INCIDENT SOLAR RADIATION * 

C * ............................................................... 

C 

QAZ=(QDIR*SINANGLE + QINDIR*DPHI)*ROCAN*DZ*QZ 

IF(J .EQ. l)HEAT=O.O 

HEAT=HEAT+2.0*QAZ 

RETURN 

END 

C 

e 
e 

A10 

e 

* 

e 

e 

f? 

* 

C******************************************~.~****************** 

C************************************************************** 

C * 

C SUBROUTINE GASPROP * 

C * 

C properties and convective coefficients for the gas coolant * 

C This program calculates the effective thermal * 

C The gas coolant is assumed to be a mixture of He and Xe. * 

C Properties for these two constituents are included as * 

C curve fits with the corresponding coefficients included * 

C as data within the subroutine. * 

C * 

C******************************************~********~********** 

C************************************************************** 

C 

C************************************************************** 

C * 

C REFERENCES * 

C M. R. Vanco, ttAnalytical Comparison of Relative Heat * 

C Gases and Their Binary Mixtures", Lewis Research Center* 

C * 

C Transfer Coefficients and Pressure Drops of Inert * 

C A. J. Chapman ttFundamentals of Heat Transfer" * 

C Macmillan, 1987, p 329. * 

C * 

C************************************************************** 

C 

C************************************************************** 

C * 

C EQUATIONS * 

C * 

C K1 = K11 + K12*(TGASP - K13)**K14 * 

c K2 = K21 + K22*(TGASP - K23)**K24 * 

C VISCl = VISCll + VISC12*(TGASP - VISC13)**VISC14 * 

C VISC2 = VISC21 + VISC22*(TGASP - VISC23)**VISC24 * 

C * 

C * 

C************************************************************** 

C 

C 

SUBROUTINE GASPROP(RICAN,TlJ,DZ) 

REAL K,K1,K2,Kll,K12,K13,Kl4,K2l,K22,K23,K24,KGAS, 

+MW,MWl,MW2,NTU,NU 

C 

COMMON/BLOKB/F, GASMIX 

COMMON/BLOKD/CPGAS,DELTDZ,DELTP,GASMAS,H,NTU,TGASIN,TGASP, 

, THMELT 

C 

DATA GAMMA,PI,RUNIV,MWl,MW2/1.66,3.14159,8.31434,4.,131./ 

DATA Kll,K12,K13,K14/0.158OE-3,0.9086E-6,333.333,0.794/ 
All 

0 

I 

0 

e 

e 

0 

0 

0 

e 

DATA K21,K22,K23,K24/0.6404E-5I0.0454E-6,3333.333,0.800/ 

DATA VISC11,VISC12,VISC13,vISc14/ 

/.2158E-4,.1022E-6,333.3,.825/ 

/.253E-4,.2452E-6,333.333,.77/ 

DATA VISC21,VISC22,VISC23,vISc24/ 

C 

F1 (XI Y) 

F2 (X1,X2 ,X3 I X4 I X5 I X6) = 

= ( 1.0+2.41* (X-Y) * (X-0 .142*Y) / (X+Y) **2 ) 

= ( (l.O+ (Xl/X2) **O. 5* (X3/X4) **O. 25) **2)/ 

/ (2.0*2.0**0.5* (l.O+X5/X6) **O. 5) 

C 

C*********************************************************** 

C * 

C FIND THE THERMAL CONDUCTIVITY &I VISCOSIT" OF PURE GASES * 

C * 

C*********************************************************** 

C 

Kl=Kll+K12* (TGASP-K13) **K14 

K2=K21+K22*(TGASP-K23)**K24 

VISC1=VISC11+VISC12*(TGASP-VISCl3)**VISCl4 

VISC2=VISC21+VISC22*(TGASP-VISC23)**VISC24 

VISCW1=VISC11+VISC12*(TlJ+THMELT-VISCl3)**VISCl4 

VISCW2=VISC21+VISC22*(TlJ+THMELT-VISC23)**VISC24 

C 

C*********************************************************** 

C * 

C FIND THE WEIGHTING PARAMETERS FROM REFERENCE * 

C * 

C*********************************************************** 

C 

PH112=F2(VISC1,VISC2,MW2,MW2fMWlfMWlfMW2) 

PHI21=PHI12*(VISC2*MWl)/(VISCl*MW2) 

PSI12=F2(K1,K2,MW1,MW2,MWl,MW2)*Fl(MwlIMW2) 

PS121=F2(K2,K1,MW2,MWlfMW21MWl)*Fl(MW2fMWl) 

C 

C*********************************************************** 

C * 

C SOLVE FOR THERMAL CONDUCTIVITY AND VISCOSlrY OF GAS MIX * 

C * 

C*********************************************************** 

C 

VISC=VISC1/(1.O+PHI12*(l.O-GASMIX)/GASMIX)+ 

+ VISC2/(l.O+PHI21*(GASMIX/(l.O-GASMIX))) 

VISCW=VISCWl/(l.O+PHIl2*(l.O-GASMIX)/GASMIX)+ 

+ VISCW2/(l.O+PHI21*(GASMIX/(l.O-GASMIX))) 

KGAS=Kl/(l.O+PSI12*(l.O-GASMIX)/GASMIX)+ 

+ K2/(l.O+PSI21*(GASMIX/(l.O-GASMIX))) 

C 

C*********************************************************** 

C * 
A12 

e 

e 

a 

0 

C SOLVE FOR SPECIFIC HEAT OF GAS MIX * 

C * 

C*********************************************************** 

C 

MW=GASMIX*MWl+ (1.0-GASMIX) *MW2 

CPGAS=GAMMA*RUNIV/ ( (GAMMA-1.0) *MW) 

C 

C 

C*********************************************************** 

C * 

C .SOLVE FOR CONVECTION COEFFICIENT * 

C USING THE SIEDER-TATE CORRELATION FOR TURBULENT FLOW * 

C AND THE NUSSELT CORRELATION FOR LAMINAR FLOW * 

C FIND THE NUMBER OF TRANSFER UNITS * 

C * 

C******************************************~**************** 

C 

PR=CPGAS*VISC/ (KGAS) 

RE=2.0*GASMAS/(PI*RICAN*VISC) 

IF(RE .LT. 10000.0) THEN 

C 

F=64.O/RE 

NU=3.66 

ELSE 

F= ( 1.82 *ALOG10 (RE) -1.64 ) ** ( -2 ) 

NU=0.027*RE**0.8*PR**O.333333*(VISC/VISC~~)**O.l4 

ENDIF 

H=KGAS*NU/ (2.O*RICAN) 

NTU=H*2.0*PI*RICAN*DZ/(GASMAS*CPGAS) 

RETURN 

END 

C 

C 
A13 

0 

a 

c*************************************************************** 

c*************************************************************** 

C * 

C SUBROUTINE CONDUCT * 

C * 

C PROGRAM DEVELOPS FINITE DIFFERENCE MATRIX * 

C FOR CONDUCTION ANALYSIS IN CYLINDRICAL COORDINATES * 

C * ................................................................ 

................................................................ 

C * 

C REFERENCE: Pantankar, S. V., ''Numerical Heat Transfer and * 

C * c*************************************************************** 

C 

C 

C 

* 

C Fluid Flow", McGraw Hill, 1980, pp. 64-74. * 

SUBROUTINE CONDUCT(BORL,SORL,QINT) 

REAL K, KE, KLIQ, KMET, KN, KP, KS, KSOL, KW, LENGTH,NTU 

COMMON/BLOKA/CVLIQ,CVMET,CVSOL,DENLIQ,DENMET,DENSOL,HFL, 

+KLIQ, KMET, KSOL, LENGTH, PI, RI , RICAN, RO, ROCAN, 

+RMELT(15,10,10),ROMELT(l5,lO,lO),VOID,WALLS 

COMMON/BLOKC/DR,DTIME,DPHI,DZ,HEAT,QAZ,QZ,TO(l5,10,10) 

COMMON/BLOKD/CPGAS,DELTDZ,DELTP,GASMAS,H,NTU,TGASIN,TGASP, 

COMMON/BLOKE/A ( 15) , B ( 15) , C ( 15) , D ( 15) , T ( 15,10,10) 

COMMON/BLOKF/I , J, L, NI , NJ, NL 

DO 50 I=l,NI 

IF (J .EQ. 1) THEN 

, THMELT 

C 

C 

AE=O. 0 

TE=O. 0 

TE=T ( I , J - 1 , L ) 

ELSE 

ENDIF 

IF (J .EQ. NJ) THEN 

C 

AW=O. 0 

TW=O. 0 

TW=T (I, J+1, L) 

ELSE 

ENDIF 

C 

C 

C************************************************************** 

C * 

C SET NODAL CONSTANTS FOR INSIDE RING * 
e 

A14 

a 

a 

e 

C * 

C************************************************************** 

C 

IF (I .EQ. 1) THEN 

R=RICAN 

a RN=RI 

VOID=( 1. 0-DENLIQ/DENSOL) * (l.O-RMELl( 2, J, L) ) 

CALL PROPERTY (CV,DEN,2,J,KN,L) 

AA=O. 5* (RN*RN-R*R) *DPHI*KMET/DZ 

AB=0.5*(RN*RN-R*R)*DPHI*KMET/DZ 

AN=KMET*KN*RN*DPHI*DZ/(O.5*KMET*DR+KN*(RI-RICAN)) 

AS=O. 0 

IF(J. GT. 1) 

IF(J .LT. NJ) 

APO~DENMET*CVMET*0.5*(RI*RI-RICAN*RICAN)*DPHI*DZ/DTIME 

TCIJzTGASIN-THMELT 

HIJ=GASMAs*CPGAS*(l.O-EXP(-NTU))/FLOAT(2*NJ) 

QIJ=QINT/(2.0*PI*R*LENGTH) 

& 

& AW=KMET*(RI-RICAN)*DZ/(O.5*(RI+RICAN)*DPHI) 

AE=KMET* (RI-RICAN) *DZ/ (0.5" (RI+RICAN) *DPHI) 

C 

C************************************************************** 

C * 

C SET NODAL CONSTANTS FOR OUTSIDE RING * 

C 

C************************************************************** 

C 

* a 

ELSEIF (I .EQ. NI) THEN 

R=ROCAN 

RS=RO 

VOID=(l.O-DENLIQ/DENSOL)*(l.O-RMELI'(NI-1,J,L)) 

CALL PROPERTY (CV,DEN,I-l,J,KS,L) 

AA=O.5*(R*R-RS*RS)*DPHI*KMET/DZ 

AB=0.5*(R*R-RS*RS)*DPHI*KMET/DZ 

AS=KMET*KS*RS*DPHI*DZ/(O.5*KMET*DR+KS*(ROCAN-R0)) 

IF(J .GT. 1) 

IF(J .LT.NJ) 

APO~DENMET*CVMET*0.5*(ROCAN*ROCAN-RO*RO)*DPHI*DZ/DTIME 

HIJ=O. 0 

QI J=QAZ 

a AN=O. 0 

& AE=KMET*(ROCAN-RO)*DZ/((.5*(ROCAN-RO)+RO)*DPHI) 

& AW=KMET* (ROCAN-RO) *DZ/ ( ( .5* (ROCAN-RO) +RO) *DPHI) 

TCIJ=O. 0 

C 

C************************************************************** 

C * 

C * 

C************************************************************** 

e C SET NODAL CONSTANTS FOR PHASE CHANGE MATERIAL * 
0 

C 

ELSE 

A15 

0 

C 

0 

e 

+ 

C * 

+ 

e C 

R=RI+(FLOAT(I)-1.5) *DR 

RN=R+0.5*DR 

RS=R-0.5*DR 

VOID= ( 1 . 0-DENLIQ/DENSOL) * (1.0-RMELT (I, J, L) ) 

CALL PROPERTY (CV, DEN, I, J, K, L) 

AA=O.S*(RN*RN-RS*RS)*DPHI*K/DZ 

AB=0.5*(RN*RN-RS*RS)*DPHI*K/DZ 

IF (I .EQ. 2) THEN 

AS=KMET*K*RI*DPHI*DZ/(O.5*KMET*DR+K*(RI-RICAN)) 

TS=T (2, J, L) +DR*AS* (T (1, J, L) -T (2, J, L) ) / 

0.125* (T (2, J, L) -2.O*T (3, J, L) +T (4, J, L) ) 

CALL PROPERTY (CV,DEN,I-l,J,KS,L) 

KP=2.O*K*KS/ (K+KS) 

AS=KP*RS*DPHI*DZ/DR 

TS=O. 5* (T (I, J, L) +T (1-1, J, L) ) + 

/ (2.O*K*RI*DPHI) + 

+ 

ELSE 

0.125* (T( I+1, J, L) -2.O*T( I, J, L) +T(I-1, J, L) ) 

ENDIF 

IF (I .EQ. NI-1) THEN 

AN=KMET*K*RO*DPHI*DZ/(O.5*KMET*DR+K*(ROCAN-RO)) 

TN=T (NI-1, J, L) +DR*AN* (T (NI, J, L) -T (NI-1, J, L) ) / 

(2. *K*RO*DPHI) +O. 125* (T(N1-1, J, L) - 

2.O*T (NI-2, J, L), +T (NI-3, J, L) ) 

(CV , DEN, I+1, J , KN, L) 

ELSE 

CALL PROPERTY 

KP=2.O*K*KN/ (K+KN) 

AN=KP*RN*DPHI*DZ/DR 

TN=O. 5* (T( I, J, L) +T (I+1, J, L) ) + 

0.125* (T (I+1, J, L) -2.O*T (I, J, L) +T (1-1, J, L) ) 

ENDIF 

c*********************************.***************************** 

C * 

C Calculate the melt front movement * 

C * c************************************************************** 

e C 

IF 

* 

(TS*TN .GE. 0.0) THEN 

OVRRLX = BORL*ABS (T (I, J, L) ) 

DMELT = OVRRLX*0.5*(TN+TS)*CV/HFL 

IF (TS .GE. 0.0) THEN 

DMELT=AMINl (DMELT, 1.0-WELT (I, J, L) ) 

ELSE 

DMELT=AMAXl(DMELT,O.O-WELT(I,J,L)) 

ENDIF 
A16 

/ 

m 

a 

0 

/ 

ELSE 

IF(T(I,J,L)*TN .GE. 0.0) THEN 

IF(TN .GE. 0.O)THEN 

ELSE 

DMELT=( 0.5*TN/ (TN-T( I, J, L) ) -MELT( I, J, L) ) *SORL 

DMELT=((O.5*TN-T(I,J,L))/ 

(TN-T( I, J, L) ) -RMELT(I, J, L) ) *SORL 

ENDIF 

IF(TN .LE. 0.O)THEN 

ELSE 

DMELT= ( .5*TS/ (TS-T ( I, J, L) ) -RMELT ( I, J, L) ) *SORL 

ELSE 

DMELT=( (0.5*TS-T (I, J, L) ) / 

(TS-T (I, J, L) +WELT( I, J, L) *SORL 

ENDIF 

ENDIF 

ENDIF 

C 

C************************************************************** 

C * 

C Check melt front movement to ensure that it meets * 

C thermodynamic limitations. * 

C * c************************************************************** 

C 

40 

C 

& 

& 

& 

& 

& 

& 

& 

& 

IF (DMELT . GE. 0.0) THEN 

IF (RMELT ( 11, J , L) 

DO 40 II=I+l,NI-l 

CONTINUE . LT . 0.99 ) DMELT=O . 0 

ENDIF 

IF(1 .EQ. 2)THEN 

IF(DMELT .LT. 0.O)THEN 

IF((RMELT(I+l,J,L) .LE. .01) .OR. 

(T(1-l,J,L) .LE. 0.0)) 

WELT( I, J, L) =MELT (I, J, L) +DMELT 

ELSE 

IF( (RMELT(I+l,J,L) .GE. 0.99) .OR. 

(T(1-l,J,L) .GE. 0.0)) 

RMELT (I, J, L) =WELT( I, J, L) +DMELT 

ENDIF 

IF(DMELT .GE. 0.O)THEN 

ELSEIF(1 .EQ. NI-1)THEN 

IF((T(NI,J,L) .GE. 0.0) .OR. 

(QAZ .GE. 0.0)) 

WELT( I, J, L) =WELT( I, J, L) +DMELT 

ELSE 

IF((T(NI,J,L) .LE. 0.0) .OR. ~ (RMELT(NI-2,J,L) .LE. .01)) 

WELT( I, J, L) =RMELT( I, J, L) +DMELT 

, 
A17 

I '* 

* 

0 

C 

C 

C 

a 

C 

0 

C 

0 

C 

ENDIF 

IF(DMELT .LT. 0.O)THEN 

& (RMELT(I+l,J,L) .LE. .01)) 

& 

ELSE 

IF( (RMELT(1-l,J,L) .LE. .01) .OR. 

RMELT (I, J, L) =MELT (I, J, L) +DMELT 

ELSE 

IF((RMELT(1-l,J,L) .GE. .99) .OR. 

& (RMELT(I+l,J,L) .GE. .99)) 

& RMELT (I, J, L) =WELT (I, J, L) +DMELT 

ENDIF 

ENDIF 

RMELT (I, J, L) =AMAX1 (0.0 ,WELT( I, J, L) ) 

RMELT( I, J, L) =AMIN1 (1.0, RMELT (I, J, L) ) 

IF (J .GT. 1) THEN 

CALL PROPERTY (CV, DEN, I, J-1, KE, L) 

KP=2.O*K*KE/ (K+KE) 

AE=KP*DR*DZ/ (R*DPHI) 

ENDIF 

IF(J .LT. NJ) THEN 

CALL PROPERTY (CV, DEN, I, J+1, KW, L) 

KP=2.0*K*KW/(K+KW) 

AW=KP*DR*DZ/(R*DPHI) 

ENDIF 

APO=DEN*CV*R*DPHI*DR*DZ/DTIME 

TCIJ=O. 0 

HIJ=O. 0 

PCMMASS=R*DPHI*DR*DENLIQ*(l.O-WALLS)*DZ 

QIJ=(ROMELT(I,J,L) -RMELT(I,J,L) ) *PCMMASS*HFL/DTIME 

ENDIF 

IF (L .EQ. 1) THEN 

AA=o. 0 

TA=O . 0 

TA=T ( I, J , L-1) 

ELSE 

ENDIF 

IF (L .EQ. NL) THEN 

AB=O. 0 

TB=O. 0 

TB=T ( I, J , L+1) 

ELSE 

ENDIF 

A (I) =-AS 
A18 

B(I)=AA+AB+AN+AS+AE+AW+APO+HIJ 

C (I) =-AN 

D(I)=APO*TO(I,J,L)+HIJ*TCIJ+QIJ+AE*TE+AW*TW+AA*TA+AB*TB 

C 

C 

50 CONTINUE 

RETURN 

END 
A19 

0 

a 

a 

c************************************************************** 

C*****************************************************~********* 

C' * 

.c SUBROUTINE PROPERTY * 

C * 

C This program calculated the effective thermal * 

C properties of the solid & liquid phases oir' the storage media* 

C * 

C************************************************************** 

C************************************************************** 

C************************************************************** 

C * 

0 

C REFERENCES Y. S. Touloukian et.al. l1Thermophysical * 

C Properties of Matter", Vol. 2, 1970. * 

C * 

C R. P. Wichner, ORNL IIBoundary Conditions and Thermophysical * 

C Properties of LiF-CaF2I8, Presented to NASA Lewis Research * 

C Center, June 10, 1987. * 

C * 

C R. I. Vachon, et.al, "Thermal Conductivity of Heterogeneous * 

C Mixtures & Lunar Soils11, NAS8-26579, October 1973. * 

C * 

C************************************************************** 

C 

C 

SUBROUTINE PROPERTY (CV, DEN, I, J, K, L) 

COMMON/BLOKA/CVLIQ, CVMET, CVSOL, DENLIQ, DENMET, DENSOL, HFL, 

+KLIQ,KMET,KSOL,LENGTH,PI,RI,RICAN,RO,ROCAN, 

,RMELT(15,10,10),ROMELT(l5,lO,lO),VOID,WALLS 

REAL K, KLIQ, KMET, KP, KPSOL, KPURE,KSOL, LENGTH 

C 

C 

C 

C************************************************************** 

C * 

C FIND THE THERMAL PROPERTIES OF THE PURE MATERIALS * 

C * 

C CALCULATE EFFECT OF VOID ON PROPERTIES OF SOLID * 

C USE CORRELATION BY MEREDITH & TOBIAS TO CORRECT CONDUCTIVITY* 

C * 

C************************************************************** 

C 

KPSOL~KSOL*(2.0*(1.0-VOID)+O.6135*VOI~**2.333-1.6*VOID** 

& 3.333)/(2.O+VOID+0.6135*V01D**2.333-0.6975*V01D**3.333) 

DENPSOL= DENLIQ 

KPURE=l. O/ (RMELT (I, J, L) /KLIQ+ (1.0-WELT (I, J, L) ) /KPSOL) 

CVPURE=RMELT(I,J,L)*CVLIQ+(l.O-RMELT(I,J,L))*CVSOL 

DENPURE=RMELT (I, J, L) *DENLIQ+ ( 1.0-WELT (I, J, L) ) *DENPSOL 

C 

a 
A2 0 

C 

C************************************************************** 

C * 

C PUT PURE STORAGE MATERIAL IN PARALLEL WITH CAN MATERIAL * 

C * 

C************************************************************** 

C 

e 

K=KPURE*(l.O-WALLS)+KMET*WALLS 

DEN=DENPURE*(l.O-WALLS)+WALLS*DENMET 

CV=(CVPURE*DENPURE*(l.O-WALLS)+WALLS*CVMET*DENMET)/DEN 

RETURN 

END 

0 C 

a 
A2 1 

0 

C************************************************************** 

C************************************************************** 

C * 

C SUBROUTINE SOLVER * 

C * 

C This program solves a tridiagonal matrix using * 

C an elimination method * 

C * 

C [A]T(I-1) + [B]T(I) + [C]T(I+l) = [D] * 

C * ............................................................... 

C************************************************************** 

C 

SUBROUTINE SOLVER 

COMMON/BLOKE/A ( 15) , B ( 15) , C ( 15) , D ( 15) , T ( 15 , 10 , 10) 

COMMON/BLOKF/I, J, L, NI, NJ, NL 

C 

C************************************************************** 

C * 

C PROCEED WITH GAUSSIAN ELIMINATION TO DEVELOP A * 

C UPPER DIAGONAL MATRIX * 

C * 

C************************************************************** 

C 

DO 1, I=l,NI-l 

R=A ( I+ 1 ) / B ( I ) 

B(I+l)=B(I+l) -R*C(I) 

D(I+l)=D(I+l) -R*D(I) 

1 CONTINUE 

C 

C************************************************************** 

C * 

C BEGIN BACK SUBSTITUTION * 

C * 

c************************************************************** 

C 

C 

T (NI , J, L) =D (NI) /B (NI) 

DO 2, I=NI-1,1,-1 

T( I , J, L)=(D( I) -C (I) *T (I+1, J, L) )/B( I) 

2 CONTINUE 

RETURN 

END 
APPENDIX B 

PROGRAM FOR THE THERMAL ANALYSIS 

OF THE 

CAVITY RECEIVER 
B1 

e 

0 

a 

e 

a 

a 

9 

e 

C******************************************************************* 

C******************************************************************* 

C******************************************************************* 

C * 

C PROGRAM SOLDTESS * 

C * 

C SOLAR DYNAMIC TRANSIENT ENERGY STORAGE SIMULATION * 

C * 

C This program is the main routine for the TESS code. * 

C Its purpose is to read C print data and to call * 

C calculational subroutines so as to complete the thermal * 

C analysis. Checks for convergence are made here. * 

C * 

C******************************************************************* 

C 

C NOMENCLATURE 

C 

C CS = Coefficient for Simpson 1/3 Integration 

C DR = Diff. element dimension in radial direction, m 

C DPHI = Diff. element dim. in azimuthal direction, rad 

C DMAX = Maximum move in melt position in one iteration 

C DMELT = Fractional movement of melt position in element 

C DTIME = Time increment, sec 

C DZ = Diff. element dimension in axial directrion, m 

C GASMAS = Coolant flow rate, kg/sec 

C LENGTH = Axial length of stack, m 

C OVRRLX = Over relaxiation factor on melt movement 

C POS = Axial position of current node, m 

C R = Radius of element (I,J), m 

C RI = Outside radius of inner can ring, m 

C RICAN = Inside radius of inner can ring, m 

C RMELT = Fraction of PCM in element melted 

C RN = Radius at outside edge of element (I,J), m 

C RO = Inside radius of outer can ring, m 

C ROCAN = Outside radius of outer can ring, m 

C RS = Radius at inside edge of node (1,J) , m 

C ROMELT = Previous fraction of PCM melted 

C T = Projected relative temperature of node 

C TCIJ = Absolute coolant temperature of gas at pos., K 

C TGAS = Relative coolant temperature at axial position 

C TGASIN = Absolute coolant temp. enterilly axial position,K 

C TGASP = Absolute coolant temperature at mid DZ, K 

C THGAS = Absolute coolant temperature at inlet, K 

C THETA = Absolute temperature of node (I,J), K 

C THETAG = Absolute temperature of gas, K 

C THMELT = Absolute melting temperature of storage media, K 

C TIME = Time into orbit run, min 

C TIMER = Time of orbit, min 

.C System dimensions 
B2 

0 

e 

a 

a 

0 

a 

0 

e 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

TO = Previous relative temperature of node (1,J) 

TP = Projected relative temperature of node (1,J) 

WALLS = Ratio of wall conduction area to storage area 

Material properties 

CPGAS = Coolant specific heat, kJ/kg K 

CVLIQ = Liquid specific heat, kJ/kg K 

CVMET = Ring specific heat, kJ/kg K 

CVPURE = Specific heat of PCM (solid & Liquid) kJ/kg K 

CVSOL = Solid specific heat, kJ/kg K 

DENLIQ = Liquid density, kg/m**3 

DENMET = Ring density, kg/m**3 

DENPURE = Density of PCM (solid & liquid) kg/m**3 

DENSOL = Solid density, kg/m**3 

F = Smooth tube adiabatic friction factor 

GASMIX = Fraction of gas 1 in coolant mix 

H = Convection coefficient of gas W/m**2 K 

HFL = Heat of fusion, kJ/kg 

K = Eff. thermal conductivity of node (I,J), W/m K 

KE = Thermal conductivity of node (I,J+l) , W/m K 

KN = Thermal conductivity of node (I+l,J), W/m K 

KP = Thermal conductivity solid with void, W/m K 

KS = Thermal conductivity of node (I-l,J), W/m K 

KW = Thermal conductivity of node (1,J-1), W/m K 

K 1 = Thermal conductivity constant (See GASPROP) 

K-2 = Thermal conductivity constant (See GASPROP) 

K-3 = Thermal conductivity constant (See GASPROP) 

K-4 = Thermal conductivity constant (See GASPROP) 

KXIQ 

KMET = Ring thermal conductivity, W/m K 

KPURE = Eff. thermal cond. of PCM (solid & liq.) W/m K 

KSOL = Solid thermal conductivity, W/m K 

MW- = Molecular weight of gas componenl; - VISC-1 = Viscosity constant (see GASPROP) 

VISC-2 = Viscosity constant (see GASPROP) 

VISC 3 = Viscosity constant (see GASPROP) 

VISC-4 = Viscosity constant (see GASPROP) 

VOID- = Void fraction in element 

= Liquid thermal conductivity, W/m K 

Indicies 

I = Counter in radial direction 

ICOUNT = Loop counter 

IFLAG = Flag to indicate insufficient convergence 

IMAX = Maximum number of loops 

J = Counter in azimuthal direction 

L = Counter in axial direction 

LIO = Counter on output file 

NORBIT = Number of orbit cycles 

NI = Number of radial nodes 

NJ = Number of azimuthal nodes 

a 
B3 

C NL = Number of axial nodes 

C NTIME = Number of time increments per orbit 

C 

C Matrix coefficients 

C A = Conduction matrix coefficient on 1-1 node 

C B = Conduction matrix coefficient on I node 

C c = Conduction matrix coefficient on 1+1 node 

C D = Conduction matrix coefficient on constant vector 

C 

C Thermal resistances 61 capacitances 

C AE = Thermal conductance of J-1 node, W/K 

C AN = Thermal conductance of 1+1 node, W/K 

C AS = Thermal conductance of 1-1 node, W/K 

C AW = Thermal conductance of J+l node, W/K 

C APO = Thermal capacitance of node (I,J), kJ/m K 

C 

C Heat fluxes 

C QAZ = Outside surface heat flux, W/m**2 

C QDIR = Direct solar heat flux, W/m**2 

C QGASP = Total heat added to gas at position DZ/2, kJ 

C QGASDZ = Total heat added to gas at postion DZ, kJ 

C QIJ = Radiation to element (1,J) , W/m*,:2 

C QIDIR = Indirect solar heat flux, W/m**2 

C QINTR = Heat input on internal surface, W/m**2 

C QZ = Integrated avg. axial multiplier on heat flux 

C QZZ = Local axial multiplier on heat flux 

C QZ- = Constant on fit of QZ 

C 

C****************~************************************************** 

C * 

PROGRAM SOLDTESS 

REAL KLIQ,KMET,KSOL,LENGTH,MIN(ll) ,NTU 

DIMENSION THETA(lO),TP(15,10) 

COMMON/BLOKA/CVLIQ, CVMET, CVSOL, DENLIQ, DENMET, DENSOL, HFL, 

COMMON/BLOKC/DTIME , DX, DY , TO (2 5,20) 

COMMON/BLOKE/A (25) , B (2 5) , C (2 5) , D (25) , Q ( 10) , QNA (2 0) , T (25,2 0) 

COMMON/BLOKF/I , J, NI , NJ 

DATA U/400.0/ 

C 

C 

, KLIQ, KMET, KSOL, LENGTH, WELT (25,20) , ROMELT (25,20) , VOID 

C c****************************************************************** 

C * 

C OPEN INPUT AND OUTPUT FILES * 

C * c****************************************************************** 

C 

OPEN (8, FILE='INPUT. TXT' ) 

OPEN (7, FILE='OUTPUT. TXT' ) 
B4 

0 

0 

II=1 

READ(8,2100,ERR=93)IMAX,NORBIT,NI,NJ,NTIME 

READ(8,2200,ERR=94)H,Hl,W,BORL,SORL 

READ(8,2200,ERR=95)CVLIQ,DENLIQ,KLIQ,LENGTH,THMELT 

READ(8,2200,ERR=96)CVSOL,DENSOL,KSOL,ERRX,HFL 

READ(8,2200,ERR=97)CVMET,DENMET,KMET,GASMIX,TIMER 

II=II+l 

MIN (11) =l. E9 

GO TO 1 

1 READ(8,2200,ERR=99,END=5)MIN(II) ,Q(II) ,DUM2,DUM3,DUM4 

C 

C*********************************************.~******************** 

C * a 

C INITIALIZE LOOP COUNTERS AND ERROR FLAGS * 

C * 

C****************************************************************** 

C 

5 ITEMP=O 

DX=Hl/FLOAT (NJ) 

NI=INT(H/DX+O. 5) 

DXX=H-DX*FLOAT(NI-1) 

DY=W/FLOAT (2*NJ) 

DTIME=60.O*TIMER/FLOAT(NTIME) 

e 

0 C 

DO 10 II=l,NI 

DO 10 JJ=l,NJ 

T( 11, JJ)=-1.0 

TP(I1, JJ) =T( 11, JJ) 

RMELT ( I I , J J ) = 0 . 0 

10 CONTINUE 

TlAVG=T (1,l) 0 

C 

c****************************************************************** 

C * 

C BEGIN CALCULATIONAL LOOPS ON LENGTH, TIVF, ANGLE * 

C * 

C****************************************************************** 

c 

DO 100 NO=l,NORBIT 

NPUT=O 

TIME=O. 0 

DO 100 N=l,NTIME 

C 

e 

IF(MIN(NPUT+l) .LE. TIME+O.lE-G)NPUT=NPUT+l 

C 

DO 20 JJ=l,NJ 

QNA(JJ) =Q(NPUT)/FLOAT(NJ) 

ROMELT (11, JJ) =MELT (11, JJ) 

TO (11, JJ) =T( 11, JJ) 

DO 20 II=l,NI 
20 

C 

30 

C 

C 

40 

C 

50 

C 

C 

C 

60 

70 

C 

80 

C 

CONTINUE 

ICOUNT=O 

IFLAG=O 

ICOUNT=ICOUNT+l 

CALL SODIUM(NPUT,THMELT,TKNA,TNA,U) 

DO 40 J=l,NJ 

CALL CONDUCT (BORL, DXX, SORL, TNA, U) 

CALL SOLVER 

CONTINUE 

DO 50 II=l,NI 

DO 50 JJ=l,NJ 

ERROR=ABS(T(II,JJ)-TP(II,JJ))/THMELT 

IF(ERR0R .GT. ERRMAX)IFLAG=2 

IF(ERR0R .GT. ERRMAX)ITEMP=ITEMP+l 

TP(I1, JJ) =T (11, JJ) 

CONTINUE 

IF(IC0UNT .GT. 1MAX)THEN 

WRITE (7,1000) ITEMP 

IFLAG=O 

ENDIF 

IF(1FLAG .GE. l)GO TO 30 

TIME=TIME+DTIME/GO.O 

IF(N0 .LT. NORB1T)GO TO 100 

WRITE (7,1050) 

WRITE (7,1250) 

WRITE (7,1100) TIME, TKNA 

DO 70 II=l,NI 

DO 60 JJ=l,NJ 

CONTINUE 

THETA(JJ)=T(II-(JJ-1) ,JJ)+THMELT 

JJJ=AMINO(II,NJ) 

WRITE(7,1200) (THETA(JJ) ,JJ=l,JJJ) 

CONTINUE 

WRITE (7,1050) 

WRITE (7,1250) 

DO 80 I=l,NI 

JJJ=AMINO (NJ, I) 

WRITE(7,1350) (RMELT(I-(JJ-l) ,JJ) ,JJ=l,JJJ) 

CONTINUE 

100 CONTINUE 

C 
B6 

0 

0 

e 

a 

................................................................ 

C * 

C INPUT AND OUTPUT FORMAT STATEMENTS * 

C * 

C*************************************************************** 

C * 

1000 FORMAT(lX,'FAILED TO CONVERGE AFTER ',16,' ITERATIONS.') 

1050 FORMAT(lX,' ') 

1100 FORMAT(lX,' Time = ', 1F6.1,' min. Sodium Temp. = ', 

1200 FORMAT( 10F7.1) 

, 1F6.1,' K.') 

1250 FORMAT(~X,'-------------------------------------------------/ ................................................. 

I +'--------------------' .................... 

+ ' , F7.4, 'KW/m2-K') 

1300 FORMAT(lX,'Avg. Gas Temp. = ',F6.1,'K Oiit. Gas temp.',FG.l,'K H = 

1350 FORMAT(lOF7.4) 

1400 FORMAT(lX,'INPUT FORMAT ERROR IN 1ST LINE') 

1401 FORMAT(lX,'INPUT FORMAT ERROR IN 2ND LINE') 

1402 FORMAT(lX,'INPUT FORMAT ERROR IN 3RD LINE') 

1403 FORMAT(lX,'INPUT FORMAT ERROR IN 4TH LINE') 

1404 FORMAT(lX,'INPUT FORMAT ERROR IN 5TH LINE') 

1405 FORMAT(lX,'INPUT FORMAT ERROR IN 6TH LINE') 

1406 FORMAT(lX,'INPUT FORMAT ERROR IN TRANSIENT SECTION') 

2100 FORMAT(I6,4110) 

2200 FORMAT(SF10.4) 

93 WRITE(*,1400) 

94 WRITE (*, 1401) 

95 WRITE(*,1402) 

96 WRITE(*,1403) 

97 WRITE(*,1404) 

98 WRITE (*, 1405) 

99 WRITE (*, 1406) 

STOP 

STOP 

STOP 

STOP 

STOP 

STOP 

STOP 

STOP 

END 
APPENDIX C 

PROGRAM FOR THE THERMAL ANALYSIS 

OF THE 

PEBBLE BED RECEIVER 
c1 

0 

0 

e 

C 

C 

0 

0 

e 

0 

e 

e 

e 

0 

PROGRAM MAIN 

C************************************::~******************* 

c********************************************************* 

C * 

C TEST COEFF * 

C * 

C * .......................................................... 

C*******************************************************~* 

C 

DIMENSION TL(15,10,20),RL(15,10,20) 

COMMON/BLOKF/I,J,L,NI,NJ,NL,ITIME,MODE 

COMMON/BLOKA/RMELT(15,lO),KLIQ,KSOL,K~AS,~ET,VOID,EPS, 

I DRM,RO,DENLIQ,CVLIQ,CVSOL,DENMET,CVMET,HFL, 

I SAPUV,PI,PR,REP,H,CVPURE,DENPURE,DENSOL, 

I ROMELT(15,lO) ,RPMELT(15,10) 

I TGAS0(50),GASMIX,TGASIN,DP,VF,DPHI,CPGAS, 

COMMON/BLOKB/DZ,DTIME,DR,QAZ,QZ,HEAT 

COMMON/BLOKD/TGAS(50),TO(l5,lO),T(l5,lO),AREA,GAS~S, 

I DENGAS, TP.( 15,lO) , TGASP (50) , TGAST (50,100) 

COMMON/BLOKG/KE(l5,10) ,KW(15,10) ,KN,KS,HWGAS(lO) ,HWSOL 

COMMON/BLOKJ/A ( 15) , B ( 15) , C ( 15) , D ( 15) 

COMMON/BLOKO/KPSOL, KPURE, K (3,3) 

COMMON/ BLOKX/FCTR 

COMMON/BLOKS/Q1DIR,QlIN,Q2DIR,Q2IN 

C 

REAL KPSOL, KSOL, KPURE, DFT, DF, KSTAG, KCOND, KRAD, KRO, FS , KWKF 

REAL KFSTAG,KFF,KWALL,KWRAD,K,KLIQ,KGAS,KMET,LENGTH 

REAL KN , KS , KE , KW 

INTEGER DTIME 

OPEN (8, FILE=' INPDATA. TXT' ) 

C 

READ(8,*)NI,NJ,NL,NTIME,DTIME,NORBIT 

READ (8, *) DP, DRM, RO, LENGTH, VF, SAPUV, SAPUVW 

READ (8, *) KLIQ, CVLIQ, DENLIQ, THMELT 

READ (8, *) KSOL, CVSOL, DENSOL, HFL 

READ (8, *) KMET, CVMET, DENMET, EPS 

READ(8,*)TGASIN,GASMAS,GASMIX 

READ(8,*)FCTR,ERRMAX,QlDIR,QlIN,Q2DIR,Q2IN,I~X 

C 

PI=3.14159 

DPHI=3.14159/NJ 

AREA=PI*RO**2 

DR=RO/(NI-1) 

DZ=LENGTH/NL 

OPEN(S,FILE='TEMP.TXT') 

DO 12 IORBIT=l,NORBIT 

DO 12 MODE=1,2 

IF (MODE .EQ. 2) NTIME=2160 

IF (MODE .EQ. 1) NTIME=3240 
c2 

a 

l 

a 

a 

a 

6 

a 

a 

7 

400 

a 

500 

10 

11 

2200 

40 

2216 

DO 12 L=l,NL 

DO 12 ITIME=l,NTIME/DTIME 

IF ((ITIME .EQ. l).AND.(MODE .EQ. 1) 

.AND. (IORBIT .EQ. 1) ) THEN 

TGAS (L) =TGASIN 

TGASO (L) =TGASIN 

DO 6 II=l,NI 

DO 6 JJ=l,NJ 

T(I1, JJ) =1218.0 

TO( 11, JJ) =1218.0 

TP(II,JJ)=1218.0 

ROMELT (11, JJ) =O. 0 

MELT( 11, JJ) =O. 0 

CONTINUE 

TGASO (L) =TGAST (L, NTIME/DTIME) 

TGAS (L) =TGASO (L) 

TGASP (L) =TGASO (L) 

DO 7 II=l,NI 

DO 7 JJ=l,NJ 

ELSEIF (ITIME .EQ. 1) THEN 

TO (11, JJ) =TL( 11, JJ, L) 

T(I1, JJ) =TO (11, JJ) 

TP( 11, JJ) =TO (11, JJ) 

ROMELT (11, JJ) =RL (11, JJ, L) 

RMELT(II,JJ)=ROMELT(II,JJ) 

CONTINUE 

ENDIF 

JKJ=O 

HEAT=O . 0 

DO 11 J=l,NJ 

CALL TGAS2 (KGAS , H , PR, REP) 

DO 10 I=l,NI 

CALL TCON 

CALL TCOEFF 

CONTINUE 

CALL TSOLVER 

CALL TMELT (THMELT) 

CONTINUE 

CALL TGASTEMP 

FORMAT (lX, 'TGAS ( 

DO 30 J=l,NJ 

DO 30 I=l,NI 

,I1, ' ) =' , F7.2) 

IF (ABS(T(1,J)-TP(1,J)) .GE. ERRMAX) THEN 

DO 40 JJ=l,NJ 

DO 40 IT=l,NI 

CONTINUE 

JKJ=JKJ+l 

WRITE(*,2216)IORBIT,MODE,L,ITIME,JKJ 

FORMAT('+',413,5X,14) 

TP( 11, JJ) =T (11, JJ) 

a 
c3 

~. 

I IF (JKJ .GT. IMAX) THEN I WRITE (5,7000) JKJ 

7000 FORMAT(lX,'FAILED TO CONVERGE AFTER',I4,' 

& ITERATIONS ' ) 

0 

30 

+ 

29 

28 

0 

700 

3500 

& 

0 

2000 

GO TO 700 

ENDIF 

GO TO 400 

ENDIF 

CONTINUE 

SUMT2=O. 0 

SUMR2=O. 0 

DO 28 I=l,NI 

SUMTl=O. 0 

SUMRl=O. 0 

DO 29 J=l,NJ 

CVPURE=RMELT (I, J) *CVLIQ+ 

SUMTl=SUMTl+(T(I,J)-TO(I,J))*CVPURE 

IF (I .EQ. 1) GO TO 29 

SUMRl=SUMRl+ (WELT (I, J) -ROMELT (I, J) ) 

(1.0-RMELT (I, J) ) *CVSOL 

CONTINUE 

RS=(NI-I) *RO/ (NI-1) 

RN=RS+RO/ (NI-1) 

VOLUME=(DPHI/2)*(RN-RS)*(RN+RS)*DZ 

PRODT=VOLUME *SUMTl 

PRODR=VOLUME*SUMRl 

SUMT2 =SUMT2 +PRODT 

SUMR2=SUMR2+PRODR 

CONTINUE 

QLAT=2*(1-VF)*DENPURE*SUMR2*HFL/DTIME 

QSENS=2*(1-VF)*DENPURE*SUMT2/DTIME 

IF (L .EQ. 1) THEN 

TGASB=TGASIN 

ELSE 

TGASB=TGAST(L-1,ITIME) 

ENDIF 

QGAS=GASMAS*CPGAS* (TGAS (L) -TGASB) 

TGASO (L) =TGAS (L) 

TGAST (L, ITIME) =TGAS (L) 

WRITE (5, *) 

WRITE (5, *) 

WRITE (5, *) 

WRITE(5,3500)IORBIT,MODE,ITIME,L 

FORMAT(1X,'IIIIII ORBIT # ',Il,' 111111 MODE ',Il,' 

WRITE (5, *) 

DO 15 I=l,NI 

111111 TIME=',Il,' 111111 L=',Il,' 111111') 

WRITE(5,2000) (T(1, J) , J=NJ, 1, -1) 

FORMAT(lX,lO(F6.l,lX)) 
a 

c4 

0 

~* 15 CONTINUE 

WRITE (5, *) 

DO 18 I=l,NI 

WRITE ( 5 ,2 10 0 ) ( RMELT ( I, J ) , J=NJ, 1 , - 1 ) 

2100 FORMAT(lX,lO(F6.3,1X)) 

18 CONTINUE 

WRITE (5, *) 

WRITE(5,8000)KGAS,H,PR,REP 

a 

8000 FORMAT(lX,'KGAS=',F9.7,3X,'H=',F6.4,3X,'PR=',F6.4, 

SI 3X,'REP=',F8.3) 

WRITE ( 5,2200) L, TGAS (L) 

HEAT=HEAT*DZ 

WRITE(5,8001)QZ,HEAT 

WRITE (5, *) 

WRITE(5,8002)QGAS,QLAT,QSENS 

a 

8001 FORMAT(lX,'HEAT INPUT FRAC. =',F8.5,3X,'HEAT=',F8.5) 

8002 FORMAT(lX,'QGAS=',F7.5,5X,'QLAT=',F7.5,5X, 

& 'QSENS=' , F7.5) 

DO 16 I=l,NI 

DO 16 J=l,NJ 

a 

TO(1, J) =T(I, J) 

TP( I, J) =T (I, J) 

TL( I, J, L) =T (I, J) 

ROMELT (I, J) =WELT (I, J) 

RL( I, J, L) =RMELT (I, J) 

16 CONTINUE 

12 CONTINUE 

5000 END 

0 

a 

a 

a 
c5 

e 

e 

0 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c********************************************************** 

C************************************************~********* 

C * 

C SUBROUTINE GASPROP * 

c * 

C This program calculates the effectiva thermal properties* 

C and convective coefficients for the gas coolant. * 

C The gas coolant is assumed to be a mixture of He and Xe.* 

C Properties of these two constituents are included as * 

C curve fits with the corresponding coefficients included * 

C as data in the subroutine. * 

C * c********************************************************** 

C********************************************************** 

C 

c********************************************************** 

C * 

C REFERENCES * 

C M. R. Vanco, ffAnalytical Comparison of Relative Heat * 

C Transfer Coefficients & Pressure Drops of Inert Gases * 

C and Their Binary Mixturesff, Lewis Research Center * 

C * 

C A. J. Chapman ffFundamentals of Heat Transfer" * 

C Macmillan, 1978, p 329. * 

C * 

C********************************************************** 

C 

C********************************************************** 

C * 

C EQUATIONS * 

C K1 = K11 + K12*(TGAS - K13)**K14 * 

C K2 = K21 + K22*(TGAS - K23)**K24 * 

C * 

C VISCl = VISCll + VISC12*(TGAS - VISC13)**VISC14 * 

C VISC2 = VISC21 + VISC22*(TGAS - VISC23)**VISC24 * 

C * 

C********************************************************** 

C 

SUBROUTINE TGAS2(KGAS,H,PR,REP) 

REAL K,Kl,K2 ,Kll,K12 ,K13,K14 ,K21,K22 ,K23 ,K24 ,MW,MWl,MW2, 

COMMON/BLOKD/TGAS ( 50 ) , TO ( 15,lO) , T ( 15,lO ) , AREA, GASMAS, , KGAS 

I TGAS0(50),GASMIX,TGASIN,DP,VF,DPHI,CPGAS, 

f DENGAS,TP(15,10) ,TGASP(5O) ,TGAST(50,100) 

COMMON/BLOKF/I,J,L,NI,NJ,NL,ITIME,MODE 

COMMON/BLOKG/KE(l5,10) ,KW(15,10) ,KN,KS,HWGAS(lO) ,HWSOL 

DATA GAMMA, PI, RUNIV, MW1, MW2/ 

DATA Kll,K12,K13,K14/0.1580E-3,0.9086E-6,333.333,0.794/ 

/1.66,3.14159,8.31434,4.0,131.0/ 
C6 

0 

0 

e 

C 

C 

C 

C 

C 

C 

DATA K21,K22,K23,K24/0.640E-5l0.0454E-6l~~~.3~3l~~800/ 

DATA V1SC11,V1SC12,V1SC13,V1SC14/ 

DATA VISC21,VISC22,VISC23,~1SC2~/ 

F1 (XI Y) = (1.0+2.41* (X-Y) * (X-O.l42*Y) / (X+Y) **2) 

/.2158E-4,.1022E-6,333.333,.825/ 

/0.253E-4,0.2452E-6,333.333,0.77/ 

F2 (Xl, X2 , X3 , X4 , X5, X6) = 

=( (1.0+(x1/x2)**0.5*(x3/x4)**0.25)**2) ’ 

/ (2.0*2.0**0.5* (l.O+X5/X6) **O. 5) 

C 

C********************************************************** 

C * 

C * C********************************************************** 

C FIND THE THERMAL CONDUCTIVITY & VISCOSITY OF PURE GASES * 

C 

Kl=Kll+K12* (TGAS (L) -K13) **K14 

K2=K21+K22*(TGAS(L)-K23) **K24 

VISC1=VISCll+VISC12*(TGAS(L)-VISCl3)**VISC14 

VISC2=VISC21+VISC22* (TGAS (L) -VISC23) *jTIISC24 

C 

c********************************************************** 

C * 

C * C********************************************************** 

* C FIND THE WEIGHTING PARAMETERS FROM REFERENCE 

C 

PHI12=F2 (V1SC1,V1SC2,MW2,MW1,MW1,MW2) 

PHI21=PHI12*(VISC2*MWl)/(VISCl*MW2) 

PS112=F2(K1,K2,MW1,MW21MWl,MW2)*Fl(MWllMW2) 

PSI21=F2 (K2,Kl,MW2,MWl,MW2,MWi) *F1 (MW2,MWl) 

C c********************************************************** 

C * 

C * C********************************************************** 

/ 

* 

* C SOLVE FOR THE THERMAL CONDUCTIVITY 

C AND VISCOSITY OF GAS MIX 

C 

VISC=VISCl/(l.O+PHI12*(l.O-GASMIX)/GASMIX)+ 

KGAS=Kl/(l.O+PSI12*(1.0-GASMIX)/GASMIX)+ 

VISC2/(l.O+PHI21*(GASMIX/(l.O-GASMIX))) 

/ K2/(l.O+PSI21*(GASMIX/(l.O-GASMIX))) 

C****t********+*****+*******f**t************hk******************* * 

* C SOLVE FOR SPECIFIC HEAT OF GASMIX 

C * ........................................................... 

C 

MW=GASMIX*MWl+ (1.0-GASMIX) *MW2 
c7 

0 

0 

0 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

-C 

C 

C 

C 

CPGAS=GAMMA*RUNIV/ ( (GAMMA-1.0) *MW) 

C 

C********************************************************** 

C * 

C CORRELATION RECOMMENDED FOR GAS FLOW IN A BED OF SPHERES* 

C 

C 

C 

0.575 

(e?) = (ej) = 2.06 (Re)- 

h m 

* 

* 

* 

C * 

C j are the Colburn j factors, and e is the void fraction * 

C * 

C REFERENCE * 

C Incropera, Frank P. , "Fundamentals of Heat and Mass * 

C Transfer", 2nd Edition, John Wiley & Sons Inc., 1981. * 

C * 

C********************************************************** 

C 

C p 292,352. * 

PR=CPGAS*VISC/(KGAS) 

REP=GASMAS*DP/(VISC*AREA) 

H= ( 2.06/VF) * ( KGAS/DP) *PR** ( 1.0/ 3 ) *REP**O .4 2 5 

DENGAS=11.22 

RETURN 

END 

C 

a 

e 
C8 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

a C 

a C 

a 

C 

C*********************************************************** 

C*********************************************************** 

C * 

C SUBROUTINE COEFF * 

C * 

C This program develops finite difference~matrix for * 

C conduction analysis in cylindrical coordinates. * 

C * 

C*********************************************************** 

C 

C*****************+********************’.k******************* 

SUBROUTINE TCOEFF 

COMMON/BLOKA/RMELT(15,10) ,KLIQ,KSOL,KGAS,KMET,VOID,EPS, 

I DRM,RO,DENLIQ,CVLIQ,CVSOL, DENMET,CVMET,HFL, I SAPUV,PI,PR,REP,H,CVPURE,DENPURE,DENSOL, I ROMELT(15,lO) ,RPMELT(15,10) 

I TGASO(SO),GASMIX,TGASIN,DP,VF,DPHI,CPGAS, 

I DENGAS,TP(15,10) ,TGASP(50) ,TGAST(50,100) 

COMMON/BLOKB/DZ,DTIME,DR,QAZ,QZ,HEAT 

COMMON/BLOKD/TGAS(50),TO(l5,lO),T(l5,lO),AREA,GASMAS, 

COMMON/BLOKF/I, J,L,NI,NJ,NL, ITIME,MODE 

COMMON/BLOKG/KE(l5,lO),~(l5,lO),KN,K~,HWGAS(l~),HWSOL 

COMMON/BLOKJ/A(15) ,B(15) ,C(15) ,D(15) 

C 

INTEGER DTIME 

REAL KN,KS,KE,KW 

C 

C 

a 

0 

C 

a 

CVPURE=RMELT(I,J)*CVLIQ+(l.O-RMELT(I,J))*CVSOL 

DENPURE=DENLIQ 

DR=RO/ (NI-1) 

RS=(NI-I) *DR 

IF(I.EQ.l) THEN 

DR=DRM 

DRS= ( (RO/ (NI-1) ) +DR) / 2 

GO TO 100 

ELSEIF (I.EQ.2) THEN 

DRN= (DR+DRM) /2 

DRS=DR 

GO TO 100 

RS=O. 0 

DRN=DR 

GO TO 100 

GO TO 90 

ELSEIF (I. EQ. NI) THEN 

ELSE 

ENDIF 

90 DRN=DR ‘‘4 
c9 

a 

DRS=DR 

RE=RS+DR/2 

RW=RE 

DV=(DPHI/2) *(RN-RS) *(RN+RS) 

PF=l-VF 

AE=KE (I, J) *DR/ (RE*DPHI) 

AW=KW(I,J) *DR/(RW*DPHI) 

IF (1.EQ.NI) THEN 

100 RN=RS+DR 

C 

C( I) =o. 0 

GO TO 250 

ENDIF 

AWALL=PI *RO/ NJ 

C(I)=-(KS*RS*DPHI/DRS) 

0 

C 

IF (I.EQ.1) GO TO 290 

C 0 

0 

250 A(I)=-(KN*RN*DPHI/DRN) 

BC=H*SAPUV*DV*TGAS(L)+PF*DENPURE*CVPURE*DV*TO(I,J)/DTIME 

B(I)=AE+AW-C(I)-A(I)+H*SAPUV*DV+PF*DENPURE*CVPURE*DV/DTIME 

D (I) =AE*T (I, J-1) +AW*T (I, J+1) +BC 

GO TO 300 

/+PF*DENPURE*HFL*DV*(ROMELT(I,J)-RMELT(I,J))/DTIME 

C 

290 C( I) =C( I) 

B (I) =AE+AW-C (I) +DENMET*CVMET*DV/DTIME+HWGAS (J) *AWALL 

CALL TSOLAR 

BC=(DENMET*CVMET*DV*TO(I,J)/DTIME)+QAZ+HWGAS(J)*AWALL* 

*TGAS (L) 

D(I)=AE*T(I,J-l)+AW*T(I,J+l)+BC 

A( I) =O. 0 

C 

300 RETURN 

END 0 

a 
c10 

C ......................................................... C C************************************~~******************* 

C C * 

C C SUBROUTINE CONDUCT * 

C C * 

e C C This program calculates the interface conductivities to* 

C C the north, south, east, and west of the element (1,J). * 

C C * 

C C********************************************************* 

C C********************************************************* 

C C 

C 

SUBROUTINE TCON 

COMMON/BLOKF/I, J, L, NI I NJ, NL, ITIME,MODE 

COMMON/BLOKA/RMELT( 15,lO) ,KLIQ,KSOL,KGAS,KMET,VOID,EPS, 

/ DRM,RO,DENLIQ,CVLIQ,CVSOL,DENMET,CVMET,HFL, 

I SAPUV,PI,PR,REP,H,CVPURE,DENPURE,DENSOL, I ROMELT(15,lO) ,RPMELT(lS,lO) 

I TGASO(5O) ,GASMIX,TGASIN,DP,VF,DPHI,CPGAS, 

I DENGAS,TP(15,10) ,TGASP(SO) ,TGAST(50,100) 

COMMON/BLOKD/TGAS (50) , TO ( 15,lO) , T ( 15,lO) , AREA, GASMAS, 

COMMON/BLOKG/KE(l5,10) ,KW(15,10) ,KN,KS,HWGAS(lO) ,HWSOL 

COMMON/BLOKO/KPSOL, KPURE, K (3,3) 

0 C 

e 

REAL KPSOL, KPURE , DFT, DF I KSTAG, KCOND, KRAD, KRO , FS , KWKF 

REAL KFSTAG,KFF,KWALL,KWRAD,K,KLIQ,KSOL,KGASt~ET,KWG~D 

REAL KN, KS , KE , KW, NUM 

e 

DO 110 II=2,3 

IF ((I.EQ.1) .AND. (II.NE.3)) GO TO 110 

IF ((I.EQ.NI).AND.(II.EQ.3)) THEN 

KN=KS 

KS=O. 0 

GO TO 155 

ENDIF 

C 

DO 120 JJ=2,3 

C 

IF (I1 .EQ. 2) THEN 

IF (JJ .EQ. 2) THEN 

ELSE 

VOID=(l.O-DENLIQ/DENSOL)*(l.O-RMELT(1,J)) 

VOID=(l.O-DENLIQ/DENSOL)*(l.O-RMELT(I,J+l)) 

ENDIF 

ELSE 

VOID=(l.O-DENLIQ/DENSOL)*(l.O-WELT(I+1,J)) 

0 ENDIF 

C 

IJ=J 
c11 

0 

0 

0 

0 

0 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

IF (II.EQ.2) IJ=J-2+JJ 

IF ((J.EQ.NJ).AND.(JJ.EQ.3)) THEN 

KW( I, J) =O. 0 

GO TO 105 

ENDIF 

C 

.......................................................... 

C * 

C CALCULATE THE EFFECT OF VOID ON CONDUCTIVITY OF SOLID * 

C USE CORRELATION BY MEREDITH & TOBIAS TO CORRECT * 

C CONDUCTIVITY * 

C * .......................................................... 

C 

KPSOL~KSOL*(2.0*(1.0-VOID)+O.6135*VOI~**2.333-1.6* 

*VOID**3.333)/ 

/(2.O+VOID+0.6135*VOID**2.333-0.6975*VOID**3.333) 

C 

C********************************************************* 

C FIND CONDUCTIVITY OF PURE MATERIAL * 

c********************************************************* 

C 

JI=I 

IF (I1 .EQ. 3) JI=I+l 

KPURE=l. O/ ( (RMELT (JI , IJ) /KLIQ) + ( ( 1. O-RMELT (JI , IJ) ) /KPSOL) ) 

HRS=O. 227 * (EPS/ (2. O-EPS) ) * (T (I, J) **3 ) *lOE-9 

HRV~.227*T(I,J)**3*1OE-9/(l+VF*(l.O-E~~)/(2*(l-VF)*EPS)) 

C 

C*******~************************************************* 

C * 

C FIND THE EFFECTIVE CONDUCTIVITY OF THE BED USING METHOD* 

C * 

C OF SUMMED RESISTANCE OVER VARIOUS HEAT TRANSFER PATHS * 

C * 

C HANDBOOK OF HEAT TRANSFER APPLICATIONS, * 

C 2nd EDITION, 1985, pp. 6-10,6-18. * 

C********************************************************* 

C * 

DFT=((KPURE-KGAS)/KPURE)**2 

DF=( (0.25*DFT) / (ALOG (KPURE/KGAS) -DFT) ) -KGAS/ (3*KPURE) 

NUM=((KPURE-KGAS)/KPURE)**2 

Dl=ALOG ( (KPURE/KGAS) -. 9128* ( (KPURE/KGAS) -1.0) ) 

D2=( (KPURE-KGAS)/KPURE) *. 0872 

PHI= ( ( 1 . O/ 12 ) *NUM/ (Dl-D2 ) ) - ( 2.0/3 ) * (KGAS/KPURE) 

D=(l/PHI)+(DP*HRS/KGAS) 

KSTAG=VF* (l+HRV*DP/KGAS) + (1-VF) / 

KCOND=(KSTAG+O.ll*PR*REP)*KGAS 

K (11, JJ) =KCOND 

/ ( ( 1.. O/D) + ( 2 . O/ 3 ) *KGAS/KPURE ) 
c12 

105 IF (II.EQ.3) GO TO 110 

120 CONTINUE 

110 CONTINUE 

a 

C C * 

C C********************************************************* 

C C FIND WALL-TO-BED CONDUCTIVITY USING THE METHOD OF * 

C C * 0 

C C "EFFECTIVE PATH LENGTH" * 

C C * 

C C IIHANDBOOK OF HEAT TRANSFER APPLICATIONS 2nd EDITION, * 

C C 1985, pp. 6-10,6-18. * 

C .......................................................... 

C C * 0 

IF (I.EQ.1) THEN 

190 KWKF=2.0*0.7+(1.0-0.7)/(DF+(KGAS/(3.0*KPURE))) 

KFSTAG=l.O/KWKF-(O.5*KGAS/KSTAG) 

KFF=l.O/KFSTAG+O.O54*PR*REP 

KWALL=KFF*KGAS/DP 

DEN=(VF* (1. O-EPS) *lOE+6)/ (2.0* (1-VF) *EPS) 

KWGRAD=(O.227*T(I, J) **3/(DEN+1.0) )/lo00 

HWGAS (J) =KWALL 

HWSOL=( (0.227* (EPS/ (2. O-EPS) ) * (T(1, J) **3*10E-6) )/1000) + 

+ KWRAD 

0 KE (I, J) =KMET 

IF (J .EQ. 1) KE(I,J)=O.O 

KW( I, J) =KMET 

IF (J .EQ. NJ) KW(I,J)=O.O 

FS=DRM/ (DRM+RO/ (NI-1.0) ) 

K(3,2) =KWKF*KGAS+O. ll*PR*REP*KGAS+IIWSOL* 

* (RO/ (NI-1.0) +DRM) /2.0 

KS=l.O/ (FS/KMET+ (1. O-FS) /K( 3,2) ) 

KN=O. 0 

GO TO 200 

0 

ENDIF 

KN=KS 

KS=l.O/ (0.5/K(2,2) +O. 5/K( 3,2) ) 

C 

0 

C 

155 IF (J.EQ.l) THEN 

KE(1, J) =O. 0 

GO TO 195 

KE ( I, J) =KW ( I, J-1) 

ELSE 

ENDIF 

IF ((KW(I,J).EQ.O) .AND. (J .EQ. NJ)) GO TO 200 

195 KW (I, J) =l. O/ (0.5/K( 2,2) +O. 5/K(2,3) ) 

200 RETURN 

END 

a 

C 

a 
e 

C13 

0 

0 

m 

a 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C********************************************************** 

C********************************************************** 

C * 

C SUBROUTINE GASPROP * 

C * 

C This program calculates the effective thermal properties* 

C and convective coefficients for the gas coolant. * 

C The gas coolant is assumed to be a mixture of He and Xe.* 

C Properties of the two constituents are included as curve* 

C fits with the corresponding coefficients included as * 

C data in the routine. * 

C * 

C********************************************************** 

C********************************************************** 

C 

C********************************************************** 

C * 

C REFERENCES * 

C M. R. Vanco, ltAnalytical Comparison of Relative Heat * 

C Gasesand Their Binary Mixtures", Lewis Research Center* 

C * 

C Transfer Coefficients and Pressure Drops of Inert * 

C A. J. Chapman llFundamentals of Heat Transfer" * 

C Macmillan, 1978, p 329. * 

C * 

C********************************************************** 

C 

C********************************************************** 

C * 

C EQUATIONS * 

C K1 = K11 + K12*(TGAS - K13)**K14 * 

C K2 = K21 + K22*(TGAS - K23)**K24 * 

C * 

C VISCl = VISCll + VISC12*(TGAS - VISC13)**VISC14 * 

C VISC2 = VISC21 + VISC22*(TGAS - VISC23)**VISC24 * 

C * 

C********************************************************** 

C 

SUBROUTINE TGAS2(KGAS,H,PR,REP) 

REAL K,K1,K2,Kll,K12,K13,Kl4,K2l,K22,K23,K24,MW,MWl,MW2, 

COMMON/BLOKD/TGAS(50),TO(l5,lO),T(l5,lO),AREA,GASMAS, 

, KGAS 

I TGAS0(50),GASMIX,TGASIN,DP,VF,DPHI,CPGAS, 

I DENGAS,TP(15,1O),TGASP(5O),TGAST(50,100) 

COMMON/BLOKF/I , J, L, NI, NJ, NL, ITIME, MODE 

COMMON/BLOKG/KE(l5,10) ,KW(15,10) ,KN,KS,HWGAS(lO) ,HWSOL 

DATA GAMMA,PI,RUNIV,MWl,M2/ 

DATA Kll,K12,K13,K14/0.15~0E-3,0.9086E-6,333.333,0.794/ 

/1.66,3.14159,8.31434,4.0,131.0/ 
C14 

e 

a 

0 

e 

0 

a 

0 

DATA K21,K22,K23,K24/0.640E-5,0.0454E-6,333.333,0.800/ 

DATA VISC11,VISC12,VISC13,vISc14/ 

DATA VISC21,VISC22,VISC23,vIsc24/ 

F1 (XI Y)=( 1.0+2.41* (X-Y) * (X-O.l42*Y) / (X+Y) **2) 

F2 (Xl,X2 I X3 ,X4 I X5, X6) = 

/.2158E-4,.1022E-6,333.333,.825/ 

/0.253E-4,0.2452E-6,333.333,0.77/ 

= ( (l.O+ (Xl/X2) **O. 5* (X3/X4) **O. 25) **2)/' 

/ (2.0*2.0**0.5* (l.O+X5/X6) **O. 5) 

C 

........................................................... 

C * 

C FIND THE THERMAL CONDUCTIVITY C VISCOSITY OF PURE GASES * 

C * 

C********************************************************** 

C 

Kl=Kll+K12* (TGAS (L) -K13) **K14 

K2=K21+K22* (TGAS (L) -K23) **K24 

VISCl=VISCll+VISC12* (TGAS (L) -VISC13) **VISC14 

VISC2=VISC2l+VISC22*(TGAS(L)-VISC23)**VISC24 

C 

C********************************************************** 

C * 

C FIND THE WEIGHTING PARAMETERS FROM REFERENCE * 

C * 

C********************************************************** 

C 

PH112=F2(VISC1,VISC2~MW21MWlIMWl~MW2) 

PHI21=PHI12*(VISC2*MWl)/(VISCl*MW2) 

PS112=F2(K1,K2,MW1,MW21MWl,MW2)*Fl(MWlIMW2) 

PSI21=F2 (K2, Kl,MW2 ,MW1 I MW2 ,MW1) *F1 (MW2 ,MW1) 

C 

C*****************************************************~**** 

C * 

C SOLVE FOR THE THERMAL CONDUCTIVITY & VISCOSITY OF MIX * 

C * 

C********************************************************** 

C 

/ VISC2/(l.O+PHI21*(GASMIX/(l.O-GASMIX))) 

/ K2/(l.O+PSI21*(GASMIX/(l.O-GASMIX))) 

C * 

C SOLVE FOR SPECIFIC HEAT OF GASMIX * 

C * 

C********************************************************** 

C 

VISC=VISC1/(1.O+PHI12*(l.O-GASMIX)/GASMIX)+ 

KGAS=Kl/(l.O+PSI12*(l.O-GASMIX)/GASMIX)+ 

C***********X*************************'k******************* 

MW=GASMIX*MWl+(l.O-GASMIX)*MW2 

CPGAS=GAMMA*RUNIV/ ( (GAMMA-1.0) *MW) 
C15 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

0 

C 

0 

C 

C********************************************************** 

C * 

C CORRELATION RECOMMENDED FOR GAS FLOW IN A BED OF SPHERES* 

C 

C 

C 

0.575 

(ej) = (ej) = 2.06 (Re)- 

h m 

* * 

* 

C * 

C j are the Colburn j factors, and e is the void fraction * 

C * 

C REFERENCE * 

C Incropera, Frank P. , "Fundamentals of Heat and Mass * 

C Transfer", 2nd Edition, John Wiley & Sons Inc., 1981, * 

C * 

C********************************************************** 

C 

C p 292,352. * 

PR=CPGAS*VISC/ (KGAS) 

REP=GASMAS*DP/ (VISC*AREA) 

H=(2.O6/VF)*(KGAS/DP) *PR**(1.0/3.0)*REP**O.425 

DENGAS=11.22 

RETURN 

END 

C 

C 
e 

l 

i 

0 

a 

0 

a 

a 

a 

e 

a 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C1G 

c********************************************************** 

c********************************************************** 

C * 

C SUBROUTINE GASTEMP * 

C * 

C This program solves for the temperature of the working * 

C fluid. The energy equation has been discretized using * 

C the upwind scheme to calculate the convection term. * 

C Ref: Numerical Heat Transfer And Fluid Flow It, * 

C Suhas V. Patankar. * 

C * 

C********************************************************** 

........................................................... 

C 

SUBROUTINE TGASTEMP 

COMMON/BLOKA/RMELT(15,lO),KLIQ,KSOL,KGAS,~ET,VOID,EPS, 

I DRM,RO,DENLIQ,CVLIQ,CVSOL,DENMET,CVMET,HFL, 

I ROMELT(15,lO) ,RPMELT(15,10) 

I TGASO(50) ,GASMIX,TGASIN,CP,VF,DPHT,CPGAS, 

I DENGAS,TP(15,1O),TGASP(5O),TGAST(50,100) 

I SAPUV, PI, PR, REP, H, CVPURE, DENPURE, DENSOL, 

COMMON/BLOKD/TGAS(50) ,T0(15,10) ,T(15,10) ,AREA,GASMAS, 

COMMON/BLOKB/DZ , DTIME, DR, QAZ , QZ, HEAT 

COMMON/BLOKF/I,J,L,NI,NJ,NL,ITIME,MODE 

COMMON/BLOKG/KE(15,10) ,KW(15,10) ,KN,KS,HWGAS(lO) ,HWSOL 

REAL KGAS 

INTEGER DTIME 

DR=RO/ (NI - 1 ) 

QWALLl=O.O 

QWALL2=O. 0 

QBEDl=O. 0 

QBED2=O. 0 

AWALL= (PI*RO/NJ) *DZ 

DO 11 JQ=l,NJ 

DO 11 IQ=2,NI 

RS=(NI-IQ) *DR 

RN=RS+DR 

DV=(DPHI/2)*(RN-RS)*(RN+RS)*DZ 

IF (IQ .EQ. 2) THEN 

QWALLl=QWALL1+2 *HWGAS (JQ) *AWALL*T ( 1, JQ) 

QWALL2=QWALL2+2*HWGAS(JQ)*AWALL 

ENDIF 

QBEDl=QBED1+2*H*SAPUV*DV*T(IQ,JQ) 

QBED2=QBED2+2 *H*SAPUV*DV 

11 CONTINUE 
a 

C17 

e 

C 

VOLUME=(PI*RO**2) *DZ 

AB=GASMAS*CPGAS 

B=(QWALLl+QBEDl)+DENGAS*VF*CPGAS*TGASO(L)*VOLUME/DTIME 

A=AB +(DENGAS*VF*CPGAS*VOLUME/DTIME)+QWALL2+QBED2 

IF (L.EQ.l) THEN 

TB=TGASIN 

ELSE 

TB=TGAST(L-1,ITIME) 

ENDIF 

TGAS (L) = (AB*TB+B) /A 

RETURN 

END 

a C 

a C 

a 

a 

a 
0 

C18 

e 

C C 

C C********************************************************** 

C C**************************************k******************* 

C C * 

C C SUBROUTINE MELT * 

C C * 

C C * 

C C * 

C C * 

C C * 

C C********************************************************** 

C C********************************************************** 

C C 

C 

C 

SUBROUTINE TMELT(THMELT) 

COMMON/BLOKF/I,J,L,NI,NJ,NL,ITIME,MODE 

COMMON/BLOKA/RMELT(15,lO),KLIQ,KSOL,KGAS,~ET,VOID,EPS, 

I DRM,RO,DENLIQ,CVLIQ,CVSOL,DENMET,CWET,HFL, 

I SAPUV,PI,PR,REP,H,CVPURE,DENPURE,DENSOL, 

ROMELT (15,lO) , RPMELT (15,lO) 

I TGASO(50) ,GASMIX,TGASIN,DP,VF,DPHI,CPGAS, 

I DENGAS,TP(15,10),TGASP(50),TGAST(50,100) 

I 

COMMON/BLOKD/TGAS(50),TO(l5,lO),T(l5,lO),AREA,GAS~S, 

COMMON/BLOKB/DZ,DTIME,DR,QAZ,QZ,HEAT 

COMMON/BLOKX/FCTR 

DO 10 II=l,NI 

C 

T ( 11, J) =T ( 11, J) -THMELT 

10 CONTINUE 

C 

DO 11 I=2,NI 

TN=T(I-1,J) 

IF (I .EQ. NI) THEN 

ELSE 

ENDIF 

WELT( 1, J) =l. 0 

DMELT=FCTR*T(I,J)*CVPURE/HFL 

IF ( (DMELT . LT. 0.0) .AND. (QAZ . LT. 0.0) ) THEN 

ELSEIF (I .EQ. 2) THEN 

C 

TS= (T (I, J) ) 

TS=T (I+1, J) 

C 

C 

IF (T(1,J) .LE. 0.0) RMELT(I,J)=RMELT(I,J)+DMELT 

IF ((T(1,J) .GE. 0.0) .AND. (QAZ .GE. 0.0)) THEN 

ENDIF 

RMELT (I, J) =WELT (I, J) +DMELT 

ELSEIF (T(1,J) .GE. 0.0) THEN 
c19 

RMELT (I, J) =WELT (I, J) +DMELT 

RMELT ( I, J) =O . 0 

ELSE 

ENDIF 

C 

RMELT ( I, J) =mxi ( 0.0, RMELT ( I, J) ) 

RMELT(I,J)=AMINl(l.O,RMELT(I,J)) 

11 CONTINUE 

C 

DO 13 I=l,NI 

T (I, J) =T (I, J) +THMELT 

13 CONTINUE 

RETURN 

END 

C 
t 

c2 0 

B 

B 

B 

D 

B 

B 

D 

B 

D 

............................................................... 

C*************************************************************** 

c * 

C SUBROUTINE SOLAR * 

C * 

C This program calculates the solar incidence on * 

C the outside perimeter of the containment ring * 

C q" = qI1 (direct) *cos (angle) + q" (indirect) * C * 

C * 

C*************************************************************** 

C*************************************************************** 

C 

C 

SUBROUTINE TSOLAR 

DIMENSION CS (21) 

REAL LENGTH 

COMMON/BLOKF/I I J , L, NI , NJ, NL, ITIME , MODE 

COMMON/BLOKB/DZ,DTIME,DR,QAZ,QZ,HEAT 

COMMON/BLOKS/Q1DIR,QlIN,Q2DIR,Q2IN 

DATA CS(1) ,CS(3) ,CS(5) ,CS(7) ,CS(9) ,CS(11) ,CS(13)/ 

DATA CS(15) ,CS(17) ,CS(19) ,CS(2) ,CS(4) ,CS(6) ,CS(8)/ 

DATA CS(10) ,CS(12) ,CS(14) ,CS(16) ,CS(18) ,CS(20) ,CS(21)/ 

DATA PI/3.14159/ 

DATA QZ1,QZ2,QZ3,Q24/24.28,0.08,0.75,0.26244/ 

/ 1.0,6*2.0/ 

/3*2.0,4*4.0/ 

/6*4.,1./ 

C 

C * 

C QZ USES SIMPSON'S 1/3 INTEGRATION TO FIND THE * 

C INTETRATED AVERAGE VALUE OF THE AXIAL POWER CURVE * 

C OVER THE INTERVAL DZ * 

C * 

C*************************************************************** 

C 

C***t f** ****************X*t***************~k******************* 

LENGTH=O .9 14 4 

POS=DZ*L 

XL=(POS-DZ)/LENGTH 

QZ=O. 0 

C 

DO 15 II=1,21 

IF(I1 .GT. 1) XL=XL+O.O5*DZ/LENGTH 

QZZ=QZl*XL*EXP(l.O-(XL/QZ2)**QZ3)+QZ4*(XL)**O.5 

QZ=QZ+CS(II)*(QZZ)/60.0 

15 CONTINUE 

C 

c*************************************************************** 

C * 
e 

e 

c2 1 

a 

e 

C FIND THE PROJECTION OF THE DIRECT SOLAR INCIDENCE * 

C * 

C*************************************************************** 

C 

ANGLEl=PI*FLOAT(J)/FLOAT(NJ) 

a ANGLE2=PI*FLOAT(J-l) /FLOAT (NJ) 

C 

IF (J+J-2 .GT. NJ) THEN 

ELSEIF (J+J .GT. NJ)THEN 

ELSE 

ENDIF 

SINANGLE=O.O 

SINANGLE=l. O-SIN (ANGLE2) 

SINANGLE=SIN(ANGLEl) -SIN(ANGLE2) 

C 

C*************************************************************** 

C * 

C * 

C******************************************~******************** 

C 

e C ADD DIRECT AND INDIRECT INCIDENT SOLAR RADIATION * 

IF (MODE .EQ. 1) THEN 

QDIR=QlDIR 

QINDIR=QlIN 

QDIR=Q2DIR 

QINDIR=Q2IN 

ELSE 

ENDIF 

DPHI=PI/NJ 

QAZ= (QDIR*SINANGLE + QINDIR*DPHI) *ROCP..V*QZ 

IF(J .EQ. l)HEAT=O.O 

HEAT=HEAT+2.0*QAZ 

RETURN 

e ROCAN=0.04191+0.00127 

C 

e END 

e 

e 

a 
e 

c22 

I

a 

e 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C********************************************************** 

C********************************************************** 

C SUBROUTINE SOLVER * 

C This program solves a tridiagonal matrix using an * 

C elimination method * 

C * 

C 

C 

[A]T(I-1) + [B]T(I) + [C]T(I+l) = [D] * 

* 

C*************************************',.~******************* 

C********************************************************** 

C 

SUBROUTINE TSOLVER 

C0MM0N/BL0KD/TGAS(50),T0(15,10),T(15,10),AREA,GASMAS 

I TGASO(50),GASMIX,TGASINIDP,VF,DPHI,CPGAS, I DENGAS,TP(15,1O),TGASP(50),TGAST(50,100) 

COMMON/BLOKJ/A ( 15) , B ( 15) , C ( 15) , D ( 15) 

COMMON/BLOKF/I,J,L,NI,NJ,NL,ITIME,MODE 

C 

C********************************************************** 

C * 

C PROCEED WITH GAUSSIAN ELIMINATION TO DEVELOP * 

C AN UPPER DIAGONAL MATRIX * 

C * 

C********************************************************** 

C 

DO 1 I=l,NI-l 

R=A ( I+ 1 ) / B ( I ) 

B (I+1) =B (I+1) -R*C (I) 

D(I+l)=D(I+l) -R*D(I) 

1 CONTINUE 

C 

C * 

C BEGIN BACK SUBSTITUTION * 

C * 

C************************************************************ 

C 

C***********X*************************kk******************* 

T (NI , J) =D (NI) /B (NI) 

DO 2 I=NI-1,1,-1 

T( I, J)=( D(1) -C(I) *T( I+1, J) )/B(I) 

2 CONTINUE 

RETURN 

END 

e 
lo 

e 

APPENDIX D 

PROGRAM FOR THE THERMAL ANALYSIS 

OF THE 

CAVITY RECEIVER 

a 

e 
D1 

e 

0 

0 

0 

e 

program main 

common/bloka/a (15) , b (15) , c (15) ; d (15) 

common/blokb/ni, ne (10) , j , t (15,lO) 

common/blokc/thmelt,tgasin,gasmas,cpgas,ntu,deltp, 

common/blokd/tgasp,gasmix,ro(l5),ht,w,dxx,p 

common/bloke/rmelt(l5,lO),denl~q,densol,walls,cvmet, 

+ denmet,kmet,kliq,ksol,cvliq,cvsol 

common/blokf/dz,dtime,rpo(lO),qaz,hfl,romelt(l5~,lO),cl,wf, 

+ to(15,10),rpme1t(15,10),dx,dy,th,qt,1e(15), 

+ ap00(15,lO),apon(l5,10) 

common/blokg/k,kn,ks,ke,kw,den,cv 

dimension cs(2l),q(20),tp(l5,lO),thgas(20,20), 

real ntu,length,kl~q,ksol,kmet,k,kn,ks,ke,kw,den,cv,le 

data cs(1) ,cs(3) ,cs(5) ,cs(7) ,cs(9) ,cs(ll) ,cs(13)/ 

data cs(15) ,cs(17) ,cs(19) ,cs(2) ,cs(4) ,cs(6) ,cs(8)/ 

C 

+ deltdz 

+ theta (20,lO) , thn (15,lO) , rmel(15,lO) , nn (10) 

C 

/l. 0,6*2.0/ 

/3*2.0,4*4.0/ 

data cs(l0) ,cs(12) ,cs(14) ,cs(16) ,cs(18) ,cs(20) ,cs(21)/ 

/6*4. , 1 ./ 

open (8,file='ginput.txtg) 

open (7,file='goutput.txt1) 

read(8,*)imax,norbit,nj,nl,ntime 

read(8,*)length,ht,w,th,sl 

read (8, *) gasmas, borl, cl, wf, tr 

read(8,*)cvliq,denliq,kliq,walls,thmelt 

read(8,*)cvsol,densol,ksol,errmax,hfl 

read(8,*)cvmet,denmet,kmet,gasmix,timer 

read(8, *) qzl, qz2, qz3, qZ4, qz5,qz6 

read (8, *) qd, qn, tgas, cor, ft 

p = 3.14 

ni = nj 

tgasp = tgas 

dy = ht/nj 

dz = length/float(nl) 

dtime = 60.0*timer/float(ntime) 

qtime = 54.0*60.0/dtime 

do 2 j = 1,111 

ro(j) = ni*dy 

ne(j) = ni-j+2 

C 

C**************************************************************** 

C 

2 continue 

C 

do 100 1 = 1,nl 

w = p*((nl-l+.5)*dz/sl+tr)/cor 
D2 

IO 

0 

0 

0 

C 

C 

C 

dx = w/nj 

dxx = (dx*dx+dy*dy) **O. 5 

do 10 jj = 1,nj 

do 10 ii = 1,15 

if (ii .gt. ne(]])) then 

else 

endif 

t(ii,jj) = 0.0 

t(ii,jj) = -O.OOOO~ 

tp(ii, jj) = t(ii, jj) 

rmelt(ii,jj) = 0.0 

apon(ii,jj) = 0.0 

pos = l*dz 

xl = pos/length 

qzz = qzl+qz2*xl+qz3*xl*xl+qz4*(x1)*+;.o+qz5*(xl)**4.O 

qz = qzz-qzzp 

qzzp = qzz 

10 continue 

+ +qz6* (xl) **5.0 

do 100 no = 1,norbit 

do 100 n = 1,ntime 

if (1 .eq. 1) thgas(n,no) = tgas 

do 30 jj = 1,nj 

do 30 ii = 1,15 

C 

C 

if (ii.gt.ne(1-j)) then 

t(ii,jj) = 0.0 

apon(ii,jj) = 0.0 

endif 

C 

C 

e 

0 

to(ii,jj) = t(ii,jj) 

romelt(ii,jj) = rmelt(ii,jj) 

apoo(ii,jj) = apon(ii,jj) 

30 continue 

if (n.gt.qtime) then 

else 

endif 

q(n) =‘qn 

q(n) = qd 

qaz = q(n) *qz/nj 

icount = o 

35 iflag = o 

e 
D3 

a 

e 

e 

a 

a 

e 

a 

C 

C 

C 

C 

C 

C 

C 

C 

24 

80 

40 

50 

icount = icount+l 

tgasin = thgas(n,no) 

ta=O. 0 

tt=O. 0 

do 24 j = 1,nj 

tt = tt + ne(j)-l 

le(j) = ro(j)-(j-l)*dy 

f = (1-exp(-cl*le(j))) 

ta = ta+f 

continue 

qt = qaz* (nj-ta) * (1-wf) 

qt = qt/tt 

qlat = 0.0 

qsens = 0.0 

qrad = 0.0 

tac = 0 

heat = 0.0 

hd = 0.0 

do 40 j = 1,nj 

call gasprop(dz,nj,h,l) 

heat = heat + qaz*dz 

hd = hd+h 

call conduct(borl,i,nj,icount,ft,qlat,qsens,qrad,tac) 

call heatgas(nj) 

call solver (i) 

continue 

hd = hd/nj 

do 50 jj = 1,nj 

do 50 ii = l,ne(]j) 

error = abs(t(ii,jj)-tp(ii,jj)) 

if (error .gt. errmax) iflag = 2 

tp(ii,jj) = t(ii,jj) 

rpmelt(ii,jj) = rmelt(ii,jj) 

continue 

qgas = gasmas*cpgas*deltdz 

qtotal = qlat+qsens+qgas 

qerror = abs(l.O-qtotal/heat) 

if (n .gt. qtime) then 

if (qerror .gt. 0.Ol)iflag = 2 
D4 

a 

e 

0 

else 

endif 

if (qerror .gt. 0.001)iflag = 2 

C 

. tgasp = thgas (n, no) +deltp 

C 

if (iflag .gt. 1) then 

if (icount .ge. imax) then 

write (7,1000) icount 

else 

endif 

goto 35 

endif 

thgas(n,no) = thgas(n,no)+deltdz 

do 47 j =l,nj 

if (no .It. norbit) goto 100 

write (7,1050) 

write (7,1250) 

write(7,1100)l,no,n 

write (7,1150) heat, qz 

write(7,1151) qtotal 

nr = 0.0 

do 83 j = 1,nj 

C 

C 

47 continue 

C 

1151 format(ix,'qtotal = ',f10.8) 

C 

0 nn( j ) =ne (j ) +j-1 

83 continue 

C 

do 34 j = 1,nj 

nr = amaxl(nr,nn(j)) 

34 continue 

C 

do 36 j = 1,nj 

do 37 i = 1,nr 

C 

if (i.gt.ne(j)) then 

rmelt(i,j) = 0.0 

theta(i,j) = 0.0 

theta(i,j) = t(i,j)+thmelt 

else 

endif 

C 

37 continue 

36 continue 

C 

do 85 j = 1,nj 
D5 

a 

do 85 i = 1,nr 

if (j.ne.1) then 

thn( (i+j-l) ,j) = theta(i, j) 

rmel( (i+j-l), j) = rmelt(i, j) 

else 

thn(i,j) = theta(i, j) 

rmel(i, j) = rmelt(i, j) 

C 

endif 

85 continue 

C 

C 

do 70 i = nr,l,-1 

write (7,1200) i, (thn( i, j ) , j=1, nj ) 

write(7,1300)tgasp,thgas(n,no) ,hd 

write (7,1050) 

write (7,1250) 

do 82 i = nr,l,-1 

write(7,1350) i, (rmel(i, j), j=l,nj) 

70 continue 

C 

C 

82 continue 

100 continue 

0 C 

................................................................ 

C * 

C input and output format statements * 

C * ................................................................ 

1000 format(lx,'failed to converge after',i6,'iterations1) 

1050 format (lx, '1 1250 forrnat(lx,'----------------------------------------------- ............................................... 

I 

1100 format(ix,'axial position = ',il,' orbit = ',il,' orbit 

1150 format(lx,'net solar heat input = ',lfl0.8,'qz = ',lf9.4) 

1200 format(lx, 'i = ',i2,20f7.1) 

1300 format(lx,'avg gas temp = ',f6.l,'k out gas temp I, 

1350 format(lx,'i = ',i2,20f7.4) 

---','--------------------I --- 

& position = ',i2) 

,f6.l,'kh = ',f7.4,'kw/rn2-k1) 

close (8) 

close (7) 

stop 

end 
DG 

C 

C 

a 

e 

C 

C 

C 

C 

subroutine gasprop(dz,nj,h,l) 

real k,kl,k2,kll,k12,k13,k141k211k221k231k241kgasl 

common/b1okb/ni,ne(10),j,t(15,10) 

common/blokc/thmelt,tgasin,gasmasmaslcpgas,ntu,deltpl 

common/blokd/tgasp,gasmix,ro(l5),ht,w,dxx,p 

data gamma,p~,run~v,mw1,mw2/1.66,3.14159,8.31434,4,131/ 

data k11,k12,k13,k14/0.1580e-3,0.9086e-6,333.333,0.794/ 

data k21,k22,k23,k24/0.6404e-5,0.0454e-6,333.333,0.800/ 

data viscll,viscl2,viscl3,viScl4/ 

data visc21,visc22,visc23,viSc24/ 

+ mw,mwl,mw2,ntu,nu 

+ deltdz 

/.2158e-4,.1022e-6,333.3,.825/ 

/.253e-4,.2452e-6,333.333,.77/ 

fl (x,y) = (1+2.41* (x-y) * (x-O.i42*y)/ (x+y) **2.0) 

f2 (xl, x2, x3 I x4 I x5, x6) = 

(2.0*2.0**. 5* (l+x5/x6) **O. 5) 

= ( (1+ (xl/x2) **. 5* (x3/x4) ** .25) **2.0) / 

+ 

kl = kll+k12*(tgasp-k13)**k14 

k2 = k21+k22*(tgasp-k23)**k24 

viscl = viscl~+viscl2*(tgasp-visci3)**visci4 

visc2 = visc2~+visc22*(tgasp-visc23)**~isc~~ 

viscwl = viscll+viscl2*(t(l,j)+thmelt-viscl3)**viscl4 

viscw2 = visc2l+visc22*(t(l,j)+thmelt-visc23)**visc24 

phi12 = f2(viscl,visc2,mw2,mwl,mwl,mw2) 

phi21 = phil2*(visc2*mwl)/(viscl*mw2) 

psi12 = f2(kl,k2,mwl,mw2,mwl,mw2)*fl(mwl,mw2) 

psi21 = f 2 (k2 I kl I mw2 I mwl I mw2 I mwl) *f 1 (mw2 ,mwl) 

visc = viscl/ (l+phil2* (1 - gasmix) /gasmix) + 

+ visc2/(l+phi21*(gasmix/(l - gasmix))) 

viscw = viscwl/(l+phil2*(l - gasmix)/sasmix)+ 

+ viscw2/(l+phi2l*(gasmix/(l - gasmix))) 

kgas = kl/(l+psil2*(1 - gasmix)/gasmix)+ 

+ k2/(1 +psi21*(gasmix/(l - gasmix))) 

mw = gasmix*mwl+(l-gasmix)*mw2 

cpgas = gamma*runiv/((gamma-l)*mw) 

pr = cpgas*visc/ (kgas) 

dt = 2.0*ht*w/(dxx*n-j) 

re = 2.0*gasmas/((p*.5*dt)*visc) 

nut = .027*re**.8*pr**.333333*(visc/viscw)**.l4 

if ((l.eq.l).and.(re.lt.2lOO.O))then 

nu = .44* (dz/ (re*pr*dt) ) ** (-. 66) 
D7 

0 

0 

elseif (re .It. 2100.0) then 

elseif ((re.gt.2100.0).and.(re.lt.lOOuO.O)) then 

else 

endif 

h = kgas*nu/dt 

ntu = h*dxx*dz*nj/(gasmas*cpgas) 

if (l.eq.1) ntu = .2337841 

if (l.eq.2) ntu = .5330746 

if (l.eq.3) ntu = .4739715 

if (l.eq.4) ntu = .4153427 

if (l.eq.5) ntu = .41995574 

return 

end 

nu = 2.35 

nu = amaxl(2.35,nut) 

nu = nut 

C 

C 

a 
D8 

0 

e 

0 

a 

0 

0 

subroutine conduct 

+ (borl, i, nj , icount, ft I qlat I qsens, qrad, tac) 

c 

common/bloka/a (15) , b (15) , c (15) , d( 15) 

common/blokb/ni,~ne(lO),~,t(l5,lO) 

common/blokc/thmelt,tgasin,gasmas,cpgas,nt~,del~p, 

common/blokd/tgasp,gasmix,ro(l5),ht,w,dxx,p 

common/bloke/rmelt(l5,lO),denliq,densol,walls,cvmet, 

common/blokf/dz,dtime,rpo(lO)Iqaz,hfl,romelt(l5,lO),c~,~f, 

common/blokg/k,kn,ks,ke,kw,den,cv 

real k,kn,ks,ke,kw,kmet,ntu,ksol,kliq,le 

+ deltdz 

+ denmet,kmet,kliq,ksol,cvliq,cvsol 

+ t0(15,l0),rpmelt(l5,10),dx,dy,th,qt,le(l5), 

+ apoo (15,lO) , apon (15,lO) 

C 

ne(j) = 200 

i=1 

as = 0 

C 

1 if (j .eq. 1) then 

aw = 0 

tw = 0 

tw = t(i,j-1) 

elseif (i. eq. 1) then 

else 

tw = t(i+l,j-l) 

endif 

C 

if (j .eq. nj) then 

ae = 0.0 

te = 0.0 

elseif (i. eq. 1) then 

te = t(i,j+i) 

else 

te = t(i-i,j+i) 

endif 

if (i .ge. 2) then 

C 

r = (float(i)+j-2.5)*dy 

rn = r + .5*dy 

rs = r - .5*dy 

endif 

C 

................................................................ * 

C set nodal constants for the wall ................................................................ 

C 

if (i .eq. 1) then 

an = kmet*dxx/th 
D9 

0 

I* 

a 

a 

0 

a 

0 

0 

a 

if (j .ne. 1) aw = kmet*th/dxx 

if (j .ne. nj) ae = kmet*th/dxx 

apon(i,j) = denmet*cvmet*th*dxx/dtime 

apo = apon(i,j) 

dapo = 0.0 

dtapo = 0.0 

tcij = tgasin-thmelt 

hij = gasmas*cpgas* (1-exp (-ntu) ) / (dz*nj ) 

abw = wf * (exp (-cl*le (j ) ) ) 

qr = qaz*abw 

qij = qr 

goto 40 

endif 

c; 

if (rpo(j) .gt. rn) then 

endif 

if ((j .eq. nj).and.(i .eq, 2)) then 

endif 

if (i .gt. (ni-j+2)) then 

if (rpo(j) .le. (rn+dy)) ne(]) = i+l 

C 

if (rpo(j).le.rn) ne(j) = i 

C 

rmelt(i,j) = 1.0 

c; 

if (i .eq. ne(j)) then 

elseif (i .eq. (ne(])-1)) then 

else 

endif 

an = 0.0 

an = kliq*dx/ ( .5* (rpo ( j ) -rn) + .5*dy) 

an = kliq*dx/dy 

C 

if (i .eq. (ni-j+3)) then 

call property (cv,den,i-l,j,ks) 

if (i.eq.ne(j)) then 

else 

endif 

as = kliq*dx/( .5*(rpo(j)-rs+dy)) 

as = kliq*dx/dy 

as = kliq*ks*dx/(kliq*.5*dy+l:s*(rpo(j) -rs) *.5) 

as = kliq*ks*dx/((kliq+ks)*.5*dy) 

elseif (i .eq. ne(])) then 

else 

endif 

C 
D10 

0 

I. 

0 

e 

0 

0 

if (j .ne. nj) then 

if (rpo(j+l) .ge.rn) then 

if (i.eq. ne(])) then 

else 

ae = kliq* ( (rpo ( j ) +rn) * .5-rs) /dx 

ae = (kliq*dy)/dx 

endif 

elseif(rpo(j+l).le.rs) then 

elseif (i.eq.ne(j)) then 

else 

endif 

ae = 0.0 

ae = (kliq*( (rpo(j)+rpo(j+l) )*.5-rs) )/dx 

ae = kliq*((rn+rpo(j+l))*.5-rs)/dx 

endif 

if (j .ne.l) then 

C 

if (rpo (j -1) . ge. rn) then 

i’f (i.eq. ne(])) then 

else 

aw = kliq* ( (rpo (j ) +rn) *. 5-rs) /dx 

aw = (kliq*dy)/dx 

endif 

elseif(rpo(j-1) .le.rs) then 

elseif (i.eq.ne(j)) then 

else 

aw = 0.0 

aw = (kliq* ( (rpo (j ) +rpo( j-1) ) *. 5-rs) ) /dx 

aw = kliq* ( (rn+rpo ( j -1) ) *. 5-rs) /dx 

endif 

endif 

if (i .eq. ne(j)) then 

C 

apon(i,j) = denliq*cvliq*(l.O - walls)*(rpo(j)-rs)* 

apon(i,j) = denliq*cvliq*(l.O - walls)*dy*dx/dtime 

* dx/dt ime 

else 

endif 

if (apoo(i,j).eq.O.O) apoo(i,j) = apon(i,j) 

C 

if (qaz .gt. 0.0) then 

apo = apoo(i,j) 

dapo = apon(i,j) - apoo(i,j) 

dtapo = 0.0 

apo = apon(i,j) 

dapo = 0.0 

else 

dtapo = apoo(i,]) - apon(i,j) 

endif 

tcij = 0.0 
D11 

0 

0 

0 

a 

0 

a 

a 

hij = 0.0 

if (i .eq. ne(])) then 

else 

abm = 1 - exp(-cl*(rpo(j)-rs)) 

abm = (-exp (-cl* (rpo (j ) -rs) ) +exp (-cl* (rpo ( j ) -rn) ) ) 

endif 

qr = qaz*abm + qt 

qij = qr 

goto 40 

endif 

C 

................................................................ 

................................................................ 

C set nodal constants for the inner pcm * 

23 call property (cv,den,i,j,k) 

C 

if (i .eq. 2) then 

else 

apon(i,j) = den*cv*.5*dx*dy/dtime 

pcmmass = .5*dx*dy*densol*(l-walls~ 

apon(i,j) = den*cv*dx*dy/dtime 

pcmmass = dx*dy*densol*(l-walls) 

endif 

if (apoo(i,j) .eq. 0.0) apoo(i,j) = apon(i,j) 

if (qaz .gt. 0.0) then 

apo = apon(i,j) 

dapo = 0.0 

C 

dtapo = apoo(i,j) - apon(i,j) 

apo = apoo(i,j) 

dapo = apon(i,j) - apoo(i,j) 

dtapo = 0.0 

else 

endif 

tcij = 0.0 

hij = 0.0 

dmelt = t(i,j)*cv*borl/hfl 

if (i.eq. ne(])) then 

elseif (i.eq. (ni-j+2) ) then 

C 

C 

an = 0.0 

if ( (i+l) . eq. ne (j ) ) then 

else 

endif 

call property(cv,den,i+l,j,kn) 

an = k*kn*dx/((k+kn)*.5*dy) 

an = k*kliq*dx/(k*(rpo(j)-rn)*.5+kliq*.5*dy) 

an = k*kliq*dx/((k+kliq)*.5*dy) 

else 

0 
D12 

a 

0 

0 

a 

0 

endif 

if (i .eq. 2) then 

else 

C 

as = kmet*dxx/th 

call property(cv,den,i-l,j,ks) 

as = k*ks*dx/ ( (k+ks) *. 5*dy) 

endif 

if (j .ne.nj) then 

C 

if (i. eq. 2) then 

else 

ae = 0.0 

call property(cv,den,i-l,j+l,ke) 

ae = k*ke*dy/ ( (k+ke) *. 5*dx) 

endif 

endif 

-if (j.ne.1) then 

endif 

C 

call property(cv,den,i+l,j-1,kw) 

aw = k*kw*dy/ ( (k+kw) *. 5*dx) 

C 

................................................................ 

if (i.eq.2) then 

if (i.eq. ne(])) then 

if (dmel t . gt . 0 . 0) then 

if ((t(i-l,j).ge.O.O).or.(rmelt(i+i,j).ge.i.O) 

& .or.(qaz.gt.O.O)) 

& rmelt(i, j) = rmelt(i,j)+dmelt 

else 

if ((t(i-l,j).le.o.o).or.(qaz.lt.O.O)) 

& rmelt(i,j) = rmelt(i,j)+dmelt 

endif 

if (dmelt.gt.O.0) then 

else 

if (rmelt(i+l,j) .ge. 1.0) 

& ' rmelt(i,j) = rmelt(i,j)+dmelt 

& rmelt (i, j) = rmelt (i, j) +dmelt 

else 

if ((t(i-1,j) .le.O.O) .or. (qaz.lt.O.0)) 

endif 

endif 

if (dmelt .gt. 0.0) then 

else 

elseif (i.eq.ne(j)) then 

if (qaz .gt. 0.0) rmelt(i,j) = rmelt(i,j)+dmelt 

if ((rmelt(i-1,j) .le. O.O).or.(qaz .It. 0.0)) 

& rmelt(i, j) = rmelt(i,.j)+dmelt 

endif 
D13 

0 

e 

a 

d 

a 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

else 

if (dmelt.gt.O.0) then 

if (rmelt(i+l,j) .ge. 1.0) 

E4 rmelt(i,j) = rmelt(i,j)+dmelt 

if ((rmelt(i-1,j) .le. 0.0) .or. (qaz .It. 0.0)) 

& rmelt (i, j) = rmelt (1, j) +dmelt 

else 

endif 

endif 

if (rmelt(i, j) .gt.l.o) rmelt(i, j) = 1.0 

if (rmelt(i,j) .It. 0.0) rmelt(i,j) = 0.0 

qij = (romelt(i,j)-rmelt(i,j)) *pcmmass*hfl/dtime 

qlat = qlat-qij *dz 

abm = (-exp (-cl* (rpo (j ) -rs) ) +exp (-cl* (rpo (j ) -rn) ) ) 

qr = qaz*abm + qt 

qij = qij + qr 

if (i .le. (ni-j+2)) then 

if (i.eq.2) then 

odmelt = (rmelt(i,j)-rpmelt(i,j))*.5*dy* 

odmel t = ( rmel t ( i , j ) -rpmelt ( i , j ) ) *dy* 

* (densol/denliq-1) 

else 

* (densol/denliq-1) 

endif 

ro(j) = ro(j) + odmelt 

endif 

if (rmelt( (ni-j+2) ,j)-0.0.1e.0.000001) then 

endif 

ro(j) = ni*dy 

ne(j) = ni-j+2 

odmelt = 0.0 

40 a(i) = -as 

b(i) = an+as+ae+aw+apo+dapo+hij 

c(i) = -an 

d (i) = (apo+dtapo) *to (i, j ) +qij+ae*te+aw*tw+hij *tcij 

delapo = apoo(i,j) - apon(i,j) 

if (delapo.ge.O.0) then 

else 

qmove = delapo*(O.O-to(i,j))*dz 

qmove = delapo*(O.O-t(i,j))*dz 

endif 

qsens = qsens+(apo+dapo)*dz*t(i,j)-(a~o+dtapo)*dz*to(i,j) 

qrad = qrad + qr*dz 

if (i.ne.1) tac = tac+abm 

if (i .eq. ne(j)) goto 60 

50 i = i+l 

a 
a 

a 

D14 

got0 1 

60 return 

end 
D15 

a 

a 

subroutine property(cv,den,i,j,k) 

common/bloke/rmelt(l5,lO),denl~q,densol,walls,cvmet,denmet, 

real kliq,ksol,kmet,k,kpure 

kpure = 1.0 / (rmelt(i,j)/kliq + (1-rmelt(i,j))/ksol) 

cvpure = rmelt(i,j)*cvliq + (1-rmelt(i,j))*cvsol 

denpur = rmelt(i,j)*denliq + (1-rmelt(i,j))*densol 

C 

+ kmet,kliq,ksol,cvliq,cvsol 

C 

c******************************************~k******************* 

C put storage material in parallel with can material * 

C*************************************************************** 

k = kpure*(l-walls)+kmet*walls 

den = denpur*(l-walls)+denmet*walls 

cv = (cvpure*denpur*(l-walls)+walls*cvmet*denmet)/den 

return 

end 

a 
C 

a 

C 

C a 

C 

0 

a 

D1G 

subroutine solver(i) 

common/bloka/a (15) , b (15) , c (15) , d (15) 

common/blokb/ni, ne (10) , j , t (15, IO) 

m = ne(j) 

do 1 i = 1,m-1 

r = a(i+l)/b(i) 

b(i+l) = b(i+l)-r*c(i) 

d(i+l) = d(i+l)-r*d(i) 

1 continue 

do 2 i = m-l,l,-1 

t (i, j) = (d( i) -c( i) *t (i+l, j) )/b( i) 

2 continue 

return 

end 
D17 

0 

0 

e 

subroutine heatgas(n1) 

real ntu 

common/blokb/ni,ne(lO),~,t(15,10) 

common/blokc/thmelt,tgasin,gasmas,cpgas,ntu,deltp, 

+ deltdz 

if (j .eq. 1) then 

qgasp = 0.0 

qgasdz = 0.0 

C 

endif 

temp=(t(l,j)+thmelt) -tgasin 

qgasp=qgasp+gasmas*cpgas*temp*(i-exp(-~.5*ntu))/(nj) 

qgasdz=qgasdz+gasmas*cpgas*temp*(l-exp(-ntu))/(nj) 

if (j .eq. nj) then 

endif 

10 continue 

return 

end 

C 

C 

deltp=qgasp/(gasmas*cpgas) 

deltdz=qgasdz/(gasmas*cpgas) 

C 

a 

0 
NASA 

1. Aewrt No. 

Report Documentation Page 

2. Government Accession No. 3. Recipient's Catalog No. . 

NASA CR-185117 

4. Title and Subtitle 5. Repon Date 

Thermal Evaluation of Advanced Solar Dynamic 

Heat Receiver Performance 

7. Author@) 

Roger A. Crane 

March 1989 

8. Performing Organization Reporl No. 

None 

10. Work Unit No. 

506-41-31 

9. Performing Organizathn Name and Address 

University of South Florida 

Tampa, Florida 33620 

2. Sponsoring Agency Name and Address 

Project Manager, Miles 0. Dustin, Power Technology Division, 

NASA Lewis Research Center. 

~ 

11. Contract or Grant No. 

NAG3-85 1 

13. Type of Report and Period Covered 

Contractor Report 

6. Abstract 

Lewis Research Center 

Cleveland, Ohio 44135-3191 

The work described herein addresses the thermal performance of a variety 

of concepts for thermal energy storage as applied to solar dynamic applications. 

It is recognized that designs providing large thermal gradients or large temperature swings during orbit are susceptible to early mechanical failure. 

incorporating heat pipe technology may encounter operational limitations over 

sufficiently large ranges. 

the relative merits of the basic concepts are to be compared. 

effect of thermal enhancement and metal utilization as applied to each design provides a partial characterization of the performance improvements to be achieved 

by developing these technologies. 

Concepts 

By reviewing the thermal performance of basic designs 

In addition the 

14. Sponsoring Agency Code 

7. Key Words (Suggested by Author(s)) 18. Distribution Statement 

9. Security Classif. (of this report) 

Unclassified 

Solar Dynamic 

Heat Receiver 

20. Security Classil. (01 this page) 21. No of pages 22. Price' 

Unclassified 140 A0 7 

, . .. Unclassified-Unlimited I 

a 

4 

4 

(I 

'For sale by the National Technical Information Service, Springfield, Virginia 22161 NASA FORM lOZ@ OCT M