Ava.

tr 5.1.:

I
c

in.“ .fiufinwwzwu
: .
«Wk

 

.I z N
. #3..) a

r
r. «3% N
V dihumflofl...muum§ :uufib .
1 civfiflrl. FM?
2 .1 .. H v
.. Sukiyaki
v.0 (2?... 44
.i.

F.

.....$

.9 .r ..
r.‘ a? .

bu. pa.) 3.

 

 

77-33715
2 LIBRARY
2003 Michigan State

,I, ,1 ,,_ W Universi
691 W w W

This is to certify that the
dissertation entitled

MATHEMATICAL MODELS TO PREDICT THE PERFORMANCE
OF INSULATING PACKAGES AND THEIR PRACTICAL USES

presented by
Seung-Jin Choi

has been accepted towards fulfillment
of the requirements for the

Ph.D. degree in School of Packaging

find '13 ummutv/

' Majoi‘ Professfir’s Signature

i2//§/0’{

Date

 

 

MSU is an Affirmative Action/Equal Opportunity Institution

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

w (fig! y'alggoa

JAN 2 93111; 1 a
on ‘ 1

I.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6/01 c;/CIRC/DateDue.p65-p.15

MATHEMATICAL MODELS TO PREDICT THE PERFORMANCE OF
INSULATING PACKAGES AND THEIR PRACTICAL USES

By

Seung-Jin Choi

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY
School of Packaging

2004

ABSTRACT

MATHEMATICAL MODELS TO PREDICT THE PERFORMANCE OF
INSULATING PACKAGES AND THEIR PRACTICAL USES

By

Seung-Jin Choi

Thermal insulation is used in variety of applications to protect temperature
sensitive products from thermal damage. Several factors affect the performance
of insulation packages. These factors include the heat transfer through the
packaging material (conduction, convection and radiation), the geometry of the
insulating package and the contact resistance between the product and the
package. In this study, a comprehensive model, which includes all of these
factors, was developed to predict the performance of the insulating package.

First, a computer model was developed to represent the heat transfer
through multi-Iayered walls and an equation was derived for the calculation of the
thermal resistance of the multi-Iayered wall.

Thereafter, three mathematical models, which included the geometry of
the insulating package, contact between the insulating package and product and
contact between the package and the ground, were developed to investigate the
heat transfer through the entire insulating package. Thermal resistances of entire
packages calculated by these models were compared to ice-melt test data for
evaluation. The model, which included the geometry of the insulating package
and contact between the insulating package and ice container, successfully

predicted experimental data within 20% error.

Based on this model, a simplified equation was derived for the practical
prediction of the thermal resistance of the entire insulating package. This
equation showed dramatically improved results compared to the conventional
equafions.

Examples of the utilization of this model are also included.

ACKNOWLEDGEMENTS

I still remember the clear blue sky and unseasonably cool, fresh day of
July, 2000 when I first stepped out of the plane into the Lansing airport. I also
remember my first class PKG323 at 8:00am in the Brody Hall complex. It was a
cloudy day in late August. Since then, there have been a lot of people who have
supported my efforts in earning my doctoral degree.

First of all, I would like to express my gratitude to Dr. Burgess, my advisor
and committee chair, for his advice and support throughout my doctoral program.
With his great insight as an engineer, he always guided me to an answer
whenever I was in need of one. With his great personality, he always listened to
me and guided me in a respectful and encouraging way. He was both a mentor
and colleague to me. I would also like to thank my committee members, Dr.

Singh, Dr. Clarke and Dr. Feeny, who challenged me to be a scholar and a
researcher.

I owe my thanks to many friends who provided me with every kind of
support, from encouraging remarks to wonderful meals. I would especially like to
thank Dr. Simon B. Shin, a true scholar and my father-in-law. Whenever I visit his
study, I’m always challenged by full of books and drafts he is working on. I would
also like to thank my parents and my mother-in-law for their thoughtful support
throughout earning my degree. Finally, I would like to thank my wife, Dr. Nary
Shin, who shared her experience and insight as a colleague in academia, and to
my daughter, Alexis, for the joy that she has brought and continues to bring to my
life. To everyone who has made it possible to complete my dissertation, I extend
my sincere gratitude and love.

TABLE OF CONTENTS

LIST OF TABLES ................................................................................................ vii

LIST OF FIGURES ............................................................................................. viii
CHAPTER 1

INTRODUCTION ................................................................................................... 1
CHAPTER 2

LITERATURE REVIEW ......................................................................................... 7

2.1 Insulating packages ...................................................................................... 7

2.2 Distribution environments and their simulations ......................................... 11

2.3 Test methods for thermal insulation quality of packages ............................ 13

2.4 Simulation of the thermal insulation quality of packages ............................ 14
CHAPTER 3

THEORETICAL DEVELOPMENT ....................................................................... 17

3.1 British thermal unit (BTU) ........................................................................... 17

3.2 Basic theories of heat transfer in insulating package ................................. 17

3.2.1 Conduction heat transfer .................................................................... 17

3.2.1.1 Heat transfer through a solid wall .............................................. 18

3.2.1.2 Calculation of effective area for conduction through the walls of

the insulating package ................................................................ 20

3.2.2 Convection heat transfer .................................................................... 22

3.2.2.1 Natural convection in unconfined spaces .................................. 25

3.2.2.2 Natural convection in enclosed spaces ..................................... 27

3.2.2.3 Heat transfer balance by convection in enclosed spaces .......... 30

3.2.3 Radiation heat transfer ....................................................................... 32

3.2.4 Combined radiation and convection ................................................... 33

3.3 Mathematical models for heat transfer through multi-layered walls ............ 35

3.4 Mathematical models for heat transfer through entire insulating package..40
3.4.1 Mathematical model considering the geometry of insulating package41
3.4.2 Mathematical model considering the contact between the Ice container

and the inside of the insulating package ............................................. 45

3.4.3 Mathematical model considering all contact areas ............................. 52
CHAPTER 4

EXPERIMENTAL DESIGN .................................................................................. 57

4.1 Materials and methods ............................................................................... 57

CHAPTER 5
RESULTS AND DISCUSSION ............................................................................ 62
5.1 Characteristics of the heat transfer in multi-layered walls ........................... 62
5.2 Simplification of the calculation of the system R- value of multi- layered
walls ............................................................................................................ 63
5.3 Evaluation of the simplified equation using computer program .................. 69
5.4 Simulation of the ice-melt test using mathematical models ........................ 76
5.4.1 Construction of mathematical model considering the geometry of
insulating package .............................................................................. 76
5.4.2 Construction of mathematical model considering the contact between
the ice container and the insulating package ...................................... 77
5.4.3 Construction of mathematical model of the insulating package
considering all contacts ....................................................................... 77
5.4.4 Comparison of models with experimental data .................................. 78
5.5 Construction of simplified model for practical use ...................................... 81
CHAPTER 6
CONCLUSIONS .................................................................................................. 84
6.1 Practical use of model equations ................................................................ 84
6.2 Example ..................................................................................................... 87
6.3 Factors affecting on package construction ................................................. 92
6.3.1 The effect of wall thickness on the performance of the insulating
package .............................................................................................. 92
6.3.2 The effect of aluminum foil on the performance of the insulating
package .............................................................................................. 94
6.3.3 The effect of the surface area on the performance of the insulating
package .............................................................................................. 96
APPENDICES ..................................................................................................... 98
APPENDIX A. Basic programs for simulation of the insulating package .......... 99
APPENDIX B. Thickness of solid layers of the multi-layer wall ...................... 133
APPENDIX C. The thickness of the air space between solid layers ............... 136
APPENDIX D. The emissivities of the 'between-the-layer' air space .............. 139
APPENDIX E. R-values and parameters generated by computer program ....143
REFERENCES .................................................................................................. 146

vi

LIST OF TABLES

Table 1. Various distribution conditions and their temperatures ............................ 2
Table 2. Thermal conductivities of various insulation materials ............................ 4
Table 3. The extremes of the thermal cycles during transportation ..................... 12

Table 4. Comparison of geometric average and arithmetic average with

effective area in equation (12) ............................................................... 23
Table 5. Constants for calculating the Nusselt number ....................................... 27
Table 6. Summary specifications for insulation materials .................................... 58
Table 7. Predicted R-values of three mathematical models ................................ 61

Table 8. Comparison of R-values in equation (87) with the computer model ...... 71
Table 9. Comparison of R-values in equation (88) with the computer model ...... 73
Table 10. Thermal resistance per inch for various insulating materials ............... 74
Table 11. Summary of the experimental data and areas according to

Figure 13. .............................................................................................. 79
Table 12. Comparison of simulation results of three computer models with the

experimental data .................................................................................. 80
Table 13. Comparison of the simulation results with equation (93) with

experimental data. ............................................ 83

Table 14. System R-values with various wall thickness and foil locations ........... 95

vii

LIST OF FIGURES

Figure 1. Expanded polystyrene (EPS) foam container ........................................ 7
Figure 2. Molded polyurethane container .............................................................. 8
Figure 3. Vacuum insulation panels (VIP) container ............................................. 9
Figure 4. Gas-filled panel (GFP) container: (A) GFP in a corrugated box, .......... 10
Figure 5. Diagram for conduction heat transfer through a wall ............................ 18
Figure 6. Relation between temperature and the physical properties of air ........ 28

Figure 7. Diagram of the enclosed air space between the surfaces a and b ....... 31
Figure 8. Diagram of heat transfer through multi-layered walls in contact with
inside and outside air .......................................................................... 36
Figure 9. Iteration procedure to solve non-linear problems. ................................ 38
Figure 10. Diagram of the model considering the geometry of the insulating
package. ............................................................................................. 42
Figure 11. Diagram of thermal contact resistance between surface a and b: ...... 46

Figure 12. Diagram of the model considering the contact between the Ice

container and the inside of the insulating package ............................. 48
Figure 13. Diagram of the model considering all contact areas .......................... 53
Figure 14. The procedure for the regular ice melt test ........................................ 60

Figure 15. The Grashof number with the distance between the layers and

the temperature difference between the layers ................................... 64
Figure 16. Resistance of the air space with various emissivities ......................... 67
Figure 17. Diagram of typical insulating package ................................................ 84
Figure 18. Temperature profile after the product melts ....................................... 91

viii

CHAPTER 1

INTRODUCTION

Temperature sensitive products such as biological materials,
pharmaceuticals, industrial adhesives, gyroscopes, blood, frozen foods, fresh-
products and dairy products should be shipped in temperature controlled
containers which keep the temperature constant during transportation and keep
products from thermal damage caused by melting, thawing or freezing.

In general, a number of steps are involved in the distribution process of
the product. According to ASTM 04169-98, there are 3 or 5 handling operations
in the overall distribution procedure (ASTM 1999a). Even the simplest case of
distribution process has handling operations: 1) packages are loaded on a truck
at the manufacturing site, 2) packages are taken from the truck into the
warehouse of the distributor, 3) packages are loaded on a vehicle from the
warehouse to the customer. In this case, the temperature may not be controlled
and fluctuations in the temperature could include very extreme cases. Various
distribution conditions and their temperatures are summarized in Table 1.

In most cases, thermal insulation is used to protect the products from
thermal damage during distribution of products. Thermal insulation is a procedure,
which can retard heat transfer and eventually conserve energy by reducing heat
loss or gain of the product (ASHRAE, 1993). Due to the complexity of the nature

of heat transfer, several factors can affect the insulating ability of the containers.

Table 1. Various distribution conditions and their temperatures (Panyarjun, 2002)

 

 

Conditions Temperatures
Freezer -18 to —35 °C
Refrigerator 1 to 4 °C

Freezer truck, rail car and airplane

Refrigerated truck, rail car and
airplane

Outside environment
(Depending on locations and
seasons)

—18 °C or below
(~29 °C or below for ice cream)

1to4°C

Based on ASTM 044332-89

20: 2 °C for Temperature high humidity
40 i 2 °C for Tropical condition
20 i 2 °C for Desert condition

Extremes, based on FedEx assuming:

-50 °C for carrier vehicles and open dock
areas during the winter in northern
climates

60 °C for closed, parked carrier vehicles
during summer in southern climates

Extremes, based on US, military
assuming:

-50 °C lowest for worldwide
70 °C highest for worldwide

 

Conduction, convection and radiation are involved in the heat transfer
through the insulating package. Heat transfer phenomena were studied
intensively during last century (Holman, 1986; Kreith, 1973; McCabe et al., 1993).
A material’s resistance to heat penetration is a very crucial factor. There are
several reports on thermal conductivities of the insulating materials. A short table
of thermal conductivities for common insulating materials used in packaging is
shown in Table 2. The geometry of the package can also affect the insulating
capacity of the container. The size, shape, thickness and structure of the wall of
the insulation package affect its insulating ability (Burgess, 2004). Generally,
larger containers have more surface area and absorb heat faster than smaller
containers. On the other hand, the container with larger interior volume can hold
more products and ice, which can hold the temperature longer. As the container
gets larger, the volume increase is faster than the surface area increases. For
this reason, large containers have more potential to hold the temperature longer
than small containers. The shape of the container is also an important factor
when designing the insulating package. For a given volume, one should design
the package with the smallest surface area and maximum distance from its
surface to center of the product. The structure and the thickness of the wall affect

the insulating capacity.

Table 2. Thermal conductivities of various insulation materials
(Kositruangchai, 2003)

 

Thermal Conductivity at 70°F

 

Matenals (BTU/hrftz'W)
Air 0.015
Water 0.346
Corrugated board 0.035
EPS foam 0.022
Polyurethane 0.018

 

It is obvious that thicker walls generate better insulating capability. One
should consider that multi-layered materials make better insulators than one thick
material with the same thickness. A thin blanket of air trapped between the liner
and box produces a significant resistance to heat penetration. This is the same
principle behind double and triple pane windows or wearing layered clothes in the
winter instead of one thick coat. Use of aluminum foil can also improve the
capacity of the insulating package. As mentioned above, one of the major heat
transfer mechanisms of the insulating package is infrared radiation. Most
insulating materials absorb about 95% of the infrared radiation. On the other
hand, aluminum absorbs only 5% of the radiation. For this reason, the addition of
aluminum foil can dramatically reduce the infrared radiation, which results
substantial improvement of the insulating ability of the package. Sealing the
package also plays important role for insulating ability of the package. Air
currents can flow in and out through very small openings and can carry enough

heat with them to render even the best insulator ineffective. For example,

corrugated boxes lined with EPS sheets on all 6 faces cannot insulate as well as
molded EPS container because there are always gaps along the edges where
the sheets meet.

Because of the variety of factors affecting the performance of a thermally
insulated package, a comprehensive model, which can represent these factors
as much as possible, is needed for designing efficient insulating packages.
Proper simulation model can reduce time-consuming efforts of preliminary
specifications and subsequent validation tests. It also eliminates over-packaging
and the resultant unnecessary costs. However, there were very few attempts to
predict the ability of insulation packages (Burgess, 1999; Kositruangchai, 2003;
Stavish, 1984). Moreover, these models were over simplified or missed very
crucial factors of heat transfer through insulating packages.

The objective of this study was to develop an effective model which could
successfully predict the complexity of heat transfer phenomena of insulating
packages. Simulation using mathematical models is very important for designing
new insulating packages. Using proper mathematical models, one can predict the
performance of the insulating package without constructing actual insulating
packages. This procedure will save time, efforts and cost by minimizing the
construction of actual packages and subsequent performance tests.

To accomplish this objective, three types of heat transfer phenomena
(conduction, convection, and radiation) were analyzed. Using the results of the
analysis, a computer model for one-dimensional heat transfer through a multi-

layered wall structure was developed. Thereafter, more complicated computer

models, which considered two-dimensional geometries and contact resistance,
were developed and compared with experimental data. Eventually, a simplified
equation was derived for a practical prediction of the insulating capacity of

packages.

CHAPTER 2

LITERATURE REVIEW

2.1 Insulating packages

The purpose of using insulation materials (those with low thermal
conductivities or high thermal resistances) is to conserve energy by reducing
heat loss or gain, to control surface temperature, to facilitate temperature control
of a chemical process in products, to prevent vapor condensation at the surface,
or to reduce temperature variations within a products (ASHRAE, 1993). There
are several insulating packages that are used for temperature sensitive products.

Expanded polystyrene foam (EPS foam), molded polyurethane, vacuum
insulation panels (VIP), gas-filled panels (GFP), and corrugated with liners or
blankets are the most widely used insulating materials. EPS foam is the most

common packaging foam used for low temperature distribution (See Figure 1).

 

Figure 1. Expanded polystyrene (EPS) foam container

The air entrapped in EPS foam gives excellent insulating capability since
air has very low heat conductivity. EPS foam is a relatively inert material and
acceptable for use in food packaging for meat and produce. EPS foam is
lightweight, inexpensive, reusable and stackable (Burgess, 2004; Hernandez,
2000). The major drawbacks of EPS foam are its difficulty of disposal and
difficulty in surface printing (Sasaki and Kato, 1999).

A molded polyurethane container is designed with a layer of high
performance urethane foam injected between layers of corrugated fiberboard
(See Figure 2). In general, urethane has higher insulating capability comparing to
EPS foam. For this reason, urethane requires less refrigerant than EPS foam. In
addition, polyurethane, which is injected between two corrugated board boxes,
increases the strength of the box and stability (Panyarjun, 2002). However, the

cost of urethane foam is higher than that of EPS foam.

 

Figure 2. Molded polyurethane container

Vacuum insulation panels (VIP) are advanced thermal insulation materials
that significantly outperform insulation materials such as closed-cell foams, foam
beads or fiber blankets (See Figure 3). It is made by placing open-cell
polyurethane planks in an insulating chamber and removing all the gases within
the insulating chamber (Burgess, 2004). The plank is wrapped in aluminum foil to
maintain the vacuum and the foil is protected by polyethylene film or any other
plastic polymers. In general, VlP’s are placed within another container such as a
corrugated box or protective covering during transportation. DuPont-Tejin Films,

Dow Chemical and Glacier Bay are among the major manufacturers of VIP.

 

Figure 3. Vacuum insulation panels (VIP) container

Gas-filled panels (GFP) are also used for insulating packages. It is
composed of exterior films made from HDPE and interior films with low-emissivity
surfaces (Griffith et al., 1991). Metalized films are used for the interior films. In

this case, low-conductivity gas-filled cavities and a series of low-emissivity

honeycomb type interior film layers minimize radiation, convection and
conduction. A typical thermal resistance of these materials is about 5.2
hrftz'oF/BTUin (Griffith et al., 1991). The structure of the GFP is shown in Figure 4.
Cargotech Airliner® by Cargo Technology is commercially used for thermal

insulation of single parcel distribution of temperature sensitive products.

 

Inside Air

=I Exterior Film: HDPE

— Interior Film: Metalized Film

Figure 4. Gas-filled panel (GFP) container: (A) GFP in a corrugated box,
(B) Diagram of the structure of the GFP

1O

Several modified corrugated boards were also developed to improve
thermal insulation ability of the corrugated board. For example, modified
corrugated board lined EPS foam was developed (Sasaki and Kato, 1999). This
material consists of a core layer, which is an EPS foamed sheet, covered on both
sides with a paperboard liner. EPS foam boards are used to improve heat
insulation by stopping the movement of air and decrease the absorption of water.
This material is collapsible, light and can be recycled separately from paperboard.
The major drawback of these materials is that the insulating ability is not much

improved comparing to that of corrugated board (Sasaki and Kato, 1999).

2.2 Distribution environments and their simulations

The environment of distribution, especially the temperature, affects the
performance of the insulating package. In the business of temperature sensitive
products, products pass through a variety of distribution processes as they are
delivered from the manufacturer to the customer. During these processes, the
product is exposed to several thermal cycles. ISTA 3G was designed to simulate
these temperature conditions (ISTA, 1999). According to this procedure, the
cycle of atmospheric conditions for 24 hour, 48 hour and 72 hour delivery are
summarized in Table 3 (Panyarjun, 2002). Using these conditions, individual

packaged products are tested for protection against temperature extremes.

11

Table 3. The extremes of the thermal cycles during transportation

 

24 hour domestic express small package freight transportation

 

 

 

 

 

 

 

Winter Profile Winter Profile
(Cold shipping and cold receiving) (Hot shipping and hot receiving)
Temperature Hours Temperature Hours
18 °C (65 °F) 0-4 22 °C (72 °F) 0-4
-10 °C (14 °F) 4-6 35 °C (95 °F) 4-6
10 °C (50 °F) 6-18 30 °C (86 °F) 6-18
-10 °C (14 °F) 18-24 35 °C (95 °F) 18-24

48 hour domestic express small package freight transportation
Winter Profile Winter Profile
Cold shippigg and cold receiving Hot shipping and hot receiving

Temperature Hours Temperature Hours
18 °C (65 °F) 0-4 22 °C (72 °F) 0-4
-10 °C (14 °F) 4-6 35 °C (95 °F) 4-6
10 °C (50 °F) 6-42 30 °C (86 °F) 6-42
-10 °C (14 °F) 42-48 35 °C (95 °F) 42-48

 

72 hour international expedited airfreight transportation

 

 

 

 

 

 

Cold shipping and cold receiving Hot shipping and hot receiving
Temperature Hours Temperature Hours
18 °C (65 °F) 0-4 22 0C (72 °F) 0-4
-10 °C (14 °F) 4-10 35 °C (95 °F) 4-10
10 °C (50 °F) 10-66 30 °C (86 °F) 10-66
-10 °C (14 °F) 66-72 35 °C (95 °F) 66-72
Cold shipping and cold receiving Hot shipping and cold receiving
Temperature Hours Temperature Hours
18 °C (65 °F) 0-4 22 °C (72 °F) 0-4
-10 °C (14 °F) 4-10 35 °C (95 °F) 4-10
10 °C (50 °F) 10-42 30 °C (86 °F) 10-42
22 °C (72 °F) 42-66 35 °C (95 °F) 42-66

35 0C (95 °F) 66-72 -10 °C (14 °F) 66—72

 

12

2.3 Test methods for thermal insulation quality of packages

Several test methods were proposed to estimate the thermal insulation
quality of packages.

ASTM 03103-92 is the standard test method for thermal insulation quality
of packages (ASTM 1999b). This test method covers the determination of the
thermal insulation quality of a package and its enclosed packaging for a
particular temperature difference between the packaged item and the outside
environment. It could be used for testing packages with and without various
internal refrigerants and with or without interior packaging. However, this
complicated procedure requires thermocouples to determine the temperature
profiles and R-value of the insulating package.

A more practical and simple procedure to determine R-value was
proposed by Burgess (Burgess, 1999). This procedure, so called the ‘ice-melt
test’, is based on the principle that 1 lb of ice at its melting point requires 144
BTU of latent heat to melt. In this method, an insulating package with a known
quantity of ice inside is stored for a designated period in a constant temperature
environment. By comparing the remaining quantity of ice to the original one, the
ice melt rate is determined. Using this rate, the rate of heat flow can be obtained
and subsequently the thermal resistance can be determined. This procedure
requires no special equipment or temperature sensors, but is still a very effective
method to determine the insulating ability of a package. In this work, the ice melt
test method was adopted for its simplicity, easy of testing, and low cost

performance (do not need any machine or sensors).

13

In general, the thermal resistance of packages is measured by an “R-
value” (R stands for resistance). The R-value is the reciprocal of the thermal

conductivity and is measured by following equation (Burgess, 1999):

AxAT _AxAT
meltrate x latentheat Q

 

R — value =

(1)

where R-value is the thermal resistance of container wall in ft2 -"F - hr/BTU , A is
the inside surface area of package in frz, AT is the temperature difference
between outside air and refrigerant used in °F, and Q is the heat transfer rate in
BTU/h. The melt rate is the rate that ice melts and is equal to the weight of the ice
melted divided by the melt time (lb/hr). When regular ice is used, the latent heat

of the ice is 144 BTU/lb.

2.4 Simulation of the thermal insulation quality of packages

As mentioned above, a comprehensive model is very important for
designing efficient insulating packages. Proper simulation models can reduce
time-consuming efforts of preliminary specification, fabrication and subsequent
validation tests. It also eliminates over-packaging and the resultant unnecessary
costs. However, there were very few attempts to predict the ability of insulation
packages (Burgess, 1999; Kositruangchai, 2003; Stavish, 1984). Moreover,
these models were over simplified or missed very crucial factors of heat transfer
through the insulating packages. Stavish proposed a simple model to predict the

heat flow through the package using the concept of total thermal resistance

14

(Stavish, 1984). The total thermal resistance, R,, is the summation of the R-values
of the material and the air film resistances on the inside and outside of the

insulating package. This can be expressed as

Rt = Ri + Rmalerral + R0 (2)

where R, is the air film resistance on the inside surface of the package, Rmm, is
the thermal resistance of the insulating material and R0 is the air film resistance of
the outside surface of the package. However, this equation overlooked the effect
of radiation and oversimplified the thermal resistance of the inside and outside of
the box.

Burgess (1999) also proposed a model to predict the R-values of the
insulating package. In his report, the system R-value could be predicted using

following fitted equation

Rm," = 3.9th +1.5np + 3.2nf (i20% accuracy) (3)

where th is the average wall thickness in inches, np is the number of plain
surfaces, and nf is the number of aluminum foil surfaces. In this equation, the
system R-value is separated into three parts: the effects of conduction,
convection and radiation. Conduction depends mainly on the wall construction
(th). Convection refers to heat transfer between air and surfaces so the R-value
depends upon the number of surfaces in contact with air (np). Radiation refers to
the emission or reflection of infrared waves and can be a significant factor. The

number of aluminum foil surfaces (nf) is therefore important to R-value.

15

Comparing to the study by Stavish (1984), equation (3) has improved
the concept of the thermal resistance by including radiation as a key variable.
However, this equation also oversimplified the complexity of the heat transfer
phenomena through the insulating package. For example, the equation is linear
but radiation is known to be a strong function of temperature. Equation (3) uses
an average of thermal conductivity (k) of 0.022 Btu/h'fi. °F but k varies for different
materials (See Table 2). Moreover, the thermal conductivity of materials at lower
temperatures should have higher R-value due to lower conductivities at lower
temperatures. For these reasons, equation (3) has a temperature limitation for
applications of various types of insulating packages. In fact, this equation was
adopted to the study by Kositruangchai (2003) to predict various conditions of

insulating package but was not successful to predict the R-value in many cases.

16

CHAPTER 3

THEORETICAL DEVELOPMENT

3.1 British thermal unit (BTU)

The British thermal unit, abbreviated by BTU, is the unit for measuring
heat quantity in the English unit measurement system. One BTU is defined as the
amount of heat required to raise the temperature of one pound of water at its
maximum density, which occurs at a temperature of 39.1 °F, by 1°F (Kreith,
1973). One BTU is approximately equivalent to the following: 251.9 calories;

1055 joules; 107.5 kilogram-meters; 0.0002928 kilowatt-hours.

3.2 Basic theories of heat transfer in insulating package

Heat transfer is the transport of thermal energy from one region to another,
normally from the hotter to the colder. The net rate of heat flow is in the direction
of decreasing temperature (Hagen, 1999). Heat transfer can occur by three

distinct modes: conduction, convection and radiation.

3.2.1 Conduction heat transfer

Conduction is the transfer of thermal energy in a solid or a liquid from
higher temperature to an adjacent point of lower temperature. The heat flux, q, is
the heat transfer rate per unit area. It is proportional to the temperature difference

since the driving force of the heat flux is the temperature difference between two

17

bodies from one region to another region (See Figure 5). The proportionality

coefficient, k, is called the thermal conductivity (Middleman, 1998).

 

 

 

 

all!” W“ g; ”1*”
TI ”T2... '4'
k .
W6!-
k .
<—>
Ax

Figure 5. Diagram for conduction heat transfer through a wall

3.2.1.1 Heat transfer through a solid wall

In Figure 5, conduction heat transfer through a wall at steady state in one
dimension follows equation (4). Based on Fourier’s law of heat conduction (Kreith,
1973), the heat flows from higher temperature T1 through the surface of wall to

surface T2 (lower temperature).

q:

g _k_A_T=k(_T1;7:2 (4,

= Ax Ax
where q is the heat flux in BTU/hftz, Q is the heat transfer rate in BTU/h, A is the
cross-sectional area in fiz, k is thermal conductivity in Btu/hft°F, AT is the
temperature difference in °F, and Ax is the thickness of the material in ft. T, and
T2 are the surface temperatures of the each side of the wall.

Thermal conductivity, k, is the rate of heat transmitted through a unit area
of material in a unit time through its total thickness with a unit of temperature
difference between the two opposite sides (Turner, 1981) and indicates how well
a material conducts thermal energy (McCabe et al., 1993). Thermal conductivity
of metals has a broad range of values, about 10 BTU/frh' °F for stainless steel to
240 BTU/frh' °F for silver. Water has a thermal conductivity of 0.3-0.4 BTU/frh' °F.
Gases have the lowest thermal conductivities; for air, k is about 0.016 BTU/frh' °F.
Solids having low k values are used as insulators to minimize the rate of heat
flow. Examples are polystyrene foam, urethane foam, and fiberglass. These
porous insulating materials act by entrapping air and thus eliminate the air-flow
within the cell, as long as the diameter is less than 4 mm (Brody and Marsh,
1997). Their k values are about equal to that for air. Thermal conductivity is a
function of temperature. In general when temperature increases, k increases.
Rearranging the equation (4):

T __ T
q 5 (5)

k

 

19

then

_ A(T1-Tz) _ T1 -T2 __A_x
Q q k

R

 

(6)

where R is the “thermal resistance” in ftz' °Fh/BTU. R is a material property that
describes the resistance to the passage of energy or heat through a material
(Turner and Malloy, 1981). The thermal resistance of a homogeneous body of
uniform cross-section is the reciprocal of conductance as shown in equation (6)
(Turner and Malloy, 1981). High thermal resistance indicates more effective
insulation (ASHRAE, 1993). There are several pertinent properties of insulation
as well as thermal conductivity of structures: the service temperature, the relative

humidity, surface emissivity, reflectivity and absorptivity, density, form and water

transmission.

3.2.1.2 Calculation of effective area for conduction through the walls of the
insulating package

In the previous section, one-dimensional conduction was discussed. In
this case, the cross sectional area of the entire material was constant. However,
real insulating packages have three-dimensional geometries. There is a
substantial surface area difference between the outside and inside of the
package due to its wall thickness. This difference affects the overall conduction
heat transfer rate significantly. For this reason, choosing the correct area for use
in the conduction equation is very important. Assuming there is a box with

outside dimensions of L, x W, x D, and the wall thickness oft, the heat transfer

20

rate through the wall in the thickness direction can be described in the following

equafion:

$194003;T (7)

where Q is the heat transfer rate in btu/h, k is the thermal conductivity in btu/hfr°F,
T is the temperature in °F, and x is the distance from the outside of the wall in ft.

In this case, the area of the box can be written as a function of x as follows:
A(x) 2 2[(Lo — 2x)(W0 _ 2x) + (Lo " 2x)(Do '- 2x) + (W0 — 2x)(Do _ 2’0] (8)
Equation (8) can be rearranged to yield

A(X) = 24(x - x1)(x - J‘2) (9)

 

where x13 =

 

L0 +W0 +00 +1J(L0 —W0)2 +(L0 -D0)2 +(W0 —D0)2
6 ’6 2

At steady state, heat transfer rate Q is constant and equation (7) can be

rearranged to the following equation:

 

To _Q r dx
IT, ”710 A(x) ‘10)

Substituting equation (9) with equation (10), and integrating gives
Toxin-=9 1 lnx2(x1_t) (11)
k 24051 “152) X102 ‘1)

21

By rearranging equation (11), the effective heat transfer area, :4, can be

expressed as

24(x1 — X2)
3‘2 (x1 - t)

x1(X2 *1)

Z:

 

(12)

Geometric averages and arithmetic averages of the inside and outside
areas of the box were investigated as candidates for the approximation of the
effective surface area in equation (12). The BASIC program AREA.BAS was
constructed to answer this question (See Appendix A). The results were shown in
Table 4. As shown in Table 4, the geometric average of the inside and outside
areas successfully represented the effective area for the calculation of
conduction heat transfer through the wall of the insulating package. The
geometric average of the two areas of the insulating package was used for the

rest of the modeling procedure.

3.2.2 Convection heat transfer

Convection is the mechanism by which thermal energy is transferred
between a solid surface and a fluid moving over the surface. Heat transfer by
convection is proportional to the temperature difference and surface area and is

known as Newton’s law of cooling (Hagen, 1999), shown in equation (13),

q = hAAT (13)

where h is heat transfer coefficient in BTU/h'fiz' °F.

22

mmocxoé \ xon 9: Co 835 can xon 65 Co 3630 cc cmmcc>m 2585:? .0
035.25 \ xon $5 Co 2365 new xon 65 cc 2333 cc cacao; 2.:ch06 .9

AN: concave 3 “3233.3 $05.22: scum Comxw .m

 

 

 

some 58 83 2mm 0.2.8 4.on 38 ~23 3: 0.2 3:
Soc 3:» Nod 38 3mm 0.8m 38 83 3: 9m oar
<85 ENS Sad «.38 News «.48 38 BE 3: car of
83 983 mood 38m 38m $8 33 BE 3: ca 0.2
ammo Ev? Sod 3&2 New? v.8: EB? 83 m2 3: emu
«Nos 9%» wood 0.2; 0.9% 35 0.0% 83 q: our 3:
33 was. «8.0 0.on ©on 3% 3% 83 3 em 0m
55 as o <52 Sew as n <80 a Scam. $2 med. e R: 9a 2
mmmcxozfimi 98.35 muss 023:0 $9.on £85 522, £80..

 

AN: concave E «Em m>=octc 53> 00.963 2585B new comco>m oEoEocm ho comtmanO .v Each

23

The heat transfer coefficient, h, is not a thermal property (Hagen, 1999). It
is determined by several factors: type of fluid (liquid or gas), flow condition
(laminar or turbulent), forced or natural convection or phase change, free-stream
velocity, surface geometry and roughness, position along the surface and
temperature dependence of fluid properties.

Rearranging equation (13),

 

T — T
q: 11 2 (14)
‘1?
then, the thermal resistance R for convection is
R--1— (15)
h

In convection heat flow, the driving force is the temperature difference
between two bodies. The thermal resistance is reciprocal of the heat transfer
coefficient as shown in equation (15) (Turner, 1981).

Convection can be classified into two categories: natural convection and
forced convection (Kreith, 1973). If air currents are resulted from buoyancy forces
generated by differences in density, which are caused by temperature gradients
in the fluid mass, then this is called natural convection. If the currents are set in
motion by the action of mechanical device such as fan, pump or agitator so that

the flow is independent of density gradients, then this is called forced convection.

24

In most cases involving the distribution of the insulating packages, the packages
are located on the floor or on pallets without any mechanical air-circulating
device. For these reasons, most of the convection heat transfers to and from
insulating packages are by natural convection.

In this study, only natural convection was considered to represent
convection heat transfer. Two types of natural convection heat transfer affect the
insulating package. Since the insulating package is located on the floor of the
warehouse or the distributing vehicle, the outside of the package contacts a very
large volume of the air in the environment. Convection in this case is
called ”natural convection in unconfined space.” At the same time, the insulating
package usually has an air space between the inside of the package and the
product. Convection in this space is called “natural convection in enclosed
space.” In next two sections, distinct characteristics of these two heat transfer

phenomena are discussed.

3.2.2.1 Natural convection in unconfined spaces

Natural convection phenomena for various geometries were studied
intensively. For an isothermal condition, the average free convection heat
transfer coefficient he for flow over plate can be calculated by the following

equation (Kreith, 1973):

k
h =—-N 16
c L u ( )

25

where k is the thermal conductivity of air in BTU/hfr"F, L is a characteristic length
dimension in ft, and Nu is the dimensionless Nusselt number. The Nusselt
number is a function of two dimensionless numbers Gr and Pr, which are the
Grashof number and the Prandtl number, respectively. General relationships
among these dimensionless numbers can be expressed by the following

equafion:

Nu =a-(Gr-Pr)m (17)

where a and m are constants. Table 5 gives these constants for free convection
heat transfer to and from plates. The Grashof number and Prandtl number are

defined by following two equations:

 

3 2
Gr=££¥£ (18)
,u
C
Pr: 2” (19)

where L is the length or height of the plate in ft, p is the density of air in lbw/fr}, p
is the viscosity of air in cp, AT is the positive temperature difference between the
wall and the air in °R, k is the thermal conductivity of air in BTU/hfl°F, C,, is the
heat capacity in BTU/lbm'0F, ,6 is the volumetric coefficient of expansion of the air
in UK and g is the acceleration due to gravity, which is 32.174 x (3600)2ft/h2. All
the air properties are evaluated at the air film temperature T] = (Tw + Tb)/2. T...

and T), represent the plate temperature and the air temperature, respectively.

26

The conductivity of air, k, the Prandtl number, Pr, and the quantity ngfl/ #2 in the
Grashof number are functions of the air film temperature. The relationship

between these parameters and the air film temperature is shown in Figure 6.

Table 5. Constants for calculating the Nusselt number (Geankoplis, 1983)

 

Physical Geometry NGNP, a m

 

Vertical plates (Vertical height L < 3ft)

<104 1.36 1/5
104~109 0.59 1/4
>109 0.13 1/4

Horizontal Plates
Lower Surface of cooled plates 105 ~ 2 x 107 0.54 1/4
2x107~3x101° 0.14 113
Upper Surface of cooled plates 105 ~1 011 0.58 1/5

 

3.2.2.2 Natural convection in enclosed spaces

Natural convection heat transfer in enclosed spaces was investigated in
various reports (Jacob, 1949; Kreith, 1973; Holman, 1986). The Nusselt number
in an enclosed air space between parallel planes can be expressed by the

following empirical formulas:

he = g - N11,; (20)

27

Conductivity, k ( BTU/h'ft'oF)

Prandtl Number, Pr

95132412 (1 flit/23). x10'6

0.026

 

0.024 -
0.022 -
0.020 -
0.018 -
0.016 5
0.014 -

 

 

Y = 0.0131 + 2.6136E-05 x - 6.1955E-09 x2

 

 

0.012
0.730

I

0

l I I T

1 00 200 300 400 500

 

0.720 -
0.710 -
0.700 -
0.690 5

0.680 -

 

0.670

Y = 0.07206 - 1.8453E-04 x + 3.5979E-07 x2 - 2.9967E-10 x3

 

 

5.000

T I

0 1 00 200 300 400 500

 

4.000 -

3.000 -

2.000 -

1.000 .

 

0.000

Y = 0.1794 + 4.1794 6'0-0095043 x

 

 

 

1 00 200 300 400 500

Temperature (0F)

Figure 6. Relation between temperature and the physical properties of air

28

where k is the thermal conductivity of air in BTU/hfroF, 6 is the thickness of the air

space between the planes on either side in ft and Nux is Nusselt number in an

enclosed air space.

For vertical planes (Jacob, 1949),

_ 0.186r61/4(-§-)‘“9 for 2000 < Gr5 < 20,000

Nua = (21)

 

L

0.0656r1/3(£)-“9 for 20,000< Gr <11 x106
5 5 6

where the Grashof number in enclosed space, Gr5, is defined by

2 1 3
Gra = p gm; 5 (22)
,u

where AT’ is the temperature difference between the planes on either side of the
enclosed space in OR and 6 is the thickness of the air space between the planes
on either side in ft. L represents the height of the air space in ft. All other symbols
of physical properties are the same with those of natural convection in
unconfined spaces. When the Grashof number is below 2000, the heat transfer is
dominated by conduction, so that Nuxis 1.0 for this region.

Nusselt number in horizontal air space with heating from below can be

also presented in the following equation (Jacob, 1949):

0.19595“4 for 104< Gr5<4x105
Nua = (23)

0.06805?3 for 4 x 105 < Grg

29

In a horizontal fluid layer with heating from above, heat transfer is by

conduction only.

3.2.2.3 Heat transfer balance by convection in enclosed spaces
Assuming two surfaces a and b form an enclosed space (See Figure 7),
the steady state heat balance between the two surfaces can be described by

following equation:

Qc,b—>a = hb Ab (Tb - Tair) = haAa (Tair '— Ta) (24)

where Qab—m is the convection heat transfer rate from surface b to a in btu/h, ha is
the convection heat transfer coefficient of the surface a in btu/ft’h °F, 11), is the
convection heat transfer coefficient for surface b in btu/fr’h °F, A0, is the area of
surface a in fr’, A), is the area of surface b in ft), Ta is the temperature of surface a
in °F, T), is the temperature of surface b in °F, and Ta), is the temperature of the air
in the space between surfaces 0 and b in °F.

Equation (24) can be rearranged in terms of Ta and T), by eliminating Tam
which results in the following equation:

_ haAa ‘thb

— T—T 25
0.1.... haAa+thb(b .) < >

 

30

Surface Surface
3 b

Air Space

 

F rl

Figure 7. Diagram of the enclosed air space between the surfaces a and b

Unfortunately, ha and h), in equation (25) are not usually known. From
equation (20), however, one can get the heat transfer coefficient, ha), for enclosed

spaces where Aa=Ab. In this case, equation (25) can be rearranged to yield

ha oh),
ha +hb

 

Qc,b-—>a = AMT}; — Ta) = habAb(Tb - Ta) (26)

If the temperatures of surface a and b are not so different, and the size
and shape of surface a and b are similar, then ha a hb. Eliminating Aa and
substituting A), with ha), one can get the relationship, ha = hb = 2hc. Using this

relationship, equation (26) can be rearranged to yield

31

Aa'Ab
Aa-I-Ab

 

Qc,b—)a = 2 ' he (Tb - Ta) (27)

In this case, heat transfer coefficient by convection, h 1.1,...“ can be

described as follows:

h, = Qc,b—)a = 2 ' hc ° Ag
8”“ Ab(Tb -T,,) Aa + A),

 

(28)

3.2.3 Radiation heat transfer

Radiation is thermal energy transfer by electromagnetic waves traveling at
the speed of light (2.998x108 mls). Radiation does not require a medium. 80 it
can transmit energy through a vacuum (Middleman, 1998). Thermal radiation lies
in the range of wavelengths from about 0.1 to 100 um (Hagen, 1999).

If there are two surfaces a and b, and the area of the surface a is smaller
than the larger surface which encloses it, radiation heat transfer rate from larger
surface b to smaller surface a, Q,,b_,a, can be described as the following
equation (Kreith, 1973):

_ Ab0'(Tb4 _Ta4)

Qr,b—)a _ _A_b_i-+ 1_€b (29)

Au Ea 5b

 

where an—m is a radiation heat transfer rate from surface b to a in BTU/h, as is the

emissivity of the surface a, ab is the emissivity of the surface b, 0' is the Stefan-

32

Boltzmann constant, which is 0.1714x10'8 BTU/hftz‘ °R, and Ta, T), are the absolute

temperatures in °R. Then, heat transfer coefficient by radiation is

9,-1.3. _ am,“ - Tab/(Tb — Ta)

 

 

 

h = _ 30
Au 5a 3b
if Aa=Ab, then
0(Tb4 - T04)/(Tb - Ta)
hr,b—>a = 1 l (31)
-—- + —— —1
8a 81,
if Aa< <Ab, then
4 4
8 0' T ——T
.3... = a ( b a ) <32)

 

(Tb T Ta)

3.2.4 Combined radiation and convection

In many practical problems, except for heat transfer at extremely high
temperatures, it is convenient to include radiation in a thermal network involving
convection and conduction. In this case, the total heat transfer coefficient where

convection, conduction and radiation occur simultaneously can be described as:

h = hc + h, (33)

where h is the effective convection heat transfer coefficient, he is the pure
convection or pure conduction heat transfer coefficient, and h, is the radiation

heat transfer coefficient. For natural convection in unconfined spaces, the overall

33

heat transfer coefficient can be written by substituting he and h, using equation

(32).

ado-(T124 T T04)

h =h +h =h +
C r a (Tb-Ta)

 

(33)

In similar manner using equation (28) and (30), the overall heat transfer

coefficient in enclosed spaces can be written for following two cases:

 

 

IfAb>Aa,
h _h '1']? _2h A0 +G(Tb4TTa4)/(TbTTa) (34)
c r c,a—>b Aa+Ab £2._I_+_1_-_8£
Aa ea 8b

IfAb=Aa,

4 4
h =2hc a—)b Aa +O(Tb Ta )/(Tb Ta) (35)

’ Aa+Ab 1 1
——+—-1
8a 8b

34

3.3 Mathematical models for heat transfer through multi-layered walls

In the previous section, basic theories of heat transfer and their application
to insulating packages were discussed. As mentioned earlier, the main objective
of this study was to develop an effective model which could successfully predict
the performance of an insulating package in distribution. The first step to
accomplish this objective is to obtain an accurate thermal resistance of the wall
of the insulating package.

If the wall consists of one material, it is very simple to calculate the
thermal resistance of the package. For example, the thermal resistance of pure
EPS is the reciprocal of its thermal conductivity, 0.022 BTU/hfr"F. However,
insulating packages used in many applications employ much more complicated
structures. Many insulating packages use a multi-layered structure with
combinations of several insulating materials. In many cases, these layered
materials are loosely fitted to each other to obtain extra thermal resistance by
entrapping air between the insulating materials. A loose-fitting EPS foam jacket
inside a corrugated box is a good example of these efforts. In this case, even
though thermal conductivities of each of the insulating materials are known, it is
still very difficult to estimate the overall thermal resistance of entire multi-layered
structure because of the complexity of the heat transfer mechanisms through the
‘between the wall’ air space.

To explain these phenomena and obtain an overall thermal resistance of
the multi-Iayered insulating materials, all of the basic theories used in the

previous part must be combined to yield a comprehensive mathematical model.

35

For these reasons, a mathematical model of the wall by itself, not the entire
package, was first developed. An effective thermal resistance of the wall was
derived and later used in the model of the whole package.

The one-dimensional diagram of a multi-Iayered wall is shown in Figure 8.
At steady state, the heat transfer rate through each solid layer is the same as the
heat transfer rate through the air space between the solid layers. Assuming an n
layered wall with n-l air spaces between them, one can establish a heat transfer

balance with following equations.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

i=1 ' i=2 i=n-l i=n
Tom 4:? Tin
hl kl h2 h2 k2 h3 hn-l kn- I h“ hn kn hn+1
T1 T2 T3 f I; T4 T2(n-l)—l T2(n-1) T2n-1 Tzn
I I I | l I
I d1 I dairl I (12 I ' dn-l ' dairn-l I dn I

Figure 8. Diagram of heat transfer through multi-layered walls in contact with
inside and outside air

36

for i=1

k
%=h1<To..—T1)=-im —T2)=h2(T2—T3) (36)
for i= 2 to 21-1
k .
hr (Tm—1) — T2i—1) = d—l-(T 21—1 - T21) = hr+1(T21 - T2141) (37)
l
for i=n
k
hn (T2014) — Tzn—i ) = d—"(Tzn—l - Tzn) = hn+l (Tzn — Tm) (38)

n

where d,- is the thickness of the ith layer of solid material, dair) is the thickness of
the i th layer of air space, k,- is the conduction heat transfer coefficient of i th layer
of solid material and h,- is the effective convection heat transfer coefficient for the i
th air space. h,- is derived from equation (35), which includes convection and
radiation.

In this system of equations, there are 2N +1 unknowns (Q, T1, T2, T2)“,
TZN) and 2N+1 equations. However, this system of equations is highly non-linear.
The convection and conduction coefficients, which are essential to solving for the
unknown temperatures, are functions of the temperature itself. Moreover,
radiation heat transfer coefficients include the terms of temperature powered by
four. To solve this non-linearity problem, an iterative procedure should be
adopted to solve this problem. The diagram of this procedure is shown in

Figure 9.

37

 

INPUT

73'", Tout, 61,8 and 6101.738

 

 

 

Initial Guess of Temperatures

Ts

 

 

 

 

Calculation of Heat Transfer
Parameters Using Temperature

k’s, h’s

 

 

 

Solve the Matrix of the System of Heat
Balance Equations

 

 

 

Obtain New Temperatures and Q

 

 

 

     

Is the difference between New and

 

 

Old Temperature is smaller than the
criteria?

No

       

 

Solution of the Temperatures and Q

 

 

 

Figure 9. Iteration procedure to solve non-linear problems.

38

Using this algorithm, the system of equations (36), (37) and (38) can be

solved in terms of the heat transfer rate, Q, to yield

Q _ Tl TTn (39)

:4“ . d,
Z—+
r=rki i=1

n-l 1

hi

In this case, the thermal resistance of the wall, Rm”, and the conductivity, km”,

of the entire wall structure can be calculated by equations (40) and (41)

respectively.
n d "—11

Rwall = Z4 + — (40)
._ k. ._ h-
[—1 I 1—1 I

n n—I
Zdi + Zdairi
kwall = T1 i=1
n n—l
di 1

25* 17.-

i=l i=1

(41)

 

39

3.4 Mathematical models for heat transfer through entire insulating
package

In the previous section, a mathematical model of one-dimensional heat
transfer through a wall was discussed. The effective thermal resistance of the
wall was obtained. If heat transfer through the insulating package is analyzed in
one-dimensional way, the model developed in previous section can be used to
explain heat transfer through the wall. In the real world, however, the wall is
three-dimensional. In the one-dimensional model, it is assumed that the heat
transfer is conducted through a series of the planes with same cross sectional
area and geometry. In real packages, however, the situation is different. If the
insulating package has a substantial wall thickness, there is a significant area
difference between the outside and inside of the wall. The difference in areas
should be considered in calculating the conduction heat transfer through the wall.
In addition, the geometry of the ice container might be totally different than the
geometry of the inside of the insulating package. This might result in different
convection and radiation heat transfer coefficients, which could affect the overall
heat transfer rate in the air space between the inside of the package and the ice.

Another factor which should be considered is contact between the
insulating package, ice container and the ground. In most cases, the insulating
packages sit on the ground and the ice container sit on the bottom of the
package.

For these reasons, three models were developed, which included more of

these factors in estimating heat transfer in real insulating packages. These

40

models are used to predict the results of ice melt tests done on various insulating
packages. For this reason, all three models consider the product inside to be ice

at its melting point in some containers.

3.4.1 Mathematical model considering the geometry of insulating package

A diagram of the model of heat transfer considering geometry differences
among the outside of the box, the inside of the box and ice container is shown in
Figure 10. This model is the simplest model of three models proposed. In this
model, contact between the ice container and the insulating package, and
between the outside of the package and the ground is ignored. Only the
convection and radiation are considered between the insulating package and the
ice container or the environment. In this model, a heat balance between the

outside of the package and the environment can be written as follows:
Q=h1Al(Too -T1) (42)

where Q is the overall heat transfer rate in BTU/h, A1 is the outside surface area
of the insulating package in ftz, and T0,, and T1 are the temperatures of the
outside air and the outside surface of the insulating package in °F, respectively.
h] is the combined convection and radiation heat transfer coefficient for the

outside surface of the package. From equation (33), h can be described as

41

A1! T1

 

Wall of the insulating package

 

I

Outside of the Box

 

 

 

 

 

 

 

Inside of the Box

 

 

 

 

 

 

 

Ice Container

 

 

 

 

 

 

 

Figure 10. Diagram of the model considering the geometry of the insulating
package.

42

slam." 41“)
(Too T Tl)

 

hl = hc,l + hr,l = hc,l + (43)

where he, 1 is the convection heat transfer coefficient for the outside surface area,
51 is the emissivity of the outside surface and a is the Boltzmann constant.
Heat transfer through the package wall is expressed as

T —T
] Ax 2 = h2A2(T2 TTice) (44)

 

Q = 11141066 41) = kIZJAIAZ

where k); is the effective conductivity of the wall in BTU/”Ffth, Ax is the thickness
of the wall in ft, A 2 is the inside surface area of the insulating package in ft", and
Tice and T2 are the temperatures of the ice container and the inside surface of the
insulating package in °F, respectively. Tice is the melting point of the ice used to
do the ice-melt tests. It is usually 32 °F. h; is the combined convection and
radiation heat transfer coefficient for the inside surface. From equation (34), h;
can be described as

A3 + 002“ - nab/(T2 — T...)

A2+A3 fli+l-82
A3 83 82

 

’12 = 2hc,2—>3 (45)

where hc,2_,3 is the convection heat transfer coefficient of the enclosed air space
between the inside of the package and ice container, and A 3 is the surface area

of the ice container in ft).

The heat transfer rate to the ice container can be described as

Q = (1343 (T2 - Tree) (45)

43

where A 3 is the surface area of the ice container in fr’ and h3 is the combined
convection and radiation heat transfer coefficient for the ice container. h3 can be
described as

A2 + om“ — Tish/(T2 — T...)
A2 '1' A3 1 :11 1T82

 

’13 = 21162—43 (47)

83 A2 82
In this system of equations, there are three unknowns: Q, T1 and T2. As
with the one dimension multi-Iayered wall, the system of equations is non-linear.
The conductivity and the convection and radiation coefficients, which are needed
to solve for the unknown temperatures, are functions of the temperature itself.
The iterative procedure in Figure 9 was introduced to solve the problem.
Solving the system of equations (42), (44) and (46) for Q in terms of T0,,

and Tice gives

1
Q=( 1 Ax 1 )(TooTTice) (48)

+ +
hlAl 1021/4142 11242

 

 

In this case, system R-value is

 

A2(Too TTice) 1 Ax 1
R t = = A2( + + )
sy‘ 3’" Q hl A1 k12 ,IAIA2 hz A2

 

 

(49)

The predicted heat transfer rate, Q, and system R—value will be
compared to experimental heat transfer rates from various insulating packages.
The degree of agreement between the two will determine whether or not the

assumptions made in the model are correct.

44

3.4.2 Mathematical model considering the contact between the Ice
container and the inside of the insulating package

The next model included contact between the ice container and the
insulating package. The heat transfer between the bottom of the ice container
and the inside of the package is treated by the concept of the contact resistance
(Cooper et al., 1969; Fried, 1963; Holman, 1986).

Heat transfer through the interface between two solid materials is limited
by the “contact resistance” at the interface. The resistance is determined mainly
by the surface roughness of the two materials (See Figure 11).

Every solid material has some roughness on its surface. Due to surface
roughness, there are two mechanisms controlling heat transfer at the interface.
One is the limited area corresponding to solid-to-solid conduction at the spots of
contact and the other is conduction through entrapped gases in the void spaces
created by the spots of contact. The major resistance of contact area is
conduction through the gases because the thermal conductivity of the gas is
quite small in comparison to that of the solid. The energy balance between two

solid surfaces, a and b, is given by the following equation

 

 

T -T T —T
Qcontacl = ka 1 a = hcontact (Ta __ Tb) = kb b 2

(50)
A Axa Axb

45

 

(A)

 

Q Q
9 Material (1 Material b +

AXa AXb

 

 

 

 

 

 

 

 

 

Figure 11. Diagram of thermal contact resistance between surface a and b:

(A) physical situation; (B) temperature profile (Holman, 1986).

46

or
T1 - T2
Qcontact = Axa 1 Axb
+ +
kaA hcontactA kbA

 

(51)

 

where 1/hcomaaA is called the thermal contact resistance and hcomac, is called the
thermal contact coefficient. However, most studies on the thermal contact
resistance dealt with the thermal conductance between the engineering materials
under high pressures (>10 atm). In insulating packages, the package sits on the
floor and the ice container sits on the bottom of the package. The pressure in this
situation is normal air pressure or 1 atm. In this case, it may not be appropriate to
use the thermal conductivities in these previous studies.

A diagram of the model of heat transfer considering the contact between
the ice container and the inside of the insulating package is shown in Figure 12.
This requires a separation of the inside surface area of the insulating package
into two parts: the area directly in contact with the ice container and the area in
contact with the air inside of the package. Therefore, the heat transfer balance is
divided into two parts explaining each situation separately.

From Figure 12, the heat balance between the bottom of the ice container

and the insulating package can be written as

 

T —T
Qcontact : k24\/ A2A4 2 Ax 4 = hcontactA4 (T4 T Tice) (52)

where Tice is the temperature of the ice container, T2 is the temperature of the

outside of the package, T4 is the temperature of the inside of the box, k24 is the

47

Outside of the Box:
ut=A1+A2

 

Inside of the Box

 

Ice Container :

A6=A4

 

Figure 12. Diagram of the model considering the contact between the Ice
container and the inside of the insulating package

48

thermal conductivity of the wall between the surface 2 and 4, and Ax is the
thickness of the wall.

Contact between the ice (or ice container) and the relatively warm air of
insulating package generates film of water, which promotes the thermal
conduction between the two surfaces. In this case, the thermal conductance
between the ice container and the inner surface of the insulating package should
be substantially bigger than the thermal conductivity of the insulating material
(hwmac, >> k24 / Ax). Then equation (52) can be rearranged to yield

_ T2 TTice .. Tz TTice _ lam/42444 53
Qcontact — Ax 1 = Ax — (T2 T Tice) ( )

Ax
+
1‘24 4244 hammer/14 k241/ 4244

 

 

The heat balance between the outside surface of the package and the

outside air is similar to the previous model.

Q = hour/10141000 ' Tout) (54)

where A0,), is the outside area of the insulating package in ftz, and Tc,o and Tow
are the temperatures of the environment and the outside of the insulating
package in °F, respectively. how is the combined convection and radiation heat
transfer coefficient for the outside surface of the package. In similar way to
equation (43), how can be described as

8100004 T Tout4)
(Too T Tout)

 

(55)

hour = hc,out

49

where hm“, is the convection heat transfer coefficient for the outside surface
area, 81 is the emissivity of the outside surface and 0' is the Boltzmann constant.

The heat balance between the two surfaces can be written as follows:

 

Q = hour Aout (Too T Tout) = Q1 + Q2 (56)
T - T
01 4131/4143 ‘Ax 3 (57)
T2 T Tice
Q2 = Qcontact = k24\1 A2 A4 T (58)

where Q, is the heat transfer rate through the area contacting the air inside the
package, Q2 is the heat transfer rate through the area contacting the ice
container, kg and k2.; are conductivities of the insulating materials in BTU/°Ffrh,
A3 is the inside area contacting the inside air in ft2 and A4 is the inside area

contacting the ice in frz. Q can also be written in the following expression:

T -T
Q1 = kI3\/AIA3 flATl = h3A3(T3 TTice) (59)

where h; is the combined convection and radiation heat transfer coefficient for
the inside surface of the package. In similar way to equation (45), h; can be

described as

 

A5 + 0(T34 TTice4)/(T3 TTice)
A3 + A5 fii+ 1-83
As 85 83

’73 = 2hc,3—>5 (60)

50

where he; _, 5 is the convection heat transfer coefficient of the enclosed air space
between the inside of the package and ice container, and A5 is the surface area
of the ice container. In this system of equations, there are three unknowns: Q, T.
and T3.

From equations (57) and (59),

Ql= 1 1 Ar (Tout—Tree) (61)

+
’73 A3 k24\/142 A4

Then, equation (56) can be rearranged to yield

 

 

k AA
Q=Ql+Q2=( 1 1 Ax + 2“"sz 4)<Tou,—T,-ce) (62)

+
h3A3 I‘m/«41143

 

Solving equations (56) and (62) for Q in terms of T00 and Tice gives

Q=(

 

1 1 )(Too _ Tice) (63)

+
hour Aout 1 + [(24 ‘/ A2 A4

1 Ax Ax

+
h3A3 I‘m/141143

 

 

In this case, system R-value is expressed as

_ Ain (Too ' Tice) 2 Am + Ain

R , _
sys em Q houtAout l + [€24 «142144
1 Ax Ax

+
h3A3 I‘m/A1143

 

 

(64)

 

51

where A," = A3 + A4 . Algorithm in Figure 9 was used to solve the non-linearity of
the equations.

The predicted heat transfer rate, Q, and system R-value will be
compared to experimental heat transfer rate from various insulating packages.
The degree of agreement between the two will determine whether or not the

assumptions made in the model are correct.

3.4.3 Mathematical model considering all contact areas

The bottom of the insulating package also contacts the ground. In this
case, heat transfer through the bottom of the insulating package is dominated by
the conduction between two surfaces rather than the convection to air, which is
assumed to occur in the previous models. A diagram of the model of the heat
transfer considering all the contact is shown in Figure 13. In addition to the
previous model, this model considers contact between the outside of the
package and the ground. Heat transfer through the outside of the package can
be divided into two parts: heat transfer by contacting the outside air and heat
transfer through contact with the ground. Using symbols in Figure 13, the heat

balance between the bottom of the box and the ground can be written as

 

 

Qcontacl = hcontact A2 (Tg - T2) + hcontactAZ (Tg " T2) (65)
T2 - T4

hcontact A2 (Tg ’ T 2 ) = k25 V A2 A5 Ax (66)
,— T3 - T6

hcontactA3 (T g T" T3) = k36 A3A6 Ax (67)

52

Outside of the Box

 

A2. T2 A3. T3 A2, T2

Inside of the Box

 

Ice Container

 

Figure 13. Diagram of the model considering all contact areas

53

where Tg is the ground temperature, T; and T3 are the temperatures of the bottom
of the package, T4 and T6 are the temperature of the inside of the box, k25 and
k36 are the thermal conductivities of the insulating wall, and Ax is the thickness of
the insulating wall. A2 and A3 are the bottom surface area of the outside of the
package. A5 is the inside surface area contacting the inside air and A6 is the
inside surface area contacting the ice.

One can assume that the ground temperature is nearly the same as the
air temperature of the environment. In addition, the conductivity of the insulating
material, k25 and k36, is very low and the thickness of the insulating package,
Ax, is quite large. Then, due to the small heat transfer coefficients, 1:25 /Ax and
k36 / Ax, the temperature of the bottom of the package should be very similar to
the ground temperature (T2 ; T3 a T8 ). In this case, equations (65), (66) and (67)

can be simplified to

Tg -T6
Ax

 

 

(68)

T — T4
Qcontact = k25 JAzAs g Ax + k36\/A3A6

The contact between the ice container and the inside of the insulating
package can be solved in similar manner with the previous model. Considering
that the ground temperature is similar to the outside air temperature (Tg .2; Too ),

the heat transfer balances for the entire model can be expressed as follows:

 

Q = M10... —T1>+k25JA2A5 T317,“ + k36JA3A6 -7"°—;xZ-"C—e (69)

54

Tl-T4

 

 

 

 

 

hlAl(Too -T1) = k14\/A1A4 Ax (70)
T,o —T4 T1 —T4
kzsxlAzAs Ax +k14JA1A4 Ax = h4,5(A4 +A5)(T4 _Tice) (71)
Too ”Tice
Q = h8A8(T4 — Tice) + k36\/A3A6 T (72)
where
4 4
h1=hc1+£“’(T°° T1 ) (73)
, (Too -T1)
4 4
h4 = 2h 4 A8 + 0'(T4 ’TI'ce )/(T4 "Tice) (74)
,5 C, ,5—)8 (A4+A5)+A8 A4 +145 i+1’84
A8 88 84
4 4
h8 = 2"! 4 5—)8 A4 +A5 + 6(T4 “Tice )/(T4 "Tice) (75)
c, 9

 

(A4+A5)+A8 .1... As 1-84
88 A4+A5 84

 

In this system of equations, there are three unknowns: Q, T1 and T4. The
algorithm in Figure 9 was used to solve the non-linearity of these equations.
Solving the system of equations (69), (70), (71) and (72) for Q in terms of

T00 and Tice gives

1 k 1/AA
Q=( +36 36

l 1 Ax

+
1 1‘25 4142/15 ’78 A8
+
1 Ar Ax

+
hlAl 16144141144

)(Too _ Tice) (76)

 

 

 

In this case, the system R-value is

55

Am

1 [‘36 x/A3A6
+
1 1 Ax

+
1 + k25JA2A5 h8A8
1 Ax Ax

+
hlAl kl4\/A1A4

where Ain=A4+A5+A6.

Rsystem = (77)

 

The predicted heat transfer rate, Q, and system R-value will be compared
to experimental heat transfer rate from various insulating packages. The degree

of agreement between the two will determine whether or not the assumptions

made in the model are correct.

56

CHAPTER 4

EXPERIMENTAL DESIGN

4.1 Materials and methods

Commercial ice (cubed), obtained from local grocery stores, was used in
ice-melt tests. Various sizes of boxes and ice containers were used in this
experiment. The geometry of the boxes and containers were summarized in
Table 6.

The ice-melt test (Burgess, 1999) was used to measure the overall heat
transfer rate and System R-value of each insulating package. The experimental
procedure is as follows:

_S£p_1: Since the ice was stored about 0 °F, it was preconditioned by
letting it sit out in the open environment to raise its temperature up to its melting
point of 32 °F. A quantity of regular ice (as much as possible to reduce the error
associated with water on ice left from preconditioning) was placed into a
container. A metallic ice container, which can interfere with the calculation of the
system R-value by providing a reflective surface not associate with the package
itself, was avoided. Plastic containers were used in all experiments. The package
was loosely closed and stored at room temperature for several hours to ensure
that all the surface of the ice was covered by water. This meant that the ice was
reached its melting point. The water was drained and discarded after

preconditioning step.

57

Table 6. Summary specifications for insulation materials

 

Experiment Dimension of each Number Outsnde d1mens10ns of Envnronmental

 

. . of ice the boxc Tern erature
number Ice container container (in x in x in) (1F)
1 6 x 5 a 1 9.25 x 9.25x 9.5 62
2 6 x 5 a 1 9.25 x 9.25x 9.5 90
3 6 x 5 a 1 9.25 x 9.25x 9.5 104
4 4.5 x 2.5 a 1 9.25 x 9.25x 9.5 62
5 6x65" 1 13.0x12.5x11.0 62
6 6x65‘1 1 13.0x12.5x11.0 90
7 6x653 1 13.0x12.5x11.0 104
6 6x5a 1 160x125x110 62
9 6 x 6.5 a 2 22.5 x 15.75 x 10.75 62
10 6 x 6.5 a 2 22.5 x 15.75 x 10.75 90
11 6 x 6.5 a 2 22.5 x 15.75 x 10.75 104
12 6 x 5 a 1 22.5 x 15.75 x 10.75 62
13 6x6x2.25b 2 26.75x11.75x5.0 62
14 6 x 6 x 2.25 b 2 26.75 x 11.75 x 5.0 90
15 6x6x2.25b 2 26.75x11.75x5.0 104
16 6 x 6 x 2.25” 1 26.75 x 11.75 x 5.0 62
17 6x58 4 43.0x12x8.75 62
16 6x53 4 43.0x12x8.75 90
19 8x5“ 4 43.0x12x8.75 104
20 9.5 x 4.5 a 1 43.0 x 12 x 6.75 62

 

a. Round plastic bucket: diameter (in) x height (in)

b. Square plastic bucket: length (in) x width (in) x height (in)

c. Airliner® insulating materials by Cargotech with the wall thickness of 1.0 inch
were used for these insulating packages.

58

$932: After the preconditioning procedure, the ice container placed back
inside the package near the center and the package was sealed with tape to
make it relatively airtight. The day and time are noted and the package is
immediately placed in a room with various temperatures (See Table 6). The
package is allowed to sit in this environment for several hours depending on the
package.

m: At the end of this period, the day and time are noted, the box is
opened, the ice container removed and the water again drained out. The drained
water is collected and weighed. A Digital mass balance (Mettler PM 2000) was
used for measuring the weight of the ice and the package.

m Thereafter, the experimental heat transfer rate, Qexp, and R-value
were calculated using equation (78) and (79) as shown below:
Qexp = MeltRate x LatentHeat (78)

_ Am XAT

Re” " Qexp

 

(79)

where Rexp is the experimental thermal resistance of insulating package in
ft2'°Fh/BTU, A in is the inside surface area of box in ftz, and AT is the temperature
difference between the outside air and the melting point of the refrigerant used.
The melt rate is the rate that ice melts in during the experiment, which is equal to

the weight of the ice melted divided by the melt time in Ib/hr and 144 BTU/lb for

59

regular ice was used as latent heat of the ice. The procedure of the ice-melt test

was summarized in Figure 14.

Put the sample box in the conditioning room
At least 24 hours
Put ice in the bucket
l Sit for several hours

Drain water out and discard

ll

Put the bucket back
into the sample box

Seal the sample box

1

Place sample box in the conditioning storage

I

 

Measure melted water and note exposure time _

Figure 14. The procedure for the regular ice melt test

60

Precondition

Procedure

These experimental data were compared with the predicted R-values of
three mathematical models, which were shown in equations (49), (64) and (77).

These three equations are summarized in Table 7.

Table 7. Predicted R-values of three mathematical models

 

Model considering geometry of the insulating package

 

R — A- (~1— + Ax + 1 )
system m hl A1 klz T1242 112/12

 

 

Model considering the contact between the ice container and the insulating
package

 

Ain + Ain

R t =
sys em hour Aout 1 + 1‘24 \/ A2 A4
1 Ax Ax

+
h3A3 kl3\/AIA3

 

 

 

 

 

Model considering all contact areas

 

 

 

 

 

R t = Ain
sys em
1 + [‘36 JA3A6
1 + 1 Ax
1 + ’62st A2145 h8A8
l Ax Ax

+
hlAl k14x/A1A4

 

 

 

61

CHAPTER 5

RESULTS AND DISCUSSION

5.1 Characteristics of the heat transfer in multi-Iayered walls

Heat transfer through multi-Iayered walls was investigated by simulation of
the model proposed in equations (36) through (41). Heat transfer through a
homogeneous wall made of single material is relatively simple. It is pure
conduction using the conductivities of various insulating materials reported in
several previous studies. However, heat transfer through walls containing air
spaces is much more complex. Conduction, convection and radiation are
involved in each of the air spaces between the layers. In equation (40) and (41),
all these heat transfer mechanisms were included.

One of the major questions about heat transfer within the airspaces in
muIti-Iayered walls is whether heat is conducted by conduction or by convection.
In fact, this is determined by the Grashof number, Gr5, based on the thickness,
6, of the air space. According to previous studies (Jacob, 1949; Kreith, 1973),
heat transfer is dominated by conduction if Gra is below 2000. Otherwise, heat
transfer is dominated by convection. According to equation (22), Grg in confined
space is a function of the distance, 6, between the solid layers and the
temperature difference between the layers, AT. The average temperature
between the layers ((Ta + Tb)/ 2 in Figure 7), also affects Gra because ngfl/ #2
is a function of the average temperature between the layers. Using this equation,

the characteristics of heat transfer of ‘between-the-wall’ air space was simulated.

62

In most packaging applications, the distance between the two solid layers
is less than 1.0 inch and the temperature difference between this thin air space is
less than 10 oF. The BASIC program NGR.BAS in Appendix A was used for
simulation. The simulation results are shown in Figure 15 as a contour plot.
Average temperatures of 40 °F, 70 °F, 100 0F were chosen to represent the low
temperature, room temperature and high temperature, respectively. When GI“; is
higher than 2000, which is the upper right part of each graph, heat transfer is
dominated by convection. When the distance between the solid walls is smaller
than 0.5 inches, heat is transferred mainly by conduction regardless of the
temperature difference or the average temperature in the simulation range. In
most cases, the thickness of the air spaces in multi-layered walls is very thin (<
0.5 in). This means that the most of the heat transfer in ‘between the wall’ air

spaces is dominated by conduction combined with radiation.

5.2 Simplification of the calculation of the system R- value of multi- layered
walls

In Equation (40), the resistance of a multi-layered wall is composed of two
main factors: one is the resistance from the solid material and one is the
resistance of air space between the walls. Assuming that heat transfer in the air
space is dominated by conduction and radiation, which is very likely in many

insulating packaging materials,

63

 

 

 

 

 

 

 

 

0.0

 

Temperature Difference between Two Layers (°F)

 

 

 

 

0.0 0.2 0.4 0.6 0.8 1.0 1.2
Thickness of the Air Space (in)

Figure 15. The Grashof number with the distance between the layers and the
temperature difference between the layers

64

 

k .
h.- = h...- +16..- = 74% hr..- (80)

71 di n—l 1
Rwall = Rsolz‘d + Rair = Zk—Jr Z— (31)
[:1 i i4. hr’.
dairi ’

From equation (80), the resistance generated by the ith air space can be

represented by the following equation:

 

 

 

 

1 . dairi
1 1 hr 1' kair
R . . = _ = = ’ 82
W" h.- kair. + h . A “”1 ( )
dai’i m hm: air

According to equation (31), the radiation heat transfer coefficient is a
function of the emissivities of both surfaces. If one or both surfaces are foiled
with aluminum, the radiation heat transfer coefficient, hm: can be reduced
dramatically. In this case, the denominator of equation (82) is dominated by 1/ hm,
and the equation becomes a linear function of the distance between the walls
over the conductivity of the air (dair/kair). For these reasons, the effect of the
emissivity on the resistance by the air space was simulated. Three cases were
simulated: both sides of the walls foiled (81 = 005,62 = 0.05), one side of the wall
foiled (6‘1 = 005,62 = 0.95), and none of the wall foiled (81 = 095,62 = 0.95) . The
BASIC program, NUSSELT.BAS, in Appendix A was used for simulation. This

program is a modified version of NGR.BAS and contained the algorithm using

65

 

equation (20), (21), (22) and (23), which calculates the heat transfer coefficients
and resistance of the air space under various conditions.

In addition to the emissivities, the distance between the two solid layers
(<1.0 in), the temperature difference between the layers (<10 °F), the average air
temperature (40 0F, 70 oF, 100 °F) were also chosen as variables. Using the
simulation results, heat resistances in three cases of the radiation were plotted
versus the distance between the walls over the conductivity of the air (dair/kair)
as an independent variable.

As shown in Figure 16, the three cases of radiation resulted in three
distinctive patterns of resistance according to the independent variable, dair/kair.
Moreover, the resistances of the air space showed a linear pattern when the
emissivity of either walls is small (one or more walls are foiled).

The resistances of the air space in these two cases were linearized as
functions of dair/kair. Linear coefficients for the cases of two-walls- foiled and
one-walI-foiled are 0.88992 and 0.80724 respectively. The coefficients for
hyperbolic type heat resistance of equation (82) for no-waII-foiled air space were
obtained from the non-linear regression function using Sigmaplot 4.0. 1.0855
and 1.1366 were obtained for 1/ hm- of the numerator and denominator,
respectively. These results were converted into the formulation of system R-
values. From equation (81), the resistance generated by the solid material can be
expressed as follows. Assuming that the conductivity of the solid material is
0.022 BTU/°Ffrh, which is usually the case with insulating materials, the

resistance of the solid material in equation (81) can be described as:

66

 

Rair

 

 
 
   
    

 

 

0 R0,, with 31:0.95 and 32:0.95 é)
v Rair with 61=O.95 and 62=0.05
CI Rair with s1=0.05 and 62=0.05 1i"
6 ‘ <> Raj, without Radiation (Control) “a
4 _
2 -
N \‘ w“ 9“ “92“"? 312:1:- 21-. TALE-it, ‘1)ng r “(z-'0
' 3‘ .,.,.“\\599:9?»09".09'09'6“,Mum...mortuolttloqflrtfi
3’: c . . . . . .
O T T I
0 2 4 e

d/k

Figure 16. Resistance of the air space with various emissivities

67

 

 

N

N N

d- Id- 1 d-
R -= —'= ——'= ———=3.79 d 63
50’“ Zk, Elk-12 =0.02212 Z ( )

1:]

where di is the thickness of the ith solid wall in inches. The heat resistance of the
air space can be classified into three cases. Assuming that the conductivity of the
air is 0.016 BTU/°Ffrh, linearized heat resistance of the air space where both of

the walls are foiled can be expressed as:

 

 

”1 dairli 0. 88992 dairl, ”1
R - = 0.88992 =4.64 dairl- 64
“"1 E, k,- =: 0.016 12 E; ’ ( )

where N1 is the number of the airspaces with both sides of the walls foiled and
dairli is the thickness of the ith both-layers-foiled air space in inches. In similar
way, linearized heat resistance of the air space where one of the walls is foiled

can be described as follows:

N2
=4.202dair2,- (85)
i=1

dair2,- 20.80724 dair2-
= 0. 80724——— '
R0“ 2: k,- i; 0.016 12

 

where N2 is the number of the airspaces with one side of the layers foiled and
dair2,- is the thickness of the ith one-wall-foiled air space in inches. Heat

resistance of the air space without any foiled walls can be expressed using non-

linear regression.

68

dair3,~

 

 

 

1.——0855
N3 N3 -
1.,0855dazr3
Rair3 = Z :21 (86)
i=11.1366+ dair3,~ il= 0.2182+dair3,-
kair,

where N3 is the number of the airspaces with non-foiled air space and dair3,~ is
the thickness of the ith non-foiled air space in inches. In this case, the R-value of

the multi-layered wall can be described as follows:

N 3 1. 086dair3,
Rm,” = 3. 79ths + 4. 64tha1 + 4. 20tha2 + 2
0.218 + dair3,-

 

(87)

where ths is the total thickness of solid layer in inches, thal is the total thickness
of both-walI-foiled air space in inches, tha2 is the total thickness of one-waIl-foiled
air space in inches, N3 is the number of the airspaces without foiled layer and

dair3i is the thickness of the ith non-foiled air space in inches.

5.3 Evaluation of the simplified equation using computer program

In the previous section, a model was constructed to represent heat
transfer through a multi-Iayered wall structure. The model was based on heat
transfer theories, so it depends on the assumptions that could generate errors.
For example, the model assumed that the conductivity of the air was constant at
0.016 BTU/’Ffrh. In the real world, however, the conductivity of air is a function of
its temperature. These kinds of errors should be minor ones, but they might
accumulate and create huge disagreements between the model and real

situations. For these reasons, evaluation and refinement of the model is

69

 

necessary. To evaluate the model, a computer program including all parameters
which could affect heat transfer was constructed. Using equation (36), (37) and
(38), which represent conduction, convection, and radiation, the BASIC program
WALLBG7.BAS in Appendix A was constructed. 96 different cases of multi-
layered walls were arbitrary chosen with various constructions (See Appendix B,
C and D). Using these data, WALLBG7.BAS created R-values, total thickness of
solid layers, total thickness of both-wall-foiled layers, total thickness of one-side-
foiled layers, the number of the airspaces without foiled layers and the thickness
of each non-foiled air space (See Appendix E). R-values generated by this
program were compared with the R—values calculated using the simplified
equation (87).

The results are shown in Table 8. The errors in the R-values between the
computer model and the simplified model were less than 20%, which means that
the simplified model successfully predicts the system R-value. Except for four
cases, however, the percentage errors, which were calculated by (R-value by
computer program - R-value by simplified model)l(R-value by computer program)
x 100, showed positive numbers. This implies that the simplified model generally
underestimates the R-value. This model was refined by fitting each coefficient in
equation (87). First, the coefficients for the linear part of the equation (87) were
fitted by multi-variable regression. The coefficients for total thickness of solid wall
(ths), total thickness of both-walI-foiled air space (thal), and total thickness of

one-side-foiled air space (tha2) were obtained by multi-variable linear regression.

70

Table 8. Comparison of R-values in equation (87) with the computer model

 

 

 

 

R by Eq. (87) R by Eq. (87)
Data Set Rwa” Rwall °/o Error Data Set Rwa” Rwall % Error

1 5.891 5.685 3.5 49 8.910 8.986 -0.9
2 11.757 11.370 3.3 50 16.154 15.384 4.8
3 12.262 11.790 3.8 51 9.234 8.986 2.7
4 13.296 12.630 5.0 52 9.561 8.986 6.0
5 21.178 19.995 5.6 53 14.931 14.201 4.9
6 4.461 3.874 13.2 54 16.695 15.804 5.3
7 9.732 8.416 13.5 55 12.697 12.262 3.4
8 10.761 8.836 17.9 56 12.697 12.262 3.4
9 10.333 8.836 14.5 57 2.796 2.635 5.8
10 10.328 8.836 14.4 58 4.982 4.678 6.1
11 9.748 7.996 18.0 59 4.982 4.678 6.1
12 8.702 7.576 12.9 60 4.982 4.678 6.1
13 5.851 5.685 2.8 61 4.982 4.678 6.1
14 11.702 11.370 2.8 62 9.426 9.164 2.8
15 12.269 11.834 3.5 63 16.154 15.384 4.8
16 13.356 12.762 4.4 64 9.645 9.709 -0.7
17 21.206 20.303 4.3 65 11.893 11.191 5.9
18 4.569 4.094 10.4 66 14.931 14.201 4.9
19 10.077 8.988 10.8 67 5.851 5.685 2.8
20 11.161 9.452 15.3 68 11.702 11.370 2.8
21 10.676 9.452 11.5 69 12.089 11.712 3.1
22 10.673 9.452 11.4 70 12.269 11.834 3.5
23 10.106 8.524 15.7 71 18.784 18.153 3.4
24 9.004 8.060 10.5 72 3.475 3.034 12.7
25 5.851 5.685 2.8 73 5.331 4.835 9.3
26 11.702 11.370 2.8 74 6.382 5.652 11.4
27 12.269 11.834 3.5 75 5.331 4.835 9.3
28 13.356 12.762 4.4 76 4.982 4.678 6.1
29 21.206 20.083 5.3 77 6.164 5.398 12.4
30 4.562 4.050 11.2 78 4.757 4.401 7.5
31 10.047 8.900 11.4 79 8.910 8.986 -0.9
32 11.033 9.188 16.7 80 16.494 15.540 5.8
33 10.571 9.276 12.2 81 11.174 10.044 10.1
34 10.488 9.144 12.8 82 10.602 9.501 10.4
35 9.846 8.260 16.1 83 16.329 14.917 8.6
36 8.738 7.708 11.8 84 20.815 18.543 10.9
37 5.851 5.685 2.8 85 12.697 12.262 3.4
38 11.702 11.370 2.8 86 15.880 16.010 -0.8
39 12.089 11.712 3.1 87 3.510 3.166 9.8
40 12.089 11.712 3.1 88 6.410 5.740 10.5
41 18.784 18.153 3.4 89 6.341 5.476 13.6
42 2.796 2.635 5.8 90 4.982 4.678 6.1
43 4.982 4.678 6.1 91 6.453 5.740 11.0
44 4.982 4.678 6.1 92 10.853 10.053 7.4
45 4.982 4.678 6.1 93 17.205 15.986 7.1
46 4.982 4.678 6.1 94 12.688 11.157 12.1
47 4.982 4.678 6.1 95 14.649 13.267 9.4
48 4.588 4.322 5.8 96 16.484 14.907 9.6

 

71

The BASIC program CHOI5.BAS in Appendix A was used for the multi-
variable linear regression. Among 96 data sets (Appendix B, C and D), 36 data
sets, which contained no non-foiled layers, were selected for linear regression.
Values of 3.89001, 5.558905 and 5.305831 were obtained for the coefficients for
total thickness of solid wall (ths), total thickness of both-walI-foiled air space (thal),
and total thickness of one-side-foiled air space (tha2), respectively. Then the non-
linear part of the equation was fitted to data for the coefficient of the thickness of
each none-foiled air space. CHOI6.BAS in Appendix A, a modified version of
CHOI5.BAS, which contained an additional loop to solve the non-linear problem,
was used for regression. Out of the 96 data sets (Appendix B, C, D), 30 data sets,
which constructed only with none foiled walls, were used for the regression. For
non-linear part of the equation (87), 1.22859 and 0.233 were obtained for the
coefficients of numerator and denominator, respectively. With these results,

equation (88) could be written as follows:

N 3 1.229dair3,
Rm,” = 3.89716 + 5.56tha1+ 5.31tha2 + Z _
0.233 + da1r3,-

i=1

 

(33)

For the verification of the model, the R-values in Appendix E were again
compared with the R-values calculated by equation (88). As shown in Table 9,
the percentage error in R-value was within 10%. Furthermore, the percentage
errors were distributed positive and negative numbers, which implies that the

refined model significantly improved the accuracy of the simulation.

72

 

Table 9. Comparison of R-values in equation (88) with the computer model

 

 

 

 

R by Eq. (88) R by Eq. (88)
Data Set Rwa” Rwall % Error Data Set Rwa” Rwall % Error

1 5.891 5.851 0.677 49 8.910 9.406 -5.565
2 11.757 11.703 0.464 50 16.154 15.991 1.012
3 12.262 12.232 0.240 51 9.234 9.406 -1.860
4 13.296 13.291 0.040 52 9.561 9.406 1.625
5 21.178 21.260 -0.389 53 14.931 14.773 1.055
6 4.461 4.473 -0.280 54 16.695 16.560 0.809
7 9.732 9.926 -2.000 55 12.697 12.699 -0.016
8 10.761 10.456 2.834 56 12.697 12.699 -0.016
9 10.333 10.456 -1.186 57 2.796 2.762 1.210
10 10.328 10.456 -1.233 58 4.982 4.916 1.331
11 9.748 9.397 3.607 59 4.982 4.916 1.331
12 8.702 8.867 -1.900 60 4.982 4.916 1.331
13 5.851 5.851 0.001 61 4.982 4.916 1.331
14 11.702 11.703 -0.002 62 9.426 9.604 -1.889
15 12.269 12.257 0.092 63 16.154 15.991 1.012
16 13.356 13.367 -0.083 64 9.645 10.221 -5.973
17 21.206 21.438 -1.093 65 11.893 11.761 1.106
18 4.569 4.600 -0.671 66 14.931 14.773 1.055
19 10.077 10.256 -1.773 67 5.851 5.851 0.001
20 11.161 10.811 3.136 68 11.702 11.703 -0.002
21 10.676 10.811 -1.261 69 12.089 12.072 0.142
22 10.673 10.811 -1.288 70 12.269 12.257 0.092
23 10.106 9.701 4.010 71 18.784 18.761 0.125
24 9.004 9.146 -1.576 72 3.475 3.414 1.752
25 5.851 5.851 0.001 73 5.331 5.237 1.762
26 11.702 11.703 -0.002 74 6.382 6.321 0.958
27 12.269 12.257 0.092 75 5.331 5.237 1.762
28 13.356 13.367 -0.083 76 4.982 4.916 1.331
29 21.206 21.311 -0.495 77 6.164 6.059 1.702
30 4.562 4.575 -0.277 78 4.757 4.679 1.636
31 10.047 10.205 -1.574 79 8.910 9.406 -5.565
32 11.033 10.659 3.390 80 16.494 16.312 1.108
33 10.571 10.709 -1.313 81 11.174 11.054 1.074
34 10.488 10.633 -1.390 82 10.602 10.395 1.955
35 9.846 9.549 3.018 83 16.329 16.109 1.350
36 8.738 8.943 —2.352 84 20.815 20.594 1.063
37 5.851 5.851 0.001 85 12.697 12.699 -0.016
38 11.702 11.703 -0.002 86 15.880 17.251 -8.637
39 12.089 12.072 0.142 87 3.510 3.490 0.551
40 12.089 12.072 0.142 88 6.410 6.372 0.598
41 18.784 18.761 0.125 89 6.341 6.220 1.912
42 2.796 2.762 1.210 90 4.982 4.916 1.331
43 4.982 4.916 1.331 91 6.453 6.372 1.257
44 4.982 4.916 1.331 92 10.853 11.000 -1.353
45 4.982 4.916 1.331 93 17.205 17.030 1.015
46 4.982 4.916 1.331 94 12.688 12.526 1.279
47 4.982 4.916 1.331 95 14.649 14.630 0.126
48 4.588 4.518 1.510 96 16.484 16.218 1.614

 

73

 

The major drawback of this equation is that thermal resistance generated
by the solid layer cannot be differentiated between each solid material. Thus,
thermal resistances per inch of each solid materials, which is summarized in

Table 10, were added to the equation (88) to yield

N N 3 1.229dair3-
R = R-thso + 5.56thal + 5.31tha2 + '
”a” Z ’ ’ 2 0.233 + dair3,-

i=1 i=1

 

(39)

where N is the number of solid layers, R,- is the thermal resistance per inch for the

solid material of the ith layer and ths, is the thickness of solid material of i th layer.

 

5.

Table 10. Thermal resistance per inch for various insulating materials

 

Thermal resistance per inch at 70°F

 

Materials
R,- (hrfth’F/BTUin)
Air 5.56
Corrugated board 2.38
EPS foam 3.78
Polyurethane 4.62
VIP 10-30

 

Using this equation, the Rm,” values of various insulating materials were
calculated. For example, the Rm," of the gas-filled panel used for ice-melt test was
calculated in the following procedure. The thickness of the gas-filled panel used
in this study was 1.0 inch. The real structure of the gas-filled panel is very
complicated honeycomb type (See Figure 4). For the convenience of the

calculation, the structure was simplified with the assumption that the gas-filled

74

panel is composed of two ‘one-wall aluminum foiled’ air spaces. The thickness of
each air space is 0.5 inches in this case. Ignoring thermal resistance generated
by very thin layers of the solid materials, the approximate value or Rm,” of the

gas-filled panel is

R...” = 5.31 x tha2 = 5.31 x 1.0 = 5.31 hrftz‘oF/BTU
where tha2 is the total thickness of the one-layer foiled air space, which is 1.0
inch. This value is similar to the R-value per inch of gas-filled panel, 5.2

hrft2'0F/BTUin in previous study (Griffith et al, 1991).

75

5.4 Simulation of the ice-melt test using mathematical models

In previous sections, one dimensional heat transfer through the insulation
package was simulated. In the real world, however, the insulating package has
three-dimensional bodies and the complexity of three.dimensional geometry
could affect the overall insulating capacity of the package. Based on one-
dimensional model, a more comprehensive approach was performed to construct
a model which could explain all of the factors affecting the insulating capacity of
the packages. These models were evaluated by comparison to the results of ice-
melt tests on insulating packages with various materials and geometries.

The major purpose of this study was to construct a simple model which
could predict the complexity of heat transfer through an insulating package. For
this reason, the simplest model was constructed based on one-dimensional
multi-wall structure model. Thereafter, two more models were constructed, which
included contact resistance. Eventually, these models were compared to the
experimental data and the simplest model which successfully explained the ice-

melt test data was selected.

5.4.1 Construction of mathematical model considering the geometry of
insulating package

The detail of the model was described in theory part of this study and in
Figure 10. BASIC program KCARGO9.BAS in Appendix A was constructed using

equations (42), (43), (44), (45), (46) and (47).

76

 

5.4.2 Construction of mathematical model considering the contact between
the ice container and the insulating package

Using equations (55), (56), (57), (58), (59) and (60), the mathematical
model considering the contact between the ice container and the inner surface of
the insulating package was constructed. The BASIC program KCARGO16.BAS
in Appendix A was constructed for computer simulation. KCARGO16.BAS is a
modified version of KCARGO9.BAS with the addition of the equation
representing contact between the ice container and the inner surface of the
package. The details of the model are described in theory part of this study and

in Figure 12.

5.4.3 Construction of mathematical model of the insulating package
considering all contacts

Using equations (69), (70), (71), (72), (73), (74) and (75), the
mathematical model considering the contact between the ice container and the
inner surfaces of the insulating package, and between the outer surfaces of the
package and the ground was constructed. The BASIC program KCARGO17.BAS
in Appendix A was constructed for computer simulation. KCARGO17.BAS is a
modified version of KCARGO16.BAS with the addition of an equation
representing contact between the bottom of the insulating package and the
ground. The details of the model are described in theory part of this study and in

Figure 13.

77

 

5.4.4 Comparison of models with experimental data

The three models proposed in the previous section were evaluated by
comparing melt rates to ice-melt test results. In addition to the experimental data
in Table 6, several previous ice-melt test data (Kositruangchai, 2003) with
various geometries and insulating materials were used for the evaluation. The
experimental data from the ice-melting test are summarized in Table 11.

Using the three BASIC programs, the overall R-values (Rsymm) were
simulated for each case. The difference between the experimental and simulation
values for each experimental data was squared and added to obtain the sum of
the squared errors (SSE) for each model. The SSE’s are shown in Table 12. The
percent errors between the experimental data and the simulation data are also
presented in the same table. The first model considering geometry (not
considering any contact) predicted the overall heat transfer rate within 36% of the
actual value and the SSE was 75.15. Except one case, the model considering
the contact between the ice container and the inner surface of the insulating
package predicted the overall heat transfer rate within 23% of the actual value.
The SSE value of this model was 44.30, which is almost half of the previous
model. This means that the consideration of the contact between the ice and the
insulating package dramatically improves the accuracy of the model. The result
of the model considering all contacts was similar to those of the model

considering the contact between the ice and the insulating package.

78

 

Table 11. Summary of the experimental data and areas according to Figure 13.

 

 

DataSet A1 A2 A3, A6 A4 A5 A8 x Qexp Rexp
1 3.035 0.396 0.196 1.675 0.169 0.851 0.063 11.09 6.06
2 3.035 0.396 0.196 1.875 0.169 0.651 0.063 21.32 6.10
3 3.035 0.398 0.196 1.875 0.169 0.651 0.063 26.73 6.03
4 3.035 0.484 0.110 1.875 0.255 0.356 0.083 5.94 11.32
5 5.024 0.779 0.349 3.490 0.453 1.464 0.063 15.49 6.31
6 5.024 0.779 0.349 3.490 0.453 1.484 0.083 33.64 7.36
7 5.024 0.779 0.349 3.490 0.453 1.464 0.083 44.63 6.92
6 5.024 0.932 0.196 3.490 0.606 0.851 0.083 11.23 11.46
9 6.172 1.763 0.698 6.120 1.259 2.967 0.063 25.66 9.44
10 8.172 1.763 0.698 6.120 1.259 2.967 0.063 56.31 8.03
11 8.172 1.763 0.698 6.120 1.259 2.967 0.083 79.02 7.36
12 8.172 2.265 0.196 6.120 1.761 0.851 0.063 15.06 16.07
13 4.656 1.037 1.146 3.113 0.530 2.302 0.083 20.13 7.40
14 4.856 1.037 1.146 3.113 0.530 2.302 0.083 57.75 4.99
15 4.856 1.037 1.146 3.113 0.530 2.302 0.083 69.69 5.13
16 5.007 1.931 0.333 3.113 1.342 0.771 0.083 13.07 11.40
17 10.267 2.167 1.396 7.629 1.451 4.667 0.063 41.04 7.66
16 10.267 2.187 1.396 7.629 1.451 4.887 0.083 66.39 7.03
19 10.267 2.187 1.396 7.629 1.451 4.887 0.063 123.04 6.13
20 10.267 3.091 0.492 7.629 2.355 1.256 0.063 19.6 15.67
21* 5.356 0.665 0.196 5.090 0.616 0.759 0.013”I 45.15 5.23
22* 5.356 0.665 0.196 5.090 0.616 0.733 0.013‘3 43.1 5.46
23* 5.358 0.665 0.196 5.090 0.616 0.720 0.013a 42.17 5.60
24* 5.266 1.113 0.196 5.000 1.054 0.733 0.013a 45.66 5.45
25* 5.266 1.113 0.196 5.000 1.054 0.707 0.013a 43.64 5.73
26* 5.266 1.113 0.196 5.000 1.054 0.694 0.013“1 42.76 5.65
27* 6.125 1.054 0.196 3.667 0.554 0.612 0.125b 17.39 10.16
26* 6.125 1.054 0.196 3.667 0.554 0.812 0.125'3 17.46 10.12
29* 6.125 1.054 0.196 3.667 0.554 0.877 0.125b 19.04 9.26
30* 5.674 1.246 0.196 4.021 0.673 0.759 0.0633b 24.66 6.25
31* 5.674 1.248 0.196 4.021 0.673 0.746 0.0833b 23.93 8.51
32* 5.674 1.248 0.196 4.021 0.873 0.746 0.0633'D 25.5 7.96
33* 5.534 0.910 0.196 3.472 0.498 0.746 01093" 15.39 10.83
34* 5.534 0.910 0.196 3.472 0.498 0.746 0.1093c 15.69 10.62
35* 5.534 0.910 0.196 3.472 0.498 0.746 01093" 15.55 10.72
36* 5.263 0.656 0.196 3.472 0.498 0.765 0.0963d 12.96 12.86
37* 5.263 0.856 0.196 3.472 0.498 0.753 0.0963d 12.96 12.86
36* 5.263 0.856 0.196 3.472 0.498 0.772 0.0963d 13.62 12.24

 

* Data by Kositruangchai 2003. All Kositruangchai 's data was performed at 72 °F.
a. Insulating material is a corrugated board

b. Insulating material is an EPS foam
c. Insulating material is a polyurethane
d. Insulating material is a VIP

79

 

Table 12. Comparison of simulation results of three computer models with the
experimental data

 

Experiment

Simulation with

Simulation with Contact Simulation with Contact

 

 

 

Geometry (ice and package) (all contact)

RexP Rsystem 3223'. SE Rsystem E623" SE Rsystem iii/gr SE
1 6.06 7.16 -18.2 1.22 7.25 -19.6 1.41 7.14 -17.8 1.16
2 6.10 6.70 -9.9 0.37 6.78 -11.3 0.47 6.69 -9.7 0.35
3 6.03 6.50 -7.7 0.22 6.58 -9.0 0.30 6.49 -7.5 0.21
4 11.32 9.14 19.3 4.75 9.83 13.1 2.20 9.71 14.2 2.59
5 8.31 7.73 7.0 0.34 7.85 5.6 0.22 7.70 7.4 0.37
6 7.36 7.23 1.8 0.02 7.34 0.2 0.00 7.21 1.9 0.02
7 6.92 7.01 -1.2 0.01 7.12 -2.9 0.04 7.00 -1.1 0.01
8 11.46 9.06 21.0 5.77 9.54 16.8 3.71 9.39 18.1 4.32
9 9.44 8.20 13.2 1.55 7.97 15.6 2.18 7.77 17.7 2.78
10 8.03 7.69 4.2 0.12 7.47 7.1 0.32 7.30 9.2 0.54
11 7.36 7.47 -1.5 0.01 7.25 1.5 0.01 7.09 3.7 0.07
12 16.07 12.06 24.9 16.06 13.49 16.0 6.63 13.29 17.3 7.74
13 7.40 6.89 3.4 0.06 6.70 6.1 0.19 6.48 9.2 0.43
14 4.99 6.53 -35.7 2.95 6.32 -31.4 2.28 6.13 -27.4 1.73
15 5.13 6.36 -28.5 1.99 6.15 -24.2 1.44 5.97 -20.6 1.03
16 11.40 10.73 2.4 0.07 12.24 -11.4 1.56 12.03 -9.4 1.07
17 7.66 7.86 ~2.6 0.04 7.48 2.3 0.03 7.27 5.0 0.15
18 7.03 7.41 -5.3 0.14 7.04 0.0 0.00 6.85 2.5 0.03
19 6.13 7.20 -17.5 1.15 6.84 -11.6 0.50 6.67 -8.8 0.29
20 15.87 11.26 29.1 21.29 12.19 23.2 13.59 11.95 24.7 15.36
21 5.23 5.87 -12.2 0.41 5.29 -1.2 0.00 4.82 7.8 0.16
22 5.48 6.00 -9.5 0.27 5.39 1.6 0.01 4.92 10.2 0.32
23 5.60 6.06 -8.3 0.22 5.44 2.8 0.02 4.97 11.3 0.40
24 5.45 6.22 -14.1 0.59 5.73 -5.0 0.08 5.23 4.1 0.05
25 5.73 6.36 -11.0 0.40 5.84 -1.9 0.01 5.33 6.9 0.16
26 5.85 6.44 -10.1 0.35 5.90 -0.9 0.00 5.39 7.8 0.21
27 10.16 8.57 15.7 2.54 9.10 10.4 1.11 8.99 11.5 1.37
28 10.12 8.57 15.3 2.41 9.10 10.0 1.03 8.99 11.2 1.28
29 9.28 8.36 9.9 0.84 8.84 4.8 0.20 8.72 6.0 0.31
30 8.25 8.01 2.9 0.06 8.76 -6.1 0.26 8.58 -3.9 0.11
31 8.51 8.06 5.2 0.20 8.83 -3.8 0.10 8.65 -1.7 0.02
32 7.98 8.06 -1.0 0.01 8.83 -10.6 0.72 8.65 -8.3 0.44
33 10.83 9.16 15.4 2.79 9.75 10.0 1.16 9.63 11.1 1.44
34 10.62 9.06 14.7 2.43 9.62 9.4 1.01 9.50 10.6 1.26
35 10.72 9.06 15.4 2.74 9.62 10.3 1.21 9.50 11.4 1.48
36 12.86 12.17 5.4 0.48 12.72 1.1 0.02 12.61 2.0 0.06
37 12.86 12.28 4.5 0.33 12.87 -0.1 0.00 12.76 0.7 0.01
38 12.24 12.21 0.2 0.00 12.78 -4.4 0.29 12.67 -3.5 0.19

SSE 75.15 44.30 49.53

 

80

 

From this result, one can conclude that the contact between the ice and
the insulating package as well as geometry difference is a very crucial factor for

simulation of the heat transfer phenomena of insulating packages.

5.5 Construction of simplified model for practical use
The computer model considering contact between the ice and the
insulating package successfully predicted the heat transfer of the insulating
packages. From this model, a simplified equation was derived for practical use.
As discussed in theory part, system R-value of the entire insulating
package was described in equation (64). ln this equation, there are four heat
transfer coefficients, k13, k24, how, hi", and all of these parameters are
temperature and geometry dependent. In packaging applications, however, these
parameters are distributed in a narrow range. Equation (89) was rearranged to
yield equation (90) for the thermal conductivities, kg and k24.
k1 1:24 _ 1

_ = _ (90)
Ax Ax Rwall

 

According to the pre-calculation (data not shown), the radiation part of
both equation and the pure convection heat transfer coefficients were also very
stable. Considering those numbers, equations (55) and (60) were modified to

represent how and hi", respectively. The modified equations are as follows:

hout = 0'312+0*912£0ut (91)

81

 

A5 0.807

 

 

 

 

 

 

h- = 0.624 + 92
A5 Sin gice
then,
Ain Ain
= + (93)
system hour Aout 1 A2
1 + Rwall Rwall
hinA3 A1143

The simulation resultscalculated using equations (89), (90), (91), (92) and

(93) were compared with the experimental data. The simulation results were also

compared with the conventional system R-value calculation proposed by Burgess

(1999), which were shown in equation (3). Table 13 summarizes the comparison

of the experimental data and the two simulation results. The simulation with the

above equations predicts the experimental data within 25%, which was

dramatically improved numbers comparing the simulation results with the

conventional model in equation (3)(< 69% error).

82

 

Table 13. Comparison of the simulation results with equation (93) with
experimental data.

 

 

 

 

Parameters Eq (93) Eq (3)
hout hin exp Rsymm % Error Rsymm % Error

1 1.02 0.51 6.06 6.50 -7.2 8.40 -38.6

2 1.02 0.51 6.10 6.50 -6.6 8.40 -37.8

3 1.02 0.51 6.03 6.50 -7.7 8.40 -39.2

4 1.00 0.22 11.32 9.04 20.1 8.40 25.8

5 1.00 0.47 8.31 7.01 15.7 8.40 -1.1

6 1.00 0.47 7.36 7.01 4.7 8.40 -14.2

7 1.00 0.47 6.92 7.01 -1.2 8.40 -21.3

8 0.97 0.27 11.46 8.56 25.4 8.40 26.7

9 0.95 0.50 9.44 7.16 24.1 8.40 11.0
10 0.95 0.50 8.03 7.16 10.8 8.40 -4.6
11 0.95 0.50 7.36 7.16 2.6 8.40 -14.1
12 0.90 0.14 16.07 11.99 25.4 8.40 47.7
13 0.94 0.60 7.40 6.35 14.2 8.40 -13.5
14 0.94 0.60 4.99 6.35 -27.4 8.40 -68.4
15 0.94 0.60 5.13 6.35 -23.8 8.40 -63.7
16 0.83 0.18 11.40 10.04 11.9 8.40 26.3
17 0.95 0.64 7.66 6.68 12.7 8.40 -9.7
18 0.95 0.64 7.03 6.68 5.0 8.40 -19.4
19 0.95 0.64 6.13 6.68 -9.0 8.40 -37.0
20 0.89 0.16 15.87 11.06 30.3 8.40 47.1
21 1.03 0.17 5.23 5.69 -8.7 5.11 2.3
22 1.03 0.17 5.48 5.80 -5.9 5.11 6.8
23 1.03 0.17 5.60 5.86 -4.7 5.11 8.8
24 0.95 0.16 5.45 6.14 -12.5 5.11 6.3
25 0.95 0.15 5.73 6.26 -9.3 5.11 10.8
26 0.95 0.15 5.85 6.33 -8.3 5.11 12.6
27 0.99 0.25 10.16 8.80 13.4 10.35 -1.9
28 0.99 0.25 10.12 8.80 13.0 10.35 -2.3
29 0.99 0.27 9.28 8.53 8.1 10.35 -11.5
30 0.95 0.20 8.25 8.53 -3.4 8.40 -1.8
31 0.95 0.20 8.51 8.61 -1.2 8.40 1.3
32 0.95 0.20 7.98 8.61 -7.8 8.40 -5.2
33 0.99 0.23 10.83 9.46 12.6 9.62 11.2
34 0.99 0.24 10.62 9.33 12.2 9.62 9.5
35 0.99 0.24 10.72 9.33 13.0 9.62 10.3
36 0.99 0.26 12.86 11.50 10.6 9.01 30.0
37 0.99 0.25 12.86 11.66 9.4 9.01 30.0
38 0.99 0.25 12.24 11.56 5.5 9.01 26.4

 

 

CHAPTER 6

CONCLUSIONS

6.1 Practical use of model equations

A mathematical model which predicts the performance of an insulating
package was developed. This model includes various factors affecting the quality
of the insulating package. Multi-layered walls, the geometry of the package and
the degree of contact between the product and insulating package are among
these factors. The system R-value for the entire insulating package was predicted
using this model. The diagram of the insulating package is shown in Figure 17. In
order for the model predictions to apply. the actual package and product must be
able to be described as in the figure. The insulating package has outside surface
area A1 + A2, inside surface area A3 + A2, and average wall thickness Ax. The
product has surface area A2 + A4 with A2 being the contact area between the

product and package.

Outside Area:
Aout=A1 +A2

Inside Area:
Aln=A3+A2

Product Area:
Aoroduct= A4+A2

 

Figure 17. Diagram of typical insulating package

84

Equations (89), (91 ), (92), (93) are used to find the system R-value for the

entire package according to the procedure described below.

Step1: Calculate the surface areas A1, A2, A3 and A4 in sq. ft. Also

calculate the average package wall thickness Ax.

Step 2: Find the outside surface emissivity 60“,, inside surface emissivity

a,” and product surface emissivity s pmdua. Use the table below.

(’18;

 

 

Material Emissivity (e )
Aluminum 0.05
Copper 0.03
EPS foam 0.95
Ice 0.96
paper 0.95
wood 0.65-0.95

 

Step 3: Find the surface heat transfer coefficient using

A4 0.807
+
A3 '1' A4 A3 1 + 1

A4 8in 8product

113 = 0.624

 

—1

 

Step 4: Calculate the wall R-value using

N N 3 1.229d ' 3.
Rwall = 212,661,- + 5.56tha1+ 5.31thaZ + Z a" 2

i=1 1.210.233 + dair3,-

 

= the number of solid layers

= the thermal resistance per inch for the solid material of the ith layer

(See Table below)

= the thickness of solid material of i th layer

85

 

 

 

 

 

 

 

thal = the total thickness of both-walI-foiled air space in inches
tha2 = the total thickness of one-wall-foiled air space in inches
N3 = the number of the airspaces without foiled layers
dair3i = the thickness of the i th non-foiled air space in inches
, Thermal resistance per inch at 70°F
Matenals
R,- (hrftz'oF/BTUin)
Air 5.56
Corrugated board 2.38
EPS foam 3.78
Polyurethane 4.62
VIP 10-30
Step 5: Calculate the system _R—value using
A- A-
R = m + m
SYSfem hout/10m 1 + A2
1 + RWG” Rwall
Step 6: Calculate the heat transfer rate using
Q = Ail! X AT
Rsystem
Q = the heat transfer rate in BTU/h
A,-,, = the inside surface area of package in fiz
AT = the temperature difference between outside air and the product in °F

Rsystem = the thermal resistance of entire package in fi2-°F - hr/ BTU

86

 

divl

hee

6.2

00

fit

The rate at which ice melts (lb/hr) is equal to the heat transfer rate Q
divided by latent heat of the ice (or gel pack). When regular ice is used, the latent

heat is 144 BTU/lb.

6.2 Example

- These equations can be used for practical situations. As an example,
consider an ordinary EPS cooler (18 in x 16 in x 12 in 0.0., 1.5 in thick) loosely
fitted in a corrugated C-flute box with foil outside. The average distance between
EPS cooler and corrugated box is 0.1 in. If 20 lbs of frozen fish (80% water) at
32 °F with dimensions of 12in x 10in x 8 in is stored in this cooler, and the cooler

is stored in an outside temperature of 100 °F, how long will the fish remain

frozen?

$2.11 Calculate the surface areas:
The insulating material is composed of 1.5 in thick EPS foam, 0.156 in thick
corrugated board and the 0.1 in thick air space between EPS cooler and
corrugated board. In this case, the outside dimensions of the entire package are

18.512 in x 16.512 in x 12.512 in. The average thickness of the insulating

material is

Ax =1.5+O.156+0.1=1.756 in

and

A0,, = (2 (18.512 x 16.512) + 2 (18.512 x 16.512) +2 (18.512 x 16.512))/144 =10.33f:2

A," = (2 (15.0 x 12.0) + 2 (12.0 x 9.0) +2 (9.0 x 15.0))/144 = 6.20 ftz

87

 

Amduc, = (2 (12.0 x 10.0) + 2 (10.0 x 8.0) +2 (8.0 x 12.0))/144 = 4.11 ft2
A2 = (12.0 x 10.0) /144 = 0.83fi2

A1= A0,“ - A2 = 10.33 — 0.83 = 9.50fz2

A3 =A,-,, - A2 = 6.20 — 0.83 = 5.37 ft2

A4 = Apmduc, - A2 = 4.11— 0.83 = 3.28 ftz

Step 2: From the table in step 2 of the procedure:

so“, = 0.05, s," = 0.95, 6,68: 0.95

Step 3: The surface heat transfer coefficients are
ho“, = 0.312 + 0.912 x 0.05 = 0.357 BTU/hft2‘°F
h,-,, = 0.624x 3.28/ (3.28+5.37) + 0.807 / (5.37/3.28 x 1/0.95 +1/O.95 —1)

= 0.691 BTU/hft2'°F

Step 4: From the equation in step 4 of the procedure

N = 2 (there are two solid layers: EPS foam and corrugated board)
R, = 2.38 for corrugated board (table in step 4)

R2 = 3.78 for EPS foam (table in step 4)

thsl = 0.156 in for the thickness of corrugated board

thsz = 1.5 in for the thickness of the EPS

thal = 0.0 (there is no both-waII-foiled air space)

thaz = 0.0 (there is no one-walI-foiled air space)

N3 = 1 there is one airspace without foil

88

(layer between EPS cooler and corrugated board)
dair3i = 0.1 in for the thickness of the non-foiled air space
Then,

Rm,” = 2.38x0.156 + 3.78xl.5 + 1.229x 0.1 /(0.233 + 0.1) = 6.41 hft2‘°F/BTU

Step 5: The system R—value is

 

 

6.20 6.20
+

Rsys’e’" = 0.357-10.33 1 + 0.83
1 6.41 6.41

+
0.691-5.37 J9.50-5.37

= 7.96 hft2'°F/BTU

 

 

Step 6: The heat transfer rate Q is

_ A,-,, ~AT _ 6.20.(100—32)
Rsystem 7'96

= 52.97 BTU/h

 

 

Q

The frozen fish is essentially a 16 lb block of ice (80% of 20 lb). The heat
required to melt 16 lbs of ice is 144 BTU/lb x 16 lbs = 2304 BTU. The time for the

frozen fish to completely thaw in the package can be calculated as

I: —1———h——- 2304BTU = 43.6h

52.9 BTU

Since the system R-value is a property of the insulating package, it is not
restricted to steady state ice melting. In an unsteady state, the system R-value
can be used to find the transient response. For example, it can be used to find

the temperature of the fish after it melts. It will gradually warm up to the outside

89

 

temperature. The heat transfer rate in unsteady state is as follows (Holman,

1986)

1 dT

00) = 41.081040» = WCp-d—t (94)

 

system

where Q(t) is the heat transfer rate as a function of time in BTU/h,t is the time in
h, A," is the inside surface area of the package in frz, T0,,(t) is the environmental
temperature as a function of time in °F, T is the temperature of the product as a
function of time, W is the weight of the product in lb, and Cp is the heat capacity
of the product in BTU/lb°F.

Once system R-value is known, one can easily integrate equation (94)
using the temperature profile of the environment and thermal properties of the
product to obtain the temperature profile of the product. For example, if the

environmental temperature, T00, is constant, integration of equation (94) results in

Air:
WC p Rsystem

 

T(t)=Too_(Too_Ti )6Xp(- t) (95)

where T, is the initial temperature of the product.

In the example discussed earlier, the fish in the package completely thaws
after 43.6 hours. Immediately after all the fish thaws, T,- is 32 °F, Wand Cp for the
product are 20.0 lb and 0.8 BTU/lb'°F, respectively. The temperature profile of the

product can be calculated using equation (95).

6.20
20 - 0.8 - 7.96

 

7(7) = 100 — (100 — 32) exp(— t) = 100 — 68 exp(—0.0486t) (96)

Equation (96) is profiled in Figure 18.

90

 

 

100 *1

m
C
1

O)
O
1

40*

 

Temperature (°F)

20-1

 

 

I I I ' l

0 20 40 60 80 100
Time (h)

 

Figure 18. Temperature profile after the product melts

The system R-value developed in this study can also be used for the
simulation of the performance of the insulating package in a temperature variant
environment such as the thermal cycling environment used in ISTA test
standards (See Table 3). Use of this model will save time, effort and cost by

minimizing the construction of actual packages and subsequent performance

tests.

91

6.3 Factors affecting on package construction
Several factors affect the thermal resistance of the insulating package.
Using the example in the previous section, these factors can be investigated for

the design of the insulating package.

6.3.1 The effect of wall thickness on the performance of the insulating F:
package
The effect of wall thickness on the performance of the insulating package

can be determined for the example discussed earlier. If the wall thickness is very

 

small, like a plastic bag, and the other factors remain the same, the effect of the
wall thickness can be estimated. Suppose that the wall is composed of EPS foam
and the average thickness of the wall is very small, 0.01 in. The inside
dimensions of the package remain the same to maintain the same interior
volume. If other conditions remain the same as the earlier example, the heat
transfer rate can be calculated by the following procedure:

M: Calculate the surface areas:
The average thickness of the insulating material is
Ax = 0.01 in
A0,, = (2 (15.01 x 12.01) + 2 (12.01 x 9.01) +2 (12.01 x 9.01))/144 = 6.22ft2
A,” = 6.20ft2, Apmduc, = 411172.212 = 0.83f12
A1: A0,“ - A2 = 6.22 -— 0.83 = 5.39fz2

A3 = 5.37 ftz, A4 = 3.28 ftz

92

Step 2: From the table in step 2 of the procedure:

80,“ = 0.05, Si" = 0.95, ewe: 0.95

Step 3: The surface heat transfer coefficients remain the same:
ho“, = 0.357 BIU/hfi2'°F

h," = 0.691 BTU/hftz'°F

Step 4: From the equation in step 4 of the procedure

N = 1 (there are two solid layers: EPS foam and corrugated board)
R, = 3.78 for EPS foam (See Table 10)

thsl = 0.01 in for the thickness of corrugated board

Then,

Rm,” = 3.78 x 0.01 = 0.038 hft2'°F/BTU

Step 5: The system R-value is

6.20 6.20

R = +
”3"?“ 0.357 - 6.22 1 + 0.83
1 0.038 0,033

+
0.691-5.37 45395.37

 

 

= 3.03 hftz'°F/BTU

 

 

Step 6: The heat transfer rate Q is

_ A," ~AT _ 6.20-(100—32)
Rsymm 3.03

 

Q = 139.4 BTU/h

93

Without wall thickness, the insulating package loses 63% of the system R-
value compared to the original example and the heat transfer rate is almost 3

times as high.

6.3.2 The effect of aluminum foil on the performance of the insulating
package

The example previously discussed is an interesting situation. In that
example, the insulating package, which has lost most of its thermal resistance of
the wall, still retains more than 37% of its system R-value. This is mainly due to
the contribution of the aluminum foil. If the effect of aluminum foil is eliminated,
which means that the ho,“ value increases up to 1.18 BIU/hft2'°F, then the
system R-value decreases 64% to 1.09 hft2'°F/BTU and heat transfer rate
dramatically increases to 388 BTU/h.

The effect of the aluminum foil is less significant when the package has
substantial thickness. When the foil effect is eliminated from the original example
with other variables remaining the same, the system R-value of the insulating
package decreases 15% to 6.80 hft2'°F/BTU. One can imply from these results
that the effect of foil is more significant when the thickness of the wall is small.
These results coincide with the prediction model by Burgess (Burgess, 1999).

One can apply aluminum foil on the outside surface, on the inside surface
or on the product. Using the previous two examples with a different thickness of

the wall, the location effect of the aluminum foil can be investigated. Alternating

94

the emissivities on three different locations, one can get the the system R-values

shown in Table 14.

Table 14. System R-values with various wall thicknesses and foil locations

 

System R-value according to the foil location

 

 

Inside of the Outside of the
Product package package
Wall thickness of 1.75 in 8.61 8.77 797

(Original example)

Wall thickness of 0.01 in
(simulation by reducing 1.10 1.10 3.03
wall thickness)

From this result, one can have a clue that the system R-value is affected
by the location of the foil. If the insulating material is thin, the effect of foil is
maximized when the foil is applied on the outer surface of the package. When
the insulating material has substantial wall thickness, foil should be located at the
inner surface of the insulating package.

In addition, proper location of aluminum foil can reduce the volume of the
insulation package significantly. With the application of the foil on the inner
surface of the package, the overall heat transfer rate of insulating package of the
original example can be reduced to 48.1 BTU/h. This value is equivalent to the
heat transfer rate of a 2.5 in thick non-foiled EPS container with the same inner
and product dimensions. In this case, one can reduce the volume of the entire

package by more than 30%.

95

6.3.3 The effect of the surface area on the performance of the lnsulatlng
package

The effect of surface area on the performance of the insulating package
can be determined for the example discussed earlier. Suppose that all
dimensions of the original example increase twice with the geometric ratio of the
entire package remaining the same. In this case, even though the system R-
value increases 60% to yield 12.85 hft"°F/BTU, the overall heat transfer rate still
increases to 131.5 BTU/h, which is approximately 2.5 times higher than original
value. This implies that the role of the surface area is crucial for the performance
of the insulating package even though the geometric ratios of all dimensions
remain the same. When all dimensions of the insulating package increase twice,
the surface areas increase four times. This faster increase of the surface can be
an explanation of the increase of the overall heat transfer rate in spite of the
increase of system R-value.

When all dimensions of the original example are reduced to 1/2 of their
original value, the system R-value decreases to 5.46 h'ftz'°F/BTU, which is 70% of
its original value. In this case, the overall heat transfer rate also decreases to
19.3 BTU/h, which is approximately 36% of its original value. The same surface
area and dimension relationship can explain this result.

This result can be applied to the fish thawing time example. When the
dimension of the product increases 2 times, the volume of the product increases
8 times. Since the density remains the same, the frozen fish in the original

example should be treated as 128 lbs of ice and the heat required to melt the fish

96

is 144 BTU/lb x 128 lbs = 18432 BTU. The time for the frozen fish to completely

thaw in the package can be calculated as

t = -1——L-1843ZBTU =140h
131.5 BTU

This is more than three times of the original thawing time. In similar way,
the thawing time of the frozen fish for 1/2 reduced dimension problem is
calculated as 15 hours, which is 34% of the original thawing time.

The effect of the surface area of the product can be also determined for
the original example. Suppose that the dimensions of the product in original
example are reduced to 1/2 of their original value with other variables remaining
the same. In this case, the system R-value of the insulating package increases
35% to yield 10.81 hft2'°F/BTU and the overall heat transfer rate decreases to
39.1 BTU/h, which is approximately 73% of its original value.

The examples discussed in this chapter are very specific cases of the
insulating package. The actual insulating package might involve more
complicated factors. However, it is obvious that the model developed in this study
can be a very useful tool to analyze and predict these factors. As discussed in
this chapter, the model can be applied to various uses of the insulating package:
ice-melt test simulation, transient thermal cycling simulation, analysis of the
factors for the construction of the package, etc. Using this model, one can get
more insightful and systematic approach to the design and analysis of the

insulating package.

97

APPENDICES

98

APPENDIX A. Basic programs for simulation of the insulating package

(AREA.BAS)

5 REM Effective Area of the heat conduction

10 OPEN "O",#1,"Area.res"

20 set=7

100 For I=1 to set

200 Read L, w, D, t

300 Ao=2*(L*w + w*o + D*L)

400 Ai=2*((L-Z*t)*(w-2*t)+(w-2*t)*(D-2*t)+(D-2*t)*(L-2*t))
500 5:5 r(((L-w)A2+(L-D)A2+(w-D)A2)/2)

600 X1= L+w+D+S)/6: X2=(L+w+D-S)/6

700 If L=W and L=D Then I=t/(24*X1*(X1-t)):GOtO 850
800 I=Log(X2*(x1-t)/(x1*(X2-t)))/(24*(x1-x2))
850 Pr1nt #1, 1/I, Sqr(Ao*Ai)/t, (Ao+Ai)/(2*t)
900 Next I

1000 Data 18, 13, 9, 2

1010 Data 100, 50, S, 1

1020 Data 24, 24, 24, 3

1030 Data 36, 12, 36, 4

1040 Data 20, 15, 10, 0.5

1050 Data 30, 20, 25, 5

1060 Data 12, 80, 15, 2

1990 CLOSE #1

2000 END

(NGR.BAS)

5 REM NGR

10 OPEN "O",#1,"NGR.res"

20 L=3

100 For Tave=40 to 100 step 30

200 For dT=1 to 10

300 For d=0.005 to 0.1 step 0.005

400 NGR=dA3*dT*(.1794+4.1794*EXP(—.0095034*TAVE))*1E6
500 Print #1, Tave, dT, d, NGR

600 Next d

700 Next dT

800 Next Tave

1990 CLOSE #1

2000 END

99

(NUSSELT.BAS)

50
100
200
1000
2000
3000
3100

REM Nusse1t And h’s
OPEN "O",#1,"Nusse1t.res"
L:
For Tave=30 to 120 step 10
For dT=1 to 10
For d=0.005 to 0.1 step 0.005
NPR=.07206-1.8453E-04*TAVE+3.5979E-07*TAVEA2—

2.9967E-10*TAVEA3

3200
4000
4100
4200
4300
4350
4400
5000
6000
7000
8000
19980
19990

29000
29010
29050
29100
29110
29150

29200
29210
29230
29270
29290

30000
30100
31050
31100
31500

32000
32100
32200
32300
32500
32700
32900

40000
40100
40200

KAIR=.0131+2.6136E-05*TAVE-6.1955E-09*TAVEA2
NGR=dA3*dT*(.1794+4.1794*EXPC-.0095034*TAVE))*1E6
NGRF=LA3*dT*(.l794+4.1794*EXP(-.0095034*TAVE))*1E6
Gosub 29000
Gosub 30000
If NGR<8000 then hb=kair/d: goto 5000
Gosub 40000
Print #1, Tave, dT, d, NGR, ha, hb, hf
Next d
Next dT
Next Tave
CLOSE #1
END

Rem *******Kreith***********************

IF NGR>=2000 THEN 29100
HA=kair/d: goto 29150
60508 29200
HA=KAIR/d*A*NGRAM
RETURN

REM convection coefficient of air for verticai plane

IF NGR>=20000 THEN GOTO 29270
A=.18*(L/d)A(-1/9):M=1/4:GOTO 29290
A=.065*(L/d)A(—1/9):M=1/3

RETURN

REM *************Free Convect-i on*****************

NPRNGR=NPR*NGRF

GOSUB 32000

HF= KAIR/L*A*NPRNGRAM
RETURN

REM convection coefficient of air for verticai plane
IF NPRNGR<1E+04 THEN GOTO 32500

IF NPRNGR>1E+09 THEN GOTO 32700
A=.59:M=1/4:GOTO 32900
A=1.36:M=1/5:GOTO 32900
A=.13:M=1/3

RETURN

Rem *************Dr . Burgess'k******‘k‘k************
Ra=NPR*NGR
If 2<L/d and L/d<10 and Npr<10 and Ra<1E10 then gosub

45000:goto 44000

40250 Ra1=Ra*Npr/(0.2+Npr)

40300 If 1<L/d and L/d<=2 and 10E-3<Npr and Npr<1E5 and Ra1>1E3
then gosub 46000zgoto 44000

100

40400 If 10<L/d and L/d<40 and 1<Npr and Npr<2E4 and 1E4<Ra and
Ra<1E7 then gosub 47000:goto 44000

40500 If 1<L/d and L/d<40 and 1<Npr and Npr<20 and 1E6<Ra and
Ra<lE9 then gosub 48000:goto 44000

40600 nb=-1:goto 44500

44000 hb=Nu*Kair/d

44500 Return

45000 Nu=0.22*(L/d)A(-1/4)*(Npr/(Npr+0.2)*Ra)A0.28
45100 Return

46000 Nu=0.18*Ra1A0.29
46100 Return

47000 Nu=0.42*RaAO.25*NprA0.012/(L/d)A0.3
47100 Return

48000 Nu=0.046*RaAO.33
47100 Return

101

(WALLBG7.BAS)

5 REM "Dr. Burgess Version #2

100 OPEN "O",#1,"wa11bg7. res"

200 READ CRITERIA, MAXITER

300 DATA 0.005, 20

350 nset=12: nnum=1000: print #1, "990 data ";nset

375 BOLZ=1.714E- 09 : inicon=0. 015: L=3: Nw=20: Unknown=2*nw+1:
nw2=2*nw

380 DIM T(2*nw), TAiran), anw), dairan), KanL H(nw+1),
BUFFERCUnknownL E(2*NwL hrad(Nw+1L hcon(nw+1L property(nw)
390 DIM HACZL C(unknown, Unknown+1L TEMP(unknown),
dairnormaIan)

400 For zzzz=1 to nset

450 Read Nw: Unknown=2*nw+1: nw2=2*nw

640 REM reading dimensions and properties of the walIs
650 FOR I= -1 TO nw

700 READ ww, xx, propertyCiL2e(2*i-1L e(2*i)

800 d(i)=ww/12: dair(i)=xx/12

900 NEXT I

902 data 1, 1.5,5, 2, 0.95.0.95

904 data 1, 3, 5, 2, 0.95, 0.95

906 data 2, 1.5, 0.1, 2, 0.95, 0.050, 1.5, 5, 2, 0.050, 0.950
908 data 2, 1.5, 0.3, 2, 0.95, 0.05, 1.5, 5, 2, 0.050, 0.950
910 data 3, 1.5, 0.2, 2, 0.95, 0.050, 1.5, 0.5, 2, 0.050,
0.050, 1.5, 5, 2,0.950,0.95

912 data 3, 0.156, 0.1, 2, 0.95, 0.05, 0.156, 0.4, 2, 0.05,
0.05, 0.156, 5,2,0.05,0.95

914 data 5, 0.156, 0.2, 2, 0.95, 0.05, 0.156, 0.1, 2, 0.05,

0. 05, 0. 156, 0.3,2,0.05,0. 05, 0.156, 0. 7, 2, 0.05, 0.05,

0. 156, 5, 2, 0. 05,0. 95

916 data 5, 0.156, 0.3, 2, 0. 9S, 0. 05, 0.156, 0.2, 2, 0.05,
0.05, 0.156, 0. 3,2,0. 05,0. 0L 0. 156, 0. 6, 2, 0.05, 0.05,
0.156,5,2,0.05,0. 95

918 data 5, 0.156, 0.4, 2, 0.95, 0.05, 0.156, 0.3, 2, 0.05,
0.05, 0.15L 0. 2,2,0. 05,0. 05, 0.156, 0.5, 2, 0.05, 0.05,
0.156,5,2,0.05,0. 95

920 data 5, 0.156, 0.5, 2, 0.95, 0.05, 0.156, 0.3, 2, 0.05,
0.05, 0.15L 0. 2,2,0. 05.0.05, 0.156, 0.4, 2, 0.05, 0.05,
0.156,5,2,0.05,0. 5

922 data 5, L 156, 0.6, 2, 0.95, 0.05, 0.156, 0.2, 2, 0.05,
0.05, 0.156, 0.1,2,0.05,0.05, 0.156, 0.3, 2, 0.05, 0.05,

0.156, 5, 2, 0. 05,0. 5

924 data 5, 0.156, 0.7, 2, 0. 95, 0. 05, 0. 156, 0. 1, 2, 0.05,
0.05, 0.156, L 1, 2,0. 05,0. 05, 0. 156, 0. 2, 2, 0. 05, 0. 05,

0.156, 5, L 0. 05, 0. 95

1000 REM “Heading0 known Temperature*****************
1100 Tin=40 Tout-
1150 T(0)=Tout: Tair(0)=Tout: Tair(nw)=Tin

1200 REM *********initiai guess of each temperature*********
1250 DELTA=(Tout- Tin)/unknown
1260 For 1: 1 to nw

1270 T(2*1-1)=T(2*(1-1))-DELTA
1280 T(2*1)=T(2*1-1)-DELTA
1300 next 1

102

1400

FOR KK=1 T0 MAXITER

1450 GOSUB 2000

1650 FOR 3= 1 TO unknown

1700 TBUFFER = TBUFFER + ABSCCTempCJ)-BUFFER(J))/Temp(J))

1750 NEXT J

1800 FOR J=1 TO unknown: BUFFER(3)=Temp(J): NEXT J

1810 For j=1 to nw2: T(j)=Temp(j):Next j

1850 IF TBUFFER < CRITERIA THEN 1930

1900 TBUFFER=0

1910 NEXT KK

1915 REM did NOT converge

1920 PRINT "dataset ";zzzz;" has not been converged":

1930 REM converged

1965 REM ***************Ca1cu1ation of h and k*******

1970 gosub 50000

1972 next 2222

1998 CLOSE #1

1999 END

2000 REM ****************main subroutine*********************

2050 REM ***************Loo for materia] properties*********

2100 Ta1=Tout : Ta2=T(1§: E=E(1): jj=1 : gosub 30100

2170 If nw=1 then 2550

2200 For i=1 to nw-1

2300 Ta1=T(2*i) : Ta2=T(2*i+1): dair1=dair(i) :El

=E(2*i) :E2=E(2*i+1) :jj=i+1: gosub 29020

2500 Next i

2550 Ta1=Tout : Ta2=T(2*nw): E=E(2*nw):jj=nw+1: gosub 30100

2600 For i=1 to nw _ .

2620 Ta1=T(2*i-1) :Ta2=T(2*i):mater1a1=property(1):gosub
35000: k(i)=k

2640 Next i

3000 REM ***********Loop for construction of Matrix********

3100 Gosub 10000

4000 REM *****impose known variab1es and reduce matrix*******
5000 REM ***************so]ving matrix********************
5100 60508 20000

6000 REM ***************************************************

8400 QTEMP=C(unknown,unknown+1)

8500 FOR I=1 T0 unknown: Temp(1)=C(I,unknown+1): NEXT I

8600 FOR I=1 TO unknown: FOR J=1 TO unknown+1: C(I,J)=0: NEXT
3: NEXT I

8700 FOR I=1 TO nw+1: H(I)=0: NEXT I

9000 RETURN

10000 REM ***********Constructing
Matrix**************************

103

11000 c(1,1)=h(1):c(1,unknown)=1:c(1,unknown+1)=h(1)*Tout

11100 c(2,1)=-h(1)-k(1)/d(1): c(2,2)=k(1)/d(1): c(2,unknown+1)=-
h(1)*Tout

11300 If nw=1 then 12000

11400 for i= 1 to nw-1

11500 c(2*i+1, 2*i-1)=k(i)/d(i): c(2*i+1, 2*i)=-h(i+1)-

k(i)/d(i): c(2*i+1, 2*i+1)= h(i+1)

11600 c(2*(i+1L 2*i)= h(i+1): c(2*(i+1L 2*i+1)=-h(i+1)—
k(i+1)/d(i+1): c(2*(i+1L 2*(i+1))= k(i+1)/d(i+1)
11800 next i

12000 c(2*nw+1,2*nw-1)=k(nw)/d(nw):c(2*nw+1,2*nw)=-h(nw+1)-
k(nw)/d(nw):c(2*nw+1,unknown+1)=-h(nw+1)*Tin

14000 Return

20000 REM ******subroutine for so1ving Matrix***********

20010 REM: SOLVES SYSTEM OF EQNS BY GAUSSIAN ELIMINATION

20060 FOR I= 1 TO unknown : PVT=I 'e1imination using row I ...

20070 FOR K=I+1 TO unknown

20080 IF ABSCCCK, I))>ABS(C(PVT, 1)) THEN PVT=K

20090 NEXT K 'row "PVT" has 1argest e1ement in Ith co1umn

20100 FOR J: I TO unknown+1'interchan e row I and pivot row

20110 CHG=C(I, J) : C(I, J)=C(PVT, J? : C(PVT, J)=CHG

20120 NEXT J 'now row I has 1argest e1ement in Ith co1umn

20130 FOR K=1 TO unknown : IF K=I THEN 20160

20140 F=C(K, I)/C(I, I)

20150 FOR J=I+1 TO unknown+1 : C(K,J)=C(K,J)-F*C(I,J) : NEXT J

20160 NEXT K 'Ith co1umn now zeroed out (except for row I)

20170 NEXT I

20180 FOR K=1 TO unknown:C(K,unknown+1)=C(K,unknown+1)/C(K,K)
:NEXT K

20200 RETURN

29020 REM ****convection coefficient for enc1osed space******

29030 TAVE=(TA1+TA2)/2

29040 NPR=. 07206- 1. 8453E-04*TAVE+3.5979E-07*TAVEA2-2.9967E—

10*TAVEA3

29045 KAIR=.0131+2.6136E-05*TAVE-6.1955E-09*TAVEA2

29080 NGR=dair1A3*ABS(TA2-TA1)*(.1794+4.1794*EXP(-

.0095034*TAVE))*1E6

29085 IF NGR>= 2000 THEN 29100

29087 HconC g)=—kair/dair1:goto 29120

29100 GOSUB 29

29110 Hcon(j8)= —KAIR/dair1*A*NGRAM

29120 If E1- then E11=1E6: goto 29140

29130 Ell=1/ E1

29140 If E2=0 then E22= 1E6: goto 29155

29150 E22=1/E2

29155 Ta1=Ta1+460z Ta2=Ta2+460

29160 hradea)=bo1;*(TalA3fTa1A2*Ta2+Ta1*Ta2A2+Ta2A3)/(E11+E22—1)

29170 h(j3)= con(33)+hrad(JJ)

29190 RETURN

29200

REM convection coefficient of air for vertica1 p1ane

104

29210
29230
29270
29290

30100
30500
30600

IF NGR>=20000 THEN GOTO 29270
A=.18*(L/dair1)A(-1/9):M=1/4:GOTO 29290
A=.06S*(L/dair1)A(-1/9):M=1/3

RETURN

REM conductivity and FREE convection coefficient for air
TAVE=(TA1+TA2)/2
NPR=.07206-1.8453E-04*TAVE+3.5979E-07*TAVEA2-2.9967E-

10*TAVEA3

30650
30800

KAIR=.0131+2.6136E-05*TAVE-6.1955E-09*TAVEA2
NGR=LA3*ABS(TA2-TA1)*(.1794+4.1794*EXPC-

.0095034*TAVE))*1E6

30850
31050
31100
31350
31400
31500

32000
32100
32200
32300
32500
32700
32900

35000
35500
35550
35560
35570
35580
35590
35640
35650
35660
35670
35680
35690
36000

NPRNGR=NPR*NGR

GOSUB 32000

HC= KAIR/L*A*NPRNGRAM
hrad(jg)=E*bo1z*((TA2+460)A4-(TA1+460)A4)/(TA2-TA1)
293:: MW

REM convection coefficient of air for vertica1 p1ane
IF NPRNGR<1E+04 THEN GOTO 32500

IF NPRNGR>1E+09 THEN GOTO 32700

A=.59:M=1/4:GOTO 32900

A=1.36:M=1/5:GOTO 32900

A=.13:M=1/3
RETURN

REM ****conduction coefficients for each materia1s*********
TAVE=(TA1+TA2)/2
If materia1 = 0 then 35650
If materia1 = 1 then 35660
If materia1 = 2 then 35670
If materia1 = 3 then 35680
If materia1 = 4 then 35690
k=0.346: goto 36000 'water
k=.0131+2.6136E-05*TAVE-6.1955E-09*TAVEA2:goto 36900’ air
k=0.035: goto 36000 'corrugated board
k=0.022: goto 36000 'EPS foam
=0.018: goto 36000 'Po1yurathane
=0.010: oto 36000 'VIP
k=k*(.013 +2.6136E-05*TAVE-6.1955E-

09*TAVEA2)/(.0131+2.6136E-05*72-6.1955E-O9*7242)

36900

50000

52100 Ta1=Tout

52170
52200
52300

52500
52550
52600
52620

RETURN

Rem ******Ca1cu1ation of Fina1 va1ues of h's and k's******
: Ta2=T(1): E=E(1): jj=1 : gosub 30100

If nw=1 then 52550

For i=1 to nw-1

Ta1=T(2*i) : Ta2=T(2*i+1): dair1=dair(i)
:El =E(2*i) :E2=E(2*i+1) :jj=i+1: gosub 29020

Next i
: Ta2=T(2*nw):E=E(2*nw):jj=nw+1:gosub 30100

Ta1=Tout
For i=1 to nw . .
Ta1=T(2*i-1):Ta2=T(2*i):mater1a1=property(1)

105

:gosub 35000: k(i)=k
52640 Next i

52750 =QTEMP
55600 TR=0
56000 For i=1 to nw :Totwa11=Totwa11+d(i)*12

:TR=TR+d(i)/k(i) :next i
56050 For i=2 to nw :TR=TR + 1/h(i) :Next i
56100 For i=1 to nw-1

56110 If e(2*i)<0.1 and e(2*i+1)<0.1 then
Totair1=Totair1+dair(i)*12:goto 56150

56120 If e(2*i)<0.1 or e(2*i+1)<0.1 then
Totair2=Totair2+dair(i)*12: oto 56150

56130 nnorma1=nnorma +1: dairnorma1(nnorma1)=dair(i)*12

56150 Next i

56155 If nnorma1 =0 then print #1, nnum;" data ";TR;" , ";
totwa11;" , "; totair1;" , "; totair2;" , ";nnorma1: goto 56195

56160 For i=1 to nnorma1: Totair3=Totair3+dairnorma (1): next

1

56170 print #1, nnum;" data ";TR;" , "; totwa11;" , ";
totair1;" , "; totair2;" , ";nnorma1:" , ";

56175 If nnorma1=1 then 56185

56180 _ For i=1 to nnorma1-1: print #1, dairnorma1(i);" , ";:
next 1

56185 print #1, dairnorma1(nnorma1)

56195 nnum=nnum+10

56200 T1ength=Totwa11+Totair1+Totair2+Totair3 :totwa11=0:
totair1=0:totair2=0ztotair3=0znnorma1=0

56220 Rsystem=TR:TR=0

59000 Return

106

(CHOI5.BAS)

10
20
30
40

OPEN "0" , #1 , "ChO'i 5 . FES"

NX=3

: Read ND

either side of the wa11s is foi1ed

DIM X(ND,NX),Y(ND),C(NX,NX+1), Residua1(ND), nnorma1(ND),
dnorma1(ND,30) 'store data in arrays X and Y

50 FOR I=1 TO ND 'read data in as x1,x2,...y
60 READ Y(I)

70 FOR J=1 TO NX : READ X(I,J) : NEXT 3
200 NEXT I

990 data 36

1000 data 5.891172 , 1.5 , 0, O

1010 data 11.75717 , 3 , 0,0

1020 data 12.26155 , 3, 0, .1

1030 data 13.29638 , 3 , , .3

1040 data 21.17806 , 4.5 ,0, .7000001
1050 data 4.460613 , .468 ,0, .5000001
1060 data 9.731501 , .78 ,0, 1.3

1070 data 10.76058 , .78 ,0, 1.4

1080 data 10.33307 , .78 ,0, 1.4

1090 data 10.32826 , .78 ,0, 1.4

1100 data 9.748288 , .78 ,0, 1.2

1110 data 8.701786 , .78 ,0, 1.1

2000 data 5.851358 , 1.5 , 0 ,0

2010 data 11.70232 , 3 , 0 ,0

2020 data 12.26874 , 3 , .1 ,0

2030 data 13.35619 , 3 , .3 ,0

2040 data 21.20628 , 4.5 .7 ,0

2050 data 4.569338 , .468 , .5000001 ,0
2060 data 10.07744 , .78 , 1.3 ,0
2070 data 11.16092 , .78 , 1.4 ,0
2080 data 10.67635 , .78 , 1.4 ,0
2090 data 10.67344 , .78 , 1.4 ,0
2100 data 10.10647 , .78 , 1.2 ,0
2110 data 9.004424 , .78 , 1.1 ,0

3000 data 5.851358 , 1.5 , 0 , 0

3010 data 11.70232 , 3 , O , 0

3020 data 12.26874 , 3 , .1 , 0

3030 data 13.35619 , 3 , .3 , 0

3040 data 21.20628 , 4.5 , .2 , .5
3050 data 4.561962 , .468 , .4 , .1
3060 data 10.04717 , .78 , 1.1 , .2
3070 data 11.03265 , .78 , .8000001 , .6
3080 data 10.57067 , .78 , 1 , .4
3090 data 10.48752 , .78 , .7000001 , .7000001
3100 data 9.846064 , .78 , .6 , .6
3110 data 8.737776 , .78 , .3 , .8
4000 For i= 1 to NO

4100 For j= 1 to NX

4200 print #1, USING" #####.##";x(i,j);
4300 Next j

4400 print #1, USING" #####.##";y(i);
4500 Print #1,

4600

Next i

107

 

5000
5100
5110
5120
5130
5140
5150
5160
5170
5180
5190

5200
5205
5210
5220
5230
5235
5240

9980
9990

toterr=0

FOR I=1 TO NX : FOR J= 1 TO NX : C(I, J)=0

FOR K=1 TO ND : C(I, J)=C(I, J)+X(K, I) *XCK, J)
NEXT J : NEXT I

FOR I=1 TO NX 2

: NEXT K
C(I,NX+1)=0

FOR K=1 TO ND : C(I,NX+1)=C(I,NX+1)+X(K,I)*Y(K) : NEXT K
NEXT I
GOSUB 10240 'so1ve NX x NX sgstem of equations
PRINT #1, " y = a1*x1 + a2*x +a3*x3 ...
FOR I=1 TO NX : PRINT #1, " a";I;"= ";CCI,NX+1) : NEXT I
PRINT #1, : PRINT #1, " fitted y given y Residua1"
FOR I=1 TO ND
YI=0 'check fit ...
FOR J=1 TO NX : YI=YI+X(I,J)*C(J,NX+1): NEXT J
Residua1(i)=Y(i)-YI : Toterr=toterr+Residua1(i)A2
PRINT #1, USING" ###.### ";YI;Y(I);Residua1(i)
NEXT I
PRINT #1, " SSE: ";Toterr
CLOSE #1
END

10240 'matrix subroutine

10250 N=NX :
10260 PVT=I
10270 IF ABSCCCK, I))>ABS(C(PVT, I)) THEN PVT=K

FOR I=1 TO N 'zero out Ith co1umn
: FOR K=I+1 TO N 'search for pivot row first
'pivot row has

1ar est number

1028

10290 FOR J=I TO N+1
10300 CHG=C(I,J)

NEXT K .
'1nterchange row I and pivot row
: C(I,J)=C(PVT,J : C(PVL J)=CHG

10310 NEXT J

10320 FOR K=1 TO N

: IF K=I THEN 10350 'use pivot row to zero

out Ith co1umn

10330

POWS

10340
10350
10360
10370
10380

F=C(K,I)/C(I,I) 'subtract F times ror I from a11 other
FOR J=I+1 TO N+1 : C(K,J)=C(K,J)—F*C(I,J) : NEXT J

NEXT K

NEXT I 'zero out next co1umn

FOR K=1 TO N : C(K,N+1)=C(K,N+1)/C(K,K) : NEXT K

RETURN

108

(CHOIGBAS)

10 ' no foi1 at the wa11

20 OPEN "O",#1,"choi6.res"

30 NX=2 : Read ND

40 DIM X(ND, NXL Y(NDL C(Nx, NX+1L Residua1(ND), nnorma1(ND),
dnorma1(ND, 30) 'store data in arrays X and Y

45(2) PRINT #1, : PRINT #1, " beta a(1)

a

50 FOR I=1 TO ND 'read data in as x1, x2, .y

60 READ Y(I)

70 FOR J= 1 TO NX-l : READ X(I,J) : NEXT J

80 READ nnorma1 (i)

90 If nnorma1(i) =0 then 200

100 For j=1 to nnorma1(i)

110 Read dnorma1(i,j)

130 Next j

200 NEXT I

990 data 30

1000 data 5.851358 , 1.5 , 0

1010 data 11.70232 , 3 , 0

1020 data 12.08868 , 3 , 1 , .1

1030 data 12.08868 , 3 , , .1

1040 data 18.78432 , 4.5 , 2 , .1 , .5

1050 data 2.795854 , .468 , 2 , .1 , .2

1060 data 4.981812 , .78 , 4 , .1 , .2 , .1 ,
1070 data 4.981812 , .78 , 4 , .1 , .2 , .1 ,
1080 data 4.981812 , .78 , 4 , .1 , .2 , .1 ,
1090 data 4.981812 , .78 , 4 , .1 , .2 , .1 ,
1100 data 4.981812 , .78 , 4 , .1 , .2 , .1 ,
1110 data 4.587751 , .78 , 4 , .1 , .1 , .1 ,
1120 data 8 910006,1.56,9,.1,.1,.1,.1,.1,.1,.1,.1,.1
1130 data 16.15394,3.248,9,.1,.1,.1,.1,.1,.1,.1,.1,.1
1140 data 9 234108,1.S6,9,.1,.1,.1,.1,.1,.1,.1,.1,.1
1150 data 9 561256,1.56,9,.1,.1,.1,.1,.1,.1,.1,.1,.1
1160 data 14.93095,2.936 ,9,.1,.1,.1,.1,.1,.1,.1,.1,.1
1170 data 16.69534,2.936 ,9,.2,.2,.2,.2,.2,.2,.2,.2,.2
1180 data 12.69694 , 3 , 1 , 1

1190 data 12. 69694 , 3 , 1 , 1

1200 data 2. 795854 , .468 , 2 , .1 , .2

1210 data 4. 981812 , .78 , 4 , .1 , .2 , .1 ,
1220 data 4.981812 , .78 , 4 , .1 , .2 , .1 ,
1230 data 4. 981812 , .78 , 4 , .1 , .2 , .1 ,
1240 data 4. 981812 , .78 , 4 , .1 , .2 , .1 ,
1250 data 9. 426338, 1. 56, 9, .1, .2,.1,.1,.1,.1,.1,.1,.1
1260 data 16.15394, 3.248,9, .1,.1,.1,.1,.1,.1,.1,.1,.1
1270 data 9.644621,1. 56,9, .L .1,. ,.1,.1,.1,.1,.4,.4
1280 data 11. 89268, 1. 904,9, .1,. 1,. 1, .1,.1,.1,.4,.4,.2
1290 data 14.93095, 2.936, 9, .1, .1,. 1, .1,.1,.1, 1,.1,.1
2100 For beta=0. 210 to 0. 240 step 0. 001

2200 For i =1 to ND

2300 If nnorma1(i) =0 then 2700

2400 For j= —1 to nnorma1 (i)

2500 Rnorma1= Rnorma1 + dnorma1(i , j) /

(beta +dnorma1(i

. 1))

109

I—‘NNNNN

NNNN

2600
2700
2900

5000
5100
5110
5120
5130
5140
5150
5160
5200
5205
5210
5220
5235
5240
5250
###.
5260

9000

9980
9990

. Next '
x(1,nx)=Rnorma : Rnorma1=0
NEXT I
toterr=0
FOR I=1 TO NX : FOR J=1 TO NX : C(I,J)=O

FOR K=1 TO ND : C(I,J)=C(I,J)+X(K,I)*X(K,J)
NEXT J : NEXT I

FOR I=1 TO NX : C(I,NX+1)=0

FOR K=1 TO ND 2 C(I,NX+1)=C(I,NX+1)+X(K,I)*Y(K)
NEXT I

GOSUB 10240
FOR I=1 TO ND

: NEXT K

: NEXT K

'soTve NX x Nx system of equations

YI=O 'check fit ...
FOR J=1 TO NX : YI=YI+X(I,J)*C(J,NX+1): NEXT J
Residua1(i)=Y(i)-YI : Toterr=toterr+Residua1(i)A2
NEXT I
Print #1, using " #. ###"; beta;
FOR I: 1 TO NX : PRINT #1, u51ng "
#### " ;C(I, NX+1);" : NEXT I

Print #L using #.#######";toterr
Next beta

CLOSE #1
END

10240 'matrix subroutine

10250 N=NX :
10260 PVT=I
10270 IF ABSCCCK, I))>ABS(C(PVT, I)) THEN PVT=K

FOR I=1 To N 'zero out Ith co1umn
: FOR K=I+1 To N 'search for pivot row first
'pivot row has

Tar est number

1028

10290 FOR J: I TO N+1
10300 CHG=C(I, J)

NEXT K .
'1nterchange: row I and pivot row
: C(I,J)=C(PVT,J C(PVT, J)=CHG

10310 NEXTJ

10320 FOR K=1 TO N : IF K=I THEN 10350 'use pivot row to zero
out Ith coTumn

10330 F=C(K,I)/C(I,I) 'subtract F times ror I from a11 other rows
10340 FOR J=I+l TO N+1 : C(K,J)=C(K,J)-F*C(I,J) : NEXT J

10350 NEXT K

10360 NEXT I 'zero out next coTumn

10370 FOR K=1 TO N : C(K,N+1)=C(K,N+1)/C(K,K) : NEXT K

10380 RETURN

110

(KCARGOQBAS)

5 rem without bottom contact a11 convection
10 OPEN "0",#1,"kcargoO9.res"

20 READ CRITERIA, MAXITER

25 DATA 0.005,100

27 Print #1, " Dataset A1 A2 A4 A5 A7 A8 x
h1 hradl h4 hrad4 h5 hradS h8 hrad8 k14 k25
k36 Qexp Qsim ErrorAZ %error "

28 Print #1,

30 N=11: Unknown=3: inicon=0.015

50(30)IM 9.096 )TACN) , L(3 , 3), L1(3 , 3) , A(N) , KCN) , HCN) , BUFFERCN) ,
E , ra N

51 DIM HA(3,3), C(unknown,unknown+1), VVVCN), VVVICN), TEMP(N),
TEMPlCN)

53 BOLZ=1.714E-O9

54 Dataset=38z Dim QexpCDataset), $SE(dataset+1), perererataset)
55 For 1: 1 to Dataset: read Qexp(1): next i

56 data 11.09, 21.32, 26.73, 5.94, 15.49, 33.84, 44.63, 11.23,
25.66, 58.31

58 data 79.02, 15.08, 20.13, 57.75, 69.69, 13.07, 41.04, 86.39,
123.04, 19.80

5199 doa4ta 45.15, 43.10, 42.17, 45.86, 43.64, 42.76, 17.39, 17.46,

60 data 24.68, 23.93, 25.50, 15.39, 15.69, 15.55, 12.96, 12.96,
13.62
61 For zzzz=1 to Dataset

63 Read property, Nbucket

64 REM *****reading dimensions and properties of the wa11s***
65 FOR I=1 TO 3

70 FOR J=1 T0 3

80 READ x

90 L(I,J)= x/12

100 NEXT J

110 READ E(i)

120 NEXT I

123 If 2222 > 12 and 2222 < 17 then 165

125 x = (L(3,1)-L(2,1))/2

130 A(8)=(3.141592*(L(1,1)+L(1,2))/2*L(1,3)
+3.141592/4*L(1,1)A2)*Nbucket

13S A(7)=(3.141592/4*L(1,1)A2)*Nbucket: A(6)=A(7): A(3)=A(7)
140 A(5)=L(2,1)*L(2,2)—A(6)

145 A(4)=2*L(2,3)*(L(2,1)+L(2,2))+L(2,1)*L(2,2)
150 A(2)=L(3,1)*L(3,2)-A(3)

155 A(1)=2*L(3,3)*(L(3,1)+L(3,2))+L(3,1)*L(3,2)
156 A(9)=A(1)+A(2)+A(3)

157 A(10)=A(4)+A(5)+A(6)

158 A(11)=A(8)+A(6)

160 goto 200

111

165
170
175
180
182
184
186
187
188
189

190

200
210
215
220
225
230

250
255
260
270
280
290
295
300

350
380
400
410
420
450
460
470
500
510
520
530
540
550
560
570
580
590
600
610
650
660
670
710
720
730
750

x = (L(3.1)- -L(2.1))/2
A(8)=(2*L(1,3)*(L(1,1)+L(1,2))+L(1,1)*L(1,2))*Nbucket
A(7)=(L(1,1)*L(1,2))* Nbucket: A(6)=A(7): A(3)=A(7)
A(5)=L(2,1)*L(2, 2)-A(6)
A(4)=2*L(2,3)*(LL(2,1)+L(2,2))+L(2,1)*L(2, 2)
A(2)=L(391)*L(3!2)-A

A(1)=2*L(3, 3)*(L(3,1)+L(3,2))+L(3,1)*L(3, 2)
A(9)=A(1)+A A(2)+A(3)

A(10)=A(4)+A(5)+A(6)

A(11)=A(8)+A(6)

REM ***************************************************

REM **********reading known temperature**************
Tice= 32:

Read Tinfinite

Tice=Tice+460: Tinfinite=Tinfinite+460

Tg=Tinfinite

REM ******************************************************

REM *********initia1 uess of each temperature*************
DELTA=(Tinfinite—Tice /4

T(2)=Tinfinite: T(3)=Tinfinite

T(1)=Tinfinite-De1ta

T(4)=Tinfinite- 2*De1ta: T(5)=T(4): T(6)=Tice

T(7)=Tice: T(8)=Tice

For i= 1 to N. BufferCi)=T(i): TA(i)=T(i): next I

REM *******************************************************

FOR KK=1 TO MAXITER
GOSUB 2000
FOR J: 1 TO unknown
VVVCJ) = T(J)-BUFFER(J)
NEXT J
FOR J= 1 TO unknown
TBUFFER = TBUFFER + ABSCCTCJ)-BUFFER(J))/T(J))
NEXT J
FOR III= 1 TO 9
FOR J: 1 TO unknown
TA(J)=BUFFER(J)+VVV(J)*III/lO
NEXT J
GOSUB 2000
FOR J: 1 TO unknown
TERR = TERR + ABS((T(J)-TA(J))/T(J))
NEXT J
IF TERR >= TBUFFER THEN 660
FOR J: 1 TO unknown
TEMP(J)=TA(J): TEMPlCJ)=T(J)
NEXT J
TBUFFER = TERR: Q=QTEMP
TERR=0
NEXT III
FOR J: 1 TO unknown. TA(J)=TEMP(J): BUFFER(J)=TA(J): NEXT J
IF TBUFFER < CRITERIA THEN 930
TBUFFER: 0
NEXT KK

112

910

920
930 REM conver
964
965
970 gosub 50000

972 next 2222

973 Print #1,

REM did NOT converge
PRINT "not converged": GOTO 1998

"Tota1 555:";SSECdataset+1)=

ed
gosub 60008: Rem Ca1cu1ation of SSE
M ***************Ca1cu1at-ion of h and k*******

sse(dataset+1)=0

974 REM ‘k***********************DATA***************************

975 DATA 2, 1, 6, 6, 5, 0.95, 7.25, 7.25, 7.5, 0.95,

9.25, 9.25, 9.5, 0. , 62

976 DATA 2, 1, 6, 6, 5, 0.95, 7.25, 7.25, 7.5, 0.95,

9.25, 9.25, 9.5, 0.95, 90

977 DATA 2, 1, 6, 6, 5, 0.95, 7.25, 7.25, 7.5, 0.95,

9.25, 9.25, 9.5, 0.95, 104

978 DATA 2, 1, 4.5, 4 5, 2.5 0 95 7.25, 7.25, 7.5, 0.95,
9.25, 9.25, 9.5, 0.95, 62

979 DATA 2, 1, 8, 8, 6.5, 0.95, 11.0, 10.5, 9.0, 0.95,
13.0, 12.5,1 .0, 0. 5, 62

980 DATA 2, 1, 8, 8, 6.5, 0.95, 11.0, 10.5, 9.0, 0.95,
13.0, 12.5, 1 .0, 0. , 90

981 DATA 2, 1, 8, 8 6.5, 0.95, 11.0, 10.5, 9.0, 0.95,
13.0,12.5,1 .0, 0.95, 104

982 DATA 2,1,6, 6, 5.0, 0.95, 11.0, 10.5, 9.0, 0.95,
13.0,12.5,11.0, 0.95, 62

983 DATA 2, 2,8, 8, 6.5, 0.95, 20.5, 13.75, 8.75, 0.95,
22.5, 15.75, 10.75, 0.95, 62

984 DATA 2, 2, 8, , 6.5, 0.95, 20.5, 13.75, 8.75, 0.95,
22.5, 15.75, 10.75, 0.95, 90

985 DATA 2, 2, 8, 8, 6.5, 0.95, 20.5, 13.75, 8.75, 0.95,
22.5, 15.75, 10.75, 0.95, 104

986 DATA 2, 1, 6, 6, 5.0, 0 95, 20.5, 13.75, 8.75, 0.95,
22.5, 15.75, 10. 75, 0.95, 62

987 DATA 2,2,11, 7.5, 2.25, 0 95 24.75, 9.75, 3, 0.95,
26.75,11.75,5.0,0.95, 62

988 DATA 2,2, 11, 7. 5, 2.25, 0.95, 24.75, 9.75, 3, 0.95,
26. 75, 11. 75, 5.0, 0.9L 90

989 DATA 2,2, 11, 7.5, 2.25, 0.95, 24.75, 9.75, 3, 0.95,
26.7S,11.75, 5.0, 0.95, 104

990 DATA 2,1, 8, 6, 2.25, 0.95, 24.75, 9.75, 3, 0.95,
27.75, 11.75, 5.0, 0.95, 62

991 DATA 2, 4, 8, 8, 5.0, 0.95, 41,10, 6.75, 0.95, 43, 12,
8.75, 0.95, 62

992 DATA 2, 4, 8, 8, 5.0, 0.95, 41,10, 6.75, 0.95, 43, 12,
8.75, 0.95, 90

993 DATA 2, 4, 8, 8, 5.0, 0.95, 41,10, 6.75, 0.95, 43, 12,
8. 75, 0. 95, 104

994 DATA 2,1,9.5, 4.5, 5.0, 0.95, 41,10, 6.75, 0.95, 43,
12, 8. 75, 0. 95, 62

995 DATA 1,1, 6, 6, 4.3, 0.95, 13,9,14, 0.95, 13.312,
9.312, 14.312, 0.95, 72

996 DATA 1, 1, 6, 13,9,14, 0.95, 13.312,

9.312, 14.312, 0.95,

6, 4.1, 0.95,
72

113

997 DATA 1, 1, 6, 6, 4.0, 0.95, 13,9,14, 0.95, 13.312,
9. 312, 14. 312, 0.9L 72

998 DATA 1,1, 6, 6, 4. L 0.9L 15, 12 ,10, 0.95, 15.312,
12.312, 10.312, 0.95, 72

999 DATA 1, L 6, 6, 3. 9, 0. 95, 15, 12 ,10, 0.95, 15.312,
IL 312, 10.312, 0.9L 72

1000 DATA 1,1, 6, 6, 3.8, 0.95, 15, 12 ,10, 0.95, 15.312,
12. 312, 10.312, 0.9L 72

1010 DATA 0,1, 6, 6, 4.7, 0.95, 12, 9, 10, 0.95, 15, 12,
13, 0. 95, 72

1020 DATA 0, 1, 6, 6, 4 7, 0.95, 12, 9, 10, 0.95, 15, 12,
13, 0.95, 72

1030 DATA 0, 1, 6, 6, 5 2, 0.95, 12, 9, 10, 0.95, 15, 12,
13, 0.95, 72

1040 DATA 0, 1, 6, 6, 4 3, 0.95, 14,11, 8.5, 0.95, 16, 13,
10.5, 0.95, 72

1050 DATA 0, 1, 6, 6, 4.2, 0.95, 14,11, 8.5, 0.95, 16, 13,
10.5, 0.95, 72

1060 DATA 0, 1, 6, 6, 4 2, 0.95, 14,11, 8.5, 0.95, 16, 13,
10. 5, 0. 95, 72

1070 DATA 3, L 6, 6, 4.0, 0.95, 10,10, 10, 0.95, 12.624,
12. 624, 12. 624, 0. 95, 72

1080 DATA 3, L 6, 6, 4.2, 0.95, 10,10, 10, 0.95, 12.624,
12.624,1L 624, 0. 95, 72

1090 DATA 3, 1, 6, 6, 4.2, 0.95, 10,10, 10, 0.95, 12.624,
1L 624,1L 624, 0. 95, 72

1100 DATA 4, 1, 6, 6, 4.5, 0.95, 10,10, 10, 0.95, 12.312,

12. 312, 1L 312, 0. 95, 72
1110 DATA 4, L 4.25, 0.95, 10,10, 10, 0.95, 12.312,

12. 312, 12.312, 0.9L 72
1120 DATA 4,1, 44, 0.95, 10,10, 10, 0.95, 12.312,

6,
IL 312, 12.312, 0.9L 72

1997 REM *******************************************************

mm
03

1998 CLOSE #1

1999 END
2000 REM **************majn subroutine*********************
2050 REM ****************Loop for materia1 pro erties*******

2100 Ta1=Ta(1): Ta2=Tinfinite: i=3: jj=9: osu 30100
2200 Ta1=Ta(2): Ta2=Tice: i=1: jj=—11:gosu 29020

2600 Ta1=Ta(1):Ta2=Ta(2): gosub 35000: k14=k

3000 REM ************Loop for construction of Matrix**********
3100 Gosub 10000 . .

4000 REM *****1mpose known var1ab1es and reduce matr1x*******
5000 REM ******************so1ving matrix********************
5100 GOSUB 20000

6000 REM ****************************************************

8400 QTEMP=C(unknown,unknown+1)
8500 FOR I=1 To unknown: T(I)=C(I,unknown+1): NEXT I

114

8600 FOR I=1 TO unknown: FOR J=1 TO unknown+1: C(I,J)=0
: NEXT J: NEXT I

8700 FOR I=1 TO N: HCI)=0: NEXT I
9000 RETURN

10000 REM ***********Construct1ng Matr1x**********************

11000 C(1,1)= h(9)*A(9): C(1, 3)=1
11100 C(1, 4): h(9)*A(9)*T1nf1n1te

11200 C(2,1)= h(9)*A(9)+k14*SQR(A(9)*A(10))/x : C(2,2)= -
k14*SQR(A(9)*A(10) /x
11300C C(Z ,4): h(9)*A(9)*T1nf1n1te

11600 C(3, 2)- -h(11)*A(11): C(3, 3): 1
11700 C(3, 4)- -h(11)*A(11)*T1ce

13000 Return

20000 REM *************subrout1ne for so1v1ng matr1x*************
20010 REM: SOLVES SYSTEM OF EQNS BY GAUSSIAN ELIMINATION

20060 FOR I=1 TO unknown : PVT=I 'e11m1nat1on using row I ...
20070 FOR K=I+1 TO unknown

20080 IF ABSCCCK, I))>ABS(C(PVT, I)) THEN PVT=K

20090 NEXT K 'row "PVT" has '1argest e1ement 1n Ith co1umn

20100 FOR J= I TO unknown+1'1nterchan e row I and p1vot row

20110 CHG=C(I, J) : C(I, J)=C(PVT, J): C(PVT, J)=CHG

20120 NEXT 3 'now row I has 1argest e1ement 1n Ith co1umn

20130 FOR K=1 TO unknown : IF K=I THEN 20160

20140 F=C(K, I)/C(I, I)

20150 FOR J= I+1 TO unknown+1 : C(K, J) =C(K, J)- F*C(I, J) : NEXT J

20160 NEXT K 'Ith co1umn now zeroed out (except for row I)

20170 NEXT I

20180 FOR K=1 TO unknown:C(K,unknown+1)=C(K,unknown+1)/C(K,K)
: NEXT K

20200 RETURN

29020 REM ****convect1on coeff1c1ent for enc1osed space******
29050 TAVE=(TA1+TA2)/2- 460
29060 NPR=. 07206- 1. 8453E— 04*TAVE+3.5979E-07*TAVEA2
-2.9967E-10*TAVEA3
29065 KAIR=.0131+2.6136E-05*TAVE-6.1955E—09*TAVEA2
29066 L1(I,1)=ABS((L(2, 3)-L(1, 3)))/2
29067 L1(I, 2)=ABS((L(2, 3)- -L(1, 3)))/2
29068 L1(I, 3)=(ABS((L(2,1)-Nbucket*L(1,1))/2)
+ ABS((L(3, 2)- Nbucket*L(1, 2))/2))/2

29070 FORJ

29080

=1 TO 3

NGRCJ)= L1(I, J)A3*ABS(TA2- TA1)*

(.1794+4.1794*EXP(- .0095034*TAVE))*1E6

29090 IF J= 1 THEN GOSUB 29400

29095 IF j=2 then gosub 29300

29100 IF J=3 THEN GOSUB 29200

29110 HA(1,J)= KAIR/LlCI,J)*A*NGR(j)AM

115

29120
29122
29124

hcj'
1.11
29130

HA(1,2)

NEXT J
If 2222 > 12 and zzzz < 17 then 29130
If jj= 11 then

)=((HA(1,1)+HA(1,2))*3.141592/4*L(1,2)A2+HA(1,3)*3.141592*L(
*L(1!3))/A(jj): oto 29135

h(jg)=((HA(1,1 +HA(1,2))*L(1,1)*L(1,2)+
* *L(1.3)*(L(1.1)+L(1.2)))/A(jj)

29135 hradej)=bo12*(TalA3+Ta1A2*Ta2+Ta1*Ta2A2+Ta2A3)
/( 1/E(1) + l/ECZ) -1)

29140
29150

29200
29210
29215
29220
29230
29270
29290

29300
29310
29315
29320
29330
29350
29390

29400
29450
29490

30100
30500
30600

h(jj)=h(jj)+hrad(jj)
RETURN

REM convect1on coeff1c1ent of a1r for vert1ca1 p1ane
If NGRCj)>2000 then 29220
A=1: M=0.0: oto 29290
IF NGRCg)<=2000 THEN GOTO 29270
A=.1 *(L(I,3)/L1(1,j))A(-1/9):M=1/4:GOTO 29290
A=.065*(L(I,3)/L1(1,J))A(-1/9):M=1/3
RETURN

REM convection coeff1c1ent of a1r for Bottom
If NGR(j)>2000 then 29320
A=1: M=0.0: oto 29390
IF NGR(g)<=4E+0 THEN GOTO 29350
A=.1 5:M=1/4:GOTO 29390
A=.068:M=1/3
RETURN

REM convect1on coeff1c1ent of a1r for Top
A=1.0:M=0.0
RETURN

REM conduct1v1ty and FREE convect1on coeff1c1ent for a1r
TAVE=(TA1+TA2)/2-460
NPR=.07206-1.8453E-04*TAVE+3.5979E

-07*TAVEA2-2.9967E-10*TAVEA3

30650
30660
30680
30700
30800

KAIR=.0131+Z.6136E-05*TAVE-6.1955E-09*TAVEA2
L1(I,1)=(L(I,1)+L(I,2))/2:L1(1,2)=L1(1,1)
L1(I,3)=L(I,3)

FOR J=1 TO 3

NGRCJ)=L1(I,J)A3*ABS(TA2-TA1)*

(.1794+4.1794*EXP(—.0095034*TAVE))*1E6

30850
30900
31000
31050
31100
31200
31300

+HA(1,2

31350
31400
31500

NPRNGR=NPR*NGR(J)
IF J=1 THEN GOSUB 34000
If J=2 THEN GOSUB 33000
IF J=3 THEN GOSUB 32000
HA(1,J)= KAIR/LICI,J)*A*NPRNGRAM
NEXT 3
h('g)=((HA(i.1)+Ha(1.2))*L(1.1)*y§1.2)
3 2*L<1.3)*(L<1.1>+L<1.2)>)/A<3 >
hrquja)fE(1)*b01If(TA2A4-TA1A4)/%TA2-TA1)
h(JJ)= (JJ)+hrad(JJ)
RETURN

116

32000
32100
32200
32300
32500
32700
32900

33000
33100
33300
33500
33900

34000
34500
34900

35000
35500
35550
35560
35570
35580
35590
35640
35650
35660
35670
35680
35690
36000

REM convection coefficient of air for vertica1 p1ane
IF NPRNGR<1E+04 THEN GOTO 32500

IF NPRNGR>1E+09 THEN GOTO 32700

A=.59:M=1/4:GOTO 32900

A=1.36:M=1/5:GOTO 32900

A=.13:M=1/3
RETURN

REM convection coefficient of air for Bottom
IF NPRNGR<2E+07 THEN GOTO 33500
A=.14:M=1/3:GOTO 33900
A=.S4:M=1/4
RETURN

REM convection coefficient of air for Top
A=.58:M=1/3
RETURN

REM *****conduction coefficients for each materia1s******
TAVE=(TA1+TA2)/2-460

If property = 2 then 35650

If property = 1 then 35660
If property = 0 then 35670
If property = 3 then 35680
If property = 4 then 35690

k=0.346: goto 36000 'water
k=.0131+2.6136E-05*TAVE-6.19SSE-09*TAVEA2:goto 36000
k=0.035: goto 36000 'corrugated board
k=0.022: goto 36000 'EPS foam
k=0.018: goto 36000 'Po1yurathane

=0.010: goto 36000 'VIP

k=k*(.013 +2.6136E-05*TAVE-6.1955E-09*TAVEA2)

/(.0131+2.6136E-05*72-6.1955E-09*72A2)

39900

40000
40100
40200

RETURN

REM ****convection coefficient for enc1osed space******
TAVE=(TA1+TA2)/2-460
NPR=.07206-1.8453E-04*TAVE+3.5979E-07*TAVEA2-2.9967E—

10*TAVEA3

40300
40400
40500

KAIR=.0131+2.6136E-05*TAVE-6.1955E-09*TAVEA2
L1(I,1)=L(2,3)
NGR=L1(I,1)A3*ABs(TA2-TA1)*(.1794+4.1794*EXP(-

.0095034*TAVE))*1E6

40600
40700
40800
40900
41000
41025
41050
41100

50000
51100
51200
51600
51810

REM convect1on coefficient of air for Top

IF NGR<=4E+05 THEN GOTO 40800
A=.195:M=1/4:GOTO 40900
A=.068:M=1/3

H(jj)= KAIR/LlCI,1)*A*NGRAM

hrad(ja)fE(1)*b012f(TA2A4-TA1A4)/(TA2-TA1)

h(JJ)= (JJ)+hrad(JJ)

RETURN

Rem ********Ca1cu1ation of Fina1 va1ues of h's and k's*****
Ta1=Temp1(1): Ta2=Tinf1nite: 1:3: jj=9: osub 30100
Ta1=Tice: Ta2=Ta(2): i=1: jj=11: osub 2 020
Ta1=Temp1(2): Ta2=Temp(1): gosub 5000: k14=k
SSE(zzzz)=(Qexp(zzzz)-Q)A2

: pererrszzz)=(Qexp(zzzz)-Q)/Qexp(zzzz)*100

117

51890 SSE(dataset+1)=SSE(dataset+1)+SSE(zzzz)
52300 Print #1, usin " ###.#### ";
zzzz;A(1);A(Z);A(4);A );A(7);A(8);x;h(9);hrad(9);h(11);hrad(11
);k14;k25;k36;Qexp(zzzz);Q;SSE(zzzz):pererr(zzzz): Print #1,
59000 Return

60000 Rem ****************Ca1cuat-ion of SSE ******~k~k**~k*~k*~k*****
62700 Return

118

(KCARG16.BAS)
10 OPEN "O",#1,"kcargolG.res"

20 READ CRITERIA, MAXITER

25 DATA 0.005,100

27 Print #L " Dataset A1 A2 A4 A5 A7 A8 x
h1 hradl h4 hrad4 hS hradS h8 hrad8 k14 k25
k36 Qexp Qsim ErrorAZ %error "

28 Print #1,

30 N=8: Unknown=3: inicon=0.015

50 DIM T(N), TA(N), L(3,3), L1(3,2), A(N), KCN), HCN), BUFFERCN),
E(3), hrad(N)

51 DIM HA(3,2), C(unknown,unknown+1), VVVCN), VVV1(N), TEMP(N),
TEMPlCN)

53 BOLZ=1. 714E- 01

54 Dataset: 38: Dim QexpCOataset), SSECdataset+1), perererataset)
55 For 1= 1 to Dataset: read QexpCi): next 1

56 data 11.09, 21. 32, 26. 7L 5. 94, 15. 49, 33. 84, 44.63, 11.23,
25. 66, 58. 31

58 data 79.02, 15.08, 20.13, 57.75, 69.69, 13.07, 41.04, 86.39,
123.04, 19.80

5199d§1ta 45.15, 43.10, 42.17, 45.86, 43.64, 42.76, 17.39, 17.46,
6103 data 24.68, 23.93, 25.50, 15.39, 15.69, 15.55, 12.96, 12.96,

61 For zzzz=1 to Dataset

63 Read property, Nbucket

64 REM *****reading d1mensions and properties of the wa11s****
65 FOR I=1 TO 3

70 FOR J=1 TO 3

80 READ X

90 L(I,J)= x/12

100 NEXT 3

110 READ E(i)

120 NEXT I

123 If 2222 > 12 and 2222 < 17 then 165

125 x = (L(3,1)-L(2,1))/2

130 A(8)=(3.141592*(L(1,1)+L(1,2))/2*L(1,3)
+3.141592/4*L(1,1)A2)*Nbucket

135 A(7)=(3.141592/4*L(1,1)A2)*Nbucket: A(6)=A(7): A(3)=A(7)
140 A(5)=L(2,1)*L(2,2)-A(6)

14S A(4)=2*L(2,3)*(L(2,1)+L(2,2))+L(2,1)*L(2,2)
150 A(2)=L(3,1)*L(3,2)-A 3

155 A(1)=2*L(3,3)*(L(3,1)+L(3,2))+L(3,1)*L(3,2)
160 goto 200

165 = (L(3.1)-L(2.1))/2

170 A(8)=(2*L(1, 3)*(L (1,1)+L(1,2))+L(1,1)*L(1,2))*Nbucket
175 A(7)=(L(1,1)*L(1,2))*Nbucket: A(6)=A(7): A(3)=A(7)
180 A(5)=L(2,1)*L(2,2)-A(6 )

182 A(4)=2*L(2,3)*(L(2,1)+L(2,2))+L(2,1)*L(2, 2)

184 A(2)=L(3,1)*L(3,2)- A(3 )

119

186
190

200
210
215
220
225
230

250
255
270
280
295
300

350
380
400
410
420
450
460
470
500
510
520
530
540
550
560
570
580
590
600
610
650
660
670
710
720
730
750

910
920

930
964
965
970

A(1)=2*L(3, 3)*(L(3,1)+L(3,H*;: (3 ,1)*L(3,2)

REM ********************* H* ***************************

REM **********read1ng known temperature*******************
Tice=32:

Read Tinfinite _ _ . . . .

T1CE?TTEEf460: T1nf1n1te=T1nf1n1te+460

Tg=T1nf1n1te

REM ******************************************************

REM *********initia1 uess of each temperature*************
DELTA=(Tinfinite-Tice /3

T(1)=Tinfinite- De1ta

T(2)=Tinfinite- 2*De1ta

For i= 1 to unknown. BufferCi)=T(i): TA(i)=T(i): next I

REM *******************************************************

FOR KK=1 T0 MAXITER
GOSUB 2000
FOR J: 1 TO unknown
VVVCJ) = T(J)-BUFFER(J)
NEXT J
FOR 3= 1 TO unknown
TBUFFER = TBUFFER + ABSCCTCJ)-BUFFER(J))/T(J))
NEXT J
FOR III= 1 TO 9
FOR J: 1 TO unknown
TA(J)=BUFFER(J)+VVV(J)*III/lO
NEXT J
GOSUB 2000
FOR 3= 1 TO unknown
TERR = TERR + ABSCCTCJ)-TA(J))/T(J))
NEXT J
IF TERR >= TBUFFER THEN 660
FOR J: 1 TO unknown
TEMP(J)=TA(J): TEMPlCJ)=T(J)
NEXT J
TBUFFER = TERR: Q=QTEMP
TERR=0
NEXT III
FOR J=1 TO unknown:TA(J)=TEMP(J):BUFFERCJ)=TA(J):NEXT J
IF TBUFFER < CRITERIA THEN 930
TBUFFER=O
NEXT KK

REM did NOT converge
PRINT "not converged": GOTO 972

REM converged .

gosub 6000 : Rem Ca1cu1ation of SSE

REM ***************Ca1cu1ation of h and k*******
gosub 50000

972 next 2222

973 Print #1, "Total SSE=";SSE(dataset+1): sse(dataset+1)=0

120

974 REM *****************‘k**DATA*******************************

975 DATA 2, 1, 6, 6, 5, 0.95, 7.25, 7.25, 7.5, 0.95,
95 62

9.25, 9.25, 9 5, O. ,

976 DATA 2, 1, 6, 6, S, 0.95, 7.25, 7.25, 7.5, 0.95,
9.25, 9.25, 9 S, 0.95, 90

977 DATA 2, 1, 6, 6, S, 0.95, 7.25, 7.25, 7.5, 0.95,
9.25, 9.25, 9.5, 0.95, 104

978 DATA 2, 1, 4.5, 4.5, 2.5, 0.95, 7.25, 7.25, 7.5, 0.95,
9.25, 9.25, 9.5, 0.95, 62

979 DATA 2, 1, 8, 8, 6.5, 0.95, 11.0, 10.5, 9.0, 0.95,
13.0, 12.5, 1 0, 0.95, 62

980 DATA 2, 1, 8, 8, 6.5, 0.95, 11.0, 10.5, 9.0, 0.95,
13.0,12.5,11.0, 0.9 , 90

981 DATA 2,1,8, 8, 6.5, 0.95, 11.0, 10.5, 9.0, 0.95,
13.0,12.5,11.0, 0.95, 104

982 DATA 2,1, 6, 6, 5.0, 0.95, 11.0, 10.5, 9.0, 0.95,
13.0,12.5,11.0, 89 , 62

983 DATA 2,2,8, , 6.5, 0.95, 20.5, 13.75, 8.75, 0.95,
22. 5, 15. 75, 10. 7

984 DATA 2,2,8, 8, 6.5, 0.95, 20.5, 13.75, 8.75, 0.95,
22.5, 15.75, 10. 75, 0.95, 90

985 DATA 2, 2, 8, 8, 6. 5, O. 95, 20.5, 13.75, 8.75, 0.95,
22.5, 15.75, 10. 75, O. 95, 104

986 DATA 2, 1, 6, 6, 5. 0, 0. 95, 20.5, 13.75, 8.75, 0.95,
22.5, 15.75, 10. 75, O. 95, 62

987 DATA 2, 2, 11, 7.5, 2. 25, 0.95, 24.75, 9.75, 3, 0.95,
26.75,11.75,5.0, 0.95, 62

988 DATA 2, 2, 11, 7.5, 2.25, 0.95, 24.75, 9.75, 3, 0.95,
26.75,11.75, 5.0, 0.95, 90

989 DATA 2,2, 11, 7.5, 2.25, 0.95, 24.75, 9.75, 3, 0.95,
26.75.11.75, 5.0, 0.95, 104

990 DATA 2,1, , 6, 2.25, 0.95, 24.75, 9.75, 3, 0.95,
27.75, 11.75, 5. , 0.95, 62

000000900

991 DATA 2, 4, , 8, 5.0, 0.95, 41,10, 6.75, 0.95, 43, 12,
8.75, 0.95, 62
992 DATA 2, 4, , 8, 5.0, 0.95, 41,10, 6.75, 0.95, 43, 12,
8.75, 0.95, 90

993 DATA 2, 4,, 8, 5.0, 0.95, 41,10, 6.75, 0.95, 43, 12,
8. 75, 0. 95, 104

994 DATA 2,1,9.5, 4.5, 5.0, 0.95, 41,10, 6.75, 0.95, 43,
12, 8. 7S, 0. 95, 62

995 DATA 1,1, 6, 6, 4.3, 0.95, 13,9,14, 0.95, 13.312,
9.312, 14.312, 0.95, 72
996 DATA 1, 1, 6, 6, 4.1, 0.95, 13,9,14, 0.95, 13.312,
9.312, 14.312, 0.95, 72
997 DATA 1, 1, 6, 6, 4.0, 0.95, 13,9,14, 0.95, 13.312,

9. 312, 14. 312, 0. 95, 72

998 DATA 1, 1, 6, 6, 4.1, 0.95, 15, 12 ,10, 0.95, 15.312,
12.312, 10.312, 0.95, 72

999 DATA 1, 1, 6, 6, 3. 9, 0. 95, 15, 12 ,10, 0.95, 15.312,
12.312, 10.312, 0.95, 72

1000 DATA 1,1,6, 6, 3. 8, 0. 95, 15, 12 ,10, 0.95, 15.312,
12. 312, 10.312, 0. 95, 72
1010 DATA 0,1, 6, 6, 4. 7, 0.95, 12, 9, 10, 0.95, 15, 12,
13, 0.95, 72

121

110320 (I))Aél’: 07,2 1, 6, 6, 4.7, 0.95, 12, 9, 10, 0.95, 15, 12,
1030 DATA 0, 1, 6, 6, 5.2, 0.95, 12, 9, 10, 0.95, 15, 12,
13, 0.95, 72

1040 DATA 0, 1, 6, 6, 4.3, 0.95, 14,11, 8.5, 0.95, 16, 13,
10.5, 0.95, 72

1050 DATA 0, 1, 6, 6, 4 2, 0.95, 14,11, 8.5, 0.95, 16, 13,
10.5, 0.95, 72

1060 DATA 0, 1, 6, 6, 4 2, 0.95, 14,11, 8.5, 0.95, 16, 13,
10. L 0. 95, 72

1070 DATA 3, 1, 6, 6, 4.0, 0.95, 10,10, 10, 0.95, 12.624,
12. 624, 12. 624, 0. 9L 72

1080 DATA 3, 1, 6, 6, 4.2, 0.95, 10,10, 10, 0.95, 12.624,
12. 624, 12. 624, 0. 9L 72

1090 DATA 3, 1, 6, 6, 4.2, 0.95, 10,10, 10, 0.95, 12.624,
12. 624, 12. 624, 0. 9L 72

1100 DATA 4, 1, 6, 6, 4.5, 0.95, 10,10, 10, 0.95, 12.312,
IL 312, 12.312, 0.9L 72

1110 DATA 4, 1, 6, 6, 4. 25, 0.95, 10,10, 10, 0.95, 12.312,
1L 312, 12. 312, O. 95, 72

1120 DATA 4,1,6, 4.4,0.95, 10,10, 10, 0.95, 12.312,

6,
12. 312, 12. 312, O. 95, 72

1997 REM *******************************************************

1998 CLOSE #1

REM *************main subroutine*************************

REM ************Loop for mater1a1 propert1es***********
jj=1z gosub 30100

gosub 29020

5000: k14= k

REM ***********Loop for construction of Matr1x***********

REM ******1mpose known var1ab1es and reduce matr1x*******

REM *******************So1ving matrix********************

REM *****************************************************

FOR I=1 TO unknown: T(I)=C(I,unknown+1): NEXT I
FOR I=1 TO unknown: FOR J=1 TO unknown+1:C(I,J)=0: NEXT J

1999 END
2000
2050
2100 Ta1=Ta(1): Ta2=T1nf1n1te: 1= 3:
2200 Ta1=Ta(2): Ta2=T1ce: 1=2: j '=4:
2600 Ta1=Ta(1): Ta2=Ta(2): gosub
2700 Ta1=Tg: TaZ=Ta(2T gosub 35000: k25= k
2800 Ta1=Tg: Ta2=T1ce: gosub 35000: k36= k
3000
3100 Gosub 10000
4000
5000
5100 GOSUB 20000
6000
8400 QTEMP=C(unknown,unknown+1)
8500
8600
: NEXT I
8700 FOR I=1 TO N: H(I)=0: NEXT I
9000 RETURN

10000 REM ***********Constructing Matrix***********************

122

10010 radiation=bO1z/((A(4)+A(5))/A(8)*1/E(1)+1/E(2)-
1)*(A(4)+A(5))*(CTAC2)/100)A4-(Tice/100)A4)

11ooo c(1,1)= h(1)*(A(1)+A(2)+A(3)): c(1,3)=1

11688 EE%’§3: h(1)*(A(l)+A(2)+A(3))*Tinfinite
h(l)*(A(1)+;(2)+A(3))+k14*SQR((A(1)+A(2))*(A(4)+A(S)))/x+k36*SQR
1&ggg)2?g6%%£xz c(2,2)= -k14*5QR((A(1)+A(2))*(A(4)+A(5)))/x
h(1)*(AC1)+X(2)+A(3))*Tinfinite+k36*SQR(A(3)*A(6))/X*Tice

11600 c(3,1)= k14*SQR((A(1)+A(2))*(A(4)+A(5)))/X: c(3,2)= -
k14*SQR((A(1)+A(2))*(A(4)+A(5)))[x-h(4)*(A(4)+A(S))

11700 c(3,4)= -h(4)*(A(4)+A(5))*T1ce

13000 Return

20000 REM **********subroutine for so1ving matrix****************
20010 REM: SOLVES SYSTEM OF EQNS BY QAQSSIAN ELIMINATION
$8068 FOR I=1 TO unknown : PVT=I 'e11m1nat10n us1ng row I ...

07

FOR K=I+1 TO unknown

20080 IF ABSCCCK,I))>ABS(C(PVT,I)) THEN PVT=K

20090 NEXT K 'row "PVT" has 1argest e1ement in Ith co1umn
20100 FOR J=1 TO unknown+1 'interchan e row I and pivot row
20110 CHG=C(I,J) : C(I,J)=C(PVT,Jg : C(PVT,J)=CHG

20120 NEXT 3 'now row I has 1argest e1ement in Ith co1umn
20130 FOR K=1 TO unknown : IF K=I THEN 20160

20140 F=C(K,I)/C(I,I)

20150 FOR J=I+1 TO unknown+1 : C(K,J)=C(K,J)—F*C(I,J) : NEXT J
20160 NEXT K 'Ith co1umn now zeroed out (except for row I)
20170 NEXT I

20180 FOR K=1 TO unknown:C(K,unknown+1)=C(K,unknown+1)/C(K,K)

: NEXT K

20200 RETURN

29020 REM ****convection coefficient for enc1osed space******
29050 TAVE=(TA1+TA2)/2-460

23820 Ni§=.07206-1.8453E-04*TAVE+3.5979E-07*TAVEA2-2.9967E—

TAVE

29065 KAIR=.0131+2.61365-05*TAVE-6.1955E-09*TAVEA2

29066 L1(I,1)=ABS((L(2,3)-L(1,3)))

29068 L1(I,2)=(ABS((L(2,1)-Nbucket*L(1,1))/2) + ABS((L(3,2)-
Nbucket*L(1,2))/2))/2

29070 FOR J=1 T0 2

29080 NGRCJ)=L1(I,J)A3*ABS(TA2-TA1)*(.1794+4.1794*EXP(-
.0095034*TAVE))*1E6

29090 IF J=1 THEN GOSUB 29400

29100 IF J=2 THEN GOSUB 29200

29110 HA(1,3)= KAIR/L1(I,J)*A*NGR(j)AM

29120 NEXT J

29122 If 2222 > 12 and 2222 < 17 then 29130

29124 If jj= 8 then
h(jj)=(HA(1,1)*3.141592/4*L(i,1)A2+HA(i,2)*3.141592*L(i,1)*L(i,3
))/A(6j): oto 29135

29130 (

. gj =§HA(i,1)*t§i,1)*L(i,2)+HA(1,2)*2*L(T,3)
*(L(1.1 +L 1.2)))/A(JJ)

29135 hradCJ )=bo1z*((TA2/100)A4-(TA1/100)A4)
/(TA2-TA;)/3 1/gq1) + 1/E(i-1)*<A(4)+A(5))/A(81:1>
29140 h(JJ)=2*h(JJ)*A(8)/(A(8)+A(4)+A(5))+hrad(JJ)

123

29150

29200
29210
29215
29220
29230
29270
29290

29300
29305
29310
29330
29350
29390

29400

29450
29490

30100

RETURN

REM convection coefficient of air for vertica1 p1ane

If NGRCj)>2000 then 29220
A=1: M=0.0: oto 29290

IF NGRCj)<=2000 THEN GOTO 29270
A=.18*(L(I,3)/L1(i,j))A(-1/9):M=1/4:GOTO 29290
A=.065*(L(I,3)/L1(i,j))A(-1/9):M=1/3

RETURN

REM convection coefficient of air for Bottom
A1(1.J)?A(1.J)-A(1.J)
IF NGR(3)<=4E+05 THEN GOTO 29350
A=.1 5:M=1/4:GOTO 29390
A=.068:M=1/3
RETURN

REM convection coefficient of air for Top

A=1.0: =0.0
RETURN

REM subroutine for conductivity and FREE convection

coefficient for air

30500
30600

TAVE=(TA1+TA2)/Z-460
NPR=.07206-1.8453E-04*TAVE+3.5979E-07*TAVEA2-2.9967E-

10*TAVEA3

30650
30660
30680
30700
30800

KAIR=.0131+2.6136E-05*TAVE-6.1955E-09*TAVEA2
L1(I,1)=(L(I,1)+L(I,2))/2

L1(I,2)=L(I,3)

FOR J=1 TO 2
NGRCJ)=L1(I,J)A3*ABS(TA2-TA1)*(.1794+4.l794*EXP(-

.0095034*TAVE))*1E6

30850
30900
31050
31100
31200

NPRNGR=NPR*NGR(J)
IF J=1 THEN GOSUB 34000
IF J=2 THEN GOSUB 32000
HACT,J)= KAIR/LlCI,J)*A*NPRNGRAM
NEXT J

31300 h(jg)=(HA(i,1)*L(i,1)*L(i,2)+

HACi,2)*

31350
31400
31500

32000
32100
32200
32300
32500
32700
32900

33000
33100
33300
33500
33900
34000

*cgi,3)f(L(i.1)+L(i.2)))/A(jj)
hradCJfl)fE(1)*b01IfCCTA2/100)A4-(TA1/100)A4)/(TA2-TA1)
2é%3a: (JJ)+hrad(JJ)

REM convection coefficient of air for vertica1 p1ane
IF NPRNGR<1E+04 THEN GOTO 32500

IF NPRNGR>1E+09 THEN GOTO 32700

A=.59:M=1/4:GOTO 32900

A=1.36:M=1/S:GOTO 32900

A=.13:M=1/3
RETURN

REM convection coefficient of air for Bottom
IF NPRNGR<2E+07 THEN GOTO 33500
A=.14:M=1/3:GOTO 33900
A=.S4:M=1/4
RETURN
REM convection coefficient of air for Top

124

34500 A=.58:M=1/3
34900 RETURN

35000 REM *****conduction coefficients for each materia15*******
35500 TAVE=(TA1+TA2)/2—460

35550 If property = 2 then 35650

35560 If property = 1 then 35660

35570 If property = 0 then 35670

35580 If property = 3 then 35680

35590 If property = 4 then 35690

35640 k=0.346: goto 36000 'water

35650 k=0.016zgoto 36000 'cargotech

35660 k=0.035: goto 36000 'corrugated board
35670 k=0.022: goto 36000 'EPS foam

35680 k=0.018: goto 36000 'PoTyurathane
35690 k=0.010: goto 36000 'VIP

36000 k=k*(.013 +2.6136E-05*TAVE-6.1955E-09*TAVEA2)
/(.0131+2.6136E-05*72—6.19SSE-09*72A2)

39900 RETURN

40000 REM ****convection coefficient for enc10$ed space******

40100 TAVE=(TA1+TA2)/2-460

40200 NPR=.07206-1.8453E-04*TAVE+3.5979E-07*TAVEA2-2.9967E-
10*TAVEA3

40300 KAIR=.0131+2.6136E-05*TAVE-6.1955E-09*TAVEA2

40400 L1(I,1)=L(2,3)

40500 NGR=L1(I,1)A3*ABS(TA2-TA1)*(.1794+4.1794*EXP(-
.0095034*TAVE))*1E6

40600 REM convection coefficient of air for Top

40700 IF NGR<=4E+05 THEN GOTO 40800

40800 A=.195:M=1/4:GOTO 40900

40900 _A=.068:M=1/3

41000 H(33)= KAIR/L1(1,1)*A*NGRAM

41025 hradea)§E(i)*b012*((TA2/100)A4-(TA1/100)A4)/(TA2-TA1)

41050 h(JJ)= (33)

41100 RETURN

50000 Rem ******Ca1cu1ation of Fina1 va1ues of h's and k's******
51100 Ta1=Temp1(1): Ta2=Tinfinite: i=3: jj=1: osub 30100
51200 'Ta1=T1ce: Ta2=Ta(2): i=2: jj=4: gosub 29 20
51600 Ta1=Temp1(2): Ta2=Temp(1): gosub 35000: k14=k
51700 Ta1=Temp1(2): Ta2=Tg: osub 35000: k25=k
51800 Ta1=Tice: Ta2=T : gosu 35000: k36=k
51810 SSE(zzzz)=(Qexp%zzzz)-Q)A2
. ererr(zzzz)=(Qexp(zzzz)-Q)/Qexp(zzzz)*100
518 0 SSE(dataset+1)=SSE(dataset+1)+SSE(zzzz)
52300 Print #1, using " ###.#### ";
zzzz;A(1);A(2);A(4);A( );A(7);A(8):x;h(1):hrad(1);h(4):hrad(4):k
14;k25;k36;Qexp(zzzz);Q;SSE(zzzz);pererr(zzzz): Print #1,
59000 Return

60000 Rem **********Ca‘lcuat-ion of SSE ********************~k****
62700 Return

125

(KCARG17.BAS)
10 OPEN "O",#1,"kcargoOG.res"

20 READ CRITERIA, MAXITER

25 DATA 0.005.100

27 Print #1, " DataSet A1 A2 A4 A5 A7 A8 x
h1 hradl h4 hrad4 h5 hrad5 h8 hrad8 k14 k25
k36 Qexp Qsim ErrorAZ %error "

28 Print #1,

30 N= 20: Unknown=11: inicon=0.015

50 DIM T(N), TA(N), L(5,3), Ll(5,2), A(N), KCN), HCN), BUFFER(N),
E(N), hradCN)

51 DIM HA(5,2), C(unknown,Unknown+1), VVV(N), vvv1(N), TEMP(N),
TEMPlCN)

53 BOLZ=1. 714E- 09

54 Dataset=4: Dim QexpCDataset), SSECdataset+1), perererataset)
55 For i= 1 to Dataset: read QexpCi): next i

56 data 22.18, 18.86, 20.16, 29.23

61 For zzzz=1 to Dataset

63 Read property, Geometry, Nbucket

64 REM ****reading dimensions and properties of the wa115****
65 FOR I=1 TO 5

70 FOR J=1 TO 3

80 READ X

90 L(I, J): x/12

95 print #1, "L(";i;",";j; ")=";L(i,J)

100 NEXT J

110 READ E(i): print #1, "E(";i;")=";E(i)

120 NEXT I

123 If geometry=0 then 165

125 xi = (L(3,1)-L(2,1))/2:xo = (L(5,1)-L(4,1))/2

130 A(11)=(3.141592*(L(1,1)+L(1,2))/2*L(1,3)+
3.141592/4*L(1,1)A2)*Nb ucket

135 A(12)=(3.141592/4*L(1,1)A2)*Nbucket: A(13)=A(12):
A(14)=A(12)

140 A(8)=L(2,1)*L(2,2)-
A(12):A(5)=2*L(2,3)*(L(2,1)+L(2,2))+L(2,1)*L(2,2)

145 A(9)=L(3,1)*L(3,2)-
A(12):A(4)=2*L(3,3)*(L(3,1)+L(3,2))+L(3,1)*L(3,2)

150 A(7)=L(4,1)*L(4,2)-
(A(12)+A(9)):A(2)=2*L(4,3)*(L(4,1)+L(4,2))+L(4,1)*L(4,2)
155 A(10)=L(5,1)*L(5,2 -
(A(12)+A(9)):A(1)=2*L(5,3)*(L(5,1)+L(5,2))+L(5,1)*L(5,2)
160 goto 190

165 xi = (L(3,1)-L(2,1))/2:xo = (L(S,1)-L(4,1))/2

170 A(11)=(2*L(1,3)*(L(1,1)+L(1,2))+L(1,1)*L(1,2))*Nbucket
175 A(12)=(LC1,1)*L(1, 2))*Nbucket: A(13)=A(12): A(14)=A(12)
180 A(8)=L(2,1)*L(2,2

A(12): A(5L 2*L(2, 3)*(L(2,1)+L(2,2))+L(2,1)*L(2,2)

182 A(9)=L(3,1)*L(3,2-

A(12): A(4L 2*L(3, 3)*(L(3,1)+L(3,2))+L(3,1)*L(3, 2)

126

184 A(7)=L(4,1)*L(4,2
(A(12)+A(9)):A(2)=2*L(4,3)*(L(4,1)+L(4,2))+L(4, 1)*L(4, 2)
186 A(10)=L(5,1)*L(L
(A(12)+A(9)):A(1)=2*L(S, 3)*(L(5,1)+L(L 2))+L(5,1)*L(L 2)

190 REM ******************************************************

195 for i= 1 to 14: print #1, "A(";i;")=";A(i):next i

200 REM ********read1ng known temperature*********************
210 Tice=32:

215 Read Tinfinite . . . . . .

220 T1ce§T1cef460: T1nf1n1te=T1nf1n1te+460

225 Tg=T1nf1n1te

23o REM *******************************************************

250 REM *********initia1 uess of each temperature*************
255 DELTA=(Tinfinite-Tice /7

260 T(7)=Tinfinite: T(8)=Tinfinite: T(9)=Tinfinite:
T(10)=Tinfinite:

270 T(1)=Tinfinite- De1ta: T(2)=Tinfinite- 2*De1ta:
T(3)=Tinfinite- 3*De1ta:

280 T(4)=Tinfinite- 4*De1ta: T(5)=Tinfinite- 5*De1ta:
T(6)=Tinfinite- 6*De1ta

290 T(11)= 10

295 For i= 1 tON N: BufferCi)=T(i): TA(i)=T(i): print #1,
"TA(";i; ")="; Ta(iL next I

300 REM *******************************************************

350 FOR KK=1 TO MAXITER

380 GOSUB 2000

400 FOR J: 1 TO unknown

410 VVV(J) = T(J)-BUFFER(J)

420 NEXT J

450 FOR J= 1 TO unknown

460 TBUFFER = TBUFFER + ABS((T(3)-BUFFER(J))/T(J))
470 NEXT J

500 FOR III= 1 TO 9

510 FOR J: 1 TO unknown

520 TA(J)=BUFFER(J)+VVV(J)*III/10

530 NEXT 3

S40 GOSUB 2000

550 FOR J: 1 TO unknown

560 TERR = TERR + ABS((T(J)-TA(J))/T(J))
570 NEXT J

580 IF TERR >= TBUFFER THEN 660

590 FOR J= 1 TO unknown

600 TEMP(J)=TA(J): TEMP1(3)=T(J)

610 NEXT 3

650 TBUFFER = TERR: Q=QTEMP

660 TERR=0

670 NEXT III

710 FOR J=1 TO unknown: TA(J)=TEMP(J):BUFFERCJ)=TA(J): NEXT
720 IF TBUFFER < CRITERIA THEN 930

730 TBUFFER=0

750 NEXT KK
910 REM did NOT converge

127

920 PRINT "not converged": GOTO 1998

930 REM converged

964 gosub 6000 Rem Ca1cu1ation of SSE

965 REM ***************Ca1cu1ation of h and k*******
970 gosub 50000

972 next 2222

973 Print #1, "Tota1 SSE=";SSE(dataset+1): sseCdataset+1)=0

974 REM *******************DATA*******************************

1020 DATA 0,0,1, 7.625, 7.625, 9, 0.95, 19.4375,
12.3215,11.9375,0.95, l9.625,12.5,12.125,0.05, 22.438, 15.188,
15.438, 0.95, 22.75, 15.5, 15.75, 0.95, 59

1030 DATA 0,0,1, 7.625, 7.625, 9, 0.95, 19.4375,
12.3215,11.9375,0.05, 19.625,12.5,12.125,0.05, 22.438, 15.188,
15.438, 0.95, 22.75, 15.5, 15.75, 0.95, 59

1040 DATA 0,0,1, 7.625, 7.625, 9, 0.95, 19.5, 12.375,
12.00 ,0.05, 19.625,12.S,12.125,0.05, 22.438, 15.188, 15.438,
0.95, 22.75, 15.5, 15.75, 0.95, 59

1050 DATA 0,0,1, 7.625, 7.625, 9, 0.95, 22.2215,
14.9125,14.4375,0.9S, 22.400,15.1,14.750,0.05, 22.438, 15.188,
15.438, 0.95, 22.75, 15.5, 15.75, 0.95, 59

1997 REM *******************************************************

1998 CLOSE #1

1999 END

2000 REM *************main subroutine*************************
2050 REM **************Loop for materia1 properties***********
2100 Ta1=Ta(1) . Ta2=Tinfinite: i=5: jj=1 . gosub 30100
2200 Ta1=Ta(3) : Ta2=Ta(2) : i=4: 11:2 : gosub 29020
2300 Ta1=Ta(4) : Ta2=Ta(3) i=3: 11:4 : gosub 29020
2400 Ta1=Ta(6) : Ta2=Ta(5) i=2: 11:5 : gosub 29020
2500 Ta1=Tice : Ta2=Ta(6) i=1: 11:11 : gosub 29020
2550 Ta1=Ta(3) : Ta2=Ta(7) i=4: 11:7 : gosub 40000
2600 Ta1=Ta(6) : Ta2=Ta(8) : i=2: 00:8 gosub 40000
2620 Ta1=Ta(1) :Ta2=Ta(2) : gosub 35 0: k12=k

2640 Ta1=Tg :Ta2=Ta(7) : gosub 35000: kg7=k

2660 Ta1=Ta(4) :Ta2=Ta(5) : gosub 35000: k45=k

2680 Ta1=Tg :Ta2=Ta(9) : gosub 35000: k =k

2700 Ta1=Ta(9) :Ta2=Ta(8) : gosub 35000: k 9=k

2720 Ta1=Tg :Ta2=Ta(10) : gosub 35000: k 10=k

2740 Ta1=Ta(10):Ta2=Tice : gosub 35000: k Oice=k

3000 REM *********Loop for construction of Matrix************
3100 Gosub 10000

4000 REM *****impose known variab1es and reduce matrix********
5000 REM ******************So]ving matrix********************
5100 GOSUB 20000

128

6000 REM *****************************************************

8400 QTEMP=C(unknown,unknown+1)
8500 FOR I=1 TO unknown: T(I)=C(I,unknown+1): NEXT I

8600 FOR I=1 TO unknown: FOR J=1 TO unknown+1:C(I,J)=0:NEXT J:
NEXT I

8700 FOR I=1 TO N: HCI)=0: NEXT I
9000 RETURN

10000 REM ***********Constructing Matrix**********************

11000 C(1,1)= -h(1)*A(1)-k12*SQR(A(1)*A(2))/XO:
C(1,2)=k12*SQR(A(1)*A(2))/xo
11100 C(1,12)= -h(1)*A(1)*Tinfinite

11200 C(2,3)= h(7)*A(7): C(2,7)= -h(7)*A(7)—
kg7*SQR(A(7)*A(10))/xo
11 00 C(2,12)= -kg7*SQR(A(7)*A(10))/xo*Tg

11400 C(3,2)= h(2)*A(2): C(3,3)=-h(2)*A(2)-h(7)*A(7)-h(4)*A(4):
C(3,4)= h(4)*A(4):C(3,7)= h(7)*A(7)
11500 c(3,12)=0

11600 C(4,3)= h(4)*A(4): C(4,4)= -h(4)*A(4)-
k45*SQR(A(4)*A(S))/Xi: C(4,5)= k45*SQR(A(4)*A(5))/xi
11700 C(4,12)= 0

11750 C(5,6)= h(5)*A(5): C(S,5)= -h(5)*A(5)-
k45*SQR(A(4)*A(5))/xi: C(5,4)= k45*SQR(A(4)*A(5))/xi
11770 C(4,12)= 0

11800 C(6,8)= k89*SQR(A(8)*A(9))/xi: C(6,9)= -
k89*SQR(A(8)*A(9))/xi-kg9*A(9)/XO
11900 C(6,12)= -k99*A(9)/xo*Tg

12000 C(7,6)= h(8)*A(8): C(7,8)= -h(8)*A(8)-
k89*SQR(A(8)*A(9))/xi: C(7,9)= k89*SQR(A(8)*A(9))/xi
12100 C(7,12)= 0

12200 C(8,6)= h(11)*A(11): C(8,10)= k101ce*A(12)/Xi: C(8,11)= -1
12300 C(8,12)= (h(11)*A(11)+k101ce*A(12)/xi)*Tice

12400 C(9,10)= —k910*A(12)/xo -k10ice*A(12)/xi
12500 C(9,12)= -kglO*A(12)/xo*Tg -k10ice*A(12)/xi*Tice

12600 C(10,5)= h(5)*A(5): c(10,6)=-h(5)*A(5)-h(8)*A(8)-
h(11)*A(11): C(10,8)= h(8)*A(8)

12700 C(10,12)=-h(11)*A(11)*Tice

12800 c(11,1)= k12*SQR(A(1)*A(2))/xo:c(ll,2): -h(2)*A(2)-
k12*SQR(A(1)*A(2))/xo:c(11,3)=h(2)*A(2)

12900 C(11,12)= 0

13000 FOR I=1 TO unknown:

129

13100 FOR J: 1 TO unknown+1
13200 PRINT #1, C(I,J);
13300 NEXT J

13350 PRINT #1,

13400 NEXT I

13500 PRINT #1,

14000 Return

20000 REM **********subroutine for so1ving matrix**************

20010 REM: SOLVES SYSTEM OF EQNS BY GAUSSIAN ELIMINATION

20060 FOR I=1 TO unknown : PVT=I 'e1imination using row I ...

20070 FOR K=I+1 TO unknown

20080 IF ABSCCCK,I))>ABS(C(PVT,I)) THEN PVT=K

20090 NEXT K 'row "PVT" has 1argest e1ement in Ith co1umn

20100 FOR J=I TO unknown+1 'interchan e row I and pivot row

20110 CHG=C(I,J) : C(I,J)=C(PVT,J§ : C(PVT,J)=CHG

20120 NEXT J 'now row I has 1argest e1ement in Ith co1umn

20130 FOR K=1 TO unknown : IF K=I THEN 20160 'zero out

20140 F=C(K,I)/C(I,I)'amount of row I to subtract off row K:

20150 FOR J=I+1 TO unknown+1 : C(K,J)=C(K,J)-F*C(I,J) : NEXT 3

20160 NEXT K 'Ith co1umn now zeroed out (except for row I)

20170 NEXT I

20180 FOR K=1 TO unknown: C(K,unknown+1)=C(K,unknown+1)/C(K,K)
:print #1, "T(";k;")=";c(k,unknown+1):NEXT K

20200 RETURN

29020 REM ****convection coefficient for enc1osed space******

29050 TAVE=(TA1+TA2)/2-460

29060 NPR=.07206-1.8453E-04*TAVE+3.5979E—07*TAVEA2-2.9967E—

10*TAVEA3

29065 KAIR=.0131+2.6136E-05*TAVE-6.1955E—09*TAVEA2

29066 L1(I,1)=ABS((L(2,3)-L(1,3)))

29068 L1(I,2)=(ABS((L(2,1)-Nbucket*L(1,1))/2) + ABS((L(3,2)-
Nbucket*L(1,2))/2))/2

29070 FOR J=1 TO 2

29080 NGR(3)=L1(I,J)A3*ABS(TA2-TA1)*
(.1794+4.1794*EXP(-.0095034*TAVE))*1E6

29090 IF J=1 THEN GOSUB 29400

29100 IF J=2 THEN GOSUB 29200

29110 HA(1,J)= KAIR/LlCI,J)*A*NGR(j)AM

29120 NEXT J

29122 If 2222 > 12 and 2222 < 17 then 29130

29124 If jj= 8 then
h(jj)=(HA(1,1)*3.141592/4*L(i,1)A2+HA(i,2)*3.141592*L(i,1)*L(i,3
))/A(gj): oto 29135

29130 (jg =(HA(i,1)*L(i,1)*L(i,2)+

HACi.Z)* *LCi.3)*(L(i.1)+L(i.2)))/A( J)

29135 hrad(jg)§E(i)*bo1;f(TA244-TA1A4 /(TA27TA1) ._

29140 h(JJ)= (JJ)+hrad(JJ)=pr1nt #1. "h(";33;")=":h(JJ)

29150 RETURN

29200 REM convection coefficient of air for vertica1 p1ane

29210 IF NGRC§)<=20000 THEN GOTO 29270

29230 A=.1 *(LCI,3)/L1(i,j))A(-1/9):M=1/4:GOTO 29290

29270 A=.065*(L(I,3)/L1(i,j))A(-1/9):M=1/3

29290 RETURN

29300 REM convection coefficient of air for Bottom

130

29305
29310
29330
29350
29390
29400
29410
29430
29450
29490

30100

A1(1.J)?A(1.J)-A(1.J)
IF NGRC )<=4E+05 THEN GOTO 29350
A=.1 5:M=1/4:GOTO 29390
A=.068:M=1/3
RETURN
REM convection coefficient of air for Top
IF NGRC')<=4E+05 THEN GOTO 29450
A=.1 5:M=1/4:GOTO 29490
A=.068:M=1/3
RETURN

REM subroutine for conductivity and FREE convection

coefficient for air

30500 TAVE=(TA1+TA2)/2-460

30600 NPR=.07206-1.8453E-04*TAVE+3.5979E-07*TAVEA2-2.9967E-
10*TAVEA3

30650 KAIR=.0131+2.6136E-05*TAVE-6.1955E-09*TAVEA2

30660 L1(I,1)=(L(I,1)+L(I,2))/2

30680 L1(I,2)=L(I,3)

30700 FOR J=1 TO 2

30800 NGRCJ)=L1(I,J)A3*ABS(TA2-TA1)*(.1794+4.1794*EXP(-
.0095034*TAVE))*1E6

30850 NPRNGR=NPR*NGR(J)

30900 IF J=1 THEN GOSUB 34000

31050 IF J=2 THEN GOSUB 32000

31100 HA(1,J)= KAIR/LlCI,J)*A*NPRNGRAM

31200 NEXT 3

31300

HACi,

31350
31400
31500

32000
32100
32200
32300
32500
32700
32900

33000
33100
33300
33500
33900
34000
34500
34900

35000
35500
35550
35560
35570
35580
35590

h(gg)=(HA(i.1)*L(i.1)*L(i.2)+

2) *LCi.3)*(L(i.1)+L(i.2)))/A( j)
hradeg)§E(i)*bo1;f(TA2A4-TA1A4 /(TA27TA1) _.
h(JJ)= (JJ)+hrad(JJ) pr1nt #1. "h(":JJ;")=";h(JJ)
RETURN

REM convection coefficient of air for vertica1 p1ane

IF NPRNGR<1E+04 THEN GOTO 32500
IF NPRNGR>1E+09 THEN GOTO 32700
A=.59:M=1/4:GOTO 32900
A=1.36:M=1/5:GOTO 32900
A=.13:M=1/3

RETURN

REM convection coefficient of air for Bottom

IF NPRNGR<2E+07 THEN GOTO 33500
A=.14:M=1/3:GOTO 33900
A=.S4:M=1/4

RETURN

REM convection coefficient of air for Top
A=.58:M=1/3

RETURN

REM ******conduction coefficients for each materia15*****
TAVE=(TA1+TA2)/2-460
If property = 2 then 35650
If property = 1 then 35660
If property = 0 then 35670
If property = 3 then 35680
If property = 4 then 35690

131

35640 k=0.346: goto 36000 'water

35650 k=inicon:goto 36000 'air

35660 k=0.035: goto 36000 'corrugated board
35670 k=0.022: goto 36000 'EPS foam

35680 =0.018: goto 36000 'PoTyurathane

35690 k=0.010: goto 36000 'VIP

36000 k=k*(.013 +2.6136E—05*TAVE-6.19SSE-09*TAVEA2)/
(.0131+2.6136E—05*72-6.19SSE-O9*72A2)

36900 RETURN

40000 REM ****convection coefficient for enCIOsed space******

40100 TAVE=(TA1+TA2)/2-460

40200 NPR=.O7206-1.8453E-04*TAVE+3.5979E-07*TAVEA2
-2.9967E-10*TAVEA3

40300 KAIR=.0131+2.6136E-05*TAVE-6.1955E-09*TAVEA2

40400 L1(I,1)=L(2,3)

40500 NGR=L1(I,1)A3*ABS(TA2-TA1)*

(.1794+4.1794*EXP(-.0095034*TAVE))*1E6

40600 REM convection coefficient of air for Top
40700 IF NGR<=4E+05 THEN GOTO 40800
40800 A=.195:M=1/4:GOTO 40900
40900 A=.068:M=1/3
41000 Hij)= KAIR/L1(1,1)*A*NGRAM
41025 hradea);E(i)*b014f(TAZA4-TA1A4)/(TA27TA1) ..
41050 h(JJ)= (JJ)+hrad(JJ)=pr1nt #1. "h(";JJ;")=";h(JJ)
41100 RETURN
50000 Rem ****Ca1cu1ation of Fina1 va1ues of h's and k's******
51100 Ta1=Temp1(1): Ta2=Tinfinite: i=3: jj=1: osub 30100
51200 Ta1=Tice: Ta2=Ta(4): i=1: jj=8: gosub 29 20
51300 Ta1=Temp1(4): Ta2=Temp1(2): i=2: jj=4: gosub 29020
51400 Ta1=Temp1(4): Ta2=Temp1(3): i=2: gg=5: gosub 40000
51600 Ta1=Temp1(4): Ta2=Temp(1): gosub 000: k14=k
51700 Ta1=Temp1(3): Ta2=Tg: osub 35000: k25=k
51800 Ta1=Tice: Ta2=T : gosu 35000: k36=k
51810 SSEszzz)=(Qexp%zzzz)-Q)A2
: ererr(zzzz)=(Qexp(zzzz)-Q)/Qexp(zzzz)*100
518 0 SSE(dataset+1)=SSE(dataset+1)+SSE(zzzz)
52300 Print #1, usin " ###.#### ";

zzzz;A(1);A(2);A(4);A 5);A(7);A(8):x;h(1);hrad(1);h(4);hrad(4):
h(5):hradCS);h(8):hrad(8);k14;k25;k36;Qexp(zzzz);Q;SSE(zzzz):Pe
rerrszzz): Print #1,

59000

60000
62700

Return

Rem ***********Ca1cuation Of SSE *********************
chun

132

APPENDIX B. Thickness of solid layers of the multi-layer wall

 

 

 

Data The thickness of the each wall (in)
Set
2 3 4 5 6 7 8 10

1 1.500

2 3.000

3 1.500 1.500

4 1.500 1.500

5 1.500 1.500 1.500

6 0.156 0.156 0.156

7 0.156 0.156 0.156 0.156 0.156
8 0.156 0.156 0.156 0.156 0.156
9 0.156 0.156 0.156 0.156 0.156
10 0.156 0.156 0.156 0.156 0.156
11 0.156 0.156 0.156 0.156 0.156
12 0.156 0.156 0.156 0.156 0.156
13 1.500

14 3.000

15 1.500 1.500

16 1.500 1.500

17 1.500 1.500 1.500

18 0.156 0.156 0.156

19 0.156 0.156 0.156 0.156 0.156
20 0.156 0.156 0.156 0.156 0.156
21 0.156 0.156 0.156 0.156 0.156
22 0.156 0.156 0.156 0.156 0.156
23 0.156 0.156 0.156 0.156 0.156
24 0.156 0.156 0.156 0.156 0.156
25 1.500

26 3.000

27 1.500 1.500

28 1.500 1.500

29 1.500 1.500 1.500

30 0.156 0.156 0.156

31 0.156 0.156 0.156 0.156 0.156
32 0.156 0.156 0.156 0.156 0.156

 

133

APPENDIX B. (Cont’d)

 

 

 

Data The thickness of the each wall (in)
Set
2 3 4 5 6 7 8 9 10
33 0.156 0.156 0.156 0.156 0.156
34 0.156 0.156 0.156 0.156 0.156
35 0.156 0.156 0.156 0.156 0.156
36 0.156 0.156 0.156 0.156 0.156
37 1.500
38 3.000
39 1.500 1.500
40 1.500 1.500
41 1.500 1.500 1.500
42 0.156 0.156 0.156
43 0.156 0.156 0.156 0.156 0.156
44 0.156 0.156 0.156 0.156 0.156
45 0.156 0.156 0.156 0.156 0.156
46 0.156 0.156 0.156 0.156 0.156
47 0.156 0.156 0.156 0.156 0.156
48 0.156 0.156 0.156 0.156 0.156
49 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156
50 1.000 0.156 0.156 0.156 0.156 0.156 1.000 0.156 0.156 0.156
51 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156
52 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156
53 0.500 0.156 0.156 0.500 0.156 0.500 0.156 0.156 0.500 0.156
54 0.500 0.156 0.156 0.500 0.156 0.500 0.156 0.156 0.500 0.156
55 1.500 1.500
56 1.500 1.500
57 0.156 0.156 0.156
58 0.156 0.156 0.156 0.156 0.156
59 0.156 0.156 0.156 0.156 0.156
60 0.156 0.156 0.156 0.156 0.156
61 0.156 0.156 0.156 0.156 0.156
62 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156
63 1.000 0.156 0.156 0.156 0.156 0.156 1.000 0.156 0.156 0.156
64 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156

 

134

APPENDIX B. (Cont’d)

 

 

 

Data The thickness of the each wall (in)
Set

2 3 4 5 6 7 8 9 1O
65 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.500
66 0.500 0.156 0.156 0.500 0.156 0.500 0.156 0.156 0.500 0.156
67 1.500
68 3.000
69 1.500 1.500
70 1.500 1.500
71 1.500 1.500 1.500
72 0.156 0.156 0.156
73 0.156 0.156 0.156 0.156 0.156
74 0.156 0.156 0.156 0.156 0.156
75 0.156 0.156 0.156 0.156 0.156
76 0.156 0.156 0.156 0.156 0.156
77 0.156 0.156 0.156 0.156 0.156
78 0.156 0.156 0.156 0.156 0.156
79 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156
80 1.000 0.156 0.156 0.156 0.156 0.156 1.000 0.156 0.156 0.156
81 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156
82 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156
83 0.500 0.156 0.156 0.500 0.156 0.500 0.156 0.156 0.500 0.156
84 0.500 0.156 0.156 0.500 0.156 0.500 0.156 0.156 0.500 0.156
85 1.500 1.500
86 1.500 1.500
87 0.156 0.156 0.156
88 0.156 0.156 0.156 0.156 0.156
89 0.156 0.156 0.156 0.156 0.156
90 0.156 0.156 0.156 0.156 0.156
91 0.156 0.156 0.156 0.156 0.156
92 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156
93 1.000 0.156 0.156 0.156 0.156 0.156 1.000 0.156 0.156 0.156
94 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156
95 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.156 0.500
96 0.500 0.156 0.156 0.500 0.156 0.500 0.156 0.156 0.500 0.156

 

135

APPENDIX C. The thickness of the air space between solid layers

 

 

 

D N of The thickness of the each air space (in)
ata the
set walls 1 2 3 4 5 6 7 8
1 1 -
2 1 -
3 2 0.100
4 2 0.300
5 3 0.200 0.500
6 3 0.100 0.400
7 5 0.200 0.100 0.300 0.700
8 5 0.300 0.200 0.300 0.600
9 5 0.400 0.300 0.200 0.500
10 5 0.500 0.300 0.200 0.400
11 5 0.600 0.200 0.100 0.300
12 5 0.700 0.100 0.100 0.200
13 1 -
14 1 -
15 2 0.100
16 2 0.300
17 3 0.200 0.500
18 3 0.100 0.400
19 5 0.200 0.100 0.300 0.700
20 5 0.300 0.200 0.300 0.600
21 5 0.400 0.300 0.200 0.500
22 5 0.500 0.300 0.200 0.400
23 5 0.600 0.200 0.100 0.300
24 5 0.700 0.100 0.100 0.200
25 1 -
26 1 -
27 2 0.100
28 2 0.300
29 3 0.200 0.500
30 3 0.100 0.400
31 5 0.200 0.100 0.300 0.700
32 5 0.300 0.200 0.300 0.600

 

136

APPENDIX C. (Cont’d)

 

 

 

D N of The thickness of the each air space (in)

ata the

set walls 1 2 3 4 5 6 7 8 9
33 5 0.400 0.300 0.200 0.500

34 5 0.500 0.300 0.200 0.400

35 5 0.600 0.200 0.100 0.300

36 5 0.700 0.100 0.100 0.200

37 1 -

38 1 -

39 2 0.100

40 2 0.100

41 3 0.100 0.500

42 3 0.100 0.200

43 5 0.100 0.200 0.100 0.200

44 5 0.100 0.200 0.100 0.200

45 5 0.100 0.200 0.100 0.200

46 5 0.100 0.200 0.100 0.200

47 5 0.100 0.200 0.100 0.200

48 5 0.100 0.100 0.100 0.100

49 10 0.100 0.100 0.100 0.100 0.100 0.100 0.100 0.100 0.100
50 10 0.100 0.100 0.100 0.100 0.100 0.100 0.100 0.100 0.100
51 10 0.100 0.100 0.100 0.100 0.100 0.100 0.100 0.100 0.100
52 10 0.100 0.100 0.100 0.100 0.100 0.100 0.100 0.100 0.100
53 10 0.100 0.100 0.100 0.100 0.100 0.100 0.100 0.100 0.100
54 10 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200
55 2 1.000

56 2 1.000

57 3 0.100 0.200

58 5 0.100 0.200 0.100 0.200

59 5 0.100 0.200 0.100 0.200

60 5 0.100 0.200 0.100 0.200

61 5 0.100 0.200 0.100 0.200

62 10 0.100 0.200 0.100 0.100 0.100 0.100 0.100 0.100 0.100
63 10 0.100 0.100 0.100 0.100 0.100 0.100 0.100 0.100 0.100
64 10 0.100 0.100 0.100 0.100 0.100 0.100 0.100 0.400 0.400

 

137

 

APPENDIX C. (Cont’d)

 

 

 

D N of The thickness of the each air space (in)

ate the

39‘ walls 1 2 3 4 5 6 7 8 9
65 10 0.100 0.100 0.100 0.100 0.100 0.100 0.400 0.400 0.200
66 10 0.100 0.100 0.100 0.100 0.100 0.100 0.100 0.100 0.100
67 1 -

68 1 -

69 2 0.100

70 2 0.100

71 3 0.100 0.500

72 3 0.100 0.200

73 5 0.100 0.200 0.100 0.200

74 5 0.100 0.200 0.100 0.200

75 5 0.100 0.200 0.100 0.200

76 5 0.100 0.200 0.100 0.200

77 5 0.100 0.200 0.100 0.200

78 5 0.100 0.100 0.100 0.100

79 10 0.100 0.100 0.100 0.100 0.100 0.100 0.100 0.100 0.100
80 10 0.100 0.100 0.100 0.100 0.100 0.100 0.100 0.100 0.100
81 10 0.100 0.100 0.100 0.100 0.100 0.100 0.100 0.100 0.100
82 10 0.100 0.100 0.100 0.100 0.100 0.100 0.100 0.100 0.100
83 10 0.100 0.100 0.100 0.100 0.100 0.100 0.100 0.100 0.100
84 10 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200
85 2 1.000

86 2 1.000

87 3 0.100 0.200

88 5 0.100 0.200 0.100 0.200

89 5 0.100 0.200 0.100 0.200

90 5 0.100 0.200 0.100 0.200

91 5 0.100 0.200 0.100 0.200

92 10 0.100 0.200 0.100 0.100 0.100 0.100 0.100 0.100 0.100
93 10 0.100 0.100 0.100 0.100 0.100 0.100 0.100 0.100 0.100
94 10 0.100 0.100 0.100 0.100 0.100 0.100 0.100 0.400 0.400
95 10 0.100 0.100 0.100 0.100 0.100 0.100 0.400 0.400 0.200
96 10 0.100 0.100 0.100 0.100 0.100 0.100 0.100 0.100 0.100

 

138

 

 

 

 

mod mod mod mod mod mod mod mod vm
mod mod mod mod mod mod mod mod mm
mod mod mod mod mod mod mod mod Nu
mod mod mod mod mod mod mod mod 3
mod mod mod mod mod mod mod mod om
mod mod mod mod mod mod mod mod 3
mod mod mod mod 3
mod mod mod mod 5
mod mod 0..
mod mod 3
3
.2
mad mod mod mod mad mod mad mod N..
mad mod mad mod mad mod mad mod 3
mad mod mad mod mod mod mad mod or
86 mod mad mod mad mod mwd mod a
mad mod mad mod mad mod mad mod m
mad mod mad mod mad mod mad mod 5
mad mod mad mod m
mod mod mad mod m
mad mod v
mmd mod m
N
F

No 5 No 5 «w 5 «w 5 «w 5 Na 5 mm 5 mm 5 «w 5
“mm
x. o m v m N . _. Ema

08% .fi 56.85-5053. 9.: ho conEac

coma» ..fi ...o>a.-o£.coo§oa. a...“ ho 3333958 on... d x.ozwam<

139

 

 

 

 

mad mad mad mwd mad mad mod mod 3
mad mad mad mad mad mad mmd mmd D4
mad mad mad mad mad mad mad mad me
mad mad mad mad mod mad mad mad mv
mad mad mad mad mod mad mad mmd 3
mad mod mad mod mad mod mod mad m...
mad mmd mad mad Ne
mad mad mad mad 5.
mad mad ov
mmd mad mm
mm
5
mod mod mad mod mod mod mad mod on
mod mod mod mod mod mod mad mod mm
mod mod mad mod mod mod mad mod 3”
mod mod mod mod med mod mmd mod mm
mod mod mad mod mod mod mad mod Nm
mod mod mod mod mod mod mad mod 5
mod mod mad mod on
mad mod mod mod N
mod mod N
mod mod N
N
N

Na Fm Nu .w No Pm Nw aw Nu Fm Nw .w we so No aw Nw Fm
“mm
x. o m v m N _. 2mm.

.8QO .__m 165.85.52.15; 9: ho .0983:

3.2.58 .0 X_n_2m_n_n_<

140

 

 

 

 

mod mad mod mod Nu

mad mad mod mad K

mod mod on

mmd mad mo

8

no

mad mad mad mad mad mad mad mad mad mad mad mod mod mad mad mad mad mad mo

mad mad mad mad mad mad mad mod mad mad mad mad mad mod mad mod mad mad mo

mad mad mad mad mad mad mad mad mad mad mad mad mad mad mad mad mad mad 3

mad mad mod mad mod mad mad mad mad mad mad mad mad mad mad mad mad mad mo

mad mad mad mod mad mm mad mad mad mad mod mad mad mmd mad mad mad mod No

mad mad mad mad mad mad mad mad Po

mmd mad 3.0 mad mad mad mad mad co

mad mad mwd mad mad mad mad mad mm

mad mad mad mad mmd mad mad mad mm

mod mad mad mad mm

mad mad om

mad mad mm

mad mad mad mad mad mad mad mad mod mad mad mad mad mad mad mad mod mmd vm

mad mad mad mad mod mad mad mad mad mod mad mad mad mmd mmd mad mad mad mm

mmd mad mad mad mad mad mad mad mad mad mad mad mmd mad mod mad mad mad Nm

mad mad mmd mad mad mad mod mad mad mad mad mad mad mad mad mad mad mad 5

mad mod mad mad mad mad mad mad mad mad mad mmd mad mmd mad mad mad mad om

mad mad mod 36 mad mmd mad mad mod mad mad mad mad mad mad mm 86 mad av
Nw K... Na 5 Nu 5 Nu 5 Nw 5 Nu 5 Nu «m Nw rm Nw Kw

«mm

m w n o m V m N w Ema

83m .6 453.105.5053. 05 ho .383:

8.208 .0 x_ozman_<

141

 

 

 

mad mod mad mod mad mod mad mod mad mod mmd mod mad mod mad mod mod mod om
mod mad mod mad mad mod mod mod mad mod mod mad mad 36 mod mad mmd mad mm
mod mad mod 36 mad mad mad mod mad mod mad mod mod mod mad mod mad mod .3
mod mod mod mad mad mod mod mod mod mad mad mad mod mod mod mmd mad mad mm
mad mad mad mod mad mm mad mad mad med mod mod mod mod no.0 mod mod mad Nm
mod mod mod mod mod mod mod mod 5
mad mad mad mad mod mad mad mmd om
mad mod mod mod mmd mod mad mod mm
mod mod mod mod mod mod mod mod mm
mod mod mod mod 3
mod mod ow
, mad mad mm
mod mad med mad mad mad mod mad mad mod mod mod mod mad mad mod mod mod vm
mod mad mod mad mad mad mod mad mad mod mod mod mod mad mad mod mod mod mm
mod mad mod mmd mad mod mod mad mad mad mad mad mmd mad mad mod mod mod Nm
mod mad mod mod mod mod mod mod mod mod mod mod mod mod mod mod mod mod 5
mod mad mad mad mad mad mad mad mad mad mod mad mmd mad mad mad mad mad om
mad mod mad mad mad mm mad mad mad mad mad mad mad mad mod mm mad mad as
mod mad mad mad mad mad mad mad mu
mod mad mod mad mod mod mad mad t.
mad mad mad mad mod mmd mod mad on
mad mad mod mod mmd mad mod mad mu
mod mod mod mod mad mod mod mod E
mad mod mod mad mad mad mod mod MK

Nw 5 Nu». ww Nw _.w Nw 5 Nw Fm Nw Pm Nw Fw Nw 5 N» 5
“mm
m w n o m v m N _. Ema

woman Em 105.85.385.89. 05 ho .383:

95:09 .0 xazwnE<

142

APPENDIX E. R-values and parameters generated by computer program

 

Liner Part Non-Liner Part
Data

 

Set ths thaI tha2 tha3 (when N>O)

 

l 2 3 4 5 6 7 8 9
1 5.891 1.50 0.00 0.00
2 11.757 3.00 0.00 0.00
3 12.262 3.00 0.00 0.10
4 13.296 3.00 0.00 0.30
5 21.178 4.50 0.00 0.70
6 4.461 0.47 0.00 0.50
7 9.732 0.78 0.00 1.30
8 10.761 0.78 0.00 1.40
9 10.333 0.78 0.00 1.40
10 10.328 0.78 0.00 1.40
11 9.748 0.78 0.00 1.20
12 8.702 0.78 0.00 1.10
13 5.851 1.50 0.00 0.00
14 11.702 3.00 0.00 0.00
15 12.269 3.00 0.10 0.00
16 13.356 3.00 0.30 0.00
17 21.206 4.50 0.20 0.50
18 4.569 0.47 0.50 0.00
19 10.077 0.78 1.30 0.00
20 11.161 0.78 1.40 0.00
21 10.676 0.78 1.40 0.00
22 10.673 0.78 1.40 0.00
23 10.106 0.78 1.20 0.00
24 9.004 0.78 1.10 0.00
25 5.851 1.50 0.00 0.00
26 11.702 3.00 0.00 0.00
27 12.269 3.00 0.10 0.00
28 13.356 3.00 0.30 0.00
29 21.206 4.50 0.20 0.50
30 4.562 0.47 0.40 0.10
31 10.047 0.78 1.10 0.20
32 11.033 0.78 0.80 0.60

OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO 2

143

 

APPENDIX E. (Cont’d)

 

 

 

 

Liner Part Non-Liner Part
Data
Set ths thaI tha2 N ”M3 (when N>O)
2 3 4 5 6 7 8 9
33 10.571 0.78 1.00 0.40 0
34 10.488 0.78 0.70 0.70 0
35 9.846 0.78 0.60 0.60 0
36 8.738 0.78 0.30 0.80 0
37 5.851 1.50 0.00 0.00 0
38 11.702 3.00 0.00 0.00 0
39 12.089 3.00 0.00 0.00 1 0.1
40 12.089 3.00 0.00 0.00 1 0.1
41 18.784 4.50 0.00 0.00 2 0.1 0.5
42 2.796 0.47 0.00 0.00 2 0.1 0.2
43 4.982 0.78 0.00 0.00 4 0.1 0.2 0.1 0.2
44 4.982 0.78 0.00 0.00 4 0.1 0.2 0.1 0.2
45 4.982 0.78 0.00 0.00 4 0.1 0.2 0.1 0.2
46 4.982 0.78 0.00 0.00 4 0.1 0.2 0.1 0.2
47 4.982 0.78 0.00 0.00 4 0.1 0.2 0.1 0.2
48 4.588 0.78 0.00 0.00 4 0.1 0.1 0.1 0.1
49 89101.56 0.00 0.00 9 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
50 16.154 3.25 0.00 0.00 9 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
51 92341.56 0.00 0.00 9 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
52 95611.56 0.00 0.00 9 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
53 149312.94 0.00 0.00 9 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
54 16.695 2.94 0.00 0.00 9 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2
55 12.697 3.00 0.00 0.00 1 1.0
56 12.697 3.00 0.00 0.00 1 1.0
57 2.796 0.47 0.00 0.00 2 0.1 0.2
58 4.982 0.78 0.00 0.00 4 0.1 0.2 0.1 0.2
59 4.982 0.78 0.00 0.00 4 0.1 0.2 0.1 0.2
60 4.982 0.78 0.00 0.00 4 0.1 0.2 0.1 0.2
61 4.982 0.78 0.00 0.00 4 0.1 0.2 0.1 0.2
62 94261.56 0.00 0.00 9 0.1 0.2 0.1 0.1 0.1 0.1 0.1 0.1 0.1
63 16.154 3.25 0.00 0.00 9 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
64 96451.56 0.00 0.00 9 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.4 0.4

 

144

APPENDIX E. (Cont’d)

 

 

 

 

Liner Part Non-Liner Part
Data
Set ths thaI tha2 N ”’03 (When N>O)

1 2 3 4 5 6 7 8 9

65 11.893190 0.00 0.00 9 0.1 0.1 0.1 0.1 0.1 0.1 0.4 0.4 0.2
66 14.931 2.94 0.00 0.00 9 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
67 5.851 1.50 0.00 0.00 0
68 11.702 3.00 0.00 0.00 0
69 12.089 3.00 0.00 0.00 1 0.1
70 12.269 3.00 0.10 0.00 0
71 18.784 4.50 0.00 0.00 2 0.1 0.5
72 3.475 0.47 0.00 0.30 0
73 5.331 0.78 0.00 0.20 2 0.2 0.2
74 6.382 0.78 0.40 0.20 0
75 5.331 0.78 0.00 0.20 2 0.2 0.2
76 4.982 0.78 0.00 0.00 4 0.1 0.2 0.1 0.2
77 6.164 0.78 0.00 0.50 1 0.1
78 4.757 0.78 0.00 0.10 3 0.1 0.1 0.1
79 89101.56 0.00 0.00 9 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
80 16.494 3.25 0.00 0.20 7 0.1 0.1 0.1 0.1 0.1 0.1 0.1
81 11.174 1.56 0.80 0.10 0
82 106021.56 0.10 0.50 3 0.1 0.1 0.1
83 16.329 2.94 0.20 0.60 1 0.1
84 20.815 2.94 0.40 1.20 1 0.2
85 12.697 3.00 0.00 0.00 1 1.0
86 15.880 3.00 1.00 0.00 0
87 3.510 0.47 0.30 0.00 0
88 6.410 0.78 0.60 0.00 0
89 6.341 0.78 0.00 0.60 0
90 4.982 0.78 0.00 0.00 4 0.1 0.2 0.1 0.2
91 6.453 0.78 0.60 0.00 0
92 10.853156 0.40 0.30 3 0.1 0.1 0.1
93 17.205 3.25 0.30 0.30 3 0.1 0.1 0.1
94 12.688 1.56 0.00 1.00 2 0.1 0.4
95 146491.90 0.10 0.80 6 0.1 0.1 0.1 0.1 0.1 0.2
96 16.484 2.94 0.00 0.90 0

 

145

REFERENCES

ASTM. (1999a). Standard Practice for Performance Testing of Shipping
Containers and Systems. In ASTM (Eds), Selected ASTM Standards on
Packaging, 5th ed (pp187-197). West Conshohocken, PA: American Society for
Testing and Materials.

ASTM. (1999b). Standard Test Methods for Thermal Insulation Quality of

Packages. ln ASTM (Eds.), Selected ASTM Standards on Packaging, 5th ed
(pp117-120). West Conshohocken, PA: American Society for Testing and
Materials.

ASHRAE. (1993). ASHRAE Handbook. Fundamentals. Atlanta, GA: American
Society of Heating, Refrigerating and Air-Conditioning Engineers.

Brody, L. A., and Marsh, 8. K. (1997). Ih_e Wilev Encyclopedia of Packaging
Technolmy, 2nd ed. New York, NY: John Wiley and Sons.

Burgess, G. (1999). Practical Thermal Resistance and Ice Requirement

Calculations for Insulating Packages. Packaging Technology and Sciencez12, 75-
80.

Burgess, G. (2004). Lecture Notes, Distribution and Packaging Dynamics. East
Lansing, MI: Michigan State University.

Cooper, M. G., Mikic, B. B., and Yovanovich, M. M. (1969). Thermal Contact
Conductance. International Journal of Heat and Mass Transfer: 12, 279-300.

Fried, E. (1963). Thermal Joint Conductance in a Vacuum. American Society of
Mechanical Engineers: 63-AGHT-18.

Geankoplis, C. J. (1983). Transport Processes and Unit Operations, 2nd ed.
Boston, MA: Allyn and Bacon.

Griffith, B., Arasteh, D., and Selkowitz, S. (1991). Gas-filled panel high-

performance thermal insulation. Insulating Materials: Test and Application: 2,
441-454.

Hagen, K. D. (1999). Heat Transfer with Applications. Upper Saddle River, NJ:
Prentice-Hall.

Hernandez, R. J., Selke, S. E. M., and Culter, J. D. (2000). Plastics Packaging

Properties. Processing, Applications and Regplations. Cincinnati, OH: Hanser
Publishers.

146

 

Holman, J. P. (1986). Heat Transfer, 6th ed. New York, NY: McGraw-Hill.

ISTA. (1999). Development Test for Thermal Performance of Insulated Transport
Packaging. In ISTA (Ed), I_STA 1998-1999 ResoprceJBook. East Lansing, MI:
lntemational Safe Transit Association.

Jacob, M. (1949). Heat Transfer. New York, NY: John Wiley and Sons.
Kreith, F. (1973). Principles of Heat Transfer. 3rd ed. New York, NY: Intext Press.

Kositruangchai, N. (2003). Theoretical, Experimental and Computer Model for
Package R-value Using Regular Ice and Dry Ice. Masters Thesis. East Lansing,
MI: Michigan State University.

McCabe, W. L., Smith, J. C., and Harriott, P. (1993). Unit Operations of Chemical
Engineering, 5th ed. New York, NY: McGraw-Hill.

Middleman, S. (1998). An lntrod_uction to Mass and Heat Transfg: Principle of
Analysis and Design. New York, NY: John Wiley and Sons.

Panyaljun, O. (2002). Thermal Insulated Packaging Design with S-Flute
Corrugated Board. PH.D. Dissertation. East Lansing, MI: Michigan State
University.

Sasaki, H., and Kato, E. (1999). Heat Insulating Cardboard Composed of

Corrugated Foamed Polystyrene Layer, Packaging Technology and
Science:19,151-157.

Stavish, L. J. (1984). Designing Insulated Packaging for Perishable in vivo
Diagnostics. Medical Device and Diagnostic Industm: 6 (18), 105-108.

Turner, W. C., and Malloy, F. J. (1981). Thermal Insulation Handbook. New York,
NY: McGraw-Hill.

147

 

     

Iiflifiifljilflij‘iifliii’ji‘jiii

L